楊紀(jì)鵬，夏 燁，孫利民,2

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

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

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

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

1 全局迭代參數(shù)卡爾曼濾波算法

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

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

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

式中：Qk和Rk分別為系統(tǒng)過程噪聲協(xié)方差矩陣和系統(tǒng)量測噪聲協(xié)方差矩陣；δkj為kronecker函數(shù)?；诓煌牟蓸硬呗?，本文所提三種方法過程如下。

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

3)時間更新

4)量測預(yù)測

因為結(jié)構(gòu)響應(yīng)與歷史路徑有關(guān)，故第k+1時刻的觀測值預(yù)測由第0時刻的初始狀態(tài)開始計算到當(dāng)前時刻，通常結(jié)構(gòu)初始狀態(tài)響應(yīng)取為0，地震動輸入已知，基于當(dāng)前的結(jié)構(gòu)參數(shù)，可用Newmark-β，Wilson-θ法等方法計算結(jié)構(gòu)動力響應(yīng)，從而得到當(dāng)前時刻的量測預(yù)測值。

5)量測更新

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

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

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

3)時間更新

4)量測預(yù)測

式(20)量測預(yù)測時，需從0時刻計算到當(dāng)前時刻。

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

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

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

曼濾波(GIP-SSRCQKF)

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

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

3)時間更新

4)量測預(yù)測

式(29)量測預(yù)測時，需從0時刻計算到當(dāng)前時刻。

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

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

以下對三種方法進(jìn)行簡要對比與總結(jié)：

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

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

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

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

所提算法流程可總結(jié)如圖1所示。

圖1 結(jié)構(gòu)參數(shù)識別流程圖Fig. 1 Flow chart of structural parameter identification process

2 數(shù)值仿真

2.1 雙線性模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

篇幅所限，選取GIP-CKF最終識別結(jié)果重構(gòu)結(jié)構(gòu)響應(yīng)，重構(gòu)結(jié)構(gòu)響應(yīng)及恢復(fù)力滯回圖如圖7所示。

圖7 雙線性模型結(jié)構(gòu)響應(yīng)重構(gòu)(GIP-CKF)Fig. 7 Reconstruction structural response of bilinear model (GIP-CKF)

由圖7可以看出，基于識別所得參數(shù)重構(gòu)響應(yīng)與理論值吻合非常好，即識別參數(shù)非常接近真值。

2.2 Takeda模型

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

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

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

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

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

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

圖9 Takeda恢復(fù)力模型Fig. 9 Restoring force of Takeda model

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

2)當(dāng)?shù)谝淮纬^屈服位移后，結(jié)構(gòu)剛度降為k2=αk1；

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

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

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

表2 Takeda模型初始參數(shù)表Table 2 Initial value of Takeda model

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

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

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

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

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

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

篇幅所限，選取GIP-CKF最終識別結(jié)果重構(gòu)結(jié)構(gòu)響應(yīng)，重構(gòu)結(jié)構(gòu)響應(yīng)及非線性恢復(fù)力滯回圖如圖13所示。

圖13 Takeda模型結(jié)構(gòu)響應(yīng)重構(gòu)(GIP-CKF)Fig. 13 Reconstruction structural response of Takeda model (GIP-CKF)

由圖13可以看出，對于Takeda模型，本文所提方法參數(shù)識別結(jié)果較好，重構(gòu)系統(tǒng)響應(yīng)與理論值吻合良好。

3 結(jié)論

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

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

(1)本文所提出的理念方法能夠?qū)崿F(xiàn)折線型本構(gòu)模型的參數(shù)識別?；谌N不同的Sigma點采樣規(guī)則，通過全局迭代，所提方法準(zhǔn)確識別雙線性和Takeda本構(gòu)模型，并具有較強(qiáng)的魯棒性。

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

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