统计容积卡尔曼滤波器的混合试验模型更新方法

王 涛， 李 勐， 孟丽岩， 许国山， 王 贞

(1.黑龙江科技大学 建筑工程学院，哈尔滨 150022；2.哈尔滨工业大学 土木工程学院，哈尔滨 150090；3.武汉理工大学 土木与建筑学院，武汉 430070)

混合试验方法是一种结合试验子结构试验加载及数值模拟的技术，是检验大型复杂结构动力灾变性能的有效试验手段。为提高混合试验数值模型精度，模型更新混合试验方法应运而生[1]，该方法将模型更新技术与混合试验方法相结合，利用试验子结构观测数据在线识别试验子结构模型参数，实时更新数值子结构模型参数，以提高数值模型精度。

迄今为止，模型更新混合试验的研究已获得较大进展[2-3]。Zhang等[4]采用最小二乘法识别了双折线模型及Bouc-Wen模型参数，并将其应用于混合试验中。Hashemi等[5]以OpenFresco为基础，完成了单层单跨框架结构的模型更新混合试验，采用无迹卡尔曼滤波器(unscented Kalman filter,UKF)识别Bouc-Wen模型参数，并在线更新了数值子结构模型参数。王涛等[6]提出UKF模型更新混合试验方法，在线识别单自由度Bouc-Wen模型的参数，并通过对弹簧试件的单自由度体系的实际试验验证了此试验方法的有效性。梅竹等[7-8]使用UKF算法对钢筋混凝土结构材料本构模型参数进行了在线参数识别。陈再现等[9]提出基于均匀设计的模型更新混合仿真方法，在试验子结构的基础上，通过均匀设计的手段构造验算子结构样本空间，得到与试验子结构最为吻合的数值模型进行模型参数更新，有效提高了识别效率。Atasaratnam等[10]提出容积卡尔曼滤波器(cubature Kalman filter，CKF)算法，CKF基于数值积分理论，其核心步骤是根据三阶球面径向容积规则通过系统的初始矩阵、初始状态和初始点集确定2n个容积点[11]，并使这些容积点接近具有高斯噪声非线性系统的状态和协方差[12]。与UKF相比，CKF避免线性化的原理不同[13]，CKF采用三阶球面径向容积准则逼近高斯积分，具有较高的精度和稳定性，其容积点和权重仅由状态维数决定，可以提前计算和存储[14-15]，是建立在贝叶斯估计理论框架下的次优滤波，其结果满足正态分布下的最小方差理论。王涛等[16]采用CKF针对Bouc-Wen模型在线参数识别并与UKF进行对比，研究表明，CKF在运行效率上，仿真耗时少，识别结果更接近真实值，并能有效降低模型参数识别过程中的波动幅度，提高强非线性模型在线参数识别精度。然而，在模型更新混合试验中，CKF算法受制于算法初始参数的选择，初始参数选择不当，往往会降低算法精度，有时导致混合试验失败，大大限制了模型更新混合试验的工程应用。

鉴于此，本文以CKF算法为例，提出基于统计CKF混合试验模型更新方法，解决由模型更新算法初始参数选择不当对模型参数识别结果的影响。以自复位摩擦耗能支撑模型为加载对象，应用该方法对模型进行参数识别，分析在不同参数条件下统计CKF算法的识别精度；以两层带有自复位摩擦耗能支撑框架结构进行模型更新混合试验数值仿真，检验该方法的鲁棒性。

1 统计CKF模型更新算法

应用CKF进行参数识别时，需人为假定算法初始参数，初始参数的选择对参数识别结果有很大影响。目前，没有方法确定合理初始参数值，为此，笔者提出统计CKF模型更新算法，以弱化模型初始参数选择不当对识别结果的影响，提高算法对初始参数的鲁棒性。

1.1 统计CKF模型更新算法步骤

首先，建立如下系统的状态方程和观测方程

Xk=f(Xk-1,uk-1)+Vk-1

(1)

Yk=h(Xk,uk)+Wk

(2)

式中：u为输入向量；V为过程噪声向量，假定过程噪声向量是均值为0的高斯白噪声，其协方差矩阵为R；W为观测噪声向量，假定观测噪声向量是均值为0的高斯白噪声，其协方差矩阵为Q；X为系统状态向量；Y为系统观测向量；k为算法运行步数，k=1,2,…N。

(3)

(4)

ζ为容积点集，其表达式为

(5)

[e]i=

(6)

