基于Landsat-8影像的北极地区入海冰川流速监测

刘 欣，张绪冰，王 耀

（中国地质大学（武汉）地理与信息工程学院，湖北武汉430078）

0 引言

入海（marine-terminating）冰川冰流由于注入海洋，其流速变化直接影响入海冰流量和全球海平面升降，对冰川物质平衡具有重要的指示意义。北极地区是地球上主要的冰冻圈，北冰洋中各岛屿沿岸广泛分布着入海冰川，其中面积最大的格陵兰岛拥有北极地区最大的冰盖，体积约为2.99×106km3，若全部消融将导致全球总海平面上升约7.42 m［1］。研究指出，冰川流速的时空变化特征与影响机制复杂，且具有区域性［2-3］；在全球气候变化的环境下，北极地区多数冰川正面临着末端后退和运动加速的问题［4-5］。因此，获取北极地区入海冰川流速，实现其时空变化的遥感监测与规律分析，为研究入海冰川的运动特征提供关键信息，对北极地区冰冻圈的资源开发与生态环境保护具有重要意义。

遥感手段能实现大范围、长时间的冰川流速提取，可分为光学遥感监测［6-10］和雷达遥感监测［11-15］。其中由于Landsat-8相较前代卫星传感器具有更高的辐射分辨率、图像获取率及图像质量，且空间分辨率和重复周期适中，且相比于雷达影像冰川特征更加直观清晰，在基于特征追踪法［16-21］的全球冰川流速监测中表现出良好的效果［22-25］。同时，北极地区冰川流速监测有关研究主要集中在格陵兰岛沿岸［26-28］，虽然有研究单独分析了加拿大埃尔斯米尔冰帽［29］，斯瓦尔巴群岛、北地群岛和法兰士约瑟夫地群岛地区冰川流速［30］，但仍然缺少对北极地区大范围、多数量的入海冰川流速特征的分析，导致难以把握北极地区冰川流速的空间分布特征以及通过入海冰流量对海平面升降和冰川物质平衡具有重要影响的关键区域。

因此，本文基于Landsat-8卫星L1T产品全色波段数据，采用特征追踪法提取冰川流速，并提出了针对北极地区入海冰川的流速异常值剔除方法。通过提取北极地区共198条入海冰川在2017年或2018年消融期流速，以及Kangerlussuaq冰川流速在2018年3—10月消融期间的流速变化，分析了北极地区入海冰川流速的时空变化特征以及可能的影响因素，有助于辨别在北极地区开发过程中需要重点研究关注的冰川区域和时间段，也为北极地区自然生态环境保护、北极航道规划、北极冰川灾害风险评估等研究提供了基础资料。

1 研究区与数据

1.1 研究区概况

本文选取了198处北冰洋沿岸岛屿入海冰川作为研究对象（图1）。其中格陵兰冰盖是仅次于南极冰盖的世界第二大冰盖，冰下基岩呈东高西低，南高北低的地势［1］；斯瓦尔巴群岛（东北地岛）和北地群岛（十月革命岛）也存在区域开阔且规模较大的冰帽。北极地区1月份的平均气温介于-40～-20℃，最暖月8月的平均气温也只有-3℃左右。北冰洋表层环流则主要受到北大西洋暖流及其支流的影响，使沿岸气候相对于内陆更加温暖湿润。现有研究资料显示北极地区入海冰川流速差异较大，例如位于格陵兰岛东海岸Ikeq Fjord冰川流速在2010年可达30 m·d-1以上，但多数入海冰川流速低于5 m·d-1［2］。

图1 研究区域示意图Fig.1 Map showing the study area

1.2 数据源及预处理

1.2.1 实验数据

研究区L1T级别的Landsat-8全色波段影像数据均来源于美国地质调查局（USGS）（https：//earthexplorer.usgs.gov/），该级别影像已经过正射校正，因此冰川表面同名点的位置变化可以被认为是由冰川自身运动导致的。影像空间分辨率为15 m，重复周期为16 d，但由于北极地区纬度较高，对于部分区域的重访天数也有可能短于16 d［31］。格陵兰岛和北地群岛影像获取于2018年7月至9月之间，占本文考察冰川总数的90%；同期少数地区的影像由于受到云层覆盖的影响导致无法成功提取冰川流速，改用2017年7—9月期间的影像，例如斯瓦尔巴群岛、法兰士约瑟夫地群岛和德文岛冰川影像。根据2017年影像提取的冰川运动速度较慢，流速随时间变化幅度较小，本文结果与前人研究也具有较高的一致性［32-33］，因此认为上述198条北极地区入海冰川流速具有可比性。

此外，本文对位于格陵兰岛东海岸68°40′N的Kangerlussuaq冰川连续提取了2018年3—10月的11组表面流速（表1），2018-08-18和2018-09-03冰川流速结果还将与同期的验证数据进行对比和方法的精度验证。

表1 用于Kangerlussuaq冰川流速提取的影像对信息Table 1 Parameters of Landsat data used to retrieve Kangerlussuaq Glacier velocities

对于上述所有Landsat-8全色波段影像，统一采用3×3像素的高斯高通滤波进行预处理，得到增强冰川表面纹理特征的新影像以提高特征追踪中的正确匹配率［22，34］。1.2.2其他数据

本文选用全球陆地冰测量计划（Global Land Ice Measurements from Space，GLIMS）中由雷达影像提取的Kangerlussuaq冰川流速产品作为验证数据；此外，由于数据质量限制，Jakobshavn冰川前缘流速结果缺失，因此将同样采用光学影像提取的2018年8月Jakobshavn冰川流速产品作为补充数据［35］。

2 冰川流速提取

2.1 特征追踪法

本文采用基于图像互相关算法的特征追踪法［9，36-37］对冰川表面同名点进行追踪，并提取冰川表面特征的偏移量。如图2所示，（x，y）为参考模板的像元位置，（x-u，y-v）为搜索模板的像元位置，f和t分别为参考模板和搜索模板的像元值，fˉ和tˉ分别为参考模板和搜索模板的像元平均值。本文选择的归一化互相关系数ρ的计算公式［13］为

图2 特征追踪方法的原理Fig.2 The principle of feature tracking method

本文根据冰川规模由大至小，选择41×41，31×31或21×21像素尺寸的参考模板［34］。北极地区入海冰川最大流速约为35 m·d-1［38］，基于影像间隔时间天数T和影像空间分辨率15 m，将搜索范围设置为（35×T）/15个像元，以保证冰川表面特征的最大像元偏移量在搜索范围之内。每隔4个像元计算像元偏移量，因此所得原始流速结果的分辨率为60 m。

2.2 异常值剔除

根据冰川运动的应力传递规则，冰川表面流速变化通常是连续的，但由于①影像中的云层遮挡，②冰川表面特征在研究期间经历了旋转、剪切、碎裂等较大变化，③流动冰川体和静止基岩之间的剪切运动等原因，部分匹配结果将出现错误。得到冰川表面特征的偏移量之后，根据影像对的空间分辨率和间隔天数计算冰川流速。本文以如图3所示3×3像素的冰川流速像元为单元剔除流速结果中的异常值。其中P表示中心像元，Q表示邻域像元Q1～Q8的集合，O为P、Q之和。本文依据前人研究以及北极地区入海冰川流速实际情况设计了如下针对北极地区入海冰川的流速异常值剔除方法［38-40］：

图3 速度筛选单元示意图Fig.3 Filter unit for outlier eliminating

①输入分辨率为60 m的冰川流速图像，像元值代表冰川流速，单位为m·d-1。

②遍历60 m分辨率的冰川流速图像，当标准差σ（Q）＜0.002 m·d-1且Qˉ-3σ（Q）＜P＜Qˉ+3σ（Q），中心像元P被剔除。

③第二次遍历冰川流速图像，如果σ（O）＞0.002 m·d-1，中心像元P被剔除。

④以3×3像元的中位数作为新的像元值进行像元聚合，图像分辨率由60 m变为180 m；一方面避免数据冗余，另一方面可以有效抑制前一步骤产生的“孤岛”状异常值（见3.1节）。

⑤第三次遍历冰川流速图像，此时根据各像元的邻域像元是否为空，分为以下3种情况：

a.当集合Q为空，中心像元被剔除；

b.当集合Q仅有一个值Qi，|P-Qi|＞1 m·d-1时，中心像元被剔除；

c.当集合Q有2个或2个以上的有效值，当标准差σ（Q）＜0.006 m·d-1且Qˉ-3σ（Q）＜P＜Qˉ+3σ（Q），中心像元被剔除。

⑥第四次遍历冰川流速图像，如果σ（O）＞0.006 m·d-1，中心像元被剔除。

异常值剔除之后，采用克里金插值法得到200 m分辨率的冰川流速图像。

3 实验结果

3.1 异常值剔除结果

