基于多時相Sentinel-1A的沼澤濕地水面時空動態變化監測

韋嫦， 付波霖， 覃嬌玲， 王雅南， 陳智瀚， 劉兵

(桂林理工大學測繪地理信息學院，桂林 541000)

0 引言

濕地是位于陸生生態系統和水生生態系統之間的至關重要的過渡地帶，它有著貯蓄水量、調節氣候、保護環境和保護生物多樣性等重要作用，是世界上最有價值的自然資源之一[1-2]。水是濕地三要素之一，是維持濕地的主要因素[3]。由于人類活動過程中長期忽視對地球環境的保護，地球生態環境逐漸惡化，水循環系統遭到破壞，全球濕地萎縮和功能退化成為突出問題[4]。因此，為了防止濕地進一步流失，更好地管理和保護現有的濕地資源，有必要研究濕地水域變化，監測其年內、年際變化。

目前國內外濕地水體遙感提取方式主要有以下3種： ①基于傳統光學遙感圖像的遙感技術，通過構造波段間的光譜關系提取水體[5]。如甄佳寧[6]使用Landsat系列衛星影像對長春濕地采取不同的水體指數和閾值分割法提取水體信息； Jin等[7]利用Landsat時間序列影像監測德爾馬瓦半島的濕地水位淹沒動態，增進對濕地的了解并確保生態系統服務的可持續性。不足之處是光學遙感圖像的遙感技術受云雨天氣等條件限制，對監測濕地水域年內變化影響頗大[8]。②基于微波遙感圖像的遙感技術，通過合成孔徑雷達(synthetic aperture Radar，SAR)后向散射系數(σ0)、干涉相干系數(μ0)，構建SAR影像數據集提取水體信息。徐怡波[9]基于多時相ENVISAT衛星影像分析水體的后向散射系數特征，采用閾值分割的方法提取水體，實現對水面變化信息的監測； 付波霖等[10]利用多頻率、多時相SAR影像，采用合成孔徑雷達差分干涉測量技術(interferometric synthetic aperture Radar，InSAR)，構建沼澤水位變化監測模型，探究了不同時相沼澤濕地干涉相干性的差異。不足之處是使用單一的σ0或單一的μ0方式有一定的局限性。③基于微波和光學遙感圖像的融合，利用2種傳感器的互補性提取水體信息[11]。關韻桐[12]利用多時相Sentinel-1 雷達影像為主數據源，以Sentinel-2和Landsat8 OLI 光學影像為輔數據源，將光學影像計算出的歸一化植被指數(normalized difference vegetation index, NDVI)分類，根據雷達影像數據分析了σ0與雷達參數之間的關系，利用反演模型對大山包高原濕地進行土壤水分反演，獲取了土壤的水分分布狀況，根據實測數據驗證了反演結果的準確性； Mizuochia等[13]使用AMSR系列、MODIS和Landsat衛星數據，開發并評估了基于查找表的時空數據融合方法，并將其應用于半納米比亞的干旱季節性濕地，結果表明時間序列圖能準確地描述濕地的季節變化和年際變化。綜上，目前缺乏整合長時間序列的微波遙感影像的后向散射系數影像和干涉相干系數影像方法，構建多源遙感數據集來檢測濕地水域面積動態變化方面的研究。

本文以東北扎龍國家級自然保護區作為研究區，以2018—2019年Sentinel-1A衛星影像數據作為主數據源，輔以Sentinel-2A衛星影像數據，提取保護區年內和年際的SAR后向散射系數和干涉相干系數，構建濕地水面積提取SAR影像數據集，利用機器學習算法構建沼澤濕地水面動態變化監測模型，實現扎龍濕地水域面積的動態變化監測，探究濕地水域年內和年際變化規律。

1 研究區概況和數據源

1.1 研究區概況

扎龍濕地(圖1)位于黑龍江省松嫩平原西部的烏裕爾河下游，地理坐標為E123°47′～124°37′，N46°52′～47°32′，邊界平面呈現不規則橄欖形，濕地南北約80.6 km，東西約58.0 km，總面積約為2 100 km2，是面積排名亞洲第一、世界第四的國家濕地生態系統保護區，也是世界上面積最大的蘆葦濕地,它于1987年內被國務院批準為國家級自然保護區[14-15]。多年來扎龍濕地面臨水資源短缺、水位下降、資源枯竭、面積分解等問題，嚴重影響了濕地的生態環境健康，甚至影響了當地經濟和社會的可持續發展。

圖1 研究區地理位置

1.2 數據源及數據預處理

本論文使用的數據主要為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所示。

表1 影像數據具體信息

2 研究方法

本文主要利用SNAP軟件對多時序Sentinel-1衛星數據進行一系列預處理，主要包括熱噪聲去除、輻射定標(將數字像素值轉換為后向散射系數)、條帶修復(去除條帶噪聲)、濾波處理(Define Lee濾波)、多視處理(距離向視數設為4)、斜距轉地距及地理編碼獲取UTM WGS 1984投影坐標系的σ0數據，如圖2(a)所示； 經過應用軌道文件(獲取準確的衛星位置和速度信息)、影像配準(將主影像和輔影像進行配準獲取2個時相影像相融合的新影像)、增強光譜多樣性(校正條帶的距離和方位向)、干涉條紋圖的相干系數提取(對原始圖像參數賦予相干系數)、條帶修復、多視處理(距離向視數設為4)、地形校正和影像裁剪獲取μ0數據，如圖2(b)所示。本文選擇每個月靠近月中的VV極化Sentinel-

(a) σ0(b) μ0

1A/B影像提取研究區σ0數據，選擇每個月2景VV極化的Sentinel-1A/B影像計算逐月的μ0數據。

本文總體技術方案如圖3所示。第一部分為預處理模塊，對原始Sentinel-1A衛星影像數據經過熱噪聲去除、輻射定標等預處理過程，獲取研究區σ0影像； 經過應用軌道文件、影像配準等預處理過程獲取研究區μ0影像； 經過波段運算方法加權研究區σ0影像和μ0影像獲取新的加權影像。第二部分為數據處理模塊，將第一部分預處理模塊得到的3種類型影像對照彩色光學影像找到水體大致位置，使用閾值分割法確定最佳閾值分割參數提取濕地水體，逐月計算提取的濕地水體面積，獲取水位動態信息。第三部分為數據分析模塊，對所得影像數據基礎建立濕地水面積提取SAR影像數據集，使用隨機森林算法對多源遙感數據集進行濕地水體和非水體分類，計算水體面積，獲取水位動態信息。第四部分為精度驗證模塊，對所的結果進行定量分析，探究濕地水域年內年際變化規律，驗證精度。

2.1 基于σ0的濕地水體提取

光在水體中傳播時會發生散射，入射能量與表面相互作用后，再輻射向各個方向，其中散射光沿入射方向相反方向的散射的過程稱為后向散射，它也是雷達接收機接收到的散射，通常以σ0表示[16]。

σ0是水體的固有光學量，水體在微波范圍內的σ0較低，雷達波的吸收相對完整，在雷達圖像上呈現較暗的黑色區域，水體和非水體邊界明顯，可通過閾值分割等多種方法輕易將水體和其他地物區分開。其中，σ0是體散射函數β(θ)的散射角θ對后向半球的積分，在數學上定義表達式為：

(1)

式中λ(θ)為σ0與體散射函數關系的常數，該常數用于此研究領域近似恒定。

本文在ENVI中實現閾值分割濕地水體提取。閾值分割是一種操作簡易且效率高的圖像分割方法，該方法使用一個或多個閾值來分割圖像的灰度級，并且在相同閾值范圍內的灰度值將被識別為相同類型的地物。本文將σ0設為F(x,y)，以某個灰度范圍[k1,k2]作為閾值范圍進行分割，分割后生成的圖像為f(x,y)，在閾值范圍內的像素全部更改為1，其他像素均更改為0，即劃分灰度級后的圖像為二值圖，在數學上定義為：

