融合多源环境变量与光谱信息的土壤重金属含量反演方法

摘要

本发明公开了融合多源环境变量与光谱信息的土壤重金属含量反演方法，方法包括：首先获取土壤样本数据、混合光谱数据和环境变量数据；然后根据所述土壤样本数据、混合光谱数据和环境变量数据，构建重金属反演模型；最后根据所述重金属反演模型对目标区域进行土壤重金属含量反演，确定所述目标区域的重金属分布图。本发明提高了预测精度且普适性强，可广泛应用于环境监测技术领域。

权利要求

1.融合多源环境变量与光谱信息的土壤重金属含量反演方法，其特征在于，包括：获取土壤样本数据、混合光谱数据和环境变量数据；

根据所述土壤样本数据、混合光谱数据和环境变量数据，构建重金属反演模型；

根据所述重金属反演模型对目标区域进行土壤重金属含量反演，确定所述目标区域的重金属分布图；

所述方法还包括从所述环境变量数据中确定重要环境变量的过程，具体包括：将所述环境变量数据作为因子，将土壤重金属浓度作为变量，通过地理探测器模型中的分异及因子探测器进行重要性分析，获取每个环境变量对于重金属浓度分布的影响力，并依据影响程度从高到低进行排序，得到第一序列；

将所述环境变量数据作为因子，将土壤重金属浓度作为变量，通过随机森林重要性分析的方法进行重要性分析，获取每个环境变量对于重金属浓度分布的影响力，并依据影响程度从高到低进行排序，得到第二序列；

从所述第一序列和所述第二序列中选取目标个数的环境变量作为重要环境变量。

2.根据权利要求1所述的融合多源环境变量与光谱信息的土壤重金属含量反演方法，其特征在于，所述获取土壤样本数据，包括：采用五点采样法采集样本农田的表层土壤样本，并利用GPS记录所述样本农田的中心点的经纬度坐标；其中，所述表层土壤样本均匀分布于所述样本农田；

将采集到的表层土壤样本置于干燥通风处进行自然风干，并对表层土壤样本进行杂质剔除，然后将表层土壤样本进行研磨，并利用100目筛对研磨后的样本进行过筛；

采用HF‑HNO3‑HClO4对样本进行消煮，然后利用等离子体质谱仪测定重金属镉含量；利用碳氮元素分析仪测定土壤有机碳含量；用pH计测定土壤酸碱度。

3.根据权利要求1所述的融合多源环境变量与光谱信息的土壤重金属含量反演方法，其特征在于，所述获取混合光谱数据，包括：利用调查区域数字高程模型以及地面控制点经纬度坐标对遥感影像进行正射校正，获取原始光谱数据；

对原始光谱数据进行辐射定标以及大气校正，将光谱数据单位转换至地表反射率；

提取采样点对应位置遥感影像像元的原始光谱数据，利用一阶微分、倒数对数的光谱变换方法对原始光谱数据进行光谱变换，分别得到一阶微分变换结果和倒数对数变换结果；

分别计算一阶微分变换结果、倒数对数变换结果、原始光谱数据各个波段与重金属浓度之间的相关系数；

对比同一波段不同处理方式下光谱数据与重金属浓度之间相关性的大小，选用每一波段中相关系数最大的光谱数据组成混合光谱数据。

4.根据权利要求1所述的融合多源环境变量与光谱信息的土壤重金属含量反演方法，其特征在于，所述获取环境变量数据，包括：根据调查区域的数字高程模型，通过ArcGIS软件生成表征调查区域地形的栅格数据，通过SAGA GIS软件生成地形湿度指数数据；其中，所述栅格数据包括坡度、坡向、平面曲率、剖面曲率；

若调查区域的路网和水系的矢量数据的质量不满足预设要求，则对所述矢量数据进行勾画修正，然后根据调查区域的范围生成栅格数据，并利用ArcGIS软件计算栅格数据中每个像元离最近的河流或者道路的距离，并将该距离的值赋值给对应的像元，生成用来表征交通、水流影响程度的栅格数据；

利用监督分类方法对遥感影像进行分类，将遥感影像分为不同地物类型，其中，地物类型包括农田、人造地表、水体及森林；提取调查区域人造地表，生成表征人类活动的栅格数据；

若调查区域周边存在重金属污染源，则根据调查区域的范围生成栅格数据，并基于污染源位置，将栅格数据中每个像元离污染源的距离赋值给该像元，生成表征离污染源距离的栅格数据。

