黑河中游是甘肃省中部河西走廊平原的主要部分^［1］，位于祁连山区的山前平原，地形平坦、阳光充足、引水方便^［2］。该区域气候干旱、降水量少、蒸散量大^［3］。因此，黑河中游绿洲以人工绿洲为主，农作物种植结构具有高度异质性和多样性^［4］。2002年开始玉米成为黑河中游甘州区、临泽县、高台县种植面积最大的农作物，生长期需水量最大、蒸散量最高^［5］，准确获取玉米种植面积信息对该区域作物种植构调整、水资源合理规划具有重要参考意义。

随着卫星遥感的快速发展，遥感技术以其独特的大面积观测能力在植被面积监测方面取得了巨大成就，利用可监测农业信息的遥感卫星（例如Landsat卫星、SPOT卫星、Terra和Aqua卫星、高分1号卫星以及Sentinel-2卫星等）数据获取大范围农作物种植面积信息已成为一种重要手段^［6-7］。过去几年农作物分类提取研究中，卫星影像数据源以Landsat数据和MODIS数据为主^［8-9］，但是受限于数据空间分辨率较低，即便Landsat数据具有较高的空间分辨率30 m，也存在混合像元的问题，仍难以满足复杂异质性条件下的农作物种植面积制图精细制图的需求。Sentinel-2是高分辨率多光谱成像卫星，包括在同一轨道的两颗极轨卫星S2A和S2B，两颗卫星互补，具有宽幅盖（290 km）、高重访周期（5天）、高空间分辨率（10 m）优势，可支持监测植被生长季节内变化^［10］，大大减小混淆像元问题。因此本研究探讨高度异质性植被景观条件下，高分辨率的Sentinel-2影像在农作物种植面积提取研究中的应用潜力。

近几年，利用机器学习方法提取农作物空间分布的研究成为了热点。其中，基于逐像元的监督分类方法应用最为广泛，如决策树算法、最大似然分类、支持向量机以及随机森林（Random Forest，RF）等^［11-12］，同时K-means聚类、迭代自组织非监督分类算法在无标签样本的大尺度区域分类中也具有一定优势^［13-14］。已有研究表明，与最大似然分类器、支持向量机等作物类型提取方法相比，随机森林算法的精度更高^［15］。故本文选择利用随机森林算法对具有高度异质性的黑河中游绿洲开展玉米种植面积提取研究。

本文采用了两种不同的玉米种植面积提取方案，研究适合于复杂异质性植被景观条件下农作物种植面积提取的更优方法，并分析影响农作物种植面积提取精度的主要因素。此外，还比较了基于单幅影像与多时相影像农作物种植面积提取结果差异，探讨多时间信息量对分类精度的影响。

选取黑河中游张掖地区的甘州区、临泽县、高台县作为研究区（图1）。研究区是甘肃中部河西走廊平原的主要部分，地势平坦，为黑河绿洲农业集中区，气候干旱、年日照时数高、降水量少、蒸发量大的水热条件决定了该区域农业多为单季种植，农作物生长期主要集中在在4~9月间。玉米是该区域种植面积最大、生长期需水量最大的农作物。本文选取该研究区开展玉米种植面积研究，主要基于以下3点考虑：①该绿洲农作物分布具有高度异质性和多样性，可为异质性区域作物提取提供参考；②该区域玉米种植收益为农民主要收入来源之一，可为种植结构提供参考；③该区域玉米种植生长期需水量最大，可为该区域水资源供应调整提供参考。

