一种利用土壤高光谱反射率计算土壤全盐含量的方法与流程

本发明涉及土壤成份测试技术领域，具体涉及一种利用土壤高光谱反射率计算土壤全盐含量的方法。

背景技术：

我国有5亿亩盐碱土地，且分布广泛，合理利用盐渍化耕地、提升盐碱地开发治理与管控水平，对保障区域资源安全、粮食安全、生态安全具有重大战略意义。而土壤盐分信息监测和获取是合理制定土壤盐碱化改良措施的前提。传统田间采样及室内分析方法比较费时、费力，特别是对于深层土壤样品的采集和大范围快速测量更加不易，在土壤盐渍化大范围快速监测中受到限制。由于不同盐分含量的土壤表现出不同的光谱信息特征，而地物高光谱传感器能获取监测对象纳米级的连续光谱信息，因此可以通过对含盐土壤的高光谱定量监测和分析获取其盐分含量。

技术实现要素：

针对现有技术的上述不足，本发明提供了一种可以对区域土壤盐分进行快速监测和诊断的利用土壤高光谱反射率计算土壤全盐含量的方法。

为达到上述发明目的，本发明所采用的技术方案为：

提供一种利用土壤高光谱反射率计算土壤全盐含量的方法，其包括以下步骤：

s1：在待测目标含盐土壤上选取若干测试点，在每个测试点上选取若干数据采集点，并采集每个数据采集点上的高光谱反射率；

s2：将每一个测试点上获得的若干高光谱反射率进行算术平均，绘制成一条高光谱反射曲线，即获得若干条高光谱反射曲线；

s3：采集每个测试点设定范围内的土样，并将其制备成标准土壤浸提液，测定标准土壤浸提液的电导率和全盐含量，建立电导率和全盐含量的关系式：

salt＝αec1:5+β，

其中，ec1:5为标准土壤浸提液的电导率(ds/m)，salt为土壤全盐量(％)，α和β均为关系系数；

s4：对每条高光谱反射曲线进行去包络线处理，获得若干条高光谱吸收特征曲线，读取每条高光谱吸收特征曲线中的两个吸收波谷的波长分布范围；

s5：采用一元多次回归方程，拟合第i条高光谱吸收特征曲线的两个吸收波谷的曲线表达式；

s6：通过两个吸收波谷的曲线表达式，计算第i条高光谱反射曲线对应的第i个测试点土壤的高光谱吸收指数；

s7：重复步骤s5-s6，获得其余测试点土壤的高光谱吸收指数；

s8：利用所有测试点土壤的高光谱吸收指数建立土壤全盐量和高光谱吸收指数的回归关系式，计算待测目标含盐土壤的土壤全盐量。

进一步地，步骤s1包括：

s11：在待测目标含盐土壤上选取若干个测试点，并以每个测试点为中心建立若干个1m²的正方形测试区域，并选取测试区域内的四个角点和中心点作为五个采集点；

s12：利用便携式地物光谱仪对每一个采集点进行三次土壤高光谱反射率测量，得到15组高光谱反射率的测量数据；

进一步地，步骤s3包括：

s31：以每个测试点为中心的1m²正方形区域内，选择五个点个点进行土壤的采集，并将采集的五份土壤混合成一个采集土样；

s32：将采集土样置于干燥通风处自然风干，并剔除采集土样内的杂物，然后搅匀样品，研磨，并通过2mm孔径筛筛除大颗粒物；

s33：将经过步骤s32处理后采集土样按照水土比例5：1进行混合，得到标准土壤浸提液；

s34：采用电导仪和烘干法分别测定标准土壤浸提液的电导率和全盐含量，建立电导率和全盐含量的关系式。

进一步地，步骤s5中两个波谷曲线的表达式fi,1(λ)和fi,2(λ)分别为：

其中，fi,1(λ)和fi,2(λ)分别为第i个土壤样品的高光谱吸收特征曲线的第一波谷曲线和第二波谷曲线的回归方程，和b^i,1均为第i个土壤样品的第一个波谷曲线的回归系数，和b^i,2为第i个土壤样品的第二个波谷曲线的回归系数；λi,1和λi,2为第一个波谷曲线和第二个波谷曲线的波长，m和k分别为第一波谷曲线和第二波谷曲线的回归方程的阶数。

