基于非均勻地殼模型的格林函數法反演2001年昆侖山口西Ms8.1地震的同震滑移*

周 宇 聶兆生 賈治革 陳正松 楊少敏

(中國地震局地震研究所(地震大地測量實驗室)，武漢 430071)

基于非均勻地殼模型的格林函數法反演2001年昆侖山口西Ms8.1地震的同震滑移*

周 宇 聶兆生 賈治革 陳正松 楊少敏

(中國地震局地震研究所(地震大地測量實驗室)，武漢 430071)

建立基于地殼不均勻的三維有限元模型，將昆侖山口西Ms8.1地震產生的長約400 km的破裂帶在走向上分為16段，生成各個破裂段滑移相對于地表位移的格林函數，并以GPS觀測的水平位移為約束，用阻尼最小二乘法反演破裂滑移分布，結果表明反演結果與地質考察的地表破裂數據接近。用反演結果重構的地表位移場與觀測值一致，并體現了斷裂帶南北兩側不同的形變模式。

有限元;最小二乘法;昆侖山口西Ms8.1地震;同震滑移;GPS

1 引言

昆侖山口西Ms8.1地震發生在東昆侖斷裂帶上，該斷裂帶是青藏高原一條古老的縫合帶，其左旋走滑非常活躍。震后，美國哈佛大學、美國地質調查局(USGS)、日本東京大學地震研究所和中國地震局地球物理研究所先后給出了此次地震的震中位置及矩張量震源機制解，結果略有差異，但多認為此次主震的破裂方式以走滑破裂為主，兼有少量的逆沖滑動分量［1，2］。野外地質考察發現，地震沿東昆侖斷裂帶西段形成了長約350～420 km的地表破裂帶，位錯以走滑為主，最大左旋水平錯距 6～7 m［3-5］。

文獻［1，2］運用全球地震波資料反演了該地震的時空破裂過程。文獻［1］的反演結果為:斷層的破裂長度約380 km、寬約30 km、最大滑動量5.8 m、平均錯動(左旋走滑)1.8 m;文獻［2］假定斷層破裂深度為40 km，得到的斷層最大滑動量為2.2 m，平均滑動量僅1.2 m。顯然這些結果與地表破裂的實測滑動量［3-5］差別較大。

為使反演結果更加真實地表現地表破裂現象及斷層的活動狀況，文獻［6～10］分別以大地測量結果為約束，計算了該地震的同震位移，并討論了地表形變特征;單新建［11］則利用InSAR資料給出了衛星視線向的同震形變場。無論是用GPS還是InSAR，地表形變的觀測結果都表明，東昆侖斷裂帶南北兩側的形變顯示出明顯的非對稱性，斷層南側的同震位移和水平位移梯度的南北向分量都明顯大于斷層北側，表明斷層南北兩側介質的彈性力學性質有明顯的差異。這些反演結果與實地考察結果具有較好的一致性，地震波速資料也證實了反演結果的可信性［12-14］。

但這些研究使用的是Okada模型的解析公式，而Okada模型基于均勻彈性半空間，所以無法解釋南北盤形變特征的明顯差異。要解決這一問題，必須加入介質在橫向和垂向上的不均勻性，然而目前的位錯理論并沒有非均勻介質條件下的解析解，所以采用數值模擬方法就成為必然。王輝等［15］用三維有限元數值方法，考慮介質的不均勻性模擬了昆侖山口西地震的同震形變場，但該研究的目的在于重構空間連續的形變場，將野外地質勘查獲得的地表走滑量作為地震破裂分布。然而，并非所有的破裂都會到達地表，而地表的破裂不一定和深部的破裂一致。Masterlark［16］發現，在板塊俯沖帶同震及震后形變的數值模擬中，在非均勻地殼模型中用基于Okada模型反演得到的滑動分布來模擬地表形變場，可能會導致不可忽視的誤差。為解決這個問題，本文提出用三維有限元建立非均勻的介質模型，計算各個破裂段對于地表位移的格林函數，然后用GPS觀測數據反演各個破裂段的滑移。

