基于地震数据线性正演表达的波形成像(二):地震数据的波形成像

陈生昌,鲁方正,刘亚楠,周华敏

(1.浙江大学地球科学学院,浙江杭州310027;2.长江水利委员会长江科学院,湖北武汉430010)

在波形成像论文的第一篇中,对于给定的地震波运动学光滑模型,将扰动法应用于波动方程,得到了描述非均匀体扰动波场的波动方程。依据地震波长与非均匀体大小或非均匀体物理性质空间变化特征尺度之间的相对大小关系,将非均匀体划分为在空间上具有有限延续度可产生散射波的局部散射体或空间上具有一定延续度可产生反射波的反射体,并将相应的扰动波场的波动方程转化为散射波动方程和反射波动方程。对于散射波,利用小扰动假定和Born近似,得到基于散射体物性参数相对扰动的散射数据线性正演表达式;对于反射波,利用波阻抗相对扰动的小值假定、一次波近似和高频近似(地震波长相对于非均匀体大小或其物理性质空间变化特征尺度为较小量),得到基于反射体波阻抗相对扰动的反射数据线性正演表达式。为了描述反射体边界对反射波的作用和反射体边界上物理量在反射数据偏移成像中的广泛应用,在高频近似条件(光滑模型的空间变化相对于地震波长尺度可视为缓慢变化甚至不变化)下定义反射体波阻抗相对扰动沿入射波传播方向的方向导数为反射体边界局部反射系数,推导出了基于反射体局部反射系数的反射数据线性正演表达式,并且给出了标量波、声波和弹性波方程下散射数据和反射数据的线性正演表达式。

本文是波形成像论文的第二篇,将第一篇推导出的地震一次散射数据和一次反射数据线性正演表达式作为地震数据线性反演的正演方程,利用线性反演方法提出包括散射数据波形成像和反射数据波形成像的地震数据波形成像方法理论。首先根据第一部分论述中提出的地震数据线性正演表达式,提出用于散射数据波形成像和反射数据波形成像的地震波传播通用方程、散射通用方程和反射通用方程,再利用线性反演理论建立地震数据波形成像的基本方法理论框架,最后在波形成像的线性反演中利用地震波传播算子的伴随算子和最小二乘逆算子分别提出波形成像的波形偏移方法和最小二乘波形偏移方法。对于地震反射数据在两种不同线性正演表达下的波形成像,首先给出面向地层波阻抗相对扰动的岩性成像和反射面局部反射系数的构造成像,然后针对具体的标量波散射数据、声波反射数据和弹性波反射数据分别给出它们的波形偏移和最小二乘波形偏移具体计算公式,最后进行合成数据的波形成像数值试验,试验结果表明波形偏移结果相对于逆时偏移结果具有更加准确的相位和更高的分辨率。

1 地震数据波形成像的方法理论框架

地震数据的正演表达式是地震数据反演成像的基础,不同的地震数据正演表达式会导致不同的反演目标和反演方法。当前地震数据波动方程偏移方法缺少与之对应的严谨的地震数据正演方程,使得偏移方法中描述波场传播的方程与描述波场成像的成像公式不匹配,皆因波场传播方程和成像公式是由不同条件下的波动方程导出的。偏移成像公式应该和产生散射波或反射波的虚源有关,即应该和散射波或反射波与入射波之间所满足的数学物理方程有关。本文将描述散射波或反射波与入射波之间关系的数学物理方程称之为散射(波)方程或反射(波)方程。在第一部分论述中关于地震数据线性正演表达的基础上,首先建立与地震数据波形成像有关的地震波传播、散射、反射的通用方程,然后在此基础上利用反演理论构建地震数据波形成像的方法理论框架。

1.1 用于散射数据波形成像的地震波传播与散射通用方程

当前的地震数据波动方程偏移成像不区分散射和反射,将无穷大平边界平面波的镜面反射方程(第一部分论述中的(6)式)作为描述散射波与入射波之间关系的散射方程,用于散射数据偏移成像。第一部分论述中将相对于地震波长尺度可近似为点源的虚源产生的地震波定义为散射波,因此我们认为在散射波偏移成像中采用第一部分论述中的(6)式描述散射不合适,它不能准确地描述局部的散射现象和散射波与入射波之间的相位变化。

根据第一部分论述中在光滑模型下建立的散射数据线性正演表达公式,我们建立用于描述地震波传播与一次散射的通用波动方程,即:

L(mb)ui(x,t;xs)=f(x,t)

(1)

L(mb)us(x,t;xs)=-δLs(mb,δm)ui(x,t;xs)

(2)

式中:mb为已知的光滑模型;L(mb)为光滑模型下的通用波动算子;ui(x,t;xs)为光滑模型下的背景波场或入射波场;f(x,t)为震源函数;xs为震源位置;us(x,t;xs)为一次散射波场;δLs(mb,δm)为散射算子,它的具体表达式与δm/mb和波动算子L(mb)的具体形式有关;δm为模型扰动,其空间尺寸大小或其物理性质空间变化特征尺度相对于地震波长为较小量,满足散射体条件。波动方程(1)描述了光滑模型下的入射波场及其传播,它所对应的Green函数也是光滑模型下散射波的波场传播算子。波动方程(2)不仅描述了散射波场与入射波场之间应满足的数学物理关系,也描述了散射波的产生(即产生散射波的虚源)与散射波的传播。我们将(1)式(或其对应的齐次方程)和(2)式统称为用于散射数据波形成像的地震波传播与散射通用方程,它们是由相同条件下的波动方程推导而来的,因此二者相互匹配,由此可以进一步推导出一次散射波数据的线性正演表达,该方程也是利用一次散射波对模型相对扰动δm/mb或模型扰动δm进行反演成像的基础方程。

1.2 用于反射数据波形成像的地震波传播与反射通用方程

在当前的地震数据波动方程逆时偏移方法中,将无穷大平边界平面波的镜面反射方程(第一部分论述中的(6)式)作为一般情况下反射波的反射方程,并用于反射数据偏移成像。这样的镜面反射方程不能准确地描述非平边界、非平面波情况下的反射以及反射波与入射波之间的相位差异。根据第一部分论述中在光滑模型下利用高频近似和局部平面波建立的反射数据线性正演表达公式,可以得到地震波传播与反射的通用方程。对于基于反射体波阻抗相对扰动产生的一次反射,可以得到:

L(mb)ui(x,t;xs)=f(x,t)

(3)

L(mb)ur(x,t;xs)=-δLr(mb,Ir,σ)ui(x,t;xs)

(4)

式中:ur(x,t;xs)表示一次反射波场;δLr(mb,Ir,σ)表示随角度σ变化的反射算子,它的表达式与波阻抗相对扰动Ir和波动算子L(mb)的关系式有关;σ表示入射方向与反射方向之间的开角。(3)式描述了光滑模型下的入射波场及其传播规律,所对应的Green函数也是光滑模型下反射波的波场传播算子。(4)式不仅反映了反射波场与入射波场之间应满足的数学物理关系,也描述了基于相对扰动Ir的反射波产生(即产生反射波的虚源)与传播规律。我们将(3)式(或其对应的齐次方程)和(4)式称为基于波阻抗相对扰动Ir的反射数据波形成像地震波传播与反射通用方程,它们均由相同条件下的波动方程推导而来,因此二者是相互匹配的,由此可以推导出基于波阻抗相对扰动Ir的一次反射波数据的线性正演表达公式,该方程也是利用一次反射波数据对波阻抗相对扰动Ir进行反演成像的基础方程。

考虑基于反射体边界局部反射系数产生的反射,推导得到如下的地震波传播与反射的通用方程:

L(mb)ui(x,t;xs)=f(x,t)

(5)

L(mb)ur(x,t;xs)=-δLr(mb,r,σ)ui(x,t;xs)

(6)

式中:r表示反射体边界局部反射系数;δLr(mb,r,σ)表示基于r随角度σ变化的反射算子,它的具体表达式与r和波动算子L(mb)的具体形式有关。(6)式不仅反映反射波场与入射波场之间应满足的数学物理关系,也描述了基于局部反射系数的反射波的产生与传播规律。我们将(5)式(或其对应的齐次方程)和(6)式称为利用基于局部反射系数的反射数据波形成像的地震波传播与反射通用方程,它们也是由相同条件下的波动方程推导而来,因此二者是相互匹配的,它们的结合可作为基于反射体边界局部反射系数产生的一次反射波数据的线性正演表达,由此可以推导出利用一次反射波对反射体边界局部反射系数r进行反演成像的基础方程。

由上述分析可知,对于反射波和散射波而言,它们具有相同的入射方程和传播方程(或传播算子),却具有不同的散射方程和反射方程。

1.3 地震数据波形成像的方法理论框架

根据上文中的地震波传播通用方程、散射通用方程和反射通用方程,再利用(1)式所对应的Green函数,可得到下述通用的散射数据或反射数据线性积分表达式:

(7)

式中:u(xg,t;xs)代表散射数据或反射数据;xg代表检波点坐标;Ω表示散射波或反射波虚源的空间分布范围;g(xg,x,t)为Green函数,代表从Ω内x到xg的波场传播算子;*t代表时间褶积;s(x,t;xs)为(2)式、(4)式或(6)式右端的散射波或反射波虚源项,该虚源项包含了光滑模型下的入射波场ui(x,t;xs)、产生散射或反射的模型参数和散射算子或反射算子。

