CN103926619B - 一种三维vsp数据的逆时偏移方法 - Google Patents

Abstract

本方法公开一种三维VSP数据的逆时偏移方法，其为基于三维声波方程进行三维VSP资料逆时偏移成像的方法。本发明可充分利用VSP数据中丰富的波场信息和逆时偏移方法对复杂波场的成像能力，实现对井周复杂构造的准确成像，为复杂油气藏勘探提供有利指导。该方法充分考虑了地下实际存在的不同波场信息(多次波，绕射波等)，对波动方程不做任何简化。该方法利用优化的有限差分系数对偏微分方程进行中心差分求解，与传统相比，精度更高，计算更稳定。将该方法应用于三维VSP地震资料处理，对井旁小断层成像、储层预测、微幅构造落实以及隐蔽性油气藏勘探与开发都具有重要意义。

Description

一种三维VSP数据的逆时偏移方法

技术领域

本发明属于地球物理勘探领域，涉及一种三维VSP数据的逆时偏移方法，具体是采用逆时偏移方法对三维VSP地震资料进行成像处理，从而实现地下构造的准确成像。

背景技术

在地震资料处理中，偏移始终是处理环节中的关键部分，偏移的准确与否也影响着地震解释等后续工作的准确程度。为了提高地下复杂地质体的成像效果，能够准确识别地下储层，地震资料处理技术需要不断的发展以应对实际需求。逆时偏移是各种现行偏移算法中成像精度最高的偏移方法，方法本身基于双程波的波动方程，几乎不对方程做任何假设，因此具有较强的适应性。并且理论上能够对地下存在的各种波准确成像，对复杂构造也有较好的成像能力，是地震偏移理论中的一个重要研究领域，并也处于不断的发展当中。逆时偏移的概念最早在1983年被提出，当时Whitmore(1983)、Baysal等(1983)阐述了不同的逆时偏移概念；Loewenthal和Mufti(1983)将逆时偏移应用于空间-频率域；Levin等(1984)概括了逆时偏移的基本原理及实现方法；Hildebrand等(1987)将逆时偏移方法应用于波阻抗成像，取得了较好的效果；Zhu和Lines(1994)比较了逆时偏移与克希霍夫积分偏移的成像效果，得出前者对Marmous模型成像精度更高，但耗时较多的结论；Wu等(1996)对三维高阶有限差分法逆时偏移技术进行了研究；随着地球物理学家们从不同的角度对其进行着研究，逆时偏移理论也在不断的被完善。Faqi Liu等(2007)介绍了一种新的成像条件，通过分解波场，使成像过程中能够进行去噪，但增加了计算，并且会损失部分有效信息；Sandip等(2008)对逆时偏移的成像条件做了比较系统和全面的介绍；Yu Zhang等(2009)提出拉普拉斯算子滤波，相当于对成像波场的角度域衰减；Robert等(2009)将随机边界思想引入逆时偏移很好的解决了波场存储问题；陈生昌等(2012)提出了一种三维VSP数据高效偏移成像的超道集方法，极大地提高了计算效率；刘守伟和王华忠(2012)通过修改波场延拓公式，提出了上下行反射波联合成像方法，获得更宽的成像范围和更高的成像精度；Jinhai Zhang(2013)对差分系数进行优化，取得了较好的模拟精度。贺茜君和杨顶辉等(2014)通过波场模拟计算，对几种数值模拟方法的数值频散进行了细致的对比分析。

当今油气勘探目标日趋复杂，传统的一些构造简单的油气藏逐渐减少，有效解决各种隐蔽性油藏、陡倾角构造、横向变速地层、薄储层及小断层等非均质复杂油气藏问题当前也是未来一段时间勘探地球物理的研究目标。上述提及的逆时偏移方法虽从不同角度对理论及关键技术进行研究与分析，但对三维VSP数据，几乎没有逆时偏移成像的研究。

发明内容

针对上述技术难题，本发明人经过多年研究，针对传统偏移方法成像VSP数据的不足之处，经过多次设计和研究，提出一种三维VSP数据的逆时偏移方法。

依据本发明的技术方案，三维VSP数据的逆时偏移方法包括以下步骤：

1)在井中垂向等间隔布置检波器，地表等间隔平行设置炮线，井中检波器记录人工震源激发的每炮反射地震数据并将其记录到磁带上，利用常规预处理方法完成叠前地震数据的噪音衰减、振幅补偿、反褶积等处理，利用克希霍夫叠前深度偏移方法获得深度域层速度模型；

