波动方程偏移直接产生角道集方法

摘要

波动方程偏移直接产生角道集方法，应用于地震勘探中地震资料处理，是用于形成地震偏移的角道集的波动方程叠前深度偏移方法。该方法通过计算震源出射波在地下成像点的入射角度，直接由炮域偏移结果生成地下反射界面成像的角道集，且补偿了观测系统非均匀覆盖对幅值的改变。角道集直接对应于目的层处物性参数变化，为叠前反演提供了准确的角度域地震反射强度。这一方法的核心是直接用单程波算法确定震源波场在成像点的入射角度，避免了射线法的不稳定也与波动方程偏移方法更匹配。入射角度计算仅利用单程波的简谐波场，与偏移计算相比计算量可以忽略。该方法能应用于二维和三维波动方程叠前深度偏移处理，对油气、矿产资源勘探有重要应用价值。

说明书

技术领域

本发明属于地震勘探中反射地震资料处理技术领域，涉及地震资料处理过程中的叠前偏移成像和叠前反演技术范畴，是一种针对波动方程叠前深度偏移生成叠前偏移的角道集的方法。

背景技术

地震资料处理流程中，叠前偏移成像是最关键的环节，它的主要任务是将地表记录的反射地震信号“聚焦”到地下发生反射的真实位置。对偏移结果的一个较普遍的应用是将不同炮或不同偏移距的偏移结果叠加，得到叠加剖面，从而获得地下构造的图像；而更深入的应用是由不同炮或不同偏移距的偏移结果得到地下反射点对不同角度入射波的反射强度，据此直接识别地下介质的含油气情况。

现行的叠前时间偏移和积分法叠前深度偏移可容易得到分偏移距的叠前偏移结果，从而简单得到地下反射点对不同偏移距入射波场的反射强度，即得到随偏移距变化的共反射点(CRP)道集。速度变化简单时偏移距与入射角有较简单的转换关系，但复杂速度构造下偏移距与入射角的对应关系变得复杂，由偏移距变换到入射角已变得很困难。

当地下构造的上覆地层速度变化较剧烈时，基于金字塔面(pyramid)旅行时假设的叠前时间偏移和基于渐进(射线)理论的积分法叠前深度偏移就不能正确反应实际地震波的传播过程，因此不能准确地将反射地震信号“聚焦”到地下发生反射的真实位置。这样，即使得到随偏移距变化的共反射点(CRP)道集，由于不是来之相同的反射点，实际上是没有意义的。

波动方程叠前深度偏移是用较准确的波动理论来描述地震波的传播过程，能正确模拟复杂上覆地层导致的多次到达等复杂的波传播现象和地震波传播的弥散现象，因此即可对复杂构造正确成像，又有望得到正确反应地下构造反射强度的成像幅值。但服务于复杂构造成像的波动方程叠前深度偏移多是通过分炮的炮域偏移实现的；由炮域偏移结果得到叠加剖面是简单的，但分炮的炮域偏移结果不能如积分法那样简单地得到随偏移距变化的共反射点(CRP)道集。既然，真正反应地下构造物性参数变化的是随入射角变化的反射强度(角道集)，人们就希望能由炮域偏移结果直接得到角道集。

已发展了两类由炮域波动方程深度偏移方法生成角道集的方法。一是在应用成像条件前对炮点的入射波场和检波点的反传波场进行局部平面波分解，对每对平面波分别成像，二是在炮域偏移的成像条件中引入“局部偏移距”参量，形成多维成像体，通过倾斜叠加等角度变换方法得到角道集。第一种方法需对逐个成像点进行局部平面波分解，将单点的成像转换成一个成像矩阵，三维偏移计算时计算量难以承受；第二种方法尽管计算量较第一种有所减少，但引入“局部偏移距”参量仍导致成像计算的内存需求聚增。两种方法均不能考虑炮点覆盖不均匀导致的成像幅值误差，且偏移计算的空间采样控制了这两类方法所生成的角道集的有效频率上限。

一个自然的想法是直接计算炮点的出射波对成像点的入射角度，然后由炮域偏移结果直接通过部分叠加得到角道集。但现行的计算入射角的方法是基于射线追踪算法，这一高频近似方法对速度模型较敏感，复杂速度变化情况下难以得到正确的入射角度，因此这一方法在实际应用中遇到许多困难。

