地震震源研究对于判定发震构造、分析发震机理、研究地震动力学等有着重要的意义(倪四道, 2008; 王卫民等, 2008; 许力生和陈运泰, 2004; 赵翠萍等, 2009; Ye et al., 2016).震源研究的一个重要手段是分析地震波形数据，利用振幅获取地震震级大小，利用波形反演获取点源机制解，以及基于有限断层模型反演震源破裂过程等(韩立波和蒋长胜, 2012; 刘成利等, 2015; 王卫民等, 2008; 张勇等, 2008; Hao et al., 2013; Yue et al., 2013).

在震源研究中，震源附近的区域宽频带地震数据尤为重要，特别是应急地震学中，需要实时获取震源机制解等(崇加军等, 2010; 韩立波和蒋长胜, 2012).限幅的地震数据有无法用面波振幅评价震级，无法用长周期波形获取震源机制等不足.而区域地震台网面临着近震中距地震数据限幅的问题，例如，对于2013年芦山M_W6.6地震，震中距远至~370 km，仍存在数据限幅的问题(Zhang et al., 2016).

针对这一问题，Zhang等(2016)发展了一种利用改进的POCS法恢复地震数据的方法，其基本原理是利用合理的约束求解最优解集.针对地震数据限幅问题，依据两个基本约束(限幅数据在时间域的振幅须超越限幅阈值；频率域的高能量频谱可信)交替地在时间域和空间域不断迭代，同时缓慢降低频域阈值，最终可以有效地修复限幅数据.数值实验结果表明该方法精度非常高(Zhang等2016)，即使限幅发生在峰值振幅30%时仍能获得可靠的恢复结果，同时利用该方法恢复了2013年芦山M_W6.6地震的区域宽频带地震记录.本研究利用恢复的区域宽频带地震波形确定芦山地震的面波震级、点源机制解和震源破裂过程，用于评估该方法在实际应用中的可行性.

Russell(2006)提出了一种在区域和远震距离测量面波震级的时域方法.震级公式是

(1)

式中, 为滤波器的拐角频率，a_b, Δ和T分别为峰值振幅(测量单位为mm)，震中距(单位为°)和周期，该震级公式也被应用于核爆震级的测定(范娜等, 2013).本研究采用Russell(2006)的方法确定面波震级，由于垂向和切向分量具有瑞利波，且垂向分量通常具有较高的信噪比，因此利用震中距为99~300 km的区域宽带波形的25个垂向分量来确定震级(图 1a).在去除仪器响应后，将地震观测转换为位移记录，并使用零相位的三阶巴特沃斯滤波器对波形进行滤波，使用18个滤波频带，滤波周期自8 s起，间隔1 s至25 s.用群速度5.5~1.8 km·s^-1的时间窗得到瑞利波波形，测量零到峰值的振幅后，根据该公式得到震级，取18个频带中的最大震级作为面波震级.基于恢复前数据(含限幅数据)和恢复后数据确定的震级分别表示为M_S_clipped和M_S_restored.

图 1

Fig. 1

图 1 用于计算震级的台站分布以及震级差异.(a)台站分布图，五角星和三角形分别表示震中位置以及台站分布; (b)由恢复前和恢复后数据所获取的震级之差随震中距的变化图 Fig. 1 The distribution of stations used to investigate the magnitude and the difference of magnitudes.(a)The distribution of stations. The star and triangles show the location of epicenter and the stations, respectively; (b) The difference between the magnitudes based on clipped data and that based on restored data with epicentral distances

震级结果如表 1所示，恢复后震级和恢复前震级之差如图 1b所示，垂向位移记录的峰值位于周期12~25 s之间.结果显示，除了最近的四个台站(99~107 km)，面波震级差异较大，恢复后震级比恢复前震级大0.13~0.45级.其余震级差异较小，而当震中距超过250 km之后，恢复前后震级几乎无差异.基于恢复前数据获取的面波震级为7.01，由恢复后数据获得的震级为7.06.

表 1 (Table 1)