5.根据权利要求1所述的融合多源环境变量与光谱信息的土壤重金属含量反演方法，其特征在于，所述根据所述土壤样本数据、混合光谱数据和环境变量数据，构建重金属反演模型，包括：根据采样点的经纬度坐标，提取采样点对应位置的影响因子数据，构建数据集；

基于十折交叉验证方法，将所述数据集随机分为10等份，选用其中9份样本作为训练集，利用偏最小二乘法、Bp神经网络、支持向量机、随机森林四种建模方法，以混合光谱数据及重要影响因子数据为自变量，以土壤重金属浓度数据为因变量训练模型；

选用另外1份样本为验证集，对四种模型进行精度评价，计算模型的决定系数、平均绝对误差、均方根误差；

重复对模型进行训练，获取满足训练条件的模型构建得到重金属反演模型。

6.根据权利要求1所述的融合多源环境变量与光谱信息的土壤重金属含量反演方法，其特征在于，所述根据所述重金属反演模型对目标区域进行土壤重金属含量反演，确定所述目标区域的重金属分布图，包括：根据混合光谱筛选结果，生成目标区域范围内混合光谱栅格数据；

根据目标区域范围生成栅格数据，以混合光谱栅格数据以及栅格形式的重要影响因子为自变量，选用所述重金属反演模型反演栅格中所有像元的重金属浓度，得到所述目标区域的重金属浓度空间分布图；

提取所述目标区域的农田区域，根据所述农田区域对所述重金属浓度空间分布图进行裁剪，得到所述目标区域的重金属分布图。

7.一种电子设备，其特征在于，包括处理器以及存储器；

所述存储器用于存储程序；

所述处理器执行所述程序实现如权利要求1至6中任一项所述的方法。

8.一种计算机可读存储介质，其特征在于，所述存储介质存储有程序，所述程序被处理器执行实现如权利要求1至6中任一项所述的方法。

说明书

融合多源环境变量与光谱信息的土壤重金属含量反演方法

技术领域

[0001] 本发明涉及环境监测技术领域，尤其是融合多源环境变量与光谱信息的土壤重金属含量反演方法。

背景技术

[0002] 随着农用化学物质使用量不断增加，土壤重金属污染问题逐步凸显。

[0003] 目前针对大范围土壤重金属含量估算的研究，大多单一使用光谱信息，即将与土壤重金属含量相关的特征波段及其组合指数与实测重金属含量进行统计建模，限制了模型的预测精度和普适性。

[0004] 少量的知识应用会针对土壤重金属污染机理，引入环境因子辅助反演模型的建立。这种方式相比传统建模方法能够对光谱信息进行补充，一定程度上提高了模型精度。然而，这些研究尚未对辅助环境因子数据进行深入分析，引入了较多的不相关环境因子，造成了数据冗余，导致模型精度仍受到限制。

发明内容

[0005] 有鉴于此，本发明实施例提供一种预测精度高且普适性强的，融合多源环境变量与光谱信息的土壤重金属含量反演方法。

[0006] 本发明的一方面提供了一种融合多源环境变量与光谱信息的土壤重金属含量反演方法，包括：

[0007] 获取土壤样本数据、混合光谱数据和环境变量数据；

[0008] 根据所述土壤样本数据、混合光谱数据和环境变量数据，构建重金属反演模型；

[0009] 根据所述重金属反演模型对目标区域进行土壤重金属含量反演，确定所述目标区域的重金属分布图。

[0010] 可选地，所述获取土壤样本数据，包括：

[0011] 采用五点采样法采集样本农田的表层土壤样本，并利用GPS记录所述样本农田的中心点的经纬度坐标；其中，所述表层土壤样本均匀分布于所述样本农田；

[0012] 将采集到的表层土壤样本置于干燥通风处进行自然风干，并对表层土壤样本进行杂质剔除，然后将表层土壤样本进行研磨，并利用100目筛对研磨后的样本进行过筛；

[0013] 采用HF‑HNO3‑HClO4对样本进行消煮，然后利用等离子体质谱仪测定重金属镉含量；利用碳氮元素分析仪测定土壤有机碳含量；用pH计测定土壤酸碱度。

[0014] 可选地，所述获取混合光谱数据，包括：

[0015] 利用调查区域数字高程模型以及地面控制点经纬度坐标对遥感影像进行正射校正，获取原始光谱数据；

