基于GBDT和Google Earth Engine的冬小麦种植结构提取

张海洋， 张 瑶， 田泽众， 吴江梅， 李民赞， 刘凯迪

中国农业大学智慧农业系统集成研究教育部重点实验室， 北京 100083

引 言

农作物种植结构是某一地区或生产单位农作物种植类型、 种植面积以及空间分布的综合反映， 是表征农业生产资源利用科学性、 合理性的重要指标。 冬小麦是世界三大粮食作物之一， 实时、 准确地获取冬小麦精确的种植结构信息对作物产量估算、 粮食政策制定、 粮食贸易流通和国家粮食安全保障具有十分重要意义[1-2]。

传统的逐级上报和抽样调查方法， 具有一定的主观因素， 耗时费力、 成本高以及时效性差。 20世纪70年代以来， 随着RS 和GIS 技术的快速发展以及影像数据源的不断丰富， 光学遥感观测技术在大尺度作物种植面积和空间分布信息的监测上已被广泛应用。 当前主要使用AVHRR， MODIS， Landsat (TM， ETM+和OLI)， SPOT等影像为数据源， 然而中国农业区的耕地经营分散、 农业景观破碎和作物种植结构复杂， 致使上述影像难以有效捕捉田块尺度的冬小麦等作物对象， 大大降低了农作物的分类准确性。 与上述卫星传感器相比， Sentinel-2卫星传感器在保证相对较高的空间分辨率和时间分辨率的同时还提供了丰富的红边信息， 可为冬小麦种植结构研究提供有力的数据支撑[3]。 谷歌地球引擎(Google Earth Engine， GEE)是由谷歌、 卡内基梅隆大学和美国地质调查局联合开发的基于云计算的开源地理信息处理平台[4]。 该平台向用户免费提供海量的卫星影像和其他地球观测数据， 而且借助Google的高性能集群服务器实现影像在线可视化处理， 具有强大的数据处理能力， 可大大提高工作效率。 目前， GEE 被成功应用于作物长势监测、 土地利用类型分类等方面[5]。

近年来， 随着机器学习分类算法的发展， 利用机器学习分类算法进行作物种植结构提取得到了国内外学者的广泛青睐， 如朴素贝叶斯、 支持向量机、 随机森林等分类器。 其中随机森林分类器是一种基于Bagging思想的并行集成基学习器， 该方法训练和预测速度快、 计算成本低， 被广泛应用于地面作物分类识别。 Jin等[6]等利用随机森林分类器构建了作物/非作物和玉米/非玉米识别模型， 作物/非作物识别准确率为85%， 玉米/非玉米在坦桑尼亚的准确率为79%， 在肯尼亚为63%。 Chong等[7]使用随机森林分类器进行黑龙江作物分布制图， 总体分类准确率为89.75%。 RF分类器是通过减少模型方差提高性能， 容易陷入过拟合， 而且， 当训练样本较少时， 地物分类效果会显著降低。 与随机森林分类器相比， GBDT分类器是基于权值的基学习器的串行集成模型， 通过不断减少模型在训练过程中产生的残差实现样本数据的分类， 利用该方式能大大增强对数据的拟合能力以及提高分类准确性， 可为提升地物分类精度提供一种有益方法。

利用以上机器学习方法提取农作物种植结构时， 多特征提取是目前的研究热点之一。 针对中国复杂的农作物种植结构， 传统的仅依赖光谱波段信息或光谱指数信息往往难以保证可靠的地物识别精度。 Shetty等[8]利用Landsat8 OLI遥感卫星的光谱波段(Blue， Green， Red和NIR)特征和归一化植被指数(NDVI)特征， 对研究区进行地物分类， 但由于采用的特征变量较为单一， 地物识别总体分类准确率较低。 针对农作物分类中依赖有限特征无法取得较高的分类精度问题， 学者们尝试引入纹理特征输入机器学习模型， 已被证实在提高作物分类精度方面具有一定优势， 对进一步提高农作物空间分布提取结果精度发挥了重要作用。 Khosravi等[9]提取加拿大温尼伯市附近的地物信息， 在光谱特征、 植被指数特征和偏振特征的基础上加上纹理特征， 可以显著提高每个分类类别的分类精度(Kappa系数超过85%)。 此外， 一些研究表明地形特征也可有效提高土地利用分类精度[10]。 目前， 研究者们主要采用机器学习方法， 提取研究区多特征信息进行大尺度作物分类研究， 但针对小尺度、 高精度的田块尺度的作物分类效果仍有待进一步探索。 此外， 冬小麦在不同生长生育期有不同的形态特征， 同期生长的其他作物也随之变化[11]。 如何找到冬小麦种植结构提取的最优生育物候期， 进一步提高冬小麦识别精度值得进一步分析与探究。

基于此， 利用GEE云平台， 选用10 m分辨率的Sentinel-2遥感数据， 充分利用研究区光谱波段特征、 光谱指数特征、 纹理特征和地形特征等多维特征， 探究冬小麦识别的最佳生育期， 对比分析不同分类模型在田块尺度下的冬小麦识别性能， 旨在提出一种基于GEE云平台的冬小麦种植结构提取方法。

1 研究区概况

研究区位于河南省新乡市封丘县陈固镇(如图1所示)， 覆盖范围处于35°5′39.82″N—35°11′22.43″和114°15′45.62″E—114°23′44.77″E之间。 全镇总面积6.1千公顷， 耕地面积约为4.53千公顷， 下辖23个行政村。 该地区地势平坦， 土壤的质地主要有壤土和粘土， 且壤土占比居多。 该区域属暖温带大陆性季风气候， 夏季炎热多雨， 冬季寒冷干燥， 年平均气温14.5 ℃， 1月和7月平均气温分别约为-1.0和27.2 ℃， 年降雨量在600 mm左右。

图1 研究区域位置图和地面样本点分布示意图Fig.1 Location of the study area and the distributions of ground sample sites

该研究区是重要的粮食产地， 以冬小麦、 夏玉米一年二熟的作物轮作模式为主， 春季的主要农作物为冬小麦、 金银花和大蒜等。 此外， 亦有少许农户种植梨、 苹果等果树和林地。 经过实地走访调查， 该地区的冬小麦、 大蒜和金银花的生育期如表1所示， 其中冬小麦的生育期有播种期、 出苗期、 分蘖期、 越冬期、 返青期、 起身拔节期、 抽穗扬花期和灌浆乳熟期， 播种期约为每年的10月上旬， 收获期为下一年的6月上旬； 该地区种植的大蒜为秋播大蒜， 其生长周期主要有播种期、 萌芽期、 幼苗期、 花芽鳞芽分化期、 抽薹期和成熟期， 种植期为每年的9月下旬或10月上旬， 收蒜头期为第二年的5月下旬和6月上旬； 金银花的生长发育期主要为萌芽期、 新梢生长期、 现蕾期、 开花期、 缓慢生长期和越冬期， 金银花种植时间不固定， 一般在2月进行插种。 金银花的开花期为每年的5月中旬至9月下旬， 可在开花期内进行4次采收。

表1 研究区内冬小麦、 大蒜和金银花的生长发育时期Table 1 Growth and development periods of winter wheat, garlic and honeysuckle in the study area

依据冬小麦生长特点， 起身拔节期、 抽穗扬花期和灌浆乳熟期是冬小麦的关键物候期， 也是对冬小麦进行遥感识别的最佳时期， 因此， 分别对这三个时期的Sentinel-2影像进行提取， 以识别冬小麦的空间分布， 并探寻哪一个时期是冬小麦种植结构的最佳提取时期。

2 实验部分

2.1 数据来源

研究所用的数据包括Sentinel-2影像数据、 SRTM高程数据和样本数据， 详细介绍如下：

2.1.1 Sentinel-2卫星影像

选用的遥感数据为Sentinel-2卫星影像数据。 Sentinel-2卫星主要由Sentinel-2A和Sentinel-2B两颗卫星组成， 拥有较高的时间分辨率(重访周期为5 d)和空间分辨率(最高为10 m)。 Sentinel-2卫星搭载的多光谱传感器主要有13个光谱波段， 包括10， 20和60 m三种空间分辨率， 其中10 m分辨率的波段有红波段(B4)、 绿波段(B3)、 蓝波段(B2)和宽近红外波段(B8)； 20 m分辨率的波段为3个红边波段(B5、 B6和B7)、 窄近红外波段(B8A)、 短波红外(B11和B12)； 60 m分辨率的波段有海岸波段(B1)、 水蒸汽波段(B9)和卷云波段(B10)。 GEE平台提供了Level-1C 和Level-2A 两种处理级别的Sentinel-2 MSI数据， Level-1C 是经过了辐射定标、 几何校正(包括空间配准和正射校正)的大气顶层反射率(top of atmosphere reflectance， TOA)产品， Level-2A是在Level-1C的基础上经过大气校正的地表反射率(surface reflectance， SR) 产品， 使用的是Level-2A级的产品。

基于GEE平台首先获取研究区冬小麦三个关键生育期内的所有Sentinel-2 L2A级MSI影像； 然后利用属性“CLOUDY_PIXEL_PERCENTAGE”筛选出云量低于10%的影像， 接着对影像进行取中值处理， 继而将各生育期内的影像合成一幅影像， 最后对各生育期内的影像应用三次卷积法将低分辨率的波段重采样为10 m。

2.1.2 其他辅助数据

SRTM (shuttle radar topography mission)是美国太空总署(NASA)和国防部国家测绘局(NIMA)以及德国与意大利航天机构合作， 由美国发射的“奋进”号航天飞机上搭载SRTM系统完成测量。 SRTM V3.0是2016年最新发布的全球高程数据集， 其中SRTMGL1为全世界1弧秒数据， 采用WGS-84投影， 空间分辨率约为30 m， 在GEE平台上通过搜索“SRTMGL1_003”获取[24]。

2.1.3 样本数据

依据研究区实际地物类型的分布情况， 将研究区分为冬小麦(winter wheat， Ww)、 建筑物(building， Bu)、 其他植被(other vegetation， Ov)和道路(road， Ro)4类地物， 其他植被主要为金银花、 大蒜和少量种植的果树与林地等。 样本数据的采集方式是目视解译标定和实地采样， 在采集过程中保证了样本均匀分布和随机性。 其中， 建筑用地和道路是利用Google Earth高分遥感影像进行目视解译标定的， 建筑用地240个， 道路190个。 冬小麦和其他植被是于2021年6月1日至6月15日利用野外调查方式进行标定的， 冬小麦的标定样本个数为260个， 其他植被共获得个200实地样本点， 包含金银花123个、 大蒜53个和果树与林地24个， 每个采样点均利用定位设备获取其经维度坐标并记录相应的植被类型等信息， 且选取样本点对应种植植被的面积均超过0.2公顷， 样本点具体分布如图1所示。 定位设备是采用千寻位置网络有限公司的千寻星矩SR6网络RTK接收机和千寻知寸技术服务(Qianxunxingju SR6, Qianxun Spatial Intelligence Inc., Shanghai, China)。

2.2 分类特征构建

特征变量的选取是地物分类识别的重要前提之一。 优化特征变量及其组合可以有效提高遥感地物解译的分类精度[12]。

从研究区的分布特征出发， 充分考虑冬小麦种植区域与其他土地类型的显着差异， 如明显不同于其他地物的光谱特征和纹理特征， 以及地形特征等。 基于此， 选取有益于提高区域冬小麦提取精度的光谱、 纹理、 地形等特征波段， 构建基于Google Earth Engine的冬小麦种植结构识别系统。

影像的光谱特征是地物遥感解译的关键特征变量。 本研究选择光谱波段特征为Sentinel-2影像输出的波段。 光谱指数特征有NDVI， NDWI， NDBI， CIRE和S2REP， 其中NDVI反映植被生长状态和植被茂密程度的指标； NDWI可有效抑制其他类型的地表类型而凸显水体信息； NDBI将城镇灰度值增高， 降低其他地类值； CIRE和S2REP是Sentinel-2影像特有的红边波段指数特征， 红边波段数据可为农作物类型识别提供了强有力的数据支持。

