小波包变换(彭玉华,1999; 刘明才,2005)推广了小波分析,通过一组低、高通正交滤波器H和G,对给定信号进行高频和低频分量多层次分解,使信号的时频特征分解的更加精细,适合于描述和刻画非平稳信号。

令正交小波基的滤波器系数分别为{h_k}_k∈Z和{g_k}_k∈Z,并将尺度函数φ(t)改记为u₀(t), 小波函数φ(t)改记为u₁(t), 由关于φ(t)和φ(t)的两尺度方程(1)式,引进小波包的概念(2)式

{u₀(t)=2^1/2∑_k∈Zh_ku₀(2t-k)

u₁(t)=2^1/2∑_k∈Zg_ku₁(2t-k).(1)

{u_2n(t)=2^1/2∑_k∈Zh_ku_n(2t-k)

u_2n+1(t)=2^1/2∑_k∈Zg_ku_n(2t-k).(2)

式(2)所定义的函数集合{u_n(t); n∈Z⁺}为由u₀=φ 所确定的小波包(正交小波包)。

定义闭包空间Uⁿ_j=span{2^j/2u_n(2^jt-k)}_k∈Z^-, n∈Z⁺, j∈Z。 设gⁿ_j(t)∈Uⁿ_j, gⁿ_j(t)可表示为

gⁿ_j(t)=2^j/2∑_k∈Zdⁿ_jku_n(2^jt-k).(3)

则小波包分解的分解算法为

{d²ⁿ_jl=∑_k∈Zh^-_k-2ldⁿ_j+1,k

d²ⁿ⁺¹_jl=∑_k∈Zg^-_k-2ldⁿ_j+1,k.(4)

重构算法为

dⁿ_j+1,l=∑_k∈Z(h_l-2kd²ⁿ_jk+g_l-2kd²ⁿ⁺¹_jk).(5)

对应L²(R)空间的子空间的每一种划分形式就能得到一种小波包基,在众多的小波包基中,笔者采用最优基,可在所需研究的尺度上提供较好的时频分辨率。

给定一个序列的代价函数,在小波库{2^j/2u_n(2^jt-k); n∈Z⁺, k, j∈Z}的所有小波包基中寻找代价函数最小的基即为最优基。

代价函数可定义为任何关于序列的实函数M,而实函数M的选取,最好是那些能够测得系数集中度的可加性代价函数。当序列中的系数大致在相同尺度时,M会大,反之,M会小。对实数序列附加信息成本函数M采用香农熵,则有

M[d_jl(i)]=-∑_id_jl(i)²log[d_jl(i)²],

且0lg0=0.(6)

对小波包分解的每个节点计算成本函数,通过比较某一节点与其子节点的值,获取熵值最低的基。由最小熵确定的最优基,对给定母小波和成本函数提供了最佳选择。

基于小波的降噪算法大多遵循以下步骤:对信号变换,修改小波系数以减少或消除难以确定的系数,重构系数序列获得所需信号。假设含噪有限长信号可表示为

x(i)=s(i)+z(i)i=0,1,…,N.(7)

其中,x(i)是记录信号; s(i)是感兴趣的信号; z(i)是噪声。 去噪的最终目标是从x(i)中找到s(i)的一个估计信号s^{^}(i),并尽可能与之类似。具体计算方法如下:

(1)小波包变换:计算小波包变换X=Wx,得到与基函数相关的小波包系数序列d_jl,该系数序列也被噪声污染,可表示为

d_jl|_x=d_jl|_s+d_jl|_z.(8)

式中,d_jl|_x表示记录信号的系数集合; d_jl|_s表示感兴趣信号的系数集合; d_jl|_z表示噪声相关的系数集合。

(2)小波包系数阈值处理:感兴趣信号的能量集中在小波域的少量系数上,其取值必然大于能量分散的噪声的小波包系数值。假设每个小波包系数都受到污染,对它们进行阈值处理,处理的主要方法有硬阈值和软阈值(Donoho,1995):

硬阈值:S^{^}=T_h(X,t)={X |X|>=t

0 |X|<t.(9)

软阈值: S^{^}=T_h(X,t)

={sgn(X)(|X|-t)|X|>=t

0 |X|<t. (10)

阈值处理以最小化估计系数与期望系数之间的均方差最小化为目的。从L₂范数误差最小的观点出发,硬阈值处理方法的估算系数d_jl|_s优于软阈值方法,但其重构信号的光滑性不如软阈值方法。

(3)小波包逆变换:计算小波包逆变换s^{^}=W^-1S^{^}得到恢复信号s^{^}。

用合成信号模拟震源P波和S波的相继到达情况,假设独立点震源的地震波强度随时间满足指数衰减规律(Aki,Richards,1980; Slawmomir,et al,1994)

m(t,ω,Q,t₀)=A₀exp[-πω(t-t₀)/Q].(12)

这里,Q为常数,取30; ω为频率; t₀为波的初始到达时间。依据这个关系,可以构造一个包含多个频率的合成信号

f(t)={0(0≤t<t₀₁)

∑³₁m_isinω_i(t-t_0i)(t₀₁≤t<t₀₂)

∑⁶₁m_isinω_i(t-t_0i)(t₀₂≤t<10).(13)

其中,m_i=m(t,A_0i,ω_i,t_0i),信号中的基本参数按表(1)给出,采样频率100 Hz,信号终止时刻为10 s。合成的原始信号及其被信噪比为10 dB的高斯白噪声污染后的信号如图1所示。

图1 原始合成信号及其加噪后信号(a)具有指数衰减规律的原始合成信号;(b)加入高斯白噪声后的信号(信噪比为10 dB)
Fig.1 Original synthetic signal and its contaminated signal(a)Original synthetic signal with exponential decay; (b)Signal contaminated by Gaussian white noise(SNR is 10 dB)

表1 模拟地震波到达的合成信号取值参数表
Tab.1 Parameter values of the synthetic signal simulating seismic wave

笔者利用最优小波包基变换,采用节点相关阈值去噪,去噪信号与加噪前的原始信号非常相似,不但信噪比非常高,第一脉冲到达前也没有产生人为引入的信号,而且第一对脉的到时几乎没有失真,去噪效果最好(图2c)。

图2 加噪信号的去噪结果
Fig.2 De-noising results of the contaminated signal

图3a为SDAES 数字声发射仪采集的近距离微震实验信号,提取信号前2 s(1 250个样本)的数据作为噪声。

去噪结果表明:3种方法对第一脉冲峰值没有产生时间延迟或提前; 采用小波变换默认阈值去噪,信噪比改善不明显,频谱分析表明噪声信号频带与微震信号频带相近,该方法去噪原理类似低通滤波器,因此去噪效果不好(图3b); 采用小波变换分层阈值去噪,信噪比有所提高,但结果仍不够理想(图3c); 最优小波包基变换采用节点阈值去噪,信噪比有很大的提高,虽然结果的幅值较其他方法要低,但不影响到时的取值(图3d),即该方法对低频噪声的去噪效果明显优于以小波变换为基础的其他去噪方法。

图3 实验信号及去噪结果(a)声发射仪采集的微震实验信号;(b)多尺度小波变换默认阈值去噪结果;(c)小波变换分层相关阈值去噪结果;(d)最优小波包基变换节点相关阈值去噪结果
Fig.3 Experimemtal signal and its de-noising signal(a)Experimental micro-seismic signal collected by SDAES digital acoustic emission instrument; (b)Signal denoised by multi-scale WT using a global threshold;(c)Signal denoised by WT using a level-dependent threshold;(d)Signal denoised by the optimum WPT using a node-dependent threshold

图4a为2007年3月25发生在日本能登半岛地震的记录波形图,数据源于NIED(K-NET)强震观测台站。地震图信噪比不足以精确确定到达时间,应用上述3种方法消噪,采用震前10 s(1000个样本)记录的环境噪声数据来描述。

类似前面的结果,最优小波包基变换采用节点相关阈值去噪在信噪比上优于小波变换为基础的去噪结果(图4b～d)。

图5是图4去噪结果的第一到达脉冲的局部放大图。由图5可见:尽管3种方法得到的第一脉冲的峰值幅值不同,但位置一致,而第一脉冲的到达时间存在明显差异。根据上述合成信号和实验信号的分析结果以及地震波形图的频谱分析推断:应用小波变换方法去噪,信号的到时有可能存在人为引入的信号波动误差; 采用最优小波包基变换节点阈值,去噪后信号第一脉冲的到时位置更准确。

分析:小波包变换将信号在低频部分和高频部分都进行了精细分解; 采用最优基能够根据该信号的特征自适应地选择相应的频带,使之与信号频谱相匹配; 采用节点相关阈值消噪使信号的去噪效果最佳。这表明本文采用的去噪方法为后续分析提供的信号更可靠。