紧致交错网格优化差分系数二维声波方程数值模拟

2019-10-08蔡文杰桂志先

石油地球物理勘探 2019年5期

汪勇王鹏蔡文杰桂志先

(油气资源与勘探技术教育部重点实验室(长江大学)，湖北武汉 430100)

0 引言

地震波场数值模拟是应用数值方法模拟地震波在地下介质中的传播过程，得到在地面或地下各观测点地震记录的一种方法，是地球物理勘探和地震学的重要基础。随着有限差分数值模拟技术的发展，针对实际应用中提高有限差分计算效率[1]、模拟精度[2]、算法稳定性[3-4]，处理复杂介质[5-7]、吸收边界条件[8-10]和压制数值频散[11-13]等多方面需求，已研发了多种实用技术，并取得相应研究成果。

提高差分格式数值计算精度最直接的方法就是在模拟时增加(加密)网格节点数，但随之也增加了运算量和存储空间。紧致差分方法能较好地解决这个矛盾；同时紧致差分是一种隐式差分格式，稳定性较强。这些优势使其成为目前颇受关注的有限差分方法之一，并被广泛应用于声波、弹性波和复杂介质等地震波场数值模拟中[14-22]，提高了地震波场数值模拟的精度和计算效率。

而Madariaga[23]提出的交错网格差分格式既可提高数值模拟的局部精度，还能加快收敛速度。将交错网格与紧致差分技术结合，可进一步提高数值模拟的精度和压制数值频散的能力[24-26]。

为使差分格式在较大波数范围内减少频散误差，Kim等[27]基于频散关系保持(Dispersion relation preserving，DRP)的思路优化了紧致差分格式，指明采用这种优化差分系数方法可使差分算子最大程度地逼近空间偏导数算子。Liu等[28]基于频散关系提出了优化的时空域有限差分系数，在不增加计算成本的情况下可显著提高模拟精度。Zhang等[29-30]应用最大范数的目标函数及模拟退火算法求解目标函数，优化了一阶和二阶常规中心差分格式的差分系数。Liu[31-32]利用最小平方法优化了二阶导数中心差分和一阶导数交错差分系数。Yu等[33]基于DRP思想，优化得到5阶精度的组合型紧致差分系数。Ren等[34]以优化后时空域差分格式进行声波和弹性波方程的数值模拟。Yang等[35-36]利用极值逼近(Minimax approximation)算法优化了交错网格差分系数，当波数较大时可提高模拟精度。

本文首先讨论了一阶导数的2N阶紧致交错差分格式的建立方法；基于频散关系保持的思想，通过最小平方法建立了最小数值波数误差的目标函数，利用拉格朗日乘数法求解目标函数，得到优化后的4～10阶精度差分系数；然后分析和对比优化前、后的紧致交错网格差分格式的模拟精度、数值频散和声波方程的稳定性条件，并用优化后的紧致交错差分格式对一阶速度—应力声波方程组进行数值模拟。研究结果表明，采用优化差分算子能有效压制数值频散，提高差分近似导数的精度。

1 二维声波方程的高阶差分近似

Madariaga[23]早年提出的波动方程交错网格数值模拟方法的差分精度为O(Δt2+Δx2)； Levander[37]应用交错网格方法求解P-SV波方程，将差分精度提高到O(Δt2+Δx4)；董良国等[3]提出了一阶弹性波方程交错网格高阶差分解法，其差分精度达到O(Δt2M+Δx2N)。现今，交错网格有限差分方法已成为一种高效且应用广泛的数值模拟方法。

根据弹性力学分析，二维情况下一阶应力—速度声波方程组(假设体力为零)可表示为

(1)

式中：P为应力(标量)；Vx和Vz分别表示x和z方向质点振动速度分量；ρ为介质密度；v为地震波速度。

1.1 时间2M阶近似

用交错网格法进行声波方程数值模拟时，速度和应力分量分别是在t和t+Δt/2时刻计算的，其中Δt为时间步长。用泰勒公式可得速度分量Vx的2M阶精度时间递推式

(2)

当M=1时，即可得到时间二阶精度差分近似

(3)

同理可得速度分量Vz和应力分量P的时间二阶精度差分近似

(4)

(5)

利用应力—速度方程组，将式(3)～式(5)中的变量对时间的导数转换为对空间的导数，则可得

(6)

1.2 空间2N阶紧致交错差分近似

紧致交错差分格式最早由Nagarajan等[38]提出，用于大涡模拟(Large eddy simulation)问题，该差分格式利用4个网格节点可得到6阶空间精度。随后，Boersma[39]提出最高可达12阶空间精度的紧致交错差分格式，用于可压缩流体的Navier-Stokes方程的数值模拟。类似于常规交错网格数值模拟，紧致交错网格法求取变量的导数时，也是在相应的变量网格点之间的半程上进行的。

一阶导数的2N阶精度紧致交错差分格式为

(7)

式中：h为空间网格尺寸；a和bn是差分系数。利用式(7)近似计算式(6)中变量对空间的一阶导数，确保2N阶差分精度的关键是确定差分系数a和bn。利用泰勒级数展开和待定系数法，即用

(8)

可求得紧致交错网格2N阶空间差分精度的差分系数。如当差分精度2N=10时，上述方程组共有5个方程，可确定a和bn(n=1～4)共5个差分系数。

此外，从上述公式可看出，常规交错差分格式(Staggered finite difference，SD)求取中心点处的导数值时，仅用到周围网格点的函数值，而紧致交错差分格式(Staggered compact finite difference，SCD)却还用到相邻点导数值。紧致交错网格法只需利用2N-2个节点即可达到2N阶空间差分精度，而常规交错网格法则须使用2N个节点，因此SCD方法可节省存储空间。表1列出4～10阶泰勒级数展开得到的差分系数。

表1 紧致交错差分格式一阶导数的差分系数

1.3 差分格式的矩阵形式

(9)

其中

向前差分为

向后差分为

据式(9)，由应力(标量)场P求取其空间一阶导数

(10)

上述各式中差分系数a和bn(n=1～4)由表1给出，可求得一阶导数的4～10阶空间差分精度近似值。

2 紧致交错格式差分系数优化及分析

在地震波场数值模拟中，为得到清晰的波场记录，需要高精度的空间离散差分算子。差分算子的精度取决于差分系数和差分阶数，差分阶数越高，差分算子越精确，但其运算量也会随之增加。优化差分系数可使差分算子最大程度地逼近空间偏导数算子，本节对一阶导数的紧致交错差分系数进行优化，进一步提高紧致交错差分格式的分辨率，并分析优化后的差分格式的精度、频散和稳定性。

2.1 优化差分系数

首先对一阶导数的紧致交错差分格式(式(7))进行数值频散分析。令

(11)

式中：h为网格尺寸；k为真波数。若数值模拟时不存在误差，则B=(ik)A。但在差分计算中，可能会产生不同结果，即产生数值波数(又称修正波数)。可令B=(ik′)A,k′为修正波数。

将式(11)代入式(7)，并用欧拉公式化简可得

(12)

式中：ω=kh;ω′=k′h。

在理想情况下，若不存在数值频散，则ω′=ω。它们的差异越大，说明该方法的数值频散越严重；反之则说明该方法能更好地压制数值频散。依据修正波数尽可能在大范围内接近真波数，以4点紧致交错差分格式为例说明优化的思路和方法。

由式(7)可得

(13)

为了满足2阶和4阶精度泰勒公式截断误差的要求，式中差分系数a、b1和b2须满足下列方程组

(14)

即为式(8)表示的前两个方程。

为了确定上述方程中的未知差分系数，按照Tam等[40]的思路和方法，在某个选定波数范围内，确定式(13)中3个未知差分系数，使得修正波数尽可能地接近真波数，故定义误差函数为

(15)

(16)

其中

(17)

由式(14)和式(16)中的3个方程，可确定3个优化后差分系数，即a=0.185715、b1=0.985714、b2=0.128572。将优化后的差分系数代入式(13)，则得到优化后的紧致交错差分格式(Optimized staggered compact finite difference，OSCD)。

