基于加速度频响函数小波分解的模型修正方法

2020-07-29彭珍瑞曹明明刘满东

航空学报 2020年7期

彭珍瑞，曹明明，刘满东

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

在航空、航天、机械、土木结构和武器系统的研制过程中，有限元仿真和试验分析已成为2种不可或缺的重要手段[1]。但是在建模过程中的许多不确定性因素会导致有限元模型与试验模型存在一定的差异，因此需要利用试验模型的测试数据对有限元模型进行修正，使修正后的有限元模型能较准确地反映结构的动态特性[2]。其中，基于模态参数和基于频响函数(FRFs)的模型修正方法较为常见[3]。相比较而言，基于频响函数的模型修正方法直接利用测量得到的频响函数进行模型修正不需要模态识别，从而避免了识别误差的引入，具有一定的优越性[4]。文献[5-6]分别从理论上介绍了设计参数型频响函数修正法在工程应用中的优势及局限。殷红等[7]利用位移频响函数结合结构损伤识别，得到了与结构实际参数相符的模型修正结果。Esfandiari等[8]采用基于频响函数的模型修正方法，实现了对一标准桁架结构的损伤识别。杨海峰等[9]使用加速度频响函数对一简支梁进行了有限元模型修正。虽然基于频响函数的模型修正方法有一定的优势，但因其包含大量频率点信息，所以频率点的选择对于修正结果有重要影响，目前尚无统一方法。Kwon和Lin[10]引入测量依赖性指标和参数依赖性指标来选择最优频率点。张勇等[11]利用频响函数截断奇异值分解进行模型修正避开了频率点的选择难题。

目前，作为一种可以替代有限元模型进行快速分析的代理模型技术，引起了众多学者的关注。代理模型通过建立一个显式的函数来代替模型参数与结构响应之间复杂的隐式关系，近年来逐渐被应用到结构模型修正与确认中[12]。常用的代理模型有响应面法(RSM)、径向基函数(RBF)、神经网络(NN)、支持向量回归机(SVR)、Kriging模型等[13]。其中Kriging模型是基于Kriging插值技术的一种等效模型，与其他代理模型不同的是，Kriging模型除了能给出对未知函数的预估值，还能得到预估值的误差估计[14]，只用少量的样本就可以较准确地描述结构输入与输出间的关系，在航天飞行器设计和结构优化设计等领域应用较为广泛。文献[15-17]将Kriging方法应用于频响函数模型修正中，减少了求解时间，提高了频响函数模型修正的效率。

在上述背景下，为保留利用频响函数进行模型修正的优点，同时避免频率点的选择难题和提高模型修正效率，本文提出了一种基于加速度频响函数小波分解的模型修正方法。将计算得到的加速度频响函数进行小波分解，用小波分解得到的第1层幅值较大的小波系数来表征频响函数作为Kriging模型的输出，以模型待修正参数作为Kriging模型的输入，并通过混合灰狼算法对Kriging模型相关系数进行寻优来构建精确的Kriging代理模型。然后，以实测的加速度频响函数小波分解得到的第1层小波系数与代理模型输出小波系数的差值最小作为目标函数，通过水循环算法进行待修正参数的求解。最后，通过一空间桁架结构证明了本文方法的有效性，特别是对于随机噪声的鲁棒性。

1 加速度频响函数小波分解

1.1 加速度频响函数

频域内，一个n自由度阻尼系统在简谐激励作用下的运动方程可以表示为

(-ω2M+iωC+K)X(ω)=F(ω)

(1)

式中：M、K和C分别表示质量矩阵、刚度矩阵和阻尼矩阵；ω为激励频率；F(ω)为简谐激励；稳态响应X(ω)可表示为

X(ω)=H(ω)F(ω)

(2)

其中：H(ω)为加速度频响函数，可表示为

(3)

1.2 小波变换

