2017, Vol. 32
2. 煤炭资源与安全开采国家重点实验室(中国矿业大学), 北京 100083
3. 湖南科技大学地质系, 湘潭 411201
2. State Key Laboratory of Coal Resources and Safe Mining, China University of Mining & Technology, Beijing 100083, China
3. Department of Geology of Hunan University of Science and Technology, Xiangtan 411201, China
随着油气勘探、开发程度的不断深入，油气工业界对高精度成像和储层描述、地震解释提出了更高的要求，急需地震勘探方法技术，尤其是精确的反演成像方法技术的进步.地震波反演成像包含两个层次，即反射系数估计(仅仅定位反射系数出现的位置的叠前偏移成像)和全波形反演(Full Waveform Inversion，FWI)估计速度、密度、各向异性参数和吸收衰减等弹性参数.
Tarantola(1984)的工作是地震反演的奠基之作.成功地把偏移理论归纳为一个数学上的局部最优化问题，即把观测数据和模拟数据之间的残差指向最小的最小二乘问题.通过从震源激发的入射波场和反向传播的残差波场之间的互相关而建立残差函数的梯度，开启了地震波形反演的大门.当存在强烈的非线性时，反演可能收敛到一个局部极小值(Gauthier et al., 1986)而失败，因此，初始模型必须位于全局极小值的附近.
20世纪80年代，地震波反演处于蓬勃发展时期，多位学者(Tarantola, 1984；Mora, 1987, 1989)对FWI理论做了非常扎实的理论研究.但受当时野外地震观测技术和高性能计算机技术的局限，无法将提出的理论方法用于工业界.
近几年，随着高性能计算机技术的飞速进步以及全方位、宽频带、高密度的地震数据野外采集技术的飞速发展，以高精度、多参数建模为优点的反演估计全波数或宽波数带的弹性参数为目标的FWI方法再次成为当前勘探地震学的研究热点.理论上，FWI已被证明是一种高精度速度建模的有效手段(Tarantola，1984; Tarantola et al., 1986；Pratt, 1990, 1999; Pratt et al., 1998；Virieux and Operto, 2009；Virieux et al., 2017).但是，在数学上它是一个高度病态的非线性问题(Tarantola，1984)，需要解决其多解性及收敛性问题.从地球物理角度看，它涉及模型的参数化、误差泛函的建立、数据预处理、波场的数值模拟、子波的估计等研究内容(Tarantola，1984; Tarantola et al., 1986；Pratt et al., 1998; Pratt，1999；许琨和王妙月，2001；Virieux and Operto, 2009；卞爱飞等，2010；廖建平等, 2011a, b；Virieux et al., 2017).
Pratt开拓了在频率-空间域使用直接法进行二维全波形正反演方法研究(Pratt，1990；Pratt et al., 1998；Pratt，1999；Sirgue and Pratt, 2004).Charl-Hyun等(1996)利用9点法进行二维声波正演，提高了计算精度；Shin(1998)推广到25点声波正演，几何数量级地增加了内存需求，降低了数值模拟的计算速度；Štekl等(1998)实现了二维弹性波9点法正演；Min等(2000)实现了弹性波25点法正演；Min等(2003, 2004)讨论了自由表面等边界条件的频率空间域有限元正演；Hustedt等(2004)实现了交错网格和PML的二维频率空间域正演.综合比较计算精度、内存需求和计算速度，大部分研究人员采用9点法进行二维频率-空间域正反演.Operto等(2006)、Shin和Min(2006)、Bednar等(2007)将二维频率-空间域全波形反演发展到黏声波、黏弹性波和各向异性介质，广泛应用于井间地震、VSP、海洋地震、地面地震资料中.
刘璐等(2013)基于修正拟牛顿公式的全波形反演；魏哲枫等(2014)基于非规则网格声波正演的时间域全波形反演；刘玉柱等(2015)基于Born敏感核函数的VTI介质多参数全波形反演；张文生等(2015)频率多尺度全波形速度反演；王本锋等(2016)基于T-matrix的非线性参数估计方法；崔永福等(2016)全波形反演在缝洞型储层速度建模中的应用；罗静蕊等(2016)地震包络反演对局部极小值的抑制特性研究；郭雪豹等(2016)基于频域衰减的时域全波形反演；陈生昌和陈国新(2016)时间二阶积分波场的全波形反演；胡勇等(2017)基于精确震源函数的解调包络多尺度全波形反演；桂生等(2017)简化混合域全波形反演多GPU加速策略.
波形反演既可以在时间域实现(Gauthier等，1986；Mora，1987)又可以在频率域实现(Pratt and Worthington, 1990；Pratt et al., 1998).在反演策略上，由于时间域FWI可以灵活地对数据进行必要的预处理，以及灵活选取所需的特征波等特点而受到关注.但时间域同时反演所有频率成分，增大了反问题的非线性，因此，Bunks等(1995)提出了时间域的多尺度反演策略，通过对资料的处理将问题分解为不同的空间尺度，增加了数值计算时的稳定性，降低了非线性程度.
由于时间域的数据量大，需要大量的计算机内存进行存储，而当前的计算机集群往往不能满足要求，特别是在三维情况下，这一问题显得尤为突出.频率域反演使用几个离散频率的傅里叶级数取代庞大的时间采样点信息，大量节省了计算机内存资源.其次，频率域迭代算法从低频成分开始反演，然后到高频的反演，充分利用低频分量的线性性，而降低了全波形反演的非线性，相比较于时间域反演而言，频率域更利于实现这一算法.尤其是对衰减系数等物理参数的反演，频率域反演具有时间域反演无法比拟的优越性.多尺度策略能够减轻全波形反演的非线性.频率域能够容易地实现基于任意带宽和采样间隔的频率组的多尺度全波形反演，而时间域对于任意观测系统都可以灵活运用时窗.
本文讨论了二维频率空间域直接法声波波动方程全波形速度反演和二维时间空间域声波波动方程全波形速度反演，并对它们进行了理论模型的数值试验，通过对数值试验结果的分析，最后得出了相应的结论.
1 二维时间空间域和频率空间域声波全波形反演理论二维常密度频率空间域声波方程用矩阵形式表达为(Pratt等，1998)：
|
(1) |
式中ω为角频率，u(ω)是傅里叶变换后的离散波场，网格点的个数为l，S(ω)是一个具有l×l个元素的阻抗矩阵，其元素为复数，f(ω)为震源项，u(ω)和f(ω)是l×1的列矢量.矩阵S(ω)为高度稀疏、非对称、非正定，对它直接求逆非常困难.通常采用直接法因式分解求解该稀疏矩阵.我们采用Charl-Hyun等(1996)的二阶有限差分格式，进行二维频率空间域直接法声波全波形反演.
假设模型参数的个数为m，根据m×l的列矢量参数集p定义模型，通常m < l，我们从参数集p中计算l×l的阻抗矩阵S(ω).
假设具有n个试验观测d，则方程(1) 中的每一项都依赖于角频率ω，我们通过减少观测数据和模拟数据之间的残差函数，而不断迭代更新模型，实现反演.认为存在一个足够接近于全局解的初始模型p(0)，给定初始模型，通过正演模拟得到模拟波场值.
定义剩余残差δd为初始模型响应和观测数据之间的差，即：
|
(2) |
式中，下标i代表检波器的个数，我们寻求最小化波场残差的平方和，公式为
|
(3) |
即误差函数的l2范数，其中的T为矩阵转置，*为复共轭.详细的理论推导请见参考文献(Pratt et al., 1998).
本文使用L-BFGS算法(即限定内存的BFGS算法)来求解模型更新量，进行二维时间空间域声波全波形反演.L-BFGS方法只需要存储有限个数的模型量，在每次迭代中只需要输入目标函数的梯度，通过保存最近几次的曲率信息来更新近似矩阵的这种方法在实践中是很有效的.虽然L-BFGS算法是线性收敛，但是每次迭代的开销非常小，因此, L-BFGS算法运算速度很快，而且由于每一步迭代都能保证近似海森矩阵正定，因此，算法的稳健性很强.
我们使用经典的伴随状态法计算目标函数的梯度.这种方法使用炮点正向传播入射波场与检波点逆时传播残差波场的互相关生成梯度，避免了直接计算Frechet导数.
2 二维时间空间域和频率空间域声波全波形反演的数值试验图 1是Marmousi模型的真实速度模型.图 2是对真实速度模型图 1进行平滑而得到的模型作为初始速度模型.
|
图 1 真实的Marmousi速度模型 Figure 1 The real Marmousi velocity model |
|
图 2 平滑后的Marmousi模型作为反演的初始速度模型 Figure 2 The starting Marmousi velocity model |
我们使用该平滑速度模型作为初始速度模型，分别进行二维时间空间域声波全波形反演和二维频率空间域直接法声波全波形反演数值试验.
为了验证二维时间空间域和频率空间域直接法声波全波形反演的优缺点，我们在作数值试验时，都采用相同的观测系统和参数.观测系统如下所示：纵向为279个采样点，横向为1359个采样点，纵向和横向的网格间距都为10 m，一共有130炮，第一炮的位置为纵向10 m，横向10 m，炮点间距为纵向0 m，横向为10 m，检波点数为1358，检波点的初始位置为纵向10 m和横向10 m，检波点间距为纵向0 m，横向10 m.
我们进行二维时间空间域声波全波形反演(使用L-BFGS算法)数值试验时，采用多尺度的策略，从低频到高频进行反演，将低频的反演结果作为高频反演时的初始速度模型，依次进行主频为5 Hz、8 Hz和12 Hz的波形反演数值试验.数值试验反演得到的速度分别为图 3、图 4和图 5.我们从反演所恢复出的速度模型中可以看出，随着反演的频率增大，反演所得速度模型的复杂构造细节愈发清晰可辨，说明反演的分辨率越高.
|
图 3 主频为5 Hz的二维时间空间域声波全波形反演后得到的速度模型 Figure 3 The inverted result of Marmousi velocity using 5 Hz in time-space domain |
|
图 4 主频为8 Hz的二维时间空间域声波全波形反演后得到的速度模型 Figure 4 The inverted result of Marmousi velocity using 8 Hz in time-space domain |
|
图 5 主频为12 Hz的二维时间空间域声波全波形反演后得到的速度模型 Figure 5 The inverted result of Marmousi velocity using 12 Hz in time-space domain |
二维频率空间域直接法声波全波形反演数值试验时，我们对该初始速度模型分别进行5 Hz、8 Hz和12 Hz的二维频率空间域直接法声波全波形反演数值试验.反演所得到的速度分别展示在图 6、图 7和图 8.我们从反演所恢复出的速度模型中可以看出，随着反演的频率增大，反演所得的模型的复杂构造细节愈发清晰可辨，说明反演的分辨率越高.
|
图 6 5 Hz的二维频率空间域声波全波形反演后得到的速度模型 Figure 6 The inverted result of Marmousi velocity using 5 Hz in frequency-space domain |
|
图 7 8 Hz的二维频率空间域声波全波形反演后得到的速度模型 Figure 7 The inverted result of Marmousi velocity using 8 Hz in frequency-space domain |
|
图 8 12 Hz的二维频率空间域声波全波形反演后得到的速度模型 Figure 8 The inverted result of Marmousi velocity using 12 Hz in frequency-space domain |
我们将二维时间空间域和频率空间域声波全波形反演所恢复的速度模型进行对比，分析二维时间空间域和频率空间域直接法声波全波形反演的优缺点.
我们使用14个计算节点(130个进程)进行二维时间空间域声波全波形反演数值试验，全波形反演该模型所需要的总时间为119172 s，需要的内存为9918 MB.使用5 Hz、8 Hz和12 Hz三组频率，每组频率都进行25次迭代.
二维频率空间域直接法声波全波形反演数值试验时，我们使用14个计算节点，共使用130个进程，全波形反演该模型所需要的总计算时间为68510 s，需要的内存为6470 MB.使用5 Hz、8 Hz和12 Hz三组频率，每组频率都进行25次迭代.
二维频率空间域使用更多的计算节点时，其加速有限(主要用于矩阵的LU分解，炮点计算波场为线性关系)，二维时间空间域更灵活，多炮同时计算，可以多倍提高计算效率.
使用相同的频率进行反演时，二维时间空间域声波全波形反演所得到的速度的分辨率比二维频率空间域直接法声波全波形反演数值试验时反演得到的速度的分辨率更高.二维时间空间域声波全波形反演，使用12 Hz反演时，所得到的速度和真实的理论模型相比较，二者之间的差别不是很大，复杂构造的形态一致，特别是在浅层，二者几乎完全相同.在中深层处，二者之间的差异稍微大些，这主要是由于我们所使用的频率仅仅为12 Hz.数值试验的结果证明了本文反演方法的分辨率较高.
3 结论 3.1通过二维时间空间域和频率空间域声波全波形反演Marmousi理论模型的数值试验，我们可以得出如下的结论和认识：
(1) 使用更多的节点进行二维频率空间域直接法声波全波形反演时，其加速有限(正演模拟的计算量主要用于矩阵的LU分解，炮点计算波场为线性关系).
(2) 时间空间域更灵活，多炮同时计算，可以多倍提高计算效率.
(3) 二维声波全波形反演时，频率空间域的计算速度远快于时间空间域，需要的内存也比时间空间域少.
3.2总之，二维声波全波形反演时，相比较于时间空间域方法，频率空间域直接法声波全波形反演具有计算速度快和节省计算机内存需求的优势.
致谢 感谢国家留学基金委和湖南省地方合作项目对本项研究的资助.| [] | Bednar J B, Shin C, Pyun S. 2007. Comparison of waveform inversion, part 2:Phase approach[J]. Geophysical Prospecting, 55(4): 465–475. DOI:10.1111/gpr.2007.55.issue-4 |
| [] | Bian A F, Yu W H, Zhou H W. 2010. Progress in the frequency-domain full waveform inversion method[J]. Progress in Geophys. (in Chinese), 25(3): 982–993. DOI:10.3969/j.issn.1004-2903.2010.03.037 |
| [] | Bunks C, Saleck F M, Zaleski S, et al. 1995. Multiscale seismic waveform inversion[J]. Geophysics, 60(5): 1457–1473. DOI:10.1190/1.1443880 |
| [] | Charl-Hyun J, 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 |
| [] | Chen S C, Chen G X. 2016. Full waveform inversion of the second-order time integral wavefield[J]. Chinese Journal of Geophysics (in Chinese), 59(10): 3765–3776. DOI:10.6038/cjg20161021 |
| [] | Cui Y F, Peng G X, Wu G C, et al. 2016. Application of full waveform inversion velocity model-building technology for the fractured-vuggy reservoir[J]. Chinese Journal of Geophysics (in Chinese), 59(7): 2713–2725. DOI:10.6038/cjg20160734 |
| [] | Gauthier O, Virieux J, Tarantola A. 1986. Two-dimensional nonlinear inversion of seismic waveforms:Numerical results[J]. Geophysics, 51(7): 1387–1403. DOI:10.1190/1.1442188 |
| [] | Gui S, Liu H, Zhang Y J. 2017. The acceleration strategy of a simplified hybrid domain full waveform inversion on multi-GPUs[J]. Chinese Journal of Geophysics (in Chinese), 60(2): 665–677. DOI:10.6038/cjg20170220 |
| [] | Guo X B, Liu H, Shi Y. 2016. Time domain full waveform inversion based on frequency attenuation[J]. Chinese Journal of Geophysics (in Chinese), 59(10): 3777–3787. DOI:10.6038/cjg20161022 |
| [] | Hu Y, Han L G, Xu Z, et al. 2017. Demodulation envelope multi-scale full waveform inversion based on precise seismic source function[J]. Chinese Journal of Geophysics (in Chinese), 60(3): 1088–1105. DOI:10.6038/cjg20170321 |
| [] | 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 |
| [] | 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. DOI:10.3969/j.issn.1004-2903.2011.04.029 |
| [] | 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. DOI:10.3969/j.issn.1004-2903.2011.05.023 |
| [] | Liu L, Liu H, Zhang H, et al. 2013. Full waveform inversion based on modified quasi-Newton equation[J]. Chinese Journal of Geophysics (in Chinese), 56(7): 2447–2451. DOI:10.6038/cjg20130730 |
| [] | Liu Y Z, Wang G Y, Yang J Z, et al. 2015. Multi-parameter full-waveform inversion for VTI media based on Born sensitivity kernels[J]. Chinese Journal of Geophysics (in Chinese), 58(4): 1305–1316. DOI:10.6038/cjg20150418 |
| [] | Luo J R, Wu R S, Gao J H. 2016. Local minima reduction of seismic envelope inversion[J]. Chinese Journal of Geophysics (in Chinese), 59(7): 2510–2518. DOI:10.6038/cjg20160716 |
| [] | Mora P. 1987. Nonlinear two-dimensional elastic inversion of multioffset seismic data[J]. Geophysics, 52(9): 1211–1228. DOI:10.1190/1.1442384 |
| [] | Mora P. 1989. Inversion=migration+tomography[J]. Geophysics, 54(12): 1575–1586. DOI:10.1190/1.1442625 |
| [] | 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(B9): B09306. DOI:10.1029/2005JB003835 |
| [] | Pratt R G. 1990. Inverse theory applied to multi-source cross-hole tomography.Part 2:Elastic wave-equation method[J]. Geophysical Prospecting, 38(3): 311–329. DOI:10.1111/gpr.1990.38.issue-3 |
| [] | 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 |
| [] | Pratt R G, Shin C, Hick G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[J]. Geophysical Journal International, 133(2): 341–362. DOI:10.1046/j.1365-246X.1998.00498.x |
| [] | Pratt R G, Worthington M H. 1990. Inverse theory applied to multi-source cross-hole tomography.Part 1:Acoustic wave-equation method[J]. Geophysical Prospecting, 38(3): 287–310. DOI:10.1111/gpr.1990.38.issue-3 |
| [] | Shin C, Min D J. 2006. Waveform inversion using a logarithmic wavefield[J]. Geophysics, 71(3): R31–R42. DOI:10.1190/1.2194523 |
| [] | 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 |
| [] | Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation[J]. Geophysics, 49(8): 1259–1266. DOI:10.1190/1.1441754 |
| [] | Virieux J, Asnaashari A, Brossier R, et al. 2017. An introduction to full waveform inversion[A].//Grechka V, Wapenaar K eds. Encyclopedia of Exploration Geophysics[M]. Tulsa, OK:SEG, R1-1-R1-40. |
| [] | Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics[J]. Geophysics, 74(6): WCC1–WCC26. DOI:10.1190/1.3238367 |
| [] | Wang B F, Wu R S, Chen X H, et al. 2016. Non-linear parameter estimation method based on T-matrix[J]. Chinese Journal of Geophysics (in Chinese), 59(6): 2257–2265. DOI:10.6038/cjg20160628 |
| [] | Wei Z F, Gao H W, Zhang J F. 2014. Time-domain full waveform inversion based on an irregular-grid acoustic modeling method[J]. Chinese Journal of Geophysics (in Chinese), 57(2): 586–594. DOI:10.6038/cjg20140222 |
| [] | Xu K, Wang M Y. 2001. Finite element inversion of the coefficients of acoustic equation in frequency domain[J]. Chinese Journal of Geophysics (in Chinese), 44(6): 852–864. DOI:10.3321/j.issn:0001-5733.2001.06.015 |
| [] | Zhang W S, Luo J, Teng J W. 2015. Frequency multiscale full-waveform velocity inversion[J]. Chinese Journal of Geophysics (in Chinese), 58(1): 216–228. DOI:10.6038/cjg20150119 |
| [] | 卞爱飞, 於文辉, 周华伟. 2010. 频率域全波形反演方法研究进展[J]. 地球物理学进展, 25(3): 982–993. DOI:10.3969/j.issn.1004-2903.2010.03.037 |
| [] | 陈生昌, 陈国新. 2016. 时间二阶积分波场的全波形反演[J]. 地球物理学报, 59(10): 3765–3776. DOI:10.6038/cjg20161021 |
| [] | 崔永福, 彭更新, 吴国忱, 等. 2016. 全波形反演在缝洞型储层速度建模中的应用[J]. 地球物理学报, 59(7): 2713–2725. DOI:10.6038/cjg20160734 |
| [] | 桂生, 刘洪, 张玉洁. 2017. 简化混合域全波形反演多GPU加速策略[J]. 地球物理学报, 60(2): 665–677. DOI:10.6038/cjg20170220 |
| [] | 郭雪豹, 刘洪, 石颖. 2016. 基于频域衰减的时域全波形反演[J]. 地球物理学报, 59(10): 3777–3787. DOI:10.6038/cjg20161022 |
| [] | 胡勇, 韩立国, 许卓, 等. 2017. 基于精确震源函数的解调包络多尺度全波形反演[J]. 地球物理学报, 60(3): 1088–1105. DOI:10.6038/cjg20170321 |
| [] | 廖建平, 刘和秀, 王华忠, 等. 2011a. 快速高精度的频率空间域声波数值模拟方法研究[J]. 地球物理学进展, 26(4): 1359–1363. DOI:10.3969/j.issn.1004-2903.2011.04.029 |
| [] | 廖建平, 刘和秀, 王华忠, 等. 2011b. 高分辨率的频率空间域声波全波形速度反演-理论模型[J]. 地球物理学进展, 26(5): 1690–1695. DOI:10.3969/j.issn.1004-2903.2011.05.023 |
| [] | 刘玉柱, 王光银, 杨积忠, 等. 2015. 基于Born敏感核函数的VTI介质多参数全波形反演[J]. 地球物理学报, 58(4): 1305–1316. DOI:10.6038/cjg20150418 |
| [] | 刘璐, 刘洪, 张衡, 等. 2013. 基于修正拟牛顿公式的全波形反演[J]. 地球物理学报, 56(7): 2447–2451. DOI:10.6038/cjg20130730 |
| [] | 罗静蕊, 吴如山, 高静怀. 2016. 地震包络反演对局部极小值的抑制特性[J]. 地球物理学报, 59(7): 2510–2518. DOI:10.6038/cjg20160716 |
| [] | 王本锋, 吴如山, 陈小宏, 等. 2016. 基于T-matrix的非线性参数估计方法[J]. 地球物理学报, 59(6): 2257–2265. DOI:10.6038/cjg20160628 |
| [] | 魏哲枫, 高红伟, 张剑锋. 2014. 基于非规则网格声波正演的时间域全波形反演[J]. 地球物理学报, 57(2): 586–594. DOI:10.6038/cjg20140222 |
| [] | 许琨, 王妙月. 2001. 声波方程频率域有限元参数反演[J]. 地球物理学报, 44(6): 852–864. DOI:10.3321/j.issn:0001-5733.2001.06.015 |
| [] | 张文生, 罗嘉, 滕吉文. 2015. 频率多尺度全波形速度反演[J]. 地球物理学报, 58(1): 216–228. DOI:10.6038/cjg20150119 |