二维时间空间域和频率空间域声波全波形速度反演方法的对比研究

1、地球物理学进展 2017,05(32),2029-2034 d01:10.6038/pg20170523二维时间空间域和频率空间域声波全波 形速度反演方法的对比研究廖建平 刘和秀 戴世鑫 赵延林andrew hursthouse湖南科技大学页岩气资源利用湖南省重点实验室煤 炭资源与安全开采国家重点实验室（中国矿业大学） 湖南科技大学地质系导出/参考文献关注分享收藏打印 摘要：木文将二维时间空间域和频率空间域声波全波形速度反演方法分别应用到 marmousi模型，进行数值试验两种方法均采用相同的观测系统和其他的参数， 理论模型的数值试验结果证实了：使用较多的计算集群的cpu进行二维频率空间 域直

2、接法声波全波形反演时，其加速有限（正演数值模拟的计算量主要用于稀 疏矩阵的lu分解，炮点计算波场时为线性关系）.二维时间空间域声波全波形反 演计算时更灵活，多炮同时计算时，可以多倍提高其计算效率;二维声波全波形 速度反演吋，直接法求解频率空间域的计算速度远快于吋间空间域，所需要的 计算机内存也比时间空间域少二维声波全波形速度反演时，相比较于时间空间 域的方法，频率空间域直接法声波全波形反演具有计算速度快和节省计算机内 存需求的优势.关键词:吋间-空间域;频率-空间域;二维声波波动方程;全波形速度反演;作者简介：廖建平，男，1976年生，博士，副教授，主要从事复杂介质屮地震 波传播及数值模拟、走

3、时反演、二维和三维频率-空间域、时间-空间域声波和弹 性波全波形反演方法研究.（e-ma订：bobl jp163. com）收稿日期:2017-02-22 基金：国家自然科学基金（41274126, 51274097）research on comparisons of 2d acoustic wave full waveform velocity inversion in time-space domain and frequency-space domainliao jian-ping liu he-xiu dai shi-xin zhaoyantin rew hursthousehuna

4、n provincial key laboratory of shale gasresource utilization, hunan university ofscience and technology;abstract：the 2 i) timc-spacc domain and frequcncy-space domain acoustic wave full waveform velocity inversion method ( fwi) are applied to marmousi model respectively. these two methods both use a

5、 same geometry and other parameters in the numerical experiments. it is confirmed that when using more computer nodes for direct fwt in frequency-space domain, its acceleration rate is limited ( the forweird modeling calculation is mainly used for the lu decomposition of the impedance matrix, shot p

6、oint calculating the wavefield is linear relationship) 2 d time-space domain fwi is more flexible, when more shots are calculated at the same time, it can be many times to improve calculation efficiency. the calculation speed of 2 d frequenc-space domei in direct fwi is much fas ter tha n that of th

7、e 2 d time-space domain fwi, the memory requirements are much smaller than that of the 2 d time-space domain fwi. compared with the method of time-space domain fwi, frequencyspace domain method of direct fwi has the advantages of high compuling speed and saving computer memory requirements.keyword：t

8、ime-space domain; frequency-space domain; 2d acoustic wave equation; ful1 waveform velocity inversion;received： 2017-02-220引言随着油气勘探、开发程度的不断深入，油气工业界对高精度成像和储层描述、地 震解释提出了更高的要求，急需地震勘探方法技术，尤其是精确的反演成像方 法技术的进步地震波反演成像包含两个层次，即反肘系数估计(仅仅定位反射 系数出现的位置的叠前偏移成像)和全波形反演(full waveform inversion, fwt)估计速度、密度、各向异性参数和吸收衰

9、减等弹性参数.tarantola (1984)的工作是地震反演的奠基之作.成功地把偏移理论归纳为一 个数学上的局部最优化问题，即把观测数据和模拟数据z间的残差指向最小的 最小二乘问题通过从震源激发的入射波场和反向传播的残差波场之间的互相关 而建立残差函数的梯度，开启了地震波形反演的大门当存在强烈的非线性时， 反演可能收敛到一个局部极小值(gauthier et al., 1986)而失败，因此，初 始模型必须位于全局极小值的附近.20世纪80年代，地震波反演处于蓬勃发展时期，多位学者(tarantola,1984;mora, 1987, 1989)对fwi理论做了非常扎实的理论研究但受当时野外

