基于改進IMK恢復力模型的鋼筋混凝土柱參數(shù)識別與應用

郭玉榮 龍沐恩

摘? ?要：提出了一種利用鋼筋混凝土柱擬靜力試驗數(shù)據(jù)識別改進IMK模型骨架曲線參數(shù)，進而提高鋼筋混凝土框架結構非線性模擬精度的方法. 通過引入可抗差的基于奇異值分解的無跡卡爾曼濾波算法（抗差SVD-UKF算法），抑制觀測值粗差對參數(shù)識別的影響，采用粒子群算法對初始協(xié)方差矩陣、過程噪聲矩陣和測量噪聲矩陣進行自動尋優(yōu)，在MATLAB中實現(xiàn)了柱滯回特征正負向對稱與非對稱兩種情況下改進IMK恢復力模型骨架曲線參數(shù)的識別. 鋼筋混凝土柱實測滯回曲線的模型骨架曲線參數(shù)識別結果及其在框架結構非線性模擬中的應用結果驗證了本文方法的有效性.

關鍵詞：恢復力模型;滯回特征;參數(shù)識別;抗差SVD-UKF算法;粒子群算法

中圖分類號：TU317? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 文獻標志碼：A

文章編號：1674—2974（2021）01—0126—09

Abstract：A method for identifying the backbone curve parameters of the modified Ibarra-Medina-Krawinkler （IMK） model by using quasi-static test data of reinforced concrete columns and thus improving the simulation accuracy of reinforced concrete frame structures is proposed in this paper. In this method， a robust unscented Kalman filtering algorithm based on singular value decomposition（robust SVD-UKF algorithm） is introduced to suppress the influence of gross error of the observation on the parameter identification， and the particle swarm optimization algorithm is adopted to automatically optimize the initial covariance matrix， the process and measurement noise matrices. The identification of backbone curve parameters of the modified IMK model is realized using MATLAB，in which the symmetric and asymmetric hysteresis behavior of the columns in the positive and negative direction is considered. The effectiveness of the proposed method is verified by model backbone curve parameter identification based on the measured hysteretic curves of reinforced concrete columns and its application in the nonlinear simulation of frame structures.

Key words：hysteretic model;hysteretic behavior;parameter identification;robust SVD-UKF algorithm;particle swarm optimization algorithm

鋼筋混凝土框架結構的地震響應混合模擬[1] 及其抗倒塌性能的分析，需要可有效模擬鋼筋混凝土構件滯回特征[2]的恢復力模型及準確的模型參數(shù). 塑性鉸模型是框架結構非線性模擬常采用的一種模型，它不僅反映構件的力學特征，還與構件的材料、約束狀況及空間布局密切相關. 幾十年以來，塑性鉸模型已有了飛速發(fā)展，Clough等[3]開發(fā)了雙線性模型;Wen[4]提出了光滑的塑性鉸模型;Takeda等[5]開發(fā)了三線性塑性鉸模型. 但是上述常見模型不能充分考慮構件在循環(huán)往復荷載作用下的剛度和強度退化，影響了整體結構模擬分析的精確性. Ibarra等[6-7]開發(fā)了復雜的塑性鉸模型即改進的IMK （Ibarra-Medina-Krawinkler）模型. 改進IMK模型作為一種以三折線為骨架曲線的塑性鉸模型，引入了基于能量耗散的退化參數(shù)β，考慮了構件在往復荷載作用下的多種剛度和強度退化，相對于其他塑性鉸模型，能夠較好地對鋼筋混凝土梁柱的滯回特征進行有效模擬. 但根據(jù)模型經(jīng)驗公式直接計算的改進IMK模型骨架曲線參數(shù)可能會存在明顯的誤差，從而影響框架結構地震響應模擬精度或導致結構抗倒塌性能分析產(chǎn)生較大誤差，因此，在選擇恢復力模型之后，模型參數(shù)的準確度成為影響結構體系模擬精度的主要因素之一. 為此，本文提出一種結合擬靜力試驗的混合模擬方法，先從框架結構中選取部分關鍵構件進行擬靜力試驗，然后基于構件實測滯回曲線對構件恢復力模型骨架曲線參數(shù)進行識別，最后將識別的參數(shù)用于更新結構中相同構件的恢復力模型參數(shù)，并進行整體結構非線性數(shù)值模擬. 在這一混合模擬方法中，模型參數(shù)識別是最關鍵的環(huán)節(jié)之一.

目前在土木工程領域廣泛使用的參數(shù)識別方法主要有最小二乘估計法（Least Square Estimation，LSE）、擴展卡爾曼濾波法（Extended Kalman Filter，EKF）和無跡卡爾曼濾波法（Unscented Kalman Filter，UKF）[8]，其中UKF算法以其高效性成為非線性參數(shù)識別的常用方法. 但UKF方法在實際運用時仍存在一些局限性，系統(tǒng)狀態(tài)先驗信息矩陣會因為觀測粗差、系統(tǒng)噪聲等不確定的因素失去正定性，中斷濾波進程，觀測粗差等因素還會影響濾波過程，影響收斂速度和識別結果;此外，在UKF方法中初始協(xié)方差矩陣及過程噪聲和測量噪聲矩陣的確定較為繁瑣，但它們的取值卻又決定著參數(shù)識別的準確性. 基于此，本文運用可抗差的基于奇異值分解的無跡卡爾曼濾波算法（robust unscented Kalman filtering algorithm based on singular value decomposition，抗差SVD-UKF算法）[9]來解決狀態(tài)先驗信息矩陣失去正定性的問題并抑制粗差等對濾波過程的影響，并通過粒子群算法對初始協(xié)方差矩陣及過程噪聲矩陣和測量噪聲矩陣進行自動尋優(yōu)，采用MATLAB實現(xiàn)了改進IMK恢復力模型骨架曲線參數(shù)的識別. 通過基于陸新征等的框架柱擬靜力試驗數(shù)據(jù)[10]的模型參數(shù)識別以及對其整體框架結構擬靜力試驗結果[11]的模擬對比，來驗證本文方法的有效性.

1? ?改進的IMK恢復力模型

改進的IMK恢復力模型是一種用塑性鉸來模擬梁柱構件非線性行為的恢復力模型，它采用三折線骨架曲線（如圖1所示），其骨架曲線形狀由5個參數(shù)EIy、θcap，p l、θpc、My、Mc /My確定，參數(shù)的計算如式（1）～式（5）[12]所示，其循環(huán)能量耗散能力λ的計算如公式（6）所示.

式中：EIy為割線模量;θcap，p l為塑性轉角;θpc為峰值后轉角;My為屈服彎矩值;Mc /My為屈服后硬化剛度;P/Ag? f ′ c、ν均為軸壓比;？準y和ky為由構件基本信息計算出的參數(shù)，在此不作展開;δ′= d′/d，d′為受壓區(qū)邊緣到受壓鋼筋中心的距離;s為塑性鉸區(qū)的箍筋間距;asl為縱向鋼筋滑移系數(shù)（考慮取1，不考慮取0）;sn為鋼筋屈曲系數(shù)，sn=（s/db）（fy /100）0.5，db為縱筋直徑，fy為縱筋屈服強度;f ′ c為混凝土軸心抗壓強度，單位為MPa;cunits是單位換算量，當f ′ c的單位為MPa時cunits取1;ρsh為柱塑性鉸區(qū)的橫向鋼筋面積比;ρ為縱筋配筋率;b和d分別為柱橫截面的寬度和高度.

2? ?無跡卡爾曼濾波算法

2.1? ?標準UKF算法

2.1.1? ?時間預測

利用第K-1步狀態(tài)向量構造一組總數(shù)為2L+1的sigma點集，L為狀態(tài)量的維度，X在本文中為待識別的參數(shù) K1、K2、K3、 fy、 fp和恢復力r組成的狀態(tài)向量：

