基于D-InSAR 技術(shù)的伊朗巴姆地震地表形變監(jiān)測(cè)1

寇 程 柯長(zhǎng)青

（南京大學(xué)地理信息科學(xué)系，南京 210046）

引言

傳統(tǒng)的地表形變監(jiān)測(cè)方法一般有水準(zhǔn)測(cè)量、三角測(cè)量、三邊測(cè)距（查顯杰，2007）、GPS測(cè)量（劉大杰等，1999）等，但是這些測(cè)量方法都具有范圍小、難度大、周期長(zhǎng)、成本高的缺陷。合成孔徑雷達(dá)干涉測(cè)量（Interferometric Synthetic Aperture Radar, InSAR）技術(shù)是以合成孔徑雷達(dá)（Synthetic Aperture Radar, SAR）復(fù)數(shù)據(jù)提取的相位信息為數(shù)據(jù)源獲取地表三維信息和變化信息的一項(xiàng)技術(shù)（王超等，2002），差分合成孔徑雷達(dá)干涉測(cè)量（D-InSAR）技術(shù)是InSAR技術(shù)的延伸，差分是指在研究區(qū)形變前后的SAR影像產(chǎn)生的干涉圖中去除地球曲率、地形因素等產(chǎn)生的相位，從而得到地表形變信息。它是迄今為止獨(dú)一無(wú)二的基于面觀測(cè)的形變遙感監(jiān)測(cè)手段（劉國(guó)祥，2006）。與常規(guī)方法相比，D-InSAR技術(shù)監(jiān)測(cè)地表形變具有大面積、連續(xù)、快速、準(zhǔn)確的優(yōu)勢(shì)，同時(shí)觀測(cè)無(wú)需布設(shè)地面控制點(diǎn)，從而降低了對(duì)危險(xiǎn)地區(qū)（比如火山地區(qū)）觀測(cè)的危險(xiǎn)性（馬超等，2004）。

國(guó)內(nèi)外學(xué)者對(duì)D-InSAR技術(shù)監(jiān)測(cè)地表形變監(jiān)測(cè)進(jìn)行了廣泛研究。D-InSAR技術(shù)被用于臺(tái)灣集集地區(qū)地震（劉國(guó)祥等，2002）、昆侖山口西地震（單新建等，2004）、張北地震（張紅等，2002）、汶川地震（孫建寶等，2008）、Landers地震（Zebker等，1994）和Hyogoken-Nanbu地震（Ito等，2003）同震地表位移監(jiān)測(cè)和同震形變場(chǎng)參數(shù)提取。這些研究證明了 D-InSAR技術(shù)在地表形變監(jiān)測(cè)方面的可行性，其監(jiān)測(cè)精度可達(dá)到厘米級(jí)。

2003年12月26日，伊朗巴姆地震發(fā)生后，許多研究者對(duì)巴姆地震進(jìn)行了研究（季靈運(yùn)等，2009；凌勇等，2006；羅揚(yáng)等，2006；王志勇等，2008），得到了巴姆地震LOS方向上的同震形變場(chǎng)。本文采用二軌法利用 ENVI SARscape模塊，基于伊朗巴姆地震前后的兩景SAR影像和SRTM DEM數(shù)據(jù)進(jìn)行差分干涉處理，并利用雷達(dá)波束的局部入射角得到垂直方向上的同震形變場(chǎng)。最后結(jié)合GIS三維空間分析技術(shù)分析了地震造成的形變特征，并探討了去相關(guān)效應(yīng)產(chǎn)生的原因。

1 研究區(qū)地理背景及數(shù)據(jù)源

巴姆市位于伊朗克爾曼省首府克爾曼市東南 195km處，地處盧特沙漠的邊緣，海拔1000m左右（張保淑，2009），地理位置見(jiàn)圖1。這里干燥少雨，多風(fēng)（Bam (Iran), 2003）。2003年12月26日，世界時(shí)01:56:52，巴姆地區(qū)發(fā)生了里氏6.7級(jí)地震，震中位于29.01°N，58.26°E（李毅等，2008），距巴姆城市約15km。巴姆古城在這次地震中幾乎毀于一旦（舒寧，2003）。

