CN105260568B - 基于离散型卡尔曼滤波的超高层建筑风荷载反分析方法 - Google Patents

Abstract

本发明公开了一种基于离散型卡尔曼滤波的超高层建筑风荷载反分析方法，该方法包括以下步骤：（1）基于实测的前几阶模态振型将结构有限测试楼层的风致响应分量（位移或速度）转化为模态风致响应；（2）利用离散型卡尔曼滤波估计未知的模态风致响应分量；（3）通过估计的模态风致响应（位移、速度、加速度）识别结构的模态风荷载；（4）利用模态振型矩阵的广义逆获得结构任意楼层的风荷载时程。该方法可以解决超高层建筑风致响应测点不足的问题，本方法在结构模态参数误差、模态截断误差及测量噪声影响下的风荷载反演结果仍然能够满足工程需要。研究技术为超高层建筑抗风设计及相关研究提供有用的工具及依据。

Description

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

技术领域

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

背景技术

近年来在我国东南沿海地区兴建了大量三百米以上的超高层建筑,由于这类结构的自振频率较低、阻尼较小，同台风动荷载的主要频率段比较接近,在强/台风作用下的风致响应比较大，风荷载及风振响应是其安全性及适用性设计的首要控制指标。一般的来说，超高层建筑的风荷载很难实现现场测量，风洞试验技术尽管能够测试结构风荷载，但由于其手段较为复杂、试验技术及分析方法尚不完善，风荷载评估结果难以准确的反映真实状态，目前对于风与结构的相互作用机理(特别是横风向及扭转)还存在着很多疑问。近年来，国内外兴起了一种基于结构响应反演动态荷载的研究方法，考虑到目前结构动力响应的实时测量技术比较成熟，量测精度远高于荷载的量测精度。因此，根据结构动力特性的测试与分析结果，以实测结构风致响应为基础，以反分析方法为手段来获取结构动态风荷载，成为近代兴起的间接量测风荷载的一种新途径。

近年来，国内外学者在结构动态荷载的反分析方面已做出了有益的探索与尝试。如Liu等(2000)利用反分析算子尝试识别了悬挑板上的稳态动荷载。Ma以及Ho等人(2003,2004)基于数值分析评估了悬臂梁上的动荷载，并将其扩展到非线性结构系统。Liu andShepard(2005)发展了一种动态荷载的频域反演方法。Lu和Law(2006)提出了基于动力响应识别结构动态激励的方法。陈建云等人(2006)基于最小二乘法开展了未知输入下的复合反演研究。在建筑结构风荷载反分析方面，Kang和Lo(2002)对典型高塔进行了风荷载反演分析。Law等人(2005)提出了一种评估风荷载的反分析算法，并以桅杆为对象进行了数值验证。Hwang及Ahsan Kareem等(2009,2010)基于混凝土烟囱风洞试验中测得的结构动力响应评估了横风向荷载，并对比研究了不同响应类型对烟囱风荷载分析结果的影响。

以上研究现状表明，基于实测结构动力响应的动荷载反分析方法是一种获取结构动态风荷载的有效手段，它对全面了解结构风荷载作用机理及风振响应规律具有重要的意义。但目前来说对建筑结构风荷载的反分析研究较少，对于超高层建筑的风荷载反演分析研究更是有限，因此,本发明提出一种基于离散型卡尔曼滤波的超高层建筑风荷载反分析方法。

发明内容

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

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

1)利用有限单元法获取超高层建筑的质量矩阵M、刚度矩阵K和阻尼矩阵C，超高层建筑的层数为n；

2)输入实测的P个楼层的风致响应分量，所述风致响应为位移或速度响应分量中的一种；

3)根据结构前q阶模态振型，将实测的风致响应分量转化为模态风致响应，并在模态空间中构造离散化的状态方程和测量方程；

3.1)将实测的风致响应分量转化到模态空间；

y_p×1＝Φ_p×q·U_q×1(1≤p≤n,1≤q≤n)

式中y_p×1为p个楼层风致响应，Φ_p×q为对应最高阶为q阶的模态振型矩阵，U_q×1为前q阶模态响应；

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

式中为U_q×1的估计值；

准确模态位移与估计模态位移之间误差向量可以用下式表示:

上式中估计误差的减小可通过增加实测楼层的数目，使其大于结构振动主要控制模态数来实现。为了确定结构振动的主要控制模态数目,本发明将基于POD方法首先获取位移响应协方差矩阵的特征值λ_i(i＝1,2,…n)，进而计算出前q阶模态对结构振动的贡献比例：

本发明取θ超过99％时所对应的q值作为结构振动主要控制模态数目。当位移传感器的数目p≥q时，模态位移之间的误差将会满足计算精度的要求。

3.2)将质量矩阵M、刚度矩阵K和阻尼矩阵C按质量归一化进行模态转化；

式中Φ_i、M_i、K_i分别为第i阶按质量规准化的振型向量、模态质量和模态刚度；

3.3)构造离散化的状态方程和测量方程；

动力方程可如下解耦为：

式中f_i、C_i分别为第i阶按质量规准化的模态荷载、模态阻尼；U_i分别为第i阶模态加速度、模态速度、模态位移，应用Taylor展开式有：

式中Δt是采样间隔，j_i(k-1)和分别是k-1时间点处的模态加速度的一阶导数和二阶导数，

令

则X_i(k)用离散型的状态方程可以如下表示：

其中A表示状态传递矩阵，B表示噪声矩阵，如下

离散型的测量方程可以表示如下：

Z_i(k)＝H·X_i(k)+W_i(k)

其中Z_i(t)表示测量的模态响应，W_i(k)表示测量噪声矢量，H表示测量矩阵，取决于测量响应的类型，如果测量响应是位移，H＝[1 0 0 0]，如果测量响应是速度，H＝[0 1 00]。

假设过程噪声W_i(k)是零均值的白噪声，其协方差可以分别用Q_i、R_i表示如下：

E[W_i(k)]＝0,Cov[W_i(k),W_i(m)]＝E[W_i(k)W_i(m)^T]＝R_iδ_km

这里δ_km是柯氏；

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

其中，

P_i(k/k-1)＝A P_i(k-1)A^T+BQ_iB^T；

G_i(k)＝P_i(k/k-1)H^T(HP_i(k/k-1)H^T+R_i)^-1；

P_i(k)＝(I-G_i(k)H)P_i(k/k-1)；

式中，G_i(k)是时间点k处的卡尔曼滤波增益，P_i(k)表示滤波误差协方差矩阵；I为单位矩阵，Q_i为噪声协方差，A表示状态传递矩阵，B表示噪声矩阵，

5)根据预测的模态响应，估计模态风荷载，进而利用模态振型矩阵的广义逆获得结构任意楼层的风荷载时程；

其中则可估计结构风致外荷载

式中Φ_n×q为振型矩阵Φ_n×n的前q列，为估计的模态荷载向量

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

附图说明

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

图1是本发明实施例中平均风速剖面及湍流剖面；

图2是本发明实施例中参考坐标轴；

图3是本发明实施例中90°风向下第80层X向速度响应时程；

图4是本发明实施例中90°风向下第80层X向加速度响应时程；

图5是本发明实施例中90°风向下结构80层X向速度响应功率谱；

图6是本发明实施例中90°风向下结构80层X向加速度响应功率谱；

图7是本发明实施例中90°风向下结构底部X向总风力时程；

图8是本发明实施例中90°风向下结构底部Y向总风力时程；

图9是本发明实施例中90°风向角下结构底部X向总风力功率谱

图10是本发明实施例中90°风向角下结构底部Y向总风力功率谱；

图11是本发明实施例中自振频率误差±10％时，位移反演X向基底总风力功率谱(90°风向)；

图12是本发明实施例中自振频率误差±10％时，位移反演Y向基底总风力功率谱(90°风向)；

图13是本发明实施例中阻尼比误差±10％时，位移反演X向基底总风力功率谱(90°风向)；

图14是本发明实施例中阻尼比误差±10％时，位移反演Y向基底总风力功率谱(90°风向)；

图15是本发明实施例中90°风向下2％噪声水平时X向准确荷载和反演荷载的相关图；

图16是本发明实施例中90°风向下2％噪声水平时Y向准确荷载和反演荷载的相关图；

