格陵蘭冰蓋次表面湖多源遙感監測
——以格陵蘭中西部流域為例

李嵐靜， 陳卓奇,2,3*， 鄭雷,2,3， 程曉,2,3

1 中山大學測繪科學與技術學院, 廣東珠海 519082 2 南方海洋科學與工程廣東省實驗室(珠海), 廣東珠海 519082 3 中國高校極地聯合研究中心, 北京 100875

0 引言

格陵蘭冰蓋是全球第二大陸地冰蓋，其完全消融將引起全球海平面上升7.4 m(Morlighem et al., 2017).過去幾十年，北極是全球氣溫上升最快的地區，增溫幅度達到1.2 ℃/10a，是全球平均增溫速度的2倍以上.快速上升的氣溫增強了格陵蘭冰蓋表面消融，促進了冰面湖的形成，加強了冰蓋表面徑流，加劇了冰蓋的物質損失(Leeson et al., 2015; Chen et al., 2006; Trusel et al., 2018).研究表明，格陵蘭冰蓋物質損失平均速率已從1991—2001年間的34(-6～74)Gt/a大幅增至2002—2011年間的215(157～274)Gt/a，物質損耗的加速度為21.9±1 Gt/a(Stocker, 2014)，其中50.3%的物質損失由冰蓋表面消融貢獻(The IMBIE Team, 2020).因此，格陵蘭冰蓋表面消融是格陵蘭冰蓋物質損失最重要的影響因素之一.

格陵蘭冰蓋表面消融加劇直接導致了格陵蘭冰面湖的增加(Howat et al., 2013).格陵蘭冰面湖數量、面積和深度提取、形成機制、生長周期和排水事件監測成為了國際研究熱點(Pope et al., 2016; Williamson et al., 2018; Yang et al., 2019).但近期研究表明格陵蘭大量冰面湖在冬季不會完全凍結.Koenig等(2014, 2015)基于現場觀測首次發現了格陵蘭冰蓋次表面湖的存在，并利用IceBridge資料證實了格陵蘭冰蓋存在大量的次表面湖.這些次表面湖在冬季保持液態，溫度在0 ℃左右，比冰蓋表面溫度高，成為冰蓋內部的一個熱庫(Koenig et al., 2015).滯留在冰蓋次表層的融水對格陵蘭冰蓋消融、穩定性和物質損失會產生重要的影響(Harper et al., 2012; Machguth et al., 2016; MacFerrin et al., 2019)，因此，對格陵蘭次表面湖的研究對于理解格陵蘭冰蓋物質平衡具有非常重要的價值.

與格陵蘭冰面湖不同，格陵蘭冰蓋次表面湖被掩埋在冰雪層的內部，難以用可見光遙感手段進行監測.Koenig等(2015)利用IceBridge計劃獲得的機載雪雷達數據(2～6.5 GHz)，提取了2009—2012年格陵蘭冰蓋的次表面湖分布.雪雷達具有很好的穿透性，可以準確的提取次表面湖的位置，但不能識別次表面湖的形狀和面積等信息，且僅能獲得飛機航線經過區域的次表面湖，難以獲得全格陵蘭冰蓋次表面湖的分布.

另一種可用于監測格陵蘭冰蓋次表面湖的方法是合成孔徑雷達技術(Synthetic Aperture Radar，SAR).有研究表明C波段SAR信號可以穿透雪層1～10 m(Rignot et al., 2001).由于冰層和水體的后向散射系數顯著不同(Fahnestock et al., 1993).科學家提出可以利用Sentinel-1衛星C波段后向散射系數監測格陵蘭冰蓋次表面湖的分布(Miles et al., 2017; Schr?der et al., 2020; Dunmire et al., 2021).這類算法的理論基礎是HH和HV極化的后向散射系數分別代表了次表面湖的面散射和冰體的體散射信號，通過對HH和HV極化信號的差異進行分離實現對格陵蘭冰蓋次表面湖的監測(Schr?der et al., 2020).

基于SAR數據的格陵蘭冰蓋次表面湖監測技術可以實現對次表面湖位置、形狀、面積的監測.但Sentinel-1衛星C波段后向散射系數易受到地面幾何屬性的影響(如表面粗糙度、雪粒徑、冰層內部結構等)(Fahnestock et al.,1993).這些問題給格陵蘭冰蓋次表面湖監測帶來了很大的困難(Schr?der et al., 2020).現有研究大多采用監督式學習的方法提取格陵蘭冰蓋次表面湖(Schr?der et al., 2020; Dunmire et al., 2021).該類方法的優勢是算法簡單，但其算法精度嚴重依賴于樣本庫的選擇且常存在眾多人為設定的閾值.但受制于HH、HV后向散射系數的特點和監督式學習方法的特性，利用監督式學習實現格陵蘭冰蓋全域次表面湖的監測仍存在很大的不確定性.本研究的目的是利用無監督式學習的方法，以格陵蘭中西部流域為例，提取夏季出現過但冬季未完全凍結的水體.

基于此，本文旨在利用Landsat8和Sentinel-1的影像，在提取出夏季冰面湖范圍的基礎上，發展基于Rosin圖像分割算法的格陵蘭冰蓋次表面湖提取方法，對方法有效性進行驗證后，應用該方法提取2018—2019年冬季格陵蘭中西部流域次表面湖.

1 研究區域

參考Rignot和Mouginot(2012)的格陵蘭流域劃分結果，本研究選取格陵蘭中西部流域為研究區(70°N—75°N, 40°W—50°W)(圖1).此流域平均高程為2427.1 m，存在廣泛而活躍的冰川水文系統(Miles et al., 2017)，并且有著大量的沿岸冰面湖(Lewis and Smith, 2009).

圖1 研究區范圍采用基于Landsat-8的格陵蘭真彩色合成影像(Chen et al., 2020).格陵蘭流域的劃分采用Rignot和Mouginot(2012)流域分割方法.Fig.1 Study areaGreenland true color synthetic image based on Landsat-8 is obtained from Chen et al., (2020). Greenland watershed is obtained from Rignot and Mouginot (2012).

2 數據與研究方法

2.1 Sentinel-1

Sentinel-1由處于同一軌道平面的兩顆攜帶C波段合成孔徑雷達極軌衛星Sentinel-1A和Sentinel-1B組成，Sentinel-1A和Sentinel-1B分別于2014年4月3日和2016年4月15日發射升空.單顆衛星重訪周期為12天，雙星重訪周期提高至6天，更有利于對研究區的連續觀測(Schr?der et al., 2020)，在Sentinel-1B發射前Sentinel-1未實現格陵蘭地區影像的全覆蓋.

作為主動微波遙感衛星，Sentinel-1在對格陵蘭的觀測中有其獨特的優勢.其一，由于其接收的是傳感器發射后經地物反射后的后向散射信息，不受黑夜、云的影響，故可以全天時、全天候對地實施觀測；其二，與光學影像相比，Sentinel-1搭載的C波段SAR具有一定的地表穿透能力，可以穿透覆蓋在次表面湖之上的冰或雪，觀測到次表面的水體.

Sentinel-1 SAR包括四種工作模式，分別為干涉寬帶模式(Interferometric Wide swath Mode，IW)、條帶模式(Strip Map Mode，SM)、超幅寬模式(Extra-Wide Swath Mode，EWS)以及波模式(Wave Mode)，IW模式以5 m×20 m的空間分辨率獲取250 km的帶狀影像，是Sentinel-1在陸地上的主要采集模式.EWS模式以犧牲空間分辨率為代價提供了非常大的區域覆蓋，由于次表面湖面積較小，故在此選擇IW模式，在此模式下，衛星具有雙極化的能力(HH+HV或VV+VH)，次表面湖在不同極化方式中的表現形式也不盡相同.同極化信號對表散射比較敏感，而交叉極化對體散射敏感，可視為是體散射的表達 (Schr?der et al., 2020)，且交叉極化信噪比較大，受表面粗糙度的影響較小(Miles et al., 2017)，反映了由于液態水的存在而導致的微波后向散射的變化(Dunmire et al., 2021)，因此本研究采用交叉極化數據開展工作.

歐洲航天局(European Space Agency，ESA)為Sentinel-1提供4種數據產品，本文采用一級GRD(Ground Range Detected)數據產品.

2.2 夏季冰面湖提取方法

已有研究表明，部分格陵蘭冰蓋次表面湖是由夏季冰面湖發育而來.以夏季冰面湖分布為先驗知識，可以大幅提高冬季次表面湖提取的精度和效率.基于此，本研究開展了格陵蘭冰蓋夏季冰面湖分布提取研究，采用2018年5月至9月Landsat8影像共108幅.

在格陵蘭冰面湖提取過程中，夏季冰面湖和融化的冰、積雪、粒雪的近紅外的反射率都很低，無法區分，故采用改進的歸一化水體指數NDWIice，利用藍光波段和紅光波段的歸一化比值，以識別這種位于格陵蘭冰面上獨特環境的水體特征(Yang and Smith, 2013).

(1)

式(1)中的BLUE和RED對應Landsat8的第二和第四波段.

本文取0.5作為NDWIice指數的全局經驗閾值,將整個研究區二值化為水區域和非水區域，認為小于0.5為非水區域(Dunmire et al., 2021).

加入歸一化積雪指數NDSI,利用綠光波段和短波紅外波段的歸一化比值可以消除云和裸巖造成的冰面湖提取誤差(Moussavi et al., 2020).

(2)

式(2)中的GREEN和SWIR對應Landsat8的第三波段和第六波段.

本文取0.85作為NDSI的經驗閾值(Dunmire et al., 2021),去除NDSI小于0.85范圍.

2.3 次表面湖探測方法

由于格陵蘭冰蓋次表面湖被掩埋在冰雪層的內部，可見光遙感手段僅能監測格陵蘭冰面湖的變化，而無法對次表面的情況進行分析，在冬季無法獲得次表面有效信息(圖2).在冬季，冰層和水體由于介電常數不同，兩者的后向散射系數顯著不同(圖3)，以此為依據可利用SAR影像進行次表面湖的提取.

圖2 冰面湖探測實例Fig.2 Example of the detection of supraglacial lakes

圖3 格陵蘭中西部流域水體、冰面直方圖Fig.3 Histogram of water and ice in the Central West basin of Greenland

本研究提取次表面湖的流程如圖4所示.首先，進行夏季冰面湖提取.對2018年5月至9月Landsat8地表反射率影像數據集進行云掩模，利用NDWIice和NDSI提取每幅影像冰面湖區域.將每幅影像的提取結果進行最大值合成，得到該年格陵蘭中西部流域的最大水體范圍，作為冬季可能存在次表面湖的區域，對最大水體范圍中每個冰面湖轉為矢量數據并進行編碼，考慮到冰流速的逐年增加(Schr?der et al., 2020)，對矢量數據進行形態學擴張以確保可以完全提取冬季次表面湖范圍.

圖4 次表面湖提取流程圖(a) 夏季冰面湖可見光影像； (b) 夏季冰面湖范圍； (c) 冬季次表面湖可見光影像； (d) 冬季次表面湖SAR影像； (e) Rosin閾值分割方法； (f) 冬季次表面湖提取結果.Fig.4 Flow chart of the detection of buried lakes(a) Visible image of a supraglacial lake in summer; (b) Range of a supraglacial lake in summer; (c) Visible image of a buried lake in winter; (d) SAR image of a buried lake in winter; (e) Rosin Threshold Segmentation method; (f) Extraction of a buried lake in winter.

其次，對Sentinel-1的SAR GRD數據集進行預處理，預處理包括：熱噪聲去除、輻射校正和地形校正，并按照月份進行合成，利用擴張后的冰面湖矢量提取其對應的SAR影像.

最后，提取格陵蘭冰蓋次表面湖.由于次表面湖經過消融、掩埋和凍結三個過程，冬季所能提取到的水體范圍要小于夏季最大水體范圍.由于次表面湖與冰面背景的后向散射差異，冰面湖矢量所對應其冬季SAR影像應存在兩種情況：其一，影像僅含有冰面背景，影像直方圖為正態分布；其二，影像同時含有冰面背景和次表面湖，影像中水體占比較小時，直方圖呈現偏峰分布(圖5a)，影像中水體占比較大時，直方圖呈現雙峰分布(圖5b).研究區夏季最大水體范圍掩膜內的SAR影像直方圖以偏峰分布為主，故文章采用基于Rosin閾值分割的方法從雙峰分布或偏峰分布中動態提取次表面湖(Rosin, 2001).

