張雪萍, 彭珍瑞, 張亞峰

(蘭州交通大學(xué) 機電工程學(xué)院，蘭州 730070)

1 引 言

在工程領(lǐng)域中，一個能夠準(zhǔn)確表征結(jié)構(gòu)行為的精確有限元模型非常重要[1]。有限元模型修正已成為振動領(lǐng)域研究的熱點之一[2]。但在實際工程應(yīng)用中，由于結(jié)構(gòu)幾何尺寸、材料變異性和模型誤差等因素的影響，不確定性廣泛存在于實際結(jié)構(gòu)的參數(shù)和響應(yīng)中[3]。

為了使動力學(xué)模型能夠盡可能地反映實際情況，國內(nèi)外開始運用隨機系統(tǒng)理論建立不確定性模型，在有限元模型修正過程中結(jié)合概率方法進(jìn)行分析是十分必要的[4]。不確定性模型修正方法中基于貝葉斯統(tǒng)計理論的研究是一個熱點。Beck等[5]將貝葉斯理論引入到模型修正中，并明確了修正的基本思路。韓芳等[6]提出了一種基于信息融合和貝葉斯理論的模型修正方法。但是，貝葉斯方法存在無法直接積分求取參數(shù)的似然函數(shù)和后驗概率公式的積分常數(shù)等缺點。因此，貝葉斯方法的參數(shù)后驗概率密度常用馬爾科夫鏈蒙特卡羅MCMC(Markov Chain Monte Carlo)方法來推斷。Beck等[7]又提出了一種基于MH(Metropolis -Hastings)算法的自適應(yīng)MCMC方法。Green[8]提出了融合改進(jìn)退火算法的MCMC方法，避免了目標(biāo)函數(shù)陷入局部最優(yōu)。Lam等[9]使用多個CPU生成多個馬爾可夫鏈(而不是單個鏈)的并行MCMC方法來提高M(jìn)CMC方法的計算效率。蔣偉等[10]在標(biāo)準(zhǔn)MCMC方法的基礎(chǔ)上，引入差分進(jìn)化算法，為提高不確定性模型修正中的計算精度提供了一種新手段。然而，在參數(shù)維度較高的情況下，雖然MCMC算法采樣所得的樣本具有一定的收斂性，但馬爾可夫鏈出現(xiàn)大量的停滯階段，樣本的接受率很低[11]。彭珍瑞等[12]引入了布谷鳥算法中新鳥巢更新的思想改進(jìn)MCMC算法。以上文獻(xiàn)大多采用結(jié)構(gòu)的響應(yīng)為頻率和位移等進(jìn)行模型修正，而在工程實踐中會出現(xiàn)選擇的響應(yīng)對局部狀態(tài)參數(shù)不敏感的情況，從而在修正過程中不能得到準(zhǔn)確、合理的應(yīng)力分布。因此，在有限元模型修正方法中，選擇一種能夠同時反映結(jié)構(gòu)局部狀態(tài)與全局特征的響應(yīng)至關(guān)重要。同時，在有限元模型修正中，由于所用代理模型精度不高也會使修正結(jié)果不理想。

基于上述背景，本文引入花朵授粉算法FPA(Flower Pollination Algorithm)與標(biāo)準(zhǔn)MCMC算法融合，解決在有效提高新樣本接受率的同時降低修正參數(shù)的相對誤差；且由于應(yīng)變模態(tài)包含結(jié)構(gòu)全局的頻率信息和能夠表征結(jié)構(gòu)局部狀態(tài)的應(yīng)變振型信息，因此選用應(yīng)變模態(tài)作為結(jié)構(gòu)響應(yīng)進(jìn)行模型修正，并使用蝙蝠算法BA(Bat Algorithm)確定Kriging模型相關(guān)系數(shù)，建立精確的代理模型，提高結(jié)構(gòu)響應(yīng)的計算效率。最后，利用三自由度質(zhì)量-彈簧系統(tǒng)和三維桁架結(jié)構(gòu)兩個數(shù)值算例驗證所提方法的有效性。

2 基于應(yīng)變模態(tài)的貝葉斯方法

