基于子帶干涉測量技術(shù)的巴基斯坦地震形變獲取研究

庾露， 單新建， 宋小剛， 屈春燕

1 中國地震局地質(zhì)研究所， 地震動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室， 北京　100029 2 廣西壯族自治區(qū)遙感中心， 南寧　530023

基于子帶干涉測量技術(shù)的巴基斯坦地震形變獲取研究

庾露1,2， 單新建1*， 宋小剛1， 屈春燕1

1 中國地震局地質(zhì)研究所， 地震動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室， 北京100029 2 廣西壯族自治區(qū)遙感中心， 南寧530023

摘要2013年9月24日發(fā)生在巴基斯坦俾路支省(Balochistan)境內(nèi)阿瓦蘭縣(Awaran)的MW7.7級地震，在地表產(chǎn)生了最大達(dá)10 m的滑動(dòng)量.利用TerraSAR-X短波雷達(dá)數(shù)據(jù)獲取的InSAR同震形變場產(chǎn)生了密集且大范圍的干涉條紋，給后續(xù)的相位解纏帶來困難.而子帶干涉法是一種無需或只需進(jìn)行少量相位解纏，即可獲得絕對相位差的新方法.其主要思路是通過縮減帶寬以增長波長，從而減少干涉條紋數(shù)，降低解纏難度或不需解纏直接得到絕對相位差.但由于帶寬的縮減，導(dǎo)致噪聲的增大和旁瓣帶來的額外干擾，使干涉圖質(zhì)量下降，因此在子帶干涉參數(shù)選取、噪聲濾波以及處理流程等方面需要特殊處理，特別是子帶的中心頻率和帶寬的選取會(huì)很大程度地影響測量精度.首先選取典型DEM實(shí)驗(yàn)區(qū)，以干涉圖相干性和誤差為評價(jià)指標(biāo)，利用逐步參數(shù)選取法，研究相關(guān)參數(shù)對子帶干涉測量的影響，制定最優(yōu)的參數(shù)方案，認(rèn)識參數(shù)選取的原則和方法.在此基礎(chǔ)上，將子帶干涉應(yīng)用于巴基斯坦地震的同震形變場獲取.最后，將子帶干涉、Landsat 8光學(xué)影像的交叉頻譜相關(guān)法、offset-tracking、常規(guī)DInSAR獲取的同震形變場進(jìn)行比較，并與模型擬合的形變場進(jìn)行對比分析.結(jié)果表明，子帶干涉雖然會(huì)受失相干的影響，其提取的形變場范圍相較于Landsat 8和offset-tracking有所缺失，但在共同覆蓋的區(qū)域其精度和噪聲水平更優(yōu)，相比較于常規(guī)DInSAR，更適用于條紋密集和形變量大的地區(qū).

關(guān)鍵詞子帶干涉； 巴基斯坦地震； 同震形變

1引言

常規(guī)合成孔徑雷達(dá)干涉測量技術(shù)(Interferometry Synthetic Aperture Radar，InSAR)在地形測繪(杜亞男等，2015；占文俊等，2015)和形變監(jiān)測方面(劉云華等，2014;李永生等，2015；單新建等，2015)的發(fā)展已經(jīng)非常成熟，其成功應(yīng)用的關(guān)鍵取決于數(shù)據(jù)處理中的一些關(guān)鍵步驟，相位解纏便是其中之一，解纏結(jié)果的質(zhì)量直接影響測量的精度(廖明生和林琿，2003).針對解纏問題，眾多學(xué)者提出了不同的解纏算法，例如：Costantini(1998)提出的最小費(fèi)用流(minimum cost flow, MCF)法；Goldstein等(1998)提出的枝切法.這些解纏算法的提出，極大地推動(dòng)了InSAR測量技術(shù)的發(fā)展.但在一些地表環(huán)境復(fù)雜地區(qū)，大部分的解纏方法存在一定的局限性.例如：一方面由于解纏算法的不同，會(huì)造成InSAR的結(jié)果各異，從而帶來了不確定性；另一方面，對于地形復(fù)雜或相位梯度過大的區(qū)域，干涉條紋會(huì)過于密集，再加上去相干因素的影響，濾波會(huì)出現(xiàn)條紋混疊現(xiàn)象，導(dǎo)致后續(xù)錯(cuò)誤的解纏結(jié)果.如何在應(yīng)用中解決這些問題，已成為當(dāng)今研究的熱點(diǎn)問題之一，在這一需求下，子帶干涉成為全球關(guān)注的焦點(diǎn).子帶干涉可以增加有效波長，大大減少條紋數(shù)量，從而降低解纏難度或直接得到絕對相位.該方法最早由Madsen等(1993)針對當(dāng)時(shí)NASA DC-8空基SAR數(shù)據(jù)提取地形提出，其基本原理是利用帶通濾波器，在距離向的頻率域上將主副兩景單視復(fù)數(shù)影像(single look complex，SLC)分別分割成上下兩個(gè)子帶，再通過上-上子帶、下-下子帶干涉處理后，再次進(jìn)行上下子帶二次干涉，從而模擬獲得較長的波長.在其后的一段時(shí)間內(nèi)，受制于天基雷達(dá)波帶寬的限制，這一技術(shù)一直未能得到廣泛的應(yīng)用.隨著新一代高分辨率，寬帶寬的雷達(dá)衛(wèi)星如TerraSAR、CosmoSkyMed的相繼發(fā)射并投入使用，雷達(dá)波帶寬已達(dá)到300 MHz，已可以提取出信噪比很高的子帶.在這一背景下，又有一批學(xué)者重新對該方法進(jìn)行了研究.如Derauw等(2010)分析了子帶干涉法下地面散射體的相干性；Brcic等(2008)通過子帶干涉法利用TerraSAR-X數(shù)據(jù)，提取了南美Salar de Uyuni鹽湖區(qū)域的DEM，并與常規(guī)InSAR的結(jié)果進(jìn)行了對比.目前對于子帶干涉技術(shù)的應(yīng)用，主要集中在提取DEM數(shù)據(jù)上，而對于地表形變的提取則少有涉及.

而對于本文所研究的2013年發(fā)生在巴基斯坦境內(nèi)的MW7.7級地震的同震地表形變場提取.目前使用遙感手段，主要是利用Landsat 8光學(xué)衛(wèi)星震前和震后影像，基于交叉頻譜相關(guān)(cross-spectrum correlation)算法(Leprince et al.，2007a)獲得(馮光財(cái)?shù)龋?015；Avouac et al.，2014).雷達(dá)影像方面，主要使用了TerraSAR-X StripMap模式和ScanSAR模式像對提取同震形變場.其中StripMap模式像對覆蓋了形變中心，但是因?yàn)闂l紋分布密集和較大范圍失相干等原因，無法通過常規(guī)DInSAR提取出與真實(shí)情況相符的形變場，一些學(xué)者通常采用offset-tracking等方法獲得(Jolivet et al.，2014).另一組TerraSAR-X ScanSAR模式的像對其覆蓋范圍遠(yuǎn)離形變中心區(qū)域，形變量的顯著下降，Yague-Martinez等(2014)采取常規(guī)DInSAR法成功提取了遠(yuǎn)離形變中心東北部的同震形變場，并與Avouac等(2014)的同震模型進(jìn)行了對比，二者結(jié)果吻合較高.

本文在前期研究的成果上，首先通過選取典型DEM示范區(qū)，研究和討論子帶干涉處理過程中參數(shù)的選取問題，并嘗試總結(jié)出最優(yōu)化的參數(shù)方案.然后在此基礎(chǔ)上，將子帶干涉法應(yīng)用于2013年巴基斯坦地震形變中心附近(StripMap模式像對所覆蓋的范圍)同震形變場的獲取，力求通過該方法，更準(zhǔn)確地提取干涉條紋密集區(qū)的形變量.

2子帶干涉測量法

2.1基本原理

根據(jù)雷達(dá)干涉測量原理(李德仁等，2000)，干涉相位φ為

(1)

其中fc是雷達(dá)傳感器的中心頻率，ΔR是斜距差，c是光速.

為了計(jì)算高程，首先需要去除平地相位，去平后的干涉相位φ′與高程差Δh之間的關(guān)系(王超等，2002)為

(2)

其中R為主影像斜距，θ為入射角，B⊥為垂直基線距，λ為波長(λ=c/fc).

當(dāng)式(2)中φ′=2π時(shí)的高程差為Hamb=λRsinθ/(2B⊥)稱為模糊高.模糊高是指一個(gè)去平后的干涉相位周期所能代表的最大高程差.對于相同的高程差，越小的模糊高會(huì)產(chǎn)生越多的干涉條紋.在地形起伏劇烈的區(qū)域，較小的模糊高會(huì)產(chǎn)生過于密集的干涉條紋，進(jìn)而影響濾波和相位解纏的結(jié)果.由式(2)可知，模糊高與頻率成反比，與垂直基線距成反比.由于一組干涉像對其基線距是已定的，如需減少條紋數(shù)增加模糊高，只能通過增長波長實(shí)現(xiàn).

而對于地表形變的探測，λ/2稱為模糊形變周期，是一個(gè)去平并去地形相位后的干涉相位周期所能代表的最大形變差.與模糊高類似，模糊形變周期與波長成正比，對于大形變的區(qū)域，過小的模糊形變周期也會(huì)產(chǎn)生過于密集的干涉條紋，如需減少條紋數(shù)，也只有通過增長波長實(shí)現(xiàn).

子帶干涉法可通過模擬的手段獲得較長的波長.其原理是：

1) 首先利用帶通濾波器，將主影像和輔影像在距離向頻率域上，分割出上下兩個(gè)互不重疊的子帶像對，如圖1所示，其中上子帶和下子帶的距離向中心頻率分別為fup=fc+f0和flo=fc-f0，fc為原中心頻率，f0為偏移量，B為全頻帶寬，b為子帶帶寬.

圖1　距離向頻譜分割Fig.1　Range spectrum divided into upper sub-band and lower sub-band

2) 對上子帶像對和下子帶像對分別進(jìn)行干涉處理，根據(jù)式(1)相應(yīng)地可得到上子帶干涉相位φup和下子帶干涉相位φlo，公式為

(3)

(4)

3) 再將上子帶干涉相位與下子帶干涉相位進(jìn)行差分，得到子帶干涉相位Δφ，公式為

(5)

由圖1可知，f0的取值范圍在B以內(nèi)，按照一般SAR傳感器的設(shè)計(jì)B?fc，因此f0?fc.根據(jù)式(5)以及波長和頻率的關(guān)系(λ=c/f)可知，經(jīng)過子帶干涉處理后，極大地增長了波長，從而在高程計(jì)算或地表形變提取中極大地增加模糊高或模糊形變周期.相應(yīng)地就能減少甚至消除干涉圖中的條紋數(shù)量，使得相位解纏變得簡單甚至無需解纏而得到絕對相位差.

2.2最優(yōu)參數(shù)選取

選取用于提取DEM數(shù)據(jù)的影像對為TerraSAR-X Spotlight模式，覆蓋范圍為澳大利亞中部的北領(lǐng)地南部，影像覆蓋澳大利亞烏魯汝-卡塔楚塔國家公園(Uluru-Kata Tjuta National Park)中心和周邊地區(qū)，影像中心是一巨型的獨(dú)塊砂巖巖體，又被稱為艾爾斯巖(Young et al.,2002).ASTER GDEM 30 m分辨率的高程數(shù)據(jù)顯示該地區(qū)最大高程差約為357 m，該地區(qū)除艾爾斯巖體外，周邊是較為平坦的戈壁，像對獲取時(shí)間間隔僅11天，可認(rèn)為期間不存在地表形變，干涉相位主要由地形相位和平地相位組成，影像相干性非常好適合做處理流程和參數(shù)選取方面的討論.因此，采用逐步參數(shù)篩選方法，在DEM提取過程中，針對處理流程和參數(shù)的選取進(jìn)行討論，得出最優(yōu)化的方案.

2.2.1處理流程

圖2給出了子帶干涉法的基本處理流程：

圖2　子帶干涉法處理流程Fig.2　Processing flow of sub-band InSAR

1) 距離向頻譜分割：首先將主輔影像分別在距離向頻率域上分割出互不重疊的上下子帶，如果子帶之間頻率出現(xiàn)重疊，則重疊部分的能量會(huì)對結(jié)果產(chǎn)生嚴(yán)重的干擾.在帶通濾波器的選取上，為了集中主瓣能量，并能夠控制旁瓣能量快速衰減，本文采用了基于Kaiser窗函數(shù)的FIR濾波器(Kaiser，1974).

2) 子帶配準(zhǔn)：重采樣后的輔影像其頻譜會(huì)產(chǎn)生一定的改變，如果首先進(jìn)行配準(zhǔn)再進(jìn)行距離向頻譜分割，則有可能在某個(gè)頻段上引入不必要的誤差，表現(xiàn)在干涉圖背景上出現(xiàn)了一些周期性的線性條紋.因此，子帶配準(zhǔn)應(yīng)在距離向頻率分割之后進(jìn)行.

3) 子帶干涉：分別對上子帶和下子帶進(jìn)行干涉和去除地形相位等處理，得到上下子帶去平干涉相位后，再進(jìn)行一次干涉得到子帶干涉相位.

4) 相位平滑降噪：由于子帶帶寬只是全帶寬的一部分，而帶寬與分辨率成正比(Cumming and Wong，2005)，因此子帶干涉后的結(jié)果會(huì)降低分辨率，引入額外的噪聲并放大原有的噪聲，需要對其進(jìn)行一定程度的平滑降噪處理.本文采用了基于Goldstein的自適應(yīng)濾波器(Goldstein and Werner，1998)，采用兩次迭代以進(jìn)一步減小噪聲.

5) 相位解纏：對于一些高程差或形變量很大的區(qū)域，子帶干涉所產(chǎn)生的模糊高或模糊形變周期有可能仍無法完全涵蓋最大高程差或最大形變量并形成相位纏繞，因此仍需進(jìn)行解纏以求出絕對相位差.反之，如果模糊高或模糊形變周期已涵蓋最大高程差或最大形變量，子帶干涉相位差就是絕對值，這一步則可以省略.

2.2.2參數(shù)選取

子帶的中心頻率會(huì)影響波長，子帶帶寬的劃分則會(huì)影響精度，由于已有ASTER GDEM等數(shù)字高程模型作為先驗(yàn)參考，因此可設(shè)計(jì)一系列配置方案，研究這兩個(gè)參數(shù)對高程測量結(jié)果的影響程度，選取線性誤差和統(tǒng)計(jì)誤差作為指標(biāo)來對不同參數(shù)產(chǎn)生的結(jié)果進(jìn)行評價(jià)，從而實(shí)現(xiàn)最優(yōu)化參數(shù)的選取.如表1和圖3所示，本文以0.05B為間隔，從0.1B開始逐步增大子帶帶寬，fc為原SLC的中心頻率.同時(shí)，為了降低譜間干擾的影響，所有方案如圖3所示，均將上下子帶帶寬的間隔設(shè)置最大，子帶中心頻率均處于子帶帶寬中心.

線性誤差對比分析：首先將上述10組方案的干涉圖進(jìn)行平滑降噪，然后計(jì)算得到高程，并與ASTER GDEM高程數(shù)據(jù)做差分，得到高程殘差圖(圖4).可看出從方案1至方案10，在高程殘差圖上出現(xiàn)線性趨勢條紋的程度越為嚴(yán)重.分析原因，是因?yàn)閺姆桨?至方案10的子帶帶寬是不斷增大的，上下子帶帶寬的間隔則逐步縮小.帶通濾波器自身存在不可避免的旁瓣能量泄漏，會(huì)隨著間隔的逐步縮小對子帶干涉相位產(chǎn)生越為嚴(yán)重的譜間干擾，即出現(xiàn)距離向的線性趨勢條紋.而將這些相位條紋轉(zhuǎn)換為高程后，也成為高程上的誤差，反映在殘差圖上就是呈線性趨勢性的條紋誤差.

表1　子帶中心頻率和帶寬參數(shù)配置方案

統(tǒng)計(jì)誤差對比分析：我們采用了均方根誤差(root mean square error, RMSE)來檢測方案1至方案10計(jì)算得到的高程值相對于ASTER GDEM高程值的偏離程度.由方案1至方案10計(jì)算得到的RMSE分別為{47.5844, 29.6204, 23.0572, 24.3553, 23.3479, 25.7188, 23.8768, 28.9921, 32.4175, 37.7373}(m)，與各自對應(yīng)的子帶帶寬可繪制出如圖5所示的關(guān)系曲線.如該曲線所示，RMSE隨子帶帶寬的增大呈現(xiàn)先下降后上升的趨勢.其中最小值出現(xiàn)在0.2B的位置，對應(yīng)方案3.而在此之前方案1和方案2的RMSE均較大，是因?yàn)槠鋵?yīng)的子帶帶寬0.1B和0.15B過小，分辨率損失過大所帶來的相位模糊造成的，而旁瓣能量泄漏此時(shí)造成的影響在此處非常輕微，基本可以忽略.當(dāng)子帶帶寬大于0.2B且繼續(xù)增大時(shí)，對應(yīng)的RMSE在0.25B～0.35B處出現(xiàn)緩慢抬升，而當(dāng)>0.35B時(shí)則出現(xiàn)較快抬升.分析原因，雖然子帶帶寬的不斷增大會(huì)對分辨率以及干涉質(zhì)量(相干性)有所提升，但是也會(huì)使譜間干擾問題不斷加重，給測量結(jié)果帶來的負(fù)面影響甚至超過了分辨率提升帶來的增益.

圖3　子帶中心頻率和帶寬配置示意圖Fig.3　Sketch of central frequency and bandwidth of sub-band

圖4　方案1至方案10的高程殘差對比Fig.4　Comparison of elevation residuals among scheme 1 to scheme 10

圖5　方案1至方案10的均方根誤差對比Fig.5　Comparison of RMSE among scheme 1 to scheme 10

綜合上述分析，由于子帶帶寬和分辨率成正比，當(dāng)設(shè)置較小帶寬以避免譜間干擾影響的同時(shí)，會(huì)出現(xiàn)因分辨率降低過大而造成的模糊和噪聲放大等誤差，反之如果設(shè)置了較大的帶寬雖然能提高子帶影像的分辨率，卻會(huì)帶來譜間干擾的嚴(yán)重影響.因此，可看出譜間干擾的減少和分辨率的提高是“此消彼長”的相互矛盾的關(guān)系，無法同時(shí)兼顧.通過在真實(shí)數(shù)據(jù)基礎(chǔ)上的一系列實(shí)驗(yàn)，可得到子帶帶寬為0.2～0.3B，且上下子帶帶寬間隔最大時(shí)，分辨率可獲得較可靠的精度保證，同時(shí)受譜間干擾的影響也較為輕微，是理想的子帶帶寬選擇區(qū)間.

