样条小波的构造方法及其应用

周义明　张　超　张化民

摘 要：从样条函数定义入手,给出了构造样条小波的具体方法。为了解决小波分析对奇异信号的奇异点定位困难的问题,重点介绍了基于余弦基展开式的零点对称的二进样条小波、零点反对称的二进样条小波构造方法,并提出了具体的实现算法及计算过程。通过与正交样条小波的比较,对称样条小波能对阶跃信号和脉冲信号奇异点进行更好的定位,实验结论也进一步验证了该样条小波对分析奇异信号具有较大优势。

关键词：样条函数；构造方法；二进样条小波；信号定位；奇异性

中图分类号：TN911 文献标识码：B

文章编号：1004-373X(2009)01-036-05

Construction Methods and Application Based on Spline Wavelets

ZHOU Yiming1,ZHANG Chao2,ZHANG Huamin3

(1.Beijing Institute of Petro-chemical Technology,Beijing,102617,China;

2.College of Electronics and Information Engineering,Tongji University,Shanghai,201804,China;

3.Hunan Tax College,Changsha,410116,China)

Abstract：Based on cardinal spline function,the contribution gives the construction ways.In order to overcome the matter about locating singularity position,the paper introduces the ways of symmetry dyadic spline wavelets on zero,anti-symmetry dyadic spline wavelets on zero,at the same time,puts forward to the transform arithmetic.Analyzing the real fault step signals and pulse signals by means of the constructed spline wavelets,the paper proves that the location of symmetry(anti-symmetry) dyadic spline wavelet is better than orthorhombic wavelet,experiment testify the result is correct.

Keywords：spline function;construction method;dyadic spline wavelets;signal location;singularity

0 引 言

小波变换因为同时具有时-频分析能力,成为了近代信号处理的主要工具。样条小波作为一种小波函数重要分支,由于具有更加优良的时频特性,而被广泛应用。然而,国内外的小波专著主要是对小波函数构造理论的推导及证明,需要较强的数学理论知识,对于大多数工程技术人员而言,主要关心的是小波函数构造的方法、结论、算法实现及应用范围。基于此,从样条函数入手,结合对小波正交性、正则性、消失矩、紧支性和对称性要求的不同,构造出了不同性质的样条小波、实现算法和应用范围。

全文的主体内容大致分为3个方面：

(1) 引入基本样条函数的定义及其性质；

(2) 利用小波构造理论给出了二进样条小波(反对称和对称)、正交样条小波的构造过程和结论,同时给出了实现算法；

(3) 结合作者在实际项目中的应用,选用不同的 3阶二进样条小波对突变信号进行了分析与比较,给出了反应该突变信号的特征。

1 B样条函数

样条函数广泛应用于曲线、曲面拟合等领域,小波分析出现后,利用样条函数构造样条小波取得了巨大成就。

1.1 B样条函数的定义

设有节点序列{ti} n+k,称函数:

B i,k(x)=(t i+k-ti)［ti,t i+1,…,t i+k］(t-x) k-1+(1)

为关于节点序列{ti}的第i个k阶(k-1次)B样条函数［1］。

说明：k阶半截幂函数f(t)=(t-x) k-1+(x为参数)的k阶差商乘上一个数(t i+k-ti)就是k阶(k- 1次)B样条函数。

由式(1)可得,当k=1时,第i个1阶(k-1次) B样条函数为:

B i,1(t)=1, ti≤t＜t i+1

0,其他(2)

由于i只是反应了样条函数在t轴的平移,不妨将B样条函数设为以0为对称中心,且节点距离为1,并重记为:

N1(t)=1, -1/2≤t＜1/2

0，其他(3)

1.2 B样条函数性质

对于m阶B样条函数Nm有以下性质［2,3］：

(1) 样条函数的支集：supp Nm=［-m/2,m/2］；

(2) 对于每个f∈C,有:

∫∞ -∞f(x)Nm(x)dx=∫ 1/2 -1/2…∫ 1/2 -1/2f(x1+x2+…

+xm)dx1dx2…dxm

(3) 对称性：Nm(-t)=Nm(t)；

(4) Nm(t)=tm-1N m-1(t)+m-tm-1N m-1(t-1)；

(5) 对所有t,∑∞k=-∞Nm(t-k)=1；

(6) Nm(t)=N m-1*N1(t)；

结合式(2)和性质(4),(6)可以求出任意阶B样条函数表达式。

2 样条小波

由样条函数可构造多种性质的样条小波,结合信号分析的目的和要求,主要从正交性、消失矩、紧支性、对称性等方面考虑,通过小波构造理论来构造不同性能的小波母函数。

2.1 二进样条小波

函数Ψ(t)∈L1(R)∩L2(R)称为一个二进小波,若存在常数0＜A≤B＜∞使得对于笑亍蔙-{0}有:

A≤∑j∈ZΨ(2jω)2≤B(4)

二进小波变换只是将连续小波变换尺度S以{2j}采样,平移参数不变,为了适宜对信号进行边沿检测或突变检测,往往小波的对称性作为构造重点,并以卷积形式定义：{W 2jf(t)} j∈Z叫做f的二进小波变换,其中:

W 2jf(t)=f*Ψ 2j(t)=12j∫Rf(x)Ψ(t-x2j)dx(5)

在式(5)中,根据卷积定理有:

W 2j(ω)=(ω)(2jω)(6)

2.1.1 二进样条小波的构造

在二进样条小波构造之前,先给出母二进小波φ(t)满足的条件：

(1) φ(t)是偶函数,即φ(-t)=φ(t)；

(2) limω→0 (ω)=1, limω→∞ (ω)=0；

(3) (2ω)=H(ω)(ω),且H(ω)满足：H(ω)∈L2(［-π,π］),H(ω)≤1,H(ω)又称为母二进小波φ(t)的生成元。

二进小波变换主要是对突变信号进行边沿和突变点分析,信号有2种性质的突变：阶跃突变；脉冲突变。因此在构造样条小波时,就对应构造出了两种不同性质的二进样条小波：反对称二进小波；对称二进小波。

(1) 反对称二进小波构造

设φ(t)是一个二进小波,H(ω)为其生成元,N是任意固定的正整数,(a1,a2,…,aN)是满足∑Nk=1ak＜1的任意N维实数向量,则o(2ω)=Go(ω)(ω)是一个定义在L2(R)中的二进小波函数,其中:

Go(ω)=isgn(ω)［1+∑Nk=1(akcos kω)］1-H(ω)2(7)

当N=1,ak=0时,有:

G0o(ω)=isgn(ω)1-H(ω)2(8)

由式(7)可知,Go(ω)是一个余弦基的展开式,可将其认为是一个对称函数的准Fourier级数展开式,当选取不同的N和ak时可构造出不同滤波系数的二进小波。

对于Go(ω)可得以下一些性质：

① 由于式(7)中k∈Z,必有Go(ω)∈ L2(［-π,π］)；

② Go(ω)是奇函数；

③ 记Go(ω)=∑n∈Zgne -inω,则有gn=-gn,g00=0；

④不同N,ak构造的反对称二进小波滤波系数之间有关系：

gn=g0n+12∑Nk=1ak(g0 n-k+g0 n+k)(9)

(2) 对称二进小波构造

与反对称二进小波构造条件一样,只是取：

Ge(ω)=［1+∑Nk=1(akcos kω)］1-H(ω)2(10)

当N=1,ak=0时，有：

G0e(ω)=1-H(ω)2(11)

相应的性质如下：

① Ge(ω)∈L2(［-π,π］)；

② Ge(ω)是偶函数；

③ 记Ge(ω)=∑n∈Zgne -inω,则有gn=-gn；

④ 不同N,ak构造的对称二进小波滤波系数之间的关系：

gn=g0n+12∑Nk=1ak(g0 n-k+g0 n+k)(12)

以3阶中心B样条函数作为母二进小波φ(t)：

根据样条函数的定义和性质,由Nm(t)=N m-1砃1(t)及卷积定理可得：

(2ω)=m(2ω)=(sin ωω)m=

(cosω2)mm(ω)=H(ω)m(ω)(13)

因此：

H(ω)=(cosω2)m(14)

现分析其是否满足母二进小波的3个条件：

① 由于φ(t)是以零为中心的B样条函数,所以其为偶函数；

② limω→0(2ω)= limω→0sin ωωm=1,同理 limω→∞ (2ω)=0；

③ 由式(14)有：当m=2k,即m为偶数时,H(ω)满足条件,可作为生成元；当m=2k+1,即m为奇数时,不是以2π为周期的,若在式(14)右边乘以一指数因子以保证其是以2π为周期的,则时域的φ(t)将会产生零点偏移,破坏了φ(t)的偶对称性,Mallat构造的3阶二进样条小波［4］就属于此类。为解决这一矛盾,许传祥等通过扩展φ(t)的时间轴使得其频域轴压缩,从而构造了极好的二进样条小波,并给出了3阶和4阶二进样条小波的滤波器系数［5］。

图1为二者用3阶样条函数构造的小波函数对一阶跃信号的分析,一共进行了9层小波分析,其中 图1(a)为被分析信号；图1(b)为由膨胀 φ(t)的时间轴构造的小波对信号的分析； 图1(c)为由Mallat构造的小波对信号的分析。由图可见,膨胀φ(t)的时间轴构造的小波分解系数不发生时间的偏移,而由Mallat构造的小波随着分解层数的增加逐步偏离突变点,但紧支性方面对称性样条小波不如Mallat二进小波。

图1 不同样条小波对突变信号的分析

2.1.2 二进小波变换的实现

由于二进样条小波采用余弦基构造,而且其变换通过卷积形式,因此不符合MRA理论,不能采用快速Mallat算法实现。à trous算法(多孔算法)能实现这一二进小波变换。

设原始信号为{s(n)} n∈［0,N),可以证明必存在f(t),使得:

s(n)=＜f(t),φ(t-n)>(15)

不妨设f(n)=Sd1f,f(n)为t=n时刻的值,对于采样较密集的信号可直接设s(n)=Sd1f,则对每一个尺度2j,Sd 2jf可分解为Sd 2 j+1f和Wd 2 j+1f,程序公式如下：

while(j≤J)

Sd 2 j+1f=Sd 2jf砲n

Wd 2 j+1f=Sd 2jf砱n

j=j+1

end

分解后的信号长度为JN个细节系数,N个概貌系数。

2.2 正交样条小波

正交小波由于在尺度和平移方面都是正交的,不含冗余项,被广泛应用于图像或数据压缩,因而在构造过程主要考虑正交性。

2.2.1 正交样条小波的构造

正交小波的构造主要采用MRA理论,在构造之前先给出几个定义或定理：

定理1 MRA理论是指满足以下条件的序列闭子个空间{Vj} j∈Z：

(1) …糣2糣1糣0糣 -1糣 -2肌；

(2) ∩j∈ZVj={0}; ∪j∈ZVj=L2(R)；

(3) f(t)∈Vj趂(2t)∈V j+1,j∈Z；

(4) 存在φ(t)∈V0,使得V0= span{φ(t-n),n∈Z}。

由条件(4)可以得到一个引理：

若满足：

∑(ω+2nπ)2=1(16)

则序列{φ(t-n),n∈Z}必为V0空间的一组正交基。

定理2 若{Vj} j∈Z满足MRA理论,尺度函数序列{φ(t-n),n∈Z}为V0空间的正交基,则有二尺度差分方程如下：

尺度:

φ(t)=2∑n∈Zcnφ(2t-n)(17)

小波:

Ψ(t)=2∑n∈Z(-1)nc 1-nφ(2t-n)(18)

可将小波离散表示为：

Ψ j,k(t)=2jΨ(2j-k)(19)

可以证明Ψ(t)的Fourier变换满足：

(ω)=-e -iω/2(ω2+π)(ω2)(20)

其中：

H(ω)=2 -1/2∑cne -inω(21)

根据式(3)和1.2节中性质(6),若取φ(t)=Nm(t)则有：

当m=1时,N1(t)具有正交性,可以构造正交小波基,即Haar小波；

当m≥2时,Nm(t-n)是Riesz基的,并不具备平移正交性,因此在构造小波之前必须对其进行正交化处理,记：

p(ω)=sin(ω/2)ω/2p(22)

Fp(ω)=p(ω+2nπ)2〗 -1/2(23)

#p(ω)=p(ω)Fp(ω)(24)

