CN107144879B - 一种基于自适应滤波与小波变换结合的地震波降噪方法 - Google Patents

Abstract

本发明提供一种基于自适应滤波与小波变换结合的地震波降噪方法，解决如何从混合信号中提取出利于分析解释的有效信号，压制噪声，即凸显有效信号,尽量去除干扰信号，实现信号与噪声分离的目的。首先考虑使用时间延迟估计的方法快速找到有用信号，此时的有用信号中仍然包含有部分频段的噪声；接着利用小波去噪的方法实现与有用信号频谱重合的部分噪声，实现精细去噪的目的，最终得到干净有效的地震波数据，可以反映地下地质结构真实情况，对地质情况做出正确合理的检测。

Description

一种基于自适应滤波与小波变换结合的地震波降噪方法

技术领域

本发明涉及基于自适应滤波与小波变换的地震波降噪方法，属于地震波数据处理技术领域。

背景技术

路基坍塌、地面断裂等地质灾害对人民的生命安全造成较大的威胁。因此，预防路基桥基地质下层结构灾害，对已修建的工程定期进行安全性地质检测显得异常重要。同时随着国家对重大基础设施的投入增大，需要开展大量的地质勘探工作，加之我国本属于地质灾害高发区、群发区，给工程地震仪设备的研发和应用提供了广阔的市场。

利用检波器采集地震波数据的过程中，除了携带地质信息的有效波之外，还无可避免地夹杂许多类型的干扰波，其中包括检波器自身硬件设备工作过程中产生的干扰与噪声，以及周围环境带来的噪声。检波器采集数据的实验中，未产生振动信号时，也就检波器在非工作状态下也会产生噪声，而由加速度传感器输出的信号被淹没在噪声信号中，使信号的辨认产生困难，要复原携带的有用信号，必须去除信号中叠加的噪声和干扰成分。这些噪声多为高斯白噪声，其频谱分布为在整个频域范围内，必然与有用信号的频谱有重叠部分，在此不能采用基于频率的滤波处理方法进行降噪。

现有技术中，对地震波进行滤波的方法主要有：(1)反Q滤波法：模拟地震波正向传播的反过程，通过补充振幅来消除地震波在传播过程中的振幅衰减和相位畸变，但在补偿过程中，随着旅行时间和频率的增大，振幅补偿项越来越大，即在补偿过程中放大了噪声。(2)反滤波，也叫反褶积，是滤波的一种逆过程，将地震子波压缩为震源脉冲形状，形成理想的地震记录，以抵消大地滤波作用；常用的反褶积算法是脉冲反褶积，但是脉冲脉冲反褶积不能区分信号和噪声，中考虑总的频谱，可以明显提高垂直分辨率，但提高分辨率的同时会降低信噪比。(3)短时傅里叶变换，短时傅里叶变换在全时频域上的分辨率是一致的，具有分辨率单一的特点，但需要的计算量非常大，计算效率低。

发明内容

针对现有技术存在的上述不足，本发明的目的是提供一种基于自适应滤波器的地震波降噪方法，解决如何从混合信号中提取出利于分析解释的有效信号，压制噪声，即凸显有效信号,尽量去除干扰信号，实现信号与噪声分离的目的。

为实现上述目的，本发明采用如下技术手段：一种基于自适应滤波与小波变换结合的地震波降噪方法，其特征在于：包括以下步骤：

1)自适应延时估计：采用自适应滤波器对检波器的输出信号进行滤波，在准则下自动调整滤波器的抽头系数，使得经滤波后的输出信号在该准则下与相应的延迟信号相等，所述准则为最小均方误差准则，而时延估值由滤波器的抽头系数加以处理得到；

2)利用小波阈值降噪处理方法：通过将信号分解为不同时段和频带的不同系数，让有用信号集中在小波域中一些较大的小波系数，而干扰信号则分布于整个小波域，从而有效区分有用信号和噪声，处理有用信号频谱重合的部分噪声，实现精细去噪的目的。

进一步，具体步骤：

步骤1：设检波器的输出信号为x_j(k)，检波器输出信号中的高斯白噪声信号为n_j(k)，j∈{1,2}；其中，x₁(k)表示基准道检波器的输出信号，x₂(k)表示当前道检波器的输出信号；

