基于隨機子空間方法的懸臂結構損傷識別研究

唐盛華， 方 志， 張國剛

(1. 湘潭大學 土木工程與力學學院，湖南 湘潭 411105； 2. 湖南大學 土木工程學院，長沙 410082；3. 湖南省交通規劃勘察設計院，長沙 410008)

損傷識別是結構健康監測系統的重要組成部分，損傷識別的目的在于通過現場測試數據，分析判斷結構的狀態，從而區別結構損傷與非損傷狀態，達到對結構進行損傷預警、損傷定位及損傷程度估計等健康監測的目的。如何有效地提取對結構損傷敏感的特征參數是損傷識別的關鍵，損傷一般會引起結構剛度的下降，導致結構模態參數發生變化，如頻率的降低。由于結構動力參數的測試不影響結構的正常使用，如橋梁結構模態參數測試時可以不封閉交通，因而，采用結構動力參數進行損傷識別非常普遍。健康監測系統能實時采集結構的各種動力響應信號，使基于動力參數的損傷識別方法也更加方便。

土木工程結構由于其質量和剛度一般很大，進行激勵非常困難，故利用環境激勵下的響應信號進行損傷識別是一個重要研究方向。隨機子空間方法基于結構受未知白噪聲激勵的假定，非常適合土木工程結構的動力特性描述。自Overschee等[1]提出隨機子空間方法(SSI)以來，國內外的學者在隨機子空間方法識別結構模態參數方面取得了較大進展[2-9]，SSI方法模態參數識別過程需輸入的參數量少，不需要迭代計算，無收斂性問題。亦有學者將隨機子空間方法應用于結構損傷識別中，取得了一定的研究成果，Basseville等[10]基于子空間方法，采用χ2檢驗識別結構損傷；Yan等[11-12]利用SSI方法獲得系統未損傷狀態下的卡爾曼濾波模型，由模型輸出殘差的大小來判斷結構有無損傷，此外，還利用子空間夾角進行損傷檢測；林友勤等[13-14]由系統矩陣能控標準形構造結構損傷指標主要用來區別結構的不同狀態；Kim等[15]使用觀測標準型對一框架結構進行損傷識別。對于實際工程結構，目前研究者多注重結構有無損傷的判斷，鮮見使用隨機子空間方法進行損傷定位的試驗研究。

本文基于隨機子空間方法，提出了一種損傷特征指標構造方法。對振動測量信號，利用數據驅動隨機子空間識別方法得到隨機狀態空間模型。通過兩個變換矩陣T、T1將狀態空間矩陣A轉化為特定形式，由向量馬氏距離構造損傷指標，通過矩陣T變換后的A矩陣對結構進行有無損傷及相對損傷程度的判斷，通過矩陣T1變換后的A矩陣進行損傷定位。通過一三自由度懸臂結構仿真算例和一座模型斜拉橋橋塔損傷試驗，表明該方法的有效性。

1 數據驅動隨機子空間模型

對m維的時間序列yk=y(kΔt) (k=1,2,…,N)，N為總采樣點數， Δt為采樣時間間隔，則隨機狀態空間模型為

(1)

式中：A∈Rn×n為離散狀態空間矩陣，C∈Rm×n為離散輸出矩陣(m為測點數、n為系統階次)，xk為n維狀態向量，wk∈Rn、vk∈Rm分別為建模誤差噪聲、測量噪聲。

輸出數據構造的Hankel矩陣為：

(2)

式中： Hankel矩陣包含2i個“塊行”，每一塊行由m行組成，j為矩陣的列數。

數據驅動隨機子空間識別方法采用投影法進行計算：

(3)

投影矩陣Pi可表示成可觀矩陣Qi與卡爾曼濾波狀態向量的乘積：

(4)

Oi=[CCA…CAi-1]T

(5)

將隨機狀態空間模型(1)通過Kalman狀態濾波，得到新息狀態模型

(6)

新息狀態模型(6)由三元組(A,C,K)完全確定，由線性系統的理論可知，具有相同輸入-輸出關系的系統，狀態方程與狀態變量的選取有關，對于式(6)所表達的系統，考慮滿秩線性變換

(7)

將式(7)代入式(6)得

(8)

A2=TAT-1,C2=CT-1,K2=TK

(9)

則稱(A,C,K)、(A2,C2,K2)是等價的。如果兩個系統等價，則A、A2相似，故它們有相同的特征值。由于T矩陣無窮多，因此，便產生了狀態模型多樣性的問題。

2 變換矩陣

設由兩個時刻未損傷狀態數據得到的系統矩陣為(A,C,K)、(A2,C2,K2)，對矩陣A、A2進行特征值分解

(10)

式中：Ψ,Ψ2∈Cn×n是特征向量矩陣，Λ,Λ2∈Cn×n是特征值矩陣。

因均為無損傷狀態下的結果，故兩個系統等價，即存在非奇異變換矩陣T，使得式(11)成立

(11)

因Λ=Λ2，于是

(12)

即

(13)

故通過變換矩陣T可以將其它狀態的A矩陣轉化為某一特定狀態。

為了進行損傷定位，構造變換矩陣[16-17]

(14)

可見，T1為式(5)中Oi矩陣的前兩個塊行組成的矩陣。

假設：

(15)

則有：

(16)

于是

(17)

(18)

即通過變換矩陣T1使得矩陣C化為單位矩陣和零矩陣的組合，變換后的A矩陣前m行由零矩陣和單位矩陣組成，后m行與測點相對應，嘗試用來檢測損傷位置。

3 損傷指標

采用距離判別兩個向量之間的接近程度，設兩個樣本總體U,Z均含N個樣本(矩陣中一行為一個樣本，一個樣本有n′個元素)

(19)

(20)

向量Ui到Z的馬氏距離為

(21)

樣本總體U到Z的馬氏距離為

(22)

為了方便損傷指標的說明，采用下標0,r,t分別表示基準狀態、參考狀態和檢驗狀態，其中前兩個狀態為未損傷狀態。

3.1 矩陣T變換后A矩陣損傷指標

取基準狀態一個樣本計算A矩陣， 分解得Ψ矩陣，再取各狀態N個樣本分別計算A矩陣， 分解到Ψ2矩陣， 由式(13)、式(11)可求得變換后的A矩陣， 取A矩陣第一行元素構造損傷指標，如式(23)。

(23)

(24)

3.2 矩陣T1變換后A矩陣損傷指標

取變換后A矩陣的后m行矩陣[CA2Q1CA2Q2]， 記為α矩陣構造損傷指標， 第i行元素X(X=0,r,t)狀態的樣本為

(25)

故第i個測點的損傷指標為

(26)

全部測點的損傷指標為

IT1=[I1I2…Im]T

(27)

3.3 模型殘差損傷指標

模型殘差損傷指標[11]用于與本文損傷指標進行對比，計算方法如下

(28)

(29)

(30)

X(X=0,r,t)狀態的樣本為

eX=[eX1eX2…eXN]T

(31)

這樣可以得到預測誤差序列，根據卡爾曼濾波的特點，剛開始ek將出現明顯的振蕩，但很快將趨向于穩定，故文中采用去除前40%數據的方法按下式計算損傷指標

(32)

4 損傷指標應用

4.1 數值算例

如圖1所示的三自由度懸臂結構，材料彈性模量為Ec=3.2×104MPa， 柱的截面尺寸為50 mm×50 mm，采用柱彈性模量降低模擬損傷，柱底施加白噪聲激勵。用有限元法模擬結構響應，按100 Hz采樣頻率提取各節點加速度信號，1號測點響應如圖2，各狀態下測點均采集10 240個數據。

圖1 三自由度模型(cm)Fig.1 Three degrees of freedom model (cm)

圖2 測點1加速度時程Fig.2 Acceleration time history of measuring point 1

4.1.1 模型殘差損傷指標Ie

分析時數據Hankel矩陣的列數取4 096，塊行數取50，階數取6。下柱損傷10%時，損傷指標如圖3，圖中前一半數據為參考狀態樣本到基準狀態樣本的距離Dr， 后一半為檢驗狀態樣本到基準狀態樣本的距離Dt。下柱各損傷狀態下的損傷指標如圖4，損傷程度和指標值在損傷程度小于50%時基本呈線性關系。

圖3 下柱損傷10% Ie指標Fig.3 Ie index of bottom column damage 10%

圖4 下柱損傷Ie指標Fig.4 Ie index of bottom column damage

4.1.2 矩陣T變換后A矩陣損傷指標IT

參數設置同4.1.1，共分析15個樣本。各柱分別損傷10%時，損傷指標如圖5，可見各樣本值基本一致，損傷后馬氏距離明顯變大，各柱損傷IT指標均能正確進行損傷判斷。下柱各損傷狀態下的損傷指標如圖6，損傷程度和指標值基本呈線性關系，識別效果優于模型殘差指標。

圖5 各柱損傷10% IT指標Fig.5 IT index of each column damage 10%

圖6 下柱損傷IT指標Fig.6 IT index of bottom column damage

4.1.3 矩陣T1變換后A矩陣損傷指標IT1

各柱分別損傷10%時，1～3測點的損傷指標如圖7，1測點損傷指標明顯比2、3測點大，與下柱損傷吻合(圖7(a))，中柱損傷時1、2測點指標值較接近，大于3測點(圖7(b))，上柱損傷時2測點指標值最大，3測點大于1測點，與損傷情況較為吻合(圖7(c))。下柱各損傷狀態下的損傷指標如圖8，指標能較好識別出損傷位置。

圖7 各柱損傷10% IT1指標Fig.7 IT1 index of each column damage 10%

圖8 下柱損傷IT1指標Fig.8 IT1 index of bottom column damage

4.2 模型斜拉橋橋塔損傷試驗

模型橋按1∶15的幾何相似比制作[18]，主梁采用C50混凝土，索塔、橋墩采用C30混凝土，斜拉索采用高強鋼絲，總體布置見圖9。上塔柱截面如圖10，上斜塔柱頂與上塔柱無連接，為懸臂構件，故上塔柱在橫橋向相當于懸臂柱。