以Kangerlussuaq冰川为例，异常值剔除过程如图4所示，其中图中黑色部分代表被剔除的像元，包括海洋区域和异常值。图4（a）为特征追踪后得到的原始冰川流速图，其局部细节如图5（a）所示，原始流速图中有许多不符合冰川运动规律的异常值，需要被剔除。

图4 （b）为步骤③之后的结果，在剔除了原始图像中大部分离散异常值的同时也残留下部分区块异常值，这部分区域通常由多个异常值组成，空缺值将其与周围像元隔开，状似“孤岛”。

图4 （c）是步骤④之后的结果，像元聚合降低了图像分辨率，同时填补了一些较小的空缺值，减少了表示“孤岛”异常值的像元数量，并重新建立起了“孤岛”与周围像元间的联系，并通过进一步筛除有效抑制了“孤岛”异常值得到图4（d）的结果，最后通过克里金插值得到图4（e）的结果。

图4 冰川流速异常值剔除流程（黑色区域为空值）Fig.4 Filtering process of glacier velocity：original velocity（a），first elimination（b），pixel aggregation（c），second elimination（d），and interpolation（e）

异常值剔除处理前后对比如图5所示，冰川体表面的异常值已基本被剔除，已经符合后续分析的要求。

图5 异常值剔除前后冰川流速图像细节变化展示（黑色区域为空值）Fig.5 Velocity field before and after outlier filtering in detail：original velocity（a），before interpolation（b），and after interpolation（c）

3.2 方法验证

本文选择Kangerlussuaq冰川进行方法验证。验证数据为Joughin等［41］采用2018年8月16日和同年9月7日、间隔22 d的TerraSAR-X/TanDEM-X雷达影像采用InSAR方法提取的冰川流速；本文结果采用2018年8月18日和同年9月3日、间隔16 d的Landsat-8光学影像得到。

首先将100 m空间分辨率的验证数据通过双线性插值法重采样为200 m分辨率，如图6所示，沿Kangerlussuaq冰川长达20 km的中流线A—B上等间隔选取100个采样点对比验证数据和本文实验结果的冰川流速值。从图7可以发现，实验结果的变化趋势和验证数据基本吻合，即冰川流速从入海口向内陆在逐渐减慢，在中点以后流速趋于平缓。两者差值在-0.40～0.40 m·d-1之间，平均差值约为0.21 m·d-1，仅占A—B流速最大值22.64 m·d-1的0.93%。鉴于两种流速结果由独立的方法和数据类型得到，认为该差异属于合理范围。

图6 Kangerlussuaq冰川中流线和控制点的位置Fig.6 Center streamline and GCPs of Kangerlussuaq Glacier

图7 实验结果与验证数据的比较Fig.7 Comparison between experimental results and validation data

根据基岩裸露区域，分别选取了29处控制点进行精度验证。流速的均方根误差约为0.0858 m·d-1，仅占流速最大值的0.38%，表明流速提取结果是可信的［18］。

3.3 北极地区冰川流速分析

3.3.1 空间分布特征

本文提取了北极地区198条入海冰川在2017年或2018年消融期的表面流速，并对冰川入海口处2 000 m中流线的流速平均值（后文简称前缘流速）进行统计和比较分析。本文根据北极地区入海冰川实际流速情况定义前缘流速为0～5 m·d-1的入海冰川为低速运动冰川，5～10 m·d-1为中低速运动冰川，10～20 m·d-1为中高速运动冰川，20～35 m·d-1为高速运动冰川。

如图8所示，格陵兰岛北海岸入海冰川均为低速运动冰川，格陵兰岛东、西海岸考察冰川数量接近，但东海岸前缘流速为5 m·d-1以上的入海冰川约占总数的46.92%，而格陵兰岛西海岸以低速运动冰川为主，占总数的73.17%。同时东海岸考察入海冰川的平均前缘流速6.13 m·d-1高于西海岸的4.14 m·d-1，最大前缘流速达到31.62 m·d-1，高于西海岸唯一高速运动冰川Jakobshavn的前缘流速（26.33 m·d-1）。

图8 格陵兰岛2018年消融期入海冰川前缘流速Fig.8 Front velocities of Greenland marine-terminating glaciers during 2018 ablation season

如图9所示，格陵兰岛以外北极地区入海冰川前缘流速均在10 m·d-1以下，且2条中低速运动冰川分别位于斯瓦尔巴群岛（东北地岛）和北地群岛（十月革命岛）的冰帽入海口处，前缘流速分别为7.96 m·d-1和7.61 m·d-1。