2.2? ?抗差SVD-UKF算法

在觀測粗差、計算機截斷誤差等因素的干擾下，狀態(tài)協(xié)方差矩陣可能會失去正定性，導致無損變換（Unscented Transformation，UT）中的Cholesky分解失效. 本文采用SVD分解生成sigma點集，之后在UKF算法中通過抗差因子調節(jié)噪聲矩陣來抑制觀測粗差等因素對濾波過程的影響，該算法的流程圖如圖2所示.

2.2.1? ?抗差模型

抗差因子將觀測值作為考量因素引入濾波過程中，起到動態(tài)調整測量噪聲的作用. 采用觀測殘差構建抗差因子，抗差因子θk的表達式與IGGIII函數(shù)的形式相似，其表達式如下：

式中：k0、k1為閾值參數(shù)，k0通常取1.0 ～ 2.0，k1取3.0 ～8.5，本文中k0取1，k1取3;vk為當前步觀測值與估計值的殘差; vk為當前步觀測值與估計值殘差的均值;svk、σk分別為當前步觀測值與估計值殘差的標準化量和觀測值與估計值殘差的樣本標準差.

2.2.2? ?SVD分解與觀測協(xié)方差矩陣的調整

為了解決在異常情況下狀態(tài)協(xié)方差矩陣失去正定性而無法使用Cholesky分解的問題，引入SVD分解計算sigma點集，計算公式如下：

2.2.3? ?粒子群算法的自動尋優(yōu)

粒子群算法是模擬鳥群隨機搜尋食物的捕食行為，鳥群通過自身經(jīng)驗和種群之間的交流調整自己的搜尋路徑，從而找到食物最多的地點，其中每只鳥的位置為自變量組合，每次到達地點的食物密度即為函數(shù)值. 每次搜尋都會根據(jù)自身經(jīng)驗（自身歷史搜尋的最優(yōu)地點）和種群交流（種群歷史搜尋的最優(yōu)地點）調整自身的搜尋方向和速度，這稱為跟蹤極值，從而找到最優(yōu)解. 粒子群算法的流程圖如圖3所示.

本文中定義粒子的位置為初始協(xié)方差矩陣的每個元素以及過程噪聲和測量噪聲矩陣的每個元素，即相當于構造了一個N維的空間，并讓鳥群在這個N維空間中搜尋使適應度函數(shù)取值最小的坐標解，適應度函數(shù)為在每一個粒子自變量組合經(jīng)過SVD-UKF算法識別參數(shù)后，由識別參數(shù)計算所得的每個恢復力值與相應試驗值的殘差平方和，殘差平方和越大表示識別的結果越不可靠，反之則識別結果越可靠，也即粒子自身包含的自變量組合的取值越合理. 狀態(tài)轉移方程是表示粒子由當前時刻至下一時刻的位置變化方式，狀態(tài)轉移方程如下：

2.2.4? ?抗差SVD-UKF算法抗差效果分析

基于改進IMK三折線模型的MATLAB仿真，以表1參數(shù)為仿真參數(shù)值，通過輸入位移峰值依次為4、10 、15 、20 、25 、30 、40 、45 、55 mm以及各級位移峰值循環(huán)兩圈的加載序列，計算得到恢復力序列，之后為檢驗抗差效果將兩類異常觀測引入仿真結果.

1）單點異常：N（0，2Z（i）2）（均值為0，方差為2Z（i）2）的單點隨機誤差，其中Z（i）為當前恢復力值.

2）連續(xù)段異常：N（0，2Z（i）2）的連續(xù)段隨機誤差.

仿真結果加入異常后的觀測值如圖4所示，之后分別采用未抗差的UKF算法及抗差的SVD-UKF算法對加入異常的觀測恢復力值進行濾波，濾波過程如圖5所示.

圖6為在引入單點異常后截取的2 000～4 000加載步的濾波過程圖. 由圖6可見，抗差算法的濾波過程跟真實值基本吻合，而未抗差算法的濾波過程受到單點異常影響而出現(xiàn)了偏離.

