一種自動識別結(jié)構(gòu)模態(tài)參數(shù)的隨機子空間方法張小寧 段忠東

張小寧 段忠東

摘要： 對結(jié)構(gòu)進行實時在線監(jiān)測的要求提出了對結(jié)構(gòu)模態(tài)參數(shù)進行自動識別的需求。目前發(fā)展的結(jié)構(gòu)模態(tài)參數(shù)識別方法均需要人工干預(yù)，為實現(xiàn)無人值守的結(jié)構(gòu)實時監(jiān)測目標，試圖發(fā)展一種結(jié)構(gòu)模態(tài)參數(shù)自動識別方法。隨機子空間法唯一需要確定的參數(shù)是“系統(tǒng)的階次”，因此，提出了一種基于頻率穩(wěn)定性和振型穩(wěn)定性自動判別系統(tǒng)階次的方法，基于此，建立了基于隨機子空間法的模態(tài)參數(shù)自動識別方法；通過兩個橋梁算例，對該方法的適用性和魯棒性進行了驗證。關(guān)鍵詞： 模態(tài)參數(shù)識別； 隨機子空間法； 斜拉橋； 健康監(jiān)測； 系統(tǒng)階次

中圖分類號：TU311.3； U441+.3文獻標志碼： A文章編號： 10044523（2017）04054207

DOI：10.16385/j.cnki.issn.10044523.2017.04.003

引言

隨著超高層建筑和大跨橋梁的興建，健康監(jiān)測系統(tǒng)越來越多的應(yīng)用于此類大型土木工程結(jié)構(gòu)中。這些結(jié)構(gòu)的損傷診斷和受損結(jié)構(gòu)評定、實時監(jiān)測和安全預(yù)警、有限元模型修正等成為了健康檢測系統(tǒng)的重點，而結(jié)構(gòu)的模態(tài)參數(shù)識別是需要首先解決的問題[13]。為了能夠及時掌握結(jié)構(gòu)的運行狀況，及時發(fā)現(xiàn)結(jié)構(gòu)損傷，提前做出安全預(yù)警，因此需要對結(jié)構(gòu)進行實時在線監(jiān)測，而這也對結(jié)構(gòu)模態(tài)參數(shù)自動識別提出了要求。

隨機子空間法（SSI）是近年來發(fā)展起來的模態(tài)參數(shù)識別的方法[45]。傳統(tǒng)的模態(tài)參數(shù)識別方法需耗費大量工作來確定具有最小量參數(shù)的模型即規(guī)范模型，而隨機子空間法是基于隨機狀態(tài)空間模型，只需確定“系統(tǒng)的階次”這一個參數(shù)；傳統(tǒng)方法都需要進行迭代運算，往往產(chǎn)生發(fā)散或者收斂緩慢等問題，而隨機子空間法由于算法中的QR分解，矩陣正交投影，奇異值分解等運算不需進行迭代運算，故隨機子空間法不存在收斂問題；傳統(tǒng)方法往往需要計算高階次的模型，而隨機子空間法通過協(xié)方差運算或矩陣正交投影使系統(tǒng)階次明顯降低，從而運算更加快捷[67]。

1基本原理〖2〗1.1傳統(tǒng)SSI定階方法隨機子空間法經(jīng)典的確定系統(tǒng)階次的方法主要有奇異值跳躍法和穩(wěn)定圖法。奇異值跳躍法是根據(jù)奇異值的跳躍性來確定系統(tǒng)的階次，奇異值的跳躍點就是真實非零奇異值與零奇異值的分界點，系統(tǒng)的階次就是奇異值跳躍點之前所有奇異值個數(shù)的一半。穩(wěn)定圖法假定系統(tǒng)有多個階次，計算每一階次系統(tǒng)的模態(tài)參數(shù)，對應(yīng)于某一階模態(tài)，將相鄰階次的模態(tài)參數(shù)作比較，如果頻率、阻尼比和振型的差異小于提前預(yù)設(shè)的限值，則該點就稱為穩(wěn)定點，穩(wěn)定點組成穩(wěn)定軸，從而組成了穩(wěn)定圖。穩(wěn)定軸對應(yīng)的階次為系統(tǒng)的階次。

