专利名称：波动方程偏移直接产生角道集方法
技术领域：
本发明属于地震勘探中反射地震资料处理技术领域，涉及地震资料处理过程中的 叠前偏移成像和叠前反演技术范畴，是一种针对波动方程叠前深度偏移生成叠前偏移的角 道集的方法。
背景技术：
地震资料处理流程中，叠前偏移成像是最关键的环节，它的主要任务是将地表记 录的反射地震信号“聚焦”到地下发生反射的真实位置。对偏移结果的一个较普遍的应用 是将不同炮或不同偏移距的偏移结果叠加，得到叠加剖面，从而获得地下构造的图像；而更 深入的应用是由不同炮或不同偏移距的偏移结果得到地下反射点对不同角度入射波的反 射强度，据此直接识别地下介质的含油气情况。现行的叠前时间偏移和积分法叠前深度偏移可容易得到分偏移距的叠前偏移结 果，从而简单得到地下反射点对不同偏移距入射波场的反射强度，即得到随偏移距变化的 共反射点(CRP)道集。速度变化简单时偏移距与入射角有较简单的转换关系，但复杂速度 构造下偏移距与入射角的对应关系变得复杂，由偏移距变换到入射角已变得很困难。当地下构造的上覆地层速度变化较剧烈时，基于金字塔面(pyramid)旅行时假设 的叠前时间偏移和基于渐进(射线)理论的积分法叠前深度偏移就不能正确反应实际地 震波的传播过程，因此不能准确地将反射地震信号“聚焦”到地下发生反射的真实位置。这 样，即使得到随偏移距变化的共反射点(CRP)道集，由于不是来之相同的反射点，实际上是 没有意义的。波动方程叠前深度偏移是用较准确的波动理论来描述地震波的传播过程，能正确 模拟复杂上覆地层导致的多次到达等复杂的波传播现象和地震波传播的弥散现象，因此即 可对复杂构造正确成像，又有望得到正确反应地下构造反射强度的成像幅值。但服务于复 杂构造成像的波动方程叠前深度偏移多是通过分炮的炮域偏移实现的；由炮域偏移结果得 到叠加剖面是简单的，但分炮的炮域偏移结果不能如积分法那样简单地得到随偏移距变化 的共反射点(CRP)道集。既然，真正反应地下构造物性参数变化的是随入射角变化的反射 强度(角道集)，人们就希望能由炮域偏移结果直接得到角道集。已发展了两类由炮域波动方程深度偏移方法生成角道集的方法。一是在应用成像 条件前对炮点的入射波场和检波点的反传波场进行局部平面波分解，对每对平面波分别成 像，二是在炮域偏移的成像条件中引入“局部偏移距”参量，形成多维成像体，通过倾斜叠加 等角度变换方法得到角道集。第一种方法需对逐个成像点进行局部平面波分解，将单点的 成像转换成一个成像矩阵，三维偏移计算时计算量难以承受；第二种方法尽管计算量较第 一种有所减少，但引入“局部偏移距”参量仍导致成像计算的内存需求聚增。两种方法均不 能考虑炮点覆盖不均勻导致的成像幅值误差，且偏移计算的空间采样控制了这两类方法所 生成的角道集的有效频率上限。—个自然的想法是直接计算炮点的出射波对成像点的入射角度，然后由炮域偏移结果直接通过部分叠加得到角道集。但现行的计算入射角的方法是基于射线追踪算法，这 一高频近似方法对速度模型较敏感，复杂速度变化情况下难以得到正确的入射角度，因此 这一方法在实际应用中遇到许多困难。由波动方程叠前深度偏移生成叠前偏移的角道集仍是一个需要深入研究的问题， 只有准确获得角道集，才能较全面地发挥波动方程叠前深度偏移的优势，从而为复杂上覆 地层下直接识别储层的含流体和含油气奠定基础，降低复杂构造情况下油气储层的勘探风险。
发明内容
本发明的目的是提供一种由炮域波动方程叠前深度偏移产出叠前偏移的角道集 的方法，能在复杂上覆地层情况下得到直接反应储层物性参数变化的角道集，且能补偿观 测系统非均勻覆盖对幅值的改变；产出的角道集为叠前反演技术提供了准确的随入射角变 化的地震反射强度。与常规的炮域波动方程叠前深度偏移方法相比，这一方法在形成角道 集时不增加太多计算量；它克服了现行的各类波动方程偏移生成角道集技术存在的计算量 和内存需求巨大的双重问题。这一方法将三维地震资料的波动方程偏移生成角道集以及基 于角道集的反演推进到工业化应用阶段。本发明采用的技术方案是波动方程深度偏移直接产生角道集方法，具体步骤包 括(1)用单边或双边采集方式，由检波器获取人工震源激发的反射地震信号，记录到 磁带上；(2)从磁带上读取反射地震信号，按叠前炮集排序存放到计算机中；(3)利用波动方程偏移的单程波算法分别将炮点的震源波场和检波点的接收波 场沿深度轴进行波场深度延拓，对延拓后的炮点和检波点波场应用反褶积成像条件进行成 像，具体是在频率空间域，用延拓后的炮点波场除检波点波场，并将不同频率的复波场值 的实部加权累加，除法计算采用可保证计算结果稳定的算法；(4)利用单程波方法计算炮点的出射波对地下各成像点的入射角度，具体是对 炮点的主频分量进行深度延拓，用两种算法对延拓结果进行计算得到两个空间场，计算两 个场的商并进行中值滤波、平滑和归一化，由反余弦求得炮点出射波对地下各成像点的入 射角度；(5)依据得到的各炮点出射波的入射角度，将由反褶积成像条件得到的各个单炮 偏移结果按角度累加并补偿非均勻覆盖对幅值的影响，形成角道集；(6)对角道集做剩余动校，将出现明显拉伸部分的对应数值置为零，用奇异值分解 方法压制随机噪音，即得到反应介质界面物理参数变化的叠前偏移的角道集。(7)由角道集中读取不同空间位置处成像幅值随角度的变化，可由反演方法获得 地层的岩石物理参数，岩石物理参数图像即能指示地下构造的含油气情况。所述的对延拓后的炮点和检波点波场应用反褶积成像条件进行成像是这样实现 的令Ttl是地震记录的时长，则频率采样Δ ω =231/%(弧度/秒)；设ΔΖ是深度方向的 空间采样，η是深度上采样点的顺序数，χ和y是水平位置坐标，m是频率采样点的顺序数； 若炮点和检波点波场经深度延拓后在η Δ ζ深度上分别为PD(x，y，η Δ ζ，πιΔ ω)和Pu(χ，y，
5
A(x, y, ηΔζ, ι Δ ω) = Pd (χ，y，η Δ ζ, m Δ ω) {PD (χ，y，η Δ ζ, m Δ ω)} *式中上标*是复共额，炮点和检波点波场Pd和Pu的单位依据检波器类型可以是 速度或压力，频率的单位是弧度/秒，深度和水平坐标的单位为米。令Aci(IIAz)是A(x，y， ηΔζ,πιΔ ω)在深度 η Δ ζ 上的平均值，当 A (χ，y，η Δ ζ，m Δ ω)彡 1· 2AQ (η Δ ζ)，计算Β(Χ，y, mAc) = ^ ("，兄 ―幼 ^^ (X,少，^
Α(χ,γ,ηΑζ, ηΑω)当 A(x，y，nAz，mA ω) < 1· 2Α0 (η Δ ζ)，计算取,，—▲’-<(2-口)(1+ 0-妁，其
Α(χ,ν,ηΑζ, πΑω) 中α =-
τ令地震资料有效频带的下、上限分别是ω，这里I1和I2是正整数，
则空间各位置处反褶积成像条件的成像幅值Φ(Χ，Υ，ηΔζ)为
kφ(χ,少，《Δζ) = ^ Κο{βΒ{χ,γ,nlSz,πι ^ω)}
m-l、若j是单位虚数，上式中做三维偏移计算时β = jmAco，二维偏移计算时
β = ^ηιΑω ；符号Re代表求复波场值的实部。 所述的利用单程波方法计算炮点的出射波对地下各成像点的入射角度是这样实 现的令X和y是水平坐标，ζ是深度坐标，震源的主频是ω p (弧度/秒)，震源的空间位置 为(&，^，0)，介质的速度为(0^^，2)；利用波动方程偏移的单程波算法，在频率《￡)对空间 脉冲δ (x-xs, y-ys)进行深度延拓，得到以Δζ为间距的各深度离散点上的复波场值Ρ(χ， y，nAz，ωρ, xs, ys)；分别计算
I0 (x，y, ηΔζ) = Re {P (χ, y, ηΔζ, ωρ, xs, ys)}
Z10,乃nlSz) = Re<jy c(x，少p〔x,少,nAz,ω ,xs,) ！
c0P
I2{x,y,nlsz)=
2
-Γ/, (χ, y, nlSz + Δζ) - /, (χ, y, ηΛζ — Δζ)1
3Αζ上式中Xs，[和(^是参数，故在左端的函数中将其忽略；符号I I代表求绝对值， η是深度上采样点的顺序数；速度c(X，y，z)的单位是米/秒，频率的单位是弧度/秒，深度 和水平坐标的单位是米，空间脉冲的单位是压强。若I2(x，y，nAz) > 1. 05I0 (χ, y, η Δ ζ), 则令12&，7，1^2) =0。令 φ。(η Δ ζ)是 IQ(X，y，nAz)在深度 η Δ ζ 上的平均值，当 I。(χ， y, ηΔζ)彡 1· 2Φ0(ηΔζ)，计算
---- [/, (χ, γ, nlSz + 2Δζ) - /, (χ, γ, ηΑζ - 2Δζ)]
12Δζ φ、χ,γ,η^ζ)..
I2(x,y,nAz)
I0(x,y,nAz)当 I0 (χ, y, ηΔζ) < 1. 2 Φ。(η Δ ζ)，计算
φ{χ, y, πΔζ) =，)(2 - (1 + (1 - )，其中《 = '。
对空间各点的
戎 ΟΔζ)武 ΟΔζ)
Φ (χ, γ,ηΔζ)值进行如下处理(1)沿深度方向进行中值滤波；(2)进行空间三维平滑；(3)对每一深度ηΔζ，选取最大值，用这一值除该深度的每一点的数值，实现该深 度上的最大值归一化。令上述处理后的值为i(x，y，ηΔζ)，则可得震源对地下各成像点的入射角度 θ(单位为度)为—θ = acos(<i) (χ, y, ηΔζ))式中(x，y，ηΔζ)即是成像点的空间坐标。。所述的依据得到的各炮点的出射波的入射角度，将由反褶积成像条件得到的各个 单炮偏移结果按角度累加并补偿非均勻覆盖对幅值的影响，形成角道集是这样实现的定 义角度采样间距Δ θ (度)，对每一成像点定义一个一维数组，数组长度小于90/θ Δ ；对全 部炮循环，对每炮的单炮偏移结果，读取由步骤(4)得到的该炮点的出射波对各成像点的 入射角度9，计算1 = 1竹肚{0/^ θ }，将该成像点的幅值累加到对应数组的第1个元素 上，同时统计累加到该元素上的炮数；完成对全部炮循环，将同一水平位置的不同深度点的 成像幅值数组按深度、角度形成二位数组，即得到未经非均勻覆盖补偿的角道集；用统计得 到的累加炮数去除各成像点的不同角度的数值，就得到完成非均勻覆盖补偿的角道集。本发明的波动方程深度偏移直接产生角道集方法，能直接由分炮的炮域偏移结 果，经局部叠加得到地下反射点的叠前偏移的角道集。本发明的波动方程深度偏移直接产生角道集方法，能剔除复杂上覆地层的影响， 得到直接反应目的层处物性参数变化的叠前偏移角道集。本发明的波动方程深度偏移直接产生角道集方法，能补偿观测系统非均勻覆盖对 幅值的改变。本发明的波动方程深度偏移直接产生角道集方法，与常规的炮域波动方程叠前深 度偏移方法相比，不增加太多计算量。本发明的波动方程深度偏移直接产生角道集方法，可应用于二维和三维波动方程 叠前深度偏移处理。本发明的波动方程深度偏移直接产生角道集方法，能用于反演地层的岩石物理参 数，直接得到地下构造的含油气情况。本发明的具体实现原理如下本发明的核心思想是直接用单程波方法计算炮点的出射波在地下成像点的入射 角度，避免了射线法对速度模型的不稳定也与波动方程偏移方法更匹配。入射角度计算仅 利用单程波的简谐波场，与偏移算法相比计算量可以忽略。计算入射角度的原理是，由单程波方程可知，下行的频率域波场户+(^，2;劝满足 「■91 δΡ+(χ,γ,ζ·ω) ~
L0052」---= ~]kzP (χ,γ, ζ; ω)(1)
οζ 其中j是单位虚数，kz是垂直波数，ω是角频率，χ和y是水平位置坐标，ζ是深度坐标，由频散关系得 [_] 式中c (X，y，ζ)是空间各位置处的介质波速，角度θ即是空间位置(X，y，ζ)处波 前面法线与深度轴ζ的夹角。将(2)式代入(1)式并做离散的时域付氏反变换，且令时间 t = 0，可得
(3)
中符号Re 代表求复波场值的实部。若令
)，其中
f(")是地震子波的谱，则(3)式中
对应于震源点置于(xs，ys，0)处时
ω>0
的脉冲响应，而波前面法线就是(Xs，ys，0)处炮点的出射波在(x，y，z)点的入射波方向。脉冲响应
在空间中仅存在一个波前面，对没有波前面的空间
ω>0
位置，不能用波前面来确定炮点入射波的方向。若仅对主频分量计算脉冲响应，即仅计算 Re
则波前面将分布于整个区域，易于求得所有成像点的入射波方向，而此
时计算脉冲响应的计算量也大幅减少。仅考虑主频分量，由(3)式可求得 式(4)中的分母将存在多个零点，此零点处将不能得到正确的入射波角度；既然 在一个波长范围内入射波的角度不会有显著变换，可用波前面的波锋和波谷处的入射波角 度插值得到邻近位置的入射波角度。我们了设计如下平滑算法，实现这一思想。首先计算 用四阶差分算法计算(4)式中的一阶导数，令计算脉冲响应的深度方向的空间采 样为Δζ，可得
式(5)和(7)中采用绝对值，是因为对下行波而言cos θ总是正。由⑷式可知， I2(x，y，z)应该大于Itl(^yd)。但当偏移噪音较强时(即在入射波较弱的阴影区)，有时 并不满足这一条件；这表明该处不存在准确的入射波，因此不能得到正确的入射角度。在计 算中，若I2(x，y，z) > 1.051(10^，7，2)，可令120^，7，2) = 0，这样避免了阴影区角度计算的困难。 为避免式(4)中除法计算的不稳定，我们发展了高精度的稳定的相除算法，这-算法的核心是利用级数展开近似式 当χ趋于零，式⑶仍可得到稳定的结果，而当|x|≤1.2时，式⑶的近似有很 好的精度。利用式(8)，我们发展了稳定的相除算法， 式中Φ。(ζ)是I。(x，y，z)在深度z上的平均值。对(9)式的结果，先做沿深度方向的中值滤波，然后进行空间三维平滑，即可由波 锋和波谷处的入射波角度插值得到邻近位置的入射角度；为避免(7)式差分等计算产生的 误差，将对每一深度上平滑处理后的Φ (χ,γ,ζ)进行归一化处理，保证垂直入射的存在。这 样就可通过单个频率的简谐波场的深度延拓，求得炮点出射波对地下各点的入射角度
θ = acos(i(x，y, z)) (10)式中(x，y，z)即是地下各点的空间坐标。本发明对炮域偏移采用反褶积成像条件，而不是采用相关成像条件；反褶积成像 条件可补偿地震波的扩散效应，获得正确的幅值。我们再次利用式(8)来解决反褶积成像条件中除法的稳定问题。炮域偏移的反褶 积成像条件可表达为 式中符号Re代表求复波场值的实部，Δ ω是频率采样，m是频率采样点的顺序数， PuU, y，z，mΔ ω)和PD(X，y，z，mΔ ω)分别是深度延拓后深度ζ上频率域的炮点和检波点 波场J1和I2是正整数，I1 Δ ω和I2 Δ ω分别对应地震资料有效频带的下、上限；引入参数 β是为了保持成像的波形与检波点波场相一致，若j是单位虚数，做三维偏移计算时β = jm Δ ω，二维偏移计算时β= √jm△ω。式(11)中的除法可进一步表达为 式中上标*代表复共额。令式(12)中实数的分母为六0^，7，2，11^(0)，由⑶式
可得
Α{χ,γ,ζ,ηιω)
1
Α(χ,γ,ζ,ιηω)
Α(χ, γ, ζ, ιη ) > 1(ζ)
Λ(ζ)
2-
Α(χ,γ,ζ,Μω)Λω.
1 +
f
(13)式中Atl(Z)是A(X，y，z，mco)在深度ζ上的平均值。将式(13)代入(12)再代入 (11)，可得稳定的反褶积成像条件。反褶积成像条件得到的偏移幅值，反应了成像点处炮点波场的入射角度所对应的 正确反射强度。在形成角道集时，尽管角度的采样是相同的，但不同角度上累加的炮数是不 同的，角度越大，累加的数量越多；而炮点覆盖越密，累加的数量也就越多。为获得正确地反 应成像点物理参数对比的成像幅值，我们在累加过程中统计各角度间距内累加的炮数，用 这一数量去除累加结果。这样即考虑了成像结果在角度域分布的不均勻，也补偿了炮点覆 盖不均勻的影响。本发明的有益效果该方法可应用于复杂上覆地层下油气储层的地震偏移成像， 能生成直接反应目的层物理参数变化的叠前偏移的角道集。该角道集能更好地服务于叠前 反演等油气和流体直接检测技术(避免了各种基于测井资料的补偿)，降低了无井或少井 地区油气勘探的风险。该方法将波动方程叠前深度偏移由构造成像推进到岩性识别，拓展 和提升了波动方程叠前深度偏移的应用范围和效果。该方法对存在复杂上覆层的大量区域 的油气、矿产资源勘探有重要应用价值。
图1是大庆LMD地区开发地震的典型单炮记录，采集是16线，双边观测，道间距 40m ；本次采用波动方程偏移直接生成角道集方法共使用了 3600炮，满覆盖11. 2平方公里。图2a是叠前深度偏移使用的三维速度模型的沿测线322的二维切片，图2b是过 ⑶P400的沿垂直测线方向的二维速度切片。图3是线号322和⑶P号400处由所有单炮偏移结果构成的共成像点道集，单炮 偏移采用了本发明的反褶积成像条件。图4a和图4b是炮点在线号318和⑶P号415处的炮点波场的入射角分布图，入 射角计算采用了本发明的新方法；图4a是沿测线方向的二维角度切片，图4b是沿垂直测线 方向的二维角度切片。图5是线号318和⑶P号415处的未进行非均勻覆盖补偿的角道集，这一道集是 依据各炮点在线号318和CDP号415处的入射角度，由各单炮偏移结果按角度做局部叠加 得到。图6是对图5角道集应用非均勻补偿后的结果，对比图5和图6可知，小角度的幅 值得到增强。图7是对图6角道集首先应用剩余动校，再将出现明显拉伸部分的对应数值置为 零，然后用奇异值分解方法压制随机噪音得到的最终角道集。图8a是图7角道集的局部放大图，图8b是在两个典型层位处的随角度变化的幅 值对比。图8a的局部放大图对应了有测井资料的地层部分。图8b的对比图中的实线是由图7读出的深度970m和1030m处幅值随角度的变化，虚线是依据该位置(线号318和⑶P 号415)处声波和密度测井资料读取的P波、S波速度以及密度，由理论计算得到的随角度 变化的反射强度。从图8b中对比可见两者拟合很好，这表面本发明求得的角道集不需“校 正”即可较好地拟合反射界面的入射角与反射幅值关系(AVA)特征。图9a是线号323和⑶P号360的井位处的角道集，图9b是角道集的局部放大，图 9c是两个典型层位(深度840m和975m)处幅值随角度变化与理论AVA特征的对比。各图 的意义与图7、图8a和图8b相同。图IOa是线号363和⑶P号448的井位处的角道集，图IOb是角道集的局部放大， 图IOc是两个典型层位(深度995m和1095m)处幅值随角度变化与理论AVA特征的对比。 各图的意义与图7、图8a和图8b相同。图Ila和图lib是本次偏移处理得到的大庆LMD区块偏移叠加剖面。图Ila是沿 测线322的二维切片，图lib是过⑶P400的沿垂直测线方向的二维切片。图Ila和图lib 是由各成像点的角道集沿角度方向叠加得到的。
具体实施例方式实施例1 波动方程深度偏移直接产生角道集方法，针对大庆LMD地区为例，具体 为以下步骤(1)用双边采集方式，由检波器获取人工震源激发的反射地震信号，记录到磁带 上，每炮布设16条采集线，采集线上的道间距为40m，每炮共含16X168 = 2688道；共激发 和记录3600炮，图1是典型的单炮记录。(2)从磁带上读取反射地震信号，按叠前炮集排序存放到PC-Cluster计算机中， 偏移计算将通过读取每个单炮记录，对炮循环来进行。(3)依据图2a和2b的偏移速度模型，利用波动方程偏移的单程波算法分别将炮点 的震源波场(即空间和时间脉冲函数)和炮记录(接收到的地震信号)沿深度轴进行波场 延拓，对延拓后的炮点和检波点波场应用反褶积成像条件进行成像。令Ttl是地震记录的时长，则频率采样Δω = 231/% (弧度/秒)；设Δ Z是深度 方向的空间采样，η是深度上采样点的顺序数，χ和y是水平位置坐标，m是频率采样点的 顺序数；若炮点和检波点波场经深度延拓后在η Δ ζ深度上分别为Pd (X，y，η Δ z，m Δ ω)和 Ρυ(χ，y，ηΔζ, ι Δ ω)，计算A (χ, y,nAz,mA ω) = Pd (χ，y，η Δ ζ, m Δ ω) {PD (χ，y，η Δ ζ, m Δ ω)} * 式中上标 * 是复共额，炮点和检波点波场Pd和Pu的单位依据检波器类型可以是速度或压力，频率的单 位是弧度/秒，深度和水平坐标的单位为米。令Atl (ηΔζ)是A(x，y，ηΔζ, mA ω)在深度 ηΔζ 上的平均值，当 A(x，y，nAz，mA ω)彡 1. 2AQ (η Δ ζ)，计算 当 A(x，y，ηΔζ, ι Δ ω) <1. 2Α0(ηΔζ)，计算Β(χ,y,Mz,>ηΑω) = ^(Χ,少, ^比(Χ，兄成^ (2-α)( + (1 -^)，其
A0(^Az)\'
中
11
^=“令地震资料有效频带的下、上限分别是I1A ω和I2 Δ ω，这里I1和I2是正整数，
则空间各位置处反褶积成像条件的成像幅值Φ(Χ，Υ，ηΔζ)为
hΦΟ，y, nAz) = Z Re {々5(χ, γ, ηΑζ, πι ^ ))
m=l\若j是单位虚数，上式中做三维偏移计算时β = >Δω，二维偏移计算时 β = ^)ηιΑω ；符号Re代表求复波场值的实部。将成像结果φ (χ,γ,ηΔζ)存储在计算机的硬盘中。图3是线号322和⑶P号400 处由所有单炮偏移结果构成的共成像点道集；所有单炮偏移采用的参数均是Δζ = 5m，深 度上的采样点总数η = 1000, Δω= 1. 40 (弧度/秒)，I1 = 13，I2 = 337。(4)对每个炮点(震源)，依据图2a和2b的偏移速度模型，利用单程波方法对频 率为25Hz的空间脉冲进行深度延拓，取ΔΖ = 5m，求得炮点的出射波对地下各成像点的入 射角度。计算过程是令χ和y是水平坐标，Z是深度坐标，震源的主频是ω p (弧度/秒)， 震源的空间位置为(Xs，^，0)，介质的速度为(3(1 y, ζ)；利用波动方程偏移的单程波算法， 在频率ωρ对空间脉冲δ (x-xs,y-ys)进行深度延拓，得到以Δ ζ为间距的各深度离散点上 的复波场值？0^，7，“2， p,xs,ys)；分别计算I0 (X，y, ηΔζ) = Re {P (χ, y, ηΔζ, ωρ, xs, ys)}/,(χ,少，nAz) = Rei) C“,y,nAZ) y,n^z,ω ,xs
ωΡJ
2I2 (χ, y, nAz)=——[/, (χ, y, nAz +Az)- Ix (χ, γ, πΔζ - Δζ)]
3Αζ——-[Z1 (χ, y, nlSz + 2Δζ) - /, (χ, γ, ηΑζ - 2Δζ)]
12Δζ上式中Xs，[和(^是参数，故在左端的函数中将其忽略；符号I I代表求绝对值， η是深度上采样点的顺序数；速度c(X，y，z)的单位是米/秒，频率的单位是弧度/秒，深度 和水平坐标的单位是米，空间脉冲的单位是压强。若I2(x，y，nAz) > 1. 05I0 (χ, y, η Δ ζ), 则令12&，7，1^2) =0。令 φ。(η Δ ζ)是 IQ(X，y，nAz)在深度 η Δ ζ 上的平均值，当 I。(χ， y, ηΔζ)彡 1· 2Φ0(ηΔζ)，计算φ(χ,γ,ηΑζ) = —-—
I0(x,y,nAz)当I0 (χ, y, η Δ ζ) < 1. 2 Φ。(η Δ ζ)，计算敝,_)=^(7，)(2一。)(1+ (1 —吶，其中0 = ^^对空间各点的Φ (χ, γ,ηΔζ)值进行如下处理1)沿深度方向进行中值滤波；2)进行空间三维平滑；
3)对每一深度η Δ ζ，选取最大值，用这一值除该深度的每一点的数值，实现该深 度上的最大值归一化。令上述处理后的值为i(x，y, ηΔζ)，则可得震源对地下各成像点的入射角度 θ(单位为度)为—θ = acos(<i) (χ, y, ηΔζ))式中(x，y，ηΔζ)即是成像点的空间坐标。图4a和图4b是计算得到的某炮点出射波的空间入射角分布图。(5)读取各炮的单炮偏移结果，根据计算得到的该炮点出射波的入射角度值，将偏 移结果按角度累加并补偿非均勻覆盖对幅值的影响，形成角道集，角度采样间距取为Δ θ =2° 0计算过程是定义角度采样间距Δ θ (度)，对每一成像点定义一个一维数组，数 组长度小于90/Δ θ ；对全部炮循环，对每炮的单炮偏移结果，读取由步骤(4)得到的该炮 点的出射波对各成像点的入射角度9，计算1 = 1竹#{9/^ Θ}，将该成像点的幅值累加 到对应数组的第1个元素上，同时统计累加到该元素上的炮数；完成对全部炮循环，将同一 水平位置的不同深度点的成像幅值数组按深度、角度形成二位数组，即得到未经非均勻覆 盖补偿的角道集；用统计得到的累加炮数去除各成像点的不同角度的数值，就得到完成非 均勻覆盖补偿的角道集。图5是某水平位置处未进行非均勻覆盖补偿的角道集。图6是某水平位置处完成 非均勻覆盖补偿的角道集。(6)对全部角道集循环，做剩余动校，将出现明显拉伸部分的对应数值置为零，用 奇异值分解方法压制随机噪音，图7是完成后续处理的角道集。(7)本领域技术人员能由角道集中读取不同空间位置处幅值随角度的变化，并由 反演方法获得地层的岩石物理参数，岩石物理参数图像即能指示地下构造的含油气情况。 参阅图8a和8b，图9a、9b和9c，图10a、10b和10c，和图Ila和lib。
权利要求
一种波动方程深度偏移直接产生角道集方法，其特征在于采用以下步骤A)用单边或双边采集方式，由检波器获取人工震源激发的反射地震信号，记录到磁带上；B)从磁带上读取反射地震信号，按叠前炮集排序存放到计算机中；C)利用波动方程偏移的单程波算法分别将炮点的震源波场和检波点的接收波场沿深度轴进行波场深度延拓，对延拓后的炮点和检波点波场应用反褶积成像条件进行成像；D)利用单程波方法计算炮点的出射波对地下各成像点的入射角度；E)依据步骤D得到的各炮点出射波的入射角度，将由步骤C得到的各个单炮偏移结果按角度累加并补偿非均匀覆盖对幅值的影响，形成角道集；F)对角道集做剩余动校，将出现明显拉伸部分的对应数值置为零，用奇异值分解方法压制随机噪音，即得到反应介质界面物理参数变化的叠前偏移的角道集；G)由角道集中读取不同空间位置处成像幅值随角度的变化，可由反演方法获得地层的岩石物理参数，岩石物理参数图像即能指示地下构造的含油气情况。
2.根据权力要求1所述的一种波动方程深度偏移直接产生角道集方法，其特征在于 所述的对延拓后的炮点和检波点波场应用反褶积成像条件进行成像是这样实现的令Ttl是 地震记录的时长，则频率采样Δ ω =231/%(弧度/秒)；设Δζ是深度方向的空间采样，η 是深度上采样点的顺序数，χ和y是水平位置坐标，m是频率采样点的顺序数；若炮点和检波 点波场经深度延拓后在ηΔζ深度上分别为PD(x，y，nAz，mA ω)和Pu (x，y，η Δ z，m Δ ω)， 计算A(x, y, η Δ z, mA ω) = PD(x, y, ηΔζ, mA ω) {PD(x, y, η Δ ζ, mA ω)}*式中上标*是复共额，炮点和检波点波场Pd和Pu的单位依据检波器类型可以是速度或 压力，频率的单位是弧度/秒，深度和水平坐标的单位为米。令Aci(HAz)是A(x，y，ηΔζ, mA ω)在深度ηΔζ上的平均值，当A(x，y，nAz，mA ω)彡1· 2Α。(η Δ ζ)，计算Pu (χ，y, ηΑζ, ηιΑω) { Pd (χ, γ, ηΑζ, τηΑω)}Β(χ,γ ηΑζ,τηΑω)=---—Α(χ,γ,ηΑζ,τηΑω)当 A(x，y，nAz，mA ω) < 1· 2Α。(η Δ ζ)，计算Α(χ,γ,ηΑζ,τηΑω)α —-Ao(^Az)令地震资料有效频带的下、上限分别是I1 Δ ω和I2 Δ ω，这里I1和I2是正整数，则空间各位置处反褶积成像条件的成像幅值Φ (x, y，ηΔζ)为 丨1φ(χ, γ, nAz) =βΒ(χ, y, nAz, πι ^ )]若j是单位虚数，上式中做三维偏移计算时β = jmA ω，二维偏移计算时/J = TI^; 符号Re代表求复波场值的实部。
3.根据权力要求1所述的一种波动方程深度偏移直接产生角道集方法，其特征在于 所述的利用单程波方法计算炮点的出射波对地下各成像点的入射角度是这样实现的令 X和y是水平坐标，Z是深度坐标，震源的主频是ωρ(弧度/秒)，震源的空间位置为(xs， ^，0)，介质的速度为(3 0^，y, ζ)；利用波动方程偏移的单程波算法，在频率ωρ对空间脉冲δ (x-xs, y-ys)进行深度延拓，得到以Δζ为间距的各深度离散点上的复波场值P(X，y， ηΔζ, ωρ,χ3, ys)；分别计算I0 (X，y, ηΔζ) = Re {P (χ, y, ηΔζ, ωρ, xs, ys)}Z1 (X’y,nAz) = Rejy c(x，:，—)尸(χ’^nhZ,ωρ,χχ’兄.)j 2I2 (χ, γ, ηΑζ) =-[I1 (χ，y, ηΑζ + Δζ) - Ix (χ, γ, ηΑζ - Δζ)]3Δζ———[Z1 (χ, γ, ηΑζ + 2Δζ) - Z1 (χ, γ, ηΑζ - 2Αζ)] \2Αζ上式中Xs，[和《5是参数，故在左端的函数中将其忽略；符号I I代表求绝对值，η是 深度上采样点的顺序数；速度c (X，y, ζ)的单位是米/秒，频率的单位是弧度/秒，深度和 水平坐标的单位是米，空间脉冲的单位是压强。gI2(x，y，nAz) > 1. 05、(x，y，η Δ z)，则 令 I2(x，y，nAz) = 0。令 Φ。(ηΔζ)是 I。(x，y，η Δ ζ)在深度 η Δ ζ 上的平均值，当 IQ(x，y， ηΔζ) ^ 1. 2Φ。(ηΔζ)，计算ntsz) = —-τ-Α (χ. ν. ηΑζ)Ov- ^yj/当 I0 (χ, γ,ηΔζ) < 1.2Φ。(ηΔζ)，计算树χ，少,^Δζ) = 7 舡少,——Ω) + (1 _ …道中口 = W，—对空间各点的Φ (χ，y，ηΔζ)值进行如下处理(1)沿深度方向进行中值滤波；(2)进行空间三维平滑；(3)对每一深度ηΔ ζ，选取最大值，用这一值除该深度的每一点的数值，实现该深度上 的最大值归一化。令上述处理后的值为i(X，y，nAz)，则可得震源对地下各成像点的入射角度θ (单位 为度)为θ = acos (Φ (χ, y, η Δ ζ))式中(x，y，nAz)即是成像点的空间坐标。
4.根据权力要求1所述的一种波动方程深度偏移直接产生角道集方法，其特征在于 依据步骤D得到的各炮点出射波的入射角度，将由步骤C得到的各个单炮偏移结果按角 度累加并补偿非均勻覆盖对幅值的影响，形成角道集是这样实现的定义角度采样间距 Δ θ (度)，对每一成像点定义一个一维数组，数组长度小于90/Δ θ ；对全部炮循环，对每 炮的单炮偏移结果，读取由步骤D得到的该炮点的出射波对各成像点的入射角度θ，计算 1 = l+int{ θ / Δ θ }，将该成像点的幅值累加到对应数组的第1个元素上，同时统计累加到 该元素上的炮数；完成对全部炮循环，将同一水平位置的不同深度点的成像幅值数组按深 度、角度形成二位数组，即得到未经非均勻覆盖补偿的角道集；用统计得到的累加炮数去除 各成像点的不同角度的数值，就得到完成非均勻覆盖补偿的角道集。
全文摘要
波动方程偏移直接产生角道集方法，应用于地震勘探中地震资料处理，是用于形成地震偏移的角道集的波动方程叠前深度偏移方法。该方法通过计算震源出射波在地下成像点的入射角度，直接由炮域偏移结果生成地下反射界面成像的角道集，且补偿了观测系统非均匀覆盖对幅值的改变。角道集直接对应于目的层处物性参数变化，为叠前反演提供了准确的角度域地震反射强度。这一方法的核心是直接用单程波算法确定震源波场在成像点的入射角度，避免了射线法的不稳定也与波动方程偏移方法更匹配。入射角度计算仅利用单程波的简谐波场，与偏移计算相比计算量可以忽略。该方法能应用于二维和三维波动方程叠前深度偏移处理，对油气、矿产资源勘探有重要应用价值。
文档编号G01V1/28GK101923175SQ20091023821
公开日2010年12月22日 申请日期2009年11月17日 优先权日2009年11月17日
发明者刘礼农, 张剑锋, 张辉 申请人:中国科学院地质与地球物理研究所