光谱色差与主分量特征融合的石窟表面风化智能评估

曹赤鹏，王慧琴*，王 可，王 展，张 刚，马 涛

（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%。光谱色差与主分量融合数据方法对于石窟表面复杂风化区域的风化类型及风化程度可以达到精细化表征，其评估结果更加准确。