一种基于集合经验模态分解的耕地土壤肥力水平监测方法

技术领域

本发明涉及土壤肥力水平监测技术领域，特别是涉及一种基于集合经验模态分解的耕地土壤肥力水平监测方法。

背景技术

传统非遥感领域的土壤肥力水平监测方法普遍采用试验田的手段，通过定期进行土壤采样化验及农作物产量测算等途径对土壤肥力水平变化进行监测。这种方法耗费的人力物力资源多，且监测结果受空间和时间的约束较大，难以实现大范围动态监测。

相比之下，遥感技术在便捷性及时效性等方面都有着极大的提高。王之杰等利用不同氮素水平下冬小麦冠层叶片的光谱反射特点差异，对冬小麦冠层氮素分布对其籽粒蛋白质含量的影响进行了研究；田玉珂等利用TM遥感数据对土壤肥力各项指标的预测模型进行蒙特卡罗模拟修正，获得了有机质、全氮、碱解氮、有效磷以及速效钾的最佳修正光谱参数；王磊等利用无人机成像高光谱数据对土壤有机质(SOM)的敏感光谱特征进行了分析，提出了一种土壤SOM快速、准确、无损估测的监测方法，为遥感领域精准农业的发展以及成像高光谱技术的应用提供了一定的参考作用。

目前，这些基于遥感技术的土壤肥力水平监测方法虽然能够满足很多情况下的需求，但其仍存在着几个层面上的缺陷与不足：一是土壤中部分速效养分的含量差异在植被光谱特征上的显著性不足，导致其最终的反演精度受限；二是土壤表层与深层养分差异较大，且深层土壤的光谱数据在遥感影像中获取困难，导致数据准确性不足等问题；另外，一般农作物长势的遥感反演结果均包含着多种胁迫因素的共同影响，如何对土壤肥力胁迫进行有效而准确的甄别也是遥感技术在土壤肥力领域需要解决的重要难题。而农作物长势受到多种胁迫因素的影响较为复杂，单一胁迫作用甄别结果容易受到其他胁迫因素的影响；此外，影响土壤肥力水平的养分指标较多，遥感光谱数据对单一养分指标含量差异的响应不够敏感。

发明内容

本发明的目的是提供一种基于集合经验模态分解的耕地土壤肥力水平监测方法，在EEMD分解过程中对不同于土壤肥力时间尺度的胁迫作用进行的提取，并将复杂的土壤肥力系统用一个综合评价数值进行表达，使遥感光谱数据对于其数值差异的响应程度更为明显。

为实现上述目的，本发明提供了如下方案：

一种基于集合经验模态分解的耕地土壤肥力水平监测方法，包括如下步骤：

步骤1，通过高分一号遥感探测器采集耕地作物的光谱影像作为数据源；

步骤2，对数据源的光谱影像进行图像预处理，并将归一化植被指数NDVI作为耕地作物长势的反演参数，提取NDVI长时间序列；

步骤3，基于集合经验模态分解EEMD方法，对NDVI长时间序列进行不同时间尺度的分解，得到多个IMF分量，并将影响耕地作物长势的多种胁迫因素分解到不同的IMF分量中；

步骤4，结合野外实测数据和历史农情数据，对不同的IMF分量进行甄别，确定代表土壤肥力胁迫作用的组分序列，实现土壤肥力胁迫作用从多种复合胁迫作用结果中的甄别提取。

进一步的，所述方法还包括：结合测得的耕地土壤养分指标数据，从统计学角度对甄别提取结果进行验证，验证所述方法对土壤肥力胁迫作用进行甄别提取的可行性与准确性。

进一步的，所述步骤3，基于集合经验模态分解EEMD方法，对NDVI长时间序列进行不同时间尺度的分解，得到多个IMF分量，并将影响耕地作物长势的多种胁迫因素分解到不同的IMF分量中，具体包括：

采用集合经验模态分解EEMD将NDVI长时间序列按照年内、年间、年际三个时间尺度进行分解，得到多个IMF分量；

将水分、病虫害、土壤肥力以及重金属污染多种胁迫因素分解到不同的IMF分量中。

