[发明专利]一种基于高时空分辨率遥感数据的植被物候提取方法

权利要求书

1.一种基于高时空分辨率遥感数据的植被物候提取方法，其特征在于，包括如下步骤：

植被NDVI生长曲线生成步骤S110：获取目标区域的高时空分辨率的序列遥感影像，对所述遥感影像进行云检测预处理，生成植被NDVI生长曲线；

NDVI生长曲线优化步骤S120：设计噪点去除方法对NDVI生长曲线进行优化，生成稳定、光滑、可靠的植被NDVI生长曲线；

植被物候期特征点提取步骤S130：

利用光滑NDVI曲线NDVI_s、超级光滑NDVI曲线NDVI_ss和NDVI_s曲线变化率SLP_NDVI_s，基于最大变化斜率法提取到植被物候期的三个特征点，对目标区域地表植被进行物候提取，所述植被物候期的三个特征点，分别为植被生长季初期SOS，生长季峰期POS和生长季末期EOS，具体的，年内NDVI时间序列最大值对应的时间点DOY(Day of Year)作为POS,处于POS之前，NDVI值上升最快的时间点被认为是SOS，处于POS之后，NDVI值下降最快的时间点被认为是EOS；

目标像元最优物候期生成步骤S140：去除伪物候点，生成目标像元的最优物候期，并最终完成目标区域的地表植被物候提取；

步骤S110具体包括：

利用FMASK云检测算法，对所述序列遥感影像中的云像元进行检测过滤；利用过滤后像元的NDVI公式，即公式(1)，计算出每个像元的NDVI时间序列，考虑目标年份冬季观测质量较差，将目标年份两端分辨延长三个月观测，形成长时间序列NDVI生长曲线，其中，所述NDVI公式为，

其中，ρ_nir为近红外波段反射率，ρ_red为红光波段反射率；

在步骤S120中，所述设计噪点去除方法对NDVI生长曲线进行优化包括：

获取每20天内的NDVI时间序列，按NDVI值从小到大排序，将有序NDVI序列值的80％作为每20天的NDVI，生成每20天的NDVI时间序列；对所述的每20天的NDVI时间序列进行噪点过滤，同时满足以下两个噪点过滤公式的NDVI值，被认为是云检测算法遗漏的云/云阴影像元，所述噪点过滤公式包括公式(2)、(3)，

Min(NDVI_t-1-NDVI_t,NDVI_t+1-NDVI_t)/Max(NDVI_t-1-NDVI_t,NDVI_t+1-NDVI_t)＞0.45(2)

Max(NDVI_t-1-NDVI_t,NDVI_t+1-NDVI_t)＞NDVI_max*0.3 (3)，

其中，NDVI_t为目标像元目标时刻t的像元NDVI值，NDVI_t-1为前一时刻t-1像元NDVI值，NDVI_t+1为后一个时刻t+1的像元NDVI值，NDVI_max为年NDVI时间序列的最大值；

在步骤S120中，所述生成稳定、光滑、可靠的植被NDVI生长曲线包括：利用线性插值算法，对过滤的NDVI时间序列中的缺值进行插补，得到逐日的NDVI时间序列；利用Savitzky-Golay算法，对插补后的NDVI时间序列采用两个不同时间长度进行滤波，得到两个光滑的NDVI曲线；计算利用短时间长度光滑的NDVI曲线变化率，以目标时间t为中心，基于相邻8天时间窗口，利用如下公式(4)、(5)计算每日NDVI的变化斜率，

ncov_t＝COV((DOY_t-8,...DOY_t-1,DOY_t,DOY_t+1,..DOY_t+8),(NDVI_t-8,...NDVI_t-1,NDVI_t,NDVI_t+1,..NDVI_t+8))(4)

SLP_t＝ncov_t[0,1]/ncov_t[0,0] (5)

其中(DOY_t-8,...DOY_t-1,DOY_t,DOY_t+1,..DOY_t+8)为时间窗口内的DOY序列作为自变量，DOY表示年内日序数，(NDVI_t-8,...NDVI_t-1,NDVI_t,NDVI_t+1,..NDVI_t+8)为时间窗口内的NDVI序列，作为因变量，COV代表计算两个变量的协方差矩阵，ncov_t[0,1]为两个变量的协方差，ncov_t[0,0]为自变量方差，SLP_t代表目标时间点的NDVI变化率；裁剪两端多余时间点，利用较短时间长度滤波的作为光滑NDVI曲线NDVI_s，利用较长时间长度滤波的作为超级光滑NDVI曲线NDVI_ss，SLP_NDVI_s为NDVI_s曲线变化率；

