CN107894612A - 一种q吸收衰减补偿的声波阻抗反演方法及系统 - Google Patents
一种q吸收衰减补偿的声波阻抗反演方法及系统 Download PDFInfo
- Publication number
- CN107894612A CN107894612A CN201710994606.5A CN201710994606A CN107894612A CN 107894612 A CN107894612 A CN 107894612A CN 201710994606 A CN201710994606 A CN 201710994606A CN 107894612 A CN107894612 A CN 107894612A
- Authority
- CN
- China
- Prior art keywords
- absorption
- attenuations
- compensation
- data
- impedance
- 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
- 238000010521 absorption reaction Methods 0.000 title claims abstract description 75
- 238000000034 method Methods 0.000 title claims abstract description 51
- 230000010354 integration Effects 0.000 claims abstract description 9
- 239000011159 matrix material Substances 0.000 claims description 21
- 238000001228 spectrum Methods 0.000 claims description 15
- 238000005070 sampling Methods 0.000 claims description 6
- 230000017105 transposition Effects 0.000 claims description 5
- 235000013399 edible fruits Nutrition 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 4
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 230000000694 effects Effects 0.000 abstract description 23
- 238000001914 filtration Methods 0.000 abstract description 20
- 239000006185 dispersion Substances 0.000 abstract description 13
- 238000009825 accumulation Methods 0.000 abstract description 4
- 230000015572 biosynthetic process Effects 0.000 description 9
- 238000005755 formation reaction Methods 0.000 description 5
- 238000012545 processing Methods 0.000 description 5
- 230000006854 communication Effects 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 239000000203 mixture Substances 0.000 description 4
- 238000003786 synthesis reaction Methods 0.000 description 4
- 230000000644 propagated effect Effects 0.000 description 3
- 238000002310 reflectometry Methods 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 230000007423 decrease Effects 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 241001269238 Data Species 0.000 description 1
- 208000027697 autoimmune lymphoproliferative syndrome due to CTLA4 haploinsuffiency Diseases 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000004615 ingredient Substances 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000007787 solid Substances 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/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
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吸收衰减频散因子集成到褶积模型中，通过迭代反演的方法，从吸收衰减的非稳态地震记录数据中，直接反演得到高分辨率阻抗结果，本发明方法避免了传统递推过程中存在的误差累积效应，避免了传统反Q滤波方法存在的数值不稳定问题，也即本发明Q值补偿技术是无条件稳定的，其可补偿比传统技术更高的频率成分，进一步提高地震勘探数据的分辨率。
Description
技术领域
本发明属于地震勘探技术中地震信号处理领域，具体涉及一种Q吸收衰减补偿的声波阻抗反演方法及系统。
背景技术
叠后波阻抗反演技术是油气地震勘探数据处理过程中较为重要的一步，技术的输入是叠后地震数据和该数据对应的地震子波，技术处理的输出为反演的地层的阻抗，也即速度和密度的乘积。阻抗的差异是油气检测的重要方法之一。上下地层的阻抗差异决定了地层分界面处的反射系数。地震勘探过程中，波动在地层中进行传播，由于地层介质的非完全弹性(也即粘弹性)，波在传播过程中发生吸收衰减与频散。地震波的吸收衰减指的是，波在地层传播过程中振幅逐渐变弱，一部分波的能量转换为热能损耗掉。高频成分相比低频成分吸收衰减的多。地震波的频散指的是波在地下介质传播过程中，不同频率成分的波传播速度不同，也即速度频散现象。地层的粘弹性特征改变了在地层中传播的波的形态和特征，衰减了地震波的能量，拉伸和改变了在地层中传播的波，降低了地震勘探的分辨率。地球物理学家引入品质因子Q来刻画地层的粘弹性滤波效应，Q值的定义为Q＝2πE/△E，其中E为地震波的能量，△E为地震波传播一个波长长度时能量的损失。
为了进一步阐述本发明所解决的科学问题，这里给出图1(a)-(d)，假定有图1(a)地下地层阻抗模型，如果波在地下传播过程中没有发生吸收衰减与频散，利用褶积模型，可以模拟得到图1(b)所示的不含Q吸收衰减与频散效应的数据。可以看出浅中深信号的能量较为一致，反映了地层分界面的反射系数信息。如果波在地下传播过程中发生吸收衰减与频散，利用非稳态Q吸收衰减频散褶积模型，可以模拟得到图1(c)(图1(c)中细粗线不含Q吸收衰减效应，黑粗线含Q吸收衰减效应)所示的经过Q吸收衰减与频散效应的数据。可以看出浅层信号能量较强，中深信号的能量较弱。如果直接从图1(c)中数据中进行反演，得到的将是如图1(d)所示的阻抗结果，相比真实声波阻抗模型图1(a)，可以看出此时反演的声波阻抗结果分辨率较低，Q吸收衰减滤波效应较为明显。由于此时反演并未妥善考虑地层的Q滤波效应，因此反演得到的阻抗结果分辨率较低。
为了恢复补偿地层的Q滤波效应，提高地震数据的分辨率，地球物理学家提出了反Q滤波技术。反Q滤波包括振幅补偿和相位补偿。其中，反Q滤波的相位补偿是无条件稳定的，因为相位补偿仅是对信号Fourier频谱的相位进行校正，根据Fourier变换的时移特性，相位校正在时间域体现为信号位置上的前后变化，相位补偿并不改变信号的振幅能量。而反Q滤波的振幅补偿则是存在固有的数值不稳定问题。为了便于阐述本发明旨在解决的第二个科学问题，这里给出了图2(a)-(d)，其中图2(a)代表了Q吸收衰减滤波后得到的信号，可以看出波的能量随着传播时间的增大，波的振幅上存在衰减，波形上也逐渐的不再对称，即波形存在相位频散。图2(b)示意了反Q滤波振幅补偿算子的数值变化，注意图中是以量级的形式进行显示。从图2(b)可以看出，反Q滤波振幅补偿算子增长迅速，成指数型增长，很快增长至10的20次幂以上，进而反Q补偿后得到的信号如图2(c)所示，这里采用了真实的Q值，可以清晰地看出反Q滤波存在的数值不稳定问题。当然为了改进反Q滤波的数值不稳定问题，许多研究人员提出了不同的技术方法，多数采取的策略是限制了反Q补偿的高频成分，补偿结果如图2(d)所示，可以看出增益控制反Q滤波对强烈Q吸收衰减的信号没能得到恰当的补偿，高频信号成分未得到恢复，也即传统增益控制反Q滤波方法并不能补偿强烈Q衰减的高频成分，提高地震数据分辨率的能力有限。
现有技术的技术方案
为了克服和避免传统反Q滤波存在的数值不稳定问题。Chai于2014年提出一种稳健的高分辨率非稳态稀疏反射系数反演方法(Nonstationary Sparse ReflectivityInversion，NSRI)。非稳态稀疏反射系数反演方法NSRI通过将Q吸收衰减因子集成到传统稳态褶积模型中，对待求解的反射系数序列施加L1范数稀疏约束，并通过基追踪迭代求解的算法，巧妙地避免了反Q滤波存在的数值不稳定问题。为了更好地阐释NSRI方法和本发明旨在解决的关键问题，这里给出图3，其中图3(a)为图1(a)所对应的真实反射系数模型，图3(b)是NSRI方法从图1(c)中直接反演得到的反射系数。可以看出NSRI方法可以直接从Q吸收衰减的数据中反演得到高分辨率的反射系数序列，而巧妙地避免了反Q滤波振幅补偿存在的数值不稳定问题。
现有技术的缺点
对比图3(b)和图3(a)，可以看出NSRI反演的反射系数和真实的反射系数结果近乎一样，但是由反演的反射系数和初始阻抗模型来递推得到声波阻抗结果时，不可避免地产生了递推公式误差累积放大效应。从图3(c)递推得到的声波阻抗结果可以看出，由反射系数递推得到的阻抗结果出现了“挂面条”、“条带状”现象，同真实阻抗模型相比存在较大的误差。这就是传统技术存在的问题。
发明内容
为了克服上述现有技术的不足，本发明提供了一种Q吸收衰减补偿的声波阻抗反演方法及系统，其旨在从如图1(c)所示的含Q吸收衰减效应的地震数据中直接反演得到高分辨的声波阻抗结果(如图3(d)所示)，避免反Q滤波存在的固有数值不稳定问题，避免递推公式所造成的误差累积放大效应。
根据本发明的其中一方面，本发明为解决其技术问题，提供了一种Q吸收衰减补偿的声波阻抗反演方法，包含如下步骤：
(1)获取数据s、W、F-1、FQ以及D；
(2)采用LSQR算法或下述第一算法求解模型s＝WF-1FQDx得到x，求解时正演算子A＝WF-1FQD，观测数据向量为b＝s；
(3)根据sDQ＝WDx得到反Q补偿后的地震记录数据sDQ作为反演结果；
第一算法：
a)获取下述参数：正演算子A＝WF-1FQD，观测数据向量s，x的初始值向量x0，最大迭代次数lmax；
b)初始化：迭代变量l＝0,将x0赋值给y0；
c)重复进行下述各d至f步骤lmax次；
d)计算上标H代表共轭转置；
e)将xl+1更新为yl+1；
f)将l更新为l+1；
g)得出解x＝xlmax；
其中，W是地震子波所构建的褶积矩阵，代表了反Fourier变换矩阵，F-1∈CN×N，N代表了时间离散采样点数，是品质因子Q所构建的吸收衰减矩阵，FQ∈CN×N，N×N表示行数和列数均为N的矩阵，D中其他元素为0，s是实际获取的地震记录s(t)的向量表达式， 为输入地震子波w(t)的Fourier变换频谱，ω＝2πf为角频率，r(z)为在深度z处的反射系数，α(ω,τz,Qz)为Q吸收衰减因子，e为指数函数，t代表时间，τz为传播旅行时，其定义为x是与深度有关的积分变量，vr是与参考频率ωr有关的相速度，Qeff(z)是一种等效Q值，其定义为：
在本发明的Q吸收衰减补偿的声波阻抗反演方法中，还包括根据求解阻抗矩阵的步骤。
在本发明的Q吸收衰减补偿的声波阻抗反演方法中，步骤(2)中在采用LSQR算法或下述第一算法进行求解时，x的初始值向量x0全部设为零值或者为先验背景场对应的值。
根据本发明的另一方面，本发明为解决其技术问题，还提供了一种Q吸收衰减补偿的声波阻抗反演系统，包含如下模块：
数据获取模块，用于获取数据s、W、F-1、FQ以及D；
算法求解模块，用于采用LSQR算法或下述第一算法求解模型s＝WF-1FQDx得到x，求解时正演算子A＝WF-1FQD，观测数据向量为b＝s；
地震记录计算模块，用于根据sDQ＝WDx得到反Q补偿后的地震记录数据sDQ作为反演结果；
第一算法：
a)获取下述参数：正演算子A＝WF-1FQD，观测数据向量s，x的初始值向量x0，最大迭代次数lmax；
b)初始化：迭代变量l＝0,将x0赋值给y0；
c)重复进行下述各d至f步骤lmax次；
d)计算上标H代表共轭转置；
e)将xl+1更新为yl+1；
f)将l更新为l+1；
g)得出解x＝xlmax；
其中，W是地震子波所构建的褶积矩阵，代表了反Fourier变换矩阵，F-1∈CN×N，N代表了时间离散采样点数，是品质因子Q所构建的吸收衰减矩阵，FQ∈CN×N，N×N表示行数和列数均为N的矩阵，D中其他元素为0，s是实际获取的地震记录s(t)的向量表达式， 为输入地震子波w(t)的Fourier变换频谱，ω＝2πf为角频率，r(z)为在深度z处的反射系数，α(ω,tz,Qz)为Q吸收衰减因子，e为指数函数，t代表时间，τz为传播旅行时，其定义为x是与深度有关的积分变量，vr是与参考频率ωr有关的相速度，Qeff(z)是一种等效Q值，其定义为：
在本发明的Q吸收衰减补偿的声波阻抗反演系统中，还包括阻抗求解模块，用于根据求解阻抗矩阵
在本发明的Q吸收衰减补偿的声波阻抗反演系统中，算法求解模块中在采用LSQR算法或下述第一算法进行求解时，x的初始值向量x0全部设为零值或者为先验背景场对应的值。
本发明具体公开了一种Q吸收衰减补偿的声波阻抗反演方法，阻抗反演是地震油气资源勘探的重要方法之一，本发明将地层Q吸收衰减频散因子集成到褶积模型中，通过迭代反演的方法，从吸收衰减的非稳态地震记录数据中，直接反演得到高分辨率阻抗结果，本发明方法避免了传统递推过程中存在的误差累积效应，避免了传统反Q滤波方法存在的数值不稳定问题，也即本发明Q值补偿技术是无条件稳定的；其可补偿比传统技术更高的频率成分，进一步提高了地震勘探数据的分辨率。
附图说明
下面将结合附图及实施例对本发明作进一步说明，附图中：
图1(a)为Marmousi2模型的真实声波阻抗模型图；
图1(b)为由Marmousi2声波阻抗模型合成的不含Q吸收衰减的地震数据图；
图1(c)为由Marmousi2声波阻抗模型合成的含Q吸收衰减的地震数据图；
图1(d)为从图1(c)中的含Q吸收衰减效应的地震数据中直接反演到的声波阻抗图；
图2(a)为用常Q＝25合成的包含Q吸收衰减效应的地震记录数据图；
图2(b)为进行反Q滤波时振幅补偿算子随频率f和时间t值的变化图；
图2(c)为对图2(a)中地震数据进行精确的反Q滤波的结果图；
图2(d)为对图2(a)中地震数据进行增益控制的反Q滤波图；
图3(a)从Marmousi2真实声波阻抗模型导出的真实反射系数模型图；
图3(b)利用NSRI方法从图1(c)中的合成的含Q吸收衰减的地震数据反演的反射系数结果图；
图3(c)为从图3(b)中NSRI反演的反射系数采用递归公式恢复的声波阻抗结果图；
图3(d)为Marmousi2模型的真实声波阻抗模型图；
图4本发明的第一算法的流程图；
图5(a)为Marmousi2模型的真实声波阻抗模型图；
图5(b)为本发明从图1(c)中含Q吸收衰减效应的地震数据直接反演到的相对声波阻抗图；
图5(c)代表一个初始声波阻抗模型；
图5(d)为本发明从图1(c)含Q吸收衰减效应的地震数据中以图5(c)中声波阻抗模型为初始值反演到的绝对声波阻抗结果图；
图6(a)为由Marmousi2声波阻抗模型合成的不含Q吸收衰减的地震数据图；
图6(b)由图5(b)经本发明反演的相对声波阻抗模型重构的Q吸收衰减补偿的数据图；
图6(c)为图6(a)和图6(b)的比对图；
图7是图6(c)中数据频率域的振幅谱图；
图8(a)为Colony实际地震数据图；
图8(b)为本发明采用一个无穷大Q模型时反演的相对阻抗结果图；
图8(c)为本发明采用一个尝试性Q模型时反演的相对阻抗结果图；
图8(d)由图8(c)本发明反演的相对声波阻抗模型重构的Q吸收衰减补偿的数据图；
图9中画出了对此Colony实际地震数据进行Q吸收衰减补偿前后的频率域振幅谱对比图；
图10(a)为Erskine三维实际地震数据中的一条线的数据图；
图10(b)为商业软件基于模型的阻抗反演模块反演的绝对声波阻抗结果图；
图10(c)为本发明采用一个无穷大Q模型时反演的相对阻抗结果图；
图10(d)为本发明采用一个尝试性Q模型时反演的相对阻抗结果图；
图10(e)为在反演过程中最终的数据匹配误差图；
图10(f)为本发明进行Q吸收衰减补偿后的地震数据图；
图11(a)为整个Erskine三维实际地震数据图；
图11(b)为本发明对于整个Erskine三维实际地震数据采用一个无穷大Q模型时反演的相对阻抗结果图；
图11(c)为本发明对于整个Erskine三维实际地震数据采用一个尝试性Q模型时反演的相对阻抗结果图；
图11(d)为本发明对于整个Erskine三维实际地震数据进行Q吸收衰减补偿后的结果图。
名词解释
地震品质因子Q：英文为Quality factor，用于表述波动在地下介质传播过程中所发生的吸收衰减与频散效应，其降低了地震勘探采集数据的分辨率；
Q吸收衰减补偿：是指在地震勘探数据处理过程中，恢复补偿地层品质因子Q值所引起的吸收衰减与频散效应，旨在提高地震勘探数据处理的分辨率；
声波阻抗(Impedance)：在地震勘探数据处理领域中，阻抗I指的是地下地层介质速度v和密度ρ的乘积，也即I＝vρ；
反射系数：由于地层介质属性的差异，波动在地下介质传播过程中会发生反射，两地层分界面处的反射系数r定义为：v2、ρ2分别为反射界面处上层地层介质速度和密度，v1、ρ1分别为反射界面处下层地层介质速度和密度；
LSQR：Least-squares QR，最小二乘QR分解算法。
具体实施方式
为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。其中下述大写加粗的字母表示向量或者矩阵。
模拟Q滤波后的地震数据
基于一维波动方程，模拟Q吸收衰减滤波后的地震数据s(t)的公式为：
其中，为输入地震子波w(t)的Fourier变换频谱，ω＝2πf为角频率，r(z)为在深度z处的反射系数，α(ω,τz,Qz)为Q吸收衰减因子，e为指数函数，t代表时间，τz为传播旅行时，其定义为：
x是与深度有关的积分变量，vr是与参考频率ωr有关的相速度。根据常Q模型理论，也即假定Q值在地震数据频带范围内是不随频率变化的。根据Kolsky-Futterman常Q模型，吸收衰减因子α(ω,τz,Qz)定义表达式为：
其中，Qeff(z)是一种等效Q值，其定义为：
方程1可以进一步离散，改写成矩阵向量相乘的形式：
s＝WF-1FQr (5)
其中，s和r是地震记录s(t)和反射系数序列r(t)的向量表达式，W是地震子波所构建的褶积矩阵，代表了反Fourier变换矩阵，N代表了时间离散采样点数，是品质因子Q所构建的吸收衰减矩阵，FQ∈CN ×N。当Q→+∞，α(ω,τ,Q)→1，FQ演变为正Fourier变换矩阵且F∈CN×N。
反射系数与阻抗的关系
在地震波垂直入射的情况下，有两种不同的公式可将阻抗I和反射系数r联系在一起，一种是：
其中，Ij指的是第j地层的阻抗值，也即向量I的第j个元素，rj指的是第j地层分界面的反射系数值，也即向量r的第j个元素。从方程6，我们可以看出反射系数r和阻抗I的关系是非线性的。
另一种将阻抗I和反射系数r联系在一起的公式为：
从上式7可以看出，反射系数r和阻抗I的关系仍然是非线性的，但是反射系数r和阻抗log(I)的关系是线性的。由方程6式，在反演得到反射系数序列r后，我们可以通过下面的递推公式8来获得拟得到的声波阻抗序列，
同样，根据方程7式，在反演得到反射系数序列r后，我们也可以通过下面的递推公式9来获得拟得到的声波阻抗序列，
当第一个地层的阻抗值I1已知时，公式8和9可以用于递推得到所有地层所对应的阻抗序列I1，I2，…，IN。当反射系数的绝对值|rj|≤0.3时，6式和7式近乎等价一致，8式和9式近乎等价一致。在小反射系数近似的情况下，我们采用公式7将阻抗I映射到反射系数r，也即
r＝Dx (10)
其中，这样我们可以通过线性反演理论先恢复x＝log(I)，再由log(I)恢复得到阻抗I。
稳定Q补偿的声波阻抗反演
将方程10代入到公式5，我们得到：
s＝WF-1FQDx (11)
不同于NSRI，本发明这里旨在求解得到的x不再是稀疏的，也即每一个采样点上都会存在数值。所以我们这里不需要再对解施加稀疏约束。另外，同NSRI一样的是，本发明方法巧妙地避免了传统反Q补偿存在的数值不稳定问题，这是通过将Q吸收衰减因子集成到正演模拟算子A＝WF-1FQD，并采用迭代求解的方法得到x。不同于传统方法，本发明无需设计成指数增长的Q补偿算子。所以本发明Q补偿技术是无条件稳定的。当Q→+∞时，本发明即为传统的阻抗反演技术。
本发明采用较为成熟稳定的LSQR算法或下述第一算法求解公式11，求解时正演算子A＝WF-1FQD，观测数据向量b＝s，然后得到x＝log(I)，再由log(I)做指数变换恢复得到阻抗I，也即I＝ex＝elog(I)。在反演得到x后，我们通过计算
sDQ＝WDx (12)
得到反Q补偿后的地震记录数据sDQ。如果LSQR算法或第一算法的初始值全部设为零值，本发明届时反演得到的是相对声波阻抗值。如果LSQR算法或第一算法的初始值设为先验背景场的数据，本发明届时反演得到的是绝对声波阻抗值。
参考图4，第一算法具体的包含下述步骤：
h)获取下述参数：正演算子A＝WF-1FQD中各参数，观测数据向量b＝s，x初始值向量x0，最大迭代次数lmax。
i)初始化:迭代变量l＝0,将x0赋值给y0。
j)重复进行下述各d至f步骤lmax次。
k)计算上标H代表共轭转置。
l)将xl+1更新为yl+1。
m)将l更新为l+1。
n)得出解x＝xlmax。
Marmousi2模型数值例子Marmousi模型是在地震勘探数据处理方法研究中应用较多的例子。该模型的处理结果主要在图5到图7中进行展示。在此模型实验过程中，我们采用了主频为40Hz的Ricker子波作为震源子波，时间采样间隔4毫秒。图5(a)为真实的阻抗模型，在模型顶部水层的Q值设定为无穷大，其余地层的Q值设定为30。图5(b)是在初始模型全为0值的情况下，本发明反演得到的相对阻抗结果，可以看出本发明反演结果分辨率较高，同真实模型构造基本一致；在给定图5(c)一维初始模型的情况下，本发明反演得到了图5(d)所示的结果，可以看出本发明反演结果分辨率和精度均较高，同真实模型近乎一致。我们可以看出本发明反演结果对此复杂的地质模型反演出了其主要构造，分辨率较高。针对图5(c)误差较大的初始模型，本发明反演结果是较为精确、可以接受的。好的初始阻抗模型将会产生更好的反演结果。
图6(a)是真实不含Q吸收衰减效应的模拟数据，图6(b)是本发明技术进行稳定Q值补偿后得到的结果。图6(c)是将图6(a)和图6(b)中的结果以波形曲线的方式画在一起，便于对比，灰色粗实线代表的是合成的真实不含Q吸收衰减效应的地震数据中的一道，黑色粗实线代表的是合成的含Q吸收衰减效应的地震数据中的一道，黑色粗虚线代表的是本发明稳定Q吸收衰减补偿的地震数据中的一道。图7给出了Q值补偿前后数据Fourier振幅谱的对比，其中，黑色实线代表的是合成的真实不含Q吸收衰减效应的地震数据的频率域振幅谱，灰色实线代表的是合成的含Q吸收衰减效应的地震数据的频率域振幅谱，黑色虚线代表的是本发明稳定Q吸收衰减补偿的地震数据的频率域振幅谱。从图6(a)-(c)和图7可以看出，本发明实现了稳定的Q值补偿，避免了传统反Q补偿方法存在的固有数值不稳定问题；高频成分得以较好地补偿；补偿结果在时间域波形上和频率域频谱上均和真实不含Q吸收衰减的结果较为一致，表明了本发明方法的可行性、有效性和优越性。
Colony实际地震数据例子
紧接着我们采用图8(a)所示的实际数据来测试本发明方法，其中叠画了三条商业软件拾取的层位，在图8(a)的右上角叠画了本发明反演时所采用的子波，其为商业软件子波提取模块的输出。图8(b)是未进行Q值补偿得到的相对阻抗结果，也即采用一个无穷大Q模型时反演的相对阻抗结果，其中叠画了一条实际测井曲线，并标注了气藏的位置。图8(c)是本发明考虑Q值补偿，反演得到的相对阻抗结果，也即采用一个尝试性Q模型时反演的相对阻抗结果。图8(d)是本发明Q值补偿后得到的地震数据。在图8(b)和图8(c)中叠画了一条来自实际钻井得到的声波阻抗曲线。图8(a)-(d)中叠画了3条由商业软件拾取的层位。从图8(b)和图8(c)可以看出，本发明反演结果同实际钻探井曲线较为一致，精确地定位出了气藏所在的位置。图9给出了本发明Q值补偿前后数据Fourier振幅谱，其中点线代表补偿前数据的振幅谱，实线代表本发明Q吸收衰减补偿后数据的振幅谱。从图8(a)-(d)和图9可以看出本发明方法，一方面反演出了高分辨率阻抗结果可用于油气矿藏检测，另一方面对地震数据实现了稳定的Q值补偿。
Erskine三维实际地震数据例子
最后采用图10(a)和图11(a)所示的另一实际地震数据来测试本发明方法的有效性，其中10(a)叠画了五条商业软件拾取的层位，在图10(a)的右上角叠画了本发明反演时所采用的子波，其为商业软件子波提取模块的输出。图10(b)给出了某商业软件处理得到的阻抗结果。图10(c)是本发明在未进行Q值补偿的情况下得到的相对阻抗结果，也即并未进行Q吸收衰减补偿情况下得到的相对阻抗结果，其中叠画了一条实际测井曲线。图10(d)是本发明在进行Q值补偿的同时反演得到的相对阻抗结果，也即采用一个尝试性Q模型时反演的相对阻抗结果。图10(e)是本发明迭代反演过程最后产生的数据误差。图10(f)是本发明进行稳定Q值补偿后得到的地震数据。图11(a)-(d)给出了本发明对此3维数据的处理结果，图11(a)为整个Erskine三维实际地震数据图；图11(b)为本发明对于整个Erskine三维实际地震数据采用一个无穷大Q模型时反演的相对阻抗结果(也即并未进行Q吸收衰减补偿)图；图11(c)为本发明对于整个Erskine三维实际地震数据采用一个尝试性Q模型时反演的相对阻抗结果(也即进行了Q吸收衰减补偿)图；图11(d)为本发明对于整个Erskine三维实际地震数据进行Q吸收衰减补偿后的结果图。在图10(a)-(f)中，叠画了5条由商业软件拾取的层位。在图10(b)至图10(d)中，叠画了一条来自实际钻探井位的曲线。从图10(c)和图10(d)可以看出本发明反演结果同实际钻探井位曲线较为一致。从图10(a)-(f)和图11(a)-(d)同样可以看出本发明方法，一方面反演出了高分辨率高精度阻抗结果可用于油气矿藏检测，另一方面对地震数据实现了稳定的Q值补偿。本发明方法处理结果地震同相轴的横向连续性较好，不存在“挂面条”等现象。证实了本发明方法的有效性、优越性。
上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。
Claims (6)
1.一种Q吸收衰减补偿的声波阻抗反演方法，其特征在于，包含如下步骤：
(1)获取数据s、W、F-1、FQ以及D；
(2)采用LSQR算法或下述第一算法求解模型s＝WF-1FQDx得到x，求解时正演算子A＝WF- 1FQD，观测数据向量为b＝s；
(3)根据sDQ＝WDx得到反Q补偿后的地震记录数据sDQ作为反演结果；
第一算法：
a)获取下述参数：正演算子A＝WF-1FQD，观测数据向量s，x的初始值向量x0，最大迭代次数lmax；
b)初始化：迭代变量l＝0,将x0赋值给y0；
c)重复进行下述各d至f步骤lmax次；
d)计算上标H代表共轭转置；
e)将xl+1更新为yl+1；
f)将l更新为l+1；
g)得出解
其中，W是地震子波所构建的褶积矩阵，代表了反Fourier变换矩阵，F-1∈CN×N，N代表了时间离散采样点数，是品质因子Q所构建的吸收衰减矩阵，FQ∈CN×N，N×N表示行数和列数均为N的矩阵，D中其他元素为0，s是实际获取的地震记录s(t)的向量表达式， 为输入地震子波w(t)的Fourier变换频谱，ω＝2πf为角频率，r(z)为在深度z处的反射系数，α(ω,τz,Qz)为Q吸收衰减因子，e为指数函数，t代表时间，τz为传播旅行时，其定义为x是与深度有关的积分变量，vr是与参考频率ωr有关的相速度，Qeff(z)是一种等效Q值，其定义为：
2.根据权利要求1所述的Q吸收衰减补偿的声波阻抗反演方法，其特征在于，还包括根据求解阻抗矩阵的步骤。
3.根据权利要求1所述的Q吸收衰减补偿的声波阻抗反演方法，其特征在于，步骤(2)中在采用LSQR算法或下述第一算法进行求解时，x的初始值向量x0全部设为零值或者为先验背景场对应的值。
4.一种Q吸收衰减补偿的声波阻抗反演系统，其特征在于，包含如下模块：
数据获取模块，用于获取数据s、W、F-1、FQ以及D；
算法求解模块，用于采用LSQR算法或下述第一算法求解模型s＝WF-1FQDx得到x，求解时正演算子A＝WF-1FQD，观测数据向量为b＝s；
地震记录计算模块，用于根据sDQ＝WDx得到反Q补偿后的地震记录数据sDQ作为反演结果；
第一算法：
a)获取下述参数：正演算子A＝WF-1FQD，观测数据向量s，x的初始值向量x0，最大迭代次数lmax；
b)初始化：迭代变量l＝0,将x0赋值给y0；
c)重复进行下述各d至f步骤lmax次；
d)计算上标H代表共轭转置；
e)将xl+1更新为yl+1；
f)将l更新为l+1；
g)得出解
其中，W是地震子波所构建的褶积矩阵，代表了反Fourier变换矩阵，F-1∈CN×N，N代表了时间离散采样点数，是品质因子Q所构建的吸收衰减矩阵，FQ∈CN×N，N×N表示行数和列数均为N的矩阵，D中其他元素为0，s是实际获取的地震记录s(t)的向量表达式， 为输入地震子波w(t)的Fourier变换频谱，ω＝2πf为角频率，r(z)为在深度z处的反射系数，α(ω,τz,Qz)为Q吸收衰减因子，e为指数函数，t代表时间，τz为传播旅行时，其定义为x是与深度有关的积分变量，vr是与参考频率ωr有关的相速度，Qeff(z)是一种等效Q值，其定义为：
5.根据权利要求4所述的Q吸收衰减补偿的声波阻抗反演系统，其特征在于，还包括阻抗求解模块，用于根据求解阻抗矩阵
6.根据权利要求4所述的Q吸收衰减补偿的声波阻抗反演系统，其特征在于，算法求解模块中在采用LSQR算法或下述第一算法进行求解时，x的初始值向量x0全部设为零值或者为先验背景场对应的值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710994606.5A CN107894612B (zh) | 2017-10-23 | 2017-10-23 | 一种q吸收衰减补偿的声波阻抗反演方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710994606.5A CN107894612B (zh) | 2017-10-23 | 2017-10-23 | 一种q吸收衰减补偿的声波阻抗反演方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107894612A true CN107894612A (zh) | 2018-04-10 |
| CN107894612B CN107894612B (zh) | 2019-05-31 |
Family
ID=61802842
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710994606.5A Expired - Fee Related CN107894612B (zh) | 2017-10-23 | 2017-10-23 | 一种q吸收衰减补偿的声波阻抗反演方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107894612B (zh) |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109143340A (zh) * | 2018-08-20 | 2019-01-04 | 中国海洋石油集团有限公司 | 一种基于常q模型的粘弹介质地震波模拟方法及系统 |
| CN110873900A (zh) * | 2018-09-04 | 2020-03-10 | 中国石油化工股份有限公司 | 频率域叠前地震道q补偿方法及系统 |
| CN110988991A (zh) * | 2019-12-16 | 2020-04-10 | 中国石油大学(北京) | 一种弹性参数反演方法、装置及系统 |
| CN111708081A (zh) * | 2020-05-29 | 2020-09-25 | 成都理工大学 | 考虑衰减频散的深度域地震记录合成方法 |
| CN112147687A (zh) * | 2019-06-28 | 2020-12-29 | 中国石油化工股份有限公司 | 一种储层含气性预测方法及预测系统 |
| CN112485826A (zh) * | 2020-11-12 | 2021-03-12 | 中国地质大学（武汉） | 绝对波阻抗反演成像方法、装置、设备及存储介质 |
| CN112596107A (zh) * | 2020-12-11 | 2021-04-02 | 中国科学院地质与地球物理研究所 | 弹性参数反演方法及装置 |
| CN113820747A (zh) * | 2020-06-18 | 2021-12-21 | 中国石油化工股份有限公司 | Q值计算与反射系数反演方法、装置、电子设备及介质 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103293551A (zh) * | 2013-05-24 | 2013-09-11 | 中国石油天然气集团公司 | 一种基于模型约束的阻抗反演方法及系统 |
| CN104280767A (zh) * | 2013-07-12 | 2015-01-14 | 中国石油天然气集团公司 | 一种基于柯西分布的稀疏脉冲反演方法 |
| CN105589833A (zh) * | 2014-10-23 | 2016-05-18 | 陕西中浩源水电工程有限公司 | 基于lsqr法频率域波形反演的存储方法 |
| CN106291682A (zh) * | 2015-06-01 | 2017-01-04 | 中国石油化工股份有限公司 | 一种基于基追踪方法的叠后声波阻抗反演方法 |
| CN106291677A (zh) * | 2015-05-22 | 2017-01-04 | 中国石油化工股份有限公司 | 一种基于匹配追踪方法的叠后声波阻抗反演方法 |
-
2017
- 2017-10-23 CN CN201710994606.5A patent/CN107894612B/zh not_active Expired - Fee Related
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103293551A (zh) * | 2013-05-24 | 2013-09-11 | 中国石油天然气集团公司 | 一种基于模型约束的阻抗反演方法及系统 |
| CN104280767A (zh) * | 2013-07-12 | 2015-01-14 | 中国石油天然气集团公司 | 一种基于柯西分布的稀疏脉冲反演方法 |
| CN105589833A (zh) * | 2014-10-23 | 2016-05-18 | 陕西中浩源水电工程有限公司 | 基于lsqr法频率域波形反演的存储方法 |
| CN106291677A (zh) * | 2015-05-22 | 2017-01-04 | 中国石油化工股份有限公司 | 一种基于匹配追踪方法的叠后声波阻抗反演方法 |
| CN106291682A (zh) * | 2015-06-01 | 2017-01-04 | 中国石油化工股份有限公司 | 一种基于基追踪方法的叠后声波阻抗反演方法 |
Non-Patent Citations (1)
| Title |
|---|
| 柴新涛 等: "基于LSQR算法的谱反演方法研究", 《石油物探》 * |
Cited By (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109143340A (zh) * | 2018-08-20 | 2019-01-04 | 中国海洋石油集团有限公司 | 一种基于常q模型的粘弹介质地震波模拟方法及系统 |
| CN109143340B (zh) * | 2018-08-20 | 2020-03-10 | 中国海洋石油集团有限公司 | 一种基于常q模型的粘弹介质地震波模拟方法及系统 |
| CN110873900A (zh) * | 2018-09-04 | 2020-03-10 | 中国石油化工股份有限公司 | 频率域叠前地震道q补偿方法及系统 |
| CN110873900B (zh) * | 2018-09-04 | 2021-07-27 | 中国石油化工股份有限公司 | 频率域叠前地震道q补偿方法及系统 |
| CN112147687A (zh) * | 2019-06-28 | 2020-12-29 | 中国石油化工股份有限公司 | 一种储层含气性预测方法及预测系统 |
| CN110988991A (zh) * | 2019-12-16 | 2020-04-10 | 中国石油大学(北京) | 一种弹性参数反演方法、装置及系统 |
| CN110988991B (zh) * | 2019-12-16 | 2021-06-25 | 中国石油大学(北京) | 一种弹性参数反演方法、装置及系统 |
| CN111708081A (zh) * | 2020-05-29 | 2020-09-25 | 成都理工大学 | 考虑衰减频散的深度域地震记录合成方法 |
| CN111708081B (zh) * | 2020-05-29 | 2022-04-15 | 成都理工大学 | 考虑衰减频散的深度域地震记录合成方法 |
| CN113820747A (zh) * | 2020-06-18 | 2021-12-21 | 中国石油化工股份有限公司 | Q值计算与反射系数反演方法、装置、电子设备及介质 |
| CN112485826A (zh) * | 2020-11-12 | 2021-03-12 | 中国地质大学（武汉） | 绝对波阻抗反演成像方法、装置、设备及存储介质 |
| CN112596107A (zh) * | 2020-12-11 | 2021-04-02 | 中国科学院地质与地球物理研究所 | 弹性参数反演方法及装置 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107894612B (zh) | 2019-05-31 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107894612B (zh) | 一种q吸收衰减补偿的声波阻抗反演方法及系统 | |
| Virieux et al. | An overview of full-waveform inversion in exploration geophysics | |
| Borisov et al. | Application of 2D full-waveform inversion on exploration land data | |
| CN101334483B (zh) | 一种在地震数据处理中衰减瑞雷波散射噪声的方法 | |
| Pan et al. | Love-wave waveform inversion in time domain for shallow shear-wave velocity | |
| CN103792573B (zh) | 一种基于频谱融合的地震波阻抗反演方法 | |
| CN107015274A (zh) | 一种缺失地震勘探数据恢复重构方法 | |
| Ha et al. | Laplace-domain full-waveform inversion of seismic data lacking low-frequency information | |
| Wang et al. | High-resolution wave-equation amplitude-variation-with-ray-parameter (AVP) imaging with sparseness constraints | |
| CN108037531A (zh) | 一种基于广义全变分正则化的地震反演方法及系统 | |
| CN104122581B (zh) | 一种叠后声波阻抗反演方法 | |
| CN104516018A (zh) | 一种地球物理勘探中岩性约束下的孔隙度反演方法 | |
| CN106291682B (zh) | 一种基于基追踪方法的叠后声波阻抗反演方法 | |
| CN111025387B (zh) | 一种页岩储层的叠前地震多参数反演方法 | |
| CN105089652A (zh) | 一种拟声波曲线重构与稀疏脉冲联合反演方法 | |
| CN107356967A (zh) | 一种压制地震资料强屏蔽干扰的稀疏优化方法 | |
| CN104237945A (zh) | 一种地震资料自适应高分辨处理方法 | |
| CN107132579A (zh) | 一种保地层结构的地震波衰减补偿方法 | |
| CN107817526A (zh) | 叠前地震道集分段式振幅能量补偿方法及系统 | |
| CN101201409A (zh) | 一种地震数据变相位校正方法 | |
| CN113031068B (zh) | 一种基于反射系数精确式的基追踪叠前地震反演方法 | |
| CN109212599A (zh) | 一种地震数据的叠前同步反演方法 | |
| CN113945982A (zh) | 用于去除低频和低波数噪声以生成增强图像的方法和系统 | |
| CN107621654A (zh) | 一种基于最大相关熵的地震叠后波阻抗反演方法 | |
| CN107179547A (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 | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190531 Termination date: 20201023 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |