基于动态规划提取信号小波脊和瞬时频率

王  超，任伟新

(中南大学 土木建筑学院，湖南 长沙，410075)

摘  要：提出一种基于动态规划提取信号小波脊和瞬时频率的方法，其基本思路是：对信号进行连续复Morlet小波变换，由变换得到的小波系数的局部模极大值初步提取其小波脊；为降低噪音影响，在初步提取的各小波脊附近选取部分小波系数，通过施加罚函数平滑噪音干扰引起的小波脊变化的不连续性，将小波脊的提取问题转变为最优化问题，采用动态规划方法计算得到新的小波脊；根据小波尺度与频率的关系由提取的小波脊识别出信号的瞬时频率。将提出的方法运用于含噪调频信号进行数值模拟分析和实测索冲击响应信号分析。研究结果表明，基于连续小波变化的模极大值可以有效提取信号小波脊和瞬时频率；采用施加罚函数的方法可有效降低噪音的影响；基于动态规划的方法可有效提高计算效率。

关键词：小波变换；小波脊；瞬时频率；动态规划

中图分类号：TN911.6         文献标识码：A         文章编号：1672-7207(2008)06-1331-06

Wavelet ridge and instantaneous frequency extraction based on dynamic optimization

WANG Chao, REN Wei-xin

(School of Civil and Architectural Engineering, Central South University, Changsha 410075, China)

Abstract: A new method of extracting signal wavelet ridge and instantaneous frequency based on dynamic optimization was presented. Its basic ideal was as follows: First, wavelet ridges were extracted initially based on maximum value of modulus of wavelet coefficients, then some coefficients near the initially wavelet ridges were chosen and penalty function was introduced. So the problem was transformed to an optimization problem. Finally, new wavelet ridges were extracted by using dynamic optimization method. According to the relation between wavelet scale and frequency of signal, the instantaneous frequency of signal was obtained. A frequency modulated signal with noise and impulse response signal of cable tested in laboratory were simulated and analyzed using the method. The results show that the present algorithm can be used to extract wavelet ridges and instantaneous frequency of signal. It can decrease influence of noise validly. In addition, the algorithm has high efficiency in computing.

Key words: wavelet transform; wavelet ridge; instantaneous frequency; dynamic optimization

信号的特征分析与提取是信号分析与处理中的一类重要问题，其应用范围非常广泛，例如雷达声纳探测、语音信号处理、机械故障诊断、结构健康监测和损伤识别等^[1-5]。目前，已有多种方法对信号进行时频分析，如双线性时频分布(Wigner-Vill分布以及其统一形式—— 广义双线性时频分布，即能量化Cohen类时频分布)和线性时频分布如Gabor变换和小波变换^[6-7]。

1992年，Delprat等^[8]提出一种新的算法提取信号特征。其方法是基于Gabor变换或小波变换后的相位信息，在时间尺度空间通过迭代的方法提取小波脊线，从而可以获得信号的时频信息。但是，噪声会对小波变换的相位产生较大影响。Carmona等^[9]提出了施加罚函数的算法来提取小波脊，它基于小波变换系数的模信息。小波系数的模通常具有比相位更强的抗噪性，然而，这种方法只适用于单一分量的信号。对于多分量信号，Carmona等^[10]提出基于随机走动的爬山算法，该算法较复杂且计算较耗时。Helene等^[11-12]提出了更简单的方法提取小波脊，他们基于小波系数模的局部极大值。Haase等^[13]提出了一个新的根据偏微分方程直接积分的方法来提取小波变换的小波脊。

总的来说，主要采用2种方法提取小波脊：一种是基于小波系数的相位信息，一种是基于模信息。这里提出一种简单方法提取多分量信号小波脊。在初步提取的小波脊线附近选取部分小波系数，通过引入Carmona等^[9]提出的罚函数，将小波脊线提取问题转换为最优化问题，用动态规划的方法重新提取各小波脊线，以降低噪音的影响，提高计算效率。

1  连续小波变换与小波脊线

连续小波变换是一种线性变换，为信号与小波基函数的内积。小波母函数是满足允许性条件的平方可积函数，小波基函数则通过将小波母函数进行膨胀和平移得到，定义为：

信号x(t)的连续小波变换定义为：

这里采用复Morlet小波进行分析，其定义为：

然而，这样的表示并不唯一，为保证表示的唯一性，通常引入其对应的解析信号：

