曹鵬輝,呂書強,侯妙樂,趙林毅,汪萬福
(1. 北京建筑大學測繪與城市空間信息學院,北京 100044; 2. 北京市建筑遺產精細重構與健康監測重點實驗室,北京 100044;3. 敦煌研究院保護研究所,甘肅酒泉 736200; 4. 國家古代壁畫與土遺址保護工程技術研究中心,甘肅酒泉 736200)
中國古代壁畫存世數量巨大、色彩瑰麗、內容豐富,具有極高的藝術和研究價值,是中國文化遺產的重要組成部分。然而,由于存世年代久遠,自然環境以及人為等因素導致壁畫出現各種病害。由于壁畫的不可再生性,所以如何對壁畫進行無損留存與保護成為現在研究的熱點之一。而壁畫的線狀特征決定著畫面的整體空間的劃分、畫面主題定位及物體間的位置關系等各元素。它們構成了每一幅畫面的雛形,也是將畫面的生命力和感染力傾注在畫面上的奠基者[1],在壁畫的保護與修復中占據著重要地位。
近年來高光譜技術快速發展[2-4],因其分析手段無損、波段數目多、波段覆蓋范圍廣、波段分布連續、光譜分辨率高,以及“圖譜合一”等特點,可以更好地表達彩繪文物的完整信息,所以被廣泛應用在彩繪文物的數字化保護中。例如在信息挖掘方面,郭新蕾等[5]利用主成分分析技術提取了古畫頭冠周圍涂抹痕跡。Herens等[6]利用高光譜技術對凱斯·凡·東根的一幅畫作進行顏料分析,發現了一個隱藏起來的女性肖像。吳太夏等[7]通過光譜分析,成功地提取了中國黑墨水的特征,在墓葬的兩個重疊木片的隱藏面上發現了60多個漢字。史寧昌等[8]利用短波紅外波段,對部分故宮館藏書畫文物進行分析,在增強印章的文字信息、發現書畫的涂改涂抹痕跡和提取底稿線等方面得到了一定進展。在顏料識別方面,王功明等[9]采用光譜表示顏料的顏色,然后借助獨立成分分析,設計了一種分離混合顏料的方法,采用蒙賽爾色卡的光譜信息進行模擬實驗,恢復出原始顏料的種類及比例。Rohani等[10]利用高光譜稀疏模型對埃及發掘的一幅古畫進行了顏料解混,并識別了不同顏料的種類。目前,常用于增強線狀特征方法有通過對高光譜數據進行主成分分析(Principle Component Analysis,PCA)或最小噪聲分離(Minimum Noise Fraction,MNF)變換后找尋線狀特征清晰的特征波段[11-12],利用模糊聚類算法對數字影像的模糊邊緣進行增強[13],利用Laplacian算子進行邊緣增強[14]等方法。但是,利用MNF及PCA方法在選擇特征波段的時候需要人為干預,影響增強結果。其他方法增強的對象限于真彩色數字影像,其波段有限,無法完全留存壁畫信息,增強結果會有一定的信息缺失。
小波變換也被頻繁應用于圖像增強、圖像融合、人臉識別等領域[15-19],但小波變換處理的圖像大多為可見光范圍內的數字影像,無法有效利用近紅外波段線狀特征信息豐富的優點。
利用高光譜影像覆蓋波長寬、信息豐富的優點,結合MNF變換與Haar小波變換,提出了一種壁畫線狀特征的增強方法,以青海省瞿曇寺的壁畫局部為研究實例進行了線狀特征增強,并與PCA線狀特征增強進行了效果對比。
瞿曇寺位于青海省海東市樂都區南21 km處,始建于明洪武二十五年(1392年),1982年被評為第二批全國重點文物保護單位。該寺現存明清兩代壁畫1 338 m2[20],為明清兩代宮廷畫師所作,是瞿曇寺的藝術瑰寶,具有極高的藝術價值。目前,由于存世久遠,部分壁畫存在明顯的空鼓、起甲、顏料層脫落、褪色等病害,對其保護修復亟待進行。因此,利用高光譜MNF變換與Haar小波變換結合,選擇瞿曇寺的東西回廊部分壁畫進行線狀特征的增強處理,期待為壁畫的保護修復提供更豐富更直觀的參考信息。
1.2.1數據獲取 所用儀器為VNIR400H型高光譜成像儀,光譜范圍為400~1 000 nm,光譜分辨率為2.8 nm,相機為內置掃描方式,波段數為1 040個,成像畫幅為1 392×1 000個像元。拍攝時間為2018年8月28日下午,現場無直射自然光,環境較暗,采用人工光源。儀器距壁畫109 cm,光圈為4.0,曝光時間為80 ms。主要實驗區域的正射影像如圖1所示。紅色區域為研究區域。

