三维保幅叠前时间偏移方法

摘要

三维保幅叠前时间偏移方法，应用于地震勘探中反射地震资料处理，是针对三维采集地震数据的叠前偏移成像方法。采取的技术方案包括三方面：一是采用基于单程波方程和深度偏移的反褶积成像条件得到的权系数进行偏移幅值计算；二是根据地下构造的最大倾角确定时变的三维偏移孔径；三是通过对偏移孔径边缘引入带状衰减，自动压制偏移噪音。本发明通过改进偏移孔径、给出保幅的权系数、在成像过程中压制偏移噪音，实现了高效、高信噪比、保幅偏移成像。本发明生成的保幅共反射点道集可更好地服务于叠前反演等油气和流体检测技术。本发明对油气、矿产资源勘探有重要应用价值。

说明书

技术领域

本发明属于油气田地震勘探中反射地震资料处理技术领域，尤其涉及地震资料处理过程中的叠前偏移成像技术范畴，是三维地震资料偏移成像的三维保幅叠前时间偏移方法。

背景技术

地震勘探中反射地震资料处理过程中，叠前偏移成像是关键的环节，而叠前时间偏移方法是叠前偏移成像中的一种重要方法。叠前时间偏移方法可对一类断层较为复杂但速度横向变化不是很剧烈的地质构造较好成像。与叠前深度偏移方法相比，除具有较高的计算效率外，其主要的优点是只需使用叠加(均方根)速度；这样可简单地通过速度扫描等方式得到恰当的偏移速度模型，回避了使用叠前深度偏移方法面临的一个主要困难：速度建模。因此，叠前时间偏移方法已成为地震勘探领域广泛应用的关键技术。

影响叠前时间偏移成像效果的因素包括：偏移速度、地震波走时计算方法、偏移计算的偏移孔径、计算偏移幅值时的权系数、偏移算法实现流程。走时计算与偏移速度共同决定了反射波能否正确归位，偏移孔径及其应用方式决定了偏移噪音和偏移算法的计算量，权系数决定了成像幅值能否正确反应地下实际界面的物性参数变化，偏移算法实现流程对偏移的计算效率和存储需求有重要影响。对偏移方法而言，成像效果、计算效率和存储需求是评价偏移方法的重要指标。

偏移孔径对叠前时间偏移是重要的，较好的偏移孔径对压制偏移噪音和降低偏移计算量有更重要作用。小的偏移孔径可减少偏移计算量，但存在着不能对陡倾角构造正确成像的风险；过大的孔径又带来了偏移噪音和较大的计算量。对三维叠前时间偏移方法而言，由于采集数据在沿测线和垂直测线方向有不同的空间采样率，偏移孔径的选取较二维偏移更加复杂。

现行的叠前时间偏移方法是以双平方根方程为基础的，其成像条件实际上使用了叠前深度偏移的相关成像条件，因此成像时没有补偿地震波传播的几何扩散效应。为获得保幅的共反射点(CRP)道集，以服务于叠前反演等利用叠前信息进行油气和流体检测，需发展保幅，即幅值能正确反应地下界面物性参数变化的叠前时间偏移方法。

三维观测已成为海上和陆上地震勘探的主要观测方式，三维观测和三维偏移可对陡倾角构造正确成像，避免了二维观测和偏移时侧反射难以正确归位的问题。随着采集和计算设备的发展，三维采集和偏移处理已变的更容易实现，因此必须发展相应的三维叠前时间偏移技术。

发明内容

针对现有情况，本发明的目的是：提供一种三维保幅叠前时间偏移方法，通过改进偏移孔径、给出保幅的权系数、在成像过程中压制偏移噪音，实现高效、高信噪比、保幅偏移成像。这一方法可获得保幅的共反射点(CRP)道集，从而更好地服务于叠前反演等利用叠前信息进行油气和流体检测等技术。

本发明采取以下技术方案：

(1)用多条拖缆或多条测线记录人工震源激发的反射地震信号，记录到磁带上；

(2)从磁带上读取地震资料，对叠前地震资料进行常规的压制噪音处理，针对几个共中心点位置，抽取共中心点道集，对这些道集作常规的动校正(NMO)速度拾取，对所得结果做横向平均，作为初始的、横向均匀的偏移速度场，将叠前地震资料按偏移距排序、分组，将不同组地震资料放置在集群计算机的不同计算节点上；

(3)根据地震资料的偏移距、方位角的变化范围和地下构造的最大倾角以及初始的、横向均匀的偏移速度场，构建时变的三维偏移孔径的偏移孔径表；

(4)根据偏移孔径表确定对应的压制偏移噪音的衰减系数表；

(5)采用基于单程波方程和深度偏移的反褶积成像条件得到的权系数进行偏移幅值计算；

(6)在每个计算节点，对全部地震道循环，对每一地震道，由偏移孔径表确定有效的成像区域和区域中各离散点偏移计算的起始成像时间，从起始成像时间开始计算各成像点的偏移幅值，对处于三维偏移孔径的边缘区中的成像点，查衰减系数表拾取衰减系数并乘以偏移幅值得到最终的幅值，将幅值累加到存放偏移结果的数组中对应的偏移距上；

(7)收集各计算节点的偏移计算结果，对成像区域的每一水平位置把偏移计算结果按时间深度、偏移距排序，形成各水平位置处的共反射点道集；

(8)依据共反射点道集中的同相轴是否平直，在偏移过程中确定偏移速度场；

(9)利用偏移过程中确定的偏移速度场做偏移计算，对形成的全部共反射点道集做剩余动校正，将出现明显拉伸和噪音部分的对应数值置为零(即切除)，将不同偏移距的数值叠加，形成偏移叠加剖面；

(10)通过显示软件将偏移叠加剖面数据转换为地下反射构造的剖面图像，剖面图像将指示地下构造的形态、断点、断距大小以及界面对地震波的反射强度。

所述的根据地震资料的偏移距、方位角的变化范围和地下构造的最大倾角，以及初始的横向均匀的偏移速度场、构建时变的三维偏移孔径的偏移孔径表是这样实现的：用一组数据(x_i，y_i，2T_i)(i＝1，k)描述三维偏移孔径，其中k是有效成像区域(成像空间在水平面的投影)中总的离散点数，(x_i，y_i)是成像区域中某离散点与地震道中心点的两个水平坐标差，而T_i则是该点偏移计算的起始成像时间(单程旅行时)；令地震道的炮点坐标是(x_s，y_s)，检波点坐标是(x_r，y_r)，引入时变的偏移孔径，就是在该地震道的偏移计算中，对成像区域中水平坐标为(x_i+0.5(x_s+x_r)，y_i+0.5(y_s+y_r))的点(i＝1，k)，从T＝T_i开始进行偏移计算。

由构造沿测线和垂直测线方向的两个最大倾角，θ_x和θ_y，构建如下不等式方程组：

式中是地震道的方位角，满足上两式的最小角度α和β，就是保证最大倾角为θ_x和θ_y的构造得以正确成像时，沿方位角方向和垂直方位角方向的两个最大成像角度。

对不同时间深度T，成像区域可看作以方位角方向为长轴方向、垂直方位角方向为短轴方向、中心在地震道的中心点上的椭圆面，由α、β和初始的、横向均匀的偏移速度场可求得其长、短半轴分别为

$a=\frac{{\bar{V}}_{\mathrm{rms}}\left(T\right)·T}{\sin \left(2\alpha \right)}\sqrt{1+{\cos }^2\left(2\alpha \right)+{\lambda }^2{\sin }^2\left(2\alpha \right)-2\cos \left(2\alpha \right){[1+{\lambda }^2{\sin }^2\left(2\alpha \right)]}^{1/2}}$

$b={\bar{V}}_{\mathrm{rms}}\left(T\right)·T\tan \beta$

式中参数λ定义为是地震道的半偏移距，时间深度T是单程旅行时。

对时间深度T和T+ΔT(ΔT是时间深度方向上成像的采样间距)，分别计算对应的长、短半轴，令这两组长、短半轴形成的两个椭圆有共同的中心点(0，0)，这两个椭圆形成的条带就是起始成像时间为T的成像区域了；将落在这个条带中的离散点坐标x＝mΔx，y＝nΔy和T构成一个数据对，其中Δx和Δy是偏移成像的水平坐标采样。对小于指定时间深度的时间深度采样点循环(大于该时间深度，偏移孔径将不随时间变化)，对得到的全部(x，y，2T)分别按x和y的大小排序，统计总点数k，就可以得到时变的三维偏移孔径。

地震道的方位角和偏移距决定了偏移孔径，需根据全部地震道的方位角和偏移距变化范围，依据角度和偏移距采样间距(为减少存储量，可取较大的间距)，计算随方位角和偏移距变化的三维偏移孔径，构建偏移孔径表；偏移孔径表是一个三维矩阵，一个维数对应偏移距，另一维数对应方位角，第三维是三维偏移孔径的三个坐标数据(x_i，y_i，2T_i)和n_i共4个数随i增加排列，这里n_i是对应衰减系数的一组整数；这一维的第一个元素是该偏移孔径的总点数k。实际上，该矩阵中存放的是坐标(x_i，y_i，2T_i)对应的离散点序号，因此偏移孔径表是一整数矩阵。

方位角的采样间距不是针对实际角度值，而是针对方位角的正弦值。拾取偏移孔径时，由地震道的偏移距和炮检点在垂直侧线方向的坐标差计算方位角的正弦值，对偏移距和方位角正弦值按采样间距取整，直接到表中拾取。

所述的根据偏移孔径表确定对应的压制偏移噪音的衰减系数表是这样实现的：每一偏移孔径对应的衰减系数数组由如下步骤得到：

1)定义三维矩阵(x，y，T)，依据从偏移孔径得到的(x_i，y_i)点的起始成像时间T_b，计算系数并将其存入对应(x_i，y_i，T)位置中，式中n是人为定义的边缘区包含的点数，ΔT是时间深度方向成像的采样间距(单程旅行时)，对偏移孔径中全部k个点循环；

2)对每个时间深度T，由下式计算a和b：

$a=\frac{{\bar{V}}_{\mathrm{rms}}\left(T\right)·T}{\sin \left(2\alpha \right)}\sqrt{1+{\cos }^2\left(2\alpha \right)+{\lambda }^2{\sin }^2\left(2\alpha \right)-2\cos \left(2\alpha \right){[1+{\lambda }^2{\sin }^2\left(2\alpha \right)]}^{1/2}}$

$b={\bar{V}}_{\mathrm{rms}}\left(T\right)·T\tan \beta$

式中各参数与计算偏移孔径时的相同，由最大的水平方向成像采样间距Δy，计算进而计算其中是地震道的方位角。对f＞1的(x_i，y_i，T)，若对应的(x_i，y_i，T)位置中已存储了c_i，不做操作，否则计算将c₂存入对应的位置中；

3)对三维矩阵(x，y，T)的各个元素循环，将在前两个步骤中没赋值的点赋值1.0；对矩阵元素做空间高斯平滑，平滑后仅保留小于1.0的元素，由这些元素形成衰减系数数组。偏移孔径数组的每组(x_i，y_i，2T_i)在衰减系数数组中对应一个序列c_j(j＝1，n_i)，其中n_i是矩阵(x，y，T)中(x_i，y_i)处、沿T_b起小于1.0的元素个数，c_j是这些小于1.0的元素值，如有关偏移孔径表中讨论的那样，n_i将存放在偏移孔径表中。

衰减系数表是与偏移孔径表对应，偏移孔径表的每个偏移孔径在衰减系数表中对应一个衰减系数数组。衰减系数表也是一个三维矩阵，一个维数对应偏移距，另一维数对应方位角，第三维是序列c_j(j＝1，n_i)随i增加排列，其中i＝1，k，k是偏移孔径表中相同偏移距和方位角处的第一个元素值。拾取衰减系数时，对地震道的偏移距和方位角(其正弦值)按采样间距取整，直接到表中拾取。地震道的方位角和偏移距决定了偏移孔径表，也决定了衰减系数表，两者采用同样方法存放和拾取。

所述的采用基于单程波方程和深度偏移的反褶积成像条件得到的权系数进行偏移幅值计算是这样实现的：定义

$p_s=\frac{{(x-x_s)}^2+{(y-y_s)}^2}{V_{\mathrm{rms}}^2{T}^2},$$p_r=\frac{{(x-x_r)}^2+{(y-y_r)}^2}{V_{\mathrm{rms}}^2{T}^2}$

式中(x_s，y_s)是地震道的炮点坐标，(x_r，y_r)是检波点坐标，(x，y，T)是成像点的水平和时间深度(单程旅行时)坐标，V_rms是成像点处的偏移速度，则炮点至成像点再到检波点的地震波走时权系数为

对地震道应用FFT(快速傅立叶变换)，在有效频带内对每个频率分量乘-jω，其中j是单位虚数，ω是频率；令地震道记录地震信号的时间采样为Δt，在频率域将频率上限取为2/Δt，将超过实际频率上限部分的值置为零；应用FFT做傅立叶反变换，此时得到的地震信号的时间采样变为0.25Δt；对4τ/Δt-0.001取整得到整数i，由变换后的时域数据上i和i+1两点的值做线性插值，得到时间τ处的值a，则偏移幅值取为

通过上述方案描述可知，本发明的核心有三点：一是从三维叠前深度偏移方法出发分析和研究三维叠前时间偏移，基于单程波理论和稳相点原理，得到三维非均匀介质中的地震波走时和幅值，然后根据叠前深度偏移的反褶积成像条件，得到了保幅的权系数；二是从三维偏移脉冲响应的倾角分析出发，给出由沿测线和垂直测线方向的两个最大成像角度(构造倾角)决定的时变三维偏移孔径；三是通过对偏移孔径边缘引入带状衰减，自动压制偏移噪音。其具体实现原理如下：

1.地震波走时、幅值和权系数

假设三维非均匀介质可近似为层状介质，在波数一频率域，单个炮点或检波点的地震波场的深度延拓可用旅行时(时间深度)表示为：

$\tilde{P}(p_x,p_y,\omega ,T=\underset{i=1}{\overset{n}{\Sigma }}{\mathrm{\Delta T}}_i)=f\left(\omega \right)\exp (-\mathrm{j\omega }\underset{i=1}{\overset{n}{\Sigma }}{\mathrm{\Delta T}}_i\sqrt{1-v_i^2(p_x^2+p_y^2)})---\left(1\right)$

式中ΔT_i是各层介质用单程旅行时表达的厚度，n是某一深度包含的层数，是该实际深度的时间深度(单程旅行时)，v_i是各层介质的速度，p_x和p_y是射线参数，j是单位虚数，ω是频率，f(ω)是炮点或检波点的时域信号的傅立叶变换。

(1)式中指数函数的指数项中的累加部分可近似表达为：

$\underset{i=1}{\overset{n}{\Sigma }}\left({\mathrm{\Delta T}}_i\sqrt{1-v_i^2(p_x^2+p_y^2)}\right)\approx \sqrt{1-V_{\mathrm{rms}}^2(p_x^2+p_y^2)}·\left(\underset{i=1}{\overset{n}{\Sigma }}{\mathrm{\Delta T}}_i\right)---\left(2\right)$

对(2)式两边做Taylor展开，并截取到第二项，可得：

$V_{\mathrm{rms}}^2=\frac{\underset{i=1}{\overset{n}{\Sigma }}{v}_i^2{\mathrm{\Delta T}}_i}{\underset{i=1}{\overset{n}{\Sigma }}{\mathrm{\Delta T}}_i}---\left(3\right)$

V_rms即是叠前时间偏移所需要的偏移速度，即均方根速度。

将(2)式代入(1)式并做空间傅立叶反变换，得空间一频率域的波场为：

$P(x,y,\omega ,T)=\frac{\omega }{{4\pi }^2}\int \int f\left(\omega \right)\exp [-\mathrm{j\omega }(T\sqrt{1-V_{\mathrm{rms}}^2(p_x^2+p_y^2)}-p_{x}x-p_{y}y)]{\mathrm{dp}}_x{\mathrm{dp}}_y---\left(4\right)$

式(4)可用稳相点原理求得渐进解为：

$P(x,y,\omega ,T)=f\left(\omega \right)\frac{\omega }{2\pi }\exp (-j\frac{\pi }{2}){\left|Q(p_x^0,p_y^0)\right|}^{-1/2}\exp [-j\mathrm{\omega \phi }(p_x^0,p_y^0)]---\left(5\right)$

式(5)中，定义

$\phi (p_x,p_y)=T\sqrt{1-V_{\mathrm{rms}}^2(p_x^2+p_y^2)}-p_{x}x-p_{y}y$

$Q(p_x,p_y)=\left|\left(\begin{array}{cc}\frac{{\partial }^2\phi }{{\partial p}_x^2}& \frac{{\partial }^2\phi }{{\partial p}_x{\partial p}_y}\\ \frac{{\partial }^2\phi }{{\partial p}_x{\partial p}_y}& \frac{{\partial }^2\phi }{{\partial p}_y^2}\end{array}\right)\right|$

而是与的零点。求得和代入(5)式，可解得三维非均匀介质中地震波的走时和振幅值为

$t=T\sqrt{1+p_0}---\left(6\right)$

$A=\frac{1}{{\mathrm{TV}}_{\mathrm{rms}}^2(1+p_0)}---\left(7\right)$

式中

$p_0=\frac{x^2+y^2}{V_{\mathrm{rms}}^2{T}^2},$

既然已基于相移法和稳相点原理得到了较准确的走时和振幅值，就可利用波动方程深度偏移的反褶积成像条件进行保幅成像。若设震源是一时间脉冲，则对一个地震道，有成像结果：

$I(x,y,T)=\frac{A_r}{A_s}\int f\left(\omega \right)\mathrm{\omega exp}(-j\frac{\pi }{2})\exp (-\mathrm{j\omega }(t_s+t_r))\mathrm{d\omega }---\left(8\right)$

$=(A_r/A_s)F^{\prime }(t_s+t_r)$

式中F′(t)是f(ω)对应的时域函数的一阶导数，A_s和t_s是炮点(x_s，y_s，0)至成像点(x，y，T)的幅值和走时，可通过计算由(6)和(7)式求得，式中的V_rms是成像点的偏移速度；A_r和t_r是由同样方法求得的检波点至成像点的幅值和走时。实际应用中不需分别计算走时和幅值，直接计算

$p_s=\frac{{(x-x_s)}^2+{(y-y_s)}^2}{V_{\mathrm{rms}}^2{T}^2},$$p_r=\frac{{(x-x_r)}^2+{(y-y_r)}^2}{V_{\mathrm{rms}}^2{T}^2}$

进而计算总走时和权系数为

$t_s+t_r=(\sqrt{1+p_s}+\sqrt{1+p_r})T$

$\frac{A_r}{A_s}=\frac{1+p_s}{{1+p}_r}$

式(8)成像条件表明：对任一地震道和任一成像点，对地震信号求一阶导数，依据成像点处的偏移速度，计算总走时t_r+t_s和权系数A_r/A_s，在地震信号的一阶导数上拾取t_s+t_r时刻的值并乘上权系数，即得到该地震道在该成像点的偏移幅值，所有地震道偏移幅值的累加，即得到偏移成像的结果。

不同于忽略权系数的现行方法，(8)式的权系数保证了偏移成像中正确补偿了地震波传播的几何扩散效应。

