基于地面高光谱的土壤有机质含量定量反演方法与流程

本发明涉及土壤有机质含量检测技术，具体是基于地面高光谱的土壤有机质含量定量反演方法。

背景技术：

土壤有机质含量是土壤肥力评价的重要指标，土壤有机质含量的检测对于土壤有机碳储量的估算具有重要意义。现有的土壤有机质检测方法主要包括土地统计学、土壤-景观关系理论及遥感反演三种方法，其中，遥感反演方法利用遥感数据建立多因素模型，其理论依据是在可见光和近红外波段中土壤有机质含量具有特定的光谱特性，即土壤有机质含量越高，反射率越低。遥感反演方法的显著优势在于其能够开展区域尺度土壤有机质时空变化研究，因此，其成为土壤有机质检测的发展趋势。然而，目前遥感反演方法只能定性地分析土壤有机质的成份，无法实现对土壤有机质的精确定量分析，给土地质量调查带来诸多限制。

技术实现要素：

本发明的目的在于克服现有遥感反演方法无法实现对土壤有机质的精确定量分析的问题，提供了一种基于地面高光谱的土壤有机质含量定量反演方法，其能快速的分析土壤有机质含量，低成本，且能提升分析精度。

本发明的目的主要通过以下技术方案实现：基于地面高光谱的土壤有机质含量定量反演方法，包括以下步骤：

步骤1、采集土壤样品；

步骤2、采用地物光谱仪采集土壤样品的高光谱数据，并采用平滑曲线方法进行去噪预处理；

步骤3、对土壤样品的有机质含量进行化学测试分析，以获得对应土壤样本有机质含量数据；

步骤4、将经过预处理的土壤高光谱曲线数据与化学测试分析后的土壤有机质化学分析数据进行相关性分析，选择相关性系数较大的前若干个反射率光谱作为特征光谱；

步骤5、将高光谱数据中的特征光谱采用动态参数多元回归法得到土壤有机质拟合模型，用以拟合整个区域的土壤有机质含量。本发明通过预处理和化学测试分析供后续与光谱相关性分析使用。本发明在具体实施时，步骤2和步骤3顺序可互换，也可同步进行。

进一步的，所述步骤1具体包括如下步骤：

步骤1.1、规划野外土壤样品数据采样点，其中，采样点的数量为多个；

步骤1.2、采集预定采样点的近地表土壤样品，混合得到混合土壤样本；其中，近地表土壤指地表以下0～10cm处的土壤。

进一步的，所述步骤1还包括将采集的土壤样品进行干燥、过筛及拌匀处理，其中，过筛时采用-20目尼龙筛过筛。本发明通过干燥、过筛及拌匀处理，能提高化学分析的精度和高光谱采样的准确性，能消除土壤湿度、颗粒度的影响。

进一步的，所述步骤2的土壤样品光谱采样使用多点多次采集方法；所述平滑曲线方法为：以平滑点为中心，前后开窗口取若干相邻数据进行平均或拟合，进而求得平滑点的最优估计值，其公式如下：

式中，和xi分别为平滑前后光谱数据向量的一个元素，wj是窗口平滑时长为2r+1的权重因子，r为移动窗口的长度。土壤样品光谱采样使用多点多次采集方法，能消除光谱误差和噪声。进行移动窗口平滑时，其平滑窗口的大小影响平滑效果。

进一步的，所述步骤2的土壤样品光谱采样时，将加工完毕的土壤样品置于实验室暗室，将土壤样品均匀平铺于0.5m╳0.5m的黑布之上，使用高光谱光谱仪以横向间隔0.1m采集土壤高光谱曲线数据，每个土壤样品采集20～25条高光谱曲线数据；其中，实验室暗室墙面反光率小于60％、地面反光率小于20％。

进一步的，所述步骤4中的土壤高光谱曲线数据在与化学测试分析后的土壤有机质化学分析数据进行相关性分析前，对原始反射率光谱依次进行一阶微分、倒数一阶微分、对数一阶微分、倒对数一阶微分和平方根一阶微分变换变换；所述步骤4中的相关性系数计算公式如下：

式中，i为土壤样本序号，取值为1到n；n为土壤样本总数；j为波段序号；rj为土壤光谱在j波段的反射率与有机质含量的相关系数，取值介于-1到1之间；xij为第i个土壤样本在j波段处的反射率；为n个土壤样本在该波段处的光谱平均值；yi为第i个土壤样本有机质的含量；表示所有土壤样品有机质含量的平均值。

进一步的，所述动态参数多元回归法采用向前选择和向后剔除的组合，根据全部自变量对因变量的贡献值大小来动态的选择优先进入方程的光谱参数自变量，即通过偏f检验来筛选光谱参数变量，f统计量包括了fin和fout，使用的公式为：

其中，x1，...，xi，...，xk为k个光谱参数自变量；

在引入新的第k+1个光谱参数自变量时，剔除第i个光谱参数自变量xi时，

式中，sse为因变量的预测误差平方和，n为样本个数；

在动态选择光谱参数xi+1进入fin时，按α＝0.05的水平进行显著性检验，从fout中剔除光谱参数自变量时，按p＝0.1的检验概率值进行，偏f检验拒绝域的临界值为fin>fout；

xi+1引入，进行偏f检验，通过则进行下一步新光谱参数自变量xi+1的引进；不通过则把xi从回归模型中剔除；

若新自变量xi+1引入时，原被剔除的xi-1自变量又满足α＝0.05，则通过检验，重新选入；

所有模型外部自变量x都不能通过偏f值检验时结束。如此，本发明使用动态参数多元线性回归法，按规则对输入的自变量进行动态筛选，找出最优光谱参数组合。限定偏f检验拒绝域的临界值为fin>fout，以避免变量x的进出循环。

进一步的，所述动态参数多元回归法得到的研究区有机质含量反演模型为：

y＝β0+β1x1++β2x2+...+βqxq+ε(3)

其中，y是拟合的有机质含量数值，x1、x2、...、xq分别表示不同波长的反射率，q为选择的特征光谱数量，β0,β1,β2,...，βq称为回归参数，ε为随机误差项。本发明应用时，在研究区任意一点采集地面高光谱数据，输入该点的高光谱敏感波段的反射率数据，就可以反演出同一点的有机质含量。

综上所述，本发明与现有技术相比具有以下有益效果：(1)本发明通过对采集的土壤样品采用地物光谱仪采集土壤样品的高光谱数据，并采用平滑曲线方法进行去噪预处理，并对土壤样品的有机质含量进行化学测试分析，以获得对应土壤样本有机质含量数据，再进行相关性分析，然后将特征光谱采用动态参数多元回归法得到土壤有机质拟合模型，用以拟合整个区域的土壤有机质含量，能快速的分析土壤有机质含量，低成本，且能提升分析精度。

(2)本发明的动态参数多元线性回归法不仅能通过变量选择来减小共线性的影响，防止过度拟合，同时还能考虑到光谱参数自变量之间存在的相互传递的相关关系，从而提高模型的精度。

附图说明

此处所说明的附图用来提供对本发明实施例的进一步理解，构成本申请的一部分，并不构成对本发明实施例的限定。在附图中：

图1为本发明一个具体实施例的流程图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施例和附图，对本发明作进一步的详细说明，本发明的示意性实施方式及其说明仅用于解释本发明，并不作为对本发明的限定。

实施例：

如图1所示，基于地面高光谱的土壤有机质含量定量反演方法，包括以下步骤：步骤1、采集土壤样品；步骤2、采用地物光谱仪采集土壤样品的高光谱数据，并采用平滑曲线方法进行去噪预处理；步骤3、对土壤样品的有机质含量进行化学测试分析，以获得对应土壤样本有机质含量数据；步骤4、将经过预处理的土壤高光谱曲线数据与化学测试分析后的土壤有机质化学分析数据进行相关性分析，选择相关性系数较大的前若干个反射率光谱作为特征光谱；步骤5、将高光谱数据中的特征光谱采用动态参数多元回归法得到土壤有机质拟合模型，用以拟合整个区域的土壤有机质含量。

本实施例的步骤1具体包括如下步骤：步骤1.1、规划野外土壤样品数据采样点，其中，采样点的数量为多个；步骤1.2、采集预定采样点的近地表土壤样品，混合得到混合土壤样本；其中，近地表土壤指地表以下0～10cm处的土壤。本实施例在采集土壤样品时，先规划野外土壤样品数据采样点，以采集该地具有代表性的土壤样本。本实施例优选采用区域网格法对一定范围内的土地进行样本采样，同时考虑实地的地形地貌、土地类型、土地利用现状等因素，兼顾到空间分布的均匀性，以确保采样点具备该地土壤类型的代表性。本实施例的野外土壤样品采集方法为：采集预定采样点的近地表土壤样品，若干份土壤样品混合形成一份表征该采样点土地特性的混合土壤样本。采样过程中应去除杂物，所述的杂质包括杂草、草根、砾石、砖块、肥料团块等非土壤成分。优选地，为增加土地样品的代表性，采样时应在采样点周围25m～30m范围内不同的田块、不同的农作物类型，多点采集3-7个子样组合为一个样品。将采集的各分样点土壤充分混合后，留取1.0kg～2kg装入样品袋，样品袋选用可密封的塑料袋。

本实施例的步骤1还包括将采集的土壤样品进行干燥、过筛及拌匀处理，其中，过筛时采用-20目尼龙筛过筛。本实施例在对采集的土壤样品进行干燥时，在日光下干燥，在样品干燥过程中揉搓样品，以免胶结，干燥后的样品在过筛前用木槌轻轻敲打，以便使土壤样品恢复至自然粒级状态。干燥后的土壤样品过筛时，采用尼龙筛过筛，将筛下粒级样品拌匀，取450g～600g装瓶。本实施例的土壤样品加工场地需干净、通风、无污染，每加工完一个土壤样品后需对加工工具进行全面清扫，防止样品间的相互污染。样品瓶使用聚乙烯塑料瓶，在装样前对样品瓶进行清洗，烘干，保持清洁。为避免土壤样品之间相互污染，样品的采集、运输、保管、晾晒、加工应由不同人员分别进行，并且不能在同时或同一场地进行，土壤采样用的新布袋及使用过的样品布袋都要经过洗涤才能继续使用。如采集的土壤样湿度较大，需套上塑料袋使湿样品分开。土壤样品加工完毕后，填写送样单，与标签对应，及时送予实验室分析测试。

本实施例的步骤2的土壤样品光谱采样使用多点多次采集方法；所述平滑曲线方法为：以平滑点为中心，前后开窗口取若干相邻数据进行平均或拟合，进而求得平滑点的最优估计值，其公式如下：

式中，和xi分别为平滑前后光谱数据向量的一个元素，wj是窗口平滑时长为2r+1的权重因子，r为移动窗口的长度。

本实施例的步骤2的土壤样品采用多点多次采集方法，其具体实施步骤为：将加工完毕的土壤样品置于实验室暗室(内部光源可关闭、墙面反光率小于60％、地面反光率小于20％)，将土壤样品均匀平铺于0.5m╳0.5m的黑布之上，使用高光谱光谱仪以横向间隔0.1m采集土壤高光谱曲线数据，每个土壤样品采集20～25条高光谱曲线数据。高光谱数据采集使用美国asd公司生产的fieldspecpro便携式地物光谱仪，asd便携式光谱仪获得的光谱测量数据，是由光源、样品、传感器、测量环境组成的复杂数据。asd野外光谱分析仪fieldspecpro是一种测量可见光到近红外波段地物波谱的有效工具。它能快速扫描地物，光纤探头能在毫秒内得到地物单一光谱，高光谱地面采样就非常便捷。高光谱采样的原理是采集样品对象的反射率，墙面反射率和地面反射率选取不合适，会形成光谱噪声，影响高光谱采样效果，本实施例应用时限定墙面反光率小于60％、地面反光率小于20％，以避免强光影响采样效果。

本实施例的化学分析方法采用现有的氧化还原法(vol)方法实现，本实施例进行化学分析时，取70g～100g样品(潮湿的土壤样品在加工前于45℃以下的烘箱中烘干后称重使用)用无污染的行星球磨机粉碎，从加工后的试样中分取30克试样装玻璃瓶于45℃烘箱中，烘两小时后送原子荧光组做有机质的取样做化学测试分析。

本实施例的步骤4中的土壤高光谱曲线数据在与化学测试分析后的土壤有机质化学分析数据进行相关性分析前，对原始反射率光谱依次进行一阶微分、倒数一阶微分、对数一阶微分、倒对数一阶微分和平方根一阶微分变换变换；所述步骤4中的相关性系数计算公式如下：

