京津冀地区冬小麦播种面积快速提取技术研究

吴风华 郎婷婷 江东 丁方宇 张紫萍

摘要：冬小麦是我国的主导粮食作物之一，及时掌握冬小麦种植信息对农业部门制定政策、调整农业结构具有十分重要的意义。以京津冀地区为例，综合考虑遥感数据源成本、数据处理时间以及时间分辨率等要素，提出了一种以MODIS-NDVI产品为主要数据源，提取冬小麦播种面积的技术方法。通过对冬小麦种植区的实地考察与采样，准确获取了反映冬小麦生长周期的NDVI时间序列曲线，并将采样点NDVI时间序列曲线与物候期相结合，识别冬小麦的典型物候期，以决策树分类方法设定阈值，逐级快速剔除了非耕地区和林地区，实现了冬小麦种植区的提取。提取的冬小麦种植区与实地采样点的匹配度达到93.33%，与统计年鉴公布数据的比较结果表明分类精度达到 91.84%，为快速提取大面积农作物种植信息提供了一种可操作的技术手段。

关键词：京津冀;MODIS-NDVI;冬小麦;播种面积;决策树分类;面积提取

中图分类号： S127  文献标志码： A

文章编号：1002-1302（2019）16-0224-06

收稿日期：2018-04-19

基金项目：国家重点研发计划（编号：2016YFC0503507）;西藏自治区重大科技专项（编号：Z2016C01G01）。

作者简介：吴风华（1972—），女，湖南宁乡人，硕士，副教授，主要从事地理信息系统的教学与研究工作。

通信作者：江 東，博士，研究员，研究方向为资源环境遥感监测。

农作物种植面积信息的准确获取是提高作物产量预测和确保粮食安全问题的关键，为政府相关管理部门和机构提供了重要的农业生产依据。遥感技术具有时效性高、覆盖率广、观测周期短等优点，在农业监测体系中逐渐取代传统的人工核查方法，推动了农业现代化发展[1]。

不同的遥感影像数据具有不同的特点，高空间分辨率影像虽然能够较好地识别地物的纹理、形状等特征[2]，但其单景影像幅宽、时间分辨率有一定的局限性，加之云、大气等噪声污染，难以获得覆盖整个研究区的影像数据[3]。然而，基于中低分辨率的MODIS数据具有覆盖范围广，数据容易获取，以及空间、光谱、时间分辨率独特组合等特点，在大尺度面积测算中到了广泛的应用[4-5]。基于MODIS数据的遥感技术已广泛应用于农作物监测中。许青云等为获取陕西省农作物种植模式和农作物空间分布信息，以MODIS遥感影像为数据源，获取多年NDVI时间序列数据，结合农作物物候信息，识别了农作物种植模式并获得了多种农作物空间分布信息[6];郭昱杉等以黄河三角洲为研究区，选取MODIS数据，利用Hants滤波重建NDVI时间序列曲线，实现农作物空间信息分布提取，空间信息分布的总体匹配精度达86.9%[7];Pan等以MODIS遥感影像为数据源，在农作物种植面积估算中，考虑到混合像元在面积提取中带来的误差，便引入农作物物候指数比例（CPPI），以通州和沭阳为研究区，构建了植被指数（EVI）与冬小麦播种面积的表达式，将单个像元的均方根误差由15%降到了5%[8]。

基于遥感技术的农作物分类方法目前有多种，支持向量机（SVM）、神经网络、构建决策树等分类方法都已被应用到遥感分类中，如平跃鹏等利用平滑后的NDVI时间序列数据以支持向量机（SVM）对农用地逐级分层分类，得到农作物分布信息[9];李祚泳基于已知类别的遥感图像样本集的先验知识，利用B-P神经网络方法对样本反复训练，实现了对4种农作物的分类[10];潘琛等以ETM+多光谱影像作为数据源，将NDVI、GVI等10种植被指数应用于影像分类中，构建分类决策树，实现了对江苏省徐州市景观格局的分类，并且总体分类精度有所提高[11];刘焕军等对覆盖黑龙江省虎林市的影像构建的NDVI时间序列曲线，建立分类决策树，成功提取了研究区范围内的土地覆被信息，且分类精度均高于最大似然法、马氏距离法等分类方法[12];本研究拟采用决策树分类法实现对京津冀地区冬小麦播种面积的提取，此方法简单易行，便于理解，能处理多种属性数据的分类，在确定了各个属性的分类标准后生成分类规则，完成地物分类。

京津冀地区作为中国人口聚集的“首都圈”，承载1亿多人口的粮食需求，且河北省是中国冬小麦主产区之一[13]，准确获取京津冀地区的主要农作物种植面积对保障粮食安全和经济稳定具有非常重要的意义[14]。本研究以MODIS-NDVI为数据源，辅以土地利用图，结合实地采样数据，通过分析样本数据的NDVI时间序列曲线，并与农作物物候特征结合，建立决策树设定阈值，快速提取了京津冀地区冬小麦的播种面积，以现场调查数据和统计数据对提取结果进行了校验。

1 材料与方法

1.1 研究区概况

京津冀位于我国华北平原中部，地理位置为113°27′～119°50′E，36°05′～42°40′N，区域气候属于温带湿润半干旱大陆性季风气候，月均气温在-12～32 ℃之间，降水主要集中在6—9月，日降水量最大可达360 mm。京津冀总面积为 21.6万km2，其中耕地面积为6.54万km2，农作物种植制度以一年两熟为主，兼有一年一熟和两年三熟等多种种植形式，农作物大致可分为谷物、豆类、薯类、棉花，其中谷物主要为小麦、玉米、稻谷，粮食作物在该地区占有重要地位。

1.2 总体技术路线

本研究的技术路线（图1）分为4个主要步骤：（1）野外采样调研。根据土地利用现状图提供的耕地区，制定野外采样调研路线，在野外采样时，以Landsat8 OLI影像作为采样底图，将GPS定位仪连接到ArcGIS 10.1中，记录典型地物的坐标点。（2）NDVI时间序列数据处理。将下载的MOD09Q1影像进行预处理，并对处理后得到的波段影像进行计算，得每幅影像的NDVI数据，利用研究区边界对影像进行裁剪，得到京津冀地区范围内的NDVI数据，通过S-G滤波重构NDVI时间序列数据。（3）建模。对一部分野外采样点的NDVI时间序列数据进行分析，提取典型物候期，通过建立决策树进行分类，获得京津冀地区冬小麦播种面积。（4）精度验证。利用统计年鉴公布的数据对提取的冬小麦播种面积进行定量分析，另一部分野外采样点作为冬小麦种植区的空间验证。

1.3 数据处理

1.3.1 MODIS-NDVI产品预处理 为获取冬小麦播种面积，选取由NASA网站（http：//modis.gsfc.nasa.gov）提供的MOD09Q1数据产品，时间为2016年第313天至2017年第313天，时间分辨率为8 d，共计47个时相，空间分辨率为 250 m，覆盖研究区行列号为h26v04、h26v05、h27v04、h27v05。首先使用MODIS产品批处理工具MRT对下载的影像进行拼接、重采样、投影，统一在WGS-84地理坐标系下，采用Albers进行等积投影转换。批处理后的影像在ArcGIS 10.1平台下以京津冀边界矢量图进行掩膜处理，得到研究区范围内的影像。经过拼接、投影、裁剪等预处理后的MODIS影像包含近红外波段（b2）和红波段（b1），在R语言中，对预处理后的图像进行波段计算，得到每个时相的NDVI数据（图2）。

1.3.2 辅助数据 本研究所用的辅助数据为由中科院地理所资源利用与环境修复重点实验室提供的2015年土地利用现状图、京津冀边界矢量图，及在美国USGS网站（http：//glovis.usgs.gov/）下载的研究区范围内的Landsat8 OLI影像。其中，京津冀边界矢量图作为MODIS影像的裁剪范围，土地利用现状图作为农作物种植信息分布图，结合Landsat8 OLI影像来寻找研究区内冬小麦种植区，为样本采集区的选定做准备。将统计年鉴中的冬小麦播种面积数据作为冬小麦播种面积提取精度的验证。