设基准道检波器接收到的携带地质信息的有用信号为s₁(k)，当前道检波器接收到的携带地质信息的有用信号为延迟信号s₂(k)，当前道检波器与基准道检波器的时延为D，则延迟信号s₂(k)＝s₁(k-D)；

基准道检波器的输出信号可表示为：x₁(k)＝s₁(k)+n₁(k)；当前道检波器的输出信号可表示为：x₂(k)＝s₂(k)+n₂(k)；

步骤2：对上述信号进行离散化采样得到：基准道检波器的输出信号x₁(k)＝s₁(k)+n₁(k)；

当前道检波器的输出信号x₂(k)＝s₂(k)+n₂(k)；延迟信号s₂(k)＝s₁(k-D)；其中，k表示采样时刻；

根据插值定理得到有用信号s₁(k)与延迟信号s₂(k)之间的关系为：

其中i表示自适应滤波器的抽头数；

步骤3：令sin c(i-D)＝h(i-D)，其中，h(i-D)表示自适应滤波器的抽头系数，则有：

其中，p为滤波器抽头长度，从而使得自适应滤波器为可变抽头滤波器；

步骤4：采用最小均方算法对自适应滤波器的抽头系数h(i-D)进行实时调整，从而根据公式计算出时延D；

步骤5：自适应滤波器根据时延D，过滤掉检波器输出信号x_j(k)中的高斯白噪声信号n_j(k)，从而得到携带地质信息的有用信号s_j(k)；

步骤6：令s_j(k)为含噪地震信号f_j(k)，对含噪地震信号f_j(k)作离散小波分解变换得到：f_j(k)＝s′_j(k)+n′_j(k)，其中，k表示采样时刻，k＝{0,1,2,...,N-1}，N为采样点数，s′_j(k)为地震信号，n′_j(k)为噪声信号；地震信号s′_j(k)对应的小波系数为u_j,k，f_j(k)对应的小波系数为w_j,k；

步骤7：采用无偏估计计算风险阈值λ，然后根据风险阈值λ对小波系数w_j,k进行阈值滤波处理：将大于风险阈值λ的w_j,k作为估计小波系数

步骤8：判断是否成立，ε＝10^-10，若否，回到步骤7；若是，则进入步骤9；

步骤9：利用估计小波系数进行小波重构，得到估计信号将估计信号作为去噪后的地震信号。

相比现有技术，本发明具有如下有益效果：

1、本发明针对传统的滤波方法不能有效去除噪声的特点，将自适应和小波阈值去噪联合起来，不仅将采集到的地震波信息中的杂波去除了，还将检波器中与有用信号特征相近的噪声去除了，该实现方法简单，计算量小，在放大信号的同时去除了噪声，提高信噪比的同时还提高了分辨率。

2、在检波器采集数据的实验中，未产生振动信号时，也就是检波器在非工作状态下也会产生噪声，在没有包含有效波的条件下，研究只有单纯的噪声的特性。同时考虑每道检波器的波形具有时间延迟的特征，且其波形具有一定的规律。通过以上特点，首先考虑使用时间延迟估计的方法快速找到有用信号，此时的有用信号中仍然包含有部分频段的噪声；接着利用小波去噪的方法实现与有用信号频谱重合的部分噪声，实现精细去噪的目的，最终得到干净有效的地震波数据，可以映地下地质结构真实情况，对地质情况做出正确合理的检测。

附图说明

图1为自适应时延估计方法的原理框图；

图2为本发明基于自适应滤波器的地震波降噪方法流程图；

图3为小波去噪流程图；

图4为原始地震波信号波形图；

图5为单一去噪地震波信号波形图；

图6为混合去噪地震波信号波形图。

具体实施方式

下面结合附图和具体实施方式对本发明的技术方案作进一步的说明。

一种基于自适应滤波与小波变换结合的地震波降噪方法，包括以下步骤：

1.自适应时延估计方法：