式中，[e]i为第i个容积点。

(7)

(8)

(9)

(10)

重采样

(11)

(12)

(13)

(14)

(15)

(16)

校正状态量

(17)

式中,yk为第k步的试验观测值。

校正协方差阵

(18)

对算法运行M次的状态量识别值进行统计，将统计后的状态量识别值均值作为最终识别结果

(19)

图1 统计CKF模型更新算法流程图
Fig.1 Statistical CKF model update algorithm flow chart

2 自复位摩擦耗能支撑模型参数识别

2.1 经改进旗形模型表达式

为验证基于统计CKF算法的有效性，以自复位摩擦耗能支撑为研究对象，采用经改进的旗形模型(modified flag-shaped model,MFS)[17]模拟自复位摩擦耗能支撑的滞回性能，对MFS模型进行模型参数在线识别。MFS模型由线性部分、双线性弹性部分、具有可滑移弹塑性部分组成，自复位摩擦耗能支撑机理图如图2所示。初始阶段刚度K为模型线性、双线性、具有可滑移弹塑性三部分的刚度之和，即K=K1+K2+K3。当摩擦装置启动后，支撑刚度减小至预应力筋刚度K1，卸载初始阶段的卸载刚度等于K1+K3，当力到达激活点b时，刚度变成K1+K2。

(a) 线性部分

(b) 双线性部分

(c) 具有滑移区弹塑性部分

(d) 整体图2 自复位摩擦耗能支撑机理图Fig.2 Self-centering energy dissipation brace mechanics diagram

MFS模型恢复力F可表达为

F=K1x+K2R(x)+K3z

(20)

式中，R(x)为由双线性弹性模型提供的恢复力，kN。

R(x)具体表达式可由式(21)给出

R(x)=x[1-H(x-b)-H(-x-b)]+

b[H(x-b)-H(x-b)]

(21)

式中：H(x)为海维赛德阶跃函数；b为激活位移，mm；z为模型滞变位移，mm。

z微分表达式由式(22)给出

(22)

K1、K2、K3由式(23)～(25)给出

K1=AK

(23)

K2=(1-A)(1-Q)K

(24)

K3=(1-A)QK

(25)

式中，A为激活刚度比，其表达式为

A=(1-Q)b

(26)

Q为耗能率，是能量耗散的重要指标，其表达式为

Q=K3/(K2+K3)

(27)

本文设定MFS模型参数分别为：耗能率Q=0.9、激活位移b=1.4 mm、初始刚度K=300 kN/mm、激活刚度比A=0.9。

2.2 统计CKF算法初始参数

首先，假定状态量为MFS模型滞变位移z、加载位移x、及模型参数b、K、Q，观测量为MFS模型恢复力，并根据式(22)、(20)建立算法状态方程及观测方程，算法状态方程和观测方程如式(28)～(29)所示。

(28)

式中,x为加载位移。

(29)

在数值积分法得到的状态方程的离散形式中，加入过程噪声Vs，过程噪声Vs和观测噪声Ws分别为均值为零的高斯白噪声，过程噪声Vs和观测噪声Ws的协方差矩阵分别为Rs、Qs，具体表达方式如式(30)～(31)所示。

Rs=10-25×I5

(30)

式中,I5为5×5的单位矩阵。

Qs=10-8×I1

(31)

(32)

式中，dt为数值积分步长。

(33)

(34)

作用于模型的加载方式选用低周期往复位移加载，加载方式如图3所示。

图3 加载制度Fig.3 Loading protocol

2.3 不同初始参数下统计CKF方法识别MFS模型参数

CKF算法识别模型参数的可靠程度往往受制于算法初始参数影响，针对这一问题，本节以MFS模型为例，研究模型使用不同初始参数条件下的模型参数识别结果，不同初始参数条件包括：不同初始协方差矩阵、初始参数预估值以及观测噪声，以验证统计CKF模型更新算法的识别精度。

2.3.1 初始协方差矩阵

(a) 滞变位移识别值

(b) 加载位移识别值

(c) b识别值

(d)K识别值

(e) Q识别值

(f) 滞回曲线图4 不同初始协方差矩阵的参数识别值结果Fig.4 Results of parameter identification values of different initial covariance matrices

本次识别过程中，CKF单次识别平均耗时4.01 s，使用统计CKF算法总耗时119.52 s。由图4(a)～(e)可以看出，使用不同的初始协方差矩阵导致每次运行CKF算法识别参数的结果并不唯一，与模型真实参数值差异较大。从图4(f)可以看出单次使用CKF算法得到的滞回曲线不能反映模型的真实性能，甚至没有意义。采用统计CKF算法识别模型参数精度接近真实值，收敛性较好，参数识别效果得到改善，恢复力预测精度明显提高。此外，通过耗能的角度分析，单次使用CKF算法得到的模型滞回曲线最大耗能值与真实耗能值相对误差为86.26%，统计CKF算法得到的滞回曲线耗能值与真实耗能值相对误差15.41%，耗能值相对误差结果降低70.85%；单次使用CKF算法得到的模型滞回曲线最小耗能值与真实耗能值的相对误差为10.12%，相比使用统计方法得到的耗能值相对误差降低了5.29%。虽然单次使用CKF算法得到的耗能值相对误差更低，但是考虑到CKF是一种随机算法，每次应用算法得到的结果均存在随机性，相比滞回曲线最大耗能值相对误差，统计得到的模型滞回曲线有明显改善，统计后的结果更具有可靠性，通过统计生成的滞回曲线与真实值差异较小，可以证实所提出方法的有效性。为了能定量评价统计CKF算法识别精度，定义均方根误差(RMSE)为

(35)

式中：xref为模型参数真实值；xsim为模型参数识别值。

使用不同初始协方差矩阵，模型参数识别值的均方根误差如图5。

图5 参数识别值均方根误差图Fig.5 Root mean square error diagram of parameter identification values

通过均方根误差分析可以看出，相比CKF算法识别模型参数值最大均方根误差，统计CKF算法得到的模型参数识别值均方根误差分别降低了18.51%、73.78%、237.64%。初始误差协方差矩阵不同，单次使用CKF算法得到的模型参数识别结果与真实值均方根误差较高，统计CKF算法识别模型参数识别值结果均方根误差更低，而且改善效果明显，更能有效的反映结构加载后的真实性能，具有较强的鲁棒性。不同统计次数下，得到的模型滞回曲线结果如图6所示。

(a) 滞回曲线

(b) 相对误差图6 不同统计次数得到的模型滞回曲线及相对误差图Fig.6 The hysteretic curves and relative error graphs of different statistical times

为了检验统计次数对模型更新精度影响，分别对比统计10～50次得到的模型滞回曲线，如图6(a)所示。由图6(a)可以看出，随着统计次数的增加，统计50次得到的模型滞回曲线较统计10～40次得到的模型滞回曲线更接近模型真实滞回形状。由图6(b)可以看出，从耗能的角度分析，统计50次得到的模型滞回曲线耗能值与真实耗能值相对误差为11.56%，相比统计10～40次得到的滞回曲线耗能值相对误差结果分别降低了32.67%、25.63%、3.85%、2.32%，可以发现当统计次数为30次时，改善效果最明显。充分说明统计方法得到的参数识别结果可以提高算法鲁棒性，优化算法随机误差。

2.3.2 初始参数预估值

本次识别过程中，CKF单次识别平均耗时3.95 s，使用统计CKF算法总耗时118.24 s。从图7(a)～(e)可以看出，当初始参数预估值随机时，每次运行CKF算法得到的参数识别结果与真实值具有很大差异，统计CKF算法识别模型参数结果可以较快地收敛至真实值。从耗能的角度看，图7(f)单次使用CKF识别滞回曲线最大耗能值与真实耗能值最大相对误差为132.37%，而统计过后的滞回曲线耗能值与真实耗能值相对误差14.76%，耗能值相对误差结果降低117.61%。单次使用CKF算法得到得模型滞回曲线最小耗能值与真实耗能值相对误差为10.45%，相比使用统计CKF算法降低4.31%。使用不同初始参数预估值，模型参数识别值的均方根误差如图8所示。

(a) 滞变位移识别值

(b) 加载位移识别值

(c) b识别值

(d) K识别值

(e) Q识别值

(f) 滞回曲线图7 不同初始参数预估值的参数识别值结果Fig.7 Results of parameter identification values estimated by different initial parameters

图8 参数识别值均方根误差图Fig.8 Root mean square error diagram of parameter identification values

通过误差分析可以得出，当使用不同初始参数预估值时，相比单次使用CKF算法识别模型参数最大均方根误差，统计CKF算法识别模型参数值均方根误差分别降低了69.68%、107.09%、244.91%，有效提高了算法识别精度。不同统计次数下，得到的模型滞回曲线结果如图9所示。