本文采用ENVISAT衛(wèi)星C波段的兩景ASAR Image SLC數(shù)據(jù)，數(shù)據(jù)的部分參數(shù)見(jiàn)表1。利用重采樣后的SRTM DEM數(shù)據(jù)進(jìn)行地形相位的去除，格網(wǎng)大小為25m×25m，UTM投影40分帶，參考橢球體采用WGS84。研究區(qū)SAR影像見(jiàn)圖1。從圖1可以看出，數(shù)據(jù)覆蓋范圍大部分地方地形平坦，西南角有部分覆蓋了山區(qū)。

表1 ENVISAT ASAR數(shù)據(jù)部分參數(shù)Table 1 Some parameters of ENVISAT ASAR data

2 D-InSAR技術(shù)原理

Zebker等（1994）、舒寧（2003）、王超等（2002）對(duì)D-InSAR技術(shù)測(cè)量形變的基本原理已經(jīng)有了非常詳細(xì)論述，本文不再贅述。能夠進(jìn)行干涉處理的 SAR影像必須處理成 SLC（Single Look Complex）格式，SLC影像上的每個(gè)象元可表示為：

式中，A (m, n) 是影像的強(qiáng)度信息即振幅；φ (m, n) 是相位信息。

干涉處理就是將兩景配準(zhǔn)后的影像逐象元進(jìn)行復(fù)共軛相乘，記兩景影像象元分別是（何秀鳳等，2012)：

逐象元進(jìn)行復(fù)共軛相乘后得到：

因此，干涉相位為：

這也是最后干涉圖上每個(gè)象元的干涉相位信息，φ (m, n) 的取值范圍為 (-π, π)。但是，地球表面是一個(gè)具有一定曲率的球面，微波信號(hào)在大氣中傳播時(shí)還會(huì)受到大氣的影響，而且雷達(dá)影像有一定的噪音，所以一般認(rèn)為干涉相位由以下幾部分組成（劉國(guó)祥等，2002）：

式中，φf(shuō)lat_Earth是地球曲率產(chǎn)生的相位信息，也就是高度不變的平地在干涉圖中表示出來(lái)的干涉條紋隨距離向和方位向的變化而呈周期性變化的現(xiàn)象（王超等，2002）；φtopo是地形產(chǎn)生的相位信息，是由于地形起伏造成的，可用于建立DEM；φdisp是由地表形變產(chǎn)生的相位信息，反映了地表形變信息；φpath是雷達(dá)信號(hào)在傳播過(guò)程中受到大氣影響產(chǎn)生的相位信息；φnoise是雷達(dá)影像噪聲產(chǎn)生的相位信息；n·2π是相位的整周數(shù)，需要通過(guò)相位解纏獲得。

3 實(shí)驗(yàn)驗(yàn)證

3.1 實(shí)驗(yàn)步驟

采用D-InSAR技術(shù)進(jìn)行地表形變監(jiān)測(cè)主要包括以下幾個(gè)步驟：基線估計(jì)、影像配準(zhǔn)、主干涉圖生成、平地相位及地形相位的消除和相位解纏形成形變圖。數(shù)據(jù)處理流程見(jiàn)圖 2（邵蕓等，2002）。

圖2 二軌法差分干涉測(cè)量流程圖（邵蕓等，2002）Fig. 2 The data processing procedure for two-pass differential radar interferometry (from Shao et al., 2002)

