九寨沟M7.0地震多个震源机制中心解的确定

马付红，李莹甄

(1.四川省地震局西昌地震中心站，四川 西昌 615000；2.防灾科技学院，河北 燕郊 065201)

2017年8月8日21时19分，在四川省阿坝藏族羌族自治州九寨沟县发生M7.0地震，地震震中位于(33.20°N,103.82°E),震源深度20 km(据中国地震台网)。此次地震造成25人死亡，525人受伤，6人失联， 176 492人受灾，73 671间房屋受损；其最大烈度达Ⅸ度，Ⅵ度及以上总面积为18 295 km2(李志强等，2017)。大地震发生后人们除了关注其发震位置和破坏程度外，还关注其发震的物理过程，错动方式，发震断层面和破裂方向等，而这些物理量的确定可以用震源机制解进行描述。九寨沟M7.0地震发生后，多家研究机构和个人对其用不同的资料和方法得出了不同的震源机制，不同震源机制之间存在一定的差别，为此有必要进行多个震源机制中心解的确定，以便更准确的确定出此次地震的震源机制最优解。最优的震源机制解，将对后续的应力分析，应力场研究，地震应力前兆分析，地震静态应力触发(万永革，2001)，大地震对周围断层影响(万永革等，2009)等研究提供帮助。多个震源机制中心解是指对于几个震源机制测量结果，确定一个震源机制，使所有震源机制与其旋转角平方和最小(万永革，2019)。目前计算同一地震多个震源机制中心解的方法除了有Levenberg─Marquadt算法(Levenberg，1944；Marquardt，1963；万永革，2019)，还可以用旋切法、牛顿-拉普森法(Press et al,1992；万永革，2019)、二分法等。求震源机制中心解首先要找到表达两个震源机制差别的方法，本文用两个震源机制空间旋转角差别的方法来表示。表达两个震源机制解空间旋转角的差别是非线性问题，引入Levenberg─Marquadt算法将非线性问题转化成线性问题进行迭代求解(Levenberg，1944；Marquardt，1963；万永革，2019)。

1 发震背景

九寨沟M7.0地震发生在塔藏断裂、虎牙断裂、岷江断裂和雪山梁子断裂所围限的区域(赵博等，2018)。塔藏断裂为全新世左旋走滑断裂，虎牙断裂为全新世逆断裂，岷江断裂为西倾的全新世逆断裂，岷江断裂和塔藏断裂分别为巴颜喀拉地块的东、北边界。本次地震发震断层以左旋走滑运动为主，是巴颜喀拉块体边界断裂持续活动的结果，主要表现为印支板块向北挤压，导致川西北地块挤出(薄景山等，2018)。结合该区域地质构造背景和震源机制解，其发震断层为塔藏断裂南段，并与虎牙断裂北端隐性趋于贯通(姚鑫等，2017；易桂喜等，2017)。

2 方法原理

获取同一地震多个震源机制的中心解，要求所有测量震源机制与中心解的空间旋转角平方和最小。在建立直角坐标系后最小空间旋转角可以将两个震源机制进行比较，把他们当中的一个震源机制的B、P、T轴通过坐标旋转到另一个震源机制的B、P、T轴上，在空间上旋转的角度值大小代表了两个震源机制之间的差异(万永革，2019)。所有测量震源机制与所求中心解的最小空间旋转角平方和最小的表达式为：

(1)

在(1)式中(θci)min表示最小空间旋转角，采用Levenberg-Marquardt方法进行迭代求解(Levenberg，1944；Marquardt，1963；万永革，2019)。进行迭代求解时选定的初始解接近于估计的解，这样估计解就可以表达为初始解附近的泰勒展开，忽略二阶及以上导数项可以得：

(2)

在(2)式中，vi是第i个震源机制与所求中心解之间最小空间旋转角的平方和。当震源机制总数为N时，将上式写成矩阵表达式得：

(3)

将(3)式左边的第一个雅克比矩阵设为J，第二个矩阵(在初始解基础上改变的参数量)设为x，同时把右边矩阵设为d，则变量x表示为：

x=(JTJ)-1JTd

(4)

对(4)式进行求解，则可得到改变量x(即Δφ，Δδ，Δλ)。令

φ0=φ0+Δφ，δ0=δ0+Δδ,λ0=λ0+Δλ

(5)

把(5)式代入(3)式，反复计算，直到解的改变量很小，然后终止迭代，这样就可以得到非线性问题求解后的最优解。为了让求解快速收敛的同时，防止雅克比矩阵产生奇异，采用Levenberg-Marquardt算法时加上一个可变的阻尼，把(4)式改为：

x=(JTJ+κI)-1JTd

(6)

