湿地是位于陆生生态系统和水生生态系统之间的至关重要的过渡地带,它有着贮蓄水量、调节气候、保护环境和保护生物多样性等重要作用,是世界上最有价值的自然资源之一^[1-2]。水是湿地三要素之一,是维持湿地的主要因素^[3]。由于人类活动过程中长期忽视对地球环境的保护,地球生态环境逐渐恶化,水循环系统遭到破坏,全球湿地萎缩和功能退化成为突出问题^[4]。因此,为了防止湿地进一步流失,更好地管理和保护现有的湿地资源,有必要研究湿地水域变化,监测其年内、年际变化。

目前国内外湿地水体遥感提取方式主要有以下3种: ①基于传统光学遥感图像的遥感技术,通过构造波段间的光谱关系提取水体^[5]。如甄佳宁^[6]使用Landsat系列卫星影像对长春湿地采取不同的水体指数和阈值分割法提取水体信息; Jin等^[7]利用Landsat时间序列影像监测德尔马瓦半岛的湿地水位淹没动态,增进对湿地的了解并确保生态系统服务的可持续性。不足之处是光学遥感图像的遥感技术受云雨天气等条件限制,对监测湿地水域年内变化影响颇大^[8]。②基于微波遥感图像的遥感技术,通过合成孔径雷达(synthetic aperture Radar,SAR)后向散射系数(σ⁰)、干涉相干系数(μ⁰),构建SAR影像数据集提取水体信息。徐怡波^[9]基于多时相ENVISAT卫星影像分析水体的后向散射系数特征,采用阈值分割的方法提取水体,实现对水面变化信息的监测; 付波霖等^[10]利用多频率、多时相SAR影像,采用合成孔径雷达差分干涉测量技术(interferometric synthetic aperture Radar,InSAR),构建沼泽水位变化监测模型,探究了不同时相沼泽湿地干涉相干性的差异。不足之处是使用单一的σ⁰或单一的μ⁰方式有一定的局限性。③基于微波和光学遥感图像的融合,利用2种传感器的互补性提取水体信息^[11]。关韵桐^[12]利用多时相Sentinel-1 雷达影像为主数据源,以Sentinel-2和Landsat8 OLI 光学影像为辅数据源,将光学影像计算出的归一化植被指数(normalized difference vegetation index, NDVI)分类,根据雷达影像数据分析了σ⁰与雷达参数之间的关系,利用反演模型对大山包高原湿地进行土壤水分反演,获取了土壤的水分分布状况,根据实测数据验证了反演结果的准确性; Mizuochia等^[13]使用AMSR系列、MODIS和Landsat卫星数据,开发并评估了基于查找表的时空数据融合方法,并将其应用于半纳米比亚的干旱季节性湿地,结果表明时间序列图能准确地描述湿地的季节变化和年际变化。综上,目前缺乏整合长时间序列的微波遥感影像的后向散射系数影像和干涉相干系数影像方法,构建多源遥感数据集来检测湿地水域面积动态变化方面的研究。

本文以东北扎龙国家级自然保护区作为研究区,以2018—2019年Sentinel-1A卫星影像数据作为主数据源,辅以Sentinel-2A卫星影像数据,提取保护区年内和年际的SAR后向散射系数和干涉相干系数,构建湿地水面积提取SAR影像数据集,利用机器学习算法构建沼泽湿地水面动态变化监测模型,实现扎龙湿地水域面积的动态变化监测,探究湿地水域年内和年际变化规律。

扎龙湿地(图1)位于黑龙江省松嫩平原西部的乌裕尔河下游,地理坐标为E123°47'~124°37',N46°52'~47°32',边界平面呈现不规则橄榄形,湿地南北约80.6 km,东西约58.0 km,总面积约为2 100 km²,是面积排名亚洲第一、世界第四的国家湿地生态系统保护区,也是世界上面积最大的芦苇湿地,它于1987年内被国务院批准为国家级自然保护区^[14-15]。多年来扎龙湿地面临水资源短缺、水位下降、资源枯竭、面积分解等问题,严重影响了湿地的生态环境健康,甚至影响了当地经济和社会的可持续发展。