2.1 應(yīng)變模態(tài)

應(yīng)變是由位移經(jīng)過一次求導(dǎo)得來的，對于結(jié)構(gòu)整體而言，應(yīng)變與位移之間的關(guān)系可表示為[13]

ε=BTφs

(1)

式中ε為結(jié)構(gòu)所有點的應(yīng)變,B為結(jié)構(gòu)整體幾何矩陣,T為局部坐標(biāo)和總體坐標(biāo)之間的轉(zhuǎn)換矩陣,φs=φYφTFej ω t為總體坐標(biāo)的節(jié)點位移向量。

根據(jù)有限元分析模型，由模態(tài)疊加法可得

ε=BTφYφTFejω t=ΨεYφTFejω t=

(2)

2.2 貝葉斯方法

基于貝葉斯統(tǒng)計理論實現(xiàn)有限元模型修正，依據(jù)參數(shù)θ的先驗信息結(jié)合觀測數(shù)據(jù)D(模態(tài)信息)的似然函數(shù)來估計參數(shù)的后驗概率分布。其簡化公式可表示為[5]

P(θ|D)=a·P(D|θ)P(θ)

(3)

假設(shè)模型誤差向量為e，則在模型修正中，對于待修正參數(shù)向量θ得到以下關(guān)系，

(4，5)

(6，7)

P(D|θ)=P(ω|θ)·P(Ψε|θ)

(8)

(9)

(10)

為減少可行解主觀成分，在貝葉斯方法中引入熵。對模型修正這個反問題，應(yīng)充分利用最大熵原理作為先驗信息和約束條件。引入熵H表達(dá)式[7]

(11)

(12)

(13)

3 Kriging代理模型的構(gòu)造

3.1 Kriging模型

Kriging模型由一個隨機過程和一個線性回歸模型構(gòu)成，其模型可表示為[14]

f(xi)=fT(xi)β+z(xi)

(14)

式中β=[β1β2…βp]T是回歸模型系數(shù)向量,f(xi)=[f1(xi)f2(xi)…fp(xi)]T是變量x的多項式函數(shù),z(x)是協(xié)方差非零的隨機分布，服從正態(tài)分布N(0,σ2)。利用最小二乘估計可得β和σ2的估計值，分別為

(15)

(16)

式中F為樣本點向量矩陣，Y為樣本點響應(yīng)列向量，R為空間矩陣，其中元素Ri j=R(xi,xj) (i,j=1,2,…,d)，d為樣本數(shù)。β和σ2均為θk的函數(shù)，則在Kriging模型中唯一決定預(yù)測響應(yīng)精度的未知數(shù)是θk。

3.2 BA算法

BA算法從蝙蝠的行為和回聲定位能力出發(fā)，利用脈沖發(fā)射率r和脈沖強度A的多樣性來探測、捕食獵物和尋找棲息地，與其他算法相比，BA在準(zhǔn)確性和有效性方面遠(yuǎn)優(yōu)于其他算法，且沒有許多參數(shù)要進(jìn)行調(diào)整。因此，采用蝙蝠算法對其進(jìn)行尋優(yōu)以提高模型精度。式(17～19)給出了BA算法的主要方程[15]為

fi=fmin+τ(fmax-fmin)

(17)

(18)

(19)

(20)

(21)

3.3 Kriging模型的構(gòu)造

在建立Kriging模型時，需要先設(shè)定相關(guān)系數(shù)θk，其直接影響Kriging模型的預(yù)測精度。采用拉丁超立方抽樣方法抽取初始待修正參數(shù)區(qū)間內(nèi)的樣本，將其分為訓(xùn)練集和測試集(訓(xùn)練集和測試集數(shù)據(jù)互不相同)，并利用有限元方法計算所對應(yīng)的響應(yīng)值。然后，以訓(xùn)練集作為Kriging模型的輸入，對應(yīng)的結(jié)構(gòu)應(yīng)變模態(tài)頻率和振型作為響應(yīng)來構(gòu)建Kriging模型；以測試集對應(yīng)的試驗響應(yīng)值與Kriging模型預(yù)測響應(yīng)值之間的均方根誤差RMSE作為BA的目標(biāo)函數(shù)即式(22)，尋求最優(yōu)相關(guān)系數(shù)。最后，用尋得的最優(yōu)相關(guān)系數(shù)構(gòu)建Kriging模型。

(22)

4 改進(jìn)的MCMC算法

4.1 MH算法

MH算法是一種基于模擬的MCMC技術(shù)，它的一個重要應(yīng)用是從給定概率分布中抽樣，使馬爾可夫鏈穩(wěn)態(tài)于該概率密度。即構(gòu)造一個Markov鏈，通過抽樣獲取后驗分布P(θ|D)的樣本數(shù)據(jù)，根據(jù)大數(shù)定理可獲得修正后參數(shù)θ的均值和方差。在模型修正中MH算法的步驟如下[16,17]。

(1) 選擇具有物理意義的初始值θ0。

(2) 利用當(dāng)前θt的值，依據(jù)建議分布q(·|θt)產(chǎn)生一個候選值θ′。

(3) 根據(jù)接受概率函數(shù)[16](23)計算接受率

(23)

(4) 從均勻分布u(0,1)中產(chǎn)生隨機變量u。

(5) 當(dāng)α(θt,θ′)>u時接受θ′，使θt +1=θ′，否則取θt +1=θt。

(6) 重復(fù)步驟(2～5)，直至收斂。

4.2 基于FPA的MCMC方法

目前，MCMC方法已成為處理復(fù)雜高維積分運算的有效工具，但隨著參數(shù)維度增加，其后驗分布更復(fù)雜，標(biāo)準(zhǔn)MH算法的樣本接受率會隨之降低，致使采樣出現(xiàn)停滯現(xiàn)象。

FPA具有異花授粉(全局尋優(yōu))和自花授粉(局部尋優(yōu))兩大特性。異花授粉過程可以逃離任何局部景觀，進(jìn)而探索更大的空間；自花授粉過程使相似的解更頻繁地選擇，從而保證更快地收斂，其收斂速度本質(zhì)上是指數(shù)的，因此尋優(yōu)效率更高[18]。本文采用FPA改進(jìn)的MCMC抽樣來提高采樣接受率及樣本多樣性。即在N條滿足條件的馬爾可夫鏈中,采用FPA的異花授粉和自花授粉再一次尋優(yōu)來更新候選參數(shù)。用待修正參數(shù)θ的目標(biāo)函數(shù)J值來度量個體所處位置的優(yōu)劣，將個體的優(yōu)勝劣汰過程類比為搜索和優(yōu)化過程中用較好的可行解取代較差可行解的迭代過程，因此，MCMC 算法的接受率有了很大程度提高的同時增加了樣本多樣性，提高了尋優(yōu)效率，且保證樣本收斂在合理的范圍。

(24)

(25)

(26)

5 模型修正過程

根據(jù)上述理論，結(jié)合先驗信息和已構(gòu)造的Kriging代理模型產(chǎn)生的預(yù)測響應(yīng)值，由貝葉斯統(tǒng)計理論推導(dǎo)后驗概率密度函數(shù)。首先，用MCMC算法迭代求解N條馬爾可夫鏈和每次迭代所得的最優(yōu)值；然后，用FPA-MCMC算法迭代求解得到l條馬爾可夫鏈，即可得到待修正參數(shù)的后驗概率密度；最后，由后驗分布統(tǒng)計特征獲得修正后的參數(shù)值，完成模型修正，修正流程如圖1所示。

6 數(shù)值算例

6.1 三自由度彈簧-質(zhì)量系統(tǒng)

