基于系統響應瞬時特性的非線性系統識別

秦 安 壯， 楊 志 勛， 張 明 杰， 武 文 華， 張 文 首*

( 1.大連理工大學 工程力學系, 遼寧 大連 116024；2.大連理工大學 土木工程學院, 遼寧 大連 116024 )

基于系統響應瞬時特性的非線性系統識別

秦 安 壯1， 楊 志 勛1， 張 明 杰2， 武 文 華1， 張 文 首*1

( 1.大連理工大學 工程力學系, 遼寧 大連 116024；2.大連理工大學 土木工程學院, 遼寧 大連 116024 )

在Krylov-Bogoliubov-Mitropolski(KBM)法的基礎上，提出一種基于系統響應瞬時特性的非線性系統識別方法．該方法通過建立系統響應瞬時特性與系統參數之間的函數關系，從而一次性識別出所有系統參數．采用歸一化Hilbert變換(normalized Hilbert transform，NHT)和廣義過零(generalized zero-crossing，GZC)法求解信號瞬時振幅和瞬時頻率，通過算例驗證了兩種方法的效果．以Duffing方程和Vanderpol方程兩類非線性振動系統為例，驗證了所提系統識別方法的精度．算例表明，即使在系統響應受到較大噪聲污染時，該方法也有很好的識別精度．

瞬時頻率；歸一化Hilbert變換；廣義過零法；非線性系統識別

0 引 言

工程界中普遍存在的振動系統多為非線性系統，非線性系統的參數識別一直是工程界中的熱點和難點．近年來，許多學者提出了多種不同的方法，Kerschen等在文獻[1]中綜述了現有的方法．在現有的方法中，一種常用于弱非線性系統的方法是通過建立系統參數與系統響應瞬時特性(包括瞬時振幅和瞬時頻率)之間的關系函數，從而識別出系統參數．因此，準確地計算系統響應瞬時特性對準確識別系統參數具有直接的影響．求解信號瞬時特性常用的方法有小波變換和Hilbert-Huang 變換兩種，這兩種方法都在非線性系統識別中得到了廣泛的應用．Staszewski[2]提出了基于小波脊和小波骨架的非線性系統識別方法，并通過數值算例驗證了該方法的精度．Kijewski-Correa[3]介紹了基于morlet小波變換的系統識別方法在土木工程領域的應用．Feldman[4]通過Hilbert變換獲得系統響應瞬時特性的估計公式，進而識別了非線性系統的類型．Huang等[5]提出了一種改進的Hilbert變換方法，即歸一化Hilbert 變換(normalized Hilbert transform，NHT)，從而更準確地計算信號的瞬時頻率．Pai[6]將Hilbert-Huang變換應用到多種類型的非線性振動識別中，均取得了較好的識別精度．Feldman[4]綜述了Hilbert變換在機械振動系統識別領域中的應用．

本文在Krylov-Bogoliubov-Mitropolski(KBM)法[7]的基礎上，基于系統響應提出系統參數識別方法．此方法可一次性識別系統參數．將此方法應用到Duffing方程和Vanderpol方程兩類非線性振動系統中，以驗證在噪聲干擾情況下該方法的精度．

1 基本理論

1.1 歸一化Hilbert變換

非線性系統響應時程具有非平穩、非線性的特性．瞬時頻率是分析非平穩、非線性信號強有力的工具．早在1946年，Gabor就提出通過Hilbert變換構造解析信號從而求瞬時頻率的方法．對于任意單頻率分量信號

y(t)=Q(t)cos[φ(t)]

(1)

式中：y(t)為一多分量信號；Q(t)為瞬時振幅；φ(t)=2πf0t+θ(t)為瞬時相位角．

y(t)的解析信號定義為

Y(t)=y(t)+iq[y(t)]

(2)

H[y(t)]=H{Q(t)cos[φ(t)]}=

Q(t)H{cos[φ(t)]}

(3)

至此，y(t)的瞬時振幅和瞬時頻率分別表示為

(4)

(5)

然而，根據Bedrosian定理[8]，只有在Q(t)和cos[φ(t)]的頻譜完全分離時，式(3)才能嚴格成立．Nuttall等[9]還指出，調幅或調頻信號的Hilbert 變換并不等同于其正交分量，并給出了兩者之間的整體誤差：

(6)

其中F(f)是Q(t)eiθ(t)的傅里葉譜．

可見，式(4)和(5)不能準確求解信號的瞬時振幅和瞬時頻率．Huang提出了將單頻率分量信號歸一化再進行Hilbert變換的方法，如式(7)所示．

(7)

其中y1(t)是y(t)歸一化得到的信號；e1(t)是y(t) 的經驗包絡，由y(t)的極大值點三次樣條插值得到．歸一化的過程往往需要重復若干次，以確保歸一化后的信號幅值不超過1，即

(8)

定義y(t)的調頻部分

F(t)=yn(t)=cos[φ(t)]

(9)

顯然F(t)保留了所有與y(t)頻率相關的信息．由于F(t)的幅值等于1，對F(t)進行Hilbert變換不再受Bedrosian定理限制．此外，Huang對F(t)與其Hilbert變換之間的誤差給出了一種更為實用的表達式：

(10)

其中H[F(t)]為F(t)的Hilbert變換．

Huang把這種先歸一化再進行Hilbert變換的方法稱為歸一化Hilbert變換(NHT)，并在文獻[5]中驗證了這種方法的優越性．

1.2 廣義過零法

過零法是一種常用的計算信號局部頻率的方法，其基本思想是根據連續過零點之間的時間間隔來確定局部頻率，因為該方法計算的頻率在連續過零點之間為一個定值．Huang等[5]在過零法的基礎上提出將過零點和極值點視為控制點，從而將時間分辨率提高到四分之一波形，并把改進的方法命名為廣義過零(generalized zero-crossing, GZC)法．Huang 將兩個連續的同類型過零點或同類型極值點間的時間間隔定義為一個整周期，信號中任意一點同時屬于4個整周期，定義為T4j，j=1,2,3,4；將兩個連續的不同類型過零點或不同類型極值點間的時間間隔定義為一個半周期，信號中任意一點同時屬于兩個半周期，定義為T2j，j=1,2；將任意過零點與其相鄰極值點間的時間間隔定義為一個四分之一周期，信號中任意一點只屬于一個四分之一周期，定義為T1．根據三類周期的局部程度，分別賦予T1、T2j、T4j權重因子4、2、1．至此，任意點的角頻率可表示為

(11)

GZC通過測量控制點間的時間間隔來計算局部頻率，其物理意義非常明確[5]．此外，該方法不涉及任何形式的變換和微分，因此具有極好的魯棒性．由于實際采樣信號總是離散的，采樣信號不能總是準確地采集到極值點和過零點．因此本文對上述方法進行改進：對于過零點，通過對最接近其的兩個采樣點線性插值得到；對于極值點，通過對最接近其的7個采樣點二次拋物線擬合得到．

1.3 仿真信號分析

為驗證NHT和GZC兩種方法的可靠性，對以下調幅調頻信號進行分析：

x(t)=e-0.15t[sin(5πt)+0.1cos(5πt)]

(12)

信號時域波形示于圖1，采樣頻率和采樣時長分別為100 Hz和15 s．由Hilbert變換求得的瞬時振幅示于圖1，NHT和GZC求得的瞬時頻率示于圖2．為了使結果更為清晰，圖2只顯示了第4 s到第8 s的瞬時頻率結果．可見NHT方法的結果與實際瞬時頻率吻合良好，GZC方法的結果也能較好反映實際瞬時頻率的趨勢．

圖1 振動響應

圖2 原信號瞬時頻率

為考察噪聲對本文方法識別精度的影響，在x(t)中加入零均值白噪聲，并定義噪聲比例如下：

(13)

在x(t)加入5%的白噪聲后，求得的瞬時頻率示于圖3．可見NHT方法結果出現較大波動，GZC方法結果與加入噪聲前的結果基本相同，說明GZC方法抗噪性能較好．這是由于GZC方法不涉及任何形式的變換和求導，通過直接測量零點和極值點間的時間間隔獲得信號的局部頻率，因此抗噪性能良好[5]．

圖3 噪聲污染信號瞬時頻率

2 方法介紹

對于單自由度弱非線性動力系統，其自由振動方程可表示為

α．．+ω20α=εf(α，α．)

(14)

α．

和

α．．

分別是系統的位移、速度和加速度；

ω

0

是系統線性化方程的振動角頻率；

ε

是小參數；

f

是關于

α

和

α．

的非線性方程．

根據KBM法，式(14)的解可近似表示為[10]

α(t)=Q(t)cos[φ(t)]

(15)

式中：Q(t)為瞬時振幅；φ(t)=ω0t+θ(t)為瞬時相位角，φ(t)對時間的一階導數為瞬時頻率；Q(t)和θ(t)均為隨時間慢變的函數．

式(15)對時間求導得

α．

(

t

)=

Q．

(

t

)cos[

φ

(

t

)]-

Q

(

t

)

θ．

(

t

)× sin[

φ

(

t

)]-

Q

(

t

)

ω

0

sin[

φ

(

t

)]

(16)

由于Q(t)和θ(t)均為隨時間慢變的函數，式(16)右端前兩項將遠小于第3項，因此可以認為

α．

(

t

)=-

Q

(

t

)

ω

0

sin[

φ

(

t

)]

(17)

Q．

(

t

)cos[

φ

(

t

)]-

Q

(

t

)

θ．

(

t

)sin[

φ

(

t

)]=0

(18)

類似地，可以得到

α．．

-Q(t)ω20cos[φ(t)]-Q．(t)ω0sin[φ(t)]-

(

t

)=

Q

(

t

)

ω

0

θ．

(

t

)cos[

φ

(

t

)]

(19)

將式(17)和(19)代入式(14)，得

Q．

(

t

)

ω

0

sin[

φ

(

t

)]+

Q

(

t

)

ω

0

θ．

(

t

)cos[

φ

(

t

)]= -

εf

{

Q

(

t

)cos[

φ

(

t

)],-

Q

(

t

)

ω

0

sin[

φ

(

t

)]}

(20)

聯立式(18)和(20)，得

Q．(t)=-εω0f{Q(t)cos[φ(t)]，-Q(t)ω0sin[φ(t)]}sin[φ(t)]θ．(t)=-εQ(t)ω0f{Q(t)cos[φ(t)]，-Q(t)ω0sin[φ(t)]}cos[φ(t)]

(21)

Q．

(

t

)和

θ．

(

t

)在一個周期內進行平均，有

Q．(t)=-ε2πω0∫2π0f{Q(t)cos[φ(t)]，-Q(t)ω0sin[φ(t)]}sin[φ(t)]dφθ．(t)=-ε2πQ(t)ω0∫2π0f{Q(t)cos[φ(t)]，-Q(t)ω0sin[φ(t)]}cos[φ(t)]dφ

(22)

由于φ(t)=ω0t+θ(t)，從而

Q．(t)=-ε2πω0∫2π0f{Q(t)cos[φ(t)]，-Q(t)ω0sin[φ(t)]}sin[φ(t)]dφθ．(t)=ω0-ε2πQ(t)ω0∫2π0f{Q(t)cos[φ(t)]，-Q(t)ω0sin[φ(t)]}cos[φ(t)]dφ

(23)

對于一般的工程問題，可以通過適當假定非線性項的形式，代入式(23)從而建立系統響應瞬時特性與系統參數之間的函數關系．下面將通過算例進一步介紹該方法的流程．

3 數值算例

3.1 Duffing非線性振動系統

考慮下式所示的Duffing方程：

α．．+ω20α+k3α3+c1α．=0

(24)

即相當于式(14)中

(25)

將式(25)代入式(23)可得

(26)