小波变换是信号处理领域的一种重要研究方法，被用于图像压缩、信号滤波、数据融合和特征提取等很多方面。Donoho[18]提出了小波软、硬阈值消噪法。该方法计算量小且容易实现，得到了广泛的关注和应用。

一维含噪信号可以表示为

y(t)=s(t)+o(t)t=0,1,…,T-1

(4)

式中：y(t)为含噪信号；s(t)为有用信号；o(t)为噪声信号；t为采样点；T为采样点数，其小波变换为

(5)

其中：λjk为含噪信号的小波系数，ψ(·)为小波函数。因为小波变换是线性变换，所以式(5)可以变为

λjk=λjk(s)+λjk(o)

(6)

其中：λjk(s)和λjk(o)分别为有用信号和噪声信号的小波系数。

小波变换特别是正交小波变换具有很强的去数据相关性，它能够在小波域中使信号的能量集中在一些大的小波系数中；而噪声的能量却分布于整个小波域内。因此，经小波分解后，有用信号的小波系数幅值要大于噪声信号的系数幅值。可以认为，幅值比较大的小波系数一般以有用信号为主，而幅值比较小的系数在很大程度上是噪声信号[19]。对比多次试验效果，设定小波基为db4，当无噪声或者弱噪声(信噪比为30 dB、20 dB)干扰时分解层数设为5，当强噪声(信噪比为10 dB、5 dB)干扰时分解层数设为4，用加速度频响函数经小波分解后的第1层幅值较大的小波系数作为加速度频响函数的特征量进行模型修正。

2 基于Kriging模型的模型修正理论

2.1 Kriging模型的建立

Kriging模型是一个基于随机过程的代理模型，最初在地质统计学中提出。模型包含了线性回归部分和非参数部分：

(7)

Cov [z(xi)，z(xj)]=σ2R(xi,xj)

(8)

其中：σ2为z(x)的方差；xi和xj为样本中任意2个样本点；R(xi,xj)是相关函数，表示样本点间的空间相关性。选计算效率较好的高斯函数为相关函数：

(9)

(10)

(11)

其中：n为试验点数；σ2和β均为θk的函数，那么Kriging模型中唯一未知数即为相关系数θk，其决定着预测响应精度，本文采用混合灰狼算法对相关系数θk进行寻优。

(12)

式中：rT(x0)=[R(x0,x1)R(x0,x2)…R(x0,xn)]。

2.2 混合灰狼算法

混合灰狼算法(DE-GWO)是由Zhu等[20]提出，该算法将差分进化算法(DE)强大的搜索能力融入灰狼算法(GWO)中，使GWO能够跳出停滞区，降低了GWO寻得局部最优值的可能性，加快了收敛速度，优于粒子群(PSO)算法、差分进化算法(DE)和灰狼算法(GWO)。算法流程如下：

1) 参数初始化。设置种群规模Npop、最大迭代次数MaxIt、自变量维数nVar、缩放因子下上界beta_min和beta_max，交叉概率pCR，参数下上界lb和ub。

2) 对种群个体进行DE变异操作，产生中间体；并进行竞争选择操作形成父代种群、子代种群和变异种群个体。

3) 社会等级。计算种群中每个灰狼个体的适应度值，并依据适应度值的大小进行排序，确定灰狼种群中的社会等级，最优解Xα为头狼，第2和第3最优解Xβ和Xδ为狼群中第2和第3等级狼，其余的候选解X为狼群中最低等级狼。

4) 搜索包围猎物。灰狼搜索猎物过程中，灰狼接近并包围猎物行为可以表示为

D=|C·Xp(t)-X(t)|

(13)

X(t+1)=Xp(t)-A·D

(14)

式中：D为灰狼与猎物之间的距离；Xp(t)和X(t)分别是第t次迭代猎物和灰狼的位置向量；C和A为系数向量，且C=2r1，A=2ar2-a，其中r1和r2为[0，1]之间的随机向量，a是一个动态向量，随着迭代次数的增加从2线性递减至0。

