引言

广义而言，一切由地球三维非均匀体扰动所引起地震波的变化均可称为地震波散射；狭义上，地震波散射现象通常是指摒除广义概念中可以用几何光学（射线）理论表征的由大尺度非均匀体扰动引起的地震波旅行时和振幅变化的部分，由剩余地球三维非均匀体扰动所引起的地震波场畸变现象（吴如山等，1993；李灿苹等，2005）。而实际地下地质构造十分复杂，各种尺度的非均质体共同存在，地面检波器接收到的往往是多波互相干涉的复杂地震波场信息，地震散射波作为地震波的一种，它是由入射波与地下的非均匀体相互作用而产生的波，因此地面接收到的散射波中携带了大量源于散射体的构造和岩性信息，具有对地下缝、洞等小尺度非均质体进行精细刻画的潜力。而对裂缝-孔洞型等非均匀体散射体进行散射波场正演模拟，分析此类介质的散射波场特征是利用散射波分析研究地下散射体的基础。

对散射波场的正演模拟是识别散射波的基础。Wu和Huang （1992）采用相位屏算子，模拟了二维垂直变背景情况下的散射场；Wu等（1995）引入De Wolf近似和相屏算子，计算了三维常背景情况下的背向散射场；符力耘等（1998）利用扰动理论，建立了用均匀介质格林函数作为基本解的体积分方程，给出了配置法求解体积分方程的数值方法；田丽花（2007）采用相位移法在f-k域对地震散射波场进行了波动方程正演模拟；刘铁华（2012）设计了一种基于微扰论的f-k域积分法，在散射场的二次震源和空间能量衰减处理两方面进行了改进，提高了计算精度和效率；Eaton和David （1999）采用波恩（Born）近似和射线理论近似方法，计算了格林函数的三维弹性波弱散射场；奚先（2018，2020）采用卷积神经网络（convolution neuralnetwork，缩写为CNN）对散射波进行识别，大致识别出测试模型中散射点的准确位置及复杂偏移剖面中的各散射体的位置。

上述研究主要是对一次散射波进行模拟，本文则推导多级散射波和各级散射波微分方程，并基于此采用高阶有限差分对能够产生多级散射波的散射体模型进行正演模拟，分析多级散射波的地震响应特征，以推断地下散射体的分布情况和性质。

1.   基础理论

根据介质分解理论，将地下介质分解为背景介质和扰动介质（图1）。震源激发产生的入射波在背景介质中的波场称为背景波场，背景介质所对应的速度称为背景速度$ {c_0} $。入射波与扰动介质相互作用所产生的波场称为散射波场$ {u_{\rm{s}}} （ {r， \omega } ） $，散射波场可以看作是观测波场与背景波场之差（符力耘等，1998；雷蕾等，2011），即$ {u_{\rm{s}}} （ {r， \omega } ） ＝ u （ {r， \omega } ） － {u_0} （ {r， \omega } ） $。

图  1  地震散射波产生原理示意图

Figure  1.  Schematic diagram of generating seismic scattered wave

下载: 全尺寸图片 幻灯片

波恩近似理论。波恩近似利用波场的振幅，当介质产生的散射波场远小于背景波场时，即$ {u_{\rm{s}}} （ {r， \omega } ） \ll u （ {r， \omega } ） $，用背景波场代替总波场将方程线性化，波恩近似适用于弱散射当围岩和散射体的速度差与围岩速度的比值$ \Delta V ＜ 15\% $的小扰动情况，或传播距离短的情况（尹军杰等，2005）。

地下介质可以划分为背景介质和扰动介质，在研究散射特征中，扰动介质相较于周围介质其体积较小，地震波在穿越背景介质遇到小扰动体时，散射体会作为一次源激发散射，相应称其为散射波。散射波响应过程如图2a所示，当散射体仅为一个时，称之为一级散射，其中S₁表示一次散射波。当地下介质散射体（扰动体）增多后，不同位置的散射体会作为二次源激发其它一次散射源产生的散射波，如图2b所示，散射体$ {V_1} $ 激发的一次散射波在散射体$ {V_2} $ 处发生二次散射（Alkhalifah，2015，2016），称为二级散射；而在散射体$ {V_2} $ 产生的二次散射传播到散射体$ {V_1} $ （也可以是散射体$ {V_3} $，为了方便描述未展示）处再次产生三次散射，称之为三级散射。由此可知，当地下存在大量散射体时，检波器最终接收到的散射波信息数量庞大，非常复杂，是地震波在各个散射体之间错综复杂的传播和被激发所导致的。因此，为了研究这些散射波的特征，将其进行分级处理，散射波对应的级数即为入射波经过几次散射体后所产生的。

