基于DInSAR技術(shù)與斷層自動剖分方法反演斷層滑動參數(shù)

陳丹蕾， 劉國祥， 王曉文， 王蕾， 蒲慧龍

(1.四川省第三測繪工程院，成都　610500； 2.西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，成都　610031)

?

基于DInSAR技術(shù)與斷層自動剖分方法反演斷層滑動參數(shù)

陳丹蕾1， 劉國祥2， 王曉文2， 王蕾1， 蒲慧龍1

(1.四川省第三測繪工程院，成都610500； 2.西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，成都610031)

摘要：以合成孔徑雷達差分干涉測量(DInSAR)技術(shù)獲得的同震形變數(shù)據(jù)為約束，基于彈性位錯模型反演獲得發(fā)震斷層的幾何、運動參數(shù)(也稱滑動參數(shù))，對于了解震區(qū)斷層活動特征和評估潛在災(zāi)害風(fēng)險具有重要意義?，F(xiàn)有的斷層滑動參數(shù)反演多是基于矩形位錯元的均勻剖分或人為設(shè)定剖分，所得結(jié)果易出現(xiàn)“偽滑動”和滑動分布過于平滑的問題，不足以精確地反映斷層面上的滑動細節(jié)。為此，引入斷層自動剖分方法。該方法兼顧了模型殘差最小化與平滑最優(yōu)原則，在形變數(shù)據(jù)的有效約束下獲得可靠解。以巴姆地震為例反演獲得斷層的幾何參數(shù)與滑動參數(shù)。實驗結(jié)果表明，對于單斷層的滑動參數(shù)反演，應(yīng)用斷層自動剖分方法可獲得較可靠的結(jié)果。

關(guān)鍵詞：斷層自動剖分；DInSAR； 滑動參數(shù)反演； 混合采樣

0引言

利用合成孔徑雷達干涉測量技術(shù)(interferometricsyntheticapertureRadar,InSAR)可以快速獲取震區(qū)地表連續(xù)覆蓋的同震形變信息[1-4]，通過在InSAR形變數(shù)據(jù)與地下斷層參數(shù)之間建立模型并進行解算，可獲得發(fā)震斷層的幾何參數(shù)和滑動參數(shù)，這對于研究地震的發(fā)震機理和了解震區(qū)的斷層活動特征具有重要意義。斷層的滑動參數(shù)反映了斷層的運動情況，為了構(gòu)建精細的斷層滑動模型，可將斷層剖分為許多個小位錯元，并細化至求解每個位錯元上的滑動值。常規(guī)的滑動參數(shù)反演大都基于矩形面源模型[2-4]，采用均勻剖分[5]和人為設(shè)定剖分[6]。然而，矩形位錯元在模擬復(fù)雜幾何模型時易在斷層連接處出現(xiàn)間隙和重疊現(xiàn)象，并且很難實現(xiàn)連續(xù)光滑曲面斷層的模擬。若使用三角形位錯元則可避免上述問題，實現(xiàn)任意形狀曲面的模擬，適用于復(fù)雜斷層[7-8]。在對斷層面進行離散剖分時，若采用均勻剖分或人為設(shè)定剖分反演得到的滑動分布常會出現(xiàn)滑動分布過于平滑，無法表現(xiàn)真實的精細分布特征； 或者在某些位置上的滑動分布太過細化，呈現(xiàn)出一些基于數(shù)據(jù)本身無法完全確定的滑動(偽滑動)特征[7]。

針對以上問題，引入斷層自動剖分方法[9]，基于尺寸可變的三角形位錯單元，以反演模型分辨率最優(yōu)為原則迭代得到最佳斷層位錯元尺寸，將滑動合理地分布到所有位錯元上，并以巴姆地震為例，應(yīng)用斷層自動剖分方法獲得斷層的滑動分布解。

1斷層自動剖分原理

1.1模型參數(shù)求解

依據(jù)地表形變數(shù)據(jù)與斷層參數(shù)之間的關(guān)系，斷層參數(shù)的反演可分為幾何參數(shù)的非線性反演和滑動參數(shù)的線性反演，統(tǒng)一用式(1)表示，即

