基于地震波方程的地震数据波形偏移与最小二乘波形偏移方法

陈生昌

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

地震数据偏移始于20世纪70年代初期[1-4],是当前地震数据处理中最重要的方法技术之一,也是获得地下构造三维图像的主要技术手段[5-8]。地震数据偏移的目的是利用数学处理的方法将在地表观测到的地震波场送回到产生它的位置上去,以得到反映地下产生地震散射波或反射波的非均匀体的空间位置的偏移成像波场[9-13]。散射波偏移获得散射体的位置,反射波偏移获得反射体边界的位置,即反射面的位置,因此地震数据偏移是一种构造成像。经过40多年的发展,地震数据偏移方法已从叠后发展到了叠前,从二维发展到了三维,从时间偏移发展到了深度偏移,从各向同性介质发展到了各向异性介质,从标量声波方程发展到了矢量弹性波方程[14-15],从不需迭代的常规偏移发展到了基于迭代的最小二乘偏移[16-21]。

当前的地震数据偏移和最小二乘偏移不区分散射和反射,在偏移中主要用反射的概念,在最小二乘偏移中主要应用散射的概念。偏移的目的是对反射面空间位置进行成像,但缺乏与之对应的数学物理正演方程。此外,虽然在偏移的波场传播中使用了波动方程,但由于其成像公式的原因,我们认为当前的地震数据偏移不是基于真正的波动方程方法。最小二乘偏移利用散射理论构建其对应的数学物理正演问题,但其结果是非均匀体的速度相对扰动,由于地震数据的频带有限性,该结果不能有效地刻画非均匀体边界位置,所以最小二乘偏移不是对偏移的完善而是一种线性反演。

随着油气地震勘探目标由构造油气藏向岩性油气藏和页岩油气藏的发展,希望地震数据偏移结果不仅能实现对地下非均匀体构造的成像,还能反映地下非均匀体岩性的变化,偏移结果可用于岩性解释。为了适应这种发展,本文首先简要回顾和评述了当前主要面向构造成像的地震数据偏移、地震数据偏移的照明补偿和地震数据最小二乘偏移方法,然后将地震数据偏移(包括最小二乘偏移)定义为基于地震波方程的波形线性反演。根据地下非均匀体的大小和速度变化与地震波长之间的相对关系,将非均匀体划分为散射体和反射体,产生相应的散射波和反射波,利用散射理论和高频近似分别推导出描述散射波运动的散射波方程和描述反射波运动的反射波方程,并将得到的散射波方程和反射波方程作为地震散射数据和反射数据的正演方程,最后利用线性反演方法,针对散射数据提出基于散射波方程的波形偏移与最小二乘波形偏移,针对反射数据提出基于反射波方程的波形偏移与最小二乘波形偏移。

1 地震数据偏移

在地震勘探中,地震散射波和反射波的物理传播过程可直观地表述为,地表震源激发的地震波向下传播,即入射波向下传播,入射波与地下非均匀体作用形成二次源,产生散射波和/或反射波,散射波和/或反射波向上传播到地表被检波器接收。在非均匀体处,入射波与散射波或反射波具有相同的旅行时。对于这样的地震波物理传播过程,在形式上可用BERKHOUT于1982年提出的“WRW”概念模型表示[22]。

基于上述地震波传播物理过程和BERKHOUT“WRW”概念模型,当前的地震数据偏移分两步实现,一是基于波动方程的震源波场正向传播和记录波场反向传播;二是成像,即根据入射波与散射波或反射波在非均匀体处具有相同旅行时的特征,利用正向传播得到的入射波场和反向传播得到的散射波场或反射波场实现对非均匀体空间位置的成像,即实现对非均匀体的构造成像。

为了实现偏移中的波场传播需要给定一个速度模型,通常称为偏移速度模型,也叫宏观速度模型。在该速度模型下,所谓的非均匀体是地下真实速度模型相对于宏观速度模型的差异,并认为用于偏移的地震数据(散射波数据和反射波数据)是由非均匀体产生的。因此,对于用于偏移的宏观速度模型,它不能产生可观测到的散射波数据和反射波数据,即宏观速度模型必须是光滑变化的,再者宏观速度模型必须保证入射波、散射波和反射波的旅行时准确。对于散射波和反射波,宏观速度模型控制波的传播,而非均匀体控制散射波的散射和反射波的反射。

