CN111175827A - 一种增强地震勘探信号的高性能时频域滤波方法 - Google Patents

Abstract

一种增强地震勘探信号的高性能时频域滤波方法，先对含噪的地震信号做同步挤压小波变换得到高分辨率的时频平面，再在所得的时频平面上采用Catte模型进行各向异性扩散滤波，然后对滤波后的时频分布做同步挤压小波重构，最后得到有效信号的时域估计即为最终的滤波结果，并对滤波结果进行分析；本发明充分利用了同步挤压时频分布高分辨、高聚集特性和Catte模型各向异性扩散滤波能够兼顾有效信号特征保持与随机噪声消减的特性，提出了一种增强地震勘探信号的高性能时频域滤波方法；经验证，该方法在地震勘探信号的增强以及随机噪声的压制方面具有非常优越的性能。

Description

一种增强地震勘探信号的高性能时频域滤波方法

技术领域

本发明涉及滤波技术，特别涉及一种在高分辨率时频域进行各向异性扩散滤波以增强地震勘探信号的高性能时频域滤波方法。

背景技术

为了得到同相轴信息更加清晰的地震记录，必须要对采集到的地震勘探资料进行处理。对于从易勘探开发地区获取的地震资料，其中噪声干扰较少且情况不太复杂，应用一般的滤波方法即可达到一定的处理要求，但是如要获得同相轴信息非常清晰的处理结果，仍需寻找更先进的滤波方法。随着勘探开发在深度和难度上的增加，采集到的地震勘探资料中噪声情况较为复杂且程度较为严重。其中，随机噪声是影响地震资料质量的重要干扰类型之一，它能够使有效同相轴信息模糊、截断甚至湮没其中。如果能够很好地压制这些噪声干扰，将会使有效同相轴得到大幅度的增强，从而大大提升地震勘探资料的质量。

截至今日，国内外尚没有在高分辨率时频变换域——同步挤压变换域应用各向异性扩散滤波增强地震勘探信号的方法。现有研究成果报道中，有多种时频域滤波方法(如小波滤波、S变换滤波、时频峰值滤波(time-frequency peak filtering,TFPF)等)对地震资料的处理研究，以及一些图像处理方法(如高斯滤波、各向异性扩散滤波、基于曲波、轮廓波变换的滤波等)对地震图像的处理研究。这些研究虽已取得了显著的效果，但仍存在提升的空间。传统的时频滤波方法消减地震勘探随机噪声时，会受到各时频分布自身的限制(如时窗选取问题)而使滤波效果达到一定程度就很难再提高了。使用图像处理方法时，如果将地震数据转换成一定格式的图片进行处理，则在抽取单道信号进行波形及频谱的对比分析时，将涉及到不同格式数据间的转换问题，很多时候难以解决；更多方法需要通过对地震数据补零以构成行、列中元素个数都为2的整数次幂的矩阵进行处理，处理完后再从滤波结果中截取出对应于原数据的矩阵，这就造成整个处理过程不够快捷方便。

目前，在众多的时频滤波方法中，时频峰值滤波(TFPF)具有非常好的滤波性能。该方法在高斯白噪声条件下(信噪比可低至-9dB)，能够对随时间呈线性变化的信号进行无偏估计。其实现过程是先对含噪信号进行调频，得到原信号的解析信号，然后利用维格纳——威尔分布(Wigner-Ville distribution，PWVD)或其加窗形式，叫做伪维格纳——威尔分布(pseudo Wigner-Ville distribution，PWVD)对该解析信号做时频变换，再在时频平面上进行峰值搜索，所得峰值的集合即为有效信号的估值。该方法在地震勘探信号处理中的应用已非常成熟，并且已取得了显著的效果。由于地震勘探信号具有非线性、非平稳性，在TFPF方法中宜采用PWVD，这就使得方法本身存在一定的局限性：窗长的选择需权衡有效信号的保持与随机噪声的消减。因此，该方法的滤波效果也受到了一定的限制。而已有的基于同步挤压时频变换的滤波方法主要是在时频平面上进行区域截取(如带通滤波)或划分(如阈值划分)，然后对保留部分的时频分布进行同步挤压重构得到有效信号的估值。此种方法可以有效地去除与保留部分不同频的干扰，但是对保留部分本身的增强效果比较有限，因此滤波效果仍存在提升的空间。

