基于隨機子空間方法的懸臂結(jié)構(gòu)損傷識別研究

唐盛華， 方 志， 張國剛

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

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

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

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

1 數(shù)據(jù)驅(qū)動隨機子空間模型

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

(1)

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

輸出數(shù)據(jù)構(gòu)造的Hankel矩陣為：

(2)

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

數(shù)據(jù)驅(qū)動隨機子空間識別方法采用投影法進行計算：

(3)

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

(4)

Oi=[CCA…CAi-1]T

(5)

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

(6)

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

(7)

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

(8)

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

(9)

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

2 變換矩陣

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

(10)

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

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

(11)

因Λ=Λ2，于是

(12)

即

(13)

故通過變換矩陣T可以將其它狀態(tài)的A矩陣轉(zhuǎn)化為某一特定狀態(tài)。

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

(14)

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

假設(shè)：

(15)

則有：

(16)

于是

(17)

(18)

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

3 損傷指標(biāo)

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

(19)

(20)

向量Ui到Z的馬氏距離為

(21)

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

(22)

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

3.1 矩陣T變換后A矩陣損傷指標(biāo)

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

(23)

(24)

3.2 矩陣T1變換后A矩陣損傷指標(biāo)

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

(25)

故第i個測點的損傷指標(biāo)為

(26)

全部測點的損傷指標(biāo)為

IT1=[I1I2…Im]T

(27)

3.3 模型殘差損傷指標(biāo)

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

(28)

(29)

(30)

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

eX=[eX1eX2…eXN]T

(31)

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

(32)

4 損傷指標(biāo)應(yīng)用

4.1 數(shù)值算例

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

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

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

4.1.1 模型殘差損傷指標(biāo)Ie

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

圖3 下柱損傷10% Ie指標(biāo)Fig.3 Ie index of bottom column damage 10%

圖4 下柱損傷Ie指標(biāo)Fig.4 Ie index of bottom column damage

4.1.2 矩陣T變換后A矩陣損傷指標(biāo)IT

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

圖5 各柱損傷10% IT指標(biāo)Fig.5 IT index of each column damage 10%

圖6 下柱損傷IT指標(biāo)Fig.6 IT index of bottom column damage

4.1.3 矩陣T1變換后A矩陣損傷指標(biāo)IT1

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

圖7 各柱損傷10% IT1指標(biāo)Fig.7 IT1 index of each column damage 10%

圖8 下柱損傷IT1指標(biāo)Fig.8 IT1 index of bottom column damage

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

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

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

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

對模型橋一側(cè)橋塔進行了損傷試驗，該塔柱預(yù)應(yīng)力大小為315 kN，損傷程度通過在塔頂施加不同大小的水平荷載進行控制，加載裝置如圖11，通過拉力計控制荷載大小，水平力距橋面6.55 m，分三個工況加載，荷載大小分別為8 kN、12 kN、16.5 kN，荷載-位移曲線如圖12。裂縫如圖13所示，裂縫主要集中在橋面以上3.5 m范圍內(nèi)，工況一～工況三的最大裂縫寬度(卸載后未閉合裂縫寬度)分別為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 模態(tài)參數(shù)識別

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

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

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

圖15 測點7響應(yīng)信號Fig.15 Response signal of measuring point 7

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

4.2.2 模型殘差損傷指標(biāo)Ie

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

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

圖18 Ie指標(biāo)Fig.18 Ie index

4.2.3 矩陣T變換后A矩陣損傷指標(biāo)IT

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

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

圖20 IT指標(biāo)Fig.20 IT index

4.2.4 矩陣T1變換后A矩陣損傷指標(biāo)IT1

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

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

圖22 工況一IT1指標(biāo)Fig.22 IT1 index of case 1

圖23 IT1指標(biāo)Fig.23 IT1 index

5 結(jié) 論

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

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