基于非均匀地壳模型的格林函数法反演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