盛艳芳, 买买提·沙吾提, 何旭刚, 李荣鹏
(1.新疆大学地理与遥感科学学院,新疆 乌鲁木齐 830017; 2.新疆大学新疆绿洲生态重点实验室,新疆 乌鲁木齐 830017; 3.新疆大学智慧城市与环境建模自治区普通高校重点实验室,新疆 乌鲁木齐 830017)
关键字:Google earth engine (GEE); 特征优化; J-M距离; 特征集
新疆是世界六大果品生产带之一和久负盛名的“瓜果之乡”[1]。自2018年开始以“南疆特色林果提质增效,助力脱贫攻坚”举措[2]为基础,全力推进林业扶贫工作,南疆林果业规模快速扩张,林果面积每年增加百万亩,形成环塔里木盆地林果主产区等林果基地,在推进新疆脱贫攻坚方面具有重要作用。然而,南疆林果业种植面积大、地块较为破碎、果粮间作种植模式普遍存在,传统的果园种植面积监测方式已经难以满足当前林果业快速发展的实际需求。
近年来,随着遥感(多光谱、中高分辨率影像)的广泛使用,以及机器学习算法、谷歌地球引擎(Google earth engine,GEE)云计算平台不断发展,遥感技术在林果业信息提取、长势监测和估产等方面发挥了重要作用[3-6]。利用 Sentinel[7]、HJ[8]、GF[9]系列等卫星影像,采用随机森林(RF)[10]、支持向量机(SVM)[11]、BP神经网络(BP)[12]等分类技术对小麦、玉米等常见的农作物进行种植面积信息提取[13]。RF分类在调查果园分布的应用和研究中广泛使用[14-17]。另外,大部分研究主要集中在江西[18]、浙江[19]、福建[20]等丘陵地区,并且研究对象也仅局限于柑橘等芸香科果树。目前对于西北干旱区果树,尤其是新疆地区果园的遥感监测缺乏深入研究。
基于现状,还需要加强对果树遥感分类的研究,不断增加遥感影像、改进面积与分布信息的提取方法以适应技术发展要求,实现快速提取果园面积空间分布状况的目的。本研究通过中高分辨率卫星影像,充分挖掘多源卫星遥感特征量,发展基于GEE平台和机器学习的果园智能识别技术,突破高精度果园空间分布遥感制图的技术瓶颈,为果园种植面积空间分布调查和定期动态更新提供技术支持。
渭干河-库车河三角洲绿洲(渭-库绿洲)是目前新疆重要的林果产业基地,位于新疆塔里木盆地(图1),是典型的冲积扇绿洲,属暖温带大陆性干旱气候,平均气温10.0 ℃,日照2 845~2 977 h/a,光热资源充裕,年降水量59 mm。相对充足的水分、热量,形成得天独厚的光、热等自然资源,有利于果树生长。
图1 研究区示意图
遥感数据源于GEE平台,包括Sentinel-2多光谱影像和Sentinel-1雷达影像(表1)。预处理包括:轨道参数标定、热噪声去除、辐射定标、地形校正、后向散射系数提取等,通过GEE内嵌算法统一坐标系确保不同数据源之间的几何配准精度,空间分辨率统一重采样至10 m。使用平台JavaScript语言构建模型进行果园信息提取识别。
表1 本研究使用的遥感数据
野外数据:分别在2021年5月26日至6月3日、2022年6月28日至7月7日进行 2次野外实地考察,结合研究区域地物类型以及相关国家标准(GB/T 21010-2017)确定研究区地物为:果园、建设用地、林地、草地、水域、耕地、其他地物,共获得典型地物样点920个。此外,使用Google earth影像结合野外获得的典型地物样本点进行目视解译来扩充样本数量,使样本点均匀分布在整个绿洲内部,最终样本总数为1 139个,训练集和验证集比例为7∶3。
本研究基于GEE平台,使用Sentinel-1和Sentinel-2数据,计算得到6类特征,使用3种特征优化方法结合随机森林分类(Random forest,RF),确定提取果园信息的最优特征集(图2)。
图2 技术路线图
1.3.1 分类特征参数和特征优化方法 通过计算得到光谱特征、植被指数特征、红边特征、纹理特征、极化特征和物候特征构建分类特征集[21],共包含42个子特征(表2)。
表2 分类特征参数
1.3.2 特征集优化
1.3.2.1 基于J-M距离的特征集优化 J-M距离衡量不同地物类型在波段间分离能力,对各地物类型的可分性进行定量分析,并据此进行特征集优化[22]。表达式为:
JM=2(1-e-B)
式中B表示某一特征的巴氏距离(Bhattacharyya distance)。不同种类样本的巴氏距离同时样本满足正态分布的时为
式中ek表示某类特征的均值,δk2表示某类特征的方差。
1.3.2.2 基于随机森林属性重要度的特征集优化 随机森林分类器(RF)是一种非参数分类器,由一组决策树组成,具有便捷的特征选择算法[23]。利用RF中的属性重要度,进行有效数据降维,根据得分由高至低对特征集进行优化。
1.3.2.3 分类方案制定和精度验证 本研究利用GEE平台中最优棵树选择,最终确定待生成的决策树(Ntree)数目值为100。为了说明有利于果园信息提取的最佳特征,设计了23种试验方案(表3)。
表3 分类方案
1.3.3 精度验证 研究结果表明,总体精度(OA)与Kappa系数在遥感分类精度验证中是最基本、最能衡量分类结果的指标[24]。其中OA是指被正确分类的像元之和占所有地类像元总数的百分比,Kappa系数表明分类结果与验证样本数据之间的一致性。因此本研究使用OA、Kappa系数作为总体分类精度评价体系,使用统计年鉴数据对提取的果园种植面积信息进行验证。
为了确定果园信息提取的最佳时间窗口,使用Sentinel-2数据,综合研究区域典型地物样本点建立归一化植被指数(NDVI)时间序列曲线,以及使用Timesat3.3软件提取植被物候特征。
NDVI时间序列曲线,反映了植被与非植被在2020年NDVI变化特征[25]。如图3所示,非植被(建设用地、裸地、水域)的NDVI值全年较低。3月植被(果园、林地、耕地、草地)NDVI最小(NDVI<0.1);3-6月植被NDVI开始快速上升,其中果园上升速度最快,NDVI最高,尤其在5月远高于其他地物(NDVI>0.35);6-8月,植被NDVI逐渐增长至峰值,即生长旺盛期;8月裸地、建设用地、草地、果园、林地、耕地NDVI均开始下降,果园与林地NDVI相近。
图3 2020年归一化植被指数(NDVI)时间序列曲线
植被物候信息如图4A所示,5月果园NDVI高于所有其他植被,结合图4B可知果园生长季在5月(a点)开始,此时果园在所有植被中NDVI最大(NDVI>0.4);而后所有植被的NDVI继续增加至生长季中期(6-7月),7月中旬所有植被的NDVI达到最大值,即各植被均到达生长旺盛期(NDVI>0.4);8月开始所有植被的NDVI逐渐下降,10月生长季结束,NDVI迅速下降。
A图为4种主要植被物候曲线图;B图为果园物候信息点图。物候特征a~k见表2。
为了探讨雷达数据在果园信息提取中的作用,基于2020年全年12期的Sentinel-1雷达遥感影像,由VV、VH极化数据算出典型地物的后向散射系数,得到典型地物后向散射特征时间序列(图5)。类型不同的地物,散射特征随着生长发育期的变化表现出不同的散射机制[26]。
图5 2020年VV与VH极化后向散射系数时间序列曲线
从时间尺度,以及不同极化方式上来看,VV后向散射系数变化特征相比于VH显著。果树的萌芽期和展叶期在4-5月,在VV极化中果园与耕地区分不明显,趋势几乎一致,在VH极化中5月果园略高于耕地。本研究果园中的果树与林地中的树木存在一定的相似性,极化数据在不同类型作物的识别中具有较好的识别效果[27];林地在VV极化方式下,1-3月与果园的后向散射系数趋势一致,在4月开始有所区分,林地后向散射系数高于果园,峰值在6月达到-6.99 dB,在VH极化中,至8月前都有较好的分离性,因此提取果园信息的最佳时间点在5-7月。
最后,综合Sentinel-1极化影像后向散射特征和Sentinel-2多光谱影像计算得到的NDVI时间序列与物候特征,可以确定5月为果园信息提取的最佳时间窗口。
2.2.1 分类方案特征集优选 所有特征排列组合形成不同的分类方案,并通过筛选得到17个特征方案(表3),使用RF分类(图6)。在原始特征组合而成的分类方案中,所有方案总体精度均达到了80%以上,总体精度排列前三的方案分别是G12、G4、G17。
G1、G2、G3、G4、G5、G6、G7、G8、G9、G10、G11、G12、G13、G14、G15、G16、G17见表3。
2.2.2 基于J-M距离与属性重要度的特征集优化 获取J-M距离和属性重要度对各特征的评分,依据2种评价方法,依次加入特征,观察果园的分类总体精度变化(图7)。2种方法表现出一些类似的特点,随着所用特征数逐渐增加,总体精度表现为快速增长,当特征数达到10即B6处时,总体精度开始趋于稳定,特征数达到40即B8_contrast时能达到最优精度约89%。
根据图7,按照6个特征类别进行分类统计(表4)。通过对比,在两种优化的方法中,纹理特征的可分离性,对分类贡献度均较高,纹理作为物体表面的一种基本属性广泛存在自然界中,是描述和识别物体的一种极为重要的特征[28],果园因其均匀、粗糙等纹理特征明显不同于其他地类,尤其是与林地间的差别较大。
A:J-M距离优化方法;B:属性重要度优化方法。a1:B8_asm;a2:B8_ent;a3:NDVIre3;a4:B8_idm;a5:a;a6:B8_sent;a7:e;a8:EVI;a9:d;a10:B6;a11:b;a12:B11;a13:f;a14:NDWI;a15:VV;a16:B8a;a17:B7;a18:B8;a19:NDVIre1;a20:B12;a21:VH;a22:NDVI;a23:NDre2;a24:NDre1;a25:B5;a26:B3;a27:B4;a28:B1;a29:B9;a30:NDVIre2;a31:j;a32:g;a33:B2;a34:i;a35:k;a36:h;a37:c;a38:B8_var;a39:B8_diss;a40:B8_contrast;a41:B10;a42:B8_corr。a、b、c、d、e、f、g、h、i、j、k、NDVIre1、NDVIre2、NDVIre3、NDre1、NDre2、B1、B2、B3、B4、B5、B6、B7、B8、B8A、B9、B10、B11、B12、B8_asm、B8_contrast、B8_corr、B8_var、B8_idm、B8_diss、B8_sent、B8_ent、NDVI、EVI、NDWI、VV、VH见表2。
表4 特征优化表
2种方法中排名靠前的物候特征是a、b、e,即生长季开始、生长季结束、生长季中期,它们可以直接反映植被的生长信息。所以纹理特征与物候特征相比其他特征更能凸显不同植被间的差异,更好突出果园与其他植被之间的差异[29]。
按照优化后的特征集,对原始特征组合中精度较高的方案G4、G12、G17依次优化。如表5所示,其中精度最高的方案G17JM,总体精度比优化前的方案高2.05个百分点,其次方案G12JM的总体精度比G12高1.16个百分点;方案G4、G4JM、G4VIP总体精度均未达到90%,因此参与的特征数较少时,识别效果较差。
表5 分类精度对比
总体而言,使用J-M距离优化的方案总体精度均较高,相比于原始特征组合方案有1~2个百分点的提升,比使用属性重要度优化的高1~3个百分点,优化效果较好。
精度最高的3种方案如图8所示,从空间分布上来看,果园主要分布在库车市和新和县;沿河流与农田交错分布,呈片状、由上至下扇形分布在绿洲中上部。
A、B、C分别为G12、G17JM、G17VIP方案提取的果园信息空间分布。G12、G17JM、G17VIP见表3。
为了比较分类效果,选取经实地考察并包括果园周围常见地类(农田、林地、建设用地等)的局部区域与高分影像(1 m)对比分析(图9)。G12方案与G17VIP方案识别效果相似,可以完整识别果园地块,但部分其他用地被错分为建设用地,林地识别不完整,草地提取范围过大;G17JM方案识别效果更精细,果园地块周围的耕地、草地识别完整,可以将果园与林地准确分开。分类识别的果园种植面积为66 921.62 hm2,统计年鉴[30]中库车、沙雅、新和果园总面积为81 066.67 hm2,面积精度为82.55%。故J-M距离优化特征的方案G17JM总体精度最高,识别效果最好。
A为高分原始影像;B、C、D分别为G12、G17JM、G17VIP方案提取的果园信息空间分布细节图。G12、G17JM、G17VIP见表3。
利用GEE平台强大的数据运算和存储能力,使用Sentinel-1/2遥感数据和RF分类算法,在果园信息提取中对特征进行优化的研究相对较少,大多研究集中在农作物监测方面[31],随着GEE平台的发展,逐渐出现不同的研究方向如湖泊面积测量[32]、沙漠化监测[33]、黄河大尺度面积制图[34]等,均取得较好的效果。渭库绿洲作为典型干旱区绿洲,是新疆重要的林果产业基地。种植的人工林区域(果园)与自然林地存在区别,果园地块破碎但边界整齐,纹理规律,经实地调查,以及综合采集多种果园地物,可以发现不同类型果树物候期存在细小差别,茂盛期基本一致,可以利用茂盛期相似的纹理特征、极化特征等,对果园信息进行提取。赵安周等[35]、徐超等[36]、张祯祺等[37]、宁晓刚等[38]从中尺度上进行的动态监测研究结果,以及本试验为实现大面积果园动态监测所做的研究结果,都证明了果园信息泛化、统一大面积提取的可行性。但是根据以往针对果园的研究中[14-16],大多使用传统的监督分类方法进行面积信息提取,在今后的研究中可以考虑更多方法的对比研究,如面向对象结合决策树、面向对象结合随机森林等方法;此外,Sentinel-2 数据空间分辨率为10 m,针对研究区域内相对较小的破碎果园地块无法识别,可以尝试精度更高的遥感影像,如GF_2等,得到更高的制图精度。
从研究结果看,对比3种优化方法,J-M距离优化方法可以有效提高识别效果,这与宁晓刚等[38]的研究结果一致,充分说明使用J-M距离是进行特征优化时优于属性重要度优化方法。此外,通过总体精度变化图发现,适用于果园信息提取的最优特征集应最少包含10个子特征,当特征数达到10时,总体精度开始趋于稳定,分类可以取得较好的效果。
研究结果表明,使用RF识别果园时物候、纹理、极化等特征结合J-M距离进行特征集优化可以降低数据量、提高计算效率,使用得到的最优特征集进行分类,总体精度与面积精度均高于80%,效果较好。本研究基本都在GEE平台中进行并获得了较好的分类总体精度,证明了使用GEE实时获取渭库绿洲果园面积动态变化的可行性,为果园面积动态监测提供强有力的基础。
本研究结论如下:(1)本研究共计包括23种分类方案,对比方案精度,所有精度均到达80%以上,其中精度最高的G17JM方案,其总体精度为91.25%,Kappa系数为0.89,面积精度为82.55%;(2)综合使用NDVI时间序列曲线、物候曲线和极化时序曲线确定的窗口为5月,即5月为果园信息提取的最佳时间;(3)综合对比3种优化方法,J-M距离为最佳优化方法,得到的优化方案(G17JM)总体精度较高,且效果最好。提取果园信息的最优特征集,具体为:B8_asm、B8_ent、B8_idm、NDVIre3、B6、B7、a、e、b、EVI、B11、B8A、B8、VV。