湯寶平，章國穩(wěn)，陳 卓

（1.重慶大學(xué) 機(jī)械傳動國家重點(diǎn)實(shí)驗(yàn)室，重慶 400030；2.重慶交通科研設(shè)計(jì)研究院 橋梁結(jié)構(gòu)動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，重慶 400067）

隨機(jī)子空間方法（SSI）［1-2］是近年來發(fā)展起來的一種行之有效的環(huán)境激勵模態(tài)參數(shù)識別方法，該方法直接工作于時域數(shù)據(jù)，沒有頻率分辨率誤差的問題，不但能準(zhǔn)確識別系統(tǒng)的頻率，而且能很好的識別系統(tǒng)的模態(tài)振型和阻尼。在識別過程中結(jié)構(gòu)模型的定階是最關(guān)鍵的環(huán)節(jié)之一，常見的做法是先對系統(tǒng)階次進(jìn)行過估計(jì)，然后結(jié)合穩(wěn)定圖進(jìn)行結(jié)果選取。目前穩(wěn)態(tài)圖中模態(tài)的選擇多數(shù)是通過人工完成，這不僅增加使用者的工作量，不適用于在線分析情況，而且由于使用者在認(rèn)識上存在著差異（易受虛假模態(tài)干擾），使得模態(tài)參數(shù)識別結(jié)果帶有一定的主觀性。因此為隨機(jī)子空間算法引入一種模態(tài)自動選取方法是一項(xiàng)亟待解決的工作。

文獻(xiàn)［3］提出借助模糊C均值算法對p_LSCF結(jié)果進(jìn)行自動拾取，文獻(xiàn)［4］提出一種改進(jìn)的譜系聚類算法對LSCF結(jié)果進(jìn)行自動選擇，都取得了良好的效果。本文以基于協(xié)方差驅(qū)動的隨機(jī)子空間算法進(jìn)行參數(shù)估計(jì)，針對其計(jì)算結(jié)果中大量虛假模態(tài)影響結(jié)果拾取的問題，提出一種能夠衡量模態(tài)可靠性的指標(biāo)稱之為模態(tài)相似指數(shù)，將其結(jié)合模態(tài)能量以剔除計(jì)算結(jié)果中由噪聲、模態(tài)過估計(jì)等因素引起的虛假模態(tài)；以頻率、阻尼比、模態(tài)振型、模態(tài)能量為聚類因子計(jì)算結(jié)果中各模態(tài)之間的相似性，采用譜系聚類法根據(jù)模態(tài)之間的相似性將計(jì)算結(jié)果分成若干類，提取元素個數(shù)大于一定值的類作為拾取結(jié)果。通過一個數(shù)值仿真以及實(shí)例分析驗(yàn)證本文方法可以實(shí)現(xiàn)系統(tǒng)物理模態(tài)地自動拾取。

1 隨機(jī)子空間虛假模態(tài)剔除

1.1 基于協(xié)方差的隨機(jī)子空間系統(tǒng)識別算法［1］

n自由度系統(tǒng)，其離散時間狀態(tài)方程為：

其中：A∈R2n×2n：狀態(tài)矩陣，B∈R2n×m：輸入矩陣，C∈Rl×2n：輸出矩陣，D∈Rl×m：直饋矩陣，m：輸入個數(shù)，l：輸出個數(shù)，Δt：采樣間隔。

定義輸出協(xié)方差矩陣Ri：

定義狀態(tài)-輸出協(xié)方差矩陣G：

可證明如下關(guān)系成立：

（1）構(gòu)造Toeplitz矩陣

（2）矩陣塊分解

將式（4）代入式（5）得到：

其中：Oi∈Ril×N為擴(kuò)展可觀測矩陣，Ti∈RN×li為擴(kuò)展可控矩陣，N為系統(tǒng)階次。

設(shè)W1和W2為兩個可逆加權(quán)矩陣［2，5］，對加權(quán) Toeplitz矩陣進(jìn)行SVD分解，可以得到：

其中：U1∈Rli×N，S1∈RN×N，V1∈Rli×N。

結(jié)合式（6）和式（7）可得：

其中：（·）+表示矩陣的偽逆。

（3）模態(tài)參數(shù)的識別

由式（6）可以看出，C為矩陣Oi的前l(fā)行，G為Ti的后l列，用Matlab語言表達(dá)如下：

定義Oi的兩個子矩陣T1和T2：

由式（10）可以看出：

對矩陣A進(jìn)行特征值分解：

其中：Λ∈CN×N是一個對角矩陣，由離散時間系統(tǒng)的極點(diǎn)λi組成。系統(tǒng)參數(shù)可以由狀態(tài)矩陣A、輸出矩陣C得到［1］。

1.2 利用模態(tài)相似指數(shù)剔除虛假模態(tài)

定義矩陣B：

可以看出，對于系統(tǒng)特征值a（系統(tǒng)的真實(shí)極點(diǎn)）下式成立：

可以得到狀態(tài)矩陣A新的估計(jì)方式：

可以看出通過式（11）和式（15）得到的控制矩陣A具有相同的系統(tǒng)特征值，即兩種方法所得結(jié)果中物理模態(tài)將會成對出現(xiàn)。而對于由噪聲、階次過估計(jì)等因素引起的虛假特征值則不具備以上性質(zhì)，即兩種方法所得結(jié)果中虛假模態(tài)不會成對出現(xiàn)。

以一數(shù)值仿真驗(yàn)證上述結(jié)論。一兩自由度線性時不變系統(tǒng)在隨機(jī)激勵下振動，固有頻率分別是13.78，18.38，阻尼比分別是 0.002 8、0.003 8，分別用式（11）和式（15）對響應(yīng)信號進(jìn)行極點(diǎn)估計(jì)（階數(shù)假設(shè)為20）。兩種算法結(jié)果的極點(diǎn)分布情況如圖1所示，可以看出，兩種方法結(jié)果對真實(shí)極點(diǎn)（物理模態(tài)）的估計(jì)是一致的（成對地出現(xiàn)），說明兩種算法能有效地識別出了系統(tǒng)的物理模態(tài)；但是它們得到的虛假極點(diǎn)（虛假模態(tài)）則是不同的，不會成對出現(xiàn)。

圖1 極點(diǎn)分布圖（‘o’：式（11）結(jié)果，‘* ’：式（15）結(jié)果）Fig.1 The distribution of poles（‘o’：the results of Eqs.（11），‘* ’：the results of Eqs.（15））

定義模態(tài)相似系數(shù)：

其中：fm、ξm為通過式（11）得到的頻率和阻尼比，fn、ξn為通過式（15）得到的頻率和阻尼比，df、dξ分別表示頻率、阻尼比的容差，本文選取為 0.01、0.05，Wf、Wξ為頻率、阻尼比在計(jì)算相似系數(shù)中的權(quán)值，本文取為0.7，0.3。

于是可以根據(jù)計(jì)算結(jié)果中的模態(tài)相似系數(shù)剔除虛假模態(tài)，基本步驟如下：

（1）在計(jì)算過程中分別用式（11）和式（15）計(jì)算系統(tǒng)的狀態(tài)矩陣得到兩組結(jié)果（頻率f、阻尼比ξ）。

（2）對于通過式（11）計(jì)算得到的每個模態(tài)m都可以在通過式（15）計(jì)算的結(jié)果中找到頻率最相近的模態(tài)n，并用式（16）計(jì)算相應(yīng)相似指數(shù)rm。如果模態(tài)m，n之間的頻率、阻尼比都在容差之內(nèi)，那么rm將小于1，因此一般設(shè)定閾值為1。如果rm小于閾值，則可認(rèn)為模態(tài)m與n屬于同一模態(tài)，即模態(tài)m可以在式（15）計(jì)算的結(jié)果中找到相同的模態(tài)，認(rèn)為其為真實(shí)模態(tài)，反之，則將其作為虛假模態(tài)予以剔除。

