楊　珍，張永志，吳　然，張文軍，葉　凱

(1.長安大學(xué) 地質(zhì)工程與測繪學(xué)院，陜西 西安 710000；2.甘肅省測繪地理信息局 地圖院,甘肅 蘭州 730000)

?

基于D-InSAR技術(shù)的當(dāng)雄地震形變場提取研究

楊珍1，張永志1，吳然1，張文軍2，葉凱1

(1.長安大學(xué) 地質(zhì)工程與測繪學(xué)院，陜西 西安 710000；2.甘肅省測繪地理信息局 地圖院,甘肅 蘭州 730000)

摘要：隨著D-InSAR技術(shù)在同震形變量測方面的優(yōu)勢日益突出，文中以2008年10月6日當(dāng)雄Mw6.3級地震為對象，采用ENVISAT搭載的ASAR傳感器c波段在圖像模式下的5景l(fā)evel 0級影像，以90 m分辨率的SRTM DEM作為外部DEM，使用GAMMA軟件采用二軌差分干涉測量的方法獲取當(dāng)雄地震同震形變場，并對其進行分析，確定震中位置為90.372°E，29.734°N，得到同震垂直形變約為33 cm，形變量精度達到厘米級。

關(guān)鍵詞：D-InSAR；地震；ENVISAT；二軌差分；地震形變場

2008-10-06T16 p0 min，在西藏自治區(qū)當(dāng)雄縣發(fā)生了Mw6.3級地震，監(jiān)測震間微小形變場，研究斷裂活動引起地表位移的時空演化特征，捕獲可能的震前形變異常，為地震的監(jiān)測預(yù)測提供可靠依據(jù)，推進地震預(yù)報工作是地震研究迫切需要解決的問題。

合成孔徑雷達差分干涉測量(D-InSAR)技術(shù)是雷達干涉測量應(yīng)用的一個拓展，雷達干涉圖的差分可以監(jiān)測雷達視線方向厘米級或更微小的地球表面形變，以揭示許多地球物理現(xiàn)象，如地震形變、火山活動、冰川漂移、地面沉降以及山體滑坡等[1]。

近年來，國內(nèi)外很多學(xué)者利用D-InSAR技術(shù)研究地震形變，許才軍等利用InSAR數(shù)據(jù)提取及分析汶川地震形變場[2]，單新建等研究了汶川地震前的垂直形變場變化特征[3]。但是以往的研究大多針對的是大地震，本文選取強震級別的當(dāng)雄地震為研究對象，基于ENVISAT ASAR數(shù)據(jù)，采用D-InSAR方法，提取及分析當(dāng)雄Mw6.3級地震形變場。

1D-InSAR基本原理

InSAR測量的相位由5部分組成[4]：

(1)

等式右端的5個分量分別為地形相位、所求形變相位、平地相位(參考相位)、大氣相位、噪聲引起相位。地表變化前后的2幅影像的干涉圖中同時包含地形信息和形變信息，由式(1)知差分干涉的目的是使用一定的方法去掉一些相位信息，最終獲得形變相位的過程[5]。

以二軌法為例，將2幅SAR影像進行干涉生成干涉圖，生成的干涉圖包含地形變化信息，隨后再利用DEM數(shù)據(jù)進行相位模擬，模擬得到地形信息的條紋圖，再將SAR影像干涉圖與DEM模擬的條紋圖疊加去除地表形變，得到地表形變信息。

2基于GAMMA軟件的二軌法數(shù)據(jù)處理流程

本文采用GAMMA軟件來實現(xiàn)二軌法數(shù)據(jù)處理，具體的實現(xiàn)流程如圖1所示，首先在聚焦模塊下，利用儀器的參數(shù)文件及精密軌道文件對數(shù)據(jù)進行聚焦得到SLC文件，然后與常規(guī)差分處理流程一樣，包括影像的配準(zhǔn)、濾波、干涉圖生成、基線估算、外部相位模擬、差分處理、相位解纏、形變量計算、地理編碼等處理步驟。

圖1　GAMMA軟件二軌法數(shù)據(jù)處理流程

3數(shù)據(jù)處理

3.1數(shù)據(jù)源

