基于Sigma點全局迭代參數卡爾曼濾波的折線型本構模型參數識別

楊紀鵬，夏 燁，孫利民,2

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

近年來，地震災害頻繁發生，土木工程結構在強震作用下因出現損傷而表現出剛度下降甚至倒塌破壞，因而地震下結構參數識別得到越來越多的關注，并發展出一些參數識別的方法。如最小二乘法(least sqaure method, LSM)[1?2]，擴展卡爾曼濾波方法(extended Kalman filter, EKF)[3?7]，聯合擴展卡爾曼粒子濾波-最小二乘算法[8]，中心差分卡爾曼濾波方法[9]，H∞濾波[10]，序貫蒙特卡羅方法[11]，無跡卡爾曼濾波(UKF)[12?14]等。除了選取合適的識別方法，合適的非線性模型選取也至關重要。根據非線性滯回模型的表達形式和形狀特征，這些模型基本上可以分為兩類，折線型本構模型和光滑型本構模型。

基于分段線性行為的模型稱為折線型本構模型，此類模型通常由構件或結構的實際行為驅動，如初始彈性階段，由于構件開裂、屈服、剛度和強度退化以及裂縫和間隙的閉合等引起的分段非線性行為特征。代表性模型如雙線性模型[15]、Clough型[16]、Takeda模型[17]、曲哲模型[18]、陸新征-曲哲模型[19]等。基于連續微分方程表示的非線性本構模型稱為光滑型本構模型，如Bouc-Wen模型[20?21]。連續型非線性模型因其數學表達式清晰，其求導及積分容易實現，相關研究比較多。折線型非線性模型表達式隨歷史路徑、當前狀態量變化而不斷變化，求偏導時其Jacobian矩陣很難求取，故相關研究較少。但折線型模型控制參數少，且具有明確的物理意義，基于識別的參數可用于區域結構震后快速評估[22]，因而對該類模型進行參數識別具有重要意義。相關研究如趙昕、李杰[23]提出加權全局迭代參數卡爾曼濾波算法，分別進行了單自由度線性系統和隨動強化雙線性結構模型的參數識別，該算法基于EKF，需要人為設定加權系數，當噪聲水平較高時，加權數選取不當將導致濾波結果發散。Lee、Yun[24]提出改進的序列擴展卡爾曼濾波算法進行鋼筋混凝土橋墩Takeda模型參數識別，計算過程中，狀態轉移矩陣采用有限差分近似得到，該方法計算過程復雜，且需要已知初始剛度和阻尼。張肖雄、賀佳[25]通過在觀測方程中引入投影矩陣，提出基于擴展卡爾曼濾波的結構參數和荷載識別的方法。在數值算例中使用兩折線分段線性模型，但在參數識別時假定劃分線性區間的臨界層間位移已知。而且該方法用于地震荷載識別時，結構相對地面的相對加速度是不可測的。

基于以上折線型本構模型參數識別存在的問題，本文提出基于Sigma點全局迭代參數卡爾曼濾波算法進行折線型模型參數識別。該方法利用Sigma點卡爾曼濾波不需要計算Jacobian矩陣的特點，嘗試折線型本構模型的參數識別。由于折線型非線性系統響應與歷史路徑有關，下一時刻系統響應不僅與當前時刻狀態量有關，還與所選非線性系統參數有關，即狀態向量中的系統響應量是沒有意義的，由此提出參數卡爾曼濾波算法，使狀態向量中不包含系統響應，降低狀態向量維度，減少計算量；在時間更新時，狀態量均由Sigma點直接給出，量測更新時，需由初始時刻計算到當前時刻。由于系統對不同參數的敏感程度不一樣，各參數收斂速度有所差異，使用全局迭代技術，即以上一輪迭代的終值作為下一輪迭代的初值；為擴大搜索范圍以獲得全局最優解，下一輪迭代時所用的協方差矩陣仍采用開始時刻設定的初始值。為驗證所提方法的有效性，將隔震支座系統、橋墩分別簡化為單自由度雙線性模型及單自由度Takeda模型進行參數識別。目前Sigma點采樣策略主要包括對稱采樣、單形采樣、3階矩偏度采樣及高斯分布4階矩對稱采樣等，其中普遍使用的是對稱采樣策略[26?29]。本文基于廣泛使用的UKF[12]、CKF[30]以及具有更高精度的SSRCQK[31]三種方法的采樣規則，分別對兩種非線性模型進行識別，并對三種方法的識別過程及結果進行對比分析。

1 全局迭代參數卡爾曼濾波算法

不同于由連續微分方程表示的非線性系統，折線型本構模型非線性系統，其響應與歷史路徑有關，數學表達式也更加復雜，使用傳統的EKF進行參數識別時，很難求解狀態轉移矩陣及量測矩陣的Jacobian矩陣。EKF的基本思想是利用泰勒展開，將非線性方程直接線性化，線性化的系統模型和系統實際的非線性模型會有差別，非線性越強，差別就會越大。UKF以Unscented變換來近似計算非線性系統狀態的后驗均值和協方差，從而避免對于非線性函數本身的近似，能以至少二階精度逼近任何非線性系統，較EKF更適合于強非線性系統，并且不需要計算Jacobian矩陣，如何選擇這幾個點去近似后驗均值和協方差的方法就叫Unscented變換，這幾個點選擇之后對均值和方差的影響不大，被選擇的點就稱為Sigma點，至于如何選擇這幾個點，已有較多的研究成果。本文選擇其中3種進行折線型本構模型識別。

基于Sigma點卡爾曼濾波的優點，本文提出基于Sigma點全局迭代參數卡爾曼濾波算法。該算法識別流程與Kalman濾波相同，而且只需一次采樣。由于生成Sigma點的規則不同，分別選取具有代表性的UKF、CKF以及具有5階精度的SSRCQKF驗證本文所提算法，即全局迭代參數無跡卡爾曼濾波(GIP-UKF)，全局迭代參數容積卡爾曼濾波(GIP-CKF)，全局迭代參數球面單純形徑向容積正交卡爾曼濾波(GIP-SSRCQKF)，以下對所提方法進行介紹。考慮如下離散的非線性模型：

式中：f(·)為非線性狀態函數，本文所提的方法中狀態向量僅包含待識別的結構參數，狀態函數進行狀態預測時，基于當前時刻的結構參數由起始時刻計算到當前時刻；h(·)為非線性量測函數；Xk+1為n維系統狀態向量；Zk+1為l維量測向量；uk為l維系統輸入向量，u0:k表示起始時刻到k時刻的輸入；wk為n維系統過程噪聲，vk+1為l維量測噪聲，wk和vk+1均為互不相關的零均值高斯白噪聲。且有：

式中：Qk和Rk分別為系統過程噪聲協方差矩陣和系統量測噪聲協方差矩陣；δkj為kronecker函數。基于不同的采樣策略，本文所提三種方法過程如下。

1.1 全局迭代參數無跡卡爾曼濾波(GIP-UKF)1)系統初始化

3)時間更新

4)量測預測

因為結構響應與歷史路徑有關，故第k+1時刻的觀測值預測由第0時刻的初始狀態開始計算到當前時刻，通常結構初始狀態響應取為0，地震動輸入已知，基于當前的結構參數，可用Newmark-β，Wilson-θ法等方法計算結構動力響應，從而得到當前時刻的量測預測值。

5)量測更新

1.2 全局迭代參數容積卡爾曼濾波(GIP-CKF)

1)系統初始化，與式(3)相同

2)選擇采樣策略，生成Sigma點，即容積點

3)時間更新

4)量測預測

式(20)量測預測時，需從0時刻計算到當前時刻。

5)量測更新，同式(12)～式(14)

式中，n為狀態向量維數，具有2n列的容積點集ξ如式(24)所示，ξi表示容積點集的第i列。

1.3 全局迭代參數球面單純形徑向容積正交卡爾

曼濾波(GIP-SSRCQKF)

1)系統初始化，同式(3)

2)選擇采樣策略，生成Sigma點

3)時間更新

4)量測預測

式(29)量測預測時，需從0時刻計算到當前時刻。

5)量測更新，同式(12)～式(14)

式(25)中，[a?a]i中下標i 表示矩陣的第i 列，矩陣中的元素計算如式(33)所示：

以下對三種方法進行簡要對比與總結：

1)三種采樣方法均為對稱采樣；

2) UKF使用最為廣泛，但缺乏嚴格的數學推導，其具有2n+1個Sigma點，其中一個為中心點，該中心點為上一時刻狀態向量的量測更新值，且具有較大的權重值。UKF進行參數識別時需要人為定義一些參數，這些參數對識別結果有較大影響，該方法通常具有3階精度。

3) CKF將難以處理的積分分解為球面積分和徑向積分，并用三階容積法則進行逼近。CKF具有2n個等權值的Cubature點(即Sigma點)，但其不包含初始原點。與UKF相比，其數學推導過程更為嚴謹，而且無須額外定義相關參數，其權重取值僅與狀態向量的維數有關，該方法通常具有3階精度。

4) SSRCQKF使用球面單純形徑向容積正交卡爾曼濾波代替球面積分，具有更高的精度，本文選取具有5階精度的采樣方法，對于n維狀態向量，其Sigma點數為4n+4。該方法計算更加耗時，通常識別效果也更好。

所提算法流程可總結如圖1所示。

圖1 結構參數識別流程圖Fig. 1 Flow chart of structural parameter identification process

2 數值仿真

2.1 雙線性模型

首先，考慮基礎隔震系統，將其簡化為單自由度結構(圖2)，假定在地震作用下隔振系統恢復力符合雙線性模型(圖3)，m為已知的集中質量，k1為初始剛度，c為阻尼系數。地震激勵選用EL Centro地震波，采樣頻率100 Hz，持續時間10 s(圖4)。

圖2 隔震支座系統模型簡化Fig. 2 Simplified model of isolation bearing system

圖3 雙線性恢復力模型Fig. 3 Bilinear restoring force model

圖4 地震波時程圖Fig. 4 Seismic wave time history

地震激勵下非線性恢復力模型動力方程為：

式中，x 為相對地面位移，x˙為相對地面速度，x¨為相對地面加速度，xg為地震動加速度。本算例中，所用參數取值如下，m=1×104kg ，k1=1×106N/m ，c=1×104N·s/m ，Xy=0.03 m，α=0.2，g(x)由下式給出：

待識別的狀態向量為X=[k1,c,Xy,α]T，Xy為屈服位移，α為屈服比。參數識別時，結構初始速度和位移取值為0，待識別參數初值取為真值的一半，GIP-UKF算法中αu=1/2；初始誤差協方差矩陣P0及系統過程噪聲協方差矩陣Q0對角線元素如表1所示。

表1 雙線性模型初始參數表Table 1 Initial value of bilinear model

觀測量為絕對加速度，為模擬實際監測中噪聲的干擾，在監測數據中加入均方根(RMS) 5%的高斯白噪聲，土木工程中結構基頻通常較低，故在參數識別時使用低通濾波器去除信號中的高頻噪聲，濾波截斷頻率為20 Hz，觀測誤差協方差矩陣R=1×10?2。參數收斂過程如圖5所示。

由圖5可以看出，對于雙線性模型：1) GIPSSRCQKF和GIP-UKF在第一次局部循環中識別結果已經收斂到真值附近；2) GIP-CKF第一次局部循環識別效果相對較差，但剛度參數收斂到真值附近，而阻尼參數、屈服位移以及屈服比則收斂于局部最優。

圖5 雙線性模型第一次迭代收斂結果Fig. 5 Convergence results of the first iteration of bilinear model

根據識別流程圖，以上一次迭代的終值作為下一次迭代的初值；為獲得全局最優解，協方差矩陣仍采用初始值。進行10次全局迭代，觀察整體收斂效果以及穩定性。全局參數收斂過程如圖6所示。

圖6 雙線性模型全局迭代收斂結果Fig. 6 Global iterative convergence results of bilinear model

由圖6可以看出：1)通過引入全局迭代，三種方法最終都能收斂到真值，其中GIP-SSRCQKF識別效果最好，且在第一步局部迭代就收斂到真值，而其他兩種方法也僅需兩次全局迭代收斂到真值；2) GIP-UKF和GIP-CKF兩種方法屈服比最終收斂結果略差；3)一次局部循環GIP-SSRCQKF所需時間最長，但該方法收斂到真值所需的全局迭代次數最少，故該方法是最優的。

篇幅所限，選取GIP-CKF最終識別結果重構結構響應，重構結構響應及恢復力滯回圖如圖7所示。

圖7 雙線性模型結構響應重構(GIP-CKF)Fig. 7 Reconstruction structural response of bilinear model (GIP-CKF)

由圖7可以看出，基于識別所得參數重構響應與理論值吻合非常好，即識別參數非常接近真值。

2.2 Takeda模型

當強震發生時，鋼筋混凝土(RC)橋墩上產生的力可能超過其屈服承載力，并導致橋墩發生較大的非彈性變形和損壞。由于橋梁結構的地震反應受受損鋼筋混凝土橋墩滯回性能的影響很大，因此需要建立可靠的滯回性能模型。改進的Takeda模型[33]能有效地再現鋼筋混凝土構件在有限參數下的復雜非線性滯回性能，因此在混凝土結構中得到了廣泛的應用。

橋墩在地震作用下的響應可以簡化為一單自由度系統(圖8)，假定在地震作用下橋墩恢復力符合Takeda模型(圖9)，地震下其動力方程如下：

圖8 橋墩模型簡化Fig. 8 Simplified model of bridge pier

式中：x、x˙、x¨分別表示相對地面位移、速度和加速度；m 為已知的集中質量；c 為阻尼系數；f(x,x˙)為非線性恢復力，地震激勵選用EL Centro地震波，采樣頻率100 Hz，持續時間10 s(圖10)。

圖10 地震波時程圖Fig. 10 Seismic wave time history

圖9中，k1為初始剛度，Xy為屈服位移，α為屈服比。Takeda模型規則簡要介紹如下：

圖9 Takeda恢復力模型Fig. 9 Restoring force of Takeda model

1)結構初始剛度為k1，若最大位移沒有超過屈服位移Xy(或?Xy)，則結構處于線性狀態，其剛度保持為k1；

2)當第一次超過屈服位移后，結構剛度降為k2=αk1；

3)線段2～3代表非線性卸載，其卸載剛度定義為k3=β1k1(φx/Dm)β0，β1為剛度折減系數，當卸載曲線在外圈時β1=1，當卸載曲線在內圈時β1=0.9，β0通常取值為0.4；Dm代表所在方向的最大位移，φx代表按照剛度k1卸載時與x軸的交點；

4)在反向加載的情況下，它指向上一個循環的歷史峰值，并且斜率小于初始剛度。如果在上一個周期中不存在這樣一個點，繼續搜索到上一個歷史峰值，直到找到這樣一個點。如果在相反方向上沒有超過屈服點，則直接指向屈服點。

單自由度Takeda模型中待識別的狀態向量為，X=[k1,c,Xy,α]T，本算例中，所用參數取值如下，m=1×104kg ，k1=1×106N/m ，c=1×104N·s/m ，Xy=0.015 m，α=0.1。參數識別時，結構初始速度和位移取值為0，待識別參數初值取為真值的一半，GIP-UKF算法中αu=1/20；初始誤差協方差矩陣P0及系統過程噪聲協方差矩陣Q0對角線元素如表2所示。

表2 Takeda模型初始參數表Table 2 Initial value of Takeda model