1.3 利用模態(tài)能量剔除虛假模態(tài)

由自然激勵技術(shù)法（NExT）［6］可知，白噪聲環(huán)境激勵下結(jié)構(gòu)響應(yīng)的自相關(guān)函數(shù)和脈沖響應(yīng)函數(shù)具有相似表達(dá)式：

它由類似于系統(tǒng)各階模態(tài)的自由響應(yīng)信號疊加而成。由式（2）對輸出協(xié)方差矩陣Ri的定義可得各輸出信號能量之和如下式所示：

其中diag（·）表示用矩陣的主對角元素組成行向量，（·）T表示轉(zhuǎn)置。

將式（4）代入式（18）可得：

將式（12）代入式（19）可得：

結(jié)合式（12）和式（20）可得各階模態(tài)所貢獻(xiàn)的能量為：

其中：M=Cψi（ψ-1）iG，ψi為矩陣 ψ 的第i列，（ψ-1）i為矩陣 ψ-1的第i行，（·）*表示共軛。

由輸出信號的組成（如式（17）所示）可知，結(jié)構(gòu)的輸出能量應(yīng)該是由各階物理模態(tài)的能量之和，由噪聲和模態(tài)過估計(jì)等因素帶來的虛假模態(tài)對其貢獻(xiàn)應(yīng)該為零。因此，可以利用各階模態(tài)能量貢獻(xiàn)來判斷其是否屬于系統(tǒng)的物理模態(tài)。

2 模態(tài)參數(shù)自動拾取

基本步驟如圖2所示，先利用隨機(jī)子空間算法估計(jì)系統(tǒng)的模態(tài)參數(shù)并且根據(jù)本文提出的方法得到模態(tài)相似指數(shù)以及相應(yīng)的模態(tài)能量，利用各模態(tài)的相似指數(shù)以及模態(tài)能量剔除部分虛假模態(tài)，對剩下的結(jié)果采用譜系聚類法進(jìn)行聚類，根據(jù)類元素個數(shù)完成模態(tài)的自動拾取。

圖2 自動識別算法流程圖Fig.2 The flowchart of automatic identification algorithm

2.1 剔除部分虛假模態(tài)

為了減少計(jì)算結(jié)果中虛假模態(tài)對后續(xù)聚類分析的影響，在識別模態(tài)參數(shù)的過程中采用式（16）計(jì)算各模態(tài)相似指數(shù)，設(shè)定相似指數(shù)閾值剔除部分虛假模態(tài)，又由第1.3節(jié)分析可知，可以通過模態(tài)能量來判斷虛假模態(tài)，對各假設(shè)模型階數(shù)下計(jì)算出來的能量進(jìn)行排序，將能量最大的前Np個模態(tài)保留用以進(jìn)行后續(xù)分析，以進(jìn)一步剔除計(jì)算結(jié)果中的虛假模態(tài)。

2.2 計(jì)算結(jié)果的聚類及物理模態(tài)的選取

傳統(tǒng)的方法是從穩(wěn)態(tài)圖中手動選擇極點(diǎn)，由于疏忽或者經(jīng)驗(yàn)不足，選擇不恰當(dāng)?shù)臉O點(diǎn)將會造成識別結(jié)果的不準(zhǔn)確。由于模糊C均值算法存在著一些問題，如需要人為指定模態(tài)數(shù)、一些參數(shù)的最佳值難以設(shè)置等［3，7］，因此本文借助譜系聚類法（hierarchical clustering methods）對剩下的計(jì)算結(jié)果進(jìn)行自動選擇，將所有的計(jì)算結(jié)果組成一數(shù)據(jù)集，將距離在一定范圍內(nèi)的數(shù)據(jù)進(jìn)行聚類，認(rèn)為同一類中的數(shù)據(jù)屬于同一模態(tài)，最后選擇元素個數(shù)大于一定值的類作為識別結(jié)果。該方法在進(jìn)行聚類之前需要建立距離矩陣，因此需要定義一個反映模態(tài)之間距離的統(tǒng)計(jì)量，以頻率f、阻尼比ξ、模態(tài)振型ψ和模態(tài)能量P作為為聚類因子，定義模態(tài)i，j之間的距離dij如下：