10、地 震观测技术和高性能计算机技术的局限，无法将提岀的理论方法用于工业界.近几年，随着高性能计算机技术的飞速进步以及全方位、宽频带、高密度的地震 数据野外采集技术的飞速发展，以高精度、多参数建模为优点的反演估计全波数 或宽波数带的弹性参数为口标的fwi方法再次成为当前勘探地震学的研究热点. 理论上，fw1已被证明是一种高精度速度建模的有效手段(tarantola, 1984;tarantola et al., 1986;pratt, 1990, 1999;pratt et al., 1998;virieux and operto, 2009;virieux et al., 2017).但是，在数

11、学上它是一个高度病态 的非线性问题(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 ff拓了在频率-空间域使用直接法进行二维全波形正反演方法研究 (prat

12、t, 1990;pratt et al., 1998;pratt, 1999;sirguc and pratt, 2004)chari-hyun等(1996)利用9点法进行二维声波正演，提高了计算精 度;shin仃998)推广到25点声波正演，儿何数量级地增加了内存需求，降低了 数值模拟的计算速度；tekl等(1998)实现了二维弹性波9点法正演;min等 (2000)实现了弹性波25点法正演;min等(2003, 2004)讨论了自由表面等边 界条件的频率空间域有限元正演jiustcdt等(2004)实现了交错网格和pml的 二维频率空间域正演.综合比较计算精度、内存需求和计算速度，大部分研

13、究人 员釆用9点法进行二维频率-空间域正反演.operto等(2006)、shin和min (2006)、bednar等(2007)将二维频率-空间域全波形反演发展到黏声波、黏弹性波和各向异性介质，广泛应用于井间地震、vsp、海洋地震、地面地震资料 中.刘璐等(2013)基于修正拟牛顿公式的全波形反演;魏哲枫等(2014)基于非规 则网格声波正演的时间域全波形反演;刘玉柱等(2015)基于born敏感核函数 的vti介质多参数全波形反演;张文生等(2015)频率多尺度全波形速度反演； 王本锋等(2016)基于t-matrix的非线性参数估计方法;崔永福等(2016)全 波形反演在缝洞型储层速度

14、建模中的应用；罗静蕊等(2016)地震包络反演对局 部极小值的抑制特性研究;郭雪豹等(2016)基于频域衰减的时域全波形反演； 陈生昌和陈国新(2016)时间二阶积分波场的全波形反演；胡勇等(2017)基于 精确震源函数的解调包络多尺度全波形反演;桂生等(2017)简化混合域全波形 反演多gpu加速策略.波形反演既可以在时间域实现(gauthier等，1986;mora, 1987) 乂可以在频率 域实现(pratt and worthington, 1990;pratt et al. , 1998) 在反演策略上, 由于吋间域fwt可以灵活地对数据进行必要的预处理，以及灵活选取所需的特 征波

15、等特点而受到关注.但时间域同时反演所有频率成分，增大了反问题的非线 性，因此，bunks等(1995)提出了时间域的多尺度反演策略，通过对资料的处 理将问题分解为不同的空间尺度，增加了数值计算时的稳定性，降低了非线性 程度.由于吋间域的数据量大，需要大量的计算机内存进行存储，而当前的计算机集 群往往不能满足要求，特别是在三维情况下，这一问题显得尤为突出频率域反 演使用几个离散频率的傅里叶级数取代庞大的时间采样点信息，大量节省了计 算机内存资源其次，频率域迭代算法从低频成分开始反演，然后到高频的反演, 充分利用低频分量的线性性，而降低了全波形反演的非线性，相比较于时间域 反演而言，频率域更利于实