应用图像处理方法中的各向异性扩散滤波对地震图像的增强，也只是沿着有效同相轴的梯度方向进行滤波处理，而地震勘探资料中的同相轴繁多且方向有所差异，因此不能在特定区域内进行集中滤波，所能达到的效果也受到一定的限制。相当多的图像处理方法，诸如基于曲波变换的滤波方法对图像进行处理时需将矩阵行、列中的元素个数凑成2的整数次幂，要么是分块截取数据，要么就是补零凑数，都会使处理过程不够简洁。其实，很多图像处理方法对于一定格式的图片文件具有较好的处理效果，如果将地震数据转换成图片格式进行滤波处理，后续在对滤波结果的分析中，需抽取单道信号进行波形对比、频谱对比，那么将会存在有关数据格式转换的一些问题，大多数情况下不能转换成原本的地震数据格式进行分析，因此非常不便。

发明内容

为了克服上述现有方法的不足，本发明的目的是提供一种在高分辨率时频域进行各向异性扩散滤波的方法，即在同步挤压时频域采用Catte模型滤波来增强地震勘探信号的高性能滤波方法，更好地增强有效信号与更好地消减噪声相辅相成，以获得同相轴清晰度更高的地震资料，为后续的资料解释等工作提供良好的前提。

为了达到上述目的，本发明的技术方案为：

一种增强地震勘探信号的高性能时频域滤波方法，包括以下步骤：

步骤一、对含噪的地震信号做同步挤压小波变换，得到高分辨率时频分布：

设含噪信号为：

s(t)＝x(t)+n(t) (1)

x(t)为有效信号，n(t)是加性随机噪声；

对s(t)做连续小波变换，参照公式(2)如下：

其中，a是尺度因子，b是平移因子，ψ(t)为小波基函数，

为ψ(t)的共轭函数，根据帕塞瓦尔定理，公式(2)的频域表达式为：

其中，ξ是角频率，

是s(t)的傅里叶变换，ψ(aξ)为ψ(t)的傅里叶变换；通过该式就将时间域的褶积形式变换为频率域的乘积形式，通过对小波系数求导，可初步估计瞬时频率，即

其中，i是虚部单位；这样就把时间—尺度平面(b,a)转换到时间—频率平面(b,ω_s(a,b))，频率ω和尺度a均被离散化，因此仅在离散点a_i{a_i-a_i-1＝(Δa)_i}计算小波系数W_s(a,b)；同步挤压小波变换中频率重排的实现过程为：假设信号的采样频率为cf，小波变换时所取尺度的个数为N，那么重排后的时频谱上，离散的频率值就为ω_l＝l*cf/N,l∈[1,N]；小波系数的同步挤压值T_s(ω_l,b)可通过挤压任一中心频率ω_l附近区间

的值来获得，即

其中，a_k为离散的尺度，且Δa_k＝a_k-a_k-1，T_s(ω_l,b)即为挤压重排后的高分辨、高聚集性时频分布；

步骤二、在高分辨率时频平面上进行Catte模型各向异性扩散滤波，Catte模型为：

其中，U(x,y,t)表示t时刻(x,y)处的像素值，U₀(x,y)表示原始图像，▽·()是散度算子，▽是梯度算子；G_σ是尺度因子为σ的高斯核函数，C是扩散系数，通常选取非负单调递减函数。Perona和Malik提出了两个扩散系数函数：

和

式中k是一个常数，可根据处理的对象选取合适的整数；|▽G_σ*U(x,y,t)|是尺度因子为σ时常系数热扩散方程的梯度模；

