基于频响函数提升小波总能量的桁架结构随机模型修正

2021-07-01赵永鹏彭珍瑞

计算力学学报 2021年3期

赵永鹏，殷红，彭珍瑞

(兰州交通大学机电工程学院，兰州 730070)

1 引言

随着工程技术领域的不断发展，有限元技术的应用更加广泛，但是各种不确定性因素，如建模时的简化及安装误差等对结果的影响很大，因此建立一个可靠有效的有限元模型就变得尤为重要。当前，通过模型修正来获得准确的有限元模型(FEM)是主流的方法。模型修正方法分为确定性和不确定性两种[1]，确定性方法在近年来有了广泛发展。但是，由于各种不确定性的情况，出现的误差较多。所以，在模型修正中充分考虑不确定性的影响，从而构建一种随机性的模型修正方法势在必行。

参数不确定性的模型修正方法指通过多次模态实验得到相应数据，计算出具有统计特征的有限元模型，来表明实际结构的静动力学行为、机理和演化规律等。方圣恩等[2]充分考虑了修正参数的不确定性，通过计算参数的统计特征量对一组钢板进行了修正。Fang等[3]通过Hermite多项式的显式多项式混沌展开，把不确定的参数和响应通过构造一个随机响应面模型(RSM)对一组金属板进行了修正。Zhai等[4]以模态参数为响应，基于改进的响应面模型和先进的蒙特卡洛(MC)方法，对一航空发动机定子系统(壳体)进行了修正。Bi等[5]构建了一个完全嵌入Bhattacharyya距离的近似贝叶斯模型的修正框架，对一质量弹簧系统进行了模型修正，证明了方法的有效性。Deng等[6]通过构建代理模型、设置校准参数以及将距离函数作为参数相关性的衡量标准等步骤进行了随机模型修正。Khodaparast等[7]利用模态响应(如固有频率和模态振型)中得到的可变性来估计随机修正参数的第一和第二统计矩，通过一个三自由度的质量弹簧系统对所提方法进行了验证。

通过结构模态参数进行修正将不可避免地在识别的过程中引入一些误差，实际工程中由于实验模态参数的识别而出现的误差有可能会在某些特定情况下比建模精度低引起的误差更大。而基于频响函数(FRFs)的模型修正方法因为可以略去结构模态参数识别而得到广泛应用，该方法减少了各种误差的出现，适于对模态分布相对比较密的结构进行修正。李伟明等[8]利用频响函数为输出对一个二维桁架进行了修正。王巨涛等[9]利用加速度频响函数曲线对应频率点处的响应值之差构造目标函数,对GARTEUR飞机模型进行了修正。但是频响函数包含的数据信息量很大，怎样选出好的频率点数据以得到期望的修正效果，现在还没有形成统一的标准。因此，可以考虑引入一个特征向量对频响函数进行表征，以避免上述问题。张勇等[10]提取频响函数的特征向量对一扭力梁模型进行了修正。上述研究大都应用于确定性的模型修正，但鲜见应用于不确定性模型修正。提升小波变换LWT(Lifting Wavelet Transform)属于第二代小波变换，与传统小波变换相比，提升小波变换不仅拥有多分辨率的特点，而且具有计算时间少，整个流程相对简单以及变换后数据长度不变等优势[11]。宋萌萌等[12]选择提升小波变换对得到的初始信号进行降低噪声的处理，得到的结果满足要求。袁海英等[13]利用提升小波变换对振动信号进行时频特性分析和信息预处理，通过预测器和更新器的设计取代了小波基函数选取过程。提升小波变换不再使用傅立叶变换，也没有利用函数的平移或伸缩，而是在时域中就可以完成原始信号不同频带上的分离。经过变换后各层分量的小波能量对信号的变化有很高的敏感性，且对噪声污染、小波种类选择及小波层数具有鲁棒性。因此，选择提升小波变换对频响函数进行分解，提取总能量值来表征频响函数可以最大程度上保证频响函数的特征不发生变化，而且也减少了频响函数信息过多引起的各种处理困难。