這兩種方法都需要人工參與，且在人工識別過程中容易受到噪聲信號的干擾，很容易確定出錯誤的系統(tǒng)階次[8]。

1.2SSI方法系統(tǒng)階次自動確定原理

系統(tǒng)階次自動判別方法首先通過求取系統(tǒng)從1到nmax每個階次的模態(tài)參數(shù)，用求得的模態(tài)參數(shù)構(gòu)造下三角矩陣，然后判定頻率的穩(wěn)定性——模態(tài)置信因子（MAF）和振型的穩(wěn)定性——模態(tài)保證準則（MAC），當(dāng)MAF·∑MAC最大時，該階次就是系統(tǒng)的最佳階次，從而實現(xiàn)了自動定階。

1.2.1求取下三角矩陣

假定系統(tǒng)的階次依次取1，2，3，4，…，nmax，nmax的取值至少要大于系統(tǒng)的真實階次。求得每個階次對應(yīng)的系統(tǒng)狀態(tài)空間模型，此時系統(tǒng)矩陣的階數(shù)為2，4，6，8，…，2nmax。

隨機子空間法是基于隨機系統(tǒng)狀態(tài)空間模型的方法[10]，隨機系統(tǒng)狀態(tài)空間方程如下所示xk+1=Axk+wk

yk=Cxk+vk （1）式中A為系統(tǒng)的狀態(tài)矩陣；C為系統(tǒng)輸出矩陣；wk為k時刻的過程噪聲；vk為k時刻的測量噪聲；x（t）為系統(tǒng)的狀態(tài)向量；y（t）為系統(tǒng)的輸出向量。

本文采用的是基于數(shù)據(jù)驅(qū)動的隨機子空間方法（SSIDATA），首先將輸出數(shù)據(jù)直接組成Hankel矩陣；然后對Hankel矩陣進行QR分解，得到投影矩陣；對投影矩陣進行奇異值分解（SVD），得到擴展的可觀測矩陣和系統(tǒng)狀態(tài)的卡爾曼濾波；最后由卡爾曼濾波序列和系統(tǒng)輸出采用最小二乘法求得系統(tǒng)的狀態(tài)矩陣A和輸出矩陣C。對系統(tǒng)狀態(tài)矩陣A進行特征值分解A=ΦΛΦ-1（2）式中Λ為包含復(fù)特征值的對角矩陣； Φ為特征向量組成的矩陣。

上式中，Λ=diag[λi]，其特征值兩兩共軛。令λi=ai±jbi （3）可算得系統(tǒng)的振動頻率fi和阻尼比εi：fi = a2i + b2i /（2π）（4）

εi = -ai a2i + b2i （5）將所求模態(tài)參數(shù)構(gòu)造下三角矩陣，故最終求得了3個nmax階的下三角方陣Fnmax=f1，100…0

f2，1f2，20…0

f3，1f3，2f3，3…0

……………

fnmax，1fnmax，2fnmax，3…fnmax，nmax（6）式中Fnmax為 振動頻率的下三角矩陣。fa，b為 第a階次的第b階振動頻率。

同時也可以得到模態(tài)保證準則（MAC）和阻尼比的下三角矩陣。

1.2.2計算模態(tài)置信因子（MAF）

第4期張小寧，等：一種自動識別結(jié)構(gòu)模態(tài)參數(shù)的隨機子空間方法振 動 工 程 學(xué) 報第30卷歐氏距離也被稱為歐幾里得度量，是比較常用的對距離的定義，它表示n維空間中兩點之間的距離，對于n維空間，歐氏距離可以表示為d（x1，x2，…，xn）=[（x1-x、1）2+

（x2-x、2）2+…+（xn-x、n）2]12（7）式中d（x1，x2，…，xn）為n維空間的歐氏距離；xn為所求階次的第n階模態(tài)參數(shù)（頻率、阻尼比、MAC）。

在這里，定義一種專門針對于振動頻率穩(wěn)定性的因子——模態(tài)置信因子（MAF），模態(tài)置信因子在數(shù)值上等于歐氏距離的倒數(shù)（d（x1，x2，…，xn））-1：MAFa，b=[（xa，1-xb，1）2+（xa，2-xb，2）2+

…+（xa，n-xb，n）2]-12（8）式中MAFa，b為a階次模態(tài)（頻率）和b階次模態(tài)之間的模態(tài)置信因子；xa，b為a階次的第b階振動頻率或阻尼比；

因為隨機子空間法所求出的各階次的模態(tài)參數(shù)之間的差異先是隨著系統(tǒng)階次的升高而減小，差異達到最小后又會隨著系統(tǒng)階次的升高而增大，歐氏距離表示了相鄰階次模態(tài)參數(shù)之間的差異，歐氏距離和系統(tǒng)階次之間的關(guān)系如圖1所示。所以，MAFa，b越?。W氏距離越大）表示該階模態(tài)的置信度越低，MAFa，b越大（歐氏距離越?。┍硎驹撾A模態(tài)的置信度越高，越接近系統(tǒng)真實模態(tài)。圖1歐氏距離隨系統(tǒng)階次的變化

Fig.1The variation of Euclidean distance with the order of system

實際應(yīng)用中，往往只取系統(tǒng)的前幾階模態(tài)參數(shù)，式（8）中n的取值為所取模態(tài)參數(shù)的最高階數(shù)。計算出nmax-1個階次的歐氏距離d1，d2，…，dnmax-1，求出對應(yīng)的nmax-1個模態(tài)置信因子MAF2，1，MAF3，2，…，MAFnmax，nmax-1，然后選取最大的m個模態(tài)置信因子MAFi，i-1，MAFj，j-1，MAFk，k-1，…，其所對應(yīng)的階次就是系統(tǒng)的近似階次，建議m取至少大于3。

1.2.3模態(tài)保證準則（MAC）

在上一步驟中，選取了模態(tài)置信度最高的m個階次，這m個階次可以看作是頻率穩(wěn)定的階次，但它們都是系統(tǒng)的近似階次，可能含有噪聲模態(tài)的階次。為了避免噪聲模態(tài)，還應(yīng)該對這m個頻率穩(wěn)定的階次進行振型的判定，看其振型是否穩(wěn)定。最后選取這m個階次中頻率穩(wěn)定性和振型穩(wěn)定性最高的階次作為系統(tǒng)的最佳階次。振型穩(wěn)定性用模態(tài)保證準則（Modal Assurance Criterion）來定義MACa，b=φTa，jφb，j2（φTa，jφa，j）（φTb，jφb，j）（9）式中MACa，b為 兩組振型（理論和試驗、完好和損傷）的相關(guān)性振型；φa，j為 a階次第j階模態(tài)的振型向量；φb，j為 b階次第j階模態(tài)的振型向量。

MAC的值處在0和1之間。當(dāng)MAC=0時表示這兩個模態(tài)之間完全無關(guān)；當(dāng)MAC=1表示這兩個模態(tài)之間完全相關(guān)。MAF表示某階次前n階的頻率穩(wěn)定性，∑MAC表示該階次前n階振型模態(tài)保證準則的和?！芃AC越高則該階次的振型穩(wěn)定性越高。

1.3自動定階方法具體步驟

系統(tǒng)階次自動判別的具體步驟如下：

（1） 假定系統(tǒng)的階次依次取1，2，3，4，…，nmax，nmax的取值至少要大于系統(tǒng)的真實階次。求得每個階次對應(yīng)的系統(tǒng)狀態(tài)空間模型，此時系統(tǒng)矩陣的階數(shù)為2，4，6，8，…，2nmax，將所求模態(tài)參數(shù)依次寫入下三角矩陣。

（2） 計算相鄰階次（頻率矩陣相鄰行空間）之間的歐式距離di，從而求得每一階模態(tài)參數(shù)的模態(tài)置信因子（MAF），選取模態(tài)置信因子（MAF）最大的m個MAFi，i-1，MAFj，j-1，MAFk，k-1，…對應(yīng)的階次作為系統(tǒng)近似階次（第一次篩選）。

（3） 計算第二步中的m個系統(tǒng)近似階次的模態(tài)保證準則MACi，MACj， MACk，…。