其中：df、dξ、dψ、dp分別表示頻率、阻尼比、模態(tài)振型、模態(tài)能量的容差，本文選取為 0.01、0.05、0.02、0.1；Wf、Wξ、Wψ、Wp表示頻率、阻尼比、模態(tài)振型、模態(tài)能量在計(jì)算模態(tài)距離中的權(quán)重，它們之和為1，本文選取為0.25、0.25、0.25、0.25。

在距離矩陣建立完畢之后，便可以設(shè)定距離閾值對模態(tài)進(jìn)行聚類。由式（22）對距離的定義可以看出，如果模態(tài)i，j之間的頻率、阻尼比、模態(tài)振型、模態(tài)能量比值都在容差之內(nèi)，那么dij將小于1，因此一般設(shè)定距離閾值為1。當(dāng)聚類完成之后，統(tǒng)計(jì)每個聚類中包含的模態(tài)個數(shù)，認(rèn)為聚類數(shù)目大于閾值Nm的聚類為有效聚類，提取離類中心最近的數(shù)據(jù)為識別結(jié)果。

3 數(shù)值仿真

以一懸臂梁模型作為仿真實(shí)驗(yàn)件，結(jié)構(gòu)如圖3所示，材料屬性：梁長度l=1m，橫截面積A=1×10-4m2，體密度D=7830 kg/m3，彈性模量E=2.068 ×1011，泊松比K=0.3。

圖3 懸臂梁模型Fig.3 The model of a cantilever beam

對梁施加白噪聲激勵，同時在各測點(diǎn)處（x1、x2、x3、x4、x5）測試位移響應(yīng)。采樣頻率為1 000 Hz，采樣時間為10 s，給每個通道加10%的隨機(jī)噪聲。

利用基于協(xié)方差驅(qū)動的隨機(jī)子空間算法對數(shù)據(jù)進(jìn)行分析，圖4是用傳統(tǒng)方法得到的穩(wěn)定圖（其中‘s’表示穩(wěn)定點(diǎn)，‘v’表示頻率和振型穩(wěn)定的點(diǎn)，‘d’表示頻率和阻尼比穩(wěn)定的點(diǎn)，‘f’表示頻率穩(wěn)定的點(diǎn)。），圖5為利用模態(tài)相似指數(shù)及模態(tài)能量（Np=5）剔除虛假模態(tài)后的穩(wěn)定圖，可以看出剔除虛假模態(tài)后所得到的穩(wěn)定圖相比原始穩(wěn)定圖少了許多虛假模態(tài)干擾，有利于物理模態(tài)的拾取，接著利用譜系聚類法自動選定模態(tài)（Nm=20），表1是本文方法識別的結(jié)果與懸臂梁前5階理論值之間的對比（頻率、阻尼比以及對應(yīng)振型之間的MAC值），可以看出本文提出方法所識別的結(jié)果與理論值誤差很?。ü逃蓄l率的誤差在0.72%以內(nèi)，阻尼比誤差在6%以內(nèi)，對應(yīng)振型之間的MAC值為1）。數(shù)值仿真實(shí)驗(yàn)說明了本文提出的方法可以在不需要人為參與的情況下正確地識別出物理模態(tài)。

表1 識別結(jié)果與理論值的比較Tab.1 The compration between identified result and the theoretical value

4 實(shí)例分析

為了進(jìn)一步驗(yàn)證本文方法的有效性，將該方法應(yīng)用于重慶華福橋模態(tài)參數(shù)識別實(shí)驗(yàn)，橋梁在過往行人、車輛以及自然環(huán)境中的風(fēng)流動等隨機(jī)激勵下振動，在橋梁豎向上共設(shè)置了11個測點(diǎn)，由加速度傳感器采集，每個測點(diǎn)的數(shù)據(jù)記錄長度為4 k，采樣頻率為200 Hz，圖6為橋梁立面圖。

