CN103149586A - 一种倾斜层状粘弹性介质中波场正演模拟方法 - Google Patents
一种倾斜层状粘弹性介质中波场正演模拟方法 Download PDFInfo
- Publication number
- CN103149586A CN103149586A CN2013100443842A CN201310044384A CN103149586A CN 103149586 A CN103149586 A CN 103149586A CN 2013100443842 A CN2013100443842 A CN 2013100443842A CN 201310044384 A CN201310044384 A CN 201310044384A CN 103149586 A CN103149586 A CN 103149586A
- Authority
- CN
- China
- Prior art keywords
- wave
- ripple
- coordinate system
- reflection
- interface
- 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
Images
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种倾斜层状粘弹性介质中波场正演模拟方法,该方法针对地质模型为含倾斜界面的粘弹性层状介质,将点源激发的球面波分解为平面波,利用矢量波动方程不依赖于坐标系这一性质,在固定坐标系和随倾斜界面变化的动坐标系下交替地讨论波的传播。首先由传统的反射和透射系数推导出单个倾斜界面上平面波的反射和透射系数,得到该界面上的反射及透射波;然后递推求得平面波经过多个倾斜界面时的反射和透射系数。最后由传播后的各平面波合并得到点源激发波场,合并时采用了快速的振荡积分计算方法和一条适合粘弹性介质的积分路径。
Description
技术领域
本发明属于地震勘探技术领域,涉及一种模拟方法,尤其是一种倾斜层状粘弹性介质中波场正演模拟方法,其是一种针对含倾斜界面的层状粘弹性介质中地震波场高效精确的正演模拟技术。
背景技术
波场模拟在油气勘探中占有重要地位,尽管已存在多种类型的波场模拟方法,然而,针对要解决的问题的要求,发展波场模拟理论及方法一直是很活跃的研究领域。另外,在油气勘探等领域,为了解释地震图,或进行数据反演,常常需要仅模拟指定类型的波,譬如仅模拟一次反射波,或仅模拟转换波等。常用的波场模拟方法有多种类型,如基于射线理论的方法,数值解方法,广义反射率法以及其推广而得到的反/透射矩阵法等。这些方法各有其优缺点。
射线方法能适用于较复杂的介质模型,可根据事先给定的射线码模拟感兴趣的类型的波,且计算速度快;不足之处是,它只是波动方程的高频近似,用于模拟复杂储层的地震响应具有局限性。采用数值法求解波动方程,是复杂介质中波场模拟的有效途径。常用的数值方法有:有限差分法,有限元法,伪谱法和谱元法等,这类方法可模拟地震波全波场,精度高。但计算量巨大,不能像射线类方法,有选择性的仅模拟我们感性趣的类型的波。
层状介质中波传播理论及波场模拟技术一直受到人们的重视。Aki和Müller等人发展的反射率法,能快速精确地模拟水平层状粘弹性介质的全波场或指定类型的波。Berkhout提出了波场外推法可用于带倾斜界面的弹性介质波场模拟。Chen等人发展的全局广义反/透射矩阵法和边界元法,能模拟带不规则界面的层状弹性介质中全波场。但这些方法不适用于含倾斜界面的层状粘弹性介质。
发明内容
针对上述问题,本发明的目的在于针对含倾斜波阻抗界面的层状粘弹性介质,提出一种倾斜层状粘弹性介质中波场正演模拟方法,该方法可以选择性模拟感兴趣的各种类型的波,如一次反射P波、S波和转换波,且模拟的波场非常精确。
本发明的目的是通过以下技术方案来解决的:
该种倾斜层状粘弹性介质中波场正演模拟方法,包括以下步骤:
1)将点源激励的球面波分解为平面波;
2)将平面波在倾斜层状粘弹性介质中传播,遇到界面时进行如下处理:以该界面为轴,界面与全局坐标系Z轴的交点为原点建立局部坐标系,垂直于界面向下,将全局坐标系中的入射波转换到局部坐标系中;在局部坐标系下求得反射波和透射波;再将局部坐标系下的反射波和透射波转换到全局坐标系下,得到全局坐标系下的反射系数和透射系数,以及反射波和透射波;在每个界面上根据需要选取波的类型,P波入射时能够选取的波包括反射PP波、反射PS波、透射PP波、透射PS波;S波入射时能够选取的波包括反射SP波、反射SS波、透射SP波、透射SS波;从而模拟所选取的波类型,得到模拟的平面波;
3)由步骤2)模拟得到的平面波根据下式合成点源激发的波场:
式中,Γ为积分路径;N为平面波的个数;Ri(p)为平面波在多个界面上的透射和反射系数综合;j是虚数单位,ω是角频率;p1i、q1i分别为第i个界面反射回第1层的平面波的水平和垂直慢度;p、q分别为点源处平面波的水平和垂直慢度,x和z表示位置坐标。
进一步的,上述步骤2)中,模拟震源为P波时的一次反射P波,设震源和检波器均在第一层,将震源处入射波逐层向下传播的递推步骤为:
④反射波φ′i,i逐层向上递推到第一层,遇到界面时执行①-③步,但只取透射波;
⑤透射波φi+1继续向下递推,遇到下一个界面时回到步骤①。
进一步,在上述步骤2)中,模拟其它类型的波采用与模拟震源为P波时的一次反射P波的相同方法和步骤进行。
在上述步骤3)中,式中积分方法采用Frazer和Gettrust提出的GFM积分法。
本发明具有以下几点有益效果:
(1)本发明适用于带倾斜界面的层状粘弹性介质,各界面的倾斜方向和倾斜角度可以不同,只需满足界面不相交的条件。
(2)本发明能正确模拟粘弹性介质中波的衰减特征,相对于有限差分法,具有较高的计算效率,更适合层较厚、且层间参数变化剧烈情况。
(3)本发明非常灵活,可以分别模拟各种波,如一次反射P波、SV波、SH波和PS、SP转换波,也可模拟我们感兴趣的多次波,只需将希望模拟的波类型和相应的层次以编码形式给定,入射波入射角超过临界角时会出现首波。
(4)本发明模拟的波场非常接近解析解,可为数值方法解提供一个对比参照。所以本发明在地震资料解释和反演中有重要应用价值和应用前景。
附图说明
图1是本发明针对的倾斜层状粘弹性介质模型;
图2是波入射到界面上时的反射和透射示意图;
图3是本发明中平面波合成点源激励波场时的积分路径示意图。
具体实施方式
本发明针对的地质模型为含倾斜界面的粘弹性层状介质,将点源激发的球面波分解为平面波,利用矢量波动方程不依赖于坐标系这一性质,在固定坐标系和随倾斜界面变化的动坐标系下交替地讨论波的传播。首先由传统的反射和透射系数推导出单个倾斜界面上平面波的反射和透射系数,得到该界面上的反射及透射波;然后递推求得平面波经过多个倾斜界面时的反射和透射系数。最后由传播后的各平面波合并得到点源激发波场,合并时采用了快速的振荡积分计算方法和一条适合粘弹性介质的积分路径。
下面结合附图和实施例对本发明进行详细的描述。
首先本发明的主要符号定义如表1所示:
表1本发明中主要符号约定
图1为倾斜层状粘弹性介质模型示意图,各层介质的密度为ρi(i=1,2,…,N-1)、品质因子为Qi、P波速度为αi、S波速度为βi,各倾斜界面的倾角为θi,各倾斜界面与垂直方向的Z轴交点为zi。点源在第一次层,检波器可以在任意位置。以水平方向为X轴,垂直向下方向为Z轴,建立全局坐标系XOZ。
1)将点源激励的球面波分解为平面波
根据Weyl积分,均匀介质中点源激励的球面波表示为一系列平面波的叠加,二维情况下每个平面波的表达式为:
φ=exp[jω(px+qz)]
上式中省略了时间因子exp(-jωt),j是虚数单位,ω是角频率,p是水平慢度,q是垂直慢度,x和z表示位置坐标。
2)将平面波在倾斜层状粘弹性介质中传播,遇到界面时进行如下处理:以该界面为轴,界面与全局坐标系Z轴的交点为原点建立局部坐标系,垂直于界面向下,将全局坐标系中的入射波转换到局部坐标系中;在局部坐标系下求得反射波和透射波;再将局部坐标系下的反射波和透射波转换到全局坐标系下,得到全局坐标系下的反射系数和透射系数,以及反射波和透射波。在每个界面上可根据需要选取波的类型,P波入射时能够选取的波包括反射PP波、反射PS波、透射PP波、透射PS波;S波入射时能够选取的波包括反射SP波、反射SS波、透射SP波、透射SS波;从而仅模拟我们感兴趣的波类型。例如需模拟震源为P波时的一次反射P波,设震源和检波器均在第一层,将震源处入射波逐层向下传播的递推步骤为:
④反射波φ′i,i逐层向上递推到第一层,遇到界面时执行①-③步,但只取透射波(上行波);
⑤透射波φi+1继续向下递推,遇到下一个界面时回到步骤①。
3)由平面波合成点源激发的波场;
其中,Ri(p)为平面波在多个界面上的透射和反射系数综合,如模拟一次P-P反射波时 p和q为震源处入射平面波的水平和垂直慢度;p1i和q1i为平面波经第i个界面反射后到达检波器所在层时的水平和垂直慢度,它们是入射波水平慢度p的函数。Γ为积分路径,本发明选择如图2所示的积分路径,1/αi和1/βi是P波和S波的慢度,路径在第一象限中的拐点B的实部应稍大于max(1/βi),对单个接收器,θ=tan-1(xzr-zs|)时积分收敛最快(x、zr、zs分别是偏移距和接收器与源的z轴坐标),如是多个接收器同时接收,可取某个折中值。线段OB的倾角应使得OB段非常靠近第一象限内最小极点的下方。积分方法采用Frazer和Gettrust提出GFM(Generalization ofFilon’s Method)积分法。
步骤2)中入射波在界面上产生的反射波和透射波的具体计算方法为:
如图2,设有一水平慢度为p1的P波φ1从原点出发在介质I中向下传播,振幅为1,则该平面波可写为:
φ1=exp[jω(p1x+q1z)],
将入射P波由全局坐标系转换到局部坐标系,得到局部坐标系下的入射P波为:
其中:
A1=exp[jω(p1x0+q1z0)],
相对于局部坐标系,界面是水平的,这种情况的反射和透射系数计算方法可参照相关文献(Muller,1985)。将所有反射波和透射波由局部坐标系转换到全局坐标系,整理变换后各种波的表达式则得到全局坐标系下各反射波和透射波为:
式中:
对于波由下方向上入射到界面的情况,设全局坐标系下有一P波φ′2在介质II中向上传播:
φ′2=exp[jω(p2x-q2z)].
则在界面处会产生反射回介质II的下行P波φ2和SV波Ψ2,以及透射到介质I中的上行P波φ′1和SV波Ψ′1。按照P波从界面上方入射时的处理流程:先将入射波φ′2转换到局部坐标系下,在局部坐标系里求得各反射波和透射波的传统反射和透射系数;再将局部坐标系下各反射波和透射波转换到全局坐标系,得到全局坐标系下反射波和透射波的表达式分别为:
其中:
A′=exp[jω(p2x0-q2z0)]
同理,对于SV波从倾斜界面上方(或下方)入射的情况,其求解与P波入射时相同。
上述步骤3)中GFM积分方法为:
为了利用GFM积分法,先将步骤3)中积分式改写成如下形式(这里只给出求和号里的一项):
P=∫Γf(p)exp[sg(p)]dp
其中:
GFM梯形法则为:
这里b-a等于积分步长,对任意函数h,δh=h(b)-h(a)。
Claims (4)
1.一种倾斜层状粘弹性介质中波场正演模拟方法,其特征在于,包括以下步骤:
1)将点源激励的球面波分解为平面波;
2)将平面波在倾斜层状粘弹性介质中传播,遇到界面时进行如下处理:以该界面为轴,界面与全局坐标系Z轴的交点为原点建立局部坐标系,垂直于界面向下,将全局坐标系中的入射波转换到局部坐标系中;在局部坐标系下求得反射波和透射波;再将局部坐标系下的反射波和透射波转换到全局坐标系下,得到全局坐标系下的反射系数和透射系数,以及反射波和透射波;在每个界面上根据需要选取波的类型,P波入射时能够选取的波包括反射PP波、反射PS波、透射PP波、透射PS波;S波入射时能够选取的波包括反射SP波、反射SS波、透射SP波、透射SS波;从而模拟所选取的波类型,得到模拟的平面波;
3)由步骤2)模拟得到的平面波根据下式合成点源激发的波场:
式中,Γ为积分路径;N为平面波的个数;Ri(p)为平面波在多个界面上的透射和反射系数综合;j是虚数单位,ω是角频率;p1i、q1i分别为第i个界面反射回第1层的平面波的水平和垂直慢度;p、q分别为点源处平面波的水平和垂直慢度,x和z表示位置坐标。
2.根据权利要求1所述的倾斜层状粘弹性介质中波场正演模拟方法,其特征在于,步骤2)中,模拟震源为P波时的一次反射P波,设震源和检波器均在第一层,将震源处入射波逐层向下传播的递推步骤为:
④反射波φ′i,i逐层向上递推到第一层,遇到界面时执行①-③步,但只取透射波;
⑤透射波φi+1继续向下递推,遇到下一个界面时回到步骤①。
3.根据权利要求1所述的倾斜层状粘弹性介质中波场正演模拟方法,其特征在于,步骤2)中,模拟其它类型的波采用与模拟震源为P波时的一次反射P波的相同方法和步骤进行。
4.根据权利要求1所述的倾斜层状粘弹性介质中波场正演模拟方法,其特征在于,步骤3)中,式中积分方法采用Frazer和Gettrust提出的GFM积分法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310044384.2A CN103149586B (zh) | 2013-02-04 | 2013-02-04 | 一种倾斜层状粘弹性介质中波场正演模拟方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310044384.2A CN103149586B (zh) | 2013-02-04 | 2013-02-04 | 一种倾斜层状粘弹性介质中波场正演模拟方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN103149586A true CN103149586A (zh) | 2013-06-12 |
CN103149586B CN103149586B (zh) | 2016-01-20 |
Family
ID=48547768
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201310044384.2A Active CN103149586B (zh) | 2013-02-04 | 2013-02-04 | 一种倾斜层状粘弹性介质中波场正演模拟方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN103149586B (zh) |
Cited By (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104570072A (zh) * | 2013-10-16 | 2015-04-29 | 中国石油化工股份有限公司 | 一种粘弹性介质中的球面pp波反射系数建模方法 |
CN104732093A (zh) * | 2015-03-30 | 2015-06-24 | 西安交通大学 | 一种基于弥散黏滞性波动方程的fct-fdm正演模拟方法 |
CN104950332A (zh) * | 2015-06-18 | 2015-09-30 | 河海大学 | 一种弹性多层介质中平面波反射系数的计算方法 |
CN106662664A (zh) * | 2014-06-17 | 2017-05-10 | 埃克森美孚上游研究公司 | 快速粘声波和粘弹性全波场反演 |
CN108445534A (zh) * | 2018-03-12 | 2018-08-24 | 中海石油(中国)有限公司 | 一种弥散粘滞介质中波场模拟方法 |
CN108802819A (zh) * | 2018-06-26 | 2018-11-13 | 西安交通大学 | 一种深度均匀采样梯形网格有限差分地震波场模拟方法 |
CN109164488A (zh) * | 2018-10-10 | 2019-01-08 | 西安交通大学 | 一种梯形网格有限差分地震波场模拟方法 |
CN110427680A (zh) * | 2019-07-26 | 2019-11-08 | 中国水利水电科学研究院 | 获取地震波斜入射作用下边坡地震动力放大效应的方法 |
CN110515125A (zh) * | 2019-08-30 | 2019-11-29 | 长安大学 | 一种岩体中柱面波引起的地表振动速度确定方法 |
CN110954959A (zh) * | 2018-09-27 | 2020-04-03 | 中国石油化工股份有限公司 | 球面透射波特性分析方法及计算机可读存储介质 |
CN110988991A (zh) * | 2019-12-16 | 2020-04-10 | 中国石油大学(北京) | 一种弹性参数反演方法、装置及*** |
CN112505751A (zh) * | 2019-09-16 | 2021-03-16 | 中国石油化工股份有限公司 | 球面波ps反射系数计算方法及*** |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20060037401A1 (en) * | 2003-08-12 | 2006-02-23 | Northwestern University | Scanning near field ultrasound holography |
CN101206266A (zh) * | 2007-12-10 | 2008-06-25 | 吴以雄 | 地球物理勘探方法及设备 |
CN102096100A (zh) * | 2010-11-24 | 2011-06-15 | 中国石油天然气集团公司 | 测井曲线与地震记录全波匹配方法及装置 |
CN102590863A (zh) * | 2012-01-29 | 2012-07-18 | 中国石油天然气股份有限公司 | 确定反射系数频散、反射能量及识别气层的方法及装置 |
-
2013
- 2013-02-04 CN CN201310044384.2A patent/CN103149586B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20060037401A1 (en) * | 2003-08-12 | 2006-02-23 | Northwestern University | Scanning near field ultrasound holography |
CN101206266A (zh) * | 2007-12-10 | 2008-06-25 | 吴以雄 | 地球物理勘探方法及设备 |
CN102096100A (zh) * | 2010-11-24 | 2011-06-15 | 中国石油天然气集团公司 | 测井曲线与地震记录全波匹配方法及装置 |
CN102590863A (zh) * | 2012-01-29 | 2012-07-18 | 中国石油天然气股份有限公司 | 确定反射系数频散、反射能量及识别气层的方法及装置 |
Non-Patent Citations (2)
Title |
---|
汪超等: "倾斜层状粘弹性介质中的地展记录合成", 《中国地球物理学会第二十四届年会论文集》, 30 April 2008 (2008-04-30), pages 416 * |
牛滨华等: "地震波长正演模拟中的有限单元方法", 《长春地质学院学报》, vol. 19, no. 2, 30 April 1989 (1989-04-30), pages 213 - 220 * |
Cited By (16)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104570072A (zh) * | 2013-10-16 | 2015-04-29 | 中国石油化工股份有限公司 | 一种粘弹性介质中的球面pp波反射系数建模方法 |
CN104570072B (zh) * | 2013-10-16 | 2017-10-27 | 中国石油化工股份有限公司 | 一种粘弹性介质中的球面pp波反射系数建模方法 |
CN106662664A (zh) * | 2014-06-17 | 2017-05-10 | 埃克森美孚上游研究公司 | 快速粘声波和粘弹性全波场反演 |
CN104732093A (zh) * | 2015-03-30 | 2015-06-24 | 西安交通大学 | 一种基于弥散黏滞性波动方程的fct-fdm正演模拟方法 |
CN104732093B (zh) * | 2015-03-30 | 2018-04-17 | 西安交通大学 | 一种基于弥散黏滞性波动方程的fct‑fdm正演模拟方法 |
CN104950332A (zh) * | 2015-06-18 | 2015-09-30 | 河海大学 | 一种弹性多层介质中平面波反射系数的计算方法 |
CN104950332B (zh) * | 2015-06-18 | 2018-02-02 | 河海大学 | 一种弹性多层介质中平面波反射系数的计算方法 |
CN108445534A (zh) * | 2018-03-12 | 2018-08-24 | 中海石油(中国)有限公司 | 一种弥散粘滞介质中波场模拟方法 |
CN108802819A (zh) * | 2018-06-26 | 2018-11-13 | 西安交通大学 | 一种深度均匀采样梯形网格有限差分地震波场模拟方法 |
CN110954959A (zh) * | 2018-09-27 | 2020-04-03 | 中国石油化工股份有限公司 | 球面透射波特性分析方法及计算机可读存储介质 |
CN110954959B (zh) * | 2018-09-27 | 2021-11-05 | 中国石油化工股份有限公司 | 球面透射波特性分析方法及计算机可读存储介质 |
CN109164488A (zh) * | 2018-10-10 | 2019-01-08 | 西安交通大学 | 一种梯形网格有限差分地震波场模拟方法 |
CN110427680A (zh) * | 2019-07-26 | 2019-11-08 | 中国水利水电科学研究院 | 获取地震波斜入射作用下边坡地震动力放大效应的方法 |
CN110515125A (zh) * | 2019-08-30 | 2019-11-29 | 长安大学 | 一种岩体中柱面波引起的地表振动速度确定方法 |
CN112505751A (zh) * | 2019-09-16 | 2021-03-16 | 中国石油化工股份有限公司 | 球面波ps反射系数计算方法及*** |
CN110988991A (zh) * | 2019-12-16 | 2020-04-10 | 中国石油大学(北京) | 一种弹性参数反演方法、装置及*** |
Also Published As
Publication number | Publication date |
---|---|
CN103149586B (zh) | 2016-01-20 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN103149586B (zh) | 一种倾斜层状粘弹性介质中波场正演模拟方法 | |
Operto et al. | On the role of density and attenuation in three-dimensional multiparameter viscoacoustic VTI frequency-domain FWI: An OBC case study from the North Sea | |
Wang et al. | Vector-based elastic reverse time migration | |
EP3259619B1 (en) | Method and system of processing seismic data by providing surface aperture common image gathers | |
Du et al. | Polarity reversal correction for elastic reverse time migration | |
EP2715405B1 (en) | Method of processing seismic data by providing surface offset common image gathers | |
EP2689273B1 (en) | System and method for seismic data modeling and migration | |
US20170299745A1 (en) | Prestack egs migration method for seismic wave multi-component data | |
CN102053261B (zh) | 一种地震数据处理方法 | |
CN105182408A (zh) | 一种合成地震记录的制作方法和装置 | |
CN102890290A (zh) | 一种起伏地表条件下的叠前深度偏移方法 | |
GB2499096A (en) | Simultaneous joint estimation of P-P and P-S residual statics | |
CN104808243A (zh) | 一种叠前地震贝叶斯反演方法和装置 | |
CN107817526A (zh) | 叠前地震道集分段式振幅能量补偿方法及*** | |
CN104237937A (zh) | 叠前地震反演方法及其*** | |
CN102540253B (zh) | 一种陡倾角地层及裂缝地震偏移方法和装置 | |
CN102053262B (zh) | 获得地震转换波的方位速度的方法及处理地震数据的方法 | |
CN104345343A (zh) | 一种复杂海底相关的层间多次波预测方法 | |
CN105242313A (zh) | 一种弹性波逆时偏移极性反转的校正方法及*** | |
CN106353798A (zh) | 多分量联合高斯束叠前逆时偏移成像方法 | |
CN107340537A (zh) | 一种p-sv转换波叠前逆时深度偏移的方法 | |
CN105572735A (zh) | 一种提高叠前深度偏移成像精度的方法及装置 | |
CN106662665A (zh) | 用于更快速的交错网格处理的重新排序的插值和卷积 | |
CN105527648B (zh) | 用于各向异性参数反演的敏感度矩阵的计算方法及*** | |
Cormier et al. | Calculation of strong ground motion due to an extended earthquake source in a laterally varying structure |
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 |