用于地震反射数据偏移的反射波方程

陈生昌

(浙江大学地球科学学院,浙江杭州310027)

地震勘探中的反射波是最主要的地震波类型,反射数据是地震观测数据的主要成分,反射数据的偏移成像是地震数据处理的重要方法技术,也是获取地下三维构造图像的重要手段[1]。但随着油气勘探开发目标的日趋复杂化和精细化(如构造岩性油气藏和非常规油气藏的勘探开发)以及宽方位、宽频带、高密度数据采集方法技术的应用,人们希望偏移方法能更好地利用数据的波形信息,获得具有更高分辨率和能更准确地反映地下岩性变化的偏移结果,不仅满足复杂构造精细成像的要求,还能开展岩性分析研究,以满足地震数据岩性处理解释的需要。

当前的地震反射数据波动方程逆时偏移成像需要给定一个光滑的宏观模型,也称为偏移模型,利用偏移模型下的波动方程分别对震源波场和记录波场进行顺时外推和逆时外推,然后再利用基于时间一致性成像原理[2-6]建立的成像公式对外推得到的震源波场和记录波场进行成像,得到反映反射面空间位置和反射面反射系数的偏移成像结果。本文将偏移定义为地震反射数据的波形线性反演问题,因为光滑偏移模型下的波动方程只能描述入射波与反射波的传播,而不能描述反射波的反射,光滑偏移模型下的波动方程不能作为表达地震反射数据的正演方程。也正是缺乏针对反射数据的数学物理正演方程(即描述反射波运动(传播与反射)的反射波方程),致使偏移中的波场外推与波场成像之间相互独立,缺乏内在联系。当前偏移方法中的成像公式不是由偏移模型下描述入射波与反射波之间关系的数学物理方程推导出来的,仅在旅行时方面满足偏移速度模型下的波动方程,也未考虑入射波场和反射波场之间的波形变化,仅适合平面波入射到平反射面及其产生的反射平面波的特殊情况。由于缺乏与之相应的正演方程,因此我们认为当前的地震反射数据波动方程逆时偏移成像方法不是一种真正意义上完全基于地震反射波方程的方法,不能有效地利用数据的波形信息。

当前的地震反射数据最小二乘偏移将偏移定义为线性反演问题,并将利用地震波传播的散射理论和Born近似构建的散射方程作为反射数据的正演方程[7],该正演方程中的模型参数是偏移模型的相对扰动(偏移模型的相对扰动在理论上是一个阶跃函数),而不是反射面的空间位置或反射面反射系数(反射面的空间位置或反射系数在理论上是一个δ函数)。最小二乘偏移中反射数据正演方程的不恰当以及地震反射数据的带限性质致使最小二乘偏移不能给出有关反射面的空间位置并得到高保真的反射系数以及高分辨的偏移成像结果[8],因此我们认为当前反射数据最小二乘偏移只是一种求解偏移速度相对扰动的线性反演,而不是对偏移的改进。

波动方程可以根据地下模型参数的变化情况来描述全部波型(如直达波、透射波、散射波、绕射波和反射波等)的产生与传播,在地震波场的传播规律、特征研究和地震波的全波场反演中得到了广泛应用[9-11]。但在上述所谓的波动方程偏移和最小二乘偏移中,不区分散射与反射(在偏移中将散射和反射均当反射考虑,而在最小二乘偏移中将散射和反射均当散射考虑)是不恰当的。散射是一种局部性的波动行为,而反射在一定程度上可视为散射的叠加,因此散射波没有特定的方向性,而反射波具有特定的方向性,散射和反射应该区别对待。对于单纯的散射数据或反射数据反演/偏移,应该有其特定的散射波方程或反射波方程,即应该在给定的偏移模型下,分别构建描述散射波运动的散射波方程和描述反射波运动的反射波方程[12]。

为了推导与真振幅Kirchhoff积分偏移(偏移/反演)成像方法相对应的与角度有关的反射系数的地震反射数据正演计算公式,BLEISTEIN等[8,13-14]将描述绕射波场传播的惠更斯原理和Kirchhoff积分公式应用于反射波场,得到了基于反射系数而不是速度扰动的Kirchhoff积分形式的地震反射数据正演计算公式。该反射数据正演公式推导中利用了有关反射面上的反射波场等于入射波场乘以反射系数的Kirchhoff近似,因而仅适合入射平面波和平反射面的情况。

为了解决反射数据偏移和最小二乘偏移方法存在的正演方程问题,给有效利用地震反射数据波形信息的偏移方法研究提供相应的正演方程,以满足地下复杂构造油气藏、岩性油气藏和非常规油气藏勘探开发对高精度偏移成像的要求,必须研究偏移模型下适合高精度偏移成像的反射数据正演方程。陈生昌等[15-19]在地震数据正演方程及其偏移成像方法研究方面进行过一些初步的探索研究,本文在前期工作的基础上,对适合反射数据高精度偏移成像的正演方程进行系统研究。在给定光滑偏移模型的前提下,依据地震波长与地下非均匀体的大小及其物理性质空间变化快慢之间的相对关系,在地震波长相对于地下非均匀体的大小及其物理性质空间变化满足高频近似的条件下,将比波长大的地下非均匀体近似为反射体,产生具有局部方向性的反射波;从散射理论中的波动方程出发,基于上述反射体和反射波概念,推导出基于物性参数相对变化的反射波方程;再针对反射面反射率在地震反射中的广泛应用,利用定义物性参数相对变化沿入射波传播方向的方向导数为反射体上反射面的反射率,推导出基于反射率变化的反射波方程;最后根据地下介质的不同,分别推导出宏观偏移模型下基于物性参数相对变化和反射率变化的标量反射波方程、声反射波方程和弹性反射波方程。

1 波动方程与地震波的散射和反射

在应力与应变间满足线性近似的假设下,地震震源激发的地震波在地下的运动可以用非线性微分算子方程描述:

(1)

式中:L为与地下物性参数模型m(也称为弹性参数模型)有关的二阶偏微分算子;u为地震波在地下运动状态的状态变量,也称为地震波场;f为激发地震波的震源函数。方程(1)也称为非齐次波动方程,或辐射方程。如果方程(1)的右端项等于0(无源项),则得到齐次波动方程,即

(2)

如果地下物性参数模型m仅为速度v(x),则方程(1)为非齐次标量波动方程,即

(3)

如果地下物性参数模型m为速度v(x)和密度ρ(x)的函数,则方程(1)为非齐次声波动方程,即

(4)

如果地下模型为各向同性弹性介质,其物性参数模型m为拉梅系数λ(x)、μ(x)和密度ρ(x)的表达式,则方程(1)为非齐次弹性波动方程[20],即

(5)

式中:u(x,t;xs)为矢量波场;f(t,xs)为矢量震源函数;L是一个3×3偏微分算子,即

波动方程(1)和波动方程(2)能描述震源和m变化产生的各种波,如直达波、散射波、反射波、透射波和转换波(对于弹性波动方程)等,因此也称为通用波动方程,是地震波场模拟、地震数据全波形反演和波场外推的基础方程。如果地下模型为光滑的宏观模型,则波动方程(1)和波动方程(2)就不能描述散射波和反射波。为了研究光滑宏观模型下地震反射波的波场模拟、反演和成像方法,特别是反演和成像方法,如反射波形反演、反射波阻抗反演、反射数据的偏移成像(包括最小二乘偏移成像)方法等,我们需要描述光滑宏观模型下的地下反射波运动的反射波方程,即反射数据的正演方程。

地震反射波在地下运动的物理过程直观上可表述为:地表震源激发的入射波场向下传播时,入射波场与反射体作用形成产生反射波场的二次震源,反射波场向上传播到地表被检波器接收。在这个直观的表述中,反射波的传播和反射相互独立,但在描述地震波场运动状态的通用波动方程(1)中,地震反射波的传播与反射交织在一起。这是因为方程(1)中的模型m不仅控制了地震波传播,也控制了地震波的反射。

为了清楚地描述地震反射波的反射与传播状态,根据上述对地震反射波的直观认识,我们需要分别描述反射波场传播和反射的数学物理方程,也就是要在模型m中将传播和反射部分分开,得到控制反射波传播的传播模型mp和控制反射波反射的反射模型mr,在形式上将m表示为m=mp+mr。这样的反射波数学物理描述对于反射数据的偏移成像方法研究很有帮助。

为了认识传播模型mp和反射模型mr,我们将地震波传播的散射理论应用于方程(1)。令m=mb+δm,其中,mb表示光滑宏观(背景)模型,δm为体现m中快速变化的模型扰动;相应的波场为u=ub+δu,其中,δu为波场扰动;ub表示背景波场。有

(6)

将模型和波场扰动表达式代入方程(1),并结合方程(6),有:

(7)

式中:L(mb)为背景模型下的波动算子;δL(mb,δm)为背景模型mb和扰动模型δm产生的扰动算子,有δL(mb,δm)=L(mb+δm)-L(mb);u为与模型m对应的全波波场。方程(7)也称为波动方程(1)的扰动形式,它在数学上与方程(1)一样可以描述模型m变化产生的各类波型,如散射波和反射波等。在δm≪mb的一阶Born近似下,方程(7)可近似为描述一次扰动波场的方程:

(8)

式中:u1代表一次扰动波场。下面针对δm的尺寸及其空间变化快慢与地震波长之间的相对关系来研究方程(7)和方程(8)所描述的地震波的性质。

若方程(7)和方程(8)中的δm对应的不均匀体的空间尺寸与地震波长相当,或者不均匀体物理性质的空间变化相对于地震波长尺度为快速变化,我们认为这样的不均匀体是散射体,其右端的源项就相当于点源,相应地产生没有特定方向性(在各向同性介质中)的散射波,方程(7)和方程(8)就分别称为全散射波方程和一次散射波方程,有

若方程(7)和方程(8)中的δm所对应的不均匀体的空间尺寸比地震波长大,或者不均匀体物理性质的空间变化相对于地震波长尺度为缓慢变化甚至常数,即在地震波长相对于地下非均匀体的大小及其物理性质的空间变化满足高频近似的条件下,我们可以将这样的不均匀体近似为反射体,其右端源项就相当于一个在空间上具有一定延续度的局部平面波源,相应地产生具有方向性的局部平面波性质的反射波,因此方程(7)和方程(8)就分别退化为全反射波方程和一次反射波方程,有

在光学中,所谓的反射是平面波在无穷大平面(镜面)上的反射,是由镜面两侧介质的速度变化产生的[21]。入射波与反射波之间的运动学特征满足Snell定律,而它们之间的动力学特征满足反射波的振幅等于镜面反射系数乘以入射波振幅的规律。在地震勘探中,地下介质主要由块状体和层状体组成。对于相对于地震波长具有很大空间延续度的层状介质产生的地震波,可以仿照光学中的镜面反射波进行研究,而对于块状体产生的地震波,就必须根据空间尺寸与地震波长之间的相对关系来决定它是散射体产生的散射波还是反射体产生的反射波。本文所讨论的地震反射是高频近似条件下对地震散射的近似,即在地震波长比较短的情况下δm的空间分布和其物理性质的空间变化相对于地震波而言可视为反射体,它与入射波场作用产生反射波场。

根据上述有关反射体和反射波的认识和定义,我们在高频近似条件下,假定背景模型mb光滑变化,即方程(6)中的波场ub不包含反射波场,则地面观测到的地震反射都是由δm产生的。与mb对应的波动算子L(mb)控制和描述了地震波场的传播,我们称光滑的背景模型mb为地震波的传播模型。与δm有关的扰动算子δLr(mb,δm,σ)控制和描述了地震波的反射,我们称δm和mb共同组成了反射波的反射模型。在高频近似条件下,方程(11)和方程(12)不仅描述了反射波的反射与传播,还描述了入射波与反射波之间的数学物理关系,是光滑宏观模型下描述反射波运动的抽象反射波方程。相对于虚点源产生的散射场,空间上具有一定延续度的局部平面波虚源产生的反射波场可视为积分波场。借助Green函数将方程(11)和方程(12)分别写为下述积分形式的反射数据正演方程,即

式中:Ω为产生反射波的局部平面波虚源的分布空间;g为背景模型下的Green函数。

2 地震反射波方程

2.1 标量反射波方程