根据上述地震数据偏移两步实现法,基于宏观速度模型v0(x)的地震数据波动方程逆时偏移计算公式如下。

1) 震源波场的正向传播:

(1)

2) 基于伴随运算的记录波场反向传播(即记录波场的逆时外推):

(2)

3) 基于时间一致性成像原理的构造成像:

(3)

或

(4)

对地震数据进行偏移处理时,我们需要一个尽可能准确的宏观速度模型v0(x),以保证在对非均匀体进行构造成像时,正向传播得到的震源波场us(x,t;xs)和反向传播得到记录波场ur(x,t;xs)有相同的旅行时,否则就不能实现对非均匀体构造的正确成像。公式(1)描述了入射波场的传播,相当于给定速度模型和震源函数的波场模拟。公式(2)以伴随运算方式描述了散射波场或反射波场的反向传播,这种伴随传播是对逆传播的近似,相对于逆传播计算效率高、稳定性好,但重建出的地下散射波场或反射波场会受到地震波照明不均匀的影响,导致后续的偏移成像结果产生成像阴影,即地震波照明弱的地方偏移成像结果较差。成像公式(3)和公式(4)将散射波场或反射波场与入射波场联系起来,并以此得到对非均匀体的成像。如果非均匀体是散射体,则成像结果是散射体空间位置的反映,如果非均匀体是反射体,则成像结果是反射体边界空间位置的反映。根据物理上的直观认识,偏移中公式(3)得到的偏移成像结果被认为是地下反射系数的反映,实际上公式(3)所表示的散射波场或反射波场与入射波场间的关系并不满足它们之间的数学物理关系,即公式(3)不是由散射波场或反射波场和入射波场所满足的地震波方程推导出来的。公式(4)是一个计算稳定的构造成像公式,在当前的偏移中被广泛应用,相对于公式(3),由公式(4)得到的偏移成像结果会受到震源地震波在地下照明不均匀和震源子波的影响,造成分辨率降低和成像阴影,即震源地震波照明弱的地方偏移成像结果较差的偏移成像阴影。

当前地震数据偏移之所以由两步组成,而且成像公式也不是由入射波和反射波(散射波)所满足数学物理方程推导出来,是因为缺乏与地震数据偏移所对应的数学物理正演方程。正因如此,当前的地震数据偏移主要是利用地震数据的相位信息,而不能正确地利用数据的波形信息。

2 地震数据偏移阴影的照明补偿

根据波场传播的Kirchhoff-Rayleigh积分Ⅱ,波场的正向传播可表示为下述矩阵方程形式:

(5)

式中:V为已知波场;L表示波场传播矩阵;U表示V经传播后得到的波场。对于由波场U经逆传播得到波场V,在最小二乘意义下有:

(6)

其中,L-1为L的最小二乘逆,即:

(7)

式中:L*称为L的伴随矩阵。公式(6)的波场传播称为波场逆传播。如果将L*作为L-1的近似,用于公式(6)的波场传播,即:

(8)

公式(8)所表示的波场传播(也称为波场的伴随传播,或逆时外推)就是上节逆时偏移中公式(2)的基于伴随运算的波场反向传播。公式(8)的伴随传播相较于公式(6)的逆传播具有计算量小和稳定性好的优势,但(L*L)-1的缺失使得经公式(8)传播得到的波场Va相对于公式(6)得到的波场V保真性差、分辨率低。在给定速度模型下,波场传播矩阵L是地下波场照明效应的反映,(L*L)是地下波场照明强度的反映。

虽然成像公式(3)相较于成像公式(4)存在计算稳定性方面的问题,但它不会受到震源地震波照明不均匀的影响,而且还能消除数据中的震源子波效应,有利于得到高分辨率的偏移结果。

在上节偏移中,如果利用从检波点方向地震波对地下的照明强度对应用公式(2)重建的地下波场进行补偿,同时应用成像公式(3)进行成像或利用震源方向地震波对地下的照明强度对成像(公式(4))的成像结果进行补偿,则我们就可以实现对地震数据偏移成像阴影的照明补偿[23]。经照明补偿得到的偏移结果在深度方向上具有更好的均衡性,分辨率也会得到一定程度的提高。偏移照明补偿对于盐下构造、高速推覆体下构造、深层和超深层构造成像具有重要意义。

