基于Sentinel-2遥感影像的农作物分类与适宜性评价

2023-11-16赵孟辰阿里木江卡斯木

西南大学学报（自然科学版） 2023年11期

赵孟辰，阿里木江·卡斯木

新疆师范大学地理科学与旅游学院，乌鲁木齐 830054

农业产业结构调整优化是当前我国建设现代化农业,保障农民持续增收,实现农业可持续发展的必然选择[1]．农作物种植结构是指一个地区主要农作物类型和其空间分布状况,及时有效地获取农业种植结构信息对农作物产量监测、农业用水精细化管理、区域农业、可持续化发展等具有重要意义．以往的研究依据统计部门的逐级上报和抽样估测数据集获取农作物种植结构信息只能停留在量的层面,难以实现农作物空间分布可视化,且此方法存在费时、费力、时间分辨率低和主观性强等缺点．遥感技术具有大范围数字成图、高时效性和成本低等优势,目前被广泛用于资源普查,并逐步成为获取农业种植结构信息的重要方式之一．

国内外众多学者对基于遥感影像的农作物识别、分类及种植区域适宜性进行了大量研究．Johnson[2]采用多尺度分割算法进行分割,通过SVM(Support Vector Machine,SVM)分类模型,计算不同尺度下影像对象隶属度及隶属度类别;Shackelford等[3]以IKONOS 影像为实验数据,采用面向对象和基于像素的方法分别进行分类,结果证明面向对象的方法分类精度更高;Strom[4]分别利用面向对象和传统目视解译的方法对SPOT(Systeme Probatoire d’Observation de la Terre)影像进行分类,结果证明面向对象分类方法的精度高于传统目视解译方法;Zhang等[5]基于最近邻分类算法、最大似然分类算法及SVM算法,提出一种面向像斑的空间特征分类方法．陈杰等[6]在SVM算法的基础上,结合粗糙集粒度的计算方法对影像进行分类;顾海燕等[7]验证了随机森林算法在运行速度和分类准确性上有明显优势．

不同特色经济作物的生长发育与气候、地形、土壤等自然环境条件密切相关,在不同自然环境条件下的特色经济作物产量及品质有明显差异,因此如何根据各地区的自然条件,研究区域特色经济作物用地布局,因地制宜地发展特色经济作物,实现区域土地资源合理利用具有重要意义．传统的土地适宜性评价方法计算和处理效率低而且精度差,利用GIS(Geographic Information Systems,地理信息系统)技术进行土地适宜性评价,不仅可以将评价单元的空间数据与评价指标数据结合起来,还可对这些数据进行分析和处理,从而提高土地适宜性评价的效率和精度．龚珍等[8]、吴波等[9]分别使用GIS技术在我国不同地区针对茶树、苹果树等多种不同品种的植物进行了农业气候适宜性区划研究;贾艳红等[10]以立地条件和土壤性状中海拔、土壤质地等6个因子为指标,建立了罗汉果种植的土地适宜性评价体系,借助层次分析法确定评价因子权重等级标准,划分土地适宜性等级,并利用GIS技术实现了桂林市罗汉果种植土地适宜性综合评价．

本文基于Sentinel-2影像数据,以河北省张家口市沽源县为研究区,选择胡萝卜、西蓝花、玉米、马铃薯、油麦菜、生菜、白菜、野燕麦这8类研究区内主要种植的农作物,分别采用随机森林分类、支持向量机(SVM)分类和决策树分类方法,对研究区内主要农作物的分布情况进行遥感分类,并将不同方法获得的研究区内农作物分类精度进行对比．通过遥感技术分析研究区内主要农作物的种植情况,从土地和气候两个方面分析河北省张家口市沽源县主要农作物的生长适宜性,为政府相关部门了解农作物种植情况提供辅助信息,同时为进一步科学规划农作物种植结构提供参考．

1 研究区概况与数据预处理

1.1 研究区概况