（4） 最后計算MAF·∑MAC，該值最大時對應(yīng)的系統(tǒng)階次就是系統(tǒng)的最佳階次。求得系統(tǒng)最佳階次后，頻率下三角矩陣和阻尼比下三角矩陣對應(yīng)行的頻率值和阻尼比值以及對應(yīng)的振型就是系統(tǒng)的真實模態(tài)參數(shù)（第二次篩選）。

通過MAF及MAF·∑MAC對系統(tǒng)階次進行兩次篩選之后，保證了系統(tǒng)最佳階次的唯一性。本文將通過工程實驗對系統(tǒng)最佳階次的正確性進行驗證。

2算例〖2/3〗2.1三跨連續(xù)梁橋模型模態(tài)參數(shù)識別該三跨連續(xù)梁橋模型是在實驗室制作而成，圖2是實物圖。該模型尺寸為3800 mm×200 mm×30 mm，橋總長為3800 mm，橋面寬為200 mm，橋面厚為30 mm。三跨的跨度分別為1000，1700，1000 mm。圖2三跨連續(xù)梁橋?qū)嵨飯D

Fig.2A threespan continuous beam bridge

本次試驗數(shù)據(jù)采集系統(tǒng)使用北京東方振動與噪聲技術(shù)研究所的整套采集儀器，加速度器的布置如圖3所示。試驗用加速度傳感器可用頻率范圍為0.5 Hz～8 kHz，采樣頻率為512 Hz，采樣時長為30 s。用小鐵錘對模型沿著A，B，C，D，E共5個點依次敲擊。圖3加速度計布置位置（單位：mm）

Fig.3The position of acceleration sensorJ（Unit：mm）

2.1.1與傳統(tǒng)SSI方法識別結(jié)果比較

應(yīng)用本文提出的自動識別方法與傳統(tǒng)方法的識別結(jié)果比較，這三種方法所求得的前4階頻率如表1所示。

表 1隨機子空間法三種方法識別結(jié)果

Tab.1The results of three SSI methods

階次SSI自動識別方法SSI穩(wěn)定圖法SSI奇異值跳躍法126.36426.36326.362257.54157.76957.770364.87964.83264.832490.54490.41390.437

由表1知，SSI模態(tài)參數(shù)自動識別方法的識別結(jié)果與傳統(tǒng)的SSI方法識別的結(jié)果非常接近，它們之間的偏差不到1%。

2.1.2與其他模態(tài)參數(shù)識別方法識別結(jié)果比較

將SSI模態(tài)參數(shù)自動識別方法、頻域分解法（FDD）、希爾伯特黃變換（HHT）、ANSYS有限元模態(tài)分析的結(jié)果匯于表2。表2各方法頻率識別結(jié)果

Tab.2The results of various related methods

階次ANSYSSSI自動識別FDDHHT頻率偏差/%頻率偏差/%頻率偏差/%125.00026.3635.45226.3805.52026.4115.644254.68457.7795.66058.1306.30158.2416.505364.40464.8570.70365.0000.92064.3210.129492.82390.4812.52390.7502.23392.9810.170

由表2可知，SSI模態(tài)參數(shù)自動識別方法、FDD和HHT這三種方法識別的振動頻率非常接近，第1和2階頻率與ANSYS有限元結(jié)果相差5.5%左右，第3和4階頻率與ANSYS有限元結(jié)果最大僅相差2.5%。SSI模態(tài)參數(shù)自動識別方法識別結(jié)果與FDD和HHT識別結(jié)果僅相差1%以內(nèi)。

2.2濱州黃河大橋縮尺模型

濱州黃河大橋?qū)嶒炇铱s尺模型根據(jù)剛度相似準則確定構(gòu)件尺寸，選用1/40 的縮尺比例?？s尺模型橋面長15.2 m，橋面寬0.82 m，中塔高3.1 m，邊塔高1.9 m，橋面板最大板厚20 mm，斜腹板厚6 mm。主梁材料采用鋁合金，拉索采用拔絲鋼筋。主梁、橋塔和斜拉索的EA和EI剛度相似比誤差都控制在5%以內(nèi)。

