鄧振鴻， 張保強(qiáng)， 蘇國強(qiáng)， 郭勤濤

(1.廈門大學(xué)航空航天學(xué)院 廈門，361005) (2.南京航空航天大學(xué)機(jī)電學(xué)院 南京，210016)

引 言

在許多工程問題中，所建立的仿真模型的準(zhǔn)確度對于結(jié)構(gòu)的設(shè)計(jì)分析至關(guān)重要，模型修正作為一種提高仿真模型可信度的手段，旨在借助試驗(yàn)數(shù)據(jù)來校準(zhǔn)仿真模型[1]。實(shí)際工程中，結(jié)構(gòu)設(shè)備在制造加工誤差、運(yùn)輸損傷及測試過程中的主客觀因素等的影響下，使得系統(tǒng)中存在諸多不確定性，從而造成多次試驗(yàn)的結(jié)果存在一定的分散性。如何準(zhǔn)確定位這些分散性即不確定性的來源并將不確定因素識別出來，已經(jīng)成為研究人員關(guān)注的重點(diǎn)。

不確定性模型修正技術(shù)是在傳統(tǒng)的確定性模型修正的基礎(chǔ)上，結(jié)合概率或非概率等統(tǒng)計(jì)學(xué)手段，對仿真模型進(jìn)行修正。Beck等[2]提出了貝葉斯框架下的不確定性模型修正。Collins等[3]在處理振動測量數(shù)據(jù)的隨機(jī)誤差時(shí)，引入線性靈敏度方法，以測量數(shù)據(jù)的協(xié)方差矩陣的逆作為權(quán)重，通過加權(quán)最小二乘法來獲得未知參數(shù)的統(tǒng)計(jì)量。文獻(xiàn)[4-6]基于文獻(xiàn)[2]的框架，發(fā)展了其他模型修正方法，并應(yīng)用在各種線性和非線性動力學(xué)系統(tǒng)中。文獻(xiàn)[7-8]發(fā)展了隨機(jī)模型修正理論，提出了梯度回歸方法，并應(yīng)用到一系列物理結(jié)構(gòu)中。文獻(xiàn)[9]比較了貝葉斯方法和攝動法，并將兩種不確定性修正方法應(yīng)用到DLR AIRMOD結(jié)構(gòu)上加以對比驗(yàn)證。方圣恩等[10]利用隨機(jī)響應(yīng)面模型的快速運(yùn)算特性來反演得到參數(shù)的均值和方差。陳喆等[11]研究了考慮試驗(yàn)?zāi)B(tài)不確定性的有限元模型修正方法。

