一种基于高光谱成像土遗址文物表面含水量检测方法

摘要

本发明提供一种基于高光谱成像土遗址文物表面含水量检测方法。所述方法包括以下步骤：在土遗址文物周边采集样本，外加水并混合均匀，将其配制为具有不同外加含水量的梯度标准样本；采集所述梯度标准样本的可见光‑近红外高光谱成像数据；对所述可见光‑近红外高光谱成像数据进行特征波段提取得到特征光谱数据；构建梯度标准样本的特征光谱数据和外加含水量之间的关系，建立土遗址文物表面含水量预测模型，向所述预测模型输入土遗址文物待测样本的可见光‑近红外高光谱成像数据，得到土遗址文物待测样本的预测含水量。

说明书

技术领域

本发明属于文物检测技术领域，具体涉及保护文物含水量的无损检测，尤其涉及一种基于高光谱成像土遗址文物表面含水量检测方法。

背景技术

我国拥有大量粘土类遗址和遗迹，是我国优秀文化遗产中重要的组成部分。土遗址是人类历史上以土为主要材料而建造的具有历史、艺术、科学、社会及文化价值的文化遗产。这些文化遗产是我国悠久历史及灿烂文化的重要载体，是中华文明的实物象征。许多土遗址文物由于受到自然侵蚀和环境变迁等因素，其保护现状不容乐观。在众多的病害中，盐分和水分的破坏是最难治理和最为严重的，也是国内外文物保护研究工作的难点和热点。

在地下水和大气含水的交互作用下，土遗址表面水分成为导致其破坏的潜在威胁，因此，土遗址表面水分含量的获取对于干旱或者潮湿地区的土遗址保护具有重要作用。传统的土壤含水量测试方法会破坏土遗址表面信息，造成环境污染，费时费力，分辨率小，不适用于土遗址文物保护领域。目前高光谱反演土壤水分的模型主要有偏最小二乘回归、神经网络法、逐步回归方法和多元线性回归方法，由于其模型的变量众多，模型运算效率较低，稳定性较差，难以快速准确测定土遗址表面含水量的应用。同时，可见光-近红外光谱的反射率与样品的水分、盐分、有机质、颜色、氧化物含量等很多性质有密切关联，土壤的含水量预测模型很难应用于土遗址表面含水量的测定，因此需要针对特定土遗址表面的材料进行详细研究。

现有资料表明，在文物保护领域尚无通过采用高光谱成像结合化学计量学的方法来建立土遗址表面含水量合理预测模型的研究以及相关报道。

发明内容

针对现有技术中针对土遗址文物表面含水量测定存在的短板，本发明提供一种基于高光谱成像土遗址文物表面含水量检测方法，该方法作为针对土遗址文物表面含水量测定的有效技术，提高了土遗址表面含水量测定的准确性。

为此，本发明提供一种基于高光谱成像土遗址文物表面含水量检测方法。所述方法包括以下步骤：

步骤(i)：在土遗址文物周边采集样本，外加水并混合均匀，将其配制为具有不同外加含水量的梯度标准样本；

步骤(ii)：采集所述梯度标准样本的可见光-近红外高光谱成像数据；

步骤(iii)：对所述可见光-近红外高光谱成像数据进行特征波段提取得到特征光谱数据；

步骤(iv)：构建梯度标准样本的特征光谱数据和外加含水量之间的关系，建立土遗址文物表面含水量预测模型，向所述预测模型输入土遗址文物待测样本的可见光-近红外高光谱成像数据，得到土遗址文物待测样本的预测含水量。

较佳地，可见光高光谱成像的特征光谱数据为400-1000nm波段的光谱数据；近红外高光谱成像的特征光谱数据为1000-2500nm波段的光谱数据。

较佳地，所述方法还包括：对所述可见光-近红外高光谱成像数据进行特征波段提取之前，对可见光-近红外高光谱成像光谱数据进行预处理。

较佳地，采用Savitzky-Golay平滑、标准正态变化、一阶导数和二阶导数的一种或多种的组合对可见光-近红外高光谱成像光谱数据进行预处理；优选地，采用Savitzky-Golay平滑与标准正态变化的组合对标准样本的可见光高光谱成像数据尤其是反射率数据进行预处理；采用Savitzky-Golay平滑与二阶导数的组合对标准样本的近红外高光谱成像数据尤其是反射率数据进行预处理。

较佳地，在每个标准样本中选定区域，对可见光-近红外高光谱的光谱曲线进行算术平均，获得用于预处理的每个标准样本的实际反射率光谱数据。

较佳地，采用连续投影算法和主成分分析法的一种或多种的组合对可见光-近红外高光谱成像数据进行特征波段提取。

