2017, Vol. 60
2. 同济大学海洋与地球科学学院波现象与反演成像研究组, 上海 200092
2. Wave Phenomena and Inversion Imaging Research Group(WPI), School of Ocean & Earth Science, Tongji University, Shanghai 200092, China
油气地震勘探的最终目标是储层.实现储层描述有两条重要路线，一是偏移速度估计、保真成像及角度道集提取、AVA(Amplitude.vs.Angle，振幅随角度变化关系)反演、储层描述，这是目前方法技术相对成熟且常用的路线；另一种由叠前地震数据经FWI(Full Wave Inversion，全波形反演)直接反演宽波数带的弹性参数场，用于储层描述.上述两条路线的核心都是地震波反演成像、最终目标是建立全频带的速度场，通过精确的、精细的速度场建立达到储层描述的目的.第一条路线是石油地震勘探学家发展的比较实用化的路线，主要是以NMO(Normal Move Out，动校正)叠加速度分析、PSTM(Pre-stack Time Migration，叠前时间偏移)偏移速度分析、PSDM(Pre-stack Depth Migration，叠前深度偏移)偏移速度分析和透射波及反射波层析估计背景速度，结合偏移来估计速度场的高波数变化成分(定位高波数成分出现的空间位置).目前这依然是石油工业界的核心技术路线，某些成熟的技术已广泛应用于石油工业界.第二条路线是理论地震学家(Tarentola, 1984)发展的一条基于数学反演理论框架的路线，主要是以FWI和LSRTM(Least Square Reverse Time Migration，最小二乘逆时偏移)相结合，FWI估计背景速度，LSRTM估计速度场的高波数变化成分，这在理论上是逼近全波数带反演的较优路线.目前尚处于理论研究和探索应用阶段.本文的研究内容属于第一条路线中的反射波层析估计背景速度.
对现有主流的层析成像方法进行分类，根据执行域可以分为成像域和数据域层析成像；根据正问题可以分为基于波动理论和基于射线理论的层析成像；根据利用数据属性可以分为波形和走时层析成像；根据利用数据类型主要可以分为透射波和反射波层析成像.根据正问题及执行域，分别总结和评价基于波动理论和基于射线理论的、在成像域和数据域实现的速度层析反演方法特点，整体可以概括为：(1) 基于射线理论的层析偏移速度反演方法(成像域)发展较为成熟，在实际应用中已经实现自动化层析(层析控制点的自动拾取、RMO(Residue Move Out，剩余时差)的自动计算、利用某些相似性准则自动筛选数据等).然而，射线理论的固有问题(阴影区、多路径问题等)限定了该方法的反演精度和在复杂速度模型中的适用性，因而反演效果仍然有待改善(这也是本文考虑的方向)；(2) 波动方程偏移速度分析(WEMVA)方法以波动方程作为正演算子，克服了射线理论的高频近似假设；其提取的成像道集(扩展成像条件)可实现自动化反演，可适用复杂模型的速度估计.其反演精度理论上优于射线类层析方法.但是，相较于射线类层析，WEMVA方法的计算量巨大，更存在地下照明不均衡及泛函梯度假象等问题(冯波，2015)，影响反演的收敛速度甚至不收敛，在实际应用中较少推广，性价比不高.(3) 数据域射线层析与成像域射线层析的理论问题一样，基于高频射线理论，对初始模型的精度要求较低，但反演精度不高，只能反演大尺度的速度异常，对低速异常体不敏感.为改善层析矩阵的稀疏性，胖射线层析(Vasco et al., 1995)将射线加宽，考虑到波传播是有一定宽度的.但这仅仅是一种感性认识，并没有严格的理论基础，不能定量描述层析正问题.Yomogida(1992)提出的菲涅尔带反演方法，在改善射线层析矩阵稀疏性的同时，兼顾了地震波在菲涅尔体内的有限频效应(刘玉柱，2010).不管是射线层析还是菲涅尔体层析，考虑到走时扰动对速度扰动的非线性程度较弱且对初始模型精度的要求较低，这类方法更多的还是利用地震数据的走时信息进行速度反演.(4) 经典的数据域波动方程层析反演即为全波形反演(FWI).由于采用波动方程作为正演算子，因而可以同时预测地震数据的运动学和动力学信息，反演时可灵活利用不同的地震属性(如走时、相位、振幅或波形).因此，FWI在实用化过程中发展出了多种解决方案，如相位/振幅目标泛函(Bedner et al., 2007)；瞬时相位/振幅/包络目标泛函(Bozda et al., 2011；Wu et al., 2014)，波场能量目标泛函(Jeon et al., 2014)；走时误差泛函(Luo and Schuster, 1991).然而，由于叠前数据(特别是陆上数据)的不完备(缺低频和大偏移距，信噪比低等)、正演模型的不完善(正演算子无法客观描述地震波在实际介质中的传播过程)、初始模型不准确(周期跳问题)、地震子波问题(空变)等因素影响，FWI常常难以反演出高分辨率的速度场，甚至是不收敛.FWI在国外海上资料探区的成功案例尚且不多，相对而言，国内陆上资料特别是复杂地表复杂构造低信噪比资料的FWI应用更需要一个过程.
上述前两类方法属于成像域速度反演方法，主要利用反射数据的走时信息，对初始模型的依赖性较低.该方法仅能反演速度模型的低波数成分，因而主要用于背景速度建模，为偏移成像服务及为更高精度的反演方法(如FWI)提供初始模型.后两类方法属于数据域速度反演方法，可选择利用走时、相位、振幅或波形属性，理论反演精度较高.但受到叠前数据复杂波场及低信噪比的影响，难以准确高效地实现全自动拾取，实际应用较为困难.从实用化的角度，本文关注成像域走时层析速度反演方法.
常规的成像域走时层析速度反演的基本过程是首先自动拾取成像道集ADCIGs(Angle Domain Common Image Gathers，角度域共成像点道集)或ODCIGs(Offset Domain Common Image Gathers，偏移距域共成像点道集)的剩余时差(RMO)，然后通过射线追踪计算射线路径并构造相应的层析线性方程，求解方程实现剩余时差的反投影从而估计速度更新量.通常而言，层析速度反演一般在偏移距域或角度域成像道集中实现.相对于偏移距域，角度域层析更有优势(射线出发角确定，无需费时且不稳定的两点射线追踪匹配地表偏移距，且避免了多路径等问题).基于射线理论的层析速度反演方法技术相对较为成熟，在实际应用中已经基本实现自动化(包括剖面上自动判断/拾取层析控制点、自动计算界面倾角、成像道集上自动筛选数据并计算RMO等).然而，射线理论的固有问题限定了该方法的反演精度、稳定性及适用范围，因而有必要进一步发展更加适合复杂速度模型的射线束层析反演方法.
与成像方法类似，在层析速度反演中也有介于常规射线理论和波动理论之间的方法，如胖射线层析(Vasco et al., 1995)、菲涅尔体层析(刘玉柱，2010)、高斯束层析(Popov et al., 2008, 2010)、波路径层析(Bakker et al., 2015)等.胖射线层析从射线层析向波动层析走近了一步，但其理论基础不够严格，仅仅是对波动理论感性认识的结果.基于菲涅尔体的波形层析方法是近年来发展较快的重要方法.它利用射线周围的波动性质来同时拟合走时和波场的振幅，其本质是有限频层析的一种.Semtchenok等(2009)最早提出了高斯束层析方法，利用高斯束展布范围内的射线合并建立一个层析方程，对于一个炮检对通过建立多个层析方程来减少层析方程的病态性，在一定程度上提高了反演精度和分辨率，且兼容了传统射线层析计算速度快、相对稳定的优点.但是Semtchenok等利用高斯束的局部波场振幅作为加权系数的推导过程的理论基础不够严格.鉴于高斯束传播算子在叠前深度偏移领域中具有较强的实用性(蔡杰雄等，2012；黄建平，2016；栗学磊和毛伟健，2016；刘玉柱，2014；岳玉波等，2012)，且能够高效地提取三维共方位-反射角度成像道集(Cai et al., 2013)，本文从高斯束偏移角度道集出发，在波动方程的一阶Born近似和Rytov近似下，推导成像域反射波走时层析方程及其敏感度核函数(以下简称核函数)，从而发展了一种新的成像域的波动方程线性化反演方法，并利用高斯束传播算子高效计算走时层析核函数，提供了一种新的成像域射线束层析反演思路.基于高斯束算子的成像域走时层析方法，与高斯束偏移相结合，可形成具有典型特征波成像特点、适应低信噪比数据的成像与建模工具.
2 成像域走时层析核函数的推导从角度域高斯束偏移(GBM)成像条件出发：
|
(1) |
(1) 式是频率域高斯束成像条件，其中pS和pR分别代表从炮点和检波点传播到成像点x的慢度矢量，表示波场在该点的传播方向. IGBM(x, θ, φ, ω)是频率域表示的角度成像结果(共方位-反射角成像道集)；S(x, pS, ω; xS)表示炮点出发的下行高斯束波场，R(x, pR, ω; xS)表示检波点出发的上行高斯束波场；x=(x, y, z)表示成像点坐标，θ, φ分别表示成像点的反射张角及方位角，*表示共轭.我们所熟知的单程波成像条件(在波场传播时没有显式地得到慢度矢量p)与高斯束成像条件类似，都是激励时间成像条件.
在波动方程的一阶Born近似下，波场U可以分解为背景波场U0和扰动波场ΔU，
|
(2) |
因此从成像条件(1) 式可以近似得到扰动像：
|
(3) |
其中，ΔS, ΔR分别为一阶Born近似散射场，S0, R0分别为炮点和检波点出发的背景波场.该成像条件表明，成像点x处像的扰动来自炮点端和检波点端两个分支的影响.
进一步地，根据Woodward(1992)和Jocker等(2006)，一阶Born近似散射场ΔS与ΔR可以表示为：
|
(4a) |
|
(4b) |
(4) 式中的VS和VR分别表示从炮点和从检波点到成像点的Born波路径；k0=ω/v0表示背景模型波数；Δv为待反演的速度扰动；G(x; p, ω, y)表示由高斯束表达的从点y到点x的格林函数，其下标S和R分别表示对应炮点波场和检波点波场; S0(y；pS, ω, xS)和R0(y; pR, ω, xS)分别表示在背景速度模型中从炮点和检波点传播到空间任意一点y的波场.
如果将(4) 式代入(3) 式，形式上可以给出像的扰动(左端项)与速度扰动(右端项)的关系式.但实际操作时，是不是能直接利用该式进行层析反演呢？显然不能！对比数据域层析反演可以分析：数据域反演利用的是正演数据与实测数据的差在某种范数下(一般是二范数)最小作为误差泛函，其实测数据是客观的，可直接用来做逼近标准；而如果直接利用(3) 式进行成像域反演，则由于客观上无法得到真实像IGBM(x, θ, φ, ω)从而无法得到扰动像ΔIGBM(x, θ, φ, ω)，因此其本质是利用一个未知的中间量来估计另一未知量Δv.直接利用像的扰动这个概念来进行反演缺乏严格的判断标准.从另一个角度分析，像的扰动是一个综合概念，其实质包括了走时(位置)扰动和振幅扰动.实际计算时需要将像的扰动退化为走时扰动(深度扰动)或振幅扰动，从而分别建立与速度扰动的关系式.由于像域的振幅扰动影响因素太复杂(抛开地表一致性因素、震源及检波器的耦合效应等因素，即使是在准确模型下的保幅成像，角度道集仍然存在AVA效应.因此角度道集上能用来层析反演的仅仅是其同相轴的剩余走时差(RMO)，不包括其振幅)，实际层析反演一般退化为仅利用走时扰动.
考虑退化到仅利用走时扰动进行层析反演，进一步地，将式(3) 重写为：
|
(5) |
整理得
|
(6) |
在波动方程的Rytov近似下，波场可以表示为u=(A0+ΔA)ei(φ0+Δφ).成像值是两个波场相关得到，因此同样可以表示I=(A0+ΔA)ei(φ0+Δφ)，I0(x, θ, φ, ω)=A0eiφ0，由于ΔA≪A0，则式(6) 可以进一步表示成：
|
(7) |
由于Δφ趋于零，则式(7) 左端项近似等于iΔφ，两边取虚部得：
|
(8) |
S0(x, pS, ω; xS)和R0(x, pR, ω; xS)分别表示在背景速度模型中从炮点和检波点传播到成像点x的波场.
进一步地，根据Spetzler和Snieder(2001, 2004)，单频相位扰动与单频走时扰动有近似关系：
|
(9) |
同时将(4) 式及(9) 式代入(8) 式整理，最终得到成像域单频走时扰动与速度扰动的关系式：
|
(10) |
其中，KF为单频走时层析核函数，其两个分支分别表示为：
|
(11a) |
|
(11b) |
需要说明的是，在背景模型中的波场可以表示成子波与格林函数的乘积，因此在推导得到(11) 式的过程中约去了子波项(假设子波不变)，仅剩下格林函数项.
由于实际操作时，走时扰动Δt是在时空域(角度域成像道集)测量得到的，与频率无关.因此，我们定义带限地震信号的走时扰动可以用单频走时扰动加权叠加(刘玉柱，2010)得到
|
(12) |
W(ω)表示归一化的加权函数.本文采用高斯函数
|
(13) |
最终得到成像域带限走时扰动与速度扰动的关系式：
|
(14) |
其中，Δt为高斯束偏移的走时扰动，也就是成像道集的RMO; KT为带限走时层析核函数，其两个分支分别表示为：
|
(15a) |
|
(15b) |
到此为止，基于高斯束角度道集，在波动方程的一阶Born近似和Rytov近似下导出了成像域带限走时层析方程(14) 及其对应的核函数表达式(15).从式(15) 可以看出，该核函数的计算主要是背景波场中格林函数的计算.利用高斯束积分计算格林函数是一种精度较高且计算量较小的实用化计算方式(以下简称利用高斯束传播算子计算成像域走时层析核函数的方法为高斯束层析方法).
成像域走时层析反演的实际操作流程可以完全借鉴常规射线层析反演的框架，两者的区别仅仅是在核函数的表达与计算上.关于单频和带限层析核函数的计算及特征分析将在第3节中通过数值试验分析总结.
需要说明的是，前述所得到的成像域带限走时层析方程及其核函数虽然是从高斯束偏移角度道集出发推导的，但其同样适用于其他偏移方法所提取的角度道集.
3 高斯束计算层析核函数及其精度分析成像域走时层析核函数的计算关键是背景波场格林函数的计算.具体计算有多种方法：全波场有限差分法(Zhao et al., 2005)可以处理复杂速度模型但十分耗时，且全波场的模拟包含了所有波至使得难以区分感兴趣的特征波场；单程波算子(Xie and Yang, 2007, 2008)具有较高的效率，适合处理大规模数据且可以模拟散射、绕射等波现象，其最大的缺点是对陡倾角的限制，不能模拟大角度及回转波场.高斯束是波动方程在射线中心坐标系下的高频近似解.与常规射线相比，高斯束克服了焦散问题及两点射线追踪问题；与全波场有限差分相比，高斯束非常高效和灵活；与单程波相比，高斯束没有倾角限制且可以模拟回转波.高斯束的缺点是仅能在光滑背景模型中传播，而成像域层析成像的目的就是建立适合于叠前深度偏移成像的光滑背景速度场.从这个角度考虑，高斯束算子是一个很好的选择.同时，由于高斯束传播算子已成功应用于叠前深度偏移成像方法(蔡杰雄等，2012)，能高效地提取三维方位-反射角度道集(Cai et al., 2013).因此，进一步地将高斯束算子应用到成像域走时层析成像中以高效计算格林函数，可以更好地将偏移成像与层析成像结合起来形成迭代反演流程，从而真正地体现偏移成像与速度建模一体化的思想.
3.1 高斯束计算格林函数前面式(11) 已给出单频走时层析核函数的表达式，计算该核函数的关键是计算背景模型格林函数.利用高斯束计算格林函数需要讨论其精度，并通过与常速情况下的格林函数解析解对比分析.
二维单频高斯束可以表达为：
|
(16) |
其中，s是中心射线传播弧长，n是垂直于中心射线的长度，τ是中心射线的走时，ω是圆频率，v是中心射线上垂点的速度，P和Q是动力学射线追踪求得的复数参数.通过高斯束积分得到光滑背景速度场中的格林函数(erveny et al., 1982)：
|
(17) |
其中，p为高斯束的初始慢度矢量；u(y, x, p, ω)为以直角坐标系参量表示的高斯束.通过(17) 式利用高斯束积分表达的格林函数受到频率及高斯束最大积分角度(高斯束的个数)的影响，下面通过一个常速介质模型数值试验来对比分析高斯束计算格林函数的精度.均匀介质中二维格林函数的解析表达式为：
|
(18) |
式中的R表示传播距离.设计一个模型：假设常速度为2 km·s-1.首先分析频率对计算格林函数的影响.图 1a，1b分别为计算频率40 Hz的高斯束积分所计算格林函数同格林函数解析解的实部、虚部的对比.图 1c，1d分别为计算频率20 Hz的高斯束积分所计算格林函数同格林函数解析解的实部、虚部的对比.图 1e，1f分别为40 Hz及20 Hz高斯束计算格林函数(实部)与解析解的相对误差.图中可以看到，频率为40 Hz时，当传播距离大于200 m时，高斯束积分所计算格林函数就已经很好地近似了格林函数解析解，其误差值基本为0；而当频率为20 Hz时，只有当传播距离大于350 m时，高斯束积分所计算格林函数才能很好地近似格林函数解析解，误差值基本为0.这一点验证了高斯束的高频假设.
|
图 1 不同频率高斯束计算格林函数与解析格林函数对比 (a)实部(40 Hz)；(b)虚部(40 Hz)；(c)实部(20 Hz)；(d)虚部(20 Hz)；(e)误差值(40 Hz)；(f)误差值(20 Hz). Fig. 1 Comparison of 20 and 40Hz Green′s function calculated using GB and analytical methods (a) and (c) Real part; (b) and (d) Imaginary part; (e) and (f) Errors. |
高斯束计算格林函数的精度与频率有关，同时和参与积分的高斯束个数有关，即(17) 式中的θ积分.取同上的模型，频率40 Hz，分别试验最大积分角度30°和60°所计算的格林函数(仅对比实部)与解析解的对比.
从图中可以看出，最大积分角度为30°时(图 2a)，高斯束计算格林函数与解析格林函数的一致性较差，直到传播距离大于200 m之后误差才基本为0；而当最大积分角度为60°时，两者的误差震荡较小，且在传播距离150 m以后即基本为0.可以看出，积分角度对高斯束计算格林函数的精度有影响，但影响较小.积分角度越大，参与的高斯束越多，所计算的格林函数越准确，越稳定，但这仅仅体现在震源附近，远场的拟合基本不受积分角度的影响.考虑到计算量的增加，实际计算时并非单纯地增加高斯束个数，仅需要在一定程度上满足精度要求，同时平衡计算效率.
|
图 2 不同最大张角(a)30°和(b)60°高斯束积分计算格林函数与解析格林函数对比 Fig. 2 Comparison of Green′s function calculated using GB with maximum open-angles 30° (a) and 60° (b) and analytical methods |
将3.1节计算的格林函数代入走时层析核函数(11) 式中即可实现高斯束计算层析核函数，并在常速情况下与解析解对比分析.
设计一个模型：假设常速3000 m·s-1, 传播距离1000 m，考虑单频20 Hz和40 Hz的核函数计算.为了简化式(11) 层析核函数的表达，省略格林函数对频率和空间坐标的显式表达，记GR=G0(x; ω, y)，GS=G0(y; ω, xS)，GS, R=G0(x; ω, xS)，将单频走时层析核函数简写为
|
(19) |
此时，
|
其中Ry, R表示从任意一点y到接收点R的传播距离, RS, R表示从炮点S到接收点R的传播距离, RS, y表示从炮点S到任意一点y的传播距离.则
|
(20) |
记时差
|
(21) |
利用该式作为常速情况下单频走时层析核函数的解析解，对比高斯束计算结果(图 3).
|
图 3 常速情况下单频高斯束计算核函数与解析解对比 (a) 20 Hz解析核函数；(b) 20 Hz高斯束计算核函数；(c) 40 Hz解析核函数；(d) 40 Hz高斯束计算核函数. Fig. 3 Comparison of monotonic sensitivity kernel calculated using GB and analytical methods with constant velocity (a) Analytical kernel and (b) GB kernel with 20 Hz; (c) Analytical kernel and (d) GB kernel with 40 Hz. |
为直观分析两者的差别，抽取图 4中横坐标x=1 km处的核函数纵截面数值进行对比.从图中可以看出，不同频率高斯束计算的核函数与解析解基本一致，完全满足后续走时层析速度反演的需求.
|
图 4 不同频率(a)20 Hz，(b)40 Hz高斯束计算核函数与解析解在x=1 km处纵截面的精度对比 蓝线表示解析核函数，红线表示高斯束计算核函数. Fig. 4 Comparison of accuracy between different frequency sensitivity kernels calculated using the GB method and the analytical methods at x=1000 m:(a) 20 Hz; (b) 40 Hz Blue line stands for kernel by analytical method and red line for GB. |
进一步地，根据(15) 式对频率加权积分可计算带限走时层析核函数，积分频段为
图 5和图 6可以看出，利用高斯束计算带限走时层析核函数与解析解基本一致.
|
图 5 高斯束计算的带限走时层析核函数与解析解对比 (a)解析带限走时层析核函数；(b)高斯束求得带限走时层析核函数. Fig. 5 Comparison of finite-frequency sensitivity kernel calculated using analytical methods (a) and GB (b) |
|
图 6 高斯束计算带限走时层析核函数与解析解在x=1000 m处纵截面的精度对比 蓝线表示解析核函数，红线表示高斯束计算核函数. Fig. 6 Comparison of accuracy between finite-frequency sensitivity kernels calculated using the GB method and the analytical methods at x=1000 m Blue line stands for kernel by analytical method and red line for GB. |
下面给出常速度梯度情况下(图 7)及层状介质情况下(图 8)的核函数计算，两种核函数的计算均通过高斯束积分得到，积分频段同上.
|
图 7 (a)常梯度速度模型和(b)常梯度模型带限核函数.震源位于(1 km, 0)，接收点位于(5.2 km, 0.9 km) Fig. 7 Finite-frequency sensitivity kernel of constant gradient velocity model. Source is at (1 km, 0) and receiver is at (5.2 km, 0.9 km) |
|
图 8 (a)层状速度模型和(b)层状模型带限核函数.震源位于(1 km, 0)，接收点位于(4.9 km, 3.1 km) Fig. 8 Finite-frequency sensitivity kernel of layered velocity model. Source is at (1 km, 0) and receiver is at (4.9 km, 3.1 km) |
从上面两节高斯束计算格林函数及进一步计算走时层析核函数的精度分析可以得到如下认识：
(1) 单频高斯束计算格林函数：除了在震源附近出现不稳定现象以外，在远离震源一定距离后的计算值基本与解析解一致；频率越高，不稳定区域越小；最大积分角度对高斯束计算格林函数的精度影响较小.
(2) 单频高斯束计算核函数：除了在核函数的中心线附近有较小误差外，不同频率高斯束计算的波路径与解析解基本一致.
(3) 高斯束计算带限核函数：除了在核函数的中心线附近有较小误差外，其他区域计算值与解析解基本一致，完全满足后续走时层析成像的需求.
(4) 不管是单频还是带限的走时层析核函数，其计算的核函数总体表现为中心值比较小，远离中心射线后逐渐增大，在某点达到最大后又逐渐减小直至为0.形状表现为双驼峰型.
4 成像域高斯束层析反演技术的具体实现与常规成像域射线层析方法一样，在成像域执行高斯束层析反演的基本思想是将成像道集的剩余时差沿着波传播的特定路径进行反投影，得到速度模型更新量，完成速度模型的一次迭代更新.研究共成像点道集剩余时间差与剩余深度差之间的关系，并以此关系为基础建立层析方程组，进而迭代求取剩余速度.
4.1 成像域高斯束层析反问题的建立当初始速度模型与真实模型接近时(低波数成分)，对于共成像点道集满足如图 10a所示的关系.图中，S是成像道集某道对应的炮点，R是对应的检波点，A是偏移速度的成像点，A′是正确速度时的成像位置，实线SAR是射线在偏移速度中的传播路径，虚线SA′R是射线在假设真速度中的传播路径，距离AB是拾取的剩余深度差，ψ是地层倾角，θ，φ分别是偏移入射角和理论入射角，当传播路径比较长，速度误差比较小时，我们近似认为θ=φ，粗线是共成像点位置.
|
图 9 带限核函数纵截面分布图 (a) 图 7b在x=3 km处的纵截面；(b) 图 8b在x=2.5 km处的纵截面. Fig. 9 Vertical profiles of finite-frequency sensitivity kernel (a) Cut slice of x=3 km in Fig. 7b; (b) Cut slice of x=2.5 km in Fig. 8b. |
|
图 10 (a)复杂介质传播路径和(b)成像点局部传播路径 Fig. 10 Ray path in complex media (a) and local ray path near image point (b) |
在反射点局部，可以假定慢度s是常数，又因为在远离反射点的地方两条传播路径比较接近，这样在反射点局部(如图 10b)，走时扰动量等于两条路径之差乘以慢度s.即
|
(22) |
即公式(14) 中的走时扰动，在成像道集中表现为剩余走时差RMO.其中
|
(23) |
将式(23) 代入式(22) 得到剩余深度差与走时残差的关系式：
|
(24) |
将式(24) 求得的RMO代入走时层析反演方程即可实现成像域层析反问题的建立.由此也可看出，成像域的层析反演需要在偏移剖面上自动拾取反射轴倾角ψ，以及成像道集上的剩余深度差Δz.
|
(25) |
A即是上述核函数计算得到的敏感度矩阵，Δm即是待反演的速度更新量.然而实际求解时，矩阵A一般是病态的，这是因为数据往往是有限的，并不足以约束所有的模型分量，破坏了反演的稳定性.因此，必须引入额外的信息对反问题进行正则化，即加入一些先验约束，使得解向所约束的方向发展.
从反演的角度，成像域层析归结为如下泛函的求极值问题：
|
(26) |
改造(26) 式的泛函，得到如下新的误差泛函：
|
(27) |
其中εd2, ε2C1, ε2C2为较小的正实数；mr为给定的先验模型，一般取mr=0或mr=mn，mn为前一次反演结果；εd2‖m－mr‖22为阻尼项，这项的引入使得每次迭代的模型更新量较小.ε2C1‖D1(m－mr)‖22和ε2C2‖D2(m－mr)‖22分别约束速度模型在横向和纵向上的光滑程度，D1、D2分别是对B样条基函数系数更新量在横向和纵向上的一阶差分算子.
对于更改后的误差泛函(27)，每次迭代需要求解的层析方程变为：
|
(28) |
利用LSQR方法求解矩阵方程组(28)，该方法是一种迭代的方法，可以在最小二乘意义下高效地求解大规模稀疏矩阵.至此，建立了带有正则化项的成像域高斯束层析成像反问题.
4.2 高斯束层析反演流程成像域高斯束层析反演的流程同常规射线层析反演的流程基本一致.
(1) 基于初始速度模型进行高斯束叠前深度偏移(GBM)，并输出(方位)反射角度道集(ADCIGs)；
(2) 在深度偏移剖面上全自动解释反射界面，拾取层析控制点；
(3) 参数化速度模型；
(4) 在界面约束下，在ADCIGs上全自动拾取相应的同相轴的剩余深度差/时间差；
(5) 构建Frechet导数矩阵；
(6) 建立走时误差泛函(可以考虑各种约束)；
(7) 利用并行化的LSQR算法求解法方程，得到速度更新量；
(8) 更新速度模型；
(9) 重新进行GBM.视ADCIGs道集是否拉平决定是否进行下一轮的速度反演.
5 理论模型及实际数据应用该模型断层发育，地层高陡，速度横向变化较大(图 11a).横纵向采样点为731×550，采样间隔均为10 m.利用声波正演得到叠前炮记录(图 11b)，炮间距和道间距都是10 m，中心放炮，每炮361道.利用等梯度速度模型作为层析初始模型(图 12a).对初始模型进行高斯束叠前深度偏移，输出偏移剖面(图 12b)及角度道集(图 13b)，可以看出初始模型对应的偏移剖面上的各层位成像深度并不准确，绕射波也没有完全收敛.图 13是CDP400位置处提取的初始模型角度道集和层析后角度道集对比，由于速度偏低，初始角度道集同相轴上翘；经过层析反演更新速度模型，偏移后角度道集得到了拉平.图 14是高斯束层析经过五次迭代后得到的速度模型及其相应的高斯束偏移剖面，图 15是常规射线层析经过八次迭代得到的速度模型及其高斯束偏移剖面.与常规射线层析对比可知，本文发展的高斯束层析方法能够提供更加丰富的速度信息，迭代次数也小于常规射线层析.更新后偏移提取的角度道集更加接近真实速度模型下偏移提取的角度道集，反射界面归位到正确的深度位置，绕射波收敛，断层位置聚焦更好，说明高斯束层析对复杂构造模型的速度反演精度优于常规射线层析方法.
|
图 11 某地区复杂速度模型(a)及其相应炮集记录(b) Fig. 11 (a) Complex velocity model in an area and (b) shot record |
|
图 12 初始速度模型(a)及其高斯束叠前深度偏移剖面(b) Fig. 12 (a) Initial velocity model and (b) PSDM profile by GBM |
|
图 13 CDP 400角度道集 (a)初始角度道集；(b)正确速度模型提取角度道集；(c)高斯束层析提取角度道集；(d)常规射线层析提取角度道集. Fig. 13 Angle CIGs at CDP=400 (a) Initial CIG; CIG extracted from; (b) Correct velocity model; (c) GB tomography; (d) Ray-based tomography. |
|
图 14 高斯束层析更新速度场(a)及其层析后高斯束偏移剖面(b) Fig. 14 (a) Velocity field by GB tomography and (b) GBM profile |
|
图 15 常规射线层析更新速度场(a)及其层析后高斯束偏移剖面(b) Fig. 15 (a) Velocity field by ray-based tomography and (b) GBM profile |
实际资料来自中国南方复杂山前带某工区.该区地形起伏剧烈(图 16中的曲线代表高程)，悬崖峭壁密布，山高谷深，碳酸盐岩大面积裸露，喀斯特岩溶地貌较发育，形成岭谷相间地形地貌.地震数据信噪比非常低.受地表地震地质条件及复杂构造的影响，深部构造的勘探工作困难重重，特别是在深部构造高隆起部位通常反射能量弱、且信噪比低，波组连续性差，造成落实构造形态难度大.
|
图 16 中国南方某山前带工区实际资料某单炮 Fig. 16 Shot record of a piedmont area in South China |
本文的研究工作是在前期处理人员进行了预处理、叠前时间偏移及深度域初始建模及偏移成像的基础上开展的，主要通过高斯束层析反演后的速度模型进行高斯束叠前深度偏移来体现本文方法的实际应用效果.
从图 17层析反演的速度模型上看，高斯束层析结果体现了速度模型的更多细节信息.鉴于高斯束偏移的抗噪性较强，进一步地，通过高斯束偏移提高成像质量.
|
图 17 某实际数据深度域初始建模(a)，常规射线层析建模(b)及高斯束层析建模(c)结果对比 Fig. 17 Comparison of (a) initial velocity model of field data, velocity model by (b) ray-based tomography and (c) GB tomography |
从图 18可以看出，同样利用高斯束偏移方法，高斯束层析成像的速度模型对应的偏移结果信噪比更高，同相轴更连续，整体成像质量更高.而从图 19的不同偏移算法处理可以看出，高斯束偏移结果信噪比明显高于Kirchhoff偏移结果，整体成像质量有较明显的提升，对后续的构造落实有较好的指导意义.从该实际资料处理可以认识到，结合高斯束层析与高斯束偏移技术的处理方式，将能更好地适应山前带探区等复杂构造低信噪比资料，提供更高质量的成像结果.
|
图 18 实际数据分别利用射线层析模型(a)及高斯束层析模型(b)的高斯束偏移结果 Fig. 18 GBM profiles by (a) ray-based tomography and (b) GB tomography |
|
图 19 实际数据利用高斯束层析模型分别进行(a)Kirchhoff偏移和(b)高斯束偏移结果 Fig. 19 Profiles of field data migrated by Kirchhoff (a) and GBM (b) |
本文从角度域高斯束偏移成像条件出发，在波动方程的一阶Born近似和Rytov近似下，推导给出了成像域走时层析成像方程及其核函数表达式，包括单频及带限核函数.利用带限层析核函数替代常规射线层析核函数(该核函数为常数1)，可以改进层析反演精度，加快反演收敛.进一步地，利用高斯束传播算子计算该核函数既能保证计算精度同时兼顾射线层析的高效性.
成像域速度反演方法主要利用反射数据的走时信息，对初始模型的依赖性较低.该方法仅能反演速度模型的低波数成分，因而主要用于背景速度建模，为偏移成像服务及为更高精度的反演方法(如FWI)提供初始模型.本文发展的高斯束层析成像方法利用高斯束传播算子计算成像域走时层析核函数，提供了一种新的成像域波动方程线性化近似层析反演方法.基于高斯束算子的成像域走时层析方法，与高斯束偏移技术相结合，可形成具有典型特征波成像特点的、适应低信噪比数据的成像与建模工具，真正体现偏移成像与速度建模一体化的处理思想.
| Bakker P, Gerritsen S, Cao Q. 2015. 3D RTM-based wave path tomography tested at a realistic scale.//77th Conference & Exhibition-Workshops, EAGE. Extended Abstracts, WS05-C02. | |
| Bednar J B, Shin C, Pyun S. 2007. Comparison of waveform inversion, part 2:phase approach. Geophysical Prospecting, 55(4): 465-475. DOI:10.1111/gpr.2007.55.issue-4 | |
| Bozda E, Trampert J, Tromp J. 2011. Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements. Geophysical Journal International, 185(2): 845-870. DOI:10.1111/gji.2011.185.issue-2 | |
| Cai J X, Fang W B, Wang H Z. 2012. Realization and application of Gaussian beam depth migration. Geophysical Prospecting for Petroleum (in Chinese), 51(5): 119-125. | |
| Cai J X, Fang W B, Wang H Z. 2013. Azimuth-opening angle domain imaging in 3D Gaussian beam depth migration. Journal of Geophysics and Engineering, 10(2): 025013. DOI:10.1088/1742-2132/10/2/025013 | |
| erveny V, Popov M M, Psencik 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 | |
| Feng B. 2015. Data domain wave equation tomography for velocity inversion (in Chinese). Shanghai: Tongji University. | |
| Huang J P, Yang J D, Li Z C, et al. 2016. An amplitude-preserved Gaussian beam migration based on wave field approximation in effective vicinity under irregular topographical conditions. Chinese Journal of Geophysics (in Chinese), 59(6): 2245-2256. DOI:10.6038/cjg20160627 | |
| Jeon S, Ha W, Kwon J, et al. 2014. Full waveform inversion using the energy objective function.//76th EAGE Conference and Exhibition. Expanded Abstracts, Tu E106 02. | |
| Jocker J, Spetzler J, Smeulders D, et al. 2006. Validation of first-order diffraction theory for the traveltimes and amplitudes of propagating waves. Geophysics, 71(6): T167-T177. DOI:10.1190/1.2358412 | |
| Li X L, Mao W J. 2016. Multimode and multicomponent Gaussian beam prestack depth migration. Chinese Journal of Geophysics (in Chinese), 59(8): 2989-3005. DOI:10.6038/cjg20160822 | |
| Liu Y Z. 2010. Theory and applications of Fresnel volume seismic tomography (in Chinese). Shanghai: Tongji University. | |
| Liu Y Z, Xie C, Yang J Z. 2014. Gaussian beam first-arrival waveform inversion based on Born wavepath. Chinese Journal of Geophysics (in Chinese), 57(9): 2900-2909. DOI:10.6038/cjg20140915 | |
| Luo Y, Schuster G T. 1991. Wave-equation traveltime inversion. Geophysics, 56(5): 645-653. DOI:10.1190/1.1443081 | |
| Popov M M, Semtchenok N M, Popov P M, et al. 2008. Reverse time migration with Gaussian beams and velocity analysis applications.//70th EAGE Conference and Exhibition incorporating SPE EUROPEC 2008. Extended Abstracts, F048. | |
| Popov M M, Semtchenok N M, Popov P M, et al. 2010. Depth migration by the Gaussian beam summation method. Geophysics, 75(2): S81-S93. DOI:10.1190/1.3361651 | |
| Semtchenok N M, Popov M M, Verdel A R. 2009. Gaussian beam tomography.//71st EAGE Conference & Exhibition. Amsterdam:EAGE. | |
| Spetzler G, Snieder R. 2001. The effect of small-scale heterogeneity on the arrival time of waves. Geophysical Journal International, 145(3): 786-796. DOI:10.1046/j.1365-246x.2001.01438.x | |
| Spetzler G, Snieder R. 2004. The Fresnel volume and transmitted waves. Geophysics, 69(3): 653-663. DOI:10.1190/1.1759451 | |
| Tarantola A. 1984. Inversion of reflection data in the acoustic approximation. Geophysics, 19(8): 1259-1266. | |
| Vasco D W, Peterson J E Jr, Majer E L. 1995. Beyond ray tomography:Wavepaths and Fresnel volumes. Geophysics, 60(6): 1790-1804. DOI:10.1190/1.1443912 | |
| Woodward M J. 1992. Wave-equation tomography. Geophysics, 57(1): 15-26. DOI:10.1190/1.1443179 | |
| Wu R S, Hu C H, Wang B F. 2014. Nonlinear sensitivity operator and inverse thin-slab propagator for tomographic waveform inversion.//84th Annual International Meeting, SEG. Expanded Abstracts, 928-933. | |
| Xie X B, Yang H. 2007. A migration-velocity updating method based on the shot-index common image gather and finite-frequency sensitivity kernel.//77th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 2767-2771. | |
| Xie X B, Yang H. 2008. The finite-frequency sensitivity kernel for migration residual moveout and its applications in migration velocity analysis. Geophysics, 73(6): S241-S249. DOI:10.1190/1.2993536 | |
| Yomogida K. 1992. Fresnel zone inversion for lateral heterogeneities in the earth. Pure and Applied Geophysics, 138(3): 391-406. DOI:10.1007/BF00876879 | |
| 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 | |
| Zhao L, Jordan T H, Olsen K B, et al. 2005. Frechet kernels for imaging regional earth structure based on three-dimensional reference models. Bulletin of the Seismological Society of America, 95(6): 2066-2080. DOI:10.1785/0120050081 | |
| 蔡杰雄, 方伍宝, 王华忠. 2012. 高斯束深度偏移的实现与应用研究. 石油物探, 51(5): 119-125. | |
| 冯波. 2015. 数据域波动方程层析速度反演方法研究[博士论文]. 上海: 同济大学. | |
| 黄建平, 杨继东, 李振春, 等. 2016. 基于有效邻域波场近似的起伏地表保幅高斯束偏移. 地球物理学报, 59(6): 2245–2256. DOI:10.6038/cjg20160627 | |
| 栗学磊, 毛伟建. 2016. 多波多分量高斯束叠前深度偏移. 地球物理学报, 59(8): 2989–3005. DOI:10.6038/cjg20160822 | |
| 刘玉柱. 2010. 菲涅尔体地震层析成像理论与应用研究[博士论文]. 上海: 同济大学. | |
| 刘玉柱, 谢春, 杨积忠. 2014. 基于Born波路径的高斯束初至波波形反演. 地球物理学报, 57(9): 2900–2909. DOI:10.6038/cjg20140915 | |
| 岳玉波, 李振春, 钱忠平, 等. 2012. 复杂地表条件下保幅高斯束偏移. 地球物理学报, 55(4): 1376–1383. DOI:10.6038/j.issn.0001-5733.2012.04.033 | |