CN114266972A - 基于遥感植被指数提取植被物候指标的方法及装置 - Google Patents

Abstract

本发明提出了一种基于遥感植被指数提取植被物候指标的方法及装置。从以下几个方面改进了植被物候指标提取精度：(1)采用权重调整函数使植被指数摆脱了多种光学噪点(如云、雪、气溶胶、阴影等)的干扰；(2)采用粗略拟合与精细拟合相结合，使植被物候指标提取在多生长季或植被类型发生变化的区域表现优异。粗略拟合用于捕捉植被周期性信号，并划分生长周期；精细拟合用于在每个生长周期中捕捉植被生长周期关键点的迅速变化；(3)融合了多种权重调整、粗略拟合、精细拟合以及植被物候提取方法，对提取的植被物候指标的不确定性进行了量化。本发明可应用于森林、草地、农田等生态系统生长周期关键性转折点的动态监测，指导生产实践。

Description

基于遥感植被指数提取植被物候指标的方法及装置

技术领域

本发明涉及遥感和植被生态学领域，具体涉及一种基于遥感植被指数提取植被物候指标的方法及装置。

背景技术

遥感植被指数数据为全球大范围、长时间植被活动和植被物候指标监测提供了一种有效手段。然而，这些数据往往受到云、雪、气溶胶、阴影等噪点的污染，使得物候数据的准确提取具有挑战性。即使是众多卫星产品中质量较好的MODIS遥感植被指数，其中未被污染的观测仍仅有44％。如何有效地进行植被指数除噪，是植被物候指标提取要面临的主要难题。

目前国际上最流行的植被物候指标提取软件有TIMESAT、greenbrown和phenopix。TIMESAT最早在2002年提出，是使用最广泛的植被物候指标提取软件之一。但是其方法单一，生长周期划分采用傅里叶变化识别的静态生长周期；精细拟合仅包含AG和Zhang，由于二者是分段函数，返青期与枯黄期的连接点往往不连续，植被物候指标提取仅采用单一的阈值法。

近20年涌现了诸多生长周期、精细拟合和植被物候指标提取方法，但由于TIMESAT是非开源软件，这些方法无法在TIMESAT中实现。phenopix是2016年提出的另一个通用的植被物候指标提取软件，phenopix继承了greenbrown的基础逻辑，融合了近20年最先进的植被指数除噪和植被物候指标提取方法。但phenopix主要为高频相机监测植被绿度指数而设计，没有处理遥感指数中常见的污染的模块。

发明内容

本发明的目的在于：(1)融合近几十年最先进的权重调整方案、生长周期划分、精细拟合和植被物候指标提取方法，提高植被物候指标的提取精度；(2)在进行物候指标提取时采用多种函数，以对提取物候指标的不确定性进行量化；(3)提供灵活的输入，可以支持不均匀间隔的时间序列以及包含缺失值的时间序列。

为实现上述技术的目的，本发明采用如下技术方案：

一种基于遥感植被指数提取植被物候指标的方法，包括以下步骤：

步骤1：获取输入数据，包括：植被指数时间序列、日期信息和质量控制波段；

步骤2：根据所述质量控制波段，赋予所述植被指数时间序列中不同质量的数据点不同的权重；

步骤3：识别所述遥感植被指数时间序列的钉值，将钉值看作缺失值，对缺失值采用相邻日期线性插值进行填补，得到预处理后的植被指数时间序列；

步骤4：对所述预处理后的植被指数时间序列进行粗略拟合，得到粗略拟合后的遥感植被指数时间序列；

步骤5：对所述粗略拟合后的植被指数时间序列中的每个数据点调整权重；

步骤6：重复运行步骤4粗略拟合和步骤5调整权重，直至最大循环次数，输出粗略拟合除噪后的植被指数序列；

步骤7：根据步骤6输出的粗略拟合除噪后的植被指数序列，以两个谷值确定一个生长周期为原则，划分生长周期；

步骤8：在每个生长周期内，对粗略拟合除噪后的植被指数序列进行精细拟合，得到精细拟合后的遥感植被指数时间序列；

步骤9：对精细拟合后的遥感植被指数时间序列中的每个数据点调整权重；

步骤10：重复运行步骤8精细拟合和步骤9调整权重，直至最大循环次数，输出精细拟合除噪后的植被指数序列；

步骤11：根据步骤9输出的精细拟合除噪后的植被指数序列，提取物候指标。

