CN105466661A - 基于改进卡尔曼滤波的超高层建筑风荷载反分析方法 - Google Patents

Abstract

本发明公开了一种基于改进卡尔曼滤波的超高层建筑风荷载反分析方法，该方法利用有限测试楼层的风致响应识别结构动态风荷载，所属于建筑结构风荷载反分析领域，通过在模态空间基于动力微分方程构造离散***状态方程，基于改进型的卡尔曼滤波估计响应和荷载，进而改进风荷载反演技术。本发明方法可以解决超高层建筑风致响应测点不足的问题，并且在结构模态参数误差、模态截断及测量噪声等因素影响下的风荷载识别精度能够满足实际工程需要，本发明方法为超高层建筑抗风设计及相关研究提供了有用的工具及依据。

Description

基于改进卡尔曼滤波的超高层建筑风荷载反分析方法

技术领域

本发明涉及建筑结构风荷载识别技术领域，尤其涉及一种基于改进卡尔曼滤波的超高层建筑风荷载反分析方法。

背景技术

近年来在我国东南沿海地区兴建了大量三百米以上的超高层建筑,由于这类结构的自振频率较低、阻尼较小，同台风动荷载的主要频率段比较接近,在强/台风作用下的风致响应比较大，风荷载及风振响应是其安全性及适用性设计的首要控制指标。一般来讲，高层建筑风荷载很难由现场测量确定，风洞试验技术尽管能够测试结构风荷载，但由于受到试验模拟技术及分析方法的限制，其评估结果难以完全准确的反映真实情况。考虑到目前结构动力响应的实测技术比较成熟，结构加速度及位移响应的量测精度远高于风荷载的量测精度。因此，基于实测风致响应反向识别结构动态风荷载，成为间接量测风荷载的一种新途径。

利用结构响应评估外荷载属于典型的反分析问题，Sanchez和Benaroya(2014)对该问题的研究进展进行了详细地介绍，并讨论了多种荷载识别技术的适用性及局限性。而对于结构风荷载的反分析研究也取得了一些成果。如陈隽及李杰(1997)基于荷载归一化统计平均法开展了高层建筑风荷载的反演分析。Kang和Lo(2002)对典型高塔进行了风荷载反演分析。Law等人(2005)提出了一种评估风荷载的反分析算法，并以桅杆为对象进行了数值验证。Lourens等人(2008)基于一种动力反分析算子识别了80m高塔的时变风荷载。Hwang及Kareem等(2009,2010)基于混凝土烟囱的气弹风洞试验评估了结构横风向荷载，并对比研究了多种因素对烟囱风荷载分析结果的影响。以往的研究表明，基于结构风致响应识别动态风荷载是一种可行且有效的方法，它对进一步理解风与结构的相互作用机理及风振响应规律具有重要的意义。但目前来说对建筑结构风荷载的反分析研究较少，对于高层建筑风荷载反演分析更是有限，因此很有必要加强这方面的研究工作。

发明内容

本发明要解决的技术问题在于针对现有技术中的缺陷，提供一种基于改进卡尔曼滤波的超高层建筑风荷载反分析方法。

本发明解决其技术问题所采用的技术方案是：一种基于改进卡尔曼滤波的超高层建筑风荷载反分析方法，包括以下步骤：

1)由有限元计算或质量统计可得到结构质量矩阵M，其中超高层建筑的层数为n；

基于现场实测获得结构前q阶结构自振频率ω_i、阻尼比ξ_i，综合应用现场实测和有限元分析获得结构模态振型Φ_n×q；

2)输入实测的p个楼层的风致响应分量，根据实测的结构前q阶模态振型，将测试的风致响应分量转化为模态风致响应；所述风致响应分量为位移、速度或加速度响应中的一种；

在结构动态响应实测时，若只测得结构部分楼层的风致响应(假设为p层加速度响应)和前q阶模态振型，由于超高层建筑的风致振动往往以前几阶模态为主，则结构风致响应可近似分解为：

${\stackrel{\·\·}{y}}_{p\×1}={\mathrm{\Φ}}_{p\×q}{\stackrel{\·\·}{U}}_{q\×1},(1\≤q\≤p\≤n)$

其中，为实测的p层加速度响应；为前q阶模态加速度向量；Φ_p×q由Φ_n×n中与p个实测响应楼层对应的行、前q列形成的子振型矩阵；Φ_n×n为按质量规准化的模态振型矩阵；确定结构振动的主要控制模态数目q采用以下方法：基于POD方法首先获取加速度响应协方差矩阵的特征值进而计算出前q阶模态对结构振动的贡献比例：

$\mathrm{\θ}=\frac{\underset{i=1}{\overset{q}{\Σ}}{\mathrm{\λ}}_i}{\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\λ}}_i}(1\≤q\≤n)$

取θ超过预设值时所对应的q值作为结构振动主要控制模态数目。

由广义逆矩阵Φ_p×q ⁺，结构实测的模态加速度响应可近似表示为：

${\widehat{\stackrel{\·\·}{U}}}_{q\×1}={\left({\mathrm{\Φ}}_{p\×q}\right)}^{+}{\stackrel{\·\·}{y}}_{p\×1}$

式中：(Φ_p×q)⁺为Φ_p×q的广义逆矩阵；

3)根据结构动力微分方程，在模态空间中构造离散化的状态方程和观测方程；

风荷载作用下有n个楼层的超高层建筑的动力微分方程可表示如下：

$M\stackrel{\·\·}{y}+C\stackrel{\·}{y}+Ky=F$

式中y、和是位移、速度和加速度向量；F为外荷载；C和K为阻尼矩阵和刚度矩阵；

动力方程可按如下解耦：

${\mathrm{\Φ}}_i^{T}F=f_i=M_i{\stackrel{\·\·}{U}}_i+C_i{\stackrel{\·}{U}}_i+K_i{U}_i$