与此同时， 也有研究表明， 由于遥感影像上的同一光谱实际上可能是不同的地物， 而同一地物在遥感影像上也可能具有不同的光谱特征。 因此， 选择单一光谱特征进行土地类型遥感分类可能会导致遥感解译识别结果出现部分误差和遗漏。 纹理特征代表图像灰度的空间变化和重复， 或图像中重复的局部图案和排列规则， 可以在一定程度上提高遥感的分类精度。 选用3×3移动窗口， 利用灰度共生矩阵(gray level co-occurrence matrix， GLCM)计算纹理特征。 由于近红外波段对植被更敏感， 因此选取Sentinel-2影像的近红外波段(B8)计算影像的纹理特征， 共得到18个纹理特征参数。 从该地区冬小麦的纹理特征出发， 综合考虑纹理参数之间的相关性、 差异性和冗余性， 从对比度、 相关性、 熵等方面选取了最常见的4种纹理参数进行构造特征参数并训练分类器， 以减少过多纹理特征之间的数据重叠和冗余。 选取的纹理特征包括角二阶距(ASM)、 对比度(CONTRAST)、 相关性(CORR)、 信息熵(ENT)。

此外， 使用GEE自带的空间分辨率为30 m的地形数据SRTMGL1_003， 通过ee.Terrain.products(input)函数计算高程和坡度两个参数。 然后将它们作为两个独立的特征带添加到合成的多波段影像中， 用于地物的遥感识别。 该数据空间分辨率为30 m， 利用三次卷积内插法将数据重采样为10m分辨率。

选取14个光谱特征(9个光谱波段特征和5个光谱指数特征)、 4个纹理特征、 2个地形特征构建冬小麦种植区域识别特征参数， 进而实现对研究区冬小麦的有效识别。 各特征变量的名称、 描述及其计算公式如表2所示。

表2 选择的特征变量及其计算公式Table 2 Selected characteristic variables and their calculation formulas

2.3 分类方法

梯度提升决策树(gradient boosting decision tree， GBDT)分类器是由Friedman提出的一种基于Boosting策略的集成机器学习算法， 其基本思想是训练多个弱分类器不断提升性能， 除了第1棵决策树采用原始预测指标生成外， 每一轮迭代中的目标都是令当前学习器的损失函数最小化， 即令损失函数总是沿着其梯度方向下降， 通过不断迭代使最终残差趋近于0， 将所有树的结果累加起来便可得到最终的预测结果。 其中， 弱分类器一般选择CART TREE(即分类回归树)。 GBDT分类器的具体实现过程如下：

(1)初始化弱学习器

(1)

式(1)中，e为经最小二乘法划分结点后的叶子结点取值。

(2)对每棵树m=1, 2, …,M， 对每个样本i=1, 2, …,N， 计算负梯度。 即残差

(2)

式(2)中，f(xi)为弱学习器的预测值，yi为弱学习器的真实值。

将得到的残差作为新样本真实值， 数据(xi,rim),i=1, 2, …,N作为下棵树的训练数据， 得到新的回归树fm(x)， 其对应的叶子结点区域为Rjm,j=1, 2, …,J为回归树t的叶子结点个数， 对叶子区域j=1, 2, …,J， 计算最佳拟合值， 有

(3)

更新强学习器， 则有

(4)

得到最终学习器

(5)

2.4 种植结构提取流程

基于GEE云平台的冬小麦种植结构提取详细过程如下， 其具体流程如图2所示。

图2 冬小麦种植结构提取流程图Fig.2 Extraction flowchart of winter wheat planting structure

(1)分别获取2021年冬小麦起身拔节期、 抽穗扬花期和灌浆乳熟期等三个生育期内所有云量低于10%的Sentinel-2 L2A级影像， 然后对各生育期内的影像进行影像裁剪、 镶嵌、 去云、 求中值和重采样等操作。 获取研究区的DEM高程数据， 并进行裁剪和镶嵌等处理。

(2)获取研究区冬小麦、 建筑物、 其他植被以及道路等4类地物的地面样本点数据， 记录样本点的地物种类和位置。 并将各类地物样本点按照8: 2随机划分为训练集和测试集。

