龙门石窟奉先寺的红外成像监测与分析——基于MATLAB图像处理技术

刘逸堃，高东亮，马朝龙，范子龙，李心坚

(1． 北京大学考古文博学院，北京 100871； 2． 龙门石窟研究院，河南洛阳 471023)

0 引 言

奉先寺是龙门石窟最具代表性的佛龛之一。保护工作者于1971年起对奉先寺内的九尊佛像进行了抢险加固，缓解了裂隙发育的同时也避免了潜在的坍塌危险[1]。然而，奉先寺石刻的风化仍在持续进行，并且不断受到由降雨造成的渗水的破坏，危害着石刻的安全及其艺术价值[2-5]。随着文物预防性保护工作的全面开展，对于奉先寺保存状况的长期监测及安全性评估显得尤为重要[6]。

红外成像技术可以直观且无损地呈现大面积岩体的红外辐射图像。目前已有利用该技术进行历史建筑湿热相关的分析研究[7-9]，以及其他行业基于红外图像数据分析的研究[10-12]；也有学者尝试将红外成像应用于奉先寺渗水的单次探测中，通过温度差异来反映渗水的情况[13]。而当收集的红外图像数据量十分庞大时，如何从中全面获取与石窟保存情况的信息，则是亟待解决的研究问题。基于MATLAB语言的图像处理技术具有极高的分析深度、复杂度和编程自由度，可以胜任绝大部分图像大数据处理任务，已广泛使用于诸多其他行业的红外图像处理中[14]。本研究旨在将该技术和大数据分析结合，应用于石窟寺红外监测数据的挖掘与分析。

因此，为了实现对部分受到风化及渗水影响的区域进行全面监测与分析，本研究利用固定架设的红外成像设备，对奉先寺部分区域进行长期定间隔红外图像的摄影，并借助MATLAB语言的图像处理技术，对收集到的大量红外图像进行数据统计分析。

1 研究方法

1．1 红外成像监测

1．1．1监测对象 选择奉先寺石造像阿难的头部区域以及中央的卢舍那大佛南侧部分岩体作为本次监测的对象区域。其中，阿难头像(后文简称“头像”)左上部的裂隙处曾经历过灌锚补加固，灌浆材料为呋喃树脂改性环氧树脂浆液[1]。目前可观察到裂隙处有部分表面脱落的现象(图1a)。选择该区域监测，以长期观察阿难头像加固后的保存情况以及可能存在的崩落风险；岩体壁面(图1b)经历过浆砌石砖的砌补[1]，可直接观察到砌补砖体的痕迹；卢舍那大佛南侧肩部上方的岩体存在多处裂隙渗水[2，5，13](图1c)，受渗水径流影响的岩体表面发白。通过红外成像观测该区域，可识别并统计渗水发生的时间规律。

图1 监测对象区域Fig.1 Monitoring areas

1．1．2监测方法与数据获取 红外成像设备架设于奉先寺南部力士雕像腿部附近(图2)，对图1区域进行长期等间隔红外成像图像的探测，采样间隔为30 min。设备的测温灵敏度为0.1 ℃，生成的红外图像分辨率为640×480。本研究使用2017年全年的红外成像图像进行分析，经过对无效、重复的图像进行筛选后，共获取有效图像13 900张。

图2 红外成像监测设备Fig.2 Infrared imager

1．2 MATLAB语言进行图像处理

1．2．1技术特征和优势 长期监测获取的图像数据量极其庞大，而仅选取小部分图像进行观察和分析，无法得出统计规律或趋势；同时，传统的用于分析数值数据(如温湿度)的统计工具也不能应用在图像数据中。因此，若要利用图像来长期跟踪石窟保存环境以及病害的发展规律，必须同时解决上述问题。MATLAB语言让用户可以借助编程的方法进行图像处理。它将图像信息表示为多维数组，使得用户可以通过编写程序，将各种数学、统计方法直接运用在图像数据中。例如对图中特定点数值的提取，借助信号处理实现图像滤波，利用机器学习对图像内容的识别、分割等。这不仅实现了大量图像数据的同时处理，也极大地提高了分析的复杂度、深度和自由度。MATLAB作为前沿数学与计算机技术的平台，还为未来文物监测数据的挖掘和利用提供了很大的研究空间。

