1 引言

随着地震波探测技术的高速发展，传统的层状构造地球模型正在经历一场变革，地球被证明从地壳、地幔、到地核处处都呈现出多尺度的横向非均匀性^[1-2].而非均匀性往往意味着更多的地质构造运动或者矿产资源的存在.对勘探地球物理来说，这些不同尺度的横向非均匀性与油气及金属矿产的分布紧密相关^[3-5].对于土木工程来说，考虑传播介质非均匀性对地震波波形的影响，可以更真实地给定结构动力分析的地震动输入波形^[6-7].地球岩石层中的非均匀性在地震图上的表现是P尾波、S尾波以及Lg尾波的存在.高质量的高频数字地震资料的日益增加，使得越来越多的散射理论在地震尾波的衰减测定以及地球内部构造勘探中得到验证与应用^[8-13].Aki等^[14-15]首先用原生波在横向非均匀性介质中传播时的单次回返散射模型来研究区域地震尾波，这种模型适用于台站与震源距离较小的情况，二者可简化为在同一点处重合.Sato^[16]推广了Aki等人的理论，给出了观测点与震源不在同一点上的单次散射模型.随着地震波传播距离的增加，多次散射逐渐占主导作用，单次散射模型不再适用.Gao^[17-19]发展了Aki等人的单次散射理论，研究了多次散射对区域性地震尾波的影响.然而，在Gao^[18]的多次散射公式中低估了高阶散射项(两次以及两次以上)的作用^[20-21].曾健和聂用安^[22-23]在Gao的多次散射理论基础上，提出了震源距不为零的多次散射模型，得到了二维、三维无限介质中地震尾波功率谱的积分表达式.在静态能量传输理论的框架内，Wu^[24-25]研究了地震波的多次散射问题，并分离了固有衰减和散射衰减的影响.20世纪90年代初，Zeng等^[20]在Wu的静态能量传输理论中引入时间项，推导了散射波能量基本方程，给出了高阶散射项的公式表达.与此同时，Hoshiba^[21]基于能量守恒原理，采用MonteCarlo方法得到了高阶散射项的扩展系数.自此以后，在能量传输理论框架内，地震波多次散射的作用逐渐被越来越多的地震学者^[26-29]所认识，并应用于地震波的尾波分析.

然而，为了数学表达上的简化，上述学者的散射模型都是基于地震波各向同性散射假定.虽然，在地震台站距离震源较近时，这种模型可以有效地模拟地震尾波能量密度包络衰减的一致性，但是不能准确模拟地震能量密度包络的展宽现象^[30].为了更真实地模拟地震波的散射过程，越来越多的地震学者开始研究地震波的各向异性散射理论，Sato^[31]首先提出了单次各向异性散射模型，并应用于较小震源距的地震波能量密度包络的计算.之后，基于能量传输理论，Sato^[32-33]建立了二维、三维无限空间中的多次各向异性散射理论，并给出了地震波能量密度的积分方程.Hoshiba^[34]基于能量传输方程，利用MonteCarlo方法求解了地震波各向异性散射的能量传输方程.Gusev和Abubakirov^[35]指出利用反映强前散射过程的散射系数，通过MonteCarlo方法合成地震波的能量密度包络可以模拟地震波包络出现的展宽现象，但其结果在不同能量密度包络之间的衰减一致性的模拟上不是很理想.

为了同时模拟不同震源距S波尾波能量密度包络的衰减一致性和S 波包络的展宽现象，本文在Sato的多次各向异性散射理论框架内，通过在方向性散射系数的球谐函数展开式中引入高阶项来表征地震波不同强度的散射过程.并引入特征时间的概念来定义震源处的辐射脉冲宽度.基于此，把任意给定频率的地震波脉冲的能量谱密度引入到Sato的地震波能量密度积分方程中.然后利用离散波数方法得到了不同散射模式下的S波能量密度的包络曲线.讨论了特征时间对能量密度包络的影响，对比分析了多次前、后散射以及各向同性散射情况下的能量密度包络的差异.最后，从包络宽度与震源距的变化结果出发，研究了不同散射模式对S波能量密度包络展宽的影响.并初步探讨了能量密度包络展宽现象的形成机理.此外，以美国2008 年Wells地震的余震为例，简单讨论了多次前散射模式的有效性与实用性.