进一步的，采用集合经验模态分解EEMD将NDVI长时间序列按照年内、年间、年际三个时间尺度进行分解，得到多个IMF分量，其中，所述集合经验模态分解EEMD进行分解的过程包括：

对一组NDVI长时间序列的原始序列，加入高斯白噪声；

采用三次样条插值法得到加入高斯白噪声后的原始序列的上下包络线及均值；

减去原始序列的均值，得到新序列；

基于新序列得到一个IMF分量；

重复上述过程，进行多次分解，得到多个IMF分量。

进一步的，所述步骤4中，结合野外实测数据和历史农情数据，对不同的IMF分量进行甄别，确定代表土壤肥力胁迫作用的组分序列，实现土壤肥力胁迫作用从多种复合胁迫作用结果中的甄别提取，具体包括：

提取出代表水分、病虫害这类瞬时胁迫作用的年内波动组分序列；

提取代表耕地作物固有生长趋势的年间波动组分序列；

提取代表土壤肥力以及重金属污染这类环境胁迫作用的年际波动组分序列。

进一步的，所述结合测得的耕地土壤养分指标数据，从统计学角度对甄别提取结果进行验证，验证所述方法对土壤肥力胁迫作用进行甄别提取的可行性与准确性，具体包括：

在SPSS软件环境中对有机质、全氮、有效磷、速效钾四项养分指标数据测定结果进行相关性检验，得到相关性矩；

根据相关型矩阵对四项养分指标数据进行主成分分析，得到土壤养分总体水平关于四项养分指标的线性表达式；

将甄别提取结果与所述土壤养分总体水平关于四项养分指标的线性表达式采用最小二乘法进行线性拟合，基于拟合结果验证所述方法对土壤肥力胁迫作用进行甄别提取的可行性与准确性。

进一步的，根据相关型矩阵对四项养分指标数据进行主成分分析，得到土壤养分总体水平关于四项养分指标的线性表达式，具体包括：

根据相关型矩阵对四项养分指标数据进行主成分分析，选择前三项主成分作为评价土壤肥力水平的一组新变量，如下式所示：

F＝0.5759F

式中，F表示土壤养分总体水平；F

基于新变量新变量关于四项养分指标的线性表达式，得到土壤养分总体水平关于四种土壤养分指标的线性表达式；

新变量关于四项养分指标的线性表达式如下所示：

F

F

F

式中，X

土壤养分总体水平关于四种土壤养分指标的线性表达式如下：

F＝0.359X

进一步的，所述耕地作物选择冬小麦。

根据本发明提供的具体实施例，本发明公开了以下技术效果：本发明提供的基于集合经验模态分解的耕地土壤肥力水平监测方法，通过集合经验模态分解方法—EEMD，对冬小麦NDVI长时间序列进行分解，获取年内、年间、年际三个时间尺度的胁迫作用曲线，并结合野外实测数据及统计学方法确定了代表土壤肥力胁迫作用的组分，实现土壤肥力胁迫的甄别与提取；其中，EEMD方法依照序列本身的时间尺度进行分解，无需借助任何基函数，相较于传统时域分析方法，EEMD在处理非线性非平稳序列过程中具有明显的优势，理论上能够实现任何信号的分解，同时，不同胁迫作用对于冬小麦长势影响也存在时间尺度上的差异，因此，以EEMD方法为主体的土壤肥力胁迫作用甄别提取方法能够有效获取单一胁迫因素对于冬小麦长势的影响，对土壤肥力水平的监测有很大的帮助。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例基于集合经验模态分解的耕地土壤肥力水平监测方法的流程示意图；

图2(a)为本发明实施例加入高斯白噪声的原始序列曲线图；

图2(b)为本发明实施例三次样条插值法得到的序列上下包络线及均值曲线图；

图2(c)为本发明实施例原始序列减去均值后的新序列的曲线图；

图2(d)为本发明实施例得到的第一个IMF分量的曲线图；

图2(e)为本发明实施例EEMD单次分解结果的曲线图；

图3为本发明实施例年内波动组分示意图；