河北省张家口市沽源县位于河北省西北部坝上地区,与内蒙古自治区接壤．全县地势南高北低、东高西低,平均海拔1 536 m,其中北部、西部为波状高原,东部为低山丘陵,南部为中低山地,全县总土地面积为3 654 km2,耕地面积15.5万hm2．沽源县“冬无严寒、夏无酷热”,为冷温带半干旱大陆性季风气候[11],年平均气温1.4 ℃,≥10 ℃积温1 900 ℃左右,年平均降水量409.9 mm,年蒸发量1 787.5 mm[12]．土壤类型有灰色森林土、栗钙土、草甸土、沼泽土、盐土,其中栗钙土为主要类型,占全县土壤分布面积的88%．植被类型为干草原稀树草原,主要由耐寒旱生多年生草本植物组成．森林覆盖率为35.64%,有杨、榆、桦、落叶松和云杉,部分地带有杞柳、水策棠梨、枸杞等;粮食作物有小麦、莜麦、马铃薯等;经济作物有胡麻;蓄牧业主要发展大牲畜,马、牛的头数逐年增加,沽源牧场是河北省3大牧场之一．

1.2 数据来源与预处理

1.2.1 Sentinel-2 影像

哨兵2号(Sentinel-2)数据是欧洲“哥白尼计划”的重要组成部分．Sentinel-2卫星搭载的MSI成像仪能够覆盖13个光谱波段,幅宽为290 km,空间分辨率分别为10 m(4个波段)、20 m(6个波段)、60 m(3个波段)．在红边范围内有3个波段,可以有效监测植被健康信息,也能够用于森林监测、气候变化、土地利用变化、产量预测、山体滑坡等应用,极大地提高了陆地监测水平．在光学数据中,Sentinel-2数据是唯一一个在红边范围含有3个波段的数据,这对监测植被健康信息非常有效．本文使用的Sentinel-2多光谱数据下载自欧空局(European Space Agency,ESA)数据中心(https：//scihub.copernicus.eu/dhus/),为表观反射率数据,对该数据的预处理过程如下：

1) 辐射定标和大气校正：使用SNAP(Sentinels Application Platform)软件中的Sen2Cor插件对Sentinel-2数据完成辐射定标和大气校正,得到地表反射率数据,采用双线性内插法将所有波段的分辨率重采样为10 m,最后将校正结果转换为ENVI格式．

2) 波段合成经SNAP软件输出的ENVI(The Environment for Visualizing Images)格式Sentinel-2影像数据存储结果为一个波段对应一个图像,本文采用波段合成工具将其顺序调整后合成多波段图像,同时将坐标系设置为UTM(Universal Transverse Mercator Projection-通用横轴墨卡托投影)投影,基准面为WGS(World Geodetic System,世界大地坐标系)-84．

3) 图像镶嵌和裁剪：图像裁剪实验区所使用的Sentinel-2影像数据单景即可覆盖,本文通过河北省张家口市沽源县矢量数据在GIS软件中对所使用的Sentinel-2数据进行批量裁剪,由此得到实验区影像．

1.2.2 野外调查数据

为保证利用遥感影像进行数据标注的准确性,本文对研究区展开野外实地数据调查,通过实地走访和拍照记录等方式进行采集．根据研究区规模,为了在调查过程中使用训练样本进行更高精度、更符合实际情况的分类,记录了实地调查样本点的经、纬度,对实测样点的农作物进行拍照留档,以期更加准确地把握各类主要农作物的特征．采样过程利用手持GPS(Global Positioning System,全球定位系统)设备记录采样点位置信息,同时将采样点在奥维地图卫星影像数据中打点标记,以采样点编号进行备注．将样本数据按照1.5∶1的比例随机分为训练样本数据集和验证样本数据集(图1)．

图1 河北省张家口市沽源县作物样本点分布图

2 研究方法

2.1 遥感影像分类算法

2.1.1 支持向量机法

