CN102890290A - 一种起伏地表条件下的叠前深度偏移方法 - Google Patents

Abstract

一种起伏地表条件下的叠前深度偏移方法，包括：步骤1：以最高高程为基准面z_beg＝0，在基准面与起伏地表之间填充近地表速度；步骤2：假设在最初基准面z_beg=0时初始波场，将z_beg=z=0时对应的波场加入到初始波场；步骤3：在单炮偏移孔径范围内将波场划分成多个小波束窗并变换到小波束域，分别将下行波场和上行波场沿深度方向延拓至深度z+Δz；步骤4：对延拓后的波场进行滤波处理；步骤5：根据互相关成像条件，应用下行波场和上行波场进行互相关成像，得到深度层z+Δz的成像结果；步骤6：判断深度z+Δz层是否存在波场，若是，将相应波场叠加到所述延拓波场中，生成深度z+Δz层的开始延拓波场；步骤7：重复步骤3至步骤6，逐步延拓并成像直至最大深度层z_max，得到单炮的成像结果。

Description

一种起伏地表条件下的叠前深度偏移方法

技术领域

本发明是关于石油勘探中地震数据处理领域的叠前深度偏移成像技术，具体的讲是关于一种起伏地表条件下的叠前深度偏移方法。

背景技术

随着地震勘探技术的不断发展，我国油气勘探的战略重点正逐渐转向复杂地质条件的区域。但是对于复杂地质条件的区域，比如中国西部山区这样地表条件复杂的地区，基于水平叠加的常规地震数据处理遇到了很大的困难。常规的高程静校正可以解决起伏地表的影响，但是这种方法的实现是基于地表一致性假设，只有在地表起伏不大、低速带横向变化缓慢，地下浅、中、深层的反射经过低速带时几乎垂直入射至地表的情况下才适用，当遇到高程较大，横向变化剧烈地区的时候，高程静校正的效果就大打折扣了。

针对高程基准面校正所带来的问题，许多学者做了大量的工作，Berryhill在1979年首先提出波动方程基准面校正的概念。当时他用Kirchhoff积分法在零炮检距理论记录上，把一个已知的任意形状的观测面上的波场延拓到另一个特定的基准面上，从而完成波动方程基准面校正。几年后，他又将这个思路扩展到叠前（Berryhill,1984）。基于这个思路，Berryhill（1984,1986）、Wibbins（1984）、Yilmaz和Lucas（1986）、Malloyetal（1990）、Schneider et al（1995）先后就波动方程基准面校正并结合层替换做了许多有益的尝试，使很多实际问题在某种程度上得到解决。

基于射线理论的成像方法可以直接从起伏的地表面进行波场偏移成像，对复杂地表构造具有较好的适应性。Jager等在2003年提出了起伏地表条件下的真振幅Kirchhoff积分偏移方法。Gray在2005年提出了一种适用于复杂地表条件下的高斯束偏移方法，然而，当地表高程及速度变化较大时，静校正对波场造成的畸变会对后续的偏移成像特别是近地表的成像造成不利影响。此外，当遇到速度场的强横向变化时，基于射线理论的偏移方法成像效果会显著下降。

为了解决地表起伏变化剧烈对地下构造成像的影响，Reshef在1991年提出直接从非水平观测面开始的“逐步外推，逐步累加”的波场外推方式来实现深度偏移。为解决高程静校正时移所带来的误差，将非水平观测面变为水平观测面，以便使用常规的偏移算子进行波场延拓，方便地做偏移，Beasley和Lynn在1992年提出了“零速层”的概念，在起伏面和基准面间插入了一个速度很小的虚拟层，也能够消除起伏地表的影响，不过在起伏较大的情况下计算不稳定，效果较差。Yang等1999年应用的基于波动方程深度外推算子“直接下延”法和何英等2002年应用的“波场上延”法都能够较好的消除起伏地表的影响，然而由于使用全局的参考速度，使得偏移精度不够，浅层起伏地表附近的构造成像被部分噪音所掩盖。Shragge在2008年提出了黎曼坐标系下进行波场延拓，将坐标系定义成与起伏面相似的情况，使得传统偏移方法可以被应用，但是由于需要对数据进行插值等规则化处理往往需要巨大的额外计算量。

发明内容

本发明是关于一种起伏地表条件下的叠前深度偏移方法，以提高起伏地表地区的叠前深度偏移精度，改善成像质量。

为了实现上述目的，本发明提供一种起伏地表条件下的叠前深度偏移方法，所述的方法包括如下步骤：

步骤1：以最高高程为基准面z_beg=0，在基准面与起伏地表之间填充近地表速度；步骤2：假设在最初基准面z_beg=0时初始波场u(x,z_beg)=0，记录波场为u_e(x,z)，将z_beg=z=0时对应的波场加入到初始波场u(x,z_beg)=u_e(x,z)；步骤3：在单炮偏移孔径范围内将波场u(x,z)划分成多个小波束窗，将波场u(x,z)变换到小波束域，分别将下行波场和上行波场沿深度方向延拓Δz至深度z+Δz，并变回到空间域，延拓后的下行波场为

和延拓后的上行波场为

步骤4：对波场下行波场和上行波场

进行滤波处理；步骤5：根据互相关成像条件对所述下行波场

和上行波场

进行互相关成像，得到深度层z+Δz的成像结果；步骤6：判断深度z+Δz层是否存在所述波场u_e(x,z+Δz)，如果存在，将相应波场叠加到所述延拓波场中，生成深度z+Δz层的开始延拓波场u(x,z+Δz)=u_con(x,z+Δz)+u_e(x,z+Δz)；步骤7：重复所述步骤3至步骤6，逐步延拓并成像直至最大深度层z_max，得到单炮的成像结果，最终的偏移结果为所有单炮成像结果的叠加。

进一步地，所述的步骤3包括：

将每一深度层z处的波场分解为：

$u(x,z,w)=\underset{m}{\mathrm{\Σ}}\underset{n}{\mathrm{\Σ}}\hat{u}({\stackrel{\‾}{x}}_n,{\stackrel{\‾}{\mathrm{\ξ}}}_m;z)b_{\mathrm{mn}}\left(x\right)---\left(1\right)$

将（1）式代入频率-空间域波动方程 $[{\∂}_x^2+{\∂}_y^2+\frac{w^2}{v^2(x,y)}]u(x,z,w)=0$ 中得到标量声波方程 $\underset{n}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}\hat{u}({\stackrel{\‾}{x}}_n,{\stackrel{\‾}{\mathrm{\ξ}}}_m;z)[{\∂}_x^2+{\∂}_z^2+\frac{w^2}{v^2(x,z)}]b_{\mathrm{mm}}\left(x\right)=0;$

将小波束的演化表示为 $a_{\mathrm{mn}}\left(x\right)=e^{\±i{A}_n\mathrm{\Δz}}{b}_{\mathrm{mn}}\left(x\right);$

根据波场延拓算子A_n将所述波场u(x,z)延拓至下一深度层z+Δz，生成延拓后的波场u_con(x,z+Δz)；

其中，w是频率；v(x,y)代表层速度；u(x,z,w)是地表(x,z)处接收的地震记录；b_mn(x)代表空间位置波数为的小波束；是一系列关于深度z的系数，a_mn是小波束b_mn经过介质传播后的波束。

进一步地，波场延拓算子

其中， $\mathrm{\Δ}{k}_n\left(x\right)=w((1/v(x,z))-(1/v_0({\stackrel{\‾}{x}}_n,z)))$ 是局部扰动。

进一步地，所述的步骤4包括：对延拓波场进行如下滤波处理：u_cor(x,z+Δz)=u_cor(x,z+Δz)filter(x,z+Δz)，滤波函数filter(x,z+Δz)为：

$\mathrm{filter}(x,z+\mathrm{\Δz})=\left\{\begin{array}{c}0,(x,z+\mathrm{\Δz})\∈\mathrm{Fictitious\; layer}\\ 1,(x,z+\mathrm{\Δz})\∈\mathrm{Others}\end{array}\right..$

进一步地，所述的步骤5包括：

根据 $R(x,z+\mathrm{\Δz})=\frac{1}{2\mathrm{\π}}\text{\∫}{u}_{\mathrm{con}}^r(x,z+\mathrm{\Δz})u_{\mathrm{con}}^{s*}(x,z+\mathrm{\Δz})\mathrm{dw}$ 对震源波场

和记录波场

进行互相关成像生成z+Δz层的成像结果，其中，

是震源波场的复共轭。

不发明实施例的有益效果在于：提高了起伏地表地区的叠前深度偏移精度，改善了成像质量，尤其是在浅层起伏面附近和速度变化剧烈地区的改善尤为明显。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。在附图中：

图1为本发明实施例起伏地表条件下的叠前深度偏移方法流程图；

图2a是起伏地表的高程示意图；

图2b是速度-深度模型示意图；

图2c是模拟的单炮记录示意图；

图3a是15Hz的小波束示意图；

图3b是40Hz的小波束示意图；

图4a是基于分布傅立叶算子的起伏地表叠前深度偏移示意图；

图4b是基于局部窗化小波束算子的起伏地表叠前深度偏移示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚明白，下面结合附图对本发明实施例做进一步详细说明。在此，本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。

如图1所示，本发明提供一种起伏地表条件下的叠前深度偏移方法，该方法包括如下步骤：

步骤S101：以最高高程为基准面z_beg＝0，在基准面与起伏地表之间填充近地表速度。

步骤S102：假设在最初基准面z_beg＝0的时候，初始波场u(x,z_beg)=0，实际检波器和震源分布在起伏面上，记录波场为u_e(x,z)，当z_beg＝z=0时，将对应的波场加入到初始波场u(x,z_beg)=u_e(x,z)。

步骤S103：确定小波束窗的长度，在单炮偏移孔径范围内划分成多个小波束窗，将波场u(x,z)变换到小波束域，小波束窗内的平均速度为波场延拓算子中的背景参考速度，分别将下行波场和上行波场沿深度方向延拓Δz至深度z+Δz，并变回到空间域，延拓后的下行波场为

和延拓后的上行波场为

下行波场

又称为震源波场，上行波场又称为记录波场。

二维声波情况下，频率-空间域波动方程可以写成

$[{\∂}_x^2+{\∂}_y^2+\frac{w^2}{v^2(x,y)}]u(x,z,w)=0---\left(2\right)$

其中，w是频率；v(x,y)代表层速度；u(x,z,w)是地表(x,z)处接收的地震记录。对于小波束的偏移方法，在每一个深度层上，将波场分解到小波域，在进行波场延拓，在每一个深度层z处，波场可以分解为

$u(x,z,w)=\underset{m}{\mathrm{\Σ}}\underset{n}{\mathrm{\Σ}}\hat{u}({\stackrel{\‾}{x}}_n,{\stackrel{\‾}{\mathrm{\ξ}}}_m;z)b_{\mathrm{mn}}\left(x\right)---\left(1\right)$

其中，b_mn(x)代表空间位置

波数为

的小波束。

将式(1)代入式(2)后得到如下的标量声波方程

$\underset{n}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}\hat{u}({\stackrel{\‾}{x}}_n,{\stackrel{\‾}{\mathrm{\ξ}}}_m;z)[{\∂}_x^2+{\∂}_z^2+\frac{w^2}{v^2(x,z)}]b_{\mathrm{mn}}\left(x\right)=0---\left(3\right)$

式中

是一系列关于深度z的系数，这些系数是常数。传播的影响也就包含在小波束的演化中。

对于单程波近似，可以将小波束的演化表示为下式：

$a_{\mathrm{mn}}\left(x\right)=e^{\±i{A}_n\mathrm{\Δz}}{b}_{\mathrm{mn}}\left(x\right)---\left(4\right)$

a_mn是小波束b_mn经过介质传播后的波束，A_n就是波场延拓的平方根算子。在局部扰动近似中，将在每个窗

内引入局部参考速度

基于这个参考速度，在每个点处计算其速度扰动。由于局部参考速度对于横向变速的适应性，使得局部点处的速度扰动相对来说也比较小，可以在每个窗内进行相移校正。

其中， $\mathrm{\Δ}{k}_n\left(x\right)=w((1/v(x,z))-(1/v_0({\stackrel{\‾}{x}}_n,z)))$ 是局部扰动。波场延拓算子由I和II项构成，这里我们需要注意的是背景参考速度已经不是传统的全局参考速度，而是局部窗内的参考速度这样大大的减少了计算误差。通过(5)式得波场延拓算子，将波场u(x,z)延拓至下一个深度层z+Δz，延拓后的波场为u_con(x,z+Δz)。

