付康胜，吕中荣，汪利

中山大学航空航天学院，广东广州510006

在我国经济发展较好的沿海城市，为了提高土地利用效率，超高层建筑往往被当作首选目标。而，在沿海地区台风频发，风荷载对建筑的影响甚是严重。超高层建筑柔性较高，低阶固有频率往往很低，同时风荷载也属于低频荷载，建筑在风荷载作用下很容易发生共振。玻璃幕墙在强风作用下容易发生损毁，高层建筑在受到绕流风作用时，建筑边缘处易发生旋涡的脱落和气流的再附，从而引起横风向的振动，这对建筑的舒适性影响尤其严重［1］。因此，需要对作用于建筑上的风荷载进行测量，在设计时使其固有频率尽可能地避开风荷载频率。考虑到建筑表面的复杂性，基于现有技术很难把荷载直接精确地测量出来。近年来有学者提出基于动力响应反演风荷载，动力响应数据的测量精度也远高于荷载的测量精度。所以，利用实测动力响应数据间接识别风荷载能获得更加接近真实风荷载时程的数据。

国内外的学者对结构荷载反演进行了一些研究。谭栋等［2］利用新兴群智能算法对作用于梁上的冲击荷载进行识别，发现该方法抗噪性能较好。陈隽等［3］在结构物理参数未知前提下，根据风荷载作用特点，利用荷载归一化平均法，获得作用在结构上的风荷载时程。方明新［4］基于Kalman 滤波理论和超高层建筑的风振特性，研究了利用有限的风响应数据反演超高层脉动风荷载的时程。Hwang 等［5］对高耸混凝土烟囱进行了风洞试验，利用卡尔曼滤波方法和有限的动力响应测量数据反演风荷载。Azam 等［6］提出一种对偶卡尔曼滤波方法识别作用在米兰摩天大楼的冲击荷载。Lourens 等［7］提出了一种增广卡尔曼滤波方法识别作用在梁上的冲击荷载，得到较好的识别结果。郅伦海等［8］基于遗忘因子最小二乘法和离散卡尔曼滤波方法识别高层建筑的脉动风荷载，并用风洞实验验证该方法的可靠性。叶静娴等［9］基于Newmark 方法和最小二乘法构建目标函数，反演出作用在海洋平台的波浪荷载。李正农等［10］提出一种基于现场实测高层建筑加速度响应反演脉动风荷载的方法，并将脉动风荷载与平均风荷载叠加，利用Newmark 方法进行时程分析，得到的加速度与实测加速度拟合得很好。

以上的工作主要基于卡尔曼滤波方法对荷载进行反演，需要假设荷载为状态变量，给定其状态演化误差的协方差。不合理的协方差将得不到满意的识别结果。本文提出一种逐步最小二乘算法反演风荷载，它基于精细积分法和最小二乘法构建目标函数，逐个时间步识别风荷载。算例表明，该方法较增广卡尔曼滤波方法精度高。

1 基于Davenport谱的风荷载

任意高度的风速通常可以表述为平均风速v 和脉动风速vf的和［11］，即

同样地，对应的风压w 可以表示成平均风压wˉ和脉动风压wf的叠加，即

其中平均风压表示为

式中w0为建筑物当地基本风压；μs为结构体型系数；μγ为重现期调整系数，模拟作用于高层建筑的风荷载时取1.1；βz为风压高度变化系数。

脉动风为随机过程，通常利用统计的方法进行描述。用于描述风荷载的风速谱有多种，我国现在规范使用的是Davenport 谱。Davenport 风速谱描述的风速功率谱Sv( f)不随高度变化，其表达式为

其中脉动风速的方差为

脉动风压的方差为

因此，脉动风压谱与脉动风速谱关系为

由脉动风压谱通过谐波叠加法模拟脉动风压时程，即

φj为[0，2π]之间的一个随机数，A 为风压作用面积，Δf表示频率步长。

2 精细积分法求解结构动力响应

对于一般的多自由度系统，其动力学方程为

在进行荷载反演之前，需要建立观测值与荷载值的关系。本文采用精细积分法进行数值求解并建立观测加速度与荷载的关系。实际中，在求解结构动力响应时，自由度数目往往会特别大，而真实响应往往只与少数低阶模态相关，因此可以采用振型叠加法对式（10）进行解耦和降维。工程应用中，各节点的质量和刚度难以准确确定，但系统前几阶固有频率、阻尼比、振型可以较为容易地识别。比如，结构的固有频率可以通过对实际测量点的时域数据进行傅里叶变换，获得对应的自功率谱，再通过拾取峰值点来获得；结构振型则可以通过互功率谱与自功率谱的比值确定［10］。取前几阶振型，形成矩阵Φ，对应的频率矩阵为Ω2，根据振型的定义，有

位移yi+1与模态坐标下的位移ui+1满足

根据质量归一化准则，假设阻尼为比例阻尼，有

基于此，振动方程（10）可以简化为

其中对角矩阵Γ 的第i 个对角分量为2ξiωi，ξi表示第i阶模态阻尼比，ωi表示第i阶固有频率，u(t)为模态响应。

以模态坐标下的位移、速度构建状态变量

由振动方程（14）可得连续状态方程

其中矩阵AC和BC为

本文仅观测模态坐标下的加速度时程，观测方程为

其中对应的矩阵G和矩阵J为

由于状态方程（16）为线性定常方程，根据精细积分法，其解可表达为

根据式（20）的结果，构建离散状态方程

其中状态转移矩阵A和B矩阵为

同样，离散的观测方程为