图9 北极地区其他入海冰川前缘流速空间分布Fig.9 Other marine-terminating glaciers’front velocities in Arctic：Franz Josef Land（a），Svalbard（b），Severnaya Zemlya（c），and Devon（d）

综上所述，北极地区入海冰川平均流速由大到小依次为：格陵兰岛（5.00 m·d-1）、北地群岛（3.07 m·d-1）、斯瓦尔巴群岛（2.22 m·d-1）、法兰士约瑟夫地群岛（1.80 m·d-1）、德文岛（0.92 m·d-1）。格陵兰岛部分冰川，以及北地群岛、斯瓦尔巴群岛的冰帽冰川由于流速较快，对海洋环境影响较大，需要进行重点研究。

3.3.2 时间变化特征

除了对北极地区198条入海冰川流速的空间分布特征进行分析以外，本文还提取了Kangerlussuaq冰川2018年3—10月的13组流速结果，并沿冰川中流线A—B以200 m等间隔设置了150个采样点对流速值进行提取，同时基于表1的控制点流速的平均值对冰川流速进行微调，并分析其时空变化特征。

如图10所示，Kangerlussuaq冰川中流线在不同时期的流速空间变化趋势较为一致，流速从入海口向内陆逐渐减慢。距离入海口A点0～10 km中流线减速幅度相对于距离入海口A点10～30 km更为剧烈。2018年6—7月期间的前缘流速（距离入海口A点2 km范围内11个采样点的平均流速）和整体流速（中流线150个采样点的平均流速）最快，分别为22.87 m·d-1和11.39 m·d-1；2018年8—9月期间的前缘流速和整体流速最慢，分别为21.02 m·d-1和10.02 m·d-1。

图10 Kangerlussuaq冰川中流线流速随时间变化Fig.10 The temporal changes of center streamline velocities of Kangerlussuaq Glacier

如图11和图12所示，Kangerlussuaq冰川前缘流速和整体流速在研究期间的变化基本一致（实验序号见表1），3—6月期间流速变化不明显，在6—7月冰川出现了第一个和最高的流速峰值，随后在8—9月出现了流速低谷，而9—10月重新攀上第二个流速峰值，前缘流速和整体流速分别为22.18 m·d-1和11.15 m·d-1，低于6—7月数值。10月以后冰川呈现减慢的趋势。由此可知，Kangerlussuaq冰川流速变化较为复杂，除了冰川融水对冰川运动具有润滑加速的作用，冰川前缘崩解导致的冰川物质平衡变化也可能对流速产生影响［42］。

图11 Kangerlussuaq冰川前缘流速随时间变化Fig.11 Front velocity change of Kangerlussuaq Glacier between March and October in 2018

图12 Kangerlussuaq冰川整体流速随时间变化Fig.12 Average velocity change of Kangerlussuaq Glacier between March and October in 2018

4 讨论

4.1 北极地区入海冰川流速时空变化

基于2017/2018年消融期的多景Landsat-8遥感影像，本文获取了北极地区198条入海冰川流速及其空间分布特征：整体而言，格陵兰岛沿岸冰川流速较快，欧亚大陆北部的斯瓦尔巴群岛、北地群岛、法兰士约瑟夫地群岛冰川，美洲北部的德文岛冰川流速较慢。其中格陵兰岛海岸的流速差异与近年的研究结果具有较高的一致性［2-3］。

对于格陵兰岛，其北海岸入海冰川主要注入北冰洋，前缘流速均在5 m·d-1以下，高纬度及严寒的气候条件导致该海岸区域平均流速最低，约为1.99 m·d-1。西海岸入海冰川主要注入巴芬湾，以低速运动冰川为主，占总数的73.17%；东海岸入海冰川主要注入格陵兰海和大西洋，低速运动冰川占总数的53.09%，中低速、中高速和高速运动冰川数量将近考察冰川总数的一半。格陵兰岛东海岸入海冰川平均流速约为6.13 m·d-1，高于西海岸的4.14 m·d-1。东西海岸的流速差异一方面受到海流作用的影响，如格陵兰岛以东的海域由于受到北大西洋暖流的影响，一定程度上促进了东海岸入海冰川的消融和运动；另一方面，由于东海岸基岩地势高于西海岸，与海平面更大的地势差也会促进冰川体的滑动。

