引言

通过理论分析和数值计算得到的理论地震图是地震学家们了解震源机制和地球介质性质的重要基础。自Lamb （1904）开展了自由界面对于施加在其上的力的响应的经典工作以来，经过一个多世纪的发展，特别是伴随着近50年来计算技术的迅猛发展，合成地震图技术已经取得了长足的进步。20世纪50年代以前，地震学家和数学家们研究介质响应多是采取渐近分析手段得到近似解（Lapwood，1948；Pekeris，1955），而从50年代计算机引入到地震学领域中以来，各种算法不断涌现。根据解析程度的大小，大致可分成半解析-半数值方法和数值方法两类。对于半解析-半数值方法，主要有基于Cagniard-de Hoop技术的广义射线方法（Helmberger，1973，1974；Johnson，1974）和基于矩阵运算的反射率方法（Fuchs，Müller，1971；Kennett，Kerry，1979；Luco，Apsel，1983；Chen，1999）。对于数值方法，各种不同的计算方式层出不穷，比如有限差分法（Madariaga，1976；Levander，1988；Graves，1996；Zhang，Chen，2006）、谱元法（Priolo，Seriani，1991；Komatitsch et al，2001；Capdeville et al，2003；Komatitsch，Tromp，2003；Tromp，Komatitsch，2008）等。因数值方法对介质的适用性更广，越来越受到地震学家们的重视。但是在很多情况下，如果介质可用一维结构近似，半解析-半数值的方法，特别是能够合成完整地震波场的反射率方法，由于其计算效率方面的优势仍然会得到广泛的应用。

在反射率方法中，基于离散波数法的广义反透射系数法（Apsel，Luco，1983；Yao，Harkri-der，1983；Chen，1999）是目前应用最广泛的方法。但是，当震源和接收点的深度相近或相同时，由于被积函数收敛得非常慢，仅仅用离散波数法而不采取其它特殊技术时计算效率很低。Apsel和Luco （1983）及Hisada （1994，1995）采用了半解析的方法，通过将积分拆成两部分的方式巧妙地解决了被积函数的慢收敛问题，但缺点是数学处理非常繁琐，且在实际的数值计算中需要特殊的处理技术。基于重复平均法的原则，张海明等（2001）提出了一种计算浅源情况下的数值积分方法，即峰谷平均法。与半解析方法相比，该方法不仅数学处理上非常简单，而且易于数值实现。

运用峰谷平均法，需要在波数k大于阈值k_c之后构造由积分值变化曲线的波峰和波谷组成的振荡序列，而这个k_c在之前的研究（张海明等，2001）中是由经验公式给出的。但是，一方面，经验公式缺乏理论依据；另一方面，在有些情况下，根据经验公式确定的k_c比实际的积分限阈值大得多，导致计算效率不高。因此有必要思考新的求解方法。在本文的研究中，我们发现波数积分的积分限阈值k_c与反透射系数行列式的零点有关，拟通过求解频率域中瑞雷函数的零点，找到可以精确地确定积分限阈值的方案，以期通过精确确定的积分阈值能够以更高的效率计算理论地震图。

1.   峰谷平均法

根据广义反透射系数法，分层半空间中任意点源产生的位移在频率域中可表示为如下类型的波数积分（张海明等，2001）：

式中，r为震中距，z和z_s分别为震源和接收点的深度，k为水平波数，ω为角频率，J_n（kr）为n阶贝塞尔函数，F（r，z，z_s，ω）是包含衰减因子${{\rm{e}}^{ {\text{－}}c(k{\text{，}} \omega) |z{\text{－}} {z_{\rm{s}}}|}}$的项，其中c（k，ω）为与波数k和频率ω有关的实部为正的复数因子。当震源和接收点的深度相近或相同时，衰减因子作用很小，在介质响应函数和贝塞尔函数的共同作用下，整个被积函数衰减得很慢。为了得到精确的结果而不得不取很大的积分上限，这样就使得数值积分的效率很低。但是，当波数k大于某一阈值k_c时，积分值随着积分上限的增大而在某一极限值附近振荡，且其包络线为单调光滑的曲线。将波峰与波谷值抽取出来组成一个慢收敛序列，利用重复平均法便可快速求出其收敛值（Zhang et al，2003）。因此，重要的一步便是求出阈值k_c。以往的工作中，阈值k_c由经验公式给出（张海明等，2001），

式中，v_min为介质模型的最小速度，ω和ω_I分别为角频率和虚频率。然而为了尽量照顾到所有情况，通过经验公式求出的k_c可能会取得过大，从而牺牲了计算效率。

2.   反透射系数行列式零点的求解

