橋梁時變系統的張量子空間識別方法研究及試驗驗證

張二華，吳 滌，劉 昊，單德山

(1.四川省公路規劃勘察設計研究院有限公司,成都 610031；2.西南交通大學 土木工程學院，成都 610031)

近年來，隨著結構健康診斷技術的蓬勃發展并廣泛應用于大跨度橋梁，基于健康監測實測信號的橋梁系統識別研究方興未艾[1]。基于實測信號的系統識別方法易于從橋梁各種激勵環境下高精度獲得橋梁的模態信息，已成為橋梁狀態評估的關鍵[2]。其中子空間系統識別方法因人為干預少，識別結果魯棒性好，在橋梁結構系統識別中日益受到重視。目前應用于橋梁系統識別的子空間方法主要有：隨機子空間方法、遞歸子空間方法、確定子空間方法。國內外學者對3類方法均開展了廣泛深入的研究。

Zarbaf等[3]用隨機子空間方法與層次聚類算法，對某斜拉橋斜拉索力進行了估計；徐健[4]將隨機子空間方法與經驗模態分解方法相結合，對某大跨斜拉橋進行了系統識別；Bui等[5]分別采用能量譜密度、協方差驅動隨機子空間及數據驅動隨機子空間的方法，并結合峰值拾取法和穩定圖，對某混凝土簡支梁進行了系統研究。Khan等[6]采用遞歸隨機子空間方法，研究了某大跨斜拉橋溫度對模態參數影響規律；Huang等[7]討論了遞歸隨機子空間方法中遺忘因子、系統矩陣維度大小等參數對地震激勵下結構模態參數識別的影響。秦世強等[8]將確定性隨機子空間方法用于菜園壩長江大橋的模態參數識別；陳棟軍[9]采用改進后的確定-隨機子空間方法識別了某混凝土斜拉橋的模態參數。Huang等[10]采用確定子空間方法識別了地震激勵下多輸出系統的模態參數。

綜上所述，隨機子空間方法常用結構一次完整動力響應來識別模態參數，而忽略非線性結構模態參數的時變特性。遞歸隨機子空間算法，能反映結構模態參數的時變規律，但識別結果精度較低，難以用于大跨復雜橋梁的時變模態參數識別。確定性隨機子空間識別精度較高，但對實際橋梁激勵和響應數據的完備性要求較高[11]。并且，子空間系統識別方法的現有改進多集中于算法的運算效率、數據前處理技術、模態自動識別方面，尚未看到有從子空間方法的基礎理論——矩陣計算方面加以改進和突破的研究文獻。

針對矩陣運算的子空間系統識別理論與方法的缺陷，本文基于短時滑窗子空間系統識別框架，在原有2維Hankel矩陣基礎上，引入時間尺度，建立3維時變Hankel張量；結合張量平行因子分解理論，提出基于張量的平行系統矩陣估計方法；結合穩定圖方法從模態參數識別結果中篩選真實模態信息，提出一種新的張量子空間系統識別理論與方法。基于曲線斜拉模型橋、斜拉橋振動臺試驗的實測數據，以計算耗時、穩定圖穩定軸及虛假模態為評判指標，通過與已有短時滑窗子空間系統識別方法對比，驗證了本文所提方法較現有方法的優越性與工程實用性。

1 張量子空間系統識別理論

1.1 張量子空間系統識別數學模型

1.1.1 滑窗子空間系統識別框架

根據空間狀態控制理論[12]，連續確定性系統的狀態空間模型可寫為如式(1)和式(2)所示

(1)

y(t)=Ccx(t)+Dcu(t)

(2)

實際運用中，系統采樣數據均為離散形式，將連續狀態空間模型轉變為離散狀態空間模型，式(1)～式(2)改寫為

xk+1=Axk+Buk+wk

(3)

yk=Cxk+Duk+vk

(4)

式中:xk=x(kΔt);A,B,C,D為連續系統矩陣Ac,Bc,Cc,Dc的離散形式;uk,yk為離散輸入和輸出矩陣;wk為過程噪聲;vk為測量噪聲。

