基于自動改進Kriging 方法的導彈艙體結構可靠性分析

秦 強，趙朋飛，張文偉，馮蘊雯

（1.航天科工防御技術研究試驗中心，北京，100854；2.西北工業大學航空學院，西安，710072）

0 引 言

導彈的艙體結構和其他彈上設備一樣，是導彈系統的重要組成部分[1,2]。各個艙段、空氣舵和燃氣舵等各種翼面、各艙段之間和艙段與翼面之間的連接件構成了艙體主要結構。在執行機動動作過程中，導彈艙體結構會因受到發動機推力、自身重力、振動應力和慣性力等載荷的共同作用導致強度失效[3]。傳統的結構設計法通過引入安全系數以表征導彈工作中的各種不確定因素對結構的影響。然而這種方法在很大程度上依賴于設計者的經驗，具有一定的不確切性或盲目性。在導彈結構上采用此方法不僅不能深入分析影響導彈結構可靠性的不確定因素而且不能完全杜絕結構失效，相反，會造成導彈結構重量的增加，導致導彈射程減小。因此，深入研究導彈艙體結構關鍵設計參數對其影響程度并且保證導彈結構安全可靠是導彈設計過程中需要迫切解決的關鍵問題。

對于導彈結構，結構可靠性方法的顯著特點是能夠減少質量，并能降低成本和提高性能。近年來，眾多學者開展了導彈結構可靠性分析研究。文獻[3]在確定性分析的基礎上，研究了某型導彈艙段連接結構強度的可靠性；文獻[4]和文獻[5]采用Patran/Nastran 軟件建立了導彈吊掛結構參數化模型，并對其開展了可靠性即靈敏度分析；文獻[6]將某導彈空氣舵肋梁結構的尺寸以及彈性模量作為隨機變量，在將最大變形作為輸出響應量的基礎上構建極限狀態方程，進而采用自主學習Kriging 方法得到了結構的失效概率。通過以上研究可知，由有限元方法得到的結構最大應力或最大變形等響應量是艙體材料性能、外部載荷以及結構尺寸等參數的隱式表達。然而，若分析導彈艙體結構可靠性，則應構造一定的函數表達式，建立輸出響應量與輸入參數之間的顯式關系，進而得到結構的極限狀態方程。開展艙體結構有限元分析時，可以選擇

ABAQUS 、 ANSYS 、 MSC.PATRAN/NASTRAN 、HYPERWORKS 等有限元分析軟件[7]。其中，ANSYS Workbench 的六西格瑪（Six Sigma）分析模塊內嵌了多種代理模型可以建立艙體結構輸出響應量與輸入參數之間的函數關系，包括人工神經網絡、Kriging 模型、完全二次多項式響應面函數等。文獻[8]～[11]分別利用該模塊對圓盤造球機、甘蔗切割器上臂軸、門式卸車機以及機床主軸等結構進行了可靠性分析，充分表明了該模塊在可靠性分析中的有效性。以上文獻在建立響應面時用到了神經網絡模型與二次多項式響應面方法，這些方法只能通過已經生成的樣本集得到相應的響應面函數，而且僅能得到預測點的預測值[12]。與神經網絡與二次響應面法不同的是，Kriging 模型在提供預測點的預測值的同時能夠得到該點的方差[13]。近年來，國內外學者依據Kriging 模型這一特性在其主動學習方面做了大量的研究工作，并應用于結構可靠性分析領域[14]。在通過Kriging 方法擬合響應面的過程中，根據Kriging 模型的特性可以得到對響應面擬合性能產生最大影響的點，并將其作為新的設計點加入到原始的訓練集中，這樣通過一定的規則反復迭代不斷地擴充訓練樣本集，達到增強Kriging 模型擬合性能的目的，直到滿足收斂要求。這種主動學習的Kriging 方法在 ANSYS Workbench 中稱為自動改進 Kriging（Auto-Refinement Kriging，AR-Kriging）方法，其自帶的AR-Kriging 分析模塊為結構可靠性分析提供了更為便捷有效的手段。

本文利用參數化建模方法在ANSYS Workbench中構建導彈艙體結構模型，然后采用最優空間填充設計抽樣獲取多組樣本點并通過有限元分析得到各樣本點對應的響應量，結合AR-Kriging 方法生成響應面并開展參數靈敏度分析以探究結構尺寸、材料性能、外部載荷等參數對艙體結構的影響程度，進而通過六西格瑪分析評估艙體結構可靠性。

1 Kriging 模型基本理論

Kriging 模型是一種插值模型，通過參數線性回歸模型疊加一個非參數隨機過程而成[12～14]。模型表達式為

式中 Γ(β , x) 為多項式回歸模型；β 為回歸系數向量，β=[β1,...,βp]T；fT( x) 為變量 x 的多項式，fT( x) =[f1(x),f2(x),...,fp(x)]T； z ( x) 為均值為零、方差為σ2的隨機高斯過程。在設計空間不同位置處，這些隨機變量之間的相關性通過協方差表述為

式中 ? ( xi, xj; θ )為 xi與xj的相關函數。在經典Kriging模型中，常見的相關函數模型有指數模型、高斯模型、線性模型、樣條函數模型等，而最為常用的函數模型是高斯模型：

式中 θ 為參數向量，θ=[θ1, θ2,···, θm]T；m 為輸入向量第m 維元素；M 為輸入向量的總維度。

式中 F 為N0×1 的單位矩陣；N0表示訓練樣本點個數。由式（1）～（3）可知，由回歸系數向量β、隨機過程的方差σ2以及參數向量θ 可以完全定義一個Kriging 模型。而由式（4）、式（5）可知，回歸系數向量β 與隨機過程的方差 σ2依賴于參數向量 θ。因此，在構造Kriging 模型時，首先應根據樣本點求出參數向量θ，這一步驟可以通過最大似然估計實現，即：

此時，對于預測點x 的預測值?( )G x 的均值和方差分別為

2 艙體結構及參數化建模

在ANSYS Workbench 中構建艙體結構，在其前后端面各有4 個孔，用于鏈接雙頭螺栓；在艙體中部開有1 個矩形孔，如圖1 所示。在設置邊界條件時，固定后端面，而施加一定的外載荷與前端面處。通過確定性分析得到艙體結構的應力云圖和應變云圖，如圖2所示。

圖1 導彈艙體結構示意Fig.1 Diagram of a Missile Cabin Structure

圖2 艙體結構受載云圖Fig.2 Stress and Deformation Nephograms

在鑄造加工艙體結構過程中，由于儀器設備與量具的精度以及工藝人員的技術水平等外部因素的影響，即使是同一批次艙體結構，在加工成型后仍會存在一定的尺寸差異。因此在對艙體結構開展可靠性分析時將艙體關鍵尺寸參數考慮為隨機變量，包括艙體長度Y1、中部開孔長度Y2、寬度H 以及內徑X 等。此外，在考慮艙體材料性能以及加載不均勻性對其結構功能的影響情況下，將材料的彈性模量E 以及加載力F亦作為隨機變量。一般情況下，產品尺寸被認為服從正態分布，而在實際情況中，由于裝配關系，產品的尺寸存在上限值及下限值，因此屬于截斷正態分布；通過加載系統可以將外載荷很好地控制在上、下限范圍內，因此認為其也服從截斷正態分布；一般認為材料的彈性模量服從正態分布。本文涉及的6 個輸入隨機變量的數字特征如表1 所示。

表1 輸入變量數字特征Tab.1 Statistical Characteristics of Random Variables

根據以上輸入變量的數字特征，在 ANSYS Workbench 中進行設置，即可獲得導彈艙體結構參數化模型。

3 基于AR-Kriging 的艙體結構可靠性分析

3.1 隨機樣本抽樣

在ANSYSY Workbench 中，可以選擇的抽樣方法包括拉丁超立方抽樣（LHS）、最優空間填充（OSF）方法、中心復合設計（CCD）抽樣以及Box-Behnken抽樣等。其中，CCD 以及BOX-Behnken 方法適用于構造經典響應面函數，并且產生的隨機樣本數與輸入變量的個數有關。但是LHS 和OSF 方法產生的樣本數與輸入變量的個數無關，可以在Workbench 中任意設置，而且OSF 是一種優化的LHS 方法，在設計空間中通過OSF 產生的樣本的均勻性要優于LHS 產生的樣本[15]。因此，本文根據OSF 抽樣出200 個隨機樣本，這樣不僅能夠獲得足夠多的樣本訓練Kriging 模型，而且可以在后續自動改進過程中生成更優的樣本點。

3.2 采用自動改進Kriging 方法構造響應面

