基于近似似然的頻響函數不確定性模型修正*

鄧振鴻， 張保強， 蘇國強， 郭勤濤

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

引 言

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

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

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

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

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

對于多自由度振動系統，其動力學方程為

(1)

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

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

(2)

其中：ω為頻率。

定義系統的動剛度矩陣

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

(3)

動剛度矩陣的逆矩陣即為系統的位移頻響函數矩陣

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

(4)

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

同理，可以得到系統的速度頻響函數和加速度頻響函數。

當系統參數θ存在不確定性的情況下，試驗測試的頻響函數和仿真計算頻響函數有如下關系

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

(5)

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

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

2.1 貝葉斯理論

(6)

因此有

(7)

可以看出，當先驗信息給定時，似然函數的值越大，參數的后驗概率也越大，意味著該參數接近目標值的可能性越大。

2.2 適用于頻響函數模型修正的似然函數

(8)

其中：e為實測值和仿真值的殘差；Σ為協方差矩陣。

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

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

2.3 DREAM算法

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

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

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

1) 結合測試得到的頻響函數數據，統計出所關注頻率下的所有試驗頻響函數的上、下限以及中值；

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

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

4) DREAM算法抽樣產生參數值，并代入模型計算頻響；

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

6) 當計算次數小于等于設定的迭代次數時，回到步驟4，否則停止迭代計算，統計出所有使似然值取非零值的參數樣本。

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

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

3 算例驗證

3.1 三自由度算例

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

圖2 三自由度系統Fig.2 3DOF vibration system

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

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

表1 參數修正前后值及誤差對比

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

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

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

3.2 H型非對稱梁算例

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

表2 材料參數及其不確定性

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

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

表 3給出了參數修正前后的上、下限及誤差對比，參數E1的誤差由[-8.47%,10.42%]降低至[0.20%,0.79%]，參數E2的誤差由[-8.47%,10.42%]降低至[0.29%,-0.11%]，參數ζ的誤差由[-28.57%,25.00%]降低至[0.81%,0.41%]，表明參數的修正精度較高。圖 6為修正前后的FRF對比圖，分別給出了修正前后FRF邊界以及目標值?？梢钥闯鱿啾扔谛拚埃拚蟮腇RF邊界縮小，并與目標值的邊界相重合。圖 7給出了其他響應點(單元10，11，12交點處節點)修正后FRF的預測對比，修正后的FRF邊界與目標值的邊界也基本吻合。進一步對所有試驗測試的頻響函數增加2%水平的高斯白噪聲，以考慮試驗噪聲對修正結果的影響。參數修正結果見表4，可以看出，相比于修正前，不確定性參數的誤差下降明顯，修正后不確定性參數也接近于真實值，說明在一定水平噪聲情況下，仍具有不錯的修正效果。綜上所述，筆者所提

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

圖7 其他響應點修正前后FRF預測對比圖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 參數修正前后值及誤差對比

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

4 結 論

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

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

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

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