步骤S104：为了消除起伏面与基准面之间的噪音，引入滤波函数，对波场下行波场和上行波场

进行滤波处理。

在波场延拓过程中，基准面与起伏面之间的虚拟层中存在部分偏移噪音，为了消除这种噪音，本发明对波场进行了(6)式所示滤波，

u_cor(x,z+Δz)=u_cor(x,z+Δz)filter(x,z+Δz)        (6)

其中，滤波函数为

$\mathrm{filter}(x,z+\mathrm{\Δz})=\left\{\begin{array}{c}0,(x,z+\mathrm{\Δz})\∈\mathrm{Fictitious\; layer}\\ 1,(x,z+\mathrm{\Δz})\∈\mathrm{Others}\end{array}\right.---\left(7\right).$

步骤S105：应用互相关成像条件对步骤(4)得到的下行波场

和上行波场

进行互相关成像，得到深度层z+Δz的成像结果。

通过小波束波场延拓算子沿深度方向的延拓后，得到了第z+Δz层的下行波场

和上行波场

通过互相关成像条件(8)式，得到该层的成像结果，

$R(x,z+\mathrm{\Δz})=\frac{1}{2\mathrm{\π}}\text{\∫}{u}_{\mathrm{con}}^r(x,z+\mathrm{\Δz};w)u_{\mathrm{con}}^{s*}(x,z+\mathrm{\Δz};w)\mathrm{dw}---\left(8\right)$

其中，

是下行波场的复共轭。

步骤S106：在深度z+Δz层，判断该层是否存在波场u_e(x,z+Δz)，如果存在，将相应波场叠加到延拓波场中，得到z+Δz层的开始延拓波场u(x,z+Δz)=u_con(x,z+Δz)+u_e(x,z+Δz)。

在起伏地表情况下，震源和检波器分布在具有不同高程在起伏面上，基于小波束算子从基准面开始向下延拓波场，在每个点位置处的波场u(x,z+Δz)由两部分构成，一部分是上层延拓下来的波场u_con(x,z+Δz)，，另一部分是该位置处本身存在的波场u_e(x,z+Δz)，第z+Δz层起始延拓波场为

u(x,z+Δz)=u_con(x,z+Δz)+u_e(x,z+Δz)           (9)

步骤S107：循环步骤103至步骤S 106，逐步延拓并成像直至最大深度层z_max，得到单炮的成像结果，最终的偏移结果是所有单炮成像结果的叠加。

下面通过加拿大逆掩断层起伏地表模型来测试本发明的正确性：

图2a是该起伏地表模型的高程，最大高程接近2km。图2b是模型的偏移速度场，该模型的地表起伏较大，并且存在许多断层和褶皱构造。图2c是原始单炮记录，由于起伏地表的影响，同相轴有明显的扭曲。模拟的激发方式是中间放炮，最大道数是480道，正演炮数据共有277炮，每道2000个采样点，4毫秒采样，道间距为15米，最大偏移距为3600米，最小偏移距为15米，炮间距为90米。

在实现起伏地表条件下的偏移之前，首先要实现小波束在起伏地表中的传播，应用本发明所述的方法，分别可以得到不同频率下小波束在起伏地表情况下的传播，图3a是15Hz小波束的传播示意图，图3(b)是40Hz小波束的传播示意图。小波束的起始位置位于起伏面上，而不是在基准面上。图4a是基于传统算子（分布傅立叶）得到的277炮叠前深度偏移结果，图4b是基于局部窗化小波束算子的起伏地表叠前深度偏移方法的叠前深度偏移结果，对比可以看出，明显地在整体效果上有改善，尤其是在浅层起伏地表附近和速度变化剧烈地区，成像效果有了明显的提高。

本发明的起伏地表条件下的叠前深度偏移方法提高了起伏地表地区的叠前深度偏移精度，改善了成像质量，尤其是在浅层起伏面附近和速度变化剧烈地区的改善尤为明显，是一种具有较高应用价值的起伏地表地区叠前深度偏移方法。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种起伏地表条件下的叠前深度偏移方法，其特征在于，所述的方法包括如下步骤：

步骤1：以最高高程为基准面z_beg=0，在基准面与起伏地表之间填充近地表速度；

步骤2：假设在最初基准面z_beg=0时初始波场u(x,z_beg)=0，记录波场为u_e(x,z)，将z_beg=z=0时对应的波场加入到初始波场u(x,z_beg)=u_e(x,z)；

步骤3：在单炮偏移孔径范围内将波场u(x,z)划分成多个小波束窗，将波场u(x,z)变换到小波束域，分别将下行波场和上行波场沿深度方向延拓Δz至深度z+Δz，并变回到空间域，延拓后的下行波场为

和延拓后的上行波场为

步骤4：对波场下行波场

和上行波场

进行滤波处理；

步骤5：根据互相关成像条件对所述下行波场

和上行波场进行互相关成像，得到深度层z+Δz的成像结果；

步骤6：判断深度z+Δz层是否存在波场u_e(x,z+Δz)，如果存在，将相应波场叠加到所述延拓波场中，生成深度z+Δz层的开始延拓波场u(x,z+Δz)=u_con(x,z+Δz)+u_e(x,z+Δz)；

步骤7：重复所述步骤3至步骤6，逐步延拓并成像直至最大深度层z_max，得到单炮的成像结果，最终的偏移结果为所有单炮成像结果的叠加。

2.根据权利要求1所述的起伏地表条件下的叠前深度偏移方法，其特征在于，所述的步骤3包括：

将每一深度层z处的波场分解为：

$u(x,z,w)=\underset{m}{\mathrm{\Σ}}\underset{n}{\mathrm{\Σ}}\hat{u}({\stackrel{\‾}{x}}_n,{\stackrel{\‾}{\mathrm{\ξ}}}_m;z)b_{\mathrm{mn}}\left(x\right)---\left(1\right)$

将（1）式代入频率-空间域波动方程 $[{\∂}_x^2+{\∂}_y^2+\frac{w^2}{v^2(x,y)}]u(x,z,w)=0$ 中得到标量声波方程 $\underset{n}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}\hat{u}({\stackrel{\‾}{x}}_n,{\stackrel{\‾}{\mathrm{\ξ}}}_m;z)[{\∂}_x^2+{\∂}_z^2+\frac{w^2}{v^2(x,z)}]b_{\mathrm{mn}}\left(x\right)=0;$