由波动方程叠前深度偏移生成叠前偏移的角道集仍是一个需要深入研究的问题，只有准确获得角道集，才能较全面地发挥波动方程叠前深度偏移的优势，从而为复杂上覆地层下直接识别储层的含流体和含油气奠定基础，降低复杂构造情况下油气储层的勘探风险。

发明内容

本发明的目的是：提供一种由炮域波动方程叠前深度偏移产出叠前偏移的角道集的方法，能在复杂上覆地层情况下得到直接反应储层物性参数变化的角道集，且能补偿观测系统非均匀覆盖对幅值的改变；产出的角道集为叠前反演技术提供了准确的随入射角变化的地震反射强度。与常规的炮域波动方程叠前深度偏移方法相比，这一方法在形成角道集时不增加太多计算量；它克服了现行的各类波动方程偏移生成角道集技术存在的计算量和内存需求巨大的双重问题。这一方法将三维地震资料的波动方程偏移生成角道集以及基于角道集的反演推进到工业化应用阶段。

本发明采用的技术方案是：波动方程深度偏移直接产生角道集方法，具体步骤包括：

(1)用单边或双边采集方式，由检波器获取人工震源激发的反射地震信号，记录到磁带上；

(2)从磁带上读取反射地震信号，按叠前炮集排序存放到计算机中；

(3)利用波动方程偏移的单程波算法分别将炮点的震源波场和检波点的接收波场沿深度轴进行波场深度延拓，对延拓后的炮点和检波点波场应用反褶积成像条件进行成像，具体是：在频率空间域，用延拓后的炮点波场除检波点波场，并将不同频率的复波场值的实部加权累加，除法计算采用可保证计算结果稳定的算法；

(4)利用单程波方法计算炮点的出射波对地下各成像点的入射角度，具体是：对炮点的主频分量进行深度延拓，用两种算法对延拓结果进行计算得到两个空间场，计算两个场的商并进行中值滤波、平滑和归一化，由反余弦求得炮点出射波对地下各成像点的入射角度；

(5)依据得到的各炮点出射波的入射角度，将由反褶积成像条件得到的各个单炮偏移结果按角度累加并补偿非均匀覆盖对幅值的影响，形成角道集；

(6)对角道集做剩余动校，将出现明显拉伸部分的对应数值置为零，用奇异值分解方法压制随机噪音，即得到反应介质界面物理参数变化的叠前偏移的角道集。

(7)由角道集中读取不同空间位置处成像幅值随角度的变化，可由反演方法获得地层的岩石物理参数，岩石物理参数图像即能指示地下构造的含油气情况。

所述的对延拓后的炮点和检波点波场应用反褶积成像条件进行成像是这样实现的：令T₀是地震记录的时长，则频率采样Δω＝2π/T₀(弧度/秒)；设Δz是深度方向的空间采样，n是深度上采样点的顺序数，x和y是水平位置坐标，m是频率采样点的顺序数；若炮点和检波点波场经深度延拓后在nΔz深度上分别为P_D(x，y，nΔz，mΔω)和P_U(x，y，nΔz，mΔω)，计算

A(x，y，nΔz，mΔω)＝P_D(x，y，nΔz，mΔω){P_D(x，y，nΔz，mΔω)}^*

式中上标*是复共额，炮点和检波点波场P_D和P_U的单位依据检波器类型可以是速度或压力，频率的单位是弧度/秒，深度和水平坐标的单位为米。令A₀(nΔz)是A(x，y，nΔz，mΔω)在深度nΔz上的平均值，当A(x，y，nΔz，mΔω)≥1.2A₀(nΔz)，计算

$B(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega })=\frac{P_U(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega }){\left\{P_D(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega })\right\}}^{*}}{A(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega })}$

当A(x，y，nΔz，mΔω)＜1.2A₀(nΔz)，计算

$B(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega })=\frac{P_U(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega }){\left\{P_D(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega })\right\}}^{*}}{A_0\left(\mathrm{n\Delta z}\right)}(2-a)(1+{(1-a)}^2),$其中$a=\frac{A(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega })}{A_0\left(\mathrm{n\Delta z}\right)}$

令地震资料有效频带的下、上限分别是l₁Δω和l₂Δω，这里l₁和l₂是正整数，则空间各位置处反褶积成像条件的成像幅值φ(x，y，nΔz)为

