基于秩分析的結構損傷識別研究

楊秋偉, 周衛東, 梁超鋒

(1.紹興文理學院 土木工程系,浙江 紹興 312000;2.華匯工程設計集團股份有限公司,浙江 紹興 312000)

重大工程結構在服役期間，由于外力碰撞、環境腐蝕、材料老化等因素的影響，將不可避免的出現損傷。結構的局部損傷將可能導致結構整體的迅速破壞而釀成重大的工程事故。為了確保結構的安全性能，必須對結構損傷的發生、損傷的位置和程度作出及時的判斷。近年來，關于結構損傷識別方面的研究已取得了長足的進展[1-4]。目前大部分損傷識別方法均是基于結構有限元模型的模型修正方法，其原理是通過不斷修正結構的有限元模型使其和測量所得的反應數據相匹配，模型的修正量即為結構的損傷量。模型修正方法通常可以分為四類：矩陣優化方法[5-6]，靈敏度方法[7-9]，特征結構分派方法[10-11]，最小秩擾動方法[12-16]。一般情況下，結構損傷均發生于結構中的少數區域，而最小秩方法中對更新矩陣的秩的限制最符合這種實際情況，因此，最小秩方法有著其獨特的優越性[16]。

本文研究了結構損傷與系統矩陣秩變化之間的對應關系，并聯合利用最小秩和柔度擾動新方法提出了一種很精確的損傷識別方法。和已有的最小秩方法相比，本文方法的創新性在于：

(1) 已有的最小秩方法均屬于動力方法，都利用了結構的振動響應參數；而本文方法則屬于靜力方法，利用了結構在靜力荷載作用下的反應參數。動力方法和靜力方法各有優缺點，Wang等[17]指出了動力方法中所存在的幾個必須解決的問題：①結構的動力反應數據除與結構的剛度有關以外，還與結構的質量和阻尼有關。而結構的損傷多數情況下只造成結構剛度的損失，所以，許多動力方法都忽略了結構阻尼和質量變化對動力測試數據的影響，這顯然與實際工況有偏差；②對于巨大的土木工程結構，難以測量出精確的振動數據；③高階的模態往往對損傷更為敏感，然而實踐中卻只能測量出低價模態數據。相比而言，結構的靜力測試數據只和結構的剛度有關，且一般可以很精確的測量得到，所以靜力方法在土木工程損傷識別領域有著可觀的應用前景。當然，靜力方法也有其相應的缺點：比如需要進行額外的靜力加載工作，某些情況下可能會對結構造成二次損傷，不利于在線測量數據等等。因此，靜力方法和動力方法都各有優缺點，應根據具體的工程實際情況來選用。

(2) 已有的最小秩方法均是求解出剛度矩陣擾動的最小秩解，而本文方法則是直接求出柔度矩陣擾動的限定秩解，它是一種理論上的精確解。本文的研究還表明：損傷前后結構的靜力位移差向量之間必然存在著某種線性相關性，所以所施加的靜力荷載組數并非越多越好，而是只要大于損傷前后柔度變化矩陣的秩即可。文中以一個桁架結構為例對所提方法進行了驗證，結果表明了所提方法的可行性。

1 基于秩分析的損傷識別

1.1 柔度改變量的秩分析

本節首先分析由于結構損傷所導致的柔度改變量的秩變化。結構損傷前后，其剛度和柔度矩陣必然滿足以下關系式：

FK=FdKd=I

(1)

式中F和K是完好結構的柔度矩陣和剛度矩陣(n×n維);Fd和Kd是結構損傷后的柔度矩陣和剛度矩陣;I是n×n維單位矩陣。

一般而言，結構損傷將導致結構剛度減小而柔度增大。令ΔF和ΔK為結構損傷前后柔度和剛度的變化量，即有：

Fd=F+ΔF

(2)

Kd=K-ΔK

(3)

將方程(2)和(3)代入(1)，整理可得

ΔFKd=FΔK

(4)

顯然，方程(4)中矩陣Kd和F均為對稱且滿秩的矩陣。因此，根據矩陣理論，由方程(4)必然有ΔF的秩和ΔK的秩相等，即

rank(ΔF)=rank(ΔK)

(5)

由有限元理論，結構損傷前后剛度變化矩陣ΔK為

(0≤αi≤1)

(6)

式中;Ki是第i個單元剛度矩陣，αi為其相應的損傷參數，N為單元總數。