选取2019年4月1日~9月30日研究区一个完整生长季上的428幅Sentinel-2 MSI（Multi-spectral instrument，多光谱成像仪）Level-1C产品，该产品提供大气顶层反射率影像，已经过辐射校正和几何校正。预处理过程包括重投影、重采样以及云掩膜。首先，将影像数据集重投影为WGS 1984投影坐标系；其次，Sentinel-2影像包括13个波段，其中B2/B3/B4/B8波段空间分辨率为10 m，B5/B6/B7/B8A/B11/B12波段空间分辨率为20 m，B1/B9/B10波段空间分辨率为60 m，利用三次卷积插值法对所有较低空间分辨率波段重采样，重采样后的空间分辨率为10 m。利用ESA提供的Sentinel-2影像QA60波段对428幅影像质量进行评估，掩去影像中云和卷云像元，最终获得6幅月度中值合成影像。

于2019年10月22日~11月1日对研究区农作物分布开展了野外实地考察，主要记录了玉米和非玉米植被（森林、麦子、辣椒、燕麦、苜蓿、油白菜、洋葱、西兰花、笋、土豆、甜菜等）。受采样时间限制，非玉米植被样本的采集只对植被类别显示明显的区域进行了记录。共采集了780个样本点，其中玉米样本点623个，非玉米植被样本点157个点。考虑到非玉米植被采样点数量较少，且黑河中游大棚种植基本上均为蔬菜种植，因此在Google Earth卫星影像手动采集了241个大棚种植点，且假设这些点为非玉米植被样本。另外在Google Earth还采集了林地样本38个，草地/林地样本333个，建设用地样本529个以及其他土地覆盖类型（例如裸地、稀疏植被、水体等）样本667个。除此之外，为避免在Google Earth上采集水体样本点造成地理位置偏差，在Sentinel-2影像采集了336个水体样点。因此研究区内总共采样点为2 896个（详见表1），其中玉米623个，非玉米植被731个，非植被1 542个。

随机森林是2001年 Breiman提出的一种集成学习方法，被广泛应用于统计分类和非参数回归问题。该算法主要过程是：①利用bootstrap重抽样方法对样本数据集随机采样，每次生成包含m个样本的训练样本集，采样共进行n次；②对每个训练样本集随机选择特征集或特征集的线性组合进行决策树建模，生成n个分类器，即n棵树；③通过投票从n棵决策树中选择最优的分类结果。随机森林算法能评估特征变量的重要性，具有快速、易于参数化和较强鲁棒性的特点，因此在具有许多不同特性的大规模数据集的分类运用中精度较高且能处理具有高维和噪声的输入数据，目前该方法已在土地覆被分类和农作物制图应用研究中已表现出卓越的性能^［16］。

研究构建两种不同基于随机森林算法的玉米种植面积提取方案，具体技术路线如图2所示。方案一，直接提取法，即基于多时相影像和3类样本（玉米、非玉米植被、非植被）构建RF分类模型（简称为RF-0），直接提取研究区玉米种植面积。方案二，两步提取法，首先基于多时相S2影像和训练样本，构建提取植被区的RF分类模型（简称为RF-1），得到植被区月度时间序列影像；然后，基于植被区月度时间序列影像和玉米/非玉米植被样本，构建提取玉米种植面积的RF分类模型（简称为RF-2）。方案一是目前农作物分类提取中常用的方法^［17］，方案二将植被区独立后再对玉米种植面积进行提取。本文旨在通过比较两种方案得出适合于黑河中游玉米种植面积提取的方法。

为了研究Sentinel-2影像中各个波段对植被信息提取的贡献性，选用了全部可用的13个光谱波段，包括可见光、近红外、短波红外波段。另外，还考虑了有利于植被区提取或与植被生长状况相关的17个植被指数。其中，NDVI（Normalized Difference Vegetation Index，归一化植被指数）用于定量光合能力、水分胁迫和植被生产力^［18］，同时，时序NDVI能够观测不同土地覆被类型的物候信息^［19］和植被特征^［20］；NDWI（Normalized Difference Water Index，归一化水指数）是对地表水分含量敏感的植被指数^［21］；EVI（Enhanced Vegetation Index，增强植被指数）可降低土壤反射率和大气污染，并增强对高覆盖植被区的敏感性^［22］；GCVI（Green Chlorophyll Vegetation Index，叶绿素指数）的空间分布对研究区肥料（氮）等施用的响应具有重要意义^［23］；STI（Soil Tillage Index，土壤耕作指数）和NDTI（Normalized Difference Tillage Index，归一化耕作指数）是识别耕作土地的有效指数^［24］，其他指数也经常被用于植被信息的提取。

