CN108107475A - 一种基于经验小波变换和多阈值函数的井中微地震去噪方法 - Google Patents
一种基于经验小波变换和多阈值函数的井中微地震去噪方法 Download PDFInfo
- Publication number
- CN108107475A CN108107475A CN201810177571.0A CN201810177571A CN108107475A CN 108107475 A CN108107475 A CN 108107475A CN 201810177571 A CN201810177571 A CN 201810177571A CN 108107475 A CN108107475 A CN 108107475A
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- mtd
- mtr
- signal
- 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
- 238000000034 method Methods 0.000 title claims abstract description 58
- 230000009466 transformation Effects 0.000 title claims abstract description 17
- 238000001228 spectrum Methods 0.000 claims abstract description 21
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 11
- 238000004364 calculation method Methods 0.000 claims description 7
- 230000003044 adaptive effect Effects 0.000 abstract description 13
- 230000000694 effects Effects 0.000 abstract description 10
- 238000006243 chemical reaction Methods 0.000 abstract description 4
- 238000003384 imaging method Methods 0.000 abstract description 4
- 238000002474 experimental method Methods 0.000 abstract description 3
- 238000012545 processing Methods 0.000 abstract description 3
- 238000000354 decomposition reaction Methods 0.000 abstract description 2
- 239000004615 ingredient Substances 0.000 abstract 3
- 238000000605 extraction Methods 0.000 abstract 1
- 230000005764 inhibitory process Effects 0.000 abstract 1
- 238000000926 separation method Methods 0.000 abstract 1
- 238000005204 segregation Methods 0.000 description 6
- 238000012544 monitoring process Methods 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 3
- 238000002592 echocardiography Methods 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 239000004215 Carbon black (E152) Substances 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 229930195733 hydrocarbon Natural products 0.000 description 1
- 150000002430 hydrocarbons Chemical class 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003204 osmotic effect Effects 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000010008 shearing Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000001131 transforming 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
-
- 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
- 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
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和小波变换WT已经广泛地应用于地震去噪，但是CEEMD缺少理论基础，WT的自适应性较弱。所以去噪效果较差。本文首次将经验小波变换结合多种阈值函数用于微地震去噪。EWT通过频谱分割建立自适应小波滤波器组提取被测信号的不同频率块。在EWT中有四种频谱分割方法，实验发现adaptive算法可以很好地将微地震数据的有效信号和噪声分离，EWT变换后，通过分析各模块的频谱和能量将信号分为两种成分。我们将硬阈值函数应用于包含较多的有效信号的成分，将改进的阈值函数应用于包含较少的有效信号的成分。将提出方法与CEEMD和WT在模拟信号与实际信号中进行对比，证明了提出方法的有效性。
Description
技术领域
本发明涉及地震勘探中的地震信号处理方法，尤其涉及一种基于经验小波变换和多阈值函数的井中微地震去噪方法。
背景技术
微地震监测是提高非常规油气藏天然渗透率的重要手段。微地震波在水力压裂过程中产生，具有幅度小、频率高、信噪比低的特点。由微地震信号的性质可知对微地震信号进行处理有一定的难度，不能直接应用于定位油气藏，需要先对微地震信号做去噪处理。过去30年来，水力压裂和采矿诱发微震事件的监测应用迅速增加。因此如何能有效地抑制噪声并恢复有效信号很重要。微地震监测分为井中监测和地面监测两种[4]。我们处理的是井中微地震信号。
微地震数据的准确分析和反演是在噪声衰减法的基础上发展起来的，现存的随机噪声抑制和反射事件恢复方法很多。一些方法需要将被测信号变换到其他域，如f-x反卷积、小波变换、τ-p变换和剪切波变换(shearlet)。f-x反卷积是一种预测误差滤波的方法，该方法基于在f-x域的AR模型抑制随机噪声。为了提高微地震信号的幅度，在时频域中提出了小波变换和阈值结合的方法。τ-p变换被应用于地面微地震信号去噪。Shearlet变换将信号变换到shearlet域，是一种新型的多尺度多方向的分析方法。一些去噪方法需要提取被测信号的主要成分，如经验模态分解(EMD)和追踪去噪方法。EMD根据被测信号的主要信息将信号分割。为了改进EMD，提出了EEMD和CEEMD。对于多组分微震数据引入了追踪去噪方法。然而，这些方法在弱反射事件中，特别是当SNR低于给定阈值时表现不佳。
发明内容
本发明设计开发了一种基于经验小波变换和多阈值函数的井中微地震去噪方法，通过经验小波算法抑制井中微地震的噪声，能有效地将微地震中的信噪分离，多种阈值函数的使用提高了去噪效果使得微地震事件变得更清晰。
本发明提供的技术方案为：
一种基于经验小波变换和多阈值函数的井中微地震去噪方法，包括：
步骤一：计算井中微地震波信号的快速傅里叶变换，得到被测信号的频谱，根据频谱确定极小值，并找到每个极小值对应的频率；
步骤二：首先设置一个初始边界，然后在初始边界周围找到一个区间，最后在区间找到最小的极小值对应的频率作为边界；
步骤三：根据所述边界，定义尺度函数和经验小波函数；
步骤四：对信号与尺度函数和经验小波函数进行卷积运算，获得地震信号的不同模块；
步骤五：计算每个模块能量值和平均模块能量值，对模块能量值大于平均值两倍的模块采用硬阈值函数滤噪，其余模块采用改进的阈值函数进行滤噪。
优选的是，所述步骤二包括：
定义0和π作为第一个边界和最后一个边界，频谱在[0,π]内被分为K个连续的部分，初始边界为ak,k＝1,2···,K-1ak∈(0,ω)，ω为信号频率的最大值，窗口的计算公式为：
其中，round(*)是四舍五入函数，lower和upper分别是区间的上限和下限；
所述边界的计算公式为：
bk-1`＝fmin+lower-2，k＝2,3···K-1
bk-1为第K-1个边界，fmin为最小的极小值对应的频率，所述最小的极小值在原始频谱上计算。
优选的是，所述尺度函数公式为：
所述经验小波函数公式为：
其中，N为频谱上最大值的个数，每个最大值对应的频率为ωi(i＝1,2,3···,N)；α(γ,ωi)＝β[(1/2γωi)-(1-γ)ωi]，γ是确保两个连续的变换不重叠的参数，β(x)是满足下一条件的任意函数
β(x)＝x4(35-84x+70x2-20x3)。
优选的是，所述尺度函数与被测信号内积得到尺度系数
经验小波函数与被测信号内积得到经验小波系数。
其中，W(k,t)被定义为在第t时刻第k个滤波器。
优选的是，硬阈值函数表达式
改进阈值函数表达式
其中，Wk,t是是第k个模块的经验小波系数；是阈值之后的系数；λk是与噪声方差相关的阈值；sgn(*)代表符号函数。
优选的是，能量计算公式为：
其中，是第k个模块的能量。
本发明所述的有益效果
本发明中提出了一种基于EWT的新的微地震去噪方法。EWT结合了CEEMD和WT的优点，是一种提取信号中不同成分的有效方法。由于井中微地震信号有低信噪比和高频率的特点，我们将EWT和多种阈值函数相结合保留信号的保幅度。在模拟信号和实际信号的实验结果中可以看出，对比CEEMD和WT，EWT算法在噪声抑制和有效信号保留方面表现更好。
附图说明
图1为本发明所述的基于经验小波变换和多阈值函数的井中微地震去噪方法的流程图。
图2为本发明所述的纯净信号、含噪信号和实际信号的频谱图。
图3为本发明所述的adaptive分隔方法的边界频谱图。
图4为本发明所述的adaptive分隔方法的滤波器组频谱图。
图5为本发明所述的Adaptive分隔方法下的IMF频谱图。
图6为本发明所述的Locmax分隔方法下的IMF频谱图。
图7为本发明所述的Locmaxmin分隔方法下的IMF频谱图。
图8为本发明所述的Locmaxminf分隔方法下的IMF频谱图。
图9为本发明所述的不同分割方法的去噪结果波形图。
图10为本发明所述的不同分割方法的去噪结果波形放大图。
图11为本发明所述的不同阈值的去噪结果。
图12为本发明所述的三种方法在第6道信号的去噪结果波形图。
图13为本发明所述的三种方法在第6道信号的去噪结果波形放大图。
图14为本发明所述的三种方法在第6道信号的去噪结果频谱图。
图15为本发明所述的三种方法的去噪结果纯净信号图。
图16为本发明所述的三种方法的去噪结果含噪信号图。
图17为本发明所述的CEEMD方法去噪结果图。
图18为本发明所述的CEEMD残差
图19为本发明所述的db5小波去噪结果。
图20为本发明所述的db5小波残差图。
图21为本发明所述的EWT去噪结果。
图22为本发明所述的EWT残差图。
图23为发明所述的三种方法的去噪结果原始信号。
图24为本发明所述的CEEMD去噪结果。
图25为本发明所述的CEEMD残差图。
图26为本发明所述的db5小波去噪结果。
图27为本发明所述的db5小波残差图。
图28为本发明所述的EWT去噪结果。
图29为本发明所述的EWT残差图。
具体实施方式
下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。
如图1所示，本发明提供的的基于经验小波变换和多阈值函数的井中微地震去噪方法，包括：
步骤一、计算被测地震波信号f(t)的快速傅里叶变换，得到被测信号的频谱F(ω)；
步骤二、适当地分隔频谱，将边界设为Ω＝{ωK}k＝1,2,···K，假设频谱在[0,π]内被分为K个连续的部分：
频谱分割方法是一个重要的步骤，它为EWT自适应的分析信号提供了一个理论支持。在EWT中有四种频谱分割方法，分别为‘locmax’、‘localmaxmin’、localmaxminf’和‘adaptive’。adaptive方法可以将边界设在100Hz和300Hz附近，而井中微地震信号的主频范围在100Hz到300Hz之间，所以adaptive方法更适用于井中微地震去噪。
adaptive方法的主要思想是：首先设置一个初始边界，然后在初始边界周围找到一个区间，最后在区间找到最小的极小值对应的频率作为边界。最小的极小值在原始频谱上计算。因为0和π通常作为第一个边界和最后一个边界，我们需要找到另外K-1个边界。假设初始边界为ak,k＝1,2···,K-1ak∈(0,ω)ω为信号频率的最大值，窗口的计算公式为：
其中，round(*)是四舍五入函数，lower和upper分别是区间的上限和下限。
边界的计算公式为：
bk-1`＝fmin+lower-2，k＝2,3···K-1
其中，bk-1为第K-1个边界，fmin为最小的极小值对应的频率。
步骤三：基于步骤二中的边界，应用Littlewood-Paley和Meyer小波结构，得到经验小波，尺度函数公式为：
经验小波函数公式为：
其中，α(γ,ωi)＝β[(1/2γωi)-(1-γ)ωi]，γ是确保两个连续的变换不重叠的参数，β(x)是满足下一条件的任意函数
满足上述条件的函数有很多，在文献中应用最多的是
β(x)＝x4(35-84x+70x2-20x3)
步骤四：根据尺度函数和经验小波函数提取出信号的不同模块。尺度函数与被测信号内积得到尺度系数
经验小波函数与被测信号内积得到经验小波系数。
其中，W(k,t)被定义为在第t时刻第k个滤波器。
步骤五、采用硬阈值函数和改进阈值函数结合的方式进一步去噪
我们假设Wk,t是是第k个模块的经验小波系数；是阈值之后的系数；λk是与噪声方差相关的阈值；sgn(*)代表符号函数。硬阈值函数和改进阈值函数的表达式如下
硬阈值函数表达式
改进阈值函数表达式
其中，
硬阈值函数已经广泛地应用于去噪并且可以很好地保留信号幅度。但是硬阈值函数有一些缺点，例如：在重构数据时在阈值点出现不连续的现象；改进的阈值函数改善硬阈值函数的缺点，但是对于有效信号幅度的保留劣于硬阈值函数，因为不同模块中含有的信息不同，应用单一的阈值函数达不到最好的去噪效果，模块能量值越大包含的有效信号越多，因此对于有效信号多的模块采用硬阈值函数去噪，其他模块采用改进的阈值函数进一步去噪。
如图2-8所示，为了证明adaptive方法的可行性，我们用加入-5dB噪声，主频为200Hz的雷克子波模拟井中微地震信号。图2展示了纯净信号、含噪信号和实际信号的频谱，可以看出信号的主频在100Hz到300Hz之间。我们可以观察到adaptive方法可以将边界设在100Hz和300Hz附近，即将微地震信号的主频(有效信号)分到一个模块中。
如图9-10中展示了基于四种分割方法的去噪结果。从去噪结果中可以看出adaptive方法在噪声抑制和有效信号保留方面优于其他三种方法。
计算每个模块能量值和平均模块能量值，对模块能量值大于平均值两倍的模块采用硬阈值函数滤噪，其余模块采用改进的阈值函数进行滤噪。实施以硬阈值函数和改进阈值函数的方法为例做进一步说明众所周知，模块的能量值越大包含的有效信号越多。根据表1，我们可以发现IMF2中含有的有效信号最多，IMF1和IMF3中也含有少量的有效信号。因为不同模块中含有的信息不同，应用单一的阈值函数达不到最好的去噪效果。我们在IMF2中应用硬阈值函数保留有效信号的幅度，在IMF1和IMF3在应用改进的阈值函数进一步去噪。
表1中应用的能量公式如下。
其中，是第k个模块的能量。
| IMF1 | IMF2 | IMF3 | |
| 信号1 | 0.0027 | 0.0061 | 0.0037 |
| 信号2 | 0.0014 | 0.0095 | 0.0022 |
| 信号3 | 0.0027 | 0.0078 | 0.0038 |
| 信号4 | 0.0019 | 0.0070 | 0.0040 |
表1 每个模块的能量
对于每个IMF，阈值λk定义为
其中，ck是通过实验得到的主要参数；σk是噪声标准差，N是离散信号采样点的数目。
因为有效信号的幅度要高于噪声幅度，所述我们使用计算公式估计噪声标准差，标准差估算公式为：
其中，median(*)代表取中值。
图11是不同阈值函数的去噪结果。明显地，我们提出的方法比单独使用硬阈值函数和改进阈值函数去噪效果好。
实验结果
如图12-14，为了验证EWT算法的可行性和有效性，我们模拟了15道主频为200Hz的微地震信号，其波速分别为400和500m/s，采样频率为1024Hz，信噪比为-5dB，检波器之间的距离为5m。EWT算法与WT和CEEMD对比如下。对于单道微地震信号，三种方法的对比结果展示在中。提出的方法能更好地在波峰和波谷处保留有效信号。与此同时，EWT算法抑制噪声的效果更好。在频谱方面，EWT算法的去噪结果与纯净信号在低频和高频处更接近。
如图14-21，三种方法应用于交叉信道。EWT得到了更满意地输出，微地震事件变得更清晰。CEEMD和WT结果中依然含有不同程度的噪声。
信噪比SNR和均方误差MSE是评判去噪效果的两个重要参数。在表2中，可以看出EWT表现效果最好。SNR和MSE的表达式如下
Fk,t是纯净信号，是去噪后的信号。
表2 不同方法去噪后SNR和MES
为了验证提出方法的实用性，将方法应用于实际的井中微地震信号。实际信号来源于中国某地区。我们选取其中的15道进行实验。图24-29中展示了三种方法的对比结果。
本文中提出了一种基于EWT的新的微地震去噪方法。EWT结合了CEEMD和WT的优点，是一种提取信号中不同成分的有效方法。由于井中微地震信号有低信噪比和高频率的特点，我们将EWT和多种阈值函数相结合保留信号的饱幅度。在模拟信号和实际信号的实验结果中可以看出，对比CEEMD和WT，EWT算法在噪声抑制和有效信号保留方面表现更好。
图24是原始实际信号，其主频在100Hz到300Hz之间。为了展示EWT算法的优越性，我们将三种方法在区域A和B中作比较。区域A中包含了较多的噪声，区域B中包含较多的有效信号。图25和图26分别是CEEMD和WT的去噪结果，可以看出结构中依然有噪声残留。从残差图中可以看出，CEEMD和WT有信号的丢失，而提出方法的效果最好。通过对比，提出方法能更好地提高信噪比并维持信号的连续性。
尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。
Claims (6)
1.一种基于经验小波变换和多阈值函数的井中微地震去噪方法，其特征在于，包括：
步骤一：计算井中微地震波信号的快速傅里叶变换，得到被测信号的频谱，根据频谱确定极小值，并找到每个极小值对应的频率；
步骤二：首先设置一个初始边界，然后在初始边界周围找到一个区间，最后在区间找到最小的极小值对应的频率作为边界；
步骤三：根据所述边界，定义尺度函数和经验小波函数；
步骤四：对信号与尺度函数和经验小波函数进行卷积运算，获得地震信号的不同模块；
步骤五：计算每个模块能量值和平均模块能量值，对模块能量值大于平均值两倍的模块采用硬阈值函数滤噪，其余模块采用改进的阈值函数进行滤噪。
2.根据权利要求1所述的基于经验小波变换和多阈值函数的井中微地震去噪方法，其特征在于，所述步骤二包括：
定义0和π作为第一个边界和最后一个边界，频谱在[0,π]内被分为K个连续的部分，初始边界为ak,k＝1,2…,K-1；ak∈(0,ω)，ω为信号频率的最大值，窗口的计算公式为：
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mi>l</mi>
<mi>o</mi>
<mi>w</mi>
<mi>e</mi>
<mi>r</mi>
<mo>=</mo>
<msub>
<mi>a</mi>
<mi>k</mi>
</msub>
<mo>-</mo>
<mi>r</mi>
<mi>o</mi>
<mi>u</mi>
<mi>n</mi>
<mi>d</mi>
<mo>(</mo>
<mfrac>
<mrow>
<msub>
<mi>a</mi>
<mi>k</mi>
</msub>
<mo>-</mo>
<msub>
<mi>a</mi>
<mrow>
<mi>k</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
</mrow>
<mn>2</mn>
</mfrac>
<mo>)</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>u</mi>
<mi>p</mi>
<mi>p</mi>
<mi>e</mi>
<mi>r</mi>
<mo>=</mo>
<msub>
<mi>a</mi>
<mi>k</mi>
</msub>
<mo>+</mo>
<mi>r</mi>
<mi>o</mi>
<mi>u</mi>
<mi>n</mi>
<mi>d</mi>
<mo>(</mo>
<mfrac>
<mrow>
<msub>
<mi>a</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>-</mo>
<msub>
<mi>a</mi>
<mi>k</mi>
</msub>
</mrow>
<mn>2</mn>
</mfrac>
<mo>)</mo>
</mtd>
</mtr>
</mtable>
</mfenced>
其中，round(*)是四舍五入函数，lower和upper分别是区间的上限和下限；
所述边界的计算公式为：
bk-1`＝fmin+lower-2，k＝2,3…K-1
bk-1为第K-1个边界，fmin为最小的极小值对应的频率，所述最小的极小值在原始频谱上计算。
3.根据权利要求2所述的基于经验小波变换和多阈值函数的井中微地震去噪方法，其特征在于，所述尺度函数公式为：
所述经验小波函数公式为：
<mrow>
<msub>
<mi>&psi;</mi>
<mi>k</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>&omega;</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mn>1</mn>
<mo>,</mo>
</mrow>
</mtd>
<mtd>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>+</mo>
<mi>&gamma;</mi>
<mo>)</mo>
<msup>
<mi>&omega;</mi>
<mi>k</mi>
</msup>
<mo>&le;</mo>
<mo>|</mo>
<mi>&omega;</mi>
<mo>|</mo>
<mo>&le;</mo>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<mi>&gamma;</mi>
<mo>)</mo>
<msup>
<mi>&omega;</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>cos</mi>
<mo>&lsqb;</mo>
<mfrac>
<mi>&pi;</mi>
<mn>2</mn>
</mfrac>
<mi>&alpha;</mi>
<mrow>
<mo>(</mo>
<mi>&gamma;</mi>
<mo>,</mo>
<msup>
<mi>&omega;</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>,</mo>
</mrow>
</mtd>
<mtd>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<mi>&gamma;</mi>
<mo>)</mo>
<msup>
<mi>&omega;</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>&le;</mo>
<mo>|</mo>
<mi>&omega;</mi>
<mo>|</mo>
<mo>&le;</mo>
<mo>(</mo>
<mn>1</mn>
<mo>+</mo>
<mi>&gamma;</mi>
<mo>)</mo>
<msup>
<mi>&omega;</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>sin</mi>
<mo>&lsqb;</mo>
<mfrac>
<mi>&pi;</mi>
<mn>2</mn>
</mfrac>
<mi>&alpha;</mi>
<mrow>
<mo>(</mo>
<mi>&gamma;</mi>
<mo>,</mo>
<msup>
<mi>&omega;</mi>
<mi>k</mi>
</msup>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>,</mo>
</mrow>
</mtd>
<mtd>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<mi>&gamma;</mi>
<mo>)</mo>
<msup>
<mi>&omega;</mi>
<mi>k</mi>
</msup>
<mo>&le;</mo>
<mo>|</mo>
<mi>&omega;</mi>
<mo>|</mo>
<mo>&le;</mo>
<mo>(</mo>
<mn>1</mn>
<mo>+</mo>
<mi>&gamma;</mi>
<mo>)</mo>
<msup>
<mi>&omega;</mi>
<mi>k</mi>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mn>0</mn>
<mo>,</mo>
</mrow>
</mtd>
<mtd>
<mrow>
<mi>o</mi>
<mi>t</mi>
<mi>h</mi>
<mi>e</mi>
<mi>r</mi>
<mi>w</mi>
<mi>i</mi>
<mi>s</mi>
<mi>e</mi>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
其中，N为频谱上最大值的个数，每个最大值对应的频率为ωi(i＝1,2,3…,N)；γ是确保两个连续的变换不重叠的参数，β(x)是满足下一条件的任意函数
<mrow>
<mi>&beta;</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mrow>
<mi>x</mi>
<mo>&le;</mo>
<mn>0</mn>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mrow>
<mi>x</mi>
<mo>&GreaterEqual;</mo>
<mn>1</mn>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>&beta;</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>&beta;</mi>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mn>1</mn>
</mrow>
</mtd>
<mtd>
<mrow>
<mi>x</mi>
<mo>&Element;</mo>
<mrow>
<mo>(</mo>
<mn>0</mn>
<mo>,</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>;</mo>
</mrow>
β(x)＝x4(35-84x+70x2-20x3)。
4.根据权利要求1所述的基于经验小波变换和多阈值函数的井中微地震去噪方法，其特征在于，所述尺度函数与被测信号内积得到尺度系数
经验小波函数与被测信号内积得到经验小波系数。
<mrow>
<mi>W</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>,</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mo><</mo>
<mi>f</mi>
<mo>,</mo>
<msub>
<mi>&psi;</mi>
<mi>k</mi>
</msub>
<mo>></mo>
<mo>=</mo>
<mo>&Integral;</mo>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mi>&tau;</mi>
<mo>)</mo>
</mrow>
<mo>&CenterDot;</mo>
<mover>
<mrow>
<msub>
<mi>&psi;</mi>
<mi>k</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>&tau;</mi>
<mo>-</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>&OverBar;</mo>
</mover>
<mo>&CenterDot;</mo>
<mi>d</mi>
<mi>&tau;</mi>
</mrow>
其中，W(k,t)被定义为在第t时刻第k个滤波器。
5.根据权利要求4所述的基于经验小波变换和多阈值函数的井中微地震去噪方法，其特征在于，硬阈值函数表达式
<mrow>
<msub>
<mover>
<mi>W</mi>
<mo>^</mo>
</mover>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
</msub>
<mo>=</mo>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<msub>
<mi>W</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
</msub>
</mtd>
<mtd>
<mrow>
<mo>|</mo>
<msub>
<mi>W</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
</msub>
<mo>|</mo>
<mo>&GreaterEqual;</mo>
<msub>
<mi>&lambda;</mi>
<mi>k</mi>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mrow>
<mo>|</mo>
<msub>
<mi>W</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
</msub>
<mo>|</mo>
<mo><</mo>
<msub>
<mi>&lambda;</mi>
<mi>k</mi>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
改进阈值函数表达式
<mrow>
<msub>
<mover>
<mi>W</mi>
<mo>^</mo>
</mover>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
</msub>
<mo>=</mo>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<mi>&mu;</mi>
<mo>)</mo>
<msub>
<mi>W</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
</msub>
<mo>+</mo>
<mi>&mu;</mi>
<mi>sgn</mi>
<mo>(</mo>
<msub>
<mi>W</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
</msub>
<mo>)</mo>
<mo>(</mo>
<mo>|</mo>
<msub>
<mi>W</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
</msub>
<mo>|</mo>
<mo>-</mo>
<mi>&lambda;</mi>
<mo>)</mo>
</mrow>
</mtd>
<mtd>
<mrow>
<mo>|</mo>
<msub>
<mi>W</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
</msub>
<mo>|</mo>
<mo>&GreaterEqual;</mo>
<msub>
<mi>&lambda;</mi>
<mi>k</mi>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mrow>
<mo>|</mo>
<msub>
<mi>W</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
</msub>
<mo>|</mo>
<mo><</mo>
<msub>
<mi>&lambda;</mi>
<mi>k</mi>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
其中，Wk,t是是第k个模块的经验小波系数；是阈值之后的系数；λk是与噪声方差相关的阈值；sgn(*)代表符号函数。
6.根据权利要求1所述的基于经验小波变换和多阈值函数的井中微地震去噪方法，其特征在于，能量计算公式为：
<mrow>
<msubsup>
<mi>E</mi>
<mi>k</mi>
<mn>2</mn>
</msubsup>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>n</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>N</mi>
</munderover>
<mo>|</mo>
<msub>
<mi>W</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
</msub>
<msup>
<mo>|</mo>
<mn>2</mn>
</msup>
</mrow>
其中，是第k个模块的能量。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810177571.0A CN108107475B (zh) | 2018-03-05 | 2018-03-05 | 基于经验小波变换和多阈值函数的井中微地震去噪方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810177571.0A CN108107475B (zh) | 2018-03-05 | 2018-03-05 | 基于经验小波变换和多阈值函数的井中微地震去噪方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108107475A true CN108107475A (zh) | 2018-06-01 |
| CN108107475B CN108107475B (zh) | 2019-07-19 |
Family
ID=62205938
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810177571.0A Expired - Fee Related CN108107475B (zh) | 2018-03-05 | 2018-03-05 | 基于经验小波变换和多阈值函数的井中微地震去噪方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108107475B (zh) |
Cited By (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108956764A (zh) * | 2018-06-06 | 2018-12-07 | 西安理工大学 | 一种爆炸复合管结合状态定量识别方法 |
| CN109270345A (zh) * | 2018-07-20 | 2019-01-25 | 国网江苏省电力有限公司常州供电分公司 | 一种电网谐波信号的检测方法 |
| CN109520611A (zh) * | 2018-11-08 | 2019-03-26 | 温州大学 | 一种地震模拟振动台工况的监测方法 |
| CN110671613A (zh) * | 2019-10-15 | 2020-01-10 | 重庆邮电大学 | 基于改进经验小波变换的流体管道泄漏信号时延估计方法 |
| CN111238362A (zh) * | 2020-01-20 | 2020-06-05 | 西安理工大学 | 一种深孔钻削中孔圆度误差实时监测的方法 |
| CN111830559A (zh) * | 2019-04-19 | 2020-10-27 | 中国石油天然气股份有限公司 | 一种地震资料拓频处理的方法及装置 |
| CN112329667A (zh) * | 2020-11-11 | 2021-02-05 | 杭州电子科技大学 | 基于经验小波变换和多尺度熵的水质数据清洗方法 |
| CN112558147A (zh) * | 2019-09-25 | 2021-03-26 | 中国石油化工股份有限公司 | 一种井中微地震数据的偏振分析方法及系统 |
| CN112633225A (zh) * | 2020-12-31 | 2021-04-09 | 矿冶科技集团有限公司 | 矿用微震信号滤波方法 |
| CN112740262A (zh) * | 2018-09-14 | 2021-04-30 | 莱卡微系统Cms有限责任公司 | 使用局部信噪比的去卷积装置和方法 |
| CN112764108A (zh) * | 2019-11-05 | 2021-05-07 | 防灾科技学院 | 一种基于改进经验小波变换的新型地震资料噪声压制算法 |
| CN113158797A (zh) * | 2021-03-16 | 2021-07-23 | 南方科技大学 | 一种微震资料去噪方法、系统、电子设备及存储介质 |
| CN114114406A (zh) * | 2020-08-26 | 2022-03-01 | 中国石油化工股份有限公司 | 一种储层渗透率估计方法及装置 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20150160358A1 (en) * | 2013-12-11 | 2015-06-11 | Raytheon Bbn Technologies Corp. | Systems and methods for combining triaxial geophone data for localizing nearby transient seismic sources |
| CN107450097A (zh) * | 2017-07-21 | 2017-12-08 | 杰奥世博（北京）技术有限公司 | 基于多子波的地震数据频谱均等分频方法 |
| CN107479093A (zh) * | 2017-09-18 | 2017-12-15 | 中南大学 | 一种基于势函数的微地震事件去噪和聚类方法 |
| CN107589454A (zh) * | 2017-07-25 | 2018-01-16 | 西安交通大学 | 一种基于vmd‑tfpf压制地震勘探随机噪声方法 |
-
2018
- 2018-03-05 CN CN201810177571.0A patent/CN108107475B/zh not_active Expired - Fee Related
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20150160358A1 (en) * | 2013-12-11 | 2015-06-11 | Raytheon Bbn Technologies Corp. | Systems and methods for combining triaxial geophone data for localizing nearby transient seismic sources |
| CN107450097A (zh) * | 2017-07-21 | 2017-12-08 | 杰奥世博（北京）技术有限公司 | 基于多子波的地震数据频谱均等分频方法 |
| CN107589454A (zh) * | 2017-07-25 | 2018-01-16 | 西安交通大学 | 一种基于vmd‑tfpf压制地震勘探随机噪声方法 |
| CN107479093A (zh) * | 2017-09-18 | 2017-12-15 | 中南大学 | 一种基于势函数的微地震事件去噪和聚类方法 |
Non-Patent Citations (1)
| Title |
|---|
| 陈学军 等: ""基于经验小波变换的振动信号分析"", 《太阳能学报》 * |
Cited By (21)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108956764B (zh) * | 2018-06-06 | 2021-01-15 | 西安理工大学 | 一种爆炸复合管结合状态定量识别方法 |
| CN108956764A (zh) * | 2018-06-06 | 2018-12-07 | 西安理工大学 | 一种爆炸复合管结合状态定量识别方法 |
| CN109270345A (zh) * | 2018-07-20 | 2019-01-25 | 国网江苏省电力有限公司常州供电分公司 | 一种电网谐波信号的检测方法 |
| CN112740262B (zh) * | 2018-09-14 | 2024-04-02 | 莱卡微系统Cms有限责任公司 | 使用局部信噪比的去卷积装置和方法 |
| CN112740262A (zh) * | 2018-09-14 | 2021-04-30 | 莱卡微系统Cms有限责任公司 | 使用局部信噪比的去卷积装置和方法 |
| CN109520611A (zh) * | 2018-11-08 | 2019-03-26 | 温州大学 | 一种地震模拟振动台工况的监测方法 |
| CN111830559A (zh) * | 2019-04-19 | 2020-10-27 | 中国石油天然气股份有限公司 | 一种地震资料拓频处理的方法及装置 |
| CN112558147B (zh) * | 2019-09-25 | 2023-10-10 | 中国石油化工股份有限公司 | 一种井中微地震数据的偏振分析方法及系统 |
| CN112558147A (zh) * | 2019-09-25 | 2021-03-26 | 中国石油化工股份有限公司 | 一种井中微地震数据的偏振分析方法及系统 |
| CN110671613A (zh) * | 2019-10-15 | 2020-01-10 | 重庆邮电大学 | 基于改进经验小波变换的流体管道泄漏信号时延估计方法 |
| CN112764108B (zh) * | 2019-11-05 | 2023-12-15 | 防灾科技学院 | 一种基于改进经验小波变换的新型地震资料噪声压制算法 |
| CN112764108A (zh) * | 2019-11-05 | 2021-05-07 | 防灾科技学院 | 一种基于改进经验小波变换的新型地震资料噪声压制算法 |
| CN111238362B (zh) * | 2020-01-20 | 2021-10-26 | 西安理工大学 | 一种深孔钻削中孔圆度误差实时监测的方法 |
| CN111238362A (zh) * | 2020-01-20 | 2020-06-05 | 西安理工大学 | 一种深孔钻削中孔圆度误差实时监测的方法 |
| CN114114406A (zh) * | 2020-08-26 | 2022-03-01 | 中国石油化工股份有限公司 | 一种储层渗透率估计方法及装置 |
| CN114114406B (zh) * | 2020-08-26 | 2023-09-01 | 中国石油化工股份有限公司 | 一种储层渗透率估计方法及装置 |
| CN112329667A (zh) * | 2020-11-11 | 2021-02-05 | 杭州电子科技大学 | 基于经验小波变换和多尺度熵的水质数据清洗方法 |
| CN112633225B (zh) * | 2020-12-31 | 2023-07-18 | 矿冶科技集团有限公司 | 矿用微震信号滤波方法 |
| CN112633225A (zh) * | 2020-12-31 | 2021-04-09 | 矿冶科技集团有限公司 | 矿用微震信号滤波方法 |
| CN113158797B (zh) * | 2021-03-16 | 2023-08-01 | 南方科技大学 | 一种微震资料去噪方法、系统、电子设备及存储介质 |
| CN113158797A (zh) * | 2021-03-16 | 2021-07-23 | 南方科技大学 | 一种微震资料去噪方法、系统、电子设备及存储介质 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108107475B (zh) | 2019-07-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108107475B (zh) | 基于经验小波变换和多阈值函数的井中微地震去噪方法 | |
| Daubechies et al. | Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool | |
| Pinnegar et al. | The S-transform with windows of arbitrary and varying shape | |
| CN110454687A (zh) | 一种基于改进vmd的管道多点泄漏定位方法 | |
| CN105116442B (zh) | 岩性油气藏弱反射地震信号的重构方法 | |
| CN112083495B (zh) | 基于变分模态分解的同步压缩小波变换提高分辨率方法 | |
| CN106842321B (zh) | 地震数据重建方法和装置 | |
| CN111505716B (zh) | 一种基于时间同步抽取广义Chirplet变换的地震时频分析方法 | |
| Bao et al. | EMD-based extraction of modulated cavitation noise | |
| CN107957566A (zh) | 基于频率选择奇异谱分析的磁共振测深信号提取方法 | |
| CN107144879A (zh) | 一种基于自适应滤波与小波变换结合的地震波降噪方法 | |
| CN105467442B (zh) | 全局优化的时变稀疏反褶积方法及装置 | |
| Leblanc et al. | Denoising of aeromagnetic data via the wavelet transform | |
| CN106053080A (zh) | 基于能量切片小波变换的滚动轴承故障特征提取方法 | |
| CN106872985B (zh) | 基于改进短时傅里叶变换的火箭弹炮口速度测量方法 | |
| Faltermeier et al. | Weighted sliding empirical mode decomposition | |
| Huerta-Lopez et al. | Time-frequency analysis of earthquake records | |
| CN106249282B (zh) | 一种适用于avaf反演的频率域地震道集生成方法 | |
| Thiebaut et al. | Observation of a fast continental shelf wave generated by a storm impacting Newfoundland using wavelet and cross-wavelet analyses | |
| CN110244360A (zh) | 基于有效频率波数域去混叠的地震数据分离方法及系统 | |
| CN107831535B (zh) | 时变相位分解与重构方法 | |
| CN105005073A (zh) | 基于局部相似度和评价反馈的时变子波提取方法 | |
| CN109959964B (zh) | 一种高铁震源地震信号的宽频背景噪声压制方法 | |
| CN112213785B (zh) | 一种基于特征增强去噪网络的地震资料沙漠噪声抑制方法 | |
| Abratkiewicz et al. | Micro-Doppler signal decomposition using second-order vertical synchrosqueezing |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190719 Termination date: 20200305 |