地震正演模拟是探索地震波在介质中传播规律的重要手段, 基于波动方程的数值模拟方法是一种地震正演模拟方法^[1]。目前波动方程法主要包括伪谱法^[2]、有限元法^[3]、边界元法^[4]、谱元法^[5]以及有限差分法。有限元法和有限差分方法相比, 二者精度相同时, 有限元法耗时长^[6]。有限差分法因其具有简单灵活、计算效率高以及占用内存小等优势被广泛应用于地震波场数值模拟^[7]。近年来，随着人们对数值模拟结果精度的要求不断提升，传统有限差分数值模拟方法的不足逐步凸显:一是数值计算时所需的网格点数多; 二是具有严重的数值频散^[8]。相较传统的差分格式, 紧致差分格式在相同的计算网格时, 精度和分辨率高, 且稳定性高^[9-10]。紧致差分格式的发展可以追溯到1989年DENNIS等针对Navier-Stokes方程首次提出了空间四阶紧致格式^[11];1992年LELE对Pade格式进行了研究, 提出了求解一阶导数和二阶导数的对称型紧致差分格式^[12], 该格式精度可达到谱方法的精度; 1998年CHU等提出了三点六阶超紧致差分(combined compact difference 6, CCD6)格式用于对流扩散方程^[13];林东等在三点CCD6格式的基础上提出了组合型超紧致差分格式, 其精度最高可达到十阶^[14]。在紧致差分格式中引入迎风机制可有效抑制非物理振荡, 理论上可以提高紧致差分格式的稳定性。1994年傅德薰在紧致差分格式中引入迎风机制^[15];1997年FU等又提出了更适合于多尺度复杂流场计算的五点五阶精度的迎风紧致格式^[16];2002年王书强等在地震波场数值模拟时, 首次将常规紧致差分格式应用于弹性波方程的数值模拟^[17], 截至目前，超紧致差分格式的地震波场数值模拟应用未见文献报道。

本文首先讨论了五点八阶超紧致差分(combined compact difference 8, CCD8)格式, 然后通过引入迎风机制, 在五点八阶超紧致差分格式的基础上, 构造了五点七阶迎风超紧致差分(upwind combined compact difference 7, UCCD7)格式, 并将这两种格式应用于位移场空间导数的求取, 实现了二维均匀介质模型、水平层状介质模型以及Marmousi模型的地震波场数值模拟。

1 超紧致差分格式及其迎风型方法原理

早期的紧致差分格式是基于Hermite多项式构造形成的, 1992年LELE对Hermite公式进行了扩展, 构造了紧致差分格式^[12]:

$ \begin{array}{l} {\beta _1}{{f'}_{i - 2}} + {\beta _2}{{f'}_{i - 1}} + {{f'}_i} + {\beta _2}{{f'}_{i + 1}} + {\beta _1}{{f'}_{i + 2}} = \\ c\frac{{{f_{i + 3}} - {f_{i - 3}}}}{{6h}} + b\frac{{{f_{i + 2}} - {f_{i - 2}}}}{{4h}} + a\frac{{{f_{i + 1}} - {f_{i - 1}}}}{{2h}}\\ {\beta _1}{{f''}_{i - 2}} + {\beta _2}{{f''}_{i - 1}} + {{f''}_i} + {\beta _2}{{f''}_{i + 1}} + {\beta _1}{{f''}_{i + 2}} = \\ c\frac{{{f_{i + 3}} - 2{f_i} + {f_{i - 3}}}}{{9{h^2}}} + b\frac{{{f_{i + 2}} - 2{f_i} + {f_{i - 2}}}}{{4{h^2}}} + \\ a\frac{{{f_{i + 1}} - 2{f_i} + {f_{i - 1}}}}{{{h^2}}} \end{array} $

(1)

式中:h为网格间距; $ a, b, c, {\beta _1}, {\beta _2}$均为差分系数; $ {f_i}$为i节点的函数值; ${{f'}_i} $和$ {{f''}_i}$分别表示i节点的一阶和二阶导数; $ {f_{i + 1}}, {f_{i + 2}}, {f_{i + 3}}, {f_{i - 1}}, {f_{i - 2}}, {f_{i - 3}}$分别表示i节点依次向前3个节点和依次向后3个节点的函数值; $ {{f'}_{i + 1}}, {{f'}_{i + 2}}, {{f'}_{i - 1}}, {{f'}_{i - 2}}$分别表示i节点依次向前两个节点和依次向后两个节点的一阶导数; $ {{f''}_{i + 1}}, {{f''}_{i + 2}}, {{f''}_{i - 1}}, {{f''}_{i - 2}}$分别表示i节点依次向前两个节点和依次向后两个节点的二阶导数。

对上述差分格式进行泰勒展开并求解差分系数, 可以得到一、二阶导数的普通八阶的紧致差分(compact difference 8, CD8)格式。

一阶导数普通八阶紧致差分格式:

$ \begin{array}{l} O\left( {{h^8}} \right) = \frac{3}{8}{{f'}_{i - 1}} + {{f'}_i} + \frac{3}{8}{{f'}_{i + 1}} - \frac{1}{{480h}}\left( {{f_{i - 3}} - } \right.\\ \left. {{f_{i + 3}}} \right) + \frac{1}{{20h}}\left( {{f_{i - 2}} - {f_{i + 2}}} \right) + \frac{{75}}{{96h}}\left( {{f_{i - 1}} - {f_{i + 1}}} \right) \end{array} $

(2)

二阶导数普通八阶紧致差分格式:

$ \begin{array}{l} O\left( {{h^8}} \right) = \frac{9}{{38}}{{f''}_{i - 1}} + {{f''}_i} + \frac{9}{{38}}{{f''}_{i + 1}} + \frac{1}{{2\;676{h^2}}}\left( {{f_{i - 3}} + } \right.\\ \left. {{f_{i + 3}}} \right) - \frac{{51}}{{1\;520{h^2}}}\left( {{f_{i - 2}} + {f_{i + 2}}} \right) - \frac{{147}}{{152{h^2}}}\left( {{f_{i - 1}} + } \right.\\ \left. {{f_{i + 1}}} \right) + \frac{{1\;017\;169}}{{508\;440{h^2}}}{f_i} \end{array} $

(3)

从公式(2)和公式(3)可知, 普通紧致差分格式若要达到八阶精度, 需用到7个节点。此外, 紧致差分格式也适用于交错网格, 其中八阶紧致交错网格差分(compact staggered difference 8, CSD8)格式可表示为:

$ \begin{array}{l} {a_1}\left( {{{f'}_{i - 1}} + {{f'}_{i + 1}}} \right) + {{f'}_i} = \frac{1}{{\Delta x}}\sum\limits_{n = 1}^3 {{b_n}} \cdot \\ \left\{ {f\left[ {{x_i} + \frac{{\left( {2n - 1} \right)\Delta x}}{2}} \right] - f\left[ {{x_i} - \frac{{\left( {2n - 1} \right)\Delta x}}{2}} \right]} \right\}\\ {a_1} = \frac{{25}}{{118}},{b_1} = \frac{{2\;675}}{{2832}},{b_2} = \frac{{925}}{{5664}},{b_3} = - \frac{{16}}{{28\;320}} \end{array} $

(4)

式中: $ {a_1}, {b_1}, {b_2}, {b_3}$是CSD8格式的差分系数; x_i, Δx分别为x方向上的第i个节点、x方向上的空间步长。为了达到2n+2阶精度, 只需2n个节点函数值。

超紧致有限差分格式在紧致差分格式基础上发展而来, 1998年KRISHNAN提出了五点CCD8格式^[18]:

$ \left\{ \begin{array}{l} {{f'}_i} + \frac{{17}}{{36}}\left( {{{f'}_{i + 1}} + {{f'}_{i - 1}}} \right) - \frac{1}{{12}}h\left( {{{f''}_{i + 1}} + {{f''}_{i - 1}}} \right)\\ \;\;\; = \frac{{107}}{{108h}}\left( {{f_{i + 1}} - {f_{i - 1}}} \right) - \frac{1}{{108h}}\left( {{f_{i + 2}} - {f_{i - 2}}} \right)\\ h{{f''}_i} - \frac{h}{6}\left( {{{f''}_{i + 1}} + {{f''}_{i - 1}}} \right) + \frac{{23}}{{18}}\left( {{{f'}_{i + 1}} - {{f'}_{i - 1}}} \right)\\ \;\;\;\;\frac{{88}}{{27h}}\left( {{f_{i + 1}} + {f_{i - 1}}} \right) - \frac{{13}}{{2h}}{f_i} - \frac{1}{{108h}}\left( {{f_{i + 2}} - {f_{i - 2}}} \right) \end{array} \right. $

(5)

将迎风机制引入五点CCD8格式, 可得到如下的五点七阶迎风超紧致差分(UCCD7)格式, 其精度可达七阶^[19]:

$ \left\{ \begin{array}{l} {{f'}_i} + \frac{{17}}{{18}}{{f'}_{i - 1}} + h\left( { - \frac{1}{{46}}{{f''}_{i + 1}} - \frac{{17}}{{46}}{{f''}_i} + \frac{{10}}{{69}}{{f''}_{i - 1}}} \right) = \\ \;\;\;\;\frac{1}{h}\left( { - \frac{{59}}{{276}}{f_{i + 1}} - \frac{{832}}{{379}}{f_{i - 1}} - \frac{3}{{514}}{f_{i + 2}} - \frac{7}{{552}}{f_{i - 2}} + } \right.\\ \;\;\;\;\left. {\frac{{221}}{{92}}{f_i}} \right)\\ {{f''}_i} - \frac{1}{3}{{f''}_{i - 1}} + \frac{1}{h}\left( { - \frac{{20}}{9}{{f'}_{i + 1}} - 2{{f'}_i} + \frac{1}{3}{{f'}_{i - 1}}} \right) = \\ \frac{1}{{{h^2}}}\left( {\frac{{23}}{{18}}{f_{i + 1}} + \frac{{283}}{{54}}{f_{i - 1}} + \frac{1}{{108}}{f_{i + 2}} - \frac{1}{{36}}{f_{i - 2}} - \frac{{13}}{2}{f_i}} \right) \end{array} \right. $

(6)

根据超紧致差分的原理, 将超紧致差分方法用于求解声波方程初值问题^[20]。在完全弹性介质中声波方程可表示为:

$ \frac{{{\partial ^2}u}}{{\partial {t^2}}} = {v^2}\left( {\frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}u}}{{\partial {z^2}}}} \right) $

(7)

式中: $u\left( {x, z} \right) $为二维介质的质点在x方向、z方向的位移; ${\partial ^2}u/\partial {t^2} $ 为$u\left( {x, z} \right) $对时间t的二阶导数; ${\partial ^2}u/\partial {x^2},{\partial ^2}u/\partial {z^2}$分别为$u\left( {x, z} \right) $在x方向和z方向的二阶导数; v为声波速度。

对(7)式采用四阶中心差分对时间导数离散可得:

$ \begin{array}{*{20}{c}} {{u^{n + 1}} = 2{u^n} - {u^{n - 1}} + {{\left( {\Delta tv} \right)}^2}\left[ {{{\left( {\frac{{{\partial ^2}u}}{{\partial {x^2}}}} \right)}^n} + {{\left( {\frac{{{\partial ^2}u}}{{\partial {z^2}}}} \right)}^n}} \right] + }\\ {\frac{{{{\left( {\Delta tv} \right)}^4}}}{{12}}\left[ {{{\left( {\frac{{{\partial ^4}u}}{{\partial {x^4}}}} \right)}^n} + {{\left( {\frac{{{\partial ^4}u}}{{\partial {z^4}}}} \right)}^n} + 2{{\left( {\frac{{{\partial ^4}u}}{{\partial {x^2}\partial {z^2}}}} \right)}^n}} \right]} \end{array} $

(8)

式中: ${u^{n + 1}}, {u^n}, {u^{n - 1}} $分别为n+1, n和n－1时刻位移; Δt为时间步长。

以z方向为例, 说明利用五点CCD8格式求解公式(8)中u对z方向的空间偏导数的过程, x方向的求解过程不再赘述。假设U为位移场, 其大小为m×n, 则求$\partial u/\partial z $和${\partial ^2}u/\partial {z^2} $的过程可以用矩阵表示为AF=BU。其中A为公式(5)左端的差分系数矩阵, 大小为2m×2m; B为公式(5)右端的差分系数矩阵, 大小为2m×m; F为待求位移场空间一阶和二阶导数矩阵, 大小为2m×n。矩阵分别为:

(9)

$ \mathit{\boldsymbol{F}} = \left[ {\begin{array}{*{20}{c}} {{{\left( {\frac{{\partial u}}{{\partial z}}} \right)}_{1,1}}}&{{{\left( {\frac{{\partial u}}{{\partial z}}} \right)}_{1,2}}}& \cdots &{{{\left( {\frac{{\partial u}}{{\partial z}}} \right)}_{1,n - 1}}}&{{{\left( {\frac{{\partial u}}{{\partial z}}} \right)}_{1,n}}}\\ {{{\left( {\frac{{{\partial ^2}u}}{{\partial {z^2}}}} \right)}_{1,1}}}&{{{\left( {\frac{{{\partial ^2}u}}{{\partial {z^2}}}} \right)}_{1,2}}}& \cdots &{{{\left( {\frac{{{\partial ^2}u}}{{\partial {z^2}}}} \right)}_{1,n - 1}}}&{{{\left( {\frac{{{\partial ^2}u}}{{\partial {z^2}}}} \right)}_{1,n}}}\\ \vdots&\vdots &{}& \vdots&\vdots \\ {{{\left( {\frac{{\partial u}}{{\partial z}}} \right)}_{m,1}}}&{{{\left( {\frac{{\partial u}}{{\partial z}}} \right)}_{m,2}}}& \cdots &{{{\left( {\frac{{\partial u}}{{\partial z}}} \right)}_{m,n - 1}}}&{{{\left( {\frac{{\partial u}}{{\partial z}}} \right)}_{m,n}}}\\ {{{\left( {\frac{{{\partial ^2}u}}{{\partial {z^2}}}} \right)}_{m,1}}}&{{{\left( {\frac{{{\partial ^2}u}}{{\partial {z^2}}}} \right)}_{m,2}}}& \cdots &{{{\left( {\frac{{{\partial ^2}u}}{{\partial {z^2}}}} \right)}_{m,n - 1}}}&{{{\left( {\frac{{{\partial ^2}u}}{{\partial {z^2}}}} \right)}_{m,n}}} \end{array}} \right] $

(10)

(11)

根据AF=BU, 可得F=A^－1BU, 其中F的奇数行矩阵元素为$ \partial u/\partial z$, 偶数行矩阵元素为$ {\partial ^2}u/\partial {z^2}$。同理, 将位移场转置, 可求得$ \partial u/\partial x$和${\partial ^2}u/\partial {x^2} $。利用UCCD7格式求公式(8)中的空间偏导数的方法与此相同, 仅有差分系数矩阵不同, 故不再赘述。

以声波方程为例, 说明了五点CCD8格式或五点UCCD7格式在地震声波波场数值模拟中的应用方法, 该方法同样可以用于弹性波方程的数值模拟。各向同性完全弹性介质中的二维弹性波方程可以表示为:

$ \left\{ \begin{array}{l} \frac{{{\partial ^2}u}}{{\partial {t^2}}} = v_{\rm{P}}^2\left( {\frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}w}}{{\partial x\partial z}}} \right) + v_{\rm{S}}^2\left( {\frac{{{\partial ^2}u}}{{\partial {z^2}}} - \frac{{{\partial ^2}w}}{{\partial x\partial z}}} \right)\\ \frac{{{\partial ^2}w}}{{\partial {t^2}}} = v_{\rm{P}}^2\left( {\frac{{{\partial ^2}w}}{{\partial {z^2}}} + \frac{{{\partial ^2}u}}{{\partial x\partial z}}} \right) + v_{\rm{S}}^2\left( {\frac{{{\partial ^2}w}}{{\partial {x^2}}} - \frac{{{\partial ^2}u}}{{\partial x\partial z}}} \right) \end{array} \right. $

(12)

式中:$u\left( {x, z} \right) $和$w\left( {x, z} \right) $分别表示x方向和z方向的质点位移; $ {v_{\rm{P}}}$和${v_{\rm{S}}} $分别表示纵、横波速度。u和w的时间四阶精度差分格式为^[21]:

$ \begin{array}{l} {u^{n + 1}} = 2{u^n} - {u^{n - 1}} + {\left( {\Delta t} \right)^2}\left\{ {v_{\rm{P}}^2\left[ {{{\left( {\frac{{{\partial ^2}u}}{{\partial {x^2}}}} \right)}^n} + } \right.} \right.\\ \left. {\left. {{{\left( {\frac{{{\partial ^2}w}}{{\partial x\partial z}}} \right)}^n}} \right]} \right\} + {\left( {\Delta t} \right)^2}\left\{ {v_{\rm{S}}^2\left[ {{{\left( {\frac{{{\partial ^2}u}}{{\partial {x^2}}}} \right)}^n} - {{\left( {\frac{{{\partial ^2}w}}{{\partial x\partial z}}} \right)}^n}} \right]} \right\} + \\ \frac{{{{\left( {\Delta tv} \right)}^4}}}{{12}}\left\{ {v_{\rm{P}}^4\left[ {{{\left( {\frac{{{\partial ^4}u}}{{\partial {x^4}}}} \right)}^n} + {{\left( {\frac{{{\partial ^4}w}}{{\partial {x^3}\partial z}}} \right)}^n} + {{\left( {\frac{{{\partial ^4}u}}{{\partial {x^2}\partial {z^2}}}} \right)}^n} + } \right.} \right.\\ \left. {\left. {{{\left( {\frac{{{\partial ^4}w}}{{\partial x\partial {z^3}}}} \right)}^n}} \right]} \right\} + \frac{{{{\left( {\Delta tv} \right)}^4}}}{{12}}\left\{ {v_{\rm{S}}^4\left[ {{{\left( {\frac{{{\partial ^4}u}}{{\partial {z^4}}}} \right)}^n} - {{\left( {\frac{{{\partial ^4}w}}{{\partial {x^3}\partial z}}} \right)}^n} + } \right.} \right.\\ \left. {\left. {{{\left( {\frac{{{\partial ^4}u}}{{\partial {x^2}\partial {z^2}}}} \right)}^n} - {{\left( {\frac{{{\partial ^4}w}}{{\partial x\partial {z^3}}}} \right)}^n}} \right]} \right\} \end{array} $

(13)

$ \begin{array}{l} {w^{n + 1}} = 2{w^n} - {w^{n - 1}} + {\left( {\Delta t} \right)^2}\left\{ {v_{\rm{P}}^2\left[ {{{\left( {\frac{{{\partial ^2}w}}{{\partial {z^2}}}} \right)}^n} + } \right.} \right.\\ \left. {\left. {{{\left( {\frac{{{\partial ^2}u}}{{\partial x\partial z}}} \right)}^n}} \right]} \right\} + {\left( {\Delta t} \right)^2}\left\{ {v_{\rm{S}}^2\left[ {{{\left( {\frac{{{\partial ^2}w}}{{\partial {x^2}}}} \right)}^n} - {{\left( {\frac{{{\partial ^2}u}}{{\partial x\partial z}}} \right)}^n}} \right]} \right\} + \\ \frac{{{{\left( {\Delta tv} \right)}^4}}}{{12}}\left\{ {v_{\rm{P}}^4\left[ {{{\left( {\frac{{{\partial ^4}w}}{{\partial {z^4}}}} \right)}^n} + {{\left( {\frac{{{\partial ^4}u}}{{\partial {x^3}\partial z}}} \right)}^n} + {{\left( {\frac{{{\partial ^4}w}}{{\partial {x^2}\partial {z^2}}}} \right)}^n} + } \right.} \right.\\ \left. {\left. {{{\left( {\frac{{{\partial ^4}u}}{{\partial x\partial {z^3}}}} \right)}^n}} \right]} \right\} + \frac{{{{\left( {\Delta tv} \right)}^4}}}{{12}}\left\{ {v_{\rm{S}}^4\left[ {{{\left( {\frac{{{\partial ^4}w}}{{\partial {x^4}}}} \right)}^n} - {{\left( {\frac{{{\partial ^4}u}}{{\partial {x^3}\partial z}}} \right)}^n} + } \right.} \right.\\ \left. {\left. {{{\left( {\frac{{{\partial ^4}w}}{{\partial {x^2}\partial {z^2}}}} \right)}^n} - {{\left( {\frac{{{\partial ^4}u}}{{\partial x\partial {z^3}}}} \right)}^n}} \right]} \right\} \end{array} $

(14)

由公式(13)和公式(14)可以看出, 弹性波方程与声波方程的差分格式相比, 增加了许多二阶和四阶混合偏导数。对于这些混合偏导数, 均可以按不同方向进行多次导数求取, 求取顺序对结果无影响, 求取方法不再赘述。

2 数值频散、稳定性及精度分析

对波动方程的时间和空间偏导数的离散化造成了数值频散, 它使相速度变成了网格间距的函数。数值模拟时, 如果空间网格长度过大, 会造成较大的求解误差, 产生数值频散, 降低模拟精度^[22-23]。压制数值频散是采用有限差分方法时面临的关键问题之一^[24]。频散关系分析既是判断数值模拟所用方法优劣的重要依据, 也是确定空间网格大小的重要依据^[25], 五点CCD8的二阶导数数值波数可以由公式(15)表示:

$ \begin{array}{l} {\left( {{{k''}_x}h} \right)^2} = - \left[ { - 13/2 + 2\left( {17/36} \right)\left( {88/27} \right) + 2461/} \right.\\ \left. {972 + 2\left( {88/27 - 221/72 - 23/972} \right)\cos {k_x}h} \right]/\left[ {1 - } \right.\\ 17/108 - 23/108 + 2\left( { - 1/6 + 17/36} \right)\cos {k_x}h + \\ \left. {2\left( { - 17/216 + 23/216} \right)\cos 2{k_x}h} \right] - \left[ {2\left( { - 1/108} \right) + } \right.\\ \left. {374/243 - 1391/1944} \right)\cos 2{k_x}h + 4\left( { - 17/3888 + } \right.\\ \left. {\left. {23/1944} \right)\cos {k_x}h\cos 2{k_x}h} \right]/\left[ {1 - 17/108 - 23/108 + } \right.\\ 2\left( { - 1/6 + 17/36} \right)\cos {k_x}h + 2\left( { - 17/216 + } \right.\\ \left. {\left. {23/216} \right)\cos 2{k_x}h} \right] \end{array} $

(15)

式中:h为空间网格大小; k为真波数; k_x为x方向的真波数; θ为波的传播方向与x轴的夹角; ${k_x}h = kh{\rm{cos}}\theta ;{{k''}_x} $为x方向的数值波数。

同理可得z方向的数值波数${{k''}_z} $, 将${{k''}_x} $和k″₂代入声波方程的差分格式(8), 利用欧拉公式, 化简可得频散关系:

$ \begin{array}{l} 2\cos \left( {kv'\Delta t} \right) = 2 - {\left( {v\Delta \frac{t}{h}} \right)^2}\left[ {{{\left( {{{k''}_x}h} \right)}^2} + {{\left( {{{k''}_z}h} \right)}^2}} \right] + \\ {\left( {v\Delta \frac{t}{h}} \right)^4} + \left[ {{{\left( {k''h} \right)}^4} + 2{{\left( {{{k''}_x}h \cdot {{k''}_z}h} \right)}^4}} \right]/12 \end{array} $

(16)

式中:v′为数值波速; v是真波速; vΔt/h称为courant数。

数值波速与真波速的比值γ可表示为:

$ \gamma = \frac{{v'}}{v} = \frac{{kv'\Delta t}}{{kv\Delta t}} = \frac{{kv'\Delta t}}{{k\frac{{v\Delta t}}{h}h}} = \frac{{kv'\Delta t}}{{\frac{{v\Delta t}}{h}kh}} $

(17)

将简谐波解$u_{i, j}^n = A{\rm{exp}}[{\rm{I}}(ih{k_x} + jh{k_z} - kv\prime n\Delta t)] $代入公式(5)得到如下公式:

$ \left\{ \begin{array}{l} {{\varphi '}_x} + \frac{{17}}{{18}}{{\varphi '}_x} \cdot \exp \left( { - {\rm{I}}{\varphi _x}} \right) - \left[ { - \frac{1}{{46}}{{\left( {{{\varphi ''}_x}} \right)}^2} \cdot } \right.\\ \;\;\left. {\exp \left( {{\rm{I}}{\varphi _x}} \right) - \frac{{17}}{{46}}{{\left( {{{\varphi ''}_x}} \right)}^2} + \frac{{10}}{{69}}{{\left( {{{\varphi ''}_x}} \right)}^2} \cdot \exp \left( { - {\rm{I}}{\varphi _x}} \right)} \right]\\ \;\; = - \frac{{59}}{{276}}\exp \left( {{\rm{I}}{\varphi _x}} \right) - \frac{{832}}{{379}}\exp \left( { - {\rm{I}}{\varphi _x}} \right) - \\ \frac{3}{{514}}\exp \left( {{\rm{2I}}{\varphi _x}} \right) - \frac{7}{{552}}\exp \left( { - {\rm{2I}}{\varphi _x}} \right) + \frac{{221}}{{92}}\\ \;\;\;\;\;\;\frac{1}{3}{\left( {{{\varphi ''}_x}} \right)^2}\exp \left( {{\rm{I}}{\varphi _x}} \right) - {\left( {{{\varphi ''}_x}} \right)^2} + \left[ \begin{array}{l} - \frac{{20}}{9}{{\varphi '}_x} \cdot \\ \exp \left( {{\rm{I}}{\varphi _x}} \right) - \end{array} \right.\\ \;\;\;\;\;\;\left. {2{{\varphi '}_x} + \frac{1}{3}{{\varphi '}_x} \cdot \exp \left( {{\rm{I}}{\varphi _x}} \right)} \right] = \frac{{23}}{{18}}\exp \left( {{\rm{I}}{\varphi _x}} \right) + \\ \;\;\;\;\;\;\;\frac{{283}}{{54}}\exp \left( { - {\rm{I}}{\varphi _x}} \right) + \frac{1}{{108}}\exp \left( {2{\rm{I}}{\varphi _x}} \right)\\ \;\;\;\;\;\;\;\frac{1}{{36}}\exp \left( { - 2{\rm{I}}{\varphi _x}} \right) - \frac{{13}}{2} \end{array} \right. $

(18)

为方便表达, 引入中间变量$ {\varphi _x}, {{\varphi '}_x}, {{\varphi ''}_x}。{\varphi _x} = h{k_x}, {{\varphi '}_x} = h{{k'}_x}, {{\varphi ''}_x} = h{{k''}_x}, {\rm{I}} = \sqrt { - 1} $。数值求解公式(18)得到UCCD7的二阶导数数值波数, 按照CCD8数值波速与真波速的比值求取方法可求得UCCD7的数值波速与真波度的比值。

将多种有限差分方法的速度比进行比较, 结果如图 1所示。图 1a, 图 1b和图 1c展示了不同的θ条件下4种常见网格差分格式的速度比曲线, 其中α=vΔt/h; 图 1d中增加了八阶紧致交错网格格式的速度比曲线; 图 1e增加了五点七阶迎风超紧致差分格式的速度比曲线。考虑到参与对比的方法的稳定性, 此处courant数均为0.25, 横坐标φ∈[0, π], 为波数与空间步长的乘积, 纵坐标γ为速度比。速度比为1表示数值波速与理论波速一致, 说明该方法能很好地压制数值频散, 反之, 则说明该方法的数值频散严重。

由图 1可以看出, 五点CCD8格式的速度比曲线偏离1时, 其横坐标的数值最大, 所以五点CCD8格式压制数值频散的能力强。图 1e说明将五点CCD8格式引入迎风机制后的五点UCCD7格式, 压制频散效果并无明显降低。五点UCCD7, 五点CCD8, 三点CCD6, CD6, CD8和CSD8格式在偏离1时所对应的横坐标数值分别为1.57, 1.69, 1.12, 0.98, 1.55和1.61, 这6种方法在保证无数值频散时所需的最低网格点数分别为4.00, 3.70, 5.60, 6.40, 4.05和3.90个。五点CCD8格式的数值频散的压制效果最好, CD6格式的数值频散的压制效果最差。

稳定性条件是影响差分方法计算效率的重要因素。本文使用Fourier方法^[26-27]进行稳定性分析, 略去推导过程, 直接给出声波方程时间四阶精度的五点CCD8和五点UCCD7格式的稳定性条件α_CCD8和α_UCCD7分别为:

$ {\alpha _{{\rm{CCD8}}}} = {v_{\max }}\Delta \frac{t}{h} < 0.74,{\alpha _{{\rm{UCCD7}}}} = {v_{\max }}\Delta \frac{t}{h} < 0.80 $

(19)

式中:v_max是指正演模拟时, 所建立模型的最大声波速度, 与五点CCD8格式相比, 引入迎风机制的五点UCCD7格式的稳定性有了一定的提高, 可使用略大的时间网格, 在一定程度上减少计算时间, 提高计算效率。

截断误差决定了有限差分格式和地震波场数值模拟精度, 将文中所提到的几种方法的二阶导数差分格式截断误差进行对比。由Taylor级数展开^[28]可以计算得到CD6、三点CCD6、CD8和五点CCD8格式的截断误差分别为－4.2163×10^-4, －4.9603×10^-5, 2.7000×10^－5和2.2046×10^-6, 即五点CCD8格式的截断误差最小, 数值模拟精度最大。

3 PML边界条件

数值模拟时, 由于模型大小的限制, 必然会存在人工边界, 而人工边界的处理直接影响到数值模拟的精度及计算效率, 若不对其进行处理则会干扰正常的地震波场。因此, 本文在模型试算时, 采用完全匹配层(perfectly matched layer, 简称PML)对人工边界进行处理, 其基本思想是:在吸收边界区域匹配与计算区域相同的波阻抗, 使入射波无反射的进入吸收层, 吸收层为衰减介质, 可使得入射波迅速衰减直至消除^[29]。理论上, PML边界能够吸收任何入射角度和频率的入射波, 实践也证明PML吸收边界条件应用效果良好, 可应用于声波和弹性波的研究^[30-31], 下式为声波方程的PML边界条件的控制方程:

$ \left\{ \begin{array}{l} u = {u_1} + {u_2}\\ \frac{{\partial {u_1}}}{{\partial t}} + d\left( x \right){u_1} = v_p^2\frac{{\partial J}}{{\partial x}}\\ \frac{{\partial {u_2}}}{{\partial t}} + d\left( z \right){u_2} = v_p^2\frac{{\partial L}}{{\partial z}}\\ \frac{{\partial J}}{{\partial t}} + d\left( x \right)J = \frac{{\partial {u_1}}}{{\partial x}} + \frac{{\partial {u_2}}}{{\partial x}}\\ \frac{{\partial L}}{{\partial t}} + d\left( z \right)L = \frac{{\partial {u_1}}}{{\partial z}} + \frac{{\partial {u_2}}}{{\partial z}} \end{array} \right. $

(20)

式中:u₁, u₂, J和L均为中间变量; d(x), d(z)分别为x, z方向的衰减系数, 用于衰减这两个方向的波场。由文献[32]可知, 由于d(x)和d(z)的衰减函数和层数影响了边界反射的衰减效果, 因此余弦型吸收衰减函数效果较好, 所以本文在模型试算中采用如下的余弦型吸收衰减函数:

$ \begin{array}{l} {\alpha _i} = \lambda \left\{ {1 - \cos \left[ {\frac{{{\rm{ \mathsf{ π} }}\left( {{P_{{\rm{ML}}}} - i} \right)}}{{2{P_{{\rm{ML}}}}}}} \right]} \right\}\\ i = 0,1,2, \cdots ,{P_{{\rm{ML}}}} \end{array} $

(21)

式中:α_i代表吸收衰减因子; λ为衰减幅度因子, λ=500;P_ML为吸收边界的层数, P_ML=20。

4 模型试算 4.1 声波方程均匀介质模型波场模拟

均匀介质模型的长度及深度均为4800m, 声波速度为4000m/s, 激发位置位于模型中心, 震源为20Hz雷克子波, 设置纵、横向空间步长均为48m, 时间步长为1ms。

图 2为400ms时刻不同模拟方法得到的波场快照。可以看出, 采用CD6和CD8格式模拟得到的波场快照存在严重的频散; 采用三点CCD6格式得到的结果仍有轻微的频散; 采用CSD8、五点CCD8和五点UCCD7格式模拟得到的波场清晰, 无数值频散。

4.2 弹性波方程均匀介质模型波场模拟

采用五点CCD8格式及五点UCCD7格式进行弹性波方程数值模拟。模型长度和深度均为4000m, 纵横向空间步长均为10m, 介质为均匀介质, 纵波速度4000m/s, 横波速度2310m/s, 模型中心处激发, 震源为20Hz雷克子波。参考五点CCD8格式与五点UCCD7格式的最大courant数, 设置时间步长为1ms。图 3为分别采用两种格式在400ms时刻的波场快照, 其中图 3a和图 3b为采用五点CCD8格式所得到的结果, 图 3c和图 3d为采用五点UCCD7格式得到的结果。波场均清晰, 无数值频散, 纵波和横波位移在水平和垂向记录上的特征明显, 说明采用五点CCD8格式与五点UCCD7格式均可有效而清晰地模拟弹性波在各向同性介质模型中的传播。

4.3 声波方程水平层状均匀介质模型数值模拟

水平层状均匀介质模型包括两层, 各层的厚度均为1000m, 长度与深度均为2000m。第1层和第2层的地震波速分别为4000m/s和4500m/s。在地表x=1000m处激发, 震源为20Hz雷克子波, 时间步长Δt=0.5ms, 空间步长Δx=Δz=20m, 地表接收的单炮记录长度为1000ms。图 4为采用五点CCD8及UCCD7格式对水平层状介质模型进行声波方程数值模拟得到的结果, 地震记录清晰, 无数值频散和不稳定数值结果, 地震记录中的直达波和反射波显示清楚, 说明采用五点CCD8格式与五点UCCD7格式可清晰有效地模拟声波在多层各向同性介质模型中的传播。

4.4 声波方程Marmousi模型数值模拟

基于二维Marmousi纵波速度模型进行数值模拟。在地表中心处(1250m, 0)激发震源为30Hz雷克子波, 时间步长与空间步长分别为Δt=0.2ms, Δx=Δz=10m, 采样时间为2000ms。图 5为对Marmousi模型采用五点CCD8格式及五点UCCD7格式进行声波方程模拟得到的地面地震记录, 地震记录平滑清晰, 边界吸收效果好, 无明显边界反射和数值频散, 验证了五点CCD8格式和五点UCCD7格式对复杂模型的适用性。

5 结论与认识

本文重点研究了五点八阶超紧致有限差分(CCD8)和五点七阶迎风超紧致有限差分(UCCD7)两种格式, 并将其应用于声波和弹性波方程的数值模拟, 得到3点认识:

1) 与传统的有限差分格式相比, 五点CCD8格式仅需5个节点就能使二阶空间偏导数的差分精度达到八阶, 该格式具有低数值频散、高模拟精度及高稳定性的特点;

2) 引入迎风机制, 优化五点CCD8格式, 得到了五点UCCD7格式, 提高了声波方程差分格式的稳定性, 同时保持了较好的低数值频散的优势;

3) 模型试算结果说明, 五点CCD8格式与五点UCCD7格式不仅适用于声波方程的波场模拟, 还可以进一步推广应用于二维或三维的各向异性介质和粘弹介质等复杂介质的弹性波方程数值模拟。

本文迭代求取了空间四阶导数, 这导致了四阶导数的模拟精度降低。在今后的研究中, 可建立一种新的组合型紧致差分格式, 一次性求取四阶导数及混合偏导数, 以提高数值模拟的精度。