一种联合土壤室内外光谱的重金属含量反演方法与流程

本发明涉及环境监测领域，特别涉及一种联合土壤室内外光谱进行重金属含量反演的方法。

背景技术：

利用常规的实地采样方式获取大范围高密度土壤重金属污染信息时，通常面临成本耗费巨大、效率较低的难题，且难以保证常态化快速高效动态监测得以实施，无法满足国家在广域地理空间快速高效查清重金属污染空间格局、发现重金属污染时空演化规律与形成机制、控制污染范围扩大和转移、合理规划农业生产、减少污染产品损害国民健康的重大需求。

目前，基于光谱分析技术的土壤重金属含量反演方法得到快速的发展，已经成为土壤重金属污染调查的主要手段。其主要原理是：利用地物光谱仪观测土壤样品室内光谱数据，探索土壤重金属含量与土壤全波段光谱(350–2500nm)间的响应关系，构建土壤重金属含量高光谱反演模型。该方法相对传统地统计分析方法具有花费较少、效率较高的优势，但由于室内、室外大气等环境条件不同导致的土壤光谱差异，使得室内光谱难以有效直接应用于野外大范围土壤重金属污染调查的条件，亟需发明一种可在野外直接利用土壤室外光谱大范围高效反演土壤重金属含量的方法。

技术实现要素：

为了有效解决传统基于土壤室内光谱反演重金属含量的局限性，以实现野外土壤重金属污染调查过程中直接利用土壤室外光谱大范围高效反演土壤重金属含量的目的，本发明提供了一种联合土壤室内外光谱的重金属含量反演方法。

为了实现上述技术目的，本发明的技术方案是，

一种联合土壤室内外光谱的重金属含量反演方法，包括以下步骤：

步骤1：在需要进行土壤重金属含量检测的污染土壤调查区采集土壤样本，并在室外采集土壤样本的室外光谱数据，然后将土壤样本送至实验室内进行重金属含量检测，并在室内采集土壤样本的室内光谱数据，最后基于采样点土壤重金属含量实测数据与土壤室内光谱数据，构建基于室内光谱的土壤重金属含量反演的全要素主成分逐步回归模型，并提取特征波段主成分；

步骤2：选取合适的转换样本集构建土壤室内、室外光谱的关联转换模型，对全部采样土壤室外光谱数据做光谱转换及系列光谱预处理；

步骤3：对步骤1中提取的特征波段主成分与步骤2中转换处理后的室外土壤光谱主成分做相关性分析，提取显著相关主成分作为解释变量，以步骤1中测定的土壤重金属含量作为因变量，构建基于土壤室外光谱的重金属含量反演的回归模型。

所述的一种联合土壤室内外光谱的重金属含量反演方法，步骤1中所述的土壤重金属含量及其室内光谱数据通过以下步骤获得：

在需要进行土壤重金属含量检测的污染土壤调查区采用格网对室外污染土壤研究靶区进行格网界定，并在格网之间形成的交点处设置土壤采样点采集不少于40个土壤样本；在实验室内对采集到的土壤样本采用三酸消化——原子吸收分光光度法化学分析测定土壤样品重金属含量，同时，使用地物光谱仪测定土壤室内光谱数据；

步骤1中所述的采样点土壤室外光谱数据通过以下步骤获得：

步骤1进行野外土壤样本采集的同时，运用地物光谱仪监测采集土壤样本的室外光谱数据。

所述的一种联合土壤室内外光谱的重金属含量反演方法，步骤1中所述的构建基于土壤室内光谱的重金属含量反演的全要素主成分逐步回归模型并提取特征波段主成分，包括以下步骤：

步骤1)：采用光谱重采样、标准正态变换、一阶/二阶微分的预处理方法对土壤室内光谱数据进行预处理，构建基于土壤室内光谱的重金属含量的全要素主成分逐步回归模型；

步骤2)：将步骤1)模型的特征波段主成分值提取，作为重金属含量的土壤室外光谱回归建模自变量筛选的判定因子。

所述的一种联合土壤室内外光谱的重金属含量反演方法，步骤1)中，对土壤室内光谱数据进行光谱重采样的预处理包括以下步骤：

以预设长度为间隔单元，采用以下公式对测定的土壤室内光谱数据进行光谱重采样：

r＝[r(λ1)+···+r(λi)+···+r(λn)]/n