为了检验统计次数对模型更新精度影响，分别对比统计10～50次模型滞回曲线，如图9(a)所示。由图9(a)可以看出，随着统计次数的增加，统计50次得到的模型滞回曲线较统计10～40次得到的模型滞回曲线更接近模型真实滞回形状。由图9(b)可以看出，从耗能的角度分析，统计50次得到的模型滞回曲线耗能值与真实耗能值相对误差为12.44%，相比统计10～40次得到的滞回曲线耗能值相对误差结果分别降低了73.29%、28.55%、2.32%、1.77%。充分证明了统计方法得到的模型滞回曲线更接近真实值，可以更好地模拟结构的真实反映。

(a) 滞回曲线

(b) 相对误差图9 不同统计次数得到的模型滞回曲线及相对误差图Fig.9 The hysteretic curves and relative error graphs of different statistical times

2.3.3 观测噪声

观测噪声对模型参数的识别结果具有很大的影响，这种影响有时会使信号变形、失真，有的严重导致数据不可用。为了研究不同观测噪声对CKF算法识别模型参数的影响，本次模拟选用不同观测噪声对MFS模型参数进行识别，观测噪声Ws选用仿真观测噪声乘随机系数，运行算法过程中算法初始参数预估值、初始误差协方差矩阵不变，参数识别结果如图10所示，其中单次识别值为单次使用CKF算法识别结果，统计识别值为CKF算法运行30次得到的参数识别值均值。

本次识别过程中，CKF单次识别平均耗时3.87 s，使用统计CKF算法总耗时120.35 s。在图10中可以发现，使用不同观测噪声时，每次运行CKF算法识别模型参数的结果并不唯一，较模型参数真实值存在很大差异。从耗能的角度出发，单次使用CKF算法识别模

(a) 滞变位移识别值

(b) 加载位移识别值

(c) b识别值

(d) K识别值

(e) Q识别值

(f) 滞回曲线图10 不同观测噪声的参数识别值结果Fig.10 Results of parameter identification values of different observation noises

型滞回曲线最大耗能与真实耗能值相对误差为77.96%，而统计得到的滞回曲线耗能与真实耗能值相对误差16.52%，耗能值相对误差结果降低61.44%，单次使用CKF算法得到得模型滞回曲线最小耗能值与真实耗能值相对误差为14.72%，相比使用统计CKF算法降低1.8%。经过统计后的参数识别值与真实值相比明显误差更小，吻合程度更好，稳定性高。使用不同观测噪声，模型参数识别值的均方根误差如图11所示。

图11 参数识别值均方根误差图Fig.11 Root mean square error diagram of parameter identification values

通过误差分析可以得出，相比单次使用CKF算法识别模型参数最大均方根误差，统计CKF算法均方根误差分别降低了18.39%、14.03%、136.05%。充分证明了当观测噪声不同时，统计CKF方法识别参数可以有效降低观测噪声对试验的影响，从而保证了试验的鲁棒性。不同统计次数下，得到的模型滞回曲线如图12所示。

为了检验统计次数对模型更新精度影响，分别对比统计10～50次模型滞回曲线，如图12(a)所示。由图12(a)可以看出，随着统计次数的增加，统计50次得到的模型滞回曲线比统计10～40次得到的模型滞回曲线更能反映模型的真实滞回现象。由图12(b)可以看出，从耗能的角度分析，统计50次得到的模型滞回曲线耗能值与真实耗能值相对误差为14.26%，相比统计10～40次得到的滞回曲线耗能值相对误差结果分别降低了53.22%、27.54%、2.16%、1.36%。证明统计方法得到的模型滞回曲线与真实滞回曲线更为接近，可以更好地模拟结构的真实反映。

(a) 滞回曲线

(b) 相对误差图12 不同统计次数得到的模型滞回曲线及相对误差图Fig.12 The hysteretic curves and relative error graphs of different statistical times

可以发现，统计30次较统计40次、50次相对误差较小，可以得到相对较满意的结果，因此后文数值模拟部分统计次数均采用统计30次的运行结果，以验证统计方法的有效性。

3 基于统计CKF模型参数在线更新混合试验数值模拟

3.1 数值模拟原理

为了进一步检验统计CKF方法模型更新混合试验的精度及鲁棒性，针对两层带有自复位摩擦耗能支撑结构采用MATLAB2018b编程进行混合试验数值模拟。首先，选用下层自复位摩擦耗能支撑作为试验子结构，上层自复位摩擦耗能支撑作为数值子结构，通过数值