式中：M_i、K_i分别为第i阶按质量规准化的模态质量和模态刚度，其中M_i＝1(i＝1,2,…,q),；f_i、C_i分别为第i阶按质量规准化的模态荷载，其中C_i＝2ξ_iω_i；U_i分别为第i阶模态加速度、速度、位移。将上式转换为状态空间形式：

${\stackrel{\·}{X}}_i\left(t\right)=A_i{X}_i\left(t\right)+B_i{f}_i$

其中:

${\begin{array}{ccc}{X}_i\left(t\right)={\left[\begin{array}{cc}{U}_i& {\stackrel{\·}{U}}_i\end{array}\right]}^T,& A_i=\left[\begin{array}{cc}0& 1\\ -K_i& -C_i\end{array}\right],& B_i=\left[\begin{array}{cc}0& 1\end{array}\right]\end{array}}^T$

将上式离散化，得到离散***方程：

X_i(k+1)＝Ψ_i(k+1/k)X_i(k)+Γ_i(k+1/k)f_i(k)

其中：Ψ_i(k+1/k)＝exp(A_iΔt)

${\mathrm{\Γ}}_i(k+1/k)={\∫}_{k\mathrm{\Δ}t}^{(k+1)\mathrm{\Δ}t}\exp \{A_i\[(k+1)\mathrm{\Δ}t-\mathrm{\τ}\]\}{B}_{i}d\mathrm{\τ}$

式中：X_i(k)和f_i(k)分别kΔt时刻的状态向量及模态荷载；Δt为采样间隔；Ψ_i(k+1/k)为kΔt时刻至(k+1)Δt时刻的一步转移矩阵；Γ_i(k+1/k)为***噪声驱动矩阵。

由Ψ_i(k+1/k)和Γ_i(k+1/k)表达式可知，当Δt取值一定时，Ψ_i(k+1/k)和Γ_i(k+1/k)都为常量矩阵(即不随时间变化)，所以离散***方程可以简写为：

X_i(k+1)＝Ψ_iX_i(k)+Γ_if_i(k)

上式即为基于结构动力微分方程构造的离散型***状态方程。

由式 ${\stackrel{\·}{X}}_i\left(t\right)=A_i{X}_i\left(t\right)+B_i{f}_i$ 可知:

$X_i\left(k\right)={\left[\begin{array}{cc}{U}_i\left(k\right)& {\stackrel{\·}{U}}_i\left(k\right)\end{array}\right]}^T$

式中：U_i(k)和分别kΔt时刻的模态位移及速度。因此，***观测方程可以写为如下的形式：

Z_i(k)＝H_iX_i(k)+D_if_i(k)+ε_i(k)

其中：Z_i(k)为kΔt时刻的响应观测值(为脉动响应)；D_i为***矩阵；ε_i(k)为观测噪声；H_i为观测矩阵，随着输入响应类型不同而发生变化。

当响应观测值Z_i(k)的类型为位移时：

H_i＝[10],D_i＝0

当响应观测值Z_i(k)的类型为速度时：

H_i＝[01],D_i＝0

当响应观测值Z_i(k)的类型为加速度时：

H_i＝[-K_i-C_i],D_i＝1

关于模态荷载和观测噪声可做如下假设：

E[f_i(k)]＝0

E[f_i(k)f_i ^T(j)]＝Q_i(k)δ_kj

E[ε_i(k)]＝0

E[ε_i(k)ε_i ^T(j)]＝R_i(k)δ_kj

式中:E[·]为求随机变量的期望；Q_i(k)是f_i(k)的方差强度矩阵，为对称非负定矩阵；R_i(k)是ε_i(k)的方差强度矩阵，为对称正定矩阵；δ_kj是Kronecker-δ函数。***的过程噪声f_i(k)和观测噪声ε_i(j)不相关，即：

E[f_i(k)ε_i ^T(j)]＝0

为了构造符合卡尔曼滤波经典理论的离散型观测方程，可令：

V_i(k)＝D_if_i(k)+ε_i(k)

则***观测方程Z_i(k)＝H_iX_i(k)+D_if_i(k)+ε_i(k)最终可以写为如下的形式：

Z_i(k)＝H_iX_i(k)+V_i(k)

则由f_i(k)和ε_i(k)的性质，可得到关于V_i(k)的如下关系：

E[V_i(k)]＝E[D_if_i(k)+ε_i(k)]

＝D_iE[f_i(k)]+E[ε_i(k)]；

＝0

$\begin{array}{c}E\[V_i\left(k\right){V_i}^T\left(j\right)\]=E\[D_i{f}_i\left(k\right){f_i}^T\left(k\right){D_i}^T+D_i{f}_i\left(k\right){{\mathrm{\ϵ}}_i}^T\left(k\right)+{\mathrm{\ϵ}}_i\left(k\right){f_i}^T\left(k\right){D_i}^T+{\mathrm{\ϵ}}_i\left(k\right){{\mathrm{\ϵ}}_i}^T\left(j\right)\]\\ =\[D_i{Q}_i\left(k\right)D_i^T+R_i\left(k\right)\]{\mathrm{\δ}}_{kj}\\ =r_i\left(k\right){\mathrm{\δ}}_{kj}\end{array}$

由上式可知：V_i(k)为零均值的白噪声随机过程，r_i(k)是V_i(k)的方差强度矩阵，

而： $\begin{array}{c}E\[f_i\left(k\right){V_i}^T\left(j\right)\]=E\[f_i\left(k\right){f_i}^T\left(j\right){D_i}^T+f_i\left(k\right){{\mathrm{\ϵ}}_i}^T\left(j\right)\]\\ =Q_i\left(k\right)D_i^T{\mathrm{\δ}}_{kj}\\ =S_i\left(k\right){\mathrm{\δ}}_{kj}\end{array}$

