基于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)针对特定的应用场景，本文仅研究了单自由度折线型本构模型参数识别。对于多自由度结构系统将更加复杂，而且折线型非线性模型中待识别参数均具有明确的物理意义与取值范围，比如屈服位移、屈服比等参数不可能为负值，而在迭代过程中不能保证这些参数非负，因而会造成算法中断，需引入其他约束方法进行进一步研究。