，为信号的Hilbert变换。假定  A_z(t)≥0和，上式表示将唯一，由此引入瞬时频率定义^[6]：

对于式(1)定义的单一分量信号x(t)，其相应的复Morlet小波变换根据平稳相位原理可以得到^[8]：

由式(8)可以看出，当小波尺度a满足条件：

             (9)

时，小波系数模取得局部极大值，相应的点(b, a_r(b))称为小波脊点。连接时间-尺度平面上相应脊点构成小波脊线。沿着小波脊线，信号能量(与小波变换系数模的平方对应)主要集中在时间-尺度平面上脊线附近。提取出小波脊，由尺度与频率的关系式(4)可以识别出信号瞬时频率。

对于多分量信号x(t)，通常可将其表示成单一分量信号的叠加：

通过选择合适的小波参数，可以将不同的单一信号分量在时间-尺度平面分离开来，从而分别提取各自小波脊，相应地得到信号瞬时频率。

2  小波参数的选取

对于复Morlet小波变换，其相应的时间分辨率和频率分辨率分别为^[14]：

时间和频率分辨率满足不确定原理，增大小波中心频率f_c可以提高频率分辨率，但同时会降低时间分辨率，增加端点效应的影响。

为了获得最佳的时频分辨率，LIN^[15]提出利用小波熵来确定小波参数。假定小波系数为WT_x(a_i, b), i=1, 2,…, m，可以得到：

小波熵定义为：

小波熵反映了小波系数的稀疏性，最优小波系数应该具有最小的小波熵。因此，可以根据小波熵来选择合适的小波基。

3  基于动态规划提取小波脊

先根据小波系数模信息初步提取小波脊，然后，在其附近频率范围内选取相应的小波系数，引入Carmona等^[9]提出的罚函数来降低噪声的影响，将小波脊提取问题转化为一个最优化问题，用动态规划方法重新提取小波脊线。

假定小波脊线为参数曲线(b, )，相应的罚函数可以表示为：

其中：用于平衡小波变换系数幅值和尺度参数变化梯度的正常数。参数的选择依赖于小波变换系数。当取0时，可直接根据小波系数模极大值来提取小波脊。因为噪声使局部极大值发生偏离，导致频率发生跳跃而使尺度变化梯度增加，从而使罚函数值增加。增大可以调节由于噪声引起的极大值和真实极大值偏离而导致的小波脊变化不连续，使提取的小波脊线变光滑。因此，当小波脊线取小波系数模极大值且变化曲线比较平滑时罚函数应具有较小值，可以通过寻找使罚函数取最小值的参数曲线来提取小波脊线。

在离散情况下，罚函数可以表示为：

其中：N为信号长度；k为小波离散平移参数；为离散尺度参数。求罚函数最小值可通过动态规划法求全局最优路径来求得。

假定前i步最优子路径已经确定，终点通过点(, i)(其中，)，其相应的罚函数为F(m, i)，下一步最优路径通过点(，i+1)(其中，)，则相应的罚函数应为：

要确定最优路径i+1步对应的点，需要将这一点所有可取候选脊点与前i步最优子路径组成的i+1步路径进行比较，具有最小罚函数值的点保留下来，从而得到前i+1步最优子路径。

用同样方法继续寻找下一步最优路径至信号全长，从而可以提取出全局最优路径，得到小波脊线。具体算法如下：

a. 计算小波系数，画出小波量图，从小波量图中根据局部模极大值在每条脊上初步找到1个初始脊点(b₀ , a₀)位置(如图1中圆圈点)。

b. 对每个初始脊点，在其局部尺度范围(a₀－d_a，a₀ +d_a)内寻找下一个脊点(b₁ , a₁)。然后，将新的脊点(b₁ , a₁)作为初始脊点继续提取下一个脊点至信号两端，由此初步得到小波脊线(如图1中实线所示)。

实线—初始小波脊；虚线—小波系数选取范围

图1  信号小波量图和小波脊线

Fig.1  Signal wavelet scalogram and wavelet ridge

c. 在初步提取的各小波脊线附近尺度范围内分别提取相应的小波系数(见图1中虚线)作为下一步提取小波脊的候选系数。

d. 施加罚函数，用动态规划的方法，通过求罚函数最小值来重新提取出各小波脊线。

4  数值算例

如下式所示的调频信号：

