張 淼

(長(zhǎng)春工程學(xué)院理學(xué)院，長(zhǎng)春 130012)

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

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

1 基本理論[1]

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

(1)

(2)

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

(K+λiM)vi=0。

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

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

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

VTCV=diag(c1，…，cN)。

考慮阻尼時(shí)的系統(tǒng)極點(diǎn)及復(fù)模態(tài)對(duì)(si，ui)(i=1，2，…，2N)滿足方程

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

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

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

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

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

其中

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

(3)

其中

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

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

(4)

(5)

其中wj是殘差rj的權(quán)因子。權(quán)因子可以依據(jù)測(cè)量的隨機(jī)性質(zhì)定義，也可以把相關(guān)測(cè)量的方差的倒數(shù)作為權(quán)因子。

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

其中

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

(6)

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

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

(7)

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

(8)

文獻(xiàn)[6]中目標(biāo)函數(shù)f(θ)：Rn→R是一個(gè)普通的加權(quán)最小二乘問題：

(9)

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

(10)

(11)

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

(12)

r：Rn→Rm，θ∈Rn

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

(13)

(14)

式中W是權(quán)矩陣。

(15)

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

(16)

(17)

(18)

(19)

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

(20)

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

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

(21)

所以加權(quán)特征值殘差可以表示為

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

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

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

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

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

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

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

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

(22)

(23)

并認(rèn)為靈敏度矩陣

(24)

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

SΔp=ΔR，

(25)

式中Δp=(Δp1，Δp2，…，Δpm)T為待解修正向量；ΔR=f(p(u))-f(p(o))為殘差向量。求解方程(25)，即可解出設(shè)計(jì)參數(shù)的修正向量Δp。

上文中指出，當(dāng)修正對(duì)象不同時(shí)，靈敏度矩陣S的計(jì)算技術(shù)也必然不同。例如當(dāng)f為實(shí)頻率時(shí)，其泰勒展開式為

式中λ=(λ1，λ2，…，λN)T，則靈敏度矩陣方程為

SΔp=ΔR，

其中

(26)

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

(27)

其中λ(p(u))為實(shí)頻率的實(shí)驗(yàn)測(cè)量值，而λ(p(o))為實(shí)頻率的解析計(jì)算值。

當(dāng)f為復(fù)頻率時(shí)，只需將式(26)～(27)中的λ=(λ1，λ2，…，λN)T換成s=(s1，s2，…，s2N)T；當(dāng)f為實(shí)模態(tài)時(shí)，只需將式(26)～(27)中的λ=(λ1，λ2，…，λN)T換成[v]=(v1，v2，…，vN)T；當(dāng)f為復(fù)模態(tài)時(shí)，只需將式(26)～(27)中的λ=(λ1，λ2，…，λN)T換成[u]=(u1，u2，…，u2N)T，其中[u]和[v]表示由相應(yīng)的復(fù)模態(tài)和實(shí)模態(tài)向量構(gòu)成的矩陣。如果在實(shí)測(cè)數(shù)據(jù)不完備的情況下，還可以將N或2N換成實(shí)測(cè)數(shù)據(jù)的個(gè)數(shù)t，此時(shí)方程(25)是超定的。

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

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

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

Δp=S-1ΔR，

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

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

當(dāng)n

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

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

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

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

(28)

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

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

(29)

式中Wp和WE為加權(quán)矩陣?？傻玫?/p>

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

(30)

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

p1≤Δp≤p2，

(31)

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

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

pk+1=pk+Δpk，

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

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

類似地，用式(30)構(gòu)造的迭代公式為

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

4 靈敏度矩陣的計(jì)算

一般地，若φ=(φ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對(duì)第j個(gè)設(shè)計(jì)參數(shù)gj(j=1，2，…，m)的一階靈敏度為

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

(32)

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

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

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

(33)

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

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

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

其中

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

其中

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

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

其中

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

其中

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

5 結(jié)語(yǔ)

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

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