利用断层自动剖分技术反演汶川MW79地震滑动分布
毕研磊++张永志++曹海坤++杨九元+杨珍
摘要：基于GPS数据，利用断层自动剖分技术来反演汶川MW79同震滑动分布。反演结果显示：（1）断层自动剖分的滑动分布反演结果在各个位错元之间表现出良好的平滑性。（2）此次地震是一次逆冲型地震，并带有少量的右旋走滑分量。其中映秀―北川断裂滑动最大值为15 m，平均滑动量261 m；灌县―江油断裂滑动的最大值为8 m，平均滑动量为158 m。反演得到此次地震的地震矩M0为9153 8×1020 N·m，矩震级MW为791。（3）正演结果与实测数据对比显示，两者位移场在运动趋势上基本保持一致，拟合残差的均方误差在合理范围内，验证了断层自动剖分技术的可行性。
关键词：断层自动剖分；同震滑动分布；位错元；地震矩；GPS数据
中图分类号：P315241文献标识码：A文章编号：1000-0666（2017）02-0211-05
0引言
[KG（0.25mm]断层的滑动分布反演是基于一定的模型得到断层面上细部的运动特征，利用大地测量观测数据（GPS、InSAR以及重力数据等）可分析断层的运动规律。利用反演研究断层的运动特征对认识大陆内部变形以及地震孕育过程有重要意义（季灵运等，2015）。以往在反演断层滑动分布时，是通过将断层剖分成若干个大小相等的矩形位错元，基于Okada矩形位错理论反演得到每个位错元的滑动量，最终得到整条断层的精细滑动分布（Okada，1992）。但是这种断层剖分經常会带来位错元平滑尺度以及数据解析能力的问题。本文采用一种全新的断层剖分方法——自动剖分技术（Barnhart，Lohman，2010），该方法依据断层几何参数将整条断层剖分成大小可变的三角位错元，利用三角位错模型求解断层滑动分布，以期能更为准确地反映数据的解析能力，减少数据不能约束的伪滑动，从而得到一个较为合理的断层模型。最后以龙门山断裂带上汶川MW79地震为例，反演了该地震2条破裂带的滑动分布。[KG）]
[JP2]1汶川地震破裂带及GPS同震位移场[JP]
2008年发生的汶川MW79地震引起2条地表破裂带，分别产生于映秀―北川断裂带和灌县―江油断裂带上，其中映秀―北川断裂为主发断裂（王卫民等，2008；张培震等，2008）。地震发生后，中国地震局以及美国地质调查局（United States Geological Survey，简称USGS）第一时间对该地震进行了反演分析，显示此次地震为逆冲型地震，同时带有少量的右旋走滑运动，并给出了此次地震的地震矩以及矩震级。之后许才军等（2009）和张国宏等（2010）对该地震作了进一步的反演研究，得到了该断裂带的精细滑动分布，并对此次地震的震源机制作了细致分析。
本文利用断层周围的213个GPS观测站点数据进行研究，其中包含了同震近场及远场位移数据，点位的密度比较大，能够充分展示断层周围地表形变的变化特征（国家重大科学工程“中国地壳运动观测网络”项目组，2008）。汶川地区龙门山断裂带及GPS同震位移场分布，如图1所示。从图中可见，GPS位移场以映秀―北川断裂为中心相向运动，整体表现为龙门山区域向四川盆地的逆冲推覆运动。同时，沿断层两侧有较为稠密的GPS点，这为反演该断层的精细滑动分布提供了较好的数据基础。
2断层自动剖分的原理与方法
首先，我们建立断层中每个位错元（三角位错）滑动与地表点位移之间的格林函数：[KH*1]
[WTHX]Gm=d[WT][JY]（1）[KH*1D]
其中：[WTHX]G表示建立断层运动与地表位移变化的格林函数；m表示三角位错元的滑动量；d[WT]表示含有噪声的地表观测数据，常用的地表位移变化数据有GPS、InSAR等观测数据。
目前有很多方法可求解滑动量[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平滑矩阵，用以克服临近断层之间滑动的剧烈梯度变化；λ是平滑因子（即正则化的权重参数），其确定有很多种方法，本文采用jRi的方法来确定λ的大小（温扬茂等，2012b；陈丹蕾，2014），该方法既能减小数据噪声的影响，又可防止模型过度拟合的现象，其公式如下：[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为观测值的方差-协方差矩阵。
由上述jRi获取的方法可知，当jRi值越小时，模拟的残差值以及数据的噪声越小，此时的断层剖分也越合理。因此取jRi的最小值作为平滑因子λ，并结合公式（1）、（2）来求解断层的滑动参数[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]
毕研磊等：利用断层自动剖分技术反演汶川MW79地震滑动分布
[BG）F][BW）]
当平滑因子确定以后，广义逆矩阵[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]
进而可得出模型分辨率矩阵[WTHX]R为：[KH*1]
[WTHX]R=G[WTBX]-g[WTHX]G[JY]（5）[KH*1]
[JP2]对于同震滑动来说，模型分辨率矩阵[WTHX]R[WT]中的元素能够反映位错元的尺寸大小。当位错模型较为合理时，模型分辨率矩阵的行元素呈现出中心位于相应断层片的狄拉克函数（温扬茂等，2012a；陈丹蕾，2014）。当位错元模型不接近真实情况时，各个位错元的滑动将会产生剧烈的梯度变化，某一位错元的滑动就会覆盖到相邻或者更远的位错元上。利用模型分辨率矩阵[WTHX]R[WT]，我们能确定哪些位错是超定的，哪些是欠定的，并求出一个新的断层滑动分布，以更好地反映数据的约束能力。[JP]
最后的反演目标既要使每一个位错元都受到数据的有力约束，也要使断层的平滑性达到最好。本文采用Gauss函数来拟合模型分辨率矩阵与位错元尺寸之间的关系，以求得位错的最佳尺寸（Barnhart，Lohman，2010）。然后通过比较当前位错尺寸与最佳尺寸，选择是否需要调节，最后根据新产生的位错数量与上次位错数量的差异率，来判别是否终止迭代。[HJ1.8mm]
[BT（12-*1]3汶川地震滑动分布反演
31断层参数反演[BT）]
[JP4]本文采用Okada矩形位错理论反演汶川地震2条断层的参数，如断层端点的位置、长度、宽度、走向、倾角、上顶埋深以及滑动角等。其计算公式如下：[KH*2][JP]
[WTHX]l+v=G0m0[WT][JY]（6）[KH*2D]
式中：[WTHX]l[WT]表示GPS观测的同震位移数据；[WTHX]v表示观测误差；G0为发震断层与地表位移变化关系的格林函数；m0[WT]为断层参数。
由于此次地震造成断裂带破裂至地表，本文根据地质考察及其他学者反演的结果（王敏，2009），假设断层的上顶埋深为0 km，宽度为40 km，倾角为43°。利用断层周边213个GPS点的观测数据，采用粒子群算法反演汶川地震2条断裂带的其余参数（张永志等，2006；王帅等，2014），结果列于表1。
32断层滑动分布反演
汶川地震在地表形成的2条破裂带中，映秀―北川破裂长约250 km，灌县―江油破裂长约70 km（张国宏等，2010）。根据上述断层几何参数的反演结果，对映秀―北川断裂带沿走向方向取274 km，其中断层NE端至青川县附近，SW端至映秀镇附近，包含了主断层的大部分破裂带。对于灌县―江油断裂带，同样沿走向方向选取115 km，几乎包含了该断裂的全部破裂部分。
本文采用GPS同震观测数据对汶川地震2条破裂断层进行滑动分布反演。首先根据断层几何参数对2条断层进行初始剖分，并计算模型分辨率矩阵；然后利用Gauss函数拟合模型分辨率矩阵与三角位错元尺寸的关系，以求得位错元的最佳尺寸；通过比较当前位错尺寸与最佳尺寸，对三角元进行优化剖分，得到一组新的三角位错元，并在此基础上重新计算模型分辨率矩阵；重复上述过程，直至2次剖分的位错元个数的差异率在允许范围内，则迭代停止，此时剖分的每一个位错元都能满足最佳尺寸的要求，本文所采用的容许差异率为01。经过一系列的计算，最终将映秀―北川断裂剖分成301个大小不等的三角位错元，灌县―江油段断裂剖分成99个三角位错元，反演得到的滑动分布如图2所示，图中箭头表示位错元区域滑动方向，从图中可知：
（1）映秀―北川主破裂带的滑动分布主要集中在3个区域，分别为青川縣、北川以及汶川映秀镇附近，滑动区主要分布在地表以下0～20 km的区域内。对于灌县―江油断裂，滑动几乎遍布整[JP2]条断裂带，滑动区域主要分布在地表以下0～20 km范围内。[JP]
（2）汶川地震为一逆冲型地震，并带有少量的右旋走滑运动。在映秀―北川主断裂带，从SW端向NE端的3个主要滑动分布区，断层的运动形式也逐渐发生变化，在汶川映秀附近断层的滑动最大值达到了12 m，运动形式以逆冲为主，并带有少量的右旋走滑。沿着断层走向到达北川附近，又出现一个滑动分布集中区，最大滑动量达15 m，右旋走滑分量也有所增加。当到达青川县附近时，运动形式主要变为右旋走滑，并带有少量的逆冲运动，滑动量的最大值也达到了11 m，反演得到主断裂带的平均滑动量为261 m。从主断裂的反演结果来看，断层的最大滑动量略大于前人的反演结果。张国宏等（2010）反演得到主断裂的最大滑动量为10 m左右，而许才军等（2009）反演得到的最大滑动量为12 m左右，这可能是由于采用的位错模型不同或者剖分形式不同造成的。对于灌县―江油断裂带，运动形式主要表现为逆冲运动，滑动的最大值为8 m左右，其平均滑动量为158 m。
（3）利用反演得到的平均滑动量，并根据Kanamori计算震级的经验公式，得到此次地震的地震矩M0为9153 8×1020 N·m，矩震级MW为791，这与USGS（2008）所提供的结果较为一致。
33正演模拟
为了检验该模型的效果，进行了正演计算，并与实测的GPS位移场进行对比，得到的残差如图3所示。从图中可以看出，该方法正演的结果与实测结果运动趋势基本保持一致，整体呈现出以主破裂带为中心的相向运动。模拟残差值的均方误差在东西、南北方向分别为43 cm、52 cm，远场相对于近场的模拟效果更好，断裂带的南段比北段模拟效果更优。其一方面原因可能为本文采用单一断层面进行剖分反演断层滑动分布，而龙门山断裂带地质构造比较复杂，主断裂也是由许多复杂的小断裂组成，把主断裂考虑为一条单一的断层与实际的破裂状况有一定差异，这样在断层模型与实际模型较为一致的区域，正演效果会更优，当在断层模型与实际模型有所差异的区域，正演效果会有所下降，正如本文中断裂带南段比北段模拟效果更优，这也间接说明在断裂带南段的断层模型更接近于实际的断层模型。另一方面原因可能为断层运动所引起的地表形变不仅受到位错的影响，还受到一些地壳区域动力不均匀的影响，离断层越近受到区域动力的影响越大，所以断层远场区域的模拟效果较近场区域会更好，其次也可能由于龙门山断裂带西南端与鲜水河断裂、安宁河断裂等复杂断裂带斜交等因素造成的。
4结论
本文基于GPS数据，利用断层自动剖分技术反演得到汶川地震2条破裂带的滑动分布，经分析研究，可得出以下主要结论：
（1）此次地震以逆冲为主，并带有少量右旋走滑的运动；
（2）从剖分效果来看，映秀―北川断裂带共剖分了301个三角元，灌县―江油断裂剖分了99个三角元，断裂带均剖分的较为细致，能够反映断层中的细部滑动特征，且各个位错元之间的滑动量表现出良好的平滑性；
（3）反演得到映秀―北川断裂最大滑动量为15 m，平均滑动量为261 m，灌县―江油断裂中最大滑动量达8 m，平均滑动量为158 m；
[JP2]（4）计算得到此次地震的地震矩M0为9153 8[JP]×1020 N·m、矩震级MW为791，与USGS所提供的结果较为一致；[JP2]
（5）由残差图看出，正演的结果在运动趋势上与实测结果基本保持一致，各个点位的残差较小，其中在远场区域正演的结果比近场区域会更好。[JP]
综上所述，基于断层自动剖分的三角位错模型能较好地解决汶川地震所引起的同震滑动分布，具有一定的参考价值。
[HTK]本人在撰写论文时得到张永志教授、姜永涛博士、曹海坤和杨九元的帮助，以及王帅博士对程序的支持，在此向他们表示衷心感谢。
参考文献：
陈丹蕾2014基于雷达差分干涉测量与断层自动剖分方法的地震震源参数反演[D].成都：西南交通大学
国家重大科学工程“中国地壳运动观测网络”项目组2008GPS测定的2008年汶川MS80地震的同震位移场[J].中国科学：地球科学，38（10）：1195-1206
季灵运，刘立炜，郝明2015利用InSAR技术研究滇西南镇康—永德地区现今地壳形变特征[J].地震研究，38（1）：84-89
王敏2009基于GPS同震位移场约束反演2008年512汶川大地震破裂空间分布[J].地球物理学报，52（10）：2519-2526
王帅，张永志，姜永涛，等2014维多样性的动态权重粒子群算法反演断层滑动速率[J].地球物理学进展，（4）：1766-1771
王卫民，赵连锋，李娟，等2008四川汶川80级地震震源过程[J].地球物理学报，51（5）：1403-1410
温扬茂，何平，许才军，等2012a联合Envisat和ALOS卫星影像确定L′Aquila地震震源机制[J].地球物理学报，55（1）：53-65
温扬茂，许才军，刘洋，等2012b利用断层自动剖分技术的2008年青海大柴旦MW63级地震InSAR反演研究[J].武漢大学学报（信息科学版），37（4）：458-462
许才军，刘洋，温扬茂2009利用GPS资料反演汶川MW79级地震滑动分布[J].测绘学报，38（3）：195-201，215
张国宏，屈春燕，汪驰升，等2010基于GPS和InSAR反演汶川MW79地震断层滑动分布[J].大地测量与地球动力学，30（4）：19-24
张培震，徐锡伟，闻学泽，等20082008年汶川80级地震发震断裂的滑动速率、复发周期和构造成因[J].地球物理学报，51（4）：1066-1073
张永志，王卫东，魏玉明，等2006利用GPS资料反演祁连山断层的三维滑动速率[J].大地测量与地球动力学，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
摘要：基于GPS数据，利用断层自动剖分技术来反演汶川MW79同震滑动分布。反演结果显示：（1）断层自动剖分的滑动分布反演结果在各个位错元之间表现出良好的平滑性。（2）此次地震是一次逆冲型地震，并带有少量的右旋走滑分量。其中映秀―北川断裂滑动最大值为15 m，平均滑动量261 m；灌县―江油断裂滑动的最大值为8 m，平均滑动量为158 m。反演得到此次地震的地震矩M0为9153 8×1020 N·m，矩震级MW为791。（3）正演结果与实测数据对比显示，两者位移场在运动趋势上基本保持一致，拟合残差的均方误差在合理范围内，验证了断层自动剖分技术的可行性。
关键词：断层自动剖分；同震滑动分布；位错元；地震矩；GPS数据
中图分类号：P315241文献标识码：A文章编号：1000-0666（2017）02-0211-05
0引言
[KG（0.25mm]断层的滑动分布反演是基于一定的模型得到断层面上细部的运动特征，利用大地测量观测数据（GPS、InSAR以及重力数据等）可分析断层的运动规律。利用反演研究断层的运动特征对认识大陆内部变形以及地震孕育过程有重要意义（季灵运等，2015）。以往在反演断层滑动分布时，是通过将断层剖分成若干个大小相等的矩形位错元，基于Okada矩形位错理论反演得到每个位错元的滑动量，最终得到整条断层的精细滑动分布（Okada，1992）。但是这种断层剖分經常会带来位错元平滑尺度以及数据解析能力的问题。本文采用一种全新的断层剖分方法——自动剖分技术（Barnhart，Lohman，2010），该方法依据断层几何参数将整条断层剖分成大小可变的三角位错元，利用三角位错模型求解断层滑动分布，以期能更为准确地反映数据的解析能力，减少数据不能约束的伪滑动，从而得到一个较为合理的断层模型。最后以龙门山断裂带上汶川MW79地震为例，反演了该地震2条破裂带的滑动分布。[KG）]
[JP2]1汶川地震破裂带及GPS同震位移场[JP]
2008年发生的汶川MW79地震引起2条地表破裂带，分别产生于映秀―北川断裂带和灌县―江油断裂带上，其中映秀―北川断裂为主发断裂（王卫民等，2008；张培震等，2008）。地震发生后，中国地震局以及美国地质调查局（United States Geological Survey，简称USGS）第一时间对该地震进行了反演分析，显示此次地震为逆冲型地震，同时带有少量的右旋走滑运动，并给出了此次地震的地震矩以及矩震级。之后许才军等（2009）和张国宏等（2010）对该地震作了进一步的反演研究，得到了该断裂带的精细滑动分布，并对此次地震的震源机制作了细致分析。
本文利用断层周围的213个GPS观测站点数据进行研究，其中包含了同震近场及远场位移数据，点位的密度比较大，能够充分展示断层周围地表形变的变化特征（国家重大科学工程“中国地壳运动观测网络”项目组，2008）。汶川地区龙门山断裂带及GPS同震位移场分布，如图1所示。从图中可见，GPS位移场以映秀―北川断裂为中心相向运动，整体表现为龙门山区域向四川盆地的逆冲推覆运动。同时，沿断层两侧有较为稠密的GPS点，这为反演该断层的精细滑动分布提供了较好的数据基础。
2断层自动剖分的原理与方法
首先，我们建立断层中每个位错元（三角位错）滑动与地表点位移之间的格林函数：[KH*1]
[WTHX]Gm=d[WT][JY]（1）[KH*1D]
其中：[WTHX]G表示建立断层运动与地表位移变化的格林函数；m表示三角位错元的滑动量；d[WT]表示含有噪声的地表观测数据，常用的地表位移变化数据有GPS、InSAR等观测数据。
目前有很多方法可求解滑动量[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平滑矩阵，用以克服临近断层之间滑动的剧烈梯度变化；λ是平滑因子（即正则化的权重参数），其确定有很多种方法，本文采用jRi的方法来确定λ的大小（温扬茂等，2012b；陈丹蕾，2014），该方法既能减小数据噪声的影响，又可防止模型过度拟合的现象，其公式如下：[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为观测值的方差-协方差矩阵。
由上述jRi获取的方法可知，当jRi值越小时，模拟的残差值以及数据的噪声越小，此时的断层剖分也越合理。因此取jRi的最小值作为平滑因子λ，并结合公式（1）、（2）来求解断层的滑动参数[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]
毕研磊等：利用断层自动剖分技术反演汶川MW79地震滑动分布
[BG）F][BW）]
当平滑因子确定以后，广义逆矩阵[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]
进而可得出模型分辨率矩阵[WTHX]R为：[KH*1]
[WTHX]R=G[WTBX]-g[WTHX]G[JY]（5）[KH*1]
[JP2]对于同震滑动来说，模型分辨率矩阵[WTHX]R[WT]中的元素能够反映位错元的尺寸大小。当位错模型较为合理时，模型分辨率矩阵的行元素呈现出中心位于相应断层片的狄拉克函数（温扬茂等，2012a；陈丹蕾，2014）。当位错元模型不接近真实情况时，各个位错元的滑动将会产生剧烈的梯度变化，某一位错元的滑动就会覆盖到相邻或者更远的位错元上。利用模型分辨率矩阵[WTHX]R[WT]，我们能确定哪些位错是超定的，哪些是欠定的，并求出一个新的断层滑动分布，以更好地反映数据的约束能力。[JP]
最后的反演目标既要使每一个位错元都受到数据的有力约束，也要使断层的平滑性达到最好。本文采用Gauss函数来拟合模型分辨率矩阵与位错元尺寸之间的关系，以求得位错的最佳尺寸（Barnhart，Lohman，2010）。然后通过比较当前位错尺寸与最佳尺寸，选择是否需要调节，最后根据新产生的位错数量与上次位错数量的差异率，来判别是否终止迭代。[HJ1.8mm]
[BT（12-*1]3汶川地震滑动分布反演
31断层参数反演[BT）]
[JP4]本文采用Okada矩形位错理论反演汶川地震2条断层的参数，如断层端点的位置、长度、宽度、走向、倾角、上顶埋深以及滑动角等。其计算公式如下：[KH*2][JP]
[WTHX]l+v=G0m0[WT][JY]（6）[KH*2D]
式中：[WTHX]l[WT]表示GPS观测的同震位移数据；[WTHX]v表示观测误差；G0为发震断层与地表位移变化关系的格林函数；m0[WT]为断层参数。
由于此次地震造成断裂带破裂至地表，本文根据地质考察及其他学者反演的结果（王敏，2009），假设断层的上顶埋深为0 km，宽度为40 km，倾角为43°。利用断层周边213个GPS点的观测数据，采用粒子群算法反演汶川地震2条断裂带的其余参数（张永志等，2006；王帅等，2014），结果列于表1。
32断层滑动分布反演
汶川地震在地表形成的2条破裂带中，映秀―北川破裂长约250 km，灌县―江油破裂长约70 km（张国宏等，2010）。根据上述断层几何参数的反演结果，对映秀―北川断裂带沿走向方向取274 km，其中断层NE端至青川县附近，SW端至映秀镇附近，包含了主断层的大部分破裂带。对于灌县―江油断裂带，同样沿走向方向选取115 km，几乎包含了该断裂的全部破裂部分。
本文采用GPS同震观测数据对汶川地震2条破裂断层进行滑动分布反演。首先根据断层几何参数对2条断层进行初始剖分，并计算模型分辨率矩阵；然后利用Gauss函数拟合模型分辨率矩阵与三角位错元尺寸的关系，以求得位错元的最佳尺寸；通过比较当前位错尺寸与最佳尺寸，对三角元进行优化剖分，得到一组新的三角位错元，并在此基础上重新计算模型分辨率矩阵；重复上述过程，直至2次剖分的位错元个数的差异率在允许范围内，则迭代停止，此时剖分的每一个位错元都能满足最佳尺寸的要求，本文所采用的容许差异率为01。经过一系列的计算，最终将映秀―北川断裂剖分成301个大小不等的三角位错元，灌县―江油段断裂剖分成99个三角位错元，反演得到的滑动分布如图2所示，图中箭头表示位错元区域滑动方向，从图中可知：
（1）映秀―北川主破裂带的滑动分布主要集中在3个区域，分别为青川縣、北川以及汶川映秀镇附近，滑动区主要分布在地表以下0～20 km的区域内。对于灌县―江油断裂，滑动几乎遍布整[JP2]条断裂带，滑动区域主要分布在地表以下0～20 km范围内。[JP]
（2）汶川地震为一逆冲型地震，并带有少量的右旋走滑运动。在映秀―北川主断裂带，从SW端向NE端的3个主要滑动分布区，断层的运动形式也逐渐发生变化，在汶川映秀附近断层的滑动最大值达到了12 m，运动形式以逆冲为主，并带有少量的右旋走滑。沿着断层走向到达北川附近，又出现一个滑动分布集中区，最大滑动量达15 m，右旋走滑分量也有所增加。当到达青川县附近时，运动形式主要变为右旋走滑，并带有少量的逆冲运动，滑动量的最大值也达到了11 m，反演得到主断裂带的平均滑动量为261 m。从主断裂的反演结果来看，断层的最大滑动量略大于前人的反演结果。张国宏等（2010）反演得到主断裂的最大滑动量为10 m左右，而许才军等（2009）反演得到的最大滑动量为12 m左右，这可能是由于采用的位错模型不同或者剖分形式不同造成的。对于灌县―江油断裂带，运动形式主要表现为逆冲运动，滑动的最大值为8 m左右，其平均滑动量为158 m。
（3）利用反演得到的平均滑动量，并根据Kanamori计算震级的经验公式，得到此次地震的地震矩M0为9153 8×1020 N·m，矩震级MW为791，这与USGS（2008）所提供的结果较为一致。
33正演模拟
为了检验该模型的效果，进行了正演计算，并与实测的GPS位移场进行对比，得到的残差如图3所示。从图中可以看出，该方法正演的结果与实测结果运动趋势基本保持一致，整体呈现出以主破裂带为中心的相向运动。模拟残差值的均方误差在东西、南北方向分别为43 cm、52 cm，远场相对于近场的模拟效果更好，断裂带的南段比北段模拟效果更优。其一方面原因可能为本文采用单一断层面进行剖分反演断层滑动分布，而龙门山断裂带地质构造比较复杂，主断裂也是由许多复杂的小断裂组成，把主断裂考虑为一条单一的断层与实际的破裂状况有一定差异，这样在断层模型与实际模型较为一致的区域，正演效果会更优，当在断层模型与实际模型有所差异的区域，正演效果会有所下降，正如本文中断裂带南段比北段模拟效果更优，这也间接说明在断裂带南段的断层模型更接近于实际的断层模型。另一方面原因可能为断层运动所引起的地表形变不仅受到位错的影响，还受到一些地壳区域动力不均匀的影响，离断层越近受到区域动力的影响越大，所以断层远场区域的模拟效果较近场区域会更好，其次也可能由于龙门山断裂带西南端与鲜水河断裂、安宁河断裂等复杂断裂带斜交等因素造成的。
4结论
本文基于GPS数据，利用断层自动剖分技术反演得到汶川地震2条破裂带的滑动分布，经分析研究，可得出以下主要结论：
（1）此次地震以逆冲为主，并带有少量右旋走滑的运动；
（2）从剖分效果来看，映秀―北川断裂带共剖分了301个三角元，灌县―江油断裂剖分了99个三角元，断裂带均剖分的较为细致，能够反映断层中的细部滑动特征，且各个位错元之间的滑动量表现出良好的平滑性；
（3）反演得到映秀―北川断裂最大滑动量为15 m，平均滑动量为261 m，灌县―江油断裂中最大滑动量达8 m，平均滑动量为158 m；
[JP2]（4）计算得到此次地震的地震矩M0为9153 8[JP]×1020 N·m、矩震级MW为791，与USGS所提供的结果较为一致；[JP2]
（5）由残差图看出，正演的结果在运动趋势上与实测结果基本保持一致，各个点位的残差较小，其中在远场区域正演的结果比近场区域会更好。[JP]
综上所述，基于断层自动剖分的三角位错模型能较好地解决汶川地震所引起的同震滑动分布，具有一定的参考价值。
[HTK]本人在撰写论文时得到张永志教授、姜永涛博士、曹海坤和杨九元的帮助，以及王帅博士对程序的支持，在此向他们表示衷心感谢。
参考文献：
陈丹蕾2014基于雷达差分干涉测量与断层自动剖分方法的地震震源参数反演[D].成都：西南交通大学
国家重大科学工程“中国地壳运动观测网络”项目组2008GPS测定的2008年汶川MS80地震的同震位移场[J].中国科学：地球科学，38（10）：1195-1206
季灵运，刘立炜，郝明2015利用InSAR技术研究滇西南镇康—永德地区现今地壳形变特征[J].地震研究，38（1）：84-89
王敏2009基于GPS同震位移场约束反演2008年512汶川大地震破裂空间分布[J].地球物理学报，52（10）：2519-2526
王帅，张永志，姜永涛，等2014维多样性的动态权重粒子群算法反演断层滑动速率[J].地球物理学进展，（4）：1766-1771
王卫民，赵连锋，李娟，等2008四川汶川80级地震震源过程[J].地球物理学报，51（5）：1403-1410
温扬茂，何平，许才军，等2012a联合Envisat和ALOS卫星影像确定L′Aquila地震震源机制[J].地球物理学报，55（1）：53-65
温扬茂，许才军，刘洋，等2012b利用断层自动剖分技术的2008年青海大柴旦MW63级地震InSAR反演研究[J].武漢大学学报（信息科学版），37（4）：458-462
许才军，刘洋，温扬茂2009利用GPS资料反演汶川MW79级地震滑动分布[J].测绘学报，38（3）：195-201，215
张国宏，屈春燕，汪驰升，等2010基于GPS和InSAR反演汶川MW79地震断层滑动分布[J].大地测量与地球动力学，30（4）：19-24
张培震，徐锡伟，闻学泽，等20082008年汶川80级地震发震断裂的滑动速率、复发周期和构造成因[J].地球物理学报，51（4）：1066-1073
张永志，王卫东，魏玉明，等2006利用GPS资料反演祁连山断层的三维滑动速率[J].大地测量与地球动力学，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