图  2  多级散射原理示意图

（a） 一级散射；（b） 二级散射；（c） 三级散射

Figure  2.  Schematic diagram of multi-order scattering principle

（a） The first-order scattering；（b） The second-order scattering；（c） The third-order scattering

下载: 全尺寸图片 幻灯片

2.   地震多级散射波方程推导

2.1   基于波恩近似的散射波场分级

李普曼-施温格（Lippmann-Swinger）方程（Wu，Huang，1992；Wu et al，1995）的总波场计算公式为：

式中：$ \omega $为角频率，${k_0} ＝ {\omega /{C_0}}$为背景介质中波数，$ r ＝ （ {x， {\textit{z}}} ） $为散射点位置，$ {r_1} ＝ （ {{x_1}， {{\textit{z}}_1}} ） $为积分空间的任一点源，$ G $为背景介质中的格林函数，$ \varepsilon （ {{r_1}} ） $为背景速度上叠加的一个扰动量。当${u_{\rm{s}}} （ {r， \omega } ） {\text{≤}} u （ {r， \omega } ） $时，上式右边被积函数中的u用$ {u_0} $代替，则波恩一级近似表达式为（Born，Wolf，1999）

若用$ {u_1} $代替右边被积函数中$ u $，则可得到波恩二级近似表达式

将式（2）带入式（3）可得：

式中：等式右端第二项积分是关于散射体的一级散射，第三项是关于散射体的二级散射，依次进行替换迭代，每项与前一项的迭代关系式为

为方便表示，多级散射的波恩近似表达式则为：

即

则有

式中，上标1，2，···，N为散射波场级数。

2.2   多级散射波方程推导

在均匀介质中，二维标量波方程（无震源项）为

式中：$C （ r ） $为波在均匀介质中的传播速度，${\nabla ^2} ＝ \dfrac{{{\partial ^2}}}{{\partial {x^2}}} ＋ \dfrac{{{\partial ^2}}}{{\partial {{\textit{z}}^2}}}$为拉普拉斯算子，$ t $为时间。为方便推导，将式（9）转换至频率域，表示为

背景波场满足标量波方程

将式（8）和式（11）带入式（10），可得到频率域多级散射波方程

随后，根据介质扰动原理及波恩近似理论，对多级散射波方程进行分级处理。

式（8）中取${u_{\rm{s}}} （ {r， \omega } ） = u_{\rm{s}}^1 （ {r， \omega } ） $，将其与式（11）一同代入式（10）后，结合波恩弱散射近似，可得到一级散射波方程

同理，取${u_{\rm{s}}} （ {r， \omega } ） ＝ u_{\rm{s}}^1 （ {r， \omega } ） ＋ u_{\rm{s}}^2 （ {r， \omega } ） $，将其与式（11）和式（13）一同代入式（10）后，结合波恩弱散射近似，可得到二级散射波方程

取${u_{\rm{s}}} （ {r， \omega } ） ＝ {u_{\rm{s}}^1} （ {r， \omega } ） ＋ {u_{\rm{s}}^2} （ {r， \omega } ） ＋ {u_{\rm{s}}^3} （ {r， \omega } ） $，将其与式（11）、式（13）及式（14）代入式（10）后，结合波恩弱散射近似，整理得到三级散射波方程

以此类推，可以得到，每一级独立的散射波方程

式中，$ N （ {N {\text{≥}} 1} ） $表示散射波的级数，当$ N ＝ 1 $时，等式右端为背景场。

将式（16）转换到时间域，得到时间域各级散射波方程：

式中等号右端整体可视为虚震源，其中二阶导数项${{{\partial ^2}{u^{N － 1}}}}/{{\partial {t^2}}}$无法直接计算获得，本文以空间导数来替代时间导数项。

