曹鵬輝，呂書強，侯妙樂，趙林毅，汪萬福

(1. 北京建筑大學(xué)測繪與城市空間信息學(xué)院，北京 100044； 2. 北京市建筑遺產(chǎn)精細(xì)重構(gòu)與健康監(jiān)測重點實驗室，北京 100044；3. 敦煌研究院保護(hù)研究所，甘肅酒泉 736200； 4. 國家古代壁畫與土遺址保護(hù)工程技術(shù)研究中心，甘肅酒泉 736200)

0 引 言

中國古代壁畫存世數(shù)量巨大、色彩瑰麗、內(nèi)容豐富，具有極高的藝術(shù)和研究價值，是中國文化遺產(chǎn)的重要組成部分。然而，由于存世年代久遠(yuǎn)，自然環(huán)境以及人為等因素導(dǎo)致壁畫出現(xiàn)各種病害。由于壁畫的不可再生性，所以如何對壁畫進(jìn)行無損留存與保護(hù)成為現(xiàn)在研究的熱點之一。而壁畫的線狀特征決定著畫面的整體空間的劃分、畫面主題定位及物體間的位置關(guān)系等各元素。它們構(gòu)成了每一幅畫面的雛形，也是將畫面的生命力和感染力傾注在畫面上的奠基者[1]，在壁畫的保護(hù)與修復(fù)中占據(jù)著重要地位。

近年來高光譜技術(shù)快速發(fā)展[2-4]，因其分析手段無損、波段數(shù)目多、波段覆蓋范圍廣、波段分布連續(xù)、光譜分辨率高，以及“圖譜合一”等特點，可以更好地表達(dá)彩繪文物的完整信息，所以被廣泛應(yīng)用在彩繪文物的數(shù)字化保護(hù)中。例如在信息挖掘方面，郭新蕾等[5]利用主成分分析技術(shù)提取了古畫頭冠周圍涂抹痕跡。Herens等[6]利用高光譜技術(shù)對凱斯·凡·東根的一幅畫作進(jìn)行顏料分析，發(fā)現(xiàn)了一個隱藏起來的女性肖像。吳太夏等[7]通過光譜分析，成功地提取了中國黑墨水的特征，在墓葬的兩個重疊木片的隱藏面上發(fā)現(xiàn)了60多個漢字。史寧昌等[8]利用短波紅外波段，對部分故宮館藏書畫文物進(jìn)行分析，在增強印章的文字信息、發(fā)現(xiàn)書畫的涂改涂抹痕跡和提取底稿線等方面得到了一定進(jìn)展。在顏料識別方面，王功明等[9]采用光譜表示顏料的顏色，然后借助獨立成分分析，設(shè)計了一種分離混合顏料的方法，采用蒙賽爾色卡的光譜信息進(jìn)行模擬實驗，恢復(fù)出原始顏料的種類及比例。Rohani等[10]利用高光譜稀疏模型對埃及發(fā)掘的一幅古畫進(jìn)行了顏料解混，并識別了不同顏料的種類。目前，常用于增強線狀特征方法有通過對高光譜數(shù)據(jù)進(jìn)行主成分分析(Principle Component Analysis，PCA)或最小噪聲分離(Minimum Noise Fraction，MNF)變換后找尋線狀特征清晰的特征波段[11-12]，利用模糊聚類算法對數(shù)字影像的模糊邊緣進(jìn)行增強[13]，利用Laplacian算子進(jìn)行邊緣增強[14]等方法。但是，利用MNF及PCA方法在選擇特征波段的時候需要人為干預(yù)，影響增強結(jié)果。其他方法增強的對象限于真彩色數(shù)字影像，其波段有限，無法完全留存壁畫信息，增強結(jié)果會有一定的信息缺失。

小波變換也被頻繁應(yīng)用于圖像增強、圖像融合、人臉識別等領(lǐng)域[15-19]，但小波變換處理的圖像大多為可見光范圍內(nèi)的數(shù)字影像，無法有效利用近紅外波段線狀特征信息豐富的優(yōu)點。

利用高光譜影像覆蓋波長寬、信息豐富的優(yōu)點，結(jié)合MNF變換與Haar小波變換，提出了一種壁畫線狀特征的增強方法，以青海省瞿曇寺的壁畫局部為研究實例進(jìn)行了線狀特征增強，并與PCA線狀特征增強進(jìn)行了效果對比。

1 研究對象與數(shù)據(jù)預(yù)處理

