地球物理学进展  2014, Vol. 29 Issue (1): 204-208   PDF    
基于多尺度迭代求解三维频率-空间域声波方程的全波形速度反演研究-理论模型
廖建平1,2, 刘和秀1, 郑桂娟3, 戴世鑫1, 彭叶辉1, 杨天春1, 王齐仁1    
1. 湖南科技大学土木工程学院地质系, 湘潭 411201 2. 煤炭资源与安全开采国家重点实验室(中国矿业大学), 北京 100083 3. 中国石油东方地球物理公司研究院, 涿州 072751
摘要:使用广义最小残量方法迭代求解三维频率-空间域声波方程, 反演时使用多尺度、多重网格的策略, 探讨了如何快速实现高分辨率的三维频率-空间域迭代法声波全波形速度反演.通过对理论模型进行三维频率-空间域迭代法声波全波形反演数值试验, 证实该方法的计算速度快、计算效率高, 反演所得速度的分辨率高.从而为基于多尺度迭代求解三维频率-空间域声波方程的全波形速度反演成像打下方法基础.
关键词多尺度     频率-空间域     迭代法     三维声波方程     全波形速度反演    
Research on multiscale iterative three dimensional acoustic wave full waveform velocity inversion in frequency-space domain-theoretical model
LIAO Jian-ping1,2, LIU He-xiu1, ZHENG Gui-juan3, DAI Shi-xin1, PENG Ye-hui1, YANG Tian-chun1, WANG Qi-ren1    
1. School of Civil Engineering of Hunan University of Science and Technology, Xiangtan 411201, China;
2. State Key Laboratory of Coal Resources and Safe Mining, China University of Mining & Technology, Beijing 100083, China;
3. Geophysical Research Institute of Bureau of Geophysical Prospecting, CNPC, Zhuozhou 072751, China
Abstract: Using an effective GMRES iterative method to solve two way acoustic wave equation in three-dimensional frequency-space domain, using multiscale and multigrid strategy in inversion procedure, we discuss how to quickly realize a high precision and high resolution 3D frequency-space domain iterative full waveform velocity inversion(fwi). We use this three-dimensional frequency-space domain acoustic wave full waveform velocity inversion method on a synthetic data. Numerical experiments demonstrate that this method has a high computational efficiency, the spatial resolution of inversion results is high. Which providing a good foundation for the three-dimensional frequency-space domain iterative acoustic wave full waveform velocity inversion on real field data.
Key words: multiscale     frequency-space domain     iterative method     three dimensional acoustic wave equation     full waveform velocity inversion    

0 引 言

随着地下勘探目标体变得越来越隐蔽,油气勘探变得愈发困难.因此,对高精度反演成像的需求也日益迫切.旅行时层析成像利用高频假设来近似波动方程,从而得到射线方程,并通过求解射线方程来获得地震波传播的路径以及相应的走时.因为利用了高频近似和只用走时信息来反演地下结构,旅行时层析成像反演结果的分辨率不高;而全波形反演使用高精度的数值算法求解双程波方程,利用波场的全部信息.因此,全波形反演比旅行时层析的计算成本高昂.但是全波形反演方法的空间分辨率远好于旅行时层析反演方法(Pratt and Worthington,1990).

Pratt开拓了在频率-空间域使用直接法进行二维全波形正反演方法研究(Pratt and Worthington,1990Pratt et al.,1998Pratt,1999Pratt and Shipp,1999Sirgue and Pratt,2004).(Operto et al.,2004Operto et al.,2006Shin and Min,2006Shin et al.,2007Bendar et al.,2007)等将二维频率-空间域全波形反演发展到粘声波、粘弹性波和各向异性介质,广泛应用于井间地震、VSP、海洋地震、地面地震资料中.国外的二维频率-空间域全波形反演已经成功地应用于油气工业生产中.

