基于本征正交分解和卡爾曼濾波的結構損傷識別

楊少沖, 姚 遠, 張 凱, 宋康寧, 劉家亮, 靳佳林, 馬連華, 方有亮

(1. 河北大學 建筑工程學院,河北 保定 071002; 2. 河北大學 質量技術監督學院,河北 保定 071002)

由于各種環境和力學因素的影響,結構易產生一定程度的損傷。損傷傳統上被定義為系統幾何或材料特性的變化[1],對結構的安全性、可靠性和使用壽命產生不利影響。結構振動損傷識別技術是近年來結構健康監測領域重要的研究方向之一[2],也是結構健康監測技術的核心[3],如何實現實時、快速的在線損傷識別,是目前損傷識別技術的關鍵。從已有文獻看,用于土木工程領域的在線損傷識別方法已取得長足的進展,但還處于起步階段,有著較大的進步空間[4]。

本征正交分解(proper orthogonal decomposition,POD)方法是一種較為高效的模型降階技術,該方法能在減少模型自由度的同時,對其主要特征信息進行提取,保證數值模型具有足夠高的精度。該方法通過試驗或數值模擬得到的響應數據,構造用于模型降階的基向量,可以快速、準確地反映出研究對象的主要特征,實現模型降階。由于POD方法能夠在保持高精度的同時從少量的POD模態中提取結構特征,它已被廣泛用于模型降階中[5]。近年來,POD模型降階方法在航空航天[6-7]、計算流體力學[8]、材料力學[9]等領域都得到了廣泛應用,在近期的研究中,Lu等[10]全面綜述了POD方法在不同領域的應用,并對該方法的發展前景進行了展望。在土木工程領域中,Eftekhar Azam等[11]研究了基于POD方法對動力結構系統降階建模的性能,通過對快照矩陣進行奇異值分解(singular value decomposition,SVD)構建降階模型,結果表明該方法得到的降階模型可以達到較高的精度,且優于模態分析。

除在模型降階方面有較高的效率和精度外,基于POD方法得到的結構本征正交模態(proper orthogonal modes, POMs),可以用于表征結構損傷的位置和程度,Galvanetto等[12]較早地利用該方法對簡支梁進行了數值模擬,成功識別出損傷的程度和位置,并在之后的研究中僅利用傳感器數據成功對損傷進行識別[13],而不需要較為精確的數值模型。在國內的損傷識別領域中,楊斌等[14]將奇異值分解與POD方法相結合,通過結構各單元的應力分布進行損傷識別分析,成功識別出結構的損傷及其位置;王丹生等[15]在早期利用POD方法對簡支梁結構的損傷位置進行了識別,驗證了該方法的有效性,并在近期的研究中,將POD方法與粒子群優化算法相結合[16],對剪切型建筑物進行損傷檢測,有效地識別出不同程度的損傷。

而在利用POD方法進行損傷識別分析時,需要利用一段時間的響應構造快照矩陣進行分析,當出現損傷時,需要進行重復的訓練階段[17],因此要對傳統方法進行改進,使其能夠更加快速、準確地對損傷進行識別。卡爾曼濾波法(Kalman filter,KF)方法是一種以最小均方差為準則的狀態最優估計法[18],根據不斷采集的新數據對上一步識別結果進行修正,廣泛應用于結構的各項參數識別中,它能夠對系統下一步的走向做出有根據的預測,可用于在線識別。黃進鵬等[19]利用擴展卡爾曼濾波法(extended Kalman filter,EKF)對車橋耦合作用下的橋梁結構進行分析,成功識別出橋梁結構的損傷,并有一定的抗噪能力。Liu等[20]提出了一種基于未知激勵的EKF方法的遞歸估計,在識別結構系統參數的同時,成功對未知外部激勵進行了識別。然而,基于KF方法的損傷識別是一個典型的反問題,在處理大型結構時,狀態向量的維數過大往往會對KF方法的計算效率產生不利的影響,也可能會導致KF方法收斂困難、識別結果不準確[21],因此考慮利用POD方法對結構進行降階分析,可以很好地解決KF方法在處理大型結構時收斂困難的問題。

1 基于POD的動態模型降階

本文提出的基于動態響應數據驅動的模型降階主要是根據獲取的系統動力學響應構造快照矩陣,選擇一組最能夠表征系統特征的標準正交基,以獲取結構的主要特征。

當結構受到外部載荷的影響時,其動力響應可以用下面的運動方程來描述

(1)

在利用POD方法進行模型降階時,我們考慮位移響應向量X∈Rm,R為實數的集合,m為位移響應向量X的維數。對于線性系統,位移響應可以寫成如下的變換形式

(2)

式中:yi為模態坐標下的位移響應,排列在列向量y中;{φi},i=1,…,m為一組標準正交基,稱為POMs;φ為坐標轉換矩陣,由所選取的正交基組成,表示為

φ=[φ1φ2…φm]

(3)

通常,我們僅需要保留前幾階POMs,就可以跟蹤系統的演化,因此定義了式(2)的簡化表示

(4)

式中:考慮l

首先,需要為向量X提供一個子空間,計算n個時刻的位移u(k)=u(tk)(k=1,…,n),并將其收集到快照矩陣U中

U=[u(1)u(2)…u(n)]

(5)

對快照矩陣U進行奇異值分解,可以得到

U=L·S·RT

(6)

(7)

當滿足p≥99%時,即可認為前l階POMs可以表示原始系統的大部分特征。以式(7)作為準則,選取此時的l值作為保留的特征模態階數,將快照矩陣奇異值分解后得到的左奇異矩陣的前l階列向量組成降階矩陣。

在研究過程中,發現較大的奇異值所對應的特征向量具有相對較多的能量,由于本文所研究的對象為線彈性模型,而線性系統的狀態可以用單個POM來表征,因此在對線彈性模型進行損傷識別分析時,僅需要選擇最大奇異值Sii所對應的POM即可。求得結構在無損條件下模型的一階POM與結構發生損傷時模型的一階POM,并繪制圖像,圖像中損傷結構一階POM發生突變的位置即為結構中損傷發生的位置,同時隨著損傷程度的增大,其圖像在損傷位置處的突變值也越大。

對于線性結構,通過一階POM便可以對結構的損傷進行識別。然而基于SVD的POD方法會對損傷識別的實時性產生影響,因此本文的研究重點是開發一種用于在線實時更新結構POMs的方法,從而對結構前期的微小損傷以及不同程度的損傷進行快速識別。

2 基于POD-KF的結構損傷識別方法

KF方法是一種基于系統狀態方程,通過系統輸入和輸出的觀測數據,以最小均方誤差為最佳估計準則,對系統的最優狀態進行估計的算法,廣泛用于結構各項參數的快速識別。

KF方法的核心思想是利用前一時刻的估計值和現在時刻的測量值不斷對狀態向量進行更新,來獲取狀態向量的最優估計值。在計算時首先定義系統方程式(8)和觀測方程式(9)

xk=Ak-1xk-1+Bk-1uk-1+wk-1

(8)

yk=Hkxk+vk

(9)

(10)

式中:xk為系統在k時刻的狀態向量;uk為已知的系統輸入值;Ak為系統的狀態轉移矩陣;Bk為已知輸入的位置矩陣;wk為系統噪聲;yk為系統在k時刻的觀測值;Hk為系統觀測轉移矩陣;vk為系統的觀測噪聲。

之后便可以根據上述系統方程和觀測方程得到5個KF方法的基本方程。首先是狀態預測過程,即根據第k時刻狀態對k+1時刻的狀態進行預測

(11)

(12)

接下來進入更新階段,定義實際的測量值與通過觀測方程得到的估計值之間的誤差為

(13)

定義此刻測量值與預測值之間誤差的協方差矩陣,并寫出更新矯正后的狀態估計方程

(14)

(15)

式中,R為觀測噪聲的協方差矩陣。則后驗誤差協方差矩陣可以表示為

(16)

式中,卡爾曼增益Kk作為最優加權值,其作用是使后驗誤差協方差矩陣取值最小,通過計算可以表示為

(17)

通過將以上公式循環迭代求解,就可以通過給定的初值x0和初始誤差協方差P0,根據已知的觀測值進行遞推,對系統的狀態進行最優估計。

本文考慮利用KF方法對基于POD的損傷識別方法進行改進,當出現損傷時,直接利用KF方法對所保留的投影子空間(即保留的POMs)進行在線更新,以達到實時定位損傷的目的。

為了考慮降階模型子空間隨時間的變化,將式(4)所保留的POMs以向量的形式排列如下

(18)

式中,Πl,k表示的POMs可以隨時間改變,考慮到因損傷導致的系統演化,子空間隨時間的演化可以表示為

(19)

與傳統KF方法一樣,式(19)需要用一個觀測方程來補充,并根據測量值對子空間進行更新,設置測量值