圖6 華福橋橋梁立面Fig.6 The bridge elevation of hua fu bridge

對采集到的信號采用隨機(jī)子空間進(jìn)行參數(shù)識別，首先根據(jù)模態(tài)相似指數(shù)與模態(tài)能量（Np=15）對計(jì)算結(jié)果進(jìn)行虛假模態(tài)剔除，接著用本文提出的自動識別算法進(jìn)行模態(tài)自動拾?。∟m=20），將結(jié)果與原始隨機(jī)子空間方法識別方法結(jié)合穩(wěn)定圖人工選取的結(jié)果和特征系統(tǒng)實(shí)現(xiàn)算法（ERA）所得到的結(jié)果對比如表2所示。

結(jié)合圖7（剔除虛假模態(tài)后的穩(wěn)定圖）及表2可以看出，穩(wěn)定圖中穩(wěn)定性較好的數(shù)據(jù)已全部被選取，由于本文數(shù)據(jù)已經(jīng)過多種方法分析，原始算法及ERA結(jié)果的選取有了一些參考，因此通過它們得到的結(jié)果是可靠的，從表2可知，本文算法的自動識別結(jié)果與其他兩種算法的結(jié)果一致（ERA算法雖然沒有得到第1階和第5階模態(tài)，但其他相應(yīng)的模態(tài)結(jié)果與本文結(jié)果也是非常接近的），并且其大部分結(jié)果都最靠近三者相應(yīng)結(jié)果的平均值，以上進(jìn)一步說明了本文所提出方法得到的結(jié)果是可靠的。

圖7 本文方法的穩(wěn)定圖Fig.7 The stabilization diagram obtained by improved method

表2 華福橋識別結(jié)果對比Tab.2 Compration of the identification results of hua fu bridge

5 結(jié)論

（1）提出一種模態(tài)可靠性的衡量指標(biāo)稱之為模態(tài)相似指數(shù)，并結(jié)合各階模態(tài)的能量有效地剔除計(jì)算結(jié)果中的虛假模態(tài)。

（2）引入譜系聚類法進(jìn)行模態(tài)結(jié)果拾取，根據(jù)計(jì)算結(jié)果之間的相似性將計(jì)算結(jié)果分成若干類，提取類元素大于閾值的類作為拾取結(jié)果，實(shí)現(xiàn)物理模態(tài)地自動拾取，提高了識別效率，免除人為因素的干擾。

（3）通過數(shù)值仿真和實(shí)例分析驗(yàn)證本文方法的有效性。

［1］ Peeters B，Roeck G de.Reference-based stochastic subspace identification for output-only modal analysis［J］.Mechanical Systems and Signal Processing，1999，（13）：855-878.

［2］ Peter V O.Bart D M.Subspace identification for linear systems：theory-implementation-applications［M］.Dordrecht，the Netherlands：Kluwer Academic Publishers，1996.

［3］姜金輝，陳國平，張 方，等.模糊聚類法在試驗(yàn)?zāi)B(tài)參數(shù)識別分析中的應(yīng)用［J］.南京航空航天大學(xué)學(xué)報，2009，（3）：344-347.

［4］ Verboven P，Cauberghe B，Parloo E，et al.User-assisting tools for a fast frequency-domain modal parameter estimation method［J］.Mechanical System and Signal Processing，2004，（18）：759-780.

［5］Hermans L，Van der Auweraer H.Modal testing and analysis of structures under operational conditions： industrial applications［J］.Mechanical Systems and Signal Processing.1999，13（2）：193-216.

［6］王 濟(jì)，胡 曉.Matlab在振動信號處理中的應(yīng)用［M］.北京：中國水利水電出版社，2006.

［7］高新波.模糊聚類分析及應(yīng)用［M］.西安：西安電子科技大學(xué)出版社，2004.