基于模态参数灵敏度的模型修正方法研究

张 淼

(长春工程学院理学院，长春 130012)

结构分析的基本目标之一是公式化能由实测数据检验的结构有限元分析模型。这个模型常常不能一开始就产生可靠的基频和模态振型，为此引入实测数据来调整此分析模型，直到分析与实验数据得到匹配。修正后的有限元模型能够更好地反映结构的动力特性。为此，需要从不完全的实测模态参数数据中建立刚度和质量，甚至阻尼矩阵，而从实验中测量出来的数据在数量上必须超过在结构中识别出来的参数的数量，因此尽可能多地应用实测数据是十分重要的。如果只用频率数据，可识别参数的数量就十分有限，那么测量的模态数据也应该包含在修正过程之中。数据和建模精度都是修正过程中需考虑的重要因素，当然模态测量的精度不如频率的测量精度，因此任何一个修正程序如果使用最小二乘法，就必须包含数据精度的权重。

传统的模型修正方法大致分为两种，即矩阵元素型和设计参数型，灵敏度分析在有限元模型修正过程中常常被使用，有时甚至成为决定模型修正成败的关键因素。近年来由于灵敏度分析以及非线性优化技术的不断进步，使得设计参数型模型修正方法的发展成效更为显著一些。基于灵敏度的模型修正方法，是目前模型修正领域中发展较为活跃的技术之一。它的应用经常融合了残差型和灵敏度矩阵方程型的目标函数，本文的目的是试图从数学角度揭示这两种模型修正技术之间的内在联系，并充分展示它们的技术轮廓。

1 基本理论[1]

描述自由度为N的线性阻尼离散系统的自由振动方程为

(1)

(2)

令λi=-ωi2，则式(2)可化为

(K+λiM)vi=0。

可以计算得到系统的实频率λ1，λ2，…，λN及实模态v1，v2，…，vN。记V=[v1，…，vN]为无阻尼规范化振型矩阵，那么此时模态质量和模态刚度矩阵分别为

VTMV=E，VTKV=diag(-λ1，…，-λN)。

对于经典阻尼系统，模态阻尼矩阵为

VTCV=diag(c1，…，cN)。

考虑阻尼时的系统极点及复模态对(si，ui)(i=1，2，…，2N)满足方程

对于N自由度振动系统，特征方程det[s2M+sC+K]=0有2N个呈复共轭对出现的特征值s1，s2，…，s2N(其中si+1为si的共轭(i=1，3，…，2N-1))，称为系统的极点或复频率。这些频率对应着一组呈复共轭对出现特征向量ui∈CN称为与si相对应的第i个模态向量。将u1，u2，…，u2N(其中ui+1为ui的共轭(i=1，3，…，2N-1))，称为复模态。复模态和复频率统称为复模态参数。复模态的正交规范化的形式有很多，对称系统常用的形式为

ΦTAΦ=E，ΦTBΦ=diag(-λ1，…，-λ2N)。

其中状态向量矩阵为Φ=[φ1，φ2，…，φ2N]，状态向量为φi=[uiλiui]T(i=1，2，…，2N)，且

2 基于残差型目标函数的模型修正方法

文献[2]中，为了修正有限元模型参数，基于某些修正参数，把测量得到的模态数据作为参考数据，通过解一个非线性最小二乘问题来优化一个价值函数，它是由实测和计算得到的解析频率之差的平方和构成的。特别指出有限元模型修正是一个非线性参数估计问题，因为模态数据是未知结构性质的非线性函数：

其中

例如若选择板的长度作为修正参数，为了解非线性最小二乘问题，需要确定板长度的上下限，然后作为初始值，通过增加回归阶数的回归算法用来解插值多项式的回归系数，即通过已知的模态频率来拟合一条插值曲线，例如线性回归拟合。基于初始计算参数值，用修正参数L的已知值和计算出来的有限元频率进行一个高阶拟合，优化价值函数，有限元频率用插值多项式替换，当参数修正后的频率残差变得比使用者定义的准则还小的时候，令修正算法终止，从而价值函数有如式(3)的形式：

(3)

其中

aij是插值多项式的系数。

