引言

重复地震事件组作为一种特殊的地震序列，是同一断层或者断层的一部分在一定时间内的重复破裂。同一重复地震事件序列内的不同地震事件具有高度重叠的震源破裂区及相近的震级和震源机制解，因此，同一地震台站记录到的波形高度相似。Chen和Lapusta （2009）认为重复地震是大尺度蠕变区内凹凸体的破裂，通过在速度-状态摩擦本构定律基础上构建断层数值模型，可成功地解释重复地震的周期与地震矩之间的关系。利用重复地震检测地球内部介质随时间的变化和估算震源滑动速率已被广泛应用。通过分析重复地震事件序列的波形特征变化，可以得到地球内部介质的物理性质随时间的变化，进而分析该变化与地震孕震环境的关系，表明重复地震事件具有重要的实际应用价值（Schaff，Beroza，2004；Peng，Ben-Zion，2006；Rubinstein et al，2007； Pacheco et al，2017 ）重复地震还为直接获取断层在地表以下不同深度的滑动速率提供了新的途径（Nadeau，McEvilly，1999；Igarashi et al，2003 ；周龙泉等，2007；Kato et al，2012 ；李乐等，2013，2015）。但判定重复地震序列中的不同地震事件难度较大，依据重复地震的定义，是否具有高度重叠的震源破裂区域是判定重复地震的根本准则，实际操作中可通过高精度地震定位和波形记录互相关两种方法进行判定。然而，依据现有实际地震台网密度和模型结构，地震定位方法的精度还不足以有效地判定重复地震事件；除此之外，重复地震序列还要求事件对具有高度相似的震源机制解，而地震定位并不能约束震源机制解，因此，波形互相关方法是实际应用中更为有效的重复地震判定方法。

由于地球内部普遍存在尺度不一的非均匀体，其散射作用能够改造地震波传播的波场。因此，即使两次地震事件具有相同的震源机制解，若震源位置存在偏差，对应的地震波形记录的相似性也会受到非均匀体散射的影响，震源间距越大，波形的相似度则越低。据此，Geller和Mueller （1980）提出了1/4优势波长准则，即两次地震事件的波形若高度相似，则事件的震源间距应小于1/4优势波长。Baisch等（2008）通过有限差分数值模拟二维非均匀介质内地震波的传播，进而研究波形相似性与震源距离之间的关系，得出当截取波形窗口长度为2.8倍P波与S波到时差时，互相关系数阀值可定为0.95，并最终支持1/4优势波长准则。这一结果定量地研究了地震对空间距离与互相关系数变化的关系，但其研究时间和空间尺度均为二维人工勘探尺度，且未进一步探讨非均匀体性质的影响。鉴于此，本文拟结合宁夏区域地震台网分布和实际的地震分布情况，通过谱元法模拟三维结构下地震波的传播，分析在不同台站上观测到的天然地震对的波形互相关系数的分布情况，并进一步研究其波型的互相关系数与非均匀体性质的关系，为重复地震的检测提供参考依据。

1.   宁夏区域台网分布与重复地震检测

宁夏位于我国强震密集的南北地震带北段，地处青藏地块、阿拉善地块以及鄂尔多斯地块交会区。2008年实行数字化观测以来，宁夏及邻区的数字地震台网布局趋于合理化，随着台站密度不断增大，对区域内中小地震的监测能力逐渐提高，2008—2016年区域内已监测到M_L≥0.1地震事件3 868次，具体分布如图1所示，其中M_L≥4.0地震12次。

图  1  宁夏区域地震台网及2008—2016年地震震中分布

Figure  1.  Location of Ningxia regional seismic network and epicentral distribution from 2008 to 2016

下载: 全尺寸图片 幻灯片

宁夏区域地震台网平均台间距为45 km（图1），可较好地监测M_L≥0.1地震，但是将其分布密度用于判定重复地震，其精度无论是距离绝对定位还是相对定位均相差甚远。假设重复地震震级为M_L2.0，其破裂空间尺度仅有几十至几百米，而宁夏区域地震台网进行绝对定位的误差大于1 km （金春华等，2015a，b），即使利用相对定位法，其精度也不足以保证重复地震对判定的可靠性。由于大多数断层的滑动速率均小于10 mm/a，因此，在几年或十几年的时间尺度内，可利用的重复地震主要集中在M_L1.0—3.0。

宁夏区域地震台网的密度在全国具有较好的代表性，依此台网密度分布，利用地震定位方法并不能有效地检测绝大多数M_L2.0左右的重复地震，但通过波形互相关法却可在该情况下有效地检测重复地震。如图2所示，为使用波形互相关检测到的一组波形相似事件组，该事件组由6次M_L1.4地震组成，相邻地震事件之间的时间间隔从十几天至一年多不等。这一事件组内的垂向和南北向速度分量的波形互相关系数均较高，且除事件2010-08-28外，南北向分量间的互相关系数均比垂直向分量的更大。若使用绝对地震定位方法，由于存在较大误差，约66次地震事件落在图2对应的同一位置。然而，除图2中的6次事件，其它均为位置相近的地震事件，并非重复地震。

图  2  5—10 Hz带通滤波后的地震波南北向（a）和垂向（b）速度分量的波形相似序列

图中各波形右上角数字为地震事件与参考事件（2010-06-13）的波形互相关系数

Figure  2.  High similarity waveforms with the north-south components （a） and the vertical components （b） after 5−10 Hz band-pass filtering

The numbers in top right corners give the cross-correlation coefficients between the associated seismogram and template seismogram （2010-06-13）

下载: 全尺寸图片 幻灯片

波形互相关方法比地震定位方法更适用于判定重复地震，较高的互相关系数是判断重复地震的关键指标，但是否为重复地震，主要取决于地震的破裂尺度以及地震对之间的真实空间间隔。因此，有必要对地震波传播过程进行数值模拟，进而对地震对空间距离和互相关系数进行定量研究，从而获取可靠的重复地震判定结果。

2.   谱元法数值模拟地震波传播

数值模拟地震波传播过程是定量地研究地震学最有效的方法，而计算机的快速发展使得计算量庞大的三维地震波数值模拟成为可能。谱元法将有限元与谱方法相结合，在保证较高计算精度的同时大量减少计算，适用于处理复杂界面模型下的弹性波数值模拟（Komatitsch，Vilotte，1998；Komatitsch，Tromp，1999）。但其网格划分比较复杂。然而随着各种网格划分软件的成熟应用，网格生成也更加简单化，从而使谱元法的应用更具可操作性.美国普林斯顿大学Jeroen Tromp研究组开发的SPECFEM3D_Cartesian程序包在区域地震波模拟中被广泛应用（Komatitsch，Vilotte，1998）。

不同地震事件对的地震波互相关系数是判断重复地震的重要依据。总体来说，事件对的间距越小，其对应地震波的相似性越高。但是，具体的互相关系数与地震对位置之间的关系，还取决于地震台站位置、震源信息（震源破裂过程和震源机制解等）的相似程度及岩石圈内非均匀体的分布等。本文主要讨论三维速度结构下相邻地震激发地震波的互相关系数分布情况，进而为重复地震判定标准提供依据。

岩石圈内部非均匀体能够在地震波传播过程中对P波和S波进行改造，同时产生尾波，进而影响地震对相应的波形互相关系数。由于地球内部的各种动力学过程，致使岩石圈内存在各种尺度的非均匀体，包括地壳内部的流体、矿物差异、岩浆侵入、剥蚀以及变质作用，断层与褶皱和板块之间的相互挤压碰撞作用所导致的岩石破碎和沉积层序改变等。大量非均匀体构成复杂的速度及密度结构，因此，通常用统计模型来描述非均匀体分布。具体可以通过功率谱函数或自相关函数（auto-correlation function，简写为ACF）来描述非均匀体的性质及其分布，通常所用的模型有高斯、指数以及卡尔曼自相关函数模型（Sato et al，2012 ）。已有研究表明，利用指数自相关函数可以较好地描述岩石圈内非均匀体模型，即