本文构建了一种以提升小波总能量为响应，通过筛选响应样本反求出待修正参数的随机模型修正方法。首先对加速度频响函数进行提升小波变换，提取对应的提升小波总能量作为输出，以待修正参数作为输入，构造一个多项式响应面代理模型代替有限元模型；然后利用蒙特卡洛抽样得到大量的响应样本，并通过设置阈值对响应样本进行优化筛选；最后构造一个优化反问题，利用布谷鸟优化算法反求出输入即待修正参数，计算修正后参数的统计特征。分别对平面桁架和三维钢桁架结构模型进行修正，证明了所构建方法的有效性。

2 频响函数提升小波总能量

2.1 加速度频响函数

在具有阻尼的系统中，通过一个简谐激励作用所得的动力学方程为

(1)

在频域内的输入输出关系可以表示为

A(ω)=H(ω)F(ω)

(2)

式中A(ω)为稳态响应，H(ω)为加速度频响函数，F(ω)为简谐激励，ω为激励频率。

相应的加速度频响函数可以表示为

(3)

2.2 提升小波变换

对于得到的频响函数，选择提升小波变换对其能量值进行提取。因为常见的小波滤波器的输出是浮点数，所以对变换后的值选择压缩时需要量化才能得出对应的整数。Swelden[14]选择了根据提升方式改进小波特性得到小波基函数的新模式即提升小波变换。提升小波变换共分为三个过程: 分裂、预测和更新。

(1) 分裂。根据原始信号sj奇偶性的不同可以得到两个没有相关性的子集。所得子集的长度是初始信号的1/2。如式(4)所示：

sj=(ej - 1,oj - 1)

(4)

式中ej - 1={ej - 1,k=sj,2k}为偶序列，oj - 1={oj - 1,k=sj,2k + 1}为奇序列。

(2) 预测。根据第一步提取的偶序列ej - 1来预测奇序列oj - 1。如式(5)所示，用实际值oj - 1与预测值ej - 1的差值dj - 1来表示其接近程度，称作细节系数，同时也属于原始信号的高频分量。

dj - 1=oj - 1-P(ej - 1)

(5)

式中P为预测函数，用相邻数据的平均值来表示：

pk(ej - 1)=(ej - 1,k+ej - 1,k + 1)/2

(6)

(3) 更新。这一步的目的是要使第二步产生的某些特征(如均值等)与原始信号的特征保持不变。可以通过算子U来实现：

sj - 1=ej - 1+U(dj - 1)

(7)

式中sj - 1为原始信号的低频部分，U为更新算子函数，可以表示为Uk(dj - 1)=dj - 1,k/2。

通过提升小波变换之后，可以将原始信号sj变为高频分量dj - 1和低频分量sj - 1；将低频部分sj - 1继续进行相同的步骤，又可以分解为二层高频dj - 2和二层低频sj - 2；以此类推，可以根据要求进行n层的提升小波分解，原始信号就可以表示为{sj - n,dj - n,…,dj - 1}。其中sj - n为信号的低频部分，{dj - n,…,dj - 1}为信号的高频部分。

对频响函数进行提升小波变换后，即可得出各层的高频部分和低频部分。设Ej d为原始信号sj在尺度n上的高频分量总能量：

(8)

相应的低频分量能量Ej s也可以得到，信号的提升小波总能量为各个尺度下的小波能量的总和[15]。

Ez=E1d+E2d+…+Ej d+Ej s

(9)

3 响应面模型

对于一些大型复杂结构，由于其内部存在很多的零部件，以及很多边界尺寸数据的不确定性，所以在有限元建模的过程中存在困难，从而建立的有限元模型变得十分复杂，在模型修正过程中可能会出现计算无法收敛或者所用时间很长等问题。选择一个代理模型来代替这个复杂的有限元模型，可以方便计算，提高效率。多项式响应面代理模型具有建立相对简单且预测精度高的优点，因此本文选择建立响应面模型来代替有限元模型。这个代理模型将一个物理系统的输入参数x和输出响应y通过数学多项式结合起来：

y=f(x1,x2,…，xn)+ε

(10)

式中f为输入与输出间的映射，ε为建模误差，n为输入参数的数量。

对于一般性的工程问题，大多采用二阶多项式模型作为一种基本的形式，即

(11)

式中bo,bi,bi i和bi j为回归系数，xi和xj为输入参数。

响应面构建完成后，需对其进行精度检验，判断是否满足后续计算要求，常用的检验方法为计算真实值与预测值间的均方根误差，即

(12)

因此,根据提取的输出与输入参数就可以构建二阶响应面代理模型，然后再对其进行精度检验。

4 响应样本获取

4.1 蒙特卡洛抽样

蒙特卡洛方法提供了一种可替代分析数学评估统计量在随机样本中行为的方法。根据样本的概率分布，用计算机实现统计模拟或抽样，属于一种随机抽样方法。常见的数值计算方法由于需要多次迭代，导致最后误差累积很大，影响结果，而蒙特卡洛抽样虽然也有迭代，但是没有像其他方法那样频繁迭代，产生的误差一般可以忽略不计，对最后的抽样结果也是基本没有影响。因此，在已知响应概率分布的情况下，运用蒙特卡洛抽样可以得到大量的响应样本，认为这些样本就是试验值，然后用于随机模型修正过程。

4.2 响应样本筛选

通过蒙特卡洛抽样得到大量响应样本，由于是一种随机的抽样方法，所以可能会存在抽取的样本波动范围过大，以至于超出了实际试验所能测得的范围，样本称之为失真样本。如果将这些数据也代入修正过程，就失去了实际的意义，可能导致修正的结果很不理想。因此，可以通过设定一个阈值来对抽样数据进行一次优化筛选，使得筛选后的数据与实际情况更加接近。本文通过计算抽样响应值与真实响应值间的均方根误差进行样本的筛选，

(13)

为了保证所抽取样本的真实性并且选择有效的样本，按照至少保留80%左右样本的准则进行筛选。计算所有的RMSE值，根据该准则设定一个阈值，当RMSE值比阈值小时保留，反之则淘汰。

5 模型修正

根据上述理论，通过响应面模型代替有限元模型后，以蒙特卡洛抽样后优选的响应值与响应面模型预测响应值差值最小来构造目标函数。

(14)

建立目标函数之后，模型修正变成对目标函数求解的问题。由于布谷鸟优化算法具有全局搜索能力强和参数少等优势，因此本文选择该算法进行寻优以得到修正后的参数值。

综合频响函数提升小波总能量的提取，响应面模型的建立，响应样本的筛选，目标函数的求解等步骤，可得随机模型修正流程如图1所示。

图1 模型修正流程

6 数值算例

6.1 平面桁架

如图2所示平面桁架，该桁架结构共有14个节点和25个自由度。杆单元的弹性模量为210 GPa，材料密度为7850 kg/m3。将该平面桁架结构作为试验模型，选择弹性模量和密度为待修正参数。对参数进行偏离，其中密度减小10%，弹性模量增加10%作为相应的有限元模型列入表1。选择第11节点y方向为激励点位置，第3节点y方向为测点位置。

图2 平面桁架模型

表1 试验模型与有限元模型参数

选择拉丁超立方法获取样本点。在待修正参数的±20%区间内抽取40个样本点。对于抽样所得40组样本，选择前20组为样本集拟合响应面模型，后20组为测试集，用来验证代理模型的精度。通过回归拟合出响应面模型如图3所示。

图3 平面桁架的响应面模型

假设试验值均服从正态分布。依据抽样所得样本点，可以得到其加速度频响函数值，对其进行三层的提升小波变换，提取相应的提升小波总能量。根据蒙特卡洛方法完成抽样，共抽取200个样本。对于抽取的样本通过设定阈值进行筛选，按照至少保留80%数据筛选过后，代入修正过程进行求解。

对构建的响应面模型进行精度检验后，所得均方根误差RMSE值如图4所示。可以看出，RMSE值最大不超过0.06。由此可得，所构建的响应面模型精度满足计算要求，该模型有效。

图4 响应面模型RMSE值

通过上述方法进行模型修正，修正结果列入表2。由表2可知，修正误差小于2.1%，达到了较高的修正精度。

为了进一步说明方法的有效性，修正前后的频响函数曲线如图5所示。可以看出，修正后的频响函数曲线与真实值的曲线重合度高，证明用提升小波总能量来表征频响函数作为结构响应进行模型修正是可行的。

图5 修正前后频响函数曲线

6.2 三维桁架

如图6所示三维桁架。该桁架共有28个节点，66个单元，48个自由度。总长2.8 m，宽0.39 m，高0.27 m。桁架共四个支座固定(节点编号为1,8,9,16)，节点铰接。桁架的材料参数弹性模量为190 GPa，密度为7800 kg/m3。尺寸参数上弦杆和下弦杆横截面积为85.5 mm2，直腹杆横截面积为141.0 mm2，斜腹杆横截面积为45.0 mm2。结构的激励位置和响应位置分别如图6所示。

图6 三维桁架模型

6.2.1 试验设计与样本选取

选择如图6所示三维桁架结构作为本算例的试验模型，相应的有限元模型通过对待修正参数进行一定范围的偏离得到。选择材料参数的弹性模量和密度以及尺寸参数的上弦杆横截面积共三个参数为待修正参数。对参数进行偏离，其中密度减小10%，弹性模量增加10%，上弦杆横截面积减少10%，有限元模型参数列入表3。

表3 试验模型与有限元模型参数

建立响应面模型，样本点的选择很重要。常用的试验设计方法有中心复合法、全因子法和拉丁超立方法等。本文选择拉丁超立方法获取样本点，该方法采样比较均匀且样本具有代表性，符合模型建立的要求。在三个待修正参数的±20%区间内利用拉丁超立方抽样抽取40个样本点。

假设试验值均服从正态分布。根据蒙特卡洛方法完成抽样，共抽取2000个样本。对于抽取的样本通过设定阈值进行筛选，按照至少保留80%数据(即1600个)的要求筛选后，代入修正过程进行求解。

6.2.2 提升小波总能量的提取

依据拉丁超立方抽样所得样本点，可以得到其加速度频响函数值，如图7所示。以试验值的频响函数曲线为例，对其进行三层的提升小波变换，可以求得各层高频分量和低频分量如图8所示，然后提取相应的提升小波能量，将提取的各层能量求和，就可以得到提升小波总能量。

图7 试验模型频响函数曲线

图8 三层提升小波变换分量

6.2.3 响应面代理模型的构造及检验

对于抽样所得40组样本，选择前20组为样本集作为响应面模型的输入x，后20组为测试集，用来验证代理模型的精度。将20组输入与其所对应的输出(即频响函数提升小波总能量)通过回归拟合求出响应面模型的各个参数，得到代理模型表达式为

y=18.2007-0.0011x1-0.0091x2-215.08x3+

0.0013x1x3+0.043x2x3

(15)

将后20组样本代入计算得到相应的输出，然后对构建的响应面模型进行精度检验。所得均方根误差RMSE值如图9所示。均方根误差越接近0表明响应面精度越高。可以看出，RMSE值不超过0.07，所构建的响应面模型精度满足计算要求，该模型有效。

图9 响应面模型RMSE值

6.2.4 模型修正

通过蒙特卡洛抽样得到2000个响应样本如图10 所示，通过多次预实验发现设定阈值为0.092时满足筛选准则要求。对所有样本排序后进行筛选，经过筛选后共有1643个样本符合要求。如图11所示为筛选后的样本。将这些样本依次代入目标函数中，通过布谷鸟算法对目标函数求得最小值，可得到对应的1643组输入即为修正后的参数。算出修正后各参数的均值，修正前后的参数均值列入表4。

图10 蒙特卡洛抽取的响应值

图11 优化筛选的响应值

表4 修正前后参数均值

由表4可知，修正后参数均值的误差相比于修正前有了明显的减小，说明修正效果良好。修正前后的提升小波总能量值列入表5。由表5可知，提升小波总能量的误差比修正前也有了大幅的减小，表明该方法完成随机模型修正是可行的。采用修正后的参数均值可得其加速度频响函数如图12 所示。图中修正后的频响函数曲线与真实值的曲线基本重合，相较于修正前的曲线，修正效果达到了预期的目的。

表5 修正前后提升小波总能量

图12 修正前后频响函数曲线

图13和图14分别给出修正前后频响函数实部和虚部曲线，可以看出，修正值与真实值的实部虚部曲线基本重合，进一步说明修正结果良好。

图13 频响函数实部曲线

图14 频响函数虚部曲线

为了进一步验证所提方法的有效性，以平面桁架为例，分别将频响函数和提升小波总能量作为响应进行模型修正,结果列入表6，所有结果均是在同一台计算机上运行20组。可以看出，将提升小波总能量作为响应进行模型修正所用时间更短；在修正精度方面，弹性模量修正精度有了很大提高，虽然材料密度修正精度略有降低，但是总体来说效率有了一定提高，计算量也有所减小，能够达到随机模型修正的目的。