2000年以来, 随着国内外老油田的勘探程度逐步加深, 油气地震勘探的目标发生了显著的变化, 主要是由以横向缓变的层状介质为主的构造油气藏逐渐转向复杂、小尺度的构造油藏以及隐蔽油气藏。该类油气藏的勘探原则上需要一套新的思路, 使得地震波偏移反演成像的结果满足岩性储层识别与解释的需求^[1-3]。

地震勘探所面对的地下介质可以抽象为横向缓慢变化的广大沉积层中分布着由于火山活动、后期地质构造运动和其它地球动力运动所产生的小尺度速度异常体。断层、裂缝、地层尖灭、粗糙界面、孔洞等(大)小尺度的地质异常体是常见的油气运移通道和(或)储集体, 对它们的刻画和描述是油气勘探的重要目标^[4-5]。在常规的地震处理方法中, 大部分算法都侧重于利用来自层状界面的反射数据, 而来自反射界面以外的同相轴(比如绕射波)通常都被忽视。然而, 绕射波数据中也包含着地下介质的丰富信息, 对绕射波成像非常必要。早在20世纪50年代, 诸多学者就认识到了绕射波在地震勘探中的重要性^[6-7], KREY^[6]指出绕射波是地下非均质地质体最直接的响应, 使用绕射波信息可以提高地震资料处理的分辨率。但是由于绕射波能量往往比反射波能量弱很多, 从而导致了即使将其偏移归位, 通常也会被反射体所掩盖^[8]。因此必须考虑如何单独对绕射波数据进行偏移成像, 突出地下介质中的小尺度地质异常体。

根据绕射波和反射波在运动学和动力学特征上的差异, 可以分别在数据域和成像域中利用两者的特征差异分离反射波与绕射波, 从而实现对绕射波场单独成像的目的^[9-10]。因此, 绕射波成像方法主要分为:数据域绕射波成像方法与成像域绕射波成像方法。数据域中的分离方法有:基于绕射波旅行时方程构建绕射波时间剖面的D-Section分离方法^[11-12]; 聚焦-反聚焦(Focusing-Defocusing)方法^[8], 其利用反射波具有固定镜像虚震源而绕射波在伪深度剖面上发散的特点, 将聚焦的反射波切除再反聚焦从而在全波场中去除反射波; 平面波域中的绕射波分离方法^[13], 其基本思想是在平面波入射的情况下, 地下层状反射界面的地震记录为线形, 而地下绕射点的地震记录为双曲形, 利用滤波的方法滤去线性的反射波从而得到主要含绕射波成分的剖面; 另一些数据域分离方法如CRS技术分离方法^[14]及Multi-focusing方法^[15], 其实质为建立零偏移距剖面的方法。数据域中的分离方法大多基于射线理论, 利用反射波与绕射波在传播旅行时上的不同特征进行分离, 然而射线理论本身并不足以表达复杂地质条件下出现的波现象, 因而这些方法在复杂介质情况下的适用性降低。成像域中的绕射波分离方法主要是在角度域成像中进行。根据反射能量集中在第一菲涅尔带而绕射能量较发散的特点, MOSER等^[16]和KOZLOV^[17]提出了修改Kirchhoff偏移积分公式中的加权函数从而压制反射波能量的反稳相绕射波成像方法, 但由于其使用的是积分类偏移方法实现角度域成像, 绕射波的收敛和聚焦效果都不如波动方程类的偏移方法。基于反射波和绕射波在倾角成像道集上的特征差异, 许多学者提出在倾角成像道集中利用平面波滤波等方法分离绕射波并进行速度分析^[18-21]。

