CN104280768A - 一种适用于逆时偏移的吸收边界条件方法 - Google Patents

Abstract

本发明是一种适用于逆时偏移的吸收边界条件方法，叠前与深度域建模后炮点波场模拟，沿时间方向二阶差分，每一个时刻t对计算区域四个水平方向分别增加n层的计算层；采用Higdon法对每个计算层循环计算，更新波场；对计算区域顶面和底面方向用二阶混合差分近似吸收边界条件公式计算，采用相同的步骤进行检波点波场模拟；利用互相关成像对炮点波场与检波点波场互相关计算得到一炮偏移结果，并完成所有炮的地震数据偏移计算，叠加形成剖面。本发明能有效处理入射角度较大的入射波的反射，在没有增加过多计算的情况下有效消除虚假反射导致的成像噪音，计算效率高、吸收效果好、以及易于实现。

Description

一种适用于逆时偏移的吸收边界条件方法

技术领域

本发明涉及反射波地震数据处理中一种地震波偏移成像技术，具体的说是一种针对地震波逆时偏移(RTM)成像的吸收边界条件方法。

背景技术

逆时偏移(RTM)能够很好的解决速度急剧变化、以及构造异常复杂的地区的成像问题，目前，该技术在国际上的主要地震探测市场，已经成为偏移成像方法的重点，也是各个地球物理公司重点发展的一种针对复杂构造成像的特色技术。逆时偏移(RTM)基于全声波方程，无成像倾角限制，能够有效地成像回转波、棱镜波、振幅海面反射(Ghosts)、以及多次波，适合各向异性偏移。运用合适的吸收边界条件对于逆时偏移(RTM)至关重要，否则将影响成像的质量和计算的效率。因此，有必要设计一种适用于逆时偏移(RTM)的吸收边界条件，消除虚假反射导致的成像噪音。

Clayton and Engquist(1977)(also Engquist and Majda,1977)提出了基于声波方程以及弹性波方程旁轴近似的吸收边界条件，另外一种类似的、同时得到广泛应用的方法是由Reynolds(1978)提出的。这类方法的不足之处在于不能很好的处理入射角度较大的入射波的反射问题。Berenger(1994)设计了一种有效的吸收边界条件，即最佳匹配层(PML)，但是其不足之处在于计算量相对来说较大，而逆时偏移(RTM)本身就是一个计算量巨大的过程。Higdon(1986,1987)给出了一种方法，直接运用差分近似来消除反射波振幅，其理论与有效吸收非零角度的入射波的解析条件具有一致性。该方法计算效率不高，吸收效果不好，不易于实现。

发明内容

本发明目的是提供一种计算效率高，吸收效果好，易于实现的适用于逆时偏移的吸收边界条件方法。

本发明通过以下具体步骤实现：

1）采集地震数据并完成叠前处理与深度域建模；

2）进行炮点有限差分波场模拟，沿时间方向的二阶差分，在每一个时刻t，针对计算区域的四个水平方向，分别增加n层的计算层；同时对于每一个计算层，循环使用基于Higdon法的一阶混合差分近似吸收边界条件公式计算，更新波场U1、U2、U3；

所述的增加n层的计算层，n取值范围是5<n<10。

所述的增加n层的计算层，n取值范围是最佳取值为8，即增加8个计算层的时候，在不增加过多计算量的前提下，吸收效果最好。

所述的更新波场U1、U2、U3分别是延拓过程中前一时刻、当前时刻、以及下一时刻的波场。

3）针对计算区域的顶面和底面两个方向，分别增加n层的计算层，同时对于每一个计算层，循环使用基于Higdon法的二阶混合差分近似吸收边界条件公式计算，更新波场U1、U2、U3；

所述的n取值范围和更新波场U1、U2、U3与步骤2）一致。

4）按照步骤2）和步骤3）进行检波点有限差分波场模拟；

5）利用互相关成像条件，对步骤2）、步骤3）得到的炮点波场与步骤4）得到的检波点波场进行互相关计算得到一炮偏移结果，完成所有炮的地震数据偏移计算，将结果叠加形成逆时偏移成像剖面。