G(m)=d ,

(1)

式中：G為格林函數(shù)矩陣，由均勻彈性半空間位錯模型[10]獲得； m為斷層參數(shù)； d為地表形變觀測值。

在反演求解參數(shù)的過程中常采用正則化算法，即通過引入附加信息來解決方程病態(tài)問題或避免過度擬合現(xiàn)象的產(chǎn)生[11]，這些附加信息通常是平滑約束信息或解空間的邊界約束條件。本文采用高階Tikhonov正則化算法求解m，權(quán)衡考慮了誤差最小化和解的平滑性，通過附加正則化條件的殘差二范式，即

(2)

來獲得殘差最小時對應(yīng)的模型解m。式中：L為拉普拉斯平滑矩陣，用以克服鄰近位錯元上的劇烈滑動梯度變化，避免解振蕩的產(chǎn)生； λ是正則化權(quán)參數(shù)，若其值過大，會濾掉很多有用的觀測值，造成模型的過度平滑; 其值過小，則會混入過多噪聲，致使擬合模型中混入不可靠信息，偏離真實情況。一個合適的正則化權(quán)參數(shù)既能平滑噪聲，又能避免模型的過度擬合。本文采用jRi的方法[9]獲得正則化權(quán)參數(shù)，jRi的值反映了反演質(zhì)量的優(yōu)劣,其值越小，代表當(dāng)下斷層剖分越合理，取jRi最小值賦給λ。jRi的表達式為

(3)

式中： Cd為觀測值的方差—協(xié)方差矩陣；n為觀測數(shù)。求得λ值后便可計算廣義逆矩陣G-g，即

(4)

在已知G，λ和L的前提下，可利用正則化方法以最小二乘為約束，求得m的初始分布模型。

1.2自動剖分

斷層自動剖分方法是利用模型分辨率矩陣擬合產(chǎn)生所有位錯元的最佳尺寸，從而對初始位錯元的大小進行調(diào)節(jié)，使滑動在位錯元上合理分布的過程。自動剖分流程如圖1所示。

圖1　斷層自動剖分流程

斷層自動剖分的具體步驟包括數(shù)據(jù)的加載、斷層的初始剖分、正則化格林函數(shù)和模型分辨率矩陣的計算、最佳位錯尺寸的確定、斷層優(yōu)化剖分和迭代收斂檢驗。

1.2.1模型分辨率矩陣的計算

模型分辨率矩陣R是與位錯中心位置和地表形變觀測數(shù)據(jù)質(zhì)量有關(guān)的一個參數(shù)矩陣，它的行元素代表凈滑動值在相關(guān)斷層片中的平滑程度[12]，可以由

R=G-gG

(5)

予以確定。根據(jù)矩陣R可以判斷位錯元上的滑動是偽滑動，還是不足以表現(xiàn)細節(jié)的粗略滑動分布，繼而通過產(chǎn)生新的位錯剖分來更好地反映數(shù)據(jù)的解析能力。

1.2.2最佳位錯尺寸的確定

初始滑動分布確定后，需要基于數(shù)據(jù)的解析能力，對位錯尺寸進行可變剖分。為了使斷層模型中每個位錯元上的滑動都是基于數(shù)據(jù)的有力約束，需要在位錯尺寸與滑動分布的平滑尺度之間建立函數(shù)對應(yīng)關(guān)系，通過高斯曲線擬合，建立模型分辨率矩陣R的第i行元素(表示第i個位錯元上的滑動在所有位錯元上的平滑程度)與位錯尺寸之間的高斯函數(shù)關(guān)系，如圖2所示。

圖2　高斯擬合獲得最佳位錯尺寸

擬合得到的高斯函數(shù)的寬度設(shè)定為平滑尺度，將其大小作為位錯的尺寸，此時滑動在各個位錯元上的分布接近狄克拉函數(shù)，對應(yīng)的斷層剖分是最佳剖分。比較平滑尺度與位錯尺寸，如果平滑尺度接近當(dāng)前位錯元的大小，該位錯元的剖分就比較適宜； 如果平滑尺度大于或小于位錯尺寸，那么位錯元的剖分就不正確，需要對位錯尺寸進行調(diào)節(jié)。

