基于Kriging模型梯度解析解的改進(jìn)一次二階矩方法

李寶玉，張磊剛，裘群海，余雄慶

1. 南京航空航天大學(xué) 航空宇航學(xué)院，南京 210016 2. 中國(guó)運(yùn)載火箭技術(shù)研究院，北京 100076

一次二階矩(First Order and Second Moment, FOSM)方法是結(jié)構(gòu)可靠性分析中一種簡(jiǎn)單易操作的方法，作為一種實(shí)用方法在工程問題中得到了廣泛應(yīng)用[1-3]。FOSM方法的基本思想就是將非線性的功能函數(shù)進(jìn)行線性化，然后通過輸入變量的一階矩和二階矩來計(jì)算線性化后的功能函數(shù)的一階矩和二階矩，使用其求得結(jié)構(gòu)可靠度指標(biāo)，進(jìn)而近似得到功能函數(shù)的失效概率。隨著方法精度的研究與發(fā)展，F(xiàn)OSM方法包括均值一次二階矩(Mean Value FOSM，MVFOSM)方法和改進(jìn)一次二階矩(Advanced FOSM, AFOSM)方法。MVFOSM方法和AFOSM方法的區(qū)別在于二者線性化的點(diǎn)是不同的，前者是在輸入變量的均值點(diǎn)處對(duì)功能函數(shù)進(jìn)行線性化，而后者則是在對(duì)失效概率貢獻(xiàn)最大的點(diǎn)，即最可能失效點(diǎn)—設(shè)計(jì)點(diǎn)處線性化。由于設(shè)計(jì)點(diǎn)是對(duì)失效概率貢獻(xiàn)最大的點(diǎn)，因此在設(shè)計(jì)點(diǎn)處線性展開比在均值點(diǎn)處線性展開對(duì)失效概率的近似具有更高的精度[3]。

值得注意的是，不論是AFOSM方法還是MVFOSM方法，對(duì)于非線性功能函數(shù)問題，2種方法均需要將非線性功能函數(shù)進(jìn)行線性展開，因此均需要求解功能函數(shù)的導(dǎo)函數(shù)，獲取功能函數(shù)對(duì)各輸入變量的梯度信息。對(duì)于顯式功能函數(shù)，其導(dǎo)函數(shù)求解較為容易，而對(duì)于復(fù)雜工程結(jié)構(gòu)可靠性分析中常見的隱式功能函數(shù)問題，其導(dǎo)函數(shù)則較難求解。有限差分法是求解隱函數(shù)的導(dǎo)函數(shù)的通用方法，但對(duì)于高度非線性問題，有限差分法的步長(zhǎng)較難確定，而且計(jì)算梯度信息所花費(fèi)的計(jì)算量相當(dāng)大，特別地，對(duì)于以有限元形式表征的功能函數(shù)而言計(jì)算代價(jià)在工程問題中難以接受。

鑒于AFOSM方法的工程適用性，同時(shí)考慮到Kriging代理模型技術(shù)在結(jié)構(gòu)可靠性工程及優(yōu)化領(lǐng)域的應(yīng)用廣泛性，本文研究并提出了一種基于Kriging代理模型的AFOSM可靠性分析方法。在眾多代理模型中，Kriging模型是一種應(yīng)用較為廣泛且精度較高的代理模型[4-9]。例如：Simpson等[6]將Kriging與響應(yīng)面法的計(jì)算精度與計(jì)算效率進(jìn)行了對(duì)比，并將其應(yīng)用于航天飛機(jī)的設(shè)計(jì)中；Sakata等[7]改進(jìn)和發(fā)展了抽樣方法，并將Kriging模型應(yīng)用于薄壁管和機(jī)翼等結(jié)構(gòu)的優(yōu)化設(shè)計(jì)中；Lucifredi等[8]將動(dòng)態(tài)Kriging模型與神經(jīng)網(wǎng)絡(luò)進(jìn)行了對(duì)比，然后將其應(yīng)用于水電力系統(tǒng)的預(yù)測(cè)維護(hù)中。Kriging模型可以表示為張量積基函數(shù)的形式[10-11]，因此可以通過張量積基函數(shù)來解析地推導(dǎo)功能函數(shù)對(duì)輸入變量導(dǎo)函數(shù)的表達(dá)式，進(jìn)而利用較少的訓(xùn)練樣本得到較高精度的Kriging代理模型后，直接得到功能函數(shù)對(duì)輸入變量梯度信息的近似解，也即以較高的效率為AFOSM提供了較高精度的梯度信息，很好地解決了隱式問題梯度信息難以獲取的難題。