将公式(5)所得的时频分布T_s(ω_l,b)看作是一幅图像，对其做各向异性扩散滤波，代入公式(6)中，得：

由此得到滤波后的时频平面，则有效信号的能量聚集点将会得到显著增强，而随机噪声则被很大程度地压制，特别是高频随机噪声被大幅消减；如果含噪信号中随机噪声较强，在滤波后的时频平面上仍会有处在中频的一部分噪声残留下来，因此为了得到更好的消噪效果，此种情况下需借助时频域带通滤波；

步骤三、对滤波后的时频分布做同步挤压小波重构，即可得到滤波后的时域波形：

由公式(7)可得经Catte模型滤波后的时频分布，记作T_s′(ω_l,b)，采用如下所示的同步挤压重构公式：

其中，

Re表示取实部操作；x(b)即为时频面滤波后重构出来的时域波形，即公式(1)中有效信号x(t)的估计。

本发明首次提出在同步挤压时频变换域上采用Catte各向异性扩散滤波模型对地震勘探信号进行增强的方法，同时也是对地震勘探随机噪声有力压制的方法。同步挤压小波变换是目前时频分辨率与时频聚集性相当高的时频分析方法，通过该变换得到的时频分布中有效信号的能量被高度地聚集在一起，随机噪声能量则散落在时频平面内。将所得的时频平面看作是一幅图像，对该图像进行Catte模型各向异性扩散滤波，则有效信号的能量聚集点会显著增强，同时随机噪声能量则被很大程度地消减。对滤波后的时频分布，采用同步挤压小波变换重构过程即可得到有效信号的估计。

该方法充分利用了同步挤压时频变换的高分辨、高聚集特性以及Catte模型各向异性扩散滤波能够兼顾特征保持与噪声消除的能力，以期为地震勘探资料中有效信号的增强以及随机噪声的消减提供更先进的滤波技术。通过对模拟地震信号及人工合成地震记录进行测试，验证了本发明所提出方法的优越性，后续对实际地震资料的滤波处理也同样得到了相当理想的结果。

附图说明

图1为本发明中的同步挤压小波变换时频域各向异性扩散滤波方法流程图。

图2为单个Ricker子波波形图，包括不含噪波形和加噪波形。图2(a)为不含噪的Ricker子波；图2(b)为信噪比大约为15dB的Ricker子波。

图3为单个子波同步挤压时频分布图，包括不含噪、含噪子波时频图以及同步挤压时频域Catte模型各向异性扩散滤波后的时频分布图。图3(a)为不含噪Ricker子波同步挤压时频分布图；图3(b)为信噪比大约为15dB的Ricker子波同步挤压时频分布图；图3(c)为含噪Ricker子波同步挤压时频域各向异性扩散滤波时频分布图。

图4为两分量Ricker子波波形图，包括不含噪波形和加噪波形图。图4(a)为不含噪的两分量Ricker子波；图4(b)为信噪比大约为5dB的两分量Ricker子波。

图5为两分量Ricker子波同步挤压时频分布图，包括不含噪、含噪子波时频图以及同步挤压时频域Catte模型各向异性扩散带通滤波后的时频分布。图5(a)为不含噪两分量Ricker子波同步挤压时频分布图；图5(b)为信噪比大约为5dB的两分量Ricker子波同步挤压时频分布图；图5(c)为含噪两分量Ricker子波同步挤压时频域各向异性扩散带通滤波时频分布图。

图6为时频峰值滤波(TFPF)和本发明方法对两分量Ricker子波的滤波波形对比图及其局部放大图。图6(a)为含噪波形、不含噪波形、TFPF波形和新方法滤波波形对比图；图6(b)为第一个子波主信号段四种波形对比图；图6(c)为第二个子波主信号段四种波形对比图；图6(d)为第一个子波非主信号段四种波形对比图；图6(e)为第二个子波非主信号段四种波形对比图。

