2017, Vol. 60
2. 同济大学海洋与地球科学学院波现象与反演成像研究组, 上海 200092
2. Wave Phenomena and Inversion Imaging Group(WPI), Tongji University, Shanghai 200092, China
无论从偏移成像还是储层分析的角度来看，地下介质的背景速度模型对于地震勘探处理均十分重要.偏移成像从叠后走向叠前，从时间域走到深度域，对深度域速度模型质量的要求越来越高.储层分析中的AVA(振幅随角度变化)、PVA(相位随角度变化)等反演技术对地下介质的像的要求不再局限于构造特征，也包括振幅、相位等信息，从而也对速度模型有更高的要求.
随着勘探技术的发展，目前分析类的速度建模方法，如深度聚焦速度分析(Doherty and Claerbout, 1976; Faye and Jeannot, 1986)、剩余曲率速度分析(Al-Yahya，1989)，已经难以满足生产需求.基于反演的速度建模方法具有更完善的理论，可以得到更加准确的背景速度.根据使用的数据，反演方法可分为透射波反演(包括直达波、折射波、回转波等)和反射波反演，本文即利用反射波反演速度模型.此外，从正演问题的角度来看，速度反演方法主要分为射线方法和波动方程法.射线类方法可在成像域或数据域实现，如层析偏移速度分析(Tomographic Migration Velocity Analysis，T-MVA)(Stock，1992；Woodward et al., 2008)和斜率层析(Farra and Madariaga, 1988；Biondi，1992；Billette and Lambaré，1998)等.射线类层析具有应用灵活、效率高的特点，同时也受困于射线正演的阴影区、焦散问题以及反演过程的强病态性等问题.波动方程速度反演方法也可分为成像域和数据域两类，成像域波动方程层析也称为波动方程偏移速度分析(Wave-Equation Migration Velocity Analysis，WEMVA)，以成像域数据的某些特征——例如地下偏移距道集或时移成像道集的聚焦(Shen et al., 2005；Yang and Sava, 2011；等)、地表偏移距成像道集或地下角度道集同相轴的拉平等(Symes and Carazzone, 1991；Shen，2005；等)——为准则建立目标函数并通过波动方程优化，从而更新背景速度.全波形反演(Full Waveform Inversion，FWI)(Tarantola，1984；Pratt，1999；杨午阳等，2013)是一种经典的数据域波动方程反演方法，试图利用全部观测数据直接反演地下介质参数，但由于FWI具有非常强的非线性性，且实际应用中受到诸多因素制约(如实际数据的缺失，震源子波未知，波动方程无法完全描述波传播，等等)，从而FWI的反演过程难以收敛.Xu等(2012)提出了数据域的反射波全波形反演方法(Reflection Full Waveform Inversion，RFWI)，仅利用反射数据建立反演问题，大大减轻了波形反演的非线性性，Ma和Hale(2013)、付继有等(2015)把RFWI进一步简化，利用反射数据的走时差作为目标函数反演背景模型，进一步降低了反演的非线性性.波动方程类的速度反演方法在理论上更加完善，但存在计算效率低下、应用不灵活的缺点.高斯波包(Ralston，1983)是一种描述局部波场传播的方法，其理论介于射线和波动方程之间，可视为时空域的高斯束(Červeny et al., 1982；Popov，1982).高斯波包在一定程度上克服了射线方法的诸多缺点，如射线阴影区、焦散问题等，并且计算效率与射线方法相当.与波动方法相比，高斯波包具有更高的计算效率和更好的灵活性，并且在计算高频波场时精度可媲美于波动方法.本文把高斯波包应用到反射层析中，在数据域发展出一种介于射线反射层析和波动反射层析之间的背景速度反演方法.
本文利用高斯波包模拟反射波场，通过高斯波包与观测波场的互相关函数来描述走时差，并以走时差的二范数作为目标函数.在此基础上，结合一阶Born近似推导出目标函数的梯度计算公式.同时给出利用高斯波包模拟反射数据的具体策略，以及走时残差的测量方式.最后对高斯波包反射反演方法进行数值实验，在此过程中讨论了梯度的计算策略并分析了如何保证模拟反射同相轴与观测反射同相轴来自同一反射层，以及在计算走时的过程中如何克服高斯波包横向不匹配等具体问题和解决方案.数值实验反演出了可靠的背景速度模型，证实了本文方法的可靠性.
2 理论地震波在光滑的背景介质中传播，遇到扰动介质时产生扰动波场，例如波的反射或散射等.按照对波场传播的影响方式不同，本文把速度模型分成两部分——低波数的背景成分mb和高波数的扰动成分mp，即
|
(1) |
其中⊕表示某种形式的组合.前者影响波传播的透射过程，也称为背景介质，后者影响波传播的散射(或反射)过程，也称作扰动介质.相应的波场也可分成两部分——背景波场ub和扰动波场up，即
|
(2) |
必须强调的是，背景模型和扰动模型、背景波场和扰动波场之间不存在定量的界线，(1) 和(2) 式仅仅是一种形式表达.
2.1 扰动波场模拟本文中的扰动模型定义为慢度平方差：
|
(3) |
其中v0和v分别是背景速度和扰动后的速度.令扰动模型的空间分布为Gabor函数，满足如下公式：
|
(4) |
其中
|
图 1 Gabor函数描述的扰动模型(a)，入射波波前面到达扰动模型(b)，背景波场经过扰动模型散射之后产生扰动波场(c) 图中蓝色震荡图案表示散射体，红色圆弧是背景波场的波前面，红色震荡图案是散射波场(引自Klime，2012). Fig. 1 (a) Gabor-perturbation. (b) Wave-front of incident wave is arriving the Gabor-perturbation. (c) Gaussian packet scattered by Gabor-perturbation The blue pattern is perturbed model. The red arc is wave-front of the background wave field. The red oscillation pattern is scattered wave field (Klime, 2012). |
|
(5) |
其中
任意空间分布的扰动模型可以分解成Gabor函数的线性表达：
|
(6) |
在一阶Born近似下上述扰动模型产生的扰动波场为扰动高斯波包的线性叠加：
|
(7) |
其中uigp是由扰动模型mpiG产生的扰动波场，可通过高斯波包传播算子描述.
2.2 目标函数与梯度高斯波包反射层析的目标函数为反射走时残差的二范数：
|
(8) |
其优化过程可通过局部线性化近似来逼近.目标函数(8) 在局部范围内被看作二次函数，从而目标函数的优化可通过梯度类迭代方法实现，例如最速下降法、共轭梯度法、拟牛顿法等，其中目标函数的梯度是核心问题.
实际应用中扰动模型一般利用偏移成像等方法计算，所以扰动模型mp是背景模型mb的函数，即
|
(9) |
模拟的反射数据可看作背景速度模型和扰动参数模型的函数：
|
(10) |
结合(8) 式可知目标函数作为背景速度和扰动参数两者的函数，对背景速度的梯度应为目标函数对背景速度的全导数，此时利用链式法则可导出梯度的计算形式：
|
(11) |
Ma和Hale(2013)指出反射波场的运动学变量对背景介质mb的敏感程度远大于mp，所以在计算目标函数梯度时忽略目标函数对扰动介质mp的偏导数.此时梯度可近似为
|
(12) |
附录推导出了目标函数对背景速度模型的偏导数，梯度的具体计算公式为：
|
(13) |
其中sadj是伴随源，具体表达式和物理意义请参考附录中(A23) 式.ub是从震源出发的正传背景波场，可利用高斯波包的角度积分方法模拟(李辉等，2012)；ugp是炮点背景波场入射到扰动介质上产生的扰动波场，可通过扰动高斯波包理论计算(Klime，2012；李辉和王华忠，2015)；G(x rec , －t; x, 0)*sadj是以sadj为伴随源逆时传播的背景伴随波场，可对sadj进行特征高斯波包分解后利用高斯波包实现反传(李辉等，2014)；ugp(x, －t, xrec, pi)*sadj是背景伴随波场逆时反传遇到扰动介质时散射出的扰动波场，同样通过扰动高斯波包理论计算.
在目标函数对背景速度模型的梯度表达式(13) 中，ügp和üb分别是正传的扰动波场和反传的背景波场对时间的二阶导数.具体的数值算法中求解波场对时间的二阶导数需已知多个时间片的波场，从而需要额外的内存和计算量.为避免存储多个时刻的波场，结合褶积和微分的性质，我们把(13) 式中波场的时间微分修改成伴随源对时间的微分：
|
(14) |
分析表达式(14)，目标函数的梯度可分成两部分，分别为：
|
(15) |
和
|
(16) |
公式(15) 中方括号内为检波点端反传的完整的背景伴随波场，ugp为炮点波场遇到扰动介质后产生的扰动高斯波包，两者相关即为检波点端的波路径.(16) 式中方括号内为检波点波场反传时遇到扰动介质后产生的扰动高斯波包，ub是炮端完整的背景波场，两者相关即为炮端的波路径.
3 实现策略高斯波包层析可归纳为图 2所示流程.高斯波包可灵活地描述背景波场和扰动波场的传播，同时也由此给实际算法的实现带来不便，所以这里详细讨论层析流程中“模拟反射波”和“测量走时差”的具体策略.
|
图 2 高斯波包反射层析流程图 Fig. 2 Flow chart of reflection tomography approach with Gaussian packet |
根据观测反射数据和背景速度场模拟扰动高斯波包可分成两步：首先估计出扰动介质参数，然后利用扰动参数模拟反射波.
根据偏移和反偏移可知，无论背景速度正确与否，作为运动学不变量，偏移+反偏移后的数据和原始数据的走时相等(Jin，1999；Guillaume et al., 2008)，从而无法分辨出观测数据与模拟数据的走时差.如图 3所示，对于单个高斯波包描述的观测扰动波场(图 3a)，分别在正确的背景速度(3000 m·s-1的常速模型)和错误的背景速度(2800 m·s-1的常速模型)中进行成像，得到准确位置和错误位置的高波数模型扰动——图 3b和图 3c中的Gabor-扰动介质参数.以错误位置的Gabor-扰动介质作为扰动模型，在错误的背景速度(2800 m·s-1)中模拟扰动波场，最终的模拟结果和观测数据有相同的走时特征(图 3d).幸运的是，利用叠前地震数据的冗余性可检测到背景速度错误引起的走时差异.在每一个共炮道集中选取小偏移距数据对应的高斯波包，用来对扰动模型成像，得到描述扰动模型的Gabor函数，进而模拟出共炮道集中的扰动波场.上述过程可总结为：
|
图 3 单个高斯波包描述的观测波场与模拟的高斯波包对比 (a)对观测的反射数据做特征高斯波包分解后的一个高斯波包；(b)和(c)特征高斯波包在正确和错误背景速度中分别估计的Gabor-扰动模型；(d)错误的Gabor-扰动模型(c)在错误背景模型中产生的扰动波场. Fig. 3 Comparison of simulated Gaussian packet and decomposed Gaussian packet from observed data (a) One Gaussian packet used in characteristic Gaussian packet decomposition of observed reflection data. (b) and (c) The perturbed models described with Gabor function estimated from characteristic Gaussian packet in the correct and incorrect velocity models, respectively. (d) The perturbed wave simulated from incorrect perturbed model (c) in incorrect background model. |
Ⅰ对于每一个共炮道集，从分解的高斯波包中挑选出小偏移距的高斯波包；
Ⅱ把挑选出的高斯波包投影到地下介质，通过反演成像得到描述模型扰动成分的Gabor函数；
Ⅲ计算背景波场，结合扰动模型模拟出相应的扰动高斯波包；
Ⅳ叠加扰动高斯波包得到模拟的反射波.
以一个单反射层模型为例，反射层之上的速度为常数3000 m·s-1，如图 4所示，水平反射层的深度为2000 m，观测的扰动波场为水平反射层的反射波.模型参数和观测系统参数列于表 1中，有限差分法数值求解声波方程模拟的共炮点道集作为观测数据，如图 5a所示.把观测到的共炮道集进行特征高斯波包分解，取偏移距最小的高斯波包，图 5b是图 5a中共炮道集对应的最小偏移距高斯波包.在已知的错误背景速度中，一个特征高斯波包可得到一个Gabor函数表达的扰动介质的像，从而本例中一个共炮道集可得到一个Gabor函数.所有共炮道集对扰动介质成像并叠加得到最终的高波数扰动介质，如图 6b所示.至此，利用已知的错误背景速度和成像得到的扰动介质模拟出了反射波场，图 7a展示了上述策略模拟的一个单炮道集.图 7b对比了观测和模拟的反射数据，可发现两者的走时在小偏移距处相差很小甚至相等，大偏移距部分的走时相差较大.这是由于模拟反射波时使用的扰动介质参数是利用小偏移距的观测波场成像得到，从而正演的小偏移距反射数据与观测数据的走时相似.上述扰动波场模拟过程中不需对扰动介质做Gabor分解，因为通过高斯波包成像得到的扰动介质本身即为Gabor函数.
|
图 4 单反射层模型 Fig. 4 Velocity model with a single reflective layer |
|
表 1 模型和观测系统参数 Table 1 The parameters of model and geometry |
|
图 5 (a)有限差分法求解波动方程模拟的单炮道集，作为观测数据(切除了直达波)；(b)单炮道集进行特征高斯波包分解后的最小偏移距特征高斯波包.炮点位于(1000 m，10 m) Fig. 5 (a) Common shot gather with direct wave muted that simulated using finite difference algorithm, which is used as observed data. (b) Characteristic Gaussian packet with minimum offset of one common shot gather. The coordinate of shot is (1000 m, 10 m) |
|
图 6 错误背景速度中(2500 m·s-1的常速模型)的成像结果 (a) 图 5b中高斯波包的扰动模型成像结果；(b)所有炮道集小偏移距高斯波包的扰动模型成像结果. Fig. 6 Images calculated with incorrect background velocity (With constant-velocity model of 2500 m·s-1) (a) Image of perturbed model calculated with Gaussian packet in Fig. 5b; (b) Image of perturbed model calculated with minimum-offset characteristic Gaussian packets of all common shot gathers. |
|
图 7 (a)高斯波包法模拟的共炮道集；(b) “观测”的和模拟的反射数据叠加对比 Fig. 7 (a) Common shot gather simulated with Gaussian packet; (b) Superimposing comparison of "observed" and simulated reflection data |
|
图 8 高斯波包中心点坐标的时间分量(实心圆点和十字)和插值后的走时随偏移距变化的曲线 (a)分解观测数据得到的特征高斯波包；(b)描述模拟数据的高斯波包；(c)观测数据和模拟数据的走时曲线对比，其中十字和圆点曲线分别为观测数据和模拟数据的走时；(d)走时残差随偏移距变化的曲线. Fig. 8 Time component (solid circle and cross) of Gaussian packets′ central point and curve of interpolated travel time varying with offset (a) Characteristic Gaussian packets from decomposition of observed data; (b) Gaussian packets for describing simulated reflection data; (c) Comparison of travel time curves for observed data and simulated data, where the cross-curve and circle-curve stand for observed data and simulated data, respectively; (d) Curve of travel-time residuals varying with offset. |
在利用高斯波包模拟扰动波场的过程中，观测的反射数据需要进行特征高斯波包分解，特征高斯波包中心点的时间分量即为反射波场到达检波点处的走时.同时，模拟的反射波走时可直接从高斯波包的时间方向中心点坐标得到.
仍然以图 4中的单反射层模型为例，对图 5所示共炮道集进行特征高斯波包分解，特征高斯波包的时间中心点(图 8a中的圆点)插值可得观测范围内所有坐标处的走时(图 8a中的曲线).在错误速度场中模拟的反射同相轴由传播至地表的扰动高斯波包叠加组成，这些高斯波包中心点坐标的时间分量如图 8b中十字所示，图 8b中的曲线由十字点插值得到.图 8c对比了观测数据和模拟数据的走时随偏移距变化的曲线，最终反射同相轴的走时残差可通过对两条曲线作差得到，结果如图 8d所示.同前文中的结论一致，大偏移距的数据残差比小偏移距的数据残差更大.
本文走时测量策略有两个主要优点.一为充分利用了高斯波包的优势——应用灵活，方便得到数据的走时信息.更重要的是可以保证用于测量走时差的高斯波包是同一反射层反射的信号，不会出现不同反射层的反射数据之间匹配错误的情况，即“周期跳”的表现形式之一.
4 数值实验 4.1 单层模型实验利用一个单反射层模型测试反演方法，分析反演中的每个环节，以清晰地解释反演背景速度的过程.真实模型的第一层为3000 m·s-1的常速介质中加入一个高速的高斯异常体，如图 9所示.具体的模型参数及观测系统参数如表 2所示.
|
图 9 单反射层速度模型 Fig. 9 Velocity field of single-reflector model |
|
|
表 2 单反射层模型和观测系统参数 Table 2 Parameters of single-reflector model and geometry |
初始模型是速度为3000 m·s-1的常速模型，与真实模型第一层的差异仅为缺少高斯异常体.利用初始背景速度模型，每炮数据可得到一个Gabor函数表达的扰动模型的像，所有Gabor-扰动模型叠加得到完整的扰动介质的像，如图 10b所示.在x等于7500 m附近的扰动模型存在一个凸起，因为在其上方的背景速度模型偏小.由于观测系统的限制，扰动模型在右侧x大于11000 m的范围内没有成像.
|
图 10 利用初始速度进行成像得到的扰动模型成分 Fig. 10 Perturbation model imaged with initial velocity |
相对于波动方程方法，高斯波包的优点之一是应用灵活，可选择部分Gabor-扰动模型计算扰动波场，进而计算梯度.为了分析高斯波包反演背景速度的具体过程，这里计算了一个Gabor函数描述的扰动模型对梯度的贡献.中心点坐标x分量为3600 m的Gabor-扰动模型计算的Frechet导数形态如图 11a所示，可分成炮端和检端的Frechet导数两部分.炮端的Frechet导数是炮点背景波场与一个由Gabor-扰动模型散射的伴随高斯波包互相关结果，从地下反射层向地表反传的扰动高斯波包宽度逐渐增大，与从炮点正传的格林函数相关得到类似于“Banana-doughnut”(Marquering et al., 1999)的波路径.检端的Frechet导数是正传的扰动高斯波包和伴随波场的互相关结果.
|
图 11 (a)一个扰动高斯波包计算出的波路径；(b)第一次迭代的梯度 Fig. 11 (a) Wave paths calculated by one scattered Gaussian packet; (b) Gradient of the first iteration |
从常速初始模型出发，第一次迭代的梯度如图 11b所示，从中可看出在正确的位置上存在极值.梯度中的两个翼状异常是由扰动介质凸起部分引起.经过10次迭代更新的背景速度如图 12a所示，其中显示了同真实背景速度模型类似的高速异常体.高速异常体的纵向分辨率比较低，原因之一是界面附近反演精度较低，以致于和真实的高速异常体连为一体；另外一个原因是观测系统，由于地表激发地表观测，使得波传播的方向多为纵向，缺少横向传播的波场.图 12b是10次迭代的收敛曲线，纵轴所示的误差是走时差的二范数，从中可看出反演问题的收敛性态也比较理想.利用反演的背景速度和观测数据，可对模型的扰动成分进行成像，如图 12c所示.对比图 10b和图 12c，反演的背景速度模型对应的扰动模型没有凸起部分，证明更新后的背景速度比较准确.
|
图 12 (a) 10次迭代后反演的背景速度模型；(b)迭代收敛曲线；(c)反演后的背景速度对扰动模型的成像结果 Fig. 12 (a) Inverted background velocity model after 10 iterations; (b) Iteration convergence curve; (c) Perturbed model imaged with the inverted background velocity field |
本文使用的局部Sigsbee2A模型如图 13所示，模型中包含断层、尖灭等地质构造，层位密布在整个地下空间中，且高、低速度相间分布.实验使用的模型参数和观测系统参数见表 3.以真实模型为基础，利用有限差分法求解常密度声波方程，可得观测系统中每一炮数据，将这些共炮点道集作为观测数据.图 14展示了震源位于(1000 m，10 m)和(6000 m，10 m)的两个共炮点道集，其中直达波已经切除.对“观测数据”进行特征高斯波包分解是利用高斯波包法进行扰动模型成像和背景速度反演的前提条件.由于背景速度反演中只使用部分反射层位即可，所以进行特征高斯波包分解时不需拾取所有反射数据同相轴，而仅仅分解比较清楚且同相性较好的数据.图 15给出了图 14中两个共炮道集利用特征高斯波包分解再重构后的数据，从中可看出特征高斯波包分解时只关注了部分同相轴，而非追求对所有数据的分解.
|
图 13 局部Sigsbee2A模型 Fig. 13 Part of Sigsbee2A velocity model |
|
|
表 3 局部Sigsbee2A模型和观测系统参数 Table 3 Parameters of local Sigsbee2A model and geometry |
|
图 14 有限差分法求解波动方程模拟的单炮数据，作为观测数据(切除了直达波).炮道集的震源分别位于(a)(1000 m，10 m)，(b)(6000 m，10 m) Fig. 14 The "observed" common shot gathers simulated using finite difference scheme with shot located at (a) (1000 m, 10 m), and (b) (6000 m, 10 m), respectively. |
|
图 15 观测数据进行特征高斯波包分解，然后重构得到的单炮数据，炮点坐标分别为(a)(1000 m，10 m)，(b)(6000 m，10 m) Fig. 15 Reconstructed common shot gathers with characteristic Gaussian packets. Shot site is at (a) (1000 m, 10 m) and (b) (6000 m, 10 m), respectively. |
实验中的初始速度模型是沿深度方向变化的常梯度速度场，速度变化范围是2000~ 5500 m·s-1，如图 16a所示.图 16b给出了第一次迭代时目标函数的梯度，梯度中显示出速度包含突变成分，这是因为模拟的反射波场在反射界面处存在突变.注意，背景速度模型不正确，从而成像得到的高波数模型的空间位置也不正确，进而梯度中速度出现突变成分的位置也是错误的.所以每次迭代中均需要处理错误的速度突变，对更新后的速度场或直接对梯度进行平滑可擦除高波数成分.本文利用地质构造约束平滑策略(李辉和王华忠，2015)处理梯度，平滑后的速度场如图 16c所示.
|
图 16 (a)常梯度模型作为初始速度模型；(b)初始背景速度计算的梯度场；(c)地质构造约束平滑后的梯度场 Fig. 16 (a) Initial velocity model with constant-gradient of velocity; (b) Gradient field calculated with initial background velocity model; (c) Gradient field after smoothing under geologic constraints |
本次实验进行了20次迭代，收敛曲线如图 17所示.每次迭代的平滑程度不同，初始几次迭代中平滑系数较大，目的是擦除掉更新的错误高波数成分，以保证背景速度收敛到正确的结果.随迭代的进行，若干次迭代后降低一次平滑程度，这也是迭代曲线中出现跳跃的原因.最终反演的背景速度模型展示在图 18a中，图 18b是反演的速度与初始速度模型之差，即反演的更新量.从图 18中的灰度图像和相应的色标中仅能看出更新后的速度相对于初始速度更加接近真实速度.为了进一步分析反演的背景速度模型，这里对初始速度、真实速度、反演速度模型分别抽出三道，详细对比三者随深度变化曲线.如图 19所示，其中图 19a、19b、19c分别是横向位置x=1000 m、3000 m、5000 m，每张图中实线、点线、虚线分别是真实速度、初始速度和反演的速度随深度变化情况.
|
图 17 Sigsbee2A模型反演的迭代收敛曲线 Fig. 17 Iteration convergence curve of inversion with Sigsbee2A model |
|
图 18 (a) 20次迭代后反演的背景速度模型；(b)反演模型与初始模型之差 Fig. 18 (a) Background velocity model of inversion after 20 iterations; (b) Difference between inverted velocity model and the initial velocity model |
|
图 19 初始背景速度场、反演背景速度场和真实背景速度场的抽线对比 (a) x=1000 m；(b) x=3000 m；(c) x=5000 m. Fig. 19 Comparison of initial, inverted, and real background velocity models in depth direction |
分别把初始背景速度、反演速度和真实速度应用于扰动模型成像中，所有炮道集的计算结果叠加得到图 20所示的成像剖面.对比不同速度对应的叠加剖面，初始速度的像与真实速度差异较大.首先，扰动介质参数的空间位置有误，此外，扰动模型的同相轴变粗，分辨率降低.反演的背景速度对应的像整体改善比较明显，扰动模型的空间位置和分辨率都与真实速度的成像结果十分接近，只是右下角区域存在差异.可见反演的背景速度除右下角区域外，其他位置均比较准确.这是因为右下角区域被反射波覆盖的次数和角度范围较小.图 21对比了初始背景速度、反演的背景速度、真实背景速度成像得到的炮域共成像点道集(李辉等，2017)，真实速度场与反演速度场计算的道集中同相轴均已拉平，即不同炮数据计算的扰动模型深度相同，说明两者对应的背景速度比较接近.而错误的初始背景速度对应的道集中的同相轴是倾斜的.
|
图 20 不同速度场成像得到的高波数模型剖面. (a)初始背景速度；(b)更新的背景速度；(c)真实背景速度. Fig. 20 Profiles of perturbed model calculated with (a) initial velocity model, (b) inverted velocity model, and (c) real velocity model |
|
图 21 不同速度模型成像得到的炮域共成像点道集(李辉等，2017) (a, b, c)分别为初始背景速度模型成像得到的道集，横向坐标为1000 m，3500 m，6000 m；(d, e, f)分别为反演的背景速度模型成像得到的道集，横向坐标为1000 m，3500 m，6000 m；(g, h, i)分别为真实背景速度模型成像得到的道集，横向坐标为1000 m，3500 m，6000 m. Fig. 21 Shot-index common image gathers (Li et al., 2017)imaged with varied velocity models (a, b, c) Initial velocity model, and x=1000 m, 3500 m, and 6000 m; (d, e, f) Inverted velocity model, and x=1000 m, 3500 m, 6000 m. (g, h, i) Correct velocity model, and x=1000 m, 3500 m, 6000 m. |
本文提出了一种利用反射波走时反演背景速度模型的方法.以反射波走时残差的二范数为目标函数，在一阶Born近似的基础上推导出了目标函数关于背景速度模型的梯度表达式，梯度公式中的波场均利用高斯波包作为正演工具计算得到.文中详细讨论了高斯波包模拟反射波的策略和反射波走时数据的提取技术，从而计算出目标函数的梯度用于更新背景速度模型.我们通过一个高斯异常体模型和局部Sigsbee2A模型测试了本文高斯波包反演背景速度模型的方法，收敛曲线和最终反演的速度模型均从数值实验的角度证明了方法的有效性与实用性.
本文梯度的推导过程与波动方程类似，但也与之有明显的差异.目标函数中走时差的表达不是利用整道数据的互相关函数，而是高斯波包与观测波场的相关结果，从而可自然地利用高斯波包描述目标函数的梯度，而非在计算波场时生硬地用高斯波包代替波动方程.
特征高斯波包分解时需要拾取观测数据中同相轴的时空位置和斜率.人工拾取的可靠程度较高，但需要大量的人工工作，而自动拾取可能引入没有物理意义的噪声.另外，无论人工还是计算机自动拾取，在复杂的数据中均不可能拾取出所有同相轴.幸运的是，反射层析中不须拾取所有的反射数据，仅利用部分反射同相轴即可.这也给高斯波包反射层析的实用化提供了一个关于特征高斯波包分解的策略——利用计算机自动拾取反射同相轴时严格约束拾取结果，保证拾取结果宁缺毋滥.事实上，文中Sigsbee2A模型的数值实验中即抛掉了一些数据，只利用部分同相轴反演背景速度模型.
作为描述波场传播的工具，高斯波包具有高效、灵活的特点，高斯波包反演背景速度模型的方法也有同样的优点.前文指出对观测数据进行高斯波包特征分解可以仅选择部分反射同相轴，从而只有部分数据用于反演，这也体现出了高斯波包的灵活性.高斯波包灵活性的另一点体现于走时差的测量方式：扰动高斯波包在时间方向的中心点坐标可视为模拟数据的走时，结合特征高斯波包描述的走时信息，可利用高斯波包的中心点坐标计算出走时残差.
观测数据中同相轴较多时，波动方程层析的问题之一是如何正确匹配数据.例如互相关方法测量走时残差(Van Leeuwen，2010)，观测数据与模拟数据中互相关的子波可能是来自不同反射层的反射信号，这可被认为是周期跳的一种表现形式.高斯波包反演方法中不存在类似问题，文中详细讨论了走时差测量的过程，该过程确保了相互作差的走时属于来自同一反射层的信号.
附录 目标函数对背景速度模型的偏导数
在推导目标函数对背景速度模型的导数之前，首先讨论扰动高斯波包(Klime，2012；李辉和王华忠，2015)对速度模型的Frechet导数.背景波场ub在背景速度模型mb中传播，和扰动介质mp无关；扰动波场up由扰动介质散射出，在背景速度模型中传播.一阶Born近似下，扰动波场可通过求解如下时间域波动方程得到：
|
(A1) |
其中y是地下介质中的空间坐标.扰动波场可通过求解波动方程的格林函数法得到
|
(A2) |
(A2) 式中的积分范围是地下介质的物理空间，“*”是褶积符号，褶积运算的自变量为t，下同.根据高斯波包模拟扰动波场的理论可知，当扰动介质参数的空间分布为Gabor函数时，相应的扰动波场是一个高斯波包，即
|
(A3) |
其中g(y)是Gabor函数；ugp是高斯波包.ugp除了用高斯波包计算之外，还可根据格林函数法写成类似于(A2) 式的积分形式：
|
(A4) |
背景速度模型的微小扰动δmb引发背景波场的微小扰动δub，根据一阶Born近似可写出两者之间的关系：
|
(A5) |
背景速度的微小扰动除引起背景波场扰动之外，也会对扰动高斯波包产生影响，把扰动高斯波包ugp的扰动量写成δugp，扰动后的背景速度模型、背景波场和扰动波场之间仍然满足(A3) 式，此时有
|
(A6) |
把(A5) 式代入(A6) 式，整理后可得
|
(A7) |
把(A3) 式代入(A7) 式，扰动高斯波包的扰动量与背景速度模型的扰动量的关系变成
|
(A8) |
(A8) 式可看作关于δugp的波动方程，等号右端项是波动方程的源，从而可通过格林函数法求解δugp.根据一阶Born近似忽略(A8) 式中微小扰动的两次项，得到
|
(A9) |
(A9) 式两边同时除以背景速度模型的微小扰动量δmb(x)，其中x是地下空间中任一点的坐标，结合褶积运算的代数性质，得到扰动高斯波包对背景速度模型的导数：
|
(A10) |
注意到(A10) 式中褶积运算的自变量为时间t，(A10) 式等价于
|
(A11) |
参考一阶Born近似下波场扰动量的积分表达式(A2)，(A11) 式中的积分可看作格林函数
|
(A12) |
把(A12) 式代入(A11) 式，得到
|
(A13) |
(A13) 式即为扰动高斯波包对背景速度模型的Frechet导数，可分成两项，第一项是正传波场经Gabor-扰动模型散射后产生扰动高斯波包的二阶导数与检波点格林函数的褶积，第二项是背景正传波场的二阶导数和检波点扰动高斯波包的褶积，检波点扰动高斯波包是以检波点为源的格林函数逆时传播经Gabor-扰动模型散射后产生的逆时传播的高斯波包.
扰动高斯波包与观测的扰动波场的互相关函数为
|
(A14) |
其中ugp(xrec, t, xsrc, pi)是第i个Gabor-扰动模型产生的扰动高斯波包.不考虑子波相位因素，函数f极大化时(A14) 式中的时移量等于走时残差，此时互相关函数的导数为零，即有
|
(A15) |
目标函数对背景速度模型的偏导数为
|
(A16) |
(A15) 式中定义走时残差的函数
|
(A17) |
根据函数
|
(A18) |
|
(A19) |
把(A18) 和(A19) 式代入(A17) 式中，得到走时残差对背景速度模型的导数：
|
(A20) |
上述导数代入目标函数对背景速度模型的偏导数中，得到
|
(A21) |
将扰动高斯波包对背景速度模型的偏导数的表达式(A13) 代入(A21) 式中，可得目标函数对背景速度模型偏导数的最终形式：
|
(A22) |
其中
|
(A23) |
考虑如下积分变换公式(Tarantola，1987)：
|
(A24) |
(A22) 式所示的目标函数对背景速度模型的偏导数变成
|
(A25) |
其中格林函数G(xrec, -t; x, 0) 和散射高斯波包ugp(x, -t, xrec, pi)均为反传波场，所以(A23) 式可看作伴随源，产生反传播的波场.
致谢感谢中国石化石油物探技术研究院(南京)对WPI研究组及本项研究的支持.
Al-Yahya K.
1989. Velocity analysis by iterative profile migration. Geophysics, 54(6): 718-729.
DOI:10.1190/1.1442699 |
|
Billette F, Lambaré G.
1998. Velocity macro-model estimation from seismic reflection data by stereotomography. Geophysical Journal International, 135(2): 671-690.
DOI:10.1046/j.1365-246X.1998.00632.x |
|
Biondi B.
1992. Velocity estimation by beam stack. Geophysics, 57(8): 1034-1047.
DOI:10.1190/1.1443315 |
|
Č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 |
|
Doherty S M, Claerbout J F.
1976. Structure independent velocity estimation. Geophysics, 41(5): 850-881.
DOI:10.1190/1.1440669 |
|
Farra V, Madariaga R.
1988. Non-linear reflection tomography. Geophysical Journal International, 95(1): 135-147.
DOI:10.1111/j.1365-246X.1988.tb00456.x |
|
Faye J P, Jeannot J P. 1986. Prestack migration velocities from focusing depth analysis.//56th Annual International Meeting, SEG Technical Program Expanded Abstracts. Houston, Texas:SEG, 438-440.
|
|
Fu J Y, Li Z C, Yang G Q, et al.
2015. Wave equation reflection travel-time inversion in acoustic media. Oil Geophysical Prospecting , 50(6): 1134-1140.
|
|
Guillaume P, Lambare G, Leblanc O, et al. 2008. Kinematic invariants:An efficient and flexible approach for velocity model building.//78th Annual International Meeting, SEG Technical Program Expanded Abstracts. Las Vegas, Nevada:SEG, 3687-3692.
|
|
Jin S D. 1999. Constrained non-linear velocity inversion of seismic reflection data.//69th Annual International Meeting, SEG Technical Program Expanded Abstracts. SEG, 1259-1262.
|
|
Klime L.
2012. Sensitivity of seismic waves to structure. Studia Geophysica et Geodaetica, 56(2): 483-520.
DOI:10.1007/s11200-010-0088-5 |
|
Li H, Feng B, Wang H Z.
2012. A new wave field modeling method by using Gaussian Packets. Geophysical Prospecting for Petroleum , 51(4): 327-337.
|
|
Li H, Wang H, Feng B, et al.
2014. Efficient pre-stack depth migration method using characteristic Gaussian packets. Chinese Journal of Geophysics , 57(7): 2258-2268.
DOI:10.6038/cjg20140720 |
|
Li H, Wang H Z.
2015. A scattering wave modeling method using Gaussian beam and Gaussian packet in the Gabor domain. Chinese Journal of Geophysics , 58(4): 1317-1332.
DOI:10.6038/cjg20150419 |
|
Li H, Yin J F, Wang H Z, et al.
2017. A study of scatterer imaging method based on Gaussian packet. Oil Geophysical Prospecting , 52(2): 273-282.
|
|
Ma Y, Hale D.
2013. Wave-equation reflection traveltime inversion with dynamic warping and full-waveform inversion. Geophysics, 78(6): R223-R233.
DOI:10.1190/geo2013-0004.1 |
|
Marquering H, Dahlen F A, Nolet G.
1999. Three-dimensional sensitivity kernels for finite-frequency traveltimes:THE banana-doughnut paradox. Geophysical Journal International, 137(3): 805-815.
DOI:10.1046/j.1365-246x.1999.00837.x |
|
Popov M M.
1982. A new method of computation of wave fields using Gaussian beams. Wave Motion, 4(1): 85-97.
DOI:10.1016/0165-2125(82)90016-6 |
|
Pratt R G.
1999. Seismic waveform inversion in the frequency domain, Part 1:Theory and verification in a physical scale model. Geophysics, 64(3): 888-901.
DOI:10.1190/1.1444597 |
|
Ralston J. 1983. Gaussian beams and the propagation of singularities.//Littman W ed. Studies in Partial Differential Equations. MAA Studies in Mathematics. Washington, DC:Mathematics Association of America, 206-248.
|
|
Shen P. 2005. Wave equation migration velocity analysis by differential semblance optimization[Ph. D. thesis]. Houston:Rice University.
|
|
Shen P, Symes W W, Morton S, et al. 2005. Differential semblance velocity analysis via shot profile migration.//75th Annual International Meeting, SEG Technical Program Expanded Abstracts. Houston, Texas:SEG, 2249-2252.
|
|
Stork C.
1992. Reflection tomography in the postmigrated domain. Geophysics, 57: 680-692.
DOI:10.1190/1.1443282 |
|
Symes W W, Carazzone J J.
1991. Velocity inversion by differential semblance optimization. Geophysics, 56(5): 654-663.
DOI:10.1190/1.1443082 |
|
Tarantola A.
1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49(8): 1259-1266.
DOI:10.1190/1.1441754 |
|
Van Leeuwen T. 2010. Correlation-based seismic velocity inversion[Ph. D. thesis]. Delft:Delft University of Technology.
|
|
Woodward M J, Nichols D, Zdraveva O, et al.
2008. A decade of tomography. Geophysics, 73(5): VE5-VE11.
DOI:10.1190/1.2969907 |
|
Xu S, Wang D, Chen F, et al. 2012. Full waveform inversion for reflected seismic data. 74th Annual International Meeting, EAGE.W24.
|
|
Yang T N, Sava P. 2011. Image-domain waveform inversion with the adjoint-state method. CWP Annual Report. 23-39.
|
|
Yang W Y, Wang X W, Yong X S, et al.
2013. The review of seismic full waveform inversion method. Progress in Geophysics , 28(2): 766-776.
DOI:10.6038/pg20130225 |
|
付继有, 李振春, 杨国权, 等.
2015. 声介质波动方程反射走时反演方法. 石油地球物理勘探, 50(6): 1134–1140.
|
|
李辉, 冯波, 王华忠.
2012. 波场模拟的高斯波包叠加方法. 石油物探, 51(4): 327–337.
|
|
李辉, 王华忠, 冯波, 等.
2014. 特征高斯波包叠前深度偏移方法. 地球物理学报, 57(7): 2258–2268.
DOI:10.6038/cjg20140720 |
|
李辉, 王华忠.
2015. 基于高斯束与高斯波包的Gabor框架散射波模拟方法. 地球物理学报, 58(4): 1317–1332.
DOI:10.6038/cjg20150419 |
|
李辉, 殷俊锋, 王华忠, 等.
2017. 高斯波包散射体成像方法. 石油地球物理勘探, 52(2): 273–282.
|
|
杨午阳, 王西文, 雍学善, 等.
2013. 地震全波形反演方法研究综述. 地球物理学进展, 28(2): 766–776.
DOI:10.6038/pg20130225 |
|