高斯過程模型的全局靈敏度分析的參數(shù)選擇及采樣方法

萬華平，任偉 新，王寧波

（1.合肥工業(yè)大學土木與水利工程學院，安徽合肥230009；2.中南大學土木工程學院，湖南長沙410004）

高斯過程模型的全局靈敏度分析的參數(shù)選擇及采樣方法

萬華平1，2，任偉 新1，2，王寧波2

（1.合肥工業(yè)大學土木與水利工程學院，安徽合肥230009；2.中南大學土木工程學院，湖南長沙410004）

提出了基于全局靈敏度分析的有限元模型修正參數(shù)選擇方法，考慮了參數(shù)整個變化空間的作用及參數(shù)間的相互作用，具有適用于參數(shù)不確定性大和無模型限制等優(yōu)勢。全局靈敏度分析常采用蒙托卡羅方法計算靈敏度指標，因此足夠大的采樣次數(shù)是獲得可靠靈敏度指標的前提，但是同時會造成計算成本的增加。為此，采用高斯過程模型取代耗時的有限元模型用于降低計算成本，同時探討了拉丁超立方抽樣、Halton序列和Sobol序列3種空間采樣方法用于全局靈敏度分析的計算效率，旨在選擇一種高效的采樣方法。最后，一桁架人行橋?qū)嵗炞C了有限元模型修正參數(shù)選擇和采樣選擇方法。

參數(shù)選擇；有限元模型修正；全局靈敏度分析；高斯過程模型；拉丁超立方抽樣

引 言

在工程領(lǐng)域，有限元模型被廣泛地應用于基于模型的研究工作，一個能夠準確表征結(jié)構(gòu)行為的高精度有限元模型至關(guān)重要，而有限元模型修正是實現(xiàn)高精度有限元模型的方法之一。有限元模型修正包括確定目標函數(shù)、參數(shù)選擇和數(shù)學優(yōu)化估計修正參數(shù)這3個重要步驟［1］。復雜土木工程結(jié)構(gòu)包含的參數(shù)較多，將所有參數(shù)用于模型修正是不現(xiàn)實的。如果不敏感的參數(shù)未被排除在修正參數(shù)之外，修正結(jié)果很可能會失去物理意義，因為自動的數(shù)學優(yōu)化過程將會導致不敏感參數(shù)的大幅度修正。此外，剔除非敏感性參數(shù)會帶來模型降階，從而大大減輕了計算負擔。因此，參數(shù)選擇是結(jié)構(gòu)有限元模型修正的關(guān)鍵問題之一。

目前廣泛用于有限元模型修正參數(shù)篩選的方法為局部靈敏度分析［2］，它采用直接求導法或者有限差分法計算參數(shù)的局部靈敏度。局部靈敏度分析反映的是單個參數(shù)在名義值處的局部梯度信息，不能度量參數(shù)在整個取值范圍的作用和參數(shù)間的共同作用。因此，局部靈敏度分析不適用于參數(shù)變化范圍較大的情況。全局靈敏度分析克服了局部靈敏度分析的不足，而且該方法不受模型限制，也就是說它適用于任何模型，比如高度非線性模型等。全局靈敏度分析已廣泛的應用于各種研究領(lǐng)域，比如結(jié)構(gòu)動力分析［3］，管路系統(tǒng)設(shè)計［4］和結(jié)構(gòu)可靠度分析［5］等。

但是，全局靈敏度分析計算花費高，因為蒙托卡羅方法需要大量采樣以確保靈敏度指標的收斂。土木結(jié)構(gòu)復雜，其通常模擬為高精度有限元模型，這樣使得蠻力蒙托卡羅采樣法（brute-force MCS，即直接基于有限元模型）不切實際。因此，降低全局靈敏度分析應用于復雜土木結(jié)構(gòu)的計算花費顯得尤為重要。為此，本文采用高斯過程模型取代耗時的高精度有限元模型以降低全局靈敏度分析的計算成本。高斯過程模型具有評估預報的不確定性，能很好地處理高維、小樣本和強非線性等復雜問題［6］。近年來，高斯過程模型已經(jīng)在工程領(lǐng)域得到廣泛的應用［5，7-11］。另一方面，選擇一種高效的采樣方法也是降低全局靈敏度分析計算花費的有效途徑?？臻g采樣方法因其具有樣本良好的分布均勻性而受到研究者的青睞。本文對比了拉丁超立方抽樣、Halton序列和Sobol序列3種常用的空間采樣方法在全局靈敏度分析中的計算效率，旨在為選擇一種高效的采樣方法提供參考。

