李蘆鈺,牛 蕓
(大連理工大學 土木工程學院,遼寧 大連 116023)
基于小波變換對非線性系統進行參數識別已經成為土木工程界研究的熱點問題之一。小波分析是一種信號的時間-頻率分析方法,它具有多分辨率分析的特點,而且在時、頻兩域都具有表征信號局部特征的能力,很適合分析時變信號,國內外學者對小波分析方法在系統識別中的應用進行了大量的研究。Staszewski[1]以Morlet小波作為分析工具,通過提取小波脊初步建立了一種非線性阻尼及剛度的識別方法,但它僅能對系統的非線性進行定性分析。Ta等[2]完善了這種方法,進行了仿真驗證并將其用于辨識一個固支梁的非線性阻尼及剛度。任宜春[3]用復Morlet小波函數對弱Duffing系統的有阻尼自由振動響應進行了辨識,得到系統的固有頻率、阻尼系數和非線性系數。伊廷華等[4]針對小波變換中遇到的邊端效應問題,提出了基于自回歸滑動平均模型(ARMA)的“預測延拓”方法,并以美國土木工程師學會(ASCE)提供的Benchmark模型為例進行了數值模擬,驗證了方法的有效性,但該方法只能對平穩信號進行預測。代煜等[5]利用脊上連續小波變換系數的幅度和相位,從結構系統的自由衰減響應中辨識了弱非線性阻尼和剛度,提出了一種檢測小波脊的新方法以消除連續小波變換幅度極值的頻移,為抑制邊界效應,提出利用最小二乘擬合誤差選擇小波函數參數的優化算法。王超等[6]采用復Morlet小波對非線性結構自由響應信號進行連續小波變換,為了降低噪聲的影響,采用一個基于奇異值分解(SVD)的方法對識別的結構進行處理,從而識別了待辨識的參數。史治宇等[7]首先運用Daubechies小波識別了時變系統的物理參數,文中借助小波尺度函數的正交性將物理空間下的二階微分方程轉化為線性代數方程組進行求解,識別算法中需要同時計算一階和二階小波連接系數,難度大,精度受限。許鑫等[8]用小波狀態空間法實現了從振動微分方程到線性代數方程組的兩次降階,解決了時變系統在自由振動狀態下的參數識別問題,許鑫等[9]用此方法解決了時變系統受迫振動下的參數識別問題。
在利用小波分析方法進行參數識別的過程中,不可避免的會出現邊端效應問題,進而會影響識別結果的準確性。在取較少采樣點的情況下,同樣也會影響參數識別結果的準確性。信號延拓是處理這兩個問題的有效方法之一,Torrence等[10]采用零延拓,未考慮信號的原有特征,效果較差;Kijewski等[11]采用對稱延拓,但該方法沒有考慮信號的連續性,存在一定的局限性;伊廷華等[4]提出了基于自回歸滑動平均模型(ARMA)的“預測延拓”方法,但只對平穩信號有效。相關研究指出,不能任意對信號進行延拓,應保留原信號的頻率和帶寬特性,否則會影響到真實信號[12]。
本文介紹了基于復Morlet小波變換的結構非線性振動模型參數識別算法,針對小波變換過程中出現的邊端效應及較少采樣點情況下參數識別精度低等問題,提出了采用BP神經網絡對原始信號進行預測延拓,并針對兩種非線性振動模型進行了數值模擬與分析。
在小波分析中,主要討論的函數空間是 L2(R)。L2(R)是指R上平方可積函數構成的函數空間,即:

如果 f(t)∈L2(R),則稱 f(t)為能量有限信號。L2(R)也被稱為能量有限的信號空間。如果ψ(t)∈L2(R),其傅里葉變換 ψ^(ω)滿足容許條件(Admissible Condition):

即Cψ有界,則稱ψ為一個基小波或母小波(Mother Wavelet)。將母小波經過伸縮和平移后,就可以得到一個小波序列:

其中,a,b∈R,且 a>0,稱 a為伸縮因子,b為平移因子。
定義下式:

為關于基小波ψ的連續小波變換。其中,ψ*(·)表示復共軛運算。從式(4)可知變換后的函數是二維的,即小波變換把原來的一維時域信號映射到為二維“時間—尺度”域上,以便分析信號的時—頻特性。定義以下變換:

為關于基小波ψ的小波逆變換。小波逆變換是把二維信號重構回原來的一維信號。
目前常用的小波函數有:Haar小波,Morlet小波,Daubechies小波,Cgau小波,復Morlet小波等。但是復解析小波變換將Hilbert變換與小波分析緊密結合在一起,具有很好的自適應分析能力[13]。本文采用復Mor-let小波變換進行參數識別,它的表達式為:

其中,fc為小波的中心頻率,fb為帶寬參數。嚴格來講,因復Morlet小波的傅里葉變換不滿足容許性條件,因此它不是一個真正的小波。但當時,復 Mor-let近似滿足容許條件,可以作為小波使用[14]。復Mor-let小波在時、頻域均具有很好的局部性,信號的頻率f與尺度a之間的關系為:

其中,fs為信號的采樣頻率。
對于有如下表達形式的信號:

其中,φ(t)=ω0t+β(t)。如果式(8)中信號的頻率變化率相對于振幅變化率來說大很多,則信號從整體來看是一個幅度變化平緩的信號即為漸進信號。對于弱Duffing系統的自由衰減響應而言,在阻尼較小的情況下可以看成漸進信號。漸進信號x(t)的復Morlet小波變換系數可近似表示為[15]:

大多數情況下是通過小波變換系數的模|Wx(a,b)|的最大值提取小波脊,但是尺度決定了分析窗的頻帶,則小波變換存在頻移現象,此現象會對參數識別精度有影響,代煜等[5]提出了從|Wx(a,b)中提取小波脊線以消除這種現象的影響。本文采用此方法,則小波脊線上的小波變換系數表示為:

其中,a=a(b)是b的函數。則根據式(10)可以得出:

其中,A(b)、φ(b)分別為信號的瞬時振幅和相位,則信號的瞬時圓頻率為:

則根據式(7)和式(11)可求出系統的瞬時頻率和瞬時振幅。
對于弱Duffing系統的立方剛度的自由振動方程為:

由多尺度法求得方程的一階近似解為:

其中,a0為初始振幅,β0為初相位。則系統的瞬時振幅和瞬時圓頻率分別為:


再由式(17)利用最小二乘法擬合出從而可得出參數ω0和α的識別值。
BP神經網絡是一種多層前饋神經網絡,該網絡的主要特點是信號前向傳遞,誤差反向傳播。在前向傳遞中,輸入信號從輸入層經隱含層逐層處理,直至輸出層。每一層的神經元只影響下一層神經元狀態。如果輸出層得不到期望輸出,則轉入反向傳播,根據預測誤差調整網絡權值和閾值,從而使BP神經網絡預測輸出不斷逼近期望輸出[16]。BP神經網絡的拓撲結構如圖1所示。

圖1 BP神經網絡Fig.1 BP neural network
在圖1中,x1,x2,…,xn是 BP神經網絡的輸入值,y1,…,ym是 BP神經網絡的預測值,ωij和 ωjk是神經網絡的權值。從圖1中可看出,BP神經網絡可看成一個非線性函數,網絡輸入值和預測值分別為函數的自變量和因變量。BP神經網絡預測前先要訓練網絡,通過訓練使網絡具有聯想和預測能力。訓練網絡的過程包括:網絡初始化、隱含層輸出計算、輸出層輸出計算、誤差計算、權值更新、閾值更新,最后判斷算法迭代是否結束,若沒結束,返回第二步[16]。
基于BP神經網絡對信號進行預測分為以下三步:BP神經網絡構建、BP神經網絡訓練和BP神經網絡預測。本文對于信號進行預測,是利用前兩時刻的值預測下一時刻的值,則選用BP神經網絡的結構為2-5-1,即輸入層有2個節點,隱含層有5個節點,輸出層有1個節點。本文先用已有的樣本點進行網絡訓練,再用訓練好的網絡進行預測。
本文首先采用弱Duffing非線性振動模型:·u·+ω20u+2μu2·u+αu3=0,ω0=100,μ=0.008,α=3,采樣頻率為1 000 Hz,振動時間為0~0.6 s,初始振幅為100 mm,用Simulink搭建了模型,然后采用四階龍格—庫塔法進行計算,振動響應如圖2實線所示。為了識別參數,采用復Morlet小波函數ψ(t)對響應信號和預測后的信號進行小波變換,根據上面提到的方法提取小波脊線,然后得到系統的瞬時振幅和瞬時圓頻率,進而識別未知參數。

