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
随着地下勘探目标体变得越来越隐蔽,油气勘探变得愈发困难.因此,对高精度反演成像的需求也日益迫切.旅行时层析成像利用高频假设来近似波动方程,从而得到射线方程,并通过求解射线方程来获得地震波传播的路径以及相应的走时.因为利用了高频近似和只用走时信息来反演地下结构,旅行时层析成像反演结果的分辨率不高;而全波形反演使用高精度的数值算法求解双程波方程,利用波场的全部信息.因此,全波形反演比旅行时层析的计算成本高昂.但是全波形反演方法的空间分辨率远好于旅行时层析反演方法(Pratt and Worthington,1990).
Pratt开拓了在频率-空间域使用直接法进行二维全波形正反演方法研究(Pratt and Worthington,1990;Pratt et al.,1998;Pratt,1999;Pratt and Shipp,1999;Sirgue and Pratt,2004).(Operto et al.,2004;Operto et al.,2006;Shin and Min,2006;Shin et al.,2007;Bendar et al.,2007)等将二维频率-空间域全波形反演发展到粘声波、粘弹性波和各向异性介质,广泛应用于井间地震、VSP、海洋地震、地面地震资料中.国外的二维频率-空间域全波形反演已经成功地应用于油气工业生产中.
二维频率-空间域波动方程正反演数值计算时,最有效的方法是使用嵌套剖分进行直接法LU分解大规模稀疏系数矩阵,使用计算代价相对便宜的追赶法,进行基于单频共享存储的LU分解因子的多炮快速正反演计算,详见文献(廖建平,2011a,b,c).
现今,油气工业界采集的地震资料大部分为三维资料.而且,随着计算机硬件的快速发展(包括相对廉价的集群的出现),使得很有必要进行三维全波形正反演研究.
三维全波形反演的主要困难在于:三维时间域波形反演数值计算时,需要很多CPU时间;三维频率-空间域全波形反演时,需要求解大规模稀疏系数矩阵,从而导致计算量庞大.
我们知道:求解线性方程组的经典的方法一般分为两类:直接法和迭代法.一般来说,直接法对于阶数比较低的方程组(少于20000至30000个未知数)比较有效;而迭代法对于比较大的方程组更有效.在实际计算中,几十万甚至几百万个未知数的方程组并不少见,譬如,较大尺寸的三维频率-空间域声波全波形正反演数值计算即为如此.在这些情况下,迭代法有无可比拟的优势.
直接法求解三维频率-空间域声波方程时,我们以立方体的模型为例,假设N是每一个方向的网格点数,则形成的稀疏系数矩阵的元素个数为N6,该矩阵的运算次数为N6,所需要的计算机存储空间为N4(Ben-Hadj-Ali et al.,2007;tekl and Pain,2002;tekl et al.,2007;Warner et al.,2007).因此,现在的计算机内存无法容纳被操作的稀疏系数矩阵,这给直接法进行三维频率-空间域声波全波形正反演 数值计算带来很大的挑战.三维时的嵌套剖分法很复杂,对于多带宽的三维稀疏系数矩阵,它的计算效率并不高.而且,一旦矩阵被因子分解,其追赶法的计算代价也很大.所以,对于现实尺寸大小的模型(比如1000*1000*1000个网格点),三维频率-空间域直接法声波全波形反演的计算代价变得不切实际地大,远远超出当今计算机的硬件承受能力.
而迭代法对于比较大的方程组更有效,迭代法具有无可比拟的优势.另外,使用迭代法可以根据不同的精度要求选择终止时间,因此,迭代法比较灵活.
本文详细讨论了基于多尺度迭代求解三维频率-空间域声波方程的全波形速度反演,实现了稳健快速的三维频率-空间域迭代法声波全波形反演.理论模型的数值试验证实该方法计算效率高,为三维频率-空间域迭代法声波全波形反演实际地震资料打下了方法基础.
频率-空间域全波形反演与其它反演一样,首先需要讨论如何求解正演问题.
三维频率-空间域声波方程为(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.,2007;Warner et al.,2007;Sirgue et al.,2007;tekl et al.,2008;Warner et al,2008a,2008b),其思想表述为:使用预条件矩阵(该预条件矩阵能够快速产生解)而近似求解矩阵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,2004;tekl et al.,2007).
三维频率-空间域迭代法声波全波形反演的缺点是:在 稀疏的频率域很难实现时窗,而时间域可以任意的截取时窗.
使用广义最小残量方法进行三维频率-空间域迭代法声波全波形反演的计算流程如下所示:
1)建立初始速度模型;2)选择震源;3)选择频率(从第一个频率到最后一个频率),进行反演;4)正演模拟多炮;5)计算数据残差、反向传播残差;6)计算梯度、求解步长;7)更新模型,判断是否满足精度要求,如果满足精度要求则结束程序,否则,转到第3步,进行下一个频率的迭代反演;8)结束程序.
图 1是用于测试本反演方法的理论模型:模型的浅部含有一个高速异常体.模型的尺寸大小是X=160,Y=160,Z=50个网格点,网格点的间距为12.5 m.震源和检波器都设置在Z=10的深度.我们将真实模型的一维平均作为反演的初始速度模型.
由于低频分量的非线性没有高频的强烈,我们将多尺度的策略耦合在反演中,按照从低频到高频的顺序,进行反演,把低频的反演结果作为高频的初始速度模型,从而把非线性问题化为逐步线性问题,尽可能减少三维频率-空间域声波全波形速度反演解的非唯一性(即减轻全波形反演的强烈非线性).
我们选择初始频率为4 Hz,频率间隔为2 Hz,最大频率为20 Hz,每个频率迭代10次(能够确保本迭代法全波形反演收敛).
图 2是真实速度模型的一个垂直切片,图 3为本反演方法反演所得的该垂直切片.对比这两个垂直切片,我们可以看出异常体的空间位置、形状和大小都被很好地反演出来,很细微的速度变化都被成功反演恢复.由于速度的增大以及覆盖次数减少,所以,分辨率随着深度的增加而降低.真实速度模型和反演所恢复的速度二者之间的数值误差小于4%,由此可以看出全波形反演的分辨率高.
图 4为初始模型底部的双程旅行时差,图 5为模型底部的双程旅行时差的反演结果,反演结果证实全波形反演减少了双程旅行时差,证明本方法的分辨率高.
为了验证本反演方法的分辨率高,我们将图 1所示的理论模型改变为包含两个高速异常体,观测系统与数值试验1相同,进行本反演方法的数值试验.
图 6为真实的速度模型的切片,其中,两个高速异常体间隔比较大.图 7为图 6所示的两个高速异常体反演所恢复的速度模型的切片.
图 8为真实的速度模型的切片,其中,两个高速异常体的间隔小,部分重叠.图 9为图 8所示的异常体反演所恢复的速度模型的切片.异常体的归位正确,异常体的细微特征都反演恢复得很好,其数值精确性近似于96%.因此,证实本反演方法的分辨率高.
整个反演过程中,计算机内存需求量为240兆,在4核奔4处理器的计算机上面,反演的时间为16小时.该模型的尺寸不太大,对于大尺寸的模型,我们可以使用大规模的计算机机群并行求解三维频率-空间域迭代法声波全波形速度反演.
通过理论模型的三维频率-空间域迭代法声波全波形反演数值试验,我们可以得出如下的结论和认识:
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. |