利用光学和SAR遥感数据的若尔盖湿地植被分类与变化监测
2023-07-13明義森刘启航柏荷黄昌
明義森,刘启航,柏荷,黄昌
西北大学 陕西省地表系统与环境承载力重点实验室,西安 710127
1 引言
湿地含有丰富的有机碳资源,与其他土壤碳库相比具有面积小但储量大的特点(Davidson and Janssens,2006)。植被是湿地生态系统的重要组成部分之一,是湿地碳固持的重要角色,其生长变化能反应湿地系统发展的成熟程度(田应兵,2005)。因此,植被的分类与制图是监测湿地生态系统变化的关键环节(van der Putten,1997)。湿地系统在空间分布上具有不均匀性和复杂性的特点,随着旱季与雨季的交替,湿地具有季节性变化或者年变化的动态特征,所以及时、准确的监测湿地植被类型的分布和变化是一个挑战。
若尔盖地区发育有丰富的湿地,其泥炭资源丰富、碳储量巨大(田应兵,2005)。在若尔盖地区,土壤碳储量与含水量分布规律总体一致,含水量越高,土壤有机碳储量越高(Lee 等,2009),同时,土壤含水量也影响着不同类型植被的生长。因此,若尔盖湿地植被的空间分布一定程度上能反映其碳固持能力的空间分异,对于研究中国甚至全球碳循环有非常重要的意义。
遥感技术为湿地植被的分类与变化监测提供了先进、高效的手段。光学遥感数据直观、易获取和处理,能有效识别湿地植被(童庆禧 等,1997)。Sun 等(2021)基于Sentinel-2 的时间序列影像对湿地植被进行分类,实现了较高的分类精度。但是,受云等天气条件限制,Sentinel-2影像的时间序列往往不够完整。受影像质量限制,张磊等(2019)通过综合一年内的多幅Sentinel-2 图像,成功提取了黄河三角洲湿地信息。合成孔径雷达SAR(Synthetic Aperture Radar),作为先进的微波遥感手段,具有穿透性强、全天时、全天候的特点,有完整的时间序列数据。SAR 已经被证实能对湿地植被的分类以及变化监测提供帮助(White 等,2015),如Tsyganskaya 等(2018)使用Sentinel-1 SAR 时间序列数据对湿地植被分类。但是SAR 遥感图像仅能反映地物的后向散射特征,无法准确实现地物对象的识别和分类(Bourgeau-Chavez等,2016)。因此,结合光学遥感与SAR 遥感,有望发挥各自的优势,更好地对湿地植被进行分类与制图。目前,利用机器学习算法结合光学与雷达遥感数据在遥感图像分类方面已得到广泛应用,其中随机森林算法被认为是分类精度最高的算法之一(Tian 等,2016;Franklin 等,2018;刘瑞清 等,2021)。Ruiz 等(2021)结合Sentienl-1和Sentinel-2 数据使用随机森林算法成功实现亚热带湿地植被物种的分类,Slagter 等(2020)从Sentinel-1 SAR 和Sentinel-2 多光谱的时间序列数据提取特征对南非Lucia 湿地植被进行分类,但没有充分发挥Sentinel-1 SAR 数据多时相的特点。动态时间归整算法DTW(Dynamic Time Warping)是一种时间序列分析算法,可以计算两个时间序列间的相似性以辅助分类,特别适合对植被及农作物等具有物候特征地物的分类。
训练样本是遥感图像监督分类的关键,决定了分类质量的好坏。由于大部分实地观测样本是在特定时间获得的,而湿地植被可能会随着时间变化而发生改变,这对需要基于多期遥感影像进行分类的植被变化监测带来了困难。Huang等(2020)提出了一种自动训练样本迁移的方法,通过测量参考光谱与目标光谱的光谱相似度和光谱距离对1990年—2010 年的训练样本像素的变化状态进行检测和识别,取得很好的效果。样本迁移的理论使得我们能在训练样本有限的情况下稳定分类,减少了样本采集方面付出的精力、成本,为快速的湿地植被分类制图,尤其是时间和空间尺度上训练样本不可得条件下的植被分类,提供了有效方法(Gong 等,2019)。Yan和Niu(2021)利用样本迁移方法对青藏高原湿地类型进行分类,表明样本集的迁移可以快速获取准确的湿地样本集,为时间序列、多时相的湿地制图的研究奠定了基础。
因此,本文基于Sentinel-1 SAR 影像的时间序列数据与Sentinel-2 光学影像数据,充分发挥光学影像多光谱的优势和SAR 不受天气条件限制便于构建完整时间序列的优势,利用2020 年实地获得的高分辨率无人机样本,并使用样本迁移方法将其迁移至2017年,借助Google Earth Engine(GEE)平台,通过随机森林算法对若尔盖湿地2017 年和2020 年湿地植被进行分类,并对期间的植被类型变化进行检测和评估。
2 研究区和数据
2.1 研究区
若尔盖湿地(33°00'N—34°10'N,101°40'E—103°20'E)是中国3 大湿地之一,横跨甘肃省和四川省,包含四川省若尔盖县、红原县、阿坝县、松潘县以及甘肃省的玛曲县、碌曲县(刘红玉和白云芳,2006)(图1),是中国特有的高原湿地和高寒湿地生态系统的典型代表。研究区内湿地面积约4400 km2。若尔盖高原湿地平均海拔3400—3900 m,素有黄河“蓄水池”之称,是黄河上游地区重要的水源涵养地。若尔盖湿地植被主要可以分为3 类:沼泽植被、沼泽化草甸植被、草甸植被(田应兵,2005)(表1)。
表1 若尔盖湿地植被类型及主要种类Table 1 Vegetation types and main species of Zoige Wetland
图1 研究区范围(湿地范围来自于Mao等(2020))Fig.1 Study area(Wetland extent was from Mao et al.(2020))
2.2 数据
2.2.1 无人机野外实地调查数据
在2020 年7 月利用大疆Mavic 2 无人机在若尔盖湿地内选取典型样区进行航摄,现场天气状况良好,晴空无云,风速较小。飞机飞行高度控制在100 m 左右。将获得的无人机航空图像利用Pixel4D软件进行处理,最终生成0.5 m分辨率,含红、绿、蓝3个可见光波段的正射影像。结合现场观察与照片,利用目视解译对这些正射影像中的植被类型进行判读,共获得4428 个样本点,其中草甸植被1337 个、沼泽化草甸植被1296 个,沼泽植被1795 个。将其随机划分,60%作为训练样本集,40%作为测试集。
2.2.2 哨兵1号卫星影像
哨兵1 号卫星(Sentinel-1)搭载了C 波段合成孔径雷达,其入射角稳定,便于基于同一轨道数据开展多时相的研究。在GEE 平台上获取了2017年和2020年相同轨道高分辨率IW模式降轨双极化(VV+VH)数据,共104 幅。这些图像已经进行过预处理(包括热噪声去除、辐射校准等),并且在GEE平台上使用3×3 Sigma Lee滤波器(Lee等,2009)进行平滑处理。
2.2.3 哨兵2号卫星影像
在GEE 平台上筛选出哨兵2 号卫星(Sentinel-2)在2017年和2020年云量比例小于8%的全波段幅影像,这些影像是经过正射校正和几何精校正的产品。利用Sentinel-2的质量检查波段进行去云处理,分别合成2017年和2020年无云全波段中值影像。
3 研究方法
基于无人机影像目视解译得到训练样本集,提取训练样本集在2017 年和2020 年的光谱信息和时间序列后向散射特征,在光谱相似性样本迁移的理论基础上,将2020 年的植被样本迁移到2017 年。首先,选择Sentinel-2 的植被红边波段、近红外波段、短波红外波段作为光谱特征;然后,基于DTW 算法计算3 种若尔盖湿地植被在一年内Sentinel-1 影像VV、VH 波段的时间序列相似度作为微波特征集。根据Mao 等(2020)发布的中国湿地地图确定研究区湿地范围,利用随机森林算法对2017 年和2020 年的湿地植被进行分类并用混淆矩阵、Kappa 系数等评价分类精度。最后,对2017 年—2020 年湿地植被的变化进行监测。流程图如图2所示。
图2 研究方法流程图Fig.2 Work flow of this study
3.1 构建光学—微波特征集
3.1.1 构建光学特征集
Amani等(2017)发现红边波段对于湿地植被分类的精度至关重要,这主要是因为该波段的反射率与植被生物参数,如叶绿素含量和植被的水分含量等(Mutanga 等,2012)关系密切。此外,短波红外波段、近红外波段也能有效提高湿地植被分类精度,常用于湿地植被分类(Bajgain 等,2015;Mobasheri 和Amani,2016)。因此本文选取了Sentinel-2 的1 个近红外波段(B8)、2 个短波红外波段(B11、B12)以及4 个植被红边波段(B5、B6、B7、B8A)共7个波段作为光学特征集。
3.1.2 构建微波特征集
由于若尔盖湿地植被一般从4月中下旬开始返青,10 月份进入普遍枯黄期。因此本文筛选出2020 年4 月—11 月的全部Sentinel-1 影像,提取3 种植被类型样本点的后向散射强度平均值,通过其在一年内不同时间VV、VH 波段变化表示3种植被类型的物候特征(图3)。沼泽化草甸植被VV、VH 波段的后向散射强度在一年内全部时期均大于草甸植被且变化趋势相同。与草甸植被和沼泽化草甸植被相比,沼泽植被的后向散射强度在一年内的波动较大,这与它的生理环境有关。沼泽植被处于常年积水或季节性积水的环境,当沼泽植被冠层下水分的增加,它的后向散射强度也会增加(Kasischke 等,1997)。而植被冠层下水分因降雨影响会有较大的波动,这使得沼泽植被年内后向散射强度变化起伏较大。因此,由SAR 影像构成的物候特征能在一定程度上反映不同植被类型之间的差异,提高植被类别间的可分性。本研究以样本点时间变化曲线为参考,利用DTW 算法计算待分类像元与样本曲线之间的距离。
图3 3种植被类型一年内后向散射强度时间序列变化Fig.3 Time series changes of backscattering intensity for three types of vegetation within one year
DTW 算法是将一个参考的时间序列与一个未知的时间序列进行比较,计算两个时间序列之间的距离(Maus等,2016),一般使用欧氏距离度量两个时间序列之间的相似性,计算公式为
式中,ui∈U(i=1,2,…,n),vj∈V(j=1,2,…,m),U={u1,u2,…,un}和V={v1,v2,…,vm}分别是长度为n和m的时间序列,然后用n×m的矩阵Dbase=(dbase(ui,vj))n×m存储。最后计算矩阵D的最小距离递归和得到DTW距离,公式为
3.2 样本迁移
由于野外无人机观测实验开展于2020 年,无法获得2017 年的训练样本。因此,借鉴Huang 等(2020)的样本迁移方法,将2020年的训练样本迁移至2017 年使用,作为其分类的训练集和验证的测试集。图像的光谱特征是分类的主要依据,当两个样本的光谱变化较小时,可以认为样本的类型没有发生变化,光谱的变化程度用参考样本与目标样本间的欧氏距离和角距离度量。在光谱空间中,每一个像素对应一个多维光谱向量,即计算两个向量间的欧氏距离和角距离。欧氏距离ED(Euclidean Distance)越大代表差异程度越大,如式(3)所示,但ED 对亮度值比较敏感。光谱角距离(SAD)是指两个光谱向量间的夹角的余弦(式(4))。SAD 越大代表相似度越高。当SAD 等于1时代表两个样本完全没有差异。与欧氏距离不同,角距离对亮度不敏感,即当亮度值增加或者减少时,角距离保持不变,可以突出光谱的特征。
式中,X,Y为参考光谱与目标光谱,N为总波段数。
3.3 随机森林算法
基于2020 年的原始训练样本和迁移得到的2017 年训练样本,本研究使用随机森林算法分别基于2017 年和2020 年所构建的光学—微波特征集数据进行植被分类,获得2017 年和2020 年若尔盖地区湿地植被分类结果。随机森林是基于bagging框架的决策树模型的一种机器学习算法。与单一分类器相比,它更准确,对噪声的反应更好(Breiman,2001)。随机森林中的每个分类器是通过随机有放回地从训练集中抽取训练样本进行训练而生成的。
随机森林算法采用信息熵和信息增益作为确定特征选取顺序的依据(Hssina 等,2014)。信息熵指的是信息的混乱程度,即熵的值越大,信息越混乱。信息增益Gain(p,T)计算所有样本数据的类的混合程度以及决策树的任意位置,它定义了检验特征T和位置P的增益。我们可以使用Gain(p,T)对全部特征进行排序,并构建决策树,其中每个节点的位置是根据尚未考虑的特征中具有最高信息增益的那个特征决定的。一个随机森林由一组分类器组成,每个分类器对数据做出投票,最后将得票最高的类别分配。以无人机影像获得的测试集为参考,利用混淆矩阵、总体精度、用户精度、生产者精度、Kappa系数来评价分类结果的精度(吉海彦 等,2019)。
4 结果与讨论
4.1 2020年植被分类结果及精度评价
2020 年若尔盖湿地植被分类结果如图4 所示,其中河流湖泊结果是从Mao 等(2020)数据集中获取。湿地植被在研究区分布广泛,沼泽植被大面积生长在研究区的东部,沼泽化草甸植被和草甸植被呈现普遍的分布状态。在近河床区域主要生长着沼泽化草甸植被和草甸植被,而沼泽植被较少,与野外实地观察结果一致如图4(Ⅰ)所示。在图4(Ⅱ),3种植被分布从边缘到中心呈带状分布,草甸植被、沼泽化草甸植被分布在边缘而湿地植被生长在中心区域。这可能与局部水文条件、泥炭层厚度有关;离湖泊近的区域沼泽植被面积较大,而沼泽化草甸植被、草甸植被一般离湖泊较远(图4(Ⅲ))。
图4 2020年若尔盖植被分类图Fig.4 Ruoergai vegetation classification map in 2020
表2 混淆矩阵可知,2020 年的植被分类总精度达到了97.43%,Kappa系数为0.96。其中沼泽化草甸植被的分类精度较差,与草甸植被与沼泽植被都出现了一定程度的混分。这是由沼泽化草甸的特性导致的。沼泽化草甸作为过渡性植被类型,它的光谱特征与其他两种植被相似,因此容易出现混分。
表2 2020年植被分类精度评价Table 2 Evaluation of vegetation classification accuracy in 2020
4.2 采用样本迁移的2017年植被分类结果
根据样本迁移原理,以2020 年影像为参考,2017年影像为目标,计算参考光谱与目标光谱之间距离以评价参考年份与目标年份样本的相似性。本研究中,样本迁移分为两个方面,首先,基于光谱相似性的光学特征样本迁移,计算参考光谱与目标光谱间的欧氏距离和角距离评价参考年份与目标年份样本的光谱相似度(图5(a),图5(b))。其次,基于DTW 的微波特征样本迁移,计算目标年份与参考年份VV、VH 波段时间序列数据的DTW 距离以评价样本点在目标年份与参考年份的相似度(图5(c),图5(d))。最后,将两个迁移结果取交集,得到2017年的植被样本集。
图5 目标样本与参考样本相似性评价Fig.5 Evaluation of the similarity between the target sample and the reference sample
本研究参考Huang等(2020)的方法,结合研究区训练样本的统计数据,配合Google Earth 高分辨率影像及无人机影像等,通过选定适当的阈值保留光谱变化较小的样本进行迁移。最终确定,对于样本的光学特征ED<0.40 且SAD>0.93,同时微波特征满足VV 波段DTW 距离<30 且VH 波段DTW<35样本点认定为满足迁移条件,认定为可迁移样本,共3480 个。其中有1078 个草甸植被样本、793 个沼泽化草甸样本和1609 个沼泽植被样本。将所有可迁移样本随机划分,其中60%作为训练集,40%作为验证集。
样本迁移能快速解决多时相分类中训练样本有限的问题,减少采集样本的成本与精力。图6以核密度曲线的形式展示了样本迁移前后样本在光学、微波特征的分布变化,发现样本数据在全部特征的分布中都有一定程度的集中。在3种植被类型中,沼泽植被核密度曲线变化最大,特别是在红边波段与基于DTW 计算的VV 波段;而草甸植被和沼泽化草甸植被变化不明显。样本迁移后去除了光谱变化较大,具有差异明显的样本点,保证了样本的纯净度。
图6 样本迁移前后样本核密度曲线变化Fig.6 The change of kernel density estimation after the migration of training samples
利用随机森林方法,基于迁移样本中的训练集进行2017 年若尔盖植被类型进行分类,混淆矩阵如表3 所示,总体精度达96.84%,Kappa 系数为0.95。分类结果如图7,植被类型分布规律总体与2020年相似。
表3 2017年植被分类精度评价Table 3 Evaluation of vegetation classification accuracy in 2017
图7 2017年若尔盖植被分类图Fig.7 Ruoergai vegetation classification map in 2017
4.3 2017年—2020年若尔盖湿地植被变化评估
沼泽植被处于常年积水或季节性积水中,它是3种植被类型中土壤含水率最高的,与另外两种植被相比,它的变化对湿地初级生产力、固碳能力有重要影响。对沼泽植被变化的评估能够反映若尔盖湿地的退化、恢复情况。将沼泽植被演变为沼泽化草甸植被、草甸植被的过程定义为退化演替过程;而草甸植被、沼泽化草甸演变为沼泽植被的过程定义为恢复演替过程(唐艳梅 等,2021)。
依据演替方向进行分类得到2017 年—2020 年植被演替变化图(图8)。研究区内大部分湿地植被保持原本类型没有发生变化,尤其是在河流附近(图8(Ⅰ)),草甸、沼泽化草甸植被类型变化很少,这可能由于河流有较好排水性,保持附近土壤含水量在一定水平,使得植被没有发生演替;恢复演替主要发生在大面积生长着沼泽植被的附近,尤其是沼泽植被分布的边缘地区(图8(Ⅱ));湖泊周围有大面积的恢复演替发生(图8(Ⅲ)),这可能与湖泊水位变化有关;而退化演替在研究区内普遍发生,这可能与降雨、降水年变化有关。
图8 沼泽植被变化分布((Ⅰ)、(Ⅱ)、(Ⅲ)是沼泽植被变化代表性区域。左下角直方图表示发生退化演替与恢复演替的面积。其中红色代表恢复演替,黄色代表退化演替)Fig.8 Distribution of wetland vegetation change((Ⅰ),(Ⅱ),and(Ⅲ)are the zoom-in areas of swamp vegetation change.The histogram in the lower left corner shows the areas of degeneration and restoration.The red represents the restoration,and the yellow represents the degeneration)
对退化、恢复演替区域的面积进行统计发现(图8 中直方图),植被退化演替区域约100.17 km2(占湿地总面积的2%),而被恢复演替区域面积约为331.89 km2(占湿地总面积的7%)。若尔盖湿地植被在近年来没有发生较大变化,发生变化的区域中,恢复演替面积大于退化演替面积,说明若尔盖湿地系统近年来呈向好趋势。
5 结论
本文通过从2020 年实地获得的高分辨率无人机数据中提取训练样本,并将其迁移至2017 年,利用随机森林方法结合Sentienl-2 光学数据与Sentinel-1 SAR 数据对若尔盖湿地植被进行分类,最终分别得到2020 年和2017 年的湿地植被分类图,并进一步开展植被变化监测和湿地植被退化/恢复分析,得到以下结论:
(1)研究发现3种湿地植被的后向散射值一年内变化曲线具有可分性,说明SAR 特征的引入有利3 种植被类型的识别。结合Sentinel-2 光学影像和Sentinel-1 SAR 对若尔盖湿地主要的3 种湿地植被进行分类,可以发挥前者多光谱的优势和后者时间序列完整的特点,分类总体精度达到97.43%,Kappa 系数为0.96。(2)引入样本迁移思想,解决了在变化监测研究中特定时刻获取的样本数据无法直接应用到其他时刻的问题。本研究通过将2020 年实地采集的样本迁移至2017 年,解决了历史时期实地样本不可得的问题,顺利实现了2017 年的植被分类过程。(3)通过对2017 年—2020 年若尔盖湿地植被类型变化进行分析,发现在河流附近生长的草甸、沼泽化草甸植被变化较少,而在沼泽植被分布的边缘地区是恢复演替发生的主要区域。总体上,若尔盖湿地近年来植被类型变化不大,变化的部分主要以恢复演替为主。
需要注意的是,虽然从训练样本的生产、方法的使用与结果的验证等方面进行了严格的质量控制和评价,本文的分类结果不可避免仍存在一定的不确定性。比如,在通过高分辨率的无人机影像进行目视解译获得训练样本的时候,可能引入一定的误差;在样本迁移设置迁移规则的时候,也有可能对迁移后的样本带入一些误差。这类误差都将直接给最终的分类结果带来不确定性。此外,为了充分使用优质遥感数据,本研究选用最高分辨率可达10 m 的Sentinel-1和Sentinel-2数据,这也使得我们的研究时段局限于2017年—2020年。事实上,本文所使用的方法亦可应用于其它光学和SAR 数据,如Landsat和JERS,以实现更长时间尺度的湿地植被变化监测。
志 谢感谢欧洲航天局(ESA)免费公开的哨兵1号和哨兵2号数据,使得本研究得以顺利开展。