摘要
地震波传播正演模拟一直是地球物理与石油勘探领域的一个研究热点。本文首先提出了求解波动方程的一种新的数值方法——Runge-Kutta(龙格-库塔)方法。该方法的基本思想是:将二阶波动方程改写为一阶偏微分方程组,在空间方向上使用网格点上的位移、粒子速度及它们的梯度值的组合来近似逼近空间高阶偏导数,使一阶偏微分方程组转化为一个半离散的常微分方程组,然后对这个半离散的常微分方程组采用三阶或四阶Runge-Kutta方法进行时间推进计算,从而获得了求解波动方程的Runge-Kutta方法。本文还针对Runge-Kutta方法存在的不足进行改进,提出了一种加权Runge-Kutta方法,并将这两种方法统称为“Runge-Kutta型方法”。
本文从理论分析和数值计算两个方面对Runge-Kutta型方法进行了系统的研究,主要包括以下研究内容:分别给出了求解波动方程的一维、二维和三维Runge-Kutta方法的具体实现步骤,并使用Fourier分析方法推导了一维、二维和三维Runge-Kutta方法的稳定性条件;对一维、二维和三维Runge-Kutta方法分别进行了误差和频散分析,并与经典的Lax-Wendroff修正方法和交错网格方法进行了数值精度和频散误差的比较;对不同介质中二维和三维的声波和弹性波的传播进行了数值模拟,并给出了模拟结果;给出了加权Runge-Kutta方法的实现步骤以及取不同加权参数时其稳定性条件;并用加权Runge-Kutta方法进行了不同介质中地震波传播的数值模拟,给出了模拟结果,同时分析了计算效率和存储效率。
理论分析和数值实验结果表明,Runge-Kutta方法具有稳定性好,数值精度高,频散误差小的优点,并且在粗网格条件下能有效压制数值频散;加权Runge-Kutta方法可以进一步节省存储空间,更加有效地压制数值频散,从而可以通过使用更大的空间和时间步长以获得更快的计算速度和更小的存储量需求。因此,Runge-Kutta方法及其加权方法将在地球物理与石油勘探领域具有巨大的应用潜力。
The forward modeling research of seismic wave propagation is one of the hotspot in the fields of Geophysics and petroleum exploration. In this thesis, we present a new numerical method, which are called Runge-Kutta method, to solve wave equations. We transform the second-order wave equations into a system of first-order partial differential equations and use not only the displacement and the particle-velocity of the grids, but also their gradients to approximate the high-order spatial derivatives, resulting that we get semi-discrete ordinary differential equations (ODEs). We then apply the Runge-Kutta method of order three or four to solve the semi-discrete ODEs and obtain our Runge-Kutta method for solving wave equations. In order to improve the Runge-Kutta method, we propose a weighted Runge-Kutta method. And we identically call the two methods“Runge-Kutta type method”.
In this thesis, we systemically investigate the Runge-Kutta type method including theoretical analyses and numerical computations. The main research contents are listed as follows. Firstly, the construction of Runge-Kutta methods to solve 1D, 2D, and 3D wave equations are respectively described, and the stability conditions of the 1D, 2D, and 3D Runge-Kutta methods are derived by Fourier analysis method. Secondly, the numerical error and dispersion of the 1D, 2D, and 3D Runge-Kutta methods are respectively analyzed, and the numerical accuracy and dispersion error of the Runge-Kutta method are also compared with those of the classic Lax-Wendroff correction method and the staggered-grid method. Thirdly, 2D and 3D numerical simulation results of acoustic and elastic wave propagation in different media are presented. Fourthly, the implementation steps and stability conditions of the weighted Runge-Kutta method with different weights are given. Finally, numerical simulation results of the seismic wave propagation are shown by the weighted Runge-Kutta method, meanwhile, the computation and storage efficiencies are investigated.
The results of the theoretical analyses and numerical computations indicate that the Runge-Kutta method has good stability, high numerical accuracy, and small dispersion error, and can effectively suppress the numerical dispersion on coarse grids. And the weighted Runge-Kutta method can further save the storage and suppress the numerical dispersion, so the coarser spatial and temporal grids can be used to acquire faster computational speed and smaller storage. Therefore, the Runge-Kutta method and its weighted method have great potentiality of applying in Geophysics and petroleum exploration.
引文
[1]王一新,王家林,万明浩,等.石油综合地球物理方法与应用.北京:石油工业出版社, 1995.
[2]杜世通.地震波动力学.北京:石油大学出版社, 1996.
[3]牟永光,裴正林.三维复杂介质地震数值模拟.北京:石油工业出版社, 2005.
[4]佘德平.波场数值模拟技术.勘探地球物理进展, 2004, 27:16~21.
[5] Booth D C, Crampin S. The anisotropic reflectivity technique: theory. Geophys J R Astr Soc, 1983, 72:755-765.
[6] Booth D C, Crampin S. The anisotropic reflectivity technique: anomalous arrives from an anisotropic upper mantle. Geophys J R Astr Soc, 1983, 72:767-782.
[7] Fu L Y. Numerical study of generalized Lipmann-Schwinger integral equation including surface topography. Geophysics, 2003, 68:665-671.
[8] Zhou H, Chen X F. The localized boundary integral equation discrete wavenumber method for simulating P-SV wave scattering by an irregular topography. Bull Seism Soc Am, 2008, 98:265-279.
[9]冯英杰,杨长春,吴萍.地震波有限差分模拟综述.地球物理学进展. 2007, 22:487-491.
[10] Eriksson K, Johnson C. Adaptive finite element methods for parabolic problems I: A linear model problem. SIAM J Numer Anal, 1991, 28:43-77.
[11]邵秀民,蓝志凌.流体饱和多孔介质波动方程的有限元解法.地球物理学报, 2000, 43:264-278.
[12]杨顶辉.双相各向异性介质中弹性波方程的有限元解法及波场模拟.地球物理学报, 2002, 45:575-583.
[13] Ichimura T, Hori M, Kuwamoto H. Earthquake motion simulation with multiscale finite-element analysis on hybrid grid. Bull Seism Soc Am, 2007, 97:1133-1143.
[14] Huang B S. A program for two-dimensional seismic wave propagation by the pseudospectrum method. Comput Geosci, 1992, 18:289-307.
[15]刘洋,李承楚.双相各向异性介质中弹性波传播伪谱法数值模拟研究.地震学报, 2000, 22:132-138.
[16]赵志新,徐纪人,堀内茂木.错格实数傅立叶变换微分算子及其在非均匀介质波动传播研究中的应用.地球物理学报, 2003, 46:234-240.
[17] Yang D H, Liu E, Zhang Z J, Teng J. Finite-difference modelling in two-dimensional anisotropic media using a flux-corrected transport technique. Geophys J Int, 2002, 148:320-328.
[18]张剑峰.各向异性介质中的弹性波数值模拟.固体力学学报, 2002, 21:234-242.
[19]祝贺君,张伟,陈晓非.二维各向异性介质中地震波场的高阶同位网格有限差分模拟.地球物理学报, 2009, 52:1536-1546.
[20]张永刚.地震波场数值模拟方法.石油物探, 2003, 42:143-148.
[21]王书强.弹性波方程的紧致差分方法及波场模拟[硕士学位论文].北京:清华大学数学系, 2002.
[22]宋国杰.三维弹性波方程的改进近似解析离散化方法及波场模拟[硕士学位论文].北京:清华大学数学系, 2008.
[23] Alterman Z, Karal F C. Propagation of elastic waves in layered media by finite-difference methods. Bull Seism Soc Am, 1968, 58:367-398.
[24] Boore D M. Finite-difference methods for seismic wave propagation in heterogeneous materials in methods in computational physics. B A Bolt, ed, Academic Press, inc, 11, 1972.
[25] Kelly K R, Ward R W, Treitel S, Alford R M. Synthetic seismograms– a finite-difference approach. Geophysics, 1976, 41:2-27.
[26]李胜军,孙成禹,高建虎,等.地震波数值模拟中的频散压制方法分析.石油物探, 2008, 47:444-448.
[27]吴国忱,王华忠.波场模拟中的数值频散分析与校正策略.地球物理学进展, 2005, 20:58-65.
[28]董良国,李培明.地震波传播数值模拟中的频散问题.天然气工业, 2004, 24:53-56.
[29] Sei A, Symes W. Dispersion analysis of numerical wave propagation. J Sci Comput, 1994, 10:1-27.
[30] Adam Y. Highly accurate compact implicit schemes and boundary conditions. J Comput Phys, 1977, 24:10-22.
[31] Dablain M A. The application of high-order differencing to the scalar wave equation. Geophysics, 1986, 51:54-66.
[32] Peyret R, Taylor T D. Computational schemes for fluid flow. New York: Springer-Verlag, Springer series in computational physics, 1993.
[33] Igel H, Weber M. Sh-wave propagation in the whole mantle using high-order finite differences. Geophys Res Lett, 1995, 22:731-734.
[34] Wang S Q, Yang D H, Yang K D. Compact finite difference scheme for elastic equations. J. Tsinghua Univ. (Sci. & Tech.), 2002, 42:1128-1131.
[35]杨宽德,杨顶辉,王书强.基于BISQ方程的波场模拟.地球物理学报, 2002, 45:853- 861.
[36] Lele S K. Compact finite difference schemes with spectral-like resolution. J Comput Phys, 1992, 103:16-42.
[37] Igel H, Mora P, Riollet B. Anisotropic wave propagation through finite-difference grids. Geophysics, 1995, 60:1203-1216.
[38] Blanch J O, Robertsson A. A modified lax-wendroff correction for wave propagation in media described by zener elements. Geophys J Int, 1997, 131:381-386.
[39] Madariaga R. Dynamics of an expanding circular fault. Bull Seism Soc Am, 1976, 65:163-182.
[40] Virieux J. Sh-wave propagation in heterogeneous media: Velocity-stress finite-difference method. Geophysics, 1984, 49:1933-1957.
[41] Virieux J. P-sv wave propagation in heterogeneous media: Velocity-stress finite-difference method. Geophysics, 1986, 51:889-901.
[42] Levander R. Fourth-order finite-difference p-sv seismograms. Geophysics, 1988, 53:1425-1436.
[43] Faria E L, Stoffa P L. Finite-difference modeling in transversely isotropic media. Geophysics, 1994, 59:282-289.
[44] Graves R W. Simulating seismic wave propagation in 3d elastic media using staggered-grid finite differences. Bull Seism Soc Am, 1996, 86:1091-1106.
[45] Moczo P, Kristek J, Vavrycuk V, Archuleta R J, Halada L. 3d heterogeneous staggered-grid finite-difference modeling of seismic motion with volume harmonic and arithmetic averaging of elastic moduli and densities. Bull Seism Soc Am, 2002, 92:3042-3066.
[46] Liu Y, Sen M K. An implicit staggered-grid finite-difference method for seismic modeling. Geophys J Int, 2009, 179: 459-474.
[47] Yang D H, Peng M J, Lu M, Terlaky T. Optimal nearly analytic discrete approximation to the scalar wave equation. Bull Seism Soc Am, 2006, 96:1114-1130.
[48] Moczo P, Kristek J, Halada L. 3d fourth-order staggered-grid finite-difference schemes: stability and grid dispersion. Bull Seism Soc Am, 2000, 90:587-603.
[49]李荣华,冯国忱.微分方程数值解法.高等教育出版社, 1997.
[50] Jiang G S, Shu C W. Efficient implementation of weighted eno schemes. J Comput. Phys, 1996, 126:202-228.
[51] Kondoh Y. On thoughs analysis of numerical schemes for simulation using a kernel optimum nearly-analytical discrimination (KOND) method. J phys Soc Japan, 1991, 60:2851-2861.
[52] Konddoh Y, Hosaka Y, Ishii K. Kernel optimum nearly analytical discretization algorithm applied to parabolic and hyperbolic equations. Comput Math Appl, 1994, 27:59-90.
[53]杨顶辉,藤吉文,张中杰.三分量地震波场的近似解析离散模拟技术.地球物理学报, 1996, 39(增刊):283-291.
[54] Yang D H, Teng J W, Zhang Z J, Liu E. A nearly-analytic discrete method for acoustic and elastic wave equations in anisotropic media. Bull Seism Soc Am, 2003, 93:882-890.
[55] Yang D H, Peng J M, Lu M, Terlaky T. A nearly analytical discrete method for wave-field simulations in 2d porous media. Commun Comput Phys, 2006, 1(3):530-549.
[56] Yang D H, Lu M, Wu R S, Peng J M. An optimal nearly-analytic discrete method for 2d acoustic and elastic wave equations. Bull Seism Soc Am, 2004, 94:1982-1991.
[57]卢明.改进的近似解析离散化方法及弹性波波场模拟[硕士学位论文].北京:清华大学数学系, 2004.
[58] Yang D H, Song G J, Lu M. Optimally accurate nearly analytic discrete scheme for wave-field simulation in 3D anisotropic media. Bull Seism Soc Am, 2007, 97:1557~1569.
[59] Yang D H, Song G J, Chen S, Hou B Y. An improved nearly analytical discrete method: an efficient tool to simulate the seismic response of 2-D porous structures. J geophys Eng, 2007, 4:40-52.
[60] Yang D H, Liu E R, Song G J, Wang N. Elastic wave modeling method based on the displacement-velocity fields: an improving nearly-analytic discrete approximation. J Seismol, 2009, 13:209-217.
[61]王磊,杨顶辉,邓小英.非均匀介质中地震波应力场的WNAD方法及其数值模拟.地球物理学报, 2009, 52:1526-1535.
[62] Biot M A. Theory of propagation of elastic waves in a fluid-saturated porous solid, 1. low-frequency range, 2. higher frequency range. J Acoust Soc Am, 1956, 28:168-191.
[63] Dvorkin J, Nur A. Dynamic poroelasticity: A unified model with the squirt an the biot mechanisms. Geophysics, 1993, 58:524-533.
[64] Parra J O. The transversely isotropic poroelastic wave equation including the biot and the squirt mechanisms: Theory and application. Geophysics, 1997, 62:309-318.
[65] Yang D H, Zhang Z J. Effects of the biot and the squirt flow coupling interaction on anisotropic elastic waves. Chinese Science Bulletion, 2000, 45:2130-2139.
[66] Richetmyer R D, Morton K W. Difference methods for initial value problems. Interscience, New York, 2002.
[67] Aki K, Richards P G. Quantitative seismology (second edition). California: University Science Books. 2002: 63-66.
[68] Yang D H, Wang S Q, Zhang Z J, et al. N-times absorbing boundary conditions for compact finite-difference modeling of acoustic and elastic wave propagation in the 2D TI medium. Bull Seism Soc Am, 2003, 93:2389-2401.