圖2 原始信號和預測信號Fig.2 Original signal and predicted signal
4.1.1 邊端效應
對圖2的振動響應原始信號和預測信號分別進行小波變換,小波變換量圖分別見圖3、4,瞬時振幅與時間的關系及瞬時圓頻率與時間的關系分別見圖5、6,瞬時振幅與瞬時圓頻率之間的關系見圖7,小波系數平方的倒數與時間的關系見圖8。圖中實線表示原始信號經復Morlet小波變換得出的各種關系圖,虛線表示預測信號經復Morlet小波變換得出的各種關系圖。從圖5~8可以看出對利用BP神經網絡預測后的信號進行復Morlet小波變換能很好的解決邊端效應。由圖7可識別出ω0和α,由圖8可識別出μ,原始信號和預測后的信號的參數識別結果見表1。從表1可看出:直接用0~0.6 s的數據進行參數識別,預測后的信號參數識別值更接近真實值,誤差比直接用原始信號小得多;取0.2~0.5 s時間段的數據進行參數識別時,用預測后的信號誤差比用原始信號的誤差稍微減小。

圖3 原始信號小波變換量圖Fig.3Wavelet scalogram of the original signal

圖4 預測信號小波變換量圖Fig.4Wavelet scalogramof the predicted signal

圖5 瞬時振幅與時間的關系曲線Fig.5 The relation curve of instantaneous amplitude and time

圖6 瞬時圓頻率與時間的關系曲線Fig.6 The relation curve of instantaneous angular frequency and time

圖7 瞬時圓頻率與瞬時振幅的關系曲線Fig.7 The relation curve of instantaneous angular frequency and instantaneous amplitude

圖8小波系數平方的倒數與時間的關系曲線Fig.8 The relation curve of inverse square of the wavelet coefficients and time

表1 參數識別值及其誤差Tab.1 Parameter identification values and their error
另外,因實際采集的信號會被噪聲污染,所以也考慮了在仿真數據中疊加一定的白噪聲的情況下該方法的有效性。加入白噪聲的幅度由信噪比(SNR)來控制,信噪比定義為信號的均方差與噪聲的均方差之比。由于神經網絡預測對于噪聲比較敏感,所以本文先對加入白噪聲的信號進行了濾波處理,然后用BP神經網絡對處理后的信號進行預測延拓,最后用復morlet小波識別出未知參數。表2給出了不同噪聲水平下參數識別的相對誤差。定義相對誤差為:

從表2中可以看出,在信噪比達到5的情況下,參數識別值的相對誤差仍在10%以內,這表明該方法對噪聲具有一定的魯棒性。

表2 不同噪聲水平下參數識別的相對誤差Tab.2 The relative error of param eter identification under different noise levels
4.1.2 較少采樣點
對上述系統取0~0.2 s的振動響應如圖9實線所示,然后進行復Morlet小波變換,小波變換量圖見圖10,各變量的關系如圖12~15實線所示。從圖中的結果可以看出因采樣點太少及邊端效應的影響,無法對參數進行識別。采用BP神經網絡對信號進行預測延拓后,再進行復Morlet小波變換,小波變換量圖見圖11,各變量的關系如圖12~15虛線所示,取圖14、15的0.15~0.3 s時間段的數據進行參數識別,可以很好識別出未知參數,識別結果見表3。由此可知,在采樣點較少的情況下,利用BP神經網絡對信號進行延拓,是很有效的方法。

圖9 采樣點較少情況下的原始信號和預測信號Fig.9 Original signal and predicted signal in the less sampling points

圖10 原始信號的小波量圖Fig.10 Wavelet scalogram of the original signal

圖11 預測信號的小波量圖Fig.11 Wavelet scalogram of the predicted signal

圖12 瞬時振幅與時間的關系曲線Fig.12 The relation curve of instantaneous amplitude and time

圖13 瞬時圓頻率與時間的關系曲線Fig.13 The relation curve of instantaneous angular frequency and time

圖14 瞬時圓頻率與瞬時振幅的關系曲線Fig.14 The relation curve of instantaneous angular frequency and instantaneous amplitude

表3 參數識別值及其誤差Tab.3 Parameter identification values and their error
為了驗證本文的算法對于不同非線性振動模型的有效性,下面考慮第二個非線性振動模型為:

其利用多尺度法解得一階近似解為:

其中,a0為初始振幅,β0為初相位。同理系統的瞬時振幅和瞬時圓頻率分別為:


則有 ln(A(t))=ln(a0)-δt即(a0)-δt,從而可知在小波脊線上|的對數與時間t呈直線的關系,則直線斜率的絕對值為參數δ
的識別值。接著,由式(7)可得:ω(t)再由式(24)利用最小二乘法擬合出從而可得出參數ω0和ε的識別值。
對于上述模型取 δ=3.4,ε=3,ω0=120,同樣采樣頻率為1 000 Hz,振動時間為0~1.0 s,初始振幅為100 mm,用Simulink搭建模型,然后采用四階龍格—庫塔法進行計算,振動響應如圖16實線所示。為了識別參數,同樣采用復Morlet小波函數2對響應信號和預測后的信號進行小波變換,根據上面提到的方法提取小波脊線,然后得到系統的瞬時振幅和瞬時圓頻率,進而識別未知參數。

圖15 小波系數平方的倒數與時間的關系曲線Fig.15 The relation curve of inverse square of the wavelet coefficients and time

圖16 原始信號和預測信號Fig.16 Original signal and predicted signal

圖17 原始信號小波變換量圖Fig.17Wavelet scalogram of the original signal
4.2.1 邊端效應
對圖16的振動響應原始信號和預測信號分別進行小波變換,小波變換量圖分別見圖17、18,瞬時振幅與時間的關系及瞬時圓頻率與時間的關系分別見圖19、20,瞬時振幅與瞬時圓頻率之間的關系見圖21,小波系數的對數與時間的關系見圖22。圖中實線表示原始信號經復Morlet小波變換得出的各種關系圖,虛線表示預測信號經復Morlet小波變換得出的各種關系圖。從圖19~22可以看出對利用BP神經網絡預測后的信號進行復Morlet小波變換能很好的解決邊端效應。由圖21可識別出ω0和ε,由圖22可識別出δ,原始信號和預測后的信號的參數識別結果見表4。從表4可看出:直接用0~1.0 s的數據進行參數識別,預測后的信號參數識別值更接近真實值,誤差比直接用原始信號小;取0.2~0.8 s時間段的數據進行參數識別時,從圖19~22可以看出,在時間段0.2~0.8 s上預測前后數據是完全重合的,故在這個時間段上參數識別的值是一樣的。

圖18 預測信號小波變換量圖Fig.18Wavelet scalogram of the predicted signal

圖19 瞬時振幅與時間的關系曲線Fig.19 The relation curve of instantaneous amplitude and time

圖20 瞬時圓頻率與時間的關系曲線Fig.20 The relation curve of instantaneous angular frequency and time

圖21 瞬時圓頻率與瞬時振幅的關系曲線Fig.21 The relation curve of instantaneous angular frequency and instantaneous amplitude

圖22 小波系數的對數與時間的關系曲線Fig.22 The relation curve of logarithm of the wavelet coefficients and time

圖23 采樣點較少情況下的原始信號和預測信號Fig.23 Original signal and predicted signal in the less sampling points

表4 參數識別值及其誤差Tab.4 Parameter identification values and their error
同樣,此算例也考慮了噪聲影響的情況,表5給出了不同噪聲水平下參數識別的相對誤差。從表中可以看出信噪比達到4時,參數識別值的相對誤差在10%以內,這表明該方法對噪聲依然具有較好的魯棒性。

表5 不同噪聲水平下參數識別的相對誤差Tab.5 The relative error of parameter identification under different noise levels
4.2.2 較少采樣點
對算例二取0~0.3 s的振動響應如圖23實線所示,然后進行復Morlet小波變換,小波變換量圖見圖24,各變量的關系如圖26~29實線所示。從圖中的結果可以看出因采樣點太少及邊端效應的影響,無法對參數進行識別。采用BP神經網絡對信號進行預測延拓后,再進行復Morlet小波變換,小波變換量圖見圖25,各變量的關系如圖26~29虛線所示,取圖28~29的0.2~1.1 s時間段的數據進行參數識別,可以很好識別出未知參數,識別結果見表6。由此可知,對于此模型在采樣點較少的情況下,利用BP神經網絡對信號進行延拓,也是很有效的方法。

圖24 原始信號的小波量圖Fig.24Wavelet scalogram of the original signal

圖25 預測信號的小波量圖Fig.25Wavelet scalogram of the predicted signal

圖26 瞬時振幅與時間的關系曲線Fig.26 The relation curve of instantaneous amplitude and time

圖27 瞬時圓頻率與時間的關系曲線Fig.27 The relation curve of instantaneous angular frequency and time

圖28 瞬時圓頻率與瞬時振幅的關系曲線Fig.28 The relation curve of instantaneous angular frequency and instantaneous amplitude

圖29 小波系數的對數與時間的關系曲線Fig.29 The relation curve of logarithm of the wavelet coefficients and time

