CN102141633A - 各向异性三维叠前时间偏移方法 - Google Patents

Abstract

各向异性三维叠前时间偏移方法，应用于地震勘探中反射地震资料处理，是针对三维采集地震资料的叠前偏移成像方法。该方法考虑地球介质速度各向异性对地震波传播的走时和幅值的影响，能在偏移过程中自主决定偏移速度和各向异性参数，因而能得到归位准确、保幅的偏移图像。该方法根据地下构造时变的倾角确定各向异性介质中的三维时变偏移孔径，且能在偏移过程中同时压制偏移噪音。该方法通过在集群计算机的计算节点间合理分配地震资料，实现高效的并行计算。该方法的核心是应用深度偏移的单程波算子和稳相点原理求得各向异性介质中地震波的走时、幅值、成像权系数以及入射角度。该方法对油气、矿产资源勘探有重要应用价值。

Description

各向异性三维叠前时间偏移方法

技术领域

本发明属于油气田和矿产资源地震勘探中反射地震资料处理技术领域，涉及地震资料处理过程中的叠前偏移成像技术范畴，是三维地震资料偏移成像的各向异性三维叠前时间偏移方法。

背景技术

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

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

实验室测试及场地数据表明，许多沉积岩具有各向异性；在偏移算法中忽略介质的各向异性，将导致错误的偏移归位和丢失(或摸糊)某些 地质构造。实际地层的速度各向异性，可用横向各向同性较好地近似。已有一些叠前时间偏移方法在计算地震波走时考虑各向异性，但尚没有准确考虑各向异性对地震波幅值的影响，以及如何在各向异性情况下实现保幅成像。

既然在计算地震波走时与幅值时需考虑介质的各向异性，必须研究如何在叠前时间偏移方法中引入各向异性参数，以及如何估算与偏移算法对应的各向异性参数。类似于各向同性介质中使用叠加速度进行偏移计算，采用叠加(等效)各向异性参数是一个主流的选择，这一方法使得各向异性参数建模变得更容易。

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

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

偏移孔径对叠前时间偏移是重要的，较好的偏移孔径对压制偏移噪音和降低偏移计算量有更重要作用。小的偏移孔径可减少偏移计算量，但存在着不能对陡倾角构造正确成像的风险；过大的孔径又带来了偏移噪音和较大的计算量。偏移孔径的选取与地震波的入射角度有关，地球介质各向异性改变了地震波的传播路径，相应地改变了入射角度，因此在确定偏移孔径时需考虑介质的各向异性。对三维地震资料而言， 由于采集信号在沿测线和垂直测线方向有不同的空间采样率，地震资料存在方位角的变化，偏移孔径的选取较二维情况变得更复杂。

发明内容

本发明的目的是：提供一种各向异性介质中的三维叠前时间偏移方法，它通过在偏移过程中自主决定偏移速度和各向异性参数，正确考虑了地球介质各向异性对地震波走时和幅值的影响；通过给出保幅成像的权系数，实现了保幅偏移；通过给出考虑各向异性的三维时变偏移孔径和对偏移孔径边缘区的幅值施加衰减，实现了陡倾角构造成像和压制偏移噪音；通过在并行计算的各个进程之间合理分配叠前地震资料，实现了高加速比的并行计算。

本发明采用的技术方案是：各向异性三维叠前时间偏移方法，具体步骤包括：

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

(2)从磁带上读取地震信号，对叠前地震资料进行常规的压制噪音处理。针对部分共中心点，抽取共中心点道集，对抽取的道集作常规的NMO动校正速度拾取，对所得结果做横向平均，作为初始偏移速度；

(3)将叠前地震资料按偏移距大小排序，对地震资料做偏移孔径相关的相等计算量分组，将不同组地震资料存放到集群计算机的不同计算节点上；

(4)利用初始偏移速度，对已存放到集群计算机各个计算节点上的地震资料，应用时变偏移孔径的三维各向同性叠前时间偏移方法，进行并行的各向同性偏移计算，收集各计算节点的偏移结果，形成共反射点道集；

(5)依据形成的共反射点道集，确定各向异性偏移使用的偏移速度和各向异性参数；

(6)根据地下构造在不同时间深度的最大倾角、偏移速度、各向异性参数和地震道的偏移距和方位角，确定时变的三维偏移孔径和对应于偏移孔径的、用于压制偏移噪音的衰减系数；

(7)采用基于单程波算子和稳相点原理的查表法求得各向异性介质中的地震波走时和保幅成像的权系数；

(8)在每个计算节点，对地震资料的所有地震道循环。对每一地震道，由偏移孔径确定有效的成像区域和区域中各成像点偏移计算的起始成像时间，利用保幅成像的权系数，从起始成像时间开始计算各成像点的偏移幅值。对偏移孔径边缘区域的成像点，进一步用衰减系数来衰减偏移幅值。将偏移幅值累加到存放偏移结果的数组中对应的偏移距上；

(9)收集各计算节点的偏移结果，对形成的全部共反射点道集做常规的各向同性剩余动校正，将出现明显拉伸和噪音部分的对应数值置为零，将不同偏移距的偏移结果叠加，形成偏移叠加剖面；

(10)通过显示软件将偏移叠加剖面数值转换为地下反射构造的剖面图像，剖面图像将指示地下构造的形态、断裂部位、断距大小和地层沉积样式及地层的波阻抗特征，用于确定地下生、储油构造和识别油气储层。

所述的将叠前地震资料按偏移距大小排序，对地震资料做偏移孔径相关的相等计算量分组，将不同组地震资料存放到集群计算机的不同计算节点上是这样实现的：设地下构造在沿测线和垂直测线方向的两个最大倾角为θ_x和θ_y，已得到的横向均匀的初始偏移速度为 计算两个长度值

$a\left(T\right)=\frac{{\stackrel{\‾}{V}}_{\mathrm{rms}}\left(T\right)\·T}{\sin \left(2{\mathrm{\θ}}_x\right)}\sqrt{1+{\cos }^2\left({2\mathrm{\θ}}_x\right)+{\mathrm{\λ}}^2{\sin }^2\left({2\mathrm{\θ}}_x\right)-2\cos \left({2\mathrm{\θ}}_x\right){[1+{\mathrm{\λ}}^2{\sin }^2\left({2\mathrm{\θ}}_x\right)]}^{1/2}}$

$b\left(T\right)={\stackrel{\‾}{V}}_{\mathrm{rms}}\left(T\right)\·T{\tan \mathrm{\θ}}_y$

式中T是用单程旅行时表达的时间深度，单位是秒，h是地震道的半偏移距，单位是米， 是无量纲参数。时间深度T处的成像区域可看作以a(T)和b(T)为长、短半轴的椭圆面。令偏移计算的起始时间深度为T₀，最大时间深度为T₂，而大于T₁偏移孔径将不随时间深度变化，则偏移距为2h的地震道的相对偏移计算量可近似为

$G\left(h\right)=a\left(T_1\right)\·b\left(T_1\right)(T_2-T_1)(1/\mathrm{\ΔT})+\underset{j=1}{\overset{n}{\mathrm{\Σ}}}a(T_0+\mathrm{j\ΔT})\·b(T_0+\mathrm{j\ΔT})$

式中n＝(T₁-T₀)/ΔT是一整数，ΔT是时间深度方向上成像的采样间距，a(T₁)和b(T₁)等可由前一个公式计算得到。设地震资料的最大半偏移距为h_max，最小半偏移距为h_min，偏移距采样间距为Δh，可用整数i＝(h-h_min)/Δh+1作为偏移距的索引；统计有相同i的地震道，可得到偏移距满足h_min+(i-1)Δh≤h＜h_min+iΔh的地震道道数m_i(i＝1，l)，其中l＝(h_max-h_min)/Δh+1。记所有m_i中最大的为m_j，可得无量纲的权系数ρ_i＝G(h_min+(i-1)Δh)/G(h_min+(j-1)Δh)。

如集群计算机共有k个核，即希望同时有k个进程并行计算，则计算 

可按如下方式确定各进程中包含的地震资料的偏移距索引：对进程1，求满足下式的n₁，

$\underset{i=1}{\overset{n_1}{\mathrm{\Σ}}}{\mathrm{\ρ}}_i{m}_i\≤m_0\≤\underset{i=1}{\overset{n_1+1}{\mathrm{\Σ}}}{\mathrm{\ρ}}_i{m}_i$

若

$m_0-\underset{i=1}{\overset{n_1}{\mathrm{\Σ}}}{\mathrm{\ρ}}_i{m}_i>\underset{i=1}{\overset{n_1+1}{\mathrm{\Σ}}}{\mathrm{\ρ}}_i{m}_i-m_0$

取n₁＝n₁+1，索引数为i＝1，n₁的地震道将分配给进程1；对进程2，求满 足下式的n₂，

$\underset{i=n_1+1}{\overset{n_2}{\mathrm{\Σ}}}{\mathrm{\ρ}}_i{m}_i\≤m_0\≤\underset{i=n_1+1}{\overset{n_2+1}{\mathrm{\Σ}}}{\mathrm{\ρ}}_i{m}_i$

若

$m_0-\underset{i=n_1+1}{\overset{n_2}{\mathrm{\Σ}}}{\mathrm{\ρ}}_i{m}_i>\underset{i=n_1+1}{\overset{n_2+1}{\mathrm{\Σ}}}{\mathrm{\ρ}}_i{m}_i-m_0$

取n₂＝n₂+1，索引数为i＝n₁+1，n₂的地震道将分配给进程2；如此类推，完成全部k个进程的地震资料分配，这样可使得各进程的偏移计算量近似相同。一般集群计算机的一个计算节点可同时进行几个进程的计算，可将相关进程的地震资料一起存放到这个计算节点上。

所述的依据形成的共反射点道集，确定各向异性偏移使用的偏移速度和各向异性参数是这样实现的：对共反射点道集，利用初始偏移速度做反动校，再做动校正得到新的速度，对这一速度做空间平滑处理，作为各向异性偏移的初始速度V_nmo。对反动校后的共反射点道集，再次利用下式进行各向异性动校正

${\mathrm{\τ}}^2=\frac{3+4\mathrm{\η}}{4(1+\mathrm{\η})}H\left(2h\right)+\frac{1}{4(1+\mathrm{\η})}\sqrt{H^2\left(2h\right)+16\mathrm{\η}(1+\mathrm{\η}\frac{{\left(2T\right)}^2{\left(2h\right)}^2}{(1+2\mathrm{\η})V_{\mathrm{nmo}}^2})}$

$H\left(2h\right)={\left(2T\right)}^2+\frac{{\left(2h\right)}^2}{(1+2\mathrm{\η})V_{\mathrm{nmo}}^2}$

式中h是地震道的半偏移距，单位是米，T是用单程旅行时表达的时间深度，单位是秒，η是无量纲的各向异性参数。可根据叠加能量最大和η的空间联续性，即η的横向变化不应过大，确定各空间点的各向异性参数η；再分别取速度为0.95V_nmo和1.05V_nmo进行上式各向异性动校正，(0.95V_nmo，η)、(V_nmo，η)和(1.05V_nmo，η)三组参数中使叠加能量最大、同相轴最平直的那组即是该点最终的偏移速度和各向异性参数，对全部各向异性参数还需进行空间平滑处理。

所述的根据地下构造在不同时间深度的最大倾角、偏移速度、各 向异性参数和地震道的偏移距和方位角，确定时变的三维偏移孔径和对应于偏移孔径的、用于压制偏移噪音的衰减系数是这样实现的：用一组整数 

(i＝1，k)描述三维偏移孔径，其中k是有效成像区域中总的离散点数。令Δx和Δy是偏移成像结果在两个水平坐标方向的空间采样，ΔT是时间深度方向的采样，则 

和 

是成像区域中某离散点与地震道中心点的沿两个水平坐标方向的距离，而 是该点偏移计算的起始成像时间；令地震道中心点坐标为 

和 引入时变的偏移孔径，就是在该地震道的偏移计算中，对成像区域中水平坐标为 

和 

的点(i＝1，k)，仅从 

开始进行偏移计算。

在几个关键深度，定义地质构造在沿测线和垂直测线方向的最大倾角为 和 

构建如下不等式方程组：

式中 

是地震道的方位角。满足上两式的最小角度αⁱ和βⁱ，就是保证该深度最大倾角为 

和 

的地质构造得以正确成像时，在沿方位角方向和垂直方位角方向的两个最大成像角度。

对αⁱ和βⁱ应用分段线性插值，可得到随时间深度变化的、沿方位角和垂直方位角方向的两个最大成像角度α(T)和β(T)。对不同时间深度T，偏移孔径所定义的成像区域，可看作以方位角方向为长轴方向、垂直方位角方向为短轴方向、中心在地震道的中心点上的椭圆面，其长、短半轴a(T)和b(T)分别为：

a(T)＝a₀+Δa

$a_0=\frac{{\stackrel{\‾}{V}}_{\mathrm{nmo}}\left(T\right)\·T}{\sin \left[2\mathrm{\α}\left(T\right)\right]}\sqrt{1+{\cos }^2\left[2\mathrm{\α}\left(T\right)\right]+{\mathrm{\λ}}^2{\sin }^2\left[2\mathrm{\α}\left(T\right)\right]-2\cos \left[2\mathrm{\α}\left(T\right)\right]{\{1+{\mathrm{\λ}}^2{\sin }^2[2\mathrm{\α}\left(T\right)\left]\right\}}^{1/2}}$

$\mathrm{\Δa}=\frac{{5a}_0{\stackrel{\‾}{\mathrm{\η}}}_{\mathrm{eff}}\left(T\right)({\tilde{a}}^2-{\mathrm{\λ}}^2-1)}{2{\sin }^2\left[2\mathrm{\α}\left(T\right)\right]({\tilde{a}}^2-{\mathrm{\λ}}^2+1)-4}$

$b\left(T\right)={\stackrel{\‾}{V}}_{\mathrm{nmo}}\left(T\right)\·T\sin \left[\mathrm{\β}\left(T\right)\right]{\left({\cos }^2\right[\mathrm{\β}\left(T\right)]-\frac{{\stackrel{\‾}{\mathrm{\η}}}_{\mathrm{eff}}\left(T\right){\sin }^2\left[\mathrm{\β}\left(T\right)\right]{\{1+{3\cos }^2[\mathrm{\β}\left(T\right)\left]\right\}}^2}{{\sin }^2\left[\mathrm{\β}\left(T\right)\right]{\{1+2{\cos }^2[\mathrm{\β}\left(T\right)\left]\right\}}^2+{4\cos }^6\left[\mathrm{\β}\left(T\right)\right]})}^{-1/2}$

式中h是地震道的半偏移距，T是用单程旅行时表达的时间深度，a₀是忽略各向异性时的长轴值，Δa是考虑各向异性时长轴的增加量， 是无量纲参数， 

是无量纲参数， 

是不同时间深度上的平均各向异性参数，无量纲， 

是不同时间深度上的平均偏移速度。

设偏移计算的起始时间深度为T₀，最大时间深度为T₂，对T＝T₀，T₂循环，若a(T)≥a(T-ΔT)，a(T)不变，否则a(T)＝a(T-ΔT)；对a(T)作时间深度方向的高斯平滑，即得随时间深度变化的长半轴，同理可得随时间深度变化的短半轴。

对T＝T₀，T₂循环，对时间深度T和T+ΔT，由其对应的两组长、短半轴形成有共同的中心点(0，0)的两个椭圆，这两个椭圆形成的条带就是起始成像时间为T的成像区域了；记录落在各个条带中的离散点坐标 和 

对应的三个整数 

分别按 

和 的大小排序，统计总点数k，就可以得到偏移距为2h、方位角为 地震道的三维时变偏移孔径。

对偏移孔径边缘区的偏移幅值进行平滑的衰减，使其与偏移孔径以外的零值不产生突变，可有效压制偏移噪音。为实现平滑的衰减，可在深度T方向和椭圆边界的法线方向联合进行系数由1到0.2的随正弦函数变化的衰减。

沿T方向衰减得到系数为

其中m₁是事 先定义的沿垂直方向的衰减点数；沿椭圆边界的法线方向衰减得到系数为

$c_2(k,j,n_T^i)=0.8\sin (\frac{1/r^2-f}{1/r^2-1}\·\frac{\mathrm{\π}}{2})+\mathrm{0.2,1}/r^2\≥f(k,j)>1$ 和 $n_T^i=T_0/\mathrm{\ΔT},T_2/\mathrm{\ΔT}$

其中 

和 

是前面得到的椭圆长、短半轴，m₂是事先定义的衰减点数，T₀和T₂分别是偏移计算的起始和最大时间深度。若偏移孔径边缘区中相同的点同时被赋予衰减系数c₁和c₂，则最终的衰减系数定义为c＝(c₁·c₂)^e，其中e＝0.5+(c₁-0.2)/1.6；其它点的衰减系数c或取c₁或取c₂。

对得到的衰减系数作空间高斯平滑，平滑时边界点的值可作为1.0。平滑时计算各空间点的拉普拉斯算子值，即计算无量纲系数

φ＝c(k+1，j，i)+c(k-1，j，i)+c(k，j+1，i)

+c(k，j-1，i)+c(k，j，i+1)+c(k，j，i-1)-6c(k，j，i)

若φ大于指定值，则对该点重新平滑。

时变的三维偏移孔径和对应的衰减系数是事先计算好，存放在指定的表中，计算过程中直接到表中拾取相应的孔径参数和衰减系数。

所述的采用基于单程波算子和稳相点原理的查表法求得各向异性介质中的地震波走时和保幅成像的权系数是这样实现的：定义

$p_0=\frac{{(x-x_0)}^2+{(y-y_0)}^2}{V_{\mathrm{nmo}}^2{T}^2},$ $p_1=\frac{p_0}{p_0+1},\mathrm{\ξ}=\sqrt{\frac{p_0^3+{6p}_0^2+{9p}_0+4-{\mathrm{\η}}_{\mathrm{eff}}({6p}_0^2-12p_0)}{p_0^3+{6p}_0^2+{9p}_0+4+{\mathrm{\η}}_{\mathrm{eff}}({2p}_0^3+{10p}_0^2+44p_0)}}$

式中(x₀，y₀)是地震道的炮点或检波点坐标，(x，y，T)是成像点的水平和时间深度坐标，T是单程旅行时，V_nmo和η_eff是成像点处的偏移速度和各向异性参数，p₀、p₁和ξ是求得的三个无量纲参数。以 

和η_eff为参数建立二维表，表中存两个无量纲数值：

$\widetilde{\mathrm{\τ}}=\sqrt{1-\frac{{\mathrm{\ξ}}^2{p}_1}{1-2{\mathrm{\η}}_{\mathrm{eff}}{\mathrm{\ξ}}^2{p}_1}}+\mathrm{\ξ}\sqrt{p_0{p}_1},$ $A=\frac{{(1-2{\mathrm{\η}}_{\mathrm{eff}}{\mathrm{\ξ}}^2{p}_1)}^2[1-({2\mathrm{\η}}_{\mathrm{eff}}+1){\mathrm{\ξ}}^2{p}_1]}{\sqrt{1+4{\mathrm{\η}}_{\mathrm{eff}}{\mathrm{\ξ}}^2{p}_1-6{\mathrm{\η}}_{\mathrm{eff}}({2\mathrm{\η}}_{\mathrm{eff}}+1){\mathrm{\ξ}}^4{p}_1^2}}$

对每一成像点，由地震道的炮点坐标、成像点坐标和偏移速度计算 

将 

和η_eff按表的间距取整，拾取相邻的4个点的值，由双线性插值得到较准确的 和A_s；再由地震道的检波点坐标、成像点坐标和偏移速度计算 

依据 

和η_eff，同样可从表中得到较准确的 

和A_r；地震波从炮点到成像点、再反射到检波点的走时可简单得到为

$\mathrm{\τ}=({\widetilde{\mathrm{\τ}}}_s+{\widetilde{\mathrm{\τ}}}_r)T$

而计算偏移幅值的保幅成像的无量纲权系数为A_r/A_s。

