王志偉，岳廣陽，吳曉東，張 文，王普昶，宋雪蓮，吳佳海

1. 中國科學(xué)院西北生態(tài)環(huán)境資源研究院，冰凍圈科學(xué)國家重點實驗室，青藏高原冰凍圈觀測研究站，甘肅 蘭州 730020 2. 貴州省農(nóng)業(yè)科學(xué)院草業(yè)研究所，貴州 貴陽 550006 3. Department of Geological Sciences, University of Texas at San Antonio, San Antonio 78249, USA

引 言

相比可見光、 紅外和近紅外光譜波段的遙感技術(shù)，具有全天時、 全天候優(yōu)勢的合成孔徑雷達(dá)干涉測量(interferometry synthetic aperture radar, InSAR)方法可以方便、 快捷的實現(xiàn)大范圍、 無接觸、 面狀的，精度為mm級的地表形變監(jiān)測。該技術(shù)方法已成熟地應(yīng)用于如地下水開采、 礦區(qū)沉降、 石油挖掘等人為因素引起[1]的，以及冰川退化、 凍土消融等全球氣候變暖等自然因素導(dǎo)致[2]的眾多地表變化研究中。特別是隨著SBAS-InSAR(small baseline subset InSAR)技術(shù)的發(fā)展，地表形變的反演精度和準(zhǔn)確度得到了進(jìn)一步提升，促使InSAR技術(shù)的應(yīng)用得到更加廣泛的推廣[3]。

多年凍土是冰凍圈的重要組成部分，指溫度能夠維持在零攝氏度以下狀態(tài)兩年及兩年以上的近地表土壤或巖石層。近百年來全球氣溫持續(xù)升高，多年凍土也開始逐步退化，眾多學(xué)者[4]指出這種變暖現(xiàn)象在21世紀(jì)可能仍將持續(xù)。活動層是多年凍土區(qū)夏季融化、 冬季凍結(jié)的土壤或巖石層，位于多年凍土之上。多年凍土退化會導(dǎo)致活動層厚度增加，從而改變土壤中的水、 熱環(huán)境。多年凍土的凍脹融沉過程不僅毀壞鐵路、 公路，也會破壞建筑工程結(jié)構(gòu)，甚至造成極大的生態(tài)、 地質(zhì)災(zāi)害。環(huán)境變暖多年凍土熱融效應(yīng)增加后，會促使大團(tuán)聚體破碎成小團(tuán)聚體，釋放大量有機(jī)碳、 硝態(tài)氮等物質(zhì)，進(jìn)而改變土壤生物和微生物環(huán)境，引起高寒生態(tài)環(huán)境的惡化。同時因氣候變暖土壤水、 熱狀況發(fā)生改變，同樣會使植被條件發(fā)生變化，進(jìn)而影響到地表的一系列特征，如反照率、 降水的滲透速度、 土壤中的蒸騰和蒸散、 以及土壤侵蝕等，從而打亂水文和氣候系統(tǒng)的循環(huán)速率。上述近地表過程的變化，會在大氣圈的下墊面底部對環(huán)境產(chǎn)生反作用，甚至影響多年凍土周邊地區(qū)乃至全球的氣候變化，進(jìn)而影響到生態(tài)系統(tǒng)的發(fā)生、 發(fā)展進(jìn)程。

新疆地處我國西北干旱區(qū)生態(tài)脆弱帶，是我國五大畜牧業(yè)生產(chǎn)基地之一，生態(tài)資源豐富。冰凍圈要素中的冰川、 凍土是區(qū)域內(nèi)水循環(huán)過程中不可或缺的重要因子。新疆地區(qū)的多年凍土分布可稱為“兩環(huán)”模式，多年凍土分布于阿爾泰山、 天山和昆侖山構(gòu)成的“兩環(huán)”脊線上，“兩環(huán)”內(nèi)部的準(zhǔn)格爾盆地和塔里木盆地為季節(jié)凍土區(qū)。其中，天山深居內(nèi)陸，分布于其上的多年凍土具有明顯的垂直地帶性[5]，為探討新疆多年凍土區(qū)地表形變提供了天然的研究場地。

