基于序貫設(shè)計(jì)和高斯過程模型的結(jié)構(gòu)動(dòng)力不確定性量化方法

萬華平，張梓楠，周家偉，任偉新

(1.浙江大學(xué) 建筑工程學(xué)院，浙江 杭州 310058；2.浙江大學(xué)平衡建筑研究中心，浙江 杭州 310028；3.浙江大學(xué)建筑設(shè)計(jì)研究院有限公司，浙江 杭州 310028；4.深圳大學(xué) 土木與交通工程學(xué)院，廣東 深圳 518060)

結(jié)構(gòu)動(dòng)力特性對(duì)結(jié)構(gòu)動(dòng)力設(shè)計(jì)和振動(dòng)分析有重要意義，它反映了結(jié)構(gòu)系統(tǒng)剛度、質(zhì)量分布及邊界條件等信息.由于構(gòu)件制造誤差、老化損傷及自身固有的隨機(jī)性等因素影響，結(jié)構(gòu)參數(shù)具有不確定性，從而導(dǎo)致結(jié)構(gòu)動(dòng)力特性具有隨機(jī)性[1].為了提供更準(zhǔn)確的結(jié)構(gòu)動(dòng)力信息[2-3]，須定量計(jì)算結(jié)構(gòu)參數(shù)傳遞到結(jié)構(gòu)動(dòng)力的不確定性.蒙特卡羅模擬法（Monte Carlo simulation,MCS）是結(jié)構(gòu)動(dòng)力不確定性量化常用的方法，其通過大量的有限元模型計(jì)算得到動(dòng)力響應(yīng)的統(tǒng)計(jì)特性[4]，不足之處在于計(jì)算成本過高.為了克服MCS方法的不足，代理模型方法構(gòu)建用于替代結(jié)構(gòu)有限元模型的近似數(shù)學(xué)模型.不確定性量化的代理模型包括響應(yīng)面[5-6]、多項(xiàng)式混沌展開[7]、高斯過程模型（Gaussian process model,GPM）[8]等.高斯過程模型是一種非參數(shù)概率模型，廣泛用于不確定性量化[9-11]，具有模擬靈活性和預(yù)測(cè)不確定性定量等優(yōu)點(diǎn).

代理模型的預(yù)測(cè)精度受試驗(yàn)設(shè)計(jì)方法的影響.Santner等[12-13]研究不同試驗(yàn)設(shè)計(jì)方法，比較不同的試驗(yàn)設(shè)計(jì)和代理模型的擬合效果.對(duì)于代理模型建模，空間填充設(shè)計(jì)是較好的選擇，包括拉丁超立方設(shè)計(jì)[14-15]、Sobol序列抽樣[16]、Hammersley序列抽樣[17]等，這些能夠較好地均勻填充參數(shù)變化空間.采用傳統(tǒng)試驗(yàn)設(shè)計(jì)構(gòu)建GPM，是事先設(shè)定樣本的數(shù)量，一次生成所有訓(xùn)練樣本，會(huì)因?yàn)檫^度采樣導(dǎo)致計(jì)算資源浪費(fèi).序貫設(shè)計(jì)[18]是把樣本選擇和代理模型結(jié)合起來，基于填充準(zhǔn)則依次選擇最佳樣本點(diǎn)，自適應(yīng)地更新當(dāng)前代理模型，以較少樣本建立一個(gè)高精度的代理模型.

本研究提出基于序貫設(shè)計(jì)和高斯過程模型的結(jié)構(gòu)動(dòng)力不確定性量化方法.基于GPM的預(yù)測(cè)方差（mean squared error,MSE）信息，通過最大化預(yù)測(cè)方差（maximizing mean squared error,MMSE）[19]依次選擇最佳設(shè)計(jì)點(diǎn)，逐步自適應(yīng)地更新GPM.基于建立的自適應(yīng)GPM，動(dòng)力特性統(tǒng)計(jì)矩的復(fù)雜高維積分可以轉(zhuǎn)化為簡(jiǎn)單一維積分，進(jìn)而可以解析計(jì)算出結(jié)構(gòu)動(dòng)力特性統(tǒng)計(jì)矩.

1 高斯過程模型與序貫設(shè)計(jì)