1.2.3迭代終止條件

以相鄰2次迭代所得位錯元總數(shù)的差異百分數(shù)作為標準，即如果在某次迭代后，新增加或減少的位錯元數(shù)目的百分比在容許范圍內(nèi)，則終止迭代，否則將繼續(xù)進行迭代計算，直至滿足迭代條件為止。

2同震形變場獲取與數(shù)據(jù)降采樣

2.1同震形變場的獲取

利用SAR影像對的干涉處理結(jié)果可獲取巴姆地震的同震形變信息。本實驗選用2幅覆蓋同一震區(qū)范圍的降軌影像組成干涉對，其參數(shù)如表1所示。

表1　干涉對參數(shù)及基線估計

由主、從影像對干涉處理及差分干涉處理得到震區(qū)干涉相位圖及同震形變圖(圖3和圖4)。

圖3　巴姆地震干涉相位圖

圖4　巴姆地震同震形變圖

圖4的雷達視線向形變圖中一個顏色周期代表了2π的相位變化，對應(yīng)的地表形變值為2.8 cm，整個震區(qū)的形變范圍在-16.82～30.33 cm之間。地表形變集中分布在巴姆城市的東北側(cè)和東南側(cè)，東北側(cè)以地表沉降為主，東南側(cè)以地面抬升為主。西北側(cè)有2～4 cm的下降，西南側(cè)有少量抬升，其值在2 cm左右。

2.2數(shù)據(jù)降采樣

為保證反演迭代收斂和提高解算速率，需要對數(shù)量龐大的形變數(shù)據(jù)進行降采樣。鑒于四叉樹降采樣法可以較好地保留原始數(shù)據(jù)的分布特征[3-4]，本文首先利用四叉樹采樣方法對原始形變數(shù)據(jù)進行降采樣，其結(jié)果如圖5所示。

(a) 四叉樹采樣后的形變場 (b) 采樣點分布 (c) 殘差分布

圖5四叉樹采樣結(jié)果

Fig.5Results of quadtree resampling

由圖5可以看出，四叉樹采樣方法雖然可以較好地保持原始數(shù)據(jù)的分布特征，但在影像西南側(cè)受山體疊掩引起的失相干部分也提取到了相當(dāng)數(shù)量的采樣點，若讓這些數(shù)據(jù)參與后續(xù)的斷層參數(shù)反演將會影響反演精度?；谝陨峡紤]，本文采取混合采樣方法，對西南側(cè)受山體疊掩影響的遠場區(qū)采用均勻采樣，而對形變比較劇烈的近場區(qū)采用四叉樹降采樣。遠場區(qū)和近場區(qū)是以整體形變量絕對值的平均作為閾值進行劃分的。這樣的采樣方式既可以滿足保留原始數(shù)據(jù)特征的要求，又可以避免遠場區(qū)精度較差數(shù)據(jù)的干擾。

采用混合采樣方法，將形變數(shù)據(jù)由6 236 235個降采樣至1 736個，這些樣點密集地分布在原始形變圖中形變較大的區(qū)域，而其他區(qū)域的數(shù)據(jù)分布均勻且相對稀少，采樣點分布情況如圖6(b)所示。降采樣前后，殘差較大的地方出現(xiàn)在圖像的頂部位置，誤差約為4～5 cm，由于其范圍較小不會對反演結(jié)果產(chǎn)生太大影響，故可予以忽略； 其他地方數(shù)據(jù)的誤差都在0值附近，經(jīng)計算殘差的均方根誤差為1.06 cm，殘差分布如圖6(c)所示。

(a) 混合采樣后的形變場(b) 采樣點分布 (c) 殘差分布

圖6混合采樣結(jié)果

Fig.6Results of mixed resampling

3反演實驗

3.1斷層幾何參數(shù)的非線性反演