5) 追捕猎物。当灰狼搜索到猎物所在位置时，首先，头狼Xα带领狼群对猎物进行包围；然后，头狼Xα带领Xβ狼和Xδ狼对猎物进行攻击捕捉。在灰狼群体中，Xα、Xβ、Xδ狼距离猎物位置最近，因此可以通过这三者的位置来计算灰狼个体向猎物移动的位置。并根据式(15)计算出种群中最优个体Xα、Xβ和Xδ与其他灰狼个体之间的距离，再根据式(16)计算出其余灰狼的移动方向。最后由式(17)更新灰狼位置。

(15)

(16)

(17)

式中：X1(t)、X2(t)和X3(t)代表最优个体Xα、Xβ和Xδ与其他灰狼个体X指导后更新的位置；X(t+1)为子代灰狼的最终寻优位置。

6) 对种群个体进行交叉、选择操作保留优良成分并产生新子代个体，计算个体的适应度值。

7) 更新灰狼Xα、Xβ和Xδ及其他对应的位置信息。

8) 判断是否达到最大迭代次数MaxIt，如果是则停止并输出当前最优解，否则返回3)继续执行。

2.3 混合灰狼算法优化Kriging模型

Kriging模型中的相关系数θk决定着预测响应值精度，只有当所建立Kriging模型精度足够高时，代理模型对修正结果的误差影响才会最小。所以相关系数的确定对于构建代理模型至关重要。本文采用拉丁超立方抽样方法抽取一定数量的样本点，然后把抽取的样本点分成训练集和测试集。以训练集作为Kriging模型的输入变量，以加速度频响函数小波分解得到的第1层小波系数作为响应值来构建Kriging模型，以测试集Kriging模型的均方误差(MSE)均值平均值作为混合灰狼算法的目标函数，寻得最优相关系数。最后以训练集作为样本点建立Kriging模型。

2.4 模型修正

当Kriging模型构造完成时，模型修正归根结底是一优化问题，在目标函数最小的条件下，通过求解该优化问题，得到设计参数的修正值。当激励点和响应点确定时，目标函数为

(18)

在优化求解时，除了灵敏度方法外，智能算法也是很简便的。水循环算法(WCA)[21]的基本概念和思想来源于自然，并基于现实世界中的水循环过程，具有良好的随机搜索能力、鲁棒性和快速收敛等优点。故本文使用水循环算法进行优化求解待修正参数，模型修正流程如图1所示。

图1 模型修正流程图Fig.1 Flow chart of model updating

3 数值算例

选择一空间桁架结构来验证本文方法。该结构包括28个节点，66个单元和48个自由度。桁架约束条件为4个支座固定(节点编号:1、8、9、16)，节点铰接，如图2所示，桁架结构参数如表1所示。

表1 三维桁架结构参数Table 1 Structure parameters of three-dimensional truss

根据样本点参数，运用动力学方法计算其质量矩阵、刚度矩阵并且添加一定的比例阻尼，将其代入式(3)中计算出其加速度频响函数矩阵。由于本文所用方法间接利用频响函数，避开了频率点选择问题，所以在频率点的选择问题上，只需保证频响函数在感兴趣频带内有足够数量的共振峰。所选结构的激励点和响应点如图2所示。

3.1 待修正参数选取与试验设计

根据表1将各杆的弹性模量E和密度ρ、上弦杆横截面积A1、下弦杆横截面积A2、直腹杆横截面积A3、斜腹杆横截面积A4，这6个参数作为初选待修正参数。以表1中的各参数数值作为“试验模型”参数，然后对“试验模型”参数进行偏移得到“有限元模型”参数。其中，弹性模量E增加10%，密度ρ减小10%，上弦杆横截面积A1增加10%，下弦杆横截面积A2减小10%，直腹杆横截面积A3增加10%，斜腹杆横截面积A4减小10%，得到对应的有限元模型参数值如表2所示。

表2 试验模型和有限元模型参数Table 2 Parameters of test model and finite element model

