基于模態(tài)參數(shù)靈敏度的模型修正方法研究

張 淼

(長春工程學院理學院，長春 130012)

結構分析的基本目標之一是公式化能由實測數(shù)據(jù)檢驗的結構有限元分析模型。這個模型常常不能一開始就產生可靠的基頻和模態(tài)振型，為此引入實測數(shù)據(jù)來調整此分析模型，直到分析與實驗數(shù)據(jù)得到匹配。修正后的有限元模型能夠更好地反映結構的動力特性。為此，需要從不完全的實測模態(tài)參數(shù)數(shù)據(jù)中建立剛度和質量，甚至阻尼矩陣，而從實驗中測量出來的數(shù)據(jù)在數(shù)量上必須超過在結構中識別出來的參數(shù)的數(shù)量，因此盡可能多地應用實測數(shù)據(jù)是十分重要的。如果只用頻率數(shù)據(jù)，可識別參數(shù)的數(shù)量就十分有限，那么測量的模態(tài)數(shù)據(jù)也應該包含在修正過程之中。數(shù)據(jù)和建模精度都是修正過程中需考慮的重要因素，當然模態(tài)測量的精度不如頻率的測量精度，因此任何一個修正程序如果使用最小二乘法，就必須包含數(shù)據(jù)精度的權重。

