Kriging與改進(jìn)一次二階矩融合的可靠性分析方法*

袁修開,孔沖沖,顧 健

(廈門大學(xué) 航空航天學(xué)院， 福建 廈門 361005)

在結(jié)構(gòu)可靠性分析中，失效概率的求解方法可分為三類：近似解析法、數(shù)字模擬法、代理模型方法[1]。近似解析方法包括改進(jìn)一次二階矩(Advanced First Order Second Moment, AFOSM)[1-2]、均值一次二階矩(Mean Value First Order Second Moment, MVFOSM)和R-F法[3]等。近似解析方法計算量少，但是在處理復(fù)雜非線性問題時，其精度難以保證。數(shù)值模擬法包括蒙特卡洛仿真(Monte Carlo Simulation， MCS)[4]、線抽樣(Line Sampling, LS)[5-6]法、重要抽樣(Importance Sampling, IS)[7]和子集模擬(Subset Simulation, SS)[8-9]等。該類方法為了保證計算精度，需調(diào)用結(jié)構(gòu)功能函數(shù)的計算次數(shù)較多，計算代價大，且在處理包含有限元模型的隱式極限狀態(tài)問題時，計算效率低下，是影響該類方法應(yīng)用的重要因素。代理模型方法包括響應(yīng)面法[7-11]、神經(jīng)網(wǎng)絡(luò)法[12-13]、支持向量機(jī)法(Support Vector Machine, SVM)[14-16]等。該類方法專門針對隱式極限狀態(tài)問題，能夠顯著提高分析效率，因而在工程中得到廣泛應(yīng)用，其中Kriging方法作為一種典型的連續(xù)插值迭代方法，以其精確的插值技術(shù)備受人們的關(guān)注[17]。本文著眼于采用該技術(shù)來進(jìn)一步提高改進(jìn)一次二階矩方法的分析效率和適用范圍。

改進(jìn)一次二階矩法通過將非線性功能函數(shù)線性展開，然后用線性功能函數(shù)的失效概率來近似原非線性功能函數(shù)的失效概率[1]。和均值一次二階矩法相比，AFOSM在設(shè)計點(diǎn)(Most Probable Point, MPP)處而非均值點(diǎn)處將功能函數(shù)線性化，提高了分析計算的精度。呂震宙等[1]對改進(jìn)一次二階矩的理論介紹和推導(dǎo)過程做了詳細(xì)的描述。改進(jìn)一次二階矩亦在眾多工程領(lǐng)域得到了廣泛的應(yīng)用，比如：葛耀君等[18]在橋梁顫振可靠性評估中，使用改進(jìn)一次二階矩方法計算了小失效概率條件下的可靠度；侯曉亮等[19]在評價軟土基坑支護(hù)設(shè)計中的抗隆起穩(wěn)定性時，采用改進(jìn)的一次二階矩可靠度計算方法評價基坑抗隆起穩(wěn)定性；鄭財?shù)萚20]在研究三軸數(shù)控機(jī)床的運(yùn)行誤差時，用改進(jìn)一次二階矩法對其可靠性及靈敏度進(jìn)行了分析計算；曾照輝等[21]在研究動力渦輪工作時的可靠性時，采用參數(shù)化建模通過混合模擬(有限元、響應(yīng)面、改進(jìn)一次二階矩法三者結(jié)合)的方法對其進(jìn)行可靠性分析等。改進(jìn)一次二階矩針對變量維數(shù)小，非線性程度不大的小失效概率問題來說，具有很高的計算效率和良好的分析效果。然而在改進(jìn)一次二階矩中求解“設(shè)計點(diǎn)”，及其包含的偏導(dǎo)數(shù)的計算是難點(diǎn)，目前常用的方法為有限差分法，但是其對于隱式極限狀態(tài)函數(shù)(需有限元分析)，有限差分的步長很難確定，且增加了功能函數(shù)的計算次數(shù)，增大了計算量、效率低。

Kriging代理模型不僅具有良好的擬合效果和局部估計的特點(diǎn)，而且具有較好的連續(xù)性和可導(dǎo)性，故具有廣泛的應(yīng)用前景。韓忠華[22]在Kriging模型及代理優(yōu)化算法研究進(jìn)展中，對Kriging方法的背景、意義、理論以及發(fā)展現(xiàn)狀做了詳細(xì)的描述；聶雪媛等[23]在研究飛行器結(jié)構(gòu)剛度氣動優(yōu)化設(shè)計中采用Kriging方法建立代理模型，該方法能夠處理復(fù)雜目標(biāo)的全局優(yōu)化問題；黃曉旭等[24]提出一種將Kriging模型與子集模擬方法結(jié)合的可靠性分析方法用于解決小失效概率的工程結(jié)構(gòu)問題；陳立立等[25]采用Kriging代理模型進(jìn)一步驗證自由變形技術(shù)在RAE2822翼型優(yōu)化設(shè)計中的應(yīng)用；韓少強(qiáng)等[26]提出了一種將梯度信息與Kriging模型構(gòu)建相結(jié)合的方法用于氣動反設(shè)計研究。

基于改進(jìn)一次二階矩的實(shí)用性及Kriging模型的效率和良好的可導(dǎo)性，本文提出一種Kriging與改進(jìn)一次二階矩融合的可靠性分析方法。所提方法在改進(jìn)一次二階矩迭代計算設(shè)計點(diǎn)的過程中借助Kriging模型計算迭代點(diǎn)的偏導(dǎo)數(shù)值，并將迭代過程中的迭代點(diǎn)用于更新Kriging模型，兩者有機(jī)融合來求解失效概率。最后通過結(jié)合數(shù)值算例與工程算例驗證該方法的可行性和高效性。

1 Kriging方法

(1)

其中：F(β,x)為線性回歸模型；z(x)為一隨機(jī)過程；f(x)=[f1(x),f2(x),…,fp(x)]Τ(p為基函數(shù)的數(shù)目)為輸入向量x的多項式基函數(shù)，提供模擬的全局近似[29]；β=[β1,β2,…,βp]Τ為回歸系數(shù)列向量。

函數(shù)z(x)為高斯隨機(jī)過程[30]，作為局部近似[29]，其均值為0，方差為σ2，協(xié)方差滿足下式的特征：

Cov[z(x(i)),z(x(j))]=σ2R(γ,x(i),x(j))

(2)

式中：γ=[γ1,γ2,…,γn]為一個相關(guān)函數(shù)中參數(shù)向量；x(i),x(j)是試驗樣本X中任意兩個樣本點(diǎn)；R(γ,x(i),x(j))為相關(guān)函數(shù)，相應(yīng)類型有高斯函數(shù)、指數(shù)函數(shù)、冪函數(shù)、樣條函數(shù)等[31-33]。這里采用常用的高斯相關(guān)函數(shù)，其數(shù)學(xué)表達(dá)式為：

(3)

當(dāng)獲得樣本輸入X={x(1),…,x(N0)}(N0為初始樣本量)和輸出Y={y(1),…,y(N0)}后,可以計算得到式(1)中的回歸系數(shù)向量β*和式(2)中過程方差σ2[32]：

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

(4)

(5)

相關(guān)函數(shù)中的參數(shù)向量可以通過極大似然估計得到：

(6)

(7)

式中，r(x)=[R(γ,x(1),x),…,R(γ,x(N0),x)]T為輸入量x與已知樣本點(diǎn)X的相關(guān)函數(shù)，r*可以通過計算Rr*=Y-Fβ*得到。對于Kriging模型的詳細(xì)理論分析可以參照文獻(xiàn)[32-33]。

Kriging模型對輸入量和響應(yīng)量的關(guān)系可達(dá)到較高的擬合精度，且在工程中已經(jīng)得到了廣泛的應(yīng)用[23-26]。

2 基于Kriging方法與改進(jìn)一次二階矩的融合方法

所提方法采用Kriging方法高效計算偏導(dǎo)數(shù)，用于進(jìn)一步提高改進(jìn)一次二階矩的效率。由于Kriging模型中的回歸函數(shù)和相關(guān)函數(shù)都是簡單的函數(shù)體，在模型建立好后，即可簡便計算出相應(yīng)的偏導(dǎo)數(shù)，因此可將之用于改進(jìn)一次二階矩的求解過程中。同時，由于一般情況下樣本量越多，Kriging模型構(gòu)建越精確。為了充分利用已有信息，所提方法將改進(jìn)一次二階矩求解進(jìn)程中的近似設(shè)計點(diǎn)信息運(yùn)用于Kriging模型的更新構(gòu)建中，提高Kriging模型的精度，進(jìn)而將Kriging方法與改進(jìn)一次二階矩有機(jī)融合。

步驟1：構(gòu)建Kriging模型。

采用蒙特卡洛法隨機(jī)抽取N0樣本點(diǎn)X={x(1),x(2),…,x(N0)}(含均值點(diǎn))，代入結(jié)構(gòu)功能函數(shù)中計算響應(yīng)值Y={y(1),…,y(N0)}，然后根據(jù)樣本點(diǎn)構(gòu)建初始的Kriging模型。

當(dāng)回歸模型采用一階多項式，則輸入向量x的基函數(shù)向量為f(x)=[1,x]，隨機(jī)過程采用高斯過程，如式(3)所示，由第1節(jié)式(4)可知樣本量、回歸模型、隨機(jī)過程確定后，便得到回歸系數(shù)向量β*和r*，于是得到極限狀態(tài)函數(shù)g(x)初始的Kriging代理模型

(8)

步驟2：將Kriging導(dǎo)數(shù)信息運(yùn)用于AFOSM。

在原有改進(jìn)一次二階矩方法中，對原極限狀態(tài)函數(shù)在設(shè)計點(diǎn)x*處線性化：

(9)

一般采用迭代的算法來獲得[1]設(shè)計點(diǎn)，即設(shè)置初始設(shè)計點(diǎn)值，如設(shè)為均值點(diǎn)x*(0)=μ，而后在迭代過程中逐步對設(shè)計點(diǎn)值進(jìn)行更新。此外，在計算導(dǎo)數(shù)的時候，需要較多的計算代價，尤其是在變量維數(shù)較多的情形下。在本文所提方法中，采用構(gòu)建的Kriging代理模型來獲取近似導(dǎo)數(shù)信息，進(jìn)一步提升AFOSM的效率。即運(yùn)用式(8)所示的代理模型計算導(dǎo)數(shù)

(10)

(11)

由此，在后續(xù)改進(jìn)一次二階的迭代求解過程中，導(dǎo)數(shù)的計算完全由構(gòu)建的Kriging模型來計算。

步驟3：更新Kriging模型及設(shè)計點(diǎn)直至收斂。

將迭代過程中的中間迭代點(diǎn)添加到構(gòu)建初始Kriging代理模型的樣本點(diǎn)中，對前面建立的Kriging代理模型進(jìn)行更新和修正，提高Kriging模型的精度，從而提供更為準(zhǔn)確的導(dǎo)數(shù)信息，為所提融合方法的最終收斂提供保障。做法如下：

重復(fù)將步驟2中得到的設(shè)計點(diǎn)值應(yīng)用到步驟1中，更新Kriging模型的迭代過程，直到前后兩次的可靠度指標(biāo)的相對誤差滿足設(shè)定的精度要求為止，最后對應(yīng)的可靠度指標(biāo)ψ和失效概率Pf為：

(12)

Pf=Φ(-ψ)

(13)

本文所提方法的流程圖如圖1所示。

圖1 Kriging與改進(jìn)一次二階矩融合方法流程圖Fig.1 Flow chart of Kriging and AFOSM

3 算例分析

為了驗證所提方法的適用性，該節(jié)結(jié)合數(shù)值算例和屏蔽閘閥結(jié)構(gòu)的工程算例進(jìn)行說明。采用本文所提融合方法對各例進(jìn)行可靠性分析，并與其他分析方法進(jìn)行比較，包括：AFOSM，即常規(guī)改進(jìn)一次二階矩方法[1]；Kriging+MCS，即建立足夠精度的Kriging模型后再使用MCS方法[28]；MCS，即直接使用蒙特卡洛法。