信号采样频率为40 Hz，共采样1 024个点，对信号施加均值为零的高斯白噪声得到含噪信号y(t)，信噪比为11 dB(定义为)。其中：A为信号幅值；为噪声信号方差)。选用复Morlet小波，小波中心频率根据小波熵选取。模拟信号小波熵与中心频率曲线关系见图2。由图2可选取中心频率为3 Hz。

图2  小波熵与中心频率关系

Fig.2  Relationship between wavelet entropy and central frequency curve

对信号进行连续小波变换，计算其相应的小波量图，如图3所示。由图3可以看出，信号存在2条小波脊线，用本文方法对小波脊线进行提取，根据模极大值初步提取结果如图4所示(尺度已由式(4)转化为信号瞬时频率)。由于噪音的影响，小波脊线存在许多局部极值，变化不光滑。在初步提取小波脊附件选取部分小波系数，施加罚函数，用动态规划算法重新提取小波脊，识别的瞬时频率结果见图5。

图3  信号小波量图

Fig.3  Signal wavelet scalogram

1—第1分量；2—第2分量

图4  直接根据模极大值提取的瞬时频率

Fig.4  Extracted instantaneous frequency based on maximum value of modulus

1—第1分量(实线为提取值，虚线为理论值)；2—第2分量(实线为提取值，虚线为理论值)

图5  动态规划法提取的瞬时频率

Fig.5  Extracted instantaneous frequency based on dynamic optimization

从图5可以看出，通过施加罚函数，降低了噪音的影响，所提取的小波脊线(信号瞬时频率)变得光滑，更符合实际结果。为与信号瞬时频率理论值进行比较，在图中同时画出了理论瞬时频率变化曲线(图中虚线所示)。信号2个分量理论瞬时频率为：

考虑识别瞬时频率与理论结果的差别，其对应的误差百分比曲线如图6所示。可以看到，第1个分量识别结果主要在频率变化极值点附近相对误差稍大(图6(a))，最大误差在4. 5 s处为2.2%，其他时刻相对误差基本小于2%；第2个分量识别结果在频率突变点和端点附件相对误差较大(图6(b))，突变点处与小波时间分辨率有关，端点处的误差主要是由于小波变换的端点效应的影响所致。其他大部分时刻相对误差也小于2%。

(a) 第1分量相对误差；(b) 第2分量相对误差

图6  瞬时频率误差

Fig.6  Errors of instantaneous frequency

5  实验验证

将所提出的方法应用于实验结构自由响应信号，实验结构如图7所示。采用1股7φ5预应力钢绞线拉索，拉索一端锚固，一端施加变化拉力，拉力从t =    0 s 时的34 kN线性变化到t = 25 s 时的67 kN。拉力的变化将改变索的刚度，从而导致其瞬时频率的变化。用力锤对拉索施加冲击荷载，采集其加速度响应，采样频率为1.2 kHz。主要考虑结构的基频，该结构基频小于50 Hz。因此，以150 Hz采样频率对信号进行低通滤波和重采样，重采样的信号如图8所示。

图7  实验模型图

Fig.7  Experimental model

图8  冲击响应信号

Fig.8  Impulse respond signal

对于索结构，结构瞬时自振频率可以直接根据拉力计算。最简单的计算公式可表示为：

式中：T 为索的拉力；＝1.1 kg/m，为拉索的线密度；L＝3.84 m，为索的长度。因为实际边界条件不确定及公式存在近似性，拉力与频率的关系与式(18)有差别，这里通过实测固有频率与拉力的关系对式(18)进行修正，修正结果为：

运用提出的方法识别结构的瞬时频率，结果如图9所示。

由图9可见，识别的结果与直接通过索拉力计算的结果非常接近。在终点时刻相对误差较大，主要是由于冲击响应到最后信号振幅衰减比较大，能量比较低，因而信噪比小，噪音的影响比较大。

实线—瞬时频率；虚线—用拉力计算所得频率

图9  识别的瞬时频率和用拉力计算的频率

Fig.9  Identified instantaneous frequency and results calculated from tension forces

6  结  论

a. 提出了一种根据连续复小波变换小波系数的模极大值来提取信号小波脊线和瞬时频率的方法。直接根据模极大值提取小波脊线易受噪音干扰而产生波动，为降低噪音影响，采用施加罚函数的方法来平滑噪音干扰引起的小波脊变化的不连续，将小波脊的提取问题转变为最优化问题，采用动态规划的方法来提取小波脊线。