。

(2)

2.2 基于μ0的濕地水體提取

μ0是描述2幅單視復數影像相干程度高低的參數，即地物2次拍攝的差異，該系數可以表明某些理化特征，因此可以用作輔助地物判別的獨立參數。由于水體的μ0值為0，故在圖上水體呈現暗黑色，極易于區分水體和其他地物。對于單視復數影像，使用SNAP的相干變化檢測處理方法來生成地理編碼的相干系數圖，干涉測量的相干系數介于0～1之間(0表示不相干，1表示完全相干)，假定s1和s2分別為同一分辨率的一定時間間隔內獲取的2幅影像，則μ0在數學上定義表達式為：

，

(3)

式中：μ0為相干系數；s*為s的共軛復數；s1和s2為經過預處理的相干系數影像。

本文在ENVI軟件中對μ0數據進行閾值分割，將圖像二值化，輸出水體的柵格文件，計算水體面積。由于濕地水域的相干系數值會明顯低于周圍地物，水體在μ0數據上表現出的灰度值明顯偏低，因此很好選定最佳閾值參數對水體進行提取。

2.3 基于加權σ0和μ0的濕地水體提取

單一的σ0和μ0數據提取水體都分別具有各自優勢，同時也存在各自的局限性。本文利用波段運算將二者分別賦予權重疊加融合，便可取得繼承二者特性的加權σ0和μ0數據的新影像數據。權重的賦予是根據奧維地圖2018年、2019年的高空間分辨率歷史光學遙感影像來驗證單一的σ0和單一的μ0提取水體精度進行權重劃分，如σ0影像提取水體效果更佳，則權重愈大，反之愈小。

目前，加權σ0和μ0數據可以通過ENVI進行波段之間的運算實現。波段運算是遙感信息提取的重要方法之一，是基于研究對象的多光譜波段反射特性，選擇合適的波段進行運算，目的是增加對象與背景的差異[17]。

本文在ENVI軟件中對新影像采用閾值分割法提取濕地水體，通過計算像元數和像元大小的關系即可獲知濕地水域面積，如表2所示。

表2 2018年加權σ0和μ0數據波段運算公式及閾值參數

2.4 基于隨機森林算法的濕地水體提取

隨機森林是一種新穎的、極其靈活的機器學習算法，通過集成學習的思想將多棵決策樹集成，其中每棵決策樹都是一個分類器，即對于一個輸入樣本，N棵樹都將具有N個分類結果。隨機森林將集成的全部分類結果進行隨機而有放回的抽樣，通過投票機制，最終輸出投票最多的類別，該分類方法在目前所有分類的算法中，具有非常高的準確率[18]。

本文將經過處理的Sentinel-1A構建濕地水面積提取SAR影像數據集，其數據構成如表3所示，波段組成示例圖如圖4所示。每個數據集由5個波段組成，對其進行濕地水體提取，采用ENVI的拓展模塊隨機森林分類器繪制訓練樣本和預測變量的隨機子集，由于決策樹數目過少會導致分類誤差大，過大則會導致隨機森林構建時間過長，故本文生成1 000棵決策樹形成整體分類，表3中展示2019年10月SAR影像數據集對照奧維地圖彩色光學影像選取的地表真實感興趣區域(region of interest，ROI)樣本信息，選取水體和非水體樣本數各100個(70%為訓練樣本，30%為驗證樣本)。

表3 2019年10月多源數據集參數及選取ROI樣本信息

圖4 扎龍濕地彩色影像示例

2.5 基于Sentinel-2的濕地水體提取

自動化水體提取指數(automated water extraction index，AWEI)是一種可以提高包括陰影和暗表面在內的區域的分類精度的指數[19]。由于AWEI方程通過應用強制非水像素低于0和水像素高于0的系數來提高水和非水像素的可分離性，可以采用默認的0閾值來自動客觀地提取水體，而不采用人為選取的主觀最佳閾值[20]。表達式為：