假定方程(1)中的模型m为速度模型v(x),并有v(x)=vb(x)+δv(x),其中,vb(x)为光滑背景速度模型,δv(x)为速度扰动。如果速度扰动体δv(x)的体积大小及其速度空间变化满足散射体条件,则由方程(6)所对应的Green函数和齐次波动方程(2),得到全散射和一次散射方程(9)、方程(10)的积分形式:

假定地震波长相对于速度扰动体δv(x)(或av(y))的大小及其空间变化快慢满足高频近似条件,即av(y)可近似为反射体,相应的散射波退化为反射波。则方程(15)和方程(16)分别退化为与方程(13)和方程(14)对应的反射波方程,即

(19)

公式(17)和公式(18)中的波场和Green函数具有局部平面波性质。将方程(17)和方程(18)写为微分形式:

方程(17)和方程(18)分别称为基于速度相对扰动的全反射数据和一次反射数据的积分形式正演方程。

在地震勘探中,反射面的反射率不仅可以表征反射面的空间位置,还能反映反射面上的物性变化。为了得到基于反射率的反射波方程,我们定义反射率为速度相对扰动av(x)沿入射波传播方向I的方向导数。在δv(x)的空间变化满足高频近似(速度相对扰动av(x)的空间变化相对于地震波的波长尺度可视为缓慢变化甚至常数)的条件下,即在反射波近似条件下,由速度相对扰动av(x)沿入射波传播方向I的方向导数所定义的反射率r(x,σ)在波数域的定义公式为:

(22)

根据反射率的定义,反射率r(x,σ)在空间域有如下定义公式:

(23)

公式(22)和公式(23)定义的反射率r(x,σ)是地下速度模型空间变化的反映,也和角度σ有关,是一个标量。不同于根据反射波场振幅与入射波场振幅之比定义的无量纲的反射系数,本文定义的r(x,σ)具有长度倒数的量纲,它是高频近似条件下具有局部平面波性质的入射波作用于局部反射面的反射率,即反射体边界局部切平面上的反射率。虽然BLEISTEIN等[8]和STOLT等[22]分别提出了类似的利用方向导数定义的反射率,但BLEISTEIN等定义的反射率是从反演成像分辨率的角度考虑的,STOLT等定义的反射率是点反射率不是局部反射面反射率,且两者的方向导数都是反射面的法向导数而不是入射波传播方向的方向导数。

下面用一个两层一面的分层速度模型来说明公式(22)和公式(23)定义的反射率与反射系数之间的关系。假定模型上层速度为v1,下层速度为v2,速度分界面深度为z1,I为垂直向下的入射波传播方向,即σ=0,如图1所示。

图1 两层一面的分层速度模型

取vb=v1,则速度扰动δv(z)=(v2-v1)H(z-z1)。其中,H(z)为阶跃函数。假定入射波的传播方向I为垂直方向,根据反射率定义公式(23),图1 中速度模型的反射率为:

(24)

它是一个δ函数。在z=z1速度分界面上,对于法向入射地震波的反射率r,有r|z=z1≈4(v2-v1)/(v2+v1)=4R,其中,R为地震波法向入射时速度分界面的反射系数。

由上述分析可知,公式(22)和公式(23)定义的反射率与反射系数虽然有不同的物理量纲,但在数值上与反射系数成正比关系。

利用反射率定义公式(22),我们可以得到基于反射率变化的全反射波和一次反射波的反射波方程:

方程(26)所表示的一次反射波方程的右端项不同于文献[12]中基于SYMES等[23]的反射波方程而提出的一次反射波方程的右端项，文献[12]中的一次反射波方程在当前的反射波全波形反演、最小二乘偏移中得到了广泛应用，它的右端项是由无量纲的反射系数和入射波场的乘积组成，其左、右两端的物理量纲不匹配。因此，我们认为文献[12]中的一次反射波方程是错误的。

由方程(25)和方程(26)可以得到基于反射率变化的全反射数据和一次反射数据的积分形式正演方程:

式中:xg代表检波点位置坐标。将上述一次反射数据正演方程变换到频率域,有

(29)

公式(29)与BLEISTEIN等[8,24]利用Kirchhoff近似推导出的反射数据正演公式相似,但BLEISTEIN等给出的公式是基于无量纲的反射系数推导的,而公式(29)是基于具有长度倒数量纲的反射率推导的。

