E-mail Alert Rss
 

物探与化探, 2018, 42(5): 977-989 doi: 10.11720/wtyht.2018.1137

方法研究·信息处理·仪器研制

基于方位—反射角度道集的高斯束层析速度建模方法及实现

蔡杰雄

中国石油化工股份有限公司石油物探技术研究院,江苏 南京 211103

Method and application of Gaussian beam velocity tomography based on azimuth-reflection angle domain common imaging gathers

CAI Jie-Xiong

SINOPEC Geophysical Research Institute,Nanjing 211103,China

收稿日期: 2017-08-24   修回日期: 2018-05-17   网络出版日期: 2018-10-05

基金资助: 国家科技重大专项.  2016ZX05014-001-002

Received: 2017-08-24   Revised: 2018-05-17   Online: 2018-10-05

作者简介 About authors

蔡杰雄(1983-),男,同济大学固体地球物理学博士毕业,副研究员,现在中石化石油物探技术研究院从事地震反演成像研究工作。Email:caijx.swty@sinopec.com 。

摘要

层析反演是速度建模中最重要的方法之一,结合偏移成像在成像域进行波动方程线性化走时层析速度建模是当前比较实用有效且精度较高的技术组合。文中首先给出了高斯束偏移提取方位—反射角度道集的方法,之后从高斯束偏移角度道集出发,在波动方程的一阶Born近似和Rytov 近似下,推导了成像域波动方程线性化走时层析方程及其显式表达的层析核函数,并利用高斯束传播算子计算该核函数。基于高斯束传播算子的偏移成像与层析成像相结合进行深度域速度建模迭代及偏移成像,体现了速度建模与成像一体化的思想。数值计算及实际数据应用证明了基于高斯束传播算子的层析成像与偏移成像方法的有效性。

关键词: 高斯束 ; 走时层析 ; 核函数 ; 偏移成像 ; 角度道集

Abstract

Tomography is one of the most important velocity building methods.Travel time tomography in image domain, implemented with migration, is widely used in velocity model building currently. The authors first introduced the method to output the azimuth-reflection angle gathers, and then, under the assumption of the first-order Born and Rytov approximation of wave equation,started with the imaging condition of Gaussian Beam Migration to derive the linear relation between traveltime perturbation and velocity perturbation in the image domain,with which the authors constructed the explicit expression of kernel function for the wave equation traveltime tomography and established the traveltime tomography equation.The key to computing the kernel is how to compute the Green function in the background model.Making use of the Gaussian beam propagation operator to compute the kernel function can be flexible and efficient.Together with the implementation of Gaussian beam propagation operator in migration,the authors truly realized the integrated technological process of velocity building and migration.Numerical tests and field data application demonstrate that the Gaussian-beam-propagator based traveltime tomography in image domain is effective.

Keywords: Gaussian beam ; traveltime tomography ; kernel function ; migration ; angle gathers