2 基本理论 2.1 Sato的各向异性散射理论

Sato^[33]基于能量传输理论，在三维随机无限介质内建立了表征地震波各向异性散射过程的数学模型.该模型假定传播介质中的不均匀体为随机均匀分布的散射点，散射点的随机分布，使得地震波与散射点相互作用产生非相干波场，从而使得地震波散射能量可以叠加.并假定散射模式是轴对称的各向异性散射，即散射强度仅与平面散射角ψ 相关，地震波散射的各向异性特征由方向性散射系数g(ψ)来表征，总的散射系数g₀ 是g(ψ)在立体角方向上积分的平均值

(1)

式中$\hat{q}$= (1，θ，φ)是球坐标系下的单位向量，图 1是球坐标系下多次各向异性散射过程的示意图.

2.2 介质的非均匀性与散射系数的关系

地球介质中随机分布的非均匀性一般用介质特征参数扰动量的自相关函数和谱密度函数来表征.例如，介质中地震波速的非均匀性可表示为

(2)

式中V₀ = 〈V(x)〉是介质中的平均波速，ξ(x)是x处速度的微小扰动，〈ξ(x)＞= 0，〈…〉表示整体均值.ξ(x)是一个关于随机变量x的均匀的、各向同性的随机函数.这样速度扰动的自相关函数可以定义为波速扰动的均值:

(3)

自相关函数做空间傅里叶变换就可以得到速度扰动的谱密度函数.根据物理学中关于方向性散射系数的定义^[2]:针对单个散射体，它是微分散射截面的函数.对于介质内随机分布的多个散射体，它是微分散射截面均值的函数，反映传播介质内的散射场的平均特性.而微分散射截面的均值是介质特征参数扰动量的谱密度函数^{[2， 36]}.因此，在统计意义上，方向性散射系数反映了介质的非均匀性质.作为方向性散射系数的特例，在均匀介质中，g(ψ)是一常数.

2.3 方向性散射系数的展开式

应用各向异性散射理论求解地震波能量密度的关键是方向性散射系数的球谐函数展开式，并且展开项数的多少决定了表征地震波散射强度的完备性.本文首先从球函数的基本概念出发，重新推导了方向性散射系数的球谐函数展开式.在物理学中，球函数常常采用复数形式表达^[37]:

(4)

式中P_n^m(cosθ)是m阶伴随Legendre方程的解，球坐标0≤θ ≤ π，0≤φ ≤2π.

球函数的正交性表示为

(5)

式中模平方

(6)

则球面上任意函数f(θ，φ)的球函数展开式为

(7)

展开系数:

(8)

为了便于利用Sato^[33]的地震波能量密度积分方程，需要对方向性散射系数做无量纲化处理:

(9)

针对轴对称的各向异性散射问题，定义在球坐标下的方向性散射系数g(ψ)仅与平面散射角ψ 有关，ψ =arccos($\hat{r}$″·$\hat{q}$′)，这里$\hat{r}$″ = (1，θ″，φ″)，$\hat{q}$′ =(1，θq′，φq′)，散射角ψ 相当于(7)式中的θ ，并且g(ψ)与另外一个角度无关，所以g(ψ)球函数展开式可简化为

(10)

(11)

根据Legendre多项式的加法定理

(12)

(13)

(14)

(15)

(13) 式与Sato^[33]的(17)式(参见附录中(A1)式)是一致的，但这里的表达较为简洁，数值求解更直观.当n=0时，g(ψ)=1，对应于各向同性散射模式，随着n的增加，g(ψ)表征的散射模式趋于完备，可以准确描述在某一方向散射较强的情况.

2.4 能量密度积分方程的修正