式中：S_i(k)是f_i(k)和V_i(k)的协方差强度矩阵，且

综上所诉：f_i(k)和V_i(k)都为零均值的白噪声随机过程，且f_i(k)和V_i(k)相关。因此***状态方程及观测方程满足白噪声相关条件下经典的离散型卡尔曼滤波方程。

4)对应不同的响应类型，基于改进卡尔曼滤波理论，利用已测部分楼层的风致响应估计结构未知风致响应分量；

将状态方程、观测方程及相关假设代入卡尔曼滤波基本方程，即可得到kΔt时刻的X_i(k)的估计

${\hat{X}}_i\left(k\right):$

${\hat{X}}_i(k/k-1)={\mathrm{\ψ}}_i{\hat{X}}_i(k-1)+J_i(k-1)\[Z_i(k-1)-H_i{\hat{X}}_i(k-1)\]$

$\hat{K}\left(k\right)={\hat{X}}_i(k/k-1)+G_i\left(k\right)\[Z_i\left(k\right)-H_i{\hat{X}}_i(k/k-1)\]$

J_i(k-1)＝Γ_iQ_i(k-1)D_i ^T[D_iQ_i(k-1)D_i ^T+R_i(k-1)]^-1

P_i(k/k-1)＝[Ψ_i-J_i(k-1)H_i]P_i(k-1)[Ψ_i-J_i(k-1)H_i]^T+

Γ_iQ_i(k-1)Γ_i ^T-J_i(k-1)D_iQ_i(k-1)Γ_i ^T

$G_i\left(k\right)=P_i(k/k-1)H_i^T{\[H_i{P}_i(k/k-1)H_i^T+D_i{Q}_i\left(k\right){D_i}^T+R_i\left(k\right)\]}^{-1}$

P_i(k)＝[I-G_i(k)H_i]P_i(k/k-1)

式中G_i(k)为最优卡尔曼滤波增益；***状态向量的估计；J_i(k)是状态一步预测增益矩阵，P_i(k/k-1)是一步预测误差方差矩阵，P_i(k)是估计误差方差矩阵。

***初值选取如下：

${\hat{X}}_i\left(0\right)=E\[X_i\left(0\right)\]$

$P_i\left(0\right)=E\{\[X_i\left(0\right)-{\hat{X}}_i\left(0\right)\]{\[X_i\left(0\right)-{\hat{X}}_i\left(0\right)\]}^T\}$

在荷载估计之前外部荷载及测量噪声是未知的，可以先假定荷载协方差矩阵Q_i(k)为单位矩阵。通过改变观测噪声协方差矩阵R_i(k)的值(一般可取10^-4-10^-8)，可以估计最优卡尔曼滤波增益G_i(k)。

根据上诉方法，只要给定初值和P_i(0)，根据kΔt时刻的量测值Z_i(k)，就可递推得到kΔt时刻的***状态向量的估计 ${\hat{X}}_i\left(k\right),\left(k=1,2,\mathrm{...}\right).$

5)根据预测的模态响应，估计模态风荷载，进而获得结构任意楼层的风荷载时程；

具体如下：根据步骤4)中求得的***状态向量估计将其离散型***状态方程即可得：

${\mathrm{\Γ}}_i{\hat{f}}_i\left(k\right)={\hat{X}}_i(k+1)-{\mathrm{\Ψ}}_i{\hat{X}}_i\left(k\right)$

求得:

$\hat{f}\left(k\right)={{\mathrm{\Γ}}_i}^{+}\[{\hat{X}}_i(k+1)-{\mathrm{\ψ}}_i{\hat{X}}_i\left(k\right)\]$

式中:为Γ_i的广义逆。

基于上述方法依次可求得结构前q阶模态荷载估计将前q阶估计模态荷载组成向量:

${\hat{f}}_{q\×1}={\left[\begin{array}{cccc}{\hat{f}}_1& {\hat{f}}_2& \mathrm{...}& {\hat{f}}_q\end{array}\right]}^T$

当超高层建筑的风振分析只考虑前q阶模态时，即可求得结构脉动风荷载向量的估计值

式中为前q列对应的子矩阵；由于模态坐标转换理论可知，振型矩阵关于质量矩阵正交，即：

(Φ_n×n)^TMΦ_n×n＝I

式中I为n×n维单位矩阵，则可由下式求出：

本发明产生的有益效果是：本发明提出的风荷载反向识别新技术，能够利用超高层建筑有限测试楼层的部分风致响应分量(特别是可利用加速度响应分量)以及结构前若干阶主要控制模态，准确的识别出结构任意楼层的脉动风荷载及未测风致响应。该方法计算收敛速度快，抗噪声能力强，识别结果对结构参数误差及模态截断误差的敏感性较小。本发明对进一步理解风与结构的相互作用机理及改进现有的风荷载理论模型具有重要的意义。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1为平均风速剖面及湍流剖面；

图2为参考坐标轴；

图3为0°风向下第55层X向位移响应时程；

图4为0°风向下第55层X向速度响应时程；

图5为0°风向下结构55层X向位移响应功率谱；

图6为0°风向下结构55层X向速度响应功率谱；

图7为90°风向下结构35层X向风荷载时程；

图8为90°风向下结构35层Y向风荷载时程；

图9为90°风向下结构35层X向风荷载功率谱；

图10为90°风向下结构35层Y向风荷载功率谱；

图11为0°风向下结构底部X向总风力时程；

图12为0°风向下结构底部Y向总风力时程；

图13为0°风向角下结构底部X向总风力功率谱

图14为0°风向角下结构底部Y向总风力功率谱；

图15为0°风向下不同响应类型反演基底X向总风力功率谱；

图16为0°风向下不同响应类型反演基底Y向总风力功率谱；

图17自振频率误差10％时，加速度反演X向基底总风力功率谱(0°风向)；

图18自振频率误差10％时，位移反演X向基底总风力功率谱(0°风向)；