1．2．2处理方法 本研究全部数据分析工作均借助MATLAB语言实现。主要流程包括：将缺失、模糊、重复等图像查找并清除；建立图像色彩和温度值的映射，将图像转化为表示温度值的数组；根据图像拍摄时间建立数据索引；根据研究内容需要编写程序，包括提取特定坐标的数据列，利用统计工具进行分析，数据可视化等。

2 结果与分析

2．1 日变化特征的时序分析

本节选取了奉先寺红外图像中最具代表性的两组来展现晴天和雨天的日变化的特征，并提取图像中特定区域的数据，绘制时间序列曲线进行分析。分别选取并编号了若干区域(图3)：头像左上部分区域(A)及其相邻区域(B)，岩体砌补区域(C，包括岩体、砌补砖体和灌浆材料3个子区域)，用于观察各区域温度的差异；两个裂隙渗水点(D、E)、表面径流区域(F)及其周围岩体，用于通过温差来观察渗水的情况。

图3 特定区域的编号与选取位置示意图Fig.3 Numbers and locations of specific areas

为确保数据的代表性，后续分析中提取各区域温度值时，均取图3中矩形区域内全部像素点的平均值。同时，引入气候环境相同的万佛洞气象站的环境温度tenv和降雨量P做数据对比。

2．1．1晴天日变化 图4～5分别为2017年4月11～12日(晴)各区域的温度曲线和红外图像。夜间至日出之前，头像A区域的表面温度低于B区域0.5 ℃左右(后文分别用tA、tB表示)；灌浆材料、砌补砖体也略低于岩体(0.2～1 ℃不等)。日出后，阳光直接照射头像及岩体表面，各区域快速升温，于7∶00时达到同一水平； 8∶00时，可观察到图5a上各处表面温度出现显著差异。其中，tA明显高于tB，达到了二者之间当日的最大正温差1.8 ℃(记为Δt+)；随后，阳光直射逐渐消失，tA、tB之差缓慢减小，并于18∶00时持平。在这期间，图像中央岩体砌补区域C中，灌浆材料温度高于周围的砌补砖体和岩体0.1～1 ℃不等；接下来，夜间的环境温度持续下降，tA逐渐低于tB，并于次日5∶00达到最大负温差1.2 ℃(记为Δt-)。此间，C区域内的灌浆材料温度略低于砌补砖体和岩体(0.1～0.8 ℃不等)。

图4 2017年4月11～12日各区域表面温度变化曲线Fig.4 Variation of surface temperatures in different areas from 11 to 12 April，2017

图5 2017年4月11～12日不同时刻的红外图像Fig.5 Infrared images at different moments from 11 to 12 April，2017

相似的日变化出现在整个监测周期内的绝大多数晴天当中。可以看出，区域A与其相邻区域B经历着不同的日温度变化，致使日间tA>tB，而夜间tA

据此，结合奉先寺的保护修复情况，以及表面热平衡方程[15]

Rn=H+G+LE

(1)

和热通量等式[16]

(2)

可对A、B之间温差的成因予以阐释：日间，A、B表面吸收的太阳净辐射基本相同(RnA=RnB)；因A、B属同种材质，理化性能相似，故近似认为两者潜热相同(LEA=LEB)；tA与环境温度tenv的差异比tB稍大，因此A向大气耗散的感热更多(HA>HB)；综上由式1可推得，A、B表面所吸收的热通量GAtB，说明相比于B，A表面吸收的热通量更多地贮存了下来，而非继续向岩体深层传播。参考式2可进一步推断，这是由于灌浆材料以及裂隙内空气的平均热导率γ低于头像本体，阻隔了A向更深层传热GA↓的路径(GA↓

对于岩体砌补区域C，因为灌浆材料、砌补砖体和岩体具有不同的热传导性能，使得三者的热平衡状态有差异，进而导致了红外图像上表面温度的略微差异[7]。

2．1．2雨天日变化 图6～7为2017年6月4～6日降雨期间各区域的温度曲线以及红外图像。 在降雨最初的20 h内，图像上观测不到明显异常，渗水尚未形成，D、E、F以及岩体的表面温度相近(后文分别用tDE、tF、tr表示)；6月4日20∶00时，D、E处开始观测到明显的渗水，此时tDE明显低于tF、tr1.5 ℃；7 h后，图像上F处开始观测到自岩体顶部垂直向下的表面径流，同时D、E处渗水波及的区域也变广，此时tF降低至与tDE相同水平，并在随后约30 h内持续低于tr0.7～1.8 ℃；6月6日7∶00时(此时降雨已结束12 h)，天气转晴，阳光开始直射岩体表面，各区域温度骤升。随后径流停止，F区域表面水迹开始消退，tF逐渐回归至tr相同水平。但区域D、E依旧持续渗水，tDE与tr温差保持在2 ℃左右，最终本次渗水持续了近14 d。

图6 2017年6月4～6日各区域温度变化曲线Fig.6 Variation of surface temperatures in different areas from 4 to 6 June， 2017

图7 2017年6月4～6日不同时刻的红外图像Fig.7 Infrared images at different moments from 4 to 6 June， 2017

当降雨量充足时，即可观测到相似的水系活动模式：降雨持续一段时间后，先由D、E处产生渗水，而后F处出现径流；径流大多会在短时间内结束，并随着表面残余水分的蒸发，逐渐淡出图像；D、E处的渗水会持续相当长的时间。图1中其他泛白的岩体区域未观察到显著的渗水或径流。

关注降雨时期的A、B区域可以发现，渗水并未波及头像，其热平衡状态也未在雨期出现过异常。因此排除了头像处热导率差异与渗水的关联。

表1～2统计了2017年内红外成像监测区域出现渗水的情况，并通过相关分析探究了影响渗水和径流的因素。可以看出，渗水次数和出现延迟的间隔存在季节差异：冬春季渗水次数少，延迟较大(降雨后15～22 h)；夏秋季渗水次数多，延迟较小(降雨后3～9 h)；渗水与径流的持续时长之间显著正相关(0.81)，两者分别与降雨持续时长(0.74、0.67)和降雨总量(0.51、0.52)正相关。说明降雨总量越大，持续越久，渗水和径流就越明显；渗水延迟间隔、最大小时降雨量，与渗水、径流持续时长之间无显著相关关系。

表1 2017年红外监测区域的渗水统计Table 1 Statistics of water seepage in the infrared monitoring areas in 2017

表2 渗水统计变量的相关系数Table 2 Correlation coefficient of statistical variables of water seepage

2．2 阿难头像温差的统计分析

2．2．1阿难头像温差与热应力的关系 由2．1节可知，阿难头像A、B交界面处的温差具有日周期变化的特征。这种变化会带来差异性的膨胀和收缩，使得交界面产生同样周期变化的热应力[17]。这种应力的变化使得头像裂隙处存在风化破坏的风险。A、B交界面温差在日间达到最大值Δt+，夜间温度梯度反向，达到负最大值Δt-。参考热应力模型(式3)，可将该过程中应力的变化近似等效为块体A固定于B之中，A升温或降温Δt=|Δt+-Δt-|后的应力变化(假设两者在该过程中均不发生形变)[18-19]：

(3)

式中，σt为岩体沿x、y、z方向的热应力；τt为沿xy、yz、xz面的剪应力；E为弹性模量；α为膨胀系数；Δt为A温度的变化量；μ为泊松比。可以看出，式中温度的变化量Δt与岩体的热应力σt存在线性关系，故通过分析A、B每日最大温差Δt+、Δt-可以间接观测和评估温差可能引起的风化。

2．2．2温差的年变化差异 利用MATLAB语言编程对2017年内每日A、B间的最大温差进行极值的识别和提取，最终获得全年每日的Δt+、Δt-数据，并绘制其时序变化图(图8)。

图8 全年A、B区域间每日正负最大温差变化图Fig.8 Variation of the maximum daily temperature differences between Region A and B

可以看出，全年之内Δt+普遍大于Δt-；春冬季节Δt+和Δt-相对高于夏秋季节，且Δt+>4 ℃以及Δt-<-1.3 ℃的极大温差均出现在春冬季节，其中4月Δt+一度达到5.8 ℃；图像的时序变化规律不显著，不能观察到诸如季节性周期等周期性变化的存在；对比Δt+、Δt-与日降雨量P发现，降雨时期的温差明显小于非降雨时期。图中灰色矩形为监测设备数据缺失的时段。

图9为最大温差的频率直方图，带宽为0.3 ℃。可以看出Δt+呈现出双峰分布，分别出现在2.05 ℃和0.55 ℃附近，这两个峰值意味着在不同天气条件下，头像表面是否有阳光直射，使得温差分布出现了系统性差异。有73 d的Δt+集中在1.6～2.5 ℃，102 d集中在0.1～1.0 ℃，另有39 d高于2.8 ℃；Δt-则呈现集中的单峰分布，说明全年夜间温差的形成不存在系统性差异，183 d Δt-都集中在-0.7～-0.1 ℃。可以认为，相较于多云和阴雨天气，有阳光直射的晴天，头像的温差应力更显著。

图9 全年A、B区域间每日正负最大温差的频率直方图Fig.9 Histogram of the maximum daily temperature differences between Region A and B

2．2．3温差的时间段分布差异 图10～11分别展示了最大温差Δt在每日不同时段的散点图及其频率直方图，其中散点图通过不同图例进行了季节划分。

图10 不同时段最大温差的散点图Fig.10 Scatter plot of the maximum temperature differences at different time periods

图11 不同时段最大温差的频率直方图Fig.11 Histogram of the maximum temperature differences at different time periods

结合以上两张图可以看出，Δt+存在两种模式：1)一部分Δt+集中出现在8∶30～10∶30(161 d)，并且Δt+>1.6 ℃的大温差全部出现在该时段内，说明该模式为晴天的阳光直射，是温差应力的主要来源；2)其余Δt+介于 0～1 ℃之间，平缓分布于8∶30～16∶30，为多云或阴雨天的温差应力。