①基線估計(jì)。干涉像對(duì)的基線大小直接影響著干涉圖的質(zhì)量，基線過(guò)大則無(wú)法產(chǎn)生干涉。因此在進(jìn)行干涉處理前，首先要進(jìn)行基線估計(jì)，以確定像對(duì)能否產(chǎn)生干涉。一般認(rèn)為ENVISAT的C波段數(shù)據(jù)臨界基線大約為1100m，最佳基線應(yīng)小于200m，而用于地面變化監(jiān)測(cè)的基線大小應(yīng)小于 70m（舒寧，2003），本文所用數(shù)據(jù)的垂直基線估計(jì)結(jié)果為 9.1m，可以產(chǎn)生干涉。②影像配準(zhǔn)。進(jìn)行干涉的兩幅SAR影像之間必須進(jìn)行精確地配準(zhǔn)，一般認(rèn)為配準(zhǔn)精度應(yīng)達(dá)到亞像元級(jí)（Ito等，2003）。本文利用SAR數(shù)據(jù)自帶的衛(wèi)星軌道參數(shù)，進(jìn)行自動(dòng)精確配準(zhǔn)（Klees等，1999）。配準(zhǔn)精度可以通過(guò)相干系數(shù)圖（圖 3）間接地進(jìn)行檢查。③主干涉圖生成。配準(zhǔn)后的主、副影像相應(yīng)象元進(jìn)行復(fù)共軛相乘，得到主干涉圖（圖4）。④平地相位及地形相位的去除。在主干涉圖上去除平地相位和地形相位以得到差分干涉圖（圖5），其中，平地相位可通過(guò)對(duì)干涉條紋乘以復(fù)相位函數(shù)來(lái)去除（王超等，2002）。對(duì)于地形相位的消除，本文則是利用該地區(qū)的SRTM DEM數(shù)據(jù)和進(jìn)行干涉的兩景SAR影像的天線姿態(tài)參數(shù)生成地形的模擬干涉圖，從而得到地形相位。⑤相位解纏形成形變圖。最后對(duì)差分干涉圖采用最小費(fèi)用流（Minimum Cost Flow, MCF）方法進(jìn)行相位解纏得到形變圖（圖6）。

3.2 實(shí)驗(yàn)結(jié)果

圖4是兩景SAR影像形成的主干涉圖，干涉條紋從黑色到白色的一個(gè)顏色周期代表了相位差從-π到π的變化，可以看出在巴姆城市的東側(cè)形成2個(gè)明顯的花瓣?duì)罡缮鏃l紋，周圍的干涉條紋主要是由于地球曲率和地形造成的。圖5是去除了平地相位和地形相位的差分干涉圖，同樣，從黑色到白色的一個(gè)顏色周期代表了相位差從-π到π的變化，可以看到在巴姆城市周圍形成了4個(gè)花瓣?duì)睿ɑ蚝麪睿┑母缮鏃l紋，東側(cè)的條紋多且密集，而西側(cè)的條紋少且稀疏。北側(cè)兩組條紋的顏色變化從里到外是由白到黑，而南側(cè)條紋的顏色變化從里到外是由黑到白，說(shuō)明巴姆城市南北兩側(cè)的形變方向不同，這一點(diǎn)從最后得到的同震形變場(chǎng)也可以看到。

在干涉圖上可以看到，一些地區(qū)（巴姆市及其東南側(cè)）干涉條紋十分雜亂，這種現(xiàn)象稱為去相關(guān)現(xiàn)象，發(fā)生去相關(guān)現(xiàn)象的斑塊與城市區(qū)域和植被區(qū)域吻合得很好，所以可以初步推斷去相關(guān)現(xiàn)象的發(fā)生是由植被的變化和城市建筑的毀壞造成的。另外在巴姆市的西側(cè)有條帶狀的去相關(guān)區(qū)域，通過(guò)光學(xué)遙感影像可以看出這里是河流，所以可以初步認(rèn)定這里的去相關(guān)現(xiàn)象是由河流造成的。發(fā)生破壞的城市、植被和河流的散射特性極不穩(wěn)定，在兩景影像獲取的70天之內(nèi)發(fā)生了極大地變化，所以這些地方的干涉效果比較差，去相關(guān)現(xiàn)象十分嚴(yán)重。去相關(guān)現(xiàn)象也可以從兩景影像的相干系數(shù)圖中看出，兩景影像成像之間散射特性的改變會(huì)降低兩景影像之間的相關(guān)性，去相關(guān)現(xiàn)象越嚴(yán)重的地方，相關(guān)系數(shù)越低，如圖3中的黑色地方（相干系數(shù)接近于0）。相干系數(shù)圖也可以用來(lái)評(píng)價(jià)干涉質(zhì)量（舒寧，2003）。

