姚金玺,王 浪,李建忠,张 焜,张 志※
(1. 中国地质大学(武汉)地球物理与空间信息学院,武汉 430074;2. 青海省青藏高原北部地质过程与矿产资源重点实验室,西宁 810300;3. 中国地质调查局乌鲁木齐自然资源综合调查中心,乌鲁木齐 830057;4. 中国地质调查局应用地质研究中心,成都 610036)
中国农业区的土地利用类型分布分散、农业景观破碎和作物种植结构复杂等特点,给土地利用类型的遥感解译带来了持续性挑战。青海省诺木洪河一带是全国第二大枸杞种植区,作物种植面积快速调查对调控农产品市场、辅助决策和保障农业可持续发展具有重要意义。
种植区地物精准、快速分类使用传统方法费时费力,而遥感影像具有快速性和周期性的特点,利用遥感技术进行土地利用分类是一种重要的方法,且前人已经开展了大量研究。但单一传感器难以全面、准确地反映复杂的地表覆盖类型。众多地物分类是利用遥感影像不同地类的光学波段反射率或植被指数在空间域或时间域上的变化特征,其原理便决定至少需要研究区一景的清晰遥感影像。光学遥感数据受光谱分辨率、空间分辨率、时间分辨率及光谱波段的不同而产生“同物异谱”或“同谱异物”现象,会极大影响地物分类精度。Wu等利用植被指数MNDWI、NDBBI去除农田中的水体、建设用地和裸地后,针对农作物不同的生命周期使用NDVI差异性研究农田利用变化,以此为废弃农田的开采提供参考。而上述是针对不同的阶段使用光谱指数差异进行研究,若只是一景影像效果则会不尽人意。纹理特征不仅反映图像的灰度统计信息,而且反映图像的空间分布信息和结构信息。Dihkan等构建光谱特征与纹理特征进行地物分类,其分类结果进一步提高,得到了高精度的地物类型制图。Zhang等针对光谱和纹理特征空间进行进一步的分析,利用均质性指数和分裂性指数分析确定了最佳分类特征和分类器,为单一或复杂种植结构区域的农作物制图参考。然而,光谱和纹理数据均会受到云雨等因素的影响,进而导致影像质量和地物分类精度降低。
合成孔径雷达(Synthetic Aperture Radar, SAR)具有全天时、全天候的监测优点而逐渐在地物分类中占据重要作用。但SAR数据的信号会受到相干斑噪声的干扰,影响对地面目标的精确识别,Xu等使用多时相或多极化的SAR数据可以减少噪声的干扰,利用微波数据结合光学数据的方式可以进一步提高地物的分类精度。Li等使用涵盖2017年旱季的Sentinel-1A(S1A)数据的时间序列,在马朗和Lampung两个地区进行了蔬菜农作物分类,以此证明了S1A双极化SAR时间序列数据在印度尼西亚蔬菜农作物分类中的应用潜力。Swadhina等整合Sentinel-1合成孔径雷达(SAR)数据和Sentinel-2光学数据,构建光谱指数与雷达纹理特征空间,并以10 m空间分辨率绘制了农田专题图。
前人研究的对象多为光谱、纹理和极化特征的单独分析与部分特征的组合,利用多源数据、多特征共同组合分类的研究较少。故本文针对青海省诺木洪地区,基于GEE云平台,利用Sentinel-1雷达数据和Sentinel-2多光谱数据,构建多特征空间并利用特征优选算法、随机森林算法对研究区地表覆盖物进行分类,综合分析评价不同特征的重要程度、相关度以及不同特征组合对分类结果的影响。此外,本文进一步综合多源数据(Sentinel-1、Sentinel-2和GF-2)协同分类以分析不同数据源的纹理特征对分类结果产生的影响。
研究区位于青海省海西蒙古族藏族自治州都兰县诺木洪乡(图1),平均海拔2 775 m,属高原干旱大陆性气候,区内年均温1.2~4.3 ℃,降水量17.8~177.5 mm,集中在6-9月,表现出雨热同季。
图1 研究区示意图Fig.1 Schematic diagram of the study area
研究区属柴达木盆地南缘绿洲农业区,其光照资源丰富、光温生产潜力大,光热条件有利于农作物生长,主要种植春小麦、枸杞等农作物。研究区以荒漠、湿地与农业用地为主,自然植被类型有梭梭、柽柳、白刺和芦苇等。在山麓洪积扇、冲积扇区及河流沿岸等含水率高的地区,植被长势良好。
Google Earth Engine(GEE)由谷歌、卡内基梅隆大学和美国地质调查局联合开发,其中的公共数据包括了近40年的全球卫星影像以及专题图,每日更新影像量高达4000景;GEE平台还提供超过800 种功能函数,从简单的数学函数到机器学习、图像处理等操作。近年来,GEE平台逐渐被专家学者用于地表覆盖物分类中,这极大地改变了传统的数据下载、存储与处理方式,有效提高了研究效率与精度。
此次数据源包括Sentinel-1、Sentinel-2和GF-2数据,Sentinel-1雷达数据在GEE平台上经过热噪声去除、辐射校准和地形校正预处理。为了确保数据可靠性及降低噪声影响,采用改进的LEE滤波处理以去除噪声,针对数据取均值操作以获取Sentinel-1在2020年1月至12月的VV、VH波段后向散射系数。Sentinel-2数据所具有的红边波段对地表植被覆盖信息提取具有较大的应用潜力。在GEE平台对Sentinel-2数据进行了大气校正、几何校正以及利用“QA60”波段进行去云操作以减轻云污染的预处理操作,选取植被生长较好的2020年6-9月时间段内符合需求的12景影像,以中值合成方法重构最小云量合成影像。
国产高分二号卫星(GF-2)包括一个全色波段和四个多光谱波段,分辨率优于1 m,同时还具有高辐射精度、高定位精度和快速姿态机动能力等特点。本次选取了两景GF-2影像,对其进行辐射定标、大气校正、正射校正、多光谱与使用Gram-Schmidt(G-S)方法进行全色融合操作、裁剪以及影像配准等预处理。时间是2020年7月25日,一方面此时遥感影像云量比较少,另一方面是仅此时间段植被发育较好,纹理信息丰富,便于多源数据协同分类研究。
在2020年8月8日至10日、2021年7月17至24日两次开展野外调查工作,对研究区内的植被类型及分布进行调查。根据繁茂程度、长势、实地了解的大致空间分布以及在不同遥感合成图像上的表现将枸杞地分别定义为早、中、晚期。晚期枸杞地因枸杞树多为低矮幼苗在影像中较多混杂地表土壤信息,在真彩色中往往颜色发白。中期枸杞地的枸杞树长势要比晚期枸杞地好,枝叶相对繁茂,在真彩色影像中呈深灰色或褐色,而在假彩色合成图像中部分呈现红色。而早期枸杞地多为枝繁叶茂的盛产期果树,冠层覆盖度大在真彩色影像中呈黑色且假彩色合成图像中多呈现深红色。早期枸杞地位于北部与湿地相接壤,中期枸杞地位于早期枸杞地的南部,枸杞种植区多有建筑物的出现,季节性河床中植被发育较少,防护林分布于居民地及道路周边。研究区主要地物类型野外照片如图2。
图2 野外调查实地照片Fig.2 Field survey field photos
研究流程如图3所示,主要内容包括:数据获取与预处理、分类前的准备工作(样本采集与构建特征空间)、分类特征重要性分析、协同分类精度与结果的分析与评价,最终得到了诺木洪地区地物空间分布精细制图,以分析评价多源、多特征的分类方法的可靠性与科学性。
图3 技术路线图Fig.3 Technical route map
依据实地调查数据可知,枸杞是区域内主要植被,耕地种植和草地植被覆盖度较低。故本文的分类系统主要包括8类:早期枸杞地、中期枸杞地、晚期枸杞地、防护林、裸地、建筑物、湿地和道路。样本的获取主要包括在GEE平台上选取与线下实地调查划定,根据地物的空间分布而进行选取样本点。本文共选择训练样本736个,包括野外采样点300个,对比影像选点436个,验证样本736个,包括野外采样点326个,对比影像选点410个,保证了样本均匀分布的同时也使分类效果达到最佳状态。不同地物类别及各类选取点数如表1所示。
表1 样本数据表Table 1 Sample data table
特征选择在遥感分类的过程中尤为重要,科学的选择和使用特征可以有效地提高遥感分类的精度。本文在构建特征空间时选择了波段与植被指数特征、纹理特征和极化特征。
2.2.1 波段与植被指数特征
基于预处理后的S2光学影像集合,提取原始13个波段、12个植被指数,其中包括7个红边指数,分别为:归一化植被指数 NDVI、归一化水体指数 NDWI、归一化差异耕作指数 NDTI、归一化建筑指数NDBI、修正型归一化水体指数MNDWI和对植被非常敏感的7个红边指数,所用植被指数见表2。为保证空间分辨率的一致性,其中原始波段采用最邻近插值法将Band5、Band6、Band7、Band8A、Band11、Band12由原来的20 m重采样为10 m,将 Band1、Band9、Band10由原来的60 m重采样为10 m。
2.2.2 纹理特征
由于光学影像存在“同物异谱,异物同谱”的现象,仅仅依靠光学特征会导致部分地物类型误分。纹理信息是图像在空间上以一定的形式变化而产生的信息,此次引入Haralick等提出的灰度共生矩阵(Gray Level Co-occurrence Matrix, GLCM)原理提取纹理信息。由于近红外波段针对植被的反映更加敏感,本文使用GEE提供的glcm Texture(size, kernel, average)函数,计算核大小为3×3的近红外波段纹理特征。GEE平台可计算18个纹理特征,为了减少数据冗余与及其运算速率,故选用二阶矩(asm)、对比度(con)、相关性(cor)、方差(var)、逆差距(idm)和熵(ent)6个纹理特征。
预处理之后的GF-2影像,通过野外实地调查和目视选择控制点对影像数据进行精确配准,并对GF影像进行主成分分析,以第一主成分作为纹理提取与GEE平台筛选、计算的相对应的6个纹理特征图像,并将其进行重采样至空间分辨率为10 m与Sentinel-1、Sentinel-2波段进行叠加分析,与训练样本一起输入RF分类器进行分类。
表2 文中使用的植被指数Table 2 Vegetation index used in the article
2.2.3 极化特征
本文选取了多时相的Sentinel-1雷达影像以讨论多时相雷达遥感数据对于地物分类精度的影响。在GEE上获取了研究区作物生长期内VV极化和VH极化2020年全年的时间序列影像,计算出各地物类型样本点的后向散射系数均值,研究不同作物后向散射系数随时间变化规律,并以此辅助选择研究时间范围。
为了分析多源数据和多特征组合对地物分类的影响,探讨地物分类最佳试验方案,本文针对Sentinel-1、Sentinel-2和GF-2数据,对波段和植被指数特征、纹理特征、极化特征进行4种组合方案(表3),
表3 试验方案Table 3 Experimental scheme
在现实情况下,选用多种分类特征会导致数据冗余,通常需要选择影响较大的特征以压缩数据,故本文采用随机森林的方法进行特征筛选。
为评价特征的贡献大小,本文采用基尼系数(GINI Index)作为评价指标。其GINI系数的计算公式为
式中表示类别个数,p表示节点中类别所占的比例。特征X在节点时的重要性为:
式中GI是分枝前节点的GINI系数,GI和GI是分枝后两个新节点的GINI系数。特征X在第棵树的重要性为
式中为特征X在决策树中出现的节点集合。假设随机森林中树的数目为,则特征X的特征重要性评分为
随机森林(RF)是由BREIMAN提出的基于决策树的串联集成机器的学习分类器,其原理是利用Bootstrap重采样方法从训练样本随机有放回的抽取大约2/3的样本并重复次,并为每个训练样本分别生成决策树,留下约1/3的训练样本作为Out-of-bag数据在内部交叉检验,评估随机森林的分类精度,并利用GINI系数确定决策树中每个节点的分裂条件,最后完成随机森林的构建。基于此算法有较好的分类性能,故本文采用随机森林分类器进行研究区地物分类。
基于野外验证及遥感影像对比得到地物样本点,此次采用混淆矩阵法对分类结果精度进行评价,混淆矩阵也称之为误差矩阵,用行列的矩阵表示。主要包括生产者精度(Producer Accuracy, PA)、用户精度(User accuracy, UA)、总体精度(Overall Accuracy, OA)和Kappa系数(Kappa Coefficient, KC),这些指标从不同角度对分类结果进行评价。
不同地物的光谱特征主要包括Sentinel-2多光谱数据的原始13个波段特征、12个植被指数(包括7个红边指数)。根据划分的8种典型地物样本点计算其光谱特征均值,结果如图4a所示。8种地物在原始波段中B11、B12和B4相较于其他波段的地物光谱区分度要好,而整体来看植被指数和红边指数的地物区分度要高。
对不同地物的原始波段、红边指数和植被指数(除红边指数)求差值,经统计得到3类的平均值、最小值和最大值分布图(图4b),植被指数在分类中发挥的作用要优于使用原始波段进行分类,尤其是CIre、IRECI、MTCI、NDVI、NDVIre1、NDVIre2和PSSRA植被指数,它们的光谱地物区分度高,由此可知原始波段和植被指数特征的地物光谱区分度顺序为:红边指数>植被指数(除红边指数)> 原始波段。
图4 地物光谱特征曲线与统计图Fig.4 Characteristic curve and statistical graph of ground feature spectrum
为了消除纹理特征单位或数量级之间存在的差异,对计算的纹理特征值进行了归一化处理,图5中对应的Snetinel-2和GF-2提取的纹理特征显示,Sentinel-2提取的纹理特征除了二阶矩之外,其余的各类地物纹理特征区分度较为明显,而GF-2中各地物的纹理特征仅协同性、对比度和二阶矩有一定的区分效果之外,其余的纹理特征地物区分度并不好,分析原因是GF-2原始影像空间分辨率较高,提取的纹理特征会较大成都上突出地物的具体细节纹理特征,但地物的边缘纹理特征区分度则会减弱。
图5 不同地物纹理特征曲线Fig.5 Texture characteristic curves of different ground features
为更好的了解地物在全年的变化状态,本文在GEE平台上获取了研究区的各类地物一年内的VV极化和VH极化时间序列数据,计算各种地物在不同时期的后向散射均值,以此作为极化特征进行对地物类型进行分类。
由图6不同地物极化特征变化曲线可知,多时相的Sentinel-1雷达数据可以明显区分不同地物以获取地物信息,能够作为不同地物分类的重要依据。除此之外各地物在6-9月的时间段内极化特征区分比较明显,由此也为GF-2影像在时间上的筛选提供了可靠数据依据。
图6 不同地物极化特征曲线Fig.6 Polarization characteristic curves of different features
本文通过提取波段特征、植被指数特征、纹理特征和极化特征构成特征空间,共选择33个特征参与地物分类。针对不同数据源和特征组合,基于随机森林的特征重要性排序及其对总体精度、Kappa系数的影响如图7所示。对于4种分类方案而言,其特征对总体精度和Kappa系数的影响均呈现出一个普遍的趋势:前期有一个比较迅速的上升期(黑色垂直虚线),中间逐渐变化趋势平缓,然后有一个最优的节点(黑色垂直实线,即在特征压缩的情况下保持分类结果最精准),最后再增加分类特征会导致分类结果精度下降的情况。由于本文研究重点是多源数据以及多特征的协同分类,因此,在根据特征重要性选择特征时,文中选择最优节点前的特征(优选特征集)进行地物分类。
特征优选方案1选择的特征包括:11个植被指数(7个红边指数)和10个原始波段;方案2选择的特征包括:7个植被指数(4个红边指数)、7个原始波段和3个纹理特征;方案3选择的特征包括:5个植被指数(2个红边指数)、6个原始波段、3个纹理特征和1个极化特征;方案4选择的特征包括:12个植被指数(7个红边指数)、13个原始波段、4个纹理特征和2个极化特征。根据4种方案选择的特征可知光谱特征,尤其是红边指数在地物分类中的作用较大,这也和上文中光谱特征分析相对应,纹理特征的选取也和纹理特征分析得出的结论相同,即地物的特征区分度大的特征在分类中占据重要地位。
图7 不同特征优选方案分类精度随特征变化图Fig.7 The classification accuracy of different feature selection schemes varies with features
表4和图8显示,整体上4种方案的分类结果均较好,这是由于结合了对植被分类影响比较重要的光谱指数,尤其是加入了对植被较为敏感的红边指数。
表4 不同试验方案的分类精度Table 4 Classification accuracy of different experimental schemes
方案1利用原始波段和植被指数对研究区主要地物进行分类的总体精度是95.91%,Kappa系数是0.951 1,即说明Sentinel-2的光谱信息对于该研究区的地物分类有比较好的效果,但有个别地物的PA和UA处于90%附近,考虑有误分的情况发生。
因此,方案2在方案1的基础上添加地物纹理信息,其结果总体精度达到96.57%,Kappa系数是0.956 1,分类精度均有所提高。方案3基于方案2加入了多时相的极化特征,其分类结果精度更高,总体精度为97.62%,Kappa系数达到了0.971 6。综上知,不同地物对原始波段与植被指数特征、纹理特征和极化特征有着不同的响应,利用分类结果图和分类精度表综合分析知:整体上利用原始波段以及植被指数(尤其是加入了红边植被指数)地物分类的效果已经较好,即各地类对光谱原始波段和植被指数特征的反应比较强的。加入纹理特征后各个时期的枸杞地分类精度均有提高,分析原因是枸杞地的纹理信息发育较好,地块边界纹理相对清晰。防护林和建筑物由于其在遥感影像上以线状或块状形式表现,故加入纹理特征后精度也有所提高。但是裸地、湿地和道路光谱特征差异大但纹理相似,加入纹理特征后会误分而导致精度下降。对于极化特征而言,各时期的枸杞地块对其的敏感程度比较低,利用极化特征可以反应地物的外观结构特征信息,防护林、裸地、湿地和建筑物由于具有一定的高度区分性,故加入极化特征后精度得到提高。
方案4在方案3的基础上改变提取纹理特征的数据源,其分类结果相较于方案3精度有所降低,分类总体精度为96.67%,Kappa系数达到了0.960 2。根据PA和UA显示尤其是早期枸杞地、中期枸杞地、湿地和道路的分类精度有所降低,分析原因是这几类地物的边缘纹理信息比较丰富,利用GF-2数据源提取的纹理信息会降低边缘信息的提取,进而导致分类精度降低。除此之外,后者纹理特征发生变化改变了原来的特征空间数据集,其也是分类精度降低的原因之一。
为了进一步比较各个方案的地物分类差异性,获取了分类结果的局部放大图(图9),从局部图可以显示出方案3分类精度优于其他3种方案,其湿地的提取效果明显要优越很多,防护林的提取效果基本近似,建筑物的提取也有优势,但是针对晚期枸杞地会有部分误分为裸地,其结果显示和精度验证表4也相互对应。
图8 不同方案分类整体图Fig.8 Overall diagram of different scheme classification
图9 不同方案分类局部图Fig.9 Different scheme classification partial map
特征重要性对于地物的分类尤为关键,为了分析最佳分类方案中参与地物分类的特征之间的相关度,本文采用“距离”指标进行描述,即各个特征重要性评分之间的差别。通过距离的大小反映特征之间的紧密度,距离越小越紧密,反之越松散。通过图10显示波段和植被指数特征、纹理特征和极化特征的内部特征关系紧密,而三组特征之间的大部分特征关系相对较弱。
依据最优方案3的分类结果统计了研究区内主要的8种地类的面积分布(图11),早期枸杞地、中期枸杞地、晚期枸杞地、防护林、裸地、建筑物、湿地和道路面积分别为35.257、32.058、8.332、4.507、54.834、0.637、31.024和3.497 km²。研究区主要地类便是早期枸杞地、中期枸杞地、裸地和湿地,晚期枸杞地是最近1~2 a刚种植的,考虑死亡等因素,其分布范围不大。
图10 最优方案特征相关性Fig.10 Feature correlation of optimal scheme
图11 最优方案地物面积Fig.11 Optimal scheme ground object area
利用多种遥感数据源和特征对地物类型进行分类,取长补短,充分利用各种数据的优点是得到高分类精度的前提。故本文为分析多源遥感与多特征协同会如何影响分类结果而设计了不同方案,表4中方案3(数据源:S2和S1;特征:波段、植被指数、纹理和极化特征)与方案1(数据源:S2;特征:波段和植被指数特征)的生产者精度(PA)和用户精度(UA)对比显示:道路和建筑物由于光谱反射率较高,故在加入极化和纹理特征后,其PA和UA精度变化比较小,且野外验证发现晚期枸杞地种植的枸杞高度较低,纹理特征发育比较弱,故晚期枸杞地的PA和UA精度变化也不大。而由于其他几种地物类型的纹理特征发育较好,长时序的植被生长期遥感影响特征不同,故在加入纹理特征和多时相极化特征之后分类精度有了较大的提升,早期枸杞地UA提升了0.78百分点,PA提升了1.17百分点,中期枸杞地UA和PA分别提升了4.82百分点和0.96百分点,防护林、湿地和裸地的分类精度也相应有所提升。
随着遥感影像空间分辨率的提高,纹理特征被广泛用于遥感分类由于不同地类空间尺度不同,纹理对不同地物分类的影响程度也有所差异。为了讨论不同空间分辨率提取的纹理特征参与分类对分类精度的影响,本文提取纹理特征的数据源空间分辨率分别10和1 m,为保证空间分辨率一致以进行特征组合分析,将1 m高分融合图像提取的纹理特征进行重采样至10 m,但是其纹理特征效果仍优于S2提取的纹理特征。将原纹理特征更换为高空间分辨率影像提取纹理特征不仅改变了纹理特征本身,而且还对特征空间数据集产生了一定的影响。由图7及表4可以看出,由于特征空间数据集和纹理特征的改变导致边缘纹理信息强的道路和建筑物精度有所降低,而湿地整体纹理特征不明显,在高空间分辨率数据源下提取的纹理特征破碎度比较高,湿地地块容易形成混合像元,导致其PA和UA均降低了2.87百分点和2.91百分点。早期枸杞地和中期枸杞地对于低空间分辨率提取的纹理特征,由于防护林和枸杞块间隙作为边缘纹理,而在高空间分辨率提取的纹理特征中便会缺失这些信息,导致精度的下降。晚期枸杞地发育比较弱,植株较小,防护林呈带状分布,相对较窄,二者在高空间分辨率提取的纹理特征会更明显,故二者PA和UA分别提升了0.87、2.52个百分点和0.83、0.85个百分点。
加入不同空间分辨率影像提取的纹理特征对地物的分类精度影响不同,在低空间分辨率数据源提取纹理特征进行分类下的各地类的分类精度大于高空间分辨率数据源下的。不同地物类型受纹理特征的影响程度不同,并与空间分辨率有着密切关系,因此纹理特征对多光谱遥感分类的影响并不是分辨率越高分类精度越高,在适宜的分辨率下进行分类识别才能达到更好的效果。另外,同一种植被在不同的生长期间数据对其有一定程度的差异,不同物候期植被的纹理特征不尽相同,与其他地类的区分度也有差异,纹理特征在不同物候期对植被分类的影响需进一步研究。
1)基于随机森林的特征重要性评估方法可以有效地提取出对地物分类影响比较大的特征,极大地压缩了数据量,减少数据冗余,提高数据处理效率。针对最优方案3,其特征由33个减少到15个,并且也是在所有特征按重要性排序组合中分类精度最高的。此外,光谱特征、纹理特征和极化特征的内部特征关系紧密,而3组特征之间的大部分特征关系相对较弱。
2)通过4种分类方案的对比分析,植被指数(尤其是红边指数)在所有特征中对分类结果的影响比较大,纹理特征次之,极化特征的影响相对较小。结果符合研究区地物特征,其研究区主要地物的光谱特征差异比较明显,而纹理特征由于部分植被分布比较散乱,进而导致纹理信息不丰富,极化特征虽然可以有效地区分光谱特征相似的地类,但是对于光谱差异较大的地类,极化特征的作用并不明显。
3)协同Sentinel-2与Sentinel-1数据源,在基于光谱特征分类的基础上加入纹理特征和极化特征可以有效地提高地物的分类精度,其总体精度提高了1.71个百分点。协同三种数据源,利用全特征空间,使用Sentinel-2提取纹理特征的各地类分类精度优于使用GF-2提取纹理特征参与分类的情况,不同空间分辨率图像提取的纹理特征对特征优选方案及分类结果有着不同影响,在适宜的分辨率下提取纹理特征参与分类才能达到更好的效果。