基于室内标样高光谱特征的土壤重金属Cd含量反演建模及其光谱响应特征波段识别方法与流程

本发明涉及土壤重金属含量高光谱反演领域，特别涉及基于室内标样高光谱特征的土壤重金属cd含量反演建模及其光谱响应特征波段识别的方法。

背景技术：

随着城市化、工业化和农业集约化的发展，农用地土壤重金属污染物因残留时间长、隐蔽性强、迁移性小、毒性大、不易被微生物降解、可以通过食物链富集等特点，逐渐受到国内外研究者的普遍关注。然而，因大范围高密采样成本耗费巨大、常态化快速高效动态监测难以实施，大范围开展重金属污染治理至今仍然存在“重金属污染安全隐患点不清、污染空间格局不明”两大突出问题，已成为摆在土壤重金属污染全面与深度治理面前的一道“理论与方法”鸿沟，无法满足国家在广域地理空间快速高效查清重金属污染空间格局、发现流域重金属污染时空演化规律与形成机制、控制污染范围扩大和转移、合理规划农业生产、减少污染产品损害国民健康的重大需求。

与此同时，新兴的遥感探测技术具有能通过非接触的方式对目标物体进行快速、动态、大范围监测的特征，高光谱更是因其多且连续的光谱波段优势逐渐被广泛应用于土壤重金属污染监测研究中，为定量反演土壤重金属含量提供了新的契机，使得大范围、高时空分辨率、高精度快速探测土壤重金属污染靶区的含量空间分布成为可能。

近年来，国内外围绕土壤重金属含量高光谱反演开展了大量的研究工作。截止目前，大量研究已证实：当土壤重金属含量达到一定数值时，在地面可见光、近红外和热红外波段光谱测试中可发现明显的重金属光谱响应特征；在土壤重金属含量偏低时，仍可通过地面土壤有机质、铁猛氧化物和粘土矿物的光谱特征间接推算重金属含量。然而，土壤成分复杂，重金属含量甚微，有效排除土壤水分、有机质、铁锰氧化物等成分光谱干扰，建立重金属含量反演模型，准确识别微量重金属元素光谱响应特征波段，是亟待解决的关键技术难题。

技术实现要素：

为了有效排除土壤水分、有机质、铁锰氧化物等成分光谱干扰，准确识别微量重金属元素光谱响应特征波段，本发明提供一种基于室内标样高光谱特征的土壤重金属cd含量反演建模及其光谱响应特征波段识别的方法。

为了实现上述技术目的，本发明的技术方案是，一种基于室内标样高光谱特征的土壤重金属cd含量反演建模及其光谱响应特征波段识别的方法，包括以下步骤：

步骤1：土壤标样制作及其准确重金属cd含量与土壤光谱数据获取；

步骤2：土壤重金属cd含量分析与光谱数据预处理；

步骤3：基于步骤2转换得到的光谱数据的土壤重金属cd含量反演建模，实现光谱响应特征波段识别。

所述的方法，步骤1中所述的土壤标样的制作及准确其重金属cd含量的获取，包括以下步骤：

步骤1)：充分考虑土壤重金属cd污染含量实际变化范围及其光谱响应效果，人为设计不同cd污染含量(>50等级)(单位：mg/kg)下的标准土壤样本；

步骤2)：采用格网采集方式以1m×1m分辨率采集研究区未经重金属污染的背景土壤样本若干(>50份)，留作标样添加土壤底样；

步骤3)：将步骤2)采集到的背景土壤经实验室烘干、研磨处理，准确称量得到标样添加土壤底样，每份100g；

步骤4)：取步骤3)处理后的底样一份，采用三酸消化——原子吸收分光光度法实验室化学分析测定背景土壤重金属cd含量(单位：mg/kg)；

步骤5)：计算达到步骤1)中设定的cd污染含量值需向100g底样中添加的1000mg/kg的cd标准溶液，计算公式如下：

c＝100(a-s)/ρ