需要说明的是： 4阶精度OSCD格式在式(14)基础上，需加入一个表示DRP要求的积分方程，所以不能严格满足空间6阶精度泰勒公式截断误差的要求，只能达到空间4阶精度近似；而常规SCD格式的3个方程即能满足空间6阶精度截断误差的要求。因此，若要达到2N阶空间差分精度，SCD方法只需利用2N-2个节点，而OSCD方法需使用2N个节点，与SD方法使用的节点数相同。

采用同样方法，对式(7)表示的6～10阶OSCD格式做差分系数优化，优化后差分系数列于表2。

表2 优化的紧致交错网格一阶导数的差分系数

本文方法与文献[33]方法基本一样，不同之处在于：文献中原差分格式是3点6阶组合型紧致差分格式，为了加入代表DRP思路的积分方程，在原差分格式基础上构造了新格式，增加了1个网格节点，降低了1阶精度，得到的是4点5阶迎风型组合紧致差分格式；本文的优化格式与原紧致交错差分格式相比，在使用同样网格节点进行差分计算时，降低了2阶精度。

本文求取优化差分系数的方法本质上是根据泰勒级数展开式，而文献[35]和[36]则是根据函数逼近的思路。文献[35]首先根据紧致交错差分格式的频散关系和极小极大原则建立了误差函数，在所有逼近格式中寻找误差最小且满足误差标准的三角函数多项式，求取过程中使用的是列梅兹(Remez)迭代算法。文献[35]指出该方法求得的差分系数在中等和大波数条件下比泰勒级数得到的要好，但算法复杂，实现难度较大，这也限制了该方法的应用。此外，该方法得到的差分格式是在给定的误差标准条件下得到的，但不能明确差分格式的近似精度，即通常意义上的2N阶。

2.2 频散分析

对优化前、后一阶导数的紧致交错差分格式进行数值频散分析，通过数值波数与真波数的比值评判优化效果并确定适用的空间网格尺寸。

据式(12)，修正波数与真波数之比可定义为

(18)

在理想情况下，若不存在数值频散则波数比R恒等于1。R偏离1越大，说明该方法的数值频散越严重；反之则说明该方法能较好地压制数值频散。利用表1和表2中不同精度的差分系数可求取优化前后的波数比曲线。为了比较，计算得到常规交错差分格式(SD)的波数比，计算公式见文献[35]，不再赘述公式。

图1为3种差分格式在不同差分网格节点个数(即差分算子长度)情况下的波数比曲线，如图1a为4个差分节点，根据前文分析，SD、SCD和OSCD格式分别可达到4阶、6阶和4阶差分精度。

此外，从图中可见：低阶差分精度时，三者波数比曲线差别最大；随差分精度增加，差异变小。显然，OSCD格式在压制数值频散方面的优势主要体现在低阶差分精度时。因此，可用较低(4或6)阶OSCD格式达到高阶SCD或SD格式的效果。

图1 不同网格节点数三种差分格式的波数比曲线

2.3 精度分析

为了比较SD、SCD和OSCD三种格式的近似精度，对比分析其截断误差。从表3可见：在达到相同空间近似精度的条件下，4和6阶精度的OSCD格式具有更小的截断误差，8和10阶时误差介于SD和SCD之间，所以低阶差分时OSCD格式的精度优势更明显。

表3 不同差分格式的一阶导数截断误差主项系数比较

利用一维平面简谐波初值问题，比较SD、SCD和OSCD格式的数值模拟精度。一维平面谐波初值问题可表示为

(19)

该偏微分方程的达朗贝尔解，即精确的解析解为

(20)

图2 三种差分格式的一维声波方程数值模拟快照

通过定量计算，可得到该时刻、4000～5000m之间，OSCD、SCD和SD格式数值解与解析解之间的相对误差依次为0.33%、0.77%和0.97%，这正是由于它们在计算空间导数时存在不同截断误差造成的，这也与表3的理论分析结果相一致。相对误差定义为

(21)

2.4 稳定性分析

按照董良国等[3]和杜启振等[4]使用的Fourier方法对差分格式进行稳定性分析，略去推导过程，直接给出式(22)表示的二维声波一阶速度—应力方程组(2M，2N)阶差分精度的紧致交错网格差分格式的稳定性条件。

(22)

其中

表4 二维声波方程不同差分格式的Courant数

3 模型试算

采用OSCD差分格式对简单模型和复杂的Marmousi模型进行二维声波方程进行数值模拟，检验该方法数值模拟的适用性。

