基于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）北極地區入海冰川流速空間分布差異是由多種因素綜合導致的結果，包括冰川規模、冰床地勢、冰川消融狀況、海流作用等，影響機制復雜，有待于后續結合冰床地勢、冰川形態、海水溫度等相關數據資料進行相關性研究。