利用反透射系数确定理论地震图计算过程中的积分限阈值

吴永祺, 张海明

吴永祺,张海明. 2018. 利用反透射系数确定理论地震图计算过程中的积分限阈值. 地震学报,40(6):719−727. doi:10.11939/jass.20180013. DOI: 10.11939/jass.20180013
引用本文: 吴永祺,张海明. 2018. 利用反透射系数确定理论地震图计算过程中的积分限阈值. 地震学报,40(6):719−727. doi:10.11939/jass.20180013. DOI: 10.11939/jass.20180013
Wu Y Q,Zhang H M. 2018. Determining the threshold value in upper limit of wavenumber integration by reflection-transmission coefficient in theoretical seismograms calculation. Acta Seismologica Sinica40(6):719−727. doi:10.11939/jass.20180013. DOI: 10.11939/jass.20180013
Citation: Wu Y Q,Zhang H M. 2018. Determining the threshold value in upper limit of wavenumber integration by reflection-transmission coefficient in theoretical seismograms calculation. Acta Seismologica Sinica40(6):719−727. doi:10.11939/jass.20180013. DOI: 10.11939/jass.20180013

利用反透射系数确定理论地震图计算过程中的积分限阈值

基金项目: 国家自然科学基金(41674050)资助
详细信息
    通讯作者:

    张海明: e-mail:zhanghm@pku.edu.cn

  • 中图分类号: P315.01

