2019, Vol. 62
2. 山东科技大学山东省沉积成矿作用与沉积矿产重点实验室, 山东青岛 266590
2. Shandong Provincial Key Laboratory of Depositional Mineralization & Sedimentary Mineral, Shandong University of Science and Technology, Qingdao Shandong 266590, China
地震波场的正演模拟是勘探地震学的重要研究课题，其不但可以用来认识地下介质的波场传播规律，进而指导地震资料的采集、处理和解释等环节，还可以利用其共轭算子(或者逆算子)来进行地下介质成像和参数反演运算(Tarantola, 1984; Bleistein, 1987; Beylkin and Burridege, 1990; Chapman, 2004).目前，地震勘探规模越来越庞大，所面临的地下地质构造也愈加的复杂，在这种背景下发展具备高精度和高效率的波场正演方法具有重要的理论意义和应用价值.
现有的地震波场正演模拟方法主要包括两类：一类是波动方程的数值求解方法，包括有限差分法、有限元法、频谱法等(Tessmer, 2000; Virieux et al., 2011).该类方法直接对波动方程进行数值求解，可以计算直达波、折射波、一次反射波和多次波等全部的波场信息，具有很高的波场模拟精度.但是该类方法的计算效率较低，特别是在模拟高频地震波场或者模型速度较低时，为保证计算过程的稳定性而采用的小计算网格往往会导致计算量大幅度增长.因此，较低的计算效率严重地制约了波动方程数值解法在三维大规模地震波场正演时的应用.另一类是射线类方法(Červeny, 1983, 2005; 李辉等, 2012)，该类方法以波动方程的高频近似解为基础，通过射线追踪来逐步求取每个射线节点的走时和振幅，进而进行一次反射地震波场的合成.常规的射线类算法(如两点射线追踪)，计算效率较高，但是其对模型的构造细节不敏感，且存在波场焦散区、阴影区等固有缺陷，致使其波场模拟精度较低.此外，在常规的射线方法中，地下的反射界面是通过离散的样点并结合样条函数进行描述的，当地下构造复杂时，准确的界面描述也是一项很大的难题.
射线Born近似或Kirchhoff近似也是常用的地震波场正演方法，具有优于常规射线方法的波场模拟精度(Hudson and Heritage, 1981; Cohen et al., 1986; Jaramillo and Bleistein, 1999; Bleistein et al., 2001)，但是该类方法依赖常规的射线追踪来求取地下走时和振幅信息，因此仍然受到常规射线方法局限性的制约，而且在模拟多次走时波场时存在较大的实现难度.针对上述问题，Huang等(2016)提出了一种基于高斯束的地震波场Born正演方法，该方法使用高斯束取代常规的地震射线进行地下散射点走时和振幅信息的计算，在理论上克服了常规射线方法的弊端.他们在实现过程中使用最小虚值走时点来构建等时线，进而通过等时线叠加求取模拟波场.由于高斯束表征的格林函数是一系列高斯束的叠加，仅仅选择最小虚值走时点显然无法代表全部的波场信息，致使文中给出的高斯束Born正演模拟结果同波动方程有限差分正演结果存在较大的差异.此外，他们所提出的算法是逐道进行的，且在频率域中实现，因此计算效率也不存在优势.
本文提出了一种新的声波介质高斯束Born正演模拟算法.该方法将Born正演过程分解为两个基本环节：首先，我们利用高斯束的走时和振幅信息，将地下散射点处的反射率函数映射为地表稀疏束中心位置处的局部平面波；然后，我们通过逆倾斜叠加将局部平面波累加到接收点波场中，得到最终的波场正演结果.在实现过程中，我们使用了Hill(2001)和Gray(2005)所提出的最速下降法来减少循环迭代次数，在保证模拟精度的同时，显著地提高了计算效率.此外，我们在平面波合成时，给出了一种wavelet-bank算法，相对于传统的wave-packet算法(Červeny, 1983)，计算效率提高了数倍.在下文中，我们首先对相关的理论公式进行推导，接下来对具体的计算实现细节进行探讨，最后通过两个数值模型对波场模拟精度和计算效率进行验证.
1 方法原理 1.1 基于高斯束的Born正演公式在各向同性声波介质中，假设震源和接收点的空间位置矢量分别为xs和xr，那么由震源xs激发，途径地下散射点x，最后到达接收点xr的一次散射地震波场u(xr, xs, ω)可以利用Born正演公式(Tarantola, 1984; Beydoun and Mendes, 1989)来表示：
|
(1) |
式中，D为地下散射点的集合，G(x, xs, ω)为由震源下行格林函数，G(x, xr, ω)为接收点上行格林函数，F(ω)为频率域的震源函数，m(x)≈2c1(x)/c03(x)为地下散射点x处速度扰动c1(x)所引起的地震反射率，c0(x)为该点的背景速度.为了保证格林函数的正则性，并且使正演过程具备处理多次走时波场的能力，我们使用高斯束来计算式(1)中的格林函数，其表达形式为一系列具有不同出射方向的高斯束的叠加(Hill，2001)：
|
(2) |
其中，p =(px, py, pz)为中心射线在出射点x′处的初始射线参数，AGB和τGB分别为高斯束的复值振幅和走时，并且具有如下的形式：
|
(3) |
式中，下标R和I分别代表复数参量的实部和虚部.将式(2)、(3)代入式(1)，可以得到以高斯束表征的Born正演公式：
|
(4) |
其中，AR(xr, xs)和AI(xr, xs)为地下散射点处震源和接收点高斯束振幅乘积的实部和虚部，τR(xr, xs)和τI(xr, xs)为相应的高斯束走时和的实部和虚部，具有如下形式：
|
(5) |
直接利用式(4)进行波场的计算，需要在每个接收点位置计算高斯束.为提高计算效率，我们选择在稀疏的地表束中心位置计算高斯束，并且通过引入一个相移校正量exp[－iωpL·(xr－ L)]，来近似计算接收点上行格林函数(Gray and Bleistein, 2009；岳玉波等，2012)
|
(6) |
并且为了减小上述近似的误差，我们在式(4)中引入下述的单位分解公式(岳玉波，2011)
得到最终的Born正演公式：
|
(8) |
其中，
|
(9) |
式(8)的作用是将每个束中心处所合成的局部平面波累加到临近的接收点波场中，该式实际上是一个局部逆倾斜叠加过程，因此可以在频率波数域进行高效计算.式(9)的作用是通过高斯束的振幅和走时信息，将地下的反射率映射为束中心处的局部平面波.该式是本文所提出的Born正演方法的核心，其计算过程决定了最终的波场模拟精度和计算效率，因此我们接下来对其进行重点讨论.
1.2 多波至模式直接求解式(9)，需要对震源和束中心高斯束的所有组合对进行循环计算.在三维情况下，若psx, psy, pLx, pLy的采样数分别为Np，那么需要Np4次的束循环运算(我们称该实现方式为全波至模式)，计算量庞大.为了减少计算量，我们应用Hill(2001)，Gray(2005)在高斯束偏移中应用的最速下降法对式(9)的计算进行简化，具体的实现过程如下：
(1) 首先，我们通过下式将震源和束中心射线参数转化为中心点和偏移距射线参数
|
(10) |
得到以中心点和偏移距射线参数索引的高斯束组合；
(2) 接下来，对于中心点射线参数为pm0的高斯束组合，按照ph进行循环找到使地下散射点处虚值走时最小的高斯束对.假设此时的偏移距射线参数为ph0，则该高斯束对中震源高斯束的射线参数为ps=(pm0－ ph0)/2，束中心高斯束的射线参数为pL=(pm0+ph0)/2，计算(11)式中的积分，并且将计算结果累加到局部平面波波场中
|
(11) |
式中，
|
(12) |
其中QS,QL分别为震源和束中心高斯束所求取的2×2复值动力学射线追踪参数矩阵(Hill, 2001; Gray and Bleistein, 2009)，c0(x), c0(L)分别为震源和束中心处的速度.
(3) 最后，重复步骤(2)直到所有的中心射线参数处理完成，得到对应不同的接收点射线参数的局部平面波U(L, xs, pL, ω).在计算式(8)后，即可得到最终的Born正演结果.
上述实现过程可以称为多波至模式，其可以对地下的大部分波场传播路径进行处理(即使存在多次走时波场)，但是所需的束循环次数却由Np4减少为4Np2次.此外，由于高斯束所携带的走时和振幅信息是平滑变化的，求取最小虚部走时的束循环过程可以在稀疏的粗网格上高效的进行.多次波至模式具有比全波至模式高得多的计算效率，而计算精度却损失不大，我们将在模型试算中对此进行验证.
1.3 Wavelet-bank算法在将束循环过程进行简化后，我们接下来探讨如何高效的实现式(11)所描述的局部平面波的合成过程.由于式(11)的计算位于整个正演模拟过程的核心位置，其求解过程在很大程度上决定了整个高斯束Born正演的精度和效率.接下来，我们首先对现有的频率域算法和wave-packet算法进行简要的介绍，然后给出一种以wavelet-bank方式实现的时间域高效算法.为了公式表达的简洁，我们将分别用AR, AI来代表式(11)中振幅项的实部和虚部，用τR, τI来代表实部和虚部走时.
(1) 频率域算法
该算法直接在频率域进行计算，然后利用反傅氏变换求取时间域的局部平面波.由于没有使用任何近似，因此其计算精度是最高的.然而，由于我们需要对每个频率成分进行一次式(11)的循环求解，因此计算量巨大.
(2) Wave-packet算法
由于频率域算法效率很低，Červeny (1983)提出了一种wave-packet算法，该方法将频率域的振荡积分转换为了时间域褶积的形式，我们在此直接给出式(11)应用wave-packet算法后的表达公式
|
(13) |
其中，
(3) Wavelet-bank算法
我们在此提出一种基于wavelet-bank方式的时间域高效算法.首先，将式(11)通过反傅氏变换转化为时间域：
|
(14) |
其中，
|
(15) |
为应用虚部走时对子波频谱进行衰减后的时间域震源子波，f H(t, τI)为其希尔伯特变换.式(14)只需要将反射率m(x)映射到时间域地震道对应τR的唯一时间样点处，但是每完成一个点的映射我们就需要一次褶积运算，计算效率很低.然而，我们发现式(15)中与虚部走时相关的振幅衰减项只会影响子波的振幅，不会改变子波的相位，因此，我们可以根据所允许的最大虚部走时τIMAX=5/ωr，创建一个规则采样的虚部走时序列：
|
(16) |
其中，N为离散样点数，ΔτI=τIMAX/N为采样间隔.然后计算出对应每个序列点τI(n)的震源子波：
|
(17) |
这样，我们便可以利用线性插值将式(15)表示为：
|
(18) |
其中，τI(n)=int(τI/ΔτI)×ΔτI，int为取整函数.此时，式(14)便可以表示为：
|
(19) |
其中，Rn(t)，In(t)为：
|
(20) |
|
(21) |
W为线性加权系数，根据虚部走时的大小，其取值如下：
|
(22) |
Wavelet-bank算法的实现过程可以大致概括为：首先，对于地下每一个散射点，根据其对应的实部和虚部走时计算式(20)，从而将反射率加权映射到对应的Rn(t)，In(t)中；然后，在所有散射点计算完成后，计算式(19)即可得到最终模拟到的局部平面波u(L, xs, pL, t).随着虚部走时的增大，式(15)中的震源子波是单调且平滑衰减的，因此我们不需要选择太大的虚部走时离散样点数N.在实际计算中，我们发现N=10即可完全满足精度的要求.
应用wavelet-bank算法，我们只需对地下每个散射点进行两次数据累加运算，此后的褶积运算可以在频率域内高效进行.因此，wavelet-bank算法的计算效率要远高于频率域算法和wave-packet算法.在接下来的模型试算中，我们同样对此进行测试和验证.
2 模型试算在本节中，我们将通过两个模型数据进行高斯束Born正演的测试.首先使用一个层状介质模型来对比wavelet-bank算法和现有方法的效率和精度.接下来，使用Marmousi模型测试高斯束Born正演在复杂模型中的应用效果.我们分别计算对应全波至和多波至模式的波场正演结果，并且同具备时间2阶、空间8阶精度的交错网格波动方程有限差分正演结果进行对比.
2.1 层状介质模型首先，使用一个层状模型测试wavelet-bank算法的计算效率和精度.图 1a展示了所使用的速度模型，该模型横向采样点数为1000，采样间隔为10 m，纵向采样点数为550，采样间隔为5m.图 1b为对速度场进行平滑后计算得到的反射率剖面.我们将炮点放在模型表面水平位置5.0 km处并使用主频为20 Hz的雷克子波作为震源，240个接收点均匀的分布在炮点两侧，间隔为20 m(如图 1a所示)，每个接收道道长为2 s，时间采样间隔为2 ms.我们利用多波至模式进行高斯束Born正演计算，并且在合成局部平面波时分别使用了频率域算法、wave-packet算法和wavelet-bank算法.我们将频率域计算结果作为基准，来比较另外两种算法的精度和效率.
|
图 1 层状介质模型 (a)速度场，地表黑色线段标识了接收点位置，炮点位于接收道中心；(b)反射率剖面. Fig. 1 Multi-layer model (a) The velocity field, where the black line segment on the surface marks the recording receivers and the shot is at the center of the receivers; (b) The corresponding reflectivity profile. |
图 2a、2b、2c分别为应用频率域算法、wave-packet算法和wavelet-bank算法计算得到的单炮正演记录，可以看到三者几乎没有任何的区别.为了进一步对比精度，我们将频率域算法结果d0作为精确解，然后使用函数ε=‖d－d0‖2/‖d0‖2来计算wave-packet或wavelet-bank算法结果d相对于d0的误差.在对比精度的同时，我们还比较了三种方法的计算时间，最终的时间和精度对比结果如表 1所示，可以看到wavelet-bank和wave-packet算法的计算精度类似，同频率域算法的相对误差均在1%以内，但是wavelet-bank算法的计算效率却比其他两种方法高得多，大约是wave-packet算法的7倍，是频率域算法的100倍以上.因此，wavelet-bank算法在保证计算精度的同时，极大提高了计算效率.
|
图 2 应用不同算法得到的单炮记录 (a)频率域算法；(b) Wave-packet算法；(c) Wavelet-bank算法. Fig. 2 Shot records simulated with different modeling algorithms (a) Frequency approach; (b) Wave-packet approach; (c) Wavelet-bank approach. |
|
表 1 不同算法的精度和效率对比 Table 1 The comparison of accuracy and efficiency between different modeling algorithms |
在进行波场模拟时，我们采取了高斯束偏移所使用的束宽度、射线参数采样间隔等参数选取准则(Hill, 2001; 岳玉波, 2001).同高斯束偏移类似，高斯束Born正演方法的波场模拟效果对于上述参数的选取并不敏感，具有较高的冗余度.
2.2 Marmousi模型接下来，使用Marmousi模型验证高斯束Born正演在复杂模型中的应用效果.图 3展示了Marmousi模型的速度场，模型横向采样点数为369，采样间隔为20 m，纵向采样点数为375，采样间隔为8 m，相应的反射率剖面如图 4所示.在正演时我们同样使用20 Hz的雷克子波作为震源函数，将炮点放置在横向位置为4.4 km处的地表位置(如图 3所示).共有201个接收道，道间隔为20 m，且沿炮点对称分布.可以看到，该炮对地下照明的区域恰好为模型中构造最复杂的部分.
|
图 3 Marmousi模型速度场 地表黑色线段标识了接收点位置，炮点位于接收道中心 Fig. 3 The velocity field of Marmousi model Where the black line segment on the surface marks the receivers and the shot is at the center of the receivers |
|
图 4 Marmousi模型反射率剖面 Fig. 4 The reflectivity profile of Marmousi model |
我们分别使用全波至和多波至模式并结合wavelet-bank算法进行单炮记录的计算，所得到的结果如图 5a和5b所示，可以看到两者的模拟效果基本相同.为了更好的验证高斯束Born正演的波场模拟精度，我们使用有限差分法进行了相同位置处单炮记录的计算，所得到的结果如图 5c所示.可以看到，高斯束Born正演结果同有限差分正演结果的波场整体分布特征相似，证明了本文方法在进行复杂模型正演时的有效性.然而，由于Born正演只能模拟一次散射波场，而有限差分法模拟了地下全部的波场信息，因此两者在局部细节上存在一定的差异.
|
图 5 应用不同正演方法得到的单炮记录 (a)全波至模式高斯束Born正演；(b)多波至模式高斯束Born正演；(c)有限差分法；(d)全波至模式(图a)同多波至模式(图b)结果之差 Fig. 5 The simulated shot gathers using different modeling methods (a) Full arrival Gaussian beam Born modeling; (b) Most arrival Gaussian beam Born modeling; (c) Finite-difference modeling; (d) The difference between (a) and (b) |
对三者的计算时间进行统计，全波至模式为3.89 s、多波至模式为0.67 s、有限差分法为31.29 s，多波至模式具有明显的效率优势.为对比全波至和多波至的差异，我们计算了两者的残差(图 5d)，并求得两者的相对误差为9.7%.相对于全波至模式，多波至模式结果精度略有损失，但是两者的差异主要在于最速下降法近似对波场整体背景振幅的改变，同时考虑到多波至模式的具有接近于全波至模式6倍的计算效率，因此我们认为上述精度损失完全可以接受.
为了进一步验证高斯束Born正演对于复杂模型的模拟效果，我们应用多波至模式模拟了总共135炮的地震记录，并且使用高斯束偏移对模拟记录进行了成像处理，所得到的偏移剖面如图 6所示.可以看到，最终的偏移结果准确的恢复了该模型复杂的构造形态，进一步验证了本文所提出的高斯束Born正演算法的正确性和有效性.
|
图 6 多波至模式模拟数据的高斯束偏移结果 Fig. 6 Gaussian beam migration image produced by the simulated data with most arrival mode |
本文提出了一种兼具模拟精度和计算效率的声波介质高斯束Born正演方法.该方法使用高斯束作为波场传播算子，不但有效的克服了常规射线方法所具有的阴影区、焦散区等固有缺陷，还具备了模拟多次走时波场的能力，可以实现复杂介质一阶散射地震波场的精确模拟.与此同时，本文还着重探讨了该方法的具体实现过程，给出了一套切实可行的实现方案，在保证波场模拟精度的前提下，极大的提高了计算效率.接下来，我们将进一步探讨弹性各向同性介质中地震波场的Born正演方法.此外，我们还将以本文算法为基础构建相应的共轭算子，进而开展最小二乘高斯束偏移技术的相关研究.
Beydoun W B, Mendes M. 1989. Elastic ray-Born L2-migration/inversion. Geophysical Journal International, 97(1): 151-160. DOI:10.1111/gji.1989.97.issue-1 |
Beylkin G, Burridege R. 1990. Linearized inverse scattering problems in acoustics and elasticity. Wave Motion, 12(1): 15-52. DOI:10.1016/0165-2125(90)90017-X |
Bleistein N. 1987. On the imaging of reflectors in the earth. Geophysics, 52(7): 931-942. DOI:10.1190/1.1442363 |
Bleistein N, Cohen J K, Stockwell J W Jr. 2001. Mathematics of Multidimensional Seismic Imaging, Migration, and Inversion. New York: Springer.
|
Červeny V. 1983. Synthetic body wave seismograms for laterally varying layered structures by the Gaussian beam method. Geophysical Journal International, 73(2): 389-426. |
Červeny V. 2005. Seismic Ray Theory. UK: Cambridge University Press.
|
Chapman C. 2004. Fundamentals of Seismic Wave Propagation. UK: Cambridge University Press.
|
Cohen J K, Hagin F G, Bleistein N. 1986. Three-dimensional Born inversion with an arbitrary reference. Geophysics, 51(8): 1552-1558. DOI:10.1190/1.1442205 |
Gray S H. 2005. Gaussian beam migration of common-shot records. Geophysics, 70(4): S71-S77. DOI:10.1190/1.1988186 |
Gray S H, Bleistein N. 2009. True-amplitude Gaussian-beam migration. Geophysics, 74(2): S11-S23. DOI:10.1190/1.3052116 |
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, 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 H, Feng B, Wang H Z. 2012. A new wave field modeling method by using Gaussian Packets. Geophysical Prospecting for Petroleum (in Chinese), 51(4): 327-337. |
Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49(8): 1259-1266. DOI:10.1190/1.1441754 |
Tessmer E. 2000. Seismic finite-difference modeling with spatially varying time steps. Geophysics, 65(4): 1290-1293. DOI:10.1190/1.1444820 |
Virieux J, Calandra H, Plessix R É. 2011. A review of the spectral, pseudo-spectral, finite-difference and finite-element modelling techniques for geophysical imaging. Geophysical Prospecting, 59(5): 794-813. DOI:10.1111/gpr.2011.59.issue-5 |
Yue Y B. 2011. Study on Gaussian beam migration in complex medium[Ph. D. thesis] (in Chinese). Qingdao:China University of Petroleum.
|
Yue Y B, Li Z C, Qian Z P, et al. 2012. Amplitude-preserved Gaussian beam migration under complex topographic conditions. Chinese Journal of Geophysics (in Chinese), 55(4): 1376-1383. DOI:10.6038/j.issn.0001-5733.2012.04.033 |
李辉, 冯波, 王华忠. 2012. 波场模拟的高斯波包叠加方法. 石油物探, 51(4): 327-337. DOI:10.3969/j.issn.1000-1441.2012.04.002 |
岳玉波. 2011. 复杂介质高斯束偏移成像方法研究[博士论文]. 青岛: 中国石油大学.
|
岳玉波, 李振春, 钱忠平, 等. 2012. 复杂地表条件下保幅高斯束偏移. 地球物理学报, 55(4): 1376-1383. DOI:10.6038/j.issn.0001-5733.2012.04.033 |