張齊民， 閆世勇， 呂明陽， 張 露， 劉 廣

（1.自然資源部國土環(huán)境與災(zāi)害監(jiān)測重點(diǎn)實(shí)驗(yàn)室，江蘇徐州 221116； 2.中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院，江蘇徐州 221116；3.中國科學(xué)院空天信息創(chuàng)新研究院，北京 100094）

0 引言

近年來，全球氣候變暖的進(jìn)程日益加劇，作為氣候變化研究的關(guān)鍵指標(biāo)，山地冰川變化受到全世界人們的廣泛關(guān)注。運(yùn)動是冰川的重要特性之一，和氣候因素、地形地勢等存在直接關(guān)系，與冰川物質(zhì)平衡變化密切相關(guān)［1-2］。相比于兩極冰川，“世界第三極”—青藏高原地區(qū)的冰川對于氣候變化的響應(yīng)更為敏感［3-4］。以青藏高原為主的高海拔地區(qū)被稱為亞洲高山區(qū)，其冰川的動態(tài)變化對該區(qū)域物質(zhì)平衡及能量循環(huán)具有重要的影響，研究該地區(qū)冰川流速時(shí)空變化對于認(rèn)識冰川動力學(xué)特征及其對氣候變化的響應(yīng)具有重要的意義［5-6］。

山地冰川所處區(qū)域氣候條件惡劣，傳統(tǒng)實(shí)地冰川運(yùn)動測量方法難以廣泛開展。遙感技術(shù)監(jiān)測范圍廣、獲取數(shù)據(jù)快、投入成本低，主要有光學(xué)遙感和雷達(dá)遙感兩種觀測方式。光學(xué)影像存檔數(shù)據(jù)多，覆蓋范圍廣，但光照、云霧等因素對其適用性有較大的限制，破壞了冰川監(jiān)測的完整性［4］；相比較而言，合成孔徑雷達(dá)（SAR）能全天時(shí)全天候?qū)ρ芯繀^(qū)進(jìn)行監(jiān)測，不受天氣條件影響且擁有較高的時(shí)空分辨率，在冰川活動性監(jiān)測中有較強(qiáng)的適用性，已成為研究冰川活動性重要的手段之一［7-9］。目前，亞洲高山區(qū)冰川運(yùn)動信息的獲取絕大多數(shù)采用國外遙感衛(wèi)星數(shù)據(jù)，我國自主研制的遙感衛(wèi)星數(shù)據(jù)在冰川運(yùn)動監(jiān)測中的應(yīng)用較少，導(dǎo)致其在冰川活動性研究中科學(xué)價(jià)值無法充分體現(xiàn)。

基于SAR 影像強(qiáng)度信息的偏移量跟蹤算法（pixel tracking，PT）是研究山地冰川運(yùn)動的主要手段之一［10-12］，該方法無須進(jìn)行復(fù)雜的相位解纏工作，對數(shù)據(jù)的時(shí)空基線要求較低，且能直接得到山地冰川在雷達(dá)方位向和距離向上的二維位移分布，對于流速較快的山地冰川運(yùn)動而言具有較好的適用性［13-15］。PT 方法基于歸一化互相關(guān)匹配算法（normalized cross correlation，NCC），計(jì)算精度較高且對高斯白噪聲具有較強(qiáng)的抵抗作用，但其巨大的運(yùn)算量導(dǎo)致效率低下，難以滿足廣域長周期冰川運(yùn)動研究的需要［16］。另外，亞洲高山區(qū)海拔較高，地形起伏較大，導(dǎo)致偏移量跟蹤監(jiān)測結(jié)果中存在大量與地形相關(guān)的偏移量，雷達(dá)圖像斜距投影的幾何特點(diǎn)使得地形偏差效應(yīng)表現(xiàn)的更為明顯［4，17-18］。

為了充分發(fā)揮出我國遙感數(shù)據(jù)在冰川運(yùn)動研究中的優(yōu)勢，及時(shí)準(zhǔn)確地分析亞洲高山區(qū)典型冰川的活動性，本文選用2019—2020年獲取的高分三號（GF-3）遙感影像，基于并行化偏移量跟蹤技術(shù)，通過全局形變擬合估計(jì)、地形幾何畸變校正和降噪濾波等處理后，提取了亞洲高山區(qū)典型冰川表面運(yùn)動特征，為分析青藏高原地區(qū)冰川活動性與氣候變化之間復(fù)雜的響應(yīng)關(guān)系提供了更好的數(shù)據(jù)支持。同時(shí)，基于非冰川區(qū)穩(wěn)定的假設(shè)，通過開展非冰川區(qū)偏移量殘差統(tǒng)計(jì)分析，并對比GF-3 和Landsat-8 影像提取的冰川流速結(jié)果，評估我國遙感衛(wèi)星數(shù)據(jù)在山地冰川運(yùn)動監(jiān)測中的適用性和可靠性。