本论文使用的数据主要为Sentinel-1A和Sentinel-2A数据,其中以Sentinel-1A的VV极化数据为主要实验数据,Sentinel-2A多光谱影像辅以验证水体提取精度,在欧洲航天局官网(https://scihub.copernicus.eu/)中公开获取。具体获取了东北扎龙国家级自然保护区2018—2019年Sentinel-1A数据逐月的IW模式的单视复数影像,在研究区和时间段内,总共获取了44景Sentinel-1A影像,为了验证本论文实验精度,获取了4景空间分辨率为10 m的Sentinel-2A影像,具体影像日期如表1所示。

本文主要利用SNAP软件对多时序Sentinel-1卫星数据进行一系列预处理,主要包括热噪声去除、辐射定标(将数字像素值转换为后向散射系数)、条带修复(去除条带噪声)、滤波处理(Define Lee滤波)、多视处理(距离向视数设为4)、斜距转地距及地理编码获取UTM WGS 1984投影坐标系的σ⁰数据,如图2(a)所示; 经过应用轨道文件(获取准确的卫星位置和速度信息)、影像配准(将主影像和辅影像进行配准获取2个时相影像相融合的新影像)、增强光谱多样性(校正条带的距离和方位向)、干涉条纹图的相干系数提取(对原始图像参数赋予相干系数)、条带修复、多视处理(距离向视数设为4)、地形校正和影像裁剪获取μ⁰数据,如图2(b)所示。本文选择每个月靠近月中的VV极化Sentinel-1A/B影像提取研究区σ⁰数据,选择每个月2景VV极化的Sentinel-1A/B影像计算逐月的μ⁰数据。

本文总体技术方案如图3所示。第一部分为预处理模块,对原始Sentinel-1A卫星影像数据经过热噪声去除、辐射定标等预处理过程,获取研究区σ⁰影像; 经过应用轨道文件、影像配准等预处理过程获取研究区μ⁰影像; 经过波段运算方法加权研究区σ⁰影像和μ⁰影像获取新的加权影像。第二部分为数据处理模块,将第一部分预处理模块得到的3种类型影像对照彩色光学影像找到水体大致位置,使用阈值分割法确定最佳阈值分割参数提取湿地水体,逐月计算提取的湿地水体面积,获取水位动态信息。第三部分为数据分析模块,对所得影像数据基础建立湿地水面积提取SAR影像数据集,使用随机森林算法对多源遥感数据集进行湿地水体和非水体分类,计算水体面积,获取水位动态信息。第四部分为精度验证模块,对所的结果进行定量分析,探究湿地水域年内年际变化规律,验证精度。

光在水体中传播时会发生散射,入射能量与表面相互作用后,再辐射向各个方向,其中散射光沿入射方向相反方向的散射的过程称为后向散射,它也是雷达接收机接收到的散射,通常以σ⁰表示^[16]。

σ⁰是水体的固有光学量,水体在微波范围内的σ⁰较低,雷达波的吸收相对完整,在雷达图像上呈现较暗的黑色区域,水体和非水体边界明显,可通过阈值分割等多种方法轻易将水体和其他地物区分开。其中,σ⁰是体散射函数β(θ)的散射角θ对后向半球的积分,在数学上定义表达式为:

式中λ(θ)为σ⁰与体散射函数关系的常数,该常数用于此研究领域近似恒定。

本文在ENVI中实现阈值分割湿地水体提取。阈值分割是一种操作简易且效率高的图像分割方法,该方法使用一个或多个阈值来分割图像的灰度级,并且在相同阈值范围内的灰度值将被识别为相同类型的地物。本文将σ⁰设为F(x,y),以某个灰度范围[k1,k2]作为阈值范围进行分割,分割后生成的图像为f(x,y),在阈值范围内的像素全部更改为1,其他像素均更改为0,即划分灰度级后的图像为二值图,在数学上定义为:

μ⁰是描述2幅单视复数影像相干程度高低的参数,即地物2次拍摄的差异,该系数可以表明某些理化特征,因此可以用作辅助地物判别的独立参数。由于水体的μ⁰值为0,故在图上水体呈现暗黑色,极易于区分水体和其他地物。对于单视复数影像,使用SNAP的相干变化检测处理方法来生成地理编码的相干系数图,干涉测量的相干系数介于0~1之间(0表示不相干,1表示完全相干),假定s₁和s₂分别为同一分辨率的一定时间间隔内获取的2幅影像,则μ⁰在数学上定义表达式为:

式中:μ0为相干系数; s^*为s的共轭复数; s₁和s₂为经过预处理的相干系数影像。

本文在ENVI软件中对μ⁰数据进行阈值分割,将图像二值化,输出水体的栅格文件,计算水体面积。由于湿地水域的相干系数值会明显低于周围地物,水体在μ⁰数据上表现出的灰度值明显偏低,因此很好选定最佳阈值参数对水体进行提取。

单一的σ⁰和μ⁰数据提取水体都分别具有各自优势,同时也存在各自的局限性。本文利用波段运算将二者分别赋予权重叠加融合,便可取得继承二者特性的加权σ⁰和μ⁰数据的新影像数据。权重的赋予是根据奥维地图2018年、2019年的高空间分辨率历史光学遥感影像来验证单一的σ⁰和单一的μ⁰提取水体精度进行权重划分,如σ⁰影像提取水体效果更佳,则权重愈大,反之愈小。

目前,加权σ⁰和μ⁰数据可以通过ENVI进行波段之间的运算实现。波段运算是遥感信息提取的重要方法之一,是基于研究对象的多光谱波段反射特性,选择合适的波段进行运算,目的是增加对象与背景的差异^[17]。

本文在ENVI软件中对新影像采用阈值分割法提取湿地水体,通过计算像元数和像元大小的关系即可获知湿地水域面积,如表2所示。

随机森林是一种新颖的、极其灵活的机器学习算法,通过集成学习的思想将多棵决策树集成,其中每棵决策树都是一个分类器,即对于一个输入样本,N棵树都将具有N个分类结果。随机森林将集成的全部分类结果进行随机而有放回的抽样,通过投票机制,最终输出投票最多的类别,该分类方法在目前所有分类的算法中,具有非常高的准确率^[18]。

本文将经过处理的Sentinel-1A构建湿地水面积提取SAR影像数据集,其数据构成如表3所示,波段组成示例图如图4所示。每个数据集由5个波段组成,对其进行湿地水体提取,采用ENVI的拓展模块随机森林分类器绘制训练样本和预测变量的随机子集,由于决策树数目过少会导致分类误差大,过大则会导致随机森林构建时间过长,故本文生成1 000棵决策树形成整体分类,表3中展示2019年10月SAR影像数据集对照奥维地图彩色光学影像选取的地表真实感兴趣区域(region of interest,ROI)样本信息,选取水体和非水体样本数各100个(70%为训练样本,30%为验证样本)。

自动化水体提取指数(automated water extraction index,AWEI)是一种可以提高包括阴影和暗表面在内的区域的分类精度的指数^[19]。由于AWEI方程通过应用强制非水像素低于0和水像素高于0的系数来提高水和非水像素的可分离性,可以采用默认的0阈值来自动客观地提取水体,而不采用人为选取的主观最佳阈值^[20]。表达式为:

式中: ρ_G为绿光波段; ρ_B为蓝光波段; ρ_N为近红外波段; ρ_M为中红外波段; ρ_S为热红外波段; AWEInsh为没有阴影下提取的自动化水体提取指数; AWEIsh为阴影下提取的自动化水体提取指数。

本文在ENVI中通过波段运算计算光学影像数据获取AWEI指数影像,将其作为本论文湿地水体提取基准,通过阈值分割获取水体信息,分析验证前文4种基于不同处理方法得到的影像数据的水体提取精度。

本文选取了2019年的3月、8月、10月和12月的Sentinel-2A卫星的影像进行水体提取方法对比实验,通过计算Sentine-2A影像数据获取4个月份的AWEI指数,以研究区AWEI数据作为本文湿地水体提取基准,分析验证前文4种基于不同处理方法得到的影像数据的水体提取精度。结果图如图5所示。

表4为2019年4期不同方法水体提取面积对比。可以得知,随机森林算法的分类精度最高; 其次为基于加权σ⁰与μ⁰影像使用阈值分割方法的分类精度; 再是基于σ⁰影像使用阈值分割方法的分类精度; 最后是基于μ⁰影像使用阈值分割方法的分类精度。

通过阈值分割法对σ⁰数据进行水体提取,阈值设定范围在0~0.01。利用2018—2019年逐月的Sentinel-1A影像提取的扎龙湿地水域面积如表4所示,变化趋势如图6所示。

结合图6(a)及表5可知,2018年1月,扎龙湿地水域面积为1 503.50 km²,2018年2—5月,湿地面积没有明显的变化,基本在1 400 km²左右,呈轻微波动性状态。从6月份开始,湿地水域面积发生巨大变化,6月份面积陡降至年内最低谷770.51 km²,较上个月下降约640 km²,下降了45.39%。7—11月份整体上稍有回升,平均面积约为1 027 km²,其中,7月份水域面积增幅较小,仅增至839.16 km²。8月份增幅较大,增至1 098.05 km²,9月份湿地水域面积呈现出略微的下降,降至1 059.18 km²,10月份水域面积有些许回升,升至1 177.69 km²,11月份又降至964.1 km²,12月份水域面积快速回升至1 609.32 km²,这个阶段为2018年湿地水域面积增长速度最快的时期,达到年内最高峰,相比于年内最低的6月份,增加了108.86%。

2019年1月,扎龙湿地水域面积为1 361.28 km²,2月份面积相对稳定,有少许回升达1 401.55 km²,3月份持续上升至1 598.70 km²。4—7月水域面积急剧下降,7月份降至年内最低谷,仅有768.11 km²,与年内最高峰的3月份相比,7月份下降830.59 km²,下降了51.95%,如图7(a)所示。8—12月整体呈高速度增长趋势,10月之时水域面积已达1 590.14 km²,11月份与10月相比较为稳定,,12月份呈现出略微下降趋势,降至1 541.74 km²,如图7(b)所示。

结合图6(b)及表5可知,2018年湿地水域面积与2019年相比,1月份2018年多出了142.22 km²,而与之相反,2月份2019年多出了89.76 km²,同样的3月份中2019年的明显较多,多出了155.92 km²,4月份则表现出相对持平状态,仅差4.48 km²。5月份、7月份和8月份2018年比2019年分别多出244.59 km²,71.05 km²和103.9 km²,10月份呈现反转,2019年多出了412.45 km²。11月份出现了年际差别最大的峰值,2019年水域面积多出了631.8 km²,如图8所示。最后12月份时2018年比2019年稍多177.23 km²。

据有效数据,2018年扎龙湿地月平均水体面积为1 269.03 km²,2019年扎龙湿地月平均水体面积为1 334.65 km²。总体上看,2019年扎龙湿地水量比2018年更为充足。可以看出,虽然扎龙湿地水域面积有明显的季节性变化,但每年都有其变化的特征,会随着降雨量、温度及人为补水等多种因素而出现独特的现象。

通过阈值分割法对μ⁰数据进行水体提取,阈值设定范围在0~1。获取2018—2019年扎龙湿地逐月水域面积如表6所示,变化趋势如图9所示。

结合图9(a)及表6可知,2018年1—5月扎龙湿地水域面积呈“W”型波动起伏变化,大致在1 300~1 500 km²范围间浮动,平均面积约为1 400 km²。7月则呈现出2018年年内最低谷,降至756.98 km²,与5月相比,湿地水域面积发生巨大变化,减少了约660 km²。8—10月水域面积持续缓慢回升,其中8月份回升幅度最大,升至998.82 km², 11月份面积降至1 027.91 km²,12月份面积急速上升至1 649.70 km²,为年内水体面积最高峰,与面积最少的7月份相差了892.72 km²,增加了117.93%。

2019年1月,扎龙湿地水域面积为1 335.07 km²,3月份面积呈现上升趋势,升至1 561.70 km²。4—7月份总体面积呈下降趋势,4月份降至1 325.67 km²,7月份降至2019年年内最低,仅有871.37 km²。8—11月扎龙湿地水域面积则呈现出高速状态,8月份湿地水域面积缓慢回升至990.74 km²,10月份湿地水域面积呈现急速上升状态,升至1 649.70 km²,为年内水域面积最高峰,同时年内水域面积最高峰的10月份与最低谷的7月份相比,面积增加了747.63 km²,增加了85.80%。11月份与10月相比较为稳定,仅呈现出略微下降趋势,降至1 575.05 km²,12月份持续下降至1 409.57 km²。

结合图9(b)及表6可知,2018年湿地水域面积与2019年相比,1月份2018年多出了154.97 km²,3月份中2019年多出了108.97 km²,4月份则表现出相对持平状态。7月份2019年多出了114.39 km²,8月份两年年表现出相对持平状态,2018年仅比2019年多出8.08 km²。10月份呈现较大的反转,2019年比2018年湿地水域面积多出了482.11 km²。11月份出现了2 a之间差别最大的峰值,2019年比多出547.14 km²。12月份产生反转,由于2018年12月份水域面积增加较大而2019年12月份水域面积又略有减少,导致2018年比2019年多出了240.13 km²。

总体上看,2019年扎龙湿地水量比2018年更为充足,计算有效数据,2019年扎龙湿地水域面积比2018年多出约863.32 km²,根据调查,了解到2019年秋季东北的部分地区连续受暴雨侵袭而致洪水灾害,故2019年的秋季水量异乎2018年水量。

将通过加权σ⁰和μ⁰影像数据得出新的加权影像数据使用阈值分割法提取水体,利用2018年3月、8月、10月、12月及2019年同期的影像提取的湿地水域面积如表7所示,变化趋势如图10所示。

由图10及表7可知,2018年3月扎龙湿地水域面积为1 461.76 km²。8月份面积下降至1 052.46 km²,与3月份对比差值巨大,达409.3 km²。到10月份面积升至1 167.33 km²,12月升到了高峰,达1 628.46 km²。2019年3月扎龙湿地水域面积为1 566.53 km²,8月份降至996.37 km²,10月升至1 613.01 km²,与8月差值616.64 km²,12月份降至1 434.72 km²。

整体上看,扎龙湿地在10月份的2个时期水域面积表现出了较为明显的差异,2019年比2018年多出了约为445 km²。其次变化较大的为12月,自10月水域面积骤增后2019年12月呈现出来了较大的降幅,降至1 434.72 km²。

选中2018年3月、8月、10月、12月及2019年同期的影像,通过随机森林算法对多源遥感数据集提取的湿地水域面积如表8所示。

由表8可知,2018年3月扎龙湿地水域面积为1 461.64 km²,2018年8月水域降至954.93 km²,落差达506.71 km²。2018年10月扎龙湿地水域面积稍有回升,升至1 187.23 km²,2018年12月升到了高峰,比2018年春季的3月仍多出136.79 km²,达1 598.43 km²,且此时与较为干旱的8月相比多出了643.5 km²。2019年3月扎龙湿地水域面积为1 575.27 km²,2019年8月水域面积又降到低谷,降至979.91 km²,下降了37.79%。2019年10月,水域面积回升至丰水期,升至1 610.59 km²,12月不升反降,面积降至1 407.98 km²。2 a皆呈现出的丰水期-枯水期-丰水期的变化趋势,其中10月份差异最大,2019年10月比2018年10月多出了约425 km²。

本论文利用Sentinel-1A的σ⁰和μ⁰影像,通过阈值分割法和随机森林算法提取扎龙湿地2018—2019年逐月的水体信息,并对结果进行年内年际变化定量分析,获取具体研究结论如下:

1)通过使用Sentinel-2A的AWEI指数数据的辅助验证,得出结论: RF算法的分类精度最高,平均差值绝对值为6.69km²; 其次为基于加权σ⁰与μ⁰影像使用阈值分割方法的分类精度,平均差值绝对值为11.2 km²; 再是基于σ⁰影像使用阈值分割方法的分类精度,平均差值绝对值为12.87 km²; 最后是基于μ⁰影像使用阈值分割方法的分类精度,该方法精度最低,平均差值绝对值为13.07 km²。

2)扎龙湿地年内变化: 水域面积变化规律和季节有很大关系,春季水域面积在1 300~1 600 km²浮动,夏季在700~1 000 km²浮动,秋季在900~1 200 km²浮动,冬季在1 300~1 600 km²浮动。

3)扎龙湿地年际变化: 总体上看,2019年扎龙湿地水量比2018年更充足,因2019年秋季受暴雨侵袭,故10月和11月份扎龙湿地水域面积比2018年多出1 044.25 km²。扎龙湿地水体面积每年皆整体上呈现丰水-枯水-丰水的变化规律,其中不同年份又各自具有独特的水体积变化现象。