2)应用如下的计算公式作为逆时偏移中的震源子波：

$f\left(t\right)=[1-2{\left(\mathrm{\π}{f}_{p}t\right)}^2]e^{-{\left(\mathrm{\π}{f}_{p}t\right)}^2}$

式中，f_p为主频，单位为Hz，实际计算时，分析地震资料的频带宽度，如果地震资料的最高频率是f_max，则其主频f_p取为f_max/2.0；t为时间，单位是s；

3)读取震源位置信息、速度信息，由成像精度要求及输入的震源子波的主频确定偏移网格大小，考虑下面的波动方程求解的稳定性条件，

$\frac{\mathrm{v\&Delta;t}}{\mathrm{\&Delta;d}}<1/\sqrt{2\underset{m=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{b}_{2m-1}}$

式中，空间网格步长Δx＝Δz＝Δd，N为差分阶数，b为差分系数，N₁为不超过N的最大奇数；Δt为时间步长，单位为s。由上式可知，已知网格大小和读入的速度模型可计算所需要的时间步长，速度确定后，当网格变小时，时间步长也要相应变小。

4)根据速度模型确定粗密网格区域划分。在数值模拟计算中，采用粗密变网格划分方法，对于低速带区以及构造复杂的区域，利用加密网格计算，其他地区采用粗网格进行计算，兼顾考虑了计算效率和计算精度；

5)基于三维声波方程，选择利用全局优化方法确定的优化的有限差分系数，沿时间正向进行波场数值模拟，在不引入额外计算的前提下，与同阶的原始差分系数相比，优化的差分系数可有效的降低频散，提高地震波数值模拟计算精度；

6)为提高计算精度，设置完全匹配层边界条件，并保存每一时刻有效波场与完全匹配层之间交界处波场信息，节省大部分存储空间；

7)读入检波点位置信息和叠前单炮地震数据；

8)同地震波场的正向传播，从三维声波方程出发，采用优化的有限差分系数和完全匹配层边界条件，沿时间轴反传检波点波场，在反传过程中拾取每时刻保存的震源波场信息，用其替代有效波场与完全匹配层边界交界处的反传波场信息；

9)对同时刻的震源波场和检波点波场采用互相关成像条件，并对所有时刻的成像值累加；

10)对成像结果采用拉普拉斯方法进行低频噪音的压制处理，得到高信噪比和高分辨率的井周围3D VSP地震成像数据体，进而提高储层描述的精度，对更精细地刻画井周围裂缝等非均质地层及其深部构造有着重要的意义，推动油藏地球物理的发展。

其中，步骤2)所设计的震源子波为雷克子波。步骤3)所给出的确定偏移网格大小及时间步长，是指依据所需要的成像精度(即所需要识别的构造尺度)以及满足地震波采样的原则确定偏移的网格大小，并根据有限差分的稳定性条件由网格大小和读入的速度模型信息计算得到所需要的时间步长。步骤4)给出的根据速度模型确定粗密网格区域划分，是指首先对速度模型进行扫描，将速度模型划分出低速区域与非低速区，对低速部分进行网格加密。并进一步分析地下大致构造形态，对形态复杂的区域同样进行加密网格。其余区域采用粗网格剖分，最后依据各区域的网格大小及所对应的速度由稳定性条件确定各个区域所需要的时间步长。

步骤5)所述的波场模拟过程中，从三维声波方程出发，采用高阶有限差分数值算法，得到离散的有限差分方程，对三维地震波场，计算每个时刻沿时间正向的波场，在粗密网格区域，采用优化的差分系数进行波场计算，而在过渡区域，采用过渡区域的差分系数进行计算。在每时刻的计算过程中，需要保存有效区域与边界区域上六个面的波场信息。在步骤8)的波场反传过程中，在每时刻反传时，用之前保存的相应时刻的六个面上的波场信息进行替换，以保证波场反传的精度。

优选地，步骤5)、8)所采用的三维声波方程为：

$\frac{{\&PartialD;}^{2}w}{\&PartialD;t^2}=v^2(\frac{{\&PartialD;}^{2}w}{\&PartialD;x^2}+\frac{{\&PartialD;}^{2}w}{\&PartialD;y^2}+\frac{{\&PartialD;}^{2}w}{\&PartialD;z^2})$