图19阻尼比误差10％时，加速度反演X向基底总风力功率谱(0°风向)；

图20阻尼比误差10％时，位移反演X向基底总风力功率谱(0°风向)；

图21为不同噪声水平下的加速度反演X向荷载功率谱(0°风向)；

图22为不同噪声水平下的加速度反演Y向荷载功率谱(0°风向)；

图23为本发明的方法流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

如图23所示，基于改进卡尔曼滤波的超高层建筑脉动风荷载反分析方法，包括以下步骤：

1)由有限元计算或质量统计可得到结构质量矩阵M，其中超高层建筑的层数为n；

基于现场实测获得结构前q阶结构自振频率ω_i、阻尼比ξ_i，综合应用现场实测和有限元分析获得结构模态振型Φ_n×q；

2)输入实测的p个楼层的风致响应分量，根据实测的结构前q阶模态振型，将测试的风致响应分量转化为模态风致响应；所述风致响应分量为位移、速度或加速度响应中的一种；

在结构动态响应实测时，若只测得结构部分楼层的风致响应(假设为p层加速度响应)和前q阶模态振型，由于超高层建筑的风致振动往往以前几阶模态为主，则结构风致响应可近似分解为：

${\stackrel{\·\·}{y}}_{p\×1}={\mathrm{\Φ}}_{p\×q}{\stackrel{\·\·}{U}}_{q\×1},(1\≤q\≤p\≤n)---\left(1\right)$

式中：为实测的p层加速度响应；为前q阶模态加速度向量；Φ_p×q由Φ_n×n中与p个实测响应楼层对应的行、前q列形成的子振型矩阵；Φ_n×n为按质量规准化的模态振型矩阵；确定结构振动的主要控制模态数目q采用以下方法：基于POD方法首先获取加速度响应协方差矩阵的特征值λ_i(i＝1,2,…n)，进而计算出前q阶模态对结构振动的贡献比例：

$\mathrm{\θ}=\frac{\underset{i=1}{\overset{q}{\Σ}}{\mathrm{\λ}}_i}{\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\λ}}_i}(1\≤q\≤n)---\left(2\right)$

取θ超过预设值时所对应的q值作为结构振动主要控制模态数目。

由广义逆矩阵Φ_p×q ⁺，结构实测的模态加速度响应可近似表示为：

${\widehat{\stackrel{\·\·}{U}}}_{q\×1}={\left({\mathrm{\Φ}}_{p\×q}\right)}^{+}{\stackrel{\·\·}{y}}_{p\×1}---\left(3\right)$

式中：(Φ_p×q)⁺为Φ_p×q的广义逆矩阵。

3)根据结构动力微分方程，在模态空间中构造离散化的状态方程和测量方程；

风荷载作用下有n个楼层的超高层建筑的动力微分方程可表示如下：

$M\stackrel{\·\·}{y}+C\stackrel{\·}{y}+Ky=F---\left(4\right)$

式中y、和是位移、速度和加速度向量；F为外荷载；C和K为阻尼矩阵和刚度矩阵。

动力方程(4)可按如下解耦：

${\mathrm{\Φ}}_i^{T}F=f_i=M_i{\stackrel{\·\·}{U}}_i+C_i{\stackrel{\·}{U}}_i+K_i{U}_i---\left(5\right)$

式中：M_i、K_i分别为第i阶按质量规准化的模态质量和模态刚度，其中M_i＝1(i＝1,2,…,q),f_i、C_i分别为第i阶按质量规准化的模态荷载，其中C_i＝2ξ_iω_i；U_i分别为第i阶模态加速度、速度、位移。将上式转换为状态空间形式：

${\stackrel{\·}{X}}_i\left(t\right)=A_i{X}_i\left(t\right)+B_i{f}_i---\left(6\right)$

其中:

${\begin{array}{ccc}{X}_i\left(t\right)={\left[\begin{array}{cc}{U}_i& {\stackrel{\·}{U}}_i\end{array}\right]}^T,& A_i=\left[\begin{array}{cc}0& 1\\ -K_i& -C_i\end{array}\right],& B_i=\left[\begin{array}{cc}0& 1\end{array}\right]\end{array}}^T$

将上式离散化，得到离散***方程：

X_i(k+1)＝Ψ_i(k+1/k)X_i(k)+Γ_i(k+1/k)f_i(k)(7)

其中：

Ψ_i(k+1/k)＝exp(A_iΔt)(8a)

${\mathrm{\Γ}}_i(k+1/k)={\∫}_{k\mathrm{\Δ}t}^{(k+1)\mathrm{\Δ}t}\exp \{A_i\[(k+1)\mathrm{\Δ}t-\mathrm{\τ}\]\}{B}_{i}d\mathrm{\τ}---\left(8b\right)$

式中：X_i(k)和f_i(k)分别kΔt时刻的状态向量及模态荷载；Δt为采样间隔；Ψ_i(k+1/k)为kΔt时刻至(k+1)Δt时刻的一步转移矩阵；Γ_i(k+1/k)为***噪声驱动矩阵。

由式(8)可知，当Δt取值一定时，Ψ_i(k+1/k)和Γ_i(k+1/k)都为常量矩阵(即不随时间变化)，所以式(7)可以简写为：

X_i(k+1)＝Ψ_iX_i(k)+Γ_if_i(k)(9)

上式即为基于结构动力微分方程构造的离散型***状态方程。

由式(6)可知:

$X_i\left(k\right)={\left[\begin{array}{cc}{U}_i\left(k\right)& {\stackrel{\·}{U}}_i\left(k\right)\end{array}\right]}^T---\left(10\right)$

式中：U_i(k)和分别kΔt时刻的模态位移及速度。因此，***观测方程可以写为如下的形式：

Z_i(k)＝H_iX_i(k)+D_if_i(k)+ε_i(k)(11)

