畢研磊++張永志++曹海坤++楊九元+楊珍

摘要：基于GPS數(shù)據(jù)，利用斷層自動剖分技術(shù)來反演汶川MW79同震滑動分布。反演結(jié)果顯示：（1）斷層自動剖分的滑動分布反演結(jié)果在各個位錯元之間表現(xiàn)出良好的平滑性。（2）此次地震是一次逆沖型地震，并帶有少量的右旋走滑分量。其中映秀―北川斷裂滑動最大值為15 m，平均滑動量261 m；灌縣―江油斷裂滑動的最大值為8 m，平均滑動量為158 m。反演得到此次地震的地震矩M0為9153 8×1020 N·m，矩震級MW為791。（3）正演結(jié)果與實(shí)測數(shù)據(jù)對比顯示，兩者位移場在運(yùn)動趨勢上基本保持一致，擬合殘差的均方誤差在合理范圍內(nèi)，驗(yàn)證了斷層自動剖分技術(shù)的可行性。

關(guān)鍵詞：斷層自動剖分；同震滑動分布；位錯元；地震矩；GPS數(shù)據(jù)

中圖分類號：P315241文獻(xiàn)標(biāo)識碼：A文章編號：1000-0666（2017）02-0211-05

0引言

[KG（0.25mm]斷層的滑動分布反演是基于一定的模型得到斷層面上細(xì)部的運(yùn)動特征，利用大地測量觀測數(shù)據(jù)（GPS、InSAR以及重力數(shù)據(jù)等）可分析斷層的運(yùn)動規(guī)律。利用反演研究斷層的運(yùn)動特征對認(rèn)識大陸內(nèi)部變形以及地震孕育過程有重要意義（季靈運(yùn)等，2015）。以往在反演斷層滑動分布時，是通過將斷層剖分成若干個大小相等的矩形位錯元，基于Okada矩形位錯理論反演得到每個位錯元的滑動量，最終得到整條斷層的精細(xì)滑動分布（Okada，1992）。但是這種斷層剖分經(jīng)常會帶來位錯元平滑尺度以及數(shù)據(jù)解析能力的問題。本文采用一種全新的斷層剖分方法——自動剖分技術(shù)（Barnhart，Lohman，2010），該方法依據(jù)斷層幾何參數(shù)將整條斷層剖分成大小可變的三角位錯元，利用三角位錯模型求解斷層滑動分布，以期能更為準(zhǔn)確地反映數(shù)據(jù)的解析能力，減少數(shù)據(jù)不能約束的偽滑動，從而得到一個較為合理的斷層模型。最后以龍門山斷裂帶上汶川MW79地震為例，反演了該地震2條破裂帶的滑動分布。[KG）]

[JP2]1汶川地震破裂帶及GPS同震位移場[JP]

2008年發(fā)生的汶川MW79地震引起2條地表破裂帶，分別產(chǎn)生于映秀―北川斷裂帶和灌縣―江油斷裂帶上，其中映秀―北川斷裂為主發(fā)斷裂（王衛(wèi)民等，2008；張培震等，2008）。地震發(fā)生后，中國地震局以及美國地質(zhì)調(diào)查局（United States Geological Survey，簡稱USGS）第一時間對該地震進(jìn)行了反演分析，顯示此次地震為逆沖型地震，同時帶有少量的右旋走滑運(yùn)動，并給出了此次地震的地震矩以及矩震級。之后許才軍等（2009）和張國宏等（2010）對該地震作了進(jìn)一步的反演研究，得到了該斷裂帶的精細(xì)滑動分布，并對此次地震的震源機(jī)制作了細(xì)致分析。

本文利用斷層周圍的213個GPS觀測站點(diǎn)數(shù)據(jù)進(jìn)行研究，其中包含了同震近場及遠(yuǎn)場位移數(shù)據(jù)，點(diǎn)位的密度比較大，能夠充分展示斷層周圍地表形變的變化特征（國家重大科學(xué)工程“中國地殼運(yùn)動觀測網(wǎng)絡(luò)”項(xiàng)目組，2008）。汶川地區(qū)龍門山斷裂帶及GPS同震位移場分布，如圖1所示。從圖中可見，GPS位移場以映秀―北川斷裂為中心相向運(yùn)動，整體表現(xiàn)為龍門山區(qū)域向四川盆地的逆沖推覆運(yùn)動。同時，沿?cái)鄬觾蓚?cè)有較為稠密的GPS點(diǎn)，這為反演該斷層的精細(xì)滑動分布提供了較好的數(shù)據(jù)基礎(chǔ)。

2斷層自動剖分的原理與方法

首先，我們建立斷層中每個位錯元（三角位錯）滑動與地表點(diǎn)位移之間的格林函數(shù)：[KH*1]

[WTHX]Gm=d[WT][JY]（1）[KH*1D]

其中：[WTHX]G表示建立斷層運(yùn)動與地表位移變化的格林函數(shù)；m表示三角位錯元的滑動量；d[WT]表示含有噪聲的地表觀測數(shù)據(jù)，常用的地表位移變化數(shù)據(jù)有GPS、InSAR等觀測數(shù)據(jù)。

目前有很多方法可求解滑動量[WTHX]m[WT]，正則化算法是在反演問題中很常用的一種方法。本文利用Tikhonow正則化算法來求解線性反演問題（Funning et al，2005）。該算法引入2個附加信息：一是各個位錯元滑動量的光滑性約束信息，二是解空間的邊界約束條件。這2個附加信息不僅考慮了模擬殘差最小的情況，同時也保證了各個位錯元之間的平滑性問題（陳丹蕾，2014）。[KH*1]

min=[JB（=][KG-*3][JB（[][WTHX][KG*3]G[WTBX]λ[WTHX]L[JB）]]m-[JB（[]d0[JB）]][JB）=][JY]（2）[KH*1D]

其中：L[WT]是Laplace平滑矩陣，用以克服臨近斷層之間滑動的劇烈梯度變化；λ是平滑因子（即正則化的權(quán)重參數(shù)），其確定有很多種方法，本文采用jRi的方法來確定λ的大?。?fù)P茂等，2012b；陳丹蕾，2014），該方法既能減小數(shù)據(jù)噪聲的影響，又可防止模型過度擬合的現(xiàn)象，其公式如下：[KH*1]

jRi=[SX（]（[WTHX]d-Gm）2+∑[DD（]2n+3i=1[DD）]GC[WTBX]GPS[WTHX]G[WTBX]Tn[SX）][JY]（3）[KH*1D]

其中：[WTHX]C[WTBX]GPS為觀測值的方差-協(xié)方差矩陣。

由上述jRi獲取的方法可知，當(dāng)jRi值越小時，模擬的殘差值以及數(shù)據(jù)的噪聲越小，此時的斷層剖分也越合理。因此取jRi的最小值作為平滑因子λ，并結(jié)合公式（1）、（2）來求解斷層的滑動參數(shù)[WTHX]m[WT]。

[BW（S][BG（；N][BHDWG1*2，WK15mmZQ，WK140mm，WK15mmYQW][HT5"][CM（22mm]地震研究[CM）]40卷[BG）F][BW）]

[BW（D][BG（；N][BHDWG1*2，WKZQ0W][HT5"]第2期[JZ]

畢研磊等：利用斷層自動剖分技術(shù)反演汶川MW79地震滑動分布

[BG）F][BW）]

當(dāng)平滑因子確定以后，廣義逆矩陣[WTHX]G[WTBX]-g求取公式為：[KH*1]

[WTHX]G[WTBX]-g=[JB<3（][KG-*3][JB（[][KG*3]G[WTBX]λ[WTHX]L[WTBX][JB）]]T[KG-*3]*[KG-*4][JB（[][KG*3]G[WTBX]λ[WTHX]L[JB）]][JB>3）]-1G[WTBX]T[JY]（4）[KH*1D]

進(jìn)而可得出模型分辨率矩陣[WTHX]R為：[KH*1]

[WTHX]R=G[WTBX]-g[WTHX]G[JY]（5）[KH*1]

[JP2]對于同震滑動來說，模型分辨率矩陣[WTHX]R[WT]中的元素能夠反映位錯元的尺寸大小。當(dāng)位錯模型較為合理時，模型分辨率矩陣的行元素呈現(xiàn)出中心位于相應(yīng)斷層片的狄拉克函數(shù)（溫?fù)P茂等，2012a；陳丹蕾，2014）。當(dāng)位錯元模型不接近真實(shí)情況時，各個位錯元的滑動將會產(chǎn)生劇烈的梯度變化，某一位錯元的滑動就會覆蓋到相鄰或者更遠(yuǎn)的位錯元上。利用模型分辨率矩陣[WTHX]R[WT]，我們能確定哪些位錯是超定的，哪些是欠定的，并求出一個新的斷層滑動分布，以更好地反映數(shù)據(jù)的約束能力。[JP]

最后的反演目標(biāo)既要使每一個位錯元都受到數(shù)據(jù)的有力約束，也要使斷層的平滑性達(dá)到最好。本文采用Gauss函數(shù)來擬合模型分辨率矩陣與位錯元尺寸之間的關(guān)系，以求得位錯的最佳尺寸（Barnhart，Lohman，2010）。然后通過比較當(dāng)前位錯尺寸與最佳尺寸，選擇是否需要調(diào)節(jié)，最后根據(jù)新產(chǎn)生的位錯數(shù)量與上次位錯數(shù)量的差異率，來判別是否終止迭代。[HJ1.8mm]

[BT（12-*1]3汶川地震滑動分布反演

31斷層參數(shù)反演[BT）]

[JP4]本文采用Okada矩形位錯理論反演汶川地震2條斷層的參數(shù)，如斷層端點(diǎn)的位置、長度、寬度、走向、傾角、上頂埋深以及滑動角等。其計(jì)算公式如下：[KH*2][JP]

[WTHX]l+v=G0m0[WT][JY]（6）[KH*2D]

式中：[WTHX]l[WT]表示GPS觀測的同震位移數(shù)據(jù)；[WTHX]v表示觀測誤差；G0為發(fā)震斷層與地表位移變化關(guān)系的格林函數(shù)；m0[WT]為斷層參數(shù)。

由于此次地震造成斷裂帶破裂至地表，本文根據(jù)地質(zhì)考察及其他學(xué)者反演的結(jié)果（王敏，2009），假設(shè)斷層的上頂埋深為0 km，寬度為40 km，傾角為43°。利用斷層周邊213個GPS點(diǎn)的觀測數(shù)據(jù)，采用粒子群算法反演汶川地震2條斷裂帶的其余參數(shù)（張永志等，2006；王帥等，2014），結(jié)果列于表1。

32斷層滑動分布反演

汶川地震在地表形成的2條破裂帶中，映秀―北川破裂長約250 km，灌縣―江油破裂長約70 km（張國宏等，2010）。根據(jù)上述斷層幾何參數(shù)的反演結(jié)果，對映秀―北川斷裂帶沿走向方向取274 km，其中斷層NE端至青川縣附近，SW端至映秀鎮(zhèn)附近，包含了主斷層的大部分破裂帶。對于灌縣―江油斷裂帶，同樣沿走向方向選取115 km，幾乎包含了該斷裂的全部破裂部分。

本文采用GPS同震觀測數(shù)據(jù)對汶川地震2條破裂斷層進(jìn)行滑動分布反演。首先根據(jù)斷層幾何參數(shù)對2條斷層進(jìn)行初始剖分，并計(jì)算模型分辨率矩陣；然后利用Gauss函數(shù)擬合模型分辨率矩陣與三角位錯元尺寸的關(guān)系，以求得位錯元的最佳尺寸；通過比較當(dāng)前位錯尺寸與最佳尺寸，對三角元進(jìn)行優(yōu)化剖分，得到一組新的三角位錯元，并在此基礎(chǔ)上重新計(jì)算模型分辨率矩陣；重復(fù)上述過程，直至2次剖分的位錯元個數(shù)的差異率在允許范圍內(nèi)，則迭代停止，此時剖分的每一個位錯元都能滿足最佳尺寸的要求，本文所采用的容許差異率為01。經(jīng)過一系列的計(jì)算，最終將映秀―北川斷裂剖分成301個大小不等的三角位錯元，灌縣―江油段斷裂剖分成99個三角位錯元，反演得到的滑動分布如圖2所示，圖中箭頭表示位錯元區(qū)域滑動方向，從圖中可知：

（1）映秀―北川主破裂帶的滑動分布主要集中在3個區(qū)域，分別為青川縣、北川以及汶川映秀鎮(zhèn)附近，滑動區(qū)主要分布在地表以下0～20 km的區(qū)域內(nèi)。對于灌縣―江油斷裂，滑動幾乎遍布整[JP2]條斷裂帶，滑動區(qū)域主要分布在地表以下0～20 km范圍內(nèi)。[JP]

（2）汶川地震為一逆沖型地震，并帶有少量的右旋走滑運(yùn)動。在映秀―北川主斷裂帶，從SW端向NE端的3個主要滑動分布區(qū)，斷層的運(yùn)動形式也逐漸發(fā)生變化，在汶川映秀附近斷層的滑動最大值達(dá)到了12 m，運(yùn)動形式以逆沖為主，并帶有少量的右旋走滑。沿著斷層走向到達(dá)北川附近，又出現(xiàn)一個滑動分布集中區(qū)，最大滑動量達(dá)15 m，右旋走滑分量也有所增加。當(dāng)?shù)竭_(dá)青川縣附近時，運(yùn)動形式主要變?yōu)橛倚呋?，并帶有少量的逆沖運(yùn)動，滑動量的最大值也達(dá)到了11 m，反演得到主斷裂帶的平均滑動量為261 m。從主斷裂的反演結(jié)果來看，斷層的最大滑動量略大于前人的反演結(jié)果。張國宏等（2010）反演得到主斷裂的最大滑動量為10 m左右，而許才軍等（2009）反演得到的最大滑動量為12 m左右，這可能是由于采用的位錯模型不同或者剖分形式不同造成的。對于灌縣―江油斷裂帶，運(yùn)動形式主要表現(xiàn)為逆沖運(yùn)動，滑動的最大值為8 m左右，其平均滑動量為158 m。

