扩展的近似解析离散化方法及弹性波方程数值模拟

2017-11-01段焱文安一凡王笑丛桂志先

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

关键词：步长差分解析

汪勇段焱文安一凡王笑丛桂志先

(①长江大学地球物理与石油资源学院，湖北武汉 430100； ②油气资源与勘探技术教育部重点实验室(长江大学)，湖北武汉 430100； ③东方地球物理公司国际勘探事业部，河北涿州 072751)

扩展的近似解析离散化方法及弹性波方程数值模拟

汪勇*①②段焱文①②安一凡①②王笑丛③桂志先①②

汪勇，段焱文，安一凡，王笑丛，桂志先.扩展的近似解析离散化方法及弹性波方程数值模拟.石油地球物理勘探，2017,52(5):928-940,955.

地震波场数值模拟是研究地震波理论、偏移成像和地震反演等工作的基础，提高数值模拟的精度具有重要意义。在前人的研究基础上，提出了一种扩展的近似解析离散化数值模拟方法，并从理论上对该方法的精度、数值频散、稳定性和计算效率进行了分析。该方法在时间差分上比扩展前提高了一阶精度、误差最大降低了88%。与其他近似解析离散化类方法一样，具有算子半径小和适应粗网格步长的优势，在最小主波长内仅需使用5.9个网格点。与四阶Lax-Wendroff修正格式和交错网格有限差分格式的频散曲线和模拟结果对比，验证了该方法能更好地压制数值频散。二维各向同性均匀介质和水平层状介质模型数值模拟的地震弹性波场特征清晰准确，说明了该方法的实用性。

扩展最优近似解析离散化弹性波方程数值模拟稳定性条件数值频散

1 引言

地震波场数值模拟是由已知的岩石结构和物理参数建立地震地质模型，依据不同的算法模拟地震波在地下介质中的传播过程，计算在地面或地下各观测点地震记录的一种方法。目前，常用的数值模拟方法主要有射线追踪法[1]和波动方程法，其中波动方程法有伪谱法[2，3]、有限元法[4]、边界元法[5]、谱元法[6]和有限差分法。20世纪 60 年代，Alterman等[7]首先将有限差分方法用于地震波场的数值模拟，该方法现已成为地震波场数值模拟中应用最为广泛的一种方法，差分格式也由早期的低阶差分发展到高阶差分，由常规网格发展到交错网格、旋转交错网格、可变网格等[8-16]。近似解析离散化方法(Nearly-Analytic Discrete Method，简称NADM)是20世纪90年代发展起来的一种有限差分数值模拟方法[17，18]，其基本思想是在时间上采用泰勒公式展开，在空间上利用截断的泰勒公式构造高阶插值函数来逼近空间偏导数，与其他有限差分法的根本区别是它利用空间节点的位移和梯度共同逼近变量的空间偏导数，从而提高计算精度，在大网格下可有效地压制数值频散，提高计算效率。杨顶辉及其研究团队对该方法进行了深入的研究[19-25]。

本文在前人研究的基础上，对优化的近似解析离散化方法进行了扩展，即将位移场的泰勒展开由5项扩展为6项，将时间差分精度提高到4阶，进而提高了数值模拟精度。

2 扩展的近似解析离散方法原理

近似解析离散化类数值模拟方法的基本思想是：在时间上采用泰勒公式展开，在空间上利用截断的泰勒公式构造高阶插值函数来逼近空间偏导数，从而得到原函数的一个最优近似。该方法的特点是在求解偏微分方程时包含质点位移和速度的各阶空间偏导数，应用了原函数、各阶偏导数之间的相互联系，能有效减少离散过程中地震信息的丢失、提高数值计算的精度。以二维声波方程为例讨论该类方法数值模拟的基本思想。

二维声波波动方程为

(1)

设满足声波方程的解为U，定义

(2)

(3)

式中m、k分别表示S对x和t的偏导数的阶次。

利用截断的泰勒公式表示n+1时刻的U，可以得到

(4)

式中：i和j分别为x和z方向空间网格点序号；n为时间序号。P及其时间的偏导数可以由波动方程式(1)转化为U和S对空间的高阶偏导，即

(5)

(6)

(7)

根据式(3)，如果式(6)中S对空间的二阶偏导数可用U表示的一阶向后时间差分近似，即

(8)

则称为原始的近似解析离散化方法。

(9)

(10)

文献[24]指出，由于采用式(8)表示的一阶向后差分格式求取S，降低了整体数值模拟精度，导致了NADM方法精度仅为O(Δt2+Δx4+Δz4)。

为了提高时间层推进精度，对S同样用泰勒公式展开，则可以得到速度场S的时间层推进差分格式

(11)

式(4)～式(7)以及式(11)构成了改进的近似解析离散化方法(Improving Nearly-Analytic Discrete Method，简称INADM)的差分格式。该方法与原方法相比，在增加了存储量和计算量的情况下，提高了数模拟值时间精度。

此外，如果利用截断的泰勒公式表示n-1时刻的U，可以得到

(12)

将式(4)和式(12)相加，可以得到最优近似解析离散化方法(Optimum Nearly-Analytic Discrete Method，简称ONADM)的差分格式

(13)

ONADM方法的核心是消除了差分格式式(4)中的速度项，避免了原方法中降低时间精度的因素，从而提高了整体差分精度。文献[21]指出ONADM的精度为O(Δt2+Δx4+Δz4)。与INADM方法相比，该方法由于不涉及速度项，所以在计算效率和存储空间方面得到了提高，但其缺点是不能进行带速度项的波动方程数值模拟，如黏弹介质波动方程模拟等。

INADM和ONADM方法有各自的优缺点，ONADM的精度比INADM高，但INADM具有更广泛的适用性，能够对黏弹介质等复杂介质波动方程进行数值模拟。所以本文选择对INADM方法进行修正或扩展，提高INADM中位移和速度项的差分精度是本文的研究目的。

(14)

因为按照INADM方法计算时，原本就需要计算速度场S，所以式(14)中的速度高阶空间导数(≤5次)是容易直接求取的，这就为提高泰勒展开精度提供了基础，进而U和S的时间层推进公式可以扩展为

(15)

(16)

式(15)和式(16)即为扩展的近似解析离散化方法(Expand Nearly-Analytic Discrete Method，简称ENADM)的差分公式。式中P的各阶时间导数分别用式(5)～式(7)和式(14)计算，其中涉及到的空间高阶偏导数则与其他NAD类方法一样，通过质点的位移场和速度场及其梯度逼近计算。在上述差分公式推导中，虽然介质的地震波速度v是常数，但在数值模拟中可以随空间位置而改变，即上述方法依然可以用于层状介质和非均匀介质等复杂介质的数值模拟研究。

3 ENADM方法分析

3.1 精度分析

NAD类有限差分法中的空间导数是利用待求空间网格点周围点的位移及梯度共同逼近， Yang等[19]指出，它们的空间精度为四阶，所以ONADM、INADM和ENADM的空间精度相同，在后面的频散关系讨论中，也可以看出其频散曲线基本一致。ONADM采用式(13)进行位移的时间层推进，由于消掉了速度项，所以它的时间精度为四阶。INADM的位移和速度分别采用式(4)和式(11)进行时间层推进，位移的时间精度为四阶，而速度为三阶，所以整体时间精度应为三阶。ENADM在INADM基础上进行了扩展，其位移和速度按式(15)和式(16)分别进行计算，所以位移和速度的时间精度分别为五阶和四阶，整体的时间精度为四阶，比INADM提高一阶时间精度，所以ENADM的差分精度为O(Δt4+Δx4+Δz4)。

应用二维平面谐波初值问题，比较四阶Lax-Wendroff修正格式(简称LWC)、最优近似解析离散化方法(ONADM)、改进的近似解析离散化方法(INADM)和扩展的近似解析离散化方法(ENADM)数值模拟精度。二维平面波初值问题可以表示为

(17)

式中：θ0是初始时刻波阵面法线方向(即传播方向)与x轴的夹角；f0是平面简谐波的频率。其准确解析解为

(18)

LWC差分格式[26]为

(19)

二维波场数值模拟中，假设平面波频率f0=20Hz,θ0=π/4。设置均匀介质模型，波速为4000m/s，模型长度和深度均为2000m，纵、横向网格间距相同。在不同空间步长和时间步长条件下，计算前述四种数值模拟结果的相对误差。相对误差定义为

(20)

表1 四种方法在不同情况下的最大相对误差(%)比较

图1 四种差分方法的相对误差曲线(a)Δx=10m，Δt=0.5ms； (b)Δx=15m，Δt=0.5ms； (c)Δx=20m，Δt=1.0ms； (d)Δx=25m，Δt=2.0ms

从表1和图1可以看出，四种数值模拟方法相对误差均随时间步长或空间步长的增大而增加。从误差统计来看，ENADM具有最高的数值精度，INADM和ONADM次之，而LWC的数值精度最低、误差增长最快。当空间步长为25m，时间步长为2ms时，ENADM最大相对误差仅为0.45%，而扩展前的INADM却有3.74%，其相对误差最大降低了88%，能更好地用于大尺度模型的长时间数值模拟。

3.2 频散分析

数值模拟中由于求解的误差而导致数值波速与频率有关，即产生所谓的数值频散现象，这不是地震波特征真实的反映，所以频散关系分析是判断一种数值模拟方法优劣的重要依据，其结果决定了该方法能否或在何种条件下适用于地震波正反演计算。选取LWC、ONADM、INADM、ENADM和交错网格有限差分(Staggered-grid Finite Difference Method，简称SFDM)五种方法进行二维频散关系的比较和分析，以说明ENADM方法在压制数值频散方面的优越性。

ENADM方法频散关系(推导过程见附录A)为

φxsinφx+4(cosφz-1)+φzsinφz]-

3[φzsinφz+2(cosφz-1)]-2(cosφxcosφz-

cosφx-cosφz+1)}

(21)

确定φ后，代入上述频散关系，可以解得对应的ω。定义数值波速与真实速度的比值为

(22)

在理想情况下，如果不存在数值频散则速度比γ恒等于1。γ偏离1越大，说明该方法的数值频散越严重，反之则说明该方法能更好地压制数值频散。

同理可以得到其他四种差分方法的频散关系。INADM频散关系为

φxsinφx+4(cosφz-1)+φzsinφz]-

(23)

ONADM频散关系为

exp(-iω)=2-exp(iω)+α2[4cosφx+Δxkxsinφx+

4cosφz+Δxkzsinφz-8]+

(24)

LWC频散关系为

cosω-1

(25)

SFDM频散关系[27]为

(26)

图2为五种方法在不同的α和θ时的频散关系曲线。

图2 五种差分方法在不同参数时的频散曲线

取φ∈[0,π]作为横坐标，它是波数与空间步长的乘积。当空间步长确定后，φ的增加表示频率随之增加。另一方面，单位波长内采样点数M=2π/φ，所以横坐标也可以看作M由∞逐渐减小至2。从图中的频散曲线可以看出：①随着空间采样点数的减小，5种方法的频散现象逐步加剧，而NAD类方法的数值频散比LWC和SFDM方法要小，其频散曲线更趋近于1，说明了这类方法在压制数值频散方面的更具优势；②假设数值速度在理论速度的99.9%以内表示不存在频散，则NAD类、SFDM和LWC对应的最小的φ分别为0.337、0.274和0.131，所以它们在最小主波长内需要使用的网格点数分别为5.9、7.3和15.3个；③由于三种NAD类方法在求取空间偏导数时所用差分格式和精度相同，所以它们的频散曲线基本重合。在θ=π/4,α=0.45时，ENADM方法的速度比最接近1，说明该方法在大Courant数和对角线方向传播时，它能更好地压制频散；④不同的θ，即波在不同的传播方向上，速度比不同，说明存在一定的数值各向异性，α越大越明显。当α=0.45时，ENADM、INADM、ONADM、LWC和SFDM在三个方向上的速度比最大相差分别为8.2%、8.4%、8.3%、20.7%和9.2%，说明NAD类方法对数值各向异性也有很好的压制，且ENADM效果最好。

图3 二维声波均匀模型不同方法800 ms时刻模拟的波场快照

图4比较了三种方法在(x=300m，z=1500m)处接收到的波形记录。图4a中的黑色实线为ENADM数值模拟结果，红色虚线为精确解析解，二者基本重合，模拟结果非常理想。与同是四阶精度的LWC和SFDM相比，ENADM的质点振动曲线没有LWC和SFDM中的“拖尾”现象，在粗网格数值模拟时具有更高的计算精度。需要说明的是，数值模拟中三种方法均采用了相同的时间和空间网格，但SFDM采用一阶速度—应力方程，模拟结果是应力分量地震记录； ENADM和LWC采用二阶位移方程，模拟结果是位移地震记录，导致了三种方法所得地震记录振幅在数量级上存在一定的差别。为了比较，图4中的三个震动曲线经过了振幅均一化处理。

图4 在(x=300m，z=1500m)处三种方法波形记录

3.3 稳定性分析

稳定性条件是有限差分数值模拟中一个非常重要的问题，如果不满足稳定性条件，则前面各层的舍入误差将会影响到后面各层的值，后面的误差的影响越来越大，会使有限差分的结果完全错误。本文采用Fourier方法[28]对ENADM的差分格式进行稳定性分析。记

(27)

(28)

式中vmax为波的最大传播速度。

3.4 计算效率分析

对上面的二维均匀模型进行声波数值模拟，采样时间为1s。LWC和三种NAD类方法使用的空间网格间距Δx和时间网格间距Δt的选择满足无数值频散和稳定性的要求。四种方法的计算效率、使用的数组大小、个数如表2所示。

表2 不同方法计算效率比较对比

从表2可看出，由于ENADM算法提高了时间精度，导致计算时所需要的内存和时间增加，与扩展前的INADM相比，其占用内存增加了28.6%，计算时间增加了28.4%，即ENADM方法在增加存储和计算量的条件下提高了数值模拟的精度。同时，因为计算了速度分量，所以很容易进行含速度项的波动方程数值模拟研究，在实际数值模拟中，可以根据介质类型和精度要求，选择合适的数值模拟方法。

4 模型试算

采用ENADM方法对二维各向同性完全弹性介质的弹性波方程进行数值模拟，波动方程为

(29)

式中：vP和vS分别为纵、横波波速；ux和uz分别为弹性波在x和z方向的位移，对其进行泰勒公式展开可得

(30)

4.1 均匀介质模型

均匀模型介质为泊松体，尺寸为2000m×2000m，纵波速度为4000m/s，横波速度为2310m/s。由于最小速度为2310m/s，按频散分析的结果，每个波长内需采样5.9个点，可取的最大空间步长为9.8m，为确保精度，设定纵、横向空间间隔均为8m。由于最大速度为4000m/s，按稳定性条件，设置时间步长为0.8ms。在模型中心处激发40Hz的Ricker子波等能量震源。图5为220ms时刻的波场快照，波场非常清晰，没有出现数值频散现象。图5中外层为纵波波前，内层为横波波前。由图可以看出，采用等能量震源时：①图5a显示的位移水平分量上，当x为1000m时，由于纵波质点垂向振动，横波质点横向振动，所以纵波位移为0，横波最大，且纵波波前相位以x=1000m左右对称。当z=1000m时，则纵波最大，而横波为0，且横波波前相位以z=1000m上下对称。②图5b显示的位移垂直分量上，纵波和横波波前以x=1000m和z=1000m为对称轴，相位相反。

图5 均匀介质模型弹性波方程ENADM数值模拟波场快照

4.2 水平层状介质模型

设置两层水平层状介质模型每层厚度均为750m，模型长度为1500m。第一层纵、横波速度分别为3800m/s和2600m/s；第二层纵、横波速度分别为4500m/s和2900m/s。在地面中间处(x=750m，z=0m)激发20Hz雷克子波震源，震源加载到x和z方向位移上。纵、横向空间步长均为15m，时间步长为1ms，采用PML人工边界条件[29，30]。图6为地面接收到的地震记录，长度为1s，记录的纵横波波场清晰，没有数值频散。图中ZP和ZS分别表示直达纵波和直达横波，PP表示反射纵波，SS表示反射横波，PS和SP表示在反射界面形成的转换波，它们达到检波点的时间相同，相互干涉叠加。