1 改進(jìn)一次二階矩方法

(1)

(2)

設(shè)功能函數(shù)的均值和標(biāo)準(zhǔn)差分別為μg和σg，上述線性極限狀態(tài)方程的可靠度指標(biāo)βR和失效概率Pf可以由式(3)和式(4)精確求解[3]：

(3)

Pf=Φ(-βR)

(4)

式中:Φ(·)為標(biāo)準(zhǔn)正態(tài)分布的累積分布函數(shù)。

對(duì)于輸入變量為非正態(tài)分布的情況，Rackwitz和Fiessler[12]提出了一種等價(jià)正態(tài)變量算法，簡(jiǎn)稱為R-F法，將非正態(tài)變量轉(zhuǎn)化為等價(jià)的正態(tài)變量，然后再采用AFOSM法求解可靠度指標(biāo)，進(jìn)而得到失效概率。

2 Kriging代理模型的梯度解析解

Kriging代理模型是一種半?yún)?shù)插值技術(shù)，Sacks等[4]將Kriging模型表達(dá)為線性回歸部分和非參數(shù)部分2部分組成，即

g(x)=fT(x)β+Z(x)

(5)

式中：f(x)=[f1(x)f2(x)…fm(x)]T為輸入向量x的多項(xiàng)式基函數(shù)，可以是0階、一階或二階多項(xiàng)式，提供模型模擬的全局性近似；β=[β1β2…βm]T是需要確定的回歸系數(shù)，m為線性回歸部分基函數(shù)的個(gè)數(shù)；Z(x)為隨機(jī)過程實(shí)現(xiàn)的非參數(shù)部分，提供線性回歸部分fT(x)β的偏差，服從N(0,σ2)的正態(tài)分布，協(xié)方差為

Cov[Z(xi),Z(xj)]=σ2R(xi,xj)

i,j=1,2,…,N

(6)

式中：xi與xj為第i和第j個(gè)訓(xùn)練樣本點(diǎn)；N表示訓(xùn)練樣本點(diǎn)個(gè)數(shù)；σ2是隨機(jī)分布誤差的方差；R(xi,xj)為表征訓(xùn)練樣本點(diǎn)之間的空間相關(guān)性的相關(guān)函數(shù)，表達(dá)為

(7)

式中：n是輸入向量x的維數(shù)；xik和xjk分別為向量xi和xj的第k個(gè)分量；τk是相關(guān)性參數(shù)；d決定了函數(shù)在第k個(gè)坐標(biāo)方向的光滑性。Kriging模型中的一元相關(guān)函數(shù)有多種形式[4-5]。

未知參數(shù)β和σ2可以通過式(8)和式(9)進(jìn)行估計(jì)：

(8)

(9)

(10)

Kriging模型對(duì)于結(jié)構(gòu)輸入輸出關(guān)系具有較高的擬合精度，已經(jīng)被廣泛運(yùn)用于結(jié)構(gòu)可靠性分析與優(yōu)化設(shè)計(jì)中。關(guān)于模型詳細(xì)的描述，可以參閱文獻(xiàn)[4]。

Kriging模型表達(dá)為張量積基函數(shù)為[5,10]

(11)

式中：λi(i=0,1,…,N)為常系數(shù)。在Kriging模型的多種相關(guān)函數(shù)中，Gaussian相關(guān)函數(shù)被廣泛選擇使用，即hik(xk)=exp(-τk(xk-xik)2)，τk為變量xk對(duì)應(yīng)的相關(guān)參數(shù)。式(11)中，Kriging函數(shù)的多項(xiàng)式基函數(shù)取0階，因此，本文所采用的Kriging模型解析表達(dá)式為