研究區(qū)域當(dāng)雄地形起伏在4 000～6 000 m，地形高低錯落，溝谷縱橫。本文獲得了5景ENVISAT ASAR level 0 級數(shù)據(jù)，天線文件和ASAR儀器參數(shù)文件(ESA提供)及精密軌道數(shù)據(jù)文件(DELFT提供)，數(shù)據(jù)情況見表1。

表1　ENVISAT ASAR數(shù)據(jù)列表

DEM數(shù)據(jù)為SRTM DEM，分辨率30 m，平面基準(zhǔn)WGS-84，高程基準(zhǔn)EGM96，標(biāo)稱精度±10 m。

3.2數(shù)據(jù)選擇

選擇適合干涉的SAR影像是進行任何干涉處理的首要步驟，也是非常關(guān)鍵的步驟。數(shù)據(jù)源的好壞直接影響到后續(xù)步驟，為了得到最好的結(jié)果，數(shù)據(jù)選擇需要考慮的參數(shù)包括傳感器類型、視角、幾何基線/空間基線、時間基線、影像獲取時刻、相干性、氣象條件。

將5景數(shù)據(jù)處理成SLC數(shù)據(jù)后，對數(shù)據(jù)進行濾波處理(前置濾波)，配準(zhǔn)以與地震發(fā)生時間最近的20080921.slc影像為配準(zhǔn)主影像，將其余四景影像配準(zhǔn)至20080921，配準(zhǔn)方法采用精密軌道數(shù)據(jù)和強度互相關(guān)及偏移量跟蹤法，配準(zhǔn)精度可達到1/10像元。再根據(jù)表2，干涉成像。基線估計、差分處理，得到10幅干涉圖，時空基線的結(jié)果見表2。

表2　當(dāng)雄Envisat ASAR干涉組合

根據(jù)衛(wèi)星影像的覆蓋范圍，InSAR數(shù)據(jù)處理對空間基線、時間基線的要求，差分干涉的效果，故選取20080921～20081026和20080921～20090104，如圖2所示，進行后續(xù)操作，最后選取相干性好的一幅來提取同震形變場。

圖2　差分處理后的干涉圖

3.3濾波處理

在影像進行精確配準(zhǔn)前，需要進行頻率域的濾波處理，即前置濾波，前置濾波是在方位向和距離向分別進行的。一般而言，如果頻譜偏移過大，方位向濾波應(yīng)該在配準(zhǔn)前進行，可提高配準(zhǔn)精度，而距離向濾波則必須在配準(zhǔn)后的主從影像之間才能進行。方位向和距離向前置濾波的具體流程相類似，首先都是對主從影像在方位向和距離向進行頻率域分析，然后通過帶通濾波器將兩者的頻譜非重疊部分濾除。

圖2得到差分干涉圖去除地形相位及平地效應(yīng)，但其他因素導(dǎo)致的隨機噪聲也會影響干涉圖的質(zhì)量，對影像也要進行濾波(后置濾波)。Adaptive濾波運算速度快、效率高、效果好、占用內(nèi)存小，對噪聲較小的區(qū)域能夠得到正確的解纏值，且能夠較多保存邊緣信息，對于高分辨率數(shù)據(jù)處理效果非常好，如TerraSAR-X或COSMO-SkyMed；Boxcar濾波是使用局部干涉條紋的頻率來優(yōu)化濾波器，該方法盡可能的保留微小的干涉條紋；Goldstein自適用濾波方法的濾波器是可變的，提高干涉條紋的清晰度，減少由空間基線或時間基線引起的失相干的噪聲，是目前最流行的方法。

本文采取了Goldstein自適用濾波的方法，且進行兩次，第一次濾波窗口64×64，第二次濾波窗口32×32。

3.4相位解纏與去除基線誤差

相位解纏方法很多，大致分為兩類：①基于路徑的相位解纏算法；②基于最小范數(shù)的相位解纏算法。最常用效果較好的方法有區(qū)域增長法、最小費用流法、枝切法、最小二乘法。本文通過比較，最后采用最小費用流法。在進行相位解纏時避開低相干性區(qū)域，生成有效掩膜文件，提取相干系數(shù)較高的點，建立Delaunay網(wǎng)，使纏繞相位與解纏相位的差異最小化，然后將最小化問題轉(zhuǎn)換為最小費用流問題。