将小波束的演化表示为 $a_{\mathrm{mn}}\left(x\right)=e^{\±i{A}_n\mathrm{\Δz}}{b}_{\mathrm{mn}}\left(x\right);$

根据波场延拓算子A_n将所述波场u(x,z)延拓至下一深度层z+Δz，生成延拓后的波场u_con(x,z+Δz)；

其中，w是频率；v(x,y)代表层速度；u(x,z,w)是地表(x,z)处接收的地震记录；b_mn(x)代表空间位置

波数为

的小波束；是一系列关于深度z的系数，a_mn是小波束b_mn经过介质传播后的波束。

3.根据权利要求2所述的起伏地表条件下的叠前深度偏移方法，其特征在于，所述的波场延拓算子其中， $\mathrm{\Δ}{k}_n\left(x\right)=w((1/v(x,z))-(1/v_0({\stackrel{\‾}{x}}_n,z)))$ 是局部扰动。

4.根据权利要求1所述的起伏地表条件下的叠前深度偏移方法，其特征在于，所述的步骤4包括：

对延拓波场进行如下滤波处理：u_cor(x,z+Δz)=u_cor(x,z+Δz)filter(x,z+Δz)，滤波函数filter(x,z+Δz)为：

$\mathrm{filter}(x,z+\mathrm{\Δz})=\left\{\begin{array}{c}0,(x,z+\mathrm{\Δz})\∈\mathrm{Fictitious\; layer}\\ 1,(x,z+\mathrm{\Δz})\∈\mathrm{Others}\end{array}\right..$

5.根据权利要求1所述的起伏地表条件下的叠前深度偏移方法，其特征在于，所述的步骤5包括：

根据 $R(x,z+\mathrm{\Δz})=\frac{1}{2\mathrm{\π}}\text{\∫}{u}_{\mathrm{con}}^r(x,z+\mathrm{\Δz})u_{\mathrm{con}}^{s*}(x,z+\mathrm{\Δz})\mathrm{dw}$ 对震源波场

和记录波场

进行互相关成像生成z+Δz层的成像结果，其中，

是震源波场的复共轭。

Images