1 高斯過程模型

高斯過程模型是基于貝葉斯思想構(gòu)建的：高斯先驗用來定義模型輸出，通過對訓練樣本集的最大似然估計得到預報值的后驗高斯分布，高斯過程模型完全由它的均值函數(shù)和協(xié)方差函數(shù)完全決定。對于均值函數(shù)一般采用零均值函數(shù)，因為缺乏對潛在函數(shù)的總體趨勢先驗知識［12］，并且零均值函數(shù)可以簡化高斯過程模型的推導；對于協(xié)方差函數(shù)使用平方指數(shù)函數(shù)，該協(xié)方差函數(shù)具有使擬合的函數(shù)光滑且無限可導的優(yōu)勢

假設(shè)有n個觀測值的訓練樣本集D=（X，T），其中X=［XT1，XT2，…，XTn］T，T=［t1，t2，…，tn］T。非觀測點x*的預報值t*亦為高斯，其均值和方差如下

式中 （·）T表示轉(zhuǎn)置運算，C*=［C（X*，X1），C（X*，X2），…，C（X*，Xn）］T，C=C（X，X）和α=C-1T。

高斯過程模型的超參數(shù)Θ通過邊緣釋然函數(shù)最大化求得，即負對數(shù)邊緣釋然函數(shù)最小化

式中L（Θ）為負對數(shù)邊緣釋然函數(shù)，它的表達式及其對超參數(shù)的梯度如下［13］：

式中 |·|和tr（·）分別表示求行列式運算和對矩陣求跡運算。

式（4）的優(yōu)化問題可以通過共軛梯度優(yōu)化算法求解，但容易陷入局部最優(yōu)，為此多初始值策略（Multi-Starting Point Strategy）用于防止優(yōu)化過程的局部最優(yōu)問題［11］。在決定使用構(gòu)建的高斯過程模型代替有限元模型之前，留一法交叉驗證（Leave-One-Out Cross-Validation）用來驗證高斯過程模型的精度。

2 全局靈敏度分析

全局靈敏度分析的基本思想是將模型輸出的總方差分解為由各參數(shù)及其相互作用貢獻的子方差之和。具有d個獨立隨機參數(shù)的模型y（x）可分解為按維數(shù)遞增形式的主效應和交互效應之和［14］

E（y）代表模型y（x） 的均值，zi（xi）表示參數(shù)xi的主效應，zi，j（xi，j）為參數(shù)xi和xj間的相互作用，以此類推可得高階項。

式（7）右邊除E（y）之外的2d-1項互相獨立［14］，因 此 可以 獲 得如 下 表達 式［15］

其中，

靈敏度指標定義于子方差與總方差的比值，如下

評價參數(shù)xi對模型輸出作用大小的指標為一階靈敏度指標Si和總靈敏度指標STi，STi由下式求得

由式（8）～（10）可知，求出Vi，V-i和V便可獲得參數(shù)xi的一階靈敏度指標Si和總靈敏度指標STi。方差項Vi，V-i和V可通過蒙托卡羅方法求?。?/p>

式中X（M1）和X（M2）為兩個獨立的N×d維抽樣矩陣，（X（M2）-im，X（M1）im）表示用數(shù)組X（M1）的第i列向量替代數(shù)組X（M2）的第i列向量后形成的數(shù)組，同理（X（M1）-im，X（M2）im）表示用數(shù)組X（M2）的第i列向量替代數(shù)組X（M1）的第i列向量后形成的數(shù)組。