表 1 基于恢复前和恢复后数据的震级 Table 1 The magnitudes based on clipped data and restored data 台站名 震中距 (km) 最大振幅周期(s) 恢复后震级 MS_restored 恢复前震级 MS_clipped ΔMS XJI 99 25.0 7.5582 7.2694 0.2888 CD2 99 25.0 7.3969 6.9481 0.4488 JYA 106 25.0 7.3994 7.2707 0.1287 YGD 107 25.0 7.5845 7.4078 0.1767 SMI 135 23.0 7.1512 7.1392 0.0120 JJS 167 13.0 6.9923 6.9799 0.0124 MBI 170 17.0 6.9552 6.8947 0.0605 MEK 193 17.0 7.0839 7.0936 -0.0097 YJI 194 16.0 6.8597 6.8680 -0.0083 AXI 202 19.0 7.0058 6.9728 0.0330 JLO 204 23.0 6.9286 6.9470 -0.0184 ZJG 230 16.0 6.9272 6.9393 -0.0121 LBO 231 25.0 7.2607 7.2192 0.0415 MNI 232 20.0 6.8479 6.8646 -0.0167 HWS 249 15.0 7.0333 7.0397 -0.0064 ROC 257 15.0 7.0599 7.0585 0.0014 YAJ 270 25.0 7.0619 7.0675 -0.0056 XSB 275 21.0 6.8928 6.8920 0.0008 BYD 277 22.0 6.8436 6.8274 0.0162 PWU 277 23.0 6.7735 6.7691 0.0044 JLI 277 16.0 6.9233 6.9382 -0.0149 YYC 281 21.0 6.8114 6.8102 0.0012 LZH 282 12.0 7.0510 7.0617 -0.0107 XCO 289 13.0 6.9893 6.9904 -0.0011 RTA 291 21.0 7.0672 7.0649 0.0023

表 1 基于恢复前和恢复后数据的震级 Table 1 The magnitudes based on clipped data and restored data