傳統(tǒng)的模型修正方法大致分為兩種，即矩陣元素型和設計參數(shù)型，靈敏度分析在有限元模型修正過程中常常被使用，有時甚至成為決定模型修正成敗的關鍵因素。近年來由于靈敏度分析以及非線性優(yōu)化技術的不斷進步，使得設計參數(shù)型模型修正方法的發(fā)展成效更為顯著一些?；陟`敏度的模型修正方法，是目前模型修正領域中發(fā)展較為活躍的技術之一。它的應用經常融合了殘差型和靈敏度矩陣方程型的目標函數(shù)，本文的目的是試圖從數(shù)學角度揭示這兩種模型修正技術之間的內在聯(lián)系，并充分展示它們的技術輪廓。

1 基本理論[1]

描述自由度為N的線性阻尼離散系統(tǒng)的自由振動方程為

(1)

(2)

令λi=-ωi2，則式(2)可化為

(K+λiM)vi=0。

可以計算得到系統(tǒng)的實頻率λ1，λ2，…，λN及實模態(tài)v1，v2，…，vN。記V=[v1，…，vN]為無阻尼規(guī)范化振型矩陣，那么此時模態(tài)質量和模態(tài)剛度矩陣分別為

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

對于經典阻尼系統(tǒng)，模態(tài)阻尼矩陣為

VTCV=diag(c1，…，cN)。

考慮阻尼時的系統(tǒng)極點及復模態(tài)對(si，ui)(i=1，2，…，2N)滿足方程

對于N自由度振動系統(tǒng)，特征方程det[s2M+sC+K]=0有2N個呈復共軛對出現(xiàn)的特征值s1，s2，…，s2N(其中si+1為si的共軛(i=1，3，…，2N-1))，稱為系統(tǒng)的極點或復頻率。這些頻率對應著一組呈復共軛對出現(xiàn)特征向量ui∈CN稱為與si相對應的第i個模態(tài)向量。將u1，u2，…，u2N(其中ui+1為ui的共軛(i=1，3，…，2N-1))，稱為復模態(tài)。復模態(tài)和復頻率統(tǒng)稱為復模態(tài)參數(shù)。復模態(tài)的正交規(guī)范化的形式有很多，對稱系統(tǒng)常用的形式為

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

其中狀態(tài)向量矩陣為Φ=[φ1，φ2，…，φ2N]，狀態(tài)向量為φi=[uiλiui]T(i=1，2，…，2N)，且

2 基于殘差型目標函數(shù)的模型修正方法

文獻[2]中，為了修正有限元模型參數(shù)，基于某些修正參數(shù)，把測量得到的模態(tài)數(shù)據(jù)作為參考數(shù)據(jù)，通過解一個非線性最小二乘問題來優(yōu)化一個價值函數(shù)，它是由實測和計算得到的解析頻率之差的平方和構成的。特別指出有限元模型修正是一個非線性參數(shù)估計問題，因為模態(tài)數(shù)據(jù)是未知結構性質的非線性函數(shù)：

其中

例如若選擇板的長度作為修正參數(shù)，為了解非線性最小二乘問題，需要確定板長度的上下限，然后作為初始值，通過增加回歸階數(shù)的回歸算法用來解插值多項式的回歸系數(shù)，即通過已知的模態(tài)頻率來擬合一條插值曲線，例如線性回歸擬合?；诔跏加嬎銋?shù)值，用修正參數(shù)L的已知值和計算出來的有限元頻率進行一個高階擬合，優(yōu)化價值函數(shù)，有限元頻率用插值多項式替換，當參數(shù)修正后的頻率殘差變得比使用者定義的準則還小的時候，令修正算法終止，從而價值函數(shù)有如式(3)的形式：

(3)

其中

aij是插值多項式的系數(shù)。

最小二乘問題允許考慮殘差的權重，依據(jù)這些殘差的重要性和噪音的數(shù)量，加權重來區(qū)分不同的數(shù)據(jù)系列能使方法變得更實用，但同時要求工程師能提供正確的權重，例如如式(4)模型

(4)

(5)

其中wj是殘差rj的權因子。權因子可以依據(jù)測量的隨機性質定義，也可以把相關測量的方差的倒數(shù)作為權因子。

文獻[2]還提出了另一種價值函數(shù)的定義方法

其中

文獻[3]中提出了有限元模型修正的多目標優(yōu)化模型，一般的公式如式(6)：

(6)

式中：θ是修正參數(shù)向量；Fi(θ)是目標函數(shù)，它是一個數(shù)量值函數(shù)；n是子目標函數(shù)的個數(shù)；wi是相應的子目標函數(shù)的權因子；gj(θ)與hk(θ)是等式和不等式約束；θL與θU是修正參數(shù)的上下限，權因子在平衡子目標函數(shù)的重要性方面具有重要的作用。它的選擇依賴于使用者的經驗和專業(yè)技能。在此文獻中用遺傳算法求解多目標規(guī)劃問題。

文獻[4]嘗試了幾種不同的權和的表達形式，并建議最平衡的一種可以選擇為模態(tài)和頻率如式(7)形式的加權殘差：

(7)

文獻[5]在基于靈敏度分析的模型修正程序中給出的目標函數(shù)是一種聯(lián)合了模態(tài)和頻率的優(yōu)化問題：

(8)

文獻[6]中目標函數(shù)f(θ)：Rn→R是一個普通的加權最小二乘問題：

(9)

rj(θ)：Rn→R(j=1，…，m)是殘差向量，它是關于修正變量θ∈Rn的非線性函數(shù)；r(θ)=[r1(θ)，r2(θ)，…，rm(θ)]T；W=diag(w1，…，wj，…，wm)，W∈Rm×m。殘差向量的第一部分包含無阻尼固有頻率的解析值及實測值之差：

(10)

(11)

文獻[7]用兩個模態(tài)參數(shù)，即模態(tài)頻率和模態(tài)振型來組成殘差向量，構建目標函數(shù)，一般來說它是非線性最小二乘問題：

(12)

r：Rn→Rm，θ∈Rn

式中：θ代表未知的修正參數(shù)組成的向量；r代表由頻率殘差rf和模態(tài)振型殘差rm組成的殘差向量：

(13)

(14)

式中W是權矩陣。

(15)

式中rj(θ)代表測量與解析模態(tài)參數(shù)之差，它是一個關于優(yōu)化變量θ∈Rn的非線性函數(shù)。為了獲得唯一解，m(殘差的個數(shù))要多于未知數(shù)θ的元個數(shù)n，殘差向量r(θ)：Rn→Rm，分成頻率殘差rf，模態(tài)振型的殘差rs和一個模態(tài)質量的殘差向量rm：

(16)

(17)

(18)

(19)

文獻[9]中指出基于靈敏度的有限元模型修正方法是最頻繁用于模型修正的方法。模型修正的優(yōu)化問題為：

(20)

這種方法直接比較頻率與模態(tài)振型：

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

(21)

所以加權特征值殘差可以表示為

在模型修正中，模態(tài)振型都是質量規(guī)范化后的，當環(huán)境振動用于提取模態(tài)參數(shù)，實驗模態(tài)保持未縮放的，因此模態(tài)振型殘差定義為，

這個因子是可變的，給出頻率與模態(tài)振型的重要性。實驗頻率一般是最精確的實驗數(shù)據(jù)，而實驗模態(tài)振型具有更多的噪音，會減少優(yōu)化問題的穩(wěn)定性，但能提供更多局部的信息，因此一般適當加權是很必要的。

在文獻[5]～[9]中都采用的是解非線性最小二乘優(yōu)化問題的信賴域牛頓高斯算法，其優(yōu)化模型分別為本文中的公式(8)～(21)，在每步迭代中搜索步長限定在一個信賴域內，避免了不希望出現(xiàn)的大步長。f(θ)可以用一個二次型來近似，它來自于泰勒級數(shù)截斷：

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

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

式中J(θ)是由ri(θ)的一階偏導數(shù)構成的靈敏度矩陣。信賴域牛頓方法是一種迭代方法，是將高斯牛頓方法用于信賴域中，來提高收斂性，在信賴域方法中，算法構成了一個價值函數(shù)mk，它的行為類似于目標函數(shù)f，在當前點θk附近，并圍繞θk得到一個區(qū)域，這個模型是可信的，盡管海森陣可以用J(θ)TJ(θ)近似。標準的高斯牛頓方法是相當穩(wěn)定的，而且通過在信賴域中實施，可以加快收斂。

3 基于靈敏度矩陣方程的模型修正方法

將g個修正對象f1，f2，…，fg(可以包括實或復的模態(tài)參數(shù)，質量、剛度矩陣，反共振頻率，模態(tài)的MAC值，頻響函數(shù)等，也可以是這些對象的組合)均看作是m個設計參數(shù)p=(p1，p2，…，pm)T的函數(shù)f(p)，在設計參數(shù)(可以取為材料常數(shù)，幾何參數(shù)或邊界條件等)的初始取值位置作一階泰勒展開，使問題線性化：

(22)

(23)

并認為靈敏度矩陣

(24)

為已知或者可求，則式(23)還可表示為

SΔp=ΔR，

(25)

式中Δp=(Δp1，Δp2，…，Δpm)T為待解修正向量；ΔR=f(p(u))-f(p(o))為殘差向量。求解方程(25)，即可解出設計參數(shù)的修正向量Δ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為實模態(tài)時，只需將式(26)～(27)中的λ=(λ1，λ2，…，λN)T換成[v]=(v1，v2，…，vN)T；當f為復模態(tài)時，只需將式(26)～(27)中的λ=(λ1，λ2，…，λN)T換成[u]=(u1，u2，…，u2N)T，其中[u]和[v]表示由相應的復模態(tài)和實模態(tài)向量構成的矩陣。如果在實測數(shù)據(jù)不完備的情況下，還可以將N或2N換成實測數(shù)據(jù)的個數(shù)t，此時方程(25)是超定的。

在許多實際問題中所遇到的方程組Ax=b的系數(shù)矩陣往往是奇異的或是m×n任意矩陣(一般m≠n)，顯然不存在通常的逆矩陣A-1來求解方程組，這就促使人們推廣逆矩陣的概念，引進某種具有普通逆矩陣類似的性質的矩陣G，使得方程組的解仍可表示為x=Gb。而廣義逆理論能夠把相容線性方程組的一般解、極小范數(shù)解以及矛盾方程最小二乘解、極小范數(shù)最小二乘解(最佳逼近解)全部統(tǒng)一起來。

1)當方程組相容時，若系數(shù)矩陣A∈Cn×n，且非奇異，則有唯一解x=A-1b。

若靈敏度矩陣方程(25)中的系數(shù)矩陣的秩n與設計參數(shù)的個數(shù)m相同，則可以直接解得：

Δp=S-1ΔR，

在很多情況下靈敏度矩陣方程中的系數(shù)矩陣的秩n與設計參數(shù)的個數(shù)m并不相同，當n>m時，引入廣義逆來求解得[11]：

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

當n

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

從靈敏度矩陣方程(25)出發(fā)，看成是最小化殘差的優(yōu)化模型，令目標函數(shù)為[12]

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

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

(28)

然而通常由于靈敏度矩陣S會接近奇異，可以同時引入了兩種加權矩陣Wp和WE，來克服因靈敏度矩陣S病態(tài)而導致方程組的病態(tài)問題，Wp和WE反映了分析者的工程經驗和判斷，對于可靠性較高的參數(shù)可以賦以較大的加權系數(shù)，反之亦然。擁有較大的加權系數(shù)的參數(shù)將對修正結果產生較大的影響，而有較大的加權系數(shù)的參數(shù)將具有較小的修正量。為此構建目標函數(shù)的優(yōu)化問題為[13-14]

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

(29)

式中Wp和WE為加權矩陣?？傻玫?/p>

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

(30)

為了使修正后的參數(shù)具有物理意義，引入設計參數(shù)修正向量的上、下限值p1、p2，然后令：

p1≤Δp≤p2，

(31)

作為式(29)的約束條件，那么由(29)和(31)構建的數(shù)學模型變成了求解一個約束優(yōu)化問題。還可以將這一約束優(yōu)化問題轉化為二次規(guī)劃問題來求解，即

在求解出Δp后，記為Δp0，假設修正參數(shù)的初始值為p0，構造迭代形式，使在每一步迭代過程中，都將使用更新后的修正參數(shù)，即

pk+1=pk+Δpk，

并將其代入有限元模型中計算新的解析解，進而得到新的殘差向量，重復迭代過程直到殘差向量為0或者小于某一設定閾值。例如在用式(28)得到修正向量和靈敏度矩陣后，在下一步的計算中，使用p0+Δp0代替p0，即可建立起使用第k步設計參數(shù)值表示第k+1步設計參數(shù)值的迭代公式

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的函數(shù)，稱φ=(φ1(g)， … ，φN(g))T為多元向量值函數(shù)。如果φ=(φ1(g)， … ，φN(g))T第i維分量φi(g1，g2，…，gm)的梯度向量為

那么，它的梯度矩陣為

如果向量φ=(φ1(g)， … ，φN(g))T對第j個設計參數(shù)gj(j=1，2，…，m)的一階靈敏度為

那么多元向量值函數(shù)φ=(φ1(g)， … ，φN(g))T梯度矩陣還可簡寫為

(32)

由式(32)可知，梯度矩陣仍然可保持向量的形式，盡管它實質上是一個矩陣。

如果向量φ=(φ1(g)， … ，φN(g))T先對第j個設計參數(shù)gj再對第l個設計參數(shù)gl(j，l=1，2，…，m)的二階靈敏度為

多元向量值函數(shù)φ=(φ1(g)， … ，φN(g))T的海森陣為

(33)

觀察式(33)，多元向量值函數(shù)的海森陣只是一個矩陣形式而已，它已不能被認為是通常意義下的矩陣，因為它的每個元素仍是向量，但把它寫成矩陣形式的好處在于可方便用于向量值函數(shù)在某點處作泰勒展開，因此仍可沿用海森陣的名稱。

為實現(xiàn)快速而穩(wěn)定的計算模態(tài)參數(shù)的靈敏度，目前常用的方法是模態(tài)展開法[15-17]，其思想基礎是先構建特征空間的基底，因為模態(tài)靈敏度也是隸屬于該特征空間的向量，因此必然可表示為此特征空間基底的某一線性組合形式。本文僅以對稱單頻系統(tǒng)為例，給出實、復模態(tài)參數(shù)的靈敏度算法。

首先是實頻率的一階靈敏度算法[18]

其中

實模態(tài)的一階靈敏度算法：

其中

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

其次是復頻率的一階靈敏度算法[19]

其中

復模態(tài)的一階靈敏度算法：

其中

有關模態(tài)的MAC值的靈敏度矩陣[20]、頻響函數(shù)靈敏度矩陣[21-23]及反共振頻率靈敏度矩陣[20]的計算方法請參見文獻。

5 結語

綜上所述，基于靈敏度矩陣方程所建立的模型修正方法的一般步驟如下：1)針對結構確定修正對象f，取得設計參數(shù)p的初始值p(o)；2)利用結構識別及參數(shù)估計方法，采集修正對象f的實測數(shù)據(jù)；3)利用模型縮聚技術和擴展技術對修正對象的實測數(shù)據(jù)進行數(shù)據(jù)處理，獲得修正對象的實驗測量值f(p(u))；4)根據(jù)含有設計參數(shù)的圖紙或基線模型建立結構的有限元初始模型；5)利用有限元初始模型計算修正對象的解析值f(p(o))；6)計算修正對象的靈敏度矩陣S(p(o))；7)建立靈敏度矩陣方程(見式(25))；8)利用廣義逆技術，非線性最小二乘技術或正則化方法求解靈敏度矩陣方程，解出設計參數(shù)的修正量Δp；9)獲得設計參數(shù)的修正值p(o)+Δp；10)返回步4)獲得有限元修正模型，直至有限元修正模型滿足一定準則，結束并輸出有限元修正模型。

而基于殘差型目標函數(shù)所建立的模型修正方法的步驟的前6步與上述方法一致，7)是建立殘差型目標函數(shù)(見式(5)～(21))；8)利用高斯牛頓算法求解非線性最小二乘優(yōu)化模型，最后解出修正參數(shù)的穩(wěn)定解p*。