(12)

分別對(duì)式(12)中等號(hào)兩邊輸入變量xq求導(dǎo)數(shù)，推導(dǎo)如下：

(13)

3 基于Kriging解析解的可靠性分析方法

本文提出一種基于Kriging解析解的AFOSM可靠性分析方法，傳統(tǒng)AFOSM方法中的導(dǎo)函數(shù)梯度信息由Kriging模型解析求得。值得注意的是，要采用改進(jìn)的一次二階矩方法求解可靠度指標(biāo)和失效概率，必須先知道設(shè)計(jì)點(diǎn)。顯然對(duì)于一個(gè)給定的非線性功能函數(shù)，其設(shè)計(jì)點(diǎn)是不可能預(yù)先知道的，這就必須采用迭代或者直接尋優(yōu)的方法來進(jìn)行問題的求解，一般是通過迭代前后兩次的可靠度指標(biāo)的相對(duì)誤差滿足設(shè)定的精度要求為止。對(duì)于成熟的AFOSM的迭代求解步驟的詳細(xì)信息可查閱文獻(xiàn)[3]。本文所提方法的簡(jiǎn)要步驟如下文所述，流程圖如圖1所示。

1) 鑒于拉丁超立方抽樣(Latin Hypercube Sampling, LHS)[13]樣本具有低偏差特性，采用LHS方法產(chǎn)生N個(gè)訓(xùn)練樣本點(diǎn)構(gòu)建Kriging模型。

圖1 可靠性分析流程圖Fig.1 Flow chart of reliability analysis

4) 采用式(4)求解結(jié)構(gòu)失效概率。

4 算例分析

本節(jié)將所提的基于Kriging解析解的AFOSM法應(yīng)用于數(shù)值及工程算例中以驗(yàn)證所提方法的正確性。4.1節(jié)利用簡(jiǎn)單的數(shù)值算例來驗(yàn)證所推導(dǎo)的Kriging梯度解析函數(shù)的正確性。4.2節(jié)將所提的基于Kriging解析解的AFOSM法應(yīng)用于某航空發(fā)動(dòng)機(jī)渦輪盤模型的可靠性分析中，以驗(yàn)證所提可靠性分析方法的正確性與高效性。在此基礎(chǔ)上，4.3節(jié)將其應(yīng)用于以有限元隱式形式表征的某導(dǎo)彈舵面結(jié)構(gòu)的可靠性分析中，進(jìn)一步驗(yàn)證本文所提方法的正確性與工程適用性。

4.1 數(shù)值算例

考慮2個(gè)非線性功能函數(shù)g1(x)和g2(x)，均包含4個(gè)基本隨機(jī)變量xi(i=1,2,3,4)，且均服從正態(tài)分布xi～N(1,0.12)。功能函數(shù)表示為

(14)

易得2個(gè)功能函數(shù)對(duì)各輸入變量偏導(dǎo)數(shù)在其均值點(diǎn)處函數(shù)值的精確解分別為

(15)

運(yùn)用式(13)所推導(dǎo)的Kriging梯度解析函數(shù)進(jìn)行上述導(dǎo)數(shù)值求解，分別抽取一定數(shù)量的訓(xùn)練樣本建立2個(gè)功能函數(shù)的Kriging代理模型，代理模型建立后，將輸入變量均值點(diǎn)代入可直接得到相應(yīng)的梯度信息如表1所示。

由表1所示的結(jié)果可以看出，本文所推導(dǎo)的Kriging梯度解析解具有較高的精度，通過對(duì)比驗(yàn)證了所推導(dǎo)的基于Kriging代理模型的梯度解析公式的正確性。同時(shí)可看出，借助于Kriging代理模型的優(yōu)勢(shì)，僅需要少量的訓(xùn)練樣本即可以獲得高精度的結(jié)果，計(jì)算代價(jià)小。