$\phi (x,y,\mathrm{n\Delta z})=\underset{m=l_1}{\overset{l_2}{\Sigma }}\mathrm{Re}\left\{\mathrm{\beta B}(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega })\right\}$

若j是单位虚数，上式中做三维偏移计算时β＝jmΔω，二维偏移计算时$\beta =\sqrt{\mathrm{jm\Delta \omega }};$符号Re代表求复波场值的实部。

所述的利用单程波方法计算炮点的出射波对地下各成像点的入射角度是这样实现的：令x和y是水平坐标，z是深度坐标，震源的主频是ω_p(弧度/秒)，震源的空间位置为(x_s，y_s，0)，介质的速度为c(x，y，z)；利用波动方程偏移的单程波算法，在频率ω_p对空间脉冲δ(x-x_s，y-y_s)进行深度延拓，得到以Δz为间距的各深度离散点上的复波场值P(x，y，nΔz，ω_p，x_s，y_s)；分别计算

I₀(x，y，nΔz)＝|Re{P(x，y，nΔz，ω_p，x_s，y_s)}|

$I_1(x,y,\mathrm{n\Delta z})=\mathrm{Re}\left\{j\frac{c(x,y,\mathrm{n\Delta z})}{{\omega }_p}P(x,y,\mathrm{n\Delta z},{\omega }_p,x_s,y_s)\right\}$

$I_2(x,y,\mathrm{n\Delta z})=\left|\frac{2}{3\mathrm{\Delta z}}\right[I_1(x,y,\mathrm{n\Delta z}+\mathrm{\Delta z})-I_1(x,y,\mathrm{n\Delta z}-\mathrm{\Delta z})]$

$-\frac{1}{12\mathrm{\Delta z}}[I_1(x,y,\mathrm{n\Delta z}+2\mathrm{\Delta z})-I_1(x,y,\mathrm{n\Delta z}-2\mathrm{\Delta z})]|$

上式中x_s，y_s和ω_p是参数，故在左端的函数中将其忽略；符号||代表求绝对值，n是深度上采样点的顺序数；速度c(x，y，z)的单位是米/秒，频率的单位是弧度/秒，深度和水平坐标的单位是米，空间脉冲的单位是压强。若I₂(x，y，nΔz)＞1.05I₀(x，y，nΔz)，则令I₂(x，y，nΔz)＝0。令φ₀(nΔz)是I₀(x，y，nΔz)在深度nΔz上的平均值，当I₀(x，y，nΔz)≥1.2φ₀(nΔz)，计算

$\phi (x,y,\mathrm{n\Delta z})=\frac{I_2(x,y,\mathrm{n\Delta z})}{I_0(x,y,\mathrm{n\Delta z})}$

当I₀(x，y，nΔz)＜1.2φ₀(nΔz)，计算

$\phi (x,y,\mathrm{n\Delta z})=\frac{I_2(x,y,\mathrm{n\Delta z})}{{\phi }_0\left(\mathrm{n\Delta z}\right)}(2-a)(1+{(1-a)}^2),$其中$a=\frac{I_0(x,y,\mathrm{n\Delta z})}{{\phi }_0\left(\mathrm{n\Delta z}\right)}$对空间各点的φ(x，y，nΔz)值进行如下处理：

(1)沿深度方向进行中值滤波；

(2)进行空间三维平滑；

(3)对每一深度nΔz，选取最大值，用这一值除该深度的每一点的数值，实现该深度上的最大值归一化。

令上述处理后的值为φ(x，y，nΔz)，则可得震源对地下各成像点的入射角度θ(单位为度)为：

θ＝acos(φ(x，y，nΔz))

式中(x，y，nΔz)即是成像点的空间坐标。。

所述的依据得到的各炮点的出射波的入射角度，将由反褶积成像条件得到的各个单炮偏移结果按角度累加并补偿非均匀覆盖对幅值的影响，形成角道集是这样实现的：定义角度采样间距Δθ(度)，对每一成像点定义一个一维数组，数组长度小于90/θΔ；对全部炮循环，对每炮的单炮偏移结果，读取由步骤(4)得到的该炮点的出射波对各成像点的入射角度θ，计算l＝1+int{θ/Δθ}，将该成像点的幅值累加到对应数组的第l个元素上，同时统计累加到该元素上的炮数；完成对全部炮循环，将同一水平位置的不同深度点的成像幅值数组按深度、角度形成二位数组，即得到未经非均匀覆盖补偿的角道集；用统计得到的累加炮数去除各成像点的不同角度的数值，就得到完成非均匀覆盖补偿的角道集。