3巴基斯坦地震子帶干涉形變場獲取

3.1震區(qū)概況與數(shù)據(jù)選取

2013年9月24日，在巴基斯坦俾路支省(Balochistan)境內(nèi)的阿瓦蘭縣(Awaran)北北東66 km處發(fā)生了一次震級為MW7.7的地震，時(shí)隔4天后，又一個(gè)震級為MW6.6級的地震在該地區(qū)發(fā)生.馮光財(cái)?shù)?2015)和Avouac等(2014)，均采用交叉頻譜相關(guān)(cross-spectrum correlation)算法(Leprince et al.，2007b)，利用Landsat 8光學(xué)衛(wèi)星震前和震后影像，計(jì)算像對的偏移量獲取了此次地震大范圍的地表形變場.這些研究顯示，此次地震形變沿著地表破裂跡線呈弧形分布，地震類型為左旋走滑兼少許逆沖，破裂跡線北部造成的形變總體上大于南部.本文選取的像對為TerraSAR-X StripMap模式，影像覆蓋此次地震的形變中心，參數(shù)如表2所示.ASTER GDEM資料顯示，該區(qū)域在破裂跡線北部山嶺分布密集，多數(shù)山峰在800～1000 m以上，地形起伏情況相較于南部更為復(fù)雜.圖6b是二軌法差分干涉相位與強(qiáng)度圖的疊加，由于X波段的模糊形變周期只有1.6 cm，對形變和地形都十分敏感，易產(chǎn)生密集的干涉條紋或失相干.干涉條紋在分布形態(tài)上，破裂跡線南北兩側(cè)各有不同.南側(cè)以大量密集且不連續(xù)的干涉條紋為主(圖7-2)；在影像中部跡線南側(cè)位置，存在一處大面積的失相干帶，由于此處是居民地，可能受地震影響建筑物大面積倒塌，或是人類活動(dòng)改變了地表地面狀況而造成失相干.跡線北側(cè)則因過大的形變和復(fù)雜的地形影響造成了大范圍的失相干，只有影像西北角殘存了部分干涉條紋(圖7-1).

表2　阿瓦蘭地區(qū)TerraSAR-X StripMap影像參數(shù)

3.2子帶干涉形變場提取與驗(yàn)證

由Landsat 8提取的形變場顯示(圖6a)，本文選用的TerraSAR-X影像覆蓋范圍內(nèi)的地距向形變量為-5.64～2.88 m，對應(yīng)的斜距向形變約為-4.69～2.4 m.使用子帶干涉法需要模擬產(chǎn)生至少18.76 m的波長，即子帶中心頻率約為0.16B=16 MHz，才不會(huì)產(chǎn)生相位纏繞.由3.3.1節(jié)分析可知，這樣的參數(shù)配置會(huì)使上下子帶的中心頻率過于接近，帶來譜間干擾的影響，從而使結(jié)果存在較為嚴(yán)重的誤差.另一方面，常規(guī)二軌法差分干涉相位顯示(圖6b)，緊鄰破裂跡線區(qū)域由于形變量過大和地形復(fù)雜等因素，導(dǎo)致失相干.這些在常規(guī)InSAR中失相干區(qū)域的區(qū)域，其干涉條紋不可識別，干涉相位被嚴(yán)重的噪聲淹沒，即使經(jīng)過子帶干涉處理，其結(jié)果也是不正確的，而實(shí)際有干涉條紋覆蓋的區(qū)域，其斜距向形變量約為-4.0～2.4 m.綜合上述分析，由于此次地震產(chǎn)生的地表形變量遠(yuǎn)大于一般子帶干涉所能模擬的模糊形變周期，因此在子帶中心頻率和帶寬的配置上，應(yīng)以減少干涉條紋數(shù)為主要目標(biāo).

圖6　(a) Landsat 8提取的2013年巴基斯坦地震地表形變(馮光財(cái)?shù)龋?015)，已轉(zhuǎn)換為雷達(dá)地距方向；正值為接近衛(wèi)星方向; (b) TerraSAR-X StripMap常規(guī)DInSAR干涉相位與強(qiáng)度圖疊加；其覆蓋區(qū)域?yàn)閳D6a中綠色矩形框所示.圖6a和圖6b中的黑色實(shí)線為地表破裂跡線Fig.6　(a) The surface deformation of 2013 Pakistan earthquake from Landsat 8 satellite images (Feng, et al, 2015), which is transformed into the ground-range direction. Positive values mean close to the satellite; (b) The stacking chart of conventional DinSAR interferometric phase and magnitude of TerraSAR-X StripMap images. The green rectangular of Fig.6a indicates its coverage region. The black solid lines in Fig.6a and Fig.6b are surface rupture zone

圖7　圖6b的局部放大，其覆蓋區(qū)域?yàn)閳D6b中紅色矩形框所示Fig.7　The partial enlargement of Fig.6b. The red rectangular of Fig.6b indicates its coverage region

通過上文的一系列實(shí)驗(yàn)，已經(jīng)獲得了子帶參數(shù)配置的最優(yōu)方案.該方案雖然基于地形提取實(shí)驗(yàn)，但實(shí)驗(yàn)中討論的分辨率縮減和譜間干擾這兩方面問題實(shí)際是由圖像自身性質(zhì)變化造成的，與試驗(yàn)區(qū)的選取以及子帶干涉的用途無關(guān).換句話說，子帶干涉無論用途為何，其最終目的都是為了減少干涉條紋數(shù)以降低解纏難度，因此上文的最優(yōu)參數(shù)配置方案，也同樣適用于子帶干涉形變場的提取.綜合上述分析，本文將采用上文中子帶中心頻率為0.8B，子帶帶寬為0.2B的最優(yōu)方案，模擬產(chǎn)生的波長為3.75 m，模糊形變周期為1.87 m，提取該震例的形變場.

圖8a是經(jīng)過子帶干涉處理、5(距離向)×6(方位向)的多視處理，以及3次迭代自適應(yīng)濾波處理(濾波窗口分別為128、96、64)后得到的干涉相位，其干涉條紋分布在破裂跡線南側(cè)十分稀疏，相位在大部分地區(qū)處在一個(gè)形變周期范圍內(nèi)，只在靠近破裂跡線附近出現(xiàn)了少量的纏繞現(xiàn)象；北側(cè)的干涉條紋出現(xiàn)在影像西北角，并出現(xiàn)了明顯的相位纏繞.參考Landsat 8提取的形變場(圖6a)，破裂跡線南北兩側(cè)的形變量由內(nèi)向外逐步遞減，且兩側(cè)的形變運(yùn)動(dòng)方向相反，因此需要對破裂跡線南北兩側(cè)的相位分別解纏.圖8b是分別對斷層兩側(cè)子帶干涉相位采用MCF方法解纏，計(jì)算地距向形變，并拼接得到的結(jié)果，紅色圓點(diǎn)所在位置為解纏起始點(diǎn).