本发明的二阶吸收边界条件公式，将该公式作用于波场表达式，使得地震波在边界得到有效吸收，特别是解决入射角度较大的入射波的反射问题的能力要明显优于传统方法，能够有效消除虚假反射导致的逆时偏移成像噪音问题。

本发明针对不同的计算方向，采用不同阶数的吸收边界条件公式，顶面和底面采用精度较高的二阶吸收边界条件公式，水平方向采用精度较低、效率较高的一阶吸收边界条件公式来更新波场U1、U2、U3。

在有限差分波场模拟时，本发明二阶吸收边界条件公式用于波场表达式，使得地震波在边界得到有效吸收，特别是解决入射角度较大的入射波的反射问题的能力要明显优于传统方法，能够有效消除虚假反射导致的逆时偏移成像噪音问题。同时，本发明针对不同的计算方向，采用不同的吸收边界条件公式，顶面和底面采用精度较高的二阶吸收边界条件公式，水平方向采用精度较低、效率较高的一阶吸收边界条件公式来更新波场U1、U2、U3，在保证计算效率的同时，最大程度消除边界入射波的反射问题，为实现高精度的三维地震波逆时偏移(RTM)成像铺平道路。

图5、图6为某三维实验数据Kirchhoff叠前深度偏移结果与逆时偏移(RTM)结果的深度切片对比，其中逆时偏移(RTM)结果使用的是本发明提出的吸收边界条件，可以看到逆时偏移(RTM)成像效果更精确，特别是对于高速体边界的成像效果，Kirchhoff叠前深度偏移结果并不是非常清晰，而逆时偏移成像结果则非常精确的刻画了高速体的边界，本发明在资料的计算中发挥了重要作用，可以应用于逆时偏移实际生产项目。

附图说明

图1常速介质中三维合成数据单炮逆时偏移(RTM)结果，(a)传统的吸收边界条件，(b)本发明提出的吸收边界条件；

图2常速介质中三维波场传播x-z方向波场快照，(a)没有应用吸收边界条件，(b)传统的吸收边界条件，(c)本发明提出的吸收边界条件；

图3常速介质中三维合成数据单炮逆时偏移(RTM)计算过程中炮点波场快照，(a)传统的吸收边界条件，(b)本发明提出的吸收边界条件；

图4检波点有限差分波场模拟波场快照

图5中亚地区某三维实际数据Kirchhoff叠前深度偏移结果深度切片；

图6中亚地区某三维实际数据逆时偏移(RTM)结果的深度切片。

具体实施方式

结合图例说明具体实施方式：

本发明在Higdon(1986,1987)给出的吸收边界条件的基础上，发展了一种适用于逆时偏移(RTM)的吸收边界条件方法，该方法计算效率高，吸收效果好，易于实现，能够很好的适用于逆时偏移(RTM)。

本发明具体实施方式为：

1）采集地震数据并完成叠前处理与深度域建模。

2）进行炮点有限差分波场模拟，沿时间方向的二阶差分，在每一个时刻t，针对计算区域的四个水平方向，分别增加n层的计算层；同时对于每一个计算层，循环使用基于Higdon法的一阶混合差分近似吸收边界条件公式计算，更新波场U1、U2、U3。

所述的增加n层的计算层，n取值范围是5<n<10。

所述的增加n层的计算层，n取值范围是最佳取值为8，即增加8个计算层的时候，在不增加过多计算量的前提下，吸收效果最好。

所述的更新波场U1、U2、U3分别是延拓过程中前一时刻、当前时刻、以及下一时刻的波场。

波动方程的表达式为：

$(\frac{\&PartialD;}{{\&PartialD;x}^2}+\frac{\&PartialD;}{{\&PartialD;y}^2}+\frac{\&PartialD;}{{\&PartialD;z}^2}-\frac{1}{V^2}\frac{\&PartialD;}{{\&PartialD;t}^2})u=0$

Higdon(1986,1987)给出的吸收边界条件在x方向、向外传播的方程如下：

