基于改進確定-隨機子空間算法的橋梁結構模態參數識別

陳永高，鐘振宇,2，何 杰

(1. 浙江工業職業技術學院 建筑工程學院，浙江 紹興 312000； 2. 浙江大學 建筑工程學院，杭州 310058；3. 中鐵二院工程集團有限責任公司，成都 610031)

隨機子空間算法[1](Stochastic Subspace Identification,SSI)作為現階段最為廣泛使用的模態參數識別算法之一，能有效識別出環境激勵[2]下橋梁結構的固有頻率、阻尼比及模態振型。隨著該算法的不斷推廣，其缺陷也逐漸凸顯；雖然已有不少學者[3]對SSI存在的缺陷提出了相應的改進算法，但由于該識別系統的輸入信號為環境激勵(信號未知)，以致在模態參數的整個求解過程中均無法建立明確的激勵-系統-響應模型，導致系統特征值解空間的擴大，使得最終的識別結果中存在模態密集和虛假模態的現象，降低了識別結果的精確度。

確定-隨機子空間系統識別(Combined Determine-Stochastic Subspace Identification，CDSI)算法[4]是在SSI算法理論基礎上建立的，CDSI是利用結構的輸入信號(激勵)和輸出信號(響應)作為系統的“輸入”進行參數識別，能夠有效地彌補SSI算法的不足之處，提高參數識別的精確度。同時由于CDSI算法的特點在于能夠考慮結構的輸入、輸出及噪聲，所以得到了國外不少學者的青睞。Belleri等[5]將CDSI運用于混凝土框架結構，實現了地震激勵下不同損傷狀態下模態參數的準確識別；Bakir[6]將CDSI運用于識別地震激勵下的框架結構，識別結果表明該算法不僅能有效識別出結構的固有頻率，識別的模態振型同樣具有可靠性；Huang等[7]將CDSI運用于多層鋼框架振動臺試驗中，并通過對模態參數變化情況的辨識實現了對結構損傷狀態的判定。

雖然SSI算法被廣泛運用于橋梁結構的模態參數識別中，但CDSI算法卻幾乎沒被運用于識別橋梁結構的模態參數。主要原因歸于：CDSI的核心算法是線性算法[8]，并不適用于像橋梁結構這種非線性時變結構?；诖?，提出將滑窗技術[9]嵌套于CDSI算法中，以實現對橋梁結構進行時變模態參數跟蹤識別；同時為了避免人為參與真假模態的區分，提出了將模糊C均值聚類算法(Fuzzy C-Means，FCM)[10]運用于模態參數結果的智能化篩選，以完成模態參數的在線跟蹤智能化識別。

1 CDSI算法基本原理

圖1是CDSI算法的基本流程圖，其中：uk為系統的輸入信號；B為輸入矩陣；wk為環境因素帶來的各種過程噪聲；xk為系統狀態向量；A為狀態矩陣；C為輸出矩陣；vk為因傳感器帶來的測量噪聲；D為中途傳遞矩陣；yk為系統的輸出數據。

圖1 確定-隨機系統框圖Fig.1 Combined determine-stochastic system flowchart

得到確定-隨機狀態方程如式(1)所示

(1)

步驟1根據輸入信號和輸出信號建立Hankel矩陣，其中：Up作為過去的輸入數據；Uf作為將來的輸出數據。

(2)

步驟2確定卡爾曼濾波狀態序列。

Kk-1=(G-APk-1CT)(Λ0-CPk-1CT)-1

Pk=APk-1AT+(G-APk-1CT)×

(Λ0-CPk-1CT)-1(G-APk-1CT)T

(3)

式中：K為Kalman濾波增益矩陣；P為向前狀態矩陣。

(4)

得到卡爾曼濾波狀態序列為

(5)

步驟3建立含有狀態矩陣(A)、輸入矩陣(B)、輸出矩陣(C)及中途傳遞矩陣(D)的線性方程式。

(6)

步驟4對矩陣A進行特征值分解。

(7)

步驟5基于式(8)求解結構的固有頻率值(f)、阻尼比(w)及模態振型(Φ)。

(8)

基于步驟1～步驟5便能實現橋梁結構的模態參數識別，其中就如何確定狀態序列、矩陣A的特征值分解以及各步驟中各參數的具體含義可見文獻[11]。

2 滑窗技術的運用

2.1 加窗原理分析

對響應信號進行模態參數識別時需控制信號的長短以提高參數識別的效率，實際運用中可通過對信號作加窗截斷處理以控制信號的長短。對信號作加窗處理時，因窗函數會直接響應信號的頻譜，且不同的窗口形式產生的頻譜泄漏情況不同，以致最終的頻率分辨效果不同。基于此，如何選擇窗函數的種類、窗口大小及滑窗步長顯得尤為重要，以下將依次詳細論述如何確定這三項參數。

2.2 窗函數選擇

現階段常用的窗函數[12]包括矩形窗、漢寧窗、海明窗、平頂窗、高斯窗等。對比各窗函數的優缺點，最終選擇采用矩形窗作為本文響應信號的加窗函數。原因如下：矩形窗作為時間變量的零次冪窗，具有主瓣集中的優點；同時相比其它窗函數，矩形窗能精確識別出頻率結果，但識別的幅值結果精確度不高。在實際運用中，如果僅需精確地識別出主瓣的頻率，而無需在乎幅值的識別精度，比如識別結構的自振頻率，則可以優先選擇矩形函數對信號進行加窗處理。矩形窗窗函數的表達式為

w(t)=1

(9)

2.3 窗口大小確定

每次輸入系統的初始運算數據點個數-窗口大小會在一定程度上影響識別結果的精度。當窗口太小時，會因輸入信號不足導致無法識別出準確的穩定圖，且容易造成模態遺漏現象。當窗口太大時，會因輸入信號過多導致穩定圖中存在過多的虛假模態，且會降低參數識別的效率。為了確定窗口的大小，葉飛等[13]認為應該綜合考慮結構自身的特點、信號的采樣頻率及信號含噪程度等因素。但筆者認為根據這些因素來綜合考慮窗口的大小不具可行性，首先是因為涉及的因素太多，其次是因為各結構都具有自身的獨有特性?；诖?，可以借鑒穩定圖定階的原理來定義窗口的大小，具體分析步驟如下：

步驟1基于隨機子空間法中的穩定圖定階[14]確定結構的真實階次為N。

步驟2假定初始窗口數據長度為L1，對應窗口時間為t1s，采樣頻率為cf，三者之間滿足L1=t1×cf。

步驟3假定前后窗口時間差均為Ts，即第i個窗口對應的數據長度為Li=[t1+(i-1)×T]×cf。

步驟5分析BFi數值大小隨信號長度增加的變化情況，可知該百分比呈現先逐漸增加最后保持平穩的變化趨勢，如圖2所示。認為當百分比超過90%時對應的信號長度為合適的窗口大小。

圖2 穩定點百分比趨勢圖Fig.2 Percentage trend chart of stable points

2.4 窗口步長

窗口大小確定后，則需進一步確定窗口的滑動步長，即每次窗口向前推進的數據點數，實現對輸入信號的實時更新。假定信號每次更新的數據量為V，則各窗口間數據的更替、系統的輸入及輸出間的關系圖如圖3所示。

圖3 滑窗步長示意Fig.3 Schematic diagram of step length

根據Hankel矩陣的定義[15]可知，過去的輸入信號長度等于將來的輸出信號長度，基于此，將定義為1/2窗口對應的的數據量，即窗口間的重疊部分為1/2窗口。

3 FCM聚類算法的運用

無論是SSI算法還是CDSI算法均涉及人為篩選穩定圖中的有效模態[16]，為了避免因人為參與有效模態辨識而帶來的主觀性，提出了基于FCM聚類算法的有效模態智能篩選法。以下將首先介紹FCM聚類算法的基本原理及步驟，其次介紹如何將FCM算法和CDSI算法進行結合以實現有效模態的智能化辨識。

3.1 FCM聚類算法的基本原理

FCM聚類算法[17]相比K-means聚類和譜系聚類等而言，其優點在于：該算法的聚類原理是依據數據點的隸屬度來辨識其分類情況，是一種模糊性聚類算法[18]；且能實現不同維度數據間的聚類。

利用FCM算法將數據集X分為c類，假定X由N組樣本構成，即X={x1,x2,…，xN}；其中每個樣本由n維向量構成，即xj={xj1,xj2,…,xjn}。定義數據集X的隸屬度矩陣為U=[uij]，U為C×N的矩陣，其中C為分類個數，uij為X中的第j個樣本歸屬于第i個類的隸屬度。uij需滿足條件

(10)

其中,uij的計算公式為

(11)

式中：m為模糊指數，m∈[1,∞]；dij為數據對象xj和第i個聚類中心(vi)間的歐氏距離,dij∈‖xj-vi‖；{v1,v2,…，vi}為聚類中心集合，其計算公式為

(12)

FCM算法的目標函數值為J(u,v)，其迭代過程的目的即是找尋到使J(u,v)函數為最小值的隸屬度矩陣和聚類中心。

(13)

FCM聚類算法的具體步驟如下：

步驟1確定分類數目c和最大迭代次數m(模糊指數)。

步驟2從數據集X中隨機選取一組數據為初始聚類的中心點。

步驟3根據式(11)求解隸屬度矩陣U。

步驟4根據式(12)求解出聚類中心v。

步驟5根據式(13)找尋使J(u,v)函數為最小值的隸屬度矩陣和聚類中心，并完成聚類。

3.2 FCM與CDSI的結合

為了實現穩定圖中有效模態的智能化辨識，提出將FCM算法與CDSI算法進行結合以得到新的參數識別算法(C-CDSI)，該算法的具體步驟如下所示：

步驟1利用CDSI算法識別得到各窗口對應的參數結果Xi={Fi,wi,Φi}，其中F為頻率值，w為阻尼比，Φ為模態振型，各參數結果均為m×n的矩陣，其中m數值大小等于1/2系統真實階次。

Fi=fi(m,n)
wi=ξi(m,n)
Φi=ψi(m,n)

(14)

步驟2頻率聚類。

對Xi={Fi,wi,Φi}和Xj={Fj,wj,Φj}分別在頻率、阻尼比及振型三項參數上進行FCM聚類處理——①確定聚類數據集，假定對f1(m,n)與f2(m,n)中的兩兩頻率模態進行聚類；②確定聚類分類數，因每次聚類的數據集只有2階頻率模態，則設定分類數c=2；③確定最大迭代次數，因參數結果矩陣的行數均為1/2系統真實階次(N)，則設定最大迭代次數m=N/2；④首先對f1中第1階頻率結果f1(m,1)和f2中的第1階頻率結果f2(m,1)進行聚類，當兩者為同一類時，則取平均值作為新的第1階頻率結果，即f1-2(m,1)=(f1(m,1)+f2(m,1))/2；依次類推完成兩頻率矩陣中所有模態間的聚類，并將所有聚類后的新模態與未聚類的模態構成新的頻率矩陣f1-2。

步驟3阻尼比聚類。

為了驗證頻率矩陣f1-2中聚類的2階模態具有真實性，則對ξ1(m,n)和ξ2(m,n)中相應的2階模態進行阻尼比聚類辨識——①假定f1-2中f1(m,1)和f2(m,1)是同類項，則基于步驟2的算法流程完成對ξ1(m,1)和ξ2(m,1)的聚類；②當兩者為同類項時，取平均值作為新的1階阻尼比結果,ξ1-2(m,1)=(ξ1(m,1)+ξ2(m,1))/2；③當兩者不為同類項時，則將f1-2(m,1)拆分為最初的f1(m,1)和f2(m,1)，并重新構建新的頻率矩陣ff1-2；④將所有聚類后的阻尼比與未聚類的阻尼比構成新的阻尼比矩陣ξ1-2。

步驟4模態振型聚類。

為驗證ff1-2中頻率聚類和ξ1-2中阻尼比聚類的結果具有真實性，則可對ψ1(m,n)和ψ2(m,n)中相應2階模態進行模態振型聚類——①假定ff1-2中f1(m,1)和f2(m,1)是同類項，則基于步驟2中頻率聚類的步驟對ψ1(m,1)和ψ2(m,1)進行振型聚類辨識；②當兩者為同類項時，取平均值作為新的1階模態振型結果：ψ1-2(m,1)=(ψ1(m,1)+ψ2(m,1))/2；③當兩者不為同類項時，則將f1-2(m,1)分為最初的f1(m,1)和f2(m,1)，并重新構建新的頻率矩陣fff1-2；同時將ξ1-2(m,1)分為最初的ξ1(m,1)和ξ2(m,1)，并重新構建新的阻尼比矩陣ξξ1-2；④最后將所有聚類后的振型與未聚類的振型構成新的振型矩陣ψ1-2。

步驟5將fff1-2,ξξ1-2和ψ1-2組建為新的一組參數結果X1-2={fff1-2,ξ1-2,ψ1-2}，并基于步驟2～步驟4完成X1-2與X3間的同類模態聚類，得到新的聚類結果X2-3；以此類推，依次完成所有窗口間的參數結果聚類，假定窗口總數為N，則最后的聚類參數結果為XN-1-N。

步驟6為篩選出XN-1-N中有效模態(頻率、阻尼比及振型均穩定的模態)，可統計XN-1-N中各聚類模態是由多少窗口中相應模態聚類而成的；并將聚類窗口數大于0.8N的模態作為有效模態繪制于穩定圖中。

綜合上述可知，可將滑窗技術、FCM聚類算法及CDSI識別算法融匯到一起構建新的模態參數在線智能化識別算法(SC-CDSI)，該算法的具體流程圖如圖4所示。

圖4 SC-CDSI識別算法流程圖Fig.4 Flowchart of SC-CDSI

4 振動臺試驗橋

為了驗證前文所提改進確定-隨機子空間算法(SC-CDSI)的可靠性，以某大型斜拉橋振動臺試驗為背景，分析其在地震激勵下的模態參數變化情況。

4.1 工程概況

該試驗橋為斜拉橋，按照1∶20的比例對實際橋梁進行縮小，實驗室中橋梁跨度為(6.5+19+6.5) m。圖5為該橋梁結構的橋型布置圖、主梁的平面圖，其中索塔和橋墩采用M15微?；炷聊M；普通鋼筋采用直徑6 mm的圓鋼模擬；抗剪鋼筋采用10#鋼絲模擬；主梁和鋼箱梁均采用厚度為5 mm的鋼板模擬；拉索采用直徑為10 mm的鋼絲繩模擬。

4.2 加速度傳感器位置

該振動臺試驗，在基礎底部共安裝4處加速度傳感器以監測結構的輸入信號；在主梁上共安裝11處加速度傳感器以收集結構的響應信號。圖6為傳感器的布置圖。

4.3 地震輸入

輸入激勵的各項參數見表1所示，可知峰值加速度均為0.1g，地震波的方向均為橫橋向，波形分別為Chichi波和場地波。圖7為C1工況下Chichi地震波的加速度時程曲線，傳感器的采樣頻率均為256 Hz。

圖5 橋梁模型立面圖及平面圖(cm)Fig.5 Plans and elevations of bridge model(cm)

圖6 加速度傳感器布置圖(cm)Fig.6 Arrangement of acceleration sensors(cm)

表1 試驗工況表Tab.1 Working conditions of test table