[0016] 对原始光谱数据进行辐射定标以及大气校正，将光谱数据单位转换至地表反射率；

[0017] 提取采样点对应位置遥感影像像元的原始光谱数据，利用利用一阶微分、倒数对数的光谱变换方法对原始光谱数据进行光谱变换，分别得到一阶微分变换结果和倒数对数变换结果；

[0018] 分别计算一阶微分变换结果、倒数对数变换结果、原始光谱数据各个波段与重金属浓度之间的相关系数；

[0019] 对比同一波段不同处理方式下光谱数据与重金属浓度之间相关性的大小，选用每一波段中相关系数最大的光谱数据组成混合光谱数据。

[0020] 可选地，所述获取环境变量数据，包括：

[0021] 根据调查区域的数字高程模型，通过ArcGIS软件生成表征调查区域地形的栅格数据，通过SAGA GIS软件生成地形湿度指数数据；其中，所述栅格数据包括坡度、坡向、平面曲率、剖面曲率；

[0022] 若调查区域的路网和水系的矢量数据的质量不满足预设要求，则对所述矢量数据进行勾画修正，然后根据调查区域的范围生成栅格数据，并利用ArcGIS软件计算栅格数据中每个像元离最近的河流或者道路的距离，并将该距离的值赋值给对应的像元，生成用来表征交通、水流影响程度的栅格数据；

[0023] 利用监督分类方法对遥感影像进行分类，将遥感影像分为不同地物类型，其中，地物类型包括但不限于农田、人造地表、水体及森林；提取调查区域人造地表，生成表征人类活动的栅格数据；

[0024] 若调查区域周边存在重金属污染源，则根据调查区域的范围生成栅格数据，并基于污染源位置，将栅格数据中每个像元离污染源的距离赋值给该像元，生成表征离污染源距离的栅格数据。

[0025] 可选地，所述方法还包括从所述环境变量数据中确定重要环境变量的过程，具体包括：

[0026] 将所述环境变量数据作为因子，将土壤重金属浓度作为变量，通过地理探测器模型中的分异及因子探测器进行重要性分析，获取每个环境变量对于重金属浓度分布的影响力，并依据影响程度从高到低进行排序，得到第一序列；

[0027] 将所述环境变量数据作为因子，将土壤重金属浓度作为变量，通过随机森林重要性分析的方法进行重要性分析，获取每个环境变量对于重金属浓度分布的影响力，并依据影响程度从高到低进行排序，得到第二序列；

[0028] 从所述第一序列和所述第二序列中选取目标个数的环境变量作为重要环境变量。

[0029] 可选地，所述根据所述土壤样本数据、混合光谱数据和环境变量数据，构建重金属反演模型，包括：

[0030] 根据采样点的经纬度坐标，提取采样点对应位置的影响因子数据，构建数据集；

[0031] 基于十折交叉验证方法，将所述数据集随机分为10等份，选用其中9份样本作为训练集，利用偏最小二乘法、Bp神经网络、支持向量机、随机森林四种建模方法，以混合光谱数据及重要影响因子数据为自变量，以土壤重金属浓度数据为因变量训练模型；

[0032] 选用另外1份样本为验证集，对四种模型进行精度评价，计算模型的决定系数、平均绝对误差、均方根误差；

[0033] 重复对模型进行训练，获取满足训练条件的模型构建得到重金属反演模型。

[0034] 可选地，所述根据所述重金属反演模型对目标区域进行土壤重金属含量反演，确定所述目标区域的重金属分布图，包括：

[0035] 根据混合光谱筛选结果，生成目标区域范围内混合光谱栅格数据；

[0036] 根据目标区域范围生成栅格数据，以混合光谱栅格数据以及栅格形式的重要影响因子为自变量，选用所述重金属反演模型反演栅格中所有像元的重金属浓度，得到所述目标区域的重金属浓度空间分布图；

[0037] 提取所述目标区域的农田区域，根据所述农田区域对所述重金属浓度空间分布图进行裁剪，得到所述目标区域的重金属分布图。

[0038] 本发明实施例的另一方面还提供了一种电子设备，包括处理器以及存储器；

[0039] 所述存储器用于存储程序；

[0040] 所述处理器执行所述程序实现如前面所述的方法。

[0041] 本发明实施例的另一方面还提供了一种计算机可读存储介质，所述存储介质存储有程序，所述程序被处理器执行实现如前面所述的方法。