在广义反透射系数法中（Chen，1999），为描述各层之间地震波场振幅的关系，引入了修正反透射系数。由于P-SV波和SH波的修正反透射系数和广义反透射系数具有相同的形式，且我们所研究的瑞雷波为P波和SV波耦合所形成，故下文的讨论只针对P-SV波进行。对于一个N层的分层介质，P-SV波的修正反透射系数满足[Chen （1999）中式（41）]

式中，j表示层数，取值为1，2，···，N−1，R为修正反射系数，T为修正透射系数。当j＝N时，修正反透射系数满足[Chen （1999）中式（42）]

进一步，为描述各层介质中地震波场振幅对某一层中地震波场振幅的影响，引入广义反透射系数，其表达式为:

当j＝1，2，···，s－1时，满足[Chen （1999）中式（45）]

当j＝s，s＋1，···，N时，满足[Chen （1999）中式（46）]

式中，s为源所在层数，${\hat{\boldsymbol R}}$为广义反射系数，${\hat{\boldsymbol T}}$为广义透射系数，I为单位矩阵。在式（5）中存在矩阵求逆的情况，若相应的求逆矩阵存在行列式等于零的情况，此时该项便存在奇异性。式（5）中，${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$，${\hat{\boldsymbol T}}_{\rm{u}}^{\left( {{j}} \right)}$，${\hat{\boldsymbol T}}_{\rm{d}}^{\left( {{j}} \right)}$由于存在矩阵求逆，因此，这三项均有可能存在某个波数使得分母为零。下文将对这三项作进一步讨论，以寻找最适于用来求解波数积分限阈值k_c的项。

以一个三层介质模型为例（Chen，1993），如表1所示。震源深度z为29.95 km，接收点深度z_s为30.05 km。分别求出不同频率下${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$，${\hat{\boldsymbol T}}_{\rm{u}}^{\left( {j} \right)}$，${\hat{\boldsymbol T}}_{\rm{d}}^{\left( {j} \right)}$各项中含有负一次幂矩阵的行列式的零点，即分别求出${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$中所含${\boldsymbol{E}}_{21}^{\left( 1 \right)}$，${\hat{\boldsymbol T}}_{\rm{u}}^{\left( {j} \right)}$中所含${\boldsymbol{I}} {\text{－}} {\boldsymbol{R}}_{{\rm{du}}}^{\left( {j} \right)}{\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( {{j} {\text{－}} 1} \right)}$和${\hat{\boldsymbol T}}_{\rm{d}}^{\left( {j} \right)}$中所含${\boldsymbol{I}}{\text{－}} {\boldsymbol{R}}_{{\rm{ud}}}^{\left( {j} \right)}{\hat{\boldsymbol R}}_{{\rm{du}}}^{\left( {{j} {\text{＋}} 1} \right)}$这三项的行列式所对应的零点。在这个三层模型中，由于s＝1，N＝3，则

表  1  不含低速层的三层地壳模型参数 （Chen，1993）

Table  1.  Parameters of the three-layer crust model without low-velocity-layers（after Chen，1993）

层序	层底深度/km	ρ/（g·cm−3）	vS/（km·s−1）	vP/（km·s−1）
1	20	2.8	3.50	6.0
2	30	2.9	3.65	6.3
3	45	3.1	3.90	6.7
4	∞	3.3	4.70	8.2

下载: 导出CSV 

| 显示表格

第三层广义透射系数${\hat{\boldsymbol T}}$对应的式子中并无矩阵求逆，故只有第一层和第二层的广义透射系数，即${\hat{\boldsymbol T}}_{\rm{u}}^{\left( {\rm{1}} \right)}$和${\hat{\boldsymbol T}}_{\rm{d}}^{\left( {\rm{2}} \right)}$，有贡献。由${\hat{\boldsymbol T}}_{\rm{u}}^{\left( {\rm{1}} \right)}$和${\hat{\boldsymbol T}}_{\rm{d}}^{\left( {\rm{2}} \right)}$这两项求得的零点有多个，只取其零点的最大值，分别记作k_t1和k_t2，同时将通过${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$求得的零点记为k_r，则可以得到在不同频率下，各项零点与积分值的关系，如图1所示。

图  1  不同频率下反透射系数的零点与位移的关系（位移U_r为归一化后的结果）

图（a）−（d）分别对应f＝0.1，0.4，3，5 Hz的不含低速层的三层地壳模型；图（e）和（f）分别对应f＝1，5 Hz的含低速层的五层地壳模型

Figure  1.  The relationship between zeros of R/T coefficient and the displacement in different frequencies with the normalized displacement U_r

Figs. （a） to （d） correspond to the three-layer crust model without low-velocity-layers when f＝0.1，0.4，3 and 5 Hz；Figs. （e） and （f） correspond to the five-layer crust model with low-velocity-layers when f＝1 Hz and 5 Hz

下载: 全尺寸图片 幻灯片

从图1中可以看出，对于不同的频率取值，通过${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$求得的零点k_r均出现在积分值产生振荡行为的位置附近。再观察由透射系数求得的k_t1和k_t2的情况：在f＝0.1 Hz的情况下（图1a），k_t1的值与k_r相近，而k_t2处于积分值产生振荡行为之前，k_t2不可作为积分限阈值；在f＝0.4 Hz的情况下（图1b），k_t2能够满足作为积分限阈值的条件，但是k_t2并不适用；在f＝3 Hz的情况下（图1c），k_t2的值与k_r相近，而k_t1处于积分值产生振荡行为之前，k_t1不可作为积分限阈值；在f＝5 Hz的情况下（图1d），k_t2适用，k_t1不适用。再考虑一个含低速层的五层地壳模型（何耀锋等，2006），如表2所示，其结果如图1e，f所示。在f＝1 Hz的情况下（图1e），k_t2适用，而k_t1不适用；在f＝5 Hz的情况下（图1f），k_t1可作为积分限阈值，而k_t2远大于k_t1和k_r，若取作积分限阈值，则计算效率相对较低，因此作为积分限阈值不合适。在柱坐标系下，位移可以分解为三个分量，分别为径向分量U_r，切向分量U_f和垂直分量U_z。图1只显示了U_r的实部，对于U_r的虚部和另外两个方向位移分量的实部和虚部均可得到相似的结果。可以看出在所有情况中，通过反射系数${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$求得的零点k_r均满足其作为积分限阈值的要求，而由透射系数求得的零点只在某些情况下满足。因此我们可以认为通过${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$获得的零点k_r作为积分限阈值更为合理。由此我们便找到了求解积分限阈值的另一个方法，即求解广义反射系数${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$中使逆矩阵的行列式为零的点k_r，将k_r作为积分限阈值。

表  2  含低速层的五层地壳模型的参数（何耀锋等，2006）

Table  2.  Parameters of the five-layer crust model with low-velocity-layers（after He et al, 2006）

层序	层底深度/km	ρ/（g·cm−3）	vS/（km·s−1）	vP/（km·s−1）
1	20	2.8	3.50	6.0
2	30	2.9	3.65	6.3
3	45	3.5	2.50	4.9
4	70	4.0	3.80	6.5
5	90	3.2	4.50	7.6
6	∞	2.4	5.60	9.3

下载: 导出CSV 

| 显示表格

对${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$进一步进行研究，得到各项的含义为

带入式（5）中第一个子式进行计算，可得瑞雷函数的表达式为

式中，$x {\text{＝}} {({\beta ^{\left( 1 \right)}})^2}{\left( {k/\omega } \right)^2}$，$m {\text{＝}} {( {{\beta ^{(1)}}})^2}\Bigr/{( {{\alpha ^{(1)}}} )^2}$。当m的取值范围为（0，0.5）时，该瑞雷函数有且仅有一个大于1的实数根（Liu et al，2016）。因此广义反射系数${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$求逆矩阵的行列式的零点可以通过求解自由表面瑞雷函数的零点来获得。

3.   应用求零点法求解积分限阈值

3.1   数值算例

仍然采用表1所示的三层地壳模型。将各参数带入式（8），求出瑞雷函数的零点。同时，画出频率域下三个位移分量积分值随积分上限的变化曲线，得到瑞雷函数零点与积分值的关系，如图2所示。可以清楚地看出，积分值的震荡行为产生于瑞雷函数零点附近。因此有理由认为，通过求瑞雷函数零点来求解k_c的方式是正确的。至此，本文便可以给出基于理论分析求解k_c的方法。在计算含有衰减因子的积分时，若需求得积分限阈值k_c，即可通过求取广义反透射系数求逆矩阵行列式的零点的方式求得，并且在实际计算中，为了保证构成慢收敛序列的最初几个波峰和波谷值的稳定，k_c的取值将会稍大于零点值。我们将这种通过求瑞雷函数零点来寻求积分限阈值的方法定义为求零点法。

图  2  f＝2 Hz时波速递增的三层地壳模型下位移三分量的实部、虚部图（位移U为归一化后的结果）

图（a），（b）分别为径向位移U_r的实部和虚部；图（c），（d）分别为切向位移U_f的实部和虚部；图（e），（f）分别为垂向位移U_z的实部和虚部

Figure  2.  The real part and image part of three components of displacement in the three-layer crust model without low-velocity-layers with the normalized displacement U when f＝2 Hz

Figs. （a） and （b） are the real part and image part of the radial direction displacement （U_r），respectively；Figs. （c） and （d） are the real part and image part of the tangential displacement （U_f），respectively；Figs. （e） and （f） are the real part and image part of the vertical displacement （U_z），respectively

下载: 全尺寸图片 幻灯片

3.2   求零点法与经验公式法的对比

与之前通过经验公式找积分限阈值k_c的方法相比，通过求零点法来寻求k_c具有一定的优势，其寻求的k_c更为准确，而且针对某些特殊的情形，求零点法比经验公式法更实用，它可以更精确地找到临界点，并在保证计算精度的前提下具有更高的计算效率。

仍然考虑一个如表1所示的三层地壳模型，震源仍置于地下深度为29.95 km处，接收点仍置于地下深度为30.05 km处，分别通过求零点法和经验公式法求出积分限阈值。为了区别两种方法求出的积分限阈值，我们将通过经验公式求得的积分限阈值记为k_e，将通过求零点法求得的积分限阈值记作k_c。结果如图3所示，可以清楚地看到：两种方法在求得积分限阈值后，积分值均具有震荡行为；由求零点法得到的积分限阈值k_c更接近积分值开始产生震荡的位置，在运算的过程中具有更高的效率；而由经验公式法得到的积分限阈值k_e比k_c大，意味着需要更大的计算量。

图  3  波速递增的三层地壳模型在不同频率下k_e与k_c的关系（位移U_r为归一化后的结果）

Figure  3.  The relationship between k_c and k_e at different frequencies in the three-layer crust model without low-velocity-layers with the normalized displacement U_r

下载: 全尺寸图片 幻灯片

若考虑如表2所示存在低速层的情况，得到的结果如图4所示。可以看出，在存在低速层的情况下，由经验公式法得到的积分限阈值k_e不得不取得很大，说明在该情况下，经验公式法的适用性较差。同时，也可以看到由求零点法得到的积分限阈值k_c仍然出现在产生震荡行为的积分值附近，与本文结果十分符合。

图  4  f＝5 Hz时含低速层的五层地壳模型下k_e与k_c的关系（位移U_r为归一化后的结果）

Figure  4.  The relationship between k_e and k_c in the five-layer crust model with low-velocity-layers while the displace-ment has been normalized when f＝5 Hz

下载: 全尺寸图片 幻灯片

3.3   与Zhang等 （2003） 结果的对比

为了进一步检验结果的正确性，采用求零点法和经验公式法（Zhang et al，2003）分别求得理论地震图，并对结果进行对比。仍然考虑如表2所示的含低速层的5层地壳模型，震源和接收点均位于地表，震中距为200 km，结果如图5所示（图中的位移为归一化后所得）。可以看到，两种方法得到的结果相同。但是利用求零点法计算理论地震图需要12.74 s，所用计算机处理器为双核i3CPU，内存为4G；而用原有的经验公式法计算理论地震图则需要14.693 s，求零点法的加速效率为13.3%。该结果充分证明了求零点法的正确性，同时也说明该方法的确能够提升理论地震图的计算效率。

图  5  利用求零点法和经验公式法分别获得的理论地震图径向位移U_r （a）、切向位移U_f （b）和垂向位移U_z （c）三分量归一化图

Figure  5.  Seismogram calculated by for-zero method and empirical-formula method respectively to obtain the radial direction displacement U_r，the tangential displacement U_f，and the vertical displace-ment U_z，which are normalized

下载: 全尺寸图片 幻灯片

4.   讨论与结论

通过比较不同深度的积分图像（图6）可以看出：当震源与接收点的相对深度越小，即震源与接收点垂直方向相对距离越小，积分收敛得越慢；当相对深度越大，即震源与接收点垂向相对距离越大，积分收敛则越快。因此，通过求零点法寻找积分限阈值的方法适用于处理震源与接收点深度相近的情况。

图  6  波速递增的三层地壳模型下位移在不同震源与接收点相对深度h时的变化图（位移U_r为归一化后的结果）

Figure  6.  Relative depth between source and receiver in the three-layer crust model without low-velocity-layers while the displacements U_r have been normalized

下载: 全尺寸图片 幻灯片

在利用峰谷平均法计算合成理论地震图的波数积分时，找到一个合适的k_c，是此方法应用的前提。本文从广义反透射系数法出发，发现由广义反射系数${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$中求逆矩阵行列式的零点作为积分限阈值最合适，并由该行列式得到了分层介质中自由表面的瑞雷函数。同时，通过对瑞雷函数零点的求解，从图像上寻求了k_c与瑞雷函数零点的关系。在不同频率、不同层数的情况下，采用本文的方法均能很好地解决求解积分限阈值的问题。同时，只要震源和接收点的相对深度不大，也即垂向相对距离使得积分值的收敛足够慢，则无论震源和接收点是否在同一层，皆不会影响求零点法的应用。与前人所用的由经验公式求解积分限阈值的方法相比，求零点法能够获取更小的积分限阈值，这样使得我们能够以更高的效率计算理论地震图；同时对于一些用经验公式求解不适用的特殊情形，例如存在低速层的情况，使用求零点法得到的阈值计算效率更高。