式中，r代表两点之间的距离，ε为密度或P波、S波速度扰动的均方根，a为自相关尺度。对于某些地区来说，S波扰动可以高达10%（Capon，1974；Wu，Aki，1985；Yoshimotoet al，1997 ；Sato et al，2012 ），而其自相关尺度在不同地区则介于0.1—10 km间变化。

本研究需先生成完全随机分布的均匀体后，再产生符合指数自相关函数的非均匀体分布统计模型。综合Sato等（2012）的研究结果，本文非均匀体模型采用指数自相关函数统计模型。图3为满足自相关尺度a＝2 km以及P波速度扰动均方根为4%的一次模型实现，其最大速度扰动均方根可达10%。一般来说，如果P波速度存在扰动，S波速度和密度也存在一定变化，更为复杂的情况甚至包括各向异性扰动。但三者扰动之间的互相关系数较为复杂，尤其是密度变化与P波速度变化的关系更为复杂。因此，假定S波速度扰动为P波速度扰动的2倍，而密度扰动设定为0。

图  3  P波速度非均匀体模型

Figure  3.  The heterogeneities model of P-wave velocity

下载: 全尺寸图片 幻灯片

3.   合成地震图检测重复地震判断标准

地震对之间的波形互相关系数取决于速度模型以及地震事件震中位置、台站位置等信息。本文选用简单的双层模型叠加指数自相关函数的非均匀体为结构模型，其顶层为30 km地壳，底层为地幔，其速度和密度参考Dziewonski和Anderson （1981）的初步地球参考模型中莫霍面两侧的地壳与地幔结构，同时加入衰减计算，具体数据列于表1。

表  1  本文所用的地震波速度、密度及衰减模型结构

Table  1.  Velocity，density and attenuation model structure of seismic waves used in this study

	厚度/km	vP/（m·s−1）	vS/（m·s−1）	密度/（g·cm−3）	Qμ
地壳	30	6 800	3 900	2 900	80
地幔	10	8 110	4 490	3 380	80

下载: 导出CSV 

| 显示表格

谱元法模拟区域大小为60 km×60 km×40 km （厚度），相应的3个方向上的单元数为128，128，80，因此，其水平和竖直方向的单元尺度为0.47 km和0.5 km。由于背景速度最低为地壳S波速度3 700 m/s，因而，能够分辨的最高频率大概为10 Hz，接近实际数据中一般采用的最高频率（图2）。同时，岩石圈内的弹性波衰减也对地震波传播产生影响，地震学中使用品质因子Q值来描述地震波衰减。本模型将区域内的品质因子Q_μ统一设置为80，而压缩能量品质因子Q_κ设置为9 999。不同地区的Q值差异较大，例如，盆地的衰减一般都远高于正常地壳厚度的衰减（Sams et al，1997 ）。因此，宁夏盆地内部台站可能具有较强衰减。简单将模型Q_μ统一设置为80；除地表为自由边界条件外，其它5个界面的边界条件为完全吸收边界条件，可几乎完全地吸收向边界外传播的地震波（Xie et al，2014 ），因此，边界条件的影响可忽略不计。

地震震源的破裂时间设定为0.1 s （10 Hz）的点源，参考地震的水平位置及其深度固定在模型内x＝10 km和y＝30 km （图4）和z＝10 km。为研究区域内地震波互相关系数的分布情况，如图4所示，在水平面内沿x和y方向每千米设置一个地震台站。对于M_L1.0—5.0的中小地震来说，其破裂尺度约为几十米至几千米。由于地震数量与震级遵循Gutenberg-Richter定律（Gutenberg，Richter，1994），因此，M_L1.0—3.0左右的小地震发生频率较高，相应的重复地震复发周期也较短，故常被用来估算断层的滑动速率。

宁夏地区小震频发，比较适用短周期地震波，故经5—10 Hz带通滤波后，进行互相关计算来判定重复地震，结果如图2所示。Baisch等（2008）提出1/4优势波长准则，图4a给出了水平距离为97.5 m （地壳中周期为0.1 s的S波波长的1/4）的两个地震事件所对应的互相关系数分布情况。互相关系数的计算公式为

This page contains the following errors:

Below is a rendering of the page up to the first error.

This page contains the following errors:

Below is a rendering of the page up to the first error.