1 研究區(qū)概況

亞洲高山區(qū)地域遼闊，受到地形因素、冰川自身屬性與規(guī)模等影響，不同區(qū)域冰川末端的進(jìn)退、面積變化及物質(zhì)平衡等狀況有所不同，不同地理單元的冰川運(yùn)動分布具有不同特征，間接反映出不同區(qū)域氣候的差異性［19］。本文選取了天山地區(qū)南以內(nèi)里切克冰川（South Inylchek）、喀喇昆侖地區(qū)卻哥隆瑪冰川（Chogolungma）、巴爾托洛冰川（Baltoro）及橫斷山區(qū)的雅弄冰川（Yanong）進(jìn)行綜合分析（圖1）。天山山脈是世界七大山系之一，極高的地形地貌和充足的降水為冰川發(fā)育提供了有利的條件，是我國最大的現(xiàn)代冰川作用區(qū)之一［20］；在全球變暖的大背景下，多數(shù)區(qū)域的冰川處于緩慢消融狀態(tài)，而喀喇昆侖地區(qū)的冰川面積幾乎不變，被稱為“喀喇昆侖異?！爆F(xiàn)象［3］；橫斷山區(qū)位于青藏高原東南部，受海洋性季風(fēng)氣候的影響，雨季較長，該區(qū)域冰川物質(zhì)平衡線較低，運(yùn)動速度較快［21］。上述典型冰川所處區(qū)域空間跨度大、氣候差異顯著，且所選冰川具有一定的代表性，能夠較好地反映出廣域范圍內(nèi)冰川運(yùn)動與氣候變化之間的內(nèi)在關(guān)聯(lián)。本文根據(jù)典型冰川所在區(qū)域的氣候特點(diǎn)，綜合可利用的雷達(dá)遙感數(shù)據(jù)對冰川運(yùn)動狀況展開分析研究，其中冰川邊界信息參照中國第二次冰川編目數(shù)據(jù)和RGI 6.0數(shù)據(jù)庫。

圖1 研究區(qū)冰川示意圖Fig.1 Overview of the study area

2 數(shù)據(jù)源與研究方法

2.1 研究區(qū)數(shù)據(jù)

本研究選用2019—2020 年GF-3 衛(wèi)星獲取的FSI 模式SAR 數(shù)據(jù)開展亞洲高山區(qū)典型冰川表面運(yùn)動分析工作。GF-3 是我國自主研制的首顆C 波段全極化星載SAR 衛(wèi)星，最高分辨率可達(dá)1 m，重訪周期為29 天。該衛(wèi)星具有12 種常規(guī)成像模式，可滿足不同尺度冰川運(yùn)動特征提取的需要［15］。詳細(xì)數(shù)據(jù)信息如表1 所示，由于研究區(qū)冰川數(shù)目較多，為方便起見，僅用相應(yīng)一條主體冰川代表各個區(qū)域。此外，為更好地反映典型冰川運(yùn)動與氣候變化之間的響應(yīng)關(guān)系，圖2 給出研究區(qū)在研究周期內(nèi)的氣象資料，源自國家氣象科學(xué)數(shù)據(jù)中心（http：//data. cma. cn）提供的中國地面氣溫、降水月值0.5°×0.5°格點(diǎn)數(shù)據(jù)集，產(chǎn)品已經(jīng)過交叉驗(yàn)證和誤差分析，質(zhì)量狀況良好［22-23］，可輔助于冰川運(yùn)動的后續(xù)分析工作。

表1 研究區(qū)GF-3 SAR影像基本信息Table 1 Basic information of GF-3 satellite SAR image in the study area

圖2 研究區(qū)氣象資料：南以內(nèi)里切克冰川（a）；卻哥隆瑪冰川（b）；巴爾托洛冰川（c）；雅弄冰川（d）Fig.2 Temperature and precipitation data of glacier region：South Inylchek gacier（a）；Chogolungma glacier（b）；Baltoro glacier（c）；Yanong glacier（d）

2.2 偏移量跟蹤算法