一般情況下，結構發生損傷僅限于局部少數單元，由方程(6)可知ΔK一般為虧損矩陣，且當發生損傷的單元數目改變時，ΔK的秩也將相應改變(例如：對于桁架結構(可參考算例)，若1個單元發生損傷，則對應rank(ΔK)=1；若2個單元發生損傷，則對應rank(ΔK)=2，其它具體的工程結構亦可類似分析)。再根據方程(5)可知，ΔF一般也是虧損矩陣，且ΔF的秩也將隨著損傷單元的數目的變化而變化(仍以桁架結構為例，若1個單元發生損傷，則對應rank(ΔF)=1；若2個單元發生損傷，則對應rank(ΔF)=2，其它具體結構亦可類似分析)。

1.2 柔度擾動量的限定秩解

本節從結構的靜力響應方程出發，推導損傷前后結構柔度擾動的限定秩解。對于未損傷結構，在已知的外荷載li(n維列向量)的作用下，可以測量或者通過計算得到結構相應的靜力位移ui，即

Kui=li

(7)

方程(7)可以改寫為

ui=K-1li=Fli

(8)

(9)

方程(9)減(7)可得損傷前后靜力位移差Δui為

Δui=ΔFli

(10)

若在結構上分別加載了n個線性無關的荷載向量li(i=1～n)，則可得n個方程，可以將其組合為如下的矩陣方程

ΔU=ΔFL

(11)

式中:矩陣ΔU=[Δu1,Δu2,…,Δun]，矩陣L=[l1,l2,…,ln]。顯然，由于所加的荷載向量是線性無關的，故矩陣L是滿秩矩陣，根據矩陣理論，有

rank(ΔU)=rank(ΔF)

(12)

