CN106803059A - 一种遥感植被指数时间序列森林监测方法 - Google Patents

Abstract

一种遥感植被指数时间序列森林监测方法，包括：获取landsat8遥感植被指数(NDVI)时序数据，并生成研究区域以日为时间步长的跨度数年的时序数据集；基于landsat8内遥感数据的质量文件，利用云雾质量和时序差分法进行噪声数据的识别和剔除；建立每个有效点的基于高斯方程的初始遥感植被指数时序拟合图；针对时序拟合图上各个取景日对应的原始值与曲线拟合值的差值设置权重，迭代计算高斯模型，直到残差达到目标阈值或迭代次数到达上限，得到最终的植被指数时序拟合图；从植被指数时序拟合曲线提取森林年度监测指标。

Description

一种遥感植被指数时间序列森林监测方法

技术领域

本发明属于遥感影像信息处理技术领域，涉及一种遥感植被指数时间序列森林监测方法。

背景技术

遥感卫星拍摄的某一时刻的地表信息被记录为遥感影像，同一地区的多个不同时刻的影像便是遥感时序影像。遥感时序数据具有非常重要的研究价值，例如它可以用来研究植被与物候的关系、植被的变化以及监测森林生长情况等。然而受到光照强度、大气湿度、云层厚度以及天气等多重因素影响，多源噪音干扰导致遥感数据质量下降，难以准确地提取地表信息，对遥感的研究造成了很大的困难。

归一化差分植被指数(NDVI)能够部分消除辐射误差，被广泛应用于监测植被生长状态、植被覆盖度，是一种最常用的描述地表特征的定量指标。NDVI时序也是研究遥感植被时序的常用数据。

遥感数据时序拟合法是目前常用的一种遥感影像降噪方法。经典的遥感时序方法有：Savitzky-Golay(S-G)滤波法、傅里叶滤波法和高斯滤波法。其中S-G滤波法能保持信号的形状宽度不变，较好地体现植被的物候规律，但其参数选择受制于人的主观性。傅里叶滤波法将时间曲线表达为一系列余(正)弦波的线性叠加，通过筛选重要波段信息而使拟合结果真实可靠，但该方法对于有效时序较稀疏的情形适应性不够好。高斯滤波法非线性的高斯函数以最小二乘法拟合，克服了主观性，在长时间范围内的信息提取中优势明显，但该方法对原始数据点的噪声较敏感。

本方法基于先验知识对原始多年度遥感时间序列植被指数进行有效预处理及高斯迭代拟合优化，并以提取的生长参数监测森林。

发明内容

本发明所要解决的问题是：提出一种遥感植被指数时间序列森林监测方法，能够有效过滤遥感时序数据中的大量噪声，使多年度植被指数序列的拟合结果值更接近真实，用于监测森林的多年度变化。

本发明提供的技术方案是：一种遥感植被指数时间序列森林监测方法，其特征在于包括以下步骤：

1.获取landsat8遥感植被指数(NDVI)时序数据，并生成研究区域以日为时间步长的跨度数年的时序数据集；

2.基于landsat8内遥感数据的质量文件，利用云雾质量和时序差分法进行噪声数据的识别和剔除；

3.建立每个有效点的基于高斯方程的初始遥感植被指数时序拟合图；

4.针对时序拟合图上各个取景日对应的原始值与曲线拟合值的差值设置权重，迭代计算高斯模型，直到残差达到目标阈值或迭代次数到达上限，得到最终的植被指数时序拟合图；

5.从植被指数时序拟合曲线提取森林年度监测指标。

步骤1的实现方式是：根据L景landsat8遥感影像数据文件合成的归一化差分植被指数(NDVI)时序数据，生成研究区域内每个植被像元以日为步长的NDVI多年度时序数据集，其中L为大于1的自然数；

遥感数据采集的等间隔时间序列T，统一以年内的第几天计，跨年度时分别按各自年度内计天数；可记为向量T＝＜t₁,t₂,…,t_L＞，其中t_i为第i个遥感数据采集的时刻；同时同一地区的多景遥感影像构成的遥感影像序列可相应表示为其中为t_i时刻采集到的一景遥感影像；

NDVI值可由遥感原始数据的近红外波段NIR和红波段R根据公式1计算得到。

