基于EM算法和模態形式的狀態空間模型自降階工作模態分析

施袁鋒，朱正言，陳 鵬，戴靠山

(1. 四川大學建筑與環境學院深地科學與工程教育部重點實驗室，四川，成都 610065；2. 四川大學破壞力學與工程防災減災四川省重點實驗室，四川，成都 610065；3. 浙江大學建筑工程學院，浙江，杭州 310058)

模態分析是近年結構健康監測領域的熱點之一[1 ? 2]，其獲取的結構模態參數在既有結構的動力特性評價、結構損傷診斷、結構振動控制等領域具有重要應用。工作模態分析只需測量結構在環境激勵下的響應信號便可進行模態參數識別，具有操作簡單可行、試驗經濟等特點，因而在工程領域中得到了廣泛應用[3 ? 5]。近年來，很多學者分別在時域內和頻域內提出了各種工作模態分析方法[3 ? 7]。

對于時域內的模態分析方法如自回歸移動平均時序法(ARMA)、隨機子空間識別法(SSI)、特征系統實現算法(ERA)，模型階次的選擇或確定直接影響著模態參數識別的維度、準確性和穩定性[8]。環境激勵下的結構實測振動響應，因環境非平穩激勵、結構非線性行為、測量噪聲和建模誤差的影響，給模態分析中模型階次的確定帶來了實際困難。目前，基于狀態空間模型的模型階次的確定通常利用Hankel矩陣的奇異值分解并結合穩定圖的方法[9]。由于穩定圖法不能完全排除非平穩激勵、噪聲等干擾的影響，學者們提出了一些解決這一影響的改進方法。常軍等[10]提出了兩階段穩定圖的隨機子空間結構模態參數識別方法。易偉建和劉翔[11]提出了以奇異值分解后的殘差期望為依據的模型定階方法。王樹青等[12]從奇異值相對變化率、崔偉成等[13]從擬合誤差最小化原則、朱銳等[14]從奇異值百分比變化增量、廖聿宸等[15]從基于奇異值的函數擬合分別提出了相應的模型定階方法。由于建立Hankel矩陣本身的影響，上述基于奇異值分解和穩定圖的模型定階方法還是存在較大主觀性。Cara等[16]于是從一個新的角度出發，首先建立了一種模態形式的狀態空間模型，并提出以模態貢獻為指標的狀態空間模型階次選擇方法，以仿真和實際結構進行了驗證。由于該方法直接以對角化系統狀態矩陣出發，整個模態形式的狀態方程的模型參數和狀態變量為復數型，在一定程度上給后續計算帶來不便。

隨機子空間識別法在工作模態分析中已有廣泛應用[3]，但在受系統初值影響的振動響應以及環境激勵與測量噪聲相關的情況下存在一定的缺陷。針對此，有學者提出了考慮初值影響和考慮環境激勵和測量噪聲相關的建模，并采用期望最大(EM)迭代算法實現極大似然估計的方法[17 ? 19]。EM算法[20]原理簡單，不必對目標函數進行復雜微分計算，可直接通過迭代漸近得到極大似然值，因而在基于狀態空間模型的系統識別中也得到了實際應用。Cara等[17]、Matarazzo和Pakzad[18]使用隨機狀態空間建模并以EM算法進行模態參數的估計，但都忽略了實際環境激勵與測量噪聲的相關性影響。利用Gibson和Ninness[21]建立的系統噪聲和觀測噪聲間相關的狀態空間模型EM解法，張康和施袁鋒[19]考慮了完整的狀態空間參數化建模，獲得了更準確的識別結果。張靜等[22]應用EM算法實現對狀態空間模型的辨識，并利用AIC準則確定狀態空間階次，通過數值分析驗證了該方法的有效性。對于維度高的動力系統，EM算法也存在著如收斂速度較慢的問題，因此合理的模型降階將有助于EM的計算效率。

本文從上述基于狀態空間模型的工作模態分析中的不足出發，通過建立一種特殊的模態形式狀態空間模型，提出了一種結合EM算法和模態貢獻的自降階工作模態分析方法。該特殊的模態形式狀態空間模型的特點是所有的模型參數和狀態變量都為實數。通過將EM迭代過程中估計的模型參數轉換成模態形式后，引入表征各模態對于總響應的重要性的模態貢獻指標，并設定模態貢獻閾值來實現對狀態空間模型的自動降階。同時結合頻譜分析和阻尼比閾值剔除識別結果的虛假模態。最后通過數值模擬和實測數據分析以驗證方法的適用性和有效性。

1 動力系統建模

1.1 狀態空間模型

記結構在環境激勵下的動力響應觀測數據為Y，具體可表示成：

式中：yk∈Rp×1表示系統p個測量通道在k離散時刻的響應向量；Yl∈R1×N為第l個通道的數據總時長為N的響應向量。觀測響應數據yk可通過隨機狀態空間模型來模擬：

式中：xk∈Rn×1為系統狀態向量；n為狀態向量的維度，也即模型階次；A∈Rn×n為系統的狀態轉換矩陣；C∈Rp×n為系統的輸出矩陣，描述內部狀態如何轉換成輸出yk；wk∈Rn×1為系統噪聲；vk∈Rp×1為觀測噪聲。假設系統噪聲和觀測噪聲是具有相關性的高斯白噪聲，其聯合概率分布為：

式中：～表示概率分布服從； N(α,β)為高斯概率分布函數， α 為均值，β 為協方差；矩陣Q為系統噪聲wk自身的協方差；矩陣R為觀測噪聲vk自身的協方差；矩陣S為系統噪聲和觀測噪聲兩者的協方差。假定系統初始狀態x1也滿足高斯分布：

式 中：μ1∈Rn×1為 系 統 初 始 狀 態x1的 均 值；P1∈Rn×n為x1的協方差。狀態空間模型式(2)的參數化形式可表示成：

式中：θ表示模型所有自由參數組成的向量； ?為定義符號，表示由右側各項矩陣中自由參數組成。

式(2)為一般形式的隨機狀態空間模型，其固有模態頻率和阻尼比與矩陣A的復數特征值有關，即：

直接基于式(2)的工作模態分析為通過實測響應Y，采用合理參數估計方法確定模型階次n和θ的值，進而求解系統相應的模態參數。從實測數據建立動力系統模型的關鍵步驟之一是確定合理的模型階次，這對模型估計的準確性和可靠性，以 及計算效率都具有重要意義。

1.2 模態形式狀態空間模型

一般形式的狀態空間模型，通過相似變換，系統仍具有相同的模態信息和輸出性質。先介紹通過相似變換得到一種特殊的模態形式狀態空間模型，這種模態形式將用于模型降階的過程中。

引入一可逆模態形式變換矩陣T，將狀態向量xk變成：

代入式(2)，可得一特殊模態形式的狀態空間方程：

此時，模態形式的狀態空間模型(8)的參數化形式可寫成：

1.3 模態貢獻

假設模型參數已知，利用觀測數據Y，用卡爾曼濾波可以估計系統狀態和輸出。當使用模態形式的狀態空間模型式(8)時，其卡爾曼濾波為：

因此，利用卡爾曼濾波后，觀測數據Y可表示為：

式中：

模態貢獻指標的大小反映了相應系統特征值的重要性。如果ci越小，說明該模態越有可能是虛假模態，故由此可判斷系統模型階次可能過高。相關學者在研究模態貢獻時，將各階模態在所有輸出通道上的貢獻進行平均處理，以此定義模態的全局貢獻[16]。但在分析中發現，經過平均處理會削弱真實模態的貢獻，本文為了避免某些貢獻較小的真實模態因平均處理被剔除，定義模態貢獻為每個模態在各輸出通道中的最大貢獻值，用以代替平均值表征模態重要性，從而提高對真實模態的識別能力。

2 基于EM算法和模態形式的狀態空間模型自降階方法

由于極大似然估計方法具有漸近無偏估計、漸近一致性和漸近有效性等優異特性，本文采用極大似然法來實現隨機狀態空間模型參數 θ的估計，進而獲取結構的工作模態參數。然而模型參數的極大似然值的求解是一復雜的非線性優化問題，直接求解往往具有較大困難，特別是對于高維參數估計問題。EM迭代算法的提出，提高了一些問題的極大似然估計的實際應用性，其中包括狀態空間模型的參數估計問題[22]。前文已指出模型階次確定的重要性，結合模態形式狀態空間模型，在這里描述本文基于EM算法的模型降階識別的方法。

2.1 EM算法

式中：似然函數L(θ|Y)=lnp(Y|θ) ，其中p(Y|θ)是測量數據Y在模型 θ下的概率密度函數。引入X為未觀測值，令系統完整數據Z=(X,Y)。由貝葉斯公式可得：

式(20)兩邊同時取對數，有：

2.2 模型階次確定

由1.3節得到，模態響應直接反映了系統的固有特性和激勵對系統輸出的影響，而模態貢獻反映了各個模態的重要性。因此提出在EM迭代過程中在滿足數據擬合條件下，剔除模態貢獻小的模型參數進行模型的降階。顯然，更高階數的狀態空間模型包含更多的模態信息，可以更好地擬合響應數據，但也帶來了過于擬合的情況，由此獲得的模態信息中必然包括結構的真實模態和相當數量的可忽略的虛假模態，同時會增加計算時間。為篩選出真實可靠的模態信息，需要把其中模態貢獻指標較小的模態剔除，引入一個貢獻閾值 [c]用于刪除那部分貢獻很小的模態，實現模型的降階。

關于EM何時執行降階和收斂，以似然函數值相對變化率建立相關的判別式 ρ作為判據指標，并確定參數 ρ0作為降階的判斷閾值和 ρg作為全局收斂的閾值：

由于式(25)～式(33)為全參數模態形式EM迭代結果，因此在EM迭代過程中滿足降階條件時，需將模型參數轉換到模態形式后計算各模態貢獻，當模態貢獻中的最小值小于設定的閾值[c]并且擬合殘差滿足白噪聲要求時，便刪除與之對應的模態信息，實現降階。當最小貢獻對應的模態為復模態，需要去除參數矩陣中對應的2次模型階次；當最小貢獻的模態對應為實模態，則去 除參數矩陣中對應的1次模型階次。

2.3 模型估計和自降階過程

基于EM算法和模態形式的狀態空間模型自降階工作模態分析的過程見圖1，具體執行步驟為：

圖1 自降階工作模態分析流程圖Fig.1 Flowchart of operational modal analysis with auto model order reduction

1) 根據觀測響應數據Y，使用隨機子空間方法(SSI)確定一個相對較大的模型階次初值n0和模型參數的估計初值 θ0。

4) 當最小模態貢獻指標超過貢獻閾值 [c]時，此時得到算法建議的最終模型階數n，EM算法迭代繼續至滿足EM全局收斂條件 ρg時停止。

此外，執行降階EM算法需要設定幾個參數，這里介紹相關參數的設置原則：

1) 狀態空間模型階次初值n0：為了使更多的真實模態信息被包含，同時狀態空間階次又不過分大，檢驗用SSI方法確定模型階次和估計參數初值后得到響應估計殘差ε ，當殘差 ε的頻譜圖中無明顯孤立峰值時可認為識別系統在整體上符合白噪聲假設條件。此時可認為選取的模型階次初值n0大小是合適的，否則，增大模型階次初值。

2) 模態貢獻閾值 [c]：使用SSI方法試算得到模態信息，然后計算出模態響應和各模態的模態貢獻；結合頻譜圖初步確定結構模態，特別是具有獨立峰值但峰值低而窄的可能結構模態，在結構模態貢獻中找到最小的模態貢獻；確定比最小結構模態貢獻稍小的值作為模態貢獻閾值，建議以1‰為基數。

3) 降階條件 ρ0：通常EM算法收斂閾值 ρ0小于 10?5即可，所以每次降階前確保EM達到一定程度收斂，這可根據模型階次的大小進行相應調整，狀態空間階數較小時 ρ0可以更小一點，狀態空間階數較大時 ρ0可適當放寬。

4) 全局收斂標準 ρg： ρg為不能降階時EM算法的最終迭代收斂標準，其可在 ρ0的基礎上進行適當降低，一般降低一個數量級即可。

5) 虛假結構模態剔除：除了在模型降階過程中模態貢獻小的虛假模態之外，最終模型的模態信息中還會存在一些非結構模態的虛假模態。為了對這些非結構的虛假模態加以判別，還需結合輸出響應的頻譜分析，相互佐證以提高識別結果的可信度。同時將模態的阻尼比信息考慮進來，一般而言，環境激勵下的建筑結構的阻尼比通常小于5%，并且考慮到以往的研究提到阻尼比的識別沒有特別可靠的方法，所以針對結構阻尼比判定可以適當放寬[23]?；诖耍O定結構阻尼比閾值為 ζc=10%，當識別模型某個模態的阻尼比超過該閾值便視為虛假模態并予以剔除。也可以通過分析模態參數的不確定性水平[5]，把不確定性水平大的模態認定為虛假模態予以剔除。

3 應用舉例

下面以數值模擬和現場實測的結構振動響應分析來對本文提出的方法進行驗證說明。

3.1 6自由度結構數值模擬分析

一個6層剪切型結構簡圖如圖2所示，各質點質量為mi=1t ，i=1,2,···,6；層間剛度分別為k1=1300 ，k2=1200 ，k3=1100 ，k4=1000，k5=900 ，k6=800，單位kN/m。結構的阻尼形式假設為Rayleigh阻尼，其第1階和第6階模態對應阻尼比均為5%。所有自由度上施加相同強度的白噪聲作為模擬結構的環境激勵。采集所有自由度上的加速度響應，采樣頻率fs=50Hz，時長為100 s。此外在結構的輸出響應中加入相同強度的白噪聲以模擬傳感器采集數據時的觀測噪聲。

圖2 6自由度結構模型簡圖Fig.2 Schematic diagram of 6-DOF structural model

圖3給出了采集加速度響應的奇異值頻譜圖，通過其峰值基本可以看出結構的前5階模態位置，但第6階模態并不明顯?；诓杉?層加速度振動響應信號，下面利用本文方法進行工作模態分析。

圖3 6自由度加速度響應奇異值頻譜圖Fig.3 Singular value spectrum of measured accelerations of 6-DOF structure

此6自由度結構對應于狀態空間模型的理論階次為12。因為一般情況下模型階次事先未知，首先假定一個較高模型階次初值n0=24，此模型階次得到的估計殘差一定能滿足白噪聲條件。EM降 階 算 法 中 設 定 為 ρ0=1×10?6，ρg=1×10?7，[c]=0.05。圖4顯示了EM降階識別的過程，圖中給出了每一次降階對應的似然函數值的變化以及根據收斂標準得到的最終模型階數和總迭代次數。最終的模型階次確定為n=12，與真實的模型階次相同。

圖4 6自由度結構EM自降階迭代過程Fig.4 EM iteration process with auto order reduction of 6-DOF structure

最終識別的模態參數識別結果見表1。對于此數值模擬結構，通過本文識別方法確定了最小的模型階次后，識別出來的模態全部為結構真實模態，無虛假模態的存在。盡管圖3中后幾階頻率峰值不明顯，由表1的結果可以看出本文方法依然獲得了非常準確的模態參數識別結果。表1中同時給出了每個模態所對應的模態貢獻值，可以看出一階頻率對應的通道最大貢獻是所有模態中最小的，只有0.063，表明一階模態對結構的總響應的貢獻不大，為非優勢模態，但實際上該模態實際上是結構的真實模態，不應該被剔除，說明合理設置模態貢獻閾值是非常關鍵的一步。通常，我們需要在一定程度上借助頻譜圖的重要模態確定可行的模態貢獻閾值。圖5給出了模態振型的識別結果，可見所有振型都吻合得很好。綜合表1和圖5，帶有降階EM算法對結構響應進行了有效處理，并獲得了非常好的結構模態信息。