(3)提取各生育期内合成影像的光谱波段特征、 光谱指数特征、 纹理特征以及地形特征等分类特征变量， 将各生育期的所有分类特征变量以及地面样本点数据的训练集输入GBDT分类器中， 设定不同生育期GBDT分类器树的个数， 得到地物分类结果。

(4)将地面样本点测试集的特征变量代入训练后的GBDT分类器中， 获得不同生育期内研究区地物的分类准确率， 然后利用最优的地物分类结果提取研究区冬小麦的种植结构。

2.5 分类精度评价指标

为有效评估分类算法对研究区不同覆盖类型的提取精度， 通过构建混淆矩阵(confusion matrix)对分类结果进行精度度量。 基于混淆矩阵， 可以计算用户精度(user accuracy， UA)、 生产者精度(producer accuracy， PA)、 总体分类准确率(overall accuracy， OA)和Kappa系数(Kappa)等分类性能评价指标。 用户精度反映各地物被正确预测的比例， 生产者精度表示各地物被正确识别为该地物的比例， 总体分类准确率和Kappa系数是根据所有待评估地物类别的漏分和错分情况给出的一种更为全面的分类准确性评估指标， 反映提取结果与真实地物空间分布的一致性， 这四种度量指标的计算公式如式(6)—式(9)

(6)

(7)

(8)

(9)

式中，N为测试集的数量；m是混淆矩阵的行列数， 代表地物的类别数；xii为混淆矩阵对角线上的值， 表示被分到正确类别的样本数；x+i为混淆矩阵第i列元素相加之和， 表示第i类的真实样本数；xi+为混淆矩阵第i行元素相加之和， 表示预测为第i类的样本数。

3 结果与讨论

3.1 冬小麦不同生育期的地物分类结果

冬小麦不同生育期内， 研究区地物物候特征、 空间分布特征和光谱特征不同， 冬小麦不同生育期内的地物识别效果也会有差异。 利用GBDT分类器对研究区冬小麦起身拔节、 抽穗扬花和灌浆乳熟等三个时期内的地物类型进行分类识别， 分别设置GBDT分类器树的数量为1 500， 1 500和1 800棵， 地物识别效果的混淆矩阵图如图3所示。 由图3可知， 186个测试样本中， GBDT分类器在冬小麦起身拔节期的测试样本中有158个样本被分类正确。 其中， 建筑物分类准确性最高， 其他植被和道路容易混淆(分别占其他植被观测样本总数的11.54%、 道路观测样本总数的6.82%)。

图3 冬小麦不同关键物候期地物分类结果的混淆矩阵图(a)： 起身拔节期的地物分类混淆矩阵结果； (b)： 抽穗扬花期的地物分类混淆矩阵结果； (c)灌浆乳熟期的地物分类混淆矩阵结果Fig.3 Confusion matrixes of ground classification results in different key phenological periods of winter wheat(a)： Standing and jointing stage; (b)： Heading and flowering stage; (c)： Grain-filling and milky stage

抽穗扬花期有181个测试样本被正确分类， 建筑物和冬小麦的分类精度较高， 其他植被易被误判为道路(占其他植被观测样本总数的7.69%)； 灌浆乳熟期有144个测试样本被分类正确， 在冬小麦3个关键生育期内分类精度最差。

为了对冬小麦3个生育期内地物识别效果进行更深入的对比分析， 采用指标UA， PA， OA和Kappa系数对此3种分类结果进行评价， 如表3所示。 由表3可知， 在冬小麦三个关键生育期中， 抽穗扬花期的地物识别精度最高， OA和Kappa系数分别为97.01%和95.92%， 比起身拔节期分别高出2.40和3.24个百分点。 灌浆乳熟期分类效果最差， OA仅为86.23%， Kappa系数为81.33%。 说明在小田块条件下， 基于GBDT分类器可以在冬小麦抽穗扬花期有效提取土地覆盖信息， 具有很好的地物分类识别效果。

表3 冬小麦不同关键物候期地物分类结果Table 3 Classification results of ground in different key phenological stages of winter wheat

对于UA和PA， 抽穗扬花期的四种地物分类结果均最高(起身拔节期和抽穗扬花期的建筑物分类效果一致， UA均为97.78%， PA均为100%)。 就冬小麦的三个关键生育期而言， 建筑物和冬小麦的UA和PA分类结果均超过90%。 在灌浆乳熟期， 由于该时期作物叶片受脱落酸的影响， 叶绿素会分解， 胡萝卜素， 类胡萝卜素， 花青素等大量合成， 其他植被和道路易相互混淆， UA和PA的值最低(其他植被的UA和PA分别为66.67%和69.23%， 道路的UA和PA分别为79.07%和86.23%)， 起身拔节期减少了其他植被和道路之间的相互混淆， 其他植被和道路的UA和PA值均超过灌浆乳熟期(其他植被的UA和PA均为84.62%， 道路的UA和PA分别为93.02%和90.91%)。

通过直观目视解译并对比原始影像可发现， 冬小麦的三个关键生育期分类制图结果整体效果均较好(图4)。 从图4中可以看出， 冬小麦和道路覆盖区域轮廓清晰、 形状基本一致； 建筑物分布连续、 边界分明； 其他植被提取较为完整。

图4 研究区地物分类结果(a)： 研究区的Sentinel-2影像(RGB)， 影像时间是2021-05-02; (b)： 冬小麦起身拔节期的地物识别结果； (c)： 冬小麦抽穗扬花期的地物识别结果； (d)： 冬小麦灌浆乳熟期的地物识别结果Fig.4 Results of ground classification in the study area(a)： Sentinel-2 image (RGB) in the study area, the image time is 2021-05-02; (b)： Ground recognition results of winter wheat in standing and jointing stage; (c)： Results of ground feature recognition in heading and flowering stage of winter wheat; (d)： Ground recognition results of winter wheat in grain-filling and milky stage

为了进一步对比冬小麦的三个关键生育期内地物识别效果， 选择了研究区中3个代表性的位置区域进行了8倍放大， 如图5所示。 从图5中可以看出， 起身拔节期的建筑物分类结果较实际建筑物偏大， 也易将冬小麦地块错分为道路； 灌浆乳熟期的建筑物容易漏分， 尽管抽穗扬花期的分类结果与实际地物会稍稍有些偏差， 但相较其他两个时期， 该时期的识别效果表现最好。

图5 在冬小麦三个关键生育期内分类结果的局部放大图Fig.5 Partial enlarged views of classification results in three key growth stages of winter wheat

3.2 特征贡献率分析

GBDT分类器可以对特征变量的重要程度和贡献率进行分析评价， 图6是GBDT分类器在冬小麦3个关键生育期的特征重要性分布， 重要性得分越高， 说明该特征变量对分类结果的影响和贡献就越大。 由图6可知， 在冬小麦的起身拔节期和抽穗扬花期， 纹理特征(CONTRAST)、 红边指数特征(S2REP)、 红边波段特征(B7)以及蓝波段特征(B2)的排序靠前。 而在灌浆乳熟期， 短波红外波段特征(B12)、 建筑特征(NDBI)和纹理特征(CONTRAST)的贡献率较大。 这可能是因为在起身拔节期和抽穗扬花期， 研究区的农作物处于生长茂盛期， 而红边波段(B7)和红边指数(S2REP)特征在这些生育期对于农作物的探测和识别较敏感， 因此这些特征对分类结果的影响较高； 灌浆乳熟期的农作物逐渐成熟， 植被覆盖度和叶绿素含量降低， 该时期的短波红外波段特征(B12)和建筑特征(NDBI)重要程度占比较高。 其中Sentinel-2 的B12特征的重要性得分最高， 对研究区土地利用分类贡献度最大， 这是因为该时期其他作物(金银花、 果树和林地)处于生长旺季， 含水量高， 冬小麦作物处于生长后期， 干物质增加， 叶绿素和水分含量会降低， 短波红外波段对叶片含水量反映敏感， 能够很好地区分出其他植被和冬小麦。

进门后，川矢笑容可掬地对百里香叽哩呱啦地说了一通，还比比划划，夸张地晃着大拇指，百里香一句也没听懂。身着鬼子军服的庄翻译说：“川矢队长说，他非常欣赏中国民间的美食文化，也非常喜欢百里香肠铺的香肠，愿拜阁下为师，学习中国厨艺，让中国民间的饮食文化融入大东亚共荣圈，使阁下成为大东亚乃至世界顶级的厨艺大师！冒昧相求，万望阁下恩准！”

图6 各个冬小麦关键生育期GBDT分类器中20个变量特征重要性得分(a)： 起身拔节期特征变量重要性得分； (b)： 抽穗扬花期特征变量重要性得分； (c)： 灌浆乳熟期特征变量重要性得分Fig.6 Feature importance scores of 20 variables in GBDT classifier at key growth stages of winter wheat(a)： Standing and jointing stage; (b)： Heading and flowering stage; (c)： Grain-filling and milky stage

在冬小麦三个关键生育期内， 纹理特征(CONTRAST)对于冬小麦、 建筑物和道路等形状轮廓明显的地方的提取发挥了重要的作用， 因此CONTRAST的分类贡献率较高。 地形特征(ASPECT和SLOPE)的重要性得分均较低， 原因可能是该研究区地势较为平坦， 而坡度和坡向是基于地形分析得到， 因此参与分类的贡献度较低。

3.3 GBDT分类器树的数量对结果的影响

GBDT分类器训练过程中通过多轮迭代， 每轮迭代产生一个弱分类器， 每个弱分类器是在上一轮弱分类器的残差基础上进行训练， 属于串行生成， 弱分类器(树)的数量会极大影响GBDT分类器的计算速度和计算成本。 因此需要找到最佳的树的数量， 既能保证GBDT分类器的分类精度， 又能保证分类器的运行速度和计算成本。

为了探寻冬小麦三个关键生育期内利用GBDT分类器对地物进行分类时最优树的数量， 选取树的数量范围为100～2 000， 步长设置为100， 以总体分类准确率作为评价指标， 地物分类结果如图7所示。 从图7中可知， 随着树的数量递增， 冬小麦抽穗扬花期的地物识别精度逐渐升高， 当树的数量升高至1 500棵时， 地物识别效果达到最高， 为0.970 1， 之后随着树的数量增加， 研究区地物识别精度保持不变； 在灌浆乳熟期， GBDT分类器书的数量达到1 800棵时， 地物分类结果达到最高， 为0.862 3， 然后随着树的数量的增加， 地物识别结果和抽穗扬花期的一致， 保持不变； 在起身拔节期， 伴随GBDT分类器树的数量增加， 总体分类准确率在树的个数为1 400， 1 500和1 600棵时会出现一个峰谷。 因此， 利用GBDT分类器对地物进行分类时冬小麦起身拔节期、 抽穗扬花期和灌浆乳熟期最佳树的数量分别是1 400， 1 500和1 800棵。

图7 GBDT分类器树的数量不同时的地物分类结果Fig.7 Ground feature classification results with different numbers of GBDT classifier trees

3.4 不同分类方法的结果对比

为了验证GBDT分类器的分类效果和优适性， 选取了随机森林(random forest， RF)、 CART(classification and regression tree)[31]和朴素贝叶斯(Naive Bayesian， NB)[32]等分类器对研究区冬小麦抽穗扬花期的地物分类结果作为对比分析。 此三种分类器均在GEE环境下运行， 使用地物样本点的训练集进行训练， 将测试集代入训练好的分类器中进行验证和精度评价， 经过大量测试， 随机森林树的个数设置为500棵， CART和朴素贝叶斯分类器的参数为默认参数， 此三种分类算法分类结果的混淆矩阵如图8所示。

由图8可知， RF， CART和NB分类器在冬小麦灌浆乳熟期分别有160， 152和142个测试样本被识别正确。 RF， CART和分类器在建筑物测试样本中分别有43， 41和44个分类正确， 在其他植被测试样本中分别有24， 23和11个分类正确， 在道路测试样本中分别有40， 35和33个分类正确， 在冬小麦测试样本中的分类正确个数一致， 均为53个。 其中， CART和NB分类器的其他植被和道路容易混淆(分别占其他植被观测总数的11.54%和23.08%、 道路观测总数的18.18%和18.18%)， NB分类器易将其他植被错分为建筑物(占其他植被观测总数的30.77%)。