优选地，步骤2中，遥感植被指数时间序列根据质量控制波段划分为三类，良好、一般和较差；这三类数据的权重分别赋值为最大权重、中等权重和最低权重。

优选地，步骤5和步骤9中，通过调整权重以降低噪点干扰。本发明提供了三种权重调整方案：wBisquare、wTSM、wChen。来调整时间序列中每个点的权重，增加被低估点的权重以逼近上包络线，降低离群点的权重以减少离群点在拟合过程中的干扰，以减轻噪点污染。

优选地，步骤3具体包括：

步骤3.1：若一个数据点的y_i与y_mov,i的绝对差大于整个时间序列标准差的2倍，那么这个点被标记为一个钉值，将所述钉值作为缺失值；其中y_i是输入植被指数时间序列的第i个观测值，y_mov,i是相应的3点滑动平均值；

步骤3.2：利用相邻日期线性插值对缺失值进行填补；

步骤3.3：确定填补后的植被指数时间序列的有效值域：根据质量较好的植被指数观测值，确定植被指数的有效值域(lower andupperboundary，简称ylu)。

优选地，步骤4中，采用的粗略拟合函数包括：加权Savitzky-Golay、加权HANTs和加权Whittaker。

优选地，步骤7具体包括：

步骤7.1：根据粗略拟合除噪后的植被指数序列，检测局部最大值或最小值，并将其作为潜在峰值或谷值；

步骤7.2：判断所述潜在峰值或谷值中是否存在由于局部细微波动引起的虚假峰值或谷值，若存在，剔除该虚假峰值或谷值；

步骤7.3：通过两个波谷定义了一个生长周期，生长周期的编号标记为year_k，其中year为生长高峰所在的年份，k为这一年生长季节的顺序。

优选地，步骤7.2具体包括：

步骤7.2.1：若max(h_peak,L,h_peak,R)<r_maxA或min(h_peak,L,h_peak,R)<r_minA，则认为是虚假值，并进行剔除，h_peak,L、h_peak,R的表达式为：

其中，y_peak为生长周期的峰值；y_trough,L、y_trough,R为生长周期峰值左侧与右侧的谷值；h_peak,L、h_peak,R分别为生长周期峰值到左右谷值的高度差；A为输入植被指数时间序列y的幅值；r_max和r_min是两个待定参数，一般情况下，相对较小的阈值(r_max＝0.2和r_min＝0.05)已经可以消除大多数虚假的峰值和谷值；相反，更大的阈值可能会掩盖住真实的生长周期峰值；

步骤7.2.2：若max(h_trough,L,h_trough,R)<r_maxA，亦认为是虚假谷值，并进行剔除，h_trough,L和h_trough,R的表达式为：

其中，y_trough为生长周期的谷值；y_peak,L、y_peak,R为生长周期谷值左侧与右侧的峰值；h_trough,L和h_trough,R分别为生长周期谷值到左右峰值的高度差；

步骤7.2.3：若相邻两个峰或谷的时间差小于len_min，则认为是虚假峰值或谷值，len_min为设定的最小时间差。

优选地，步骤8中，采用的精细拟合函数包括：非对称高斯分布函数AG和多种双侧逻辑斯蒂函数；

多种双侧逻辑斯蒂函数包括Zhang、Beck、Elmore和Gu。

优选地，步骤8还包括：

根据精细拟合函数中各参数的生态意义确定各个参数的有效值域，在参数求解过程中进行约束；

在每个生长周期精细拟合过程中，将该生长季节前后的质量前n的观测加入进来，使不同生长周期之间过渡平滑。

优选地，步骤11中，采用的植被物候指标提取方法包括：阈值法TRS、最大斜率法、曲率曲线法和Gu法，在阈值法TRS中，生长季节开始和结束分别定义为植被指数超过阈值和返回值低于阈值的第一天，阈值通常设置为0.1、0.2或0.5。

为了实现上述的一种基于遥感植被指数提取植被物候指标的方法，本发明还提供了一种基于遥感植被指数提取植被物候指标的装置，包括以下模块：

获取模块，用于获取输入数据，包括：植被指数时间序列、日期信息和质量控制波段；

权重初始化模块，用于根据所述质量控制波段，赋予所述植被指数时间序列中不同质量的数据点不同的权重；

预处理模块，用于识别所述遥感植被指数时间序列的钉值，将钉值看作缺失值，对缺失值采用相邻日期线性插值进行填补，得到预处理后的植被指数时间序列；

粗略拟合模块，用于对所述预处理后的植被指数时间序列进行粗略拟合，得到粗略拟合后的遥感植被指数时间序列；

权重调整模块，用于对所述粗略拟合后的植被指数时间序列中的每个数据点调整权重；

循环模块，用于重复运行粗略拟合模块和权重调整模块以调整权重，直至最大循环次数，输出粗略拟合除噪后的植被指数序列；

生长周期划分模块，用于根据输出的粗略拟合除噪后的植被指数序列，以两个谷值确定一个生长周期为原则，划分生长周期；

精细拟合模块，用于在每个生长周期内，对粗略拟合除噪后的植被指数序列进行精细拟合，得到精细拟合后的遥感植被指数时间序列；

所述权重调整模块，还用于对精细拟合后的遥感植被指数时间序列中的每个数据点调整权重；

所述循环模块，还用于重复运行精细拟合模块和权重调整模块以调整权重，直至最大循环次数，输出精细拟合除噪后的植被指数序列；

物候指标提取模块，用于根据输出的精细拟合除噪后的植被指数序列，提取物候指标。

本发明提供的技术方案具有以下有益效果：

本发明从以下几个方面提高了植被物候指标提取精度：(1)采用权重调整函数使植被指数摆脱了多种光学噪点(如云、雪、气溶胶、阴影等)的干扰；(2)采用粗略拟合与精细拟合相结合，使植被物候指标提取在多生长季或植被类型发生变化的区域表现优异。粗略拟合用于捕捉植被周期性信号，并划分生长周期；精细拟合用于在每个生长周期中捕捉植被生长周期关键点的迅速变化；(3)融合了多种权重调整、粗略拟合、精细拟合以及植被物候提取方法，对提取的植被物候指标的不确定性进行了量化。本发明可应用于森林、草地、农田等生态系统生长周期关键性转折点的动态监测，指导生产实践。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明实施例中一种基于遥感植被指数提取植被物候指标的方法的流程图；

图2是本发明实施例中生长周期的划分方法；

图3是本发明实施例中在2010年1月至2016年12月期间，使用MOD13A1.006步长为16日的EVI为输入，在“CA-NS6”站点进行粗略拟合和生长周期划分的结果；

图4是本发明实施例中逻辑斯蒂函数表3中参数k的值域说明；

图5是在2010年1月至2016年12月期间，使用MOD13A1.006步长为16日EVI作为输入，在“CA-NS6”站点通过五种精细拟合方法进行精度拟合的结果；

图6是本发明实施例中没有参数边界约束的情况下精细拟合的结果；

图7是本发明实施例中使用2010-2016年在“CA-NS6”站点MOD13A1.006的16日EVI作为输入，本发明方法与TIMESAT、phenopix之间的比较结果；

图8是本发明实施例中植被物候提取方法的示意图；

图9是本发明实施例中对“CA-NS6”站点利用AG法进行精细拟合得到的植被物候指标；

图10是本发明实施例中2015年河南省第1个生长周期的空间分布；

图11是本发明实施例中2015年河南省第2个生长周期的空间分布；

图12是本发明实施例中2015年河南省第3个生长周期的空间分布。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

如图1所示，本实施例提供了一种基于遥感植被指数提取植被物候指标的方法，主要包括以下步骤：

步骤一：准备输入数据和预处理

本实施例中使用MOD13A1.006增强型植被指数(enhanced vegetation index，简称EVI)来阐述本发明的使用方法。MOD13A1.006 EVI是地表植被动态研究中最常用到的植被指数之一，其时空分辨率为500m和16天。采用CA-NS6站点2010-2016年MOD13A1.006 EVI数据用于后续分析。该站点位于-98.9644°E，55.9167°N，是典型的寒带针叶林。CA-NS6站点具有明显的周期性，且频繁受到云、雪的污染。该站点的植被指数数据可以有效地测试本发明除噪的能力。

为提取植被物候指标，本发明需要植被指数时间序列、相应的观测日期、以及权重。输入数据已经包括植被指数时间序列、日期信息、质量控制波段。权重则根据质量控制波段获得。根据质量控制波段，将不同质量的数据点赋值以不同的权重。质量较好的点，赋予其较高的权重；反之则赋予其较低的权重。权重的值域在[0,1]。

具体地，输入的植被指数时间序列根据其质量波段QC划分为三类，良好、一般、较差；这三类数据的权重分别赋值为最大权重(w_max)、中等权重(w_mid)、最低权重(w_min)。w_max、w_mid、w_min一般为1.0、0.5和0.2。被云、雪污染的植被指数认为其质量较差，赋值为最低权重w_min，但w_min应该大于0。

其次，对输入的植被指数时间序列进行如下预处理：

(1)识别钉值；如果一个点的y_i与y_mov,i的绝对差大于整个时间序列标准差的2倍，那么这个点被标记为一个“钉值”，其中y_i是输入植被指数时间序列的第i个观测值，y_mov,i是相应的3点滑动平均值。

(2)缺失值插补；利用相邻日期线性插值对缺失值进行填补，所有的“钉值”采用与缺失值相同的方法进行填补。

(3)确定植被指数有效值域；根据质量较好的植被指数观测值，确定植被指数的有效值域(lower and upper boundary，简称ylu)。

步骤二：粗略拟合

对于预处理后得到的植被指数时间序列进行粗略拟合与生长周期划分。三种粗略拟合函数如表1所示，其中，wSG(加权Savitzky-Golay)适用于植被生长周期信号波动较大的植被；wHANTs(加权HANTs)适用于具有多个生长周期，且周期信号稳定的植被；wWHIT(加权Whittaker)不仅适用于具有多个生长周期的植被，还适用于生长周期动态变化的情况。相对于另外两种粗略拟合函数，对于异常值wWHIT具有更好的稳定性。即使植被指数有超过30％的点被污染，wWHIT也能捕获植被的季节信号。在使用粗略拟合时，用户应仔细选择参数。

表1三种粗略拟合函数

其中nf、frame、λ分别是wHANTS(加权HANTs)、wSG(加权Savitzky-Golay)、wWHIT(加权Whittaker)的参数；nf为傅里叶变换中的谐波数；frame为wSG中滑动窗口的宽度；λ为wWHIT的粗糙度惩罚参数。wSG与wWHIT的最优参数比wHANTs的更加稳定。(1)frame为nptperyear*2/5(nptperyear为每年的点数)；(2)nf取决于季节的数量(nseason)，理想情况下应该为nseason的偶数倍(即2×nseason×m，m∈正整数)；(3)λ随输入的植被指数的时间尺度变化而变化，对于日尺度、8日尺度、16日尺度，λ分别为10000、15、5。

图3显示了本发明对数据进行处理后返回的粗略拟合的结果。水平的虚线是植被指数时间序列y的上下边界(ylu)；较细的折线为原始EVI时间序列；“iter1”和“iter2”是使用wWHIT进行粗略拟合的第一次迭代与第二次迭代的结果。圆点对应生长周期划分模块中检测到的峰(谷)值。图3的顶部显示了粗略拟合的性能指标，包括决定系数(R²)、Nash-Sutcliffe模型效率系数(NSE)、观测和模拟植被指数的变异系数(CV_obs、CV_sim)。其中，粗略拟合使用了更加稳定的wWHIT函数。从图3中可以看出，上下边界ylu可以有效地过滤掉非生长周期的异常值，进而得到更加合理的植被生长周期背景值。

在一些实施例中，粗略拟合函数还可以采用wHANTS或wSG，可以根据实际需求进行选择。

步骤三：生长周期划分

对粗略拟合得到的结果进行生长周期划分。图2是生长周期划分方法的具体说明。图3是该方法在CA-NS6站点的EVI数据上的具体应用。从图3可以看出，本发明可以有效地摆脱光学噪点的干扰(如图3中的2015年冬末)，稳健地捕捉到植被周期性信号，合理并正确地划分生长周期。

步骤四：精细拟合

基于划分好的生长周期，在每个生长周期中，对植被指数进行精细拟合，重建植被的日尺度时间序列，为后续的植被物候指标提取奠定好数据基础。本实施例中，给出五种精细拟合函数，如表2，具体包括非对称高斯分布函数(AG)和多种双侧逻辑斯蒂函数(包括Zhang、Beck、Elmore和Gu)。

表2提供的五种精细拟合函数

值得注意的是，在精细拟合过程中，需要对参数边界进行约束，表3为对精细拟合函数中参数的约束范围；

表3对精细拟合函数中参数的约束范围

pop代表代表生长周期峰值出现的时间；y_min(y_max)是植被指数时间序列的最小(最大)值；rsp(rap)是返青期和枯黄期的生长速率；t_min(t_max)是生长周期的开始(结束)时间；k₀是rsp的估计值；T_half是生长周期长度的一半；A是生长周期的振幅；k₀＝4/T_half*2.67，T_half＝(t_max-t_min)/2，A＝y_max-y_min，Δy＝0.1A，Δt＝T_half/2。

本实施例中，用最具代表性的Zhang逻辑斯蒂函数推导出参数k₀的范围。其中，y为植被指数；mn与mx是y的最小值与最大值；rsp是返青期的生长速率；sos是生长周期开始的时间；t是y的对应日期。由于返青期和枯黄期植被生长曲线近似对称，此处仅选取返青期(t≤t₀)的时间序列进行分析。其中，t₀为生长周期返青期和枯黄期的转折点。为了简化方程(A1)，将y归一化为p，p＝(y-mn)/(mx-mn)，将(A1)重新整理可得：

p’是p的一阶导数，易知p在t＝sos处得到最大值，其中p’(t＝sos)＝rsp/4。为求得参数k₀的范围，问题可以转化为使用p估计参数rsp的合理范围。rsp是返青期的生长速率；sos是生长周期开始的时间；t是p的对应日期。

如图4(a)，令p₁＝p，p₂＝1-p₁，因此p₂＝1-p，且

则O点(sos，0.5)的生长率可由斜率k(x)得出估计值：

其中，rsp是返青期的生长速率；sos是生长周期开始的时间；t₁与t₂是p₁与p₂的对应日期。由于p’(t＝sos)＝rsp/4，因此rsp＝4k(0)。需要用归一化p估计k(0)。如果仅是相邻的两个观测点，其被污染的概率较高，其一阶导也可能存在问题。因此，本发明通过建立平均增长率与最大增长率之间的关系来估计rsp。由于p在头、尾趋近于无穷，这有碍于平均增长率的估计。采用k(0.99)来粗略估计，其中k(0.99)≈1/T_half(T_half为生长周期长度的一半)。根据极限定理，可以得出k(0)/k(0.99)≈2.67(如图4(b))。最后，通过公式(A6)与(A7)得到rsp的估计值，即表4中的参数的范围。

p'(t＝sos)＝k(0)＝rsp/4＝2.67k(0.99) (A6)

图5为AG、Beck、Elmore、Gu、Zhang在CA-NS6站点精细拟合的结果。垂直虚线代表生长周期的开始(或结束)。需要注意，上一个生长周期的结束时间通常与下一个生长季节的开始时间重叠。参数边界约束有利于提高精细拟合的稳定性。如果不进行参数边界约束，精细拟合很容易变得不收敛(如图6中的AG)，精细拟合函数参数的物候意义将丧失。

精细拟合是粗略拟合、生长周期划分等方法表现好坏的集中体现，它也直接决定着后续植被物候指标提取的精度。图7进一步探讨了本发明与现有方法(TIMESAT、phenopix)在精细拟合上的差异。其中，TIMESAT使用了非对称高斯(AG)；phenopix由于没有AG函数，使用了Beck函数(见表2)。结果表明，三种方法在大多数情况下都能成功地重建植被指数。但当输入时间序列受到污染或在非生长周期时，重建效果会发生很大变化。phenopix由于没有权重调整模块，所以对异常值更加敏感，如图7中的2014年秋季。TIMESAT的切断函数以及全局函数让重建的时间序列稳定且平滑，但TIMESAT与phenopix的背景值在不同年份之间存在很大的波动，如图7中的2013-2014年和2015-2016年。当植被类型未发生变化时，认为植被非生长季背景值近似可看作固定值，不随年份变化而改变。不可靠的非生长季背景值可能会导致植被物候指标提取过程中出现明显的误差，尤其是TRS法。这方面，本发明在每个生长周期精细拟合过程中，将该生长季节前后的n个质量较好的观测也加入进来，使不同生长周期之间的过渡平滑，更合理地确定非生长季背景值，取得了更好的表现。

步骤五：植被物候指标提取

最后，本发明利用精细拟合重建的日尺度时间序列，提取植被物候指标。图8是植被物候提取方法的示意图。包含如下方法：(a)阈值法(TRS)。其中f(t)为原始植被指数时间序列，f(t)_ratio＝(f(t)-max)/(max-min)，f(t)_ratio为归一化植被指数时间序列，阈值在[0,1]。生长周期开始(SOS)和生长周期结束(EOS)定义为植被指数超过或低于阈值的第一天；(b)最大斜率法(DER)。f(t)'是f(t)的一阶导数，SOS(EOS)是生长(衰老)速度最快的一天，峰值(POS)是f(t)的最大值的位置；(c)曲率曲线法(Inflexion)。K'是曲率的变化率；曲率变化率K'最大的两个局部极大值定义为Greenup(返青期)和Maturity(成熟期)，K'最小的两个局部极小值定义为Senescence(衰老期)和Dormancy(休眠期)；(d)Gu法(Gu)：baseline(基线)和maxline(最大线)分别对应f(t)的上下边界，f(t)'的最大值和最小值定义了两条切线(恢复线和衰老线)到f(t)的斜率，恢复线和衰老线与baseline、maxline的交点定义了4个物候指标：上升期(UD)、稳定期(SD)、下降期(DD₀)以及衰退日期(RD)。为了考虑夏季植被指数的微弱下降，Gu法采用SD和DD₀之间的植被指数数据拟合出一条平坦直线，该直线与衰老线的交点定义了最终调整的下降期(DD)。

图9给出了在CA-NS6站点2010-2016年每个生长周期中，使用本发明提供的四种物候提取方法提取的物候指标。其中，AG精细拟合函数重建的日尺度EVI时间序列结果。f(t)为精细拟合平滑的EVI日尺度时间序列(精细拟合列)；f(t)’为f(t)的一阶导数(DER列，长虚线)；K'为曲率变化率(拐点列，虚线)。第一列是精细拟合的结果。其余4列分别对应TRS法、DER法、Inflexion法、Gu法提取的物候指标。其中TRS5为TRS法的缩写，对应的阈值为0.5。图9左列为原EVI和重构EVI，精细拟合指标为确定相关性(R²)、Nash-Sutcliffe模型效率系数(NSE)、均方根误差(RMSE)。右边四列为提取的植被物候指标，分别对应阈值法(TRS)、最大斜率法(DER)、曲率曲线法(Inflexion)和Gu法(Gu)四种植被物候提取方法。

表4本发明与TIMESAT、phenopix之间的比较结果

图10-12进一步展示了本发明在区域尺度上提取物候指标的可靠性。选择的研究区域为河南省，研究时段为2015年，输入的植被指数数据为MOD13A1.006步长为16天的EVI植被指数序列。

图10为第1个生长周期的空间分布。左侧为各方法生长周期开始的时间，右侧为生长周期结束的时间。西部地区生长周期开始时间集中在4月，结束时间在6月；东部地区开始时间在1至2月，结束时间在5月。从图中可以看出，各方法提取的西部、南部与北部地区第一个生长周期开始的时间较为吻合，大约在4月，其中阈值法(TRS)与最大斜率法(DER)的提取结果最为相似；图11为第2个生长周期的空间分布。第2个生长周期大约在7至9月。其中，Gu法(Gu)与曲率曲线法(Inflexion)提取的东部的生长周期时间略晚与阈值法(TRS)与最大斜率法(DER)提取的结果，东部的生长周期多为夏季玉米的生长周期，除此之外，四种方法均未从河南省西部地区提取出植被物候指标，可能原因是河南省西部多为森林、山脉，并不存在具有此生长周期的作物；图12为第3个生长周期的空间分布。第三个植被生长周期仅出现在河南省北部及东部部分地区，时间为11月至次年1月。

总体而言，本发明在对河南省2015年的植被物候指标提取中展现了良好的效果，各方法提取的物候指标较为吻合，体现了本发明的可靠性。

在一些实施例中，还提供了一种基于遥感植被指数提取植被物候指标的装置，包括以下模块：

获取模块，用于获取输入数据，包括：植被指数时间序列、日期信息和质量控制波段；

权重初始化模块，用于根据所述质量控制波段，赋予所述植被指数时间序列中不同质量的数据点不同的权重；

预处理模块，用于识别所述遥感植被指数时间序列的钉值，将钉值看作缺失值，对缺失值采用相邻日期线性插值进行填补，得到预处理后的植被指数时间序列；

粗略拟合模块，用于对所述预处理后的植被指数时间序列进行粗略拟合，得到粗略拟合后的遥感植被指数时间序列；

权重调整模块，用于对所述粗略拟合后的植被指数时间序列中的每个数据点调整权重；

