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

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

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

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

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

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

1 動力系統(tǒng)建模

1.1 狀態(tài)空間模型

記結(jié)構(gòu)在環(huán)境激勵下的動力響應(yīng)觀測數(shù)據(jù)為Y，具體可表示成：

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

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

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

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

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

式(2)為一般形式的隨機(jī)狀態(tài)空間模型，其固有模態(tài)頻率和阻尼比與矩陣A的復(fù)數(shù)特征值有關(guān)，即：

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

1.2 模態(tài)形式狀態(tài)空間模型

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

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

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

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

1.3 模態(tài)貢獻(xiàn)

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

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

式中：

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

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

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

2.1 EM算法

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

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

2.2 模型階次確定

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

關(guān)于EM何時執(zhí)行降階和收斂，以似然函數(shù)值相對變化率建立相關(guān)的判別式 ρ作為判據(jù)指標(biāo)，并確定參數(shù) ρ0作為降階的判斷閾值和 ρg作為全局收斂的閾值：

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

2.3 模型估計和自降階過程

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

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

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

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

此外，執(zhí)行降階EM算法需要設(shè)定幾個參數(shù)，這里介紹相關(guān)參數(shù)的設(shè)置原則：

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

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

3) 降階條件 ρ0：通常EM算法收斂閾值 ρ0小于 10?5即可，所以每次降階前確保EM達(dá)到一定程度收斂，這可根據(jù)模型階次的大小進(jìn)行相應(yīng)調(diào)整，狀態(tài)空間階數(shù)較小時 ρ0可以更小一點，狀態(tài)空間階數(shù)較大時 ρ0可適當(dāng)放寬。

4) 全局收斂標(biāo)準(zhǔn) ρg： ρg為不能降階時EM算法的最終迭代收斂標(biāo)準(zhǔn)，其可在 ρ0的基礎(chǔ)上進(jìn)行適當(dāng)降低，一般降低一個數(shù)量級即可。

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

3 應(yīng)用舉例

下面以數(shù)值模擬和現(xiàn)場實測的結(jié)構(gòu)振動響應(yīng)分析來對本文提出的方法進(jìn)行驗證說明。

3.1 6自由度結(jié)構(gòu)數(shù)值模擬分析

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

圖2 6自由度結(jié)構(gòu)模型簡圖Fig.2 Schematic diagram of 6-DOF structural model

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

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

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

圖4 6自由度結(jié)構(gòu)EM自降階迭代過程Fig.4 EM iteration process with auto order reduction of 6-DOF structure

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

表1 6自由度模態(tài)分析結(jié)果Table 1 Modal analysis results of 6-DOF structure

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

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

圖6 6自由度結(jié)構(gòu)質(zhì)點2處(10 s～20 s)加速度響應(yīng)分析Fig.6 Analysis of acceleration response at mass No.2 of 6-DOF structure (10 s-20 s)

圖7 6自由度結(jié)構(gòu)質(zhì)點2處(10 s～20 s)各階模態(tài)響應(yīng)及貢獻(xiàn)Fig.7 Modal response and contribution ratio at mass No.2 of 6-DOF structure (10 s-20 s)

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

3.2 廣州塔實測振動分析

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

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

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

圖9給出了部分傳感器(傳感器6、12和16)在18:00～19:00的加速度響應(yīng)曲線，所有傳感器獲得的加速度響應(yīng)的奇異值頻譜圖如圖10所示?？梢园l(fā)現(xiàn)3 Hz～5 Hz結(jié)構(gòu)的模態(tài)信息成分非常復(fù)雜，不易辨認(rèn)。

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

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

圖10 廣州塔加速度響應(yīng)奇異值頻譜圖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 廣州塔識別模型的響應(yīng)估計殘差奇異值頻譜圖Fig.12 Singular value spectrum of residuals of the identified model for Canton Tower

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

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

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

圖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 結(jié)論

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

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

(2) 提出了表征各階模態(tài)對總響應(yīng)的重要性的模態(tài)貢獻(xiàn)指標(biāo)，以此作為模型降階標(biāo)準(zhǔn)較之振型參與系數(shù)等指標(biāo)更具有物理意義和有效性。

(3) 提出的方法中一些閾值的選擇，仍帶有一定的經(jīng)驗性，需結(jié)合數(shù)據(jù)和結(jié)構(gòu)特征進(jìn)行調(diào)整。

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

(5) 通過模擬結(jié)構(gòu)數(shù)據(jù)和廣州塔實測數(shù)據(jù)的分析，表明本文的方法具有良好的準(zhǔn)確率和穩(wěn)健性。