基于系統(tǒng)響應(yīng)瞬時(shí)特性的非線性系統(tǒng)識別

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

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

基于系統(tǒng)響應(yīng)瞬時(shí)特性的非線性系統(tǒng)識別

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

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

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

瞬時(shí)頻率；歸一化Hilbert變換；廣義過零法；非線性系統(tǒng)識別

0 引 言

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

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

1 基本理論

1.1 歸一化Hilbert變換

非線性系統(tǒng)響應(yīng)時(shí)程具有非平穩(wěn)、非線性的特性．瞬時(shí)頻率是分析非平穩(wěn)、非線性信號強(qiáng)有力的工具．早在1946年，Gabor就提出通過Hilbert變換構(gòu)造解析信號從而求瞬時(shí)頻率的方法．對于任意單頻率分量信號

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

(1)

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

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)的瞬時(shí)振幅和瞬時(shí)頻率分別表示為

(4)

(5)

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

(6)

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

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

(7)

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

(8)

定義y(t)的調(diào)頻部分

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

(9)

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

(10)

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

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

1.2 廣義過零法

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

(11)

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

1.3 仿真信號分析

為驗(yàn)證NHT和GZC兩種方法的可靠性，對以下調(diào)幅調(diào)頻信號進(jìn)行分析：

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

(12)

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

圖1 振動響應(yīng)

圖2 原信號瞬時(shí)頻率

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

(13)

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

圖3 噪聲污染信號瞬時(shí)頻率

2 方法介紹

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

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

(14)

α．

和

α．．

分別是系統(tǒng)的位移、速度和加速度；

ω

0

是系統(tǒng)線性化方程的振動角頻率；

ε

是小參數(shù)；

f

是關(guān)于

α

和

α．

的非線性方程．

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

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

(15)

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

式(15)對時(shí)間求導(dǎo)得

α．

(

t

)=

Q．

(

t

)cos[

φ

(

t

)]-

Q

(

t

)

θ．

(

t

)× sin[

φ

(

t

)]-

Q

(

t

)

ω

0

sin[

φ

(

t

)]

(16)

由于Q(t)和θ(t)均為隨時(shí)間慢變的函數(shù)，式(16)右端前兩項(xiàng)將遠(yuǎn)小于第3項(xiàng)，因此可以認(rèn)為

α．

(

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)

聯(lián)立式(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

)在一個(gè)周期內(nèi)進(jìn)行平均，有

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)

對于一般的工程問題，可以通過適當(dāng)假定非線性項(xiàng)的形式，代入式(23)從而建立系統(tǒng)響應(yīng)瞬時(shí)特性與系統(tǒng)參數(shù)之間的函數(shù)關(guān)系．下面將通過算例進(jìn)一步介紹該方法的流程．

3 數(shù)值算例

3.1 Duffing非線性振動系統(tǒng)

考慮下式所示的Duffing方程：

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

(24)

即相當(dāng)于式(14)中

(25)

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

(26)

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

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

圖4 Duffing非線性振動系統(tǒng)的振動響應(yīng)

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

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

圖5 Duffing非線性振動系統(tǒng)的功率譜密度

圖6 Duffing非線性振動系統(tǒng)的瞬時(shí)角頻率

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

表1 Duffing非線性振動系統(tǒng)的識別結(jié)果

3.2 Vanderpol非線性振動系統(tǒng)

如下式所示的Vanderpol方程：

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

(27)

即相當(dāng)于式(14)中

(28)

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

(29)

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

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

圖7 Vanderpol非線性振動系統(tǒng)的振動響應(yīng)

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

圖8 Vanderpol非線性振動系統(tǒng)的功率譜密度

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

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

圖9 Vanderpol非線性振動系統(tǒng)的瞬時(shí)角頻率

表2 Vanderpol非線性振動系統(tǒng)的識別結(jié)果

4 結(jié) 語

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

[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ì)躍波,秦樹人,湯寶平. 零相位數(shù)字濾波器[J]. 重慶大學(xué)學(xué)報(bào)(自然科學(xué)版), 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．

國家自然科學(xué)基金資助項(xiàng)目(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