基于Kriging模型梯度解析解的改進一次二階矩方法

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

1. 南京航空航天大學 航空宇航學院，南京 210016 2. 中國運載火箭技術研究院，北京 100076

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

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

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

1 改進一次二階矩方法

(1)

(2)

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

(3)

Pf=Φ(-βR)

(4)

式中:Φ(·)為標準正態分布的累積分布函數。

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

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

Kriging代理模型是一種半參數插值技術，Sacks等[4]將Kriging模型表達為線性回歸部分和非參數部分2部分組成，即

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

(5)

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

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

i,j=1,2,…,N

(6)

式中：xi與xj為第i和第j個訓練樣本點；N表示訓練樣本點個數；σ2是隨機分布誤差的方差；R(xi,xj)為表征訓練樣本點之間的空間相關性的相關函數，表達為

(7)

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

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

(8)

(9)

(10)

Kriging模型對于結構輸入輸出關系具有較高的擬合精度，已經被廣泛運用于結構可靠性分析與優化設計中。關于模型詳細的描述，可以參閱文獻[4]。

Kriging模型表達為張量積基函數為[5,10]

(11)

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

(12)

分別對式(12)中等號兩邊輸入變量xq求導數，推導如下：

(13)

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

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

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

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

4) 采用式(4)求解結構失效概率。

4 算例分析

本節將所提的基于Kriging解析解的AFOSM法應用于數值及工程算例中以驗證所提方法的正確性。4.1節利用簡單的數值算例來驗證所推導的Kriging梯度解析函數的正確性。4.2節將所提的基于Kriging解析解的AFOSM法應用于某航空發動機渦輪盤模型的可靠性分析中，以驗證所提可靠性分析方法的正確性與高效性。在此基礎上，4.3節將其應用于以有限元隱式形式表征的某導彈舵面結構的可靠性分析中，進一步驗證本文所提方法的正確性與工程適用性。

4.1 數值算例

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

(14)

易得2個功能函數對各輸入變量偏導數在其均值點處函數值的精確解分別為

(15)

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

由表1所示的結果可以看出，本文所推導的Kriging梯度解析解具有較高的精度，通過對比驗證了所推導的基于Kriging代理模型的梯度解析公式的正確性。同時可看出，借助于Kriging代理模型的優勢，僅需要少量的訓練樣本即可以獲得高精度的結果，計算代價小。

表1 數值算例的Kriging梯度解析結果

Table 1 Results of Kriging based gradient analysis of numerical example

項目導函數值誤差/%訓練樣本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 航空發動機渦輪盤算例

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

圖2 某型航空發動機渦輪盤示意圖[14]Fig.2 Schematic diagram of an aero-engine turbine disc[14]

表2 某型航空發動機渦輪盤模型輸入變量分布信息

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

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

表3 某型航空發動機渦輪盤模型可靠性分析結果

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

4.3 飛行器舵面算例

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

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

(16)

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

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

圖3 某型飛行器舵面結構示意圖[15-16]Fig.3 Schematic diagram of an aircraft rudder structure[15-16]

表4 某型飛行器舵面結構輸入變量分布信息

Table 4 Input variable distribution information of an aircraft rudder structure

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

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

表5 某型飛行器舵面結構可靠性分析結果

Table 5 Reliability analysis results of an aircraft rudder structure

方法失效概率計算量Kriging_ANA_AFOSM0.008980Kriging_MC0.0092200

5 結 論

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

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

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