統計容積卡爾曼濾波器的混合試驗模型更新方法

王 濤， 李 勐， 孟麗巖， 許國山， 王 貞

(1.黑龍江科技大學 建筑工程學院，哈爾濱 150022；2.哈爾濱工業大學 土木工程學院，哈爾濱 150090；3.武漢理工大學 土木與建筑學院，武漢 430070)

混合試驗方法是一種結合試驗子結構試驗加載及數值模擬的技術，是檢驗大型復雜結構動力災變性能的有效試驗手段。為提高混合試驗數值模型精度，模型更新混合試驗方法應運而生[1]，該方法將模型更新技術與混合試驗方法相結合，利用試驗子結構觀測數據在線識別試驗子結構模型參數，實時更新數值子結構模型參數，以提高數值模型精度。

迄今為止，模型更新混合試驗的研究已獲得較大進展[2-3]。Zhang等[4]采用最小二乘法識別了雙折線模型及Bouc-Wen模型參數，并將其應用于混合試驗中。Hashemi等[5]以OpenFresco為基礎，完成了單層單跨框架結構的模型更新混合試驗，采用無跡卡爾曼濾波器(unscented Kalman filter,UKF)識別Bouc-Wen模型參數，并在線更新了數值子結構模型參數。王濤等[6]提出UKF模型更新混合試驗方法，在線識別單自由度Bouc-Wen模型的參數，并通過對彈簧試件的單自由度體系的實際試驗驗證了此試驗方法的有效性。梅竹等[7-8]使用UKF算法對鋼筋混凝土結構材料本構模型參數進行了在線參數識別。陳再現等[9]提出基于均勻設計的模型更新混合仿真方法，在試驗子結構的基礎上，通過均勻設計的手段構造驗算子結構樣本空間，得到與試驗子結構最為吻合的數值模型進行模型參數更新，有效提高了識別效率。Atasaratnam等[10]提出容積卡爾曼濾波器(cubature Kalman filter，CKF)算法，CKF基于數值積分理論，其核心步驟是根據三階球面徑向容積規則通過系統的初始矩陣、初始狀態和初始點集確定2n個容積點[11]，并使這些容積點接近具有高斯噪聲非線性系統的狀態和協方差[12]。與UKF相比，CKF避免線性化的原理不同[13]，CKF采用三階球面徑向容積準則逼近高斯積分，具有較高的精度和穩定性，其容積點和權重僅由狀態維數決定，可以提前計算和存儲[14-15]，是建立在貝葉斯估計理論框架下的次優濾波，其結果滿足正態分布下的最小方差理論。王濤等[16]采用CKF針對Bouc-Wen模型在線參數識別并與UKF進行對比，研究表明，CKF在運行效率上，仿真耗時少，識別結果更接近真實值，并能有效降低模型參數識別過程中的波動幅度，提高強非線性模型在線參數識別精度。然而，在模型更新混合試驗中，CKF算法受制于算法初始參數的選擇，初始參數選擇不當，往往會降低算法精度，有時導致混合試驗失敗，大大限制了模型更新混合試驗的工程應用。

鑒于此，本文以CKF算法為例，提出基于統計CKF混合試驗模型更新方法，解決由模型更新算法初始參數選擇不當對模型參數識別結果的影響。以自復位摩擦耗能支撐模型為加載對象，應用該方法對模型進行參數識別，分析在不同參數條件下統計CKF算法的識別精度；以兩層帶有自復位摩擦耗能支撐框架結構進行模型更新混合試驗數值仿真，檢驗該方法的魯棒性。

1 統計CKF模型更新算法

應用CKF進行參數識別時，需人為假定算法初始參數，初始參數的選擇對參數識別結果有很大影響。目前，沒有方法確定合理初始參數值，為此，筆者提出統計CKF模型更新算法，以弱化模型初始參數選擇不當對識別結果的影響，提高算法對初始參數的魯棒性。

1.1 統計CKF模型更新算法步驟

首先，建立如下系統的狀態方程和觀測方程

Xk=f(Xk-1,uk-1)+Vk-1

(1)

Yk=h(Xk,uk)+Wk

(2)

式中：u為輸入向量；V為過程噪聲向量，假定過程噪聲向量是均值為0的高斯白噪聲，其協方差矩陣為R；W為觀測噪聲向量，假定觀測噪聲向量是均值為0的高斯白噪聲，其協方差矩陣為Q；X為系統狀態向量；Y為系統觀測向量；k為算法運行步數，k=1,2,…N。

(3)

(4)

ζ為容積點集，其表達式為

(5)

[e]i=

(6)

式中，[e]i為第i個容積點。

(7)

(8)

(9)

(10)

重采樣

(11)

(12)

(13)

(14)

(15)

(16)

校正狀態量

(17)

式中,yk為第k步的試驗觀測值。

校正協方差陣

(18)

對算法運行M次的狀態量識別值進行統計，將統計后的狀態量識別值均值作為最終識別結果

(19)

圖1 統計CKF模型更新算法流程圖
Fig.1 Statistical CKF model update algorithm flow chart

2 自復位摩擦耗能支撐模型參數識別

2.1 經改進旗形模型表達式

為驗證基于統計CKF算法的有效性，以自復位摩擦耗能支撐為研究對象，采用經改進的旗形模型(modified flag-shaped model,MFS)[17]模擬自復位摩擦耗能支撐的滯回性能，對MFS模型進行模型參數在線識別。MFS模型由線性部分、雙線性彈性部分、具有可滑移彈塑性部分組成，自復位摩擦耗能支撐機理圖如圖2所示。初始階段剛度K為模型線性、雙線性、具有可滑移彈塑性三部分的剛度之和，即K=K1+K2+K3。當摩擦裝置啟動后，支撐剛度減小至預應力筋剛度K1，卸載初始階段的卸載剛度等于K1+K3，當力到達激活點b時，剛度變成K1+K2。

(a) 線性部分

(b) 雙線性部分

(c) 具有滑移區彈塑性部分

(d) 整體圖2 自復位摩擦耗能支撐機理圖Fig.2 Self-centering energy dissipation brace mechanics diagram

MFS模型恢復力F可表達為

F=K1x+K2R(x)+K3z

(20)

式中，R(x)為由雙線性彈性模型提供的恢復力，kN。

R(x)具體表達式可由式(21)給出

R(x)=x[1-H(x-b)-H(-x-b)]+

b[H(x-b)-H(x-b)]

(21)

式中：H(x)為海維賽德階躍函數；b為激活位移，mm；z為模型滯變位移，mm。

z微分表達式由式(22)給出

(22)

K1、K2、K3由式(23)～(25)給出

K1=AK

(23)

K2=(1-A)(1-Q)K

(24)

K3=(1-A)QK

(25)

式中，A為激活剛度比，其表達式為

A=(1-Q)b

(26)

Q為耗能率，是能量耗散的重要指標，其表達式為

Q=K3/(K2+K3)

(27)

本文設定MFS模型參數分別為：耗能率Q=0.9、激活位移b=1.4 mm、初始剛度K=300 kN/mm、激活剛度比A=0.9。

2.2 統計CKF算法初始參數

首先，假定狀態量為MFS模型滯變位移z、加載位移x、及模型參數b、K、Q，觀測量為MFS模型恢復力，并根據式(22)、(20)建立算法狀態方程及觀測方程，算法狀態方程和觀測方程如式(28)～(29)所示。

(28)

式中,x為加載位移。

(29)

在數值積分法得到的狀態方程的離散形式中，加入過程噪聲Vs，過程噪聲Vs和觀測噪聲Ws分別為均值為零的高斯白噪聲，過程噪聲Vs和觀測噪聲Ws的協方差矩陣分別為Rs、Qs，具體表達方式如式(30)～(31)所示。

Rs=10-25×I5

(30)

式中,I5為5×5的單位矩陣。

Qs=10-8×I1

(31)

(32)

式中，dt為數值積分步長。

(33)

(34)

作用于模型的加載方式選用低周期往復位移加載，加載方式如圖3所示。

圖3 加載制度Fig.3 Loading protocol

2.3 不同初始參數下統計CKF方法識別MFS模型參數

CKF算法識別模型參數的可靠程度往往受制于算法初始參數影響，針對這一問題，本節以MFS模型為例，研究模型使用不同初始參數條件下的模型參數識別結果，不同初始參數條件包括：不同初始協方差矩陣、初始參數預估值以及觀測噪聲，以驗證統計CKF模型更新算法的識別精度。

2.3.1 初始協方差矩陣

(a) 滯變位移識別值

(b) 加載位移識別值

(c) b識別值

(d)K識別值

(e) Q識別值

(f) 滯回曲線圖4 不同初始協方差矩陣的參數識別值結果Fig.4 Results of parameter identification values of different initial covariance matrices

本次識別過程中，CKF單次識別平均耗時4.01 s，使用統計CKF算法總耗時119.52 s。由圖4(a)～(e)可以看出，使用不同的初始協方差矩陣導致每次運行CKF算法識別參數的結果并不唯一，與模型真實參數值差異較大。從圖4(f)可以看出單次使用CKF算法得到的滯回曲線不能反映模型的真實性能，甚至沒有意義。采用統計CKF算法識別模型參數精度接近真實值，收斂性較好，參數識別效果得到改善，恢復力預測精度明顯提高。此外，通過耗能的角度分析，單次使用CKF算法得到的模型滯回曲線最大耗能值與真實耗能值相對誤差為86.26%，統計CKF算法得到的滯回曲線耗能值與真實耗能值相對誤差15.41%，耗能值相對誤差結果降低70.85%；單次使用CKF算法得到的模型滯回曲線最小耗能值與真實耗能值的相對誤差為10.12%，相比使用統計方法得到的耗能值相對誤差降低了5.29%。雖然單次使用CKF算法得到的耗能值相對誤差更低，但是考慮到CKF是一種隨機算法，每次應用算法得到的結果均存在隨機性，相比滯回曲線最大耗能值相對誤差，統計得到的模型滯回曲線有明顯改善，統計后的結果更具有可靠性，通過統計生成的滯回曲線與真實值差異較小，可以證實所提出方法的有效性。為了能定量評價統計CKF算法識別精度，定義均方根誤差(RMSE)為

(35)

式中：xref為模型參數真實值；xsim為模型參數識別值。

使用不同初始協方差矩陣，模型參數識別值的均方根誤差如圖5。

圖5 參數識別值均方根誤差圖Fig.5 Root mean square error diagram of parameter identification values

通過均方根誤差分析可以看出，相比CKF算法識別模型參數值最大均方根誤差，統計CKF算法得到的模型參數識別值均方根誤差分別降低了18.51%、73.78%、237.64%。初始誤差協方差矩陣不同，單次使用CKF算法得到的模型參數識別結果與真實值均方根誤差較高，統計CKF算法識別模型參數識別值結果均方根誤差更低，而且改善效果明顯，更能有效的反映結構加載后的真實性能，具有較強的魯棒性。不同統計次數下，得到的模型滯回曲線結果如圖6所示。

(a) 滯回曲線

(b) 相對誤差圖6 不同統計次數得到的模型滯回曲線及相對誤差圖Fig.6 The hysteretic curves and relative error graphs of different statistical times

為了檢驗統計次數對模型更新精度影響，分別對比統計10～50次得到的模型滯回曲線，如圖6(a)所示。由圖6(a)可以看出，隨著統計次數的增加，統計50次得到的模型滯回曲線較統計10～40次得到的模型滯回曲線更接近模型真實滯回形狀。由圖6(b)可以看出，從耗能的角度分析，統計50次得到的模型滯回曲線耗能值與真實耗能值相對誤差為11.56%，相比統計10～40次得到的滯回曲線耗能值相對誤差結果分別降低了32.67%、25.63%、3.85%、2.32%，可以發現當統計次數為30次時，改善效果最明顯。充分說明統計方法得到的參數識別結果可以提高算法魯棒性，優化算法隨機誤差。

2.3.2 初始參數預估值

本次識別過程中，CKF單次識別平均耗時3.95 s，使用統計CKF算法總耗時118.24 s。從圖7(a)～(e)可以看出，當初始參數預估值隨機時，每次運行CKF算法得到的參數識別結果與真實值具有很大差異，統計CKF算法識別模型參數結果可以較快地收斂至真實值。從耗能的角度看，圖7(f)單次使用CKF識別滯回曲線最大耗能值與真實耗能值最大相對誤差為132.37%，而統計過后的滯回曲線耗能值與真實耗能值相對誤差14.76%，耗能值相對誤差結果降低117.61%。單次使用CKF算法得到得模型滯回曲線最小耗能值與真實耗能值相對誤差為10.45%，相比使用統計CKF算法降低4.31%。使用不同初始參數預估值，模型參數識別值的均方根誤差如圖8所示。

(a) 滯變位移識別值

(b) 加載位移識別值

(c) b識別值

(d) K識別值

(e) Q識別值

(f) 滯回曲線圖7 不同初始參數預估值的參數識別值結果Fig.7 Results of parameter identification values estimated by different initial parameters

圖8 參數識別值均方根誤差圖Fig.8 Root mean square error diagram of parameter identification values

通過誤差分析可以得出，當使用不同初始參數預估值時，相比單次使用CKF算法識別模型參數最大均方根誤差，統計CKF算法識別模型參數值均方根誤差分別降低了69.68%、107.09%、244.91%，有效提高了算法識別精度。不同統計次數下，得到的模型滯回曲線結果如圖9所示。

為了檢驗統計次數對模型更新精度影響，分別對比統計10～50次模型滯回曲線，如圖9(a)所示。由圖9(a)可以看出，隨著統計次數的增加，統計50次得到的模型滯回曲線較統計10～40次得到的模型滯回曲線更接近模型真實滯回形狀。由圖9(b)可以看出，從耗能的角度分析，統計50次得到的模型滯回曲線耗能值與真實耗能值相對誤差為12.44%，相比統計10～40次得到的滯回曲線耗能值相對誤差結果分別降低了73.29%、28.55%、2.32%、1.77%。充分證明了統計方法得到的模型滯回曲線更接近真實值，可以更好地模擬結構的真實反映。

(a) 滯回曲線

(b) 相對誤差圖9 不同統計次數得到的模型滯回曲線及相對誤差圖Fig.9 The hysteretic curves and relative error graphs of different statistical times

2.3.3 觀測噪聲

觀測噪聲對模型參數的識別結果具有很大的影響，這種影響有時會使信號變形、失真，有的嚴重導致數據不可用。為了研究不同觀測噪聲對CKF算法識別模型參數的影響，本次模擬選用不同觀測噪聲對MFS模型參數進行識別，觀測噪聲Ws選用仿真觀測噪聲乘隨機系數，運行算法過程中算法初始參數預估值、初始誤差協方差矩陣不變，參數識別結果如圖10所示，其中單次識別值為單次使用CKF算法識別結果，統計識別值為CKF算法運行30次得到的參數識別值均值。

本次識別過程中，CKF單次識別平均耗時3.87 s，使用統計CKF算法總耗時120.35 s。在圖10中可以發現，使用不同觀測噪聲時，每次運行CKF算法識別模型參數的結果并不唯一，較模型參數真實值存在很大差異。從耗能的角度出發，單次使用CKF算法識別模

(a) 滯變位移識別值

(b) 加載位移識別值

(c) b識別值

(d) K識別值

(e) Q識別值

(f) 滯回曲線圖10 不同觀測噪聲的參數識別值結果Fig.10 Results of parameter identification values of different observation noises

型滯回曲線最大耗能與真實耗能值相對誤差為77.96%，而統計得到的滯回曲線耗能與真實耗能值相對誤差16.52%，耗能值相對誤差結果降低61.44%，單次使用CKF算法得到得模型滯回曲線最小耗能值與真實耗能值相對誤差為14.72%，相比使用統計CKF算法降低1.8%。經過統計后的參數識別值與真實值相比明顯誤差更小，吻合程度更好，穩定性高。使用不同觀測噪聲，模型參數識別值的均方根誤差如圖11所示。

圖11 參數識別值均方根誤差圖Fig.11 Root mean square error diagram of parameter identification values

通過誤差分析可以得出，相比單次使用CKF算法識別模型參數最大均方根誤差，統計CKF算法均方根誤差分別降低了18.39%、14.03%、136.05%。充分證明了當觀測噪聲不同時，統計CKF方法識別參數可以有效降低觀測噪聲對試驗的影響，從而保證了試驗的魯棒性。不同統計次數下，得到的模型滯回曲線如圖12所示。

為了檢驗統計次數對模型更新精度影響，分別對比統計10～50次模型滯回曲線，如圖12(a)所示。由圖12(a)可以看出，隨著統計次數的增加，統計50次得到的模型滯回曲線比統計10～40次得到的模型滯回曲線更能反映模型的真實滯回現象。由圖12(b)可以看出，從耗能的角度分析，統計50次得到的模型滯回曲線耗能值與真實耗能值相對誤差為14.26%，相比統計10～40次得到的滯回曲線耗能值相對誤差結果分別降低了53.22%、27.54%、2.16%、1.36%。證明統計方法得到的模型滯回曲線與真實滯回曲線更為接近，可以更好地模擬結構的真實反映。

(a) 滯回曲線

(b) 相對誤差圖12 不同統計次數得到的模型滯回曲線及相對誤差圖Fig.12 The hysteretic curves and relative error graphs of different statistical times

可以發現，統計30次較統計40次、50次相對誤差較小，可以得到相對較滿意的結果，因此后文數值模擬部分統計次數均采用統計30次的運行結果，以驗證統計方法的有效性。

3 基于統計CKF模型參數在線更新混合試驗數值模擬

3.1 數值模擬原理

為了進一步檢驗統計CKF方法模型更新混合試驗的精度及魯棒性，針對兩層帶有自復位摩擦耗能支撐結構采用MATLAB2018b編程進行混合試驗數值模擬。首先，選用下層自復位摩擦耗能支撐作為試驗子結構，上層自復位摩擦耗能支撐作為數值子結構，通過數值

