基于NDVI时序数据的华北地区耕地物候参数时空变化特征*

查 燕，宋 茜，卫 炜，陈仲新，杨 鹏※

(1.中国农业科学院农业资源与农业区划研究所/农业部农业遥感重点实验室，北京 100081；2.农业部规划设计研究院，农业部耕地利用遥感重点实验室，北京 100125)

0 引言

植被物候的动态变化是生物圈对全球气候变化以及陆地水文循环机制变化的动态响应，与陆地生态系统初级生产力和碳循环息息相关，因而成为全球变化的研究热点之一，在全球变化、生态学和农业生产应用等领域发挥了重要作用[1-3]。农作物物候期数据是重要的农业信息，是农业生产、田间管理、计划决策等的重要依据，也是作物模拟模型的重要参数[4-5]。农业物候期信息可以通过田间观测法、积温法和遥感法等多种途径获得[6]。由于卫星遥感具有宏观、高效和便捷的特点，能够提供大范围的多时相重复观测数据，为群落尺度和区域尺度的物候研究提供了有利条件，逐渐成为物候信息获取的重要手段[7]。目前，归一化植被指数(normalized difference vegetation index，NDVI)是应用最为广泛的植被指数，它对植被生长状况、植被覆盖类型等信息非常敏感，能以最直观的形式反映作物从播种、出苗、抽穗、成熟和收割等物理过程，可用于跟踪作物的季节性动态变化[8-9]。

对于大区域而言，SPOT/VGT NDVI 长时间序列数据对植被生长期特征监测有较好地反应[10]，其逐旬的时间分辨率可以弥补空间分辨率的不足，可有效地对农作物长势进行动态监测[11]，因此，SPOT/VGT NDVI长时间序列数据已经在作物生长过程监测、作物面积提取、作物类型识别和作物种植制度判定等方面得到了广泛应用[12-14]。利用SPOT/VGT 的旬合成NDVI数据对区域尺度的耕地复种指数研究较多，如范锦龙等[15]利用谐函数分析法(HANTS)提取了1999～2002年全国耕地复种指数； 唐鹏钦等[16]利用小波变换提取了2007年华北平原的耕地复种指数； 丁明军等[17]利用Savitzky-Golay滤波法研究了1999～2013年中国耕地复种指数的时空变化。对于耕地典型物候期的识别，李正国等[18-19]采用非对称高斯函数法(Asymmetric Gaussian，AG)对华北地区(1998～2007年)和东北地区(1998～2009年)耕地物候期进行了识别； 卫炜等[20]基于2005年SPOT/VGT数据，提取了我国北方耕地物候信息。作为地表物候的重要组成部分，不同于草原、灌丛和森林等自然植被在年内完成一个生长过程，耕地由于受人类活动影响大和多熟种植制度的存在，耕地物候比自然植被物候更为复杂[9]。然而，目前耕地物候的研究都只是停留在某一年或过去一个时期上，很少进行长时间的动态监测。因此开展基于时序植被指数提取耕地生长季特征研究，能更好地反映长时间区域尺度上作物物候期特征的时空差异。

文章以我国华北地区5省(市)为研究区域，基于1km×1km 分辨率的SPOT/VGT逐旬NDVI数据，利用非对称高斯函数拟合法对NDVI数据进行平滑去噪，重建作物生长植被指数曲线，通过比例阈值法提取了1999年和2013年华北地区耕地物候参数(生长季开始期和结束期)，并用波峰阈值法对耕地种植制度进行了识别。研究1999～2013年华北地区物候参数和种植制度的时空变化过程，为中国农业政策制定及管理提供科学依据。

1 研究区域与数据

1.1 研究区域概况

研究区域位于我国粮食主产区之一的华北平原地区，主要包括北京市、天津市、河北省、河南省和山东省，耕地面积约1 700万hm2。大部分区域属于暖温带大陆性季风气候，气候条件良好，年平均温度13.0℃，年均总辐射量5 100～5 300MJ/m2，>0℃年积温4 200～5 500℃，无霜期为170～220d，热量条件可以满足一年一熟和一年两熟农作需要。年降水量为500～900mm，主要集中在夏季的7～9月份，季节分配不均。主要农作物包括冬小麦、玉米、棉花和大豆。

1.2 数据来源

该文所采用的遥感数据是比利时佛莱芒技术研究所(Flemish Institute for Technological Research Vito)免费提供的SPOT VGT 逐旬 NDVI 最大合成数据(http：//www.vito-eodata.be)，空间分辨率为1km，时间范围为1999年和2013年全年，共72景影像。该数据已经完成了一系列几何校正、辐射校正、地图投影、状态标识以及大气校正等处理。该数据集包含了各波段的反射率、NDVI和太阳方位角等多个数据层，利用VITO提供的处理工具VGTExtract提取其中的NDVI数据层，并将影像由经纬度转换为双标准纬线Albers投影，并通过与区域行政边界掩膜处理，得到研究区域内各时相的遥感影像。数据产品中NDVI的DN值范围为0～250，通过公式NDVI=DN×0.004-0.1转换成标准植被指数值。

耕地数据来自中国科学院遥感应用研究所等7个单位2005年共同完成的中国1： 25万土地覆盖遥感调查与监测数据库。该数据库内容包括森林、草地、农田、聚落、湿地与水体、荒漠等6个一级类型和25个二级类型，其中农田数据可分为水田、水浇地和旱地，空间分辨率为100m。研究将华北地区的农田数据提取出来作为耕地图层，转换其地图投影为阿尔伯斯等面积圆锥投影并将空间分辨率重采样到1 km使之与VGT-S10 NDVI数据相匹配，从而得到研究区域的耕地分布数据(图1)。

图1 研究区域耕地分布

2 研究方法

2.1 NDVI时序数据预处理

尽管VGT-S10 NDVI产品经过最大值合成处理后能够消除部分云、气溶胶、太阳高度角等干扰因素的影响，但仍然残留了许多噪声，不能清晰地反映出耕地像元作物的生长过程，因此必须对NDVI时间序列数据进行去噪和平滑等处理。

该文采用非对称高斯函数来分段模拟植被的生长过程，最后通过平滑连接各段高斯拟合曲线实现时间序列的重构[21]。该方法是一种由局部最优化拟合到全局拟合的方法，具有一定的灵活性，避免了整体数据对局部拟合的干扰，使得拟合后的NDVI曲线可以较好描述NDVI时序数据中复杂的和微小的变化，更加接近真实情况[22]。具体步骤如下：

首先，提取原始时序数据曲线中的谷值和峰值，采用高斯函数分别拟合曲线的左右部分。局部拟合函数为：

f(t)=f(t；c1，c2，a1，…，a5)=c1+c2g(t；a1，…，a5)

(1)

式(1)中，g(t；a1，…，a5)为高斯函数；c1和c2为控制曲线的基准和幅度；a1是对应时间变量t的峰值或谷值的位置参数；a2、a3和a4、a5分别控制曲线右半部分和左半部分的宽度和平度。

局部拟合函数可以很好地描述最大值与最小值间隔内的NDVI数据变化曲线。但是在曲线突出部，拟合效果不是很好。在此将变化曲线划分成左边谷值间隔区、中部峰值间隔区与右边谷值间隔区，分别以不同的局部拟合函数来描述，最后再利用各局部拟合函数构建整体拟合函数，从而较好地描述整个植被生长期的NDVI变化过程。整体拟合函数为：

(2)

式(2),[tL，tR]区间是NDVI数据中待拟合部分的变化区间；fL(t)、fC(t)和fR(t)分别代表[tL，tR]区间内左边谷值、中间峰值及右边谷值对应的局部拟合函数；α(t)和β(t)为介于0～1的剪切系数。

2.2 作物物候特征参数提取

经过去噪和平滑处理的NDVI时间序列曲线可以很好地反映作物生长的年内动态变化特征，因此根据NDVI时间序列曲线能提取出作物主要的物候信息，如作物的出苗期或返青期、抽穗期和收获期等特征参数。动态阈值法是目前提取作物物候特征参数常用的方法，其考虑了NDVI季节变化幅度，在像元尺度上对阈值动态设定，从而一定程度上消除了土壤背景值和植被类型的影响。

在我国北方地区的耕地生长季开始前和结束后，可能存在一些影响物候参数准确提取的因素，如冬小麦种植区小麦冬前峰的影响。为避免这些因素的干扰，该文根据先验知识对NDVI时间序列数据进行“掐头去尾”处理，即将1月、11月和12月的数据从参与分析计算的数据中加以剔除。另外由于耕地物候不同于自然植被物候，其受人类耕作活动的影响较大，参考我国北方地区的作物物候历将生长季开始和结束的阈值分别设置为10%和50%[20]。因此，当NDVI时间序列曲线达到的上升阶段振幅的10%时即认为是生长季开始。类似地，定义NDVI时间序列曲线达到下降阶段振幅的50%时为生长季结束。

2.3 作物种植制度提取

每种作物生长周期都对应着NDVI时间序列曲线上的一个波峰，多熟种植制度下会有多个生长周期，波峰的个数即代表了生长周期的个数。其中一年一季作物的NDVI在一年内形成明显的单峰曲线，一年两季作物的NDVI形成双峰曲线，一年三季作物的NDVI形成三峰曲线。因此，可以通过NDVI时间序列曲线上波峰的数目来确定作物种植制度[23]。然而在实际的NDVI时间序列曲线上除了表征作物生长周期的有效波峰外，还存在着一部分由噪声引起的干扰波峰。这些干扰波峰使得拟合曲线上的波峰个数比实际偏多，导致对生长周期个数的误判。前人研究表明，对于农作物而言其生长季的NDVI峰值通常会达到0.4以上[4， 7]。因此，该文选择0.4作为波峰阈值对探测到的峰值进行判定取舍，只有峰值大于0.4的波峰才被认为是能够真实反映作物生长周期的有效波峰，每个耕地像元的种植制度由该像元NDVI时间序列曲线上有效波峰的个数来确定。

3 结果与分析

3.1 农作物典型物候期的空间分布特征

3.1.1 作物返青/出苗期

图2 华北地区作物返青/出苗期空间分布

图2a和2b分别描述了1999年华北地区第一季作物返青/出苗期和第二季作物出苗期空间分布。从分布特征看，第一季、第二季作物生长开始期都存在明显的空间差异，具有从南向北逐渐推迟的空间特征。第一季作物返青/出苗期最早发生在2月上旬，主要是河南南部少数地区； 大部分一年两熟区，如河南中部、北部，河北中部、南部的低平原区、太行山和燕山山前平原区，以及山东西部、西北部平原区，作物出苗期或返青期多开始于3月上旬至4月上旬。河北省北部和东部一年一熟区的作物生长季开始约在4月下旬至6月上旬，而山东省沂蒙山地丘陵区和鲁东丘陵区，作物返青/出苗期多开始于5月上旬至6月下旬，也有少部分区域作物返青/出苗期开始较早，大约在3月上旬至4月下旬。北京和天津地区作物返青/出苗期多发生在4月上旬至6月上旬。第二季作物的出苗期多在6 月下旬至7月上旬，河南南部少数地区发生5 月下旬至6 月上旬。

图3 夏玉米生长季开始期估算值与地面站点观测数据对比散点

2013年华北地区第一季、第二季作物返青/出苗期空间分布与1999年相比，河南大部分地区、山东鲁西、鲁西北部平原区，河北南部地区，第一季作物返青/出苗期明显提前，多发生在2月上旬至3月上旬。河北北部、北京、天津等地区，第一季作物返青/出苗期变化不大。而第二季作物出苗期明显提前的是河南南部和中部地区，由1999年6月下旬7月上旬，提前至6月上旬。有研究表明，华北平原近50年来，年平均气温整体呈现显著上升趋势，尤其春季和冬季增温对年均温增温的贡献较大[24]。华北平原夏玉米播种至出苗期的平均温度在近30 年的上升幅度为0.51℃/10年(P<0.05)，其中河北中部的石家庄、保定、廊坊和唐山该期的平均温度增加显著，平均增幅达1.04℃/10年(P<0.05)[25]。因此，华北地区作物生长开始期提前可能与气候变暖关系密切。

图4 华北地区作物成熟期空间分布

图5 华北地区农作物种植制度的空间分布

图6 华北地区农作物种植制度变化的空间分布 图7 1999～2013年华北地区各省冬小麦播种面积与一年两熟制耕地像元数对比

冬小麦夏玉米是华北地区一年两熟制的典型作物，而夏玉米是最主要的第二季作物。利用2013年SPOT/VGT数据提取的第二季作物生长开始期与25个地面站点夏玉米观测数据对比如图3所示，二者相关系数为0.40，平均绝对误差(Mean Absolute Error)和均方根误差(Root Mean Square Error)分别为5.2和7.09，说明遥感提取结果与实际情况比较吻合，可以用来分析典型作物物候期的空间变化。

3.1.2 作物成熟期

华北地区作物成熟期空间分布也存在明显的地域差异，河南除豫西和豫南的山地丘陵区以外，大部分地区作物成熟比较早。河北中南部的太行山和燕山山前平原区，山东的鲁西和鲁西北平原区等，由于第一季作物生长期开始早，因此其成熟期结束也较早，主要在5月下旬至6月下旬。第一季作物成熟早为第二季作物完成整个生长周期提供了足够的热量资源保障。而河北的中部和北部、山东中部和东部，以及北京市和天津市等地区，作物返青/出苗期较晚，其生长结束期也就较晚，多发生在9月上旬至10月中上旬(图4a，图4c)。华北地区第二季作物成熟期比较集中，基本都集中在9月上旬至10月上旬，南部地区总体上要早于北部地区(图4b，图4d)。华北地区作物成熟期空间分布差异与该区域耕地第一季作物的生长开始期差异大有关外，还与不同区域的作物种类和作物生长周期不一样有关[4]。

3.2 农作物种植制度的空间分布特征

图5是利用NDVI数据分别提取的1999年和2013年华北地区耕地种植制度空间分布图。从图5可以看出，华北地区耕地种植制度以一年两熟为主，一年一熟区主要分布在河北省北部和中部地区，山东省中部的沂蒙山地丘陵区和东部的鲁东丘陵区，以及河南省西部和南部的山地丘陵区。一年两熟区主要分布在河北省南部、山东省西部、西北部的平原区，以及河南省北部地区。北京市和天津市耕地主要为一年一熟，两熟制区域面积很小。华北地区作物种植制度空间分布格局和该区域自然地理分布特征具有很好的一致性，因为耕地种植制度在很大程度上受水、热、光和土等自然资源共同作用和影响[18]。

由于受农业技术、社会经济发展，以及全球气候变化等因素的综合影响，经过15年的耕种(1999～2013年)，华北地区农作物种植制度的空间格局发生了很大的变化(图6)，河北省中部和南部，以及河南省南部有部分区域由一年两熟变成了一年一熟； 而河南西部、山东中部和东部地区，有零星的区域由一熟制变成了两熟制； 华北地区北部受热量资源制约，仍旧保持一年一熟制不变。与1999年相比， 2013年两熟制种植面积下降了21.1%，而一熟制种植面积增加了38.7%。杨婷等[12]利用遥感数据研究中国近些年熟制变化结果表明，中国两熟制区中2 819个像元变回了一熟制，这种现象在河北和山东较明显。

华北地区绝大多数一年两熟区，冬小麦是最主要的第一季作物。因此，为了验证上述种植制度空间提取结果的可靠性，该文利用《中国农业统计资料》华北地区各省15年(1999～2013年)的冬小麦播种面积与各省一年两熟制耕地像元的数量进行对比。从图7可以看出，冬小麦播种面积与一年两熟区具有极显著的线性相关关系，R2为0.97，说明遥感提取结果符合实际情况，可以用来分析区域种植制度的空间变化。

4 讨论

基于NOAA/AVHRR、SPOT/VGT和EOS/MODIS等中高分辨率NDVI时间序列数据已经在作物生长过程监测、作物类型识别、作物种植制度判定等方面得到了广泛应用。NOAA/AVHRR的NDVI数据产品其空间分辨率通常为8km，由于我国华北地区耕地种植制度复杂多样， 8km 的空间分辨率在一定程度上影响了耕地种植制度和物候信息结果的准确性[4]。MODIS提供了500m和1km两种相对较高空间分辨率的NDVI数据产品，但是其时间分辨率为16d。SPOT/VGT的VGT-S10 NDVI产品具有适中的1km的空间分辨率和10d的时间分辨率，相对上述两种数据而言更加适合我国北方地区的耕地物候研究[20]。

由于小麦冬前峰的存在可能使得重构的NDVI拟合曲线形状发生变化，从而导致提取的生长季开始期比实际提前，而生长季结束期比实际推迟。在该研究中根据先验知识将耕地生长季之外的1月、11月和12月的NDVI时间序列数据加以剔除，这种“掐头去尾”处理，能够避免小麦冬前峰的影响，进而准确获取耕地物候信息。

耕地种植制度与物候分布格局和自然地理环境紧密相关，尤其是不同区域的温度、降水和光照等气候资源的禀赋和匹配程度对耕地种植结构和作物布局有很大影响[4， 26]。该文仅对物候期和种植制度时空变化特征进行了分析，未对这些影响因素进行深入探讨。因此，未来计划通过分析区域内主要农业气候资源，如≥0℃积温、平均气温、温度生长期天数、降水等的时空变化特征，进一步完善耕地物候期对农业气候资源特征的响应机制研究。

5 结论

该文基于SPOT/VGT数据提取了我国华北地区耕地典型物候期的空间格局，并对该区域种植制度进行了识别。研究结果表明，华北地区作物生长开始期和结束期都存在明显的空间差异， 15年来河北北部、北京、天津等地区，第一季作物返青/出苗期变化不大。而河南南部和中部地区， 2013年第二季作物出苗期明显提前，由1999年6月下旬7月上旬，提前至6月上旬。我国华北地区种植制度仍以一年二熟制为主，华北地区北部受热量资源制约，仍旧保持一年一熟制不变。与1999年相比， 2013年两熟制种植面积下降了21.1%，而一熟制种植面积增加了38.7%。总之，华北地区耕地物候的时空变化与种植制度密切相关，同时也受到自然资源和人类活动的共同影响。

[1] 侯学会， 牛铮，高帅，等.基于SPOT-VGT NDVI 时间序列的农牧交错带植被物候监测.农业工程学报， 2013, 29(1): 142～150

[2] Cleland E E，Chuine I，Menzel A，et al.Shifting plant phenology in response to global change.Trends in Ecology and Evolution， 2007， 22： 357～365

[3] Weiss J L，Gutzler D S，Coonrod J E A，et al.Seasonal and inter-annual relationships between vegetation and climate in central New Mexico，USA.Journal of Arid Environments， 2004, 57(4): 507～534

[4] 吴文斌， 杨鹏，唐华俊，等．基于NDVI 数据的华北地区耕地物候空间格局研究．中国农业科学， 2009, 42(2): 552～560

[5] 辛景峰， 宇振荣，Driessen P M.利用NOAA NDVI数据集监测冬小麦生育期的研究.遥感学报， 2001, 6(5): 442～447

[6] 张峰， 吴炳方，刘成林，等.利用时序植被指数监测作物物候的方法研究.农业工程学报， 2004, 20(1): 155～159

[7] Wu Wenbin，Yang Peng，Tang Huajun，et al.Characterizing spatial patterns of phenology in cropland of China based on remotely sensed data.Agricultural Sciences in China， 2010, 9(1): 101～112

[8] 张峰， 吴炳方，刘成林，等.区域作物生长过程的遥感提取方法，遥感学报， 2004, 8(6): 515～528

[9] 李正国， 唐华俊，杨鹏，等．植被物候特征的遥感提取与农业应用综述.中国农业资源与区划， 2012, 33(5): 20～28

[10]Chuai X W，Huang X J，Wang W J，et al.NDVI，temperature and precipitation changes and their relationships with different vegetation types during 1998～2007 in Inner Mongolia，China.International Journal of Climatology， 2013， 33： 1696～1706

[11]李杨， 江南，侍昊，等.基于SPOT/VGT NDVI 的大区域农作物空间分布.农业工程学报， 2010, 26(12): 242～247

[12]杨婷， 赵文利，王哲怡，等.基于遥感影像 NDVI 数据的中国种植制度分布变化.中国农业科学， 2015, 48(10): 1915～1925

[13]欧立业， 罗烈琴，易明华.基于 NDVI时间序列数据的江西省水稻种植制度变化研究.测绘与空间地理信息， 2016, 39(4): 63～66

[14]Li Zhengguo，Tang Huajun，Yang Peng，et al.Spatio-temporal responses of cropland phenophases to climate change in Northeast China.Journal of Geographical Sciences， 2012, 22(1): 29～45

[15]范锦龙， 吴炳方.复种指数遥感监测方法.遥感学报， 2004, 8(6): 628～636

[16]唐鹏钦， 吴文斌，姚艳敏，等.基于小波变换的华北平原耕地复种指数提取.农业工程学报， 2011, 27(7): 220～225

[17]丁明军， 陈倩，辛良杰，等.1999～2013年中国耕地复种指数的时空演变格局.地理学报， 2015, 70(7): 1080～1090

[18]李正国， 杨鹏，周清波，等.基于时序植被指数的华北地区作物物候期/种植制度的时空格局特征.生态学报， 2009, 29(11): 6216～6225

[19]李正国， 唐华俊，杨鹏，等．基于时序植被指数的东北地区耕地生长季特征识别与应用研究．北京大学学报：自然科学版， 2011, 47(5): 882～892

[20]卫炜， 吴文斌，李正国，等.基于SPOT/VGT数据的中国北方耕地物候提取研究.中国农业资源与区划， 2016, 37(4): 77～86

[21]Jönsson P，Eklundh L.Seasonality extraction and noise removal by function fitting to time-series of satellite sensor data，IEEE Transactions of Geoscience and Remote Sensing.2002, 40(8): 1824～1832

[22]Jönsson P，Eklundh L.TIMESAT a program for analyzing time-series of satellite sensor data.Computers & Geosciences， 2004， 30： 833～845

[23]闫慧敏， 曹明奎，刘纪远，等.基于多时相遥感信息的中国农业种植制度空间格局研究.农业工程学报， 2005, 21(4): 85～90

[24]阿多， 熊凯，赵文吉，等.1960～2013年华北平原气候变化时空特征及其对太阳活动和大气环境变化的响应.地理科学， 2016, 36(10): 1555～1564

[25]胡洵， 王靖.华北平原夏玉米各生育阶段农业气候要素变化特征.干旱地区农业研究， 2015, 33(4): 251～258

[26]李正国， 唐华俊，杨鹏，等.东北三省耕地物候期对热量资源变化的响应.地理学报， 2011, 66(7): 928～939