表1 數(shù)值算例的Kriging梯度解析結(jié)果

Table 1 Results of Kriging based gradient analysis of numerical example

項(xiàng)目導(dǎo)函數(shù)值誤差/%訓(xùn)練樣本Ng1(x)/x11.0000g1(x)/x23.0000g1(x)/x35.0000g1(x)/x47.000020g2(x)/x12.0030.15g2(x)/x25.9990.017g2(x)/x39.9970.03g2(x)/x414.0020.01425

4.2 航空發(fā)動(dòng)機(jī)渦輪盤算例

本算例分別抽取50、100、300個(gè)訓(xùn)練樣本點(diǎn)建立Kriging代理模型，獲取相應(yīng)的梯度解析結(jié)果，進(jìn)而采用AFOSM法進(jìn)行渦輪盤可靠性分析。從表3可以看出，3種情況下，采用本文所提的基于Kriging解析解的AFOSM方法的可靠性分析結(jié)果均與Monte Carlo(MC)方法比較吻合，驗(yàn)證了該方法的正確性。隨著訓(xùn)練樣本的增加，Kriging模型的精度也隨之增加，相應(yīng)的梯度解析解精度也隨之增加。然而，隨著訓(xùn)練樣本數(shù)量的增加，可靠性分析結(jié)果的精度提高程度很小，可見Kriging模型本身已經(jīng)具有較高的精度，然而其可靠性分析結(jié)果的結(jié)果受限于AFOSM方法本身的精度。

圖2 某型航空發(fā)動(dòng)機(jī)渦輪盤示意圖[14]Fig.2 Schematic diagram of an aero-engine turbine disc[14]

表2 某型航空發(fā)動(dòng)機(jī)渦輪盤模型輸入變量分布信息

Table 2 Input variable distribution information of an aero-engine turbine disc

輸入變量分布類型均值變異系數(shù)ρ/(kg·m-3)正態(tài)8 2400.1C/(kg·m)正態(tài)5.670.1A/m2正態(tài)6.2×10-30.1J/m4正態(tài)1.22×10-40.1n/(r·s-1)正態(tài)2000.05

表3 某型航空發(fā)動(dòng)機(jī)渦輪盤模型可靠性分析結(jié)果

表3中所示方法的計(jì)算量也顯示出基于Kriging解析解的AFOSM法的計(jì)算效率，針對(duì)一般工程問題，只需少量的計(jì)算成本即可得到較高精度的可靠性分析結(jié)果，很大程度降低了計(jì)算代價(jià)。

4.3 飛行器舵面算例

圖3為某型飛行器舵面結(jié)構(gòu)示意圖[15-16]，其由5條桁條、6條翼肋、蒙皮以及用于和機(jī)體相連接的固定端所組成。各個(gè)部件的材料均為某鈦合金材料，材料的楊氏模量為E，泊松比為ν。翼肋與桁條的厚度均為trs，蒙皮的厚度為tsk。在飛行器的飛行過程中，舵面用于調(diào)節(jié)以及穩(wěn)定飛行姿態(tài)，其所承受的載荷主要為氣動(dòng)載荷，大小為P，垂直作用于上蒙皮。舵面結(jié)構(gòu)材料的彈性模量、泊松比，桁條、翼肋、蒙皮的厚度以及所承受的氣動(dòng)載荷均為服從正態(tài)分布的隨機(jī)變量，詳細(xì)分布信息見表4。以舵面結(jié)構(gòu)在承受氣動(dòng)載荷情況下的最大變形不超過δc=8 mm為門限值，建立結(jié)構(gòu)系統(tǒng)的功能函數(shù)為

g(E,ν,P,tsk,trs)=δc-max(δ)

(16)

