2012, Vol. 55
2. 中国科学院研究生院, 北京 100049
2. 2. Graduate University of Chinese Academy of Sciences, Beijing 100049, China
全波形反演是近些年来地震勘探领域比较热门的研究课题，虽然目前没能大规模地应用到实际当中，但它表现的巨大潜力激起了越来越多的研究者的兴趣.全波形反演试图利用地震记录中的所有信息[1]，它是一种基于全波场模拟的数据拟合过程，可以从地震记录中定量提取信息[2]，它可在时间域或频率域实现.时间域全波形反演最先由Tarantola[3]在1984年提出，它使用显式有限差分格式；频率域全波形反演最先由Pratt和Worthington[4-6]在1990年提出，它使用隐式有限差分格式.不管是在时间域还是频率域，正演都是基础.在数值计算过程中，频率域全波形反演在多炮情况下计算更加高效，这正是由于频率域正演过程中直接法求解大型稀疏方程组时，不同炮(即方程组中的不同右端项)使用同一LU 分解因子，而时间域不同炮的波场必须分别求解.频率域全波形反演的过程中需进行多次频率域正演，如每次迭代过程中求目标函数时需一次正演求波场，在快速计算梯度时需要一次正演计算虚源的波场，求步长时需要多次正演求取不同模型上的波场[7].
频率域正演最早由Lysmer等[8]在1972 年提出，使用的是有限元方法；后Marfurt和Shin[9-11]对Lysmer等[8]的工作进一步发展，Marfurt详细比较了时间域和频率域正演使用有限元方法和有限差分方法的精度和效率[9]；Pratt等在1990年详细给出了频率域声波方程和弹性波方程差分格式的推导过程、阻抗矩阵的构造及求解方法，指出频率域正演适合多炮计算，便于模拟衰减效应(因为衰减系数可表示为频率的函数)，以及可灵活选择频率成分进行反演和层析成像等优点[4-6]；Jo等[12]在1996年提出了频率域常密度标量波动方程的最优9 点差分格式，将最小波长内所需剖分网格点数减小到3.3； Stekl等将Jo的最优9点格式推广到了变密度声波方程及黏弹介质中的波动方程[13-14]；Shin等[15]构造了频率域常密度标量波动方程的25点格式，将最小波长内所需的网格点数减少到2.5，但此格式仍是二阶精度格式；Hustedt等[16]分析了频率域变密度声波方程正演中使用的不同有限差分方法，并研究了交错网格上导数算子的四阶精度有限差分格式，优化基于直角坐标系，此四阶格式涉及13 个点，纵横轴各7个点，给正演过程中求解方程组带来过大的带宽.
为了得到优化的常密度标量波动方程四阶精度有限差分格式，以适应高精度正演模拟的需要，本文在经典的四阶精度9点格式的基础上，结合0°坐标系和45°旋转坐标系，构造了四阶精度17点格式，并在均匀模型、层状模型以及Overthrust模型上对二阶精度最优9点格式，四阶精度9 点格式和四阶精度17点格式进行了对比研究.
2 传统的最优9点格式传统的最优9点格式在频率域全波形反演的正演过程中广泛使用.它将均匀各向同性介质中频率域标量波动方程中拉普拉斯算子项Δ2P表示成0°坐标系和45°坐标系下的加权，并将质量加速度项ω2P/v2 也用Δ2P差分格式中对应点加权，用最优化方法寻找到最优加权系数，这样构造的最优9 点格式在保证所需精度的前提下大大减少了一个波长内所需的剖分网格点数，从而减少了计算量.均匀各向同性介质频率域标量波动方程为
|
(1) |
其中P为纵波波场，ω 为角频率，v为速度.
将0°坐标系和45°坐标系下的两个五点二阶中心有限差分格式结合起来应用于(1)式中的拉普拉斯算子，并将质量加速度项也表示成拉普拉斯算子中相应点的结合，则得到一个二阶精度的9 点格式差分方程:
|
(2) |
Pm，n为[xm，zn]= [x0 + (m-1)Δx，z0 + (n-1)Δz]点处的波场值，这里取Δ = Δx= Δz.
式(2)中各个系数可通过最优化方法求得，分别为a=0.5461，c=0.6248，d=0.09381.最终得到的最优9点格式可将G(G为一个波长内所需的剖分网格点数)减小到3.3.
3 四阶精度9点格式和17点格式传统的最优9点格式只具有二阶精度，不能满足高精度成像的需要.下文通过对经典的四阶精度9点格式进行优化，构造了新的四阶精度17点格式.
将四阶精度有限差分格式应用于(1)式中拉普拉斯算子，则得到四阶精度9点格式:
|
(3) |
将拉普拉斯算子表示成0°坐标系和45°坐标系下四阶差分格式(如图 1a和1b)的组合，并将质量加速度项表示成相应点的结合，则(1)式的四阶精度17点格式(如图 1c)差分方程可写成
|
(4) |
|
图 1 拉普算子的四阶差分格式 (a)0°坐标系9点四阶差分格式;(b)45°坐标系;(c)结合(a)和(b)的17点格式. Fig. 1 Fourth order finite-difference stars for Laplacian operator (a)Fourth order 9-point difference star;(b)45° rotated star;(c)17-point star combining (a) and (b). |
其中b+4c+4d+4e+4f=1.
将平面简谐波的表达式P(x，z，ω)=P0e-i(kxx+kzz) 代入式(4)，得到频散方程:
|
(5) |
其中
结合相速度定义vph =ω/k可得
|
(6) |
系数a，b，c，d，e可用最优化方法求得，目标函数定义为:
|
(7) |
其中k* =1/G，传播角度θ 的范围选为[0，π/4]是因为相速度值是关于π/4对称的；1/G的范围最大选到0.4是因为优化的四阶精度格式应该允许更大的网格间距，最终求得的优化系数分别为:a=1.0673，b= 0.8875，c= 0.0251，d= 0.0237，e=-0.0204.
图 2a、2b和2c分别给出了二阶精度最优9 点格式、四阶精度9点格式和四阶精度17点格式的频散曲线(图中0°，15°，30°，45°表示传播角度).可看出四阶9点格式由于没有优化，频散明显大于最优9点格式和17点格式，当G≥5时，相速度误差才可控制在±1%以内；而17 点格式的频散情况整体上又比最优9点格式好，当G≥2.56时，相速度误差都可控制在±1%以内.
|
图 2 三种有限差分格式的频散曲线 (a)最优9点格式；(b)四阶精度9点格式；(c)四阶精度17点格式. Fig. 2 Dispersion curves of 3 finite-difference schemes (a) Optimal 9-point scheme； (b) Fourth-order 9-point scheme； (c) 17-point scheme. |
取模型为nz×nx=100×100 的网格，波速为3000m/s，震源子波为主频30Hz的雷克子波(其频谱如图 3所示)，置于第2层的中间位置，接收点置于第24 层，模型示意图如图 4.分别考虑网格间距Δ=10m和Δ=20m，用四阶精度9点格式、二阶精度最优9点格式、四阶精度17 点格式做正演模拟(正演模拟流程如图 5)，将此三种格式在这两种网格间距情况下得到的地震记录与解析解作对比(如图 6、图 7所示)，并抽取Δ=20m 时，三种格式得到的单道记录与解析解进一步作对比(如图 8所示).
|
图 3 震源子波频谱 Fig. 3 Spectrum of source wavelet |
|
图 4 均匀模型示意图 Fig. 4 Sketch map of homogeneous model |
|
图 5 频率域正演模拟流程图 Fig. 5 Flow chart for frequency domain forward modeling |
|
图 6 均匀模型，△ = 10 m时，三种格式地震记录与解析解对比图 (a)四阶精度9点格式与解析方法；(b)最优9点格式与解析方法；(c) 17点格式与解析方法. Fig. 6 Comparison of seismograms obtained by 3 numerical schemes and analytical solution in homogeneous media when grid interval is 10 m (a) Fourth-order 9-point scheme and analytical method； (b) Optimal 9-point scheme and analytical method； (c) 17-point scheme and analytical method. |
|
图 7 均匀模型，△= 20m时，三种格式地震记录与解析解对比图 (a)四阶精度9点格式与解析方法；(b)最优9点格式与解析方法；(c) 17点格式与解析方法. Fig. 7 Comparison of seismograms obtained by 3 numerica schemes and analytical solution in homogeneous media when grid interal is 20 m (a) Fourth-order 9-point scheme and analytical method； (b) Optimal 9-point scheme and analytical method； (c) 17-point scheme and analytical method. |
|
图 8 均匀模型，△=20 m时，三种格式单道地震记录与解析解对比图 (a)四阶精度9点格式与解析方法；(b)最优9点格式与解析方法；(c)17点格式与解析方法. Fig. 8 Comparison of single traces obtained by 3 numerical schemes and analytical solution in homogeneous media when grid interval is 20 m (a) Fourth-order 9-point scheme and analytical method; (b) Optimal 9-point scheme and analytical method; (c) 17-point scheme and analytical method. |
图 6中Δ=10 m，可发现三种格式得到的地震记录与解析解地震记录能很好地重合，因为此时网格间距较小，三种格式均无明显频散.计算时对震源子波的采样频率实际为160Hz，那么可采到的频率成分为0~80Hz，从震源子波的频谱图可看出高于80Hz的频率成分振幅趋于0，所以0~80 Hz已涵盖绝大部分有效频率成分.此时最小波长内网格点数为3000/80/10，即3.75，最优9点格式和17点格式所要求的最小波长内网格点数至少为3.3和2.56(只要保证最小波长内网格点数达到要求值，所有波长成分的波频散误差都能控制在要求范围内)，所以频散误差满足要求；而四阶精度9 点格式要求最小波长内的网格点数至少为5，只有频率为0~60 Hz的波成分频散误差能满足要求，但从震源子波的频谱图可看出，高于60 Hz的频率成分振幅已经比较小，实际贡献不大，所以从地震记录图中较难辨出这小部分频率成分的频散.数值解和解析解的吻合同时也验证了解的正确性.
图 7中Δ=20 m，可看出三种格式地震记录与解析解地震记录有不同程度的出入.从Alford等[17]对频散在波形上的表现特征的描述(延迟的、变宽的、带有震荡的尾巴)可知，三种格式得到的地震记录均出现了频散.因为此时网格间距增大，每个波长内的网格点数均减少，最小波长内的网格点数减小到1.875，与三种格式的要求都有一定差距，导致了明显的频散.但从图 7的综合地震记录和图 8的单道地震记录都可看出，17点格式的频散程度明显轻于四阶精度9点格式和二阶精度最优9点格式.
4.2 层状模型为了更好地检验算法的有效性，这里再研究一个简单两层模型，网格大小，震源子波，震源和检波点位置均与简单模型中相同，速度分界面置于第50层(如图 9所示).这里令界面上层速度V1 为3000m/s，界面下层速度V2 为6000 m/s，同样分别考虑网格间距Δ=10m 和Δ=20m，将三种格式在此两种网格间距情况下得到的地震记录与解析方法得到的地震记录作对比(如图 10、图 11所示)，并进一步比较Δ=20m 时的单道记录(如图 12所示).
|
图 9 层状模型示意图 Fig. 9 Sketch map ol layered model |
|
图 10 层状模型，Δ= 10m时，三种格式与解析方法得到的地震记录对比图 (a)四阶9点格式与解析方法；(b)最优9点格式与解析方法；(c)17点格式与解析方法. Fig. 10 Comparison of seismograms obtained by 3 numerical schemes and analytical method on layered model when grid interval is 10 m (a) Fourth-order 9-point scheme and analytical method； (b) Optimal 9-point scheme and analytical method； (c) 17-point scheme and analytical method. |
|
图 11 层状模型△ = 20 m时，三种格式与解析方法得到的地震记录对比图 (a)四阶9点格式与解析方法；(b)最优9点格式与解析方法；(c)17点格式与解析方法. Fig. 11 Comparison of seismograms obtained by 3 numerical schemes and analytical method on layered model when grid interval is 20 m (a) Fourth-order 9-point scheme and analytical method； (b) Cptimal 9-point scheme and analytical method； (c) 17-point scheme and analytical method. |
|
图 12 层状模型△= 20 m时，三种格式与解析方法得到的单道地震记录对比图 (a)四阶9点格式与解析方法；(b)最优9点格式与解析方法；(c)17点格式与解析方法. Fig. 12 Comparison of single traces obtained by 3 numerical schemes and analytical method on layered model when grid interval is 20 m (a) Fourth-order 9-point scheme and analytical method； (b ) Cptimal 9-point scheme and analytical method； (c) 17-point scheme and analytical method. |
从图中可看出，与均匀模型中情况相似，网格间距较小时，三种格式均无明显频散，他们得到的地震记录与解析方法得到的地震记录基本重合；而网格间距较大时，三种格式均出现了明显频散，但从比较中同样可看出，17点格式频散情况好于四阶精度9点格式和二阶精度最优9点格式.
4.3 Overthrust速度模型最后研究一个复杂模型，这里取二维Overthrust模型的一部分，模型网格为nz×nx=150×300，所用震源子波与简单模型中相同，将其置于第10层的(3/8)×nx处，检波器点置于第2 层，如图 13.最小速度为2674m/s，最大速度为6000m/s，平均速度为4292m/s.取网格间距Δ=25m，分别用四阶精度9点格式、二阶精度最优9 点格式、四阶精度17 点格式做正演模拟，相应的综合地震记录如图 14.可以发现三幅地震记录上都存在明显的频散，但17点格式的地震记录与其他两种格式的地震记录相比相对清晰干净，频散程度比其他两者轻.
|
图 13 Overthrust 模型 Fig. 13 Overthrutt model |
|
图 14 Overthrust模型综合地震记录图(△=25m) (a)四阶精度9点格式；(b)最优9点格式；(c) 17点格式. Fig. 14 Seismograms in Overthrurt model (△=25 m) (a) Fourth-order 9-point scheme； (b) Optimal 9-point scheme； (c) 17-point scheme |
实际的地下介质均为复杂介质，17点格式在复杂Overthrust模型上的正演效果进一步证明了算法的可行性.
实现高精度正演模拟的另外一条途径是通过加密网格利用低阶格式来实现，但这种方法相比于高阶格式而言，其存储要求增加且计算效率降低[18-19].具体来讲，对本文提出的格式，如果考虑同样的网格，则三种格式所形成阻抗矩阵的阶数相同，显然，由于最优9点格式形成的阻抗矩阵带宽是其他两种格式的1/2，在对其进行LU 分解时，存储和计算效率方面都会占优势.但最优9 点格式只具有二阶精度，若要达到与17 点格式相当的精度，则网格剖分必须大大加密以减小网格间距，这样所形成的阻抗矩阵的阶数、带宽会随之大大增加，相应会带来更大的存储量，计算效率也会明显降低.
5 结 论通过本文的分析可知，要发展高精度正演，提高差分格式精度的阶数并不能保证频散情况比二阶精度差分格式要好，因此必须要使用各种方法对差分格式进行优化.本文利用Jo等构造最优9点格式的思想，将拉普拉斯算子四阶差分格式表示成两个坐标系的加权平均，并将质量加速项也表示成相应点的加权平均，求得最优化系数，构造了一个优化的四阶精度17点格式，将最小波长内的网格点数减小到2.56.并在均匀模型、层状模型和Overthrust模型上进行了频率域正演模拟，比较了传统的二阶精度最优9点格式、未优化的四阶精度9 点格式和优化了的四阶精度17点格式的正演效果，从地震记录图上均可看出17点格式的频散情况好于其他两种格式，简单模型上的数值实验明显体现17点格式在克服频散误差方面的优越性，而复杂模型上的数值实验则进一步承认了算法的可行性.
| [1] | Pratt R G. Seismic waveform inversion in the frequency domain, Part 1: Theory and verification in a physical scale model. Geophysics , 1999, 64(3): 888-901. DOI:10.1190/1.1444597 |
| [2] | Virieux J, Operto S. An overview of full-waveform inversion in exploration geophysics. Geophysics , 2009, 74(6): WCC1-WCC26. DOI:10.1190/1.3238367 |
| [3] | Tarantola A. Inversion of seismic reflection data in the acoustic approximation. Geophysics , 1984, 49(8): 1259-1266. DOI:10.1190/1.1441754 |
| [4] | Pratt R G. Frequency-domain elastic wave modeling by finite differences: A tool for cross-hole seismic imaging. Geophysics , 1990, 55(5): 626-632. DOI:10.1190/1.1442874 |
| [5] | Pratt R G, Wothington M H. Inverse theory applied to multi-source cross-hole tomography, Part 1: Acoustic wave equation method. Geophysical Prospecting , 1990, 38(3): 287-310. DOI:10.1111/gpr.1990.38.issue-3 |
| [6] | Pratt R G. Inverse theory applied to multi-source cross-hole tomography, Part 2: Elastic wave-equation method. Geophysical Prospecting , 1990, 38(3): 311-329. DOI:10.1111/gpr.1990.38.issue-3 |
| [7] | Pratt R G. Guass-Newton and full Newton methods in frequency domain seismic waveform inversion. Geophysics , 1998, 173: 922-931. |
| [8] | Lysmer Drake. Methods in computational physics, Volume 11: Seismology: Surface waves and earth oscillations: A finite-element method for seismology. MA: Academic Press Inc. , 1972: 181-216. |
| [9] | Marfurt K J. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations. Geophysics , 1984, 49(5): 533-549. DOI:10.1190/1.1441689 |
| [10] | Marfurt K J, Shin C S. The future of iterative modeling in geophysical exploration. //Eisner E, ed. Handbook of Geophysical Exploration: I Seismic Exploration, Volume 21:Supercomputers in seismic exploration. New York: Pergamon Press, 1989: 203-228. |
| [11] | Shin C S. Nonlinear Elastic wave Inversion by Blocky Parameterization. Tulsa: University of Tulsa, 1988 . |
| [12] | Jo C H, Suh J H, Shin C S. An optimal nine-point finite difference frequency-space 2-D acoustic wave extrapolator. Geophysics , 1996, 61(2): 529-537. DOI:10.1190/1.1443979 |
| [13] | tekl I. Frequency Domain Seismic Forward Modeling: A Tool for Waveform Inversion. London: Imperial College London, 1997 . |
| [14] | tekl I, Pratt R G. Accurate visco-elastic modeling by frequency-domain finite differences using rotated operators. Geophysics , 1998, 63(5): 1779-1794. DOI:10.1190/1.1444472 |
| [15] | Shin C, Sohn H. A frequency-space 2-D scalar wave extrapolator using extended 25-point finite-difference operator. Geophysics , 1998, 63(1): 289-296. DOI:10.1190/1.1444323 |
| [16] | Hustedt B, Operto S, Virieux J. Mixed-grid and staggered-grid finite-difference methods for frequency-domain acoustic wave modeling. Geohpys. J. Int. , 2004, 157(3): 1269-1296. |
| [17] | 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 |
| [18] | Chen J B. Lax-Wendroff and Nyström methods for seismic modeling. Geophysical Prospecting , 2009, 57(6): 931-941. DOI:10.1111/gpr.2009.57.issue-6 |
| [19] | Chen J B, Liu G F, Liu H. Seismic imaging based on spectral differentiation matrix and GPU implementation. Journal of Applied Geophysics , 2012, 79: 1-5. DOI:10.1016/j.jappgeo.2012.01.002 |