2.3   多级散射波方程的有限差分离散

采用规则网格有限差分法进行离散，时间二阶导数的离散格式为

式中，$ u $为压力场，$ \Delta t $为时间步长。空间二阶导数的空间2N阶离散格式为

式中, $ {\omega _n} $表示二阶导数的$ 2N $阶精度差分系数（赵茂强，2010；刘庆敏，2007）。规则网格有限差分方程的稳定性条件为（梁展源，2016；李青阳等，2018）

式中，$ C （ {x， {\textit{z}}} ） $为速度，$ \Delta t $为时间步长，$\Delta x 和 \Delta {\textit{z}}$为空间步长，$ {\omega _n} $表示差分系数。最后，加入完美匹配层（perfectly matched layer，缩写为PML）边界条件（何燕，2008；吴国忱等，2014，2020），对多级散射波方程采用有限差分法进行正演模拟。

3.   数值模拟结果及分析

3.1   两点强散射体模型

为验证多级和各级散射波方程的正演模拟精度，设计一简单两点散射模型（图3），模型大小为2 000 m×2 000 m，背景速度为2 500 m/s，在深度1 000 m和1 400 m处分别设置一个点散射体，散射体的速度为2 000 m/s，纵、横向网格间距均为10 m，震源采用主频为25 Hz的雷克子波，炮点位于地表1 000 m处，采样时间点数为2 000，采样间隔为0.8 ms，全排列接收。

图  3  两点强散射体模型

Figure  3.  A model with two strong scatterers

下载: 全尺寸图片 幻灯片

两点强散射模型对应的三个时刻的多级散射波波场和每一级独立的散射波波场如图4所示。图4a中二级散射波和三级散射波的能量弱于一级散射波，为方便展示，对t＝960 ms和t＝1200 ms的波场快照添加增益，因此整体未加色标，图4b−d可直观地看出各级散射波方程所模拟的波场快照仅包含自身级数的散射波，不包含其它级的散射干扰。各级散射波的传播均符合惠更斯-菲涅尔（Huygens-Fresnel）原理，且随着散射波的级数增加，其能量衰减。

图  4  不同级数的散射波波场快照

（a） 多级散射波；（b） 一级散射波；（c） 二级散射波；（d） 三级散射波

Figure  4.  Snapshots of scattered wavefields at different orders

（a） The multi-order scattering；（b）The first-order scattering；（c） The second-order scatterings；（d） The third-order scattering

下载: 全尺寸图片 幻灯片

为讨论各级散射波方程稳定性，抽取各级散射波第100道记录与参考地震记录进行对比，结果如图5所示。图中参考记录为观测记录和背景速度正演模拟记录之差，多级散射记录为本文方程正演所得记录，可以看出：两条记录曲线在相位和振幅上完全吻合，验证了多级散射波方程的准确性。各级散射记录与参考记录第100道的波形对比如图6a所示，其中，t为0.832—0.960 s时，参考记录中的一级散射信息与一级散射记录的相位信息一致，振幅信息存在较小的误差，该误差主要是速度模型中扰动体与围岩的速度差未满足波恩弱散射近似（图6b）；而t为1.240—1.360 s时的误差与一级散射信息的误差来源相同（图6c）。对扰动体与围岩的速度差需满足波恩弱散射近似的速度模型做出如下分析。

图  5  多级散射记录与参考记录第100道的波形对比

Figure  5.  Comparison of the 100th record of the multi-order scattering records with that of the reference records

下载: 全尺寸图片 幻灯片

图  6  各级散射记录与参考记录第100道的波形对比

Figure  6.  Comparison of the 100th record of each scattering records and the reference records

（a） t＝0—1.6 s；（b） t＝0.832—0.960 s；（c） t＝1.240—1.360 s

下载: 全尺寸图片 幻灯片

3.2   两点弱散射体模型

设计两点弱散射模型，设散射体的速度为2 400 m/s，其余参数均与两点强散射模型相同。