总体来看，双层模型中地面模拟接收到的地震记录非常清晰，没有数值频散和不稳定数值结果，地震记录中的直达波、反射波和转换波显示清楚，说明ENADM算法可以有效地模拟弹性波在多层各向同性介质模型中的传播。

图6 水平层状介质模型弹性波方程ENADM数值模拟地面地震记录

5 结论与建议

本文在前人的研究基础上，提出并推导了二维各向同性介质弹性波方程的扩展近似解析离散化数值模拟方法。通过理论分析和模型试算，进行了该方法的模拟精度、频散关系分析、稳定性研究及计算效率分析。研究结果表明，扩展的近似解析离散化数值模拟方法具有模拟精度高、适用粗网格计算的优势，在大尺度模型波场模拟时具有明显的优势，在研究深部地球物理和石油地球物理勘探等方面有广泛的应用前景。本文仅将该方法应用于二维各向同性介质弹性波模拟，但该方法能进一步推广到三维、各向异性、黏弹及双相介质等复杂介质的模拟，也可以采用同样的思路将其他NAD类方法(如加权近似解析离散)进行扩展，为研究地震波在复杂介质中的传播规律提供了一种新的方法和思路。

[1] Cerveny V.Seismic rays and ray intensities in inhomogeneous anisotropic media.Geophyscial Journal International，1972，29(1)：1-33.

[2] Kosloff D D，Baysal E.Forward modeling by a Fourier method.Geophysics，1982，47(10)：1402-1412.

[3] 唐怀谷，何兵寿.一阶声波方程时间四阶精度差分格式的伪谱法求解.石油地球物理勘探，2017，52(1)：71-80. Tang Huaigu，He Bingshou.Pseudo spectrum method of first-order acoustic wave equation finite-difference schemes with fourth-order time difference accuracy.OGP，2017，52(1)：71-80.

[4] Drake L A.Rayleigh waves at a continental boundary by the finite element method.Bulletin of the Seismological Society of America，1972，62(5)：1259-1268.

[5] Bouchon M，Coutant O.Calculation of synthetic seismograms in a laterally varying medium by the bounda-ry element-discrete wavenumber method.Bulletin of the Seismological Society of America，1994，84(6)：1869-1881.

[6] Komatissch D，Barnes C，Tromp J.Simulation of anisotropic wave propagation based upon a spectral element method.Geophysics，2000，65(4)：1251-1260.

[7] Alterman Z，Karal F C.Propagation of seismic wave in layered media by finite difference methods.Bulletin of the Seismological Society of America，1968，58(1)：367-398.

[8] Kelly K R，Ward R W，Treitel S et al.Synthetic seismograms：a finite-difference approach.Geophysics，1976，41(1)：2-27.

[9] Virieux J.P-SV wave propagation in heterogeneous media：Velocity-stress finite-difference method.Geo-physics，1986，51(4)：889-901.

[10] Graves R W.Simulation seismic wave propagation in 3D elastic media using staggered-grid finite differences.Bulletin of the Seismological Society of America，1996，86(4)：1091-1106.

[11] 董良国，马在田，曹景忠等.一阶弹性波方程交错网格高阶差分解法.地球物理学报，2000，43(3)：411-419. Dong Liangguo，Ma Zaitian，Cao Jingzhong et al.The staggered-grid high-order difference method of first-order elastic equation.Chinese Journal of Geophysics，2000，43(3)：411-419.

[12] 董良国，马在田，曹景忠.一阶弹性波方程交错网格高阶差分解法稳定性研究.地球物理学报，2000，43(6)：856-864. Dong Liangguo，Ma Zaitian，Cao Jingzhong.The stability study of the staggered-grid high-order difference method of first-order elastic equation.Chinese Journal of Geophysics，2000，43(6)：856-864.

[13] 董良国.弹性波数值模拟中的吸收边界条件.石油地球物理勘探，1999，34(1)：46-56. Dong Liangguo.Absorptive boundary condition in elastic-wave numerical modeling.OGP，1999，34(1)：46-56.

[14] 王书强，杨顶辉，杨宽德.弹性波方程的紧致差分方法.清华大学学报(自然科学版)，2002，42(8)：1128-1131. Wang Shuqiang，Yang Dinghui，Yang Kuande.Com-pact finite difference scheme for elastic equations.Journal of Tsinghua University(Science and Techno-logy)，2002，42(8)：1128-1131.

