基于D-InSAR 技術的伊朗巴姆地震地表形變監測1

寇 程 柯長青

（南京大學地理信息科學系，南京 210046）

引言

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

國內外學者對D-InSAR技術監測地表形變監測進行了廣泛研究。D-InSAR技術被用于臺灣集集地區地震（劉國祥等，2002）、昆侖山口西地震（單新建等，2004）、張北地震（張紅等，2002）、汶川地震（孫建寶等，2008）、Landers地震（Zebker等，1994）和Hyogoken-Nanbu地震（Ito等，2003）同震地表位移監測和同震形變場參數提取。這些研究證明了 D-InSAR技術在地表形變監測方面的可行性，其監測精度可達到厘米級。

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

1 研究區地理背景及數據源

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

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

表1 ENVISAT ASAR數據部分參數Table 1 Some parameters of ENVISAT ASAR data

2 D-InSAR技術原理

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

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

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

逐象元進行復共軛相乘后得到：

因此，干涉相位為：

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

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

3 實驗驗證

3.1 實驗步驟

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

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

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

3.2 實驗結果

圖4是兩景SAR影像形成的主干涉圖，干涉條紋從黑色到白色的一個顏色周期代表了相位差從-π到π的變化，可以看出在巴姆城市的東側形成2個明顯的花瓣狀干涉條紋，周圍的干涉條紋主要是由于地球曲率和地形造成的。圖5是去除了平地相位和地形相位的差分干涉圖，同樣，從黑色到白色的一個顏色周期代表了相位差從-π到π的變化，可以看到在巴姆城市周圍形成了4個花瓣狀（或蝴蝶狀）的干涉條紋，東側的條紋多且密集，而西側的條紋少且稀疏。北側兩組條紋的顏色變化從里到外是由白到黑，而南側條紋的顏色變化從里到外是由黑到白，說明巴姆城市南北兩側的形變方向不同，這一點從最后得到的同震形變場也可以看到。

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

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

圖3 相干系數圖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 地表同震形變場分析

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

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

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

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

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

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

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

4 結論與討論

本文利用D-InSAR技術對伊朗巴姆地區2003年12月26日6.7級地震造成的地表形變進行了研究，得到了該地區垂直向的三維同震形變場（圖6）。

（1）巴姆地震造成的地表形變主要分布在巴姆市的東側，垂直形變量在-18.3cm到33.0cm之間。這個結果與李毅等（2008）、王志勇等（2008）得到的結論接近，但是存在幾毫米到幾厘米的誤差。主要是因為他們得到的是雷達視線方向LOS上的形變，而本文將這個形變量轉換成垂直方向的形變。從三維同震形變場的圖上可以看到巴姆地區的形變主要表現為：北部地塊發生沉降，最大沉降地點位于巴姆城市東北側，距震中約 15km；南部地塊被抬升，最大抬升地點位于巴姆城市東南側，距震中約12km。震中附近形變量只有1—2cm。另外，地震還在巴姆城市南側造成了長約10km、南北走向的斷層。

（2）在巴姆地區的某些區域發生了去相關效應，這些現象主要發生在巴姆城市、河流和植被茂盛的區域。這主要是河流和植被的散射特性不穩定，在相隔2個多月的時間里，植被和河流散射特性已經發生了較大的變化，從而產生了去相關現象。巴姆城市地區則由于建筑物被地震嚴重地損毀，散射特性發生了極大變化，所以產生了去相關效應。一般認為，沙漠地區比密林地區的干涉條件好，干燥環境比潮濕環境的干涉效果好，波長較長的雷達圖像會得到較好效果的干涉圖（舒寧，2003）。

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

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

查顯杰，2007. InSAR形變測量及其在地震學中應用的研究. 合肥：中國科學技術大學.

何秀鳳，何敏，2012. InSAR對地觀測數據處理方法. 北京：科學出版社.

季靈運，許建東，2009. 利用D-InSAR和AZO技術獲取Bam地震同震三維形變場. 大地測量地球動力學，29（6）40—44.

李毅，柳林濤，薛懷平等，2008. 用D-InSAR研究巴姆地震形變場. 大地測量與地球動力學，28（1）：32—35.

劉大杰，胡叢瑋，1999. 應用GPS監測城市地表形變的初步分析. 地殼形變與地震，19（1）：37—42.

劉國祥，2006. 利用雷達干涉技術監測區域地表形變. 北京：測繪出版社.

劉國祥，丁曉利，李志偉等，2002. ERS衛星雷達干涉測量：1999年臺灣集集大地震震前和同震地表位移. 地球物理學報，45（增刊）：165—174.

凌勇，曾祺明，羅揚，趙永紅，2006. 巴姆地震變形場和應力場：Ⅰ. 用差分干涉雷達和Okada方法求解. 巖石學報，22（9）：2367—2374.

羅揚，施旭，趙永紅，2006. 巴姆地震變形場和應力場：Ⅱ. 用FEPG有限元方法求解. 巖石學報，22（9）：2375—2380.

馬超，單新建，2004. 星載合成孔徑雷達差分干涉測量（D-InSAR）技術在形變監測中的應用概述. 中國地震，20（4）：410—418.

單新建，柳稼航，馬超，2004. 2001年昆侖山口西8.1級地震同震形變場特征的初步分析. 地震學報，26（5）：474—480.

邵蕓，譚衢霖，劉浩等，2002. 雷達差分干涉測量技術及其應用研究——以1997年西藏瑪尼7.9級地震區域形變測量為例. 地球物理學報，45（增刊）：205—213.

舒寧，2003. 雷達影像干涉測量原理. 武漢：武漢大學出版社.

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

王超，張紅，劉智，2002. 星載合成孔徑雷達干涉測量. 北京：科學出版社.

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

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

張紅，王超，劉智，2002. 獲取張北地震同震形變場的差分干涉測量技術. 中國圖象圖形學報，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.