利用(7)式左端项的地震数据u(xg,t;xs),估计右端积分中的虚源项s(x,t;xs),即为求解第一类线性积分方程,也称为反源问题(inverse source problem)[1]。将(7)式代表的积分方程离散为矩阵方程形式,如下:

d=Lm

(8)

式中:d代表共炮点道集地震数据(即一次散射数据或一次反射数据)矩阵;L代表由波场传播算子组成的波场传播矩阵;m表示待求的地下虚源矩阵。

对于矩阵方程(8),有下述的最小二乘解:

m=(L*L)-1L*d=L-gd

(9)

式中:L-g代表波场传播矩阵L的最小二乘逆矩阵;L*代表矩阵L的伴随矩阵。如果(L*L)为酉矩阵或忽略(L*L)-1的作用,即用L*作为L的逆矩阵,则可得到矩阵方程(8)的近似解ma,即:

ma=L*d

(10)

对于大规模反演成像而言,利用最小二乘方法求解(9)式存在逆矩阵计算困难和不稳定性问题,但得到的解保真性和分辨率高。(10)式的求解方法称为反投影法或反传播法,具有计算量小、稳定的特点,所得到的解保真性和分辨率不及最小二乘法所得结果。对于求解(9)式中存在的逆矩阵计算困难和解不稳定问题,均可由迭代方法解决[2]。

求解(9)式或(10)式均可以得到产生散射波或反射波的虚源s(x,t;xs)的解估计。为了得到产生散射波的散射体物性参数扰动或产生反射波的反射体波阻抗相对扰动或反射体边界局部反射系数的成像结果,我们还需要利用波动方程(1)计算虚源s(x,t;xs)中的入射波场ui(x,t;xs),然后再从s(x,t;xs)的解估计中消除ui(x,t;xs),同时还要考虑散射算子或反射算子的具体形式,并将其从s(x,t;xs)的解估计中消除。这些操作类似于常规逆时偏移中的波场成像运算,但在常规逆时偏移方法中不区分散射算子和反射算子之间的差异,并利用无穷大平界面平面波的镜面反射算子R(σ)I近似本文中的散射算子和反射算子,R(σ)为随角度σ变化的反射系数[3],I为单位算子。因此,偏移成像公式应该与散射或反射虚源中的散射算子或反射算子(或散射或反射虚源的表达式)有关。

根据上述内容可得到如下的地震数据波形成像方法理论框架:

1) 利用观测到的地震数据(散射数据或反射数据)求解反源问题,得到对地下散射波或反射波虚源的解估计;

2) 利用光滑模型下的波动方程(1)求解地下入射波场;

3) 利用地下入射波场、散射算子或反射算子的具体形式,根据散射虚源或反射虚源的解估计得到散射体物性参数扰动或反射体波阻抗相对扰动或反射体边界局部反射系数的成像。

在波形成像方法理论框架的反源问题求解中,利用(10)式得到虚源的解估计,对应的地震数据波形成像方法称为地震数据波形偏移方法。在波形成像方法理论框架的反源问题求解中,利用(9)式得到虚源的解估计,则对应的地震数据波形成像方法称为地震数据最小二乘波形偏移方法。

利用上述最小二乘波形偏移方法,在理论上可以得到散射体物性参数扰动、反射体波阻抗相对扰动、反射体边界局部反射系数的高分辨率和高保真成像结果。但由于实际地震数据的频带局限性(特别是缺少低频成分)和观测范围的有限性(缺失大偏移距观测),导致得到的最小二乘波形偏移结果通常缺少低波数成分,仅仅是一种对散射体物性参数扰动或反射体波阻抗相对扰动的视反演,因此需要在最小二乘波形偏移中加入有关散射体物性参数扰动或反射体波阻抗相对扰动低波数成分的先验信息,才能实现对散射体物性参数扰动值或反射体波阻抗相对扰动值的高分辨率和高保真反演。

2 标量波散射数据的波形成像

根据上节给出的地震数据波形成像方法理论框架,对于不同波动方程所描述的地震数据线性正演表达式,给出不同的地震波传播方程和散射或反射方程,进而可以得到不同的地震数据波形成像公式。在给定的光滑速度模型下,对于标量波动方程所描述的散射线性数据正演表达式,存在下述具体的地震波传播方程和散射方程,即:

(11)

(12)

(13)

式中:g(xg,x,t)表示光滑速度模型下方程(11)的Green函数。

将上节的地震数据波形成像方法理论框架及其波形偏移的伴随求解方法以及(10)式应用于求解线性积分方程(13),再经过一些简单的推导可得到下述求取速度相对扰动av(x)近似估计的标量波散射数据波形偏移计算公式。

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

(14)

2) 基于边值条件和波场逆时传播的地下散射虚源近似反演:

(15)

3) 对av(x)的近似估计:

(16)

上述公式中:us(xg,t;xs)表示共炮点道集散射数据;T表示地震数据的最大记录时间。由于我们将产生散射波的散射体相对于地震波长近似为一个点源体,因此利用(16)式得到的av(x)的近似估计可以实现对散射体空间位置的成像[4]。

(16)式相较于常规逆时偏移成像公式的主要区别是多了对震源波场的时间二阶导数,因此(16)式得到的成像结果相对于常规逆时偏移结果存在180°的相位差异(即反极性),而且分辨率也会更高。

为了得到具有高保真性和高分辨率的av(x)反演结果,根据上节的地震数据波形偏移方法理论框架及最小二乘波形偏移的最小二乘公式,再结合上述的散射数据波形偏移计算公式,可通过迭代方式求解av(x)的散射数据最小二乘波形偏移计算公式。

① 初始化。

入射波场计算:利用公式(14)计算入射波场。

② 迭代过程

计算数据残差:

(17)

数据残差的逆时外推:

(18)

(19)

k=k+1,迭代结束。

修正公式(19)中的αk为第k次迭代步长因子(下同),可采用简单的线性搜索法获得[5],也可通过震源方向和检波点方向的地震波双向照明强度近似求解[6]得到。

3 声波反射数据的波形成像

根据前文所述的反射数据波形成像的地震波传播与反射通用方程,针对速度、密度变化的声波动方程所描述的地震反射数据,在给定光滑速度vb(x)、密度ρb(x)的模型下,对于基于波阻抗相对扰动产生的反射,有下述具体的地震波传播方程和反射方程:

(20)

(21)

(22)

式中:g(xg,x,t)表示光滑速度、密度模型下声波方程(20)的Green函数。

如果不考虑密度的变化,即δρ(x)=0,则与角度有关的声波阻抗相对扰动简化为与角度无关的速度相对扰动,即Ir(x,σ)=av(x),方程(20)至方程(22)简化为下述描述基于速度相对扰动产生的反射所对应的方程:

(20-1)

(21-1)

(22-1)

对于基于反射体边界局部反射系数产生的反射,可以得到具体的地震波传播方程、反射方程和反射数据积分表达式,如下:

(23)

(24)

(25)

式中:r(x,σ)表示反射体边界的局部反射系数,是声波阻抗相对扰动Ir(x,σ)沿入射波传播方向的方向导数,r(x,σ)=(iω)/[vb(x)]Ir(x,σ)=(iω)/[vb(x)][av+aρ(1+cosσ)]。方程(24)的右端项是产生反射波的虚源,包含局部反射系数、反射算子和入射波场,其中反射算子由光滑速度、密度模型和时间一阶导数组成,即-1/[ρb(x)vb(x)]∂/∂t。

对于方程(23)至方程(25),如果不考虑密度的变化,则δρ(x)=0,可简化为如下方程:

(23-1)

(24-1)

(25-1)

式中:r(x)代表反射体边界的局部反射系数,是速度相对扰动沿入射波传播方向的方向导数,r(x)=iω/vb(x)av(x)。

将上文的地震数据波形成像方法理论框架及波形偏移的伴随求解方法公式(10)分别应用于线性积分方程(22)和方程(25),再经过一些简单的推导可以得到求取声波阻抗相对扰动Ir(x,σ)和反射体边界局部反射系数r(x,σ)近似估计的声波反射数据波形偏移计算公式,如下。

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

(26)

2) 基于边值条件和波场逆时传播的地下反射波虚源近似反演:

(27)

3)对波阻抗相对扰动Ir(x,σ)和局部反射系数r(x,σ)的近似估计:

(28)

(29)

式中:ur(xg,t;xs)为共炮点道集反射数据。(28)式和(29)式相较于常规逆时偏移成像公式的主要区别在于多了时间二阶导数和时间一阶导数,因此它们得到的成像结果相对于常规逆时偏移结果分别存在180°和90°的相位差异,此外分辨率有所提高。

为了得到σ角度域的Ir(x,σ)和r(x,σ),在进行(28)式和(29)式右边的运算前需要将波场ur(x,t;xs)和ui(x,t;xs)转换到相应的角度域[7-8]。由(28)式可知,利用角度域的共成像点道集Ir(x,σ)数据和Ir(x,σ)的表达式可以进行类似于基于Zoeppritz方程的振幅随入射角变化(AVA)速度扰动、密度扰动和波阻抗扰动分析和反演。但是波阻抗相对扰动Ir(x,σ)不同于Zoeppritz方程中随入射角变化的反射系数,在地震波长相对于反射体的大小或其物理性质空间变化的特征尺度满足高频近似条件下,其表达式是由波动方程导出的,而Zoeppritz方程是由波动方程在平边界和平面波条件下导出的;再者Ir(x,σ)中的角度σ是根据局部切平面和局部平面波的概念定义的,是x点处的局部角度,不同于Zoeppritz方程中假定无穷大平边界和平面波条件下定义的角度。由于本文定义的局部反射系数与入射波的局部波数有关(即与频率有关),因此我们目前还不能推导出r(x,σ)与速度变化、密度变化、波阻抗变化和角度之间的具体函数关系,因此利用波形偏移得到的角度域共成像点道集r(x,σ)难以进行AVA速度分析、密度分析和反演。

如果不考虑密度的变化,即δρ(x)=0,将上文中的地震数据波形成像方法理论框架及其波形偏移的伴随求解方法公式(10)应用于线性积分方程(22-1)和方程(25-1),再经过一些简单的推导可分别得到常密度声波反射数据波形偏移计算公式、速度相对扰动av(x)和反射体边界局部反射系数r(x)近似估计公式。

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

(26-1)

2) 基于边值条件和波场逆时传播的地下反射波虚源近似反演:

(27-1)

3) 对速度相对扰动av(x)和局部反射系数r(x)的近似估计:

(28-1)

(29-1)

由于标量波方程(或常密度声波方程)所描述的反射波与入射角无关[9],所以标量波散射数据的波形成像公式(14)式至(16)式和常密度声波反射数据的波形成像速度扰动公式(26-1)式至(28-1)式在形式上完全一致。

1) 初始化。

入射波场计算:利用(26)式计算入射波场。

2) 迭代过程。

计算数据残差:

(30)

数据残差的逆时外推:

(31)

(32)

k=k+1,迭代结束。

如果不考虑密度的变化,即δρ(x)=0,可以得到如下的速度相对扰动av(x)的反射数据最小二乘波形偏移。

1) 初始化。

入射波场计算:利用(26-1)式计算入射波场。

2) 迭代过程。

计算数据残差:

(33)

数据残差的逆时外推:

(34)

(35)

k=k+1,迭代结束。

1) 初始化。

入射波场计算:利用(26)式计算入射波场。

2) 迭代过程。

计算数据残差:

(36)

数据残差的逆时外推:

(37)

(38)

k=k+1,迭代结束。

如果不考虑密度的变化,即δρ(x)=0,得到下述局部反射系数下r(x)反射数据最小二乘波形偏移成像结果,r(x)是速度相对扰动沿入射波传播方向的方向导数。

1) 初始化。

令k=0,rk(x)=r0(x)=0。

入射波场计算:利用(26-1)式计算入射波场。

2) 迭代过程。

计算数据残差:

(39)

数据残差的逆时外推:

(40)

对rk(x)的修正:

(41)

k=k+1,迭代结束。

4 弹性波反射数据的波形成像

对于纵、横波速度和密度变化的各向同性弹性波动方程所描述的地震反射数据,假定地震波在给定的纵、横波速度和密度光滑模型αb(x),βb(x),ρb(x)下不发生散射、反射和波型转换,利用反射数据波形成像的地震波传播与反射通用方程,可以得到如下的弹性入射波方程和一次反射波方程:

Lbui(x,t;xs)=f(t,xs)

(42)

Lbur(x,t;xs)=ΔLrui(x,t;xs)

(43)

(44)

式中:g(xg,x,t)为光滑模型下弹性波方程(42)所对应的矢量Green函数。

由于笛卡尔坐标系下的弹性波场中P波和S波是耦合的,在弹性波反射数据的逆时偏移成像中,为了避免弹性波场的分量对分量成像产生的交叉噪声以及得到易于物理解释的成像结果,需要在成像前将弹性波场解耦为P/S波型波场[10-15]。本文将笛卡尔坐标系下的弹性波场转换成以P、SH和SV波极化方向为坐标系下的P/S波型弹性波场,实现弹性波场的P波和S波解耦。

将光滑模型下的弹性波传播算子Lb和反射算子ΔLr转换成与P/S波型对应的算子,根据波形成像论文的第一篇,本文(43)式可表示为:

(45)

(46)

(47)

(48)

(49)

(50)

(51)

(52)

(53)

P波入射、P波反射的波场局部反射系数为:

(54)

SH波入射、SH波反射的波场局部反射系数为:

(55)

SV波入射、SV波反射的波场局部反射系数为:

(56)

P波入射、SV波反射的波场局部反射系数为:

(57)

SV波入射、P波反射的波场局部反射系数为:

(58)

由上述的rmn(m,n=p,sh,sv)表达式可知,rsvp,rpsv是关于σ的奇函数,所以rsvp,rpsv在σ=0的两侧也会出现极性反转。

利用上述的P波、SH波和SV波极化方向单位矢量和P波、SH波和SV波标量波场(本文的P波、SH波和SV波分量),还可以得到P、SH和SV波的矢量波场[12-13,17]。对入射波场ui(x,t;xs)而言,P波的矢量波场为:

(59)

SH波的矢量波场为:

(60)

SV波的矢量波场为:

(61)

对于反射波场ur(x,t;xs),P波的矢量波场为:

(62)

SH波的矢量波场为:

(63)

SV波的矢量波场为:

(64)

1) 地下入射弹性波场的计算:

Lbui(x,t;xs)=f(t,xs)

(65)

2) 基于边值条件和波场逆时传播的地下反射波虚源近似反演:

(66)

式中:ur(xg,t;xs)代表弹性反射波的炮集数据。

(67)

(68)

(69)

(70)

1) 初始化。

入射波场计算:利用公式(65)计算入射波场。

反射波场逆时外推:

(71)

2) 迭代过程。

计算分量位波场残差:

(72)

(73)

k=k+1,迭代结束。

1) 初始化。

入射波场计算:利用公式(65)计算入射波场。

反射波场逆时外推:利用公式(71)进行反射波场逆时外推。

2) 迭代过程。

计算分量位波场残差:

(74)

(75)

k=k+1,迭代结束。

(76)

(77)

(78)

(79)

由上述的标量波、声波和弹性波数据的波形成像计算公式可知,波形成像方法中的波形偏移、最小二乘波形偏移的计算过程与逆时偏移、最小二乘逆时偏移的计算过程相同,计算量相当。对于地震数据的波形成像处理而言,如果波形偏移成像处理的目标是散射体,则将地震数据视为散射数据进行本文提出的散射数据波形偏移,可以得到散射体的波形偏移成像结果;如果偏移成像处理的目标是反射体(或反射面),则将地震数据视为反射数据进行本文提出的反射数据波形偏移,可以得到反射体(或反射面)的波形偏移成像结果。

上述提出的反射数据波形成像不仅能实现反射界面的准确成像,还能得到地层波阻抗相对变化(岩性)的成像,综合应用波阻抗相对扰动的角度域共成像点道集数据和波阻抗相对扰动的表达式还可以进行类似于基于Zoeppritz方程的AVA分析和反演。

5 数值试验

为了验证上述波形成像方法理论,我们首先采用不考虑密度变化的标量波动方程合成地震数据进行标量波地震数据的波形成像试验,然后再采用考虑密度变化的声波动方程合成地震数据进行声波地震数据的波形成像试验。

5.1 标量波地震数据的波形成像试验

5.1.1 Sigsbee2A模型试验

采用SEG提供的Sigsbee2A模型数据验证本文提出的散射数据波形偏移方法和反射数据局部反射系数波形偏移方法效果。图1a为Sigsbee2A的地层模型,该模型中不仅包含高速盐体、精细的地层,还包含许多局部高速散射体和底部的水平高速反射面,十分适合用于验证散射数据波形偏移和反射数据波形偏移方法的效果。图1b为用于波形偏移试验的光滑速度模型,模型大小为1067×1201个网格点,水平和垂直方向的网格间距分别为22.86m和7.62m。该模拟数据的观测系统为海上拖缆观测系统,共有500炮,炮间距为45.72m,拖缆长度为7932.42m,检波点间距为22.86m,各炮的道数不均,最少为57道,最多为348道,时间采样间隔为8ms,时间采样点数为1500。

图1 Sigsbee2A模型

分别应用常规逆时偏移方法(波形成像第一部分论述中的(1)式至(3)式)、本文提出的散射数据波形偏移方法计算公式(本文(14)式至(16)式)和反射数据局部反射系数波形偏移方法计算公式(本文(26-1)式至(29-1)式)对Sigbee2A模型数据进行试验。图2 为应用常规逆时偏移方法、反射数据波形偏移方法和散射数据波形偏移方法得到的Sigsbee2A模型偏移结果。

图2 Sigbees2A模型偏移结果

由图2中的偏移结果可以看出,3种偏移方法都取得了不错的成像结果。对比图2中的3种偏移结果可知,对于高速散射体的成像,散射数据波形偏移方法得到的散射体偏移结果要优于常规偏移方法和反射波形偏移方法的偏移成像结果;对于反射面的成像,反射数据波形偏移方法得到的反射面的偏移成像结果要优于常规偏移方法和散射波形偏移方法的偏移成像结果。

为了仔细对比常规逆时偏移、反射数据波形偏移和散射数据波形偏移对散射体和平反射面的成像结果的差异,我们分别对图1所示的红色方框中的高速散射体和蓝色方框中的高速平反射面的成像结果进行对比。图3为图1红色方框中的高速散射体模型和分别利用逆时偏移方法、反射数据波形偏移方法、散射数据波形偏移方法得到的成像结果。由图3可以清楚地看出,散射数据波形偏移方法(其偏移结果见图3d)对高速散射体成像(成像位置和零相位)准确;逆时偏移方法(其偏移结果见图3b)对高速散射体成像位置准确,但出现了180°的相位差(反极性),将高速散射体成像为低速散射体;反射数据波形偏移方法(其偏移结果见图3c)对高速散射体成像出现了90°的相位差,致使成像位置变浅,有1/4波长的深度误差,高速散射体模型位于成像结果的波峰至波谷过渡带的中间位置。由图3还可以看出,散射数据波形偏移结果的分辨率最高,而逆时偏移结果的分辨率最低。图3的试验结果与本文的理论(见波形成像第一部分论述中逆时偏移方法中的(3)式和本文的(16)式(29-1)式)相符,波形成像第一部分论述中的(3)式与本文的(16)式相差一个时间二阶导数,即产生180°的相位差,本文的(29-1)式与(16)式相差一个负号和时间一阶导数,即产生-90°的相位差。本文的(16)式和(29-1)式中入射波场的时间导数有助于提高偏移结果的分辨率。

图4为图1蓝色方框中的高速平反射面模型和分别利用逆时偏移方法、反射数据波形偏移方法、散射数据波形偏移方法得到的成像结果。由图4可以清楚地看出,采用反射数据波形偏移方法(其偏移结果见图4c)得到的高速平反射面成像结果(成像位置和零相位)准确;采用逆时偏移方法(其偏移结果见)图4b 得到的高速平反射面成像结果存在90°的相位差,致使成像位置变浅,并存在1/4波长的深度误差,高速平反射面模型位于成像结果的波峰至波谷过渡带的中间位置;采用散射数据波形偏移方法(其偏移结果见图4d)对高速平反射面成像出现了-90°的相位差,致使成像位置变深,并存在1/4波长的深度误差,高速平反射面模型位于成像结果的波谷至波峰过渡带的中间位置。图4中的试验结果与本文的理论相符,波形成像第一部分论述中的(3)式与本文的(28-1)式相差一个一阶时间导数,即产生90°的相位差,本文的(16)式与(29-1)式相差一个负号和时间一阶导数,即产生-90°的相位差。

图4 高速平反射面模型(a)、逆时偏移(RTM)方法(b)、反射数据波形偏移(RWFM)方法(c)和散射数据波形偏移(SWFM)方法(d)成像结果

对比上述偏移成像结果可知,将本文提出的散射数据波形偏移方法应用于散射体成像、反射数据波形偏移方法应用于反射面成像,均可以取得理想的效果,它们的成像结果相对于逆时偏移结果不仅分辨率高,而且位置准确、相位(极性)正确。散射数据波形偏移方法有利于散射体的准确成像,反射数据波形偏移方法有利于反射面的准确成像。

5.1.2 反射数据的局部反射系数最小二乘波形偏移试验

为了验证本文提出的反射数据局部反射系数最小二乘波形偏移方法,我们截取Sigsbee2A速度模型(图5a)的一部分进行试验。对于图5a的速度模型采用标量波动方程的有限差分方法合成地震数据。本次试验合成了300炮地震数据,炮点范围是0～6km,炮间距20m,为中间放炮两边接收,每炮由601道接收,道间距10m,时间采样间隔是0.5ms,时间采样点数为3000,震源为主频18Hz的Ricker子波。由于图5a的速度模型为层状结构,因此合成地震数据主要为反射数据。图5b为用于最小二乘波形偏移试验的光滑速度模型,它是对图5a的速度模型进行光滑处理得到的。

在试验中,我们首先应用反射数据局部反射系数波形偏移方法对合成的炮道集数据进行局部反射系数波形偏移,图5c为得到的波形偏移结果。由图5c可以看出,将局部反射系数波形偏移方法应用于该模型可获得很好的成像结果。图5d为应用本文提出的反射数据局部反射系数最小二乘波形偏移方法得到的偏移结果,共进行了35次迭代。对比图5c和图5d可以看出,图5d的最小二乘波形偏移结果的分辨率相较于图5c的波形偏移结果明显得到了提高,并且层位面与断层面的成像更清晰。由图5e所示的模型真实垂向反射率可以看出,图5d的最小二乘波形偏移结果的振幅保真性相较于图5c的波形偏移结果也存在明显提升。

图5 部分Sigsbee2A模型及其波形偏移结果、最小二乘波形偏移结果和模型的理论垂向反射率

5.2 声波地震数据的波形成像

为了验证声波地震数据的波形成像方法,采用考虑密度变化的声波动方程合成地震数据进行声波阻抗相对扰动的声波反射数据波形偏移试验。