式中，i为土壤样本序号，取值为1到n；n为土壤样本总数；j为波段序号；rj为土壤光谱在j波段的反射率与有机质含量的相关系数，取值介于-1到1之间；xij为第i个土壤样本在j波段处的反射率；为n个土壤样本在该波段处的光谱平均值；yi为第i个土壤样本有机质的含量；表示所有土壤样品有机质含量的平均值。本实施例优选地选择相关性系数最大的前两个反射率光谱作为特征光谱。

发明人对不同光谱分辨率、不同变换形式的土壤反射率光谱曲线与有机质相关性统计发现：土壤原始反射率光谱与有机质呈负相关，且在不同光谱分辨率下，土壤原始反射率光谱与有机质的相关系数的绝对值均在0.4-0.6之间。对原始反射率光谱进行二阶微分、倒数、倒数二阶微分、对数、对数二阶微分、倒对数、倒对数二阶微分、平方根和平方根二阶微分变换后，土壤原始反射率光谱与有机质的相关系数只有微弱的增大；而对原始反射率光谱进行一阶微分、倒数一阶微分、对数一阶微分、倒对数一阶微分和平方根一阶微分变换后，土壤原始反射率光谱与有机质的相关系数的绝对值能提高到0.68-0.8左右。相关系数较高的波段主要位于1200～1250nm、1500～1800nm，而土壤原始反射率光谱曲线在相应的波段亦有许多微弱或小的吸收峰。根据平方根一阶微分在不同光谱分辨率下与有机质的相关性可知，适当改变光谱分辨率可以增大变换形式与有机质之间的相关性。针对有机质，适合进行反演的光谱分辨率以及变换形式可为：一阶微分、倒数一阶微分、对数一阶微分、倒对数一阶微分和平方根一阶微分。

本实施例采用的多元回归法是指在相关变量中将一个变量视为因变量，其他一个或多个变量视为自变量，建立多个变量之间线性或非线性数学模型数量关系式并利用样本数据进行分析的统计分析方法。多元回归模型的数学形式为设因变量为y，影响因变量的k个自变量分别为x1,x2,...，xk,假设每一个自变量对因变量y的影响都是线性的，也就是说，在其他自变量不变的情况下，y的均值随着自变量xk的变化均匀变化，这时我们把y＝β0+β1x1++β2x2+...+βkxk+ε称为总体回归模型，把β0,β1,β2,...，βk称为回归参数，ε为随机误差项。回归分析的基本任务是：任务1、利用样本数据对模型参数作出估计；任务2、对模型参数进行假设检验；任务3、应用回归模型对因变量(被解释变量)作出预测。

为了保证多元回归分析的参数估计、统计检验以及置信区间估计的有效性，与一元线性回归分析类似，我们需要对总体回归模型及数据作一些基本假定。

假定1：随机误差项ε的概率分布具有零均值，即e(ε)＝0。

假定2：随机误差项ε的概率分布对于不同的自变量表现值而言，具有同方差。即ε的方差不随着第i个土壤样本在j波段处的反射率xij的变化而变化，d(ε)＝σ²。

假定3：随机误差项ε不存在自相关，即cov(εi,εj)＝0。

假定4：εi与任一解释变量xi不相关，可以表示为cov(εi,xi)＝0。

假定5：解释变量x之间不存在完全共线性。以上假定1～4与一元回归分析的假定是相同的。假定5是针对解释变量而言，在一元回归分析中，由于只有一个解释变量，因此这一点是不需要的。在模型和数据满足上述假定时，对上述总体回归模型两边取期望，可得到：

e(y|x1，x2，...，xk)＝β0+β1x1+β2x2+...+βkxk

上式称为总体回归方程(populationregressionequation，pre)或总体回归函数(populationregressionfunction，prf)，e(y|x1，x2，...，xk)表示在给定自变量xi的条件下观察值y的条件均值。在实际问题中，总体参数β0,β1,β2,...，βk往往是未知的，我们需要根据样本观察值给出总体参数的相应的估计值此时，

该式称为样本回归方程(sampleregressionequation，sre)或样本回归函数(sampleregressionfunction，srf)，也就是e(y|x1，x2，...，xk)的点估计值。

