馬永亮 曲先強 崔洪斌 石德新

（哈爾濱工程大學(xué)多體船技術(shù)國防重點學(xué)科實驗室 哈爾濱 150001）

在結(jié)構(gòu)可靠性分析中最重要的內(nèi)容是求解結(jié)構(gòu)的失效概率.一般來說，除了簡單的極限狀態(tài)函數(shù)外，結(jié)構(gòu)失效概率的計算都比較困難.尋求高效、精確的計算方法是結(jié)構(gòu)可靠性研究的一個重要內(nèi)容.Monte Carlo方法具有模擬的收斂速度與隨機變量的個數(shù)無關(guān)，極限狀態(tài)函數(shù)的復(fù)雜程度與模擬過程無關(guān)的優(yōu)點［1－2］，是一種應(yīng)用廣泛的數(shù)值計算方法.直接 Monte Carlo法（crude monte carlo method）的誤差與無偏統(tǒng)計量的方差成正比，與模擬次數(shù)的平方根成反比.所以需要采用縮減方差抽樣方法來提高計算效率［1－4］.Harbitz提出了一種特殊的Monte Carlo法，稱為β球法［5－6］.β球法在安全域剔除了一個球形或超球形的抽樣區(qū)域，特別適用于極限狀態(tài)函數(shù)為凸曲面的情況.孫海虹［7］、張曉軍［8］、吳亞舸［9］采用β球法進行了結(jié)構(gòu)可靠性研究.雖然β球法得到了廣泛的應(yīng)用，但這種方法一般需要在求解之前給出β球的半徑.所以，這種方法不能獨立使用.基于這一點本文根據(jù)β球法自身求解特點提出了一種自適應(yīng)抽樣算法來確定初始β球的大小，使β球法不再依賴于其他方法，并將此方法應(yīng)用于潛水器耐壓圓柱殼結(jié)構(gòu)可靠性分析中.

1 一般β球方法

1.1 一般β球方法的基本原理

β球方法必須滿足以下的假設(shè)：（1）基本隨機變量服從正態(tài)分布；（2）可靠性指標(biāo)β的估計值必須事先知道.

β球?qū)⒖臻g分為了2個部分｜X｜≤β和｜X｜＞β，其中｜X｜≤β為中心在原點的m維超球體.根據(jù)全概率公式，失效概率Pf可以表示為

式中：χm為自由度為m的χ2分布函數(shù).

失效概率Pf的一個無偏估計表示為

式中：k為樣本的大小；I［G（Yi）］為示性函數(shù)，在安全域I＝0，在失效域I＝1；Yi為標(biāo)準(zhǔn)正態(tài)空間的一組隨機變量.

1.2 β球方法抽樣策略

基本隨機向量X抽樣可以表示為

式中：Ri為Xij的模，有；αij為隨機方向矢量，滿足

由于Ri和αij之間是相互獨立的，所以可以對Ri和αij分別進行抽樣，具體的抽樣方法見文獻［5，6，10］.

1.3 樣本變換及其改進

從式（2）可以看出，失效概率的求解只與示性函數(shù)有關(guān).只要求解出示性函數(shù)的結(jié)果就可以了.所以可以將相互獨立的標(biāo)準(zhǔn)正態(tài)分布樣本變換為原來的分布，再代入原來的極限狀態(tài)方程求解示性函數(shù).這種變換方法是一般變換方法的逆變換.

1.3.1 一般變換方法

任意相關(guān)隨機變量必須經(jīng)過2步轉(zhuǎn)化才能變成獨立標(biāo)準(zhǔn)正態(tài)分布.第一步是將相關(guān)非正態(tài)隨機變量轉(zhuǎn)化為相關(guān)正態(tài)隨機變量；第二步是將相關(guān)正態(tài)隨機變量轉(zhuǎn)化為獨立正態(tài)隨機變量.

將非正態(tài)隨機變量轉(zhuǎn)化為正態(tài)隨機變量最常用的方法是當(dāng)量正態(tài)化方法和映射變換方法.通過當(dāng)量正態(tài)化方法進行變換后的相關(guān)系數(shù)可以根據(jù)NATAF變換得到［11］.變換后的相關(guān)系數(shù)表示為

式中系數(shù)F已經(jīng)由Liu給出［12］.

相關(guān)正態(tài)隨機變量轉(zhuǎn)化為獨立的標(biāo)準(zhǔn)隨機變量可以采用Choleshy分解方法進行.通過這些方法就可以將任意相關(guān)隨機變量轉(zhuǎn)化為獨立的標(biāo)準(zhǔn)正態(tài)分布的隨機變量.變換流程圖見圖1.上述方法的逆變換也是成立的.文中將上述方法的逆變換用于結(jié)構(gòu)可靠性分析的β球方法中.

圖1 變換過程

1.3.2 逆變換的驗證

有2個隨機變量R和S，R服從正態(tài)分布，S服從極值I型分布，相關(guān)系數(shù)為ρRS＝0.5，其中μR＝100，μS＝50，δR＝0.12，δS＝0.15.采用 Monte－Carlo方法，使用標(biāo)準(zhǔn)正態(tài)分布的隨機數(shù)經(jīng)過逆變換來得到R和S的分布以及相關(guān)系數(shù).樣本數(shù)為40 000時的模擬結(jié)果見圖2和圖3.模擬得到的相關(guān)系數(shù)為0.491 3.可以看出本文采用的逆變換是準(zhǔn)確可靠的.

圖2 隨機變量R的模擬結(jié)果

圖3 隨機變量S的模擬結(jié)果

1.4 擴展的一般β球方法的計算過程

擴展的一般β球方法的計算流程圖見圖4.圖4中Yi為標(biāo)準(zhǔn)正態(tài)空間的一組隨機變量；Xi為一組非正態(tài)分布隨機變量；G（X）為極限狀態(tài)方程.

圖4 擴展的一般β球方法計算流程圖

2 自適應(yīng)抽樣算法的提出

β球方法的使用過程中需要對各隨機變量進行抽樣，根據(jù)抽樣計算結(jié)構(gòu)的失效概率.可以利用抽樣結(jié)果估算結(jié)構(gòu)可靠性指標(biāo)，以這個估算的可靠性指標(biāo)作為初始β球的大小，然后再根據(jù)擴展的一般β球方法求解結(jié)構(gòu)的失效概率.

2.1 自適應(yīng)抽樣算法

以包含兩個標(biāo)準(zhǔn)正態(tài)分布隨機變量的極限狀態(tài)方程為例進行說明所提出的自適應(yīng)抽樣算法，對于不服從正態(tài)分布的情況按照1.3節(jié)的變換方法變換為標(biāo)準(zhǔn)正態(tài)分布.算法的原理如圖5所示.具體的步驟描述為［13］：

步驟1假設(shè)抽樣范圍為β0和β1之間，β0取得相對較小一點，β1取得相對較大一點.

步驟2抽取樣本，并代入極限狀態(tài)方程中判斷極限狀態(tài)方程G（X）是否小于0，如果小于0，計為P1，停止抽樣.

步驟3求取P1點處的方向余弦，P1點的方向余弦和P2點的方向余弦相等，進行變換，根據(jù)極限狀態(tài)方程采用二分法求解出β2.

步驟4以β0和β2為抽樣區(qū)間重復(fù)步驟2，步驟3，在P2點根據(jù)重要抽樣策略進行抽樣，求出β3.

步驟5重復(fù)步驟2～4，3～5次就可以得到合適的β值.

圖5 自適應(yīng)抽樣算法原理示意圖

2.2 自適應(yīng)抽樣算法的驗證

2.2.1 獨立標(biāo)準(zhǔn)正態(tài)分布

以文獻［14］中的一個例題為例進行可靠性指標(biāo)的估算.已知非線性極限狀態(tài)方程0.1（x1－式中：x1，x2都服從標(biāo)準(zhǔn)正態(tài)分布.假設(shè)β在1和6之間，經(jīng)過5次計算，得到的結(jié)果見圖6.由圖6可見兩者十分接近.為了進一步驗證計算結(jié)果，極限狀態(tài)函數(shù)和抽樣區(qū)域見圖7.

2.2.2 獨立非標(biāo)準(zhǔn)正態(tài)分布

圖6 正態(tài)分布的估算結(jié)果

圖7 計算結(jié)果的比較

以文獻［4］中的一個例題為例進行可靠性指標(biāo)的估算.已知非線性極限狀態(tài)方程567fr－0.5H2＝0.f服從正態(tài)分布，μf＝0.6，δf＝0.131；r服從正態(tài)分布，μr＝2.18，δr＝0.03；H服從對數(shù)正態(tài)分布，μH＝32.8，δH＝0.03.假設(shè)β在1～6之間，經(jīng)過6次計算，得到的結(jié)果如圖8所示.圖中虛線為驗算點法的解，可以看出兩者十分接近.

圖8 非正態(tài)分布的估算結(jié)果

2.2.3 相關(guān)非正態(tài)分布

以文獻［4］中的一個例題為例進行可靠性指標(biāo)的估算.已知極限狀態(tài)方程為X1X2－130＝0.X1服從對數(shù)正態(tài)分布，μX1＝38.0，δX1＝0.1；X2服從正態(tài)分布，μX2＝7.0，δX2＝0.15；相關(guān)系數(shù)ρX1X2＝0.5.假設(shè)β在1和6之間，經(jīng)過6次計算，得到的結(jié)果和驗算點法的比較見圖9.

3 自適應(yīng)β球方法的提出

3.1 自適應(yīng)β球方法

圖9 相關(guān)非正態(tài)分布的估算結(jié)果

對于一個含有任意個非正態(tài)分布隨機變量的極限狀態(tài)方程，首先采用自適應(yīng)抽樣算法估算出初始β球的大小，再根據(jù)擴展的一般β球方法來計算結(jié)構(gòu)的可靠性指標(biāo).

3.2 計算實例

我國現(xiàn)行的《潛水系統(tǒng)和潛水器入級與建造規(guī)范》，給出了耐壓環(huán)肋圓柱殼結(jié)構(gòu)的5種校核公式.大量的潛水器結(jié)構(gòu)可靠性研究文獻將這5種校核公式作為失效模式的極限狀態(tài)方程.具體說明以及隨機變量的分布和特征參數(shù)見文獻［10］.對于潛水器這樣安全要求比較高的結(jié)構(gòu)，一般認(rèn)為只要有一個失效模式出現(xiàn)就認(rèn)為整個結(jié)構(gòu)失效，所以潛水器結(jié)構(gòu)系統(tǒng)可靠性是一個串聯(lián)模式的可靠性問題.

3.2.1 單個失效模式的可靠性計算

采用提出的自適應(yīng)β球方法對五種失效模式的可靠性分別進行了計算，全部的計算結(jié)果見表1，其中G3的計算結(jié)果和驗算點法的比較見圖10.

表1 計算結(jié)果

圖10 G3的計算結(jié)果

3.2.2 系統(tǒng)可靠性計算

采用自適應(yīng)抽樣算法計算每一個失效模式的可靠性指標(biāo)βi（i＝1，2，3，4，5），然后取min（βi）作為β球的大小，采用擴展的一般β球法進行計算，示性函數(shù)I［g（X）］表達(dá)為

式中：I［gXi（X）］為第i個失效狀態(tài)的示性函數(shù).系統(tǒng)可靠性的計算結(jié)果如圖11所示.

圖11 系統(tǒng)可靠性計算結(jié)果

將本文的計算結(jié)果和重要樣本法（important sample method）、PCM 法［15］以及 Drezner積分方法（Drezner integral method）［16］的結(jié)果進行了比較.結(jié)果表明本文的計算結(jié)果和Drezner積分方法以及PCM法接近.Drezner積分方法是一種特殊的Gauss－Hermite求積法，在理論上只存在積分誤差.PCM方法的計算結(jié)果比本文方法略小，文獻［17］指出PCM方法在計算串聯(lián)結(jié)構(gòu)系統(tǒng)可靠性時往往給出較保守的計算結(jié)果，這和本文的計算結(jié)果基本一致.本文方法和重要樣本法有一定的差異，這主要是因為重要樣本法收斂較慢，在相同的樣本大?。╯ample size）下，本文方法已收斂，而重要樣本法還沒有收斂.這主要歸功于自適應(yīng)抽樣算法提供了一個合理的初始β球尺寸.