1.3.3 采样数据 以2015年土地利用现状图作为农作物种植分布参考图，选定河北省中部及中南部的冬小麦种植集中区作为样本采集区。对下载好的Landsat8 OLI影像进行预处理，将4、3、2波段融合，便于目视解译地物，将融合后的影像作为采样底图加载到ArcGIS 10.1中，同时将GPS定位仪连接到ArcGIS 10.1中，开始定位采样。由于温度变化、冬小麦播种品种、起始种植时间等因素的影响，导致冬小麦物候期在地理位置分布上存在差异[15]，为了得到准确可靠的冬小麦物候期NDVI特征值来进行冬小麦播种面积提取，本研究选定的冬小麦样本区起于北京南部，贯穿河北保定、石家庄，至邢台北部，覆盖河北省6个县、16个村，采集冬小麦样本个数为215个，均匀分布于冬小麦种植区，另采集到28个建筑用地样本，19个林地样本，樣本分布如图3所示。

2 研究方法

2.1 NDVI时间序列重构

由于MODIS数据受外界条件的干扰，导致数据存在严重的异常值，为了削弱异常值，对原始的NDVI时序曲线进行 S-G滤波平滑，以便能较好地反映研究区内各类地物的NDVI时序曲线走向。S-G滤波是由Savitzky和Golay提出的一种简化的最小二乘卷积法，用来平滑和计算一组连续值或者光谱值导数，卷积法基于多项式给定的权重进行滑动平均滤波[16]，即S-G滤波是一种权重滑动平均滤波，权重由一个滤波窗口内做最小二乘拟合的多项式次数来确定。S-G滤波基本公式：

Y*j=∑i=mi=-mCiYj+1N。

式中：Y表示NDVI原始值，Y*表示NDVI拟合值，j表示原始NDVI数组的系数，Ci表示第i个NDVI值滤波时的系数，m指窗口的宽度，N指卷积数目，等于滑动数组宽度（2m+1），滑动数组包含2m+1个点。根据NDVI观测值调整2个参数，第1个是滑动窗口大小m，通常情况，m越大曲线越平滑，但可能将峰值过于平滑;第2个参数是平滑多项式次数d，通常在2～4范围内设置，d值越小，曲线越平滑，但可能引入偏差，d值较高时，减少过滤偏差，但去噪效果低。

S-G滤波平滑适用于连续、等距的一系列数据，MODIS-NDVI数据符合这一特征，该滤波法能够有效去除遥感影像中存在的云、大气等噪声，清晰地描述NDVI时间序列的变化趋势特征和局部突变信息，理论简单易于实现[17]。本研究使用R语言对47时相的MODIS-NDVI进行平滑滤波，结果如图4所示。

2.2 冬小麦物候特征提取

不同农作物都有不同的物候特征，物候特征被作为区分绿色植被的重要依据之一[18]。华北平原冬小麦播种时期大概在9—10月，物候期可大致分为播种期、出苗期、休眠期、返青期、开花期、成熟期，华北平原1981—2009年的年均物候期及标准差如表1所示[19]。

为了提取冬小麦与其他地物有明显差异的物候期，选取部分采样点作为样本参考点，其中，冬小麦分布范围大于 250 m×250 m的点定义为纯净冬小麦样本参考点，居民地及道路定义为建筑用地样本参考点，由于实地采样中的林地范围较小，定义林地分布范围最接近250 m×250 m的点为林地样本参考点。本研究分别选取了20个冬小麦采样点、10个建筑用地采样点及5个林地采样点作为样本参考点进行分析（图5）。

