光譜色差與主分量特征融合的石窟表面風化智能評估

曹赤鵬，王慧琴*，王 可，王 展，張 剛，馬 濤

（1.西安建筑科技大學 信息與控制工程學院，陜西 西安710055；2.陜西省文物保護研究院，陜西 西安710075）

1 引 言

在大型石質文物表面風化等級評估方法研究中，不同風化類型及風化程度之間的表征問題一直是研究的熱點。傳統的風化程度檢測主要依靠人工采集風化表面巖體的物理化學參數，通過計算各參數與風化類型及程度之間的相關性對不同風化類型及風化程度進行表征［1］。采集過程工作量巨大，并且采樣點的規劃極易受人為主觀因素影響，難以滿足大型石質文物表面復雜風化區域的精細化表征。

光譜成像技術［2-4］迅速發展，廣泛應用于壁畫顏料識別［5］、油畫顏料信息提取［6］、石質文物表面風化程度檢測［7］等工作。在石質文物表面風化程度檢測方法研究中，國內外眾多學者將多光譜成像技術和光譜分析技術相結合，通過采集石質文物表面風化離散光譜通道的光譜信息和空間信息對不同風化類型及風化程度區域進行稀疏表征。反射光譜不僅可以反映物質的狀態特征，還反映了表面顏色、表面粒度等物理特征［8］。柴勃隆［9］等將多光譜成像技術與石窟考古研究相結合，利用壁畫表面不同物質反射光譜特征的不同，成功對莫高窟壁畫在可見光范圍內肉眼無法辨識的彩色區域進行了修復。周霄［10］等利用光譜成像技術采集了云岡石窟表面風化區域的光譜反射率數據，利用光譜反射率幅值的差異對風化程度進行表征，有效區分出強風化和弱風化區域。但反射光譜曲線特征變化較為平緩，光譜反射率數據存在較大冗余，不能實現數據的快速處理，并且對于石窟表面復雜風化區域，光譜反射率幅值特征較為單一，難以達到精細化表征，需要探索一種有效的表征方法。

顏色是對可見光特征的表征，可見光的波長決定了人類眼睛所感知的顏色。對于石窟而言，不同風化類型及風化程度區域的顏色由其表面的反射光譜決定。不同風化區域的顏色差異真實地反映了石窟表面不同風化類型及風化程度之間的差別。主成分分析方法［11］是對多波段數據的線性變換，將原始光譜數據轉換到一個數據維度較少的新空間，通過計算使光譜數據的差異達到最大，對于增強石窟表面風化光譜信息含量、隔離噪聲及減少數據維數非常有效［12］。基于此，本文提出了一種基于光譜色差與主分量特征融合的石窟表面風化智能評估方法。通過重建石窟表面風化光譜反射率，計算多光譜成像數據中每個像素點與基準點的色差；并對石窟表面原始多光譜成像數據進行主成分分析，提取主分量特征數據，將光譜色差數據與主分量特征數據進行融合，建立不同風化類型及風化程度表征數據庫；利用隨機森林算法（Random Forest，RF）對石窟表面復雜風化區域每個像素點進行智能評估，為大型石質文物表面復雜風化區域精細化評估提供可行性方案。

2 石窟表面不同風化區域光譜色差特性分析

石窟表面在長期風化作用下其表面物質、表面顏色發生改變。表面巖體中的鈣、鎂等氫氧化物在空氣中二氧化碳的長期作用下，使得巖體表面顏色變為白色，風化程度越深顏色表征越為明顯［13］。

以黃色砂礫巖巖體石窟為例，利用多光譜成像系統對石窟表面風化區域進行多光譜成像數據采集，如圖1所示為石窟表面不同風化類型及風化程度純凈風化區域的多光譜成像平均反射光譜數據。隨鹽析風化程度的加深，其所對應的反射光譜的幅值越高；隨積塵風化的加深，其所對應的反射光譜幅值越低。反射光譜不僅反映石窟表面巖體物質成分的改變，同時還反映不同風化類型及風化程度區域間顏色的差異。為了準確分析石窟表面風化類型及風化程度的改變對顏色變化的影響，分別對不同的風化類型及風化程度區域的反射光譜進行CIE1976L*a*b*、CIE1931XYZ以及RGB顏色空間的線性變換，通過各參數值的變化來反映石窟表面不同風化類型及程度區域間表面顏色的變化情況。

