王涛，刘立波，张鹏*，王晓丽

1.宁夏大学信息工程学院，银川 750021

2.中国农业科学院农业信息研究所，北京 100081

3.国家农业科学数据中心，北京 100081

4.中国农业科学院国家南繁研究院，海南三亚 572024

引 言

估产数据集是表征农作物长势和产量的重要信息，被广泛用于作物长势分析、产量预测等研究。然而，传统的估产数据集制作方法采用人工实地测量土壤墒情、作物干重、气象条件等参数，从而构建数据集[1-3]。该方法速度慢、工作量大，且过度依赖于实测数据。

近年来，遥感影像因其低成本、覆盖面广等优势，常被用于提取归一化植被指数（Normalized Difference Vegetation Index，NDVI）和增强型植被指数（Enhanced Vegetation Index，EVI）。如 Quarmby[4-6]等利用红外波段和近红外波段计算得到了NDVI，Bolton等[7-8]利用红外、近红外和蓝光波段计算得到了EVI；此外，满卫东等[9]基于AVHRR影像制作了基于辽宁省1982至2009年AVHRR归一化植被指数数据集（https://www.osgeo.cn/data/wef57）。更进一步，Aerial Intelligence公司2017年发布的美国冬小麦产量预测数据集（https://aerialintel.blob.core.windows.net/recruiting/datasets/wheat-2013-supervised.csv 和 https://aerialintel.blob.core.windows.net/recruiting/datasets/wheat-2014-supervised.csv）不仅提取遥感影像中的NDVI和EVI，而且加入了气象、地理位置等特征，提升了估产特征的多样性。但在这类数据集中，NDVI、EVI等手工特征仅用少许的波段计算得到，忽略了其余波段重要信息，且特征的选取过度依赖人工经验，具有一定局限性。

因此，针对以上问题，本文摒弃手工特征制作的方法，制作了一种基于MODIS高光谱遥感影像的、多波段和多时相融合的宁夏县级尺度枸杞估产数据集，用于卷积神经网络特征的自动提取，简化特征提取操作的同时，进一步增强了特征丰富度。

1 数据采集和处理方法

1.1 数据来源

本文采用的实验数据由宁夏回族自治区遥感影像、枸杞种植区域矢量图和年际枸杞产量3类数据组成。其中，遥感影像采用 MODIS高光谱影像数据，来源于 EARTHDATA网站（https://search.earthdata.nasa.gov/），行列号为h26v04，时间范围为枸杞生长季内的每年第97天至第297天，包括MOD09A1、MOD13A1、MYD11A2、MCD15A2H 4种类型的MODIS产品数据（表1）；枸杞种植区域矢量图由宁夏农林科学院研究人员实地记录枸杞种植区域经纬度制成，为shp文件；年际枸杞产量数据来源于宁夏回族自治区统计局，包括2010-2019年宁夏16县（县级市区）枸杞种植面积和实际产量。其中，2019年产量数据如表2所示，枸杞种植总面积为27960公顷，总产量为94843吨，平均产量为3.39吨/公顷。其中，同心县、红寺堡区等地区枸杞产量均高于9000吨，属于高产地区；金凤区、平罗县等地区产量均在1000吨左右，属于中产地区；其余为低产地区。

表1 影像及矢量数据表Table 1 The table of remote sensing images and vector data

表2 2019年产量数据表Table 2 The table of yield data in 2019

1.2 数据处理

1.2.1 重投影与重采样

为了保证遥感影像数据空间位置的一致性，首先，利用 MODIS影像重投影工具（MODIS Reprojection Tool， MRT）工具将MODIS遥感影像和枸杞种植区域掩膜数据重投影为基于WGS-84椭球体的UTM投影；然后，将MYD11A2影像和枸杞矢量数据重采样为500 m，使MYD11A2影像和枸杞矢量数据与其余MODIS数据产品的空间分辨率相互统一。

1.2.2 时间序列补充

由于MOD13A1的时间分辨率为16天，其余MODIS数据产品为8天（表1），为了保证时间序列的完整性，根据式（1）采用上下影像求平均的方法对枸杞生长季内缺失影像进行补充，如将MOD13A1第97天和第113天影像的均值作为第105天的影像数据。

式中，Ii、Ii-8和Ii+8分别为MOD13A1第i天、第i-8天和第i+8天影像数据，i的取值范围为[105, 289]，时间间隔为8天。

1.2.3 波段融合和时间序列融合

为了提高波段信息的丰富度，分别提取了同一景MOD09A1、MOD13A1等遥感影像中的band1-band7、NDVI、EVI等13个波段（表1），并对其进行了融合。随后，为了进一步融合影像不同时相上的枸杞生长信息，对波段融合结果在时间维度上进行融合，最终形成了10张波段数为13，时相为26的时间序列影像。

1.2.4 枸杞种植区域提取

遥感图像是基于像素点的，包括了各种地面覆盖类型，为了使各县年鉴统计产量与其遥感影像相对应，利用了不同县域枸杞种植矢量图提取上述融合数据中的种植区域，以中宁县枸杞种植区域提取为例。首先，采用ArcGIS软件从宁夏行政区划矢量图提取出中宁县行政区划矢量图；然后，根据中宁县行政区划矢量图左上角和右下角经纬度裁剪得到中宁县融合影像（图1）；接着，调用GDAL（Geospatial Data Abstraction Library）库中的warp函数以实现枸杞中宁枸杞种植区域矢量图（图2）对图1的裁剪，从而得到中宁市枸杞种植区域图（图3）；最后，以此类推，提取出其余县枸杞种植区域。