本研究结合冬小麦物候期，对经S-G滤波后样本点的NDVI时间序列参考曲线进行分析（图5）：该地区冬小麦在第281天开始播种，第289天后出苗，此时为9月底10月初，树木趋于衰落，NDVI值下降，冬小麦生长，NDVI值上升;随气温逐渐降低直到休眠期前，冬小麦NDVI值达到一个小波峰，而林地树叶凋落，NDVI持续最低;在进入返青期后（3月初），冬小麦迅速生长，NDVI值不断上升，开花期前后一直处于波峰状态，期间树木生长，NDVI值处于上升状态，但NDVI值一直低于冬小麦;在冬小麦到达成熟期（5月底），NDVI值逐渐跌入波谷，直至收割阶段NDVI值达到最低，由于树木一直处于生长状态，其NDVI值逐渐上升，此阶段两者NDVI值有明显的区别。建筑用地NDVI值一直处于低值状态，由于建筑用地内种有树木花草等绿植，受其影响，在8月份有一个波峰状态，但整体NDVI值都偏低。为了对冬小麦、林地、建筑用地进行分类，提取NDVI值有明显差异的时期，选取开花期（第124天）、休眠期前（第336天）、成熟期后（第161天）作为提取冬小麦的典型物候期。

2.3 基于决策树分类提取冬小麦

决策树分类是一种基于空间数据挖掘和知识发现的分类法，通过对训练样本进行归纳学习生成决策规则，实现对数据的分类，分类样本不需要满足正态分布，但要严格“非參”[20-21]。决策树的基本组成为根节点、节点、分支和叶节点，在分类过程中更直观、更清晰，便于理解，能够提取最有价值的特征参数，且计算效率高，常被应用于遥感影像分类中[22]。根据上述NDVI时间序列参考曲线分析提取的典型冬小麦物候期，结合图5设定阈值提取冬小麦种植区：

（1）在第124天，根据冬小麦样本点的NDVI最大值和最小值设定阈值，区分绿色植被与建筑用地区，提取含有冬小麦的绿色植被区;

（2）在第336天，根据冬小麦样本点的NDVI最小值设定阈值，区分冬小麦与其他绿色植被，高于设定的阈值，则为含有冬小麦的区域;

（3）在第161天，根据冬小麦样本点的NDVI最大值设定阈值，区分冬小麦丰度高低的区域，低于设定的阈值，为冬小麦丰度高的区域，高于设定的阈值，则为冬小麦丰度低的区域，在此研究中被视为林地，不被作为冬小麦提取。在NDVI时间序列参考曲线中，1时相到27时相期间（如图5-a和图5-b所示），林地的一条最大NDVI样本曲线与冬小麦的最小NDVI样本曲线趋势基本一致，此后便有较大差异，可能是两者互相影响，出现混合像元，造成NDVI曲线趋势相似，也可能是林地套种冬小麦（图3），为了达到快速分类提取的目的，统一将其视为林地。

3 结果与分析

3.1 冬小麦面积提取结果及验证

利用MODIS-NDVI数据，基于决策树分类方法提取了含有冬小麦的像元，像元分布如图6，研究区内冬小麦播种面积提取结果为220.587 6万hm2。

本研究精度验证从2个方面进行分析——定量分析和空间分析。

定量分析是将研究区范围内冬小麦提取面积与统计年鉴数据报告进行对比，统计年鉴中研究区内冬小麦播种面积为240.179 2万hm2（由于北京、天津统计年鉴数据未更新至2016年，且天津市统计年鉴的小麦面积未将冬小麦与春小麦分离，故将2012—2014年冬小麦平均播种面积作为参考面积），实际提取冬小麦播种面积为220.587 6万hm2，精度达到91.84%，各市的提取面积及精度如表2所示。

空间分析是通过野外采样点来进行验证。冬小麦野外采样点为215个，其中20个采样点已作为冬小麦样本参考点，剩余195个采样点作为精度验证数据。分类结果表明，182个点被正确分类，空间匹配精度达到93.33%。

3.2 结果分析