图8 不同分类器地物分类结果的混淆矩阵图(a)： RF分类器的地物分类混淆矩阵结果； (b)： CART分类器的地物分类混淆矩阵结果； (c)： NB分类器的地物分类混淆矩阵结果Fig.8 Confusion matrixes of ground classification results of different classifiers(a)： RF classifier; (b)： CART classifier; (c)： NB classifier

为了评定分类器的分类性能， 利用指标UA， PA， OA和Kappa对此3种模型分类结果进行进一步的分析和评价， 加上本工作的GBDT分类器， 4种分类分类器的分类结果如表4所示。 由表4可知， 在冬小麦的抽穗扬花期， GBDT分类器的分类精度最高， 总体分类准确率比随机森林分类器和CART分类器分别提高了1.20%和5.99%， Kappa系数比随机森林分类器和CART分类器分别提高了1.61%和8.04%， 朴素贝叶斯分类器的识别效果最差， 总体分类准确率和Kappa系数分别为84.43%和78.69%。

表4 不同分类模型在冬小麦抽穗扬花期的地物分类结果Table 4 Ground classification results of different classification models at heading and flowering stage of winter wheat

对于UA， 建筑物、 冬小麦和道路类别在4种分类器均获得了可观的结果(UA值均超过80%)。 GBDT在其他作物类别的值最高(95.83%)， 但CART分类器和NB分类器的UA值较低(CART分类器和NB分类器分别为67.65%和57.89%)， 究其原因可能是CART分类器尽管方法简单， 生成的分类树很直观， 但容易过拟合， 导致泛化能力不强。 NB分类器虽然有着坚实的数学基础， 但不能学习特征间的相互作用， 而且该分类器通过先验和数据来决定后验的概率从而决定分类， 所以分类决策存在一定的错误率。

图9 不同分类器的地物分类结果(a)： RF分类器的分类结果； (b)： CART分类器的分类结果； (c)： NB分类器的分类结果Fig.9 Ground classification results of different classifiers(a)： RF classifier; (b)： CART classifier; (c)： NB classifier

尽管已有10 m分辨率全国土地利用产品[13]， 但该产品是以全国范围尺度为研究区域， 存在多样性、 复杂性以及影像的时相无法统一的问题， 此外， 研究区域越大其局部区域纹理细节越难兼顾， 降低了作物提取精度。 对比其他学者的地物分类研究成果， 在空间分辨率一致的情况下， 本研究提出的方法能够实现略高于Belgiu等[14]提出TWDTW(time-weighted dynamic time warping)算法分类精度， 远高于Kristof等[15]利用RF分类器进行地物识别的准确度。

但本研究也存在一定的缺陷， 首先本研究是在降低地物分类类别的基础上提高冬小麦的分类准确率， 其次本研究的研究区域范围较小。 下一步研究中将增加作物分类种类和扩大研究区域对本研究提出的方法进行进一步的验证。

表5 其他冬小麦种植结构提取研究结果Table 5 Extraction results of planting structure of other winter wheat

4 结 论

针对中国农业耕地多呈现规模小、 破碎分散的特点， 基于GEE云平台开发了一种冬小麦种植结构提取模型， 能够实现10m空间分辨率的精细提取。

(1)对比冬小麦起身拔节期、 抽穗扬花期和灌浆乳熟期等三个关键生育期， 抽穗扬花期的分类效果最好(OA: 97.01%, Kappa: 95.52%)更接近研究区地表真实情况。

(2)在所有的光谱指数特征、 光谱波段特征、 纹理特征和地形特征中， 纹理特征中的影响最高， 其次是光谱波段和光谱指数特征， 地形特征的贡献率最少。

(3)在冬小麦的抽穗扬花期， 与RF， CART， NB等分类器相比， GBDT分类器的分类效果最佳， 其次是RF分类器， NB的地物识别效果最差。