[15] Saenger E H，Gold N，Shapiro S A.Modeling the propa-gation of elastic waves using a modified finite-difference grid.Wave Motion，2000，31(1)：77-92.

[16] 徐文才，杨国权，李振春等.横向各向同性介质拟声波一阶速度—应力方程.石油地球物理勘探，2016，51(1)：87-96. Xu Wencai，Yang Guoquan，Li Zhenchun et al.First order velocity-stress equation in TI media.OGP，2016，51(1)：87-96.

[17] Kondoh Y，Hosaka Y，Ishii K.Kerrel optimum nearly analytical discrimination algorithm applied to parabo-lic and hyperbolic equations.Computers & Mathema-tics with Applications，1994，27(3)：59-90.

[18] 杨顶辉，滕吉文，张中杰.三分量地震波场的近似解析离散模拟技术.地球物理学报，1996，39(增刊)：283-291. Yang Dinghui，Teng Jiwen，Zhang Zhongjie.Nearly-analytical discretization modeling technique of 3-component seismic wave-fields.Chinese Journal of Geophysics，1996，39(S)：283-291.

[19] Yang D H，Teng J W，Zhang Z J et al.A nearly-analytic discrete method for acoustic and elastic wave equation.Bulletin of the Seismological Society of America，2003，93(2)：882-890.

[20] Yang D H，Lu M，Wu R S et al.An optimal nearly-ana-lytic discrete method for 2D acoustic and elastic wave equations.Bulletin of the Seismological Society of America，2004，94(5)：1982-1991.

[21] Yang D H，Song G J，Lu M.Optimally accurate nearly analytic discrete scheme for wave-field simulation in 3D anisotropic media.Bulletin of the Seismological Society of America，2007，97(5)：1557-1569.

[22] Yang D H，Song G J，Chen S et al.An improved nearly analytical discrete method：an efficient tool to simulate the seismic response of 2-D porous structures.Journal of Geophysics & Engineering，2007，4(1)：40-52.

[23] 卢明.改进的近似解析离散化方法及弹性波波场模拟[学位论文].北京：清华大学数学科学系，2004.

[24] 宋国杰.三维弹性波方程的改进近似解析离散化方法及波场模拟[学位论文].北京：清华大学数学科学系，2008.

[25] Tong P，Yang D H，Hua B L.High accuracy wave simulation-revised derivation，numerical analysis and testing of a nearly analytic integration discrete method for solving acoustic equation.International Journal of Solid and Structures，2011，48(1)：56-70.

[26] Lukacova M M，Warnecke G.Lax-Wendroff type se-cond order evolution Galerkin methods for multidimensional hyperbolic systems.East-West Journal of Numerical Mathematics，2000，8(2)：127-152.

[27] 董良国，李培明.地震波传播数值模拟中的频散问题.天然气工业，2004，24(6)：53-56. Dong Liangguo，Li Peiming.Dispersive problem in seismic wave propagation numerical modeling.Natural Gas Industry，2004，24(6)：53-56.

[28] 张文生.科学计算中的偏微分方程有限差分法.北京：高等教育出版社，2006.

[29] 高正辉，孙建国，孙章庆等.基于完全匹配层构建新方法的2.5维声波近似方程数值模拟.石油地球物理勘探，2016，51(6)：1128-1133. Gao Zhenghui，Sun Jianguo，Sun Zhangqing et al.2.5D acoustic approximate equation numerical simulation with a new construction method of the perfectly matched layer.OGP，2016，51(6)：1128-1133.

[30] 汪勇，段焱文，王婷等.优化近似解析离散化方法的二维弹性波波场分离模拟.石油地球物理勘探，2017，52(3)：458-467. Wang Yong，Duan Yanwen，Wang Ting et al.Numerical simulation of elastic wave separation in 2D isotropic medium with the optimal nearly-analytic discretization.OGP，2017，52(3)：458-467.

附录AENADM频散关系

将式(5)～式(7)和式(4)代入式(15),可以得到ENADM的差分格式为