区域地震资料被广泛用于反演中小地震震源机制解(韩立波和蒋长胜, 2012; 罗艳等, 2010; 易桂喜等, 2019; Hao and Yao, 2012; Zhao and Helmberger, 1994; Zhu and Helmberger, 1996).其具体流程包括数据处理、格林函数计算以及波形反演三个部分.本项研究中将区域宽频带地震记录去仪器响应变换到速度记录，并进行带通滤波，滤波器频带0.02~0.06 Hz.基于Crust1.0地壳速度模型(表 2)(https://igppweb.ucsd.edu/~gabi/crust1.html)，使用广义反射透射系数矩阵和离散波数方法计算理论地震图(Yao and Harkrider, 1983).基于模拟退火反演方法(姚振兴和纪晨, 1997; 张霖斌等, 1997)，通过拟合长周期波形获取双力偶点源机制解，包括走向、倾角、滑动角、深度以及矩张量，误差函数定义为

表 2 (Table 2)

表 2 地壳模型 Table 2 Crust model P波速度VP (km·s-1) S波速度VS (km·s-1) 密度ρ (g·cm-3) 层厚度 (km) 2.50 1.20 2.10 1.0 4.00 2.10 2.40 1.0 6.10 3.50 2.75 20.0 6.30 3.60 2.80 20.0 7.20 4.00 3.10 4.0 8.08 4.47 3.38

表 2 地壳模型 Table 2 Crust model

(2)

其中，Err为误差函数，N_s表示台站个数，N_ci表示第i个台站使用到的分量个数，S_{(ij, m)}和O_{(ij, m)}分别为为第i个台站第j分量的观测速度记录以及理论合成记录，mpoint为地震记录的采样点数.

选取震中距小于400 km的区域宽频带地震记录中没有限幅的数据，同时兼顾方位角覆盖(图 2a中的黑色三角形)，利用这些数据反演芦山地震的震源机制解，将之定义为参考机制解，得到的结果为走向、倾向和滑动角分别为213°、45°和94°，震源深度15 km，地震矩张量7.5×10¹⁹N·m(表 3)，波形拟合结果如图 3所示.选取震中距150 km至250 km间七个数据限幅的台站(图 2a中的红色三角形)，分别使用恢复前的数据和恢复后的数据进行反演，结果分别记作恢复前结果和恢复后结果(表 3)，图 2b显示了恢复前结果和恢复后结果的误差函数及震源机制随深度的变化关系，蓝线串起的沙滩球表示恢复前数据的结果，红线串起的沙滩球表示恢复后数据的结果.未填充颜色的大沙滩球中的黑线、蓝线和红线分别表示参考机制解、恢复前机制解和恢复后机制解的一个节面.相应的波形拟合结果如图 4和图 5所示.

图 2

Fig. 2

图 2 台站分布与点源机制解比较.(a)台站分布图；红色五角星表示了震中位置，黑色三角形以及黑色台站名表示用来获取参考机制解的台站分布，红色三角形以及红色台站名表示使用到的含有限幅数据的台站分布.(b)基于恢复前和恢复后数据，震源机制与误差函数随深度的变化；红色的沙滩球表示恢复后数据的结果，蓝色的沙滩球表示恢复前数据的结果；未填充颜色的大沙滩球中的黑线、蓝线和红线分别表示参考机制解、恢复前机制解和恢复后机制解的一个节面 Fig. 2 The distribution of stations and comparison of focal mechanisms. (a) The distribution of stations. The red star and black triangles show the location of epicenter and the stations used to get the reference mechanism, respectively. The red triangles indicate the used stations with clipped data. (b) The focal mechanisms and misfits in different depths based on the clipped and restored data. The beach balls connected by the red and blue line represent the results of restored data and clipped data, respectively. The black, blue and red line in the large beach ball represent a nodal plane of the reference mechanism, clipped mechanism and restored mechanism, respectively

表 3 (Table 3)

表 3 点源机制解反演比较 Table 3 Comparison of focal mechanisms 走向 (°) 倾角 (°) 滑动角 (°) 深度 (km) 矩张量 (N·m) 恢复前结果 203 50 81 14 8.1×1019 恢复后结果 210 45 90 14 8.3×1019 参考机制解 213 45 94 15 7.5×1019

表 3 点源机制解反演比较 Table 3 Comparison of focal mechanisms

图 3

Fig. 3

图 3 参考机制解的波形比较.细线和粗虚线分别表示实际记录和理论波形，每道数据左边为台站名，N、E和Z分别表示北向、东向和垂向分量，左边数字上下分别表示方位角和震中距(单位km)，右边的数字表示实际记录的峰值(单位cm·s-1)，波形使用实际记录的峰值进行了归一化 Fig. 3 Comparison of waveforms predicted using reference mechanism. Comparison of three-component velocity records (thin solid line) and synthetic seismograms (thick dashed line) predicted using reference mechanism. The station name is indicated at the left side of each seismogram. N, E and Z indicate the North-south component, East-west component and Vertical component, respectively. The value above the beginning of each trace is the source azimuth in degrees, and below is the epicentral distance in kilometers. The peak velocity of the observation in cm/s is indicated above the end of each trace, which is used to normalize both data and the corresponding synthetic seismogram

图 4

Fig. 4

图 4 恢复后数据在最佳机制解下的波形比较.细线和粗虚线分别表示实际记录和理论波形，每道数据左边为台站名，N、E和Z分别表示北向、东向和垂向分量，左边数字上下分别表示方位角和震中距(单位km)，右边的数字表示实际记录的峰值(单位cm·s-1)，波形使用实际记录的峰值进行了归一化 Fig. 4 Comparison of waveforms predicted by preferred mechanism based on restored data. Comparison of three-component velocity records (thin solid line) and synthetic seismograms (thick dashed line) predicted using reference mechanism. The station name is indicated at the left side of each seismogram. N, E and Z indicate the North-south component, East-west component and Vertical component, respectively. The value above the beginning of each trace is the source azimuth in degrees, and below is the epicentral distance in kilometers. The peak velocity of the observation in cm·s-1 is indicated above the end of each trace, which is used to normalize both data and the corresponding synthetic seismogram

图 5

Fig. 5

图 5 恢复前数据在最佳机制解下的波形比较.细线和粗虚线分别表示实际记录和理论波形，每道数据左边为台站名，N、E和Z分别表示北向、东向和垂向分量，左边数字上下分别表示方位角和震中距(单位km)，右边的数字表示实际记录的峰值(单位cm·s-1)，波形使用实际记录的峰值进行了归一化 Fig. 5 Comparison of waveforms predicted by preferred mechanism based on clipped data. Comparison of three-component velocity records (thin solid line) and synthetic seismograms (thick dashed line) predicted using reference mechanism. The station name is indicated at the left side of each seismogram. N, E and Z indicate the North-south component, East-west component and Vertical component, respectively. The value above the beginning of each trace is the source azimuth in degrees, and below is the epicentral distance in kilometers. The peak velocity of the observation in cm·s-1 is indicated above the end of each trace, which is used to normalize both data and the corresponding synthetic seismogram

利用恢复前地震数据得到的机制解中，滑动角差别最大，与参考机制解相比，差别达到了13°；而利用恢复后的地震数据得到的机制解，这一差别仅为4°.相应的震源深度和标量地震矩的差别均不大.综上所述，使用震中距较近台站的恢复后地震数据可以获取稳定可靠的震源机制解，对发震构造性质的判定具有重要意义，且数据恢复时间近实时，故可应用于地震应急中，近实时获取震源机制解.

Hao等(2013)采用远场地震波形联合近场强地面运动数据，反演获得了芦山地震的破裂过程.本研究中采用恢复后的地震数据，选取震中距较近的6个台站的三分量记录(震中距81~134 km)，联合远场体波与面波波形，反演芦山地震的破裂过程.比较两个破裂过程，评价利用恢复后数据获取破裂过程结果的有效性.采用和Hao等(2013)完全相同的反演方法、反演参数范围、断层几何形态(走向210°，倾角40°，长度42 km，宽度38 km，子断层尺度2 km×2 km)以及震源位置(30.3°E，103°N，12.5 km).将恢复后的地震数据，去仪器响应至速度记录，进行带通滤波，滤波频带为0.02~0.5 Hz，联合远场体波(0.0033~1 Hz)与面波波形(4~6 mHz)反演震源破裂过程.

震源破裂过程结果比较以及7个区域台站的分布如图 6所示，主要震源参数的比较如表 4所示，波形拟合如图 7—9所示.结果显示两者滑动量较大的区域，位置较为一致，位于深度约6 km至21 km，震中两侧10 km范围内.本文得到的主要震源破裂过程参数(如破裂速度2.1 km·s^-1，平均滑动量0.4 m等)和Hao等(2013)得到的结果都十分接近，验证了恢复数据的可用性，而直接使用限幅数据做震源破裂过程反演得到的结果与两者存在显著差异(图 6).也可以注意到，恢复后数据得到的破裂过程，细节上与Hao等(2013)的结果存在一些差异，这主要是由于使用的数据集不同所造成的.Hao等(2013)使用更多且更近的波形数据，因此可以提供更多细节信息；本文所用台站的震中距多在100 km左右，且有四个台站位于四川盆地内部(图 6)，由于受到盆地复杂三维结构的影响，地震波形复杂，难以很好地拟合(图 9).盆地里的台站波形拟合较差与Hao等(2013)的结果是一致的.

图 6

Fig. 6

图 6 破裂过程、台站分布与地震矩释放率.(a)强地面运动加远场波形数据获取的破裂过程(Hao et al., 2013), 颜色、箭头和等值线分别表示滑动量的大小、滑动方向以及破裂传播时间.(b)恢复后的区域宽频带数据加远场数据获取的破裂过程.(c)恢复前的区域宽频带数据加远场数据获取的破裂过程.(d)台站分布图.红色五角星表示了震中位置，黑色三角形表示台站分布，黑色矩形框为断层在地表的投影.(e)地震矩释放率随时间变化的比较 Fig. 6 Rupture process, distribution of stations and moment rate. (a) Cross sections of slip distributions of the model constrained by the strong motion and teleseismic data. The red star indicates the hypocenter. The color and the white arrows denote the slip amplitude and direction, respectively. The contours show the rupture initiation time. (b) Cross sections of slip distributions of the model constrained by the restored regional broadband waveforms and teleseismic data. (c) Cross sections of slip distributions of the model constrained by the clipped regional broadband waveforms and teleseismic data. (d) The distribution of stations. The red star and black triangles show the location of epicenter and the used stations, respectively. The black box shows the surface projection of the fault plane used in this study. (e) Comparison of moment rate functions

表 4 (Table 4)

表 4 震源参数对比 Table 4 The comparison of source parameters 参数 Hao等(2013) 本文 标量地震矩 8.9×1018Nm 8.4×1018Nm 矩心深度 12.5 km 13.6 km 峰值滑动深度 15.1 km 13.8 km 主要滑动范围 ~5-20 km ~6-21 km 峰值滑动量 1.2 m 1.1 m 平均滑动量 0.4 m 0.4 m 破裂持续时间 10.3 s 9.6 s 破裂速度 2.0 km·s-1 2.1 km·s-1

表 4 震源参数对比 Table 4 The comparison of source parameters

图 7

Fig. 7

图 7 远场体波波形比较.细线和粗虚线分别表示实际记录和理论波形，每道数据左边为台站名，P和SH分别表示垂向P波速度记录和切向S波位移记录，左边数字上下分别表示方位角和震中距(单位°)，右边数字表示实际记录的峰值(单位为cm·s-1或cm)，波形使用实际记录的峰值进行了归一化 Fig. 7 Comparison of teleseismic body waveforms. Comparison of teleseimic body waveforms (thin solid line) and synthetic seismograms (thick dashed line). The station name is indicated at the left side of each seismogram. P and SH indicate the Vertical P wave component and Tangential SH component, respectively. The value above the beginning of each trace is the source azimuth in degrees, and below is the epicentral distance in degrees. The peak velocity of the observation in cm·s-1 is indicated above the end of each trace, which is used to normalize both data and the corresponding synthetic seismogram

图 8

Fig. 8

图 8 远场面波波形比较.细线和粗虚线分别表示实际记录和理论波形，每道数据左边为台站名，UD和SH分别表示垂向记录和切向记录，左边数字上下分别表示方位角和震中距(单位°)，右边的数字表示实际记录的峰值(单位mm)，波形使用实际记录的峰值进行了归一化 Fig. 8 Comparison of teleseismic surface waveforms. Comparison of teleseimic surface waveforms (thin solid line) and synthetic seismograms (thick dashed line). The station name is indicated at the left side of each seismogram. UD and SH indicate the Vertical component and Tangential component, respectively. The value above the beginning of each trace is the source azimuth in degrees, and below is the epicentral distance in degrees. The peak velocity of the observation in mm is indicated above the end of each trace, which is used to normalize both data and the corresponding synthetic seismogram

图 9

Fig. 9

图 9 区域宽频带波形记录比较.细线和粗虚线分别表示实际记录和理论波形，每一道数据左边为台站名，N、E和Z分别表示北向、东向和垂向分量，左边数字上下分别表示方位角和震中距(单位°)，右边的数字表示实际记录的峰值(单位为cm·s-1)，波形使用实际记录的峰值进行了归一化 Fig. 9 Comparison of regional broadband waveforms. Comparison of three-component velocity records (thin solid line) and synthetic seismograms (thick dashed line). The station name is indicated at the left side of each seismogram. N, E and Z indicate the North-south component, East-west component and Vertical component, respectively. The value above the beginning of each trace is the source azimuth in degrees, and below is the epicentral distance in kilometers. The peak velocity of the observation in cm·s-1 is indicated above the end of each trace, which is used to normalize both data and the corresponding synthetic seismogram

采用芦山地震震中距位于99~300 km台站的恢复后数据测定面波震级.结果显示，100 km震中距左右的台站(99~107 km)修复后震级与修复前震级存在显著差异，对于更远的台站(震中距大于135 km)，两者的差异小于0.061级，差异不大.同时，由于严重限幅的地震数据无法修复，因此限制了利用修复后数据对震级的评价.

对于同一个地震的点源机制，使用不同的资料以及不同的反演方法获得的结果可能存在较大差异，本研究中探讨的是由于直接使用限幅数据所产生的误差.以芦山地震为例，在使用相同的速度模型和反演方法条件下，误用限幅数据获得的机制解与参考机制解相比，滑动角误差可达13°，而使用恢复后的地震数据将误差减小到4°以内，可以认为恢复后机制解与参考机制解基本是一致的.由此可见，利用恢复后的地震数据可以获得更为准确的震源机制解，所采用的地震数据恢复方法(Zhang et al., 2016)可在地震应急中利用近场限幅波形记录获取较为准确的点源机制解.

对于有限断层反演震源破裂过程，近场记录以及关键台站如靠近节平面的台站，往往起着关键的作用，如日本熊本地震中，关键强地面运动台站的波形拟合支持了破裂发生在两条断层上的结论(Hao et al., 2017)，玉树地震破裂方向上台站记录的分析揭示了玉树地震发生了超剪切破裂(Wang and Mori, 2012).对于中国大陆地区强地面运动台站分布并不均匀，多数地震区尚未覆盖，而宽频带地震台站分布较为均匀，多数实际震例都有近场的宽频带地震记录，但是存在数据限幅问题.本研究中使用恢复后的区域宽频带数据进行联合反演，无论在破裂区形状还是得到的主要震源参数和利用强地面运动数据进行联合反演获取的结果都十分接近，验证了恢复数据的可用性.因此用该方法恢复近场以及关键台站的波形记录有望在有限断层反演中发挥作用.