图17是本发明实施例中90°风向下5％噪声水平时X向准确荷载和反演荷载的相关图；

图18是本发明实施例中90°风向下5％噪声水平时Y向准确荷载和反演荷载的相关图；

图19是本发明实施例的方法流程图。

具体实施方式

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

如图19所示，基于离散型卡尔曼滤波的超高层建筑风荷载反分析方法,包括以下步骤：

1)利用有限单元法获取超高层建筑的质量矩阵M、刚度矩阵K和阻尼矩阵C，超高层建筑的层数为n；

2)输入实测的P个楼层的风致响应分量，所述风致响应为位移或速度响应分量中的一种；

3)根据结构前q阶模态振型，将实测的风致响应分量转化为模态风致响应，并在模态空间中构造离散化的状态方程和测量方程。

3.1)将实测的风致响应分量转化到模态空间；

y_p×1＝Φ_p×q·U_q×1(1≤p≤n,1≤q≤n)

式中y_p×1为p个楼层风致响应，Φ_p×q为对应最高阶为q阶的模态振型矩阵，U_q×1为前q阶模态响应。

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

式中为U_q×1的估计值。准确模态位移与估计模态位移之间误差向量可以用下式表示:

上式中估计误差的减小可通过增加实测楼层的数目，使其大于结构振动主要控制模态数来实现。为了确定结构振动的主要控制模态数目,本发明将基于POD方法首先获取位移响应协方差矩阵的特征值λ_i(i＝1,2,…n)，进而计算出前q阶模态对结构振动的贡献比例：

本发明取θ超过99％时所对应的q值作为结构振动主要控制模态数目。当位移传感器的数目p≥q时，模态位移之间的误差将会满足计算精度的要求。

3.2)将质量矩阵M、刚度矩阵K和阻尼矩阵C按质量归一化进行模态转化；

式中Φ_i、M_i、K_i分别为第i阶按质量规准化的振型向量、模态质量和模态刚度；

3.3)构造离散化的状态方程和测量方程；

动力方程可如下解耦为：

式中f_i、C_i分别为第i阶按质量规准化的模态荷载、模态阻尼；U_i分别为第i阶模态加速度、模态速度、模态位移，应用Taylor展开式有：

式中Δt是采样间隔，j_i(k-1)和分别是k-1时间点处的模态加速度的一阶导数和二阶导数，

令

则X_i(k)用离散型的状态方程可以如下表示：

其中A表示状态传递矩阵，B表示噪声矩阵，如下

离散型的测量方程可以表示如下：

Z_i(k)＝H·X_i(k)+W_i(k)

其中Z_i(t)表示测量的模态响应，W_i(k)表示测量噪声矢量，H表示测量矩阵，取决于测量响应的类型，如果测量响应是位移，H＝[1 0 0 0]，如果测量响应是速度，H＝[0 1 00]。

假设过程噪声W_i(k)是零均值的白噪声，其协方差可以分别用Q_i、R_i表示如下：

E[W_i(k)]＝0,Cov[W_i(k),W_i(m)]＝E[W_i(k)W_i(m)^T]＝R_iδ_km

这里δ_km是柯氏δ。

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

其中，

P_i(k/k-1)＝A P_i(k-1)A^T+BQ_iB^T；

G_i(k)＝P_i(k/k-1)H^T(HP_i(k/k-1)H^T+R_i)^-1；

P_i(k)＝(I-G_i(k)H)P_i(k/k-1)；

式中，G_i(k)是时间点k处的卡尔曼滤波增益，P_i(k)表示滤波误差协方差矩阵；I为单位矩阵，Q_i为噪声协方差，A表示状态传递矩阵，B表示噪声矩阵，

5)根据预测的模态响应，估计模态风荷载，进而利用模态振型矩阵的广义逆获得结构任意楼层的风荷载时程：

其中

则可估计结构风致外荷载

式中，Φ_n×q为振型矩阵Φ_n×n的前q列，为估计的模态荷载向量

下面结合实例对该发明做进一步详细的说明，该实例用于验证本发明的有效性。

验证实例：基于香港某超高层建筑的风洞试验数据进行风荷载反演分析