Determining the threshold value in upper limit of wavenumber integration by reflection-transmission coefficient in theoretical seismograms calculation

  • 摘要: 本文尝试基于理论分析来求解积分限阈值kc,即通过反透射系数来确定kc。根据计算理论地震图的广义反透射系数法,在反透射系数中进行求逆运算的矩阵行列式的零点将会使得被积函数产生较大的变化,通过具体实例显示出自由界面处的反射系数中含有的零点可以作为合适的kc。多种情况的实例显示,通过反射系数来确定kc具有较好的普适性。通过与经验公式的计算结果进行比对,表明根据本文方案所确定的kc在保证准确性的同时可以有效地提高计算效率。
    Abstract: In multi-layered half-space, the displacement field is expressed as the integration of wavenumber k using reflection-transmission (R/T) coefficient when the seismogram is calculated in numerical ways in the frequency domain. Therefore, some special methods were introduced to solve those kinds of integration as to accelerate computation. For example, when the depth of source is equal or close to that of receiver, the peak-trough averaging method (PTAM) is applied. To apply PTAM, however, one must determine the threshold value called kc after which the integration oscillates regularly. In previous studies, kc was estimated empirically without theoretical support. In this study, a scheme based on theoretical analysis to determine kc is proposed, and kc is related to the R/T coefficients. According to the generalized R/T coefficient method, violent variation of integrand will occur when the determinant of relevant R/T coefficient vanishes. We show by examples that for a given frequency, root of the determinant of the reflection coefficient on the free surface is a proper value for kc. Compared with the method calculated with empirical formula, the new method to determine kc will be more efficient for a given precision to prove the validity.
  • 通过理论分析和数值计算得到的理论地震图是地震学家们了解震源机制和地球介质性质的重要基础。自Lamb (1904)开展了自由界面对于施加在其上的力的响应的经典工作以来,经过一个多世纪的发展,特别是伴随着近50年来计算技术的迅猛发展,合成地震图技术已经取得了长足的进步。20世纪50年代以前,地震学家和数学家们研究介质响应多是采取渐近分析手段得到近似解(Lapwood,1948Pekeris,1955),而从50年代计算机引入到地震学领域中以来,各种算法不断涌现。根据解析程度的大小,大致可分成半解析-半数值方法和数值方法两类。对于半解析-半数值方法,主要有基于Cagniard-de Hoop技术的广义射线方法(Helmberger,19731974Johnson,1974)和基于矩阵运算的反射率方法(Fuchs,Müller,1971Kennett,Kerry,1979Luco,Apsel,1983Chen,1999)。对于数值方法,各种不同的计算方式层出不穷,比如有限差分法(Madariaga,1976Levander,1988Graves,1996Zhang,Chen,2006)、谱元法(Priolo,Seriani,1991Komatitsch et al,2001Capdeville et al,2003Komatitsch,Tromp,2003Tromp,Komatitsch,2008)等。因数值方法对介质的适用性更广,越来越受到地震学家们的重视。但是在很多情况下,如果介质可用一维结构近似,半解析-半数值的方法,特别是能够合成完整地震波场的反射率方法,由于其计算效率方面的优势仍然会得到广泛的应用。

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

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

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

    $$ Y\left( {r{\text{,}}z{\text{,}}{{{z}}_{\rm{s}}}{\text{,}}\omega } \right) {\text{=}} \mathop \int \nolimits_0^{ + \infty } F\left( {r{\text{,}}z{\text{,}}{{{z}}_{\rm{s}}}{\text{,}}\omega } \right){{\rm{J}}_n}\left( {kr} \right){\rm{d}}k{\text{,}} $$ (1)

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

    $$ {{{k}}_{\rm{c}}} {\text{=}} 1.5\frac{{\sqrt {{\omega ^2} {\text{+}} \omega _{\rm{I}}^2} }}{{{v_{{\rm{min}}}}}}{\text{,}} $$ (2)

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

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

    $$ \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{T}}_{\rm{d}}^{{{(j)}}}} & {{\boldsymbol{R}}_{{\rm{ud}}}^{{{(j)}}}}\\ {{\boldsymbol{R}}_{{\rm{du}}}^{{{(j)}}}} & {{\boldsymbol{T}}_{\rm{u}}^{{{(j)}}}} \end{array}} \right]{\rm{ {\text{=}} }}{\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{E}}_{{\rm{11}}}^{{{(j {\text{+}} 1)}}}} & { {\text{-}} {\boldsymbol{E}}_{{\rm{12}}}^{{{(j)}}}}\\ {{\boldsymbol{E}}_{{\rm{21}}}^{{{(j {\text{+}} 1)}}}} & { {\text{-}} {\boldsymbol{E}}_{{\rm{22}}}^{{{(j)}}}} \end{array}} \right]^{{\rm{ {\text{-}} 1}}}}\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{E}}_{{\rm{11}}}^{{{(j)}}}} & { {\text{-}} {\boldsymbol{E}}_{{\rm{12}}}^{{{(j {\text{+}} 1)}}}}\\ {{\boldsymbol{E}}_{{\rm{21}}}^{{{(j)}}}} & { {\text{-}} {\boldsymbol{E}}_{{\rm{22}}}^{{{(j {\text{+}} 1)}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{K}}_{\rm{d}}^{\left( {{j}} \right)}{\rm{(}}{{{z}}^{{{(j)}}}}{\rm{)}}} & {\rm{0}}\\ {\rm{0}} & {{\boldsymbol{K}}_{\rm{u}}^{\left( {{{j {\text{+}} 1}}} \right)}{\rm{(}}{{{z}}^{{{(j)}}}}{\rm{)}}} \end{array}} \right], $$ (3)

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

    $$ \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{T}}_{\rm{d}}^{{{(N)}}}}\\ {{\boldsymbol{R}}_{{\rm{du}}}^{{{(N)}}}} \end{array}} \right]{\rm{ {\text{=}} }}{\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{E}}_{{\rm{11}}}^{{{(N {\text{+}} 1)}}}} & { {\text{-}} {\boldsymbol{E}}_{{\rm{12}}}^{{{(N)}}}}\\ {{\boldsymbol{E}}_{{\rm{21}}}^{{{(N {\text{+}} 1)}}}} & { {\text{-}} {\boldsymbol{E}}_{{\rm{22}}}^{{{(N)}}}} \end{array}} \right]^{ {\text{-}} {\rm{1}}}}\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{E}}_{{\rm{11}}}^{{{(N)}}}{\boldsymbol{K}}_{\rm{d}}^{\left( {{N}} \right)}{\rm{(}}{{{z}}^{{{(N)}}}}{\rm{)}}}\\ {{\boldsymbol{E}}_{{\rm{21}}}^{{{(N)}}}{\boldsymbol{K}}_{\rm{d}}^{\left( {{N}} \right)}{\rm{(}}{{{z}}^{{{(N)}}}}{\rm{)}}} \end{array}} \right]{\text{,}} $$ (4)

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

    j=1,2,···,s-1时,满足[Chen (1999)中式(45)]

    $$ \left\{ \begin{aligned} &{\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}{\rm{ {\text{=}} }} {\text{-}}{\boldsymbol{E}}{_{{\rm{21}}}^{\left( {\rm{1}} \right)^{\text{-}1}}}{\boldsymbol{E}}_{{\rm{22}}}^{\left( {\rm{1}} \right)}{\boldsymbol{K}}_{\rm{u}}^{\left( {\rm{1}} \right)}{\rm{(}}{{{z}}^{{\rm{(0)}}}}{\rm{)}}{\text{,}}\\ &{\hat{\boldsymbol T}}_{\rm{u}}^{\left( {{j}} \right)}{\rm{ {\text{=}} }}{[{\boldsymbol{I}} {\text{-}} {\boldsymbol{R}}_{{\rm{du}}}^{\left( {{j}} \right)}{\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( {{{j}} {\text{-}} 1} \right)}]^{ {\text{-}}1}}{\boldsymbol{T}}_{\rm{u}}^{\left( {{j}} \right)}{\text{,}}\\ &{\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( {{j}} \right)} {\text{=}} {\boldsymbol{R}}_{{\rm{ud}}}^{\left( {{j}} \right)} {\text{+}}{\boldsymbol{T}}_{\rm{d}}^{\left( {{j}} \right)}{\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( {{{j}} {\text{-}} 1} \right)}{\hat{\boldsymbol T}}_{\rm{u}}^{\left( {{j}} \right)}{\text{;}} \end{aligned} \right. $$ (5)

    jss+1,···,N时,满足[Chen (1999)中式(46)]

    $$ \left\{ \begin{aligned} &{\hat{\boldsymbol R}}_{{\rm{du}}}^{\left( {{{N}} {\text{+}} 1} \right)}{\rm{ {\text{=}} }}0{\text{,}}\\ &{\hat{\boldsymbol T}}_{\rm{d}}^{\left( {{j}} \right)}{\rm{ {\text{=}} }}{[{\boldsymbol{I}} {\text{-}} {\boldsymbol{R}}_{{\rm{ud}}}^{\left( {{j}} \right)}{\hat{\boldsymbol R}}_{{\rm{du}}}^{\left( {{{j}} {\text{+}} 1} \right)}]^{ {\text{-}} 1}}{\boldsymbol{T}}_{\rm{d}}^{\left( {{j}} \right)}{\text{,}}\\ &{\hat{\boldsymbol R}}_{{\rm{du}}}^{\left( {{j}} \right)} {\text{=}} {\boldsymbol{R}}_{{\rm{du}}}^{\left( {{j}} \right)} {\text{+}} {\boldsymbol{T}}_{\rm{u}}^{\left( {{j}} \right)}{\hat{\boldsymbol R}}_{{\rm{du}}}^{\left( {{{j}} {\text{+}} 1} \right)}{\hat{\boldsymbol T}}_{\rm{d}}^{\left( {{j}} \right)}{\text{,}} \end{aligned} \right. $$ (6)

    式中,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)}$由于存在矩阵求逆,因此,这三项均有可能存在某个波数使得分母为零。下文将对这三项作进一步讨论,以寻找最适于用来求解波数积分限阈值kc的项。

    以一个三层介质模型为例(Chen,1993),如表1所示。震源深度z为29.95 km,接收点深度zs为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 R}}_{{\rm{du}}}^{\left( 4 \right)}{\rm{ {\text{=}} }}0{\text{,}} {\hat{\boldsymbol T}}_{\rm{d}}^{\left( {\rm{3}} \right)} {\text{=}} {\boldsymbol{T}}_{\rm{d}}^{\left( {\rm{3}} \right)}{\text{.}} $$ (7)

    第三层广义透射系数${\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)}$这两项求得的零点有多个,只取其零点的最大值,分别记作kt1kt2,同时将通过${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$求得的零点记为kr,则可以得到在不同频率下,各项零点与积分值的关系,如图1所示。

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

    表  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)}$进一步进行研究,得到各项的含义为

    $$ {\boldsymbol{E}}_{{\rm{21}}}^{{\rm{(1)}}} {\text{=}} \frac{1}{\omega }\left[ {\begin{array}{*{20}{c}} { {\text{-}} 2{\alpha ^{\left( 1 \right)}}{\mu ^{\left( 1 \right)}}k{\gamma ^{\left( 1 \right)}}} & { {\text{-}} {\beta ^{\left( 1 \right)}}{\mu ^{\left( 1 \right)}}{\chi ^{\left( 1 \right)}}}\\ { {\text{-}} {\alpha ^{\left( 1 \right)}}{\mu ^{\left( 1 \right)}}{\chi ^{\left( 1 \right)}}} & { {\text{-}} 2{\beta ^{\left( 1 \right)}}{\mu ^{\left( 1 \right)}}k{\nu ^{\left( 1 \right)}}} \end{array}} \right] {\text{,}} $$ (8)
    $$ {\boldsymbol{E}}_{{\rm{22}}}^{{\rm{(1)}}} {\text{=}} \frac{1}{\omega }\left[ {\begin{array}{*{20}{c}} {2{\alpha ^{\left( 1 \right)}}{\mu ^{\left( 1 \right)}}k{\gamma ^{\left( 1 \right)}}} & {{\beta ^{\left( 1 \right)}}{\mu ^{\left( 1 \right)}}{\chi ^{\left( 1 \right)}}}\\ { {\text{-}} {\alpha ^{\left( 1 \right)}}{\mu ^{\left( 1 \right)}}{\chi ^{\left( 1 \right)}}} & {{\text{-}} 2{\beta ^{\left( 1 \right)}}{\mu ^{\left( 1 \right)}}k{\nu ^{\left( 1 \right)}}} \end{array}} \right] {\text{,}} $$ (9)
    $$ {\boldsymbol{K}}_u^{\left( 1 \right)} {\text{=}} \left[ {\begin{array}{*{20}{c}} {{{\rm e}^{ {\text{-}} {\gamma ^{\left( 1 \right)}}\left( {z {\text{-}} {z^{\left( 0 \right)}}} \right)}}} & 0\\ 0 & {{{\rm e}^{{\text{-}} {\nu ^{\left( 1 \right)}}\left( {z {\text{-}} {z^{\left( 0 \right)}}} \right)}}} \end{array}} \right]{\text{,}} $$ (10)
    $$ {\nu ^{\left( 1 \right)}} {\text{=}} \sqrt {{k^2}{\text{-}} \frac{\omega }{{{{\left(\, {{\beta ^{\left( 1 \right)}}} \right)}^2}}}} {\text{,}}{\gamma ^{\left( 1 \right)}} {\text{=}} \sqrt {{k^2} {\text{-}} \frac{\omega }{{{{\left( {{\alpha ^{\left( 1 \right)}}} \right)}^2}}}}{\text{,}}(R( {{\nu ^{\left( 1 \right)}}} ){\text{>}}0 {\text{,}}R( {{\gamma ^{\left( 1 \right)}}}) {\text{>}} 0) {\text{,}} $$ (11)

    带入式(5)中第一个子式进行计算,可得瑞雷函数的表达式为

    $$ \begin{align} R\left( k \right) &{\text{=}} 4{k^2}\sqrt {\displaystyle\frac{{{{({\alpha ^{\left( 1 \right)}})}^2}{k^2} {\text{-}} {\omega ^2}}}{{{{({\alpha ^{\left( 1 \right)}})}^2}}}} {(\,{\beta ^{\left( 1 \right)}})^4}\sqrt {\displaystyle\frac{{{{(\,{\beta ^{\left( 1 \right)}})}^2}{k^2} {\text{-}} {\omega ^2}}}{{{{(\,{\beta ^{\left( 1 \right)}})}^2}}}} {\text{-}} 4{({\beta ^{\left( 1 \right)}})^4}{k^4}{\text{+}}\\ & \ 4{(\,{\beta ^{\left( 1 \right)}})^2}{k^2}{\omega ^2} {\text{-}} {\omega ^4} {\text{=}} {\left( {1 {\text{-}} 2x} \right)^2} {\text{-}} 4x\sqrt {m {\text{-}} x} \sqrt {1 {\text{-}} x}\;\;{\text{,}} \end{align} $$ (12)

    式中,$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)}$求逆矩阵的行列式的零点可以通过求解自由表面瑞雷函数的零点来获得。

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

    图  2  f=2 Hz时波速递增的三层地壳模型下位移三分量的实部、虚部图(位移U为归一化后的结果)
    图(a),(b)分别为径向位移Ur的实部和虚部;图(c),(d)分别为切向位移Uf的实部和虚部;图(e),(f)分别为垂向位移Uz的实部和虚部
    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 (Ur),respectively;Figs. (c) and (d) are the real part and image part of the tangential displacement (Uf),respectively;Figs. (e) and (f) are the real part and image part of the vertical displacement (Uz),respectively

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

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

    图  3  波速递增的三层地壳模型在不同频率下kekc的关系(位移Ur为归一化后的结果)
    Figure  3.  The relationship between kc and ke at different frequencies in the three-layer crust model without low-velocity-layers with the normalized displacement Ur

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

    图  4  f=5 Hz时含低速层的五层地壳模型下kekc的关系(位移Ur为归一化后的结果)
    Figure  4.  The relationship between ke and kc in the five-layer crust model with low-velocity-layers while the displace-ment has been normalized when f=5 Hz

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

    图  5  利用求零点法和经验公式法分别获得的理论地震图径向位移Ur (a)、切向位移Uf (b)和垂向位移Uz (c)三分量归一化图
    Figure  5.  Seismogram calculated by for-zero method and empirical-formula method respectively to obtain the radial direction displacement Ur,the tangential displacement Uf,and the vertical displace-ment Uz,which are normalized

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

    图  6  波速递增的三层地壳模型下位移在不同震源与接收点相对深度h时的变化图(位移Ur为归一化后的结果)
    Figure  6.  Relative depth between source and receiver in the three-layer crust model without low-velocity-layers while the displacements Ur have been normalized

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

  • 图  1   不同频率下反透射系数的零点与位移的关系(位移Ur为归一化后的结果)

    图(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 Ur

    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

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

    图(a),(b)分别为径向位移Ur的实部和虚部;图(c),(d)分别为切向位移Uf的实部和虚部;图(e),(f)分别为垂向位移Uz的实部和虚部

    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 (Ur),respectively;Figs. (c) and (d) are the real part and image part of the tangential displacement (Uf),respectively;Figs. (e) and (f) are the real part and image part of the vertical displacement (Uz),respectively

    图  3   波速递增的三层地壳模型在不同频率下kekc的关系(位移Ur为归一化后的结果)

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

    图  4   f=5 Hz时含低速层的五层地壳模型下kekc的关系(位移Ur为归一化后的结果)

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

    图  5   利用求零点法和经验公式法分别获得的理论地震图径向位移Ur (a)、切向位移Uf (b)和垂向位移Uz (c)三分量归一化图

    Figure  5.   Seismogram calculated by for-zero method and empirical-formula method respectively to obtain the radial direction displacement Ur,the tangential displacement Uf,and the vertical displace-ment Uz,which are normalized

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

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

    表  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

    表  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
  • 何耀锋,陈蔚天,陈晓非. 2006. 利用广义反射-透射系数方法求解含低速层水平层状介质模型中面波频散曲线问题[J]. 地球物理学报,49(4):1074–1081. doi: 10.3321/j.issn:0001-5733.2006.04.020

    He Y F,Chen W T,Chen X F. 2006. Normal mode computation by the generalized reflection-transmission coefficient method in planar layered half space[J]. Chinese Journal of Geophysics,49(4):1074–1081 (in Chinese).

    张海明,陈晓非,张似洪. 2001. 峰谷平均法及其在计算浅源合成地震图中的应用[J]. 地球物理学报,44(6):805–813. doi: 10.3321/j.issn:0001-5733.2001.06.010

    Zhang H M,Chen X F,Zhang S H. 2001. Peak-trough averaging method and its application to calculation of synthetic seismograms with shallow focuses[J]. Chinese Journal of Geophysics,44(6):805–813 (in Chinese).

    Apsel R J,Luco J E. 1983. On the Green’s functions for a layered half-space. Part II[J]. Bull Seismol Soc Am,73(4):931–951.

    Capdeville Y,Chaljub E,Montagner J P. 2003. Coupling the spectral element method with a modal solution for elastic wave propagation in global earth models[J]. Geophys J Int,152(1):34–67. doi: 10.1046/j.1365-246X.2003.01808.x

    Chen X F. 1993. A systematic and efficient method of computing normal modes for multilayered half-space[J]. Geophys J Int,115(2):391–409. doi: 10.1111/gji.1993.115.issue-2

    Chen X F. 1999. Seismogram synthesis in multi-layered half-space Part I. Theoretical formulations[J]. Earthquake Research in China,13(2):149–174.

    Fuchs K,Müller G. 1971. Computation of synthetic seismograms with the reflectivity method and comparison with observations[J]. Geophys J Int,23(4):417–433. doi: 10.1111/j.1365-246X.1971.tb01834.x

    Graves R W. 1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences[J]. Bull Seismol Soc Am,86(4):1091–1106.

    Helmberger D V. 1973. On the structure of the low velocity zone[J]. Geophys J Int,34(2):251–263. doi: 10.1111/j.1365-246X.1973.tb02395.x

    Helmberger D V. 1974. Generalized ray theory for shear dislocations[J]. Bull Seismol Soc Am,64(1):45–64.

    Hisada Y. 1994. An efficient method for computing Green’s functions for a layered half-space with sources and receivers at close depths[J]. Bull Seismol Soc Am,84(5):1457–1472.

    Hisada Y. 1995. An efficient method for computing Green’s functions for a layered half-space with sources and receivers at close depths (Part 2)[J]. Bull Seismol Soc Am,85(4):1080–1093.

    Johnson L R. 1974. Green’s function for Lamb’s problem[J]. Geophys J Int,37(1):99–131. doi: 10.1111/j.1365-246X.1974.tb02446.x

    Kennett B L N,Kerry N J. 1979. Seismic waves in a stratified half space[J]. Geophys J Int,57(3):557–583. doi: 10.1111/gji.1979.57.issue-3

    Komatitsch D,Martin R,Tromp J,Taylor M A,Wingate B A. 2001. Wave propagation in 2-D elastic media using a spectral element method with triangles and quadrangles[J]. J Comput Acoust,9(2):703–718. doi: 10.1142/S0218396X01000796

    Komatitsch D,Tromp J. 2003. A perfectly matched layer absorbing boundary condition for the second-order seismic wave equa-tion[J]. Geophys J Int,154(1):146–153. doi: 10.1046/j.1365-246X.2003.01950.x

    Lamb H. 1904. On the propagation of tremors over the surface of an elastic solid[J]. Philos Trans Roy Soc A,203(359/371):1–42.

    Lapwood E R. 1948. Convection of a fluid in a porous medium[J]. Math Proc Cambridge Philos Soc,44(4):508–521. doi: 10.1017/S030500410002452X

    Levander A R. 1988. Fourth-order finite-difference P-SV seismograms[J]. Geophysics,53(11):1425–1436. doi: 10.1190/1.1442422

    Liu T S,Feng X,Zhang H M. 2016. On Rayleigh wave in half-space:An asymptotic approach to study the Rayleigh function and its relation to the Rayleigh wave[J]. Geophys J Int,206(2):1179–1193. doi: 10.1093/gji/ggw189

    Luco J E,Apsel R J. 1983. On the Green’s functions for a layered half-space. Part I[J]. Bull Seismol Soc Am,73(4):909–929.

    Madariaga R. 1976. Dynamics of an expanding circular fault[J]. Bull Seismol Soc Am,66(3):639–666.

    Pekeris C L. 1955. The seismic surface pulse[J]. Proc Natl Acad Sci USA,41(7):469–480. doi: 10.1073/pnas.41.7.469

    Priolo E, Seriani G. 1991. A numerical investigation of Chebyshev spectral element method for acoustic wave propa-gation[C]//Proceedings of the 13th IMACS World Congress on Computational and Applied Mathematics. Dublin: Criterion Press: 551–556.

    Tromp J,Komattisch D,Liu Q. 2008. Spectral-element and adjoint methods in seismology[J]. Commun Comput Phys,3(1):1–32.

    Yao Z X,Harkrider D G. 1983. A generalized reflection-transmission coefficient matrix and discrete wavenumber method for syn-thetic seismograms[J]. Bull Seismol Soc Am,73(6A):1685–1699.

    Zhang H M,Chen X F. 2001. Self-adaptive filon’s integration method and its application to computing synthetic seismograms[J]. Chinese Phys Lett,18(3):313–315. doi: 10.1088/0256-307X/18/3/301

    Zhang H M,Chen X F,Chang S. 2003. An efficient numerical method for computing synthetic seismograms for a layered half-space with sources and receivers at close or same depths[J]. Pure Appl Geophys,160(3/4):467–486.

    Zhang W,Chen X F. 2006. Traction image method for irregular free surface boundaries in finite difference seismic wave simul-ation[J]. Geophys J Int,167(1):337–353. doi: 10.1111/gji.2006.167.issue-1

图(6)  /  表(2)
计量
  • 文章访问数:  1727
  • HTML全文浏览量:  937
  • PDF下载量:  55
  • 被引次数: 0
出版历程
  • 收稿日期:  2018-01-14
  • 修回日期:  2018-07-08
  • 网络出版日期:  2018-11-04
  • 发布日期:  2018-10-31

目录

/

返回文章
返回