基于梯度投影取樣的改進Kriging方法

趙維濤，唐 愷，祁武超

(沈陽航空航天大學 遼寧省飛行器復合材料結(jié)構(gòu)分析與仿真重點實驗室，沈陽 110136)

?

基于梯度投影取樣的改進Kriging方法

趙維濤，唐 愷，祁武超

(沈陽航空航天大學 遼寧省飛行器復合材料結(jié)構(gòu)分析與仿真重點實驗室，沈陽 110136)

在結(jié)構(gòu)可靠性分析過程中，Kriging方法是工程實際中常用方法之一。但樣本點選取對Kriging模型預測精度具有一定影響，樣本點過多計算效率低下，樣本點過少難以保證預測精度?；诮Y(jié)構(gòu)可靠性分析特點，充分利用梯度投影技術(shù)的優(yōu)勢，結(jié)合樣本點篩選技術(shù)，合理選取用于構(gòu)建Kriging模型的樣本點。數(shù)值算例表明，方法求得的設計點、可靠性指標和失效概率均具有滿意的計算精度，且方法所需樣本點數(shù)量明顯少于其它方法，計算效率較高。

結(jié)構(gòu)可靠性；Kriging模型；梯度投影；樣本篩選

在實際工程中，對于某些大型復雜結(jié)構(gòu)，輸入變量和輸出變量之間往往具有較高的非線性關系，很難得到輸入變量和輸出變量之間的明確表達式，這給結(jié)構(gòu)可靠性分析帶來一定的困難。雖然Monte Carlo Simulation(MSC)可以解決此類問題并能夠得到較精準的解，但該方法需要調(diào)用有限元次數(shù)過多，計算效率低下。針對上述問題，國內(nèi)外學者提出了很多改進方法，例如Response Surface Method(RSM)[1-2]和Kriging[3-5]方法等。RSM通常使用一個簡單的函數(shù)代替設計點附近的極限狀態(tài)函數(shù)，但RSM的計算精度對響應函數(shù)形式和樣本點選取具有較大的敏感性[6-8]。Kriging方法相對而言比較靈活，在解決具有高度非線性極限狀態(tài)函數(shù)的可靠性分析問題時，亦能取得良好的擬合效果。Kaymaz[5]指出在結(jié)構(gòu)可靠性分析過程中，如果Kriging模型建立適當，其結(jié)果要優(yōu)于經(jīng)典RSM。張崎和李興斯[9]使用全局網(wǎng)格型和拉丁超立方抽樣對結(jié)構(gòu)進行了基于Kriging模型的可靠性分析，這種全局模型可以將所有樣本點的響應值都模擬出來，相當于得到了極限狀態(tài)函數(shù)的顯式表達式。劉瞻等[10]提出了基于重要抽樣的Kriging方法，并將該方法用于結(jié)構(gòu)的可靠性分析之中。

Kriging模型是由一個參數(shù)模型和一個非參數(shù)隨機過程組成的半?yún)?shù)模型，已有研究[11-14]表明“樣本點的選擇對Kriging模型預測精度具有一定的影響，樣本點過多計算效率低下，樣本點過少難以保證預測精度”。Gomes等[15]指出在結(jié)構(gòu)可靠度計算時，梯度投影技術(shù)是最優(yōu)的樣本點選取方法之一。本文將結(jié)合結(jié)構(gòu)可靠性分析的特點，充分利用梯度投影技術(shù)的優(yōu)勢，合理選取用于構(gòu)建Kriging模型的樣本點，在保證計算精度的前提下，提高結(jié)構(gòu)可靠性分析的計算效率。

1 經(jīng)典Kriging方法簡介

已知有m個n維樣本點[X1,X2,…,Xm]及其響應值Y={y1,y2,…,m}T，其中Xi={xi1,xi2,…,xin}T，i=1,2,…,m。Kriging模型[11-12]如式(1)所示。

(1)

式中：F為回歸函數(shù)矩陣，β為回歸系數(shù)，ε為隨機誤差，下標p為回歸函數(shù)的項數(shù)。

回歸函數(shù)提供模擬的全局近似，即響應值的期望。在經(jīng)典Kriging模型中，常用的回歸函數(shù)為多項式函數(shù)。通常情況下，回歸函數(shù)的形式對擬合精度不起決定性的作用。當回歸函數(shù)取為線性多項式時有p=n+1，此時回歸函數(shù)矩陣F的任意一行，如式(2)所示。

(2)

