基于樹形馬氏鏈模型的可靠性分析方法

張洪銘，顧曉輝,*，邸憶,2

1. 南京理工大學 機械工程學院，南京 210094 2. 武昌理工學院 信息工程學院，武漢 430223

工程實踐中的系統(tǒng)可靠性分析過程具有設計變量多、低失效率和極限狀態(tài)函數(shù)形式復雜甚至隱式的特點。失效概率(Failure Probability)Pf的準確計算是機械、電子系統(tǒng)可靠性分析的關鍵，Pf的本質(zhì)是設計變量的聯(lián)合概率密度函數(shù)在失效區(qū)域上的積分：

(1)

式中：X為各設計變量x構成的基本隨機向量；G(x)為系統(tǒng)的極限狀態(tài)函數(shù)(Limit State Function，LSF)，G(x)≤0表征失效區(qū)域；fX(x)為設計變量的聯(lián)合概率密度函數(shù)。前述工程實踐問題的特點使得直接處理式(1)的難度極大，難以直接積分獲得精確的系統(tǒng)失效概率，因此針對強非線性問題的高效模型是目前工程可靠性分析研究的關鍵。

目前的4類系統(tǒng)可靠性分析方法[1]中，以蒙特卡羅模擬(Monte Carlo Simulation, MCS)為基礎的數(shù)字模擬算法原理簡潔、適應性強，對問題的維度不敏感，適合用于處理考慮不確定性的系統(tǒng)可靠性分析問題。進一步地，針對原始MCS方法在實際應用中效率低下的問題，重要抽樣(Importance Sampling，IS)法[2-3]利用方便采樣的重要抽樣分布，在靠近失效域邊界的設計點處進行抽樣，使得對失效概率貢獻大的區(qū)域內(nèi)的樣本有較大概率被抽到，增大失效樣本的比重從而提升抽樣效率。確定設計點位置和構造最優(yōu)的重要抽樣分布是提升IS算法性能的關鍵，為此學者們進行了廣泛的研究。馬紀明等[4]提出的條件遞歸方法能在系統(tǒng)的LSF為隱式時快速找到設計點；為了避免重要抽樣分布的抽樣中心落入局部最優(yōu)區(qū)域，文獻[5]借助模擬退火算法尋找全局最優(yōu)設計點。上述文獻中較少提及重要抽樣分布的構造過程，實際上，不合適的重要抽樣分布容易導致估計結果的方差很大[6]。為了避免因主觀選擇重要抽樣分布類型而導致的算法性能下降，同時減少尋找設計點的計算開銷，自適應重要抽樣法應運而生。侯本偉等[7]采用基于演化過程和交叉熵模型的多準則迭代方法來計算重要抽樣函數(shù)，使預抽樣次數(shù)下降為原來的1/40；陳向前等[8]提出的主動導引技術，使抽樣中心能快速自適應地逼近設計點；唐承等[9]提出的粒子群序列二次規(guī)劃混合算法有較強的全局搜索能力；Au和Beck[10]將馬爾可夫鏈與IS算法結合，構造收斂規(guī)則，使馬氏鏈收斂后的穩(wěn)定分布與最優(yōu)重要抽樣分布相同或近似，用收斂后的狀態(tài)點進行核密度估計從而得到最優(yōu)的重要抽樣分布密度函數(shù)。該方法中馬氏鏈的收斂過程即是狀態(tài)點在失效域中的自適應探索過程，適用于失效域形狀復雜的情況。戴鴻哲等在Au工作的基礎上進行了深入研究，其先在文獻[11]中借助快速高斯變換大幅降低核密度估計的計算復雜度，有效縮短了計算耗時；而后在文獻[12]中以馬氏鏈收斂后的狀態(tài)點作為支持向量機的數(shù)據(jù)樣本，擬合出近似極限狀態(tài)函數(shù)，有效解決了失效域復雜、極限狀態(tài)函數(shù)隱式條件下的可靠性計算問題。

在前述可靠性分析研究的基礎上，為了進一步提高處理強非線性問題的效率，提出了基于樹形馬氏鏈重要抽樣(Tree Markov Chain Importance Sampling, TMCIS)模型的系統(tǒng)可靠性分析方法。文中提出的樹形馬爾可夫鏈算法是對原始馬氏鏈的改進，新的抽樣機制使其具有局部多鏈并行和自適應樣本轉移的特性，能夠從抽樣起點自動向失效域邊界搜索，并在關鍵樣本密集區(qū)域內(nèi)靈活地生成支鏈來提升抽樣效率。由樹形馬氏鏈生成的樣本，能夠最大程度地反映出失效分布的特征信息，使得核密度估計的結果能最大程度地近似于真實失效分布。文末給出了2個強非線性數(shù)值算例和BAT子彈藥彈翼的工程算例，計算結果顯示本文提出的算法在提高抽樣效率及計算結果準確度的同時具有較好的穩(wěn)定性。

1 樹形馬氏鏈模型

1.1 Metropolis-Hasting算法

滿足了細致平穩(wěn)條件(Detailed-balance Condition)的非周期馬氏鏈會收斂到一個平穩(wěn)分布，分布形式與轉移矩陣或轉移函數(shù)有關。MCMC(Markov Chain Monte Carlo)方法利用了馬氏鏈的這一特性，即通過構造馬氏鏈的轉移矩陣或函數(shù)，使得鏈的平穩(wěn)分布類型為目標分布類型或與目標分布類型十分近似，則馬氏鏈收斂后的狀態(tài)點即可看作從目標分布中抽得的樣本點。M-H算法(Metropolis-Hasting Method)[13]是MCMC方法的常見形式，其算法步驟為

步驟1給定馬氏鏈的初始狀態(tài)X0=x0，其中x0∈R，初始化計數(shù)變量t、設置總樣本點數(shù)N，這里的t,N∈N。

步驟2對t=1,2,…,N循環(huán)以下過程進行采樣：

1)第t次抽樣，馬氏鏈狀態(tài)為Xt=xt，依建議分布產(chǎn)生候選點y～q(x|xt),

2)從均勻分布抽得一個隨機數(shù)：u～U(0,1)，計算步驟1)中候選點y的接受率：

(2)

式中:p(·)為先驗分布；q(·|·)為建議分布；α為接受率。

3)如果α(xt,y)>u，則接受轉移，令xt+1=y；否則拒絕y，令xt+1=xt。

分析上述算法，在馬氏鏈滿足細致平穩(wěn)條件的前提下，若建議分布為對稱分布，亦即q(xt|y)=q(y|xt)，則式(2)中的接受率可以進一步寫為：α(xt|y)=min{1,p(y)/p(xt)}，可見當候選點處的目標分布概率密度越大時，候選點越容易被選中；同時，接受率的形式保證了概率密度較小的樣本點也有被選中的可能，從而防止馬氏鏈“陷入”目標分布密度函數(shù)的局部峰值區(qū)域，這增強了算法對狀態(tài)空間的遍歷性能，使抽到的樣本點更加貼合目標分布的分布特點。

原始M-H算法由于引入了馬氏鏈，使得抽樣過程具有自適應特性，但當狀態(tài)空間維度較高時，概率密度稀疏處的樣本很容易被拒絕，且M-H算法每次只決定一個點被接受與否，使抽樣效率很低；同時，馬氏鏈中相鄰狀態(tài)的產(chǎn)生方式?jīng)Q定了相鄰樣本點間是相關的，為了滿足樣本的獨立性要求，普遍的做法是從已經(jīng)收斂的馬氏鏈中，每隔r個點選取一個點形成新的鏈，若總樣本點數(shù)為N，則原始馬氏鏈的抽樣次數(shù)至少應為[N(r+1)-r]，即原始馬氏鏈的抽樣次數(shù)會成倍增長，需要更多的時間和計算成本來完成抽樣，導致抽樣效率進一步下降。

1.2 樹形馬氏鏈算法

針對上述不足，樹形馬氏鏈算法充分利用相鄰狀態(tài)點間的關系信息，對單次抽樣步的候選點數(shù)量進行合理擴充，允許當前狀態(tài)點可以產(chǎn)生多個候選點，使得單次抽樣步中被接受候選點的數(shù)量的平均值大于1。樹形馬氏鏈的抽樣機制，使得抽樣結果在生成關系上呈現(xiàn)出以抽樣起點為根節(jié)點的“樹狀”結構，有效樣本密集的區(qū)域會有多個候選點被接受為樣本點，這是樹結構中產(chǎn)生“枝杈”的直接原因。樹形馬氏鏈中的支鏈增強了算法對設計域充分搜索的能力，提高了抽樣效率。為了完整、充分地表達算法內(nèi)容，提出了基本向量和跳轉向量的概念來實現(xiàn)樹形馬氏鏈算法模型的抽樣機制，概念的具體定義為：

定義1(基本向量) 指馬氏鏈的當前狀態(tài)點依據(jù)建議分布產(chǎn)生候選點(該候選點稱為基本候選點)后，由當前狀態(tài)點指向基本候選點的向量。

定義2(跳轉向量) 由基本向量在狀態(tài)空間中經(jīng)過若干次變換后得到的向量。跳轉向量必須以當前狀態(tài)點為起點，則其末端指向的點稱為擴充候選點，跳轉向量與基本向量的L2范數(shù)相等。

給出樹形馬氏鏈抽樣步驟如下：

步驟2在當前樣本量的值小于N時，令t=1,2,…；循環(huán)下列步驟采樣：

2)依據(jù)步驟1)中產(chǎn)生的成對父子點對，計算父點的基本向量：

(3)

3)計算每個候選點的接受率：

(4)

樹形馬氏鏈算法對由基本向量生成跳轉向量的具體實現(xiàn)方法不做強制要求(滿足定義2的內(nèi)容即可)，以適應不同實際情況下的具體問題。本文給出了設計變量為二維隨機變量時，采用快速生成跳轉向量并得到候選點的方法?；鞠蛄康男D動作依靠變換矩陣實現(xiàn)，二維時的變換矩陣為

(5)

(6)

圖1為設計變量為二維時，樹形馬氏鏈抽樣機制示意圖，圖中2條極限狀態(tài)函數(shù)曲線將設計域分割成4個區(qū)域，左下方區(qū)域為非失效區(qū)域，候選點的產(chǎn)生、接受或拒絕狀態(tài)如圖所示；圖1中的抽樣過程所得的點按照生成關系和接受、拒絕狀態(tài)可以表達成圖2所示的樹形結構，這是樹形馬氏鏈算法的命名由來。

圖1 TMC抽樣機制Fig.1 Mechanism of TMC sampling

圖2 TMC樣本數(shù)據(jù)結構Fig.2 Data structure of TMC samples

2 核密度估計

2.1 重要抽樣法

原始IS算法是對MCS的改進，是一種典型的方差縮減技術。MCS算法的抽樣中心在遠離失效域邊界的均值點處，在低失效率問題下MCS算法抽到失效樣本的概率很小，因而需要大量抽樣次數(shù)才能有少數(shù)失效樣本被抽中。為了減少抽樣次數(shù)，IS算法利用測度變換增大極限事件的發(fā)生概率，亦即使對失效率貢獻大的區(qū)域以較大概率被抽中。設計點是失效域附近密度函數(shù)值最大的點，因此IS算法選擇設計點為抽樣中心來提高抽樣效率?；谏鲜鯥S方法的原理，可將式(1)中失效概率的積分區(qū)域從失效域改寫到整個設計空間：

(7)

式中：Θ∈RN為設計空間；IG(x)為示性函數(shù)，表征樣本是否在失效域中；當G(x)≤0時，有IG(x)=1；當G(x)>0時，IG(x)=0。引入的hX(x)是關于隨機變量的重要抽樣密度函數(shù)，hX(x)需要滿足歸一性、非負性且其分布類型應當便于采樣?？紤]式(7)的統(tǒng)計意義可以得到失效概率的期望為

(8)

表明重要抽樣法求得的失效概率的估計值為無偏估計，此處可得最優(yōu)重要抽樣分布為

(9)

(10)

由于樣本方差依概率收斂于母體方差，即有

(11)

(12)

2.2 重要抽樣分布的核密度估計

使用樹形馬氏鏈抽樣得到表征失效域的樣本點之后，需要通過樣本點得到重要抽樣分布，本文采用描述能力較強的自適應核密度估計(Kernel Density Estimation, KDE)算法，來構建重要抽樣分布密度函數(shù)。核密度估計[14-15]的本質(zhì)是用一系列核概率密度函數(shù)的組合來產(chǎn)生合適的目標分布，形式為[16]

(13)

式中：θ(j)(j=1,2,…,M)為樹形馬氏鏈生成的樣本點；這里的n為設計變量的維數(shù)；ω為窗口寬度；λi為局部帶寬因子，用來調(diào)整窗口寬度的值；K(·)為核概率密度函數(shù)，通常取正態(tài)分布類型。

窗口寬度參數(shù)ω的本質(zhì)是各個核概率密度函數(shù)進行組合時的權重，其決定了最終的核函數(shù)的平滑性。過大的ω值會使核函數(shù)過于“光滑”而丟失局部特性；過小的ω值會導致核函數(shù)尾部因為疊加效應而產(chǎn)生波動(或稱為噪聲)。文獻[17]中給出了較好的值選擇：

(14)

式中:Md為不同樣本的個數(shù)(Md

(15)

式中：α∈[0,1]，通常取α=0.5[18]；當α=0時，自適應核密度估計變?yōu)楣潭ù翱趯挾群嗣芏裙烙嫛?/p>

2.3 TMCIS可靠性分析方法流程

提出了基于TMCIS模型的系統(tǒng)可靠性分析方法，分析流程如圖3所示。其基本思想是：通過樹形馬氏鏈，快速生成盡可能服從真實目標分布的隨機樣本后，對這些樣本進行核密度估計得到重要抽樣分布函數(shù)，從獲得的分布函數(shù)進行抽樣，得到隨機的設計狀態(tài)樣本，判斷每個樣本下系統(tǒng)對應的工作狀態(tài)，來評估系統(tǒng)的失效概率并計算該次評估的方差。

圖3 TMCIS可靠性分析方法流程Fig.3 Flowchart of TMCIS reliability analysis method

3 算例分析

本節(jié)通過2個數(shù)值算例和1個工程算例，對所提算法的適應性、計算效率和準確度進行檢驗。選用了其他算法與本文方法進行對比分析，其中包括：原始蒙特卡羅(MC)算法、一次二階矩法(FOSM)及改進一次二階矩法(AFOSM)、原始重要抽樣法(IS)、Metropolis-Hasting算法(M-H)。各算例結果如表1所示，N為各方法的樣本總量和TMCIS方法中樹形馬氏鏈的抽樣數(shù)，M為TM-CIS方法的重要抽樣數(shù)；表中數(shù)據(jù)除AFOSM的

表1 算例1各方法計算結果Table 1 Calculation results of methods in Example 1

結果為公式推導所得，其余均為50次計算結果均值。本文提出的TMCIS方法在各算例中均使用高斯核進行自適應核密度估計。

算例1選取非線性程度較強的LSF檢驗TMCIS方法對強非線性多維問題的處理能力。所選極限狀態(tài)函數(shù)為

g(θ)=3-Y+(4X)4

(16)

這里θ=(X,Y)為隨機向量，2個隨機變量X、Y獨立同分布且均服從標準正態(tài)分布，系統(tǒng)有一個設計點位于(0,3)處。TMC算法中m=4，跳轉向量采用1.2節(jié)中變換矩陣法生成。TMCIS模型對本例抽樣結果如圖4所示，樣本點自適應地從起點轉移至設計點周圍，圖中帶數(shù)字標簽的等高線為KDE方法得到的重要抽樣分布(重要抽樣樣本量與TMC樣本量皆為800)，KDE樣本為該分布生成的重要抽樣樣本。

本算例使用了4種方法分別進行了計算，結果如表1所列，以MC算法的失效率估計值作為準確解，其中ε是其他3種算法計算結果與準確解間的相對誤差，CALL為調(diào)用極限狀態(tài)函數(shù)的次數(shù)。

圖4 算例1的抽樣結果Fig.4 Sampling results of Example 1

表1數(shù)據(jù)表明TMCIS算法在少樣本量下計算結果的相對誤差是AFOSM算法的1/5，這是由于AFOSM算法本質(zhì)上是將非線性LSF在設計點處進行泰勒級數(shù)展開，僅取線性部分來近似計算失效概率。因此，當LSF的非線性程度較高時，AFOSM算法會產(chǎn)生較大誤差，而TMC的抽樣機制使得其抽得的樣本盡可能地反映LSF的非線性信息，因此，相較而言TMCIS算法更適合處理非線性問題且能夠在較少的樣本量下得到較高準確度的計算結果。

算例2某串聯(lián)系統(tǒng)[10]的極限狀態(tài)函數(shù)為

g(x)=min(g1,g2)

(17)

式中:串聯(lián)的2個子系統(tǒng)極限狀態(tài)函數(shù)分別為

表2中給出了上述各算法在不同樣本量時的計算結果，圖6反映了表2中計算結果隨樣本量的變化規(guī)律。圖6中水平虛線處的縱坐標為Pf的精確值，虛線上下兩側的點線分別是誤差為±5%的位置，從圖中可以看出算法的準確度和穩(wěn)定性方面TMCIS算法最佳，IS算法次之，MC算法最差。并且TMCIS算法相較于IS算法在樣本量較小時即可得到高準確度的結果。另外，從表2中可以發(fā)現(xiàn)當N/M的值接近1且N+M的值在800～10 000之間時TMCIS算法性能最優(yōu)。

圖5 不同方法抽樣的結果Fig.5 Sampling results by different methods

表2 不同樣本量下3種算法計算失效概率結果Table 2 Calculation results of failure probability by three methods under different sample sizes

算法Pf/10-5(Cov)N/104=0.020.040.080.10.20.40.81MC7.90 (0.21)7.50 (0.30)IS7.69 (0.37)5.87 (0.24)6.81 (0.21)6.50 (0.19)6.72 (0.11)7.04 (0.17)6.83 (0.23)6.59 (0.18)TMCIS(M=800)6.84(M=200)6.70 (0.24)6.77 (0.11)6.73 (0.23)6.96 (0.12)6.79 (0.34)7.21 (0.19)6.98 (0.25)算法Pf/10-5N/104=24810204080100MC7.95(0.19)7.76(0.12)7.22(0.25)7.29(0.16)7.05(0.21)7.15(0.22)7.23(0.15)6.97(0.14)IS6.82(0.09)6.80(0.13)6.79(0.14)TMCIS(M=800)6.95(0.17)6.90(0.23)7.02(0.22)

圖6 3種方法計算結果隨樣本量的變化趨勢Fig.6 Variation of calculation results with sample sizes of three methods

算例3BAT子彈藥(Brainpower Anti-Tank submunition)是一種智能的末制導反坦克彈藥，采用聲波/毫米波復合導引手段來追蹤、定位目標[19]。BAT子彈藥的彈體中部有4個與彈體縱軸垂直、呈十字形分布的彈翼，4片彈翼的末端安裝了聲信號傳感器，形成一個平面四元聲傳感器陣列，如圖7所示。

由于要保證傳感器陣列的測量精度，則各傳感器之間的相對距離應保持在誤差范圍內(nèi)。為了滿足減輕總重和傳感器布線的要求，彈翼框架的主體結構為由28根桿件和16片板殼剛性連接構成的三段盒式結構[20]。所有桿件的橫截面積A相同，T為板的厚度；以L1、L2、L3分別表示X、Y、Z方向上桿的長度；桿和板的彈性模量為E，泊松比為常量取0.3；P為加在節(jié)點上的外載荷。翼盒結構和外載荷加載如圖8所示。

圖7 BAT子彈藥外觀結構Fig.7 Structure appearance of BAT submunition

取L1、L2、L3、A、E、P、T為隨機變量，各隨機變量的分布類型和分布參數(shù)如表3所列。距離聲傳感器最近的16號節(jié)點在Y方向上的最大允許位移為4.5 mm，以16號節(jié)點的Y向位移建立極限狀態(tài)函數(shù)：

g(X)16#,Y=[δ]-δmax(X)

(18)

式中：X為由各隨機變量構成的隨機向量，即有：X=[L1,L2,L3,A,E,P,T]T;δmax(X)為各基本隨機變量為X狀態(tài)時16號節(jié)點的最大位移；TMC算法的m=8，8個子點中包含一個基本候選點和7個擴充候選點。不同于算例1、2使用變

圖8 翼盒結構及外在加載示意圖Fig.8 Schematic of wing box structure and external loading

表3 翼盒結構隨機變量分布信息

Table 3 Random variable distribution information of wing box structure

隨機變量均值變異系數(shù)分布類型L1/m0.60.10正態(tài)分布L2/m0.20.10正態(tài)分布L3/m0.40.10正態(tài)分布A/m20.00010.10正態(tài)分布E/GPa710.10正態(tài)分布P/N15000.10正態(tài)分布T/m0.0030.10正態(tài)分布

換矩陣生成跳轉向量，本算例依靠基本向量的坐標輪轉生成跳轉向量。坐標輪轉矩陣為

7個跳轉向量依下列遞推式生成：

表4中給出了各方法的失效概率計算結果[21]以及算得該結果所用的樣本點數(shù)。從表中數(shù)據(jù)可以看出對非線性部分忽略最嚴重的AFOSM算法的計算結果誤差最大。文中提出的TMCIS算法相較于文獻[22]的算法和MCS算法，用更少的樣本點數(shù)得到了相對準確的結果。

表4 算例3失效概率計算結果[22]

4 結 論

1) 文中通過分析算法原理，指出了FOSM算法及AFOSM算法處理強非線性問題時誤差大的根源，同時針對原始MCMC算法單次迭代中，僅抉擇一個候選狀態(tài)接受與否帶來的效率問題，提出了TMCIS模型。TMCIS的多候選點抽樣機制保證了模型抽得的樣本能充分自適應地探索設計域，在失效概率密度大的區(qū)域進行局部多鏈并行抽樣，使得模型搜索效率提高的同時加快收斂速度。

2) 為了提高模型的穩(wěn)定性、魯棒性以更加適應帶有強非線性特點的工程問題，TMCIS模型將TMC抽樣算法與自適應核密估計結合，使得到的重要抽樣分布與真實失效分布達到最大程度的近似。文末算例顯示出了所提算法在不同樣本量下，對超低失效率、串聯(lián)系統(tǒng)問題的計算結果均能收斂到真值的范圍內(nèi)，驗證了TMCIS算法的穩(wěn)定性和可信度。

3) TMC抽樣算法作為TMCIS模型的核心部分，其所具有的多候選狀態(tài)抽樣機制，在提高遍歷性能以解決非線性問題的同時，能夠進一步利用多個候選狀態(tài)中攜帶的額外失效域信息，使得核密度估計的結果更加接近真實失效分布。這是文中所列算例中基于TMCIS模型的可靠性分析方法能以較少的樣本點得到準確失效率的原因。