在Sato的散射理论中，由于假定震源激发的地震波脉冲为理想脉冲，导致在计算有限频率带宽时，地震波能量密度包络往往出现Gibbs现象.而实际上震源处辐射波脉冲往往是有一定宽度的，本文利用Sato^[2]在抛物线近似方法中提出的特征时间t_M来表征震源处初始脉冲宽度，并基于Sato的各向异性散射模型，在地震波能量密度积分方程中引入初始地震波脉冲能量谱密度的频域表达^[2]I(ω)，对原有的积分方程进行修正.I(ω)中t_M 的引入，可以有效地避免初始理想脉冲引起的地震波能量密度包络中Gibbs现象的出现，相关的数值验证见第4 节.修正之后的无量纲化地震波能量密度表达式如下:

(16)

(17)

(18)

(19)

(20)

E_S(ρ，τ)是Sato的能量密度积分方程.式(18)中，t_M 是与震源相关的常数，通过与实际台站记录的拟合计算，一般地t_M =0.0~0.5.并且取值与震源距的大小无关，它表征震源处地震波能量脉冲的初始宽度.此处t_M 与Sato^[2]的物理含义是不同的，Sato定义的t_M 是一个随震源距递增的量，表征地震波能量脉冲随着震源距的增加而逐渐展宽.(16)、(17)式中ρ = |ρ| 是无量纲化的震源距，τ 是无量纲化时间，f(ρ，τ;$\hat{q}$)是ρ 处散射点在τ 时刻沿方向$\hat{q}$的无量纲化能量密度，f₀(k，s)的含义及表达式见附录中的(A2)式.无量纲化参数可以由(20)式还原，其中V₀ 是地震波传播的平均速度，b是固有吸收因子，W是震源处的辐射能量.

3 数值计算方法

为了求解地震波能量密度积分方程，按照地震波能量的构成，我们可以把无量纲化能量密度分解为直达波、单次散射项、多次散射项(二次散射以上)三部分:

(21)

式中

单次散射项E¹(ρ，τ)的求解直接利用Sato^[33]推导的积分方程，见附录(A6)式.

关于多次散射项E^M(ρ，τ)，本文采用了Zeng^[38]的离散波数法对积分方程(19)式进行波数域上的离散.其中积分函数f₀(k，s)可以由附录(A2)式得到.

求得E_S(ρ，τ)以后，采用傅里叶变换得到其频域结果E_S(ρ，ω)，乘以初始地震波脉冲能量谱密度的频域表达I(ω)，获得地震波修正的能量密度频域表达.然后按照(16)式做傅里叶逆变换，并通过(20)式还原到实际的地震波能量密度E(r，t).

4 数值计算实例

以指数型自相关函数表征的随机介质为例，本文采用Hoshiba^[34]给出的方向性散射系数表达，来讨论不同散射模式对S波能量密度包络的影响.方向性散射系数表示如下:

(22)

式中ν>-1，当ν>0时，前散射强于后散射，称之为前散射模式;反之ν ＜ 0 时，后散射强于前散射，称之为后散射模式，ν = 0 对应于各向同性散射模式.在图 2中分别取ν=-0.9，-0.7，0，3.0，9.0来表征不同的散射模式.在散射平面上，横轴表示g(ψ)的水平分量g(ψ)cosψ，纵轴表示g(ψ)竖向分量g(ψ)sinψ.其中针对ν= -0.9、-0.7、3.0、9.0的强散射情况，需要在g(ψ)的球函数展开式中包含高阶项，在附录(A3)式中由于高斯型超几何函数的存在，使得高阶的格林函数的计算很不稳定，本文采用了Zhang^[39]给出的计算高斯型超几何函数的数值方法，并通过格林函数的递推公式进行高阶项的计算，可以将方向性散射系数的球函数展开至10阶以上.一般地，n≤8就能满足实际地震波散射强度的需要，这一点将做另文讨论.在本算例中取n=6.

首先，取震源距r=50、100、200、400km;根据S波波速、散射系数和固有吸收系数的一般取值范围^[20]，地震波平均传播速度V₀ 取为3.5km/s;总散射系数g₀ 取为0.01km^-1;固有吸收系数ηi取为0.01km^-1，并表示为ηi=b/V₀;震源辐射能量W取为1.针对图 2中给出的不同散射模式，利用修正的能量密度积分方程合成S 波的能量密度包络曲线，如图 3所示.从计算结果可以看出，当震源距较小时，无论是多次前散射模式、各向同性散射还是后散射模式都能预测到能量密度包络尾波段的同一衰减形式.然而随着震源距的增加，如r=400km 时，各向同性散射、后散射模式合成的能量密度包络尾波段衰减趋势逐渐出现偏差，而前散射模式的结果则能准确呈现不同震源距尾波段的衰减一致性.

其次，基于图 2给出的散射模式，在图 4中讨论了t_M 取值对S 波能量密度包络的影响.不考虑t_M时，即t_M =0.0，合成的地震波能量密度包络在S波到时附近出现高度的振荡现象，以致于S波能量密度包络的峰值不稳定.而当t_M = 0.1 时，则能合成平滑的地震波能量密度包络.因此，在地震波能量密度积分表达中引入t_M ，可以有效地消去能量密度包络的Gibbs现象.

然后，使用相同的输入参数:r，V₀，g₀，η_i，我们还讨论了图 2所示的不同散射模式对S波能量密度包络的影响，从图 5 中给出的地震波能量密度包络曲线可以看出，当震源距较小时，比如r≤100km时，各向同性散射模式与各向异性散射模式的结果十分吻合.但随着震源距的增加，从图 5右侧部分结果可以看到，在S波到达以及到达以后的早期阶段，前散射模式合成的地震波能量密度包络与各向同性散射的结果相比较大，这点与Sato^[33]的结果是一致的，并且随着前散射强度的加强，能量密度包络峰值逐渐增加.对于相同的时间段，如图 5 左侧部分所示，同样随着震源距的增加，各向同性散射模式合成的地震波能量密度与后散射模式之间结果的差别逐渐增大，并且各向同性散射的结果大于后散射的结果.因此，当震源距较大时，如r>100km 时，各向异性散射的结果与各向同性的结果会出现较大差异，尤其是强散射情况，此时，需要考虑地震波各向异性散射过程的影响.

此外，在实际观测的地震图记录中，随着震源距的增加，地震波能量包络会出现逐渐展宽的现象，之前关于这方面的研究一直采用抛物线近似的方法，但该方法完全忽略了后散射模式的影响，不能准确地模拟地震波能量密度包络的尾波部分.本文的多次各向异性散射模型基于能量传输理论建立，可以同时考虑地震波的前、后散射过程，在图 6中给出了不同散射模式下S波能量密度包络的宽度t_P、t_Q 随震源距的变化情况，t_P、t_Q 的物理意义如图 7所示.从图 6的数值结果可以看出，在多次前散射模式下，能量密度脉冲宽度t_P、t_Q 随着震源距的增加依次递增，从而说明多次前散射模式的结果能够反映S波能量密度包络的展宽现象，而其他散射模式则不能说明这一点.由此，S波能量密度包络随震源距的增加而出现的展宽现象可解释为地震波与介质中不均匀体的多次前散射作用所引起.这点与Gusev 和Abubakirov^[35]结论是一致的.

最后，以2008年美国内华达州的Wells地震的一次余震为例，该余震发生于当地时间2008年2月21日，震级4.84，震源深度10.56km.取其中三个台站BMN、K10、K11的竖向地震波速度记录，并利用地震波处理软件SAC2000 扣除记录中台站仪器响应的影响，在8~16Hz频率范围内对记录进行滤波，最后得到对应的地震波竖向速度包络曲线，如图 8中灰色线所示，其中图 8 中的黑色细实线表示灰色线的平滑处理结果.利用本文提出的修正的能量密度积分方程，基于图 2 给出的不同散射模式，由(16)式求得S波能量密度的时空分布，取能量密度的平方根拟合台站记录的速度包络.得到t_M = 0.09，g₀ =0.005，以及ν=9.0时的地震波前散射模式.拟合得到的S 波能量密度的平方根结果如图 8右侧的黑色粗实线所示，图中竖向虚线表示S波到时所在位置，从S波到时算起，对照实际台站记录的平滑处理结果，多次前散射模式可以得到与其一致的能量密度的平方根曲线.因此，与各向同性散射、多次后散射模式相比，多次前散射模式可以拟合高频的S波台站记录.该散射模式为观测地震图的解释提供了理论依据.