图4为本发明实施例年间波动组分示意图；

图5(a)为本发明实施例季节平均模型示意图；

图5(b)为本发明实施例季节项分量示意图；

图6为本发明实施例年际波动组分示意图；

图7(a)为本发明实施例大厂县2015年降水量示意图；

图7(b)为本发明实施例大厂县2016年降水量示意图；

图7(c)为本发明实施例大厂县2017年降水量示意图；

图7(d)为本发明实施例大厂县2018年降水量示意图；

图7(e)为本发明实施例大厂县2019年降水量示意图；

图8(a)为本发明实施例陈府镇土壤肥力胁迫作用水平曲线图；

图8(b)为本发明实施例东马各庄土壤肥力胁迫作用水平曲线图；

图8(c)为本发明实施例兰庄户村土壤肥力胁迫作用水平曲线图；

图8(d)为本发明实施例梁庄村土壤肥力胁迫作用水平曲线图；

图8(e)为本发明实施例荣马坊村土壤肥力胁迫作用水平曲线图；

图8(f)为本发明实施例王必屯村土壤肥力胁迫作用水平曲线图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的是提供一种基于集合经验模态分解的耕地土壤肥力水平监测方法，在EEMD分解过程中对不同于土壤肥力时间尺度的胁迫作用进行的提取，并将复杂的土壤肥力系统用一个综合评价数值进行表达，使遥感光谱数据对于其数值差异的响应程度更为明显。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

如图1所示，本发明实施例提供的基于集合经验模态分解的耕地土壤肥力水平监测方法，包括如下步骤：

步骤1，通过高分一号遥感探测器采集耕地作物的光谱影像作为数据源；

步骤2，对数据源的光谱影像进行图像预处理，并将归一化植被指数NDVI作为耕地作物长势的反演参数，提取NDVI长时间序列；

步骤3，基于集合经验模态分解EEMD方法，对NDVI长时间序列进行不同时间尺度的分解，得到多个IMF分量，并将影响耕地作物长势的多种胁迫因素分解到不同的IMF分量中；具体包括：

采用集合经验模态分解EEMD将NDVI长时间序列按照年内、年间、年际三个时间尺度进行分解，得到多个IMF分量；

将水分、病虫害、土壤肥力以及重金属污染多种胁迫因素分解到不同的IMF分量中。

其中，所述集合经验模态分解EEMD进行分解的过程包括：

对一组NDVI长时间序列的原始序列，加入高斯白噪声；

采用三次样条插值法得到加入高斯白噪声后的原始序列的上下包络线及均值；

减去原始序列的均值，得到新序列；

基于新序列得到一个IMF分量；

重复上述过程，进行多次分解，得到多个IMF分量。

步骤4，结合野外实测数据和历史农情数据，对不同的IMF分量进行甄别，确定代表土壤肥力胁迫作用的组分序列，实现土壤肥力胁迫作用从多种复合胁迫作用结果中的甄别提取，具体包括：

提取出代表水分、病虫害这类瞬时胁迫作用的年内波动组分序列；

提取代表耕地作物固有生长趋势的年间波动组分序列；

提取代表土壤肥力以及重金属污染这类环境胁迫作用的年际波动组分序列。

所述方法还包括：结合测得的耕地土壤养分指标数据，从统计学角度对甄别提取结果进行验证，验证所述方法对土壤肥力胁迫作用进行甄别提取的可行性与准确性，具体包括：

在SPSS软件环境中对有机质、全氮、有效磷、速效钾四项养分指标数据测定结果进行相关性检验，得到相关性矩；

根据相关型矩阵对四项养分指标数据进行主成分分析，得到土壤养分总体水平关于四项养分指标的线性表达式；

将甄别提取结果与所述土壤养分总体水平关于四项养分指标的线性表达式采用最小二乘法进行线性拟合，基于拟合结果验证所述方法对土壤肥力胁迫作用进行甄别提取的可行性与准确性。

