基于长时间序列MODIS影像的黑龙江省主要农作物种植结构变化研究

李 晖，陈雨琪，张紫滟，武洪峰，刘占宇

(1. 杭州师范大学遥感与地球科学研究院，浙江 杭州 311121； 2. 黑龙江省农垦科学院情报研究所，黑龙江 哈尔滨 154005)

0 引言

随着人口的快速增长和城市化趋势的不断增强，世界的农作物分布格局都在发生大小不一的变化.农作物的种植结构信息、空间分布格局和农业精细化管理对于保障国家粮食安全、应对全球农情变化以及农业可持续发展等具有重要意义.农作物种植结构是空间分布格局的重要组成部分，它用来描述一块区域或生产组织内农作物的组成和布局，即主要农作物类型(种什么)和空间分布(种在哪里)[1-3].快速、准确地获取耕作区域农作物的空间分布信息和时空动态变化信息不仅是监测农作物种植面积、预测区域农作物产量、保障区域粮食安全平衡的主要数据来源，也是进行农作物种植结构布局优化的重要依据[4].同时，作为可持续化发展的重要形式，农业土地的开垦利用、作物种植更替以及与其他地物的相互转换，必然会对整个生态环境包括大气圈、水圈和生物圈产生不同程度的影响[5].遥感技术具有覆盖范围大、重访周期短、现势性强和效率高等优点，比传统的人工调查更能准确直观地了解农作物的分布特点，因此早已广泛应用于大尺度范围的作物空间分布信息提取中[6].20世纪70年代，美国、日本、澳大利亚等国家逐渐开始对水稻、大小麦、大豆等农作物空间分布进行遥感监测[7-8]，20世纪90年代，我国也逐步建立了专业的农作物遥感监测平台.随着MODIS(传感器)到Landsat(平台)、TM(传感器)、SPOT(平台)，再到国产高分系列卫星的普及应用，时间和空间分辨率不断提升，为农作物种植信息的精准提取带来了便利.

在农作物的遥感识别提取方面，众多学者进行了深入研究.在方法上可以分为3类：1)基于单一影像识别法.此方法适用于作物种植结构比较简单的区域，通过单幅影像或两到三幅影像对作物的种植信息进行提取，其重点是找到待分作物的光谱信息特征与其他地物的区别.郑长春等[9]基于SPOT-5数据，利用简单决策树结合归一化植被指数NDVI、光谱波段和背景差异等特征量，对黑龙江省852农场主要作物种植信息进行提取.Mathur等[10-11]利用IRS-1D数据，采用支持向量机(SVM)的方法结合光谱波段分析，对印度旁遮普地区的棉花和水稻种植信息进行提取.2)基于时间序列影像识别法.农业土地系统是多种农作物通过一种特定的规律种植并组合形成的空间表达，由于不同作物生长发育时间和阶段不同，单一遥感影像难以覆盖所有作物生长发育期和光谱识别的最佳时相[3].因此，利用时间序列遥感影像表现物候变化规律从而进行农作物遥感提取成为目前最为主流的技术方法[10-12].黄青等[13]基于MODIS-NDVI时间序列数据集，通过分析东北三省主要作物的物候历及时序变化特征，构建作物遥感提取模型并通过调整阈值的方法对东北三省主要作物种植信息进行提取.Zhang等[14]根据时序MODIS影像，结合地表温度产品LST(land surface temperature)，对中国和印度2000—2015年水稻种植面积进行提取，并对两国16年水稻面积的变化趋势进行分析.Atzberger等[15]基于AVHRR影像，通过神经网络算法自动计算NDVI时间序列曲线和端元丰度之间的非线性关系，对意大利托斯卡尼区域进行作物种植信息自动提取.3)基于遥感影像数据与统计数据结合识别法.遥感技术为有效提取作物信息发挥了很大作用，但是由于遥感数据的获取方式、大气干扰、混合像元、分辨率等因素的影响，只采用遥感技术对大尺度区域进行作物提取仍然有一定的缺陷[16-17].Leff等[18]采集已有全球范围的遥感耕地数据，结合国家、省级和县等多个行政单位的作物统计资料，通过非耕地区域掩膜、亚像素面积计算等方法，对20世纪90年代全球5分栅格尺度下的18种农作物进行空间分布信息提取.Fischer等[19]利用GAEZ模型，根据全球耕地分布图、人口密度图等信息，采用交叉信息熵理论和方法将作物统计信息分成5分栅格尺度像元，得到全球农业生态区23类农作物的空间分布特征.

针对目前农作物空间分布信息遥感提取的研究现状，对于大面积耕作地区多年份、多作物的研究较少，本文以整个黑龙江省为研究区，农作物种植面积空间分布的快速提取为主要目标，根据不同作物的波谱特性和物候特征，选取 MODIS-NDVI时间序列数据集的合适时相及特征参量，对黑龙江省主要农作物进行识别提取，并对3个年份的农作物空间分布图进行分析，为研究区农业遥感监测和作物种植结构调整等决策信息提供有效可行的参考和借鉴.

1 研究区概况

黑龙江省位于我国东北部，总面积约47.3万km2，东北和西南部(三江平原、松嫩平原等地区)地势较低(图1a)，有利于农作物的种植.以大陆性季风气候为主，春冬季低温少雨，夏季温热多雨(图1b)，≥10 ℃有效积温为1 600～3 300 ℃(图1c).第一积温带分布在西南部地区，热量充足，降水充沛，适合种植水稻等作物，而第六积温带纬度较高，水热条件不足，作物呈零散种植分布.

黑龙江省作为我国重要的商品粮基地，农业种植条件优越，全省耕地近160 000 km2，居全国第一位.无霜期140～180 d，主要农作物类型为水稻、玉米、大豆和大小麦等.2016年黑龙江省利用118 040 km2耕地面积生产了6. 057 8×107t粮食[20]，为全国提供了1/4的商品粮.

(a)数字高程图

(b)月平均最高、最低气温和降水量

(c)积温分布图

2 数据与方法

2.1 数据来源

从美国地质勘探局数据中心(USGS EROS，https://ladsweb.modaps.eosdis.nasa.gov/)获取的2001、2010和2017年46景8天合成MODIS地表反射率产品(MOD09Q1)，空间分辨率为250 m.利用MODIS产品批处理工具MRT(MODIS reprojection tool)，将投影统一转换成Albers等积投影，地理坐标系为WGS-84.以研究区范围为裁剪框对转换后的数据进行裁剪，并利用红光波段数据和近红外波段反射率数据计算归一化植被指数NDVI，获得3个年份的研究区MODIS-NDVI时间序列数据集.

辅助数据为黑龙江省2001、2010和2017年土地覆盖类型产品MCD12Q1，精度验证数据来源于黑龙江省历年农业统计年鉴.

2.2 研究方法

2.2.1 NDVI时间序列数据集

归一化植被指数NDVI是目前应用最广泛的植被指数[21].完整的NDVI时间序列数据集能够较好地反映农作物从播种—生长—成熟—收割的全过程.由于农作物的种类、生长特性和播种收割时间不同，表现出的光谱差异也不同.利用ENVI 5.2处理得到NDVI时间序列数据集.计算公式为

(1)

式中：ρnir为近红外波段的反射率，ρred为红波段的反射率.

2.2.2 MODIS-NDVI时间序列数据集重构

由于云、气溶胶以及地物的双向性反射等因素会对地表覆盖及陆地生态系统的监测产生极大的影响[3]，因此需要对植被指数时间序列集进行平滑处理，即重构时序植被指数.

目前常用的植被指数时间序列滤波去除噪声分析方法主要有Savizky-Golay滤波(S-G)、非对称高斯滤波(A-G)、双逻辑调和函数滤波(D-L)和傅里叶谐波(Hants)等.Jönsson等[22]的研究表明Hants滤波在重现植被物候特征时具有较好的拟合度,平滑时间序列数据集的效果较好.其计算公式为

(2)

式中：谐波的余项A0等于序列的平均值；Aj为各谐波的振幅；ω=2jπ/N为各谐波的频率；N为序列的长度；θj为各谐波的初相位；m为谐波的个数，m=N-1，这些正弦函数叠加构成傅里叶序列.

2.2.3 农作物遥感提取模型

在作物的整个生长周期内，不同作物的光谱曲线特征都有一定程度的差异性，根据不同类型作物独特的物候历特征，抓住不同作物在不同时间段生长周期不同的物候差异，可以实现作物种植类型提取和空间分布提取[23].

根据对黑龙江省主要农作物(玉米、水稻、大豆和大小麦)的物候历(表1)及植被指数与物候期对应关系的分析，利用MODIS-NDVI时序数据获取不同农作物在种植空间信息上的关键时间点，构建研究区农作物遥感提取模型(表2)，其中T1—T13为不同农作物在不同生长发育期内的NDVI阈值.

表1 黑龙江省主要农作物(玉米、水稻、大豆和大小麦)物候特征Tab.1 Phenological characteristics of main crops (corn, rice, soybean and wheat) in Heilongjiang Province

表2 黑龙江省玉米、水稻、大豆和大小麦遥感提取模型Tab.2 Remote sensing extraction model for corn, rice, soybean and wheat in Heilongjiang Province

以大豆为例进行说明.黑龙江省大豆种植时间大约在5月中上旬，5—6月出苗、开花，7月中下旬结荚，8—9月从鼓粒到成熟.随着大豆的出苗、开花、结荚、成熟和收割，反映在遥感影像上的NDVI值随着时间的变化先逐步增大，达到一个峰值后再逐渐减小.根据这一特性，大豆种植面积信息的识别要同时满足表2中序号为3的所有条件，为了保证结果的准确性，避免受到因各地种植习惯不同、气候干扰、成熟和收割时间差异带来的影响，模型中某一NDVI值的选择包含了前后两个时相的“或”.设置这些阈值参数并带入上述模型后，能对农作物的空间分布进行提取.

3 结果与分析

3.1 作物种植结构空间分布研究

根据表2作物空间分布遥感提取模型，得到黑龙江省主要农作物2001、2010和2017年空间分布图(图2).

2001年

2010年

2017年

图2 黑龙江省2001、2010和2017年主要农作物种植结构空间分布图Fig.2 Spatial distribution map of main crop planting structures in Heilongjiang Province in 2001, 2010 and 2017

由图2可见，自2001年开始，除了大小麦的种植面积不断减少外，其他3种作物种植面积不断扩张，且耕作范围有向三江平原逐步扩张的趋势.从行政区域上看，2001年作物种植主要分布在哈尔滨、齐齐哈尔和绥化，水稻种植区主要在南部，玉米和大豆主要分布于西南部和东北部；2010年，中东部的鹤岗、佳木斯和双鸭山种植面积也逐步扩大，多为玉米和大豆的混种模式；2017年，受到市场、政策等多方面因素的影响，种植结构进一步调整，大豆种植区主要分布于黑河和齐齐哈尔，佳木斯种植面积扩大且结构变化较为明显，作物类型以玉米、水稻为主.

3.2 作物种植结构时空变化分析

2001至2017年，黑龙江省主要农作物种植面积显著增长，种植结构发生了明显变化.为了探究种植结构调整的空间分布格局，将整个研究区主要作物种植空间信息分布图以市为行政单位进行划分，同时按经度和纬度两个方向进行分析，以获取其在时空变化上的规律.具体规则为：分市划分时将4种作物的种植面积总和计为单位“1”，统计每个市范围内每种作物所占的百分比(图3)；在经纬度上，每隔0.5°统计各作物种植面积(图4).

2001年

2010年

2017年

图3 黑龙江省主要农作物种植面积分市统计Fig.3 Statistical chart of Heilongjiang Province’s main crop planting surface integral city

由图3可知，除了大小麦的种植面积明显减少，玉米和水稻一直处于平稳增长的状态，而伊春和黑河两市大豆种植占比显著升高.玉米在近20年中的种植面积逐渐增大，并大幅向东部扩张；水稻面积主要变化集中在黑龙江省南部和三江平原地区，原因是那里水热条件好，积温充足，更有利于水稻种植；大豆种植区主要向北纬48°—51°、东经124°—126°迁移；大小麦的分布在经纬上都呈现出少量分布到零星分布的趋势(图4).

3.3 精度评价分析

由于本次进行作物提取选用的MODIS影像空间分辨率为250 m,不可避免地存在混合像元.但由于黑龙江省多以大面积成片种植为主，纯净像元较多，因此仅通过历年农业统计年鉴就可以对提取结果进行精度评价.表3显示，利用MODIS-NDVI时间序列数据集构建的黑龙江省主要农作物遥感提取模型最高精度为92.1%，最低精度77.0%，平均精度87.4%.

表3 黑龙江省主要农作物遥感提取模型精度验证Tab.3 Precision verification of remote sensing extraction model of main crops in Heilongjiang Province

主要从以下几个方面对误差精度进行评价分析：1)在构建黑龙江省主要农作物种植面积遥感提取模型中，阈值T1—T13主要基于各农场分布及选取的感兴趣区域(ROI，region of interest)设定，且模型选取的时间节点并非完全准确，有可能对提取结果产生影响.2)本次作物提取使用的MOD09Q1产品，时间分辨率为8 d，不能完全对应研究区农作物完整生长周期的关键时间节点，可能对作物面积信息提取的结果产生一定程度误差.3)本研究中使用的MODIS影像空间分辨率为250 m，即一个像元内存在像元“同物异谱”“同谱异物”的现象，对中小散户混合种植信息提取不够精确，从而影响作物的提取面积和精度.