對于流速較快的山地冰川而言，偏移量跟蹤算法具有較強(qiáng)的適用性。該方法以歸一化互相關(guān)算法為基礎(chǔ)，利用主輔影像間非冰川區(qū)穩(wěn)定特征點(diǎn)的亞像素級匹配結(jié)果，精確解算二者之間的幾何變換模型，進(jìn)而估算冰川運(yùn)動分布。針對傳統(tǒng)的偏移量跟蹤算法計(jì)算效率較低的現(xiàn)狀，本研究基于CUDA架構(gòu)高性能計(jì)算實(shí)現(xiàn)了聯(lián)合GPU 多運(yùn)算核心和CPU 多線程的并行化偏移量算法，極大地提高了運(yùn)行效率；此外，對于成像過程中由于衛(wèi)星姿態(tài)及位置差異引入的全局性偏移量、冰川分布區(qū)高程變化引起的地形相關(guān)殘差以及無關(guān)噪聲等進(jìn)行了相應(yīng)的補(bǔ)償，確保了冰川表面運(yùn)動估計(jì)具有較高的準(zhǔn)確性，其相應(yīng)的數(shù)據(jù)處理流程如圖3所示。

圖3 并行化偏移量跟蹤算法處理流程Fig.3 Flowchart of parallel program design of PT

歸一化互相關(guān)算法是一種經(jīng)典的統(tǒng)計(jì)匹配算法，通過計(jì)算模板圖像和匹配圖像之間的歸一化互相關(guān)系數(shù)來決定各像素的匹配程度［24］。在初始匹配的過程中，首先在主影像上選取模板窗口，而后模板窗口在輔影像的搜索窗口內(nèi)不斷滑動，遍歷計(jì)算模板窗口和搜索窗口的歸一化互相關(guān)系數(shù)，得到相關(guān)矩陣并以其最大值對應(yīng)的位置為最佳匹配窗口，得到相應(yīng)的位移量，進(jìn)而完成初級配準(zhǔn)。為使配準(zhǔn)精度能夠達(dá)到亞像素級，還需要對影像和相關(guān)矩陣進(jìn)行過采樣處理。通過上述步驟能夠獲取研究區(qū)內(nèi)各個模板的總體偏移量，其中包含冰川運(yùn)動相關(guān)偏移量、因軌道差異而引入的全局偏移量、地形起伏帶來的偏移量及噪聲偏移量等信息?；诒O(jiān)測周期內(nèi)非冰川區(qū)未發(fā)生移動的假設(shè)，精準(zhǔn)地估算出冰川運(yùn)動無關(guān)偏移量并對其進(jìn)行補(bǔ)償，進(jìn)而得到可靠的冰川表面流速結(jié)果［25-26］。

在衛(wèi)星成像時(shí)，由于傳感器姿態(tài)和軌道差異導(dǎo)致的主輔影像中同名點(diǎn)偏移一般可以通過多項(xiàng)式擬合進(jìn)行估計(jì)與補(bǔ)償［27］。本研究采用隨機(jī)抽樣一致性算法（random sample consensus，RANSAC），通過不斷的迭代來逐次確定符合變換模型的匹配點(diǎn)，并重新優(yōu)化計(jì)算變換參數(shù)，最終在滿足閾值條件時(shí)終止迭代計(jì)算［28］，準(zhǔn)確擬合并去除衛(wèi)星軌道和姿態(tài)等差異導(dǎo)致的全局形變，從而降低錯誤匹配點(diǎn)對估計(jì)全局變換參數(shù)的干擾［4，16］；此外，成像位置差異會導(dǎo)致偏移量跟蹤結(jié)果中存在與地形相關(guān)的偏移量，一般來說，當(dāng)兩景影像垂直基線小于300 m 時(shí)，無需進(jìn)行復(fù)雜的地形校正［14，29］。然而，受到GF-3 軌道穩(wěn)定性、定位精度等方面的制約，本研究所選用的影像對垂直基線較長（表1），研究區(qū)地勢起伏較大，亟需進(jìn)行地形補(bǔ)償工作。距離向地形起伏導(dǎo)致的偏移量與雷達(dá)成像的垂直基線長度及所在區(qū)域地形高差成正比，而方位向上地形偏移量則主要與所在區(qū)域的地形起伏、成像時(shí)主輔影像之間的軌道交角及衛(wèi)星下視角有關(guān)［30］?；谕獠緿EM 數(shù)據(jù)和兩景SAR 影像的空間基線信息，能夠準(zhǔn)確地獲取地形起伏在各方向上產(chǎn)生的偏移，對結(jié)果進(jìn)行地形效應(yīng)補(bǔ)償，從而準(zhǔn)確地估計(jì)冰川表面運(yùn)動。另外，研究區(qū)冰川表面積雪覆蓋容易導(dǎo)致冰川運(yùn)動場中存在大量的噪聲，在保證冰川邊緣位移信息的完整性和結(jié)果準(zhǔn)確性的前提下，本文采取自適應(yīng)中值濾波（rankedorder based adaptive median filter，RAMF）來削弱噪聲的影響［16］。借助濾波窗口可變特點(diǎn)，RAMF 既能保證冰川偏移量信息的精度，又能適當(dāng)?shù)叵魅踉肼暎诒ㄟ\(yùn)動信息提取中具有較高的實(shí)用價(jià)值。