对地震道的时域信号求一阶导数，通过插值计算拾取时间τ处的值，将该值乘上权系数A_r/A_s，即得到该地震道在该成像点的偏移幅值。插值计算可采用加密重采样的插值技术，以提高插值计算的精度和效率。

本发明的各向异性三维叠前时间偏移方法，能在偏移过程中自主决定偏移速度和各向异性参数，准确获得各向异性介质中地震波的走时和幅值。

本发明的各向异性三维叠前时间偏移方法，通过求得各向异性介质中保幅成像的权系数，能在成像过程中正确补偿地震波传播的几何扩散效应，得到保幅的共反射点道集。

本发明的各向异性三维叠前时间偏移方法，可依据地质构造在沿测线和垂直测线这两个方向的随深度变化的最大倾角，确定考虑介质各向异性的时变的三维偏移孔径，从而可在保证陡倾角构造得到正确成像的前提下，降低偏移噪音和提高偏移方法的计算效率。

本发明的各向异性三维叠前时间偏移方法，通过对偏移孔径边缘区的偏移幅值进行平滑衰减，实现了在成像过程中自动压制偏移噪音， 从而得到更高信噪比的偏移成像结果。

本发明的各向异性三维叠前时间偏移方法，适宜于并行实现，通过在并行计算的各个进程之间合理分配叠前地震资料，实现了高加速比的并行计算。

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

本发明的核心思想有四点：一是从三维叠前深度偏移方法出发分析和研究叠前时间偏移，基于单程波理论和稳相点原理，得到三维各向异性介质中的地震波走时、幅值和等效偏移参数，同时根据叠前深度偏移的反褶积成像条件，得到了保幅成像的权系数；二是从各向异性介质中地震波入射角度和三维空间脉冲响应的倾角分析出发，给出由地质构造随深度变化的最大倾角决定的三维时变偏移孔径，并通过对偏移孔径边缘的偏移幅值进行平滑的衰减来压制偏移噪音；三是结合各向同性偏移和各向异性动校正，获取各向异性偏移的偏移速度和各向异性参数；四是通过对地震资料做偏移孔径相关的等偏移计算量分组，实现在并行计算的各个进程之间最优分配地震资料，获得高加速比的并行计算。其具体实现原理如下：

1.各向异性介质中地震波走时、幅值和权系数

近似三维非均匀介质为层状介质，用横向各向同性近似地球介质的各向异性。由深度偏移的相移法可得，在波数-频率域，单个炮点或检波点的地震波场的深度延拓可用旅行时(时间深度)表示为：

$\tilde{P}(p_x,p_y,\mathrm{\ω},T=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ΔT}}_i)=f\left(\mathrm{\ω}\right)\exp (-\mathrm{j\ω}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ΔT}}_i\sqrt{1-\frac{v_{\mathrm{nmo}}^2\left(i\right)(p_x^2+p_y^2)}{1-{2\mathrm{\η}}_i{v}_{\mathrm{nmo}}^2\left(i\right)(p_x^2+p_y^2)}})---\left(1\right)$

式中ΔT_i是各层介质用单程旅行时表达的厚度，n是某一深度包含的层数， 

是该实际深度的时间深度(单程旅行时)，v_nmo(i)是各层 介质的各向异性动校正速度，η_i是各层介质的各向异性参数，p_x和p_y是射线参数，j是单位虚数，ω是频率，f(ω)是炮点或检波点的时域信号的傅立叶变换。式(1)默认横向各向同性的对称轴是旅行时方向，由于实际偏移实践中旅行时方向并不总是垂直向下，因此式(1)并不限定各向异性介质是垂直对称轴的横向各向同性介质，即不限定为VTI各向异性。

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

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left({\mathrm{\ΔT}}_i\sqrt{1-\frac{v_{\mathrm{nmo}}^2\left(i\right)(p_x^2+p_y^2)}{1-2{\mathrm{\η}}_i{v}_{\mathrm{nmo}}^2\left(i\right)(p_x^2+p_y^2)}}\right)\≈\sqrt{1-\frac{V_{\mathrm{nmo}}^2(p_x^2+p_y^2)}{1-{2\mathrm{\η}}_{\mathrm{eff}}{V}_{\mathrm{nmo}}^2(p_x^2+p_y^2)}}\·\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ΔT}}_i\right)---\left(2\right)$

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

$V_{\mathrm{nmo}}^2=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{v}_{\mathrm{nmo}}^2\left(i\right){\mathrm{\ΔT}}_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ΔT}}_i}---\left(3\right)$

${\mathrm{\η}}_{\mathrm{eff}}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{v}_{\mathrm{nmo}}^4\left(i\right){\mathrm{\η}}_i{\mathrm{\ΔT}}_i}{V_{\mathrm{nmo}}^4\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ΔT}}_i}+\frac{1}{8}\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{v}_{\mathrm{nmo}}^4\left(i\right){\mathrm{\ΔT}}_i}{V_{\mathrm{nmo}}^4\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ΔT}}_i}-\frac{1}{8}---\left(4\right)$

V_nmo和η_eff即是各向异性介质叠前时间偏移所使要的两个等效偏移参数，分别是均方根动校正速度和无量纲的均方根各向异性参数。

将(2)式的近似相移量代入(1)式并作空间付氏反变换，得空间-频率域波场为

$P(x,y,\mathrm{\ω},T)=\frac{\mathrm{\ω}}{4{\mathrm{\π}}^2}\∫\∫f\left(\mathrm{\ω}\right)\exp [-\mathrm{j\ω}(\sqrt{1-\frac{V_{\mathrm{nmo}}^2(p_x^2+p_y^2)}{1-2{\mathrm{\η}}_{\mathrm{eff}}{V}_{\mathrm{nmo}}^2(p_x^2+p_y^2)}}\·T-p_{x}x-p_{y}y)]{\mathrm{dp}}_x{\mathrm{dp}}_y---\left(5\right)$

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

$P(x,y,\mathrm{\ω},T)=f\left(\mathrm{\ω}\right)\frac{\mathrm{\ω}}{2\mathrm{\π}}\exp (-j\frac{\mathrm{\π}}{2}){\left|Q(p_x^0,p_y^0)\right|}^{-1/2}\exp [-\mathrm{j\ω\φ}(p_x^0,p_y^0)]---\left(6\right)$

式(6)中，定义

$\mathrm{\φ}(p_x,p_y)=\sqrt{1-\frac{V_{\mathrm{nmo}}^2(p_x^2+p_y^2)}{1-2{\mathrm{\η}}_{\mathrm{eff}}{V}_{\mathrm{nmo}}^2(p_x^2+p_y^2)}}\·T-p_{x}x-p_{y}y$

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

而 

是 

与 

的零点。求得 

和 

代入(6)式，可解得地震波传播的走时和振幅值。

定义：

$p_0=\frac{x^2+y^2}{V_{\mathrm{nmo}}^2{T}^2},$ $p_1=\frac{p_0}{p_0+1},$ $\mathrm{\ξ}=\sqrt{\frac{p_0^3+{6p}_0^2+{9p}_0+4-{\mathrm{\η}}_{\mathrm{eff}}({6p}_0^2-12p_0)}{p_0^3+{6p}_0^2+{9p}_0+4+{\mathrm{\η}}_{\mathrm{eff}}({2p}_0^3+{10p}_0^2+44p_0)}}$

可解得地震波的走时和振幅值为：

$\mathrm{\τ}=(\sqrt{1-\frac{{\mathrm{\ξ}}^2{p}_1}{1-2{\mathrm{\η}}_{\mathrm{eff}}{\mathrm{\ξ}}^2{p}_1}}+\mathrm{\ξ}\sqrt{p_0{p}_1})T---\left(7\right)$

$A=\frac{{(1-2{\mathrm{\η}}_{\mathrm{eff}}{\mathrm{\ξ}}^2{p}_1)}^2[1-({2\mathrm{\η}}_{\mathrm{eff}}+1){\mathrm{\ξ}}^2{p}_1]}{{\mathrm{TV}}_{\mathrm{nmo}}^2\sqrt{1+4{\mathrm{\η}}_{\mathrm{eff}}{\mathrm{\ξ}}^2{p}_1-6{\mathrm{\η}}_{\mathrm{eff}}({2\mathrm{\η}}_{\mathrm{eff}}+1){\mathrm{\ξ}}^4{p}_1^2}}---\left(8\right)$

尽管(7)和(8)式的走时和振幅值是基于(1)式的层状介质假设推出的，通过允许两个等效偏移参数V_nmo和η_eff横向变化，式(7)和(8)的解实际上代表了三维非均匀各向异性介质中地震波的走时和振幅值。这一点与在波动方程叠前深度偏移中，先由相移法得到单程波方程及其解，然后通过允许速度横向变化来处理介质横向非均匀是一致的。

实际计算中式(7)和(8)的走时和振幅值可应用表的技术简便求得：以 

和η_eff作为两个维构建一个二维表，表中的元素为由 和η_eff计算得到的两个值

和 

实际计算时根据具体的 

和η_eff，通过双线性插值由表中拾取对应的参数值，然后将走时和幅值参数代入(7)和(8)式，即可快速求得准确的走时和振幅值。

既然已基于单程波的相移法和稳相点原理得到了三维非均匀各向异性介质中地震波较准确的走时和振幅值，就可利用波动方程叠前深度偏移的反褶积成像条件保幅成像。设震源是一时间脉冲，则对单个地震道，有成像结果：

$I(x,y,T)=\frac{A_r}{A_s}\∫f\left(\mathrm{\ω}\right)\mathrm{\ωexp}(-j\frac{\mathrm{\π}}{2})\exp (-\mathrm{j\ω}(t_s+t_r))\mathrm{d\ω}$ (9)

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

式中F′(t)是f(ω)对应的时域函数的一阶导数，A_s和t_s是炮点(x_s，y_s，0)至成像点(x，y，T)的幅值和走时，它们可通过先计算 

直接由前述的查表法求得，其中V_nmo是成像点的偏移速度；A_r和t_r是由同样方法求得的检波点至成像点的幅值和走时。实际应用中不必分别计算走时和幅值，直接计算权系数A_r/A_s和总走时即可。

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

不同于忽略成像权系数的现行方法，(9)式的成像权系数保证了偏移成像中正确补偿了地震波在非均匀各向异性介质中传播的几何扩散效应。

式(9)中的F′(t)是离散函数，拾取其t_s+t_r时刻的值需应用插值技术。采用加密重采样的插值技术，可提高插值计算的计算效率和精度。式(9)中使用的时间深度T是单程旅行时，为与现行的以双程旅行时表达的成像结果一致，(9)式的成像结果可用I(x，y，2T)表达。

2.各向异性介质中时变三维偏移孔径

对单个地震道偏移计算而言，(9)式的计算在多大的成像空间，即(x，y，T)范围内进行，就是偏移计算的偏移孔径。偏移孔径对叠前偏移是重要的：较小的孔径可减少偏移计算量，但存在着不能对陡倾角构造正确成像的风险，过大的孔径又带来了偏移噪音和较大的计算量。偏移孔径的选取与地震波的入射角度有关，地球介质各向异性改变了地震波的传播路径，相应地改变了地震波入射角度，因此在确定偏移孔径时需考虑介质的各向异性。

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

以沿方位角方向和垂直方位角方向作为坐标轴看单个地震道的三维偏移脉冲响应，可将三维偏移脉冲响应看作一组由椭圆以方位角方向为轴旋转形成的一组椭球面。本发明根据地质构造的倾角确定偏移孔径，就是仅把倾角小于最大倾角的椭球面部分作为(9)式的计算区域。对不同时间深度T，切除了大于最大倾角部分的三维偏移脉冲响应，可看作以方位角方向为长轴方向、垂直方位角方向为短轴方向、中心在该地震道的中心点上的椭圆面。因此，本发明的偏移孔径是与地震道的偏移距、方位角和地下构造的倾角，以及偏移速度和各向异性参数有关。

在几个关键深度，定义构造沿测线和垂直测线方向的两个最大倾角 和 

对方位角为 

的地震道，测线与成像椭圆面的长轴方向有角 度为 

的夹角，要保证最大倾角为 

和 

的地质构造得以正确成像，需重新计算该深度沿方位角方向和垂直方位角方向的最大成像角度。为此构建如下不等式方程组：

满足上两式的最小角度αⁱ和βⁱ，就是保证该深度最大倾角为 

和 

的地质构造得以正确成像时，沿方位角方向和垂直方位角方向的两个最大成像角度。对αⁱ和βⁱ应用分段线性插值，可得到随时间深度变化的、沿方位角方向和垂直方位角方向的两个最大成像角度α(T)和β(T)。

可根据角度α(T)和β(T)分别计算不同时间深度上成像椭圆的长、短轴。不失一般性，假设椭圆长轴方向为x轴方向；定义入射波与界面法线的夹角为γ，则椭圆长轴顶点处入射和反射地震波的角度为α-γ和α+γ。由(6)式可解得的射线参数为

$p_0^x=\frac{\mathrm{x\ξ}}{V_{\mathrm{nmo}}\sqrt{x^2+V_{\mathrm{nmo}}^2{T}^2}}---\left(12\right)$

式中V_nmo是各向异性的偏移速度，T是时间深度，x是沿长轴方向的坐标，ξ与(7)式中使用的相同，但其中的 

根据射线参数的定义，若忽略偏移速度和层速度的差异，可得：

$\sin (\mathrm{\α}+\mathrm{\γ})=\frac{(a+h){\mathrm{\ξ}}_1}{\sqrt{{\left(V_{\mathrm{nmo}}T\right)}^2+{(a+h)}^2}}---\left(13\right)$

$\sin (\mathrm{\α}-\mathrm{\γ})=\frac{(a-h){\mathrm{\ξ}}_2}{\sqrt{{\left(V_{\mathrm{nmo}}T\right)}^2+{(a-h)}^2}}---\left(14\right)$

式中a是椭圆的长半轴，h是地震道的半偏移距，ξ₁和ξ₂分别是取 

和 

时得到的ξ。直接求解(13)和(14) 式的非线性方程是困难的，由于η_eff是远小于1的正数，可近似

a＝a₀(η_eff＝0)+Δa(η_eff)，Δa/a□1      (15)

式中a₀(η_eff＝0)是忽略各向异性时的长轴值，Δa是考虑各向异性时长轴的增加量。令η_eff＝0，可由(13)和(14)式解得

$a_0=\frac{V_{\mathrm{nmo}}\·T}{\sin \left(2\mathrm{\α}\right)}\sqrt{1+{\cos }^2\left(2\mathrm{\α}\right)+{\mathrm{\λ}}^2{\sin }^2\left(2\mathrm{\α}\right)-2\cos \left(2\mathrm{\α}\right){[1+{\mathrm{\λ}}^2{\sin }^2\left(2\mathrm{\α}\right)]}^{1/2}}---\left(16\right)$

式中λ＝h/(V_nmo·T)是无量纲参数。将(15)和(16)式代入(13)和(14)式中，忽略Δa/a和η_eff的高级项，可解得

$\mathrm{\Δa}=\frac{{5a}_0{\mathrm{\η}}_{\mathrm{eff}}({\tilde{a}}^2-{\mathrm{\λ}}^2-1)}{{2\sin }^2\left(2\mathrm{\α}\right)({\tilde{a}}^2-{\mathrm{\λ}}^2+1)-4}---\left(17\right)$

式中 

是无量纲参数。

令h＝0，可知γ＝0，此时(13)和(14)式变的相同，可直接由该式解得椭圆的短半轴b为

$b=V_{\mathrm{nmo}}\·T\sin \mathrm{\β}{[{\cos }^2\mathrm{\β}-{\mathrm{\η}}_{\mathrm{eff}}\left(\frac{{\sin }^2\mathrm{\β}{(1+3{\cos }^2\mathrm{\β})}^2}{{\sin }^2\mathrm{\β}{(1+2{\cos }^2\mathrm{\β})}^2+4{\cos }^6\mathrm{\β}}\right)]}^{-1/2}---\left(18\right)$

实际应用中，为简化计算，将对有相同偏移距和方位角的地震道使用相同的偏移孔径，因此(16)-(18)式中的V_nmo和η_eff将使用不同时间深度上的平均值；这避免了偏移孔径的横向变化。对地下构造的倾角、偏移速度和各向异性参数存在明显变化的偏移工区，可采用分区域方法处理。

由于有限采集缆长、有效偏移孔径等限制因素，偏移的累加计算将产生偏移噪音。若在单个地震道的偏移幅值累加前引入沿偏移孔径边缘的衰减，可有效实现偏移噪音压制。单个地震道的成像空间与偏移孔径对应，为实现平滑的衰减，可在深度T和椭圆边界的法线方向联合进行系数由1到0.2的随正弦函数变化的衰减。

依据从偏移孔径得到的(x_i，y_i)点的起始成像时间T_b，可得沿T方向衰减的系数：

$c_1=0.8\sin (\frac{T-T_b}{\mathrm{m\&Delta;T}}\&CenterDot;\frac{\mathrm{\&pi;}}{2})+0.2,(T<T_b+\mathrm{m\&Delta;T})---\left(19\right)$

式中m是定义的垂直方向衰减点数，ΔT是时间深度(单程旅行时)方向的采样。

沿椭圆边界的法线方向的衰减可根据由下式计算得到的椭圆方程值f进行，

式中 是地震道的方位角。f＝1表明该点落在椭圆边界上，f越小该点就更靠近中心。令沿椭圆短轴方向缩进m个点时衰减系数为1(即内部不衰减)，计算 

可得沿椭圆边界的法线方向衰减的系数为

$c_2=0.8\sin (\frac{f/r^2-1/r^2}{1-1/r^2}\&CenterDot;\frac{\mathrm{\&pi;}}{2})+0.2---\left(21\right)$

若偏移孔径边缘区中相同的点同时有衰减系数c₁和c₂，为避免双重衰减，可定义系数为(c₁·c₂)^e，其中e＝0.5+(c₁-0.2)/1.6。

为避免衰减系数的突变引入的新噪音，需进一步对得到的衰减系数作空间高斯平滑。当遇到不衰减的区域时，可令对应的系数为1.0。平滑时计算各空间点的拉普拉斯算子值，若该值大于指定值，则对该点重新平滑。

3.各向异性介质偏移参数建模

式(3)和(4)表明，各向异性叠前时间偏移利用V_nmo和η_eff这两个偏移参数，它们分别是均方根动校正速度和均方根各向异性参数。由于各向异性参数η_eff导致的走时变化较偏移速度V_nmo的影响是小量，我们 发展了结合各向同性偏移和各向异性动校正的各向异性偏移参数建模方法。

首先选取部分共中心点，对这些点的共中心点道集作常规的动校正(NMO)速度拾取，对所得结果做横向平均，可作为各向同性偏移的初始偏移速度。基于这一速度做各向同性偏移，生成对共反射点道集，该道集较常规的CMP道集更接近真正的“共反射点”道集。利用横向均匀的初始偏移速度对共反射点道集做反动校，再做动校正可得到新的偏移速度，对这一速度做空间平滑处理，可作为各向异性偏移的初始速度V_nmo。

地球介质各向异性或使得长偏移距地震道不能被拉平，或使的将同相轴拉平的偏移速度场存在突变。因此，当各向异性存在时，平滑的V_nmo不能将反动校后的共反射点道集拉平，而这一剩余量恰好可用于拾取各向异性参数η_eff。对反动校后的道集，利用式(22)做各向异性动校正

${\mathrm{\&tau;}}^2=\frac{3+4\mathrm{\&eta;}}{4(1+\mathrm{\&eta;})}H\left(2h\right)+\frac{1}{4(1+\mathrm{\&eta;})}\sqrt{H^2\left(2h\right)+16\mathrm{\&eta;}(1+\mathrm{\&eta;})\frac{{\left(2T\right)}^2{\left(2h\right)}^2}{(1+2\mathrm{\&eta;})V_{\mathrm{nmo}}^2}}$ (22)

$H\left(2h\right)={\left(2T\right)}^2+\frac{{\left(2h\right)}^2}{(1+2\mathrm{\&eta;})V_{\mathrm{nmo}}^2}$

式中h是地震道的半偏移距，T是用单程旅行时表达的时间深度，V_nmo是已知的偏移速度，η是各向异性参数。根据叠加能量最大和η的空间联续性(η的横向变化不应过大)，可确定各向异性参数η，它即可作为η_eff。

为避免V_nmo的误差影响η_eff的拾取，再分别取速度为0.95V_nmo和1.05V_nmo进行式(22)各向异性动校正；叠加能量最大、同相轴最平直的速度和各向异性参数即是最终的各向异性偏移的偏移速度和各向异性参数，对各向异性参数还需进行空间平滑处理。

4.偏移孔径相关的等计算量分组

三维叠前时间偏移一般是利用集群计算机，由多个计算节点并行计算来完成的。并行计算是通过将叠前地震资料分组，由每个计算节点独立处理一组地震资料来进行的。因此，叠前地震资料分组是否恰当，能否使各个计算节点同时完成计算任务是提高并行计算效率的关键。叠前时间偏移计算中相同偏移距数据的计算结果可相互叠加，而不同偏移距数据的偏移结果各自独立，因此共偏移距数据是叠前时间偏移数据分组的基础。但同样的数据量，不同偏移距地震资料的偏移计算量是不同的，因此合理的分组既要考虑数据量，也要考虑地震资料的偏移距。

导致不同偏移距地震资料的偏移计算量不同的关键是偏移孔径，偏移距越大偏移孔径也越大，相应的偏移计算量也就越大。为简化计算，在求取估算偏移计算量使用的偏移孔径时，忽略各向异性和方位角的影响，仅在一个深度定义最大成像角度(但引入时变偏移孔径的临界时间深度，即大于此值偏移孔径将不随时间深度变化)；这样，可由(16)和(18)式直接得到不同时间深度上偏移孔径对应的成像椭圆的长、短半轴为

$a\left(T\right)=\frac{{\stackrel{\&OverBar;}{V}}_{\mathrm{rms}}\left(T\right)\&CenterDot;T}{\sin \left(2{\mathrm{\&theta;}}_x\right)}\sqrt{1+{\cos }^2\left({2\mathrm{\&theta;}}_x\right)+{\mathrm{\&lambda;}}^2{\sin }^2\left({2\mathrm{\&theta;}}_x\right)-2\cos \left({2\mathrm{\&theta;}}_x\right){[1+{\mathrm{\&lambda;}}^2{\sin }^2\left({2\mathrm{\&theta;}}_x\right)]}^{1/2}}---\left(23\right)$

$b\left(T\right)={\stackrel{\&OverBar;}{V}}_{\mathrm{rms}}\left(T\right)\&CenterDot;T{\tan \mathrm{\&theta;}}_y---\left(24\right)$

式中θ_x和θ_y是地下构造在沿测线和垂直测线方向的最大倾角， 是初始的、横向均匀的偏移速度，h是地震道的半偏移距， 

是无量纲参数。

令偏移计算的起始时间深度(单程旅行时)为T₀，最大时间深度为T₂，时变偏移孔径的临界时间深度是T₁，则偏移距为2h的地震道的相对偏移计算量可近似为

$G\left(h\right)=a\left(T_1\right)\&CenterDot;b\left(T_1\right)(T_2-T_1)(1/\mathrm{\&Delta;T})+\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}a(T_0+\mathrm{j\&Delta;T})\&CenterDot;b(T_0+\mathrm{j\&Delta;T})---\left(25\right)$

式中n＝(T₁-T₀)/ΔT是一整数，ΔT是时间深度方向上成像的采样间距(单程旅行时)。

根据G(h)，可得到不同偏移距地震资料的与计算量有关的权系数，将权系数与数据量相乘，得到的新的数据量就同时考虑了偏移距和数据量对偏移计算量的影响；对这一新数据量按等数据量分组，就可以实现各个计算节点同时完成计算任务。

本发明的有益效果：该方法可应用于多条拖缆或多条测线记录的三维地震资料，生成具有较高信噪比、能正确刻画地下三维陡倾角构造的偏移图像。偏移图像能指示地下构造的形态、断裂部位和地层沉积样式，对确定有利生、储油构造和识别有利油气储层有重要应用价值。该方法可正确处理地球介质各向异性对地震波传播的走时和幅值的影响，能在偏移过程中决定偏移速度和各向异性参数，因而可得到归位更准确且保幅的偏移图像；它对陡倾角构造有更好的成像效果。该方法生成的保幅共反射点(CRP)道集可更好地服务于叠前反演等油气和流体检测技术。

附图说明

图1是渤海区域某区块常规噪音压制处理后的典型单炮记录，采集是每炮3条拖缆。

图2是偏移速度在沿测线切面上的等值线图。图中数字是偏移速度值。

图3是各向异性参数在沿测线切面上的等值线图。图中数字是各向异性参数值。

图4是偏移距为1845m、方位角为30度地震道的偏移孔径在沿测 线和垂直测线切面上的起始成像时间(双程旅行时)曲线。

图5是图4偏移孔径的孔径边缘区局部的衰减系数等值线图。图中数字是衰减系数值。

图6是渤海区域某区块地下潜山构造偏移结果沿测线方向的剖面图像。

图7是各计算节点上地震资料的偏移距范围和地震道道数表。

具体实施方式

实施例1：各向异性三维叠前时间偏移方法，针对渤海区域某区块为例，具体为以下步骤：

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

(2)从磁带上读取地震信号，对叠前地震资料进行常规的压制噪音处理。针对部分共中心点，抽取共中心点道集，对抽取的道集作常规的NMO动校正速度拾取，对所得结果做横向平均，作为初始偏移速度，记为 

(3)将叠前地震资料按偏移距大小排序，对地震资料做偏移孔径相关的相等计算量分组，将不同组地震资料存放到集群计算机的不同计算节点上。所使用的PC集群计算机共有35个节点，70个核，表1给出了各计算节点上地震资料的偏移距范围和地震道道数。

设地下构造在沿测线和垂直测线方向的两个最大倾角为θ_x和θ_y，已得到的横向均匀的初始偏移速度为 

计算两个长度值

$a\left(T\right)=\frac{{\stackrel{\&OverBar;}{V}}_{\mathrm{rms}}\left(T\right)\&CenterDot;T}{\sin \left(2{\mathrm{\&theta;}}_x\right)}\sqrt{1+{\cos }^2\left({2\mathrm{\&theta;}}_x\right)+{\mathrm{\&lambda;}}^2{\sin }^2\left({2\mathrm{\&theta;}}_x\right)-2\cos \left({2\mathrm{\&theta;}}_x\right){[1+{\mathrm{\&lambda;}}^2{\sin }^2\left({2\mathrm{\&theta;}}_x\right)]}^{1/2}}$

$b\left(T\right)={\stackrel{\&OverBar;}{V}}_{\mathrm{rms}}\left(T\right)\&CenterDot;T{\tan \mathrm{\&theta;}}_y$

式中T是用单程旅行时表达的时间深度，单位是秒，h是地震道的半偏移距，单位是米， 

是无量纲参数。时间深度T处的成像区域可看作以a(T)和b(T)为长、短半轴的椭圆面。令偏移计算的起始时间深度为T₀，最大时间深度为T₂，而大于T₁偏移孔径将不随时间深度变化，则偏移距为2h的地震道的相对偏移计算量可近似为

$G\left(h\right)=a\left(T_1\right)\&CenterDot;b\left(T_1\right)(T_2-T_1)(1/\mathrm{\&Delta;T})+\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}a(T_0+\mathrm{j\&Delta;T})\&CenterDot;b(T_0+\mathrm{j\&Delta;T})$

式中n＝(T₁-T₀)/ΔT是一整数，ΔT是时间深度方向上成像的采样间距，a(T₁)和b(T₁)等可由前一个公式计算得到。设地震资料的最大半偏移距为h_max，最小半偏移距为h_min，偏移距采样间距为Δh，可用整数i＝(h-h_min)/Δh+1作为偏移距的索引；统计有相同i的地震道，可得到偏移距满足h_min+(i-1)Δh≤h＜h_min+iΔh的地震道道数m_i(i＝1，l)，其中l＝(h_max-h_min)/Δh+1。记所有m_i中最大的为m_j，可得无量纲的权系数ρ_i＝G(h_min+(i-1)Δh)/G(h_min+(j-1)Δh)。

如集群计算机共有k个核，即希望同时有k个进程并行计算，则计算 

可按如下方式确定各进程中包含的地震资料的偏移距索引：对进程1，求满足下式的n₁，