3 地震数据最小二乘偏移

常规偏移中记录波场的地下重建是借助波场传播算子的伴随完成的,得到的重建波场会受到地下速度变化导致的地震波照明不均匀的影响,再者,成像公式(3)也不是基于控制地震波的数学物理方程得到的。虽然利用上节的照明补偿可以在一定程度上改进照明不均匀导致的偏移结果的成像阴影,但为了获得不受地震波照明不均匀影响的高保真而且可以反映地下岩性变化的偏移结果,人们将地震数据偏移定义为一个基于波动方程的线性反演问题,并严格求解。为了构建该反问题,我们首先要在波动方程的基础上建立其对应的正演问题。

假定速度模型可表示为:

(9)

其中,v0(x)为变化缓慢的背景模型,也可视为前述地震数据偏移中的宏观速度模型;δv(x)为变化剧烈的扰动模型,对应的地震波可表示为:

(10)

其中,us(x,t)为背景模型下震源激发产生的波场;δu(x,t)为扰动模型产生的扰动波场。根据散射理论和Born近似,有:

(11)

公式(11)就是当前地震数据最小二乘偏移的正演方程,式中的m(x)就是待求的最小二乘偏移结果。在该方程中扰动波场δu(x,t;xs)与m(x)呈线性关系。利用Green函数,方程(11)可写为下述积分方程:

(12)

式中“*”代表时间褶积。利用最小二乘法求解积分方程(公式(12)),得到m1(y,t;xs)的解估计,有:

(13)

式中:M1为对应m1(y,t;xs)的矩阵;G-1为对应Green函数g(x,t;y)的矩阵G的最小二乘逆矩阵,有G-1=(G*G)-1G*,G*为G的伴随矩阵;ΔU为对应δu(x,t;xs)的矩阵。公式(13)的运算即是波场的逆传播。利用反演得到的m1(y,t;xs)和背景模型下的模拟波场∂2us(y,t;xs)/∂t2,可以实现对m(y)的求解,即:

(14)

公式(13)和公式(14)虽然在理论上给出了m(y)的结果,但存在一些数值计算方面的难题:一是最小二乘逆矩阵G-1中逆矩阵(G*G)-1的大规模计算问题;二是公式(14)右端中的波场相除运算会出现计算不稳定问题。针对这些计算方面的问题,常用迭代的方法求解m(y)。用G的伴随矩阵G*代替式(13)中的G-1,即:

(15)

公式(15)的运算是波场的伴随传播,可以利用波场逆时外推的方法实现,如上述公式(2)的计算。同时将公式(14)改写为:

(16)

应用Landweber迭代方法[24],可以得到下述基于迭代的最小二乘偏移。

1) 初始化。

2) 迭代过程。

式中α为步长因子,可用不同方法求取,如简单的线性搜索法,也可用(G*G)-1的近似值,如检波点方向的地震波单向照明强度。

上述最小二乘偏移实质上是对速度相对扰动的线性反演,不是对速度扰动体边界空间位置和边界反射系数的成像。它只是在迭代中利用了波场传播算子的伴随近似波场传播算子的逆,即利用波场逆时外推的方式进行地下记录波场的重建,这一步运算与逆时偏移是相同的。最小二乘偏移是在地震数据所满足的地震波方程(即公式(12))的基础上,利用最小二乘反演方法推导得到,因此不需要利用偏移方法中所谓的成像原理,这也是公式(14)与公式(3)的区别。我们认为,最小二乘偏移不是对当前偏移方法的改进,而是利用地震数据的散射理论表达重新定义的一种地震数据线性反演方法,这是因为一方面当前由波场外推和成像两个步骤构成的逆时偏移方法缺乏与之对应的数学物理正演方程,另一方面带限地震数据的最小二乘偏移结果不能清楚地反映速度扰动体的边界位置。