其中，r表示高光谱曲线经过重采样后的反射率值；n为间隔单元范围内所包含的光谱波段数，r(λi)表示第i个光谱波段λi的反射率；

步骤1)中，对土壤室内光谱数据进行标准正态变换的预处理包括以下步骤：

基于重采样后的室内光谱，基于加权平均化方法采用以下公式进行标准正态变量变换以完成归一化处理：

其中，n为间隔单元范围内所包含的光谱波段数；表示第i个样本n个波段光谱反射率的平均值；xi,j表示第i个样本的第j个波段的光谱反射率值；δ为自定义的偏移量；wi表示第i个样本的标准差；x′i,j表示第i个样本的第j个波段标准正态变换完成归一化处理的结果值；

步骤1)中，对土壤室内光谱数据进行一阶/二阶微分的预处理包括以下步骤：

基于一阶/二阶微分采用以下公式计算各样本光谱数据以去除干扰：

r′(λi)＝[r(λi)-r(λi-1)]/2Δλ

r″(λi)＝[r(λi)-2r(λi-1)+r(λi-2)]/(2Δλ)²

其中，r(λi)为第i个光谱波段λi，r′(λi)和r″(λi)分别为光谱波段λi的一阶/二阶微分光谱值，Δλ为光谱波段λi-1与λi之间的间隔单元。

所述的一种联合土壤室内外光谱的重金属含量反演方法，步骤1)中，构建基于室内光谱的土壤重金属含量的全要素主成分逐步回归模型包括以下步骤：

以经过预处理后的土壤室内光谱数据为输入自变量，步骤1中测定的土壤样品重金属含量为因变量，构建全要素主成分逐步回归模型：

Cj＝v1jX1+v2jX2+…vijXi…+vmjXm

Y＝a0+a1C1+a2C2+…anCn

其中，Cj代表主成分变换后的第j个新成分，Xi代表第i个波段的光谱反射率；m代表主成分变换前的光谱变量数，vij代表m个原始变量的相关系数矩阵中第i个特征值对应的特征向量，Y代表土壤重金属含量估算值，an代表第n个主成分变量对应的重金属含量反演模型系数。

所述的一种联合土壤室内外光谱的重金属含量反演方法，所述的步骤2包括以下步骤：

步骤①：删除室内、室外土壤光谱水汽影响波段，获得有效土壤光谱数据；

步骤②：利用Kennard-Stone算法通过计算采样土壤室内光谱之间的欧氏距离，选择转换集样品：

其中，rik与rjk分别代表样本i和样本j的第k个波段的光谱反射率值；dij代表采样土壤样本i和样本j之间的欧氏距离，p为样品光谱的波长个数；

步骤③：采用直接校正算法构建步骤2)选取的转换集样品的室内、室外光谱转换模型：

其中，E为残差矩阵，XQ为分析样本室内光谱值，Xq为分析样本室外光谱值，，为分析样本室内中心化处理后光谱值，为分析样本室外中心化处理后光谱值；

步骤④：采用步骤③构建的土壤室内外光谱转换模型将土壤样品室外光谱进行光谱转换，得到转换后的土壤室外光谱数据；

步骤⑤：采用光谱重采样—标准正态变换—一阶/二阶微分系列光谱预处理方法对步骤④转换后的土壤室外光谱进行系列光谱数据预处理。

所述的一种联合土壤室内外光谱的重金属含量反演方法，步骤3包括以下步骤：

步骤一：采用主成分分析方法提取转换处理后的土壤室外光谱主成分；

步骤二：采用Pearson相关性分析法计算步骤一提取的光谱主成分数据与重金属含量土壤室内光谱反演模型选入的主成分判定因子的相关系数；

步骤三：分析步骤二相关系数计算结果，选取与判定因子相关系数最大的光谱主成分数据作为重金属含量反演模型输入自变量，以步骤1测得的土壤样品重金属含量为因变量，构建线性回归模型。