當系統處于環境激勵下，其激勵源往往不可測，通常將環境激勵假定為與隨機噪聲影響相似，則含有uk項與噪聲項可合并,隨機子空間系統狀態空間模型可寫為

xk+1=Axk+wk

(5)

yk=Cxk+vk

(6)

無論是確定子空間識別或隨機子空間系統識別問題，均可通過將輸入-輸出數據或輸出數據組成特定的Hankel矩陣，計算該矩陣的行空間投影，對投影進行奇異值分解(singular value decomposition,SVD)，從而得到可觀測矩陣和狀態序列的卡爾曼濾波估計，進而確定系統矩陣以及噪聲協方差矩陣。然后經特征值分解，獲得系統的模態參數，其求解過程可參考Van Overschee等的研究。

因確定子空間與隨機子空間系統識別方法均基于線性系統假設，近年來，Moaven等[13]將滑窗技術引入到上述方法中，以解決線性子空間系統識別方法在追蹤復雜非線性時變系統模態參數方面的不足。其核心思想為：基于短時平穩性假設，對測試數據進行分段，各窗口內將結構視為時不變系統；然后對每個窗口內動力數據重復使用子空間識別方法求解結構模態參數，達到對非線性系統時變模態參數的追蹤。算法流程圖,如圖1所示。

圖1 滑窗子空間系統識別算法流程圖Fig.1 The sliding window based subspace system identification flow chart

由圖1可知：上述方法未從根本上對基于矩陣運算的子空間系統識別理論和方法進行改進，而是不同時段內對方法的重復使用，因此其耗時為N次子空間求解過程花費時間的總和。

1.1.2 張量子空間數學模型

為保持和張量表達形式一致，滑窗子空間中t時刻的輸入輸出Hankel矩陣重新寫為U2i,j,t、Y2i,j,t，其中第一下標2i表示t時刻Hankel矩陣的行數，第二下標j表示Hankel矩陣的列數，第三下標t表示當前的時刻。假設t+1時刻新的數據列為ut+1和yt+1，形成如式(7)的行向量空間

(7)

將式(7)行向量空間附于t時刻Hankel矩陣的最右一列，同時將原U2i,j,t、Y2i,j,t最左一列剔除，保持t+1時刻的Hankel矩陣與t時刻Hankel矩陣維度相同，形成t+1時刻的輸入輸出Hankel矩陣，如式(8)所示

(8)

假設信號時間窗口為L，可將信號時程為T的時變振動信號劃分為N個時間窗，則t=1,2,…,N，按信號采集時間先后，構建時變3維Hankel張量，3維張量沿時間維度的切片為

(9)

基于上述時變3維Hankel張量，對系統矩陣進行估計，首先需要解決的問題是剖析其數學本質，即建立對應的數學模型進行求解。為此，以下將首先通過理論推導分析非線性系統的系統矩陣的時變性，在此基礎上，建立系統矩陣估計的數學模型。

(10)

式中:矩陣U∈RM×M,V∈RN×N為正交矩陣;UsU(1:M,1:r);VsV(1:N,1:r);Σs=diag(σ1,…,σr),σ1≥σ2≥…≥σr≥σr+1=…σM=0。

由式(10)經推導，可得到

(11)

在X中加入白噪聲項Ψ，則可得到受噪聲污染的實測信號為

Y=X+Ψ

(12)

式中，Ψ為測量噪聲分量。

當數據長度足夠大時，存在如式(13)關系

(13)

進一步可得到式(14)

(14)

(15)

式中,S=diag(s1,…,sM),si>0。si表達式為

(16)

由式(15)和式(16)可知，含噪數據組成矩陣Y經SVD后左奇異矩陣U和對角陣S包含結構時變動力指紋信息，對于時變的非線性系統，左奇異矩陣U和對角陣S集中反映了系統時變的特性

由張量展開和高階張量SVD，時變Hankel張量沿時間維度展開的切片可類似看作某時段如式(12)含噪信號Y，則由式(9)～式(15)可得，張量子空間系統矩陣求解的數學模型為

H(k)=UkSkVT+Rk

(17)

式中:Rk為擬合的殘差平方和;Uk為酋矩陣;Sk為非負對角陣;k為沿時間軸展開的3維Hankel張量H“切片”下角標。

式(17)在形式上與張量平行因子分解理論中的PARAFAC模型改進形式PARAFAC2模型保持一致。因此，以下將引入改進的張量平行因子分解理論，求解上述數學問題。

1.2 基于張量分解的系統矩陣估計

1.2.1 PARAFAC2模型簡介

Harshman[14]在1972年提出了PARAFAC2模型，其數學表達為：定義張量的展開式Xk,k=1,…,K為矩陣集合，其中，Xk是Ik×J的矩陣，Ik的值隨著k的改變而改變，假設PARAFAC2分解的維度為R，此時，PARAFAC2具有

Xk≈UkSkVT,k=1,…,K

(18)

式中:Uk為Ik×R的矩陣;Sk為R×R的對角矩陣;k=1,…,K;V為J×R的因子矩陣。

PARAFAC2分解示意圖,如圖2所示。

圖2 PARAFAC2分解示意圖Fig.2 Illustration of PARAFAC2

Xk≈QkHSkVT,k=1,…,K

(19)

(20)

1.2.2 系統矩陣估計

基于約束的PARAFAC2模型，式(17)可視為廣義的主成分分析問題，其問題可表達為

σk(Q1,…,Qk,H,V,S1,…,Sk)=

(21)

首先求解Qk，式(21)的最小化問題，等價于式(22)最大化問題

(22)

(23)

則列正交的Qk為

(24)

通過H,S1,…,Sk,V，最小化式(21)，轉變為下述問題

σk(H,V,S1,…,Sk|Q1,…,Qk)=

(25)

綜上所述，利用PARAFAC2模型平行求解張量系統矩陣過程可歸納為：

步驟4由式(25)判斷求解過程是否達到收斂準則，若σold,k-σnew,k>εσold,k，ε為收斂準則(經試算，本文中ε取值為0.1)，則回到步驟3，否則，進入下一步。

至此，基于PARAFAC2模型，包含結構時變系統參數的系統矩陣Uk和Sk已求解得到。

則每個時間窗內的擴展觀測矩陣可按式(26)進行求解

(26)

式中:Wk為權重影響因子;Γk為擴展可觀測矩陣。由擴展觀測矩陣的移位不變性，可計算獲得系統矩陣A,C，經特征值分解運算，即可獲得系統的頻率、阻尼比和振型，具體計算過程可參考Van Overschee等的研究，由于篇幅所限，不在贅述。

1.3 求解過程對比分析

本文所提出的張量子空間系統識別方法與兩種滑窗子空間系統(確定系統或隨機系統)識別方法的求解過程流程圖,如圖3～圖4所示。

圖3 滑窗子空間系統識別方法流程示意圖Fig.3 Schematic diagram of the sliding window-based subspace system identification method

圖4 張量子空間系統識別方法求解過程示意圖Fig.4 Schematic diagram of the tensor subspace system identification method

由圖3～圖4可知：假設測試信號劃分時間窗口數為K，穩定圖中最大系統階次為N，對于滑窗子空間系統識別方法，每一窗口內數據均需進行一次完整的基于矩陣運算的子空間系統識別過程，則需進行的SVD次數總共為K·N/2；而本文中張量子空間系統識別方法，通過K個時間窗口構建時變3維Hankel張量，在系統識別過程中，僅在穩定圖分析過程中，每一階次進行兩次SVD，系統識別全過程SVD次數總共為2·N/2=N次，隨著K增大(K為大于2的整數)，SVD次數將成倍減少。而已有研究表明[16]，子空間系統識別過程中SVD是耗時嚴重的步驟之一，且SVD過程極易受到噪聲影響導致奇異值發生變化，致使最后模態參數識別結果不準確，在穩定圖上表現為穩定點散亂，穩定軸扭曲。