幾何參數(shù)的反演是滑動參數(shù)反演的基礎(chǔ)。本文對巴姆地震斷層的幾何參數(shù)反演是基于Okada[13]彈性半空間位錯理論，利用Levenberg-Marquardt(LM)最小二乘優(yōu)化算法迭代產(chǎn)生最優(yōu)解。Okada彈性半空間位錯理論為模型提供了聯(lián)系待求參數(shù)與InSAR形變數(shù)據(jù)之間的格林函數(shù)，而LM優(yōu)化算法是在最小二乘算法基礎(chǔ)上附加阻尼因子，從而在解空間中搜索得到殘差最小時的幾何參數(shù)最優(yōu)解，即尋求真實形變量與模型解算的模擬形變量之間的殘差最小時對應(yīng)的模型參數(shù)解。得到的幾何參數(shù)作為滑動分布反演部分的輸入?yún)?shù)，結(jié)果如表2所示。

表2　斷層幾何參數(shù)反演結(jié)果

根據(jù)所得斷層幾何參數(shù)反演結(jié)果，并結(jié)合已有資料，初步確定巴姆地震是由隱伏單斷層的右旋走滑運動(滑動角接近180°)為主引起的。為進一步求得斷層面上的精細滑動分布情況，需要將斷層面離散剖分為三角形位錯元，并求得每個位錯元上的滑動分布情況，本文采用斷層自動剖分方法實現(xiàn)。

3.2斷層滑動參數(shù)的反演

基于斷層自動剖分的反演過程分為2部分： 一是基于均勻彈性半空間三角位錯模型獲得格林函數(shù)表達式，并通過正則化優(yōu)化算法求得滑動解； 另一部分是斷層自動剖分過程，通過高斯擬合計算獲得所有位錯元的最佳尺寸以及滑動在所有位錯元的平滑程度情況，從而對當(dāng)前的位錯元尺寸和平滑尺度進行調(diào)節(jié)，迭代直至滿足終止條件為止，產(chǎn)生的滑動分布如圖7和圖8所示。

圖7　斷層滑動平面分布

圖8　斷層滑動的地下空間分布

利用斷層自動剖分方法將引發(fā)巴姆地震的右旋走滑斷層剖分為151個大小不一的三角形位錯元，每個位錯元上對應(yīng)不同的滑動值?？傮w滑動值分布在0.2～2.35 m之間，集中分布于地下0.5～11 km范圍內(nèi)，最大滑動值為2.35 m，出現(xiàn)在地下4.8 km處，平均滑動值為1.8 m。尺寸較小的三角形位錯元分布在距離地表較近的地方，該區(qū)域滑動變化梯度較大，較小尺寸的位錯元能夠精細地反映滑動的變化情況； 而尺寸較大的三角形位錯元分布在離地表較遠的斷層深處，滑動變化比較平緩，這反映了三角形可變剖分模型的平滑特征。

3.3結(jié)果分析與評價

為了評價反演所得發(fā)震斷層參數(shù)的精度，需要通過正演解算獲得模擬形變圖，即由斷層參數(shù)和格林函數(shù)求得模擬形變數(shù)據(jù)，并計算模擬形變量與真實形變之間的殘差情況，以此為依據(jù)分析反演效果。模型殘差可由式(6)獲得，即

ε=G(m)-d，

(6)

式中ε表示殘差。將斷層參數(shù)作為已知數(shù)據(jù)，正演獲得巴姆震區(qū)的模擬形變?nèi)鐖D9所示。

圖9　LOS向模擬形變

模擬形變與原形變干涉分布特征相似，形變集中區(qū)域可分為4個象限，形變分布情況與真實形變較吻合，但東南側(cè)形變在25～30 cm之間的區(qū)域比真實形變范圍略大，分析是由模擬數(shù)據(jù)內(nèi)插至原始InSAR數(shù)據(jù)量級的過程中產(chǎn)生的誤差所致?？傮w模擬形變范圍為-16.45～34.16 cm之間，與真實形變范圍接近，殘差的均方根誤差值為0.97 cm。如圖10所示，由斷層自動剖分方法反演所得的模型殘差分布于-3.2～6.4 cm之間，整體殘差以0為中心呈正態(tài)分布。

圖10　反演結(jié)果的誤差分布直方圖

4結(jié)論