$\underset{i=1}{\overset{p}{\mathrm{\&Pi;}}}{B}_i=\underset{i=1}{\overset{p}{\mathrm{\&Pi;}}}(\frac{\&PartialD;}{\&PartialD;x}+\frac{\cos {\mathrm{\&theta;}}_1}{c}\frac{\&PartialD;}{\&PartialD;t}+{\mathrm{\&epsiv;}}_i)u=0$

其中，B是吸收边界公式，该公式能够很好的吸收入射角度为θ₁,θ₂,…,θ_p的平面波的任意线性组合。ε_i是针对低频分量的吸收衰减因子。将方程(2)写成有限差分的形式，可以得到如下的方程：

$B_i=\frac{I-D^{-1}}{\mathrm{\&Delta;x}}((1-a)I+a{K}^{-1})$

$+\frac{\cos {\mathrm{\&theta;}}_i}{c}\left(\frac{I-K^{-1}}{\mathrm{\&Delta;t}}\right)((1-b)I+{\mathrm{bD}}^{-1})+{\mathrm{\&epsiv;}}_i$

其中，系数a和系数b是空间与时间差分的加权系数。公式I、D和K的定义表达式为：

Iuⁿ(i,j,k)=uⁿ(i,j,k)

Duⁿ(i,j,k)=uⁿ(i+1,j,k)

Kuⁿ(i,j,k)=uⁿ⁺¹(i,j,k)

方程(3)能够写成如下形式：

B_i=I-α_iK^-1-β_iD^-1-γ_iD^-1K^-1

其中：

${\mathrm{\&alpha;}}_i=\frac{(a-g_i(1-b))}{(a-1-g_i(1-b)-{\mathrm{\&epsiv;}}_i\mathrm{\&Delta;l})}$

${\mathrm{\&beta;}}_i=\frac{(a-1+g_{i}b)}{(a-1-g_i(1-b)-{\mathrm{\&epsiv;}}_i\mathrm{\&Delta;l})}$

${\mathrm{\&gamma;}}_i=\frac{(-a-g_{i}b)}{(a-1-g_i(1-b)-{\mathrm{\&epsiv;}}_i\mathrm{\&Delta;}l)}$

参数g_i的定义如下：

$g_i=\frac{\cos {\mathrm{\&theta;}}_i}{c}\frac{\mathrm{\&Delta;l}}{\mathrm{\&Delta;t}}$

对于一阶(p=1)，波场u在第n个时刻的第m个点的波场值可以由沿着边界方向的相邻点得到：

uⁿ(m,j,k)=α₁u^n-1(m,j,k)+β₁uⁿ(m-1,j,k)

+γ₁u^n-1(m-1,j,k)

上式即对于水平方向应用的一阶吸收边界条件公式。

3）针对计算区域的顶面和底面两个方向，分别增加n层的计算层，同时对于每一个计算层，循环使用基于Higdon法的二阶混合差分近似吸收边界条件公式计算，更新波场U1、U2、U3。

所述的n取值范围和更新波场U1、U2、U3与步骤2）一致。

图1为常速介质中三维合成数据单炮逆时偏移(RTM)结果，其中a使用了传统的吸收边界条件，b使用了本发明提出的吸收边界条件，可以看到由于顶面吸收效果不理想，导致实际反射层的上部出现了一个影子，而这是不应该出现的噪音，b的吸收效果相对于a较好。可以看出，顶面需要采用精度高于一阶的二阶吸收边界条件公式。

对于二阶吸收边界来说，公式变成如下的形式：

B=B₁B₂

=I-(α₁+α₂)K^-1+α₁α₂K^-2-(β₁+β₂)D^-1

+(α₁β₂+α₂β₁-γ₁-γ₂)D^-1K^-1

+(α₁γ₂+α₂γ₁)D^-1K^-2+β₁β₂D^-2

+(β₁γ₂+β₂γ₁)D^-2K^-1+γ₁γ₂D^-2K^-2

将该公式直接应用于波场表达式可以得到如下公式：

uⁿ(m,j,k)=(α₁+α₂)u^n-1(m,j,k)+α₁α₂u^n-2(m,j,k)

-(β₁+β₂)uⁿ(m-1,j,k)

