李 书,陈 益
(北京航空航天大学飞机设计研究所,北京100083)
摘要:建立了超声探测缺陷回波的数学模型,讨论了信号奇异性同其小波变换之间的关系以及通过小波变换模极大值精确重构原信号的原理和方法,利用Mallat的交替投影算法对仿真的超声信号进行了精确重构和对实际检测到的超声信号进行了消噪处理。结果表明,利用小波变换模极大值重构信号的交替投影算法来重构超声信号,重构精度高,实现速度快,用于处理染噪信号,消噪效果好,是一种较为理想的处理超声信号的方法。
关键词:小波变换;模极大值;超声信号;交替投影算法;奇异性
中图分类号: TG115.28 文献标识码: A 文章编号: 1004-4523(2006)02-0206-06
引 言
超声检测是目前国内外最为实用有效、应用最为广泛的无损检测技术,人们更多的是利用现代数字信号处理技术来对超声信号进行处理,而小波变换正是一种理想的处理超声信号这种时变非稳态脉冲信号的方法[1]。
在信号处理中,人们总希望能提取反映信号特征的信息,而信号的剧烈变化点(奇异点)处通常携带有反映信号特征的重要信息,所以人们关注能否用信号的奇异点表征和重构原信号。最初Logan证明了窄带信号的过零点可以完整地表征原信号[2],但重构不稳健,难以实际应用。此后有很多学者都在这方面展开了研究工作[3～5]。而信号小波变换的奇异点(极值点和过零点)和信号的剧烈变化处有着密切的联系。由信号二进小波变换的奇异点在多尺度上的表示来表征和重构信号,是小波变换的一个重要应用领域[6～9]。Mallat采用信号二进小波变换模极大值来表征信号,并且提出了用信号二进小波变换模极大值重构信号的交替投影算法[7],实验证明,该算法具有良好的逼近特性,能精确地重构原信号。
由于其在信号的数据压缩和信号去噪等方面的良好应用前景,也引起了工程界的广泛关注[9～11]。本文建立了超声回波信号的数学模型,讨论了信号奇异性同其小波变换之间的关系以及通过小波变换模极大值精确重构原信号的原理和方法,应用Mallat的交替投影算法对超声信号进行了精确的重构,并利用交替投影算法对两个实际的超声检测信号进行了消噪处理,消噪效果良好。
1 超声回波数学建模
在宽带超声检测中,超声回波信号通常是一个被探头中心频率调制的宽带信号,超声缺陷回波的数学模型可建立如下[12]
式中 f0为探头的中心频率;B0确定f(t)的带宽;超声脉冲信号的功率谱通常被建模为Gaussian函数[13],考虑包络h(t)为Gaussian函数时,式(1)变为
即超声缺陷回波的数学模型是Gaussian函数经过调制得到的。考虑μ(f)exp(iθ(f))为噪声n(t)的频域建模N(f)时,则系统接受到的信号频域表达式为
2 小波变换检测信号奇异性
信号中所包含的信息主要体现在信号的瞬变点或瞬变的区域中,信号的瞬变程度常用信号的奇异性来描述,而Lipschitz指数(李氏指数)则是描述信号奇异性的重要指数,其定义如下
对于任意给定信号f(t),若存在常数K>0及n=[α]阶的多项式pt0(t),使得对于任意t
称f(t)在t0处具有李氏指数α。
李氏指数α和信号f(t)在t0或在区间[a,b]上的可微性有关,用来衡量信号f(t)在此处的可微程度。
若f(t)在此处的导数阶次越高,相应的α越大,反映在信号的特征上,f(t)在此处越平滑;若f(t)在此处的李氏指数小于1,则信号在此处是不可微的,也就是奇异的。
设ψ(t)∈L2(R)(L2(R)表示平方可积的实数空间,即能量有限的信号空间),其傅立叶变换为当满足相容性条件
时,称ψ(t)为一个基本小波或母小波。将母函数ψ(t)经过伸缩和平移后,就可以得到一个小波序列。对于连续的情况,小波序列为
式中 a为伸缩因子,b为平移因子。
用 表示卷积操作,则对于任意的函数f(t)∈L2(R),它的连续小波变换为
假定θ(t)是一个低通函数,且高阶可微,令 ,则
所以,是带通函数,它们可以作为小波母函数使用,定义
则对于任意的函数f(t)∈L2(R),有
式(10),(11)两式说明对信号平滑(即通过一个低通滤波器)后求一阶或二阶导数等效于直接用该平滑函数(即低通函数θ(t))的一阶或二阶导数来滤波。此等效过程还可以推广到θ(t)的更高阶导数,如信号f(t)通过θ(t)后求p阶导数等效于直接让f(t)通过
从数学的角度,一个函数的一阶导数等于零的点对应于该函数的极值点,而二阶导数等于零的点对应于该函数的拐点,即转折点。这样,如果用于小波变换的小波函数是来自于某一个低通函数的一阶导数或二阶导数,那么小波变换的结果将体现出信号的极值点或转折点。信号中常见的瞬变有两种,一是边缘的突变,这相当于在该处叠加了一个阶跃信号;另一个是峰值的突变,这相当于在该处叠加了一个冲激信号。这两种情况分别对应了信号的极值点和转折点,它们统称为信号的奇异点。显然,这些奇异点一般都可以在其小波变换的幅值中反映出来,即或是对应小波变换的过零点,或是对应小波变换的峰值点。具体的说,用ψ(1)(t)对f(t)做小波变换,其等于零的点反映了f(t)的极值,因此可以实现极值点的检测;用ψ(2)(t)对f(t)做小波变换,其等于零的点反映了f(t)的转折点,从而可实现转折点的检测。
3 小波变换模极大值与李氏指数
设f(t)∈L2(R),并且f(t)在t0处具有李氏指数α≤n,n∈Z+为小波ψ(t)的消失矩的阶数,则存在常数A使得[7,14,15]
成立。考虑f(t)在t0处有一个奇异点的情况,很容易想象,f(t)在t0处的奇异点将不会影响到整个尺度时间平面上的小波变换,而主要影响该平面上围绕t0的一个小的范围。该范围成为t0的影响锥。假定所使用的小波ψ(t)具有紧支撑,支撑范围是[C,C],那么ψa,b(t)的支撑范围是[t-Ca,t+Ca],所谓影响锥,是指尺度时间平面上使得t0包含在ψa,b(t)范围内所有点的集合。所以,t0的影响锥为
将式(13)带入式(12)可得
在二进制小波变换中,令a=2j,对式(14)两边取以2为底的对数,有
由上式可知,对于任意的信号f(t),有
1)若t0处的李氏指数α>0,那么小波变换的模极大值随着尺度j的增大而增大;
2)若t0处的李氏指数α<0,那么小波变换的模极大值随着尺度j的增大而减小;
3)若t0处的李氏指数α=0,那么小波变换的模极大值不随尺度j的变化而变化。
4 小波变换模极大值重构信号
设信号f(t)∈L2(R),其二进小波变换为二进小波变换在处取得局部模极大值 为了通过小波变换模极大值重构原信号f(t),假定有一个信号集合h(t),该集合中信号的小波变换和f(t)的小波变换具有相同的模极大值。显然,希望在某一准则下在h(t)中选取一个信号来最佳地逼近f(t)。记h(t)的小波变换为则h(t)必须满足以下两个条件:
1） 对应每一个尺度j,在所有的模极大值横坐标{tjn}n∈Z处,都应有
2） 对应每一个尺度j,的局部模极大值都应位于模极大值横坐标处。
对于条件1,由于对于任意的点t0有
所以条件1等价于
设U为函数族所张成的L2(R)的闭子空间,所以函数h(t)满足式(18)的充分必要条件是h(t)在U上的正交投影等于f(t)在U上的正交投影。
记O为U在L2(R)中的正交补空间,即
取g(t)∈O,则有
由此可见,只有O={0}时,才有h(t)=f(t),而在一般情况下,O≠{0},所以条件1下h(t)并不能唯一确定f(t)。
对于条件2,这是一个很难实现的要求。要求的局部极值点都位于模极大值横坐标处,Mallat提出了一个近似的方法,其思路是:条件1确定了处取模极大值,但不强求在其他的横坐标处没有模极大值,把条件2修改为要求 在其他点上的平均值尽可能为最小模极大值点数的多少取决于其振荡情况,为了使以外有尽可能少的模极大值点,除了要求最小外,还需要求{W2jh(t)}j∈Z导数的能量也最小。为此,引入Sobolev范数[16]
该范数将两个最小都包括了进去。通过式(21)求解出使得 ‖h‖ 2为最小的W2jh(t),就可以得到f(t)的最佳逼近h(t),采用交替投影方法求解。
5 信号重构的交替投影算法
令V是L2(R)上所有信号的二进小波变换所组成的空间,令K是序列gj(t)所组成的空间,gj(t)满足
显然,对比式(21)可知,序列gj(t)应是某一个序列的小波变换。令Γ是空间K上的一个闭包,其中的元素gj(t)满足
这样,满足条件1的小波变换应是空间
中的元素。则现在只需求Λ中的一个元素,使其Sobolev范数为最小即可。采用V,Γ之间的交替投影可以实现这一目标。
显然Γ中的元素即为前述的h(t),求h(t)的过程,即是令式(22)的Sobolev范数为最小。记投影算子
PV可以将K中的任一序列投影到空间V。
投影算子PΓ的作用是将序列X投影到空间Γ。
由于X∈K,所以在其中的序列都定义了Sobolev范数。经PΓ投影后,X的极值点及其大小将和 W2jf(tjn) 相同,并最接近hj(t)。所谓最接近,即是令
并且使其在Sobolev范数意义上为最小,即
为最小,此时,应该保证在各个尺度j,式(27)中的每一项都为最小,即
设的两个相邻的模极大值点,由于,因而
因此,在区间[tn,tn+1]上,式(28)最小等价于
为最小。式中εj(t)满足式(29)的约束条件。
式(26)～(30)构成了一个有约束的最小化问题,对式(30)可通过求解如下的Euler方程给出
在区间[tn,tn+1]上,式(31)的解为
式中 常数α和β由式(29)的约束条件给出。
实现PΓ时,先通过式(32)和(29)求出εj(t),再由式(26)得到
令gj(t)=0开始迭代,则交替投影可求得hj(t)。
令为空间V,Γ之间的交替投影算子,为P是的n次迭代,则对任一序列可以证明有
因此,空间V,Γ之间的交替投影收敛到对空间Λ的正交投影。迭代开始时,选gj(t)=0,则交替投影后收敛到Λ中的元素的范数接近于0,做到了使Sobolev范数最小。空间V,Γ之间的交替投影如图1所示。
6 仿真实验结果分析
仿真实验采样频率为50 MHz,采样点数为1 024个,小波函数采用sym8小波,分解尺度为6。本次实验仿真了两个具有不同中心频率和幅值的超声缺陷回波信号,并根据本文所述算法重构原信号,其结果如图所示。
图2给出了仿真的超声回波原信号,图3,图4分别是小波变换的近似系数和细节系数,图5从原信号的小波变换中提取出了各分解尺度下的模极大值,图6给出了15次迭代和20次迭代后重构出的信号。由图6可以看出,15次迭代后,重构信号信噪比为35. 142, 20次迭代后,重构信号信噪比为37.101 5,该方法重构超声信号精度非常高。15次迭代用时3.58 s,20次迭代用时4.27 s,可以看出,该方法用来重构超声信号运算量不大,实现方便。
图7和图8是进行超声探伤时超声探伤仪所接收到的两个超声回波信号波形(为了说明问题,特意在接收到的回波信号中加上了高斯白噪声)。探伤采用自发自收的A扫描方式,用超声探伤仪激发并接收中心频率为2 MHz的宽带、窄脉冲超声信号,该信号经采样频率为10 MHz的信号采集板数字化后存于微机,采样点数4 096个。图7的回波信号中含有一个缺陷回波,它是由直径为1 mm、深25 mm的圆柱形横通孔产生的;图8的回波信号中不含任何缺陷回波。从图7和图8来看,缺陷回波完全淹没于强噪声中,很难判断缺陷的存在与否。
图9和图10是采用交替投影算法对原信号进行消噪处理后的信号波形,图中最左边的回波为发射波,最右边的回波为底面回波,中间的回波为缺陷回波。由经消噪处理后的波形可以看出,缺陷的存在与否一目了然,采用交替投影算法对超声回波信号进行处理,消噪效果较为理想。
7 结 论
超声回波信号的处理是超声检测的难点。小波变换以其特有的“变焦距”特性,非常适于用来分析超声信号。在使用小波变换分析超声信号时,利用小波变换的奇异点在多尺度上的表示来表征和重构信号,是小波变换的一个重要应用领域,在信号的数据压缩和信号消噪等方面有着良好的应用前景。本文讨论了利用小波变换模极大值精确重构原信号的交替投影算法,并对超声信号进行了仿真实验。结果表明,这种方法用来重构超声信号,重构精度高,实现速度快,用于处理染噪信号,消噪效果好。而且这种方法思路清晰,原理简单,实现方便,在工程上有着很高的应用价值。
参考文献:
[1] 崔锦泰.小波分析导论[M].西安:西安交通大学出版社,1995.
[2] Logan B. Information in the zero-crossings of bandpass signals[J].Bell System Tech J.,1977,(56):510.
[3] Curtis S, Oppenheim A. Reconstruction of multidi-mensional signals from zero-crossings[J].J Opt SocAmer,1987,(4):221.
[4] Zeevi Y Y,Rotem D.Image reconstruction from zero-crossings[J].IEEE Trans ASSP,1986,(34):1 269.
[5] Sanz J,Huang T.Image representation by sign infor-mation[J].IEEE Trans PAMI,1992,(11):729.
[6] Mallat S.Zero-crossings of a wavelet transform[J].IEEE Trans Information Theory,1991,37(4):1 019—1 033.
[7] Mallat S, Zhang S. Characterization of signal frommultiscale edges [J]. IEEE Trans Pattern Analysisand Machine Intelligence,1992,14(7):710—732.
[8] Liu G Z,Zhang Z M.Iterative shaping reconstructionalgorithm based on the modulus maxima of signals'dyadic wavelet transform[J]. Progress in Natural Sci-ence,2000,10(8):660—664.
[9] Cvetkovic Z, Vetterli M. Discrete-time wavelet ex-treme representation: Design and consistent recon-struction[J].IEEE Trans Signal Processing,1995,43(3):681—693.
[10] Berman Z,Baras J.Properties of the multiscale maximaand zero-crossing representation[J].IEEE Trans SP,1993,41(12):3 216.
[11] Cetin A, Ansari R. Signal recovery from wavelettransform maxima[J].IEEE Trans SP,1994,42(1):194.
[12] Gustafsson M G, Stepinski T. Split spectrum algo-rithms rely on instantaneous phase information-a geo-metrical approach [J]. IEEE Trans UFFC, 1993, 40(6):659—665.
[13]张广明.超声无损检测中的时频分析理论及应用研究[D].西安:西安交通大学,1999
[14] Mallat S,Hwang W L.Singularity detection and pro-cessing with wavelets[J].TEEE Trans Inf Theory,1992,38(2):617—643.
[15] Mallat S.A Wavelet Tour of Signal Processing[M].San Diego,CA:Academic Press,1997
[16] Mazja V G.Sobolev sPACes[M].Springer-Verlag,1985
