基于移動窗卡爾曼濾波算法的結構響應重構

2021-11-17 12:25:42張笑華吳志彪吳圣斌黃梅萍

振動與沖擊 2021年21期

張笑華，吳志彪，吳圣斌，黃梅萍

(福州大學土木工程學院，福州 350108)

結構健康監測的首要關鍵工作就是通過傳感器系統正確地測量結構的動力響應。然而由于經濟條件、儀器設備和現場環境等因素的制約，安裝在結構上的傳感器數目始終有限，許多位置并未安裝傳感器，從而無法為結構健康監測提供完整有效的信息。因此，利用有限位置上的測量響應來重構未測量位置的結構響應具有重要意義。

目前，已有的響應重構方法主要是基于傳遞的概念，即通過在時域、頻域或者小波域里構造傳遞矩陣和傳遞函數來實現響應重構[1-3]。這些方法主要考慮了單種響應的重構。已有的研究結果表明：用多類型的測量信息來識別結構參數，其精度高于用單一信息識別的結果，更有利于結構分析、監測和安全評估[4]。因此，Zhang等[5]基于模態理論研究了位移和應變響應重構的方法。Zhu等[6]擴展了該方法，用經典卡爾曼濾波(Kalman filter, KF)算法，考慮傳感器位置的優化，重構未測點加速度、位移和應變三種響應。

KF算法是一種以最小均方誤差為最佳準則，利用測量響應對系統狀態向量進行最優估計的迭代算法。因其優越的性能，該算法也被廣泛應用于土木工程領域中，包括線性和非線性結構參數識別[7-8]，結構損傷識別[9]，有限元模型修正[10]等。近年來，也有不少學者將KF算法應用于結構響應的重構。Papadimitriou等[11]利用KF算法估計應力功力譜，用于識別因疲勞而引起的結構損傷累積。Jo等[12]融合加速度和應變測量響應，利用KF算法估計結構的位移響應。Zhang等[13]基于擴展KF算法，通過有限的傳感器測量響應，首先估計外部激勵，然后重構未測量位置的結構響應。萬志敏[14]基于貝葉斯理論，研究了類卡爾曼濾波算法的外激勵識別和結構響應重構方法。董康立等[15]分別將逐步消去法和螢火蟲算法與類卡爾曼濾波算法相結合，提出了面向多類型傳感器優化布置的結構響應重構方法。

目前大多數基于KF算法的結構響應重構方法，需預先假定測量與過程噪聲方差已知且為常數。然而，大部分情況下，它們是未知且是時變的。不精確的噪聲方差可能導致較大誤差的狀態向量估計，甚至可能導致濾波發散。Fuad等[16-17]研究了過程和測量噪聲方差未知情況下的KF算法。Lai等[18]利用移動窗技術，通過測量響應和系統實際輸出之間的殘差，實時估計過程噪聲方差和測量噪聲方差，并用于KF算法以識別結構剛度的變化。Zhang等[19]改進了Lai等的測量噪聲估計方法，研究了測量噪聲方差未知下的結構響應重構，但該方法需先假定過程噪聲方差已知。

基于以上背景，考慮測量和過程噪聲方差均未知的情況，研究有限測點下基于移動窗卡爾曼濾波(moving-window Kalman filter, MWKF)算法的結構響應重構。利用有限測點的位移和應變測量響應，先采用移動窗技術估計結構的測量與過程噪聲方差，然后基于KF算法重構結構中未安裝傳感器位置處位移與應變響應，最后以一個平面單跨框架結構為例進行數值模擬和試驗分析，驗證所提方法的有效性和可行性。

1 MWKF算法

線性系統的結構動力學方程可寫成連續的狀態空間方程。由于實測數據都是離散的，且實際工程的信號測量中，總是存在系統的不確定性，若將這些不確定性分為測量噪聲和過程噪聲，則連續的狀態空間方程經過離散后，可得到下述的離散時間隨機狀態空間方程

zk+1=Azk+Buk+wk

(1)

yk=Czk+Duk+vk

(2)

其中：

(3)

(4)

(5)

式中：下標k表示第k時刻；z為離散狀態向量；q是模態坐標向量，y是離散觀測向量；A、B分別是離散狀態矩陣和外部激勵的輸入矩陣；C，D分別代表測量輸出矩陣和傳遞矩陣。Δt是采樣時間間隔；I是單位矩陣；u是外部激勵向量；w表示過程噪聲向量；v表示測量噪聲向量；w和v假定為均值為零且互不相關的高斯白噪聲，它們的協方差矩陣分別為Q和R。Φh是質量歸一化后的位移振型矩陣；Bu表示激勵位置矩陣；ξ是模態阻尼系數矩陣；ω0是模態頻率矩陣。

KF算法是無偏遞歸算法，可最優估計未知狀態向量。經典KF算法為

(6)

(7)

先驗誤差方差Pk|k-1

Pk|k-1=APk-1AT+Q

(8)

狀態向量估計均方差Pk

Pk=[I-KkC]Pk|k-1

(9)

卡爾曼濾波增益Kk

Kk=Pk|k-1CT[CPk|k-1CT+R]-1

(10)

算法中的R和Q通常需要通過大量實際實驗的統計特性或在算法實際操作之前憑經驗評估。在土木工程領域，它們常被假定為常數且為恒定。然而，事實上它們經常是時變的，尤其是在結構長期的運營服役期間。因此，本文基于移動窗技術，首先估計噪聲方差陣R和Q。

定義測量噪聲σk為實際響應與測量響應之間的差值。由于實際響應難以獲得，這里用測量響應的加權平均來代替實際響應，具體計算如下

yk，smooth=λ1yk-1，m+λ2yk，m+λ3yk+1，m

(11)

式中：λ1,λ2,λ3為加權系數；下標‘m’表示測量的。用yk，smooth代替實際響應計算測量噪聲

σk=yk，m-yk，smooth

(12)

(13)

其中

(14)

式中，Nr表示移動窗長度。

(15)

其中γi為加權系數，計算如下

(16)

式中，tr(·)表示計算矩陣的跡。測量噪聲方差的估計僅需要通過測量數據來獲得而不涉及任何結構特性。

測量殘值即新息(innovation)ek定義為測量估計值與實際測量值之差，它表示存在于測量響應中而濾波器無法預測的那部分信息，即

(17)

其對應的自相關系數Hτ定義如下

(18)

式中，Nq為移動窗長度。測量殘差的自相關系Hτ與先驗誤差方差之間有如下關系

當τ=0時

Hτ=CPk|k-1CT+R

(19)

當τ>0時

Hτ=

CA(I-KkC)τ-1A[Pk|k-1CT-Kk(CPk|k-1CT+R)]

(20)

(21)

式中：a為控制收斂速率的系數矩陣；H0為測量殘值的方差矩陣。

2 結構位移和應變響應重構

考慮應變與位移兩種響應，觀測向量為

(22)

因此，式(2)中的C和D矩陣可寫為

(23)

式中:Ψh=BsdΦh為應變模態矩陣;Bsd是應變-位移轉換矩陣，與有限元劃分的單元類型形函數有關。這里將輸出矩陣C分為Cm和Ce兩部分，分別對應于測量的位置和需響應重構的位置。

由于同時考慮了應變和位移兩種響應，二者在數值上有量級差異，引起對應輸出矩陣C的高度病態，導致卡爾曼增益矩陣計算結果不準確。可采用應變和位移對應的測量噪聲標準差來規則化輸出矩陣Cm

(24)

其中，

(25)

(26)

未測點的響應時程按下式計算

(27)

實際上，并非所有的振型對結構的整體動力響應具有同等重要的貢獻，因此在使用基于MWKF算法重構未測點結構響應時，可根據實際情況僅考慮對動態響應貢獻最大的前幾模態振型，忽略高階模態的影響。

3 數值模擬分析

采用三層單跨框架為數值算例驗證所提方法的有效性和可靠性。該框架每層高度為0.6 m，跨度為0.6 m，梁與柱的橫截面尺寸分別為50 mm×8.8 mm和50 mm×5 mm，彈性模量為206 GPa，密度為7 850 kg/m3，泊松比為0.3。同時也建立了相應的有限元模型，模型共有53個節點和54個單元，其中節點1和節點20為固定端，模型共有153個自由度，如圖1所示。根據模態振型對減小響應重構誤差的貢獻情況，選取了前6階模態為目標模態進行未測點的響應重構[20]。前6階的模態頻率分別為3.78 Hz、11.21 Hz、17.17 Hz、46.06 Hz、54.48 Hz和66.17 Hz。根據文獻[21]的方法得到傳感器位置，如圖1所示。

外荷載以隨機激勵為例，利用狀態空間方程并考慮過程噪聲的影響計算結構動力時程響應。該隨機激勵的頻帶寬度為0.5～50 Hz，最大幅值為100 N，施加于節點38-x向。假設應變對應的過程噪聲方差陣為Qs=1×10-15·I，位移對應的過程噪聲方差陣Qd=1×100·I。提取對應有安裝傳感器位置處的動力響應并人為加入3%白噪聲后作為測量響應。

圖1 框架結構傳感器布置Fig.1 Sensor placement of the frame model

圖2是測量噪聲方差時程曲線。圖2(a)對應于節點19處位移測量響應，圖2(b)對應于單元7處應變測量響應。由圖2可知，估計的測量噪聲隨著時間的推移在整個過程中不斷變化且會在某一值上下波動。由于用測量響應的加權平均值代替了實際響應，因此實際的測量噪聲方差和估計的測量噪聲方差均值之間稍有差異。圖3為節點19處位移響應對應的過程噪聲方差時程與單元7處的應變響應對應的過程噪聲方差時程。由圖3可知，估計的過程噪聲方差也隨著時間的變化而變化，應變對應的過程噪聲方差趨于1×10-15，位移對應的過程噪聲方差趨于1×100。不論是應變還是位移，過程噪聲方差都趨于預設值。