本文以MCS直接方法計算的值作為精確值(通過大樣本數(shù)來計算得到)，以AFOSM 和Kriging+MCS兩種方法來與本文方法進(jìn)行對比。需要指出的是，對于可靠性計算效率的評定，一般通過極限狀態(tài)的計算(調(diào)用)次數(shù)來衡量。本文構(gòu)建Kriging模型所采用的樣本點(diǎn)是通過計算極限狀態(tài)函數(shù)得到，所以樣本點(diǎn)數(shù)也表示了極限狀態(tài)函數(shù)的計算次數(shù)。為了考察本文所提方法中初始Kriging模型構(gòu)建的隨機(jī)性，采用不同數(shù)目的樣本點(diǎn)構(gòu)建初始Kriging模型，然后進(jìn)行可靠性分析。在以下三個算例中用本文方法1和本文方法2表示不同樣本點(diǎn)下本文所提方法計算效果。

算例1非線性極限狀態(tài)函數(shù)g1(x)為：

g1(x)=x1-x2x3

(14)

其中：變量x1,x2,x3均服從正態(tài)分布，其分布信息分別為x1～N(7000,14002)，x2～N(187 500,28 1252)，x3～N(0.024,0.001 442)。

首先驗證用Kriging方法計算改進(jìn)一次二階矩偏導(dǎo)數(shù)的可行性。本文方法采用蒙特卡洛隨機(jī)抽取5個樣本點(diǎn)和均值點(diǎn)組成的初始樣本點(diǎn)構(gòu)建初始Kriging模型，由構(gòu)建的Kriging模型計算第一次迭代點(diǎn)即均值處偏導(dǎo)數(shù)信息，同時使用AFOSM中有限差分(其步長為0.001)計算偏導(dǎo)數(shù)，以下算例與此相同，計算結(jié)果如表1所示。

從表1可以看出，在均值點(diǎn)處，根據(jù)初始Kriging方法得到的偏導(dǎo)數(shù)與理論計算出來的偏導(dǎo)數(shù)的相對誤差在3%以內(nèi)，已具有足夠的精度。

采用本文方法及其他方法進(jìn)行可靠性分析的結(jié)果如表2所示。以MCS方法的結(jié)果(視作精確值)作為參照，各方法的失效概率與精確值都比較接近，誤差均小于1%。從計算效率上講，AFOSM計算極限狀態(tài)函數(shù)26次，Kriging+MCS方法中使用200個樣本點(diǎn)構(gòu)建了Kriging代理模型，而后用蒙特卡洛抽樣方法計算失效概率。本文方法收斂準(zhǔn)則為|ψi-ψi+1|/ψi≤5%，以下兩個算例中設(shè)計點(diǎn)的確定方法相同。本文方法1使用6個初始樣本點(diǎn)構(gòu)建初始Kriging模型，經(jīng)過兩次迭代修正便可計算出極限狀態(tài)函數(shù)的失效概率，計算次數(shù)為6+2。為了檢驗方法穩(wěn)健性，使用不同的初始樣本點(diǎn)進(jìn)行分析，亦列在表中，記為本文方法2，計算次數(shù)為4+2。表2中本文方法1和本文方法2計算結(jié)果一致，與MCS計算結(jié)果在誤差允許范圍內(nèi)，由此可以看出本文所提方法的效率和適用性。

表1 初始Kriging模型計算極限狀態(tài)函數(shù)的偏導(dǎo)數(shù)(均值點(diǎn))Tab.1 Partial derivative of the limit state function by the initial Kriging model (at the mean point)

注：?g1/?xi為極限狀態(tài)函數(shù)對變量的偏導(dǎo)數(shù)值。

表2 算例1不同方法得到的可靠性分析結(jié)果Tab.2 Reliability analysis results obtained by different methods in example 1

算例2非線性極限狀態(tài)函數(shù)g2(x)為：

(15)

其中：基本隨機(jī)變量x1,x2,x3,x4均服從正態(tài)分布，分別為x1～N(83.5,0.122)，x2～N(83.5,0.122)，x3～N(83.5,0.122)，x4～N(150,0.252)。

采用本文融合方法計算均值點(diǎn)處的偏導(dǎo)數(shù)值如表3所示，可靠性分析結(jié)果如表4所示。

從表3可以看出，在均值點(diǎn)處，根據(jù)初始Kriging方法得到的偏導(dǎo)數(shù)與理論計算出來的偏導(dǎo)數(shù)值相對誤差在7%以內(nèi)，與算例1相比相對誤差值偏大，原因是極限狀態(tài)函數(shù)的非線性提高，隨機(jī)變量的維度增多。

表3 基于初始Kriging模型計算極限狀態(tài)函數(shù)的偏導(dǎo)數(shù)(均值點(diǎn))Tab.3 Partial derivative of the limit state function by the initial Kriging model (at the mean point)

表4 算例2不同方法得到的可靠性分析結(jié)果Tab.4 Reliability analysis results obtained by different methods in example 2

同樣從表4可以看出，以MCS方法的結(jié)果(視作精確值)作為參照，其他三種方法的失效概率與精確值都比較接近，但是從計算極限狀態(tài)函數(shù)的次數(shù)來說，對于常規(guī)改進(jìn)一次二階矩調(diào)用次數(shù)最多為25次，Kriging+MCS方法采用200個樣本構(gòu)建Kriging模型。而本文方法采用不同的初始樣本數(shù)來進(jìn)行分析，結(jié)果亦列在表中本文方法1和本文方法2?？梢钥闯霰疚姆椒ㄓ嬎愦螖?shù)最少，效率最高。

算例3閘板可靠性分析一屏蔽閘閥的結(jié)構(gòu)如圖2所示，建立相應(yīng)的有限元來分析接觸面(閘板與密封座)最大應(yīng)力值。該結(jié)構(gòu)的相關(guān)參數(shù)有：閘板的傾斜角θ，材料彈性模型E，閘板與密封座的摩擦系數(shù)f，閘板上端的均布載荷p。假定各參數(shù)均服從正態(tài)分布，分布參數(shù)為θ～N(5，0.12)，p～N(140,102)，E～N(197 000,10002)，f～N(0.24,0.012)。以結(jié)構(gòu)的最大應(yīng)力不超過給定極限應(yīng)力，建立結(jié)構(gòu)的極限狀態(tài)函數(shù)，即

g(θ,E,f,p)=[σ]-σmax

(16)

其中：[σ]=310 MPa為極限應(yīng)力；σmax為結(jié)構(gòu)中閥板和密封座接觸面的最大應(yīng)力值可以表達(dá)為：

σmax=ANSYS(θ,E,f,p)

(17)

其中，ANSYS(θ,E,f,p)表示由調(diào)用ANSYS有限元分析得到。

(a) 正視圖(a) Front view

(b) 有限元網(wǎng)格劃分(b) Finite element meshing圖2 屏蔽閘閥結(jié)構(gòu)的三維有限元模型Fig.2 The three-dimension finite element model of shielded gate valve structure

表5給出了各方法的計算結(jié)果?？梢钥闯?，對于隱式函數(shù)，常規(guī)的AFOSM方法無法計算結(jié)構(gòu)的失效概率，從側(cè)面也反映了本文所提方法的優(yōu)勢。Kriging+MCS使用300個樣本點(diǎn)構(gòu)建了Kriging代理模型，然后通過MCS計算失效概率。本文方法計算兩次，一次調(diào)用6+6次結(jié)構(gòu)極限狀態(tài)函數(shù)，另一次調(diào)用5+4次，計算得到的失效概率值相近。同樣從該例可看出本文所提方法的效率和對解決工程結(jié)構(gòu)問題的適用性。

表5 屏蔽閘閥的可靠性分析結(jié)果Tab.5 Reliability analysis results of shielded gate valves

4 結(jié)論

本文提出了將Kriging方法和改進(jìn)一次二階相融合的可靠性分析方法。該方法針對AFOSM在處理復(fù)雜非線性函數(shù)或包含有限元分析的隱式極限狀態(tài)函數(shù)的可靠性問題時，采用Kriging方法求解極限狀態(tài)函數(shù)的偏導(dǎo)數(shù)，為AFSOM迭代計算設(shè)計點(diǎn)時提供偏導(dǎo)數(shù)信息，同時將AFSOM迭代產(chǎn)生的中間迭代點(diǎn)用于更新Kriging模型進(jìn)一步提高偏導(dǎo)數(shù)的計算精度。

通過兩個數(shù)值算例可以看出，所提方法中運(yùn)用Kriging得到求解梯度信息與理論計算的基本一致。最終計算的結(jié)果亦能滿足精度要求，而在效率上與傳統(tǒng)的AFSOM相比有了進(jìn)一步的提升。所提方法拓寬了AFOSM使用的范圍，可以用于解決復(fù)雜的工程可靠性問題，且提高了其求解效率。