式中：w为波场值，是空间坐标(x，y，z)和时间t的函数，t为时间，单位是s，v是相应空间位置的速度，单位为m/s。

更优选地，步骤5)、8)的波场数值模拟过程中采用优化的高阶有限差分系数。

进一步地，步骤5)、8)的波场数值模拟过程中依据步骤3)的粗密网格区域划分，对粗密网格部分采用常规的差分计算式，对粗密网格交叉部分采用过渡差分计算式。

更近一步地，步骤6)、8)的波场数值模拟过程中，保存每一时刻的有效波场与完全匹配层之间交界处波场信息，并在步骤7)中的每时刻进行应用。

与传统的VSP数据成像方法相比，本发明公开的三维VSP逆时偏移方法具有如下优势：

1)本方法对三维声波方程不做任何近似，直接求解，因此对于地下实际介质具有更强的适应性，可充分利用实际存在的各种地震波信息，实现高信息量地震数据成像。

2)本发明在逆时偏移的波场模拟过程中采用优化差分系数，在相同差分阶数的情况下，减少了采用常规差分系数时出现的频散，明显提高数值模拟精度。

3)本方法采用的变网格变时间步长可明显改善由全局细网格小时间步长带来的计算效率低的状况。

4)传统逆时偏移算法应用完全匹配层边界条件时，需要保存每一时刻的全部波场信息，随着所处理数据尺度的增大，这种存储方式将使计算机硬件无法负荷，与其有所不同，本发明通过记录每一时刻有效区域与边界区域交界部分的波场信息来减少存储，并将每一时刻保存的波场信息用于震源波场反传，且几乎不损失成像精度。

附图说明

图1-1为均匀介质中采用传统差分系数正演剖面；图1-2为采用优化差分系数正演剖面；两图中除差分系数不同外，其余参数全部相同；

图2为变网格差分原理示意图；

图3为变时间步长原理示意图；

图4为三维VSP逆时偏移方法流程图；

图5为三维盐丘模型示意图；

图6-1为三维盐丘模型的一水平切片示意图；

图6-2为成像体所对应的水平切片示意图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。另外地，不应当将本发明的保护范围仅仅限制至下述具体步骤或具体参数。本发明不局限于模型数据，可针对实际资料进行处理，具有广泛的适应性。

本发明公开的一种三维VSP数据的逆时偏移方法包括以下步骤：

1)在井中垂向等间隔布置检波器，地表等间隔平行设置炮线，井中检波器记录人工震源激发的每炮反射地震数据，并将其记录到磁带上，利用常规预处理方法完成叠前地震数据的噪音衰减、振幅补偿、反褶积等处理，利用克希霍夫叠前深度偏移方法获得深度域层速度模型；

2)应用如下的计算公式作为逆时偏移中的震源子波：

$f\left(t\right)=[1-2{\left(\mathrm{\&pi;}{f}_{p}t\right)}^2]e^{-{\left(\mathrm{\&pi;}{f}_{p}t\right)}^2}$

式中，f_p为主频，单位为Hz，实际计算时，分析地震资料的频带宽度，如果地震资料的最高频率是f_max，则其主频f_p取为f_max/2.0；t为时间，单位是s；

3)读取震源位置信息、速度信息，由成像精度要求及输入的震源子波的主频确定偏移网格大小，考虑下面的波动方程求解的稳定性条件，

$\frac{\mathrm{v\&Delta;t}}{\mathrm{\&Delta;d}}<1/\sqrt{2\underset{m=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{b}_{2m-1}}$

式中，空间网格步长Δx＝Δz＝Δd，N为差分阶数，b为差分系数，N₁为不超过N的最大奇数；Δt为时间步长，单位是s。由上式可知，已知网格大小和读入的速度模型可计算所需要的时间步长，速度确定后，当网格变小时，时间步长也要相应变小。

4)根据速度模型确定粗密网格区域划分。在数值模拟计算中，采用粗密变网格划分方法，对于低速带区以及构造复杂的区域，利用加密网格计算，其他地区采用粗网格进行计算，兼顾考虑了计算效率和计算精度；

5)基于三维声波方程，选择利用全局优化方法确定的优化的有限差分系数，沿时间正向进行波场数值模拟，在不引入额外计算的前提下，与同阶的原始差分系数相比，优化的差分系数可有效的降低频散，提高地震波数值模拟计算精度；

