基于Landsat-8時間序列影像分析西昆侖山地區冰川滑移特征

鄭　茜　孫建寶　張　永

1　中國地震局地質研究所地震動力學國家重點實驗室,北京市華嚴里甲1號,100029

基于Landsat-8時間序列影像分析西昆侖山地區冰川滑移特征

鄭茜1孫建寶1張永1

1中國地震局地質研究所地震動力學國家重點實驗室,北京市華嚴里甲1號,100029

摘要：利用圖像亞像元互相關分析方法處理了Landsat-8衛星獲取的時間序列影像數據，得到中國青藏高原西北部地區西昆侖峰區冰川勻速滑移的時空演化過程。利用亞像元影像互相關技術對Landsat-8光學影像精確配準，配準精度達到0.01像元，即該光學影像的水平形變監測精度達到0.15 m。通過對2013-07～2014-08的15景Landsat-8影像進行互相關和形變時間序列反演分析，獲得了西昆侖峰區兩條冰川的滑動位移場和速度場。研究表明，該區域的冰川基本處于勻速滑移狀態(無明顯加速和減速現象)；同時也驗證了Landsat-8光學影像在監測較大地表位移和地殼形變事件(如沙丘移動、地震、滑坡、火山等)上的應用潛力。

關鍵詞：西昆侖峰區；Landsat-8；配準；冰川滑移；亞像元圖像互相關

依據中國冰川編目，西昆侖山共有冰川5 485條，冰川面積8 817.78 km2，約占整個昆侖山冰川總面積的3/4以上[1]。來自昆侖山冰川的冰川水是塔里木盆地的主要水源，為生活在這一流域的近630萬人口提供飲用水[2-3]。分析研究這些高山冰川的位移場和它們的動態時空演化過程，對了解青藏高原地區的氣候變化、冰川活動對人類活動的影響等有重要意義。

目前，常用的形變監測技術包括差分全球定位系統技術(DGPS)、干涉合成孔徑雷達技術 (InSAR)和亞像元光學影像互相關技術等。前兩種方法雖然可以獲得較高的形變觀測精度，但在監測高山冰川方面具有較大的局限性[4]。在過去的幾十年里，美國的Landsat-1～7系列衛星所獲得的數據廣泛用于地球科學研究和人類日常生活的各個領域。特別是在冰川研究領域，通過分析數10 a尺度的Landsat數據，可以直觀地測量冰川的演進速度。然而，由于軌道控制方面不夠精準，應用Landsat-1～7系列衛星數據需要手工選擇地面控制點，以實現精確的地面定位，形變的量算也只能通過肉眼對圖像進行比較獲得。Landsat-8解決了上述問題，為冰川科學研究提供了非常有價值的數據源。

在本研究中，我們通過亞像元級的圖像匹配技術對軌道和幀號為145/35的15景Landsat-8影像數據進行互相關計算，并對其組成的時間序列進行SVD反演，獲得了西昆侖地區目標冰川的二維位移場及速度場。

1數據源

研究區為西昆侖峰區,位于青藏高原西北部的新疆自治區境內，是我國冰川分布最密集的區域之一(圖1(a))。

圖1　研究區域Fig.1　The studied area

Landsat-8圖像覆蓋范圍為185 km×185 km的區域，因此一景圖像可以覆蓋多條冰川。在一次圖像互相關形變計算中，可以得到圖像內所有冰川的水平位移場，因而具有很高的空間監測效率。同時，16 d的固定重訪周期和不間斷的數據獲取，使得該數據的時間分辨率比其他衛星數據源具有更大優勢。圖1(b)為2013-09-11獲取的Landsat-8影像；圖1(c)顯示了本研究中的目標冰川分布，即圖1(b)中矩形框所示的范圍。本研究將分析該景圖像范圍內的形變信息。

我們從USGS網站下載了2013-07～2014-08共15景覆蓋西昆侖峰區的Landsat-8影像，構成形變研究的時間序列，數據信息如表1所示。為了提高估算冰川速度場的精度，我們只提取空間分辨率為15 m的全色波段圖像進行圖像互相關計算。同時，選擇云量在30%以下的影像數據，以確保圖像可以有效地精確配準。

表1　研究區 Landsat-8影像數據

2數據分析方法

從USGS網站下載的Landsat-8影像為使用SRTM DEM數據正射糾正過的圖像。只有經過正射校正的圖像，排除成像幾何畸變，才能互相進行精確配準并分析形變。