圖1 西回廊十五區正射影像圖Fig.1 Orthophoto of the 15th district of west corridor
1.2.2預處理 高光譜成像儀采集的原始數據是目標物反射的輻射亮度,其中包含壁畫本身的輻射信息和噪聲。另外考慮到地面高光譜成像的物距只有1 m左右,大氣輻射傳輸對輻射亮度的影響可以忽略。因此,只需要對高光譜原始數據進行反射率重建,如式(1)。
(1)
式中,Rref為反射率重建后數據;Rdate為高光譜原始影像數據;Rwhite為同等環境下白板原始影像數據;Rdark為暗電流噪聲數據。
白板使用反射率為99%的標準白板,為減少燈光和儀器拍攝角度對白板數據的影響,采集白板的環境與采集壁畫時的環境和儀器參數保持一致。暗電流數據在關閉光源,蓋上鏡頭蓋時進行數據獲取,反映了儀器在沒有輻射能量輸入時儀器的噪聲數據。
針對壁畫線狀特征經常會褪色而難以辨認的問題,綜合利用MNF變換能分離噪聲、小波變換能分離低頻和高頻信號的特點,設計了壁畫線狀特征增強的流程(圖2)。對反射率重建后的高光譜影像進行MNF變換,計算MNF變換后的前m個波段的平均梯度(Average Gradient,AG),選擇AG最高的波段作為最優波段,該波段線狀特征最豐富。同時對MNF逆變換的影像進行RGB組合并轉化為灰度圖像。利用Haar小波對最優波段和MNF逆變換合成的灰度圖像進行小波分解,將兩者獲得的低頻信號相融合,并與MNF逆變換的灰度影像分解的高頻信號進行Haar小波逆變換,得到線狀特征明顯的增強影像。

圖2 線狀特征增強流程圖Fig.2 Flow chart of the algorithm proposed
MNF變換[21]是一種高光譜數據常用的降維方法,本質是進行兩次主成分變換。第一次變換分離出原始數據中的噪聲,如式(2),使得變換后的噪聲數據具有最小的方差且各波段不相關。第二次變換是對噪聲數據進行標準主成分變換,如式(3)。
X=S+N
(2)
式中,X為反射率重建后數據;S為無噪聲數據;N為噪聲數據。
Y=U×X
(3)
式中,Y為經過MNF變換后數據;U為N的協方差的逆矩陣與X的協方差矩陣相乘的特征向量。
原始高光譜數據通過MNF變換以后,獲得各波段不相關且信噪比由高到低排列的數據。
平均梯度的是描述圖像在垂直和水平方向上灰度的平均變化率[22]。由于線狀特征的灰度值與其他物質的灰度值差距較大,所以若影像中線狀特征越豐富、清晰,該影像的平均梯度就越大。利用這種方法對降維后的高光譜數據進行最優波段選擇,選擇出線狀特征豐富并且清晰的波段,如式(4)。


Haar小波是小波分析中具有緊支撐的正交小波函數。Haar小波函數定義如式(5),其尺度函數如式(6)。