本发明的技术效果在于，基于土壤室内、室外光谱内在的关联一致性，顾及地理环境要素对于土壤光谱反射率值的变化影响，本发明基于土壤室内、室外光谱的关联转换模型构建实现了联合土壤室内外光谱的重金属含量反演。具体地，本发明提出，基于土壤室内光谱的系列光谱预处理，构建土壤重金属含量反演的全要素主成分逐步回归模型，并提取特征波段主成分；采用Kennard-Stone算法选择合适的转换集构建土壤室内、室外光谱的关联转换模型，并同样对转换后的土壤室外光谱做系列光谱预处理；最后，融合重金属含量的土壤室内光谱反演模型提取的特征波段主成分，基于转换处理后的土壤室外光谱数据，构建重金属含量高光谱反演模型。这是当前土壤重金属含量高光谱遥感反演多局限于采样土壤室内光谱特征研究，难以有效直接应用于野外大范围土壤重金属污染调查的条件下，实现野外直接利用土壤室外光谱的大范围高效土壤重金属含量反演。

附图说明

图1为本发明室内外光谱联合转换关键技术示意图；

图2中，(a)为基于土壤室内光谱的重金属含量全要素主成分回归模型结果示例图，(b)为基于模型结果的特征波段主成分提取结果示例图。

图3中，(a)为联合室内外光谱的土壤重金属含量反演模型结果示例图，(b)为基于模型结果的重金属含量的土壤光谱响应特征波段识别结果示例图。

具体实施方式

下面是对本发明一个优选实施例，结合附图进行的详细说明。

1、构建基于土壤室内光谱的重金属含量反演的全要素主成分逐步回归模型并提取特征波段主成分，本发明采用的土壤重金属含量室内光谱反演特征波段主成分提取的步骤包括：

步骤1)：在土壤调查区采用格网采集方式进行采样，以1m×1m分辨率对污染土壤研究靶区进行网格界定即画网格，并在格网之间形成的交点处设置采样点采集土壤样本；采用ASD地物光谱仪采集土壤采样点室外光谱数据(350-2500nm)；然后将采集的土壤样本带入实验室内，采用九点采样法采集上述采样点土壤样品，采用三酸消化——原子吸收(石墨炉)分光光度法实验室化学分析测定土壤样品重金属含量，同样运用ASD地物光谱仪测定土壤样本室内光谱数据(350–2500nm)；这样就得到了同一土壤样本在室外监测和室内监测不同环境条件下获得的光谱数据，也就是后面转换用到的两组光谱。

步骤2)：以10nm为间隔单元对步骤1)测定的土壤室内光谱数据进行光谱重采样；具体的计算公式如下：

r＝[r(λ1)+…+r(λi)+…+r(λn)]/n

式中，r表示高光谱曲线经过重采样后的反射率值；n为10nm间隔单元范围内所包含的光谱波段数，r(λi)表示第i个光谱波段λi的反射率；

步骤3)：基于重采样后的光谱曲线，具体采用加权平均化方法对其进行标准正态变量变换，以完成采集样本光谱曲线的归一化处理。其具体的计算公式如下：

其中，n为光谱波段的数量；表示第i个样本n个波段光谱反射率的平均值；xi,j表示第i个样本的第j个波段的光谱反射率值；δ为自定义的偏移量；wi表示第i个样本的标准差；x′i,j表示第i个样本的第j个波段标准正态变换完成归一化处理的结果值；

步骤4)：采用一阶/二阶微分处理计算各样本光谱数据，以去除部分线性或类线性光谱背景值对重金属光谱的干扰(如：噪声光谱)，达到增强光谱特征差异、提取光谱特征吸收带的目的。具体计算公式如下：

r′(λi)＝[r(λi)-r(λi-1)]/2Δλ

r″(λi)＝[r(λi)-2r(λi-1)+r(λi-2)]/(2Δλ)²

其中，r(λi)为第i个光谱波段λi，r′(λi)和r″(λi)分别为光谱波段λi的一阶/二阶微分光谱值，Δλ为光谱波段λi-1与λi之间的间隔单元10nm；

步骤5)：以上述预处理后的土壤室内光谱数据为输入自变量，步骤1)中测定的土壤样品重金属含量为因变量，构建全要素主成分逐步回归模型。具体计算公式如下：

Cj＝v1jX1+v2jX2+…vijXi…+vmjXm

Y＝a0+a1C1+a2C2+…anCn

式中，Cj代表主成分变换后的第j个新成分，Xi代表第i个波段的光谱反射率；m代表主成分变换前的光谱变量数，vij代表m个原始变量的相关系数矩阵中第i个特征值对应的特征向量，Y代表土壤重金属含量估算值，an代表第n个主成分变量对应的重金属含量反演模型系数；