$\underset{i=1}{\overset{n_1}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{m}_i\&le;m_0\&le;\underset{i=1}{\overset{n_1+1}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{m}_i$

若

$m_0-\underset{i=1}{\overset{n_1}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{m}_i>\underset{i=1}{\overset{n_1+1}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{m}_i-m_0$

取n₁＝n₁+1，索引数为i＝1，n₁的地震道将分配给进程1；对进程2，求满足下式的n₂，

$\underset{i=n_1+1}{\overset{n_2}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{m}_i\&le;m_0\&le;\underset{i=n_1+1}{\overset{n_2+1}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{m}_i$

若

$m_0-\underset{i=n_1+1}{\overset{n_2}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{m}_i>\underset{i=n_1+1}{\overset{n_2+1}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{m}_i-m_0$

取n₂＝n₂+1，索引数为i＝n₁+1，n₂的地震道将分配给进程2；如此类推，完成全部k个进程的地震资料分配，这样可使得各进程的偏移计算量近似相同。一般集群计算机的一个计算节点可同时进行几个进程的计算，可将相关进程的地震资料一起存放到这个计算节点上。

(4)利用初始偏移速度，对已存放到集群计算机各个计算节点上的地震资料，应用时变偏移孔径的三维各向同性叠前时间偏移方法，进行并行的各向同性偏移计算，收集各计算节点的偏移结果，形成共反射点道集。

(5)依据形成的共反射点道集，确定各向异性偏移使用的偏移速度和各向异性参数。图2是得到的偏移速度在沿测线切面上的等值线图，图3是得到的各向异性参数在沿测线切面上的等值线图。

对共反射点道集，利用初始偏移速度做反动校，再做动校正得到新的速度，对这一速度做空间平滑处理，作为各向异性偏移的初始速度V_nmo。对反动校后的共反射点道集，再次利用下式进行各向异性动校正

${\mathrm{\&tau;}}^2=\frac{3+4\mathrm{\&eta;}}{4(1+\mathrm{\&eta;})}H\left(2h\right)+\frac{1}{4(1+\mathrm{\&eta;})}\sqrt{H^2\left(2h\right)+16\mathrm{\&eta;}(1+\mathrm{\&eta;}\frac{{\left(2T\right)}^2{\left(2h\right)}^2}{(1+2\mathrm{\&eta;})V_{\mathrm{nmo}}^2})}$

$H\left(2h\right)={\left(2T\right)}^2+\frac{{\left(2h\right)}^2}{(1+2\mathrm{\&eta;})V_{\mathrm{nmo}}^2}$

式中h是地震道的半偏移距，单位是米，T是用单程旅行时表达的时间深度，单位是秒，η是无量纲的各向异性参数。可根据叠加能量最大和η的空间联续性，即η的横向变化不应过大，确定各空间点的各向异性参数η；再分别取速度为0.95V_nmo和1.05V_nmo进行上式各向异性动校正，(0.95V_nmo，η)、(V_nmo，η)和(1.05V_nmo，η)三组参数中使叠加能量最大、同相轴最平直的那组即是该点最终的偏移速度和各向异性参数，对全部各向异性参数还需进行空间平滑处理。

(6)根据地下构造在不同时间深度的最大倾角、偏移速度、各向异性参数和地震道的偏移距和方位角，确定时变的三维偏移孔径和对应于偏移孔径的、用于压制偏移噪音的衰减系数。定义不同时间深度处构造沿测线和垂直测线方向的最大倾角为：双程旅行时1s处25度和10度，2s处45度和20度，2.2s处以上0度和0度；用双程旅行时表达的成像起始时间T₀＝0.4s。图4是偏移距为1845m、方位角为30度地震道的偏移孔径在沿测线和垂直测线切面上的双程旅行时起始成像时间曲线，图5是偏移孔径边缘区局部的压制偏移噪音的衰减系数的等值线图。

用一组整数 

描述三维偏移孔径，其中k是有效成像区域中总的离散点数。令Δx和Δy是偏移成像结果在两个水平坐标方向的空间采样，ΔT是时间深度方向的采样，则 

和 是成像区域中某离散点与地震道中心点的沿两个水平坐标方向的距离，而 

是该点偏移计算的起始成像时间；令地震道中心点坐标为 

和 

引入时变的偏移孔径，就是在该地震道的偏移计算中，对成像区域中水平坐标为 

和 

的点(i＝1，k)，仅从 

开始进行偏移计算。

在几个关键深度，定义地质构造在沿测线和垂直测线方向的最大 倾角为 

和 构建如下不等式方程组：

式中 是地震道的方位角。满足上两式的最小角度αⁱ和βⁱ，就是保证该深度最大倾角为 

和 的地质构造得以正确成像时，在沿方位角方向和垂直方位角方向的两个最大成像角度。

对αⁱ和βⁱ应用分段线性插值，可得到随时间深度变化的、沿方位角和垂直方位角方向的两个最大成像角度α(T)和β(T)。对不同时间深度T，偏移孔径所定义的成像区域，可看作以方位角方向为长轴方向、垂直方位角方向为短轴方向、中心在地震道的中心点上的椭圆面，其长、短半轴a(T)和b(T)分别为：

a(T)＝a₀+Δa

$a_0=\frac{{\stackrel{\&OverBar;}{V}}_{\mathrm{nmo}}\left(T\right)\&CenterDot;T}{\sin \left[2\mathrm{\&alpha;}\left(T\right)\right]}\sqrt{1+{\cos }^2\left[2\mathrm{\&alpha;}\left(T\right)\right]+{\mathrm{\&lambda;}}^2{\sin }^2\left[2\mathrm{\&alpha;}\left(T\right)\right]-2\cos \left[2\mathrm{\&alpha;}\left(T\right)\right]{\{1+{\mathrm{\&lambda;}}^2{\sin }^2[2\mathrm{\&alpha;}\left(T\right)\left]\right\}}^{1/2}}$

$\mathrm{\&Delta;a}=\frac{{5a}_0{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_{\mathrm{eff}}\left(T\right)({\tilde{a}}^2-{\mathrm{\&lambda;}}^2-1)}{2{\sin }^2\left[2\mathrm{\&alpha;}\left(T\right)\right]({\tilde{a}}^2-{\mathrm{\&lambda;}}^2+1)-4}$

$b\left(T\right)={\stackrel{\&OverBar;}{V}}_{\mathrm{nmo}}\left(T\right)\&CenterDot;T\sin \left[\mathrm{\&beta;}\left(T\right)\right]{\left({\cos }^2\right[\mathrm{\&beta;}\left(T\right)]-\frac{{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_{\mathrm{eff}}\left(T\right){\sin }^2\left[\mathrm{\&beta;}\left(T\right)\right]{\{1+{3\cos }^2[\mathrm{\&beta;}\left(T\right)\left]\right\}}^2}{{\sin }^2\left[\mathrm{\&beta;}\left(T\right)\right]{\{1+2{\cos }^2[\mathrm{\&beta;}\left(T\right)\left]\right\}}^2+{4\cos }^6\left[\mathrm{\&beta;}\left(T\right)\right]})}^{-1/2}$

式中h是地震道的半偏移距，T是用单程旅行时表达的时间深度，a₀是忽略各向异性时的长轴值，Δa是考虑各向异性时长轴的增加量， 

是无量纲参数， 

是无量纲参数， 

是不同时间深度上的平均各向异性参数，无量纲， 

是不同时间深度上的平均偏移速度。

设偏移计算的起始时间深度为T₀，最大时间深度为T₂，对T＝T₀，T₂循 环，若a(T)≥a(T-ΔT)，a(T)不变，否则a(T)＝a(T-ΔT)；对a(T)作时间深度方向的高斯平滑，即得随时间深度变化的长半轴，同理可得随时间深度变化的短半轴。

对T＝T₀，T₂循环，对时间深度T和T+ΔT，由其对应的两组长、短半轴形成有共同的中心点(0，0)的两个椭圆，这两个椭圆形成的条带就是起始成像时间为T的成像区域了；记录落在各个条带中的离散点坐标 

和 

对应的三个整数 

分别按 和 的大小排序，统计总点数k，就可以得到偏移距为2h、方位角为 

地震道的三维时变偏移孔径。

对偏移孔径边缘区的偏移幅值进行平滑的衰减，使其与偏移孔径以外的零值不产生突变，可有效压制偏移噪音。为实现平滑的衰减，可在深度T方向和椭圆边界的法线方向联合进行系数由1到0.2的随正弦函数变化的衰减。

沿T方向衰减得到系数为

(j＝1，m₁和i＝1，k)，其中m₁是事先定义的沿垂直方向的衰减点数；沿椭圆边界的法线方向衰减得到系数为

$c_2(k,j,n_T^i)=0.8\sin (\frac{1/r^2-f}{1/r^2-1}\&CenterDot;\frac{\mathrm{\&pi;}}{2})+\mathrm{0.2,1}/r^2\&GreaterEqual;f(k,j)>1$ 和 $n_T^i=T_0/\mathrm{\&Delta;T},T_2/\mathrm{\&Delta;T}$

其中 

和 是前面得到的椭圆长、短半轴，m₂是事先定义的衰减点数，T₀和T₂分别是偏移计算的起始和最大时间深度。若偏移孔径边缘区中相同的点同时被赋予衰减系数c₁和c₂，则最终的衰减系数定义为c＝(c₁·c₂)^e，其中e＝0.5+(c₁-0.2)/1.6；其它点的衰减系数c或取c₁或取c₂。

对得到的衰减系数作空间高斯平滑，平滑时边界点的值可作为1.0。平滑时计算各空间点的拉普拉斯算子值，即计算无量纲系数

φ＝c(k+1，j，i)+c(k-1，j，i)+c(k，j+1，i)

+c(k，j-1，i)+c(k，j，i+1)+c(k，j，i-1)-6c(k，j，i)

若φ大于指定值，则对该点重新平滑。

时变的三维偏移孔径和对应的衰减系数是事先计算好，存放在指定的表中，计算过程中直接到表中拾取相应的孔径参数和衰减系数。

(7)采用基于单程波算子和稳相点原理的查表法求得各向异性介质中的地震波走时和保幅成像的权系数。

定义

$p_0=\frac{{(x-x_0)}^2+{(y-y_0)}^2}{V_{\mathrm{nmo}}^2{T}^2},$ $p_1=\frac{p_0}{p_0+1},$ $\mathrm{\&xi;}=\sqrt{\frac{p_0^3+{6p}_0^2+{9p}_0+4-{\mathrm{\&eta;}}_{\mathrm{eff}}({6p}_0^2-12p_0)}{p_0^3+{6p}_0^2+{9p}_0+4+{\mathrm{\&eta;}}_{\mathrm{eff}}({2p}_0^3+{10p}_0^2+44p_0)}}$

式中(x₀，y₀)是地震道的炮点或检波点坐标，(x，y，T)是成像点的水平和时间深度坐标，T是单程旅行时，V_nmo和η_eff是成像点处的偏移速度和各向异性参数，p₀、p₁和ξ是求得的三个无量纲参数。以 

和η_eff为参数建立二维表，表中存两个无量纲数值：

$\widetilde{\mathrm{\&tau;}}=\sqrt{1-\frac{{\mathrm{\&xi;}}^2{p}_1}{1-2{\mathrm{\&eta;}}_{\mathrm{eff}}{\mathrm{\&xi;}}^2{p}_1}}+\mathrm{\&xi;}\sqrt{p_0{p}_1},$ $A=\frac{{(1-2{\mathrm{\&eta;}}_{\mathrm{eff}}{\mathrm{\&xi;}}^2{p}_1)}^2[1-({2\mathrm{\&eta;}}_{\mathrm{eff}}+1){\mathrm{\&xi;}}^2{p}_1]}{\sqrt{1+4{\mathrm{\&eta;}}_{\mathrm{eff}}{\mathrm{\&xi;}}^2{p}_1-6{\mathrm{\&eta;}}_{\mathrm{eff}}({2\mathrm{\&eta;}}_{\mathrm{eff}}+1){\mathrm{\&xi;}}^4{p}_1^2}}$

对每一成像点，由地震道的炮点坐标、成像点坐标和偏移速度计算 

将 

和η_eff按表的间距取整，拾取相邻的4个点的值，由双线性插值得到较准确的 

和A_s；再由地震道的检波点坐标、成像点坐标和偏移速度计算 

依据 

和η_eff，同样可从表中得到较准确的 

和A_r；地震波从炮点到成像点、再反射到检波点的走时可简单得到为

$\mathrm{\&tau;}=({\widetilde{\mathrm{\&tau;}}}_s+{\widetilde{\mathrm{\&tau;}}}_r)T$

而计算偏移幅值的保幅成像的无量纲权系数为A_r/A_s。

对地震道的时域信号求一阶导数，通过插值计算拾取时间τ处的值，将该值乘上权系数A_r/A_s，即得到该地震道在该成像点的偏移幅值。插值计算可采用加密重采样的插值技术，以提高插值计算的精度和效率。

(8)在每个计算节点，对地震资料的所有地震道循环。对每一地震道，由偏移孔径确定有效的成像区域和区域中各成像点偏移计算的起始成像时间，利用保幅成像的权系数，从起始成像时间开始计算各成像点的偏移幅值。对偏移孔径边缘区域的成像点，进一步用衰减系数来衰减偏移幅值。将偏移幅值累加到存放偏移结果的数组中对应的偏移距上。具体是，若三维成像空间大小为(n_x，n_y，n_T)，则偏移结果将依据每一地震道的半偏移距h，累加到该数组第4维的m₀＝(h-h_min)/Δh+1上，其中h_min为地震资料的最小半偏移距，这里定义Δh＝12.5m。

(9)收集各计算节点的偏移结果，对形成的全部共反射点道集做常规的各向同性剩余动校正，将出现明显拉伸和噪音部分的对应数值置为零，将不同偏移距的偏移结果叠加，形成偏移叠加剖面。

(10)通过显示软件将偏移叠加剖面数值转换为地下反射构造的剖面图像，剖面图像将指示地下构造的形态、断裂部位、断距大小和地层沉积样式及地层的波阻抗特征，用于确定地下生、储油构造和识别油气储层。图6是用本发明方法得到的渤海区域某区块地下潜山构造沿测线方向的剖面图像，潜山侧翼、顶部、断层和地层沉积样式得到很好地刻画。

Claims (5)

1.一种各向异性三维叠前时间偏移方法，其特征在于：采用以下步骤：A)用拖缆或测线记录人工震源激发的反射地震信号，记录到磁带上；B)从磁带上读取地震信号，对叠前地震资料进行常规的压制噪音处理，针对部分共中心点，抽取共中心点道集，对抽取的道集作常规的NMO动校正速度拾取，对所得结果做横向平均，作为初始偏移速度；C)将叠前地震资料按偏移距大小排序，对地震资料做偏移孔径相关的相等计算量分组，将不同组地震资料存放到集群计算机的不同计算节点上；D)利用初始偏移速度，对已存放到集群计算机各个计算节点上的地震资料，应用时变偏移孔径的三维各向同性叠前时间偏移方法，进行并行的各向同性偏移计算，收集各计算节点的偏移结果，形成共反射点道集；E)依据形成的共反射点道集，确定各向异性偏移使用的偏移速度和各向异性参数；F)根据地下构造在不同时间深度的最大倾角、偏移速度、各向异性参数和地震道的偏移距和方位角，确定时变的三维偏移孔径和对应于偏移孔径的、用于压制偏移噪音的衰减系数；G)采用基于单程波算子和稳相点原理的查表法求得各向异性介质中的地震波走时和保幅成像的权系数；H)在每个计算节点，对地震资料的所有地震道循环，对每一地震道，由偏移孔径确定有效的成像区域和区域中各成像点偏移计算的起始成像时间，利用保幅成像的权系数，从起始成像时间开始计算各成像点的偏移幅值，对偏移孔径边缘区域的成像点，进一步用衰减系数来衰减偏移幅值，将偏移幅值累加到存放偏移结果的数组中对应的偏移距上；I)收集各计算节点的偏移结果，对形成的全部共反射点道集做常规的各向同性剩余动校正，将出现明显拉伸和噪音部分的对应数值置为零，将不同偏移距的偏移结果叠加，形成偏移叠加剖面；J)通过显示软件将偏移叠加剖面数值转换为地下反射构造的剖面图像，剖面图像将指示地下构造的形态、断裂部位、断距大小和地层沉积样式及地层的波阻抗特征，用于确定地下生、储油构造和识别油气储层。