圖5 次表面湖SAR影像直方圖(a) 水體占比較小時SAR影像直方圖； (b) 水體占比較大時SAR影像直方圖.Fig.5 Histogram of SAR image of buried lakes(a) Histogram of SAR image with small water occupation; (b) Histogram of SAR image with large water occupation.

Rosin閾值分割算法適用于研究區存在一個主要地物類型，除此類型之外還會存在一個相對次要地物類型的情況.后一類地物可能會產生可辨別的峰或不可辨別的峰，需要與主峰合理地分開，以免其包含的信息被主峰淹沒(Rosin, 2001).其主要思想是從直方圖的最大值點引出一條直線，連接向左尋找到的第一個零值點，從直方圖中各點向該直線做垂線，距離最大的點所對應的灰度值即為所求的閾值(Rosin, 2001)(圖4e).利用上述方法遍歷所有水體矢量，提取格陵蘭中西部流域冬季次表面湖.

2.4 方法驗證

方法驗證主要分析了利用監督分類方法提取次表面湖的精度及存在的問題，并對自動提取算法的提取結果進行驗證(圖6).

圖6 監督分類法提取次表面湖Fig.6 Buried lakes extracted by supervised classification

〗首先，利用研究區冬季合成Sentinel-1 HV波段數據，結合研究區高程數據進行監督分類，提取冬季格陵蘭中西部流域次表面湖.高程數據采用來源于美國國家地理空間情報局(National Geospatial-Intelligence Agency，NGA)和美國國家科學基金會(National Science Foundation，NSF)的ArcticDEM，分辨率為2 m(Morin et al., 2016).

由表1可得，研究區監督分類的精度隨樣本的數量的減少快速下降，故監督分類方法需在樣本量合適的情況下才可以得到較好的分類精度，但其分類精度嚴重依賴于樣本的選取及樣本選取的質量.

表1 監督分類精度評定Table 1 Accuracy assessment of supervise classification

其次，在格陵蘭中西部流域選取了10個含有次表面湖的測試區域，利用冬季SAR影像，對測試區域內的次表面湖進行手動勾繪，并與提取的次表面湖結果進行對比，圖7a—j代表所選擇的10個測試區域，分析其總體精度(Overall Accuracy，OA)及Kappa系數(表2).

表2 次表面湖提取精度評定Table 2 Accuracy evaluation of buried lakes extraction

10個測試區域的平均總體精度約為99.1%，平均Kappa系數約為0.85.10個測試區域Kappa系數大多處于0.6～0.8和高于0.8兩個區間，屬于高度的一致性和幾乎完全一致這兩個類別，少部分低于0.6，屬于中等的一致性.驗證了利用該方法提取格陵蘭中西部流域次表面湖的可行性.

3 結果

3.1 格陵蘭中西部流域夏季冰面湖空間分布

通過NDWIice和NDSI提取格陵蘭中西部流域融化季節(5—9月)的冰面湖(圖8)，主要分布在沿岸區域，圖8a—f是格陵蘭中西部流域夏季冰面湖提取結果的樣區，左側代表樣區夏季冰面湖可見光影像，右側代表夏季冰面湖的提取結果.2018年夏季探測到冰面湖共352個，總面積為102.28 km2，每個湖的平均面積為0.29 km2，面積主要集中在0.05～0.1 km2，約占總面積分布的38.9%，高程主要集中在900～1400 m，約占總高程分布的87.2%(圖9).

圖8 格陵蘭中西部流域冰面湖結果Fig.8 Supraglacial lakes in the Central West basin of Greenland

圖9 格陵蘭中西部流域冰面湖分布(a) 冰面湖面積分布范圍； (b) 冰面湖高程分布范圍.Fig.9 The distribution of supraglacial lakes in the Central West basin of Greenland(a) The area distribution of the supraglacial lakes; (b) The elevation distribution of the supraglacial lakes.

3.2 格陵蘭中西部流域冬季次表面湖分布

