土壤有机质测定方法及测定系统

1.本发明涉及一种土壤中有机质的测定技术领域，具体涉及土壤有机质测定方法及测定系统。


背景技术：

2.土壤是农业生产中重要的自然资源，与人类生产活动有着极为密切的关系。农作物的生长离不开土壤，土壤为农作物提供各种营养成分，但土壤中的养分分布不均匀且处于动态之中，另外，我国目前的施肥技术大多还停留在经验施肥阶段，生产中存在着较多肥料施用不合理及肥料利用率低等现象，直接影响了作物的生长发育，同时，对环境也造成了一定的污染。因此，采用先进的农业生产技术准确及时的获取土壤中的养分含量信息，可为肥料合理利用和作物精确栽培管理提供依据和参考。
3.全氮、有机质作为土壤主要养分含量，对作物生长具有非常重要的作用。具体来讲，全氮是衡量土壤氮素供应状况和评价作物的氮素营养状况的重要指标，其丰缺直接影响作物的产量；有机质的含量大小是评估土壤肥力水平的一项重要指标，其对全球碳平衡过程有着至关重要的作用。利用传统的实验室化验分析方法提取土壤全氮和有机质含量的过程操作繁琐，既费时又费力，不能准确及时应用到实际生产中，并且存在一定程度的污染和浪费。高光谱遥感技术作为重要的监测手段之一，可以快速、无损地收集土壤信息，获取的数据量大，同时包含多方面的信息，能够对多种土壤养分及影响因素同时进行分析研究，构建不同的养分估测模型，对土壤信息快速准确获取具有重要意义。


技术实现要素：

4.为了解决上述技术问题，本发明提出了土壤有机质测定方法，包括如下步骤：
5.步骤一，土壤样品采集及土壤样品高光谱测定，获得土壤样本的反射光谱数据；
6.步骤二，数据分析，剔除反射光谱数据中的异常值，对剔除所述异常值后的土壤样本的反射光谱进行变换，对建模需要的特征波段进行选择；
7.步骤三，构建特征波段的反射光谱数据与有机质含量之间的估测模型。
8.进一步地，采用垂线图示法剔除所述异常值，利用界限参数m，构建垂线图的垂线延伸的边缘：
9.m＝med h(pi,pj),pi≤m≤pj；
10.式中：h(pi,pj)＝{p1,p2,
…
,pn}为按升序排列的样本序列，med为求取中位数函数，m为样本序列的中位数。
11.进一步地，当pi≠pj时，
[0012][0013]
当pi＝pj＝m时，n1＜n2＜
…
＜nk表示与中位数相等的样本序号，表示与中位数相等的样本序号，
[0014][0015]
进一步地，所述步骤二中，基于最小二乘的savitzky-golay平滑滤波去除光谱噪声，设光谱曲线为t次多项式，其中t为自然数，其公式为：
[0016]
y＝a0+a1i+a2i2+a3i3+...+a
tit
；
[0017]
其中y表示光谱数据值，a0、a1、a2、a3、ai为多项式系数，i表示光谱波长，t为自然数。
[0018]
进一步地，采用一阶微分变换进行光谱变换，一阶微分变换的方程式为：
[0019][0020]
式中，λi是波长值，ρ
′
(λi)是λi处的一阶微分光谱数据，ρ(λ)是λi处的原始光谱数据，λ是范围内的所有波长值,
[0021]
λi＝(λ
i-1
+λ
i+1
)/2,δλ＝λ
i+1-λ
i-1
。
[0022]
进一步地，所述对建模需要的特征波段进行选择，包括：进行相关分析后，采用交叉校验均方根误差作为适应度函数进行优化目标；
[0023]
所述相关分析的计算公式如下:
[0024][0025]
式中，ri为土壤全氮含量与原始光谱波段的单相关系数，n为土壤样本的个数，i为波段序号，r
ni
表示第n个土壤样本第z个波段的原始光谱或其变换形式光谱值，为n个土壤样本在i波段原始光谱或其变换形式光谱值的平均值，tn为第n个土壤样品的全氮含量，为实测的n个土壤样品全氮含量的平均值。
[0026]
进一步地，所述步骤三中，估测模型的构建算法如下：
[0027]
将n个土壤样品m个组分的土壤有机质含量浓度矩阵y和n个土壤样品p个波长点处吸光度矩阵x，分解成特征向量形式:
[0028]
y＝uq+f；
[0029]
x＝tp+e；
[0030]
其中u和t分别为n行d列的有机质含量浓度特征因子矩阵和吸光度特征因子矩阵，q为d
×
m阶浓度载荷矩阵，p为d
×
p阶吸光度载荷矩阵，f和e分别为n
×
m，n
×
p阶浓度残差阵和吸光度残差阵；
[0031]
根据特征向量的相关性分解y和x，建立估测模型，吸光度向量为x，待测样本有机质浓度y为:
[0032]
y＝xt+bq；
[0033]
b＝(u-ed)/t；
[0034]
其中ed为随机误差阵，b为d维对角回归系数矩阵。
[0035]
本发明还提出了土壤有机质测定系统，其用于实现土壤有机质测定方法，包括：样品采集单元、高光谱数据测量单元、数据库、处理器和模型构建单元；
[0036]
所述样品采集单元，用于依据土壤类型和位置采集土壤样品，并采用gps记录采样位置的信息；
[0037]
所述高光谱数据测量单元，用于对采集的土壤样本进行反射光谱测量，获得土壤样本的反射光谱数据，并存储在所述数据库中；
[0038]
所述处理器，用于剔除反射光谱数据中的异常值，对剔除所述异常值后的土壤样本的反射光谱进行变换，对建模需要的特征波段进行选择；
[0039]
所述模型构建单元；用于构建特征波段的反射光谱数据与有机质含量之间的估测模型。
[0040]
本发明应用时，能快速的分析土壤有机质含量，低成本，且能提升分析精度，对于减少数据量，尤其是去除冗余信息起到了很好的作用，能够快速估算土壤中的有机质含量。
附图说明
[0041]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0042]
附图1是本发明的土壤有机质测定方法的流程示意图；
[0043]
附图2是本发明的土壤样品采集数据示意图；
[0044]
附图3是本发明的垂线图的组成示意图；
[0045]
附图4为剔除异常极值前后的数据分布对比情况。
具体实施方式
[0046]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0047]
参考图1为本发明的土壤有机质测定方法的流程图，其包括如下步骤：
[0048]
步骤1、样品采集及土壤高光谱测定。
[0049]
样品采集，将采集的新鲜土壤样本平面铺开，放置在通风处，让土壤样本自然风干。风干之后，用木棒压磨并除去砂砾等杂质物。将每个土壤样本平均分为多份，用于高光谱测试。
[0050]
在优选实施例中，采集土壤样品均为0～20cm的表层土壤，野外采集，从不同地点收集到16个土类1581个土壤样品。如图2所示，为土壤样品采集数据示意图。对采集的土壤样本进行风干、研磨并通过2mm孔筛。
[0051]
土壤高光谱测定，采用地物光谱仪进行土壤样本的光谱反射率数据测量。在优选实施例中，采用美国asd公司的fieldspecprofr型光谱仪对采集的1581个土壤样本进行土壤测试，其波长范围是350～2500nm，采样间隔为1.4nm(350～1000nm区间)和2nm(1000～
2500nm区间)，重采样间隔为1nm，输出波段数为2150。光谱测量在一个能控制光照条件的暗室内进行。光源是功率为50w的卤素灯，距土壤样品表面70cm，天顶角30
°
。待测土壤样本分别放置在直径10cm、深1.5cm的盛样器皿内，土样表面刮平。传感器探头位于土壤样本表面垂直上方15cm处。采用25
°
视场角探头，探头接受土壤光谱的区域为直径6.7cm的圆，远小于盛样器皿面积，保证探头接收的均为土壤的反射光谱。测试之前先进行白板校正。
[0052]
每个土样采集10条光谱曲线，将采集到10条光谱曲线进行算术平均计算后得到该土样的实际反射光谱数据。测试完毕后，导出土壤样本的原始光谱数据，对每一份光谱曲线数据导出至excel表格中保存，利用matlab求出多个样本曲线的平均值作为最后样本的光谱反射率数据进行保存。
[0053]
步骤2、数据分析，具体包括：
[0054]
s21、反射光谱数据中的异常值。由于异常值的存在会影响实验结果的精度，因此数据分析之前首先需要对异常值进行剔除。
[0055]
本实施例中采用垂线图示法。垂线图的组成如图3所示。垂线图的中心位置为中数，垂线图的中间块体长度表示垂线图的间距(d)，两端分别是上边界数(h
上
)和下边界数(h
下
)，垂线图两端的垂线向上所延伸到的位置为最大值，向下所延伸到的位置为最小值。垂线图示法定义的异常值是指样本数据中大于h
上
+1.5d和小于h
下
－1.5d的数值。其中h
上
+1.5d和h
下
－1.5d为异常值截止线，h
上
+3d和h
下
－3d为极端值截止线，介于异常值截止线与极端值截止线之间的异常值叫边缘异常值，位于极端值截距线以外的异常值叫异常极值。
[0056]
利用界限参数m，构建垂线图的垂线延伸的边缘。界限参数m的表达形式为：
[0057]
m＝med h(pi,pj),pi≤m≤pjꢀꢀꢀ
(1)；
[0058]
式中：h(pi,pj)＝{p1,p2,
…
,pn}为按升序排列的样本序列，med为求取中位数函数，m为样本序列的中位数。
[0059]
当pi≠pj时，
[0060][0061]
当pi＝pj＝m时，n1＜n2＜
…
＜nk表示与中位数相等的样本序号，即
[0062][0063][0064]
由式(1)-(3)可知，界限参数m的分布区间为[-1，1]。当数据序列向右偏时，界限参数m可定义为正值；数据序列向左偏时，界限参数m可定义为负值；当数据序列对称时，界限参数m为0。
[0065]
充分利用界限参数m的特性，将其引入垂线图的垂线延伸的边缘，垂线图的垂线延伸的边缘区间为[h
下-fm(m)
·
d，h
上
+fn(m)
·
d],d＝h
上
－h
下
；
[0066]
其中，fm(m)，fn(m)是界限参数m的线性函数。
[0067]
fm(m)＝1.5+am,fn(m)＝1.5+bm,其中a和b为线性系数。
[0068]
根据垂线图示法对样本数据的数值状态进行标注，区分出边缘异常值和异常极值，并根据数据分析的精确程度，进一步地选择是剔除边缘异常值还是异常极值。参考图4
为剔除异常极值前后的数据分布对比情况。
[0069]
s22、对剔除异常极值或剔除边缘异常值的土壤样本反射光谱进行变换。
[0070]
为了有效消除在光谱测量时各种非目标因素对土壤光谱的影响，并且在一定程度上削弱土壤类型等对土壤养分高光谱模型的影响，放大原始光谱数据中的细微信息。
[0071]
本发明将对土壤样本光谱曲线进行土壤光谱变换，具体地用基于最小二乘的savitzky-golay平滑滤波去除光谱噪声。
[0072]
基于最小二乘的savitzky-golay平滑滤波方法使用简化的最小二乘拟合卷积方法。假设光谱曲线为t次多项式，其中t为自然数，其公式为：
[0073]
y＝a0+a1i+a2i2+a3i3+...+a
tit
；
[0074]
其中y表示光谱数据值，a0、a1、a2、a3、ai为多项式系数，i表示光谱波长，t为自然数。
[0075]
基于最小二乘的savitzky-golay平滑滤波方法的推导过程中，最后简化的一般最小二乘卷积方程为：
[0076][0077]
其中，ci是窗口平滑时第i个光谱值系数，n是卷积中点值个数，y是原始光谱数据值，y是平滑后的光谱数据值，n＝2m+1为平滑窗口大小，对处于光谱两端的2m个(开始m个和最后m个)点，利用公式计算其平滑值。(i，j)为平滑窗口坐标值。
[0078]
采用一阶微分变换进行光谱变换。一阶微分变换可以提高光谱数据与有机质含量之间的相关系数，快速确定光谱曲线的拐点及最大、最小反射率的波长位置。
[0079]
一阶微分变换的方程式为：
[0080][0081]
式中，λi是波长值，ρ
′
(λi)是λi处的一阶微分光谱数据，ρ(λ)是λi处的原始光谱数据，λ是范围内的所有波长值,
[0082]
λi＝(λ
i-1
+λ
i+1
)/2,δλ＝λ
i+1-λ
i-1
。
[0083]
s23、对建模需要的特征波段进行选择。
[0084]
进行光谱变换后的光谱数据中，有效信息较弱，会存在大量的冗余信息。在具体的建模分析中，需要从光谱数据中挑选出部分特征波段，减少建模的工作量，提供工作效率。
[0085]
首先，进行相关分析，相关分析是分别对土壤样本的全氮和有机质含量与进行变换后的光谱进行各个波段的相关分析，计算每个波段与选定的土壤样本的有机质含量的相关系数。在变换后的光谱与波段的相关系数曲线中，选取相关系数明显的波谷和波峰作为敏感波段。相关系数的计算公式如下:
[0086][0087]
式中，ri为土壤有机质含量与光谱波段的单相关系数，n为土壤样本的个数，i为波
段序号。r
ni
表示第n个土壤样本第z个波段的光谱值，为n个土壤样本在i波段光谱值的平均值，tn为第n个土壤样品的有机质含量，为实测的n个土壤样品有机质含量的平均值。
[0088]
其次，优化目标算法。在进行土壤养分信息高光谱的特征波段选择时，采用pls交叉校验均方根误差rmsecv作为适应度函数，当最优土壤有机质含量与光谱波段的单相关系数达到给定的阈值或者优化次数达到预设的代数时，终止运算。
[0089]
适应度函数由偏最小二乘回归中的交叉性均方根误差rmsecv和决定系数rz决定。由于优化次数越多，波段之间选中的次数差异会增大，频率高的波长数目会明显减少且更稳定，对每次优化目标运行出来的结果进行对比分析，最后选择运行结果最好的一次作为入选的特征波段。在优选实施例中，优化目标算法需要运行10次来显示运行结果。
[0090]
步骤3、构建特征波段的反射光谱数据与有机质含量之间的估测模型。
[0091]
构建高光谱数据与有机质含量之间的估测模型，精简光谱自变量和去除冗余噪声，能够对土壤有机质含量信息进行有效筛选，避免了高光谱建模中的过拟合问题，最后建立的模型性最强的综合变量进行有很好的稳健性。
[0092]
具体地，构建算法如下：
[0093]
将n个土壤样品m个组分的土壤有机质含量浓度矩阵y和n个土壤样品p个波长点处吸光度矩阵x，分解成特征向量形式:
[0094]
y＝uq+f；
[0095]
x＝tp+e；
[0096]
其中u和t分别为n行d列的有机质含量浓度特征因子矩阵和吸光度特征因子矩阵，q为d
×
m阶浓度载荷矩阵，p为d
×
p阶吸光度载荷矩阵，f和e分别为n
×
m，n
×
p阶浓度残差阵和吸光度残差阵。
[0097]
根据特征向量的相关性分解y和x，建立估测模型:
[0098]
y＝xt+bq；
[0099]
b＝(u-ed)/t；
[0100]
其中ed为随机误差矩阵，b为d维对角回归系数矩阵。
[0101]
本发明还提出了一种系统用于实现所述的土壤有机质测定方法，包括：样品采集单元、高光谱数据测量单元、数据库、处理器和模型构建单元；
[0102]
样品采集单元，用于依据土壤类型和位置采集土壤样品，并采用gps记录采样位置的位置信息。
[0103]
所述高光谱数据测量单元，用于对采集的土壤样本进行反射光谱测量，获得土壤样本的反射光谱数据，并存储在所述数据库中；
[0104]
优选地，采用地物光谱仪进行土壤样本的光谱反射率数据测量。地物波谱仪是测量地表植被、农作物、土壤、岩石、水体等地物光谱的光电仪器，量测的物理内容有亮度、照度反射率、分光反射率、光谱分布等项目,是提取地物光谱特征信息的重要工具。
[0105]
优选地，采用specview软件对高光谱数据测量单元进行镜头校准-白板校准-大气校正，得到反射光谱数据。其中，specview软件针对机载高光谱成像系统，用来实现对光谱的控制和校正。
[0106]
所述处理器，用于剔除反射光谱数据中的异常值，对剔除所述异常值后的土壤样
本的反射光谱进行变换，对建模需要的特征波段进行选择；
[0107]
所述模型构建单元，用于构建特征波段的反射光谱数据与有机质含量之间的估测模型。
[0108]
本发明应用时，能快速的分析土壤有机质含量，低成本，且能提升分析精度，对于减少数据量，尤其是去除冗余信息起到了很好的作用，能够快速估算土壤中的有机质含量。
[0109]
以上所描述的装置实施例仅仅是示意性的，其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。本领域普通技术人员在不付出创造性的劳动的情况下，即可以理解并实施。
[0110]
通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到各实施方式可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件。基于这样的理解，上述技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在计算机可读存储介质中，如rom/ram、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行各个实施例或者实施例的某些部分所述的方法。
[0111]
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。