本文共选取30个波段/指数，得到的月度合成影像共有30×6个波段/指数，因此，模型输入共有180个特征。

在随机森林算法中，将一列确定性伪随机数添加到样本数据集合中，这些数字是双精度浮点数，范围是［ 0.0， 1.0），然后按照伪随机数小于0.7和大于等于0.7划分训练样本和验证样本，即RF-0、RF-1、RF-2模型的训练样本和验证样本个数分别为2016/869、2069/816、900/394。对给定的训练数据集，使用字段随机分割方法对模型参数进行评估，在RF-0/RF-1/RF-2分类器中设置模型参数树的个数为200，随机特征变量数为总特征数的平方根以及最小叶节点数为10。

构建混淆矩阵，定义总体精度（OA，Overall Accuracy）、F1_Score以及Kappa系数3个评价指数：

其中：N表示训练样本的个数，n表示分类类别数，i表示某一具体类别。OA指在验证样本集上分类器预测正确的概率；F1_Score由用户精度和生产者精度计算得到，是用户精度和生产者精度调和均值，用于表征模型输出的好坏，取值范围位于0~1之间，其值越接近1表示模型输出越好；Kappa系数是检测验证样本与预测结果吻合度的指标，其值的大小通常位于0 ~ 1之间，越接近1表示吻合度越高^［25］。

首先选择基于Sentinel-2的所有180个特征作为输入参数，分别训练两种研究方案下的随机森林分类模型，方案一中RF-0模型总体精度为93.06%、Kappa系数为0.88；方案二中RF-1模型与RF-2模型的总体精度分别为97.68%和89.78%，Kappa系数分别为0.95和0.79，表明通过两种方案构建的随机森林模型精度均较高。

进一步利用验证样本对所构建的模型进行精度验证，为了使两分类方案间具有可比性，对方案一中的混淆矩阵进行了处理，结果如表2所示。两种方案提取的植被区OA、Kappa和F1_Score基本相同，OA为97%，Kappa为0.935左右，以及植被类型F1_Score为0.97，表明两种研究方案均能较准确地提取植被区。在玉米种植面积提取中，方案一中OA为84.48%，Kappa为0.68，玉米F1_Score为0.82；方案二中OA为85.03%，Kappa为0.70，玉米F1_Score为0.83。与方案一相比，方案二中RF-1首先将植被区独立出来，RF-2接着在植被区中分类出玉米作物，玉米错分为非植被类型的样本数明显减少，OA、Kappa系数、玉米F1_Score均明显高于方案一，此方案能有效降低分类过程中的总体类究区玉米种植面积提取研究。基于方案二的研究区玉米种植面积分布结果如图3所示。

研究区玉米种植面积提取过程受复杂因素影响，例如研究区作物灌溉方式大多为渠灌，受开渠时间等影响，地理位置相邻的大面积农作物灌溉时间相似，这部分农作物物候信息、生长健康状况相似，可能监测到的植被指数差异小，导致错分；玉米种植地边缘的杂草、套种玉米以及玉米地生长的其他果树、枣树等，这些因素均可能造成玉米与其他植被间光谱差异变小，导致分类困难；此外，山前平原土壤侵蚀严重，土质较差，而距离黑河较近的耕地土质较好，土质之间的差异导致玉米类内光谱差异变大从而导致玉米错分。