该超高层建筑及基本信息如下：高420米,地面以上共88层，该塔楼结构的平面布置为方形，底部尺寸为57m×57m，到顶部逐渐变化为39m×39m，高宽比大约为8，属于典型的风敏感性结构。风洞试验在香港城市大学的边界层风洞试验室进行，模拟边界层流场，风剖面指数α＝0.28，试验时，利用挡板、尖塔等模拟装置在风洞中形成规定的风剖面(如图1所示)。风洞试验段的是4米宽×2米高,模型与实物在外形上保持几何相似，缩尺比为1：400，周边环境模型比例也为1：400。试验共进行24个风向(0°～360°,间隔15°)的结构表面风压的测量，采样频率为600Hz，试验风向角与参考坐标轴定义如图2所示。

实例分析时利用风洞试验的测试数据,结合我国现行荷载规范的有关规定，确定了对应于50年重现期(基本风速为59.5m/s)结构各层风荷载时程(在对比分析时作为准确荷载),再运用结构动力分析得到风荷载作用下结构各层加速度、速度和位移响应，这些结果将用于风荷载反演分析时的对比研究。

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

本次分析将分别选取位移和速度响应作为已测的风致响应分量反演得到结构脉动风荷载。基于模态参与系数公式求得该高层结构两个方向前4阶模态响应的能量贡献率均超过了99％，因此反演分析时选择的结构自由度数为4个。输入相应响应的楼层分别为:第24,39,64和88层。结构的质量矩阵、刚度矩阵为已知，结构阻尼矩阵选择瑞雷阻尼模型，阻尼比取5％。

基于所选择的四个楼层的位移响应，利用本发明方法可以估计任意楼层的风致响应(假设P_i(0/-1)＝10⁶I)。图3、图4给出了风向角为90°时，80层X向的速度和加速度响应，作为比较，给出对应的基于风洞试验计算的准确响应。由图可见,估计得到的速度、加速度响应与准确响应时程吻合的非常好。图5、图6给出了90°风向角下,结构80层X向准确响应和反演响应的功率谱密度。由图可知，位移和速度响应的反演功率谱与准确功率谱在整个频率段均符合的非常好，这意味着本发明提出的反演方法能够准确的预测到结构未知响应分量。此外，两个明显的峰值出现在图5、图6中0.143Hz和0.342Hz处，恰好分别对应于结构X向的两个自振频率。

b.第二步：根据预测的模态响应，估计模态风荷载，进而得到结构各个楼层的脉动风荷载。

利用本发明反演方法并结合预测的风致响应识别了结构各层动态风荷载，并对各楼层反演风荷载沿建筑高度进行积分，获得了结构底部总风力的识别结果。图7～图10给出了90°风向角下分别基于位移及速度响应识别的风荷载结果。作为对比，图中同时给出了结构底部准确总体风荷载的变化曲线。由图可见，在时域及频域内，反演的结构底部总风力均与相应的原始结果吻合良好，这验证了本发明反分析方法的准确性及可靠性。

此外，表1、表2分别给出0°、90°、180°和270°风向角时，不同响应类型反演的结构X向和Y向风荷载与准确值的标准差及其差别。由表可以看出，两种响应反演的风荷载非常接近准确值，进一步说明了本发明方法的准确性。c.第三步：评估结构模态参数误差、模态截断及测量噪声对反演结果的影响。

实例分析将通过对结构模态参数大小进行人为增加(或减小)10％来考察结构模态参数误差对风荷载反演结果的影响。图11和图12给出了自振频率误差为±10％时，基于位移响应反演的结构X向和Y向基底总风力谱对比。由图可知，自振频率误差对结构风荷载反演结果有一定影响，但基本满足工程需要。图13和图14分别给出了阻尼比误差为±10％时，基于位移响应反演的X向和Y向结构基底风荷载谱与准确荷载谱的对比结果。由图可知，反分析结果对阻尼比误差不敏感。

表1 位移反演的风荷载与准确荷载值比较

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

表2 速度反演的风荷载与准确荷载值比较

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

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

表3 不同模态数情况下位移响应反演的风荷载与准确荷载根方差比较(风向为90°)

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

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

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