Δt-则大多分布在夜间1∶00～6∶00，以4∶00～5∶00最多(78 d)，说明最大负温差多出现在日出之前；由图10得知，温差的时间段分布没有明显的季节差异(冬春季小部分Δt-出现在了19∶00～24∶00，而夏秋季无此现象，具体原因有待进一步研究)。

2．3 红外图像的整体统计分析

本节借助MATLAB图像处理技术，对全年的奉先寺13 900张红外图像进行逐像素的垂直叠加，计算平均值和去趋势标准差两个统计量，并分别绘制与原始红外图像空间位置一一对应的矩阵的热图，从而实现对表面温度空间层面的统计分析。

图12中，由于水将部分太阳辐射吸收，并以潜热LE的形式蒸发耗散，使得渗水的区域的平均温度低于其他区域，并且温度越低，说明渗水的累积时长或岩体的潮湿时间越长。可以看到，D、E区域的温度平均值低于周围岩体约1 ℃，并且渗水源头处在图中展现出明显的黑点(低于周围岩体约2 ℃)；岩体仅F处观察到明显的径流痕迹，说明目前该区域受表面径流影响最严重。因此，该统计方法实现对不同区域渗水情况进行可视化的半定量估计，并准确寻找渗水点的位置。A、B、C区域之间的温度均值没有明显的差异，故无法通过该方法来观测温差应力。

