CN108445539B - 一种消除地震子波旁瓣干扰的方法、设备及系统 - Google Patents
一种消除地震子波旁瓣干扰的方法、设备及系统 Download PDFInfo
- Publication number
- CN108445539B CN108445539B CN201810186640.4A CN201810186640A CN108445539B CN 108445539 B CN108445539 B CN 108445539B CN 201810186640 A CN201810186640 A CN 201810186640A CN 108445539 B CN108445539 B CN 108445539B
- Authority
- CN
- China
- Prior art keywords
- spectrum
- amplitude
- butterworth
- broadband
- sub
- 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 59
- 238000001228 spectrum Methods 0.000 claims abstract description 197
- 230000003044 adaptive effect Effects 0.000 claims abstract description 32
- 238000004422 calculation algorithm Methods 0.000 claims description 33
- 230000003595 spectral effect Effects 0.000 claims description 22
- 238000009499 grossing Methods 0.000 claims description 16
- 230000006870 function Effects 0.000 claims description 10
- 238000004590 computer program Methods 0.000 claims description 4
- 230000006978 adaptation Effects 0.000 abstract 1
- 230000008030 elimination Effects 0.000 abstract 1
- 238000003379 elimination reaction Methods 0.000 abstract 1
- 238000010586 diagram Methods 0.000 description 35
- 238000012545 processing Methods 0.000 description 16
- 230000008569 process Effects 0.000 description 11
- 230000004044 response Effects 0.000 description 10
- 230000006872 improvement Effects 0.000 description 9
- 238000005311 autocorrelation function Methods 0.000 description 6
- 239000010410 layer Substances 0.000 description 5
- 238000013508 migration Methods 0.000 description 5
- 239000004576 sand Substances 0.000 description 5
- 238000012935 Averaging Methods 0.000 description 4
- 230000005012 migration Effects 0.000 description 4
- 230000008859 change Effects 0.000 description 3
- 238000013213 extrapolation Methods 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 238000004321 preservation Methods 0.000 description 2
- 238000003860 storage Methods 0.000 description 2
- 241001123248 Arma Species 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000005314 correlation function Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000008021 deposition Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000005055 memory storage Effects 0.000 description 1
- 230000035772 mutation Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000001373 regressive effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 239000002356 single layer Substances 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
- 230000007480 spreading Effects 0.000 description 1
- 238000003892 spreading Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 239000013598 vector 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
- G01V1/364—Seismic filtering
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/30—Noise handling
- G01V2210/32—Noise reduction
- G01V2210/324—Filtering
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
本申请实施方式公开了一种消除地震子波旁瓣干扰的方法、设备及系统，其中，消除地震子波旁瓣干扰的方法包括：根据地震数据获得地震记录的振幅谱；根据所述地震记录的振幅谱获得振幅谱的振幅包络；构建频率域宽带Butterworth子波谱，并根据所述频率域宽带Butterworth子波谱获得宽带Butterworth子波谱包络；根据所述振幅谱的振幅包络和所述宽带Butterworth子波谱包络获得自适应加权因子；利用所述自适应加权因子将所述地震记录的振幅谱沿着所述宽带Butterworth子波谱包络进行载波调制，消除地震子波旁瓣干扰。
Description
技术领域
本申请涉及地震数据技术领域，特别涉及一种消除地震子波旁瓣干扰的方法、设备及系统。
背景技术
随着勘探程度的提高，我国已经从寻找构造性油气藏转变为寻找隐蔽性油气藏。其中，由沉积环境直接形成油气聚集的岩性油气藏发育着很多的薄层或薄互层储层，它们在不同岩性薄层间隔出现，具有相邻层地震波速度不同、以致地震记录中薄互层单层顶、底反射系数符号相异、走时差小等特点，由于地震勘探中激发子波频带有限，薄互层有效反射波呈现为子波相互叠加、相互干涉的结果。这种干涉叠加削弱薄互层界面的真实反射系数使地震记录成为一组整体地震响应，难以分辨。因此，砂泥岩薄互层成为地震资料处理和储层预测的重点与难点之一，需要发展一种识别薄互层的具有保幅性的高分辨率处理方法。能够综合考虑子波频率变化、补偿和保护地震子波高频信息，从而提高薄互层成像精度。
发明内容
本申请实施方式的目的是提供一种消除地震子波旁瓣干扰的方法、装置及系统，使用零相位宽带Butterworth子波约束谱消除子波旁瓣干扰影响，提高地震数据薄互层分辨率，解决如何提高地震图像质量和地震资料分辨率的技术问题。
为实现上述目的，本申请实施方式提供一种消除地震子波旁瓣干扰的方法，包括：
根据地震数据获得地震记录的振幅谱；
根据所述地震记录的振幅谱获得振幅谱的振幅包络；
构建频率域宽带Butterworth子波谱，并根据所述频率域宽带Butterworth子波谱获得宽带Butterworth子波谱包络；
根据所述振幅谱的振幅包络和所述宽带Butterworth子波谱包络获得自适应加权因子；
利用所述自适应加权因子将所述地震记录的振幅谱沿着所述宽带Butterworth子波谱包络进行载波调制，消除地震子波旁瓣干扰。
优选地，所述地震记录的振幅谱通过Burg最大熵谱估计方法获得。
优选地，所述振幅谱的振幅包络通过多点滑动平均方法获得。
优选地，所述频率域宽带Butterworth子波谱的表达式为：
其中，Pf、Qf均为频率参数；D(f)表示频率域宽带Butterworth子波谱。
优选地，所述自适应加权因子通过NLMS自适应滤波器算法获得。
优选地，所述宽带Butterworth子波谱包络通过基于Burg最大熵谱的振幅多点平滑及谱平滑因子获得。
为实现上述目的，本申请实施方式还提供一种消除地震子波旁瓣干扰的装置，包括：
振幅谱获得单元，用于根据地震数据获得地震记录的振幅谱；
振幅包络获得单元，用于根据所述地震记录的振幅谱获得振幅谱的振幅包络；
宽带Butterworth子波谱包络获得单元，用于构建频率域宽带Butterworth子波谱，并根据所述频率域宽带Butterworth子波谱获得宽带Butterworth子波谱包络；
自适应加权因子获得单元，用于根据所述振幅谱的振幅包络和所述宽带Butterworth子波谱包络获得自适应加权因子；
载波调制单元，用于利用所述自适应加权因子将所述地震记录的振幅谱沿着所述宽带Butterworth子波谱包络进行载波调制，消除地震子波旁瓣干扰。
优选地，所述振幅谱获得单元通过Burg最大熵谱估计方法获得地震记录的振幅谱。
优选地，所述振幅包络获得单元通过多点滑动平均方法获得振幅谱的振幅包络。
优选地，所述宽带Butterworth子波谱包络获得单元获得的频率域宽带Butterworth子波谱的表达式为：
其中，Pf、Qf均为频率参数；D(f)表示频率域宽带Butterworth子波谱。
优选地，所述自适应加权因子获得单元通过NLMS自适应滤波器算法获得自适应加权因子。
优选地，所述宽带Butterworth子波谱包络获得单元通过基于Burg最大熵谱的振幅多点平滑及谱平滑因子获得所述宽带Butterworth子波谱包络。
为实现上述目的，本申请实施方式还提供一种消除地震子波旁瓣干扰的系统，所述系统包括：存储器和处理器，所述存储器中存储计算机程序，所述计算机程序被所述处理器执行时，实现以下功能：
根据地震数据获得地震记录的振幅谱；
根据所述地震记录的振幅谱获得振幅谱的振幅包络；
构建频率域宽带Butterworth子波谱，并根据所述频率域宽带Butterworth子波谱获得宽带Butterworth子波谱包络；
根据所述振幅谱的振幅包络和所述宽带Butterworth子波谱包络获得自适应加权因子；
利用所述自适应加权因子将所述地震记录的振幅谱沿着所述宽带Butterworth子波谱包络进行载波调制，消除地震子波旁瓣干扰。
由上可见，与现有技术相比较，本技术方案基于零相位宽带Butterworth子波约束谱的自适应高分辨率处理地震数据，可以实现在压制面波噪声的同时能够有效地消除由于子波旁瓣造成的分辨率降低的问题，并且，在提高地震资料信噪比的同时也很好地提高了地震资料的分辨率。由于采用Burg最大熵谱估计和沃什变换算法使地震资料处理速度有很大的提高。
附图说明
为了更清楚地说明本申请实施方式或现有技术中的技术方案，下面将对实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本申请实施方式提供一种消除地震子波旁瓣干扰的方法流程图；
图2为Ricker子波滤波器响应曲线图；
图3为带通子波滤波器响应曲线图；
图4为宽带Butterworth子波滤波器响应曲线图；
图5为10-40Hz宽带Butterworth子波的曲线图；
图6为10-40Hz宽带Butterworth子波的频谱的曲线图；
图7a为宽带PSD的真实功率谱密度示意图；
图7b为窄带PSD的真实功率谱密度示意图；
图8a为Burg法(宽带)重叠的结果示意图；
图8b为Burg法(宽带)平均的结果示意图；
图9a为Burg法(窄带)重叠的结果示意图；
图9b为Burg法(窄带)平均的结果示意图；
图10a为西北新疆某区二维地震叠后处理剖面示意图；
图10b为本技术方案处理之后的西北新疆某区二维地震剖面示意图；
图10c为未使用本方案处理的频谱示意图；
图10d为使用本方案处理的频谱示意图；
图11a为西南四川某区三维地震叠后处理剖面示意图；
图11b为本技术方案处理之后的西南四川某区三维地震剖面示意图；
图11c为未使用本方案处理的西南四川某区三维地震频谱示意图；
图11d为使用本方案处理的西南四川某区三维地震频谱示意图；
图12a为国外哈萨克斯坦某油田三维叠前时间剖面未使用本方案处理的地震剖面示意图；
图12b为国外哈萨克斯坦某油田三维叠前时间剖面使用本方案处理的地震剖面示意图；
图12c为国外哈萨克斯坦某油田三维叠前时间剖面未使用本方案处理的频谱示意图；
图12d为国外哈萨克斯坦某油田三维叠前时间剖面使用本方案处理的频谱示意图；
图13a为中东部冀东某区三维叠前时间偏移原始地震剖面示意图；
图13b为中东部冀东某区三维叠前时间偏移CGG提频处理之后的地震剖面示意图；
图13c为中东部冀东某区三维叠前时间偏移剖面使用本方案处理的地震剖面示意图；
图14为本申请实施方式提供一种消除地震子波旁瓣干扰的装置功能框图；
图15为本申请实施例提出的一种消除地震子波旁瓣干扰的系统示意图。
具体实施方式
为了使本技术领域的人员更好地理解本申请中的技术方案，下面将结合本申请实施方式中的附图，对本申请实施方式中的技术方案进行清楚、完整地描述，显然，所描述的实施方式仅仅是本申请一部分实施方式，而不是全部的实施方式。基于本申请中的实施方式，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施方式，都应当属于本申请保护的范围。
为了拓宽地震记录优势频率的带宽，并有效减少子波旁瓣，突出高频弱反射信息，达到改善地震图像的质量，提高地震资料分辨率的目的，使用宽带Butterworth子波约束自适应谱“载波调制”的方案，达到消除子波旁瓣干扰影响的同时提高地震资料分辨率处理的技术效果。
地震数据处理的理想目标是得到反射系数，即子波是冲激函数。但由于受原始数据的频带和信噪比的限制，最终只能得到具有带限子波的记录。常用的子波是带通子波和Ricker子波。前者延续时间长，旁瓣波形复杂，后者旁瓣幅度比较大，都不是冲激函数的很好的近似。为此，设计了一种新的子波，它由不同宽度的Ricker子波合成，称为宽带Butterworth子波。它的主瓣比较窄、旁瓣幅度小、波形简单，在同样的主瓣宽度的情况下，振幅谱的峰值频率比较低，因而在保真度和信噪比方面都优于常用的带通子波和Ricker子波。
如图2所示，为Ricker子波滤波器响应曲线图。Ricker子波的表达式为：
其振幅谱为：
式中，g为Ricker子波表征频率的一个参数，等于峰值频率。
如图3所示，带通子波滤波器响应曲线图。带通子波的表达式为：
B(t)＝sinπ(fh-fl)tcosπ(fh+fl)t (3)
其振幅谱为
式中，fl和fh分别为低截频和高截频。
时间域宽带Butterworth子波表达式：
频率域宽带Butterworth子波表达式：
高分辨率宽带Ricker子波期望输出算子中，Pf、Qf为其表征频率参数，改变Pf影响谱的左支(低频)，改变Qf影响谱的右支(高频)。
如图4所示，为宽带Butterworth子波滤波器响应曲线图。可以看出相比带通滤波器，宽带Butterworth滤波器低频响应较陡，有利于地震记录中面波的衰减，高频响应较缓且有较宽的通放带，有利于地震记录的高分辨率处理。
如图5所示，为10-40Hz宽带Butterworth子波的曲线图。如图6所示，为10-40Hz宽带Butterworth子波的频谱的曲线图。宽带Butterworth子波最大的特点是最小相位和物理上可实现的。Butterworth子波在时间为零处开始，其在通频带内最平直，需要4个参数限定，从通带的上限频率和下限频率开始。然后是所需的上限和下限的截止陡度，即通带外的振幅响应的陡度。
基于上述描述，本申请实施方式提供一种消除地震子波旁瓣干扰的方法，如图1所示。包括：
步骤101)：根据地震数据获得地震记录的振幅谱。
在本实施例中，假设地震数据序列为x(n)，对于M阶AR(Auto Regressive model)模型有：
式中，为序列x(n)的估值，e(n)为估计误差，k表示AR模型的阶数。
在本实施例中，所述地震记录的振幅谱通过Burg最大熵谱估计方法获得。地震记录的振幅谱Gx(fi)的表达式为：
式中，Gx(fi)为频谱；σ2为线性预测误差的方差；ak为M阶AR模型的系数；f为频率。
经典谱估计是用有限长信号进行的线性估计，具有分辨率低和旁瓣“泄漏”的问题。其根本原因是自相关函数加窗，这样克服这些问题必须对自相关函数进行外推。现代谱估计是针对经典谱估计方差性能较差、分辨率较低的缺点提出并逐渐发展起来的，属于非线性估计，具有较高的频率分辨率。现代谱估计分为参数模型谱估计和非参数模型谱估计。而参数模型谱估计主要有AR模型、MA模型、ARMA模型等，其中AR模型应用最多。
AR模型是最主要的三种参数模型之一，其参数估计可归结为求解一组线性方程。要利用AR模型进行功率谱估计，必须由Yule-Walker方程求得AR模型的参数。Burg提出直接建立在信号上的AR模型系数求解方法，以最大熵作为自相关函数外推的准则，其合理性在于对自相关函数的/约束最少，因而时间序列的随机性最大，熵最大，功率谱最平坦。
研究最大熵谱估计时，Levinson递推一直受制于反射系数Km的求出。而Burg算法秉承着使前、后向预测误差平均功率最小的基本思想，不直接估计AR模型的参数，而是先估计反射系数Km，再利用Levinson关系式求得AR模型的参数，继而得到功率谱估计。
在信息论中，对在观测到一个以概率Pk发生的事件X＝xk后能够得到的信息，称之为信息量I(xk)。信息量的平均称为熵：
Burg的熵连续型随机变量的熵：
H(X)＝∫p(x)I(x)dx＝-∫p(x)ln p(x)dx＝E{I(x)} (7)
Burg根据连续随机变量的熵定义功率谱S(ω)的熵为：
利用给定的2m+1个样本自相关函数k＝0,±1,…,±m估计功率谱时，使得功率谱熵H[S(ω)]最大，这就是Burg最大熵谱估计方法。
求功率谱密度S(ω)，使得S(ω)在约束条件：
之下，能够使其谱熵H[S(ω)]为最大。这一具有约束条件的优化问题可以用拉格朗日乘子法求解。
构造目标函数：
式中λk为拉格朗日乘子。求则有：
作变量代换αk＝λ-k，则式(11)可以写作：
根据Fejer-Riesz定理知，Burg最大熵功率谱与AR功率谱等价。
为了实现最大熵谱估计，需要确定AR参数的阶数和系数。使用线性预测方法递推计算不同阶数的预测器系数，然后比较各预测器的预测误差功率。一种求解AR参数方法是要计算各阶情况下的预测系数和误差功率，进而选择出比较合理的模型阶数，计算量很大，因此可以考虑从低阶模型的参数递推出高阶模型的参数，实现这一递推计算的基础就是Levinson递推。
最大熵方法同时使用前向预测和后向预测。所谓前向预测就是利用给定的m个数据y(n-m)，…，y(n-1)预测y(n)的值，并称之为m阶前向线性预测。
y(n)的m阶前向线性预测值记为定义为：
式中αm(i)表示m阶预测器的第i个系数。
利用给定的m个数据y(n-m+1)，…，y(n)预测y(n-m)的值称之为m阶后向线性预测，定义为：
式中，表示是αm(i)的复数共轭。
下面根据最小均方误差准则设计前向和后向预测滤波器。前、后向预测误差分别定义为：
根据正交性原理知，为了使预测值是y(n)的线性均方估计，则前向预测误差f(n)必须与已知的数据y(n-m)，…，y(n-1)正交，即有：
E{f(n)y*(n-k)}＝0,1≤k≤m (17)
将式(15)代入式(17)，立即得到一组法方程：
式中，Rx(k)＝E{y(n)y*(n-k)}是{y(n)}的自相关函数。
定义前向预测均方误差，即前向预测均方误差的输出功率为：
由式(18)和(19)得到：
求解式(20)，即可直接得到m阶前向预测滤波器的系数αm(1)，…，αm(m)。
正交性原理为：
由正交性原理得到：
由预测误差均方值得到：
由正交性原理和预测误差均方值得到：
已知第k阶满足：
目标是求第k+1阶：
式(25)向目标方程(26)推进(加一行加一列)将条件向目标转化，扩充条件方程组：
其中，
式(27)先进行行倒序得到式(28)：
式(28)再进行列倒序得到：
行列倒序后的矩阵为式(30)：
利用共轭对称性R*(-τ)＝R(τ)，得到扩充方程组的矩阵式(31)：
行列倒序和共轭对称的方程组(扩充方程组)：
与目标方程组(26)相比，缺少一个自由度αk+1,k+1，所以利用共轭对称性质创造一个条件：
根据Dk+γk+1Pk＝0，得：
Pk+1＝Pk+γk+1D*k＝(1-|γk+1|2)Pk (36)
αk+1,i＝αk,i+γk+1α*k,k+1-i (37)
式中，γk称为反射系数。于是，可由1阶推出2阶，2阶推出3阶……
Burg算法的基本思想是使前、后向预测误差的平均功率为最小。m阶前、后向预测误差为：
将式(37)分别代入式(39)和式(40)，经整理后，即得到前、后向预测误差的阶数递推公式：
fm+1(n)＝fm(n)+γm+1gm(n-1) (41)
定义阶前、后向预测误差的平均功率：
将阶数递推公式(41)和(42)代入式(43)，并令则：
式中，m＝1，2，....
归纳起来，计算前向预测滤波器系数的Burg算法如下：
步骤1)：计算预测误差功率的初始值和前、后向预测误差的初始值。
f0(n)＝g0(n)＝x(n) (46)
并令m＝1。
步骤2)：求反射系数。
步骤3)：计算前向预测滤波器系数。
步骤4)：计算预测误差功率。
Pk+1＝(1-|γk+1|2)Pk (49)
步骤5)：计算滤波器输出。
fk+1(n)＝fk(n)+γk+1gk(n-1) (50)
步骤6)：令计算m←m+1，并重复步骤2至步骤5，直到预测误差功率Pk+1不再明显减小。
如图7a所示，为宽带PSD的真实功率谱密度示意图。如图7b所示，为窄带PSD的真实功率谱密度示意图。如图8a所示，为Burg法(宽带)重叠的结果示意图。如图8b所示，为Burg法(宽带)平均的结果示意图。如图9a所示，为Burg法(窄带)重叠的结果示意图。如图9b所示，为Burg法(窄带)平均的结果示意图。选择两种AR(4)过程，分别为宽带过程和窄带过程。观测数据分别由下面两式表达：
x(n)＝1.352x(n-1)-1.338x(n-2)+0.662x(n-3)-0.240x(n-4)+w(n)；
x(n)＝2.760x(n-1)-3.809x(n-2)+2.654x(n-3)-0.924x(n-4)+w(n)；
其中，w(n)是一高斯白噪声，用Burg算法估计模型参数、阶数和功率谱。每个采样样本的长度为256，每次产生50个采样序列，分别估计出功率谱密度(PSD)，然后将这50个结果重叠画在一幅图中以观察方差的大小，同时将这50个结果的均值(实线)和x(n)的理论功率谱密度(虚线)画在另一幅图中以观察偏差的大小。
通过使用Burg算法仿真50次得出的功率谱图，可以看出Burg算法的鲁棒性，通过将估计的结果同实际的功率谱相比较，可以看出Burg算法的精度和分辨率都很高，除个别频率点外，估计值与真实值的偏差几乎为0，进一步说明其利用谱熵最大化准则进行自相关函数外推的优越性。
步骤102)：根据所述地震记录的振幅谱获得振幅谱的振幅包络。
在本实施例中，振幅谱的振幅包络的表达式为：
式中，E(fi)是振幅谱平滑包络线；fi为频率；L平滑算子长度。
在本实施例中，振幅谱的振幅包络通过多点滑动平均方法获得。滑动平均法就是相当于有一个固定长度为L的滑动窗口，沿离散时间序列滑动。每滑动一个采样间隔，窗口前面进入一个新的数据，窗口后面去掉一个旧的数据，这样在窗口中始终有L个“最新”的数据。只要每次在滑动后把窗口中的L个数据进行算术平均，就可得到一组经过滑动平均滤波的新序列，其表达式为：
式中，L为滑动窗口的宽度。
步骤103)：构建频率域宽带Butterworth子波谱，并根据所述频率域宽带Butterworth子波谱获得宽带Butterworth子波谱包络。
在本实施例中，所述频率域宽带Butterworth子波谱的表达式为：
其中，Pf、Qf均为频率参数；D(f)表示频率域宽带Butterworth子波谱。
所述宽带Butterworth子波谱包络通过基于Burg最大熵谱的振幅多点平滑及谱平滑因子获得。
步骤104)：根据所述振幅谱的振幅包络和所述宽带Butterworth子波谱包络获得自适应加权因子。
在本实施例中，自适应加权因子通过NLMS自适应滤波器算法获得。NLMS自适应滤波器算法是在LMS算法的基础上发展而来的，解决了LMS算法因为输入信号的突然改变而造成权系数的突变，从而影响系统稳定的问题。NLMS自适应滤波器算法的主要原理是基于最小均方误差准则，算法通过在滤波过程中不断调节数字滤波器的权系数使得系统的期望信号d(n)与输出信号y(n)之间的均方误差E{e2(n)}最小的算法。NLMS算法的数学表达式如下：
y(n)＝wT(n)x(n) (52)
e(n)＝d(n)-y(n) (53)
权系数的表达式，如下：
上式中，步长因子μ的取值为0＜μ＜2，NLMS自适应滤波器算法对输入信号进行了归一化，因而该算法相对于LMS算法有很好的鲁棒性，解决了输入信号突变带来的系统不稳定为问题。NLMS自适应滤波器算法就是在LMS算法基础上采用了可变因子，归一化输入向量，使瞬时误差最小化。
由于NLMS自适应滤波器算法的收敛条件与输入信号的特征值是无关的，同时增大了算法的动态输入范围，计算量又与LMS算法相差不多，且加之其具有变步长的特点使得NLMS自适应滤波器算法具有更快的收敛速度。
步骤105)：利用所述自适应加权因子将所述地震记录的振幅谱沿着所述宽带Butterworth子波谱包络进行载波调制，消除地震子波旁瓣干扰。
在本实施例中，采用沃尔什变换算法对地震数据进行FFT正变换，获得频域的地震数据。使用自适应加权因子w(n)实现频域的地震数据的“载波调制”。从频率域地震数据的实部和虚部获得的载波调制包络和振幅包络，将载波调制包络和振幅包络中各个频点的数值对应的自适应加权因子相乘，得到载波调制后的地震数据结果。采用沃尔什变换算法对载波调制后的地震数据结果实现IFFT变换，将频域数据转换到时间域数据，输出最终处理结果，该结果就是消除地震子波旁瓣干扰的时间域数据。
如图10a所示，为西北新疆某区二维地震叠后处理剖面示意图。如图10b所示，为本技术方案处理之后的地震剖面示意图。可以看出剖面的纵向分辨率有显著的提高。在图10b的剖面中上部可以看分辨出相对于图10a呈现更多的地质层位，并且各波组特征更加明显，更利于各地层层位横向连续追踪对比。在图10b的剖面下部所表现的绕射波特征更加明显，分辨率更高。如图10c所示，为未使用本方案处理的频谱示意图。如图10d所示，为使用本方案处理的频谱示意图。在图10c、图10d中，黑色曲线为整个剖面频谱，红色为选中区域频谱。由图10c所示，剖面的频带较窄，主频在10-30hz，处理后剖面频带较宽，主频在8-50hz，从频谱分析的结果可以看出。图10d所示，经本方案处理之后，可以将原始资料频带拓宽20Hz以上。使用本方案可以提高构造成图和砂体描述的精度能够加深对砂体纵横向展布规律的认识，同时可有效提高后续地震反演的精度，减少其多解性。
如图11a所示，为西南四川某区三维地震叠后处理剖面示意图。如图11b所示，为本技术方案处理之后的西南四川某区三维地震剖面示意图。可以看出地震剖面的分辨率有显著的提高，同时有效波组更加突出。各地层层位横向连续可等时追踪对比。如图11c所示，为未使用本方案处理的频谱示意图。如图11d所示，为使用本方案处理的频谱示意图。由图11c可知，拓频前剖面的频带很窄，低频及高频信息能量弱。由图11d可知，拓频后剖面频带得以较大幅度的拓宽，低频及高频信息能量得以明显提升，使地震波频带变宽，时间域分辨率提高，含较多低频成分的资料其反演结果分辨率较高、成像较好。从而提高了储层预测解释的正确性。
如图12a所示，为国外哈萨克斯坦某油田三维叠前时间剖面未使用本方案处理的地震剖面示意图。如图12b所示，为国外哈萨克斯坦某油田三维叠前时间剖面使用本方案处理的地震剖面示意图。可以看出提高分辨率的同时突出有效波组，地层反射续至减少，分辨率高。如图12c所示，为国外哈萨克斯坦某油田三维叠前时间剖面未使用本方案处理的频谱示意图。由图12c可知，拓频前剖面的频带很窄(主频40Hz-100Hz)，低频及高频信息能量弱。如图12d所示，为国外哈萨克斯坦某油田三维叠前时间剖面使用本方案处理的频谱示意图。由图12d可知，拓频后剖面频带(主频18Hz-110Hz)得以较大幅度的拓宽，低频及高频信息能量得以明显提升，有利于提高砂体的描述精度，开展沉积微相及砂体分布的研究。
如图13a所示，为中东部冀东某区三维叠前时间偏移原始地震剖面示意图。如图13b所示，为中东部冀东某区三维叠前时间偏移CGG提频处理之后的地震剖面示意图。由图13b可知，法国CGG公司地震处理软件反Q滤波提频后的剖面，此方法不保幅，波组续至多。如图13c所示，为中东部冀东某区三维叠前时间偏移剖面使用本方案处理的地震剖面示意图。由图13c可知，可以看出提高分辨率的同时突出有效波组，且具保幅特性。在剖面中下目的层部位，法国CGG公司软件提频后波组特征及振幅特性被破坏。本技术方案提频结果分辨率提高的同时振幅特性得到很好的保持，为下一步准确预测岩性油气藏提供一套可靠的基础数据。
如图14所示，为本申请实施例提出的一种消除地震子波旁瓣干扰的装置功能框图。包括包括：
振幅谱获得单元11，用于根据地震数据获得地震记录的振幅谱；
振幅包络获得单元12，用于根据所述地震记录的振幅谱获得振幅谱的振幅包络；
宽带Butterworth子波谱包络获得单元13，用于构建频率域宽带Butterworth子波谱，并根据所述频率域宽带Butterworth子波谱获得宽带Butterworth子波谱包络；
自适应加权因子获得单元14，用于根据所述振幅谱的振幅包络和所述宽带Butterworth子波谱包络获得自适应加权因子；
载波调制单元15，用于利用所述自适应加权因子将所述地震记录的振幅谱沿着所述宽带Butterworth子波谱包络进行载波调制，消除地震子波旁瓣干扰。
在本实施例中，所述振幅谱获得单元通过Burg最大熵谱估计方法获得地震记录的振幅谱。
在本实施例中，所述振幅包络获得单元通过多点滑动平均方法获得振幅谱的振幅包络。
在本实施例中，所述宽带Butterworth子波谱包络获得单元获得的频率域宽带Butterworth子波谱的表达式为：
其中，Pf、Qf均为频率参数；D(f)表示频率域宽带Butterworth子波谱。
在本实施例中，所述自适应加权因子获得单元通过NLMS自适应滤波器算法获得自适应加权因子。
在本实施例中，所述宽带Butterworth子波谱包络获得单元通过基于Burg最大熵谱的振幅多点平滑及谱平滑因子获得所述宽带Butterworth子波谱包络。
如图15所示，为本申请实施例提出的一种消除地震子波旁瓣干扰的系统示意图。所述系统包括：存储器a和处理器b，所述存储器a中存储计算机程序，所述计算机程序被所述处理器b执行时，实现以下功能：
根据地震数据获得地震记录的振幅谱；
根据所述地震记录的振幅谱获得振幅谱的振幅包络；
构建频率域宽带Butterworth子波谱，并根据所述频率域宽带Butterworth子波谱获得宽带Butterworth子波谱包络；
根据所述振幅谱的振幅包络和所述宽带Butterworth子波谱包络获得自适应加权因子；
利用所述自适应加权因子将所述地震记录的振幅谱沿着所述宽带Butterworth子波谱包络进行载波调制，消除地震子波旁瓣干扰。
在本实施方式中，所述存储器包括但不限于随机存取存储器(Random AccessMemory,RAM)、只读存储器(Read-Only Memory,ROM)、缓存(Cache)、硬盘(Hard DiskDrive,HDD)或者存储卡(Memory Card)。
在本实施方式中，所述处理器可以按任何适当的方式实现。例如，所述处理器可以采取例如微处理器或处理器以及存储可由该(微)处理器执行的计算机可读程序代码(例如软件或固件)的计算机可读介质、逻辑门、开关、专用集成电路(Application SpecificIntegrated Circuit，ASIC)、可编程逻辑控制器和嵌入微控制器的形式等等。
本说明书实施方式提供的消除地震子波旁瓣干扰的系统，其存储器和处理器实现的具体功能，可以与本说明书中的前述实施方式相对照解释，并能够达到前述实施方式的技术效果，这里便不再赘述。
在20世纪90年代，对于一个技术的改进可以很明显地区分是硬件上的改进(例如，对二极管、晶体管、开关等电路结构的改进)还是软件上的改进(对于方法流程的改进)。然而，随着技术的发展，当今的很多方法流程的改进已经可以视为硬件电路结构的直接改进。设计人员几乎都通过将改进的方法流程编程到硬件电路中来得到相应的硬件电路结构。因此，不能说一个方法流程的改进就不能用硬件实体模块来实现。例如，可编程逻辑器件(Programmable Logic Device,PLD)(例如现场可编程门阵列(Field Programmable GateArray，FPGA))就是这样一种集成电路，其逻辑功能由用户对器件编程来确定。由设计人员自行编程来把一个数字系统“集成”在一片PLD上，而不需要请芯片制造厂商来设计和制作专用的集成电路芯片。而且，如今，取代手工地制作集成电路芯片，这种编程也多半改用“逻辑编译器(logic compiler)”软件来实现，它与程序开发撰写时所用的软件编译器相类似，而要编译之前的原始代码也得用特定的编程语言来撰写，此称之为硬件描述语言(Hardware Description Language，HDL)，而HDL也并非仅有一种，而是有许多种，如ABEL(Advanced Boolean Expression Language)、AHDL(Altera Hardware DescriptionLanguage)、Confluence、CUPL(Cornell University Programming Language)、HDCal、JHDL(Java Hardware Description Language)、Lava、Lola、MyHDL、PALASM、RHDL(RubyHardware Description Language)等，目前最普遍使用的是VHDL(Very-High-SpeedIntegrated Circuit Hardware Description Language)与Verilog2。本领域技术人员也应该清楚，只需要将方法流程用上述几种硬件描述语言稍作逻辑编程并编程到集成电路中，就可以很容易得到实现该逻辑方法流程的硬件电路。
本领域技术人员也知道，除了以纯计算机可读程序代码方式实现客户端、服务器以外，完全可以通过将方法步骤进行逻辑编程来使得客户端、服务器以逻辑门、开关、专用集成电路、可编程逻辑控制器和嵌入微控制器等的形式来实现相同功能。因此这种客户端、服务器可以被认为是一种硬件部件，而对其内包括的用于实现各种功能的装置也可以视为硬件部件内的结构。或者甚至，可以将用于实现各种功能的装置视为既可以是实现方法的软件模块又可以是硬件部件内的结构。
通过以上的实施方式的描述可知，本领域的技术人员可以清楚地了解到本申请可借助软件加必需的通用硬件平台的方式来实现。基于这样的理解，本申请的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本申请各个实施方式或者实施方式的某些部分所述的方法。
本说明书中的各个实施方式均采用递进的方式描述，各个实施方式之间相同相似的部分互相参见即可，每个实施方式重点说明的都是与其他实施方式的不同之处。尤其，针对客户端的实施方式来说，均可以参照前述方法的实施方式的介绍对照解释。
本申请可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本申请，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。
虽然通过实施方式描绘了本申请，本领域普通技术人员知道，本申请有许多变形和变化而不脱离本申请的精神，希望所附的权利要求包括这些变形和变化而不脱离本申请的精神。
Claims (9)
1.一种消除地震子波旁瓣干扰的方法，其特征在于，包括：
根据地震数据获得地震记录的振幅谱；
根据所述地震记录的振幅谱获得振幅谱的振幅包络；
构建频率域宽带Butterworth子波谱，并根据所述频率域宽带Butterworth子波谱获得宽带Butterworth子波谱包络；
根据所述振幅谱的振幅包络和所述宽带Butterworth子波谱包络获得自适应加权因子；
利用所述自适应加权因子将所述地震记录的振幅谱沿着所述宽带Butterworth子波谱包络进行载波调制，消除地震子波旁瓣干扰；
其中，所述频率域宽带Butterworth子波谱的表达式为：
其中，Pf、Qf均为频率参数；D(f)表示频率域宽带Butterworth子波谱；
所述宽带Butterworth子波谱包络通过基于Burg最大熵谱的振幅多点平滑及谱平滑因子获得。
2.如权利要求1所述的方法，其特征在于，所述地震记录的振幅谱通过Burg最大熵谱估计方法获得。
3.如权利要求1所述的方法，其特征在于，所述振幅谱的振幅包络通过多点滑动平均方法获得。
4.如权利要求1所述的方法，其特征在于，所述自适应加权因子通过NLMS自适应滤波器算法获得。
5.一种消除地震子波旁瓣干扰的装置，其特征在于，包括：
振幅谱获得单元，用于根据地震数据获得地震记录的振幅谱；
振幅包络获得单元，用于根据所述地震记录的振幅谱获得振幅谱的振幅包络；
宽带Butterworth子波谱包络获得单元，用于构建频率域宽带Butterworth子波谱，并根据所述频率域宽带Butterworth子波谱获得宽带Butterworth子波谱包络；
自适应加权因子获得单元，用于根据所述振幅谱的振幅包络和所述宽带Butterworth子波谱包络获得自适应加权因子；
载波调制单元，用于利用所述自适应加权因子将所述地震记录的振幅谱沿着所述宽带Butterworth子波谱包络进行载波调制，消除地震子波旁瓣干扰；
其中，所述宽带Butterworth子波谱包络获得单元获得的频率域宽带Butterworth子波谱的表达式为：
其中，Pf、Qf均为频率参数；D(f)表示频率域宽带Butterworth子波谱；
所述宽带Butterworth子波谱包络获得单元通过基于Burg最大熵谱的振幅多点平滑及谱平滑因子获得所述宽带Butterworth子波谱包络。
6.如权利要求5所述的装置，其特征在于，所述振幅谱获得单元通过Burg最大熵谱估计方法获得地震记录的振幅谱。
7.如权利要求5所述的装置，其特征在于，所述振幅包络获得单元通过多点滑动平均方法获得振幅谱的振幅包络。
8.如权利要求5所述的装置，其特征在于，所述自适应加权因子获得单元通过NLMS自适应滤波器算法获得自适应加权因子。
9.一种消除地震子波旁瓣干扰的系统，其特征在于，所述系统包括：存储器和处理器，所述存储器中存储计算机程序，所述计算机程序被所述处理器执行时，实现以下功能：
根据地震数据获得地震记录的振幅谱；
根据所述地震记录的振幅谱获得振幅谱的振幅包络；
构建频率域宽带Butterworth子波谱，并根据所述频率域宽带Butterworth子波谱获得宽带Butterworth子波谱包络；
根据所述振幅谱的振幅包络和所述宽带Butterworth子波谱包络获得自适应加权因子；
利用所述自适应加权因子将所述地震记录的振幅谱沿着所述宽带Butterworth子波谱包络进行载波调制，消除地震子波旁瓣干扰；
其中，所述频率域宽带Butterworth子波谱的表达式为：
其中，Pf、Qf均为频率参数；D(f)表示频率域宽带Butterworth子波谱；
所述宽带Butterworth子波谱包络通过基于Burg最大熵谱的振幅多点平滑及谱平滑因子获得。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810186640.4A CN108445539B (zh) | 2018-03-07 | 2018-03-07 | 一种消除地震子波旁瓣干扰的方法、设备及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810186640.4A CN108445539B (zh) | 2018-03-07 | 2018-03-07 | 一种消除地震子波旁瓣干扰的方法、设备及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108445539A CN108445539A (zh) | 2018-08-24 |
| CN108445539B true CN108445539B (zh) | 2019-08-30 |
Family
ID=63193520
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810186640.4A Active CN108445539B (zh) | 2018-03-07 | 2018-03-07 | 一种消除地震子波旁瓣干扰的方法、设备及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108445539B (zh) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113933902B (zh) * | 2020-07-10 | 2024-07-02 | 中国石油化工股份有限公司 | 基于立体空间子波的高频拓展方法 |
| CN113359187B (zh) * | 2021-06-16 | 2022-08-02 | 王仰华 | 地震数据的子波旁瓣消除方法 |
| CN113420452B (zh) * | 2021-06-30 | 2022-04-01 | 中国工程物理研究院总体工程研究所 | 地基微振动设计载荷确定方法 |
| CN114384580B (zh) * | 2021-12-31 | 2023-05-02 | 同济大学 | 一种基于可控震源的理想子波定制方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| GB2359363B (en) * | 2000-02-15 | 2002-04-03 | Geco Prakla | Processing simultaneous vibratory seismic data |
| CN103364827A (zh) * | 2012-03-30 | 2013-10-23 | 中国石油化工股份有限公司 | 基于双参数目标寻优的自适应谱模拟反褶积方法 |
| CN104597502A (zh) * | 2014-12-08 | 2015-05-06 | 翟明岳 | 一种新的石油地震勘探数据去噪方法 |
| CN107179550A (zh) * | 2017-07-05 | 2017-09-19 | 西安交通大学 | 一种数据驱动的地震信号零相位反褶积方法 |
-
2018
- 2018-03-07 CN CN201810186640.4A patent/CN108445539B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| GB2359363B (en) * | 2000-02-15 | 2002-04-03 | Geco Prakla | Processing simultaneous vibratory seismic data |
| CN103364827A (zh) * | 2012-03-30 | 2013-10-23 | 中国石油化工股份有限公司 | 基于双参数目标寻优的自适应谱模拟反褶积方法 |
| CN104597502A (zh) * | 2014-12-08 | 2015-05-06 | 翟明岳 | 一种新的石油地震勘探数据去噪方法 |
| CN107179550A (zh) * | 2017-07-05 | 2017-09-19 | 西安交通大学 | 一种数据驱动的地震信号零相位反褶积方法 |
Non-Patent Citations (2)
| Title |
|---|
| 《几种常用解析子波的特性分析》;张海燕等;《石油地球物理勘探》;20071231;第42卷(第6期);全文 |
| 《谱模拟融合自适应分段和局部相似度的时变子波提取方法研究》;王蓉蓉;《中国石油大学硕士学位论文》;20160531;全文 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108445539A (zh) | 2018-08-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108445539B (zh) | 一种消除地震子波旁瓣干扰的方法、设备及系统 | |
| CN103995289B (zh) | 基于时频谱模拟的时变混合相位地震子波提取方法 | |
| Liu et al. | Random noise attenuation using f-x regularized nonstationary autoregression | |
| CN108549100B (zh) | 基于非线性高次拓频的时间域多尺度全波形反演方法 | |
| CN108919347A (zh) | 基于vmd的地震信号随机噪声压制方法 | |
| Ghaderpour et al. | Antileakage least-squares spectral analysis for seismic data regularization and random noise attenuation | |
| Li et al. | Multidimensional seismic data reconstruction using frequency-domain adaptive prediction-error filter | |
| CN103645507B (zh) | 地震记录的处理方法 | |
| CN111159888B (zh) | 一种基于互相关函数的协方差矩阵稀疏迭代时延估计方法 | |
| CN108318921B (zh) | 一种基于横向约束的快速地震随机反演方法 | |
| CN107179550B (zh) | 一种数据驱动的地震信号零相位反褶积方法 | |
| CN114200525B (zh) | 一种自适应的多道奇异谱分析地震数据去噪方法 | |
| Liu et al. | Irregularly sampled seismic data reconstruction using multiscale multidirectional adaptive prediction-error filter | |
| CN111399057B (zh) | 一种基于非凸稀疏约束的地震资料噪声压制方法 | |
| Kimiaefar et al. | Random noise attenuation by Wiener-ANFIS filtering | |
| CN112147678B (zh) | 一种高精度的时变子波反演方法 | |
| CN111474582A (zh) | 生成高精度时频谱的精准s变换方法 | |
| CN103308945A (zh) | 一种陆地勘探初至前噪声的模拟产生与预测方法 | |
| CN114779332B (zh) | 地震数据沉积背景去除方法、装置及电子设备 | |
| CN113341463B (zh) | 一种叠前地震资料非平稳盲反褶积方法及相关组件 | |
| Tang et al. | Non‐minimum phase seismic wavelet reconstruction based on higher order statistics | |
| CN104360322B (zh) | 基于四阶非对称乘积型核函数的qfm信号参数估计方法 | |
| DONG et al. | Fast implementation technique of adaptive Kalman filtering deconvolution via dyadic wavelet transform | |
| Wu et al. | A new method for high resolution well-control processing of post-stack seismic data | |
| CN109085649A (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 |