[0042] 本发明实施例的另一方面还提供了一种计算机程序产品，包括计算机程序，所述计算机程序被处理器执行时实现如前面所述的方法。

[0043] 本发明的实施例首先获取土壤样本数据、混合光谱数据和环境变量数据；然后根据所述土壤样本数据、混合光谱数据和环境变量数据，构建重金属反演模型；最后根据所述重金属反演模型对目标区域进行土壤重金属含量反演，确定所述目标区域的重金属分布图。本发明提高了预测精度且普适性强。

附图说明

[0044] 为了更清楚地说明本申请实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

[0045] 图1为本发明实施例提供的整体步骤流程图。

具体实施方式

[0046] 为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。

[0047] 针对现有技术存在的问题，本发明的一方面提供了一种融合多源环境变量与光谱信息的土壤重金属含量反演方法，包括：

[0048] 获取土壤样本数据、混合光谱数据和环境变量数据；

[0049] 根据所述土壤样本数据、混合光谱数据和环境变量数据，构建重金属反演模型；

[0050] 根据所述重金属反演模型对目标区域进行土壤重金属含量反演，确定所述目标区域的重金属分布图。

[0051] 可选地，所述获取土壤样本数据，包括：

[0052] 采用五点采样法采集样本农田的表层土壤样本，并利用GPS记录所述样本农田的中心点的经纬度坐标；其中，所述表层土壤样本均匀分布于所述样本农田；

[0053] 将采集到的表层土壤样本置于干燥通风处进行自然风干，并对表层土壤样本进行杂质剔除，然后将表层土壤样本进行研磨，并利用100目筛对研磨后的样本进行过筛；

[0054] 采用HF‑HNO3‑HClO4对样本进行消煮，然后利用等离子体质谱仪测定重金属镉含量；利用碳氮元素分析仪测定土壤有机碳含量；用pH计测定土壤酸碱度。

[0055] 可选地，所述获取混合光谱数据，包括：

[0056] 利用调查区域数字高程模型以及地面控制点经纬度坐标对遥感影像进行正射校正，获取原始光谱数据；

[0057] 对原始光谱数据进行辐射定标以及大气校正，将光谱数据单位转换至地表反射率；

[0058] 提取采样点对应位置遥感影像像元的原始光谱数据，利用利用一阶微分、倒数对数的光谱变换方法对原始光谱数据进行光谱变换，分别得到一阶微分变换结果和倒数对数变换结果；

[0059] 分别计算一阶微分变换结果、倒数对数变换结果、原始光谱数据各个波段与重金属浓度之间的相关系数；

[0060] 对比同一波段不同处理方式下光谱数据与重金属浓度之间相关性的大小，选用每一波段中相关系数最大的光谱数据组成混合光谱数据。

[0061] 可选地，所述获取环境变量数据，包括：

[0062] 根据调查区域的数字高程模型，通过ArcGIS软件生成表征调查区域地形的栅格数据，通过SAGA GIS软件生成地形湿度指数数据；其中，所述栅格数据包括坡度、坡向、平面曲率、剖面曲率；

[0063] 若调查区域的路网和水系的矢量数据的质量不满足预设要求，则对所述矢量数据进行勾画修正，然后根据调查区域的范围生成栅格数据，并利用ArcGIS软件计算栅格数据中每个像元离最近的河流或者道路的距离，并将该距离的值赋值给对应的像元，生成用来表征交通、水流影响程度的栅格数据；

[0064] 利用监督分类方法对遥感影像进行分类，将遥感影像分为不同地物类型，其中，地物类型包括但不限于农田、人造地表、水体及森林；提取调查区域人造地表，生成表征人类活动的栅格数据；

[0065] 若调查区域周边存在重金属污染源，则根据调查区域的范围生成栅格数据，并基于污染源位置，将栅格数据中每个像元离污染源的距离赋值给该像元，生成表征离污染源距离的栅格数据。

[0066] 可选地，所述方法还包括从所述环境变量数据中确定重要环境变量的过程，具体包括：

[0067] 将所述环境变量数据作为因子，将土壤重金属浓度作为变量，通过地理探测器模型中的分异及因子探测器进行重要性分析，获取每个环境变量对于重金属浓度分布的影响力，并依据影响程度从高到低进行排序，得到第一序列；

[0068] 将所述环境变量数据作为因子，将土壤重金属浓度作为变量，通过随机森林重要性分析的方法进行重要性分析，获取每个环境变量对于重金属浓度分布的影响力，并依据影响程度从高到低进行排序，得到第二序列；

