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

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

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

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

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

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

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

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

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

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

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

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

（2）矩陣塊分解

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

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

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

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

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

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

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

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

定義Oi的兩個(gè)子矩陣T1和T2：

由式（10）可以看出：

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

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

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

定義矩陣B：

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

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

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

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

（2）對(duì)于通過(guò)式（11）計(jì)算得到的每個(gè)模態(tài)m都可以在通過(guò)式（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)

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

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

其中diag（·）表示用矩陣的主對(duì)角元素組成行向量，（·）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行，（·）*表示共軛。

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

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

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

圖2 自動(dòng)識(shí)別算法流程圖Fig.2 The flowchart of automatic identification algorithm

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

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

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

傳統(tǒng)的方法是從穩(wěn)態(tài)圖中手動(dòng)選擇極點(diǎn)，由于疏忽或者經(jīng)驗(yàn)不足，選擇不恰當(dāng)?shù)臉O點(diǎn)將會(huì)造成識(shí)別結(jié)果的不準(zhǔn)確。由于模糊C均值算法存在著一些問(wèn)題，如需要人為指定模態(tài)數(shù)、一些參數(shù)的最佳值難以設(shè)置等［3，7］，因此本文借助譜系聚類法（hierarchical clustering methods）對(duì)剩下的計(jì)算結(jié)果進(jìn)行自動(dòng)選擇，將所有的計(jì)算結(jié)果組成一數(shù)據(jù)集，將距離在一定范圍內(nèi)的數(shù)據(jù)進(jìn)行聚類，認(rèn)為同一類中的數(shù)據(jù)屬于同一模態(tài)，最后選擇元素個(gè)數(shù)大于一定值的類作為識(shí)別結(jié)果。該方法在進(jìn)行聚類之前需要建立距離矩陣，因此需要定義一個(gè)反映模態(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è)定距離閾值對(duì)模態(tài)進(jìn)行聚類。由式（22）對(duì)距離的定義可以看出，如果模態(tài)i，j之間的頻率、阻尼比、模態(tài)振型、模態(tài)能量比值都在容差之內(nèi)，那么dij將小于1，因此一般設(shè)定距離閾值為1。當(dāng)聚類完成之后，統(tǒng)計(jì)每個(gè)聚類中包含的模態(tài)個(gè)數(shù)，認(rèn)為聚類數(shù)目大于閾值Nm的聚類為有效聚類，提取離類中心最近的數(shù)據(jù)為識(shí)別結(jié)果。

3 數(shù)值仿真

以一懸臂梁模型作為仿真實(shí)驗(yàn)件，結(jié)構(gòu)如圖3所示，材料屬性：梁長(zhǎng)度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

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

利用基于協(xié)方差驅(qū)動(dòng)的隨機(jī)子空間算法對(duì)數(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)的拾取，接著利用譜系聚類法自動(dòng)選定模態(tài)（Nm=20），表1是本文方法識(shí)別的結(jié)果與懸臂梁前5階理論值之間的對(duì)比（頻率、阻尼比以及對(duì)應(yīng)振型之間的MAC值），可以看出本文提出方法所識(shí)別的結(jié)果與理論值誤差很小（固有頻率的誤差在0.72%以內(nèi)，阻尼比誤差在6%以內(nèi)，對(duì)應(yīng)振型之間的MAC值為1）。數(shù)值仿真實(shí)驗(yàn)說(shuō)明了本文提出的方法可以在不需要人為參與的情況下正確地識(shí)別出物理模態(tài)。

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

4 實(shí)例分析

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

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

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

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

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

表2 華福橋識(shí)別結(jié)果對(duì)比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)地自動(dòng)拾取，提高了識(shí)別效率，免除人為因素的干擾。

（3）通過(guò)數(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］姜金輝，陳國(guó)平，張 方，等.模糊聚類法在試驗(yàn)?zāi)B(tài)參數(shù)識(shí)別分析中的應(yīng)用［J］.南京航空航天大學(xué)學(xué)報(bào)，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在振動(dòng)信號(hào)處理中的應(yīng)用［M］.北京：中國(guó)水利水電出版社，2006.

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