本发明的波动方程深度偏移直接产生角道集方法，能直接由分炮的炮域偏移结果，经局部叠加得到地下反射点的叠前偏移的角道集。

本发明的波动方程深度偏移直接产生角道集方法，能剔除复杂上覆地层的影响，得到直接反应目的层处物性参数变化的叠前偏移角道集。

本发明的波动方程深度偏移直接产生角道集方法，能补偿观测系统非均匀覆盖对幅值的改变。

本发明的波动方程深度偏移直接产生角道集方法，与常规的炮域波动方程叠前深度偏移方法相比，不增加太多计算量。

本发明的波动方程深度偏移直接产生角道集方法，可应用于二维和三维波动方程叠前深度偏移处理。

本发明的波动方程深度偏移直接产生角道集方法，能用于反演地层的岩石物理参数，直接得到地下构造的含油气情况。

本发明的具体实现原理如下：

本发明的核心思想是直接用单程波方法计算炮点的出射波在地下成像点的入射角度，避免了射线法对速度模型的不稳定也与波动方程偏移方法更匹配。入射角度计算仅利用单程波的简谐波场，与偏移算法相比计算量可以忽略。

计算入射角度的原理是，由单程波方程可知，下行的频率域波场满足

$\frac{\partial {\tilde{P}}^{+}(x,y,z;\omega )}{\partial z}=-j{k}_z{\tilde{P}}^{+}(x,y,z;\omega )---\left(1\right)$

其中j是单位虚数，k_z是垂直波数，ω是角频率，x和y是水平位置坐标，z是深度坐标，由频散关系得

$k_z=\frac{\omega }{c(x,y,z)}\cos \theta ---\left(2\right)$

式中c(x，y，z)是空间各位置处的介质波速，角度θ即是空间位置(x，y，z)处波前面法线与深度轴z的夹角。将(2)式代入(1)式并做离散的时域付氏反变换，且令时间t＝0，可得

$\frac{\partial \left[\underset{\omega >0}{\Sigma }\mathrm{Re}\left(j\frac{c(x,y,z)}{\omega }{\tilde{P}}^{+}(x,y,z,\omega )\right)\right]}{\partial z}=\cos \theta \left[\underset{\omega >0}{\Sigma }\mathrm{Re}\left({\tilde{P}}^{+}(x,y,z,\omega )\right)\right]---\left(3\right)$式中符号Re  代表求复波场值的实部。若令${\tilde{P}}^{+}(x,y,z=0,\omega )=\delta (x-x_s,y-y_s)f\left(\omega \right),$其中f(ω)是地震子波的谱，则(3)式中对应于震源点置于(x_s，y_s，0)处时的脉冲响应，而波前面法线就是(x_s，y_s，0)处炮点的出射波在(x，y，z)点的入射波方向。

脉冲响应在空间中仅存在一个波前面，对没有波前面的空间位置，不能用波前面来确定炮点入射波的方向。若仅对主频分量计算脉冲响应，即仅计算Re则波前面将分布于整个区域，易于求得所有成像点的入射波方向，而此时计算脉冲响应的计算量也大幅减少。仅考虑主频分量，由(3)式可求得

$\cos \theta =\frac{\partial \left[\mathrm{Re}\left(j\frac{c(x,y,z)}{{\omega }_p}{\tilde{P}}^{+}(x,y,z,{\omega }_p)\right)\right]/\partial z}{\mathrm{Re}\left({\tilde{P}}^{+}(x,y,z,{\omega }_p)\right)}---\left(4\right)$

式(4)中的分母将存在多个零点，此零点处将不能得到正确的入射波角度；既然在一个波长范围内入射波的角度不会有显著变换，可用波前面的波锋和波谷处的入射波角度插值得到邻近位置的入射波角度。我们了设计如下平滑算法，实现这一思想。首先计算

$I_0(x,y,z)=\left|\mathrm{Re}\left({\tilde{P}}^{+}(x,y,z,{\omega }_p)\right)\right|---\left(5\right)$

$I_1(x,y,z)=\mathrm{Re}\left(j\frac{c(x,y,z)}{{\omega }_p}{\tilde{P}}^{+}(x,y,z,{\omega }_p)\right)---\left(6\right)$

用四阶差分算法计算(4)式中的一阶导数，令计算脉冲响应的深度方向的空间采样为Δz，可得