[0069] 从所述第一序列和所述第二序列中选取目标个数的环境变量作为重要环境变量。

[0070] 可选地，所述根据所述土壤样本数据、混合光谱数据和环境变量数据，构建重金属反演模型，包括：

[0071] 根据采样点的经纬度坐标，提取采样点对应位置的影响因子数据，构建数据集；

[0072] 基于十折交叉验证方法，将所述数据集随机分为10等份，选用其中9份样本作为训练集，利用偏最小二乘法、Bp神经网络、支持向量机、随机森林四种建模方法，以混合光谱数据及重要影响因子数据为自变量，以土壤重金属浓度数据为因变量训练模型；

[0073] 选用另外1份样本为验证集，对四种模型进行精度评价，计算模型的决定系数、平均绝对误差、均方根误差；

[0074] 重复对模型进行训练，获取满足训练条件的模型构建得到重金属反演模型。

[0075] 可选地，所述根据所述重金属反演模型对目标区域进行土壤重金属含量反演，确定所述目标区域的重金属分布图，包括：

[0076] 根据混合光谱筛选结果，生成目标区域范围内混合光谱栅格数据；

[0077] 根据目标区域范围生成栅格数据，以混合光谱栅格数据以及栅格形式的重要影响因子为自变量，选用所述重金属反演模型反演栅格中所有像元的重金属浓度，得到所述目标区域的重金属浓度空间分布图；

[0078] 提取所述目标区域的农田区域，根据所述农田区域对所述重金属浓度空间分布图进行裁剪，得到所述目标区域的重金属分布图。

[0079] 本发明实施例的另一方面还提供了一种电子设备，包括处理器以及存储器；

[0080] 所述存储器用于存储程序；

[0081] 所述处理器执行所述程序实现如前面所述的方法。

[0082] 本发明实施例的另一方面还提供了一种计算机可读存储介质，所述存储介质存储有程序，所述程序被处理器执行实现如前面所述的方法。

[0083] 本发明实施例的另一方面还提供了一种计算机程序产品，包括计算机程序，所述计算机程序被处理器执行时实现如前面所述的方法。

[0084] 下面结合说明书附图，对本发明的具体实现过程进行详细描述：

[0085] 如图1所示，本发明包括以下步骤：步骤一：农田采样点土壤样本数据获取。步骤二：混合光谱数据获取。步骤三：影响土壤镉浓度空间分布的环境变量数据生产。步骤四：重要环境变量确定。步骤五：最优重金属反演模型的建立。步骤六：区域尺度农田土壤重金属含量反演。

[0086] 下面具体描述各个步骤的实现过程：

[0087] 步骤一：农田采样点土壤样本数据获取。具体细化步骤如下Step1‑3：

[0088] Step1：样本采集：采用五点采样法采集农田0‑20cm表层土壤，并利用GPS记录中心点的经纬度坐标。采样点应均匀分布于调查区域。

[0089] Step2：样本预处理：将样本置于干燥通风处使其自然风干，并在室内将样本内的稻茬、杂草等杂质剔除。之后将样本进行研磨，过100目筛。

[0090] Step3：样本成分检测：采用HF‑HNO3‑HClO4消煮，并利用电感耦合等离子体质谱仪(ICP‑MS，Thermo Fisher，USA)进行测定重金属镉(Cd)含量；利用碳氮元素分析仪测定土壤有机碳含量；利用pH计测定土壤酸碱度，土水比为1:2.5。

[0091] 步骤二：混合光谱数据获取。具体细化步骤如下,包括Step1‑5：

[0092] Step1：光谱预处理：对原始光谱数据进行辐射定标以及大气校正，将光谱数据单位转换至地表反射率。利用调查区域DEM以及地面控制点经纬度坐标对遥感影像进行正射校正，提高遥感影像几何精度。其中地面控制点坐标需要在土壤采集时利用GPS实地记录。

[0093] Step2：光谱数据提取：根据步骤一Step1所采集的采样点经纬度坐标，提取采样点对应位置遥感影像像元的原始光谱数据。

[0094] Step3：光谱变换：利用利用一阶微分、倒数对数的光谱变换方法对原始光谱数据进行光谱变换，分别得到一阶微分和倒数对数的变换结果。其中一阶微分的变换公式为倒数对数的变换公式为 式中Dj、Lj代表光谱变换结果。Bj代表第j个波段的原始光谱反射率，Dj表示第j个波段一阶微分后的光谱数据，gj为波段之间窗口宽度。Lj表示第j个波段倒数对数变换后的光谱数据。