用ANSYS對濱州黃河大橋的縮尺模型進行了建模，主梁采用BEAM4單元，彈性模量為6.86×104 N/mm2，泊松比取0.34，密度為3.2×10-3 kg/mm3；主塔、邊塔和橫隔梁采用BEAM188單元，彈性模量為6.86×107 N/mm2，泊松比取0.34，密度為2.7×10-3 kg/mm3；斜拉索采用LINK8單元，彈性模量為20×104 N/mm2，泊松比取0.3，密度為7.85×10-3 kg/mm3。ANSYS有限元模型如圖4所示。圖4縮尺模型的ANSYS有限元模型

Fig.4The ANSYS finite element model2.2.1振動試驗過程描述

本次試驗的加速度傳感器可用頻率范圍為0.5 Hz～8 kHz，采樣頻率為256 Hz，采樣時長為50 s。由于該橋沿主塔嚴格對稱，其振型也是沿主塔的主軸線對稱或反對稱。所以將加速度計布置在橋的左半側(cè)或右半側(cè)求得半橋的振型即可。本實驗將加速度計布置在橋的左半側(cè)，第1個加速度計布置在橋左側(cè)邊緣0.25 m處，加速度計之間間隔0.45 m，第16個加速度計布置在中塔左側(cè)0.5 m處，如圖5所示。圖5加速度計布置圖（單位：m）

Fig.5The position of acceleration sensor（Unit：m）

通過人在橋上走動作為橋的外部激勵方式。本次試驗共采集到215組數(shù)據(jù)。

2.2.2模態(tài)參數(shù)識別

將采集到的數(shù)據(jù)用本文所提出的SSI模態(tài)參數(shù)自動識別法進行識別，求得的系統(tǒng)階次為20，如圖6所示。圖 6系統(tǒng)的最佳階次

Fig.6The optimal order of the system

用SSI模態(tài)參數(shù)自動識別法對每一組數(shù)據(jù)進行識別，表3是一組數(shù)據(jù)所求系統(tǒng)的前6階模態(tài)與ANSYS有限元結(jié)果的對比。

SSI模態(tài)參數(shù)自動識別法識別了該橋的左半橋的振型，與ANSYS識別的前6階振型對比如圖7所示。圖7前6階振型

Fig.7The first six mode shapes

表3SSI自動識別法與ANSYS有限元頻率識別結(jié)果

Tab.3The results of SSI automatic identification method and ANSYS finite element for frequency

階次SSI自動

識別法/HzANSYS

理論值/Hz振型描述頻率值

偏差/%14.0714.110一階反對稱0.97528.6299.540一階對稱10.557310.68410.995二階反對稱2.911411.89612.822二階對稱7.784515.26915.494三階反對稱1.452616.32217.286三階對稱5.906

由表3可知，SSI模態(tài)參數(shù)自動識別方法識別的濱州黃河大橋縮尺模型前6階頻率與ANSYS模態(tài)分析結(jié)果相比較，相差均在10%以內(nèi)。由圖7可知，SSI模態(tài)參數(shù)自動識別方法識別的濱州黃河大橋縮尺模型前6階振型與ANSYS模態(tài)分析結(jié)果在走勢和形狀上比較一致。可見本文提出的SSI模態(tài)參數(shù)自動識別方法可以很好地識別出該模型的模態(tài)參數(shù)。

3結(jié)論

本文提出了基于隨機子空間法的模態(tài)參數(shù)識別方法SSI模態(tài)參數(shù)自動識別法， 將SSI模態(tài)參數(shù)自動識別法識別結(jié)果與ANSYS模態(tài)分析的結(jié)果以及其他模態(tài)參數(shù)的識別結(jié)果進行了對比：

（1） 本文提出了一種系統(tǒng)階次自動判別方法，并對該方法從理論上進行了推導(dǎo)，給出了自動定階方法計算流程，論證了系統(tǒng)最佳階次識別結(jié)果的唯一性，通過實驗驗證了識別結(jié)果的正確性。

（2）相比較于奇異值跳躍法和穩(wěn)定圖法，自動判別方法克服了它們的缺點，由于不需要人工參與，該方法實現(xiàn)了系統(tǒng)階次的自動識別，其識別結(jié)果是不受人的主觀判斷的影響。

（3） 提出了SSI模態(tài)參數(shù)自動識別方法，對3個不同類型結(jié)構(gòu)的模態(tài)參數(shù)進行識別，該方法識別出的振動頻率與ANSYS有限元結(jié)果、傳統(tǒng)SSI法、FDD和HHT識別結(jié)果很接近，表明SSI模態(tài)參數(shù)自動識別法可以很好地識別出這3種結(jié)構(gòu)的振動頻率和振型，且具有較高的魯棒性。因此SSI模態(tài)參數(shù)自動識別法可以作為一種模態(tài)參數(shù)自動識別方法運用于結(jié)構(gòu)實時在線參數(shù)識別。

總之，通過這兩個算例的計算，本文提出的SSI模態(tài)參數(shù)自動識別法能很好地識別出結(jié)構(gòu)的頻率和振型，表現(xiàn)出了良好地適用性和魯棒性。

參考文獻：

[1]常軍， 張啟偉， 孫利民. 隨機子空間產(chǎn)生虛假模態(tài)及模態(tài)遺漏的原因分析[J].工程力學(xué)， 2007， 24（11）： 57—62.

CHANG Jun， ZHANG Qiwei， SUN Limin.Analysis how stochastic subspace identification brings false modes and mode absence[J]. Engineering Mechanics， 2007， 24（11）： 57—62.

[2]Lennart L . System Identification： Theory for the User[M]. 北京： 清華大學(xué)出版社， 2002： 219—221.

[3]Herman V A， Bart P. International research projects on structural health monitoring： an overview[J]. Structural Health Monitoring， 2003， 2（4）： 341—358.

[4]Bart P. System identification and damage detection in civil engineering[D]. Belgium： Department of Civil Engineering， K.U. Leuven， 2000： 123—127.

[5]Bart P， Roeck D. Referencebased stochastic subspace identification for outputonly modal analysis [J]. Mechanical Systems and Signal Processing， 1999， 13（6）： 855—878.

[6]單德山， 徐敏. 數(shù)據(jù)驅(qū)動隨機子空間算法的橋梁運營模態(tài)分析[J]. 橋梁建設(shè)， 2011， 6： 16—21.

SHAN Deshan， XU Min. Operational modal analysis of bridge structure based on datadriven stochastic subspace identification algorithm[J].Bridge Construction， 2011， 6： 16—21.

[7]Bart P， Roeck D. Stochastic system identification for operational modal analysis： a review[J]. Journal of Dynamic Systems， Measurement and Control， 2001， 123（4）： 659—667.

[8]Filipe M， Cunha A， Caetano E. Online automatic identification of the modal parameters of a long span arch bridge[J]. Mechanical Systems and Signal Processing， 2009， 23 （5）： 316—329.

[9]Andersen P， Brincker R， Goursat M， et al. Automated modal parameter estimation for operational modal analysis of large systems[C]. Proceedings of the 2nd International Operational Modal Analysis Conference， Copenhagen： 2007： 299—308.

[10]肖祥， 任偉新. 實時工作模態(tài)參數(shù)數(shù)據(jù)驅(qū)動隨機子空間識別[J]. 振動與沖擊， 2009， 28（8）： 148—153.

XIAO Xiang，REN Weixin.Improved datadriven stochastic subspace identification of online operational modal parameters[J]. Journal of Vibration and Shock， 2009， 28（8）： 148—153.

A stochastic subspace identification method for automotive identification

of structural modal parameters

ZHANG Xiaoning1， DUAN Zhongdong2

（1.The Ningxia Hui Autonomous Region Electric Power Design Institute，Yinchuan 750000， China；

2. Harbin Institute of Technology（Shenzhen）， Shenzhen 518055， China）

Abstract： Realtime online monitoring of the structure put forward the demand for the automatic identification of the structural modal parameters. The current methods of the structural modal parameter identification require manual intervention， in order to realize realtime monitoring of the structure， this article attempts to develop an automatic structural modal parameters identification method. Stochastic subspace method only needs to determine one parameter， which is "system order". This paper proposes a method for automatic identification system order based on the frequency stability and vibration mode stability， based on those， we established an automatic identification method for structural modal parameters. Through identification two bridge examples， the reliability and applicability of the method was verified.Key words： modal parameters identification； stochastic subspace method； cablestayed bridge； health monitoring； system order作者簡介：張小寧（1988—），男，碩士研究生。電話：18295390955。Email： zxn880625@163.com