舒寧（2003）提出變形監(jiān)測(cè)采用雷達(dá)干涉測(cè)量方法時(shí)必須滿足2個(gè)條件：一是2次觀測(cè)期間的變化梯度必須小于一定閾值；二是2次觀測(cè)期間，一個(gè)象元內(nèi)的所有地物散射體的特性應(yīng)具有相似性，特別是表面位置變化的誤差應(yīng)只是雷達(dá)波長(zhǎng)的10%—20%的量級(jí)。去相關(guān)效應(yīng)與干涉的兩景SAR影像之間發(fā)生變化的程度等有關(guān)，大氣影響、人類活動(dòng)等都有可能引起去相關(guān)效應(yīng)（Klees等，1999)。因此，雷達(dá)干涉測(cè)量監(jiān)測(cè)地表形變的方法適用于地表變化在一定范圍內(nèi)、干燥、少植被的地區(qū)。

圖3 相干系數(shù)圖Fig. 3 Coherence coefficient image

圖4 主干涉圖Fig. 4 The major interferogram

圖5 差分干涉圖Fig. 5 The differettntial interferogram

圖6 同震形變圖（單位：cm）Fig. 6 Coseismic deformation map (unit: cm)

3.3 地表同震形變場(chǎng)分析

從干涉圖直接解纏得到的形變是沿著視線方向（Line of Sight, LOS）的，為了將形變轉(zhuǎn)換到垂直方向上，假設(shè)LOS方向的形變量是垂直方向形變量的分量，根據(jù)公式（6）轉(zhuǎn)換：

式中，D⊥是垂直形變量；DLOS是LOS方向的形變量；α是雷達(dá)波束的局部入射角。

差分干涉圖經(jīng)過(guò)相位解纏和形變值的轉(zhuǎn)換得到垂直方向的同震形變場(chǎng)，如圖6所示。圖上的線條是形變等值線。震中的位置距離形變最大點(diǎn)距離分別是15km（最大沉降點(diǎn)）和12km（最大抬升點(diǎn)），而震中附近的形變只有1—2cm。

從圖6所示的形變圖可以看出，在遠(yuǎn)離巴姆城市的地區(qū)形變量基本在0cm左右，說(shuō)明在長(zhǎng)達(dá)70天的時(shí)間里，巴姆地區(qū)除了地震造成了嚴(yán)重形變外，幾乎沒(méi)有大的形變發(fā)生，所以得到的形變可以認(rèn)為基本是由地震造成的。地震造成的地表形變主要發(fā)生在巴姆城市的東南側(cè)和東北側(cè)。其中最大的地形抬升發(fā)生在巴姆城市東南側(cè)，如圖中白點(diǎn)所在位置，形變量達(dá)到 33.0cm，這點(diǎn)的高程約為 1042m，最大的地形沉降發(fā)生在巴姆城市東北側(cè)，如圖中黑點(diǎn)所在位置，沉降量達(dá)到 18.3cm，這點(diǎn)的高程約為 1021m。在巴姆城市西南側(cè)有少量的上升，上升量在2cm左右；巴姆城市西北側(cè)有少量的下沉，沉降量為2—4cm左右。在圖中黑色粗線所示位置可以明顯地看到一條斷層特征。斷層呈現(xiàn)南北走向，長(zhǎng)約10km。

在形變圖上過(guò)最大形變點(diǎn)分別作2條東西走向和南北走向的剖面線（剖面線A、B、C、D），如圖7所示。

圖7 形變剖面圖Fig. 7 The profiles (A, B, C and D) of the deformation map

圖7 a和b是過(guò)沉降最大點(diǎn)的剖面線，圖7c和d是過(guò)抬升最大點(diǎn)的剖面線，其中，圖7a、c為東走向剖面圖，圖7b、d為南北走向剖面圖。從這4個(gè)剖面圖可以看出，巴姆地區(qū)的形變基本表現(xiàn)為巴姆城市北部地區(qū)的地塊有所沉降，而巴姆城市南部地區(qū)的地塊被抬升。

4 結(jié)論與討論

本文利用D-InSAR技術(shù)對(duì)伊朗巴姆地區(qū)2003年12月26日6.7級(jí)地震造成的地表形變進(jìn)行了研究，得到了該地區(qū)垂直向的三維同震形變場(chǎng)（圖6）。

