CN114624766B - 基于行波分离的弹性波最小二乘逆时偏移梯度求取方法 - Google Patents
基于行波分离的弹性波最小二乘逆时偏移梯度求取方法 Download PDFInfo
- Publication number
- CN114624766B CN114624766B CN202210525461.5A CN202210525461A CN114624766B CN 114624766 B CN114624766 B CN 114624766B CN 202210525461 A CN202210525461 A CN 202210525461A CN 114624766 B CN114624766 B CN 114624766B
- Authority
- CN
- China
- Prior art keywords
- wave
- traveling wave
- field
- wave field
- formula
- 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.)
- Active
Links
- 238000013508 migration Methods 0.000 title claims abstract description 32
- 230000005012 migration Effects 0.000 title claims abstract description 32
- 238000000034 method Methods 0.000 title claims abstract description 25
- 238000000926 separation method Methods 0.000 title claims abstract description 24
- 239000013598 vector Substances 0.000 claims description 19
- 238000003384 imaging method Methods 0.000 claims description 9
- 238000000354 decomposition reaction Methods 0.000 claims description 6
- 230000005540 biological transmission Effects 0.000 claims description 5
- 230000001902 propagating effect Effects 0.000 claims description 3
- 238000011144 upstream manufacturing Methods 0.000 claims 2
- 238000002474 experimental method Methods 0.000 abstract description 5
- 238000010586 diagram Methods 0.000 description 7
- 238000004364 calculation method Methods 0.000 description 3
- 238000001514 detection method Methods 0.000 description 3
- 238000004422 calculation algorithm Methods 0.000 description 2
- 230000002194 synthesizing effect Effects 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000005422 blasting Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000003607 modifier Substances 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 238000004088 simulation Methods 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
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/301—Analysis for determining seismic cross-sections or geostructures
-
- 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/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- 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/64—Geostructures, e.g. in 3D data cubes
- G01V2210/642—Faults
-
- 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
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/70—Other details related to processing
- G01V2210/74—Visualisation of seismic data
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Remote Sensing (AREA)
- General Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Geophysics (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Mathematical Physics (AREA)
- Geology (AREA)
- Environmental & Geological Engineering (AREA)
- Acoustics & Sound (AREA)
- Mathematical Optimization (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Databases & Information Systems (AREA)
- Algebra (AREA)
- Operations Research (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
本发明涉及一种基于行波分离的弹性波最小二乘逆时偏移梯度求取方法,属于地震速度建模领域,首先,进行线性化正演,得到正传波场,通过伴随方程得到反传波场,然后利用Poynting矢量行波分离方法将正传波场分为正传波场上行波、下行波、左行波和右行波,将反传波场分为反传波场上行波、下行波、左行波和右行波,最后将分离后的正、反传波场带入梯度公式,得到梯度剖面。模型实验结果表明本方法得到的梯度剖面可有效去除低波数噪声,提高梯度剖面的精度,从而加快弹性波方程最小二乘逆时偏移的收敛速度。
Description
技术领域
本发明属于地震速度建模领域,具体涉及一种基于Poynting矢量行波分离的弹性波最小二乘逆时偏移的梯度求取方法。
背景技术
最小二乘逆时偏移作为一种高精度的地震速度建模算法,其每次迭代过程中,首先基于初始反射系数模型和背景速度,采用双程波方程有限差分数值模拟技术进行线性化正演得到模拟地震记录,然后将模拟地震记录与实测地震记录的差值放入L2范数中建立目标泛函,基于梯度类方法进行迭代反演使得目标泛函最小,从而得到模型修改量完成一次迭代。整个最小二乘逆时偏移通常需要几十次甚至上百次的迭代反演。
最小二乘逆时偏移可以将逆时偏移的成像问题作为最小二乘反演问题求解,通过迭代算法求解一个与观测地震数据最佳匹配的成像结果,相对于逆时偏移而言,其具有压制偏移噪音、改善深部成像、提高同相轴的连续性和均衡性以及改善振幅保真性的优势,因此,最小二乘逆时偏移在未来的实际生产中有着巨大的应用潜力。
常规最小二乘逆时偏移算法求取梯度剖面时,并未对正传波场和反传波场进行方向上的区分,求取梯度时不同方向波场的互相关计算会产生大量的低波数噪声,其会严重降低反演收敛速度和成像精度,限制了该方法在实际数据处理中的应用。
发明内容
本发明要解决的技术问题在于提供一种基于行波分离的弹性波最小二乘逆时偏移梯度求取方法。首先,进行线性化正演,得到正传波场,通过伴随方程得到反传波场,然后利用光流矢量行波分离方法将正传波场分为正传波场上行波、下行波、左行波和右行波,将反传波场分为反传波场上行波、下行波、左行波和右行波,最后将正传波场上行波和反传波场下行波、正传波场下行波和反传波场上行波、正传波场左行波和反传波场右行波、正传波场右行波和反传波场左行波分别进行互相关成像并相加,得到梯度剖面。模型实验结果表明本方法得到的梯度剖面可有效去除低波数噪声,提高梯度剖面的精度,从而加快最小二乘逆时偏移的收敛速度。
本发明采取以下技术方案:
一种基于行波分离的弹性波最小二乘逆时偏移梯度求取方法,所述方法具体包括以下步骤:
(1)给定背景横波速度V x (x, z)、背景纵波速度V z (x, z)、背景密度D(x, z)和初始反射系数模型M(x, z),其中x、z表示空间位置坐标,x=1, 2, 3, … , Nx, z=1, 2, 3, …, Nz,Nx、Nz分别代表模型横向和纵向网格点总数;给定雷克子波W(t),t表示时间,基于背景速度、背景密度和初始反射系数模型,应用弹性波方程进行线性化正演,得到正传波场τ xx (x, z, t)、τ xz (x, z, t)、τ zz (x, z, t)、v x (x, z, t)、v z (x, z, t)和合成地震记录Cal(x,z, t);
(2)获得实际地震记录Obs(x, z, t),将其与合成地震记录的差值作为逆时扰动带入伴随方程得到反传波场Φ xx (x, z, t)、Φ xz (x, z, t)、Φ zz (x, z, t)、ϕ x (x, z, t)、ϕ z (x, z, t);
(3)利用Poynting矢量行波分离公式(1),分别带入正传波场τ xx (x, z, t)、τ xz (x,z, t)、τ zz (x, z, t)、v x (x, z, t)、v z (x, z, t)和反传波场Φ xx (x, z, t)、Φ xz (x, z, t)、Φ zz (x, z, t)、ϕ x (x, z, t)、ϕ z (x, z, t),计算其Poynting矢量的水平分量和垂直分量和,以正传波场为例,所述的公式(1)为:;
(4)利用Poynting矢量方向指示公式(2)将正传波场v x (x, z, t)分解为正传波场上行波、下行波、左行波和右行波,将正传波场v z (x, z, t)分解为正传波场上行波、下行波、左行波和右行波,将反传波场ϕ x (x, z, t)分解为反传波场上行波、下行波、左行波和右行波,将反传波场ϕ z (x, z, t)分解为反传波场上行波、下行波、左行波和右行波;以正传波场v x (x, z, t)为例,所述的公式(2)为:;
弹性波方程最小二乘逆时偏移每次迭代中的梯度剖面求取均需要重复以上5个步骤。
本发明与现有技术相比的有益效果:
本发明提出的一种基于Poynting矢量行波分离的弹性波方程最小二乘逆时偏移梯度求取方法,首先,进行线性化正演,得到正传波场,通过伴随方程得到反传波场,然后利用Poynting矢量行波分离方法将正传波场分为正传波场上行波、下行波、左行波和右行波,将反传波场分为反传波场上行波、下行波、左行波和右行波,最后将分离后的正、反传波场带入梯度公式,得到梯度剖面。模型实验结果表明本方法得到的梯度剖面可有效去除低波数噪声,提高梯度剖面的精度,从而加快弹性波最小二乘逆时偏移的收敛速度,节约时间,进而推动最小二乘逆时偏移的工业化进程。
附图说明
图1为弹性波最小二乘逆时偏移常规梯度求取方法流程图;
图2为一种基于Poynting矢量行波分离的弹性波最小二乘逆时偏移梯度求取方法流程图;
图3为横波速度模型的水平层状模型图;
图4为平滑后横波速度模型的速度模型图;
图5为真实反射系数模型图;
图6为正传波场的波场快照图;
图7为正传波场行波分离后的波场快照图;
图8为反传波场的波场快照图;
图9为反传波场行波分离后的波场快照图;
图10为常规梯度剖面图;
图11为基于Poynting矢量行波分离后得到的梯度剖面图;
图12为常规弹性波最小二乘逆时偏移的第10次迭代反演结果图;
图13为基于Poynting矢量行波分离的弹性波最小二乘逆时偏移第10次迭代反演结果图;
图14为抽道对比图。
具体实施方式
下面通过实施例结合附图来对本发明的技术方案做进一步解释,但本发明的保护范围不受实施例任何形式上的限制。
本发明通过将Poynting矢量公式引入弹性波最小二乘逆时偏移的行波分离中,将正传波场和反传波场进行方向上的区分后再进行互相关成像,(具体流程见图2所示)。实验证明,该方法可获得高质量梯度剖面,减少迭代次数,节约计算时间,从而显著提升最小二乘逆时偏移的效果。
实施例1
本发明以水平层状模型为例(如图3所示)进行阐述,该模型横向长度为2000m,纵向深度为1000m,模型x、z方向的网格步长均为5m。
基于该模型,一共进行1次放炮,将炮点位置设置为(1000,0),设置400个检波点,检波点以5m为间隔均匀的分布在模型的正上方,检波点深度均为0m。
下面详细阐述本发明具体实施方式,本发明方法的具体流程见图2所示。
(1)给定横向长度均为2000m,纵向深度均为1000m的水平模型:横波速度模型V x 、纵波速度模型V z 、密度模型D。将这些模型网格化,模型x、z方向的网格步长均为5m,得到网格大小均为400×200的模型,将其平滑处理得到背景横波速度模型V x (x, z)、背景纵波速度模型V z (x, z)、背景密度模型D(x, z)(如图4),计算真实反射系数模型M(x, z)(如图5),并设反射系数模型初始值为0值,其中x、z表示网格位置坐标(x=1, 2, 3, … ,400,z=1,2, 3, … , 200);给定雷克子波W(t),t表示时间,基于背景速度、背景密度和初始反射系数模型,应用弹性波方程进行线性化正演,得到正传波场τ xx (x, z, t)、τ xz (x, z, t)、τ zz (x, z, t)、v x (x, z, t)、v z (x, z, t)和合成地震记录Cal(x, z, t);
(2)给定实际地震记录Obs(x, z, t),将其与合成地震记录的差值作为逆时扰动带入伴随方程得到反传波场Φ xx (x, z, t)、Φ xz (x, z, t)、Φ zz (x, z, t)、ϕ x (x, z, t)、ϕ z (x, z, t);
(3)利用Poynting矢量行波分离公式(1),分别带入正传波场τ xx (x, z, t)、τ xz (x,z, t)、τ zz (x, z, t)、v x (x, z, t)、v z (x, z, t)和反传波场Φ xx (x, z, t)、Φ xz (x, z,t)、Φ zz (x, z, t)、ϕ x (x, z, t)、ϕ z (x, z, t),计算其Poynting矢量的水平分量和垂直分量和,以正传波场为例,所述的公式(1)为:;反传波场的公式与正传波场公式一样
(4)利用Poynting矢量方向指示公式(2)将正传波场v x (x, z, t)(如图6)分解为正传波场上行波、下行波、左行波和右行波(如图7),将正传波场v z (x, z, t)分解为正传波场上行波、下行波、左行波和右行波,将反传波场ϕ x (x, z, t)(如图8)分解为反传波场上行波、下行波、左行波和右行波(如图9),将反传波场ϕ z (x, z, t)分解为反传波场上行波、下行波、左行波和右行波;所述的公式(2)为(以正传波场v x (x, z, t)为例):;
反传波场的公式与正传波场一样
为了说明本发明所述方法的有效性,这里与常规梯度剖面进行对比。常规梯度剖面的求解方法如图1所示,如图10为常规方法所求的梯度剖面,图11为本发明方法基于Poynting矢量行波分离方法求取的梯度剖面。对比两图可以发现,图11的低频噪声得到消除。图12为常规弹性波最小二乘逆时偏移的第10次迭代反演结果,图13为本发明方法基于Poynting矢量行波分离的最小二乘逆时偏移第10次迭代反演结果,对图12、图13以及真实反射系数模型图5的水平方向1000m位置处进行抽道,抽道图如图14所示,可以看出,本发明方法基于Poynting矢量行波分离的最小二乘逆时偏移结果与真实反射系数吻合度更高。模型实验结果表明,本方法可以降低梯度的低波数噪声,提高梯度的精度,将求取的高精度的梯度带入弹性波最小二乘逆时偏移的运算中,可以加快剖面收敛,节约计算时间,从而推动最小二乘逆时偏移的工业化进程。
Claims (1)
1.一种基于行波分离的弹性波最小二乘逆时偏移梯度求取方法,其特征在于所述方法具体包括以下步骤:
(1)给定背景横波速度Vx(x,z)、背景纵波速度Vz(x,z)、背景密度D(x,z)和初始反射系数模型M(x,z),其中x、z表示空间位置坐标,x=1,2,3,…,Nx,z=1,2,3,…,Nz,Nx、Nz分别代表模型横向和纵向网格点总数;给定雷克子波W(t),t表示时间,基于背景速度、背景密度和初始反射系数模型,应用弹性波方程进行线性化正演,得到正传波场τxx(x,z,t)、τxz(x,z,t)、τzz(x,z,t)、vx(x,z,t)、vz(x,z,t)和合成地震记录Cal(x,z,t);
(2)获得实际地震记录Obs(x,z,t),将其与合成地震记录的差值作为逆时扰动带入伴随方程得到反传波场Φxx(x,z,t)、Φxz(x,z,t)、Φzz(x,z,t)、φx(x,z,t)、φz(x,z,t);
(3)利用Poynting矢量行波分离公式(1),分别带入正传波场τxx(x,z,t)、τxz(x,z,t)、τzz(x,z,t)、vx(x,z,t)、vz(x,z,t)和反传波场Φxx(x,z,t)、Φxz(x,z,t)、Φzz(x,z,t)、φx(x,z,t)、φz(x,z,t),计算其Poynting矢量的水平分量和垂直分量和对于正传波场,所述的公式(1)为:反传波场的公式(1)能够类推得到;
(4)利用Poynting矢量方向指示公式(2)将正传波场vx(x,z,t)分解为正传波场上行波下行波左行波和右行波将正传波场vz(x,z,t)分解为正传波场上行波下行波左行波 和右行波将反传波场φx(x,z,t)分解为反传波场上行波下行波左行波和右行波将反传波场φz(x,z,t)分解为反传波场上行波下行波左行波和右行波 对于正传波场vx(x,z,t),所述的公式(2)为:反传波场的公式(2)能够类推得到;
(5)利用梯度公式(3)将行波分离后的正传波场和反传波场进行互相关成像,其中λ0为拉梅常数,公式(3)为
弹性波方程最小二乘逆时偏移每次迭代中的梯度剖面求取均需要重复以上5个步骤。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210525461.5A CN114624766B (zh) | 2022-05-16 | 2022-05-16 | 基于行波分离的弹性波最小二乘逆时偏移梯度求取方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210525461.5A CN114624766B (zh) | 2022-05-16 | 2022-05-16 | 基于行波分离的弹性波最小二乘逆时偏移梯度求取方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114624766A CN114624766A (zh) | 2022-06-14 |
CN114624766B true CN114624766B (zh) | 2022-08-02 |
Family
ID=81907138
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202210525461.5A Active CN114624766B (zh) | 2022-05-16 | 2022-05-16 | 基于行波分离的弹性波最小二乘逆时偏移梯度求取方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114624766B (zh) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117741746A (zh) * | 2023-12-07 | 2024-03-22 | 中国地震局地质研究所 | 断层介质成像方法、装置和电子设备 |
Family Cites Families (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105137486B (zh) * | 2015-09-01 | 2017-10-20 | 中国科学院地质与地球物理研究所 | 各向异性介质中弹性波逆时偏移成像方法及其装置 |
CN105467444B (zh) * | 2015-12-10 | 2017-11-21 | 中国石油天然气集团公司 | 一种弹性波全波形反演方法及装置 |
CN106526674B (zh) * | 2016-11-14 | 2020-03-31 | 中国石油化工股份有限公司 | 一种三维全波形反演能量加权梯度预处理方法 |
US10295685B2 (en) * | 2017-04-06 | 2019-05-21 | Saudi Arabian Oil Company | Generating common image gather using wave-field separation |
CN107272058B (zh) * | 2017-07-05 | 2020-02-11 | 中国科学院地质与地球物理研究所 | 成像方法、成像装置以及计算机存储介质 |
CN107153216B (zh) * | 2017-07-05 | 2019-05-07 | 中国科学院地质与地球物理研究所 | 确定地震波场的坡印廷矢量的方法、装置以及计算机存储介质 |
CN110531410B (zh) * | 2019-09-23 | 2020-04-14 | 中国海洋大学 | 一种基于直达波场的最小二乘逆时偏移梯度预条件方法 |
CN113031062B (zh) * | 2021-04-09 | 2022-01-28 | 中国海洋大学 | 一种基于波场分离的相关加权逆时偏移成像方法 |
-
2022
- 2022-05-16 CN CN202210525461.5A patent/CN114624766B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN114624766A (zh) | 2022-06-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108873066B (zh) | 弹性介质波动方程反射波旅行时反演方法 | |
CN104122585B (zh) | 基于弹性波场矢量分解与低秩分解的地震正演模拟方法 | |
CN110058307B (zh) | 一种基于快速拟牛顿法的全波形反演方法 | |
CN110187382B (zh) | 一种回折波和反射波波动方程旅行时反演方法 | |
CN105388520B (zh) | 一种地震资料叠前逆时偏移成像方法 | |
CN110007340B (zh) | 基于角度域直接包络反演的盐丘速度密度估计方法 | |
CN111948708B (zh) | 一种浸入边界起伏地表地震波场正演模拟方法 | |
CN111239819B (zh) | 一种基于地震道属性分析的带极性直接包络反演方法 | |
CN114624766B (zh) | 基于行波分离的弹性波最小二乘逆时偏移梯度求取方法 | |
CN110888159B (zh) | 基于角度分解与波场分离的弹性波全波形反演方法 | |
Li et al. | An immersed boundary method with iterative symmetric interpolation for irregular surface topography in seismic wavefield modelling | |
CN111665556B (zh) | 地层声波传播速度模型构建方法 | |
US11199641B2 (en) | Seismic modeling | |
CN114460640A (zh) | 有限差分模拟弹性波全波形反演方法和装置 | |
CN114924313B (zh) | 基于行波分离的声波最小二乘逆时偏移梯度求取方法 | |
CN107807392A (zh) | 一种自适应抗频散的分块时空双变逆时偏移方法 | |
CN102313902B (zh) | 一种基于切比雪夫展开的广义屏叠前深度偏移方法 | |
CN113126151B (zh) | 一种基于纯波延拓方程的弹性反射波旅行时反演方法 | |
CN111665550A (zh) | 地下介质密度信息反演方法 | |
US10451756B2 (en) | Apparatus and method for calculating efficient 3D traveltime by using coarse-grid mesh for shallow depth source | |
Liu et al. | Using pseudo-spectral method on curved grids for SH-wave modeling of irregular free-surface | |
CN114325829B (zh) | 一种基于双差思想的全波形反演方法 | |
CN115903042A (zh) | 一种基于构造整形正则化的波形反演方法及其设备 | |
CN108563802B (zh) | 一种提高地震转换波数值模拟精度的方法 | |
CN113960663A (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 |