二维频率-空间域波动方程正反演数值计算时,最有效的方法是使用嵌套剖分进行直接法LU分解大规模稀疏系数矩阵,使用计算代价相对便宜的追赶法,进行基于单频共享存储的LU分解因子的多炮快速正反演计算,详见文献(廖建平,2011ab,c).

现今,油气工业界采集的地震资料大部分为三维资料.而且,随着计算机硬件的快速发展(包括相对廉价的集群的出现),使得很有必要进行三维全波形正反演研究.

三维全波形反演的主要困难在于:三维时间域波形反演数值计算时,需要很多CPU时间;三维频率-空间域全波形反演时,需要求解大规模稀疏系数矩阵,从而导致计算量庞大.

我们知道:求解线性方程组的经典的方法一般分为两类:直接法和迭代法.一般来说,直接法对于阶数比较低的方程组(少于20000至30000个未知数)比较有效;而迭代法对于比较大的方程组更有效.在实际计算中,几十万甚至几百万个未知数的方程组并不少见,譬如,较大尺寸的三维频率-空间域声波全波形正反演数值计算即为如此.在这些情况下,迭代法有无可比拟的优势.

直接法求解三维频率-空间域声波方程时,我们以立方体的模型为例,假设N是每一个方向的网格点数,则形成的稀疏系数矩阵的元素个数为N6,该矩阵的运算次数为N6,所需要的计算机存储空间为N4(Ben-Hadj-Ali et al.,2007tekl and Pain,2002tekl et al.,2007Warner et al.,2007).因此,现在的计算机内存无法容纳被操作的稀疏系数矩阵,这给直接法进行三维频率-空间域声波全波形正反演 数值计算带来很大的挑战.三维时的嵌套剖分法很复杂,对于多带宽的三维稀疏系数矩阵,它的计算效率并不高.而且,一旦矩阵被因子分解,其追赶法的计算代价也很大.所以,对于现实尺寸大小的模型(比如1000*1000*1000个网格点),三维频率-空间域直接法声波全波形反演的计算代价变得不切实际地大,远远超出当今计算机的硬件承受能力.

而迭代法对于比较大的方程组更有效,迭代法具有无可比拟的优势.另外,使用迭代法可以根据不同的精度要求选择终止时间,因此,迭代法比较灵活.

本文详细讨论了基于多尺度迭代求解三维频率-空间域声波方程的全波形速度反演,实现了稳健快速的三维频率-空间域迭代法声波全波形反演.理论模型的数值试验证实该方法计算效率高,为三维频率-空间域迭代法声波全波形反演实际地震资料打下了方法基础.

1 三维频率-空间域迭代法声波全波形速度反演理论
1.1 如何求解正演

频率-空间域全波形反演与其它反演一样,首先需要讨论如何求解正演问题.

三维频率-空间域声波方程为(Ben-Hadj-Ali et al.,2007)

式中ω为角频率,k(x,y,z)为体积模量,P(x,y,z,ω)为依赖频率的压力波场,ρ(x,y,z)为介质的密度,S(x,y,z,ω)为震源.

方程(1)写成矩阵的形式是AX=B(Ben-Hadj-Ali et al.,2007),其中A是大规模稀疏系数矩阵,取决于模型的速度、密度等,X是一个列向量,它包含每个网格点的单个频率的波场,B是一个类似的列向量,它包含每个网格点的震源.

采用tekl和 Pain (2002)的27点法的二阶有限差分格式,单位波长内采样5个网格点,其数字误差小于1%,能够确保产生高精度的数值模拟结果.

使用迭代法进行三维频率-空间域声波方程正演数值模拟(tekl et al.,2007Warner et al.,2007Sirgue et al.,2007tekl et al.,2008Warner et al,2008a2008b),其思想表述为:使用预条件矩阵(该预条件矩阵能够快速产生解)而近似求解矩阵A.这些近似的解再被矩阵A去获得一个有效的震源,然后,从震源B中减去该有效的震源,而获得剩余残差,将该残差作为震源,而使用在下一次迭代中,一直重复该过程,直到残差小于一个事先给定的小值.

迭代法求解的关键点是选择合适的预条件矩阵,和用于计算迭代更新残差的方法,以得到快速的收敛速度.我们使用以下的预条件:首先是在上、下、前、后和左、右这6个垂直方向,使用单程波传播算子,把能量传播若干个网格点.该预条件子稳定,可以得到一个快速收敛的解.然后,双程波方程迭代使用该近似解,从而获得双程波方程的解;第二个预条件是我们使用多重网格方法.首先求解粗网格上的波动方程,粗网格的解作为细网格的解的预条件子.

在数学上,广义最小残量方法(GMRES)是一个求解线性方程组数值解的迭代方法.该方法利用在Krylov子空间中有着最小残量的向量来逼近解.Arnoldi迭代方法被用来求解该向量.GMRES方法由Yousef Saad和Martin H. Schultz在1986年提出(Saad,2003).

本文使用广义最小残量方法(GMRES)(Saad,2003)迭代求解三维频率-空间域声波方程,只需要存储和使用前面的结果,从而减少了计算机内存使用量.迭代法反演所需的时间正比于炮数.

与直接法相比较(比如LU分解法),该方法的优点是内存需求减少了几个数量级.

与时间域方法相比较, 该方法的优点是:仅仅需要很少的几个离散频率进行层析成像(Sirgue and Pratt,2004),则对于一炮,频率-空间域方法所需要的计算时间,通常远少于时间域的方法所需要的时间,而且所需要的内存要少很多(Sirgue and Pratt,2004tekl et al.,2007).

三维频率-空间域迭代法声波全波形反演的缺点是:在 稀疏的频率域很难实现时窗,而时间域可以任意的截取时窗.

1.2 计算流程

使用广义最小残量方法进行三维频率-空间域迭代法声波全波形反演的计算流程如下所示:

1)建立初始速度模型;2)选择震源;3)选择频率(从第一个频率到最后一个频率),进行反演;4)正演模拟多炮;5)计算数据残差、反向传播残差;6)计算梯度、求解步长;7)更新模型,判断是否满足精度要求,如果满足精度要求则结束程序,否则,转到第3步,进行下一个频率的迭代反演;8)结束程序.

2 三维频率-空间域迭代法声波全波形速度反演数值试验

图 1是用于测试本反演方法的理论模型:模型的浅部含有一个高速异常体.模型的尺寸大小是X=160,Y=160,Z=50个网格点,网格点的间距为12.5 m.震源和检波器都设置在Z=10的深度.我们将真实模型的一维平均作为反演的初始速度模型.

图 1 160*160*50个网格点的真实速度模型 Fig. 1 True velocity model with 160*160*50 grids

由于低频分量的非线性没有高频的强烈,我们将多尺度的策略耦合在反演中,按照从低频到高频的顺序,进行反演,把低频的反演结果作为高频的初始速度模型,从而把非线性问题化为逐步线性问题,尽可能减少三维频率-空间域声波全波形速度反演解的非唯一性(即减轻全波形反演的强烈非线性).

我们选择初始频率为4 Hz,频率间隔为2 Hz,最大频率为20 Hz,每个频率迭代10次(能够确保本迭代法全波形反演收敛).

图 2是真实速度模型的一个垂直切片,图 3为本反演方法反演所得的该垂直切片.对比这两个垂直切片,我们可以看出异常体的空间位置、形状和大小都被很好地反演出来,很细微的速度变化都被成功反演恢复.由于速度的增大以及覆盖次数减少,所以,分辨率随着深度的增加而降低.真实速度模型和反演所恢复的速度二者之间的数值误差小于4%,由此可以看出全波形反演的分辨率高.

图 2 真实模型1的垂直切片 Fig. 2 Vertical slice of true model 1

