基于ANN多源遥感数据融合的耕地种植强度 估算方法研究

徐 猛，陶建斌，吴琪凡

（地理过程分析与模拟湖北省重点实验室/华中师范大学城市与环境科学学院，湖北武汉430079）

0 引言

农业生产是人类社会存在和发展的基石，耕地是进行农业生产的基本资源和条件。随着我国城市化迅速发展，每年约有20万hm2耕地被转化为建设用地［1-3］。同时，随着人口持续增长、饮食结构改变和非食物用途（如生物能源）农产品消耗的增加，人类对农产品的需求不断增长。2015年末全国人口为13.75亿，预计到2030年全国人口将增加至14.5亿［4］，增长5.45%，而同时期粮食需求将增长47.19%，达到9.17亿t［5］。如何缓解日益加剧的人地矛盾，是国际社会普遍关注的焦点问题。在耕地面积减少、粮食需求持续增长和部分耕地质量恶化等多重压力下［6］，耕地集约利用成为保障国家粮食安全的重要选项。同时，由于城市化快速发展对农业劳动力的虹吸效应，以及市场经济的发展对农民种粮收益的冲击，部分地区农作物种植频率呈现下降趋势（表现为耕地撂荒和冬闲田）。因此，耕地的集约利用成为解决未来世界粮食安全的重要途径，引起学界和业界的广泛关注。

当前，国内外研究者多数将耕地的集约利用定义为种植频率（单季和双季等）［7-13］或者复种指数［14-20］，也有部分研究者从投入产出、生产管理和技术角度来对耕地利用强度进行定义［21-27］。Kehoe等［23］从输入、输出、系统3个测度总结了现有的对Agricultural Land-use Intensity（简称ALI）的定义。Jiang等［24］认为ALI通常有3种测度：种植频率、农业产出、农业投入。然而，无论是种植频率还是复种指数，都是对种植模式的硬性划分（如单季或双季）。一方面，尽管各种复种指数定义不尽相同，但本质上都是种植频率在空间尺度上的聚合或是时间尺度上的综合。另一方面，基于统计数据得到的一系列耕地指标缺乏详细的空间分布信息，忽略了行政区内部耕地集约利用水平的空间异质性。因此现有评价体系存在时、空分辨率上的矛盾，不能满足耕地集约化利用的精细评价的要求。

在国外，很早就有了Cropping Intensity或ALI的概念［28-30］，但却有许多不同的定义。在对Cropping Intensity进行系统梳理的基础上，文章定义种植强度指数对耕地集约化利用进行精确表征，通过构建人工神经网络（Artificial Neural Network，ANN）对时间序列MODIS数据与种植强度间的映射关系进行表达，以湖北省为研究区进行种植强度估算。该文研究方法可获取高时空分辨率的种植强度数据集，对推动作物制图的进一步发展，提高耕地集约化精细监测具有重要意义，为全球变化过程的分析与模拟提供理论、方法和数据支撑。

1 研究区及数据

1.1 研究区

湖北省地处我国中部（东经108°21′42″～116°07′50″、北纬29°01′53″～33°6′47″），面积18.59万km2，约占全国国土面积的1.97%，耕地面积占到全国耕地面积的3.89%。湖北省东、西、北三面环山，中部低平。除部分山区外，大部分地区属亚热带季风气候，光热充足，降水丰沛。优越的自然环境和悠久的种植传统使湖北省成为著名的鱼米之乡，其粮油生产在全国占有重要地位。

1.2 数据来源及数据预处理

该文采用的主要数据包括：（1）MOD13Q1（MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid）250 m分辨率16 d合成的植被指数数据［31］；（2）Landsat 8遥感影像［32］。另外还采用了其他数据作辅助，包括Google Earth历史影像，GlobeLand30 —2010全球30 m空间分辨率土地覆盖数据［33］，中国1∶400万矢量数据集。

MOD13Q1是陆地植被指数数据，具有250 m的空间分辨率和16 d的时间分辨率，在地表植被监测领域的应用十分广泛。MOD13Q1数据包括增强型植被指数（Enhanced Vegetation Index，EVI）、归一化植被指数（Normal Differential Vegetation Index，NDVI）、像素质量图层（Pixel Reliability，PR）等12个波段。该文采用2015年MOD13Q1数据，研究区所涉及图幅编号为h27v05、h27v06、h28v05、h28v06。高分辨率遥感数据采用了3期Landsat 8多光谱遥感数据（表1），空间分辨率为30 m，重返周期为16 d。该研究中，MOD13Q1数据用于神经网络的输入特征，Landsat 8遥感数据主要用于训练样本和验证样本的制作。

图1 研究区地理区位及地物覆盖类型（GlobeLand30—2010）Fig.1 The location of the study area and the land-cover types（GlobeLand30 —2010）

表1 Landsat 8实验数据详细信息Table 1 Landsat 8 experimental data details

数据预处理包括3个部分：（1）利用MODIS批量处理工具（MRT）对MOD13Q1数据进行预处理，重采样至240 m，投影坐标选择WGS_1984_UTM_Zone_49N，提取出EVI图层，并以Landsat 8影像为基准对MOD13Q1数据进行配准。（2）MOD13Q1数据受到云、雪、阴影等诸多因素的影响，原始EVI数据存在很多噪声，需要通过数学方法对原始EVI数据进行重建。TIMESAT是用于平滑长时间序列EVI数据的有效工具［34］，提供多种数据重建方法，该文采用Savitsky-Golay算法（简称S-G滤波算法），这是一种基于最小二乘法的移动窗口加权平均的数据重建方法［35］。S-G滤波算法在去除噪声的同时，能够清晰的保留局部细节，并且对EVI曲线的峰值和峰宽保真度较高［36］。（3）提取耕地掩膜，对GlobalLand 30 —2010年数据产品进行投影变换、裁剪、重分类等处理，提取湖北省的耕地掩膜。

2 研究方法

利用神经网络估算耕地种植强度的基本假设是：多时相遥感影像的EVI数据反映了耕地的利用状态，耕地的不同利用状态与耕地种植强度之间存在映射关系。该文将MOD13Q1时间序列EVI数据和基于Landsat 8数据获取的样本区种植强度分别作为神经网络的输入和输出层，训练网络模型，利用该网络模型估算整体耕地种植强度。

2.1 种植强度指数

种植强度指数是指在亚像元尺度下综合表征农作物的空间分布和种植频率的指数，其公式为：

式（1）中，CII为种植强度指数，Dcrop为1个MODIS像元内Landsat 8遥感影像上双季作物的像元个数，Scrop为1个MODIS像元内Landsat 8遥感影像上单季作物的像元个数，n为1个MODIS像元内Landsat 8遥感影像上的像元总数。

CII为MODIS像元内的种植频率的均值，值域范围为0～1（图2）。在地块完整的双季作物种植区，其种植强度为1或接近于1。在地块完整的单季作物种植区，其种植强度为0.5或接近0.5。非耕地区域的种植强度则为0（忽略研究区内少量的三季作物）。

图2 耕地种植强度示意图Fig.2 Diagram of cropping intensity

2.2 样本种植强度

Landsat 8遥感影像用于样本区耕地种植频率的分类。在种植频率分类结果的基础上采用ArcGIS的分区统计与重采样工具实现种植频率的概化工作。

湖北省耕地种植模式主要以单季和双季为主［37］。单双季作物的物候期存在较大差异，单季作物生长起始期较晚，大概在每年4月初。双季作物第1个生长季约在每年1月中下旬开始，在4月中下旬结束，EVI曲线在第129、145 d处明显下降（图3）。以2015年为例，在1—4月的Landsat 8真彩色合成图像上，单季作物种植区呈现暗灰色，双季作物种植区呈现暗绿色（图4）。选取MOD13Q1的001、145这2期EVI影像，将001、145、001分别赋予红、绿、蓝通道，能够明显地凸显出双季作物种植区，根据加色法原理，双季作物种植区呈现出品红色（图5）。借此结合Google Earth高清影像以辅助判读Landsat 8遥感影像的非监督分类结果。样本区的选择应包括耕地、水体、天然植被、建设用地、裸地等多种地物类型，且面积不应小于研究区面积的15%，以增强模型的鲁棒性和泛化性。该文采用非监督分类ISODATA分类方法对样本区Landsat 8遥感影像进行分类，初始类数量为20类，迭代次数为15次。然后将初始分类结果重分类为双季作物、单季作物和其它3类。借助ArcGIS创建240 m的格网，分区统计每个网格单元内种植频率的平均值，得到每个网格单元的种植强度，值域为0～1（图6）。

图3 湖北省耕地单季作物与双季作物时间序列EVI曲线Fig.3 EVI profiles of single cropping crops and double cropping crops of cropland in Hubei Province

图4 2015年湖北省耕地单季作物与双季作物典型分布区Fig.4 Typical distribution areas of single cropping crops and double cropping crops of cropland in Hubei Province in 2015

图5 2015年湖北省耕地MOD13Q1 EVI假彩色合成Fig.5 False color composite of MOD13Q1 EVI of cropland in Hubei Province in 2015

图6 基于Landsat 8影像的训练样本区数据：a. 种植频率；b. 耕地种植强度Fig.6 Cropping frequency and intensity of training sample area based on Landsat 8

2.3 模型训练

人工神经网络是一种模拟生物神经网络，进行分布式并行信息处理的数学模型，它具有自主学习的能力，通过不断改变神经节点之间的连接权重，完成由输入数据到目标数据的复杂映射［38］。BP神经网络又称为前馈型神经网络，是人工神经网络中应用最为广泛的模型之一，其由3部分组成：输入层、隐含层（可以是单层，也可以是多层）、输出层。BP神经网络通过误差的反向传播，不断调整输入层与隐含层、隐含层与输出层之间的权重和误差，使网络输出与实际输出的误差最小化［39］。该文利用MATLAB设计BP神经网络自主学习样本区23期EVI数据集与耕地种植强度间的映射关系，样本区23期EVI作为输入层，耕地种植强度作为输出层。在隐含层数量及隐含层神经元个数的选择上，该文参考Kolomogorov定理。Hecht-Nielsen指出，映射网络存在Kolomogorov定理，即对于任意连续函数f ：［0，1］n→Rm都可以用1个3层前向神经网络实现，且中间层神经元个数不多于2n+1个［40］。因此该文设计的BP神经网络为3层：输入层、输出层和隐含层。关于神经元个数的选择，在遵循Kolomogorov定理的基础上，通过多次试验，最终将神经元个数定为15，即构成1个23-15-1的BP神经网络模型，见图7。神经网络其他参数设置如下：训练函数TRAINLM，学习函数LEARNGDM，隐含层与输出层传递函数均为LOGSIG。

图7 BP神经网络模型Fig.7 Back propagation neural network model

3 结果与分析

3.1 样本精度验证

基于神经网络估算的耕地种植强度其准确性受样本精度影响，为了利用验证样本对估算结果进行精度评价，需要保证样本精度。该文从Google Earth历史影像（2015年3月6日、2015年3月25日、2015年3月27日、2015年10月11日）上选取了其他、单季作物、双季作物3类共计1 218个真实样本（图8），以验证样本精度。结果表明，基于Landsat 8影像得到的样本总体精度达到了93.51%，Kappa系数为0.914（表2），证明了样本具有很高的可信度。

图8 2015年湖北省耕地真实样本空间分布Fig.8 Spatial distribution of truth sample of cropland in Hubei Province in 2015

表2 样本分类结果混淆矩阵Table 2 Confusion matrix of sample classification results

3.2 种植强度精度验证

该文通过种植强度指数和人工神经网络估算方法对验证样本区进行处理，制作验证样本区耕地种植强度，并与BP神经网络提取的种植强度结果作对比。从目视结果来看，Landsat 8遥感影像上的种植频率分布与BP神经网络提取的耕地种植强度呈现出较高的一致性（图9）。

