2016, Vol. 31
随着地震勘探目标日益复杂化和勘探要求日益精细化，全波形反演以其高精度、多参数建模的能力已成为目前研究的热点.而频率域正演是频率域全波形反演的基础(Virieux and Operto，2009;刘璐等，2013).相对时间域正演而言，频率域正演是按频率片进行求解，其计算误差分配到每一个网格点上，不存在累积误差，每个频率都是单独计算，因此特别适用于多炮并行计算.在FWI反演梯度的构造中，频率域不需要进行波场存储(Symes and Dussaud，2007; Clapp，2009)，节省大量内存存储空间.对于黏弹性介质，在进行频率域正演时，其衰减系数可以表示为频率的函数，因此频率域正演模拟比时间域正演模拟更加方便.但是巨大的稀疏矩阵的存储和大规模阻抗矩阵方程组求解是制约频率域有限差分正演发展的一大难题(任浩然等，2009).Lysmer和Drake(1972)首先提出了使用有限元法来实现频率域正演模拟，在此基础上，Marfurt和Shin对其进一步发展并用于地震波反演成像.Pratt(1990)将有限差分算法用于频率域正演中进行层析成像和地震反演的研究，并且Pratt和Worthington(1990)提出了频率-空间域五点差分格式，但其频散现象严重，在保持相速度误差为1%内，每个波长至少需要13个网格点，计算效率和正演模拟精度低，并未得到广泛的应用.为了改善这个问题，Jo等(1996)提出了优化9点有限差分格式，保持相速度误差在1%内，每个波长内所需网格点减小到4个，明显压制了频散现象，提高了计算的精度.Shin和Sohn(1998)基于优化9点差分格式的基础上，将网格点扩充到25点，引入0°，26.6°，45°和63.4°旋转坐标系来近似拉普拉斯算子，进而构造25点差分格式，此差分格式每个波长内只需要2.5个网格点将误差控制在1%内.Min等(2000)将加权平均算子的思想引入到25点频率域正演中，构造加权平均25点差分格式，并应用于弹性波频率域正演.Hustedt等(2004)结合优化九点差分算法和交错网格算法提出了混合九点算法，此算法适合于变密度有限差分正演.在国内，曹书红和陈景波(2012)在常规四阶九点差分格式的基础上，引入了45°旋转坐标系推导出优化17点差分格式，该方法每个波长内所需2.56个网格点即可控制误差在1%内.无论25点差分格式还是17点差分格式，他们的系数矩阵的带宽相对九点差分格式而言，增加了一倍，在运用直接法求解大型稀疏矩阵方程组时，明显增加了LU分解的时间(Hustedt，2004)，计算效率下降，不利于应用在大规模的全波形反演中.针对稀疏矩阵的带宽以及大型阻抗矩阵求解问题，刘璐等(2013)提出了优化15点差分格式，此方法在减小系数矩阵带宽的同时，把每个波长内所需网格点数减少到2.97个，进而提高了计算效率.
在实际资料处理中，地质模型比较复杂，往往出现纵横采样间隔不相等的情况，在这种情况下，上述所提出的旋转坐标系差分方法不再适用.因此，针对上述问题，Chen(2012)利用平均导数法(ADM)提出了基于平均导数的二阶九点差分格式，将每个波长内所需网格点降低到4个，但是其精度相对较低.因此在基于平均导数法的基础上，张衡等(2014)提出了一种显式的四阶25点差分格式(简称ADM25点法)，将每个波长所需网格点降低到2.78个，唐祥德等(2015)提出显式的四阶精度17点差分格式，通过使用模拟退火法进行差分系数优化，使得每个波长仅需2.4网格点，这两种方法都进一步提高了正演的精度.但是无论ADM九点差分格式还是ADM17点差分格式亦或是ADM25点差分格式都不能同时兼顾效率和精度.
针对上述提到的效率与精度兼顾的问题，本文将平均导数法与15点差分格式相结合，推导出一种适用于纵横向采样间隔不相等情况下的15点差分格式，其基本思想：在纵向上采用二阶差分格式，横向上采用显式的四阶差分格式，此差分格式使得非零元素尽可能的分布在主对角线两侧，相比于ADM25点与ADM17点差分格式，大大减小了系数矩阵的带宽，节约了正演求解的时间，提高了计算效率，并且横向采用显式的四阶差分有效的压制了由角度各向异性引起的数值频散.本文通过最小二乘法求得优化系数，经频散分析，每个波长所需网格点数减小到2.85个，相比于ADM9点算法以及常规四阶九点算法，明显提高了正演的精度.为了吸收衰减人为边界反射，本文引入最佳匹配层边界条件(PML)进行吸收衰减，确保了正演算法的正确性.最后进行了数值模拟，通过正反演模型测试表明了ADM15点算法能适用于不同纵横向采样间隔的正演模拟，而且横向上采用显式四阶差分对于横向速度变化剧烈的地区反演效果较好.
1 基于平均导数法的15点差分格式构造二维笛卡尔坐标系中，均匀介质的频率-空间域的标量波波动方程为
|
(1) |
式中：P是介质中每一点的波场值，w为角频率，v为介质速度，S(w)为震源项.xs和zs是震源的坐标.
在实际应用中，经常会出现水平方向采样间隔不等于深度方向采样间隔的情况，此时只适用于等间隔采样的旋转坐标系下的混合网格就不再适用.因此本文基于Chen等(2012)平均导数的思想提出了基于平均导数的15点差分格式(如图 1).
|
图 1 本文差分格式及系数加权示意图 Figure 1 Finite difference computational grid and coefficients of 15-point scheme |
基于平均导数的15点差分格式是由拉普拉斯算子、质量加速度项以及优化差分系数构成.其构造过程如下，空间导数项：在x方向的空间导数的四阶差分近似的波场值用其正交方向上3个网格点的波场值进行加权平均，在z方向的空间导数的二阶差分近似的波场值用其正交方向上5个网格点的波场值进行加权平均.而质量加速度项用其本身与周围的14个点进行加权平均求取.由此可得平均导数的15点有限差分格式为
|
(2) |
式中Δx、Δz分别为原坐标系中水平方向与深度方向的空间采样间隔.并且：
|
|
其中∂1，∂2，β1，β2，β3，b1，b2，b3，b4，b5，b6为加权优化系数，满足：
|
将平面波P(x，z，w)=P0e-i(kxx+kzz)带入(2)式的ADM15点差分格式，可以求得相速度的频散关系式为
|
(3) |
vph表示相速度，v表示实际速度. 其中A=cos(kxΔx)，B=cos(2kxΔx)，C=cos(kzΔz)，
当Δx≥Δz;即
|
(4) |
本文利用最小二乘法，将频散曲线与1作差，进而求得ADM15点差分格式的最优化系数.在求取系数时，1/G的取值范围为[0.001，0.4]，间隔为0.0001，由于差分格式的不对称性，传播角度的取值范围为[0，π/2]间隔为1°，最终采用最小二乘法求取不同r=ΔxΔz值下的优化系数如表 1所示.
|
表 1
当Δx≥Δz时，不同的  |
本文通过常规九点四阶差分算法与ADM15点差分算法的频散曲线图(如图 2所示)对比可知，以相速度误差控制在1%内为标准，常规9点四阶差分格式每个波长内最少需要5个网格点，才能满足精度要求.而本文的ADM15点差分格式仅需要2.85个网格点，就能达到相同的精度，与ADM25点四阶差分格式所需2.78个网格点相近，因此ADM15点差分格式压制频散的效果明显.横向上采用四阶差分对ADM15点差分格式的角度各向异性有了较好的改善.从内存存储量来比较，如果在同一精度下，常规四阶九点差分格式需要nx×nz个网格点，ADM25点四阶差分格式则需要
|
图 2 不同 |
大型稀疏线性方程组的求解和存储是制约频率域正演的难题.稀疏线性方程组在求解时，系数矩阵的带宽是影响矩阵LU分解效率的重要因素.若模型的网格点为nz×nx，常规四阶九点差分格式非零元素的带宽为4×nx，ADM25点四阶差分格式的带宽为4×nx+5，ADM九点差分格式的带宽为2×nx+3，而本文所采用的ADM15点差分格式的带宽为2×nx+5，约为常规九点差分格式和ADM25点差分格式的一半，因此在进行LU分解时，节省了计算时间.综上所述，ADM15点差分格式在保证较高精度的前提下，降低计算成本，提高了计算效率.
由于平均导数15点差分格式不具有几何对称性，所以我们仍需讨论Δx＜Δz的情形.
当Δx≥Δz;即 
|
(5) |
采用上述同样的方法对差分系数进行优化，得到不同r值下的优化系数(如表 2所示.)
|
|
表 2
当Δx≤Δz时，不同的 |
在有限差分正演模拟中，由于计算机硬件和计算时间的限制，人们往往对无限介质加以截断，截断后在正演的边界处就会产生了毫无意义的人为边界反射，严重的干扰了正演的结果.因此必须构造边界条件来吸收衰减边界反射，确保正演模拟的准确性.本文采用由Berenger(1994)提出的完全匹配层(PML)吸收边界条件，此吸收边界被公认为目前效果最好的吸收边界条件.完全匹配层介质中的波动方程可以看作是常规波动方程的推广，在该介质中，地震波传播时振幅随指数衰减；对于弹性介质参数相同而衰减系数不同的PML，它们的波阻抗将完全匹配，因此理论上入射波将无反射的传播，反射主要由数值离散引起.本文基于二维频率-空间域标量波方程，构建了ADM15点差分格式的PML方程.
二维频率-空间域标量波PML波动方程(任浩然等，2009)表示为
|
(6) |
其中ex和ez分别是x和z方向上的衰减项，其表达式为： 
|
其中f0为震源主频，a0为经验值(吴国忱等，2007)，li为PML区域计算点与PML边界的距离，LPML为PML边界的总厚度.
将PML波动方程带入ADM15点差分方程，得ADM15点PML差分方程为
|
(7) |
其中：
|
对于均匀介质模型我们采用单道记录与解析解对比验证本文方法的精度，首先我们采用Alford(1974)提出的均匀介质模型的解析解表达式为
|
F和F-1分别为傅里叶变换与傅立叶反变换，H0(2)为第二类零阶汉克函数，f(t)为震源函数，本文采用雷克子波为震源子波.
观测系统如图 3所示，本文采用均匀介质模型速度为2000 m/s，模型网格点数为201×201，震源采用雷克子波，其主频为30 Hz，模拟频率范围1~75 Hz，震源位于模型的中心位置如下图所示(图 3)，检波点距震源位置右侧60个网格点，边界条件采用完全匹配层边界条件(PML).
|
图 3 均匀速度模型 Figure 3 Homogeneous velocity model |
为了进行精度和计算效率的对比测试，在相同的测试条件下，我们分别采用常规四阶九点差分格式、ADM 9点差分格式、ADM 25点差分格式与本文ADM15点差分格式进行频率-空间域的正演模拟.首先我们采用了纵横采样间隔比为2，横向采样间隔8 m，纵向采样间隔为4 m的模型进行数值模拟.由图 4a、b、c、d、e对比可知本文方法与ADM25点法数值模拟所得单道记录与解析解单道记录能很好地重合，而常规四阶九点格式和ADM 9点格式所得单道记录数值频散严重，与解析解误差较大.通过以上对比可以看出本文方法在精度上的优势.
|
图 4 5种差分格式的单道记录 (a)解析解的单道记录；(b)常规四阶九点格式单道记录;(c)ADM9点格式单道记录; (d)ADM15点单道记录;(e)ADM25点单道记录. Figure 4 The result of single trace record with different finite scheme (a)The anlytical sloution;(b)The conventional forth-order 9-point scheme;(c)9-point scheme based on average-derivative method; (d)15-point scheme based on average-derivative method;(e)25-point scheme based on average-derivative method |
在相同的测试条件下，仍采用纵横空间采样间隔比为2，对上述均匀介质模型进行计算效率的测试.经测试对比可知，在保持同等计算精度下，常规四阶九点格式单炮数值模拟所用时间为2240 s，ADM9点单炮数值模拟所用时间为520 s，ADM25点格式单炮数值模拟所用时间为786 s，本文ADM15点差分格式单炮数值模拟所用时间为489 s.因此在保持相同精度下，本文方法计算效率高于ADM9点差分格式和ADM25点格式，并明显高于常规九点格式，这是由于ADM15点差分格式在精度得到较大的提高，而其系数矩阵的带宽并没有明显增加，降低了计算成本.综合计算效率与计算精度的分析，ADM15点差分格式在不明显增加计算量的前提下，较好的压制频散，保证了正演模拟的精度，不失为一种较好的数值模拟方法.
3.2 凹陷模型测试为了进一步验证本文ADM15点差分方法适用于纵横采样间隔不相等的情况，本文采用复杂的速度模型来进行数值模拟测试，我们采用横向网格间距为12.5m，纵向网格间距为4 m.横纵向采样间隔比为3.125.速度模型网格点数为301×301，数值模拟时，雷克子波主频取20 Hz，震源位置如图所示，PML层的厚度为40，传播时间为2.78 s，时间采样间隔为0.0054 s.
|
图 5 洼陷速度模型 Figure 5 The hollow velocity model |
本文分别采用适用于不等间隔采样的常规四阶9点差分格式、ADM9点差分格式、与本文推导的ADM15点差分格式进行频率-空间域数值模拟.为了验证其正演模拟的正确性，本文采用了时间域12阶有限差分数值模拟进行对比(图 6d所示)，通过图 6c、d地震记录的对比可以看出，本文ADM15点差分格式与时间域12阶有限差分所得正演记录的同相轴位置有很好的一致性，并且数值模拟精度也相当.从图 6a、b、c的地震记录中所标识的区域，可以明显的看出常规九点格式数值模拟地震记录(图 6a)与ADM九点格式数值模拟地震记录(图 6b)数值频散严重，而本文ADM15点格式波场特征保持良好，频散压制效果显著.本文引入的PML边界吸收条件很好的压制了人工边界反射(图 6a、b、c所示)，确保了频率-空间域有限差分数值模拟的正确性.
|
图 6 洼陷模型正演记录 (a)常规九点格式得到的地震记录；(b)ADM九点格式得到的地震记录；(c)ADM15点格式得到的地震记录；(d)时间域12阶差分精度得到的地震记录. Figure 6 Forward modeling record of hollow model (a)Seismographic record computed with conventional nine-point scheme;(b)Seismographic record computed with ADM nine-point scheme;(c)Seismographic record computed withADM 15-point scheme;(d)Seismographic record computed with time domain 12-order finite difference scheme. |
为了进一步验证本文方法的实用性，对其进行了Marmousi模型全波形反演测试.由于Marmousi较大，对其选取了一部分进行了反演.初始模型为其真实模型横向与纵向大尺度平滑后的速度.本文参考Sirgue和Pratt(2004)等人对频率选取策略，分别选取以下频率：1.95 Hz、3.26 Hz、4.56 Hz、6.56 Hz、9.11 Hz、12.36 Hz、15.63 Hz、18.88 Hz、22.78 Hz、26.69 Hz、30.60 Hz进行了反演.震源是主频为20 Hz的雷克子波，反演策略为由低频到高频的多尺度反演.
观测系统设为地表水平方向(30 m，20 m)处开始放炮，炮间距120 m，炮数为60，道间距30 m，道数240，正演模拟采用的横纵向采样间隔分别为30 m，20 m.横纵项间隔比为1.5，反演方式：前一次反演的最终速度为下一次反演的初始速度.最终反演结果如图 7所示.
|
图 7 速度(a)真实模型;(b)初始模型;(c)22.78 Hz反演结果;(d)30.6 Hz点反演结果 Figure 7 Velocity(a)true model;(b)initial model;(c)22.78 Hz inversion result;(d)30.06 Hz inversion result |
所得最终反演速度场与真实速度进行对比，如图 8a、b所示，尽管初始速度与真实速度差距较大，但是全波形反演测试所得的反演速度场与实际速度场的匹配效果良好，由于在波形反演中，梯度算子对深层速度不敏感，因而造成图 8b的反演速度不如图 8a中的反演速度效果好.对于横向速度变化较剧烈的区域，由于ADM15点差分格式横向上采用四阶差分，反演速度场的匹配程度较高.
|
图 8 反演结果与真实速度对比 (a)(z=700 m);(b)(z=980 m). Figure 8 Comparison of inverted velocity and true velocity |
本文将平均导数法与15点差分格式相结合，推导出一种适用于纵横方向采样间隔不相等情况下的频率-空间域正演模拟算法，拓宽了原有方法的适用范围.相比于常规四阶九点差分格式每个波长内所需5个网格点与ADM9点差分格式每个波长所需4个网格点而言，本文通过频散分析优化加权系数，使得每个波长内仅需2.85个网格点就能使误差控制在1%内，明显的提高了数值模拟精度.本文采用的ADM15点差分格式将系数矩阵的带宽尽可能的聚焦在主对角线的两侧，减少了LU分解所用时间，因此在不明显提高计算量的情况下提高了计算效率.最后进行了频率域全波形正反演测试，进一步验证了本文方法正演模拟的高效性，并且横向采用显式四阶差分对横向速度变化剧烈的地区反演效果好，具有较好的实用性.
致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！| [] | Alford R M, KellyK R, Boore D M .1974. Accuracy of finite-difference modeling of the acoustic wave equation[J]. Geophysics, 39 (6) : 834–842. DOI:10.1190/1.1440470 |
| [] | Berenger J P .1994. A perfectly matched layer for the absorption of electromagnetic waves[J]. Journal of Computational Physics, 114 (2) : 185–200. DOI:10.1006/jcph.1994.1159 |
| [] | Cao S H, Chen J B .2012. A 17-point scheme and its numerical implementation for high-accuracy modeling of frequency-domain acoustic equation[J]. Chinese J. Geophys.(in Chinese), 55 (10) : 3440–3449. DOI:10.6038/j.issn.0001-5733.2012.10.027 |
| [] | Chen J B .2012. An average-derivative optimal scheme for frequency-domain scalar wave equation[J]. Geophysics, 77 (6) : T201–T210. DOI:10.1190/geo2011-0389.1 |
| [] | Hustedt B, Operto S, Virieux J .2004. Mixed-grid and staggered-grid finite-difference methods for frequency-domain acoustic wave modelling[J]. Geophys. J.Int., 157 (3) : 1269–1296. DOI:10.1111/gji.2004.157.issue-3 |
| [] | Jo C H, Shin C, Suh J H .1996. An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolator[J]. Geophysics, 61 (2) : 529–537. DOI:10.1190/1.1443979 |
| [] | Liu L, Liu H, Liu H W .2013. Optimal 15-point finite difference forward modeling in frequency-space domain[J]. Chinese J. Geophys. (in Chinese), 56 (2) : 644–652. DOI:10.6038/cjg20130228 |
| [] | Lysmer J, Drake LA. 1972. Afiniteelementmethod for seismology[A].//Bolt B A. Methods inComputational Physics.Seismology:Surface Waves and Earth Oscillations[M]. NewYork:Academic Press, 11:181-216. |
| [] | MarfurtK J .1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic waveequations[J]. Geophysics, 49 (5) : 533–549. DOI:10.1190/1.1441689 |
| [] | Min D J, Shin C, Kwon B D, et al .2000. Improved frequency-domain elastic wave modeling using weighted-averaging difference operators[J]. Geophysics, 65 (3) : 884–895. DOI:10.1190/1.1444785 |
| [] | Pratt R G .1990. Frequency-domain elastic wave modeling by finite differences:A tool for crosshole seismic imaging[J]. Geophysics, 55 (5) : 626–632. DOI:10.1190/1.1442874 |
| [] | Pratt R G .1999. Seismic waveform inversion in the frequency domain, Part 1:Theory and verification in a physical scale model[J]. Geophysics, 64 (3) : 888–901. DOI:10.1190/1.1444597 |
| [] | Ren H R, Wang H Z, Gong T .2009. Seismic modeling of scalar seismic wave propagation with finite-difference scheme in frequency-space domain[J]. Geophysical Prospecting for Petroleum (inChinese), 48 (1) : 20–26. |
| [] | Shin C, Sohn H .1998. A frequency-space 2-D scalar wave extrapolator using extended 25-point finite-difference operator[J]. Geophysics, 63 (1) : 289–296. DOI:10.1190/1.1444323 |
| [] | Sirgue L, Pratt R G .2004. Efficient waveform inversion and imaging:A strategy for selecting temporal frequencies[J]. Geophysics, 69 (1) : 231–248. DOI:10.1190/1.1649391 |
| [] | Symes W W, Dussaud E. 2007. Optimal scaling for prestack depth migration[C].//SEG Technical Program Expanded Abstracts,1770-1774,doi:10.1190/1.2792835. |
| [] | Tang X D, Liu H, ZhangH .2015. Frequency-space domain acoustic modeling based on an average-derivative and GPU implementation[J]. Chinese J. Geophys. (in Chinese), 58 (4) : 1341–1354. DOI:10.6038/cjg20150421 |
| [] | Virieux J, Operto S .2009. An overview of full-waveform inversion in exploration geophysics[J]. Geophysics, 74 (6) . |
| [] | Yin W, Yin X Y, Wu G C, et al .2006. The method of finite difference of high precision elastic wave equations in the frequency domain and wave-field simulation[J]. Chinese J.Geophys.(in Chinese), 49 (2) : 561–568. DOI:10.3321/j.issn:0001-5733.2006.02.032 |
| [] | Zhang H, Liu H, Liu L, et al .2014. Frequency domain acoustic equation high-order modeling based on an average-derivative method[J]. Chinese J.Geophys.(in Chinese), 57 (5) : 1599–1611. DOI:10.6038/cjg20140523 |
| [] | 曹书红, 陈景波.2012. 声波方程频率域高精度正演的17点格式及数值实现[J]. 地球物理学报, 55 (10) : 3440–3449. DOI:10.6038/j.issn.0001-5733.2012.10.027 |
| [] | 刘璐, 刘洪, 刘红伟.2013. 优化15点频率-空间域有限差分正演模拟[J]. 地球物理学报, 56 (2) : 644–652. DOI:10.6038/cjg20130228 |
| [] | 任浩然, 王华忠, 龚婷.2009. 标量地震波频率-空间域有限差分法数值模拟[J]. 石油物探, 48 (1) : 20–26. |
| [] | 唐祥德, 刘洪, 张衡.2015. 基于加权平均导数的频率-空间域正演模拟及GPU实现[J]. 地球物理学报 (4) : 1341–1354. DOI:10.6038/cjg20150421 |
| [] | 殷文, 印兴耀, 吴国忱, 等.2006. 高精度频率域弹性波方程有限差分方法及波场模拟[J]. 地球物理学报, 49 (2) : 561–568. DOI:10.3321/j.issn:0001-5733.2006.02.032 |
| [] | 张衡, 刘洪, 刘璐, 等.2014. 基于平均导数方法的声波方程频率域高阶正演[J]. 地球物理学报, 57 (5) : 1599–1611. DOI:10.6038/cjg20140523 |