图7为TFPF和本发明方法对两分量Ricker子波滤波后的频谱对比图及其局部放大图。图7(a)为含噪信号、不含噪信号、TFPF信号和新方法滤波信号的频谱对比图；图7(b)为中低频、中频部分频谱对比图；图7(c)为高频部分频谱对比图。

图8为TFPF和本发明方法对信噪比大约为-4dB人工合成地震记录的滤波实验结果图。图8(a)为不含噪的合成地震记录；图8(b)为信噪比大约为-4dB的合成地震记录；图8(c)为本发明方法滤波记录；图8(d)为TFPF记录；图8(e)为本发明方法滤波后的差记录；图8(f)为TFPF差记录。

图9为TFPF和本发明方法对信噪比为-12dB人工合成地震记录的滤波实验结果图。图9(a)为不含噪的合成地震记录；图9(b)为信噪比大约为-12dB的合成地震记录；图9(c)为本发明方法滤波记录；图9(d)为TFPF记录；图9(e)为第23道波形对比；图9(f)为第一个子波波形对比；图9(g)为第二个子波波形对比。

图10为TFPF和本发明方法对实际地震记录的滤波结果图。图10(a)为原始记录；图10(b)为本发明方法滤波记录；图10(c)为TFPF记录。

图11为实际地震记录中单道信号的同步挤压时频分布以及同步挤压时频域Catte模型各向异性扩散滤波后的时频分布图。图11(a)为第49道信号的同步挤压时频分布；图11(b)为第49道信号的时频域各向异性扩散滤波时频分布；图11(c)为第67道信号的同步挤压时频分布；图11(d)为第67道信号的时频域各向异性扩散滤波时频分布。

图12为实际地震记录及其滤波结果中第49道信号的频谱对比图。图12(a)为整道信号滤波前后的频谱对比；图12(b)为该道信号频谱对比的局部放大图。

具体实施方式

下面结合附图对本发明的技术方案做详细叙述。

参照图1，一种增强地震勘探信号的高性能时频域滤波方法，包括以下步骤：

步骤一、对含噪的地震信号做同步挤压小波变换，得到高分辨率时频分布：

设含噪信号为：

s(t)＝x(t)+n(t) (1)

x(t)为有效信号，n(t)是加性随机噪声；

对s(t)做连续小波变换，参照公式(2)如下：

其中，a是尺度因子，b是平移因子，ψ(t)为小波基函数，

为ψ(t)的共轭函数；根据帕塞瓦尔定理，公式(2)的频域表达式为：

其中，ξ是角频率，

是s(t)的傅里叶变换，ψ(aξ)为ψ(t)的傅里叶变换；通过该式就将时间域的褶积形式变换为频率域的乘积形式，通过对小波系数求导可初步估计瞬时频率，即

其中，i是虚部单位；这样就把时间—尺度平面(b,a)转换到时间—频率平面(b,ω_s(a,b))，频率ω和尺度a均被离散化，因此仅在离散点a_i{a_i-a_i-1＝(Δa)_i}计算小波系数W_s(a,b)；同步挤压小波变换中频率重排的实现过程为：假设信号的采样频率为cf，小波变换时所取尺度的个数为N，那么重排后的时频谱上，离散的频率值就为ω_l＝l*cf/N,l∈[1,N]；小波系数的同步挤压值T_s(ω_l,b)可通过挤压任一中心频率ω_l附近区间

的值来获得，即

其中，a_k为离散的尺度，且Δa_k＝a_k-a_k-1，T_s(ω_l,b)即为挤压重排后的高分辨、高聚集性时频分布。

所得的时频分布如图3、图5所示。图3(a)、(b)分别为单个Ricker子波(主频为25Hz)不含噪(图2(a)所示)与含噪时(图2(b)所示)的同步挤压时频分布；图5(a)、(b)分别为两分量Ricker子波(主频分别为20Hz和15Hz)不含噪(图4(a)所示)与含噪时(图4(b)所示)的同步挤压时频分布；可以看出，无论随机噪声是较弱还是较强，时频平面上有效信号的能量聚集度都很高，频率分量曲线都很精细，只是噪声较强时时频曲线受到一定程度的干扰而不清晰。对比图3(b)和图5(b)，可以看出，后者中高频部分随机噪声散布得较多。

