格陵兰冰盖次表面湖多源遥感监测——以格陵兰中西部流域为例

李岚静， 陈卓奇,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/).