Research on bandwidth of ground motion simulated using FK approach
-
摘要: 地震动合成方法的应用前景与其所能表达的有效频带范围紧密相关。本文研究基于频率波数域格林函数的地震动合成方法（FK法），分析频带范围的关键影响要素及处理措施。在介绍计算原理的基础上，概括FK法合成地震动的主要影响要素，分析地壳速度模型对子源地震矩、破裂时间和传播时间的影响。随后，分析格林函数传播宽频带地震波的能力，提出格林函数计算中的频带影响参数的建议值。为分析FK法的震源模型辐射宽频带地震波的能力，重点考察控制子源错动过程的震源时间函数和上升时间。研究认为，FK法具有合成宽频带地震动的能力，这需要采用合理的参数和模型。Abstract: The application prospect of ground motion simulation approach is closely related to the effective bandwidth of synthetics. Hence we studied the simulation approach based on the frequency-wavenumber Green’s function （FK approach）, and analyzed the key factors affecting the bandwidth and the treatment measures. On the basis of the calculation principle introduction, we summarized the main factors influencing the ground motion simulated by FK approach, and analyzed the influences of the crustal velocity model on sub-source seismic moment, rupture time and propagation time. Then, we analyzed the ability of the Green’s function to propagate broadband seismic waves, and put forward the recommended values of parameters affecting bandwidth in the calculation of Green’s function. The source time function and the rise time controlling the sub-source slip process are investigated to analyze the ability of source model in FK approach to radiate broadband seismic waves. The results show that the FK approach has the ability to simulate broadband ground motion, and this requires reasonable parameters and models.
-
引言
在地震破裂过程研究中，破裂模型反演结果的可靠性依赖于断层参数的准确性。研究人员一般通过反演震源机制解获得地震断层参数。在震源机制反演研究方面已积累了大量方法和技术，例如P波初动、P/S振幅比（梁尚鸿等，1984）、面波波形反演（Patton，1980）、体波波形反演（Wallace，Helmberger，1982）、近震全波场反演（Dreger，Helmberger，1993）、CAP （cut and paste）方法（Zhao，Helmberger，1994）和W震相方法（Kanamori，1993）等。对于地震破裂模型研究，则一般采用远场地震波、近场强震地震波、全球定位系统（global positioning system，缩写为GPS）、合成孔径雷达干涉 （interferometric synthetic aperture radar，缩写为InSAR）等数据进行反演。尽管这些反演工作在大地震研究中已经比较成熟，但在中等强度地震的应用中仍然存在一些困难。由于震级相对较小，中等地震的地震波和形变资料信号偏弱，信噪比不高，可能会对破裂模型反演结果造成影响。另一方面，对我国而言，中等强度的地震足以造成破坏，导致人员伤亡和财产损失，同样需要对震源特征开展细致的研究工作。
2017年精河MS6.6地震即是一次典型的中等强度地震。据中国地震台网中心（China Earthquake Networks Center，缩写为CENC）测定，北京时间2017年8月9日上午7时27分52秒，在新疆精河县发生MS6.6地震，震中位置为（44.27°N，82.89°E），震源深度为11 km。截至8月20日23时59分，共记录到324次余震，其中最大震级MS4.7。地震导致36人受伤，未造成人员死亡（阿里木江·亚力昆等，2017；常想德等，2017）。此次地震发生于天山北脉博罗科努山北沿的库松木楔克山前断裂的东段附近。天山南北两侧的构造主要受印度板块与欧亚板块的挤压控制。根据GPS测定结果，天山南北向的缩短速率约为20 mm/a （Abdrakhmatov et al，1996），地壳持续的缩短与增厚使天山不断向南北两侧盆地推覆扩展。库松木楔克山前断裂发育于天山推覆体前缘的逆断层（陈建波等，2007）。根据野外考察结果，本次地震未观测到明显的地表破裂。对本次地震序列重定位的研究（白兰淑等，2017；姜祥华等，2017；刘建明等，2017；徐志国等，2019；翟亮等，2019）认为，库松木楔克山前断裂可能为本次地震的发震构造（图1）。何骁慧等（2020）结合区域地质资料、卫星影像等，也认为本次地震的发震断层为位于库松木楔克山前断裂以北的精河南断层。精河地震附近区域地震活动性较强，自1900年以来，在距本次地震震中200 km范围内，一共发生过6次M≥6.0地震。在精河县境内，最近发生的中等强度地震为2018年10月16日的MS5.4地震（44.19°N，82.53°E）。
![]()
图 1 2017年8月9日精河MS6.6地震及其附近地区的构造示意图Figure 1. Geological background of the Jinghe MS6.6 earthquake on August 9，2017对于精河地震，数据质量可能影响了其破裂过程的反演结果。由于地震震级不大，断层尺度较小，远震地震波可能难以分辨较小空间尺度内的破裂过程；同时，近场强震台站方位覆盖较为有限，且区域沉积层较厚，较强的近地表或场地效应难以在反演中较好地考虑；InSAR数据记录的地表位移较小，造成观测位移的相对误差较大等。已有的研究中，不同破裂模型也存在一定差异。张勇等（2017）使用远场地震波形资料反演得到了本次地震的破裂时空过程，单新建等（2017）、刘传金等（2018）、施贺青等（2019）和Gong等（2019）使用InSAR同震形变数据对本次地震发震断层的滑动分布进行了反演；Zhang等（2020）分别使用远场地震波形反投影和远震、近场强震、InSAR数据联合反演得到了本次地震的破裂特征。在InSAR反演结果中，刘传金等（2018）、Gong等（2019）的破裂区域相对更浅；张勇等（2017）的远震反演结果，其破裂区域距离震源较近；Zhang等（2020）的联合反演结果中破裂则更倾向于向深处扩展。为确证此次地震的更细节的破裂行为，本文发展了一套确定断层面参数的流程，采用多种类型数据联合反演本次地震的破裂模型，以期获得更为详细的、准确的震源破裂细节。
精河MS6.6地震后，一些机构发布了其测定的震源位置和矩张量解（表1）。不同结果中的震中位置较为接近，震源深度位于11—24 km不等。根据本次地震发震构造研究结果（徐志国等，2019；何骁慧等，2020），真实的发震断层面应为近南倾的节面I，然而不同结构给出的南倾节面参数在一定程度上存在差异，无从确定哪个结果更为准确、更适合于破裂过程联合反演。由于不准确的断层参数可能对反演结果造成较大影响，因此本文试图通过密集分布的InSAR同震形变数据来确定更为准确的断层几何模型。
表 1 研究机构发布的2017年精河MS6.6地震震源定位结果和矩张量解Table 1. Epicenter locations and focal mechanism solutions of the 2017 Jinghe MS6.6 earthquake released by different research institutes研究机构 北纬/° 东经/° 震源深度/km 节面Ⅰ 节面Ⅱ 走向/° 倾角/° 滑动角/° 走向/° 倾角/° 滑动角/° CENC 44.270 82.890 11 76 44 80 269 47 99 USGS 44.302 82.832 20 92 60 92 269 30 87 GFZ 44.330 82.840 24 85 47 81 277 44 99 注：GFZ为德国地学中心German Research Centre for Geosciences的缩写. 1. 数据与方法
我们从美国地震学研究联合会（Incorporated Research Institutions for Seismology，缩写为IRIS）数据中心下载了全球地震台网震中距介于40°—90°范围内的远震宽频带垂直向P波记录。为保证良好的方位角覆盖，本文以最小方位角间隔5°对台站进行进一步挑选，最终挑选出27个台站的地震记录（图1左上插图），之后对这些地震记录以0.01—0.4 Hz的频带进行滤波，并将波形重采样为4 sps。
从中国地震局工程力学研究所国家强震动台网中心获取了32个近场强震台站的三分量地震波数据。为了尽可能减弱浅层速度结构对地震波的放大效应，我们使用较低的0.04—0.1 Hz频带对其进行滤波，同样将波形重采样至4 sps。经过多次试反演，最终从32个台站中选出11个波形拟合较好的台站，用于最终的联合反演（图1左下插图）。
利用升降轨道哨兵-1 （Sentinel-1）雷达数据，两个同震干涉像对（Track85和Track63）被处理用于约束断层几何参数和滑动分布（表2）。采用Feng等（2016）发展的自动处理程序gInSAR处理和分析InSAR数据。数据处理中，采用30 m分辨率的SRTM （Shuttle Radar Topography Mission）地形数据用于地形相位的去除和数据配准。为减轻反演迭代中的计算压力，采用四叉树算法对原始干涉图进行降分辨率重采样，将Track85和Track63的采样点分别降为3 285个和4 271个（表2）。
表 2 本文所用雷达数据信息Table 2. Information of the radar data used in this study轨道信息 轨道方向 时间信息 空间基线/m 时间基线/d 采样点个数 震前图 震后图 T85 升轨 2017-08-08 2017-08-14 − 92.2 37 3 285 T63 降轨 2017-08-07 2017-08-13 −119.6 37 4 271 对于地震波记录，基于CRUST1.0模型（Laske et al，2013），采用Wang （1999）的Qseis程序计算格林函数。对于InSAR记录，采用基于弹性位错理论的Okada （1985）方法计算地表形变。
2. 断层模型
有限断层破裂过程反演中，断层面为通过震源（点源）的平面断层。需将其沿走向和倾向划分为若干个正方形的子断层，震源位于其中一个子断层的中心。考虑到地震发生在隐伏断层上，断层几何参数不能很好地确定。为获得合理的断层参数，基于InSAR同震形变数据设计了一套确定断层面走向、倾角和震源深度的格点搜索流程。基于InSAR形变数据，采用不同的断层几何参数和位置参数进行静态滑动分布反演，根据残差最小的反演结果所对应的走向、倾角和震源位置，构建更合理的断层几何模型，并得到震源在断层面上的相对位置。
在以上搜索流程中，我们将对断层面位置参数的搜索转化为对震源深度的搜索。首先，设置一个较大且能较好地覆盖潜在滑动区域的断层面；其次，将断层面的上边界固定于地表，断层面沿倾向方向的移动等同于改变震源所在子断层在倾向方向上的序数，这样震源深度即可由震源到断层面上边界的距离乘以断层面倾角的正弦值得到。通过这种做法，只需对震源所在子断层在倾向方向上的序数进行搜索，即实现了对震源深度的搜索，同时也完成了对断层面沿倾向方向位置的搜索。
我们将三维搜索空间中的任意一个格点（走向、倾角、序数）记为FP （fault-plane parameters）。理想情况下，存在一个全局最优点，具有最小残差，三维搜索空间中的其它点的残差均随着最小残差点的距离单调递增，即最小残差点是一个稳定的点。
InSAR是一种相对测量技术，即任意干涉像对中理论上只能得到相对形变。通常如精河MS6.6地震这样的小形变事件，可以假定远场形变为零，进而得到近场绝对形变量。但由于大气和电离层的影响，得到的形变场中仍然可能含有一定的零漂和线漂成分。考虑一般情况，本文尝试从反演流程中自动估计数据偏移，实现在反演中进行去漂移处理。本文在Zhang等（2012）的InSAR同震滑动反演方法的基础上，采用线性拟合的方式去除数据中的漂移。InSAR观测数据可以拆分为去漂移后的数据和漂移项，即
$$ {{\boldsymbol{Q}}_{{\rm{ob}}}}{\text{＝}}{{\boldsymbol{Q}}_0} {\text{＋}} {{\boldsymbol{Q}}_{{\rm{drift}}}}{\text{，}} $$ (1) 式中，Qob为观测数据，Q0为无漂移数据，Qdrift为漂移大小。Q0可以表示为
$$ {{{\boldsymbol{Q}}_{{0}}}} {\text{＝}}{\boldsymbol{B}} {\boldsymbol{f}} {\text{，}} $$ (2) 式中，B为由断层面上各子断层单位滑动形成的地表形变矩阵，f为包含子断层滑动大小的未知参数。假设漂移是线性的，则Qdrift可以表示为
$$ {{\boldsymbol{Q}}_{{\rm{drift}}}}{\text{＝}}{{\boldsymbol{s}}_0} {\text{＋}} {\boldsymbol{X}}{{\boldsymbol{s}}_x} {\text{＋}} {\boldsymbol{Y}}{{\boldsymbol{s}}_y}{\text{，}} $$ (3) 式中：X和Y分别为平面坐标系中InSAR数据覆盖区域内各点的坐标值；s0，sx和sy分别为拟合参数。则考虑漂移的InSAR同震滑动反演方程即可列为
$$ {{{\boldsymbol{Q}}_{{\rm{ob}}}}}{\text{＝}}\left[\!\!\!\! {\begin{array}{*{20}{c}} {\boldsymbol{B}}&{\boldsymbol{I}}&{\begin{array}{*{20}{c}} \!\!\!\!{\boldsymbol{X}}&{\boldsymbol{Y}} \end{array}} \end{array}} \!\!\!\!\!\!\!\!\right]\left[\!\!\!\!\!\!\!\! {\begin{array}{*{20}{c}} {\boldsymbol{f}}\\ {{{\boldsymbol{s}}_0}}\\ {\begin{array}{*{20}{c}} {{{\boldsymbol{s}}_{{x}}}}\\ {{{\boldsymbol{s}}_{{y}}}} \end{array}} \end{array}} \!\!\!\!\!\!\!\!\!\!\right]{\text{，}} $$ 式中，I为元素全为1的列向量。基于该方程，我们可以反演出InSAR数据的漂移项，最终得到去漂移后的InSAR数据。
在去漂移的InSAR反演中，需要给定断层面参数（initial fault-plane parameters，缩写为IFP）。通过给定不同的IFP并对InSAR数据进行去漂移处理后看出，不同的IFP得到的无漂移InSAR数据Q0存在明显的区别。可能的原因是，精河MS6.6地震破裂产生的地表形变量较小，接近InSAR数据的固有误差水平（约0.1 m），即InSAR数据的相对误差较大，使得IFP较小的变动也会导致去漂移结果的明显变化。为解决此问题，我们尝试了多组IFP，得到多组去漂移的InSAR数据，然后对每一组处理后的InSAR数据进行三维格点搜索，得到多个FP值，记为SFP （searched fault-plane parameters），再根据得到的SFP点的分布确定合适的、适用于最终联合反演的SFP。图2给出了基于InSAR数据滑动分布反演的三维格点搜索流程。
![]()
图 2 基于InSAR滑动分布反演的格点搜索流程示意图Figure 2. Flow chart of the grid search based on the InSAR slip inversions反演结果显示，精河地震在断层面上约21 km×14 km范围内产生滑动，为避免边缘效应，我们将断层面尺度适当拓展为30 km×30 km，并沿走向方向和倾向方向分别划分为15个子断层，即每个子断层大小为2 km×2 km。断层面滑动角在90°±45°之间变化。经过若干次反演尝试，我们将震源固定在沿走向方向的第12个子断层上，并采用中国地震台网中心的定位结果（44.27° N，82.89° E）作为震中位置。将准备的256个IFP作为输入：走向取值分别为75°，80°，85°，90°，95°，100°，105°和110°；倾角取值分别为37°，40°，42°，45°，47°，50°，52°和55°；震源在断层倾向方向上的序数分别为8，9，10，11。采用这256组初始断层面参数对原始InSAR数据分别进行去漂移反演，得到256组去漂移的InSAR数据，随后对256组去漂移的InSAR数据进行基于InSAR数据滑动分布反演的三维格点搜索：各组数据IFP走向搜索范围为±16°，间隔1°；倾角搜索范围为±16°，间隔1°；序数搜索范围为±4，间隔为1。搜索完成后，得到256个SFP点的分布。从所有输出结果为该SFP点的IFP点中寻找残差最小的IFP点，并采用该IFP点去漂移得到的InSAR数据用于联合反演。
图3给出了搜索得到的256个SFP点的个数和密度分布，其中每一点的密度计算方法为该点SFP数量加上边与之相接的四个点SFP数量之和的四分之一，再加上顶点与之相接的四个点SFP数量之和的八分之一。从三维搜索空间的三个侧面均能看到，SFP点存在一定的聚集。通过计算得到的256个SFP点在走向、倾角和序数的平均值为[95，47，10]，并将其作为搜索得到的断层面参数。为考察该点是否稳定，我们比较了该点附近256个SFP点的平均残差。结果显示，从三个侧面来看，[95，47，10]在搜索空间中均为其附近的最小残差点，结果具有稳定性。然后，我们找到输出结果为SFP＝[95，47，10]的最小残差的IFP点，即IFP＝[100，40，8]，并将该IFP点的去漂移InSAR数据用于接下来的联合反演。
![]()
图 3 搜索得到的256个SFP点的分布图（a）—（c）分别为倾角与走向、震源所在子断层沿倾向方向序数（简称序数）与走向、序数与倾角剖面所显示的SPF三维空间分布，圆圈表示SFP点，其尺寸对应残差大小；图（d）—（f）分别为与图（a）—（c）对应的SFP分布密度；图（g）—（i）分别为已选取的SFP＝[95，47，10]为中心，搜索空间中的走向截面、倾角截面和序数截面的平均残差Figure 3. Distribution of the searched 256 SFPsFigs.（a）− （c） show the SFP distribution in the cross sections of strike-dip，ordinal number of subfault along dip the epicenter located in （“ordinal”for short）-strike，and ordinal-dip，respectively，where the circles indicate the SFPs. The size of the circle denotesthe residuals，and circle color represents the number of SFPs. Figs.（d）−（f） show the density distribution of SFP correspon-ding to Figs.（a）− （c），respectively. Figs.（g）−（i） show the average residuals in the search space （selected SPF＝[95，47，10] centered） in the cross section of strike，dip and ordinal，respectively3. 破裂过程联合反演
基于搜索流程给出的结果，我们进一步明确联合反演的断层模型参数：断层面走向为95°，倾角为47°，为一个30 km×30 km的矩形，并将其划分为225个2 km×2 km的子断层。震源位于沿走向方向的第12个子断层也是沿倾向方向的第10个子断层的中心点。反演中滑动角允许在90°±45°之间变化。最大破裂速度限定在2.8 km/s以内，允许每个子断层上的最长破裂持续时间为4 s。
我们采用远震数据、近场强震数据和InSAR数据，基于Zhang等（2012）的方法进行震源破裂过程的联合反演。反演中另行添加了空间光滑约束、时间光滑约束（Horikawa，2001；Yagi et al，2004）和最小标量地震矩约束（Hartzell，Iida，1990），采用非负共轭梯度法（Ward，Barrientos，1986）求解反演方程组。经尝试，给定强震数据3倍于远震数据的权重，并赋予InSAR数据和地震波数据同等权重。同时，为获得更好的波形拟合结果，允许地震波形在时间上进行平移：当首次联合反演完成后，对于远震波形，在时间域上对每个台站的观测波形进行不超过5 s的平移，使得观测波形与合成波形的相关系数达到最大；对于强震波形，将同一台站的三分量波形同时进行不超过10 s的平移，使得三分量的观测波形与合成波形的平均相关系数最大。平移完成后，再次联合反演，确定最终的反演结果。
图4展示了联合反演得到的滑动分布和震源时间函数。反演结果显示，精河MS6.6地震为一次单侧破裂事件，滑动集中分布于震源以西5—15 km，从地表沿断层面向下15—25 km，考虑倾角因素，对应为11—18 km的深度范围。最大滑动量约为0.8 m，所处位置在震源以西8 km，深度与震源相近。滑动整体没有出露地表。断层面上的平均滑动角为106°。震源时间函数显示本次地震只有一次主事件，破裂整体持续约9 s：从破裂开始地震矩快速释放，第2.5 s时达到顶峰，随后平缓下行至第9 s。从破裂过程的快照，能够清楚地看到破裂向西传播的过程。反演得到的标量地震矩为3.6×1018 N·m，对应矩震级为MW6.3，与中国地震台网中心、USGS和德国地球科学研究中心（German Research Center for Geosciences，缩写为GFZ）的结果一致。图5显示了地震波形和InSAR数据的拟合情况。
![]()
图 4 基于远震资料、近场强震资料和InSAR资料的地震破裂过程联合反演结果（a） 滑动分布和震源时间函数；（b） 子断层上的子震源时间函数；（c） 每2 s积累的滑动量分布快照Figure 4. Joint inversion results of rupture process based on teleseismic，near-field strong-motion and InSAR data（a） Slip distribution and source time function；（b） Subfault source time functions；（c） Snapshots ofthe slip accumulated at each 2-second interval![]()
图 5 联合反演的观测与合成资料比较（a） 观测（黑线）与合成（红线）远震地震波；（b） 观测（黑线）与合成（红线）近场强震地震波；（c） 观测与合成InSAR数据Figure 5. Comparison of the observed and synthetic data of the joint inversion（a） Observed （in black） and synthetic （in red） teleseismic data； （b） Observed （in black） and synthetic （in red） strong-motion data；（c） Observed and synthetic InSAR data图6给出了地震破裂过程时空图像。主要考虑矩心的破裂速度，由于反演结果的不确定性，为了去除破裂滑动中的非主要信息，在计算每1 s内的矩心位置时，我们将1 s内滑动量小于整个破裂过程最大累积滑动量15%的子断层剔除。同时，考虑到反演结果的不确定性所导致的矩心瞬时速度的不稳定，采用了平均破裂速度（自破裂发生起矩心的累积位移除以已经发生的破裂时间），并计算以每1 s内地震矩大小为权重的加权平均破裂速度。平均破裂速度计算结果显示，破裂速度从第2 s的1.0 km/s缓慢增加至第6 s的2.1 km/s，后缓慢降至第8 s的1.9 km/s；从加权平均破裂速度来看，破裂速度从第2 s的1.0 km/s增加至第5 s的2.6 km/s，后降至第8 s的2.4 km/s。综合来看，本次地震的破裂速度介于2.1—2.6 km/s。
![]()
图 6 矩心破裂传播速度（a） 矩心的位置、时间和平均移动速度，图中时间为破裂开始后经过的时间 ；（b） 矩心的位置、时间和地震矩加权平均移动速度；（c） 破裂传播速度-时间曲线Figure 6. Rupture velocity of the moment centroid（a） The location，time and average velocity of the centroid，the time represents the duration after the rupture；（b） The location，time and moment-weighted average velocity of the centroid；（c） Velocity-time curves4. 讨论
对于精河MS6.6地震这样的中等强度地震，破裂模型反演结果可能会受到数据的影响。因此，采用大地测量数据确定断层面参数，并综合远震地震数据、近场强震数据和大地测量数据进行联合反演，有助于得到更好的破裂模型。
本文也尝试采用三维格点搜索流程对真实的发震断层进行判别（图7）。在对北倾节面的FP进行搜索后，同样得到了256个SFP点的分布。北倾节面对应的SFP相对残差范围为0.081—0.136，平均值为0.096；而南倾节面SFP的相对残差范围为0.080—0.148，平均值为0.098。在北倾节面对应的SFP中选取[261，38，8]作为联合反演使用的参数，该SFP点对应的最小残差为0.094，同时南倾节面中选取的SFP对应的最小残差为0.090。从残差分布范围来看，北倾节面与南倾节面重合度较高；北倾节面平均残差略低于南倾节面，但二者差别较小。综合来看，我们采用的搜索流程难以对发震断层进行有效地判别，其原因是此次地震破裂位置较深，造成的有效地表形变较小，导致InSAR数据的信噪比较低。基于距离破裂区域较远且具有一定误差水平的数据，在区分断层的倾向上存在困难。
![]()
图 7 南倾节面与北倾节面对应的SFP残差分布Figure 7. Residuals of the SFPs for south-dipping and north-dipping nodal planes需要说明的是，一些因素可能会影响联合反演得到的破裂模型。首先，区域沉积层对近场强震台站记录的地震波形有较大影响，导致强震波形拟合结果相较于远震拟合更差。虽然本文采用的CRUST1.0模型包含区域沉积层速度结构信息，但由于该模型的最小分辨尺度相对本次地震破裂尺度大得多，其包含的沉积层速度信息可能与实际速度结构存在一定差别，因而在计算理论格林函数时，我们对区域沉积层速度结构的考虑有可能不够全面。若能获得震中附近区域更为精确的速度结构信息，则有可能提高强震资料的拟合。其次，强震台站大部分位于震中南侧，方位覆盖较差，这样的台站布局可能对反演结果产生一定影响。但从反演结果来看，虽然方位角分布不均匀，但强震数据反演得到的破裂模型与其它数据的反演结果具有一致性；同时，强震数据有助于更好地约束地震破裂位置和位错大小，在时空上的分辨能力比远震数据更强。因此，即使强震台站分布不均匀，强震数据在破裂过程反演中仍然具有一定的优势。最后，震源定位的结果直接影响破裂模型中断层面上滑动区域的位置和破裂速度的计算结果。因此，若能采用其它方法对震源进行更精确的重定位，则有可能提高反演得到破裂模型的准确程度。
为考察在联合反演中不同类型的数据所起的作用，我们分别采用远震数据、强震数据和InSAR数据进行了单独反演。反演结果如图8所示。其中，强震数据反演、InSAR反演和联合反演显示出较为一致的破裂方向性，即破裂从震源向西传播了10—15 km，而远震反演结果中的滑动集中在震源附近，未显示明显的破裂方向性，这与张勇等（2017）的结果相似。这表明，对于中等强度地震产生的破裂，远震P波（本文中为震中距40°以上）可能难以有效地分辨破裂方向性，而距离较近（本文中为震中距3°以内）的强震台站和InSAR数据能够较好地分辨。由于在反演中，InSAR数据对破裂区域的绝对位置的约束力较好，因此我们倾向于认为本次地震破裂向西传播。虽然在本震例中远震数据对于破裂方向性难以分辨，但由于远震地震波穿过地幔，有助于约束震级，另外远震地震波对于断层的不规则性和弯曲不敏感，能够更好地约束断层宏观破裂行为，因此在联合反演中仍起到一定作用。联合反演中，经过多次尝试，我们将远震数据与强震数据的权重比设置为1 ∶ 3。
![]()
图 8 不同类型数据反演得到的滑动分布和震源时间函数（a） 远震、近场强震和InSAR数据联合反演；（b） 远震数据反演；（c） InSAR数据反演；（d） 近场强震数据反演Figure 8. Slip distributions and source time functions of inversions with different data（a） Joint inversion of teleseismic，strong-motion and InSAR data；（b） Inversion of teleseismic data；（c） Inversion of InSAR data；（d） Inversion of strong-motion data我们将联合反演结果与张勇等（2017）的远震快速反演结果、Gong等（2019）的InSAR反演结果、Zhang等（2020）的联合反演结果进行比较。从反演得到的矩震级看，本文结果、张勇等（2017）、Gong等（2019）、Zhang等（2020）依次为MW6.3，MW6.3，MW6.2，MW6.3，结果基本一致；最大滑动量分别为0.8 m，1.1 m，0.3 m和0.6 m，其中远震反演结果最大，InSAR反演结果最小，联合反演结果居中，与我们进行的不同类型数据反演结果一致；从破裂方向性上看，除张勇等（2017）的远震反演结果外，其余反演结果均显示出较为明显的西向破裂，显示远震数据对本次地震破裂方向性的分辨能力较差。此外，Zhang等（2020）的联合反演结果显示破裂向深部扩展，而本文联合反演和Gong等（2019）显示破裂向浅部扩展，这可能是Zhang等（2020）基于北倾节面反演所致，但各反演结果中破裂未出露地表的主要特征是一致的。
在联合反演中，我们在初次反演后，将观测地震波形在时间上进行平移，使得合成波形与相关波形之间具有最大的相关系数，并在此基础上再次进行反演得到最终结果。这一操作的目的是只拟合波形信息，而不考虑到时误差。这有助于减小理论到时的不准确性对反演结果的干扰，也是震源波形反演中常见的操作（Yagi et al，2004）。对于远震数据，由于地震波传播距离较远，到时通常存在最大1—3 s的不确定性，因此在远震数据反演中允许波形时间平移是较为合理的。对于强震数据，由于精河地震发生区域沉积层较厚，浅层速度结构对强震波形有较大的影响，不但影响到时的准确性，还导致强震数据质量偏低，这也是我们在32个强震台站中最终只挑选出11个台站进行反演的原因之一。本文中，我们允许强震波形的时间平移，即选择放弃使用到时信息，以提高波形的拟合程度。同时注意到，在联合反演中InSAR数据对破裂绝对位置的约束较好，部分补偿了强震数据到时信息缺失导致的对破裂绝对位置分辨能力下降的问题。因此，允许强震波形时间平移，能尽可能地降低到时信息不准确造成的负面影响。
精河MS6.6地震主震破裂分布与余震分布可能存在一定的互补关系。图1显示了主震滑动分布和余震分布在地表上的投影。我们共选取2017年8月14日16:00之前中国地震台网中心记录到的209次余震事件，其位置已由白兰淑等（2017）采用双差定位法（Waldhauser，Ellsworth，2000）进行重定位确定。从地表上看，大部分余震集中在破裂区域的西边，小部分余震集中在破裂区域的西南角。余震分布和主震破裂分布的互补性说明主震对余震具有很强的触发作用。主震破裂后，在破裂区域周围引起较强的应力扰动和库仑应力增加，可能是震后余震密集发生的主要原因（Harris，1998；Kilb et al，2002）。
5. 结论
本文中，我们利用InSAR同震形变数据，设计了一套基于滑动分布反演的格点搜索流程，得到了2017年8月9日新疆精河MS6.6地震断层面的走向、倾角和震源深度；在此基础上，利用远震资料、近场强震资料和InSAR同震形变资料，联合反演得到了此次地震的震源破裂时空过程。根据搜索和反演结果，本次地震断层面走向为95°，倾角为47°，震源深度为14 km。整个地震主要以单侧破裂为主，破裂过程持续约9 s，破裂朝西以2.1—2.6 km/s的速度传播，最终在震源以西5—15 km、沿倾向15—25 km范围形成了一个主要的滑动区域，最大滑动量约0.8 m，平均滑动角为106°。整个破裂过程释放的标量地震矩为3.6×1018 N·m，对应矩震级为MW6.3。由于此次地震的破裂深度主要集中在10 km以下，未来可能需要关注震源区浅部0—10 km处发生潜在地震的可能。
本文研究所用的远震数据下载自美国地震学研究联合会IRIS，强震资料由中国地震局工程力学研究所国家强震动台网中心提供，InSAR资料来自于欧洲空间局哨兵-1雷达数据。作者在此对以上机构一并表示感谢。
-
![]()
图 2 不同破裂速度计算的破裂时间。星号表示破裂起始点
（a） 同一破裂速度；（b） 深度相关破裂速度；（c） 深度相关破裂速度的分布；（d） 图（b）相对图（a）的差值
Figure 2. Rupture time calculated by different rupture velocity vr
（a） The same vr ；（b） Depth-related vr；（c） Distribution of depth-related vr； （d） Difference between Fig.（b） and Fig.（a）
![]()
图 3 不同剪切波速计算的S波的传播时间
（a） 同一剪切波速；（b） 深度相关剪切波速；（c） 图（b） 相对图 （a） 的差值
Figure 3. Propagation time of S wave calculated by different
$v_{\rm{S}} $ （a） The same $v_{\rm{S}} $；（b） Depth-related $v_{\rm{S}} $；（c） Difference between Fig.（b） and Fig.（a）
![]()
图 4 上升时间为0 s （a）和1.5 s （b） 时奈奎斯特频率对合成地震动加速度幅值谱FAS的影响
Figure 4. Influence of Nyquist frequency on FAS of synthetic ground motion with rise time of 0 s （a） and 1.5 s （b）
![]()
图 5 最大波数（a）和波数间隔（b）对合成地震动加速度幅值谱FAS和位移时程的影响
Figure 5. Influence of the maximum wavenumber （a） and the wavenumber interval （b） on FAS and displacement time history of synthetic ground motion
![]()
图 6 品质因子对合成地震动加速度幅值谱FAS的影响
Figure 6. Influence of quality factor on FAS of synthetic ground motion
![]()
图 7 五种震源时间函数的错动速率（a）及其加速度幅值谱FAS （b）
Figure 7. Slip rate （a） and its FAS （b） for five source time functions
![]()
图 8 震源时间函数对合成地震动加速度幅值谱FAS的影响
Figure 8. Influence of source time function on FAS of synthetic ground motion
![]()
图 9 上升时间（a）和破裂速度（b）对合成地震动加速度幅值谱FAS的影响
Figure 9. Influence of rise time （a） and rupture velocity （b） on FAS of synthetic ground motion
表 1 本文采用的地壳速度结构模型
Table 1 Crustal velocity structure model adopted in this study
厚度/km vP/（km·s−1） vS/（km·s−1） ρ/（g·cm−3） QP QS 3 4.88 2.86 2.55 300 150 5 5.80 3.40 2.70 400 200 14 6.04 3.55 2.75 400 200 21 6.82 3.98 2.90 600 300 25 7.61 4.45 3.10 800 400 - 8.08 4.47 3.38 1 000 500 
下载: 导出CSV
表 2 不同dt时FK法合成地震动的计算量
Table 2 Calculation cost of ground motion simulation using FK approach for different dt
算例 子源数 dt/s 地表某点地震动的
计算耗时/s地震动场（网格2×2 km2） 耗时/h 存储空间/GB 范围/km2 本文算例
MW6.516×8 0.02 1.4 0.7 0.23 100×70 16×8 0.01 2.8 1.4 0.46 100×70 汶川地震
MW7.964×32 0.02 22.4 126.3 2.56 400×200 64×32 0.01 44.8 252.6 5.12 400×200
下载: 导出CSV
表 3 不同dt，dk和kmax时FK法格林函数的计算量
Table 3 Calculation cost of Green’s function in FK approach with different dt, dk and kmax
dt/s dk kmax 计算耗时（源深6 km）/s 150个地表点60组（源深差0.5 km） 30 km处一点 150个地表点 耗时/h 存储空间/GB 0.02 0.1 15 8 239 3.98 3.09 0.01 0.1 15 34 954 15.90 6.18 0.02 0.05 15 17 483 8.05 3.09 0.02 0.3 15 3 81 1.35 3.09 0.02 0.1 5 8 236 3.93 3.09 0.02 0.1 20 8 240 4.00 3.09 注：150个地表点的震中距设为1，2，···，150 km
下载: 导出CSV
-
曹泽林,陶夏新. 2018. 基于频率波数域格林函数的宽频带地震动合成方法综述[J]. 地震工程与工程振动,38(5):35–42. Cao Z L,Tao X X. 2018. Review on broadband ground motion simulation based on frequency-wavenumber Green’s function[J]. Earthquake Engineering and Engineering Dynamics,38(5):35–42 (in Chinese).
姜伟,陶夏新,赵凯. 2017. 基于NGA数据的震源模型全局参数定标律的统计[J]. 地震工程学报,39(2):221–226. doi: 10.3969/j.issn.1000-0844.2017.02.0221 Jiang W,Tao X X,Zhao K. 2017. Scaling laws of the global parameters of source models from NGA data[J]. China Earthquake Engineering Journal,39(2):221–226 (in Chinese).
孙晓丹,陶夏新. 2012. 宽频带地震动混合模拟方法综述[J]. 地震学报,34(4):571–577. doi: 10.3969/j.issn.0253-3782.2012.04.013 Sun X D,Tao X X. 2012. Hybrid simulation of broadband ground motion:Overview[J]. Acta Seismologica Sinica,34(4):571–577 (in Chinese).
Aki K, Richards P G. 2002. Quantitative Seismology[M]. 2nd ed. Sausalito, California: University Science Books: 37−62.
Beresnev I A. 2002. Source parameters observable from the corner frequency of earthquake spectra[J]. Bull Seismol Soc Am,92(5):2047–2048. doi: 10.1785/0120010266
Brune J N. 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes[J]. J Geophys Res,75(26):4997–5009. doi: 10.1029/JB075i026p04997
Cao Z L,Tao X X,Tao Z R,Tang A P. 2019. Kinematic source modeling for the synthesis of broadband ground motion using the f-k approach[J]. Bull Seismol Soc Am,109(5):1738–1757. doi: 10.1785/0120180294
Day S M. 1982. Three-dimensional simulation of spontaneous rupture:The effect of nonuniform prestress[J]. Bull Seismol Soc Am,72(6A):1881–1902.
Dreger D,Tinti E,Cirella A. 2007. Slip velocity function parameterization for broadband ground motion simulation[J]. Seismol Res Lett,78(2):308.
Douglas J,Aochi H. 2008. A survey of techniques for predicting earthquake ground motions for engineering purposes[J]. Surv Geophys,29(3):187–220. doi: 10.1007/s10712-008-9046-y
Fortuño C,de la Llera J C,Wicks C W,Abell J A. 2014. Synthetic hybrid broadband seismograms based on InSAR coseismic displacements[J]. Bull Seismol Soc Am,104(6):2735–2754. doi: 10.1785/0120130293
Frankel A. 2009. A constant stress-drop model for producing broadband synthetic seismograms:Comparison with the Next Generation Attenuation relations[J]. Bull Seismol Soc Am,99(2A):664–680. doi: 10.1785/0120080079
Graves R,Pitarka A. 2010. Broadband ground-motion simulation using a hybrid approach[J]. Bull Seismol Soc Am,100(5A):2095–2123. doi: 10.1785/0120100057
Hao J L,Ji C,Wang W M,Yao Z X. 2013. Rupture history of the 2013 MW6.6 Lushan earthquake constrained with local strong motion and teleseismic body and surface waves[J]. Geophys Res Lett,40(20):5371–5376. doi: 10.1002/2013GL056876
Hartzell S. 1978. Earthquake aftershocks as Green’s functions[J]. Geophys Res Lett,5(1):1–4. doi: 10.1029/GL005i001p00001
Hartzell S,Guatteri M,Mai P M,Liu P C,Fisk M. 2005. Calculation of broadband time histories of ground motion,Part Ⅱ:Kinematic and dynamic modeling using theoretical Green’s functions and comparison with the 1994 Northridge earthquake[J]. Bull Seismol Soc Am,95(2):614–645. doi: 10.1785/0120040136
Hartzell S,Liu P,Mendoza C,Ji C,Larson K M. 2007. Stability and uncertainty of finite-fault slip inversions:Application to the 2004 Parkfield,California,earthquake[J]. Bull Seismol Soc Am,97(6):1911–1934. doi: 10.1785/0120070080
Hartzell S,Frankel A,Liu P,Zeng Y,Rahman S. 2011. Model and parametric uncertainty in source-based kinematic models of earthquake ground motion[J]. Bull Seismol Soc Am,101(5):2431–2452. doi: 10.1785/0120110028
Irikura K,Miyake H. 2011. Recipe for predicting strong ground motion from crustal earthquake scenarios[J]. Pure Appl Geophys,168(1/2):85–104.
Kennett B L N. 2009. Seismic Wave Propagation in Stratified Media[M]. 2nd ed. Canberra: ANU Press: 59–100.
Kieling K,Wang R,Hainzl S. 2014. Broadband ground-motion simulation using energy-constrained rise-time scaling[J]. Bull Seismol Soc Am,104(6):2683–2697. doi: 10.1785/0120140063
Konno K,Ohmachi T. 1998. Ground-motion characteristics estimated from spectral ratio between horizontal and vertical components of microtremor[J]. Bull Seismol Soc Am,88(1):228–241. doi: 10.1785/BSSA0880010228
Liu P,Archuleta R J,Hartzell S H. 2006. Prediction of broadband ground-motion time histories:Hybrid low/high-frequency method with correlated random source parameters[J]. Bull Seismol Soc Am,96(6):2118–2130. doi: 10.1785/0120060036
Motazedian D,Atkinson G M. 2005. Stochastic finite-fault modeling based on a dynamic corner frequency[J]. Bull Seismol Soc Am,95(3):995–1010.
Sun X,Hartzell S,Rezaeian S. 2015. Ground-motion simulation for the 23 August 2011,Mineral,Virginia,earthquake using physics-based and stochastic broadband methods[J]. Bull Seismol Soc Am,105(5):2641–2661. doi: 10.1785/0120140311
Tinti E,Fukuyama E,Piatanesi A,Cocco M. 2005. A kinematic source-time function compatible with earthquake dynamics[J]. Bull Seismol Soc Am,95(4):1211–1223. doi: 10.1785/0120040177
Wang R. 1999. A simple orthonormalization method for stable and efficient computation of Green’s functions[J]. Bull Seismol Soc Am,89(3):733–741.
Zhu L,Rivera L A. 2002. A note on the dynamic and static displacements from a point source in multilayered media[J]. Geophys J Int,148(3):619–627. doi: 10.1046/j.1365-246X.2002.01610.x
-
期刊类型引用(1)
1. 郑兴群，陶正如，白凯. 面向地震动估计需求的区域传播介质参数. 地震地质. 2024(05): 1091-1105 . 
百度学术
其他类型引用(4)
计量
- 文章访问数: 394
- HTML全文浏览量: 226
- PDF下载量: 92
- 被引次数: 5
