基于MODIS-EVI的内蒙古沿黄平原区作物种植结构分析

2021-04-28贾博中白燕英魏占民张志胤

灌溉排水学报 2021年4期

贾博中，白燕英*，魏占民，闫东，张志胤

（1.内蒙古农业大学水利与土木建筑工程学院，呼和浩特010018；2.内蒙古自治区土壤肥料和节水农业工作站，呼和浩特010010）

0 引言

【研究意义】我国是人口大国，农业的发展是国民经济和社会稳定发展的重要基础和保障[1]。快速、高效地获取农作物种植结构信息对农业部门的生产管理、农产品产量估计以及农作物种植结构优化和调整都具有重要的意义[2-3]。遥感技术相比于费时费力的传统人工调查，有着成本低、观测周期短以及覆盖面积广等优势，广泛应用于作物种植结构提取，并且正逐渐成为一种热门手段[4-6]。因此，怎样应用遥感影像和采用适合的分类方法，来取得较高精度的分类结果，成为低成本较高精度农业遥感研究的一个重要方向。【研究进展】目前基于多时相作物分类方法充分地利用了作物季节特征，能清楚地反映出不同农作物随着时间的变化趋势，能有效地减少农作物误分现象，是当前用遥感进行农作物种植结构提取的主流方法[7-11]。杨小唤等[12]基于时间序列的MODIS 影像，通过分析各类作物的归一化植被指数NDVI值在生育期的生长变化规律，提取出北京市冬小麦和夏玉米等农作物。MAUS 等[13]在DTW 算法的基础上引入了时间权重因子，利用MODIS 时间序列影像，基于时间加权的动态弯曲方法在巴西中部地区提取作物的空间分布情况。陈颖姝等[14]对Landsat8-OLI 和MODIS影像的时间序列综合利用，提取出感兴趣区域，再利用BP 神经网络的方法对提取数据进行监督分类，进而获取研究区农作物的空间分布信息。韩宇平等[15]利用Landsat8 和MODIS 影像的NDVI时间序列数据，采用监督分类和非监督分类2 种方法对灌区不同水源的灌溉面积进行分类。由此可见，虽然目前学者对利用卫星遥感数据提取空间分布信息开展了广泛研究，但如何利用大尺度中低分辨率在大区域内进行种植结构分析仍值得研究。

【切入点】目前，MOD13Q1 影像能够有效地构建NDVI和EVI时间序列。由于在农作物生育高峰期或高作物覆盖下，NDVI的增长速度没有EVI快，导致生长或衰落时农作物的EVI值相比NDVI值的变化更为明显。甚至有些情况下当作物开始衰落时，EVI值显示为降低，而NDVI值仍表现为升高，不能敏锐地感应出农作物长势变化，在监测农作物长势变化方面EVI更敏感[16]。【拟解决的关键问题】本研究采用2019年的MCD12Q1 影像和2019年4―9月12 个时相的MOD13Q1 影像构建EVI时间序列，分析内蒙古沿黄平原区小麦、玉米、葵花、西葫芦、番茄和苜蓿6 种作物的生育期和增强型植被指数EVI的时间序列变化特征，提取出主要农作物的时间变化曲线，通过6 种主要农作物的EVI时序曲线，获取作物空间分布和种植面积等信息。

1 材料与方法

1.1 研究区概况

研究区位于内蒙古自治区中西部的内蒙古沿黄平原区地处黄河中上游，研究区疆域辽阔地势平坦（图1）。平原区属于典型的温带大陆性气候，降水量少而蒸发能力大，全年日照时间长。地理坐标为北纬39°01′30″―41°20′23″，东经105°22′09″―112°00′14″。内蒙古沿黄平原区拥有着得天独厚的灌溉条件，适宜多种作物生长，是我国重要的粮食生产基地，主要的粮食作物有小麦和玉米，主要的经济作物有各类瓜果蔬菜和葵花，种植结构较为复杂，在自治区经济发展中有着举足轻重的作用。

图1 研究区地理位置Fig.1 Geographical location of the study area

1.2 研究区数据获取与分析

1.2.1 MODIS 遥感数据获取

本研究所用的MCD12Q1 和MOD13Q1 影像时间分辨率分别为1 a 和16 d，空间分辨率分别为500 m和250 m，均下载自ladsweb.modaps.eosdis.nasa.gov。中分辨率成像光谱仪（MODIS）是AQUA 卫星和TERRA 卫星上装载的主要传感器之一，在2 颗卫星的相互配合下可以得出36 个波段的观测信息，并可重复观测整个地球表面。

1.2.2 MODIS 遥感数据预处理

MODIS 产品都为正弦投影坐标系，通常使用的地理坐标系为WGS-84，而在正弦投影下中国地区形变较大，需要坐标转换后使用，所以对MODIS 影像数据预处理的流程是：重投影-重采样-镶嵌-裁剪。

数据预处理使用The Environment for Visualizing Images（ENVI）软件。首先对MCD12Q1 和MOD13Q1影像进行重投影，将图像转换为WGS-84 地理坐标系；其次进行重采样，将MCD12Q1 和MOD13Q1 影像的像素分别设为500 m 和250 m；然后将覆盖研究区需要的h26v04 和h26v05 二部分图像镶嵌；最后根据内蒙古沿黄平原区的矢量数据文件对测区进行裁剪。

1.2.3 野外采样数据

为了更好地验证农作物分类结果，2019年4―9月对内蒙古沿黄平原区内的农作物种植情况进行了多次实地采样，共收集到533 个样本点，其中小麦64 个，玉米129 个，葵花72 个，西葫芦79 个，番茄75 个，苜蓿62 个和其他作物52 个。MODIS 影像空间分辨率较低，在野外采集数据时需要用Landsat8-OLI 影像中分辨率为15 m 的全色波段和30 m 的多光谱波段，配合Google Earth 地图，通过目视解译来辅助野外布点，保证每一个采样点都设置在面积相对较大的农田上，且平均分散在研究区内。在野外实地采样时用GPS 仪打点并拍照记录，尽量选取农作物种类一致、覆盖度相对较高、长势均匀且面积大于250 m×250 m 的区域作为采样点。

表1 主要农作物生育期Table 1 Main crop growth period

1.3 农作物生育期

内蒙古沿黄平原区农作物的生育期为4月初到9月底，本文主要对研究区2019年种植面积相对较大的小麦、玉米、葵花、西葫芦、番茄和苜蓿的种植情况进行遥感获取，剩余农作物归为其他作物一类，6 种主要农作物生育期详见表1。

1.4 混淆矩阵的计算方法

混淆矩阵可以进一步得到农作物总体分类精度（PC）和Kappa 系数（k），计算式为：

式中：n为混淆矩阵的总行列数；Pkk为混淆矩阵的对角线元素；P为总像元个数。

式中：N为总像元个数；r为混淆矩阵的总行列数；Xii为混淆矩阵的对角线元素；Xi+为混淆矩阵各行之和；X+i为混淆矩阵各列之和。

2 结果与分析

2.1 土地利用分类

MCD12Q1 影像基本的土地覆盖分类确定了17 种土地覆盖类别，其中包括11 种自然植被类别，3 种土地类别以及3 种非植被土地类别，这些是通过MCD12Q1影像本身自带的决策树分类方法得出（图2）。本研究应用的是全球植被分类（IGBP），提取其中DN 值为12 的农用地在ENVI软件下进行掩膜处理，需要在内蒙古沿黄平原区内掩膜掉非农用地部分，掩膜后剩余的部分作为下一步作物提取分类的研究区。

图2 MCD12Q1 土地利用分类Fig.2 MCD12Q1 Land use classification

完成分类后统计各类所占比例，得出内蒙古沿黄平原区总面积约为3.86×104km2；提取出农用地面积约为1.03×104km2，占研究区总面积的26.55%，把这部分得到的分类结果作为农用地研究区。

2.2 不同农作物EVI 值时间特征分析