目前，利用InSAR技術(shù)針對天山地區(qū)的研究多集中于對冰川的討論[2]，較少應(yīng)用到對多年凍土的分析中，本工作試圖利用SBAS-InSAR方法和39景ENVISAT影像(覆蓋時間從2003年6月17日到2010年6月15日)對天山多年凍土區(qū)地表變形規(guī)律做出部分探討，以期為后續(xù)的冰凍圈相關(guān)科學(xué)研究提供數(shù)據(jù)基礎(chǔ)和方法借鑒。

1 研究區(qū)概括

研究區(qū)位于天山山脈中段，位于北緯43.4°—44.7°，東經(jīng)83.9°—85.6°之間(圖1所示)，區(qū)域內(nèi)北部主要為平原地區(qū)，南部基本為山區(qū)。天山山脈深居大陸內(nèi)陸，屬于典型的大陸性山地氣候。根據(jù)早期研究可知[5]，天山山脈具有典型的凍土分布特征，區(qū)域內(nèi)高寒草甸和沼澤草甸分布廣泛。

圖1 研究區(qū)概括(a): 研究區(qū)地理位置；(b): 高程圖Fig.1 Location of study areas(a): Distribution of study area; (b): Digital elevation model

2003年至2010年期間，區(qū)域內(nèi)年均溫和年降水分別為(-2.7±8.27) ℃和(340±85) mL，如圖2所示(依據(jù)《中國區(qū)域高時空分辨率地面氣象要素驅(qū)動數(shù)據(jù)集》[6-7]計算獲取)，其基本分布規(guī)律為北部平原地區(qū)溫度高、 降水少，南部山地溫度低、 降水多，以上雨、 熱分布特點為多年凍土的發(fā)育提供了十分有利的條件。以3 000 m海拔為界，研究區(qū)的山地和平原地區(qū)2003年至2010年間的年均溫和年降水分別為(-10.87±2.98) ℃、 (420±38) mL和(0.54±7.34) ℃，(310±77) mL，具有明顯的差異。

圖2 2003年—2010年間多年平均溫度(℃)(a)和降水(mm·yr-1)(b)分布圖Fig.2 The spatial distribution of annual average temperature (℃) (a) and annual average precipitation (mm·yr-1) (b) from 2003 to 2010

2 地表形變反演研究算法

2.1 研究數(shù)據(jù)

如表1所示，本研究中涉及的ASAR (advanced synthetic aperture radar)數(shù)據(jù)主要有39景。經(jīng)初始處理[8]后，因基線和多普勒質(zhì)心差的原因，其中6景影像沒有參與配對。該ENVISAT ASAR數(shù)據(jù)集屬于單視復(fù)數(shù)影像(single look complex, SLC)，為C波段、 IS2模式遙感數(shù)據(jù)，入射角23°。研究數(shù)據(jù)獲取方式為降軌，數(shù)據(jù)來源于歐空局(European Space Agency, ESA)。方位向分辨率和距離向分辨率分別為22.6和4.05 m。本文中所用DEM資料為STRM V4版本數(shù)據(jù)，該數(shù)據(jù)在赤道的分辨率約為90 m[9]。《中國區(qū)域高時空分辨率地面氣象要素驅(qū)動數(shù)據(jù)集》[6-7]下載于寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心，覆蓋時間范圍為2003年1月1日到2010年12月31日，包括月、 日、 年三種時間尺度數(shù)據(jù)，其空間分辨率為0.1°。此外，文中所有地圖投影方式均為WGS84坐標(biāo)系。

表1 本研究中ENVISAT ASAR影像集的詳細(xì)信息(連接圖生成時主影像為20090106)Table 1 ENVISAT ASAR datasets products applied in this study (The image of 20090106 was used as the super master-imagine generating a connected graph)

2.2 算法原理介紹

星載InSAR技術(shù)在地表形變監(jiān)測領(lǐng)域應(yīng)用已經(jīng)十分成熟，其工作原理主要是集合成孔徑雷達(dá)技術(shù)與干涉測量技術(shù)于一體，計算距離信息的經(jīng)典干涉相位計算公式如式(1)[3]

φint=φflat+φtopo+φdef+φatmo+φnoise

(1)