3 結(jié)果與分析

3.1 偏移量跟蹤算法優(yōu)化分析

3.1.1 并行化偏移量跟蹤算法處理分析

本研究數(shù)據(jù)處理平臺內(nèi)存為128 G，CPU 主頻2.1 GHz，顯卡內(nèi)存12 GB，包含3840 個計(jì)算核心（CUDA），峰值浮點(diǎn)性能為12T Flops，CUDA計(jì)算能力為6.1，滿足GPU 開發(fā)及大數(shù)據(jù)高性能計(jì)算的需要。為保證冰川表面運(yùn)動采樣點(diǎn)的可靠性，偏移量跟蹤算法中各參數(shù)的選擇和設(shè)定需要充分考慮研究區(qū)冰川運(yùn)動特征、規(guī)模類型及其地表后向散射的穩(wěn)定性?；贕PU 并行化偏移量跟蹤算法與傳統(tǒng)偏移量跟蹤算法的參數(shù)設(shè)置及計(jì)算效率對比如表2所示，各冰川區(qū)圖像尺寸及窗口大小均表示為長×寬，單位為像素。結(jié)果顯示各影像對并行化加速比均超過20，極大地提高了計(jì)算效率，為大規(guī)模監(jiān)測冰川運(yùn)動時(shí)空變化特征提供了良好的硬件和軟件保障。

表2 并行化PT算法加速效率對比參數(shù)設(shè)置Table 2 Efficiency comparison between traditional and parallel PT algorithm

3.1.2 冰川相關(guān)偏移量準(zhǔn)確估計(jì)

采用RANSAC 算法對軌道相關(guān)的偏移量進(jìn)行全局?jǐn)M合估計(jì)，基于主輔影像的空間基線信息及外部DEM 對地形起伏引入的地形偏移量進(jìn)行精確補(bǔ)償，最后對上述結(jié)果進(jìn)行濾波處理，完整地提取出亞洲高山區(qū)典型冰川表面運(yùn)動分布特征。經(jīng)過SAR影像高精度配準(zhǔn)和全局形變估計(jì)與補(bǔ)償之后，殘余的誤差主要為地形效應(yīng)造成的偏差（圖4）。黑色輪廓為冰川矢量邊界，圖4（a1）～（d1）為距離向結(jié)果，圖4（a2）～（d2）為方位向結(jié)果，圖示中正負(fù)號表示冰川運(yùn)動方向。

圖4 冰川區(qū)地形相關(guān)位移：南以內(nèi)里切克冰川（a）；卻哥隆瑪冰川（b）；巴爾托洛冰川（c）；雅弄冰川（d）Fig.4 The topography-related displacement of glacier region：South Inylchek gacier（a）；Chogolungma glacier（b）；Baltoro glacier（c）；Yanong glacier（d）

據(jù)圖4 可知，不同方向上地形相關(guān)偏移量差異性較大，距離向殘余偏移量明顯高于方位向，且隨垂直基線的增大而增大，地形起伏嚴(yán)重的區(qū)域，造成的非冰川運(yùn)動相關(guān)位移可達(dá)像素級。為進(jìn)一步分析地形起伏對初始偏移量結(jié)果的影響，沿各區(qū)域非冰川區(qū)取一條剖面線（圖1 中黃色虛線，方向?yàn)?→2），繪制距離向和方位向的地形效應(yīng)偏差曲線如圖5 所示。圖5（a1）～（d1）為距離向結(jié)果，圖5（a2）～（d2）為方位向結(jié)果，正負(fù)號表示方向。

圖5 地形效應(yīng)校正結(jié)果：南以內(nèi)里切克冰川（a）；卻哥隆瑪冰川（b）；巴爾托洛冰川（c）；雅弄冰川（d）Fig.5 Correction results of topographic effects：South Inylchek gacier（a）；Chogolungma glacier（b）；Baltoro glacier（c）；Yanong glacier（d）

結(jié)果顯示，在地形校正之前，由于研究區(qū)地勢起伏過大且空間基線較長，偏移量結(jié)果受地形效應(yīng)誤差影響較為嚴(yán)重。距離向上地形偏移量與垂直基線及地勢高低起伏成正比，在地勢復(fù)雜的區(qū)域不可忽略其影響。南以內(nèi)里切克冰川及卻哥隆瑪冰川區(qū)域，距離向上由地形效應(yīng)所引起的誤差在2 m左右，而巴爾托洛冰川及雅弄冰川區(qū)，由于周圍地形起伏過大，在距離向上誤差可達(dá)10 m 左右，基于外部DEM 和垂直基線數(shù)據(jù)對距離向地形偏差進(jìn)行精確擬合后，非冰川區(qū)位移量結(jié)果整體趨于0，校正效果顯著；方位向上各研究區(qū)由地形因素所帶來的殘差效應(yīng)在0.5 m 左右，地形偏移量量級較小，對其進(jìn)行地形校正之后整體結(jié)果更接近真實(shí)值。由于影響方位向地形偏移量大小的軌道交角及衛(wèi)星入射角量級較小，一般情況下可以忽略該方向的地形效應(yīng)殘差影響［31］。

