CN109212602A - 一种改进地震数据分辨率的反射系数反演方法 - Google Patents
一种改进地震数据分辨率的反射系数反演方法 Download PDFInfo
- Publication number
- CN109212602A CN109212602A CN201811031301.5A CN201811031301A CN109212602A CN 109212602 A CN109212602 A CN 109212602A CN 201811031301 A CN201811031301 A CN 201811031301A CN 109212602 A CN109212602 A CN 109212602A
- Authority
- CN
- China
- Prior art keywords
- formula
- reflection coefficient
- low
- matrix
- seismic
- 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
- 238000000034 method Methods 0.000 title claims abstract description 34
- 239000011159 matrix material Substances 0.000 claims abstract description 25
- 238000001228 spectrum Methods 0.000 claims abstract description 4
- 239000002131 composite material Substances 0.000 description 3
- 230000011514 reflex Effects 0.000 description 3
- 230000015572 biosynthetic process Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 239000004615 ingredient Substances 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 238000002310 reflectometry Methods 0.000 description 1
- 239000011435 rock Substances 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/30—Analysis
- G01V1/301—Analysis for determining seismic cross-sections or geostructures
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
本发明公开了一种改进地震数据分辨率的反射系数反演方法。本发明的技术要点是,它是一种基于L1范数的反演方法,用于反演地震反射系数;先在此基础上引入低频模型来控制反射系数的幅值变化范围;再引入光滑矩阵F;比较在引入低通滤波器F前后地震反射系数的变化;在构建低通滤波器时,首先在频率域构建其频谱,之后变换到时间域,得到滤波器函数F(t),最后构建滤波器矩阵;求解引入低频模型来控制反射系数后的线性方程式,即可求出式中的反射系数数组R。本发明的方法通过使用低通滤波器,改进算法,将反射系数反演的频带宽度提至更高,因而使地震数据分辨率也得以改进。
Description
技术领域
本发明属于地震勘探中地震反演和油藏描述技术领域。具体涉及一种改进地震数据分辨率的反射系数反演方法。
背景技术
使用地震波反演来指导油藏描述主要基于两种原理:(1)基于波动方程,通过地震波场正演生成的地震合成记录和野外观测的地震记录的差异,来不断迭代更新油藏速度模型,直到差异最小时,我们认为此时的反演得到的模型接近真实模型,是我们的最终解;(2)基于褶积模型,即我们认为地震反射记录是地震反射系数和地震子波褶积得到的,在从地震和测井数据出发,估计出地震子波之后,可以使用子波和反射记录推算出地震反射系数,结合由地震、测井和地震解释结果得到的低频速度模型,即可推算出地层速度或波阻抗。使用上面两种方法得到的地层速度或波阻抗模型,我们即可对油藏进行定性或定量的描述,以确定岩石岩性和流体性质。
由于上述的第一种方法的输入数据是叠前数据,而且反演过程包括波动方程的正演过程,计算量巨大,尤其是在三维情况下,目前计算机的计算能力远远不够,所以第一种方法多在二维条件下应用,虽然此方法是目前研究的热点之一,但是实际应用的规模还是相当有限。相比较而言,第二种方法是在实际油藏描述中使用最为广泛的反演方法,这种方法基于褶积模型,本质上是一维反演技术,所以输入的地震数据必须是角道集数据或叠后数据,因此数据量相对较少,计算速度快,大型三维处理方式成熟。
第二种方法的原理简单,但是由于地震数据的频带宽度有限,基于常规反演方法得到的速度或波阻抗剖面光滑,不利于油藏描述。为了提高地震反演的频带宽度,许多学者提出了稀疏反演概念,以反演地层反射系数,提高分辨率。本发明使用稀疏反演的方法提高频率带宽,并且针对稀疏反演方法的弊端,提出低频滤波器,改进算法,将频带宽度提至更高。
发明内容
本发明的目的在于针对现有稀疏反演方法的弊端,提供一种改进地震数据分辨率的反射系数反演方法。
本发明的目的是通过如下的技术方案来实现:该改进地震数据分辨率的反射系数反演方法,是一种基于L1范数的反演方法,用于反演地震反射系数,具体实现是求取式(1)的最小值对应的解:
式(1)中,S0是观测地震道数据数组,W是地震子波矩阵,R是反射系数数组,λ1是L1范数规则化加权因子;式(1)中第一项是观测数据和合成数据的残差L2范数的平方,第二项是反射系数的L1范数的加权;
它还包括如下步骤:
(一)在式(1)中,引入低频模型来控制反射系数的幅值变化范围,得式(2):
式(2)中,λ2是低频控制加权因子,RLFM是低频模型的反射系数数组,L是下三角单位矩阵,完成对反射系数的积分,其形式如式(3):
(二)在式(2)中,引入光滑矩阵F,对LR项进行光滑,得式(4):
式(4)中,光滑矩阵F是一个低通滤波器,滤波范围设为0-6HZ;
(三)比较式(2)和式(4)在引入低通滤波器F前后地震反射系数的变化;由于式(2)和式(4)中含有地震子波矩阵,因此,需要在反演前使用地震数据和测井数据进行标定,估计出合理的地震子波,之后利用地震子波构建地震子波矩阵如式(5):
式(5)中,w1、w2、w3…是估计得到地震子波的第一个、第二个、第三个…的振幅值;
(四)在构建低通滤波器时,首先在频率域构建其频谱,之后变换到时间域,得到滤波器函数F(t),最后构建滤波器矩阵如式(6):
式(6)中,f1、f2、f3…是滤波器函数F(t)第一个、第二个、第三个…的振幅值;
(五)λ1、λ2通过测试得到,介于0.1到0.01之间;最后解式(4)的最终形式的线性方程如式(7):
[[GTG]-1GTG+λ1I]R=[GTG]-1GTS (7);
式(7)中:
I是单位矩阵;
使用求解线性方程的方法,即可求出式(7)中的反射系数数组R。
本发明的方法通过使用低通滤波器,改进算法,将反射系数反演的频带宽度提至更高,因而使地震数据分辨率也得以改进。
附图说明
图1和图2是本发明实施例未使用低通滤波器时应用式(2)的求解结果曲线图。图1中,实线曲线即为没有使用低频滤波器的反射系数曲线,虚线曲线为低频反射系数的积分曲线,点线曲线为反演到反射系数积分曲线;图2中,虚线曲线为地震道信号曲线,即实际记录曲线,实线曲线为合成记录曲线。
图3和图4是本发明实施例使用低通滤波器时应用式(7)的求解结果曲线图。图3中,实线曲线即为使用低频滤波器的反射系数曲线,虚线曲线为低频反射系数的积分曲线,点线曲线为反演到反射系数积分曲线;图4中,虚线曲线为地震道信号曲线,即实际记录曲线,实线曲线为合成记录曲线。
图5和图6分别是本发明实施例F低通滤波器在频率域和时间域的示意图。
图7是本发明实施例未使用低通滤波器的反射系数剖面图。
图8是本发明实施例使用低通滤波器的反射系数剖面图。
具体实施方式
下面结合附图和实施例对本发明作进一步详细的描述。
本发明方法是一种基于L1范数的反演方法,用于反演地震反射系数,具体实现是求取式(1)的最小值对应的解:
式(1)中,S0是观测地震道数据数组,W是地震子波矩阵,R是反射系数数组,λ1是L1范数规则化加权因子;式(1)中第一项是观测数据和合成数据的残差L2范数的平方,第二项是反射系数的L1范数的加权。
求解上式的最小值的解,即可得到尖脉冲反射系数数组,但是本质上不是真正的尖脉冲,而是具有一定宽度的尖脉冲,式(1)通常可以反演出反射系数对应的时间位置,但是反射系数的振幅和宽度没有得到有效控制,因此,本发明引入低频模型来控制反射系数的幅值变化范围,得式(2):
式(2)中,λ2是低频控制加权因子,RLFM是低频模型的反射系数数组(LFM表示LowFrequency Model),L是下三角单位矩阵,完成对反射系数的积分,其形式如式(3):
基于式(2),如果理想情况下R的解是一系列尖脉冲,LR的结果就是阶梯状的积分结果,而LRLFM是低频分量平滑的曲线,LR与LRLFM比较时,两者差异较大,这样会影响式(2)的第一项值,为了减小这两项的差异,在式(2)中,引入光滑矩阵F,对LR项进行光滑,取LR的低频分量,以便和LRLFM接近,得式(4):
式(4)中,光滑矩阵F是一个低通滤波器,滤波范围设为0-6HZ,滤掉高频的成分,这样,FLR与LRLFM具备相同的频带范围,有利于减小式(2)第一项的值,提高反射系数R的尖锐度,提高地震资料反演频带范围。式(4)是本发明的主要控制式,能提高反射系数反演的频带宽度。
接下来说明如何求解式(4)。比较式(2)和式(4)在引入低通滤波器F前后地震反射系数的变化;由于式(2)和式(4)中含有地震子波矩阵,因此,需要在反演前使用地震数据和测井数据进行标定,估计出合理的地震子波,之后利用地震子波构建地震子波矩阵如式(5):
式(5)中,w1、w2、w3…是估计得到地震子波的第一个、第二个、第三个…的振幅值。
在构建低通滤波器时,首先在频率域构建其频谱,如图5所示,之后变换到时间域,如图6所示,得到滤波器函数F(t),最后构建滤波器矩阵如式(6):
式(6)中,f1、f2、f3…是滤波器函数F(t)第一个、第二个、第三个…的振幅值。
对于λ1、λ2,主要通过测试得到,一般介于0.1到0.01之间;得最后解式(4)的最终形式的线性方程如式(7):
[[GTG]-1GTG+λ1I]R=[GTG]-1GTS (7);
式(7)中:
I是单位矩阵;
使用求解线性方程的方法,即可求出式(7)中的反射系数数组R。
下面是运用本发明方法的实例:
我们首先提取野外地震记录道来做反演,地震道信号如图2中的虚线曲线所示,我们使用式(2)的求解结果见图1和图2;图1中的实线曲线即为没有使用低通滤波器的反射系数曲线,虚线曲线为低频反射系数的积分曲线,点线曲线为反演得到反射系数积分曲线;图2中的实线曲线为合成记录。
图3、图4中的各个曲线和图1、图2中的一一对应,不同之处是图3、图4中我们使用了低通滤波器。比较图1、图2、图3、图4中的反射系数数组(实线曲线),我们不难发现,图3、图4中的反射系数更尖锐,说明图3、图4的反射系数数组的频带宽度更宽,更利于油藏描述。
图7和图8是在应用低通滤波器前后的反射系数剖面结果,这是一个实际数据的处理结果。我们比较图7和图8中矩形框里面的信号,不难发现,图8的反射系数更为集中,表示使用低通滤波器后图8的反射系数系列的频带宽度更宽,由此可见本发明方法的效果显著。
最后应说明的是:上述仅用以说明本发明而并非限制本发明所描述的技术方案;尽管本说明书对本发明已进行了详细的说明,但是,本领域的技术人员仍然可以对本发明进行修改或等同替换,一切不脱离本发明的精神和范围的技术方案及其改进,其均应涵盖在本发明的权利要求范围中。
Claims (1)
1.一种改进地震数据分辨率的反射系数反演方法,是一种基于L1范数的反演方法,用于反演地震反射系数,具体实现是求取式(1)的最小值对应的解:
式(1)中,S0是观测地震道数据数组,W是地震子波矩阵,R是反射系数数组,λ1是L1范数规则化加权因子;式(1)中第一项是观测数据和合成数据的残差L2范数的平方,第二项是反射系数的L1范数的加权;
其特征在于包括如下步骤:
(一)在式(1)中,引入低频模型来控制反射系数的幅值变化范围,得式(2):
式(2)中,λ2是低频控制加权因子,RLFM是低频模型的反射系数数组,L是下三角单位矩阵,完成对反射系数的积分,其形式如式(3):
(二)在式(2)中,引入光滑矩阵F,对LR项进行光滑,得式(4):
式(4)中,光滑矩阵F是一个低通滤波器,滤波范围设为0-6HZ;
(三)比较式(2)和式(4)在引入低通滤波器F前后地震反射系数的变化;由于式(2)和式(4)中含有地震子波矩阵,因此,需要在反演前使用地震数据和测井数据进行标定,估计出合理的地震子波,之后利用地震子波构建地震子波矩阵如式(5):
式(5)中,w1、w2、w3…是估计得到地震子波的第一个、第二个、第三个…的振幅值;
(四)在构建低通滤波器时,首先在频率域构建其频谱,之后变换到时间域,得到滤波器函数F(t),最后构建滤波器矩阵如式(6):
式(6)中,f1、f2、f3…是滤波器函数F(t)第一个、第二个、第三个…的振幅值;
(五)λ1、λ2通过测试得到,介于0.1到0.01之间;最后解式(4)的最终形式的线性方程如式(7):
[[GTG]-1GTG+λ1I]R=[GTG]-1GTS (7);
式(7)中:
I是单位矩阵;
使用求解线性方程的方法,即可求出式(7)中的反射系数数组R。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811031301.5A CN109212602B (zh) | 2018-09-05 | 2018-09-05 | 一种改进地震数据分辨率的反射系数反演方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811031301.5A CN109212602B (zh) | 2018-09-05 | 2018-09-05 | 一种改进地震数据分辨率的反射系数反演方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109212602A true CN109212602A (zh) | 2019-01-15 |
CN109212602B CN109212602B (zh) | 2019-11-08 |
Family
ID=64986348
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201811031301.5A Active CN109212602B (zh) | 2018-09-05 | 2018-09-05 | 一种改进地震数据分辨率的反射系数反演方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109212602B (zh) |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN1797037A (zh) * | 2004-12-29 | 2006-07-05 | 中国石油天然气集团公司 | 一种地震波波阻抗反演的方法 |
CN103487835A (zh) * | 2012-06-12 | 2014-01-01 | 中国石油化工股份有限公司 | 一种基于模型约束的多分辨率波阻抗反演方法 |
WO2016044538A1 (en) * | 2014-09-19 | 2016-03-24 | Conocophillips Company | Bandwidth extension beyond the vibrator sweep signal via a constrained simultaneous multiple vibrator inversion |
CN105467451A (zh) * | 2016-01-13 | 2016-04-06 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 基于全变差最小化约束的地震反射系数反演方法 |
CN106569262A (zh) * | 2015-10-12 | 2017-04-19 | 中国石油化工股份有限公司 | 低频地震数据缺失下的背景速度模型重构方法 |
CN107589448A (zh) * | 2017-07-13 | 2018-01-16 | 西安交通大学 | 一种多道地震记录反射系数序列同时反演方法 |
-
2018
- 2018-09-05 CN CN201811031301.5A patent/CN109212602B/zh active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN1797037A (zh) * | 2004-12-29 | 2006-07-05 | 中国石油天然气集团公司 | 一种地震波波阻抗反演的方法 |
CN103487835A (zh) * | 2012-06-12 | 2014-01-01 | 中国石油化工股份有限公司 | 一种基于模型约束的多分辨率波阻抗反演方法 |
WO2016044538A1 (en) * | 2014-09-19 | 2016-03-24 | Conocophillips Company | Bandwidth extension beyond the vibrator sweep signal via a constrained simultaneous multiple vibrator inversion |
CN106569262A (zh) * | 2015-10-12 | 2017-04-19 | 中国石油化工股份有限公司 | 低频地震数据缺失下的背景速度模型重构方法 |
CN105467451A (zh) * | 2016-01-13 | 2016-04-06 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 基于全变差最小化约束的地震反射系数反演方法 |
CN107589448A (zh) * | 2017-07-13 | 2018-01-16 | 西安交通大学 | 一种多道地震记录反射系数序列同时反演方法 |
Non-Patent Citations (2)
Title |
---|
张丰麒 等: "基于低频软约束的叠前AVA稀疏层反演", 《石油地球物理勘探》 * |
韩立国 等: "基于压缩感知和稀疏反演的地震数据低频补偿", 《吉林大学学报(地球科学版)》 * |
Also Published As
Publication number | Publication date |
---|---|
CN109212602B (zh) | 2019-11-08 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Chen et al. | Iterative deblending of simultaneous-source seismic data using seislet-domain shaping regularization | |
Shi et al. | Reverse time migration of 3D vertical seismic profile data | |
CN109669212B (zh) | 地震数据处理方法、地层品质因子估算方法与装置 | |
CN106033124B (zh) | 一种基于随机最优化的多震源粘声最小二乘逆时偏移方法 | |
CN108108331B (zh) | 一种基于拟空间域弹性波方程的有限差分计算方法 | |
CN105652321A (zh) | 一种粘声各向异性最小二乘逆时偏移成像方法 | |
Wang et al. | Seismic, waveform modeling and tomography | |
Sun et al. | Deep learning for low-frequency extrapolation of multicomponent data in elastic FWI | |
CN111596366A (zh) | 一种基于地震信号优化处理的波阻抗反演方法 | |
CN104965222B (zh) | 三维纵波阻抗全波形反演方法及装置 | |
Zhang et al. | 3D simultaneous seismic data reconstruction and noise suppression based on the curvelet transform | |
CN110703331A (zh) | 一种基于常q粘滞声波方程的衰减补偿逆时偏移实现方法 | |
CN110542928A (zh) | 基于vti各向异性传播矩阵的地震响应模拟方法 | |
CN107179550A (zh) | 一种数据驱动的地震信号零相位反褶积方法 | |
CN111077567B (zh) | 一种基于矩阵乘法的双程波叠前深度偏移的方法 | |
CN110687597B (zh) | 一种基于联合字典的波阻抗反演方法 | |
Bai et al. | Least-squares Gaussian beam transform for seismic noise attenuation | |
Chen et al. | 3-D seismic diffraction separation and imaging using the local rank-reduction method | |
CN106574980A (zh) | 用于地下地质体的岩石性质估计的系统和方法 | |
Yang et al. | A seismic interpolation and denoising method with curvelet transform matching filter | |
CN109212602B (zh) | 一种改进地震数据分辨率的反射系数反演方法 | |
Huang et al. | Pure qP-wave least-squares reverse time migration in vertically transverse isotropic media and its application to field data | |
Zhang et al. | Deep-learning for accelerating prestack correlative least-squares reverse time migration | |
CN114740528A (zh) | 一种超微分拉普拉斯块约束的叠前多波联合反演方法 | |
Wellington et al. | Preconditioning full-waveform inversion with efficient local correlation operators |
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 |