上述不確定性模型修正研究均是以結(jié)構(gòu)振動的模態(tài)參數(shù)作為輸出響應(yīng)量進(jìn)行的，相比于模態(tài)參數(shù)，結(jié)構(gòu)振動的頻率響應(yīng)量，如頻響函數(shù)(frequency response function, 簡稱FRF)則包含了結(jié)構(gòu)更充分的信息，能夠更加準(zhǔn)確地反應(yīng)結(jié)構(gòu)的動力學(xué)特性。同時(shí)，采用頻響函數(shù)進(jìn)行模型修正可以有效避免由模態(tài)分析所引入的誤差。文獻(xiàn)[12-15]介紹了頻響函數(shù)驅(qū)動的模型修正研究，但均屬于確定性修正范疇內(nèi)的工作。Mares等[16]基于貝葉斯方法以三自由度質(zhì)量彈簧系統(tǒng)頻響函數(shù)的不確定性修正為例，比較了不同頻率點(diǎn)的選擇策略，發(fā)現(xiàn)靠近共振峰處的頻率點(diǎn)應(yīng)該被舍去，但仍沒有得到很理想的修正效果。Arora等[17]針對動力學(xué)系統(tǒng)中存在的密集模態(tài)和阻尼問題，提出了一種基于復(fù)頻響的參數(shù)不確定性修正方法。Hegde等[18]將單自由度二階系統(tǒng)頻響函數(shù)的幅值和相位分離，構(gòu)建出參數(shù)與響應(yīng)之間的線性模型，并采用極大似然估計(jì)方法反演出參數(shù)的均值和方差。Vakilzadhe等[19]在阻尼平衡模型校準(zhǔn)技術(shù)的基礎(chǔ)上，引入Bootstrating方法對試驗(yàn)頻響函數(shù)重復(fù)采樣，提出了一種新的隨機(jī)有限元模型修正框架，可用于處理含有噪聲的頻響函數(shù)模型修正問題。Machado等[20]將隨機(jī)場理論同基于靈敏度的頻響函數(shù)模型修正方法結(jié)合起來，對梁材料隨機(jī)場特性進(jìn)行估計(jì)。Zhang等[21]同時(shí)考慮了參數(shù)和模型形式引起的不確定性，基于貝葉斯方法對重復(fù)測量的頻響函數(shù)進(jìn)行了不確定性校準(zhǔn)研究。曹詩澤[22]對基于頻響函數(shù)的結(jié)構(gòu)有限元模型修正的不確定性量化方法進(jìn)行了研究，基于頻率響應(yīng)函數(shù)的概率模型，提出了頻響函數(shù)驅(qū)動的貝葉斯模型修正方法，并采用漸進(jìn)馬爾科夫蒙特卡洛算法進(jìn)行求解待修正參數(shù)的最優(yōu)解及后驗(yàn)概率密度函數(shù)。

傳統(tǒng)的貝葉斯模型修正通常假設(shè)不確定性參數(shù)服從正態(tài)分布，而實(shí)際工程中由于試驗(yàn)數(shù)據(jù)有限，其樣本的分布特性并不明顯，從而影響似然函數(shù)的選擇。因此，筆者在貝葉斯框架下，引入近似貝葉斯計(jì)算(approximate Bayesian computation, 簡稱ABC)理論[23]，采用區(qū)間作為參數(shù)的統(tǒng)計(jì)量，提出了一種適用于頻響函數(shù)的似然函數(shù)，使用DREAM算法對不確定性參數(shù)進(jìn)行識別，統(tǒng)計(jì)出參數(shù)的后驗(yàn)范圍，并分別采用了三自由度數(shù)值算例和H型非對稱梁的有限元模型修正算例加以驗(yàn)證。

1 基于頻響函數(shù)的不確定性模型修正理論

對于多自由度振動系統(tǒng)，其動力學(xué)方程為

(1)

其中：M，K，C分別為其質(zhì)量、剛度和阻尼矩陣；F(t)為激振力。

對式(1)進(jìn)行傅里葉變換得到

(2)

其中：ω為頻率。

定義系統(tǒng)的動剛度矩陣

Z(ω)=K-ω2M+jωC

(3)

動剛度矩陣的逆矩陣即為系統(tǒng)的位移頻響函數(shù)矩陣

H(ω)=Z(ω)-1=(K-ω2M+jωC)-1

(4)

其中：H(ω)為在頻率ω下的位移頻響函數(shù)矩陣，其元素Hij(ω)代表第j個(gè)自由度激勵(lì)下第i個(gè)自由度的頻率響應(yīng)。

同理，可以得到系統(tǒng)的速度頻響函數(shù)和加速度頻響函數(shù)。

當(dāng)系統(tǒng)參數(shù)θ存在不確定性的情況下，試驗(yàn)測試的頻響函數(shù)和仿真計(jì)算頻響函數(shù)有如下關(guān)系

He(ω)=Ha(ω;θ)+e

(5)

其中：He(ω)為試驗(yàn)頻響函數(shù)；Ha(ω;θ)為仿真頻響函數(shù)；e為二者的誤差，它是一個(gè)隨機(jī)變量。

2 貝葉斯框架下的模型修正

2.1 貝葉斯理論

(6)

因此有

(7)