圖8　(a) 子帶干涉相位; (b) 分別對破裂跡線南北兩側(cè)的干涉相位解纏，并計(jì)算得到的地距向形變場；紅色圓點(diǎn)為解纏起始點(diǎn).圖8a和圖8b中的黑色實(shí)線為地表破裂跡線Fig.8　(a) Sub-band interferometric phase; (b) Phase unwrapping from north and south side of surface rupture zone respectively, and obtaining deformation field in slant direction. Red point is the beginning of phase unwrapping. Note: The black solid lines in Fig.8a and Fig.8b are surface rupture zone

圖9　不同手段獲得的同震形變場及殘差分布圖(a) 模型擬合得到的地表形變場； (b) 子帶干涉法提取的地表形變場，(低相干區(qū)已被掩膜去除)； (c) Landsat 8提取的地表形變場； (d) TerraSAR-X基于offset-tracking算法提取的地表形變場； (e) 常規(guī)DInSAR提取的地表形變場，(低相干區(qū)已被掩膜去除)； (f—i) 子帶干涉法、Landsat 8、offset-tracking法和常規(guī)DInSAR各自的的殘差分布.Fig.9　The coseismic deformation field and residuals distribution map obtained by different approaches.(a) Model prediction, (b) sub-band InSAR (low coherence area has been removed by mask), (c) Landsat 8, (d) TerraSAR-X based on offset-tracking, (e) conventional DInSAR (low coherence area has been removed by mask); (f—i) The residuals distribution of sub-band interferometric, Landsat 8, offset-tracking and conventional DInSAR.

圖10　剖面地距向形變量Fig.10　The deformation value in the profile ground-range direction

4形變場比較與驗(yàn)證

由于阿瓦蘭地震發(fā)生在偏遠(yuǎn)的山區(qū)，GPS臺(tái)站分布稀疏，地面調(diào)查資料缺乏，因此本文使用的基準(zhǔn)參考形變場來自馮光財(cái)?shù)?2015)的研究成果，其主要基于Landsat 8、USGS(美國地質(zhì)調(diào)查局)的觀測結(jié)果，經(jīng)Okada彈性半空間形變模型(Okada，1985)反演出此次地震的斷層滑動(dòng)分布，最后正演得到(圖9a).本文將這一結(jié)果分別與子帶干涉法、Landsat 8光學(xué)影像的交叉頻譜相關(guān)法、offset-tracking法和常規(guī)DInSAR各自提取的地表形變場(圖9b—e)進(jìn)行差分，得到各自的殘差圖(圖9f—i)，并分別在圖9a—e的相同位置從南至北繪制了3條橫跨斷層的剖面p1，p2，p3(圖10).

對比分析圖9和圖10，可看出以上4種方法均不同程度地獲得了此次地震的地表形變場.其中：

1) offset-tracking法和交叉頻譜相關(guān)法，都是一種基于強(qiáng)度信息的子像素精密匹配算法，適用于大形變量的探測(Gray et al., 1998)，offset-tracking法還可避免InSAR測量中低相干或失相干而造成無法解纏等問題(劉云華等，2012；李佳等，2013).如圖9c和圖9d所示，上述兩種方法相比子帶干涉法(圖9b)和常規(guī)DInSAR(圖9e)，均獲得了更完整的地表形變場和更清晰的地表破裂跡線，且跡線南北兩側(cè)的運(yùn)動(dòng)方向也與模型擬合的結(jié)果(圖9a)相吻合.但這兩種方法會(huì)受到地形、陰影和數(shù)據(jù)質(zhì)量等因素的影響，其精度一般要比InSAR低，offset-tracking法的理論精度在分米級(Strozzi et al.，2002)，而交叉頻譜相干法的理論最高精度為像元大小的1/50(Leprince et al.，2007a)，另外二者均存在細(xì)節(jié)分辨力偏低和噪聲偏大等問題.圖9c和圖9d顯示，在破裂跡線南側(cè)，存在大量的噪聲異常值；與之對應(yīng)的殘差分布(圖9g和圖9h)上，可看到北側(cè)在緊鄰跡線的區(qū)域形變值偏大；在遠(yuǎn)離破裂跡線的北側(cè)區(qū)域，Landsat 8表現(xiàn)出總體吻合局部偏大的特點(diǎn)，offset-tracking則呈現(xiàn)由近到遠(yuǎn)先偏大后偏小的趨勢.從剖面中可看出，上述2種方法均存在形變值在各處均呈現(xiàn)波動(dòng)較大的現(xiàn)象，說明受噪聲干擾較為嚴(yán)重.

2) 子帶干涉法和常規(guī)DInSAR法，受到影像中2條失相干帶的影響，無法計(jì)算出正確的相位，進(jìn)而在最終結(jié)果中(圖9b和圖9e)該處的形變值為幅值抖動(dòng)劇烈的噪聲，故本文已將其掩膜去除.如圖9f、圖9i和圖10所示，2條失相干帶將影像分為南部、中部和西北部3塊相對獨(dú)立的區(qū)域(分別在圖9b和圖9e中以I、II和III標(biāo)出)，其中：

I號區(qū)域，位于影像南部是遠(yuǎn)離形變中心的區(qū)域，上述2種方法在此處獲得的形變場均與模型擬合結(jié)果擬合較好，其平均誤差分別在5 cm和8 cm的范圍內(nèi)；剖面線(圖10)顯示與模型擬合結(jié)果相吻合，且噪聲相較于offset-tracking法和交叉頻譜相關(guān)法少，形變趨勢變化平穩(wěn).