上述地震数据的两步偏移法及其成像公式是一种直观的方法,没有考虑宏观速度模型下散射波和反射波的运动方程,在地震数据偏移和最小二乘偏移中也不区分散射波和反射波。偏移时将散射波视为反射波,统一用反射波的偏移方法进行偏移。在最小二乘偏移中将反射波视为散射波,统一用散射理论构建最小二乘偏移方法。当前偏移是对产生反射波的反射体边界空间位置进行成像,得到的结果理论上是一个δ函数;而最小二乘偏移是对产生散射波的非均匀体的速度扰动进行反演,得到的结果理论上是一个阶跃函数。由于实际地震数据是缺失低频和高频成分的带限数据,因此所得到的最小二乘偏移结果不是一个阶跃函数。我们认为,实际地震数据最小二乘偏移结果可以准确地反映局部散射体的空间位置,但不能准确反映非局部反射层的边界空间位置。

4 地震散射波方程与反射波方程

为了使得到的地震数据偏移结果不仅能实现对地下非均匀体的构造成像,还能反映地下非均匀体岩性的变化并用于岩性解释,我们需要重新定义地震数据偏移,将其定义为基于地震数据正演方程的波形线性反演,因而必须要有合适的描述地震波运动(如传播和散射或反射等)的数学物理方程。

对于产生散射波和/或反射波的非均匀体,根据其体积大小,以地震波长为尺度,将空间尺寸与地震波长相当的非均匀体视为产生散射波的散射体,而将空间尺寸比地震波长大的非均匀体视为产生反射波的反射体。此外,我们还可以根据速度扰动δv(x)的变化,以地震波长为尺度将变化快的δv(x)视为产生散射波的散射体,将变化相对较慢的δv(x)视为产生反射波的反射体。上述划分相当于散射体在高频近似条件下退化为反射体。在我们的研究中,散射体对于地震波视为空间中的一个点,反射体对于地震波视为在空间上具有一定延续度的几何体,因此散射体产生的散射波没有方向性,而反射体产生的反射波具有方向性。

根据散射理论,对于相对于地震波长具有局部性和变化快的δv(x)产生的一次散射波,在Born近似下,可用下述方程描述:

(17)

其中,av(x)=[2δv(x)]/[v0(x)],称为速度相对扰动;us(x,t;xs)为背景速度模型v0(x)下的背景波场,有:

(18)

公式(17)称为入射波作用下散射体产生的散射波的散射波方程,其右端项称为产生散射波的二次源。

由于us(x,t;xs)在非源处满足齐次波动方程,因此公式(17)可写为:

(19)

利用Green函数将公式(19)写为积分形式:

(20)

如果地震波长相对于非均匀体的尺寸比较小或/和δv(x)的空间变化相对于地震波长尺度缓慢,即在av(y)相对地震波长尺度具有变化缓慢的高频近似条件下,则可将非均匀体视为反射体,相应地产生具有方向性的反射波,这种具有方向性的反射波可近似为局部平面波。相应的一次散射波方程(20)退化为一次反射波方程:

(21)

(22)

将方程(22)变换到时间域,并写为微分形式,可得到反射波动方程(23),即:

(23)

式中:vr(x,σ)称为与角度σ有关的速度相对扰动,有vr(x,σ)=av(x)cosσ。方程(23)为高频近似下基于速度相对扰动的反射波方程。

在地震勘探中,反射体产生的反射波直观上可视为由反射面产生,而且反射面上的反射率也是一个常用的概念。反射率不仅能刻画反射面的空间位置,还能反映反射面上速度扰动的变化。我们将反射波方程(23)中的速度相对扰动vr(x,σ)沿入射波传播方向的方向导数定义为反射率。在vr(x,σ)相对地震波长满足高频近似条件的情形下,有下述反射率定义式:

(24)

利用反射率定义公式(24),可得到高频近似下基于反射率变化的反射波方程:

(25)

5 地震散射数据的波形偏移

上节推导出的散射波方程(17)和背景速度下的入射波方程(18)是进行散射波模拟和反演成像的基础方程。所谓地震散射数据的波形偏移是在给定背景速度模型v0(x)、震源函数f(t)和散射地震数据δu(xg,t;xs)的条件下,根据散射数据的正演方程(17),利用反演的方法实现对散射体空间位置的成像。

利用Green函数,散射数据的正演方程(17)在频率域可写为积分形式:

(26)

(27)

(28)

(29)

