基于馬氏距離累積量和EMD的結構損傷識別兩步法

陳 闖， 俞 鵬， 王銀輝

(1.浙江大學 寧波理工學院，寧波 315100；2.浙江公路技師學院，杭州 310023)

結構損傷識別一直是結構安全狀態評估的熱點問題，隨著結構健康監測技術的發展和應用，大量、不同類型的結構監測信息為準確、有效的識別結構損傷提供了數據支撐[1]。然而，對于整體結構而言，損傷時刻和損傷位置往往具有隨機性，尤其對于早期、微小損傷以及遠離監測位置的損傷，監測數據中包含的損傷信息少且信噪比低，這極大的增加了結構損傷識別的難度[2]。

目前，國內外學者基于結構監測數據已提出了一系列的損傷識別方法，主要是包括時域、頻域以及時頻域的損傷識別方法。基于馬氏距離的時域損傷識別方法以其較好的適用性及敏感性在結構損傷識別領域已有較多應用。Mosavi等[3]以馬氏距離為損傷位置評判指標提出一種在環境激勵下鋼梁橋損傷位置的識別方法，工字鋼模型試驗驗證了該方法對損傷具有很好的敏感性。曹軍宏等[4]提出一種基于獨立分量分析結合馬氏距離判斷結構損傷的方法，鋼框架模型的沖擊試驗結果驗證了該方法的有效性。駱志高等[5]利用不同程度裂紋對應的馬氏距離作為標準值，通過對比測試組馬氏距離與標準值的差異判斷損傷程度。Zhou等[6]通過馬氏距離的最大值評判強迫振動下結構的損傷程度，并與頻率響應函數為指標的損傷識別方法進行比較，證明了該方法的損傷識別能力。劉綱等[7]利用馬氏距離作為損傷特征值用于評價自回歸模型的損傷敏感向量，結果表明該方法具有較好的抗噪性。王清華等[8]采用電磁干擾技術結合馬氏距離指標對木結構進行了損傷程度評估。相比于直接利用采集的數據構造馬氏距離，Szymon等[9]利用結構的加速度響應構造的Hank矩陣計算馬氏距離，并與子空間損傷識別方法進行了比較，證明了兩種方法的高效性，并能夠降低風、溫度等周圍環境對識別效果的影響。羅鈞等[10]利用自回歸滑動平均模型的系數向量構造馬氏距離對剪切型框架結構進行了損傷定位研究。

綜上所述，馬氏距離作為結構損傷評判指標的適用性及有效性已被證實。然而，對于結構健康監測中信噪比低、包含的損傷信息較少情況下的損傷識別研究較少。針對以上情況，本文基于馬氏距離累積量和經驗模態分解提出了一種結構損傷識別定位的“兩步法”，數值模擬和模型試驗結果驗證了該方法的有效性。

1 結構損傷識別“兩步法”的基本原理

1.1 馬氏距離基本原理

1936年印度統計學家Mahalanobis提出了馬氏距離的基本原理[11]，其主要目的是比較兩樣本(集合)的差異性。馬氏距離的特點是能夠考慮樣本中元素的大小和特性之間相關性，同時不區分數據類型，其計算公式為

(1)

1.2 經驗模態分解(EMD)基本原理

(2)

式中，ci為第i階IMF，rn(t)為殘差項。

經驗模態分解能夠將非平穩信號按照頻率高低排列成各階本征模態函數，結構損傷導致的剛度降低將引起結構響應中某些頻帶能量的改變，從而造成頻段能量的轉移，這種頻帶能量的轉移對結構的損傷具有較好敏感性[13-14]。

1.3 損傷識別向量的構造

結構發生損傷通常導致結構動力參數的改變，然而，較小的參數變化很難有效識別。根據結構損傷的“單一向”特性(即一旦損傷發生，產生的結構性裂縫只能一直存在且持續增加)，可通過損傷信息的“疊加累積”形成一種“放大效應”，使微小損傷有效識別，且信噪比低的信號有效提取損傷信息。本文通過將監測數據樣本的馬氏距離進行累積形成損傷信息的“放大效應”，從而實現損傷信息的累積。

假設監測數據樣本長度為l，等分成k組，每組數據長度為j，即l=j×k。定義第i組馬氏距離累積量(Mahalanobis distance cumulates, MDC)為

(3)

式中，xn為待測樣本中的各類數據，μ為參考樣本的均值估計，S為參考樣本的協方差估計；MDCi為第i組馬氏距離累積量。損傷識別向量Rk為各組MDC值組成的向量，表示為

(4)

1.4 結構損傷識別的“兩步法”

本文提出的結構損傷識別“兩步法”主要是針對由于損傷較小或傳感器測點距離損傷較遠而導致的監測信息中包含的損傷信息較少且信噪比低的情況，首先需要對原始的監測數據進一步處理，即對原始監測數據進行經驗模態分解，進而在更精細化的頻段內進行損傷識別。具體的步驟為

步驟1 利用損傷前后直接監測的物理參量(加速度、應變等)作為參考樣本和待測樣本計算MDC值構造損傷識別向量，并以參考樣本損傷識別向量MDC值的均值估計作為閾值，即

(5)

式中，mean(Rk)為Rk向量的均值估計。當采集到的物理參量信噪比較高，信號中包含了較多的損傷信息時，往往MDC值的變化已經很明顯，則僅需第一步即可識別出裂縫引起的剛度降低。

步驟2 對原始結構監測數據進行EMD分解，獲得各階IMF，以損傷前后的各階IMF作為參考樣本和待測樣本構造損傷識別向量，并建立各階IMF樣本的MDC值概率密度函數，以參考樣本MDC值概率密度函數95%置信區間的上限值作為損傷識別的閾值，并通過待測樣本與閾值之間的對比關系進一步識別結構損傷。

本文提出的結構損傷識別兩步法，在數值模擬過程中，微小的剛度變化經過第一步方法即能夠識別出損傷；然而，在模型試驗和實際工程中由于噪聲、溫度等影響，為損傷的精確識別增加了難度，距離損傷位置較遠處的傳感器或者損傷較小的情況下的MDC值變化不明顯，因此，在直接利用監測數據無法明確識別出損傷的情況下，則進入第二步基于原始監測數據EMD分解后，利用各階IMF進行更精細化的損傷識別。

2 數值模擬算例

2.1 有限元模型

采用有限元分析軟件ANSYS14.0建立跨徑為10 m的簡支梁模型，主梁采用BEAM4單元模擬，截面高和寬分別為0.6 m和0.3 m，慣性矩Izz和慣性矩Iyy分別為0.001 35 m4和0.005 4 m4，彈性模量為2.06 GPa，密度為7 850 kg/m3，泊松比為0.3。主梁共劃分成25個單元和26個節點，每個單元長度為0.4 m。在節點1和節點26處設置支座，簡支梁單元和節點編號如圖1所示。

圖1 簡支梁有限元模型(cm)

在主梁13號節點施加F=100×Sin(2t)的正弦力激勵，其振幅為100 kN，荷載激勵時間步長為0.1 s，數據采集時長為500 s。目前，通常采用降低單元剛度的形式來模擬單元的損傷，本文通過降低9號單元剛度10%和20%兩種工況模擬結構損傷。圖2為3號節點和10號節點損傷前后的加速度時程曲線，可以看出，損傷前后加速度變化不大，單從加速度時程曲線上很難識別出結構損傷。

(a) 3號節點加速度時程曲線

(b) 10號節點加速度時程曲線

圖2 損傷前后簡支梁節點加速度時程曲線

Fig.2 Acceleration time-history curves of nodes of the simple support beam before and after damage

2.2 損傷位置識別

按照上文1.2節的方法對加速度數據構造馬氏距離累積量MDC值，首先，分別以健康狀態下的加速度數據作為參考樣本，損傷后的加速度數據作為待測樣本，將數據按照j=100進行分段，計算參考樣本和待測樣本的MDC值變化如圖3所示。