可以看出，當(dāng)先驗(yàn)信息給定時(shí)，似然函數(shù)的值越大，參數(shù)的后驗(yàn)概率也越大，意味著該參數(shù)接近目標(biāo)值的可能性越大。

2.2 適用于頻響函數(shù)模型修正的似然函數(shù)

(8)

其中：e為實(shí)測值和仿真值的殘差；Σ為協(xié)方差矩陣。

然而在實(shí)際情況中，受限于試驗(yàn)樣本數(shù)量，其所服從的分布并不總能滿足正態(tài)分布。另外，當(dāng)在正態(tài)分布假設(shè)下選用頻響函數(shù)數(shù)據(jù)直接進(jìn)行修正時(shí)，必須選擇少量頻率點(diǎn)下的響應(yīng)進(jìn)行修正，否則可能出現(xiàn)響應(yīng)數(shù)據(jù)的協(xié)方差矩陣奇異的情況，而不同的頻率點(diǎn)選取策略則會對修正結(jié)果產(chǎn)生不同的影響[16]。近似貝葉斯估計(jì)是一種似然自由貝葉斯技術(shù)[23],該方法不假定數(shù)據(jù)服從某個(gè)特定的分布類型，而是通過設(shè)定一個(gè)容忍誤差ε，并計(jì)算某個(gè)參數(shù)下模型數(shù)據(jù)和觀測數(shù)據(jù)之間的距離d，當(dāng)d<ε時(shí)，就接受該參數(shù)，否則拒絕該參數(shù)?；谠摾碚摚P者采用一種近似似然函數(shù)，以實(shí)現(xiàn)頻響函數(shù)的不確定性模型修正。具體的似然函數(shù)表達(dá)式為

式(9)表明，當(dāng)試驗(yàn)與仿真之間的誤差滿足試驗(yàn)本身的變異水平時(shí)，令似然函數(shù)值等于一個(gè)相對較大的數(shù)，而不滿足時(shí)其值取零，從而實(shí)現(xiàn)對參數(shù)的接受和拒絕。同時(shí)為了保證繼續(xù)搜尋其他滿足收斂指標(biāo)的參數(shù)而不停留在該處，故將似然函數(shù)設(shè)置成某個(gè)較大數(shù)附近的隨機(jī)數(shù)，而非固定的值。

2.3 DREAM算法

對于參數(shù)后驗(yàn)概率的求解，由于Y關(guān)于θ的函數(shù)通常相對復(fù)雜或是以隱式存在，很難直接利用表達(dá)式進(jìn)行積分計(jì)算，因此對于后驗(yàn)分布的求解一般借助于馬爾科夫鏈蒙特卡洛(Markov chain Monte Carlo,簡稱MCMC)方法進(jìn)行抽樣。傳統(tǒng)的MCMC抽樣主要包括Metropolis算法、Metropolis-Hastings算法以及Gibbs抽樣等，然而上述算法普遍存在計(jì)算效率低、拒絕率高及依賴建議分布嚴(yán)重的問題。后續(xù)發(fā)展起來的延遲拒絕MCMC算法能夠自適應(yīng)地改變建議分布，加快其收斂效率，但可能造成得到的后驗(yàn)分布與真實(shí)分布存在偏差。

DREAM算法[24]是一種改進(jìn)的MCMC算法，其全稱為差分進(jìn)化自適應(yīng)Metropolis算法。DREAM算法具備多鏈并行抽樣的能力，每條并行鏈在產(chǎn)生不同初始樣本后，通過差分進(jìn)化方程來產(chǎn)生新樣本，在搜索過程中能自適應(yīng)地調(diào)整步長和方向，從而對多個(gè)全局最優(yōu)區(qū)域進(jìn)行有效搜索，有利于樣本從局部最優(yōu)點(diǎn)跳到全局最優(yōu)點(diǎn)，擁有高效率的計(jì)算能力[25-26]。筆者采用 DREAM-Matlab工具箱[27]求解不確定性參數(shù)的后驗(yàn)分布。