為了評估測量和過程噪聲未知條件下，基于KF算法重構響應的精度，定義了相對百分比誤差，計算如下

(a) 節點19

(b) 單元7圖2 測量噪聲方差Fig.2 Measurement noise covariance

(a) 節點19

(b) 單元7圖3 過程噪聲方差Fig.3 Process noise covariance

(28)

式中:std指標準差;下標j表示的是第j個單元或節點。

圖4顯示了使用KF算法(R與Q均已知)與使用MWKF算法(R與Q均未知)重構未測點響應得到的相對百分比誤差對比。其中，R與Q均已知時，位移響應重構的相對百分比誤差大部分都在2%以內，平均相對百分比誤差0.81%；應變響應重構的相對百分比誤差大部分都在8%以內，平均相對百分比誤差為2.08%；R與Q均未知時，位移響應重構的相對百分比誤差在3%以內，平均相對百分比誤差為1.49%；應變重構的相對百分比誤差都在10%以內，平均相對百分比誤差為2.93%。兩種情況得到的響應重構結果精度相差不大。圖5為使用MWKF算法得到的單元4應變重構響應時程、節點19-x向位移重構響應時程與相應計算響應時程的對比。由圖5可知，重構響應時程與有限元軟件計算得到的響應時程吻合良好。數值算例的結果表明，在測量和過程噪聲方差未知的情況下，基于MWKF算法，可對測量和過程噪聲方差進行實時估計，并用于結構未測點的響應重構，能取得較好的重構精度。

(a) 應變

(b) 位移圖4 相對百分誤差Fig.4 Relative percentage errors

(a) 單元4應變響應

(b) 節點19-x向位移響應圖5 計算響應與重構響應時程對比Fig.5 Comparison of reconstructed and real time history responses

4 實驗室框架試驗驗證

根據數值算例中框架結構的幾何尺寸和物理參數，在實驗室制作了相應的單跨三層鋼框架，模型如圖6所示。框架的材料選用Q235鋼，框架底座通過兩個L型角鋼錨固在固定于地面的工字型鋼梁底座上。本次試驗中應變片的型號為BX120-5AA；位移傳感器采樣sensorpart FT 25 RA型激光位移計和KEYENCE IL-300型激光位移傳感器；數據采集系統的型號為JM5959采集系統。

圖6 實驗室框架結構Fig.6 Experimental model

文章提出的基于MWKF算法的結構響應重構方法是基于有限元模型計算的振型向量。因此，首先須對有限元模型進行修正，確保有限元模型與實際結構的動力特性基本一致。為此，本次試驗首先進行了結構模態試驗，以獲得結構的模態頻率和振型。表1列出了前6階的頻率值和振型的MAC(modal assurance criterion)值。從表中可以看到，頻率的測量值與計算值對比，最大誤差為1.59%，其余誤差基本都低于或者接近1%。前6階模態振型的MAC值均與1.0非常接近。這些對比結果表明有限元模型與實驗室模型的動力特性有很高的相似度，因此可以用有限元模型來預測實際結構的動態特性，不需要進行有限元模型修正。

表1 前6階頻率值及MAC值

用于響應重構試驗的傳感器位置如圖7所示。除了12個測量點外，在節點11和26額外安裝了2個激光位移計，在單元4、18、28、39和47增加貼了5個應變片，用于檢驗重構響應的精度。本文提出的響應重構方法與外荷載的類型無關。本試驗施加的外荷載以錘擊激勵為例，施加在節點35-x位置上。為了驗證不同的采樣頻率對重構方法的影響，本試驗考慮了兩種情況下的響應重構：

工況1：采樣頻率200 Hz；

工況2：采樣頻率500 Hz。

響應重構時考慮前6階模態，加權因子、移動窗口大小與數值模擬分析中的取值相同。采集的動力響應首先用70 Hz的截止頻率進行濾波，去除高頻響應的影響，然后結合測量的位移和應變響應，對測量和過程噪聲進行實時估計，并用于結構未測點的位移和應變響應重構。

圖7 框架結構傳感器布置Fig.7 Sensor placement of frame structure

將測量響應替代式(28)中的計算響應y，計算相對百分比誤差。表2列出了兩個工況的相對百分比誤差。工況1中，應變和位移響應重構最大相對百分比誤差分別為14.72%和9.75%；工況2中，應變和位移響應重構最大相對百分比誤差分別為15.33%和11.15%。兩個工況的相對百分比誤差大部分都低于15%，采樣頻率不論是為500 Hz還是200 Hz，二者重構結果精度相差不大。圖8和圖9分別為工況1估計的測量噪聲方差時程和估計的過程噪聲方差時程，對應于單元7的應變和節點32的位移。圖10和圖11分別為工況1在時域和頻域內單元4的應變重構響應、節點11-x位移重構與測量響應之間的對比，無論是時域內還是頻域內，重構響應與測量響應之間的誤差都很小，吻合良好。