圖1 石窟表面不同風化類型及風化程度反射光譜表征Fig.1 Reflection spectrum characterization of different weathering types and degrees on the surface of grotto

表1為石窟表面風化區域的CIE1976L*a*b*顏色空間、CIE1931XYZ顏色空間以及R，G，B各參數隨不同風化類型及風化程度的變化情況，其中 變 量Wstrong，Wweak，Wslightly，Wdust，Wbenchmark分 別表示強鹽析風化，弱鹽析風化，微鹽析風化，積塵風化，基準點。鹽析風化類型中風化程度1強于風化程度2，積塵風化類型中風化程度2強于風化程度1。通過對比表1中各參數值可以發現，隨鹽析風化程度的加深，石窟表面巖體中鈣、鎂等的氫氧化物溶液與空氣中的二氧化碳反應生成白色的碳酸鈣或碳酸鎂等的結晶，導致表面巖體顏色變白，使得巖體表面亮度逐漸增加，三種顏色空間的各參數值呈現逐步上升趨勢。對于積塵風化區域，由于石窟表面酥化，產生的碎屑摻雜著塵土覆蓋于石窟表面，使得積塵風化區域的明亮度變低，相對于鹽析風化區域其RGB各通道值最低。不同的風化類型及風化程度區域間其表面顏色的差異可以用光譜色差來表征。以基準點L*a*b*值為標準，計算不同風化類型及風化程度間的光譜色差值，用ΔE表示。在國際顏色標準CIE1976L*a*b*的評價體系中，色差的判斷標準如表2所示。不同風化類型及風化程度區域與基準間的光譜色差值如表3所示。

表1 不同風化類型及風化程度區域的L*a*b、XYZ、RGB變化表Tab.1 Changes of L*a*b，XYZ and RGB in different weathering types and weathering degree areas

表2 色差判斷標準T ab.2 Color difference judgment standard

通過表3可以看出，由于石窟表面不同風化程度區域表面顏色的不同，導致不同風化程度區域與基準間的色差存在較大差異。鹽析風化作用越強與基準間的色差值越高，強鹽析風化與弱鹽析風化之間色差為6.38、弱鹽析風化與微鹽析風化之間色差為9.69、微鹽析風化與積塵風化之間色差為17.88，在色差判斷標準中其差異均處于很大或者非常大水平。根據以上分析，光譜色差可以對石窟表面不同風化類型及風化程度進行有效表征，為光譜色差在石窟表面風化智能評估提供理論支撐與數據基礎。

表3 不同風化類型及風化程度與基準的色差T ab.3 Color difference between different weathering types and weathering degrees and reference points

3 相關原理

3.1 基于多光譜成像數據的石窟表面光譜反射率重建

石窟表面風化多光譜成像數據由多波段組成，具有較強的離散特性，使得光譜色差計算時存在較大誤差。光譜反射率重建可以消除多光譜成像數據的離散特性，重建后的光譜信息可以充分表征石窟不同風化區域表面物質及顏色差異，精確的重現不同風化類型及風化程度區域的顏色信息。

在前期光譜重建理論研究與實驗基礎之上，本文采用多核支持向量回歸方法，重建石窟表面多光譜成像的光譜反射率數據。與偽逆法、多項式法、單核支持向量回歸等傳統重建方法相比，多核支持向量回歸將局部核函數和全局核函數相結合，通過局部核函數對測試點較近函數值的強學習能力與全局核函數對測試點較遠函數值的強泛化能力，能夠有效解決傳統重建模型泛化能力差、重建精度低的問題，并在光譜重建顏色復原實驗中獲得較高的精度［14］。本文采用局部核函數Poly核函數與全局核函數Cauchy核函數相結合的方法，引入到支持向量回歸中，利用多核支持向量回歸對石窟表面不同風化類型及風化程度的光譜反射率進行重建。回歸函數為：