至此，如果可以得到系統振動響應時程的瞬時振幅及瞬時頻率，即可根據式(26)最小二乘擬合識別系統參數ω0、k3、c1．

作為算例，取ω0=45 rad/s，k3=10，c1=0.1，采用四階龍格-庫塔法計算其在初始振幅為5 mm，初始速度為零時的振動響應，如圖4所示．

圖4 Duffing非線性振動系統的振動響應

圖5為該響應時程的功率譜密度，其中縱坐標是以10為底的對數坐標．可見響應包含7、21和35 Hz 3個頻率帶，高階成分能量很低．根據文獻[11-12]，EMD分解很難篩選出能量較低的頻率成分．因此，本例將振動響應時程通過一個零相位數字濾波器[13]從而獲得單頻率分量信號．

通過Hilbert變換得到的瞬時振幅示于圖4中，NHT和GZC方法求得的瞬時角頻率示于圖6中．在使用NHT方法求瞬時角頻率時，信號兩端分別用與其第一個和最后一個周期相同的波形延長，從而減輕端部效應．可見，兩種方法求得的瞬時角頻率趨勢十分相似，NHT方法的結果波動更大，這與Huang在文獻[5]中的結果相吻合．

圖5 Duffing非線性振動系統的功率譜密度

圖6 Duffing非線性振動系統的瞬時角頻率

將瞬時振幅和瞬時角頻率結果代入式(26)，通過最小二乘擬合識別系統參數，各參數識別結果示于表1．可見各參數識別結果具有很好的精度，即使在信號受到10%噪聲污染時，各參數識別結果誤差都在4%以內．因此，本文方法對該類系統識別精度較高，且具有良好的抗噪性．

表1 Duffing非線性振動系統的識別結果

3.2 Vanderpol非線性振動系統

如下式所示的Vanderpol方程：

α．．+ω20α+c1α．+c2α2α．=0

(27)

即相當于式(14)中

(28)

將式(28)代入式(23)可得

(29)

至此，如果可以得到系統振動響應時程的瞬時振幅及瞬時頻率，即可根據式(29)最小二乘擬合識別系統參數ω0、c1、c2．

作為算例，取ω0=45 rad/s，c1=0.1，c2=0.1，采用四階龍格-庫塔法計算其在初始振幅為5 mm，初始速度為零時的振動響應，結果如圖7所示．

圖7 Vanderpol非線性振動系統的振動響應

圖8為該響應時程的功率譜密度，其中縱坐標是以10為底的對數坐標．可見響應包含7、21和35 Hz 3個頻率帶，高階成分能量很低．同樣，本例將振動響應時程通過一個零相位數字濾波器從而獲得單頻率分量信號．

圖8 Vanderpol非線性振動系統的功率譜密度

通過Hilbert變換得到的瞬時振幅示于圖7中，NHT和GZC方法求得的瞬時角頻率示于圖9中．在使用NHT方法求瞬時角頻率時，信號兩端分別用與第一個和最后一個波形相同的正弦波延長，從而減輕了端部效應．可見，兩種方法求得的瞬時角頻率結果吻合良好．

將瞬時振幅和瞬時角頻率結果代入式(29)，通過最小二乘擬合識別系統參數，各參數識別結果示于表2．可見各參數識別結果具有很好的精度，即使在信號受到10%噪聲污染時，各參數識別結果誤差都在3%以內．因此，本文方法對該類系統識別精度較高，且具有良好的抗噪性．

圖9 Vanderpol非線性振動系統的瞬時角頻率

表2 Vanderpol非線性振動系統的識別結果

4 結 語

本文在KBM法的基礎上，提出了一種基于系統響應瞬時特性的非線性系統識別方法；采用歸一化Hilbert變換和廣義過零法求解信號瞬時振幅和瞬時頻率，通過算例驗證了兩種方法的效果；將本文提出的系統識別方法應用到Duffing方程和Vanderpol方程兩類非線性振動系統，發現即使在系統響應受到較大噪聲污染時，本文方法也有很好的識別精度．

[1] KERSCHEN G, GOLINVAL J C. Generation of accurate finite element models of nonlinear systems — application to an aeroplane-like structure [J]. Nonlinear Dynamics, 2005, 39(1/2):129-142.

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

[3] KIJEWSKI-CORREA T L. Full-scale measurements and system identification:A time-frequency perspective [D]. Notre Dame:University of Notre Dame, 2003.

[4] FELDMAN M. Non-linear system vibration analysis using Hilbert transform — I. Free vibration analysis method ′Freevib′ [J]. Mechanical Systems and Signal Processing, 1994, 8(2):119-127.

[5] HUANG N E, WU Zhaohua, LONG S R,etal. On instantaneous frequency [J]. Advances in Adaptive Data Analysis, 2009, 1(2):177-229.

[6] PAI P F. Nonlinear vibration characterization by signal decomposition [J]. Journal of Sound and Vibration, 2007, 307(3/4/5):527-544.

[7] BOGOLIUBOV N N, MITROPOLSKI Y A. Asymptotic Methods in the Theory of Non-Linear Oscillations [M]. New York:Gordon and Breach, 1961.

[8] BEDROSIAN E, RICE S O. Output properties of Volterra systems (nonlinear systems with memory) driven by harmonic and Gaussian inputs [J]. Proceedings of the IEEE, 1971, 59(12):1688-1707.

[9] NUTTALL A H, BEDROSIAN E. On the quadrature approximation to the Hilbert transform of modulated signals [J]. Proceedings of the IEEE, 1966, 54(10):1458-1459.

[10] SCHMIDT G, TONDL A. Non-Linear Vibrations [M]. Cambridge:Cambridge University Press, 2009.

[11] PENG Z K, TSE P W, CHU F L. An improved Hilbert-Huang transform and its application in vibration signal analysis [J]. Journal of Sound and Vibration, 2005, 286(1/2):187-205.

[12] PENG Z K, TSE P W, CHU F L. A comparison study of improved Hilbert-Huang transform and wavelet transform:Application to fault diagnosis for rolling bearing [J]. Mechanical Systems and Signal Processing, 2005, 19(5):974-988.

[13] 紀躍波,秦樹人,湯寶平. 零相位數字濾波器[J]. 重慶大學學報(自然科學版), 2000, 23(6):4-7.

JI Yuebo, QIN Shuren, TANG Baoping. Digital filtering with zero phase error [J]. Journal of Chongqing University (Natural Science Edition), 2000, 23(6):4-7. (in Chinese)

Nonlinear system identification based on instantaneous characteristics of dynamic response

QIN Anzhuang1, YANG Zhixun1, ZHANG Mingjie2, WU Wenhua1, ZHANG Wenshou*1

( 1.Department of Engineering Mechanics, Dalian University of Technology, Dalian 116024, China；2.School of Civil Engineering, Dalian University of Technology, Dalian 116024, China )

In the light of the Krylov-Bogoliubov-Mitropolski (KBM) method, a nonlinear system identification method is developed based on the instantaneous characteristics of the dynamic response of the system. The method identifies all the system parameters through establishing function relationship connecting system transient response characteristics with system parameters. The normalized Hilbert transform (NHT) and the generalized zero-crossing(GZC) method are introduced to calculate the instantaneous amplitude and frequency of the dynamic response. An example is applied to verify the efficiencies of these two methods. The proposed system identification method is applied to the nonlinear vibration systems of the Duffing and Vanderpol equations. Experimental results show that the method has good identification accuracy even when the dynamic responses are largely polluted by noises.

instantaneous frequency; normalized Hilbert transform; generalized zero-crossing method; nonlinear system identification

2016-06-02；

2017-03-16．

國家自然科學基金資助項目(11572072).

秦安壯(1991-)，男，碩士生，E-mail:zuoyetingfeng@126.com；張文首*(1963-)，男，博士，教授，E-mail:wszhang@dlut.edu.cn.

1000-8608(2017)03-0221-06

TN911.6；O32

A

10.7511/dllgxb201703001