基于空谱特征联合的烟熏壁画线条增强方法

毛锦程，吕书强，侯妙乐，汪万福

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）方法相比，本方法能较好区分壁画线条与其他图案背景，增强线条与其他图案的对比度；与数字影像相比，本次提出的方法能更好地复原壁画色彩，使画面更富有层次感，人物形象更加鲜活，提升壁画的视觉效果。本研究能够为壁画的保护与修复工作提供借鉴与参考。