其中：R(x i)為波長為λi處的不同風化類型及風化程度的光譜反射率重建結果；K cp(x i，y)表示重建光譜反射率的核函數；i表示需要重建的光譜響應值個數，取值范圍為(i=1，2，???，n)；x i表示第i個測試樣本光譜響應數據；y表示訓練樣本光譜反射率數據；b表示測量誤差。通過多核核函數對重建算法進行優化，多核核函數構造如式（2）所示：

其中：x表示石窟表面待重建區域不同風化類型及風化程度的多光譜相機響應輸入值；σ表示柯西核函數的寬度，作為各支持向量之間的關聯系數；d表示多項式核函數階數。則石窟表面風化多光譜成像數據中每個像素點的光譜反射率重建結果可以表示為(R1(x)，R2(x)，???，Ri(x))T。

3.2 石窟表面風化像素級光譜色差計算

利用支持向量回歸對石窟表面不同風化類型及風化程度區域每個像素點的光譜反射率進行重建；將重建后的光譜反射率，從光譜空間轉換到CIE1976L*a*b*顏色空間，計算每個像素點的L*，a*，b*參數值。在CIE1976L*a*b*顏色空間中，顏色通過三種對比的刺激值進行表示，分別為黑-白、綠-紅、藍-黃。L*表示黑/白刺激值，當L*為0時，表示黑色，當L*為100時表示白色；a*表示綠/紅刺激值，當a*為正值時表示紅的程度，當a*為負值時表示綠的程度；b*表示藍/黃刺激值，當b*為正值時表示黃的程度，當b*為負值時表示藍的程度。CIE1976L*a*b*顏色空間的各參數值是由CIE1931XYZ顏色空間通過線性變換而來。石窟表面風化多光譜成像數據計算X、Y、Z三刺激值的計算過程如式（3）～式（7）所示：

其中：X，Y，Z分別是CIE1931XYZ顏色空間的三刺激值；X n，Y n，Z n表示標準白板的三刺激值，并通過Y n值為100進行標準化處理。通過顏色空間的線性變換得到了石窟表面不同風化類型及風化程度區域內每個像素點的L*，a*，b*參數值。并以基準點為計算基準對石窟表面不同風化類型及風化程度區域的色差進行計算，用ΔE i表示，則計算公式如下：

其中：L*0，a*0，b*0表示基準點的CIE1976L*a*b*值；L*i，a*i，b*i表示石窟表面不同風化類型及風化程度區域每個像素點的CIE1976L*a*b*值。

3.3 石窟表面風化多光譜成像數據主成分分析

對石窟表面多光譜成像數據進行降維，可以有效降低多光譜成像數據的冗余，提高智能評估方法的計算效率。光譜成像數據降維的常用方法主要分為非線性和線性降維，其中非線性降維以核主成分分析（Kernel Principal Component Analysis，KPCA）［15］、核獨立成分分析（Kernel In?dependent Component Analysis，KICA）［16］為主，通過核函數完成非線性映射，以實現對非線性數據降維，同時保留了原始光譜數據的主要信息［17］；線性降維以主成分分析（Principal Compo?nent Analysis，PCA）［18］、獨立成分分析（Indepen?dent Component Analysis，ICA）［19］為主，通過線性空間變換將高維數據空間投影到低維空間，即保留原始數據的大部分方差信息又消除了原始數據成分間的相關分量［20］。由于獨立成分分析方法適用于數據的解混，且石窟表面風化光譜數據變化較為平緩，不同風化類型及風化程度間的反射光譜特性曲線線性可分，所以本文采用主成分分析方法對石窟表面多光譜成像數據進行線性降維處理。首先以統計計算的方式獲得石窟表面風化多光譜成像數據中各波段間的相關性，從中提取特征向量和特征值，再以二者為基礎進行正交變換。主成分分析方法對石窟表面風化多光譜成像數據處理的關鍵在于計算轉換矩陣。通過轉換矩陣將石窟表面風化多光譜成像數據從高維空間映射到低維空間，用較少的獨立的主分量圖像全面綜合反映有用信息［21］。特征值的大小也反映著數據的變化程度，特征值越大說明該像素點的反射光譜響應數據變化程度越大，對于分類的評估結果影響程度越高，在多光譜成像數據中越重要。將石窟表面多光譜成像數據按相同像素點位置的相機響應值作為向量，以此建立石窟表面風化多光譜成像數據矩陣。原始的石窟表面風化多光譜成像數據由16個通道組成，每個通道的圖片大小為1 276×1 028，將16個通道的光譜成像數據表示為16個大小為1×1 311 728的向量集合X=[x1，x2???x i]T，每個向量x i有1 311 728個反射光譜數據構成，x i=[x i1，x i2，???，x in]，并 對 每 一 列 光 譜 數 據 進 行 零 均值化處理，即每一列的反射光譜數據減去該行數據的均值，然后計算x ij的協方差矩陣Cov，由協方差矩陣求其特征值λn及其特征向量，計算公式為：