最小二乘问题允许考虑残差的权重，依据这些残差的重要性和噪音的数量，加权重来区分不同的数据系列能使方法变得更实用，但同时要求工程师能提供正确的权重，例如如式(4)模型

(4)

(5)

其中wj是残差rj的权因子。权因子可以依据测量的随机性质定义，也可以把相关测量的方差的倒数作为权因子。

文献[2]还提出了另一种价值函数的定义方法

其中

文献[3]中提出了有限元模型修正的多目标优化模型，一般的公式如式(6)：

(6)

式中：θ是修正参数向量；Fi(θ)是目标函数，它是一个数量值函数；n是子目标函数的个数；wi是相应的子目标函数的权因子；gj(θ)与hk(θ)是等式和不等式约束；θL与θU是修正参数的上下限，权因子在平衡子目标函数的重要性方面具有重要的作用。它的选择依赖于使用者的经验和专业技能。在此文献中用遗传算法求解多目标规划问题。

文献[4]尝试了几种不同的权和的表达形式，并建议最平衡的一种可以选择为模态和频率如式(7)形式的加权残差：

(7)

文献[5]在基于灵敏度分析的模型修正程序中给出的目标函数是一种联合了模态和频率的优化问题：

(8)

文献[6]中目标函数f(θ)：Rn→R是一个普通的加权最小二乘问题：

(9)

rj(θ)：Rn→R(j=1，…，m)是残差向量，它是关于修正变量θ∈Rn的非线性函数；r(θ)=[r1(θ)，r2(θ)，…，rm(θ)]T；W=diag(w1，…，wj，…，wm)，W∈Rm×m。残差向量的第一部分包含无阻尼固有频率的解析值及实测值之差：

(10)

(11)

文献[7]用两个模态参数，即模态频率和模态振型来组成残差向量，构建目标函数，一般来说它是非线性最小二乘问题：

(12)

r：Rn→Rm，θ∈Rn

式中：θ代表未知的修正参数组成的向量；r代表由频率残差rf和模态振型残差rm组成的残差向量：

(13)

(14)

式中W是权矩阵。

(15)

式中rj(θ)代表测量与解析模态参数之差，它是一个关于优化变量θ∈Rn的非线性函数。为了获得唯一解，m(残差的个数)要多于未知数θ的元个数n，残差向量r(θ)：Rn→Rm，分成频率残差rf，模态振型的残差rs和一个模态质量的残差向量rm：

(16)

(17)

(18)

(19)

文献[9]中指出基于灵敏度的有限元模型修正方法是最频繁用于模型修正的方法。模型修正的优化问题为：

(20)

这种方法直接比较频率与模态振型：

rf：Rn→Rmf，rs：Rn→Rms，

(21)

所以加权特征值残差可以表示为

在模型修正中，模态振型都是质量规范化后的，当环境振动用于提取模态参数，实验模态保持未缩放的，因此模态振型残差定义为，

这个因子是可变的，给出频率与模态振型的重要性。实验频率一般是最精确的实验数据，而实验模态振型具有更多的噪音，会减少优化问题的稳定性，但能提供更多局部的信息，因此一般适当加权是很必要的。

在文献[5]～[9]中都采用的是解非线性最小二乘优化问题的信赖域牛顿高斯算法，其优化模型分别为本文中的公式(8)～(21)，在每步迭代中搜索步长限定在一个信赖域内，避免了不希望出现的大步长。f(θ)可以用一个二次型来近似，它来自于泰勒级数截断：

minm(z)=fs+{2fs]z(‖z‖≤Δu)，

式中：z是步长向量(Δθ)；{fs}和[2fs]是f(θ)的梯度及海森矩阵。{fs}={ri(θ)}=J(θ)Tr(θ)，

式中J(θ)是由ri(θ)的一阶偏导数构成的灵敏度矩阵。信赖域牛顿方法是一种迭代方法，是将高斯牛顿方法用于信赖域中，来提高收敛性，在信赖域方法中，算法构成了一个价值函数mk，它的行为类似于目标函数f，在当前点θk附近，并围绕θk得到一个区域，这个模型是可信的，尽管海森阵可以用J(θ)TJ(θ)近似。标准的高斯牛顿方法是相当稳定的，而且通过在信赖域中实施，可以加快收敛。

3 基于灵敏度矩阵方程的模型修正方法

将g个修正对象f1，f2，…，fg(可以包括实或复的模态参数，质量、刚度矩阵，反共振频率，模态的MAC值，频响函数等，也可以是这些对象的组合)均看作是m个设计参数p=(p1，p2，…，pm)T的函数f(p)，在设计参数(可以取为材料常数，几何参数或边界条件等)的初始取值位置作一阶泰勒展开，使问题线性化：

(22)

(23)

并认为灵敏度矩阵

(24)

为已知或者可求，则式(23)还可表示为

SΔp=ΔR，

(25)

式中Δp=(Δp1，Δp2，…，Δpm)T为待解修正向量；ΔR=f(p(u))-f(p(o))为残差向量。求解方程(25)，即可解出设计参数的修正向量Δp。

上文中指出，当修正对象不同时，灵敏度矩阵S的计算技术也必然不同。例如当f为实频率时，其泰勒展开式为

式中λ=(λ1，λ2，…，λN)T，则灵敏度矩阵方程为

SΔp=ΔR，

其中

(26)

ΔR=λ(p(u))-λ(p(o))，

(27)

其中λ(p(u))为实频率的实验测量值，而λ(p(o))为实频率的解析计算值。

当f为复频率时，只需将式(26)～(27)中的λ=(λ1，λ2，…，λN)T换成s=(s1，s2，…，s2N)T；当f为实模态时，只需将式(26)～(27)中的λ=(λ1，λ2，…，λN)T换成[v]=(v1，v2，…，vN)T；当f为复模态时，只需将式(26)～(27)中的λ=(λ1，λ2，…，λN)T换成[u]=(u1，u2，…，u2N)T，其中[u]和[v]表示由相应的复模态和实模态向量构成的矩阵。如果在实测数据不完备的情况下，还可以将N或2N换成实测数据的个数t，此时方程(25)是超定的。

在许多实际问题中所遇到的方程组Ax=b的系数矩阵往往是奇异的或是m×n任意矩阵(一般m≠n)，显然不存在通常的逆矩阵A-1来求解方程组，这就促使人们推广逆矩阵的概念，引进某种具有普通逆矩阵类似的性质的矩阵G，使得方程组的解仍可表示为x=Gb。而广义逆理论能够把相容线性方程组的一般解、极小范数解以及矛盾方程最小二乘解、极小范数最小二乘解(最佳逼近解)全部统一起来。

1)当方程组相容时，若系数矩阵A∈Cn×n，且非奇异，则有唯一解x=A-1b。

若灵敏度矩阵方程(25)中的系数矩阵的秩n与设计参数的个数m相同，则可以直接解得：

Δp=S-1ΔR，

在很多情况下灵敏度矩阵方程中的系数矩阵的秩n与设计参数的个数m并不相同，当n>m时，引入广义逆来求解得[11]：

Δp=(STS)-1STΔR；

当n

Δp=ST(SST)-1ΔR。

从灵敏度矩阵方程(25)出发，看成是最小化残差的优化模型，令目标函数为[12]

minJ(Δp)=(SΔp-ΔR)TWE(SΔp-ΔR)，

Δp=(STWES)-1SWEΔR。

(28)

然而通常由于灵敏度矩阵S会接近奇异，可以同时引入了两种加权矩阵Wp和WE，来克服因灵敏度矩阵S病态而导致方程组的病态问题，Wp和WE反映了分析者的工程经验和判断，对于可靠性较高的参数可以赋以较大的加权系数，反之亦然。拥有较大的加权系数的参数将对修正结果产生较大的影响，而有较大的加权系数的参数将具有较小的修正量。为此构建目标函数的优化问题为[13-14]

minJ(Δp)=ΔpTWpΔp+(SΔp-ΔR)TWE(SΔp-ΔR)

(29)

式中Wp和WE为加权矩阵。可得到

Δp=Wp-1ST(WE-1+SWp-1ST)ΔR。

(30)

为了使修正后的参数具有物理意义，引入设计参数修正向量的上、下限值p1、p2，然后令：

p1≤Δp≤p2，

(31)

作为式(29)的约束条件，那么由(29)和(31)构建的数学模型变成了求解一个约束优化问题。还可以将这一约束优化问题转化为二次规划问题来求解，即