较佳地，采用偏最小二乘回归和/或支持向量回归的方法建立土遗址文物表面含水量预测模型；优选地，采用偏最小二乘回归对基于可见光高光谱成像的特征光谱数据建立预测模型，采用支持向量回归对基于近红外高光谱成像的特征光谱数据建立预测模型。

较佳地，所述土遗址文物为粘土基文物。包括但不限于砂岩、石质等一些便于选取周边模拟样本且存在水潜在破坏的硅酸盐质文物。

较佳地，外加含水量＝外加水质量/样本的绝干质量。

较佳地，梯度标准样本的外加含水量分布在0-0.45之间。

本发明所述基于高光谱成像土遗址文物表面含水量检测方法，可针对不同性质、不同地区的土遗址样本建立具有针对性的专用模型，可推广性强。而且，所述方法采用与土遗址文物组成基本一致的样本，进一步提高了预测结果的适配性。经过验证，该方法的模型稳定、预测针对性强且预测精度高。

附图说明

图1是敦煌土样本的可见光高光谱成像(400-1000nm)光谱数据；

图2是敦煌土样本的近红外高光谱成像(1000-2500nm)光谱数据；

图3是采用Savitzky-Golay平滑与标准正态变换组合的方法获取的敦煌土样本的可见光高光谱成像(400-1000nm)光谱预处理数据；

图4是采用Savitzky-Golay平滑与二阶导数组合的方法获取的敦煌土样本的近红外高光谱成像(1000-2500nm)光谱数据；

图5是采用连续投影算法提取的敦煌土样本的可见光高光谱成像(400-1000nm)光谱预处理数据的特征波段；

图6是采用连续投影算法提取的敦煌土样本的近红外高光谱成像(1000-2500nm)光谱预处理数据的特征波段；

图7是连续投影算法的运行情况示意图；

图8是采用偏最小二乘回归对基于可见光高光谱成像(400-1000nm)的特征光谱数据建立预测模型的算法示意图；

图9是采用支持向量回归对基于近红外高光谱成像(1000-2500nm)的特征光谱数据建立预测模型的算法示意图；

图10是基于可见光高光谱成像(400-1000nm)敦煌土样本含水量预测结果；

图11是基于近红外高光谱成像(1000-2500nm)敦煌土样本含水量预测结果。

具体实施方式

通过下述实施方式进一步说明本发明，应理解，下述实施方式仅用于说明本发明，而非限制本发明。以下示例性说明本发明所述基于高光谱成像土遗址文物表面含水量检测方法。

采集样本。在土遗址文物周边采集与土遗址文物组成(基本)一致的样本。土遗址文物周边取样前最好除去植物枝干、砂砾等对样本成分干扰较大的杂质。同时尽量保持样本的平整度与土遗址文物表面基本一致，如此便于后期数据的采集。在样本中混入额外的水，从而配置多个具有不同外加含水量的样本，即为梯度标准样本。这里的梯度指的是不同梯度样本的外加含水量呈现梯度分布。外加含水量＝外加水质量(g)/样本的绝干质量(g)。该样本的绝干质量可为样本在105℃的烘干质量。一些技术方案中，梯度标准样本的外加含水量分布在0-0.45之间。可以根据实际需要选择合适的标准样本数量。例如，该标准样本的数量可为100个。为保证样本与水充分混合，将烘干的土遗址文物样本置于密封的样品袋中，搅拌并喷洒一定质量的去离子水。

采集所述梯度标准样本的可见光-近红外高光谱成像数据。试验中采用上述外加水后的湿样进行可见光-近红外高光谱成像。作为优选，对所述标准样本的可见光-近红外高光谱成像数据(光谱数据)进行算术平均，获得标准样本的实际反射率光谱数据。具体地，对所述梯度标准样本进行可见光-近红外高光谱成像数据的采集，每个标准样本中间位置选定感兴趣区域，并对光谱曲线进行算术平均，获得每个标准样本的实际反射率光谱数据。例如，采用高光谱成像系统数据分析软件HSI Analyzer来导出选择光谱数据的算术平均值。

采用Savitzky-Golay平滑(卷积平衡算法)、标准正态变化、一阶导数和二阶导数等一种或多种的组合对可见光-近红外高光谱成像光谱数据进行预处理。预处理的目的是降低外界因素的干扰以及提高模型的可靠性和精准性。作为优选，采用Savitzky-Golay平滑与标准正态变化的组合对基于可见光高光谱成像(400-1000nm)的标准样本反射率数据进行预处理；采用Savitzky-Golay平滑与二阶导数的组合对基于近红外高光谱成像(1000-2500nm)的标准样本反射率数据进行预处理。