对于多元回归方程，在模型和数据满足前文所述的基本假定的前提下，参数估计可以通过最小二乘估计来得到，同样假设即

根据高等数学知识，q分别对

对求偏导数，令其等于0，得到

求解上式(中的方程组，即可得到参数的估计值

本实施例的动态参数多元回归法采用向前选择和向后剔除的组合，根据全部自变量对因变量的贡献值大小来动态的选择优先进入方程的光谱参数自变量，即通过偏f检验来筛选光谱参数变量，f统计量包括了fin和fout，使用的公式为：

其中，x1，...，xi，...，xk为k个光谱参数自变量；

在引入新的第k+1个光谱参数自变量时，剔除第i个光谱参数自变量xi时，

式中，sse为因变量的预测误差平方和，n为样本个数；

在动态选择光谱参数xi+1进入fin时，按α＝0.05的水平进行显著性检验，从fout中剔除光谱参数自变量时，按p＝0.1的检验概率值进行，偏f检验拒绝域的临界值为fin>fout；其中，α、p为概率统计中的显著性检验参数。

xi+1引入，进行偏f检验，通过则进行下一步新光谱参数自变量xi+1的引进；不通过则把xi从回归模型中剔除；

若新自变量xi+1引入时，原被剔除的xi-1自变量又满足α＝0.05，则通过检验，重新选入；

所有模型外部自变量x都不能通过偏f值检验时结束。

本实施例的动态参数多元回归法得到的研究区有机质含量反演模型为：

y＝β0+β1x1++β2x2+...+βqxq+ε(5)

其中，y是拟合的有机质含量数值，x1、x2、...、xq分别表示不同波长的反射率，q为选择的特征光谱数量，β0,β1,β2,...，βq称为回归参数，ε为随机误差项。

本实施例通过步骤4选择相关性系数较大的前8个反射率光谱作为特征光谱，式(3)中43.588、13.4449、3.4161、0.6043、4.2284、0.7607、0.6529、4.9248及0.2831总共9个数值由步骤5中构建反演模型得到的，ε为0。面积得到的研究区有机质含量反演模型为：

y＝43.588+13.4449x1+3.4161x2+0.6043x3-4.2284x4+0.7607x5

+0.6529x6+4.9248x7-0.2831x8。

运用野外实地采集的数据验证本实施例反演模型的精度，以验证模型对区域内有机质含量反演的有效性，并对该反演模型进行评价，其具体实施步骤如下：

模型的验证优选主要通过交叉验证均方根误差rmsecv指标对比，rmsecv可以近似看做模型泛化误差的无偏估计，模型优选的目标通常是使rmsecv最小化的过程，即选择泛化误差最小的模型作为优选的结果。

ym，yp分别为测量值和预测值，ncv交叉验证的样本个数。

模型的建模精度，即模型对于建模样本的模拟效果，主要通过建模均方根误差rmsec及建模决定系数r2评价；

nc为建模样本个数。

模型的预测精度，即模型对于预测样本的模拟效果，主要通过预测均方根误差rmsep、预测决定系数r2评价。rmsec、rmsep越小，r2和rp2越大说明建模精度和预测精度越好。

为实测样本的平均值。rpd为样本值标准差与估算均方根误差比，主要用来评价模型的预测能力，该值越大模型的预测能力越好，rpd>2时可以认为该模型有较好的预测能力，当1.4<rpd<2时则认为可对样品进行粗略估计，当rpd<1.4则说明模型预测能力较差，不能应用到样品预测中。

rpd＝stdev(ym)/rmsep(9)。

常规的土地质量调查需要地球化学分析，传统的土壤质量检测方法需投入大量的人力物力对研究区进行地球化学采样。以1:5万地球化学采样为例：采样密度为4-8个点/km2，每个土地采样点需要耗时一小时以上，每个土壤样点全分析费用为320元/个，进行野外实地监测通常需要长时间和大范围的密集监测。传统土壤质量检测方法劳动强度大、效率低下、周期漫长、成本耗高昂，这导致基础数据的获得非常困难。传统的土壤质量检测方法只能说哪些波段与有机质含量高低有关(定性)，或是给出含量的相对多少(不准确)。本实施例应用时，只需要对地面样点采集高光谱数据，就可以反演该点的有机质含量，可以比较准确计算出含量多少(定量)。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。