其中：Z_i(k)为kΔt时刻的响应观测值(为脉动响应)；D_i为***矩阵；ε_i(k)为观测噪声；H_i为观测矩阵，随着输入响应类型不同而发生变化。

当响应观测值Z_i(k)的类型为位移时：

H_i＝[10],D_i＝0(12a)

当响应观测值Z_i(k)的类型为速度时：

H_i＝[01],D_i＝0(12b)

当响应观测值Z_i(k)的类型为加速度时：

H_i＝[-K_i-C_i],D_i＝1(12c)

关于模态荷载和测量噪声可做如下假设：

E[f_i(k)]＝0(13a)

E[f_i(k)f_i ^T(j)]＝Q_i(k)δ_kj(13b)

E[ε_i(k)]＝0(14a)

E[ε_i(k)ε_i ^T(j)]＝R_i(k)δ_kj(14b)

式中:E[·]为求随机变量的期望；Q_i(k)是f_i(k)的方差强度矩阵，为对称非负定矩阵；R_i(k)是ε_i(k)的方差强度矩阵，为对称正定矩阵；δ_kj是Kronecker-δ函数。***的过程噪声f_i(k)和观测噪声ε_i(j)不相关，即：

E[f_i(k)ε_i ^T(j)]＝0(15)

为了构造符合卡尔曼滤波经典理论的离散型观测方程，可令：

V_i(k)＝D_if_i(k)+ε_i(k)(16)

则***观测方程(11)可以写为如下的形式：

Z_i(k)＝H_iX_i(k)+V_i(k)(17)

则由f_i(k)和ε_i(k)的性质，可得到关于V_i(k)的如下关系：

E[V_i(k)]＝E[D_if_i(k)+ε_i(k)]

＝D_iE[f_i(k)]+E[ε_i(k)](18a)

＝0

$\begin{array}{c}E\[V_i\left(k\right){V_i}^T\left(j\right)\]=E\[D_i{f}_i\left(k\right){f_i}^T\left(k\right){D_i}^T+D_i{f}_i\left(k\right){{\mathrm{\ϵ}}_i}^T\left(k\right)+{\mathrm{\ϵ}}_i\left(k\right){f_i}^T\left(k\right){D_i}^T+{\mathrm{\ϵ}}_i\left(k\right){{\mathrm{\ϵ}}_i}^T\left(j\right)\]\\ =\[D_i{Q}_i\left(k\right)D_i^T+R_i\left(k\right)\]{\mathrm{\δ}}_{kj}\\ =r_i\left(k\right){\mathrm{\δ}}_{kj}\end{array}---\left(18b\right)$

由上式可知：V_i(k)为零均值的白噪声随机过程，r_i(k)是V_i(k)的方

差强度矩阵， $r_i\left(k\right)=D_i{Q}_i\left(k\right)D_i^T+R_i\left(k\right).$

而：

$\begin{array}{c}E\[f_i\left(k\right){V_i}^T\left(j\right)\]=E\[f_i\left(k\right){f_i}^T\left(j\right){D_i}^T+f_i\left(k\right){{\mathrm{\ϵ}}_i}^T\left(j\right)\]\\ =Q_i\left(k\right)D_i^T{\mathrm{\δ}}_{kj}\\ =S_i\left(k\right){\mathrm{\δ}}_{kj}\end{array}---\left(19\right)$

式中：S_i(k)是f_i(k)和V_i(k)的协方差强度矩阵，且

综上所诉：f_i(k)和V_i(k)都为零均值的白噪声随机过程，且f_i(k)和V_i(k)相关。因此***状态方程(9)及观测方程(17)满足白噪声相关条件下经典的离散型卡尔曼滤波方程。

4)基于改进卡尔曼滤波理论，利用已测部分楼层的风致响应估计结构未知风致响应分量；

将方程(9)、(17)及相关假设代入卡尔曼滤波基本方程，即可得到kΔt时刻的X_i(k)的估计

${\hat{X}}_i(k/k-1)={\mathrm{\Ψ}}_i{\hat{X}}_i(k-1)+J_i(k-1)\[Z_i(k-1)-H_i{\hat{X}}_i(k-1)\]---\left(20a\right)$

${\hat{X}}_i\left(k\right)={\hat{X}}_i(k/k-1)+G_i\left(k\right)\[Z_i\left(k\right)-H_i{\hat{X}}_i(k/k-1)\]---\left(20b\right)$

J_i(k-1)＝Γ_iQ_i(k-1)D_i ^T[D_iQ_i(k-1)D_i ^T+R_i(k-1)]^-1(20c)

P_i(k/k-1)＝[Ψ_i-J_i(k-1)H_i]P_i(k-1)[Ψ_i-J_i(k-1)H_i]^T+

(20d)

Γ_iQ_i(k-1)Γ_i ^T-J_i(k-1)D_iQ_i(k-1)Γ_i ^T

$G_i\left(k\right)=P_i(k/k-1)H_i^T{\[H_i{P}_i(k/k-1)H_i^T+D_i{Q}_i\left(k\right){D_i}^T+R_i\left(k\right)\]}^{-1}---\left(20e\right)$

P_i(k)＝[I-G_i(k)H_i]P_i(k/k-1)(20f)

式中G_i(k)为最优卡尔曼滤波增益；***状态向量的估计；J_i(k)是状态一步预测增益矩阵，P_i(k/k-1)是一步预测误差方差矩阵，P_i(k)是估计误差方差矩阵。

***初值选取如下：

${\hat{X}}_i\left(0\right)=E\[X_i\left(0\right)\]---\left(21a\right)$

$P_i\left(0\right)=E\{\[X_i\left(0\right)-{\hat{X}}_i\left(0\right)\]{\[X_i\left(0\right)-{\hat{X}}_i\left(0\right)\]}^T\}---\left(21b\right)$

