王 婷，段焱文，蔡伟祥，易院平，汪 勇

(1. 长江大学 油气资源与勘探技术教育部重点实验室，湖北 武汉 430100； 2. 长江大学 地球物理与石油资源学院，湖北 武汉 430100； 3. 中石化重庆涪陵页岩气勘探开发有限公司，重庆 涪陵 408014; 4. 中冶集团武汉勘察研究院有限公司，湖北 武汉 430000)

0 前 言

波动方程数值求解法是一种首先离散网格化计算区域，然后通过数值求解描述地震波传播的波动方程，最终得到波动方程较高精度的近似解析解，从而模拟地震波传播的方法。有限差分法是一种快速有效求解波动方程的数值模拟方法，可以方便灵活地应用于复杂的地质模型，并且运动速度不受影响，边界处理简单。但是，当在一个波长内空间采样点较少时，就会产生严重的网格频散现象。谭未一等[1]采用规则网格有限差分方法模拟了粘弹VTI介质中的波场，董良国等[2-4]将交错网格与高阶差分法有效结合，求解了横向同性介质中一阶速度—应力弹性波方程，并对其稳定性进行了探究。霍凤斌等[5]对各向异性介质弹性波进行了高阶交错网格有限差分数值模拟， Saenger等[6-7]首先提出了旋转交错网格有限差分算法，旋转交错网格与常规交错网格有限差分不同，旋转交错网格有限差分是将所有物理量仅定义在两个不同的位置。若是将应力(或者质点振动速度)定义于离散网格单元的中心，那么就将质点振动速度(或者应力)定义于离散网格单元的顶点，并且可将所有弹性模量都定义于相同位置(与应力所在位置相同)。因此，在模拟非均匀介质时，弹性常数就不需要进行平均处理，在模拟各向异性介质时，弹性常数就不需要进行内插处理[8]，Saenger等对此已经进行了有力证明。严红勇等[9]随后在粘弹TTI介质中进行了旋转交错网格高阶有限差分数值模拟，进一步探究了旋转交错网格的特性。通量传输校正技术(Flux-corrected transport method，FCT)是一种在增加少量计算量的前提下，有效压制数值模拟中造成数值频散的技术。Book等[10]首次将流体动力学中的通量校正传输方法应用于声波方程求解，杨顶辉等[11]将通量校正传输方法与正演模拟相结合，应用到弹性波方程正演模拟。均取得了较好的效果。张省等[12]在VTI介质中将常规交错网格与FCT技术相结合进行了数值模拟，有效压制了数值频散。本文在前人研究的基础上，首次将通量传输校正技术与旋转交错网格相结合，实现了弹性波方程旋转交错网格高阶有限差分数值模拟，提高了数值模拟稳定性，并利用FCT技术有效压制了粗网格情况下造成的数值频散，提高了数值模拟精度。

1 方法原理

1.1 一阶速度—应力弹性波方程

在二维TI介质中，用速度—应力表示的弹性波方程(假设体力为零)[5]为：

(1)

其中，Vi为i(i=x或z)方向速度分量，τij为ij(ij=xx,zz或xz)方向应力，ρ为密度，Cij为介质的弹性常数。在各向同性情况下，C11=C33=λ+2μ，C13=λ，C44=μ。

1.2 时间上2M阶差分近似

在本文中，使用旋转交错网格法来数值求解一阶弹性波方程时，速度是在t+△t/2时刻，应力是在 时刻进行计算的，我们可以利用泰勒公式展开式来提高时间差分精度，得到2M阶精度的时间差分近似[4]：

(2)

其中，∂2m-1Vx/∂t2m-1为x方向速度对时间的2m-1阶偏导数。

将公式(2)中速度对时间的奇数阶高阶导数用公式(1)中应力对空间的导数表示，同理也可以将应力对时间的奇数阶高阶导数用速度对空间的导数表示。

1.3 旋转交错网格有限差分格式

在常规交错网格中，若将速度分量(或应力分量)定义在不同整网格点，应力分量(或速度分量)定义在不同半网格点，则此时弹性模量应该定义在所有的半网格点和整网格点上，而在实际情况中，通常只将弹性模量定义在半网格点上，故需要对弹性模量进行插值。当弹性模量变化很大时，普通交错网格就容易出现计算不稳定的问题。

在旋转交错网格中，若将应力(或者质点振动速度)位于离散网格单元的中心，质点振动速度(或者应力)位于离散网格单元的顶点，则此时所有弹性模量都可定义于相同网格上(与应力所在位置相同)。因此，在对非均匀介质进行数值模拟时，就不需要对弹性常数进行平均处理，在对各向异性介质进行数值模拟时，就不需要对弹性常数进行内插处理。图1对二维情况下的常规交错网格有限差分算法和旋转交错网格有限差分算法的空间位置进行了解释。

(a)常规交错网格 (b)旋转交错网格图1 二维网格单元

旋转交错网格有限差分算法分两步完成，首先是计算沿着网格对角线方向相关物理量的空间导数，然后线性叠加这些沿对角线方向计算得到的相关物理量结果，从而得到沿坐标轴方向的空间导数[7]。

将旋转交错网格沿坐标轴方向的空间导数与弹性波波动方程相结合，可以得到一阶弹性波旋转交错网格差分格式，对此，Saenger进行了公式推导[6-7]。

1.4 旋转交错网格与常规交错网格的稳定性条件

分析一阶速度—应力弹性波方程组的时间二阶，空间2N阶精度旋转交错网格[10]和常规交错网格[3]有限差分方程，可得到两者数值解的稳定性条件(见表1)。

表1稳定性条件

1.5 通量传输校正技术

在有限差分法中，数值频散包括空间频散和时间频散，但空间频散占其主要方面。旋转交错网格属于有限差分法，虽然与常规网格相比，数值频散得到了改善，精度得到了提高，但仍存在有限差分法固有的数值频散的问题，网格越大，数值频散越严重。所以，提高差分方程阶数来提高波场模拟精度和减小空间网格间距来压制数值频散是一种常规方法。但是，随着网格间距的缩小，运算效率会降低。通量传输校正技术(FCT)适用于高阶差分并有较高的计算精度,在采样间隔较大时，可以在增加少量计算量的前提下有效压制数值模拟中造成的数值频散。通量传输校正技术虽然能有效的压制数值频散，但存在较难选取合适参数的问题，从而影响校正效果。作者将采用改进后的FCT与旋转交错网格相结合，借助旋转交错网格来实现对弹性波方程的有限差分数值模拟，在增加少量计算机内存的情况下有效压制粗网格造成的数值频散。

FCT校正[12]是在假设所有的极值点都是由数值频散引起的前提下，对参与计算的所有网格点进行扩散通量校正处理，再对非局部极值点进行逆扩散通量校正来补偿损耗。FCT旋转交错网格有限差分方法的实现分两个阶段:旋转交错网格有限差分阶段和通量校正阶段。通量校正阶段又分为两步：平滑方程的解和反漫射矫正。

图2 旋转交错网格FCT

其具体步骤如下：

1)用旋转交错网格有限差分格式表示弹性波波动方程。求解k+1时刻和k时刻的速度分量。

2)以速度vx为例，计算k时刻的弥散通量P和Q。

3)运用k+1时刻计算得到的弥散通量对所有网格点(包括非局部极值点)上的速度vx进行修正。

4)反弥散通量的计算

a 计算k+1时刻的弥散通量P和Q；

c 进行反弥散通量处理来得到消除数值频散后的正确波场值。

反弥散通量的计算不是对所有的网格点进行处理，仅针对进行漫射通量校正前未产生频散的网格点(即非局部极值点)，通过修正被平滑过的非局部极值点方程的解，恢复不需要进行漫反射位置(未产生数值频散位置)的真实波场。具体计算公式可参考文献[12]。

2 模型试算

2.1 均匀介质VTI介质

设计一个均匀各向同性介质模型，其大小为2 000 m×2 000 m，纵波速度vp=3 000 m/s，横波速度vs=1 732 m/s，密度ρ=1 500 kg/m3，空间网格间距8 m和空间网格间距10 m，时间步长为0.5 ms，在模型(1 000 m，1 000 m)处使用震源信号主频为30 Hz的Ricker子波激发。采用精度为o(△t2+△x4)的旋转交错网格进行数值模拟，得到X分量在t=0.2 s时的波场快照，如图3所示，得到相应时间切片相同位置x分量记录振幅如图4所示。

由图4(a)和(b)可以看出，在网格大小8 m的情况下，P波波场无明显的数值频散，而S波由于速度较低，其波场模拟结果数值频散较为严重，存在明显的"拖尾"现象。在网格大小10 m的情况下，P波波场有轻微的数值频散，而S波波场数值频散比小网格严重。可见，空间网格较小的情况下，频散得到一定的改善。

比较图4(a)、(c)、(d)可得，利用通量传输校正技术后，数值频散(尤其是横波的数值频散)均得到了有效地压制，分辨率提高。能够清晰地反映通量校正传输技术的校正结果与漫射因子和反射因子有关，比较图4(c)、(d)、(f)，可得出压制效果(d)>(f)>(c)，分析可得漫射因子的取值没有特定标准，和数值频散的大小有关，当数值频散较为严重时，漫射因子取较大值时，频散压制效果较为出色。当漫射因子取值达到合适值时，频散压制效果不再随着漫射因子的变化而发生显著变化。