(20)

首先進入預測階段

(21)

(22)

對卡爾曼增益進行計算

(23)

對子空間進行更新并對其濾波誤差協方差矩陣進行求解

(24)

(25)

通過對上述公式進行迭代求解,便可以實時對結構的POMs進行在線更新,以達到識別損傷的目的。

3 一階POM的損傷敏感性分析

本章首先對POM的損傷敏感性進行分析,考慮了一個受地震加速度激勵的六層剪切型框架,如圖1所示。各層質量為:m1～m6=400 kg;層間剛度分別為:k1=120 kN/m,k2=k3=k4=100 kN/m,k5=k6=70 kN/m;層間阻尼分別為:c1=6 kNs/m,c2=c3=c4=5 kNs/m,c5=c6=3.5 kNs/m。地震動選擇EL-Centro波,作用于框架底部,峰值加速度為PGA=0.5g,采樣時間間隔為0.02 s,所收集位移響應的總時間為55 s。

圖1 六層剪切型框架計算模型Fig.1 Calculation model of six-story shear frame

首先對該結構的動響應進行求解,利用Newmark-β法計算框架結構在地震作用下的位移和速度響應,并選取其位移響應作為觀測值,同時為了更貼近實際情況,在位移響應中加入了5%的高斯噪聲。之后對該剪切型框架進行POD分析:利用獲取的位移響應構造快照矩陣,對其進行奇異值分解得到結構的POMs,如圖2所示,從圖2可以看出POMs所包含的能量隨著階數的增大而急劇減小,同時從表1的結果中發現,一階POM的能量占比已經達到了99%,可以用于表示結構的主要特征并保證降階模型有足夠的精度。將全階模型的結構響應與一階POM構成的降階模型的結構響應進行對比,如圖3和圖4所示,可以看出降階模型和全階模型的速度、位移響應大致吻合,再次證明了一階POM可以揭示出該系統最本質的特征。

表1 未損傷結構和損傷結構本征值及能量占比

圖2 未損傷結構和損傷結構的本征值(POVs)Fig.2 POVs of undamaged and damaged structures

圖3 全階模型與降階模型三層、四層間位移響應Fig.3 Displacement response at the inter-story between the third and fourth floors of the full-order model and the reduced-order model

圖4 全階模型與降階模型三層、四層間速度響應Fig.4 Velocity response at the inter-story between the third and fourth floors of the full-order model and the reduced-order model

圖5 位于第三層和第四層間不同程度損傷結構的一階POMFig.5 The first-order POM of structures with different degrees of damage located at the inter-story between the third and fourth floors

同樣,將損傷施加到其余樓層位置,以驗證該方法對不同位置損傷進行識別的可行性。圖6為1～6層分別發生50%損傷的一階POM圖,從圖6可以看出:對于不同的損傷位置,一階POM對損傷均有著較高的敏感性;當損傷程度相同的情況下,框架底部及其相鄰層發生損傷時其識別效果最明顯,這是由于其距載荷施加的位置較近,其他樓層識別效果雖然相對較弱,但同樣可以對損傷進行識別。

圖6 d=50%,未損傷結構與不同位置損傷結構的一階POMFig.6 d=50%,the first-order POM for the undamaged structure, and for the damaged located at a varying inter-story level

綜上所述,對于剪切型框架結構,傳統POD方法得到的一階POM在處理損傷問題時有著較高的敏感性,而本研究的目的是開發一種能夠實時跟蹤損傷演化的方法,接下來將利用POD-KF方法對數值模型進行損傷識別分析,驗證其有效性。

4 POD-KF方法的數值驗證

本章利用所提出的POD-KF方法對結構損傷進行識別,將加入高斯噪聲的位移響應作為KF方法的觀測值進行識別。在進行濾波前,需要對算法的協方差矩陣進行設置:初始協方差矩陣Ξ0=1×10-12I;系統噪聲協方差矩陣WΠ=1×10-10I;觀測噪聲協方差矩陣Vy=1×10-10I。上述參數被稱作卡爾曼濾波調節參數,通常需要在進行濾波之前進行設置與調整,這些參數應不斷調整到最優值,以允許在濾波的同時對變化的參數做出快速反應[22]。本文所提出的方法能夠利用結構響應對結構的時變參數進行追蹤,因此調試的根本宗旨是能盡量利用測量方程,也就是使卡爾曼增益在迭代過程中始終保持較大值。

模型選用第3章的剪切型框架模型,而損傷的設置與上章有所不同,本文所分析的是時變損傷(即在分析的某一時刻出現一定程度的損傷),首先考慮單損傷工況,損傷出現時間設置在第8 s,損傷程度為50%,損傷位置位于第三層和第四層之間;利用POD-KF方法對其進分析,識別結果如圖7所示。

圖7 d=50%,POD-KF法一階POM時程圖Fig.7 d=50%, the first-order POM time-history diagram of the POD-KF method

圖7的結果顯示,當損傷出現在第8 s時,結構的一階POM圖像開始發生變化,但是在第10 s左右便可以很快收斂到穩定值,并且在損傷位置處發生突變。為了更加直觀地看到損傷識別結果,圖8對比了POD-KF方法得到的損傷與未損傷結構的一階POM時程變化,灰色的球型線表示未損傷結構,而黑色的方形線表示損傷結構,可以看出當損傷發生時,方形線的時程圖像會在損傷位置出現明顯的突變。

圖8 d=50%,POD-KF方法損傷與未損傷結構一階POMFig.8 d=50%, the first-order POM of damaged and undamaged structures by POD-KF method

正如第2章提出的方法,這里對損傷結構的POMs的求解,并不需要像傳統POD方法一樣,需要對損傷結構模型進行重新訓練獲得,而是利用KF方法不斷進行更新,以更好地與實際測量的結構響應相匹配,大大縮短計算時間。由于當損傷出現時,沒有經過重新訓練,因此子空間更新的最終結果會與傳統POD方法有一定的差異,如圖9所示。

圖9 d=50%,POD-KF方法一階POM與POD方法對比Fig.9 d=50%, comparison of the first-order POM of POD-KF method and POD method

圖9將傳統POD方法得到的結果與通過KF方法更新后達到穩定的POM進行了對比,并在表2中列出了兩者間的誤差,誤差較小,說明利用KF方法更新的POM具有較高的精度。同時為了驗證該方法對不同損傷程度與不同損傷位置的可行性,對利用一階POM構成的降階模型進行分析,考慮降階剛度κ11的變化,從圖10的結果中我們可以看出,當損傷位置從1層變到6層,損傷程度從0.05%增加到50%時,降階剛度κ11呈現線性的變化,這為KF方法的應用提供了前提。

表2 POD-KF方法與傳統POD方法對比

圖10 降階剛度矩陣對不同位置和損傷程度的變化Fig.10 Variation of reduced order stiffness matrix on different positions and damage degrees

而對于大型結構,大損傷對系統參數會產生較為明顯的影響,小損傷則影響較小,再加上噪聲的存在,導致目前大多數損傷識別方法無法對小損傷進行識別。然而如果不對小損傷加以重視,其會在極短的時間內發展為大損傷,在人們還未察覺的情況下對整個結構進行破壞,因此探究本文提出的方法對小損傷的敏感性,無論是在理論和實際工程中,都是一個值得研究的內容。

仍采用第3章提出的框架模型,設置損傷位于第三層和第四層之間,損傷出現時間設置在第8 s,為了驗證該方法對小程度損傷的敏感性,損傷程度為5%,利用POD-KF方法對其進行損傷識別分析,識別結果如圖11所示。

圖11 d=5%,POD-KF方法損傷與未損傷結構一階POMFig.11 d=5%, the first-order POM of damaged and undamaged structures by POD-KF method

結果顯示,當損傷程度較小時,本文提出的POD-KF的方法仍可以較為高效地對其進行識別,從圖像中可以看出,利用KF方法進行更新后,仍然能夠很快收斂到穩定值,并將達到穩定后損傷結構與未損傷結構的一階POM進行對比,結果如圖12所示,可以看出當出現小程度損傷時,一階POM的曲線在損傷處會發生突變,具有較好的識別效果,綜上所述,本文采用的POD-KF方法即使是處理較小程度的損傷時仍具有較好的識別效果,證實了該方法識別結構早期微小損傷在理論上的可能性。

圖12 d=5%,損傷結構與未損傷結構的一階POM對比Fig.12 d=5%, comparison of the first-order POM between damaged and undamaged structures

然而在實際工程中,結構損傷往往不是單獨存在的,而是多個不同的損傷同時存在,這對結構的威脅更大,因此接下來將考慮多工況損傷,損傷出現時間同樣設置在第8 s,損傷程度為50%,但是損傷位置設置為兩位置損傷,分別在第二層和第五層。利用POD-KF方法對其進分析,識別結果如圖13所示。