6)为提高计算精度，设置完全匹配层边界条件，并保存每一时刻有效波场与完全匹配层之间交界处波场信息，节省大部分存储空间；

7)读入检波点位置信息和叠前单炮地震数据；

8)同地震波场的正向传播，从三维声波方程出发，采用优化的有限差分系数和完全匹配层边界条件，沿时间轴反传检波点波场，在反传过程中拾取每时刻保存的震源波场信息，用其替代有效波场与完全匹配层边界交界处的反传波场信息；

9)对同时刻的震源波场和检波点波场采用互相关成像条件，并对所有时刻的成像值累加。

10)对成像结果采用拉普拉斯方法进行低频噪音的压制处理，得到高信噪比和高分辨率的井周围3D VSP地震成像数据体，进而提高储层描述的精度，对更精细地刻画井周围裂缝等非均质地层及其深部构造有着重要的意义，推动油藏地球物理的发展。

在本发明中，步骤(3)所述的确定偏移网格大小及时间步长，是指根据输入的速度模型、稳定性条件以及实际对构造的成像效果是否精细决定偏移的网格大小及时间步长。

在本发明中，步骤(4)所述的根据速度模型确定粗密网格区域划分，是指根据速度模型判断地下的大致构造形态，对于形态复杂的构造采用细密网格和小时间步长，其他区域采用粗网格和大时间步长。

步骤(5)所述的波场模拟过程中，在每一时刻，对于粗网格和密网格区域采用优化的有限差分系数，对于粗密网格交界的部分采用过渡区域差分系数，并保存有效区域与边界之间交界部分的波场信息。并应用到步骤(8)的震源波场反传中。

在本发明所提供的三维VSP逆时偏移方法中，其波场计算方法如下：

三维声波方程的形式如下：

$\frac{{\&PartialD;}^{2}w}{\&PartialD;t^2}=v^2(\frac{{\&PartialD;}^{2}w}{\&PartialD;x^2}+\frac{{\&PartialD;}^{2}w}{\&PartialD;y^2}+\frac{{\&PartialD;}^{2}w}{\&PartialD;z^2})$

式中，w为波场值，是空间坐标(x，y，z)和时间t的函数，t为时间，单位是s；v是相应空间位置的速度，单位为m/s。对二阶偏导数采用2N阶中心有限差分法进行展开得到：

$\frac{{\&PartialD;}^{2}w}{\&PartialD;x^2}\&ap;\frac{1}{\mathrm{\&Delta;}{x}^2}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{a}_i*[w(x+\mathrm{i\&Delta;x})+w(x-\mathrm{i\&Delta;x})-2w\left(x\right)]$

式中，Δx为网格大小，a_i为有限差分系数。

差分系数需求解以下矩阵：

式中，a₁、a₂和a_N等为差分系数，N为差分阶数。

求解得：

$a_i=\frac{{(-1)}^{i+1}}{i^2}\frac{{\mathrm{\&Pi;}}_{\underset{j\&NotEqual;i}{j=1}}^N{j}^2}{{\mathrm{\&Pi;}}_{\underset{j\&NotEqual;i}{j=1}}^N|i^2-j^2|}i=\mathrm{1,2}\&CenterDot;\&CenterDot;\&CenterDot;N$

而优化的有限差分系数通过最优化以下目标函数求得：

$E(k_x^{\max },T)\&equiv;\underset{0\&le;k_x\&le;k_x}{{\max }_{\max }}|-k_x^2{\mathrm{\&Delta;}}^2-\underset{n=-N-2}{\overset{N/2}{\mathrm{\&Sigma;}}}{b}_n\cos \left(n{k}_x\mathrm{\&Delta;}\right)|\&le;T$

式中，k_x为与空间偏导数相对应的解析波数，为最大波数，T为误差限值，Δ为网格步长，N为阶数，b_n为优化系数。

优化的有限差分系数与传统差分系数相比，可明显改善波场模拟过程中的数值频散，提高数值模拟的精度，如图1-1和图1-2所示。

经上述差分方式展开后，在计算过程中，Δx、Δz空间网格步长与Δt时间步长需满足稳定性条件，其形式如下：