2.2 声反射波方程

为了考虑地下密度变化对地震波运动的影响,假定方程(1)中的模型m为包含速度v(x)和密度ρ(x)的声学介质,并有v(x)=vb(x)+δv(x),ρ(x)=ρb(x)+δρ(x),其中,vb(x)和ρb(x)为光滑背景模型,δv(x)和δρ(x)为扰动模型。假定地震波长相对于δv(x)和δρ(x)的空间体积大小和空间变化快慢满足高频近似条件,即与δv(x)和δρ(x)对应的非均匀体可近似为反射体,则根据方程(11)和方程(12),进行与文献[16]相类似的推导可得到下述声波反射波方程:

(32)

方程(30)和方程(31)中的av+aρ(1+cosσ)是一个与σ有关的声波阻抗相对扰动,令Ir(x,σ)=av+aρ(1+cosσ)。如果取σ=0,即沿反射面法向入射,则有Ir(x,σ=0)=av+2aρ=2δv(x)/vb(x)+2δρ(x)/ρb(x)≈2δI(x)/Ib(x),即小角度反射主要受阻抗变化的控制。其中,Ib(x)=vb(x)ρb(x),δI(x)=δ[v(x)ρ(x)]=ρ(x)δv(x)+v(x)δρ(x)。如果取σ=π,则有Ir(x,σ=π)=av=2δv(x)/vb(x),即大角度反射主要受速度变化的控制。因此,在高频近似和vb(x),ρb(x)为光滑变化的条件下,我们分别称方程(30)和方程(31)为基于声波阻抗相对扰动的全反射波和一次反射波的声波反射波方程。

由方程(30)和方程(31)可以得到基于声波阻抗相对扰动的全反射数据和一次反射数据的积分形式正演方程:

(33)

(34)

定义背景模型vb(x)、ρb(x)中入射波传播方向的局部波数为ki,有ki=ω/vb(x)。为了得到基于反射率的反射波动方程,我们定义反射率为声波阻抗相对扰动Ir(x,σ)沿入射波传播方向I的方向导数。在δv(x)、δρ(x)的空间变化满足高频近似(即声波阻抗相对扰动Ir(x,σ)的空间变化相对于地震波的波长尺度可视为缓慢变化甚至常数)的条件下,声波阻抗相对扰动Ir(x,σ)沿入射波传播方向I的方向导数所定义的反射率r(x,σ)在波数域有如下定义公式:

(35)

根据反射率的定义,反射率r(x,σ)在空间域有如下定义公式:

(36)

公式(35)和公式(36)所定义的反射率是一个与反射开角σ有关的变量,是介质声波阻抗空间变化的反映,具有长度倒数的量纲。

图2为两层一面分层声阻抗模型。对于法向入射地震波在z=z1声阻抗分界面上的反射率r,有:r|z=z1≈4(ρ2v2-ρ1v1)/(ρ2v2+ρ1v1)=4Ra。其中,Ra为地震波法向入射时声阻抗分界面的反射系数。因此,公式(35)和公式(36)定义的反射率与声阻抗分界面上的反射系数也具有正比关系。

利用公式(35)给出的反射率定义和2.1节的推导,我们可以得到高频近似条件下基于反射率变化的全反射波和一次反射波的声反射波方程:

由方程(37)和方程(38)可以得到基于反射率变化的全反射数据和一次反射数据的积分形式正演方程:

(39)

(40)

图2 两层一面分层声阻抗模型

2.3 弹性反射波方程

(43)

ΔLr为弹性波反射算子,有

(44)

在高频近似和αb(x)、βb(x)、ρb(x)为光滑变化的条件下,我们分别称方程(41)和方程(42)为基于弹性参数相对扰动的全反射波和一次反射波的弹性反射波方程。

为了将弹性矢量波场u=(ux,uy,uz)转换为纵波场和横波场ups=(up,ush,usv),得到P波和S波型的反射波方程,对方程(41)和方程(42)应用Helmholtz变换和对S波的旋转[22]得到:

(47)

公式(47)中各元素代表了不同类型入射波与不同类型反射波间的反射作用。在高频近似条件下,这些元素在频率域有[22,25]:

在法向入射时,σ=0,有

(48)

其中,

(49)

(49)式中的各元素称为弹性波阻抗相对扰动,有

将(48)式代入方程(45)和方程(46),有

(50)

(51)

方程(50)和方程(51)分别称为高频近似条件下基于弹性波阻抗相对扰动的P波和S波型的全反射波和一次反射波的弹性反射波方程。

由方程(50)和方程(51)可以得到基于弹性波阻抗相对扰动的P波和S波型的全反射数据和一次反射数据的积分形式正演方程:

(54)

对于SH波入射、SH波反射的反射率,有

(55)

对于SV波入射、SV波反射的反射率,有

(56)

对于P波入射、SV波反射的反射率,有

(57)

对于SV波入射、P波反射的反射率,有

(58)

图3为两层一面分层弹性阻抗模型,对于法向入射的不同波型地震波,在z=z1弹性波阻抗分界面上的反射率r,分别有

(59)

式中:Rpp为P波法向入射时弹性阻抗分界面的反射系数。

(60)

式中:Rshsh为SH波法向入射时弹性阻抗分界面的反射系数。

(61)

图3 两层一面分层弹性阻抗模型

式中:Rsvsv为SV波法向入射时弹性阻抗分界面的反射系数。

由上述的简单弹性分层模型可知,我们定义的反射率与弹性阻抗分界面上的反射系数也同样具有正比关系。

(62)

其中,

利用(62)式可得到基于反射率变化的全反射波和一次反射波的弹性反射波方程,即

(63)

(64)

由方程(63)和方程(64)可以得到基于反射率变化的P波和S波型的弹性波全反射数据和一次反射数据的积分形式正演方程:

在上述推导出的标量波、声波和弹性波的反射波方程中,基于波阻抗相对变化的反射波方程的右端项是由波场的时间二阶导数与波阻抗相对变化的乘积组成,而基于反射率的反射波动方程的右端项是由波场的时间一阶导数与反射率的乘积组成。与平面波无穷大平边界的无量纲反射系数不同,本文基于地震波长相对于地下非均匀体大小满足高频近似和局部平面波而定义的反射率是一个具有长度倒数量纲的变量,是真实介质与背景介质之间模型参数相对扰动沿入射波传播方向的方向导数,它不仅依赖于背景介质,也与入射波传播方向有关。Zoeppritz方程描述了平面波无穷大平边界条件下反射系数与入射方向和介质模型参数之间的定量关系,本文推导的基于反射率的反射波方程建立了高频近似条件下反射波与入射波和反射率之间的数学物理关系,但目前我们还不能写出倾斜入射情况下反射率与模型参数相对扰动和入射波传播方向之间的具体函数关系。

3 结论

光滑宏观模型下的波动方程只能描述反射波的传播,而不能描述反射波的反射,致使当前地震反射数据的波动方程偏移缺少与之对应的反射数据正演方程。根据地下非均匀体的大小及其物理性质空间变化快慢与地震波长之间的相对关系,把非均匀体划分为产生散射波的散射体和产生反射波的反射体,本文视反射体和反射波为散射体和散射波的高频近似,并认为反射波是一种具有方向性的局部平面波。利用散射理论的波动方程和本文所得到的有关反射体、反射波的认识与定义,在给定的光滑宏观模型下,首先推导出基于地下模型参数相对扰动的反射波方程,然后将地下模型参数相对扰动沿入射波传播方向的方向导数定义为地下反射率分布,再推导出基于反射率变化的反射波方程。根据地下介质的不同,分别推导出了光滑背景模型下基于模型参数相对变化和反射率变化的标量反射波方程、声反射波方程和弹性反射波方程。本文在高频近似条件下定义的反射率反映了地下模型参数相对扰动沿入射波传播方向的变化,与基于波场振幅变化定义的无量纲的反射系数不同,本文的反射率具有长度倒数的量纲,它是高频近似条件下,局部平面波的入射波作用于反射体边界的局部切平面的反射率。从简单的分层介质模型可知,本文定义的反射率与反射系数在数值上具有正比关系。本文推导出的反射波方程可作为反射数据的正演方程,是开展地震反射数据反演和偏移方法研究的基础方程,特别是开展可有效利用反射数据波形信息的偏移方法研究的基础方程。