(8)式中的F′(t)是离散函数，拾取其t_s+t_r时刻的值一般需应用插值技术。为减少插值的计算量并提高插值精度，本发明采用重采样插值技术：在频率域将最大频率f_max＝0.5/Δt扩大4倍，把增加的频率部分的频率分量置为零，然后用FFT做傅立叶反变换；此时时间采样变为0.25Δt，在这一采样下对4(t_s+t_r)/Δt取整，由相邻两点值做线性插值，即可得到很准确的F′(t_s+_tr)值。

(8)式中使用的时间深度T是单程旅行时，为与现行的以双程旅行时表达的成像结果一致，(8)式的成像结果可用I(x，y，2T)表达。

2.时变的三维偏移孔径

对单个地震道偏移计算而言，(8)式的计算在多大的成像空间，即(x，y，T)范围内进行，就是偏移计算的偏移孔径。偏移孔径对叠前偏移是重要的：较小的孔径可减少偏移计算量，但存在着不能对陡倾角构造正确成像的风险，过大的孔径又带来了偏移噪音和较大的计算量。为此，我们发明了利用构造在两个方向的最大倾角，自适应确定时变的三维偏移孔径的方法。

时变的三维偏移孔径是用各水平坐标点(x，y)上起始成像时间T_b来描述的，T_b大于成像的最大时间深度的区域就是不需要成像的区域。

本发明的偏移孔径是与地震道的偏移距和方位角有关的。首先由构造沿测线和垂直测线方向的两个最大倾角，θ_x和θ_y，构建如下不等式方程组：

式中是该地震道的方位角，满足(9)和(10)式的最小α和β，就是保证两个最大倾角为θ_x和θ_y的构造得以正确成像时，沿方位角方向和垂直方位角方向的两个最大成像角度。

以沿方位角方向和垂直方位角方向作为坐标轴看单个地震道的三维偏移脉冲响应，可将三维偏移脉冲响应看作一组椭圆以方位角方向为轴旋转形成的一组椭球面。考虑成像(即构造)倾角的偏移孔径，就是仅把倾角小于最大倾角的椭球面部分作为(8)式的计算区域。对不同时间深度T，切除了大于最大倾角部分的三维偏移脉冲响应，可看作以方位角方向为长轴方向、垂直方位角方向为短轴方向、中心在该地震道的中心点上的椭圆面，利用如下几何关系：长轴是偏移脉冲响应椭球面在沿方位角方向切面上椭圆的一个弦，该弦与椭圆交点处，椭圆切线的倾角为α；可求得弦的长度。同理可求得垂直方位角方向切面上圆的弦的长度，此时圆的切线的倾角为β。这样可求得计算区域椭圆面的长、短半轴分别为

$a=\frac{{\bar{V}}_{\mathrm{rms}}\left(T\right)·T}{\sin \left(2\alpha \right)}\sqrt{1+{\cos }^2\left(2\alpha \right)+{\lambda }^2{\sin }^2\left(2\alpha \right)-2\cos \left(2\alpha \right){[1+{\lambda }^2{\sin }^2\left(2\alpha \right)]}^{1/2}}---\left(11\right)$

$b={\bar{V}}_{\mathrm{rms}}\left(T\right)·T\tan \beta ---\left(12\right)$

式中是不同时间深度上的平均均方根速度，参数λ定义为是地震道的半偏移距，时间深度T是单程旅行时。

对时间深度T和T+ΔT(ΔT是时间深度方向上成像的采样间距)，分别计算对应的长、短半轴，这两个椭圆形成的条带就是起始成像时间为T的成像区域，对小于指定时间深度的时间深度采样点循环(大于该时间深度，偏移孔径将不随时间变化)，就可以得到时变的三维偏移孔径。

为与以双程旅行时表达的成像结果一致，偏移孔径中的起始成像时间用2Tb表达。

3.偏移噪音压制

由于有限采集缆长、有效偏移孔径等限制因素，偏移的累加计算将产生偏移噪音，本发明采用输入道成像方式，通过在单个地震道的偏移幅值累加前引入沿偏移孔径边缘的衰减，来实现偏移噪音压制。单个地震道的成像空间与偏移孔径对应，衰减系数可分别沿T增加方向和不同时间深度的椭圆边界的法线方向施加，具体方法如下：

1)定义三维矩阵(x，y，T)，对成像区域的每个水平坐标点，依据从偏移孔径得到的(x_i，y_i)点的起始成像时间T_b，计算系数

并将其存入对应(x_i，y_i，T)位置中，式中n是人为定义的边缘区包含的点数，一般取为20，ΔT是时间深度方向成像的采样；

2)对每个时间深度T，由(11)和(12)式计算a和b，再由最大的水平方向成像采样间距Δy，计算进而计算

其中是地震道的方位角。对f＞1的(x_i，y_i，T_i)，若对应的(x_i，y_i，T)位置中已存储了c₁，不做操作，否则计算将c₂存入对应的位置中；

