龍門石窟奉先寺的紅外成像監測與分析
——基于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作為裂隙滲水點，滲水關系更復雜隱蔽，因此可優先將治理的重點集中在該區域的短時表面徑流中。相關性分析發現，滲水和徑流的持續時長受降雨總量和降雨時長的直接影響。