誤差后截取的5 000～7 000加載步的濾波過程圖. 由圖7可知，抗差算法的濾波過程跟真實值基本吻合，而未抗差算法的濾波過程受到了連續(xù)段異常的影響而出現(xiàn)了較大波動.

可見，在兩種異常的影響下抗差SVD-UKF算法擁有較強的魯棒性與吻合程度.

3? ?低周反復荷載下鋼筋混凝土柱的骨架曲線參數(shù)識別

3.1? ?試驗概況

本文試驗實測數(shù)據(jù)采自陸新征等[10-11]鋼筋混凝土框架結構擬靜力倒塌試驗研究，試驗包括一榀3層3跨整體框架的擬靜力試驗和一些關鍵構件的擬靜力試驗. 關鍵構件包括底層邊柱和中柱各兩個試件，層一的邊節(jié)點和中節(jié)點各一個試件，本節(jié)參數(shù)識別選用了其中邊柱A和中柱C的試驗數(shù)據(jù). 邊柱A和中柱C的幾何尺寸和配筋如圖8所示，其中中柱角筋直徑為10 mm，鋼筋材料參數(shù)如表2所示，混凝土強度和對柱子施加的軸力如表3所示. 柱子擬靜力試驗的加載路徑如下.

邊柱A試驗中首先施加柱端豎向荷載至試驗軸壓比，然后以水平荷載控制加載10 kN/級至30 kN，各級荷載循環(huán)一圈. 此后以水平加載點位移控制加載，所加水平位移峰值分別為10、15、20、25、30、37.5、45、55 mm，各級位移循環(huán)兩圈.

中柱C試驗中首先施加柱端豎向荷載至試驗軸壓比，之后以水平荷載控制加載10 kN/級至40 kN，各級荷載循環(huán)一圈. 此后以水平加載點位移控制加載，所加水平位移峰值依次為10、15、20、25、30、37.5、45、55 mm，各級位移循環(huán)兩圈.

3.2? ?基于柱試驗結果的改進IMK對稱模型參數(shù)識別

本文采用抗差SVD-UKF算法，在MATLAB中編寫了識別改進IMK對稱模型骨架曲線參數(shù)的程序，識別參數(shù)包括彈性剛度K1、硬化剛度K2、退化剛度K3、屈服力fy和極限力fp這5個參數(shù). 在后續(xù)的滯回曲線對比中，試驗結果為柱子擬靜力試驗的實測滯回曲線;初始結果是指通過公式（1）～（6）計算得出改進IMK模型參數(shù)值后，依據(jù)柱子的位移加載路徑計算得到的滯回曲線;以初始結果中計算得到的模型參數(shù)為參數(shù)識別的初始值，識別出5個待識別參數(shù)后，依據(jù)柱子的位移加載路徑計算得到的滯回曲線稱為識別結果.

邊柱A的改進IMK模型骨架曲線參數(shù)的初始值和識別值對比如表4所示. 用初始值和識別值分別計算出柱子的滯回曲線，邊柱A計算滯回曲線和試驗滯回曲線的對比分別如圖9和圖10所示，與初始結果相比，識別結果與試驗結果吻合較好，也即通過參數(shù)識別提高了模型參數(shù)的準確度.

圖11為邊柱A的試驗結果、識別結果和初始結果的滯回環(huán)面積對比. 由圖11的滯回環(huán)面積數(shù)據(jù)計算得到邊柱A的識別結果滯回環(huán)總面積與試驗結果滯回環(huán)總面積相差8.15%，邊柱A初始結果滯回環(huán)總面積與試驗結果滯回環(huán)總面積相差15.77%. 從滯回環(huán)總面積差來看，與初始參數(shù)相比，采用識別參數(shù)可以更好地模擬柱構件的滯回曲線.

3.3? ?基于柱試驗結果的改進IMK非對稱模型參數(shù)識別

在上述識別方法中，通過修改MATLAB中的程序，實現(xiàn)了改進IMK非對稱模型骨架曲線參數(shù)的識別.

