CN108845352B - 基于vmd近似熵与多层感知机的沙漠地震信号去噪方法 - Google Patents
基于vmd近似熵与多层感知机的沙漠地震信号去噪方法 Download PDFInfo
- Publication number
- CN108845352B CN108845352B CN201810682626.3A CN201810682626A CN108845352B CN 108845352 B CN108845352 B CN 108845352B CN 201810682626 A CN201810682626 A CN 201810682626A CN 108845352 B CN108845352 B CN 108845352B
- Authority
- CN
- China
- Prior art keywords
- noise
- component
- signal
- desert
- mlp
- 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
- 238000000034 method Methods 0.000 title claims abstract description 56
- 230000008447 perception Effects 0.000 title claims abstract description 42
- 230000008901 benefit Effects 0.000 claims abstract description 68
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 37
- 230000008569 process Effects 0.000 claims description 18
- 230000006870 function Effects 0.000 claims description 17
- 238000005457 optimization Methods 0.000 claims description 17
- 230000000694 effects Effects 0.000 claims description 14
- 238000001228 spectrum Methods 0.000 claims description 13
- 238000012545 processing Methods 0.000 claims description 12
- 238000013528 artificial neural network Methods 0.000 claims description 11
- 238000004422 calculation algorithm Methods 0.000 claims description 11
- 238000012549 training Methods 0.000 claims description 9
- 230000005540 biological transmission Effects 0.000 claims description 6
- 238000000605 extraction Methods 0.000 claims description 6
- 230000009466 transformation Effects 0.000 claims description 5
- 230000015572 biosynthetic process Effects 0.000 claims description 4
- 230000001537 neural effect Effects 0.000 claims description 4
- 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 description 3
- 238000007664 blowing Methods 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 claims description 3
- 238000013329 compounding Methods 0.000 claims description 3
- 238000001914 filtration Methods 0.000 claims description 3
- 238000002156 mixing Methods 0.000 claims description 3
- 238000004800 variational method Methods 0.000 claims description 3
- 235000013399 edible fruits Nutrition 0.000 claims description 2
- 230000004069 differentiation Effects 0.000 claims 1
- 230000007547 defect Effects 0.000 description 3
- 239000000284 extract Substances 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 238000000926 separation method Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 1
- 238000005452 bending Methods 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 239000004615 ingredient Substances 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 238000003786 synthesis reaction 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/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/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
本发明涉及一种基于VMD近似熵与多层感知机的沙漠地震信号去噪方法，属于地球物理技术领域。将二维的沙漠地震记录进行变分模态分解，得到一系列的本征模态分量，计算各本征模态分量的近似熵，将所有的本征模态分量分为有效信号优势分量与噪声优势分量，通过有效信号相关性构造一个特征量，将特征量输入多层感知机进行分类，将多层感知机分类器判定为有效信号部分保留噪声部分去除，噪声优势分量则通过自相关系数与多层感知机结合进行去噪，最后通过重构实现沙漠地震信号去噪。本发明实现了强噪低信噪比条件下的噪声去除，具有速度快、准确性高及较强的抗干扰能力，可以实现低频低信噪比下的沙漠地震信号去噪。
Description
技术领域
本发明属于地球物理技术领域，尤其是指沙漠地震信号去噪方法。
背景技术
地震勘探是油气勘探的一种主要方式，是人们获得地质信息的主要手段。地震勘探利用地震波在不同地质具有不同的反射、散射和透射等来对地下的情况进行勘探。近年来，油气资源的需求量不断增长，易开采、易勘探的资源不断减少，从而对地震资料的质量提出了更高的要求。沙漠地区是目前主要的地震勘探地点之一，所以沙漠地区噪声压制是需要我们研究的一个重要课题。我国沙漠地带地震勘探随机噪声的复杂性、随机性给地震资料的处理带来很大困难。沙漠地区随机噪声主要由风吹地表所致的噪声和远场近场人文噪声组成。其中，风吹地表所致的噪声是由风对地表作用力使地表发生形变而产生弹性波引起的。近场人文噪声是由人类活动引起的噪声，远场人文噪声是由不同方向、不同大小、不同位置的噪声源共同激发的。沙漠地带随机噪声最主要的成分是由风吹地表所致的噪声，该噪声主要集中在低频段，所以沙漠地带随机噪声具有低频、弱相似性等性质，且二维的沙漠地震记录具有低信噪比的特点。
传统方法通常通过频率、振幅、波形变化等来区分有效信号与噪声，但在沙漠地震勘探中，有效信号与噪声都是低频信号，且它们的振幅都较小，所以我们需要寻找其他有效信号与噪声的区别来区分信号。fx反褶积滤波器是应用最为广泛的叠后消除随机干扰的方法之一，对于一个频率一定，道间距一定的二维的沙漠地震记录，在频率域可以用前一道沙漠地震勘探的记录来近似表示后一道沙漠地震勘探的记录，在频率域计算一个预测算子进行褶积运算，达到消除随机干扰提高信噪比的目的，但是这种方法容易产生假轴，并且对直轴较为有效但对弯曲轴效果不佳；近年来，学者将小波变换与熵值相结合进行去噪，主要应用于心电信号、故障检测等方面。小波熵去噪主要分为两类：将二维的沙漠地震记录进行小波变换，选择含有有效信号的系数，依据小波熵自适应计算全局阈值压制噪声；另一类是熵值判断做预处理，首先对分解后的小波系数进行熵值计算，根据熵值设置阈值，将小波系数分为含有有效信号与不含有有效信号的部分，再将含有有效信号的部分进行硬阈值处理达到噪声压制的目的，但是小波变换存在无法处理高维数据的缺陷，所以学者们开始使用其他分频的方法对二维的沙漠地震记录进行处理，变分模态分解VMD是一种新型分频方法，它利用迭代搜索变分模型最优解的方式来提取每个VIMF分量，通过变分模态分解，可以将一个复杂的二维沙漠地震记录分解为K个VIMF分量信号，K为任意正整数，要得到好的分解效果,K的选取十分重要，K值若选取合适，可以避免待测信号由于模态混叠现象而给分解带来的不利干扰，变分模态分解改善了经典模态分解对噪声敏感，容易产生严重混叠现象的缺点，也解决了传统小波变换在处理高维数据的缺陷，具有分解能力强、抗噪性能好和处理速度快等特点，但利用变分模态分解与普通硬阈值结合的方法处理沙漠噪声，在有效信号与噪声频率相近情况下无法有效实现信噪分离。
发明内容
本发明提供一种基于VMD近似熵与多层感知机的沙漠地震信号去噪方法，以解决存在的极低信噪比下的沙漠噪声去噪不准确，沙漠噪声在低频与有效信号重叠、导致信噪分离困难的问题。
本发明采取的技术方案是，包括下列步骤：
1)二维的沙漠地震记录获取：
沙漠地震勘探利用与地质构造走向垂直的方向安排测线，在一条测线上等距放置N个检波器，可布置m条测线，m、N为任意正整数，通过在中间点放炮，测线接收地震波信号的形式，将一条测线上的检波器接收到的地震波信号，用放大器处理后，即可得到一道沙漠地震勘探的记录，m条测线得到的m道记录即可组成一幅二维的沙漠地震记录Z；
2)二维的沙漠地震记录的分频处理：
利用变分模态分解VMD的可变尺度分解对二维的沙漠地震记录进行处理，变分模态分解的过程是将一个复杂的二维的沙漠地震记录分解为K个VIMF分量，K为任意正整数，变分模态分解的基本过程如下：
(1)希尔伯特变换，对每一个VIMF分量进行希尔伯特变换，得到单边频谱如下：
式中，uk(t)是第k个VIMF分量，k为任意正整数，uk(t)可看作谐波信号：其中，Ak(t)是uk(t)的瞬时幅值，且Ak(t)≥0；为uk(t)的相位，且δ(t)为冲激函数；t是时间；*为卷积运算；
(2)频率混合，给上述单边频谱信号混合一个预先估计的中心频率混合后，每一个uk(t)的频谱就将移动到它相应的基频带上，即：
(3)带宽估计，计算公式(2)所得到的解调信号的梯度的二范数的平方：
其中，为函数对时间t的偏导数；
(4)建立最优化模型，通过上述步骤可得到变分约束模型如下：
式中，K为VIMF分量总数；{uk(t)}＝{u1(t),u2(t),…,uk(t)}为模态的集合，{ωk(t)}＝{ω1(t),ω2(t),…,ωk(t)}表示{uk(t)}的频率中心，x(t)为VIMF分量{uk(t)}的总和；
(5)转化为非约束问题，引入二次惩罚因子α与拉格朗日乘子λ(t)，构成表达式：
{uk(t)}、{ωk(t)}和λ(t)初始化值均为0，通过交替方向乘子算法迭代求取式(4)的最优解，也可以得到模态分量{uk(t)}与相应的中心频率{ωk(t)}如下所示：
式中：为当前余项的维纳滤波；为模态功率谱的中心；为{uk(t)}的傅里叶变换，n为任意正整数，ω为频率，i为任意正整数；
变分模态分解对二维的沙漠地震记录分解的具体过程如下：设置n＝n+1，k＝k+1，根据公式(6)(7)更新{uk(t)}和{ωk(t)}，直到k＝K结束循环；根据式(8)更新λ(t)，即：
式中γ代表噪声容限参数，一般取γ＝0，一直循环上述操作，ε为误差，取任意正整数，直到满足公式(9)：
即得到由变分模态分解后的一系列优化的本征模态分量；
3)根据二维的沙漠地震记录的特性，设计变分模态近似熵特征提取：
(1)由步骤2)得到的一系列优化的本征模态分量，计算所有这些优化的本征模态分量的近似熵系数，由于近似熵是度量信号随机性的量，所以近似熵较大的分量被判定为噪声优势分量，相对地近似熵较小的分量被判定为有效信号优势分量，可以将信噪进行一个初步分离：
式中，E[]表示有效信号优势分量，F[]表示噪声优势分量，IMF1表示一个优化的本征模态分量，IMF2也表示一个优化的本征模态分量，J1表示IMF1的近似熵，J2表示IMF2的近似熵，y表示阈值；
由公式(10)得到有效信号优势分量E[]与噪声优势分量F[]，为了得到更好的去噪效果，分别处理有效信号优势分量与噪声优势分量，有效信号优势分量根据信号互相关性构造更精确的特征量v；虽然噪声优势分量中大部分是噪声，但为了得到更好的保留有效信号的效果，同样处理噪声优势分量，将噪声优势分量中无法被近似熵J1，J2区分的有效信号提取出来，所以噪声优势分量根据信号自相关性构造自相关系数ESIMF2作为噪声优势分量的信号特征量同样进行分类去噪：
(2)对于有效信号优势分量E[]来说，选取E[]中任一优化的本征模态分量IMF1，B是以IMF1为中心的一个s×s的领域窗口，s为任意正整数，令V为领域窗口B内优化的本征模态分量的平方和，即：
V＝∑B2 (11)
V即为窗口B内中心点与周围点的互相关性系数，由于噪声具有随机性而有效信号具有规律性，所以有效信号点与周围点之间的互相关性强于噪声点，则含有有效信号的窗口求出的V大于只含有噪声的窗口所求出的V，由于在二维的沙漠地震记录中，所有的B数值都是非常小的，所以对于有效信号的V与噪声的V之间的差异也是较小的，通过下列公式(12)，可以增加有效信号与噪声之间的差异，构成更精确的特征量，即定义一个更精确的特征量v：
将构造好的特征量v输入多层感知机分类器中进行分类，通过多层感知机训练数据，将输出与期望值做对比得到误差(p-q)，再根据误差(p-q)构造损失函数Hp，并且通过优化方法更新神经网络的权重Wcd，最后通过不断的训练与优化，得到训练好的神经网络分类器，再利用神经网络分类器对需要测试的二维的沙漠地震记录进行分类，将判定为有效信号部分v2保留，噪声部分v1去除,如下式：
最后得到一个去除噪声后的本征模态分量
(3)对于噪声优势分量F[]，计算自相关系数作为噪声优势分量的信号特征量同样输入多层感知机分类器中进行分类去噪，自相关系数定义如下：
有效信号自相关性是大于噪声的自相关性的，所以在噪声优势分量中，计算出的自相关性系数大的部分将判定为有效信号，而自相关性系数小的部分则判定为噪声，将噪声优势分量的信号特征量输入多层感知机中训练分类，如下式：
最后达到去除噪声的目的，得到一个去除噪声后的本征模态分量
(4)最后通过步骤3)中的(1)(2)处理后得到的去除噪声后的本征模态分量与进行重构如下：
得到去噪后的二维的沙漠地震记录
本发明所述步骤3)中(1)的变分模态近似熵的分量判定过程如下：
近似熵的基本算法：确定维数g，用时间序列信号e构造一组g维向量：
e(a)＝{ua,ua+1,…,ua+g-1} (17)
将以上构造向量中任意两向量之间的距离定义为：
计算以上构造向量中任意两向量之间的关联程度Ca(g,r)：
式中，he(·)为Heaciside函数；r为相似容限；h为[0,g-1]区域内任意正整数；l为任意正整数；计算矢量序列{e(l)}的平均自相关程度Φ(g,r)：
即可求解近似熵J1：
J1＝Φ(g,r)-Φ(g+1,r) (21)
J2求法与J1相同；
本发明所述步骤3)中(2)的基于相关性的特征提取与多层感知机的使用包括：
多层感知机是由一个输入层、多层隐藏层和一个输出层组成的，它重复两个过程：输入信号的前向传输和误差反向传播，多层感知机的模型公式如下：
隐藏层信号：Oc＝f(∑Wcd×Xc-θd)；
输出层信号：Ye＝f(∑Tcd×Oc-θd)；
在此公式中，Xc是输入信号，Ye是输出信号，Wcd是权重系数，θd表示神经单元d的阈值，误差函数Ep定义为：
在此公式中，p是输出，q是期望值。
本发明利用有效信号的特征：相比于噪声，有效信号随机性更小，所以使用近似熵对变分模态分解后的一系列优化的本征模态分量进行分量判定；并且有效信号的互相关性强于噪声，所以构造了一个基于互相关性的更精确的特征量，将该特征量输入多层感知机进行分类，通过分类实现去噪；同样有效信号的自相关性也是强于噪声的，所以对于噪声优势分量，采用自相关性系数作为信号特征量输入多层感知机分类去噪。
本发明利用变分模态分解VMD的特点将输入的二维的沙漠地震记录先进行分频处理，通过噪声随机性强于有效信号，将一系列优化的本征模态分量分开，分开后分别处理有效信号优势分量与噪声优势分量。有效信号优势分量通过基于互相关性的特征量输入多层感知机分类，对分类结果进行保留与舍去；噪声优势分量通过自相关系数与多层感知机结合保留有效信号压制噪声。通过去除噪声后的本征模态分量进行重构达到沙漠噪声去噪的目的。
本发明使用多层感知机代替传统硬阈值来处理二维的沙漠地震记录，多层感知机是一种前向网络，由感知机推广而来，最主要的特点是有多个神经元层，因此也可以称为深度神经网络，它可以捕捉隐藏在任何物理过程背后的复杂底层机制，是解决分类问题的有力工具。本发明用多层感知机对有效信号与噪声进行分类。根据分类的结果我们可以去除噪声，保留有效信号，达到去噪的效果。多层感知机是由一个输入层、多层隐藏层和一个输出层组成，没有规定固定的隐层数量，因此可以根据各自的需求选择合适的隐层层数。它重复两个过程：输入信号的前向传输和误差反向传播。在前向传输中，输入信号通过输入层向前传播，并逐层处理直到输出层。然后将输出与期望值进行比较，将期望值与输出结果之间的差值称为误差；在误差反向传播中，利用误差来构造损失函数，并且通过优化方法更新神经网络的链路权重，以最小化损失函数。然后停止训练。由于存在多层的网络结构，因此无法直接对中间的隐层利用损失来进行参数更新，但本发明通过损失从顶层到底层的反向传播来实现参数估计。此时，经过训练的神经网络可以输出具有最小误差的信息，因此利用多层感知机进行分类去噪解决了传统硬阈值在信噪频带重叠部分无法分离的问题。
本发明比现有方法检测的更为准确，能够实现极低信噪比下的沙漠噪声去噪。利用传统的硬阈值在信噪重叠的频带无法实现有效分离，本发明解决了沙漠噪声在低频与有效信号重叠的困难，通过多层感知机分类实现信噪分离。二维的沙漠地震记录去噪是地震勘探十分重要的一步，为精确探明后续地质结构处理提供了重要依据，有利于提供更好的资源开发设计。
本发明的优点：与传统使用变分模态分量系数、分解后分量的振幅、频率等作为信号特征量相比，本发明提出了利用近似熵对信噪进行分离，它更能体现出信号与噪声的差异。对于后续多层感知机分类来说，相比于传统特征量，本发明构造基于互相关性的特征量更能表明二维的沙漠地震记录的特征，针对二维的沙漠地震记录幅值小，互相关性强的特点构造出合适的特征量作为多层感知机的输入；同时利用多层感知机进行分类达到去噪目的，避免了传统阈值的设置。不会存在阈值选择过小带来的将噪声误判为有效信号、阈值选择过大带来的漏判有效信号的情况。也解决了传统阈值在频带重合处无法实现信噪分离的缺陷。与传统方法相比，本发明的方法使用于低频低信噪比的情况，具有更高的准确率，能够更有效的保留有效信号、去除噪声。
附图说明
图1(a)是模拟的两个主频为30Hz、采样频率为512Hz、振幅分别为0.5、0.2的雷克子波、一共40道的沙漠地震记录；
图1(b)是向纯净雷克子波中加入高斯白噪声，合成信噪比为-9.0873的含噪二维的沙漠地震记录；从图中看到信号被大量噪声淹没，并且沙漠噪声具有弱相似性，与有效信号也十分相似，人工识别有效信号是十分困难；
图1(c)是本发明去噪的效果图；
图1(d)是变分模态硬阈值去噪的效果图；
图1(e)是fx反褶积滤波器去噪的效果图；
图2(a)是将图1(b)中一道沙漠地震勘探的记录经过变分模态分解后的图像；
图2(b)是将图1(b)中一道沙漠地震勘探的记录经过变分模态分解后的频谱图；
图3是单道信号利用本发明方法、变分模态硬阈值、fx反褶积滤波器去噪与原始信号的对比图。
具体实施方式
包括以下步骤：
1)二维的沙漠地震记录获取：
沙漠地震勘探利用与地质构造走向垂直的方向安排测线，在一条测线上等距放置N个检波器，可布置m条测线，m、N为任意正整数，通过在中间点放炮，测线接收地震波信号的形式，将一条测线上的检波器接收到的地震波信号，用放大器处理后，即可得到一道沙漠地震勘探的记录，m条测线得到的m道记录即可组成一幅二维的沙漠地震记录Z；
2)二维的沙漠地震记录的分频处理：
利用变分模态分解VMD的可变尺度分解对二维的沙漠地震记录进行处理，变分模态分解的过程是将一个复杂的二维的沙漠地震记录分解为K个VIMF分量，K为任意正整数，变分模态分解的基本过程如下：
(1)希尔伯特变换，对每一个VIMF分量进行希尔伯特变换，得到单边频谱如下：
式中，uk(t)是第k个VIMF分量，k为任意正整数，uk(t)可看作谐波信号：其中，Ak(t)是uk(t)的瞬时幅值，且Ak(t)≥0；为uk(t)的相位，且δ(t)为冲激函数；t是时间；*为卷积运算；
(2)频率混合，给上述单边频谱信号混合一个预先估计的中心频率混合后，每一个uk(t)的频谱就将移动到它相应的基频带上，即：
(3)带宽估计，计算公式(2)所得到的解调信号的梯度的二范数的平方：
其中，为函数对时间t的偏导数；
(4)建立最优化模型，通过上述步骤可得到变分约束模型如下：
式中，K为VIMF分量总数；{uk(t)}＝{u1(t),u2(t),…,uk(t)}为模态的集合，{ωk(t)}＝{ω1(t),ω2(t),…,ωk(t)}表示{uk(t)}的频率中心，x(t)为VIMF分量{uk(t)}的总和；
(5)转化为非约束问题，引入二次惩罚因子α与拉格朗日乘子λ(t)，构成表达式：
{uk(t)}、{ωk(t)}和λ(t)初始化值均为0，通过交替方向乘子算法迭代求取式(4)的最优解，也可以得到模态分量{uk(t)}与相应的中心频率{ωk(t)}如下所示：
式中：为当前余项的维纳滤波；为模态功率谱的中心；为{uk(t)}的傅里叶变换，n为任意正整数，ω为频率，i为任意正整数；
变分模态分解对二维的沙漠地震记录分解的具体过程如下：设置n＝n+1，k＝k+1，根据公式(6)(7)更新{uk(t)}和{ωk(t)}，直到k＝K结束循环；根据式(8)更新λ(t)，即：
式中γ代表噪声容限参数，一般取γ＝0，一直循环上述操作，ε为误差，可取任意正整数，直到满足公式(9)：
即可得到由变分模态分解后的一系列优化的本征模态分量；
3)根据二维的沙漠地震记录的特性，设计变分模态近似熵特征提取：
(1)由步骤2)得到的一系列优化的本征模态分量，计算所有这些优化的本征模态分量的近似熵系数，由于近似熵是度量信号随机性的量，所以近似熵较大的分量被判定为噪声优势分量，相对地近似熵较小的分量被判定为有效信号优势分量，可以将信噪进行一个初步分离：
式中，E[]表示有效信号优势分量，F[]表示噪声优势分量，IMF1表示一个优化的本征模态分量，IMF2也表示一个优化的本征模态分量，J1表示IMF1的近似熵，J2表示IMF2的近似熵，y表示阈值；
由公式(10)得到有效信号优势分量E[]与噪声优势分量F[]，为了得到更好的去噪效果，分别处理有效信号优势分量与噪声优势分量，有效信号优势分量根据信号互相关性构造更精确的特征量v；虽然噪声优势分量中大部分是噪声，但为了得到更好的保留有效信号的效果，同样处理噪声优势分量，将噪声优势分量中无法被近似熵J1，J2区分的有效信号提取出来，所以噪声优势分量根据信号自相关性构造自相关系数作为噪声优势分量的信号特征量同样进行分类去噪：
(2)对于有效信号优势分量E[]来说，选取E[]中任一优化的本征模态分量IMF1，B是以IMF1为中心的一个s×s的领域窗口，s为任意正整数，令V为领域窗口B内优化的本征模态分量的平方和，即：
V＝∑B2 (11)
V即为窗口B内中心点与周围点的互相关性系数，由于噪声具有随机性而有效信号具有规律性，所以有效信号点与周围点之间的互相关性强于噪声点，则含有有效信号的窗口求出的V大于只含有噪声的窗口所求出的V，由于在二维的沙漠地震记录中，所有的B数值都是非常小的，所以对于有效信号的V与噪声的V之间的差异也是较小的，通过下列公式(12)，可以增加有效信号与噪声之间的差异，构成更精确的特征量，即定义一个更精确的特征量v：
将构造好的特征量v输入多层感知机分类器中进行分类，通过多层感知机训练数据，将输出与期望值做对比得到误差(p-q)，再根据误差(p-q)构造损失函数Hp，并且通过优化方法更新神经网络的权重Wcd，最后通过不断的训练与优化，得到训练好的神经网络分类器，再利用神经网络分类器对需要测试的二维的沙漠地震记录进行分类，将判定为有效信号部分v2保留，噪声部分v1去除,如下式：
最后得到一个去除噪声后的本征模态分量
(3)对于噪声优势分量F[]，计算自相关系数作为噪声优势分量的信号特征量同样输入多层感知机分类器中进行分类去噪，自相关系数定义如下：
有效信号自相关性是大于噪声的自相关性的，所以在噪声优势分量中，计算出的自相关性系数大的部分将判定为有效信号，而自相关性系数小的部分则判定为噪声，将噪声优势分量的信号特征量输入多层感知机中训练分类，如下式：
最后达到去除噪声的目的，得到一个去除噪声后的本征模态分量
(4)最后通过步骤3)中的(1)(2)处理后得到的去除噪声后的本征模态分量与进行重构如下：
得到去噪后的二维的沙漠地震记录
本发明所述步骤3)中(1)的变分模态近似熵的分量判定过程如下：
近似熵的基本算法：确定维数g，用时间序列信号e构造一组g维向量：
e(a)＝{ua,ua+1,…,ua+g-1} (17)
将以上构造向量中任意两向量之间的距离定义为：
计算以上构造向量中任意两向量之间的关联程度Ca(g,r)：
式中，he(·)为Heaciside函数；r为相似容限；h为[0,g-1]区域内任意正整数；l为任意正整数；计算矢量序列{e(l)}的平均自相关程度Φ(g,r)：
即可求解近似熵J1：
J1＝Φ(g,r)-Φ(g+1,r) (21)
J2求法与J1相同；
本发明所述步骤3)中(2)的基于相关性的特征提取与多层感知机的使用包括：
多层感知机是由一个输入层、多层隐藏层和一个输出层组成的，它重复两个过程：输入信号的前向传输和误差反向传播，多层感知机的模型公式如下：
隐藏层信号：Oc＝f(∑Wcd×Xc-θd)；
输出层信号：Ye＝f(∑Tcd×Oc-θd)；
在此公式中，Xc是输入信号，Ye是输出信号，Wcd是权重系数，θd表示神经单元d的阈值，误差函数Ep定义为：
在此公式中，p是输出，q是期望值。
下边通过模拟例来进一步说明本发明的效果。
合成记录
模拟40道二维的沙漠地震记录，主频为30Hz，如图1(a)所示。每道采样频率为512Hz。幅度分别为0.5、0.2的雷克子波，图1(b)为向纯净雷克子波中加入高斯白噪声，合成信噪比为-9.0873dB的含噪数据。在变分模态分解中，k选取为6。多层感知机选择4个隐藏层。
为了说明本发明方法的有效性与优越性，将本发明方法与变分模态硬阈值算法和fx反褶积滤波器的方法进行对比，去噪结果如图1(c)(d)(e)所示。从图1(c)(d)可以看到，使用变分模态硬阈值算法随机噪声去除不够彻底，大量数据依旧被噪声干扰呈弱相似弯曲状，而本发明的方法在没有有效信号的区域将噪声压制的十分彻底，将同相轴清晰连续的恢复出来；由图1(c)(e)可以看到，只使用fx反褶积滤波器的结果同样不能很好的压制沙漠噪声，并且存在着一些细小毛刺，而本发明去除噪声十分彻底，细小毛刺也得到压制，可以得到更好的去噪效果。
选取图1(b)中的一道沙漠地震勘探的记录详细说明变分模态分解与分量判定过程，将此一道沙漠地震勘探的记录通过变分模态分解后，分解K为6，α为1000，如图2(a)(b)所示，画出分解后的VIMF分量与其相对于的频谱。计算每一个VIMF分量的近似熵，通过判定，将第1 3 4 5分量判定为有效信号优势分量，2与6分量判定为噪声优势分量。接着对两种分量分别处理。
为了更好地说明本发明方法的优越性通过观察单道信号的波形和信噪比来进一步分析，从去噪结果中抽取第5道来进行波形对比，如图3所示。可以看到，本发明方法最接近原始纯净信号的波形，尤其在波峰、波谷处，本发明方法将有效信号幅度保持的很好，且在不含有有效信号的位置噪声压制的十分干净。而硬阈值和fx反褶积滤波器的方法，噪声压制不彻底，且有效信号幅值丢失较多。
接下来，通过信噪比和均方根误差计算来说明方法的有效性。若信号信噪比越高、均方误差越小，表示去噪信号越接近原始信号，效果越好。信噪比与均方根误差的公式如下。表1、2为三种不同的方法去噪前后的信噪比与均方根误差对比：
其中为纯净沙漠地震信号，为含噪沙漠地震信号，R为信号长度。
表1在不同噪声下的信噪比对比
表2在不同噪声下的均方根误差对比
从表中我们可以看出沙漠地震记录的信噪比非常低，同样可以从表中看到，本发明的方法明显优于另外两种的方法，去噪后信噪比大幅度提高、均方误差更小，对原始纯净信号进行了更为精准的估计。
Claims (3)
1.一种基于VMD近似熵与多层感知机的沙漠地震信号去噪方法，其特征在于，包括下列步骤：
1)二维的沙漠地震记录获取：
沙漠地震勘探利用与地质构造走向垂直的方向安排测线，在一条测线上等距放置N个检波器，布置m条测线，m、N为任意正整数，通过在中间点放炮，测线接收地震波信号的形式，将一条测线上的检波器接收到的地震波信号，用放大器处理后，即可得到一道沙漠地震勘探的记录，m条测线得到的m道记录即可组成一幅二维的沙漠地震记录Z；
2)二维的沙漠地震记录的分频处理：
利用变分模态分解VMD的可变尺度分解对二维的沙漠地震记录进行处理，变分模态分解的过程是将一个复杂的二维的沙漠地震记录分解为K个VIMF分量，K为任意正整数，变分模态分解的基本过程如下：
(1)希尔伯特变换，对每一个VIMF分量进行希尔伯特变换，得到单边频谱如下：
式中，uk(t)是第k个VIMF分量，k为任意正整数，uk(t)看作谐波信号：其中，Ak(t)是uk(t)的瞬时幅值，且Ak(t)≥0；为uk(t)的相位，且δ(t)为冲激函数；t是时间；*为卷积运算；
(2)频率混合，给上述单边频谱信号混合一个预先估计的中心频率混合后，每一个uk(t)的频谱就将移动到它相应的基频带上，即：
(3)带宽估计，计算公式(2)所得到的解调信号的梯度的二范数的平方：
其中，为函数对时间t的偏导数；
(4)建立最优化模型，通过上述步骤可得到变分约束模型如下：
式中，K为VIMF分量总数；{uk(t)}＝{u1(t),u2(t),…,uk(t)}为模态的集合，{ωk(t)}＝{ω1(t),ω2(t),…,ωk(t)}表示{uk(t)}的频率中心，x(t)为VIMF分量{uk(t)}的总和；
(5)转化为非约束问题，引入二次惩罚因子α与拉格朗日乘子λ(t)，构成表达式：
{uk(t)}、{ωk(t)}和λ(t)初始化值均为0，通过交替方向乘子算法迭代求取式(4)的最优解，也可以得到模态分量{uk(t)}与相应的中心频率{ωk(t)}如下所示：
式中：为当前余项的维纳滤波；为模态功率谱的中心；为{uk(t)}的傅里叶变换，n为任意正整数，ω为频率，i为任意正整数；
变分模态分解对二维的沙漠地震记录分解的具体过程如下：设置n＝n+1，k＝k+1，根据公式(6)(7)更新{uk(t)}和{ωk(t)}，直到k＝K结束循环；根据式(8)更新λ(t)，即：
式中γ代表噪声容限参数，取γ＝0，一直循环上述操作，ε为误差，取任意正整数，直到满足公式(9)：
即得到由变分模态分解后的一系列优化的本征模态分量；
3)根据二维的沙漠地震记录的特性，设计变分模态近似熵特征提取：
(1)由步骤2)得到的一系列优化的本征模态分量，计算所有这些优化的本征模态分量的近似熵系数，由于近似熵是度量信号随机性的量，所以近似熵较大的分量被判定为噪声优势分量，相对地近似熵较小的分量被判定为有效信号优势分量，将信噪进行一个初步分离：
式中，E[]表示有效信号优势分量，F[]表示噪声优势分量，IMF1表示一个优化的本征模态分量，IMF2也表示一个优化的本征模态分量，J1表示IMF1的近似熵，J2表示IMF2的近似熵，y表示阈值；
由公式(10)得到有效信号优势分量E[]与噪声优势分量F[]，为了得到更好的去噪效果，分别处理有效信号优势分量与噪声优势分量，有效信号优势分量根据信号互相关性构造更精确的特征量v；虽然噪声优势分量中大部分是噪声，但为了得到更好的保留有效信号的效果，同样处理噪声优势分量，将噪声优势分量中无法被近似熵J1，J2区分的有效信号提取出来，所以噪声优势分量根据信号自相关性构造自相关系数作为噪声优势分量的信号特征量同样进行分类去噪：
(2)对于有效信号优势分量E[]来说，选取E[]中任一优化的本征模态分量IMF1，B是以IMF1为中心的一个s×s的领域窗口，s为任意正整数，令V为领域窗口B内优化的本征模态分量的平方和，即：
V＝∑B2 (11)
V即为窗口B内中心点与周围点的互相关性系数，由于噪声具有随机性而有效信号具有规律性，所以有效信号点与周围点之间的互相关性强于噪声点，则含有有效信号的窗口求出的V大于只含有噪声的窗口所求出的V，由于在二维的沙漠地震记录中，所有的B数值都是非常小的，所以对于有效信号的V与噪声的V之间的差异也是较小的，通过下列公式(12)，可以增加有效信号与噪声之间的差异，构成更精确的特征量，即定义一个更精确的特征量v：
将构造好的特征量v输入多层感知机分类器中进行分类，通过多层感知机训练数据，将输出与期望值做对比得到误差(p-q)，再根据误差(p-q)构造损失函数Hp，并且通过优化方法更新神经网络的权重Wcd，最后通过不断的训练与优化，得到训练好的神经网络分类器，再利用神经网络分类器对需要测试的二维的沙漠地震记录进行分类，将判定为有效信号部分v2保留，噪声部分v1去除,如下式：
最后得到一个去除噪声后的本征模态分量
(3)对于噪声优势分量F[]，计算自相关系数作为噪声优势分量的信号特征量同样输入多层感知机分类器中进行分类去噪，自相关系数定义如下：
有效信号自相关性是大于噪声的自相关性的，所以在噪声优势分量中，计算出的自相关性系数大的部分将判定为有效信号，而自相关性系数小的部分则判定为噪声，将噪声优势分量的信号特征量输入多层感知机中训练分类，如下式：
最后达到去除噪声的目的，得到一个去除噪声后的本征模态分量
(4)最后通过步骤3)中的(1)(2)处理后得到的去除噪声后的本征模态分量与进行重构如下：
得到去噪后的二维的沙漠地震记录
2.根据权利要求1所述的一种基于VMD近似熵与多层感知机的沙漠地震信号去噪方法，其特征在于：所述步骤3)中(1)的变分模态近似熵的分量判定过程如下：
近似熵的基本算法：确定维数g，用时间序列信号e构造一组g维向量：
e(a)＝{ua,ua+1,…,ua+g-1} (17)
将以上构造向量中任意两向量之间的距离定义为：
计算以上构造向量中任意两向量之间的关联程度Ca(g,r)：
式中，he(·)为Heaciside函数；r为相似容限；h为[0,g-1]区域内任意正整数；l为任意正整数；计算矢量序列{e(l)}的平均自相关程度Φ(g,r)：
求解近似熵J1：
J1＝Φ(g,r)-Φ(g+1,r) (21)
J2求法与J1相同。
3.根据权利要求1所述的一种基于VMD近似熵与多层感知机的沙漠地震信号去噪方法，其特征在于：所述步骤3)中(2)的基于相关性的特征提取与多层感知机的使用包括：
多层感知机是由一个输入层、多层隐藏层和一个输出层组成的，它重复两个过程：输入信号的前向传输和误差反向传播，多层感知机的模型公式如下：
隐藏层信号：Oc＝f(∑Wcd×Xc-θd)；
输出层信号：Ye＝f(∑Tcd×Oc-θd)；
在此公式中，Xc是输入信号，Ye是输出信号，Wcd是权重系数，θd表示神经单元d的阈值，误差函数Ep定义为：
在此公式中，p是输出，q是期望值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810682626.3A CN108845352B (zh) | 2018-06-27 | 2018-06-27 | 基于vmd近似熵与多层感知机的沙漠地震信号去噪方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810682626.3A CN108845352B (zh) | 2018-06-27 | 2018-06-27 | 基于vmd近似熵与多层感知机的沙漠地震信号去噪方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108845352A CN108845352A (zh) | 2018-11-20 |
| CN108845352B true CN108845352B (zh) | 2019-09-27 |
Family
ID=64200612
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810682626.3A Expired - Fee Related CN108845352B (zh) | 2018-06-27 | 2018-06-27 | 基于vmd近似熵与多层感知机的沙漠地震信号去噪方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108845352B (zh) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110109186A (zh) * | 2019-04-18 | 2019-08-09 | 河海大学 | 一种探地雷达数据三维时频分析方法 |
| CN110208856B (zh) * | 2019-06-05 | 2020-12-25 | 吉林大学 | 一种基于流形分区2d-vmd的沙漠复杂噪声压制方法 |
| CN110208375B (zh) * | 2019-06-13 | 2021-12-14 | 石家庄铁道大学 | 一种锚杆锚固缺陷的检测方法及终端设备 |
| CN110186682B (zh) * | 2019-07-08 | 2021-03-23 | 石家庄铁道大学 | 基于分数阶变分模态分解的滚动轴承故障诊断方法 |
| CN111723701B (zh) * | 2020-06-08 | 2022-05-20 | 西安交通大学 | 一种水中目标识别方法 |
| CN113093282A (zh) * | 2021-04-18 | 2021-07-09 | 吉林大学 | 一种基于几何模态特征并行网络的沙漠数据消噪方法 |
| CN113109873B (zh) * | 2021-04-20 | 2022-11-29 | 吉林大学 | 一种基于秩残差约束的沙漠地区地震信号噪声抑制方法 |
| CN113484913B (zh) * | 2021-08-16 | 2023-06-16 | 成都理工大学 | 一种多粒度特征融合卷积神经网络的地震资料去噪方法 |
| CN114236610B (zh) * | 2021-12-24 | 2023-03-31 | 西安交通大学 | 基于深度卷积网络先验的迭代地震资料解混方法及系统 |
| CN114545504B (zh) * | 2022-02-18 | 2023-05-09 | 吉林大学 | 一种基于深度学习的沙漠地震勘探随机噪声建模方法 |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8050414B2 (en) * | 2008-10-16 | 2011-11-01 | Gas Technology Institute | Robust pipe-strike pulse detector |
| CN103543469B (zh) * | 2012-07-17 | 2016-12-21 | 中国石油化工股份有限公司 | 一种基于小波变换的小尺度阈值去噪方法 |
| CN105093317A (zh) * | 2014-05-14 | 2015-11-25 | 中国石油化工股份有限公司 | 地面阵列式微地震数据独立分量分离去噪方法 |
| CN105513048A (zh) * | 2015-11-24 | 2016-04-20 | 西安电子科技大学昆山创新研究院 | 基于子带信息熵度量的图像质量评价方法 |
| CN106896407B (zh) * | 2017-03-28 | 2018-07-13 | 吉林大学 | 一种基于近似负熵的微地震信号初至拾取方法 |
| CN108132486A (zh) * | 2017-11-28 | 2018-06-08 | 吉林大学 | 一种重磁梯度与地震数据联合界面反演的优化模拟退火法 |
-
2018
- 2018-06-27 CN CN201810682626.3A patent/CN108845352B/zh not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| CN108845352A (zh) | 2018-11-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108845352B (zh) | 基于vmd近似熵与多层感知机的沙漠地震信号去噪方法 | |
| Adler et al. | Deep learning for seismic inverse problems: Toward the acceleration of geophysical analysis workflows | |
| Dong et al. | Denoising the optical fiber seismic data by using convolutional adversarial network based on loss balance | |
| Pinnegar et al. | The S-transform with windows of arbitrary and varying shape | |
| CN108037526B (zh) | 基于全波波场vsp/rvsp地震资料的逆时偏移方法 | |
| CN107144879B (zh) | 一种基于自适应滤波与小波变换结合的地震波降噪方法 | |
| Van den Ende et al. | A self-supervised deep learning approach for blind denoising and waveform coherence enhancement in distributed acoustic sensing data | |
| CN104216014B (zh) | 一种地震信号分频处理方法 | |
| CN108415077A (zh) | 新的边缘检测低序级断层识别方法 | |
| Sun et al. | Deep learning for low-frequency extrapolation of multicomponent data in elastic FWI | |
| CN102939546A (zh) | 用于在地震处理中的局部属性匹配的系统和方法 | |
| CN110501742A (zh) | 一种采用Boosting集成学习算法对地震事件进行区分的方法 | |
| Liu et al. | Irregularly sampled seismic data reconstruction using multiscale multidirectional adaptive prediction-error filter | |
| CN109946739A (zh) | 一种基于压缩感知理论的地震剖面增强方法 | |
| CN109507726A (zh) | 时间域弹性波多参数全波形的反演方法及系统 | |
| CN110208856A (zh) | 一种基于流形分区2d-vmd的沙漠复杂噪声压制方法 | |
| CN115905805A (zh) | 一种基于全局信息判别gan的das数据多尺度噪声消减方法 | |
| CN111948711A (zh) | 一种利用深度学习法提取地震数据低频部分的方法和系统 | |
| CN108226997A (zh) | 一种基于叠前地震数据的地震相划分方法 | |
| CN109164492B (zh) | 一种提取套管井地层声波速度的方法 | |
| Liu et al. | Unsupervised deep learning for ground roll and scattered noise attenuation | |
| Wang et al. | Low-frequency extrapolation of prestack viscoacoustic seismic data based on dense convolutional network | |
| CN113740906A (zh) | 一种水下垂直缆地震波干涉成像方法及装置 | |
| Wang et al. | Intelligent shot gather reconstruction using residual learning networks | |
| Rovetta et al. | A density-based spatial clustering application for a fully automatic picking of surface wave dispersion curves |
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: 20190927 |