1.1 高斯過程模型

高斯過程模型是基于貝葉斯理論的一種非參數(shù)概率模型，利用高斯過程先驗(yàn)，得到待預(yù)測(cè)點(diǎn)的后驗(yàn)概率分布.GPM完全由其均值函數(shù)m(x)和平方指數(shù)協(xié)方差函數(shù)C(x,x')決定[9-10]，通常均值函數(shù)采用常數(shù)形式μ，平方指數(shù)協(xié)方差函數(shù)計(jì)算式為

式中：η2為協(xié)方差函數(shù)的變化尺度；d為輸入?yún)?shù)的維度；xk為參數(shù)x的第k個(gè)元素；lk為協(xié)方差函數(shù)的變化速率；協(xié)方差函數(shù)的參數(shù) θ 為超參數(shù)，

假設(shè)有n個(gè)觀測(cè)值的樣本點(diǎn)集D={X,Y}，其中根據(jù)高斯過程先驗(yàn)，模型輸出服從高斯分布：

待預(yù)測(cè)點(diǎn)x*處的預(yù)測(cè)值y*與Y也服從高斯分布，計(jì)算式為

根據(jù)貝葉斯原理，預(yù)測(cè)值y*的后驗(yàn)分布為

將式(2)、(3)代入式(4)，得到

式中：μ*為 待預(yù)測(cè)點(diǎn)的均值e為長(zhǎng)度為n的單位列向量；α=C-1(Y-eμ*) ；C*為待預(yù)測(cè)點(diǎn)與訓(xùn)練樣本點(diǎn)之間的協(xié)方差矩陣，表示待預(yù)測(cè)點(diǎn)與訓(xùn)練樣本點(diǎn)的相關(guān)性，C*=[C(x1,x*),C(x2,x*),···,C(xn,x*)]T；C為訓(xùn)練樣本點(diǎn)之間的協(xié)方差矩陣，表示訓(xùn)練樣本點(diǎn)自身的相關(guān)性，C=C(X,X).

1.2 序貫設(shè)計(jì)

序貫設(shè)計(jì)是基于樣本填充準(zhǔn)則，依次選擇最有用的樣本點(diǎn)填充至初始樣本集中，不斷地通過新樣本點(diǎn)更新代理模型.序貫設(shè)計(jì)的關(guān)鍵是如何根據(jù)樣本點(diǎn)和當(dāng)前代理模型選擇新的樣本點(diǎn)，即如何定義樣本填充準(zhǔn)則.GPM能夠提供預(yù)測(cè)方差信息（mean squared error,MSE），MSE反映GPM與原始物理模型之間的差異.MSE越大，表明GPM的預(yù)測(cè)誤差越大.在預(yù)測(cè)方差最大處增加樣本點(diǎn)，能夠減小GPM的預(yù)測(cè)誤差.因此提出基于最大化預(yù)測(cè)方差的樣本填充準(zhǔn)則，計(jì)算式為

樣本迭代采用的停止準(zhǔn)則[20]為

式中：r為比例系數(shù)，取值為0.01%～1.00%；ΔY為輸出響應(yīng)中最大值與最小值的差值.

通過MMSE填充準(zhǔn)則的動(dòng)態(tài)序貫設(shè)計(jì)，構(gòu)建自適應(yīng)GPM，由于包含最優(yōu)設(shè)計(jì)點(diǎn)，可以有效地提高GPM的預(yù)測(cè)精度.建立自適應(yīng)GPM的具體步驟如下.

1) 獲取初始訓(xùn)練樣本點(diǎn)X；

2) 以X作為輸入值計(jì)算有限元模型，得到相應(yīng)的響應(yīng)輸出值Y；

3) 得到初始樣本集D={X,Y}，建立初始GPM；

2 動(dòng)力特性統(tǒng)計(jì)矩的解析計(jì)算

將結(jié)構(gòu)動(dòng)力特性的輸入?yún)?shù)x與輸出響應(yīng)y的關(guān)系用自適應(yīng)GPM來表達(dá)，其中輸入?yún)?shù)x的概率分布函數(shù)為p(x).根據(jù)統(tǒng)計(jì)原理，動(dòng)力特性的均值和方差表達(dá)式為