式中，c代表加入100g底样土壤中的cd标准溶液量(单位：ml)，a代表步骤1)中设计的土壤标样的重金属cd含量(单位：mg/kg)，s代表步骤4)中测得的背景土壤样本的重金属cd含量，ρ代表cd标准溶液的密度(单位：kg/m³)。

特别地，向100g土壤底样中人为添加cd标准溶液后，需经风干处理后测定土壤重金属cd含量，标样cd含量最高不超过0.1mg，为底样的10^-4数量级，可忽略不计，因此土壤总量按100g计。

步骤6)：实验室按步骤5)中计算得到的cd标准溶液量依次添加到步骤3)获取的土壤底样中，制得土壤标准样本；

步骤7)：采用三酸消化一一一原子吸收分光光度法实验室化学分析测定步骤6)中制得的土壤标样，准确获取土壤标样重金属cd含量。

所述的方法，步骤2中所述的土壤污染标样光谱与背景样本光谱作比值的光谱预处理方法，包括以下步骤：

步骤1)：将实验室人为添加cd污染标样制作的土壤标样光谱与未经人为添加重金属cd标准溶液的土壤背景样本光谱同时做10nm间隔重采样—标准正态变换(snv)—一阶/二阶微分处理，得到转换后的光谱数据；

步骤2)：将经步骤1)光谱变换后的土壤污染标样光谱与背景样本光谱作比值，得到光谱特征数据，计算公式如下：

ri＝ai/bi(i＝1,2,…,n)

式中，n代表光谱波段数，i代表第i个波段，ai代表第i个波段的步骤1)光谱转换后的土壤污染样本光谱反射率，bi代表第i个波段的步骤1)光谱转换后的土壤光谱反射率，ri代表第i个波段经比值变换后得到的光谱反射率，作为土壤重金属cd含量高光谱反演建模输入变量光谱特征数据。

附图说明

图1为实验室标样制作方法流程；

图2为光谱基于各类光谱预处理方法的光谱变换效果实例；

图3为偏最小二乘法建模因子提取过程的基本原理；

图4为偏最小二乘法建模最小平方误差理论原理。

具体实施方式

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

1、如图1流程图所示，土壤标样制作及其准确重金属cd含量与土壤光谱数据获取，具体原理方法如下：

步骤1)：充分考虑土壤重金属cd污染含量实际变化范围及其光谱响应效果，人为设计不同cd污染含量(>50等级)(单位：mg/kg)下的标准土壤样本；

步骤2)：采用格网采集方式以1m×1m分辨率采集研究区未经重金属污染的背景土壤样本若干(>50份)，留作标样添加土壤底样；

步骤3)：将步骤2)采集到的背景土壤经实验室烘干、研磨处理，准确称量得到标样添加土壤底样，每份100g；

步骤4)：取步骤3)处理后的底样一份，采用三酸消化——原子吸收分光光度法实验室化学分析测定背景土壤重金属cd含量(单位：mg/kg)；

步骤5)：计算达到步骤1)中设定的cd污染含量值需向100g底样中添加的1000mg/kg的cd标准溶液，计算公式如下：

c＝100(a-s)/ρ

式中，c代表加入100g底样土壤中的cd标准溶液量(单位：ml)，a代表步骤1)中设计的土壤标样的重金属cd含量(单位：mg/kg)，s代表步骤4)中测得的背景土壤样本的重金属cd含量，ρ代表cd标准溶液的密度(单位：kg/m³)。

特别地，向100g土壤底样中人为添加cd标准溶液后，需经风干处理后测定土壤重金属cd含量，标样cd含量最高不超过0.1mg，为底样的10^-4数量级，可忽略不计，因此土壤总量按100g计。

步骤6)：实验室按步骤5)中计算得到的cd标准溶液量依次添加到步骤3)获取的土壤底样中，制得土壤标准样本；

步骤7)：采用三酸消化一一一原子吸收分光光度法实验室化学分析测定步骤6)中制得的土壤标样，准确获取土壤标样重金属cd含量。

