潛水器結(jié)構(gòu)可靠性計(jì)算的自適應(yīng)β球方法＊

馬永亮 曲先強(qiáng) 崔洪斌 石德新

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

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

1 一般β球方法

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

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

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

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

失效概率Pf的一個(gè)無偏估計(jì)表示為

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

1.2 β球方法抽樣策略

基本隨機(jī)向量X抽樣可以表示為

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

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

1.3 樣本變換及其改進(jìn)

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

1.3.1 一般變換方法

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

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

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

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

圖1 變換過程

1.3.2 逆變換的驗(yàn)證

有2個(gè)隨機(jī)變量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)分布的隨機(jī)數(shù)經(jīng)過逆變換來得到R和S的分布以及相關(guān)系數(shù).樣本數(shù)為40 000時(shí)的模擬結(jié)果見圖2和圖3.模擬得到的相關(guān)系數(shù)為0.491 3.可以看出本文采用的逆變換是準(zhǔn)確可靠的.

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

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

1.4 擴(kuò)展的一般β球方法的計(jì)算過程

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

圖4 擴(kuò)展的一般β球方法計(jì)算流程圖

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

以文獻(xiàn)［4］中的一個(gè)例題為例進(jìn)行可靠性指標(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次計(jì)算，得到的結(jié)果如圖8所示.圖中虛線為驗(yàn)算點(diǎn)法的解，可以看出兩者十分接近.

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

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

以文獻(xiàn)［4］中的一個(gè)例題為例進(jìn)行可靠性指標(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次計(jì)算，得到的結(jié)果和驗(yàn)算點(diǎn)法的比較見圖9.

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

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

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

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

3.2 計(jì)算實(shí)例

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

3.2.1 單個(gè)失效模式的可靠性計(jì)算

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

表1 計(jì)算結(jié)果

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

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

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

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

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

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

4 結(jié)束語

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

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

［1］貢金鑫.工程結(jié)構(gòu)可靠度計(jì)算方法［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ìn)的重要抽樣法［J］.武漢交通科技大學(xué)學(xué)報(bào)，1994，18（3）：241－246.

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

［9］吳亞舸，吳劍國.潛艇耐壓圓柱殼結(jié)構(gòu)可靠性的數(shù)值模擬計(jì)算法［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.