3 采樣方法

3.1 拉丁超立方抽樣

拉丁超立方抽樣是由Mc Kay等人［16］提出的，它是一種基于分層抽樣思想的方法。拉丁超立方抽樣將變量的概率分布函數(shù)分成N（N為抽樣次數(shù)）個互不重疊的等概率子區(qū)間，然后在每個子區(qū)間內(nèi)分別進行獨立隨機采樣，以確保在每個子區(qū)間精確地采樣一次。拉丁超立方抽樣的樣本可由下式產(chǎn)生

式中π為1到N的獨立隨機排列數(shù)，U為［0，1］之間獨立于π的隨機數(shù)。

3.2 Halton序列

Halton序列［17］是vander Corput序列的擴展。它是通過將一系列整數(shù)表示成素數(shù)基的數(shù)位的形式，然后將這些數(shù)位反序排列再在前面加小數(shù)點得到的值。

設(shè)任意整數(shù)n和任意素數(shù)R，則整數(shù)n均有唯一的R進制表示

式中m=［lnn/lnR］，［·］表示取整運算。然后，將數(shù)位反序排列再在前面加小數(shù)點得到新的值如下

則d維的Halton序列可表示為

式中 （R1，R2，…，Rd）為前d個素數(shù)。

3.3 Sobol序列

Sobol序列［18］是對Halton序列的重新排列，它是基于一組直接數(shù)而構(gòu)造的。Halton序列以一組連續(xù)素數(shù)為基，而Sobol序列均以最小素數(shù)2為基，這樣可以避免Halton方法產(chǎn)生高維序列需要計算大的素數(shù)的麻煩。

任意整數(shù)n以2為基表示如下

式中m=［lnn/lnR］，［·］表示取整運算，n0，…，nm取0或1。

Sobol序列的生成需借助不可約的本原多項式

式中pi（x）為第i維的本原多項式，Si為本原多項式pi（x）的度數(shù)，系數(shù)a1，i，a2，i，…，aSi-1，i取0或1。采 用 遞推 關(guān)系 定義 一 組 正 數(shù) ｛m1，i，m2，i，…，mSi，i｝，如下

式中 ⊕表示二進制按位異或運算，即1⊕0=0⊕1=0，1⊕1=0⊕0=1；mk，i（k=1，2，…，Si）為小于2i的奇數(shù)。

定義直接數(shù)νk，i=（k=1，2，…，Si），則基于直接數(shù)的Sobol序列可由下式求得

Antonov等人［19］提出了一種基于格雷碼（Gray code）快速生成Sobol序列的方法。利用格雷碼的性質(zhì)，Sobol序列可由如下遞歸公式生成［20］

式中 下標ck表示整數(shù)k二進制表示的最右零位。比如，對于0=（0000）2，1=（0001）2，2=（0010）23= （0011）2，可得到c0=1，c1=2，c2=1，c3=3。

4 實橋算例

4.1 橋梁描述及有限元模擬

風荷橋是位于湖南省常德市白馬湖公園內(nèi)的一座鋼桁架人行橋（圖1）。該橋?qū)?.5 m，長115 m （20 m+3×25 m+20 m），橋型布置如圖2所示。該桁架人行橋上部結(jié)構(gòu)由主桁、腹桿、橫縱聯(lián)、橋面系和欄桿扶手。主桁、腹桿、橫縱聯(lián)均由空心鋼管構(gòu)成，其尺寸分別為Φ273 mm×16 mm，Φ140 mm× 12 mm，Φ114 mm×3.5 mm。構(gòu)成橋面系的實木尺寸為3.5 m（長度）×0.2 m（寬度）×0.15 m（厚度），實木之間間隙為0.02 m。欄桿立柱均勻焊接在主桁鋼管上，立柱間間距為2.2 m。

圖1 風荷橋Fig.1 Feheng Bridge