3.2 冰川運(yùn)動分布

利用2019—2020 年覆蓋研究區(qū)的GF-3 影像數(shù)據(jù)，基于改進(jìn)的偏移量跟蹤算法獲取了典型冰川表面運(yùn)動分布信息，將不同方向冰川流速合成二維平面流速圖［圖6（a）～（d）］。此外，為更好地分析冰川運(yùn)動變化特征，自冰川下游至上游沿主干取一條剖面線（圖1中紅色曲線，方向?yàn)?→2），繪制相應(yīng)區(qū)域冰川運(yùn)動變化及高程變化曲線［圖6（e）～（h）］。四個區(qū)域依次為南以內(nèi)里切克冰川、卻哥隆瑪冰川、巴爾托洛冰川及雅弄冰川。

圖6（a）、6（e）表示2019 年7 月—2020 年6 月天山地區(qū)南以內(nèi)里切克冰川及其附近冰川運(yùn)動分布情況。南以內(nèi)里切克干流冰川整體上呈現(xiàn)為中上游流速較快，下游至末端流速較為緩慢的態(tài)勢，冰川流速方向?yàn)樽詵|向西，與整體高程下降趨勢保持一致。距干流區(qū)末端8～25 km 處運(yùn)動速率在30～35 cm·d-1左右，研究周期內(nèi)累積位移量約為120 m；下游區(qū)域至冰川末梢（0～8 km）冰磧物覆蓋現(xiàn)象嚴(yán)重且處于不斷堆積的狀態(tài)，阻礙冰川繼續(xù)向前移動，加之基底地形較為平整，流速呈現(xiàn)出逐漸遞減的態(tài)勢。在距S1 點(diǎn)約8 km 處將冰川運(yùn)動整體劃分為兩部分，自上游至該處前端流速較為穩(wěn)定，而后出現(xiàn)小的流速峰值，最后流速逐漸減慢。該區(qū)域位于冰川主干與麥茲巴赫（Merzbacher）冰湖所在冰川分支的交界處，來自中上游的冰川不斷輸送至此，應(yīng)力積累導(dǎo)致冰川運(yùn)動速度有所加快，在向北轉(zhuǎn)向麥茲巴赫冰湖處集中釋放；轉(zhuǎn)彎處至冰川末端則由于麥茲巴赫冰湖的分流作用及冰磧物不斷堆積導(dǎo)致速度逐漸減小，其運(yùn)動變化情況與之前學(xué)者研究成果基本一致［32-33］。

圖6 典型冰川運(yùn)動分布：南以內(nèi)里切克冰川（a）（e）；卻哥隆瑪冰川（b）（f）；巴爾托洛冰川（c）（g）；雅弄冰川（d）（h）Fig.6 Distribution of ice motion on typical glacier：South Inylchek gacier（a）（e）；Chogolungma glacier（b）（f）；Baltoro glacier（c）（g）；Yanong glacier（d）（h）

圖6（b）、6（f）表示2019 年7 月—2020 年5 月喀喇昆侖地區(qū)卻哥隆瑪冰川及周圍冰川的運(yùn)動分布情況。卻哥隆瑪冰川整體呈西北東南走向，主干區(qū)域運(yùn)動場保持的比較完整，流速隨海拔的降低逐漸減慢，最大高差超過1 000 m，冰川上游最大運(yùn)動速率約為60 cm·d-1，末端流速約為15 cm·d-1。在距冰川末端剖線起點(diǎn)C1 約11 km、16 km 及20 km 處冰川速率明顯高于周圍區(qū)域，這三處是支流干流冰川的交匯點(diǎn)，來自西南部的支流匯入干流，使得交匯點(diǎn)處冰川速度激增。希斯帕爾冰川位于卻哥隆瑪冰川東北部，其走向規(guī)模與卻哥隆瑪冰川較為相似，但運(yùn)動情況差異性較大。希斯帕爾冰川是該區(qū)域較為典型的躍動冰川，躍動周期較長，干流中部區(qū)域曾在2015—2016年間發(fā)生劇烈的躍動現(xiàn)象，期間觀察到的最大流速約為14 m·d-1，表面特征前移量較大，其余時(shí)間均保持近乎靜止的狀態(tài)［34］。2019—2020 年研究周期內(nèi)，中部躍動區(qū)流速較小，干流兩端及支流處流速相對較大，約為15～20 cm·d-1左右。

