基于MODIS数据的漳滏河灌区作物种植结构提取

2022-12-13高叶鹏刘海新王冬利文韶鑫钱以临

江苏农业科学 2022年22期

高叶鹏，刘海新，王冬利，文韶鑫，钱以临

(河北工程大学矿业与测绘工程学院，河北邯郸 056038)

作物种植结构作为农业生产活动结果的直接体现，可以反映农作物的种植类型、种植布局、种植面积以及种植模式的综合农业信息[1]。快速、准确地获取农作物种植结构的时空动态变化信息，是研究作物长势、优化种植结构的重要依据，也是指导农业生产、估测农作物产量、实施农作物精细化管理的重要支撑[2]。因此，开展作物种植结构提取研究具有非常重要的现实意义。目前，农作物种植结构提取的研究主要考虑遥感数据源的选择、特征变量的选取以及分类方法的选用3个方面[3]。在数据源的选择方面，目前常被学者用于研究的是高时间、低空间分辨率的MODIS和中高空间分辨率、时间分辨率较低的Sentinel-2、Landsat、高分系列等数据[4-6]。在特征变量选取方面，常被国内外学者用来研究农作物分类的特征变量是纹理特征、光谱特征、植被指数等，其中植被指数作为表征地表植被生长状况的重要参数应用最广泛[7-13]。在分类方法选用方面，国内外学者较常用且分类精度较高的是随机森林、决策树、面向对象等分类方法，这些方法都较成熟，已经得到广大学者的验证与应用[14-17]。但是根据已有研究结果[18-20]，这些分类方法都存在以下问题：一是需要人工采集一定数量的训练样本；二是需要过多的人为干扰；三是自动化程度较低。因此，本研究基于遥感技术在种植结构领域的应用，选用MODIS数据作为数据源，以增强植被指数(EVI)作为主要特征变量，分析基于决策树和广义DEM模型分类方法的河北省邯郸市漳滏河灌区主要作物种植结构提取，对比评价2种分类方法的提取精度和分类结果，分析讨论广义DEM在提取作物方面的适用性以及提取优势。同时，获得区域内农作物的种植结构信息，以期为漳滏河灌区作物的精细化管理和可持续发展提供科学依据。

1 研究区概况与数据处理

1.1 研究区概况

1.1.1 地理位置漳滏河灌区位于河北省邯郸市中部，介于36°15′～36°55′N、113°50′～115°20′E之间，范围涉及磁县、永年区、曲周县、鸡泽县、临漳县、成安县、魏县、广平县、肥乡区、馆陶县、大名县、邯郸市区(邯山区、丛台区)、峰峰矿区，是全国特大型灌区之一，灌区内总耕地面积达到26多万hm2。灌区内主要取水模式为水库引水，大多来源于滏阳河上的东武仕水库和漳河上的岳城水库，另灌区内建设有总干渠、干渠、分干渠等多个灌溉渠道，可实现2个水库之间的联合调度、水源共享(图1)。

1.1.2 农作物种植情况漳滏河灌区属于温带半干旱、半湿润季风气候，气候条件适宜小麦、玉米、辣椒等农作物的种植以及苹果、梨、桃等林果的生长，且该地区冬小麦和夏玉米大多采取连种的方式，其主要农作物物候期见表1。

表1 主要农作物物候信息

1.2 数据来源及处理

