地震波正演数值模拟的计算区域大小有限, 若对计算区域的边界不做任何处理, 地震波到达计算区域边界时会产生数值噪声, 对正常传播的地震波造成干扰。为逼近无限大空间的地震波传播, 需要在有限区域的边界设计数值边界条件, 以吸收边界产生的数值噪声。

早期地震波正演模拟多采用基于旁轴近似的单程波方程吸收条件^[1-2], 该吸收条件仅对接近垂直入射至边界的地震波具有较好的衰减, 对大入射角地震波的衰减效果较差。另一种被广泛应用的吸收边界条件是海绵吸收边界条件, 该边界条件在正常运算区域外部加了一个衰减层, 使得入射至边界衰减层的地震波振幅被逐渐衰减而不产生明显反射^[3]。海绵吸收边界的吸收效果取决于衰减系数和边界宽度等参数, 为了达到较好的衰减效果人们通常要将衰减边界设置得足够宽, 以满足对地震波振幅衰减程度的要求, 这会大幅增加计算量和对计算机内存的占用^[4]。

完全匹配层吸收边界条件(PML)最早由BERENGER^[5]在研究电磁场有限域波动问题时提出, 相比传统的吸收边界条件, PML的吸收效果更好^[6-7]。随后该吸收边界被推广到声波和弹性波的数值模拟^[8-13]。常规PML吸收边界按照不同的实现方式, 可以分为分裂格式PML(SPML)和非分裂格式PML(NPML)^[14]。NPML不需要分裂波场分量, 编程实现较容易, 但实现过程中包含大量卷积, 计算过程复杂^[15-16]。为了提高NPML计算效率和简化运算过程, 前人提出了非卷积形式的NPML吸收边界^[17]。SPML通过分裂波场的传播分量来实现, 计算方程简洁, 计算量相对NPML更小, 所以目前很多正演模拟仍然采用SPML^[18-19]。

然而, SPML边界条件对于切向入射到边界的地震波吸收效果差, 当入射波频率较低时还会在边界附近产生虚假反射^[20]。针对上述问题, 前人提出一种不分裂褶积型PML边界条件(CPML)^[21-22], 其在声波和弹性波数值模拟中均取得了较好的应用效果^[23-24]。但是, CPML在推导中使用了递归卷积的方法, 导致其应用受到较多限制, 如无法适应空间网格的灵活变化, 难以用于高阶的时间离散格式等^[25]。与CPML相比, SPML限制条件少, 可以适用于各种正演模拟方法。可见, 改进SPML以提高边界的吸收效果对于地震波数值模拟在地震学不同领域的应用具有重要意义。

针对SPML对切向入射到边界的地震波吸收效果差的问题, 本文将SPML边界与海绵吸收边界进行空间衔接, 构建内部为SPML, 外部为海绵吸收边界的复合吸收边界, 使SPML未吸收完全的地震波得到二次衰减。同时针对常规海绵边界对大角度入射波吸收不足的问题, 改进了其衰减函数, 以进一步改善复合吸收边界的吸收效果。

1 SPML边界吸收能力的理论分析

伪谱法弹性波正演模拟使用有限差分求时间导数, 通过傅里叶变换在波数域求空间导数, 具有节省空间离散样点数和不产生空间数值频散的优点^[26-28]。

COLLINO等^[29]通过复坐标伸展变换推导出了各向同性介质一阶应力-速度方程的完全匹配层吸收边界条件, 其应用在伪谱法数值模拟中需转化到时间-波数域, 具体形式如下:

$ \left\{ \begin{array}{l} {\rm{i}}{k_x}{{\tilde \sigma }_{xx}} = \rho \left[ {\frac{\partial }{{\partial t}} + d\left( x \right)} \right]\tilde V_x^ \bot \\ {\rm{i}}{k_y}{{\tilde \sigma }_{xy}} = \rho \left[ {\frac{\partial }{{\partial t}} + d\left( y \right)} \right]\tilde V_x^\parallel \\ {\rm{i}}{k_x}{{\tilde \sigma }_{xy}} = \rho \left[ {\frac{\partial }{{\partial t}} + d\left( x \right)} \right]\tilde V_y^ \bot \\ {\rm{i}}{k_y}{{\tilde \sigma }_{yy}} = \rho \left[ {\frac{\partial }{{\partial t}} + d\left( y \right)} \right]\tilde V_y^\parallel \end{array} \right. $

(1)

$ \left\{ \begin{array}{l} {\rm{i}}\left[ {\lambda + 2\mu } \right]{k_x}{{\tilde V}_x} = \left[ {\frac{\partial }{{\partial t}} + d\left( x \right)} \right]\tilde \sigma _{xx}^ \bot \\ {\rm{i}}\lambda {k_y}{{\tilde V}_y} = \left[ {\frac{\partial }{{\partial t}} + d\left( x \right)} \right]\tilde \sigma _{xx}^\parallel \\ {\rm{i}}\lambda {k_x}{{\tilde V}_x} = \left[ {\frac{\partial }{{\partial t}} + d\left( x \right)} \right]\tilde \sigma _{yy}^ \bot \\ {\rm{i}}\left[ {\lambda + 2\mu } \right]{k_y}{{\tilde V}_y} = \left[ {\frac{\partial }{{\partial t}} + d\left( y \right)} \right]\tilde \sigma _{yy}^\parallel \\ {\rm{i}}\mu {k_x}{{\tilde V}_y} = \left[ {\frac{\partial }{{\partial t}} + d\left( x \right)} \right]\tilde \sigma _{xy}^ \bot \\ {\rm{i}}\mu {k_y}{{\tilde V}_x} = \left[ {\frac{\partial }{{\partial t}} + d\left( y \right)} \right]\tilde \sigma _{xy}^\parallel \end{array} \right. $

(2)

式中:k_x和k_y分别代表波矢量的x和y方向分量; λ=ρ(v_P²－2v_S²), μ=ρv_S²为拉梅常数, v_P, v_S和ρ分别为介质的纵波速度、横波速度和密度; ${\tilde \sigma _{xx}}$ , ${\tilde \sigma _{yy}}$ , ${\tilde \sigma _{xy}}$ 为波数域的水平正应力、垂直正应力和剪应力; ${\tilde V_x}$ , ${\tilde V_y}$ 分别是波数域中质点振动速度的水平和垂直分量; $\tilde V_x^{\parallel}$ , $\tilde V_y^{\parallel}$ 表示分裂后平行于边界的速度分量; $\tilde \sigma _{xx}^{\parallel}$ , $\tilde \sigma _{xy}^{\parallel}$ , $\tilde \sigma _{yy}^{\parallel}$ 表示分裂后平行于边界的应力分量; $\tilde V_x^ \bot $ , $\tilde V_y^ \bot $ 表示分裂后垂直于边界的速度分量; $\tilde \sigma _{xx}^ \bot $ , $\tilde \sigma _{xy}^ \bot $ , $\tilde \sigma _{yy}^ \bot $ 表示分裂后垂直于边界的应力分量。中间变量满足下列关系^[30]:

$ \left\{ \begin{array}{l} {{\tilde V}_x} = \tilde V_x^\parallel + \tilde V_x^ \bot \\ {{\tilde V}_y} = \tilde V_y^\parallel + \tilde V_y^ \bot \\ {{\tilde \sigma }_{xx}} = \tilde \sigma _{xx}^\parallel + \tilde \sigma _{xx}^ \bot \\ {{\tilde \sigma }_{xy}} = \tilde \sigma _{xy}^\parallel + \tilde \sigma _{xy}^ \bot \\ {{\tilde \sigma }_{yy}} = \tilde \sigma _{yy}^\parallel + \tilde \sigma _{yy}^ \bot \end{array} \right. $

(3)

d(x)和d(y)是笛卡尔坐标系沿x和y方向的衰减函数, 其表达式为^[31-32]:

$ \left\{ \begin{array}{l} d\left( x \right) = {d_{\rm{o}}}{\left( {\frac{x}{\delta }} \right)^2}\\ d\left( y \right) = {d_{\rm{o}}}{\left( {\frac{y}{\delta }} \right)^2} \end{array} \right. $

(4)

$ {d_{\rm{o}}} = \lg \frac{1}{R}\frac{{3{v_{\rm{P}}}}}{{2\delta }} $

(5)

式中:δ为SPML边界的总厚度; R为理论反射系数。

在SPML区域内, 衰减函数d(x)和d(y)控制了对边界入射波的吸收能力, 它们通过复坐标伸展变换被引入到弹性波方程。二维空间的复坐标伸展变换公式为^[33-34]:

$ \left\{ \begin{array}{l} \tilde x\left( x \right) = x - \frac{{\rm{i}}}{\omega }\int_0^x {d\left( x \right){\rm{d}}x} \\ \tilde y\left( y \right) = y - \frac{{\rm{i}}}{\omega }\int_0^y {d\left( y \right){\rm{d}}y} \end{array} \right. $

(6)

SPML对边界入射波的吸收能力可以通过平面波试探解来评价。设波动方程的平面波解为A·exp[－i(kx－ωt)], 其中, A和ω分别是波的振幅和角频率, $\mathit{\boldsymbol{k}} = {k_x}\mathit{\boldsymbol{\hat x}} + {k_y}\mathit{\boldsymbol{\hat y}},$ , $\mathit{\boldsymbol{x}} = x\mathit{\boldsymbol{\hat x}} + y\mathit{\boldsymbol{\hat y}}$ , ${\mathit{\boldsymbol{\hat x}}}$ 和 ${\mathit{\boldsymbol{\hat y}}}$ 分别代表x方向和y方向的单位向量。

平面波在SPML区域沿x方向的衰减形式为^[35]:

$ A{{\rm{e}}^{ - {\rm{i}}\left( {{k_x}\mathit{\boldsymbol{\tilde x}} + {k_y}y - \omega t} \right)}} = A{{\rm{e}}^{ - {\rm{i}}\left( {\mathit{\boldsymbol{kx}} - \omega t} \right)}}{{\rm{e}}^{ - {k_x}/\omega \int_0^x {d\left( x \right){\rm{d}}x} }} $

(7)

其中, ${\mathit{\boldsymbol{\tilde x}}}$ 代表对平面波进行x方向的衰减。由于(7)式的第2个指数项不显含时间, 可以将其记为衰减项D_x, 代表对振幅的衰减, 表达式如下:

$ {D_x} = {{\rm{e}}^{ - {k_x}/\omega \int_0^x {d\left( x \right){\rm{d}}x} }} $

(8)

衰减项D_x在正常运算区域中的值为1, 表示未对地震波进行衰减, 在SPML区域中的值按指数形式衰减。波矢量的x和y方向分量可以写为:

$ \left\{ \begin{array}{l} {k_x} = \frac{{\omega \cos \gamma }}{{{c_{\rm{p}}}}}\\ {k_y} = \frac{{\omega \sin \gamma }}{{{c_{\rm{p}}}}} \end{array} \right. $

(9)

式中:γ是波前面法向与对应边界法向(例如图 1中x轴方向)的夹角; c_p为相速度。将(9)式代入(8)式, D_x的表达式变为:

$ {D_x} = {{\rm{e}}^{ - \cos \gamma /{c_{\rm{p}}}\int_0^x {d\left( x \right){\rm{d}}x} }} $

(10)

当入射波传播方向近似平行于边界时, 波前面的法线方向n与x轴的夹角γ趋于90°(如图 1所示), D_x趋近于1, 表明SPML边界对大角度入射波几乎没有吸收。

衰减项D_x随入射角γ的变化情况如图 2所示。这里用到的地震纵波速度v_P为3000m/s。当地震波入射角γ趋于90°时, 边界所对应的衰减系数项D_x趋近于1, 导致入射到边界的地震波所受衰减很小, 因此SPML边界难以对其有效吸收。

2 复合吸收边界 2.1 复合吸收边界的结构

复合吸收边界的基本思想是在SPML吸收边界之外添加一层海绵吸收边界, 以对SPML边界未完全吸收的地震波进行二次衰减, 对其充分吸收。复合吸收边界的结构如图 3所示。

在图 3中, 整个区域由内至外依次是正常运算区域(白色)、SPML边界(灰色)和海绵吸收边界(阴影)。上、下SPML边界只衰减含有y方向偏导数的波场分量, 左、右SPML边界只衰减含有x方向偏导数的波场分量。当未被SPML边界充分吸收的地震波进入到海绵吸收边界区域时, 海绵吸收边界会进一步对其进行衰减, 实现对边界波场的二次衰减。

2.2 海绵吸收边界及其改进

由于海绵吸收边界的作用是对SPML边界未完全衰减的地震波作进一步衰减, 它应该能够吸收大角度入射波。然而, 常规海绵吸收边界^[2]对于大角度边界入射波的吸收能力也不够强, 因此有必要对常规海绵吸收边界进行改进, 以提高其对大角度边界入射波的衰减效果。

常规海绵吸收边界条件在计算区域外划定边界层, 在数值模拟过程中将质点的速度值和应力值在边界内分别乘以指数型衰减函数, 使得地震波在边界内逐渐衰减而不产生反射, 从而实现对地震波的吸收。常规海绵吸收边界的指数型衰减函数C在第i层的表达式为^[2]:

$ {C_i} = {{\rm{e}}^{ - {{\left[ {\alpha \left( {l - i} \right)} \right]}^2}}}\;\;\;1 \le i \le l $

(11)

式中:i代表海绵吸收边界内的边界层号; l为海绵吸收边界总层数; α为衰减系数。

当以大角度入射至边界的地震波传播方向与海绵边界切向接近平行时, 其主要沿着衰减系数较小的表层传播, 海绵边界的衰减效果有限, 难以在复合边界中起到对SPML边界剩余波场的充分吸收。为使海绵边界在沿边界的切线方向也具备一定的衰减能力, 我们将海绵吸收边界的衰减函数由仅沿y方向的一维变化推广为沿x和y方向的二维变化, 使指数型衰减函数C具有沿边界法线方向和切线方向同时衰减的能力。改进后的衰减函数C′在(i, j)位置(i代表海绵吸收边界内的边界层号, j代表每层吸收边界的网格点号)的表达式为:

$ \begin{array}{*{20}{c}} {{{C'}_{ij}} = {{\rm{e}}^{ - \left[ {\alpha \left( {l - i} \right)} \right]2}} - \beta \left\{ {1 - {{\rm{e}}^{ - \left[ {\alpha \left( {l - i} \right)} \right]2}}} \right\}\left[ {\sin \left( {2{\rm{ \mathit{ π} }} \cdot \frac{j}{T}} \right)} \right]}\\ {1 \le i \le l\;\;\;\;1 \le j \le {N_x}} \end{array} $

(12)

$ T \cdot \min \left( {{d_x},{d_y}} \right) \ge \lambda $

(13)

式中:α和β分别代表沿边界垂向和切向的衰减系数; N_x, N_y代表正演区域沿纵、横向的网格点数。由于地震波横向振幅变化剧烈, 会产生明显的反射噪声, 所以衰减函数沿x轴方向的变化需要足够缓慢。一种设置方法是使周期T与网格间距d_x或d_y的乘积大于地震波的一个波长λ。

图 4显示了海绵吸收边界改进前、后的衰减函数C和C′的分布。y=0是海绵边界与SPML边界的交界。图 4a和图 4b分别对应(11)式和(12)式的衰减函数, 参数N_x=601, d_x=10, d_y=10, T=100。对于l=10的边界, 对应的垂向和切向衰减系数分别为α=0.04, β=0.20。

图 4a中常规海绵边界的衰减函数在y方向上呈指数变化, 而边界内各层的衰减函数值在x方向都为常量。若地震波以大角度入射至海绵边界时(即接近平行于x轴方向), 由于其表层衰减函数值趋于1(等于1表示没有衰减), 致使传播至表层的地震波不能得到快速有效的吸收, 因此常规海绵吸收边界所采用的衰减函数在沿x方向的衰减能力不足。相比之下, 图 4b中改进型海绵边界的衰减函数值除了沿y方向呈指数变化外, 沿x方向也存在高低相间的周期变化, 谷值位置对地震波的衰减能力强于峰值。由于常规海绵吸收边界的衰减能力与峰值的衰减效果相同, 因此改进型海绵边界对大角度入射波的总体衰减能力要强于常规海绵吸收边界。

3 数值测试 3.1 改进型海绵吸收边界衰减效果测试

改进前、后海绵吸收边界对大角度入射波的衰减效果可以通过对比波场和地震记录来观察和分析。为了不产生额外干扰, 我们选择了均匀介质模型, 地震波速度v_P=3300m/s, v_S=1833m/s, 介质密度ρ=2800kg/m³, 模型的大小为6400m×1000m, 离散网格尺寸为10m×10m。数值模拟过程中, 为满足稳定性条件, 选择的时间采样间隔为0.3ms; 震源为P波源, 震源子波采用主频为7Hz的Ricker子波。海绵边界衰减函数的参数α=0.04, β=0.20。复合吸收边界由20层平行于边界的网格排列构成, SPML和海绵边界分别占用内层和外层的10层网格。震源设置在接近复合边界的位置(图 5中五角星处), 以产生大角度边界入射波。

当复合吸收边界使用图 4中所示的两种海绵吸收边界的衰减函数时, 弹性波数值模拟得到的波场快照在边界处将产生明显差异, 如图 5所示。使用常规海绵吸收边界时, 运算区域的边界处会产生明显的噪声波形, 而在改进海绵边界的快照中则较弱(图 5中箭头所示位置)。这说明, 常规海绵吸收边界对大角度边界入射波的吸收效果弱于改进型海绵吸收边界。

3.2 复合吸收边界衰减效果测试

为分析复合吸收边界对大角度入射波的衰减效果, 本节使用复合吸收边界和SPML边界分别计算地震波场, 并对两种情况所获得的地震波场和理论地震记录进行对比分析。我们选择了层状速度模型和Marmousi速度模型分别进行正演模拟和结果分析。

3.2.1 水平层状模型

速度模型中出现较强的速度跳变时, 由低速介质中震源发出的地震波入射至高速介质后的传播方向与边界法向夹角较大。图 6是一个两层速度模型, 其上层和下层的速度分别为3300m/s和1500m/s, 纵横波速度比为1.8, 介质密度ρ=2800kg/m³。模型被离散化为620×240个规则网格, 纵向和横向网格大小均为10m。为满足数值模拟稳定性条件, 时间采样间隔采用0.3ms, 震源子波采用主频为7Hz的Ricker子波。

图 7给出了层状速度模型采用SPML和复合吸收边界得到的水平分量波场快照。由图 7可见, 使用SPML边界的波场快照在边界附近产生了较强的数值噪声(图 7a中箭头所示), 产生这种现象的主要原因是地震波大角度入射时, 边界处的SPML衰减函数失效, 导致地震波未被有效衰减, 在边界处产生了数值噪声。在速度分界面处, 边界数值噪声的振幅并未出现明显变化, 表明这种干扰噪声能够传播至正常运算区域, 这对于所要模拟波场的危害很大。随着波传播时间进一步增加, 这些数值噪声将透过边界, 向正常运算区域内扩散, 并与正常地震波混合, 降低正演模拟结果的精度。

采用复合吸收边界后, 边界处的数值噪声要弱很多, 如图 7b中箭头所示。虽然在边界层内部存在少量数值噪声, 但它们被较好地控制在吸收边界内部, 并未明显出现在正常运算区域, 因而对正常传播的地震波的干扰可以忽略, 有效提高了对大角度入射波的吸收效果。

图 8是在平行于边界的地震测线上(图 6中虚线)采集到的不同边界条件水平分量地震单炮记录。图 8a显示, SPML吸收边界在较大偏移距处产生了较强的虚假波形(椭圆虚线处所示)。这些虚假波形紧贴边界传播, 并有随着偏移距和波传播时间增加而变强的趋势。这种变化趋势与地震波在远偏移距变弱的趋势相反, 导致有效地震信号最终被噪声所淹没。由于边界吸收不足进入正常运算区域的噪声不严格遵守正演模拟控制方程, 因此随着计算时间的增加, 较强的数值噪声会使正演模拟结果变得不稳定^[30]。图 8b中的地震记录并未出现明显的虚假波形(椭圆虚线处所示), 这表明复合边界对大角度入射波进行了有效吸收。这种现象在图 8放大后显得尤为突出(图 9)。可见, 复合吸收边界对大角度入射波的吸收效果要明显优于传统SPML吸收边界, 降低噪声的同时提高了计算的精度。

3.2.2 Marmousi模型

为了进一步测试复合吸收边界在复杂速度模型中的应用效果, 采用如图 10所示的Marmousi速度模型, 纵横波速度比取1.8, 介质密度ρ为2800kg/m³。整个空间离散为1149×375个均匀网格, 网格大小为20m×20m, 数值模拟采用的时间采样间隔为0.3ms, 震源子波采用主频为7Hz的Ricker子波, 炮检点位置如图 10所示。

图 11给出了Marmousi模型采用不同边界条件得到的波场快照(5600ms)。由图 11可见, 当地震波接近切向入射到边界时, 使用SPML边界的波场快照在边界处出现了数值噪声(图 11a中绿色箭头所示), 这些数值噪声向正常运算区域内扩散, 与正常地震波混合, 降低了正演模拟结果的精度。与采用SPML边界的结果不同, 采用复合吸收边界后, 边界处并未产生明显的数值噪声(图 11b中绿色箭头所示)。

吸收边界条件不同所产生的数值噪声的差异也反映在地震记录中。图 12是平行于边界地震测线(图 10中黄色虚线)的水平分量单炮地震记录。由图 12可见, SPML边界地震记录在5.5s之后存在明显的虚假波形(较强部分如图 12a中绿色椭圆处所示), 而复合吸收边界的地震记录中未出现明显的虚假波形(图 12b中绿色椭圆处所示)。

为对比地震记录的细节, 从图 12地震记录中抽取10~15km的地震道, 并以波形的形式显示, 如图 13所示。图 13中的波形对比更加明显地展示了SPML边界地震记录中的边界干扰波与模拟所产生的波场重叠, 导致部分反射波难以辨别。然而, 这种边界干扰波在使用复合吸收边界所得到的地震记录中并不明显。这表明复合吸收边界对于复杂速度模型伪谱法正演模拟产生的边界干扰波也具有较好的吸收效果。

4 结论

本文将SPML边界与海绵吸收边界相结合, 提出了一种适用于伪谱法弹性波数值模拟的复合吸收边界条件, 有效提高了对接近平行于边界传播的大角度边界入射波的吸收能力, 并得到如下结论和认识:①理论分析显示, SPML边界对大角度入射波衰减能力不足; ②将常规海绵吸收边界的衰减函数由一维变化改进为二维变化可以有效吸收平行于边界传播的地震波; ③复合吸收边界继承了SPML对小角度入射波的良好吸收效果, 且对大角度入射波的吸收效果优于SPML边界。

本文提出的复合吸收边界数值实现过程较为方便, 理论上既可以应用于伪谱法正演模拟, 也可以应用于有限差分、有限元等格式的地震波数值模拟中, 在不进行复杂数学推导的情况下提高了对边界干扰波的吸收能力。