2017, Vol. 32
2. 煤炭资源与安全开采国家重点实验室(中国矿业大学), 北京 100083
3. 湖南科技大学地质系, 湖南湘潭 411201
4. School of Science and Sport, University of the West of Scotland, Paisley PA1 2BE, UK
2. State Key Laboratory of Coal Resources and Safe Mining, China University of Mining and Technology, Beijing 100083, China
3. Department of Geology of Hunan University of Science and Technology, Hunan Xiangtan 411201, China
4. School of Science and Sport, University of the West of Scotland, Paisley PA1 2BE, UK
随着油气勘探开发的不断深入，地下勘探目标体变得越来越隐蔽，油气勘探变得越来越困难.因此，对高精度反演成像的需求也日益迫切.需要不断寻找新的方法和技术，以适应当前油气勘探开发的需求.全波形反演(FWI)就是顺应这种需求而迅速发展起来的一项新技术，并且已经成为当今地球物理界的研究热点.走时层析反演利用高频假设来近似波动方程，从而得到射线方程，并通过求解射线方程来获得地震波传播的路径以及相应的走时.走时层析反演是一种计算地震速度长波长的重要方法，由于其数值计算的相对高效性, 如今广泛应用于各种尺度，从全球研究整个地幔，到近地表环境和工程研究，以及油气勘探中(Luo and Schuster, 1991; Ramachandran et al., 2011; Baumann-Wilke et al., 2012; Zelt et al., 2013; Chen and Zelt, 2016).
因为利用了高频近似和只用走时信息来反演地下结构，走时层析反演结果的分辨率不高；全波形反演利用波场的全部信息，使用高精度的数值算法求解双程波方程.因此，全波形反演比走时层析反演的计算成本高，但是全波形反演方法的空间分辨率远好于走时层析反演方法.
Tarantola(1984)的工作是地震反演的奠基之作.20世纪80年代，地震波反演处于蓬勃发展时期，多位学者(Tarantola, 1984；Mora, 1987, 1989)对FWI理论进行了非常扎实的理论研究.但受当时野外地震观测技术和高性能计算机技术的局限，人们无法将提出的理论方法应用于工业界.
近年来，随着高性能计算机技术的飞速进步以及全方位、宽频带、高密度的地震数据野外采集技术的飞速发展，以高精度、多参数建模为优点的反演估计全波数或宽波数带的弹性参数为目标的FWI方法再次成为当前勘探地震学的研究热点.
FWI是在正则化约束下，通过减小计算模拟数据和实际观测数据之间的误差，逐步逼近真实模型，迭代更新初始模型(Tarantola，1984).理论上，FWI已经被证明为一种高精度速度建模的有效手段(Tarantola，1984；Pratt, 1990, 1999; Pratt et al., 1998；Virieux and Operto, 2009；Virieux et al., 2017).但是，在数学上它是一个高度病态的非线性问题(Tarantola，1984)，需要解决其多解性及收敛性问题.从地球物理角度看，它涉及模型的参数化、误差泛函的建立、数据预处理、波场的数值模拟、子波的估计等研究内容(Tarantola，1984；Pratt et al., 1998; Pratt，1999；Virieux and Operto, 2009；Virieux et al., 2017).FWI的理论体系完美，但它对数据、初始模型、子波、地震波正演算子是有理论假设的.尽管FWI某些理论模型试验时可以得到几乎和真实模型完全匹配的反演结果.但是，在复杂介质情况下，理论模型数据的FWI反演结果不收敛，实际地震资料的FWI愈发困难.FWI反演结果很难收敛到真实解上，核心问题在于叠前地震数据的不完备(譬如有限观测孔径、有限带宽(尤其是缺失低频)，空间不规则的观测方式，地表不一致的观测，以及观测时具有较强烈的噪声)、地震波波场(尤其是反射地震波场)与反演参数之间具有严重的非线性关系，误差泛函存在非常多的局部极值点.因此，需要很好的初始模型来降低误差泛函的非线性性.但是，在复杂介质和低信噪比情况下，难以获取较高精度的初始模型；正演模型的不完善(地震波正演模拟算子无法很好地模拟实际地震波场中各种复杂的波现象)、初始模型距离反演的“真实”参数模型太远、地震子波未知且空变(愈发加重了地震波场与反演参数之间的非线性关系)等多种因素制约.
波形反演可以在时间域实现(Gauthier et al., 1986；Mora，1987；Bunks et al., 1995；Zhou et al., 1995；许琨和王妙月，2001；卞爱飞等，2010；廖建平等, 2011a, b；Yang et al., 2016).近年来，国内的时间-空间域声波全波形反演也取得了很多进展.胡光辉等(2013)实现了3D声波FWI及应用；刘璐等(2013)基于修正拟牛顿公式的全波形反演；魏哲枫等(2014)基于非规则网格声波正演的时间域全波形反演；刘玉柱等(2015)基于Born敏感核函数的VTI介质多参数全波形反演；王本锋等(2016)基于T-matrix的非线性参数估计方法；崔永福等(2016)全波形反演在缝洞型储层速度建模中的应用；罗静蕊等(2016)地震包络反演对局部极小值的抑制特性研究；郭雪豹等(2016)基于频域衰减的时域全波形反演；陈生昌和陈国新(2016)时间二阶积分波场的全波形反演；胡勇等(2017)基于精确震源函数的解调包络多尺度全波形反演.
在反演策略上，由于时间域FWI可以灵活地对数据进行必要的预处理，以及灵活选取所需的特征波等特点而受到关注.但时间域同时反演所有频率成分，增大了反问题的非线性.因此，Bunks等(1995)提出了时间域的多尺度反演，通过对资料的处理将问题分解为不同的空间尺度，增加了问题求解的稳定性，降低了非线性程度.
时间域具有对于任意观测系统都可以灵活运用时窗，计算方便等优点.本文讨论了L-BFGS算法的二维时间-空间域声波波动方程全波形速度反演，并进行了复杂介质Marmousi理论模型的数值试验，通过对数值试验的分析，最后得出相应的结论.
1 二维时间-空间域声波全波形速度反演理论定义初始模型响应和观测数据之差为剩余残差δd，即：
|
(1) |
式中下标i代表检波器的个数，我们求解最小化波场残差的平方和为
|
(2) |
即误差函数的L2范数，其中的T为矩阵转置，*为矩阵复共轭运算.
梯度计算是全波形反演最核心的部分.如何高效率地计算梯度，决定了反演的计算量.为了避免Frechet矩阵导数的计算所需要花费的巨大代价，我们使用伴随状态法计算梯度(Plessix，2006).该方法使用炮点正向传播的波场与检波点逆时传播的残差波场的零延迟互相关而生成误差泛函的梯度，避免了直接计算Frechet导数.此外，伪保守形式的一阶速度-压力场声波波动方程的使用(胡光辉等，2013)，使梯度计算更加简单有效.
我们使用伴随状态法来计算梯度，其计算公式为：
|
(3) |
式中：mi和λ*分别表示模型参数和反向传播的残差波场.在计算炮点正向传播波场与检波点逆时传播的残差波场时，可以使用同一个正演算子进行求解，由于正演算子的对称性，在求取对角散射内核
我们采用区域分解(物理上分割模型，从而可以节约单个CPU内核的内存使用量)和炮并行(炮并行能够最大限度地加速运算)的双并行算法，进行L-BFGS算法的二维时间-空间域声波全波形速度反演.
2 二维时间-空间域声波全波形速度反演的数值试验图 1是Marmousi模型的真实速度模型.由于本文是进行理论模型的全波形反演数值试验，因此，我们对真实速度模型进行平滑而得到的速度模型，作为全波形反演时的初始速度模型.我们使用图 2所示的初始速度模型进行二维时间-空间域声波全波形速度反演数值试验.考察图 1和图 2，我们可以看出真实速度模型的构造非常复杂，初始速度模型已经比较难以看到复杂的构造细节.
|
图 1 真实的Marmousi速度模型 Figure 1 The real Marmousi velocity model |
|
图 2 Marmousi模型平滑后作为初始速度模型 Figure 2 The starting Marmousi velocity model |
为了验证二维时间-空间域声波全波形速度反演的优缺点，我们采用如下所示的观测系统进行全波形反演数值试验：纵向采样点数为279，横向采样点数为1359，纵向、横向的网格间距都为10 m；总共130炮，第1炮位于(10 m，10 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赫兹的二维时间-空间域声波全波形反演得到的速度模型 Figure 3 The inverted result of Marmousi using 5 Hz in time space domain |
|
图 4 主频为8赫兹的二维时间-空间域声波全波形反演得到的速度模型 Figure 4 The inverted result of Marmousi using 8 Hz in time space domain |
|
图 5 主频为12赫兹的二维时间-空间域声波全波形反演后得到的速度模型 Figure 5 The inverted result of Marmousi velocity using 12 Hz in time space domain |
从反演所恢复的速度模型，我们分析二维时间-空间域声波全波形反演的优缺点.二维时间-空间域声波全波形反演数值试验时：我们使用14个计算节点，总共使用130个进程，全波形反演该模型所需要的总时间为33小时6分钟，所需要的内存为9918 M.5 Hz、8 Hz和12 Hz三组频率迭代计算时，每组频率都使用25次迭代.时间-空间域计算灵活，多炮同时计算，使用更多的计算节点时，可以多倍提高计算速度.
当使用12 Hz进行二维时间-空间域声波全波形反演时，反演所恢复得到的速度和真实的理论模型相比较，二者之间的差别不太大，大致形态基本一致.特别是在浅层，二者几乎相同.在中深层处，二者之间的差异稍微大些，这主要是由于我们所使用的频率仅仅为12 Hz.数值试验的结果证明了本文反演方法的分辨率较高.
3 结论通过二维时间-空间域声波全波形反演Marmousi理论模型的数值试验，我们得出如下的结论：
(1) 时间-空间域计算灵活，适用于各种采集观测系统，多炮同时计算，使用更多的计算节点时，可以多倍提高计算效率；
(2) 二维时间-空间域声波全波形速度反演理论模型的分辨率较高.
致谢 感谢国家留学基金委和湖南省地方合作项目对本项研究的资助，感谢审稿专家和编辑的宝贵意见.| [] | Baumann-Wilke M, Bauer K, Schovsbo N H, et al. 2012. P-wave traveltime tomography for a seismic characterization of black shales at shallow depth on Bornholm, Denmark[J]. Geophysics, 77(5): EN53–EN60. DOI:10.1190/geo2011-0326.1 |
| [] | Bian A F, Yu W H, Zhou H W. 2010. Progress in the frequency-domain full waveform inversion method[J]. Progress in Geophysics, 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 |
| [] | Chen J X, Zelt C A. 2016. Application of frequency-dependent traveltime tomography and full waveform inversion to realistic near-surface seismic refraction data[J]. Journal of Environmental and Engineering Geophysics, 21(1): 1–12. DOI:10.2113/JEEG21.1.1 |
| [] | Chen S C, Chen G X. 2016. Full waveform inversion of the second-order time integral wavefield[J]. Chinese Journal of Geophysics, 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, 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 |
| [] | Guo X B, Liu H, Shi Y. 2016. Time domain full waveform inversion based on frequency attenuation[J]. Chinese Journal of Geophysics, 59(10): 3777–3787. DOI:10.6038/cjg20161022 |
| [] | Hu G H, Jia C M, Xia H R, et al. 2013. Implementation and validation of 3D acoustic full waveform inversion[J]. Geophysical Prospecting for Petroleum, 52(4): 417–425. |
| [] | 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, 60(3): 1088–1105. DOI:10.6038/cjg20170321 |
| [] | Liao J P, Liu H X, Wang H Z, et al. 2011a. High resolution acoustic wave full waveform velocity inversion in frequency space domain-theoretical model[J]. Progress in Geophysics, 26(5): 1690–1695. DOI:10.3969/j.issn.1004-2903.2011.05.023 |
| [] | Liao J P, Liu H X, Wang H Z, et al. 2011b. Study on rapid highly accurate acoustic wave numerical simulation in frequency space domain[J]. Progress in Geophysics, 26(4): 1359–1363. DOI:10.3969/j.issn.1004-2903.2011.04.029 |
| [] | Liu L, Liu H, Zhang H, et al. 2013. Full waveform inversion based on modified quasi-Newton equation[J]. Chinese Journal of Geophysics, 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, 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, 59(7): 2510–2518. DOI:10.6038/cjg20160716 |
| [] | Luo Y, Schuster G T. 1991. Wave-equation traveltime inversion[J]. Geophysics, 56(5): 645–653. DOI:10.1190/1.1443081 |
| [] | 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 |
| [] | Plessix R E. 2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[J]. Geophys. J. Int., 167(2): 495–503. DOI:10.1111/gji.2006.167.issue-2 |
| [] | 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]. Geophys. J. Int., 133(2): 341–362. DOI:10.1046/j.1365-246X.1998.00498.x |
| [] | Ramachandran K, Bellefleur G, Brent T, et al. 2011. Imaging permafrost velocity structure using high resolution 3D seismic tomography[J]. Geophysics, 76(5): B187–B198. DOI:10.1190/geo2010-0353.1 |
| [] | Virieux J, Asnaashari A, Brossier R, et al. 2017. An introduction to full waveform inversion[M].//Grechka V, Wapenaar K. Encyclopedia of Exploration Geophysics. 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, 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, 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, 44(6): 852–864. DOI:10.3321/j.issn:0001-5733.2001.06.015 |
| [] | Yang P L, Brossier R, Métivier L, et al. 2016. A review on the systematic formulation of 3-D multiparameter full waveform inversion in viscoelastic medium[J]. Geophys. J. Int., 207(1): 129–149. DOI:10.1093/gji/ggw262 |
| [] | Zelt C A, Haines S, Powers M H, et al. 2013. Blind test of methods for obtaining 2-D near-surface seismic velocity models from first-arrival traveltimes[J]. Journal of Environmental and Engineering Geophysics, 18(3): 183–194. DOI:10.2113/JEEG18.3.183 |
| [] | Zhou C X, Cai W Y, Luo Y, et al. 1995. Acoustic wave-equation traveltime and waveform inversion of crosshole seismic data[J]. Geophysics, 60(3): 765–773. DOI:10.1190/1.1443815 |
| [] | 卞爱飞, 於文辉, 周华伟. 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 |
| [] | 郭雪豹, 刘洪, 石颖. 2016. 基于频域衰减的时域全波形反演[J]. 地球物理学报, 59(10): 3777–3787. DOI:10.6038/cjg20161022 |
| [] | 胡光辉, 贾春梅, 夏洪瑞, 等. 2013. 三维声波全波形反演的实现与验证[J]. 石油物探, 52(4): 417–425. |
| [] | 胡勇, 韩立国, 许卓, 等. 2017. 基于精确震源函数的解调包络多尺度全波形反演[J]. 地球物理学报, 60(3): 1088–1105. DOI:10.6038/cjg20170321 |
| [] | 廖建平, 刘和秀, 王华忠, 等. 2011a. 高分辨率的频率空间域声波全波形速度反演-理论模型[J]. 地球物理学进展, 26(5): 1690–1695. DOI:10.3969/j.issn.1004-2903.2011.05.023 |
| [] | 廖建平, 刘和秀, 王华忠, 等. 2011b. 快速高精度的频率空间域声波数值模拟方法研究[J]. 地球物理学进展, 26(4): 1359–1363. DOI:10.3969/j.issn.1004-2903.2011.04.029 |
| [] | 刘璐, 刘洪, 张衡, 等. 2013. 基于修正拟牛顿公式的全波形反演[J]. 地球物理学报, 56(7): 2447–2451. DOI:10.6038/cjg20130730 |
| [] | 刘玉柱, 王光银, 杨积忠, 等. 2015. 基于Born敏感核函数的VTI介质多参数全波形反演[J]. 地球物理学报, 58(4): 1305–1316. DOI:10.6038/cjg20150418 |
| [] | 罗静蕊, 吴如山, 高静怀. 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 |