由式（22）和（23）可知，通过第i 时刻的位移、速度和已知的荷载时程可以求解出第i + 1 时刻的位移、速度和加速度。循环这个过程，可以获得所有离散时间点位移、速度和加速度。

3 风荷载反演

荷载反演是典型的反问题，它通过观测动力响应反求荷载时程。本文以模态坐标下的加速度为测量数据，主要考虑两种荷载反演方法-增广卡尔曼滤波方法和逐步最小二乘算法。

3.1 增广卡尔曼滤波方法

基于已推导出的离散状态方程（21），引入满足零均值高斯白噪声假定的模型噪声wi，则可以得到状态方程

令相邻时刻的荷载满足以下关联方程

ηi是一个随机过程，满足零均值高斯白噪声假定。合并式（24）和式（25）的状态向量，得到增广的状态向量

以及增广状态方程为

其中增广系统矩阵Aa和增广噪声向量ξi为

在观测方程中，取vi为第i 时刻的观测噪声，满足零均值高斯白噪声假定，则有增广观测方程

其中增广观测矩阵为

卡尔曼滤波在最小方差无偏估计下定义为最优的线性状态估计器，上文给出的离散状态方程（30）和离散观测方程（29），涉及到的噪声之间互不相关。用矩阵Q、R 和S 分别表示模型噪声方差强度矩阵、观测噪声方差强度矩阵和荷载方差强度矩阵，则应满足以下关系

其中δi-j是克罗内克函数。

增广状态变量的估计误差协方差矩阵Pi|j定义为

增广卡尔曼滤波由5个步骤组成：

1）用第i 时刻的状态向量预测第i + 1 时刻的状态向量

2）用第i 时刻的估计误差协方差矩阵Pi|i预测第i + 1时刻的估计误差协方差矩阵

其中增广状态向量的方差强度阵为

3）计算卡尔曼增益矩阵

4）利用预测的第i + 1 时刻的状态向量和第i时刻的观测加速度修正对状态向量的预测，即

5）更新估计误差协方差矩阵

3.2 逐步最小二乘算法

逐步最小二乘算法针对每一个时间步建立目标函数。首先，对于两个相邻时刻，已知第i 时刻的位移、速度和荷载和第i + 1时刻的观测加速度，反演第i + 1 时刻的位移、速度和荷载。构建的目标函数

式中λ 为权值，其作用是平衡模型误差和观测误差。

对式（40）的xi+1和fi+1求变分，获得线性方程组

其中

式中w 为权重矩阵，本文取为xi+1协方差矩阵的逆矩阵。当系统自由度较大时，Am矩阵往往会具有较高的条件数，可以对式（41）进行1-范数均衡预处理来降低待求逆矩阵的条件数［12］，利用矩阵指数的无穷积分代替其求逆过程，并且可通过精细积分法快速求解无穷积分过程［13］，该方法的精确性远高于利用matlab直接求逆。

4 数值算例

考虑图1 所示的16 自由度的剪切层模型，取各楼层具有相同的质量m1= m2= m3= …= m16=3 500 kg 和相同的刚度k1= k2= k3= …= k16=1500 kN/m。

模型的各阶模态阻尼比均取为2%。荷载的空间分布矩阵为

地面粗糙度指数α 取0.23。第10 层的风荷载Davenport功率谱取为

模拟风荷载的频率范围取为ω ∈[0，15]，参考利用谐波叠加法（9） 构建风荷载时程，具体见图2。

图1 高层建筑简化而来的剪切层模型Fig.1 Simplified shear layer model for high-rise buildings

图2 风荷载时程Fig.2 Wind load time history

基于精细积分法求解加速度时程时，取时间长度为5 s，时间间隔为0.01 s，初始位移和速度均取为0。考虑到测量加速度数据时存在误差，给仿真的加速度添加一定水平的高斯白噪声

图4为基于增广卡尔曼滤波方法反演的风荷载结果，其荷载协方差矩阵S 根据L 曲线方法已经取到最优（见图5）；由于模型精度很高，其模型噪声协方差矩阵Q 取为I × 10-20；一般来说，观测噪声协方差矩阵取决于传感器精度。

图3 模态加速度时程Fig.3 Time history of modal acceleration

图4 增广卡尔曼滤波反演风荷载结果（5%噪声水平）Fig.4 The wind load is retrieved by the augmented Kalman filter(5%noise level)

图5 L曲线法确定最优荷载协方差矩阵SFig.5 L curve method is used to determine the optimal load covariance matrix S

简化起见，本算例的观测噪声协方差矩阵直接由真实加速度和污染加速度计算出来。该方法的识别结果相对误差为10.07%，本文相对误差δ计算公式为

图6为基于最小二乘法反演的风荷载结果，可以看出，仅考虑前四阶的模态加速度数据，在5%噪声干扰下，识别的风荷载很好地逼近真实荷载，识别结果相对误差为5.66%。说明该方法抗噪性能良好。可以看出，本文所提的逐步最小二乘算法比增广卡尔曼滤波方法具有更好的识别精度和抗噪能力。

5 结 论

本文提出了一种逐步最小二乘算法反演风荷载。该方法以模态加速度为观测数据，针对每一个时间步建立目标函数，逐步识别荷载。与已有的增广卡尔曼滤波方法相比，所提方法不需要荷载的协方差信息。算例表明，两种方法在5%噪声水平下，仍然获得不错的识别结果，说明两种方法的抗噪性能都较好，但逐步最小二乘算法的识别精度优于增广卡尔曼滤波方法。

图6 逐步最小二乘法反演风荷载时程结果（5%噪声水平）Fig.6 Stepwise least square method is used to invert wind load time history(5%noise level)