由于我们的目的是得到av(y)的空间位置,根据文献[25],可得到下述求取av(y)空间位置的散射地震数据波形偏移计算公式。

1) 地下入射波场的计算(对应公式(27))。

(30)

2) 基于波场逆时传播的地下散射波场近似重建(对应积分方程(28)的近似求解)。

(31)

3) 得到对av(x)的近似反演(对应公式(29)的近似计算)。

(32)

由于我们将产生散射波的散射体近似为一个点源体,因此公式(32)所得到的av(x)的近似反演结果也同样能实现对其空间位置的成像[26]。

上述散射数据偏移计算公式完全是根据散射数据的正演表达式,利用反演方法得到,不需要应用常规偏移中必须用的成像原理(公式)。由于利用了散射数据的正演方程,所得到的偏移方法可以可靠地利用散射数据的波形信息,因此这样的偏移方法称为波形偏移方法。

6 地震散射数据的最小二乘波形偏移

上节提出的由散射数据波形偏移方法得到的av(x)虽然能有效地实现对散射体相对速度扰动空间位置的成像,但在有关相对速度扰动av(x)数值的保真度和分辨率方面与上文的最小二乘偏移方法一样会受地下速度变化导致的地震波照明不均匀的影响和公式(32)对公式(29)所做的近似的影响。为了消除这些影响,得到高分辨和高保真度的相对速度扰动av(x)反演结果,我们依照第4节的最小二乘偏移方法并结合文献[27]提出的时间二阶积分波场的全波形反演方法,提出了下述基于迭代的地震数据最小二乘波形偏移方法。

给定背景速度模型v0(x)、震源函数f(t)和散射地震数据δu(xg,t;xs),可利用下述迭代格式的最小二乘波形偏移求解av(x)。

1) 初始化。

2) 迭代过程。

计算数据残差:

数据残差的逆时外推:

(33)

(34)

k=k+1

式中:δui(xg,t;xs)表示共炮道集地震散射数据的时间二阶积分。修正公式(34)中的α为步长因子,α可用类似于上面最小二乘偏移中的方法求取。

在迭代过程中,如果将前述的照明补偿加进来则可以有效提高迭代的收敛速度。

7 地震反射数据的波形偏移

相对于散射波数据,反射波数据是地震勘探中最常见的数据,反射数据的偏移与反演方法是地震数据处理中重要的方法技术,也是获取地下三维构造图像的重要手段。上面推导出的基于速度相对扰动和反射率变化的两个反射波方程(23)、方程(25)和背景速度下的入射波方程(18)是进行反射波模拟和反演成像的基础方程。所谓地震反射数据的波形偏移就是在给定背景速度模型v0(x)、震源函数f(t)和反射地震数据ur(xg,t;xs)的条件下,根据反射数据的正演方程(反射波方程(23)或方程(25)),利用反演的方法实现对反射体速度相对扰动和反射率的近似反演。

7.1 速度相对扰动的波形偏移

基于上面推导出的基于速度相对扰动的反射波方程(23)和入射波方程(18),采用类似于推导地震散射数据波形偏移的方法,我们可以推导出地震反射数据的速度相对扰动波形偏移(近似反演)计算公式。

1) 地下入射波场的计算。

(35)

2) 基于波场逆时传播的地下反射波场近似重建。

(36)

式中:ur(xg,t;xs)为共炮道集地震反射数据。

3) 得到对速度相对扰动vr(x,σ)的近似反演结果。

(37)

速度的相对扰动是一个阶跃函数,因此在对其进行反演成像时数据的低频成分非常重要,为了凸显数据中的低频成分,根据公式(29)和文献[26]提出的时间二阶积分波场的全波形反演方法,我们将公式(36)和公式(37)分别改写公式(38)和公式(39)。

(38)

(39)

本节所提出的速度相对扰动波形偏移方法可以用于对地下反射层岩性变化进行成像[30]。

7.2 反射率的波形偏移

基于反射率变化的反射波方程(25)和入射波方程(18),采用类似于推导上述地震反射数据的速度相对扰动波形偏移方法,我们可以推导出地震反射数据的反射率波形偏移计算公式[25]。

1) 地下入射波场的计算。

(40)

2) 基于波场逆时传播的地下反射波场近似重建。