2.根据权力要求1所述的一种各向异性三维叠前时间偏移方法，其特征在于：在C步骤中，所述的将叠前地震资料按偏移距大小排序，对地震资料做偏移孔径相关的相等计算量分组，将不同组地震资料存放到集群计算机的不同计算节点上是这样实现的：设地下构造在沿测线和垂直测线方向的两个最大倾角为θ_x和θ_y，已得到的横向均匀的初始偏移速度为计算两个长度值

$a\left(T\right)=\frac{{\stackrel{\&OverBar;}{V}}_{\mathrm{rms}}\left(T\right)\&CenterDot;T}{\sin \left(2{\mathrm{\&theta;}}_x\right)}\sqrt{1+{\cos }^2\left({2\mathrm{\&theta;}}_x\right)+{\mathrm{\&lambda;}}^2{\sin }^2\left({2\mathrm{\&theta;}}_x\right)-2\cos \left({2\mathrm{\&theta;}}_x\right){[1+{\mathrm{\&lambda;}}^2{\sin }^2\left({2\mathrm{\&theta;}}_x\right)]}^{1/2}}$

$b\left(T\right)={\stackrel{\&OverBar;}{V}}_{\mathrm{rms}}\left(T\right)\&CenterDot;T{\tan \mathrm{\&theta;}}_y$

式中T是用单程旅行时表达的时间深度，单位是秒，h是地震道的半偏移距，单位是米，

是无量纲参数，时间深度T处的成像区域可看作以a(T)和b(T)为长、短半轴的椭圆面，令偏移计算的起始时间深度为T₀，最大时间深度为T₂，而大于T₁偏移孔径将不随时间深度变化，则偏移距为2h的地震道的相对偏移计算量可近似为

$G\left(h\right)=a\left(T_1\right)\&CenterDot;b\left(T_1\right)(T_2-T_1)(1/\mathrm{\&Delta;T})+\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}a(T_0+\mathrm{j\&Delta;T})\&CenterDot;b(T_0+\mathrm{j\&Delta;T})$

式中n＝(T₁-T₀)/ΔT是一整数，ΔT是时间深度方向上成像的采样间距，a(T₁)和b(T₁)等可由前一个公式计算得到，设地震资料的最大半偏移距为h_max，最小半偏移距为h_min，偏移距采样间距为Δh，可用整数i＝(h-h_min)/Δh+1作为偏移距的索引；统计有相同i的地震道，可得到偏移距满足h_min+(i-1)Δh≤h＜h_min+iΔh的地震道道数m_i(i＝1，l)，其中l＝(h_max-h_min)/Δh+1，记所有m_i中最大的为m_j，可得无量纲的权系数ρ_i＝G(h_min+(i-1)Δh)/G(h_min+(j-1)Δh)；

如集群计算机共有k个核，即希望同时有k个进程并行计算，则计算

可按如下方式确定各进程中包含的地震资料的偏移距索引：对进程1，求满足下式的n₁，