$\frac{\mathrm{v\&Delta;t}}{\mathrm{\&Delta;d}}<1/\sqrt{2\underset{m=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{b}_{2m-1}}$

式中，Δx＝Δz＝Δd，N₁为不超不N的最大奇数。由上式可知，当网格变小时，时间步长也要相应变小。

为提高计算效率，本发明采用的变网格方法如下：

假设加密区域x、z方向的网格大小分别为Δx、Δz，非加密区域x、z方向的网格大小分别为ΔX、ΔZ(ΔX与Δx之比为L)。

加密区域和非加密区域的差分格式分别为：

$\frac{{\&PartialD;}^{2}w}{\&PartialD;x^2}\&ap;\frac{1}{\mathrm{\&Delta;}{x}^2}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{b}_i[w(x+\mathrm{i\&Delta;x},z)-2w(x,\mathrm{z\&Delta;})+w(x-\mathrm{i\&Delta;x},z)]$

$\frac{{\&PartialD;}^{2}w}{\&PartialD;x^2}\&ap;\frac{1}{\mathrm{\&Delta;}{X}^2}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{b}_i[w(x+\mathrm{i\&Delta;X},z)-2w(x,\mathrm{z\&Delta;})+w(x-\mathrm{i\&Delta;X},z)]$

式中：b_i为差分系数，2N为空间差分精度。

加密区与非加密区交界部分的差分格式为

$\begin{array}{c}\frac{{\&PartialD;}^{2}w}{\&PartialD;x^2}\&ap;\frac{1}{\mathrm{\&Delta;}{x}^2}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{b}_i^{\&prime;}[w(x+\mathrm{i\&Delta;x},z)-2w(x,z)+w(x-\mathrm{i\&Delta;x},z)]\\ +\frac{1}{\mathrm{\&Delta;}{X}^2}\underset{i=n+1}{\overset{N}{\mathrm{\&Sigma;}}}{b}_i^{\&prime;}[w(x+\mathrm{i\&Delta;X},z)-2w(x,z)+w(x-\mathrm{i\&Delta;X},z)]\end{array}$

式中：bi′为过渡区差分系数，n为加密区中所在网格点与加密区边界之间的网格点数，其原理如图2所示(以N＝2，K＝4为例)。

由稳定性条件可知，密网格需小时间步长，为解决全局小时间步长带来的计算效率降低的情况，本发明采用变时间步长算法，其原理如图3所示(以大小时间步长相差2倍为例)。

图中每个切片(wf1、wf2、wf3、wf4、wf5)为对应时刻的波场，对于大时间步长区域，已知wf1和wf3依据常规大网格差分格式、ΔT得出wf5；在过渡区域，已知wf1和wf3依据过渡区差分格式、ΔT得出wf5，已知wf2和wf3依据过渡区差分格式、Δt得出wf4；已知wf2和wf3依据常规细差分格式、Δt得出wf4，已知wf3和wf4依据常规细网格差分格式、Δt得出wf5。

完全匹配层边界条件可以几乎完美的吸收由计算机有限存储所带来的边界反射，其形式如下：

$\frac{{\&PartialD;}^2{w}_x}{\&PartialD;t^2}+2d_x\left(x\right)\frac{\&PartialD;w_x}{\&PartialD;t}+d_x^2\left(x\right)w_x=v^2\frac{{\&PartialD;}^{2}w}{\&PartialD;x^2}$

$\frac{{\&PartialD;}^2{w}_z}{\&PartialD;t^2}+2d_z\left(z\right)\frac{\&PartialD;w_z}{\&PartialD;t}+d_z^2\left(z\right)w_z=v^2\frac{{\&PartialD;}^{2}w}{\&PartialD;z^2}$

式中，d_x(x)、d_z(z)为沿x方向的衰减函数，v是速度，单位为m/s。

为减少应用完全匹配层边界条件所耗费的大量存储，本发明采用保存每时刻有效区域与边界之间的交界部分，从而节省存储。

在本发明中，逆时偏移过程中包含大量的波场数值计算，有限差分算法兼顾效率与精度，是目前较为流行的方法。利用有限差分法对波动方程进行数值求解时，在阶数固定的前提下，差分系数对数值模拟的精度起着决定性作用。优化的差分系数可以在相同阶数的前提下，比利用传统差分系数获得更高的精度，且能有效减少因求解不稳定带来的频散现象。将优化差分系数引入到三维波场模拟中可以明显的提高算法的稳定性及成像精度。

计算和存储是影响逆时偏移发展的两个主要因素。在正演模拟中，为在不引入过多计算量的前提下，采用变网格变时间步长技术模拟地下小的地质体。只需在需要的区域加密网格就可得到与全局细网格相类似的响应，将其引入到三维逆时偏移中可以有效的节省大量计算，增强算法的实用性。

本发明所述的三维VSP逆时偏移方法充分考虑了VSP地震记录中丰富的波场信息，通过逆时偏移方法对复杂波场的成像能力进行偏移处理，能够充分发挥三维VSP数据高分辨率、高信噪比的优势，对井旁复杂地质体成像更清晰。下面参考附图，对本发明做出进一步说明。

参考附图1，两图显示的均为均匀介质内部震源激发后的波场快照，其中两图选用相同频率的雷克子波、网格大小、时间步长等参数，仅差分系数的选取不同，其中图1-1中使用传统差分系数，图1-2中使用优化差分系数，对比可发现图1-1中明显的数值频散在附图1-2中得到明显的减弱，优化差分系数在相同阶数的前提下可明显减弱使用传统差分系数带来的频散现象，有效提高数值模拟的精度，改善成像效果。在实际操作过程中，只需用优化差分系数取代传统差分系数即可，并不会产生额外的计算及存储。

参考附图2、3说明了变网格变时间步长的基本原理，从附图2中可以看出对于位于粗网格和密网格的区域，差分方式原理与传统格式并无区别，仅就粗密网格交界部分有所差异，过渡区域的宽度与所采用的差分阶数有关。当空间采用2N阶差分时，过渡层则为N-1层，因此在处理这N-1层时需要使用过渡区域的差分公式和差分系数，每一层的差分系数均不相同。并且对于处理不在网格点上的点时，需要插值处理。附图3表示的是变时间步长原理，同变网格类似，粗网格和密网格区域都分别采用大时间步长和小时间步长进行运算，对于过渡区域，则需要大时间步长、小时间步长交替使用，才能保证过渡平稳。

参考附图4，公开了本发明三维VSP逆时偏移方法流程图，包括下述步骤：1)设计震源子波；2)读取震源位置信息、速度信息，确定偏移网格大小及时间步长；3)根据速度模型确定粗密网格区域划分；4)采用优化的有限差分系数沿时间正向进行波场数值模拟，并设置完全匹配层边界条件；5)保存每一时刻有效波场与完全匹配层边界交界处波场信息；6)读入检波点位置信息和单炮数据；7)沿时间轴反向反传检波点波场，依旧采用优化查分系数及完全匹配层边界；8)同时利用每时刻保存的震源波场信息进行反传；9)对同时刻的震源波场和检波点波场采用互相关成像条件，并对所有时刻的成像值累加；10)对成像结果采用拉普拉斯方法进行低频噪音的压制处理，输出最终成像结果。