（1）巴姆地震造成的地表形變主要分布在巴姆市的東側(cè)，垂直形變量在-18.3cm到33.0cm之間。這個(gè)結(jié)果與李毅等（2008）、王志勇等（2008）得到的結(jié)論接近，但是存在幾毫米到幾厘米的誤差。主要是因?yàn)樗麄兊玫降氖抢走_(dá)視線方向LOS上的形變，而本文將這個(gè)形變量轉(zhuǎn)換成垂直方向的形變。從三維同震形變場(chǎng)的圖上可以看到巴姆地區(qū)的形變主要表現(xiàn)為：北部地塊發(fā)生沉降，最大沉降地點(diǎn)位于巴姆城市東北側(cè)，距震中約 15km；南部地塊被抬升，最大抬升地點(diǎn)位于巴姆城市東南側(cè)，距震中約12km。震中附近形變量只有1—2cm。另外，地震還在巴姆城市南側(cè)造成了長(zhǎng)約10km、南北走向的斷層。

（2）在巴姆地區(qū)的某些區(qū)域發(fā)生了去相關(guān)效應(yīng)，這些現(xiàn)象主要發(fā)生在巴姆城市、河流和植被茂盛的區(qū)域。這主要是河流和植被的散射特性不穩(wěn)定，在相隔2個(gè)多月的時(shí)間里，植被和河流散射特性已經(jīng)發(fā)生了較大的變化，從而產(chǎn)生了去相關(guān)現(xiàn)象。巴姆城市地區(qū)則由于建筑物被地震嚴(yán)重地?fù)p毀，散射特性發(fā)生了極大變化，所以產(chǎn)生了去相關(guān)效應(yīng)。一般認(rèn)為，沙漠地區(qū)比密林地區(qū)的干涉條件好，干燥環(huán)境比潮濕環(huán)境的干涉效果好，波長(zhǎng)較長(zhǎng)的雷達(dá)圖像會(huì)得到較好效果的干涉圖（舒寧，2003）。

（3）巴姆地區(qū)差分干涉圖是利用精度較低的SRTM DEM數(shù)據(jù)去除地形相位信息的，這樣造成的地形誤差較大，在今后的研究中，在提高地形測(cè)量精度方面還需要進(jìn)一步努力。另外，由于D-InSAR技術(shù)對(duì)地表散射體散射特性的穩(wěn)定性要求較高，所以對(duì)于植被豐富，氣候潮濕的地區(qū)，就很難產(chǎn)生較好的干涉效果，因此在去相關(guān)問(wèn)題上還需進(jìn)一步的研究。

利用雷達(dá)影像差分干涉測(cè)量得到的形變圖僅僅包含了LOS向的形變信息（張紅，2002），即使得到了垂直的形變量，也是基于LOS方向形變是垂直形變量的分量的假設(shè)。要得到地表形變的真實(shí)方向及大小，需要兩個(gè)影像對(duì)進(jìn)行D-InSAR形變分析，從而合成出真實(shí)的形變方向及大小。另外，為了對(duì)地表形變場(chǎng)認(rèn)識(shí)得更加全面充分，還需要結(jié)合GPS、GIS等手段建立形變模型，從而確定研究區(qū)的形變機(jī)制。如果能夠通過(guò)技術(shù)進(jìn)步提高雷達(dá)干涉測(cè)量的精度，并降低觀測(cè)成本，建立長(zhǎng)期的觀測(cè)機(jī)制，周期性的獲取高質(zhì)量的SAR影像用來(lái)觀測(cè)地表形變量，則可以對(duì)地震、滑坡、火山爆發(fā)、地面沉降等災(zāi)害的預(yù)報(bào)提供有價(jià)值的觀測(cè)數(shù)據(jù)，從而在災(zāi)難發(fā)生前后采取有效措施以降低損失，真正讓這項(xiàng)技術(shù)發(fā)揮價(jià)值。

查顯杰，2007. InSAR形變測(cè)量及其在地震學(xué)中應(yīng)用的研究. 合肥：中國(guó)科學(xué)技術(shù)大學(xué).

何秀鳳，何敏，2012. InSAR對(duì)地觀測(cè)數(shù)據(jù)處理方法. 北京：科學(xué)出版社.