16、现这一算法尤其是对衰减系数等物理参数的反演， 频率域反演具有时间域反演无法比拟的优越性多尺度策略能够减轻全波形反演 的非线性频率域能够容易地实现基于任意带宽和采样间隔的频率组的多尺度全 波形反演，而时间域对于任意观测系统都可以灵活运用时窗.木文讨论了二维频率空间域直接法声波波动方程全波形速度反演和二维时间空 间域声波波动方程全波形速度反演，并对它们进行了理论模型的数值试验，通 过对数值试验结果的分析，最后得出了相应的结论.1二维时间空间域和频率空间域声波全波形反演理论二维常密度频率空间域声波方程用矩阵形式表达为(pratt等，1998): 式中3为角频率，u（3）是傅里叶变换后的离散波场，网格

17、点的个数为1, s （3）是一个具有1x1个元素的阻抗矩阵，其元素为复数，f（3）为震源项，u （3）和f（3）是1x1的列矢量.矩阵s（3）为高度稀疏、非对称、非正定，对 它直接求逆非常困难通常采用直接法因式分解求解该稀疏矩阵我们采用 chan-hyun等（1996）的二阶有限差分格式，进行二维频率空间域直接法声波 全波形反演.假设模型参数的个数为m,根据nixl的列矢量参数集p定义模型，通常水1,我 们从参数集p中计算1x1的阻抗矩阵s（3）.假设具有n个试验观测d,则方程（1）屮的每一项都依赖于角频率3,我们通 过减少观测数据和模拟数据之间的残差函数，而不断迭代更新模型，实现反演. 认为