进一步地，步骤s6中的高光谱吸收指数包括吸收谷深h：

其中，hi,1和hi,2分别为第一波谷曲线和第二个波谷曲线的谷深，ρi,1,max和ρi,1,min分别为第一波谷的反射率最大值和最小值，ρi,2,max和ρi,2,min分别为第二波谷的反射率最大值和最小值，λi,1,min和λi,1,max为第一波谷的波长最小值和最大值，λi,2,min和λi,2,max为第二波谷的波长最小值和最大值；

吸收度s：

其中，si,1和si,2分别为第一波谷和第二波谷的吸收度；

吸收幅宽w：

wi,1＝λi,1,w2-λi,1,w1，wi,2＝λi,2,w2-λi,2,w1，其中，wi,1和wi,2分别为第一波谷和第二波谷的幅宽。λi,1,w1和λi,1,w2为：以为纵坐标，绘制与x轴平行的直线，该直线与第一个吸收谷曲线的两个交点的横坐标。λi,2,w1，λi,2,w2为：以为纵坐标，绘制与x轴平行的直线，该直线与第二个吸收谷曲线的两个交点的横坐标。

进一步地，步骤s8包括：

s81：建立土壤电导率与高光谱吸收指数的回归方程：ec1：5＝p1×hi,1+p2×hi,2+p3×si,1+p4×si,2+p5×wi,1+p6×wi,2+p0，其中，p0～p6为回归系数；

s82：利用步骤s3中的电导率和全盐含量的关系式，建立土壤全盐量和高光谱吸收指数的回归关系式：

salt＝αec1:5+β＝α×p1×hi,1+α×p2×hi,2+α×p3×si,1+α×p4×si,2+α×p5×wi,1+α×p6×wi,2+α×p0+β。

进一步地，测试点不少于30个。

本发明的有益效果为：本方案利用含盐土壤的高光谱信息推算土壤中全盐含量，可以用于区域土壤盐分的快速监测和诊断，用户只需要通过监测有限个样品土壤的全盐量、电导率和相应的土壤高光谱反射率信息，即可利用本方法建立该区域同类土壤的全盐量与高光谱吸收指数的关系，并且基于土壤全盐量和高光谱吸收指数的回归关系式，用户可通过便携式地物高光谱仪快速地获取该区域同类土壤的高光谱吸收指数，即可计算得到土壤全盐量，而不需要用户再逐一费时费力地钻孔取样并通过实验室方法测量其土壤全盐量了。

附图说明

图1为土壤高光谱反射率计算土壤全盐含量的方法的流程图。

图2为若干高光谱反射曲线图。

图3为连续统去除处理后的土壤样品高光谱吸收特征曲线

图4为两个吸收谷的谷深示意图。

图5为两个吸收谷的吸收度示意图。

图6为两个吸收谷的吸收幅宽示意图。

图7为第一个土壤样品高光谱吸收特征曲线图。

图8为第一个土壤样品的两个吸收谷的吸收幅宽图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图1所示，利用土壤高光谱反射率计算土壤全盐含量的方法包括以下步骤：

s1：在待测目标含盐土壤上选取至少30个测试点，在每个测试点上选取5数据采集点，并采集每个数据采集点上的高光谱反射率。

步骤s1包括以下步骤：

s11：在待测目标含盐土壤上选取至少30个测试点，并以每个测试点为中心建立30个1m²的正方形测试区域，并选取测试区域内的四个角点和中心点作为五个采集点；

s12：利用便携式地物光谱仪对每一个采集点进行三次土壤高光谱反射率测量，则每个测试点可得到15组高光谱反射率的测量数据；

本实施例数据的采集选择在天气晴朗的上午10:00至下午3：00之间进行，采集人员着深色服装携带3便携式地物光谱仪，在目标检测区域选择n＝30个测试点进行土壤高光谱反射率数据采集，每个测试点按照步骤s11的方式选择5个采集点，每个采集点位采集3次，则每个测试点可获得15组高光谱反射率数据。