該桁架橋采用ANSYS進行有限元建模，橋面板采用質(zhì)量單元（MASS21），其余構(gòu)件均采用梁單元（BEAM44）。需要指出的是，通常被認作為附屬結(jié)構(gòu)的欄桿扶手在有限元模擬時采用了梁單元而非質(zhì)量單元是因為該橋的欄桿扶手足夠粗壯，其剛度貢獻不應忽略。建立的有限元模型總共有3 035個節(jié)點和4 962個單元（4 832個梁單元，130個質(zhì)量單元），如圖3所示。

圖2 橋型布置圖（單位：mm）Fig.2 Conguration of Fenghe Bridge（Unit：mm）

圖3 有限元模型Fig.3 Three-dimensional finite element model

4.2 參數(shù)選擇及采樣方法比較

該橋?qū)崪y的前五階模態(tài)頻率為4.132，5.096，6.044，6.409，7.248 Hz；對應的初始有限元模型計算頻率為4.352，5.369，6.365，6.721，7.634 Hz。假設(shè)橋梁模態(tài)頻率的降低是由結(jié)構(gòu)的剛度折減而引起。另外鋼橋在實際服役過程中，鋼材彈性模量易受工作環(huán)境影響（比如溫度等）而發(fā)生改變。因此，本文選擇彈性模量作為修正參數(shù)，根據(jù)構(gòu)件所處的位置和作用劃分為6個子塊（如表1所示）。

表1 參數(shù)列表Tab.1 The list of parameters

采用基于頻率殘差的無量綱目標函數(shù)

式中fiGP和fiTest分別表示高斯過程預測頻率和實測頻率。由文獻［21］可知，即使實測頻率和計算頻率相差不大，修正參數(shù)的幅度也可能很大。文獻［21］中彈性模量修正量最大為30%，所以參數(shù)先驗范圍設(shè)定為名義值210 GPa的±30%改變量。

對前面敘述的3種空間采樣方法，均采用足夠大樣本數(shù)5×106計算各參數(shù)的一階靈敏度指標Si和總靈敏度指標STi，結(jié)果如圖4所示（圖中數(shù)字僅給出了有限位精度）。結(jié)果表明：

圖4 參數(shù)對目標函數(shù)的靈敏度指標Fig.4 Sensitivity indices of parameters with respect to objective function

·參數(shù)E2，E3，E6與其他參數(shù)的相互作用明顯，因為它們各自的總靈敏度指標明顯大于一階靈敏度指標；

·參數(shù)E2對目標函數(shù)起著主導作用，對目標函數(shù)的總效應占82.3%，這是因為主桁構(gòu)成了該橋的主骨架，起著重要作用，因此參數(shù)E2應該確定為修正參數(shù)；

·參數(shù)E3，E6對目標函數(shù)影響也比較顯著，所以它們均應該選擇為修正參數(shù)；

·參數(shù)E1，E4，E5應該從修正參數(shù)組移除，因為它們對目標函數(shù)的影響極其微小。

接下來可考察3種空間采樣方法用于計算靈敏度指標的效率。采樣次數(shù)依次為103，5×103，10× 103，15×103，…，100×103，共21種采樣。引入絕對誤差（Absolute Error，AE）評價靈敏度指標計算的收斂效果如下

式中SMaxi（SMaxTi）為采樣次數(shù)為5×106的一階（總）靈敏度指標設(shè)為指標，SSami（SSamTi）表示采樣次數(shù)為103，5×103，…，100×103的一階（總）靈敏度指標，d為參數(shù)個數(shù)。3種空間采樣方法效果對比如圖5所示。

圖5 采樣方法計算收斂效果比較Fig.5 Convergence comparison of sampling methods

由圖5可知，Sobol序列用于全局靈敏度分析計算效率最高，Halton序列次之，最后為拉丁超立方抽樣。它們計算效率的高低可以從其均勻性優(yōu)劣得到解釋：拉丁超立方保證了一維的均勻性，但未考慮維度間的均勻性；Sobol序列和Halton序列為2種低差異性（Low Discrepancy）的擬蒙托卡羅采樣（Quasi-Monte Carlo）方法，它們在維度間呈現(xiàn)出良好的均勻性；Sobol序列避免了Halton序列高維抽樣時的集聚效應，前者的均勻性優(yōu)于后者［22］。