格陵蘭冬季次表面湖的提取采用最大值合成的方法，取2018年12月、2019年1月、2019年2月出現過的次表面湖合成 2018—2019年冬季次表面湖(圖10)，圖10a—f是格陵蘭中西部流域冬季次表面湖提取結果的樣區，左側代表樣區冬季次表面湖SAR影像，右側代表冬季次表面湖的提取結果.2018—2019年冬季探測到次表面湖共166個，總面積為44.07 km2，每個湖的平均面積為0.27 km2，面積主要集中在0.05～0.1 km2，約占總面積分布的40.7%，高程主要集中在900～1400 m，約占總高程分布的92.8%(圖11).

圖10 格陵蘭中西部流域次表面湖結果Fig.10 Buried lakes in the Central West basin of Greenland

圖11 格陵蘭中西部流域次表面湖分布(a) 次表面湖面積分布范圍； (b) 次表面湖高程分布范圍.Fig.11 The distribution of buried lakes in the Central West basin of Greenland(a) The area distribution of the buried lakes; (b) The elevation distribution of the buried lakes.

4 討論

結果表明基于SAR數據和Rosin閾值分割方法的格陵蘭冰蓋次表面湖監測技術可以有效提取冬季格陵蘭次表面湖的范圍，證明格陵蘭中西部流域在冬季仍有大量水體存在于次表面.

此方法在提取次表面湖時，也會存在一定誤判的情況.首先冬季可能存在次表面湖的范圍依賴于夏季所提取的水體范圍，即與NDWIice指數閾值的選取有關.NDWIice指數是提取冰蓋水體較為常用的方法，但其閾值的選擇會存在較大的差異(Miles et al., 2017; Schr?der et al., 2020; Benedek and Willis, 2021).Bell以0.14作為表面水體提取的閾值(Bell et al., 2017)，而Dunmire等(2021)則利用0.5為閾值提取夏季冰面湖.

由于NDWI用來提取格陵蘭水體的閾值集中在 0.14～0.5范圍內(Miles et al., 2017; Schr?der et al., 2020; Benedek and Willis, 2021)，我們以0.05作為步長，對水體指數NDWI做敏感性實驗，分析NDWI變化對次表面湖提取結果的影響(表3).由測試結果可知，隨著NDWI取值的減少，夏季提取的水體范圍會有所增加，導致更多的淺水范圍被劃分為冬季可能存在次表面湖的情況，從而導致誤分的情況.本文僅針對較深水體進行次表面湖的提取，故選擇0.5作為水體提取的閾值.

表3 冰面湖提取閾值的敏感性分析Table 3 Sensitivity analysis of threshold of supraglacial lakes extraction

其次，Rosin閾值分割法在提取冬季次表面湖時也會存在一定誤差.第一，雖然在冬季冰面和次表面水體的后向散射系數存在可分離性，但仍有一定的交集，尤其是在近岸地區，冬季水體與冰面的后向散射系數接近.由于Rosin閾值分割法僅根據每個冰面湖矢量內后向散射系數的概率密度分布得到分割閾值，故在沿岸地區會存在將冰面誤分為次表面湖的情況.第二，若在研究區域內，冬季次表面湖面積小而冰面占比較大，其后向散射信息可能會淹沒于冰面的主峰故難以通過此方法而分離；第三，在僅有冰面的情況，研究區矢量范圍內直方圖呈現正態分布，在此情況下，Rosin閾值分割法也有可能將少數冰面誤分為水體；第四，雖然Rosin閾值分割法可以處理部分雙峰的情況，但對于研究區矢量范圍內兩種地物比例相當或者次表面湖占主體的情況下，該方法也無法準確提取次表面湖；第五，Rosin閾值分割法的精度也依賴于SAR影像的質量.由于研究區內水體矢量范圍影像直方圖以偏峰分布為主，綜上所述，利用Rosin閾值分割法提取次表面湖會存在一定的高估.

由于本文中次表面湖提取主要是依賴Rosin自適應閾值分割的方式實現，故在獲得自適應閾值后，分別對其進行增加或減少1 dB、2 dB、3 dB，對Rosin閾值選取做敏感性實驗，分析其不確定性(表4).由測試結果可知，Rosin所得到的次表面湖范圍對所使用的閾值非常敏感，提取閾值增加或減少會對提取范圍造成很大影響，例如，增加/減少3 dB會對提取結果造成317%/-481%的偏差，這正說明使用全局閾值無法在大尺度上準確提取所有次表面湖.