圖9 模型橋總體布置(mm)Fig.9 The overall arrangement of the model bridge (mm)

圖10 上塔柱截面(mm)Fig.10 The upper tower cross section (mm)

對模型橋一側橋塔進行了損傷試驗，該塔柱預應力大小為315 kN，損傷程度通過在塔頂施加不同大小的水平荷載進行控制，加載裝置如圖11，通過拉力計控制荷載大小，水平力距橋面6.55 m，分三個工況加載，荷載大小分別為8 kN、12 kN、16.5 kN，荷載-位移曲線如圖12。裂縫如圖13所示，裂縫主要集中在橋面以上3.5 m范圍內，工況一～工況三的最大裂縫寬度(卸載后未閉合裂縫寬度)分別為0.09(0) mm、0.15(0) mm、0.22 (0.06) mm。

圖11 塔損傷加載裝置(mm)Fig.11 Experimental setup of tower damage (mm)

圖12 荷載位移曲線Fig.12 Load displacement curve

圖13 裂縫圖(mm)Fig.13 Crack map (mm)

4.2.1 模態參數識別

橋塔模態測試傳感器布置如圖14，按1～4測點、4～7測點分組測試，共7個測點，為了增加橋塔橫橋向的振動量，使用力錘在塔頂(8號點)施加一定的激勵，加速傳感器采樣頻率109 Hz，每個工況進行三次測試，每次測試采集10 240個數據，7測點未損傷狀態的數據如圖15。

各工況下識別的基頻變化如圖16，三次測試樣本識別的頻率結果接近，頻率隨著損傷程度的增加而降低，工況三時基頻降低4.4%。

圖14 塔測點布置Fig.14 Tower measuring point arrangement

圖15 測點7響應信號Fig.15 Response signal of measuring point 7

圖16 一階頻率變化Fig.16 First order frequency variation

4.2.2 模型殘差損傷指標Ie

分析時數據Hankel矩陣的列數取4 096，塊行數取50，階數取6。工況一由1～4測點數據識別的損傷指標如圖17，指標值接近1，無法判別有無損傷。各工況的損傷指標如圖18，可見，兩組測點各工況下的指標值均很小，不能判斷結構有無損傷。說明對于變化幅度大的響應信號(如突加激勵影響)模型殘差損傷指標幾乎不能進行損傷識別。

圖17 工況一1～4測點Ie指標Fig.17 Ie index of 1～4 measuring points of case 1

圖18 Ie指標Fig.18 Ie index

4.2.3 矩陣T變換后A矩陣損傷指標IT

參數設置同4.2.2，共分析15個樣本。工況一由1～4測點數據識別的損傷指標如圖19，各工況的損傷指標如圖20，指標值明顯大于1，可判斷結構存在損傷，隨著損傷程度的增加，指標值增加，1～4測點損傷指標值大于4～7測點。

圖19 工況一1～4測點IT指標Fig.19 IT index of 1～4 measuring points of case 1

圖20 IT指標Fig.20 IT index

4.2.4 矩陣T1變換后A矩陣損傷指標IT1

工況一由1～4測點數據識別的損傷指標如圖21，工況一1～4測點、4～7測點兩組數據識別的損傷指標如圖22，通過4測點，將后3個測點的損傷指標與前四個測點統一起來，各工況的損傷指標如圖23，可見，損傷區域的測點損傷指標值較未損傷區域的大，可以較好的識別出上塔柱下部區域發生損傷。

圖21 工況一1～4測點IT1指標Fig.21 IT1 index of 1～4 measuring points of case 1

圖22 工況一IT1指標Fig.22 IT1 index of case 1

圖23 IT1指標Fig.23 IT1 index

5 結 論

本文對基于隨機子空間方法的損傷識別進行了試驗研究，為了克服離散狀態空間矩陣A的多樣性，構造了兩個變換矩陣T、T1，通過變換矩陣T可以將任意狀態的A矩陣轉化為某一特定狀態(損傷指標記為IT)；通過變換矩陣T1使得離散輸出矩陣C化為單位矩陣和零矩陣的組合，變換后的A矩陣前m行(m為測點數)由零矩陣和單位矩陣組成，后m行與測點相對應，可用來進行損傷定位(損傷指標記為IT1)。損傷指標采用統計方法確定，由損傷前后樣本馬氏距離計算得到。

通過一三自由度懸臂結構仿真模型和一座模型斜拉橋橋塔損傷試驗，分析結果表明，模型殘差損傷指標Ie能用于白噪聲激勵響應信號的損傷識別，無法對變化幅度大的響應信號進行損傷識別；IT指標對兩類信號均能夠較好的識別有無損傷及相對損傷程度，損傷程度與指標值基本呈線性關系；IT1指標可以對損傷進行定位。故通過IT、IT1指標的聯合應用，可對懸臂結構進行相對損傷程度和損傷位置的識別。對其它結構類型，指標的損傷識別效果有待進一步的檢驗。