(福建工程學院 土木工程學院， 福建 福州 350118)

為了保障安全性，許多重大工程結構通過安裝結構健康監(jiān)測系統(tǒng)來實時監(jiān)測其服役期間安全狀況[1-5]。結構損傷識別是結構健康監(jiān)測的重要內容，近30年，基于動力特性的損傷監(jiān)測研究有了重大的進展，其基本原理為通過結構損傷前后結構模態(tài)信息的改變進行結構損傷診斷，采用的方法主要有模態(tài)應變能法[6]、模態(tài)置信度判據(jù)法[7]、曲率模態(tài)法[8]、剛度法[9]、柔度法[10]等。其中，基于有限元模型的結構損傷識別方法主要通過結構實測振動響應數(shù)據(jù)提取的模態(tài)參數(shù)或損傷特征參數(shù)，結合有限元模型，識別結構系統(tǒng)參數(shù)，然后進行損傷識別[11]。正則化方法是目前求解動力學反問題的主流方法，如Tikhonov正則化[12]和奇異值分解法[13]，等。這些方法的原理是基于2范數(shù)的正則化算法，存在對識別結果過度光滑的效果。而結構局部損傷的特點是突變性和稀疏性，因此可將L1范數(shù)正則化方法運用到求解結構損傷識別中的動力學反問題，提高識別精度，改善求解過程的不適應性[14-19]。

為此，本文利用結構模態(tài)振型，引入L1正則化，提出一種基于模態(tài)振型和L1正則化的結構損傷識別方法。

1 結構損傷識別方法

1.1 特征方程

無損結構用特征值方程描述結構的模態(tài)特性[14]，即：

(1)

式中，K和M分別為整體剛度矩陣和質量矩陣；ωi和Φi是無損結構的第i階固有頻率和第i階振型；n為無損結構的振型總數(shù)。

當結構發(fā)生損傷時，其特征方程為：

(2)

1.2 單元損傷指標

以單元的剛度減少來定義損傷即：

(3)

(4)

1.3 單元損傷方程

(5)

又知單元剛度矩陣與整體剛度矩陣的關系，如下：

(6)

(7)

式中，N是單元總數(shù)。子矩陣ke是第e個單元對結構整體剛度矩陣的貢獻。

將式(4)代入式(7)，則：

(8)

同理，第i階振型為：

(9)

假設結構的質量矩陣在損傷前后保持不變。這種假設在大多數(shù)實際應用中是可行的。取式(5)的轉置，將式(1)(8)(9)代入，得到的損傷方程如下：

i=1,…,n且j=1,…,m

(10)

式(10)可以表示為：

[S]{Δa}={ΔR}

(11)

式中靈敏度矩陣[S]的系數(shù)為：

(12)

右側殘差向量ΔR為：

(13)

2 稀疏正則化優(yōu)化算法

設原始信號為{x}，其長度為l，假設通過線性映射[17]得到長度為k的向量{y}：

{y}=[Θ]{x}

(14)

式中，[Θ]稱為測量矩陣或感知矩陣，其維數(shù)為k×l。假設{x}中只有s項(s≈k?l)非零，其他n-s項為零，當k?l時，未知數(shù)遠多于方程數(shù)，故式(14)為欠定方程組。

根據(jù)稀疏恢復理論，{x}可以通過式(15)優(yōu)化求解：

(15)

利用l1范數(shù)代替l0范數(shù)，將式(15)的非凸組合優(yōu)化問題轉化為凸松弛問題求解：

(16)

實際應用中，多數(shù)信號近似稀疏，而且信號測量過程中由于硬件設備等原因會引入噪聲，當測量存在噪聲時，壓縮感知測量過程表示為：

{y}=[Θ]{x}+e

(17)

式中，e為測量誤差。

x可通過求解如下l1優(yōu)化問題進行重構：

(18)

由式(11)可進一步得到[15,16]：

[S]{Δa}={ΔR}={R0}-{RE}

(19)

式中，{R0}為無損結構測得的特征值，{RE}為損傷結構測得的特征值。

式(19)建立了結構單元損傷剛度變化系數(shù)向量{Δa}與測試的結構模態(tài)參數(shù)直接的關系，但{Δa}的維數(shù)遠大于測點數(shù)目，即未知數(shù)遠多于方程數(shù)，式(19)為一病態(tài)的線性方程組，考慮到結構單元損傷向量{Δa}具有稀疏性，只有少數(shù)損傷的結構單元其值為非零，因此可通過式(20)優(yōu)化求解：

(20)

將式(19)代入(20)中：

(21)

式中，R({Δa})=[S]{Δa}+{RE}。

式(21)可轉化為無約束最小化問題：

(22)

式中，拉格朗日乘子β為正則化參數(shù)且大于0。

一般采用優(yōu)化方法求解式(22)，傳統(tǒng)的模型修正方法是采用L2正則化方法，即：

(23)

采用L2正則化方法常常導致結果過于光滑，這與實際損傷的稀疏性不符，為此，本文采用式(22)作為收斂方程，且考慮到損傷值的非負性，增加限定條件0≤Δa≤1。

3 數(shù)值算例分析

3.1 簡支梁數(shù)值模擬

建立簡支梁有限元模型，簡支梁總長6 m，含15個單元，彈性模量E=32 GPa，密度ρ=2 500 kg/m3。模型如圖1所示。

為了研究方法的有效性，設置種損傷工況，工況模擬見表1。其中，結構損傷主要根據(jù)模擬剛度降低。為了研究噪聲的影響，本文在工況3中對振型添加了1%的噪聲水平，即：

z=za×(1+εR)

(24)

式中，za和z分別代表原始信號和噪聲污染后的振型；R是均值為0，偏差為1的正態(tài)分布隨機數(shù)；ε是噪聲水平為1%。

正則化參數(shù)β采用L曲線法確定，取0.1。

圖1 簡支梁模型Fig.1 Model of the simply-supported beam

表1 損傷工況模擬

3.2 計算結果

取n=4個，m=3個((n×m)

3.3 振型數(shù)對識別影響分析

為進一步研究無損結構的振型數(shù)n、損傷結構的振型數(shù)m以及二者的總振型數(shù)對損傷識別結果的影響對3種損傷工況展開分析。

圖2 3種損傷工況下的損傷識別結果Fig.2 Identified damage results of three cases

表2為無損結構的振型數(shù)n∈[1,6]，損傷結構的振型數(shù)m∈[1,6]時工況1的識別結果。從表2可以看出，針對單損傷工況，無損結構和損傷結構的振型數(shù)對識別結果均無影響，基于L1正則化的損傷識別方法均可精確地識別結構損傷。

表3為損傷工況2在不同振型數(shù)下的損傷識別結果。從表3可見，當n×m為1×1，1×2，1×3時，單元2識別值的相對誤差達到100%，識別錯誤。分析原因為n與m的值均較小時，提供信息較少，且損傷單元有2個，工況較為復雜，故識別結果較差。對于剩余振型數(shù)，單元2和單元8的識別誤差最大分別為0.5%和16.5%，均小于20.0%，且大多數(shù)情況相對誤差為0，故本文方法可以較準確地識別出損傷位置與損傷程度。

表4為損傷工況3在不同振型數(shù)下的損傷識別結果。從表4可知，當損傷結構的振型數(shù)為1或無損結構的振型數(shù)為1時，識別誤差為100%，但當n×m>4時，不再出現(xiàn)識別值為0的情況，隨著n和m的不斷增加，識別精度逐漸提高。當n×m≥6時，單元2檢測值的最大誤差為0.8%，單元8檢測值的最大誤差為0.5%，單元12檢測值的最大誤差為2%，均小于20%，處于可接受的范圍。故建議通過本文所提方法進行結構損傷識別時，n×m可取大于等于6。

表2 損傷工況1損傷識別結果

表3 損傷工況2損傷識別結果

表4 損傷工況3損傷識別結果

4 結語

1)對于單損傷工況和多損傷工況，基于模態(tài)振型和L1正則化的方法均可有效地定位結構損傷和量化損傷程度。

2)隨著所采用的結構損傷模型與有限元模型的振型數(shù)的增加，結構損傷識別精度有所提高；但較少的模態(tài)數(shù)據(jù)(模態(tài)振型數(shù)為1時)，識別效果不理想，甚至出現(xiàn)誤判現(xiàn)象。

3)損傷單元數(shù)導致所需的結構損傷模型與有限元模型的振型數(shù)均不同，總體而言，當n×m≥6，可以達到實際應用要求。