基于Gibbs抽样的结构时域载荷识别

王 婷， 万志敏， 郑伟光

(1. 南通职业大学 机械工程学院，江苏 南通 226000； 2. 华中科技大学 机械科学与工程学院，武汉 430074)

在许多实际工程问题中，精确地掌握作用于结构上的外载荷对结构优化设计、响应重构、故障诊断以及健康监测等起着非常重要的作用。然而，很多情况下的载荷根本无法直接测量，而结构响应却往往较易获得，因此载荷识别技术得以发展，主要包括频域法和时域法两大类[1]。频域法首先获得频域载荷，再通过傅里叶变换获得时域载荷，故频域法对冲击载荷识别效果不佳。时域法[2-5]直接在时域内识别载荷的时间历程，故不受载荷形式限制，且识别结果具有明确的物理意义，无变换误差，在工程中具有良好的应用前景。

另一方面，载荷识别是动力学的第二类反问题，其中的结构矩阵求逆过程往往是病态的。正则化技术被国内外学者广泛用于获得该病态问题的稳定解，而其中的关键在于正则化参数的选取。Law等[6]基于时域法识别了移动载荷，并引入Tikhonov正则化技术寻求近似解，其中正则化参数采用L-曲线法获得。Mao 等[7]建立了基于精细积分的载荷识别模型，并采用GCV准则获得了正则化参数。马超等[8-9]分别采用一种基于局部最小准则及改进的正则化技术识别了载荷。刘杰[10]基于Green核函数建立时域载荷识别模型，并比较了不同正则化参数的选择方法，包括L-曲线、GCV准则等。

然而，上述传统方法获得的正则化参数均为恒定值，往往很难适用于不同类型载荷的反演。Zhang等[11]提出了贝叶斯框架下的载荷识别频域方法，该法具有本征的自适应正则化性能。本文将该法拓展到载荷识别的时域法，将载荷与测量噪声看作是随机变量，建立载荷识别的贝叶斯模型，并采用Gibbs抽样法获得了载荷后验值。对照Tikhonov正则化的数学程式，可知本文提出的载荷识别法具有可变的正则化参数。数值算例分别采用六自由度弹簧质量模型与四边简支板模型为对象，验证了本文方法的有效性。

1 载荷识别的状态空间模型

对于线性时不变动态系统，其离散状态空间方程可以表示为

xt+1=Axt+Bft

(1)

yt=Rxt+Dft

(2)

式中：t∈[1,nt]表示离散时刻；xt为状态空间变量；ft为外载荷向量；A,B,R和D分别为系统矩阵、输入矩阵、输出影响矩阵及输出传递矩阵；yt为响应向量。

假设初始状态x1为零， 同时令H1=D,Ht=R(A)nt-1B， 将式(1)代入式(2)中，则可得矩阵形式方程如下

Y=Hf

(3)

有H=

2 基于贝叶斯理论的动载荷识别

由于实测的响应都是存在噪声的，故实际载荷识别方程应为

Y=Hf+W

(4)

式中：W为噪声向量。

本文采用贝叶斯思想来求解上述动载荷识别问题，即通过后验概率密度分布来定量未知载荷时间历程。首先定义向量θ为不确定性参数，包括载荷大小、测量噪声，那么根据贝叶斯原理，联合概率后验概率密度函数可以写成

p(θ|Y)=p(Y|θ)p(θ)/p(Y)

(5)

2.1 贝叶斯框架下的载荷识别

根据中心极限定理，测量噪声可以近似表示为均值为0的正态分布。例如，第k个响应的噪声分布为

(6)

式中： 符号～表示“分布为”。那么似然函数可以写成

(7)

假设时域载荷向量为共轭先验分布，则概率分布为

f～N(f0,Cf0)

(8)

(9)

(10)

(11)

式中： 超参数{kW,αW,U0,σU0,kf,αf}均假设已知。结合上述超先验分布，得到了一个三层贝叶斯载荷识别模型：其一为似然函数层；其二为载荷先验模型层；其三为底部超参数层。由于分层贝叶斯模型对于超参数的取值非常稳定，故对于上述共轭先验中的超参数可以取任意值。

2.2 Gibbs抽样法求解不确定参数

根据上述的似然函数以及参数先验分布，可以得到所有不确定性参数的联合后验概率密度分布如下：

(12)

其中先验分布

(13)

式(13)是由于各个参数之间相互独立。另外，由于本文采用的先验分布以及似然函数都为标准分布(正态与伽马分布)，那么各个参数的条件概率密度函数分别可以推导出为标准分布形式，具体推导过程如下：

(1) 载荷的条件概率密度函数

(14)

对式(14)右边取对数，则

(15)

对式(15)取载荷f的一阶偏导求极值，得到极值点为

(16)

对式(15)取载荷f的二阶偏导求极值，得到极值点为

(17)

(18)