中柱C的改進IMK模型骨架曲線參數(shù)初始值和識別值對比如表5所示. 用初始值和識別值分別計算出柱子的滯回曲線，計算滯回曲線和試驗滯回曲線的對比分別如圖12和圖13所示. 與初始結果相比，識別結果與試驗結果吻合較好，也即通過參數(shù)識別提高了模型參數(shù)的準確度.

圖14為中柱C試驗結果、識別結果和初始結果的滯回環(huán)面積對比. 由圖14的滯回環(huán)面積數(shù)據(jù)計算得到識別結果滯回環(huán)總面積與試驗結果滯回環(huán)總面積相差26.03%，初始結果滯回環(huán)總面積與試驗結果滯回環(huán)總面積相差31.79%. 可見，與初始參數(shù)相比，采用識別參數(shù)可更好地模擬柱構件的滯回曲線.

4? ?基于柱構件參數(shù)識別結果的應用

為檢驗基于柱實測滯回曲線的模型骨架曲線參數(shù)識別值的有效性，以下將識別的模型參數(shù)應用于整體框架結構的數(shù)值模擬，所選的框架結構來自陸新征等[10-11]完成的一榀3層3跨鋼筋混凝土框架結構擬靜力倒塌試驗. 整體框架的數(shù)值模擬采用OpenSees建模，建模流程如下.

首先計算得到改進IMK模型的輸入?yún)?shù)，然后在OpenSees中調用改進IMK模型并定義塑性鉸材料，最后將材料用于梁柱端部零長度單元的彎曲特性，其余單元均為彈性梁柱單元. 柱端部改進IMK模型參數(shù)按兩種工況取值換算，工況1采用由柱端截面的尺寸、材料、配筋、軸力等參數(shù)按照公式（1）～（6）計算得到的模型參數(shù)值;工況2采用基于柱實測試驗數(shù)據(jù)得到的參數(shù)識別值，其中中柱的參數(shù)值選用中柱C負向的參數(shù)識別值，將參數(shù)識別值換算成模型參數(shù)值后，對割線剛度EIy、塑性轉角θcap，p l和峰值后轉角θpc根據(jù)框架的實際軸力按照公式（1）（3）（4）中軸壓比的分項作了調整，考慮到屈服彎矩與軸壓比關系的復雜性，對屈服彎矩My和屈服后硬化剛度Mc /My分別按公式（2）和公式（5）計算取值，對λ取一較大值，即不考慮構件的循環(huán)能量耗散. 之后，根據(jù)擬靜力試驗設定的加載路徑，進行了兩種工況下框架結構的數(shù)值模擬. 圖15為兩種工況下基底剪力-頂點位移數(shù)值模擬滯回曲線和試驗滯回曲線的對比. 圖16為工況1、工況2與試驗結果的滯回環(huán)面積對比.

從圖15曲線對比可以看出，工況1與工況2均對試驗情況作出了較好的模擬. 由圖16的滯回環(huán)面積數(shù)據(jù)計算得到基底剪力-頂點位移識別結果的滯回環(huán)總面積與試驗結果的滯回環(huán)總面積相差23.71%，初始結果滯回環(huán)總面積與試驗結果滯回環(huán)總面積相差36.50%. 從滯回環(huán)總面積差來看，與初始參數(shù)相比，采用識別參數(shù)可以更好地模擬框架結構的整體響應.

5? ?結? ?論

1）將抗差SVD-UKF算法運用于鋼筋混凝土柱改進IMK模型的參數(shù)識別，有效抑制了觀測粗差對識別進程的影響，并通過粒子群算法對初始協(xié)方差矩陣、過程噪聲矩陣、測量噪聲矩陣進行自動尋優(yōu)，省去了大量的人力調參工作.

2）基于陸新征等[10-11]鋼筋混凝土框架柱構件擬靜力試驗實測滯回曲線數(shù)據(jù)，對柱子的改進IMK模型骨架曲線參數(shù)進行了識別，結果表明識別參數(shù)精度較初始參數(shù)精度有明顯提高.