季靈運(yùn)，許建東，2009. 利用D-InSAR和AZO技術(shù)獲取Bam地震同震三維形變場(chǎng). 大地測(cè)量地球動(dòng)力學(xué)，29（6）40—44.

李毅，柳林濤，薛懷平等，2008. 用D-InSAR研究巴姆地震形變場(chǎng). 大地測(cè)量與地球動(dòng)力學(xué)，28（1）：32—35.

劉大杰，胡叢瑋，1999. 應(yīng)用GPS監(jiān)測(cè)城市地表形變的初步分析. 地殼形變與地震，19（1）：37—42.

劉國(guó)祥，2006. 利用雷達(dá)干涉技術(shù)監(jiān)測(cè)區(qū)域地表形變. 北京：測(cè)繪出版社.

劉國(guó)祥，丁曉利，李志偉等，2002. ERS衛(wèi)星雷達(dá)干涉測(cè)量：1999年臺(tái)灣集集大地震震前和同震地表位移. 地球物理學(xué)報(bào)，45（增刊）：165—174.

凌勇，曾祺明，羅揚(yáng)，趙永紅，2006. 巴姆地震變形場(chǎng)和應(yīng)力場(chǎng)：Ⅰ. 用差分干涉雷達(dá)和Okada方法求解. 巖石學(xué)報(bào)，22（9）：2367—2374.

羅揚(yáng)，施旭，趙永紅，2006. 巴姆地震變形場(chǎng)和應(yīng)力場(chǎng)：Ⅱ. 用FEPG有限元方法求解. 巖石學(xué)報(bào)，22（9）：2375—2380.

馬超，單新建，2004. 星載合成孔徑雷達(dá)差分干涉測(cè)量（D-InSAR）技術(shù)在形變監(jiān)測(cè)中的應(yīng)用概述. 中國(guó)地震，20（4）：410—418.

單新建，柳稼航，馬超，2004. 2001年昆侖山口西8.1級(jí)地震同震形變場(chǎng)特征的初步分析. 地震學(xué)報(bào)，26（5）：474—480.

邵蕓，譚衢霖，劉浩等，2002. 雷達(dá)差分干涉測(cè)量技術(shù)及其應(yīng)用研究——以1997年西藏瑪尼7.9級(jí)地震區(qū)域形變測(cè)量為例. 地球物理學(xué)報(bào)，45（增刊）：205—213.

舒寧，2003. 雷達(dá)影像干涉測(cè)量原理. 武漢：武漢大學(xué)出版社.

孫建寶，梁芳，沈正康等，2008. 汶川 MS8.0地震 InSAR形變觀測(cè)及初步分析. 地震地質(zhì)，30（3）：789—795.

王超，張紅，劉智，2002. 星載合成孔徑雷達(dá)干涉測(cè)量. 北京：科學(xué)出版社.

王志勇，張繼賢，張永紅等，2008. 伊朗巴姆 6.5級(jí)地震同震形變場(chǎng)的獲取與解譯. 地震研究，31（1）：70—76.

張保淑，2009. 科技之謎：追尋巴姆古城堡. http://www.people.com.cn/GB/keji/1058 /2284182.html.

張紅，王超，劉智，2002. 獲取張北地震同震形變場(chǎng)的差分干涉測(cè)量技術(shù). 中國(guó)圖象圖形學(xué)報(bào)，5A（6）：497—500.

Bam (Iran), December 2003 [2009-04-26]. http://earth.esa.int/ew/earthquakes/Bam_Iran_Earth- quake.

Ito Y., Hosokawa M., 2003. A degree-of-damage estimation model of earthquake damage using interferometric SAR data. Electrical Engineering in Japan, 143 (3): 49—57.

Klees R., Massonnet D., 1999. Deformation measurements using SAR interferometry: potenital and limitation.Geologie en Mijnbouw, 77 (2): 161—176.

Zebker H.A., Rousen P., 1994. On the derivation of coseismic displacement fields using differential radar interferometry: The Landers earthquake. Geoscience and Remote Sensing Symposium, 1994. IGARSS '94.Surface and Atmospheric Remote Sensing: Technologies, Data Analysis and Interpretation. 1: 286—288.