潘 竞,艾尼玩·艾买尔,阿斯娅·曼力克,阿 仁,田 聪,马海燕,孙宗玖,郑逢令*
(1.新疆农业大学草业学院,乌鲁木齐 830052;2.新疆畜牧科学院草业研究所,乌鲁木齐 830057;3.新疆畜牧科学院天山北坡草地生态环境野外定位观测研究站,乌鲁木齐 830057)
小麦是现今世界上最主要的经济作物,在全世界其重要性仅次于稻米,同时也是我国主要的作物之一,在中国粮食生产和保障粮食安全方面占有举足轻重的地位[1]。因此开展小麦耕地的识别有助于政府有关部门指导农业生产活动,估计产量,优化调整农业结构,确保粮食安全[2]。对于农作物识别及长势信息采集的传统农业调查法缺少精确的作物空间分布信息,调查数据的质量往往易受到人为因素的影响[3],因而缺乏客观性。遥感能够在广泛的地理区域获取数据,可以为监测小麦的种植信息提供快速又准确的数据源。田间秸秆废弃物是指粮食作物收割后残留在田里的残余物料,它对耕地生态系统具有很大危害。新疆小麦、玉米、棉花等的秸秆废弃物,都是草食家畜饲草料重要的来源。正确评价可利用秸秆资源量与分布是秸秆杂物资源化综合利用的基础。在农业生产规划或管理决策过程中,草谷比是估算不同种类、不同地域大田作物秸秆资源量的重要且较快捷可行的方式[4],遥感等空间技术有快速在大范围实现作物分类并估算作物秸秆产量的潜力,但国内相关研究报道较少,国外利用遥感和GIS 等技术调查秸秆资源主要是服务于生物质能源的开发利用[5,6]。
遥感应用的各种卫星图像采集和数据挖掘方法相对比较成熟,能够有效地服务于多种农业应用场景,如产量预测、作物类型识别、杂草分类等[7~9]。近年来高分辨率卫星可以在较小的空间尺度上获得更高的空间和时间分辨率数据。遥感数据可以通过光学、雷达、红外等各种传感器采集[10],使得遥感卫星能够采集的数据类型越来越丰富。遥感图像处理技术不断有新突破,分类方法朝着自动化、智能化方向发展。Gulsen等[11]提出一种线性SVM和支持向量的混合模型可以有效地实现遥感图像分类。Bruno 等[12]使用Landsat 8 传感器的多时间光谱信息,提出了一种结合多分辨率分割和随机森林(Random Forests)分类的方法,实现巴西东南部农业地区的五种作物自动分类。与常规的地面遥感研究方法比较,谷歌地球引擎(Google Earth Engine,GEE)优点更明显,可以迅速获得大范围的高分辨率图像信号。目前Google Earth Engine云平台己经广泛应用于多个学科领域,如:植被制图与监测[13]、农业应用[14]、土地利用(覆盖研究)。
本研究以Google earth engine 为平台,采用Sentinel-1 雷达影像、Sentinel-2 多光谱影像为主要数据源,使用三种机器学习算法进行小麦作物的遥感分类,通过精度评价选择出最优的数据源以及机器学习算法。结合地面的实测产量来估算奇台县小麦秸秆资源的存量和分布。利用遥感技术可以更直观有效地对农作物种植面积和产量进行识别和估算,为行业管理部门监督指导农业活动提供基础数据。
奇台县地处中国新疆东北部,东西横间隔约150 km,南北纵线约250 km。由南向北地理环境特殊,地形地貌错综复杂,地形北部高,中部较低,呈马鞍状,属于中温带大陆性半荒漠干旱性气候,主要栽培农作物有小麦、玉米、西瓜、西葫芦和甜菜。通常在4至5月播种,8至9月收获[1],见表1。
表1 奇台县主要农作物物候期
样本数据采用户外实地调查数据,采集时间6月、7月与小麦关键物候期接近。野外实测数据中分别调查了小麦、玉米、西葫芦、西瓜、甜菜等11 种农作物田地、林地及裸地,样地数据共1 075 块。其中小麦的样本选取有302块,玉米的样本选取有356块,西葫芦的样本选取有190块。本次分类只选取在研究区中农作物种植面积较大的前五种作物作为样本点。同时在GEE 平台上使用人工目视的方法在线选取水体、建筑物、道路、撂荒地(戈壁)和林(草)地块[3],见表2。在样地块中随机选择70%的样品数据作为训练样本,另外30%为验证样本,用以判断分析及验证结论的可靠性。
表2 样本地块数据统计
1.2.1 Sentinel-1数据
本研究采用谷歌地球引擎(GEE)系统开展信息收集与管理。通过Sentinel-1工具箱对数据进行了预处理。预处理过程涉及辐射校正、消除热噪声以及地形修正。针对本次研究,使用了Sentinel-1中的干涉宽幅模式(IW)的地距多视图像方式作为数据源,该图像方式包括了VV 极化和VH 偏振的双极化模式。为了消除Sentinel-1 数据集中的斑点噪声,使用了改进的Lee 滤波器[15]来过滤反向散射时间序列。
1.2.2 Sentinel-2数据
使用GEE 平台获取已经经过大气校正的大气底层反射率的哨兵L2A 级数据,其中包括了所有10 m 和20 m 分辨率波段,波段为蓝、绿、红、红边1、红边2、红边3、红边4、近红外(NIR)、短波红外1(SWIR1)和短波红外2(SWIR2)。并在GEE 中构建了一个函数,使用QA60 波段检测和云掩蔽来获得云量较少的Sentinel-2 图像,见表3。
表3 遥感影像参数
为了防止其它土地覆盖类型对作物识别造成影响,本研究在进行作物分类时,先通过欧空局的ESRI十米年度全球农田覆盖产品数据及掩膜耕地信息,得到研究区的Sentinel-2农田图像。在时间序列范围内,最终共得到2022 年奇台县作物生育期内云量较低(云量小于20%)的研究区合成影像共6幅,其中四月到五月1幅、五月到六月1幅、六月到八月2幅、八月2幅。
1.2.3 数据特征选取
1.2.3.1光谱波段 选择了Sentinel-2卫星波段(共13个)中的9个,主要涉及用于进行土地覆盖变化观测、植被调查、陆表调查和环境监测的波段,即B2、B3、B4、B5、B6、B7、B8、B8A和B11波段。
1.2.3.2植被指数 根据研究区域植物的生态环境与物候特征,可以选择利用Sentinel-2影像的光谱反射率与相关植物指标来帮助划分,包含归一化植被指数(Normalized Difference Vegetation Index,NDVI)、增强植被指数(Enhanced Vegetation Index,EVI)、比值植被指数(Enhanced Vegetation Index,RVI)、归一化差异红边植被指数(Normalized difference red edge,NDRE)。其公式如下:
其中,G:增益系数,2.5;C1:Red波段气溶胶阻抗系数,6;C2:Blue波段气溶胶阻抗系数,7.5;L:冠层背景调整因子,1。
1.2.3.3 Sentinel-1 双极化特征 光学遥感方法在农作物的关键生长期内的参数信息获取方面依然存在明显的局限性,如可能受云雨天气影响,导致光学遥感产品数据源无法保障。与之相比,雷达遥感不受光照、云雨等天气条件的影响,合成孔径雷达成像的分辨率不受距离影响,通过对多次雷达回波信号进行处理,可生成高分辨率的成像结果。由于成像分辨率不随距离变化,因此可以实现在远距离下的高精度成像[16]。加入Sentinel-1 数据VV、VH 双极化特征用于对Sentinel-2 影像分类的补充。
首先通过GEE 平台对遥感影像进行筛查,以筛选出不到百分之十云量的合成影像。选择了欧空局的ESRI十米年度全球农田覆盖产品数据,掩膜出只包含农田的研究区域数据,然后将五种作物的样本点导入三种机器学习分类器,通过比较光谱、极化波段和植被特性对作物空间分布提取的影响,筛选出最佳的分类模型,用确定最优特征波段组合和机器学习算法来获取研究区小麦等作物的分类数据。
1.3.1 机器学习监督分类方法
决策树、随机森林和支持向量机是近年来遥感监督分类最普遍使用的机器学习方法。
1.3.1.1 决策树 决策树(Classification And Regression Tree,CART)的影像分类过程是通过使用专家知识库、数学计算和机器学习算法等方式收集信息并设定分类规则,通过这些规则进行的。这些规则有助于算法对影像进行分类,并自动识别出不同的影像类型[17]。每个节点都是基于样本属性的判断条件,通过这些条件进行分支和决策,最终形成叶子节点代表的类别划分。
1.3.1.2随机森林 随机森林分级计算是很有效的一种CART 决策树机器教学方式,其本质是对决策树分级算法的改进。该计算可以防止由于单棵决策树所造成的过模型拟合现象[18],旨在小样本容量情况下保证较好的稳定性[19]。通过分析发现,本研究设置准确率最高的颗数树的数量为110,其他参数设置为默认值。
1.3.1.3 支持向量机 支持向量机(Support Vector Machine,SVM)由Cortes 和Vapnik 在1995 年提出,基于统计学方法理论的一种算法模型[20]。SVM 中的软间隔概念可以应用于农作物分类中,解决线性不可分的问题,从而能够有效降低计算量。
1.3.2 精度评价方法
混淆矩阵是一种常用的分类结果评价方法,其主要作用是将分类结果与真实信息进行比较,并将分类精度展示在矩阵中。基于混淆矩阵,其中总体分类精度和Kappa 系数应用最广泛。计算公式如下:
其中,n 是待划分的样本类别数目;X是混淆矩阵中第i 行第j 列的元素;OA 表示总体精度,k代表Kappa 系数。
为评价不同机器学习算法在复杂土地结构制图中的特性,将通过GEE 对奇台县遥感图像光谱波段、极化波段以及对植被特征的九种综合方案分别作出精度评估。通过表4 可得知,在不同方案下的RF 分类器的总体精度和Kappa 系数都是最高,总体精度和Kappa 系数依次是0.98 和0.97。以CART 算法为例,仅采用极化波段的方法总体精度为0.73,明显低于使用其他波段特征的分类方案,但增加了植被特征之后将分类准确度提高至0.86。另外,在采用了光谱和植被特征分类方案后,CART 和SVM 分类器的方法总体精度也都有所提高,比仅采用极化波段的方案总体准确度也提升了0.12。表明在奇台县的小麦作物识别中,光谱特征和植被特征都对识别的精度有至关重要影响。尽管基于极化波段或植被特征方案的识别精度并不高,但它们对于进一步提高识别精度也有促进作用。
表4 基于不同分类方案的识别精度
将实地采集的地块数据与使用最优特征的三种分类器识别出的小麦结果进行直观对比,如图1所示,效果最好的是使用RF 分类器的图1(a),与图1(b)、图1(c)相比RF 分类器识别出的小麦地块较为完整、边缘清晰,与实地采集数据重合度最高、错分最少。而使用SVM 分类器的图1(c)效果最差,错分最多。使用RF分类器和极化波段、光谱波段以及植被特征组合对小麦识别的效果最好。
图1 三种不同分类器在相同特征条件下(极化波段+光谱波段+植被指数)的识别结果
混淆矩阵如图2 所示,可以看出用于精度验证的样本共1 664 个像元,其中被正确分类的有1 636个像元,总体效果较好。因此对于本研究区的小麦识别,使用RF分类器和极化波段、光谱波段以及植被特征组合的结果最好。
图2 最优方案分类后的混淆矩阵
依据小麦生长特点与五种作物NDVI 曲线(4 月初-8 月末)对比结果,小麦的关键物候期为四月末、六月初和七月末,分别对应小麦出苗、抽穗和灌浆乳熟期,也是对小麦进行遥感识别的最佳时期。因此,使用分类精度最高的RF 分类器和分类特征分别对这三个时期的耕地影像进行提取,识别小麦的空间分布,是判断小麦种植面积最佳提取时期的较好方式,见图3。
图3 植被指数时间序列曲线
为了对小麦3个生长期内识别效果进行更深入的对比分析,采用指标总体精度和Kappa 系数对此3 种分类结果进行评价。在小麦三个关键生长期中,抽穗期的识别精度最高,总体精度和Kappa系数分别为0.98 和0.97,比出苗期分别高出0.2 和0.1。灌浆乳熟期分类效果最差,总体精度仅为0.69,Kappa系数为0.61,见表5。说明在小麦关键物候期条件下,基于RF分类器和波段特征可以在小麦抽穗期有效提取土地覆盖信息,具有很好的地物分类识别效果。
表5 使用RF分类器对三个生长期小麦的识别精度
在2022年7月末对已收获的小麦地块的秸秆产量进行数据采集,在全县分布的4个地块中进行了实测后取平均值得到小麦秸秆每公顷产量(干重)2 700 kg。将研究区小麦种植地块识别后的数据导入Arcgis,计算小麦作物的种植面积为37 265 hm2,折算出奇台县小麦秸秆总产量100 615.5 t。
随着遥感影像分辨率的不断提高,对遥感影像处理的需求也在不断增加,特别是在进行大区域作物分类时,需要对大量的影像数据进行分析处理,这对本地电脑的硬件和性能提出了更高的要求。因此,越来越多的人选择在云平台上进行影像数据的处理。可以将GEE 平台中的遥感数据上传至Google 云平台,利用其数百万服务器进行并行处理,从而大幅缩短遥感影像处理所需的时间和数据处理成本。这一方式可有效提高遥感数据处理的效率,为相关领域的研究和应用提供了更为便捷的支持[21]。本研究基于GEE 云平台,利用Sentinel1、2 影像通过编程迅速地调用并管理这些数据,从而使原本本地计算平台需要几天甚至几周时间的数据处理分析工作能够在数小时内完成,大大地提高了效率、降低了成本。
在特征选择方面,本研究主要选择Sentinel-2光学影像的光谱波段和植被指数。而SAR 数据与光学数据相比能提供更丰富的农作物纹理信息,可以有效解决光学数据中“异物同谱”、“同物异谱”的现象,提高农作物识别精度[22]。本研究中可以明显发现雷达数据的导入对小麦分类精度有明显的提升作用。
此外,本研究在对农作物进行户外调查时,尽管获取的样本量较大且具备代表性,但实地数据空间分布不够均匀,对植被类型识别结果存在一定的影响。因此,通过适当增加农作物样方数量和农作物样方分布范围,会在一定程度上提升农作物样本的多样性、可靠性和连续性,同时将多源遥感数据融合,可以进一步提高农作物遥感分类的精确性和有效性。
本研究以Google earth engine 为平台,以奇台县域内主要农作物种植区域为研究对象。采用Sentinel-1 雷达影像、Sentinel-2 多光谱影像为主要数据源,评价三种机器学习分类器的作物分类精度,对比不同特征组合的效果,并使用地面实测数据来估算小麦秸秆产量。结论如下:对于小麦作物分类,使用随机森林分类器结合光谱波段、极化波段和植被特征识别精度最高,总体精度达到98%。对比小麦出苗期、抽穗期和灌浆乳熟期识别精度,其中以抽穗期的识别效果最佳。通过实地测产方法估算奇台县小麦秸秆总产量为100 615.5 t。