地震发生之后从IRIS分布在全球的地震观测台网获取了震中距在0° ~60° 范围内信噪比较高的宽频带及长周期三分量数字波形资料。为了获得可靠稳定的反演结果, 在反演开始之前需要对下载的原始波形资料进行预处理, 包括去除仪器响应并转化为地面运动速度值, 对数据进行0.006i7~0.01Hz频段内的带通滤波, 旋转三分量记录得到垂向、 径向和切向波形记录。最后对得到的波形资料进行质量控制, 挑选出高质量的观测数据用于反演地震矩张量。

目前大部分的震源机制反演工作利用短周期体波(P波和S波)或面波资料来开展, 然而短周期体波或面波主要在地壳内传播, 容易受到地壳结构不均匀性的影响。本研究使用长周期的W-Phase震相资料, W-Phase震相是地震所产生的非常长周期地震波, 并且以高于S波速的群速度在地球内部传播, 此快速传播的长周期地震波非常适合用来快速、 可靠地确定地震的震源参数(Kanamori et al., 2008)。W-Phase震相首次在1992年的1次尼加拉瓜地震的位移记录中被观测到, 因其具有与众不同的斜坡状长周期波形特征, 被命名为W-Phase震相(Kanamori, 1993)。长周期的W-Phase波形能量主要集中于100~1i000s之内, 以4.5~9.0km/s的相速度传播。该周期范围的地震波能量主要集中在地幔内, 因此W-Phase波在传播过程中较少受到海洋或大陆所引起的强烈浅层各向异性速度结构的影响。

由于W-phase震相具有长周期特性, 故不能直接在频率域使用地震仪的零极点信息去除原始波形资料中的仪器响应, 而是需要把地震仪仪器响应的零极点信息转换为在时间域表示的机械地震仪响应二阶微分方程式的常系数(Zhu, 2003):

z¨+2hω0z˙+ω02z=Gv¨(1)

式(1)中, ω0表示地震仪传感器的自然频率, z表示摆相对于地面的位移, h表示阻尼系数, G表示增益因子, v¨表示地表加速度。随后使用由二阶微分方程表示的时间域递归滤波器, 对原始波形资料进行处理, 去除仪器响应:

vi=2vi-1-vi-2+1G(azi+bzi-1+zi-2)v0=0, v1=0a=1+2hω0Δt+ω02Δt2b=-2(1+hω0Δt)(2)

其中 Δt表示原始资料采样间隔。

获取地面加速度之后, 使用4阶巴特沃斯滤波器对得到的波形进行带通滤波。带通滤波器的频段根据地震震级的大小进行调整(Hayes et al., 2009; Duputel et al., 2012)。最后对滤波后的加速度 v进行2次积分得到地面位移。

通过拟合给定频率范围零极点信息所表示的仪器幅度响应与二阶微分方程式(式(1))所表示的仪器幅度响应, 可以获取可靠的递归滤波器的参数(Kanamori et al., 2008)。对于部分地震台站, 有可能无法在要求的频率范围之内找到满足精度需求的递归滤波器参数。在求取递归滤波器参数过程中直接去除这类台站的波形资料。