即使地震对的空间间隔为半个优势波长（图5a），甚至大于一个优势波长（图5b），某些台站的互相关系数仍然高于0.9，并且可以更明显地看出这些具有较高互相关系数的台站位于P波辐射较强的45°和135°方向。当然，也有一部分台站的互相关系数低于0.3。虽然方位角为45°和135°方向的互相关系数整体均偏高，但两者并不对称，在图5a和5b中方位角45°附近的互相关系数增长得更加明显，这是由于震源附近散射波对方位角45°左右的波形改造较小所致。然而，图5给出的是统计模型的一次具体实现，而实际情况则取决于真实震源附近的散射体分布情况。

综上，尾波比直达波对非均匀体更为敏感，但其强度取决于非均匀体的强度，那么互相关系数也自然应当与非均匀体的强度直接相关。图4和图5中的非均匀体模型参数为a＝2 km和P波的速度拢动均方根ε＝4%，如果保持非均匀体的分布或相位信息不变，其它模型及震源参数也保持一致，其P波的ε从4%降低为2%，则其互相关系数如图6所示。通过比较图6a与图4a以及图6b与图5a，可以看出当P波的ε＝2%时，非均匀体模型所对应的互相关系数总体要比ε＝4%的模型高，其分布也与P波的辐射花样具有一定的相关性，尽管某些特殊台站的结果可能不一致。

图  6  震源间距为97.5 m （a）和195 m （b）时的波形互相关系数分布

自相关尺度a＝2 km，P波和S波的速度扰动均方根分别为2%和4%。图中震源机制的空间分布或相位信息均与图4和图5一致

Figure  6.  Waveform cross-correlation coefficients distribution with the distance separation of earthquake pairs of 97.5 m （a） and 195 m （b）

The auto-correlation length a is 2 km，the P-wave and S-wave velocity perturbations are 2% and 4%，respectively. The focal mechanism of the earthquakes used here is exactly the same as Fig.4 and Fig.5

下载: 全尺寸图片 幻灯片

图  5  地震对震源间距为195 m （a）和290 m （b）时波形互相关系数分布

非均匀体模型参数及震源机制解与图4相同

Figure  5.  Waveform cross-correlation coefficients distribution with the distance separation of earthquake pairs of 195 m （a） and 290 m （b）

The parameters of heterogeneities model and the focal mechanism of the earthquakes are exactly the same as Fig. 4

下载: 全尺寸图片 幻灯片

直达波相对尾波的振幅可对互相关系数产生较大影响，由于直达波受非均匀体的影响比尾波要弱，因此，在直达波振幅较大的台站，其互相关系数也相对较高。走滑型地震的P波辐射在水平面投影内的极小值处于走向及垂直走向的方位角方向上，而逆冲型地震的直达P波在震源附近几乎所有方位角方向上均具有较强的辐射。图7所对应的模型结构与图4a完全一致，但震源机制由走滑型地震变为逆冲型地震。对比90°方位角的震中距剖面 （图7b与图4b），其对应的直达P波相对于尾波更强，而互相关系数整体上也比走滑型地震要大，这与前述“如果直达P波相较尾波更强，则对应的互相关系数也更大”的推论一致。

图  7  震源间距97.5 m时的波形互相关系数分布（a）和垂向速度分量合成地震图（b）

非均匀体模型与图4完全相同，但震源机制为倾角45°逆冲型，选取剖面台站为y＝30 km，x在10—55 km范围内递增

Figure  7.  Waveform cross-correlation coefficients distribution （a） and a distance profile of vertical velocity seismograms （b） for an earthquake pair with a distance separation of 97.5 m

The heterogeneity model used here is exactly the same as Fig.4，but the earthquake focal mechanism is replaced by a thrust type with a dip angle of 45°，this profile is composed of synthetics from stations at y＝30 km and x ranging from 10 km to 55 km

下载: 全尺寸图片 幻灯片

图  4  走滑型地震对应的波形互相关系数分布（a）及部分台站波形合成图（b）

