基于小波變換的時變及典型非線性振動系統識別

黃東梅, 周 實, 任偉新,3

(1.中南大學 土木工程學院，長沙 410075；2. 高速鐵路建造技術國家工程實驗室，長沙 410075;3.合肥工業大學 土木與水利工程學院，合肥 230009)

非線性振動現象在控制工程、機械工程、土木工程等領域廣泛存在，非線性振動系統的識別是現代振動分析領域中研究的熱點和難點問題。總的來說，非線性振動系統的識別方法可以分為基于信號處理技術的方法、時間序列分析類方法、子空間類方法等[1-9]。近年來，時頻聯合分析被廣泛應用于線性系統模態參數識別方法中，同時也為弱非線性系統的辨識提供了新途徑。由于非線性振動的一個主要特征是系統的固有頻率和振幅有關，因此通過測量隨時間變化的信號頻率成分的變化情況可以實現對系統中非線性參數的識別。小波變換是一種信號的時間-尺度(與信號頻率相對應)分析方法，在時域和頻域同時具有表征信號局部特性的能力，非常適合進行時變和非線性振動系統的參數識別，因此，國內外對基于小波變換的非線性系統識別方法進行了一系列的研究：Staszewski[5]用Morlet小波變換對一個包含庫侖阻尼和三次方剛度的系統進行了辨識。Ta等[6]利用平均法推導了單自由度弱非線性自由振動系統的的瞬時頻率、解析幅度和解析幅度導數之間的理論函數關系，信號經Morlet小波變換后通過提取脊可以得到這三者的實際關系曲線，最后通過最小二乘曲線擬合即可估計系統的非線性阻尼和剛度系數，數值仿真表明該方法可以從弱非線性系統的自由振動響應中準確地辨識出庫侖阻尼、平方阻尼和三次方剛度。任宜春等[7]用Morlet復小波函數對弱Duffing系統的有阻尼自由振動響應進行了辨識，得到系統的固有頻率、阻尼系數和非線性系數。王超等[8]用連續小波變化的方法識別結構的非線性，基于Feldman提出的非線性結構骨架曲線概念，通過提取結構響應信號的小波脊和小波骨架，識別出結構的瞬時頻率和瞬時幅值，得到非線性結構的骨架曲線，從而識別結構的非線性。Kitad[9]對具有非線性恢復力系統的剪切型結構的方程進行離散，利用小波的時間局部性將切線剛度用離散小波尺度函數表示，阻尼系數分段線性化，使振動微分方程變為以阻尼系數和小波尺度函數為變量的代數方程組，采用最小二乘法計算得到未知向量，進而確定阻尼系數和切線剛度，進一步識別了結構的滯回曲線，最終通過對一個多自由度系統的試驗驗證了該方法的有效性。

1 時變及典型非線性振動系統及其解

1.1 時變阻尼自由振動系統

時變阻尼系統的自由振動方程可以寫為：

(1)

式中:ξ(t)為隨時間變化的阻尼比；ω0為無阻尼固有圓頻率。

在進行求解時，假定阻尼比在一個很短的時間內是慢變函數，則在這個短時間內可以看成是一個定常阻尼，此時方程(1)的解為

x(t)=A0e-ξω0tcos(ωdt+θ0)

(2)

式(1)中頻率是不隨時間變化的常數，瞬時振幅為

A(t)=A0e-ξω0t

(3)

可得

ln(A(t))=ln(A0)-ξω0t

(4)

1.2 達芬非線性振動系統

達芬有阻尼非線性自治方程為：

(5)

用Krylov-Bogoliubov平均法[10]進行求解，可得

(6)

瞬時振幅和瞬時頻率為

A(t)=A0e-ξω0t

(7a)

(7b)

達芬有阻尼非線性非自治方程為

(8)

由KBM法[10]，可得解為

(9)

式中的第一項為自由振動項，第二項為倍頻項，第三項為簡諧激勵項。其中

A(t)=A0e-ξω0t

(10a)

(10b)

式中:A0和θ0分別為式(9)第一項的初始振幅和初始相位，同式(6)。

對比達芬有阻尼非線性自由振動響應可知，其簡諧激勵響應第一項的瞬時振幅、瞬時頻率的表達式不變，而簡諧激勵振動響應為自由振動和簡諧振動的疊加。

1.3 范德波爾非線性振動系統

范德波爾非線性自治方程為：

(11)

用Krylov-Bogoliubov平均法[10]進行求解，可得

(12)

式中:A0為初始振幅；相位為常數θ0。當t=0時，若已知初始振幅為A0，則可以由式(12)得θ0=0。

瞬時振幅和瞬時頻率為

(13a)

ω(t)=ω0

(13b)