(5)

(6)
利用Haar小波對圖像進行一級分解可以得到一個低頻信號(LL)和3個不同方向的高頻信號,分別為水平方向分量(HL),豎直方向分量(LH)和對角線分量(HH)。
圖像分解后,主要信息集中在低頻信號上,而噪聲和細節信息集中在高頻信號,所以只需要對低頻信號進行增強,對高頻信號進行弱化即可以達到圖像增強的目的。
使用經過反射率重建的數據進行MNF變換,選擇信噪比較高的前m個波段,分別利用式(4)計算各波段的平均梯度,選擇平均梯度值最大的波段為最優波段。利用前m個波段信噪比高的特點進行MNF逆變換獲得低噪聲的高光譜數據。對其進行RGB組合,利用式(7)將RGB影像轉化為灰度影像。
I=R×0.299+G×0.587+B×0.114
(7)
式中,I為合成灰度影像;R為紅波段影像;G為綠波段影像;B為藍波段影像。
由于最優波段影像與MNF逆變換后合成的灰度影像的灰度區間不同,所以對最優波段影像和MNF逆變換后合成的灰度影像進行歸一化,保證其灰度區間在[0,1]內,如式(8)。
(8)
式中,sn為歸一化后的影像;s為待歸一化影像;smin為待歸一化影像s中灰度的最小值;smax為待歸一化影像s中灰度的最大值。
利用Haar小波對進行歸一化處理后的灰度影像和最優波段影像進行分解。雖然最優波段影像的低頻信號可以清晰地反映出絕大多數線狀特征信息并且沒有噪聲的影響,但單一波段無法完全保留壁畫的線狀特征,需要與MNF逆變換的灰度影像的低頻信號相融合,從而確保線狀特征的完整度。MNF逆變換影像比原始影像的噪聲更低。在經過小波分解后,高頻信號中的噪聲信息降低,突出了圖像細節信息。這不但達到削弱高頻信號的目的,也保留了圖像的細節信息。利用增強后的低頻信號與削弱后的高頻信號后進行重構,獲得增強后影像。在低頻信號融合的過程中,應按式(9)對兩幅影像的低頻信號進行加權。具體流程圖如圖3。

圖3 融合示意圖Fig.3 Fusion flow chart
L=λL1+(1-λ)L2
(9)
式中,L1、L2分別為最優波段影像與灰度影像利用Haar小波分解的低頻部分;L為低頻部分融合結果;λ為權重系數。
以瞿曇寺西回廊15區的部分壁畫為例,進行實驗并分析結果。在ENVI 5.1上對原始影像進行MNF變換、MNF逆變換、RGB組合與灰度圖像轉換。其余所有處理均使用Matlab 2013 a實現。
3.1.1MNF變換與MNF逆變換 對輻射校正后的影像進行MNF變換獲得信噪比由高到低排列的各波段不相關的數據,其特征值隨波段增加而減小,選擇特征值高的前10個波段作為MNF變換的結果。利用前10波段進行MNF逆變換,去除高光譜圖像中的噪聲。
3.1.2基于平均梯度的最優波段選擇 對MNF變換后的10個波段進行歸一化后分別計算平均梯度,結果如圖4。MNF變換后的第二波段的平均梯度明顯高于其他波段,利用目視法進行檢驗,證明該波段為線狀特征最為豐富、清晰的波段,所以選擇該波段為最優波段。

圖4 MNF變換前10波段的平均梯度Fig.4 Average gradient of the first 10 bands
3.1.3基于Haar小波變換的線狀特征增強 選擇MNF逆變換后的band 383(640 nm)作為紅光波段,band 241(550 nm)為綠色波段,band 94(460 nm)為藍色波段進行RGB組合,并利用式(7)將其轉換為灰度影像(圖5a)。
利用Haar小波對歸一化后最優波段和MNF逆變換后的灰度圖像進行分解、重構,得到線狀特征的增強影像(圖5b)。