方程(12)表明，即使測量了很多組靜力位移向量，所得到的靜力位移差向量Δui中只有少數幾個是獨立的(仍以桁架結構為例，若僅一個單元發生損傷(rank(ΔF)=1)，則即使施加了n組荷載向量并測量得到n組位移數據，最終所得的位移差矩陣ΔU的秩仍然為1，即n個Δui中只有一個是獨立的，其它的均可用其線性表示。因此，應用靜力方法進行結構損傷識別時，施加靜力荷載的組數可以不用很多，只要大于損傷前后剛度或者柔度變化矩陣的秩即可。工程實踐中，由于結構損傷前后的剛度或柔度改變量的秩事先是未知的，我們可以采取嘗試的辦法，首先施加一定組數的靜力荷載并計算相應的ΔU的秩(采用計算ΔU的奇異值來確定其秩的方法)，如果計算結果中出現了近似為零的奇異值，則表明所施加的靜力荷載組數是足夠了的。反之則繼續增加靜力荷載組數，直到出現近似為零的奇異值為止即可。由于大部分情況下結構損傷所引起的秩變化都是比較小的，因此通常所施加的靜力荷載組數都不會很多，通過上述嘗試過程很快就可以知道所施加的靜力荷載組數是否夠用。

計算柔度改變矩陣限定秩解的過程如下：假設ΔU的秩為r，則可以從ΔU中取出某個最大線性無關組ΔUr=[Δu1,Δu2,…,Δur](其中r=rank(ΔU))，該組對應的荷載向量矩陣為Lr=[l1,l2,…,lr]，則方程(11)簡化為

ΔUr=ΔFLr

(13)

類似于最小秩理論[12-16]，由方程(13)可得ΔF的限定秩解為

(14)

1.3 損傷識別

當ΔF的限定秩解計算出來以后，接下來可用文獻[9]所提的柔度擾動新方法來計算各單元損傷參數αi(i=1～N)，據此便可判斷結構中哪個單元損傷及其損傷程度。該方法的突出優點在于計算過程不需要迭代運算或高階靈敏度分析，計算量小且精度高，且對于靜定結構而言，該方法是一種精確方法。該方法的主要公式簡述如下(詳細推導過程見文獻[9])：

首先利用單元剛度矩陣的特征值分解與重新組合，可得結構損傷前后的總剛度矩陣的分解形式：

K=CPCT

(15)

Kd=CPdCT

(16)

根據柔度矩陣和剛度矩陣互逆可得

F=(C+)TP-1C+

(17)

(18)

方程(18)減去方程(17)可得

ΔF=EΔBET

(19)

式中:矩陣E稱為柔度聯系矩陣，其計算公式為：

E=(C+)T=(CCT)-1C

(20)

對角矩陣ΔB為

ΔB=diag(β1,β2,…,βN)

(21)

式中:βi為第i個單元的柔度擾動參數，它和剛度擾動參數之間的關系為

(22)

綜上所述，計算損傷參數(即剛度擾動參數)的主要步驟為：首先，根據方程(14)計算出ΔF的限定秩解，然后根據方程(19)計算出各柔度擾動參數ΔB，最后利用方程(22)計算出各剛度擾動參數αi(i=1～N)，根據所得結果便可對損傷情況作出判斷。

2 算例

以圖1所示桁架結構為例，驗證本文所提的損傷識別方法。該結構基本參數為：彈性模量E=200 GPa, 密度ρ=7.8×103kg/m3,單元長度L=1 m和橫截面面積A=7.85×10-5m2。可能的靜力加載點也見圖1中。考慮3種加載方式，方式1：F1=10 kN，F2=0，F3=0；方式2：F1=10 kN，F2=10 kN，F3=0；方式3：F1=10 kN，F2=10 kN，F3=10 kN。顯然，這3種加載方式是線性無關的。假設3種損傷情況：① 單個損傷：單元17剛度損傷20%；② 多個小損傷：單元10和15剛度損傷10%和15%；③ 多個大損傷：單元10和15剛度損傷30%和40%。

對于第1種損傷情況，不考慮測量誤差時三種加載方式下所得的靜力位移差向量列于表1中(限于篇幅，表中只給出了節點7-12所對應的數據，即對應于自由度12-23處的數據)。為說明所得三組位移差向量的線性相關性，各自由度所對應的位移差比值也同時列于表1中。由表1可見，Δu1，Δu2和Δu3成比例，說明這3個位移差向量是線性相關的，只有一個是獨立的，由這3個向量所組成的位移差矩陣ΔU的秩為1。因此，可以任意取一個位移差向量Δui(i=1,2或3)，用方程(14)來計算ΔF的限定秩解，然后應用方程(19)來計算各柔度擾動參數βi(i=1～N)，最后用方程(22)計算出各單元損傷參數αi(i=1～N)，結果見于圖2中。圖2中單元17的損傷程度計算值為α17=0.2，和假設值完全一致。這說明了在不考慮測量誤差的情況下，ΔF的限定秩解是一種理論上的精確解，而柔度擾動新方法對于靜定結構而言，亦是一種精確方法。

圖1 平面桁架結構及其靜力加載

表1 靜力位移差向量 (損傷情況1)

表2 靜力位移差向量 (損傷情況2)

圖2 不考慮測量誤差時的損傷識別結果(損傷情況1)

對于第2種損傷情況，不考慮測量誤差時三種加載方式下所得的靜力位移差向量列于表2中。由表2可見，Δu1，Δu2和Δu3并不完全成比例，為了更好的判斷這組位移差向量的線性相關性，我們通過計算這3個向量所組成的位移差矩陣ΔU的奇異值，根據所得非零奇異值的個數來判斷。表2中Δu1，Δu2和Δu3組成的矩陣ΔU的奇異值為：{1.080 4,0.027 3,0.000 0}。顯然，ΔU的秩為2，這說明Δu1，Δu2和Δu3三個向量中只有兩個是獨立的，這也說明僅兩個單元發生損傷。可從ΔU中取出兩個向量用方程(14)計算出ΔF的限定秩解，然后應用柔度擾動新方法計算出各單元損傷參數，結果示于圖3中，其中單元10和15損傷程度計算值為α10=0.1，α15=0.15，和假設值完全一致。

對于第2種損傷情況，若考慮3%的測量誤差時(測量誤差的添加方法為：在每個精確的位移值的基礎上，加上誤差水平(如3%)乘以一個[-1,1]范圍內的隨機數)，我們可以求出矩陣ΔU的三個奇異值為{1.073 6,0.028 0,0.001 3}，其中第3個奇異值很小可以近似認為等于0，仍然可以認為ΔU的秩應為2，因此取前兩個位移差向量來計算柔度擾動量的限定秩解，然后應用柔度擾動新方法計算出各單元損傷參數見圖4中。另外，為了討論誤差大小對本文方法的影響，添加10%的誤差時的損傷識別結果也同時列于圖4中。由圖4可見，3%噪聲水平下本文方法識別結果和假設值很接近，10%噪聲水平下識別精度有所降低，但仍然可以清楚判斷出損傷單元為10和15。

對于第3種損傷情況，不考慮測量誤差時ΔU的奇異值為{4.158 9,0.103 4,0.000 0}；考慮3%的測量誤差時，ΔU的奇異值為[4.084 7,0.103 4,0.016 0}，考慮10%的測量誤差時，ΔU的奇異值為{3.655 0,0.154 1,0.066 2}。以上三種情況下，均可以認為ΔU的秩應為2，因此采用兩個位移差向量來計算柔度擾動的限定秩解，最后的損傷識別計算結果都列于圖5中。由圖5可見，無誤差時損傷程度計算值為α10=0.3，α15=0.4，和假設值完全一致。3%和10%噪聲水平下也清楚的表明單元10和15發生損傷。因此，由于計算過程中對秩的限制，所提方法對測量噪聲有著較好的魯棒性。

圖5 損傷識別結果(損傷情況3)

3 結 論

本文研究了結構損傷與系統矩陣秩變化之間的對應關系，并聯合利用最小秩和柔度擾動新方法提出了一種很精確的損傷識別方法。所提方法的優勢在于：由于方法中對系統矩陣秩的限制，因此計算結果對測量誤差具有很好的魯棒性。數值算例結果表明：所提方法是合理可行的。

[1]楊秋偉, 劉濟科. 工程結構損傷識別的柔度方法研究進展[J]. 振動與沖擊, 2011, 30(12): 147-153.

YANG Qiu-wei, LIU Ji-ke. Structural damage identification by flexibility change: a review[J]. Journal of Vibration and Shock, 2011, 30(12): 147-153.

[2]程遠勝, 王真. 基于提高頻率靈敏度的結構損傷統計識別方法[J]. 計算力學學報, 2008, 25(6): 844-849.

CHENG Yuan-sheng, WANG Zhen. Statistical damage detection based on frequencies of sensitivity-enhanced structures[J]. Chinese Journal of Computational Mechanics, 2008, 25(6): 844-849.

[3]安永輝, 歐進萍. 大跨度鋼桁橋模型的精細化損傷定位模擬和試驗研究[J]. 計算力學學報, 2011, 28(5): 730-736.

AN Yong-hui, OU Jin-ping. Experimental and numerical studies on precise damage localization of a long-span steel truss bridge model[J]. Chinese Journal of Computational Mechanics, 2011, 28(5): 730-736.

[4]Farrar C R, Lieven N A J. Damage prognosis: the future of structural health monitoring[J]. Philosophical Transactions of the Royal Society A, 2007, 365: 623-632.

[5]Kammer D C. Optimum approximation for residual stiffness in linear system identification. AIAA J, 1988, 26(1): 104-112.

[6]Smith S W, Beattie C A. Secant-method adjustment for structural models[J]. AIAA J, 1991, 29(1): 119-126.

[7]Yu L, Cheng L, Yam L H, et al. Application of eigenvalue perturbation theory for detecting small structural damage using dynamic responses[J]. Composite Structures, 2007, 78: 402-409.

[8]Shi Z Y, Law S S, Zhang L M. Damage localization by directly using incomplete mode shapes[J]. Journal of Engineering Mechanics, 2000, 126(6): 656-660.

[9]Yang Q W. A new damage identification method based on structural flexibility disassembly. Journal of Vibration and Control, 2011, 17(7): 1000-1008.

[10]Zimmerman D C, Kaouk M. Eigenstructure assignment approach for structural damage detection[J]. AIAA Journal, 1992, 30(7): 1848-1855.

[11]Kiddy J, Pines D. Eigenstructure assignment technique for damage detection in rotating structures[J]. AIAA Journal, 1998, 36(9): 1680-1685.

[12]Zimmerman D C, Kaouk M.Structural damage detection using a minimum rank update theory[J]. Journal of Vibration and Acoustics, 1994, 116(2): 222-231.

[13]Zimmerman D C, Smith S W, Kim H M, et al. An experimental study of structural damage detection using incomplete measurements[J]. Journal of Vibration and Acoustics, 1996, 118(4): 543-550.

[14]Doebling S W. Minimum-Rank optimal update of elemental stiffness parameters for structural damage identification[J]. AIAA Journal, 1996, 34(12): 2615-2621.

[15]Souza K D, Epureanu B I. Damage detection in nonlinear systems using system augmentation and generalized minimum rank perturbation theory[J]. Smart Materials and Structures, 2005, 14: 989-1000.

[16]Zimmerman D C. Statistical confidence using minimum rank perturbation theory[J]. Mechanical Systems and Signal Processing, 2006, 20: 1155-1172.

[17]Wang X, Hu N, Fukunaga H, et al. Structural damage identification using static test data and changes in frequencies[J]. Engineering Structures, 2001, 23: 610-621.