隨機誤差提供模擬局部偏差的近似，即響應值的局部變化。在經(jīng)典Kriging模型中，假定隨機誤差ε1,ε2,…,εm均服從正態(tài)分布且相互不獨立[13]，即如式(3)和(4)所示。

(3)

(4)

式中：σε為隨機誤差的標準差，R(θ,Xi,Xj)為樣本點Xi和Xj之間的相關系數(shù)，θ={θ1,θ2,…,θn}T為模型參數(shù)。

在經(jīng)典Kriging模型中，相關系數(shù)一般選為高斯模型[14]，即如式(5)所示。

R(θ,Xi,Xj)=exp[-θ1(xi1-xj1)2-θ2(xi2-xj2)2-…-θn(xin-xjn)2]

(5)

在經(jīng)典Kriging模型中，式(1)中的回歸系數(shù)可以利用廣義最小二乘法得到，如式(6)所示。

β*=(FTR-1F)-1FTR-1Y

(6)

式中：R是隨機誤差的相關系數(shù)矩陣，具體形式如式(7)所示。

(7)

(8)

模型參數(shù)θ可通過如下的非線性無約束優(yōu)化問題求得，具體公式如式(9)所示。

(9)

(10)

式中：r(X*)是預測點X*與樣本點之間的相關系數(shù)向量，即如式(11)所示。

r(X*)={R(θ,X1,X*),R(θ,X2,X*),…,R(θ,Xm,X*)}T

(11)

根據(jù)式(10)可求得Kriging模型對預測點X*偏導數(shù)，如式(12)所示。

(12)

式中：Jf(X*)和Jr(X*)分別是f(X*)和r(X*)的雅可比矩陣。

2 改進Kriging方法

2.1 基本流程

本文方法的總體思路如下：首先利用少量樣本點構(gòu)建初始Kriging模型，并利用初始Kriging模型采用FORM(First Order Reliability Method)求解設計點；然后在設計點附近區(qū)域利用梯度投影技術(shù)增加新的樣本點，并對新增樣本點進行篩選；最后將篩選后的樣本點增加到之前樣本點集合中，進而更新Kriging模型以提高Kriging模型在設計點附近的預測精度。本文方法的基本流程如圖1，具體步驟如下:

圖1 計算流程圖

(1) 以均值點為中心，按照Bucher提出的實驗設計[1]確定初始樣本點，并計算樣本點對應的結(jié)構(gòu)響應值；

(2)構(gòu)建初始Kriging模型，其中回歸模型根據(jù)樣本點數(shù)量選擇線性或二次多項式模型，相關系數(shù)模型選取高斯模型；

(3)基于初始Kriging模型，使用FORM求解結(jié)構(gòu)可靠指標β及設計點XD；

(4)在XD處利用梯度投影法增加2n+1個樣本點；

(5)對新增樣本點進行篩選，然后將篩選后的樣本點增加到先前樣本點集合中，并更新已有Kriging模型；

(6)根據(jù)更新后的Kriging模型，使用FORM求解可靠指標β及設計點XD；

(7)重復步驟4到步驟6，直至滿足收斂準則(收斂準則為相鄰兩次迭代的可靠性指標相對變化小于0.1%)；

(8)根據(jù)最終的Kriging模型，利用重要抽樣法計算結(jié)構(gòu)失效概率Pf。

2.2 樣本點選取

設計點XD={xD1,xD2,…,xDn}T，則設計點XD中第j個分量xDj投影用到的向量Tj及單位向量tj分別如式(13)和(14)所示。

Tj=ej-(αTej)α, j=1,2,...,n

(13)

(14)

式中：ej為第j坐標軸方向的單位基向量，α為極限狀態(tài)函數(shù)G(X)在XD處的單位負梯度向量，即如式(15)所示。

(15)

為了使新增樣本點合理地位于真實極限狀態(tài)面附近，用到單位投影向量qj，即如式(16)所示。

(16)

式中：γ為一小數(shù)(γ=10-3)，ω是加權(quán)因子，0≤ω≤1，本文采用全向量投影即ω=1。 根據(jù)上式，新增樣本點Sj坐標如式(17)所示。

(17)

j=1,2,…,n

(18)

式中：σj是第j個隨機變量的標準差，參數(shù)k為取樣步長，文中k=1。

梯度投影法可以得到2n個樣本點，加上設計點XD共新增樣本點個數(shù)為2n+1個。令Xd=Sn+1，新增樣本點與之前樣本點總個數(shù)為m+2n+1，如式(19)所示。