自适应时延估计方法采用自适应滤波器对检波器的输出信号进行滤波，并在某种准则下自动调整滤波器的抽头系数，使得经滤波后的输出信号在该准则(例如最小均方误差准则)下与相应的延迟信号相等，而时延估值由滤波器的抽头系数加以处理得到。该滤波器除了能较准确的的时延值外，还能够去除一部分噪声，从而达到提高信噪比的目的。自适应时延估计方法的原理框图如图1所示。

设两个检波器输出的信号的时延估计问题的数学模型为：

其中，x₁(k)表示第i(i＝1,2)个检波器输出的信号，s₁(k)表示第一个检波器输出的携带地质信息的有用信号，s₂(k)＝s₁(k-D)表示相应的延迟信号，其相对于第一个检波器的时延为D。n₁(k)表示第i(i＝1,2)个检波器输出信号中的噪声。根据上述数学模型，考虑其离散形式，根据插值定理，有用信号s₁(k)与其延迟信号s₁(k-D)的关系为：

若把sin c(i-D)看成是一个有限长横向滤波器的抽头系数，那么(2)式可以写为

其中，2P+1是滤波器抽头长度，且要求P足够长，上式的约等于号才能成立。将该滤波器以可变抽头滤波器的形式实现，并采用LMS准则对抽头系数进行自适应实时调整，那么时延估计问题最终可以得到以下的稳态解：

即维纳解。其中W^*是滤波器的抽头系数矢量，是输出信号的自相关矩阵，是输出信号与期望信号的互相关矢量。其对应的频域解为：

2.小波去噪：

小波分析具有优良的时-频分析能力，通过将信号分解为不同时段和频带的不同系数，可让有用信号集中在小波域中一些较大的小波系数，而干扰信号则分布于整个小波域，从而有效区分有用信号和噪声，成为地震信号去噪的有力工具。

小波阈值噪声压制的方法是最早提出来的小波变换压制噪声的方法，它是一种实现简单而效果较好的噪声压制方法。有效信号和随机噪声在不同的尺度上进行小波分解时有不同的奇异性随着尺度的逐渐增加，对应于有效信号的小波变换模极大值也将越来越大，而噪声的小波变换模极大值幅度及模极大值稠密度则越来越小。因此当变换尺度最小时，模极大值几乎完全由噪声所控制，而变换尺度最大时，模极大值几乎完全由信号所控制。小波变换下信号与噪声奇异性截然不同的表现特征为区分信号和噪声提供了条件。

对含噪地震信号f(k)＝s(k)+n(k)k＝0,1,2,...N-1作小波变换，得到一组小波系数w_j,k；通过对w_j,k用软硬阈值函数进行阈值处理，得出估计小波系数使得尽量小，u_j,k为对f(k)做离散小波变换后，有效信号s(k)对应的小波系数；利用进行小波重构，得到估计信号即为去噪后的地震信号。

最关键的是如何选取阈值和如何进行阈值的量化，从某种程度上说，它直接关系到信号消噪的质量。采用Rigrsure阈值，利用MATLAB可以计算出合理的阈值，该阈值即是利用Stein的无偏估计求出阈值，其具体算法如下：把信号x(n)中的每一个元素取绝对值，再由小到大排序，然后将各个元素取平方，从而得到新的信号序列s(k)＝(sort(|x|))²k＝0,1...N-1，式中，sort是MATLAB中的排序命令。若取阈值为s(k)的第k个元素的平方根，即则该阈值产生的风险为根据所得到的风险曲线，记其最小风险点所对应的值为k_min，那Rigrsure阈值定义为

由于地震信号频率的时变性，以及地震信号不同频率成分中的信噪比对分辨率具有不同贡献的性质，利用地震信号剖面相邻道的有效信号具有较好相关性这一性质，并根据有效信号和噪声在小波变换下奇异性截然不同的表现特性，采用二进小波变换的方法对原始地震信号进行分解，然后在各分量记录上逐道进行分时分频相关分析，并依此对其进行加权处理，衰减噪声。

参见图2，一种基于自适应滤波器的地震波降噪方法，具体处理方法为：

步骤1：设检波器的输出信号为x_j(k)，检波器输出信号中的高斯白噪声信号为n_j(k)，j∈{1,2}；其中，x₁(k)表示基准道检波器的输出信号，x₂(k)表示当前道检波器的输出信号；

设基准道检波器接收到的携带地质信息的有用信号为s₁(k)，当前道检波器接收到的携带地质信息的有用信号为延迟信号s₂(k)，当前道检波器与基准道检波器的时延为D，则延迟信号s₂(k)＝s₁(k-D)；

基准道检波器的输出信号可表示为：x₁(k)＝s₁(k)+n₁(k)；当前道检波器的输出信号可表示为：x₂(k)＝s₂(k)+n₂(k)；

步骤2：对上述信号进行离散化采样得到：基准道检波器的输出信号x₁(k)＝s₁(k)+n₁(k)；

当前道检波器的输出信号x₂(k)＝s₂(k)+n₂(k)；延迟信号s₂(k)＝s₁(k-D)；其中，k表示采样时刻；

根据插值定理得到有用信号s₁(k)与延迟信号s₂(k)之间的关系为：

其中i表示自适应滤波器的抽头数；

步骤3：令sin c(i-D)＝h(i-D)，其中，h(i-D)表示自适应滤波器的抽头系数，则有：

其中，p为滤波器抽头长度，从而使得自适应滤波器为可变抽头滤波器；

步骤4：采用最小均方算法对自适应滤波器的抽头系数h(i-D)进行实时调整，从而根据公式计算出时延D；

步骤5：自适应滤波器根据时延D，过滤掉检波器输出信号x_j(k)中的高斯白噪声信号n_j(t)，从而得到携带地质信息的有用信号s_j(k)；

步骤6：令s_j(k)为含噪地震信号f_j(k)，对含噪地震信号f_j(k)作离散小波分解变换得到：f_j(k)＝s′_j(k)+n′_j(k)，其中，k表示采样时刻，k＝{0,1,2,...,N-1}，N为采样点数，s′_j(k)为地震信号，n′_j(k)为噪声信号；地震信号s′_j(k)对应的小波系数为u_j,k，f_j(k)对应的小波系数为w_j,k

步骤7：采用无偏估计计算风险阈值λ，然后根据风险阈值λ对小波系数w_j,k进行阈值滤波处理：将大于风险阈值λ的w_j,k作为估计小波系数

步骤8：判断是否成立，ε＝10^-10，若否，回到步骤7；若是，则进入步骤9；

步骤9：利用估计小波系数进行小波重构，得到估计信号将估计信号作为去噪后的地震信号。

进一步，步骤4中抽头系数h(i-D)按照如下步骤进行实时调整：

步骤401：选定抽头系数h(i-D)的初始值h(0)，令h(0)＝0，则由抽头系数h(i-D)组成的抽头系数矩阵H(k)的初始值H(0)＝0；

步骤402：计算当前采样时刻k的自适应滤波器的输出y(k)，根据如下公式：

y(k)＝X₁ ^T(k)H(k)，其中，X₁(k)为自适应滤波器的输入矩阵；X₁(k)为x₁(k)的矩阵形式。

步骤403：计算误差信号g(k)，g(k)＝y(k)-x₂(k)；(x₂(k)对应步骤1)x₂(t))

步骤404：计算下一采样时刻的抽头系数矩阵：H(k+1)＝H(k)+2μg(k)X₁(k)，其中μ为自适应迭代的步长；

步骤405：判断误差信号g(k)是否满足收敛条件，若是，则输出抽头系数矩阵(见步骤4)；若否，则令k＝k+1，回到步骤402。

所述步骤7中采用如下步骤计算风险阈值λ：

步骤701：把f_j(k)中每一个元素取绝对值(f_j(k)为离散序列)，再从小到大排序，然后将各个元素取平方，从而得到新的信号序列p(k)＝(sort(|k|))²，k＝{0,1,2,...,N-1}，(k指的是p(k)中的某一个具体元素的序号，如p(1)p(2))其中，sort是指MATLAB中的排序命令；

步骤702：选取信号序列p(k)中的第k个元素的平方根作为随机阈值，那么阈值的表达式选取该阈值λ_k所产生的风险曲线为：

步骤703：选取风险曲线上的最小风险点计算风险阈值λ，按如下公式：其中σ为噪声方差，k_min为最小风险点所对应的值。

图3为小波去噪步骤。将含有噪音的信号进行小波分解，得到多个频带的小波系数，然后将高低频成分加权叠加，得到噪音压制后的信号并计算噪音压制后的信号与原信号相减，得到误差曲线。最终得到干净有效的地震波数据，可以映地下地质结构真实情况，对地质情况做出正确合理的检测。如图2所示是利用小波包多阈值去噪的流程，该流程主要由：数据获得、小波变换、分层数确定、阈值选择、阈值函数、小波重构、数据导出选择等组成。其中：

(1)数据获得：在地面上布置多个检波器，而地震波的原始数据就是通过这些高精度的检波器获得。

(2)离散小波包分解：小波变换是近十几年信号处理领域研究的一个热点，许多学者将小波在理论上的研究成果应用到诸如图像压缩、特征提取、信号滤波和数据融合等方面。而且小波变换的应用领域还在不断地发展当中，而离散小波包分解是在Mallat算法的基础上进行更加精细的划分，即对高频系数也进行划分。

(3)分层数确定：对于不同点数的信号，可以分解的层数也不同，例如对于点数N＝1024，则进行分层时，最多可以分解为10层。

(4)阈值选择：小波阈值去噪最关键的是如何选取阈值和如何进行阈值的量化，从某种程度上说，它直接关系到信号消噪的质量。采用Rigrsure阈值，利用MATLAB可以计算出合理的阈值，该阈值即是利用Stein的无偏估计求出阈值，其具体算法如下：把信号f(k)中的每一个元素取绝对值，再由小到大排序，然后将各个元素取平方，从而得到新的信号序列p(k)＝(sort(|k|))² k＝0,1...N-1，式中，sort是MATLAB中的排序命令。若取阈值为p(k)的第k个元素的平方根，即则该阈值产生的风险为根据所得到的风险曲线，记其最小风险点所对应的值为k_min，那Rigrsure阈值定义为其中λ为预置阈值，σ为噪声方差。

(5)小波函数选择：对于不同的小波函数，所表现出的正交性、紧支撑性、正则性和对称性也会不同，而这些性质的不同会影响降噪的效果，因此进行小波阈值降噪时要选择合适的小波函数。

(6)阈值函数选择：阈值函数可以选择硬阈值函数与软阈值函数。

(7)小波重构：用小波分层、选择阈值进行降噪后需要将这些小波系数进行重构，得到去完随机噪声之后的信号。

图4为经检波器采集并输出的原始地震波信号，从图4中可以明显的看出，各道集中含有较多的噪声，并且携带地质信息的有用信号已经被各种类型的噪声给淹没。

针对上述情况，需要对地震波信号进行去噪，从而提高信号的信噪比，对地震波信号进行滤波处理到不同阶段的波形效果图。首先考虑使用自适应时延估计的方法快速找到有用信号，加快去噪处理的速度。该方法处理后，与图4的原始地震波信号进行对比分析可以看出，信噪比明显地提高了。

自适应时延估计方法是用于快速找到有用信号，加快去噪处理的速度，但这一步骤并没有完全去除道集中的噪声，为此需要执行进一步的去噪工作。本专利采用小波去噪的方法去除剩余噪声，实现精细去噪的目的，最终得到干净有效的地震波数据。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (3)

1.一种基于自适应滤波与小波变换结合的地震波降噪方法，其特征在于：包括以下步骤：

1)自适应延时估计：采用自适应滤波器对检波器的输出信号进行滤波，在准则下自动调整滤波器的抽头系数，使得经滤波后的输出信号在该准则下与相应的延迟信号相等，所述准则为最小均方误差准则，而时延估值由滤波器的抽头系数加以处理得到；

2)利用小波阈值降噪处理方法：通过将信号分解为不同时段和频带的不同系数，让有用信号集中在小波域中一些较大的小波系数，而干扰信号则分布于整个小波域，从而有效区分有用信号和噪声，处理有用信号频谱重合的部分噪声，实现精细去噪的目的；

具体步骤为：

步骤1：设检波器的输出信号为x_j(k)，检波器输出信号中的高斯白噪声信号为n_j(k)，j∈{1,2}；其中，x₁(k)表示基准道检波器的输出信号，x₂(k)表示当前道检波器的输出信号；

