改進的強追蹤平方根無跡卡爾曼濾波時變結構參數識別

楊紀鵬，夏 燁，閆業祥，孫利民,2

(1.同濟大學 土木工程學院，上海 200092;2.同濟大學 土木工程防災國家重點實驗室，上海 200092)

土木工程結構在強地震荷載作用下由于結構損傷而引起系統參數突變，其包括結構剛度、阻尼以及描述系統行為的其他參數，因而考慮參數時變特性的結構系統識別得到越來越多的重視。

Yang等[1]提出自適應最小二乘法，該方法在識別結構參數時需要觀測結構的位移、速度和加速度，實際中地震作用下結構的位移和速度是很難觀測到的，而加速度是容易測得的物理量。

Yang等[2-3]提出自適應擴展卡爾曼濾波法，該方法基于擴展卡爾曼濾波法(extended Kalman filter，EKF)，在狀態誤差協方差矩陣中引入自適應矩陣，追蹤突變結構參數。擴展卡爾曼濾波是用泰勒展開式去線性近似狀態估計值和量測觀測值。該理念處理非線性模型的思想清晰，濾波過程簡單有效；但該方法也有明顯的缺陷：①進行泰勒級數近似展開時需要計算雅克比矩陣(Jacobian)，但對于有些非線性系統，其Jacobian矩陣不易計算；②通常EKF僅保留一階精度，若要提高到二階精度則要計算海塞矩陣，計算過程更加復雜；③因僅保留一階精度，通常需要較高的采樣頻率才能保證濾波效果。

雷鷹等[4]將靜力凝聚法與擴展卡爾曼濾波法相結合，識別結構參數，同時定量識別結構節點損傷程度，該方法基于擴展卡爾曼濾波，對于復雜系統求解Jacobian矩陣比較困難。

Julier等[5]根據確定性采樣提出的無跡卡爾曼濾波(unscented Kalman filter，UKF)以Unscented變換來近似計算系統狀態的后驗均值和協方差，相較于直接近似非線性函數本身，前者更容易實現，UKF能以至少二階精度逼近任何非線性系統。但傳統UKF存在協方差矩陣開方時矩陣奇異的數學問題。謝強等[6-7]提出奇異值分解(singular value decomposition,SVD)計算狀態誤差協方差矩陣的平方根，但效果欠佳。

池傳國等[8]提出一種基于M估計的強跟蹤SVD-UKF算法。該算法利用M估計理論，對異常新息數據進行“篩選”，保留有用新息，剔除有害新息，有效避免由衛星信號野值引起的粗差對強跟蹤SVD-UKF算法的魯棒性影響。

王小旭等[9-10]提出強追蹤無跡卡爾曼濾波器，該方法克服了傳統UKF無法追蹤時變結構參數的缺點，但其計算過程中需要計算狀態預測誤差協方差矩陣的平方根，其數值穩定性較差。

Van der Merwe等[11]提出平方根卡爾曼濾波，該方法基于UKF的算法框架，以協方差矩陣的平方根代替協方差矩陣進行遞推運算。該算法避免了時間狀態預測時協方差矩陣開方計算，但在濾波結果更新協方差矩陣平方根時，需要進行喬里斯基分解一階更新，其對矩陣的正定性的要求與UKF一樣高，并沒有從根本上解決數值計算時不穩定的問題。

張玉峰等[12]基于平方根濾波的思想，對傳統的Sage-Husa自適應濾波算法進行了改進，并與平方根無跡卡爾曼濾波(square root unscented Kalman filter,SRUKF)結合，該算法能直接對非線性系統的狀態方差矩陣和噪聲方差矩陣的平方根進行遞推與估算，確保狀態和噪聲方差矩陣的對稱性和非負定性，但文中仍需計算量測矩陣的Jacobian矩陣，未引入強追蹤濾波因子，不能有效追蹤時變系統參數。

葉浩澤等[13]在標準的平方根UKF算法上，首先改用了球型無跡變換對權系數以及Sigma點進行計算選??；其次改進了SRUKF中平方根矩陣的分解方法；同時在預測誤差協方差矩陣中引入了自適應衰減因子。該文中也需計算觀測矩陣的Jacobian矩陣，且文中計算強追蹤濾波因子的方法并不適用于土木結構系統參數識別。

李敏等[14]提出一種改進的強跟蹤平方根無跡卡爾曼濾波(strong tracking square root unscented Kalman filter，STSRUKF)導航方法，采用一種改進的平方根分解方法，改善了濾波器的穩定性。同時，基于強跟蹤濾波器理論，引入多重自適應衰減因子調節協方差矩陣，使得濾波器具有強跟蹤能力。數值模擬表明該方法效果并不好，同時該方法需要計算量測矩陣的Jacobian矩陣，文中所使用的強追蹤濾波因子算法也不適用于土木結構系統參數識別。

基于以上問題，本文提出改進的強追蹤平方根無跡卡爾曼濾波方法(modified strong tracing square root unscented Kalman filter ，MSTSRUKF)。該方法首先改進了協方差矩陣平方根計算方法，根據協方差定義直接使用QR分解計算其平方根，避免使用喬里斯基分解一階更新，使協方差矩陣平方根計算過程無條件數值穩定；其次，引入量測矩陣Jacobian矩陣的等價形式進行計算，對于非線性量測系統，可避免求解Jacobian矩陣，減小計算工作量，使該方法更具有普適性；最后，引入強追蹤濾波因子，根據土木工程結構參數特點，調整濾波因子矩陣對應項的比值，使該方法在引起較小參數波動的情況下即可識別到時變結構參數。

1 改進的平方根強追蹤無跡卡爾曼濾波

1.1 標準的SRUKF算法

考慮如下的離散非線性系統模型

(1)

式中：f(·)為非線性狀態函數；h(·)為非線性量測函數；Xk+1為n維系統狀態向量；Zk+1為m維量測向量；uk為s維系統輸入向量；wk為n維系統過程噪聲;vk+1為m維量測噪聲;wk和vk+1均為互不相關的零均值高斯白噪聲。且有

E(wk)=0,cov(wk,wj)=Qkδkj,

E(vk)=0,cov(vk,vj)=Rkδkj,

cov(wk,vj)=0

(2)

式中：Qk和Rk分別為系統過程噪聲協方差矩陣和系統量測噪聲協方差矩陣;δkj為kronecker函數?；赨KF的平方根無跡卡爾曼濾波(SRUKF)標準流程如下。

步驟1初始化狀態X0，狀態誤差協方差矩陣平方根S0

(3)

(4)

(5)

步驟3時間預測更新

(6)

(7)

(8)

(9)

步驟4量測更新

(10)

(11)

(12)

(13)

(14)

步驟5濾波結果更新

(15)

(16)

U=KkSZ

(17)

Sk+1|k+1=cholupdate{Sk+1|k,U,-1}

(18)

上述計算過程中用到的權值參數計算如下

(19)

式中:λ=α2(n+κ)-n;α為決定Sigma點在先驗均值附近擴展程度的主要尺度因子，通常取1×10-3<α≤1;β為第二個尺度因子，用于強調后驗協方差計算的第0個Sigma點的權重;κ為三級比例因子，通常取κ=3-n;n為狀態向量的維數;qr(·)為QR分解;cholupdate(·)為喬里斯基分解一階更新;上標T為矩陣的轉置。

1.2 改進的SRUKF算法

1.2.1 改進的時間預測狀態向量誤差協方差矩陣平方根計算

在標準的SRUKF算法中，由式(8)與式(9)計算狀態向量時間預測的誤差協方差矩陣的平方根，即式(9)可等價描述為

(20)

式(9)通過cholupdate計算Sk+1|k時，要求式(20)等式右邊必須是正定的，其要求與標準的UKF中計算誤差協方差矩陣平方根一樣嚴格，實際運算中由于舍入誤差的影響，很難保證矩陣的正定性，為此根據時間預測誤差協方差矩陣的定義

(21)

Sk+1|k=rT

(22)

又由式(21)推導得到

(23)

(24)

(25)

數值分析結果表明，式(25)對計算結果影響不大。

而李敏等提出如下的時間預測誤差協方差矩陣的計算方法

(26)

(27)

數值分析結果表明，式(26)、式(27)所提出的協方差計算方法參數識別結果并不理想。故本文提出如下協方差矩陣平方根計算方法

(28)

再通過式(24)計算Sk+1|k。

1.2.2 改進的量測誤差協方差平方根計算

同理，式(12)、式(13)可改進為

(29)

(30)

1.2.3 改進的濾波更新協方差矩陣計算

計算濾波結果更新狀態向量誤差協方差矩陣平方根時，式(17)、式(18)可等價描述為

(31)

這要求式(31)等式右邊必須是正定的，由于舍入誤差的影響，很容易使該部分失去正定性，因此還需對狀態估計誤差協方差矩陣平方根的計算做如下改進。

根據狀態估計誤差協方差的定義

(32)

將式(16)代入式(32)得

(33)

同時，式(1)中的量測方程可寫為

Zk+1=Hk+1Xk+1+Vk+1

(34)

(35)

將式(34)、式(35)代入式(33)可得

(36)

式中，I為單位矩陣。

由于Vk+1與其他向量不相關，根據協方差矩陣定義，由式(36)可推導得到

(37)

根據狀態向量時間預測誤差協方差矩陣平方根計算的改進方法，可得到改進的濾波更新狀態估計協方差矩陣的平方根計算方法

(38)

(39)

1.2.4 量測矩陣的改進

對于復雜的非線性結構，若觀測向量與狀態向量之間為非線性關系，為避免計算Jacobian矩陣，可使用量測矩陣Hk+1的等價形式。根據定義，狀態預測自協方差矩陣Pk+1|k，輸出預測自協方差矩陣PZZ,k+1和輸出預測互協方差矩陣PXZ,k+1可寫為式(40)～式(42)的形式

(40)

則可求得量測矩陣的等價形式為

(43)

1.3 強追蹤濾波算法

要使濾波器式(16)成為強追蹤濾波器的一個充分條件是實時在線調整增益矩陣Kk，使得

(44)

(45)

(46)

(47)

(48)

(49)

(50)

式中:tr[·]為矩陣的跡；l≥1為弱化因子；Γk+1為實際輸出殘差序列，可以表示為

(51)

設引入強追蹤漸消因子后的時間預測誤差協方差矩陣的平方根表示為Sk+1|k,new，計算出強追蹤濾波因子后，式(28)、式(24)表示的時間預測誤差協方差矩陣的平方根計算方法可更新為

(52)

(53)

根據以上推導，MSTSRUKF算法步驟可以描述為

步驟1根據式(3)、式(4)初始化系統狀態;

步驟3時間預測更新，根據式(6)、式(7)、式(28)、式(24)計算時間預測狀態向量及協方差矩陣的平方根;

步驟4量測更新，根據式(10)、式(11)、式(29)、式(30)、式(14)計算量測更新

步驟5計算漸消因子矩陣，更新時間預測誤差協方差矩陣的平方根；根據式(43)計算量測矩陣Jacobian矩陣等價形式，根據式(46)～式(51)計算漸消因子矩陣，根據式(52)、式(53)計算引入漸消因子后的時間預測誤差協方差矩陣的平方根;

步驟6濾波更新,根據式(15)計算增益矩陣，根據式(16)計算k+1時刻狀態后驗估計值；引入漸消因子后，k+1時刻狀態預測協方差矩陣根據式(54)～式(56)計算

(54)

(55)

(56)

2 數值算例驗證

2.1 三自由度線性結構系統

首先考慮一個三自由度線性結構系統遭受地震激勵，其運動方程可以寫為

(57)

本算例中，取m1=m2=m3=1 000 kg，k1=120 kN/m，k2=120 kN/m，k3=60 kN/m。c1=c2=c3=0.6 kN·s/m，地震激勵選取El Centro地震波，采樣頻率100 Hz，地震持續時間31.2 s。

為模擬實際觀測中噪聲的干擾，在模擬的量測加速度數據中加入5%RMS(root mean square)的高斯白噪聲，然后使用低通濾波去除加速度時程中的高頻成分，濾波截斷頻率為20 Hz，為了驗證MSTSRUKF在識別線性結構參數突變時的有效性，設定當t=10 s時:結構一層剛度由k1=120 kN/m突變到k1=80 kN/m，二層剛度k2=120 kN/m突變到k2=80 kN/m，三層剛度k3=60 kN/m突變到k3=40 kN/m；結構一層阻尼由c1=0.6 kN·s/m突變到c1=0.7 kN·s/m，二層阻尼c2=0.6 kN·s/m突變到c2=0.65 kN·s/m，三層阻尼c3=0.6 kN·s/m突變到c3=0.65 kN·s/m。狀態向量初值取為X0=[0,0,0,0,0,0,90,90,50,0.3,0.3,0.3]T

結構參數識別結果,如圖1所示。結構各層位移識別結果,如圖2所示。

若以1～10 s為第一階段，10.0～31.2 s為第二階段，不同階段參數最終識別結果如表1所示。

由圖1、圖2及表1可以看出，本文提出的MSTSRUKF：①能夠快速有效地識別突變的結構參數；②由于剛度參數的數量級遠大于阻尼參數，其在協方差矩陣中起主要作用，故剛度參數識別精度高于阻尼參數，且阻尼參數震蕩幅度更大，收斂速度慢于剛度參數；③能夠準確估計出結構位移時程；④由各個階段最終識別結果看出，即使加入5%白噪聲，最終收斂精度仍然很高，剛度參數識別誤差不大于0.2%，阻尼參數識別誤差不大于4%，具有較高的魯棒性；⑤所提出算法數值穩定，不存在因矩陣奇異而計算中斷的問題。

表1 線性系統時變參數識別結果Tab.1 Time-varying parameter identification results of linear system

2.2 三自由度Bouc-Wen模型系統

考慮一個三自由度剪切型框架結構，底層采用Bouc-Wen滯回模型，受地震荷載激勵，其運動方程可寫為

(59)

(a)

本算例中，取m1=450 kg，m2=400.5 kg，m3=350.5 kg；k1=20.5 kN/m，k2=23.5 kN/m，k3=23.5 kN/m；c1=0.205 kN·s/m，c2=0.255 kN·s/m，c3=0.255 kN·s/m；α=0.2，β=2，γ=1，n=1，由于該模型參數的冗余性，n不作為待識別參數。地震激勵選取El Centro地震波，采樣頻率100 Hz，地震持續時間31.2 s。

(a)

本文中認為地震輸入已知，觀測量為各層絕對加速度。為模擬實際觀測中噪聲的干擾，在模擬的量測加速度數據中加入5%RMS的高斯白噪聲，然后使用低通濾波去除加速度時程中的高頻成分，濾波截斷頻率為20 Hz，為了驗證MSTSRUKF在識別非線性結構參數突變的有效性，設定當t=15 s時：結構一層剛度由k1=20.5 kN/m突變到k1=14 kN/m，結構二層剛度k2=23.5 kN/m突變到k2=16 kN/m，三層剛度k3=23.5 kN/m突變到k3=16 kN/m。狀態向量初值取為X0=[0,0,0,0,0,0,0,16,16,16,0.16,0.16,0.16,0.16,1.6,0.8]T

結構參數識別結果，如圖3所示。結構各層位移識別結果，如圖4所示。第一層滯回曲線及滯變位移識別結果，分別如圖5、圖6所示。

(a)

(a)

圖5 第一層滯回曲線識別結果Fig.5 Hysteric loop estimation results of 1st floor

圖6 第一層滯變位移識別結果Fig.6 Hysteric displacement estimation results of 1st floor

若以1～15 s為第一階段，15.0～31.2 s為第二階段，不同階段參數最終識別結果,如表2所示。

由圖3、圖4、圖5、圖6以及表2可以看出，本文提出的MSTSRUKF：①對于非線性系統具有良好的識別效果，能較為準確地識別各層位移；②能夠有效識別到參數突變，剛度參數在經過短暫震蕩后很快收斂到真值；③根據先驗知識調整阻尼項在漸消因子中對應μk值，阻尼項經較小幅度的震蕩后，也很快收斂到真值；④由于剛度參數的數量級遠大于阻尼參數的數量級，故剛度參數的識別精度高于阻尼參數；⑤結構參數突變后，描述Bouc-Wen模型系統的參數識別精度有所降低，最高誤差達到-11.19%；⑥所提方法數值計算穩定，對于強非線性系統仍能穩定、有效運行。

表2 Bouc-Wen模型系統時變參數識別結果Tab.2 Time-varying parameter identification results of Bouc-Wen model system

2.3 三自由度杜芬結構系統

考慮一個三自由度剪切型杜芬系統遭受地震激勵，其運動方程可以寫為

(60)

本算例中取m1=m2=m3=1 000 kg，c1=c2=c3=0.6 kN·s/m，k1=120 kN/m，k2=120 kN/m，k3=60 kN/m，kd1=200 kN/m3，kd2=200 kN/m3，kd3=-50 kN/m3。地震激勵選取El Centro地震波，采樣頻率100 Hz，地震持續時間31.2 s。

本文中認為地震輸入已知，觀測量為各層絕對加速度。為模擬實際觀測中噪聲的干擾，在模擬的量測加速度數據中加入5%RMS的高斯白噪聲，然后使用低通濾波去除加速度時程中的高頻成分，濾波截斷頻率為20 Hz，為了驗證MSTSRUKF在識別非線性結構參數突變的有效性，設定當t=15 s時，結構一層剛度由k1=120 kN/m突變到k1=100 kN/m，結構二層剛度k2=120 kN/m突變到k2=100 kN/m。狀態向量初值取為X0=[0,0,0,0,0,0,90,90,50,0,0,0,0.3,0.3,0.3]T。

結構參數識別結果,如圖7所示，結構各層位移識別結果,如圖8所示。

(a)

(a)

若以1～15 s為第一階段，15.0～31.2 s為第二階段，不同階段參數最終識別結果,如表3所示。

表3 Duffing系統時變參數識別結果Tab.3 Time-varying parameter identification results of Duffing system

由圖7、圖8及表3可以看出，本文提出的MSTSRUKF：①能夠準確估計Duffing系統位移時程；②能夠有效追蹤到結構剛度突變，并很快收斂到真值；③根據先驗知識調整阻尼項在漸消因子中對應μk值，阻尼參數微小震蕩后快速收斂到真值；④與位移三次方相關的剛度參數識別誤差略微偏大，最大誤差為5.54%；⑤計算過程數值穩定，不存在因數學計算不穩定而意外終止的問題。

3 結 論

針對土木工程結構地震作用下時變參數識別問題，提出改進的強追蹤平方根無跡卡爾曼濾波算法。算法克服了傳統UKF及SRUKF協方差矩陣平方根遞推過程中數值計算不穩定的問題，使遞推過程無條件數值穩定；其次使用量測矩陣Jacobian矩陣的等價形式，避免復雜系統求解Jacobian矩陣，簡化計算方法，擴展該方法的應用范圍；最后引入強追蹤濾波技術，實現時變結構的參數識別。數值分析結果表明：

(1)線性系統數學模型相對簡單，其時變參數識別結果優于非線性系統，因剛度參數與阻尼參數數量級相差較大，剛度參數識別結果優于阻尼參數。

(2)對于Bouc-Wen模型系統，由于控制滯回特性的參數較多，而系統對不同參數的敏感性也不同，故應科學選取敏感參數的初值及協方差矩陣中對應的元素值，避免識別過程中系統發散。

(3)對于杜芬系統，與位移三次方對應的剛度參數識別結果精度低于與一次項對應的剛度參數，說明其對系統的影響程度不如一次項對應的剛度參數。

(4)對于時變參數的線性系統和非線性系統，所提方法均能較好地識別結構參數，結構狀態預測結果在初始階段有一定誤差，待參數收斂后預測結構狀態也與真實狀態吻合良好。

(5)在數值模擬過程中，所提出的算法數值計算穩定，最終識別結果表明本文提出的方法在避免求解Jacobian矩陣后能夠準確識別結構參數。

針對時時參數突變的情況因收斂需要過程，當結構剛度、阻尼以及控制滯回特性的參數變化在合理范圍內時，所提算法能夠收斂；本算法確保了遞推過程計算的穩定性，但實際工程中觀測量數目遠小于結構的自由度數目，實際應用效果還需進一步驗證。