基于粒子滤波算法的风力发电塔地震动力响应预测

徐亚洲， 任倩倩， 于明阳， 时文浩

(西安建筑科技大学 土木工程学院，西安 710055)

随着电子计算机技术的发展，有限元数值分析模型和结构动力试验成为结构动力学中愈加成熟的领域，因此动力学模型更新和对动态系统结构性能的合理预测仍是非常重要的研究方向[1-2]。利用测量系统响应进行动力模型更新，在结构振动控制、结构健康监测以及可靠性评估等方面都有广泛的应用，但在建立结构计算模型及其激励的数学模型的过程中，可能存在由于外部原因导致的建模误差及参数不确定性影响，因此能够准确和适当地量化此类不确定性对于结构动力模型的更新至关重要[3-5]。土木工程结构中的有限元模型往往与许多来源的不确定性和建模误差有关，尤其是对于复杂的土木工程结构，在建模过程中有限元模型相比实际结构进行了大量的假设和简化，当此类有限元模型用于预测一些感兴趣的量时，这些不确定性因素会影响土木工程结构的设计和性能评估[6-8]。这些不确定性因素促使研究人员通过概率有限元模型更新方法来考虑潜在的结构不确定性。

基于经典概率的理论分析模型是目前工程中常用的不确定性分析模型[9]，但在统计信息不充分的条件下，经典概率分析往往具有一定的局限性。近年来，贝叶斯估计理论由于克服了经典概率理论中难以处理的小样本问题成为概率不确定性分析研究的新热点[9]，贝叶斯的基本思路是根据状态的先验信息和测量信息来不断递推更新状态的后验概率密度函数[10]。毛文贵等[11]运用基于遗传智能采样技术的改进贝叶斯理论，可以很好地处理轴承转子系统不平衡识别中的模型不确定性参数的影响。Pepi等[12]建立了位于特尔尼的弯曲斜拉桥有限元模型，该有限元计算模型考虑了模型不确定性以及对连接和构件的有限认识，在贝叶斯框架下进行环境振动试验以确定用于模型更新的主要动态参数，并用于评价基于有限元模型计算结果的准确性。Sedehi等[13]在贝叶斯层次模型的基础上，提出了渐近近似和最大后验估计方法，该方法可以有效解决由外界因素引起的参数不确定性的影响。

在贝叶斯理论估计计算框架中，由于粒子滤波(particle fliter,PF)方法具有较好的处理复杂系统的能力而被广泛研究及应用[14]。万春风等[15-17]基于粒子滤波算法对单自由度和多自由度结构进行参数损伤识别研究，指出了粒子滤波算法在参数识别方面相较于其他滤波方法的优越性。一些学者[18-21]采用粒子滤波方法在线预测疲劳裂纹扩展情况并更新剩余寿命，并通过试验结果验证了粒子滤波对疲劳裂纹扩展和剩余寿命在线预测的可行性和合理性。樊学平等[22-23]采用粒子滤波算法建立了随时间更新的动态模型，并结合极值监测数据实现桥梁结构可靠度指标的动态预测。

风力发电塔对动力荷载的高敏感性引起了研究者的关注[24-26]，且目前关于风力发电塔的抗震性能研究多采用有限元模拟的手段，但在构造结构数值分析模型的过程中，总是存在简化模型所带来的建模误差和参数不确定性。振动台试验作为目前获得结构动力响应分析的主要手段，被广泛应用于结构性能评价与分析中，试验结果对于校正能够准确描述结构响应的有限元数值分析模型至关重要。因此本文提出了一个用于结构动力响应预测的概率贝叶斯估计计算框架，即首先在贝叶斯估计原理基础上，建立地震激励下结构有限元分析模型。然后采用粒子滤波方法，充分利用结构先验信息并结合振动台试验观测的风力发电塔动力响应数据信息选择重要性函数，实现风力发电塔动力响应预测，并以有限元数值解与试验值之间的差值来量化有限元模型的不确定性误差。最后通过风力发电塔振动台试验结果，对计算方法的合理性及有效性进行验证。

1 粒子滤波算法

在存在结构建模不确定性的情况下，需要确定并建立合理的结构模型，并基于实测结果进行加权修正。粒子滤波算法具有显著处理参数不确定性估计的优越性，目前被广泛用于结构损伤识别和结构不确定性分析中。该算法主要是应用蒙特卡洛方法产生大量的随机样本(粒子)，然后结合观测值利用样本权重预测下一时刻的系统状态或参数概率分布，最后根据算法不断递推更新得到结构参数的最优贝叶斯估计值。

与其他滤波方法一样，粒子滤波算法需要基于估计问题构建状态空间模型，主要由状态方程、观测方程和初始状态信息组成。状态方程描述了系统状态量随时间的演变过程，观测方程描述了系统适时输出与系统当时状态的模型，函数表达式如下