[0095] Step4：相关性分析：分别计算一阶微分变换结果、倒数对数变换结果、原始光谱数据各 个波段 与重金 属浓 度之间的 相关 系数 。相关 系数的 计算 公式 为Rj是其中一种变换方式中第j个波段与重金属浓度之间的相关系数。 代表一种变换方式中第j个波段第i个样本的数值， 代表一种变换方式中第ji
个波段所有样本的平均值。C代表第i个样本的重金属浓度，代表所有样本重金属浓度的平均值。

[0096] Step5：对比同一波段不同处理方式(原始光谱、一阶微分、倒数对数)的光谱数据与重金属浓度之间相关性的大小，选用每一波段相关系数最大的光谱数据组成混合光谱数据。

[0097] 步骤三：影响土壤镉浓度空间分布的环境变量数据生产。具体细化步骤如下步骤Step1‑5：

[0098] Step1：根据调查区域的数字高程模型，基于ArcGIS软件生产表征调查区域地形的栅格数据，包括坡度、坡向、平面曲率、剖面曲率，采用SAGA GIS软件生产地形湿度指数数据。

[0099] Step2：若调查区域OpenStreetMap(OSM)路网和水系的矢量数据质量不满足要求，需要在其基础上进行人工勾画。在此基础上，根据调查区域范围生成栅格数据，并利用ArcGIS软件计算栅格中每个像元离最近的河流或者道路的距离，并将该值赋值给该像元，生成用来表征交通、水流影响程度的栅格数据。

[0100] Step3：利用监督分类方法对遥感影像进行分类，根据调查区域具体情况将遥感影像分为农田、人造地表、水体、森林等地物类型，并提取调查区域人造地表。根据调查区域范围生成栅格数据，利用数字图像处理技术，计算栅格中每个像元周边1公里范围内的人造地表面积，并将此面积赋值给该像元，生成表征人类活动的栅格数据。周围人造地表面积的详细计算过程为：

[0101] 1)、将遥感影像的分类结果二值化，人造地表区域赋值为1，非人造地表区域赋值为0；

[0102] 2)、利用MATLAB读入二值化结果，生成数值为0或1的矩阵，并读入二值化结果的地理信息。

[0103] 3)、根据二值化结果像元大小，设置对应于真实地表1公里的窗口。若每个像元的大小为10m*10m，则窗口大小设置为101*101。

[0104] 4)、以此窗口遍历二值化矩阵中的所有元素。以二值化矩阵中的元素为窗口中心，统计被窗口覆盖的数值为1的元素个数。以此来计算此元素周边人造地表面积。若数值为12
的元素个数为100，则周边人造地表面积为10000m。

[0105] 5)、生成与二值化结果同样大小的空矩阵，储存相同位置的周边人造地表面积。

[0106] 6)、根据二值化结果的地理信息，利用MATLAB将储存周围人造地表面积的矩阵转换为带有地理信息的栅格数据。

[0107] Step4：若调查区域周边存在矿业开采、污水排放、废气排放等可能的重金属污染源，需要根据调查区域范围生成栅格数据，基于污染源位置，将栅格中每个像元离污染源的距离赋值给该像元，生成表征离污染源距离的栅格数据。

[0108] Step5：根据步骤一所述的采样点土壤酸碱度、有机碳含量、重金属浓度，利用克里金插值法，获取调查区域范围内的土壤酸碱度、土壤有机碳含量以及土壤重金属浓度的栅格数据。

[0109] 步骤四：重要环境变量确定。具体细化步骤如下Step1‑3：

[0110] Step1：利用地理探测器(geographical detectors model，GDM)模型中的分异及因子探测器，以步骤三所获取的环境变量(DEM、坡度、坡向、平面曲率、剖面曲率、地形湿度指数、离道路距离、离河流距离、周围人造地表面积、离污染源距离、土壤有机碳、土壤酸碱度)作为因子X，以土壤重金属浓度为变量Y进行重要性分析，获取每个环境变量对于重金属浓度分布的影响力，并依据影响程度从高到低进行排序。

[0111] 分异及因子探测可用于探测某因子X多大程度上解释了属性Y的空间分异，用q值度量。

[0112]

[0113]

[0114] SST＝Nσ2