设基准道检波器接收到的携带地质信息的有用信号为s₁(k)，当前道检波器接收到的携带地质信息的有用信号为延迟信号s₂(k)，当前道检波器与基准道检波器的时延为D，则延迟信号s₂(k)＝s₁(k-D)；

基准道检波器的输出信号可表示为：x₁(k)＝s₁(k)+n₁(k)；当前道检波器的输出信号可表示为：x₂(k)＝s₂(k)+n₂(k)；

步骤2：对上述信号进行离散化采样得到：基准道检波器的输出信号x₁(k)＝s₁(k)+n₁(k)；

当前道检波器的输出信号x₂(k)＝s₂(k)+n₂(k)；延迟信号s₂(k)＝s₁(k-D)；其中，k表示采样时刻；

根据插值定理得到有用信号s₁(k)与延迟信号s₂(k)之间的关系为：

其中i表示自适应滤波器的抽头数；

步骤3：令sinc(i-D)＝h(i-D)，其中，h(i-D)表示自适应滤波器的抽头系数，则有：

其中，p为滤波器抽头长度，从而使得自适应滤波器为可变抽头滤波器；

步骤4：采用最小均方算法对自适应滤波器的抽头系数h(i-D)进行实时调整，从而根据公式计算出时延D；

步骤5：自适应滤波器根据时延D，过滤掉检波器输出信号x_j(k)中的高斯白噪声信号n_j(k)，从而得到携带地质信息的有用信号s_j(k)；

步骤6：令s_j(k)为含噪地震信号f_j(k)，对含噪地震信号f_j(k)作离散小波分解变换得到：f_j(k)＝s′_j(k)+n′_j(k)，其中，k表示采样时刻，k＝{0,1,2,...,N-1}，N为采样点数，s′_j(k)为地震信号，n′_j(k)为噪声信号；地震信号s′_j(k)对应的小波系数为u_j,k，f_j(k)对应的小波系数为w_j,k；

步骤7：采用无偏估计计算风险阈值λ，然后根据风险阈值λ对小波系数w_j,k进行阈值滤波处理：将大于风险阈值λ的w_j,k作为估计小波系数

步骤8：判断是否成立，ε＝10^-10，若否，回到步骤7；若是，则进入步骤9；

步骤9：利用估计小波系数进行小波重构，得到估计信号将估计信号作为去噪后的地震信号。

2.根据权利要求1所述的一种基于自适应滤波与小波变换结合的地震波降噪方法，其特征在于：步骤4中抽头系数h(i-D)按照如下步骤进行实时调整：

步骤401：选定抽头系数h(i-D)的初始值h(0)，令h(0)＝0，则由抽头系数h(i-D)组成的抽头系数矩阵H(k)的初始值H(0)＝0；

步骤402：计算当前采样时刻k的自适应滤波器的输出y(k)，根据如下公式：

y(k)＝X₁ ^T(k)H(k)，其中，X₁(k)为自适应滤波器的输入矩阵；

步骤403：计算误差信号g(k)，g(k)＝y(k)-x₂(k)；步骤404：计算下一采样时刻的抽头系数矩阵：H(k+1)＝H(k)+2μg(k)X₁(k)，其中μ为自适应迭代的步长；

步骤405：判断误差信号g(k)是否满足收敛条件，若是，则输出抽头系数矩阵；若否，则令k＝k+1，回到步骤402。

3.根据权利要求1所述的一种基于自适应滤波与小波变换结合的地震波降噪方法，其特征在于：步骤7中采用如下步骤计算风险阈值λ：

步骤701：把f_j(k)中每一个元素取绝对值，再从小到大排序，然后将各个元素取平方，从而得到新的信号序列p(k)＝(sort(|k|))²，k＝{0_,1,2,...,N-1}，其中，sort是指MATLAB中的排序命令；

步骤702：选取信号序列p(k)中的第k个元素的平方根作为随机阈值，那么阈值的表达式k∈{0,1,...,N-1}；

选取该阈值λ_k所产生的风险曲线为：

步骤703：选取风险曲线上的最小风险点计算风险阈值λ，按如下公式：其中σ为噪声方差，k_min为最小风险点所对应的值。