4 結(jié)束語

Harbitz提出的β球法是一種應(yīng)用廣泛且高效的數(shù)值計算方法，但這種方法不能獨立使用.針對這一點，本文提出了一種自適應(yīng)抽樣算法.這種算法可以根據(jù)β球方法的抽樣特點估計出β球的大小，算法的效率較高，經(jīng)過5～6步的迭代就可以得出和驗算點法相近的計算結(jié)果.包含自適應(yīng)抽樣算法的自適應(yīng)β球方法收斂速度高于重要樣本法，且計算準(zhǔn)確性較高.這種自適應(yīng)β球方法特別適用于極限狀態(tài)函數(shù)是凸曲面以及結(jié)構(gòu)系統(tǒng)可靠性的模擬，對一般的極限狀態(tài)方程也是有效的.

本文最后使用提出的方法計算了潛水器耐壓圓柱殼結(jié)構(gòu)各個失效模式的可靠性和系統(tǒng)可靠性，并和驗算點法、重要抽樣法、PCM方法以及Drezner積分方法進行了比較，驗證了本文提出的方法的計算效率和準(zhǔn)確性.

［1］貢金鑫.工程結(jié)構(gòu)可靠度計算方法［M］.大連：大連理工大學(xué)出版社，2003.

［2］肖 剛，李天柁.系統(tǒng)可靠性分析中的蒙特卡羅方法［M］.北京：科學(xué)出版社，2003.

［3］謝祚水，王自力，吳劍國.潛艇結(jié)構(gòu)分析［M］.武漢：華中科技大學(xué)出版社，2002.

［4］趙國藩，金偉良，貢金鑫.結(jié)構(gòu)可靠度理論［M］.北京：中國建筑工業(yè)出版社，2000.

［5］Harbitz A.An efficient sampling method for probability of failure calculation［J］.Structural Safety，1986（3）：109－115.

［6］余建星，郭振邦，徐 慧，等.船舶與海洋結(jié)構(gòu)物可靠性原理［M］.天津：天津大學(xué)出版社，2001.

［7］孫海虹.結(jié)構(gòu)可靠性分析改進的重要抽樣法［J］.武漢交通科技大學(xué)學(xué)報，1994，18（3）：241－246.

［8］張曉軍，常新龍，楊 青.利用改進的Monte Carlo β球法計算結(jié)構(gòu)可靠性［J］.強度與環(huán)境，2005，32（2）：52－56.

［9］吳亞舸，吳劍國.潛艇耐壓圓柱殼結(jié)構(gòu)可靠性的數(shù)值模擬計算法［J］.造船技術(shù)，2002（1）：15－18.

［10］馬永亮.考慮腐蝕影響的潛艇結(jié)構(gòu)可靠性研究［D］.哈爾濱：哈爾濱工程大學(xué)船舶工程學(xué)院，2009.

［11］Melchers R E.Structural reliability analysis and prediction（Second Edition）［M］.New York：John Wiley ＆Sons Ltd.，1999.

［12］Liu Peilin，Kiureghian A D.Multivariate distribution models with prescribed marginals and covariances［J］，Prob.Engineering Mechanics，1986（2）：105－112.

［13］Nie Jinsuo.A new directional method to assess structural system reliability in the context of performance－based design［D］.Baltimore：Johns Hopkins University，2003.

［14］Borri A，Speranzini E.Structural reliability analysis using a standard deterministic finite element code［J］.Structural Safety，1997，19（4）：361－82.

［15］Pandey M D.An effective approximation to evaluate multi－normal integrals［J］.Structural Safety，1998，20：51－67.

［16］Drezner Z.Computation of the multivariate normal integral［J］.ACM Transactions on Mathematical Software，1992，18（4）：470－480.

［17］Mori Y，Kato T.Multi－normal integrals by importance sampling for series system reliability［J］.Structural Safety，2003（25）：363－378.