PDF (10265KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

蔡杰雄. 基于方位—反射角度道集的高斯束层析速度建模方法及实现. 物探与化探[J], 2018, 42(5): 977-989 doi:10.11720/wtyht.2018.1137

CAI Jie-Xiong. Method and application of Gaussian beam velocity tomography based on azimuth-reflection angle domain common imaging gathers. Geophysical and Geochemical Exploration[J], 2018, 42(5): 977-989 doi:10.11720/wtyht.2018.1137

0 引言

地震成像所利用的数据主要是反射波。反射波场由背景速度扰动和高波数速度扰动共同决定,因此基于反射波的地震成像分为叠前深度偏移成像(确定高波数扰动)和反射波层析成像(确定背景速度扰动)两个步骤。偏移成像用于定位地下反射位置,实现了从数据空间到成像空间的映射,是正过程;层析成像用于反演传播路径上的速度场,实现了从成像空间到模型空间的映射,是反过程。当然,这两个步骤也可以耦合在一起进行数据域的迭代反演成像(如全波形反演,FWI)。相比于常规射线走时层析,结合偏移成像在成像域进行波动方程线性化走时层析速度反演是当前比较实用有效且精度较高的技术组合,两者的核心都归结于格林函数的计算。

层析速度建模的实施涉及到两个问题,即正问题所采用的理论是基于射线还是基于波动,以及所输入的数据空间(成像域或数据域)、数据属性(走时、振幅、相位或波形)和数据类型(折射波、初至波、透射波、反射波等)。对现有主流的层析成像方法进行分类,根据执行域可以分为成像域和数据域层析成像;根据正问题可以分为基于波动理论和基于射线理论的层析成像;根据利用数据属性可以分为波形和走时层析成像等;根据利用数据类型可以分为折射波、初至波、透射波和反射波层析成像等。根据正问题及执行域的不同,组合衍生出具有不同反演能力和需求的方法,其中:①基于射线理论的成像域层析速度反演方法发展较为成熟,在实际应用中已经实现自动化层析。然而,射线理论的固有问题(阴影区、多路径问题等)限定了该方法的反演精度和在复杂速度模型中的适用性,因而反演效果仍然有待改善(这也是本论文考虑的方向)。②波动方程偏移速度分析(WEMVA)方法以波动方程作为正演算子,克服了射线理论的高频近似假设,可适用复杂模型的速度估计,其反演精度理论上优于射线类层析方法,但是,WEMVA方法的计算量巨大,更存在地下照明不均衡及泛函梯度假象等问题,影响反演的收敛速度甚至不收敛,在实际应用中较少推广,性价比不高。③数据域射线层析与成像域射线层析的理论问题一样,基于高频射线理论,但反演精度不高,只能反演大尺度的速度异常。为提高反演精度及改善层析矩阵的稀疏性,出现了胖射线层析[1]、菲涅尔体层析[2,3]等。不管是射线层析还是菲涅尔体层析,考虑到走时扰动对速度扰动的非线性性程度较弱且对初始模型精度要求较低,这类方法更多的还是利用地震数据的走时信息进行速度反演。④经典的数据域波动方程层析反演即为全波形反演(FWI)。由于采用波动方程作为正演算子,可以同时预测地震数据的运动学和动力学信息,反演时可灵活利用不同的地震属性(如走时、相位、振幅或波形),因此,FWI在实用化过程中发展出了多种解决方案,如相位/振幅目标泛函[4]、瞬时相位/振幅/包络目标泛函[5,6,7]、波场能量目标泛函[8]、走时误差泛函[9]。然而,由于叠前数据(特别是陆上数据)的不完备(缺低频和大偏移距,信噪比低等)、正演模型的不完善(正演算子无法客观描述地震波在实际介质中的传播过程)、初始模型不准确(周期跳问题)、地震子波问题(空变)等因素影响,FWI常常难以反演出高分辨率的速度场,甚至是不收敛。FWI在国外海上资料探区的成功案例不多,相对而言,国内陆上资料特别是复杂地表复杂构造低信噪比资料的FWI应用更需要一个过程。

上述前两类方法属于成像域速度反演方法,主要利用反射数据的走时信息,对初始模型的依赖性较低;该方法仅能反演速度模型的低波数成分,因而主要用于背景速度建模,为偏移成像服务及为更高精度的反演方法(如FWI)提供初始模型。后两类方法属于数据域速度反演方法,可选择利用走时、相位、振幅或波形属性,理论反演精度较高,但受到叠前数据复杂波场及低信噪比的影响,难以准确高效地实现全自动拾取,实际应用较为困难。从实用化的角度,笔者关注成像域走时层析速度反演方法及技术。

成像域走时层析速度反演的输入是偏移成像道集,需要在此道集上拾取剩余时差(RMO)。目前工业界在进行层析速度反演时所输入的成像道集主要是利用Kirchhoff叠前深度偏移技术提取的偏移距域共成像点道集(ODCIGs),但Kirchhoff偏移所利用的偏移距参数是定义在数据面(地表)的炮检点距离而非地下成像点的信息,对于复杂构造区,由于忽略了不同偏移距矢量对应不同方位的波传播,且基于单旅行时的Kirchhoff偏移无法准确描述多路径从而导致成像质量不佳,成像道集产生假象。而通过计算成像点处的方位—反射角进而提取共方位—反射角度成像道集可以避免偏移距域成像道集的这些缺陷。方位—反射角度道集提供了地下成像点的初始射线追踪方向,这使得层析速度反演向地表实施射线追踪的过程相较于偏移距域更加地自然和有效,且避免了多路径等问题。角度域成像道集可以通过Kirchhoff叠前深度偏移[10,11]及高斯束偏移[12,13]等射线类偏移方法获得,也可通过波动方程偏移[14]和逆时偏移[15]等波动类偏移方法获得。通过波动类偏移方法提取三维地震数据的方位—反射角度道集的计算代价较大,因此其在工业界中并没有得到大规模的应用。考虑到射线类偏移方法在提取方位—反射角度道集上相对容易和代价较小,笔者利用高斯束叠前深度偏移提取地下局部方位—反射角度道集。

进而,基于高斯束偏移提取的角度道集,笔者在波动方程的一阶Born近似和Rytov近似下,推导成像域反射波走时层析方程及其敏感度核函数(以下简称核函数),并利用高斯束传播算子计算该核函数,从而发展一种新的成像域波动方程线性化走时层析反演方法(简称基于高斯束传播算子的成像域波动方程线性化走时层析反演方法为高斯束层析)。基于高斯束算子的成像域走时层析方法与高斯束偏移相结合,可形成具有典型特征波成像特点、适应低信噪比数据的成像与建模工具。

1 高斯束偏移提取方位—反射角成像道集

在叠前深度偏移过程中提取方位—反射角度道集是进行后续成像域高斯束层析速度反演的必备过程。计算地下成像点的方位—反射角的一种有效且高效的方法是分别估算从炮点和检波点出发到达成像点处的波场传播方向。一旦获得两者的传播方向,即可通过简单的向量代数运算得到方位角和反射角(张角)。对于射线类偏移方法而言,这个过程相对容易,只需通过旅行时场的空间导数分别计算震源波场的射线慢度ps和检波点波场射线慢度pr即可。而对高斯束函数的解析求解,即可方便高效地计算旅行时场的空间导数。如图1定义所示,方位角φ和反射张角θ可以通过如下的向量代数运算得到:

cosθ=cos(2α)=ps·pr|ps||pr|,

且,

cosφ=(z×psr)·(z×x)|z×psr||z×x|=(psr×z)·y|psr×z|

其中,

psr=ps-pr,x=(1,0,0),y=(0,1,0),z=(0,0,1)

图1

图1   定义地下成像点的方位—反射角示意

从炮点S出发和从检波点G出发到达成像点R处的射线慢度矢量分别表示为pspr,通过在成像点附近的局部空间中定义方位角(旋转角)和反射角(张角)


式(1a)中的反射张角θ定义为反射角α的两倍。式(1b)表示的方位角φ是相对于参考方向单位矢量x而言的。文中,我们定义参考方向沿着坐标轴x正方向(图1),而方位角φ对应的几何意义是局部炮检点方向矢量psr在水平面x-y上的投影沿着x轴逆时针方向的旋转角。这里我们认为方位角的参考方向固定不变,如文中的x轴方向(正东方向,当然也可以选择y轴方向,即正北方向为参考方向),该方向不随着局部反射平面[ps, pr]的变化而变化。

利用三维频率域高斯束函数[16,17]在射线中心坐标系中表达波动方程:

u(s,q1,q2)=v(s)detQ(s0)v(s0)detQ(s)·expiωτ(s)+12qTP(s)Q(s)q,

其中:ω表示圆频率,v(s)和τ(s)分别代表中心射线的速度和旅行时。qT=(q1,q2),其中q1,q2是射线中心坐标系中垂直于中心射线的两个坐标分量。矩阵P(s)和Q(s)是一个2×2的复数矩阵,沿着中心射线通过动力学射线追踪方程求解得到,表征动力学射线追踪参量。高斯束函数的具体求解在文献[16,17]中已经有详细探讨,这里不再赘述。

图2,高斯束邻域内任意一点Q的走时可由其对应的中心射线上的点R的旅行时及动力学射线追踪参量表示,进一步地、通过旅行时场的空间导数得到点Q的高斯束波场慢度矢量:

P(Q)=T(Q)x

图2

图2   三维射线中心坐标系示意

高斯束有效范围内任意一点Q的旅行时可用其对应的中心射线上点R的旅行时解析表达


分别得到炮点波场和检波点波场的慢度矢量后,根据反射角和方位角的定义式(1),即可求得成像点处的方位—反射角,从而输出方位—反射角度成像道集。详见参考文献[11],此处不再赘述。

2 成像域高斯束层析方法

进一步地,从角度域高斯束偏移(GBM)成像条件出发:

IGBM(x,θ,φ,ω)=S(x,ps,ω;xs)R*(x,pr,ω;xs),

上式是引入了高斯束传播方向p的频率域高斯束成像条件。IGBM(x, θ, φ, ω)是频率域表示的角度成像结果(共方位—反射角成像道集);S(x, ps, ω;xs)表示炮点出发的下行高斯束波场、R(x, pr, ω;xs)表示检波点出发的上行高斯束波场;x=(x, y, z)表示成像点坐标,θ, φ分别表示成像点的反射张角及方位角,*表示共轭。我们所熟知的单程波成像条件(在波场传播时没有显式地得到方向特征p)与高斯束成像条件类似,都是激励时间成像条件。

在波动方程的一阶Born近似下,波场U可以分解为背景波场U0和扰动波场ΔU:

U=U0+ΔU,

因此从成像条件式(3)可以近似得到扰动像:

ΔIGBM(x,θ,φ,ω)ΔS(x,ps,ω;xs)R0*(x,pr,ω;xs)+S0(x,ps,ω;xs)ΔR*(x,pr,ω;xs)

式中:ΔS、ΔR分别为一阶Born近似散射场,S0R0分别为炮点和检波点出发到成像点的背景波场。该成像条件表明,成像点x处像的扰动来自炮点端和检波点端两个分支的影响。

根据文献[18,19],一阶Born近似散射场ΔS与ΔR可以表示为:

ΔS(x,ps,ω;xs)=vS2k02Δv(y)v0(y)Gs(x,ps,ω;y)S0(y,ps,ω;xs)dy,

ΔR*(x,pr,ω;xs)=vR2k02Δv(y)v0(y)Gr(x,pr,ω;y)R0*(y,pr,ω;xs)dy

上面两式中的积分范围vSvR分别表示从炮点和从检波点到成像点的Born波路径;k0=ω/v0表示背景模型波数;Δv为待反演的速度扰动;S0(y, ps, ω;xs)和R0(y, pr, ω;xs)分别表示在背景速度模型中从炮点和检波点传播到空间任意一点y的波场;G(x, p, ω;y)表示由高斯束表达的从点y到点x的格林函数。

如果将式(6)代入式(5),形式上可以给出像的扰动(左端项)与速度扰动(右端项)的关系式。但实际操作时,显然是不能直接利用该式进行层析反演。对比数据域层析反演可以分析:数据域反演利用的是正演数据与实测数据的差在某种范数下(一般是二范数)最小作为误差泛函,其实测数据是客观的,可直接用来做逼近标准;而如果直接利用式(5)进行成像域反演,则由于客观上无法得到真实像IGBM(x, θ, φ, ω)从而无法得到扰动像ΔIGBM(x, θ, φ, ω),因此其本质是利用一个未知的中间量来估计另一未知量Δv。直接利用像的扰动这个概念来进行反演缺乏严格的判断标准。从另一个角度分析,像的扰动是一个综合概念,其实质包括了走时(位置)扰动和振幅扰动。实际计算时需要将像的扰动退化为走时扰动(深度扰动)或振幅扰动,从而分别建立与速度扰动的关系式。由于像域的振幅扰动影响因素太复杂,实际层析反演一般退化为仅利用走时扰动。

考虑退化到仅利用走时扰动进行层析反演,将式(5)重写为

ΔIGBM(x,θ,φ,ω)ΔS(x,ps,ω;xs)S0(x,ps,ω;xs)S0(x,ps,ω;xs)R0*(x,pr,ω;xs)+ΔR*(x,pr,ω;xs)R0*(x,pr,ω;xs)S0(x,ps,ω;xs)R0*(x,pr,ω;xs)=I0(x,θ,φ,ω)·ΔS(x,ps,ω;xs)S0(x,ps,ω;xs)+ΔR*(x,pr,ω;xs)R0*(x,pr,ω;xs),

整理得:

I(x,θ,φ,ω)I0(x,θ,φ,ω)-1ΔS(x,ps,ω;xs)S0(x,ps,ω;xs)+ΔR*(x,pr,ω;xs)R0*(x,pr,ω;xs)

在波动方程Rytov近似下,波场可以表示为u=(A0A) ei(φ0+Δφ)。成像值是两个波场相关得到,因此同样可以表示I=(A0A) ei(φ0+Δφ),I0(x, θ, φ, ω)=A0eiφ0,由于ΔAA0,则式(8)可以进一步表示成

eφ-1ΔS(x,ps,ω;xs)S0(x,ps,ω;xs)+ΔR*(x,pr,ω;xs)R0*(x,pr,ω;xs)

由于Δφ趋于零,则式(9)左端项近似等于iΔφ,两边取虚部得:

Δφ(x,θ,φ,ω)=ImΔS(x,ps,ω;xs)S0(x,ps,ω;xs)+ImΔR*(x,pr,ω;xs)R0*(x,pr,ω;xs)

根据Spetzler和Snieder[20],单频相位扰动与单频走时扰动有近似关系:

Δt(x,θ,φ,ω)=Δφ(x,θ,φ,ω)ω

同时将式(11)及式(6)代入式(10)整理,最终得到成像域单频走时扰动与速度扰动的关系式:

Δt(x,θ,φ,ω)=vSKSF(y,x,ps,ω;xs)Δv(y)dy+vRKRF(y,x,pr,ω;xs)Δv(y)dy,

其中,Δt(x, θ, φ, ω)为高斯束偏移的走时扰动,也就是成像道集的RMO。KF为单频走时层析核函数,其两个分支分别表示为:

KSF(y,x,ps,ω;xs)=Im2ωv03(y)Gs(x,ps,ω;y)G0(y,ps,ω;xs)G0(x,ps,ω;xs),

KRF(y,x,pr,ω;xs)=Im2ωv03(y)Gr(x,pr,ω;y)G0*(y,pr,ω;xs)G0*(x,pr,ω;xs)

上式可以看出,成像域走时层析核函数的表现形式与Liu[21]给出的数据域菲尼尔体走时层析核函数的表现形式类似,其本质都是有限频核函数,其计算关键是背景速度下格林函数的计算。

需要说明的是,在背景模型中的波场可以表示成子波与格林函数的乘积,因此在推导得到式(13)的过程中约去了子波项(假设子波不变),仅剩下格林函数项。

由于实际操作时,走时扰动Δt是在时空域(角度域成像道集)测量得到的,与频率无关。因此,我们定义带限地震信号的走时扰动可以用单频走时扰动加权叠加[3]得到

Δt=ω0-Δωω0+ΔωW(ω)Δt(ω)dω,

W(ω)表示归一化的加权函数,文中采用高斯函数g(ω)= 12πδexp -(ω-ω0)22δ2,以ω0为高斯函数的期望值(对称中心),Δω为高斯函数的标准差(半宽度)。则

W(ω)=g(ω)ω0-Δωω0+Δωg(ω)dω,(15)

最终得到成像域带限走时扰动与速度扰动的关系式:

Δt(x,θ,φ)=vSKST(y;x,ps;xs)Δv(y)dy+vRKRT(y;x,pr;xs)Δv(y)dy

其中,KT为带限走时层析核函数,其两个分支分别表示为

KST(y,x,ps;xs)=W(ω)Im2ωv03(y)Gs(x,ps,ω;y)G0(y,ps,ω;xs)G0(x,ps,ω;xs)dω,

KRT(y,x,pr;xs)=W(ω)Im2ωv03(y)Gr(x,pr,ω;y)G0*(y,pr,ω;xs)G0*(x,pr,ω;xs)dω

到此为止,基于高斯束角度道集、在波动方程的一阶Born近似和Rytov近似下导出了成像域带限走时层析方程式(16)及其对应的核函数表达式(17)。从式(17)可以看出,该核函数的本质是有限频核函数,其求解主要是背景波场中格林函数的计算。利用高斯束积分计算格林函数是一种精度较高且计算量较小的实用化计算方式。

下面通过数值试验分析单频和带限走时层析核函数的特征。设计背景模型参数:网格个数Nx=Ny=Nz=201,网格间距dx=dy=dz=10.0 m,速度为3 000 m/s,震源和检波点坐标分别为(500,1000,1000)和(1500,1000,1000),带限走时层析核函数的频率范围取(0,40 Hz),单频核函数计算取中心频率20 Hz。图3给出了对应的单频及带限走时层析核函数。从图3g和图3h上可以看出,核函数在炮检连线的中心线上的敏感度为零,这与常规射线层析的核函数仅分布在射线上的假设完全相反。刘玉柱[3]详细分析了这种差异,并指出了常规射线层析之所以仍然能取得较好效果的原因。

图3

图3   不同频率单频和带限走时层析核函数的数值模拟特征

a—单频10 Hz核函数xz切片;b—单频10 Hz核函数yz切片;c—单频20 Hz核函数xz切片;d—单频20 Hz核函数yz切片;e—单频30 Hz核函数xz切片;f—单频30 Hz核函数yz切片;g—带限核函数xz切片;h—带限核函数yz切片


图3反应的是单频和带限走时层析核函数。根据式(17),该核函数仅仅是成像域走时层析核函数的一个分支,两个分支(炮点到成像点及成像点到检波点)则构成了成像域走时层析反演核函数的完整形态,如图4。用此带有一定宽度的高斯束层析核函数替换常规射线层析核函数能更准确逼近波实际传播方式,提高层析反演精度和稳定性。

图4

图4   成像域高斯束走时层析反演核函数


成像域高斯束层析反演的实际操作流程及具体计算方法可以完全借鉴常规射线层析反演的框架和算法,包括成像剖面层位及成像道集RMO的自动拾取,矩阵求解,层位约束正则化方法等等。两者的区别仅仅是在核函数的表达与计算上,因此该方法的实用性较强。

需要说明的是,成像域走时层析的具体实现是与其走时误差的具体计算形式有关的。Xie[22,23,24]给出了基于炮偏移后按照炮索引排列的成像道集的剩余时差(RMO)所对应的核函数的表现形式,该核函数的两个分支并不对称,这与成像道集上每一个RMO代表一炮的成像走时误差有关。文中所推导的核函数与Xie给出的核函数表达形式一致,但由于是针对方位角度域成像道集,所提取的RMO仅与方位角及反射角有关[25],核函数则表现为关于剖面上反射轴法方向对称的两个分支,这一点并不会影响理论上的精度,但会使得层析的实现更加自然方便,与角度域的高斯束偏移更好地互为正反过程。

3 成像域高斯束层析方法的具体实现

与常规成像域射线层析方法一样,在成像域执行高斯束层析反演的基本思想是将成像道集的剩余时差(RMO)沿着波传播的特定路径进行反投影,得到速度模型更新量,完成速度模型的一次迭代更新。研究共成像点道集剩余时间差与剩余深度差之间的关系,并以此关系为基础建立层析方程组,进而迭代求取剩余速度。

3.1 成像域高斯束层析反问题的建立

当初始速度模型与真实模型接近时(低波数成分),对于共成像点道集满足如图5a所示的关系。图中,S是成像道集某道对应的炮点,R是对应的检波点,A是偏移速度的成像点,A'是正确速度时的成像位置,实线SAR是射线在偏移速度中的传播路径,虚线SA'R是射线在假设真速度中的传播路径,距离AB是拾取的剩余深度差,ψ是地层倾角,θϕ分别是偏移入射角和理论入射角,当传播路径比较长、速度误差比较小时,近似认为θ=ϕ,红线是共成像点位置。

图5

图5   高斯束层析传播路径几何示意

a—复杂介质传播路径;b—成像点局部传播路径


在反射点局部,可以假定慢度s是常数,又因为在远离反射点的地方两条传播路径比较接近,这样在反射点局部(图5b),走时扰动量等于两条路径之差乘以慢度s,即

ΔtRMO=2sΔL,

其中,ΔtRMO即式(16)中的走时扰动,在成像道集中表现为剩余走时差。式(18)中:

ΔL=Δzcosθcosψ,Δz=AB,

将式(19)代入式(18)得到剩余深度差与走时残差的关系式:

ΔtRMO=2scosθcosϕΔz

将式(20)求得的每个层析控制点的剩余走时差代入走时层析反演方程即可实现成像域层析反问题的建立。由此也可看出,成像域的层析反演需要在偏移剖面上自动拾取反射轴倾角ψ,以及成像道集上的剩余深度差Δz:

AΔm=ΔtRMO

A即是上述核函数计算得到的敏感度矩阵,Δm即是待反演的速度更新量。然而实际求解时,矩阵A一般是病态的,这是因为数据往往是有限的,并不足以约束所有的模型分量,破坏了反演的稳定性。因此,必须引入额外的信息对反问题进行正则化,即加入一些先验约束,使得解向所约束的方向发展。

从反演的角度,成像域层析归结为如下泛函的求极值问题:

S(m)=ΔtRMO22,

改造上式的泛函,得到如下新的误差泛函:

S(m)=ΔtRMO22+εd2m-mr22+εC12D1(m-mr)22+εC22D2(m-mr)22,

其中: εd2, εC12, εC22为较小的正实数;mr为给定的先验模型,一般取mr=0或mr=mn, mn为前一次反演结果; εd2m-mr22为阻尼项,这项的引入使得每次迭代的模型更新量较小。 εC12D1(m-mr) 22εC22D2(m-mr) 22分别约束速度模型在横向和纵向上的光滑程度,D1D2分别是对B样条基函数系数更新量在横向和纵向上的一阶差分算子。

对于更改后的误差泛函(22),每次迭代需要求解的层析方程变为:

AεdIεC1D1εC2D2Vm=VtRMO-εd(m-mr)-εC1D1(m-mr)-εC2D2(m-mr)

利用LSQR方法求解矩阵方程组(24),该方法是一种迭代的方法,可以在最小二乘意义下高效地求解大规模稀疏矩阵。至此,建立了带有正则化项的成像域高斯束层析成像反问题。

3.2 高斯束层析反演流程

成像域高斯束层析反演的基本流程同常规射线层析反演的流程基本一致(图6)。

图6

图6   高斯束层析反演流程


从流程图上可以看出,高斯束层析的输入要求从叠前深度偏移剖面上自动拾取控制点,并在这些控制点对应的成像道集上自动拾取剩余深度差。由于自动拾取精度往往受信噪比等因素的影响,而高斯束偏移本身具备更高成像信噪比的优势,因此更有利于高斯束层析的迭代反演。当然在实际操作中还需人工的检查与局部调整。

4 理论模型及实际数据应用

4.1 理论模型实验

该模型断层发育,地层高陡,速度横向变化较大(图7a)。横纵向采样点为731×550,横向采样间隔10 m,纵向采样间隔5 m。利用声波正演得到叠前炮记录(图7b),炮间距和道间距都是10 m,中心放炮,每炮361道。利用等梯度速度模型作为层析初始模型(图8a)。对初始模型进行高斯束叠前深度偏移,输出偏移剖面(图8b)及角度道集(图9),可以看出初始模型对应的偏移剖面上的各层位成像深度并不准确,绕射波也没有完全收敛。图9是CDP400位置处提取的初始模型角度道集和层析后角度道集对比,由于速度偏低,初始角度道集同相轴上翘;经过层析反演更新速度模型,偏移后角度道集得到了拉平。图10是高斯束层析经过5次迭代后得到的速度模型及其相应的高斯束偏移剖面,图11是常规射线层析经过8次迭代得到的速度模型及其高斯束偏移剖面。与常规射线层析对比可知,文中发展的高斯束层析方法能够提供更加丰富的速度信息,迭代次数也小于常规射线层析。更新后偏移提取的角度道集更加接近真实速度模型下偏移提取的角度道集,反射界面归位到正确的深度位置,绕射波收敛,断层位置聚焦更好,说明高斯束层析对复杂构造模型的速度反演精度优于常规射线层析方法。

图7

图7   某地区复杂速度模型(a)及其相应炮集记录(b)


图8

图8   初始速度模型(a)及其高斯束叠前深度偏移剖面(b)


图9

图9   CDP 400角度道集

a—初始角度道集;b—正确模型角度道集;c—高斯束层析角度道集;d—射线层析角度道集


图10

图10   高斯束层析更新速度场(a)及其层析后高斯束偏移剖面(b)


图11

图11   常规射线层析更新速度场(a)及其层析后高斯束偏移剖面(b)


成像域高斯束层析反演的实际操作流程及具体计算方法可以完全借鉴常规射线层析反演的框架和算法,包括反演过程中的一些正则化方法等。两者的区别仅仅是在核函数的表达与计算上。但从核函数的表达式可以看出,由于引进了有限频率段,利用高斯束积分计算格林函数从而计算有限频核函数的计算量将比仅计算高频近似下的射线层析核函数(射线长度)大得多。在相同的计算资源情况下,高斯束层析反演计算效率大概比射线层析慢一个数量级,但该值并不固定,且取决于不同的优化算法。后续在实用化过程中将着重进行计算效率的提升。

表1   计算机集群系统软硬件环境

硬件环境①处理器:Intel(R)Xeon(R)X5650 ②处理器主频:2.66 GHz
③物理内存:48 GB ④每个节点CPU数:2个(每个CPU8核)
软件环境① 操作系统:Red Hat Enterpise Linux 4-64 Update 5
② 并行计算环境:MPICH1.2.6
③ 编译器:INTEL C++ 、FORTRAN编译器 10.0.023 Linux 版(64位)
④ 配置英特尔 MKL9.1 (Math Kernel Library Cluster Edition)数学库
网络环境1000 M,4 GB光纤

新窗口打开| 下载CSV


表2   高斯束层析与射线层析计算效率对比

处理名称处理时间/min折算单核处理效率
高斯束层析110.30.766 M/h
射线层析9.998.46 M/h

新窗口打开| 下载CSV


该模型数据量1.1 G,在完全相同的机器环境条件下分别利用50个节点处理该模型数据,高斯束层析增加的计算量是射线层析的11倍左右。

4.2 实际资料应用

实际资料来自中国南方复杂山前带某工区。该区地形起伏剧烈(图12中的曲线代表高程),悬崖峭壁密布,山高谷深,碳酸盐岩大面积裸露,喀斯特岩溶地貌较发育,形成岭谷相间地形地貌。地震数据信噪比非常低。受地表地震地质条件及复杂构造的影响,深部构造的勘探工作困难重重,特别是在深部构造高部位通常反射能量弱、且信噪比低,波组连续性差,造成落实构造形态难度大。

图12

图12   中国南方某山前带工区某单炮实际资料


文中的研究工作是在前期处理人员进行了预处理、叠前时间偏移及深度域初始建模及偏移成像的基础上开展的,主要通过高斯束层析反演后的速度模型进行高斯束叠前深度偏移来体现本文方法的实际应用效果。

图13层析反演的速度模型上看,高斯束层析结果体现了速度模型的更多细节信息。鉴于高斯束偏移的抗噪性较强,可通过高斯束偏移进一步提高成像质量。

图13

图13   某实际数据深度域初始建模(a)常规射线层析建模(b)及高斯束层析建模(c)结果对比


图14可以看出,同样利用高斯束偏移方法,高斯束层析成像的速度模型对应的偏移结果信噪比更高,同相轴更连续,整体成像质量更高。而从图15的不同偏移算法处理可以看出,高斯束偏移结果信噪比明显高于Kirchhoff偏移结果,整体成像质量有较明显的提升,对后续的构造落实有较好的指导意义,这得益于高斯束偏移的实现可以采取类似于控制束偏移的方式,在数据合成及传播过程中进行筛选及优选[26]。通过该实际资料处理可以认识到,结合高斯束层析与高斯束偏移技术的处理方式,将能更好地适应山前带探区等复杂构造低信噪比资料,提供更高质量的成像结果。

图14

图14   实际数据的高斯束偏移结果

a—利用射线层析模型;b—利用高斯束层析模型


图15

图15   实际数据利用高斯束层析模型偏移结果

a—利用Kirchhoff偏移;b—利用高斯束偏移


5 结论

文中首先给出了高斯束偏移提取方位—反射角度道集的方法,进而从高斯束偏移角度道集出发,在波动方程的一阶Born近似和Rytov 近似下,推导给出了成像域走时层析成像方程及其核函数表达式,包括单频及带限核函数。利用带限层析核函数替代常规射线层析核函数(该核函数为常数1),可以改进层析反演精度,加快反演收敛,从而形成了基于高斯束算子的偏移成像与层析成像的联合迭代速度建模与偏移方法,相对于常规射线方法具有更高的精度和较强的实用性。

成像域速度反演方法主要利用反射数据的走时信息,对初始模型的依赖性较低。该方法仅能反演速度模型的低波数成分,因而主要用于背景速度建模,为偏移成像服务及更高精度的反演方法(如FWI)提供初始模型。文中发展的高斯束层析成像方法利用高斯束传播算子计算成像域走时层析核函数,提供了一种新的成像域波动方程线性化近似层析反演方法。基于高斯束算子的成像域走时层析方法,与高斯束偏移技术相结合,可形成具有典型特征波成像特点的、适应低信噪比数据的成像与建模工具,真正体现偏移成像与速度建模一体化处理思想。

The authors have declared that no competing interests exist.
作者已声明无竞争性利益关系。

参考文献

Vasco D W, Peterson J E, Majer E L .

Beond ray tomography:Wavepaths and Fresnel volumes

[J]. Geophysics, 1995,60(6):1790-1804.

DOI:10.1190/1.1443912      URL     [本文引用: 1]

Yomogida K .

Fresnel Zone Inversion for lateral heterogeneities in the earth

[J]. Pageoph, 1992,138(3):391-406.

DOI:10.1007/BF00876879      URL     [本文引用: 1]

刘玉柱 . 菲涅尔体地震层析成像理论与应用研究[D].上海:同济大学, 2011.

[本文引用: 3]

Bednar J B, Shin C, Pyun S .

Comparison of waveform inversion,part 2:phase approach

[J]. Geophysical Prospecting, 2007,55(4):465-475.

DOI:10.1111/gpr.2007.55.issue-4      URL     [本文引用: 1]

Bozdağ E, Trampert J, Tromp J .

Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements

[J]. Geophys. J. Int, 2011,185(2):845-870.

DOI:10.1111/gji.2011.185.issue-2      URL     [本文引用: 1]

Wu R S, Hu C, Wang B.

Nonlinear sensitivity operator and inverse thin-slab propagator for tomographic waveform inversion

[C]//Expanded Abstracts of the 84 th Annual SEG Meeting,Society of Exploration Geophysicists , 2014: 928-933.

[本文引用: 1]

Jeon S.

Full waveform inversion using the energy objective function

[C]//76 th EAGE Conference & Exhibition,Expanded Abstracts , 2014: 106:108.

[本文引用: 1]

Luo Y, Schuster G T .

Wave-equation traveltime inversion

[J]. Geophysics, 1991,56:645-653.

DOI:10.1190/1.1443081      URL     [本文引用: 1]

Xu S, Chauris H, Lambaré G,Noble M S.

Common angle image gather:A new strategy for imaging complex media

[C]//Expanded Abstracts of the 68 th Annual SEG Meeting,Society of Exploration Geophysicists , 1998: 1538-1541.

[本文引用: 1]

Xu S H, Chauris G Lambaré, Noble M S .

Common-angle migration:A strategy for imaging complex media

[J]. Geophysics, 2001,66(3):1877-1894.

DOI:10.1190/1.1487131      URL     [本文引用: 1]

蔡杰雄, 王华忠, 王立歆 .

基于三维高斯束算子解析的方位反射角道集提取技术研究

[J]. 石油物探, 2016,55(1):76-83.

[本文引用: 2]

李振春, 岳玉波, 郭朝斌 , .

高斯波束共角度保幅深度偏移

[J]. 石油地球物理勘探, 2010,45(3):360-365.

DOI:      Magsci     [本文引用: 1]

<P><FONT face=Verdana>高斯波束偏移是一种准确、灵活、高效的深度偏移方法,不但具有接近于波动方程偏移的成像能力,还可以对陡倾角地层以及各向异性介质进行成像。由于高斯波束偏移过程中含有地下射线的传播角度信息,因而可以利用此角度信息直接提取角度域共成像点道集,而不需要波动方程偏移中复杂的映射转换。本文首先介绍了高斯波束偏移的基本原理,然后实现了基于互相关成像条件的保幅高斯波束偏移,并提出了角度域共成像点道集的提取方法。模型和实际资料的试算分析表明,此方法是有效的。</FONT></P>

Sava P, Fomel S .

Angle-domain common image gathers by wavefield continuation methods

[J]. Geophysics, 2003,68(2):1065-1074.

DOI:10.1190/1.1581078      URL     [本文引用: 1]

Zhang Y, Xu S, Bleistein N , et al.

True-amplitude,angle domain,common-image gathers from one-way wave-equation migrations

[J]. Geophysics, 2007,72(1):S49-S58.

DOI:10.1190/1.2399371      URL     [本文引用: 1]

Xu S, Zhang Y, Tang B.

3D common image gathers from reverse time migration

[C]//Expanded Abstracts of the 80 th Annual SEG Meeting,Society of Exploration Geophysicists , 2010: 3257-3260.

[本文引用: 1]

Ross Hill N .

Gaussian beam migration

[J]. Geophysics, 1990,55(11):1416-1428.

DOI:10.1190/1.1442788      URL     [本文引用: 1]

Ross Hill N .

Prestack Gaussian beam depth migration

[J]. Geophysics, 2001,66(4):1240-1250.

DOI:10.1190/1.1487071      URL     [本文引用: 1]

Woodward M J .

Wave-equation tomography

[J]. Geophysics, 1992,57(1):15-26.

DOI:10.1190/1.1443179      URL    

Jeroen J, Spetzler J, Smeulders D , et al.

Validation of first-order diffraction theory for the traveltimes and amplitudes of propagating waves

[J]. Geophysics, 2006,71(6):167-177.

Spetzler G, Snieder R .

The fresnel volume and transmitted waves

[J]. Geophysics, 2004,69(3):653-663.

DOI:10.1190/1.1759451      URL     [本文引用: 1]

Liu Y Z, Dong L G, Wang Y W , et al.

Sensitivity kernels for seismic Fresnel volume tomography

[J]. Geophysics, 2009,74(5):U35-U46.

DOI:10.1190/1.3169600      URL     [本文引用: 1]

Xie X B, Yang H.

A migration velocity updating method based on the shot index common image gather and finite-frequency sensitivity kernel

[C]//Expanded Abstracts of the 77 th Annual SEG Meeting,Society of Exploration Geophysicists , 2007: 2767-2771.

[本文引用: 1]

Xie X B, Yang H .

The finite-frequency sensitivity kernel for migration residual moveout and its applications in migration velocity analysis

[J]. Geophysics, 2008a,73(6):S241-S249.

DOI:10.1190/1.2993536      URL     [本文引用: 1]

Xie X B, Yang H.

A wave-equation migration velocity analysis approach based on the finite-frequency sensitivity kernel

[C]//Expanded Abstracts of the 78 th Annual SEG Meeting,Society of Exploration Geophysicists , 2008b: 3093-3097.

[本文引用: 1]

蔡杰雄, 王华忠, 陈进 , .

基于高斯束传播算子的成像域走时层析成像方法

[J]. 地球物理学报, 2017,60(9):3539-3554.

[本文引用: 1]

刘少勇, 蔡杰雄, 王华忠 , .

山前带地震数据射线(束)叠前成像方法研究与应用

[J]. 石油物探, 2012,51(6):598-605.

[本文引用: 1]

/

京ICP备05055290号-3
版权所有 © 2021《物探与化探》编辑部
通讯地址:北京市学院路29号航遥中心 邮编:100083
电话:010-62060192;62060193 E-mail:whtbjb@sina.com