在荷载估计之前外部荷载及测量噪声是未知的，可以先假定荷载协方差矩阵Q_i(k)为单位矩阵。通过改变观测噪声协方差矩阵R_i(k)的值(一般可取10^-4-10^-8)，可以估计最优卡尔曼滤波增益G_i(k)。

根据上诉方法，只要给定初值和P_i(0)，根据kΔt时刻的量测值Z_i(k)，就可递推得到kΔt时刻的***状态向量的估计 ${\hat{X}}_i\left(k\right),(k=1,2,\mathrm{...}).$

5)根据预测的模态响应，估计模态风荷载，进而获得结构任意楼层的风荷载时程；

根据上步求得***状态向量估计将其代入式(9)即可得：

${\mathrm{\Γ}}_i{\hat{f}}_i\left(k\right)={\hat{X}}_i(k+1)-{\mathrm{\Ψ}}_{i^{.}}{\hat{X}}_i\left(k\right)---\left(22\right)$

求得:

${\hat{f}}_i\left(k\right)={{\mathrm{\Γ}}_i}^{+}\[{\hat{X}}_i(k+1)-{\mathrm{\Ψ}}_i{\hat{X}}_i\left(k\right)\]---\left(23\right)$

式中:为Γ_i的广义逆。

基于上述方法依次可求得结构前q阶模态荷载估计将前q阶估计模态荷载组成向量:

${\hat{f}}_{q\×1}={\left[\begin{array}{cccc}{\hat{f}}_1& {\hat{f}}_2& \mathrm{...}& {\hat{f}}_q\end{array}\right]}^T---\left(24\right)$

当超高层建筑的风振分析只考虑前q阶模态时，即可求得结构脉动风荷载向量的估计值

式中为前q列对应的子矩阵。由于模态坐标转换理论可知，振型矩阵关于质量矩阵正交，即：

(Φ_n×n)^TMΦ_n×n＝I(26)

式中I为n×n维单位矩阵，则可由下式求出：

下面结合实施例对本文做进一步详细的说明，但是此说明不会构成对本发明的限制。

实施例1：基于广州市中心某高层建筑的风洞试验，进行风荷载反分析

位于广州市中心的某高层建筑共63层，其中地面以上58层，屋顶标高为256.9m，其所处地形地貌接近于我国《建筑结构荷载规范》GB50009规定的C类粗糙度地区。该建筑平面呈正方形，尺寸为48.00m×48.00m，最大高宽比为5.20，属于风敏感性结构。

该建筑风洞试验在湖南大学边界层风洞进行，风洞的高速试验段为3.0m×2.5m(宽×高)的矩形，该试验段的风速在1.0～58.0m/s内可调。

由于该建筑所在场地为国家规范所划分的C类地貌(风剖面指数为0.22),试验时，利用挡板、尖塔等模拟装置在风洞中形成规定的风剖面(如图1)。风洞试验以主建筑物为中心，模拟半径500m范围内的主要周边建筑及地形，置于风洞试验段转盘上，进行数据测量。试验模型是用ABS板制成的刚体测压模型，具有足够的强度和刚度。模型与实物在外形上保持几何相似，缩尺比为1：300，周边环境模型比例也为1：300。

试验共进行24个风向(0°～360°,间隔15°)的结构表面风压的测量，采样频率为331Hz，采样时间为60s。试验风向角与参考坐标轴定义如图2所示。

第一步：基于改进卡尔曼滤波理论，利用已测部分楼层的风致响应估计未知结构风致响应分量。

本次分析将选取结构加速度响应作为已测的风致响应分量，考虑到加速度响应对应于风荷载中的脉动分量，反演分析时主要获得结构的脉动风荷载分量。基于模态参与系数公式求得该高层结构两个方向前6阶模态响应的能量贡献率均超过了99％，因此反演分析时选择的结构自由度数为6个。选取观测噪声协方差矩阵为：R(k)＝10^-6。输入加速度响应的楼层分别为:第10,20,30,40,50,58层。结构的质量矩阵、刚度矩阵为已知，结构阻尼矩阵选择瑞雷阻尼模型，阻尼比取5％。

基于上述反演方法利用0°及90°风向角下，6层风致响应获得了高层结构未知风致响应分量的分析结果。本实施例给出部分楼层的研究结果。

图3、图4给出了风向角为0°时，结构55层X向(横风向)加速度反演响应时程，作为对比，图中同时给出了结构的准确响应。由图可见,反演得到的位移、速度响应与准确响应时程吻合的非常好。反演的位移、速度响应的均方根为0.029m和0.027m/s，而准确位移响应及速度响应的均方根统计值分别为0.029m和0.026m/s，反演的风致响应分量误差均控制在1％以内。

功率谱是频域内的重要数值特征量，它表征了随机过程的能量分布。图5、图6给出了0°风向角下,结构55层X向(横风向)准确响应和反演响应的功率谱密度。由图可知，位移和速度响应的反演功率谱与准确功率谱在整个频率段均符合的非常好，这意味着本文提出的反演方法能够准确的预测到结构未知响应分量。此外，图中的响应功率谱在结构基频处(0.184Hz)呈现出明显的峰值，这说明结构X向振动以第一阶频率为主。

第二步：根据预测的超高层建筑的模态响应，估计模态风荷载，输出动态风荷载。

利用本发明反演方法并结合预测的风致响应识别了结构各层动态风荷载，图7～图10分别给出了90°风向下35层反演风荷载与准确荷载的时程和功率谱对比结果。由图可知,两个方向的反演风荷载时程均与原始荷载吻合较好，反演风荷载谱与原始功率谱基本一致；由表1可知，反演得到的风荷载峰值及标准差与准确荷载的统计结果差别在8％以内，研究结果表明本发明的反演算法能准确地识别出结构动态风荷载。

表190°风向下结构35层风荷载峰值及标准差(KN)

*差别＝(反演值-准确值)/准确值.