利用平方指數(shù)協(xié)方差函數(shù)的分離特性，式（6）、（7）可以轉(zhuǎn)換為

當(dāng)參數(shù)xk服從其他分布時(shí)，可根據(jù)概率相等的原則將其轉(zhuǎn)化為服從正態(tài)分布的參數(shù)，采用上述解析結(jié)果.不同概率分布轉(zhuǎn)化的表達(dá)式為

式中：Fxk為參數(shù)xk的累 積分布函 數(shù)；Φ 為標(biāo) 準(zhǔn)正態(tài)累積分布函數(shù).

3 自適應(yīng)高斯過程模型擬合過程展示

采用一維函數(shù)f1(x) 來說明本研究所提自適應(yīng)GPM的擬合過程，表達(dá)式為

初始樣本設(shè)為4個(gè)，根據(jù)MMSE填充準(zhǔn)則，在每次迭代過程中選擇最優(yōu)的設(shè)計(jì)點(diǎn)填充至初始樣本集中，并更新當(dāng)前的GPM，停止準(zhǔn)則的比例系數(shù)r=0.02%.整個(gè)迭代建模的部分動(dòng)態(tài)過程如圖1所示，可以看出，隨著迭代次數(shù)的增加，GPM的預(yù)測(cè)均值與真實(shí)函數(shù)曲線的擬合程度越來越高；GPM的預(yù)測(cè)方差逐漸減小，即圖中95%的置信區(qū)間面積不斷減小，表明所建立的GPM精度越來越高.

圖1 一維函數(shù)自適應(yīng)GPM的迭代過程Fig.1 Evolution of adaptive GPM for a one-dimensional test function

二維函數(shù)f2(x) 用來進(jìn)一步展示自適應(yīng)GPM方法的迭代擬合過程，表達(dá)式為

初始樣本設(shè)為10個(gè)，停止準(zhǔn)則的比例系數(shù)r=0.02%.在每次迭代過程中，通過MMSE選擇最佳設(shè)計(jì)點(diǎn)，迭代過程如圖2所示.圖中，實(shí)線為真實(shí)函數(shù)，虛線為自適應(yīng)GPM預(yù)測(cè)值，圓點(diǎn)為初始設(shè)計(jì)點(diǎn)，菱形為序貫設(shè)計(jì)點(diǎn)，正三角形為最優(yōu)設(shè)計(jì)點(diǎn).對(duì)于二維函數(shù)，以等高線圖呈現(xiàn)擬合程度來提高可視化.與一維函數(shù)類似，隨著迭代次數(shù)增加，GPM的預(yù)測(cè)曲線與真實(shí)函數(shù)曲線之間的擬合優(yōu)度提高.

圖2 二維函數(shù)自適應(yīng)GPM的迭代過程Fig.2 Evolution of adaptive GPM for a two-dimensional test function

為了更清晰地說明基于MMSE的樣本填充機(jī)制，迭代過程的部分預(yù)測(cè)方差如圖3所示.可以看出，隨著迭代次數(shù)的增加，GPM的預(yù)測(cè)方差明顯減小，進(jìn)一步表明基于MMSE填充準(zhǔn)則建立的自適應(yīng)GPM可以有效地提高模型預(yù)測(cè)精度.

圖3 二維函數(shù)自適應(yīng)GPM迭代過程的MSEFig.3 MSE against number of iterations for a two-dimensional test function

4 方法驗(yàn)證

將基于序貫設(shè)計(jì)的自適應(yīng)高斯過程模型方法用于一個(gè)雙層柱面網(wǎng)殼（見圖4）的動(dòng)力特性不確定性量化.該網(wǎng)殼跨度為30 m，矢高為6.8 m，長(zhǎng)度為50 m，采用兩邊支承，整個(gè)網(wǎng)殼均由截面直徑d=80 mm、薄壁厚度為2 mm的圓形薄壁鋼管組成.在ANSYS環(huán)境下建立該網(wǎng)殼的有限元模型，采用LINK8桿單元模擬全部桿件，所建有限元模型共326個(gè)節(jié)點(diǎn)，1 200個(gè)桿單元.該網(wǎng)殼的有限元模型及前4階振型如圖5所示.

圖4 雙層柱面網(wǎng)殼Fig.4 Double layer cylindrical reticulated shell

