基于空譜特征聯合的煙熏壁畫線條增強方法

毛錦程，呂書強，侯妙樂，汪萬福

1. 航天規劃設計集團有限公司，北京 100162；

2. 北京建筑大學 測繪與城市空間信息學院，北京 100044；

3. 建筑遺產精細重構與健康監測北京市重點實驗室，北京 100044；

4. 敦煌研究院保護研究所，敦煌 736200；

5. 國家古代壁畫與土遺址保護工程技術研究中心，敦煌 736200

1 引 言

古壁畫色彩豐富、風格多樣，具有極高的文化價值和研究價值。但由于壁畫長期受煙熏污染影響，部分線條信息已較難辨認，迫切需要保護修復。線條作為古壁畫中的核心元素，體現了畫家的初始想法以及壁畫圖案的繪畫結構和藝術思想（汪萬福等，2015）。此外，線條有助于文物修復人員更加清晰地感受匠師的技法和情感，對壁畫的傳承與保護具有極其重要的作用，線條的提取也是壁畫修復至關重要的一環（袁小樓，2018）。

近年來，國內外學者利用高光譜技術對壁畫線條信息的提取進行了多種嘗試，主要可以分為三類方法。第一類是利用壁畫的空間信息進行特征提取。侯妙樂等（2014）利用高光譜成像技術結合主成分分析（principal component analysis，PCA）等方法提取并增強了壁畫的底稿信息；郭新蕾等（2017）利用PCA 結合光譜分析挖掘古畫的隱藏信息和涂抹信息；Pan 等（2017）通過PCA、波段選擇、密度分割，增強了墓葬壁畫中的片狀圖案的輪廓信息。第二類是利用壁畫的光譜信息對線條的特征信息進行增強。Han 等（2015）選擇近紅外波段結合光譜角匹配算法分析底稿的顏料，與可見光影像圖像融合后得到壁畫的視覺增強影像；閆青（2019）提出了稀疏非負矩陣欠近似的方法對書畫、墓葬壁畫的待修復區域進行線稿提取。第三類是利用深度學習及其他計算機輔助技術實現線條的特征信息提取。呂書強等（2022）通過改進U-Net 網絡較好地提取了壁畫顏料層脫落病害區域；Cao 等（2020）通過建立網絡模型進行特征提取，將多種損失函數與網絡模型相結合獲取優化后的古代壁畫影像；孫振榮（2020）通過在邊緣檢測算法中融入雙向級聯網絡實現彩繪文物的線條提取，針對噪聲較多或損毀相對嚴重的文物引入了生成對抗網絡，實現了線條的虛擬修復；Pan 等（2018）利用卷積神經網絡對壁畫線稿信息進行預測，并通過檢測特定特征區域來保留壁畫的藝術風格。

但是，壁畫高光譜影像多數呈非線性，基于整幅影像進行線性降維未能充分挖掘壁畫線條的隱含信息，難以顧及壁畫的光譜信息、局部特征差異性等因素。此外，高光譜影像波段數目多，其行、列相關性存在較大差異也會導致部分線條信息難以提取。深度學習等方法所提取的線條雖然精度較高，但需要大量的訓練數據集及復雜的調整參數過程，限制了其在壁畫保護工程中的實用性。

因此，本文利用壁畫高光譜影像波段數多、圖譜合一的優勢，通過分析影像中的光譜信息和空間信息，提取出壁畫高光譜影像中潛在的線條信息。將核主成分分析（kernel principal component analysis，KPCA）算法引入分區降維中，結合二維主成分分析(two dimensional principal component analysis，2DPCA)與光譜特征分析對線條進行增強與提取，利用自適應伽馬校正的圖像增強算法實現其顏色空間的復原。因此，本文建立了一種基于空譜特征聯合的煙熏壁畫線條增強方法，實現壁畫的線條信息提取。

2 數據來源及預處理

1）數據來源

瞿曇寺建于明洪武二十五年（1392 年），地處青海省海東市樂都區瞿曇鎮，主體建筑為明代官式，1982 年被列入第二批全國重點文物保護單位（仲旭輝，2018）。隨著時間推移，壁畫受到自然或人為等多種因素的影響，病害較嚴重，部分壁畫的線條模糊不清，急需保護和修復。本研究選取瞿曇寺大黑天殿的部分壁畫，對壁畫中被煙熏覆蓋區域下的線條進行增強與提取，以期為壁畫的保護工作提供科學的參考數據。