b. 该方法具有较强的抗噪性和较高的计算效率。该方法不必事先假定信号瞬时频率的变化规律，具有自适应性。此外，该方法识别结构瞬时频率时只需对结构响应信号进行处理，不需要结构的输入信号，便于实际应用。

参考文献：

[1] 任伟新, 韩建刚, 孙增寿. 小波分析在土木工程结构中的应用[M]. 北京: 中国铁道出版社, 2006.
REN Wei-xin, HAN Jian-gang, SUN Zeng-shou. Application of wavelet analysis in civil engineering structure[M]. Beijing: China Railway Press, 2006.

[2] 孙增寿, 韩建刚, 任伟新. 基于小波分析的结构损伤检测研究进展[J]. 地震工程与工程振动, 2005, 25(2): 93-99.
SUN Zeng-shou, HAN Jian-gang, REN Wei-xin. State-of-the-art research and development of wavelete analysis based structural damage detection[J]. Earthquake Engineer and Engineering Vibration, 2005, 25(2): 93-99.

[3] 吕 琛, 王桂增, 叶 昊, 等. 基于噪声小波包络谱的主轴承磨损故障诊断[J]. 中南工业大学学报: 自然科学版, 2003, 34(4): 459-462.
L? Chen, WANG Gui-zeng, YE Hao, et al. Abrasion fault diagnosis of main bearing based on noise and wavelet envolope spectrum[J]. Journal of Central South University of Technology: Natural Science, 2003, 34(4): 459-462.

[4] 李夕兵, 张义平, 左宇军, 等. 岩石爆破振动信号的EMD 滤波与消噪[J]. 中南大学学报: 自然科学版, 2006, 37(1): 150-154.
LI Xi-bing, ZHANG Yi-ping, ZUO Yu-jun. Filtering and denoising of rock blasting vibration signal with EMD[J]. Journal of Central South University: Science and Technology, 2006, 37(1): 150-154.

[5] WANG Chao, REN Wei-xin. A wavelet ridge and SVD-based approach for instantaneous frequency identification of structure[C]//Proceedings of the Second International Conference on Structural Condition Assessment, Monitoring and Improvement. Beijing: Science Press, 2007: 803-807.

[6] Cohen L. Time-frequency analysis[M]. New Jersey: Prentice Hall, 1995.

[7] 张 晔. 信号时频分析及应用[M]. 哈尔滨: 哈尔滨工业大学出版社, 2006.
ZHANG Ye. Signal time-frequency analysis and application[M]. Harbin: Harbin Institute of Technology Press, 2006.

[8] Delprat N, Escudie B, Guillemain P, et al. Asymptotic wavelet and gabor analysis: extraction of instantaneous frequencies[J]. IEEE Transactions on Information Theory, 1992, 38(2): 644-644.

[9] Carmona R, Wen L H, Torrdsani B. Characterization of signals by the ridges of their wavelet transforms[J]. Transactions on Signal Processing, 1997, 45(10): 2586-2590.

[10] Carmona R, Wen L H, Torrdsani B. Multi-ridge detection and time-frequency reconstruction[R]. California: University of California at Irvine, 1995.

[11] Helene L, Catherine M. Ridge extraction from the scalogram of the uterine electromyogram[C]//Proceedings of the IEEE-SP International Symposium. Pittsburgh, 1998: 245-248.

[12] Heng L, Alexander N. Moire interferogram phase extraction: a ridge detection algorithm for continuous wavelet transforms[J]. Applied Optics, 2004, 43(4): 850-857.

[13] Haase M, Widjajakusuma J. Damage identification based on ridges and maxima lines of the wavelet transform[J]. International Journal of Engineering Science, 2003, 41: 1423-1443.

[14] Kijewski T, Kareem A. Wavelet transforms for system identification in civil engineering[J]. Journal of Computer-Aided Civil and Infrastructure Engineering, 2003, 18(5): 339-355.

[15] LIN Jing. Feature extraction based on Morlet wavelet and its application for mechanical fault diagnosis[J]. Journal of Sound and Vibration, 2000, 234(1): 135-148.

收稿日期：2008-01-08；修回日期：2008-03-25

基金项目：国家自然科学基金资助项目(50378021)

通信作者：任伟新(1960-)，男，湖南长沙人，博士生导师，教授，从事结构健康监测、损伤诊断、稳定与振动研究；电话：0731-2654349；E-mail: renwx@mail.csu.edu.cn