支持向量机在线性可分的多特征空间中计算出最优超平面,该平面在将被分类样本正确分开的同时使分类间隔最大化[13]．通过SVM中软间隔概念,将需要分类的样本通过引入核函数映射到高维空间中,同时将解平面由线性扩展为非线性并求得最优分类超平面[14-15],能够有效地解决农作物分类中的线性不可分问题,减少运算量,使SVM更加实用化．SVM算法通过选择恰当的函数子集、合适的对应子集的判别函数,依据有限样本构建最佳模型,使机器学习的实际风险降到最低．SVM分类器采用二次规划寻找最优解,利用核函数,即使面对高维特征也可获得较好的结果,有利于解决非线性分类问题．

SVM可以表示为二次优化问题,公式为：

(1)

(2)

(3)

在SVM分类器中,核函数的选择十分关键,常用的核函数包括4种[16],其中高斯径向基核函数能够有效解决线性不可分问题,应用较为普遍[17-18],因此本文选用高斯径向基核函数进行农作物分类．

2.1.2 决策树法

决策树法适合对无序、无规则的多种空间数据通过寻找内部暗藏规律进行分类．在决策树模型使用中需要考虑算法特性．决策树模型的优点：① 浅层的(Shallow)决策树视觉上非常直观,而且容易解释;② 对数据结构和分布不需作任何假设;③ 可以捕捉住变量间的相互作用(Interaction)．但是,决策树法也存在深层视觉上和解释上比较困难、容易过分微调样本数据而失去稳定性和抗震荡性、对样本量(Sample Size)需求比较大和处理缺失值功能非常有限等缺陷．

2.1.3 随机森林法

随机森林(Random Forest,RF)是一种特殊的套袋算法,本质是决策树法的改进．该算法特性良好,相对于单棵决策树而言泛化能力更强、稳定性更高,运算速度较快且鲁棒性较好,能够避免因单棵决策树引起的过拟合现象,在小样本容量情况下依旧能够保持较好的稳定性,具备良好的抗噪能力,容易实现,分类效果较好,在农作物种植信息提取应用中优势明显．在随机森林分类算法模型构建过程中,特征变量数量m和决策树数量n为需要设置的参数,应用不同其取值不同,计算公式为：

Model=RF_train(Index,Value)

(4)

Result=RF_class(Index,Value)

(5)

式(4)-式(5)中,Model为RF建立的模型,Index为RF输入变量的数据集,Value为分类结果隶属的代表数值,Result是RF的分类结果．

2.2 层次分析法

层次分析法是一种多层次权重分析决策方法[19],被广泛应用于农作物适宜性评价．层次分析方法是将复杂问题先看作一个系统,然后对其各组成要素进行两两分析,合并层次形成各组成要素间有序的层次,进而形成有序的层次模型[20]．然后,根据层次间的重要性对每一个层次因素给出判断矩阵,求出每一层次全部要素的相对权值,就会得到最终每个因子的权重值,否则继续进行层次划分[21]．层次分析是一种决策分析方法,也是一个量化过程,根据各因素之间的比较和计算可以得出相应的决策方案,这些方案可以被选为最佳方案或作为被选为最佳方案的依据．

3 基于Sentinel-2数据的农作物分类

3.1 基于Sentinel-2影像植被指数数据提取

利用多时相的Sentinel-2影像提取光谱特征、分割特征、颜色特征及纹理特征,构建特征数据集,这些特征可以将土地覆盖类型的光谱与空间信息进行有效地结合,从而极大地提高遥感影像中农作物的识别能力和准确性,改善作物制图的准确率．

植被指数由卫星影像数据不同波段组合而得,将光谱信息压缩为一个通道,能够反映地表植被的生长状况[22-23]．多种植被指数,如归一化建筑指数、归一化水体指数、归一化差值红边指数、简单比值指数和红边叶绿素指数已经被广泛应用于植被分类[24-25]、干旱监测[26]、全球与区域土地覆盖[27-28]、农作物识别[29-31]、作物和牧草估产等领域的研究中．红边是描述植物色素、健康状况的重要指示波段,能够有效地应用于参数计算、植被状态监测等领域．

