CN109782350B - 模式识别自适应全波形反演方法 - Google Patents

模式识别自适应全波形反演方法 Download PDF

Info

Publication number
CN109782350B
CN109782350B CN201910084956.7A CN201910084956A CN109782350B CN 109782350 B CN109782350 B CN 109782350B CN 201910084956 A CN201910084956 A CN 201910084956A CN 109782350 B CN109782350 B CN 109782350B
Authority
CN
China
Prior art keywords
obtaining
matrix
wave field
full waveform
waveform inversion
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
Application number
CN201910084956.7A
Other languages
English (en)
Other versions
CN109782350A (zh
Inventor
国九英
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Guoyou Weitai Beijing Technology Co ltd
Original Assignee
Guoyou Weitai Beijing Technology Co ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Guoyou Weitai Beijing Technology Co ltd filed Critical Guoyou Weitai Beijing Technology Co ltd
Priority to CN201910084956.7A priority Critical patent/CN109782350B/zh
Publication of CN109782350A publication Critical patent/CN109782350A/zh
Application granted granted Critical
Publication of CN109782350B publication Critical patent/CN109782350B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Abstract

本发明公开了一种模式识别自适应全波形反演方法，其采用如下步骤：步骤1、对每一炮集记录完成如下步骤：1‑1计算每一时刻的正传播波场；1‑2计算伴随源
Figure DEST_PATH_IMAGE002
，并将其反传播得到反传播数据；1‑3计算反传播波场与正传播波场的互相关得到单炮的梯度；步骤2、将所有炮求取得梯度叠加，从而得到模型空间的全局梯度；步骤3、通过最速下降法或局部反演算法得到速度模型的修改量，进而得到优化后的速度模型。本发明没有只能利用地震数据的低频成分的限制，可以利用地震记录的中频、高频成分，比如主频范围内的频率成分等，且没有常规全波形反演技术发生周期跳等导致错误求解的缺陷。

Description

模式识别自适应全波形反演方法
技术领域
本发明涉及一种模式识别自适应全波形反演方法，其属于油气勘探应用领域。
背景技术
石油天然气资源不仅对国家的经济意义而且国家的战略意义都十分重大。我国是一个能源需求大国，每年都要花去大量的资金购买国外的石油，这样除给国家造成巨大经济损失外，还给国家的能源战略提出巨大挑战。如何寻找埋藏于地下几千米到上万米深的石油天然气资源是目前世界各大石油公司面临的巨大难题和努力的目标。它要求首先要得到地下三维的图像，这个过程也就是通常所说的地球物理勘探，然后根据一些原则判断地下有无石油天然气资源的存在、预测石油天然气的储量、提供钻井井位等。众所周知，地下是看不见摸不着的，只能通过在地球表面人工产生地震波，这些地震波传入到地球内部，当这些地震波传入到地下不同深度后往回反射弹性波到地球表面，这些反射波以数字信号形式被高灵敏度专用数字记录仪记录下来。
通过人工激发并接收回来的地震数字信号利用体现地震波传播规律的波动方程来反演地下地质结构，从而得到地下的高精度成像。为达到此目的必须有高精度的地下速度模型。如何得到地下高精度的速度模型一直是人们努力的目标，但是由于其技术的高难度，一直在该领域没有得到突破性进展。
随着地球物理勘探的深入，人们在努力研究如何得到高精度的速度模型，为此提出了全波形反演这项技术，该技术在一定程度上取得了一定效果，但是该技术只能利用地震记录的低频分量来反演地下速度模型，否则就会出现周期跳这个令人十分头疼的问题，但是遗憾的是地震记录的低频成分往往受干扰波的影响比较大，信噪比比较低，比如面波或海上波浪产生的浪涌等干扰波都以强能量低频形式出现，它们影响地震信号的低频成分，而常规全波形反演技术恰恰只能利用地震波的低频成分，显然这是我们不希望的。另外地震波的低频成分不能得到速度模型的高频成分，也就是说得不到速度模型的细节变化，这样也制约了速度模型求取的精度。
发明内容
本发明所要解决的技术问题是提供了一种不限于利用低频数据并且无“周期跳”的高精度模式识别自适应全波形反演方法。
本发明采用了如下技术方案：
一种模式识别自适应全波形反演方法，其实现包括如下步骤：
步骤1、对每一炮集记录完成如下步骤：
(1)计算每一时刻的正传播波场；
(2)计算伴随源δs，并将其反传播得到反传播数据；
(3)计算反传播波场与正传播波场的互相关得到单炮的梯度；
步骤2、将所有炮求取得梯度叠加，从而得到模型空间的全局梯度；
步骤3、通过最速下降法或局部反演算法得到速度模型的修改量，进而得到优化后的速度模型。
进一步的，所述伴随源δs为如下公式(6)：
Figure BDA0001961454000000021
其中，
Figure BDA0001961454000000022
式(4)中，矩阵T是对w加权的对角矩阵；w为滤波算子；
式(6)中，wT表示与滤波算子w的互相关；
I为单位矩阵；
上角标T代表共轭转置；
P是模拟的地震数据的Toeplitz矩阵；
Ep是信号的投影预测误差滤波算子矩阵；
Bp是信号的投影预测滤波算子矩阵；
矩阵T是对w加权的对角矩阵。
本发明的有益效果如下：
本发明没有只能利用地震数据的低频成分的限制，可以利用地震记录的中频、高频成分，比如主频范围内的频率成分等，且没有常规全波形反演技术发生周期跳等导致错误求解的缺陷。本发明可以提供高精度的地下速度模型，进而可以得到地下高精度的地震成像，从而为高精度地震探测地下石油天然气及其它矿物和地热资源等提供了可靠的基础，也为石油天然气等公司减少风险和降低花费提供强有力的技术保证，具有广阔的用户市场和经济价值。
附图说明
图1为真实的速度模型。
图2为用常规全波形反演方法得到的速度模型。
图3为用模式识别自适应全波形反演方法得到的速度模型。
图4为某实际区块用模式识别自适应全波形反演方法得到的速度模型的速度优化量。
图5为用常规全波形反演方法得到的速度模型做的叠前深度偏移成像。
图6为用模式识别自适应全波形反演方法得到的速度模型做的叠前深度偏移成像。
具体实施方式
为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图1-6和具体实施方式对本发明进行进一步的详细描述。需要说明的是，在不冲突的情况下，本申请的实施例及实施例中的特征可以相互结合。
在下面的描述中，阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述的其他方式来实施，因此，本发明的保护范围并不受下面公开的具体实施例的限制。
本发明的方案可应用在石油天然气的勘探、煤田勘测、地热勘测、水文勘测，以及地震预测防灾等。现有技术不能满足高精度地球物理勘探的精度要求，所以无法高精度地发现储存在地下构造中的油气资源，尤其是当地下构造复杂时更是如此。
本实施例涉及一种模式识别自适应全波形反演方法，其原理如下：
设计目标函数g的公式(1)如下：
Figure BDA0001961454000000031
Figure BDA0001961454000000032
Figure BDA0001961454000000033
其中，P是模拟的地震数据的Toeplitz矩阵，即计算出的正传播波场；
D是输入的实际地震数据矩阵；
w是滤波算子，被应用到模拟的地震记录上从而使模拟的地震记录与实际观测到的地震记录有好的匹配；
Fp是信号的投影预测误差滤波算子矩阵；
Bp是信号的投影预测滤波算子矩阵；
ε是预白化系数，即稳定常数；
上角标T代表共轭转置；
I是单位矩阵。
公式(1)是一个线性最小平方问题，它有一个全局最小值，没有第二个局部最小值，所以没有常规全波形反演所具有的“周期跳”问题。
采用如下公式(4)替代的目标函数g，公式(4)如下所示：
Figure BDA0001961454000000041
其中，矩阵T是对w加权的对角矩阵；
通过速度模型m对公式(4)中的目标函数f最小化，得到反传播波场与正传播波场的互相关，进而得到单炮对应的梯度，即如下公式(5)：
Figure BDA0001961454000000042
其中，矩阵A是波动方程正演模拟算子的矩阵表达；
R是选择检波点位置的波场的限定算子；
u表示震源s在速度模型m中在所有点和时间产生的波场；
δS是伴随源；
其中，伴随源δS为如下公式(6)：
Figure BDA0001961454000000043
其中，wT表示与滤波算子w的互相关；
I为单位矩阵；
上角标T代表共轭转置；
f代表公式(4)。
通过上述所述方法，常规全波形反演技术固有的周期跳问题可以得到有效解决，从而全波形反演可以利用地震数据有效频带内的信息，不再有只能利用地震记录低频信息的限制。
图1为真实的速度模型，用基于此速度模型产生的地震记录做全波形速度反演。
图2为用常规全波形反演技术得到的速度模型，可以清楚的看到此模型与图1中的真实速度模型相差甚远，也就是说反演结果的精度相对不高。
图3为用模式识别自适应全波形反演得到的速度模型，可以清楚的看出它与图1中真实的速度模型很接近，也说明它的反演速度模型精度明显好于常规全波形反演的精度。
图4为某实际区块用模式识别自适应全波形反演得到的速度模型与初始速度模型相减得到的速度变化量；从图4中可以看出新速度细节明显更加清楚，也就是说速度的高频成分更加丰富。
图5为用常规全波形反演得到的速度模型做的叠前深度偏移成像，可以看出圆圈所围部分出现了一个小凹陷，实际钻井证明这个小凹陷是不存在的。
图6为用模式识别自适应全波形反演得到的速度模型做的叠前深度偏移成像，从图6中可以清楚的看出原来在图5中的圆圈所围的小凹陷不存在了，且与实际钻井深度吻合。
综上所述，说明模式识别自适应全波形反演方法的速度模型的精度明显好于常规全波形反演得到的速度模型的精度。这为后面的基于波动方程的叠前深度偏移成像提供了坚实的速度基础，从而可以大大提高地球物理勘探的精度。
以上结合附图详细说明了本发明的技术方案，本发明中的步骤可根据实际需求进行顺序调整、合并和删减。
尽管参考附图详地公开了本发明，但应理解的是，这些描述仅仅是示例性的，并非用来限制本发明的应用。本发明的保护范围由附加权利要求限定，并可包括在不脱离本发明保护范围和精神的情况下针对发明所作的各种变型、改型及等效方案。

Claims (1)

1.模式识别自适应全波形反演方法，其特征在于，其包括如下步骤：
步骤1、对每一炮集记录完成如下步骤：
(1)计算每一时刻的正传播波场；
(2)计算伴随源δs，并将其反传播得到反传播数据；
(3)计算反传播波场与正传播波场的互相关得到单炮的梯度；
步骤2、将所有炮求取得梯度叠加，从而得到模型空间的全局梯度；
步骤3、通过最速下降法或局部反演算法得到速度模型的修改量，进而得到优化后的速度模型；
所述伴随源δs为如下公式(6)：
Figure FDA0003007534070000011
其中，
Figure FDA0003007534070000012
式(4)中，w为滤波算子；
式(6)中，
上角标T代表共轭转置；
wTw表示wT与滤波算子w的互相关；
I为单位矩阵；
P是模拟的地震数据的Toeplitz矩阵；
EP是信号的投影预测误差滤波算子矩阵；
BP是信号的投影预测滤波算子矩阵；
矩阵T′是对w加权的对角矩阵；
所述投影预测滤波算子矩阵BP为如下公式(3)：
Figure FDA0003007534070000013
式(3)中，
ε是预白化系数，即稳定常数。
CN201910084956.7A 2019-01-29 2019-01-29 模式识别自适应全波形反演方法 Active CN109782350B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910084956.7A CN109782350B (zh) 2019-01-29 2019-01-29 模式识别自适应全波形反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910084956.7A CN109782350B (zh) 2019-01-29 2019-01-29 模式识别自适应全波形反演方法

Publications (2)

Publication Number Publication Date
CN109782350A CN109782350A (zh) 2019-05-21
CN109782350B true CN109782350B (zh) 2021-06-01

Family

ID=66503148

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910084956.7A Active CN109782350B (zh) 2019-01-29 2019-01-29 模式识别自适应全波形反演方法

Country Status (1)

Country Link
CN (1) CN109782350B (zh)

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9448317B2 (en) * 2010-08-19 2016-09-20 Pgs Geophysical As Method for swell noise detection and attenuation in marine seismic surveys
KR101182839B1 (ko) * 2010-08-26 2012-09-14 서울대학교산학협력단 송신원 추정을 통한 시간 영역 역시간 구조보정 방법 및 장치
EP2622457A4 (en) * 2010-09-27 2018-02-21 Exxonmobil Upstream Research Company Simultaneous source encoding and source separation as a practical solution for full wavefield inversion
GB2509223B (en) * 2013-10-29 2015-03-18 Imp Innovations Ltd Method of, and apparatus for, full waveform inversion

Also Published As

Publication number Publication date
CN109782350A (zh) 2019-05-21

Similar Documents

Publication Publication Date Title
Landro et al. The Gullfaks 4D seismic study
Graves Modeling three-dimensional site response effects in the Marina District Basin, San Francisco, California
US6188964B1 (en) Method for using global optimization to the estimation of surface-consistent residual statics
Barclay et al. Seismic inversion: Reading between the lines
CN109669212B (zh) 地震数据处理方法、地层品质因子估算方法与装置
Frankel et al. Three-dimensional simulations of ground motions in the Seattle region for earthquakes in the Seattle fault zone
Kelamis et al. Velocity-independent redatuming: A new approach to the near-surface problem in land seismic data processing
CN104516018A (zh) 一种地球物理勘探中岩性约束下的孔隙度反演方法
CN111722284B (zh) 一种基于道集数据建立速度深度模型的方法
WO2012131356A2 (en) A process for characterising the evolution of an oil or gas reservoir over time
Talukdar et al. Sub-basalt imaging of hydrocarbon-bearing Mesozoic sediments using ray-trace inversion of first-arrival seismic data and elastic finite-difference full-wave modeling along Sinor–Valod profile of Deccan Syneclise, India
Ding et al. Reverse time migration (RTM) imaging of iron-oxide deposits in the Ludvika mining area, Sweden
Mittet et al. Fast finite-difference modeling of 3-D elastic wave propagation
CN109782350B (zh) 模式识别自适应全波形反演方法
Rossi et al. Joint 3D inversion of SWD and surface seismic data
Murty et al. Velocity structure of the West-Bengal sedimentary basin, India along the Palashi-Kandi profile using a travel-time inversion of wide-angle seismic data and gravity modeling–an update
Li et al. Novel strategies for complex foothills seismic imaging—Part 1: Mega-near-surface velocity estimation
Nanda Seismic modelling and inversion
Burianyk et al. Broadside wide-angle seismic studies and three-dimensional structure of the crust in the southeast Canadian Cordillera
Li et al. Application of discontinuous-Galerkin-based acoustic FWI on land data from the mountainous region of strong surface topography variations in China
Pullammanappallil et al. Use of advanced data processing techniques in the imaging of the Coso geothermal field
OSINOWO Data workflow for generating static geological model platform for petroleum prospect re-evaluation of Wytch Farm field, Wessex Basin, South Coast, United Kingdom.
Alhajni Comparison study between the elevation and datum statics in NC 210, western Libya
Young Reverse-Time Migration of a Methane Gas Hydrate Distributed Acoustic Sensing Three-Dimensional Vertical Seismic Profile Dataset
Wei et al. Unlocking A Global Ocean Mixing Dataset: toward Standardization of Seismic-derived Ocean Mixing Rates

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