基于多时相遥感数据的我国东北地区耕地生长季特征提取与空间格局分析

基于多时相遥感数据的我国东北地区耕地生长季特征提取与空间格局分析

农业物候期作为一个重要的农业生态系统特征，不仅是重要的农业生产、管理和规划决策的重要依据，也是全球气候变化研究和耕地生态循环系统的重要组成部分。农业物候期监测主要是发现作物在形态上发生显著变化所对应的日期,其中作物的出苗期、抽穗期、成熟期是关键的物候期。物候期特征可以通过田间观测法、积温法和遥感法等多种途径获得。其中利用遥感法提取作物物候期特征的基本假设是:从多时相遥感影像中获取的NDVI数据可以反映作物生长的动态变化特征,即NDVI的时序变化对应于作物返青生长与成熟等活动过程。通常认为经函数拟合方法重构的NDVI时序变化曲线,可以准确地反映农作物的出苗、拔节、抽穗、成熟等物理过程,适合大范围、快速监测作物的关键物候期。目前,基于NOAA/AVHRR,SPOT/VGT以及EOS/MODIS等中高分辨率NDVI时序资料已经在作物生长过程监测、作物面积提取、作物类型识别和作物种植制度判定等方面得到了广泛的应用。作为我国重要的商品粮生产基地,东北地区是玉米、水稻等粮食作物和大豆、甜菜等经济作物的主产区,在我国粮食安全保障体系和农业生产中占有重要地位。在我国面临气候变暖问题的区域中,东北地区又是增温最显著的地区之一。目前已有大量关于该地区气候变化对作物生长条件、播种期和生育期、作物产量以及种植格局的影响研究,结果表明作物品种和播种期调整是适应气候变化的重要措施。传统作物物候期研究多从农业气候资源的角度出发来分析气候变化对作物播种期和品种布局的影响,主要通过比较不同品种熟型作物从播种到成熟所需的热量条件,探讨不同作物播种期和品种布局与气候条件的匹配性,而实际农业生产中由于受到社会因素及人为因素的影响,所种植作物的实际物候期与理想条件下的适宜物候期并不完全对应,从而导致带来更多不确定性因素。此外,由于传统基于站点的农业物候观测数据无法及时全面反映区域尺度上作物生长过程对外界环境条件变化(诸如区域温度、降水和日照时间等)的时空响应,因此,快速获取区域尺度上作物物候期特征的时空动态信息,纠正品种熟型不搭和播种期不适的情况,使东北作物在气候变暖背景下获得最优产量对于制定未来气候变化农业适应策略尤为重要。本文基于多时相遥感信息可以有效描述区域尺度上作物物候期特征的原理,利用1998—2009年期间的SPOT/VGTNDVI逐旬时间序列数据,采用非对称高斯函数拟合方法对数据序列进行平滑重构处理,进而根据拟合曲线变化特征定义和提取了东北地区耕地生长季特征参数,并利用农业物候观测数据对提取结果进行了对比验证研究。在此基础上重点分析了近10年来东北地区耕地典型物候期的空间格局及其时间波动特征。1研究区域和数据1.1气候、气候条件中国东北地区包括辽宁省、吉林省和黑龙江省区域面积为79.18万km2,其中耕地面积为26.44万km2,约占全国耕地总面积的16.5%。该区属寒温带、中温带湿润(半湿润)气候,大部分地区≥0℃积温为2000~4200℃·d,≥10℃积温为1600~3600℃·d,夏季平均气温为20~25℃,无霜期为140~170d,降雨量500~800mm,其中60%集中在7—9月。该地区主要作物包括大豆、玉米、小麦和水稻,种植结构的地域差异很大,品种熟型亦随纬度、地形变化而改变。一般认为,对于上述东北地区主要农作物而言,日平均气温稳定≥7℃初日是适宜播种期,≥10℃初日和终日分别是作物生长发育的起始期和终止期。目前的相关研究显示20世纪80年代以来东北地区农业物候期特征的总体变化趋势为生长季提前(表现为≥7℃和≥10℃初日提前),生长季长度增加,生长季内总积温增加≥10℃积温带北移,玉米和水稻晚熟品种的种植范围北移和东扩。1.2数据及研究方法本文所采用的遥感数据为SPOT/VGT逐旬NDVI最大值合成数据(VGT-S10),时间范围为1998年4月至2009年12月,空间分辨率为1km。该时序数据从比利时佛莱芒技术研究所(FlemishInstituteforTechnologicalResearch,Vito)的VEGET-ATION影像处理中心网站下载(http://free.vgt.vito.be)。该数据已经完成了几何校正、辐射校正、地图投影、状态标识以及大气校正等处理,其中包含了各波段的多光谱数据,NDVI数据以及太阳方位角等多个数据层。整个影像经过投影变换由经纬度转换为双标准纬线等积圆锥投影(ALBERS),通过与区域行政边界掩膜处理,得到研究区范围内NDVI时序数据。本文使用的农业气象观测资料为东北3省79个气象台站(包括黑龙江省34个、吉林省21个和辽宁省24个)1998—2009年的4种主要作物生育期资料(图1),数据项包括作物出苗期、抽穗期或抽雄期(其中不包括大豆)、成熟期以及生育期天数(出苗期与成熟期之间的天数)。该数据来源为国家气象局提供的中国农作物生长发育状况资料数据集,资料根据1991年以来中国农业气象台站上报的农业气象旬月报报文资料整理而得。研究中使用的2000年土地利用数据来自中国科学院全国土地利用遥感调查结果,该数据是由覆盖全国的分辨率为30m的LandsatTM/ETM数据,采用人工交互式图像解译制作的,数据以栅格形式存储,比例尺为1:10万,分辨率为100m,共有26类土地利用类型,其中耕地层数据可以分为水田和旱地。本研究主要提取东北地区的耕地图层用作农业用地分析。2学习方法2.1ndvi时间序列拟合利用多时相NDVI数据集提取作物物候期特征的重要基础是对数据序列进行平滑处理重构,以最大可能减弱噪声影响,减少数据序列的空值点,从而提升其年际间和区域间的可对比性。从各种数据处理方法的效果对比来看,基于谐函数或高斯函数的最小二乘拟合方法能很好地应用于具有明显植被生长与衰落等季节特征的数据序列,而对于季节特征不明显的序列则存在较大的不确定性;滤波方法实现较为简单,能较好地反映数据序列的局部变化特征,受局部误差影响,年际间差异较大。本研究采用的数据平滑方法是非对称性高斯函数拟合方法。该方法是一个从局部最优化拟合到全局拟合的方法,具有很多灵活性,使得重建的NDVI曲线可以较好描述NDVI时序数据中复杂的和微小的变化。具体的步骤如下。首先,提取原始NDVI时序曲线中的谷值和峰值,采用高斯函数分别拟合曲线的左右部分,局部拟合公式为其中,式(1)和(2)中,c1和c2为定义基线和振幅的线性参数g(t;a1,…,a5)是高斯函数,其中,a1是对应时间变量的峰值或谷值的位置参数,参数a2和a3决定函数曲线右半部分的宽度和平度,相应地,参数a4和a5决定函数曲线左半部分的宽度和平度。上述参数可以通过优化函数计算得到。在假定数据误差皆为负向偏差的前提条件下,通过与NDVI数据的上包络曲线拟合,基于最小二乘法原理构建最优化拟合函数计算方程系数。优化函数公式为为数据序列中第i点对应的时间ti与NDVI值Ii;σi为假定已知的第i个数据点的不确定性。如果该确定性未知,则设定该值为常数1。局部拟合函数可以很好地描述最大值与最小值间隔内的NDVI数据变化曲线。但是在曲线突出部,拟合效果不是很好。在此将变化曲线划分成左边谷值间隔区,中部峰值间隔区与右边谷值间隔区,分别以不同的局部拟合函数来描述,最后再利用各局部拟合函数构建整体拟合函数,从而较好地描述整个植被生长期的NDVI变化过程。简单来说,整体拟合就是将各局部拟合函数的特征加以综合的过程。整体拟合函数为式(4)中,[tL,tR]区间是NDVI数据序列整个变化区间,fL(t),fC(t),fR(t)分别代表[tL,tR]区间内左边谷值、中间峰值及右边谷值所对应的局部函数,α(t)和β(t)为位于之间的剪切系数。一般来说,将局部拟合函数整合为全局函数是该方法的重要特征,从而增加了拟合函数的弹性,可以模拟更复杂的时间序列变化特征。研究中应用非对称性高斯函数拟合方法,对1998—2009年间的NDVI时间序列数据进行了平滑处理,得到了可用于描述农作物生长特征的光滑曲线。对比平滑前后的NDVI曲线(图2)可以看出:平滑前NDVI时间序列曲线有很严重的锯齿状波动,直接用于时间序列分析比较困难;平滑后重构的NDVI曲线保持了原有曲线的基本形状,而且去除一些干扰值后,对曲线特征的刻画更为突出。总体而言,平滑处理后重构的NDVI时序曲线有利于提取研究区农作物物候特征信息。2.2动态阈值法:基于生长曲线变化的模型识别经过平滑处理的NDVI时间序列曲线可以反映出作物生长的年内动态变化特征,因此,可以根据NDVI时间序列曲线提取得到作物物候信息,如作物的出苗期、抽穗期和成熟期等。基于NDVI时序数据提取作物物候期的常用方法是NDVI阈值法即根据预先定义的NDVI参考值来确定典型物候期。该方法操作简单,适宜于小区域的作物物候研究,但是由于不同区域和不同植被类型的NDV值存在显著差异,确定适合大尺度物候期提取的单一NDVI阈值是不可行的。张峰等依据作物生长过程线的时间特征,提出利用生长曲线最大上升或下降速率来确定作物生长过程的转折点,以此反映作物的出苗和成熟期。该方法能有效反映作物生长历程的差异,但存在计算过程复杂,计算结果精度受时间分辨率影响大等缺陷,较难应用于大区域物候期特征提取。Jönsson等提出动态阈值法,其改进之处在于通过考虑NDVI季节变化幅度,实现在像元尺度上阈值的动态设定,从而在一定程度上消除土壤背景值和植被类型的影响。目前,该方法已经应用于全球很多区域,取得了较好的应用结果如Heumann等利用20%的阈值研究了非洲大陆撒哈拉和苏丹区域的植被物候变化特征,Beck等[42则利用25%的阈值监测了北半球高纬度的芬诺斯坎底亚半岛的植被物候变化,于信芳等利用20%的动态阈值监测了我国东北地区森林植被物候并进行了验证,吴文斌等亦选择20%的动态阈值提取了我国耕地生长季起始期和华北地区耕地典型物候期。本研究依据作物生长过程曲线的变化特征,在耕地作物生长的NDVI曲线上升阶段,将作物生长起始期定义为曲线初始阶段上升速率急剧增加处具体日期确定为NDVI拟合曲线增长到某一水平(整体增幅的20%)所对应的日期。类似地,定义NDV拟合曲线下降到某一水平(整体增幅的20%)所对应日期为生长季结束期。生长季高峰期定义为拟合曲线上升的极限处所对应的日期。生长季长度则为从生长季开始期至结束期所需的时间长度(图2)。理论上,对于东北地区一年一熟农作物,上述生长季开始期可对应作物的出苗期,生长季高峰期可代表作物的抽穗期或抽雄期(玉米),生长季结束期可反映作物的成熟期,而生长季长度则和作物的生育期天数相对应(表1)。3结果分析3.1主要收获季节的基本特征3.1.1出苗期和出/短时分布通过对近10年来东北地区农业气象观测资料的汇总分析,可以发现区内主要农作物的出苗期存在明显差异(图3(a))。从平均出苗期的分布区间来看,一季稻的出苗期最早,主要出现在4月下旬;大豆的出苗期最晚,多分布在5月下旬;春小麦和春玉米介于其间,分别处于5月上旬和5月中旬;整体而言,研究区主要作物的出苗期从4月下旬一直延续到6月上旬,总体时间跨度为50d左右。从平均出苗期的年际变化来看,上述4种作物的年际变化幅度基本不超过10d,总体变化趋势基本一致没有出现交叠的区间,其中一季稻的年际变化幅度最小,相对其他作物其出苗期最为平稳。3.1.2春、玉米生长规律从作物的抽穗期的多年平均水平状况来看,区内3种主要粮食作物的抽穗期具有显著的差异(图3(b))。其中春小麦的抽穗期早在6月下旬,而一季稻的抽穗期则晚至8月上旬,春玉米则介于其间主要分布于7月下旬。整体而言,区内耕地作物抽穗期时间跨度为50d左右。此外,3种作物抽穗期的年际变化趋势均不明显,总体变幅均小于10d。3.1.3春小麦生育期特点由于研究区内种植有以春小麦为主的夏收作物和以水稻、玉米和大豆为主的秋收作物,因此,区内的作物成熟期亦存在明显的夏秋两季分布的特点(图3(c))。从多年平均状况来看,除春小麦的成熟期早在7月下旬和8月上旬外,其他3种作物的成熟期集中分布在9月中下旬。整体而言,区内以的夏收作物和秋收作物的成熟期时间跨度分别为20d左右。从平均成熟期的年际变化来看,春小麦的年际变化相对较大,个别年份差异达15d以上;其他秋收作物成熟期的年际变化趋势不明显,总体变幅均小于10d。3.1.4生育期和生育期日期、年份变化作物生育期天数指作物从出苗期至成熟期所需日数,由于区内主要作物出苗期存在明显的先后次序,导致作物生育期天数也存在较大差异(图3(d))。从多年平均状况来看,春小麦的生育期最短,大约在90d左右;一季稻的生育期最长,达150d左右;玉米的生育期略长于大豆,两者皆为130d左右。从平均生育期天数的年际变化来看,一季稻和春小麦的年际变化相对较大,个别年份差异达10d左右;春玉米和大豆的年际变化趋势不明显,总体变幅均小于5d。3.2对回收地区主要农业物候观值的散点分析为了验证基于NDVI拟合曲线提取的生长季特征与农业物候观测值之间的对应关系,本文根据农业气象观测站点的位置,提取了周边地区耕地范围内的主要生长季特征,进一步与农业物候观测值进行散点分析(图4)和相关分析(表2)。3.2.1春不同作物类型种植区相关系数从提取的生长季开始期与作物出苗期的散点分析图(图4(a))可见,两者间具有明显的线性关系,整体相关系数达0.93(表2)。从不同作物类型来看,春小麦种植区域相关系数最高,达0.92;大豆种植区域两者相关系数最低,为0.69;春玉米和一季稻种植区的相关系数介于其间,约为0.82;上述分析表明基于NDVI拟合曲线提取的生长季开始期可以很好地表征不同作物的出苗期。3.2.2种植区域与春小麦产量的相关分析类似于生长季开始期与作物出苗期之间的关系,从生长季高峰期与作物抽穗期的散点分析图(图4(b))可见,两者间亦具有较明显的线性关系,整体相关系数达0.91(表2)。从不同作物类型来看,春小麦种植区域相关系数最高,达0.91;一季稻种植区域两者相关系数次之,为0.87;春玉米种植区的相关系数较低,约为0.77;上述分析表明不同作物的抽穗期可以通过提取的生长季高峰期加以表征。3.2.3种植区域与时期分布相对于其他物候期散点图的差异化分布,区内生长季结束期与作物成熟期的散点分布更为集中(图4(c)),主要分布在8月上旬和9月中下旬两个时段,前者为以春小麦为主的夏收作物种植区,后者为春玉米、一季稻和大豆为主的秋收作物种植区。两者间的线性关系亦较为显著,整体相关系数为0.85(表2)。从不同作物类型来看,春玉米种植区域相关系数最高,达0.85;一季稻种植区域两者相关系数次之,为0.77;春小麦和大豆种植区的相关系数较低,约低于0.75;上述分析表明提取的生长季结束期亦可以表征不同作物的成熟期。3.2.4种植区域与品种间的线性关系从提取的生长季长度与作物生育期天数的散点分析图(图4(d))可见,由于夏收与秋收作物的差异,大致可分为80~100d和120~160d两大区间,前者为以春小麦为主的夏收作物种植区,后者为春玉米、一季稻和大豆为主的秋收作物种植区。相对于其他散点图较显著的线性关系,两者间的线性相关关系表现明显,整体相关系数略低于0.90(表2)。从不同作物类型来看,春玉米种植区域相关系数最高,为0.74;一季稻种植区域两者相关系数次之,约为0.70;春小麦和大豆种植区的相关系数较低,均在0.50左右;上述分析表明提取的生长季长度在不同程度上可表征作物的生育期天数,其中在春玉米和一季稻种植区域内的表达效果较好。3.3主要生长季特征均值的空间分布为了获取东北地区农业物候特征的空间分布信息,本文提取了1998—2009年耕地范围内主要生长季特征均值的空间分布(图5),以期通过分析作物物候特征的空间分异为快速识别作物类型提供研究基础。3.3.1耕地生长季开始期的时空差异从该地区耕地生长季开始期空间分布特征(图5(a))来看,与≥10℃初日从南至北逐渐推迟的空间分布格局不同,生长季开始期空间分布的南北差异并不显著。其中三江平原东部与松嫩平原北部耕地生长季开始期相对较早,多分布在4月中下旬至5月上旬;黑龙江南部以及吉林和辽宁省大部地区生长季开始期均在5月中下旬至6月上旬,该研究结果与吴文斌等利用NOAA-AVHRR序列提取的全国范围内的耕地生长季开始期基本一致。由于东北地区北部耕地区域大量种植生育期较长的一季稻,其出苗期相对较早,而其他区域的玉米、春小麦和大豆出苗期较晚,导致区域内生长季开始期与≥10℃初日在空间上的对应关系不明显。上述结果亦表明东北地区作物种植类型以出苗期较晚的品种为主。3.3.2晚7月至7月东北部从该地区耕地生长季高峰期空间分布特征(图5(b))来看,松嫩平原北部区域、辽宁和吉林的东南部山区的生长季高峰期最早,多分布在6月中下旬至7月上旬,与该区域春小麦的抽穗期分布较为吻合;三江平原大部、松嫩平原中部以及辽河平原中部的生长季高峰期多发生于7月中下旬,与春玉米的抽穗、抽雄期基本一致;而三江平原的局部、松嫩平原南部、辽河平原北部以及辽宁省东南沿海地带的生长季高峰期较晚,多为8月中上旬,与一季稻的抽穗期存在较大重叠。上述结果亦表明东北地区主要粮食作物的抽穗期存在较大的空间差异,而这种空间差异主要受作物种植类型的影响。3.3.3生长季结束期为9月至20月从该地区耕地生生长季结束期空间分布特征(图5(c))来看,相对于其他生长季特征较明显的空间差异化分布,区内大部分地区生长季结束期集中在9月中下旬,与该区域主要种植作物(春玉米、大豆和一季稻)的成熟期基本一致。只有松嫩平原西北部的少量耕地生长季结束期在8月中下旬,亦与以春小麦为主的夏收作物的成熟期相吻合。生长季结束期的空间分布特征也充分反映了研究区域内秋收作物占据主导地位的种植格局。3.3.4与春玉米生育期长度的关系上文的分析表明,区内耕地生长季长度在一定程度上可以反映主要种植作物的生育期长度,从其分布特征(图5(d))来看,区内耕地生长季长度主要可分为3个分布区间:在三江平原的东部、松嫩平原东北部以及辽河平原东南部地区,生长季长度为4~5个月左右(约130~160d),与一季稻以及晚熟品种的春玉米生育期长度具有较好的对应关系;中部的松嫩平原东部、三江平原中西部以及辽河平原东部地区的耕地生长季长度为4个月左右(约110~130d),与大豆以及中熟品种的春玉米生育期长度大致相当;而松嫩平原西北部,如吉林西部和黑龙江西部等,耕地生长季长度较短,一般为3~4个月左右(约90~110d),主要与夏收的春小麦以及春玉米早熟品种的生育期长度相对应。生长季长度的空间分布特征亦进一步反映了研究区作物种植类型以长生育期的作物为主。3.4生长季特征参数变化的区域基于提取的1998—2009年逐年的东北地区生长季特征参数的空间分布信息,本文计算了期间各生长季特征参数值的最大变化幅度(图6),进而识别出物候期发生显著变动的区域,为进一步区分潜在的作物类型转换提供研究基础。3.4.1作物生育期变化的空间分布前文分析表明耕地生长季开始期与不同作物类型的出苗期均存在显著关系,因此,通过分析期间生长季开始期的变化幅度,可以初步了解研究区内作物出苗期的变化类型。近10年来的农业物候观测资料(图3(a))表明,研究区内同一作物类型出苗期的最大年际变化幅度均小于20d,旱地作物类型之间的出苗期差异也都在20d以内,而水田和旱地作物互转可使出苗期波动达30d以上,例如从春玉米转变为一季稻后,出苗期可能提前1个月以上。从耕地生长季开始期变化幅度的空间分布(图6(a))来看,大部分地区的生长季开始期最大变化幅度均在20d以内,表明这些区域作物类型总体保持不变或以旱地作物之间的转换为主;而松嫩平原的中西部地区和三江平原的部分区域生长季起始期变幅达40d以上,表明上述地区作物类型存在水田和旱地作物转换的情形。3.4.2生长季之际影响年际变化的作物时期变化幅度较作为与作物抽穗期关系紧密的生长季高峰期其年际变化幅度也能反映潜在的作物类型转换。基于农业物候观测资料的分析结果(图3(b))表明,15d以内的抽穗期变化幅度多为作物正常的年际波动20~30d的变化幅度通常源自水稻和春玉米之间的转换,超过30d的变化幅度主要由于春小麦与其他作物类型的转换。从耕地生长季高峰期变化幅度的空间分布(图6(b))来看,大部分地区的生长季高峰期最大变化幅度均在20d以内;而松嫩平原的中西部地区、辽河平原和三江平原的部分区域生长季高峰期变幅达30d以上。3.4.3小麦生育期变幅原因生长季结束期可以在一定程度上表征作物的成熟期,基于农业物候观测资料的分析结果(图3(c))表明,20d以内的成熟期变幅多源于作物的年际波动或秋收作物间类型转换,超过30d的变化幅度主要由春小麦与其他作物类型的转换引起。从耕地生长季结束期变化幅度的空间分布(图6(c))来看,大部分地区的生长季结束期最大变化幅度均在20d以内;而松嫩平原的中西部地区、辽河平原和三江平原的部分区域生长季结束期变幅达30d以上。3.4.4年际波动或秋收作物间类型转换前文分析表明提取的耕地生长季长度可以较好地表征作物生育期长度,农业物候观测资料(图3(d))表明,20d以内的变化幅度多为作物的年际波动或秋收作物间类型转换(主要为春玉米和大豆)40~50d的变化幅度通常源自水稻和春玉米或大豆之间的转换,超过50d的变化幅度主要由于春小麦与其他作物类型的转换。从耕地生长季长度变化幅度的空间分布(图6(d))来看,大部分地区的生长季长度最大变化幅度均在30d以内;而松嫩平原的中西部地区、辽河平原和三江平原的部分区域生长季长度变幅达40d以上。4作物物候期特征的宏观时空分异基于多