[X1,X2,…,Xm,S1,S2,…,S2n+1]=

(19)

2.3 樣本點篩選

若新增樣本點與之前樣本點之間距離過小，一方面會造成樣本點的浪費，一方面可能會產(chǎn)生病態(tài)的回歸函數(shù)矩陣。因此為提高計算效率和模型預測精度，本文采用如下方法對新增樣本點進行篩選。

首先將新增樣本點與之前樣本點均轉(zhuǎn)換到標準正態(tài)空間內(nèi)，如式(20)和(21)所示。

(20)

(21)

式中：Ui和Hj為標準正態(tài)空間內(nèi)的樣本點，μ和σ分別為隨機變量的均值矢量和標準差矢量。

然后計算任意新增樣本點Hj與之前所有樣本點在標準正態(tài)空間內(nèi)的距離，并判斷該距離是否過小，具體公式如式(22)所示。

‖Ui-Hj‖≤γ i=1,2,…,m，j=1,2,…,2n+1

(22)

式中：γ為一小數(shù)(γ=10-3)

若Hj滿足式(22)，則認為樣本點Sj與之前樣本點之間的距離過小，進而將樣本點Sj刪除；若不滿足，則保留樣本點Sj。

3 算例

3.1 算例1

已知結(jié)構(gòu)極限狀態(tài)函數(shù)，如式(23)所示：

G(X)=x1x2-1500

(23)

式中：x1和x2均服從正態(tài)分布，均值為μ={38，54}T，標準差為σ={3.8,2.7}T。

在標準正態(tài)空間內(nèi)(U空間)，應用本文方法所布置的樣本點如圖2。從圖2中可以看出，本文方法所布置的樣本點能夠較好地落在設計點和極限狀態(tài)曲面附近，能夠提高Kriging模型在設計點附近的預測精度。

圖2 本文方法的樣本點

數(shù)值計算結(jié)果見表1。表1中xD1和xD2為設計點坐標，NFE代表極限狀態(tài)函數(shù)的評估次數(shù)(即樣本點個數(shù))。從表1中可以看出，對比設計點坐標、可靠性指標和失效概率，本文方法均具有較好的計算精度；對比NFE(即樣本點個數(shù))，本文方法所需樣本點數(shù)量明顯少于其它方法，本文方法具有較高的計算效率。

表1 數(shù)值計算結(jié)果(算例1)

方法xD1xD2βNFEPf/10-3MCS---10663經(jīng)典RSM[10]2915350049239246583經(jīng)典Kriging[10]2924751092247513667文獻[10]2921451344247953665FORM29144514682512030-本文方法2914451468251201563

3.2 算例2

已知結(jié)構(gòu)極限狀態(tài)函數(shù)，如式(24)所示：

(24)

式中：x1至x4均服從正態(tài)分布，均值為μ={107,2.5×10-4,0.98,20}T，標準差為σ={106,2.5×10-5,0.098,8}T。

應用本文方法求得的數(shù)值計算結(jié)果見表2。從表2中可以看出，就設計點、可靠性指標和失效概率而言，本文方法均具有滿意的計算精度。本文方法共迭代4次，樣本點總數(shù)量9+9+2+1=21個。在第1次迭代過程中，本文方法共布置9個初始樣本點，在后續(xù)每次迭代過程中均利用梯度投影法新增9個樣本點并對新增樣本點進行篩選，在第3次和第4次迭代過程中分別篩選掉7個樣本點和8個樣本點。由于本文采用了樣本點的篩選技術(shù)，使得樣本點數(shù)量進一步減少。對比NFE，本文方法具有較高的計算效率。

表2 數(shù)值計算結(jié)果(算例2)

方法NFExD1/106xD2/10-4xD3xD4βNFEPf/10-3MCS-----1063328經(jīng)典RSM[10]865273862913481153248053671622957經(jīng)典Kriging[10]115340276830959721166204016723440文獻[10]105310238620964323364004291483339FORM9968525000097382342650434145-本文方法99731250780973923412104330213315

4 結(jié)論

基于結(jié)構(gòu)可靠性分析的特點，本文充分利用梯度投影技術(shù)的優(yōu)勢，合理選取用于構(gòu)建Kriging模型的樣本點，能夠以少量樣本點給出在設計點附近的極限狀態(tài)函數(shù)信息。在迭代計算過程中，為了避免樣本點的距離過近而導致的樣本點浪費現(xiàn)象發(fā)生，本文對新增樣本點進行篩選，進一步合理減少樣本點數(shù)量。