1)本文引入斷層自動剖分方法，以巴姆地震為例，對發(fā)震斷層的滑動參數(shù)進行反演，并通過正演比較模擬形變量與InSAR真實形變之間的殘差，表明利用斷層自動剖分方法對淺源單斷層進行反演可獲得較可靠滑動參數(shù)解。

2)利用斷層自動剖分方法可以兼顧殘差最小化和模型平滑性原則，避免了均勻剖分或人為設(shè)定剖分引起的“偽滑動”或滑動過于平滑的現(xiàn)象，提高解的可靠性。

3)本文仍存在一些不足之處。對于淺源單斷層引起的地震，斷層自動剖分方法可以反演獲得較可靠的滑動參數(shù)解，但對于震級較高且發(fā)震斷層構(gòu)造較復(fù)雜的地震，該方法是否適用還有待于進一步研究論證； 地震反演結(jié)果的精度評價一直是該領(lǐng)域長期存在的普遍性問題，僅通過正演模擬進行殘差比較具有一定的局限性，也有待于開展進一步研究。

參考文獻(References)：

[1]洪順英,申旭輝,單新建,等.基于升降軌ASAR的于田Ms 7.3級地震同震形變場信息提取與分析[J].國土資源遙感,2010,22(4):98-102.doi:10.6046/gtzyyg.2010.04.20.

Hong S Y,Shen X H,Shan X J,et al.The calculation and analysis of the co-seismic deformation field of Yutian Ms 7.3 earthquake basing on the ascending and descending orbit ASAR data[J].Remote Sensing for Land and Resource,2010,22(4):98-102.doi:10.6046/gtzyyg.2010.04.20.

[2]羅旭巍,孫建寶,沈正康,等.基于InSAR同震形變觀測反演2010年新西蘭南島Mw7.1 Darfield地震同震破裂分布[J].地球物理學(xué)報,2013,56(8):2613-2624.

Luo X W,Sun J B,Shen Z K,et al.Co-seismic slip distribution of 2010 Darfield, New Zealand Mw7.1 earthquake inverted using InSAR measurements[J].Chinese Journal of Geophysics,2013,56(8):2613-2624.

[3]薛蓮,孫建寶,沈正康.2010年1月12日海地Mw7.0級地震InSAR同震形變觀測及同震滑動分布反演[J].地震地質(zhì),2011,33(1):157-174.

Xue L,Sun J B,Shen Z K.InSAR coseismic deformation observation of the Jan 12th, 2010 HaiTi earthquake and its coseismic slip distribution inversion[J].Seismology and Geology,2011,33(1):157-174.

[4]周輝,馮光財,李志偉,等.利用InSAR資料反演緬甸Mw6.8地震斷層滑動分布[J].地球物理學(xué)報,2013,56(9):3011-3021.

Zhou H,Feng G C,Li Z W,et al.The fault slip distribution of the Myanmar Mw6.8 earthquake inferred from InSAR measurements[J].Chinese Journal of Geophysics,2013,56(9):3011-3021.

[5]溫揚茂,何平,許才軍,等.聯(lián)合Envisat和ALOS衛(wèi)星影像確定L’Aquila地震震源機制[J].地球物理學(xué)報,2012,55(1):53-65.

Wen Y M,He P,Xu C J,et al.Source parameters of the 2009 L’Aquila earthquake, Italy from Envisat and ALOS satellite SAR images[J].Chinese Journal of Geophysics,2012,55(1):53-65.

[6]Resor P G,Pollard D D,Wright T J,et al.Integrating high-precision aftershock locations and geodetic observations to model coseismic deformation associated with the 1995 Kozani-Grevena earthquake,Greece[J].Journal of Geophysical Research:Solid earth,2005,110(B9):B09402.

[7]Walters R J,Elliott J R,D’Agostino N,et al.The 2009 L’Aquila earthquake(central Italy):A source mechanism and implications for seismic hazard[J].Geophysical Research Letters,2009,36(17):L17312.

[8]張國宏,屈春燕,汪馳升,等.基于GPS和InSAR反演汶川Mw7.9地震斷層滑動分布[J].大地測量與地球動力學(xué),2010,30(4):19-24.