已有研究发现,红边信息能够用来有效地区分地物类别．目前,几十种植被指数已被定义,对不同地物的敏感性也具有差异．其中,归一化差分植被指数(Normalized Difference Vegetation Index,NDVI)在农作物提取和监测中应用较为广泛．考虑到Sentinel-2影像的波段特点,本文除NDVI植被指数外,选取了归一化建筑指数(Normalized Difference Building Index,NDBI)、归一化水体指数(Normalized Difference Water Index,NDWI)、归一化差值红边指数(Normalized Difference Red Edge Index,NDRE1)[34]、简单比值指数(Simple Ratio Index,SRre)[35]和红边叶绿素指数(Red-Edge Chlorophyll Index,CIred-edge)共同作为指数特征进行研究,各植被指数的计算公式见表1．

表1 植被指数及计算公式

3.2 主要农作物分类精度评价

3.2.1 精度评价标准

通过计算混淆矩阵来获取各分类方案在不同分类器下农作物分类结果的精度．混淆矩阵是衡量不同分类器分类精度最直观、最便捷的方式．本文主要利用总体精度和Kappa系数这2个评价指标来分析各方案及算法的农作物分类精度．

进行常规健康教育，健康教育6 个月后，对病人集体进行1次 60 min的骨质疏松预防知识讲解，同时为每例病人免费提供1本 2 型糖尿病病人骨质疏松预防小册。

总体精度(Overall Accuracy,OA)为被正确分类的样本个数占样本总数的比例,其计算公式为：

(6)

式(6)中,OA为总体精度;N为样本总数;r为分类类别总个数;i为第j类地物被正确分类的样本个数．

Kappa系数计算公式为：

(7)

式(7)中,N为样本点像元总数;xi+和x+i分别为第i列、第i行的和;xii为主对角线上的值．

3.2.2 精度评价

本研究所选用的分割方法为棋盘分割(Chess Board),以3个像素格为单位作为一个棋盘格,一个像素格面积为900 m2．为了比较农作物分类过程中支持向量机(SVM)、决策树、随机森林(RF)3种分类器的分类精度,优选出适合地形复杂地区农作物的分类算法,通过野外实地调查获取验证数据,采用建立混淆矩阵的方法对分类结果进行精度验证．

具体评价指标为总体分类精度(OA)和Kappa系数,其中总体精度是正确分类像元数占参与分类总像元数比例的表征;Kappa系数则用来度量分类结果与实际数据的一致性,分类结果评价精度见表2．

表2 分类结果评价精度

3.3 基于Sentinel-2影像农作物分类结果分析

根据支持向量机(SVM)分类、决策树分类、随机森林(RF)分类算法的分类结果(表3),可以看出所有分类方法中只有决策树分类法未识别出西蓝花,其余方法中SVM分类法提取的作物面积最小．

表3 不同分类结果的各农作物种植面积

3.4 干预后随机森林分类结果

将随机森林分类法得到的结果进行人为干预后得到研究区主要农作物分布情况(图2)．从图2中可以看出,8种农作物广泛分布在研究区内．经分析发现,研究区内主要农作物分布整体呈镶嵌结构,结合表3中各农作物的种植面积分析可知,在干预后随机森林分类法得到的研究区主要农作物分类中,白菜的种植面积最大,约为471.11 km2,广泛分布于研究区全区;油麦菜种植面积其次,约为451.12 km2,接着从大到小依次为玉米、野燕麦、生菜、马铃薯、胡萝卜、西蓝花;其中种植面积最小的是西蓝花,约为47.58 km2(图2)．

图2 干预后随机森林分类主要农作物提取分布情况

4 种植区域土地适宜性评价

4.1 种植区域土地适宜性指标体系结构及权重确定