3）將基于柱子試驗數(shù)據(jù)識別的骨架曲線參數(shù)用于框架結構的擬靜力加載試驗數(shù)值模擬，結果表明利用識別參數(shù)可提高整體結構的數(shù)值模擬精度.

參考文獻

[1]? ? 郭玉榮，劉鐘真. 基于OpenSees的子結構擬動力試驗方法研究[J]. 湖南大學學報（自然科學版），2015，42（9）：1—7.

GUO Y R，LIU Z Z. Research on substructure pseudo-dynamic test method based on OpenSees[J]. Journal of Hunan University （Natural Sciences），2015，42（9）：1—7. （In Chinese）

[2]? ? 過鎮(zhèn)海. 鋼筋混凝土原理[M]. 北京：清華大學出版社，1999：147—178，192—212.

GUO Z H. Reinforcement concrete theory[M]. Beijing： Tsinghua University Press，1999：147—178，192—212. （In Chinese）

[3]? ? CLOUGH R W，JOHNSTON S B. Effect of stiffness degradation on earthquake ductility requirements[C]//Proceedings of the 2nd Japan earthquake engineering symposium. Tokyo，1966.

[4]? ? WEN Y K. Method for random vibration of hysteretic systems[J]. Journal of the Engineering Mechanics Division，1976，102（2）：249—263.

[5]? ? TAKEDA T，SOZEN M A，NIELSEN N N. Reinforced concrete response to simulated earthquakes [J]. Journal of the Structural Division，1970，96（12）：2557—2573.

[6]? ? IBARRA L F，KRAWINKLER H. Global collapse of frame structures under seismic excitations [R]. California：Pacific Earthquake Engineering Research Center，2005.

[7]? ? IBARRA L F，MEDINA R A，KRAWINKLER H. Hysteretic models that incorporate strength and stiffness deterioration [J]. Earthquake Engineering & Structural Dynamics，2005，34（12）：1489—1511.

[8]? ? 王濤，吳斌. 基于約束UKF模型更新的混合試驗方法[J]. 地震工程與工程振動，2013，33（5）：100—109.

WANG T，WU B.Hybrid testing method based on model updating with constrained unscented Kalman filter[J]. Journal of Earthquake Engineering and Engineering Vibration，2013，33（5）：100—109. （In Chinese）

[9]? ? 譚興龍，王堅，李增科. 基于SVD的改進抗差UKF算法及在組合導航中的應用[J]. 控制與決策，2014，29（10）：1744—1750.

TAN X L，WANG J，LI Z K. SVD aided improved robust UKF algorithm and its application for integration navigation [J]. Control and Decision，2014，29（10）：1744—1750.（In Chinese）

[10]? 陸新征，葉列平，潘鵬，等. 鋼筋混凝土框架結構擬靜力倒塌試驗研究及數(shù)值模擬競賽Ⅱ：關鍵構件試驗[J]. 建筑結構，2012， 42（11）：23—26.

LU X Z，YE L P，PAN P，et al. Pseudo-static collapse experiments and numerical prediction competition of RC frame structure Ⅱ：key elements experiment[J]. Building Structure，2012，42（11）：23—26. （In Chinese）

[11]? 陸新征，葉列平，潘鵬，等. 鋼筋混凝土框架結構擬靜力倒塌試驗研究及數(shù)值模擬競賽Ⅰ：框架試驗[J]. 建筑結構，2012，42（11）：19—22.

LU X Z，YE L P，PAN P，et al. Pseudo-static collapse experiments and numerical prediction competition of RC frame structureⅠ：RC frame experiment[J]. Building Structure，2012，42（11）：19—22. （In Chinese）

[12]? HASELTON C B，LIEL S，TAYLOR L，et al. Beam-column element model calibrated for predicting flexural response leading to global collapse of RC frame buildings[R]. California：Pacific Earthquake Engineering Research Center，2007.