(2) 测量噪声与模型误差方差的条件概率密度函数

(19)

(3) 载荷先验均值的条件概率密度函数

(20)

(4) 载荷先验方差的条件概率密度函数

(21)

学术界现有众多的MCMC(Markov Chain Monte Carlo)算法[12-14]，由于本文推导的所有条件后验概率分布均为标准分布，如正态分布、伽马分布，故满足Gibbs抽样法的应用条件。综上，基于Gibbs抽样的载荷识别算法总结如下：

⑥返回第②步，重复迭代抽样，直到抽样到足够多的粒子。

(22)

(23)

2.3 贝叶斯方法与正则化方法比较

由上述贝叶斯方法可得结构动载荷的后验解为

(24)

可以看出式(24)与传统正则化方法有些许联系。标准Tikhonov法的实质是求解满足最小均方差的载荷，

(25)

即动载荷的正则化解形式如下：

(26)

3 数值分析

采用两种结构模型来验证本文方法的可行性，并与Tikhonov正则化方法进行对比。两种仿真结构模型分别为六自由度弹簧质量结构模型与四边简支板结构模型。

3.1 算例 1

六自由度的弹簧质量结构，如图1所示，共包含6个质量单元和11个弹簧单元，其质量和刚度参数如表1。假设结构为Rayleigh阻尼，结构的6阶模态频率分别为9.76 Hz, 53.26 Hz, 72.54 Hz, 83.51 Hz, 85.58 Hz和162.28 Hz。

图1 弹簧质量结构模型Fig.1 Mass-spring structural model

表1 结构的质量刚度值

载荷作用于质量块2上，采用如下形式

F1(t)=50(1-cos(20πt))sin(45πt)

(27)

测量响应采用质量块3、5的加速度响应，并将噪声水平为5%的测量噪声加入测量信号中。采样频率采用1 000 Hz，仿真时间为0.3 s。此时，Markov矩阵H的条件数为1.284 5E+008，非常大的条件数意味着矩阵呈病态特性。

针对上述载荷反演不适定性问题，采用Tikhonov正则化技术消除矩阵奇异性，获得载荷的正则化稳定解。首先，采用L-曲线准则获得了正则化参数为λ=2.508 6E-005。进而把正则化参数值代入Tikhonov方程(26)计算出载荷识别值，如图2所示。此时虽然载荷识别值不发散，但与真实值的误差仍然较大，说明正则化参数的选取不合适。

然后，采用另一种常用的正则化参数选择准则——GCV准则，来选取正则化参数。GCV选择参数值为λ=0.007 020 7。采用该参数，代入Tikhonov方程(26)计算出载荷识别值，如图3所示。由图可以看出，载荷识别曲线基本和真实曲线一致，整体误差较小，仅曲线拐点处误差较大，总体正则化效果要比上述L-曲线正则化方法好得多。另外，对比上述两种准则选取的参数值，发现通过GCV准则选择的参数值要远大于通过L-曲线准则选取的参数值，说明L-曲线选取的参数值偏小，导致矩阵仍然病态。

图2 载荷识别值与真实值(L-曲线)Fig.2 The actual and identified values of the load (L-curve)

图3 载荷识别值与真实值(GCV准则)Fig.3 The actual and identified values of the load (GCV criterion)

图4 载荷识别值与真实值(Gibbs抽样)Fig.4 The actual and identified values of the load (Gibbs sampling)

图5 三种正则化参数对比图Fig.5 The comparison results of the regularization parameter based on these three methods

3.2 算例 2

采用如图6所示的简支板为验证对象，力锤冲击载荷作用于薄板中心(325 mm，225 mm)处，采用两个加速度信号(位置如图6，黑色方框代表加速度传感器)作为测量信号来识别该冲击载荷。

不采用正则化技术时，冲击载荷识别曲线如图7所示，可以看出载荷识别曲线从开始阶段到结束阶段逐渐发散，表明该反求系统病态，故采用Tikhonov正则化技术来获得其稳定解。采用了L-曲线准则和GCV准则来获取正则化参数分别为0.014 07、0.890 1，那么反求的冲击载荷分别如图8和图9所示，可以看出采用了正则化技术后，识别载荷已不再发散，但由于Tikhonov正则化的存在，反求结果已经偏离了真实解，特别是GCV准则下的正则化参数过大，导致了冲击载荷峰值与真实值的误差较大。

图6 激励位置以及传感器布置Fig.6 The load and sensors layout

采用本文方法来识别该冲击载荷，超参数设置为(kW=10,αW=1)、(U0=0,σU0=1)、(kf=10,αf=0.1)，得到载荷识别值如图10所示，可以看出载荷识别值既稳定又准确，识别的曲线能够非常接近真实载荷曲线，特别是冲击峰值比传统正则化技术的识别精度要高的多。这是因为本文的载荷识别方法具有本征自适应正则化属性，其变化的正则化参数曲线如图11所示，可以看出该正则化参数曲线类似冲击载荷曲线，整体呈冲击状，仅在载荷冲击时段具有较大的值，在冲击峰值达到参数峰值。本算例再次说明了基于Gibbs抽样的动载荷识别方法相较于传统正则化方法具有自适应正则化的优越性，适合各种类型的载荷识别问题。

图7 冲击载荷识别值与真实值(无正则化)Fig.7 The actual and identified values of the impact load (No regularization)

图8 冲击载荷识别值与真实值(L-曲线)Fig.8 The actual and identified values of the impact load (L-curve)

图9 冲击载荷识别值与真实值(GCV准则)Fig.9 The actual and identified values of the impact load (GCV criterion)

图10 冲击载荷识别值与真实值(Gibbs抽样)以及局部放大图Fig.10 The actual and identified values of the impact load (Gibbs sampling)

图11 正则化参数对比图Fig.11 The comparison results of the regularization parameter based on these two methods

4 结 论

本文在贝叶斯框架下建立了动载荷识别的Gibbs抽样方法。该法具有本征的自适应正则化性能，正则化参数即为噪声方差与载荷方差的比值。数值算例验证了本文方法的有效性，并与Tikhonov正则化方法对比，结果表明本文方法的精度要高于Tikhonov正则化法。当前的载荷识别方法大都以准确的结构模型为基础，而本文的方法可以进一步考虑结构的模型误差，此为作者下一步的研究内容。

[ 1 ] INOUE H, HARRIGAN J J, REID S R. Review of inverse analysis for indirect measurement of impact force [J]. Applied Mechanics Reviews, 2001, 54(6): 503-524.

[ 2 ] 孙兴盛，刘杰，丁飞，等. 基于矩阵摄动的随机结构动态载荷识别技术[J]. 机械工程学报, 2014, 50(13): 148-156.

SUN Xingsheng, LIU jie, DING Fei, et al. Identification method of dynamic loads for stochastic structures based on matrix perturbation theory [J]. Journal of Mechanical Engineering, 2014, 50(13): 148-156.

[ 3 ] 杨智春，贾有. 动载荷的识别方法[J]. 力学进展, 2015: 29-54.

YANG Zhichun, JIA You. The identification of dynamic loads [J]. Advances in Mechanics, 2015: 29-54.

[ 4 ] 朱涛，肖守讷，阳光武. 一种新的时域动态载荷识别方法 [J]. 西南交通大学学报, 2012, 47(6): 968-973.

ZHU Tao, XIAO Shoune, YANG Guangwu. A new time domain method for force identification [J]. Journal of Southwest Jiaotong University, 2012, 47(6): 968-973.

[ 5 ] 徐菁，张方，姜金辉，等. 基于拟静态初值的载荷识别数值修正算法[J]. 振动与冲击, 2016, 35(2): 39-44.

XU Jing, ZHANG Fang, JIANG Jinhui, et al. Numerical correcting algorithm for load identification based on quasi-static initial value [J]. Journal of Vibration and Shock, 2016, 35(2): 39-44.

[ 6 ] LAW S S, FANG Y L. Moving force identification: optimal state estimation approach [J]. Journal of Sound and Vibration, 2001, 239(2): 233-254.

[ 7 ] MAO Y M, GUO X L, ZHAO Y. A state space force identification method based on Markov parameters precise computation and regularization technique [J]. Journal of Sound and Vibration, 2010, 329(15): 3008-3019.

[ 8 ] 马超，华宏星. 正则化技术在状态空间载荷识别中的应用 [J]. 振动. 测试与诊断, 2014, 34(6): 1154-1158.

MA Chao, HUA Hongxing. State space force identification based on regularization technique [J]. Journal of Vibration, Measurement & Diagnosis, 2014, 34(6): 1154-1158.

[ 9 ] 马超，华宏星. 基于改进正则化方法的状态空间载荷识别技术[J]. 振动与冲击, 2015, 34(11): 146-149.

MA Chao, HUA Hongxing. State space load identification technique based on an improved regularized method [J]. Journal of Vibration and Shock, 2015, 34(11): 146-149.

[10] 刘杰. 动态载荷识别的计算反求技术研究[D]. 长沙：湖南大学, 2011.

[11] ZHANG E, ANTONI J, FEISSEL P. Bayesian force reconstruction with an uncertain model [J]. Journal of Sound and Vibration, 2012, 331(4): 798-814.

[12] HASTINGS W K. Monte Carlo sampling methods using Markov chains and their applications [J]. Biometrika, 1970, 57(1): 97-109.

[13] GEMAN S, GEMAN D. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence ,1984(6): 721-741.

[14] GILKS W R. Introducing markov chain monte carlo, markov chain monte carlo in practice [M](eds. WR Gilks, S. Richardson and DJ Spiegelhalter), 1—19. Chapman & Hall, London, 1996.