（3）利用反演得到的平均滑動量，并根據(jù)Kanamori計(jì)算震級的經(jīng)驗(yàn)公式，得到此次地震的地震矩M0為9153 8×1020 N·m，矩震級MW為791，這與USGS（2008）所提供的結(jié)果較為一致。

33正演模擬

為了檢驗(yàn)該模型的效果，進(jìn)行了正演計(jì)算，并與實(shí)測的GPS位移場進(jìn)行對比，得到的殘差如圖3所示。從圖中可以看出，該方法正演的結(jié)果與實(shí)測結(jié)果運(yùn)動趨勢基本保持一致，整體呈現(xiàn)出以主破裂帶為中心的相向運(yùn)動。模擬殘差值的均方誤差在東西、南北方向分別為43 cm、52 cm，遠(yuǎn)場相對于近場的模擬效果更好，斷裂帶的南段比北段模擬效果更優(yōu)。其一方面原因可能為本文采用單一斷層面進(jìn)行剖分反演斷層滑動分布，而龍門山斷裂帶地質(zhì)構(gòu)造比較復(fù)雜，主斷裂也是由許多復(fù)雜的小斷裂組成，把主斷裂考慮為一條單一的斷層與實(shí)際的破裂狀況有一定差異，這樣在斷層模型與實(shí)際模型較為一致的區(qū)域，正演效果會更優(yōu)，當(dāng)在斷層模型與實(shí)際模型有所差異的區(qū)域，正演效果會有所下降，正如本文中斷裂帶南段比北段模擬效果更優(yōu)，這也間接說明在斷裂帶南段的斷層模型更接近于實(shí)際的斷層模型。另一方面原因可能為斷層運(yùn)動所引起的地表形變不僅受到位錯的影響，還受到一些地殼區(qū)域動力不均勻的影響，離斷層越近受到區(qū)域動力的影響越大，所以斷層遠(yuǎn)場區(qū)域的模擬效果較近場區(qū)域會更好，其次也可能由于龍門山斷裂帶西南端與鮮水河斷裂、安寧河斷裂等復(fù)雜斷裂帶斜交等因素造成的。

4結(jié)論

本文基于GPS數(shù)據(jù)，利用斷層自動剖分技術(shù)反演得到汶川地震2條破裂帶的滑動分布，經(jīng)分析研究，可得出以下主要結(jié)論：

（1）此次地震以逆沖為主，并帶有少量右旋走滑的運(yùn)動；

（2）從剖分效果來看，映秀―北川斷裂帶共剖分了301個三角元，灌縣―江油斷裂剖分了99個三角元，斷裂帶均剖分的較為細(xì)致，能夠反映斷層中的細(xì)部滑動特征，且各個位錯元之間的滑動量表現(xiàn)出良好的平滑性；

（3）反演得到映秀―北川斷裂最大滑動量為15 m，平均滑動量為261 m，灌縣―江油斷裂中最大滑動量達(dá)8 m，平均滑動量為158 m；

[JP2]（4）計(jì)算得到此次地震的地震矩M0為9153 8[JP]×1020 N·m、矩震級MW為791，與USGS所提供的結(jié)果較為一致；[JP2]

（5）由殘差圖看出，正演的結(jié)果在運(yùn)動趨勢上與實(shí)測結(jié)果基本保持一致，各個點(diǎn)位的殘差較小，其中在遠(yuǎn)場區(qū)域正演的結(jié)果比近場區(qū)域會更好。[JP]

綜上所述，基于斷層自動剖分的三角位錯模型能較好地解決汶川地震所引起的同震滑動分布，具有一定的參考價值。

[HTK]本人在撰寫論文時得到張永志教授、姜永濤博士、曹海坤和楊九元的幫助，以及王帥博士對程序的支持，在此向他們表示衷心感謝。

參考文獻(xiàn)：

陳丹蕾2014基于雷達(dá)差分干涉測量與斷層自動剖分方法的地震震源參數(shù)反演[D].成都：西南交通大學(xué)

國家重大科學(xué)工程“中國地殼運(yùn)動觀測網(wǎng)絡(luò)”項(xiàng)目組2008GPS測定的2008年汶川MS80地震的同震位移場[J].中國科學(xué)：地球科學(xué)，38（10）：1195-1206

季靈運(yùn)，劉立煒，郝明2015利用InSAR技術(shù)研究滇西南鎮(zhèn)康—永德地區(qū)現(xiàn)今地殼形變特征[J].地震研究，38（1）：84-89

王敏2009基于GPS同震位移場約束反演2008年512汶川大地震破裂空間分布[J].地球物理學(xué)報(bào)，52（10）：2519-2526

王帥，張永志，姜永濤，等2014維多樣性的動態(tài)權(quán)重粒子群算法反演斷層滑動速率[J].地球物理學(xué)進(jìn)展，（4）：1766-1771

王衛(wèi)民，趙連鋒，李娟，等2008四川汶川80級地震震源過程[J].地球物理學(xué)報(bào)，51（5）：1403-1410

溫?fù)P茂，何平，許才軍，等2012a聯(lián)合Envisat和ALOS衛(wèi)星影像確定L′Aquila地震震源機(jī)制[J].地球物理學(xué)報(bào)，55（1）：53-65

溫?fù)P茂，許才軍，劉洋，等2012b利用斷層自動剖分技術(shù)的2008年青海大柴旦MW63級地震InSAR反演研究[J].武漢大學(xué)學(xué)報(bào)（信息科學(xué)版），37（4）：458-462

許才軍，劉洋，溫?fù)P茂2009利用GPS資料反演汶川MW79級地震滑動分布[J].測繪學(xué)報(bào)，38（3）：195-201，215

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

張培震，徐錫偉，聞學(xué)澤，等20082008年汶川80級地震發(fā)震斷裂的滑動速率、復(fù)發(fā)周期和構(gòu)造成因[J].地球物理學(xué)報(bào)，51（4）：1066-1073

張永志，王衛(wèi)東，魏玉明，等2006利用GPS資料反演祁連山斷層的三維滑動速率[J].大地測量與地球動力學(xué)，26（1）：31-35

BARNHART W D，LOHMAN R B2010Automated Fault Model Discretization for Inversions for Coseismic Slip Distributions[J].Journal of Geophysical Research：Solid Earth，115（B10），DOI：101029/2010JB007545

FUNNING G J，PARSON B，WRIGHT T J，et al2005Surface Displacements and Source Parameters of the 2003 Bam（Iran）Earthquake from Envisat Advanced Synthetic Aperture Radar Imagery[J].Journal of Geophysical Research：Solid Earth，110（B9），DOI：101029/2004JB003338

OKADA Y1992Internal deformation due to shear and tensile faults in a half-space[J].Bulletin of the Seismological Society of America，82（2）：1018-1040

USGS2008Magnitude 79-EASTERN SICHUAN，China[EB/OL].（2008）[2016-11-01].http：//earthquakeusgsgov/eqdenter/eqinthenews/2008/us2008ryan[ZK）][HJ][FL）]

[STHZ][WT4HZ][JZ]The Wenchuan MW79 Earthquake Slip Distribution Inversion Based[JZ]on Automated Fault Discretization Method

[WT5B1][STBZ][JZ]BI Yanlei，ZHANG Yongzhi，CAO Haikun，YANG Jiuyuan，YANG Zhen

[WT5"B1X][JZ]（School of Geology Engineering and Geomatics，Changan University，Xian 710064，Shaanxi，China）

Based on GPS data，we have inversed the Wenchuan MW79 coseismic slip distributions by utilizing automated fault model discretization method The results show：（1）The slip distributions based on automated fault model discretization appear well smooth（2）It was a thrust-type earthquake with a small amount of right-lateral slip component In Yingxiu-Beichuan fracture zone，the maximal sliding reached 15 m，average magnitude is 261 m In the Guanxian-Jiangyou fracture zone，the maximal sliding reached 8 m，average sliding magnitude is 158 m Meanwhile the inversed seismic moment M0 is 9153 8×1020N·m corresponding to moment magnitude MW is 791（3）The comparison between forward modeling results and the measured results shows that the displacement fields are consistent on the trend of movement，the root mean square error（RMSE）within a reasonable range，and this also verifies the feasibility of automated fault model discretization method

[WTHZ]Keywords：[WTB1]automated fault discretization；coseismic slip distribution；dislocation element；seismic moment；GPS data