自相关尺度a＝2 km，P波和S波的速度扰动均方根分别为4%和8%。图（a）中的地震对间距为97.5 m，震源机制为走向0°的走滑型，图（b）中垂直向速度分量为5—10 Hz带通滤波后的震中距剖面，选取剖面台站为y＝30 km，x从10 km递增至55 km

Figure  4.  Waveforms cross-correlation coefficients （a） and waveform synthetics of station （b） for a strike-slip earthquake

The auto-correlation length a is 2 km，the P-wave and S-wave velocity perturbations are 4% and 8% respectively. The earth-quake pair has a distance separation of 97.5 m，and the earthquake has a strike-slip focal mechanism with a strike of 0° in Fig.（a），Fig.（b） shows a distance profile of vertical velocity seismograms after 5−10 Hz band-pass filtering，and the profile is composed of synthetics from the stations at y＝30 km and x ranging from 10 km to 55 km

下载: 全尺寸图片 幻灯片

4.   讨论与结论

重复地震研究在断层深部滑动速率估算以及地球内部介质时间变化监测方面具有重要应用，而准确地判定重复地震是深入研究重复地震的基础。本文以宁夏区域台网为例，探讨了基于稀疏台网记录的波形互相关法在判定重复地震中的优势，并进一步采用谱元法模拟三维结构下地震波的传播，分析了波形互相关系数分布情况与非均匀体性质之间的关系，本文通过理论分析及数值模拟得到如下结论：

1） 在稀疏台网情况下，例如宁夏区域台网（平均台间距45 km），无论是绝对定位还是相对定位的精度都不足以准确地判定重复地震。波形互相关法对台站分布要求较低，对地震事件对的空间间距更为敏感，因此广泛地应用于重复地震判定。

2） 波形互相关系数对地震间距的敏感性主要来源于震源附近的三维非均匀体。通过建立三维非均匀体统计模型，我们使用谱元法数值模拟三维结构下的地震波传播并分析了波形互相关系数与地震对间距的关系。研究表明，互相关系数分布受到速度模型结构（主要是震源附近三维非均匀体性质）、震源机制解和地震对间距的影响。具体来说，震源附近非均匀体强度越大，则散射效应越强，进而互相关系数对地震间距更敏感；由于尾波比直达P和S波对震源附近非均匀体更敏感，因此，互相关系数在不同台站的分布与直达波和尾波的振幅比（即震源辐射花样）有关，直达波越强，则互相关系数越高。

3） Geller和Mueller （1980）提出利用1/4优势波长准则来判定重复地震，1/4优势波长准则虽具有统计上的意义，但具体互相关系数的判定阀值受到很多因素的影响。本文研究表明，判定重复地震的波形互相关系数至少与震源附近非均匀体性质和震源机制解有关，研究地壳内的非均匀体性质及震源机制解有助于重复地震的判定。如果震源附近结构简单或者直达波较强，则互相关系数对震源间距的敏感性更低，实际判定过程中应当选取更高的互相关系数阀值或者只截取尾波进行互相关计算。

当然，实际观测还受到干扰噪音的影响。对于板内重复地震研究来说，主要选取M_L1.0—3.0的重复地震，对宁夏区域台网来说，一般都能在至少两个以上台站观测到较高信噪比的波形记录。对于大于M_L3.0的重复地震，其复发周期较长，一般很难在几年或者十几年的时间尺度内观测到；而对于震级小于M_L1.0的重复地震，其破裂尺度更小，信号更弱，导致干扰噪声影响更大，针对其的判定对台网密度和台站波形数据质量要求更高。

综上，对于平均台间距30—50 km左右的区域台网来说，利用波形互相关法可有效地判定约为M_L1.0—3.0的重复地震。但是，判定过程中的互相关系数阀值选取与台站所处的震源辐射花样及地壳非均匀体结构有关，对于直达波振幅较强的波形记录，应当选取更高的阀值（例如高于0.95）或者只截取尾波进行互相关系数的计算。

宁夏地震台网中心、甘肃地震台网中心和内蒙古地震台网中心提供了地震波形数据，美国普林斯顿大学的Jeroen Tromp研究组提供了数据处理所需的SPECFEM3D_Cartesian程序，该大学的吴文波博士对本文的撰写提供了指导，作者在此一并表示感谢。