CN117310802A - 一种深度域储层地震反演方法 - Google Patents
一种深度域储层地震反演方法 Download PDFInfo
- Publication number
- CN117310802A CN117310802A CN202311183468.4A CN202311183468A CN117310802A CN 117310802 A CN117310802 A CN 117310802A CN 202311183468 A CN202311183468 A CN 202311183468A CN 117310802 A CN117310802 A CN 117310802A
- Authority
- CN
- China
- Prior art keywords
- time domain
- seismic data
- seismic
- depth
- domain 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 50
- 238000001914 filtration Methods 0.000 claims abstract description 11
- 230000006870 function Effects 0.000 claims description 42
- 238000005070 sampling Methods 0.000 claims description 37
- 239000011159 matrix material Substances 0.000 claims description 24
- 238000006243 chemical reaction Methods 0.000 claims description 14
- 238000001228 spectrum Methods 0.000 claims description 9
- 230000004927 fusion Effects 0.000 claims description 4
- 238000012937 correction Methods 0.000 claims description 3
- 125000004122 cyclic group Chemical group 0.000 claims description 3
- 230000002194 synthesizing effect Effects 0.000 claims description 3
- 238000010586 diagram Methods 0.000 description 6
- 238000012545 processing Methods 0.000 description 6
- 238000004590 computer program Methods 0.000 description 5
- 230000008569 process Effects 0.000 description 5
- 238000011160 research Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 230000005012 migration Effects 0.000 description 2
- 238000013508 migration Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000011426 transformation method Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 238000013528 artificial neural network Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000006698 induction Effects 0.000 description 1
- 239000011229 interlayer Substances 0.000 description 1
- 239000010410 layer Substances 0.000 description 1
- 238000012417 linear regression Methods 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 230000035699 permeability Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 230000009466 transformation Effects 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/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/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- 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
- G01V1/362—Effecting static or dynamic corrections; Stacking
-
- 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
- G01V1/364—Seismic filtering
-
- 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
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
- Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping
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
技术领域
本发明涉及地震反演技术领域，具体涉及一种深度域储层地震反演方法。
背景技术
地震勘探技术中的综合分析过程经历三个发展阶段：“二维到三维、叠后到叠前、时间域到深度域”。时间域反演己经相当成熟，然而对于深度域的反演，是一个前缘课题，虽然一直有人探讨和研究，但是到目前仍然没有相对完善的理论和方法。基于褶积模型的反演，大部分停留在把时间域子波转换为深度域进行褶积，没有达到消除时深转换的最终目的。
虽然目前在深度域方面的研究还比较少，一些专门的软件也未得以开发，但近几年来国内许多学者在深度域地震数据处理方面也开始做了一定的研究。2000年张雪建等提出了深度域合成地震记录的制作方法研究；2001年林金逞等提出了应用深度域高分辨率地震反演识别低渗透薄互层储层研究；2002柴春艳等提出了深感应测井深度域反演算法及应用；2003年姚振兴等提出了用于深度域地震剖面衰减与频散补偿的反Q滤波方法；2010张静等利用多元线性回归变换方法建立波阻抗、自然伽玛、孔隙度等测井曲线与地震属性之间存在的线性变换来预测岩性和物性；2009胡中平等提出了伪深度变换的方法，这一方法有效的解决了深度域中子波随深度变化的问题；YESHPALSINGH,Repsol也是根据上述子波提取的理论在2012年对深度域反演这一课题做了较为深入的研究。
目前，深度域反演大致分为下列几类方法：1)通过数次时深转换将深度域地震记录转换到时间域再进行反演，但是几次时深转换使得有效信息大量丢失，甚至丢失反映薄层的地震信息，最终使得反演误差放大；2)忽略深度域子波的时变性，利用平均速度计算深度域子波，直接利用线性褶积模型直接进行地震记录反演，但该方法忽略深度域子波的时变性，在长时间窗内反演会造成估计的弹性参数准确度不高；3)利用滑动的窗口进行反演，在窗口内忽略子波的时变性，是方法2)的改良，使得弹性参数的准确度获得提升，但并未从根本上解决子波时变带来的反演误差大的问题；4)借助深度域偏移，将时间域合成的地震记录转换到深度域，再进行深度域反演，直接将偏移介入反演中，使得计算效率大大降低；5)利用一个等效速度将地震记录转换到满足褶积模型的“伪深度”域，进行线性褶积反演，再利用该等效速度将估计的弹性参数转换到真深度域，从而完成反演整个流程，但该方法仍然无法避免两次不同域之间的转换，使得反演误差放大；6)基于非确定性映射关系，利用统计原理和机器学习的模糊深度域反演方法，但如何提高标签数据的质量、如何扩充标签数据数量以及神经网络结构设计、效率等问题仍待解决，使得该类方法在工业界并未得到广泛应用。
发明内容
针对现有技术中的上述不足，本发明提供了一种深度域储层地震反演方法。
为了达到上述发明目的，本发明采用的技术方案为：
一种深度域储层地震反演方法，包括以下步骤：
S1、获取包含不同地震波反射特征的深度域地震数据；
S2、采用低通滤波插值方法将深度域地震数据转换为时间域地震数据；
S3、以反射系数的先验概率分布建立目标函数，以先验波阻抗约束和反射系数稀疏性约束建立约束条件，根据时间域地震数据进行目标函数求解，得到时间域地震数据对应的反射系数；
S4、根据反射系数对深度域地震数据进行反演，得到最终的地震反演结果。
进一步地，步骤S2具体包括以下步骤：
S21、采用逐道转换的方式，从需要得到的等间隔时间域地震道开始，通过时深转换速度模型计算各时间点对应的深度位置；
S22、采用低通滤波插值算子从等间隔的深度域地震道上插值得到该深度对应的振幅值。
进一步地，步骤S22中采用的低通滤波插值算子具体为：
其中，x(t)为时间域地震道上采样点在时间t的振幅值，x(nT)为等间隔时间域地震道上采样点的振幅值，n为采样序号，T为采样周期，F为采样频率。
进一步地，步骤S3中以反射系数的先验概率分布建立目标函数具体包括：
根据时间域地震数据的噪声分布将包含不同地震波反射特征的时间域地震数据进行数据融合，并基于时间域地震数据域反射系数之间的似然函数建立目标函数，表示为：
其中，J(r)为目标函数，d1,d2,d3分别为第一时间域地震数据、第二时间域地震数据和第三时间域地震数据，G1,G2,G3分别为第一时间域地震数据、第二时间域地震数据和第三时间域地震数据对应的地震子波褶积矩阵，r为反射系数，α为第二时间域地震数据的误差最小二乘拟合项的权值，β为第三时间域地震数据的误差最小二乘拟合项的权值，μ为反射系数正则化约束项的权值，ρ为先验波阻抗约束项的权值，T为转置符号，M为时间域地震数据类别数量，ri为第i类时间域地震数据对应的反射系数，δ为稀疏化参数，C为积分算子矩阵，ξ为相对波阻抗。
进一步地，步骤S3中以先验波阻抗约束和反射系数稀疏性约束建立约束条件具体包括：
以相对波阻抗与先验波阻抗之间的最小平方误差建立先验波阻抗约束，表示为：
其中，J1为先验波阻抗约束，C为积分算子矩阵，r为反射系数，ξ为相对波阻抗；
对反射系数的先验概率分布进行稀疏化，建立反射系数稀疏性约束，表示为：
其中，J2为反射系数稀疏性约束，M为时间域地震数据类别数量，ri为第i类时间域地震数据对应的反射系数，δ为稀疏化参数。
进一步地，步骤S3中根据时间域地震数据进行目标函数求解，得到时间域地震数据对应的反射系数，具体包括：
根据叠加地震速度谱模型和层位模型对时深曲线速度模型进行校正，得到校正后的空变速度模型；
根据校正后的空变速度模型的时深关系对第一时间域地震数据、第二时间域地震数据和第三时间域地震数据进行标定；
分别对第一时间域地震数据、第二时间域地震数据和第三时间域地震数据提取地震子波；
根据第一时间域地震数据、第二时间域地震数据和第三时间域地震数据对应的地震子波分别计算对应的噪声分布协方差矩阵，并转换为相对波阻抗；
根据设定的初始反射系数序列计算初始迭代矩阵向量，并通过当前迭代矩阵向量计算当前迭代反射系数；
根据相对波阻抗和当前迭代反射系数计算目标函数，直至目标函数满足设定误差阈值。
进一步地，根据叠加地震速度谱模型和层位模型对时深曲线速度模型进行校正，得到校正后的空变速度模型，具体包括：
根据叠加地震速度谱模型和层位模型对时深曲线速度模型进行校正，得到校正速度模型；
利用分层伪速度模型对校正速度模型进行二次校正，得到最终校正后的空变速度模型。
进一步地，分别对第一时间域地震数据、第二时间域地震数据和第三时间域地震数据提取地震子波，具体包括：
确定第一时间域地震数据、第二时间域地震数据和第三时间域地震数据所要提取地震子波的深度层段；
确定地震子波的采样点数范围，分别计算第一时间域地震数据、第二时间域地震数据和第三时间域地震数据的地震子波特征参数；
根据第一时间域地震数据、第二时间域地震数据和第三时间域地震数据的地震子波特征参数合成地震记录，提取得到对应的地震子波。
进一步地，确定地震子波的采样点数范围具体包括：
A1、设定初始采样区间[a1,b1]和采样精度ε，求解使得搜索函数满足以下条件的采样点数：
Γ(n)≥(b1-a1)/ε
A2、设定迭代次数k，根据采样点数计算当前迭代搜索参数：
A3、判定是否满足Γ(λk)>Γ(γk)；若是，则执行步骤A4；否则执行步骤A6；
A4、根据当前迭代搜索参数计算下次迭代搜索参数：
A5、判定迭代次数k是否满足k＝n-2；若是，则执行步骤A9；否则执行步骤A8；
A6、根据当前迭代搜索参数计算下次迭代搜索参数：
A7、判定迭代次数k是否满足k＝n-2；若是，则执行步骤A9；否则执行步骤A8；
A8、返回步骤A2将迭代次数加1；
A9、令λn＝λn-1,γn＝λn-1+1；若Γ(λn)>Γ(γn)，则将λn作为地震子波的采样点数；否则将γn作为地震子波的采样点数。
进一步地，地震子波特征参数的确定方式为：
其中，κ为地震子波特征参数，R为反射系数循环矩阵，T为转置符号，s为地震道，I为单位矩阵，tr为矩阵求迹符号，||||为范数。
本发明具有以下有益效果：
本发明通过融合包含不同地震波反射特征的深度域地震数据，并采用低通滤波插值方法将深度域地震数据转换为时间域地震数据，然后基于不同深度域地震数据与反射系数之间的似然函数确定深度域地震数据与反射系数的关联关系，最终实现深度域地震反演，能够有效提高深度域地震反演的分辨率。
附图说明
图1为本发明中一种深度域储层地震反演方法的流程示意图。
具体实施方式
下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。
如图1所述，本发明实施例提供了一种深度域储层地震反演方法，包括以下步骤S1至S4：
S1、获取包含不同地震波反射特征的深度域地震数据；
在本发明的一个可选实施例中，本实施例首先获取包含不同地震波反射特征的深度域地震数据，以相同地下地质条件在不同观测尺度下得到的三维地面地震数据、三维垂直地震剖面数据和井间地震数据来表征不同地震波反射特征的深度域地震数据，这三种不同地震波反射特征对应相同的反射系数。
本实施例通过井中接收的三维垂直地震剖面数据能够比较精确的建立地震资料与测井资料之间的时深关系，增强地面地震反演的可靠性；井间地震数据分辨率较高，具备将高精度井旁信息向空间外推的能力；两种井中地震资料都具备了介于地面地震和测井之间的研究尺度，不仅具有一定的空间连续性而且具有较高的纵向分辨率，是有效的将测井和地面地震资料分辨率衔接起来的高分辨率地球物理技术。因此本实施例充分利用包含不同地震波反射特征的深度域地震数据，减小地震反演问题固有的多解性、提高分辨率，以更好的解决地下地质问题。
S2、采用低通滤波插值方法将深度域地震数据转换为时间域地震数据；
在本发明的一个可选实施例中，步骤S2具体包括以下步骤：
S21、采用逐道转换的方式，从需要得到的等间隔时间域地震道开始，通过时深转换速度模型计算各时间点对应的深度位置；
S22、采用低通滤波插值算子从等间隔的深度域地震道上插值得到该深度对应的振幅值。
具体而言，将井间地震资料转换到时间域的优势在于可以借鉴时间域地面地震资料中已有的认识来建立井间地震的反射特征，并且可以采用时间域比较成熟的反演方法来进行处理。在进行深时转换时需要建立合理的速度模型，井间地震观测深度段之上的深时转换速度模型可以采用所在探区的平均速度，当横向速度变化比较平缓时采用多井标定的时深关系来拟合出平均速度则能够更好的保证时间域井间地震反射特征与地面地震资料相匹配；而井间地震观测深度段的速度模型则采用与井间地震剖面相对应的层析速度，并用声波测井速度所确定的时深关系对速度模型进行校正。在此过程中需要对层析速度资料边界处的奇异值进行剔除和平滑处理，以避免由于速度奇异点和层间不光滑而造成地震反射同相轴畸变和锯齿状等现象。
由于深度域井间地震剖面与期望输出的时间域剖面都要求等间隔采样，因此在采用速度模型对井间资料进行深时转换时不可避免的需要对中间过程产生的不等间隔地震道进行插值。用常规线性算子进行插值会影响原始信号的波形及频谱特征，因此采用精确的插值算子来获取等时间间隔上的合理地震振幅值是实现深时转换中波形保真的关键。根据采样定理，重采样后的信号可以表示成采样前信号同理想低通滤波器的卷积，表示为
其中，x(t)为时间域地震道上采样点在时间t的振幅值，x(nT)为等间隔时间域地震道上采样点的振幅值，n为采样序号，T为采样周期，F为采样频率。
S3、以反射系数的先验概率分布建立目标函数，以先验波阻抗约束和反射系数稀疏性约束建立约束条件，根据时间域地震数据进行目标函数求解，得到时间域地震数据对应的反射系数；
在本发明的一个可选实施例中，步骤S3中以反射系数的先验概率分布建立目标函数具体包括：
根据时间域地震数据的噪声分布将包含不同地震波反射特征的时间域地震数据进行数据融合，并基于时间域地震数据域反射系数之间的似然函数建立目标函数，表示为：
其中，J(r)为目标函数，d1,d2,d3分别为第一时间域地震数据、第二时间域地震数据和第三时间域地震数据，G1,G2,G3分别为第一时间域地震数据、第二时间域地震数据和第三时间域地震数据对应的地震子波褶积矩阵，r为反射系数，α为第二时间域地震数据的误差最小二乘拟合项的权值，β为第三时间域地震数据的误差最小二乘拟合项的权值，μ为反射系数正则化约束项的权值，ρ为先验波阻抗约束项的权值，T为转置符号，M为时间域地震数据类别数量，ri为第i类时间域地震数据对应的反射系数，δ为稀疏化参数，C为积分算子矩阵，ξ为相对波阻抗。
具体而言，由于包含不同地震波反射特征的深度域地震数据的采集、处理相互独立，数据间的差异主要表现在分辨尺度和观测噪声的不同，在贝叶斯反演理论中似然函数可以通过噪声分布特征来表示，而不同深度域地震数据之间的相互独立性确保了不同深度域地震数据d与反射系数r之间的似然函数可以通过联合概率分布来建立，因此本实施例根据时间域地震数据的噪声分布将包含不同地震波反射特征的时间域地震数据进行数据融合，并基于时间域地震数据域反射系数之间的似然函数建立目标函数，表示为：
本实施例的步骤S3中以先验波阻抗约束和反射系数稀疏性约束建立约束条件具体包括：
以相对波阻抗与先验波阻抗之间的最小平方误差建立先验波阻抗约束，表示为：
其中，J1为先验波阻抗约束，C为积分算子矩阵，r为反射系数，ξ为相对波阻抗；
对反射系数的先验概率分布进行稀疏化，建立反射系数稀疏性约束，表示为：
其中，J2为反射系数稀疏性约束，M为时间域地震数据类别数量，ri为第i类时间域地震数据对应的反射系数，δ为稀疏化参数。
本实施例的步骤S3中根据时间域地震数据进行目标函数求解，得到时间域地震数据对应的反射系数，具体包括：
根据叠加地震速度谱模型和层位模型对时深曲线速度模型进行校正，得到校正后的空变速度模型；
根据校正后的空变速度模型的时深关系对第一时间域地震数据、第二时间域地震数据和第三时间域地震数据进行标定；
分别对第一时间域地震数据、第二时间域地震数据和第三时间域地震数据提取地震子波；
根据第一时间域地震数据、第二时间域地震数据和第三时间域地震数据对应的地震子波分别计算对应的噪声分布协方差矩阵，并转换为相对波阻抗；
根据设定的初始反射系数序列计算初始迭代矩阵向量，并通过当前迭代矩阵向量计算当前迭代反射系数；
根据相对波阻抗和当前迭代反射系数计算目标函数，直至目标函数满足设定误差阈值。
本实施例根据叠加地震速度谱模型和层位模型对时深曲线速度模型进行校正，得到校正后的空变速度模型，具体包括：
根据叠加地震速度谱模型和层位模型对时深曲线速度模型进行校正，得到校正速度模型；
利用分层伪速度模型对校正速度模型进行二次校正，得到最终校正后的空变速度模型。
本实施例分别对第一时间域地震数据、第二时间域地震数据和第三时间域地震数据提取地震子波，具体包括：
确定第一时间域地震数据、第二时间域地震数据和第三时间域地震数据所要提取地震子波的深度层段；
确定地震子波的采样点数范围，分别计算第一时间域地震数据、第二时间域地震数据和第三时间域地震数据的地震子波特征参数；
根据第一时间域地震数据、第二时间域地震数据和第三时间域地震数据的地震子波特征参数合成地震记录，提取得到对应的地震子波。
其中确定地震子波的采样点数范围具体包括：
A1、设定初始采样区间[a1,b1]和采样精度ε，求解使得搜索函数满足以下条件的采样点数：
Γ(n)≥(b1-a1)/ε
A2、设定迭代次数k，根据采样点数计算当前迭代搜索参数：
A3、判定是否满足Γ(λk)>Γ(γk)；若是，则执行步骤A4；否则执行步骤A6；
A4、根据当前迭代搜索参数计算下次迭代搜索参数：
A5、判定迭代次数k是否满足k＝n-2；若是，则执行步骤A9；否则执行步骤A8；
A6、根据当前迭代搜索参数计算下次迭代搜索参数：
A7、判定迭代次数k是否满足k＝n-2；若是，则执行步骤A9；否则执行步骤A8；
A8、返回步骤A2将迭代次数加1；
A9、令λn＝λn-1,γn＝λn-1+1；若Γ(λn)>Γ(γn)，则将λn作为地震子波的采样点数；否则将γn作为地震子波的采样点数。
其中地震子波特征参数的确定方式为：
其中，κ为地震子波特征参数，R为反射系数循环矩阵，T为转置符号，s为地震道，I为单位矩阵，tr为矩阵求迹符号，||||为范数。
S4、根据反射系数对深度域地震数据进行反演，得到最终的地震反演结果。
在本发明的一个可选实施例中，本实施例根据深度域地震数据与反射系数的关联关系，对深度域地震数据进行反演，得到最终的地震反演结果。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。
Claims (10)
1.一种深度域储层地震反演方法，其特征在于，包括以下步骤：
S1、获取包含不同地震波反射特征的深度域地震数据；
S2、采用低通滤波插值方法将深度域地震数据转换为时间域地震数据；
S3、以反射系数的先验概率分布建立目标函数，以先验波阻抗约束和反射系数稀疏性约束建立约束条件，根据时间域地震数据进行目标函数求解，得到时间域地震数据对应的反射系数；
S4、根据反射系数对深度域地震数据进行反演，得到最终的地震反演结果。
2.根据权利要求1所述的一种深度域储层地震反演方法，其特征在于，步骤S2具体包括以下步骤：
S21、采用逐道转换的方式，从需要得到的等间隔时间域地震道开始，通过时深转换速度模型计算各时间点对应的深度位置；
S22、采用低通滤波插值算子从等间隔的深度域地震道上插值得到该深度对应的振幅值。
3.根据权利要求2所述的一种深度域储层地震反演方法，其特征在于，步骤S22中采用的低通滤波插值算子具体为：
其中，x(t)为时间域地震道上采样点在时间t的振幅值，x(nT)为等间隔时间域地震道上采样点的振幅值，n为采样序号，T为采样周期，F为采样频率。
4.根据权利要求1所述的一种深度域储层地震反演方法，其特征在于，步骤S3中以反射系数的先验概率分布建立目标函数具体包括：
根据时间域地震数据的噪声分布将包含不同地震波反射特征的时间域地震数据进行数据融合，并基于时间域地震数据域反射系数之间的似然函数建立目标函数，表示为：
其中，J(r)为目标函数，d1,d2,d3分别为第一时间域地震数据、第二时间域地震数据和第三时间域地震数据，G1,G2,G3分别为第一时间域地震数据、第二时间域地震数据和第三时间域地震数据对应的地震子波褶积矩阵，r为反射系数，α为第二时间域地震数据的误差最小二乘拟合项的权值，β为第三时间域地震数据的误差最小二乘拟合项的权值，μ为反射系数正则化约束项的权值，ρ为先验波阻抗约束项的权值，T为转置符号，M为时间域地震数据类别数量，ri为第i类时间域地震数据对应的反射系数，δ为稀疏化参数，C为积分算子矩阵，ξ为相对波阻抗。
5.根据权利要求1所述的一种深度域储层地震反演方法，其特征在于，步骤S3中以先验波阻抗约束和反射系数稀疏性约束建立约束条件具体包括：
以相对波阻抗与先验波阻抗之间的最小平方误差建立先验波阻抗约束，表示为：
其中，J1为先验波阻抗约束，C为积分算子矩阵，r为反射系数，ξ为相对波阻抗；
对反射系数的先验概率分布进行稀疏化，建立反射系数稀疏性约束，表示为：
其中，J2为反射系数稀疏性约束，M为时间域地震数据类别数量，ri为第i类时间域地震数据对应的反射系数，δ为稀疏化参数。
6.根据权利要求1所述的一种深度域储层地震反演方法，其特征在于，步骤S3中根据时间域地震数据进行目标函数求解，得到时间域地震数据对应的反射系数，具体包括：
根据叠加地震速度谱模型和层位模型对时深曲线速度模型进行校正，得到校正后的空变速度模型；
根据校正后的空变速度模型的时深关系对第一时间域地震数据、第二时间域地震数据和第三时间域地震数据进行标定；
分别对第一时间域地震数据、第二时间域地震数据和第三时间域地震数据提取地震子波；
根据第一时间域地震数据、第二时间域地震数据和第三时间域地震数据对应的地震子波分别计算对应的噪声分布协方差矩阵，并转换为相对波阻抗；
根据设定的初始反射系数序列计算初始迭代矩阵向量，并通过当前迭代矩阵向量计算当前迭代反射系数；
根据相对波阻抗和当前迭代反射系数计算目标函数，直至目标函数满足设定误差阈值。
7.根据权利要求6所述的一种深度域储层地震反演方法，其特征在于，根据叠加地震速度谱模型和层位模型对时深曲线速度模型进行校正，得到校正后的空变速度模型，具体包括：
根据叠加地震速度谱模型和层位模型对时深曲线速度模型进行校正，得到校正速度模型；
利用分层伪速度模型对校正速度模型进行二次校正，得到最终校正后的空变速度模型。
8.根据权利要求6所述的一种深度域储层地震反演方法，其特征在于，分别对第一时间域地震数据、第二时间域地震数据和第三时间域地震数据提取地震子波，具体包括：
确定第一时间域地震数据、第二时间域地震数据和第三时间域地震数据所要提取地震子波的深度层段；
确定地震子波的采样点数范围，分别计算第一时间域地震数据、第二时间域地震数据和第三时间域地震数据的地震子波特征参数；
根据第一时间域地震数据、第二时间域地震数据和第三时间域地震数据的地震子波特征参数合成地震记录，提取得到对应的地震子波。
9.根据权利要求8所述的一种深度域储层地震反演方法，其特征在于，确定地震子波的采样点数范围具体包括：
A1、设定初始采样区间[a1,b1]和采样精度ε，求解使得搜索函数满足以下条件的采样点数：
Γ(n)≥(b1-a1)/ε
A2、设定迭代次数k，根据采样点数计算当前迭代搜索参数：
A3、判定是否满足Γ(λk)>Γ(γk)；若是，则执行步骤A4；否则执行步骤A6；
A4、根据当前迭代搜索参数计算下次迭代搜索参数：
A5、判定迭代次数k是否满足k＝n-2；若是，则执行步骤A9；否则执行步骤A8；
A6、根据当前迭代搜索参数计算下次迭代搜索参数：
A7、判定迭代次数k是否满足k＝n-2；若是，则执行步骤A9；否则执行步骤A8；
A8、返回步骤A2将迭代次数加1；
A9、令λn＝λn-1,γn＝λn-1+1；若Γ(λn)>Γ(γn)，则将λn作为地震子波的采样点数；否则将γn作为地震子波的采样点数。
10.根据权利要求1所述的一种深度域储层地震反演方法，其特征在于，地震子波特征参数的确定方式为：
其中，κ为地震子波特征参数，R为反射系数循环矩阵，T为转置符号，s为地震道，I为单位矩阵，tr为矩阵求迹符号，||为范数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311183468.4A CN117310802B (zh) | 2023-09-13 | 2023-09-13 | 一种深度域储层地震反演方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311183468.4A CN117310802B (zh) | 2023-09-13 | 2023-09-13 | 一种深度域储层地震反演方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117310802A true CN117310802A (zh) | 2023-12-29 |
| CN117310802B CN117310802B (zh) | 2024-06-07 |
Family
ID=89296351
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311183468.4A Active CN117310802B (zh) | 2023-09-13 | 2023-09-13 | 一种深度域储层地震反演方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117310802B (zh) |
Citations (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20060291328A1 (en) * | 2004-05-21 | 2006-12-28 | Robertsson Johan O A | Interpolation and extrapolation method for seismic recordings and use of same in multiple supression |
| CN104769458A (zh) * | 2014-07-15 | 2015-07-08 | 杨顺伟 | 一种基于柯西分布的叠后波阻抗反演方法 |
| CN106094027A (zh) * | 2016-06-01 | 2016-11-09 | 中国海洋石油总公司 | 一种垂直地震剖面vsp钻前压力预测方法和系统 |
| CN107193040A (zh) * | 2017-06-27 | 2017-09-22 | 中国石油天然气股份有限公司 | 深度域合成地震记录的确定方法和装置 |
| CN107329171A (zh) * | 2017-06-07 | 2017-11-07 | 中国石油天然气股份有限公司 | 深度域储层地震反演方法及装置 |
| WO2018201114A1 (en) * | 2017-04-28 | 2018-11-01 | Pioneer Natural Resources Usa, Inc. | High resolution seismic data derived from pre-stack inversion and machine learning |
| CN110954953A (zh) * | 2019-12-16 | 2020-04-03 | 中国地质大学（武汉） | 一种基于柯西分布的叠后波阻抗反演方法 |
| CN111077571A (zh) * | 2019-12-12 | 2020-04-28 | 成都理工大学 | 一种提高分辨率的孔隙度反演方法 |
| CN111208561A (zh) * | 2020-01-07 | 2020-05-29 | 自然资源部第一海洋研究所 | 基于时变子波与曲波变换约束的地震声波阻抗反演方法 |
| CN111948712A (zh) * | 2020-08-10 | 2020-11-17 | 中海石油(中国)有限公司 | 一种基于深度域地震记录的叠前线性反演方法 |
| CN113917540A (zh) * | 2021-11-10 | 2022-01-11 | 同济大学 | 基于稀疏约束的抗假频射线束进行地震数据去噪的方法 |
| WO2023123971A1 (zh) * | 2021-12-30 | 2023-07-06 | 中国石油天然气集团有限公司 | 基于vsp的深度域地震剖面层位标定方法及装置 |
-
2023
- 2023-09-13 CN CN202311183468.4A patent/CN117310802B/zh active Active
Patent Citations (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20060291328A1 (en) * | 2004-05-21 | 2006-12-28 | Robertsson Johan O A | Interpolation and extrapolation method for seismic recordings and use of same in multiple supression |
| CN104769458A (zh) * | 2014-07-15 | 2015-07-08 | 杨顺伟 | 一种基于柯西分布的叠后波阻抗反演方法 |
| CN106094027A (zh) * | 2016-06-01 | 2016-11-09 | 中国海洋石油总公司 | 一种垂直地震剖面vsp钻前压力预测方法和系统 |
| WO2018201114A1 (en) * | 2017-04-28 | 2018-11-01 | Pioneer Natural Resources Usa, Inc. | High resolution seismic data derived from pre-stack inversion and machine learning |
| CN107329171A (zh) * | 2017-06-07 | 2017-11-07 | 中国石油天然气股份有限公司 | 深度域储层地震反演方法及装置 |
| CN107193040A (zh) * | 2017-06-27 | 2017-09-22 | 中国石油天然气股份有限公司 | 深度域合成地震记录的确定方法和装置 |
| CN111077571A (zh) * | 2019-12-12 | 2020-04-28 | 成都理工大学 | 一种提高分辨率的孔隙度反演方法 |
| CN110954953A (zh) * | 2019-12-16 | 2020-04-03 | 中国地质大学（武汉） | 一种基于柯西分布的叠后波阻抗反演方法 |
| CN111208561A (zh) * | 2020-01-07 | 2020-05-29 | 自然资源部第一海洋研究所 | 基于时变子波与曲波变换约束的地震声波阻抗反演方法 |
| CN111948712A (zh) * | 2020-08-10 | 2020-11-17 | 中海石油(中国)有限公司 | 一种基于深度域地震记录的叠前线性反演方法 |
| CN113917540A (zh) * | 2021-11-10 | 2022-01-11 | 同济大学 | 基于稀疏约束的抗假频射线束进行地震数据去噪的方法 |
| WO2023123971A1 (zh) * | 2021-12-30 | 2023-07-06 | 中国石油天然气集团有限公司 | 基于vsp的深度域地震剖面层位标定方法及装置 |
Non-Patent Citations (4)
| Title |
|---|
| LAMBERT M ET AL: "Born-type schemes for the acoustic probing of 1-D fluid media from time-harmonic planar reflection coefficients at two incidences", JOURNAL OF THE ACOUSTICAL SOCIETY OF AMERICA, vol. 99, no. 1, 1 January 1996 (1996-01-01), pages 243 - 253 * |
| LAURENCE LETKI;张荣忠;: "时间域和深度域反演声阻抗对比", 油气地球物理, no. 02, 26 April 2016 (2016-04-26) * |
| 姚振岸: "深度域地震资料反演方法研究", 中国博士学位论文全文数据库, no. 2022, 15 July 2020 (2020-07-15) * |
| 陈珂磷 等: "页岩裂缝型储层模型参数化及AVAZ 反演预测方法研究", 地球物理学进展, vol. 37, no. 6, 5 September 2022 (2022-09-05), pages 2364 - 237 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117310802B (zh) | 2024-06-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106405651B (zh) | 一种基于测井匹配的全波形反演初始速度模型构建方法 | |
| US20090119076A1 (en) | Method for Generating a 3D Earth Model | |
| CA2616379C (en) | Method for tomographic inversion by matrix transformation | |
| CN109541685B (zh) | 一种河道砂体识别方法 | |
| CN101201409B (zh) | 一种地震数据变相位校正方法 | |
| CN110895348B (zh) | 一种地震弹性阻抗低频信息提取方法、系统及存储介质 | |
| WO2005026776A1 (en) | Wide-offset-range pre-stack depth migration method for seismic exploration | |
| CN108957554B (zh) | 一种地球物理勘探中的地震反演方法 | |
| CN117310802B (zh) | 一种深度域储层地震反演方法 | |
| CN106842297A (zh) | 井约束非稳态相位校正方法 | |
| CN106842291B (zh) | 一种基于叠前地震射线阻抗反演的不整合圈闭储层岩性预测方法 | |
| CN112213776B (zh) | 叠前道集和vsp资料联合层控q模型建立方法 | |
| CN109061737B (zh) | 一种基于合成地震记录的储层预测方法及装置 | |
| CN110398773B (zh) | 一种针对部分缺失的地震数据的恢复重构方法 | |
| Li et al. | Automatic extraction of seismic data horizon across faults | |
| CN110857997A (zh) | 一种基于横向约束的分步叠前弹性参数反演方法及系统 | |
| CN113960658B (zh) | 一种基于地质地震模型的测井约束速度建模方法及装置 | |
| CN112415580B (zh) | 消除速度模型突变界面的方法、叠前深度偏移处理方法 | |
| CN112558147B (zh) | 一种井中微地震数据的偏振分析方法及系统 | |
| CN113009579B (zh) | 地震数据反演方法及装置 | |
| CN107561580A (zh) | 基于多层多井建立初始地质模型的方法 | |
| Dutta et al. | Joint 4D full-waveform inversion using enhanced template-matching objective function | |
| CN116520409A (zh) | 油藏薄小储层高分辨率地震反演识别方法 | |
| CN116699678A (zh) | 一种多信息融合的速度建模方法及装置 | |
| CN116859454A (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 |