基于K-means聚類算法的橋梁結構真實模態篩選研究

尹紅燕，劉東霞，唐 莉

(重慶交通職業學院，重慶 402247)

0 引言

隨著我國橋梁結構的飛速發展，各種大型橋梁結構應運而生，隨之而來的弊端是橋梁安全事故的發生頻率日益增長，目前如何判別橋梁結構在運營階段的健康狀態[1]已經成為了一個十分重要的研究課題。實際工程中，常通過對橋梁結構進行動力監測[2]來實現對其的安全監測，因為一旦結構發生損傷[3]，則其自身的動力特性會發生不同程度的變化。基于此，提出了通過分析橋梁結構的固有頻率、阻尼比以及模態振型[4]的變化情況來判別其健康狀態。

隨著橋梁跨度的不斷增大、新型材料的使用以及環境因素的影響導致現有的模態參數識別算法[5]已不能很好地滿足人們對結果精度的要求。隨機子空間算法(Stochastic Subspace Identification-SSI)[6-7]作為主要識別算法之一，已有不少學者研究了如何辨識穩定圖中的真假模態，包括湯寶平等[8]通過在參數識別過程中引入譜系聚類算法來完成模態參數的自動化識別；UBERTITN F等[9]利用譜系聚類分析法來辨識穩定圖中的真假模態；章國穩等[10]將聚類算法運用于模態參數結果的篩選過程中，以辨識穩定圖中的真假模態。

雖然這些算法能在不同程度上實現真實模態的智能化識別，但識別效果均不佳，其主要原因在于：在參數識別過程中忽略了結構自身真實模態[11]存在的一般規律。基于此，本研究以某試驗橋梁為研究對象，首先通過分析大量穩定圖中真假模態的存在形式來總結其存在的一般規律，其次基于K-means算法[12]提出了一種新的模態參數識別算法，最后將該算法運用于試驗橋和實際橋梁結構中以驗證算法的可信性。

1 DATA-SSI算法

隨機子空間算法之所以被廣泛運用于橋梁結構的模態參數識別中，其主要原因在于該算法相比其他參數識別算法而言有如下3方面的優點：(1)可直接將時域信號作為輸入數據；(2)不會出現頻率分辨率誤差等問題；(3)識別的固有頻率、阻尼比以及模態振型均能滿足實際工程的需求。

該算法共分兩類，分別是基于數據驅動的隨機子空間算法[13](Data driven Stochastic Subspace Identification,DATA-SSI)和基于協方差驅動的隨機子空間算法[14](Covariance driven Stochastic Subspace Identification,COV-SSI)。對比這兩種算法的基本原理和計算效率可知：

(1)DATA-SSI算法在數據縮減方面具有更高的效率，因為該算法是通過正交投影的方式對 Hankel矩陣中的數據進行縮減；而COV-SSI算法則是利用協方差的方式進行數據的縮減。

(2)DATA-SSI算法在計算效率方面具有更高的效率，因為該算法的計算過程中僅涉及一次QR分解；而COV-SSI算法過程中需對Toeplitz進行多個協方差計算。

鑒于上述原因可知：DATA-SSI算法更適應于實現模態參數的自動化識別，所以本研究以該算法為研究對象進行模態參數自動化識別研究。以下將簡單介紹基于數據驅動的隨機子空間算法的基本理論。

1.1 基本原理

DATA-SSI算法的大致計算步驟如下：

Step1：利用公式(1)建立Hankel矩陣，Yp為“過去”信息，Yf為“將來”信息；

(1)

Step2：QR分解

利用公式(2)對Hankel矩陣進行QR分解；

(2)

Step3：奇異值分解

(3)

Step4：特征值分解

利用公式(4)對系統矩陣A進行特征值分解；

A=eAcΔt，

(4)

Step5：模態參數計算

利用公式(5)計算模態參數結果；

Φi=CiΨi，

(5)

1.2 穩定圖定階

在利用DATA-SSI算法進行參數識別時，需事先明確系統的真實階次，階次過高過低都會導致識別的結果不可信，因為過高時會伴隨虛假模態的產生；過低時會伴隨模態遺漏現象的發生。基于此，可利用穩定圖法[15]對系統進行真實階次的確定，該方法的基本原理是以頻率為橫坐標，系統階次為縱坐標，將各階次對應的極點繪于同一圖中，并分析前后兩階對應各參數是否滿足公式(6)所示容差，當滿足時則可判定該階次為系統的真實階次。

(1-MAC(j,j+1))×100%<Δψ(振型容差)，

(6)

式中，j為階次數；f，ξ，ψ分別為頻率值、阻尼比及振型；MAC為模態置信因子。關于容差的具體取值，文獻[16]對其進行了深入研究，根據該文獻可將頻率容差定為1%；振型容差和阻尼比容差為5%。圖1為穩定圖定階的一般流程圖。

圖1 穩定圖定階原理Fig.1 Principle of stable diagram order determination

2 橋梁結構的真實模態

為了分析橋梁結構真實模態存在的一般規律，以某大型斜拉橋振動臺試驗為研究對象，首先將各節點各時間段內的振動信號作為DATA-SSI算法的輸入完成參數識別；并對比分析大量穩定圖中模態參數的存在形式，進而總結出穩定圖中真假模態的存在規律。

2.1 試驗橋概況

該橋梁結構為斜拉橋，其跨度為(130+380+130)m，試驗橋梁是以1∶20的比例對實橋進行比例縮放。橋型布置圖、主梁的平面圖如圖2所示。橋梁結構各部分對應的材料情況如下：

(1)采用M15微粒混凝土模擬索塔和橋墩；

(2)采用直徑為6 mm的圓鋼模擬普通鋼筋；

(3)采用10#鋼絲模擬抗剪鋼筋；

(4)采用厚度為5 mm的鋼板模擬主梁結構和鋼箱梁；

(5)采用直徑為10 mm的鋼絲繩模擬拉索。

2.2 傳感器布置

橋梁結構建立于振動臺上，試驗過程中可通過振動臺對其施加不同的激勵，并利用事先布置在主梁上的加速度傳感器采集各點對應的振動信號，傳感器布置圖如圖3所示，主梁上共設置11個傳感器。

圖2 橋梁實景圖(單位：cm)Fig.2 Real images of bridge(unit:cm)

圖3 加速度傳感器布置圖(單位：cm)Fig.3 Arrangement of acceleration sensors (unit: cm)

2.3 地震波及響應信號

輸入地震波分為兩類，分別是Chichi波和場地波(Sit波)，各地震波對應的輸入方向、輸入波形以及峰值加速度如表1所示。根據該表可知，根據地震波波形和峰值加速度的不同將振動臺輸入地震波分為6種工況。各工況下信號的采樣頻率均為256 Hz，采樣時間為240 min。圖4是G1工況下6#傳感器在最開始20 s內對應的加速度響應信號時程圖。

表1 試驗工況表

2.4 真實模態的一般規律

以每分鐘采集的響應信號為DATA-SSI算法的輸入數據，即輸入信號的矩陣大小為15 360×11(11代表共11個傳感器)，即同一地震波可識別得到720幅穩定圖，圖5僅為G1工況下隨機的3幅穩定圖。

對比分析6種工況下共1 440幅穩定圖，得到如下3點結論：

(1)同一工況的不同時間段內，穩定圖中的穩定軸(真實模態)具有穩定性，即不會因為虛假模態的存在發生較大差異；

圖4 加速度時程曲線Fig.4 Curve of acceleration vs. time

圖5 穩定圖Fig.5 Stability diagrams

(2)同一工況的不同時間段內，穩定圖中的虛假模態具有不穩定性，即各穩定圖中虛假模態的分布情況具有離散性；

(3)兩種地震波識別所得真實模態具有互補性，即如果僅對試驗橋進行一種地震波試驗，則識別的真實模態會出現模態遺漏現象。

綜合上述：對于某一固定的結構而言，其自身的動力特性參數具有時不變性，當該結構受到外界因素的影響而發生損傷時，其動力特性參數(真實模態) 會隨之發生一定的變化。

3 真實模態的篩選

人為辨識穩定圖中的真假模態不僅會耗費大量時間，還會導致辨識的結果不可信。基于此，可以利用多元統計學中的聚類算法[17-18]來從大量穩定圖中篩選真實模態和虛假模態，通過引入K-means聚類算法的方式來辨識真假模態。之所以選擇K-means聚類算法，是因為該聚類算法相比譜系聚類算法而言，不僅能實現不同維度數據的聚類，也能判別數據間是否屬于同一類，還能篩選出同類的數據點。以下將詳細介紹K-means算法的基本原理以及如何將其融入到模態參數識別的流程中。

3.1 K-means算法的基本原理

K-means算法的基本原理：按照樣本之間的距離大小將樣本集劃分為K個簇，使得各簇內的點都盡可能得緊密聚集，進而保證簇間的距離盡可能得大。假定輸入樣本集D={x1,x2,…xm}，聚類的簇樹為K，最大迭代次數為N，輸出是簇劃分C={C1,C2,…,Ck}。

該算法的詳細步驟如下：

Step1：從樣本集D中任意選擇K個樣本作為最初的K個質心向量：{u1,u2,…uk}。

Step2：對于各迭代次數n=1,2,…N。

(1)將簇劃分C初始化為Ct=?，t=1,2,…k；

(2)對于各樣本i=1,2,…,m，計算樣本xi和各個質心向量uj(j=1,2,…,k)的距離：

(7)

標記最小的dij所對應的類別為λi，并更新Cλi=Cλi∪{xi}；

(3)對于各簇樹j=1,2,…,k，重新計算Cj中所有樣本點的新質心uj，

(8)

(4)當K個質心向量均未變化時，則進行Step3。

Step3：輸出簇劃分C={C1,C2,…,Ck}。

3.2 真實模態的篩選

將K-means算法運用于辨識多幅穩定圖中的真實模態，具體步驟如下：

Step1：基于1.2節所提穩定圖定階算法確定結構的真實系統階次N；

Step2：利用1.1節所提DATA-SSI算法識別得到結構在多個時間段內的參數結果，假定參數結果共有NN組。通過不斷的反復試驗，建議NN≥300。定義每組參數結果Xi={Fi,Di,Mi}(i=1,2,…,NN)；

Fi=fi(m,n)，

Di=ξi(m,n)，

Mi=ψi(m,n,c)，

m=1,2,…,N，

n=1,2,…,j，

c=1,2,…,g，

(9)

式中，fi為第i組參數結果中的頻率二維數據矩陣(m×n)；ξi為第i組參數結果中的阻尼比二維數據矩陣(m×n)；ψi為第i組參數結果中的阻尼比三維數據矩陣(m×n×c)；m為矩陣的行數，等于系統的階次N；f為傳感器個數，即輸入信號的例數。

Step3：基于K-means算法實現第1組參數結果X1={F1,D1,M1}與第2組參數結果X2={F2,D2,M2}間的同類項聚類，步驟如下：

(1)頻率矩陣f1與頻率矩陣f2間的聚類

基于K-means算法對f1中的第i列頻率向量f1(m,i)(m=1,2,…,N)與f2中的第j列頻率向量f2(m,j)(m=1,2,…,N)進行聚類，由于頻率矩陣中存在為0的項，則需首先剔除向量中的0項，得到重構的向量f′1(m,i)(m=1,2,…,N)和f′2(m,j)(m=1,2,…,N)。

基于以上步驟實現頻率矩陣f1與f2的頻率聚類，即將所有的同類模態都聚為一類，同時將所有的聚類模態與未聚類的模態組成新的參數結果f1-2。

(2)利用阻尼比進一步驗證頻率聚類結果

假定第1組頻率參數結果的第1列模態f1(m,1)和第2組頻率參數結果的第1列模態f2(m,1)為同一模態，則基于頻率聚類的原理驗算兩者對應的阻尼比結果是否屬于同類項，即驗證ξ1(m,1)(m=1,2,…,N)和ξ2(m,1)(m=1,2,…,N)所得的簇劃分結果是否滿足95%的原則。當滿足上述條件時，則認為f1(m,1)和f2(m,1)為同類模態予以保留，反之則剔除這組聚類結果。

(3)基于穩定圖法的原理驗證步驟(1)和步驟(2)所得的同類模態

(10)

基于步驟(1)～(3)完成第1組參數結果X1={F1,D1,M1}與第2組參數結果X2={F2,D2,M2}間的所有聚類，將聚類后的模態和未聚類的所有模態整合到一起，并定義為新的參數結果X1-2。

Step4：基于Step3完成X1-2與X3間的同類模態聚類，得到新的聚類參數結果X2-3；以此類推，完成NN組參數結果的同類聚類，得到X(NN-1)-NN。

Step5：由于X(NN-1)-NN中的各階模態有的是從X1到XNN中通過聚類得到的，有的則是非聚類的模態。為了進一步篩選出最為真實的模態，可首先統計X(NN-1)-NN中各聚類模態的聚類項數(即：該聚類模態是由多少組模態聚類組成的)；其次選擇聚類項數大于0.85NN的模態，并將其繪制于最終的穩定圖中；同時為了使識別所得的模態參數包含所有聚類項，所以取同類模態的均值作為最終的穩定軸繪制于穩定圖中，進而得到最終的聚類穩定圖。

梳理上述所有步驟，可得真實模態辨識流程，如圖6所示。

圖6 真實模態辨識流程Fig.6 Flowchart of real mode identification

4 試驗橋參數識別

4.1 有限元結果

采用MIDAS軟件建立2.1節所提試驗橋梁結構的三維有限元模型，其中橋墩、索塔以及主梁對應的截面均采用纖維截面，且單元采用非線性梁柱單元(Nonlinear Beam-Column Elements)，采用桁架(Truss)單元模擬拉索；支座為盆式橡膠支座，且利用零長度單元(Zero-Length Element)模擬；混凝土的本構采用Concrete02 Material材料，鋼筋本構采用Steel02材料。全橋模型如圖7所示。

經特征值分析可獲得其固有階頻率值和模態振型圖，本研究僅僅羅列了其前7階頻率值和前4階模態振型，頻率值如表2所示，振型圖如圖8所示。

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

表2 自振頻率(單位：Hz)

圖8 前4階模態振型圖(MIDAS)Fig.8 Mode shapes of the first 4 orders (MIDAS)

4.2 識別結果

由于該橋梁結構的輸入地震波共分兩類，分別為Chichi波和場地波(Sit波)，信號采集時間為240 min，則以每分鐘對應的響應信號為參數識別的輸入數據，即每種地震波可得720幅穩定圖。基于第3節所提真實模態篩選算法對這1 440幅穩定圖進行辨識，繪制出最終的聚類穩定圖，如圖9所示。

圖9 聚類穩定圖Fig.9 Clustering stability diagram

將圖9中的各階頻率值與表2中MIDAS所得理論頻率值作對比分析，結果如表3所示。

表3 頻率對比

根據該表可得如下結論：

(1)本研究所提真假模態的一般規律具有可信性，且所提算法不僅能有效地辨識出橋梁的真實模態和虛假模態，還能有效避免模態遺漏現象的發生；

(2)識別所得頻率值與MIDAS有限元軟件所得理論值間的差值百分比在[-4.3%,4.96%]，驗證了所提算法具有可信性。

為進一步驗證識別所得模態振型同樣具有可信性，繪制了該試驗橋的前3階二維模態振型圖和三維模態振型圖，如圖10所示。將其與MIDAS有限元軟件特征值分析所得的前3階模態振型(圖8)作對比，結果表明：本研究算法識別的模態振型圖與理論模態振型圖具有較高的相似度。

圖10 前3階模態振型圖Fig.10 Mode shapes of the first 3 orders

之所以沒有對阻尼比的結果進行對比分析，主要是因為結構的阻尼比受環境因素的影響很大，且人們對其的認知水平還比較有限。同時在MIDAS模型計算中采用的假定與實際工程結構間存在一定的差異。基于此，本研究僅對比分析了試驗橋的固有頻率值和模態振型結果。

5 某大型斜拉橋參數識別

5.1 實橋工程概況

為驗證所提算法不僅能運用于試驗橋梁，還能適用于實際橋梁結構。以重慶市某大型懸浮體系斜拉橋為識別對象，其主跨為330 m，邊跨為149 m。其主梁上共設置有11處加速度傳感器，橋跨布置圖及傳感器位置圖如圖11所示，其中主梁上的圓圈為傳感器的大致布置位置。加速度信號的采樣頻率為1 s 80次，即頻率為80 Hz。圖12為主梁上某一傳感器在隨機50 s內采集到的信號時程圖。

圖11 橋跨布置圖(單位：cm)Fig.11 Layout of bridge spans(unit:cm)

圖12 加速度響應信號時程Fig.12 Time history of acceleration response signal

5.2 自振頻率理論值

在主梁的主跨跨中采用跳車激振的方式產生激勵，即通過汽車后輪越過帶有坡面的三角橫木(高度約12 cm)，利用車輪落下對結構的沖擊來激勵橋梁振動；并利用加速度傳感器收集結構自振信號，最后對脈動信號進行分析識別得到其自振頻率，跳車自振頻譜圖如圖13所示，圖中橫坐標為頻率值，縱坐標為對應的加速度值。

圖13 實測跳車自振頻譜圖Fig.13 Actual frequency spectrum of jumps off

5.3 參數識別結果

利用加速度傳感器采集斜拉橋在2019年4月份的所有加速度信號，并以每小時對應的響應信號作為參數識別的輸入數據，即4月份30天可識別得到720幅穩定圖。基于第3節所提真實模態篩選算法對這720幅穩定圖進行辨識，繪制出最終的聚類穩定圖，如圖14所示。

圖14 穩定圖Fig.14 Stability diagram

將穩定圖中前6階頻率值與跳車自振頻率譜圖中的頻率值進行對比分析，得到表4所示結果，分析表中數據可知，識別所得頻率值與跳車試驗所得頻率值間的差值百分比在[-4.7%,6.5%]，驗證了本研究算法的識別結果具有可信性。

表4 頻率結果對比

圖15 前3階模態振型圖(二維)Fig.15 Mode shapes of the first 3 orders (2D)

圖15和圖16分別是該斜拉橋的前3階二維模態振型圖和三維模態振型圖，將其與理論模態振型圖作對比分析，可知識別所得模態振型與理論振型圖具有較高的相似度。

圖16 前3階模態振型圖(三維)Fig.16 Mode shapes of the first 3 orders (3D)

6 結論

本研究首先基于某試驗橋結構分析了穩定圖中真假模態的一般規律，其次通過在參數識別的過程中引入“K-means聚類算法”來實現穩定圖中真假模態的智能化辨識,最后將所提算法運用于識別某試驗和某實橋的模態參數，并將結果與理論值作對比，結果表明：

(1)所提真假模態的存在規律具有可信性，即結構的真實模態在結構未受損傷時具有穩定性，而虛假模態則會隨著時間的變化發生變化，呈現的現象是在各穩定圖中具有不穩定性和離散性；

(2)K-means算法能夠實現多幅穩定圖中同類模態參數的聚類完成真假模態的智能化篩選；

(3)本研究算法不僅適用于試驗橋的模態參數識別，還適用于實橋結構的模態參數識別，且識別結果均具有較高的可信性。