使用如圖2所示的三自由度彈簧-質(zhì)量系統(tǒng)對所提算法的有效性進(jìn)行驗證。該系統(tǒng)的確定性參數(shù)為m1=m2=m3=1.0 kg，k3=k4=10 N/m，k6=30 N/m。假設(shè)不確定參數(shù)的試驗真實均值為k1=10 N/m，k2=14 N/m和k5=12 N/m，且服從獨立正態(tài)分布。結(jié)構(gòu)的固有頻率為ω1=3.351 Hz，ω2=6.759 Hz，ω3=9.005 Hz，假定不確定參數(shù)的初始值是k1=6 N/m，k2=16 N/m，k5=10 N/m。選用k1，k2和k5為待修正參數(shù)，表示為θ={θ1,θ2,θ3}。使用貝葉斯公式計算其后驗概率分布，用所提算法對該分布進(jìn)行采樣，迭代次數(shù)n=5000，待修正參數(shù)θi的取值范圍為[5,20],初始值θ0={6,16,10}。

圖1 模型修正流程

圖2 三自由度彈簧-質(zhì)量系統(tǒng)

分別采用標(biāo)準(zhǔn)MCMC算法與FPA-MCMC算法進(jìn)行5000次迭代后得到待修正參數(shù)的馬爾可夫鏈如圖3所示。可以看出，兩種算法得到的馬爾科夫鏈均能收斂到試驗真實值，但改進(jìn)的FPA-MCMC 算法能夠有效克服采樣停滯的缺點，其抽樣多樣性優(yōu)于標(biāo)準(zhǔn)MCMC算法。由式(23)計算得到改進(jìn)算法的樣本接受率為74.5%，與同維度下標(biāo)準(zhǔn)MCMC算法3%的接受率相比，改進(jìn)算法的采樣效果更好。其中，兩種算法迭代所得θ1的馬爾可夫鏈如圖4所示，可以看出，所提算法遍歷性遠(yuǎn)優(yōu)于標(biāo)準(zhǔn)MCMC算法；圖5為θ1的概率分布直方圖，可以看出，其后驗均值基本上與預(yù)設(shè)的真實值相吻合。

雖然初始值的選取不會影響收斂的最終結(jié)果，但對參數(shù)的選取誤差越大，則燃燒期越長，因此去掉樣本集的前300次迭代值，降低所選初始值不佳的影響。然后，檢驗后驗樣本的正態(tài)概率結(jié)果，如 圖6 所示，可以看出，改進(jìn)算法抽樣得到的樣本點基本與假設(shè)的參數(shù)分布線重合，因此樣本的正態(tài)分布假設(shè)成立。并對不確定性參數(shù)θi的置信區(qū)間做出了概率性的估計，如圖7所示，在置信橢圓內(nèi)為含95%樣本參數(shù)的區(qū)域。改進(jìn)算法經(jīng)過5000次迭代后，待修正參數(shù)基本位于置信橢圓輪廓內(nèi)。

圖3 三條馬爾科夫鏈

圖4 θ1兩種算法的馬爾科夫鏈

由后驗分布統(tǒng)計特征得到修正后的不確定性參數(shù)值和響應(yīng)值，修正結(jié)果列入表1和表2，可知待修正參數(shù)在修正后均接近假設(shè)的真實值，且修正后參數(shù)對應(yīng)響應(yīng)值的相對誤差都小于0.02%，表明所提算法修正效果較好。

圖5 θ1概率分布直方圖

圖6 參數(shù)的正態(tài)概率檢驗圖

圖7 置信橢圓的散點圖

表1 修正參數(shù)的對比

表2 修正前后頻率對比

6.2 三維桁架結(jié)構(gòu)

圖8 三維桁架模型

采用拉丁超立方抽樣抽取2500組樣本參數(shù)，分為2000組訓(xùn)練集和500組測試集，選取前6階應(yīng)變模態(tài)頻率值和第60桿單元的前3階應(yīng)變模態(tài)振型值作為目標(biāo)響應(yīng)，通過訓(xùn)練集的樣本和其響應(yīng)值來構(gòu)造Kriging模型。設(shè)置BA的種群大小為20，迭代次數(shù)為100，從區(qū)間[0.1,50]內(nèi)尋得最優(yōu)相關(guān)系數(shù)θk，其中第5階頻率的迭代過程如圖9所示。

圖9 第5階頻率迭代曲線