在Workbench 中的響應面類型中選擇“Kriging”，而在核變量類型選項中有默認選項“Variable”和“Constant”，其中，Variable 表示在式（3）中的參數向量θ 的每一分量均不同，而Constant 表示θ 的每一個分量均相同[15]，本文保持其默認選項。

在Workbench 中利用AR-Kriging 方法構造響應面時，各參數的設置如表2 所示。

表2 AR-Kriging 在Workbench 中的設置Tab.2 The setup of AR-Kriging in Workbench

在ANSYS Workbench 中，AR-Kriging 方法的原理是通過自動改進策略選擇特定的改進準則并根據自動優化程序在設計空間中增加更多的設計點以達到自動更新Kriging 模型的目的。在自動改進迭代的每一次循環中，Kriging 模型在全部參數空間范圍內對預測相對誤差進行評估。而誤差預測是一個連續可微的方程。為了得到最優的改進設計點，改進過程通過運行一個基于梯度優化的程序計算這個預測方程的最大值。如果新的候選設計點的預測精度優于設置的允許精度，那么認為該候選設計點是新的設計點。

3.2.1 雙輸出變量情況

本節將最大應力及最大變形同時作為輸出變量，分析AR-Kriging 方法在擬合響應面時的性能。采用AR-Kriging 方法對艙體結構進行分析時，輸出變量相對誤差隨改進樣本點數的變化曲線如圖3 所示。

圖3 輸出變量相對誤差隨改進樣本點數變化Fig.3 The Curve of Relative Error with Design Point

由圖3 可以看出，Kriging 模型對最大變形的擬合誤差始終能夠保持在誤差要求范圍內，但是當增加到第19 個樣本點時，該輸出變量的誤差從0.5%跳躍至4.5%，此后隨著樣本點的增加有緩慢下降的趨勢。另一方面，最大應力的擬合誤差總體上呈現出下降的趨勢，當改進樣本點增加到第78 個的時候，最大變形的擬合誤差降低到 4.6259%，滿足收斂要求，此時Workbench 停止繼續迭代。

為了對比分析AR-Kriging 模型與原始Kriging 模型的擬合精度，表3 列出了由驗證點產生的各種誤差值。其中，原始Kriging 模型是通過原始生成的200 個樣本點在Workbench 中得到的響應面模型。誤差值包括均方根誤差（RMSE）、相對均方根誤差（RRMSE）、相對最大絕對誤差（RMAE）以及相對平均絕對誤差（RAAE）。

表3 驗證點通過兩種Kriging 模型產生的誤差對比分析Tab.3 Comparison of the Errors Produced by Two Kriging Models

由表3 可知，對最大應力這一輸出變量而言，20 個驗證點通過自動改進Kriging 模型產生的RMSE、RRMSE 以及RAAE 相對于由Kriging 模型產生的相應的誤差有所減小，但是由AR-Kriging 模型產生的最大變形量的RMSE、RRMSE 以及RAAE 相比于Kriging模型均有不同程度的增加。另外通過表3 還可以看出，由自動改進Kriging 產生的RMAE 高于Kriging 產生的RMAE。以上的數據分析表明，在AR-Kriging 迭代過程中，為了降低最大應力的RMSE、RRMSE以及RAAE是以犧牲最大變形量的擬合精度為代價的，這一情況也可由圖3 中的誤差收斂曲線看出。

3.2.2 單輸出變量情況

由3.2.1 節分析可知，在對雙輸出變量同時構造響應面時，即使采用Kriging 方法就已經能夠對最大變形的擬合達到非常良好的效果，然而不論是Kriging 方法還是AR-Kriging 方法得到的關于最大變形的擬合精度都有待進一步提高。因此本節僅以最大應力為輸出變量，對比分析Kriging 方法及AR-Kriging 方法在單輸出變量情況下的擬合性能。保持AR-Kriging 參數如表2 所示，最大應力相對誤差隨改進樣本點數的變化曲線如圖4 所示。

圖4 最大應力響應面相對誤差隨改進樣本點數變化Fig.4 The Curve of Relative Error with Design Point

由圖4 可知，在增加11 個樣本點的情況下，最大應力的響應面預測誤差滿足收斂要求，此時Workbench停止繼續迭代。相比于圖3 可知，大大減小了新的樣本的加入，提高了計算效率。表4 列出了由驗證點產生的誤差值。