圖5 雙層柱面網(wǎng)殼有限元模型和前4階振型Fig.5 Finite element model and first-four-order vibration modes of double layer cylindrical reticulated shell

假定鋼管半徑、鋼材密度和彈性模量為不確定性參數(shù)（見表1），以網(wǎng)殼的前4階固有頻率為分析對(duì)象，計(jì)算不確定性參數(shù)下該結(jié)構(gòu)固有頻率的統(tǒng)計(jì)特性.基于MMSE建立自適應(yīng)GPM用于結(jié)構(gòu)固有頻率的統(tǒng)計(jì)矩計(jì)算.與傳統(tǒng)GPM方法對(duì)比，大量次數(shù)（2×104）采樣的MCS方法用于近似固有頻率統(tǒng)計(jì)矩真值.計(jì)算平臺(tái)：電腦品牌為DELL，處理器為Intel(R) Core(TM) i5-10500CPU@3.10 GHz，有限元軟件為ANSYS @19.2.自適應(yīng)GPM法、GPM法和MCS法的計(jì)算結(jié)果如表2～5所示，均采用Sobol序列抽樣獲取初始樣本點(diǎn)，初始樣本數(shù)為12.自適應(yīng)GPM的樣本填充停止準(zhǔn)則的比例系數(shù)r=0.05%.由表2～5可知，在增加了15個(gè)樣本點(diǎn)后，自適應(yīng)GPM已經(jīng)達(dá)到較高的計(jì)算精度，均值最大誤差為0.003 7%，方差最大誤差為0.561 7%，而傳統(tǒng)GPM方法在任意增加32個(gè)樣本點(diǎn)之后才達(dá)到相似的精度，表明所提基于序貫設(shè)計(jì)的自適應(yīng)高斯過程模型方法能夠在保證高精度的同時(shí)顯著降低計(jì)算成本.

表1 不確定性參數(shù)的統(tǒng)計(jì)特征Tab.1 Statistical characteristics of uncertain parameters

表2 自適應(yīng)GPM、GPM和MCS法的均值和方差計(jì)算結(jié)果對(duì)比（第1階固有頻率）Tab.2 Comparison of mean and variance of adaptive GPM,GPM and MCS (first natural frequency)

表3 自適應(yīng)GPM、GPM和MCS法的均值和方差計(jì)算結(jié)果對(duì)比（第2階固有頻率）Tab.3 Comparison of mean and variance adaptive GPM,GPM and MCS (second natural frequency)

表4 自適應(yīng)GPM、GPM和MCS法的均值和方差計(jì)算結(jié)果對(duì)比（第3階固有頻率）Tab.4 Comparison of mean and variance adaptive GPM,GPM and MCS (third natural frequency)

表5 自適應(yīng)GPM、GPM和MCS法的均值和方差計(jì)算結(jié)果對(duì)比（第4階固有頻率）Tab.5 Comparison of mean and variance adaptive GPM,GPM and MCS (fourth natural frequency)

5 結(jié)論

(1) 所提自適應(yīng)GPM方法可以通過較少樣本數(shù)量達(dá)到較高的預(yù)測(cè)精度，均值和方差的相對(duì)誤差均較低，而傳統(tǒng)GPM則需要較多樣本點(diǎn)才能達(dá)到相當(dāng)?shù)挠?jì)算精度，表明所提方法能夠明顯降低計(jì)算成本.

(2) 迭代過程中統(tǒng)計(jì)矩的相對(duì)誤差的變化表明，所提自適應(yīng)GPM法的相對(duì)誤差在幾次迭代后明顯降低，而傳統(tǒng)GPM的相對(duì)誤差在增加大量樣本后仍然較高，表明所提自適應(yīng)GPM法相對(duì)于傳統(tǒng)GPM法具有高效率的優(yōu)勢(shì).

(3) 所提自適應(yīng)高斯過程模型方法只涉及一種樣本填充準(zhǔn)則，可進(jìn)一步研究其他序貫抽樣方法，如最大熵準(zhǔn)則、整體均方誤差準(zhǔn)則，比較幾種抽樣方法的計(jì)算結(jié)果，得到不同抽樣準(zhǔn)則的使用范圍。