基于高效求矩法與最大熵原理的矩方法

姜昊銳，趙衍剛

(北京工業大學城市建設學部，北京 100124)

工程結構在使用過程中存在著諸如材料參數、幾何參數以及外部荷載作用等不確定性，日益引發了人們對基于概率設計的重視，結構可靠度分析在結構設計與使用過程中起到越來越重要的作用[1]。目前常用的可靠度分析方法有蒙特卡洛模擬法[2]、一次二階矩法(first order reliability method, FORM)[3-5]與二次二階矩法(second order reliability method, SORM)[6-8]等,但蒙特卡洛模擬方法在處理小失效概率問題時需要大量的計算,效率低下；FORM與SORM因為需要計算導數以及迭代運算，在處理復雜工程問題時精度較低。矩方法也是可靠性分析方法之一，因其方法簡單，便于操作而受到越來越多的重視。矩方法具體包含兩個步驟：首先求解功能函數的統計矩，如前四階中心矩：均值，標準差，偏度，峰度；接著依據求得的統計矩信息，選擇合適的分布模型近似真實的概率分布，最終在失效域上進行積分求解結構的失效概率與可靠指標。由上述步驟可知，矩方法不需要像蒙特卡洛模擬法進行大量運算，也無需像FORM或SORM一樣計算導數以及迭代運算，但求矩的效率與精度直接影響著矩方法的精度與效率[7-9]，如何高效且準確地求解功能函數的矩,目前仍舊是矩方法所研究的重點。

求解功能函數的統計矩，實質上就是計算多維積分。在實際工程中，由于結構的功能函數通常較為復雜且以隱式函數表示，直接進行多維積分運算往往較為困難[10]。因此，許多近似計算方法相繼被提出。Zhao等[11]提出了新點估計方法，通過Rosenblatt變換[12]和一維減維方法[13]將功能函數表示成一元函數和的形式，繼而用一維點估計進行求矩運算，很大程度上提高了計算的效率，但該方法在功能函數較為復雜的情況下計算精度較低。為了提高計算精度，Zhao等[14]提出了將Rosenblatt變換和二維減維方法[15]相結合的二維減維點估計方法，將功能函數表示成一元函數與二元函數和的形式進行求矩運算，但隨著輸入隨機變量數目的增加，該方法往往會造成“維數災難”問題。Xu等[16]提出了用高階無跡變換[17-19]計算二維減維點估計中的二維積分問題，該求矩方法可一定程度上提高計算精度與效率。為了進一步提高計算的效率，Fan等[20]提出了自適應二維減維點估計方法，通過引入交叉項判定準則，判斷二維減維后的二元函數中是否存在交叉項，若不存在交叉項，則將二元函數進一步用一元函數表示，并用一維點估計法進行求矩運算；若存在交叉項，則用二維點估計進行計算，但在功能函數所含交叉項較多的情況下，該方法依舊需要大量的運算。針對此問題，Cai等[21]提出采用二維稀疏網格法[22]計算自適應二維減維點估計方法中的二維積分，可一定程度上提高計算效率，但在輸入隨機變量數目較多的情況下計算量依舊較大。綜上所述，已有的求矩方法均不能在計算過程中兼顧精度與效率。

鑒于此，現提出一種基于快速求矩法與最大熵原理的矩方法，將自適應二維減維與高階無跡變換相結合求解結構功能函數的前四階矩，并選擇最大熵原理[11]近似真實的概率分布，最終依據最大熵分布求解結構的可靠指標。算例部分通過兩道例題，對比所提方法與幾種已有方法的計算結果，驗證所提方法在求矩運算以及可靠指標計算問題中在效率以及精度上的優勢。

1 基于矩方法的可靠度分析

1.1 問題表述

結構的功能函數[5]的表達式通常為

Z=G(X)

(1)

式(1)中：Z為功能函數；X=[X1,X2,…,Xn]為n維隨機向量，功能函數的前四階中心矩計算公式為

μG=μ1G

(2)

(3)

(4)

(5)

式中：μG、σG、α3G和α4G分別為功能函數的均值、標準差、偏度與峰度；μkG為功能函數的k階原點矩(k=1,2,3,4),具體計算表達式為

(6)

式(6)中：fx(x)為輸入隨機向量的聯合概率密度函數。由式(6)可知，求矩運算本質上是多維積分計算問題，但實際工程中功能函數較為復雜，且隨機向量的聯合概率密度函數難以直接獲得，因此直接求解并不可行。為了高效且精確地求矩，提出基于自適應二維減維與高階無跡變換的高效求矩法。

1.2 自適應二維減維

結構可靠性分析一般在標準正態空間下進行，因此可通過Rosenblatt變換或者Nataf變換[23]，將輸入的具有相關性的非正態隨機變量轉換為獨立的標準正態隨機變量U=[U1,U2,…,Un]，此時功能函數可用標準正態隨機變量[18]表示為

Z=G(X)=G[T-1(U)]=L(U)

(7)

式(7)中：T-1為Rosenblatt逆變換或Nataf逆變換；L為G與T-1的復合函數。

基于二維減維[12]方法，功能函數可進一步寫為

(8)

式(8)中：L2為n(n-1)/2個二元函數的和；L1為n個一元函數的和；L0為功能函數在均值點處的函數值，為一常數。L2、L1、L0的具體表達式為

(9)

(10)

L0=L(0,…,0,…,0)

(11)

式中：i,j=1,2,…,n,i

基于式(8)～式(11)可得，功能函數在二維減維后的k階原點矩的近似計算表達式為

(12)

(13)

(14)

式中：φ(u)為標準正態隨機變量的概率密度函數。由式(14)易知求解μk-Li為一維積分問題，可通過一維點估計法[15]進行計算，即

(15)

式(15)中：ur為積分點；Pr為對應權重；N為積分點個數。一維點估計中積分點及其對應權重的計算式[16]為

(16)

(17)

式中：xr、ωr分別為以exp(-x2)為權函數的2N-1階高斯埃爾米特積分的積分點與權重值[24]。

式(13)中含有兩個隨機變量，為二維積分問題，文獻[21]采用二維點估計進行運算，即

(18)

針對此情況，自適應二維減維點估計方法[17]提出了交叉項判定準則，通過判定準則判斷函數中是否存在交叉項，將判定為不含有交叉項的二元函數進一步用一元函數表示，對應的矩的計算可直接利用一維積分的計算結果；判定為存在交叉項的二元函數，其矩仍通過二維點估計進行計算。

判斷二元函數中是否存在交叉項的具體步驟如下。

(1)首先選取參考點Ur=(uir，uir),并帶入功能函數計算，用Lij(Ur)表示。

(2)另選取3個參考點，分別記為Ui-1=(ui-1，ujr),Uj-1=(uir，uj-1)以及Uij-1=(ui-1，uj-1)，代入式(19)計算Δ1，即

(19)

若Δ1>ε，則可判定在二元函數Lij(Ui，Uj)中Ui，Uj存在交叉項，否則執行下一步。

(3)再選取3個計算點，分別記為Ui-2=(ui-2，ujr),Uj-2=(uir，uj-2)以及Uij-2=(ui-2，uj-2)，代入式(20)計算Δ2，即

(20)

若Δ2>ε，則可判定在二維函數Lij(Ui，Uj)中Ui，Uj存在交叉項，否則Lij(Ui，Uj)中Ui，Uj不存在交叉項。

在上述步驟中ε為一特定值，作為判別交叉項存在的標準，取ε=10-6。值得注意的是，交叉項判定選取的計算點通常是進行一維點估計計算的積分點，因此不需額外增加計算量。

如果二維函數Lij(Ui，Uj)經上述判別準則判定為不存在交叉項時，則二元函數可用一元函數表示為

Lij(Ui,Uj)=Li(Ui)+Lj(Uj)-L0

(21)

繼而不含交叉項的二元函數矩的求解，可直接利用一元函數Li(Ui)、Lj(Uj)矩的計算結果進行計算，即

μ1-Lij=μ1-Li+μ1-Lj-L0

(22)

(23)

(24)

(25)

式中：M2ij、M3ij、M4ij計算表達式分別為

(26)

(27)

(28)

由式(22)～式(28)可知，此時二維積分計算可直接利用一維積分計算結果，無需額外增加計算量，因此一定程度上提高了求矩的效率。

若判定二元函數中存在交叉項，則采用二維點估計進行計算，為保證精度通常采用5點或7點估計，此時每個二元函數的求矩過程分別需計算25次或49次功能函數值。在交叉項較多的情況下，進行二維點估計求矩的二元函數數目較多，因此仍舊需要大量的運算。針對此情況，通過采用高階無跡變換方法替換原有的二維點估計，減少了二維積分中積分點的個數，相應減少了功能函數的計算次數，一定程度上提高了求矩的計算效率。

1.3 基于高階無跡變換的二維積分算法

高階無跡變換可以用于解決n維高斯權重積分問題[21]。高階無跡變換的基本思想是[18]：首先，選取一些特定的積分點(Sigma點)及其對應的權重來匹配輸入隨機變量的概率信息，通常為輸入隨機變量的各階統計矩；其次，將選取的積分點通過非線性變換，即帶入結構的功能函數計算對應的輸出狀態變量的樣本點；最后，輸出狀態變量的統計矩可以利用輸入變量的積分點對應權重和輸出狀態變量的樣本點進行估計，其實質上為一種改進的高斯積分。

由前述可知，高階無跡變換積分點總個數為N=2n2+1，在自適應二維減維后的二維積分問題中，n=2，則每個二元函數矩的計算僅需9次運算，其效率明顯高于二維點估計中的5點或7點估計。