作物经济性评价理论认为,作物在其区域能否正常生长及取得多大的经济效益取决于3个因素：① 作物生长所处的自然物理环境,包括气候条件、地形地貌、土壤质量等．② 作物生长的社会条件,包括基础设施情况、人文环境等．③ 作物本身对经济适宜性评价的影响,包括叶绿素含量、水分含量、有效氮含量等．这3个因素决定了物种的地理空间分布,本文也正是基于这个理论来进行作物种植经济适宜性评价研究．

河北省张家口市沽源县2012年土地变更调查耕地面积为1 319 km2,主要农作物有燕麦、马铃薯、亚麻、杂豆、蔬菜等,是华北最大的脱毒薯基地．根据乡镇种植意向调查资料显示,蔬菜、莜麦、甜菜、玉米、杂豆等农作物种植面积增长,零星种植的马铃薯、小麦、油料等种植面积有所下降．研究中考虑到评价因子的相关性、稳定性及各因子数据的可获得性,选取海拔、坡度、坡向、土壤类型及土壤侵蚀度作为研究区主要农作物种植土地适宜性评价因子．采用层次分析法(Analytic Hierarchy Process,简称AHP)来确定各评价因子的权重(表4)．

表4 种植区域土地适宜性评价指标及权重

4.2 种植区域土地适宜性分级结果分析

根据河北省张家口市沽源县主要农作物种植生态适宜性指标体系,将土壤侵蚀度和土壤类型作为土壤要素因子;并将海拔、坡度和坡向3个地形指标综合为地形要素因子．参照以往的研究,根据层次分析法(AHP)的具体步骤,利用和积法对判断矩阵进行归一化和正规化处理,确定各评价因子权重[32],并对各因子的影响程度进行分级．按照自然分断法,划分为3类土地适宜性等级[33]．根据《土地评价纲要》中对特定用途土地适宜性程度评价等级划分要求及河北省张家口市沽源县主要农作物种植的自然条件,将土地适宜性分为适宜、中度适宜和不适宜3个等级(表5)．

表5 主要农作物种植区域土地适宜性分级指标

4.3 种植区域土地适宜性评价

利用GIS层次分析法与土地适宜性分级指标进行评估,获得河北省张家口市沽源县主要农作物种植土地综合适宜性等级分布图(图3)．由图3分析可知,研究区内中度适宜种植当地主要农作物的用地面积最大,其次是适宜种植面积,不适宜种植地区的面积最小．从整体上看,河北省张家口市沽源县适宜和中度适宜种植用地均匀分布在研究区内,不太适宜种植用地主要分布在沽源县白土窑乡、西辛营乡、西坡、闪电河乡等地区．由此分析,对于沽源县主要农作物种植中度适宜地区,建议培肥地力,提升该区域的土地质量;在不适宜地区发展其他产业,更加合理地规划利用土地;在种植适宜区应完善现有种植园,创新种植技术和生产模式,促进当地产业发展(图3)．

图3 区域种植土地综合适宜性等级分布图

5 结论

5.1 研究区内农作物种植信息

以Sentinel-2提取的NDV,NDBI,NDWI,NDRE1,SRre,CIred-edge数据及组合为分类特征,分别采用SVM分类法、决策树分类法、随机森林分类法提取研究区内农作物种植信息．结果表明,随机森林分类法对研究区内8种主要农作物进行分类的精度最高,其总体准确率为65.10%,Kappa系数为0.587 1．

5.2 研究区内主要农作物分布情况

将随机森林分类法得到的结果进行人为干预后得到研究区主要农作物分布情况,发现8种农作物广泛分布在研究区内,整体呈镶嵌结构,其中马铃薯种植面积最大,约为22 728 hm2,广泛分布于研究区全区;生菜种植面积其次,约为21 565 hm2,接着从大到小依次为野燕麦、玉米、油麦菜、白菜、西蓝花,其中野燕麦和玉米主要分布在沽源县的中部和西部地区,有较少部分分布于以“塞北管理区—丰源店”为界的东部地区;种植面积最小的是胡萝卜,约6 826 hm2,最集中区域在沽源县南部,在东北边远地区种植较少．