2013, Vol. 56
2. 中国科学院大学, 北京 100049
2. University of Chinese Academy of Sciences, Beijing 100049, China
声波方程有限差分方法因为计算效率高、所需内存小、实现简单，而广泛应用于地震正演研究[1-9]，同时是地震逆时偏移成像技术得以快速发展的基础[10-13].网格频散是有限差分中至关重要的问题[1, 14]，直接影响着有限差分法在波动方程中的应用.网格频散是由对时间和空间偏导数的网格化造成的，相速度变成了网格间距的函数[1, 15]，这会导致地震波的数值相速度不等于地球介质的真实相速度，使得波场模拟的精度降低.一般来说，如果存在时间频散，则高频的相速度的增大；如果存在空间频散，则高频的相速度减小[15].
为了缓解或者压制网格数值频散效应，一般采用两种措施：其一是采用低阶差分格式，要求时间步长和空间步长非常小，这会造成所需内存和计算量过大而无法应用于三维的情况，特别是难以应用于实际地震观测系统的数据逆时偏移（空间采样间隔已固定）；其二是采用高阶差分格式，一般情况下，时间方向采用二阶差分格式，空间方向采用高阶差分格式或伪谱法.高阶差分格式主要是利用等波纹准则或最小误差准则优化设计差分系数实现对空间偏导数的近似.伪谱法则是通过正、反傅里叶变换来实现空间偏导数的精确求解.两种方法相比，高阶空间差分格式具有精度适中、计算效率高的特征而得到广泛应用，特别是用于地震逆时偏移成像.伪谱法计算精度高，但因计算量过大，一般常用于地震正演模拟.此外，Dablain[15]和Chen[16, 17]提出了时间方向导数的四阶有限差分格式，提高了波场模拟的精度.
在空间方向偏导数差分格式近似计算方面已开展了大量研究工作，积累了大量的研究成果.考虑到在地震波场数值模拟中，需在时间方向和空间方向同时网格化，两者共同决定着数值频散特征，进而共同影响着波场模拟的精度，使得单纯依靠空间方向偏导数计算精度的提高，并不意味着可直接导致波场模拟精度的提高.鉴于上述原因，近年来Etgen[18]提出了相速度误差最小的高斯-牛顿迭代法确定时间-空间域有限差分系数.Finkelstein提出了在时间-空间域内设计与确定有限差分格式的新方法，以期提高正演模拟计算精度.其方法原理是时间-空间域对波数进行泰勒展开、特定频率点满足相速度关系、群速度关系的组合，用以控制精度的阶数和频散[19, 20].Liu提出了在时间-空间域对波数进行泰勒展开，通过匹配系数确定有限差分系数，提高了时间-空间域有限差分算子精度的阶数.频散关系分析与正演模拟表明与传统的差分方法相比，时间-空间域对波数进行泰勒展开确定有限差分系数能够提高波场模拟精度[21-23].
从时间-空间域差分方法的实现中，可以看出有以下两点能够进一步提高，其一是波数方向上的泰勒展开近似限制了差分格式只能适于一定的波数范围，其二是利用波场特定传播方向上的频散关系来确定差分格式的系数，这不仅难以保证其他传播方向的精度，而且还存在着会导致数值各向异性出现的可能.为此，本文引入了特定波数点满足时间-空间域频散关系的方法[19, 20].通过震源频率、速度和空间间距选取合适的波数点满足频散关系来确定有限差分系数，方法本身无需在波数方向上进行泰勒展开近似，同时能够考虑波场在不同传播方向上的频散关系，这会使得新的差分格式不仅能在更大的频率范围符合波场传播规律，也会使数值各向异性问题得到处理.通过频散分析和波场正演模拟，验证了本文方法的有效性.
2 声波方程时间-空间域有限差分格式对于三维声波方程
|
(1) |
在（x，y，z）三个方向上的二阶空间偏导数，使用相同有限差分格式加以计算[20-22]，公式为
|
(2) |
其中
|
(3) |
h表示网格长度，τ表示时间步长，M表示差分格式的算子长度.
对于时间导数，采用二阶有限差分近似，得到
|
(4) |
将公式（2）和公式（4）带入公式（1），可得
|
(5) |
对上式进行时间-空间傅里叶变换，可得，
|
(6) |
令
|
(7) |
式中，k表示波数，θ表示平面波与水平面夹角，Φ表示其方位角.把kx，ky，kz代入公式（6），并定义r=vτ/h，得到
|
(8) |
将
|
(9) |
比较k2j项的系数，得到
|
(10) |
上式是待定系数（Cm，m=1，2，…，M）的线性方程组，线性方程组的系数为θ和Φ的函数，在给定方向（θ和Φ）后，即可求解出差分格式系数Cm，m=1，2，…，M.
从上述时间-空间域有限差分格式实现过程中，可以看出差分格式的精度直接取确于三角函数在kh等于0时泰勒展开的近似精度，从图 1中，可以看出在波数较小范围内，泰勒展开具有高的近似精度，随着波数增大，近似精度快速下降.此外，从图 2中可以看出该方法明显优于传统的差分格式，也确实存在数值各向异性效应（即频散关系随传播方向而变化）.
|
图 1 余弦函数不同阶数泰勒展开的误差 Fig. 1 The error of Taylor expansion for the cosine function |
|
图 2 空间域泰勒展开法和时间-空间域常规方法的频散误差图(a)空间域泰勒展开法; (b)时间-空间域常规方法.其中:M=10, 速度v=3000m/s, 空间步长h=10m, 时间步长τ=0.001s Fig. 2 Dispersion errors:(a) the conventional method; (b) the time-space domain Taylor expansion method when M=10, v=3000m/s, h=10m, τ=0.001s |
本文时间-空间域有限差分格式基于下面的认识：根据震源频率、速度和网格间距确定需要满足频散关系的波数范围占总波数的比例，即
|
(11) |
其中，f是震源频率的上限，v是波速，h是空间采样间距.假设震源最高频率40Hz，波速2000m/s，空间间距20m，则需要满足频散关系的波数范围占总波数范围的80%；其它条件不变，如果波速2500m/s，则需要满足频散关系的波数范围占总波数的64%.因此速度越高，满足频散关系的波数范围越小，需要的算子长度越短.在确定的波数范围内选取M+1个均匀分布的波数点，在这些特定的波数点满足时间-空间域频散关系[18-19]，并考虑不同的传播角度.定义a kl，mh=cos（mklh），kl（l=x，y，z）的第i个分量表示为ki，l，kx=kcosθcosΦ，ky=kcosθsinΦ，kz=ksinθ，则由公式（8）得到
|
(12) |
公式（12）是一个M+1维的线性方程组，易通过迭代法或直接的矩阵分解方法解得有限差分系数.但直接法的计算量大，因此在本文中，我们采用迭代求解方法.为此，定义（12）左边的系数矩阵为A，系数向量为c，右端项为d，则（12）的解可以通过求极小化问题
|
(13) |
来达到.容易算得函数ψ的梯度为
|
(14) |
于是我们构造迭代算法如下[24]：
|
(15) |
其中gk=AT（Ack-d），γ1和γ2为0到1之间的数.这种解法的效率高，在更大的频率范围和更广的传播角度内频散小.
4 频散关系分析 4.1 一维有限差分系数的频散分析一维声波方程有限差分的频散公式为[22]
|
(16) |
定义1-δ是频散误差，则1-δ的值越接近0，频散误差越小；否则频散误差越大.
由图 3可以看到，随着有限差分算子长度的增加，时间-空间域常规方法[22]和新方法的频散误差都减小.时间-空间域新方法的频散误差在整个kh范围内分布比较均匀，在更广的kh范围内保持较小误差，且可以根据震源频率、空间步长调整需要优化的波数的上限使得kh比较小的情况下频散误差也很小.
|
图 3 两种不同方法随算子长度(2M+1)增加的频散误差图, 其中, v=3000m/s, 时间步长Δx=10m, 空间步长τ=0.001s.(a)时间-空间域常规方法; (b)时间-空间域新方法. Fig. 3 1D dispersion error curves for different spatial sampling points 2M+1 (a) The conventional time-space domain method; (b) The new time-space domain optimizing method when v=3000 m/s, Δx=10 m and τ=0.001s. |
由图 4和5可以看出，随速度和时间步长的变化，时间-空间域常规方法和时间-空间域新方法的频散误差变化都比较小，因此这两种方法都比传统方法的稳定性好.
|
图 4 两种方法随速度变化的频散误差图, 速度分别是2500, 3500, 4500, 5500m/s, M=10, h=10m, τ=0.001s (a)时间-空间域常规方法; (b)时间-空间域新方法. Fig. 4 1D dispersion error curves for different velocities (a) The conventional time-space domain method; (b) The new time-space domain optimizing method whenM=10, h=10 m, τ=0.001s. |
|
图 5 两种方法随时间步长的变化的频散误差图 (a)时间-空间域常规方法; (b)时间-空间域新方法.时间步长分别是τ=0.0005, 0.001, 0.002, 0.003s, v=3000m/s, h=10m, M=10 Fig. 5 1D dispersion error curves for different time steps (a) The conventional time-space domain method; (b) The new time-space domain method where v=3000m/s, h=10m, M=10. |
4.2二维有限差分系数的频散分析
二维声波有限差分的频散关系为[22]
|
(17) |
其中，1-δ是频散误差，其值越接近0，频散误差越小；否则频散误差越大.
图 6和图 7分别显示了随有限差分算子长度和角度的变化，时间-空间域常规方法和时间-空间域新方法的频散误差.从图 6观察到，时间-空间域常规方法在kh比较低的时候，频散误差比较小；在kh比较大的时候，频散误差大，各向异性较强.从图 7观察到，时间-空间域新方法的频散误差在整个kh范围内分布比较均匀，在kh较大的时候，频散误差也较小，各向异性也较弱.
|
图 6 时间-空间域常规方法频散误差图 (a)M=2;(b)M=6;(c)M=10.参数是v=3000m/s, h=10m, τ=0.001s, θ=iπ/16(i=0, 1, …, 4) Fig. 6 2D dispersion error curves of the conventional time-space domain method for different spatial sampling points and angles of wave propagations where M=2, 6, 10, v=3000m/s, h=10m, τ=0.001s, θ=iπ/16(i=0, 1, …, 4) |
|
图 7 时间-空间域新方法频散误差图 (a)M=2;(b)M=6;(c)M=10.参数是v=3000m/s, h=10m, τ=0.001s, θ=iπ/16(i=0, …, 4) Fig. 7 2D dispersion error curves of the new time-space domain method for different spatial sampling points and angles of wave propagations where M=2, 6, 10, v=3000m/s, h=10m, τ=0.001s, θ=iπ/16(i=0, …, 4) |
二维声波正演使用公式（18）的震源，其中f0=150，主频约50Hz，震源的频谱如图 8所示.
|
(18) |
|
图 8 二维正演震源的频谱 Fig. 8 Spectrum of the wavele |
其中c为振幅，f0为震源频率.
图 9和图 10是时间-空间域常规方法系数和时间-空间域新方法系数的正演快照图.其中，M=10，介质速度是3000m/s，时间步长0.001s，空间步长Δx=Δz=10m，网格大小是300×300，震源在中心位置.图 9显示了时间-空间域常规方法得到的系数在500ms，900ms，1400ms时的波场快照图，从中可以看到明显的空间频散.图 10显示了时间-空间域新方法得到的系数在500ms，900ms，1400ms的波场快照图，在其他条件相同的情况下，空间频散明显减小.
|
图 9 时间-空间域常规方法正演的快照图 (a)500ms时的快照; (b)900ms时的快照; (c)1400ms时的快照.参数是速度等于3000m/s, Δx=Δz=10m, τ=1ms Fig. 9 Snapshots of wave fields using traditional time-space domain method v=3000m/s, Δx=Δz=10m, τ=1ms |
|
图 10 时间-空间域优化方法正演的快照图 (a)500ms时的快照; (b)900ms时的快照; (c)1400ms时的快照.参数是速度等于3000m/s, Δx=Δz=10m, τ=1ms. Fig. 10 Snapshots of wave fields using the proposed time-space domain method v=3000m/s, Δx=Δz=10m, τ=1ms |
图 11显示了美国勘探地球物理学会的盐丘模型，其速度从1500m/s变化到4800m/s.震源函数使用公式（18），其中震源频率f0等于45Hz.
|
图 11 BP盐丘速度模型 Fig. 11 BP saltmodel |
图 12显示了空间域泰勒展开法、时间-空间域常规方法、时间-空间域新方法的地震记录.其中，M=7，空间步长20m，时间步长1ms，对于每个速度确定有限差分系数.图 12a是空间域泰勒展开法差分格式的地震记录，频散明显；图 12b是时间-空间域常规方法的地震记录，频散也很明显；图 12c是时间-空间域新方法的地震记录，频散得到压制. Etgen[18]的方法也能达到图 12c相似的结果，但是Etgen方法使用高斯-牛顿迭代方法确定有限差分系数，对于大的速度模型的每个速度使用其方法，仅系数就需要较长的计算时间.
|
图 12 盐丘模型地震正演模型图 (a)空间域泰勒展开的地震记录; (b)时间-空间域传统方法的地震记录; (c)时间-空间域新方法的地震记录. Fig. 12 Seismograms for the BP salt model (a) the Taylor expansion method in the space domain; (b) the Taylor expansion method in the time-space domain method; (c) the new method. |
本文在时间-空间域常规方法确定声波方程有限差分系数基础上，引入了在时间-空间域内特定波数点满足频散关系为准则确定声波方程有限差分系数的方法.方法本身无需在波数方向上进行泰勒展开近似，而是通过震源频率、波速和网格间距确定波数上限，同时考虑了波场在不同传播方向上的频散关系，这会使得新的差分格式不仅能在更大范围符合波场传播规律，也会使数值各向异性问题得到处理.
通过分析时间-空间域常规方法、时间-空间域新方法在一维、二维的频散误差，以及二维的地震波数值模拟，验证了本文提出的新方法的有效性.因此，时间-空间域特定频率点满足频散关系的方法可以代替空间域泰勒展开法用于地震的正演、逆时偏移中，以提高地震波数值模拟的精度和效率.
致谢感谢审稿人对文章内容和格式提出的良好建议，感谢编辑老师对文章的排版、修正.
| [1] | de Basabe J D, Sen M K. Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equations. Geophysics , 2007, 72(6): T81-T95. DOI:10.1190/1.2785046 |
| [2] | Liu Y, Sen M K. Acoustic VTI modeling with a time-space domain dispersion-relation-based finite-difference scheme. Geophysics , 2010, 75(3): A11-A17. DOI:10.1190/1.3374477 |
| [3] | Kosloff D, Pestana R C, Tal-Ezer H. Acoustic and elastic numerical wave simulations by recursive spatial derivative operators. Geophysics , 2010, 75(6): T167-T174. DOI:10.1190/1.3485217 |
| [4] | Chu C L, Stoffa P L. Implicit finite-difference simulations of seismic wave propagation. Geophysics , 2012, 77(2): T57-T67. DOI:10.1190/geo2011-0180.1 |
| [5] | Chu C L, Stoffa P L. Determination of finite-difference weights using scaled binomial windows. Geophysics , 2012, 77(3): W17-W26. DOI:10.1190/geo2011-0336.1 |
| [6] | Alford R M, Kelly K R, Boore D M. Accuracy of finite-difference modeling of the acoustic wave equation. Geophysics , 1974, 39(6): 834-842. DOI:10.1190/1.1440470 |
| [7] | Kelly K R, Ward R W, Treitel S, et al. Synthetic seismograms: A finite-difference approach. Geophysics , 1976, 41(1): 2-27. DOI:10.1190/1.1440605 |
| [8] | 兰海强, 刘佳, 白志明. VTI介质起伏地表地震波场模拟. 地球物理学报 , 2011, 54(8): 2072–2084. Lan H Q, Liu J, Bai Z M. Wave-field simulation in VTI media with irregular free surface. Chinese J. Geophys. (in Chinese) , 2011, 54(8): 2072-2084. |
| [9] | 王珺, 杨长春, 冯英杰. 用优化通量校正传输技术压制数值模拟的频散. 勘探地球物理进展 , 2007, 30(4): 252–256. Wang J, Yang C C, Feng Y J. Elimination of numerical dispersion in finite-difference modeling by optimized flux-corrected transport. Progress in Exploration Geophysics (in Chinese) , 2007, 30(4): 252-256. |
| [10] | 刘红伟, 刘洪, 李博, 等. 起伏地表叠前逆时偏移理论及GPU加速技术. 地球物理学报 , 2011, 54(7): 1883–1892. Liu H W, Liu H, Li B, et al. Pre-stack reverse time migration for rugged topography and GPU acceleration technology. Chinese J. Geophys. (in Chinese) , 2011, 54(7): 1883-1892. |
| [11] | 李博, 李敏, 刘红伟, 等. TTI介质有限差分逆时偏移的稳定性探讨. 地球物理学报 , 2012, 55(4): 1366–1375. Li B, Li M, Liu H W, et al. Stability of reverse time migration in TTI media. Chinese J. Geophys. (in Chinese) , 2012, 55(4): 1366-1375. |
| [12] | 刘红伟, 刘洪, 邹振, 等. 地震叠前逆时偏移中的去噪与存储. 地球物理学报 , 2010, 53(9): 2171–2180. Liu H W, Liu H, Zou Z. The problems of denoise and storage in seismic reverse time migration. Chinese J. Geophys. (in Chinese) , 2010, 53(9): 2171-2180. |
| [13] | 韩令贺, 何兵寿, 张会星. VTI介质中准P波方程叠前逆时深度偏移. 地震学报 , 2012, 33(2): 209–218. Han L H, He B S, Zhang H X. Prestack reverse-time depth migration of quasi P-wave equations in VTI media. Acta Seismologica Sinica (in Chinese) , 2012, 33(2): 209-218. |
| [14] | 冯英杰, 杨长春, 吴萍. 地震波有限差分模拟综述. 地球物理学进展 , 2007, 22(2): 487–491. Feng Y J, Yang C C, Wu P. The review of the finite-difference elastic wave motion modeling. Progress in Geophysics (in Chinese) , 2007, 22(2): 487-491. |
| [15] | Dablain M A. The application of high-order differencing to the scalar wave equation. Geophysics , 1986, 51(1): 54-66. DOI:10.1190/1.1442040 |
| [16] | Chen J B. High-order time discretizations in seismic modeling. Geophysics , 2007, 72(5): SM115-SM122. DOI:10.1190/1.2750424 |
| [17] | Chen J B. A stability formula for Lax-Wendroff methods with fourth-order in time and general-order in space for the scalar wave equation. Geophysics , 2011, 76(2): T37-T42. DOI:10.1190/1.3554626 |
| [18] | Etgen J T. A tutorial on optimizing time domain finite-difference scheme: "Beyond Holberg". Stanford Exploration Project , 2007. |
| [19] | Finkelstein B, Kastner R. A comprehensive new methodology for formulating FDTD schemes with controlled order of accuracy and dispersion. IEEE Transactions on Antennas and Propagation , 2008, 56(11): 3516-3525. DOI:10.1109/TAP.2008.2005458 |
| [20] | Finkelstein B, Kastner R. Finite difference time domain dispersion reduction schemes. Journal of Computational Physics , 2006, 221(1): 422-438. |
| [21] | Liu Y, Sen M K. Advanced finite-difference methods for seismic modeling. Geohorizons , 2009, 14(2): 5-16. |
| [22] | Liu Y, Sen M K. Finite-difference modeling with adaptive variable-length spatial operators. Geophysics , 2011, 76(4): T79-T89. DOI:10.1190/1.3587223 |
| [23] | Liu Y, Sen M K. A new time-space domain high-order finite-difference method for the acoustic wave equation. Journal of Computational Physics , 2009, 228(23): 8779-8806. DOI:10.1016/j.jcp.2009.08.027 |
| [24] | 王彦飞. 反演问题的计算方法及其应用. 北京: 高等教育出版社, 2007 . Wang Y F. Computational Methods for Inverse Problems and Their Applications (in Chinese). Beijing: Higher Education Press, 2007 . |