$I_2(x,y,z)=\left|\frac{2}{3\mathrm{\Delta z}}\right[I_1(x,y,z+\mathrm{\Delta z})-I_1(x,y,z-\mathrm{\Delta z})]---\left(7\right)$

$-\frac{1}{12\mathrm{\Delta z}}[I_1(x,y,z+2\mathrm{\Delta z})-I_1(x,y,z-2\mathrm{\Delta z})]$

式(5)和(7)中采用绝对值，是因为对下行波而言cosθ总是正。由(4)式可知，I₂(x，y，z)应该大于I₀(x，y，z)。但当偏移噪音较强时(即在入射波较弱的阴影区)，有时并不满足这一条件；这表明该处不存在准确的入射波，因此不能得到正确的入射角度。在计算中，若I₂(x，y，z)＞1.05I₀(x，y，z)，可令I₂(x，y，z)＝0，这样避免了阴影区角度计算的困难。

为避免式(4)中除法计算的不稳定，我们发展了高精度的稳定的相除算法，这一算法的核心是利用级数展开近似式：

$\frac{1}{x}=\frac{1}{1-(1-x)}\approx 1+(1-x)+{(1-x)}^2+{(1-x)}^3=(2-x)(1+{(1-x)}^2)---\left(8\right)$

当x趋于零，式(8)仍可得到稳定的结果，而当|x|≤1.2时，式(8)的近似有很好的精度。利用式(8)，我们发展了稳定的相除算法，

$\phi (x,y,z)=\frac{I_2(x,y,z)}{I_0(x,y,z)}=\left(\begin{array}{c}\left(\begin{array}{cc}\frac{I_2(x,y,z)}{I_0(x,y,z)},& I_0(x,y,z)\ge 1.2{\phi }_0\left(z\right)\end{array}\right)\\ \frac{I_2(x,y,z)}{{\phi }_0\left(z\right)}(2-\frac{I_0(x,y,z)}{{\phi }_0\left(z\right)})(1+{(1-\frac{I_0(x,y,z)}{{\phi }_0\left(z\right)})}^2)\end{array}\right)---\left(9\right)$

式中φ₀(z)是I₀(x，y，z)在深度z上的平均值。

对(9)式的结果，先做沿深度方向的中值滤波，然后进行空间三维平滑，即可由波锋和波谷处的入射波角度插值得到邻近位置的入射角度；为避免(7)式差分等计算产生的误差，将对每一深度上平滑处理后的φ(x，y，z)进行归一化处理，保证垂直入射的存在。这样就可通过单个频率的简谐波场的深度延拓，求得炮点出射波对地下各点的入射角度

θ＝acos(φ(x，y，z))(10)

式中(x，y，z)即是地下各点的空间坐标。

本发明对炮域偏移采用反褶积成像条件，而不是采用相关成像条件；反褶积成像条件可补偿地震波的扩散效应，获得正确的幅值。

我们再次利用式(8)来解决反褶积成像条件中除法的稳定问题。炮域偏移的反褶积成像条件可表达为

$\phi (x,y,z)=\underset{m=l_1}{\overset{l_2}{\Sigma }}\mathrm{Re}\left(\beta \frac{P_U(x,y,z,\mathrm{m\Delta \omega })}{P_D(x,y,z,\mathrm{m\Delta \omega })}\right)---\left(11\right)$

式中符号Re代表求复波场值的实部，Δω是频率采样，m是频率采样点的顺序数，P_U(x，y，z，mΔω)和P_D(x，y，z，mΔω)分别是深度延拓后深度z上频率域的炮点和检波点波场；l₁和l₂是正整数，l₁Δω和l₂Δω分别对应地震资料有效频带的下、上限；引入参数β是为了保持成像的波形与检波点波场相一致，若j是单位虚数，做三维偏移计算时β＝jmΔω，二维偏移计算时$\beta =\sqrt{\mathrm{jm\Delta \omega }}.$

式(11)中的除法可进一步表达为

$\frac{P_U(x,y,z,\mathrm{m\Delta \omega })}{P_D(x,y,z,\mathrm{m\Delta \omega })}=\frac{P_U(x,y,z,\mathrm{m\Delta \omega }){\left\{P_D(x,y,z,\mathrm{m\Delta \omega })\right\}}^{*}}{P_D(x,y,z,\mathrm{m\Delta \omega }){\left\{P_D(x,y,z,\mathrm{m\Delta \omega })\right\}}^{*}}---\left(12\right)$