圖6（c）、6（g）展示了2019 年8 月—2020 年4 月喀喇昆侖地區(qū)巴爾托洛冰川的運(yùn)動分布情況。該冰川是典型的樹枝狀冰川，干流支流層次分明，主干呈東西走向。受到影像數(shù)據(jù)覆蓋范圍的限制，本文研究范圍主要位于巴爾托洛冰川中下游區(qū)域，即Urdukas 附近區(qū)域至冰川末端部分［27］。研究區(qū)冰川干流自東向西流速逐漸減小，遠(yuǎn)端最大流速達(dá)25 cm·d-1左右，冰川末梢最小流速低于5 cm·d-1，Wendleder、Sun 等人在該區(qū)域觀察到類似的冰川運(yùn)動情況［27，35］，同Mayer 等［36］采用GPS 進(jìn)行實(shí)地測量的數(shù)據(jù)基本吻合，佐證了本研究的結(jié)果具有較高的可靠性。由于該區(qū)域地形高差較大，北側(cè)支流平均速率高于南側(cè)，在支流冰川與干流冰川交匯處冰川運(yùn)動速率較高。此外，由于冰川中部消融區(qū)兩側(cè)山體的阻尼作用及冰川厚度向兩側(cè)減薄，流速在空間分布上呈現(xiàn)出由中間向兩側(cè)逐漸減弱的趨勢，這與多數(shù)山地冰川運(yùn)動規(guī)律特征相一致。

圖6（d）、6（h）為橫斷山區(qū)雅弄冰川上游部分2020 年1 月—2020 年12 月的運(yùn)動分布情況。冰川主冰流的速度約為35～60 cm·d-1，剖線Y1 點(diǎn)到Y(jié)2點(diǎn)運(yùn)動速度先緩慢增加，后逐漸減小，最后急劇增加。剖線末端觀察到最大運(yùn)動速度約為62 cm·d-1，在海拔4 620 m 處觀察到最小運(yùn)動速度約為34 cm·d-1。第一速度峰值出現(xiàn)在約1.3 km 處，這可能是受到支流冰川匯入的影響；第二峰值速度出現(xiàn)在約3.6 km 處，高程變化是造成此現(xiàn)象的原因之一；在Y2 點(diǎn)處觀察到流速最大值，結(jié)合光學(xué)影像分析可知，冰川干流上游被大量的積雪覆蓋，夏季溫度升高易產(chǎn)生積雪性融水，由于上游地區(qū)海拔較高，巨大的高差轉(zhuǎn)換成融水的動能，使得該區(qū)域冰川流速達(dá)到峰值，該結(jié)果同前人的研究成果基本一致［37-38］。另外，受多種失相干因素的影響，冰川上游粒雪盆區(qū)域運(yùn)動分布場受噪聲影響較為嚴(yán)重。

3.3 冰川運(yùn)動結(jié)果精度分析

冰川區(qū)海拔較高，地勢險(xiǎn)峻，氣候條件惡劣，冰川流速野外實(shí)地驗(yàn)證工作難以有效開展，本文按照如下兩種方式對流速結(jié)果的精度展開分析。

3.3.1 非冰川區(qū)殘差統(tǒng)計(jì)分析

基于非冰川區(qū)在研究周期內(nèi)保持穩(wěn)定且未發(fā)生運(yùn)動的基本假設(shè)，本研究均勻選取各區(qū)域非冰川覆蓋區(qū)域進(jìn)行殘差統(tǒng)計(jì)，進(jìn)而對流速提取精度進(jìn)行評估，結(jié)果如圖7 所示。在研究周期內(nèi)冰川運(yùn)動總體精度約為0.5 m，方差接近于0，殘差結(jié)果整體呈正態(tài)分布，表明本研究的監(jiān)測結(jié)果具有較高的準(zhǔn)確性。

圖7 非冰川區(qū)殘余誤差統(tǒng)計(jì)Fig.7 Histograms of the residual error in the non-glacial region

3.3.2 同期國外產(chǎn)品對比