为了进一步评估本发明方法的反演结果的准确性，这里对0°风向角下各楼层反演风荷载沿建筑高度进行积分，获得结构底部总风力的识别结果，如图11～图14所示。作为对比，图中同时给出了结构底部准确总体风荷载的变化曲线。由图可见，在时域及频域内，反演的结构底部总风力均与相应的原始结果吻合良好，这进一步验证了本发明提出的反分析方法的准确性及可靠性。此外，X向(横风向)功率谱有明显的峰值，表现出涡激力谱的特征，而Y向功率谱则为明显的顺风向湍流荷载谱。

结构任意一种风致响应分量(如加速度、速度、位移等)均可作为已知输入来识别结构动态风荷载，但不同的风致响应类型可能会对荷载识别精度有一定影响。考虑到现场实测中常常难以测量结构的速度响应，这里对位移响应及加速度响应的风荷载反演结果开展对比研究。图15、图16对比了0°风向下加速度及位移响应反演的结构底部总风力功率谱。由图可见,基于加速度反演的结构X向(横风向)及Y向(顺风向)风荷载功率谱与原始功率谱吻合的非常好，低频段(小于1Hz)两个方向位移反演风力谱与原始谱也基本吻合，但在高频段，位移响应的反演结果与原始结果有一定差别。总体来讲，加速度反演风荷载精度优于位移的反演精度。

第三步：检验结构模态参数误差、模态截断及多种噪声水平下对该发明方法的反演结果的影响。

准确评估结构模态参数(阻尼、频率等)对预测高层建筑风荷载及风致响应有重要的意义。以往研究结果表明，有限单元法计算的结构动力特性与实测结果往往存在着差异，而实测的结构模态参数本身也常常具有不确定性。这些计算误差及识别不确定性可能会影响风荷载估计精度。本发明将通过对结构模态参数大小进行人为增加(或减小)10％来考察结构模态参数误差对风荷载反演的影响。

图17给出了自振频率误差为±10％时，基于加速度反演的结构X向基底总风力谱对比。由图17可知，在低频段(小于0.2Hz)，自振频率误差对结构风荷载反演结果有一定影响，但基本满足工程需要；在高频段，两个方向的加速度反演荷载功率谱与原始荷载谱吻合较好。图18对比分析了自振频率误差对位移响应反演结果的影响。由图可见，自振频率误差的存在对主要频率段(0.2～1Hz)的结构反演风力谱影响较小，但在低频及高频段，反演风荷载谱与准确谱呈现出一定的差异。分析结果表明，反演结果对自振频率误差不太敏感。

图19和图20分别给出了阻尼比误差为±10％时，加速度及位移反演的结构风荷载谱与准确荷载谱的对比结果。由图可知，加速度反演的风力谱与准确谱符合较好，反演结果基本不受阻尼比误差的影响；位移响应反演的风荷载功率谱密度在高频段与原始风力谱有一定的差异。计算结果表明，与自振频率误差的影响相似，结构阻尼比误差对位移响应反演结果有一定影响,但仍然能够满足工程实际的需要。

表2分别给出了分析模态为前1阶、2阶、4阶和6阶情况下，基于加速度响应反演的结构基底风荷载结果。由表可知，反演分析时当选择的结构模态数超过4阶时，风荷载识别结果的准确度已经能够满足工程需要。

表2不同模态数下加速度反演的风荷载与准确荷载根方差比较(KN)

*差别＝(反演值-准确值)/准确值.

为了检验本发明风荷载反演方法的抗噪声能力，实例分析中将在计算得到的准确风致响应中按照下式迭加一定强度的人工噪声时程，并将含有噪声的动力响应作为输入开展风荷载反向识别。准确响应中拟加入的人工噪声模型为：

d_实测＝d_准确+E_pNoiseσ(d_准确)

式中d_实测为实测的风致响应。d_准确为准确响应。E_p代表噪声强度水平。Noise为利用MTLAB程序中的“randn”函数生成的正态分布随机序列。σ(d_准确)为准确响应的标准差。本次分析共进行了5％、10％两种噪声水平下的对比研究。图21、图22给出了不同噪声水平下，加速度反演风荷载的功率谱对比。由图可见，主要频率段的反演风力谱曲线与原始谱符合较好，本发明方法识别风荷载具有较强的抗噪声能力，在测量噪声影响下，识别风荷载的准确性仍在可接受范围。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (1)

1.一种基于改进卡尔曼滤波的超高层建筑风荷载反分析方法，其特征在于，包括以下步骤：

1)由有限元计算或质量统计可得到结构质量矩阵M，其中超高层建筑的层数为n；

基于现场实测获得结构前q阶结构自振频率ω_i、阻尼比ξ_i，综合应用现场实测和有限元分析获得结构模态振型Φ_n×q；

2)输入实测的p个楼层的风致响应分量，根据实测的结构前q阶模态振型，将测试的风致响应分量转化为模态风致响应；所述风致响应分量为位移、速度或加速度响应中的一种；

在结构动态响应实测时，若只测得结构部分楼层的风致响应(假设为p层加速度响应)和前q阶模态振型，由于超高层建筑的风致振动往往以前几阶模态为主，则结构风致响应可近似分解为：

${\stackrel{\·\·}{y}}_{p\×1}={\mathrm{\Φ}}_{p\×q}{\stackrel{\·\·}{U}}_{q\×1},(1\≤q\≤p\≤n)$

其中，为实测的p层风致响应分量；为前q阶模态向量；Φ_p×q由Φ_n×n中与p个实测响应楼层对应的行、前q列形成的子振型矩阵；Φ_n×n为按质量规准化的模态振型矩阵；

确定结构振动的主要控制模态数目q采用以下方法：基于POD方法首先获取加速度响应协方差矩阵的特征值λ_i(i＝1,2,…n)，进而计算出前q阶模态对结构振动的贡献比例：