图13 基于统计CKF模型更新混合模拟方法示意图Fig.13 Schematic diagram of hybrid simulation method based on statistical CKF model update

3.2 参数设置

设定每层支撑恢复力模型均采用本文2.1节MFS模型。试验过程中假定主体结构处于线性状态，框架结构模型参数设置为：每层结构质量Mn1=Mn2=2 000 t，框架结构的层间水平刚度刚度Kn1=Kn2=80 000 kN/m，框架结构的层间阻尼系数Cn1=Cn2=1 550 kN/(m/s)，其中n1、n2分别代表了数值子结构的一层和二层，支撑与楼面夹角均为28.81°。MFS真实模型参数分别为：激活位移b=0.001 4 m、初始刚度K=300 000 kN/m、耗能率Q=0.9、激活刚度比A=0.015。地震作用选取EI Centro(1940,NS)地震记录，地震动峰值加速度为400 cm/s2，地震记录如图14所示。支撑滞变位移利用四阶Runge-Kutta积分算法得出，积分步长为0.01 s。

图14 EI Centro地震记录Fig.14 EI Centro seismic record

假定状态量为MFS模型滞变位移z、加载位移x、及模型参数b、K、Q，观测量为MFS模型恢复力，利用识别关键参数值去更新数值子结构支撑部分参数。试验子结构的状态方程和观测方程同式(28)～(29)所示，试验结构加载速度采用式(32)位移差分方法确定。过程噪声Vs和观测噪声Ws分别为均值为零的高斯白噪声，过程噪声Vs和观测噪声Ws的协方差矩阵分别为Rs、Qs，具体表达如式(36)～(37)所示。

Rs=10-26×I5

(36)

Qs=10-10×I1

(37)

(38)

(39)

3.3 自复位摩擦耗能支撑框架混合试验数值模拟

(a) 下层支撑滞变位移识别值

(b) 下层支撑加载位移识别值

(c) b识别值

(d) K识别值

(e) Q识别值

(f) 下层支撑滞回曲线图15 MFS模型参数识别值Fig.15 MFS model parameter identification values

上层支撑的滞回曲线数值模拟结果如图16所示。

图16 上层支撑滞回曲线数值模拟结果Fig.16 Numerical simulation results of upper bracing hysteretic curve

本次模拟过程中，CKF单次识别平均耗时6.71 s，使用统计CKF算法总耗时180.92 s。从图15～16可以看出，单次使用CKF方法得到的模型参数识别值、模型滞回曲线与真实值有明显差异，基于统计CKF方法得到的模型参数识别值、模型滞回曲线较真实值吻合程度较高。从耗能的角度分析，单次使用CKF方法得到的下层支撑模型滞回曲线最大耗能值较真实耗能值相对误差18.2%，统计CKF方法得到的下层支撑模型滞回曲线耗能值与真实耗能值相对误差为1.83%，耗能值相对误差降低16.37%。单次使用CKF方法得到的上层支撑模型滞回曲线耗能值较真实耗能值相对误差为48.91%，统计CKF方法得到的上层支撑模型滞回曲线耗能值较真实耗能值相对误差为22.26%，耗能值相对误差降低了26.65%。单次使用CKF算法得到得模型滞回曲线最小耗能值与真实耗能值相对误差为18.55%，相比使用统计CKF算法降低3.71%。模型参数识别值的均方根误差如图17所示。

通过均方根误差分析，相比单次使用CKF方法识别模型参数最大均方根误差，统计CKF方法均方根误差分别降低了30.84%、43.22%、25.08%，混合试验模型参数识别精度均有很大改善。可见，统计CKF方法可以明显改善参数识别结果，更具有鲁棒性。

4 结 论

本文提出了一种基于统计CKF模型参数更新方法，通过对MFS模型进行不同初始参数条件下的参数识别，并针对两层带有自复位摩擦耗能支撑结构进行混合试验数值仿真，可以得到以下结论：

(1) 在不同初始参数情况下，采用统计CKF算法对MFS模型进行参数识别，参数识别结果及滞回曲线与MFS模型真实值吻合程度更高，具有较高的识别精度。

(2) 采用统计CKF方法对两层带有自复位摩擦耗能支撑结构进行混合试验数值仿真，可以弱化初始参数对混合试验识别结果的影响，提高混合试验模型参数在线识别精度及鲁棒性。