圖13 基于統計CKF模型更新混合模擬方法示意圖Fig.13 Schematic diagram of hybrid simulation method based on statistical CKF model update

3.2 參數設置

設定每層支撐恢復力模型均采用本文2.1節MFS模型。試驗過程中假定主體結構處于線性狀態，框架結構模型參數設置為：每層結構質量Mn1=Mn2=2 000 t，框架結構的層間水平剛度剛度Kn1=Kn2=80 000 kN/m，框架結構的層間阻尼系數Cn1=Cn2=1 550 kN/(m/s)，其中n1、n2分別代表了數值子結構的一層和二層，支撐與樓面夾角均為28.81°。MFS真實模型參數分別為：激活位移b=0.001 4 m、初始剛度K=300 000 kN/m、耗能率Q=0.9、激活剛度比A=0.015。地震作用選取EI Centro(1940,NS)地震記錄，地震動峰值加速度為400 cm/s2，地震記錄如圖14所示。支撐滯變位移利用四階Runge-Kutta積分算法得出，積分步長為0.01 s。

圖14 EI Centro地震記錄Fig.14 EI Centro seismic record

假定狀態量為MFS模型滯變位移z、加載位移x、及模型參數b、K、Q，觀測量為MFS模型恢復力，利用識別關鍵參數值去更新數值子結構支撐部分參數。試驗子結構的狀態方程和觀測方程同式(28)～(29)所示，試驗結構加載速度采用式(32)位移差分方法確定。過程噪聲Vs和觀測噪聲Ws分別為均值為零的高斯白噪聲，過程噪聲Vs和觀測噪聲Ws的協方差矩陣分別為Rs、Qs，具體表達如式(36)～(37)所示。

Rs=10-26×I5

(36)

Qs=10-10×I1

(37)

(38)

(39)

3.3 自復位摩擦耗能支撐框架混合試驗數值模擬

(a) 下層支撐滯變位移識別值

(b) 下層支撐加載位移識別值

(c) b識別值

(d) K識別值

(e) Q識別值

(f) 下層支撐滯回曲線圖15 MFS模型參數識別值Fig.15 MFS model parameter identification values

上層支撐的滯回曲線數值模擬結果如圖16所示。

圖16 上層支撐滯回曲線數值模擬結果Fig.16 Numerical simulation results of upper bracing hysteretic curve

本次模擬過程中，CKF單次識別平均耗時6.71 s，使用統計CKF算法總耗時180.92 s。從圖15～16可以看出，單次使用CKF方法得到的模型參數識別值、模型滯回曲線與真實值有明顯差異，基于統計CKF方法得到的模型參數識別值、模型滯回曲線較真實值吻合程度較高。從耗能的角度分析，單次使用CKF方法得到的下層支撐模型滯回曲線最大耗能值較真實耗能值相對誤差18.2%，統計CKF方法得到的下層支撐模型滯回曲線耗能值與真實耗能值相對誤差為1.83%，耗能值相對誤差降低16.37%。單次使用CKF方法得到的上層支撐模型滯回曲線耗能值較真實耗能值相對誤差為48.91%，統計CKF方法得到的上層支撐模型滯回曲線耗能值較真實耗能值相對誤差為22.26%，耗能值相對誤差降低了26.65%。單次使用CKF算法得到得模型滯回曲線最小耗能值與真實耗能值相對誤差為18.55%，相比使用統計CKF算法降低3.71%。模型參數識別值的均方根誤差如圖17所示。

通過均方根誤差分析，相比單次使用CKF方法識別模型參數最大均方根誤差，統計CKF方法均方根誤差分別降低了30.84%、43.22%、25.08%，混合試驗模型參數識別精度均有很大改善。可見，統計CKF方法可以明顯改善參數識別結果，更具有魯棒性。

4 結 論

本文提出了一種基于統計CKF模型參數更新方法，通過對MFS模型進行不同初始參數條件下的參數識別，并針對兩層帶有自復位摩擦耗能支撐結構進行混合試驗數值仿真，可以得到以下結論：

(1) 在不同初始參數情況下，采用統計CKF算法對MFS模型進行參數識別，參數識別結果及滯回曲線與MFS模型真實值吻合程度更高，具有較高的識別精度。

(2) 采用統計CKF方法對兩層帶有自復位摩擦耗能支撐結構進行混合試驗數值仿真，可以弱化初始參數對混合試驗識別結果的影響，提高混合試驗模型參數在線識別精度及魯棒性。