式中上标*代表复共额。令式(12)中实数的分母为A(x，y，z，mΔω)，由(8)式可得

$\frac{1}{A(x,y,z,\mathrm{m\omega })}=\left(\begin{array}{c}\left(\begin{array}{cc}\frac{1}{A(x,y,z,\mathrm{m\omega })},& A(x,y,z,\mathrm{m\omega })\ge 1.2A_0\left(z\right)\end{array}\right)\\ \frac{1}{A_0\left(z\right)}(2-\frac{A(x,y,z,\mathrm{m\omega })}{A_0\left(z\right)})(1+{(1-\frac{A(x,y,z,\mathrm{m\omega })}{A_0\left(z\right)})}^2)\end{array}\right)---\left(13\right)$

式中A₀(z)是A(x，y，z，mω)在深度z上的平均值。将式(13)代入(12)再代入(11)，可得稳定的反褶积成像条件。

反褶积成像条件得到的偏移幅值，反应了成像点处炮点波场的入射角度所对应的正确反射强度。在形成角道集时，尽管角度的采样是相同的，但不同角度上累加的炮数是不同的，角度越大，累加的数量越多；而炮点覆盖越密，累加的数量也就越多。为获得正确地反应成像点物理参数对比的成像幅值，我们在累加过程中统计各角度间距内累加的炮数，用这一数量去除累加结果。这样即考虑了成像结果在角度域分布的不均匀，也补偿了炮点覆盖不均匀的影响。

本发明的有益效果：该方法可应用于复杂上覆地层下油气储层的地震偏移成像，能生成直接反应目的层物理参数变化的叠前偏移的角道集。该角道集能更好地服务于叠前反演等油气和流体直接检测技术(避免了各种基于测井资料的补偿)，降低了无井或少井地区油气勘探的风险。该方法将波动方程叠前深度偏移由构造成像推进到岩性识别，拓展和提升了波动方程叠前深度偏移的应用范围和效果。该方法对存在复杂上覆层的大量区域的油气、矿产资源勘探有重要应用价值。

附图说明

图1是大庆LMD地区开发地震的典型单炮记录，采集是16线，双边观测，道间距40m；本次采用波动方程偏移直接生成角道集方法共使用了3600炮，满覆盖11.2平方公里。

图2a是叠前深度偏移使用的三维速度模型的沿测线322的二维切片，图2b是过CDP400的沿垂直测线方向的二维速度切片。

图3是线号322和CDP号400处由所有单炮偏移结果构成的共成像点道集，单炮偏移采用了本发明的反褶积成像条件。

图4a和图4b是炮点在线号318和CDP号415处的炮点波场的入射角分布图，入射角计算采用了本发明的新方法；图4a是沿测线方向的二维角度切片，图4b是沿垂直测线方向的二维角度切片。

图5是线号318和CDP号415处的未进行非均匀覆盖补偿的角道集，这一道集是依据各炮点在线号318和CDP号415处的入射角度，由各单炮偏移结果按角度做局部叠加得到。

图6是对图5角道集应用非均匀补偿后的结果，对比图5和图6可知，小角度的幅值得到增强。

图7是对图6角道集首先应用剩余动校，再将出现明显拉伸部分的对应数值置为零，然后用奇异值分解方法压制随机噪音得到的最终角道集。

图8a是图7角道集的局部放大图，图8b是在两个典型层位处的随角度变化的幅值对比。图8a的局部放大图对应了有测井资料的地层部分。图8b的对比图中的实线是由图7读出的深度970m和1030m处幅值随角度的变化，虚线是依据该位置(线号318和CDP号415)处声波和密度测井资料读取的P波、S波速度以及密度，由理论计算得到的随角度变化的反射强度。从图8b中对比可见两者拟合很好，这表面本发明求得的角道集不需“校正”即可较好地拟合反射界面的入射角与反射幅值关系(AVA)特征。

图9a是线号323和CDP号360的井位处的角道集，图9b是角道集的局部放大，图9c是两个典型层位(深度840m和975m)处幅值随角度变化与理论AVA特征的对比。各图的意义与图7、图8a和图8b相同。

图10a是线号363和CDP号448的井位处的角道集，图10b是角道集的局部放大，图10c是两个典型层位(深度995m和1095m)处幅值随角度变化与理论AVA特征的对比。各图的意义与图7、图8a和图8b相同。