s2：将每一个测试点上获得的15组高光谱反射率数据进行算术平均，得到该测试点的平均高光谱反射曲线，30个测试点即可获得30条高光谱反射曲线。

s3：采集每个测试点设定范围内的土样，利用采集土样制作标准土壤浸提液，并测定标准土壤浸提液的电导率ec1：5(ds.m-1)和全盐含量salt(％)，建立电导率和全盐含量的关系式为：

salt＝0.3473ec1:5-0.0185，

其中，ec1:5为标准土壤浸提液的电导率(ds/m)，salt为土壤全盐量(％)，α＝0.3473，β＝-0.0185。

步骤s3具体的步骤为：

s31：以每个测试点为中心建立1m²正方形测试区域，选择测试区域内的四个角点和中心点作为五个采集点进行土壤的采集，并将采集的五份土壤混合成一个采集土样；

s32：将采集土样置于干燥通风处自然风干，并剔除采集土样内的杂物，然后搅匀样品，研磨，并通过2mm孔径筛筛除大颗粒物；

s33：将经过步骤s32处理后采集土样按照水土比例5：1进行混合，得到标准土壤浸提液；

s34：采用电导仪和烘干法分别测定标准土壤浸提液的电导率和全盐含量，建立电导率和全盐含量的关系式。

s4：对30条高光谱反射曲线进行去包络线处理，获得30条高光谱吸收特征曲线，如图2所示。读取每条高光谱吸收特征曲线中的两个吸收波谷的波长分布范围；本实施例选第1条高光谱吸收特征曲线，如图7所示，读取第1条曲线上两个吸收波谷的波长分布范围(λi,1,min,λi,1,max)和(λi,2,min,λi,2,max)，i＝1；本实施例中第1条曲线对应第一个土壤样品的电导率和全盐量数据。

在遥感软件envi环境下，利用连续统去除功能对每个土壤样品的高光谱反射曲线进行去包络线处理(菜单spectral–mappingmethods–continuumremoval)，获得每个测试点的高光谱吸收特征曲线，并逐个找出每个测试点的高光谱吸收特征曲线的两个吸收谷的波长具体分布范围(λi,1,min,λi,1,max)和(λi,2,min,λi,2,max)。

本实施例以第一个样品土壤为例，该土壤样品的高光谱吸收特征曲线的两个吸收谷的波长分布范围分别为(1323，1653)和(1849，2142)。

s5：利用一元多次回归方程拟合两个吸收波谷的波谷曲线表达式fi,1(λ)和fi,2(λ)：

其中，如图3所示，fi,1(λ)和fi,2(λ)分别为第i个土壤样品的高光谱吸收特征曲线的第一个和第二个吸收谷曲线的回归方程，和b^i,1为第i个土壤样品的第一个波谷曲线的回归系数，和b^i,2为第i个土壤样品的第二个波谷曲线的回归系数。λi,1和λi,2为第一个波谷曲线和第二个波谷曲线的波长，m和k分别为第一波谷曲线和第二波谷曲线的回归方程的阶数。

本实施例中表达式fi,1(λ)和fi,2(λ)中的具体为：

如图7所示，表达式中各变量在本实施例的具体取值如表1所示，

表1

s6：通过两个波谷曲线的表达式fi,1(λ)和fi,2(λ)计算第1条高光谱反射曲线对应的第1个测试点的高光谱吸收指数，高光谱吸收指数包括吸收谷深h、吸收度s和吸收幅宽w；

如图4所示，吸收谷深h为：

其中，hi,1和hi,1分别为第一波谷曲线和第二个波谷曲线的谷深，ρi,1,max和ρi,1,min分别为第一波谷的反射率最大值和最小值，ρi,2,max和ρi,2,min分别为第二波谷的反射率最大值和最小值，λi,1,min和λi,1,max为第一波谷的波长最小值和最大值，λi,2,min和λi,2,max为第一波谷的波长最小值和最大值。

本实施例中两个谷深hi,1和hi,1具体为：

如图5所示，吸收度s为：

其中，si,1和si,2分别为第一波谷和第二波谷的吸收度。

本实施例中吸收度si,1和si,2具体为：

如图6所示，吸收幅宽w为：