式中d_实测为实测的风致响应。d_准确为准确响应。E_p代表噪声强度水平。Noise为利用MTLAB程序中的“randn”函数生成的正态分布随机序列。σ(d_准确)为准确响应的标准差。

图15、图16分别给出了在2％噪声水平下90°风向角时，X向和Y向的准确荷载和反演荷载的相关图，相关系数分别是0.90和0.89。图17、图18分别给出了在5％噪声水平下90°风向角时，X向和Y向的准确荷载和反演荷载的相关图，相关系数分别是0.84和0.82。由图可知，本方法识别风荷载有一定的抗噪声能力，在测量噪声影响下，识别风荷载的准确性仍在可接受范围。

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

Claims (3)

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

1)利用有限单元法获取超高层建筑的质量矩阵M、刚度矩阵K和阻尼矩阵C，超高层建筑的层数为n；

2)输入实测的p个楼层的风致响应分量，所述风致响应为位移或速度响应分量中的一种；

3)根据结构前q阶模态振型，将实测的风致响应分量转化为模态风致响应，并在模态空间中构造离散化的状态方程和测量方程；

3.1)将实测的风致响应分量转化到模态空间；

y_p×1＝Φ_p×q·U_q×1,1≤p≤n,1≤q≤n；

式中y_p×1为p个楼层风致响应，Φ_p×q为对应最高阶为q阶的模态振型矩阵，U_q×1为前q阶模态响应；

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

<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mover> <mi>U</mi> <mo>^</mo> </mover> <mrow> <mi>q</mi> <mo>&amp;times;</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>&amp;Phi;</mi> <mrow> <mi>p</mi> <mo>&amp;times;</mo> <mi>q</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>+</mo> </msup> <mo>&amp;CenterDot;</mo> <msub> <mi>y</mi> <mrow> <mi>p</mi> <mo>&amp;times;</mo> <mn>1</mn> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>U</mi> <mo>^</mo> </mover> <mrow> <mi>q</mi> <mo>&amp;times;</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>&amp;Phi;</mi> <mrow> <mi>p</mi> <mo>&amp;times;</mo> <mi>q</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>-</mo> </msup> <mo>&amp;CenterDot;</mo> <msub> <mi>y</mi> <mrow> <mi>p</mi> <mo>&amp;times;</mo> <mn>1</mn> </mrow> </msub> <mo>,</mo> <mi>i</mi> <mi>f</mi> <mi> </mi> <mi>p</mi> <mo>=</mo> <mi>q</mi> </mrow> </mtd> </mtr> </mtable> </mfenced>

式中为U_q×1的估计值；

准确模态位移与估计模态位移之间误差向量用下式表示:

<mrow> <mi>&amp;Delta;</mi> <mo>=</mo> <msub> <mi>U</mi> <mrow> <mi>q</mi> <mo>&amp;times;</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mover> <mi>U</mi> <mo>^</mo> </mover> <mrow> <mi>q</mi> <mo>&amp;times;</mo> <mn>1</mn> </mrow> </msub> </mrow>

3.2)将质量矩阵M、刚度矩阵K和阻尼矩阵C按质量归一化进行模态转化；

<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>M</mi> <mi>i</mi> </msub> <mo>=</mo> <msubsup> <mi>&amp;Phi;</mi> <mi>i</mi> <mi>T</mi> </msubsup> <msub> <mi>M&amp;Phi;</mi> <mi>i</mi> </msub> <mo>=</mo> <mn>1</mn> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>n</mi> <mo>;</mo> </mrow> </mtd> <mtd> <mrow> <msub> <mi>K</mi> <mi>i</mi> </msub> <mo>=</mo> <msubsup> <mi>&amp;Phi;</mi> <mi>i</mi> <mi>T</mi> </msubsup> <msub> <mi>K&amp;Phi;</mi> <mi>i</mi> </msub> <mo>=</mo> <msubsup> <mi>&amp;omega;</mi> <mi>i</mi> <mn>2</mn> </msubsup> </mrow> </mtd> <mtd> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>n</mi> </mrow> </mtd> </mtr> </mtable> </mfenced>

式中Φ_i、M_i、K_i分别为第i阶按质量规准化的振型向量、模态质量和模态刚度；

3.3)构造离散化的状态方程和测量方程；

