CN113917526A - 基于非分裂完全匹配层吸收边界的正演模拟方法 - Google Patents
基于非分裂完全匹配层吸收边界的正演模拟方法 Download PDFInfo
- Publication number
- CN113917526A CN113917526A CN202010661155.5A CN202010661155A CN113917526A CN 113917526 A CN113917526 A CN 113917526A CN 202010661155 A CN202010661155 A CN 202010661155A CN 113917526 A CN113917526 A CN 113917526A
- Authority
- CN
- China
- Prior art keywords
- matching layer
- absorption boundary
- equation
- absorption
- order
- 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
- 238000010521 absorption reaction Methods 0.000 title claims abstract description 95
- 238000000034 method Methods 0.000 title claims abstract description 33
- 230000000694 effects Effects 0.000 claims abstract description 14
- 230000006870 function Effects 0.000 claims description 6
- 238000006073 displacement reaction Methods 0.000 claims description 3
- 238000001914 filtration Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 238000013507 mapping Methods 0.000 claims description 2
- 238000004364 calculation method Methods 0.000 abstract description 6
- 238000009304 pastoral farming Methods 0.000 abstract description 3
- 238000004519 manufacturing process Methods 0.000 abstract description 2
- 238000004088 simulation Methods 0.000 description 11
- 238000005070 sampling Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 3
- 150000001875 compounds Chemical class 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/66—Subsurface modeling
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明提供一种基于非分裂完全匹配层吸收边界的正演模拟方法,该基于非分裂完全匹配层吸收边界的正演模拟方法包括:步骤1,建立二阶弹性波动方程;步骤2,利用迭代法代替卷积运算的思想构建适用于二阶弹性波方程的非分裂完全匹配层吸收边界;步骤3,实现非分裂完全匹配层吸收边界二阶精度弹性波正演模拟。该基于非分裂完全匹配层吸收边界的正演模拟方法与常规的分裂完全匹配层吸收边界相比,非分裂完全匹配层吸收边界的吸收效果更好,能够吸收高角度入射波和掠射波,计算速度更快、所占内存更小,降低了算法对计算机内存的需求,更易满足实际生产需求。
Description
技术领域
本发明涉及地球物理勘查技术领域,特别是涉及到一种基于非分裂卷积完全匹配层弹性波数值模拟方法。
背景技术
波动方程法是应对日益复杂的地质构造地震波传播数值模拟的最重要方法。由于无法构建无限区域的模型,使用波动方程法进行地震波正演数值模拟时,一般需要用有限的区域模拟无限区域波传播过程,即添加人工边界。为了消除地震波传播到人工边界引起的虚假反射,就需要设置吸收边界条件,使得地震波可以无反射的穿过人工边界。因此,发展形成有效而稳定的吸收边界条件一直是地震波数值模拟的一项重要研究内容。完全匹配层吸收边界是在模型区域外人为增加了一层吸收介质,当波传播到该介质中时,波的能量将被介质以指数衰减的速度迅速吸收,直至消失,而不会反弹回模型区域。理论上,通过合理的设置该吸收层的吸收系数,无论什么频率的波以任意角度从模型区域进入吸收介质层时,都不产生反射。与传统的吸收边界条件相比,完全匹配层具有更优秀的吸收效果。目前应用较为广泛的完全匹配层吸收边界为分裂完全匹配层(SPML)吸收边界。但在大角度入射的情况下,SPML吸收边界的吸收效果并不理想,虚假反射将会更加明显,对薄层介质、震源靠近模型离散的网格边缘或者大偏移距模型的数值模拟效果不理想,此外,SPML边界无法吸收掠射波,并且正演时对内存需求较大。在分裂完全匹配层(SPML)吸收边界的基础上发展起来的卷积完全匹配层(convolution perfectly matched layer,CPML)吸收边界能够吸收大角度入射波和掠射波,但目前构建的卷积完全匹配层吸收边界多适用于一阶波动方程,难以直接应用于二阶波动方程。
为解决以上技术问题,本发明提出一种基于非分裂完全匹配层吸收边界正演模拟方法。
发明内容
本发明的目的是提供一种适用于二阶弹性波方程的非分裂卷积完全匹配层吸收边界条件的数值模拟方法。
本发明的目的可通过如下技术措施来实现:基于非分裂完全匹配层吸收边界的正演模拟方法,该基于非分裂完全匹配层吸收边界的正演模拟方法包括:步骤1,建立二阶弹性波动方程;步骤2,利用迭代法代替卷积运算的思想构建适用于二阶弹性波方程的非分裂完全匹配层吸收边界;步骤3,实现非分裂完全匹配层吸收边界二阶精度弹性波正演模拟。
本发明的目的还可通过如下技术措施来实现:
在步骤1中,通过柯西方程、几何方程以及纳维尔方程推导出二阶弹性波波动方程:
其中λ和μ为介质的拉梅常数;ρ表示介质密度;u表示地震波场;ux和uz分别为水平方向和垂直方向的位移波场,t表示地震波传播时间,x,z分别表示水平方向和垂直方向。
在步骤2中,采用复坐标变换,忽略衰减因子空变影响计算出二阶弹性波方程的空间二阶偏导算子:
其中:δ(t)和H(t)分别代表狄拉克函数以及单位阶跃函数,t为时间,dx、dz为衰减透射波的衰减参数;αx、αz为滤波的频率参数;为x方向在复频域的拉伸变量,x与z方向在复频域的拉伸变量,kx,kz为衰减渐消波的尺度因子。
在步骤2中,利用迭代法代替卷积运算的思想,将二阶偏导算子中的时域卷积运算转换成迭代运算计算出时间域基于非分裂完全匹配层吸收边界的二阶弹性波方程。
在步骤3中,由连续方程的时间二阶、空间十阶有限差分计算基于非分裂完全匹配层吸收边界二阶弹性波方程的离散方程。
在步骤3中,对均匀介质模型和复杂层状模型分别开展分裂完全匹配层和非分裂完全匹配层吸收边界二阶精度弹性波方程正演模拟,通过成图验证非分裂完全匹配层吸收边界具有更好的吸收效果、更快的计算速度和对计算机内存较低的需求。
本发明的基于非分裂完全匹配层吸收边界的正演模拟方法,由连续方程的时间二阶、空间十阶有限差分计算基于非分裂完全匹配层吸收边界二阶弹性波方程的离散方程,基于卷积完全匹配层(CPML)吸收边界利用迭代运算取代卷积运算的思想推导出非分裂完全匹配层吸收边界(NCPML,new convolution perfectly matched layer,又称新卷积完全匹配层),非分裂完全匹配层吸收边界能够吸收大角度的入射波和掠射波,能够直接对二阶弹性波方程做数值模拟,减小内存使用空间,提高计算效率,更适用于实际生产需求。
附图说明
图1是本发明的基于非分裂完全匹配层吸收边界的正演模拟方法的一具体实施例流程图;
图2是本发明的具体实施例中均匀介质模型纵波源激发分裂完全匹配层(SPML)和非分裂完全匹配层(NCPML)两种不同吸收边界条件下二阶弹性波正演的水平分量波场图;(a)为SPML边界700ms时刻的波场快照;(b)为NCPML边界700ms时刻的波场快照;
图3是本发明的具体实施例中不同吸收边界条件下x=900m处700ms时刻水平分量波形对比图;
图4是本发明的具体实施例中层状速度模型的示意图;
图5是本发明的具体实施例中层状模型t=350ms时两种不同吸收边界条件下水平分量波场快照对比图;(a)层状模型SPML吸收边界条件下水平分量波场快照;(b)层状模型NCPML吸收边界条件下水平分量波场快照;
图6是本发明的具体实施例中第150道SPML吸收边界和NCPML吸收边界条件下的水平分量波形曲线与理论水平分量波形曲线对比图;
图7是本发明的具体实施例中图6在深度为400m至800m处局部放大水平分量波形曲线对比图;
图8是本发明的具体实施例中图6SPML吸收边界和NCPML吸收边界波场水平分量的数值解与理论值的误差曲线对比图;
图9是本发明的具体实施例中NCPML吸收边界的层状模型对应的地震记录的水平分量和垂直分量图,(a)水平分量地震记录,(b)垂直分量。
具体实施方式
为使本发明的上述和其他目的、特征和优点能更明显易懂,下文特举出较佳实施例,并配合附图所示,作详细说明如下。
图1所示为本发明的基于非分裂完全匹配层弹性波数值模拟方法的流程图。
第一步:建立二阶弹性波波动方程
通过柯西方程、几何方程以及纳维尔方程推导出二阶弹性波波动方程:
其中λ和μ为介质的拉梅常数;ρ表示介质密度;ux和uz分别为水平方向和垂直方向的位移波场;
第二步:利用迭代法代替卷积运算的思想构建适用于二阶弹性波方程的非分裂完全匹配层吸收边界
CPML吸收边界x方向的复拉伸变量可以表示为:
采用以上复坐标变换公式(2),将二阶弹性波方程表示为:
忽略衰减因子空变影响计算出方程(3)中空间二阶偏导算子:
其中:δ(t)和H(t)分别代表狄拉克函数以及单位阶跃函数,t为时间,dx、fz为衰减透射波的衰减参数;αx、αz为滤波的频率参数;为x方向在复频域的拉伸变量,为x与z方向在复频域的拉伸变量,kx,kz为衰减渐消波的尺度因子。
利用迭代法代替卷积运算的思想,将时域卷积运算转换成迭代运算并计算得到时间域基于NCPML的二阶弹性波方程:
1)当dx/kx+αx≠dz/kz+αz时,基于NCPML时间域的控制方程为:
2)当dx/kx+αx=dz/kz+αz时,基于NCPML时间域的控制方程为:
方程(5)(6)中:是将卷积运算转换成迭代运算的迭代项;bx、ax为转换过程中产生的两个系数, 类似。w为波场垂直分量,u为波场水平分量,t为采样时间。vp表示纵波速度、vs表示横波速度。其余符号含义与公式(4)相同。
第三步:实现非分裂完全匹配层吸收边界二阶精度弹性波正演模拟为实现正演模拟,对连续方程进行时间二阶、空间十阶差分计算,差分方程格式如下:
时间二阶精度差分格式:
空间十阶精度差分格式:
离散后的弹性波方程为:
1)当dx/kx+αx≠dz/kz+αz时,基于NCPML时间域的控制方程的离散形式为:
2)当dx/kx+αx=dz/kz+αz时,基于NCPML时间域的控制方程离散形式为:
正演模拟使用的子波主频为20Hz,时间采样间隔为1ms,采样点数为1000。均匀模型震源位于网格点(100,100),层状模型震源位于(1000m,0m)处。
通过数值模拟对比NCPML边界和SPML边界在吸收效果、运行速度以及内存占用方面的不同,并用复杂的层状模型测试NCPML边界对复杂模型的适定性。首先通过层状模型对比分析NCPML吸收边界与SPML吸收边界的吸收效果以及计算效率;然后通过层状模型验证所推导的二阶弹性波NCPML吸收边界的稳定性。
图2是本发明的一具体实施例中均匀介质模型纵波源激发的水平分量波场图,模型大小为200*200,Vp=3000m/s,Vs=2000m/s,网格间距为10m。(a)为SPML边界700ms时刻的波场快照;(b)为NCPML边界700ms时刻的波场快照;对比图2(a)和图2(b),直观上两种吸收边界对人工边界产生的反射波都具有良好的吸收效果,在边界处没有明显的虚假反射。定量分析两种吸收边界的吸收效果,图3是图2中虚线处(在x=900m的垂线上)用两种PML吸收边界数值模拟得出的波场值与理论值的对比图,其中理论值是通过增加计算区域得到的,其计算结果等同于对边界处的透射波完全吸收。通过图2和图3可以看出与SPML吸收边界相比,NCPML吸收边界更接近理论值,吸收效果更好。
图4为存在多个速度界面的复杂层状速度模型。在层状模型数值模拟中,图5(a)为层状模型SPML吸收边界t=350ms时水平分量波场快照,(b)为层状模型NCPML吸收边界t=350ms时水平分量波场快照;图6是第150道SPML吸收边界和NCPML吸收边界的水平分量波形曲线与理论水平分量波形曲线对比;由图5和图6可以看出NCPML吸收边界与SPML吸收边界的吸收效果都较好。为进一步分析两者的区别,对图6所示曲线进行局部放大显示,并分别计算两者与理论值的误差。图7是图6在深度为400m至800m处水平分量波形曲线;图8是用图6中SPML吸收边界和NCPML吸收边界水平分量波场的数值解与理论值的误差曲线;通过图7和图8可以看出NCPML吸收边界的数值更加接近理论值,其吸收效果更好。
表1 SPML边界与NCPML边界计算效率与占用内存数据对比
表1为均匀介质模型SPML边界与NCPML边界计算效率与占用内存的数据对比,通过表1的数据对比可知NCPML吸收边界相较于SPML吸收边界计算速度更快,所占内存更小,更适用于生产实践。图9是基于NCPML吸收边界的层状模型对应的地震记录,图9(a)水平分量,图9(b)垂直分量;由图9可以看出,得到的地震记录都没有明显的虚假反射波,证明了本方法对于复杂模型的数值模拟稳定性良好。
Claims (6)
1.基于非分裂完全匹配层吸收边界的正演模拟方法,其特征在于,该基于非分裂完全匹配层吸收边界的正演模拟方法包括:
步骤1,建立二阶弹性波动方程;
步骤2,利用迭代法代替卷积运算的思想构建适用于二阶弹性波方程的非分裂完全匹配层吸收边界;
步骤3,实现非分裂完全匹配层吸收边界二阶精度弹性波正演模拟。
4.根据权利要求3所述的基于非分裂完全匹配层吸收边界的正演模拟方法,其特征在于,在步骤2中,利用迭代法代替卷积运算的思想,将二阶偏导算子中的时域卷积运算转换成迭代运算计算出时间域基于非分裂完全匹配层吸收边界的二阶弹性波方程。
5.根据权利要求4所述的基于非分裂完全匹配层吸收边界的正演模拟方法,其特征在于,在步骤3中,由连续方程的时间二阶、空间十阶有限差分计算基于非分裂完全匹配层吸收边界二阶弹性波方程的离散方程。
6.根据权利要求5所述的基于非分裂完全匹配层吸收边界的正演模拟方法,其特征在于,在步骤3中,对均匀介质模型和复杂层状模型分别开展分裂完全匹配层和非分裂完全匹配层吸收边界二阶精度弹性波方程正演模拟,通过成图验证非分裂完全匹配层吸收边界具有更好的吸收效果、更快的计算速度和对计算机内存较低的需求。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010661155.5A CN113917526B (zh) | 2020-07-10 | 2020-07-10 | 基于非分裂完全匹配层吸收边界的正演模拟方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010661155.5A CN113917526B (zh) | 2020-07-10 | 2020-07-10 | 基于非分裂完全匹配层吸收边界的正演模拟方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113917526A true CN113917526A (zh) | 2022-01-11 |
CN113917526B CN113917526B (zh) | 2023-07-28 |
Family
ID=79231998
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010661155.5A Active CN113917526B (zh) | 2020-07-10 | 2020-07-10 | 基于非分裂完全匹配层吸收边界的正演模拟方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113917526B (zh) |
Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20080046185A1 (en) * | 2006-04-10 | 2008-02-21 | Tim Niebauer | Method and apparatus for processing an under-sampled chirped sinusoidal waveform using a complex-heterodyne |
CN102253415A (zh) * | 2011-04-19 | 2011-11-23 | 中国石油大学(华东) | 基于裂缝等效介质模型的地震响应模式建立方法 |
CN104459791A (zh) * | 2014-12-15 | 2015-03-25 | 中国石油大学(华东) | 一种基于波动方程的小尺度大模型正演模拟方法 |
US20160223697A1 (en) * | 2015-02-04 | 2016-08-04 | Tetyana Vdovina | Poynting Vector Minimal Reflection Boundary Conditions |
WO2016130208A1 (en) * | 2015-02-13 | 2016-08-18 | Exxonmobil Upstream Research Company | Efficient and stable absorbing boundary condition in finite-difference calculations |
CN106777472A (zh) * | 2016-11-16 | 2017-05-31 | 西安理工大学 | 基于拉盖尔多项式的减少分裂误差的完全匹配层实现方法 |
CN107479092A (zh) * | 2017-08-17 | 2017-12-15 | 电子科技大学 | 一种基于方向导数的频率域高阶声波方程正演模拟方法 |
CN107944214A (zh) * | 2017-11-27 | 2018-04-20 | 河北工业大学 | 笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法 |
CN111208563A (zh) * | 2020-02-18 | 2020-05-29 | 吉林大学 | 一种非分裂完全匹配层吸收边界方法 |
-
2020
- 2020-07-10 CN CN202010661155.5A patent/CN113917526B/zh active Active
Patent Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20080046185A1 (en) * | 2006-04-10 | 2008-02-21 | Tim Niebauer | Method and apparatus for processing an under-sampled chirped sinusoidal waveform using a complex-heterodyne |
CN102253415A (zh) * | 2011-04-19 | 2011-11-23 | 中国石油大学(华东) | 基于裂缝等效介质模型的地震响应模式建立方法 |
CN104459791A (zh) * | 2014-12-15 | 2015-03-25 | 中国石油大学(华东) | 一种基于波动方程的小尺度大模型正演模拟方法 |
US20160223697A1 (en) * | 2015-02-04 | 2016-08-04 | Tetyana Vdovina | Poynting Vector Minimal Reflection Boundary Conditions |
WO2016130208A1 (en) * | 2015-02-13 | 2016-08-18 | Exxonmobil Upstream Research Company | Efficient and stable absorbing boundary condition in finite-difference calculations |
CN106777472A (zh) * | 2016-11-16 | 2017-05-31 | 西安理工大学 | 基于拉盖尔多项式的减少分裂误差的完全匹配层实现方法 |
CN107479092A (zh) * | 2017-08-17 | 2017-12-15 | 电子科技大学 | 一种基于方向导数的频率域高阶声波方程正演模拟方法 |
CN107944214A (zh) * | 2017-11-27 | 2018-04-20 | 河北工业大学 | 笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法 |
CN111208563A (zh) * | 2020-02-18 | 2020-05-29 | 吉林大学 | 一种非分裂完全匹配层吸收边界方法 |
Also Published As
Publication number | Publication date |
---|---|
CN113917526B (zh) | 2023-07-28 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Dragna et al. | A generalized recursive convolution method for time-domain propagation in porous media | |
Fuyuki et al. | Finite difference analysis of Rayleigh wave scattering at a trench | |
CN103630933B (zh) | 基于非线性优化的时空域交错网格有限差分方法和装置 | |
CN104122585B (zh) | 基于弹性波场矢量分解与低秩分解的地震正演模拟方法 | |
CN105467444B (zh) | 一种弹性波全波形反演方法及装置 | |
CN112327358B (zh) | 一种粘滞性介质中声波地震数据正演模拟方法 | |
CA2947410A1 (en) | Fast viscoacoustic and viscoelastic full-wavefield inversion | |
AU2007201220A1 (en) | Method of evaluating the interaction between a wavefield and a solid body | |
CN107894612A (zh) | 一种q吸收衰减补偿的声波阻抗反演方法及系统 | |
CN105911584B (zh) | 一种隐式交错网格有限差分弹性波数值模拟方法及装置 | |
MX2011003850A (es) | Estimado de señal de dominio de imagen a interferencia. | |
CN109946742B (zh) | 一种TTI介质中纯qP波地震数据模拟方法 | |
Qin et al. | The implementation of an improved NPML absorbing boundary condition in elastic wave modeling | |
CN104965222A (zh) | 三维纵波阻抗全波形反演方法及装置 | |
CN115600373A (zh) | 黏滞各向异性介质qP波模拟方法、系统、设备及应用 | |
CN114895351A (zh) | 模拟地震波在任意不连续界面处传播的介质模型化的方法及装置 | |
Ha et al. | 3D Laplace-domain waveform inversion using a low-frequency time-domain modeling algorithm | |
CN104732093B (zh) | 一种基于弥散黏滞性波动方程的fct‑fdm正演模拟方法 | |
CN113917526A (zh) | 基于非分裂完全匹配层吸收边界的正演模拟方法 | |
Nakagawa et al. | Three-dimensional elastic wave scattering by a layer containing vertical periodic fractures | |
CN115201896A (zh) | 吸收衰减介质逆时偏移方法、装置、成像方法、介质 | |
WO2022254247A1 (en) | Laplace-fourier 1.5d forward modeling using an adaptive sampling technique | |
Zyserman et al. | Analysis of the numerical dispersion of waves in saturated poroelastic media | |
MX2011003852A (es) | Atributos de procesamiento de imagen con inversion de tiempo. | |
CN112649848A (zh) | 利用波动方程求解地震波阻抗的方法和装置 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |