CN107102356B - 基于ceemd的地震信号高分辨率处理方法 - Google Patents
基于ceemd的地震信号高分辨率处理方法 Download PDFInfo
- Publication number
- CN107102356B CN107102356B CN201710408170.7A CN201710408170A CN107102356B CN 107102356 B CN107102356 B CN 107102356B CN 201710408170 A CN201710408170 A CN 201710408170A CN 107102356 B CN107102356 B CN 107102356B
- Authority
- CN
- China
- Prior art keywords
- signal
- imf
- imf component
- seismic
- ceemd
- 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.)
- Expired - Fee Related
Links
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/32—Transforming one recording into another or one representation into another
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
本发明公开了一种基于CEEMD的地震信号高分辨率处理方法，先根据地震剖面的平均振幅谱确定阈值集合G，接着用CEEMD对单道地震信号进行分解得到IMF分量，利用相关系数‑阈值法确定阈值H，通过阈值H选取有效频段内的IMF分量进行叠加，再用阈值集合G对叠加后的IMF分量进行高分辨率处理，然后将处理后的IMF分量和剩余的IMF分量进行信号重构，得到高分辨率处理后的地震信号，最后逐道处理余下的地震信号得到高分辨率的地震剖面。本发明通过计算地震剖面的平均振幅谱确定阈值集合G来避免仅使用标准的地震数据引起的随机干扰；通过阈值H选取有效频段内的IMF分量进行叠加后高分辨率处理，减少噪声对地震信号分辨率的影响。
Description
技术领域
本发明涉及地震勘探领域，尤其涉及一种基于CEEMD的地震信号高分辨率处理方法。
背景技术
地震勘探是利用地下介质弹性和密度的差异，通过观测和分析人工地震产生的地震波在地下的传播规律，推断地下岩层的性质和形态的勘探方法。其中，上述的地震波在传播过程中，受大地滤波、周围噪声等影响，导致了地表接收到的地震信号分辨率降低的问题。而地震信号分辨率是地震勘探工作中获取地层细节信息的关键因素，对研究薄层或小的地质体有着重要的意义，因此针对地震信号的特性来合理有效地提高其分辨率就显得尤为重要。
目前，提高地震信号分辨率的方法：谱白化、反Q滤波、多尺度联合分析等几类。谱白化是通过展宽振幅谱来达到补偿目的的方法，可是该方法存在破坏各个频率间振幅空间关系的问题；反Q滤波可以补偿振幅衰减、频率损失和改善相位特性，但该方法的准确程度高度依赖品质因子Q值求取；多尺度联合分析方法是利用测井、井间地震等技术对地下得同一目标进行不同尺度性质的反映，通过他们之间的联合作用提高地震资料分辨率，可是该方法需要特殊的井中资料，缺乏应用的普遍性。
可见，现有的技术中尚无较好的提高地震信号分辨率的方法。为了达到对非线性、非稳态的数据进行自适应和多角度分析的目的，Huang等人提出经验模态分解(EmpiricalMode Decomposition,EMD)。EMD方法作为一种多尺度分解方法，算法直观简单，能够自适应将信号分解为一系列固有模态函数(Intrinsic Mode Function,IMF)。EMD方法对非平稳信号的分析具有很大的优势，但它的分解过程存在模态混叠现象，导致某一个IMF分量中包含不同尺度的信号，或者相似的尺度信号存在于不同的IMF中。基于以上问题，WU等人提出一种噪声辅助信号分析的EMD方法，即集合经验模态分解(Ensemble Empirical ModeDecomposition,EEMD)，利用白噪声功率谱密度在整个频域内均匀分布特征，向原始地震信号中加入一定的白噪声，使得信号在不同尺度上具有连续性，之后再进行EMD分解。该方法在继承EMD分解的优点的同时，可以有效的解决EMD方法存在的模态混叠问题，但它是利用增加平均次数来降低重构误差，这是以增加计算成本为代价，而且它不具有完备性，信号的重构过程中也存在残余噪声的影响。
Torres等人对EEMD方法作了进一步的改进，提出互补集合经验模态分解(Complementary Ensemble Empirical Mode Decomposition,CEEMD)。该方法主要是通过向原信号中添加正负白噪声分别进行EMD分解，这样就能够很好地消除重构信号过程中残余的白噪声，而且加入的噪声集合次数可以很低，可提高信号处理的计算效率。CEEMD方法在继承EMD分解的优点的同时克服其分解会引起模态混叠现象的问题，还几乎消除了EEMD中添加白噪声引起的残余噪声的影响，适合于对精确性要求较高的地震资料进行高分辨率处理。
发明内容
本发明的目的就在于提供一种解决上述问题，不破坏各个频率间振幅空间关系，减少噪声影响，能快速提高地震信号的分辨率，充分准确的展示其细节信息的基于CEEMD的地震信号高分辨率处理方法。
为了实现上述目的，本发明采用的技术方案是一种基于CEEMD的地震信号高分辨率处理方法，包括以下步骤：
(1)获取地震剖面，求取该剖面的平均振幅谱，分析平均振幅谱的有效频带范围，根据频带范围确定其低通、低截、高通、高截的阈值，排列成阈值集合G，同时可以参考随机选取的数道单道地震信号的振幅谱，对阈值G进行修正，这里的地震剖面由单道地震信号构成，所述单道地震信号为用于勘探地下岩层性质和形态的时域信号；
(2)对其中一单道地震信号进行CEEMD分解，得到多个IMF分量和剩余分量；
(3)根据CEEMD分解得到的IMF分量排列特征，将IMF分量分为有效频段内的IMF分量和剩余的IMF分量；
(4)对有效频段内的IMF分量进行叠加，并用阈值集合G对叠加后的信号进行高分辨率处理；
(5)将高分辨率处理后的IMF分量和剩余的IMF分量进行信号重构，得到高分辨率处理后的单道地震信号；
x′(t)是高分辨率处理后的单道地震信号；imf′(t)表示高分辨率处理后的有效频带内的IMF分量；imf(t)代表剩余的IMF分量；rn(t)为剩余分量；M为CEEMD分解得到的IMF分量的个数，m为有效频段内的IMF分量个数；
(6)重复步骤(2)-(5)，对剩余单道地震信号逐道进行处理，最后依据处理的顺序依次排列所有高分辨率处理后的单道地震信号，获得高分辨率的地震剖面，用于地震解释。
在依据地震剖面平均振幅谱的有效频带确定其低通、低截、高通、高截的阈值集合G时，同时参考随机选取的数道地震信号振幅谱，对阈值G进行修正。
作为优选：步骤(2)的具体方法为：
(21)向原始单道地震信号中加入n组白噪声，白噪声是以正、负对互补的方式加入，从而生成一对添加噪声后的信号：
式中：x(t)为原始单道地震信号；Ni(t)为白噪声；Mi1(t),Mi2(t)分别为加入正、负成对噪声后的信号对,得到的加噪信号个数为2n；
(22)对每个加噪信号进行EMD分解，各自得到一组IMF分量，其中第i个信号的第j个IMF分量表示为imfij(t)，并通过下式得到imfj(t)，
其中，imfj(t)表示CEEMD分解得到的第j个IMF分量；
(23)原始地震信号x(t)表示为M个IMF分量与剩余分量rn(t)之和，如下式所示
在CEEMD分解的过程中需要确定所添加高斯白噪声的幅值ε及CEEMD分解的次数N，添加的白噪声幅值ε过小或者分解的次数N过少，都不能改变极值点分布，达到均匀极值点分布的目的；当添加的白噪声幅值ε过大或者分解的次数N过多时，虽然能减少添加白噪声的影响，但是会增加运行的时间。
所以，为了实际的分解效果及计算效率，结合其他作者选择的参数与自己测试地震信号的经验，添加高斯白噪声的幅值ε＝0.1时，CEEMD分解的次数N＝100为宜；添加高斯白噪声的幅值ε＝0.02时，CEEMD分解的次数N＝500为宜。
作为优选：步骤(3)具体方法为：
(31)利用下式计算出各个IMF分量与原始地震信号的相关程度Q，即
式中Q(j)代表第j个IMF信号与原始地震信号的相关程度，
(32)通过以下公式计算阈值H，
其中M为CEEMD分解得到的IMF分量的个数；
(33)依据CEEMD分解得到的IMF分量满足从高频到低频的依次排列分布的特征，首先利用下式将IMF分量分成A、B两部分，即
QL＝maxQ(j)，j＝1，2，3，…，M
接着，利用下式将A部分的IMF分量分成A1、A2两部分，即
最终得到的B部分和A2部分内的IMF分量为有效频带内的IMF分量，A1部分的IMF分量为剩余的IMF分量。
因CEEMD方法与EMD方法一样具有二进滤波特性，分解出的IMF分量满足从高频到低频的分布，前面几个IMF分量主要包含信号中的高频成分，常常随机噪声包含其中，后面的IMF分量主要包含信号中的低频成分，一般低频成分以有用信息为主。所以在进行IMF分量分类的时候，按照其排列顺序，将其分为A部分、B部分，而不会出现混乱。
作为优选：步骤(3)具体方法为：
(31)计算出各个IMF分量与原始地震道信号的相关程度Q，即
式中Q(j)代表第j个IMF信号与原始地震信号的相关程度，
(32)通过以下公式计算出ρ值，得到阈值H；
ρH＝max(ρj),j＝2,3,…M
(33)依据CEEMD分解得到的IMF分量满足从高频到低频依次排列的特点，用阈值H对IMF分量进行分类，具体公式如下：
作为优选：步骤(4)具体为：
(41)原始地震信号记录为x(t)，变换到频率域为
其中A(f)为地震记录的振幅谱，φ(f)为相位谱，采用平滑滤波的方法，由A(f)估算出地震记录振幅谱的衰减曲线B(f),其公式为
其中H(n)为滤波器的平滑因子，采用三角滤波算子，其公式为
N为半滤波算子长，式中n为滤波算子样点序号；
(42)振幅补偿曲线C(f)为：
式中fLC、fLP、fHP、fHC、fNq分别为地震信号的低截、低通、高通、高截、Nyquist频率；T是与白噪系数U有关的稳定系数，即
(43)用地震记录的原振幅谱对振幅补偿曲线C(f)进行标定，同时对补偿的结果进行线性加权处理,具体按照下式：
K为高频增强系数，取值范围为K＞1，△f为频率采样间隔；
(44)将A′(f)变换回时域，得到振幅补偿后的最终结果为：
此为高分辨率处理的方法，先将原始信号变换到频域，接着在频域中进行高分辨率处理，最后返回到时域，得到高分辨率的地震信号。
本发明的关键点在于：
(1)通过地震剖面的平均振幅谱，确定高分辨率处理的阈值集合G；
(2)CEEMD方法对单道地震信号的分频处理更有优势；
(3)阈值H选取的多样性，通过阈值H对对IMF分量进行分类，选出地震信号有效频带内的IMF分量；
(4)只对有效频带内的IMF分量叠加后进行高分辨率处理。
与现有技术相比，本发明的优点在于：通过计算地震剖面的平均振幅谱确定阈值集合G来避免仅使用标准的地震数据引起的随机干扰；通过阈值H的甄别，选取有效频段内的IMF分量叠加后进行高分辨率处理，不破坏各个频率间振幅空间关系，同时减少噪声对地震信号分辨率的影响。
本发明基于CEEMD方法对地震信号进行高分辨率处理，CEEMD方法对单道地震信号的分频处理更有优势,有效的避免了分解过程中存在的模态混叠，几乎消除了添加白噪声引起的残余噪声的影响，同时提高了信号处理的计算效率，适合于对精确性要求较高的地震资料进行高分辨率处理。
地震剖面中单道地震信号经CEEMD分解后，利用相关系数-阈值法选出阈值H，通过阈值H甄别出有效频段内的IMF分量进行叠加，并用阈值集合G对其进行高分辨率处理，将处理后的IMF分量与剩余IMF分量重构得到高分辨率处理后的单道地震信号，最后将所有处理后的单道地震信号排列为地震剖面。整个过程具备严格的理论基础，可以有效的提高分辨率，为地震勘探开发中的构造识别提供更为精确的地震资料。
上述优点具体可以分为以下三点：
(1)因为CEEMD方法对信号有较好的分解效果，能够很好的抑制模态混叠，对噪声的中和效果较好，所以本发明向原信号中添加正负白噪声分别进行EMD分解，很好地消除重构信号过程中残余的白噪声，提高信号处理的计算效率。CEEMD方法在克服EMD分解模态混叠现象的同时，还几乎消除了EEMD中添加白噪声引起的残余噪声的影响，适合于对精确性要求较高的地震资料进行高分辨率处理。
(2)获得阈值集合G：通过计算地震剖面的平均振幅谱，同时参考随机选取的几道地震信号振幅谱情况，分析剖面振幅谱的有效频带范围，分别确定高分辨率处理时的阈值集合G，避免仅使用标准的地震数据引起的随机干扰；
(3)获得阈值H，利用相关系数-阈值法选出阈值H，通过阈值H的甄别，选取有效频段内的IMF分量叠加后进行高分辨率处理，不破坏各个频率间振幅空间关系，同时减少噪声对地震信号分辨率的影响。
附图说明
图1为本发明流程图；
图2a是实施例1中实际资料1地震剖面的平均振幅谱；
图2b为图2中实际资料1第1道地震信号的振幅谱；
图2c为图2中实际资料1第50道地震信号的振幅谱；
图2d为图2中实际资料1第210道地震信号的振幅谱；
图2e为图2中实际资料1第420道地震信号的振幅谱；
图3a为实施例1中仿真信号各组成成分的时域波形；
图3b为图3a的仿真信号；
图3c为EMD方法对图3b的分解效果图；
图3d为CEEMD方法对图3b的分解效果图；
图4a为CEEMD分解参数高斯白噪声幅值0.02，执行次数500的分解效果；
图4b为CEEMD分解参数高斯白噪声幅值0.1，执行次数100的分解效果；
图5a为实施例1中确定有效频段内的IMF分量的示意图；
图5b为实施例2中确定有效频段内的IMF分量的示意图；
图6a为实际资料2高分辨率处理前的地震剖面图；
图6b为如6a所示地区高分辨率处理后的地震剖面图；
图6c为如6a所示地区高分辨率处理前的频谱曲线图；
图6d为如6a所示地区高分辨率处理后的频谱曲线图；
图7a为实际资料3高分辨率处理前的地震剖面图；
图7b为如7a所示地区高分辨率处理后地震剖面图；
图7c为如7a所示地区高分辨率处理前的频谱曲线图；
图7d为如7a所示地区高分辨率处理后的频谱曲线图。
具体实施方式
下面将结合附图对本发明作进一步说明。
实施例1：参见图1到图5。
如图1所示的一种基于CEEMD的地震信号高分辨率处理方法，包括以下步骤：
(1)获取地震剖面，求取该剖面的平均振幅谱，分析平均振幅谱的有效频带范围，根据频带范围确定其低通、低截、高通、高截的阈值，排列成阈值集合G，同时可以参考随机选取的数道单道地震信号的振幅谱，对阈值G进行修正，这里的地震剖面由单道地震信号构成，所述单道地震信号为用于勘探地下岩层性质和形态的时域信号；
参见图2，实际资料1是采样点为3001，道数为461，采样时间为2ms的地震剖面。从图2a可以看出，实际资料剖面的平均振幅谱的频带范围为0-70Hz,从图2b、2c、2d、2e了解到随机抽取的4个单道(第1、50、210、420道)的振幅谱的频带范围也在0-70之间，说明平均振幅谱求取的合理性。但各个单道的主频、频带宽度各有差异，根本不易确定出高分辨率处理时的参数，通过剖面的平均振幅谱的图像，我们可以确定高分辨率处理时的阈值集合G为[05 55 70]。实际效果表明，处理后的地震资料的分辨率较之前有很大的提高。而且通过求取地震剖面的平均振幅谱来确定高分辨率处理的阈值，可以避免仅使用标准的地震数据引起的随机干扰。
(2)对其中一单道地震信号进行CEEMD分解，得到多个IMF分量和剩余分量，具体方法为：
(21)向原始单道地震信号中加入n组白噪声，白噪声是以正、负对互补的方式加入，从而生成一对添加噪声后的信号：
式中：x(t)为原始单道地震信号；Ni(t)为白噪声；Mi1(t),Mi2(t)分别为加入正、负成对噪声后的信号对,得到的加噪信号个数为2n；
(22)对每个加噪信号进行EMD分解，各自得到一组IMF分量，其中第i个信号的第j个IMF分量表示为imfij(t)，并通过下式得到imfj(t)，
其中，imfj(t)表示CEEMD分解得到的第j个IMF分量；
(23)原始地震信号x(t)表示为M个IMF分量与剩余分量rn(t)之和，如下式所示
分量imfj(t)包含了信号从高到低的不同频率段的成分，剩余分量rn(t)取值为平均趋势或常量。
CEEMD方法对信号有较好的分解效果，能够很好的抑制模态混叠，对噪声的中和效果较好。为了说明该方法的优势，现对一仿真信号进行CEEMD分解，仿真信号各组成成分的时域波形如图3a所示，仿真信号如图3b所示，图3c是EMD方法对该仿真信号分解的效果，图3d是CEEMD分解得到的各个IMF分量的时域波形。通过对比分析可知，EMD分解方法的模态混叠效果比较严重，没有将原始仿真信号的3个信号分解出来，CEEMD分解结果中IMF1-3分别对应组成原始仿真信号的三个子信号。可以看出CEEMD消除了模态混叠的现象，由分解得到的间歇信号为零的部分可以看出，原信号所加的白噪声基本消除，噪声中和的效果较好。而且在整个的复分解过程中，CEEMD用时较短，大大提高了运算效率。
地震信号经由CEEMD分解得到的多个IMF分量，前面几个IMF分量主要包含信号中的高频成分，随机噪声常常包含其中，后面的IMF分量主要包含信号中的低频成分，一般低频成分以有用信息为主。CEEMD分解时选择的参数不同，其分解的效果也有差异。现对实际资料1的第33道进行分解CEEMD分解，当选择高斯白噪声的幅值ε＝0.1时，执行CEEMD分解次数N＝100，得到分解之后的9个IMF分量(见图4a)，当选择高斯白噪声的幅值ε＝0.02时，执行CEEMD分解次数N＝500，得到分解之后的10个IMF分量(见图4b)。通过图4可知，因为参数的选择不同，其分解的效果存在差异。或者在CEEMD分解的过程中可以选择其他的参数数值，本发明实施例不作限制，只要能够得到CEEMD分解的最佳效果即可。
(3)依据CEEMD分解方法具有完备性特点和CEEMD分解得到的IMF分量排列特征，利用相关系数-阈值法确定阈值H，通过阈值H的甄别，将IMF分量分为有效频段内的IMF分量和剩余的IMF分量，方法示意图见图5a所示，具体方法为：
(31)利用下式计算出各个IMF分量与原始地震信号的相关程度Q，即
式中Q(j)代表第j个IMF信号与原始地震信号的相关程度，
(32)通过以下公式计算阈值H，
其中M为CEEMD分解得到的IMF分量的个数；
(33)依据CEEMD分解得到的IMF分量满足从高频到低频的依次排列分布的特征，首先利用下式将IMF分量分成A、B两部分，即
QL＝maxQ(j)，j＝1，2，3，…，M
接着，利用下式将A部分的IMF分量分成A1、A2两部分，即
最终得到的B部分和A2部分内的IMF分量为有效频带内的IMF分量，图5中的C部分所示，A1部分的IMF分量为剩余的IMF分量。
(4)对有效频段内的IMF分量进行叠加，并用阈值集合G对叠加后的信号进行高分辨率处理，具体分为以下步骤：
(41)原始地震信号记录为x(t)，变换到频率域为
其中A(f)为地震记录的振幅谱，φ(f)为相位谱，采用平滑滤波的方法，由A(f)估算出地震记录振幅谱的衰减曲线B(f),其公式为
其中H(n)为滤波器的平滑因子，采用三角滤波算子，其公式为
N为半滤波算子长，式中n为滤波算子样点序号；
(42)振幅补偿曲线C(f)为：
式中fLC、fLP、fHP、fHC、fNq分别为地震信号的低截、低通、高通、高截、Nyquist频率；T是与白噪系数U有关的稳定系数，即
(43)用地震记录的原振幅谱对振幅补偿曲线C(f)进行标定，同时对补偿的结果进行线性加权处理,具体按照下式：
K为高频增强系数，取值范围为K＞1，△f为频率采样间隔；
(44)将A′(f)变换回时域，得到振幅补偿后的最终结果为：
(5)将高分辨率处理后的IMF分量和剩余的IMF分量进行信号重构，得到高分辨率处理后的单道地震信号；
x′(t)是高分辨率处理后的单道信号；imf′(t)表示高分辨率处理后的有效频带内的IMF分量；imf(t)代表剩余的IMF分量；rn(t)为剩余分量；
(6)重复步骤(2)-(5)，对剩余单道地震信号逐道进行处理，最后依据处理的顺序依次排列所有高分辨率处理后的单道地震信号，获得高分辨率的地震剖面，用于地震解释。
实施例2：步骤(3)中，根据IMF分量的特征，利用相关系数-阈值法确定阈值H，通过阈值H的甄别，将IMF分量分为有效频段内的IMF分量和其他的IMF分量，方法示意图见图5b所示,具体方法为：
(31)计算出各个IMF分量与原始地震道信号的相关程度Q，即
式中Q(j)代表第j个IMF信号与原始地震信号的相关程度，
(32)通过以下公式计算出ρ值，得到阈值H；
ρH＝max(ρj),j＝2,3,…M
(33)依据CEEMD分解得到的IMF分量满足从高频到低频依次排列的特点，用阈值H对IMF分量进行分类，具体公式如下：
其余与实施例1相同。
实施例3：参见图6a-图6d，图6a为实际资料1高分辨率处理前的剖面图，6c为该地区高分辨率处理前的频谱曲线图，图6b为实际资料1高分辨率处理后的剖面图，6d为该地区高分辨率处理后的频谱曲线图。通过实际资料1处理前后的剖面图和频谱曲线图分析，可以看出通过本发明实施例中的技术方案，整个剖面的分辨率明显提高，在处理前的剖面图中难以直接分辨的信息变得清晰，而且区分出一些薄层地区，地震剖面的频谱带宽明显变宽，主频显著提高。
实施例4：参见图7a-图7d，图7a为实际资料2高分辨率处理前的剖面图，7c为该地区高分辨率处理前的频谱曲线图，图7b为实际资料2高分辨率处理后的剖面图，7d为该地区高分辨率处理后的频谱曲线图。通过对比处理前后的效果图可知，通过本发明实施例中的技术方案，地震资料的平面分布特征与原始资料保持一致，层间信息变得丰富，细节的刻画更加清晰。地震剖面的频谱带宽明显变宽，主频显著提高。
Claims (3)
1.一种基于CEEMD的地震信号高分辨率处理方法，其特征在于,包括以下步骤：
(1)获取地震剖面，求取该剖面的平均振幅谱，分析平均振幅谱的有效频带范围，根据频带范围确定其低通、低截、高通、高截的阈值，排列成阈值集合G，地震剖面由单道地震信号构成，所述单道地震信号为用于勘探地下岩层性质和形态的时域信号；
(2)对其中一单道地震信号进行CEEMD分解，得到多个IMF分量和剩余分量；
(3)根据CEEMD分解得到的IMF分量排列特征，利用相关系数-阈值法确定阈值H，通过阈值H的甄别，将IMF分量分为有效频段内的IMF分量和剩余的IMF分量；利用相关系数-阈值法确定阈值H能通过以下两种方案：
方案一：(31)利用下式计算出各个IMF分量与原始地震信号的相关程度Q，即
式中Q(j)代表第j个IMF信号与原始地震信号的相关程度，
(32)通过以下公式计算阈值H，
其中M为CEEMD分解得到的IMF分量的个数；
(33)依据CEEMD分解得到的IMF分量满足从高频到低频的依次排列分布的特征，首先利用下式将IMF分量分成A、B两部分，即
QL=maxQ(j)，j=1，2，3，...，M
接着，利用下式将A部分的IMF分量分成A1、A2两部分，即
最终得到的B部分和A2部分内的IMF分量为有效频带内的IMF分量，A1部分的IMF分量为剩余的IMF分量；
方案二：(31)计算出各个IMF分量与原始地震道信号的相关程度Q，即
式中Q(j)代表第j个IMF信号与原始地震信号的相关程度，
(32)通过以下公式计算出ρ值，得到阈值H；
ρH＝max(ρj),j＝2,3,…M
(33)依据CEEMD分解得到的IMF分量满足从高频到低频依次排列的特点，用阈值H对IMF分量进行分类，具体公式如下：
(4)对有效频段内的IMF分量进行叠加，并用阈值集合G对叠加后的信号进行高分辨率处理；
(5)将高分辨率处理后的IMF分量和剩余的IMF分量进行信号重构，得到高分辨率处理后的单道地震信号；
x′(t)是高分辨率处理后的单道地震信号；imf′(t)表示高分辨率处理后的有效频带内的IMF分量；imf(t)代表剩余的IMF分量；rn(t)为剩余分量；M为CEEMD分解得到的IMF分量的个数，m为有效频段内的IMF分量个数；
(6)重复步骤(2)-(5)，对剩余单道地震信号逐道进行处理，最后依据处理的顺序依次排列所有高分辨率处理后的单道地震信号，获得高分辨率的地震剖面，用于地震解释。
2.根据权利要求1所述的基于CEEMD的地震信号高分辨率处理方法，其特征在于：步骤(2)中，具体方法为：
(21)向原始单道地震信号中加入n组白噪声，白噪声是以正、负对互补的方式加入，从而生成一对添加噪声后的信号：
式中：x(t)为原始单道地震信号；Ni(t)为白噪声；Mi1(t),Mi2(t)分别为加入正、负成对噪声后的信号对,得到的加噪信号个数为2n；
(22)对每个加噪信号进行EMD分解，各自得到一组IMF分量，其中第i个信号的第j个IMF分量表示为imfij(t)，并通过下式得到imfj(t)，
其中，imfj(t)表示CEEMD分解得到的第j个IMF分量；
(23)原始地震信号x(t)表示为M个IMF分量与剩余分量rn(t)之和，如下式所示
3.根据权利要求1所述的基于CEEMD的地震信号高分辨率处理方法，其特征在于：步骤(4)具体为：
(41)原始地震信号记录为x(t)，变换到频率域为
其中A(f)为地震记录的振幅谱，φ(f)为相位谱，采用平滑滤波的方法，由A(f)估算出地震记录振幅谱的衰减曲线B(f),其公式为
其中H(n)为滤波器的平滑因子，采用三角滤波算子，其公式为
N为半滤波算子长，式中n为滤波算子样点序号；
(42)振幅补偿曲线C(f)为：
式中fLC、fLP、fHP、fHC、fNq分别为地震信号的低截、低通、高通、高截、Nyquist频率；T是与白噪系数U有关的稳定系数，即
(43)用地震记录的原振幅谱对振幅补偿曲线C(f)进行标定，同时对补偿的结果进行线性加权处理,具体按照下式：
K为高频增强系数，取值范围为K＞1，△f为频率采样间隔；
(44)将A′(f)变换回时域，得到振幅补偿后的最终结果为：
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710408170.7A CN107102356B (zh) | 2017-06-02 | 2017-06-02 | 基于ceemd的地震信号高分辨率处理方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710408170.7A CN107102356B (zh) | 2017-06-02 | 2017-06-02 | 基于ceemd的地震信号高分辨率处理方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107102356A CN107102356A (zh) | 2017-08-29 |
| CN107102356B true CN107102356B (zh) | 2019-01-11 |
Family
ID=59659977
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710408170.7A Expired - Fee Related CN107102356B (zh) | 2017-06-02 | 2017-06-02 | 基于ceemd的地震信号高分辨率处理方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107102356B (zh) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107422381B (zh) * | 2017-09-18 | 2019-07-02 | 西南石油大学 | 一种基于eemd-ica的地震低频信息流体预测方法 |
| CN109557578A (zh) * | 2017-09-27 | 2019-04-02 | 中国石油化工股份有限公司 | 一种储层含气性检测方法及装置 |
| CN108333629B (zh) * | 2018-01-24 | 2021-08-24 | 中国矿业大学 | 一种利用经验模态分解和支持向量机定量预测煤厚的方法 |
| CN108303738A (zh) * | 2018-02-05 | 2018-07-20 | 西南石油大学 | 一种基于hht-mfcc的地震声纹流体预测方法 |
| CN108596144A (zh) * | 2018-05-09 | 2018-09-28 | 大连海事大学 | 一种分解参数自适应确定的互补集成经验模态分解方法 |
| CN108983286B (zh) * | 2018-07-23 | 2019-11-15 | 中国石油大学(华东) | 一种联合ceemd与广义s变换的地震数据去噪方法 |
| CN110908001B (zh) * | 2019-12-16 | 2020-09-29 | 吉林大学 | 一种大地电磁测深信号的重构方法及系统 |
| CN112200069B (zh) * | 2020-09-30 | 2022-11-04 | 山东大学 | 时频域谱减法和经验模态分解联合的隧道滤波方法及系统 |
| CN115144647B (zh) * | 2022-08-30 | 2022-12-30 | 国网江西省电力有限公司电力科学研究院 | 一种过电压智能识别方法及系统 |
| CN119644430A (zh) * | 2025-02-18 | 2025-03-18 | 中国地震局地质研究所 | 一种地震信号去噪方法、装置、介质及设备 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| RU2563323C1 (ru) * | 2014-04-02 | 2015-09-20 | Открытое акционерное общество "Нефтяная компания "Роснефть" | Способ реконструкции тонкой структуры геологического объекта и прогноза его флюидонасыщения |
| CN105700020A (zh) * | 2016-03-23 | 2016-06-22 | 中国石油天然气集团公司 | 一种地震数据随机噪声压制方法及装置 |
-
2017
- 2017-06-02 CN CN201710408170.7A patent/CN107102356B/zh not_active Expired - Fee Related
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| RU2563323C1 (ru) * | 2014-04-02 | 2015-09-20 | Открытое акционерное общество "Нефтяная компания "Роснефть" | Способ реконструкции тонкой структуры геологического объекта и прогноза его флюидонасыщения |
| CN105700020A (zh) * | 2016-03-23 | 2016-06-22 | 中国石油天然气集团公司 | 一种地震数据随机噪声压制方法及装置 |
Non-Patent Citations (3)
| Title |
|---|
| CEEMD-TK方法及其在致密砂岩储层含气性检测中的应用;徐丹 等;《中国地球科学联合学术年会2016》;20161231;第1466-1468页 |
| 振幅谱补偿和相位校正;姚逢昌;《石油物探》;19900930;第29卷(第3期);第46-59页 |
| 提高迭后地震记录分辨率的频率振幅补偿方法;程荃 等;《物探化探计算技术》;19960831;第18卷(第3期);第238-241页 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107102356A (zh) | 2017-08-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107102356B (zh) | 基于ceemd的地震信号高分辨率处理方法 | |
| Neelamani et al. | Coherent and random noise attenuation using the curvelet transform | |
| US11880011B2 (en) | Surface wave prediction and removal from seismic data | |
| Yuan et al. | Multiscale adjoint waveform tomography for surface and body waves | |
| Serdyukov et al. | Slant f-k transform of multichannel seismic surface wave data | |
| CN102998706B (zh) | 一种衰减地震数据随机噪声的方法及系统 | |
| Parolai et al. | The importance of converted waves in comparing H/V and RSM site response estimates | |
| EP0873527A1 (en) | Time-frequency processing and analysis of seismic data using very short-time fourier transforms | |
| Liu et al. | Improving vertical resolution of vintage seismic data by a weakly supervised method based on cycle generative adversarial network | |
| Fu | Joint inversion of seismic data for acoustic impedance | |
| Chen et al. | Seismic resolution enhancement by frequency-dependent wavelet scaling | |
| CA2741543A1 (en) | Simultaneous multiple source extended inversion | |
| CN106873036A (zh) | 一种基于井震结合的去噪方法 | |
| Jiang et al. | An adaptive generalized S-transform algorithm for seismic signal analysis | |
| Zheng et al. | Microseismic event denoising via adaptive directional vector median filters | |
| CN110109179A (zh) | 频宽补偿处理方法、装置和设备 | |
| Naskar et al. | Spectral whitening based seismic data preprocessing technique to improve the quality of surface wave's velocity spectra | |
| Jiang et al. | Seismic wavefield information extraction method based on adaptive local singular value decomposition | |
| Prajapati et al. | Delineation of stratigraphic pattern using combined application of wavelet-Fourier transform and fractal dimension: A case study over Cambay Basin, India | |
| Singh et al. | Modelling discontinuous well log signal to identify lithological boundaries via wavelet analysis: an example from KTB borehole data | |
| Abbasi et al. | Using geometric mode decomposition for the background noise suppression on microseismic data | |
| Steeghs et al. | Time-frequency analysis of seismic sequences | |
| CN107678065B (zh) | 提高地震分辨率的保构造井控空间反褶积方法和装置 | |
| Huang et al. | A fast least-squares reverse time migration method using cycle-consistent generative adversarial network | |
| Zhang et al. | Multi‐model stacked structure based on particle swarm optimization for random noise attenuation of 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 | ||
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20210128 Address after: Room 205C, 2nd floor, building D, 2-2, Beijing Shichuang hi tech Development Corporation Patentee after: WOKENSHI ENERGY TECHNOLOGY (BEIJING) Co.,Ltd. Address before: Three road 610059 Sichuan city of Chengdu province Chenghua District Erxian Qiaodong No. 1 Patentee before: Chengdu University of Technology |
|
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190111 |