步骤2的实现方式是：对于每景landsat8影像可表述为R行C列的矩阵，若植被像元位于矩阵第r行c列，则可从影像序列I(T)中提取得到该植被像元的NDVI时序，表示为其中第i个分量表示位于第r行c列的植被像元在时间序列T的t_i时刻采集到的数据，根据公式1计算得到的NDVI值；针对所有的植被像元，根据其质量波段的高四位，将表现为1*1*的二进制数值对应的植被像元及其时序数据剔除；由于遥感所处的复杂噪声环境，初步去噪后的遥感观测时序向量仍需要进一步提升精度；因而还需要分析剩余数据点，根据差分法剔除时序数据中的噪声点；具体为：针对每个像元，遍历其时序数据，计算前后相邻点的植被指数绝对差并除以两点对应的间隔天数；如果连续两次的结果都大于阈值，则判断该时刻点的数据为噪声值应剔除；

遍历整个向量，判断每个NDVI时序分量是否为噪声：

then flag＝true

else flag＝false

其中，θ为阈值，序号i∈{2,...,L-1}，L为大于1的自然数，t_i是时间序列T的第i个分量，表示位于第r行c列的植被像元在t_i时刻的NDVI值，flag是噪声判别变量，flag值为真(true)时表示分量为噪声应去除，值为假(false)时表示分量应保留；

剔除噪声后，得到新的时间向量T'＝＜t'₁,t'₂,…,t'_m＞，其中m≤L且相应每个像元原始的NDVI时序长度变短，对于每个植被像元来说，m大小不一，普遍不足原序列长度的一半，必定小于L；第r行c列的植被像元的NDVI时序去噪后变更为

步骤3的实现方式是：将初步去噪后的遥感观测时序向量采用高斯模型，以最小二乘法求解得到初始的连续拟合曲线f₁(x)；

其中e为自然常数，x表示时间x∈[0,t_L]，t_L为T的最大分量，T'的每个分量为x的采样点；n为年数，对应拟合曲线f₁的高斯函数波峰个数；j为某个年度，参数a_1,j，b_1,j和c_1,j对应每个峰的振幅、中心点和半宽；g₁是求解方程的目标函数，即在采样点(T'的每个分量)处曲线f₁与去噪后的NDVI时序分量的拟合误差的差方最小；

当前的残差向量可通过计算曲线在T'的各有效观测时点上时序上的原始观测值与高斯函数初次拟合值的差得到，即

步骤4的实现方式是：基于有效观测时刻的拟合残差，并结合概率密度曲线得到修正的权重系数，继续迭代计算高斯模型，直到残差达到目标阈值或迭代次数到达上限，得到最终的植被指数时序拟合图；

由于遥感所处的复杂噪声环境，初次高斯拟合后的时序植被植被指数误差一般高于目标阈值；采用多次高斯拟合迭代修正的方式，使遥感观测时序向量逐步提升精度得到高精度的最终拟合结果其过程可以表示为序列f₁(T')，f₂(T')，…，其中下标0,1,2，…，k等表示拟合次数，f_i(T')表示第i次的拟合曲线在T'各分量时刻采样得到的向量，第i次的拟合f_i转向下一次拟合f_i+1借助拟合修正向量Y_i和Y_i+1完成；

初次拟合后，如残差D₁高于阈值，需要进行后续拟合；设计迭代调整的权重向量W₁如公式4所示，其中u₁、σ₁为向量的均值和方差；

令第r行c列的植被像元的初始拟合修正向量为下一次拟合的修正向量Y₁ ^(r，c)的修正规则如下：

其中i表示向量的第i个分量，W₁为权重向量，D₁为残差，T'为去噪后的时间序列；

设第k次迭代，相应的高斯函数拟合曲线f_k、目标函数g_k、权重向量W_k和残差D_k计算相比于初次迭代变更如下：

第k次迭代时拟合修正向量和的计算规则为：

其中i表示向量的第i个分量，W_k为权重向量，D_k为残差，T'为去噪后的时间序列；

直到残差小于阈值或达到迭代上限，则终止迭代计算，得到最终优化的植被指数拟合曲线，否则继续本步骤迭代。

步骤5的实现方式是：以NDVI年度生长期、快速生长起始日期、生长幅度作为森林监测指标；根据最终的植被指数时序拟合曲线，先保存末次迭代得出的高斯函数的参数，并从曲线提取出年度区间内的植被指数最大斜率点、最小值、最大值及对应时间，进一步计算得到年度生长期、快速生长起始日期、生长幅度；

森林一般为多年生植物，每年稳定有一个生长峰，高斯函数的参数表征了该峰的振幅、中心点位置和方差，可以保存作为长期动态监测参考；设末次迭代为第p次迭代，得到的高斯函数如公式9，其中e为自然常数，x表示时间x∈[0,t_L]，t_L为T的最大分量，T'的每个分量为x的采样点；n为年数，对应拟合曲线f_p的高斯函数波峰个数；j为某个年度，参数a_p,1，b_p,1，c_p,1，…，a_p,n，b_p,n，c_p,n对应每个峰的振幅、中心点和半宽，最终保存参数a_p,1，b_p,1，c_p,1，…，a_p,n，b_p,n，c_p,n；

通过森林像元在自然年度的生长期、快速生长起始日期、生长幅度这三个指标来监测森林的生长动态；定义年度生长期为年度内最大值点日期减去最小值点日期，定义年度快速生长起始日期为年度内曲线斜率最大点对应的日期，定义年度NDVI生长幅度为年度内最大值减去最小值；在各自然年度区间内，提取得到曲线的最小值、最大值、斜率最大点对应的年度内日期，可计算得到年度生长期、快速生长起始日期、生长幅度。

本发明的优点是：(1)具有噪声鲁棒性，可以较好地适应遥感有效数据稀疏的情况。(2)利用NDVI被噪声低估的先验知识，通过数据质量评价评估权重，迭代拟合获取稳定的拟合结果。(3)人工主观影响小，可以实现多年度区域级的参数自动提取，自动监测森林生长状况。

附图说明

图1是本发明方法的流程图。

具体实施方式

下面结合附图及实施例对本发明做进一步说明。

本发明提供一种遥感植被指数时间序列森林监测方法，图1是具体流程。以下是本发明一实施例及其步骤：

步骤1：获取landsat8遥感植被指数(NDVI)时序数据，并生成研究区域以日为时间步长的跨度数年的时序数据集；

在USGS(http://earthexplorer.usgs.gov)上获取地区条带号为119、039，空间分辨率为30米，包含浙江省杭州2014和2015两年的等间隔16天的共45景的LandSat OLI遥感数据。根据45景landsat8遥感影像数据文件合成的归一化差分植被指数(NDVI)时序数据，生成研究区域内每个植被像元以日为步长的NDVI多年度时序数据集；

遥感数据采集的等间隔时间序列T，统一以年内的第几天计，跨年度时分别按各自年度内计天数；可记为向量T＝＜t₁,t₂,…,t₄₅＞，其中t_i为第i个遥感数据采集的时刻；同时同一地区的多景遥感影像构成的遥感影像序列可相应表示为其中为t_i时刻采集到的一景遥感影像；

NDVI值可由遥感原始数据的近红外波段NIR和红波段R根据公式1计算得到。

步骤2：基于landsat8内遥感数据的质量文件，利用云雾质量和时序差分法进行噪声数据的识别和剔除；

对于每景landsat8影像可表述为R行C列的矩阵，选取位于(30^°10＇11＂N,119^°57＇14.6＂E)的森林像元，则可从影像序列I(T)中提取得到该植被像元的NDVI时序，表示为其中第i个分量表示该森林像元在时间序列T的t_i时刻采集到的数据，根据公式1计算得到的NDVI值；根据其质量波段的高四位，将表现为1*1*的二进制数值对应的像元及其时序数据剔除；由于遥感所处的复杂噪声环境，初步去噪后的遥感观测时序向量仍需要进一步提升精度；因而还需要分析剩余数据点，根据差分法剔除时序数据中的噪声点；具体为：针对每个像元，遍历其时序数据，计算前后相邻点的植被指数绝对差并除以两点对应的间隔天数；如果连续两次的结果都大于阈值，则判断该时刻点的数据为噪声值应剔除；

遍历整个向量，判断每个NDVI时序分量是否为噪声：

then flag＝true

else flag＝false

其中，0.017为阈值，序号i∈{2,...,45-1}，t_i是时间序列T的第i个分量，表示森林像元在t_i时刻的NDVI值，flag是噪声判别变量，flag值为真(true)时表示分量为噪声应去除，值为假(false)时表示分量应保留；

剔除噪声数据后，本例剩余15个有效数据，得到新的时间向量T'＝＜t'₁,t'₂,…,t'₁₅＞，相应每个像元原始的NDVI时序长度变短，变更为

步骤3：将初步去噪后的遥感观测时序向量N_T'采用高斯模型，以最小二乘法求解得到初始的连续拟合曲线f₁(x)；

其中e为自然常数，x表示时间x∈[0,t₄₅]，t₄₅为T的最大分量，T'的每个分量为x的采样点；本例中拟合曲线f₁的高斯函数波峰为2个；参数a_1,1,a_1,2对应两个峰的振幅，b_1,1,b_1,2对应两个峰的中心点，c_1,1,c_1,2对应两个峰半宽；g₁是求解方程的目标函数，即在采样点(T'的每个分量)处曲线f₁与去噪后的NDVI时序分量N_T'的拟合误差的差方最小；

当前的残差向量可通过计算曲线在T'的各有效观测时点上N_T'时序上的原始观测值与高斯函数初次拟合值的差得到，即D₁＝N_T'-f₁(T')，本例此处值为0.12。

步骤4：基于有效观测时刻的拟合残差，并结合概率密度曲线得到修正的权重系数，继续迭代计算高斯模型，直到残差达到目标阈值或迭代次数到达上限，得到最终的植被指数时序拟合图；

由于遥感所处的复杂噪声环境，初次高斯拟合后的时序植被植被指数误差一般高于目标阈值；采用多次高斯拟合迭代修正的方式，使遥感观测时序向量N_T'逐步提升精度得到高精度的最终拟合结果其过程可以表示为序列N_T'(＝f₀(T'))，f₁(T')，f₂(T')，…，其中下标0,1,2，…，k等表示拟合次数，f_i(T')表示第i次的拟合曲线在T'各分量时刻采样得到的向量，第i次的拟合f_i转向下一次拟合f_i+1借助拟合修正向量Y_i和Y_i+1完成；

初次拟合后，如残差D₁高于阈值，需要进行后续拟合；本例迭代次数上限设为50，最大误差分量绝对值阈值为0.01。本例初次迭代后的残差最大分量绝对值达到0.12，超过了阈值，因此需要继续进行下一次迭代；

设计迭代调整的权重向量W₁如公式4所示，其中u₁、σ₁为向量N_T'的均值和方差；

令森林像元的初始拟合修正向量为Y₀＝f₁(T')，下一次拟合的向量Y₁的修正规则如下：

If(Y_0,i＜f₁(T'_i))

then Y_1,i＝f₁(T'_i)-W₁(Y₀)×D_1,i

else Y_1,i＝f₁(T'_i)

其中i表示向量的第i个分量，W₁为权重向量，D₁为残差，T'为去噪后的时间序列；本例后续以Y₁作为新的高斯函数拟合采样值进行迭代；

设第k次迭代，相应的高斯函数拟合曲线f_k、目标函数g_k、权重向量W_k和残差D_k计算相比于初次迭代变更如下：

D_k＝N_T'-f_k(T') (8)

第k次迭代时拟合修正向量Y_k-1和Y_k的计算规则为：

If(Y_k-1,i＜f_k(T'_i))

then Y_k,i＝f_k(T'_i)-W_k(Y_k-i,i)×D_k,i

else Y_k,i＝f_k(T'_i)

其中i表示向量的第i个分量，W_k为权重向量，D_k为残差，T'为去噪后的时间序列；

直到残差小于阈值或达到迭代上限，则终止迭代计算，得到最终优化的植被指数拟合曲线，否则继续本步骤迭代；

本例在第14次迭代时，残差为0.008小于阈值0.01，终止迭代，得到该森林像元的植被指数拟合曲线。

步骤5：以NDVI年度生长期、快速生长起始日期、生长幅度作为森林监测指标。根据最终的植被指数时序拟合曲线，先保存末次迭代得出的高斯函数的参数，并从曲线提取出年度区间内的植被指数最大斜率点、最小值、最大值及对应时间，进一步计算得到年度生长期、快速生长起始日期、生长幅度；

森林一般为多年生植物，每年稳定有一个生长峰，高斯函数的参数表征了该峰的振幅、中心点位置和方差，可以保存作为长期动态监测参考；本例末次迭代为第14次迭代，得到的高斯函数如公式9，其中e为自然常数，x表示时间x∈[0,t₄₅]，t₄₅为T的最大分量，T'的每个分量为x的采样点；年数为2，对应拟合曲线f₁₄的高斯函数波峰个数；j为某个年度，参数a_14,1，b_14,1，c_14,1，a_14,2，b_14,2，c_14,2对应每个峰的振幅、中心点和半宽，最终保存参数a_14,1，b_14,1，c_14,1，a_14,2，b_14,2，c_14,2；本例森林像元分别得到上述参数值为：0.43137，10.02，8.5184，0.44454，33.297，12.852；

通过森林像元在自然年度的生长期、快速生长起始日期、生长幅度这三个指标来监测森林的生长动态；定义年度生长期为年度内最大值点日期减去最小值点日期，定义年度快速生长起始日期为年度内曲线斜率最大点对应的日期，定义年度NDVI生长幅度为年度内最大值减去最小值；在各自然年度区间内，提取得到曲线的最小值、最大值、斜率最大点对应的年度内日期，可计算得到年度生长期、快速生长起始日期、生长幅度；

本实施例中得到第一年度(2014年)生长期为167天、快速生长起始日期为3月6日、生长幅度为0.34045，第二年度(2015年)内生长期为188天、快速生长起始日期为3月4日、生长幅度为0.25081。

本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的保护范围不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也及于本领域技术人员根据本发明构思所能够想到的等同技术手段。

Claims (7)

1.一种遥感植被指数时间序列森林监测方法，其特征在于如下处理流程：

S1.获取landsat8遥感植被指数(NDVI)时序数据，并生成研究区域以日为时间步长的跨度数年的时序数据集；

S2.基于landsat8内遥感数据的质量文件，利用云雾质量和时序差分法进行噪声数据的识别和剔除；

S3.建立每个有效点的基于高斯方程的初始遥感植被指数时序拟合图；

S4.针对时序拟合图上各个取景日对应的原始值与曲线拟合值的差值设置权重，迭代计算高斯模型，直到残差达到目标阈值或迭代次数到达上限，得到最终的植被指数时序拟合图；

S5.从植被指数时序拟合曲线提取森林年度监测指标。

2.根据权利要求1所述的遥感植被指数时间序列森林监测方法，其特征在于：所述步骤S1的实现方式是：根据L景landsat8遥感影像数据文件合成的归一化差分植被指数(NDVI)时序数据，生成研究区域内每个植被像元以日为步长的NDVI多年度时序数据集，其中L为大于1的自然数；

遥感数据采集的等间隔时间序列T，统一以年内的第几天计，跨年度时分别按各自年度内计天数；可记为向量T＝＜t₁,t₂,…,t_L＞，其中t_i为第i个遥感数据采集的时刻；同时同一地区的多景遥感影像构成的遥感影像序列可相应表示为其中为t_i时刻采集到的一景遥感影像；

NDVI值可由遥感原始数据的近红外波段NIR和红波段R根据公式1计算得到。

$NDVI=\frac{NIR-R}{NIR+R}---\left(1\right)$

3.根据权利要求1所述的遥感植被指数时间序列森林监测方法，其特征在于：所述步骤S2的实现方式是：对于每景landsat8影像可表述为R行C列的矩阵，若植被像元位于矩阵第r行c列，则可从影像序列I(T)中提取得到该植被像元的NDVI时序，表示为其中第i个分量表示位于第r行c列的植被像元在时间序列T的t_i时刻采集到的数据，根据公式1计算得到的NDVI值；针对所有的植被像元，根据其质量波段的高四位，将表现为1*1*的二进制数值对应的植被像元及其时序数据剔除；由于遥感所处的复杂噪声环境，初步去噪后的遥感观测时序向量仍需要进一步提升精度；因而还需要分析剩余数据点，根据差分法剔除时序数据中的噪声点；具体为：针对每个像元，遍历其时序数据，计算前后相邻点的植被指数绝对差并除以两点对应的间隔天数；如果连续两次的结果都大于阈值，则判断该时刻点的数据为噪声值应剔除；

遍历整个向量，判断每个NDVI时序分量是否为噪声：

$If\left(\begin{array}{ccc}\frac{N_{t_i}^{(r,c)}-N_{t_{i-1}}^{(r,c)}}{t_i-t_{i-1}}>\mathrm{\θ}& and& \frac{N_{t_{i+1}}^{(r,c)}-N_{t_i}^{(r,c)}}{t_{i+1}-t_i}>\mathrm{\θ}\end{array}\right)$

then flag＝true

else flag＝false

其中，θ为阈值，序号i∈{2,...,L-1}，L为大于1的自然数，t_i是时间序列T的第i个分量，表示位于第r行c列的植被像元在t_i时刻的NDVI值，flag是噪声判别变量，flag值为真(true)时表示分量为噪声应去除，值为假(false)时表示分量应保留；

剔除噪声后，得到新的时间向量T'＝＜t'₁,t'₂,…,t'_m＞，其中m≤L且相应每个像元原始的NDVI时序长度变短，对于每个植被像元来说，m大小不一，普遍不足原序列长度的一半，必定小于L；第r行c列的植被像元的NDVI时序去噪后变更为

4.根据权利要求1所述的遥感植被指数时间序列森林监测方法，其特征在于：所述步骤S3的实现方式是：将初步去噪后的遥感观测时序向量采用高斯模型，以最小二乘法求解得到初始的连续拟合曲线f₁(x)；

$f_1\left(x\right)=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{a}_{1,j}\×e^{-{\left(\frac{x-b_{1,j}}{c_{1,j}}\right)}^2}---\left(2\right)$

其中e为自然常数，x表示时间x∈[0,t_L]，t_L为T的最大分量，T'的每个分量为x的采样点；n为年数，对应拟合曲线f₁的高斯函数波峰个数；j为某个年度，参数a_1,j，b_1,j和c_1,j对应每个峰的振幅、中心点和半宽；g₁是求解方程的目标函数，即在采样点(T'的每个分量)处曲线f₁与去噪后的NDVI时序分量的拟合误差的差方最小；

$g_1\left(T^{\′}\right)=\arg \underset{i=1}{\overset{m}{\min }}{(N_{{t^{\′}}_i}^{(r,c)}-f_1({T^{\′}}_i\left)\right)}^2---\left(3\right)$

当前的残差向量可通过计算曲线在T'的各有效观测时点上时序上的原始观测值与高斯函数初次拟合值的差得到，即

5.根据权利要求1所述的遥感植被指数时间序列森林监测方法，其特征在于：所述步骤S4的实现方式是：基于有效观测时刻的拟合残差，并结合概率密度曲线得到修正的权重系数，继续迭代计算高斯模型，直到残差达到目标阈值或迭代次数到达上限，得到最终的植被指数时序拟合图；

由于遥感所处的复杂噪声环境，初次高斯拟合后的时序植被植被指数误差一般高于目标阈值；采用多次高斯拟合迭代修正的方式，使遥感观测时序向量逐步提升精度得到高精度的最终拟合结果其过程可以表示为序列其中下标0,1,2，…，k等表示拟合次数，f_i(T')表示第i次的拟合曲线在T'各分量时刻采样得到的向量，第i次的拟合f_i转向下一次拟合f_i+1借助拟合修正向量Y_i和Y_i+1完成；

初次拟合后，如残差D₁高于阈值，需要进行后续拟合；设计迭代调整的权重向量W₁如公式4所示，其中u₁、σ₁为向量的均值和方差；

$W_1\left(T^{\′}\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_1}{e}^{\frac{{(T^{\′}-u_1)}^2}{2{{\mathrm{\σ}}_1}^2}}---\left(4\right)$

令第r行c列的植被像元的初始拟合修正向量为下一次拟合的修正向量Y₁ ^(r,c)的修正规则如下：

$\begin{array}{cc}If& (Y_{0,i}^{\left(r,c\right)}<f_1({T^{\&prime;}}_i\left)\right)\end{array}$

$\begin{array}{cc}then& Y_{1,i}^{\left(r,c\right)}=f_1\left({T^{\&prime;}}_i\right)-W_1\left(Y_0^{\left(r,c\right)}\right)\&times;D_{1,i}\end{array}$

$\begin{array}{cc}else& Y_{1,i}^{(r,c)}=f_1\left({T^{\&prime;}}_i\right)\end{array}$

其中i表示向量的第i个分量，W₁为权重向量，D₁为残差，T'为去噪后的时间序列；

设第k次迭代，相应的高斯函数拟合曲线f_k、目标函数g_k、权重向量W_k和残差D_k计算相比于初次迭代变更如下：

$f_k\left(x\right)=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{a}_{k,j}\&times;e^{-{\left(\frac{x-b_{k,j}}{c_{k,j}}\right)}^2}---\left(5\right)$

$g_k\left(T^{\&prime;}\right)=\underset{i=1}{\overset{m}{\arg \min }}{(Y_{k-1,i}^{(r,c)}-f_k({T^{\&prime;}}_i\left)\right)}^2---\left(6\right)$

$W_k\left(Y_{k-1}^{(r,c)}\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_k}{e}^{\frac{{(Y_{k-1}^{(r,c)}-u_k)}^2}{2{{\mathrm{\&sigma;}}_k}^2}}---\left(7\right)$

$D_k=N_{T^{\&prime;}}^{\left(r,c\right)}-f_k\left(T^{\&prime;}\right)---\left(8\right)$

第k次迭代时拟合修正向量和的计算规则为：

$\begin{array}{cc}If& (Y_{k-1,i}^{(r,c)}<f_k({T^{\&prime;}}_i\left)\right)\end{array}$

$\begin{array}{cc}then& Y_{k,i}^{(r,c)}=f_k\left({T^{\&prime;}}_i\right)-W_k\left(Y_{k-1,i}^{(r,c)}\right)\&times;D_{k,i}\end{array}$

$\begin{array}{cc}else& Y_{k,i}^{(r,c)}=f_k\left({T^{\&prime;}}_i\right)\end{array}$

其中i表示向量的第i个分量，W_k为权重向量，D_k为残差，T'为去噪后的时间序列；

直到残差小于阈值或达到迭代上限，则终止迭代计算，得到最终优化的植被指数拟合曲线，否则继续本步骤迭代。

6.根据权利要求1所述的遥感植被指数时间序列森林监测方法，其特征在于：所述步骤S5的实现方式是：以NDVI年度生长期、快速生长起始日期、生长幅度作为森林监测指标；根据最终的植被指数时序拟合曲线，先保存末次迭代得出的高斯函数的参数，并从曲线提取出年度区间内的植被指数最大斜率点、最小值、最大值及对应时间，进一步计算得到年度生长期、快速生长起始日期、生长幅度；

森林一般为多年生植物，每年稳定有一个生长峰，高斯函数的参数表征了该峰的振幅、中心点位置和方差，可以保存作为长期动态监测参考；设末次迭代为第p次迭代，得到的高斯函数如公式9，其中e为自然常数，x表示时间x∈[0,t_L]，t_L为T的最大分量，T'的每个分量为x的采样点；n为年数，对应拟合曲线f_p的高斯函数波峰个数；j为某个年度，参数a_p,1，b_p,1，c_p,1，…，a_p,n，b_p,n，c_p,n对应每个峰的振幅、中心点和半宽，最终保存参数a_p,1，b_p,1，c_p,1，…，a_p,n，b_p,_n，c_p,n；

$f_p\left(x\right)=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{a}_{p,j}\&times;e^{-{\left(\frac{x-b_{p,j}}{c_{p,j}}\right)}^2}---\left(9\right)$

通过森林像元在自然年度的生长期、快速生长起始日期、生长幅度这三个指标来监测森林的生长动态；定义年度生长期为年度内最大值点日期减去最小值点日期，定义年度快速生长起始日期为年度内曲线斜率最大点对应的日期，定义年度NDVI生长幅度为年度内最大值减去最小值；在各自然年度区间内，提取得到曲线的最小值、最大值、斜率最大点对应的年度内日期，可计算得到年度生长期、快速生长起始日期、生长幅度。

7.根据权利要求1至6任一项所述的遥感植被指数时间序列森林监测方法，其特征在于：该方法可应用于林业等部门实施区域以上尺度的林区年度动态监测。