圖5 結果對比Fig.5 Results of comparison
為同時利用MNF逆變換后合成的灰度影像與最優波段影像的低頻信號,在低頻信號融合時權重均為0.5。
如圖5,對比增強后灰度影像與合成灰度影像發現區域1和區域2中的線狀特征明顯增強,區域3中被顏料覆蓋的線狀特征通過本研究方法增強后也清晰可見。將數據歸一化后,截取區域1~3利用Sobel算子進行邊緣提取后,計算其平均梯度,結果如表1。可見增強后影像平均梯度均大于原始灰度影像,進一步證明線狀特征獲得增強。

表1 平均梯度對比結果Table 1 Results of average gradient
為了證明該方法的普適性,選擇瞿曇寺東回廊三區和西回廊十五區的另一部分壁畫進行實驗,實驗結果如圖6~7所示。

圖6 東回廊三區實驗結果Fig.6 Experimental results of the third district of east corridor

圖7 西回廊十五區另一部分實驗結果Fig.7 Experimental results of another area in the 15th district of west corridor
通過對比真彩色影像與增強后灰度影像發現東回廊三區人物的腰帶和衣服細節等區域的線狀特征明顯增強。西回廊十五區的另一部分除增強線狀特征外,還使被顏料覆蓋的文字信息清晰可見。證明該方法在壁畫的線狀特征增強和隱含信息的挖掘的方面具有良好的普適性。
文獻[11]中,提出基于PCA變換的線狀特征增強方法。利用該方法對圖1所示區域高光譜數據進行PCA,選擇最優波段,然后與圖5中所得結果進行對比,如圖8。另外,計算局部區域的平均梯度與結構相似性(Structural similarity index measurement,SSIM)進行對比。對數據進行歸一化后截取區域1~3,分別計算基于PCA增強、MNF-Haar小波結合增強兩種線狀特征增強方法與降噪后合成灰度影像的結構相似性,結果如表2。利用Sobel算子對歸一化后的區域1~3進行邊緣檢測后,分別計算兩種方法不同區域的平均梯度,結果如表2。區域1、區域2由于裂縫的影響,使用文獻[11]中方法增強后裂縫過于明顯,所以導致計算平均梯度值高于MNF-Haar小波結合增強方法。對于無裂縫的區域3,MNF-Haar小波結合增強方法的平均梯度較好,證明該方法可以更好地實現線狀特征的增強,且增強后圖像與原圖具有較高的結構相似性。

圖8 方法對比Fig.8 Methods contrast

表2 不同方法的結構相似性與平均梯度Table 2 Structural similarity index measurement and average gradient of different methods
利用高光譜MNF變換與Haar小波結合,提出了一種壁畫線狀特征的增強方法。首先,對高光譜影像進行MNF變換,選擇信息量集中的前幾個波段進行MNF逆變換實現降噪處理。其次,對重構后的影像選擇真彩色波段變換為灰度圖像,對其進行Haar小波分解。然后,對MNF變換后的影像,利用最大平均梯度法進行最優波段選擇,挑選出線狀特征信息豐富的波段進行Haar小波分解,利用分解后的低頻信號與灰度圖像分解后的低頻信號相融合,高頻信號使用重構后的降噪圖像代替。最后對優化組合的低頻和高頻信號進行Haar小波逆變換得到結果圖像。以青海省瞿曇寺壁畫局部為研究實例,利用提出的方法進行了線狀特征增強,經過與原始灰度影像、主成分分析線狀特征增強方法進行對比,結果表明MNF變換與Haar小波結合的壁畫線狀特征增強方法具有較好的效果,驗證了提出方法的有效性。研究表明,利用高光譜與小波變換相結合的方法可以有效增強壁畫模糊的線狀特征,實現隱含信息提取,能夠為壁畫修復提供線狀特征參考,使得修復有據可依。