CN104280768A - 一种适用于逆时偏移的吸收边界条件方法 - Google Patents
一种适用于逆时偏移的吸收边界条件方法 Download PDFInfo
- Publication number
- CN104280768A CN104280768A CN201310293637.XA CN201310293637A CN104280768A CN 104280768 A CN104280768 A CN 104280768A CN 201310293637 A CN201310293637 A CN 201310293637A CN 104280768 A CN104280768 A CN 104280768A
- Authority
- CN
- China
- Prior art keywords
- wave field
- boundary condition
- calculation
- absorbing boundary
- layer
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Landscapes
- Geophysics And Detection Of Objects (AREA)
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分别是延拓过程中前一时刻、当前时刻、以及下一时刻的波场。
波动方程的表达式为:
Higdon(1986,1987)给出的吸收边界条件在x方向、向外传播的方程如下:
其中,B是吸收边界公式,该公式能够很好的吸收入射角度为θ1,θ2,…,θp的平面波的任意线性组合。εi是针对低频分量的吸收衰减因子。将方程(2)写成有限差分的形式,可以得到如下的方程:
其中,系数a和系数b是空间与时间差分的加权系数。公式I、D和K的定义表达式为:
Iun(i,j,k)=un(i,j,k)
Dun(i,j,k)=un(i+1,j,k)
Kun(i,j,k)=un+1(i,j,k)
方程(3)能够写成如下形式:
Bi=I-αiK-1-βiD-1-γiD-1K-1
其中:
参数gi的定义如下:
对于一阶(p=1),波场u在第n个时刻的第m个点的波场值可以由沿着边界方向的相邻点得到:
un(m,j,k)=α1un-1(m,j,k)+β1un(m-1,j,k)
+γ1un-1(m-1,j,k)
上式即对于水平方向应用的一阶吸收边界条件公式。
3)针对计算区域的顶面和底面两个方向,分别增加n层的计算层,同时对于每一个计算层,循环使用基于Higdon法的二阶混合差分近似吸收边界条件公式计算,更新波场U1、U2、U3。
所述的n取值范围和更新波场U1、U2、U3与步骤2)一致。
图1为常速介质中三维合成数据单炮逆时偏移(RTM)结果,其中a使用了传统的吸收边界条件,b使用了本发明提出的吸收边界条件,可以看到由于顶面吸收效果不理想,导致实际反射层的上部出现了一个影子,而这是不应该出现的噪音,b的吸收效果相对于a较好。可以看出,顶面需要采用精度高于一阶的二阶吸收边界条件公式。
对于二阶吸收边界来说,公式变成如下的形式:
B=B1B2
=I-(α1+α2)K-1+α1α2K-2-(β1+β2)D-1
+(α1β2+α2β1-γ1-γ2)D-1K-1
+(α1γ2+α2γ1)D-1K-2+β1β2D-2
+(β1γ2+β2γ1)D-2K-1+γ1γ2D-2K-2
将该公式直接应用于波场表达式可以得到如下公式:
un(m,j,k)=(α1+α2)un-1(m,j,k)+α1α2un-2(m,j,k)
-(β1+β2)un(m-1,j,k)
+(α1β2+α2β1-γ1-γ2)un-1(m-1,j,k)
+(α1γ2+α2γ1)un-2(m-1,j,k)+β1β2un(m-2,j,k)
+(β1γ2+β2γ1)un-1(m-2,j,k)+γ1γ2un-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分别是延拓过程中前一时刻、当前时刻、以及下一时刻的波场。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310293637.XA CN104280768B (zh) | 2013-07-12 | 2013-07-12 | 一种适用于逆时偏移的吸收边界条件方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310293637.XA CN104280768B (zh) | 2013-07-12 | 2013-07-12 | 一种适用于逆时偏移的吸收边界条件方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104280768A true CN104280768A (zh) | 2015-01-14 |
CN104280768B CN104280768B (zh) | 2017-03-15 |
Family
ID=52255830
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201310293637.XA Active CN104280768B (zh) | 2013-07-12 | 2013-07-12 | 一种适用于逆时偏移的吸收边界条件方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104280768B (zh) |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105447225A (zh) * | 2015-11-06 | 2016-03-30 | 中国海洋大学 | 一种应用于声波有限差分数值模拟的组合吸收边界条件 |
CN106054242A (zh) * | 2016-05-04 | 2016-10-26 | 中国地质大学(北京) | 三维各向异性衰减介质波场模拟方法 |
CN108828668A (zh) * | 2018-03-27 | 2018-11-16 | 中国石油天然气集团有限公司 | 一种叠前时间偏移数据处理方法及装置 |
CN109188517A (zh) * | 2018-09-03 | 2019-01-11 | 中国海洋大学 | 基于Higdon余弦型加权的混合吸收边界条件方法 |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6446007B1 (en) * | 1998-11-05 | 2002-09-03 | Exxonmobil Upstream Research Company | Method for controlled-amplitude prestack time migration of seismic data |
US6687659B1 (en) * | 2000-03-24 | 2004-02-03 | Conocophillips Company | Method and apparatus for absorbing boundary conditions in numerical finite-difference acoustic applications |
CN102156296A (zh) * | 2011-04-19 | 2011-08-17 | 中国石油大学(华东) | 地震多分量联合弹性逆时偏移成像方法 |
CN102879816A (zh) * | 2012-07-17 | 2013-01-16 | 中国科学院地质与地球物理研究所 | 一种地震多次波偏移方法 |
CN102906599A (zh) * | 2010-02-10 | 2013-01-30 | 埃克森美孚上游研究公司 | 用于全波场反演和逆时偏移中的地下参数估计的方法 |
CN103091710A (zh) * | 2013-01-15 | 2013-05-08 | 中国石油天然气股份有限公司 | 一种逆时偏移成像方法及装置 |
-
2013
- 2013-07-12 CN CN201310293637.XA patent/CN104280768B/zh active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6446007B1 (en) * | 1998-11-05 | 2002-09-03 | Exxonmobil Upstream Research Company | Method for controlled-amplitude prestack time migration of seismic data |
US6687659B1 (en) * | 2000-03-24 | 2004-02-03 | Conocophillips Company | Method and apparatus for absorbing boundary conditions in numerical finite-difference acoustic applications |
CN102906599A (zh) * | 2010-02-10 | 2013-01-30 | 埃克森美孚上游研究公司 | 用于全波场反演和逆时偏移中的地下参数估计的方法 |
CN102156296A (zh) * | 2011-04-19 | 2011-08-17 | 中国石油大学(华东) | 地震多分量联合弹性逆时偏移成像方法 |
CN102879816A (zh) * | 2012-07-17 | 2013-01-16 | 中国科学院地质与地球物理研究所 | 一种地震多次波偏移方法 |
CN103091710A (zh) * | 2013-01-15 | 2013-05-08 | 中国石油天然气股份有限公司 | 一种逆时偏移成像方法及装置 |
Non-Patent Citations (2)
Title |
---|
任济时: ""高阶Higdon吸收边界条件的直接算法及其评价"", 《电子学报》 * |
王守东: ""声波方程完全匹配层吸收边界"", 《石油地球物理勘探》 * |
Cited By (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105447225A (zh) * | 2015-11-06 | 2016-03-30 | 中国海洋大学 | 一种应用于声波有限差分数值模拟的组合吸收边界条件 |
CN106054242A (zh) * | 2016-05-04 | 2016-10-26 | 中国地质大学(北京) | 三维各向异性衰减介质波场模拟方法 |
CN108828668A (zh) * | 2018-03-27 | 2018-11-16 | 中国石油天然气集团有限公司 | 一种叠前时间偏移数据处理方法及装置 |
CN108828668B (zh) * | 2018-03-27 | 2020-04-10 | 中国石油天然气集团有限公司 | 一种叠前时间偏移数据处理方法及装置 |
CN109188517A (zh) * | 2018-09-03 | 2019-01-11 | 中国海洋大学 | 基于Higdon余弦型加权的混合吸收边界条件方法 |
Also Published As
Publication number | Publication date |
---|---|
CN104280768B (zh) | 2017-03-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Archuleta et al. | Effects of fault finiteness on near-source ground motion | |
Du et al. | Polarity reversal correction for elastic reverse time migration | |
Yan et al. | Improving the efficiency of elastic wave-mode separation for heterogeneous tilted transverse isotropic media | |
CN103091710A (zh) | 一种逆时偏移成像方法及装置 | |
CN103513277B (zh) | 一种地震地层裂隙裂缝密度反演方法及系统 | |
CN104280768A (zh) | 一种适用于逆时偏移的吸收边界条件方法 | |
CN102636809A (zh) | 一种传播角度域共成像点道集的生成方法 | |
Yang et al. | A strong stability-preserving predictor-corrector method for the simulation of elastic wave propagation in anisotropic media | |
CN105447225B (zh) | 一种应用于声波有限差分数值模拟的组合吸收边界条件 | |
Zhou et al. | An iterative factored topography-dependent eikonal solver for anisotropic media | |
Wei‐Hong et al. | Elastic reverse time migration based on wavefield separation | |
CN102830424B (zh) | 一种检波器组合参数计算方法 | |
Heng et al. | Forward Modeling of VTI Media in the Frequency‐Space Domain Based on an Average‐Derivative Optimal Method | |
CN102998702B (zh) | 保幅平面波叠前深度偏移方法 | |
Kamal et al. | 3D basin-shape ratio effects on frequency content and spectral amplitudes of basin-generated surface waves and associated spatial ground motion amplification and differential ground motion | |
WANG et al. | Multi‐Azimuth Three‐Component Surface Seismic Modeling in Cracked Monoclinic Media | |
Wang et al. | Numerical simulation of elastic wave equation and analysis of wave field characteristics in 2-D VTI medium | |
CN102313902B (zh) | 一种基于切比雪夫展开的广义屏叠前深度偏移方法 | |
Fang et al. | An unsplit complex frequency-shifted perfectly matched layer for second-order acoustic wave equations | |
Chanu et al. | Along‐strike variation in the shear wave crustal structure of the NE Himalayan and Indo‐Burmese arc: Evidence based on surface wave dispersion analysis | |
CN104216012A (zh) | 三维Born-Kirchhoff变步长插值成像方法 | |
Gao‐Xiang et al. | A quantitative analysis method for the seismic geological complexity of near surface | |
Li et al. | Numerical simulation of seismic wave in elastic and viscoelastic TTI media | |
Zhang et al. | Frequency-space domain high-order modeling based on an average-derivative optimal method | |
Gu et al. | An excitation potential imaging condition for elastic reverse time migration |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant |