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

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

2、1+.3 文獻(xiàn)標(biāo)志碼： A 文章 編號： 10044523 （2017）04054207DOI：10.16385/ki.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)進(jìn) 行實時在線監(jiān)測，而這也對結(jié)構(gòu)模態(tài)參數(shù)自動識別提出了要 求。隨機子空間法（ SS）I 是近年來發(fā)展起

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

4、圖法。奇 異值跳躍法是根據(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 到

5、nmax 每個階次的模態(tài)參數(shù)，用求得的模態(tài)參數(shù)構(gòu)造下三角矩陣， 然后判定頻率的穩(wěn)定性一一模態(tài)置信因子(MAF )和振型的穩(wěn)定性一一模態(tài)保證準(zhǔn)則( 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+wkyk=Cxk+vk (1)式中A為系統(tǒng)的狀態(tài)

6、矩陣；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 矩陣進(jìn)行 QR 分解，得到投影矩陣；對投影矩陣進(jìn)行奇異值分解(SVD)，得到擴展的可觀測矩陣和系統(tǒng)狀態(tài)的卡 爾曼濾波；最后由卡爾曼濾波序列和系統(tǒng)輸出采用最小二乘 法求得系統(tǒng)的狀態(tài)矩陣 A和輸出矩陣C。對系統(tǒng)狀態(tài)矩陣A 進(jìn)行特征值分解 A= A-1 (2)式中A為包含復(fù)特征值的 對角矩陣；為特征向量組成的矩陣。上式中，A =diag入i,其特征值兩

7、兩共軛。 令入i=ai jbi(3) 可算得系統(tǒng)的振動頻率 fi和阻尼比 i: fi = a2i + b2i /(2 n ) (4)e i = -ai a2i + b2i (5)將所求模態(tài)參數(shù)構(gòu)造下三角矩陣，故最終求得了 3個nmax階的下三角方陣 Fnmax=f1, 1000 f2 , 1f2 , 20 0f3 , 1f3 , 2f3, 3 0fnmax, Ifnmax , 2fnmax , 3fnmax, nmax (6) 式中Fnmax為振動頻率的下三角矩陣。fa, b為第a階次的第b 階振動頻率。同時也可以得到模態(tài)保證準(zhǔn)則(MAC )和阻尼比的下三角矩陣。1.2.2 計算模態(tài)置信因子(

8、 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） 212 （7）式中 d （x1, x2,， xn）為n維空間的歐氏距離；xn為所求階次的第n階模態(tài)參 數(shù)（頻率、阻尼比、 MAC）。在這里，定義一種專門針對于振動頻率穩(wěn)定性的因子模態(tài)置信因子（MAF），模態(tài)置信因子在數(shù)值上等于歐氏距離的倒數(shù)（d （x1, x2,，

9、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)階次的升高而減小,差異達(dá)到最小后又會 隨著系統(tǒng)階次的升高而增大,歐氏距離表示了相鄰階次模態(tài) 參數(shù)之間的差異,歐氏距離和系統(tǒng)階次之間的關(guān)系如圖 1 所 示。所以, MAFa, b 越小（歐氏距離越大）表示該階模態(tài)的 置信度越低, MAFa, b 越大（歐氏距離越

10、?。┍硎驹撾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個階次的歐氏距離 di, 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)的

11、階次就是系統(tǒng)的近似階次，建 議 m 取至少大于 3。1.2.3 模態(tài)保證準(zhǔn)則( MAC)在上一步驟中， 選取了模態(tài)置信度最高的 m 個階次， 這 m 個階次可以看作是頻率穩(wěn)定的階次，但它們都是系統(tǒng)的近 似階次，可能含有噪聲模態(tài)的階次。為了避免噪聲模態(tài)，還 應(yīng)該對這 m 個頻率穩(wěn)定的階次進(jìn)行振型的判定， 看其振型是 否穩(wěn)定。 最后選取這 m 個階次中頻率穩(wěn)定性和振型穩(wěn)定性最 高的階次作為系統(tǒng)的最佳階次。振型穩(wěn)定性用模態(tài)保證準(zhǔn)則(Modal Assurance Criterion )來定義 MACa, b= $ Ta, b, j2 ( Ta, j$ a, j) ( $ Tb, j $ b, j)

12、(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)保證準(zhǔn)則的和。刀 MAC越高則該階次的振型穩(wěn) 定性越高。1.3 自動定階方法具體步驟 系統(tǒng)階次自動判別的具體步驟如下：（1）假定系統(tǒng)的階次依次取 1, 2, 3, 4,，nmax, nmax 的取值至少要大于系統(tǒng)的真實階次。求得每個階次對 應(yīng)的

13、系統(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)保證準(zhǔn) 則 MACi, MACj, MACk,。（4）最后計算MAF?刀MAC,該值最大時對應(yīng)的系統(tǒng) 階次就是系統(tǒng)的最佳階次。求得系統(tǒng)最佳階次后，頻率下三 角矩陣和阻尼比下三角矩陣

14、對應(yīng)行的頻率值和阻尼比值以 及對應(yīng)的振型就是系統(tǒng)的真實模態(tài)參數(shù)（第二次篩選） 。通過MAF及MAF?刀MAC對系統(tǒng)階次進(jìn)行兩次篩選之后，保證了系統(tǒng)最佳階次的唯一性。本文將通過工程實驗對系統(tǒng) 最佳階次的正確性進(jìn)行驗證。2算例 2/3 2.1 三跨連續(xù)梁橋模型模態(tài)參數(shù)識別該三 跨連續(xù)梁橋模型是在實驗室制作而成，圖 2 是實物圖。該模 型尺寸為 3800 mm x 200 mm x 30 mm，橋總長為 3800 mm , 橋面寬為200 mm，橋面厚為30 mm。三跨的跨度分別為 1000, 1700, 1000 mm。圖2三跨連續(xù)梁橋?qū)嵨飯DFig.2A threespan continuous

15、beam bridge 本次試驗數(shù)據(jù)采集系統(tǒng)使用北京東方振動與噪聲技術(shù) 研究所的整套采集儀器，加速度器的布置如圖 3所示。試驗 用加速度傳感器可用頻率范圍為0.5 Hz8 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é)

16、果Tab.1The results of three SSI methods階次SSI自動識別方法SSI穩(wěn)定圖法SSI奇異值跳躍法126.36426.36326.362257.54157.76957.770364.87964.83264.8 32490.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

17、 results of various related methods階次ANSYSSS自動識別FDDHHT頻率偏差/%頻率偏差/% 頻率偏差/%125.00026.3635.45226.3805.52026.4115.644254.68457.779 5.66058.1306.30158.2416.505364.40464.8570.70365.0000.920 64.3210.129492.82390.4812.52390.7502.23392.9810.170由表2可知，SSI模?B參數(shù)自動識別方法、FDD和HHT這三種方法識別的振動頻率非常接近，第 1 和 2 階頻率與 ANSYS 有

18、限元結(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ù)剛度相似準(zhǔn)則確定 構(gòu)件尺寸， 選用 1/40 的縮尺比例。 縮尺模型橋面長 15.2 m， 橋面寬0.82 m，中塔高3.1 m，邊塔高1.9 m，橋面板最大板 厚20 mm，斜腹板厚6 mm。主梁材料采用鋁合金，拉索采 用拔絲鋼筋。主梁、橋塔和斜拉索的EA和EI剛度相似比誤差都控制在 5%以內(nèi)用ANSYS對濱州黃河大橋的縮尺模型進(jìn)行了建模，主梁采用BEAM4單元，彈性模量為6.

19、86 X 104 N/mm2 ,泊松比取 0.34，密度為3.2 X 10-3 kg/mm3 ;主塔、邊塔和橫隔梁采用 BEAM188單元，彈性模量為 6.86X 107 N/mm2，泊松比取 0.34, 密度為2.7X 10-3 kg/mm3 ;斜拉索采用LINK8單元，彈性模 量為 20X 104 N/mm2 ，泊松比取 0.3，密度為 7.85X 10-3 kg/mm3。ANSYS有限元模型如圖4所示。圖4縮尺模型的 ANSYS有限元模型Fig.4The ANSYS finite element model2.2.1 振動試驗過程 描述本次試驗的加速度傳感器可用頻率范圍為0.5 Hz8

20、kHz,采樣頻率為256 Hz,采樣時長為50 s。由于該橋沿主塔嚴(yán)格 對稱,其振型也是沿主塔的主軸線對稱或反對稱。所以將加 速度計布置在橋的左半側(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ù)識別將采集到的

21、數(shù)據(jù)用本文所提出的SSI模態(tài)參數(shù)自動識別法進(jìn)行識別，求得的系統(tǒng)階次為20，如圖 6 所示。圖 6 系統(tǒng)的最佳階次Fig.6The optimal order of the system用SSI模態(tài)參數(shù)自動識別法對每一組數(shù)據(jù)進(jìn)行識別，表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 identif

22、ication methodand 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)

23、。由圖7可知，SSI模態(tài)參數(shù)自動識別方法識 別的濱州黃河大橋縮尺模型前6階振型與ANSYS模態(tài)分析結(jié)果在走勢和形狀上比較一致?？梢姳疚奶岢龅腟SI模態(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é)果進(jìn)行 了對比：（1）本文提出了一種系統(tǒng)階次自動判別方法，并對該 方法從理論上進(jìn)行了推導(dǎo)，給出了自動定階方法計算流程， 論證了系統(tǒng)最佳階次識別結(jié)果的唯一性，通過實驗驗證了識 別結(jié)果的正確性。（2）相比較于奇異值跳躍法和穩(wěn)定圖法，自動判

24、別方 法克服了它們的缺點，由于不需要人工參與，該方法實現(xiàn)了 系統(tǒng)階次的自動識別，其識別結(jié)果是不受人的主觀判斷的影 響。（3）提出了 SSI模態(tài)參數(shù)自動識別方法，對 3個不同 類型結(jié)構(gòu)的模態(tài)參數(shù)進(jìn)行識別，該方法識別出的振動頻率與 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)出了 良好地適用性和魯

25、棒性。參考文獻(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 absenceJ. Engineering Mechanics ，2007，24(11)：5762.2Lennart L . System Identification ： Theory for theUserM. 北京： 清華大學(xué)出版

26、社， 2002：219221.3 Herman V A， Bart P. International research projects onstructural health monitoring ： an overviewJ. Structural Health Monitoring ，2003，2( 4)：341 358.4 Bart P. System identification and damage detection incivil engineeringD. Belgium ： Department of Civil Engineering ， K.U. Leuven， 200

27、0：123127.5 Bart P ，Roeck D. Referencebased stochastic subspaceidentification for outputonly modal analysis J. MechanicalSystems and Signal Processing 1999,13 (6):855 878.6 ?蔚律劍 ? 徐敏 . 數(shù)據(jù)驅(qū)動隨機子空間算法的橋梁運 營模態(tài)分析 J. 橋梁建設(shè)， 2011， 6：1621.SHAN Deshan， XU Min. Operational modal analysis of bridge structure base

28、d on datadriven stochastic subspace identification algorithmJ.Bridge Construction ，2011，6：1621.7 Bart P ， Roeck D. Stochastic system identification for operational modal analysis ： a reviewJ. Journal of Dynamic Systems， Measurement and Control ， 2001 ， 123( 4)： 659 667.8 Filipe M ， Cunha A， Caetano

29、E. Online automaticidentification of the modal parameters of a long span arch bridgeJ. Mechanical Systems and Signal Processing，2009，23 (5)：316329.9 Andersen P， Brincker R， Goursat M， et al. Automated modal parameter estimation for operational modal analysis of large systemsC. Proceedings of the 2nd

30、 International Operational Modal Analysis Conference ， Copenhagen： 2007： 299 308 .10 肖祥，任偉新 . 實時工作模態(tài)參數(shù)數(shù)據(jù)驅(qū)動隨機子空間識別J.振動與沖擊， 2009,28 (8):148153.XIAO Xiang， REN Weixin.Improved datadriven stochastic subspace identification of online operational modal parametersJ. Journal of Vibration and Shock ， 2009， 28( 8)： 148153.A stochastic subspace identification method for automotive identificationof structural modal parametersZHANG Xiaoning1， DUAN Zhongdong2( 1.The Ningxia Hui Autonomous Region Electric Power Design Institute ， Yinchuan 750000， China；2. Harbin Institute of Technology ( Shenzhen)， Shenzh