进一步采用增量式学习方法^［26］分析多期月中值合成影像数据对分类准确性的影响。该方法在数据子集上训练初始模型，随后不断加入新的数据更新该模型。首先以4月中值合成影像作为数据子集训练初始模型，并对模型进行验证，计算得出玉米的Kappa；然后在此基础上加入5月份中值合成影像再次训练模型，并再次计算此模型分类中玉米的Kappa。以此类推，不断加入新的数据直到新数据为9月份中值合成影像。基于单幅影像和多时相影像分类得到的总体分类精度如图4所示。基于单幅影像的分类中，4、5月和8月的玉米提取Kappa较高。在4月和5月该地区农田上基本压了薄膜，这些压薄膜的区域基本上种植玉米，因此玉米提取精度较高，但是由于薄膜和大棚信号很相似，导致大棚也被识别为玉米；在8月份，玉米处于抽穗期，生长旺盛，光谱特征与其它地物差异明显，因此提取精度较高。而6、7、9月精度较低，可能是由于6、7月份植被覆盖度整体较低，不同作物间光谱差异较小，存在同物异谱/异物同谱现象；9月部分成熟农作物被收割，导致不同作物间植被指数或波段值差异减小，玉米种植面积提取精度稍有下降。

基于多时相影像的分类，随模型输入影像数量的增加，玉米的Kappa明显提高。在不同月份中，基于多时相影像的提取精度均高于单幅影像的提取精度，主要是因为多时相影像能在很大程度上能克服基于单幅影像玉米提取过程中的外界因素影响，并提供更丰富的光谱信息，有效减少作物分类过程的错分、漏分现象，从而提高作物分类精度。在本研究中，到8月份为止，Kappa值达到最高（0.70），比基于单幅影像的最佳分类精度高2%。

随机森林算法中共使用了180个变量，这些变量在分类过程中的重要性各不相同。对RF-1和RF-2模型中180个变量进行特征重要性降序排序，并对排名前60变量统计其对应指数的重要性以及出现的频次（图5）。RF-1模型中排名前4的指数是NDVI、STI、NDTI、RDGCVI1（Red edge Green Chlorophyll Vegetation Index 1，红边叶绿素指数1），这4个指数与B4、B5、B8、B11、B12波段相关，这些波段均可有效监测植被生长信息。其中NDVI指数排在第一位且重要性远高于其他指数或波段，说明研究区植被提取在很大程度上依赖于NDVI。与RF-1模型特征参数重要性排序相比，RF-2模型使用指数量明显增多，这是由黑河中游复杂异质性植被结构决定的，其中排名前4的指数为STI、B10、NDT I、RDNDVI1（Red Edge Normalized Difference Vegetation Index 1，红边归一化植被指数1）。STI、NDTI指数为RF-1和RF-2共有指数，两指数均由B11、B12波段计算所得，表明短波红外波段在植被信息提取中的重要性。另外，值得注意的是，B10波段属于“大气”短波红外波段，具有水汽强吸收特性，该波段对本研究中玉米识别具有重要作用，这可能与玉米是研究区蒸散量最大的农作物相关。

本文基于Sentinel-2多时相影像，通过搜集大量有利于研究的参数指数，采用随机森林算法开展了黑河中游玉米种植面积提取研究，得出以下结论：

（1）Sentinel-2高分辨率多时相影像可用于复杂异质性景观条件下农作物种植面积精细提取，有效减少混淆像元的影响。本研究中玉米种植面积提取总体精度可达85.03%，F1_Score为0.70，Kappa系数为0.83，可满足精细农业种植面积提取的需求，为该区域农作物种植结构调整、水资源合理规划提供数据支撑。

（2）与单幅影像相比，多时相影像能监测不同作物物候信息，获取更丰富的光谱信息，可有效减少作物分类过程中出现错分、漏分现象，提高作物分类精度。

（3）玉米种植面积提取过程中，将非植被区独立出来可减小植被总体类内方差，有利于作物精准提取。作物灌溉的时空差异也可能导致作物类内光谱差异变大，而类间光谱差异减小，这些光谱上的差异使得玉米真实种植面积减小，同时伪种植面积增大。