两点弱散射模型如图7a所示，其对应的多级散射记录与参考记录第100道的波形对比如图7b所示，由该图可见：两条记录曲线完全吻合，验证了多级散射波方程的准确性。两点弱散射模型对应的各级散射记录与参考记录第100道的波形对比如图8所示，当t为0.848—0.944 s时，参考记录中的一级散射信息与一级散射记录的相位和振幅信息基本一致，验证一级散射波的准确性（图8b）；t为1.296—1.318 s时，参考记录中的二级散射信息与二级散射记录的相位和振幅信息基本一致，验证二级散射波的准确性（图8c）。由于散射的波的能量逐级递减，三级散射波的能量十分微弱，在单道记录上振幅信息不明显，此处不作讨论。

图  7  两点弱散射体模型（a）及多级散射记录与参考记录第 100 道的波形对比（b）

Figure  7.  A model with two weak scatterers （a） and comparison of the 100th record of the multi-orderscattering records and with that of the reference records （b）

下载: 全尺寸图片 幻灯片

图  8  各级散射记录与参考记录第100道的波形对比

Figure  8.  Comparison of the 100th record of each scattering records and the reference records

（a） t＝0—1.6 s；（b） t＝0.848—0.944 s；（c） t＝1.296—1.318 s

下载: 全尺寸图片 幻灯片

3.3   复杂散射体模型

为了验证本文方法对复杂模型的适用性和稳定性，设计一个同时包含不同倾斜角度裂缝和不同尺度孔洞的复杂散射体模型，如图9所示。模型大小2 000 m×2 000 m，背景速度为2 500 m/s，三条裂缝从左到右的倾角依次为90°，0°和45°，速度均为2 300 m/s，裂缝下面四个不同尺度的溶洞速度为2 400 m/s，纵、横向网格间距均为10 m，震源采用主频为25 Hz的雷克子波，炮点位于地表1 000 m处，采样时间点数为2 000，采样间隔为0.8 ms，全排列接收。

图  9  复杂散射体模型

Figure  9.  A model with complex scatterers

下载: 全尺寸图片 幻灯片

复杂散射体模型对应的多级散射波地震记录如图10所示，由该图可见：多级散射记录无直达波，散射波呈现双曲形态，且多个散射体同时存在且距离相近时，各级散射波会产生干涉；同一尺度不同倾角的情况下，横向裂缝的散射能量强于纵向裂缝；当孔洞的尺度小于1/4波长时，炮记录表现为单条双曲线，随着孔洞散射体尺度的增大，其顶底对应的两条双曲线能量区分愈加明显，直至完全分离。

图  10  多级散射记录

Figure  10.  Multi-order scattering records

下载: 全尺寸图片 幻灯片

分别抽取复杂散射体模型对应的多级散射波地震记录与参考地震记录的第100道记录进行对比，如图11所示，图中两条曲线在振幅和相位上完全吻合，表明多级散射波方程对复杂模型具有适用性，并验证了多级散射波方程的准确性。

图  11  多级散射记录与参考记录的第100道波形对比

Figure  11.  Comparison of the 100th record of the multi-order scattering records with that of and the reference records

下载: 全尺寸图片 幻灯片

4.   结论

针对非均匀散射体中地震波传播正演模拟问题，本文推导了多级和各级散射波方程。通过算法分析和模型测试得到以下认识：

1） 本文推导了多级和各级散射波方程，通过有限差分的正演数值模拟，抽取其中心道与参考记录中心道进行对比，结果显示，多级散射单道记录与参考单道记录完全吻合，验证了多级散射波方程的准确性，各级散射记录在波恩弱散射近似范围内，其单道记录与参考单道记录的振幅及相位基本吻合，验证了在波恩弱散射近似下各级散射方程的准确性。

2） 各级散射波的传播符合惠更斯-菲涅尔原理，且随着散射波级数的增加，能量随之衰减。

3） 在散射能量方面，尺度相同的情况下，横向裂缝比纵向裂缝能量更强，散射体与围岩的速度差异越大，其能量也越大。

本文基于标量波方程的散射波方程推导思路理论上能够推广至各向异性介质，但具体推导内容仍需进一步分析；本文采用的背景模型基本为半空间模型，而对于复杂背景介质下的孔洞、缝洞的散射波识别、提取仍需进一步研究。

审稿专家提出有益的修改意见，作者在此表示感谢。