采用Savitzky-Golay平滑与标准正态变化的组合对基于可见光高光谱成像(400-1000nm)的标准样本反射率数据进行预处理的操作为：首先将两组高光谱成像采集的数据导入The Unscramble X 10.4软件中，其中列值为各波段的光谱反射率，行值为样本的个数。在Savitzky-Golay平滑处理方法参数中，由于此数据为线性，故选择多项式次数为1，平滑点数为9，即得到两组Savitzky-Golay平滑后数据。然后将两组预数据分别进行标准正态变换和二阶导数预处理。

其中标准正态变换原理为以下公式：

其中x

同理，用Savitzky-Golay平滑与二阶导数的组合对基于近红外高光谱成像(1000-2500nm)的标准样本反射率数据进行预处理的操作为：将基于近红外高光谱成像的Savitzky-Golay平滑后数据进行二阶导数运算，其运算公式为：

其中，原始光谱为x

上述预处理采用The Unscrambler X 10.4软件进行。具体实施方式中采用不同预处理方法建立了模型，根据决定系数和均方根误差，从中优选了较优的预处理方法。

对预处理的可见光-近红外高光谱成像数据进行特征波段提取。可采用连续投影算法和主成分分析等方法对预处理后的可见光-近红外高光谱成像数据进行特征波段提取。通过建立连续投影算法和主成分分析等特征提取方法，提取了可见光和近红外的特征波段，特征波段涵盖了高光谱数据中的有效信息，减少了光谱的变量数，有利于降低土遗址表面样本含水量预测模型的复杂度，提高了运算效率和预测精度。采用Matlab软件实现连续投影算法提取特征波段的建立，其原理是优选特征波段的数量，得到最小的均方根误差RMSE，此时模型的拟合效果最佳。连续投影算法是一种向前循环的特征波长提取方法，设矩阵X

步骤1：在迭代开始前，任意选取光谱矩阵的一列j，并把建模集的第j列赋值给x

步骤2：把剩下列向量位置的集合记为S，

步骤3：分别计算x

步骤4：记

步骤5：:记

步骤6：令n＝n+1，如果n

基于梯度标准样本建立土遗址文物表面含水量预测模型。构建梯度标准样本的特征光谱数据和外加含水量之间的关系，建立土遗址文物表面含水量预测模型，向所述预测模型输入土遗址文物待测样本的可见光-近红外高光谱成像数据，得到土遗址文物待测样本的预测含水量。作为优选，向所述预测模型输入土遗址文物待测样本的反射率光谱数据。

可以采用偏最小二乘回归和支持向量回归等方法建立土遗址表面含水量预测模型。优选地，采用偏最小二乘回归对基于可见光高光谱成像(400-1000nm)的特征光谱数据建立预测模型，采用支持向量回归对基于近红外高光谱成像(1000-2500nm)的特征光谱数据建立预测模型。上述预测模型的建立在Matlab软件中实现。

其中偏最小二乘回归算法的具体操作步骤如下：

步骤1：对光谱反射率矩阵X和含水量矩阵Y进行分解，其模型为：X＝TP+E，Y＝UQ+F，式中，T和U分别为X和Y矩阵的得分矩阵；P和Q分别为X和Y矩阵的载荷矩阵；E和F分别为X和Y矩阵的偏最小二乘回归拟合残差矩阵；