在求解出Δp后，记为Δp0，假设修正参数的初始值为p0，构造迭代形式，使在每一步迭代过程中，都将使用更新后的修正参数，即

pk+1=pk+Δpk，

并将其代入有限元模型中计算新的解析解，进而得到新的残差向量，重复迭代过程直到残差向量为0或者小于某一设定阈值。例如在用式(28)得到修正向量和灵敏度矩阵后，在下一步的计算中，使用p0+Δp0代替p0，即可建立起使用第k步设计参数值表示第k+1步设计参数值的迭代公式

pk+1=pk+(SkTWESk)-1SkWEΔRk。

类似地，用式(30)构造的迭代公式为

pk+1=pk+Wp-1SkT(WE-1+SkWp-1SkT)ΔRk。

4 灵敏度矩阵的计算

一般地，若φ=(φ1(g)， … ，φN(g))T的每一维分量皆是向量g=(g1， … ，gm)T的函数，称φ=(φ1(g)， … ，φN(g))T为多元向量值函数。如果φ=(φ1(g)， … ，φN(g))T第i维分量φi(g1，g2，…，gm)的梯度向量为

那么，它的梯度矩阵为

如果向量φ=(φ1(g)， … ，φN(g))T对第j个设计参数gj(j=1，2，…，m)的一阶灵敏度为

那么多元向量值函数φ=(φ1(g)， … ，φN(g))T梯度矩阵还可简写为

(32)

由式(32)可知，梯度矩阵仍然可保持向量的形式，尽管它实质上是一个矩阵。

如果向量φ=(φ1(g)， … ，φN(g))T先对第j个设计参数gj再对第l个设计参数gl(j，l=1，2，…，m)的二阶灵敏度为

多元向量值函数φ=(φ1(g)， … ，φN(g))T的海森阵为

(33)

观察式(33)，多元向量值函数的海森阵只是一个矩阵形式而已，它已不能被认为是通常意义下的矩阵，因为它的每个元素仍是向量，但把它写成矩阵形式的好处在于可方便用于向量值函数在某点处作泰勒展开，因此仍可沿用海森阵的名称。

为实现快速而稳定的计算模态参数的灵敏度，目前常用的方法是模态展开法[15-17]，其思想基础是先构建特征空间的基底，因为模态灵敏度也是隶属于该特征空间的向量，因此必然可表示为此特征空间基底的某一线性组合形式。本文仅以对称单频系统为例，给出实、复模态参数的灵敏度算法。

首先是实频率的一阶灵敏度算法[18]

其中

实模态的一阶灵敏度算法：

其中

(k≠i，k=1，…，N)。

其次是复频率的一阶灵敏度算法[19]

其中

复模态的一阶灵敏度算法：

其中

有关模态的MAC值的灵敏度矩阵[20]、频响函数灵敏度矩阵[21-23]及反共振频率灵敏度矩阵[20]的计算方法请参见文献。

5 结语

综上所述，基于灵敏度矩阵方程所建立的模型修正方法的一般步骤如下：1)针对结构确定修正对象f，取得设计参数p的初始值p(o)；2)利用结构识别及参数估计方法，采集修正对象f的实测数据；3)利用模型缩聚技术和扩展技术对修正对象的实测数据进行数据处理，获得修正对象的实验测量值f(p(u))；4)根据含有设计参数的图纸或基线模型建立结构的有限元初始模型；5)利用有限元初始模型计算修正对象的解析值f(p(o))；6)计算修正对象的灵敏度矩阵S(p(o))；7)建立灵敏度矩阵方程(见式(25))；8)利用广义逆技术，非线性最小二乘技术或正则化方法求解灵敏度矩阵方程，解出设计参数的修正量Δp；9)获得设计参数的修正值p(o)+Δp；10)返回步4)获得有限元修正模型，直至有限元修正模型满足一定准则，结束并输出有限元修正模型。

而基于残差型目标函数所建立的模型修正方法的步骤的前6步与上述方法一致，7)是建立残差型目标函数(见式(5)～(21))；8)利用高斯牛顿算法求解非线性最小二乘优化模型，最后解出修正参数的稳定解p*。