则可证明存在{V #0}使得{φ #p(t-n),n∈Z}为其正交基,这时再利用定理2就可以构造正交的样条小波基,即有：

(ω)=e -iω/2 m #0(ω/2+π)φ #p(ω/2)(25)

其中：m #0(ω)=m0(ω)［∑n∈Zp(ω+2nπ)2］ 1/2［∑n∈Z

p(2ω+2nπ)2］ -1/2,且m0(ω)=2 -1/2∑n∈Zhne -inω, hn=cn,∑n∈Zp(2ω+2nπ)2=-sin 2pω(2p-1)!d 2p-1dω 2p-1 cot ω［3］。

以时域表示式(24)和式(25)有：

φ #p(t)=2∑n∈Zh #nφ #p(2t-n)(26)

Ψ #p(t)=2∑n∈Z(-1)nh # 1-nφ #p(2t-n)(27)

该小波还具有对称性,但却牺牲了小波的紧支性,二阶和三阶正交样条小波的具体构造过程和其他阶正交样条小波滤波器系数在相关文章中给出［6,7］。

2.2.2 正交样条小波变换的实现

正交样条小波的构造是基于MRA理论,因此可以通过快速Mallat算法实现,其原理为：

信号的小波变换可以看作对信号分别进行高通或低通滤波,利用二尺度差分方程对小波变换进一步推导可得：

x (j)k=∑nh0(n-2k)x (j-1)n(28)

d (j)k=∑ng0(n-2k)x (j-1)n(29)

其中：h0,g0分别为式(25),式(26)中的滤波器系数h #n,(-1)nh # 1-n。

将待小波分解信号通过一个高通滤波器和一个低通滤波器进行滤波,得到一组低频信号和一组高频信号,并且继续对低频信号进行类似分解直到满足要求,每一次分解得到的低频信号和高频信号长度都是原信号长度的一半,两者之和等于原信号的长度,可看作是在滤波后进行隔点采样,分解结果既不冗余,也不损失原信号的任何信息,此即为Mallat算法,将其用电路结构图表示如图2所示。

图2 Mallat算法电路结构示意图

3 实例分析

根据分析信号的目的和要求的不同,选择的小波基和变换方法也是不一样的,基于作者在北京环铁基地做实验所采集的信号,选用不同的小波基对其进行了分析。

实验目的主要是分析牵引接触网的接地故障行波信号,测试线路长度为1 400 m,数据采样频率为 100 MHz,分为在线和离线检测两种情况：

(1) 在线情况是将线路首端接一触发恒压电源 15 V,末端通过电阻(有几种等级)接地,电源触发后,在末端采集故障行波；

(2) 离线情况是用电容球隙放电产生的突发脉冲代替(1)中恒压电源作为信号源,末端仍接地,信号采集端在首端,图3为用虚拟仪器实测的故障行波波形。

图3 两种情况下的实测故障行波

注意实测波形：由恒压电源产生的故障行波是带一定倾斜度台阶上升的,最终达到15 V；由球隙放电产生的故障行波是隔一定间距就产生相反方向的突变脉冲(不严格)信号,第一个极大脉冲是由球隙放电产生,并不是故障行波,振荡波是由球隙逐渐放电过程引起的。由于测试现场噪声较强,采样率较高,数据量较大,为了便于分析,在不破坏原始数据波形成分和形状的条件下,只对其(突变剧烈处)进行部分截取并放大,分别用不同小波基分析直流恒压源故障行波,分解层数为 9级,得到图4,其中(a)为原始信号； (b)为零点反对称小波(N=1,a1=0.28)的分析；(c)为Mallat小波的分析；(d)为零点对称小波(N=1,a1=-0.64)的分析。

图4 直流恒压源故障行波在不同小波基下的分析

由图4可见,对于分析阶跃信号,在低级数分解过程中,由于噪声的影响,无法分辨突变点,提高分解级数(第8级),可以分辨突变点,但是难以精确定位,不同小波基下的分解,定位也是完全不同的,零点反对称小波的分析最佳,对于脉冲性质的信号,零点对称小波的分析最佳,当选用正交小波基时,每级分解后数据量减半,但由于进行了二抽取,完全没有定位性,因此可总结出以下性质：

(1) 用小波对信号进行边沿或突变点检测时,宜选用二进小波函数作为分析母小波；

(2) 对阶跃信号的分析宜选用零点反对称小波,对触发信号的分析宜选用零点对称小波；

(3) 阶跃信号的突变点分别对应零点反对称小波变换系数的极大值点和零点对称小波变换系数的过零点；触发信号的突变点分别对应零点反对称小波变换系数的过零点和零点对称小波变换系数的极大值点；

(4) Mallat小波在进行突变检测时会发生定位偏移,不利于突变点的定位；

(5) 正交小波基更适宜应用于数据的压缩和重构,不适宜突变点检测。

4 结 语

根据分析信号的要求,提出了基于样条函数构造不同性质的样条小波的理论和方法,并通过对实测故障行波数据的分析,得出了分析阶跃信号和突变信号的小波基函数选择的理论依据,证实其完全符合理论的推导。

参考文献

［1］易大义,陈道琦.数值分析引论［M］.杭州:浙江大学出版社,2002.

［2］崔锦泰.多元样条理论及应用［M］.程正兴,译.西安:西安交通大学出版社,1991.

［3］崔锦泰.小波分析导论［M］.程正兴,译.西安：西安交通大学出版社,1997.

［4］Stephane Mallat,Sifen Zhou.Characterization of Signals from Multiscale Edges［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence,1992,14(7):710-732.

［5］许传祥,石青云,程民德.零对称和反对称二进小波及其在边缘检测中的应用［J］.中国图像图形学报,1996,1(1): 4-11.

［6］Anna Wang,Zhang Xinhua,Chen Yu,et al.B Spline Wavelet and SVM Threshold Based Medical Image Edge Extraction.Industrial Electronics.2007.ISIE 2007.IEEE International Symposium on.2007:1 628-1 632.

［7］Daubechies Ten.Lectures on Wavelets.Philadelphia: Society for Industrial and Applied Mathe,1992:129-150.

［8］Liu L T,Hsu H T,Gao B X.A New Family of Orthonormal Wavelet Bases［J］.Journal of Geodesy,1998,72:294-303.

作者简介

周义明 男,1976年出生,讲师。主要研究方向为信号处理。

张 超 男,硕士研究生。主要研究方向为通信与信息系统。

张化民 男,1958年出生,副教授。主要研究方向为计算机应用。