图6a为用于合成地震数据的速度模型,它是在常速为2000m/s的背景模型上叠加一些理论速度值在(-100m/s,100m/s)变化的随机层。图6b为用于合成地震数据的密度模型,它是在常密度为2000kg/m3的背景模型上叠加一些理论密度值在(-50kg/m3,50kg/m3)变化的随机层。图6c为利用图6a和图6b的速度、密度扰动、背景速度和密度得到的声波阻抗相对扰动模型。利用声波方程的有限差分方法合成地震数据,本试验合成了50炮地震数据,炮点范围为0～2km,炮间距为40m,为中间放炮两边接收,每炮由201道接收,道间距为10m,时间采样间隔为1ms,时间采样点数为3000,震源为主频20Hz的Ricker子波。由于图6的模型为层状结构,因此合成的地震数据为反射数据。用于波形成像的速度和密度模型分别为速度为2000m/s的常速模型和密度为2000kg/m3的常密度模型。

利用(26)式至(28)式所表示的波阻抗相对扰动波形偏移计算公式对合成地震反射数据进行波阻抗相对扰动的波形偏移。图7为波形偏移得到的结果,其中图7a为波形偏移的角度域叠加结果,图7b为x=1000m处的角度域共成像点道集,图7c为x=1000m处角度域共成像点道集零度角对应的波阻抗相对扰动波形偏移与波阻抗相对扰动模型值对比结果(右上)和角度域共成像点道集中60°角对应的速度相对扰动波形偏移与速度相对扰动模型值对比结果(右下)。对比图6和图7中的波形偏移结果可以看出,本文提出的声波阻抗相对扰动波形偏移方法可以很好地实现对声波阻抗相对扰动的成像。由图6c 和图7a验证了在(28)式中不进行波场角度分解而直接得到角度域的平均波阻抗相对扰动是一种快速的波阻抗相对扰动成像方法,可用来对地层波阻抗变化进行成像。由图7c中的对比曲线可以看出,波阻抗相对扰动的角度域共成像点道集中的小角度道和大角度道可分别作为波阻抗相对扰动和速度相对扰动的近似估计。图7b的波阻抗相对扰动角度域共成像点道集为进一步的角度域叠前岩性分析和反演提供了基础数据。

图6 用于波形成像试验的速度模型(a)、密度模型(b)和相对声波阻抗相对扰动模型(c)

图7 波阻抗相对扰动的波形偏移结果

上述合成数据的数值试验结果验证了本文提出的标量波和声波的散射数据波形偏移方法、反射数据波形偏移方法和最小二乘波形偏移方法的正确性和有效性。根据波形成像方法的计算公式,我们认为本文提出的波形成像方法具有很强的实用性。弹性波数据的波形成像数值试验结果将在后续的文章中展示。

6 结论

对于波形成像的光滑背景模型,在地震波不发生散射、反射和波型转换的假定下,基于地震数据的线性正演表达,本文首先给出了地震数据波形成像的地震波通用传播方程、散射数据波形成像的通用散射方程和反射数据波形成像的通用反射方程,然后利用线性反演理论建立了面向散射数据和反射数据的地震数据波形成像方法理论框架。散射数据波形成像是对散射体物性参数扰动的线性反演,适合相对于地震波长较小尺寸的散射体和绕射体的成像。反射数据波形成像是对反射体波阻抗相对扰动或反射体边界局部反射系数的线性反演,适合相对于地震波长较大尺寸的反射体物性和地层反射界面的成像。反射数据波形成像不仅能实现对反射界面位置的准确成像,还能得到对地层波阻抗相对变化的岩性成像。在波形成像中,如果利用波传播算子的伴随算子进行地震波的反传播,则波形成像退化为波形偏移,如果利用波传播算子的最小二乘逆算子进行地震波的反传播,则波形成像退化为最小二乘波形偏移。波形偏移计算量小稳定性好,但其偏移成像结果会因数据观测系统和地下复杂构造的影响在保真性和分辨率方面存在不足。最小二乘波形偏移计算量大,但能有效校正数据观测系统和地下复杂构造的影响,因此其偏移成像结果保真性好、分辨率高。我们给出了标量波散射数据、声波反射数据和弹性波反射数据波形成像的计算式,并进行了标量波地震数据和声波地震数据的波形成像数值试验,取得了与理论相一致的试验结果。我们提出的波形成像方法弥补了当前地震数据偏移成像方法理论上存在的不足,可以准确地利用地震数据的波形信息,修正当前逆时偏移方法存在的相位误差,得到的偏移结果没有相位误差,分辨率也有所提高。与常规的逆时偏移和最小二乘逆时偏移的计算相比,波形偏移和最小二乘波形偏移仅增加了对入射波场的时间导数计算。反射体波阻抗相对扰动随反射开角、岩性变化之间的定量表达式和反射数据波形成像得到的波阻抗相对扰动角度域共成像点道集为角度域成像数据的岩性分析和反演提供了方法理论和数据基础,这是我们下一步研究工作的方向。