數(shù)值算例分析表明，本文方法求得的設計點、可靠性指標和失效概率均具有滿意的計算精度；相對極限狀態(tài)函數(shù)的評估次數(shù)(即樣本點個數(shù))，本文方法所需樣本點明顯少于其它方法，計算效率較高。

[1]BUCHER C G,BOURGUND U.A fast and efficient response surface approach for structural reliability problems[J].Structural Safety,1990,7(1):57-66.

[2]KIM S H,NA S W.Response surface method using vector projected sampling points[J].Structural Safety,1997,19(1):3-19.

[3]ANTHONY A G,WATSON L T.A comparison of approximation modeling techniques:polynomial versus interpolating models[R].AIAA 98-4758,St.Louis,USA,September ,1998.

[4]TIMOTHY W S.Comparison of response surface and kriging models in the multidisciplinary design of an aerospike nozzle[R].NASA/CR-1998-206935,Hampton,USA,February,1998.

[5]KAYMAZ I.Application of Kriging method to structural reliability problems[J].Structural Safety,2005,27(2):133-151.

[6]WEITAO Z,XUEYAN S,KAI T.A response surface method based on sub-region of interest for structural reliability analysis[J].Structural Engineering and Mechanics,2016,57(4):587-602.

[7]趙維濤,邱志平.基于合理子域的改進響應面方法[J].力學學報,2014,46(3):409-416.

[8]趙維濤,邱志平.基于切平面布點的一種改進響應面方法[J].工程力學,2014,31(10):21-26.

[9]張崎,李興斯.基于Kriging模型的結(jié)構(gòu)可靠性分析[J].計算力學學報,2006,23(2):175-179.

[10]劉瞻,張建國,王燦燦,等.基于優(yōu)化Kriging模型和重要抽樣法的結(jié)構(gòu)可靠度混合算法[J].航空學報,2013,34(6):1347-1355.

[11]SOREN N L,NIELSEN H B,SONDERGAARD J.Aspects of the matlab toolbox DACE[R].IMM-REP-2002-13,Technical University of Denmark,Lyngby,Denmark,2002.

[12]SOREN N L,NIELSEN H B.A matlab Kriging toolbox[R].IMM-TR-2002-12,Technical University of Denmark,Lyngby,Denmark,2002.

[13]RASMUSSEN C E,WILLIAMS CKI.Gaussian processes for machine learning[M].MIT Press,Massachusetts,USA,2006.

[14]AU S K,BECK J L.Estimation of small failure probabilities in high dimensions by subset simulation[J].Probabilistic Engineering Mechanics,2001,16(4):263-277.

[15]GOMES H M,AWRUCH A M.Comparison of response surface and neural network with other methods for structural reliability analysis[J].Structural Safety,2004,26(1):49-67.

(責任編輯：劉劃 英文審校：趙歡)

Improved Kriging method based on gradient projection sampling

ZHAO Wei-tao,TANG Kai,QI Wu-chao

(Key Laboratory of Liaoning Province for Composite Structural Analysis of Aerocraft and Simulation, Shenyang Aerospace University,Shenyang 110136,China)

Kriging method is one of the common methods for analyzing structural reliability in engineering practice.However,the selection of samples has a great influence on the accuracy of Kriging model.A large number of samples can lead to low efficiency,and few samples can not ensure the accuracy of prediction.In this paper,an improved Kriging method was proposed based on the characteristics of structural reliability analysis,the gradient projection technique and the sample screening technique.The method selected appropriate amount of samples to construct the Kriging model.The numerical results indicate that the improved method can obtain approving accuracy for the calculation of design point,reliability index and failure probability.The number of samples selected by the improved method is remarkably less than that by other methods and the calculation efficiency is higher.

structural reliability;Kriging model;gradient projection;sample selecting

2016-05-06

國家自然科學基金(項目編號：11502149)；遼寧省自然科學基金(項目編號：201602579)；航空科學基金(項目編號：2013ZA54004)；遼寧省教育廳(項目編號：L2014072) 。

趙維濤(1977-)，男，遼寧沈陽人，副教授，主要研究方向：飛機結(jié)構(gòu)強度與可靠性，E-mail:zhwt201@126.com。

航空宇航工程

2095-1248(2016)05-0001-05

V215.7

A

10.3969/j.issn.2095-1248.2016.05.001