本研究所提模型修正方法具體包含如下步驟：

1) 結(jié)合測試得到的頻響函數(shù)數(shù)據(jù)，統(tǒng)計(jì)出所關(guān)注頻率下的所有試驗(yàn)頻響函數(shù)的上、下限以及中值；

2) 根據(jù)試驗(yàn)數(shù)據(jù)的上、下限及中值計(jì)算出數(shù)據(jù)的離散度，將其作為容忍誤差向量；

3) 設(shè)定參數(shù)的初始范圍及迭代次數(shù)；

4) DREAM算法抽樣產(chǎn)生參數(shù)值，并代入模型計(jì)算頻響；

5) 根據(jù)式(9)計(jì)算近似似然函數(shù)值L，若L>0則接受該參數(shù)，L=0則拒絕該參數(shù)；

6) 當(dāng)計(jì)算次數(shù)小于等于設(shè)定的迭代次數(shù)時(shí)，回到步驟4，否則停止迭代計(jì)算，統(tǒng)計(jì)出所有使似然值取非零值的參數(shù)樣本。

上述步驟的流程圖如圖 1所示。

圖1 修正過程流程圖Fig.1 The flow chart of model updating process

3 算例驗(yàn)證

3.1 三自由度算例

某三自由度振動系統(tǒng)見圖 2，已知m1=m2=m3=1.0 kg，k3=k4=1.0 N/m，k6=3.0 N/m。其他3個(gè)剛度系數(shù)k1，k2，k5為待識別的不確定性參數(shù)，表現(xiàn)為在區(qū)間[0.8,1.2]N/m范圍內(nèi)波動。

圖2 三自由度系統(tǒng)Fig.2 3DOF vibration system

本算例在參數(shù)的目標(biāo)區(qū)間范圍內(nèi)采用拉丁超立方抽樣100次得到的頻響函數(shù)作為目標(biāo)值，認(rèn)為參數(shù)的初始值在區(qū)間[0.4,2.4] N/m內(nèi)，采用DREAM算法對參數(shù)的后驗(yàn)分布進(jìn)行更新。圖3給出了5 000次抽樣計(jì)算的參數(shù)迭代結(jié)果，可以看出參數(shù)在約500次迭代后趨于收斂，將收斂后參數(shù)樣本的最大值和最小值統(tǒng)計(jì)出來，并與目標(biāo)值進(jìn)行比較，結(jié)果見表 1。

圖3 待修正參數(shù)迭代收斂過程Fig.3 Iterative convergence process of modified parameters

表1 參數(shù)修正前后值及誤差對比

從表1可以看出，修正后參數(shù)的上限值誤差由原來的100%減小到3%以下，下限誤差的絕對值由50%降低至不超過8.75%，可見修正效果顯著。

圖 4為修正前后頻響函數(shù)的對比，分別給出了初始參數(shù)下得到的FRF值、目標(biāo)FRF值和修正后的FRF邊界值。從圖 4可以看出，修正后FRF邊界能夠較好地包含所有的目標(biāo)FRF曲線，相比于修正前，不確定性的范圍大大縮小。

圖4 修正前后FRF對比圖Fig.4 Comparison of FRF before and after updating

3.2 H型非對稱梁算例

結(jié)合某H型非對稱梁有限元模型修正算例來說明筆者提出方法的可行性。H型非對稱梁及其尺寸如圖 5所示，將其劃分為12個(gè)單元，材料為鋁。在加工過程中，兩側(cè)的豎梁和橫梁分別選用了兩批不同的鋁材料，同時(shí)由于安裝過程中的損傷，其結(jié)構(gòu)剛度受到了影響，主要體現(xiàn)在彈性模量的不確定性中，具體的材料參數(shù)見表2。

表2 材料參數(shù)及其不確定性

圖5 非對稱H型梁結(jié)構(gòu)及其尺寸示意圖(單位：mm)Fig.5 The unsymmetrical H-shaped structure(unit:mm)