表1 6自由度模態分析結果Table 1 Modal analysis results of 6-DOF structure

圖5 6自由度振型識別結果Fig.5 Identified mode shapes of 6-DOF structure

為演示模態貢獻的概念，以下以質點2的振動響應數據為例進行說明。首先，圖6(a)展示了質點2處10 s～20 s時間段內的真實和采集的加速度響應的時程曲線。由于模擬測量噪聲較大，可見圖6(a)中真實加速度序列和測量加速度序列相差較大。使用模態形式狀態空間模型求得各階模態在質點2處的模態響應和模態貢獻，通過振型疊加得到質點2的擬合加速度響應。各階模態響應及貢獻如圖7所示，可以看出模態1和模態4的振幅都非常小，模態貢獻都小于1%，比表1中對應模態的通道最大貢獻還小很多；模態2、模態3、模態6的響應幅值較大，對應的模態貢獻也比較大，意味著這3個模態于質點2而言是優勢模態。

圖6 6自由度結構質點2處(10 s～20 s)加速度響應分析Fig.6 Analysis of acceleration response at mass No.2 of 6-DOF structure (10 s-20 s)

圖7 6自由度結構質點2處(10 s～20 s)各階模態響應及貢獻Fig.7 Modal response and contribution ratio at mass No.2 of 6-DOF structure (10 s-20 s)

將質點2處各階模態響應疊加可得到其最終擬合加速度曲線，如圖6(b)所示。從圖6(b)中看出，得到的擬合加速度曲線在幅值較高的地方與真實加速度有一定誤差，這是由于原始信號中包含了很高的噪聲成分，從而影響了對結構響應的估計，但與圖6(a)相比，擬合的加速度響應與真實響應較接近。

3.2 廣州塔實測振動分析

為了驗證本文方法對實際結構的模態分析能力，選擇廣州新電視塔(以下簡稱“廣州塔”)為應用對象。國內外已有許多學者[23 ? 26]對廣州塔在環境激勵下進行了模態分析研究。本節通過對廣州塔的環境振動數據進行模態分析，嘗試在0 Hz～5 Hz間確定一個合適的狀態空間階次并估計可靠的模態信息。