循环模块，用于重复运行粗略拟合模块和权重调整模块以调整权重，直至最大循环次数，输出粗略拟合除噪后的植被指数序列；

生长周期划分模块，用于根据输出的粗略拟合除噪后的植被指数序列，以两个谷值确定一个生长周期为原则，划分生长周期；

精细拟合模块，用于在每个生长周期内，对粗略拟合除噪后的植被指数序列进行精细拟合，得到精细拟合后的遥感植被指数时间序列；

所述权重调整模块，还用于对精细拟合后的遥感植被指数时间序列中的每个数据点调整权重；

所述循环模块，还用于重复运行精细拟合模块和权重调整模块以调整权重，直至最大循环次数，输出精细拟合除噪后的植被指数序列；

物候指标提取模块，用于根据输出的精细拟合除噪后的植被指数序列，提取物候指标。

本发明实施例提供的一种基于遥感植被指数提取植被物候指标的方法及装置从以下几个方面改进了植被物候指标提取精度：(1)采用权重调整函数使植被指数摆脱了多种光学噪点(如云、雪、气溶胶、阴影等)的干扰；(2)采用粗略拟合与精细拟合相结合，使植被物候指标提取在多生长季或植被类型发生变化的区域表现优异。粗略拟合用于捕捉植被周期性信号，并划分生长周期；精细拟合用于在每个生长周期中捕捉植被生长周期关键点的迅速变化；(3)融合了多种权重调整、粗略拟合、精细拟合以及植被物候提取方法，对提取的植被物候指标的不确定性进行了量化。本发明可应用于森林、草地、农田等生态系统生长周期关键性转折点的动态监测，指导生产实践。

需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者系统不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者系统所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括该要素的过程、方法、物品或者系统中还存在另外的相同要素。

上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。在列举了若干装置的单元权利要求中，这些装置中的若干个可以是通过同一个硬件项来具体体现。词语第一、第二、以及第三等的使用不表示任何顺序，可将这些词语解释为标识。

以上仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。

Claims (10)

1.一种基于遥感植被指数提取植被物候指标的方法，其特征在于，包括以下步骤：

步骤1：获取输入数据，包括：植被指数时间序列、日期信息和质量控制波段；

步骤2：根据所述质量控制波段，赋予所述植被指数时间序列中不同质量的数据点不同的权重；

步骤3：识别所述遥感植被指数时间序列的钉值，将钉值看作缺失值，对缺失值采用相邻日期线性插值进行填补，得到预处理后的植被指数时间序列；

步骤4：对所述预处理后的植被指数时间序列进行粗略拟合，得到粗略拟合后的遥感植被指数时间序列；

步骤5：对所述粗略拟合后的植被指数时间序列中的每个数据点调整权重；

步骤6：重复运行步骤4粗略拟合和步骤5调整权重，直至最大循环次数，输出粗略拟合除噪后的植被指数序列；

步骤7：根据步骤6输出的粗略拟合除噪后的植被指数序列，以两个谷值确定一个生长周期为原则，划分生长周期；

步骤8：在每个生长周期内，对粗略拟合除噪后的植被指数序列进行精细拟合，得到精细拟合后的遥感植被指数时间序列；

步骤9：对精细拟合后的遥感植被指数时间序列中的每个数据点调整权重；

步骤10：重复运行步骤8精细拟合和步骤9调整权重，直至最大循环次数，输出精细拟合除噪后的植被指数序列；

步骤11：根据步骤10输出的精细拟合除噪后的植被指数序列，提取物候指标。

2.如权利要求1所述的基于遥感植被指数提取植被物候指标的方法，其特征在于，步骤2中，遥感植被指数时间序列根据质量控制波段划分为三类，良好、一般和较差；这三类数据的权重分别赋值为最大权重、中等权重和最低权重。

3.如权利要求1所述的基于遥感植被指数提取植被物候指标的方法，其特征在于，步骤3具体包括：

步骤3.1：若一个数据点的y_i与y_mov,i的绝对差大于整个时间序列标准差的2倍，那么这个点被标记为一个钉值，将所述钉值作为缺失值；其中y_i是输入植被指数时间序列的第i个观测值，y_mov,i是相应的3点滑动平均值；

步骤3.2：利用相邻日期线性插值对缺失值进行填补；

步骤3.3：确定填补后的植被指数时间序列的有效值域。

4.如权利要求1所述的基于遥感植被指数提取植被物候指标的方法，其特征在于，步骤4中，采用的粗略拟合函数包括：加权Savitzky-Golay、加权HANTs和加权Whittaker。

5.如权利要求1所述的基于遥感植被指数提取植被物候指标的方法，其特征在于，步骤7具体包括：

步骤7.1：根据粗略拟合除噪后的植被指数序列，检测局部最大值或最小值，并将其作为潜在峰值或谷值；

步骤7.2：判断所述潜在峰值或谷值中是否存在由于局部细微波动引起的虚假峰值或谷值，若存在，剔除该虚假峰值或谷值；

步骤7.3：通过两个波谷定义了一个生长周期，生长周期的编号标记为year_k，其中year为生长高峰所在的年份，k为这一年生长季节的顺序。

6.如权利要求5所述的基于遥感植被指数提取植被物候指标的方法，其特征在于，步骤7.2具体包括：

步骤7.2.1：若max(h_peak,L,h_peak,R)<r_maxA或min(h_peak,L,h_peak,R)<r_minA，则认为是虚假值，并进行剔除，h_peak,L、h_peak,R的表达式为：

其中，y_peak为生长周期的峰值；y_trough,L、y_trough,R为生长周期峰值左侧与右侧的谷值；h_peak,L、h_peak,R分别为生长周期峰值到左右谷值的高度差；A为输入植被指数时间序列y的幅值；r_max和r_min是两个待定参数；

步骤7.2.2：若max(h_trough,L,h_trough,R)<r_maxA，亦认为是虚假谷值，并进行剔除，h_trough,L和h_trough,R的表达式为：

其中，y_trough为生长周期的谷值；y_peak,L、y_peak,R为生长周期谷值左侧与右侧的峰值；h_trough,L和h_trough,R分别为生长周期谷值到左右峰值的高度差；

步骤7.2.3：若相邻两个峰或谷的时间差小于len_min，则认为是虚假峰值或谷值，len_min为设定的最小时间差。

7.如权利要求1所述的基于遥感植被指数提取植被物候指标的方法，其特征在于，步骤8中，采用的精细拟合函数包括：非对称高斯分布函数AG和多种双侧逻辑斯蒂函数；

多种双侧逻辑斯蒂函数包括：Zhang、Beck、Elmore和Gu。

8.如权利要求1所述的基于遥感植被指数提取植被物候指标的方法，其特征在于，步骤8还包括：

根据精细拟合函数中各参数的生态意义确定各个参数的有效值域，在参数求解过程中进行约束；

在每个生长周期精细拟合过程中，将该生长季节前后的质量前n的观测加入进来，使不同生长周期之间过渡平滑。

9.如权利要求1所述的基于遥感植被指数提取植被物候指标的方法，其特征在于，步骤11中，采用的植被物候指标提取方法包括：阈值法、最大斜率法、曲率曲线法和Gu法。

10.一种基于遥感植被指数提取植被物候指标的装置，其特征在于，包括以下模块：

获取模块，用于获取输入数据，包括：植被指数时间序列、日期信息和质量控制波段；

权重初始化模块，用于根据所述质量控制波段，赋予所述植被指数时间序列中不同质量的数据点不同的权重；

预处理模块，用于识别所述遥感植被指数时间序列的钉值，将钉值看作缺失值，对缺失值采用相邻日期线性插值进行填补，得到预处理后的植被指数时间序列；

粗略拟合模块，用于对所述预处理后的植被指数时间序列进行粗略拟合，得到粗略拟合后的遥感植被指数时间序列；

权重调整模块，用于对所述粗略拟合后的植被指数时间序列中的每个数据点调整权重；

循环模块，用于重复运行粗略拟合模块和权重调整模块以调整权重，直至最大循环次数，输出粗略拟合除噪后的植被指数序列；

生长周期划分模块，用于根据输出的粗略拟合除噪后的植被指数序列，以两个谷值确定一个生长周期为原则，划分生长周期；

精细拟合模块，用于在每个生长周期内，对粗略拟合除噪后的植被指数序列进行精细拟合，得到精细拟合后的遥感植被指数时间序列；

所述权重调整模块，还用于对精细拟合后的遥感植被指数时间序列中的每个数据点调整权重；

所述循环模块，还用于重复运行精细拟合模块和权重调整模块以调整权重，直至最大循环次数，输出精细拟合除噪后的植被指数序列；

物候指标提取模块，用于根据输出的精细拟合除噪后的植被指数序列，提取物候指标。

Images