[0115] 式中：h为因子X的分层，即分类或分区(可采用自然断点法对因子X进行分层)。Nh和N分别为层h和全区变量Y的像元数；σh和σ为别为层h和全区变量Y值的方差。SSW和SST分别为层内方差之和和全区总方差。q值越大表示自变量X对属性Y的解释力越强，反之则越弱。极端情况下，q值为1表明因子X完全控制了Y的空间分布，q值为0则表明因子X与Y没有任何关系。

[0116] Step2：根据步骤一Step1所采集的采样点经纬度坐标，提取采样点对应位置的环境变量数据。利用随机森林(random forest，RF)重要性分析功能以环境变量为自变量X，重金属浓度为因变量Y，计算所有环境变量对于重金属浓度分布的重要性，并依据重要性从高到低进行排序。随机森林模型拟合时会通过自助抽样形成大约2/3的样本数据，剩余样本数据称为袋外样本数据(OOB)，用于评估模型整体误差以及变量的重要性。在进行重要性评价时，首先利用自助抽样的样本数据建立回归树模型，计算OOB的预测误差，然后随机置换自变量X的观测值，再次建立回归树模型，计算OOB'(随机置换变量观测后袋外数据)的预测误差，将2次袋外数据预测误差的差值进行标准化处理后，将所有回归树结果取均值即为变量的重要性。若2次袋外数据预测误差的差值较大，则说明随机置换自变量X的观测值后，随机森林的预测精度受到了较大程度的影响，此自变量X对于因变量Y的预测至关重要。本项目通过随机森林重要性排序验证地理探测器因子分析的结果，以不同的方法确定重要的驱动因子，分析驱动因子与重金属浓度之间的相关性。

[0117] Step3：为降低使用单一方法而产生的不确定性，综合选取Step1、Step2排序结果前五项的并集作为重要环境变量参与后续模型建立。

[0118] 步骤五：最优重金属反演模型的建立。具体细化步骤如下：

[0119] Step1：根据步骤一的Step1所采集的采样点经纬度坐标，提取采样点对应位置的影响因子数据。

[0120] Step2：将数据集随机分为10等份，选用其中9份样本作为训练集，利用偏最小二乘法、Bp神经网络、支持向量机、随机森林四种建模方法，以混合光谱数据及重要影响因子数据为自变量，以土壤重金属浓度数据为因变量训练模型。

[0121] Step3：选用另外1份样本为验证集，对四种模型进行精度评价，计算模型的决定系数、平均绝对误差、均方根误差。

[0122] Step4：重复Step1、Step2步骤10次，保证每份样本均作为验证集参与过模型精度评价，记录每次模型评价结果。

[0123] Step5：对10次模型评价计算得到的决定系数、平均绝对误差、均方根误差取平均值，以此作为四种模型的最终精度。对比分析四种模型的精度，选择精度最高的模型参与大面积重金属浓度反演。

[0124] 步骤六：区域尺度农田土壤重金属含量反演。具体细化步骤如下：

[0125] Step1：根据混合光谱筛选结果，生产研究区范围内混合光谱栅格数据。

[0126] Step2：根据研究区范围生成栅格数据，以混合光谱栅格数据以及栅格形式的重要影响因子为自变量，选用最优模型反演栅格中所有像元的重金属浓度，得到调查区域重金属浓度空间分布图。

[0127] Step3：利用步骤三Step 3中的监督分类结果，提取其中的农田区域。用农田区域对Step2获取的重金属分布图进行裁剪，得到调查区域的农田重金属分布图。

[0128] 可以理解的是，本发明通过融合光谱信息以及环境因子进行建模，进而基于最优模型大面积预测农田土壤重金属浓度。

[0129] 综上所述，本发明提供了一种融合光谱信息与环境变量的土壤重金属含量反演方法。该技术具备以下特点：

[0130] 1、本发明利用高空间分辨率WorldView‑3多光谱数据对土壤重金属浓度进行反演建模。WorldView‑3多光谱数据空间分辨率较高，具有16个波段，涵盖可见光、近红外和短波红外谱段，为土壤重金属含量反演提供了较多的光谱信息。

[0131] 2、本发明通过对光谱数据进行光谱变换，并利用相关性分析对光谱数据进行筛选，有效地提高了光谱数据质量。

[0132] 3、本发明引入多源环境变量数据，并对其深入分析，选择对土壤重金属分布具有重要驱动作用的环境因子数据辅助反演模型的建立，进一步提高了模型反演精度。