图11a和图11b是本次偏移处理得到的大庆LMD区块偏移叠加剖面。图11a是沿测线322的二维切片，图11b是过CDP400的沿垂直测线方向的二维切片。图11a和图11b是由各成像点的角道集沿角度方向叠加得到的。

具体实施方式

实施例1：波动方程深度偏移直接产生角道集方法，针对大庆LMD地区为例，具体为以下步骤：

(1)用双边采集方式，由检波器获取人工震源激发的反射地震信号，记录到磁带上，每炮布设16条采集线，采集线上的道间距为40m，每炮共含16×168＝2688道；共激发和记录3600炮，图1是典型的单炮记录。

(2)从磁带上读取反射地震信号，按叠前炮集排序存放到PC-Cluster计算机中，偏移计算将通过读取每个单炮记录，对炮循环来进行。

(3)依据图2a和2b的偏移速度模型，利用波动方程偏移的单程波算法分别将炮点的震源波场(即空间和时间脉冲函数)和炮记录(接收到的地震信号)沿深度轴进行波场延拓，对延拓后的炮点和检波点波场应用反褶积成像条件进行成像。

令T₀是地震记录的时长，则频率采样Δω＝2π/T₀(弧度/秒)；设Δz是深度方向的空间采样，n是深度上采样点的顺序数，x和y是水平位置坐标，m是频率采样点的顺序数；若炮点和检波点波场经深度延拓后在nΔz深度上分别为P_D(x，y，nΔz，mΔω)和P_U(x，y，nΔz，mΔω)，计算

A(x，y，nΔz，mΔω)＝P_D(x，y，nΔz，mΔω){P_D(x，y，nΔz，mΔω)}^*式中上标*是复共额，炮点和检波点波场P_D和P_U的单位依据检波器类型可以是速度或压力，频率的单位是弧度/秒，深度和水平坐标的单位为米。令A₀(nΔz)是A(x，y，nΔz，mΔω)在深度nΔz上的平均值，当A(x，y，nΔz，mΔω)≥1.2A₀(nΔz)，计算

$B(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega })=\frac{P_U(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega }){\left\{P_D(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega })\right\}}^{*}}{A(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega })}$

当A(x，y，nΔz，mΔω)＜1.2A₀(nΔz)，计算

$B(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega })=\frac{P_U(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega }){\left\{P_D(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega })\right\}}^{*}}{A_0\left(\mathrm{n\Delta z}\right)}(2-a)(1+{(1-a)}^2),$其中

$a=\frac{A(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega })}{A_0\left(\mathrm{n\Delta z}\right)}$

令地震资料有效频带的下、上限分别是l₁Δω和l₂Δω，这里l₁和l₂是正整数，则空间各位置处反褶积成像条件的成像幅值φ(x，y，nΔz)为

$\phi (x,y,\mathrm{n\Delta z})=\underset{m=l_1}{\overset{l_2}{\Sigma }}\mathrm{Re}\left\{\mathrm{\beta B}(x,y,\mathrm{n\Delta z},\mathrm{m\Delta \omega })\right\}$

若j是单位虚数，上式中做三维偏移计算时β＝jmΔω，二维偏移计算时$\beta =\sqrt{\mathrm{jm\Delta \omega }};$符号Re代表求复波场值的实部。

将成像结果φ(x，y，nΔz)存储在计算机的硬盘中。图3是线号322和CDP号400处由所有单炮偏移结果构成的共成像点道集；所有单炮偏移采用的参数均是：Δz＝5m，深度上的采样点总数n＝1000，Δω＝1.40(弧度/秒)，l₁＝13，l₂＝337。

(4)对每个炮点(震源)，依据图2a和2b的偏移速度模型，利用单程波方法对频率为25Hz的空间脉冲进行深度延拓，取Δz＝5m，求得炮点的出射波对地下各成像点的入射角度。

计算过程是：令x和y是水平坐标，z是深度坐标，震源的主频是ω_p(弧度/秒)，震源的空间位置为(x_s，y_s，0)，介质的速度为c(x，y，z)；利用波动方程偏移的单程波算法，在频率ω_p对空间脉冲δ(x-x_s，y-y_s)进行深度延拓，得到以Δz为间距的各深度离散点上的复波场值P(x，y，nΔz，ω_p，x_s，y_s)；分别计算