高階無跡變換三類積分點及其對應權重的具體表達式[13,15]如下。

第一類：

θ0=0

(29)

(30)

第二類：

(31)

(32)

(33)

式中：ej1=[0,…,1,…,0]，即只有第j1項為1。

第三類：

(34)

(35)

(36)

(37)

(38)

在計算第二類和第三類積分點時，需要引入自由參數β，其影響著高階無跡變換的計算結果的精度，文獻[16,18]中對β的取值進行過研究，在此取β=7.2，經不同例題試算結果表明，所選取的β值計算結果精度較高。

當自由參數確定后，三類積分點及權重由式(29)～式(38)進行計算，此時含有交叉項的二元函數的矩可直接由式(39)進行計算，即

(39)

式(39)中：ωi、θi分別為高階無跡變換積分點與對應權重，i=1,2,…,n。

最后，將計算所得的自適應二維減維后的一元函數、二元函數的矩代入式(12)中，功能函數的前四階原點矩μkG即可確定(k=1,2,3,4),繼而可通過式(2)～式(5)計算功能函數G(X)的前四階中心矩。

獲得功能函數的前四階中心矩后，需選擇合適的分布模型對功能函數的真實概率分布進行近似，并以此計算結構可靠指標。常用的分布模型有Pearson分布族、Burr分布族等，盡管該類分布系統非常靈活，但在不同分布類型的分界處往往難以實現分布近似。因此，采用最大熵原理近似結構功能函數的概率分布，避免了上述問題，同時該分布模型對真實概率分布的擬合效果較好，且分布參數計算方便。

1.4 基于最大熵原理的可靠指標求解

若功能函數Z服從概率密度函數為fZ(z)的連續分布，其熵定義[25-26]為

(40)

根據最大熵原理可知，最大熵分布為所有可能的分布形式中，在給定約束條件下,使熵最大的概率分布，其具有最小偏見，因此被認為是構造“最優”概率分布的一種途徑。

為使式(40)的熵取最大值，此處以計算的功能函數前四階原點矩為約束條件。在計算最大熵分布時，為求解方便，提高收斂速度，通常先將功能函數進行標準化處理為

(41)

此時對功能函數Z的分布近似轉變為對標準化變量Y的概率分布的近似。

最大熵分布的求解可通過引入拉格朗日算法[26]進行計算，表達式為

(42)

令式(42)關于拉格朗日乘子的偏導數為0，即

(43)

(44)

利用所得到的最大熵分布，在失效域上積分可計算得到結構的失效概率為

(45)

對應的結構可靠度指標計算公式[3]為

β=Φ-1(1-pf)

(46)

式(46)中：Φ-1為標準高斯隨機變量的累計分布函數的逆函數。

綜上，以自適應二維減維與高階無跡變換求得的功能函數前四階矩為約束條件，結合最大熵原理，最終便可計算得到結構的可靠指標。

1.4 算法步驟

所提出的一種新的矩方法，基于自適應二維減維與高階無跡變換高效求解功能函數的前四階中心矩，并采用最大熵原理求解結構的可靠指標。具體計算步驟如下。

(1)根據式(8)將結構功能函數進行二維減維，表示成二元函數與一元函數和的形式。

(2)一元函數的矩直接通過式(15)的一維點估計進行運算求解。

(3)引入交叉項判定準則：①通過準則判定為不含有交叉項的二元函數，利用式(21)進一步用一元函數表示，其對應的矩可利用已求得的一元函數矩的結果直接進行計算；②判定為含有交叉項的二元函數，其矩采用高階無跡變換方法進行計算。

(4)將求得的一元函數、二元函數的矩的計算結果代入式(12)中，即可計算功能函數的前四階原點矩，繼而由式(2)～式(4)計算前四階中心矩。

(5)以計算得到的前四階中心矩為約束條件，基于最大熵原理求解功能函數的最大熵分布，繼而通過最大熵分布在失效域上的積分，可最終求解結構的失效概率與可靠指標。

2 算例分析

為驗證所提方法在計算精度與效率上的優勢，選取兩道算例，將所提高效求矩方法計算得到的功能函數前四階中心矩結果，與蒙特卡洛模擬方法、一維減維點估計方法、二維減維點估計方法、自適應二維減維點估計、二維減維無跡變換及自適應二維減維稀疏網格積分方法進行了對比。此外，將基于最大熵原理計算得到的可靠指標與蒙特卡洛方法結果進行了對比。

2.1 算例1：軸向壓力環形柱可靠度分析

某根受軸向壓力的環形柱，其受力圖與平面幾何尺寸如圖1所示，其中環形柱的彈性模量為E，長度為L，內徑為D，壁厚為T，受到軸向荷載P的作用。其中，E、L、D、T、P為相互獨立的隨機變量，具體對應的概率分布及分布參數如表1所示。

