CN114063151A - 高精度叠前地震数据衰减属性提取方法及系统 - Google Patents
高精度叠前地震数据衰减属性提取方法及系统 Download PDFInfo
- Publication number
- CN114063151A CN114063151A CN202111355584.0A CN202111355584A CN114063151A CN 114063151 A CN114063151 A CN 114063151A CN 202111355584 A CN202111355584 A CN 202111355584A CN 114063151 A CN114063151 A CN 114063151A
- Authority
- CN
- China
- Prior art keywords
- matrix
- seismic
- frequency
- formula
- wavelet
- 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.)
- Pending
Links
- 238000000605 extraction Methods 0.000 title claims abstract description 18
- 239000011159 matrix material Substances 0.000 claims description 57
- 238000000034 method Methods 0.000 claims description 28
- 230000003595 spectral effect Effects 0.000 claims description 27
- 238000001228 spectrum Methods 0.000 claims description 21
- 230000015572 biosynthetic process Effects 0.000 claims description 19
- 238000009499 grossing Methods 0.000 claims description 11
- 238000000354 decomposition reaction Methods 0.000 claims description 10
- 239000000126 substance Substances 0.000 claims description 9
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims 2
- 230000017105 transposition Effects 0.000 claims 2
- 238000001914 filtration Methods 0.000 abstract description 4
- 239000010410 layer Substances 0.000 description 12
- 238000012545 processing Methods 0.000 description 6
- 238000007493 shaping process Methods 0.000 description 6
- 238000013459 approach Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 230000009466 transformation Effects 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 238000009792 diffusion process Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 239000011229 interlayer Substances 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000013075 data extraction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000012417 linear regression Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000010363 phase shift Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000013519 translation Methods 0.000 description 1
Images
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
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/51—Migration
- G01V2210/512—Pre-stack
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值反演方程求取CMP道集层位处等效Q值，最后根据等效Q与层间Q的关系求得层间Q值。本发明能避免叠后Q值的不稳定性，提高Q值估算准确性，通过反Q滤波能有效恢复实际地震信号能量，获得高分辨率和高信噪比的地震数据。
Description
技术领域
本发明涉及油气地球物理勘探技术领域，具体的说，涉及一种高精度叠前地震数据衰减属性提取方法及系统。
背景技术
与叠后地震资料不同，叠前地震数据包含更多地层信息，故从叠前地震资料提取Q值具有更高的精度。Reine等(2009)将时-频分析方法引入到Q值的提取中，比较了利用STFT、Gabor变换、S变换和小波变换所得振幅谱进行Q值反演时对薄层调谐的适应能力，通过比较发现变时窗的S变换和小波变换具有较强的鲁棒性；Wang(2016)分析了S变换的缺陷后提出了一种新窗口函数的S变换，并将其应用于Q值估算及地震数据补偿；Wang(2004，2014)利用时频函数将时频谱转变为一维数据，通过拟合实际数据曲线与理论曲线，最后进行补偿分析得到Q值，该方法较为稳定，但计算步骤较繁琐。对于叠前数据，叠前Q值估算难度相对大于叠后数据，受噪声影响大。Dasgupta和Clark(1998)先对CMP道集进行动校正，然后利用Q值和偏移距的关系拟合求出Q值，该方法首次将谱比法用于CMP道集，并取得了较好效果。Reine等(2012a，2012b)提出利用频率差与走时差坐标联合反演减小残留的谱干扰，同时将t-x域的数据转换到τ-p域，该方法提高了Q值求取的精度。但由于时频分辨率及其他因素约束，导致提取的Q值准确度较低。
Gladwin(1974)提出脉冲上升时间法，该方法利用地震波在介质中传播时振幅的变化，通过子波的最大振幅、振幅的最大斜率与Q值的关系式来求得Q值，相应的方法还有脉冲宽度法(Wright和Hoy，1981)、振幅衰减法(Tyce，1981)等。该类方法统称为时间域的Q值估计方法。由于波在传播过程中其振幅受到散射、几何扩散等因素的影响，该类估计方法的精度往往不高。除了时间域的Q值估计方法外，还有频率域的Q值估计方法，比如，Bath(1974)提出将时间域数据转换到频率域，对不同时刻数据取对数比值，再利用振幅谱对数比值与频率的关系得到Q值；Quan和Harris(1997)提出先求出震源信号的质心频率，再求出接收信号的质心频率，利用两者之差与Q值的关系求出Q值。除此之外，该类方法还有峰值频率偏移法(Zhang and Ulrych，2002)等，Tonn(1991)对比了多种Q值的估计方法，结果表明在噪声较大的情况下，每种估计方法都没有较好的效果，含有少量噪声时，谱比法的效果优于其他方法。相对于时间域Q值估算方法，频率域方法通常需要利用一个时间窗去截取地震记录，然后通过傅里叶变换得到频率域的数据，所以处理中对窗口类型及长度的要求较高，如果选取的窗口不适宜，往往会影响Q值求取的精度。
因此，现有技术中提取Q值多数从叠后地震记录进行提取，对于叠后地震记录由于叠加次数越高，低通滤波性质越明显，导致提取的Q值分辨率不能满足生产实际的需求。
发明内容
本发明针对叠后地震数据提取Q值准确度较低，提出了一种高精度叠前地震数据衰减属性提取方法及系统，该方法利用反演谱分解及整形正则化技术，使得Q值提取更准确，从而有利于辨别含油气区域。
本发明的具体技术方案如下：
根据本发明的第一技术方案，提供了一种高精度叠前地震数据衰减属性提取方法，包括以下步骤：
基于不同地震子波的主频、不同主频的复雷克子波以及对应时刻的噪声，通过如下公式(2)确定褶积模型：
其中，s(t)为地震数据，M表示地震子波的个数，fi表示i个地震子波的主频；w(t,fi)表示不同主频的复雷克子波；r(t,fi)为不同主频子波对应的复反射系数，所述反射系数至少包含子波频率和相位信息，n(t)为噪声，t表示时刻；
将所述公式(2)表示成矩阵形式如下公式(3)所示：
其中，
表示子波褶积矩阵，每个子波的主频为fi；
表示每个子波对应的复反射系数矩阵；
表示复子波库的褶积矩阵；
表示所有复反射系数组成的矩阵；
设定地震子波传播时间t后的振幅为A(t,f)，并通过如下公式(4)表示A(t，f)：
A(t,f)＝P(t)A(t0,f)e-πft/Q 公式(4)
其中，t为旅行时间，f为频率；A(t,f)为地震子波传播时间t后的振幅；A(t0,f)为地震波在t0时刻的振幅；P(t)表示与频率无关的能量损失；Q为地层的品质因子；
基于公式(4)，获得地震子波经过时间t1到达某层顶的地震波振幅以及地震子波经过时间t2到达某层底的地震波振幅，分别表示为：
其中，A(t1,f)为地震子波经过时间t1到达某层顶的地震波振幅，f为频率，P(t1)表示经过时间t1与频率无关的能量损失，A(t0,f)为地震波在t0时刻的振幅，Q为地层的品质因子；
其中，A(t2,f)为地震子波经过时间t1到达某层顶的地震波振幅，f为频率，P(t2)表示经过时间t2与频率无关的能量损失，A(t0,f)为地震波在t0时刻的振幅，Q为地层的品质因子；
将公式(5)与公式(6)取比值并在式两端取对数，并通过如下公式(7)表示：
其中，ln[r(f)]表示对数谱比值，ln(D)表示与频率无关的能量损失，f为频率，Q为地层的品质因子；
通过求取对数谱比值与频率的斜率得到Q值。
根据本发明的第二技术方案，提供了一种高精度叠前地震数据衰减属性提取系统，包括处理器，所述处理器被配置为：
基于不同地震子波的主频、不同主频的复雷克子波以及对应时刻的噪声，通过如下公式(2)确定褶积模型：
其中，s(t)为地震数据，*表示褶积运算，M表示地震子波的个数，fi表示i个地震子波的主频；w(t,fi)表示不同主频的复雷克子波；r(t,fi)为不同主频子波对应的复反射系数，所述反射系数至少包含子波频率和相位信息，n(t)为噪声，t表示时刻；
将所述公式(2)表示成矩阵形式如下公式(3)所示：
其中，
表示子波褶积矩阵，每个子波的主频为fi；
表示每个子波对应的复反射系数矩阵；
表示复子波库的褶积矩阵；
表示所有复反射系数组成的矩阵；
设定地震子波传播时间t后的振幅为A(t,f)，并通过如下公式(4)表示A(t，f)：
A(t，f)＝P(t)A(t0，f)e-πft/Q 公式(4)
其中，t为旅行时间，f为频率；A(t，f)为地震子波传播时间t后的振幅；A(t0,f)为地震波在t0时刻的振幅；P(t)表示与频率无关的能量损失；Q为地层的品质因子；
基于公式(4)，获得地震子波经过时间t1到达某层顶的地震波振幅以及地震子波经过时间t2到达某层底的地震波振幅，分别表示为：
其中，A(t1,f)为地震子波经过时间t1到达某层顶的地震波振幅，f为频率，P(t1)表示经过时间t1与频率无关的能量损失，A(t0,f)为地震波在t0时刻的振幅，Q为地层的品质因子；
其中，A(t2,f)为地震子波经过时间t1到达某层顶的地震波振幅，f为频率，P(t2)表示经过时间t2与频率无关的能量损失，A(t0,f)为地震波在t0时刻的振幅，Q为地层的品质因子；
将公式(5)与公式(6)取比值并在式两端取对数，并通过如下公式(7)表示：
其中，ln[r(f)]表示对数谱比值，ln(D)表示与频率无关的能量损失，f为频率，Q为地层的品质因子；
通过求取对数谱比值与频率的斜率得到Q值。
根据本发明实施例的高精度叠前地震数据衰减属性提取方法及系统，能避免叠后Q值的不稳定性，提高Q值估算准确性，通过反Q滤波能有效恢复实际地震信号能量，获得高分辨率和高信噪比的地震数据。
附图说明
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍。在所有附图中，类似的元件或部分一般由类似的附图标记标识。附图中，各元件或部分并不一定按照实际的比例绘制。
图1为根据本发明实施例的一种高精度叠前地震数据衰减属性提取方法的流程图；
图2为根据本发明实施例的高精度叠前地震数据衰减属性提取系统的硬件图。
具体实施方式
下面将对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定发明。
现在结合说明书附图对本发明做进一步的说明。
图1示出了根据本发明实施例的一种高精度叠前地震数据衰减属性提取方法的流程图。如图1所示，本发明实施例提供了一种高精度叠前地震数据衰减属性提取方法。
具体说来，本发明实施例基于Robinson(1967)提出的经典的褶积模型，用恒定的地震子波与反射系数褶积得到地震记录。公式为
s(t)＝w(t)*r(t)+n(t) 公式(1)
其中，*表示褶积运算；s(t)为地震数据，t表示地震波传播的旅行时间；w(t)为地震子波；r(t)为反射系数；n(t)为噪声。
发明人通过实验发现，公式(1)表示的是平稳地震信号的褶积模型，但地震波向地下传播时会发生吸收衰减，故不同时刻子波的主频和相位应该有所差异，因此该过程可表示为
其中，M表示子波的个数，fi表示不同的主频；w(t，fi)表示不同主频的复雷克子波；r(t,fi)为不同主频子波对应的复反射系数，该反射系数包含子波频率和相位信息。
将公式(2)表示成矩阵形式为
其中，
表示子波褶积矩阵，每个子波的主频为fi；
表示每个子波对应的复反射系数矩阵；
表示复子波库的褶积矩阵；
表示所有复反射系数组成的矩阵。
假设地震子波传播时间t后的振幅为A(t,f)，根据Futterman理论(Futterman,1962)与频率域相移原理(Gazdag and Sguazzero,1984)，A(t,f)可表示为：
A(t,f)＝P(t)A(t0,f)e-πft/Q 公式(4)
其中，t为旅行时间，f为频率；A(t,f)为ISD方法求取的振幅；A(t0,f)为地震波在t0时刻的振幅；P(t)表示由于几何扩散、透射反射等引起的与频率无关的能量损失；Q为地层的品质因子。
根据公式(4)可知，地震子波经过时间t1和t2到达某层顶、底的地震波振幅可表示为：
公式(5)与公式(6)取比值可得：
公式(7)两端取对数表示为：
其中，ln(D)表示与频率无关的能量损失。通过求取对数谱比值与频率的斜率即可求得Q值。
在上述求取对数谱比时，由于公式(6)中分母趋近于无穷大时，r(f)的值趋近于零，所以公式(7)通常求取的对数比值会存在不稳定性。为了使求出的Q值相对准确，减少不稳定的一个简单的方法是在分母中添加一个小值δ，则公式(6)改写为：
发明人通过实验发现由公式(8)得到的结果，容易产生振动，且不够光滑。因此，本发明实施例引入反演理论求解该方程，公式(8)可表示为：
d＝Wu 公式(9)
其中，u为待求解的谱比矩阵，d和W为反演谱分解得到的A(t2,f)和A(t1,f)矩阵。
为了使公式(9)的得到的解更稳定，下面本发明实施例将使用整形正则化来进行求解。
基于现有技术中的考虑整形算子作用的整形正则化理论。该理论可以方便地选择正则化算子，如高斯平滑算子和三角平滑算子，对于公式(9)所表示的典型线性正问题，根据Tikhonov正则化，正则化解可以表示为：
D是正则化算子和μ是一个标量扩展参数。整形正则化加入光滑算子S。一般意义上，S可以看作约束模型在可接受空间中的映射。公式(11)作为正则化算子D的定义可以写成：
S＝[I+S(WTW-I)]-1SWTor μ2DTD＝S-1-I 公式(11)
将公式(10)代入公式(9)，得到了整形正则化问题的形式解为：
本发明实施例考虑到算子W可能需要单位的转化。因此，在本发明实施例中，可以引入算子W的1/η进行缩放比例，则将公式(12)重写为：
其中，整形算子S可以是三角平滑算子，也可以是高斯平滑算子。
对于CMP道集，现有技术中普遍采用是QVO(Dasgupta and Clark，1998)法，QVO通常采用两步拟合求出Q值，其中第一步：求出对数谱比得：
其中，Δt1为零炮检距地震波通过目的层上下界面时所用的时间差；Qx为不同炮检距处的Q值。拟合出自然对数谱比值与频率f之间的斜率。
第二步：跑检距归零处理得：
其中，KD为零偏移距为x时的谱比斜率；v为波速；Q0为零偏移距的Q值；Δt0为零偏移距处地震波通过参考同相轴到目标同相轴所用的双程旅行时差。KD经最小平方线性回归即可得到参考同相轴到目标反射轴的射线平均Q值。
本发明实施例基于以上理论，将两步QVO方法改造为单步QVO方法。
具体来说，本发明实施例中，对于
其中，rNM表示第N道、频率为fM的自然对数谱比值。
利用R表示自然对数谱比值组成的向量，C为系数矩阵，U为待求参数向量，则公式(16)可表示为：
R＝CU 公式(17)
由最小二乘原理可得：
U＝(C′C)-1(C′R) 公式(18)
通过公式(18)求得U后即可得到Q值。
本发明实施例提供的方法先利用反演复谱分解求取地震波时频谱，通过整形正则化得到对数谱比值，再利用单步Q值反演方程求取CMP道集层位处等效Q值，最后根据等效Q与层间Q的关系求得层间Q值。能避免叠后Q值的不稳定性，提高Q值估算准确性，通过反Q滤波能有效恢复实际地震信号能量，获得高分辨率和高信噪比的地震数据。
图2示出了根据本发明实施例的高精度叠前地震数据衰减属性提取系统的硬件图，如图2所示，本发明实施例还提供了一种高精度叠前地震数据衰减属性提取系统，该系统200包括处理器201，所述处理器201配置为实现根据本发明各个实施例所述的方法。
需要注意的是处理器201可以是包括一个以上通用处理设备的处理设备，诸如微处理器、中央处理单元(CPU)、图形处理单元(GPU)等。更具体地，处理器201可以是复杂指令集计算(CISC)微处理器、精简指令集计算(RISC)微处理器、超长指令字(VLIW)微处理器、运行其他指令集的处理器或运行指令集的组合的处理器。处理器201还可以是一个以上专用处理设备，诸如专用集成电路(ASIC)、现场可编程门阵列(FPGA)、数字信号处理器(DSP)、片上系统(SoC)等。
以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围，其均应涵盖在本发明的权利要求和说明书的范围当中。
Claims (6)
1.高精度叠前地震数据衰减属性提取方法，其特征在于，包括以下步骤：
基于不同地震子波的主频、不同主频的复雷克子波以及对应时刻的噪声，通过如下公式(2)确定褶积模型：
其中，s(t)为地震数据，*表示褶积运算，M表示地震子波的个数，fi表示i个地震子波的主频；w(t，fi)表示不同主频的复雷克子波；r(t，fi)为不同主频子波对应的复反射系数，所述反射系数至少包含子波频率和相位信息，n(t)为噪声，t表示时刻；
将所述公式(2)表示成矩阵形式如下公式(3)所示：
其中，
表示子波褶积矩阵，每个子波的主频为fi；
表示每个子波对应的复反射系数矩阵；
表示复子波库的褶积矩阵；
表示所有复反射系数组成的矩阵；
设定地震子波传播时间t后的振幅为A(t，f)，并通过如下公式(4)表示A(t，f)：
A(t，f)＝P(t)A(t0，f)e-πft/Q 公式(4)
其中，t为旅行时间，f为频率；A(t，f)为地震子波传播时间t后的振幅；A(t0，f)为地震波在t0时刻的振幅；P(t)表示与频率无关的能量损失；Q为地层的品质因子；
基于公式(4)，获得地震子波经过时间t1到达某层顶的地震波振幅以及地震子波经过时间t2到达某层底的地震波振幅，分别表示为：
其中，A(t1，f)为地震子波经过时间t1到达某层顶的地震波振幅，f为频率，P(t1)表示经过时间t1与频率无关的能量损失，A(t0，f)为地震波在t0时刻的振幅，Q为地层的品质因子；
其中，A(t2，f)为地震子波经过时间t1到达某层顶的地震波振幅，f为频率，P(t2)表示经过时间t2与频率无关的能量损失，A(t0，f)为地震波在t0时刻的振幅，Q为地层的品质因子；
将公式(5)与公式(6)取比值并在式两端取对数，并通过如下公式(7)表示：
其中，ln[r(f)]表示对数谱比值，ln(D)表示与频率无关的能量损失，f为频率，Q为地层的品质因子；
通过求取对数谱比值与频率的斜率得到Q值。
2.根据权利要求1所述的方法，其特征在于，所述公式(6)替换为：
引入反演理论对所述公式(8)进行求解得到：
d＝Wu 公式(9)
其中，u为待求解的谱比矩阵，d和W为反演谱分解得到的A(t2，f)和A(t1，f)矩阵；
对所述公式(9)进行整形正则化求解，得到：
其中，
为求解后的谱比矩阵，D是正则化算子，μ是一个标量扩展参数，d和W为反演谱分解得到的A(t2，f)和A(t1，f)矩阵；
在所述公式(10)中加入光滑算子S，光滑算子S通过如下公式(11)表示：
S＝[I+S(WTW-I)]-1SWT or μ2DTD＝S-1-I (11)
其中，S是光滑算子，D是正则化算子，I是单位矩阵，μ是一个标量扩展参数，d和W为反演谱分解得到的A(t2，f)和A(t1，f)矩阵；
将公式(10)代入公式(9)，得到了整形正则化问题的形式解为：
其中，
为求解后的谱比矩阵，D是正则化算子，I是单位矩阵，d和W为反演谱分解得到的A(t2，f)和A(t1，f)矩阵；
基于求解后的谱比矩阵，来确定Q值。
3.根据权利要求2所述的方法，其特征在于，基于求解后的谱比矩阵，来确定Q值，包括：
基于所述求解后的谱比矩阵，获得如下公式(16)：
其中，rNM表示第N道、频率为fM的自然对数谱比值，ΔtN为时间间隔，ln(DN)为耗散因子的对数，Q为地层的品质因子；
利用R表示自然对数谱比值组成的向量，C为系数矩阵，U为待求参数向量，则所述公式(16)表示为：
R＝CU (17)
其中，R表示自然对数谱比值组成的向量，C为系数矩阵，U为待求参数向量；
基于最小二乘原理，对所述公式(17)求解得到如下公式(18)：
U＝(C′C)-1(C′R) (18)
其中，R表示自然对数谱比值组成的向量，C为系数矩阵，C′为C的转置，U为待求参数向量；
通过公式(18)求得U后得到Q值。
4.高精度叠前地震数据衰减属性提取系统，包括处理器，其特征在于，所述处理器被配置为：
基于不同地震子波的主频、不同主频的复雷克子波以及对应时刻的噪声，通过如下公式(2)确定褶积模型：
其中，s(t)为地震数据，*表示褶积运算，M表示地震子波的个数，fi表示i个地震子波的主频；w(t，fi)表示不同主频的复雷克子波；r(t，fi)为不同主频子波对应的复反射系数，所述反射系数至少包含子波频率和相位信息，n(t)为噪声，t表示时刻；
将所述公式(2)表示成矩阵形式如下公式(3)所示：
其中，
表示子波褶积矩阵，每个子波的主频为fi；
表示每个子波对应的复反射系数矩阵；
表示复子波库的褶积矩阵；
表示所有复反射系数组成的矩阵；
设定地震子波传播时间t后的振幅为A(t，f)，并通过如下公式(4)表示A(t，f)：
A(t，f)＝P(t)A(t0，f)e-πft/Q 公式(4)
其中，t为旅行时间，f为频率；A(t，f)为地震子波传播时间t后的振幅；A(t0，f)为地震波在t0时刻的振幅；P(t)表示与频率无关的能量损失；Q为地层的品质因子；
基于公式(4)，获得地震子波经过时间t1到达某层顶的地震波振幅以及地震子波经过时间t2到达某层底的地震波振幅，分别表示为：
其中，A(t1，f)为地震子波经过时间t1到达某层顶的地震波振幅，f为频率，P(t1)表示经过时间t1与频率无关的能量损失，A(t0，f)为地震波在t0时刻的振幅，Q为地层的品质因子；
其中，A(t2，f)为地震子波经过时间t1到达某层顶的地震波振幅，f为频率，P(t2)表示经过时间t2与频率无关的能量损失，A(t0，f)为地震波在t0时刻的振幅，Q为地层的品质因子；
将公式(5)与公式(6)取比值并在式两端取对数，并通过如下公式(7)表示：
其中，ln[r(f)]表示对数谱比值，ln(D)表示与频率无关的能量损失，f为频率，Q为地层的品质因子；
通过求取对数谱比值与频率的斜率得到Q值。
5.根据权利要求4所述的系统，其特征在于，所述处理器被进一步配置为：
引入反演理论对所述公式(8)进行求解得到：
d＝Wu 公式(9)
其中，u为待求解的谱比矩阵，d和W为反演谱分解得到的A(t2，f)和A(t1，f)矩阵；
对所述公式(9)进行整形正则化求解，得到：
其中，
为求解后的谱比矩阵，D是正则化算子，μ是一个标量扩展参数，d和W为反演谱分解得到的A(t2，f)和A(t1，f)矩阵；
在所述公式(10)中加入光滑算子S，光滑算子S通过如下公式(11)表示：
S＝[I+S(WTW-I)]-1SWT or μ2DTD＝S-1-I (11)
其中，S是光滑算子，D是正则化算子，I是单位矩阵，μ是一个标量扩展参数，d和W为反演谱分解得到的A(t2，f)和A(t1，f)矩阵；
将公式(10)代入公式(9)，得到了整形正则化问题的形式解为：
其中，
为求解后的谱比矩阵，D是正则化算子，I是单位矩阵，d和W为反演谱分解得到的A(t2，f)和A(t1，f)矩阵；
基于求解后的谱比矩阵，来确定Q值。
6.根据权利要求5所述的系统，其特征在于，所述处理器被进一步配置为：
基于所述求解后的谱比矩阵，获得如下公式(16)：
其中，rNM表示第N道、频率为fM的自然对数谱比值，ΔtN为时间间隔，ln(DN)为耗散因子的对数，Q为地层的品质因子；
利用R表示自然对数谱比值组成的向量，C为系数矩阵，U为待求参数向量，则所述公式(16)表示为：
R＝CU (17)
其中，R表示自然对数谱比值组成的向量，C为系数矩阵，U为待求参数向量；
基于最小二乘原理，对所述公式(17)求解得到如下公式(18)：
U＝(C′C)-1(C′R) (18)
其中，R表示自然对数谱比值组成的向量，C为系数矩阵，C′为C的转置，U为待求参数向量；
通过公式(18)求得U后得到Q值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111355584.0A CN114063151A (zh) | 2021-11-16 | 2021-11-16 | 高精度叠前地震数据衰减属性提取方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111355584.0A CN114063151A (zh) | 2021-11-16 | 2021-11-16 | 高精度叠前地震数据衰减属性提取方法及系统 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN114063151A true CN114063151A (zh) | 2022-02-18 |
Family
ID=80272865
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111355584.0A Pending CN114063151A (zh) | 2021-11-16 | 2021-11-16 | 高精度叠前地震数据衰减属性提取方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114063151A (zh) |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103984011A (zh) * | 2014-04-16 | 2014-08-13 | 孙赞东 | 一种动态q补偿偏移方法 |
| CN106291693A (zh) * | 2015-05-21 | 2017-01-04 | 中国石油化工股份有限公司 | 一种基于广义s变换的叠前q值反演方法及系统 |
| CN107300718A (zh) * | 2016-04-14 | 2017-10-27 | 中国石油天然气股份有限公司 | 一种品质因子三维衰减模型的建立方法 |
| CN109669212A (zh) * | 2017-10-13 | 2019-04-23 | 中国石油化工股份有限公司 | 地震数据处理方法、地层品质因子估算方法与装置 |
| CN110515127A (zh) * | 2019-09-26 | 2019-11-29 | 中国石油大学(北京) | 一种地震品质因子确定方法、装置、设备、介质 |
| CN112305586A (zh) * | 2019-07-29 | 2021-02-02 | 中国石油化工股份有限公司 | 非稳态地震资料时频分析方法、计算机存储介质及系统 |
| CN112578448A (zh) * | 2020-12-03 | 2021-03-30 | 成都理工大学 | 一种高精度的地层品质因子的提取方法 |
| CN112630837A (zh) * | 2021-01-08 | 2021-04-09 | 中海石油(中国)有限公司深圳分公司 | 一种地震品质因子估算方法、装置、设备及存储介质 |
-
2021
- 2021-11-16 CN CN202111355584.0A patent/CN114063151A/zh active Pending
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103984011A (zh) * | 2014-04-16 | 2014-08-13 | 孙赞东 | 一种动态q补偿偏移方法 |
| CN106291693A (zh) * | 2015-05-21 | 2017-01-04 | 中国石油化工股份有限公司 | 一种基于广义s变换的叠前q值反演方法及系统 |
| CN107300718A (zh) * | 2016-04-14 | 2017-10-27 | 中国石油天然气股份有限公司 | 一种品质因子三维衰减模型的建立方法 |
| CN109669212A (zh) * | 2017-10-13 | 2019-04-23 | 中国石油化工股份有限公司 | 地震数据处理方法、地层品质因子估算方法与装置 |
| CN112305586A (zh) * | 2019-07-29 | 2021-02-02 | 中国石油化工股份有限公司 | 非稳态地震资料时频分析方法、计算机存储介质及系统 |
| CN110515127A (zh) * | 2019-09-26 | 2019-11-29 | 中国石油大学(北京) | 一种地震品质因子确定方法、装置、设备、介质 |
| CN112578448A (zh) * | 2020-12-03 | 2021-03-30 | 成都理工大学 | 一种高精度的地层品质因子的提取方法 |
| CN112630837A (zh) * | 2021-01-08 | 2021-04-09 | 中海石油(中国)有限公司深圳分公司 | 一种地震品质因子估算方法、装置、设备及存储介质 |
Non-Patent Citations (4)
| Title |
|---|
| 刘国昌 等: "基于整形正则化和S变换的Q值估计方法", 《石油地球物理勘探》 * |
| 文莙翔: "基于稀疏反演谱分解的Q值估计方法研究", 《中国优秀硕士学位论文全文数据库 基础科学辑》 * |
| 王小杰 等: "基于叠前地震数据的地层Q值估计", 《石油地球物理勘探》 * |
| 郝亚炬: "高精度Q值反演方法研究", 《中国优秀硕士学位论文全文数据库 基础科学辑》 * |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Liu et al. | Time-frequency analysis of seismic data using local attributes | |
| RU2579164C1 (ru) | Способ обращения для определения добротности геологической среды | |
| US8755249B2 (en) | Automatic dispersion extraction of multiple time overlapped acoustic signals | |
| CN104849756B (zh) | 一种提高地震数据分辨率增强有效弱信号能量的方法 | |
| Tary et al. | Applications of high-resolution time-frequency transforms to attenuation estimation | |
| CN110187388B (zh) | 一种基于变分模态分解的稳定地震品质因子q估计方法 | |
| Wang et al. | High-resolution wave-equation amplitude-variation-with-ray-parameter (AVP) imaging with sparseness constraints | |
| Tary et al. | Attenuation estimation using high resolution time–frequency transforms | |
| CN107132579A (zh) | 一种保地层结构的地震波衰减补偿方法 | |
| Liu et al. | Seismic quality factor estimation using frequency-dependent linear fitting | |
| CN109901224B (zh) | 一种地震资料低频信号保护压制噪声方法 | |
| Jiang et al. | A combined denoising method of empirical mode decomposition and singular spectrum analysis applied to Jason altimeter waveforms: A case of the Caspian Sea | |
| CN102854530B (zh) | 基于对数时频域双曲平滑的动态反褶积方法 | |
| CN117111155B (zh) | 一种基于集成框架的微地震数据去噪方法 | |
| CN114063151A (zh) | 高精度叠前地震数据衰减属性提取方法及系统 | |
| Huang et al. | Robust time-frequency analysis of seismic data using general linear chirplet transform | |
| CN112684501B (zh) | 一种基于谱比面积的q值估计方法与应用 | |
| Chopra et al. | Coherence attribute applications on seismic data in various guises | |
| Naghadeh et al. | Source wavelet phase extraction | |
| CN114371505A (zh) | 一种基于地震分频技术的多子波反演方法及系统 | |
| CN112764095B (zh) | 基于广义能量比的vsp数据q值计算方法及系统 | |
| CN109459788B (zh) | 地层品质因子计算方法及系统 | |
| Weishi et al. | Statistical denoising of signals in the S-transform domain | |
| Wang et al. | Structure-oriented edge-preserving smoothing based on accurate estimation of orientation and edges | |
| Wang et al. | Continuous time-varying Q-factor estimation method in the time-frequency domain |
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 | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20220218 |