圖13 d=50%,多損傷工況下POD-KF方法損傷與未損傷結構一階POMFig.13 d=50%, the first-order POM of damaged and undamaged structures by POD-KF method in multi-stage damage

在圖13中對比了POD-KF方法得到的損傷與未損傷結構的一階POM時程變化,灰色的球型線表示未損傷結構,黑色的方形線表示損傷結構,當損傷出現在第8 s時,與單損傷工況結果一致,結構的一階POM圖像開始發生變化,同樣在第10 s左右便可以收斂到穩定值,并且在損傷位置處發生突變。可以看出當損傷發生時,方形線的時程圖像會在損傷位置出現明顯的突變。圖14比較了POD-KF方法得到的一階POM終值與未損傷結構的一階POM,與圖6中的結果相比較,發現當損傷出現在某兩層時,其圖像是其單層損傷圖像的疊加:設置第二層、第五層出現損傷,其圖像可以看作是圖6(b)和圖6(e)的疊加。該結果表明,所發展的POD-KF方法無論是處理單損傷工況還是多損傷工況,均可以對其進行快速、準確地識別。

圖14 d=50%,多損傷工況下損傷結構與未損傷結構一階POM對比Fig.14 d=50%, comparison of the first-order POM between undamaged and damaged structures in multi-stage damage

5 抗噪性分析

由于實際工程中更容易受到外部環境、監測設備等諸多因素的影響,所得到的響應數據往往需要考慮不同水平的噪聲,為了使該方法更好地應用于工程實際,本節在第4章單損傷工況的基礎上,添加了更高水平的噪聲,以驗證該方法在實際工程中的適用性。

選擇的工況與第4章單損傷工況相同,但是在位移響應中加入了15%的高斯噪聲,利用POD-KF方法對其進行分析,識別結果如圖15所示。

圖15 d=50%,高水平噪聲下POD-KF方法損傷與未損傷結構一階POMFig.15 d=50%, the first-order POM of damaged and undamaged structures by POD-KF method in a high level of noise

圖15對比了噪聲水平增高后POD-KF方法得到的損傷與未損傷結構的一階POM時程變化,可以看出圖像前4 s的結果相比于低水平的噪聲,黑色的方形線波動較大,這說明在該算法的初始階段,噪聲影響較大,但隨著POD-KF方法的推進,圖像逐漸收斂到穩定值,之后在第11 s左右,圖像在損傷位置處發生變化,并逐漸穩定,成功識別出結構的損傷。圖16表示了不同水平噪聲下損傷識別結果誤差,可以看出雖然有一定的差異,但是基本重合,仍然滿足損傷識別的要求。該結果說明,得益于卡爾曼濾波顯著的降噪能力,POD-KF方法在處理高水平噪聲時仍有著出色的識別能力,有著較強的抗噪能力,并且隨著時間的推進,穩定性會持續增強,這為該方法應用到實際工程中提供了前提。

圖16 d=50%,不同水平噪聲下損傷結構與未損傷結構的一階POM對比Fig.16 d=50%, comparison of the first-order POM between damaged and undamaged structures in different levels of noise

6 結 論

本研究充分利用POD和KF的優點,并將兩者有機結合,發展了一種基于子空間更新的、數據驅動的POD-KF損傷識別方法,并以六層剪切型框架為例,通過數值試驗,驗證了該方法的有效性。主要得到以下結論:

(1)利用KF方法對子空間進行更新,可以動態跟蹤反應結構損傷演化的一階POM,無論是單損傷工況還是多損傷工況,均可以實時、快速地對結構的損傷進行識別,而不需要高精度的模型,結果表明該方法可以在保持較高精度的前提下,大大減少計算量,并可以有效地抑制噪聲,即使在處理結構早期小程度的損傷時,仍具有較高的精度。

(2)KF方法在識別結構參數時需要同時對外載荷進行識別,具有很大的局限性;而本文所采用的POD-KF方法,僅利用不斷獲取的位移響應數據便可對結構的一階POM進行更新,無需對外載荷進行求解,使其能夠對受復雜載荷影響的結構進行快速識別,具有一定的創新性。

綜上所述,利用POD-KF方法對結構的一階POM進行更新具有實時識別損傷的能力,該方法可為基于響應數據驅動的動力學降階以及結構損傷識別領域提供新的理論和技術支撐。在接下來的研究中,將考慮更復雜的模型,并針對實際工程中的問題進行分析,在滿足損傷識別精度的同時,對傳感器的數量以及位置進行優化。