Landsat、Sentinel和ALOS等系列衛(wèi)星是亞洲高山區(qū)冰川運(yùn)動提取的主要數(shù)據(jù)源。其中，Landsat影像存檔數(shù)據(jù)多，覆蓋范圍廣且分辨率適中，大量研究表明該數(shù)據(jù)在山地冰川運(yùn)動提取研究中具有較強(qiáng)的適用性，是研究冰川運(yùn)動學(xué)的重要遙感數(shù)據(jù)源［35，37，39］，美國航空航天局（NASA）以其為主要數(shù)據(jù)源制作了青藏高原地區(qū)1985—2018 年逐年冰川流速產(chǎn)品數(shù)據(jù)集（https：//its-live. jpl. nasa. gov/），可為GF-3 影像提取的冰川運(yùn)動評估提供一定參考。2018 年產(chǎn)品數(shù)據(jù)集冰川流速與本文GF-3 冰川流速整體相同，空間分布趨于一致，僅在冰川支流區(qū)域及規(guī)模較小的區(qū)域吻合度較弱，該現(xiàn)象與所選研究周期和空間分辨率不完全一致密切相關(guān)。為更好地說明本文監(jiān)測結(jié)果的準(zhǔn)確性，選取準(zhǔn)同期的Landsat-8 數(shù)據(jù)對研究區(qū)冰川同步進(jìn)行流速提取工作，兩種不同數(shù)據(jù)源沿冰川主干剖線的運(yùn)動速度對比如圖8所示。各研究區(qū)不同數(shù)據(jù)源影像提取的冰川流速結(jié)果高度吻合，變化趨勢基本一致，兩者的離差遠(yuǎn)小于冰川實(shí)際運(yùn)動速率，表明本文監(jiān)測結(jié)果具有良好的精準(zhǔn)性，同時(shí)說明GF-3影像在冰川活動性監(jiān)測中具有較強(qiáng)的適用性。

圖8 多源遙感數(shù)據(jù)運(yùn)動速度對比圖：南以內(nèi)里切克冰川（a）；卻哥隆瑪冰川（b）；巴爾托洛冰川（c）；雅弄冰川（d）Fig.8 The curves of ice velocity along the longitudinal profile in Fig.1 derived from GF-3 and Landsat-8：South Inylchek glacier（a）；Chogolungma glacier（b）；Baltoro glacier（c）；Yanong glacier（d）

4 討論

4.1 冰川運(yùn)動影響因素分析

亞洲高山區(qū)幅員遼闊，地形錯綜復(fù)雜，氣候環(huán)境差異較大。本文所選的研究區(qū)遍及天山地區(qū)、喀喇昆侖地區(qū)以及橫斷山區(qū)，區(qū)域地理跨度大，在不同氣候因素的作用下冰川流速空間分布具有一定的差異性。此外，受氣候條件、冰川自身屬性及地形地貌等因素的影響，不同地理單元的冰川流速、消融特征等存在顯著性差異，相同區(qū)域不同規(guī)模和形態(tài)的冰川流速也存在明顯區(qū)別，與冰川長度、面積、基底滑移、地表坡度、物質(zhì)平衡等密切相關(guān)。

喀喇昆侖地區(qū)的卻哥隆瑪冰川及巴爾托洛冰川地理位置相近，但冰川運(yùn)動情況卻有所差異。由于數(shù)據(jù)覆蓋范圍受到限制，巴爾托洛冰川運(yùn)動圖所展示出的并不是全部主干，僅是整條干流的中下游，為此我們選用幅寬更大的Landsat-8影像同步進(jìn)行準(zhǔn)同期流速提取工作，結(jié)果表明主干區(qū)域流速分布在7～42 cm·d-1之間，與卻哥隆瑪冰川流速相比，在上游區(qū)域差異性較大。卻哥隆瑪冰川區(qū)周圍的山體起伏較大，周圍海拔較高的支流冰川匯入干流，加速了主干區(qū)域冰川的流速，流速特征改變明顯。此外，冰川坡度、坡向的不同也是影響流速分布的主要因素［40］。文中獲取了天山地區(qū)南以內(nèi)里切克冰川及其周圍的北以內(nèi)里切克冰川（North Inylchek）、凱迪冰川（Kaindy）和托木爾冰川（Tomur）流速分布，對比可知，規(guī)模較小的凱迪冰川和托木爾冰川流速較小，南以內(nèi)里切克冰川各支流冰川流速低于主干冰川，而北以內(nèi)里切克冰川則處于近乎停滯狀態(tài)。上述結(jié)果表明，同一區(qū)域內(nèi)部不同冰川流速具有差異性，一般而言，大型冰川流速比小型冰川快，小型冰川對于氣候的敏感程度較高。此外，冰川的屬性、地形走勢及高差變化等都會對冰川運(yùn)動分布產(chǎn)生影響［41-42］。

4.2 GF-3影像適用性分析