式(1)中,φflat指平地相位，φtopo指地形相位，φdef指地表形變相位；φatmo指大氣延遲相位；φnoise指噪聲引起的相位。首次將該技術(shù)應(yīng)用于厘米級地表形變的研究方法稱為雷達(dá)差分干涉技術(shù)(D-InSAR)[10-11]。在多年凍土區(qū)的研究指出早期InSAR方法的主要限制因素是時空失相關(guān)和凍土地表的非線性形變[9, 12-13]，而SBAS-InSAR方法則因奇異值分解法(singular value decomposition, SVD)可以很好地克服上述問題[14]，SBAS-InSAR方法的具體原理如式(2)所示[15]。首先，生成差分干涉圖。干涉圖的數(shù)量如式(2)所示

(N+1)/2≤M≤N(N+1)/2

(2)

式(2)中，N+1景的ASAR研究區(qū)影像獲取于具體的研究時期，M對干涉圖則需滿足特定的時空基線閾值設(shè)置規(guī)則，假設(shè)每期ASAR影像都至少有1景影像與之干涉，則如式(2)，N+1景的ASAR影像生成的M對干涉圖數(shù)量應(yīng)在(N+1)/2至N(N+1)/2之間。時空基線的閾值設(shè)置規(guī)則一般通過多普勒質(zhì)心差來控制[13]。之后，通過Ta和Tb時刻獲取的ASAR影像來計算像元上的解纏相位，假設(shè)該像元在方位向和距離向的數(shù)值為(x,r)，則生成第i景的干涉圖公式為

d(Ta，x，r)]+Δφi，topo(x，r)+Δφi,orb+Δφi，res(x，r)

(3)

式(3)中，φ是干涉相位信息；i是干涉圖的序號，i∈(1, 2, …,M)；λ是傳輸信號的中心波長(本研究中ENVISAT ASAR數(shù)據(jù)的中心波長是5.67 cm)；d(Tb，x，r)和d(Ta，x，r)是雷達(dá)視線方向(line-of-sight，LOS)的積累量，Ta和Tb時刻主要相比T0時刻而言(T0時刻是參考影像的獲取時間)；Δφi,topo是來自于DEM誤差的相位信息；Δφi,orb是來自于衛(wèi)星軌道信息誤差的相位信息；Δφi,res是來自于殘余相位數(shù)據(jù)的相位信息。

假設(shè)vj是p到q時刻的平均相位速度，則相位和時間滿足如式(4)

vj×(Tp-Tq)=Δ(φp-φq)

(4)

式(4)中，p和q是獲取ASAR影像的具體時刻，φ是相位信息?？紤]到不同的時間分段，T0到Tj時刻的第j幅干涉圖的相位信息可以通過如式(5)計算得到

(5)

式(5)中，Δφj是影響不同時間間隔的相位速度積分信息，改寫成矩陣表達(dá)式為

Av=Δφ

(6)

式(6)中，A是系數(shù)矩陣；速度v是速度矢量，可以通過系數(shù)矩陣A計算獲取。因為在SBAS-InSAR處理過程中會有多個主影像，有可能導(dǎo)致系數(shù)矩陣A秩虧。通常利用SVD方法來處理，首先會生成一個逆矩陣，然后獲得速度矢量的最小范數(shù)解，最終得到各個時間段的形變量[14]。

2.3 數(shù)據(jù)處理流程

研究中的數(shù)據(jù)處理流程如圖3所示。其中生成連接圖時，選取空間基線距小于500 m和時間基線距小于550 d的ASAR影像，最終生成M幅(本研究中M為126)差分干涉圖，圖4顯示結(jié)果為生成連接圖時的時空基線圖。

圖3 地表形變反演流程Fig.3 Processing of surface deformation calculation

圖4 時空基線圖(a): 空間基線距；(b): 時間基線距；(c): 3D空間基線距黃色、 綠色和紅色點分別代表超級主影像、 從影像及未參與配對影像Fig.4 Plots of space and time baselines(a): Plot of position baseline; (b): Plot of time baseline distance; (c): Delaunay 3D plot of position baseline distance Blue lines represent interferometric pairs；Yellow, green and red diamonds denote super-master, slave and no-matched images, respectively

