2021, Vol. 64
2. 西南交通大学高速铁路线路工程教育部重点实验室, 成都 610031;
3. 东方地球物理公司物探技术研究中心, 涿州 072751
2. MOE Key Laboratory of High-speed Railway Engineering, Southwest Jiaotong University, Chengdu 610031, China;
3. Research and Development Center, Bureau of Geophysical Prospecting(BGP), China National Petroleum Corporaton, Zhuozhou 072751, China
地震波场数值模拟是勘探地震学重要组成部分，可以揭示地震波的传播规律，有效指导地震资料的采集、处理和解释，同时还是多种地震反演方法的基础(Long et al., 2019; Wang and Zhan, 2020).传统地震正演方法通常忽略介质黏滞性对地震波传播的影响，但这一性质会显著影响地震波的动力学特征，包括能量衰减、频带变窄、相位畸变等(Walcott, 1970; White, 1975；Zong et al., 2020).因此研究快速、高精度黏声介质地震正演方法具有重要意义.
地震波场数值模拟方法主要分为两大类.一类是基于波动方程的数值解法，这类方法虽然计算效率不高但具有较高的计算精度，且可以计算全波场，因此众多学者针对这一类方法在黏声介质中的应用展开了研究(Arntsen et al., 1998; Operto et al., 2007; Abubakar and Habashy, 2013; 吴玉等，2017；Li et al., 2019; Wang et al., 2020).另一类方法是射线类地震波场数值模拟方法，这一类方法基于射线理论，具有较高的计算效率，但传统射线类正演方法存在计算阴影区，在焦散区无法获得准确的振幅并且难以处理多至走时问题(Jaramillo and Bleistein, 1999; Duquet et al., 2000; Fawcett, 2001; Chapman, 2004; Shekar and Tsvankin, 2014).
Born散射理论描述了模型扰动和散射波场之间的关系，射线Born正演即为Born散射理论与射线理论相结合的产物(Hudson, 1970; Hudson and Heritage, 1981).Coates和Chapman(1991)使用射线Born正演方法对弹性介质中的一次散射波进行正演，并在之后将相关正演结果与有限差分弹性波正演结果进行了对比(Coates and Charrette, 1993).Červený和Coppoli(1992)在已知格林函数的条件下利用射线Born正演对简单声波模型进行数值模拟.Gibson等(1993)使用射线Born正演方法对井间地震模型进行正演并与实际接收到的地震数据进行了对比.Thierr等(1999)将射线Born正演方法扩展到了三维介质中，并将其应用到了三维真振幅偏移成像中.Pollitz(2011)同样在三维地震波数值模拟中使用了射线Born正演方法.Moser(2012)分析了射线Born正演方法并强调了其针对特殊构造绕射波模拟的有效性.Šachl(2013)实现了三维弹性波射线Born正演方法，并分析了网格大小对模拟结果的影响.Sarajaervi和Keers(2018)利用等时线实现了三维时间域射线Born正演算法.
虽然射线Born正演相对于传统射线类地震正演方法的计算精度更高(Bleistein et al., 2001; 符立耘，2010)，但由于该方法采用射线类方法计算格林函数，仍然无法从根本上解决射线类正演方法所存在的问题.作为解决办法，可以使用基于高斯束叠加的格林函数代入Born正演公式中.由高斯束叠加获取的格林函数可以处理多路径及焦散问题，其有效性已在相关偏移成像方法中得到了验证(Hill, 1990, 2001; Gray and Bleistein, 2009; Liu and Palacharla, 2010;曹文俊等；2013；杨继东等，2016；韩建光等，2017；杨晶和黄建平，2018；Sun et al., 2020).Shekar和Tsvankin(2014)将基于高斯束叠加的格林函数和Kirchhoff散射理论相结合提出了一种针对各向异性衰减介质的地震波正演方法.Huang等(2016)以及岳玉波等(2019)将由高斯束叠加获取的格林函数融入了Born散射理论中，提出了一种针对声波介质的地震波正演方法，并取得了良好的数值模拟结果.
本文提出了一种针对黏声介质的二维高斯波束Born正演方法，在Born正演理论的框架下采用由高斯束叠加获得的格林函数，在高斯束传播算子中考虑了介质黏滞性信息，并将其融入了wavelet-bank平面波合成方法中.在接下来的章节中，首先推导出基于高斯束的Born正演公式，接着介绍正演方法的实现策略，最后通过两个数值模拟对本文提出方法的计算精度和计算效率进行验证.
1 基本原理 1.1 2D黏声介质中基于高斯束的格林函数计算格林函数可以由不同初射方向的高斯束叠加获得，在声波介质中，其表达式为(Červený et al., 1982)：
|
(1) |
其中G表示格林函数；x为散射点；x′为源点；p=(px, pz)为慢度，px和pz分别为慢度在水平方向和垂直方向上的分量；ω为角频率.高斯波束uGB(x, x′, p, ω)可以进一步表示为:
|
(2) |
其中s、n分别为射线中心坐标系中的坐标；p(s)、q(s)为动力学射线参数，可以通过求解动力学射线方程组获得；τ(s)为沿着射线路径的走时.
令：
|
(3) |
式(1)中格林函数可以表示为：
|
(4) |
其中：
|
(5) |
AGB(x, x′)和τGB(x, x′)分别为复振幅和复走时，下标R和I表示复数的实部和虚部.
黏声介质有多种模型，本文采用的是常Q模型，此时介质黏滞性强弱由与频率无关的品质因子Q决定.地震波在黏声介质中的传播可以看作为以一个复速度在声波介质中传播，当满足条件Q≫1, 此时的复速度V(x)可以表示为(Aki and Richards, 2002)：
|
(6) |
其中ωr为参考频率.
此时波束的复走时τGBc(x, x′)可以表示为：
|
(7) |
其中：
|
(8) |
使用τGBc(x, x′)替换掉式(4)中的τGB(x, x′)即可获得黏声介质中基于高斯束叠加的格林函数.
2D黏声介质中基于高斯束的Born正演公式
Born散射理论可以通过省略高阶项实现针对一次散射波的正演，此时的Born正演公式为(Moser，2012)：
|
(9) |
其中xs和xr分别表示源点和接收道的位置；F(ω)为震源函数；x为散射点；D为散射点的集合；c0(x)为背景速度；c1(x)为扰动速度；G(x, xs, ω)和G(x, xr, ω)分别为源点和接收点的格林函数，它们的表达式为：
|
(10) |
将式(10)代入式(9)可得黏声介质中基于高斯束的Born正演公式：
|
(11) |
其中
|
(12) |
其中pLx和pLz为中心道初始慢度在水平和垂直方向上的参数.
结合公式(13)：
|
(13) |
式(11)可以进一步转化为：
|
(14) |
其中ΔL为窗间隔；ωr为参考频率；w0为初始波束宽度；U(L, xs, psx, pLx, ω)的表达式为：
|
(15) |
式(15)的含义是将反射率信息映射为局部平面波，接着通过式(14)将局部平面波转化为窗内接收道的地震记录.式(15)为整个计算的核心也决定了正演算法的计算效率，下一小节将着重讨论式(15)的计算方法.
1.2 基于wavelet-bank的时间域实现策略直接按照式(15)在频率域中对U(L, xs, psx, pLx, ω)进行计算效率并不高，文中将给出高效的时间域计算方法.将式(5)和式(7)代入到式(15)中可得：
|
(16) |
对U(L, xs, psx, pLx, ω)进行傅里叶反变换可将其转换到时间域：
|
(17) |
其中*位褶积运算符号，f(t, τI, τQ)的表达式为：
|
(18) |
fH(t, τI, τQ)为f(t, τI, τQ)的希尔伯特变换，直接按照式(17)需要大量的褶积运算计算效率并不高，文中采取wavelet-bank方法对式(17)进行计算.
首先针对τI和τQ建立规则的序列：
|
(19) |
其中J和K分别为τI和τQ的离散点数量，根据高斯波束偏移中类似用法经验，J和K的取值当达到10时即可满足后续计算的精度要求；τIMAX为最大虚部走时，其可以通过τIMAX=5/ωr来计算；τQMAX为计算中心射线时记录的τQ最大值.在已知离散点τIj和τQk后，建立子波库：
|
(20) |
对
|
(21) |
随后通过双线性插值以及高斯束的振幅和走时将散射点的反射率信息映射到Rjk(t)和Ijk(t)中：
|
(22) |
最后通过将Rjk(t)、Ijk(t)同滤波后的子波f(t, τIj, τQk)进行褶积获得所需的局部平面波信息：
|
(23) |
整个正演算法的实现流程如图 1所示.
|
图 1 黏声介质高斯束Born正演实现流程 Fig. 1 Flow chart of the Born forward modeling for visco-acoustic media using Gaussian beam |
在本节中，将使用一个水平层状黏声模型和一个复杂黏声模型对所提出新正演方法的计算能力进行验证.
2.1 水平层状黏声模型计算图 2为水平层状黏声模型的的速度场分布和反射率信息.模型横向和纵向的网格点数均为301；横向和纵向网格间距分别为20 m和10 m.从上到下各层的声波速度分别为2000 m·s-1、2500 m·s-1、3000 m·s-1；Q值分别为60、50以及40.炮点位于模型顶部中间位置，接收道以炮点为中心向两侧均匀分布201道，道间距为20 m，正演使用的子波为雷克子波，主频为20 Hz，时间采样间隔为4 ms，时间采样点数为750.
|
图 2 水平层状黏声模型 (a)声波速度分布；(b)反射率分布. Fig. 2 Layered visco-acoustic model (a) Acoustic velocity distribution; (b) Reflectivity distribution. |
图 3为水平层状黏声模型采用不同方法获得的正演结果，图 3a为声波高斯束Born正演结果，图 3b为黏声高斯束Born正演结果.图 4进一步展示了两种方法在零偏移距处的波形对比.通过对比可以发现：没有考虑介质黏滞性的声波正演结果中与地层对应的同相轴能量更强，黏声介质高斯束Born正演方法可以很好的将介质黏滞性对地震波场的影响反映到正演结果中.在测试中，由于考虑介质黏滞性的相关计算更加耗时，黏声正演结果总共耗时2.2 s，比声波正演多消耗了约50%的时间.
|
图 3 水平层状黏声模型不同方法正演结果 (a)声波高斯束Born正演结果；(b)黏声高斯束Born正演结果. Fig. 3 Modeling results with different modeling methods for layered visco-acoustic model (a) Result of acoustic Gaussian beam Born modeling; (b) Result of visco-acoustic Gaussian beam Born modeling. |
|
图 4 零偏移距不同正演方法波形对比 虚线为图 3a中零偏移距波形，实线为图 3b中零偏移距波形. Fig. 4 Waveform comparison of zero-offset traces between different modeling methods Dashed line marks the waveform in Fig. 3a, solid line is the waveform in Fig. 3b. |
图 5为复杂黏声模型的的速度场分布和反射率分布.模型横向网格点数为250；纵向网格点数为750；横向和纵向网格间距分别为10 m和5 m.模型为常Q模型，Q的取值为50.炮点位于横向1000 m处，接收道以炮点为中心向两侧均匀分布50道，道间距为15 m，正演使用的子波为雷克子波，主频为15 Hz，时间采样间隔为4 ms，时间采样点数为750.
|
图 5 复杂黏声模型 (a)声波速度分布；(b)反射率分布. Fig. 5 Complex visco-acoustic model (a) Acoustic velocity distribution; (b) Reflectivity distribution. |
图 6为复杂黏声模型采用不同方法获得的正演结果，图 6a为黏声有限差分方法正演结果，图 6b为黏声高斯束Born正演结果.通过对比可以发现：有限差分方法可以模拟直达波等全波场，黏声高斯束Born正演方法只针对一次散射波进行模拟.两种方法得到的一次散射波场结果几乎相同.就计算时间而言，黏声有限差分正演方法共耗时101.4 s，黏声高斯束Born正演方法共耗时1.3 s，可以看出本文提出的正演方法在计算效率上具有很大的优势.
|
图 6 复杂黏声模型不同方法正演结果 (a)黏声有限差分正演结果；(b)黏声高斯束Born正演结果. Fig. 6 Forward modeling results with different methods for layered visco-acoustic model (a) Result of visco-acoustic finite-difference modeling; (b) Result of visco-acoustic Gaussian beam Born modeling. |
本文提出了一种针对二维黏声介质一次散射波场模拟的高斯束Born正演方法.该方法通过叠加不同初射方向的高斯束获取格林函数，克服了传统射线类方法计算格林函数的缺点，保证了格林函数以及波场模拟的精度.同时基于黏声高斯束所包含的信息，新的正演方法采用wavelet-bank方法快速有效的将反射率信息映射成为了窗内的局部平面波，保证了正演算法的计算效率.通过两个数值模型测试表明，本文提出的正演方法具有良好的计算精度，并且计算效率相较于黏声介质有限差分地震波正演要高很多.
Abubakar A, Habashy T M. 2013. Three-dimensional visco-acoustic modeling using a renormalized integral equation iterative solver. Journal of Computational Physics, 249: 1-12. DOI:10.1016/j.jcp.2013.04.008 |
Aki K, Richards P G. 2002. Quantitative Seismology. 2nd ed. Sausalito: University Science Books.
|
Arntsen B, Nebel A G, Amundsen L. 1998. Viscoacoustic finite-difference modeling in the frequency domain. Journal of Seismic Exploration, 7(1): 45-64. |
Bleistein N, Cohen J K, Stockwell J W Jr, et al. 2001. Mathematics of Multidimensional Seismic Imaging, Migration, and Inversion. .
|
Cao W J, Wang H Z, Li Z C, et al. 2013. Application of complex topographic Gaussian beam amplitude-preserved pre-stack depth migration in piedmont zone. Progress in Geophysics (in Chinese), 28(6): 3086-3091. DOI:10.6038/pg20130633 |
Červený V, Coppoli A D M. 1992. Ray Born synthetic seismograms for complex structures containing scatterers. Journal of Seismic Exploration, 1: 191-206. |
Červený V, Popov M M, Pšen? ík, I. 1982. Computation of wave fields in inhomogeneous media-Gaussian beam approach. Geophysical Journal International, 70(1): 109-128. DOI:10.1111/j.1365-246X.1982.tb06394.x |
Chapman C. 2004. Fundamentals of Seismic Wave Propagation. Cambridge: Cambridge University Press, 2004.
|
Coates R T, Chapman C H. 1991. Generalized Born scattering of elastic waves in 3-D media. Geophysical Journal International, 107(2): 231-263. |
Coates R T, Charrette E E. 1993. A comparison of single scattering and finite difference synthetic seismograms in realizations of 2-D elastic random media. Geophysical Journal International, 113(2): 463-482. DOI:10.1111/j.1365-246X.1993.tb00900.x |
Duquet B, Marfurt K J, Dellinger J A. 2000. Kirchhoff modeling, inversion for reflectivity, and subsurface illumination. Geophysics, 65(4): 1195-1209. DOI:10.1190/1.1444812 |
Fawcett J A. 2001. Modeling of high-frequency scattering from objects using a hybrid Kirchhoff/diffraction approach. The Journal of the Acoustical Society of America, 109(4): 1312-1319. DOI:10.1121/1.1350400 |
Fu L Y. 2010. Born-series dispersion equations and Born-Kirchhoff propagators. Chinese Journal of Geophysics (in Chinese), 53(2): 370-384. DOI:10.3969/j.issn.0001-5733.2010.02.015 |
Gibson R L Jr, Toks uz M N, Batini F. 1993. Ray-Born modelling of fracture-zone reflections in the Larderello geothermal field. Geophysical Journal International, 114(1): 81-90. DOI:10.1111/j.1365-246X.1993.tb01468.x |
Gray S H, Bleistein N. 2009. True-amplitude Gaussian-beam migration. Geophysics, 74(2): S11-S23. DOI:10.1190/1.3052116 |
Han J G, Wang Y, Chen C X, et al. 2017. Multiwave Gaussian beam prestack depth migration for 2D anisotropic media. Progress in Geophysics (in Chinese), 32(3): 1130-1139. DOI:10.6038/pg20170324 |
Hill N R. 1990. Gaussian beam migration. Geophysics, 55(11): 1416-1428. DOI:10.1190/1.1442788 |
Hill N R. 2001. Prestack Gaussian-beam depth migration. Geophysics, 66(4): 1240-1250. DOI:10.1190/1.1487071 |
Huang X G, Sun H, Sun J G. 2016. Born modeling for heterogeneous media using the Gaussian beam summation based Green's function. Journal of Applied Geophysics, 131: 191-201. DOI:10.1016/j.jappgeo.2016.06.004 |
Hudson J A. 1970. Scattered waves in the coda of P. Journal of Geophysics, 43(1): 359-374. |
Hudson J A, Heritage J R. 1981. The use of the Born approximation in seismic scattering problems. Geophysical Journal International, 66(1): 221-240. DOI:10.1111/j.1365-246X.1981.tb05954.x |
Jaramillo H H, Bleistein N. 1999. The link of Kirchhoff migration and demigration to Kirchhoff and Born modeling. Geophysics, 64(6): 1793-1805. DOI:10.1190/1.1444685 |
Li Q Q, Fu L Y, Zhou H, et al. 2019. Effective Q-compensated reverse time migration using new decoupled fractional Laplacian viscoacoustic wave equation. Geophysics, 84(2): S57-S69. DOI:10.1190/geo2017-0748.1 |
Liu J, Palacharla G. 2011. Multiarrival Kirchhoff beam migration. Geophysics, 76(5): WB109-WB118. DOI:10.1190/geo2010-0403.1 |
Long Y, Lin J, Li B, et al. 2019. Fast-AIC method for automatic first arrivals picking of microseismic event with multitrace energy stacking envelope summation. IEEE Geoscience and Remote Sensing Letters, 17(10): 1832-1836. DOI:10.1109/LGRS.2019.2952571 |
Moser T J. 2012. Review of ray-Born forward modeling for migration and diffraction analysis. Studia Geophysica et Geodaetica, 56(2): 411-432. DOI:10.1007/s11200-011-9046-0 |
Operto S, Virieux J, Amestoy P, et al. 2007. 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver:A feasibility study. Geophysics, 72(5): SM195-SM211. DOI:10.1190/1.2759835 |
Pollitz F. 2011. High-frequency Born synthetic seismograms based on coupled normal modes. Geophysical Journal International, 187(3): 1420-1442. DOI:10.1111/j.1365-246X.2011.05188.x |
Šachl L. 2013. Validation of 3D synthetic seismograms based on the ray-Born approximation. Studia Geophysica et Geodaetica, 57(1): 84-102. DOI:10.1007/s11200-012-0722-5 |
Sarajaervi M, Keers H. 2018. Computation of ray-Born seismograms using isochrons. Geophysics, 83(5): T245-T256. DOI:10.1190/geo2017-0669.1 |
Shekar B, Tsvankin I. 2014. Kirchhoff modeling for attenuative anisotropic media using Gaussian beams. Geophysics, 79(5): WB51-WB61. DOI:10.1190/geo2013-0432.1 |
Sun H, Gao C, Zhang Z H, et al. 2020. High-resolution anisotropic prestack Kirchhoff dynamic focused beam migration. IEEE Sensors Journal, 20(20): 11753-11760. DOI:10.1109/JSEN.2019.2933200 |
Thierry P, Lambaré G, Podvin P, et al. 1999. 3-D preserved amplitude prestack depth migration on a workstation. Geophysics, 64(1): 222-229. DOI:10.1190/1.1444518 |
Walcott R I. 1970. Flexural rigidity, thickness, and viscosity of the lithosphere. Journal of Geophysical Research, 75(20): 3941-3954. DOI:10.1029/JB075i020p03941 |
Wang N, Zhu T Y, Zhou H, et al. 2020. Fractional Laplacians viscoacoustic wavefield modeling with k-space based time-stepping error compensating scheme. Geophysics, 85(1): T1-T13. DOI:10.1190/geo2019-0151.1 |
Wang X, Zhan Z W. 2020. Moving from 1-D to 3-D velocity model:automated waveform-based earthquake moment tensor inversion in the Los Angeles region. Geophysical Journal International, 220(1): 218-234. DOI:10.1093/gji/ggz435 |
White J E. 1975. Computed seismic speeds and attenuation in rocks with partial gas saturation. Geophysics, 40(2): 224-232. DOI:10.1190/1.1440520 |
Wu Y, Fu L Y, Chen G X. 2017. Forward modeling and reverse time migration of viscoacoustic media using decoupled fractional Laplacians. Chinese Journal of Geophysics (in Chinese), 60(4): 1527-1537. DOI:10.6038/cjg20170425 |
Yang J, Huang J P. 2018. Time domain Gaussian beam migration method based Hanning window filtering. Progress in Geophysics (in Chinese), 33(3): 1161-1166. DOI:10.6038/pg2018BB0183 |
Yang J D, Huang J P, Li Z C, et al. 2016. Gaussian beam migration based on matching pursuit sparse decomposition. Progress in Geophysics (in Chinese), 31(3): 1237-1245. DOI:10.6038/pg20160342 |
Yue Y B, Qian Z P, Zhang X D, et al. 2019. Gaussian beam based Born modeling method for single-scattering waves in acoustic medium. Chinese Journal of Geophysics (in Chinese), 62(2): 648-656. DOI:10.6038/cjg2019M0367 |
Zong J J, Stewart R R, Dyaur N. 2020. Attenuation of rock salt:Ultrasonic lab analysis of Gulf of Mexico coastal samples. Journal of Geophysical Research:Solid Earth, 125(7): e2019JB019025. |
曹文俊, 王华忠, 李振春, 等. 2013. 复杂地表高斯束保幅叠前深度偏移在山前带地区的应用. 地球物理学进展, 28(6): 3086-3091. DOI:10.6038/pg20130633 |
符力耘. 2010. Born序列频散方程和Born-Kirchhoff传播算子. 地球物理学报, 53(2): 370-384. DOI:10.3969/j.issn.0001-5733.2010.02.015 |
韩建光, 王赟, 陈传绪, 等. 2017. 二维各向异性多波高斯束叠前深度偏移. 地球物理学进展, 32(3): 1130-1139. DOI:10.6038/pg20170324 |
吴玉, 符力耘, 陈高祥. 2017. 基于分数阶拉普拉斯算子解耦的黏声介质地震正演模拟与逆时偏移. 地球物理学报, 60(4): 1527-1537. DOI:10.6038/cjg20170425 |
杨继东, 黄建平, 李振春, 等. 2016. 基于匹配追踪稀疏分解的高斯束成像方法. 地球物理学进展, 31(3): 1237-1245. DOI:10.6038/pg20160342 |
杨晶, 黄建平. 2018. 基于汉宁窗函数滤波时间域高斯束成像方法. 地球物理学进展, 33(3): 1161-1166. DOI:10.6038/pg2018BB0183 |
岳玉波, 钱忠平, 张旭东, 等. 2019. 声波介质一次散射波场高斯束Born正演. 地球物理学报, 62(2): 648-656. DOI:10.6038/cjg2019M0367 |