由地下反射界面产生的反射波可以看成不同绕射点的绕射叠加过程, 反射波与绕射波的区别在于反射波遵循Snell定律, 具有方向性, 而绕射波无方向性。因此可以根据地震波传播过程中的传播方向区分绕射波和反射波。在波场传播过程中分解不同方向的波场, 并采用合适的成像条件提取绕射波进行成像是可行的, 但是需要在每个波场外推过程中, 进行波场分解提取对应的波场方向, 计算量巨大, 无法高效和精确地实现。为了避免波场外推中波场分解带来的计算量, 可以分别将来自某一特定方向的地震波进行成像(例如可以分别对左行波和右行波进行成像)。根据绕射波和反射波的方向差异性, 针对地下某一成像点, 反射波只在满足Snell定律的方向上才能成像, 而在其它方向上无法成像, 但是绕射波在所有方向上均可成像。基于此, 本文提出一种基于行波分解的绕射波逆时偏移成像方法, 该方法首先利用解析时间波场外推及波场分解的方法, 在逆时偏移的每一个时间片外推过程中得到震源端及检波点端的下左行波和下右行波, 然后通过修改RTM零延迟相关成像条件, 使用下左和下左行波相关、下右和下右行波相关, 得到“正倾角反射层+绕射体”以及“负倾角反射层+绕射体”两种成像结果, 最后将两种成像结果进行相关并得到最终的绕射波成像结果。该方法旨在较好地压制反射波能量, 有效增强绕射波能量, 实现绕射波成像的目的。最后利用数值实验和实际资料验证该方法的有效性。

1 行波分解理论

波场分解可通过傅里叶变换在频率-波数域实现^[22-24]。利用傅里叶变换可以将地震波场分解为上下、左右行波, 以震源端波场为例, 关于时间和空间的傅里叶变换可表示为:

$ \begin{array}{c} U_{\mathrm{s}}\left(k_{x}, k_{z}, \omega\right)= \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} u_{\mathrm{s}}(x, z, t) \mathrm{e}^{-\mathrm{i}\left(\omega t+x k_{x}+z k_{z}\right)} \cdot \\ \mathrm{d} t \mathrm{d} x \mathrm{d} z \end{array} $

(1)

式中:u_s表示震源端波场; U_s是u_s的傅里叶变换形式。将震源端波场u_s和检波点端波场u_R进行上、下行波分解^[25-26]:

$ \begin{aligned} &u_{\rm s}=u_{\rm s}^{\mathrm{up}}+u_{\mathrm{s}}^{\mathrm{down}}\\ &u_{\mathrm{R}}=u_{\mathrm{R}}^{\mathrm{up}}+u_{\mathrm{R}}^{\text {down }} \end{aligned} $

(2)

式中:上标up代表上行波, down代表下行波。在频率-波数域, 震源端和检波点端的上、下行波可定义为^[26]:

$ \begin{aligned} &U_{\rm s}^{\mathrm{up}}\left(k_{x}, k_{z}, \omega\right)=\left\{\begin{array}{c} U_{\rm s}\left(k_{x}, k_{z}, \omega\right), \omega k_{z}>0 \\ 0, 其它 \end{array}\right.\\ &U_{\rm s}^{\mathrm{down}}\left(k_{x}, k_{z}, \omega\right)=\left\{\begin{array}{c} U_{\rm s}\left(k_{x}, k_{z}, \omega\right), \omega k_{z}<0 \\ 0, 其它 \end{array}\right.\\ &U_{\rm R}^{\text {up }}\left(k_{x}, k_{z}, \omega\right)=\left\{\begin{array}{c} U_{\mathrm{R}}\left(k_{x}, k_{z}, \omega\right), \omega k_{z}<0 \\ 0, 其它 \end{array}\right.\\ &U_{\mathrm{R}}^{\text {down }}\left(k_{x}, k_{z}, \omega\right)=\left\{\begin{array}{c} U_{\mathrm{R}}\left(k_{x}, k_{z}, \omega\right), \omega k_{z}>0 \\ 0, 其它 \end{array}\right. \end{aligned} $

(3)

式中:U_s^up, U_s^down, U_R^up, U_R^down分别是u_s^up, u_s^down, u_R^up, u_R^down的傅里叶变换形式。本文震源端和检波点端的下行波方向都定义为+z方向, 即震源端正方向为波场正传的方向, 检波点端正方向为波场反传的方向, 左、右行波的方向同理。进一步地, 可在上、下行波中分解出左、右行波:

$ \begin{aligned} &u_{\mathrm{s}}^{\mathrm{up}}=u_{\mathrm{s}}^{\mathrm{right}-\mathrm{up}}+u_{\mathrm{s}}^{\mathrm{left}-\mathrm{up}}\\ &u_{\mathrm{s}}^{\text {down }}=u_{\mathrm{s}}^{\text {right-down }}+u_{\mathrm{s}}^{\text {left-down }}\\ &u_{\mathrm{R}}^{\mathrm{up}}=u_{\mathrm{R}}^{\mathrm{right}-\mathrm{up}}+u_{\mathrm{R}}^{\mathrm{left}-\mathrm{up}}\\ &u_{\mathrm{R}}^{\text {down }}=u_{\mathrm{R}}^{\text {right-down }}+u_{\mathrm{R}}^{\text {left-down }} \end{aligned} $

(4)

式中:上标right代表右行波, left代表左行波, right-down表示下右行波, u_s^right-down表示震源端波场下右行波, u_R^right-down表示检波点端波场下右行波。对应地, 下右行波和下左行波也可通过傅里叶变换得到, 以震源端为例:

$ \begin{aligned} U_{\rm s}^{\text {right-down }}\left(k_{x}, k_{z}, \omega\right) &=\left\{\begin{aligned} U_{\rm s}\left(k_{x}, k_{z}, \omega\right), \omega, k_{x}, k_{z}>0 & \;\;\;\; \text { 或} \;\;\;\; \omega, k_{x}, k_{z}<0 \\ 0, & \text { 其它 } \end{aligned}\right.\\ U_{\rm s}^{\text {left-down }}\left(k_{x}, k_{z}, \omega\right) &=\left\{\begin{aligned} U_{\rm s}\left(k_{x}, k_{z}, \omega\right), \omega, k_{z}>0, k_{x} &<0 \;\; \;\;或 \;\; \;\; \omega, k_{z} < 0, k_{x} > 0 \\ 0, & \text { 其它 } \end{aligned}\right. \end{aligned} $

(5)

式中:U_s^right-down, U_s^left-down分别是u_s^right-down, u_s^left-down的傅里叶变换形式。

为了实现上述分解, 需要对全波场进行高维傅里叶变换, 因而需要巨大的计算量及存储量。另外一种有效的方法是采用希尔伯特变换构建解析波场^[27-28], 其实部为波场本身, 虚部为波场在时间方向上的希尔伯特变换结果, 在解析波场中仅包含正频率信息, 因此可以通过空间波数的正、负来判断波场的方向。以震源端为例, 震源端的解析波场为:

$ \tilde{u}_{\mathrm{s}}(x, z, t)=u_{\mathrm{s}}(x, z, t)+\mathrm{i} u_{\mathrm{s}}^{\mathrm{H}}(x, z, t) $

(6)

$ u_{\mathrm{s}}^{\mathrm{H}}(x, z, t)=\mathrm{Hilbert}\left[u_{\mathrm{s}}(x, z, t)\right] $

(7)

式中: $\tilde u_{\rm s}$(x, z, t)表示震源端解析波场; u_s(x, z, t)表示震源端波场; u_s^H(x, z, t)表示震源端波场的希尔伯特变换。想要获得解析波场, 就需要获得实部波场和虚部波场。解析波场的实部可以通过波动方程传播过程得到, 其虚部可以通过对地下波场进行Hilbert变换得到, 但是由于Hilbert变换是关于时间方向的卷积, 这需要将整个波场存盘并进行卷积, 不利于计算。Hilbert变换为一个90°相移算子, 借助于波动方程的震源项和波场的线性关系^[29], 可以把对波场的Hilbert变换转化为震源项的Hilbert变换, 利用新的“震源”进行传播, 从而可以获得解析波场的虚部。以震源端为例, 解析波场外推表示为^[27]:

$ \left\{\begin{array}{l} \frac{1}{v^{2}(\boldsymbol{x})} \frac{\partial^{2} u_{\mathrm{s}}\left(\boldsymbol{x}, t ; \boldsymbol{s}_{ \boldsymbol{x}}\right)}{\partial t^{2}}- \boldsymbol{\nabla}^{2} u_{\mathrm{s}}\left(\boldsymbol{x}, t ; \boldsymbol{s}_{\boldsymbol{x}}\right)=f_{\mathrm{s}}\left(t, \boldsymbol{s}_{\boldsymbol{x}}\right) \\ \frac{1}{v^{2}(\boldsymbol{x})} \frac{\partial^{2} u_{\mathrm{s}}^{\mathrm{H}}\left(\boldsymbol{x}, t ; \boldsymbol{s}_{\boldsymbol{x}}\right)}{\partial t^{2}}- \boldsymbol{\nabla}^{2} u_{\mathrm{s}}^{\mathrm{H}}\left(\boldsymbol{x}, t ; \boldsymbol{s}_{\boldsymbol{x}}\right)=f_{\mathrm{s}}^{\mathrm{H}}\left(t, \boldsymbol{s}_{\boldsymbol{x}}\right) \\ \tilde{u}_{\mathrm{s}}=u_{\mathrm{s}}+\mathrm{i} u_{\mathrm{s}}^{\mathrm{H}} \end{array}\right. $

(8)

式中:f_s为地震子波震源项; v为地下介质速度; u_s为震源波场, 即解析时间波场的实部; f_s^H为f_s的Hilbert变换; u_s^H为u_s的Hilbert变换, 即解析波场的虚部; $\tilde u_{\rm s}$就是震源端解析波场; x是空间点坐标的矢量形式; s_x是震源点坐标的矢量形式。

同理, 检波点端波场的反传公式为^[27]:

$ \left\{\begin{array}{l} \frac{1}{v^{2}( \boldsymbol{x})} \frac{\partial^{2} u_{\mathrm{R}}\left(\boldsymbol{x}, t ; \boldsymbol{s}_{\boldsymbol{x}}\right)}{\partial t^{2}}- \boldsymbol{\nabla}^{2} u_{\mathrm{R}}\left( \boldsymbol{x}, t ; \boldsymbol{s}_{\boldsymbol{x}}\right)= \\ d_{\mathrm{R}}\left( \boldsymbol{g}_{ \boldsymbol{x}}, t ; \boldsymbol{s}_{ \boldsymbol{x}}\right) \\ \frac{1}{v^{2}( \boldsymbol{x})} \frac{\partial^{2} u_{\mathrm{R}}^{\mathrm{H}}\left( \boldsymbol{x}, t ; \boldsymbol{s}_{ \boldsymbol{x}}\right)}{\partial t^{2}}- \boldsymbol{\nabla}^{2} u_{\mathrm{R}}^{\mathrm{H}}\left( \boldsymbol{x}, t ; \boldsymbol{s}_{ \boldsymbol{x}}\right)= \\ d_{\mathrm{R}}^{\mathrm{H}}\left( \boldsymbol{g}_{ \boldsymbol{x}}, t ; \boldsymbol{s}_{ \boldsymbol{x}}\right) \\ \tilde{u}_{\mathrm{R}}=u_{\mathrm{R}}+\mathrm{i} u_{\mathrm{R}}^{\mathrm{H}} \end{array}\right. $

(9)

式中:d_R为地表记录的地震数据; d_R^H为d_R的Hilbert变换; u_R^H为u_R的Hilbert变换; $\tilde u_{\rm R}$为检波点端解析波场; g_x为检波点坐标的矢量形式。

得到解析时间波场后, 就可以在解析波场外推的每个时间层上进行空间傅里叶变换, 将时间-空间域的波场转换到频率-波数域, 从而利用(5)式来实现波场的下左和下右行波分量的分解了。图 1为源检端下左、下右行波的分解结果, 由图 1可以看出, 该方法可以较好地对波场进行下左和下右的方向分解。

2 基于行波分解的绕射波逆时偏移成像方法

ZHANG等^[30-31]提出在时间-空间域通过构建波场的希尔伯特变换对正、负倾角反射层分别进行成像, 从而实现分离反射波、提取绕射波并对其单独成像的目的。该方法的核心思想是利用绕射波与反射波的方向性差异, 在两次成像中, 反射层只存在单一的成像剖面中, 而绕射体在两次成像的剖面中都存在。由于该方法无法区分垂直向下传播的波场, 在下左行波和下右行波分量中均含有垂直向下传播的波场, 因此水平反射层会作为干扰留在绕射波成像结果中, 并且该方法需要存储虚部波场, 对硬盘存储和I/O均有较高的要求。本文提出通过构建解析时间波场, 在时间-波数域通过傅里叶变换进行波场方向分解, 并对分解的方向波场分别成像, 然后将成像结果进行相关, 从而实现绕射体成像的目的。该方法不需要存储虚部波场, 并且可以灵活地选择波场方向范围, 通过对波场传播角度的控制, 消除垂直向下传播的波场, 有效压制了水平反射层的干扰。

2.1 传统逆时偏移成像方法

逆时偏移(RTM)技术在20世纪80年代由MCMECHAN^[22]、BAYSAL等^[32-34]提出并取得了一系列的发展。基于波动理论的逆时偏移, 由于其不受倾角和偏移孔径的限制, 可以有效地处理纵横向剧烈变化的地球介质物性特征(如速度、密度等), 是现行偏移方法中最精确的一种成像方法。传统的RTM零延时互相关成像条件公式可以表示为:

$ I=\sum\limits_{n_{\text {shot }}} \sum\limits_{n_{t}} u_{\rm s} \cdot u_{\mathrm{R}} $

(10)

式中:n_t代表时间采样点数; n_shot为炮数; I为成像结果。直接采用(10)式作为成像条件会引起成像假象和低频噪声, 为避免成像假象和低频噪声, LIU等^[26]提出将波场进行上、下行波分解, 利用震源端和检波点端同向传播的波场进行相关成像, 成像条件变为:

$ \begin{aligned} I &=\sum\limits_{n_{\mathrm{shot}}} \sum\limits_{n_{t}}\left(u_{\mathrm{s}}^{\mathrm{down}} \cdot u_{\mathrm{R}}^{\mathrm{down}}+u_{\mathrm{s}}^{\mathrm{up}} \cdot u_{\mathrm{R}}^{\mathrm{up}}\right) \\ &=I_{1}+I_{2} \end{aligned} $

(11)

式中:I₁为震源端下行波场与检波点端下行波场相关成像结果; I₂为震源端上行波场与检波点端上行波场相关成像结果。在宏观速度纵向梯度较大, 且存在高陡构造时, I₂成像分量才有显著贡献, 在不存在高陡构造的情况下, 可以忽略震源端上行波和检波点端上行波相关成像的结果。因此, 若只采用震源端下行波和检波点端下行波相关进行成像, 则修改后的成像条件变为:

$ I=I_{1}=\sum\limits_{n_\text { shot }} \sum\limits_{n_{t}} u_{\rm s}^{\text {down }} \cdot u_{\mathrm{R}}^{\text {down }} $

(12)

2.2 绕射波逆时偏移成像方法

绕射波是由于地下地质异常体的尺度小于地震波波长而产生的, 绕射波与反射波在运动学特征上存在一定的差异, 绕射波的能量通常比反射波低一到两个量级。如图 2所示, 反射波的传播符合Snell定律, 即入射角等于反射角, 地震波入射到反射体后, 反射波只沿着符合Snell定律的方向发生反射, 因此反射波具有方向特征。根据入射波场与反射波场的方向, 可以计算反射层的倾角, 相应地, 特定的反射层的入射射线与反射射线的方向与该反射层的倾角有关。而绕射波的传播不符合Snell定律, 地震波入射到绕射体后, 绕射波会向各个方向出射, 并由地面的检波器接收, 绕射波的传播不具有方向特征。这一特征差异可被用于反射波和绕射波的分离成像。

将(12)式中的下行波分解成下左行波和下右行波, 则(12)式可以写成:

$ \begin{aligned} I_{1}=& \sum\limits_{n \text { shot }} \sum\limits_{n_{t}} u_{\rm s}^{\text {down }} \cdot u_{\mathrm{R}}^{\text {down }} \\ =& \sum\limits_{n_{\text {shot }}} \sum\limits_{n_{t}}\left(u_{\rm s}^{\text {rd }} \cdot u_{\mathrm{R}}^{\text {rd }}+u_{\rm s}^{\text {ld }} \cdot u_{\mathrm{R}}^{\text {ld }}+u_{\rm s}^{\mathrm{rd}} \cdot u_{\mathrm{R}}^{\text {ld }}+\right.\\ &\left.u_{\rm s}^{\text {ld }} \cdot u_{\mathrm{R}}^{\mathrm{rd}}\right)=I_{\mathrm{rr}}+I_{\rm ll}+I_{\mathrm{rl}}+I_{\mathrm{lr}} \end{aligned} $

(13)

式中:

$ I_{\rm rr}=\sum\limits_{n_{\rm shot}}\sum\limits_{n_{\rm t}}u^{\rm right-down}_{\rm s}·u^{\rm right-down}_{\rm R} $

(14)

$ I_{\rm ll}=\sum\limits_{n_{\rm shot}}\sum\limits_{n_{\rm t}}u^{\rm left-down}_{\rm s}·u^{\rm left-down}_{\rm R} $

(15)

$ I_{\rm rl}=\sum\limits_{n_{\rm shot}}\sum\limits_{n_{\rm t}}u^{\rm right-down}_{\rm s}·u^{\rm left-down}_{\rm R} $

(16)

$ I_{\rm lr}=\sum\limits_{n_{\rm shot}}\sum\limits_{n_{\rm t}}u^{\rm left-down}_{\rm s}·u^{\rm right-down}_{\rm R} $

(17)

其中, I_rr为震源端下右行波和检波点端下右行波相关成像结果, 如图 3a所示; I_ll为震源端下左行波和检波点端下左行波相关成像结果, 如图 3b所示; I_rl为震源端下右行波和检波点端下左行波相关成像结果, I_lr为震源端下左行波和检波点端下右行波相关成像结果, 如图 3c所示。反射层可以根据法向向量方向分为正倾角与负倾角两种反射层, 如图 4所示, 本文定义正倾角反射层为朝向地面的法向向量与向右的水平线夹角为锐角(0~90°)的反射层, 负倾角反射层为朝向地面的法向向量与向右的水平线夹角为钝角(90°~180°)的反射层。

图 5为正倾角反射层成像原理示意图。由图 5可见, 在震源端入射波场为下右行波的情况下, 由几何关系可知, 检波点端出射波场必定为下左行波(图 5中红色实线所示), 也即不可能为下右行波, 因此I_rr的成像结果中没有正倾角反射层的像; 在震源端波场为下左行波时, 检波点端出射波场可能为下左行波(图 5中绿色实线所示), 因此I_ll的成像结果中含有正倾角反射层的像。同理, I_ll的成像结果中没有负倾角反射层的像, 而I_rr的成像结果中含有负倾角反射层的像。而绕射波由于不符合Snell定律, 波场入射到绕射体上, 绕射体向各个方向出射绕射波, 因此绕射体在I_rr、I_ll中均成像。此时可以认为I_ll代表了正倾角反射层和绕射体的像, I_rr代表了负倾角反射层和绕射体的像, 由于正、负倾角的反射层都只在I_rr或I_ll其中一项中成像, 在另一项中不成像, 而绕射体在两项中均可以成像, 对两项结果进行相关, 即可以将正、负倾角的反射层的像都去除, 突出绕射体的像。具体的成像公式如下:

$ I_{\rm diffraction}=I_{\rm rr}·I_{\rm ll} $

(18)

式中:I_diffraction为绕射波成像结果。

在水平反射层或地层倾角较小的情况下, 波场分解的过程中, 由于垂直向下传播的波场既在下左行波分量中存在, 也在下右行波分量中存在, 因此水平反射层(包括小倾角反射层)既在正倾角反射层的像中存在, 也在负倾角反射层的像中存在, 利用(18)式无法压制水平反射层的能量。为了解决这一问题, 在频率-波数域提取不同象限的波场进行波场分解时, 可以在每个象限的边界处设置衰减窗函数, 这样就可以灵活地分解得到特定角度范围的下左和下右行波分量。在下左和下右行波分量中都去除垂直向下传播的波场(以及沿小角度传播的波场), 从而实现压制水平反射层(小倾角反射层)能量的目的, 而去除一定角度范围的下左行波和下右行波分量对绕射波成像的影响很小。在时间-空间域通过卷积进行希尔伯特变换从而实现波场分解的方法无法去除沿特定角度范围传播的波场, 因此选择利用傅里叶变换在频率-波数域实现波场分解对压制反射层的能量更加有效。

3 数值试验 3.1 模拟数据

首先采用“凹”形模型验证本文方法的有效性。图 6a为速度模型, 该速度模型存在一个正倾角的反射层、一个负倾角的反射层以及4个绕射点; 图 6b为常规的RTM成像结果; 图 6c和图 6d分别为I_ll与I_rr的成像结果, 分别代表了“正倾角反射层+绕射体”的像和“负倾角反射层+绕射体”的像; 图 6e为利用本文方法得到的绕射体成像结果。对比图 6c与图 6d可见, 在I_ll的成像结果中只存在正倾角的反射层与绕射体的像, 其中两个绕射体清晰可见, 但依然有两个绕射体被正倾角反射层掩盖, 难以有效识别; 在I_rr的成像结果中只存在负倾角的反射层与绕射体的像, 与I_ll的结果正好相反。因此利用公式(18)进行相关后可以有效压制反射层的像, 增强绕射体的成像结果。

进一步地, 为了验证复杂模型下本文方法的有效性, 采用Sigsbee模型进行测试。图 7a为Sigsbee偏移速度模型, 速度模型以层状反射地层为主, 并发育数组断层, 此外, 该模型还存在大量由粗糙界面引起的绕射点, 第一炮在x=0处激发, 共50炮; 图 7b为常规RTM成像结果; 图 7c和图 7d分别为I_ll与I_rr成像结果; 图 7e为采用本文方法得到的成像结果。由图 7可以看出, 常规RTM成像中, 绕射体被能量更强的反射界面掩盖, 难以识别, 而采用本文方法得到的结果, 可以有效地对绕射体进行识别, 并更加清晰地指示了断层的位置。为验证本文方法对速度误差的敏感性, 采用不同的偏移速度模型进行成像, 实验结果如图 8所示。其中图 8a为采用95%真实Sigsbee速度模型平滑之后的偏移速度模型; 图 8b为采用80%真实Sigsbee速度模型平滑之后的偏移速度模型。受偏移速度误差的影响, 常规成像结果中成像位置发生变化, 绕射体的成像不再聚焦, 但采用本文方法仍然可以有效地压制反射波, 较好地突出绕射波成像, 识别绕射体的相对位置。

3.2 实际地震数据

采用某工区二维实际地震数据测试本文方法的有效性, 该数据包含130炮, 每炮间隔150 m, 时间采样为4 s, 采样间隔为0.5 ms, 双边观测, 最大偏移距为1 976 m。图 9a为实际地震资料常规逆时偏移成像结果; 图 9b为采用本文方法得到的绕射波成像结果。在浅层地区, 反射波同相轴表现为水平层状, 上、下反射界面分界明显, 同相轴的分辨率较高。在横向位置1.0 km和1.8 km处出现了明显的断层, 该断层延伸至深度1 400 m左右的不整合面上。从图 9可以看出, 在常规成像剖面中, 断层和绕射体被掩盖掉, 无法识别; 而在本文方法的成像结果中, 几乎没有连续的反射层的像存在, 断层和绕射点可以清晰地识别。由于本文方法可以灵活地选择波场方向, 不存在垂直向下传播的波场无法区分的问题, 因此浅层地区的水平层状反射层几乎完全被去掉。实际资料应用结果说明了本文提出的绕射波成像方法的有效性。

4 结论与讨论

绕射波成像对油气勘探具有独特的价值, 本文发展了一种基于波场方向分解理论的绕射波逆时偏移成像方法。该方法基于反射波和绕射波传播方向的特征差异, 采用波场方向分解理论构建了针对绕射波成像的逆时偏移成像条件, 最终实现了绕射波成像的目的。为了在逆时偏移过程中高效便利地进行波场方向分解, 采用了基于解析时间波场方向分解的波场传播方法。基于波场方向分解成像条件获得“正倾角反射层+绕射体”及“负倾角反射层+绕射体”两种成像结果, 基于绕射波传播角度特征, 将二者进行相关处理获得绕射波成像结果。凹陷模型和复杂模型数据的测试表明本文方法能够很好地获取绕射波成像结果, 模型中断点的成像清晰。同时, 在速度存在一定误差的情况下, 仍然可对绕射点的相对位置进行成像。本文方法在实际数据中的实验结果表明该方法能够很好地对断点及潜山顶界面不连续位置进行成像, 体现出一定的实际应用潜力。

数值试验表明:本文方法在速度存在一定误差情况下仍能获取绕射点的相对位置成像结果, 但其成像位置不准确。在实际数据处理中, 为了精确定位绕射体, 仍然需要高精度的偏移速度模型。本文方法仅使用下行波进行成像, 与单程波偏移类似, 对于缝洞较为发育的探区存在的多次绕射及多次波无法很好地利用, 将该方法拓展到双程波的情况下是进一步研究的方向。此外, 在波场方向分解中为了消除水平层的影响, 去除了少量绕射波信息, 也会对绕射波成像振幅存在一定的影响。对于复杂波场, 如何准确估计其传播方向仍然值得进一步研究。

致谢: 感谢中石油勘探开发研究院及西北分院、中海油研究院和湛江分公司、中国石化石油物探技术研究院和胜利油田分公司对波现象与智能反演成像研究组(WPI)研究工作的资助与支持。