本算例中選取E1和E2作為待修正參數(shù)，另外考慮系統(tǒng)阻尼的不確定性，將模態(tài)阻尼比也作為不確定性參數(shù)同時(shí)進(jìn)行修正，采用Matlab和NASTRAN聯(lián)合編程，在參數(shù)目標(biāo)值范圍內(nèi)執(zhí)行100次拉丁超立方抽樣獲得的頻響函數(shù)作為目標(biāo)值來代替試驗(yàn)。將參數(shù)的初始區(qū)間定義在[60 000, 80 000]MPa范圍內(nèi)。在圖 5中標(biāo)出的激勵(lì)點(diǎn)處施加單位加速度激勵(lì)，選取的頻率帶寬范圍為100～800 Hz，模態(tài)阻尼比ζ的目標(biāo)區(qū)間設(shè)定為[0.028,0.032]，初始區(qū)間為[0.02,0.04]。以圖 5所示的響應(yīng)點(diǎn)處的位移頻響函數(shù)為目標(biāo)，對參數(shù)進(jìn)行修正。

表 3給出了參數(shù)修正前后的上、下限及誤差對比，參數(shù)E1的誤差由[-8.47%,10.42%]降低至[0.20%,0.79%]，參數(shù)E2的誤差由[-8.47%,10.42%]降低至[0.29%,-0.11%]，參數(shù)ζ的誤差由[-28.57%,25.00%]降低至[0.81%,0.41%]，表明參數(shù)的修正精度較高。圖 6為修正前后的FRF對比圖，分別給出了修正前后FRF邊界以及目標(biāo)值。可以看出相比于修正前，修正后的FRF邊界縮小，并與目標(biāo)值的邊界相重合。圖 7給出了其他響應(yīng)點(diǎn)(單元10，11，12交點(diǎn)處節(jié)點(diǎn))修正后FRF的預(yù)測對比，修正后的FRF邊界與目標(biāo)值的邊界也基本吻合。進(jìn)一步對所有試驗(yàn)測試的頻響函數(shù)增加2%水平的高斯白噪聲，以考慮試驗(yàn)噪聲對修正結(jié)果的影響。參數(shù)修正結(jié)果見表4，可以看出，相比于修正前，不確定性參數(shù)的誤差下降明顯，修正后不確定性參數(shù)也接近于真實(shí)值，說明在一定水平噪聲情況下，仍具有不錯(cuò)的修正效果。綜上所述，筆者所提

表4 參數(shù)修正前后值及誤差對比(2%噪聲)

圖7 其他響應(yīng)點(diǎn)修正前后FRF預(yù)測對比圖Fig.7 Comparison of FRF prediction before and after updating of other response point

圖6 修正前后FRF對比圖Fig.6 Comparison of FRF before and after updating

表3 參數(shù)修正前后值及誤差對比

方法對有限元模型的不確定性參數(shù)修正也具有適用性。

4 結(jié) 論

1) 在貝葉斯模型修正框架下，引入近似貝葉斯計(jì)算方法，提出了一種近似似然函數(shù)，適用于頻響函數(shù)不確定性模型修正，分別結(jié)合三自由度數(shù)值算例以及非對稱梁的有限元模型修正算例進(jìn)行了驗(yàn)證。

2) 修正過程中可直接應(yīng)用離散化后的頻響函數(shù)數(shù)據(jù)，不需要人為地選擇頻率點(diǎn)來避開共振峰處的響應(yīng)。

3) 在參數(shù)分布不明確的情況下，采用區(qū)間來描述參數(shù)的不確定性，并引入近似貝葉斯計(jì)算理論，使用近似似然函數(shù)對參數(shù)進(jìn)行篩選，避免了主觀假設(shè)分布類型帶來的偏差。

4) 算例結(jié)果表明，修正后的模型與目標(biāo)值更加接近，且在一定噪聲水平下仍能保持不錯(cuò)的修正效果，對于工程實(shí)際有一定的參考價(jià)值。