AWEInsh=4(ρG-ρM)-
(0.25ρN+2.75ρS)

，

(4)

AWEIsh=ρB+0.25ρG-
1.5(ρN+ρM)-0.25ρS

，

(5)

式中：ρG為綠光波段；ρB為藍光波段；ρN為近紅外波段；ρM為中紅外波段；ρS為熱紅外波段；AWEInsh為沒有陰影下提取的自動化水體提取指數；AWEIsh為陰影下提取的自動化水體提取指數。

本文在ENVI中通過波段運算計算光學影像數據獲取AWEI指數影像，將其作為本論文濕地水體提取基準，通過閾值分割獲取水體信息，分析驗證前文4種基于不同處理方法得到的影像數據的水體提取精度。

3 扎龍濕地水域面積時空動態變化定量分析

3.1 基于Sentinel-2的水域面積驗證

本文選取了2019年的3月、8月、10月和12月的Sentinel-2A衛星的影像進行水體提取方法對比實驗，通過計算Sentine-2A影像數據獲取4個月份的AWEI指數，以研究區AWEI數據作為本文濕地水體提取基準，分析驗證前文4種基于不同處理方法得到的影像數據的水體提取精度。結果圖如圖5所示。

(a) 基于AWEI數據(閾值分割法)(b) 基于σ0數據(閾值分割法)(c) 基于μ0數據(閾值分割法)(d) 基于加權σ0與μ0數據(閾值分割法)(e) 基于多源SAR影像數據集(RF分類法)

表4為2019年4期不同方法水體提取面積對比?？梢缘弥S機森林算法的分類精度最高； 其次為基于加權σ0與μ0影像使用閾值分割方法的分類精度； 再是基于σ0影像使用閾值分割方法的分類精度； 最后是基于μ0影像使用閾值分割方法的分類精度。

表4 2019年基于不同處理方法得到的數據的水體提取結果比較

3.2 基于σ0的扎龍濕地水域面積時空動態變化定量分析

通過閾值分割法對σ0數據進行水體提取，閾值設定范圍在0～0.01。利用2018—2019年逐月的Sentinel-1A影像提取的扎龍濕地水域面積如表4所示，變化趨勢如圖6所示。

(a) 水域面積年內變化曲線 (b) 水域面積年際變化對比柱狀圖

3.2.1 扎龍濕地年內水域面積動態變化分析

結合圖6(a)及表5可知，2018年1月，扎龍濕地水域面積為1 503.50 km2，2018年2—5月，濕地面積沒有明顯的變化，基本在1 400 km2左右，呈輕微波動性狀態。從6月份開始，濕地水域面積發生巨大變化，6月份面積陡降至年內最低谷770.51 km2，較上個月下降約640 km2，下降了45.39%。7—11月份整體上稍有回升，平均面積約為1 027 km2，其中，7月份水域面積增幅較小，僅增至839.16 km2。8月份增幅較大，增至1 098.05 km2，9月份濕地水域面積呈現出略微的下降，降至1 059.18 km2，10月份水域面積有些許回升，升至1 177.69 km2，11月份又降至964.1 km2，12月份水域面積快速回升至1 609.32 km2，這個階段為2018年濕地水域面積增長速度最快的時期，達到年內最高峰，相比于年內最低的6月份，增加了108.86%。

表5 2018—2019年基于σ0數據提取的扎龍濕地水域面積(閾值分割法)

2019年1月，扎龍濕地水域面積為1 361.28 km2，2月份面積相對穩定，有少許回升達1 401.55 km2，3月份持續上升至1 598.70 km2。4—7月水域面積急劇下降，7月份降至年內最低谷，僅有768.11 km2，與年內最高峰的3月份相比，7月份下降830.59 km2，下降了51.95%，如圖7(a)所示。8—12月整體呈高速度增長趨勢，10月之時水域面積已達1 590.14 km2，11月份與10月相比較為穩定，，12月份呈現出略微下降趨勢，降至1 541.74 km2，如圖7(b)所示。

