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

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

(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)雄地震同震形變場，并對其進(jìn)行分析，確定震中位置為90.372°E，29.734°N，得到同震垂直形變約為33 cm，形變量精度達(dá)到厘米級。

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

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

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

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

1D-InSAR基本原理

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

(1)

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

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

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

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

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

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

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

研究區(qū)域當(dāng)雄地形起伏在4 000～6 000 m，地形高低錯(cuò)落，溝谷縱橫。本文獲得了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影像是進(jìn)行任何干涉處理的首要步驟，也是非常關(guān)鍵的步驟。數(shù)據(jù)源的好壞直接影響到后續(xù)步驟，為了得到最好的結(jié)果，數(shù)據(jù)選擇需要考慮的參數(shù)包括傳感器類型、視角、幾何基線/空間基線、時(shí)間基線、影像獲取時(shí)刻、相干性、氣象條件。

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

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

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

圖2　差分處理后的干涉圖

3.3濾波處理

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

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

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

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

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

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

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

3.5形變量計(jì)算

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

3.6地理編碼

計(jì)算出每個(gè)像素對應(yīng)的地表幾何信息后，得到一個(gè)在影像坐標(biāo)系下的點(diǎn)陣圖，還需要把各種數(shù)據(jù)從影像坐標(biāo)系轉(zhuǎn)換到統(tǒng)一的地理坐標(biāo)系下。將雷達(dá)坐標(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形變圖進(jìn)行結(jié)果分析。

圖3　地理編碼后的形變圖

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

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

4.1剖面分析

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

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

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

4.2等直線分析

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

圖5　同震形變場剖面圖

圖6　同震形變場等值線

5結(jié)束語

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

參考文獻(xiàn)：

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

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

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

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

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

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

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

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

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

DOI:10.19349/j.cnki.issn1006-7949.2016.08.002

收稿日期：2015-09-06

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

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

中圖分類號：TP751

文獻(xiàn)標(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