2 模型的建立

昆侖山地震后，中國地震局地震研究所聯合中國地震局第二地形變監測中心迅速開展了應急觀測，獲得了40個點的同震位移數據。我們將利用這些數據，以單日時段為基本單位解算各GPS點的三維坐標。由于地震位錯以走滑為主，GPS站點的垂直位移不大，而GPS觀測在垂向上的誤差比水平分量的要高出一個數量級，所以本文拋棄垂向位移，只用水平位移。

2.1 單元劃分

考慮巖石圈縱向上的分層和斷裂帶南北兩側的分界，我們依據亞東格爾木額濟納旗地學斷面建立的巖石圈地質模型［9］，并參考此斷面得到的介質密度、速度、泊松比和莫霍面深度，用縱波速和泊松比換算出楊氏彈性模量。研究區介質的分區見圖1，參數的設置見表1，單元的劃分如圖2。

圖1 模型剖面示意圖Fig.1 Sketch of model profile

表1 介質參數Tab.1 Medium parameters

模型中我們將斷裂同一側的介質視為橫向均勻。

2.2 三維有限元網格與邊界條件

網格總體長1 500 km，寬1 000 km，深200 km，頂面的中心與斷裂帶中點重合。斷裂帶用長400 km、深14 km，完全與地面垂直的矩形面表示，走向277°(圖2)。模型的網格采用四面體單元，破裂面附近形變梯度大的地方做精細劃分，遠場位移梯度小的地方做粗疏劃分。數據試驗表明，在破裂面附近的網格單元的邊長應該比破裂面的邊長小一個數量級，才能保證有限元數值計算的最大誤差在2%以內，因此在斷裂帶附近的網格邊長設定為2～4 km，在邊界處邊長為35～50 km。網格共有節點58 349個，四面體單元313 502個(圖2)。

約束條件為:網格4個側面的邊界的水平位移為零，底面的垂向位移為零，側面和底面其他方向上的位移分量自由，地表位移的3個分量均為自由。

使用Pylith軟件進行有限元計算。

圖2 三維有限元網格的水平面Fig.2 Horizontal surface of 3D-FEM mesh

2.3 格林函數的生成

設任意地面點xi的位移矢量vi為關于斷層某個區域ξj滑移矢量uj的線性函數:

式中g(u)為格林函數。

將式(1)寫成矩陣形式，有:

由于觀測條件所限，GPS點不多，分布也不夠理想，因此沒有基于破裂模式在平面上對斷裂帶做不等長、垂向上做不等深的精細劃分，為方便起見，將破裂帶分為16個等長的分段，每段為長25 km、深14 km的矩形。并設定滑動向量中的一個分量為1，其他分量為零。用三維有限元計算正演各GPS站的位移。將2M個位移向量從左到右并列即可獲得模型的格林函數矩陣。

3 反演方法

以最小方差為反演目標函數，用廣義逆反演超定方程(2)，有:

其中CV是觀測向量V的協方差矩陣。

由于各GPS站觀測值獨立，且同一個GPS站的位移觀測值的兩個水平分量相關系數接近于零，故CV為對角陣，有:

σi是觀測向量第i個分量的標準差。為表示方便，令:

代入方程(3)有

模擬結果表明，用方程(8)反演，其結果發生了畸變，破裂面滑移的絕對值遠遠超出了正常范圍。經分析，發生畸變是因為矩陣G'TG'的最大和最小特征值相差太大，比率高達1011，導致反演的病態，所以無法用解決超定問題的方法來反演。由此本文改用阻尼最小二乘法進行反演，由于阻尼最小二乘法加入阻尼因子避免了病態問題。反演方程如下:

其中:β＞0，為阻尼因子;I為單位矩陣。殘差大小為:

其中norm()為向量的2-范數，tr()為矩陣的跡。

反演結果顯示，雖然?U相差很大，但都在GPS點上產生了幾乎一樣的位移。分析其原因，一是觀測數據不足，導致數據向量對模型向量不敏感;二是阻尼因子太小無法解決反演的病態問題。由此我們引入地震矩來選取合適的阻尼因子。地震矩公式如下:

其中Aj是第j個破裂面的面積，sj是第j個破裂面的滑動距離，μ是剪切模量。哈佛矩心張量的地震矩為5.9×1020Nm，與之對應的阻尼因子β=10-3.5。

最終，反演獲得的破裂面滑動分布的最優解如圖3所示。

圖3 基于非均勻地殼模型的同震滑移分布的反演結果Fig.3 Inversed coseismic slip distribution based on the heterogeneous crust model

4 反演結果分析

從圖3可以看到，反演得到的同震滑移以走滑為主，傾滑分量在大多數分段上小于50 cm。左旋走滑分布有三個峰值，從大到小分別分布于東經93.6°、92.6°與94.4°處，其量級分別為5.45、4.84和3.82 m，與地質考察的結果［5］基本一致。但傾滑分量在兩個分段上超過了1 m，與實際不符。

將反演獲得的滑動分布，代入公式(2)計算得到各GPS點的水平位移矢量如圖4所示。

圖4 水平位移的模擬值與GPS觀測值的比較Fig.4 Comparison between the simulated horizontal displacement field and the GPS observations

圖4中，紅色箭頭是GPS觀測得到的水平位移，綠色箭頭是模擬的水平位移，誤差橢圓置信區間為95%。其中BS33和WT02兩個點的位移超出圖幅，為表示方便，圖中比例尺所示的位移大小為真實值的一半。在大部分觀測點上，觀測值和模擬值都比較一致，在幾個離斷層較近的點上二者幾乎相同。模擬位移也很好地體現了斷裂帶南側的水平位移的南北向梯度比北側的大。

BS26、G0CQ兩個點的震前觀測時間分別在1995年和1997年［7］，而1995—2001年該地區發生5級以上地震10次，所以GPS觀測位移可能包含了這些地震引起的位移，使得觀測值遠大于模擬值。

在斷裂南側，除了近場的WT02，其他幾個點的模擬位移比實測的位移小，但在各點上二者相差的比例幾乎相同，這說明模擬結果的位移梯度或者應變場與觀測結果幾乎是相同的。引起位移不一致的原因可能是在WT02和WUDA之間存在著彈性模量較大的局部區域，導致該區域的應變較小，使得從WT02到WUDA的位移梯度變小。

處于阿爾金斷裂帶附近的幾個點上(圖4中的西北角)，模擬位移和觀測位移的大小比較接近，但是觀測位移為北東向，模擬位移為北西向，二者近于垂直。經分析可能是沒有用軟弱帶模擬阿爾金斷裂所致。

從整體上看，斷裂帶西側的模擬值與觀測值的差別明顯大于東側，從反演的角度來看，可能是斷裂帶西半段的近場沒有觀測點約束之故。

圖5 基于彈性半空間模型的同震滑移分布的反演結果Fig.5 Inversed coseismic slip distribution based on the homogeneous elastic half space model

圖6 非均勻地殼模型對于兩種不同模型的滑動反演結果的地表響應Fig.6 Ground surface response of heterogeneous crust model to two kinds of inversed slip distribution

GPS觀測不能提供同震位移場在空間上的連續分布，在反演獲得滑移分布之后帶入有限元模型正演，可以重建全區域的位移場和應變場。為了對比，本文建立了均勻彈性半空間的模型，用Okada公式計算格林函數，并用同樣的反演方法計算了同震滑移分布。雖然均勻彈性半空間模型也可以在觀測點上產生與非均勻地殼模型幾乎相同的模擬位移，然而其反演的滑移分布(圖5)，與非均勻地殼模型有明顯差別，雖然起伏趨勢基本一致，但其量值有差別，如走滑分量最大值為7.2 m，最小值為-1.5 m，傾滑分量最大值為3.6 m，并出現了右旋走滑，且傾滑分量過大，顯然不如非均勻模型更接近真實情形。將均勻彈性半空間模型的滑移反演結果代入非均勻模型中，兩種模型的位移場比較如圖6所示，在絕大多數點上，均勻彈性半空間模型的位移比非均勻模型的位移小30%以上，這說明基于地殼非均勻性反演同震滑移更加合理，而在非均勻性模型中直接使用均勻模型反演得到滑移分布是值得商榷的。