5 結(jié) 論

本文提出了基于全局靈敏度分析的有限元模型修正參數(shù)選擇方法。相比傳統(tǒng)的局部靈敏度分析方法，全局靈敏度分析方法具有量化參數(shù)整個變化范圍作用，考慮了參數(shù)間的相互作用等優(yōu)點。但是，全局靈敏度分析應用于復雜的土木結(jié)構(gòu)計算花費高。為此，一方面，本文采用高斯過程模型用于降低全局靈敏度分析方法計算成本高；另一方，比較了拉丁超立方抽樣、Halton序列、Sobol序列3種空間采樣方法用于全局靈敏度分析的計算效率，對尋找一種高效的采樣方法用于減少全局靈敏度分析計算花費具有指導意義。

［1］ Ren W X，Chen H B.Finite element model updating in structural dynamics by using the response surface method［J］.Engineering Structures，2010，32（8）：2 455—2 465.

［2］ Brownjohn J M W，Xia P Q，Hao H，et al.Civil structure condition assessment by fe model updating：Methodology and case studies［J］.Finite Elements in Analysis and Design，2001，37（10）：761—775.

［3］ 于德介，李睿.Sobol’法在非線性隔振系統(tǒng)靈敏度分析中的應用研究［J］.振動工程學報，2004，17（2）：210—213. Yu D J，Li R.Application of Sobol′method to sensitivity analysis of a non-linear passive vibration isolators ［J］.Journal of Vibration Engineering，2004，17（2）：210—213.

［4］ 陳剛，汪玉，毛為民，等.沖擊載荷作用下艦艇管路系統(tǒng)全局參數(shù)靈敏度分析［J］.振動與沖擊，2007，26 （3）：45—48. Cheng G，Wang Y，Mao Y M，et al.Global parameter sensitivity analysis of shipboard piping systems under shock loads［J］.Journal of Vibration and Shock，2006，26（3）：45—48.

［5］ Wang P，Lu Z，Tang Z.An application of the Kriging method in global sensitivity analysis with parameter uncertainty［J］.Applied Mathematical Modelling，2013，37（9）：6 543—6 555.

［6］ Simpson T W，Poplinski J D，Koch P N，et al.Metamodels for computer-based engineering design：Survey and recommendations［J］.Engineering with Computers，2001，17（2）：129—150.

［7］ 劉開云，劉保國，徐沖.基于遺傳-組合核函數(shù)高斯過程回歸算法的邊坡非線性變形時序分析智能模型［J］.巖石力學與工程學報，2009，28（10）：2 128—2 134. Liu K Y，Liu B G，Xu C.Intelligent analysis model of slope nonlinear displacement time series based on genetic-Gaussian process regression algorithm of combined kernel function［J］.Chinese Journal of Rock Mechanics and Engineering，2009，28（10）：2 128—2 134.

［8］ 劉信恩，肖世富，莫軍.高斯過程響應面法研究［J］.應用力學學報，2010，27（1）：190—195. Liu X E，Xiao S F，Mo J.A new response surface method based on Gaussian process［J］.Chinese Journal of Applied Mechanics，2010，27（1）：190—195.

［9］ 萬華平，任偉新，魏錦輝.基于高斯過程響應面的結(jié)構(gòu)有限元模型修正方法［J］.振動與沖擊，2012，31 （24）：82—87. Wan H P，Ren W X，Wei J H.Strucutral finite element model updating based on Gaussian process response surface methodology［J］.Journal of Vibration and Shock，2012，31（24）：82—87.

［10］張冬冬，郭勤濤.Kriging響應面代理模型在有限元模型確認中的應用［J］.振動與沖擊，2013，32（9）：187—191. Zhang D D，Guo Q T.Application of Kriging response surface in finite element model validation［J］.Journal of Vibration and Shock，2013，32（9）：187—191.