状态方程为

Xk=f(Xk-1)+ωk

(1)

观测方程为

Zk=h(Xk)+νk

(2)

初始状态信息为

p(X0Z0)=p(X0)

(3)

式中:Xk和Xk-1分别为k时刻和k-1时刻的状态值；Zk为k时刻的观测值；f(·)为表达系统当前时刻与前一时刻的变化关系的已知非线性函数；h(·)为状态值和观测值之间的关系非线性函数；wk为过程噪声，描述了系统的各种不确定性影响因素；vk为观测噪声，代表了观测过程中的不确定性观测误差。

1.1 贝叶斯估计理论

贝叶斯的基本步骤包括预测和更新[27]，预测为利用模型前一时刻的先验已知信息对当前时刻进行预测，更新是结合当前时刻的实际观测进行递推，从而获得待估值在当前时刻的后验概率密度，在贝叶斯更新过程中，状态变量的更新是为了更好地表示真实的结构，其中更新过程是由先验已知信息和所研究的结构测量中包含的测量信息完成的。状态量中包含了所有的不确定性，因此，假设前一时刻的后验概率p(Xk-1Z1:k-1)已知，则基于贝叶斯估计算法，可得到后验概率密度p(XkZ1:k)的估计式

状态预测方程

(4)

状态更新方程

(5)

式(4)和式(5)组成了贝叶斯估计的计算框架。首先若已知k-1时刻的状态值Xk-1，则可通过式(1)预测得到k时刻Xk的先验概率密度p(Xk-1Z1:k-1)，基于观测方程(式(2))可得到似然概率p(ZkX1:k)。在贝叶斯计算中，直接计算后验概率密度函数牵涉到复杂的积分计算，对于非线性动态系统的状态估计问题求解困难。因此，引入了基于随机抽样计算的蒙特卡洛方法来计算后验概率。

1.2 蒙特卡洛原理

基于大数定理，蒙特卡洛方法[28]将所求解的随机问题发生的概率用大量试验中发生的频率来估计，当统计数据比较大时，可以认为该随机事件发生的频率能无限接近期望值，从而可使贝叶斯后验概率中的积分运算用大量样本点的平均回报来逼近。

如对本研究来说，可以在状态的后验概率密度中大量的随机采样，则状态的后验概率密度可表达如式(6)。作为一种复杂问题有效的处理求解方式，该方法回避了工程应用中积分计算的困难，用状态发生的频率来表示概率，这在积分运算中有着非常重要的意义。

(6)

式中，δ(·)为狄拉克函数。然后待估状态的期望值计算如式(7)所示

(7)

1.3 粒子滤波算法

(8)

作为本文研究问题所要求解的对象，故不能在待求解的后验概率中采样。一般采取从已知重要性密度函数q(XkZ1:k)中选取的N个样本近似表示

(9)

其中,ω(Xk)为样本权重，即

(10)

依据式(10)，式(9)可进一步化简为

(11)

其中，归一化后的权重值为

(12)

基于马尔可夫假设与简化，归一化后粒子的权值(重要性权值)在每个时间步长上有如下递推关系

(13)

通过权重值的不断计算，这些权重被附加到每个相应的粒子上。但是这些粒子并没有被修改，只是它们的相对权重被改变了，也就是说，当所有的粒子偏离观测值时它们不会被拉回观测值，重要性采样只改变了它们的相对权重。这个缺点可能会导致单个粒子的权重占比非常大，而其他粒子的权值很小。因此，需要进行重采样步骤，定义有效粒子数量为

(14)

当有效粒子的数量小于采样数量N时，则表示粒子样本退化严重，需要进行重采样，即高重要性权值的粒子会被多次复制，低权值的粒子少量复制或被去除以保持粒子数量不变。重采样后每个后验粒子的权值相同，均为1/N。

重要性密度函数一般选取为状态转移密度函数，即

(15)

将式(15)代入式(13)，可重新得到权值更新方程的表达式

(16)

2 动力响应预测的粒子滤波实现

建立有限元模型进行风力发电塔动力响应分析的过程中，由于理想假设或简化容易导致参数不确定性或模型不确定性，继而导致有限元模拟值与实际测量值之间存在较大区别。因此本章基于第1章关于粒子滤波算法的计算原理的介绍，利用粒子滤波算法中的状态空间模型将这些不确定信息考虑进来，作为先验已知信息实现结构动力响应的预测过程。

2.1 力学模型