廣州塔傳感器布置如圖8所示，在8個不同高度總共安裝有20個水平加速度傳感器，分別沿截面的長軸和短軸方向布置。2010年1月19日—20日，基準測試工作組采集了廣州塔在環境激勵下24 h的加速度響應數據。這里隨機選用1月20日18:00～19:00間的測量數據進行分析(數據來源：http://www.zn903.com/ceyxia/benchmark/index.htm)。原始測量數據的采樣頻率為50 Hz，因此每個通道的采樣點數量N=180000?？紤]到主要考慮結構低階范圍內的模態信息，我們對原始信號進行濾波重采樣處理，以減少數據量并提高計算效率。對原始測量數據使用巴特沃斯低通濾波，截斷頻率為5 Hz，然后以10 Hz進行重采樣，采樣點降為N=36000。

圖8 廣州塔水平加速度傳感器布置簡圖Fig.8 Layout diagram of horizontal acceleration sensor placement for Canton Tower

圖9給出了部分傳感器(傳感器6、12和16)在18:00～19:00的加速度響應曲線，所有傳感器獲得的加速度響應的奇異值頻譜圖如圖10所示。可以發現3 Hz～5 Hz結構的模態信息成分非常復雜，不易辨認。

圖9 廣州塔傳感器(6、12和16)加速度響應曲線Fig.9 Acceleration responses from sensors(6, 12 and 16) on Canton Tower

針對重采樣的所有數據，應用本文方法展開模態識別。首先，設置初始狀態空間階數n0=160 ，EM算法降階條件設為ρ0=5×10?6，全局收斂率 ρg=1×10?6，設置模態貢獻閾值 [c]=0.005，帶有模型降階識別過程如圖11所示，圖中每次似然函數值突降點為進行了模型降階，最后確定的模型階數為n=119，相比于模型階次初值有較大減少。由于廣州塔結構復雜，存在較多的模態，所以需要較大模型階次得到合理的識別模型，這一點也可以從圖10中預見到。圖12展示了最終確定模型階次后的估計殘差奇異值頻譜圖。盡管圖12中也存在一些小峰值，但因為這些峰值對應的模態貢獻都非常小，可認為其為噪聲，殘差頻譜圖在整體上看可認為符合高斯白噪聲分布。綜合圖11和圖12，可以認為模型階次得到了合理的確定。

圖10 廣州塔加速度響應奇異值頻譜圖Fig.10 Singular value spectrum of acceleration responses of Canton Tower

圖11 廣州塔EM自降階迭代過程Fig.11 EM iteration process with auto order reduction of Canton Tower

圖12 廣州塔識別模型的響應估計殘差奇異值頻譜圖Fig.12 Singular value spectrum of residuals of the identified model for Canton Tower

最終得到的模型所有模態信息包含13個實模態和53個復模態，但由于激勵假設誤差、建模誤差、噪聲等影響，這些模態信息既包含結構的真實模態，又包含虛假模態。因此，需要結合結構模態的特征來對識別結果進行篩選處理。利用加速度響應奇異值頻譜圖及一般結構阻尼判據ζc=10%，對虛假模態進行剔除。將剔除虛假模態后的最終前12階模態信息與Ye等[23]使用NExTERA方法得到的結果列于表2，可以看出本文方法與NExT-ERA方法對于頻率的識別具有很好的一致性，兩者得到的模態信息極為接近。此外本算例中使用的數據信息頻寬為0 Hz～5 Hz，在這個范圍內可以獲得更多對結構響應貢獻較大的模態，這對廣州塔的模型更新具有更多應用價值。

表2 廣州塔本文方法與NExT-ERA[23]識別結果比較Table 2 Comparison of the identified results of the proposed method and NExT-ERA[23] for Canton Tower

圖13和圖14畫出了廣州塔前12階模態振型分別在x方向和y方向的分量，可以通過比較判斷出模態振型在2個方向上的耦合情況，比如模態1和模態2都在一個特定方向呈現出彎曲型特征，而在另外一個方向的振幅很小，說明模態1和模態2在相互正交的方向上耦合效應不明顯，符合結構的一階振型特征。

圖13 廣州塔前12階振型在x方向的分量(橫截面長軸)Fig.13 The x-component of the first 12 identified modes of Canton Tower (the longitudinal axis of the cross section)

圖14 廣州塔前12階振型在y方向的分量(橫截面短軸)Fig.14 The y-component of the first 12 identified modes of Canton Tower (the transverse axis of the cross section)

4 結論

本文提出了利用一種特殊的模態形式狀態空間模型，在EM迭代算法求解模型參數的極大似然估計中實現模型降階和模態參數估計的方法。重點可以概況如下：

(1) 基于模態形式狀態空間模型，可利用卡爾曼濾波獲得各階模態的模態響應量。

(2) 提出了表征各階模態對總響應的重要性的模態貢獻指標，以此作為模型降階標準較之振型參與系數等指標更具有物理意義和有效性。

(3) 提出的方法中一些閾值的選擇，仍帶有一定的經驗性，需結合數據和結構特征進行調整。

(4) 該模型降階方法的結果雖不容易遺漏真實模態，但結構虛假模態不能全部剔除。因此最后仍需要結合結構模態信息的特征對系統模型全部模態結果進一步篩選。

(5) 通過模擬結構數據和廣州塔實測數據的分析，表明本文的方法具有良好的準確率和穩健性。