觀測量為絕對加速度，為模擬實際監測中噪聲的干擾，在監測數據中加入均方根(RMS) 5%的高斯白噪聲，土木工程中結構基頻通常較低，結構損傷后剛度下降，基頻會進一步降低，在參數識別時使用低通濾波器去除信號中的高頻噪聲，濾波截斷頻率為20 Hz，觀測誤差協方差矩陣R=1×10?2。第一次濾波識別結果如圖11所示。

由圖11可以看出，對于Takeda模型：GIPSSRCQKF方法第一次局部迭代就收斂到真值附近；GIP-CKF效果略微差于前者，基本也收斂到真值附近；GIP-UKF第一次局部迭代效果最差，初始剛度收斂到真值附近，阻尼參數、屈服位移、屈服比均收斂于局部最優解。

圖11 Takeda模型第一次迭代識別結果Fig. 11 Convergence results of the first iteration of Takeda model

為驗證全局迭代的效果及算法收斂的穩定性，以上一輪迭代的終值作為下一輪迭代的初值，進行10次全局迭代，整體識別結果如圖12所示。

由圖12可以看出，即使對于滯回規則更復雜的Takeda模型，所提方法經過全局迭代后也能收斂到真值附近。其中，1) GIP-SSRCQKF方法收斂最快，結果也最好；2) GIP-UKF識別過程中，因選擇參數較多，尤其是比例縮放因子αu的選取對識別結果有較大影響，故其全局收斂過程中波動幅度較大，但最終也能收斂于真值附近；3) GIP-CKF無須定義控制參數，雖然效果略微遜于GIP-SSRCQKF，但因其計算流程簡單，仍具有很大優勢。

圖12 Takeda模型全局迭代結果Fig. 12 Global iterative convergence results of Takeda model

篇幅所限，選取GIP-CKF最終識別結果重構結構響應，重構結構響應及非線性恢復力滯回圖如圖13所示。

圖13 Takeda模型結構響應重構(GIP-CKF)Fig. 13 Reconstruction structural response of Takeda model (GIP-CKF)

由圖13可以看出，對于Takeda模型，本文所提方法參數識別結果較好，重構系統響應與理論值吻合良好。

3 結論

本文基于Sigma點卡爾曼濾提出了全局迭代參數卡爾曼濾波算法。該方法利用Sigma點卡爾曼濾波避免求解非線性系統Jacobian矩陣，從而實現折線型本構模型參數識別。通過提出參數卡爾曼濾波，降低狀態向量維度，減少計算量，節約計算時間。需要注意的是量測更新時需要從初始時刻計算到當前時刻，最后通過全局迭代，得到待識別參數全局最優解。

為驗證所提方法的有效性，分別使用GIPUKF、GIP-CKF以及GIP-SSRCQKF進行雙線性及Takeda模型的參數識別，結論如下：

(1)本文所提出的理念方法能夠實現折線型本構模型的參數識別。基于三種不同的Sigma點采樣規則，通過全局迭代，所提方法準確識別雙線性和Takeda本構模型，并具有較強的魯棒性。

(2)三種方法中，GIP-SSRCQKF識別精度最高，雖然一次局部循環所需時間最長，但收斂到全局最優解所需的全局循環次數最少，故其總體計算時間成本并不高。GIP-UKF識別效果與GIPCKF差異不大，但前者選擇參數較多，不恰當的αu取值對結果有較大影響，選用該方法時需慎重選擇參數。GIP-CKF無須選擇計算參數，計算程序最簡單，經過全局迭代后也可以收斂到全局最優解，也具有較強的應用前景。

(3)針對特定的應用場景，本文僅研究了單自由度折線型本構模型參數識別。對于多自由度結構系統將更加復雜，而且折線型非線性模型中待識別參數均具有明確的物理意義與取值范圍，比如屈服位移、屈服比等參數不可能為負值，而在迭代過程中不能保證這些參數非負，因而會造成算法中斷，需引入其他約束方法進行進一步研究。