(a)△x=△z=8，η1=0，η2=0；(b)△x=△z=10，η1=0，η2=0；(c)△x=△z=8，η1=0.01，η2=0.02；(d)△x=△z=8，η1=0.02，η2=0.03；(e)△x=△z=10，η1=0.01，η2=0.02；(f)△x=△z=10，η1=0.02，η2=0.03

图3旋转交错网格不同参数波场快照对比

(a)△x=△z=8，η1=0，η2=0；(b)△x=△z=10，η1=0，η2=0；(c)△x=△z=8，η1=0.01，η2=0.02；(d)△x=△z=8，η1=0.02，η2=0.03；(e)△x=△z=10，η1=0.01，η2=0.02；(f)△x=△z=10，η1=0.02，η2=0.03

图4旋转交错网格不同参数振幅记录对比

为了进一步探究在得到相同计算精度的情况下使用FCT和不使用FCT的运行效率和资源占用情况，本文在未使用FCT情况下缩小空间网格间距，若要达到使用FCT(△x2=△z2=10，η1=0.02,η2=0.03)的精度(模型相同)，需要将网格空间大小取到4 m，得到波形对比图如图5所示，a实线为旋转交错网格应用FCT后的数值模拟结果，b实线为旋转交错网格数值模拟结果，c虚线为常规交错网格的精确解析解。

图5 3种方法振幅记录对比

从图5可以看出，三者基本一致，旋转交错网格

应用FCT技术的模拟结果略优于旋转交错网格的数值模拟，旋转交错网格略优于常规交错网格。

资源占用量和效率如表2所示。从表2中数组个数和数组大小可以得到，使用FCT所用的存储空间(数组个数×数组大小)比未使用FCT的存储空间小，达到相同效果所需计算机运行时间也略小。因此，我们可以看到使用FCT技术确实提高了效率和节省了存储空间。

表2 运行效率和资源占用对比

2.2 水平层状介质模型

设计了水平层状介质模型(两种介质分界面位于深度1 000 m处)，模型总大小为2 000 m×2 000 m(vp1=3 000 m/s，vs1=1 732 m/s，ρ1=2 000 kg/m3，vp2=4 200 m/s，vs2=2 425 m/s，ρ2=2 500 kg/m3)，空间网格间距△x=△z=8 m，时间步长为1 ms，在模型(1 000 m，1 000 m)处使用震源信号主频为30 Hz的Ricker子波激发，得到速度X分量在t=0.5 s时的波场快照，如图6(a)所示，取漫射因子η1为0.02，反漫射因子η2为0.03，波场快照如图6(b)所示。

①为第一层介质入射S波；②为第一层介质入射P波；③为第一层介质PP透射波；④为第一层介质PP反射波；⑤为第一层介质PS透射波；⑥为第一层介质PS转换波

图6水平层状介质波场快照

抽取第125道记录振幅进行对比，如图7所示。

图7 第125道记录振幅对比(△x1=△z1=8)

观察图5和图6可得，在层状介质模型中进行正演模拟，通量传输校正技术也可以有效压制数值频散，但由于通量校正传输方法进行扩散通量校正处理是对所有网格点进行处理，造成了非局部极值点的损失，随后再对有损失的非局部极值点进行逆散通量校正处理来进行补偿损失，所以也造成了一定的振幅偏差。图7可以十分清晰的看出使用FCT后的效果得到显著提高。

3 结 论

旋转交错网格有限差分数值模拟方法是一种有效的地震波场数值模拟方法，可以将地震波在地下介质中的传播过程及规律清晰明了地给呈现出来，模拟精度较高，稳定性相较于常规网格得到了加强，实现过程较简单容易。通量传输校正技术可以有效压制数值模拟过程中由于网格差分产生的数值频散，特别是在大网格情况下很大地提高了计算精度，清晰地呈现出波在地下介质中传播过程，反射因子和漫反射因子的取值没有唯一标准，取值大小直接影响通量传输校正技术的效果，当频散较为严重(网格取值较大)时，漫射因子可以取较大值，当产生轻微频散时，可以取较小值。当取值的变化不再影响频散的压制效果时，说明取值合适。本文首次将通量传输校正技术与旋转交错网格相结合，对均匀介质和层状介质模型进行正演模拟，不仅有效压制数值频散，达到较高模拟精度，还减少了增加网格点来提高计算精度所产生的计算量，提高了运算速度。