表4 最大應力驗證點通過兩種Kriging 模型產生誤差對比分析Tab.4 Comparison of the Errors Produced by Two Kriging Modes

由表4 可知，通過AR-Kriging 方法大大提高了響應面的擬合精度，RMAE 與RAAE 分別從采用Kriging方 法 時 的 132.67% 和 96.399% 降 低 到 了 采 用AR-Kriging 方法時的28.073%和10.451%。由于RMAE與RAAE 分別能夠表明Kriging 模型局部以及全局的擬合能力，因此通過AR-Kriging 方法提高了生成的響應面的全局擬合能力。

3.3 艙體結構靈敏度分析

由結構可靠性靈敏度的定義可知，靈敏度數值的絕對值越高表明輸入變量對輸出變量的影響程度越大。通過選擇局部靈敏度（Local Sensitivity）分析6 個輸入變量對2 個輸出變量的影響程度。艙體結構靈敏度矩形圖分別如圖5、圖6 所示。

圖5 輸入變量對最大應力的靈敏度變化Fig.5 Sensitivity of the Variables to the Maximum Stress

圖6 輸入變量對最大變形量的靈敏度變化Fig.6 Sensitivity of the Variables to the Maximum Deformation

為了清晰直觀地分析艙體結構靈敏度，圖5 和圖6中相應的靈敏度數值如表5 所示。

由圖5、圖6 及表5 可知，艙體結構最大應力和最大變形量對內徑X 的變化最為敏感，而且隨著內徑X的減小，最大應力和最大變形均會減小。同理，以上兩個輸出變量會隨著艙體長度、彈性模量和外載的增加而逐漸增大。通過靈敏度值的大小可知，6 個輸入變量對 2 個輸出變量影響程度的排序均為X>E>Y1>H>Y2>F。

3.4 艙體結構可靠性分析

通過Six Sigma Analysis 模塊，可以計算出艙體結構輸出變量的累計分布函數，根據該累計分布函數可以獲得任意樣本點對應的概率值。本文采用LHS 方法生成100 000 個樣本點獲得最大應力和最大變形量的累計分布函數，分別如圖7 和圖8 所示。

圖7 最大應力累計分布函數Fig.7 Cumulative Distribution Function of Maximum Stress

圖8 最大變形累計分布函數Fig.8 Cumulative Distribution Function of Maximum Deformation

由圖7 及圖8 可以看出，兩個輸出變量柱狀圖無較大的間隙或跳躍，說明100 000 次模擬已經足夠，且柱狀圖分布情況較好。在Six Sigma Analysis 中，有兩種獲得失效概率的方法：一種是通過累計分布函數，通過該函數可以檢測產品是否滿足可靠性要求，根據圖中的黑色樣本點擬合出輸出變量的累計分布函數曲線，曲線上的每一點的值表示的是輸出變量小于該點的概率；另一種方法是從參數概率列表中獲得概率值。在參數列表中直接輸入某一固定數值，其后顯示出的概率值是輸出變量小于該固定數值的概率，這種方法得到的概率值更為準確。分別在最大應力和最大變形量參數概率列表中輸入最大許用應力358 MPa 以及最大允許變形量0.285 mm，得到相應的概率值，如表6、表7 所示。

表6 最大應力概率Tab.6 Probability List of the Maximum Stress

表7 最大變形量概率Tab.7 Probability List of the Maximum Deformation

由表6 和表7 可知，在最大許用應力確定時，該艙體結構的強度可靠度為0.976 86，而在最大允許變形量確定時，該艙體結構的剛度可靠度為0.994 37。因此，在允許值確定的情況下，以該艙體結構的強度作為衡量其可靠性的依據更加保守。

4 結 論

本文采用ANSYS Workbench 對導彈艙體結構參數化建模后，分析了其可靠性以及參數靈敏度，通過研究可知：

a）與Kriging 方法相比，采用AR-Kriging 方法可以提高響應面的全局擬合能力，尤其是在以最大應力為單輸出變量時的全局擬合能力；

b）艙體結構最大應力和最大變形量對內徑的變化最為敏感，而彈性模量、艙體長度、開孔寬度、開孔長度和外載荷對以上兩個輸出變量的影響程度逐漸遞減；

c）在最大應力和最大應變的允許值確定的情況下，以該艙體結構的強度作為衡量其可靠性的依據顯得更加保守。