由于此算例為以有限元形式表征的隱式問題，舵面結(jié)構(gòu)在氣動(dòng)載荷作用下的位移無法解析求解，本節(jié)采用有限元分析軟件MSC.Patran/Nastran對(duì)結(jié)構(gòu)的變形進(jìn)行求解，有限元模型如圖4[15-16]所示。一次響應(yīng)計(jì)算時(shí)間較長(zhǎng)，若采用Monte Carlo法進(jìn)行可靠性分析，其計(jì)算代價(jià)是無法接受的。為了結(jié)果對(duì)比，在建立Kriging模型后，采用MC方法(Kriging_MC)針對(duì)建立的代理模型進(jìn)行結(jié)構(gòu)可靠性分析，以此結(jié)果作為參照解進(jìn)行對(duì)比。

如圖5所示，抽取200個(gè)訓(xùn)練樣本建立舵面結(jié)構(gòu)功能函數(shù)的Kriging代理模型，進(jìn)而抽取106個(gè)樣本代入建立的代理模型中進(jìn)行舵面結(jié)構(gòu)的可靠性分析，以此結(jié)果作為參照解。在抽取80個(gè)訓(xùn)練樣本的情況下，采用本文所提方法的可靠性分析結(jié)果與參照解較為接近，再次驗(yàn)證了所提結(jié)構(gòu)可靠性分析方法的正確性，結(jié)果的差異性主要來源于AFOSM方法本身對(duì)于高度非線性問題的適用性。然而值得指出的是，若直接采取數(shù)值模擬法進(jìn)行有限元隱式問題的可靠性分析的計(jì)算代價(jià)是無法承受的，進(jìn)行200次有限元模型分析同樣需要花費(fèi)較長(zhǎng)時(shí)間，計(jì)算成本相對(duì)較高，本文方法在80次模型調(diào)用的情況下，所得結(jié)果精度較高，進(jìn)而進(jìn)一步驗(yàn)證了本文所提方法的工程適用性。

圖3 某型飛行器舵面結(jié)構(gòu)示意圖[15-16]Fig.3 Schematic diagram of an aircraft rudder structure[15-16]

表4 某型飛行器舵面結(jié)構(gòu)輸入變量分布信息

Table 4 Input variable distribution information of an aircraft rudder structure

輸入變量分布類型均值變異系數(shù)E/MPa正態(tài)1176000.1ν正態(tài)0.30.1P/MPa正態(tài)0.160.1tsk/mm正態(tài)10.05trs/mm正態(tài)30.05

圖4 某型飛行器舵面結(jié)構(gòu)有限元模型[15-16]Fig.4 Finite element model for an aircraft rudder structure[15-16]

表5 某型飛行器舵面結(jié)構(gòu)可靠性分析結(jié)果

Table 5 Reliability analysis results of an aircraft rudder structure

方法失效概率計(jì)算量Kriging_ANA_AFOSM0.008980Kriging_MC0.0092200

5 結(jié) 論

1) 針對(duì)工程廣泛應(yīng)用的AFOSM法對(duì)于隱式功能函數(shù)問題存在梯度信息較難求解的現(xiàn)狀，尤其是對(duì)于以有限元形式表征的工程復(fù)雜問題，梯度計(jì)算難度大、計(jì)算代價(jià)大，本文提出了基于Kriging模型梯度解析解的AFOSM可靠性分析方法，為AFOSM法提供了高精確的梯度信息，有效地解決了AFOSM方法的工程應(yīng)用難題。

2) 鑒于Kriging模型在工程應(yīng)用中的優(yōu)勢(shì)以及考慮到Kriging模型可以很好地表征為張量積基函數(shù)的形式，本文研究推導(dǎo)出了基于Kriging代理模型的功能函數(shù)對(duì)輸入變量導(dǎo)函數(shù)的解析表達(dá)式，Kriging模型建立好后可直接得到梯度解析結(jié)果。

3) 本文提出的基于Kriging代理模型梯度解析解的AFOSM法拓寬了AFOSM方法的應(yīng)用范疇，成為AFOSM方法的補(bǔ)充，為工程復(fù)雜問題的結(jié)構(gòu)可靠性分析提供了可供選擇的方法。算例驗(yàn)證了基于Kriging模型的梯度解析結(jié)果具有較高精度，因此，基于Kriging模型解析解的AFOSM方法的精度與適用性與AFOSM本身密切相關(guān)。