$\underset{i=1}{\overset{n_1}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{m}_i\&le;m_0\&le;\underset{i=1}{\overset{n_1+1}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{m}_i$

若

$m_0-\underset{i=1}{\overset{n_1}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{m}_i>\underset{i=1}{\overset{n_1+1}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{m}_i-m_0$

取n₁＝n₁+1，索引数为i＝1，n₁的地震道将分配给进程1；对进程2，求满足下式的n₂，

$\underset{i=n_1+1}{\overset{n_2}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{m}_i\&le;m_0\&le;\underset{i=n_1+1}{\overset{n_2+1}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{m}_i$

若

$m_0-\underset{i=n_1+1}{\overset{n_2}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{m}_i>\underset{i=n_1+1}{\overset{n_2+1}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{m}_i-m_0$

取n₂＝n₂+1，索引数为i＝n₁+1，n₂的地震道将分配给进程2；如此类推，完成全部k个进程的地震资料分配，这样可使得各进程的偏移计算量近似相同；一般集群计算机的一个计算节点可同时进行几个进程的计算，可将相关进程的地震资料一起存放到这个计算节点上。

3.根据权力要求1所述的一种各向异性三维叠前时间偏移方法，其特征在于：在E步骤中，所述的依据形成的共反射点道集，确定各向异性偏移使用的偏移速度和各向异性参数是这样实现的：对共反射点道集，利用初始偏移速度做反动校，再做动校正得到新的速度，对这一速度做空间平滑处理，作为各向异性偏移的初始速度V_nmo，对反动校后的共反射点道集，再次利用下式进行各向异性动校正

${\mathrm{\&tau;}}^2=\frac{3+4\mathrm{\&eta;}}{4(1+\mathrm{\&eta;})}H\left(2h\right)+\frac{1}{4(1+\mathrm{\&eta;})}\sqrt{H^2\left(2h\right)+16\mathrm{\&eta;}(1+\mathrm{\&eta;}\frac{{\left(2T\right)}^2{\left(2h\right)}^2}{(1+2\mathrm{\&eta;})V_{\mathrm{nmo}}^2})}$

$H\left(2h\right)={\left(2T\right)}^2+\frac{{\left(2h\right)}^2}{(1+2\mathrm{\&eta;})V_{\mathrm{nmo}}^2}$

式中h是地震道的半偏移距，单位是米，T是用单程旅行时表达的时间深度，单位是秒，η是无量纲的各向异性参数，可根据叠加能量最大和η的空间联续性，即η的横向变化不应过大，确定各空间点的各向异性参数η；再分别取速度为0.95V_nmo和1.05V_nmo进行上式各向异性动校正，(0.95V_nmo，η)、(V_nmo，η)和(1.05V_nmo，η)三组参数中使叠加能量最大、同相轴最平直的那组即是该点最终的偏移速度和各向异性参数，对全部各向异性参数还需进行空间平滑处理。

4.根据权力要求1所述的一种各向异性三维叠前时间偏移方法，其特征在于：在F步骤中，所述的根据地下构造在不同时间深度的最大倾角、偏移速度、各向异性参数和地震道的偏移距和方位角，确定时变的三维偏移孔径和对应于偏移孔径的、用于压制偏移噪音的衰减系数是这样实现的：用一组整数

描述三维偏移孔径，其中k是有效成像区域中总的离散点数，令Δx和Δy是偏移成像结果在两个水平坐标方向的空间采样，ΔT是时间深度方向的采样，则

和

是成像区域中某离散点与地震道中心点的沿两个水平坐标方向的距离，而

是该点偏移计算的起始成像时间；令地震道中心点坐标为

和

引入时变的偏移孔径，就是在该地震道的偏移计算中，对成像区域中水平坐标为

和

的点(i＝1，k)，仅从

开始进行偏移计算；

在几个关键深度，定义地质构造在沿测线和垂直测线方向的最大倾角为和

构建如下不等式方程组：

式中

是地震道的方位角，满足上两式的最小角度αⁱ和βⁱ，就是保证该深度最大倾角为

和

的地质构造得以正确成像时，在沿方位角方向和垂直方位角方向的两个最大成像角度；

对αⁱ和βⁱ应用分段线性插值，可得到随时间深度变化的、沿方位角和垂直方位角方向的两个最大成像角度α(T)和β(T)，对不同时间深度T，偏移孔径所定义的成像区域，可看作以方位角方向为长轴方向、垂直方位角方向为短轴方向、中心在地震道的中心点上的椭圆面，其长、短半轴a(T)和b(T)分别为：

a(T)＝a₀+Δa

$a_0=\frac{{\stackrel{\&OverBar;}{V}}_{\mathrm{nmo}}\left(T\right)\&CenterDot;T}{\sin \left[2\mathrm{\&alpha;}\left(T\right)\right]}\sqrt{1+{\cos }^2\left[2\mathrm{\&alpha;}\left(T\right)\right]+{\mathrm{\&lambda;}}^2{\sin }^2\left[2\mathrm{\&alpha;}\left(T\right)\right]-2\cos \left[2\mathrm{\&alpha;}\left(T\right)\right]{\{1+{\mathrm{\&lambda;}}^2{\sin }^2[2\mathrm{\&alpha;}\left(T\right)\left]\right\}}^{1/2}}$

$\mathrm{\&Delta;a}=\frac{{5a}_0{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_{\mathrm{eff}}\left(T\right)({\tilde{a}}^2-{\mathrm{\&lambda;}}^2-1)}{2{\sin }^2\left[2\mathrm{\&alpha;}\left(T\right)\right]({\tilde{a}}^2-{\mathrm{\&lambda;}}^2+1)-4}$

$b\left(T\right)={\stackrel{\&OverBar;}{V}}_{\mathrm{nmo}}\left(T\right)\&CenterDot;T\sin \left[\mathrm{\&beta;}\left(T\right)\right]{\left({\cos }^2\right[\mathrm{\&beta;}\left(T\right)]-\frac{{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_{\mathrm{eff}}\left(T\right){\sin }^2\left[\mathrm{\&beta;}\left(T\right)\right]{\{1+{3\cos }^2[\mathrm{\&beta;}\left(T\right)\left]\right\}}^2}{{\sin }^2\left[\mathrm{\&beta;}\left(T\right)\right]{\{1+2{\cos }^2[\mathrm{\&beta;}\left(T\right)\left]\right\}}^2+{4\cos }^6\left[\mathrm{\&beta;}\left(T\right)\right]})}^{-1/2}$

式中h是地震道的半偏移距，T是用单程旅行时表达的时间深度，a₀是忽略各向异性时的长轴值，Δa是考虑各向异性时长轴的增加量，是无量纲参数，

是无量纲参数，是不同时间深度上的平均各向异性参数，无量纲，

是不同时间深度上的平均偏移速度；

设偏移计算的起始时间深度为T₀，最大时间深度为T₂，对T＝T₀，T₂循环，若a(T)≥a(T-ΔT)，a(T)不变，否则a(T)＝a(T-ΔT)；对a(T)作时间深度方向的高斯平滑，即得随时间深度变化的长半轴，同理可得随时间深度变化的短半轴；

对T＝T₀，T₂循环，对时间深度T和T+ΔT，由其对应的两组长、短半轴形成有共同的中心点(0，0)的两个椭圆，这两个椭圆形成的条带就是起始成像时间为T的成像区域了；记录落在各个条带中的离散点坐标

和

对应的三个整数

分别按

和

的大小排序，统计总点数k，就可以得到偏移距为2h、方位角为

地震道的三维时变偏移孔径；

对偏移孔径边缘区的偏移幅值进行平滑的衰减，使其与偏移孔径以外的零值不产生突变，可有效压制偏移噪音，为实现平滑的衰减，可在深度T方向和椭圆边界的法线方向联合进行系数由1到0.2的随正弦函数变化的衰减；

沿T方向衰减得到系数为

(j＝1，m₁和i＝1，k)，其中m₁是事先定义的沿垂直方向的衰减点数；沿椭圆边界的法线方向衰减得到系数为

$c_2(k,j,n_T^i)=0.8\sin (\frac{1/r^2-f}{1/r^2-1}\&CenterDot;\frac{\mathrm{\&pi;}}{2})+\mathrm{0.2,1}/r^2\&GreaterEqual;f(k,j)>1$ 和 $n_T^i=T_0/\mathrm{\&Delta;T},T_2/\mathrm{\&Delta;T}$

其中

和

是前面得到的椭圆长、短半轴，m₂是事先定义的衰减点数，T₀和T₂分别是偏移计算的起始和最大时间深度，若偏移孔径边缘区中相同的点同时被赋予衰减系数c₁和c₂，则最终的衰减系数定义为c＝(c₁·c₂)^e，其中e＝0.5+(c₁-0.2)/1.6；其它点的衰减系数c或取c₁或取c₂；

对得到的衰减系数作空间高斯平滑，平滑时边界点的值可作为1.0，平滑时计算各空间点的拉普拉斯算子值，即计算无量纲系数

φ＝c(k+1，j，i)+c(k-1，j，i)+c(k，j+1，i)

+c(k，j-1，i)+c(k，j，i+1)+c(k，j，i-1)-6c(k，j，i)

若φ大于指定值，则对该点重新平滑；

时变的三维偏移孔径和对应的衰减系数是事先计算好，存放在指定的表中，计算过程中直接到表中拾取相应的孔径参数和衰减系数。

5.根据权力要求1所述的一种各向异性三维叠前时间偏移方法，其特征在于：在G步骤中，所述的采用基于单程波算子和稳相点原理的查表法求得各向异性介质中的地震波走时和保幅成像的权系数是这样实现的：定义

$p_0=\frac{{(x-x_0)}^2+{(y-y_0)}^2}{V_{\mathrm{nmo}}^2{T}^2},$ $p_1=\frac{p_0}{p_0+1},$ $\mathrm{\&xi;}=\sqrt{\frac{p_0^3+{6p}_0^2+{9p}_0+4-{\mathrm{\&eta;}}_{\mathrm{eff}}({6p}_0^2-12p_0)}{p_0^3+{6p}_0^2+{9p}_0+4+{\mathrm{\&eta;}}_{\mathrm{eff}}({2p}_0^3+{10p}_0^2+44p_0)}}$

式中(x₀，y₀)是地震道的炮点或检波点坐标，(x，y，T)是成像点的水平和时间深度坐标，T是单程旅行时，V_nmo和η_eff是成像点处的偏移速度和各向异性参数，p₀、p₁和ξ是求得的三个无量纲参数，以

和η_eff为参数建立二维表，表中存两个无量纲数值：

$\widetilde{\mathrm{\&tau;}}=\sqrt{1-\frac{{\mathrm{\&xi;}}^2{p}_1}{1-2{\mathrm{\&eta;}}_{\mathrm{eff}}{\mathrm{\&xi;}}^2{p}_1}}+\mathrm{\&xi;}\sqrt{p_0{p}_1},$ $A=\frac{{(1-2{\mathrm{\&eta;}}_{\mathrm{eff}}{\mathrm{\&xi;}}^2{p}_1)}^2[1-({2\mathrm{\&eta;}}_{\mathrm{eff}}+1){\mathrm{\&xi;}}^2{p}_1]}{\sqrt{1+4{\mathrm{\&eta;}}_{\mathrm{eff}}{\mathrm{\&xi;}}^2{p}_1-6{\mathrm{\&eta;}}_{\mathrm{eff}}({2\mathrm{\&eta;}}_{\mathrm{eff}}+1){\mathrm{\&xi;}}^4{p}_1^2}}$

对每一成像点，由地震道的炮点坐标、成像点坐标和偏移速度计算

将和η_eff按表的间距取整，拾取相邻的4个点的值，由双线性插值得到较准确的和A_s；再由地震道的检波点坐标、成像点坐标和偏移速度计算

依据和η_eff，同样可从表中得到较准确的

和A_r；地震波从炮点到成像点、再反射到检波点的走时可简单得到为

$\mathrm{\&tau;}=({\widetilde{\mathrm{\&tau;}}}_s+{\widetilde{\mathrm{\&tau;}}}_r)T$

而计算偏移幅值的保幅成像的无量纲权系数为A_r/A_s；

对地震道的时域信号求一阶导数，通过插值计算拾取时间τ处的值，将该值乘上权系数A_r/A_s，即得到该地震道在该成像点的偏移幅值，插值计算可采用加密重采样的插值技术，以提高插值计算的精度和效率。

Images