(a) 損傷前后3號節點MDC值變化圖

(b) 損傷前后10號節點MDC值變化圖

Fig.3 MDC value variation node 3 and node 10 before and after damage

從圖3中可以看出，當結構出現損傷后，3號和10號節點的MDC值明顯的增大，通過比較健康狀態與損傷狀態的MDC值即可識別結構是否發生了損傷。隨著單元損傷程度的加深(9號單元剛度降低10%到降低20%)，MDC值的整體增加愈加明顯。所以，當采集數據包含了較多損傷信息或者其損傷較為明顯時，其MDC值的變化較為明顯，可以直接通過對比MDC值的變化即可識別出結構的損傷。圖4為各個節點在損傷前后MDC值的均值，通過比較均值的變化可判斷結構損傷的位置。

圖4(a)表示損傷前后各個節點處MDC的均值，可以看出，結構損傷導致了MDC均值明顯增大，且隨著損傷的加深，MDC均值變化愈加明顯。圖4(b)為損傷前后各個節點相較于健康狀態MDC值均值的變化情況，與圖4(a)變化一致，發生損傷單元即9號節點處MDC均值的變化量較其他單元更為明顯，形成了明顯的峰值。因此，可通過各節點處MDC均值的變化情況來判斷損傷的位置。

2.3 抗噪性分析

噪聲是工程試驗和現場監測過程中無法避免的，主要是由于環境激勵、測試儀器與傳感器本身等因素所導致的，結構在動力測試的過程中，一般會將環境等產生的噪聲作為一種激勵結構振動的輸入，則其輸出信號中包含了原測試信號與噪聲。通常采用信號與噪聲的相對比例，即信噪比描述監測數據中噪聲對監測數據的影響。信噪比越小，說明輸出信號中的噪音越大，反之亦然。現對10號節點健康狀態、單元損傷10%與損傷20%情況下結構的加速度信號中加入5%高斯白噪聲，并對其損傷識別情況進行分析，以驗證該方法的抗噪性。圖5分別為10號節點在加入了5%的高斯白噪聲后的加速度時程損傷前后的對比圖。可以看出，通過對比損傷指標向量MDC值均值變化依然可以較好的識別出結構是否發生損傷。

(a) 各節點MDC均值圖

(b) 各節點MDC均值相對變化圖

Fig.4 MDC average values and relative variation of MDC average values of each node of the simply supported beam

3 模型試驗方法驗證

結構在實際監測環境中，往往會受到如風、溫度等外界環境的影響；加之，結構損傷具有隨機性，傳感器與損傷的相對位置往往決定了是否能夠收集到足夠的損傷信息；此外，初始損傷信噪比低，也極大的增加了損傷識別的難度。本節對工字鋼進行了損傷識別的模型試驗，從梁底向上切割工字鋼梁人為制造損傷，距離損傷遠近不同的位置分別布置了應變傳感器，以驗證監測位置距離損傷較遠、信噪比較低情況下的損傷識別效果。

3.1 工字鋼模型及損傷試驗

選取一根14#工字鋼進行模型試驗研究，其橫截面尺寸為140 mm×80 mm×5.5 mm(高度×寬度×腹板厚度)，質量密度為16.89 kg/m，標準跨徑為2.3 m，計算跨徑為2.0 m。

(a) 加入5%白噪聲后損傷前后加速度對比圖

(b) 10號節點MDC值變化圖

圖5 簡支梁10號節點抗噪性分析

Fig.5 Anti-noise analysis of the node 10 of the beam

工字鋼底部等間距布置5個斷面的應變測點，每個斷面2個應變測點，測點編號為y1～y10；布置環境溫度傳感器T1及工字鋼表面溫度傳感器T2和T3，傳感器布置及現場試驗采集設備如圖6和圖7所示。利用激振器在工字鋼接近跨中位置進行強迫振動激勵，采樣頻率50 Hz，采集健康狀態下24 h的動態應變數據。在工字鋼底部向上切割5 mm和10 mm損傷，分別定義為損傷工況a和損傷工況b，并繼續進行強迫振動，保持50 Hz的采樣頻率，分別各采集24 h應變數據。

圖6 應變片與溫度傳感器布置圖

圖7 工字鋼試驗加載及數據采集

3.2 趨勢項分析

通過24 h的動態應變數據采集，獲得在健康狀態下工字鋼底部的應變時程曲線如圖8所示，可以看出，應變時程具有明顯的趨勢項，采用滑動平均值法[15]獲得整體趨勢項及剔除趨勢項后應變。圖9為提取的應變趨勢項，與圖10測試得到溫度變化對比可以看出，其應變趨勢項與溫度變化具有較高的相關性，表1列出了應變趨勢項與溫度的相關性分析結果，可以看出，各測點的應變趨勢項與工字鋼表面溫度的相關系數均較高。因此，可以認為應變趨勢項為溫度趨勢項，在數據分析時將溫度趨勢剔除，只分析在外荷載激勵下的受迫振動數據。

圖8 應變時程曲線

3.3 損傷識別分析

以健康狀態下剔除趨勢項后的應變數據為參考樣本，損傷工況a和損傷工況b下剔除趨勢項后的應變數據計為待測樣本。參考樣本和待測樣本均按照j=10 000進行分段，分別計算損傷識別向量的MDC值，各個位置損傷前后MDC值的變化情況如圖11所示。

圖9 應變趨勢項時程曲線

圖10 大氣溫度及工字鋼表面溫度時程曲線

Fig.10 Time-history curves of the atmospheric temperature and the surface temperature of thei-beam

表1 應變趨勢項與溫度相關系數

Tab.1 Correlation coefficients between the strain trend items and the temperature

相關系數應變測點y1y3y5y7y9溫度測點T1T2T30.710.930.850.700.940.850.690.940.850.710.920.830.710.930.85

通過圖11可以看出，測點y3和y5的MDC值相較于其他傳感器變化更為明顯，MDC值均值也隨著損傷的發生而逐漸增大，并且隨著損傷的加重，其MDC均值的增幅更大，可判斷損傷位置在測點y3和y5之間，這與實際損傷的損傷位置相吻合。因此，可直接通過損傷識別過程中的第一步，即可識別出損傷及損傷位置。

3.4 經驗模態分解的應用

為了進一步驗證本文提出方法在距離損傷位置較遠處和直接利用MDC值很難判斷是否損傷發生情況下的識別效果，選取測點y1和y3兩位置處的應變結果進行分析。通過G.Rilling的EMD工具箱對各測點的應變數據進行分解得到各階的本征模態函數，測點y3健康狀態下的EMD分解后第3階本征模態函數如圖12所示。

(a) 測點y1損傷前后MDC值變化

(b) 測點y3損傷前后MDC值變化

(c) 測點y5損傷前后MDC值變化

(d) 測點y7損傷前后MDC值變化

(e) 測點y9損傷前后MDC值變化

圖11 工字鋼部分應變測點損傷前后MDC值變化

Fig.11 MDC values of partial testing strain points of thei-beam before and after damage

圖13為y3應變剔除趨勢項后信號構造了損傷識別向量MDC值的概率密度函數，分別采用Burr分布、Normal分布和Logistic分布擬合概率密度函數得到的結果如表2所示。采用極大似然法對3種概率密度函數擬合結果對比發現，Burr分布擬合的極大似然值最大，因此，可采用Burr分布對MDC值進行分布擬合，其95%置信區間為[963.3,1 045.5]，并將其置信區間上限作為損傷識別向量的閾值。

表2 概率密度函數擬合結果對比

Tab.2 Comparison results of the probability density function fitting

擬合方式Burr分布Normal分布Logistic分布極大似然值-4 375.67-4 478.33-4 429.21

圖12 健康狀態下測點y3經驗模態分解的IMF3

圖13 健康狀態下測點y3 MDC值概率密度函數擬合圖