2）數據預處理

地面高光譜成像儀采集得到的是目標物體的輻射亮度，主要包含兩部分：一是目標物本身的輻射信息；二是噪聲信息。由于在實際拍攝時距離較短，大氣輻射傳輸對輻射亮度的影響往往可以忽略不計，因此，僅需要對原始高光譜數據進行反射率校正。高光譜影像預處理包括反射率校正與影像去噪兩部分：

式中，R為反射率影像；RData為原始的高光譜影像數據，即物體本身的輻射亮度；RWhite為同等拍攝環境下的標準白板影像數據；RDark為暗電流數據，即拍攝環境中的噪聲信息。暗電流采集時，需關閉所有光源，蓋上鏡頭蓋。數據采集過程中使用的標準白板數據反射率為99%，環境與真實壁畫數據采集環境一致。

考慮到壁畫的原始高光譜數據有1040 個波段，數目較多。受儀器影響波長兩端含有大量噪聲信息，故去除噪聲信息較多的前后50 個波段。對裁剪后的940 個波段進行最小噪聲分離(minimum noise fraction，MNF)正逆變換。逆變換后的影像數據能在保留影像的空間信息和光譜信息的同時抑制噪聲信息，達到圖像去噪的目的。

3 研究方法

針對壁畫部分區域由煙熏覆蓋導致的線條信息模糊無法辨認的問題，本文提出一種基于空譜特征聯合的煙熏壁畫線條增強方法，將改進主成分分析算法與自適應伽馬校正算法相結合，對煙熏壁畫的高光譜影像中的線條進行增強，具體技術流程如圖1 所示。首先，用地面高光譜成像儀采集壁畫的原始高光譜數據并進行數據預處理，選取波段合成真彩色影像；利用多尺度模糊C 均值算法進行圖像聚類，將分類后數據與高光譜影像一一對應得到分區高光譜數據；對各分區數據進行KPCA 降維獲取前幾主成分，將各分區數據的主成分影像重構，得到壁畫整幅影像的前幾主成分數據。利用平均梯度分析選出分區KPCA 降維后的最優主成分影像。其次，為進一步分析壁畫每景影像的行、列相關性，對預處理后的高光譜影像進行2DPCA 獲取其行、列主成分，經圖像重構和平均梯度分析獲得2DPCA降維得到的最優主成分影像。最后，結合光譜特征分析與自適應伽馬校正算法得到包含壁畫色彩信息和線條增強信息的影像。

圖1 線條信息增強技術路線Fig.1 Schematic diagram of line information enhancement

3.1 多尺度模糊C 均值算法

傳統模糊C 均值（fuzzy C-means，FCM）算法通過對目標函數進行優化更新樣本數據的聚類中心和隸屬度，從而實現數據聚類（Ahmed 等，2002）。FCM 算法的目標函數可表示為

式中，U為隸屬矩陣；Ci為聚類中心；d2ij為第i個聚類中心與第j個樣本數據的歐氏距離；m為模糊加權指數；uij為隸屬度，是評估樣本屬于某個類別的程度，取值為[0,1]。

由于傳統的FCM 算法對噪聲敏感，常導致像素的空間信息缺失，存在過分割現象。在壁畫中，圖案中相鄰的像素通常屬于同一類別，因此，可利用超像素根據影像的局部空間信息，將其分割成均勻的區域。由于超像素圖像承載了圖像的空間信息，用超像素圖像中不同區域顏色像素的平均值替換原圖像中的各顏色像素，可以有效融合自適應局部空間信息和全局顏色特征。

首先，利用多尺度形態梯度重建可以選取正確的分割尺寸，提高FCM 算法的分類效果。多尺度形態梯度重建函數可定義為

式中，GMC是采用多尺度形態梯度重建函數對梯度圖像進行重構后所得到的重構圖像；A為原始影像；B為標記影像；Z為結構元素；r1為最小區域的尺寸大??；r2為最大區域的尺寸大??；GC為閉運算。由于形態梯度重建可以在保持目標形態輪廓細節的同時，消除了噪聲和無用的形態梯度細節，可有效減少影像過分割。因此，經過多尺度形態梯度重建方法處理后的影像，具有更準確的輪廓信息（Lei等，2018）。