6 結論

1)阻尼最小二乘法可以克服反演中的病態問題，而不必事先限定滑移分布的上下界。地質考察結果驗證了反演結果的可靠性。雖然沒有水準測量數據，使得斷裂傾滑分量的反演結果不夠準確，然而在觀測數據如此少的情況下能得到這樣的結果，充分證明了反演方法的有效性。

2)基于地殼橫向和縱向不均勻的三維有限元模型模擬的同震位移場與觀測值基本一致，雖然由于模型沒有考慮地殼東西向的不均勻性和其他局部特征，導致在少數觀測點上與觀測值差別略大，但還是很好地解釋了東昆侖斷裂帶南北兩側形變模式的不同。

3)在沒有地表滑移分布的情況下，如果地殼存在明顯的橫向不均勻性，將均勻彈性半空間模型的滑移反演結果作為內邊界條件代入非均勻地殼模型來模擬形變場是不恰當的。

1 Lin A，Kikuchi M and Fu B.Rupture segmentation and process of the 2001 Mw7.8 central Kunlun earthquake，China［J］.Bull Seism Soc Amer.，2003，93:2 477-2 492.

2 許力生，陳運泰.用長周期波形資料反演2001年11月14日昆侖山口地震時空破裂過程［J］.中國科學(D輯)，2004，34(3):256-264.(Xu Lisheng and Chen Yuntai.Temporal and spatial rupture process of the great Kunlun Mountain Pass earthquake of Nov.14，2001，from GDSN long period waveform data［J］.Science in China(Ser D)，2004，34(3):256-264)

3 徐錫偉，等.2001年11月14日昆侖山庫湖地震(Ms8.1)地表破裂帶的基本特征［J］.地震地質，2002，24(1):1-13.(Xu Xiwei，et al.Characteristic featur es of the surface ruptures of the HohSai Hu(Kunlunshan)earthquake(Ms 8.1)northern Tibetan Plateau，China［J］.Seismology and Geology，2002，24(1):1-13)

4 Lin Aiming，et al.Coseismic strikeslip and rupture length produced by the 2001 Ms8.1 central Kunlun earthquake［J］.Science，2002，296:2 015-2 017.

5 陳杰，等.2001年昆侖山口西Ms8.1地震地表同震位移分布特征［J］.地震地質，2004，26(3):378-392.(Chen Jie，et al.Surface rupture zones of the 2001 earthquake Ms 8.1 west of Kunlun Pass，northern Qinghai-Tibet plateau［J］.Quaternary Sciences，2003，23(6):629-639)

6 喬學軍，等.昆侖山口西Ms8.1地震的地殼變形特征［J］.大地測量與地球動力學，2002，(4):6-11.(Qiao Xuejun，et al.Characteristics of crustal deformation relating to Ms8.1 Kunlunshan earthquake［J］.Journal of Geodesy and Geodynamics，2002，(4):6-11)

7 任金衛，王敏.GPS觀測的2001年昆侖山口西Ms8.1地震地殼變形［J］.第四紀研究，2005，25(1):34-44.(Ren Jinwei and Wang Min.GPS measured crustal deformation of the Ms8.1 Kunlun earthquake on November 14th 2001 in Qinhai-Tibet Plateau［J］.Quaternary Sciences，2005，25 (1):34-44)

8 譚凱，王琪，申重陽.用大地測量數據反演2001年昆侖山地震［J］.大地測量與地球動力學，2004，(3):47-50.(Tan Kai，Wang Qi and Shen Chongyang.Using geodetic data to inverse coseismic dislocation of 2001 Kunlun earthquake［J］.Journal of Geodesy and Geodynamics，2004，(3): 47-50)