动力方程如下解耦为：

<mrow> <msub> <mi>f</mi> <mi>i</mi> </msub> <mo>=</mo> <msub> <mover> <mi>U</mi> <mo>&amp;CenterDot;&amp;CenterDot;</mo> </mover> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>C</mi> <mi>i</mi> </msub> <msub> <mover> <mi>U</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>i</mi> </msub> <mo>+</mo> <msubsup> <mi>&amp;omega;</mi> <mi>i</mi> <mn>2</mn> </msubsup> <msub> <mi>U</mi> <mi>i</mi> </msub> </mrow>

式中f_i、C_i分别为第i阶按质量规准化的模态荷载、模态阻尼；U_i分别为第i阶模态加速度、模态速度、模态位移，应用Taylor展开式有：

<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>U</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>U</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mover> <mi>U</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>+</mo> <msub> <mover> <mi>U</mi> <mo>&amp;CenterDot;&amp;CenterDot;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <msup> <mi>&amp;Delta;t</mi> <mn>2</mn> </msup> <mo>/</mo> <mn>2</mn> <mo>+</mo> <msub> <mi>j</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <msup> <mi>&amp;Delta;t</mi> <mn>3</mn> </msup> <mo>/</mo> <mn>6</mn> <mo>+</mo> <msub> <mover> <mi>j</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <msup> <mi>&amp;Delta;t</mi> <mn>4</mn> </msup> <mo>/</mo> <mn>24</mn> <mo>+</mo> <mi>o</mi> <mrow> <mo>&amp;lsqb;</mo> <mrow> <msup> <mi>&amp;Delta;t</mi> <mn>4</mn> </msup> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>U</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mover> <mi>U</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mover> <mi>U</mi> <mo>&amp;CenterDot;&amp;CenterDot;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>+</mo> <msub> <mi>j</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <msup> <mi>&amp;Delta;t</mi> <mn>2</mn> </msup> <mo>/</mo> <mn>2</mn> <mo>+</mo> <msub> <mover> <mi>j</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <msup> <mi>&amp;Delta;t</mi> <mn>3</mn> </msup> <mo>/</mo> <mn>6</mn> <mo>+</mo> <mi>o</mi> <mrow> <mo>&amp;lsqb;</mo> <mrow> <msup> <mi>&amp;Delta;t</mi> <mn>3</mn> </msup> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>U</mi> <mo>&amp;CenterDot;&amp;CenterDot;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mover> <mi>U</mi> <mo>&amp;CenterDot;&amp;CenterDot;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>j</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>+</mo> <msub> <mover> <mi>j</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <msup> <mi>&amp;Delta;t</mi> <mn>2</mn> </msup> <mo>/</mo> <mn>2</mn> <mo>+</mo> <mi>o</mi> <mrow> <mo>&amp;lsqb;</mo> <mrow> <msup> <mi>&amp;Delta;t</mi> <mn>2</mn> </msup> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>j</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>j</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mover> <mi>j</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>+</mo> <mi>o</mi> <mrow> <mo>&amp;lsqb;</mo> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced>

式中Δt是采样间隔，j_i(k-1)和分别是k-1时间点处的模态加速度的一阶导数和二阶导数，

令

则X_i(k)用离散型的状态方程如下表示：

<mrow> <msub> <mi>X</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>A</mi> <mo>&amp;CenterDot;</mo> <msub> <mi>X</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>+</mo> <mi>B</mi> <mo>&amp;CenterDot;</mo> <msub> <mover> <mi>j</mi> <mo>&amp;CenterDot;</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>

其中A表示状态传递矩阵，B表示噪声矩阵，如下