图 3 三维频率-空间域迭代法声波全波形反演所得的图 2所示的垂直切片 Fig. 3 Result of 3D frequency-space domain iterative acoustic fwi on Fig. 2

图 4为初始模型底部的双程旅行时差,图 5为模型底部的双程旅行时差的反演结果,反演结果证实全波形反演减少了双程旅行时差,证明本方法的分辨率高.

图 4 初始模型底部的双程旅行时差 Fig. 4 Two way travel time difference of initial model bottom

图 5 模型底部的双程旅行时差的反演结果 Fig. 5 Fwi result of two way travel time difference of model bottom

为了验证本反演方法的分辨率高,我们将图 1所示的理论模型改变为包含两个高速异常体,观测系统与数值试验1相同,进行本反演方法的数值试验.

图 6为真实的速度模型的切片,其中,两个高速异常体间隔比较大.图 7图 6所示的两个高速异常体反演所恢复的速度模型的切片.

图 6 真实的模型2的切片 Fig. 6 Vertical slice of true model 2

图 7图 6反演所恢复的模型的切片 Fig. 7 Result of 3D frequency-space domain iterative acoustic fwi on Fig. 6

图 8为真实的速度模型的切片,其中,两个高速异常体的间隔小,部分重叠.图 9图 8所示的异常体反演所恢复的速度模型的切片.异常体的归位正确,异常体的细微特征都反演恢复得很好,其数值精确性近似于96%.因此,证实本反演方法的分辨率高.

图 8 真实的速度模型3的切片 Fig. 8 Vertical slice of true model 3

图 9图 8反演所恢复的速度模型的切片 Fig. 9 Result of 3D frequency-space domain iterative acoustic fwi on Fig. 8

整个反演过程中,计算机内存需求量为240兆,在4核奔4处理器的计算机上面,反演的时间为16小时.该模型的尺寸不太大,对于大尺寸的模型,我们可以使用大规模的计算机机群并行求解三维频率-空间域迭代法声波全波形速度反演.

3 结 论

通过理论模型的三维频率-空间域迭代法声波全波形反演数值试验,我们可以得出如下的结论和认识:

1)理论数据数值试验所得的反演结果证实,本反演方法对理论模型的分辨率高;

2)本反演方法稳健、计算速度快、计算效率高.

3)如何使得三维频率-空间域迭代法声波全波形速度反演能够成功反演三维地震实际资料,并且能够尽快地应用于油气工业生产,这需要进行持续深入的研究.

参考文献
[1] Bendar J B, Shin C, Pyun, S. 2007. Comparison of waveform inversion, part 2: Phase approach[J].   Geophysical Prospecting, 55(4): 465-475.
[2] George A, Liu J W H. 1981. Computer Solution of Large Sparse Positive Definite Systems[M]. New York: Prentice Hall.
[3] H. Ben Hadj Ali, S. Operto, J. Virieux, et al. 2007. 3D acoustic frequency-domain full waveform inversion. SEG Conference, San Antonio, (Extended abstract) .1730-1734.
[4] Liao J P, Wang H Z, Ma Z T. 2009. 2-D elastic wave modeling with frequency-space 25-point finite-difference operators[J].   Applied Geophysics, 6(3): 259-266.
[5] Liao J P, Wang H Z, Ma Z T. 2009. Frequency-space domain two dimension elastic wave model using compressed storage format[C]. SEG-CPS Annual Meeting. Beijing, 1180.
[6] Liao J P, Liu H X, Wang H Z, et al. 2011a. Study on rapid highly accurate acoustic wave numerical simulation in frequency space domain[J]. Progress in Geophys.   (in Chinese), 26(4): 1359-1363.
[7] Liao J P, Liu H X, Wang H Z, et al. 2011b. High resolution acoustic wave full waveform velocity inversion in frequency space domain-theoretical model[J]. Progress in Geophys.   (in Chinese), 26(5): 1690-1695.
[8] Liao J P, Liu H X, Wang H Z, et al. 2011c. Two dimensional Visco-acoustic wave modeling research in frequency-space domain[J]. Progress in Geophys.   (in Chinese), 26(6): 2075-2081.
[9] Operto, S., Ravaut, C., Improta, L., et al. 2004. Quantitative imaging of complex structures from multi-fold wide-aperture seismic data: a case study[J].   Geophysical Prospecting,52:625-651.
[10] Operto S, Virieux J, Dessa J X, et al. 2006. Crustal seismic imaging from multifold ocean bottom seismometer data by frequency domain full waveform tomography: Application to the eastern Nankai trough[J]. Journal of Geophysical Research, 111, B09306, doi:10.  1029/2005JB003835.
[11] Pratt R G, Worthington M H. 1990. Inverse theory applied to multi-source cross-hole tomography. PartⅠ: Acoustic wave-equation method[J].   Geophysical Prospecting, 38(3): 287-310.
[12] Pratt R G, Shin C, Hicks. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[J].   Geophysical Journal International, 133(2): 341-362.
[13] 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.
[14] Pratt R G, Shipp R M. 1999. Seismic waveform inversion in the frequency domain, Part 2: Fault delineation in sediments using crosshole data[J].   Geophysics, 64(3): 902-914.
[15] Saad Y. 2003. Iterative Methods for Sparse Linear Systems (2nd edition) [M].   Society for Industrial and Applied Mathematics, ISBN 978-0-89871-534-7.
[16] Shin C, Min D J. 2006. Waveform inversion using a logarithmic wavefield[J].   Geophysics, 71(3): R31-R42.
[17] Shin C, Pyn S, Bendar J B. 2007. Comparison of waveform inversion, part 1: conventional wave field vs logarithmic wave field[J].   Geophysical Prospecting, 55(4): 449-464.
[18] Sirgue L, Pratt R G. 2004. Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies[J].   Geophysics, 69(1): 231-248.
[19] Sirgue L, Etgen J, Albertin U. 2007. 3D full-waveform inversion: Wide- versus narrow-azimuth acquisitions[A].// Proceedings of the 77th Annual International Meeting SEG. pp.   1760-1764.
[20] tekl I, Pain C. 2002. 3D frequency domain visco-acoustic modeling using rotated finite difference operators.   Abstract, pages C-27 EAGE.
[21] tekl I, Warner M R, Umpleby A P. 2007. 3D Frequency Domain Waveform Inversion-Synthetic Shallow Channel Example[C]. 69th EAGE Conference & Exhibition (Extended abstract).
[22] tekl I, Warner M, Umpleby A, et al. 2008. Practical 3D wavefield tomography on field datasets[C]. 70th EAGE Conference & Exhibition.
[23] Warner M, tekl I, Umpleby A. 2007. Full wavefield seismic tomography-iterative forward modelling in 3D[C].   69th EAGE Conference & Exhibition.
[24] Warner M, tekl I, Umpleby A. 2008a. Efficient and Effective 3D Wavefield Tomography[C]. 70th EAGE F023.
[25] Warner M, tekl I, Umpleby A, et al. 2008b. 3D wavefield tomography: Problems, opportunities and future directions[C].  // 70th EAGE.
[26] 廖建平, 刘和秀, 王华忠,等. 2011a. 快速高精度的频率空间域声波数值模拟方法研究[J].   地球物理学进展, 26(4): 1359-1363.
[27] 廖建平, 刘和秀, 王华忠,等. 2011b. 高分辨率的频率空间域声波全波形速度反演-理论模型[J].   地球物理学进展, 26(5): 1690-1695.
[28] 廖建平, 刘和秀, 王华忠,等. 2011c. 二维频率空间域粘声波正演模拟研究[J].   地球物理学进展, 26(6): 2075-2081.