wi,1＝λi,1,w2-λi,1,w1，wi,2＝λi,2,w2-λi,2,w1，其中，wi,1和wi,2分别为第一波谷和第二波谷的幅宽。λi,1,w1和λi,1,w2为：以为纵坐标，绘制与x轴平行的直线，该直线与第一个吸收谷曲线的两个交点的横坐标。λi,2,w1，λi,2,w2为：以为纵坐标，绘制与x轴平行的直线，该直线与第二个吸收谷曲线的两个交点的横坐标。

如图8所示，本实施例中在第1个测试点的高光谱吸收特征曲线上，针对第1个吸收谷，以即(0.916175+0.083825/2＝0.9580875)为纵坐标，绘制与x轴平行的直线，该直线与第一个吸收谷曲线的两个交点的横坐标(λi,1,w1，λi,1,w2)分别为(1401，1535)，针对第2个吸收谷，以即(0.728538+0.271462/2＝0.864269)为纵坐标，绘制与x轴平行的直线，该直线与第二个吸收谷曲线的两个交点的横坐标(λi,2,w1，λi,2,w2)分别为(1899，2001)。

第1个土壤样品的两个吸收谷的幅宽分别为：w1＝λ1,w2-λi,1,w1＝1535-1401＝134，

w2＝λ2,w2-λi,2,w1＝2001-1899＝102；

根据上述关系式可以得出，第1个土壤样品的6个吸收指数(hi,1，hi,2，si,1，si,2，wi,1，wi,2)分别为(0.083825，0.271462，9.36051，28.41955，134，102)。

如图7和图8所示，表达式中各变量在本实施例的具体取值如表2所示，

表2

s7：重复步骤s5-s6，获得其余29个测试点土壤的6个高光谱吸收指数；

s8：利用所有30个测试点高光谱吸收指数建立土壤全盐量和高光谱吸收指数的回归关系式。

土壤电导率ec1：5与高光谱吸收指数的回归方程为：ec1：5＝p1×hi,1+p2×hi,2+p3×si,1+p4×si,2+p5×wi,1+p6×wi,2+p0，其中，p0～p6为回归系数。

以土壤电导率为因变量，以每个土壤样品的6个吸收指数为自变量，通过统计分析软件dps自带的多元稳健回归估计功能建立土壤电导率和吸收指数的回归关系。具体操作为：dps软件菜单栏—多元分析—回归分析—稳健回归(m估计)，本实施例回归方程如下：

ec1：5＝16.4015×hi,1+25.2694×hi,2-0.3036×si,1-0.1769×si,2

-0.0197×wi,1+0.0265×wi,2+0.0816

即土壤电导率ec1：5与高光谱吸收指数的回归方程中，p0＝0.0816，p1＝16.4015，p2＝25.2694，p3＝-0.3036，p4＝-0.1769，p5＝-0.0197，p6＝0.0265。

土壤全盐量的含量和高光谱吸收指数的回归关系式为：

salt＝αec1:5+β＝α×p1×hi,1+α×p2×hi,2+α×p3×si,1+α×p4×si,2+α×p5×wi,1+α×p6×wi,2+α×p0+β。

本实施例回归方程如下：salt＝αec1:5+β＝0.3473×(16.4015×hi,1+25.2694×hi,2-0.3036×si,1-0.1769×si,2-0.0197×wi,1+0.0265×wi,2+0.0816)-0.0185＝5.6962×hi,1+8.7761×hi,2-0.1054×si,1-0.0614×si,2-0.0069×wi,1+0.0092×wi,2-0.0101。利用该公式计算待测目标含盐土壤的土壤全盐量。

本方案利用含盐土壤的高光谱信息推算土壤中全盐含量，可以用于区域土壤盐分的快速监测和诊断，用户只需要通过监测有限个样品土壤的全盐量、电导率和相应的土壤高光谱反射率信息，即可利用本方法建立该区域同类土壤的全盐量与高光谱吸收指数的推算关系，并且基于土壤全盐量和高光谱吸收指数的回归关系式，用户可通过便携式地物高光谱仪快速地获取该区域同类土壤的高光谱信息，即可简单推求得到土壤全盐量，而不需要用户再逐一费时费力地钻孔取样并通过实验室方法测量其土壤全盐量了。