潛水器結構可靠性計算的自適應β球方法＊

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

（哈爾濱工程大學多體船技術國防重點學科實驗室 哈爾濱 150001）

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

1 一般β球方法

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

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

β球將空間分為了2個部分｜X｜≤β和｜X｜＞β，其中｜X｜≤β為中心在原點的m維超球體.根據全概率公式，失效概率Pf可以表示為

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

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

式中：k為樣本的大?。籌［G（Yi）］為示性函數，在安全域I＝0，在失效域I＝1；Yi為標準正態空間的一組隨機變量.

1.2 β球方法抽樣策略

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

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

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

1.3 樣本變換及其改進

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

1.3.1 一般變換方法

任意相關隨機變量必須經過2步轉化才能變成獨立標準正態分布.第一步是將相關非正態隨機變量轉化為相關正態隨機變量；第二步是將相關正態隨機變量轉化為獨立正態隨機變量.

將非正態隨機變量轉化為正態隨機變量最常用的方法是當量正態化方法和映射變換方法.通過當量正態化方法進行變換后的相關系數可以根據NATAF變換得到［11］.變換后的相關系數表示為

式中系數F已經由Liu給出［12］.

相關正態隨機變量轉化為獨立的標準隨機變量可以采用Choleshy分解方法進行.通過這些方法就可以將任意相關隨機變量轉化為獨立的標準正態分布的隨機變量.變換流程圖見圖1.上述方法的逆變換也是成立的.文中將上述方法的逆變換用于結構可靠性分析的β球方法中.

圖1 變換過程

1.3.2 逆變換的驗證

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

圖2 隨機變量R的模擬結果

圖3 隨機變量S的模擬結果

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

擴展的一般β球方法的計算流程圖見圖4.圖4中Yi為標準正態空間的一組隨機變量；Xi為一組非正態分布隨機變量；G（X）為極限狀態方程.

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

2 自適應抽樣算法的提出

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

2.1 自適應抽樣算法

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

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

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

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

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

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

圖5 自適應抽樣算法原理示意圖

2.2 自適應抽樣算法的驗證

2.2.1 獨立標準正態分布

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

2.2.2 獨立非標準正態分布

圖6 正態分布的估算結果

圖7 計算結果的比較

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

圖8 非正態分布的估算結果

2.2.3 相關非正態分布

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

3 自適應β球方法的提出

3.1 自適應β球方法

圖9 相關非正態分布的估算結果

對于一個含有任意個非正態分布隨機變量的極限狀態方程，首先采用自適應抽樣算法估算出初始β球的大小，再根據擴展的一般β球方法來計算結構的可靠性指標.

3.2 計算實例

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

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

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

表1 計算結果

圖10 G3的計算結果

3.2.2 系統可靠性計算

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

式中：I［gXi（X）］為第i個失效狀態的示性函數.系統可靠性的計算結果如圖11所示.

圖11 系統可靠性計算結果

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

4 結束語

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

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

［1］貢金鑫.工程結構可靠度計算方法［M］.大連：大連理工大學出版社，2003.

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

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

［4］趙國藩，金偉良，貢金鑫.結構可靠度理論［M］.北京：中國建筑工業出版社，2000.

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

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

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

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

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

［10］馬永亮.考慮腐蝕影響的潛艇結構可靠性研究［D］.哈爾濱：哈爾濱工程大學船舶工程學院，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.