为了加快W-Phase反演地震矩张量的速度, 在正式的反演开始之前, 需要预先构建好可供反演过程使用的格林函数库。通常在使用短周期体波或面波开展震源机制反演研究时, 可以使用矩阵传播法或射线叠加的方法来计算理论格林函数, 并且不需要考虑地球重力场对理论格林函数波形的影响。但是对于长周期的W-Phase理论格林函数来说, 地球自身重力场的影响不能忽略。经过考察目前可用的理论格林函数技术程序包, 本研究选择使用基于简振正型叠加(Normal mode summation)方法的Mineos软件包^①(①https://geodynamics.org/cig/software/mineos/。)来构建格林函数库。计算格林函数库时使用一维PREM速度模型(Dziewonski et al., 1981)作为输入的地震波速度模型, 预先计算的格林函数震中距以0.1° 为间隔, 从0° 增加到90° ; 震源深度从0km增加到760km, 随着震源深度的增加, 格林函数库的深度间隔从2km逐渐增加到10km。

发震断层的几何特征与地震波的传播效应共同确定了地震观测波形的变化, 理论上这些效应均可以通过计算得出, 因此可以合成理论地震图。将理论地震图与实际地震观测波形进行对比, 通过二者之间的最佳拟合得到震源机制解。通常地震震源可以用具有6个独立分量的矩张量 Mkl表示(Jost et al., 1989)。地震矩张量 Mkl可以分解为1个各向同性分量、 1个双力偶分量和1个补偿线性偶极子分量(Knopoff et al., 1970)。对于大部分与构造相关的地震, 各向同性部分通常假设为零, 则矩张量只包含偏量成分, 震源为纯剪切位错源, 这时矩张量 Mkl只包含5个独立分量, 这是大部分矩张量反演方法对震源矩张量特征的普遍假设。虽然大部分构造地震都可以用纯偏矩张量来表示, 但部分地震震源过程中存在不可忽略的体积变化, 因此仍然有必要进行完全矩张量的反演。在使用完全矩张量的情况下, 地震台站 i对应的理论地震图可用这6个矩张量分量及其对应格林函数的线性组合来表示:

∑uwiklMkl=uwi(3)

式(3)中, M_kl表示震源矩张量的k-l分量, $u_{wi}^{kl}$表示使用M_kl=1计算的在地震台站i所产生的理论地震图(格林函数), u_wi表示在地震台站i所观测到的W-Phase震相波形。

获取观测波形与地震矩张量的线性方程组之后, 使用雅可比迭代方法求解线性方程, 从而获得所求的地震矩。在雅可比迭代反演中, 首先把格林函数组成的矩阵分解为由对角矩阵 D、严格下三角矩阵L和严格上三角矩阵U表示的等式:

∑uwikl=D-L-U(4)

将式(4)带入式(3), 整理可以得到迭代方程式:

Mkl(n+1)=D-1(L+U)Mkl(n)+D-1uwi(5)

然后给定1组初始的地震矩张量数值, 带入式(5)进行迭代, 直到得到满足收敛条件的1组地震矩张量。

为了得到可靠的结果, 在反演中需要对数据进行多次质量控制, 并利用空间网格搜索方法获取最优解。具体反演工作主要分为2步。 第一步, 使用哈佛大学全球中心矩张量解给出的震源信息作为参考开展初步反演, 反演开始之前首先对原始波形数据进行检查, 去除波形数据中最大振幅与最小振幅之间的过小和过大的反常间隔数据。之后使用剩下的数据再开展3次反演, 每次反演之后, 去除波形拟合均方差分别> 5、 3和0.6的台站数据, 然后再次使用剩下的数据进行反演。经过上述的数据筛选, 共获得56个台站、 88个分量的数据。第二步, 利用初步反演剩下的高质量观测数据, 参考哈佛大学全球中心矩张量解给出的震源参数, 对震源进行网格划分, 开展空间网格搜索, 获得最优的震源信息, 包括震源的地震矩张量解及震源三维位置信息。

图2a为通过空间网格搜索反演获得的最优矩张量解与哈佛大学全球中心矩张量解的对比。本研究获得的结果显示四川九寨沟2017年8月8日M7.0地震的标量地震距为 M0=7.04×1025dyn· cm, 折合矩震级M_W=6.5; 地震矩张量 (Mxx, Myy, Mzz, Mxy, Myz, Mxz)为(0.003i4, 0.473i8, -0.477i2, 0.268i2, 0.231i5, -0.385i5)× 10²⁶dyn· cm。 2个发震断层节面解为: 节面Ⅰ 走向152.7° , 倾角61.4° , 滑动角-4.8° ; 节面Ⅱ 走向245.0° , 倾角85.8° , 滑动角-151.3° 。从图中可以看到, 本研究获得的结果与哈佛大学全球中心矩张量解给出结果比较接近。矩震级差异为0.02, 2个可能的发震断层走向差别较小, 而滑动角、 倾向差异较大, 约20° 。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 a 2017年8月8日九寨沟M7.0地震矩张量反演结果及与哈佛大学全球中心矩张量解对比图; b 部分参加反演台站的波形拟合及台站分布图 黑线表示观测资料, 红线表示理论地震图; 台站名称, 震中距等信息标示于波形图上方; 波形图右边地图显示地震震中位置(蓝色五角星)、 当前波形对应台站的位置(红色圆圈)及所有其他台站的位置(橙色圆圈)Fig.2 Comparison of August 8, 2017 Jiuzhaigou earthquake moment tensor solutions coming from global CMT and this study(a)and comparison of the observed and synthetic records for moment tensor solution(b).

图2b为部分参加反演台站的理论波形与观测波形的对比, 其中红色为理论波形, 黑色为实际观测资料, 可见大部分台站的理论波形与实际观察资料拟合度较好。图2b中每个波形对比图右边的地图显示反演中使用的地震台站分布(橙色圆点)、 当前波形图对应的台站位置(红色圆点)及震中位置(蓝色五角星)。从图中可以看出大部分地震台站位于震中的东部和南部, 同时中亚、 中东及部分欧洲台站覆盖了西部、 西北方向的区域。总体来说, 台站的方位角覆盖较为均匀, 为反演结果的可靠性提供了坚实的数据基础。

在深度空间网格搜索反演中, 将深度网格搜索步长设置为2km, 起始深度为2km, 最大深度为24km, 共搜索11个深度位置。图3a为波形拟合均方差随震源深度的变化情况, 此次地震的最佳矩心深度为14km, 与美国地质调查局(USGS)和GCMT给出的震源深度位置比较接近。图3b显示了在深度为14km的平面上拟合均方差随空间位置的变化, 其中 “ +” 标示出最优震中的位置。

图3

Fig. 3

	Figure Option ViewDownloadNew Window
	图3 2017年8月8日九寨沟M7.0地震震源深度及震中位置空间网格搜索结果 a 深度搜索与波形拟合均方差变化图; b 经纬度搜索与波形拟合均方差分布图。 黑色十字标示最佳震中位置Fig. 3 Results of optimal centroid location grid search of August 8, 2017 Jiuzhaigou earthquake.

九寨沟M7.0地震发生之后, 国内外多家机构在第一时间发布了地震源参数, 相关研究机构及本文获得的结果均列于表1中。其中美国地质调查局(USGS)发布结果与哈佛大学全球中心矩张量解(GCMT)比较接近, 断层的性质为正斜滑断层; 中国地震局地球物理研究所(CEA-IG)和中国地震局地震预测研究所(CEA-IES)的结果与USGS和GCMT的结果差异比较大。本文获得的断层面解、 震源深度与估计震级都接近USGS和GCMT的结果。总体而言, 本次地震的走滑分量比较强, 可视为1个近走滑型地震。各个研究机构甚至同一个研究机构的反演结果存在的差异, 来自于与地震矩张量反演相关的不确定性, 其中与观测波形资料相关的不确定性包括选取的震相(P波、 SH波、 面波或W-Phase震相)、 滤波带宽以及波形资料的方位角覆盖等。另外, 体波或面波资料相对W-phase震相频率较高, 对于体波而言地壳速度模型的相对微小变化有可能会对理论地震图产生显著的影响。各相关研究机构获得的震源参数有些是使用体波或面波资料反演的, 且合成理论地图使用的地壳速度模型也可能不一致, 最终导致公布的发震断层面参数存在一定差异。

表1

Table 1

表1(Table 1)

表1 各个研究机构得到的2017年8月8日九寨沟M7.0地震震源机制解 Table 1 Focal mechanism of August 8, 2017 Jiuzhaigou earthquake published by different institutes 来源 节面 走向/(° ) 倾角/(° ) 滑动角/(° ) 震源深度/km 震级MW USGS① Ⅰ 246 57 -173 13.5 6.5 Ⅱ 153 84 -33 GCMT② Ⅰ 242 77 -168 14.9 6.5 Ⅱ 150 78 -13 CEA-IG③ Ⅰ 328 48 -11 5 6.5 Ⅱ 65 82 -137 CEA-IES④ Ⅰ 152 75 13 5 6.5 Ⅱ 59 77 164 本文 Ⅰ 245.0 85.8 -151.3 14 6.5 Ⅱ 152.7 61.4 -4.8 注 ①美国地质调查局; ②哈佛大学全球中心矩张量解; ③中国地震局地球物理研究所; ④中国地震局地震预测研究所。

表1 各个研究机构得到的2017年8月8日九寨沟M7.0地震震源机制解 Table 1 Focal mechanism of August 8, 2017 Jiuzhaigou earthquake published by different institutes

使用与研究九寨沟地震相同的方法反演精河M6.6地震。经过初步反演及波形资料筛选, 共获得60个台站、 90个分量的波形数据进行下一步的空间网格搜索反演。

图4a为本研究获得的地震矩张量与哈佛大学全球中心矩张量解的对比。结果显示新疆精河2017年8月9日M6.6地震标量地震距为M₀=4.04× 10²⁵dyn· cm, 折合矩震级M_W=6.34, 地震矩张量(Mxx, Myy, Mzz, Mxy, Myz, Mxz)为(0.311i7, -0.309i2, -0.002i5, -0.233i9, -0.110i9, -0.010i1)× 10²⁶dyn· cm。 2个发震断层节面解为: 节面Ⅰ 走向100.6° , 倾角27.5° , 滑动角114.1° ; 节面Ⅱ 走向253.9° , 倾角65.1° , 滑动角78.0° 。从图中可以看到, 本研究获得结果与哈佛大学全球中心矩张量解比较接近。矩震级差异为0.04, 可能的发震断层面走向差异较小, 但倾角和滑动角的差异较大, 约20° 。

图4

Fig. 4

	Figure Option ViewDownloadNew Window
	图4 a 2017年8月9日新疆精河M6.6地震矩张量反演结果与哈佛大学全球中心矩张量解对比图; b 部分参加反演台站的波形拟合及台站分布图 黑线表示观测资料, 红线表示理论地震图; 台站名称、 震中距等信息标示于波形图上方; 波形图右边地图显示地震震中位置(蓝色五角星)、 当前波形对应台站的位置(红色圆圈)及所有其他台站的位置(橙色圆圈)Fig. 4 Comparison of August 9, 2017 Xinjiang Jinghe earthquake moment tensor solution coming from global CMT and this study(a) and comparison of the observed and synthetic records for moment tensor solution(b).

图4b为部分参加反演台站的理论波形和观测波形的对比, 图中显示理论波形与实际观察资料拟合度较好。从波形对比图右边的地图可以看出, 除震中南部地震台站较少之外, 其他3个方向台站的方位角覆盖都较为均匀, 可以保证结果的可靠性。

在深度空间网格搜索中, 将深度网格搜索步长设置为2km, 起始深度为2km, 最大深度为34km, 共搜索16个深度位置。图5a为波形拟合均方差随震源深度的变化情况, 此次地震的最佳矩心深度为16km, 与USGS和GCMT给出的震源深度位置存在明显差异。图5b显示在深度为16km的平面上进行空间网格搜索时, 波形拟合均方差随空间位置的变化, 其中 “ +” 标示出最优震中位置。

图5

Fig. 5

	Figure Option ViewDownloadNew Window
	图5 2017年8月9日新疆精河M6.6地震震源深度及震中位置空间网格搜索结果 a 深度搜索与波形拟合均方差变化图; b 经纬度搜索与波形拟合均方差分布图。 黑色十字标示最佳震中位置Fig. 5 Results of optimal centroid location grid search of August 9, 2017 Xinjiang Jinghe earthquake.

表2列出了各相关研究机构发布的2017年8月9日新疆精河M6.6地震的震源机制解。从表中可以看出, 各相关研究机构对矩震级M_W的估计比较接近, 为6.2~6.4。但是对于震源深度的估计存在较大差异, 从约16km到约31km。各相关研究机构对发震断层2个节面的参数估计差异不大, 断层面的走向、 倾角及滑动角的差异都≤ 35° 。总体而言, 目前发布的该次地震的发震断层面都呈近EW向, 大体上与北天山的走向接近, 震源的性质与区域地质构造背景一致。新疆精河地震是发生在中下地壳的1次逆冲型破裂事件。

表2

Table 2

表2(Table 2)

表2 各个研究机构得到的2017年8月9日新疆精河M6.6地震震源机制解 Table 2 Focal mechanism of August 9, 2017 Xinjiang Jinhe earthquake published by different institutes 来源 节面 走向/(° ) 倾角/(° ) 滑动角/(° ) 震源深度/km 震级MW USGS① Ⅰ 269 30 87 30.5 6.3 Ⅱ 92 60 92 GCMT② Ⅰ 255 41 75 27.9 6.3 Ⅱ 94 50 103 CEA-IG③ Ⅰ 262 45 80 20 6.2~6.3 Ⅱ 96 46 100 CEA-IES④ Ⅰ 272 44 95 17 6.2 Ⅱ 85 56 85 本文 Ⅰ 253.9 61.5 78.0 16 6.34 Ⅱ 100.6 27.5 114.1 注 ①美国地质调查局; ②哈佛大学全球中心矩张量解; ③中国地震局地球物理研究所; ④中国地震局地震预测研究所。

表2 各个研究机构得到的2017年8月9日新疆精河M6.6地震震源机制解 Table 2 Focal mechanism of August 9, 2017 Xinjiang Jinhe earthquake published by different institutes