基于譜系聚類的隨機子空間模態參數自動識別

湯寶平，章國穩，陳 卓

（1.重慶大學 機械傳動國家重點實驗室，重慶 400030；2.重慶交通科研設計研究院 橋梁結構動力學國家重點實驗室，重慶 400067）

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

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

1 隨機子空間虛假模態剔除

1.1 基于協方差的隨機子空間系統識別算法［1］

n自由度系統，其離散時間狀態方程為：

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

定義輸出協方差矩陣Ri：

定義狀態-輸出協方差矩陣G：

可證明如下關系成立：

（1）構造Toeplitz矩陣

（2）矩陣塊分解

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

其中：Oi∈Ril×N為擴展可觀測矩陣，Ti∈RN×li為擴展可控矩陣，N為系統階次。

設W1和W2為兩個可逆加權矩陣［2，5］，對加權 Toeplitz矩陣進行SVD分解，可以得到：

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

結合式（6）和式（7）可得：

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

（3）模態參數的識別

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

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

由式（10）可以看出：

對矩陣A進行特征值分解：

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

1.2 利用模態相似指數剔除虛假模態

定義矩陣B：

可以看出，對于系統特征值a（系統的真實極點）下式成立：

可以得到狀態矩陣A新的估計方式：

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

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

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

定義模態相似系數：

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

于是可以根據計算結果中的模態相似系數剔除虛假模態，基本步驟如下：

（1）在計算過程中分別用式（11）和式（15）計算系統的狀態矩陣得到兩組結果（頻率f、阻尼比ξ）。

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

1.3 利用模態能量剔除虛假模態

由自然激勵技術法（NExT）［6］可知，白噪聲環境激勵下結構響應的自相關函數和脈沖響應函數具有相似表達式：

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

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

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

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

結合式（12）和式（20）可得各階模態所貢獻的能量為：

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

由輸出信號的組成（如式（17）所示）可知，結構的輸出能量應該是由各階物理模態的能量之和，由噪聲和模態過估計等因素帶來的虛假模態對其貢獻應該為零。因此，可以利用各階模態能量貢獻來判斷其是否屬于系統的物理模態。

2 模態參數自動拾取

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

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

2.1 剔除部分虛假模態

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

2.2 計算結果的聚類及物理模態的選取

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

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

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

3 數值仿真

以一懸臂梁模型作為仿真實驗件，結構如圖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

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

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

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

4 實例分析

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

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

對采集到的信號采用隨機子空間進行參數識別，首先根據模態相似指數與模態能量（Np=15）對計算結果進行虛假模態剔除，接著用本文提出的自動識別算法進行模態自動拾取（Nm=20），將結果與原始隨機子空間方法識別方法結合穩定圖人工選取的結果和特征系統實現算法（ERA）所得到的結果對比如表2所示。

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

圖7 本文方法的穩定圖Fig.7 The stabilization diagram obtained by improved method

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

5 結論

（1）提出一種模態可靠性的衡量指標稱之為模態相似指數，并結合各階模態的能量有效地剔除計算結果中的虛假模態。

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

（3）通過數值仿真和實例分析驗證本文方法的有效性。

［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］姜金輝，陳國平，張 方，等.模糊聚類法在試驗模態參數識別分析中的應用［J］.南京航空航天大學學報，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］王 濟，胡 曉.Matlab在振動信號處理中的應用［M］.北京：中國水利水電出版社，2006.

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