I₀(x，y，nΔz)＝|Re{P(x，y，nΔz，ω_p，x_s，y_s)}|

$I_1(x,y,\mathrm{n\Delta z})=\mathrm{Re}\left\{j\frac{c(x,y,\mathrm{n\Delta z})}{{\omega }_p}P(x,y,\mathrm{n\Delta z},{\omega }_p,x_s,y_s)\right\}$

$I_2(x,y,\mathrm{n\Delta z})=\left|\frac{2}{3\mathrm{\Delta z}}\right[I_1(x,y,\mathrm{n\Delta z}+\mathrm{\Delta z})-I_1(x,y,\mathrm{n\Delta z}-\mathrm{\Delta z})]$

$-\frac{1}{12\mathrm{\Delta z}}[I_1(x,y,\mathrm{n\Delta z}+2\mathrm{\Delta z})-I_1(x,y,\mathrm{n\Delta z}-2\mathrm{\Delta z})|$

上式中x_s，y_s和ω_p是参数，故在左端的函数中将其忽略；符号||代表求绝对值，n是深度上采样点的顺序数；速度c(x，y，z)的单位是米/秒，频率的单位是弧度/秒，深度和水平坐标的单位是米，空间脉冲的单位是压强。若I₂(x，y，nΔz)＞1.05I₀(x，y，nΔz)，则令I₂(x，y，nΔz)＝0。令φ₀(nΔz)是I₀(x，y，nΔz)在深度nΔz上的平均值，当I₀(x，y，nΔz)≥1.2φ₀(nΔz)，计算

$\phi (x,y,\mathrm{n\Delta z})=\frac{I_2(x,y,\mathrm{n\Delta z})}{I_0(x,y,\mathrm{n\Delta z})}$

当I₀(x，y，nΔz)＜1.2φ₀(nΔz)，计算

$\phi (x,y,\mathrm{n\Delta z})=\frac{I_2(x,y,\mathrm{n\Delta z})}{{\phi }_0\left(\mathrm{n\Delta z}\right)}(2-a)(1+{(1-a)}^2),$其中$a=\frac{I_0(x,y,\mathrm{n\Delta z})}{{\phi }_0\left(\mathrm{n\Delta z}\right)}$

对空间各点的φ(x，y，nΔz)值进行如下处理：

1)沿深度方向进行中值滤波；

2)进行空间三维平滑；

3)对每一深度nΔz，选取最大值，用这一值除该深度的每一点的数值，实现该深度上的最大值归一化。

令上述处理后的值为φ(x，y，nΔz)，则可得震源对地下各成像点的入射角度θ(单位为度)为：

θ＝acos(φ(x，y，nΔz))

式中(x，y，nΔz)即是成像点的空间坐标。

图4a和图4b是计算得到的某炮点出射波的空间入射角分布图。

(5)读取各炮的单炮偏移结果，根据计算得到的该炮点出射波的入射角度值，将偏移结果按角度累加并补偿非均匀覆盖对幅值的影响，形成角道集，角度采样间距取为Δθ＝2°。

计算过程是：定义角度采样间距Δθ(度)，对每一成像点定义一个一维数组，数组长度小于90/Δθ；对全部炮循环，对每炮的单炮偏移结果，读取由步骤(4)得到的该炮点的出射波对各成像点的入射角度θ，计算l＝1+int{θ/Δθ}，将该成像点的幅值累加到对应数组的第l个元素上，同时统计累加到该元素上的炮数；完成对全部炮循环，将同一水平位置的不同深度点的成像幅值数组按深度、角度形成二位数组，即得到未经非均匀覆盖补偿的角道集；用统计得到的累加炮数去除各成像点的不同角度的数值，就得到完成非均匀覆盖补偿的角道集。

图5是某水平位置处未进行非均匀覆盖补偿的角道集。图6是某水平位置处完成非均匀覆盖补偿的角道集。

(6)对全部角道集循环，做剩余动校，将出现明显拉伸部分的对应数值置为零，用奇异值分解方法压制随机噪音，图7是完成后续处理的角道集。

(7)本领域技术人员能由角道集中读取不同空间位置处幅值随角度的变化，并由反演方法获得地层的岩石物理参数，岩石物理参数图像即能指示地下构造的含油气情况。参阅图8a和8b，图9a、9b和9c，图10a、10b和10c，和图11a和11b。