表4 次表面湖提取閾值的敏感性分析Table 4 Sensitivity analysis of threshold of buried lakes extraction

最后，影像本身也會存在一定的噪聲.影像的后向散射信息與采集的季節性時機、視線方向、入射角、表面和近表面的濕度、介電常數以及表面特征的結構有關(Miles et al., 2017; Schr?der et al., 2020)，對提取的精確度造成影響.并且由于C波段微波雷達在冰雪中的穿透深度為1～10 m(Rignot and Mouginot, 2012)，我們的方法無法探測到超過此深度的水體.此外，利用此方法僅可以提取當年夏季出現過表面融水的次表面湖，無法提取格陵蘭更靠近內陸的區域全年不會暴露于地表的次表面湖.

5 總結與展望

本文提出了一種基于Rosin圖像分割算法，使用Sentinel-1數據自動提取次表面湖范圍的方法，并得到2018—2019年冬季格陵蘭冰蓋中西部流域夏季冰面湖和冬季次表面湖分布范圍，得出結論如下：

(1)本研究提出的次表面湖提取算法的總體精度達到99%以上，Kappa系數達0.85，實現了對格陵蘭中西部流域次表面湖自動提取的功能.該方法能用于整個格陵蘭各流域次表面湖的提取，深化對格陵蘭冰蓋物質平衡的理解.對比機器學習方法，該方法在提取次表面湖時具有易于實現、不需要人工采集樣本且無需人為設定大量閾值的特點.

(2)利用該方法提取2018—2019年格陵蘭中西部流域夏季冰面湖共352個，總面積為102.28 km2，冬季探測到次表面湖共166個，總面積為44.07 km2，說明43.09%的夏季冰面湖在冬季不會完全凍結.

(3)對比夏季冰面湖和冬季次表面湖，其面積和高程的總體分布趨勢一致，均呈現面積主要集中在0.05～0.1 km2，高程主要集中在900～1400 m的規律.在海拔大于1400 m的區域，夏季融化較少，冬季具有隔熱條件，更有利于物質積累；對于海拔低于900 m的區域冰面湖排水系統更發達，不易形成次表面湖.

本研究也存在很多可以改進的地方.首先，Rosin圖像分割算法在提取次表面湖時還存在一定的不確定性，主要體現在研究區內僅有冰這一類地物存在時的情況，利用Rosin閾值分割方法會將處于直方圖末端的少量像元誤分為次表面湖；在兩類地物占比接近時提取精度較低，但此種情況的占比很小.其次，該方法在提取次表面湖時，僅采用2018年12月、2019年1月、2019年2月三個月份作為2018—2019年冬季次表面湖提取的時間范圍，未能分析次表面湖隨時間的變化情況，也未以單個湖體為單位分析其面積變化及監測其排水情況，難以分析其儲量及生消規律等特征，也無法開展其變化機制的研究.最后，由于目前大量可用的SAR數據從2017年開始獲取，難以對次表面湖進行長期動態監測結果的分析.多數次表面湖是一個重新凍結的過程，其對冰蓋穩定性及物質平衡的影響應重點監測其排水過程，隨著數據積累，此項工作將是將來研究的重點.

因此，進一步研究將通過增添直方圖分類情況、增加其他分割算法以及引入同極化數據等方法優化地物類型判斷方法，對格陵蘭次表面湖的演變進行監測并研究其形成及變化機理，分析次表面湖的形成、發展與排水事件對冰蓋穩定性或物質平衡的影響.

致謝感謝中山大學蔣彌副教授對研究提出的建設性意見；感謝美國地質勘探局(United States Geological Survey，USGS)提供的Landsat8衛星數據(https:∥earthexplorer.usgs.gov/)；感謝歐洲航天局(European Space Agency，ESA)提供的Sentinel-1衛星數據(https:∥sentinels.copernicus.eu/web/sentinel/missions/sentinel1)；感謝美國明尼蘇達大學極地空間數據中心提供的ArcticDEM高程數據(https:∥www.pgc.umn.edu/data/arcticdem/).