圖像亞像元互相關分析技術，最先用于SAR圖像幅值信息的分析[5]，一方面用于計算SAR斜距向和方位向的位移場，另一方面用于SAR圖像的精確配準。其基本原理是在兩幅幅值圖像上，按照一定間隔選取一定的配準網格點，在每個網格點附近按照一定的窗口大小搜索另一幅圖像的相應位置，直到找到窗口內互相關系數最大的點，認為這兩個圖像互相關系數最大的點是地面同一個散射體目標。通過計算同一個目標在不同時間獲取的SAR圖像上的位置變化，就能知道該點地面發生的位移大小。為了實現高精度測量，SAR圖像的計算都是在單視復圖像上、雷達坐標系下進行的。而在光學圖像的分析中，圖像經過了地理編碼和正射校正，所有形變計算直接在地理坐標系下進行。采用同樣的方法，通過比較兩幅圖像的地面反射信息，可以實現亞像元級的圖像精確匹配和位移場計算。

使用SAR圖像配準方法，將兩幅光學圖像配準。一般配準精度(或稱配準誤差、偏差)在0.001～0.01像元范圍內，可以達到理想的形變監測信噪比。在圖像精確配準過程中，有必要控制獲取影像的云量。 在配準之后，把一幅圖像(稱為“輔圖像”或者“從圖像”)重采樣到另一幅圖像(稱為“主圖像”)的幾何坐標下[6]。圖像重采樣采用SINC插值方法。將任意兩幅配準和重采樣后的圖像在COSI-CORR軟件[7]中作互相關計算，可得到初步的位移結果。

2.1圖像配準

配準步驟是整個形變計算的關鍵。使用圖像亞像元互相關的方法，尋找兩幅圖像中共同的地面目標，用最小二乘的方法建立主從圖像的多項式變換關系。配準和重采樣過程消除了兩幅圖像的投影偏差，因而能夠保證獲得精確的地表形變信息[4]。另外一種配準方法是在主從圖像之間建立映射關系，即查找表。查找表記錄了由輔圖像幾何坐標到主圖像的映射信息。建立查找表的同時，考慮地形的影響，可以進一步優化圖像的配準過程。在建立了主從圖像的查找表之后，繼續用前面的配準多項式來校正查找表，并對其進行更新。基于更新后的查找表，重采樣輔圖像到參考圖像坐標系下。經過該配準步驟，主從圖像的配準精度可以進一步提升，配準偏差更小，因而更有利于形變的監測。在本研究中，我們選用獲取時間較早的圖像作為主圖像來配準獲取時間較晚的另一幅，當這幅圖像精確配準之后，也可將它作為一幅新的主圖像來配準其他從圖像。以此類推，配準所有獲得的Landsat-8數據。總的來說，較短時間間隔內的圖像對，其配準結果相對較好，信噪比較高。

2.2圖像互相關計算

地表水平位移場的獲取是通過多時相正射影像圖的像元互相關計算得到的。本研究采用的圖像互相關方法是在傅里葉域里對圖像進行反復迭代、無偏估計的過程[8]。這一過程已經在COSI-CORR 軟件中實現，它將生成兩幅互相關形變圖像，分別是地表水平位移場的東西向和南北向分量；同時給出一幅評估形變監測結果精度的信噪比圖像(圖2)。圖2給出的形變監測結果中，EW向形變不明顯，而NS向位移場中可以看到幾條明顯的白色高值區域，顯示幾條最活躍冰川形成的地面位移場。由于這些冰川主要為近南北走向，與地形梯度方向一致，因而冰川位移也以南北向為主，其NS向的形變比較顯著，這與本次觀測的結果完全一致。整個圖像互相關的信噪比SNR比較高，圖像顯示為高亮色。

圖2　主輔圖像互相關結果Fig.2　Cross-correlation products from a pair of master and slave images

2.3后處理

2.3.1條帶去除

目前，光學圖像一般通過多傳感器機械掃描得到，這是為了避免單一光電傳感器在掃描速率上的局限性[9]。例如，Landsat所使用的傳感器系統，每一個頻譜波段數據都由6臺傳感器同時掃描得到。因此，一次單相掃描能夠產生圖像上的6行，當下一次掃描時，衛星軌道適量前進，之前的6臺傳感器掃描生成下一個6行圖像(LANDSAT 用戶手冊)。然而，這些傳感器并沒有統一的增益函數，因此導致互相關結果圖上有明顯的平行軌道的條帶，并在垂直軌道方向起伏波動。例如，在ASTER和SPOT影像上出現過的類似振蕩偏差或CCD陣列的錯誤振幅[10-11]。這種因傳感器偏差導致的條帶誤差，可以通過COSI-CORR軟件的后處理工具Destrip移除。

2.3.2數據濾波和異常值消除

互相關計算得到的形變監測結果必須通過濾波來降低噪聲，但某些合理可信的結果常會被一些低相關系數的斑塊所掩蓋，必須先從形變場中移除一些過大或過小的異常值。首先計算位移值的概率密度分布，剔除低信噪比異常值，如圖3所示(藍色標識的部分認為是異常值，予以剔除)。

圖3　位移值概率分布Fig.3　The probability distribution of a displacement field

我們選用基于非局部均值算法的濾波方法[12]對結果圖像進行濾波處理，其在保證局部細節的同時降低了加性高斯白噪聲。非局部均值算法是通過逐步調試來實現數據的濾波，具體方法是：對每個像元點，將其周圍的小區域作為一個中心斑塊，斑塊大小一般為5×5或7×7像元。依據中心斑塊與其他斑塊的相似度來確定中心斑塊周圍更大范圍的搜索區域，最終以該搜索區域內所有像元的平均值作為中心像元的濾波值。

3觀測結果分析

3.1西昆侖峰區冰川位移場

在西昆侖山地區，大部分高山冰川呈南北走向，因此在互相關結果中，東西分量的水平位移信號相對較弱(圖2)。所以，我們以南北向水平位移結果為例，研究西昆侖峰區冰川滑移及其隨時間的演化。圖4(a)給出了西昆侖峰區的地表南北向水平位移場，這是由2013-09-27和2013-12-16獲取的兩幅Landsat-8正射影像圖作互相關得到的，其中一幅覆蓋目標冰川的正射影像如圖4(b)所示。該地區的位移最大值出現在圖4(a)中的灰色矩形框內。在80 d內，兩條目標冰川的滑移速率達到了1.5 m/d。此外，該地區內其他冰川上稍小的滑移速率也比較明顯。圖上紅色表示冰川向南移動，藍色表示向北移動。我們將滑移速率最大的兩條目標冰川作為研究對象，其位移場在圖4(c)中顯示。

圖4　西昆侖峰區的地表南北向水平位移場Fig.4　Horizontal displacement field of the glaciers in the west Kunlun mountains

圖5　目標冰川位移場的時間序列演化Fig.5　The time-series plot of the displacement fields of the target glaciers

3.2 時間序列分析

我們處理分析了2013-07～2014-08的15景Landsat-8影像(表1)，并且利用奇異值分解(SVD)算法計算了目標冰川的年地表形變[13]。圖5顯示西昆侖山地區的冰川在觀測時間段內，滑移量逐漸累積(南北向)。在384 d內，冰川位移隨時間呈明顯遞增趨勢。

為了確定目標冰川精確的滑移速率，我們在兩條冰川上選取6個參考點(圖4(c))，分析它們的位移(圖6)和速率(圖7)差異。由于圖像處理過程中存在一些不可避免的隨機噪聲，導致利用奇異值分解算法反演位移場時間序列的結果存在部分“負增長”值(圖6)，但依然可以明顯看到，在384 d內位移場隨時間近似呈線性遞增，同時曲線的斜率在不同的參考點上是不同的，說明各點的滑動速率并不完全一致。圖7給出了6個參考點的速率變化。可以看出，2013-11～2014-08冰川滑移速率保持在0.2～0.6 m/d，而對于每個參考點，冰川速度場呈勻速運動趨勢。因此，研究區的這兩條冰川始終保持了極為高速的滑動，其能夠被Landsat-8這種圖像分辨率的數據檢測到。與其相比，InSAR形變監測方法，由于時間去相干效應的影響，監測如此高速的地面冰川滑動存在很大的局限性，這也反映了Landsat-8數據的觀測優勢。

圖6　位移時間序列曲線Fig.6　The curves of the displacement time-series

圖7　速度時間序列曲線Fig.7　The curves of the velocity time-series

4結語

本研究首次利用Landsat-8光學圖像形變分析方法，獲得西昆侖峰區冰川的形變場時間序列演化特征。基于光學影像資料提取地表形變，應用時間序列分析的方法、原理，對Landsat-8光學影像進行形變分析，獲得西昆侖峰區冰川滑移的地表水平形變及其時間序列演化過程。結論如下：

1)光學影像互相關技術監測地表形變，特別是在監測較大形變和較快移動的地面目標時具有一定的優勢。

2)基于幅值信息的影像配準方法、亞像元影像互相關形變計算、SVD反演時間序列——這些方法能夠有效獲得研究區時空形變場，提高測量精度。

3)在384 d觀測期內，西昆侖峰區冰川位移場隨時間近似呈線性變化，滑移速率在0.2～0.6 m/d。

4)本研究的成果不僅對冰川動力學和氣候變化方面的研究有一定意義，而且驗證了Landsat-8光學影像在形變分析中的應用潛力。該方法也適用于分析大地震和滑坡等形變場。

致謝：USGS為本研究提供了所有的Landsat-8影像數據,在此表示感謝!

參考文獻

[1]焦克勤,姚檀棟,李世杰.西昆侖山32 Ka來的冰川與環境演變[J].冰川凍土,2000，22(3),250-256(Jiao Keqin,Yao Tandong,Li Shijie.Evolution of Glaciers and Environment in the West Kunlun Mountains During the Past 32 Ka[J].Journal of Glaciology and Geocryology,2000,22(3),250-256)

