洞庭湖流域植被物候特征参数遥感提取

彭检贵，宋亚斌，宁小斌，江腾宇

(国家林业和草原局中南调查规划设计院, 长沙 410014)

0 引言

植被的物候变化直接反映了周边地区的气候变化，成为植被动态模拟的关键[1-2]。同时，植被物候数据可以为植被动态监测、区域农业规划、区域资源配置与决策、区域气候变化研究等提供重要参考。

植物物候监测一般可分为传统监测方法和遥感监测方法。传统物候监测方法需要耗费大量的人力、物力，而且受时间限制，监测周期短，覆盖区域小[3-4]。而遥感物候监测方法仅需要耗费相对较低的费用而获取大尺度区域的物候参数，能及时、有效地对植物物候进行监测[4-8]。已有大量学者基于NDVI和EVI植被指数的各种遥感数据源的时间序列提出了大量物候特征参数遥感提取方法，如White等[3]将传统物候观测模型与遥感物候参数提取结合，针对特定生物群落生态系统，提出NDVI比率阈值物候参数提取方法；Reed等[4]基于NDVI时间序列数据，使用滑动平均方法计算提取了美国4种代表性植物的(荒漠灌丛、森林、草地、农作物)物候特征；Yu等[9]提出了最大变化斜率方法，通过限制不同阈值的最大变化斜率范围，计算了中亚东部区域植物生长季开始期；Zhang等[10]通过对每年NDVI时间序列进行分段Logistic函数拟合，利用拟合后的曲线的曲率变化的特点，得出植物的生长季开始期、结束期、生长季长度以及生长季峰值等。已有的研究表明，基于遥感数据的物候特征提取方法主要有经验公式法、NDVI比率阈值法、滑动平均法、阈值法、最大变化斜率法及Logistic函数拟合法等。不同的物候特征提取方法适应的地区也不一样，因此需要分析与对比。

洞庭湖流域在长江流域水系中占十分重要的比重，也是我国商品粮的主要生产基地，棉麻等经济作物也非常丰富。为此，本文以洞庭湖流域为研究区域，基于时间序列EVI数据，采用滑动平均法、动态阈值法、最大变化斜率法以及Logistic函数拟合法4种方法进行物候参数提取，并对四种物候参数提取方法进行对比分析，为洞庭湖流域植被物候提供技术方法支持。

1 研究区域与数据

洞庭湖流域即洞庭湖水系流经的广大地区，属于典型的亚热带季风气候，年均气温 16.5～17.0 ℃。流域多年平均年降水量1427 mm，多年平均年径流量为2016亿m3，约占长江流域地表水资源的21%，其比重为长江流域各水系之首。洞庭湖流域覆盖湖南省大部和湖北省、广西壮族自治区、贵州省和重庆市部分地区，流域面积约262800 km2，占长江流域面积的14%。洞庭湖平原区为湖南农业主产区，以种植粮食、棉花为主，也是中国主要淡水养殖区之一，耕地面积约占湖南全省的1/6，居住人口约占湖南全省1/9，水、气、生物等资源丰富。

本文选用已进行畸变校正和定标之后的MODIS 13Q1产品的EVI数据与MODIS 12Q1遥感数据经过数据提取、影像镶嵌、投影转换、裁剪等一系列处理后进行研究。其中，用TIMESAT软件对MODIS 13Q1产品中的EVI数据进行了S-G滤波重构处理；MODIS 12Q1为土地覆盖数据，主要包括林地、草地、耕地、水域、城市用地等类型。

2 研究方法

2.1 滑动平均法

滑动平均法主要是将植被指数时间序列数据突然增长时对应的时间代表植物开始生长的时期，是在回归滑动平均模型的基础上，将植被指数时间序列曲线与计算后的滑动平均曲线进行比较，以反映实际植被指数时间序列曲线和既定趋势之间的偏差。滑动平均模型的公式为：

Yt=(Xt+Xt-1+Xt-2+…+Xt-(w-1))/w

(1)

其中，Yt是时间t时滑动平均模型的计算值，Xt表示时间t时的标准植被指数时间序列值，表示滑动平均模型的窗口大小，即滑动平均时间间隔[11-13]。

2.2 动态阈值法

动态阈值法是在固定阈值法的基础上，由Jonsson等提出，相较于固定阈值法其考虑了植被指数时间序列值的变化幅度[14-17]，在一定程度上消除了植被类型和土壤背景值的影响。具体研究时，根据植被指数时间序列曲线的变化特点，将曲线最初阶段的植被指数值快速增加处定义为植被生长开始期，可具体确定为植被指数时间序列重构曲线的左边起始值或最小值增长到指定程度(如整体增幅的10%)，其他物候特征说明见表1。

2.3 最大变化斜率法

最大变化斜率法将植物在春季开始迅速生长的时期定义为植被生长季的开始。具体的方法步骤为：第一步，建立标准的植被指数时间序列数据；第二步，利用最大变化斜率模型计算植物的生长季开始点，并设置条件阈值剔除掉结果中不合理的计算值。最大变化斜率法的公式为：

Δyt=yt-yt-1，θt=arctg(Δyt)

Δyt+1=yt+1-yt，θt+1=arctg(Δyt+1)

Δθ(t+1)=θt+1-θt

(2)

表1 动态阈值法植被物候特征定义物候参数参数含义生长季开始时间EVI增加至拟合函数的左半部分EVI振幅的20%的时刻生长季结束时间EVI降低至拟合函数的右半部分EVI振幅的20%的时刻生长季长度开始生长时间与停止生长时间之间的时间间隔基准值拟合函数左半部分和右半部分EVI最小值的平均值生长季中期时刻EVI增至拟合函数左半部分EVI 振幅的80%的时刻和EVI降至拟合函数右半部分EVI振幅的80%的时刻的平均值生长季EVI最大值EVI的最大值生长季EVI振幅EVI 峰值与左半部分最小值和右半部分最小值的均值之间的差值生长速度拟合曲线左边振幅20%和80%之间的EVI差值与相应的时间差值的比值减缓速度拟合曲线右边振幅20%和80%之间的EVI差值与相应的时间差值的比值生长季EVI活跃累积量EVI拟合曲线与基准值以上的区域围成的面积生长季EVI总累积量EVI拟合曲线与作物开始生长到生长结束时间段内的分量

其中，yt-1指t-1时的植被指数值，yt是指t时的植被指数值，yt+1是指t+1时植被指数值，Δθt和Δθt+1指t和t+1时的角度变化值[18-22]。

2.4 Logistic函数拟合法

该方法通过曲率变化率的极值点反映植物各物候期的转换[23-26]。具体的方法步骤为：

第一步，建立标准的植被指数时间序列数据，并根据实际情况划分出植被指数值持续增长或下降的阶段。

第二步，利用Logistic函数对划分出来的各个阶段进行拟合，生成拟合曲线，拟合函数为：

(3)

其中，t为具体的时间值，y(t)为随时间t变化的植被指数拟合值，d为标准的植被指数值；a，b，c为拟合函数的参数。

第三步，计算Logistic拟合后的植被指数时间序列曲线的曲率K及变化率，将曲率变化率的极大值定义为植物生长的开始时间。曲率K及变化率K′的计算公式如下：

(4)

其中，z=e(a+bt)，α为曲线在时间t时的单位弧长转过的角度，s为单位弧长。

(5)

本文利用MATLAB软件，通过编程，实现上述4种方法，并利用4种方法对洞庭湖流域植被的物候信息进行了提取与对比分析。

3 结果与分析

3.1 不同方法下的物候特征计算结果

利用4种方法提取的2005与2015年洞庭湖流域的部分物候特征如表2，包括生长季开始(start of season， SOS)、生长季结束(end of season， EOS)以及生长季长度(length of season， LOS)，EVI最大值(max of EVI， MOE)。

表2 不同方法下的物候特征参数d提取方法SOSEOSLOSMOE2005年2015年2005年2015年2005年2015年2005年2015年动态阈值法24753213532982790.830.91最大变化斜率法33332252731922400.830.91Logistic函数拟合法1849324916366700.830.91滑动平均法51423262372761960.830.91 注:像元坐标为x=265,y=415

依据湖南气象网农业气象旬、月报上的相关实地观测资料，洞庭湖流域主要农作物为油菜、水稻、棉花和玉米，其物候分别为：油菜一般为上一年9月底10月初种植，翌年5月中下旬左右成熟收割；水稻当年4月上旬左右为返青期，双季稻7月左右收割第一季，然后种第二季，11月份收割，而单季稻于当年10月份左右收割；棉花当年4月左右为播种出苗期，9月左右开始收获，直到11月份；春玉米当年四月上旬左右为播种出苗期，6月左右为成熟期。考虑到研究区域内有一年一熟、一年两熟的作物，这次实验提取的为年内植被生长的开始期和结束期。由于实验所采用的数据分辨率为250 m，一个像元内的植被可能并不是单一的，综合对比实验提取的物候特征参数和实地观测资料，可以看出：①动态阈值法提取2005年像元的SOS值为第24天，早于4种农作物观测数据的生长开始时间；提取2015年像元的EOS值为第353天，相比其他三种方法的提取结果过晚，也晚于4种农作物观测数据生长结束时间。②最大变化斜率法提取的2005年和2015年的SOS普遍过早，且由于最大变化斜率法的模型问题，影像一个时间间隔内的斜率变化角相同，使得最终的计算结果存在误差。③Logistic函数拟合法，因研究区域内多种土地覆盖类型，针对个别土地覆盖类型如水体等，其植被指数时间序列曲线呈锯齿状重复出现，若直接采用传统的Logistic模型方法对原始的植被指数时间序列影像进行分段拟合，MATLAB在求解Logistic函数时会出现问题，导致出错，因此利用MOE前后的植被指数时间序列数据，分段用植被指数时间序列累积曲线来进行Logistic函数拟合。Logistic函数拟合法提取2005年像元的SOS值为第184天，相比其他三种方法的提取结果过晚，也晚于4种农作物观测数据的生长开始时间；提取2015年像元的EOS值为第163天，相比其他三种方法的提取结果过早，且也早于4种农作物的综合生长结束时间。④滑动平均法提取2005年和2015年像元的SOS值分别为第51天和第42天，EOS值分别为第326天和237天，对比其他三种方法，该方法的提取结果与4种农作物观测数据生长开始时间和结束时间接近，准确性较高。综合考虑到上述4种作物物候特征参数提取方法的模型原理和实验数据物候参数提取结果，滑动平均法模型简单，效果较好，更适用于整个研究区域的物候特征参数提取。

3.2 近10年洞庭湖流域植被物候特征参数变化分析

选择滑动平均法对2005年与2015年整个洞庭湖流域EVI植被指数时间序列影像进行物候特征参数提取，提取结果如图1。

从图1可看出，2005年SOS值从0～第120天均有一定区域的分布，但2015年SOS值大部分分布在0～第80天之间，集中在第25～75天之间；且2005年图中中部区域的SOS的值在第80～120天之间的，在2015年时主要分布在第40～80天之间，提前了40 d左右；2005年研究区东北部和南部区域的SOS值在第40～50天之间的，在2015年时主要分布在第60～75天之间，推迟了20 d左右；2005年和2015年，洞庭湖流域EOS值整体上都分布在第220～350天之间，但2005年研究区北部和中东部区域的EOS值在第300～350天之间的，在2015年时主要分布在第230～320天之间，提前了30～70 d左右；2005年研究区中西部区域的EOS值在第220～260天之间的，在2015年时主要分布在第300～350天之间，推迟了85 d左右； 2005年北部和东部区域的LOS值主要集中在第270～320天之间，在2015年主要分布在第200～270天之间，减短了60 d左右；中西部区域2005年LOS值集中在第200～260天之间的，在2015年分布在第240～300天之间，延长了50 d左右。

图1 2005年和2015年物候特征参数

续图1 2005年和2015年物候特征参数

利用MODIS 12Q1数据与洞庭湖流域物候特征参数数据叠加，并进行统计分析，得到了洞庭湖流域的林地、草地和耕地这三种主要植被覆盖类型的物候特征参数，见表3。

表3 不同植被覆盖类型的物候特征参数d植被覆盖类型SOS均值EOS均值LOS均值2005年2015年2005年2015年2005年2015年林地51.46652.294305.23300.362254.213248.538草地52.16853.205304.646297.831252.93245.098耕地52.85652.6311.253297.647258.847245.522

从表3可以看出，2005年林地、草地、耕地的SOS均值分别为51.466，52.168，52.856；2015年林地、草地、耕地的SOS均值分别为52.294，53.205，52.6；两年不同植被覆盖类型的SOS均值相差不远。2005年林地、草地、耕地的EOS均值分别为305.23，304.646，311.253；2015年林地、草地、耕地的EOS均值分别为300.362，297.831，297.647，所以2015年林地和草地的EOS均值比2005年提前了6 d左右，2015年耕地的EOS均值比2005年提前了14 d左右。2005年林地、草地、耕地的LOS均值分别为254.213，252.93，258.847；2015年林地、草地、耕地的LOS均值分别为248.538，245.098，245.522。分区统计结果显示，两年SOS均值变化不大，2015年EOS均值比2005年提前了6～14 d左右，而2015年LOS值比2005年减短了6～13 d左右。

这些变化引起的原因有气候变化、农作物类型发生变化、耕种方式发生变化等，但从湖南气象网了解到，农作物类型和耕种方式变化不大，气候变化是主要原因。2005年年均气温18 ℃、年降水量约1 300 mm、年日照时数1 400 h左右，2015年年均气温17.9 ℃、年降水量约1 580.5 mm、年日照时数1 180 h左右，降水量较常年偏多。且2015年相比2005年，年初气温异常偏高，出现强“暖冬”现象，清明节前后，却出现“倒春寒”现象，盛夏时气温异常偏低，立冬后又出现罕见冬汛，所以这些气候变化是造成物候特征参数值空间格局发生变化的主要原因，且气候变化对对农作物的影响较大，对林地、草地影响较小。

4 结论

以洞庭湖流域为研究区域，采用MATLAB GUI将四种物候特征参数提取方法通过编程实现，然后分析了四种物候特征参数提取方法的优劣性，最后分析了2005年和2015年物候参数的空间变化及原因。结论如下：

1)采用滑动平均法计算各类植被生长季开始时间简单、有效，并且克服了阈值法在不同土壤背景、光照以及植被类型等条件下需要设置不同阈值的问题，但该方法如果受气候事件如雪融影响，所得结果可能会早于植物实际生长季开始时间。

2)动态阈值法计算生长季开始和结束时间简单、耗时短，但需要多次实验才能确定适合区域植被覆盖类型的最佳阈值。

3)最大变化斜率法充分考虑植被指数时间序列曲线特征，将生长季开始时间评估限制在合理生长期内，以减少对稀疏植被地区不经历主要生长季开始时间事件的错误计算，且设定阈值条件可提高植物生长季开始时间的计算效率与精度，但相关阈值的设置须主观综合考虑研究区域、植被类型及背景条件的具体实际情况。阈值条件选择合理与否，直接影响计算结果及运行效率。

4)Logistic函数拟合法监测大区域植被动态变化灵活方便，结果理想；逐个处理像素时不需设置阈值或经验限制条件，因此具有适用性；由于是分段式拟合，能对一年内农田和半干旱地区常见的多熟物候行为进行监测；不受植被物候历法时期限制，为监测接近真实时间的植被物候提供了可能，但实际的植被指数时间序列曲线并非理想的规则曲线，突变点较多，拟合精度的高低直接影响植被物候参数的确定。

5)2005年与2015年洞庭湖流域物候特征参数的空间格局的变化主要是由气候变化引起，且气候变化对林地和草地的影响较小，对耕地的影响较大。

总体而言，从多时相遥感影像中获取的植被指数时间序列数据，经时间序列重构后可以较为准确地提取植物物物候特征参数，反映植物的生长过程，并可以反映出区域气候变化，适合大范围、快速监测植物的物候变化。