(a) 2019年7月(b) 2019年12月

3.2.2 扎龍濕地年際水域面積動態變化分析

結合圖6(b)及表5可知，2018年濕地水域面積與2019年相比，1月份2018年多出了142.22 km2，而與之相反，2月份2019年多出了89.76 km2，同樣的3月份中2019年的明顯較多，多出了155.92 km2，4月份則表現出相對持平狀態，僅差4.48 km2。5月份、7月份和8月份2018年比2019年分別多出244.59 km2，71.05 km2和103.9 km2，10月份呈現反轉，2019年多出了412.45 km2。11月份出現了年際差別最大的峰值，2019年水域面積多出了631.8 km2，如圖8所示。最后12月份時2018年比2019年稍多177.23 km2。

(a) 2018年11月(b) 2019年11月

據有效數據，2018年扎龍濕地月平均水體面積為1 269.03 km2，2019年扎龍濕地月平均水體面積為1 334.65 km2。總體上看，2019年扎龍濕地水量比2018年更為充足?？梢钥闯觯m然扎龍濕地水域面積有明顯的季節性變化，但每年都有其變化的特征，會隨著降雨量、溫度及人為補水等多種因素而出現獨特的現象。

3.3 基于μ0的扎龍濕地水域面積時空動態變化定量分析

通過閾值分割法對μ0數據進行水體提取，閾值設定范圍在0～1。獲取2018—2019年扎龍濕地逐月水域面積如表6所示，變化趨勢如圖9所示。

表6 2018—2019年基于μ0數據提取的扎龍濕地水體面積(閾值分割法)

(a) 水域面積年內變化曲線 (b) 水域面積年際變化對比柱狀圖

3.3.1 扎龍濕地年內水域面積動態變化分析

結合圖9(a)及表6可知，2018年1—5月扎龍濕地水域面積呈“W”型波動起伏變化，大致在1 300～1 500 km2范圍間浮動，平均面積約為1 400 km2。7月則呈現出2018年年內最低谷，降至756.98 km2，與5月相比，濕地水域面積發生巨大變化，減少了約660 km2。8—10月水域面積持續緩慢回升，其中8月份回升幅度最大，升至998.82 km2， 11月份面積降至1 027.91 km2，12月份面積急速上升至1 649.70 km2，為年內水體面積最高峰，與面積最少的7月份相差了892.72 km2，增加了117.93%。

2019年1月，扎龍濕地水域面積為1 335.07 km2，3月份面積呈現上升趨勢，升至1 561.70 km2。4—7月份總體面積呈下降趨勢，4月份降至1 325.67 km2，7月份降至2019年年內最低，僅有871.37 km2。8—11月扎龍濕地水域面積則呈現出高速狀態，8月份濕地水域面積緩慢回升至990.74 km2，10月份濕地水域面積呈現急速上升狀態，升至1 649.70 km2，為年內水域面積最高峰，同時年內水域面積最高峰的10月份與最低谷的7月份相比，面積增加了747.63 km2，增加了85.80%。11月份與10月相比較為穩定，僅呈現出略微下降趨勢，降至1 575.05 km2，12月份持續下降至1 409.57 km2。

3.3.2 扎龍濕地年際水域面積動態變化分析

結合圖9(b)及表6可知，2018年濕地水域面積與2019年相比，1月份2018年多出了154.97 km2，3月份中2019年多出了108.97 km2，4月份則表現出相對持平狀態。7月份2019年多出了114.39 km2，8月份兩年年表現出相對持平狀態，2018年僅比2019年多出8.08 km2。10月份呈現較大的反轉，2019年比2018年濕地水域面積多出了482.11 km2。11月份出現了2 a之間差別最大的峰值，2019年比多出547.14 km2。12月份產生反轉，由于2018年12月份水域面積增加較大而2019年12月份水域面積又略有減少，導致2018年比2019年多出了240.13 km2。