[0133] 相较于现有技术，本发明对光谱数据进行光谱变换，并通过相关性分析筛选出混合光谱，提高了光谱数据质量。本发明结合自然因子(地形、离河流距离、土壤酸碱度、土壤有机碳)，人为因子(周边人造地面积、离道路距离)，对光谱信息进行补充，提高了反演建模精度。本发明以周边人造地表面积来反映人类活动对此地区土壤性质的影响程度。本发明综合利用地理探测器、随机森林对环境因子重要性进行分析，选择重要影响因子参与建模。本发明的精度高体现在通过筛选混合光谱数据，结合环境影响因子共同参与建模，极大地提高了土壤重金属含量的反演精度。

[0134] 在一些可选择的实施例中，在方框图中提到的功能/操作可以不按照操作示图提到的顺序发生。例如，取决于所涉及的功能/操作，连续示出的两个方框实际上可以被大体上同时地执行或所述方框有时能以相反顺序被执行。此外，在本发明的流程图中所呈现和描述的实施例以示例的方式被提供，目的在于提供对技术更全面的理解。所公开的方法不限于本文所呈现的操作和逻辑流程。可选择的实施例是可预期的，其中各种操作的顺序被改变以及其中被描述为较大操作的一部分的子操作被独立地执行。

[0135] 此外，虽然在功能性模块的背景下描述了本发明，但应当理解的是，除非另有相反说明，所述的功能和/或特征中的一个或多个可以被集成在单个物理装置和/或软件模块中，或者一个或多个功能和/或特征可以在单独的物理装置或软件模块中被实现。还可以理解的是，有关每个模块的实际实现的详细讨论对于理解本发明是不必要的。更确切地说，考虑到在本文中公开的装置中各种功能模块的属性、功能和内部关系的情况下，在工程师的常规技术内将会了解该模块的实际实现。因此，本领域技术人员运用普通技术就能够在无需过度试验的情况下实现在权利要求书中所阐明的本发明。还可以理解的是，所公开的特定概念仅仅是说明性的，并不意在限制本发明的范围，本发明的范围由所附权利要求书及其等同方案的全部范围来决定。

[0136] 所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read‑Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。

[0137] 在流程图中表示或在此以其他方式描述的逻辑和/或步骤，例如，可以被认为是用于实现逻辑功能的可执行指令的定序列表，可以具体实现在任何计算机可读介质中，以供指令执行系统、装置或设备(如基于计算机的系统、包括处理器的系统或其他可以从指令执行系统、装置或设备取指令并执行指令的系统)使用，或结合这些指令执行系统、装置或设备而使用。就本说明书而言，“计算机可读介质”可以是任何可以包含、存储、通信、传播或传输程序以供指令执行系统、装置或设备或结合这些指令执行系统、装置或设备而使用的装置。

[0138] 计算机可读介质的更具体的示例(非穷尽性列表)包括以下：具有一个或多个布线的电连接部(电子装置)，便携式计算机盘盒(磁装置)，随机存取存储器(RAM)，只读存储器(ROM)，可擦除可编辑只读存储器(EPROM或闪速存储器)，光纤装置，以及便携式光盘只读存储器(CDROM)。另外，计算机可读介质甚至可以是可在其上打印所述程序的纸或其他合适的介质，因为可以例如通过对纸或其他介质进行光学扫描，接着进行编辑、解译或必要时以其他合适方式进行处理来以电子方式获得所述程序，然后将其存储在计算机存储器中。

[0139] 应当理解，本发明的各部分可以用硬件、软件、固件或它们的组合来实现。在上述实施方式中，多个步骤或方法可以用存储在存储器中且由合适的指令执行系统执行的软件或固件来实现。例如，如果用硬件来实现，和在另一实施方式中一样，可用本领域公知的下列技术中的任一项或他们的组合来实现：具有用于对数据信号实现逻辑功能的逻辑门电路的离散逻辑电路，具有合适的组合逻辑门电路的专用集成电路，可编程门阵列(PGA)，现场可编程门阵列(FPGA)等。

[0140] 在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

[0141] 尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。

[0142] 以上是对本发明的较佳实施进行了具体说明，但本发明并不限于所述实施例，熟悉本领域的技术人员在不违背本发明精神的前提下还可做出种种的等同变形或替换，这些等同的变形或替换均包含在本申请权利要求所限定的范围内。