参考附图5所示，三维盐丘模型(SEG/EAGE Salt Model)是SEG所发布，其本身具有复杂的构造，是墨西哥湾盐丘构造的代表，常用于正演地震资料及偏移测试。

参考附图6是利用本发明的三维VSP逆时偏移方法对三维盐丘模型测试的多炮偏移叠加结果。对于附图6中的a、b两水平切片，可以看出本发明所采用的方法成像准确且较为清晰，可为实际三维VSP勘探提供有利指导。

概述地，本方法公开的三维VSP数据的逆时偏移方法为基于三维声波方程进行三维VSP资料逆时偏移成像的方法。本发明可充分利用VSP数据中丰富的波场信息和逆时偏移方法对复杂波场的成像能力，实现对井周复杂构造的准确成像，为复杂油气藏勘探提供有利指导。该方法充分考虑了地下实际存在的不同波场信息(多次波，绕射波等)，对波动方程不做任何简化，因而算法更符合实际，从而能够实现多种波的准确成像。该方法利用优化的有限差分系数对偏微分方程进行中心差分求解，与传统相比，精度更高，计算更稳定。算法考虑了实际的网格剖分与时间步长关系，应用变网格变时间步长方法解决地下局部构造成像所带来的计算精度和计算效率问题。在边界条件方面，方法采用设置检测层的完全匹配层边界条件弥补使用随机边界给成像结果引入的随机噪音问题。将该方法应用于三维VSP地震资料处理，对井旁小断层成像、储层预测、微幅构造落实以及隐蔽性油气藏勘探与开发都具有重要意义。