根据小麦、玉米、葵花、西葫芦、番茄和苜蓿6种作物在野外采样时得到的经纬度，在ENVI软件上可以显示出每个采样点EVI的时序曲线走向，并得到每个采样点时间序列EVI值，合并起来即为EVI时间序列曲线（如图3）。同一种农作物在不同的采样点由于长势、种植方式等因素的不同导致EVI时间序列曲线也不完全相同，但其趋势基本相近。在分析研究区目标农作物时间序列EVI值特征的基础上，可以提取6 种主要农作物EVI的最大值、最小值、平均值以及峰值出现的时间这4 个特征阈值作为典型特征，用于辅助提取6 种主要农作物EVI时间序列曲线。

图3 农作物EVI 时间序列曲线Fig.3 Crop EVI time series curve

小麦的EVI值在4月下旬出苗时开始增加，在6月下旬抽穗时达到最大值，并在7月下旬成熟收获时会迅速下降。内蒙古沿黄平原区部分区域春小麦收割时间早，如果收割后温度适宜，且灌溉条件好，部分地区还会在收割后的春小麦区域重新种植其他农作物，导致有的地区9月小麦EVI时序曲线有一定幅度上升，达到了1年2 个波峰。因为峰值和峰值出现时间是区分小麦和其他作物的关键，所以在提取小麦时找到这2 个因素尤为重要。

玉米的EVI值在4—5月都处于较低范围，EVI值小于0.25，此时的玉米处于生长初期，到了6—8月，玉米的EVI指数一直呈较快速度增长，此时玉米处于快速生长期，EVI值于8月初达到峰值。

葵花4—6月EVI值始终较小，EVI值低于0.2，其在此时期处于生长初级阶段，7月初至8月底，葵花的EVI值迅速增加，此时葵花处于快速生长期，8月底葵花EVI指数达到峰值，并在9月呈下降趋势，在9月下旬进入成熟期。值得注意的是葵花是内蒙古沿黄平原区最晚收割的农作物，其峰值出现时间在农作物提取中尤为关键。

西葫芦的EVI指数在4—5月都处于较低范围，EVI值均小于0.2，此时西葫芦处于生长初期，5—6月EVI指数增速加快，到了6—7月西葫芦EVI指数增速减缓，并在7月初达到峰值。

番茄与西葫芦生长情况较为相似，在4—5月EVI指数处于较低范围，EVI值均小于0.2，5—6月EVI指数快速增加，6—7月番茄EVI指数增长速率减慢，并在7月中旬达到峰值。

苜蓿和其他农作物的EVI时间序列曲线差异明显，在研究区的地理和气候条件下，当地为了实现饲草料就近就地化解决，从而降低奶牛养殖成本增加收益，苜蓿每年收割3 次，间隔1—2 个月就会出现1次EVI值先上升再下降的情况，相应在EVI时序曲线升降较为频繁，便于区分。

2.3 作物种植结构提取

本研究根据农作物EVI时间序列曲线和农作物样本点分别列出6 种主要农作物的分类提取条件，其中EVI时间序列曲线的波形可以反映出农作物的生长变化情况。依据不同农作物生长的物候特征，利用作物分类提取条件（如图4）在MOD13Q1 影像ENVI软件中对研究区内的农作物进行分类提取。

图4 作物分类提取条件Fig.4 Crop classification extraction conditions

表2 各类农作物分类像元数及面积占比Table 2 Image number and area ratio of various crop classification

图5 主要农作物分类Fig.5 Main crop classification

按上述分类提取条件（如图4）和分类结果（如图5）可得，农用地大部分被玉米所覆盖占比达到72.29%，面积为7 412.75 km2，尤其在研究区东部的呼和浩特和包头周围地区种植意向更倾向于玉米。葵花面积为575.26 km2，占农用地总面积的5.61%，主要分布在乌拉特前旗和五原县。小麦在研究区内仅占农用地的0.97%，面积为99.47 km2，集中在杭锦后旗一带。西葫芦和番茄的面积分别为900.06 km2和836.75 km2，分别占农用地的8.78%和8.16%，在研究区分布较为广泛。苜蓿主要集中在土默特左旗，面积为183.55 km2，占农用地的1.79%。其他作物面积为246.10 km2，占农用地的2.40%，得到的统计结果与野外实地采样情况基本符合。

2.4 精度验证

根据野外采样点结合6 种主要农作物EVI时间序列曲线，在分类结果图像中随机提取每种作物50 个样本点，共计350 个点进行农作物分类结果精度验证，结果如图6 所示。

图6 农作物分类混淆矩阵热度Fig.6 Crop classification confusion matrix heat

农作物的总体分类精度为78.29%，Kappa 系数为0.747，精度较为理想。从用用精度和制图精度分析，除去其他作物的用用精度和制图精度较低，分别为68.75%和66%以外，小麦、玉米、葵花、西葫芦、番茄和苜蓿的用用精度和制图精度都较好。

因为研究区玉米种植面积超过70%，所以获取了沿黄灌区各旗县2019年玉米种植统计年鉴面积统计结果作为分类效果的依据。由统计结果可知，各旗县2019年玉米种植面积为：阿拉善10.47 km2、鄂尔多斯2 287.13 km2、巴彦淖尔2 617.53 km2、包头1 047.50 km2、呼和浩特1 949.54 km2，在沿黄灌区种植的玉米实地统计面积合计为7 912.17 km2。而从表2 得知玉米的提取面积为7 412.75 km2，相对误差为6.31%，得到了较为满意的结果。

3 讨论

作物种植结构是作物长势监测和估产分析、种植结构调整及优化、作物灌溉管理的主要依据，充分利用作物的季相节律特征是区分作物与其他绿色植被以及作物之间的关键理论依据[17]。本研究尝试在内蒙古沿黄平原区的粮食生产基地建立快速、准确并且相对有效的作物分类识别方法，使用MODIS-EVI时间序列曲线探索研究区6 种主要农作物的分布情况。研究结果表明，解译精度虽然未能达到张威等[1]、邓荣鑫等[18]、黄思宇等[19]、魏鹏飞等[20]研究结果的精度，学者使用高分辨率影像分类文献中的平均精度，但在局部区域和整体研究区上也取得了较好的分类结果。小麦、玉米、葵花、西葫芦、番茄、苜蓿和其他作物的用用精度分别为：79.59%、80%、83.67%、78.18%、75.93%、82.22%、68.75%，制图精度分别为78%、80%、82%、86%、82%、74%、66%。农作物总体分类精度达到78.29%，Kappa 系数为0.747，而且在沿黄灌区种植面积最大的玉米与实际统计相对误差为6.31%。证明了基于大区域大尺度中低分辨率的EVI时间序列分类方法在识别小麦、玉米、葵花、西葫芦、番茄和苜蓿等农作物有较强的适用性。

由于研究区内农作物种类较多、面积较广，在野外实地采样中未能将研究区内全部作物进行采样，未得到研究区内所有作物的EVI时间序列曲线，也未能够对全部作物进行分类，只提取了数量相对较多的玉米、葵花、小麦、西葫芦、番茄和苜蓿6 种农作物。在农作物分类混淆矩阵中其他作物的用用精度和制图精度仅为68.75%和66%，该方法如需进一步提高分类精度，就要对研究区内更多作物类别如：糖菜、马铃薯和大豆等进行分类提取，从而提高其他作物的用用精度和制图精度，来进一步提高农作物总体分类精度和Kappa 系数。此外，在实地采样过程中一个像元250 m×250 m 的范围内存在插花种植现象影响EVI时间序列曲线的准确性，从而影响分类精度。研究区部分农作物生育期差异较小，仅凭借主要作物生育期得到的农作物EVI时序曲线来确定作物分类提取条件再进行作物种植结构分析，不能够对这部分生育期较为相近的作物进行很好的分类提取。针对本文研究的不足之处下一步应着重在像元内插花种植现象和生育期差异较小的农作物分类提取这两部分进行深入分析。在今后的研究中，应融合中高空间分辨率、低时间分辨率的遥感影像（如高分系列卫星影像、Landsat 系列卫星影像）和低空间分辨率、高时间分辨率的遥感影像（如MODIS 卫星影像），开展多源多尺度的种植结构分析研究，进一步提高分类精度[22]。