總體上看，2019年扎龍濕地水量比2018年更為充足，計算有效數據，2019年扎龍濕地水域面積比2018年多出約863.32 km2，根據調查，了解到2019年秋季東北的部分地區連續受暴雨侵襲而致洪水災害，故2019年的秋季水量異乎2018年水量。

3.4 基于加權σ0和μ0的扎龍濕地水域面積時空動態變化定量分析

將通過加權σ0和μ0影像數據得出新的加權影像數據使用閾值分割法提取水體，利用2018年3月、8月、10月、12月及2019年同期的影像提取的濕地水域面積如表7所示，變化趨勢如圖10所示。

由圖10及表7可知，2018年3月扎龍濕地水域面積為1 461.76 km2。8月份面積下降至1 052.46 km2，與3月份對比差值巨大，達409.3 km2。到10月份面積升至1 167.33 km2，12月升到了高峰，達1 628.46 km2。2019年3月扎龍濕地水域面積為1 566.53 km2，8月份降至996.37 km2，10月升至1 613.01 km2，與8月差值616.64 km2，12月份降至1 434.72 km2。

表7 基于加權σ0與μ0數據提取的水體面積(閾值分割法)

圖10 扎龍濕地基于加權σ0與μ0數據提取的年際水體面積折線圖(閾值分割法)

整體上看，扎龍濕地在10月份的2個時期水域面積表現出了較為明顯的差異，2019年比2018年多出了約為445 km2。其次變化較大的為12月，自10月水域面積驟增后2019年12月呈現出來了較大的降幅，降至1 434.72 km2。

3.5 基于隨機森林算法的扎龍濕地水域面積時空動態變化定量分析

選中2018年3月、8月、10月、12月及2019年同期的影像，通過隨機森林算法對多源遙感數據集提取的濕地水域面積如表8所示。

表8 扎龍濕地基于隨機森林算法提取水體面積

由表8可知，2018年3月扎龍濕地水域面積為1 461.64 km2，2018年8月水域降至954.93 km2，落差達506.71 km2。2018年10月扎龍濕地水域面積稍有回升，升至1 187.23 km2，2018年12月升到了高峰，比2018年春季的3月仍多出136.79 km2，達1 598.43 km2，且此時與較為干旱的8月相比多出了643.5 km2。2019年3月扎龍濕地水域面積為1 575.27 km2，2019年8月水域面積又降到低谷，降至979.91 km2，下降了37.79%。2019年10月，水域面積回升至豐水期，升至1 610.59 km2，12月不升反降，面積降至1 407.98 km2。2 a皆呈現出的豐水期-枯水期-豐水期的變化趨勢，其中10月份差異最大，2019年10月比2018年10月多出了約425 km2。

4 結語

本論文利用Sentinel-1A的σ0和μ0影像，通過閾值分割法和隨機森林算法提取扎龍濕地2018—2019年逐月的水體信息，并對結果進行年內年際變化定量分析，獲取具體研究結論如下：

1)通過使用Sentinel-2A的AWEI指數數據的輔助驗證，得出結論： RF算法的分類精度最高，平均差值絕對值為6.69km2； 其次為基于加權σ0與μ0影像使用閾值分割方法的分類精度，平均差值絕對值為11.2 km2； 再是基于σ0影像使用閾值分割方法的分類精度，平均差值絕對值為12.87 km2； 最后是基于μ0影像使用閾值分割方法的分類精度，該方法精度最低，平均差值絕對值為13.07 km2。

2)扎龍濕地年內變化： 水域面積變化規律和季節有很大關系，春季水域面積在1 300～1 600 km2浮動，夏季在700～1 000 km2浮動，秋季在900～1 200 km2浮動，冬季在1 300～1 600 km2浮動。

3)扎龍濕地年際變化： 總體上看，2019年扎龍濕地水量比2018年更充足，因2019年秋季受暴雨侵襲，故10月和11月份扎龍濕地水域面積比2018年多出1 044.25 km2。扎龍濕地水體面積每年皆整體上呈現豐水-枯水-豐水的變化規律，其中不同年份又各自具有獨特的水體積變化現象。