據(jù)有關(guān)實驗研究證明在海拔5 000 m的高原水汽變化引起的誤差小于1 mm，當(dāng)雄海拔高，屬于干燥性氣候，水汽總含量相對較低，所以大氣誤差造成的影響本文不考慮。

采用Delft精密軌道數(shù)據(jù)，軌道精度為5～10 cm，基線誤差屬于系統(tǒng)誤差，本文采用遠場最佳二次多項式擬合的方法去除非線性軌道殘余相位信息，最后得到去除基線軌道誤差的干涉圖。

3.5形變量計算

形變量計算是將差分解纏相位轉(zhuǎn)換為視線(LOS)向形變量，LOS向形變同時包含了東西向、南北向、及垂直向3個方向的分量，經(jīng)過計算分析發(fā)現(xiàn)本次地震LOS向形變在東西向和南北向的形變量非常小，經(jīng)過處理轉(zhuǎn)換最終得到地表垂直方向形變。

3.6地理編碼

計算出每個像素對應(yīng)的地表幾何信息后，得到一個在影像坐標(biāo)系下的點陣圖，還需要把各種數(shù)據(jù)從影像坐標(biāo)系轉(zhuǎn)換到統(tǒng)一的地理坐標(biāo)系下。將雷達坐標(biāo)系下的原始數(shù)據(jù)通過一定的幾何校正方法消除由軌道、傳感器、地球模型引起的扭曲和畸變，然后變換到某種制圖參考系中，地理編碼之后影像坐標(biāo)系與提供的外部DEM坐標(biāo)系一致，故形變監(jiān)測結(jié)果為WGS-84坐標(biāo)系統(tǒng)。圖3為地理編碼后的形變圖，總體上兩幅圖的同震形變場保存高度一致，圖3(a)相干性較好，圖3(b)失相干現(xiàn)象比較嚴(yán)重，故本文選取地理編碼后的20080921～20081026形變圖進行結(jié)果分析。

圖3　地理編碼后的形變圖

4當(dāng)雄地震形變結(jié)果分析

根據(jù)圖4地震同震形變圖，從宏觀上分析，沉降區(qū)的中心位置就是震中位置，形變條紋左上側(cè)圓環(huán)沉降，區(qū)域面積大于右側(cè)的隆升區(qū)，越靠近地表破裂帶，干涉條紋越密集，地表形變量越大，越遠離地表破裂帶，干涉條紋越疏松，地表形變量越小，地表破裂帶沉降區(qū)的最大變形量約為33 cm，隆升區(qū)最大變形量約為5 cm，圓環(huán)的中心位置即為震中位置。可判斷出震中位置在90.372°E，29.734°N。與喬學(xué)軍推算的震中位置90.374°E，29.745°N非常接近[6]，相差約640 m。

4.1剖面分析

取兩條測線進行剖面分析，兩條測線的位置均過震，一條對應(yīng)雷達的視線方向且近似垂直地表破裂帶[7]，如圖4中的WE，一條近似平行于地表破裂帶，如圖4中的NS，剖面圖結(jié)果如圖5所示。

圖4　當(dāng)雄地震同震形變

從圖5(a)中可以看出，從起點至大約12.5 km處，形變值約從-0.03變化到-0.32，在大約12.5～19 km處，形變值又由約-0.32變化到0，在約19～31 km處，形變值在0.03左右波動，這說明地表變形從沉降逐漸變?yōu)槁∩?，最大沉降量約為32 cm。從圖5(b)中可以看出，剖面圖近似呈對稱狀，在距離起點10.5 km處，沉降量達最大值約33 cm，說明離震中越近，地表位錯變換越劇烈。

4.2等直線分析

對同震形變做等值線圖[8]，如圖6所示，可以看出最大形變量在32 cm左右，根據(jù)圖6的形變場特征可以看出，與圖4、圖5中的形變場保持高度一致，等值線中因部分像元無變形值出現(xiàn)“孤島”現(xiàn)象，導(dǎo)致這種情況的原因有多種，大氣改正失敗，相位解纏失敗，山體滑坡等，另一方面來說，對等值線進行分析，可評判InSAR結(jié)果的質(zhì)量。