為充分凸顯我國自主研制的遙感衛(wèi)星數(shù)據(jù)在冰川活動性研究中的優(yōu)勢及應(yīng)用價(jià)值，選取兩處典型區(qū)域?qū)F-3 和Landsat-8 影像獲取的流速結(jié)果進(jìn)行比較（圖9），圖9（a）、9（b）為南以內(nèi)里切克南側(cè)支流冰川及托木爾冰川，圖9（c）、9（d）為巴爾托洛冰川部分區(qū)域。對比顯示，不同遙感數(shù)據(jù)源獲取的相同冰川區(qū)域運(yùn)動分布在空間上高度一致，但細(xì)節(jié)方面有所不同（圖中圈選區(qū)域）。高空間分辨率的GF-3 數(shù)據(jù)提取的冰川流速分布場在細(xì)節(jié)方面更為完整，較好地兼顧了小尺度冰川運(yùn)動信息，有利于冰川躍動研究；而15 m空間分辨率的Landsat-8數(shù)據(jù)僅能獲取流速較快、規(guī)模較大的冰川主干運(yùn)動分布信息，無法完整的反映整個冰川區(qū)運(yùn)動情況。此外，基于GF-3影像獲取的運(yùn)動分布場連續(xù)性更好，冰川支流匯入后對主干速度變化體現(xiàn)更加清晰，更加符合真實(shí)情況。上述對比表明分辨率更高的GF-3 影像在冰川運(yùn)動監(jiān)測中具有顯著的優(yōu)勢，為基于我國遙感衛(wèi)星數(shù)據(jù)的冰川運(yùn)動研究提供了可靠的支撐。

圖9 不同數(shù)據(jù)源流速提取對比分析：GF-3獲取冰川流速結(jié)果（a）；Landsat-8獲取流速結(jié)果（b）；GF-3獲取冰川流速結(jié)果（c）；Landsat-8獲取流速結(jié)果（d）Fig.9 The estimated average glacier velocity field from GF-3 and Landsat-8：the velocity extracted by GF-3（a）；the velocity extracted by Landsat-8（b）；the velocity extracted by GF-3（c）；the velocity extracted by Landsat-8（d）

5 結(jié)論

基于改進(jìn)的偏移量跟蹤技術(shù)和多景GF-3影像，高效地提取了2019—2020 年間天山地區(qū)南以內(nèi)里切克冰川、喀喇昆侖地區(qū)卻哥隆瑪冰川及巴爾托洛冰川、橫斷山區(qū)雅弄冰川的表面運(yùn)動特征，結(jié)果表明：

（1）GPU 并行化偏移量跟蹤算法能有效地提高計(jì)算效率，通過彌補(bǔ)全局性位移和地形效應(yīng)校正以及濾波降噪處理后，冰川表面流速結(jié)果的精確性得以大幅提升，冰川區(qū)運(yùn)動分布形態(tài)更加直觀。

（2）對比GF-3 影像及普適性較強(qiáng)的準(zhǔn)同期Landsat-8 影像提取的冰川流速結(jié)果顯示，二者的吻合程度較高，同時(shí)基于非冰川區(qū)未發(fā)生相關(guān)位移的假設(shè)開展殘差統(tǒng)計(jì)分析，精度可達(dá)0.5 m，驗(yàn)證了本文結(jié)果的可靠性。

（3）GF-3 數(shù)據(jù)能夠完整地提取冰川表面運(yùn)動分布，憑借其高分辨率的特點(diǎn)在小尺度冰川運(yùn)動獲取過程中具有較大的優(yōu)勢，為小型冰川運(yùn)動提取與分析以及冰川躍動研究提供了良好的數(shù)據(jù)支撐。

（4）受冰川自身屬性以及外部地形地貌、氣候等因素的影響，亞洲高山區(qū)內(nèi)部各區(qū)域冰川活動性存在顯著性差異，后期將繼續(xù)搜集實(shí)測數(shù)據(jù)和冰川區(qū)氣象資料，多維多層次分析冰川的運(yùn)動規(guī)律和氣候變化之間內(nèi)在的響應(yīng)機(jī)制，同時(shí)繼續(xù)完善不同地理單元典型區(qū)域冰川的運(yùn)動速度提取，為廣域青藏高原地區(qū)冰川活動性研究奠定基礎(chǔ)。

謹(jǐn)以此文紀(jì)念在慕士塔格科考中失蹤的中國科學(xué)院青藏高原研究所林樹標(biāo)博士。

致謝：感謝中國科學(xué)院空天信息創(chuàng)新研究院（AIR）提供的GF-3 數(shù)據(jù)，國家氣象科學(xué)數(shù)據(jù)中心（CMA）提供的研究區(qū)地面氣象數(shù)據(jù)，美國航空航天局（NASA）提供的ASTER 30 m DEM 數(shù)據(jù)，美國地質(zhì)勘探局（USGS）提供的Landsat-8 影像數(shù)據(jù)，中國科學(xué)院西北生態(tài)環(huán)境資源研究院提供的中國第二次冰川編目數(shù)據(jù)集。