魏春蕊，楊成生，魏云杰，熊國華，李曉陽

（1. 長安大學地質工程與測繪學院，陜西 西安 710054；2. 中國地質環(huán)境監(jiān)測院（自然資源部地質災害技術指導中心），北京 100081）

0 引言

全球氣候變暖導致冰川和凍土消融，進而引發(fā)滑坡、泥石流等自然災害。西藏是我國山地冰川分布最廣泛的區(qū)域，據資料統(tǒng)計，20世紀以來西藏地區(qū)發(fā)生過多次冰湖潰決事件[1?2]。1991年由冰川活動誘發(fā)西藏八宿縣吉達鄉(xiāng)江查村各布馬錯冰湖潰決引發(fā)泥石流，沖毀10余處川藏公路路基、55 ha農田，沖走11戶民房、144頭牲畜、18座木橋[3]；2000年位于康馬縣薩馬達鄉(xiāng)雅拉山口北坡的沖巴嚇錯湖因后緣冰川滑坡引發(fā)涌浪潰壩洪水和泥石流，使冰湖下游的薩馬達、康馬、少崗和南尼等4個鄉(xiāng)受災[4]；2014年由于冰川泥石流導致波密縣的那隆藏布支溝冰湖潰決造成易貢鄉(xiāng)2.7ha農田和550 m鄉(xiāng)村道路被洪水沖毀等[5]。冰川運動變化的精確監(jiān)測能夠為冰川災害的預警和氣候分析提供理論依據，具有重大意義。

近年來，偏移量跟蹤技術在冰川位移的研究中得到了廣泛應用[6?7]。目前，偏移量跟蹤技術根據影像來源不同可分為兩大類：SAR偏移量技術、光學偏移量技術。由于SAR影像的獲取不受天氣條件的限制，且測量周期短，監(jiān)測范圍廣，許多國內外學者將SAR偏移量跟蹤技術應用到測量冰川流速的研究中[8?9]。然而SAR影像側視成像的特點使結果受疊掩和陰影等幾何畸變的影響較大，且計算結果為距離向和方位向形變，因此一些學者也集中研究基于光學影像的像素偏移量跟蹤技術[10?11]。例如，BERTHIER等[12]利用SPOT5光學影像獲取了阿爾卑斯地區(qū)山地冰川流速場[12]，DEHECQ等[13]利用Landsat-5和7影像，獲取了帕米爾-喀喇昆侖-喜馬拉雅(PKH)三年期間的冰川年速度場。

隨著多源數據和多平臺影像的發(fā)展，針對大量級的位移監(jiān)測，也有學者相繼提出多孔徑差分干涉測量技術(MAI)[14]、融合InSAR升降軌技術的冰川二維或三維運動監(jiān)測[15]，融合多個軌道的Offset-tracking或MAI形變測量技術[16?17]、小基線時序處理的PO-SBAS技術[18]和融合多源InSAR數據的三維形變測量[19]，但目前對于融合SAR影像與光學影像監(jiān)測多維度形變的研究較少。因此，文中選取西藏聶拉木縣希夏邦馬峰地區(qū)的大型冰川作為實驗區(qū)，基于方差分量估計，利用Sentinel-1與Landsat8兩種影像進行了冰川位移的三維分解研究，分析該方法在冰川運動監(jiān)測中的適用性和精確性。

1 研究區(qū)與數據

1.1 研究區(qū)概況

希夏邦馬峰為喜馬拉雅山脈中段的高峰之一，空間位置如圖1(a)所示，是一座完全在中國西藏自治區(qū)聶拉木縣境內的8 000 m級高峰。希夏邦瑪峰地區(qū)山地冰川持續(xù)退縮，并且冰湖加速擴張，該地區(qū)冰湖潰決等自然災害[20]增多。文中的研究對象屬于希夏邦馬峰的大型冰川，如圖1(b)所示，該冰川有3條上游支流并且在中游匯集，為冰川三維運動特征分析提供了典型實驗區(qū)域。

1.2 研究區(qū)數據

文中采用了Sentinel-1和Landsat8兩種數據源?？紤]到計算的高效率以及冰川在較長時間間隔里其位移量較明顯的因素，采用了從歐空局獲取的間隔48 d的2019年12月至2020年12月9景Sentinel-1升軌影像，詳細參數信息如表1所示，其距離向分辨率為5 m、方位向分辨率為20 m。對于光學數據，采用了2019年11月至2021年1月之間的6景美國Landsat8影像，參數信息如表2所示，Landsat8的全色影像空間分辨率為15 m。這兩種數據的影像覆蓋范圍如圖1所示。

表1 研究所采用的Sentinel-1數據參數Table 1 The Sentinel-1 data parameters used in the study

表2 研究所采用的Landsat8影像參數Table 2 The Landsat8 imaging parameters used in the study

圖1 實驗區(qū)數據覆蓋圖和冰川光學影像圖Fig.1 Data coverage map of the experimental area and glacier optical images

2 研究方法與數據處理

文中對獲取的Sentinel-1影像進行偏移量估計得到其距離向和方位向的形變量，對Landsat8影像進行偏移量估計得到東西向和南北向的形變量。為了融合這四個方向的形變量來求解冰川三維流速，操作步驟如下：(1)對距離向和方位向的形變位移結果進行投影轉換和分辨率重采樣，并將其分解至東西、南北和垂直向的位移；(2)根據觀測誤差的來源不同，與光學形變數據組成四個觀測方程，基于方差分量估計對計算的SAR偏移量數據和光學偏移量數據進行隨機模型的驗后估計，從而得到冰川的三維形變速率；(3)對于解算的結果采用同期的光學影像對進行冰川流速及其運動特征的對比分析，并通過計算穩(wěn)定區(qū)域的標準差進行精度估計。技術流程如圖2所示，主要分為三個關鍵的步驟：SAR偏移量估計、光學偏移量估計、聯(lián)合解算。

圖2 聯(lián)合解算冰川三維形變速率技術流程圖Fig.2 Technical flow chart of 3D deformation rate

2.1 SAR偏移量估計

由于SAR影像易受失相干的影響，文中采用小基線集像素跟蹤算法[21]，通過設置時空閾值得到最優(yōu)的基線組合（圖3），并利用對相干性要求較低的強度跟蹤算法進行偏移量的估算，其主要的參數設置包括配準窗口、互相關系數和步長。匹配窗口的大小對影像配準的精度非常重要，采用多尺度窗口可以提高計算結果的精度[22]，獲得更加精確的SAR影像距離向和方位向偏移結果。為了在偏移量信息和噪聲值之間達到較好的平衡，經過多次實驗確定最終的搜索窗口為128×128像素（距離向×方位向），搜索步長為5×1像素(距離向×方位向)，將相關系數閾值設置為0.2，以掩膜在失相干嚴重區(qū)域獲得的不可靠位移值。

圖3 SAR影像時空基線組合分布圖Fig.3 Space otemporal baseline combination distribution ofSAR images

2.2 光學偏移量估計

光學影像的偏移量估算結果易受失相關噪聲、軌道誤差、條帶誤差以及衛(wèi)星姿態(tài)角誤差的影響[23]，其中失相關噪聲包括由云雪、地表建筑物以及地形陰影等。為了減少失相關噪聲的影響，文中對獲取的時間跨度小且含云量少的影像計算每一影像對的太陽高度角和太陽方位角差值，選取其差值較小的影像對（表3），共計7組影像對用于影像對的匹配。

表3 研究所采用的Landsat8影像對Table 3 The Landsat8 image pairs used in the study

文中以COSI-Corr軟件為偏移量估算平臺，基于亞像素相關性匹配算法獲取影像對的東西向和南北向形變量。為了得到較好的數據處理結果，進行多次實驗選擇合適的參數：初始搜索窗口設置為64×64，最終相關搜索窗口設置為32×32，步長設置為4，迭代次數為2，掩膜閾值為0.9。處理過程中，采用了以下措施來削弱或消除誤差的影響：設置信噪比閾值>0.9去除失相關噪聲引起的誤差；利用一次多項式曲面擬合模型去除軌道誤差；利用“均值相減法”去除條帶誤差；利用改進后的“均值相減法”去除衛(wèi)星姿態(tài)角誤差[24]。

2.3 聯(lián)合解算

在多類觀測數據的聯(lián)合平差中，Helmert方差分量估計是對不同類型的觀測值進行定權的常用方法之一。它的基本思想：將觀測值按不同的觀測來源分類，確定各類觀測值的初始權，進行預平差；然后根據一定的原則計算各類觀測值的驗后單位權方差，進行觀測量方差的迭代運算，直至各類觀測值的單位權中誤差相等或各類單位權方差之比等于1。之前有學者已經用該方法探索過InSAR觀測量與GPS數據的融合[25]。同樣，文中基于方差分量估計對計算的SAR偏移量數據和光學偏移量數據進行隨機模型的驗后估計，具體步驟和原理：

首先，以Landsat8數據為參考對象，將SAR數據的坐標系統(tǒng)轉換到Landsat8數據的坐標系統(tǒng)下，并對SAR干涉對重采樣進行分辨率的統(tǒng)一，使得兩種圖像的像元能夠一一對應，便于后續(xù)的計算。

假設有N幅研究區(qū)的SAR影像組成了m個干涉對，對第r（r≤m）個影像對的相干點i來說，假設發(fā)生了勻速形變，則有：

式中：Losri、Azri——第r個干涉對得到的相干點i在Los向和方位向的形變量；

在東西、南北和垂直三個方向上的形變速率；

tr——第r個干涉對的時間間隔；

——Los向上相干點i在東西、南北和垂直方向上的投影；

——方位向上相干點i在東西、南北和垂直方向上的投影。其值根影的幾何關系（圖4）來計算。

圖4 升軌SAR影像的幾何關系(箭頭的方向為正)Fig.4 Geometric relation of rail SAR image (arrow direction is positive)

利用光學偏移量技術獲取了f幅影像對，對第k幅影像對的相干點i來說，則有：

基于最小二乘模型融合SAR干涉對觀測量和光學影像對觀測量求解三維地表形變：

式中：X——求解的三維形變速率，

L—2m個SAR偏移量測量值和2f個光學偏移量組成的觀測量，L=[Los1i···LosmiAz1i···

V——相對應的觀測殘差；

B——三個方向上的觀測值組成的矩陣。

假設觀測值的方差是已知的，則利用最小二乘平差可以計算得到三維形變速率的最優(yōu)估值：

其中，P為各個觀測量方差組成的權陣：

要得到最優(yōu)的三維形變速率估值，除函數模型外，還需要觀測量的隨機模型，也就是先驗方差，就可以在平差模型中精確定權。通常觀測值的方差往往很難精確獲取，所以基于方差分量估計進行權陣的驗后估計。由于SAR偏移量的觀測誤差受失相干、過采樣和匹配窗口的影響，光學偏移量的觀測誤差受各種噪聲、匹配窗口和步長等影響，因此根據觀測值所受觀測誤差的不同進行分組，將SAR觀測量的距離向和方位向各自分成一組，光學觀測量的東西向和南北向各自分為一組，總共四組觀測值。設每一組觀測值為Li、權重為Pi(i=1,2,3,4)，基于最小二乘方法得到第一次的估值：

假設初始權陣Pi為 單位陣，則方差分量和觀測殘差的關系：

最后計算新的權陣估計值：

3 結果與分析

3.1 冰川解算結果

3.1.1 SAR偏移量冰川位移結果

對于覆蓋研究實驗區(qū)冰川的9景Sentinel-1影像，基于小基線集的思想得到了14對SAR影像的優(yōu)化組合（圖3），利用SAR偏移量強度跟蹤技術獲取了目標冰川每對影像的距離向位移量（圖5）和方位向位移量（圖6）。從目標冰川的SAR偏移量組合對結果中可以得出，在2020年，1?3月整體形變量較小，4?8月的形變量增加，9月份形變量減緩，11?12月形變量較小，這也符合冰川季節(jié)性的變化趨勢。對于圖中形變值缺失的部分，由于SAR影像的失相干特性，在計算過程中掩膜了低相干區(qū)域的值。因此，形變值缺失的區(qū)域部分可以通過與光學偏移量的融合得到形變信息的補償。

圖5 SAR影像對之間的距離向冰川表面位移分布Fig.5 Distance between the SAR image pairs is displaced to the glacial surface

圖6 SAR影像對之間的方位向冰川表面位移分布Fig.6 Azimuthal to glacial surface displacement distribution between SAR image pairs

3.1.2 光學偏移量冰川位移結果

為了獲取目標冰川的光學偏移量觀測值，利用同時期的6景Landsat8影像，計算每一影像對的太陽高度角和太陽方位角差值，選取差值較小的影像對組合成7對影像用于偏移量計算?；贑OSI-Corr平臺的頻率域相關算法獲得初始偏移量解算結果，冰川的南北向和東西向表面位移分布分別如圖7和圖8所示。同樣，由于在誤差處理過程中，將SNR<0.9的形變值剔除而導致了局部形變值的缺失。因此，與SAR偏移量結果融合可以進行信息互補。

圖7 Landsat8影像對之間的南北向冰川表面位移分布Fig.7 Distribution of north-south glacial surface displacement between Landsat8 image pairs

圖8 Landsat8影像對之間的東西向冰川表面位移分布Fig.8 Distribution of east-west glacial surface displacement between Landsat8 image pairs

為了驗證太陽高度角和太陽方位角對結果精度的影響，選取光學偏移量計算結果中東西向和南北向穩(wěn)定區(qū)域（見圖1(b)的紅色框Roi1和Roi2）的標準差進行精度評定。通過統(tǒng)計發(fā)現標準差與太陽高度角差值和太陽方位角差值之間的相互關系如圖9所示，當太陽高度角差值和太陽方位角差值較小時，其標準差較低；當太陽高度角差值和太陽方位角差值較大時，其標準差較高。這也說明，在時間跨度較小的條件下，選擇太陽高度角差值和太陽方位角差值較小的光學影像對能夠提高偏移量估算結果的精度。

圖9 太陽高度角差值和太陽方位角差值與穩(wěn)定區(qū)域標準差的關系Fig.9 Relationship between the solar height angular difference and the solar azimuth difference and the standard deviation of the stable region

3.1.3 SAR和光學偏移量聯(lián)合解算冰川位移結果

按照文中介紹的方差分量估計聯(lián)合解算模型，將SAR與Landsat8兩種數據源提取結果進行融合，得到了冰川同期的三維形變速率（圖10）。結果顯示，該冰川在東西向的最大形變速率為21 cm/d，南北向的最大形變速率為68 cm/d，垂直向的最大形變速率為17 cm/d。總體上呈現出該冰川上游的運動速度快，中游的速度比較慢，下游的速度逐漸減小。對于冰川上游三個分支的冰川流速變化較大且不均勻，并且南北向的最大流速達到68 cm/d，主要原因有：地形起伏變化大，下坡運動增大了冰川的流速；上游為冰川消融區(qū)域，在監(jiān)測結果中存在特征點消失導致的局部信息不連續(xù)。

圖10 聯(lián)合解算的冰川表面運動速率分布Fig.10 Jointly solved distribution of glacier surface motion rates

為了更好地分析目標冰川的表面流速，將解算后的結果進行東西向、南北向和垂直向上的流速合成（圖11）?？傮w上該冰川從西南方向往東北方向流動，從上游開始消融到中游匯合，與冰川的實際運動方向具有一致性，其中流速較大的地區(qū)包括上游冰川的消融區(qū)域，坡度較大的區(qū)域和三條上游冰川支流匯合的區(qū)域。

圖11 合成后的冰川表面運動速率分布Fig.11 Distribution of glacial surface motion rates after synthesis

3.2 結果對比

為了檢驗聯(lián)合模型解算的冰川位移及其運動特征的可靠性，對聯(lián)合解算的南北向和東西向的冰川流速與同一時期光學影像對的冰川流速進行對比分析。這里選擇光學影像對的原因是，光學偏移量的結果與聯(lián)合解算后的南北和東西向為一致的方向，具有對比性，而SAR偏移量的結果為距離向和方位向，無可對比性。因此，選擇 2019年11月19日與2021年1月24日的Landsat8影像對，從表2可知該影像對的太陽高度角差值和太陽方位角差值較小，其差值分別為2.42°和7.4°，結合圖9分析得到的結果具有較高的精度。在進行誤差消除處理和去除低質量點后，對影像對進行了偏移量估計，得到該冰川南北向和東西向的表面運動速率。如圖12所示，從整體上看，冰川在南北向和東西向的運動速率達到了60 cm/d，并且運動速率較大的區(qū)域為冰川上游消融區(qū)和中游支流匯集區(qū)，與聯(lián)合結算后的結果與圖10(a)、圖10(b)較為符合。由于光學影像偏移量技術不能獲取垂直向位移，且沒有其他監(jiān)測數據可以被利用，因此我們無法對冰川的垂直向位移進行對比。

圖12 2019年11月—2021年1月光學影像對間的冰川運動速率分布Fig.12 Glacial motion rate distribution between optical image pairs from November 2019 to January 2021

為了進一步分析兩種數據監(jiān)測結果的一致性，按冰川的主流方向在上游和中游提取了這兩種監(jiān)測結果分別在南北向和東西向上冰川剖線MM′、NN′、PP′（圖1）的流速剖線圖（圖13），分別對這三條剖線進行分析：

（1）MM′剖線：在南北向上，整體趨勢具有一致性，同期光學影像對出現了因低質量點剔除的流速值見圖13(a)綠色框，但在數據融合后出現了較大的流速變化可達到60cm/d，根據地形因素分析，流速增大的位置對應于較陡的斜坡，受地形的影響較大；在東西向上，整體的變化趨勢較為一致，同樣在融合后變化較大的地方見圖13(b)綠色框，也是對應的斜坡位置。所以，融合后冰川流速的運動特征與地形起伏的變化相關。由于偏移量估算易受冰川表面特性變化的影響，對于同期影像對低質量的缺失值，通過融合后的數據得到了補償，從圖13中可以看出，融合后其對應的流速也較大。

圖13 冰川剖線MM′、NN′、PP′在南北向和東西向上的流速和高程圖Fig.13 Flow velocity and elevation of MM′, NN′, PP′

（2）NN′剖線：從圖13(c)和圖13(d)可以看出，融合后與同期影像對的結果整體趨勢上的變化是一致的，在南北向上，流速變化較為平緩；在東西向上，在坡度較陡的位置冰川流速較大。根據地形因素分析，冰川的運動與坡度的變化有關。

（3）PP′剖線：在冰川的中游區(qū)域，從圖13(e)和圖13(f)可以看出，總體趨勢具有一致性，并且地勢起伏變化小，同期影像對與融合后的結果沒有出現較大的流速變化差異，微小的差異是因為研究方法的不同。對于圖13(f)中東西向凸起的流速變化，是因為冰川的三條支流匯集增大了流速的變化。

3.3 精度分析

由于缺乏實測數據，采用先驗信息輔助進行精度分析是常用的方法，比如利用穩(wěn)定區(qū)域的標準差。由于光學偏移量的結果與聯(lián)合解算后的南北和東西方向一致，具有可對比性。因此，對于兩種不同方法的冰川流速監(jiān)測結果，我們統(tǒng)計了融合后與同期影像對在穩(wěn)定區(qū)域見圖1(b)的紅色框Roi1和Roi2的標準差（表4）。在南北方向上，Landsat8影像對在穩(wěn)定區(qū)域的統(tǒng)計結果分別為1.0 cm/d和0.77 cm/d，而融合解算結果在穩(wěn)定區(qū)域的標準差分別為0.56 cm/d和0.48 cm/d，相比于同期光學影像對結果提升了44%和38%；在東西方向上，Landsat8影像對在穩(wěn)定區(qū)域的統(tǒng)計結果分別為0.99 cm/d和0.84 cm/d，而融合解算結果在穩(wěn)定區(qū)域的標準差分別為0.91 cm/d和0.6 cm/d，相比于同期影像對結果提升了8%和29%?？傮w而言，融合后冰川位移監(jiān)測結果的不確定性小于同期光學影像對的監(jiān)測結果。

表4 融合后和同期影像對穩(wěn)定區(qū)域的標準差統(tǒng)計Table 4 Standard deviation statistics for st able regions after fusion and concurrent images

4 結論

多維、高精度的冰川運動變化監(jiān)測能夠為冰川災害的預警和氣候分析提供理論依據。因此，文中基于強度信息的SAR偏移量和基于亞像素頻率域相關性匹配的光學偏移量方法，結合方差分量估計，解算了覆蓋研究區(qū)冰川的三維流速，分析了冰川在不同方向的形變特征，并與同期光學影像對的計算結果進行對比分析，最后統(tǒng)計了穩(wěn)定區(qū)域的標準差進行精度分析。

（1）不同平臺影像偏移量融合可改善冰川位移監(jiān)測的時間分辨率，提供更多的冰川位移細節(jié)信息。

（2）實驗區(qū)冰川在融合后東西向的最大形變速率為21 cm/d，南北向的最大形變速率為68 cm/d，垂直向的最大形變速率為17 cm/d，其中流速較大的地區(qū)包括上游冰川的消融區(qū)域，坡度較大的區(qū)域和三條上游冰川支流匯合的區(qū)域。

（3）融合后的監(jiān)測結果與同期光學影像對在南北向和東西向上具有較好的一致性，并且山地冰川的運動與坡度有關，陡坡處流速增大。

（4）融合后冰川位移監(jiān)測結果的不確定性小于同期光學影像對的監(jiān)測結果，在南北方向和東西方向上監(jiān)測結果精度分別提升了38%和8%。

總體來說，同期單影像對解算方法雖然計算量小，但時間分辨率低，并且結果為二維形變。而不同平臺影像的融合解算具有獲取冰川的多維度形變信息，可進一步改善冰川監(jiān)測的時間分辨率和獲得更多形變特征的優(yōu)勢，但對于垂直向的研究及其精度分析需要進一步探索，因為垂直向的位移包含了冰川厚度的變化和沿斜坡的垂直向運動?；诓煌脚_解算多維度的形變信息研究可能是未來的發(fā)展趨勢，本文可為利用不同平臺的數據聯(lián)合監(jiān)測山地冰川的多維和高精度變化提供參考和技術支持。