1.2.1 遥感数据本研究采用从美国国家航空航天局(National Aeronautics and Space Administration，NASA)官网(https：//ladsweb.modaps.eosdis.nasa.gov/search/)获取的2018年空间分辨率为250 m，16 d合成的MOD13Q1影像，共计23个时相作为时间序列数据集，并利用MRT(MODIS Reprojection Tool)对其进行拼接、重投影、波段提取，通过裁剪以获取研究区增强植被指数(EVI)、红光波段反射率、近红外波段反射率等多个数据产品。

1.2.2 样本选择本研究结合研究区2018年谷歌地球专业版(Google Earth Pro)高分辨率历史影像，根据研究区内主要的地物类型，通过目视解译的方法，选取研究区冬小麦—夏玉米、辣椒、大棚蔬菜、其他地物(建设用地、林地等未识别地类)4类样本，并保证在选择过程中样本的随机性和均匀分布。选择的样本数据主要分为2个部分用途，一部分用于作物分类方法的时间序列构建，一部分用于作物分类结果的精度验证，具体样本选择见表2。

2 研究方法

2.1 决策树模型

遥感领域的决策树分类模型是一种自上而下且根据不同地物的特殊属性，结合先验知识以及数据的统计分析基于训练集构造的1个以二叉树和多叉树形式展现的分类规则模型[21]。决策树分类模型具有分类规则灵活、分类数据多源、运算高效等特点，在作物识别方面具有广泛应用[22-24]。本研究利用增强植被指数(EVI)、比值植被指数(RVI)、红光波段反射率、近红外波段反射率等多源数据，共同参与决策树分类规则的制定，以更准确地提取不同地物。

表2 样本个数统计

2.1.1 MODIS时间序列的构建刘俊伟等认为，在提取多种作物时，尤其是作物周期存在一定重叠的，仅采用单一特征变量来提取作物存在一定困难，可能会出现较多错分和漏分现象，提取精度会降低，因而通常会选择多个特征变量作为约束条件进行作物提取[25-26]。而作物的光谱特征和光谱指数作为表征作物生长过程中重要的指标，可以反映不同作物对波段信号反射的差异性以及作物的生长状况，故常被组合共同用于作物的提取。因此，本研究选择EVI、RVI、红光波段反射率、近红外波段反射率4个特征变量作为约束条件来提取作物，并基于MODIS时序数据利用样本分别构建EVI、RVI、ρR、ρNIR的时序曲线，根据统计的均值数据绘制其变化曲线(图2)，其中RVI由ρR、ρNIR计算得出。各个时相所对应的日期见表3。

2.1.2 基于MODIS数据决策树模型构建决策树技术路线见图3。结合表1、图2可知，区域内对作物识别能够形成干扰的其他地物主要是建筑、道路等非农业用地，故先将非农业用地从识别影像中剔除掉，观察发现非农业用地的EVI值在整个时间序列中都处于较低的水平，且在第7时相非农业用地的EVI值都小于其他作物的EVI值，经试验设置EVI(7)<0.21作为约束条件之一，可以提取出大部分非农用地。为保证非农业用地充分剔除掉，引入近红外波段反射率作为另一个约束条件，由图2-d 可知，非农业用地在近红外波段反射率第1时相其反射率远低于其他地物，且在第1时相非农业用地的最大值都小于其他作物的最小值，因而根据样本统计的非农业用地第1时相反射率最大值0.19作为阈值，形成约束条件NIR(1)<0.19。因此，利用约束条件EVI(7)<0.21、NIR(1)<0.19剔除影像内的非农业用地。

表3 MODIS数据时相

由图2-a可知，大棚蔬菜在7月下旬，即第14时相其EVI值与其他作物有明显区别，这是由于其他作物在7月下旬左右正处于生长旺盛期EVI值较高，因而可以利用该时相设置阈值提取大棚蔬菜，根据样本统计值设置约束条件为EVI(14)<0.4。另外，由图2-b可知，大棚蔬菜的RVI在第14时相其值也远小于其他作物，为避免提取出来的大棚蔬菜可能混有部分其他作物，将RVI作为另一个约束条件进一步提取大棚蔬菜，经反复验证设置约束条件RVI(14)≤5较合适。因此，最终利用规则EVI(14)<0.4、RVI(14)≤5来提取大棚蔬菜。

辣椒与冬小麦、夏玉米等作物存在周期重叠问题，单依靠EVI很难完全把辣椒提取出来，而辣椒在红光波段比其他作物有较强的反射率，因而本研究利用EVI和红光波段反射率提取该作物。由图2-a 可知，辣椒在第7时相与冬小麦—夏玉米的EVI值有明显差距，且辣椒大部分EVI值都集中在0.38附近，可以通过设置合理的EVI区间阈值将绝大部分辣椒提取出来，根据样本统计值设置区间0.35≤EVI(7)≤0.45提取辣椒较合适。另外，由图2-c可知，辣椒的红光波段反射率在春季时较冬小麦、夏玉米都高，故可通过该时期将辣椒与冬小麦—夏玉米区分开，经过试验发现设置RED(2)≥0.12作为约束条件可以进一步提取辣椒。因此，选择提取辣椒的约束条件为0.35≤EVI(7)≤0.45、RED(2)≥0.12。

冬小麦和夏玉米的种植模式为轮种，故提取其中1种作物就可以获得另一种作物的种植结构。4月上旬是冬小麦生长最旺盛的时期，该时段冬小麦的EVI值明显高于其他地物，使用第7时相冬小麦EVI最小值0.42作为阈值，大于该阈值的即为冬小麦和少许辣椒。因提取物里面还含有部分辣椒，需要进一步进行提取，因辣椒与冬小麦在第7时相红光波段反射率有明显差距，其最小值大于冬小麦的最大值，故使用辣椒红光波段反射率最小值0.08将两者分离，即约束条件RED(7)≤0.08。因此，选择EVI(7)≥0.42、RED(7)≤0.08作为提取冬小麦—夏玉米的条件。

2.2 广义DEM模型

2.2.1 广义DEM基本理论数字高程模型(DEM)是地形高程数据对地形曲面的数字化表达。通过DEM可以获取区域地形因子(坡度、坡向、坡向变率、地形起伏度等)和特征地形要素(山顶点、脊点、山脊线、山谷线等)，使其常应用于土壤、水文、地质灾害、地貌等多个领域。DEM还可以扩展为广义DEM，即在二维地理空间上具有连续变化特征的各种非高程属性，按照DEM的构建方法，使其具备类似地形上的高程属性，进而可以通过地形分析反映其空间分布特征。

根据广义DEM思想[27-28]，可以将遥感植被指数作为代表植被信息的DEM，进而以植被指数为基础计算植被指数的“地形因子”和“特征地形要素”，根据这些“地形因子”和“特征地形要素”携带的空间信息，可以对农作物进行识别和提取。

2.2.2 广义DEM模型原理和实施流程选择提取作物植被指数最大和最小时期的遥感数据做差值运算，通过差值运算增强植被信息，使得不同的地物表现出不同的差值大小，而同一地物的差值将会集中在某一数值附近，根据作物的物候信息，作物植被指数差值将会集中在一个高值区域，故可选定合适的阈值实现对作物的提取。

广义DEM模型提取作物基本操作流程见图4。以冬小麦作物提取为例，其具体实施流程如下：第一，获取研究区遥感影像，差值增强植被信息。选择冬小麦植被指数最大和最小的两期影像，并用最大和最小植被指数相邻的2个时相做最大值合成和最小值合成，合成后通过差值运算增强农作物信息，获得作物差值增强影像。第二，根据广义DEM思想，提取作物差值影像上特征要素“山顶点”作为提取作物的候选样本。其中“山顶点”提取步骤包括：利用邻域统计最大值函数统计差值影像指定邻域范围内所有输入像元值的最大值；利用邻域最大值影像与原始差值影像做差得到“山顶点”影像(山顶点即差值为0的栅格)；利用“山顶点”影像和原始差值影像提取出“山顶点”对应的植被指数值。第三，“山顶点”提取的作物样本可能存在一定数量的错误样本，须要对其进行二次筛选。因而，利用广义DEM思想中“山顶点”一定位于“山脊线”上的特点，将其作为二次筛选样本条件，将错误样本去掉。其中，提取山脊线的原理是根据地面的坡向变率及地面正负地形，取正地形上的坡向变率值较大的部分即为山脊。而坡向变率值较大部分的提取则根据获取的坡向变率统计值，选用标准方差作为分类间断点将坡向变率分为若干类，取最大的一类即为坡向变率值较大的部分。故最终将“山脊线”与“山顶点”叠加进行样本筛选，取“山顶点”对应的植被指数的最小值作为作物的提取阈值。

2.2.3 基于MODIS数据的广义DEM方法作物提取

2.2.3.1 数据选取与增强由图2-a可知，夏玉米在7月下旬生长最旺盛，EVI达到峰值，在其播种期6月和成熟期9月EVI值最小，但由于9月绿色植被生长旺盛容易对夏玉米的提取产生干扰，故本研究选择7月下旬影像作为最大值、6月中旬影像作为最小值，为增强图像，使用MODIS第14、15时相进行最大值合成，第11、12时相进行最小值合成，二者差值运算获得差值增强影像。

辣椒在5月上旬到达生长旺盛期，EVI值达到最大，故可以选择MODIS中相对应的第8时相作为最大值影像；春植辣椒一般在2月进行播种，此时对应的EVI值最小，所以可以选择MODIS中相对应的第4时相作为最小值影像。根据图像增强处理，选择MODIS第7、8时相进行最大值合成，第3、4时相进行最小值合成，并进行差值运算获得增强影像。

由图2-a还可知，在10月中旬各种作物都基本已经成熟或刚播种，所以EVI值较低，而大棚蔬菜在10月中旬EVI值较其他作物都较高，故选择10月中旬作为大棚蔬菜的最大值影像；观察EVI时序曲线在10月以后，冬小麦、辣椒的EVI值开始逐渐升高，而大棚蔬菜的EVI逐渐开始降低，在12月左右其他作物EVI值比大棚蔬菜都高，故选择12月中旬作为大棚蔬菜的最小值影像。因此，选择MODIS第18、19时相进行最大值合成，第22、23时相进行最小值合成，对图像做增强处理，进而通过差值运算获得差值影像。

2.2.3.2 EVI广义DEM“山顶点”提取以冬小麦—夏玉米提取为例，根据广义DEM“山顶点”提取原理，首先应用邻域最大值统计算法，采用3×3窗口对冬小麦—夏玉米的差值影像取邻域最大值获得最大值影像，用该影像减去原始差值影像，差值结果像元值为0的即为冬小麦—夏玉米EVI差值影像的“山顶点”，进而提取对应“山顶点”的EVI差值数值，构建冬小麦—夏玉米提取样本数据集。

2.2.3.3 二次筛选样本点利用广义DEM思想，“山顶点”必定位于“山脊线”上，采用3×3的移动窗口计算“坡向”和“坡向变率”，统计“坡向变率”的标准方差，将其作为分类间断点，提取“坡向变率”值较大部分即为山脊线，将“山脊线”与“山顶点”相交，二次筛选样本的EVI差值数值，取其最小值作为冬小麦—夏玉米的提取阈值，其样本筛选见表4。由于辣椒和冬小麦—夏玉米生长期的重叠，故在提取辣椒时应将冬小麦—夏玉米提取结果先从提取辣椒影像中掩膜去掉，再进行辣椒和大棚蔬菜的提取，进而对作物辣椒和大棚蔬菜按照相同方法和流程执行“2.2.3.2”节、“2.2.3.3”节部分的步骤，得到样本筛选表见表4。

表4 主要作物样本筛选

3 结果与分析

3.1 主要作物提取结果

本研究利用决策树和广义DEM分类方法，获取漳滏河灌区主要作物种植结构，结果见图5、图6。由图5和图6可知，冬小麦—夏玉米是漳滏河灌区主要的种植作物，其分布范围广、种植面积大，且种植区域基本固定，而辣椒、大棚蔬菜种植面积相对较小且较分散，但这2种作物主要的种植地方都较集中。从决策树分类结果来看，辣椒主要位于漳滏河灌区的东北地区，即鸡泽县、邱县和曲周县，其余地方零散分布；而广义DEM辣椒分类结果同样主要集中在这3县范围内，但其大部分都集中在鸡泽县附近。贾令晨认为，辣椒是鸡泽县的特色优势产业和支柱产业，全县常年大规模种植辣椒，已经形成了集辣椒种植、加工、销售完整的产业体系，因而鸡泽县也被誉为“中国辣椒之乡”[29]。近几年，随着辣椒产业的发展也带动了周边曲周县、邱县、永年区等地区开始种植辣椒，逐渐形成了以鸡泽县为中心，其他地区为辅的辣椒种植基地，这与本研究决策树和广义DEM的分类结果基本相似。

由图6可知，2个分类结果大棚蔬菜都主要分布在永年区和曲周县，但二者在面积和分布方面存在一定差异，其中决策树提取的大棚蔬菜都主要集中在永年区中部地区，在曲周县只有少部分，而广义DEM提取的大棚蔬菜在永年区中部地区较少，部分分布于永年区与曲周县东北交界处，绝大部分分布在曲周县中部偏西地区。

3.2 精度对比

对决策树和广义DEM方法的分类结果精度验证进行统计，结果见表5、表6。由表5、表6可知，决策树的总体分类精度为83.37%，κ系数为0.78，广义DEM的总体分类精度为85.33%，κ系数为0.80，2种方法的总体精度都达到80%以上，分类效果都较理想。但与决策树分类方法相比，广义DEM方法的总体分类精度略高于决策树，较决策树精度提升了1.96百分点，且各作物的分类精度也都有小幅度提升。因而就提取精度和使用数据量的多少而言，广义DEM在数据量较小的情况下也能够实现对作物的识别和区分。

3.3 方法对比评价

决策树分类方法在作物识别领域已经得到广泛的验证和应用，该方法对于各种尺度的范围、各种时间和空间尺度的数据都具有良好的适用性。比较决策树和广义DEM分类方法二者各有其特点和优势，决策树方法对数据量需求较大，需要人为采集大量的训练样本，且实施过程需要人为干预选取阈值，决策树分类方法可以将多种数据结合到一起，共同参与决策。相较于决策树分类方法，广义DEM分类方法对数据量需求较小，人工干预较少，可以自动选取样本，自动化程度较高，但其通常只能选取1种数据作为数据源，对于某些生长周期存在重叠、植被指数数值相近的作物，提取可能较麻烦。

4 结论与讨论

本研究利用MODIS时间序列数据采用决策树和广义DEM方法提取漳滏河灌区的主要作物，并对2种方法进行对比评价，得到以下结论。第一，广义DEM作为一种新的作物提取方法，具有很好的可性行和适用性，相对于决策树而言，不仅对数据量要求较小，且人工干预较少，自动化程度较高，有较好的提取效果。第二，广义DEM方法与决策树对漳滏河灌区的作物提取结果较一致，其中大宗作物冬小麦—夏玉米基本一致，辣椒和大棚蔬菜在面积和分布上存在一定差距，但基本与作物的主要种植区域吻合。