+(α₁β₂+α₂β₁-γ₁-γ₂)u^n-1(m-1,j,k)

+(α₁γ₂+α₂γ₁)u^n-2(m-1,j,k)+β₁β₂uⁿ(m-2,j,k)

+(β₁γ₂+β₂γ₁)u^n-1(m-2,j,k)+γ₁γ₂u^n-2(m-2,j,k)

上式即对于顶面和底面应用的二阶吸收边界条件公式。

图2为常速介质中三维波场传播x-z方向波场快照，其中a没有应用吸收边界条件，b使用了传统的吸收边界条件，c使用了本发明提出的吸收边界条件，可以看到，a中的边界反射非常明显，使用了传统的吸收边界条件以后，反射明显减弱，但是仍然能够看到反射波的存在，而使用了本发明提出的吸收边界条件以后，在同一时刻的波场快照中，已经几乎看不到反射波的存在。

图3为常速介质中三维合成数据单炮逆时偏移(RTM)计算过程中炮点波场快照，a使用了传统的吸收边界条件，b使用了本发明提出的吸收边界条件，a中可以看到比较明显的顶面反射波向下传播，而这个反射波在逆时偏移(RTM)结果中将形成噪音。而使用了本发明提出的吸收边界条件以后，在同一时刻的波场快照中，顶面反射波几乎消失。

4）按照步骤2）和步骤3）进行检波点有限差分波场模拟，图4为检波点有限差分波场模拟的波场快照结果。

5）利用互相关成像条件，对步骤2）、步骤3）得到的炮点波场与步骤4）得到的检波点波场进行互相关计算得到一炮偏移结果，完成所有炮的地震数据偏移计算，将结果叠加形成逆时偏移成像剖面。

图6为中亚地区某三维实际数据逆时偏移(RTM)结果的深度切片，使用了本发明提出的吸收边界条件。

本发明所述的一种适用于逆时偏移(RTM)的吸收边界条件方法能够解决三维复杂构造成像问题，且具有如下特点：

（1）本发明所述的方法采用了二阶的吸收边界公式，同时应用于计算边界的8层吸收区域，能够有效处理入射角度较大的入射波的反射问题，因此能够有效消除逆时偏移(RTM)剖面上虚假反射导致的成像噪音，有利于解决三维复杂构造的成像问题。

（2）本发明所述的方法只增加了8层二阶吸收边界公式的计算，计算效率高，吸收效果好，易于实现，能够很好的适用于逆时偏移(RTM)商业化软件的开发。

Claims (4)

1.一种适用于逆时偏移的吸收边界条件方法，特点是通过以下步骤实现：

1）采集地震数据并完成叠前处理与深度域建模；

2）进行炮点有限差分波场模拟，沿时间方向的二阶差分，在每一个时刻t，针对计算区域的四个水平方向，分别增加n层的计算层；同时对于每一个计算层，循环使用基于Higdon法的一阶混合差分近似吸收边界条件公式计算，更新波场U1、U2、U3；

3）针对计算区域的顶面和底面两个方向，分别增加n层的计算层，同时对于每一个计算层，循环使用基于Higdon法的二阶混合差分近似吸收边界条件公式计算，更新波场U1、U2、U3；

所述的n取值范围和更新波场U1、U2、U3与步骤2）一致；

4）按照步骤2）和步骤3）进行检波点有限差分波场模拟；

5）利用互相关成像条件，对步骤2）、步骤3）得到的炮点波场与步骤4）得到的检波点波场进行互相关计算得到一炮偏移结果，完成所有炮的地震数据偏移计算，将结果叠加形成逆时偏移成像剖面。

2.根据权利要求1的方法，特点是步骤2）所述的增加n层的计算层，n取值范围是5<n<10。

3.根据权利要求1的方法，特点是步骤2）所述的增加n层的计算层，n取值范围是最佳取值为8，在增加8个计算层的时候，在不增加过多计算量的前提下，吸收效果最好。

4.根据权利要求1的方法，特点是步骤2）所述的更新波场U1、U2、U3分别是延拓过程中前一时刻、当前时刻、以及下一时刻的波场。