1.1 研究對象

瞿曇寺位于青海省海東市樂都區(qū)南21 km處，始建于明洪武二十五年(1392年)，1982年被評為第二批全國重點文物保護(hù)單位。該寺現(xiàn)存明清兩代壁畫1 338 m2[20]，為明清兩代宮廷畫師所作，是瞿曇寺的藝術(shù)瑰寶，具有極高的藝術(shù)價值。目前，由于存世久遠(yuǎn)，部分壁畫存在明顯的空鼓、起甲、顏料層脫落、褪色等病害，對其保護(hù)修復(fù)亟待進(jìn)行。因此，利用高光譜MNF變換與Haar小波變換結(jié)合，選擇瞿曇寺的東西回廊部分壁畫進(jìn)行線狀特征的增強處理，期待為壁畫的保護(hù)修復(fù)提供更豐富更直觀的參考信息。

1.2 數(shù)據(jù)獲取與預(yù)處理

1.2.1數(shù)據(jù)獲取 所用儀器為VNIR400H型高光譜成像儀，光譜范圍為400～1 000 nm，光譜分辨率為2.8 nm，相機為內(nèi)置掃描方式，波段數(shù)為1 040個，成像畫幅為1 392×1 000個像元。拍攝時間為2018年8月28日下午，現(xiàn)場無直射自然光，環(huán)境較暗，采用人工光源。儀器距壁畫109 cm，光圈為4.0，曝光時間為80 ms。主要實驗區(qū)域的正射影像如圖1所示。紅色區(qū)域為研究區(qū)域。

圖1 西回廊十五區(qū)正射影像圖Fig.1 Orthophoto of the 15th district of west corridor

1.2.2預(yù)處理 高光譜成像儀采集的原始數(shù)據(jù)是目標(biāo)物反射的輻射亮度，其中包含壁畫本身的輻射信息和噪聲。另外考慮到地面高光譜成像的物距只有1 m左右，大氣輻射傳輸對輻射亮度的影響可以忽略。因此，只需要對高光譜原始數(shù)據(jù)進(jìn)行反射率重建，如式(1)。

(1)

式中，Rref為反射率重建后數(shù)據(jù)；Rdate為高光譜原始影像數(shù)據(jù)；Rwhite為同等環(huán)境下白板原始影像數(shù)據(jù)；Rdark為暗電流噪聲數(shù)據(jù)。

白板使用反射率為99%的標(biāo)準(zhǔn)白板，為減少燈光和儀器拍攝角度對白板數(shù)據(jù)的影響，采集白板的環(huán)境與采集壁畫時的環(huán)境和儀器參數(shù)保持一致。暗電流數(shù)據(jù)在關(guān)閉光源，蓋上鏡頭蓋時進(jìn)行數(shù)據(jù)獲取，反映了儀器在沒有輻射能量輸入時儀器的噪聲數(shù)據(jù)。

2 研究方法

針對壁畫線狀特征經(jīng)常會褪色而難以辨認(rèn)的問題，綜合利用MNF變換能分離噪聲、小波變換能分離低頻和高頻信號的特點，設(shè)計了壁畫線狀特征增強的流程(圖2)。對反射率重建后的高光譜影像進(jìn)行MNF變換，計算MNF變換后的前m個波段的平均梯度(Average Gradient，AG)，選擇AG最高的波段作為最優(yōu)波段，該波段線狀特征最豐富。同時對MNF逆變換的影像進(jìn)行RGB組合并轉(zhuǎn)化為灰度圖像。利用Haar小波對最優(yōu)波段和MNF逆變換合成的灰度圖像進(jìn)行小波分解，將兩者獲得的低頻信號相融合，并與MNF逆變換的灰度影像分解的高頻信號進(jìn)行Haar小波逆變換，得到線狀特征明顯的增強影像。

圖2 線狀特征增強流程圖Fig.2 Flow chart of the algorithm proposed

2.1 MNF變換

MNF變換[21]是一種高光譜數(shù)據(jù)常用的降維方法，本質(zhì)是進(jìn)行兩次主成分變換。第一次變換分離出原始數(shù)據(jù)中的噪聲，如式(2)，使得變換后的噪聲數(shù)據(jù)具有最小的方差且各波段不相關(guān)。第二次變換是對噪聲數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)主成分變換，如式(3)。

X=S+N

(2)

式中，X為反射率重建后數(shù)據(jù)；S為無噪聲數(shù)據(jù)；N為噪聲數(shù)據(jù)。

Y=U×X

(3)

式中，Y為經(jīng)過MNF變換后數(shù)據(jù)；U為N的協(xié)方差的逆矩陣與X的協(xié)方差矩陣相乘的特征向量。

原始高光譜數(shù)據(jù)通過MNF變換以后，獲得各波段不相關(guān)且信噪比由高到低排列的數(shù)據(jù)。

2.2 平均梯度法

平均梯度的是描述圖像在垂直和水平方向上灰度的平均變化率[22]。由于線狀特征的灰度值與其他物質(zhì)的灰度值差距較大，所以若影像中線狀特征越豐富、清晰，該影像的平均梯度就越大。利用這種方法對降維后的高光譜數(shù)據(jù)進(jìn)行最優(yōu)波段選擇，選擇出線狀特征豐富并且清晰的波段，如式(4)。

2.3 Haar小波變換

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

(5)

(6)

利用Haar小波對圖像進(jìn)行一級分解可以得到一個低頻信號(LL)和3個不同方向的高頻信號，分別為水平方向分量(HL)，豎直方向分量(LH)和對角線分量(HH)。

圖像分解后，主要信息集中在低頻信號上，而噪聲和細(xì)節(jié)信息集中在高頻信號，所以只需要對低頻信號進(jìn)行增強，對高頻信號進(jìn)行弱化即可以達(dá)到圖像增強的目的。

2.4 線狀特征增強方法

使用經(jīng)過反射率重建的數(shù)據(jù)進(jìn)行MNF變換，選擇信噪比較高的前m個波段，分別利用式(4)計算各波段的平均梯度，選擇平均梯度值最大的波段為最優(yōu)波段。利用前m個波段信噪比高的特點進(jìn)行MNF逆變換獲得低噪聲的高光譜數(shù)據(jù)。對其進(jìn)行RGB組合，利用式(7)將RGB影像轉(zhuǎn)化為灰度影像。

I=R×0.299+G×0.587+B×0.114

(7)

式中，I為合成灰度影像；R為紅波段影像；G為綠波段影像；B為藍(lán)波段影像。

由于最優(yōu)波段影像與MNF逆變換后合成的灰度影像的灰度區(qū)間不同，所以對最優(yōu)波段影像和MNF逆變換后合成的灰度影像進(jìn)行歸一化，保證其灰度區(qū)間在[0，1]內(nèi)，如式(8)。

(8)

式中，sn為歸一化后的影像；s為待歸一化影像；smin為待歸一化影像s中灰度的最小值；smax為待歸一化影像s中灰度的最大值。

利用Haar小波對進(jìn)行歸一化處理后的灰度影像和最優(yōu)波段影像進(jìn)行分解。雖然最優(yōu)波段影像的低頻信號可以清晰地反映出絕大多數(shù)線狀特征信息并且沒有噪聲的影響，但單一波段無法完全保留壁畫的線狀特征，需要與MNF逆變換的灰度影像的低頻信號相融合，從而確保線狀特征的完整度。MNF逆變換影像比原始影像的噪聲更低。在經(jīng)過小波分解后，高頻信號中的噪聲信息降低，突出了圖像細(xì)節(jié)信息。這不但達(dá)到削弱高頻信號的目的，也保留了圖像的細(xì)節(jié)信息。利用增強后的低頻信號與削弱后的高頻信號后進(jìn)行重構(gòu)，獲得增強后影像。在低頻信號融合的過程中，應(yīng)按式(9)對兩幅影像的低頻信號進(jìn)行加權(quán)。具體流程圖如圖3。

圖3 融合示意圖Fig.3 Fusion flow chart

L=λL1+(1-λ)L2

(9)

式中，L1、L2分別為最優(yōu)波段影像與灰度影像利用Haar小波分解的低頻部分；L為低頻部分融合結(jié)果；λ為權(quán)重系數(shù)。

3 結(jié)果與分析

以瞿曇寺西回廊15區(qū)的部分壁畫為例，進(jìn)行實驗并分析結(jié)果。在ENVI 5.1上對原始影像進(jìn)行MNF變換、MNF逆變換、RGB組合與灰度圖像轉(zhuǎn)換。其余所有處理均使用Matlab 2013 a實現(xiàn)。

3.1 實驗過程

3.1.1MNF變換與MNF逆變換 對輻射校正后的影像進(jìn)行MNF變換獲得信噪比由高到低排列的各波段不相關(guān)的數(shù)據(jù)，其特征值隨波段增加而減小，選擇特征值高的前10個波段作為MNF變換的結(jié)果。利用前10波段進(jìn)行MNF逆變換，去除高光譜圖像中的噪聲。

3.1.2基于平均梯度的最優(yōu)波段選擇 對MNF變換后的10個波段進(jìn)行歸一化后分別計算平均梯度，結(jié)果如圖4。MNF變換后的第二波段的平均梯度明顯高于其他波段，利用目視法進(jìn)行檢驗，證明該波段為線狀特征最為豐富、清晰的波段，所以選擇該波段為最優(yōu)波段。

圖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)為藍(lán)色波段進(jìn)行RGB組合，并利用式(7)將其轉(zhuǎn)換為灰度影像(圖5a)。

利用Haar小波對歸一化后最優(yōu)波段和MNF逆變換后的灰度圖像進(jìn)行分解、重構(gòu)，得到線狀特征的增強影像(圖5b)。

圖5 結(jié)果對比Fig.5 Results of comparison

為同時利用MNF逆變換后合成的灰度影像與最優(yōu)波段影像的低頻信號，在低頻信號融合時權(quán)重均為0.5。

3.2 結(jié)果分析

如圖5，對比增強后灰度影像與合成灰度影像發(fā)現(xiàn)區(qū)域1和區(qū)域2中的線狀特征明顯增強，區(qū)域3中被顏料覆蓋的線狀特征通過本研究方法增強后也清晰可見。將數(shù)據(jù)歸一化后，截取區(qū)域1～3利用Sobel算子進(jìn)行邊緣提取后，計算其平均梯度，結(jié)果如表1?？梢娫鰪姾笥跋衿骄荻染笥谠蓟叶扔跋?，進(jìn)一步證明線狀特征獲得增強。

表1 平均梯度對比結(jié)果Table 1 Results of average gradient

為了證明該方法的普適性，選擇瞿曇寺東回廊三區(qū)和西回廊十五區(qū)的另一部分壁畫進(jìn)行實驗，實驗結(jié)果如圖6～7所示。

圖6 東回廊三區(qū)實驗結(jié)果Fig.6 Experimental results of the third district of east corridor

圖7 西回廊十五區(qū)另一部分實驗結(jié)果Fig.7 Experimental results of another area in the 15th district of west corridor

通過對比真彩色影像與增強后灰度影像發(fā)現(xiàn)東回廊三區(qū)人物的腰帶和衣服細(xì)節(jié)等區(qū)域的線狀特征明顯增強。西回廊十五區(qū)的另一部分除增強線狀特征外，還使被顏料覆蓋的文字信息清晰可見。證明該方法在壁畫的線狀特征增強和隱含信息的挖掘的方面具有良好的普適性。

3.3 線狀特征增強對比

文獻(xiàn)[11]中，提出基于PCA變換的線狀特征增強方法。利用該方法對圖1所示區(qū)域高光譜數(shù)據(jù)進(jìn)行PCA，選擇最優(yōu)波段，然后與圖5中所得結(jié)果進(jìn)行對比，如圖8。另外，計算局部區(qū)域的平均梯度與結(jié)構(gòu)相似性(Structural similarity index measurement，SSIM)進(jìn)行對比。對數(shù)據(jù)進(jìn)行歸一化后截取區(qū)域1～3，分別計算基于PCA增強、MNF-Haar小波結(jié)合增強兩種線狀特征增強方法與降噪后合成灰度影像的結(jié)構(gòu)相似性，結(jié)果如表2。利用Sobel算子對歸一化后的區(qū)域1～3進(jìn)行邊緣檢測后，分別計算兩種方法不同區(qū)域的平均梯度，結(jié)果如表2。區(qū)域1、區(qū)域2由于裂縫的影響，使用文獻(xiàn)[11]中方法增強后裂縫過于明顯，所以導(dǎo)致計算平均梯度值高于MNF-Haar小波結(jié)合增強方法。對于無裂縫的區(qū)域3，MNF-Haar小波結(jié)合增強方法的平均梯度較好，證明該方法可以更好地實現(xiàn)線狀特征的增強，且增強后圖像與原圖具有較高的結(jié)構(gòu)相似性。

圖8 方法對比Fig.8 Methods contrast

表2 不同方法的結(jié)構(gòu)相似性與平均梯度Table 2 Structural similarity index measurement and average gradient of different methods

4 結(jié) 論

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