CN107894618B - 一种基于模型平滑算法的全波形反演梯度预处理方法 - Google Patents
一种基于模型平滑算法的全波形反演梯度预处理方法 Download PDFInfo
- Publication number
- CN107894618B CN107894618B CN201711101512.7A CN201711101512A CN107894618B CN 107894618 B CN107894618 B CN 107894618B CN 201711101512 A CN201711101512 A CN 201711101512A CN 107894618 B CN107894618 B CN 107894618B
- Authority
- CN
- China
- Prior art keywords
- gradient
- waveform inversion
- full waveform
- model
- rate pattern
- 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
- 238000000034 method Methods 0.000 title claims abstract description 36
- 238000009499 grossing Methods 0.000 title claims abstract description 17
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 15
- 238000005516 engineering process Methods 0.000 claims abstract description 10
- 238000004088 simulation Methods 0.000 claims abstract description 10
- 238000004364 calculation method Methods 0.000 claims description 2
- 230000000694 effects Effects 0.000 abstract description 6
- 238000002474 experimental method Methods 0.000 abstract description 2
- 238000012937 correction Methods 0.000 description 3
- 238000007781 pre-processing Methods 0.000 description 3
- 238000005553 drilling Methods 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 230000005540 biological transmission Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
-
- 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
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
一种基于模型平滑算法的全波形反演梯度预处理方法，属于地震波形反演领域，具体包括以下步骤：每次迭代反演时，首先将迭代得到的速度模型进行平滑，然后基于此平滑速度模型，应用声波方程有限差分模拟技术，模拟得到炮点震源子波的正时波场能量和以实际地震记录为扰动的逆时波场能量，最后应用此正时波场能量和逆时波场能量对梯度进行预处理。模型实验结果表明本方法可有效降低波场能量计算中反射波及多次波的影响，提高梯度预处理的精度，从而显著提升全波形反演的效果。
Description
技术领域
本发明属于地震波形反演领域，具体涉及一种基于模型平滑算法的全波形反演梯度预处理方法。
背景技术
上世纪八十年代，Tarantola等人首先基于最小平方理论提出了时间域全波形反演方法，该方法充分利用实测地震记录的走时、振幅以及相位等信息来重建地下介质的速度结构，理论上对于理想观测系统其精度可达到波长数量级，因此全波形反演被认为是能够进一步提升地震勘探的油气勘探能力的重要方法，有望在未来的油气勘探中发挥重要作用。
全波形反演是以理论地震记录与实际地震记录误差的最小平方泛函作为反演的目标函数，应用双程波动方程有限差分数值模拟方法分别模拟炮点震源子波的正向传播波场与接收点理论地震记录与实际地震记录残差的逆时传播波场，然后应用正向及逆时波场的时间导数值生成梯度方向并计算共轭梯度方向，最后通过梯度寻优实现速度模型参数迭代更新。具体实现流程如图1所示。
常规地震全波形反演计算时，通常情况下模型浅层地震波能量要强于深层，并且浅层的地震射线密度也要高于深层，这就使得在迭代反演的过程中浅层的速度更新一直占据主导地位，从而导致深层速度反演收敛慢且精度不高。当前主要有两种策略解决深层反演精度低下的问题，一种是基于“层剥法”的思想，由浅层到深层逐层反演，其存在着难以做到精确层位划分的问题；另一种策略则是在反演的过程中对梯度进行预处理，以校正几何扩散对梯度计算的影响。因此梯度预处理方法在实际全波形反演中得到广泛应用。
Pratt与1998年提出了应用Hessian矩阵或近似Hessian矩阵对梯度进行预处理的方法，随后诸多学者在Hessian矩阵类梯度预处理方向上做了大量的研究工作。Hessian矩阵或近似Hessian的梯度预处理方法通常都需要存储Hessian矩阵并需进行求逆运算，而Hessian矩阵或近似Hessian维数通常都非常大(为模型网格点个数的平方)，因此若应用该类方法基于较大模型进行反演，将产生巨大的内存消耗且由此带来的大型矩阵的求逆也极易产生不稳定现象，因此该类方法仅适合于小模型数据的反演。
Zhang等于2012年SEG年会上提出了一种应用地震波能量进行梯度预处理的方法，其应用地下每个网格点正时、逆时波场的地震波能量对梯度进行加权，其本质是以波场能量来作为每个点的几何扩散校正因子，以使得浅中深部的梯度值达到均衡。Zhang的方法是基于当前迭代得到的模型进行波动方程模拟，获得各点的能量值，并以其作为几何扩散及照明度不均衡的校正因子。理论上，应用透射波场的能量作为几何扩散及照明度不均衡的校正因子更为精确，但由于当前迭代得到的模型往往已具备成层性，因此基于当前迭代得到的模型模拟得到的波场中除透射波能量外还包含了大量的反射波/多次波，反射波/多次波的能量会干扰梯度预处理的精度。
本发明基于模型平滑算法对全波形反演中的梯度进行预处理(具体流程见图2所示)，可有效降低波场中的反射波/多次波的影响，提高与处理的精度，从而显著提升反演的效果。
发明内容
本发明要解决的技术问题在于提供一种基于模型平滑算法的全波形反演梯度预处理方法。在全波形反演的每次迭代时，首先将迭代得到的速度模型进行平滑，然后基于此平滑速度模型，应用声波方程有限差分模拟技术，模拟得到炮点震源子波的正时波场能量和以实际地震记录为扰动的逆时波场能量，最后应用此正时波场能量和逆时波场能量对梯度进行预处理。该方法可有效降低波场能量计算中反射波及多次波的影响，提高梯度预处理的精度，从而显著提升全波形反演的效果
本发明采取以下技术方案：
一种基于模型平滑算法的全波形反演梯度预处理方法，其特征在于它具体包括以下步骤：
(1)基于模型，模拟得到炮地震记录作为反演的观测地震记录；在全波形反演的第i(i>1)次迭代计算时，输入上次的梯度值gi-1(x)，速度模型Mn、地震子波w及观测地震记录Pobs(xr,xs)，其中x表示空间位置坐标，r表示检波点，s表示炮点，xr和xs分别表示检波点和炮点位置；
(2)应用五点三次平滑方法，将速度模型Mn平滑N次(N>＝50)，得到新的平滑速度模型
(3)基于平滑速度模型应用声波方程有限差分模拟技术，以地震子波w为正时扰动，模拟得到正时波场P(x,t,xs)，其中t表示时间；
(4)基于平滑速度模型应用声波方程有限差分模拟技术，以观测地震记录Pobs(xr,xs)为逆时扰动，逆时延拓得到逆时波场U(x,t,xs)；
(5)利用公式(1)计算基于模型的正时波场P(x,t,xs)的能量Ws(x)；所述的公式(1)为
(6)利用公式(2)计算基于模型的逆时波场U(x,t,xs)的能量Wr(x)；所述的公式(2)为
(7)应用公式(3)对上次的梯度gi-1(x)进行预处理，得到新的梯度γi(x)；所述的公式(3)为
全波形反演中每次梯度预处理时均需重复以上7个步骤。
本发明与现有技术相比的有益效果：
本发明提出的一种基于模型平滑算法的全波形反演梯度预处理方法，每次迭代反演时，首先将迭代得到的速度模型进行平滑，然后基于此平滑速度模型，应用声波方程有限差分模拟技术，模拟得到炮点震源子波的正时波场能量和以实际地震记录为扰动的逆时波场能量，最后应用此正时波场能量和逆时波场能量对梯度进行预处理。模型实验结果表明本方法可有效降低波场能量计算中反射波及多次波的影响，提高梯度预处理的精度，从而显著提升全波形反演的效果，有望在实际数据的全波形反演中发挥重要作用。
附图说明
图1基于波场能量的全波形反演流程图；
图2基于模型平滑算法的梯度预处理的全波形反演流程图；
图3OverThrust速度模型；
图4初始速度模型；
图5基于波场能量梯度预处理的反演结果；
图6基于模型平滑算法梯度预处理的反演结果；
图7基于波场能量梯度预处理的反演结果(局部)；
图8基于模型平滑算法梯度预处理的反演结果(局部)。
具体实施方式
这里采用OverThrust模型(OverThrust速度模型见图3)进行具体实施方式阐述。OverThrust模型是一个测试地震数据处理、反演方法准确性的著名模型，模型横向长10000m、纵向深2325m。
基于OverThrust模型，采用全排列的观测方式，模拟得到共521炮地震记录作为反演的观测地震记录。该炮集记录每炮为521道接收，炮间隔和道间隔均为25m，炮点深度和接收点深度也均为25m。
全波形反演的初始模型如图4所示，下面详细阐述本发明的具体实施过程：
(1)输入上次迭代的梯度值gi-1(x)(i为当前迭代次数)，速度模型Mn、地震子波w及观测地震记录Pobs(xr,xs)，其中x表示空间位置坐标，r表示检波点，s表示炮点，xr和xs分别表示检波点位置和炮点位置；
(2)应用五点三次平滑方法，将速度模型Mn平滑50次，得到新的平滑速度模型
(3)基于平滑应用声波方程有限差分模拟技术，以地震子波w为正时扰动，模速度模型拟得到正时波场P(x,t,xs)，其中t表示时间；
(4)基于平滑速度模型应用声波方程有限差分模拟技术，以观测地震记录Pobs(xr,xs)为逆时扰动，逆时延拓得到逆时波场U(x,t,xs)；
(5)利用公式(1)计算基于模型的正时波场P(x,t,xs)的能量Ws(x)；
(6)利用公式(2)计算基于模型的逆时波场U(x,t,xs)的能量Wr(x)；
(7)应用公式(3)对上次的梯度gi-1(x)进行预处理，得到新的梯度γi(x)。
全波形反演过程中每次迭代时，均需重复以上7个步骤进行梯度预处理。需要特别说明的是在首次迭代(即i＝1)时，不需要进行梯度预处理(即不需要以上步骤处理)。
为了说明本发明所述方法的有效性，这里与基于常规的梯度预处理方法进行对比。如图5为基于常规梯度预处理的全波形反演方法得到的结果，图6为基于模型平滑算法的梯度预处理全波形反演得到的结果。为了更清晰的对比两个结果，将这两个速度剖面选取同样的区域进行局部放大显示，分别如图7、8所示。对比图7和图8可知，基于模型平滑算法的梯度预处理算法，其全波形反演结果的反演精度明显高于基于常规梯度预处理算法的全波形反演。
Claims (1)
1.一种基于模型平滑算法的全波形反演梯度预处理方法，其特征在于它具体包括以下步骤：
(1)基于模型，模拟得到炮地震记录作为反演的观测地震记录；在全波形反演的第i次迭代计算时，其中i>1，输入上次的梯度值gi-1(x)，速度模型Mn、地震子波w及观测地震记录Pobs(xr,xs)，其中x表示空间位置坐标，r表示检波点，s表示炮点，xr和xs分别表示检波点和炮点位置；
(2)应用五点三次平滑方法，将速度模型Mn平滑N次，其中N≥50，得到新的平滑速度模型
(3)基于平滑速度模型应用声波方程有限差分模拟技术，以地震子波w为正时扰动，模拟得到正时波场P(x,t,xs)，其中t表示时间；
(4)基于平滑速度模型应用声波方程有限差分模拟技术，以观测地震记录Pobs(xr,xs)为逆时扰动，逆时延拓得到逆时波场U(x,t,xs)；
(5)利用公式(1)计算基于模型的正时波场P(x,t,xs)的能量Ws(x)；所述的公式(1)为
(6)利用公式(2)计算基于模型的逆时波场U(x,t,xs)的能量Wr(x)；所述的公式(2)为
(7)应用公式(3)对上次的梯度gi-1(x)进行预处理，得到新的梯度γi(x)；所述的公式(3)为
全波形反演中每次梯度预处理时均需重复以上7个步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711101512.7A CN107894618B (zh) | 2017-11-10 | 2017-11-10 | 一种基于模型平滑算法的全波形反演梯度预处理方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711101512.7A CN107894618B (zh) | 2017-11-10 | 2017-11-10 | 一种基于模型平滑算法的全波形反演梯度预处理方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107894618A CN107894618A (zh) | 2018-04-10 |
| CN107894618B true CN107894618B (zh) | 2018-08-21 |
Family
ID=61804890
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201711101512.7A Active CN107894618B (zh) | 2017-11-10 | 2017-11-10 | 一种基于模型平滑算法的全波形反演梯度预处理方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107894618B (zh) |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109975869B (zh) * | 2019-03-27 | 2020-10-27 | 中国石油大学(北京) | 一种沿地层走向光滑约束的反射波波形反演方法 |
| CN110531410B (zh) * | 2019-09-23 | 2020-04-14 | 中国海洋大学 | 一种基于直达波场的最小二乘逆时偏移梯度预条件方法 |
| CN111208568B (zh) * | 2020-01-16 | 2021-04-09 | 中国科学院地质与地球物理研究所 | 一种时间域多尺度全波形反演方法及系统 |
| US11747499B2 (en) * | 2021-01-14 | 2023-09-05 | Cgg Services Sas | Full-waveform inversion using partial time shifts and enhanced kinematic transforms |
| CN112925023B (zh) * | 2021-02-01 | 2022-03-22 | 中国石油大学(北京) | 全波场反演地震数据多次波压制方法 |
| CN113296146B (zh) * | 2021-05-19 | 2022-02-22 | 中国海洋大学 | 一种基于梯度道集相关加权的全波形反演梯度预处理方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103091711A (zh) * | 2013-01-24 | 2013-05-08 | 中国石油天然气集团公司 | 全波形反演方法及装置 |
| CN104391323A (zh) * | 2014-11-21 | 2015-03-04 | 中国石油大学(华东) | 一种利用反射波信息反演速度场中低波数成分的方法 |
| CN106526674A (zh) * | 2016-11-14 | 2017-03-22 | 中国石油化工股份有限公司 | 一种三维全波形反演能量加权梯度预处理方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9176930B2 (en) * | 2011-11-29 | 2015-11-03 | Exxonmobil Upstream Research Company | Methods for approximating hessian times vector operation in full wavefield inversion |
-
2017
- 2017-11-10 CN CN201711101512.7A patent/CN107894618B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103091711A (zh) * | 2013-01-24 | 2013-05-08 | 中国石油天然气集团公司 | 全波形反演方法及装置 |
| CN104391323A (zh) * | 2014-11-21 | 2015-03-04 | 中国石油大学(华东) | 一种利用反射波信息反演速度场中低波数成分的方法 |
| CN106526674A (zh) * | 2016-11-14 | 2017-03-22 | 中国石油化工股份有限公司 | 一种三维全波形反演能量加权梯度预处理方法 |
Non-Patent Citations (4)
| Title |
|---|
| A guided tour of multiparameter full-waveform inversion with multicomponent data: From theory to practice;S. OPERTO 等;《THE LEADING EDGE》;20130930;第32卷(第9期);第1040-1054页 * |
| An overview of full-waveform inversion in exploration geophysics;J. Virieux 等;《GEOPHYSICS》;20091130;第74卷(第6期);第WCC127-WCC152页 * |
| 一种新的用于全波形反演的能量加权梯度方法;张凯 等;《石油物探》;20150131;第54卷(第1期);第77-82页 * |
| 基于共轭梯度法的全波形反演;林朋 等;《煤田地质与勘探》;20170228;第45卷(第1期);第131-136、142页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107894618A (zh) | 2018-04-10 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107894618B (zh) | 一种基于模型平滑算法的全波形反演梯度预处理方法 | |
| CN106526674B (zh) | 一种三维全波形反演能量加权梯度预处理方法 | |
| CN103713315B (zh) | 一种地震各向异性参数全波形反演方法及装置 | |
| CN105319581B (zh) | 一种高效的时间域全波形反演方法 | |
| WO2023087451A1 (zh) | 基于观测数据自编码的多尺度无监督地震波速反演方法 | |
| CN105093278B (zh) | 基于激发主能量优化算法的全波形反演梯度算子提取方法 | |
| CN110058302A (zh) | 一种基于预条件共轭梯度加速算法的全波形反演方法 | |
| CN106932819B (zh) | 基于各向异性马尔科夫随机域的叠前地震参数反演方法 | |
| CN104570082B (zh) | 一种基于格林函数表征的全波形反演梯度算子的提取方法 | |
| CN108345031A (zh) | 一种弹性介质主动源和被动源混采地震数据全波形反演方法 | |
| CN107765308B (zh) | 基于褶积思想与精确震源的重构低频数据频域全波形反演方法 | |
| CN107505654A (zh) | 基于地震记录积分的全波形反演方法 | |
| CN113740901B (zh) | 基于复杂起伏地表的陆上地震数据全波形反演方法及装置 | |
| CN104977607B (zh) | 利用变步长网格声波波场模拟的时间域全波形反演方法 | |
| CN105549079A (zh) | 一种地球物理参数的全波形反演模型的建立方法和装置 | |
| CN108645994A (zh) | 一种基于多点地质统计学的地质随机反演方法及装置 | |
| CN108508482A (zh) | 一种地下裂缝地震散射响应特征模拟方法 | |
| CN105911584B (zh) | 一种隐式交错网格有限差分弹性波数值模拟方法及装置 | |
| CN104965222B (zh) | 三维纵波阻抗全波形反演方法及装置 | |
| CN110058307A (zh) | 一种基于快速拟牛顿法的全波形反演方法 | |
| CN109188519A (zh) | 一种极坐标下的弹性波纵横波速度反演系统及方法 | |
| CN113552625A (zh) | 一种用于常规陆域地震数据的多尺度全波形反演方法 | |
| CN110531410A (zh) | 一种基于直达波场的最小二乘逆时偏移梯度预条件方法 | |
| CN104237937A (zh) | 叠前地震反演方法及其系统 | |
| CN111290019B (zh) | 一种应用于最小二乘逆时偏移的l-bfgs初始矩阵求取方法 |
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 |