表6 參數識別值及其誤差Tab.6 Parameter identification values and their error
本文利用BP神經網絡很好地解決了信號復Morlet小波變換所產生的邊端效應,使非線性系統的參數識別精度有所提高。尤其是在采樣點較少的情況下,原始信號進行復Morlet小波變換后不能夠進行參數識別;而采用BP神經網絡進行信號的預測延拓,并且基于預測延拓的結果利用復Morlet小波變換進行參數識別,可較好地識別出結構非線性模型的未知參數。而且利用BP神經網絡預測的信號,能保持原始信號的變化趨勢,不失信號的原有特征。同時,本文也驗證了該方法對噪聲具有很好的魯棒性。訓練時所需的樣本點,原則上是越多效果越好,在本文算例中較少采樣點的情況下,用200左右的樣本點就能達到很好的效果,這為非線性系統的參數在線識別奠定了一定的基礎。
[1]Staszewski W J.Identification of non-linear systems using multi-scale ridges and skeletons of thewavelet transform[J].Journal of Sound and Vibration,1998,214(4):639-658.
[2]Ta M N,Lardiès J.Identification of weak nonlinearities on damping and stiffness by the continuous wavelet transform[J].Journal of Sound and Vibration,2006,293(1):16-37.
[3]任宜春.基于小波分析的結構參數識別方法研究[D].長沙:湖南大學,2007.
[4]伊廷華,李宏男,王國新.基于小波變換的結構模態參數識別[J].振動工程學報,2006,19(1):51-56.YI Ting-hua, LI Hong-nan, WANG Guo-xin. Structural modal parameter identification based on wavelet transform[J].Journal of Vibration Engineering,2006,19(1):51-56.
[5]代煜,孫和義,李慧鵬,等.基于小波變換的弱非線性阻尼和剛度辨識方法[J].振動與沖擊,2009,28(2):51-55.DAIYu,SUN He-yi,LI Hui-peng,et al.Identification of weak nonlin-earities in damping and stiffness based on wavelet transform [J].Journal of Vibration and Shock,2009,28(2):51-55.
[6]王超,任偉新,黃天立.基于小波的非線性結構系統識別[J].振動與沖擊,2009,28(3):10-13.WANG Chao, REN Wei-xin, HUANG Tian-li. System identification of a nonlinear structure based on wavelet transformation[J].Journal of Vibration and Shock,2009,28(3):10-13.
[7]史治宇,沈林.基于小波方法的時變動力系統參數識別[J].振動.測試與診斷,2008,28(2):108-112.SHIZhi-yu,SHEN Lin.Parameter identification of linear time-Varying dynamical system based on wavelet method[J].Journal of Vibration,Measurement&Diagnosis,2008,28(2):108-112.
[8]許鑫,史治宇,龍雙麗.基于小波狀態空間法的時變結構瞬時頻率識別[J].中國機械工程,2011,22(8):901-904.XUXin, SHI Zhi-yu, LONG Shuang-li. Instantaneous frequency identification of a time-varying structure using wavelet-based state space method[J].China Mechanical Engineering,2011,22(8):901-904.
[9]許鑫,史治宇.用于時變系統參數識別的狀態空間小波方法[J].工程力學,2011,28(3):23-28.XUXin,SHIZhi-yu.Parameter identification for time-varying system using state space and wavelet method [J].Engineering Mechanics,2011,28(3):23-28.
[10]Torrence C,Compo GP.A practical guide towaveletanalysis[J].Bulletin of the American Meteorological society,1998,79(1):61-78.
[11] Kijewski T L,Kareem A.Wavelet transforms for system identification and associated processing concerns[A].15th ASCE Engineering Mechanics Conference[C].Columbia University,New York,2002.6:1-8.
[12] Meyers S D,Kelly B G,O’Brien J J.An introduction to wavelet analysis in oceanography and meteorology:With application to the dispersion of Yanai waves[J].Monthly Weather Review,1993,121(10):2858-2866.
[13]于德介,成瓊,程軍圣.基于復解析小波變換的瞬時頻率分析方法[J].振動與沖擊,2004,23(1):108-109.YU De-jie,CHENG Qiong,CHENG Jun-sheng.Transient frequency analysis based on complex analytical wavelet transform and its application to fault diagnosis in gear drive[J].Journal of Vibration and Shock,2004,23(1):108-109.
[14]Yan B F,Miyamoto A,Brühwiler E.Wavelet transform-based modal parameter identification considering uncertainty[J].Journal of Sound and Vibration,2006,291(1-2):285-301.
[15] 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-664.
[16]史峰,王小川,郁 磊,等.Matlab神經網絡30個案例分析[M].北京:北京航空航天大學出版社,2010.