基于贝叶斯估计的动载荷识别方法

雷 博，郑 辉，2

（1.上海交通大学 振动、冲击、噪声研究所，上海 200240；2.高新船舶与深海开发装备协同创新中心，上海 200240）

载荷信息的获取是系统动力学响应计算、动态设计和故障分析的重要组成部分，在实际工程中具有重要的意义。但由于工程环境中存在诸多限制条件，使得多数情况下，载荷无法直接测量[1]。因此，建立正确有效的载荷识别方法是进行后续动力学分析和动态设计的前提。

载荷识别是指利用系统响应和系统特性对系统所受载荷进行求解，属于动力学第二类反问题。Thite等提出利用直接求逆法和最小二乘正则化方法进行设备噪声传递路径上的载荷识别，后者利用正则化方法解决了传递矩阵的病态性问题[2–3]。杨立娟等利用Tikhonov正则化方法进行悬臂壳载荷识别，并通过混合传递路径分析确定了该结构的主要振动传递路径[4]。张方等推导了基于时间有限元的动载荷识别模型，将时域动载荷识别的逆卷积关系转变为广义正交域的线性算子逆运算，大大降低了该类问题的复杂程度[5]。随着计算机技术的发展，载荷识别也可通过智能算法进行。智能算法可避免传递矩阵的病态性问题，同时可降低运算量。如Yan等利用微遗传算法对复合加筋板所受冲击载荷进行了识别，结果表明该算法具有较高的计算效率[6]。沙瑞华利用三种神经网络算法对水电机组所受载荷进行了识别，发现基于LM优化算法的BP网络对各类动载荷识别的效果优于附加动量法BP网络和RBF网络[7]。

由此可见，对动载荷识别的研究已取得很多成果，然而考虑动载荷识别过程中由于测量误差等不确定因素所带来识别误差的研究尚不多见。为了降低不确定性因素带来的动载荷识别误差，本文从系统运动方程出发，考虑测量误差这一不确定性因素，利用贝叶斯估计理论，结合马尔可夫蒙特卡罗方法，建立了基于贝叶斯估计的动载荷识别方法，并通过仿真和实验数据对比对该方法的有效性进行了验证。

1 载荷和误差参数的概率分布

1.1 贝叶斯估计理论[8–9]

传统的统计推断方法利用统计量的先验信息进行推断，而贝叶斯估计理论同时考虑统计量的先验信息和观测数据中包含的有效信息。

将系统所受载荷看作系统参数，对系统参数进行估计的贝叶斯公式为

式中M表示给定的结构模型，D表示观测数据，θ表示系统参数，p(θ|D,M)为给定模型和观测数据下参数的后验概率分布，p(D|θ,M)为给定模型在输入D下参数的似然函数（观测数据信息），p(θ|M)为给定模型下参数的先验概率分布（先验信息），通常由经验获得，p(D|M)为观测数据的边缘概率分布，与系统参数无关，起正则化作用，其计算公式为

贝叶斯估计的实质是在利用观测数据对系统参数的先验概率分布进行修正，进而得到系统参数的后验概率分布，以此作为估计推断的依据。

1.2 概率分布估计

系统的运动方程为

式中Fm=[F1，F2，…Fm]T表示系统所受载荷的频域列向量，Yn表示系统响应的频域列向量，Hn×m表示系统传递矩阵，e表示由不确定因素引起的预测误差，服从均值为零、协方差矩阵为σe2In×n的正态分布。

假设系统所受载荷的先验概率分布分别为均匀分布，即

进一步假设误差参数σe2的先验概率分布为倒伽马分布IG(α，β)，即

式中α为形状参数，β为尺度参数。

将系统所受载荷看作系统参数，由式（2）、式（4）、式（5）、式（6）可得系统所受载荷和误差参数的联合后验分布为

对传递矩阵进行奇异值分解为

式中，上标H表示共轭转置。

由式（7）、式（8）可得系统所受载荷和误差参数的边缘概率分布分别为

2 载荷抽样估计

2.1 马尔可夫蒙特卡罗方法[10]

马尔科夫蒙特卡罗方法是根据统计抽样理论获取近似解的一种统计试验方法，即构造平稳分布为目标概率分布的马尔科夫链，通过随机抽样获得容量足够大的样本。根据大数定律，每个参数的样本值分布近似于其边缘概率分布，因此以该样本的统计特征（期望）作为满足目标概率分布的参数估计值。在本文建立的方法中，其目标概率分布为载荷和误差参数的联合后验分布。

2.2 Gibbs抽样算法[11]

马尔科夫蒙特卡罗方法常用的抽样算法为Metropolis-Hastings(MH)抽样算法和Gibbs 抽样算法。

MH抽样算法在每次抽取样本后需要计算接受概率，根据接受概率的大小判断是否接受本次抽取的样本，若接受，则在本次抽取样本的基础上继续抽取，若不接受，则需要重新抽取，直至接受，抽样效率较低。而Gibbs抽样算法可看作MH抽样算法中接受概率为1的特例，即每次抽取的样本均被接受，大大提高了抽样效率。

构造目标概率分布为载荷和误差参数的联合后验概率分布的马尔可夫链，并通过Gibbs抽样算法进行抽样的具体流程如图1所示。

图1 Gibbs抽样算法流程图

3 仿真和实验验证

3.1 仿真算例

在基于本文方法进行载荷识别时，需要动力学系统的传递矩阵和输出响应。仿真算例以两端简支梁为例，通过ANSYS有限元仿真，获得相应的传递矩阵和输出响应。该梁长度为1 m，密度为2 700 kg/m3，弹性模量为71 GPa，泊松比为0.3，结构阻尼损耗因子为0.002，前3阶固有频率分别为118.95 Hz、473.8 Hz、1 058.8 Hz，仿真中载荷（F1、F2）和响应参考点（1、2、3）的信息见表1，表中位置均表示与简支梁左端的距离。

表1 载荷幅值、位置及响应参考点位置

前两次在不同位置施加单位激励力，进行谐响应分析，获得传递矩阵，而第三次施加一定大小的力，获得系统加速度响应。

考虑不确定性因素的影响，在仿真获得的传递矩阵和加速度响应中均加入信噪比为40 dB的高斯白噪声，噪声模型为

式中Nr取分布为N(0，A(ω)10-(B/20))的正态分布随机数，A(ω)表示响应或传递函数，B表示信噪比，Nu取0到1上的均匀分布随机数。

利用添加噪声后的传递矩阵和加速度响应，假设误差参数的先验分布为IG（3,1），分别通过Gibbs抽样算法和基于奇异值分解的载荷识别方法[12]，对第三次仿真施加的载荷进行识别，结果如图2所示。

以载荷频率1 200 HZ时为例，使用Gibbs抽样算法所得载荷的抽样曲线如图3所示，样本直方图如图4所示。

对比图2的结果可知，在大部分频段内，两种方法的识别结果均较好，而在梁的固有频率附近，通过Gibbs抽样算法获得的识别结果较奇异值分解法的识别结果更接近真实值，说明该方法能够降低不确定因素对识别结果的影响。

将梁模型的边界条件改为两端固支，其余参数均不变，使用Gibbs抽样算法所得载荷识别结果如图5所示。

图2 简支梁所受载荷的识别结果

图3 1 200 Hz时简支梁所受载荷抽样曲线

图4 1 200 Hz时简支梁所受载荷样本直方图

图5 固支梁所受载荷识别结果

结果表明，识别精度与两端简支边界条件下类似，在大部分频段内，载荷识别结果均较好，仅在固有频率734.41 Hz和1 425.6 Hz附近具有较大误差。因此，在固支边界条件下，本文方法同样适用。

3.2 实验验证

由于在工程应用中，系统结构复杂，且存在加工误差、安装误差等多种因素的影响，为了验证该方法在工程应用中的适用性，故实验对象采用结构相对复杂的钢制夹层梁，且安装条件为两端固支。该梁长0.468 m，宽0.03 m，高0.024 m，实验台架见图6，从左至右，加速度传感器分别安装在响应参考点1至6。

使用激振器施加一定大小的激励力，稳定后，通过LMS数采系统采集激励力信号和加速度响应信号。

图6 实验台架

选用参考点1、3、6作为本次载荷辨识的响应参考点，建立该结构的有限元模型，通过仿真所得各响应参考点的传递函数如图7所示。

图7 响应参考点传递函数

利用仿真所得的传递函数和测量所得各响应参考点的加速度响应，通过Gibbs抽样算法计算结构所受的激励力，计算所得的激励力和测量所得的激励力如图8所示。

图8 激励力实测和计算结果

从图8可以看出，基于贝叶斯估计理论，利用Gibbs抽样算法进行载荷识别，识别结果与实测结果基本相符，能够满足工程实际中的识别精度要求。

4 结语

本文从系统运动方程出发，考虑不确定性因素的影响，对载荷和误差参数服从的先验概率分布进行了假设，从而推导了载荷和误差参数的后验联合概率分布，并分别得到二者的后验边缘概率分布。基于二者的后验边缘概率分布，通过马尔科夫蒙特卡罗方法，采用Gibbs抽样算法对系统所受载荷进行了估计，并通过仿真算例和实验数据对本方法进行了验证。结果表明，该方法满足工程实际的要求，且在一定程度上降低了不确定性因素对载荷识别的影响，对提高载荷识别精度具有参考意义。

[1]胡寅寅，率志君，李玩幽，等.设备载荷识别与激励源特性的研究现状[J].噪声与振动控制，2011，31(4)：1-5.

[2]THITE A N,THOMPSON D J.The quantification of structure-borne transmission paths by inverse methods.Part 1:Improved singular value rejection methods[J].Journal of Sound&Vibration,2003,264(2):411-431.

[3]CHOI H G,THITE A N,THOMPSON D J.Comparison of methods for parameter selection in Tikhonov regularization with application to inverse force determination[J].Journal of Sound&Vibration,2007,304(3-5):894-917.

[4]杨立娟，罗新杰，苗慧慧，等.混合传递路径分析方法[J].噪声与振动控制，2016，36(6)：12-15.

[5]张方，朱德懋.动载荷识别的时间有限元模型理论及其应用[J].振动与冲击，1998，17(2)：1-4.

[6]YAN G,ZHOU L.Impact load identification of composite structure using genetic algorithms[J].Journal of Sound&Vibration,2009,319(3-5):869-884.

[7]沙瑞华.基于神经网络的水电机组动载识别研究[D].大连：大连理工大学，2005.

[8]BOLSTAD W M.Introduction to Bayesian statistics,second edition[M].2007.

[9]BERNARDO J M,SMITH A F M.Bayesian theory[M].2008.

[10]曹晨.基于MCMC方法的统计模型的参数估计[D].南京：南京航空航天大学，2007.

[11]刘书奎，吴子燕，张玉兵，等.基于Gibbs抽样的马尔科夫蒙特卡罗方法在结构物理参数识别及损伤定位中的研究[J].振动与冲击，2011，30(10)：203-207.

[12]郭荣，房怀庆，裘剡，等.基于Tikhonov正则化及奇异值分解的载荷识别方法[J].振动与冲击，2014，33(6)：53-58.