步骤2：对T和U进行线性回归分析，其表述为：U＝TB，B＝(T

步骤3：根据步骤1的公式得到预测含水量值：Y

其中支持向量回归算法的具体操作步骤如下：

步骤1：支持向量回归是通过一个非线性的函数

步骤2：在模型中引入ξ为敏感因子，假设容忍f(x)与y之间最多有ξ的偏差。由此模型转变为凸优化问题，即

步骤3：本发明的核函数选用高斯径向基(RBF)函数，具体形式为：k(x，x

此时支持向量回归的性能由误差惩罚因子C和高斯径向基函数宽度γ。误差惩罚因子C过大或者过小都会导致模型泛化能力下降，高斯径向基函数宽度γ反映支持向量之间的向关联程度，过大或者过小也会影响模型的精度。本发明采用粒子群优化算法得到最优参数：C＝11，γ＝0.0039。

一些技术方案中，将用于可见光-近红外高光谱成像采集的标准样本划分为建模样本和验证样本。所述建模样本和验证样本的数量比例可以根据需要进行调整。例如，建模样本和验证样本的数量比例为7:3。建模样本用于建立含水量预测模型。验证样本用于验证样本用于检验预测模型。作为优选，所述建模样本和验证样本的外加含水量均在标准样本的外加含水量分布区间内等数值分布。

通过以下公式根据决定系数R

上式中，y

本发明提供一种基于可见光-近红外高光谱成像土遗址表面含水量的检测方法，高效便捷，对于基于可见光-近红外高光谱成像的其他土遗址表面含水量具有类似的指导意义。

下面进一步例举实施例以详细说明本发明。同样应理解，以下实施例只用于对本发明进行进一步说明，不能理解为对本发明保护范围的限制，本领域的技术人员根据本发明的上述内容作出的一些非本质的改进和调整均属于本发明的保护范围。下述示例具体的工艺参数等也仅是合适范围中的一个示例，即本领域技术人员可以通过本文的说明做合适的范围内选择，而并非要限定于下文示例的具体数值。

以典型敦煌土作为土遗址文物样本，结合图1至图11来详细说明本发明基于可见光-近红外高光谱成像对土遗址样本含水量检测的方法。

步骤1：将敦煌土样本按梯度配制成不同外加含水量的(0-0.45)标准样本。为保证敦煌土样本与水充分混合，将105℃条件下烘干的敦煌土置于密封的样品袋中，搅拌并喷洒一定质量的去离子水。外加含水量＝外加水质量/敦煌土样本在105℃烘干质量。标准样本共有100个。将配制完成的标准样本称重记录，并放置24h备用。

步骤2：高光谱成像分为两种波段，即可见光高光谱成像(400-1000nm)和近红外高光谱成像(1000-2500nm)。分别采用可见光高光谱成像(400-1000nm)和近红外高光谱成像(1000-2500nm)对配制好的敦煌土标准样本进行数据采集。采集完成后将敦煌土样本置于105℃烘箱中烘干12h，计算样本的样本含水量真实值，并做记录。含水量真实值＝(样本质量-样本绝干质量)/样本绝干质量。

对基于可见光-近红外高光谱成像的每个敦煌土样本的数据进行分析，在每个样本中间位置选定感兴趣区域，并对光谱曲线进行算术平均，获得每个敦煌土样本的实际反射率光谱数据。基于可见光高光谱成像的敦煌土样本反射率数据如图1所示。基于近红外高光谱成像的敦煌土样本反射率数据如图2所示。

步骤3：对步骤2获得的基于可见光高光谱成像(400-1000nm)的敦煌土样本反射率数据，采用Savitzky-Golay平滑与标准正态变化组合的方法进行数据预处理，其结果如图3所示。同时，对基于近红外高光谱成像(1000-2500nm)敦煌土样本反射率数据，采用Savitzky-Golay平滑与二阶导数组合的方法进行数据预处理，其结果如图4所示。

步骤4：采用连续投影算法对步骤3获得的基于可见光高光谱成像(400-1000nm)的光谱预处理数据进行特征长提取，设定最大特征波长为18时，此时得到最小的均方根误差(RMSE＝0.0163)，其特征波长分别为414.384nm、523.410nm、629.895nm、788.908nm和932.025nm，其具体结果如图5所示。采用连续投影算法对步骤3中获得的基于近红外高光谱成像(1000-2500nm)的光谱预处理数据进行特征长提取，设定最大特征波长为8时，此时得到最小的均方根误差(RMSE＝0.0133)，其特征波长分别为1149.251nm、1233.7nm、1464.197nm、2214.801nm和2262.161nm，其具体结果如图6所示。连续投影算法在Matlab软件中进行，运行算法如图7所示。

步骤5：将用于可见光-近红外高光谱成像采集的100个敦煌土标准样本划分为建模样本和验证样本。选取70个敦煌土样本作为建模样本，剩余的30个样本作为验证样本。采用偏最小二乘回归对基于可见光高光谱成像(400-1000nm)的特征光谱数据建立预测模型，其模型算法如图8所示。采用支持向量回归对基于近红外高光谱成像(1000-2500nm)的特征光谱数据建立预测模型，其模型算法如图9所示。两种模型均采用Matlab软件运行。

步骤6：对预测模型的结果进行分析评价。基于可见光高光谱成像(400-1000nm)的特征光谱数据与已知外加含水量建立的预测模型，其建模决定系数R

模型精度由决定系数R

还采用了其他数据预处理方法结合建立不同含水量预测模型。RAW为未经任何处理的原始数据；SG为Savitzky-Golay平滑；SG+FD为Savitzky-Golay平滑和一阶导数组合；SG+SD为Savitzky-Golay平滑和二阶导数组合；SG+MSC为Savitzky-Golay平滑和多元散射校正组合；SG+SNV为Savitzky-Golay平滑和标准正态变换组合。

表1基于可见光高光谱成像采用不同预处理方法的含水量预测模型的结果

表2基于近红外高光谱成像采用不同预处理方法的含水量预测模型的结果

其中，R