上式表明在一次近似精度范圍內的計算結果對頻率無修正，即相位不隨時間變化。

范德波爾非線性非自治方程

(14)

由KBM法[10]，可得解為

(15)

式中

(16a)

φ=ω0t-θ0

(16b)

式中:A0為初始振幅。當t=0時，若振幅為A0，則此時θ0=0。

對比范德波爾非線性自由振動響應可知，其簡諧激勵響應第一項的瞬時振幅、瞬時頻率的表達式不變，而簡諧激勵振動響應為自由振動和簡諧振動的疊加。

2 小波變換法

2.1 小波變換的基本理論

設信號x(t) ∈L2(R)，基本小波或母小波函數ψ(t) ∈L2(R)，L2(R)表示平方可積的實數空間。若ψ(t)滿足小波允許條件，則信號x(t)的連續小波變換定義為：

WTx(a,b)=x(t),ψab(t)=

(17)

小波變換的頻域表示為：

(18)

任意實信號x(t)可以表示成如下時變振幅與時變相位的乘積形式：

x(t)=A(t)cosφ(t)

(19)

然而，這樣的表示并不唯一，為保證表示的唯一性，通常引入其對應的解析信號：

(20)

(21)

對于單一分量信號x(t)，其相應的復Morlet小波變換根據平穩相位原理可以得到[13]：

WTx(a,b)=x(t),ψab(t)z(t),ψab(t)=

(22)

式中：O(?) 為二階小量，可以忽略不計；φ(b)為小波變換的相位。

故小波系數的模和相位分別為：

(23a)

∠(WTx(a,b))=φ(b)

(23b)

信號的瞬時頻率為：

(24)

小波系數的模反映了信號在時頻平面上的能量密度分布，而信號的能量在時頻平面上主要集中在脊線(每一時刻的頻譜峰值連成的曲線)附近，因此可以由任一時刻的信號小波系數模的脊線點所對應的頻域來確定瞬時頻率。

由式(23a)可以看出，當小波尺度a滿足條件：

(25)

小波系數模可取得局部極大值，相應的點(b,ar(b))稱為小波脊點。連接時間-尺度平面上的相應脊點便構成小波脊線。

此時，脊線上的小波系數的模為：

(26)

由此，可由脊線上小波系數模求出瞬時振幅：

(27)

對于多分量信號x(t)，通常可將其表示成單一分量信號的疊加：

(28)

因為小波變換是線性變換，相應x(t)的小波變換可表示為：

(29)

通過選擇合適的小波參數，可以將不同的單一信號分量在時間-尺度平面分離開來，從而分別提取各自小波脊，相應地得到信號瞬時頻率。

2.2 連續Morlet小波變換

Morlet小波是一種單頻復正弦調制高斯波，也是最常用的復值小波，其時、頻兩域都具有很好的局部性，它的時域、頻域形式如下[14]：

時域ψ(t)=e-t2/2ejωct

(30a)

(30b)

ωc為中心頻率，當ωc≥5 時，可認為Ψ(ω= 0) ≈0 近似滿足容許性條件。

對于復解析小波Morlet 小波，頻率f與尺度a之間的對應關系為：

(31)

其中:fs為信號的采樣頻率，*表示卷積運算[15]。

對于復Morlet 小波變換，其相應的時間分辨率和頻率分辨率分別為[16]：

(32a)

(32b)

時間和頻率分辨率滿足不確定原理，增大小波中心頻率ωc可以提高頻率分辨率，但同時會降低時間分辨率，增加端點效應的影響。

3 數值算例

3.1 時變阻尼自由振動系統

算例：

具體步驟如下所示：

(1)用四階龍格—庫塔法計算非線性系統的振動響應，如圖1所示，圖中還給出了上下包絡線。

圖1 振動響應圖

圖2 小波色譜圖

(2)用Morlet小波函數“cmor 2 - 1”對響應進行小波變換，得到色譜圖、三維網圖和等高線圖分別如圖2～4所示，圖2中還畫出了根據小波系數模的極大值確定的脊線。

(3)根據脊線可以得到瞬時頻率如圖5所示，為一常數，即ωd=44.44 rad·s-1，可近似取ω0=ωd。

(4)由脊線上的小波系數和式(27)可以求出瞬時振幅，如圖6所示，圖中還給出了由振動曲線圖1得到的振幅包絡線。另外，把瞬時振幅曲繪于圖1中予以對比。

(5)圖7給出了瞬時振幅的對數—ln(A(t))與時間t的關系曲線，剔除邊端效應，取0.4 s～1.6 s的數據，對曲線求導可得曲線ξ(t)ω0，進而求得阻尼比隨時間變化的曲線ξ(t)，如圖8所示，圖中還給出了擬合結果，ξ(t)=2.02×10-2t，以及理論結果ξ(t)=2×10-2t。

從以上分析可以看出，小波變換方法對時變阻尼自由振動系統的頻率和瞬時阻尼比識別是有效的。

圖5 瞬時圓頻率圖

圖8 瞬時阻尼比圖

3.2 達芬有阻尼非線性簡諧激勵振動系統

算例：

具體步驟如下所示：

(1)用四階龍格-庫塔法計算非線性系統的振動響應，如圖9所示，圖中還給出了根據KBM法的式(9)得到的振動響應，兩者是非常接近的。

(2)用Morlet小波函數“cmor 2 - 1”對響應進行小波變換，得到色譜圖、三維網圖和等高線圖分別如圖10～12所示，圖10中還畫出了根據小波系數模的極大值確定的脊線。從圖中可以看出，在1.2 s附近，簡諧激勵產生的振動能量開始大于初始位移引起的自由振動產生的能量，倍頻項很小，可以忽略不計。

(3)根據圖10的色譜圖可知，振動主要由兩個頻率的振動疊加而成，對比式(9)就可以知道這兩個頻率形成的原因，因此分別在[0～100] rad·s-1和[100～200] rad·s-1兩個頻率段提取脊線，便可以得到瞬時頻率如圖13(a)和圖13(b)所示，圖13(a)中還給出了根據式(7b)得到的理論結果，圖13(b)中還給出了ω的擬合結果和實際結果。從圖中可以看出識別結果為ω=138.28 rad·s-1,與理論結果140 rad·s-1非常接近。

圖9 振動響應圖

圖12 小波等高線圖

(4)由脊線上的小波系數和式(27)可以求出在[0～100] rad·s-1頻率段的瞬時振幅-即自由振動振幅，如圖14(a)所示，圖中還給出了根據式(7a)得到的理論結果。從圖中可以看出，剔除邊端效應，0.4 s～1.6 s的結果比較接近；由脊線上的小波系數和式(27)可以求出在[100～200] rad·s-1頻率段的瞬時振幅—即受迫振動振幅，如圖14(b)所示，圖14(b)中還給出了擬合結果為2.283和根據公式(9)第3項計算的的理論結果為2.276，兩者誤差很小。

(5) 圖15給出了瞬時振幅的對數-ln(A(t))與時間t的關系曲線，由式(7a)可得ln(A(t))=ln(A0)-ξω0，近似于線性關系，剔除邊端效應，取0.4 s～1.6 s的數據，由最小二乘擬合可以得到ξω0=2.287 4和A0=19.85 mm。

圖14 瞬時振幅圖

(6) 圖16給出了瞬時頻率ω(t)與瞬時振幅A(t)平方—A2(t)的關系，剔除邊端效應，近似于線性關系，取0.4 s～1.6 s的數據，由最小二乘擬合可以得到ω0=44.32 rad·s-1和ε=3.65，同時，根據步驟(5)的結果可以得到ξ=0.051 6。

圖16 瞬時圓頻率與瞬時振幅平方的關系圖

從以上分析可以看出，小波變換方法對達芬有阻尼非線性簡諧激勵振動系統的固有頻率和阻尼比、簡諧激勵頻率和幅值、自由振動初始振幅等的識別精度是較高的，而對非線性參數的ε識別誤差較大。

3.3 范德波爾非線性自由振動系統

算例：

具體步驟如下所示：

(1) 用四階龍格-庫塔法計算非線性系統的振動響應，如圖17所示，圖中還給出了根據平均法的式(13a)得到的振動響應，從圖中可以看出，由于平均法只有一次近似精度，所以存在微小的相位差。

(2) 用Morlet小波函數“cmor 2 - 1” 對響應進行小波變換，得到譜圖、三維網圖和等高線圖分別如圖18～20所示，圖18中還畫出了根據小波系數模的極大值確定的脊線。

(3) 根據脊線便可以得到瞬時頻率如圖21所示，圖中還給出了根據式(13b)得到的理論結果。已知頻率為常數，剔除邊端效應，取0.2 s～1.8 s的數據，直接進行最小二乘曲線擬合，便可以得到ω0=44.44(rad/s)，與理論值ω0=45(rad/s)非常接近。

(4)由脊線上的小波系數和式(27)可以求出瞬時振幅，如圖22所示，圖中還給出了根據式(13a)得到的理論結果。若已知公式(13a)，剔除邊端效應，取0.2s～1.8s的數據，初始值ε、δ和A0分別取0.4、15和30，直接進行最小二乘曲線擬合，可以得到ε、δ和A0的擬合值分別為0.045、2.18和20.38。擬合結果和理論結果如圖22所示。需要注意的是，由于用式(13a)進行擬合時要把3個參數當成未知數，因此擬合結果不唯一，與初始值的取值有很大的關系，并且剔除的數據不能太多，否則難以獲得足夠的信息。若把A0當作已知，則即使ε和δ的初始值取值不同，其擬合結果也較為穩定。

對于范德波爾非線性簡諧激勵振動系統的識別，可采用3.2節的類似的方法。

圖17 振動響應圖

圖20 小波等高線圖

4 結 論

本文用Morlet小波變換方法對時變振動系統和典型非線性振動系統進行參數識別，得出以下結論：

(1) 小波變換方法對時變阻尼自由振動系統的頻率和瞬時阻尼比識別是有效的，并且可以達到較高的精度。

(2) 由于小波變換是線性變換，因此用小波變換方法可以分別獲得達芬有阻尼非線性系統的自由振動和簡諧激勵振動，并對其分別識別。固有頻率和阻尼比、簡諧激勵頻率和幅值、自由振動初始振幅等的識別精度是較高的，而對非線性參數的ε識別誤差較大。

(3) 在用小波變換方法對范德波爾有阻尼非線性自治方程進行識別時，由于邊界效應的存在，會使得識別結果不夠穩定，在應用最小二乘法進行曲線擬合時，初始值的取值對識別結果的影響很大。

(4) 鑒于以上的分析，小波變換邊界效應的處理是下一步需要研究的方向，從而提高小波變換的識別精度，特別是對初值依賴性較強的情況。

[1] 于開平, 龐世偉, 趙婕. 時變線性/非線性結構參數識別及系統辨識方法研究進展[J]. 科學通報, 2009, 54: 3147-3156.

YU Kai-ping, PANG Shi-wei, ZHAO Jie. Advances in method of time-varying linear/nonlinear structural system identification and parameter estimate[J]. Chinese Sci Bull (Chinese Ver), 2009, 54: 3147-3156.

[2] 竇蘇廣, 葉敏氣, 張偉卞. 基于增量諧波平衡的參激系統非線性識別法[J]. 力學學報, 2010,42(2):332-336.

DOU Su-guang, YE Min-qi, ZHANG Wei-bian. Nonlinearity system identification method with parametric excitation based on the incremental harmonic balance method[J]. Chinese Journal of Theoretical and Applied Mechanics, 2010,42(2):332-336.

[3] 張根輩, 臧朝平.基于振動測試的非線性參數識別方法 [J].振動與沖擊,2013,32(1):83-89.

ZHANG Gen-bei, ZANG Chao-ping.A novel method for nonlinear parametric identification based on vibration tests[J].Journal of Vibration and Shock,2013,32(1):83-89.

[4] Jang T S. A method for simultaneous identification of the full nonlinear damping and the phase shift and amplitude of the external harmonic excitation in a forced nonlinear oscillator[J]. Computers and Structures, 2013, 120: 77-85.

[5] Staszewski W J. Identification of non-linear systems using multi-scale ridges and skeletons of the wavelet transform[J]. Journal of Sound and Vibration, 1998,214(4):639-658.

[6] Ta M N, Lardies J. Identification of weak nonlinearities on damping and stiffness by the continuous wavelet transform[J]. Journal of Sound and Vibration, 2006,293(1-2):16-37.

[7] 任宜春，易偉建. 非線性系統識別的小波方法研究[J]. 振動與沖擊，2007, 26 (3) : 68-71.

REN Yi-chun, YI Wei-jian. Identification of a nonlinear system using wavelet transformation[J]. Journal of Vibration and Shock,2007,26(3):68-71.

[8] 王超，任偉新，黃天立. 基于小波的非線性結構系統識別[J]. 振動與沖擊, 2009, 28 (3) : 10-14.

WANG Chao, REN Wei-xin, HUANG Tian-li. System identification of a non linear structure based on wavelet transformation[J]. Journal of Vibration and Shock, 2009,28(3):10-14.

[9] Kitada Y. Identification of nonlinear structural dynamic systems using wavelet[J]. Journal of Engineering Mechanics, 1998, 124: 1059-1066.

[10] 劉延柱，陳立群. 非線性振動[M]. 北京：高等教育出版社，2001.

[11] 陳逢時. 子波變換理論及其在信號處理中的應用[M]. 北京: 國防工業出版社, 1998.

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

[13] 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.

[14] 任偉新，韓建剛，孫增壽. 小波分析在土木工程結構中的應用[M]. 北京：中國鐵道出版社，2006

[15] Ruzzene M, Fasana A, Garibaldi, et al. Natural frequencies and dampings identification using wavelet transform: application to real data [J]. Mechanical Systems and Signal Processing, 1997, 11(2): 207-218.

[16] Kijewski T, Kareem A. Wavelet transforms for system 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.