具体的，步骤S130包括如下步骤：

(1)提取目标像元物候期的POS：基于目标像元的NDVI_ss曲线，通过三个相邻时间点NDVI值差值对比，确定NDVI_ss曲线的局部极大值点和局部极小值点，利用NDVI_ss曲线确定局部极值点的目的是为了大致确定SOS、POS、EOS位置，同时由于过度光滑可防止提取到较多的伪物候点位置；将多个局部极大值点依据NDVI_ss曲线值进行降序排列，然后从第一个局部极大值点开始迭代提取多个POS、SOS和EOS；以目标迭代的局部极大值点为中心，依据预先设定的物候期长度，分别向前向后延长半个物候长度，基于NDVI_s曲线，在以局部极大值点为中心的一个物候期内寻找NDVI_s曲线的最大NDVI值，即NDVI_i,如果NDVI_i满足公式(6)，那么对应的DOY作为物候期的一个有效POS，

NDVI_i＞(NDVI_{s_min}+(NDVI_{s_max}-NDVI_{s_min})/2)*NDVI_i/NDVI_{s_max}， (6)

其中，NDVI_{s_min}为NDVI_s曲线的最小值，NDVI_{s_max}为NDVI_s曲线的最大值；

(2)提取目标像元物候期的SOS：基于已提取到的POS位置，如果局部最小值点NDVI_ss(i)满足公式(7)，则被认为是POS之前可靠的局部极小值点，

NDVI_ss(i)/NDVI_ss(POS)＜0.9 (7)，

其中，NDVI_ss(POS)为POS处NDVI_ss曲线对应的NDVI值，将所述局部极小值向前延长一个预先设定的物候期长度，在其对应的时间点和POS之间求取NDVI_s最小值对应的DOY；在所述DOY和POS之间，如果满足公式(8)、公式(9)，则该点被认为是SOS，

SLP_NDVI_s(SOS)＝＝Max(SLP_NDVI_s) (8)

NDVI_s(SOS)/NDVI_s(POS)＜0.95 (9)

其中，SLP_NDVI_s(SOS)代表SOS处的曲线变化率为SLP_NDVI_s曲线的最大值，且该处和POS处的NDVI比值小于0.95；

(3)提取目标像元物候期的EOS：基于已提取到的POS位置，如果局部最小值点NDVI_ss(i)满足公式(7)，则被认为是POS之后可靠的局部极小值点，

将所述局部极小值向后延长一个预先设定的物候期长度，在POS和其延长对应的时间点之间求取NDVI_s最小值对应的DOY；在POS和所述DOY之间，如果满足公式(10)和公式(11)，则该点被认为是EOS，

SLP_NDVI_s(EOS)＝＝Min(SLP_NDVI_s) (10)

NDVI_s(EOS)/NDVI_s(POS)＜0.95 (11)，

其中，SLP_NDVI_s(EOS)代表EOS处的曲线变化率为SLP_NDVI_s的最小值，且该处和POS处的NDVI比值小于0.95；

(4)循环上述步骤(1)-(3)，直至遍历完目标像元NDVI_ss曲线上获取的全部局部极大值，完成目标像元多个物候期提取；

步骤S140具体为，如果只提取到一个物候期的像元，则直接存储该目标像元的最优物候期，并记录物候期的质量信息；如果提取到多个物候期，对于三个物候点均存在且POS离180天较近的物候期直接作为最优物候期，并记录物候期的质量信息；对于连续两个不完整物候期，即第一个物候期的SOS和后一个物候期的EOS均未检测到，那么判定为噪点影响，将所述第一个物候期的SOS作为最终SOS，两个物候期POS处NDVI值最大的作为最终POS，第二个物候期的EOS作为最终EOS，生成目标像元最优物候期，并最终完成目标区域的地表植被物候提取。

2.一种存储介质，用于存储计算机可执行指令，其特征在于：

所述计算机可执行指令在被处理器执行时执行权利要求1所述的基于高时空分辨率遥感数据的植被物候提取方法。