对于全局优化问题，比较好的方法是通过试验设计(DoE)来选取一组样本点。对于Kriging模型，初始样本点数理论上不受设计空间维数的限制，且优化效率对初始样本点数的依赖也不明显[14]。为此，采用拉丁超立方抽样在上述有限元模型6个参数的±20%区间内抽取200个样本点，使所构建Kriging模型满足精度要求。根据样本点运用全局灵敏度中的Sobol法计算各参数对于响应的灵敏度值，将各参数灵敏度值归一化并转换成百分比形式，如图3所示。

由图3可以看出，各杆的弹性模量E、密度ρ和上弦杆横截面积A1灵敏度最高，而其他杆的横截面积灵敏度相对较小，所以选择E、ρ、A1作为待修正参数进行模型修正。

图3 各参数灵敏度Fig.3 Sensitivity of each parameter

3.2 Kriging模型的构造及评估

基于拉丁超立方抽样抽出的200个样本构造Kriging模型。用前150个样本数据作为训练集，后50个数据作为测试集。以训练集作为Kriging模型的输入，以训练集所对应的加速度频响函数小波分解得到的第1层的小波系数作为输出，构造初始Kriging模型。以测试集Kriging模型的均方误差(MSE)均值平均值最小作为混合灰狼算法的个体目标函数，寻得最优相关系数θk。混合灰狼算法参数设置：种群规模Npop=40，最大迭代次数MaxIt=300，自变量维数nVar=1，缩放因子上下界beta_max和beta_min分别为0.2和0.8，交叉概率pCR=0.2，参数上下界ub和lb分别为0.01和100。最终寻得最优相关系数为0.875 8，迭代曲线如图4所示。

图4 混合灰狼算法迭代曲线Fig.4 Iteration curve of DE-GWO

由图4可以看出，当迭代次数增加至250次时，迭代曲线基本收敛，寻得最优值。相关系数确定以后，就可以构造Kriging模型。为评估所构建Kriging模型的精度，训练集样本第5个小波系数真实值和Kriging模型预测值的拟合曲线如图5所示。为进一步评估Kriging模型的预测精度，测试集样本第10个小波系数真实值和Kriging模型预测值的拟合曲线如图6所示，并根据式(19) 计算出小波系数真实值和Kriging模型预测值间的均方根误差(RMSE)，当RMSE越接近0时Kriging模型预测精度越高。

(19)

由图5看出，通过训练集构建的Kriging模型预测值和真实值几乎重合，表明训练集的拟合精度很高。训练集的拟合精度并不能说明Kriging模型预测精度，必须用训练集以外的样本点来验证Kriging模型的预测精度，本文用测试集样本来验证。由图6看出，除个别样本点外，测试集样本点对应的Kriging模型预测值几乎和真实值重合，并且通过式(19)计算出的RMSE值为8.164 1×10-4，表明所构建的Kriging模型预测精度很高，可以代替理论模型进行模型修正。

图5 训练集样本第5个小波系数预测值Fig.5 Predicted value of the fifth wavelet coefficient of training set sample

图6 测试集样本第10个小波系数预测值Fig.6 Predicted value of the 10th wavelet coefficient of test set sample

3.3 模型修正

当完成Kriging模型构造以后，就可以根据式(12)通过水循环算法以式(18)作为目标函数迭代求解待修正参数。假定试验参数值在有限元参数值的±20%区间内，使用水循环算法进行迭代寻优。水循环参数设置：雨滴层数为50，河流和大海个数总和为4，极小值为1×10-5，最大迭代次数为100。修正结果如表3所示，水循环迭代曲线如图7所示。

由图7可以看出水循环算法具有良好的搜索能力和快速收敛性能，算法在10次迭代以前就已收敛。多次运行寻优程序，迭代曲线基本不变化，也表明了该算法寻优的鲁棒性。由表3可以看出，用本文所提方法进行模型修正取得了很好的修正效果。使用表3中的各参数试验值和有限元值计算出试验模型和有限元模型对应的加速度频响函数。使用所得修正参数的平均修正值计算得到修正后的加速度频响函数，称其为Kriging模型加速度频响函数。Kriging模型、试验模型及有限元模型加速度频响函数曲线如图8所示。

表3 模型修正参数及其误差Table 3 Parameters and their errors of model updating

图7 水循环迭代曲线Fig.7 Iteration curve of WCA

图8 Kriging模型、试验模型及有限元模型频响函数曲线Fig.8 FRF curves of Kriging model, test model and finite element model

由图8可以看出修正过后的结构加速度频响函数曲线与试验加速度频响函数曲线几乎重合，与有限元模型加速度频响函数曲线形状相似，只是沿坐标轴有一定的移动。

为进一步检验本文修正方法的修正效果，比较频响函数修正前后实、虚部曲线，如图9所示。由图9可以看出修正后的有限元模型无论是实部还是虚部，都能很好地与试验模型的频响函数曲线重合。这进一步证明了本文所提方法的有效性。

图9 修正前后频响函数曲线Fig.9 FRF curves before and after updating

实测试验时目标频响函数不可避免会受到随机噪声的干扰，在一些特殊环境下外界噪声干扰还很大。所以在计算的试验频响函数中加入信噪比为30、20、10和5 dB的高斯白噪声，来验证本文方法对于噪声的鲁棒性。图10为加入信噪比为20 dB和5 dB时加速度频响函数图。

图10 20 dB和5 dB信噪比下加速度频响函数Fig.10 FRF under 20 dB and 5 dB SNRs

将加入噪声的加速度频响函数按前文所提小波参数设置方法进行小波分解，得到不同信噪比下小波系数曲线和局部放大曲线，如图11所示。由图11可以看出，将加速度频响函数进行小波分解以后得到的序号较小的小波系数幅值较大，而绝大部分小波系数趋近于0。当在信号中加入随机噪声时，随着加入的高斯白噪声信噪比逐渐减小，对于幅值较大的小波系数影响较小，而对于幅值小的小波系数干扰较大。为了更直观地看出不同信噪比下噪声对于小波系数的影响，将小波系数随小波系数序号的变化分成3个部分，如各图中的3个放大曲线所示。放大曲线中左侧为加速度频响函数小波分解得到的第1层小波系数，个数较少，幅值较大；右上为加速度频响函数小波分解得到的中间靠前的小波系数，个数中等，幅值有大有小；右下为加速度频响函数小波分解得到的中间靠后的小波系数，个数较多，幅值趋近于0。

根据图11(a)～图11(d)，由其中的左侧放大曲线可以看出第1层小波系数曲线随着加入的噪声加大，加噪小波系数曲线与未加噪小波系数曲线相似，大部分重合，只是当噪声很大时略有波动。由右上放大曲线可以看出，随着噪声的加大幅值大的小波系数曲线波动较小，而幅值小的小波系数曲线波动开始加大。由右下放大曲线可以看出，随着噪声的加大对于幅值趋近于0的小波系数影响较大，噪声越大小波系数曲线波动越大。综上所述，加速度频响函数经小波分解后得到较少数目幅值较大的小波系数，绝大部分小波系数趋近于0。当加速度频响函数受到噪声干扰时，对于幅值较大的小波系数干扰较小，而对于幅值较小的小波系数干扰较大。

图11 不同信噪比下小波系数曲线和局部放大图Fig.11 Wavelet coefficient curves and local enlarged curves under different SNRs

在不同噪声水平下连续运行10次程序的模型修正结果均值如表4所示。由表可以看出/有无噪声对于模型修正结果影响很大，而当在有噪声情况下，随着噪声的加大用本文所述方法仍然能够得出较满意的模型修正误差，即使当信噪比低至5 dB时，各参数模型修正误差仍然低于4%。由此证明了本文所提模型修正方法对于噪声具有较强的鲁棒性。