步骤6)：将步骤5)中运用全要素主成分逐步回归方法建立的模型自变量，即基于土壤室内光谱的重金属含量反演模型的特征波段主成分值提取，作为重金属含量的土壤室外光谱回归建模自变量筛选的判定因子；

2、土壤室内外光谱的关联转换模型构建及土壤室外光谱的转换处理，本发明采用的土壤室内、室外光谱转换模型构建及对土壤室外光谱的转换与系列光谱预处理步骤包括：

步骤1)：去除室内、室外土壤光谱水汽影响波段(1360-1490nm和1810-1960nm)，获得有效土壤光谱数据；

步骤2)：利用Kennard-Stone算法通过计算采样土壤室内光谱之间的欧氏距离，选择转换集样品。欧式距离具体计算公式如下：

其中，rik与rjk分别代表样本i和样本j的第k个波段的光谱反射率值；dij代表采样土壤样本i和样本j之间的欧氏距离，p为样品光谱的波长个数；

Kennard-Stone算法具体实现步骤为：

(1)首先，选择两两样本之间距离最大的两个样本作为第一个和第二个转换集样品；

(2)然后，分别计算剩余的样本与已选样本之间的距离；

(3)对于每个剩余样本而言，其与已选样品之间的最短距离被选择，然后选择这些最短距离中相对最长的距离所对应的样本，作为第三个转换集样品；

(4)重复步骤(3)直至所选的转换集样品的个数等于事先确定的数目为止；

步骤3)：采用直接校正算法构建步骤2)选取的转换集样品的室内、室外光谱转换模型，具体原理如下：

(1)同一转换土壤样品室内、室外光谱数据的差异代表了在两种测量环境下土壤光谱特征的差异，两者之间存在如下的转换关系式：

XQ＝XqA+E

其中，XQ为分析样本室内光谱值，Xq为分析样本室外光谱值，A为室内、外光谱转换矩阵，E为残差矩阵；

(2)当将分析样本室内光谱值XI和室外光谱值Xi分别进行中心化处理以后，依据上述的转换方程有如下关系式：

其中，为分析样本室内中心化处理后光谱值，为分析样本室外中心化处理后光谱值，则光谱转换矩阵A的计算值为：

(3)将上式带入可得残差矩阵：

由此，土壤室内、室外光谱转换模型构建完成；

步骤4)：采用步骤3)构建的土壤室内外光谱转换模型将土壤样品室外光谱进行光谱转换，得到转换后的土壤室外光谱数据；

步骤5)：采用光谱重采样—标准正态变换—一阶/二阶微分系列光谱预处理方法对步骤4)转换后的土壤室外光谱进行系列光谱数据预处理；

3、联合土壤室内外光谱的重金属含量反演回归模型构建，本发明采用的融合重金属含量的土壤室内光谱反演模型提取的特征波段主成分，构建基于转换处理后的土壤室外光谱的重金属含量高光谱反演的回归模型的步骤包括：

步骤1)：采用主成分分析方法提取转换处理后的土壤室外光谱主成分。具体计算公式如下：

Dj＝q1jX1+q2jX2+…qijXi…+qmjXm

式中，第i个波段反射率用Xi表示；变换后成分用Dj表示，其系数qmj是m个变量相关系数矩第j个特征值对应的特征向量；

步骤2)：采用Pearson相关性分析法计算步骤1)提取的光谱主成分数据与重金属含量室内土壤光谱反演模型选入的主成分判定因子的相关系数。其计算公式为：

式中，Dj代表n个样本的第j个主成分构成的向量，Cp代表基于室内土壤光谱构建的重金属含量反演模型提取的第p个主成分变量构成的矩阵，σ代表标准差，相关系数t检验过程中0.05水平显著的波段为显著相关的主成分；

步骤3)：分析步骤2)相关系数计算结果，选取与判定因子相关系数最大的光谱主成分数据作为重金属含量反演模型输入自变量以土壤样品重金属含量为因变量，构建线性回归模型。具体计算公式如下：

Y＝b0+b1D1+b2D2+…bnDn

式中，Y代表土壤重金属含量估算值，Dj代表主成分变换后的第j个新成分，bn代表第n个主成分变量对应的重金属含量反演模型系数；

至此，联合土壤室内外光谱的重金属含量反演模型构建完成。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施方式仅限于此，对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单的推演或替换，都应当视为属于本发明由所提交的权利要求书确定专利保护范围。