［11］Wan H P，Mao Z，Todd M D，et al.Analytical uncertainty quantification for modal frequencies with structural parameter uncertainty using a Gaussian process metamodel［J］.Engineering Structures，2014，75：577—589.

［12］Neal R M.Regression and classification using Gaussian process priors.In：Bernardo J，Berger J，Dawid A，et al.Bayesian Statistics 6［M］.Oxford University Press，1998，475—501.

［13］Rasmussen C E，Williams C K I.Gaussian Processes for machine learning［M］.MIT Press，2006.

［14］Efron B，Stein C.The jackknife estimate of variance ［J］.The Annals of Statistics，1981：586—596.

［15］Sobol I M.Sensitivity estimates for nonlinear mathematical models［J］.Mathematical Modeling and Computational Experiment，1993，1（4）：407—414.

［16］Mc Kay M D，Beckman R J，Conover W J.A comparison of three methods for selecting values of input variables in the analysis of output from a computer code ［J］.Technometrics，1979，21：239—245.

［17］Halton J H.On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals［J］.Numerische Mathematik，1960，2（1）：84—90.

［18］Sobol I M.On the distribution of points in a cube and the approximate evaluation of integrals［J］.USSR Computational Mathematics and Mathematical Physics，1967，7：86—112.

［19］Antonov I A，Saleev，V M.An economic method of computing LPτ-sequences［J］.USSR Computational Mathematics and Mathematical Physics，1979，19（1），252—256.

［20］Bratley P，F(xiàn)ox B L.Algorithm 659：Implementing Sobol＇s quasirandom sequence generator［J］.ACM Transactions on Mathematical Software（TOMS），1988，14（1），88—100.

［21］Jaishi B，Ren W X.Structural finite element model updating using ambient vibration test results［J］.Journal of Structural Engineering，2005，131（4）：617—628.

［22］Cheng J，Druzdzel M J.Computational investigation of low-discrepancy sequences in simulation algorithms for bayesian networks［A］.Proceedings of the Sixteenth Conference on Uncertainty in Artificial Intelligence ［C］.Morgan Kaufmann Publishers，Stanford，CA，2000，72—81.

WAN Hua-ping1，2，REN Wei-xin1，2，WANG Ning-bo2
（1.School of Civil and Hydraulic Engineering，Hefei University of Technology，Hefei 230009，China 2.School of Civil Engineering，Central South University，Changsha 410004，China）

This paper proposes a global sensitivity analysis（GSA）approach for parameter selection in finite element modal updating.The GSA has the capability to quantify the effects of individual parameters over their entire space and interaction effects among parameters.This approach has a couple of advantages，such as capability to high degree of parameter uncertainty and model independence.However，it is still computationally intensive because a large number of model evaluations for Monte Carlo simulation（MCS）are needed to obtain a confident estimate of the sensitivity indices（SIs）.Therefore，the fast-running Gaussian process model is adopted to replace the time-consuming finite element model to alleviate the computational burden.In addition，this study explores the efficiency of three space-filling sampling methods，i.e.，Latin hypercube sampling，Halton sequence and Sobol sequence，in the SIs evaluation，which has great value for choosing an efficient sampling method.At last，a real-world truss pedestrian bridge is used to demonstrate the proposed framework.

parameter selection；finite element modal updating；global sensitivity analysis；Gaussian process model；Latin hypercube sampling

A Gaussian process model based global sensitivity analysis approach for parameter selection and sampling methods

TB123

A

1004-4523（2015）05-0714-07

10.16385/j.cnki.issn.1004-4523.2015.05.005

萬華平（1986—），男，博士，講師。電話：15556928238；E-mail：huaping.wan@gmail.com

任偉新（1960—），男，博士，教授。電話：13856966457；E-mail：renwx@hfut.edu.cn

2014-04-21

2014-12-15

國家自然科學基金資助項目（51508144）；中央高校基本科研業(yè)務費專項資金資助項目（JZ2015 HGBZ0098，JZ2015HGQC0215）