Fig.13 Probability density function fitting of MDC values of testing pointy3in health condition

對測點y1與y3的應變數據利用EMD分解得到各階IMF，再分別對各階IMF構造損傷識別向量以判斷損傷情況，因篇幅有限，只列出前四階IMF構造的MDC值概率密度函數擬合結果如圖14和圖15所示。

(a) IMF1損傷前后MDC值概率密度函數擬合

(b) IMF2損傷前后MDC值概率密度函數擬合

(c) IMF3損傷前后MDC值概率密度函數擬合

(d) IMF4損傷前后MDC值概率密度函數擬合

圖14 測點y1前四階IMF損傷前后MDC值概率密度函數擬合

Fig.14 Probability density function fitting of MDC values of first four order IMFs of the testing pointy1before and after damage

圖14僅列出了測點y1經EMD分解后前四階IMF損傷前后的MDC值概率密度函數擬合結果。未發生損傷狀態下，各階IMF損傷識別向量MDC值的概率密度函數應較為一致；一旦結構發生損傷，某些頻率IMF構造的損傷識別向量MDC值的概率密度函數面積將顯著超過閾值，在統計學角度上表明超過閾值的MDC值的概率顯著增加，超過閾值面積越大，損傷程度也越大。本例中測點y1利用分解的IMF3、IMF4、IMF6和IMF7信號構造損傷識別向量后，都較好的識別出結構是否發生損傷，其他階IMF的MDC值變化不明顯，且IMF3和IMF4相對于其他階本征模態函數而言，具有更好的識別效果，其MDC值的概率密度函數隨著損傷程度的增加變化更為明顯。這是由于結構損傷導致了結構整體剛度的改變，在輸入能量基本不變的情況下，各階IMF頻帶能量發生了重分布，各階頻率對于損傷的敏感程度也各不相同。圖15是y3應變經EMD分解后損傷前后前四階IMF的MDC值概率密度函數擬合結果，測點y3的損傷識別情況與測點y1的較為相似，其IMF3、IMF4和IMF8信號構造損傷識別向量后，可以較好的識別出結構的損傷情況，而IMF3和IMF4相較于其他幾階IMF，更好的識別出結構是否發生損傷，且其對損傷程度敏感度明顯大于其他階本征模態函數。

(a) IMF1損傷前后MDC值概率密度函數擬合

(b) IMF2損傷前后MDC值概率密度函數擬合

(c) IMF3損傷前后MDC值概率密度函數擬合

(d) IMF4損傷前后MDC值概率密度函數擬合

圖15 測點y3前四階IMF損傷前后MDC值概率密度函數擬合

Fig.15 Probability density function fitting of MDC values of first four order IMFs of the testing pointy3before and after damage

4 結 論

本文提出了一種基于馬氏距離累積量和經驗模態分解的結構損傷識別“兩步法”，首先，直接利用原始監測數據的馬氏距離累積量構造損傷識別向量，并通過馬氏距離累積量的均值作為閾值進行初步判斷；在損傷較小、信噪比低的情況下，通過EMD方法將原始監測數據按頻率高低分解成各階IMF再構造損傷識別向量，并對其進行概率密度函數擬合，以概率密度函數95%置信區間上限值作為閾值進一步識別損傷。通過有限元模型和工字鋼模型試驗對該方法進行了驗證，得出主要結論如下：

(1) 利用原始監測數據的馬氏距離損傷累積量均值作為閾值，在損傷較為明顯，監測數據中包含了較多損傷信息時可直接識別出損傷及損傷位置；隨著損傷程度的增加，損傷識別向量MDC值的變化愈加明顯。

(2) 溫度趨勢項是監測數據中最主要的影響因素，在利用該方法進行實際結構損傷識別時，必須先剔除溫度趨勢項的影響。

(3) 結構損傷造成了結構動力響應中各階頻率能量的改變，損傷識別向量MDC值的概率密度函數能夠有效反映出結構剛度降低引起的頻率能量改變，但損傷前后頻率能量的重分布模式有待于進一步研究。