圖7 加速度時程曲線(C1)Fig.7 Time-curve of acceleration(C1)

4.4 MIDAS有限元結果

采用MIDAS 2017建立該橋梁結構的仿真分析模型。橋墩、索塔以及主梁均采用非線性梁柱單元模擬，拉索采用桁架單元模擬；支座均采用盆式橡膠支座，且采用零長度單元模擬。全橋模型如圖8所示。對該橋作特征值分析，獲得其前5階頻率值，如表2所示，圖9為前4階模態振型圖。

圖8 全橋模型(MIDAS)Fig.8 Bridge model(MIDAS)

表2 自振頻率Tab.2 Natural frequency Hz

圖9 前4階模態振型圖(MIDAS)Fig.9 Modal pattern of the first four orders (MIDAS)

4.5 各窗口參數識別

根據滑窗技術和CDSI算法識別該斜拉橋在兩種地震激勵下的模態參數，設定初始窗口時間為1 s，試算出合適窗口大小為4 s；同時對應的滑窗步長為2 s，即每隔2 s更新一次窗口數據。對各窗口的信號進行參數識別，得到兩種工況下各窗口對應的穩定圖結果，圖10為C1工況下前兩幅穩定圖結果。

圖10 穩定圖(C1工況)Fig.10 Stable diagram (working condition: C1)

識別出C1和S1兩種工況下頻率值在各窗口中的具體數值并繪制前5階頻率值隨時間的趨勢圖，結果如圖11所示。分析圖中各階頻率的變化情況，可知兩種工況下橋梁結構在前5階的頻率值均處于穩定的狀態，即橋梁結構在兩種地震波作用下，其自身的固有頻率具有穩定性，也表明該斜拉橋在兩種工況下未產生損傷。

圖11 頻率走勢圖Fig.11 Frequency chart

4.6 有效模態智能辨識

基于3.2節所提FCM聚類算法對兩種工況下各窗口參數結果進行聚類分析，實現穩定圖中有效模態的智能化辨識，圖12為C1和S1工況下對應的有效穩定圖。

圖12 有效穩定圖Fig.12 Effective stable diagram

將MIDAS有限元結果、C1和S1工況下前5階頻率結果作對比分析，結果如表3所示。根據表中兩兩結果間的差值百分比，可知SC-CDSI算法所得有效模態頻率值與有限元結果間的誤差在7%以內，即所提SC-CDSI算法的頻率識別結果具有可靠性。

表3 頻率對比分析結果Tab.3 Comparative analysis results of frequency

4.7 前3階模態振型結果

為驗證SC-CDSI算法不僅能精確地篩選出穩定圖中的真實模態，還能識別出結構的模態振型。得到該斜拉橋的前3階振型圖，如圖13所示，將其與圖9中的前3階振型圖進行對比，可知，SC-CDSI算法所得模態振型與有限元結果具有90%以上的相似度，即SC-CDSI算法所得振型圖具有可靠性。

圖13 前3階振型圖Fig.13 Modal pattern of the first three orders

5 結 論

考慮到確定-隨機子空間算法并不適用于橋梁結構這種時變結構，提出了將滑窗技術運用于CDSI算法中實現對橋梁結構模態參數的實時跟蹤識別；同時考慮到需人為參與穩定圖中有效模態的辨識，提出了將FCM運用于辨識穩定圖中的真假模態，以實現模態參數結果的智能化篩選。

最后以某大型斜拉橋振動臺試驗為背景，將所提SC-CDSI算法運用于識別該橋的模態參數，并將所得結果與MIDAS有限元結果進行對比分析，結果表明，所提SC-CDSI算法不僅能精確地識別出橋梁結構的固有頻率值，還能精確地識別出結構的模態振型結果，能夠實現試驗橋梁結構模態參數的在線智能化跟蹤識別。因缺乏實際地震的輸入信號和輸出信號，所以該算法在實橋中的運用還需進一步分析論證。