4 结论

MODIS是具有中等空间分辨率、高时间分辨率且免费获取的影像数据，在大尺度区域作物信息识别和提取等方面具有明显优势.本文以2001、2010和2017年MODIS数据产品MOD09Q1为数据源，通过Hants滤波分析平滑MODIS-NDVI时间序列数据集，根据研究区主要农作物的种植区划、种植结构和物候历，分析玉米、水稻、大豆和大小麦的时间序列曲线特征，确定提取不同农作物的特征点和参量，构建农作物遥感提取模型，设定阈值，利用物候关键特征参数值对黑龙江省主要作物进行分类提取.结果显示：1)黑龙江省在近20年中主要农作物种植面积显著提高，且扩张范围主要在东北三江平原和中北部地区；2)种植结构方面，大小麦的种植面积持续减少，呈零散种植，玉米、水稻和大豆的种植面积持续增长，且玉米的增速最高；3)构建的农作物种植面积遥感提取模型的最高精度达到92.1%，最低精度77.0%，平均精度87.4%，表明通过MODIS影像可以有效对大范围耕作区域进行农作物信息提取，能够为区域农情遥感监测、作物种植结构调整等提供可行的参考和借鉴.

本研究对于大范围耕作区域农作物种植面积提取效果较好，但在模型构建、参数选择、多源数据融合等方面还存在缺陷，未来的研究方向主要有：1)由于作物是经过长期自然进化和人工选择后形成的[24]，而中国各地气候类型和地域性差异较大，农业生产利用面积广且季节性较强，散户、个体种植导致各地农业种植习惯和种植结构极其复杂，因此在作物面积提取前根据调查和统计信息及时更新不同地区的物候历和种植结构是农作物种植信息提取本地化的基础[14]；2)根据不同农作物的物候历特征及光谱特征，及时更新遥感提取模型中的关键数据阈值是保证模型精度的前提；3)由于MODIS影像为250 m中等空间分辨率，必然存在较多的混合像元，怎样在作物种植信息提取中考虑混合像元分解问题、如何与其他中高分辨率的卫星影像结合来提高提取精度都是未来研究的重点方向.