5 结论

基于Sato的多次各向异性散射理论，本文在方向性散射系数的球谐函数展开式中引入了高阶项来表征地震波的不同散射强度，在地震波能量密度积分方程中引入表征初始地震波能量密度脉冲宽度的特征时间t_M ，避免了数值求解能量密度包络时出现的Gibbs现象.针对指数型自相关函数表征的随机传播介质，利用离散波数法求解能量密度积分方程，对比分析了不同散射模式对S波能量密度包络曲线的影响，得到了以下几点结论:

(1) 在能量传输理论框架内，利用Hoshiba的方向性散射系数表达生成的前散射、各向同性散射、后散射模式，在震源距较小时，均能预测S波能量密度包络尾波段的同一衰减特性.然而，随着震源距的增加，各向同性散射、后散射模式模拟的尾波包络同一衰减逐渐出现偏差，但前散射模式的结果始终能够给出能量密度尾波段的衰减一致性.

(2) 在震源距较小(r≤100km)时，各向异性散射模式的能量密度包络与各向同性散射结果十分吻合;但随着震源距的增加，与前散射模式相比，前散射强度越强，在S波到达以及到达以后的早期阶段，能量密度包络与各向同性散射的结果相比越大.在同样的时间段，与后散射结果相比，能量密度包络与各向同性散射结果相比越小.所以当震源距较大时，应当充分考虑地震波的各向异性散射过程的影响.

(3) 与多次后散射、各向同性散射模式相比，前散射模式能够重现S波能量密度包络在传播过程中出现的展宽现象，因此，我们可以推断:S 波能量密度包络随着震源距的增加而出现的展宽现象源于地震波与介质中非均匀体的多次前散射作用.

(4) 通过与美国内华达州Wells地震余震的台站记录比较，证明了多次前散射模式可以拟合实际地震台站的记录资料.

综上所述，在能量传输理论框架内，利用修正的能量密度积分方程，基于多次前散射模式，合成的S波能量密度包络不仅能够模拟尾波段包络的衰减一致性，而且能够模拟S波能量密度包络在传播过程中出现的展宽现象.所以，多次前散射模式在实测地震图中S波的解释，以及理论地震图的合成方面，具有良好的应用价值.

附录 A Sato相关公式^{[2， 33]}

下述公式均在球坐标系下建立，参见正文图 1.无量纲化方向性散射系数的球谐函数展开式:

(A1)

针对无量纲化能量密度函数f(ρ，τ;$\hat{q}$)做空间域上的傅里叶变换，之后做时间域上的拉普拉斯变换，然后做球函数展开.在其球函数展开式中，球函数的系数$\overline{\overline{f}}$_l(k，s)可以由(A2)式得到，其中当l=0时，$\overline{\overline{f}}$₀(k，s)为能量密度积分方程中的积分核函数，

(A2)

式中格林函数$\overline{\overline{G}}$_l(k，s)的表达式如下:

(A3)

按照Wigner3-j符号可表示为三个球函数的积分:

(A4)

上述积分式需要满足:m₁ = m+m₂ 和|l-l₂| ≤l₁ |l+l₂| ，否则积分值为零.这种情况下l+l₁ +l₂ 为偶数，(A4)式的平方值为

(A5)

Sato的单次散射能量积分方程:

(A6)

这里是扁平椭球体坐标系变量，1≤v ＜ ∞，0≤u≤π，其他参数说明参见正文.

致谢

本文第一作者在美国学习期间得到了国家留学基金委“建设高水平研究型大学"项目的资助，美国密苏里科技大学土木系以及USGSEarthquakeHazardsProgram 提供了学习交流的机会，两位匿名审稿专家提出了建设性意见，使得本文对问题的讨论更加清楚，作者在此表示衷心的感谢.