步骤二、在高分辨率时频平面上进行Catte模型各向异性扩散滤波，Catte模型为：

其中，U(x,y,t)表示t时刻(x,y)处的像素值，U₀(x,y)表示原始图像，▽·()是散度算子，▽是梯度算子；G_σ是尺度因子为σ的高斯核函数，C是扩散系数，通常选取非负单调递减函数；Perona和Malik提出了两个扩散系数函数：

和

式中k是一个常数，可根据处理的对象选取合适的整数；|▽G_σ*U(x,y,t)|是尺度因子为σ时常系数热扩散方程的梯度模；

将公式(5)所得的时频分布T_s(ω_l,b)看作是一幅图像，对其做各向异性扩散滤波，代入公式(6)中，得：

由此得到滤波后的时频平面，则有效信号的能量聚集点将会得到显著增强，而随机噪声则被很大程度地压制，特别是高频随机噪声被大幅消减；如果含噪信号中随机噪声较强，在滤波后的时频平面上仍会有处在中频的一部分噪声残留下来，因此为了得到更好的消噪效果，此种情况下需借助时频域带通滤波。

如图3(c)所示的15dB单个Ricker子波同步挤压时频域Catte模型滤波后的时频分布，时频平面较为干净，无需做进一步的处理；而图4(b)所示的两分量Ricker子波信噪比相对较低，其同步挤压时频平面上(图5(b)所示)随机噪声能量较强，因此在该时频面上进行Catte模型滤波后需借助带通滤波将增强的有效信号能量分布截取出来，如图5(c)所示。

步骤三、对滤波后的时频分布做同步挤压小波重构，即可得到滤波后的时域波形：

由公式(7)可得经Catte模型滤波后的时频分布，记作T_s′(ω_l,b)，采用如下所示的同步挤压重构公式：

其中，

Re表示取实部操作；x(b)即为时频面滤波后重构出来的时域波形，即公式(1)中有效信号x(t)的估计。图6展示出了本发明提出的滤波方法与TFPF方法对两分量含噪子波处理后的波形对比。可以看出本发明方法在对有效信号的保持方面(波峰、波谷)更具优势，在随机噪声压制方面也优于TFPF。

滤波方法的性能测试：

一、该步骤需从对人工合成地震数据模型的滤波实验着手，分析滤波结果，然后对实际地震数据进行处理。

首先是对单道加噪的模拟地震信号按照步骤一、二、三进行处理，然后对加噪的合成同相轴地震记录进行处理。对合成记录的滤波处理其实也是按照单道信号的处理方法完成的。一般情况下，采用Ricker子波模拟地震反射信号，其数学表达式为：

x(t)＝[1-2(πf₀t)²]exp[-(πf₀t)²] (9)

以上的滤波实验应在不同的信噪比条件下进行，这样对方法性能的测试将会比较全面。即对于单道的模拟地震信号以及含有同相轴的合成记录，可加入不同强度的随机噪声，使其信噪比从高到低取值。

滤波实验过程中，加入随机噪声使得两分量Ricker子波信噪比大约为5dB，使得两同相轴的人工合成记录信噪比大约为-4dB和-12dB。图6(a)—(e)为两分量Ricker子波滤波前后的波形对比及其局部放大图；图8(a)—(f)展示了不含噪的两同相轴合成记录及其加噪记录(-4dB)，以及分别采用本发明方法和TFPF方法滤波后的记录及差记录；图9(a)—(d)展示了不含噪的两同相轴合成记录及其加噪记录(-12dB)，以及分别采用本发明方法和TFPF方法滤波后的记录。

这里需要说明的是，虽然实际地震信号与模拟地震信号在形态上存在一些差异，但由于本发明所设计的滤波方法参数鲁棒性好，具有普适性，因此在处理实际地震资料时仍会取得与模拟地震信号实验结论一致的效果。图10(b)、(c)分别为本发明提出的滤波方法和TFPF方法对一个实际地震记录(图10(a)所示)的处理结果。

二、对地震勘探信号及记录的滤波结果进行定性与定量分析。

对模拟地震信号和人工合成地震记录进行滤波处理后，由于模拟数据实验中有不含噪的情况作为参考标准，所以对滤波结果均能够进行定性与定量分析。而对实际地震资料的处理，由于没有不含噪的情况作参考，因此一般主要做定性分析。

定性分析主要有波形对比、谱分析、做差等方式；定量分析主要是计算滤波前后的信噪比和均方误差。将滤波前后的波形与理想波形作对比，波形越接近理想情况，说明滤波效果越好。波形对比如图6(a)—(e)以及图9(e)—(g)。从这些图中可以看出，本发明方法滤波后的波形更接近于不含噪波形，说明其对有效信号的保持和对随机噪声的消减能力均优于TFPF方法。

谱分析方法中最普遍的是采用傅里叶变换做频谱分析，主要观察中、高频随机噪声的消减情况(频谱的平坦程度)以及处于低频的有效信号的保留情况(频谱幅度的保持)。频谱对比如图7(a)—(c)以及图12(a)、(b)，可以看出，本发明方法滤波所得的信号，其频谱中、高频部分非常平坦，而TFPF方法所得的信号，其频谱一直存在抖动现象，说明前者对随机噪声的压制能力更强。而观察频谱低频部分的幅值保持(如图12(a)所示)，可以看出TFPF仍稍逊一筹。

做差法是用原始含噪记录减去滤波后的记录得到差记录，这个差记录主要是滤除的随机噪声和一小部分损失掉的有效同相轴。如果差记录中损失的有效同相轴越少，那么说明滤波方法对有效信号的保持能力越强。做差比较如图8(e)、(f)，可以看出TFPF差记录中损失的有效同相轴较为明显。图9(c)、(d)是对比滤波后同相轴的清晰度与连续性，在信噪比很低时，可以看出TFPF处理后的记录中同相轴出现了不连续的现象(图中矩形框所示)，而本发明方法处理后的记录中同相轴更清晰连贯。

图10为实际地震记录经两种方法处理后的效果对比，可以很明显地看出，本发明方法处理后的记录中，有效同相轴清晰度高、连续性好，随机噪声消减力度大，而TFPF处理后的记录中同相轴仍不够清晰。由于该实际记录含噪并不是太多(图11(a)—(d)所示的时频分布可以看出随机噪声能量很弱)，所以无需在同步挤压时频域Catte模型滤波后做进一步的处理。通过图12(a)、(b)所示的单道信号频谱对比，可以看出，本发明方法在对中、高频随机噪声的压制方面确实具有很大的优势，而TFPF方法对中、高频随机噪声的压制能力较为有限，其所得信号频谱直至高频部分仍不够平坦，存在小幅抖动。

定量分析中通过计算信噪比，并与原始信噪比作对比，可得滤波后信噪比的提升程度；计算均方误差并与原始均方误差作对比，可得滤波后均方误差的减小程度。计算信噪比和均方误差的数学公式如下(以二维地震记录为例)：

公式中，x(i,j)为不含噪记录，

为滤波后的记录，N为采样点数，T为道数。

表1是对不同信噪比的人工合成记录及其滤波结果所做的定量分析，即通过公式(10)、(11)计算出滤波前后的信噪比和均方误差。

表1不同信噪比的人工合成记录及其滤波结果评价指标

本发明提出的方法既充分利用了同步挤压时频变换的高分辨、高聚集特性，即能够将有效信号的能量高度地聚集在一起且瞬时频率曲线很精细，又利用了Catte模型各向异性扩散滤波能够兼顾有效信号特征保持与随机噪声消减的能力。该方法易于实现，运算速度快，参数鲁棒性强，适用于不同噪声强度下的地震勘探资料处理。因此，本发明所述的地震勘探信号增强方法具有良好的推广应用前景和价值。