用超像素區域各像素對應顏色層級的平均值替代原始圖像對應區域中各像素的顏色層級，減少不同顏色的數量，對具有自適應局部空間信息的影像執行FCM 聚類。

3.2 改進主成分分析算法

KPCA 算法是一種典型的非線性降維方法。根據多尺度模糊C均值算法所得到的聚類結果作為標簽數據，劃分預處理后的高光譜數據，得到分區高光譜數據。通過引入非線性映射函數將各分區高光譜數據映射到高維空間使其線性可分，再進行主成分分析降維，剔除方差較小的特征達到降維的目的（Schlkopf 等，1997）。高維特征空間中的主成分分析可表達為

式中，λi、wi分別為高維特征空間的特征值和特征向量；Φ(X)為非線性函數。由于在特征空間中直接計算點積，計算量較大，因此，引入核函數替代點積。計算核矩陣并將其歸一化，解得歸一化后矩陣的特征值與特征向量，并計算樣本在高維特征空間下的投影。

3.3 二維主成分分析

2DPCA 是一種基于圖像矩陣的特征提取算法。與主成分分析算法相比，2DPCA 直接利用圖像矩陣構建協方差矩陣，包含的結構信息更多，能更為準確地估計協方差矩陣，大大提高了計算效率（Li 等，2010）。設原始高光譜影像數據為X∈Rl×M×N，其中，l為波段數目，M和N為圖像的行、列數，將影像經過翻轉重塑為二維矩陣，轉換后的二維圖像矩陣分別為X∈和X∈RNl×M。以圖像矩陣X∈RMl×N為例，將其投影到一個N維的單位列向量上得到投影特征向量。圖像的協方差矩陣可定義為

式中，D為圖像樣本總數；為訓練樣本的均值影像；Ai為第i個波段的影像，是M×l的矩陣。

選擇協方差矩陣的前幾個特征值所對應的特征向量構成投影矩陣。將投影矩陣依次投影到前幾個特征向量上，得到對應的主成分，使用特征提取后的主成分和特征向量計算重構影像。同理，計算轉換后的另一個二維矩陣，采用平均加權的方法將行、列對應的兩幅主成分圖像進行融合，得到前幾個主成分影像。

3.4 光譜特征分析

頂點成分分析(vertex component analysis, VCA)是一種基于線性混合模型的端元提取算法，通過尋找單形體的頂點來提取端元（Nascimento 和Dias，2005）。單形體頂點對應各端元向量，迭代地將數據投影到已確定端元組成的子空間正交方向上，投影距離最大的點對應為新的端元，遍歷得到所有端元。在實際拍攝中，由于壁畫受不同顏料混合的復雜性、光照、環境及數據采集條件等多種因素的影響，像元中端元光譜通常會發生變化，產生光譜變異性，且難以兼顧像元中所有物質。因此，選擇非負約束最小二乘算法求解壁畫中各圖案區域的豐度（Das 和 Chakravortty，2021）。由于約束項為不等式，故通常情況下可轉化成一個最優化問題，得到拉格朗日表達式為

式中，A為端元矩陣；s為各端元在像元中所占比例；m∈Rl×1為某個像元矩陣。

3.5 加權分布自適應伽馬校正算法

在數據采集過程中，環境條件、光照不足或室內照明不均勻等多因素的影響，可能導致壁畫影像的亮度和對比度較低。經過圖像融合后，雖能提高整幅影像的對比度，但其邊緣亮度較低，導致部分線條無法分辨。因此，通過自適應伽馬校正函數來提高圖像較暗區域的亮度（Huang 等，2013）。

首先，利用二維直方圖生成圖像中的上下文和變量信息，使用高斯混合模型來補償圖像的灰度分布（Kim 和 Chung，2008）；通過伽馬函數來修改直方圖，使其圖像亮度保持均衡化。但由于修改后的直方圖可能會丟失一些統計信息，降低了增強效果。通過構造一個歸一化的伽馬函數來修改變換曲線，保留可用的統計直方圖。自適應伽馬校正公式如下：

式中，cdf 為累積分布函數；V為圖像的像素值；γ為自適應參數。

自適應伽馬校正算法可以避免圖像亮度的顯著衰減。利用加權分布函數能對統計直方圖進行略微修改，減少誤差。

4 結果與分析

選取瞿曇寺大黑天殿北墻部分壁畫為研究區域，對采集到的高光譜數據進行特征提取并分析結果。在ENVI5.3 上對壁畫原始影像進行MNF 正逆變換，其余所有處理均使用Matlab2016a 實現。

4.1 壁畫線條增強結果

選擇預處理后高光譜影像中的不同波段合成真彩色影像，壁畫原始影像中包含了不同圖案區域的信息，利用多尺度模糊C 均值算法對合成真彩色影像進行聚類，本次的加權指數、收斂條件和最大迭代次數分別設為2、1×10-5、50，窗口大小設為3×3。利用鄰域的局部信息得到的同質區域影像，根據對應的超像素圖像計算彩色圖像的直方圖，因為超像素圖像所劃分的區域數量遠遠小于原始影像中的像素數量。用一個區域內所有像素的平均值來代替該區域內的像素，可減少原始影像中不同顏色的數量。分割后的影像中包含了大量的小區域，且與合成真彩色影像顏色相似，分割效果較好。

由于融合多尺度分割后的分類結果可以更好地利用同質區域的空間信息，提取每個分割尺度下的同質區域的前幾主成分，得到不同分割尺度下的非線性低維特征。通過對前幾主成分影像進行平均梯度分析，結果如表1 所示。各主成分影像的信息熵相差較小。因此，選擇平均梯度最大的第五個主成分影像作為最優主成分影像，如圖2(a)所示。

表1 FCM+KPCA 主成分影像梯度分析Tab. 1 FCM + KPCA principal component image gradient analysis

圖2 影像對比Fig.2 Comparison of results

為減少高光譜影像的行、列相關性的影響，對預處理后影像進行旋轉變換后，進行二維主成分分析獲取前幾主成分影像。由于第六主成分之后均為噪聲，故對前五個主成分影像進行平均梯度計算，結果如表2 所示。第二主成分影像的平均梯度數值較高，因此，將其作為最優主成分影像，結果如圖2(b)所示。為了進一步利用高光譜影像中的光譜信息，利用VCA 對預處理后的高光譜影像進行端元提取。通過非負約束最小二乘算法得到各端元的豐度圖，其中，黑色圖案區域的線條信息更加豐富，結果如圖2(c)所示。

由圖2(c)可以看出，豐度圖中雖線條信息較為清晰，但與其背景的對比度較低，為了進一步增強線條信息，將2DPCA 和分區KPCA 的最優主成分影像歸一化后與線條豐度圖進行波段運算，得到的線狀特征增強影像，結果如圖2(e)所示。將增強影像與合成真彩色影像融合，通過加權分布自適應伽馬校正對融合影像進行色彩校正，結果如圖2(f)所示。對比圖2(e)、(d)，可以看到線狀特征增強影像中線條信息清晰，影像中人物內部的線條信息相比數字影像有較大程度的增強。融合后的影像包含了壁畫的顏色，增強了線條與其他顏色圖案的對比度，提高了壁畫的視覺效果。但影像邊緣區域增強效果略差，而經過顏色調整后，壁畫影像邊緣的亮度得到加強。

為定量分析各模塊的增強效果，將原始灰度影像與線狀特征增強影像、合成真彩色影像、線狀特征融合影像和最終影像進行了精度評價，結果如表3 所示。線狀特征增強影像的數值均高于原始灰度影像，最終影像與線狀特征融合影像的數值均高于合成真彩色影像，表明本方法中各模塊處理算法均對壁畫影像有一定的增強效果。

表3 不同模塊精度對比Tab. 3 Accuracy comparison of different modules

4.2 消融實驗

為了比較各部分影像對壁畫線條信息增強的影響，分別對FCM+KPCA、2DPCA 及光譜特征分析三部分得到的線條信息增強影像轉換到由色度（hue）、飽和度（saturation）、明度（value）組成的HSV 顏色空間得到融合圖像，結果如圖3 所示。FCM+KPCA 歸一化后的融合影像線條信息增強效果較弱，且與其背景的對比度較低。2DPCA 歸一化后的融合影像較合成真彩色影像、數字影像中的線條信息有一定的增強，但影像中人物、邊緣的深色區域顏色較深，該部分線條信息仍被遮蓋。線條豐度圖的融合影像中人物主體及邊緣的線條信息較為清晰，但因背景顏色較深，與線條的對比度較弱，視覺效果較差。

圖3 消融實驗對比結果Fig.3 Comparison results of ablation experiments

4.3 與現有方法對比

為了驗證本方法的有效性，選取瞿曇寺大黑天殿的部分壁畫，進行煙熏覆蓋下的線條增強與提取，與單尺度Retinex(single scale retinex，SSR)算法(李佳等，2020)、帶色彩恢復的多尺度Retinex(multi-scale retinex with color restoration，MSRCR)算法(Wei 和Li，2021) 和侯妙樂等（2014）方法（簡稱高光譜增強算法）進行對比，結果如圖4所示。

圖4 不同方法的影像對比結果Fig.4 Image comparison results of different methods

計算不同方法得到增強影像的平均梯度及邊緣強度進行定量分析，結果如表4 所示。圖像平均梯度和邊緣強度越大，線條信息增強效果越好。

表4 不同方法平均梯度及邊緣強度對比Tab. 4 Comparison of average gradient and edge intensity of different methods

由表3、表4 可知，SSR 算法、MSRCR 算法及高光譜增強算法得到的增強影像中平均梯度及邊緣強度的數值均低于本方法。此外，從圖4(b)可以看出，SSR 算法雖較合成真彩色影像的色彩有所增強，但其噪聲較大，且影像中脫落區域附近的線條信息增強效果不明顯；圖4 中MSRCR 算法和高光譜增強算法雖噪聲信息較少，但影像整體偏暗，影像部分區域的線條信息因圖案顏色較深或煙熏污染導致影像局部對比度較低，細節增強效果較差。高光譜增強算法的增強影像較原始灰度影像的平均梯度、邊緣強度分別提高24.44%、24.67%；SSR 算法、MSRCR 算法和本方法影像較合成真彩色影像的平均梯度分別提高 39.08%、4.85%、150.86%；SSR 算法、MSRCR 算法和本方法影像較合成真彩色影像的邊緣強度分別提高41.64%、6.65%、143.49%。因此，本方法較SSR 算法、MSRCR算法和高光譜增強算法的精度有較大提高。

為了更好地研究本方法的普適性，將瞿曇寺其他區域的壁畫按上述研究方法進行處理，與合成真彩色影像對比，結果如圖5 所示。經過本方法增強后，該區域可以更加直觀地辨認出被煙熏污染的壁畫影像，使肉眼較難辨認的模糊線條得以增強。因此，本研究提出的壁畫線條信息增強方法可以較好地實現不同壁畫的線條信息增強。

圖5 瞿曇寺大黑天殿部分區域壁畫影像對比Fig.5 Image comparison of some areas in Qutan Temple

5 結 論

針對壁畫線條信息被表面物質覆蓋辨認難的問題，提出兼顧高光譜影像空間和光譜信息的線條增強與提取的方法，實現了被煙熏污染下的壁畫線條信息增強與提取。通過多尺度模糊C 均值算法獲取標簽數據并劃分壁畫的高光譜影像，利用平均梯度分析選取分區KPCA 和2DPCA 的最優主成分影像。結合高光譜影像數據中的光譜信息利用光譜特征分析得到線條豐度圖。通過波段運算得到壁畫線條增強影像，將其進行HSV 融合得到線狀特征融合影像，利用加權分布自適應伽馬校正算法更好地復原壁畫色彩。以青海省瞿曇寺大黑天殿壁畫的高光譜影像為例進行了驗證。結果表明：與合成真彩色影像相比，影像線條增強效果明顯；與SSR 算法、MSRCR 算法及侯妙樂等（2014）方法相比，本方法能較好區分壁畫線條與其他圖案背景，增強線條與其他圖案的對比度；與數字影像相比，本次提出的方法能更好地復原壁畫色彩，使畫面更富有層次感，人物形象更加鮮活，提升壁畫的視覺效果。本研究能夠為壁畫的保護與修復工作提供借鑒與參考。