CN108535775B - 非平稳地震资料声波阻抗反演方法及装置 - Google Patents
非平稳地震资料声波阻抗反演方法及装置 Download PDFInfo
- Publication number
- CN108535775B CN108535775B CN201810299727.2A CN201810299727A CN108535775B CN 108535775 B CN108535775 B CN 108535775B CN 201810299727 A CN201810299727 A CN 201810299727A CN 108535775 B CN108535775 B CN 108535775B
- Authority
- CN
- China
- Prior art keywords
- seismic
- wave impedance
- impedance
- iteration
- given
- 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 71
- 239000011159 matrix material Substances 0.000 claims abstract description 65
- 230000006870 function Effects 0.000 claims description 58
- 230000003190 augmentative effect Effects 0.000 claims description 21
- 230000015572 biosynthetic process Effects 0.000 claims description 19
- 238000004590 computer program Methods 0.000 claims description 13
- 238000005070 sampling Methods 0.000 claims description 12
- 238000004422 calculation algorithm Methods 0.000 claims description 11
- 230000002238 attenuated effect Effects 0.000 claims description 7
- 230000001131 transforming effect Effects 0.000 claims description 6
- 230000017105 transposition Effects 0.000 claims description 5
- 230000003595 spectral effect Effects 0.000 claims description 4
- 230000003416 augmentation Effects 0.000 abstract 2
- 238000005755 formation reaction Methods 0.000 description 17
- 238000010586 diagram Methods 0.000 description 12
- 238000012545 processing Methods 0.000 description 12
- 230000008569 process Effects 0.000 description 10
- 238000003860 storage Methods 0.000 description 9
- 230000000694 effects Effects 0.000 description 8
- 238000001914 filtration Methods 0.000 description 7
- 238000005516 engineering process Methods 0.000 description 4
- 230000008859 change Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 230000008094 contradictory effect Effects 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 230000002708 enhancing effect Effects 0.000 description 1
- 230000007274 generation of a signal involved in cell-cell signaling Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 239000013067 intermediate product Substances 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 230000005055 memory storage Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 239000000047 product Substances 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 230000003068 static effect 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/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
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/622—Velocity, density or impedance
- G01V2210/6226—Impedance
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
本申请实施例提供了一种非平稳地震资料声波阻抗反演方法及装置，该方法包括：确定叠后衰减地震数据的初始地震子波、地层常Q模型和初始波阻抗模型；根据所述初始地震子波、所述地层常Q模型和所述初始波阻抗模型生成地震子波褶积矩阵，并根据所述地震子波褶积矩阵建立反演声波阻抗的目标函数；将所述目标函数变换为增广拉格朗日形式的目标函数；对于给定地震道数据，通过交替方向乘子法求解所述增广拉格朗日形式的目标函数，从而获得声波阻抗；根据求解出的声波阻抗确定给定地震道数据的绝对波阻抗。本申请实施例可以获得更为准确的非平稳地震资料的波阻抗反演结果。
Description
技术领域
本申请涉及石油地球物理勘探地震参数反演技术领域，尤其是涉及一种非平稳地震资料声波阻抗反演方法及装置。
背景技术
地震解释通常需要对地下地层的弹性、物性信息进行定性甚至定量的获取，以便达到更为精细的储层刻画描述。而地震解释人员通常希望获得的地震资料不仅仅限于高分辨率地震剖面，相应直观准确的波阻抗数据体可作为地层属性用于定量地震解释。因此，通过反射波地震资料得到波阻抗数据同样是现阶段地震资料处理技术所涵盖的范围。如何通过接收到的地震数据获取波阻抗信息也一直在为了克服各种问题而持续得到深入的研究。
现阶段针对地震数据波阻抗反演主要集中在以反射系数为间接变量，通过地震道积分和递推反演，得到相对波阻抗信息。但由于地震信号频带的缺失导致估计的波阻抗并不稳定，通过低波数的剔除可以得到相对合理的反演结果。但是相对阻抗并不能完全反映大部分地质构造信息。并且噪声的加入使得得到地震频段外的阻抗信息可信度降低。因此通过添加合理的约束或先验信息，达到地震数据直接反演波阻抗的技术是现阶段地震地质解释工作中更为认可的地层弹性参数获取手段。
然而，由于地震子波的带限滤波效应及地层Q滤波效应的存在，地震数据必须经过一系列高分辨率处理技术流程从而向理论模型靠近，而这样的处理方法在消除相应滤波效应的同时引入了新问题，尤其是不稳定现象及噪声振幅同时补偿；即处理技术中新加入的干扰使得原本具有很强多解性的波阻抗反演可能出现不合理的结果。另一方面，现有的振幅补偿和反Q滤波对输入吸收衰减参数的偏差影响无法做到自动识别及调节，残余的衰减效应在地震资料中依然能够引起波阻抗反演结果的偏差。
发明内容
本申请实施例的目的在于提供一种非平稳地震资料声波阻抗反演方法及装置，以获得更为准确的非平稳地震资料的波阻抗反演结果。
为达到上述目的，一方面，本申请实施例提供了一种非平稳地震资料声波阻抗反演方法，包括：
确定叠后衰减地震数据的初始地震子波、地层常Q模型和初始波阻抗模型；
根据所述初始地震子波、所述地层常Q模型和所述初始波阻抗模型生成地震子波褶积矩阵，并根据所述地震子波褶积矩阵建立反演声波阻抗的目标函数；
将所述目标函数变换为增广拉格朗日形式的目标函数；
对于给定地震道数据，通过交替方向乘子法求解所述增广拉格朗日形式的目标函数，从而获得声波阻抗；
根据求解出的声波阻抗确定给定地震道数据的绝对波阻抗。
本申请实施例的非平稳地震资料声波阻抗反演方法中，所述根据所述初始地震子波、所述地层常Q模型和所述初始波阻抗模型生成地震子波褶积矩阵，并根据所述地震子波褶积矩阵建立反演声波阻抗的目标函数，包括：
根据公式确定地震子波褶积矩阵；
根据公式建立反演声波阻抗的目标函数；
其中，G为地震子波褶积矩阵，M为频率域采样点数，F-1为反傅里叶变换矩阵，W为以初始地震子波的频谱值作为对角线元素得到的矩阵，A为衰减矩阵，D为差分矩阵，zr为地震声波阻抗，s为衰减地震信号且s＝Gzr+n，n为噪声，μ为正则化参数，参数p的取值范围为[0,1]。
本申请实施例的非平稳地震资料声波阻抗反演方法中，所述对于给定地震道数据，通过交替方向乘子法求解所述增广拉格朗日形式的目标函数，包括：
将给定地震道数据输入以下公式以进行迭代求解，直至达到预设的迭代停止条件；
其中，为第k次迭代得到的zr，为第k+1次迭代得到的zr，r为地层反射系数，rk为第k次迭代得到的r，rk+1为第k+1次迭代得到的r，λ为拉格朗日乘子项，λT为λ的转置，λk为第k次迭代得到的λ，λk+1第k+1次迭代得到的λ，ρ为拉格朗日乘子项更新步长，zr的迭代更新通过L-BFGS算法实现，r的迭代更新通过广义迭代阈值算法实现。
本申请实施例的非平稳地震资料声波阻抗反演方法中，所述迭代停止条件为：
j≥itmax或
其中，j为迭代次数，itmax为给定迭代次数阈值，zr为地震声波阻抗，为第j次迭代得到的zr，为第j-1次迭代得到的zr，ε为给定模型更新误差阈值。
本申请实施例的非平稳地震资料声波阻抗反演方法中，所述根据求解出的声波阻抗确定给定地震道数据的绝对波阻抗，包括：
根据公式z＝exp(2×zr)确定给定地震道数据的绝对波阻抗；
其中，z为给定地震道数据的绝对波阻抗，zr为地震声波阻抗。
另一方面，本申请实施例还提供了一种非平稳地震资料声波阻抗反演装置，包括存储器、处理器、以及存储在所述存储器上的计算机程序，所述计算机程序被所述处理器运行时执行如下步骤：
确定叠后衰减地震数据的初始地震子波、地层常Q模型和初始波阻抗模型；
根据所述初始地震子波、所述地层常Q模型和所述初始波阻抗模型生成地震子波褶积矩阵，并根据所述地震子波褶积矩阵建立反演声波阻抗的目标函数；
将所述目标函数变换为增广拉格朗日形式的目标函数；
对于给定地震道数据，通过交替方向乘子法求解所述增广拉格朗日形式的目标函数，从而获得声波阻抗；
根据求解出的声波阻抗确定给定地震道数据的绝对波阻抗。
本申请实施例的非平稳地震资料声波阻抗反演装置中，所述根据所述初始地震子波、所述地层常Q模型和所述初始波阻抗模型生成地震子波褶积矩阵，并根据所述地震子波褶积矩阵建立反演声波阻抗的目标函数，包括：
根据公式确定地震子波褶积矩阵；
根据公式建立反演声波阻抗的目标函数；
其中，G为地震子波褶积矩阵，M为频率域采样点数，F-1为反傅里叶变换矩阵，W为以初始地震子波的频谱值作为对角线元素得到的矩阵，A为衰减矩阵，D为差分矩阵，zr为地震声波阻抗，s为衰减地震信号且s＝Gzr+n，n为噪声，μ为正则化参数，参数p的取值范围为[0,1]。
本申请实施例的非平稳地震资料声波阻抗反演装置中，所述对于给定地震道数据，通过交替方向乘子法求解所述增广拉格朗日形式的目标函数，包括：
将给定地震道数据输入以下公式以进行迭代求解，直至达到预设的迭代停止条件；
其中，为第k次迭代得到的zr，为第k+1次迭代得到的zr，r为地层反射系数，rk为第k次迭代得到的r，rk+1为第k+1次迭代得到的r，λ为拉格朗日乘子项，λT为λ的转置，λk为第k次迭代得到的λ，λk+1第k+1次迭代得到的λ，ρ为拉格朗日乘子项更新步长，zr的迭代更新通过L-BFGS算法实现，r的迭代更新通过广义迭代阈值算法实现。
本申请实施例的非平稳地震资料声波阻抗反演装置中，所述迭代停止条件为：
j≥itmax或
其中，j为迭代次数，itmax为给定迭代次数阈值，zr为地震声波阻抗，为第j次迭代得到的zr，为第j-1次迭代得到的zr，ε为给定模型更新误差阈值。
本申请实施例的非平稳地震资料声波阻抗反演装置中，所述根据求解出的声波阻抗确定给定地震道数据的绝对波阻抗，包括：
根据公式z＝exp(2×zr)确定给定地震道数据的绝对波阻抗；
其中，z为给定地震道数据的绝对波阻抗，zr为地震声波阻抗。
由以上本申请实施例提供的技术方案可见，本申请实施例中，在确定叠后衰减地震数据的初始地震子波、地层常Q模型和初始波阻抗模型的基础上，根据初始地震子波、地层常Q模型和初始波阻抗模型生成地震子波褶积矩阵，并根据地震子波褶积矩阵建立反演声波阻抗的目标函数；然后将目标函数变换为增广拉格朗日形式的目标函数；对于给定地震道数据，通过交替方向乘子法求解增广拉格朗日形式的目标函数，从而获得声波阻抗；最后根据求解出的声波阻抗确定给定地震道数据的绝对波阻抗。由于该技术方案通过地震数据与波阻抗之间的映射关系建立反演方程，直接求取绝对波阻抗，抛弃了现阶段普遍使用反射系数进行波阻抗反演的策略；而且，该方法是针对叠后衰减地震资料进行的，避免了处理流程中反Q滤波方法引入新的干扰对反演结果的影响。从而使得本申请实施例可得到更为准确的非平稳地震资料的绝对波阻抗结果。
附图说明
为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。在附图中：
图1为本申请一实施例中非平稳地震资料声波阻抗反演方法的流程图；
图2a为本申请一实施例中一真实Marmousi波阻抗模型；
图2b为本申请一实施例中用于正演合成衰减的零相位Ricker地震子波；
图2c为本申请一实施例中通过Ricker地震子波、Marmouis真实波阻抗以及常Q模型合成衰减地震记录；
图2d为本申请一实施例中反演波阻抗使用的初始模型；
图2e为本申请一实施例中基于交替方向乘子法获得的叠后地震声波阻抗反演结果；
图2f为本申请一实施例中是CDP号为240处真实波阻抗，初始波阻抗和反演波阻抗的对比图；
图2g为本申请一实施例中是CDP号为800处真实波阻抗，初始波阻抗和反演波阻抗的对比图；
图3a为本申请一实施例中某油田勘探区域局部叠后地震数据；
图3b为本申请一实施例中根据测井资料中的波阻抗曲线得到剖面的初始波阻抗模型；
图3c为采用本申请实施例的非平稳地震资料声波阻抗反演方法反演得到的绝对波阻抗剖面；
图3d为本申请一实施例中CDP号为401处声波测井资料中波阻抗曲线、初始波阻抗和反演波阻抗的对比图；
图4为本申请一实施例中非平稳地震资料声波阻抗反演方法的结构框图。
具体实施方式
为了使本技术领域的人员更好地理解本申请中的技术方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都应当属于本申请保护的范围。例如在下面描述中，在第一部件上方形成第二部件，可以包括第一部件和第二部件以直接接触方式形成的实施例，还可以包括第一部件和第二部件以非直接接触方式(即第一部件和第二部件之间还可以包括额外的部件)形成的实施例等。
而且，为了便于描述，本申请一些实施例可以使用诸如“在…上方”、“在…之下”、“顶部”、“下方”等空间相对术语，以描述如实施例各附图所示的一个元件或部件与另一个(或另一些)元件或部件之间的关系。应当理解的是，除了附图中描述的方位之外，空间相对术语还旨在包括装置在使用或操作中的不同方位。例如若附图中的装置被翻转，则被描述为“在”其他元件或部件“下方”或“之下”的元件或部件，随后将被定位为“在”其他元件或部件“上方”或“之上”。
对于反射波地震资料，常规的褶积模型认为地震子波在传播过程中是稳定的，即子波的振幅与相位不随传播时间的增加而发生改变。最终接收到的为稳态地震反射信号。然而，本申请的发明人研究发现，实际中地震子波的传播由于地层对不同频率谐波衰减程度不同，导致地震子波的振幅和相位发生变化。该过程可以通过地层品质因子(Q)模型描述，当考虑地层Q滤波效应时，反射地震信号具有非平稳性质。而反射波地震信号生成过程又可以通过地震子波和反射系数的褶积运算进行刻画。同时该反射系数可以通过地层波阻抗进行运算结果。因此，非平稳地震数据可以通过绝对波阻抗和初始地震子波及Q模型进行相关运算得到，具体数学表示为：
其中，s为衰减地震信号，M为频率域采样点数，F-1为反傅里叶(Fourior)变换矩阵，A为衰减矩阵，W为以初始地震子波的频谱值作为对角线元素得到的矩阵，zr为地震声波阻抗在进行简单对数运算得到的向量表示，n为噪声，D为差分矩阵。矩阵F-1和A具体表示为：
ωi为角频率采样点，ti为时间采样点，Q为地层品质因子，频率域采样点个数M、时间域采样点个数N、及采样间隔Δt需满足采样定理。D为差分矩阵，具体形式为：
衰减矩阵A中元素值的大小同时与地层品质因子(即Q值)有关。通过非平稳叠后反射波地震数据的模型刻画及数学描述，可以建立直接反演地层声波阻抗的反演方程，即求解公式(1)中的zr。而差分矩阵D和声波阻抗zr的乘积即为地层反射系数r，运算结果即为不同时间位置的地震子波褶积矩阵。为了表述方便，公式上述公式(1)可以重新写为：
s＝Gzr+n(2)
其中，G为地震子波褶积矩阵，且由于衰减性质及地震子波的带限性质，衰减矩阵A具有低秩特征，并且差分矩阵D存在不稳定性。因此对于地震子波褶积矩阵G的逆几乎无法求取，因此基于方程(2)求解zr较为困难。现阶段针对稳态地震资料的声波阻抗反演集中在引入其他关于阻抗的约束条件，或者加入稳定因子，也可实现波阻抗反演。然而，衰减因素的存在使得加入稳定因子的策略不再适用。反演方程的病态性和绝对波阻抗的反演需求相互矛盾。因此需要更为先进的技术对波阻抗进行精确的估计。
为了从具有衰减特征的叠后反射波地震资料中准确地反演得到声波阻抗剖面，协助进一步的地震资料解释工作，并且在噪声水平较高的情况下仍然能够得到相对准确的反演结果。因此，本申请实施例从上述公式(2)出发，提出了具有较高稳定性的声波阻抗反演方案。具体如下：
由于地震子波褶积矩阵G的低秩性和不可逆性，通过公式(2)直接求解绝对波阻抗会极不稳定。并且针对波阻抗的假设条件一般很难满足。因此，本申请实施例针对反演波阻抗的中间产物施加约束，间接协助减小反演的多解性和增强其稳定性。在一般的高分辨率处理中，反射系数的先验分布通常情况下能够符合地质统计学规律。现阶段基于稀疏约束的反射系数反演在一定程度上能够获得良好的效果，因此，本申请实施例可通过将反射系数的lp范数项引入反演方程，建立关于绝对波阻抗的目标函数，通过求取目标函数的极小值最终确定绝对波阻抗。
对于求解声波阻抗zr的目标函数具体可表示为：
其中，||Dzr||p为反射系数的lp范数，μ为正则化参数，参数p的取值范围为[0,1]。对正则化参数和迭代步长的确定例如可以通过L-curve曲线捕获，或者通过局部数据试算得到。由于lp约束项的存在以及矩阵G和D的不可逆性，公式(3)难以通过部署简单的线性方法获取最终的极小值。因此，在本申请一实施例中，可通过引入交替方向乘子法对zr进行求解。该方法的主要思想是对于目标函数进行分离求解并引入自变量和中间变量的乘子项实现算法的收敛。为此，首先需要对目标函数进行一定的修正，即目标函数的增广拉格朗日(Lagrangian)形式为：
其中，λ为拉格朗日乘子项，ρ为拉格朗日乘子项更新步长。对于新的目标函数，通过交替方向乘子法进行迭代求解，最终算法迭代的形式为：
在每次算法迭代过程中更新相应的变量，首先对于zr的更新，该求解函数为严格的凸函数，因此本申请一实施例中，可采用L-BFGS(Limited-memory BFGS)算法求解；对于r的更新，由于lp范数的存在，传统的凸优化方法并不能用于求解该项，因此，本申请一实施例中，可采用广义迭代阈值算法求取相应的变量值；对于拉格朗日乘子项的更新最为简单，直接求取即可。
在得到zr后，可通过公式以下公式(6)计算得到最终的声波阻抗数值。
z＝exp(2×zr) (6)
基于以上原理说明，参考图1所示，本申请实施例的非平稳地震资料声波阻抗反演方法可以包括以下步骤：
S101、确定叠后衰减地震数据的初始地震子波、地层常Q模型和初始波阻抗模型；
在本申请一实施例中，初始波阻抗模型可以通过测井资料的声波阻抗曲线的低频成分得到，或者通过资料处理过程中的速度模型得到；根据地层常Q模型可以确定地层品质因子。
S102、根据所述初始地震子波、所述地层常Q模型和所述初始波阻抗模型生成地震子波褶积矩阵，并根据所述地震子波褶积矩阵建立反演声波阻抗的目标函数；
在本申请一实施例中，所述根据初始地震子波、地层常Q模型和初始波阻抗模型生成地震子波褶积矩阵，并根据所述地震子波褶积矩阵建立反演声波阻抗的目标函数，例如可以是首先根据公式确定地震子波褶积矩阵，然后，根据公式(3)建立反演声波阻抗的目标函数。
S103、将所述目标函数变换为增广拉格朗日形式的目标函数。
S104、对于给定地震道数据，通过交替方向乘子法求解所述增广拉格朗日形式的目标函数，从而获得声波阻抗；
在本申请一实施例中，所述对于给定地震道数据，通过交替方向乘子法求解所述增广拉格朗日形式的目标函数，从而获得声波阻抗，可以是将给定地震道数据输入上述公式(5)以进行迭代求解，直至达到预设的迭代停止条件；
其中，为第k次迭代得到的zr，为第k+1次迭代得到的zr，rk为第k次迭代得到的r，rk+1为第k+1次迭代得到的r，λT为λ的转置，λk为第k次迭代得到的λ，λk+1第k+1次迭代得到的λ。
在本申请一实施例中，所述迭代停止条件可以为：j≥itmax或否则，若j＜itmax或则j＝j+1，即继续循环迭代。
其中，j为迭代次数，itmax为给定迭代次数阈值，zr为地震声波阻抗，为第j次迭代得到的zr，为第j-1次迭代得到的zr，ε为给定模型更新误差阈值。
S105、根据求解出的声波阻抗确定给定地震道数据的绝对波阻抗。
在本申请一实施例中，可将求解出的声波阻抗输入上述公式(6)，从而可以计算出给定地震道数据的绝对波阻抗。
由此而可见，本申请实施例的非平稳地震资料声波阻抗反演方法中，以叠后衰减地震资料为输入数据，在增广拉格朗日框架下，通过分裂目标函数及引入乘子项策略，采用LBGFS算法及广义迭代阈值算法分别对变量进行求解，最终达到估计地层声波阻抗的目的，该方法通过地震数据与波阻抗之间的映射关系建立反演方程，直接求取绝对波阻抗，抛弃了现阶段普遍使用反射系数进行波阻抗反演的策略。而且，该方法是针对叠后衰减地震资料进行的，避免了处理流程中反Q滤波方法引入新的干扰对反演结果的影响，最终得到更为准确的绝对波阻抗结果。
下面介绍本申请一示例性实施例。图2a和图2b分别示出一真实Marmousi波阻抗模型及用于正演合成衰减的零相位Ricker地震子波。在图2a中，时间长度为2442ms，共1360个CDP(common depth point，共深度点)炮集；在图2b中，Ricker子波主频为30Hz，采样间隔1ms，采样点数为75。通过Ricker地震子波、Marmouis真实波阻抗以及常Q模型(Q＝50)可合成如图2c所示的衰减地震记录，衰减地震记录中加入高斯随机噪声，信噪比为15dB。但由于衰减现象的存在，深部信噪比低于15dB，因此反演较为困难。
在给定的如图2d所示的初始波阻抗模型基础上，通过实施本申请实施例的非平稳地震资料声波阻抗反演方法，可得到如图2e所示的叠后地震声波阻抗反演结果。图2f和图2g分别示出了CDP号为240、800处真实波阻抗(见图2f和图2g中的True所对应的曲线)，初始波阻抗(见图2f和图2g中的Initial所对应的曲线)和反演波阻抗(见图2f和图2g中的inverted所对应的曲线)的对比图。对比可以看出，本申请实施例的非平稳地震资料声波阻抗反演方法具有较高的准确性。
图3a示出了国内某油田勘探区域局部叠后地震数据，在给定的如图3b所示的初始波阻抗模型基础上，通过实施本申请实施例的非平稳地震资料声波阻抗反演方法，可得到如图3c所示的绝对波阻抗剖面。而图3d示出了CDP号为401处声波测井资料中波阻抗曲线(见图3d中的Log所对应的曲线)、初始波阻抗(见图3d中的Initial所对应的曲线)和反演波阻抗(见图3d中的inverted所对应的曲线)的对比图。由此可以看出，在应用实际时，采用本申请实施例的非平稳地震资料声波阻抗反演方法，同样可以获得具有较高的准确性的绝对波阻抗剖面。从而进一步验证本申请实施例的非平稳地震资料声波阻抗反演方法的技术效果。
参考图4所示，本申请实施例的一种非平稳地震资料声波阻抗反演装置，包括存储器、处理器、以及存储在所述存储器上的计算机程序，所述计算机程序被所述处理器运行时执行如下步骤：
确定叠后衰减地震数据的初始地震子波、地层常Q模型和初始波阻抗模型；
根据所述初始地震子波、所述地层常Q模型和所述初始波阻抗模型生成地震子波褶积矩阵，并根据所述地震子波褶积矩阵建立反演声波阻抗的目标函数；
将所述目标函数变换为增广拉格朗日形式的目标函数；
对于给定地震道数据，通过交替方向乘子法求解所述增广拉格朗日形式的目标函数，从而获得声波阻抗；
根据求解出的声波阻抗确定给定地震道数据的绝对波阻抗。
虽然上文描述的过程流程包括以特定顺序出现的多个操作，但是，应当清楚了解，这些过程可以包括更多或更少的操作，这些操作可以顺序执行或并行执行(例如使用并行处理器或多线程环境)。
为了描述的方便，描述以上装置时以功能分为各种单元分别描述。当然，在实施本申请时可以把各单元的功能在同一个或多个软件和/或硬件中实现。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
在一个典型的配置中，计算设备包括一个或多个处理器(CPU)、输入/输出接口、网络接口和内存。
内存可能包括计算机可读介质中的非永久性存储器，随机存取存储器(RAM)和/或非易失性内存等形式，如只读存储器(ROM)或闪存(flash RAM)。内存是计算机可读介质的示例。
计算机可读介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机的存储介质的例子包括，但不限于相变内存(PRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁带磁磁盘存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算设备访问的信息。按照本文中的界定，计算机可读介质不包括暂存电脑可读媒体(transitory media)，如调制的数据信号和载波。
还需要说明的是，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法或者设备中还存在另外的相同要素。
本领域技术人员应明白，本申请的实施例可提供为方法、系统或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本申请，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。
本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
以上所述仅为本申请的实施例而已，并不用于限制本申请。对于本领域技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本申请的权利要求范围之内。
Claims (6)
1.一种非平稳地震资料声波阻抗反演方法，其特征在于，包括：
确定叠后衰减地震数据的初始地震子波、地层常Q模型和初始波阻抗模型；
根据所述初始地震子波、所述地层常Q模型和所述初始波阻抗模型生成地震子波褶积矩阵，并根据所述地震子波褶积矩阵建立反演声波阻抗的目标函数；
将所述目标函数变换为增广拉格朗日形式的目标函数；
对于给定地震道数据，通过交替方向乘子法求解所述增广拉格朗日形式的目标函数，从而获得声波阻抗；
根据求解出的声波阻抗确定给定地震道数据的绝对波阻抗；其中，
所述根据所述初始地震子波、所述地层常Q模型和所述初始波阻抗模型生成地震子波褶积矩阵，并根据所述地震子波褶积矩阵建立反演声波阻抗的目标函数，包括：
根据公式确定地震子波褶积矩阵；
根据公式建立反演声波阻抗的目标函数；
其中，G为地震子波褶积矩阵，M为频率域采样点数，F-1为反傅里叶变换矩阵，W为以初始地震子波的频谱值作为对角线元素得到的矩阵，A为衰减矩阵，D为差分矩阵，zr为地震声波阻抗，s为衰减地震信号且s＝Gzr+n，n为噪声，μ为正则化参数，参数p的取值范围为[0,1]；
所述对于给定地震道数据，通过交替方向乘子法求解所述增广拉格朗日形式的目标函数，包括：
将给定地震道数据输入以下公式以进行迭代求解，直至达到预设的迭代停止条件；
其中，为第k次迭代得到的zr，为第k+1次迭代得到的zr，r为地层反射系数，rk为第k次迭代得到的r，rk+1为第k+1次迭代得到的r，λ为拉格朗日乘子项，λT为λ的转置，λk为第k次迭代得到的λ，λk+1第k+1次迭代得到的λ，ρ为拉格朗日乘子项更新步长，zr的迭代更新通过L-BFGS算法实现，r的迭代更新通过广义迭代阈值算法实现。
2.如权利要求1所述的非平稳地震资料声波阻抗反演方法，其特征在于，所述迭代停止条件为：
j≥itmax或
其中，j为迭代次数，itmax为给定迭代次数阈值，zr为地震声波阻抗，为第j次迭代得到的zr，为第j-1次迭代得到的zr，ε为给定模型更新误差阈值。
3.如权利要求1所述的非平稳地震资料声波阻抗反演方法，其特征在于，所述根据求解出的声波阻抗确定给定地震道数据的绝对波阻抗，包括：
根据公式z＝exp(2×zr)确定给定地震道数据的绝对波阻抗；
其中，z为给定地震道数据的绝对波阻抗，zr为地震声波阻抗。
4.一种非平稳地震资料声波阻抗反演装置，包括存储器、处理器、以及存储在所述存储器上的计算机程序，其特征在于，所述计算机程序被所述处理器运行时执行如下步骤：
确定叠后衰减地震数据的初始地震子波、地层常Q模型和初始波阻抗模型；
根据所述初始地震子波、所述地层常Q模型和所述初始波阻抗模型生成地震子波褶积矩阵，并根据所述地震子波褶积矩阵建立反演声波阻抗的目标函数；
将所述目标函数变换为增广拉格朗日形式的目标函数；
对于给定地震道数据，通过交替方向乘子法求解所述增广拉格朗日形式的目标函数，从而获得声波阻抗；
根据求解出的声波阻抗确定给定地震道数据的绝对波阻抗；其中，
所述根据所述初始地震子波、所述地层常Q模型和所述初始波阻抗模型生成地震子波褶积矩阵，并根据所述地震子波褶积矩阵建立反演声波阻抗的目标函数，包括：
根据公式确定地震子波褶积矩阵；
根据公式建立反演声波阻抗的目标函数；
其中，G为地震子波褶积矩阵，M为频率域采样点数，F-1为反傅里叶变换矩阵，W为以初始地震子波的频谱值作为对角线元素得到的矩阵，A为衰减矩阵，D为差分矩阵，zr为地震声波阻抗，s为衰减地震信号且s＝Gzr+n，n为噪声，μ为正则化参数，参数p的取值范围为[0,1]；
所述对于给定地震道数据，通过交替方向乘子法求解所述增广拉格朗日形式的目标函数，包括：
将给定地震道数据输入以下公式以进行迭代求解，直至达到预设的迭代停止条件；
其中，为第k次迭代得到的zr，为第k+1次迭代得到的zr，r为地层反射系数，rk为第k次迭代得到的r，rk+1为第k+1次迭代得到的r，λ为拉格朗日乘子项，λT为λ的转置，λk为第k次迭代得到的λ，λk+1第k+1次迭代得到的λ，ρ为拉格朗日乘子项更新步长，zr的迭代更新通过L-BFGS算法实现，r的迭代更新通过广义迭代阈值算法实现。
5.如权利要求4所述的非平稳地震资料声波阻抗反演装置，其特征在于，所述迭代停止条件为：
j≥itmax或
其中，j为迭代次数，itmax为给定迭代次数阈值，zr为地震声波阻抗，为第j次迭代得到的zr，为第j-1次迭代得到的zr，ε为给定模型更新误差阈值。
6.如权利要求4所述的非平稳地震资料声波阻抗反演装置，其特征在于，所述根据求解出的声波阻抗确定给定地震道数据的绝对波阻抗，包括：
根据公式z＝exp(2×zr)确定给定地震道数据的绝对波阻抗；
其中，z为给定地震道数据的绝对波阻抗，zr为地震声波阻抗。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810299727.2A CN108535775B (zh) | 2018-03-30 | 2018-03-30 | 非平稳地震资料声波阻抗反演方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810299727.2A CN108535775B (zh) | 2018-03-30 | 2018-03-30 | 非平稳地震资料声波阻抗反演方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108535775A CN108535775A (zh) | 2018-09-14 |
| CN108535775B true CN108535775B (zh) | 2019-08-23 |
Family
ID=63483185
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810299727.2A Active CN108535775B (zh) | 2018-03-30 | 2018-03-30 | 非平稳地震资料声波阻抗反演方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108535775B (zh) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110471113B (zh) * | 2019-08-01 | 2020-08-04 | 中国石油大学(北京) | 基于非稳态地震资料的反演动校正方法、装置及存储介质 |
| CN112415587B (zh) * | 2019-08-21 | 2024-06-18 | 中国石油化工股份有限公司 | 储层地震波衰减特性分析方法及储层反射系数的反演方法 |
| CN112649848B (zh) * | 2019-10-12 | 2024-01-23 | 中国石油化工股份有限公司 | 利用波动方程求解地震波阻抗的方法和装置 |
| CN110865409B (zh) * | 2019-12-02 | 2021-08-31 | 怀化学院 | 一种基于波阻抗低秩正则化的地震波阻抗反演方法 |
| CN112630835B (zh) * | 2020-12-03 | 2023-10-17 | 重庆三峡学院 | 一种高分辨率的叠后地震波阻抗反演的方法 |
| CN112596107B (zh) * | 2020-12-11 | 2021-09-28 | 中国科学院地质与地球物理研究所 | 弹性参数反演方法及装置 |
| CN113341463B (zh) * | 2021-06-10 | 2023-05-26 | 中国石油大学(北京) | 一种叠前地震资料非平稳盲反褶积方法及相关组件 |
| CN113640871B (zh) * | 2021-08-10 | 2023-09-01 | 成都理工大学 | 一种基于重加权l1范数稀疏约束的地震波阻抗反演方法 |
| CN113589386B (zh) * | 2021-09-15 | 2022-06-10 | 中国石油大学(北京) | 一种基于对比函数的块状声波阻抗反演方法、装置及设备 |
| CN114994758B (zh) * | 2022-08-02 | 2022-10-28 | 北京京鲁聚源能源科技有限公司 | 碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104769458A (zh) * | 2014-07-15 | 2015-07-08 | 杨顺伟 | 一种基于柯西分布的叠后波阻抗反演方法 |
| CN106443775A (zh) * | 2016-05-25 | 2017-02-22 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 高分辨率转换波裂缝预测方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US10031009B2 (en) * | 2011-08-23 | 2018-07-24 | Cidra Corporate Services, Inc. | Flow profiling techniques based on modulated magnetic-electrical impedance tomography |
-
2018
- 2018-03-30 CN CN201810299727.2A patent/CN108535775B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104769458A (zh) * | 2014-07-15 | 2015-07-08 | 杨顺伟 | 一种基于柯西分布的叠后波阻抗反演方法 |
| CN106443775A (zh) * | 2016-05-25 | 2017-02-22 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 高分辨率转换波裂缝预测方法 |
Non-Patent Citations (1)
| Title |
|---|
| 井间地震约束下的高分辨率波阻抗反演方法研究;曹丹平 等;《石油物探》;20100930;第49卷(第5期);第425-429页 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108535775A (zh) | 2018-09-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108535775B (zh) | 非平稳地震资料声波阻抗反演方法及装置 | |
| Zu et al. | A periodically varying code for improving deblending of simultaneous sources in marine acquisition | |
| Wang | Inverse Q-filter for seismic resolution enhancement | |
| US9470811B2 (en) | Creating a high resolution velocity model using seismic tomography and impedance inversion | |
| US8923093B2 (en) | Determining the quality of a seismic inversion | |
| Fomel | Local seismic attributes | |
| Yuan et al. | Stable inversion-based multitrace deabsorption method for spatial continuity preservation and weak signal compensation | |
| Edgar et al. | How reliable is statistical wavelet estimation? | |
| Wang et al. | Enhancing resolution of nonstationary seismic data by molecular-Gabor transform | |
| CN107894612B (zh) | 一种q吸收衰减补偿的声波阻抗反演方法及系统 | |
| CN107132579B (zh) | 一种保地层结构的地震波衰减补偿方法 | |
| CN110471113B (zh) | 基于非稳态地震资料的反演动校正方法、装置及存储介质 | |
| Yang et al. | Viscoacoustic least-squares reverse time migration using a time-domain complex-valued wave equation | |
| US10310117B2 (en) | Efficient seismic attribute gather generation with data synthesis and expectation method | |
| Choi et al. | Time-domain full-waveform inversion of exponentially damped wavefield using the deconvolution-based objective function | |
| CN112882099B (zh) | 一种地震频带拓宽方法、装置、介质及电子设备 | |
| CN110687597B (zh) | 一种基于联合字典的波阻抗反演方法 | |
| CN112711065A (zh) | 叠前地震反演方法及装置 | |
| CN108680957B (zh) | 基于加权的局部互相关时频域相位反演方法 | |
| Oliveira et al. | Nonlinear impedance inversion for attenuating media | |
| Wu et al. | Q-factor estimation in CMP gather and the continuous spectral ratio slope method | |
| Ma et al. | Random noise attenuation by fx spatial projection-based complex empirical mode decomposition predictive filtering | |
| CN113341463B (zh) | 一种叠前地震资料非平稳盲反褶积方法及相关组件 | |
| CN112363222A (zh) | 叠后自适应宽带约束波阻抗反演方法及装置 | |
| Zhao et al. | Absorption-constrained wavelet power spectrum inversion for enhancing resolution of nonstationary seismic data |
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 |