其中：將特征值矩陣記為Λ，由特征值所對應的特征向量組成的正交矩陣記為An，將特征值從大到小排列，將其中對應的前n個特征向量組成轉換矩陣T，將主成分集表示為Y，則主成分提取公式為：

4 實驗及結果分析

4.1 實驗部分

利用海洋光學公司的SpectroCamVIS型號多光譜成像系統對陜西省延安市清涼山萬佛寺萬佛窟表面風化巖體進行多光譜圖像數據采集。對于復雜風化區域，波長采集范圍400～940 nm，波段間隔20～60 nm，共采集16個波段，分辨率大小為1 276×1 028 pixel，如圖2所示。首先對清涼山萬佛寺萬佛窟表面風化多光譜圖像數據進行處理，將大小為1 276×1 028×3的主分量特征數據與大小為1 276×1 028×1的光譜色差特征數據進行重組，重組后的特征數據大小為1 276×1 028×4；然后按像素點位置提取該像素點的光譜色差數據和主分量特征數據，建立大小為1 311 728×4的石窟表面風化待測樣本數據庫；最后對石窟表面不同風化類型及風化程度純凈風化區域進行標記，按標記位置提取光譜色差與主分量特征融合數據，并將該像素點所對應的不同風化類型及風化程度標簽標記于特征數據矩陣最后一列，建立數據大小為34 860×5的訓練及測試樣本數據庫。

圖2 石窟表面風化區域各波段多光譜成像數據Fig.2 Multi-spectral imaging data of weathering area on the surface of grotto

為驗證光譜色差與主分量融合數據對不同風化類型及風化程度表征的適用性及有效性，利用隨機森林算法、K-近鄰算法（K Near Neigh?bor，KNN）、BP神經網絡（Back Propagation Neu?ral Networks，BP）、徑向基神經網絡（Radial Ba?sis Function，RBF）分別對光譜色差數據、主分量特征數據、光譜色差數據與主分量融合數據以及16通道原始反射光譜數據進行對比實驗。為降低不同風化類型及風化程度數據間的不平衡對模型訓練的影響，將融合表征數據庫的10%用于預測，剩余部分采用5折交叉驗證方法，輪流選取每種風化類型及風化程度數據的20%用于測試，80%用于訓練，循環5次。

通過四種實驗的預測準確率、Kappa系數及整體風化程度預測評估結果與該像素區域實際風化情況對比，對評估結果進行精度評價。石窟表面風化智能評估方法技術框圖如圖3所示。

圖3 石窟表面風化智能評估方法技術框圖Fig.3 Technical block diagram of intelligent evaluation method for grotto surface weathering

4.1.1 石窟表面風化光譜色差計算

為準確還原石窟表面不同風化區域的顏色信息，減少離散多光譜成像數據對不同風化類型及風化程度間光譜色差計算的影響，利用支持向量回歸方法對石窟表面多光譜成像數據進行光譜反射率重建。不同風化類型及風化程度區域的光譜反射率重建后的光譜特性曲線如圖4所示。

圖4 石窟表面重建光譜特性曲線Fig.4 Reconstructed spectral characteristic curve of grot?to surface

對石窟表面不同風化區域的重建光譜特征數據進行顏色空間的線性變換，計算每個像素點的L*，a*，b*參數值；以基準點L*a*b*值為標準計算石窟表面不同風化類型及風化程度區域間的光譜色差。

在國際色差評價標準下，強鹽析風化與弱鹽析風化區域光譜色差值為5.83，顏色差異較為明顯；強鹽析風化與微鹽析風化、積塵風化之間的光譜色差值均大于12，顏色差異非常大；弱鹽析風化與微鹽析風化之間的光譜色差值為6.85，顏色差異很大；弱鹽析風化與積塵風化區域光譜色差值均大于12，顏色差異非常大；微鹽析風化與積塵風化之間光譜色差值為8.33，顏色差異很大。不同風化類型及風化程度間的光譜色差計算結果如表4所示。利用光譜色差特性可以對石窟表面不同風化類型及風化程度區域進行有效表征。

表4 不同風化類型及風化程度間的色差Tab.4 Color difference between different weathering types and degrees

4.1.2 主分量數據特征提取

利用主成分分析方法對原始多光譜成像數據的16個波段圖像進行降維處理。圖5表示主成分分析方法對石窟表面風化多光譜成像數據進行處理的各波段特征值分布情況，橫軸代表波段值，縱軸代表特征值。

圖5 石窟表面風化多光譜圖像主成分分析結果Fig.5 Principal component analysis results of multi-spec?tral image of cave surface weathering

通過特征值數量隨波段的變化結果可以看出，石窟表面風化多光譜成像的數據信息主要集中在第一波段和第二、三波段，其余波段所含的信息較少。圖6（a）表示石窟表面風化多光譜成像數據經過主成分分析第一主成分分量（PC1）圖像，圖6（b）表示第二主成分分量（PC2）圖像，圖6（c）表示第三主成分分量（PC3）圖像，圖6（d）表示前三主成分分量合成的偽彩色圖像。通過對比可以看出第一主成分具有的信息最多，第三主成分具有的信息較少，噪聲占比上升。

圖6 主成分分析結果圖Fig.6 Principal component analysis results image color image image

前3個主分量貢獻率及累計貢獻率如表5所示。第一主成分分量具有的圖像信息十分豐富，其信息貢獻率達到總信息量的95.41%，第二主成分分量和第三主成分分量貢獻率較低分別為3.30%和0.53%，但包含了較多的石窟表面的巖性信息，有利于對不同風化類型及風化程度的區分識別。

表5 前3個主分量貢獻率及累計貢獻率Tab.5 Contribution rate and cumulative contribution rate of the first three principal components

前三個主成分分量累計貢獻率為99.24%。利用主成分分析方法對多光譜成像數據進行處理，即降低了數據維度又保留了石窟表面風化反射光譜的主要信息。

4.2 石窟表面風化整體評估結果分析

將石窟表面風化光譜色差數據、主分量特征數據、反射率幅值表征數據、光譜色差與主分量特征融合表征數據分為訓練、測試以及預測訓練樣本。通過RF、KNN、BP、RBF四種算法分別對不同風化類型及風化程度的石窟表面風化特征數據進行學習，并保存訓練模型，通過預測準確率及Kappa系數對四種風化程度表征方法的有效性進行對比分析。表6記錄了四種評估方法對石窟表面不同風化類型及風化程度特征數據的分類準確率實驗結果。

表6 四種評估方法準確率及Kappa系數對比Tab.6 Comparison of accuracy and kappa coefficient of four evaluation methods

實驗表明，在四種分類算法中光譜色差與主分量融合表征數據的預測準確率均較高，相比于光譜反射率幅值表征方法、主分量表征方法、色差表征方法，光譜色差與主分量特征融合的表征方法對石窟表面不同風化類型及風化程度的表征能力更強。其中在隨機森林算法中的預測準確率最高，高達99.86%，Kappa系數為0.99，與評估結果較好的16通道反射光譜幅值表征方法相比，光譜色差與主分量特征融合表征方法使評估準確率提高了1.37%。利用訓練好的隨機森林算法模型，對石窟表面整體風化類型及風化程度分布情況進行評估，石窟表面風化多光譜成像數據中，待評估區域的總像素點數為1 311 728 points。在不同風化類型及風化程度表征數據下，隨機森林算法對石窟表面每個像素點的預測評估結果如圖7所示。將不同風化類型及風化程度使用不同顏色進行可視化標記，紅色表示強鹽析風化、藍色表示弱鹽析風化、青色表示微鹽析風化、綠色表示積塵風化。

圖7 石窟表面風化整體評估結果Fig.7 Overall evaluation results of weathering on the surface of grotto

其中光譜色差與主分量特征融合表征方法對石窟表面整體風化程度的評估結果為：強鹽析風化區域像素點數為67 625 points（5.16%），弱鹽析風化區域像素點數為368 921 points（28.12%），微鹽析風化區域像素點數為359 457 points（27.40%），積塵風化區域像素點數為515 725 points（39.32%）。

對于圖像中部，弱鹽析風化與強鹽析風化混雜區域，由于表面風化類型同為鹽析風化，其物質狀態特征相似，顏色差異較小，使得主分量特征難以達到精細化表征。

將光譜色差數據與主分量特征數據進行融合表征可以有效減少同種風化不同風化程度間的誤劃分，通過實驗結果對比可以得出，光譜色差與主分量特征融合方法對石窟表面不同風化類型及風化程度的表征能力較強，其評估結果與實際風化狀況最為接近。但對于圖像底部弱鹽析風化與積塵風化混雜區域，中部以弱鹽析風化為主；對于底部突起部分的邊緣，其表面物質主要為積塵，風化類型以積塵風化為主。由于此區域的風化類型及風化程度是由多種風化因素共同作用的結果，此區域上部為大面積的鹽析風化區域，鹽析風化產物隨石窟表面的酥化脫落，覆蓋于下部突起部分的表面，使得此區域的表面的物質狀態特征與鹽析風化區域相近，表面的反射光譜特性與弱鹽析風化區域的反射光譜形態及幅值差異較小。利用隨機森林算法對下部突起部分的風化類型及風化程度進行評估時，光譜色差數據根據顏色的差異將此部分劃分為積塵風化，而主分量特征數據及其與光譜色差融合表征數據將大部分劃分為弱鹽析風化，光譜反射率幅值表征數據將此區域的風化類型進行了折中處理。對于顏色差異較大，物質狀態特征相似的復雜風化區域，光譜色差特征難以達到精細化表征。主分量特征具有石窟表面物質光譜特性及部分巖性信息，可以減小不同風化類型及風化程度區域因顏色差異而帶來的誤劃分，使得突起部分其表面風化類型主要以弱鹽析風化為主，邊緣以積塵風化為主，與實際表面物質狀態及風化狀況基本一致。光譜色差與主分量特征融合數據對于復雜風化區域不同風化類型及風化程度區域，可以達到精細化表征。

圖7（d）為只用主分量特征數據進行訓練，對石窟表面風化類型及風化程度整體表征的實驗結果。對于弱鹽析風化與微鹽析風化、弱鹽析風化與強鹽析風化混雜區域，由于石窟表面物質狀態特征相似，主分量特征對同種風化類型不同風化程度的表征能力較弱，使得評估結果與實際風化狀況差異較大。圖7（e）為光譜色差對不同風化類型及風化程度評估結果，由于色差數據較為單一，不能有效表征石窟表面風化區域的物質狀態特征。在通過隨機森林算法進行分類時，弱鹽析風化與積塵風化混雜區域、微鹽析風化與積塵風化混雜區域，出現大面積的誤劃分，石窟表面整體風化程度預測的準確率降低。16通道的光譜反射率幅值特征數據對不同風化類型及風化程度的表征能力雖然較好，但強鹽析風化與弱鹽析風化混雜區域、弱鹽析風化與積塵風化混雜區域的風化狀況與實際風化狀況及光譜色差與主分量融合表征評估結果對比分析表明，光譜色差與主分量融合的特征數據對風化狀況的表征能力更強，對復雜風化區域的評估結果更加精細、準確。

5 結 論

本文提出了一種基于光譜色差與主分量特征融合的石窟表面風化智能量化評估方法。實驗表明，光譜色差數據可以對同風化類型不同風化程度區域進行有效表征；主分量特征可以有效反映石窟表面風化的物質狀態特征及部分巖性信息。融合光譜色差與主分量特征對石窟表面復雜風化區域的風化類型及風化程度進行表征的方法，有效解決了對于復雜風化區域，單一光譜特征表征能力弱的問題，使得隨機森林算法對石窟表面整體風化程度評估的準確率提高了1.37%。光譜色差與主分量融合數據方法對于石窟表面復雜風化區域的風化類型及風化程度可以達到精細化表征，其評估結果更加準確。