引言

玉米是世界三大粮食作物之一，在农业生产中占据重要地位，客观、及时、准确提取玉米种植信息，对保障国家粮食安全和经济持续发展尤为重要^[1-4]。夏玉米精细化提取依赖于对夏玉米和大豆玉米带状复合种植的精确提取，特别是大豆玉米带状复合种植推广面积连年增长，使得对夏季作物的田间管理等有了新的要求^[5]，准确估计夏玉米及大豆玉米带状复合种植面积对农业生产有重要意义。

早期，通过人工地面抽样调查和层层上报等方式获取农作物种植面积信息，存在效用低、数据空间分布不连续的问题^[6]。遥感自诞生以来，以其快速、连续、高效的特点为监测农作物种植面积信息提供了极佳的数据源^[7-8]，遥感结合抽样调查的技术手段调查农作物面积在农作物种植面积提取方面具有潜在的应用价值^[9]。近年来，随着高时空分辨率遥感数据日渐丰富，遥感数据已经具备了明显的大数据特征，如何利用不断更新的遥感影像精确提取作物种植面积，继而服务于农业生产一直是农业遥感关注的主要问题^[10-11]。

2015年谷歌地球引擎(Google Earth Engine，GEE)公测上线，标志着利用云计算进行遥感分析的时代正式到来^[12]。GEE是一个专门处理卫星影像和其他地球观测数据的云端运算平台(https://earthengine.google.com/)，由谷歌、卡内基梅隆大学、美国地质调查局共同开发，能够支持海量的数据运算与地理信息数据可视化^[13]。此平台存储了大量的遥感数据，包括中分辨率成像光谱仪(Moderate-resolution Imaging Spectroradiometer，MODIS)、美国陆地卫星(Landsat)和哨兵卫星(Sentinel)等系列影像数据集及各类产品数据集，在相关领域的研究中提供了极大的便利。宫诏健等^[14]利用MODIS数据构建玉米归一化植被指数(normalized difference vegetation index，NDVI)时序曲线，采用动态振幅阈值等方法提取了玉米关键物候期，利用决策树分类法估计了辽宁春玉米种植面积，精度超过89.29%。张健康等^[15]和李莉等^[10]运用多时相Landsat和MODIS增强植被指数(enhanced vegetation index，EVI)遥感影像数据，采用监督分类与决策树分类相结合的方法对研究区主要农作物进行遥感解译，结果表明多时相数据分类精度高, 能较好地反映作物的分布状况。单捷等^[16]利用Radarsat-2影像和支持向量机及最大似然法对各时相水稻种植面积进行提取，精度均超过90%。谷祥辉等^[17]利用Sentinel-2时间序列数据和多种植被指数，采用多种分类方法进行作物分类研究，结果表明随机森林法结合时间序列植被指数的分类效果最佳。Huang等^[18]和Fan等^[19]基于Sentinel-2时间序列影像使用随机森林和分段逻辑函数算法进行作物提取，在准确度和时效性方面显现出较大优势。Bhogapurapu等^[20]利用GEE和Sentinel-1地距多视影像(Ground Range Detected，GRD)数据的完整极化信息，提取双极化指标，用于小麦和油菜的生长评估。叶智燕^[21]利用Sentinel-2计算的EVI时序曲线锁定冬小麦关键物候期，并结合Sentinel-1影像，极化方式为垂直-垂直极化(vertical-vertical polarization，简称为“VV”)的后向散射系数在冬小麦关键物候期的特征，实现冬小麦的种植面积提取。董心君等^[22]利用动态归一化后向散射系数(normalized difference sigma naught index，NDSI)实现水稻快速分类。夏季作物种植结构复杂，大豆玉米带状复合种植的不断推广使得遥感用于夏玉米(大豆玉米带状复合种植)的提取面临更多难点。目前，夏季作物遥感提取对夏玉米(大豆玉米带状复合种植)精细化提取研究仍较为缺乏，已有研究大多局限于中大尺度或单一指标，如MODIS、Landsat相关指标及Sentinel-1的VV波段、Sentinel-2的NDVI等，基于多源数据与多种指标的10 m级空间分辨率夏玉米(大豆玉米带状复合种植)精细化空间分布特征有待深入研究。

此文拟基于GEE云平台构建夏玉米遥感识别最佳时相、夏玉米遥感识别光学时序和夏玉米遥感识别光学与星载合成孔径雷达(synthetic aperture radar，SAR)融合时序及其相关指数如NDVI、EVI、地表水体指数(land surface water index，LSWI)，构建3种分类数据集应用于夏玉米种植区提取中，并分析各方法精度。首先，基于米级分辨率数据和带有地理信息的实地调查图片制作夏季作物分类样本点；其次，基于GEE云平台按照不同时间段统筹3种数据集及相应指数；最后，将样本和数据集作为GEE随机森林分类器的输入，导出分类结果并进行整体和样例的精度验证，比较分析不同数据集情境下夏玉米种植区域提取结果。

1 材料与方法 1.1 研究区概况

研究区为山东省济南市商河县(图 1)，地处山东北部，地理位置为116°58′~117°26′E、37°06′~37°32′N，地势平缓、河流较多，属大陆性暖温带半湿润季风气候，年均气温为12.3 ℃，年降水量大约为600 mm。商河是典型农业县，种植结构较为单一。由于前茬作物的不同，夏玉米播种时间存在差异，主要分为大蒜茬夏玉米和小麦茬夏玉米。

1.2 研究数据 1.2.1 地面数据采集与样本构建

利用外业数据采集地图软件外业精灵分别于2023年7月中旬和8月中旬在商河现场采集样本，现场采集的地面调查地块共201处。基于现场采集样本和卫星影像，制作间距大于30 m的样本点数据集，地类主要包括玉米、大豆玉米带状复合种植和其他类型，样本点总数为3 987个(图 1)。在样本点中随机抽取90%进行训练，其余样本点用于精度验证。

1.2.2 分类影像选择

本研究以Sentinel-1和Sentinel-2为数据源(表 1)，Sentinel-1目前仅有Sentinel-1A在轨工作，重访周期为12 d，搭载基于C波段雷达成像系统；Sentinel-2由Sentinel A、B两颗卫星组成星座，组成星座后重访周期仅为6 d，搭载多光谱成像仪有13个通道。本研究采集2023年6月1日—10月31日的Sentinel-2大气校正产品(Sentinel-2 L2A)和Sentinel-1 GRD产品。

1.3 研究方法 1.3.1 技术路线

为准确探究不同分类指标组合对夏玉米和大豆玉米带状复合种植识别的影响，基于GEE云平台和不同分类指标组合进行夏玉米和大豆玉米带状复合种植面积提取，并进行精度验证。选取Sentinel-1/2构建包含波段特征和相关指数的夏玉米遥感识别最佳时相、夏玉米遥感识别光学时序、夏玉米遥感识别光学与SAR融合时序特征数据集3种方案(表 2)，分析不同分类指标组合情境对夏玉米及大豆玉米带状复合种植提取精度的影响(图 2)。

1.3.2 数据集构建

Sentinel-1影像经过按照时间、研究区和工作模式的筛选处理，Sentinel-2影像经过按照时间、研究区的筛选以及利用Cloud Score+的去云处理。选取Sentinel-1的VV、垂直-水平极化(vertical-horizontal polarization，简称为“VH”)波段及SDWI^[23][公式(1)]，Sentinel-2的6个波段(蓝、绿、红、近红外、短波红外1、短波红外2)、NDWI^[24][公式(2)]和3个常用植被指数NDVI^[25][公式(3)]、EVI^[26] [公式(4)]和LSWI^[27][公式(5)]共10个指标参与监督分类(表 2)。

$ I_{\mathrm{SDW}}=\ln \left(10 \cdot \rho_{\mathrm{vv}} \cdot \rho_{\mathrm{vH}}\right)-8, $

(1)

$ I_{\mathrm{NDW}}=\frac{\rho_{\text {GREEN }}-\rho_{\mathrm{NIR}}}{\rho_{\text {GREEN }}+\rho_{\mathrm{NIR}}}, $

(2)

$ I_{\mathrm{NDV}}=\frac{\rho_{\mathrm{NIR}}-\rho_{\mathrm{RED}}}{\rho_{\mathrm{NIR}}+\rho_{\mathrm{RED}}}, $

(3)

$ I_{\mathrm{EV}}=\frac{2.5\left(\rho_{\mathrm{NRR}}-\rho_{\mathrm{RED}}\right)}{\rho_{\mathrm{NIR}}+6 \rho_{\text {RED }}-7.5 \rho_{\mathrm{BLUE}}+1}, $

(4)

$ I_{\mathrm{LSW}}=\frac{\rho_{\mathrm{NIR}}-\rho_{\mathrm{SWIR} 1}}{\rho_{\mathrm{NIR}}+\rho_{\mathrm{SWIR} 2}} \text {. } $

(5)

式中：I_SDW为Sentinel-1双极化水体指数，ρ_VV为Sentinel-1 VV波段反射率，ρ_VH为Sentinel-1 VH波段反射率，I_NDW为归一化水体指数，ρ_GREEN为Sentinel-2绿色波段反射率，ρ_NIR为Sentinel-2近红外波段反射率，I_NDV为归一化植被指数，ρ_RED为Sentinel-2红色波段反射率，I_EV为增强植被指数，ρ_BLUE为Sentinel-2蓝色波段反射率，I_LSW为地表水体指数，ρ_SWIR1和ρ_SWIR2为Sentinel-2短波红外波段反射率。

夏玉米遥感识别最佳时相一般在穗期—抽雄期(8月10—20日)^[28-29]，因此以8月中旬拼接的Sentinel-2影像为最佳时相。为了避免夏季多云和重访次数不均匀导致可用性数量不确定的问题，应为时间序列影像的融合指定适当的时间间隔，构建最佳时相、光学时序和光学与SAR融合时序及相关指数的3种数据集。时间序列步长应满足各期Sentinel-1/2影像覆盖整个研究区，经验证步长应以月为周期可满足需求，GEE云平台统计数据显示序列影像共包含5期数据(图 3)。

1.3.3 随机森林法与指标重要性评价

目前遥感监督分类常用方法有决策树法、最大似然值法、支持向量机法、随机森林法等^[30]。本研究通过GEE平台建立随机森林分类模型，进而提取研究区地类信息。随机森林法通过随机采样并放回地抽取的方式对样本进行重采样，其中约80%的样本数据作为袋内数据创建决策树，约20%的样本作为袋外数据进行模型验证，相对于其他方法，随机森林法具有更高的准确率且不会过度拟合数据^[31]。

特征重要性指的是特征指标对于分类结果的影响程度，能够判断输入指标对于结果的贡献。随机森林可以用来对特征重要性进行评估，主要分为2种方法：(1)平均不纯度减少(Mean Decrease Impurity，MDI)，通过统计节点分裂时不纯度的下降数值来衡量某个节点的重要性；(2)平均精确率减少(Mean Decrease Accuracy，MDA)，方法是随机交换袋外数据集中某个特征的值，然后重新进行预测，通过衡量分类或回归的准确度下降程度来计算该特征的重要性。本研究采用GEE随机森林法重要性评估方法对方案3中包含的相关特征重要性进行评估并统计制作柱状图，重要性越高，表明特征对结果的贡献值和重要性越大。

1.3.4 精度验证

混淆矩阵也称误差矩阵，是对实地验证的像元类型与遥感影像相对应分类结果的比较，通过混淆矩阵方法可以得到总体精度(overall accuracy，OA)[公式(6)]、生产者精度(producer accuracy，PA)[公式(7)]、用户精度(user accuracy，UA)[公式(8)]和Kappa系数(Kappa coefficient)[公式(9)]等精度衡量指标。

$ V_{\mathrm{OA}}=\frac{\sum\limits_{i=1}^n X_{i i}}{X} $

(6)

$ V_{\mathrm{PA}}=\frac{X_{i i}}{X_{i *}} $

(7)

$ V_{\mathrm{UA}}=\frac{X_{i i}}{X_{* i}} $

(8)

$ K=\frac{X \sum\limits_{i=1}^n X_{i i}-\sum\limits_{i=1}^r X_{i *} X_{* i}}{X^2-\sum\limits_{i=1}^n X_{i *}+X_{* i}} 。$

(9)

式中：n为类别总数；X_ii为正确分类的像元数量；X_i*和X_*i分别为第i行和第i列的像元数量；X为评估样本像元总数；V_OA为总体精度，指所有正确分类的土地覆盖类别的检验点数所占总抽取的检核点数的百分比；V_PA为生产者精度，表示在此次分类中，该类别的地面真实参考数据被正确分类的概率；V_UA为用户精度，表示在分类图上，落在该类别上的检验点，被正确分类为该类别的比率；K为Kappa系数，表示分类结果比随机分类好多少的指标，Kappa值的范围为0~1，Kappa值越大表示分类精度越高。

2 结果与分析 2.1 结果及各分类方法比较

基于夏玉米遥感识别最佳时相、夏玉米遥感识别光学时序、夏玉米遥感识别光学与SAR融合时序3种方案提取商河玉米和大豆玉米带状复合种植面积(图 4)。3个方案提取结果均表明夏玉米全县分布广泛，中部和西南地区城市、工业用地集中的区域玉米分布相对稀疏。大豆玉米带状复合种植区主要分布在西部和北部，分布相对集中，与实地调查及从当地农业部门获取的信息较为符合。2022年统计年鉴显示全县夏玉米播种6.20万hm²、大豆玉米带状复合种植播种620.00 hm²，占耕地总面积的84.03%。方案1—3夏玉米提取面积分别为7.33万、7.53万和6.90万hm²，大豆玉米带状复合种植面积分别为913.33、813.33和600.00 hm²。与2023年统计年鉴相比，方案3提取的玉米和大豆玉米带状复合种植面积结果与统计数据仅差1.03万hm²和20.00 hm²，精度分别达到86.96%和97.05%(表 3)。

利用验证数据集通过混淆矩阵得到分类精度评价(表 4)。所有分类方案的总体精度均大于80%，Kappa系数均大于0.77。方案1的总体精度最低，其分类特征仅包含光学卫星波段及衍生指数，OA和Kappa系数分别为83.01%、0.77。方案2在方案1的基础上加入了时序影像，OA和Kappa系数分别提高到89.44%、0.86。方案3的OA和Kappa系数最高，其加入SAR卫星波段及衍生指数，与方案1比较，OA和Kappa系数分别提高了9.5%、0.12，达到92.51%、0.89(表 4)。这表明方案3夏玉米遥感识别光学与SAR融合时序是研究区夏玉米、大豆玉米带状复合种植面积提取的最佳方案。

图 4鹰眼图区域展示3种方案分类结果与夏玉米遥感识别最佳时相RGB合成(R：B8，G：B4，B：B3)、6—8月SAR影像VV波段RGB合成(R：6月，G：7月，B：8月)和6—8月NDVI数据RGB合成(R：6月，G：7月，B：8月)进行对比分析，结果显示3种方案均包含区分夏玉米和大豆玉米带状复合种植的部分信息，时间序列数据和SAR数据挖掘了夏玉米和大豆玉米带状复合种植的其他信息，使得作物信息提取结果在图斑边缘和“椒盐现象”有一定改善，图斑边界平滑性明显提升。针对样例区域进行人工提取玉米、大豆玉米带状复合种植和其他作物面积的正确分类像元，正确分类像元验证各方案提取各作物面积的用户精度(表 5)。所有分类方案中玉米和其他的分类精度均较高，超过92%，大豆玉米带状复合种植的分类精度在时序特征和SAR特征加入后有明显提升。其中方案1的大豆玉米带状复合种植分类精度为71.81%，方案2分类精度为89.97%，方案3分类精度为95.07%。结果表明，夏玉米遥感识别最佳时相、夏玉米遥感识别光学时序、夏玉米遥感识别光学与SAR融合时序3种方案均能实现夏玉米的有效提取，且3种方案数据数量不断增加提升了夏玉米和大豆玉米带状复合种植面积提取精度。

2.2 分类贡献度

基于最优方案3，利用GEE计算12个分类特征在不同时段的重要性，并根据分类指标和时段进行统计。

按照分类特征进行特征重要性均值汇总，光谱分类特征中SWIR1、BLUE和NIR重要性相对较高，重要性分别达到322.06、316.52和316.15；衍生指数中EVI、NDWI和LSWI指数重要性相对较高，重要性分别达到312.82、310.91和310.53；光谱分类特征及其衍生指数重要性偏差不大，SAR卫星及其衍生指数重要性相对较低(图 5)。重要性均值反映的是各指标整体的贡献程度，不同分类特征重要性标准差反映各指标在不同月份贡献程度的差异性，其中标准差最小的是LSWI，其重要性指数也相对较高，表明LSWI指数在各月对分类结果的贡献程度均较大(图 6)。

按照不同时段进行特征重要性月度汇总，8月重要性为309.69。9月和10月重要性相对较低，分别为283.64和297.63(图 7)；不同分类时段内各指标重要性标准差反映各指标贡献程度的差异，不同时段重要性标准差最大月份为8月，最小为10月(图 8)。这说明8月对夏秋作物识别具有较高价值，其间月合成影像及其指数对夏季作物识别重要性的差异最大。

3 讨论

GEE云平台的发展有力推动了卫星遥感影像的数据挖掘、研究及应用，如土地利用分类和变化监测等诸多方面^{[7, 13, 31-34]}。传统遥感影像处理方法存在数据收集、存储、预处理繁琐和运算缓慢等问题，与传统方法相比，GEE云平台极大地简化了对遥感影像的预处理、筛选和运算过程。统计表明GEE云平台综合效率提升超过90%，如果将内存资源消耗、用户使用体验等因素加权在内，GEE云平台的综合效率会提升更多^[35-36]。在以往的作物面积信息提取工作中，在地面调查完成的情况下仍需要花费数天时间，占用大量存储和运算资源，本文使用GEE云平台则显著提升了效率。

夏季作物提取的研究主要集中于夏玉米，大豆玉米带状复合种植的情况较少受到关注，且现有的夏季作物提取研究多使用中低分辨率卫星数据，很少采用高分辨率图像，这使得破碎地块环境下存在大量混合像元，导致作物种植面积的错估，特别是种植面积较少的作物如大豆玉米带状复合种植的低估。本文基于高分辨率光学和SAR影像，在GEE云平台中设计实现了基于夏玉米遥感识别最佳时相、夏玉米遥感识别光学时序和夏玉米遥感识别光学与SAR融合时序的作物面积信息提取，夏玉米遥感识别最佳时相、夏玉米遥感识别光学时序和夏玉米遥感识别光学与SAR融合时序3种方案识别作物分布总体趋势一致，光学与SAR融合时序在作物面积信息提取精度上最为优异。

以往诸多农业遥感的研究均致力于探索解决作物最佳时相的选择，依靠其弱化其他因素的干扰，降低遥感信息中可能存在的误差，减少信息处理难度^[37-38]。夏玉米遥感识别最佳时相一般在穗期—抽雄期前^{[28-29, 39]}，本文以对应的8月中旬Sentinel-2 L2A与植被相关的可见光、近红外、红边波段为夏玉米遥感识别最佳时相数据^{[6, 28, 34]}，结果实现了80%以上的总体精度。Sentinel-2空间分辨率最高达10 m，相对于MODIS、Landsat8/9数据250 m和30 m的分辨率，较好分类依赖于其能够更为精细地识别破碎地块及减少混合像元。夏玉米遥感识别光学时序通过夏玉米发育期内构建每月Sentinel-2夏玉米遥感识别最佳时相用于作物面积提取，作物信息提取结果显示基于时序特征的分类结果比夏玉米遥感识别最佳时相分类结果总体精度高6.43%以上，尽管夏玉米统计值仅差3.16%(方案2)和1.69%(方案3)，但图斑存在的“椒盐现象”得到了改正，地块连续性、完整以及边界平滑性明显提升。各分类指标对分类结果重要性上相差不大，但各指标在不同时段中重要性有较大差异，这可能与作物不同发育期的作物形态等特性有关，不同时段重要性稳定的指标对作物发育早期作物信息的提取有较大贡献，通过这些指标可能实现作物在发育早期的面积提取。夏玉米遥感识别光学与SAR融合时序在夏玉米遥感识别光学时序加入Sentinel-1A GRD相关波段和SDWI指数，基于光学与SAR融合时序特征的提取结果显示，大豆玉米带状复合种植面积仅比统计值偏小20 hm²，表明具有时序特征的分类指标特别在玉米及大豆玉米带状复合种植区精准识别上具有优势。在作物种植结构复杂的夏季，“同谱异物”现象导致作物光谱信息差别不显著^[29]，无法将作物特别是玉米及大豆玉米带状复合种植区分开，SAR及相关指标为研究区的夏玉米和大豆玉米带状复合种植的区分提供额外的信息，如土壤水分、植株高度等信息，有助于提高分类模型的作物分类能力^{[22, 40-41]}。但重要性评估结果SAR波段重要性较低，其原因一方面可能是SAR数据的噪声难以避免，另一方面可能是未充分利用双极化Sentinel-1数据中包含的后向散射系数等极化信息^[20]。

主要存在2方面不足：一是按月构建时间序列可能存在多云导致的无数据现象，可能存在的数据缺失问题值得进一步研究，如Sentinel-2后续卫星如期发射，将极大改善数据不足的现象；二是分类指标选择依赖于文献报道，仅选择了2个卫星的8个波段及衍生的5个指数，未考虑其他波段、指数及其他数据的作用，如海拔、坡度、后向散射系数和目标极化分解参数等，有待进一步探索。

4 结论

基于高分辨率光学和SAR影像，在GEE云平台中设计实现了基于夏玉米遥感识别最佳时相、夏玉米遥感识别光学时序和夏玉米遥感识别光学与SAR融合时序的作物面积信息提取，光学与SAR融合时序的分类结果总体准确率(92.51%)高于其他2种方案(83.01%和89.44%)，结果表明夏玉米遥感识别最佳时相、夏玉米遥感识别光学时序和夏玉米遥感识别光学与SAR融合时序3种方案识别作物分布总体趋势一致，光学与SAR融合时序在作物面积信息提取精度上最为优异。后续拟基于GEE进一步优化多源数据融合与时序分析指标等，提高作物识别和提取能力。