3)步骤2)求得的系数是空间分布的，若将没赋值点的系数作为1.0，可对这些系数做空间高斯平滑，平滑后仅保留小于1.0的元素，由这些元素形成衰减系数。

对单个地震道，对偏移孔径内由(8)式计算得到的I(x_i，y_i，2T)，在偏移孔径边缘区，用衰减系数乘于偏移幅值I(x，y，2T)，即实现了边缘衰减，所有这样的单个地震道偏移结果相累加，即可有效地压制偏移噪音。

通过上述技术方案表明，本发明的有益效果是：1、本发明方法可应用于多条拖缆或多条测线记录的三维地震资料，通过求得保幅的权系数，能在成像过程中正确补偿地震波传播的几何扩散效应，得到保幅的共反射点道集；可依据构造沿测线和垂直测线方向的两个最大倾角，自适应确定时变的三维偏移孔径，从而能降低偏移噪音和提高偏移方法的计算效率；能在成像过程中自动压制偏移噪音，获得较高信噪比的偏移成像结果。2、本发明方法适宜于集群计算机的并行实现。本发明方法得到的偏移图像能指示地下构造形态、构造断裂部位和地层沉积样式，对确定有利生、储油构造和有利油气储层有重要应用价值。该方法生成的保幅共反射点(CRP)道集可更好地服务于叠前反演等油气和流体检测技术。

附图说明

图1是渤海区域JZ32-4地区常规噪音压制处理后的典型单炮记录，采集是每炮3条拖缆；

图2是偏移距为200m、方位角正弦值为0.14时的偏移孔径在沿测线和垂直测线切面上的起始成像时间(双程旅行时)曲线；

图3是偏移速度场在沿测线切面上的等值线图，图中数字是偏移速度值；

图4是经过剩余动校、拉伸和噪音切除后的典型共反射点道集；

图5是渤海区域JZ32-4区块地下潜山构造偏移结果沿测线方向切面的剖面图像。

具体实施方式

实施例1：三维保幅叠前时间偏移方法，针对渤海区域JZ32-4区块为例，具体为以下步骤：

(1)用多条拖缆或多条测线记录人工震源激发的反射地震信号，记录到磁带上。具体是，用海上拖缆采集方式获取人工震源激发的反射地震信号，记录到磁带上；每炮有3条拖缆，拖缆间距离100m，拖缆上检波器的间距(即道间距)3.125m，地震信号的记录时长5s，时间采样1ms，每炮共含3×1152＝3456道；共激发和记录2745炮。图1是常规噪音压制处理后的典型单炮记录。

(2)从磁带上读取地震资料，对叠前地震资料进行常规的压制噪音处理。针对几个共中心点位置，抽取共中心点道集，对这些道集作常规的动校正(NMO)速度拾取，对所得结果做横向平均，作为初始的、横向均匀的偏移速度场。将叠前地震资料按偏移距排序、分组，将不同组地震资料放置在集群计算机的不同计算节点上。具体如下，将基于NMO速度拾取得到的横向均匀的偏移速度场记为，将叠前地震资料共分了70组，每组地震资料仅包含有限偏移距范围内的地震道，第1组偏移距变化范围是150-250m，第2组是250-350m，以此类推，第70组是3500-3600m；分组要尽量使的各组地震资料的数据量接近相同。

(3)根据地震资料的偏移距、方位角的变化范围和地下构造的最大倾角以及初始的、横向均匀的偏移速度场(即)，构建时变的三维偏移孔径的偏移孔径表。具体如下，计算三维偏移孔径时令沿侧线方向构造的最大倾角θ_x＝50度，垂直侧线方向构造的最大倾角θ_y＝30度，偏移成像的水平坐标采样是Δx＝6.25m和Δy＝25m，时间深度(单程旅行时)采样ΔT＝0.5ms；偏移孔径表按偏移距从150到3600m，方位角正弦值从0到0.8计算，表中偏移距和方位角正弦值的采样间距分别是100m和0.14。拾取偏移孔径时，对地震道的偏移距和方位角正弦值按100m和0.14取整，直接到表中拾取。图2是偏移距为200m、方位角正弦值为0.14时的偏移孔径在沿测线和垂直测线切面上的起始成像时间(双程旅行时)曲线。

用一组数据(x_i，y_i，2T_i)(i＝1，k)描述三维偏移孔径，其中k是有效成像区域(成像空间在水平面的投影)中总的离散点数，(x_i，y_i)是成像区域中某离散点与地震道中心点的两个水平坐标差，而Ti是则是该点偏移计算的起始成像时间(单程旅行时)；令地震道的炮点坐标是(x_s，y_s)，检波点坐标是(x_r，y_r)，引入时变的偏移孔径，就是在该地震道的偏移计算中，对成像区域中水平坐标为(x_i+0.5(x_s+x_r)，y_i+0.5(y_s+y_r))的点(i＝1，k)，从T＝T_i开始进行偏移计算。

由构造沿测线和垂直测线方向的两个最大倾角，θ_x和θ_y，构建如下不等式方程组：

式中是地震道的方位角，满足上两式的最小角度α和β，就是保证最大倾角为θ_x和θ_y的构造得以正确成像时，沿方位角方向和垂直方位角方向的两个最大成像角度。

对不同时间深度T，成像区域可看作以方位角方向为长轴方向、垂直方位角方向为短轴方向、中心在地震道的中心点上的椭圆面，由α、β和初始的、横向均匀的偏移速度场可求得其长、短半轴分别为

$a=\frac{{\bar{V}}_{\mathrm{rms}}\left(T\right)·T}{\sin \left(2\alpha \right)}\sqrt{1+{\cos }^2\left(2\alpha \right)+{\lambda }^2{\sin }^2\left(2\alpha \right)-2\cos \left(2\alpha \right){[1+{\lambda }^2{\sin }^2\left(2\alpha \right)]}^{1/2}}$

$b={\bar{V}}_{\mathrm{rms}}\left(T\right)·T\tan \beta$

式中参数λ定义为是地震道的半偏移距，时间深度T是单程旅行时。

对时间深度T和T+ΔT(ΔT是时间深度方向上成像的采样间距)，分别计算对应的长、短半轴，令这两组长、短半轴形成的两个椭圆有共同的中心点(0，0)，这两个椭圆形成的条带就是起始成像时间为T的成像区域了；将落在这个条带中的离散点坐标x＝mΔx，y＝nΔy和T构成一个数据对，其中Δx和Δy是偏移成像的水平坐标采样。对小于指定时间深度的时间深度采样点循环(大于该时间深度，偏移孔径将不随时间变化)，对得到的全部(x，y，2T)分别按x和y的大小排序，统计总点数k，就可以得到时变的三维偏移孔径。

地震道的方位角和偏移距决定了偏移孔径，需根据全部地震道的方位角和偏移距变化范围，依据角度和偏移距采样间距(为减少存储量，可取较大的间距)，计算随方位角和偏移距变化的三维偏移孔径，构建偏移孔径表；偏移孔径表是一个三维矩阵，一个维数对应偏移距，另一维数对应方位角，第三维是三维偏移孔径的三个坐标数据(x_i，y_i，2T_i)和n_i共4个数随i增加排列，这里n_i是对应衰减系数的一组整数；这一维的第一个元素是该偏移孔径的总点数k。实际上，该矩阵中存放的是坐标(x_i，y_i，2T_i)对应的离散点序号，因此偏移孔径表是一整数矩阵。

方位角的采样间距不是针对实际角度值，而是针对方位角的正弦值。拾取偏移孔径时，由地震道的偏移距和炮检点在垂直侧线方向的坐标差计算方位角的正弦值，对偏移距和方位角正弦值按采样间距取整，直接到表中拾取。

(4)根据偏移孔径表确定对应的压制偏移噪音的衰减系数表。具体是，取边缘区包含的点数n＝40，时间深度(单程旅行时)采样ΔT＝0.5ms，最大的水平方向成像采样间距Δy＝25m；衰减系数表按偏移距从150到3600m，方位角正弦值从0到0.8计算，表中偏移距和方位角正弦值的采样间距分别是100m和0.14。拾取衰减系数时，对地震道的偏移距和方位角正弦值按100m和0.14取整，直接到表中拾取。

每一偏移孔径对应的衰减系数数组由如下步骤得到：

1)定义三维矩阵(x，y，T)，依据从偏移孔径得到的(x_i，y_i)点的起始成像时间T_b，计算系数并将其存入对应(x_i，y_i，T)位置中，式中n是人为定义的边缘区包含的点数，ΔT是时间深度方向成像的采样间距(单程旅行时)，对偏移孔径中全部k个点循环；

2)对每个时间深度T，由下式计算a和b：

$a=\frac{{\bar{V}}_{\mathrm{rms}}\left(T\right)·T}{\sin \left(2\alpha \right)}\sqrt{1+{\cos }^2\left(2\alpha \right)+{\lambda }^2{\sin }^2\left(2\alpha \right)-2\cos \left(2\alpha \right){[1+{\lambda }^2{\sin }^2\left(2\alpha \right)]}^{1/2}}$

$b={\bar{V}}_{\mathrm{rms}}\left(T\right)·T\tan \beta$

式中各参数与计算偏移孔径时的相同。由最大的水平方向成像采样间距Δy，计算进而计算其中是地震道的方位角。对f＞1的(x_i，y_i，T)，若对应的(x_i，y_i，T)位置中已存储了c₁，不做操作，否则计算将c₂存入对应的位置中；

3)对三维矩阵(x，y，T)的各个元素循环，将在前两个步骤中没赋值的点赋值1.0；对矩阵元素做空间高斯平滑，平滑后仅保留小于1.0的元素，由这些元素形成衰减系数数组。偏移孔径数组的每组(x_i，y_i，2T_i)在衰减系数数组中对应一个序列c_j(j＝1，n_i)，其中n_i是矩阵(x，y，T)中(x_i，y_i)处、沿T_b起小于1.0的元素个数，c_j是这些小于1.0的元素值，如有关偏移孔径表中讨论的那样，n_i将存放在偏移孔径表中。

衰减系数表是与偏移孔径表对应，偏移孔径表的每个偏移孔径在衰减系数表中对应一个衰减系数数组。衰减系数表也是一个三维矩阵，一个维数对应偏移距，另一维数对应方位角，第三维是序列c_j(j＝1，n_i)随i增加排列，其中i＝1，k，k是偏移孔径表中相同偏移距和方位角处的第一个元素值。拾取衰减系数时，对地震道的偏移距和方位角(其正弦值)按采样间距取整，直接到表中拾取。地震道的方位角和偏移距决定了偏移孔径表，也决定了衰减系数表，两者采用同样方法存放和拾取。

(5)采用基于单程波方程和深度偏移的反褶积成像条件得到的权系数进行偏移幅值计算。

定义

$p_s=\frac{{(x-x_s)}^2+{(y-y_s)}^2}{V_{\mathrm{rms}}^2{T}^2},$$p_r=\frac{{(x-x_r)}^2+{(y-y_r)}^2}{V_{\mathrm{rms}}^2{T}^2}$

式中(x_s，y_s)是地震道的炮点坐标，(x_r，y_r)是检波点坐标，(x，y，T)是成像点的水平和时间深度(单程旅行时)坐标，V_rms是成像点处的偏移速度，则炮点至成像点再到检波点的地震波走时权系数为

对地震道应用FFT(快速傅立叶变换)，在有效频带内对每个频率分量乘-jω，其中j是单位虚数，ω是频率；令地震道记录地震信号的引间采样为Δt，在频率域将频率上限取为2/Δt，将超过实际频率上限部分的值置为零；应用FFT做傅立叶反变换，此时得到的地震信号的时间采样变为0.25Δt；对4τ/Δt-0.001取整得到整数i，由变换后的时域数据上i和i+1两点的值做线性插值，得到时间τ处的值a，则偏移幅值取为

(6)在每个计算节点，对全部地震道循环。对每一地震道，由偏移孔径表确定有效的成像区域和区域中各离散点偏移计算的起始成像时间，从起始成像时间开始计算各成像点的偏移幅值。对处于三维偏移孔径的边缘区中的成像点，查衰减系数表拾取衰减系数并乘以偏移幅值得到最终的幅值。将幅值累加到存放偏移结果数组中对应的偏移距上。具体是，定义半偏移距间距Δh＝12.5m，计算m＝(h_max-h_min)/Δh+1，其中h_max为该节点上所有地震道的最大半偏移距，h_min为最小半偏移距，若三维成像空间大小为(n_x，n_y，n_T)，则该节点存放偏移结果的数组为(n_x，n_y，n_T，m)；偏移结果将依据每一地震道的半偏移距h，累加到该数组第4维的m₀＝(h-h_min)/Δh+1上。不同组地震资料可在集群计算机的不同计算节点上分别独立计算，而对应的偏移孔径表和衰减系数表只需存放该组地震资料对应的偏移距和方位角范围的部分。

(7)收集各计算节点的偏移计算结果，对成像区域的每一水平位置把偏移计算结果按时间深度、偏移距排序，形成各水平位置的共反射点道集。具体是，将各计算节点上的存放偏移结果的矩阵(n_x，n_y，n_T，m_i)整合为一个矩阵(n_x，n_y，n_T，∑m_i)，其中第4维的长度是各个计算节点该矩阵对应维长度的累加；对不同计算节点上地震资料的偏移距相同的部分，将把这些结果相叠加，然后存放到矩阵对应的位置上，矩阵第4维的长度也要相应减少。各个水平位置(n_xΔx，n_yΔy)处的二维数组(n_T，∑m_i)即是该水平位置的共反射点道集。

(8)依据共反射点道集中的同相轴是否平直，在偏移过程中确定偏移速度场。具体是，采用横向均匀的偏移速度场进行第一次偏移计算，对偏移生成的共反射点道集，依据做反动校，再做动校正得到新的速度，对这一速度做空间平滑处理，即得到最终的偏移速度场。图3是求得的偏移速度场在沿测线切面上的等值线图。

(9)利用偏移过程中确定的偏移速度场做偏移计算，对形成的全部共反射点道集做剩余动校，将出现明显拉伸和噪音部分的对应数值置为零(即切除)，将不同偏移距的数值叠加，形成偏移叠加剖面。图4是经过剩余动校、拉伸和噪音切除后的典型共反射点道集。

(10)通过显示软件将偏移叠加剖面数据转换为地下反射构造的剖面图像，剖面图像将指示地下构造的形态、断点、断距大小以及界面对地震波的反射强度。图5是用本发明方法得到的渤海区域JZ32-4区块地下潜山构造沿测线方向切面的剖面图像，潜山侧翼、顶部、断层、地层沉积样式得到很好的刻画。