式中κ为计算中可变的数，I为单位矩阵。把测量的各个震源机制与中心震源机制的最小空间旋转角的标准差作为计算的最终误差范围，其标准差表达式为：

(7)

在(7)式中s为所有测量震源机制与中心解的最小空间旋转角的平方和。

3 数据与处理

本文收集了九寨沟M7.0地震后美国地质勘探局(USGS)(https://earthquake.usgs.gov/earthquakes/eventpage/us2000a5x1/moment-tensor)、全球质心矩张量项目(gCMT)(https://www.globalcmt.org/)、中国地球物理研究所(IG-CEA)、中国地震台网中心(CENC)(http://www.cenc.ac.cn/eportal/fileDir/cenc/resource/cms/2018/06/2018060811421545641.pdf)、中国地震预测研究所(IEF-CEA)、易桂喜等(2017)、谢祖军等(2018)、杨宜海等(2017)、刘旭宙(2017)、皱立晔(2018)、梁姗姗等(2018)、王莹等(2017)以及韩立波等(http://www.cea-igp.ac.cn/Uploads/image/20170811/2017年8月8日四川九寨沟7-v1.0.jpg)机构和研究人员的震源机制解(见表1)。根据求中心解的原理，选择不同的初始解，将产生不同的震源机制和标准误差。本研究均以各机构和个人的震源机制作为初始解，把每次计算所得的目标函数进行比较，选标准误差最小的作为中心解。以研究机构和个人产出的震源机制解作为初始震源机制解进行计算得出的中心震源机制结果及标准差见表2。从表2可以看出不管以哪个机构或个人得到的震源机制为初始震源机制，得到的中心震源机制解的差别都很小，标准差到小数点后4位都是相同的。根据标准误差的描述以韩立波等的震源机制作为初始解，计算得到的最小三维空间旋转角的标准差最小(表2第4列14行)，其值为23.12492。通过求中心解的程序得出本研究的中心解节面Ⅰ的走向152.52°、倾角84.40°、滑动角-3.95°，节面Ⅱ的走向242.91°、倾角86.07°、滑动角-174.39°。压应力轴P走向107.78°，不确定范围为85.53°～130.53°，倾伏角6.75°，不确定范围为-9.32°～22.31°。张应力轴T走向17.64°，不确定范围为-4.61°～40.39°，倾伏角1.17°，不确定范围为-14.78°～17.20°。中间应力轴B走向为277.83，不确定范围为168.51°～409.34°，倾伏角83.15°，不确定范围为61.90°～86.05°。该地震的中心解(a)和空间三维辐射花样(b)见图1。

表1 九寨沟7.0级地震多个震源机制解

表2 九寨沟7.0级地震的中心解、标准误差和最小旋转角结果

图1 九寨沟地震的中心震源机制解(a)及空间三维辐射花样(b)(a)中的粗黑色弧线表示中心震源机制的两个节面，淡灰色细弧线覆盖区域为其不确定范围；黑色的圆点表示中心震源机制解的P轴、T轴和B轴，其周围对应颜色的封闭曲线表示其不确定性的范围，雪花、小圆圈和加号的点表示各个机构得到的震源机制解P轴、T轴和B轴，黑色细弧线表示各个机构和作者得到的震源机制节面。(b)中的压缩区域和膨胀区域分别用黑色和灰色表示。

4 结论与讨论

(1)本文收集了多家研究机构和个人用不同资料和方法得出的多个震源机制解，分别以他们的震源机制解作为初始解进行计算比较。以标准误差最小的韩立波等的震源机制为初始解，得出九寨沟M7.0地震的中心解节面Ⅰ的走向152.52°、倾角84.40°、滑动角-3.95°，节面Ⅱ的走向242.91°、倾角86.07 °、滑动角-174.39°。压应力轴P走向107.78°，不确定范围为85.53°～130.53°。张应力轴T走向17.64°，不确定范围为-4.61°～40.39°。中间应力轴B走向为277.83，不确定范围为168.51°～409.34°。(2)九寨沟M7.0地震中心解的确定，首先要找出两个震源机制解空间旋转角差别的表达方法，再引入Levenberg─Marquadt算法使非线性问题转化成线性问题进行迭代求解。该中心解不是多个震源机制的平均值，而是所有测量震源机制与中心解的最小空间旋转角的平方和最小为标准解。(3)之前机构和个人得出九寨沟M7.0地震的多个震源机制解之间存在一定的差别，用求中心解的原理，求出该地震震源机制最优解，为该区域进行应力及应力场、地震应力前兆分析、地震静态应力触发和大地震对周围断层影响等研究发挥作用。(4)后续可以用本文提及的其他求中心解的方法，求出该地震的中心解，并与本方法求得的中心解进行比较研究。

致谢：感谢防灾科技学院万永革教授提供了求中心解的程序。