如上述，已经清楚详细地描述了本发明提出的一种三维VSP数据的逆时偏移方法。尽管本发明的优选实施例详细描述并解释了本发明，但是本领域普通的技术人员可以理解，在不背离所附权利要求定义的本发明的精神和范围的情况下，可以在形式和细节中做出多种修改。

Claims (7)

1.一种三维VSP数据的逆时偏移方法，其包括以下步骤： 

1)在井中垂向等间隔布置检波器，在地表等间隔平行设置炮线，井中检波器记录人工震源激发的每炮反射地震数据并将其记录到磁带上，利用常规预处理方法完成叠前地震数据的噪音衰减、振幅补偿、反褶积处理，利用克希霍夫叠前深度偏移方法获得深度域层速度模型； 

2)应用如下的计算公式作为逆时偏移中的震源子波： 

式中，f_p为主频，单位为Hz，实际计算时，分析地震资料的频带宽度，如果地震资料的最高频率是f_max，则其主频f_p取为f_max/2.0；t为时间，单位是s； 

3)读取震源位置信息、速度信息，根据成像精度要求及输入的震源子波的主频确定偏移网格大小，考虑下面的波动方程求解的稳定性条件， 

式中，空间网格步长Δx＝Δz＝Δd，N为差分阶数，b为差分系数，N₁为不超过N的最大奇数；Δt为时间步长，单位为s，由上式可知，已知网格大小和读入的速度模型可计算所需要的时间步长，速度确定后，当网格变小时，时间步长也要相应变小； 

4)根据速度模型确定粗密网格区域划分：在数值模拟计算中，采用粗密变网格划分方法，对于低速带区以及构造复杂的区域，利用加密网格计算，其他地区采用粗网格进行计算，兼顾考虑了计算效率和计算精度； 

5)基于三维声波方程，选择利用全局优化方法确定的优化的有限差分系数，沿时间正向进行波场数值模拟，在不引入额外计算的前提下，与同阶的原始差分系数相比，优化的差分系数可有效的降低频散，提高地震波数值模拟计算精度； 

6)为提高计算精度，设置完全匹配层边界条件，并保存每一时刻有效波场与完全匹配层之间交界处波场信息，节省大部分存储空间； 

7)读入检波点位置信息和叠前单炮地震数据； 

8)同地震波场的正向传播，从三维声波方程出发，采用优化的有限差分系数和完全匹配层边界条件，沿时间轴反传检波点波场，在反传过程中拾取每时刻保存的震源波场信息，用其替代有效波场与完全匹配层边界交界处的反传波场信息； 

9)对同时刻的震源波场和检波点波场采用互相关成像条件，并对所有时刻的成像值累 加； 

10)对成像结果采用拉普拉斯方法进行低频噪音的压制处理，得到高信噪比和高分辨率的井周围3D VSP地震成像数据体，进而用于提高储层描述的精度和精细地刻画井周围裂缝的非均质地层。 

2.根据权利要求1的方法，特征是：步骤2)所设计的震源子波为雷克子波。 

3.根据权利要求1的方法，特征是：步骤4)给出的根据速度模型确定粗密网格区域划分，是指首先对速度模型进行扫描，将速度模型划分出低速区域与非低速区，对低速部分进行网格加密 ，并进一步分析地下大致构造形态，对形态复杂的区域同样进行加密网格 ，其余区域采用粗网格剖分，最后依据各区域的网格大小及所对应的速度由稳定性条件确定各个区域所需要的时间步长。 

4.根据权利要求1的方法，特征是：步骤5)和步骤8)所采用的三维声波方程为： 

式中：w为波场值，是空间坐标(x，y，z)和时间t的函数，t为时间，单位是s，v是相应空间位置的速度，单位为m/s。 

5.根据权利要求1的方法，特征是：步骤5)和步骤8)的波场数值模拟过程中采用优化的高阶有限差分系数。 

6.根据权利要求1的方法，特征是：步骤5)和步骤8)的波场数值模拟过程中依据步骤3)的粗密网格区域划分，对粗密网格部分采用常规的差分计算式，对粗密网格交叉部分采用过渡差分计算式。 

7.根据权利要求1的方法，特征是：步骤6)和步骤8)的波场数值模拟过程中，保存每一时刻的有效波场与完全匹配层之间交界处波场信息，并在步骤7)中的每时刻进行应用。