对于建筑结构，常用的力学模型分为层模型、杆模型及杆系-层间剪切模型。根据风力发电塔有限元分析方法的需要，目前有4种比较常见的建模方式。第一种为IEC规范[29]中所建议的单自由度模型简化方法，即将塔体顶部叶轮、机舱的质量及塔筒50%的质量集中到塔筒顶部；第二种为我国GB 50011—2010《建筑抗震设计规范》[30]中常用的建议多自由度模型简化方法，即将塔体顶部叶轮、机舱的质量集中到塔体顶部质量点处，各个筒段的质量分别集中到对应质量点处；第三种为简化的有限元模型[31]，即叶轮、机舱采用刚体单元，将质点位置按照刚体运动耦合到塔筒顶部中心；第四种为精细化有限元模型[32]，即将叶片按照实际结构原尺寸建模从而考虑叶片与塔体相互作用对地震作用下结构性能的影响，每个叶片单元的质量集中在一个参考点，轮毂和机舱用集中质量点表示，通过刚性连接连接到塔顶部分。本节采用的第二种建模方法，即将风力发电塔简化为多自由度层剪切结构进行分析。

以剪切型动力学模型为例，地震激励下其动力运动方程为

(17)

结构的阻尼采用瑞丽阻尼，即

C=αM+βK

(18)

式中，α，β为阻尼系数。

2.2 动力方程求解

在结构动力方程的求解过程中，Newmark-β法应用较为广泛，其主要思想是假定ti～ti+1时间段内的加速度变化呈线性规律，其假定为

(19)

式中，a，b为可调整参数。当a=1/2，b=1/6时，即为线性加速度法；当a=1/2，b=1/4时，即为平均常加速度法；当a=1/2，b=0时，即为中心差分法。本节采用线性加速度法计算。由式(19)可求得t+Δt时刻的速度和加速度

(20)

其增量表示形式为

(21)

将式(21)代入式(17)，可得到

(22)

通过求解式(22)可得到t+Δt时刻的位移，将其代入式(20)即可求得结构加速度及速度响应。

2.3 结构动力响应预测实现

基于试验观测数据，粒子滤波算法的结构动力响应预测过程的主要步骤如下。

步骤1首先需要构造状态空间模型，根据Newmark-β方法求得有限元简化模型结构的动力响应数值解作为趋势项，用来构造本文的状态方程；使用监测仪器测得结构的实际动力响应，然后利用状态量及振动台试验动力响应观测量的信息构造观测方程；初始状态先验信息基于有限元数值计算结果进行概率统计给出。

步骤2基于状态空间模型，根据初始先验信息的概率分布函数，利用蒙特卡洛抽样方法可得到状态变量的初始参数样本，即初始先验粒子{x0,i,i=1,2,3,…,N}。

步骤3在每一个时间步长上进行粒子更新，利用状态方程可获得下一个时刻的后验粒子

(i=1,2,3,…,N)

步骤4基于k时刻的实际测量值Zk，可通过似然概率密度函数更新每个后验粒子的权值，此权值衡量了该粒子产生的模拟数据接近实际观测数据的近似程度。似然概率密度函数可近似为正态分布[33]，即

(23)

式中:Zmn,i为第i个粒子的模拟观测数据；R为量测噪声方差，可通过状态量及观测值的数据信息近似估计。

步骤5根据式(12)进行权重归一化。

步骤6根据式(14)计算有效粒子的数量，当有效粒子的数量小于初始采样数量N时，则进行重采样，重采样后每个粒子的权值相同，均为1/N。

以此不断递推更新，求得结构响应的最优预测值。

具体流程图如图1所示。

图1 粒子滤波算法预测过程流程图

3 试验验证

本文开展了2 MW陆上风力发电塔振动台试验，基于状态空间模型和试验观测值对有限元计算结果的不确定性进行修正，说明本文所提用于结构动力响应预测的粒子滤波框架的合理性。

3.1 试验模型信息

根据量纲协调原理和结构动力方程[34-35]确定相似关系，按照相似原则及抗弯刚度等效设计模型。试件与原型结构的几何相似比为1/20，模型材料为Q345钢。塔呈锥形，直径从底部的200 mm线性减小到顶部的130 mm。塔筒为4段，段与段之间用10 mm厚法兰板连接，每段高度从下到上分别为1 m,0.85 m,1 m,1 m。由于薄板在轧制和焊接过程中制作困难，根据抗弯刚度等效原则，4个管段对应的管壁厚度分别为4 mm,3 mm,3 mm和3 mm。机舱和转子的总质量为157.5 kg，4段塔筒的质量分别为85 kg，58 kg，51 kg，35 kg。试验模型的示意图如图2所示。

按照抗震设计规范中输入地震波的选取原则，从太平洋地震工程研究中心[36]数据库选取了Chi-Chi波、Westmorland波、El-Centro波和Taft波作为模型结构振动台输入地震激励，将垂直于叶片的方向定义为地震激励输入的方向，图3给出了4条地震动记录水平方向加速度时程。将地震动加速度峰值(peak ground acceleration，PGA)统一调幅至0.07g，其加速度反应谱如图4所示。

(a) 试验模型

(a)

(c)

地震波的输入工况分别设置为0.07g，0.2g和0.4g。在不同幅值地震动输入前后采用PGA为0.03g的白噪声激励来识别试验模型的动力特性。本试验加速度相似比为2，频率相似比为6.324 6，故实际台面输入地震激励需按照相似关系进行调整。

在风力发电塔结构抗震设计中，塔顶的动力响应是抗震需求的重要指标，因此通过合理布置测点测量了塔顶加速度响应及位移响应。不同地震动强度幅值输入前后利用白噪声扫频，并对结果进行处理得到模型结构的动力特性。利用传递函数对试验模型的固有频率和振型进行识别，得到结构前两阶自振频率分别为3.22 Hz,22.46 Hz，结构一阶阻尼比为1.5%。

图4 地震动记录的加速度反应谱

3.2 模型动力响应预测

基于本文方法，对模型结构动力响应进行预测。

有限元数值分析：基于2.1节的介绍，本研究将风力发电塔动力系统简化为5个自由度的剪切动力学模型来描述，如图2(b)所示，根据质量点处的截面尺寸计算出截面抗弯刚度。基于MATLAB有限元软件，计算5个自由度简化模型的质量和抗弯刚度矩阵。根据结构实测的阻尼比及自振频率，采用瑞丽阻尼计算阻尼矩阵，最后利用2.2节中介绍的Newmark-β法求得有限元模型的塔顶加速度响应和位移响应。

状态空间方程建立：根据12个地震工况作用下求得的5个自由度结构的动力响应数值计算结果，选用塔顶加速度响应和位移响应为状态待估量，可估计出每个地震工况作用下的状态噪声方差；利用数值计算结果及对应工况下的振动台试验动力响应观测量信息估计出测量噪声方差；基于数值计算结果进行概率统计给出状态量的初始信息，并基于初始状态信息选取200个粒子。利用第2章中介绍的粒子滤波算法框架进行风力发电塔动力响应预测。

以振动台试验测得的观测值为参考，基于粒子滤波算法对有限元模型结果进行修正预测。定义考虑和不考虑观测值修正后的计算值与试验实测值之间差值的绝对值作为动力响应的计算偏差。则考虑观测值修正后(PF)，和不考虑观测值修正(finite element method,FEM)的加速度和位移响应的偏差最大值计算结果，如图5所示。结果表明有限元计算的动力响应偏差随地震动PGA幅值的增大而增大。以Chi-Chi地震动激励下塔顶加速度响应为例，当PGA值为140 gal时，加速度偏差为2.21 m/s2，而PGA值为800 gal时，加速度偏差为14.89 m/s2。考虑观测值修正之后有限元计算结果偏差显著减小，以塔顶位移响应为例，考虑观测值修正之后位移偏差减小近一倍。

(a) 塔顶加速度响应偏差

(b) 塔顶位移响应偏差

为进一步说明本文方法的有效性及可行性，选取了El-Centro波和Chi-Chi波激励作用下剪切动力学模型结构的加速度响应时程进行分析，图6给出了粒子滤波算法修正前后加速度时程曲线对比。同一地震激励下，随输入地震动幅值的增大，有限元计算结果相对于试验值的偏差逐渐增大。考虑观测值修正后粒子滤波预测的结果精度较好，预测值与试验实测值近似相等，具有很好的一致性，进一步验证了本文所提校正不确定有限元分析模型的概率贝叶斯框架和粒子滤波预测算法的合理性，表明粒子滤波算法可以利用有限元分析结果不确定性的所有先验信息，从而实现对包含变异性结构动力响应的稳步预测。

4 结 论

本研究基于粒子滤波算法，提出了一个用于校正不确定有限元数值模型的概率贝叶斯计算框架对有限元计算的动力响应结果进行预测。利用某2 MW风力发电塔缩尺模型的振动台试验验证了此框架的有效性和可行性。从分析结果可以得到以下结论：

(1) 有限元计算结果具有显著的不确定性，且其计算误差随地震输入幅值的增大而增大，考虑观测值修正之后预测结果不确定性显著减小。

(a)

(d)

(2) 粒子滤波模拟预测的动力响应值精度较高，模拟预测值与试验实测值的时程曲线近似相等，变化规律具有很好的一致性。

(3) 粒子滤波算法与试验观测手段相结合，可以完成高效的结构动力响应预测，继而为结构不确定性及可靠度分析提供参考。

结构不确定性分析应成为工程结构分析的一个重要部分，粒子滤波算法可以利用有限元分析结果不确定性的所有先验信息，从而实现对结构动力响应的稳步预测。