张朝元

(大理大学 数学与计算机学院， 云南 大理 671003)



求解弹性波方程的高精度ONAD方法及其波场模拟

张朝元

(大理大学 数学与计算机学院， 云南 大理 671003)

基于弹性波传播方程，发展了一种高精度低数值频散的八阶ONAD(optimal nearly-analytic discrete)方法，该方法利用八阶精度的近似解析离散算子对空间高阶偏导数进行离散，采用四阶精度的截断泰勒展开式离散时间高阶导数。八阶ONAD方法被用于模拟地震波在VTI介质模型和2个复杂层状介质模型中的传播。计算效率结果表明，该方法在运算速度和存储量上明显优越于八阶LWC方法。波场模拟结果显示，八阶ONAD方法在粗网格条件下可有效消除由速度强间断所造成的数值频散，有利于在强间断介质中使用粗网格进行波场模拟，是一种在地震勘探领域有着巨大应用潜力的数值方法。

弹性波方程；ONAD方法；数值频散；波场模拟；高精度

众所周知，正演模拟是地球物理学中开展地震波全波形反演的一个关键步骤和核心环节。如果正演数值模拟方法过多地消耗和占用计算资源，如占用计算机内存空间过多和计算运行速度过慢等，则必定会降低全波形反演速度，增加地震反演成本，严重时会导致地震反演不能成功实现。因此，研究效果好、精度高的正演地震模拟方法，成为地震波方程全波反演研究工作中非常关键的问题。

在已经发展的地震数值模拟方法[1-8]中，有限差分法因其计算机编程简单易实现且并行性强、同网格数下相比存储量小和计算速度快等优点，被有效地应用于地震数值模拟中。但传统中的有限差分方法在数值离散地震波方程时，由于数值波速严重受波数的非线性函数影响与真实波速有差异，因此在层间速度反差大或大尺度粗网格时，进行地震波数值模拟会产生严重的数值频散。为了避免这个问题，Yang等人于2003年率先通过引入近似解析离散(NAD)算子替代传统的差分算子来对波动方程中的空间高阶偏导数进行离散化[9]，从而发展了一系列具有低数值频散的数值方法[9-17]。但是这些基于NAD算子的数值方法在空间离散上仅为四阶精度，还没有更高精度的研究。为了进一步提高数值模拟精度，本文在Yang等人[10]研究基础上，基于弹性波方程，在空间上采用Tong等人发展的八阶精度NAD算子对高阶偏导数进行离散，在时间上采用四阶截断的泰勒展开式对高阶导数进行离散，得到时间精度为四阶、空间精度为八阶的ONAD方法[17]，称为八阶ONAD方法。

我们将该方法应用到弹性波在垂直对称的单层横向各向同性介质、双层均匀各向同性介质和三层均匀横向各向同性介质中的传播模拟中。数值结果表明，同具有相同精度的八阶LWC方法[18,19][精度为O(Δt4+Δx8+Δz8)]相比，本文提出的八阶ONAD方法具有计算速度快和存储量小等优势，可有效消除由粗网格条件或层状速度反差较大所引起的数值频散，能在强间断介质中很好地进行粗网格或大尺度条件下的地震波场数值模拟。

1 八阶ONAD方法的计算公式

在Yang等人[10]的研究基础上，我们发展了一种求解弹性波方程的八阶ONAD方法。设二维弹性波传播方程如下

(1)

其中:下标j=1,3;σij和ρ=ρ(x,z)分别代表应力和介质密度;fi和ui分别为i方向的震源函数分量和位移分量。在后面的波场模拟实验中，我们选择随时间变化的震源函数为

exp[-8(0.6f0t-1)2]

(2)

其中:f(t)为地震波震源主频f0的Ricker子波;f0是震源函数的频率。

(3)

其中D是三阶微分算子。

(4)

(5)

(6)

将(5)式和(6)式相加，可得

(7)

(8)

我们把用于计算二维弹性波方程(1)式的逼近方程(8)称为八阶ONAD方法。

对于八阶ONAD方法，由截断的泰勒展开式得到时间离散方程(8)，其时间误差精度为O(Δt4)；由插值函数和泰勒展开式离散表达式(8)右边的高阶空间偏导数，其逼近公式的空间误差精度[17]为O(Δx8+Δz8)。因此，二维弹性波方程求解的八阶ONAD方法的误差精度为O(Δt4+Δx8+Δz8)。

2 计算效率分析

衡量地震波模拟方法优越的重要标准就是考察该方法的计算效率，而计算效率主要包括计算耗时和计算存储量。为研究八阶ONAD方法的计算效率和模拟精度，在这里我们分别用八阶ONAD方法和具有相同精度的八阶LWC方法[精度为O(Δt4+Δx8+Δz8)]在粗网格条件下对具有垂直对称轴的横向各向同性介质(记为VTI)模型进行波场模拟，同时与八阶LWC方法在细网格下模拟的结果进行比较。

该模型的弹性系数和介质密度分别为c11=45 GPa,c13=9.6 GPa,c33=37.5 GPa,c44=12 GPa和ρ=1.0 g/cm3。模型的模拟区域为0≤(x,z)≤12 km；采用频率f0=30 Hz的Ricker震源子波函数，f1=f3=f(t)，震源位于模拟区域中央。设定时间步长为Δt=0.8 ms，空间步长为Δx=Δz=40 m,接收器位于R(6 km, 5 km)。

图1是粗网格(Δx=Δz=40 m)下分别由八阶ONAD和八阶LWC两种方法在t=0.6 s时刻产生的水平位移方向的瞬时波场快照。从图1可知，八阶ONAD方法产生的快照图(图1-A)非常清晰，无可见数值频散，可观测到清晰的尖点与三叉区，各种波形清晰分明；而八阶LWC方法却存在明显的数值频散(图1-B)，尤其在三叉区和尖点等处存在严重的数值频散。在相同Courant数下，当取空间步长为Δx=Δz=18 m时，八阶LWC方法此时能有效地消除数值频散。但在计算机环境Intel(R) Core(TM) 2Duo CPU和2G内存下，该两种方法为抑制数值频散所花费的存储量和CPU时间却明显不同。在存储量上，八阶LWC方法需要6个存储规模为667×667，但是八阶ONAD方法需要18个存储规模仅为300×300，这说明八阶ONAD方法需要的存储量是八阶LWC方法的60.69%。在CPU时间上，八阶LWC方法为消除数值频散所消耗的CPU时间为784 s，而用八阶ONAD方法模拟得到图1-A所消耗的时间仅为314 s，这表示八阶ONAD方法的计算速度约为八阶LWC方法的2.5倍。

我们以八阶LWC方法在细网格条件(Δx=Δz=18 m)下模拟所得到的结果作为精确解。图2分别给出了粗网格条件(Δx=Δz=40 m)下由八阶ONAD方法(图2-A)和八阶LWC方法(图2-B)模拟得到在接收器R(6 km, 5 km)处的波形记录，同时分别与精确解进行比较。从图2可知，八阶ONAD方法与精确解的波形基本一致，而八阶LWC方法与精确解却有着明显的误差。这意味着，同八阶LWC方法相比，八阶ONAD方法在粗网格条件下能够提供更加精确的数值解。因此，八阶ONAD方法有望在大尺度规模地震波模拟中得到推广和应用。

图1 在t=0.6 s时刻由八阶ONAD方法和八阶LWC方法产生的水平位移分量的波场快照Fig.1 Snapshots of wave fields at time t=0.6 s on the horizontal displacement component

图2 分别由精确解与八阶ONAD方法和八阶LWC方法计算得到在接收器R(6 km, 5 km)处的水平位移分量的波形记录比较Fig.2 Comparison of waveforms recorded at the receiver R(6 km, 5 km) on the horizontal displacement component by the analytic solutions

3 波场模拟

为考察本文发展的八阶ONAD方法模拟弹性波传播的效果及压制数值频散的能力，我们利用该方法来模拟弹性波在双层均匀各向同性介质(记为HI)和三层均匀横向各向同性介质(记为TI)中的传播。为体现本文提出的八阶ONAD方法的优越性，我们将数值模拟结果同具有相同精度的八阶LWC方法进行比较。

3.1 双层HI介质模型

研究地震波在层状介质模型的传播规律有利于人们对地球深层内部结构进行合理有效的了解，从而，有效模拟弹性波在地球内部介质分层界面的折射反射等具有深刻意义。其中，双层介质模型是层状介质模型中最重要的又是最简单的模型。为了考察八阶ONAD方法模拟地震波在强间断介质传播的有效性，我们进行一个强间断双层HI介质模型的数值模拟。

实验中，设0≤(x,z)≤4 km为双层HI介质模型的计算区域，z0=2.4 km为分层界面。上层介质中的弹性系数和介质密度分别为λ1=1.5 GPa,μ1=2.5 GPa,ρ1=1.5 g/cm3；下层介质中的弹性系数和介质密度分别为λ2=11.0 GPa,μ2=15.0 GPa,ρ2=2.0 g/cm3。设时间步长为Δt=1 ms，空间步长为Δx=Δz=20 m。震源位于O(2 km, 1.85 km)，震源函数为f0=22 Hz的Ricker子波，且f1(x)=f3(x)=f(t)。一个SV波的波长内包含3.2个空间采样点。

图3是八阶ONAD方法模拟获得t=0.6 s时刻在水平位移分量(图3-A)和垂直位移分量(图3-B)的瞬时波场快照。图4是八阶LWC方法模拟获得t=0.6 s时刻在水平位移分量(图4-A)和垂直位移分量(图4-B)的瞬时波场快照。从图3与图4比较可看出，八阶LWC方法取得的模拟波场快照在间断面邻域和速度较低的上层有较为严重的数值频散，而同时刻下八阶ONAD方法模拟取得的波场快照则非常清晰，没有可见的数值频散。这表明八阶ONAD方法在大尺度粗网格条件下能很好地抑制因强间断速度所导致的数值频散，有利于在强间断介质中使用大尺度粗网格进行地震波场数值模拟，从而提高模拟精度。因此，八阶ONAD方法对于二维弹性波在强间断层状介质中有着良好的地震波传播模拟效果。

3.2 三层HI介质模型

为进一步体现本文提出的八阶ONAD方法的优势，我们选择了三层TI介质模型的地震波模拟，并同八阶LWC方法的数值结果进行比较。三层TI介质模型中的模型参数如表1所示。实验的模拟区域为0≤(x,z)≤12 km，选择频率f0=20 Hz的Ricker震源子波函数，同样f1=f3=f(t)，震源位于模拟区域中央，取Δx=Δz=40 m为空间模拟步长，取Δt=2 ms为时间模拟步长。

图3 八阶ONAD方法在t=0.6 s时刻水平位移分量和垂直位移分量的波场瞬时快照Fig.3 Snapshots of wave fields on the horizontal and vertical components generated by the eighth-order ONAD method at time t=0.6 s

图4 八阶LWC方法在t=0.6 s时刻水平位移分量和垂直位移分量的波场瞬时快照Fig.4 Snapshots of wave fields on the horizontal and vertical components generated by the eighth-order LWC method at time t=0.6 s

表1 三层介质模型中的参数

图5给出了粗网格(Δx=Δz=40 m)情况下在t=1.2 s时刻分别由八阶ONAD方法(图5-A)和八阶LWC方法(图5-B)模拟获得在水平位移方向的波场瞬时快照。从图5可以观察出，这2种方法模拟的数值波场几乎相同，但是图5-B说明八阶LWC方法有着明显的数值频散现象。同样模拟条件下，我们发展的新方法获得了清晰的波场瞬时快照(图5-A)，没有看到数值频散。因此，在高采样点即粗网格条件下，八阶ONAD方法在抑制数值频散上明显优于传统的八阶LWC方法，能很好地应用于复杂介质模型地震波场数值模拟研究中。

4 结 论

弹性波在地球内部介质中的传播本质上是一个波动方程的正演模拟行为，为了更好地掌握地下信息，正演模拟技术成了计算数学和地球物理学中一个非常重要的研究项目。为了进一步提高正演数值方法的计算速度和模拟精度，在Yang等人[10]的研究基础上，本文发展了一种八阶ONAD数值模拟方法。基于二维弹性波传播方程，本方法以八阶精度的近似解析离散(NAD)算子[17]替代传统有限差分方法的差分算子对空间偏导数进行离散，以截断泰勒展开式对时间导数进行四阶离散。

图5 t=1.2 s时刻由八阶ONAD方法和八阶LWC方法产生的水平位移分量的波场快照Fig.5 Snapshots of wave fields for the horizontal displacement component at time t=1.2 s

以VTI介质模型为例进行了波场模拟以考察八阶ONAD方法的计算效率，结果表明八阶ONAD方法在粗网格条件下能够提供高精度的数值解，且存储量是八阶LWC方法的60.69%，计算速度约为传统八阶LWC方法的2.5倍。为进一步体现本文提出的八阶ONAD方法数值模拟的优越性，我们将该方法模拟二维弹性波在多层复杂介质中的传播。数值结果显示，八阶ONAD方法在高采样点即粗网格条件下能很好地抑制由强间断速度引起的数值频散结果，有利于强间断介质模型中使用高采样点进行地震波场模拟。总之，本文所提出的八阶ONAD方法是一种高效率、高精度的正演数值模拟方法。

[1] Chen K H. Propagating numerical model of elastic wave in anisotropic in homogeneous media-finite element method[J]. Symposium of 54th SEG, 1984, 54: 631-632.

[2] Dablain M A. The application of high-order differencing to scalar wave equation[J]. Geophysics, 1986, 51: 54-66.

[3] Chapman C H, Shearer P M. Ray tracing in azimuthally anisotropic media-Ⅱ: quasi shear wave coupling[J]. Geophys J, 1989, 96: 65-83.

[4] Komatitsch D, Vilotte J P. The spectral element method: an efficient tool to simulate the seismic responses of 2D and 3D geological structures[J]. Bull Seism Soc Am, 1998, 88: 368-392.

[5] 董良国,马在田,曹景忠,等.一阶弹性波方程交错网格高阶差分解法[J].地球物理学报,2000,43(3):411-419.

Dong L G, Ma Z T, Cao J Z,etal. A staggered-grid high-order difference method of one-order elastic wave equation[J]. Chinese J Geophys, 2000, 43(3): 411-419. (In Chinese)

[6] Yang D H, Liu E, Zhang Z J,etal. Finite-difference modeling in two-dimensional anisotropic media using a flux-corrected transport technique[J]. Geophys J Int, 2002, 148: 320-328.

[7] 杨顶辉.双相各向异性介质中弹性波方程的有限元解法及波场模拟[J].地球物理学报,2002,45(4):575-583.

Yang D H. Finite element method of the elastic wave equation and wave field simulation in two-phase anisotropoc media[J]. Chinese J Geophys, 2002, 45(4): 575-583. (In Chinese)

[8] Liu Y, Sen M K. A new time-space domain high-order finite-difference method for the acoustic wave equation[J]. J Comput Phys, 2009, 228: 8779-8806.

[9] Yang D H, Teng J W, Zhang Z J,etal. A nearly-analytic discrete method for acoustic and elastic wave equations in anisotropic media[J]. Bull Seism Soc Am, 2003, 93(2): 882-890.

[10] Yang D H, Peng J M, Lu M,etal. Optimal nearly-analytic discrete approximation to the scalar wave equation[J]. Bull Seism Soc Am, 2006, 96(3): 1114-1130.