在生成連接圖后，配合DEM數(shù)據(jù)，經(jīng)過干涉圖去平、 自適應(yīng)濾波、 相干系數(shù)圖生成、 相位解纏后，再對M對干涉圖進(jìn)行編輯，保留處理較優(yōu)的結(jié)果，例如圖5(a，b)所示分別為20090421和20081202干涉的結(jié)果)，對處理較差的結(jié)果(例如圖6(a，b)所示分別為20050308和20060221干涉的結(jié)果)進(jìn)行剔除處理，此時移除的干涉圖數(shù)量即為L(本研究中剔除的L干涉圖為52對)。

圖5 較優(yōu)干涉結(jié)果(a): 較優(yōu)結(jié)果；(b): 較差結(jié)果Fig.5 High quality results of interferograms(a): High quality result; (b): Low quality result

圖6 較差干涉結(jié)果(a): 20050308; (b): 20060621Fig.6 Low quality results of interferograms(a): 20050308; (b): 20060621

之后利用地面控制點，進(jìn)行軌道精煉和重去平處理，然后估算形變速率和殘余地形，經(jīng)過相關(guān)系數(shù)閾值控制和SVD方法計算，再結(jié)合空間低通濾波和時間高通濾波方法，完成最終的研究區(qū)時間序列地表形變反演，如圖7所示。

3 反演結(jié)果與分析

研究區(qū)反演結(jié)果共包括33期形變影像，除第一景20040113(參考影像)的所有數(shù)值為0外，其余32期形變結(jié)果都是相對參考影像的數(shù)值[14]。本文以其中4期形變量影像為例，包括20041228(與20040113相比，季節(jié)接近)、 20050726(與20040113相比，季節(jié)相反)、 20090106(與20040113相比，季節(jié)接近)和20100615(與20040113相比，季節(jié)相反)，如圖7(a,b,c,d)所示。根據(jù)圖7可知，研究區(qū)內(nèi)反演質(zhì)量較佳區(qū)域為平原地區(qū)，山區(qū)的地表變形反演稀少而零散。不過如有需要，可通過各種空間插值方法來獲取完整的面狀數(shù)據(jù)，從而用以評估感興趣區(qū)域。為保證更準(zhǔn)確地探索地表形變與溫度、 降水的關(guān)系，不做插值處理。此外，僅通過從影像表現(xiàn)來看，20041228和20050726在山區(qū)西部表現(xiàn)出地表抬升狀況，而在20090106和20100615兩期影像中則轉(zhuǎn)變?yōu)槌两惮F(xiàn)象，四期影像在山區(qū)東部的區(qū)域基本都呈現(xiàn)出抬升的態(tài)勢。而在平原地區(qū)，地表的沉降和抬升在4期影像中變化不明顯，大部分表現(xiàn)為抬升，在接近城市的區(qū)域多表現(xiàn)為沉降。

圖7 研究區(qū)4期示例形變結(jié)果(a): 2004/12/28; (b): 2005/07/26; (c): 2009/01/06; (d): 2010/06/15Fig.7 Deformation results of 4 example imagines in the study area(a): 2004/12/28; (b): 2005/07/26; (c): 2009/01/06; (d): 2010/06/15

在研究區(qū)監(jiān)測范圍內(nèi)，2004年至2010年的地表形變速率如圖8所示。研究區(qū)內(nèi)平原地區(qū)的形變主要表現(xiàn)為抬升，僅北方接近城市的區(qū)域表現(xiàn)出強(qiáng)烈的沉降現(xiàn)象。研究區(qū)內(nèi)山區(qū)的形變結(jié)果較為零散，大體上在西部和東部地區(qū)分別呈現(xiàn)出沉降和抬升現(xiàn)象。盡管研究區(qū)范圍內(nèi)有不同程度的沉降和抬升，但研究期間內(nèi)的形變速率基本在±5 cm·yr-1之內(nèi)，平均形變速率為(-0.07±3.38) mm·yr-1。雖然整體變形不是十分明顯，但是均值體現(xiàn)出該區(qū)域內(nèi)地表存在輕微的沉降現(xiàn)象，同全球氣候變暖導(dǎo)致凍土退化的情形一致。此外，與北方西部湖泊附近沉降速率在-6和7 cm·yr-1之間相比，位于平原區(qū)北方中部的烏蘇市地表沉降速率基本在-10 cm·yr-1以上。山區(qū)的整體形變基調(diào)以沉降為主，符合全球變暖凍土開始退化的規(guī)律。山區(qū)每一個有具體數(shù)值的形變點滿足海拔大于3 000 m的提取條件，其具體分布位置如圖9所示。利用溫度和降水?dāng)?shù)據(jù)，分析這些形變點的年變化規(guī)律如圖10所示?？傮w和不同變形速率間隔的形變樣本點基本都呈現(xiàn)出逐步變暖的趨勢，而且同樣受2008年末和2009年初的極端降溫事情影響。不過除地表形變速率小于-2.0 cm·yr-1的樣點在2010年溫度仍然有所下降外，區(qū)域形變樣本點在2010年都有繼續(xù)變暖的趨勢。研究區(qū)的山區(qū)地表形變速率小于-2.0 cm·yr-1、 在-2.0到2.0 cm·yr-1之間和大于2.0 cm·yr-1的樣本點數(shù)量依次為6 364，6 449和2 385個。

圖8 平均形變速率Fig.8 Average deformation rates

圖9 山區(qū)形變樣點空間位置分布圖Fig.9 Deformation points in mountain regions of study areas

圖10 山區(qū)形變樣點年均溫度和降水曲線圖(a): 所有樣本點；(b): 小于2 cm·yr-1樣本點；(c): -2～2 cm·yr-1樣本點；(d): 大于2 cm·yr-1樣本點Fig.10 Plots of annual precipitations and annual mean temperatures in mountain regions(a): All points；(b): Ground deformation rate <2 cm·yr-1；(c): Ground deformation rate between -2～2 cm·yr-1；(d): Ground deformation rate >2 cm·yr-1

4 結(jié) 論

近年來，因人類活動的干擾和全球氣候變暖的日益嚴(yán)重，不同區(qū)域的多年凍土開始出現(xiàn)退化現(xiàn)象，退化嚴(yán)重的地區(qū)甚至?xí)l(fā)生強(qiáng)烈的水土流失，對草場、 環(huán)境、 建筑及工程設(shè)置造成威脅[3]。特別是對于天山地區(qū)，不僅存在“全球氣候變化指示器”的多年凍土，而且還是重要的畜牧業(yè)生產(chǎn)基地。無論是在環(huán)境保護(hù)，還是人類生產(chǎn)、 生活方面，都需要著重監(jiān)測。

下列的幾方面需要重點關(guān)注：

(1)地表形變反演結(jié)果從空間來分析，雖然在平原地區(qū)表現(xiàn)優(yōu)良，但是在山區(qū)卻較為稀少和零散。即使可以后期可以進(jìn)行插值處理，卻不一定能完全滿足一些需要高精度形變數(shù)據(jù)的研究和應(yīng)用。后期可以考慮利用其他軟件和不同InSAR方法進(jìn)行反演對比。

(2)從時間角度分析，難以反演不間斷的長時間序列擬合形變結(jié)果?，F(xiàn)有常用的InSAR數(shù)據(jù)，特別是開源數(shù)據(jù)反演結(jié)果基本上一年僅存在幾期。當(dāng)利用月尺度氣象或其他科學(xué)數(shù)據(jù)集對比分析時(如本文中使用的溫度和降水?dāng)?shù)據(jù))，難以實現(xiàn)一對一匹配。如果試圖利用日或者小時尺度的數(shù)據(jù)集進(jìn)行分析時，則存在更大的科學(xué)挑戰(zhàn)，因此后期地表形變的研究需要重點關(guān)注更細(xì)時間尺度的長時間序列地表形變結(jié)果反演。

(3)凍土變化的時間滯后性也需著重考慮?，F(xiàn)有研究對時間滯后性的分析不是十分充足，但是凍土發(fā)生變化處于地表內(nèi)部，自然因素對其的影響必然存在一個積累后再變化的過程。后期針對凍土變化的研究分析需要越來越多的關(guān)注時間滯后性對其影響，衡量溫度或者降水等要素對其時間滯后尺度和周期的具體數(shù)值，正是一個亟需科學(xué)研究揭示的關(guān)鍵點。