圖1 軸向壓力環形柱示意圖

表1 各隨機變量的分布信息

此處受軸向壓力環形柱的失效模式為屈曲破壞，其對應的極限狀態函數表達式為

(47)

表2中分別列出了蒙特卡洛模擬方法，一維減維點估計方法、二維減維點估計方法，自適應二維減維點估計，二維減維無跡變換及自適應二維減維稀疏網格積分方法的求矩結果與計算次數。同時以蒙特卡洛模擬法結果為標準，不同方法計算結果的相對誤差也在表內括號里列出。其中，相對誤差計算公式為

(48)

式(48)中：e為相對誤差；rMCS為蒙特卡洛模擬方法計算結果；r為各求矩方法計算結果。

由表2可知，一維減維點估計方法計算效率較高，但偏度與峰度的計算誤差較大。所提求矩方法相較于其他方法，計算效率有一定程度的提升，僅需計算59次，且提高效率的同時，計算的精度也能保證，且誤差更小。其中，本文方法計算得到的前四階中心矩中偏度的相對誤差最大，但僅有1.31%。

表2 各求矩方法計算結果

求解得到前四階中心矩后，便可按照1.4節步驟計算對應的最大熵分布，繼而進行可靠指標計算。表3中列出了本文方法結合最大熵分布計算得到的可靠指標結果，并與蒙特卡洛方法計算結果進行對比，以驗證所求結果的精度。同樣，采用式(12)計算本文方法計算得到的可靠指標相對于蒙特卡洛模擬方法結果的相對誤差。

表3 可靠指標計算結果

從表3可知，本文方法求得的可靠指標與蒙特卡洛較為吻合，在保留三位有效數字情況下結果與蒙特卡洛方法結果基本一致，相對誤差僅有0.7%。由此可知本文方法結合最大熵原理進行可靠指標計算時精度較高，符合誤差范圍。

2.2 算例2：空間剛架可靠度分析

如圖2所示，有一個由72根桿件構成的空間桁架，其中各個桿件的橫截面面積均為34.849 cm2,F1、F2、F3、F4、F5、F6、F7、F8分別為作用在節點5、8、9、12、13、16、17、20處的水平作用力，所有作用力與構成空間剛架的桿件的彈性模量E互為獨立隨機變量，對應的概率分布與分布參數如表4所示。

圖2 空間72桿框架示意圖

表4 各隨機變量的分布信息

此處考慮該空間桁架的極限狀態為所有節點的水平位移絕對值最大值超過界限值，該失效模式所對應的結構功能函數表達式為

Z=G(X)=0.05-max|μi(x)|,i=5,6,…,20

(49)

式(49)中：0.05為沿x方向水平位移限值，m。本例題計算過程中，空間桁架的各個節點水平位移通過MATLAB軟件進行有限元分析計算。表5中列出了各個求矩方法計算得到的功能函數前四階中心矩，同算例1，表中括號內也列出了各個求矩方法計算結果相較于蒙特卡洛模擬方法結果的相對誤差。

在計算得到功能函數前四階中心矩后，依據最大熵原理計算相應的結構可靠指標并與蒙特卡洛方法對比，具體結果如表6所示。

表6 可靠指標計算結果

由表5可知，在求解隱式功能函數(14)的前四階中心矩時，本文方法計算功能函數次數較少，僅需計算211次功能函數值，且效率提升的同時也能保證計算的精度，其中本文方法計算得到的前四階中心矩中峰度結果的相對誤差最大，但也僅有3.9%。

表5 各求矩方法計算結果

最后通過最大熵分布擬合功能函數的分布函數并進行結構可靠度分析，計算得到該空間桁架的結構可靠指標為2.934，蒙特卡洛模擬法計算結果為3.051，相對誤差為3.8%，可知本文方法在效率提升的同時能保證一定的精度。

3 結論

提出了一種新的矩方法，基于自適應二維減維與高階無跡變換進行高效求矩，并結合最大熵原理計算結構的可靠指標。通過算例部分的兩道例題，對比了本文方法與既有方法的求矩結果以及可靠指標計算結果，最終能得到以下的結論。

(1)求矩方面，所提高效求矩方法在運算過程中，計算結構功能函數的次數少于其他求矩方法，因此在計算效率方面優于既有的求矩方法。通過不同的求矩結果與蒙特卡洛模擬方法結果的對比，可知所提求矩方法在提高效率的同時，在計算精度上也有明顯優勢，誤差更小。

(2)可靠指標計算方面，通過結合高效求矩法與最大熵原理計算得到可靠指標，并與蒙特卡洛方法結果對比，可知所提矩方法計算的結構可靠指標與蒙特卡洛模擬方法的結果較為吻合，相對誤差滿足工程要求。