技术领域
[0001] 本
发明涉及土壤中重金属污染检测领域，具体的说是一种基于无人机高光谱反演土壤中重金属污染监测方法。
背景技术
[0002] 土壤重金属污染是指由于人类活动，土壤中的微量金属元素在土壤中的含量超过背景值，过量沉积而引起的含量过高，统称为土壤重金属污染。土壤重金属是指由于人类活动将金属加入到土壤中，致使土壤中重金属明显高于原生含量、并造成生态环境
质量恶化的现象。农田土壤作为人类各种粮食作物赖以生存的自然环境，其生态环境问题一直是环境科学界研究的热点。随着世界工业化和城市化
进程整体快速推进，环境污染日益加剧，土壤的重金属污染甚为严重，为了更好地对土壤环境进行监测评估和控制，研究和开发有效的检测土壤重金属方面的仪器和监测方法十分必要。
[0003] 依据目前的研究现状，农田土壤重金属污染的监测仍以现场
采样实验室分析为主。但是，随着对农田土壤重金属的监测机理及其影响监测的因素不断深入的研究总结、监测分析技术不断的改进和完善，必然使土壤重金属污染监测朝着准确、快速、原位、实时、经济、定量化及易操作的现场快速监测的方向发展。
[0004] 目前对于土壤重金属含量的监测主要是通过传统的野外采样和室内化学分析方法获得，其具有测量
精度高、准确性强等优点，但相对费时费
力，而且很难获取大面积空间上连续的污染物含量分布信息。遥感技术因其多时相、大面积等特点逐渐被研究者应用于土壤性质的监测，高光谱遥感则以其多且连续的光谱波段特点被应用于监测土壤重金属含量，可以实现大范围、非侵入性和非
接触性的快速测样。但是由于重金属在土壤中的含量一般比较微量，因此其对光谱曲线的影响十分细微，直接分析光谱曲线无法得到较精确的污染情况。
发明内容
[0005] 为了克服
现有技术中的不足，本发明将提供一种基于高光谱技术的大范围、非侵入性和非接触性的快速测样的重金属污染监测方法。
[0006] 为了实现上述功能，本发明将采用以下技术方案：
[0007] 第一步，现场采样：在研究区以平均分布原则随机选取采样点，根据采样点面积、地势、土壤组成和受污染程度选择包括蛇形、梅花桩式、对
角线式、棋盘式等不同的采样方式。其中，对角线式采样法具体为对角线分5等分，以等分点为采样分点，适用于
灌溉农田土壤；梅花桩式采样法具体为：5个左右样本，适用于面积较小，地势平坦，土壤组成和受污染程度相对比较均匀的
地块；棋盘式采样法具体为：10个～20个左右样本，适宜中等面积、地势平坦、土壤不够均匀的地块和受
污泥、垃圾等固体废物污染的土壤；蛇形采样法具体为：15个左右样本，适宜于面积较大、土壤不够均匀且地势不平坦的地块，多用于农业污染型土壤。
[0008] 第二步，对样本进行预处理：
风干样本，去除土壤中贝壳、
杂草、树根等有机质，并
研磨至过2mm筛；
[0009] 第三步，采集样本的重金属含量：采用试样制备简单，可多样品、多元素同时分析，分析速度快，重现性好和非破坏性测定的
X射线荧光光谱法(XRF)对样本进行分析。具体操作为将研磨好的土壤样品装入专用塑料容器中，用
薄膜封装；装样后放入仪器配套架构中，薄膜一侧对准X射线发射窗口，每个样本测5次，每次旋转90°，多次数据取平均值。本发明使用X射线荧光分析仪(Thermo Scientific Niton XL3t GOLDD+)采集样本的多种元素含量，该仪器直接获取元素含量，一次能获得多种元素，从中选取重金属污染源，可以从中选择主要研究元素，一般来讲会选取Pb、Zn、Cu、Cr、Ni作为主要研究元素；
[0010] 第四步，采集样本的实验室高光谱反射率：将研磨好的土壤置于直径6cm高1cm的玻璃皿中，利用地物光谱仪(Spectral Evolution PSR+3500)进行重复测量，取平均得到样本在350～2500nm的光谱信息，具体操作：在暗室中进行，使用
卤素灯(全谱段
辐射源)模拟太阳灯光，控制入射角45°，光谱仪
探头与样本距离15cm，每个样本测4次，每次旋转90°，四次数据取平均值。进行白板定标后，仪器自动除白板去暗流获得样本高光谱反射率信息；
[0011] 第五步，
数据处理：为了增强光谱差异性，减少因光照条件变化引起的乘性因素影响，分解和识别重叠混合光谱，扩大样本之间的光谱特征差异，对原始光谱反射率数据进行倒数(1/R)，对数(LogR)，一阶微分(R’)，二阶微分(R”)变换,通过上述变换，能够突出光谱曲线中变化强烈的地方，从而突出光谱曲线特征；变换后利用九点加权平均法进行去噪处理，计算公式如下：
[0012] R'i＝0.04Ri-4+0.08Ri-3+0.12Ri-2+0.16Ri-1+0.20Ri
[0013] +0.16Ri+1+0.12Ri+2+0.08Ri+3+0.04Ri+4
[0014] 式中：R'i为为该样本经平滑处理后的光谱反射率值，各项系数
[0015] 为以i点为中心点的9个点的不同权重值。
[0016] 第六步，建模：使用偏最小二乘回归
算法分别将已经测得的主要研究元素的含量分别与实验室高光谱原始反射率数据、倒数变换后的原始光谱反射率数据、对数变换后的原始光谱反射率数据、一阶微分变换后的原始光谱反射率数据以及二阶微分变换变换后的原始光谱反射率数据进行相关性分析，建立待测金属含量与实验室高光谱曲线的反演模型；
[0017] 建模的具体步骤如下：
[0018] 设X0为PSR测得的高光谱曲线，共n份样本，m个波段，Y0为XRF测得的元素含量，共p个元素。
[0019] 现以一个元素为例即p＝1，假设n＝42,m＝2150,
[0020]
[0021] 1.将X0，Y0进行标准化(减均值，除标准差等)得到X(n×m)，Y(n×p)以对样本进行归一化，同一量级的两个变量的回归速度更快；
[0022] 2.设X的第一个主成分轴分量为p1，Y的第一个主成分轴分量为q1，并对两者进行单位化，(注意，此处的主成分并不是通过主成分分析(PCA)得出的，是需要求的未知向量)。
[0023] 注：由于每个样本X包含m个波段即m个特征，特征数远远大于样本数导致X矩阵为奇异矩阵即不可逆矩阵，不便于后续的变换，因此，在此引入主成分轴分量以降低特征数。
[0024] 3.令u1＝X·p1，v1＝Y·q1，即将原始数据投影至主成分轴上
[0025] 4.代价函数1：Var(u1)→max,Var(v1)→max，即主成分轴上的投影，方差最大化；代价函数2：Corr(u1,v1)→max，即两者相关性最
[0026] 注：偏最小二乘回归结合了主成分分析法和最小二乘回归，因此代价函数(cost function)也要集合两种方法。
[0027] 5.综合一下两个条件，得到优化目标
[0028]
[0029] 注：同时考虑两个代价函数在数学上变现为两者相乘
[0030] i.由此，该问题就转化为一个优化问题：将一个实际问题抽象为数学问题，方便转化为机器语言，
[0031] Maximize〈Xp1,Yq1〉
[0032] ii.Subject to:||p1||＝1,||q1||＝1
[0033] 优化问题的求解：
[0034] 1.引入拉格朗日乘子：
[0035] 分别对p1，q1求偏导：
[0036]
[0037] 目标函数 即要求λ最大
[0038] 即得到：p1是对称阵XTYYTX最大特征值对应的单位
特征向量，q1是对称阵YTXXTY最大特征值对应的单位特征向量
[0039] 2.求得了p1，q1，即可得到u1＝X·p1，v1＝Y·q1，然后对其主成分进行回归建模其中c1，d1与p1，q1不同，E和G是残差矩阵
[0040] 3.利用X的主成分对Y进行回归，即用u1对Y进行回归：Y＝u1r1T+F
[0041] 4.通过最小二乘法得到c，d，r分别是： 即得到各个投影向量
[0042] 5.已知u1＝X·p1， 其中p1是X投影出u1的方向向量，而c1则是在回归思想下根据最小二乘法求出的，两者之间一般不相同。
[0043] 6.将X中主成分u1不能解释的残差部分E当做新的X，回归方程Y中主成分u1不能解释的残差部分F当做新的Y：
[0044] 重复上述步骤进行循环回归，直到残差F达到精度要或者主成分数量已经达到上限(初始X的秩)，算法结束
[0045] 7.最后得到：
[0046]
[0048] 其中
[0049] 写成矩阵形式：
[0050] X＝UCT+E
[0051] Y＝URT+F＝XPRT+F≈XB
[0052] 即X→Y的回归方程，其中B＝PRT。
[0053] 8.从7中可以了解Y其实是多个回归的
叠加( 已经回归出Y的最主要信息)，计算模型的过程中我们得到了p和r，当使用回归模型进行预测的时候，新来一个x(新的光谱曲线)，首先计算u，得到u1＝xTp1,u2＝xTp2,u3＝xTp3...，然后代入Y＝URT即可求出预测的Y(即重金属含量)；
[0054] 建模后，对模型进行验证优化，使用留一法作为验证优化方法，如共n个样本，则进行n次回归分析，每次使用n-1个样本作为训练样本，剩下的1个作为检验样本，取n次分析得到的R2(拟合系数)和RMSE(均方根误差)的平均值 和 作为模型检验标准，以 最大， 最小为优化目标对模型的预测能力和鲁棒性进行检验和优化；
[0055] 第七步，大面积反演重金属含量：使用六旋翼无人机(Aibotix GmbH Aibot X6)搭载高光谱成像光谱仪(Headwall公司的Nano-Hyperspec VNIR高光谱成像仪)采集的研究区高光谱数据作为待测数据。将采集到的无人机高光谱数据通过最优变换方法进行数据变换，利用反演土壤重金属含量的高光谱估算模型计算待测研究区的重金属含量，得到大面积重金属含量图像，并对其中含量突出区域进行提取，突出区域根据研究区的属性如农田、养殖区、工业区、林牧区等，依照国家环境土壤质量标准文件GB15618-2009选取
阈值，以判定当地土壤是否超标以及超标面积等参量。
[0056] 第八步，多时相分析：在一个重金属防治周期或当地污染源排放周期后再次在研究区以平均分布原则随机选取采样点，在每个采样点采集多个样本，重复步骤第二步至第七步，得到最优的反演模型，对研究区进行无人机高光谱信息采集，同时进行反演获得本次待测研究区的重金属含量，将本次待测研究区的重金属含量与第八步获取的待测研究区的重金属含量进行对比，以获取防治效果判别或重金属扩散情况，为需求方提供决策依据。
[0057] 与卫星遥感、航空航天遥感相比，无人机有着独特的优势：第一，对场地要求低，作业方式灵活快捷，能快速响应拍摄任务；第二，平台构建、维护以及作业的成本极低，与卫星遥感、航空航天遥感相比几可不计；第三，因其飞行高度低，能够获取大比例尺高精度影像，在局部信息获取方面有着巨大的优势；第四，飞行高度一般低于1000米，不必
申请空域：第五，能够获取高重叠度的影像，增强了后续处理的可靠性：第六，具有便于携带，转移方便等优点。
[0058] 本发明的目的在于克服现有探测技术的不足和光谱曲线对重金属含量的响应过于微弱，建立待测金属含量与实验室高光谱曲线的反演模型，通过无人机与高光谱成像光谱仪结合采集待测光谱信息后通过模型进行大规模反演，得到大面积重金属含量图像。同时，通过多时相分析，可以对获取防治效果判别或重金属扩散情况进行判断，为需求方提供可靠的决策依据。本发明提不仅能够监测某地区的重金属含量是否超标，还能通过多时相的分析获得重金属污染的扩散情况以及防治效果。本发明具有可探测范围大、非侵入性和非接触性的快速测样等优点。
附图说明
具体实施方式
[0060] 下面将结合附图对本发明作进一步的说明：
[0061] 第一步，现场采样：在研究区以平均分布原则随机选取采样点，根据采样点面积、地势、土壤组成和受污染程度选择包括蛇形、梅花桩式、对角线式、棋盘式等不同的采样方式。其中，对角线式采样法具体为对角线分5等分，以等分点为采样分点，适用于灌溉农田土壤；梅花桩式采样法具体为：5个左右样本，适用于面积较小，地势平坦，土壤组成和受污染程度相对比较均匀的地块；棋盘式采样法具体为：10个～20个左右样本，适宜中等面积、地势平坦、土壤不够均匀的地块和受污泥、垃圾等固体废物污染的土壤；蛇形采样法具体为：15个左右样本，适宜于面积较大、土壤不够均匀且地势不平坦的地块，多用于农业污染型土壤。
[0062] 第二步，对样本进行预处理：风干样本，去除土壤中贝壳、杂草、树根等有机质，并研磨至过2mm筛；
[0063] 第三步，采集样本的重金属含量：采用试样制备简单，可多样品、多元素同时分析，分析速度快，重现性好和非破坏性测定的X射线荧光光谱法(XRF)对样本进行分析。具体操作为将研磨好的土壤样品装入专用塑料容器中，用薄膜封装；装样后放入仪器配套架构中，薄膜一侧对准X射线发射窗口，每个样本测5次，每次旋转90°，多次数据取平均值。本发明使用X射线荧光分析仪(Thermo Scientific Niton XL3t GOLDD+)采集样本的多种元素含量，该仪器直接获取元素含量，一次能获得多种元素，从中选取重金属污染源，可以从中选择主要研究元素，一般来讲会选取Pb、Zn、Cu、Cr、Ni作为主要研究元素；
[0064] 第四步，采集样本的实验室高光谱反射率：将研磨好的土壤置于直径6cm高1cm的玻璃皿中，利用地物光谱仪(Spectral Evolution PSR+3500)进行重复测量，取平均得到样本在350～2500nm的光谱信息，具体操作：在暗室中进行，使用卤素灯(全谱段辐射源)模拟太阳灯光，控制入射角45°，光谱仪探头与样本距离15cm，每个样本测4次，每次旋转90°，四次数据取平均值。进行白板定标后，仪器自动除白板去暗流获得样本高光谱反射率信息；
[0065] 第五步，用matlab进行数据处理：对原始光谱反射率数据进行倒数(1/R)，对数(LogR)，一阶微分(R’)，二阶微分(R”)变换,去噪处理，小波分析。
[0066] 第六步，建模：使用matlab建立plsr模型对高光谱数据进行反演，反演之后用matlab采用留一法作为验证优化方法，获得最优变换方法；
[0067] 第七步，大面积反演重金属含量：使用六旋翼无人机(Aibotix GmbH Aibot X6)搭载高光谱成像光谱仪(Headwall公司的Nano-Hyperspec VNIR高光谱成像仪)采集研究区高光谱数据作为待测数据。将采集到的无人机高光谱数据经matlab通过最优变换方法进行数据变换，利用反演土壤重金属含量的高光谱估算模型计算待测研究区的重金属含量，得到大面积重金属含量图像，并对其中含量突出区域进行提取，突出区域根据研究区的属性如农田、养殖区、工业区、林牧区等，依照国家环境土壤质量标准文件GB15618-2009选取阈值，以判定当地土壤是否超标以及超标面积等参量。
[0068] 第八步，多时相分析：在一个重金属防治周期或当地污染源排放周期后再次在研究区以平均分布原则随机选取采样点，在每个采样点采集多个样本，重复步骤第二步至第七步，得到最优的反演模型，对研究区进行无人机高光谱信息采集，同时进行反演获得本次待测研究区的重金属含量，将本次待测研究区的重金属含量与第八步获取的待测研究区的重金属含量进行对比，以获取防治效果判别或重金属扩散情况，为需求方提供决策依据。
[0069] 本发明的目的在于克服现有探测技术的不足和光谱曲线对重金属含量的响应过于微弱，建立待测金属含量与实验室高光谱曲线的反演模型，通过无人机与高光谱成像光谱仪结合采集待测光谱信息后通过模型进行大规模反演，得到大面积重金属含量图像。同时，通过多时相分析，可以对获取防治效果判别或重金属扩散情况进行判断，为需求方提供可靠的决策依据。本发明提不仅能够监测某地区的重金属含量是否超标，还能通过多时相的分析获得重金属污染的扩散情况以及防治效果。本发明具有可探测范围大、非侵入性和非接触性的快速测样等优点。