本领域的技术人员应当理解，现今油气矿藏资源的埋藏愈来愈深，勘探环境愈来愈复杂，勘探难度随之愈来愈大，所采集到的地震勘探资料中有效信号的识别难度也愈来愈大。为了给资料解释工作提供高质量的解释素材，在资料处理环节中必须采用性能优越的滤波方法对其中的有效信号进行恢复提取，对各种噪声进行消除，而随机噪声分布于整个资料当中，是主要存在的噪声类型之一。鉴于此，需采用能够显著增强有效信号同时大幅消减随机噪声的高性能滤波方法处理地震勘探资料，才能为后续的资料解释环节打下坚实的基础。

Claims (1)

1.一种增强地震勘探信号的高性能时频域滤波方法，其特征在于，包括以下步骤：

步骤一、对含噪的地震信号做同步挤压小波变换，得到高分辨率时频分布：

设含噪信号为：

s(t)＝x(t)+n(t) (1)

x(t)为有效信号，n(t)是加性随机噪声；

对s(t)做连续小波变换，参照公式(2)如下：

其中，a是尺度因子，b是平移因子，ψ(t)为小波基函数，

为ψ(t)的共轭函数；根据帕塞瓦尔定理，公式(2)的频域表达式为：

其中，ξ是角频率，

是s(t)的傅里叶变换，ψ(aξ)为ψ(t)的傅里叶变换；通过该式就将时间域的褶积形式变换为频率域的乘积形式，通过对小波系数求导可初步估计瞬时频率，即

其中，i是虚部单位；这样就把时间—尺度平面(b,a)转换到时间—频率平面(b,ω_s(a,b))，频率ω和尺度a均被离散化，因此仅在离散点a_i{a_i-a_i-1＝(Δa)_i}计算小波系数W_s(a,b)；同步挤压小波变换中频率重排的实现过程为：假设信号的采样频率为cf，小波变换时所取尺度的个数为N，那么重排后的时频谱上，离散的频率值就为ω_l＝l*cf/N,l∈[1,N]；小波系数的同步挤压值T_s(ω_l,b)可通过挤压任一中心频率ω_l附近区间

的值来获得，即

其中，a_k为离散的尺度，且Δa_k＝a_k-a_k-1，T_s(ω_l,b)即为挤压重排后的高分辨、高聚集性时频分布；

步骤二、在高分辨率时频平面上进行Catte模型各向异性扩散滤波，Catte模型为：

其中，U(x,y,t)表示t时刻(x,y)处的像素值，U₀(x,y)表示原始图像，

是散度算子，

是梯度算子；G_σ是尺度因子为σ的高斯核函数，C是扩散系数，通常选取非负单调递减函数；Perona和Malik提出了两个扩散系数函数：

和

式中k是一个常数，可根据处理的对象选取合适的整数；

是尺度因子为σ时常系数热扩散方程的梯度模；

将公式(5)所得的时频分布T_s(ω_l,b)看作是一幅图像，对其做各向异性扩散滤波，代入公式(6)中，得：

由此得到滤波后的时频平面，则有效信号的能量聚集点将会得到显著增强，而随机噪声则被很大程度地压制，特别是高频随机噪声被大幅消减；如果含噪信号中随机噪声较强，在滤波后的时频平面上仍会有处在中频的一部分噪声残留下来，因此为了得到更好的消噪效果，此种情况下需借助时频域带通滤波；

步骤三、对滤波后的时频分布做同步挤压小波重构，即可得到滤波后的时域波形：

由公式(7)可得经Catte模型滤波后的时频分布，记作T_s′(ω_l,b)，采用如下所示的同步挤压重构公式：

其中，

Re表示取实部操作；x(b)即为时频面滤波后重构出来的时域波形，即公式(1)中有效信号x(t)的估计。

Images