圖5　同震形變場剖面圖

圖6　同震形變場等值線

5結(jié)束語

本文通過D-InSAR的兩軌差分方法提取當(dāng)雄地震的同震形變場，得到的同震垂直形變約為33 cm，確定的震中位置：90.372°E，29.734°N。采用Goldstein自適應(yīng)濾波2次，減少斑點噪聲，相位數(shù)據(jù)連續(xù)，周期性更加明顯。但是本文沒有考慮大氣誤差對形變結(jié)果的影響，還有待進一步研究。D-InSAR技術(shù)對于跨斷層的震間形變監(jiān)測具有較大潛力，特別是對GPS測站分布稀疏的地區(qū)，利用D-InSAR獲取的震間形變場對于理解活動斷層的地殼形變特征、地球動力學(xué)問題具有重要的應(yīng)用價值。

參考文獻：

[1]廖明生，林琿.雷達干涉測量：原理與信號處理基礎(chǔ)[M].測繪出版社，2003.

[2]許才軍，江國焰，王浩，等.基于GIS的InSAR結(jié)果分析方法及在汶川Mw7.9級地震同震解釋中的應(yīng)用[J].武漢大學(xué)學(xué)報(信息科學(xué)版)，2011，36(4):379-383.

[3]單新建，宋小剛，韓宇飛，等.汶川Ms8.0地震前InSAR垂直形變場變化特征研究[J].地球物理學(xué)報，2009(11):2739-2745.

[4]吳中海，葉培盛，吳珍漢.2008年10月6日西藏當(dāng)雄Ms6.6級強震的地震烈度控震構(gòu)造和發(fā)震機理[J].地質(zhì)通報，2009，28(6):713-725.

[5]張軍，尼瑪，尚榮波.2008年西藏當(dāng)雄6.6級地震淺[J].高原地震，2009，20(4):40-44.

[6]喬學(xué)軍，王琪，楊少敏，等.2008年新疆烏恰Mw6.7地震震源機制與形變特征的InSAR研究[J].地球物理學(xué)報，2014(6):1805-1813.

[7]雷廣淵，周輝.基于SBAS技術(shù)的金屬礦山沉陷規(guī)律研究[J].測繪工程，2015，24(3):40-46.

[8]谷延超，范菇宇，范東明，等.基于不同濾波算法差異的機載LiDAR數(shù)據(jù)橋梁提取[J].測繪工程，2014，23(11):67-70.

[責(zé)任編輯：李銘娜]

DOI:10.19349/j.cnki.issn1006-7949.2016.08.002

收稿日期：2015-09-06

基金項目：國家自然科學(xué)基金資助項目(41374028)

作者簡介：楊珍(1984-)，女，博士研究生.

中圖分類號：TP751

文獻標(biāo)識碼：A

文章編號：1006-7949(2016)08-0005-06

On the extraction of surface deformation caused by Damxung earthquake based on D-InSAR technology

Yang Zhen1,Zhang Yongzhi1,Wu Ran1,Zhang Wenjun2,Ye Kai1

(1.School of Geology Engineering and Geomatics,Chang’an University,Xi’an 710000,China;2.Map Institute,Gansu Administration of Surveying and Mapping Geoinformation,Lanzhou 730000,China)

Abstract：As advantages of D-InSAR technology in coseismic deformation and measurement are becoming more and more prominent,this paper takes Damxung Mw 6.3 earthquake on October 6,2008 as the case,adopting the image of 5 scenes and level 0 of C band in the image mode of ASAR sensor being equipped on ENVISAT,using SRTM DEM with 90 m resolution ratio as the external DEM,and with GAMMA software and the approach of two-pass differential InSAR technology to obtain the coseismic displacement of Damxung earthquake.And analysis is made to determine the location of the epicenter as 90.372°E,29.734°N,gaining about 35 cm of the coseismic vertical deformation,of which the deformation precision reaches the cm-level.

Key words：D-InSAR;earthquake;ENVISAT;two-pass differential;coseismic displacement