<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>A</mi> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </mtd> <mtd> <mrow> <msup> <mi>&amp;Delta;t</mi> <mn>2</mn> </msup> <mo>/</mo> <mn>2</mn> </mrow> </mtd> <mtd> <mrow> <msup> <mi>&amp;Delta;t</mi> <mn>3</mn> </msup> <mo>/</mo> <mn>6</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </mtd> <mtd> <mrow> <msup> <mi>&amp;Delta;t</mi> <mn>2</mn> </msup> <mo>/</mo> <mn>2</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> </mtr> </mtable> </mfenced> </mrow> </mtd> <mtd> <mrow> <mi>B</mi> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mi>&amp;Delta;</mi> <msup> <mi>t</mi> <mn>4</mn> </msup> <mo>/</mo> <mn>24</mn> </mtd> </mtr> <mtr> <mtd> <mi>&amp;Delta;</mi> <msup> <mi>t</mi> <mn>3</mn> </msup> <mo>/</mo> <mn>6</mn> </mtd> </mtr> <mtr> <mtd> <mi>&amp;Delta;</mi> <msup> <mi>t</mi> <mn>2</mn> </msup> <mo>/</mo> <mn>2</mn> </mtd> </mtr> <mtr> <mtd> <mi>&amp;Delta;</mi> <mi>t</mi> </mtd> </mtr> </mtable> </mfenced> </mrow> </mtd> </mtr> </mtable> </mfenced>

离散型的测量方程表示如下：

Z_i(k)＝H·X_i(k)+W_i(k)

其中Z_i(k)表示测量的模态响应，W_i(k)表示测量噪声矢量，H表示测量矩阵，取决于测量响应的类型，如果测量响应是位移，H＝[1 0 0 0]，如果测量响应是速度，H＝[0 1 0 0]；

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

<mrow> <msub> <mover> <mi>X</mi> <mo>^</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mover> <mi>X</mi> <mo>^</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>/</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>G</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mrow> <msub> <mi>Z</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>H</mi> <msub> <mover> <mi>X</mi> <mo>^</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>/</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

其中，

<mrow> <msub> <mover> <mi>X</mi> <mo>^</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>/</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <mi>A</mi> <msub> <mover> <mi>X</mi> <mo>^</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

P_i(k/k-1)＝A P_i(k-1)A^T+BQ_iB^T；

G_i(k)＝P_i(k/k-1)H^T(HP_i(k/k-1)H^T+R_i)^-1；

P_i(k)＝(I-G_i(k)H)P_i(k/k-1)；

式中，G_i(k)是时间点k处的卡尔曼滤波增益，P_i(k)表示滤波误差协方差矩阵；I为单位矩阵，Q_i为噪声协方差，A表示状态传递矩阵，B表示噪声矩阵；

5)根据预测的模态响应，估计模态风荷载，进而利用模态振型矩阵的广义逆获得结构任意楼层的风荷载时程；

<mrow> <msub> <mover> <mi>f</mi> <mo>^</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&amp;Psi;</mi> <mi>i</mi> </msub> <mo>&amp;CenterDot;</mo> <msub> <mover> <mi>X</mi> <mo>^</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow>

其中则估计结构风致外荷载

<mrow> <mover> <mi>F</mi> <mo>^</mo> </mover> <mo>=</mo> <msup> <mrow> <mo>(</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>&amp;Phi;</mi> <mrow> <mi>n</mi> <mo>&amp;times;</mo> <mi>q</mi> </mrow> </msub> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>)</mo> </mrow> <mo>+</mo> </msup> <mover> <mi>f</mi> <mo>^</mo> </mover> </mrow>

式中Φ_n×q为振型矩阵Φ_n×n的前q列，为估计的模态荷载向量

2.根据权利要求1所述的超高层建筑风荷载反分析方法，其特征在于，所述步骤3.1)中准确模态位移与估计模态位移之间误差向量通过增加实测楼层的数目减小。

3.根据权利要求1所述的超高层建筑风荷载反分析方法，其特征在于，所述步骤3.1)中确定结构振动的主要控制模态数目采用以下方法：基于POD方法首先获取位移响应协方差矩阵的特征值λ_i,i＝1,2,…n，进而计算出前q阶模态对结构振动的贡献比例：

<mrow> <mi>&amp;theta;</mi> <mo>=</mo> <mfrac> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>q</mi> </munderover> <msub> <mi>&amp;lambda;</mi> <mi>i</mi> </msub> </mrow> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msub> <mi>&amp;lambda;</mi> <mi>i</mi> </msub> </mrow> </mfrac> <mo>,</mo> <mn>1</mn> <mo>&amp;le;</mo> <mi>q</mi> <mo>&amp;le;</mo> <mi>n</mi> </mrow>

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