綜上所述，本文所提出的張量子空間系統識別方法較滑窗子空間系統識別方法過程，計算效率更高，計算結果受噪聲影響更少，穩定軸將更加清晰可辨，更易識別到系統更多階次的模態參數信息。以下通過兩座斜拉模型橋進一步驗證上述理論分析的結論。

2 試驗驗證

2.1 曲線斜拉模型橋試驗

采用一座幾何縮尺比例為1∶20的曲線斜拉模型橋，驗證本文所提方法用于橋梁隨機子空間系統識別問題的可行性。

模型斜拉橋由5跨構成，跨度布置為(2.45+4.05+14.24+4.05+2.45)m，橋梁總長27.25 m，以中跨跨中為界分為曲線梁段和直線梁段，曲線梁段曲率半徑為27.5 m。主梁為Π型鋼梁，寬1.1 m，高0.081 m，由板厚6 mm的Q235鋼板焊接制作。斜拉索采用雙幅索面布置，每幅60對，由10 mm鋼絲繩制作。輔助墩與過渡墩由C30混凝土與Φ6 mm光面鋼筋制作完成。橋塔為鉆石型橋塔，由板厚10 mm的Q235鋼板焊接為矩形空心截面，橋塔每側鋼板中軸線位置由6 mm的Q235鋼板焊接成加緊肋。橋梁模型示意圖,如圖5所示。

圖5 曲線斜拉模型橋示意圖(m)Fig.5 A model cable-stayed bridge (m)

為獲得主梁損傷后的動力測試信號，在主梁截面中線布置941B型加速度傳感器，傳感器布置方式為：分別在兩側邊跨跨中各布置1個測點，中跨八分點各布置1個測點，共計11個測點。信號采樣頻率為256 Hz，在中跨跨中位置處施加脈沖激勵后，采集主梁自由振動信號，采集時長為40 s，如圖6所示。

穩定圖是子空間方法中用于篩選結構真實模態信息的常用方法，本文選取10 s長度信號作為初始數據，滑窗步長為0.4 s，由滑窗隨機子空間識別方法與張量子空間系統識別方法分別獲得系統穩定圖,如圖7所示，前3階識別結果的對比,如圖8所示。圖7中，“◇”穩定極軸代表結構頻率、阻尼比與振型結果均穩定，“*”極軸代表僅頻率穩定，“○”極軸代表僅頻率與阻尼比穩定。穩定圖方法以頻率、阻尼比與振型結果均穩定的極軸代表結構的真實模態結果。穩定圖穩定極軸條件及判別方法可參考Khan等的方法。本節中曲線斜拉模型橋的理論模態參數識別結果可參考文獻[17]。

(a)頻率

由圖7～圖8可知：① 由張量子空間方法所得穩定圖紅色穩定軸的數量較滑窗隨機子空間方法更多，如2.12 Hz處，滑窗隨機子空間方法并未得到“◇”穩定軸；張量子空間方法所得虛假模態散點數量更少，如在3～6 Hz內，滑窗子空間方法所得結果出現多個“*”零星散點，張量子空間方法結果零星散點較少，魯棒性更好；② 張量子空間方法識別到的頻率與振型結果與傳統子空間方法結果均極為接近，兩種方法識別出的阻尼比值差異較大，這主要是由于基于子空間的時域類系統識別方法在識別阻尼比過程中對噪聲極為敏感造成的。③ 滑窗隨機子空間方法得到上述結果用時337.49 s，而張量子空間識別方法僅用時280.31 s(計算機配置CPU均為I7-6700K，內存均為16 G)。以上結果表明，針對隨機系統識別問題，本文所建立的張量子空間方法較滑窗隨機子空間方法在大幅提高計算效率的同時，能夠從信號中挖掘更多系統模態信息，且識別精度與滑窗隨機子空間方法相同。

2.2 斜拉模型橋振動臺試驗

采用一座縮尺比例為1∶20的振動臺斜拉模型橋，驗證本文所提方法用于確定子空間系統識別問題的可行性。

斜拉模型橋由5跨構成，跨度布置為(2.9+3.6+19.0+3.6+2.9)m，橋梁總長32 m。斜拉模型橋索塔、過渡墩及輔助墩均選用M15微粒混凝土模擬原橋混凝土材料，鋼筋采用φ6圓鋼作為縱筋、10#鋼絲作為箍筋。模型橋主梁設計為矩形空心斷面，采用厚度為5 mm鋼板制作。橋塔為門式框架結構，高度為4.55 m；邊墩、輔助墩高度為1.9 m。拉索采用直徑10 mm鋼絲繩模擬。斜拉模型橋示意圖,如圖9所示。

圖9 振動臺斜拉模型橋示意圖(cm)Fig.9 A model cable-stayed bridge tested on a shaking table (cm)

主梁與橋塔分別布置941B型加速度傳感器，測點布置方式為：沿橋塔兩側塔柱間隔1 m布置縱橋向和橫橋向加速度計，每側塔柱布置5個測點；中跨主梁在1/8測點各布置1個豎向和橫向加速度計，兩邊跨跨中及各支點處均布置1個豎向和橫向加速度計；各振動臺分別布置1個縱橫向加速度計；輔助墩墩頂各布置1個縱向和橫向加速度計。信號采樣頻率256 Hz，采樣時長為47 s，橋塔地震加速度響應,如圖10所示。

本節同樣采用穩定圖作為結構真實模態的篩選方法，滑窗步長為5 s，由滑窗確定子空間識別方法與張量子空間系統識別方法分別獲得系統穩定圖，如圖11所示，前3階識別結果的對比如圖12所示。圖11中，“◇”極軸、“*”極軸、“○”極軸含義與2.1節相同，穩定圖方法以頻率、阻尼比與振型結果均穩定極軸代表結構的真實模態結果。穩定圖穩定極軸條件及判別方法可參考Huang等的研究。本節中斜拉模型橋的理論模態參數識別結果可參考文獻[18]。

(a)頻率

由圖11～圖12可知：① 針對確定性系統而言，張量子空間方法所得穩定圖中“◇”穩定極軸上的極點數量較滑窗確定子空間方法更多，如0～10 Hz內，滑窗確定子空間方法雖然也有“◇”穩定軸，但穩定極點數量較張量子空間方法少，這表明張量子空間方法求解過程魯棒性更好，結果更加穩定；② 張量子空間方法識別到的頻率與振型結果與傳統確定子空間方法結果均極為接近，兩種方法識別出的阻尼比值差異較大，其原因與2.1節相同，不再贅述；③ 滑窗確定子空間方法得到全部結果用時179.28 s，而張量子空間方法僅用時61.79 s，耗時減少約2/3(計算機配置CPU均為I7-6700K，內存均為16 G)。以上結果表明，針對確定性系統識別問題，本文所提出的張量子空間方法較滑窗確定子空間方法能大幅提高計算效率，并能從信號中挖掘更多系統模態信息，且識別精度相同。

3 結 論

在傳統2維矩陣運算的子空間系統識別理論基礎上，通過引入時間維度將傳統3維數據陣列拓展到3維，基于張量展開及平行因子分解理論，結合穩定圖方法，建立了適用于時變非線性系統的張量子空間快速系統識別理論及方法，所得主要結論如下：

(1)基于本文理論推導，在2維Hankel矩陣中引入時間維度，提出時變3維Hankel張量的概念，由張量展開和高階張量SVD運算，建立沿時間維度展開的張量子空間系統識別數學模型，其在形式上與PARAFAC2模型一致。

(2)在3維Hankel張量基礎上，基于約束的PARAFAC2模型，可快速求解時變的系統矩陣，通過與本文中傳統基于矩陣運算的滑窗子空間系統識別方法對比，驗證了方法在計算效率與結果魯棒性上的提升。

(3)通過曲線斜拉模型橋試驗和斜拉模型橋振動臺試驗驗證表明，本文所提出的張量子空間系統識別方法適用于求解確定和隨機非線性橋梁系統的瞬時模態參數追蹤識別問題。