(41)

式中:ur(xg,t;xs)为共炮道集地震反射数据。

3) 得到关于反射率r(x,σ)的近似反演结果。

(42)

反射率是速度相对扰动的方向导数,理论上在速度分界面处是一个δ函数,一个沿速度分界面的分布函数,因此即使对于带限地震反射数据,本节所得到的r(x,σ)的近似反演结果也能有效地刻画反射面的空间位置。

上述推导出的波形偏移方法不同于当前基于Claerbout成像原理的波动方程逆时偏移方法,它是一种基于反射波动方程的近似反演方法,不需要应用所谓的偏移成像原理及其成像公式,能够在波动方程意义下真实地应用地震数据的波形信息,是一种真正的波动方程偏移方法。

8 地震反射数据的最小二乘波形偏移

上节提出的反射数据波形偏移方法虽然可以实现对速度相对扰动vr(x,σ)和反射率r(x,σ)的近似反演成像,但其结果的保真度和分辨率会受到因地下速度变化引起的地震波照明不均匀和反演近似的影响。为了消除这些影响,得到对速度相对扰动vr(x,σ)和反射率r(x,σ)的高分辨和高保真反演结果,依照地震散射数据最小二乘波形偏移方法,提出地震反射数据最小二乘波形偏移方法。

8.1 速度相对扰动的最小二乘波形偏移

根据前述地震反射数据速度相对扰动波形偏移计算公式,我们提出基于迭代的地震反射数据速度相对扰动的最小二乘波形偏移。

给定背景速度模型v0(x)、震源函数f(t)和反射地震数据ur(xg,t;xs),可利用下述迭代格式的最小二乘波形偏移求解速度相对扰动vr(x,σ)。

1) 初始化。

2) 迭代过程。

计算数据残差:

(43)

数据残差时间二阶积分的逆时外推:

(44)

(45)

k=k+1

在迭代过程中,如果将照明补偿加进来则可以有效提高迭代的收敛速度。

本节所提出的速度相对扰动最小二乘波形偏移实际上是一种对地下反射层岩性变化的线性反演,所得到的反演结果相对于速度相对扰动波形偏移的结果具有更好的保真性和更高的分辨率。

8.2 反射率的最小二乘波形偏移

根据前述地震反射数据反射率波形偏移计算公式,我们提出基于迭代的地震反射数据反射率的最小二乘波形偏移[31]。

给定背景速度模型v0(x)、震源函数f(t)和反射地震数据ur(xg,t;xs),可利用下述迭代格式的最小二乘波形偏移求解反射率r(x,σ)。

1) 初始化。

令:rk(x,σ)=0,k=0,计算震源波场:

2) 迭代过程。

计算数据残差:

(46)

数据残差的逆时外推:

(47)

对rk(x,σ)的修正:

rk+1(x,σ)=rk(x,σ)+αv0(x)·

(48)

k=k+1

在迭代过程中,如果将照明补偿加进来则可以有效提高迭代的收敛速度。

本节所提出的反射率最小二乘波形偏移实际上是一种对地下反射层界面反射率变化的线性反演,所得到的反演结果相对于反射率波形偏移的结果具有更好的保真性和更高的分辨率。

9 结论

在给定偏移的宏观速度模型下,本文根据地下非均匀体大小和其速度变化与地震波长之间的相对关系,将非均匀体划分为散射体和反射体,相应地产生散射波和反射波。在散射理论的基础上,首先得到描述散射波传播与散射的散射波方程;然后再利用高频近似,得到两种描述反射波的反射波方程,即基于速度相对扰动的反射波方程和基于反射面上反射率的反射波方程。将得到的散射波方程和反射波方程作为地震散射数据和反射数据的正演方程,利用线性反演方法,分别提出地震散射数据的波形偏移与最小二乘波形偏移和地震反射数据的波形偏移与最小二乘波形偏移。地震反射数据的波形偏移与最小二乘波形偏移包括了目标分别为速度相对扰动和反射率两种不同形式。基于地震数据正演方程所提出的地震数据波形偏移与最小二乘波形偏移是对当前地震数据偏移和最小二乘偏移方法的理论完善,其偏移结果相对于当前的偏移与最小二乘偏移的结果具有更明确的数学物理意义。