1.2.5 直方图降维

鉴于影像数据集的稀疏性，无法采用端到端方式训练估产模型。因此，本文将影像无差别划分为32个像素区间，进而将影像中每个波段不同像素值映射至不同区间，以达到直方图降维的目的。

2 数据样本描述

基于MODIS的宁夏县级尺度枸杞估产数据集（2010-2019）主要由遥感影像数据和年鉴统计产量数据组成，具体数据样本描述如下。

图1 中宁县高光谱影像Figure 1 Hyperspectral image of Zhongning County

图2 中宁县枸杞种植区域矢量图Figure 2 Vector map of wolfberry planting area in Zhongning County

图3 中宁县枸杞种植区域图Figure 3 The area of interest for wolfberry planting in Zhongning County

2.1 遥感影像数据描述

遥感影像数据以MODI高光谱影像为数据源，经上述方法处理后，共形成了160个大小为32×26×13（32为像素区间数，26为时相数，13为波段数）的直方图矩阵，以“hist_年份_地名.npy”的方式进行命名，存放于“dataset/hist_data”文件夹中。其中，16个县（县级市区）（表2）中的每个地区各涵盖10年直方图矩阵数据，时间范围为2010-2019年，直方图矩阵数据示意如图4所示。

图4 直方图矩阵数据示意图Figure 4 Schematic diagram of histogram matrix data

2.2 统计产量数据描述

为了方便用户对数据的加载及使用，本文提供了与直方图矩阵年份和地名相对应一致的年鉴统计产量数据。该数据为表格数据，保存于“dataset/yield_data/ yield_data.csv”文件中，共包含7个属性列，分别为年份、省份名、省份编号、县名、县级编号、枸杞种植面积和产量，总计160条产量数据，平均每个县（县级市区）包含10条2010-2019年的枸杞产量数据，部分产量数据如表3所示。

表3 部分产量数据表Table 3 The table of partial yield data

3 数据质量控制和评估

为了验证本文数据集的可用性，本文分别从定性和定量两个角度评价本数据质量和准确性，具体如下：

3.1 定性评价

通过可视化对比高中低产量所对应的直方图降维结果（图5）可见，在高产、中产和低产中波段1、波段7、NDVI和EVI波段明显存在视觉差异，表明可以从本数据集中提取到有利的特征，用于表征枸杞的长势和产量。此外，在时间维度上，像素区间最大值基本集中于第180天左右，正值枸杞植株长势茂盛期，符合枸杞的生长规律，说明本文数据集质量较好。

图5 可视化结果图Figure 5 Visualized result graph

3.2 定量评价

为了进一步说明该数据集的准确性，本文以平均相对误差（Mean Relative Error, MRE）、均方根误差（Root Mean Square Error, RMSE）和决定系数（Coefficient of Determination, R2）为评价指标，计算公式如式（2）-式（4）所示，分别基于本文数据集和Aerial Intelligence公司发布的美国冬小麦产量预测数据集，设置了两个对比实验组，分别是：

（1）以本文数据集为数据源，构建了卷积神经网络（Convolutional Neural Networks, CNN），用于枸杞长势特征的抽取，进而采用全连接网络（Fully Connected Network, FCN）、支持向量回归（Support Vector Regression, SVR）和岭回归（Ridge Regression, RR）方法预测枸杞年际产量。卷积神经网络结构如图6所示，第一层为输入层，为本文数据集，训练集、验证集和测试集划分比率为7:2:1；第二层为特征提取层，由 6个卷积层构成，卷积核数量分别是 128、256、256、512、512、512，卷积核大小全为3×3；最后一层为全连接层，将特征图映射至2048维embedding（嵌入向量）空间，以表征枸杞的长势信息。

（2）与上述方法相比，对比实验基于冬小麦产量预测数据集，仅用单一的FCN、RR和SVR预测冬小麦产量。其中，FCN包含5个隐藏层，维度分别为1024、512、256、128和1，实验结果如表4所示。

表4 对比实验结果表Table 4 The table of comparative experiment results

由表4可知，在本文数据集上，MRE和RMSE分别为14.52%、859.23吨，且R2达到了0.83，均优于对比数据集。与2017年冬小麦产量预测比赛数据集相比，本文数据集采用波段融合和时间序列融合方法分别融合了不同波段和时相上的影像信息，增强了特征丰富度，所以基于本文数据集的CNN回归方法，在 MRE和 RMSE上分别下降了1.48%和89.14吨，R2上升了0.06，说明增加了0.06%的特征可由回归方法解释，验证了本文数据集的准确性。

其中， Ti为第i个县枸杞统计产量，Pi为第i个县产量预测值，为16县统计数据平均产量，单位为吨。

图6 卷积神经网络结构图Figure 6 Convolutional neural network structure diagram

4 数据价值

目前，公开的遥感估产数据集非常有限，且目视解译纹理特征等相关提法过度依赖于实测数据和人工经验。因此，本文构建了一种基于MODIS的宁夏县级尺度枸杞估产数据集，具有高光谱、多波段、多时相等特点，可用于卷积神经网络特征自动提取农作物长势特征，为作物产量预测研究提供数据支撑，同时为遥感多时相估产影像数据集的制作提供了一定的参考价值。