9 萬永革，等.利用GPS和水準測量資料反演2001年昆侖山口西8.1級地震的同震滑動分布［J］.地震地質，2004，26(3):393-404.(Wan Yongge，et al.Co-seismic slip distribution of the 2001 west of kunlun mountain pass earthquake inverted by GPS and levelling data［J］.Seismology and Geology，2004，26(3):393-404)

10 馬超，單新建.昆侖山Ms8.1地震震源參數的多破裂段模擬研究［J］.地球物理學報，2006，49(2):428-437.(Ma Chao and Shan Xinjian.A multi-segment analytic modeling of hypocentral geometric characteristic parameters of the Ms8.1 earthquake at the Kunlun mountains［J］.Chinese J Geophys.，2006，49(2):428-437)

11 單新建，柳稼航，馬超.2001年昆侖山口西8.1級地震同震形變場特征的初步分析［J］.地震學報，2004，26 (5):474-480.(Shan Xinjian，Liu Jiahang and Ma Chao.Preliminary analysis on characteristics of coseismic deformation associated with Ms8.1 western Kunlunshan Pass earthquake in 2001［J］.Acta Seimologica Sinica，2004，26 (5):474-480)

12 吳功建，高銳，欽范.青藏高原“亞東-格爾木地學斷面”綜合地球物理調查與研究［J］.地球物理學報，1991，34 (5):552-561.(Wu Gongjian，Gao Rui and Yu Qinfan.Integrated inverstigation of the Qinghai-Tibet plateau along the Yadong-Golmud geosciences transect［J］.Chinese J Geophys.，1991，34(5):552-562)

13 Zhu L and Helmberger D V.Moho offset across the northern margin of the Tibetan plateau［J］.Science，1998，281: 1 170-1 172.

14 Okada Y.Surface deformation due to shear and tensile faults in a half space［J］.Bull seism Soc Am.，1985，75: 1 135-1 154.

15 王輝，等.昆侖山口西Ms8.1地震同震影響場的數值模擬［J］.地震地質，2007，29(3):637-647.(Wang Hui，et al.The numerical simulation on coseismic effect of the 14 November 2001 great Kunlun earthquake，Northern Tibet，China［J］.Seismology and Geology，2007，29(3):637-647)

16 Masterlark T，Demets C and Wang H F.Homogeneousvs heterogeneous subduction zone models:Coseismic and postseismic deformation［J］.J Geophys Res.，2001，28(21): 4 047-4 050.

INVERSION OF COSEISMIC SLIP DISTRIBUTION OF 2001 Ms8.1 KUNLUN EARTHQUAKE FROM GPS OBSERVATIONS WITH GREEN’S FUNCTION METHOD BASED ON A HETEROGENEOUS CRUST MODEL

Zhou Yu，Nie Zhaosheng，Jia Zhige，Chen Zhengsong and Yang Shaomin
(Key Laboratory of Earthquake Geodesy，Institute of Seismology，CEA，Wuhan 430071)

2001 Kunlun Ms8.1 earthquake ruptured as long as about 400 km.The displacement gradient perpendicular to the strike of rupture belt is significantly greater on the south side than that on the north.This fact cannot be explained by Okada model which is based on elastically homogeneous half space.A horizontally and vertically heterogeneous 3D-FEM model is established to generate Green’s functions of slip on 16 segments of the rupture belt with respect to the ground displacement field.With the damping LS(least square)algorithm，the slip distribution of rupture is inversed from observed horizontal displacements of 40 GPS sites，and the result is in good agreement with that from geological field investigation.The ground displacement field rebuilt by the inversed result is close to that observed by GPS，demonstrating different deformation patterns on the two sides of the rupture.

finite element;least square algorithm;Kunlun Ms8.1 earthquake;coseismic slip;GPS

1671-5942(2012)03-0022-05

2012-02-24

中國地震局地震行業科研專項(201208009)

周宇，男，1983年生，碩士，研究實習員，主要從事地震形變與活動構造的數值模擬與動力學研究.E-mail:zhouyuking@foxmail.com

P315.72+5

A