[2]李成秀,楊太保,田洪陣.1990～2011 年西昆侖峰區冰川變化的遙感監測[J].地理科學進展,2013,32(4)：548-559 (Li Chengxiu,Yang Taibao,Tian Hongzhen.Variation of West Kunlun Mountains Glacier during 1990-2011[J].Progress in Geography,2013,32(4):548-559)

[3]李成秀. 昆侖山冰川和積雪變化的遙感監測[D]. 蘭州:蘭州大學,2014(Li Chengxiu. Remote Sensing Monitoring of Glacier and Snow Cover Changes in the Kunlun Mountain[D]. Lanzhou :Lanzhou University,2014)

[4]Berthier E,Vadon H,Baratoux D,et al. Surface Motion of Mountain Glaciers Derived from Satellite Optical Imagery[J].Remote Sensing of Environment,2005,95(1):14-28

[5]劉方堅. 星載雷達影像配準方法的研究與實現[D]. 北京:中國地質大學(北京),2006(Liu Fangjian. The Registration Algorihtm Research of Spaceborne InSAR Images[D]. Beijing :China University of Geosciences(Beijing),2006)

[6]Massonnet D,Feigl K. Radar Interferometry and Its Application to Changes in the Earth’s Surface[J].Reviews of Geophysics,1998,36(4): 441-500

[7]Leprince S,Berthier E,Ayoub F,et al. Monitoring Earth Surface Dynamics with Optical Imagery[J].EOS Transactions,2008,89(1):1-12

[8]Scherler D,Leprince S,Strecker M R. Glacier-Surface Velocities in Alpine Terrain from Optical Satellite Imagery -Accuracy Improvement and Quality Assessment [J].Remote Sensing of Environment,2008,112(10):3 806-3 819

[9]Horn B K P, Woodham R J. Dest-Riping Landsat MSS Images by Histogram Modification[J].Computer Graphics and Image Processing,1979(10):69-83

[10] Leprince S,Muse P,Avouac J P. In-Flight CCD Distortion Calibration for Pushbroom Satellites Based on Subpixel Correlation[J]. IEEE Transactions on Geoscience and Remote Sensing,2008,46(9):2 675-2 683

[11] Teshima Y,Iwasaki A. Correction of Attitude Fluctuation of Terra Correction of Attitude Fluctuation of Terra with Parallax Observation[J].IEEE Transactions on Geoscience and Remote Sensing,2008,46(1):222-227

[12] Ayoub F,Leprince S,Keene L.User’s Guide to COSI-CORR[Z].Co-Registration of Optically Sensed Images and Correlation,2009

[13] Berardino P,Fornaro G,Lanari R,et al. A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms[J].IEEE Transactions on Geoscience and Remote Sensing,2002,40(11):2 375-2 383

Foundation support:National Natural Science Foundation of China,No. 41374040.

About the first author:ZHENG Qian,postgraduate,majors in time-series analysis of SAR interferometry and crustal deformation application,E-mail: 359300096@qq.com.

Fast and Uniformly Slipping Western-Kunlun Glaciers from Time-series Deformation Analysis Using Periodically Captured Landsat-8 Imagery

ZHENGQian1SUNJianbao1ZHANGYong1

1State Key Laboratory of Earthquake Dynamics,Institute of Geology,CEA,A1 Huayanli,Beijing 100029,China

Abstract:We processed time-series data,using a sub-pixel cross-correlation analysis method, to obtain the temporal evolution of the surface displacements of the west Kunlun mountain glaciers in the north-western Tibetan plateau. This work co-registers periodically acquired Landsat-8 images using the sub-pixel cross-correlation technique at an accuracy of 0.01 pixels, implying a horizontal displacements detecting precision of the optical image 0.15 m. By cross-correlation analysis of the Landsat-8 images acquired from July 2013 to August 2014 and time-series inversion, this study constructs the displacement and velocity fields of two glaciers in the west Kunlun mountains. It demonstrates that the glaciers in the study area are slipping fast and uniformly (no obvious acceleration or deceleration). Based on Landsat-8 images, the method of this study can potentially be used for measuring ground deformation caused by dunes shifting, earthquakes, landslides, or volcanoes etc.

Key words:west Kunlun mountains; Landsat-8; co-registration; glacial motion; sub-pixel cross-correlation

收稿日期：2016-03-16

第一作者簡介：鄭茜，碩士生，研究方向為InSAR時間序列分析與地殼形變,E-mail: 359300096@qq.com;zhengqianabc@sina.com。

DOI：10.14075/j.jgg.2016.07.010

文章編號：1671-5942(2016)07-0604-05

中圖分類號：P228

文獻標識碼：A

項目來源：國家自然科學基金(41374040)。