利用回归分析验证BP神经网络提取的湖北省耕地种植强度的精度。图10显示了验证样本区基于Landsat 8影像的耕地种植强度与BP神经网络提取的耕地种植强度的散点图和线性关系。R2达到0.923，回归方程系数为0.948，截距项为0.023，说明二者之间的数据偏差很小，验证了BP神经网络提取耕地种植强度的可靠性。

图9 a. 基于Landsat 8影像分类结果的验证样本区种植频率 b. 基于BP神经网络提取的验证样本区耕地种植强度Fig.9 a. Cropping frequency in validation sample area b. Cropping intensity in validation sample area extracted from BP neural network

图10 验证样本区耕地种植强度与BP神经网络耕地种植强度散点图Fig.10 Scatter plot of cropping intensity in validation sample area extracted from Landsat 8 and cropping intensity in validation sample area extracted from BP neural network

3.3 湖北省耕地种植强度分析

该文利用神经网络提取了2015年湖北省耕地种植强度（图11）。基于BP神经网络提取的种植强度实现了对耕地集约利用程度的软划分。BP神经网络提取的种植强度与陶建斌等［41］提取的2015年湖北省种植频率分布基本吻合（图12）。图13为陶建斌提取的2015年湖北省耕地种植频率与该文BP神经网络提取的耕地种植强度在江汉平原地区的对比，结果显示BP神经网络提取的耕地种植强度空间分辨率更高，细节表现能力更强，层次更加分明。

图11 2015年湖北省耕地种植强度Fig.11 Cropping intensity in Hubei Province in 2015

图12 陶建斌等提取的2015年湖北省耕地种植频率［41］Fig.12 Cropping frequency from Tao Jianbin，et al in Hubei province in 2015

基于统计数据［42］计算得到2015年湖北省县级尺度的耕地复种指数（图14）。BP神经网络提取的耕地种植强度与基于统计数据得到的耕地复种指数间总体上较为吻合，复种指数较高的江汉平原地区对应着种植强度高值区，复种指数较低的西部山区对应着种植强度的低值区。然而相较之下，BP神经网络提取的耕地种植强度分辨率更高，对耕地集约利用水平的评价更为精确。

图13 2015年湖北省耕地种植频率与种植强度对比a. 该文提取的种植频率 b. 陶建斌等提取的种植强度Fig.13 Comparison of cropping frequency and cropping intensity in Hubei province in 2015

图14 2015年湖北省耕地复种指数Fig.14 Multi-cropping index in Hubei Province in 2015

4 结论与讨论

该文提出了种植强度的概念，并定义了种植强度指数以精确表征耕地的集约化利用。构建BP神经网络模拟时间序列EVI数据集与种植强度间的映射关系。估算结果精度达到92.3%，证明了BP神经网络提取耕地种植强度的可靠性。该方法充分发挥遥感技术多尺度的时空表达能力和地物时间谱特征捕获能力的优势，并结合人工神经网络自主学习能力的特点，可获得高时空融合的耕地种植强度数据集。

传统的基于统计数据的耕地集约化研究，尽管在时间尺度上能够保证良好的连续性，但基于行政区单元的统计数据往往忽略了内部的空间异质性，空间分辨率较低；而国内外现有的基于遥感数据的耕地集约化利用研究虽然保证了较高的空间精度，但在方法上大都采用了连续年份的平均；该文提出的耕地种植强度弥补了现有指标的不足，将耕地集约利用研究尺度缩小至亚像元级别，具有融合高时间分辨率和高空间分辨率的特点，能够实现对耕地种植强度的精细监测。

该文的研究重点在于方法的探究，未来将关注该方法应用于更大区域尺度耕地种植强度的时空变化以及种植强度估算中的地域分异问题。同时，该文在使用BP神经网络时未探究其作用机理，导致网络模型训练的不确定性缺乏合理解释。