图12 表面温度的平均值矩阵热图Fig.12 Heat map of the average value of surface temperatures

图13的绘制对原始数据进行了如下处理：先利用高通滤波的方法去除数据的季节性趋势，仅保留日变化周期和噪声，再计算标准差，旨在刻画表面温度的短期波动幅度[20]，以反映全年的温差风化情况。可以看出，A处的温度波动显著高于B；C区域内不同材质之间温度的波动也有差异：灌浆材料>砌补砖体>岩体。因此，该方法可以用于判断温差应力存在的位置，比较和评估不同区域温差风化的程度。D、E、F的温度波动与周围岩体十分相近，故不宜通过该方法观察渗水。

图13 表面温度的去趋势标准差矩阵热图Fig.13 Heat map of the detrended standard deviation of surface temperatures

3 讨 论

在尝试将红外成像监测和MATLAB图像处理技术结合应用于石窟寺预防性保护的过程中，发现该方法可以很好地描绘整个监测时段内石窟寺岩体的温差应力变化和渗水的发生情况：可以在观察红外图像，发现温度异常的区域后，进行各种针对特定点温度值的数学计算(如温差、温度梯度等)，并通过编程将对单个图像的处理复现在所有图像数据中，达到发掘温度应力变化规律的目的；可以快速定位岩体表面水的活跃区域、时段以及程度，统计渗水的发生规律。这避免了信息总量过于庞大而导致的冗余数据堆砌和分析困难，也为未来进一步深入挖掘图像中有价值信息提供了技术支持。

将实时监测和实际病害的发展相关联，一直是预防性保护工作的难点。本研究间接通过监测和计算日最大温差来表征温差应力的变化程度。由于石刻和灌浆材料的材质和形制独特，环境的热交换及太阳辐射也不断变化，使得阿难头像内部温度场的实际情况十分复杂，因此该方法的准确度还有待提高。若要相对准确地估计温差和风化程度之间的关系，需要进一步借助实验或热力学仿真模型进行研究。

单一的红外图像仅反映区域之间的温度差异，而在大量红外图像监测与MATLAB图像处理结合下，可以实现对不同材料热传导性能差异的表征，因此建议将该方法推广至其他复合材质文物或经保护材料干预的文物的预防性保护工作中，帮助进行保存情况的监测和保护修复效果的评估。

4 结 论

奉先寺阿难造像头部区域的裂隙处存在显著的温差应力循环，会引起温差风化。温差的成因是灌浆材料及裂隙内的空气带来的岩体内部热传导性能的差异，最终形成了交界面陡峭的温度梯度，且该温差应力存在日循环特征。

阿难头部温差应力的变化可概括为：晴天上午(约8∶30～10∶30)阳光直射带来的强烈的日温差(1.6～2.5 ℃)，这是温差应力的主要来源；多云或阴雨天时，正午前后较温和的日间温差(0.1～1.0 ℃)；夜间，温度梯度反向，并于日出之前出现反向最大温差(-0.7～-0.1 ℃)。因此，考虑到温差应力的来源，为了有效减缓其带来的破坏，可通过适当遮挡的手段来避免晴天上午阳光对阿难造像头部的短期直接照射，以显著降低日间的最大温差。

卢舍那大佛南侧岩体D、E两处渗水点长期活跃，表面径流则集中在F区域，其他区域未见明显径流。时间分布上，渗水出现于降雨后3～22 h，持续时间较长(4～33 d不等)；表面径流略晚于渗水出现，而持续时长较短。由于F区域径流的影响范围比D、E更大，且D、E作为裂隙渗水点，渗水关系更复杂隐蔽，因此可优先将治理的重点集中在该区域的短时表面径流中。相关性分析发现，渗水和径流的持续时长受降雨总量和降雨时长的直接影响。