本发明实施例中以冬小麦为目标作物，从冬小麦长势长时间序列入手，选取高分一号16m分辨率高光谱影像作为数据源，归一化植被指数NDVI作为冬小麦长势的反演参数，基于一种时域信号分析方法——集合经验模态分解EEMD，对长时间序列进行不同时间尺度的分解，将水分、病虫害、土壤肥力以及重金属污染等多种胁迫因素分解到不同频率的IMF分量中。结合相关理论以及定期开展的野外实测工作数据结果，对不同的IMF分量进行甄别，提取出代表了水分、病虫害等瞬时胁迫的年内波动组分序列、代表冬小麦固有生长趋势的年间波动组分序列以及代表土壤肥力胁迫作用的年际波动组分序列，有效实现了土壤肥力胁迫作用从多种复合胁迫作用结果中的甄别提取。结合测得的土壤养分指标数据，从统计学角度对分析结果进行验证，证明了此方法对土壤肥力胁迫进行甄别提取的可行性与准确性，并进一步确定了基于遥感影像数据参数的土壤肥力水平定量评价模型。

冬小麦长势时间序列构建：获取若干典型样区自2015年以来的高分一号卫星影像数据，选择NDVI作为冬小麦长势的反演参数，获取复合胁迫因素作用下的冬小麦长势长时间序列曲线。

EEMD分解过程实现：对于一组NDVI长时间序列，其在EEMD算法中的单次分解流程如图2(a)至图2(e)所示。其中三次样条插值法求信号上下包络线的过程可用spline函数实现，为与实际数据时间尺度相匹配，差值过程的采样频率设置为12，采样周期为1/12，即每一年内进行12次采样。

将EEMD分解过程重复N次即得到最终EEMD的分解结果。通常情况下，综合考虑对计算效率以及结果精度的影响，N的次数一般设置在50至100次之间。在这一区间内N的数值变化对结果精度的影响较小，为保证研究的科学性，本实验中设置N为100次。

在EEMD分解的整体流程中，添加高斯白噪声信号的强弱对最终结果的精度起到了决定性的影响。为选择合适的数值，研究中规定μ为加入的高斯白噪声信号标准差与原始序列间的比值，μ的数值过小则加入的噪声信号对原始序列的影响过小从而难以消除模态混叠的影响；μ过大则会影响分解结果的真实性。在一般的EEMD分解研究中，μ的数值通常设置在0.1-0.4之间，该数值的选取应结合实际研究内容考虑。

表1不同μ值下的年间波动组分相关系数

对于冬小麦而言，其不同年份间的固有生长趋势应大致接近，后文中将提到分解结果中的年间波动组分代表了冬小麦的固有生长趋势，因此通过计算比较不同μ值下各年份间年间波动组分的相关系数即可确定本实验中最佳的μ值。表1给出了不同μ值下EEMD分解结果年间波动组分各年份间的相关系数，可以看出，当μ＝0.3时，平均相关系数达到最大值0.849，因此本研究中选取μ＝0.3。

EEMD分解结果的组合：经过EEMD分解的NDVI时间序列根据频率差异被分解为两项周期大于年的组分、一项周期约等于年的组分、一项周期小于年的组分、一项无规则的组分以及一个残差项，根据研究需求对这些不同的分量进行组合得到三个不同时间尺度上的组分分量。通过后续甄别工作能够将这些组分归纳为以时间尺度分类的年际胁迫、年间胁迫以及年内胁迫(瞬时胁迫)，同时根据不用胁迫因素的实际作用类型，可对胁迫因素进行有效的甄别。

①年内波动组合的合成：EEMD分解结果中的分量IMF1以及IMF2的频率较高，二者的周期均小于1年，因此可以判定IMF1以及IMF2为信号的年内波动组分，这一部分分量代表着植物生长周期中受到的瞬时胁迫影响作用以及信号处理过程中产生的误差结果。其中，由于信号处理过程产生的误差组分为噪声信号，这一组分不应参与模型最终分解结果的合成。因此，IMF2代表了冬小麦长势长时间序列分解结果中的年内波动组分，如图3所示。

②年间波动组分的合成：年间波动组分又可称为长时间序列分解结果中的季节项，其典型特征是信号组分具有一个明显的周期。根据农作物的实际物候参数不同，这个周期也会出现一定的变化。通常情况下，这个时间周期的长度是1年。但由于实际生产过程中耕地的利用方式以及农作物的耕作方式的影响，这个年间波动组分会呈现出周期在1年左右变动的趋势。为确定研究区域内冬小麦长势的年间波动组分，在此引入季节平均模型的方法以确定冬小麦长势的年间波动组分，其流程如下：

将长时间序列所包含年份中的每年同月的数据进行平均值求取，并重复6年以获取平均季节数据；

此平均季节数据代表了冬小麦长时间序列中所有周期特征的平均波动结果，能够反映出不同年份间冬小麦生长所共有的周期特征。利用EEMD算法对平均季节数据进行分解，得到代表该基本波动的IMF分量，如图5(a)和图5(b)所示；

将原始序列EEMD分解结果中所有与上一步中得到的IMF分量周期一致的分量进行组合，即可得到冬小麦长势的年间波动组分，如图4所示。

③年际波动组分的合成：年际波动组分又称为长时间序列分解结果的趋势项，它代表了信号在一个长时间尺度上的振动趋势。对于农作物而言，趋势项代表的意义是其受到时间尺度大于季节项的胁迫因素影响时所产生的响应。因此，将振动频率小于年频率的若干IMF分量进行合成，得到冬小麦长势长时间序列分解结果中的年际波动组分，如图6所示。

土壤肥力胁迫甄别与提取：

①年内波动组分：这部分是EEMD将原始序列分解得到的高频分量。高频分量的时间分辨率较高，频率分辨率较低，能够对原始序列的高频振荡部分进行很好地表现，可以检测信号夹带的瞬时异常现象。这些信息的特征表现为年内波动杂乱，振幅规律性不强。对相邻年份同一月份的信号数值进行相关性分析能够较好地对这一特征进行验证，在此选择3月、4月、5月、6月四个冬小麦生长较为旺盛的月份进行相关性分析，结果如表2所示。可以看出，相邻年份间的同一月份数值相关性较小，表明曲线的规律性较弱，符合上述猜想。因此，可以认为冬小麦长势长时间序列的年内波动组分代表了冬小麦受各种瞬时胁迫因素胁迫的作用结果。

表2年内波动组分相关性分析结果

相对而言，水分胁迫、病虫害胁迫等胁迫对于农作物长势的影响相较于土壤肥力、重金属污染等环境因素的响应更快，影响程度更大。而年内波动组分与其他组分相比，在曲线上表现为振动频率更大且振幅更大。这一特点能够反映出水分胁迫、病虫害胁迫的胁迫信息。因此，假定信号的高频分量代表了冬小麦受水分、病虫害等瞬时胁迫作用的响应结果。

图7(a)至图7(e)及表3给出了廊坊市大厂县2015年-2019年的年降水量。其中，2018年的年降水量远高于其他年份。根据冬小麦的耕作技术，较多的降水和连续阴雨天气会对冬小麦的长势和最终产量造成一定的影响，如影响扬花授粉、容易诱发病虫害、降低小麦粒重等。此外，从图中可以看出，大厂县降水主要集中于每年的7-9月份，结合该地区冬小麦耕作“春灌1水”的灌溉方式，水分胁迫作用在一年应在上半年和下半年各出现一个峰值。代表年内波动组分的曲线中，2018年的曲线振幅明显小于其他年份，并且曲线在每年的4月份以及9月份左右出现了两个波峰。这些特性均符合水分胁迫对冬小麦长势的影响，因此，可以认为信号的高频分量也就是年内波动组分能够代表冬小麦受水分、病虫害等瞬时胁迫作用影响的结果，将这些胁迫称为年内胁迫。

表3大厂县2015-2019年降水量

②年间波动组分：年间波动组分反映了农作物的固有生长趋势，可以近似为一个稳定的周期波动信号，如图5所示。对于试验区内的冬小麦，其在一年中有一个高峰振荡，出现在每年的5月至6月左右。理论上讲，农作物在不受外界胁迫因素影响时，不同年份的年间胁迫组分应呈现出相似的趋势，但在实际研究中，耕地作物类型的变更以及实验过程中数据获取的误差会导致年间波动组分在不用年份存在一定的差异，这种差异主要体现在波峰与波谷的数值差异上，但在时间尺度上峰值对应的时间点在不同年份内大体是一致的。

③年际波动组分：根据上文的研究结果，冬小麦长势的长时间序列由年内波动组分、年间波动组分、年内波动组分构成。在剔除代表水分胁迫、病虫害胁迫等瞬时胁迫的年内波动组分以及代表冬小麦固有生长趋势的年间波动组分后，剩余的年际波动组分可代表持续的环境胁迫对农作物长势的影响，包括土壤肥力胁迫、重金属胁迫等因素。从图6中的曲线可以看出，这种胁迫是长期而平缓的。土壤肥力胁迫的作用形式主体为农作物对土壤中养分的吸收，土壤中各种养分的含量通过不同途径的生物界循环在整体上是保持一定幅度内的平衡的，这一特征符合年际波动组分曲线平稳且振幅较小的特点。因此，可以认为年际波动组分中包含着土壤肥力胁迫作用的相应结果。同时，根据野外实测工作中利用手持XRF分析仪对几个典型样区耕地土壤中重金属元素含量的测定结果，如表4所示，所有样区耕地土壤中的重金属元素含量均在正常水平范围内，可以认为该地区受到重金属胁迫作用的影响很弱。因此，年际波动组分能够较好地反应研究区域内土壤肥力胁迫作用对于农作物长势的影响。如图8(a)至图8(f)给出了表4中几个样区的土壤肥力胁迫作用水平曲线。

表4土壤重金属污染元素含量测定结果(2020.12)

分析结果验证与土壤肥力水平综合评价模型构建：选择有机质、全氮、有效磷、速效钾四项典型养分指标进行含量测定，基于测定结果构建土壤肥力综合水平的定量评价模型。

①在SPSS软件环境中对表5所示四项养分指标测定结果进行相关性检验，得到相关性矩阵如表6所示；

表5大厂县部分地区土壤养分含量测定结果

表6土壤养分指标相关性检验结果

表7主成分分析结果

②根据相关型矩阵对四项养分指标数据进行主成分分析，结果如表7所示。前三项主成分的累计贡献率打到99.16％，基本能够代表总体水平。因此选择前三项主成分作为评价土壤肥力水平的一组新变量，如下式所示：

F＝0.5759F

式中，F表示土壤养分总体水平；Δ(F)表示余项。

此外，新变量关于四项养分指标的线性表达式如下：

F

F

F

式中，F

将两组表达式相结合，即可得到土壤养分总体水平关于四种土壤养分指标的线性表达式，如下式所示：

F＝0.359X

基于上式得到几个典型样区的土壤肥力综合水平数值如表8所示。

表8大厂县样区土壤肥力水平

③与EEMD分解甄别结果进行拟合。表9列出了几个典型样区土壤肥力水平与EEMD分解的土壤肥力胁迫组分数据，利用最小二乘法对二者进行拟合，结果如表10所示。可以看出，虽然IMF

F＝0.1559+7669.2IMF-6.5*10

利用最小二乘法对数据进行线性拟合的方法，在实现过程上较为简单，且拟合结果的相关系数达到了0.9741，能够较好地反映出原始数据的变化水平。因此，此过程实现了土壤肥力水平基于遥感影像数据的定量评价。

表9拟合数据

表10几种不同拟合方式的拟合结果

本发明提供的基于集合经验模态分解的耕地土壤肥力水平监测方法，通过集合经验模态分解方法—EEMD，对冬小麦NDVI长时间序列进行分解，获取年内、年间、年际三个时间尺度的胁迫作用曲线，并结合野外实测数据及统计学方法确定了代表土壤肥力胁迫作用的组分，实现土壤肥力胁迫的甄别与提取；本方法在EEMD分解过程中对不同于土壤肥力时间尺度的胁迫作用进行提取和甄别，并将复杂的土壤肥力系统用一个综合评价数值进行表达，使遥感光谱数据对于其数值差异的响应程度更为明显。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。