$\mathrm{\θ}=\frac{\underset{i=1}{\overset{q}{\Σ}}{\mathrm{\λ}}_i}{\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\λ}}_i},(1\≤q\≤n)$

取θ超过预设值时所对应的q值作为结构振动主要控制模态数目。

由广义逆矩阵Φ_p×q ⁺，结构实测的模态加速度响应可近似表示为：

${\widehat{\stackrel{\·\·}{U}}}_{q\×1}={\left({\mathrm{\Φ}}_{p\×q}\right)}^{+}{\stackrel{\·\·}{y}}_{p\×1}$

式中：(Φ_p×q)⁺为Φ_p×q的广义逆矩阵；

3)根据结构动力微分方程，在模态空间中构造离散化的状态方程和观测方程；

风荷载作用下有n个楼层的超高层建筑的动力微分方程可表示如下：

$M\stackrel{\·\·}{y}+C\stackrel{\·}{y}+Ky=F$

式中，y、和是位移、速度和加速度向量；F为外荷载；C和K为阻尼矩阵和刚度矩阵；

动力微分方程可按如下解耦：

${\mathrm{\Φ}}_i^{T}F=f_i=M_i{\stackrel{\·\·}{U}}_i+C_i{\stackrel{\·}{U}}_i+K_i{U}_i$

式中：M_i、K_i分别为第i阶按质量规准化的模态质量和模态刚度，其中M_i＝1(i＝1,2,…,q),f_i、C_i分别为第i阶按质量规准化的模态荷载，其中C_i＝2ξ_iω_i；U_i分别为第i阶模态加速度、速度、位移

由结构动力微分方程构造的离散型***状态方程：

X_i(k+1)＝Ψ_iX_i(k)+Γ_if_i(k)；

其中，X_i(k)为kΔt时刻的状态向量，f_i(k)为kΔt时刻的模态荷载；Δt为采样间隔；Ψ_i为kΔt时刻至(k+1)Δt时刻的一步转移矩阵；Γ_i为***噪声驱动矩阵；

***观测方程为如下形式：

Z_i(k)＝H_iX_i(k)+V_i(k)

其中，Z_i(k)为kΔt时刻的响应观测值；H_i为观测矩阵；V_i(k)＝D_if_i(k)+ε_i(k)，其中，

D_i为***矩阵；ε_i(k)为观测噪声；

4)基于改进卡尔曼滤波理论，利用已测部分楼层的风致响应估计结构未知风致响应分量；

将状态方程和观测方程及相关假设代入卡尔曼滤波基本方程，即可得到kΔt时刻的X_i(k)的估计

${\hat{X}}_i(k/k-1)={\mathrm{\Ψ}}_i{\hat{X}}_i(k-1)+J_i(k-1)\[Z_i(k-1)-H_i{\hat{X}}_i(k-1)\]$

${\hat{X}}_i\left(k\right)={\hat{X}}_i(k/k-1)+G_i\left(k\right)\[Z_i\left(k\right)-H_i{\hat{X}}_i(k/k-1)\]$

J_i(k-1)＝Γ_iQ_i(k-1)D_i ^T[D_iQ_i(k-1)D_i ^T+R_i(k-1)]^-1

P_i(k/k-1)＝[Ψ_i-J_i(k-1)H_i]P_i(k-1)[Ψ_i-J_i(k-1)H_i]^T+

Γ_iQ_i(k-1)Γ_i ^T-J_i(k-1)D_iQ_i(k-1)Γ_i ^T

$G_i\left(k\right)=P_i(k/k-1)H_i^T{\[H_i{P}_i(k/k-1)H_i^T+D_i{Q}_i\left(k\right){D_i}^T+R_i\left(k\right)\]}^{-1}$

P_i(k)＝[I-G_i(k)H_i]P_i(k/k-1)

***初值选取如下：

${\hat{X}}_i\left(0\right)=E\[X_i\left(0\right)\];$

$P_i\left(0\right)=E\{\[X_i\left(0\right)-{\hat{X}}_i\left(0\right)\]{\[X_i\left(0\right)-{\hat{X}}_i\left(0\right)\]}^T\};$

其中，Q_i(k)为荷载协方差矩阵，R_i(k)为观测噪声协方差矩阵，G_i(k)为最优卡尔曼滤波增益；***状态向量的估计；J_i(k)是状态一步预测增益矩阵，P_i(k/k-1)是一步预测误差方差矩阵，P_i(k)是估计误差方差矩阵。

5)根据预测的模态响应，估计模态风荷载，进而获得结构任意楼层的风荷载时程；

具体如下：根据步骤4)中求得的***状态向量估计将其离散型***状态方程即可得：

${\mathrm{\Γ}}_i{\hat{f}}_i\left(k\right)={\hat{X}}_i(k+1)-{\mathrm{\Ψ}}_i{\hat{X}}_i\left(k\right)$

求得:

${\hat{f}}_i\left(k\right)={{\mathrm{\Γ}}_i}^{+}\[{\hat{X}}_i(k+1)-{\mathrm{\Ψ}}_i{\hat{X}}_i\left(k\right)\]$

式中:Γ_i ⁺为Γ_i的广义逆；

基于上述方法依次可求得结构前q阶模态荷载估计将前q阶估计模态荷载组成向量:

${\hat{f}}_{q\×1}={\left[\begin{array}{cccc}{\hat{f}}_1& {\hat{f}}_2& \mathrm{...}& {\hat{f}}_q\end{array}\right]}^T$

当超高层建筑的风振分析只考虑前q阶模态时，即可求得结构脉动风荷载向量的估计值

式中为前q列对应的子矩阵；由于模态坐标转换理论可知，振型矩阵关于质量矩阵正交，即：

(Φ_n×n)^TMΦ_n×n＝I

式中I为n×n维单位矩阵，则可由下式求出：