[11] Yang D H, Song G J, Chen S,etal. An improved nearly analytical discrete method: an efficient tool to simulate the seismic response of 2-D porous structures[J]. J Geophys Eng, 2007, 4: 40-52.

[12] 王磊,杨顶辉,邓小英.非均匀介质中地震波应力场的WNAD方法及其数值模拟[J].地球物理学报,2009,52(6):1526-1535.

Wang L, Yang D H, Deng X Y. A WNAD method for seismic stress-field modeling in heterogeneous media[J]. Chinese J Geophys, 2009, 52(6): 1526-1535. (In Chinese)

[13] Yang D H, Wang N, Chen S,etal. An explicit method based on the implicit Runge-Kutta algorithm for solving the wave equations[J]. Bull Seism Soc Am, 2009, 99(6): 3340-3354.

[14] Ma X, Yang D H, Liu F Q. A nearly analytic symplectically partitioned Runge-Kutta method for 2-D seismic wave equations[J]. Geophys J Int, 2011, 187: 480-496.

[15] 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 wave equation[J]. J Comput Phys, 2011, 48: 56-70.

[16] Yang D H, Tong P, Deng X Y. A central difference method with low numerical dispersion for solving the scalar wave equation[J]. Geophysical Prospecting, 2012, 60(5): 885-905.

[17] Tong P, Yang D H, Hua B L,etal. A high-order stereo-modeling method for solving wave equations[J]. Bulletin of the Seismological Society of America, 2013, 103(2): 811-833.

[18] Lax P D, Wendroff B. Difference schemes for hyperbolic equations with high order of accuracy[J]. Communication on Pure and Applied Mathematics, 1964, 17: 381-398.

[19] Blanch J O, Robertsson A. A modified Lax-Wendroff correction for wave propagation in media described by Zener elements[J]. Geophysical Journal International, 1997, 131: 381-386.

附录A：空间高阶偏导数的逼近公式

利用局部插值近似方法，可以得到位移关于空间的二阶和三阶偏导数的逼近公式，具体推导过程见参考文献[17]。这里，只给出了关于位移的高阶偏导数的计算公式

(A1)

(A2)

(A4)

(A5)

(A6)

(A7)

其中:Δx,Δz分别代表x方向和z方向的空间步长。

A high-accuracy ONAD method and its wave-field simulation for solving elastic wave equation

ZHANG Chao-yuan

CollegeofMathematicsandComputer,DaliUniversity,Dali671003,China

Based on the elastic wave equation, this paper develops an eighth-order ONAD (optimal nearly-analytic discrete) method with high accuracy and low numerical dispersion. This new method uses the nearly analytic discrete operators with eighth-accuracy to approximate the high-order derivatives in space and the truncated Taylor series expansion with fourth-order accuracy to discretize the temporal high-order derivatives. The eighth-order ONAD method is used to model the elastic wave propagations through the VTI medium and two complex layered media. The results of the computational efficiency show that this method is obviously superior to the eighth-order Lax-Wendroff Correction (LWC) method in the computation speed and the storage capacity. The wave-fields modeling results show that the eighth-order ONAD method can effectively eliminate the numerical dispersion caused by the speed of strong discontinuity in the coarse grid, and is conducive to the wave-field simulation in the discontinuous medium using a coarse grid. Therefore, the eighth-order ONAD method has great application potentiality in seismic exploration.

elastic wave equation; ONAD method; numerical dispersion; wave-field simulation; high accuracy

10.3969/j.issn.1671-9727.2016.05.13

1671-9727(2016)05-0623-07

2014-01-09。

国家自然科学基金资助项目(41230210, 41464004); Statoil Company资助项目(4502502663); 云南省教育厅科学研究基金重点项目(2013Z152)。

张朝元(1978-),男,副教授,研究方向:地震波方程的数值方法及波场模拟, E-mail:zcy_km@163.com。

P631.4; O241

A