Zhang G H,Qu C Y,Wang C S,et al.Inversion of slip distribution of 2008 Wenchuan Mw7.9 earthquake constrained jointly by InSAR and GPS measurements[J].Journal of Geodesy and Geodynamics,2010,30(4):19-24.

[9]Barnhart W D,Lohman R B.Automated fault model discretization for inversions for coseismic slip distributions[J].Journal of Geophysical Research,2010,115(B10):B10419.

[10]Meade B J.Algorithms for the calculation of exact displacements,strains, and stresses for triangular dislocation elements in a uniform elastic half space[J].Computer and Geosciences,2007,33(8):1064-1075.

[11]Golub G H,Hansen P C,O’Leary D P.Tikhonov regularization and total least squares[J].SIAM Journal on Matrix Analysis and Applications,1999,21(1):185-194.

[12]Du Y J,Aydin A,Segall P.Comparison of various inversion techniques as applied to the determination of a geophysical deformation model for the 1983 Borah Peak earthquake[J].Bulletin of the Seismological Society of America,1992,82(4):1840-1866.

[13]Okada Y.Surface deformation due to shear and tensile faults in a half-space[J].Bulletin of the Seismological Society of America,1985,75(4):1135-1154.

[14]凌勇,曾祺明,羅揚,等.巴姆地震變形場和應(yīng)力場:I.用差分干涉雷達和Okada方法求解[J].巖石學(xué)報,2006,22(9):2367-2374.

Ling Y,Zeng Q M,Luo Y,et al.Deformation and stress field of Bam earthquake:Ⅰ.InSAR and Okada calculation[J].Acta Petrologica Sinica,2006,22(9):2367-2374.

(責(zé)任編輯： 陳理)

Inversion of fault slip parameters based on DInSAR and automated fault model discretization method

CHEN Danlei1, LIU Guoxiang2, WANG Xiaowen2, WANG Lei1, PU Huilong1

(1. The Third Engineering Institution of Surveying and Mapping of Sichuan, Chengdu 610500, China; 2. Faculty of GeosciencesandEnvironmentalEngineering,SourthwestJiaotongUniversity,Chengdu610031,China)

Abstract:The inversion of geometry and motion parameters of the seismogenic fault based on elastic dislocation model and constrained by coseismic deformation data obtained by DInSAR has great significance for understanding the activity characteristics of the fault and assessing the risk of potential disasters. The existing inversion methods of fault slip parameters are mostly based on the uniform discretization with rectangular dislocation units or artificially setting discretizaiton, which will lead to “pseudo-slip” or the phenomenon that the slip distribution is too smooth to reflect the slip details on the fault plane. Therefore, the automated fault model discetization method is introduced in this paper, which takes into account the principle of the minimized residuals and the optimal smoothing scales of the model, so the reliable solution can be obtained under the effective constraints. The inversion of geometry and slip parameters of fault in Bam earthquake is taken as an example. The experimental results show that using the automated fault model discretization method to invert the motion parameters of the single fault can generate reliable results.

Keywords:automated fault model discretization； DInSAR； inversion of slip parameters； mixed resampling

doi:10.6046/gtzyyg.2016.03.05

收稿日期：2015-02-02；

修訂日期：2015-03-20

基金項目：四川省測繪地理信息局科技計劃項目 “基于Web的四川省地理國情監(jiān)測數(shù)據(jù)成果展示方法與實現(xiàn)”(編號：J2014ZC16)資助。

中圖法分類號：TP79

文獻標志碼：A

文章編號：1001-070X(2016)03-0025-06

第一作者簡介：陳丹蕾(1987-)，女，碩士研究生，助理工程師，主要從事InSAR技術(shù)應(yīng)用與遙感應(yīng)用等方面的科研工作。Email: 281560206@qq.com。

引用格式： 陳丹蕾，劉國祥，王曉文,等.基于DInSAR技術(shù)與斷層自動剖分方法反演斷層滑動參數(shù)[J].國土資源遙感,2016,28(3):25-30.(ChenDL,LiuGX,WangXW,etal.InversionoffaultslipparametersbasedonDInSARandautomatedfaultmodeldiscretizationmethod[J].RemoteSensingforLandandResources,2016,28(3):25-30.)