3.1 均匀介质模型

为比较SD、SCD和OSCD三种差分方法的差异，设置模型尺寸为3000m×3000m，纵波速度为2000m/s，模型中心加载正弦衰减子波，峰值频率为30Hz，空间网格为12m×12m，时间步长为1ms。采用三种格式的时间2阶、空间4～8阶差分精度对一阶速度—应力声波方程进行波场模拟，得到同一时刻的应力P分量的波场快照(图3)。

在时间和空间网格尺寸相同的情况下，图3a所示的4～8阶SD格式的波场快照频散非常严重；图3b的4～8阶SCD格式波场快照好于SD格式，但4阶SCD格式频散依然明显，6阶在0°和90°方向上也存在较明显频散，8阶精度时波场快照较清晰，频散得到较好压制；图3c的4～8阶OSCD格式波场快照均很清晰，几乎看不到数值频散。与图1的理论分析结果一致，4阶OSCD格式(图3c左)不仅明显好于4阶SCD格式(图3b左)，也要好于使用了相同节点数的6阶SCD格式(图3b中)，达到了8阶SCD格式(图3b右)的效果。若要SD和SCD格式达到与OSCD相同的数值模拟效果，在相同网格步长条件下，必须增加差分节点提高差分精度，所以OSCD格式能使用较少的差分节点，从而减少了运算量，提高了计算效率。

3.2 水平层状介质模型

设置四层水平层状介质模型，模型长度和深度均为2000m，各层厚度均为500m，纵波速度以500m/s间隔，从3000m/s增至4500m/s，密度由Gardner经验公式得到。在地表(1000m，0)处激发主频为20Hz的雷克子波(震源)，时间步长Δt=1ms，记录时间长度为1s。边界处理采用分裂的PML边界条件对边界反射进行吸收衰减，其控制方程见文献[7]，此处不赘述。

理论分析和均匀介质模拟都说明了常规SD格式与SCD和OSCD格式相比，不利于压制数值频散，所以下文中不再对常规SD格式进行分析，只比较SCD和OSCD的数值模拟结果。在都使用4个网格节点情况下，也就是利用4阶OSCD格式和6阶SCD格式进行数值模拟，并采用不同的网格尺寸进行计算，数值模拟地表接收的单炮记录(长度为1s)如图4所示。所有地震记录使用同样的参数进行增益处理，并提取x=1000m处质点模拟得到的不含直达波单道地震记录(图5)。

从图4和图5可看出，随着网格尺寸h的增加，数值频散总体上逐渐变得严重。图4显示6阶SCD格式和4阶OSCD格式分别最大在15m和17m时地震记录清晰，无明显数值频散和边界反射；图5中单道记录在波尾处也无抖动，SCD和OSCD分别从16m和18m开始出现“拖尾”现象。仅以该模型的单个变量而言，15m网格SCD格式的网格数约为133×133，17m网格OSCD格式的网格数为117×117，相对减少内存23%，也减少了运算时间，这对于大尺度模型的数值模拟或逆时偏移还是比较有意义的。

图3 4阶(左)、6阶(中)、8阶(右)三种格式模拟的600ms时刻波场快照

3.3 Marmousi模型

为了验证OSCD格式对复杂介质的适用性，用经典二维Marmousi纵波速度模型(图6)做数值模拟。模型的速度范围是1729～5500m/s，根据Gardner公式ρ=0.31v0.25，得到密度范围是2.00～2.67g/cm3。模型范围是501×501个网格点，空间网格尺寸取12m，时间步长取0.8ms，纵波震源位于地表中心位置，激发30Hz的Ricker子波，采样时间3s，采用20层PML吸收边界。

在上节主要对比了4阶OSCD与6阶SCD格式，本节对比6阶OSCD与8阶SCD格式数值模拟的结果。一方面验证OSCD格式对复杂介质的适用性；另一方面验证本文方法对其他精度SCD格式的优化效果。图7是分别使用优化前8阶SCD和优化后6阶OSCD格式得到的地震记录，时间精度均为2阶。为了显示清晰，对地震记录进行了瞬时自动增益控制(AGC)处理，时窗为500ms。