图6及表2结果显示，京津冀地区的冬小麦主要分布在河北省东南部，即河北省沧州市、邯郸市、衡水市、石家庄市、邢台市、保定市，且各市冬小麦面积均在30万hm2左右。表2结果显示，基于各市的冬小麦提取精度高低各不相同，这一现象是由低分辨率MODIS影像中的混合像元存在的结果造成的。在实地采样调研中，冬小麦种植区有小范围的林地和套种地存在（图3），范围为50 m×100 m、100 m×200 m不等，小范围林地的存在，使得混合像元NDVI时序曲线表现得和冬小麦NDVI时序曲线相似，导致在提取覆盖冬小麦像元的同时，也可能提取了含有林地的像元，使得提取面积大于实际播种面积;在冬小麦种植不集中的区域，由于建筑或者裸地的存在，使得含有冬小麦和建筑的混合像元的NDVI时序曲线和纯净冬小麦像元NDVI时序曲线不一致，导致提取面积小于实际播种面积，这种现象多反映在北京、廊坊、唐山等小麦种植不集中的区域。

本研究以低分辨率MODIS影像为数据源，在大尺度范围内对农作物进行识别，整体提取精度达到91.84%，精度结果表明基于MODIS-NDVI数据产品以决策树分类方法进行分类获得的冬小麦面积是可靠便捷的：数据源易获取，分类方法清晰易懂、执行效率高，分类结果精度高，为快速提取农作物面积提供了一定的依据。

4 讨论

本研究通过对冬小麦分布位置进行实地采样，得到冬小麦NDVI时间序列图，以真实数据为样本来提取影像中冬小麦像元，提取结果具有可靠性。结合作物物候期识别能够提取冬小麦的典型物候期，通过决策树分类法在典型物候期上设定阈值进行分类，达到快速提取冬小麦像元的目的。

由于卫星影像受到云、大气等噪声影响，构建出来的NDVI时间序列曲线存在一定波动，对地物分类造成影响，经S-G滤波后的NDVI时间序列曲线更为平滑，且曲线走势与农作物物候吻合，便于对地物特征进行分析，有利于遥感影像的解译。

京津冀地区冬小麦NDVI时间序列曲线与物候期表达一致，且根据典型物候期NDVI值提取冬小麦结果的方法是可取的：在出苗期后，冬小麦先生长，NDVI曲线表现上升趋势，到寒冷的12月份停止生长，达到第1个小波峰，这期间已无其他绿色植被生长，冬小麦NDVI值表现为最高，进入休眠期后，不再生长，NDVI曲线表现为下降趋势，因此在进入休眠期前设定为识别冬小麦的物候期之一，能达到区分冬小麦与其他绿色植被的目的。在返青期后，冬小麦快速生长，直到开花期前后，NDVI值到达波峰，这个时期红外波段与红波段反射率差值最大，能区分植被与非植被区域。在成熟期后，冬小麦逐渐停止生长，麦穗逐渐成熟，麦叶逐渐发黄，在近红外波段与红波段的反射率差值迅速降低，直至收割后出现裸地，NDVI值表现为裸地特征，因此选取成熟期后期阶段作为区分冬小麦与其他绿植的典型物候期。

在实地采样中，多发现在冬小麦种植区内有小范围的林地种植，从冬小麦播种到成熟期，出现林地最大NDVI时间序列曲线和冬小麦的最小NDVI时间序列曲线走势一致，而实际调研时该像元是林地和冬小麦的混合像元，可能是由于这段时期树木开始萌芽、抽枝展叶并生长，其NDVI值表现与冬小麦极为相似，造成曲线走势相似，然而在第161天后，冬小麦由于处于收割季而出现明显波谷，而混合像元的NDVI时间序列曲线没有明显波谷，整体上仍然处于上升趋势，这一现象为混合像元的分解提供了参考。

本研究下载的影像为2016—2017年一个整周年的影像，总体精度达到91.84%，但在地市级水平上，精度高低表现不同，因此为了提高提取精度，需要在冬小麦播种面积少且不集中的区域增加采样数量;同时，为验证此分类方法在京津冀地區实施的稳定性，下一步应考虑获取更长时间序列的NDVI数据，并对冬小麦的种植区变化情况做全面解析。

参考文献：

[1]江 东，丁方宇，郝蒙蒙，等. 农田遥感识别方法与应用[J]. 甘肃科学学报，2017，29（2）：43-47.

[2]单 捷，孙 玲，于 堃，等. 基于不同时相高分一号卫星影像的水稻种植面积监测研究[J]. 江苏农业科学，2017，45（22）：229-232.

[3]黄健熙，侯矞焯，武洪峰，等. 基于时间序列MODIS的农作物类型空间制图方法[J]. 农业机械学报，2017，48（10）：142-147.

[4]Wardlow B D，Egbert S L. Large-area crop mapping using time-series MODIS 250 m NDVI data：an assessment for the U.S. Central Great Plains[J]. Remote Sensing of Environment，2008，112（3）：1096-1116.

[5]黄 青，吴文斌，邓 辉，等. 2009年江苏省冬小麦和水稻种植面积信息遥感提取及长势监测[J]. 江苏农业科学，2010（6）：508-511.

[6]许青云，杨贵军，龙慧灵，等. 基于MODIS NDVI多年时序数据的农作物种植识别[J]. 农业工程学报，2016，30（11）：134-142.

[7]郭昱杉，刘庆生，刘高焕，等. 基于MODIS时序NDVI主要农作物种植信息提取研究[J]. 自然资源学报，2017，32（10）：1808-1818.

[8]Pan Y Z，Li L，Zhang S，et al. Winter wheat area estimation from MODIS-EVI time series data using the Crop Proportion Phenology Index[J]. Remote Sensing of Environment，2012，119：232-242.

[9]平跃鹏，臧淑英. 基于MODIS时间序列及物候特征的农作物分类[J]. 自然资源学报，2016，31（3）：503-511.

[10]李祚泳. 用B-P神经网络实现多波段遥感图像的监督分类[J]. 红外与毫米波学报，1998，17（2）：153-156.

[11]潘 琛，杜培军，罗 艳，等. 一种基于植被指数的遥感影像决策树分类方法[J]. 计算机应用，2009，29（3）：777-780.

[12]刘焕军，于胜男，张新乐，等. 一年一季农作物遥感分类的时效性分析[J]. 中国农业科学，2017，50（5）：830-839.

[13]黄健熙，贾世灵，马鸿元，等. 基于WOFOST模型的中国主产区冬小麦生长过程动态模拟[J]. 农业工程学报，2017，33（10）：222-228.

[14]申 健，常庆瑞，李粉玲，等. 基于时序NDVI的关中地区冬小麦种植信息遥感提取[J]. 农业机械学报，2017，48（3）：215-220.

[15]Wang J，Wang E，Feng L，et al. Phenological trends of winter wheat in response to varietal and temperature changes in the North China Plain[J]. Field Crops Research，2013，144（6）：135-144.

[16]Chen J，Jnsson P，Tamura M，et al. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter[J]. Remote Sensing of Environment，2004，91（3/4）：332-344.

[17]陈晓苗. 基于MODIS-NDVI的河北省主要农作物空间分布研究[D]. 石家庄：河北师范大学，2010.

[18]Pan Y Z，Li L，Zhang J S，et，al. Crop area estimation based on MODIS-EVI time series according to distinct characteristics of key phenology phases：a case study of winter wheat area estimation in small-scale area[J]. Journal of Remote Sensing，2011，15（3）：578-594.

[19]Xiao D，Tao F，Liu Y，et al. Observed changes in winter wheat phenology in the North China Plain for 1981—2009[J]. International Journal of Biometeorology，2013，57（2）：275-285.

[20]申文明，王文杰，罗海江，等. 多时相Landsat8 OLI影像的作物种植结构提取[J]. 遥感技术与应用，2015，30（4）：775-783.

[21]刘吉凯，钟仕全，梁文海. 基于策树分类法及其在遥感图像处理中的应用[J]. 测绘科学，2008，33（1）：208-211.

[22]潘 琛，杜培军，张海荣. 决基于决策树分类技术的遥感影像分类方法研究[J]. 遥感技术与应用，2007，22（3）：333-337.