由13景遥感影像得到的11组Kangerlussuaq冰川流速结果的分析可知：中流线流速从入海口向内陆逐渐减慢；冰川前缘流速和整体流速在2018年6—7月和9—10月均出现明显流速峰值，且6—7月的前缘流速和整体流速最高，分别为22.87 m·d-1和11.39 m·d-1，而在8—9月前缘流速和整体流速最低，分别为21.02 m·d-1和10.02 m·d-1。

4.2 冰川流速空间分布差异原因

冰川流速空间分布差异的影响因素可能有以下几点：

（1）冰床地势

冰川运动包括塑性变形和底部滑动，冰川表面或冰床高度的不同会让冰川体在重力或压力的驱使下自地面高处（或冰层厚处）流向地面低处（或冰层薄处），冰床地势高低因而会影响冰川体的流动速度。根据Morlighem等［1］构建的格陵兰岛基岩高程模型可知，格陵兰东海岸地势相对于西海岸较高，或许是导致格陵兰东海岸相较于西海岸整体平均流速更大、中高速和高速运动冰川数量更多的原因之一。此外，较为崎岖的冰床地势对冰川流动也具有明显的阻滞作用［43］。

（2）海流作用

北极地区入海冰川的主要特点是冰川前缘以及冰架部分与海水接触频繁，而且随着北极海洋温度变暖［44-46］，入海冰川前缘与海水的接触会加速冰川的底部融化并形成底部凹槽，导致冰架变薄崩解、冰川流动加速、冰川前缘后退等现象。受到北大西洋暖流的影响，格陵兰岛以东的海域以及斯瓦尔巴群岛东南部温度相对较高，促进了入海冰川的前缘消融，一定程度上影响了冰川流速的空间分布［47］。

（3）冰盖消融和冰盖表面水文系统状况

在格陵兰冰盖消融理想模式中，融水通过冰盖表面冰裂隙、冰川竖井等进入冰盖内部甚至到达冰盖底部后，能够起到润滑作用，加速冰盖运动［48］，因此区域冰川融水越丰富，冰盖表面水文系统越发育，冰川的运动速度也有可能更快。但也有部分学者认为融水对冰川的加速作用并不显著甚至实际上减缓了冰川运动速度［49-50］。

综上所述，北极地区入海冰川流速空间分布差异是不同气温降水、冰床地势、海流作用和冰盖消融状况等复杂因素综合影响的结果。此外，前缘流速较快的入海冰川更容易出现在冰盖/冰帽的溢出部分，例如集中在格陵兰岛沿岸的中高速和高速运动冰川，和位于东北地岛和十月革命岛的中低速运动冰川。

5 结论与展望

本文基于Landsat-8卫星影像数据，采用特征追踪方法提取了北极地区入海冰川在2018或2017年消融期的冰川流速，并提出了具有区域针对性的流速异常值剔除方法，得到的Kangerlussuaq冰川中流线流速与验证数据的平均偏差仅为0.21 m·d-1，仅占中流线A—B流速最大值22.64 m·d-1的0.93%；29处控制点流速的均方根误差约为0.0858 m·d-1，仅占流速最大值的0.38%，证明了本文流速监测结果可信。本文分析了北极地区198条入海冰川流速的空间分布差异和影响因素，以及Kangerlussuaq冰川流速在消融期的时空变化特征，结论如下：

（1）北极地区入海冰川流速的空间差异主要体现为冰盖/冰帽外围处入海冰川流速相对较快，如格陵兰冰盖和东北地岛及北地群岛的冰帽；而零散的小型冰川流速相对较慢。北极地区冰储量最大的格陵兰岛冰川流速空间差异则体现为东海岸入海冰川平均流速最高，西海岸冰川次之，北海岸入海冰川平均流速最低。因此一方面有必要对冰盖/冰帽边缘流速较快入海冰川进行长期、多频次的监测，并重点关注其对气候变化的响应。

（2）格陵兰岛东海岸Kangerlussuaq冰川中流线流速从入海口向内陆冰川流速逐渐减慢。2018年3月至10月期间，冰川前缘流速和整体流速在2018年6—7月和9—10月均出现明显流速峰值，而在8—9月出现低谷，变化趋势较复杂。冰川融水作用，以及随着冰川前缘崩解的冰川物质平衡变化可能是变化趋势复杂的主要原因。为了更加深入理解冰川运动机制，在今后的研究中有待加大观测密度和时间段。

（3）北极地区入海冰川流速空间分布差异是由多种因素综合导致的结果，包括冰川规模、冰床地势、冰川消融状况、海流作用等，影响机制复杂，有待于后续结合冰床地势、冰川形态、海水温度等相关数据资料进行相关性研究。