18、存在一个足够接近于全局解的初始模型p,给定初始模型，通过正演模拟 得到模拟波场值.图 1 真实的 marmousi 速度模型 fig. 1 the real marmousi velocity model下载原图cross(km)(llnplpdocl14(x)0* 3500-3000ij 2500120001500initial model图2平滑后的marmousi模型作为反演的初始速度模型fig. 2 the starting marmousi velocity model 下载原图o.l l2.2.(uj2)£dsoo 5- o 5 0 5图3主频为5 hz的二维时间空间域声

19、波全波形反演后得到的速度模型fig. 3 the inverted resuit of marmousi velocity using 5 hz in time-space domain 下载原图图4主频为8 hz的二维时间空间域声波全波形反演后得到的速度模型fig. 4 the inverled result of marmousi velocity using 8 hz in lime-space domain下载原图定义剩余残差§ d为初始模型响应和观测数据之间的差，即：式屮，下标i代表检波器的个数，我们寻求最小化波场残差的平方和，公式为£( p) =a (3)即误差

20、函数的12范数，其中的t为矩阵转置，*为复共觇详细的理论推导请见 参考文献(pratt et al., 1998)(udi)£dcr(cn图8 12贱的二维频率空间域声波全波形反演后得到的速度模型fig. 8 the inverted result of marmousi velocity using 12 ilz in frequency-space domain 下载原图本文使用l-bfgs算法（即限定内存的bfgs算法）来求解模型更新量，进行二维 时间空间域声波全波形反演.l-bfgs方法只需要存储有限个数的模型量，在每 次迭代中只需要输入目标函数的梯度，通过保存最近几次的曲率

21、信息来更新近 似矩阵的这种方法在实践屮是很有效的虽然l-bfgs算法是线性收敛，但是每 次迭代的开销非常小，因此，l-bfgs算法运算速度很快，而且由于每一步迭代 都能保证近似海森矩阵正定，因此，算法的稳健性很强.我们使用经典的伴随状态法计算目标函数的梯度这种方法使用炮点正向传播入 射波场与检波点逆吋传播残差波场的互相关生成梯度，避免了直接计算 frcchet 导数.2二维时间空间域和频率空间域声波全波形反演的数值试验图1是marmousi模型的真实速度模型图2是对真实速度模型图1进行平滑而得 到的模型作为初始速度模型.我们使用该平滑速度模型作为初始速度模型，分别进行二维时间空间域声波全 波形

22、反演和二维频率空间域直接法声波全波形反演数值试验.为了验证二维时间空间域和频率空间域直接法声波全波形反演的优缺点，我们 在作数值试验时，都采用相同的观测系统和参数观测系统如下所示:纵向为 279个采样点，横向为1359个采样点，纵向和横向的网格间距都为10 ni, 共 有130炮，第一炮的位置为纵向10 m,横向10 m,炮点间距为纵向0 m,横向 为10 m,检波点数为1358,检波点的初始位置为纵向10 m和横向10 m,检波 点间距为纵向0 m,横向10 m.我们进行二维时间空间域声波全波形反演（使用l-bfgs算法）数值试验时，采 用多尺度的策略，从低频到高频进行反演，将低频的反演结果

23、作为高频反演时 的初始速度模型，依次进行主频为5 hz、8 hz和12 hz的波形反演数值试验.数值试验反演得到的速度分别为图3、图4和图5.我们从反演所恢复出的速度模 型中可以看岀，随着反演的频率增大，反演所得速度模型的复杂构造细节愈发 清晰可辨，说明反演的分辨率越高.二维频率空间域直接法声波全波形反演数值试验时，我们对该初始速度模型分 别进行5 hz、8 hz和12 hz的二维频率空间域直接法声波全波形反演数值试验. 反演所得到的速度分别展示在图6、图7和图8.我们从反演所恢复出的速度模型 中可以看岀，随着反演的频率增大，反演所得的模型的复杂构造细节愈发清晰 可辨，说明反演的分辨率越高.我

24、们将二维时间空间域和频率空间域声波全波形反演所恢复的速度模型进行对 比，分析二维时间空间域和频率空间域直接法声波全波形反演的优缺点.我们使用14个计算节点仃30个进程）进行二维时间空间域声波全波形反演数 值试验，全波形反演该模型所需要的总吋间为119172 s,需要的内存为9918 mb. 使用5 hz、8 hz和12 hz三组频率，每组频率都进行25次迭代.二维频率空间域直接法声波全波形反演数值试验时，我们使用14个计算节点， 共使用130个进程，全波形反演该模型所需要的总计算时间为68510 s,需要的 内存为6470 mb.使用5 hz、8 hz和12 hz三组频率，每组频率都进行25次

25、迭 代.二维频率空间域使用更多的计算节点时，其加速有限（主要用于矩阵的lu分解, 炮点计算波场为线性关系），二维时间空间域更灵活，多炮同时计算，可以多 倍提高计算效率.使用相同的频率进行反演时，二维时间空间域声波全波形反演所得到的速度的 分辨率比二维频率空间域直接法声波全波形反演数值试验吋反演得到的速度的 分辨率更高二维时间空间域声波全波形反演，使用12 iiz反演时，所得到的速 度和真实的理论模型相比较，二者z间的差别不是很大，复杂构造的形态一致， 特别是在浅层，二者几乎完全相同在中深层处，二者之间的差异稍微大些，这 主要是由于我们所使用的频率仅仅为12 hz.数值试验的结果证明了本文反演方

26、 法的分辨率较高.3结论3. 1通过二维时间空间域和频率空间域声波全波形反演marmousi理论模型的数 值试验，我们可以得出如下的结论和认识：仃)使用更多的节点进行二维频率空间域直接法声波全波形反演时，其加速有 限(正演模拟的计算量主要用于矩阵的lu分解，炮点计算波场为线性关系)吋间空间域更灵活，多炮同吋计算，可以多倍提高计算效率.(3)二维声波全波形反演时，频率空间域的计算速度远快于时间空间域，需要 的内存也比时间空间域少.3.2总之，二维声波全波形反演时，相比较于时间空间域方法，频率空间域直 接法声波全波形反演具有计算速度快和节省计算机内存需求的优势.参考文献bednar j b, sh

27、in c, pyun s. 2007.comparison of waveform inversion, part 2:phase approachjgeophysical prospecting, 55 (4) :465-475.bian a f, yu w h, zhou h w. 2010. progress in the frequency-domain full waveform inversion methodj. progress in geophys. (in chinese) , 25 (3) :982-993, doi: 10. 3969/j. issn. 1004-290

28、3. 2010. 03. 037.bunks c, saleck f m, zaleski s, ct al. 1995. multiscale seismic waveform inversionjgeophysics, 60 (5) :1457-1473.chari-hyun j, shin c, suh j h.1996. an optimal 9-point, finitedifference, frequency-speice, 2d scalar wave extrapolatorj geophysics, 61(2) : 529-537.chen s c, chen g x. 2

29、016. full waveform inversion of the seco nd-order time integral wavefieldjchinese journal of geophysics (in chinese) , 59 (10) :3765-3776, doi:10. 6038/cjg20161021.cui y f, peng g x, wu g c, el al. 2016. appl ication of ful 1 waveform inversion velocity model-building technology for the fractured-vu

30、ggy reservoirjchinese journal of geophysics (in chinese) , 59(7) :2713-2725, doi:10. 6038/cjg20160734.gauthier 0, virieux j, tarantola a. 1986. two-dimensional nonlinear inversion of seismic waveforms:numerical resultsj. geophysics, 51 (7) : 1387-1403.gui s, liu ii, zhang y j. 2017. the acceleration

31、 strategy of a simplified hybrid domain full waveform inversion on multi-gpusj. 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 c

32、hinese), 59 (10) :3777-3787, doi:10. 6038/cjg20161022hu y, han l g, xu z, et al. 2017. demodulation envelope multi-scale full waveform inversion based on precise seismic source functionj. chinese journal of geophysics (in chinese) , 60 (3) :1088-1105, doi:10. 6038/cjg20170321.jllustedt b, operto s,

33、virieux j. 2004. mixed-grid and staggered-grid finite-difference methods for frequency-domain acoustic wave modellingtj.geophys. j. int., 157 (3) :1269-1296.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. progr

34、ess 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 domaintheoretical modelj. progress in geophys. (in chinese) , 26 (5) :1690-169

35、5, doi :10. 3969/j. issn. 1004-2903. 2011. 05. 023.liu l, liu ii, zhang ii, et al. 2013. full waveform inversion based on modified quasinewton cquationjchinese 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 f

36、ullwaveform inversion for vti media based on born sensitivity kernelsj.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 inversionj. chinese journal of geophysics (in chinese) , 59(7) :25

37、10-2518, doi:10. 6038/cjg20160716.mora p. 1987. nonlinear two-dimensional elastic inversion of multioffset seismic dataj. geophysics, 52 (9) :1211-1228.mora p. 1989. tnversi on二migrat ion+tomographyj. geophysics, 54 (12) : 1575-1586.0perto s, virieux j, dessa j x, et al. 2006. crus tai seismic imagi

38、ng from multifold ocean bottom seismometer data by frequency domain full waveform tomography:application to the eastern nankai troughj. 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:el

39、astic wave-equation methodj. geophysical prospecting, 38 (3) :311-329.pratt r g. 1999. seismic waveform inversion in the frequency domain, part 1:theory and verification in a physical scale modelj. geophysics, 64 (3) :888-901pratt r g, shin c, hick g j. 1998. gauss-newton and full newton methods in

40、frequency-space seismic waveform inversionj. geophysical journal international, 133 (2) :341-362.prat t r g, worthington m h. 1990. inverse theory applied to multi source cross-hole tomography. part l:acoustic waveequation methodj.geophysical prospecting, 38 (3) :287-310.shin c, min d j. 2006. wavef

41、orm inversion using a logarithmic wavefieldj. geophysics, 71 (3) :r31-r42.sirgue l, pratt r g. 2004. efficient waveform inversion and imaging:astrategy for selecting temporal frequenciesjgeophysics, 69 (1) : 231-248tarantola a. 1984. inversion of seismic reflection data in the acoustic approximation

42、j. geophysics, 49 (8) :1259-1266.virieux j, asnaashari a, brossier r, et al. 2017. an introduction to full waveform inversiona. /grechka v, wapermar k eds. encyclopedia of exploration geophysics m. tulsa, ok:seg, r1t-r1 -40.virieux j, operto s. 2009. an overview of full-waveform inversion in explora

43、tion geophysicsj. geophysics, 74 (6) :wcc1-wcc26.wang b f, wu r s, chen x h, et al. 2016. non-linear parameter estimation method based on t-matrixj. chinese journal of geophysics (in chinese), 59 (6) :2257-2265, doi:10. 6038/cjg20160628.wei z f, gao 11 w, zhang j f 2014. time-domain full waveform in

44、version based on an irregular-grid acoustic modeling methodjchinese 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 domainj. 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 inversionj.chinese journal of geo