II號區(qū)域，在影像中部，破裂跡線的南側(cè)，是形變劇烈的區(qū)域.在圖9b—d中黑色實(shí)線矩形框標(biāo)出了一處形變值急劇增大區(qū)域，對應(yīng)圖10b剖面線可看出相較于模型擬合結(jié)果，形變值發(fā)生了急劇的抬升，在殘差圖上，此處出現(xiàn)了較大的負(fù)值誤差，而常規(guī)DInSAR的結(jié)果(圖9e)在此處與模型擬合結(jié)果吻合較好.本文認(rèn)為，上述出現(xiàn)較大誤差的3種方法的結(jié)果，具有一定的一致性，且是由不同數(shù)據(jù)源和不同處理手段得到的，因此這一形變值劇烈增大的區(qū)域可能客觀存在.如果這一假設(shè)成立，則模型擬合結(jié)果此處形變值變化平緩，可能與四叉樹降采樣的密度設(shè)置有關(guān)，而DInSAR的結(jié)果也未出現(xiàn)形變值劇增的趨勢，則可能是因?yàn)閷γ芗瘲l紋做解纏時(shí)，因噪聲和低相干性的影響，造成條紋混疊并丟失導(dǎo)致.

III號區(qū)域，位于影像西北部，由模型擬合得到的形變場(圖9a)及對應(yīng)的剖面顯示(圖10)，該地區(qū)形變量約為-1.4～-3.7 m，由南至北基本呈線性遞減，遞減速率(斜率)約為0.15 m·km-1.Landsat 8和offset-tracking提取的形變場在剖面上，與該結(jié)果基本吻合，但存在振幅波動(dòng)較大的噪聲.而常規(guī)DInSAR計(jì)算得到的形變，與其他方法得到的形變差異很大，這些差異同時(shí)反映在形變速率以及形變值分布這2個(gè)方面：剖面顯示，其形變遞減速率十分平緩只有2.6 cm·km-1；形變分布(圖9e)顯示，其計(jì)算得到的地距向形變值約為-0.15～-0.36 cm，與模型擬合結(jié)果的誤差很大(圖9i).造成上述2方面誤差的原因，也與該處干涉條紋密集，解纏時(shí)條紋丟失有關(guān).對于子帶干涉，在該處只發(fā)生了一個(gè)周期的相位纏繞，可十分容易地.解纏得到更接近真實(shí)情況的絕對相位差.圖9b所示，由子帶干涉獲取的西北部的地距向形變值約為-0.87～-3.77 m與模型擬合結(jié)果相接近.殘差分布(圖9f)顯示，殘差由北至南，由形變遠(yuǎn)場至近場呈現(xiàn)遞減的趨勢；剖面(圖10)顯示，最遠(yuǎn)端與模型擬合結(jié)果的誤差最大達(dá)到0.7 m，隨著逐步靠近形變近場，誤差逐步控制在0.5～0.2 m或更低的水平，相較于offset-tracking和Landsat 8提取的形變值的平均誤差(分別為1.5 m和1.8 m)更小，且趨勢平穩(wěn)增大，噪聲水平更低.我們注意到，子帶干涉法在此處出現(xiàn)了相較于影像南部更大的誤差，其原因是因?yàn)椴捎昧顺Ｒ?guī)的多項(xiàng)式全局配準(zhǔn)策略，在處理影像局部存在的大偏移、復(fù)雜地形情況時(shí)出現(xiàn)的局部偏差，這一配準(zhǔn)處理方法上的局限，在常規(guī)InSAR中也是存在且難以避免的.

5結(jié)論

子帶干涉法作為一種傳統(tǒng)InSAR技術(shù)的衍生，其減少干涉條紋數(shù)的特性，可有效避免對條紋密集區(qū)域進(jìn)行解纏出現(xiàn)的誤差；也可避免因絕對相位差過大而待解纏區(qū)域過小，造成條紋丟失，產(chǎn)生錯(cuò)誤的解纏結(jié)果.因此可直接用于測量，也可作為一種手段用于檢驗(yàn)傳統(tǒng)InSAR解纏結(jié)果的準(zhǔn)確性.本文通過子帶干涉法以提取澳大利亞艾爾斯巖區(qū)域的DEM為基礎(chǔ)，分析了子帶中心頻率和帶寬參數(shù)的選取，并將最優(yōu)化方案應(yīng)用于2013年巴基斯坦地震的同震地表形變場提取中，通過對這實(shí)例數(shù)據(jù)進(jìn)行處理和對比分析，得出如下結(jié)論：

1) 子帶中心頻率和帶寬參數(shù)配置的不同會(huì)對最終結(jié)果產(chǎn)生較大的影響.在保證波長足夠長的前提下，應(yīng)選擇高相干點(diǎn)分布比例大的配置方案，可保證干涉圖的總體質(zhì)量；應(yīng)將上下子帶的中心頻率和帶寬的間隔設(shè)置盡量大，可減少譜間干擾產(chǎn)生的偽信號混淆測量結(jié)果以提高測量精度.

2) 為避免子帶干涉圖上產(chǎn)生斜向條紋影響測量結(jié)果，應(yīng)首先對主輔影像做距離向頻率分割，生成上下子帶影像，再分別對其進(jìn)行配準(zhǔn)和重采樣.在對子帶干涉圖進(jìn)行降噪平滑時(shí)，可將搜索窗口適當(dāng)增大，并進(jìn)行二次或更多次的迭代濾波，以達(dá)到更好的降噪效果.

3) 同常規(guī)InSAR一樣，子帶干涉法也依賴于相干性.常規(guī)InSAR中失相干的區(qū)域，在子帶干涉中該區(qū)域的干涉相位也可能是不準(zhǔn)確.因此，可利用常規(guī)InSAR的相干性分布圖，結(jié)合其他觀測數(shù)據(jù)，預(yù)先估計(jì)出影像中有效干涉相位所覆蓋的區(qū)域內(nèi)的值域范圍，以選取合適的子帶中心頻率和帶寬等參數(shù).在對常規(guī)InSAR干涉條紋密集的區(qū)域做解纏時(shí)，易發(fā)生條紋丟失使參與計(jì)算的條紋數(shù)低于實(shí)際情況，從而造成較大的誤差.而子帶干涉法能有效減少或消除干涉條紋數(shù)，降低解纏難度甚至無需解纏，可有效地避免條紋丟失問題，獲得更為準(zhǔn)確的結(jié)果.

4) 通過對2013年巴基斯坦震例進(jìn)行的研究，可看出子帶干涉法在提取大尺度地表形變上具有很大的優(yōu)勢.首先，雖然受失相干的影響其提取的形變場范圍相較于交叉頻譜相關(guān)法和offset-tracking法有所缺失，但在共同覆蓋的區(qū)域，子帶干涉法所測量的結(jié)果均顯示出了更高的精度和更低的噪聲水平.其次，在DInSAR存在干涉條紋且十分密集的區(qū)域(如影像中部區(qū)域)，子帶干涉法不但能更為準(zhǔn)確地獲得絕對相位差，而且與交叉頻譜相關(guān)法或offset-tracking法相比，其反映的形變場細(xì)節(jié)更豐富.

致謝感謝夏耶教授生前對本文完成給予的十分寶貴的幫助.

References

Avouac J P, Ayoub F, Wei S J, et al. 2014. The 2013,MW7.7 Balochistan earthquake, energetic strike-slip reactivation of a thrust fault.EarthandPlanetaryScienceLetters, 391: 128-134.

Brcic R, Eineder M, Bamler R. 2008. Absolute phase estimation from TerraSAR-X acquisitions using wideband interferometry.∥Proceedings of IEEE Radar Conference. Pasadena, CA: IEEE.

Costantini M. 1998. A novel phase unwrapping method based on network programming.IEEETransactionsonGeoscienceandRemoteSensing, 36(3): 813-821.

Cumming I G, Wong F H. 2005. Digital Processing of Synthetic Aperture Radar Data: Algorithms and Implementation. Boston: Artech House. Derauw D, Orban A, Barbier C. 2010. Wide band SAR sub-band splitting and inter-band coherence measurements.RemoteSensingLetters, 1(3): 133-140.

Du Y N, Feng G C, Li Z W et al .2015.Generation of high precision DEM from TerraSAR-X/TanDEM-X.ChineseJournalofGeophysics(in Chinese),58(9): 3089-3102,doi: 10.6038/cjg20150907.Feng G C, Xu B, Shan X J, et al. 2015. Coseismic deformation and source parameters of the 24 September 2013 Awaran, PakistanMW7.7 Earthquake derived from optical Landsat 8 satellite images.ChineseJournalofGeophysics(in Chinese), 58(5): 1634-1644, doi: 10.6038/cjg20150515.

Goldstein R M, Zebker H A, Werner C L. 1988. Satellite radar interferometry: Two-dimensional phase unwrapping.RadioScience, 23(4): 713-720.

Goldstein R M, Werner C L. 1998. Radar interferogram filtering for geophysical applications.GeophysicalResearchLetters, 25(21): 4035-4038. Gray A L, Mattar K E, Vachon P W, et al. 1998. InSAR results from the RADARSAT Antarctic Mapping Mission data: estimation of glacier motion using a simple registration procedure. 1998 IEEE International Geoscience and Remote Sensing Symposium, 1998, 3: 1638-1640.Jolivet R, Duputel Z, Riel B, et al. 2014. The 2013MW7.7 Balochistan Earthquake: Seismic potential of an accretionary wedge.BulletinoftheSeismologicalSocietyofAmerica, 104(2): 1020-1030.

Kaiser J F. 1974. Nonrecursive digital filter design using theI0-Sinh window function.∥Proceedings of the 1974 IEEE International Symposium on Circuits and Systems. New York: IEEE, 20-23. Leprince S, Barbot S, Ayoub F, et al. 2007a. Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements.IEEETransactionsonGeoscienceandRemoteSensing, 45(6): 1529-1558. Leprince S, Ayoub F, Klingert Y, et al. 2007b. Co-registration of optically sensed images and correlation (COSI-Corr): An operational methodology for ground deformation measurements.∥Proceedings of IEEE International Geoscience and Remote Sensing Symposium. Barcelona: IEEE, 1943-1946.

Li D R, Zhou Y Q, Ma H C. 2000. Principles and applications of interferometry SAR.ScienceofSurveyingandMapping(in Chinese), 25(1): 9-12.

Li J, Li Z W, Wang C C et al .2013.Using SAR offset-tracking approach to estimate surface motion of the South Inylchek Glacier in Tianshan.ChineseJournalofGeophysics,56(4): 1226-1236,doi: 10.6038/cjg20130417.

Li Y S, Feng W P, Zhang J F，et al. 2015.Coseismic slip of the 2014MW6.1 Napa, California Earthquake revealed by Sentinel-1A InSAR.ChineseJournalofGeophysics(in Chinese),58(7): 2339-2349,doi: 10.6038/cjg20150712.

Liao M S, Lin H. 2003. Synthetic Aperture Radar Interferometry: Principle and Signal Processing (in Chinese). Beijing: Surveying and Mapping Press.Liu Y H, Qu C Y, Shan X J. 2012.Two-dimensional displacement field of the Wenchuan earthquake inferred from SAR intensity offset-tacking.ChineseJournalofGeophysics(in Chinese),55(10): 3296-3306,doi: 10.6038/j.issn.0001-5733.2012.10.012.

Liu Y H, Wang C S, Shan X J, et al. 2014.Result of SAR differential interferometry for the co-seismic deformation and source parameter of theMS7.0 Lushan Earthquake.ChineseJournalofGeophysics(in Chinese), 57(8): 2495-2506,doi: 10.6038/cjg20140811.

Madsen S N, Zebker H A, Martin J. 1993. Topographic mapping using radar interferometry: Processing techniques.IEEETransactionsonGeoscienceandRemoteSensing, 31(1): 246-256.

Okada Y. 1985. Surface deformation due to shear and tensile faults in a half-space.BulletinoftheSeismologicalSocietyofAmerica, 75(4): 1135-1154.

Shan X J, Zhang G H, Wang C S, et al. 2015.Joint inversion for the spatial fault slip distribution of the 2015 NepalMW7.9 earthquake based on InSAR and GPS observations.ChineseJournalofGeophysics(in Chinese),58(11): 4266-4276,doi: 10.6038/cjg20151131.Strozzi T, Luckman A, Murray T, et al. 2002. Glacier motion estimation using SAR offset-tracking procedures.IEEETransactionsonGeoscienceandRemoteSensing, 40(11): 2384-2391. Wan Y G. 2012. Digital Signal Processing-Using MATLAB (in Chinese). 2nd ed. Beijing: Science Press.

Wang C, Zhang H, Liu Z. 2002. Spaceborne Synthetic Aperture Radar Interferometry (in Chinese). Beijing: Science Press.

Wang X P. 2010. Minimum cost flow and its ameliorative algorithm for phase unwrapping of InSAR image.ScienceofSurveyingandMapping(in Chinese), 35(4): 129-131.

Yague-Martinez N, Fielding E, Haghshenas-Haghighi M, et al. 2014. Ground displacement measurement of the 2013M7.7 and

M6.8 Balochistan earthquake with TerraSAR-X ScanSAR data.∥Proceedings of IEEE International Geoscience and Remote Sensing Symposium. Quebec City, QC: IEEE, 950-953.Young D N, Duncan N, Camacho A, et al. 2002. Ayers Rock, Northern Territory, Map Sheet GS52-8 (2nd edition) (Map). 1:250 000. Northern Territory Geological Survey. Geological Map Series Explanatory Notes.

Zhan W J, Li Z W, Wei J C, et al. 2015.A strategy for modeling and estimating atmospheric phase of SAR interferogram.ChineseJournalofGeophysics(in Chinese), 58(7): 2320-2329,doi: 10.6038/cjg20150710.

附中文參考文獻(xiàn)

杜亞男, 馮光財(cái), 李志偉等. 2015.TerraSAR-X/TanDEM-X獲取高精度數(shù)字高程模型技術(shù)研究. 地球物理學(xué)報(bào),58(9): 3089-3102,doi: 10.6038/cjg20150907.

馮光財(cái), 許兵, 單新建等. 2015. 基于Landsat 8光學(xué)影像的巴基斯坦AwaranMW7.7地震形變監(jiān)測及參數(shù)反演研究. 地球物理學(xué)報(bào), 58(5): 1634-1644, doi: 10.6038/cjg20150515.

李德仁, 周月琴, 馬洪超. 2000. 衛(wèi)星雷達(dá)干涉測量原理與應(yīng)用. 測繪科學(xué), 25(1): 9-12.

李佳, 李志偉, 汪長城等. 2013.SAR偏移量跟蹤技術(shù)估計(jì)天山南依內(nèi)里切克冰川運(yùn)動(dòng). 地球物理學(xué)報(bào),56(4): 1226-1236,doi: 10.6038/cjg20130417.

李永生, 馮萬鵬, 張景發(fā)等. 2015. 2014年美國加州納帕MW6.1地震斷層參數(shù)的 Sentinel-1A InSAR反演. 地球物理學(xué)報(bào), 58(7): 2339-2349, doi: 10.6038/cjg20150712.

廖明生, 林琿. 2003. 雷達(dá)干涉測量: 原理與信號處理基礎(chǔ). 北京: 測繪出版社.

劉云華, 屈春燕, 單新建 .2012.基于SAR影像偏移量獲取汶川地震二維形變場. 地球物理學(xué)報(bào),55(10): 3296-3306,doi: 10.6038/j.issn.0001-5733.2012.10.012.

劉云華, 汪馳升, 單新建等. 2014.蘆山MS7.0級地震InSAR形變觀測及震源參數(shù)反演. 地球物理學(xué)報(bào),57(8): 2495-2506,doi: 10.6038/cjg20140811.

單新建，張國宏，汪馳升等. 2015. 基于InSAR和GPS觀測數(shù)據(jù)的尼泊爾地震發(fā)震斷層特征參數(shù)聯(lián)合反演研究. 地球物理學(xué)報(bào), 58(11): 4266-4276, doi: 10.6038/cjg20151131.

占文俊, 李志偉, 韋建超等. 2015.一種InSAR大氣相位建模與估計(jì)方法. 地球物理學(xué)報(bào),58(7): 2320-2329,doi: 10.6038/cjg20150710.

(本文編輯劉少華)

Deformation of the 2013 PakistanMW7.7 earthquake derived from sub-band InSAR

YU Lu1,2， SHAN Xin-Jian1*， SONG Xiao-Gang1， QU Chun-Yan1

1StateKeyLaboratoryofEarthquakeDynamics,InstituteofGeology,ChinaEarthquakeAdministration,Beijing100029,China2GuangxiRemoteSensingCenter,Nanning530023,China

AbstractOn September 24, 2013, an MW7.7 earthquake strok Awaran, Balochistan Province, Pakistan, causing large strike slip up to 10 m. The coseismic deformation field, produced by TerraSAR-X short wave radar data, shows dense and extensive interference fringes which bring difficulties to thephase unwrapping. The sub-band InSAR is a new method for obtaining absolute phase, with a little or without phase unwrapping. This method shortens the bandwidth to increase the wavelength, and could reduce the number of fringes, which decreases the difficulty of phase unwrapping. However the noise and the extra interference of sidelobes will increase and the interferogram quality will decline, while the bandwidth reduces. More consideration should be taken into the parameter selection, noise filtering and processing flow, especially the selection of center frequency and bandwidth. In this paper, we first select the typical DEM experiment zone as the research area, and study the effect of parameter selection on sub-band InSAR with coherences and elevation errors as the evaluation index. Then, we apply the optimal parameter scheme to derive the coseismic deformation field of Pakistan. The coseismic deformation, extracted by sub-band interference, Landsat 8 optical images cross-spectrum correlation, offset-tracking, and conventional DInSAR separately, are analyzed and compared with the deformation field from modeling. It indicates that sub-band InSAR performs batter in precision and noise level, in spite of the slight influence of decorrelation and loss of the extracted deformation field, compared with Landsat-8 and offset tracking. Therefore this method is more suitable to be applied in areas with large deformation and dense fringes.

KeywordsSub-band InSAR; Pakistan earthquake; Coseismic deformation

基金項(xiàng)目國家自然科學(xué)基金(41461164002)、中國地震局地質(zhì)研究所所長基金(IGCEA1404)、地震動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室自主研究課題(LED2013A02)以及中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金(14CX02110A)聯(lián)合資助.

作者簡介庾露，男，1984年生，工程師，主要從事InSAR算法方面研究.E-mail:yulu_enter@163.com *通信作者單新建，男，1966 年生，研究員，主要從事地殼形變觀測與動(dòng)力學(xué)研究.E-mail:xjshan@163.com

doi:10.6038/cjg20160418 中圖分類號P315

收稿日期2015-11-11，2016-02-20收修定稿

庾露， 單新建， 宋小剛等. 2016. 基于子帶干涉測量技術(shù)的巴基斯坦地震形變獲取研究.地球物理學(xué)報(bào)，59(4):1371-1382,doi:10.6038/cjg20160418.

Yu L, Shan X J, Song X G, et al. 2016. Deformation of the 2013 PakistanMW7.7 earthquake derived from sub-band InSAR.ChineseJ.Geophys. (in Chinese),59(4):1371-1382,doi:10.6038/cjg20160418.