由優(yōu)化所得θk值構(gòu)建Kriging模型，并重新在參數(shù)取值范圍內(nèi)抽取500組樣本，根據(jù)式(22)計算得到RMSE為1.808×10-7，其中第5階頻率的擬合情況如圖10所示?？梢钥闯?，Kriging模型預(yù)測值和真實值幾乎全部重合，說明所建立的Kriging模型預(yù)測精度很高，可以替代有限元模型進(jìn)行模型修正。

圖10 模型預(yù)測

按假設(shè)試驗?zāi)Ｐ碗S機產(chǎn)生50組試驗樣本，將樣本代入結(jié)構(gòu)得到應(yīng)變模態(tài)作為仿真試驗響應(yīng)，然后根據(jù)試驗響應(yīng)的均值和方差修正有限元模型，兩種算法經(jīng)20000次采樣產(chǎn)生的馬爾可夫鏈如圖11所示?？梢钥闯觯倪M(jìn)算法產(chǎn)生的馬爾可夫鏈一直在試驗?zāi)Ｐ蛥?shù)值附近變化，其接受率與遍歷性均有所提高，改進(jìn)算法的采樣效果較好。為進(jìn)一步驗證所提方法的有效性，假定頻率和模態(tài)振型的噪聲在不同模態(tài)階數(shù)上相互獨立，對試驗頻率和模態(tài)振型加入5%的隨機噪聲來考慮試驗數(shù)據(jù)的不確定性，得到如圖12所示的馬爾科夫鏈，可以看出，待修正參數(shù)均在均值附近收斂，對噪聲具有良好的魯棒性。

圖11 兩種算法的馬爾科夫鏈

圖12 FPA-MCMC算法馬氏鏈

修正后參數(shù)的后驗正態(tài)概率結(jié)果如圖13所示?？梢钥闯?，改進(jìn)算法抽樣所得樣本點與試驗?zāi)Ｐ蛥?shù)值分布線基本重合，能夠較好地模擬后驗樣本。并對修正參數(shù)θi的置信區(qū)間做出了概率性的估計如圖14所示，在置信橢圓內(nèi)為含95%樣本參數(shù)的區(qū)域，可以看出，待修正參數(shù)基本位于置信橢圓輪廓內(nèi)的高概率區(qū)域。修正后參數(shù)的正態(tài)概率分布直方圖如圖15所示，可以看出，修正后參數(shù)的后驗均值基本上與預(yù)設(shè)的真實值相吻合，表明后驗分布的樣本符合正態(tài)分布。

表3為參數(shù)θi的修正結(jié)果，可知當(dāng)初始值選取誤差較大時,仍可得到較精確的結(jié)果。表4 和表5為修正所得應(yīng)變模態(tài)振型和頻率的結(jié)果，可知其修正值和對應(yīng)的響應(yīng)值的相對誤差均較小，因此所提算法的修正結(jié)果精度較高，模型修正效果較好。

圖13 正態(tài)概率檢驗圖

圖14 置信橢圓的散點圖

圖15 參數(shù)的概率分布直方圖

表3 修正參數(shù)對比

表4 修正前后模態(tài)振型對比

表5 修正前后頻率對比

7 結(jié) 論

(1) 在貝葉斯統(tǒng)計理論的基礎(chǔ)上，將FPA的局部優(yōu)化和全局優(yōu)化特性融入MCMC算法，保證了候選樣本的合理性，也提高了尋優(yōu)效率，且對噪聲也具有一定的魯棒性，實現(xiàn)了較高維待修正參數(shù)下的隨機有限元模型修正。

(2) 在構(gòu)造Kriging模型時，利用BA對Kriging模型的相關(guān)系數(shù)進(jìn)行尋優(yōu)，使所構(gòu)造的Kriging模型具有良好的擬合精度和預(yù)測能力，能代替有限元模型進(jìn)行迭代計算，使修正結(jié)果更可靠且提高了模型修正計算效率。

(3) 修正結(jié)果顯示，所提修正方法增加了馬爾可夫鏈的遍歷性和樣本合理性，同時大大地提高了樣本接受率，克服了傳統(tǒng)MH算法在參數(shù)維度高的情況下接受率低的缺點。