2、采用光谱重采样—标准正态变换(snv)—一阶/二阶微分—背景样比值系列光谱变换方法进行光谱数据预处理，具体原理方法如下：

步骤1)：为消除数据冗余，降低噪声，平滑波形，得到平稳的变化波谱，以10nm为间隔进行光谱数据的重采样，计算公式如下：

r＝(r1+r2+…+rn)/n

式中，r代表重采样后的光谱反射率；n代表该10nm间隔内波段数，rn代表第n个波段的光谱反射率。

步骤2)：对重采样后的光谱采用加权平均化方法进行标准正态变量变换(snv)。对于给定的样本，snv计算了所有变量的标准差，整个样本再通过该值进行归一化，计算公式如下：

式中，n是变量的个数；xi,j是第i个样本的第j个变量的值；δ是自定义的偏移量。

步骤3)：采用一阶/二阶微分处理去除部分线性或接近线性的背景值干扰，即噪声光谱对目标光谱的影响，增强光谱特征差异，提取光谱特征吸收带，光谱微分计算公式为：

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

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

式中，λi为每个波段的光谱反射率，r′(λi)和r″(λi)分别为波长λi的一阶和二阶微分光谱，δλ是波长λi-1到λi的间隔，本研究间隔经重采样变换均为10nm。

步骤4)：将经重采样—snv—一阶/二阶微分光谱变换后的土壤污染标样光谱与背景样本光谱作比值，得到光谱特征数据，计算公式如下：

ri＝ai/bi(i＝1,2,…,n)

式中，n代表光谱波段数，i代表第i个波段，ai代表第i个波段的步骤1)光谱转换后的土壤污染样本光谱反射率，bi代表第i个波段的步骤1)光谱转换后的土壤光谱反射率，ri代表第i个波段经比值变换后得到的光谱反射率，作为土壤重金属cd含量高光谱反演建模输入变量光谱特征数据。

3、基于数据预处理光谱转换得到的光谱数据采用偏最小二乘回归(plsr)方法构建土壤重金属含量反演模型，初步识别出土壤重金属含量反演的特征波段土壤重金属cd含量高光谱反演建模的步骤包括：

步骤1)：设光谱变换后得到的波段数为m，土壤样本数量为n，t，u分别为从自变量与因变量中提取的因子，称为偏最小二乘因子。从变换后得到的各波段光谱反射率(x1,x2,…,xm)中提取相互独立的成分th(h＝1,2,…)，从各土壤样本重金属cd含量值(y1,y2,…,yn)中提取相独立的成分uj(j＝1,2,…)，要求th和uj之间的相关程度达到最大。从原始变量集中提取第一对因子t1和u1的线性组合公式为：

式中，w1′＝(w1,1,w1,2,…,w1,m)′为模型效应权重；v1′＝(v1,1,v1,2,…,v1,m)′为因变量权重。

步骤2)：建立初始变量对t1的方程。

其中，α′1＝(α1,1,α1,2,…,α1,m)′；β1′＝(β1,1,β1,2,…,β1,m)′为仅有一个自变量t1时的参数向量，其中α′1为模型效应载荷量。e1和f1分别为残差阵。根据普通最小二乘法可求得参数向量α1，β1。

步骤3)：利用多元回归方法建立提取的成分与因变量土壤重金属含量值之间的回归方程。若提取的第一因子不能达到回归模型的精度要求，则利用残差阵e1和f1相应取代x0和y0，重复步骤2)，步骤3)继续提取因子，以此类推。假设最终提取了r个因子，x0和y0对r个因子的回归方程为：

把第一步分析所得的自变量中提取因子th(h＝1,2,…r)的线性组合带入因变量对r个因子建立的回归方程，即得因变量的回归方程，即为土壤重金属含量反演模型，具体公式为：

y＝α1x1+α2x2+···+αmxm

步骤3)：提取上述模型得到的与土壤重金属含量显著相关的波段(xs1,xs2,…,xsm)，即为光谱响应特征波段。