范文亮，袁 滿，劉潤宇，楊曉陽，李正良

(1.重慶大學(xué)土木工程學(xué)院，重慶400045；2.山地城鎮(zhèn)建設(shè)與新技術(shù)教育部重點實驗室(重慶大學(xué))，重慶400045)

點估計法[1?3]是隨機系統(tǒng)常用分析方法，其計算簡便有效且可直接與矩方法結(jié)合進行結(jié)構(gòu)可靠度分析[4?5]。點估計法中多維數(shù)值積分不可避免，往往遭遇維數(shù)災(zāi)難問題。高效的選點策略與降維近似是目前常用的緩解維數(shù)災(zāi)難的方法。基于稀疏點集的點估計法[6]是高效選點策略中具有相對較廣適用范圍與較高精度的方法，但對于交叉項影響較顯著的情形精度會有所減弱[6?7]。文獻[8]的比較研究表明將多維函數(shù)近似為一系列較低維數(shù)分量函數(shù)組合的降維近似模型在有效緩解維數(shù)災(zāi)難的同時可保持較好的精度。因此，近年來基于降維近似模型的統(tǒng)計矩估計引起了廣泛的關(guān)注。對于一個多變量函數(shù)，通常可有乘積降維近似[9?10]和加和降維近似[11]兩種策略。Rosenblueth[2]首先引入了單變量乘積降維近似模型計算了多變量函數(shù)的統(tǒng)計矩。多變量乘積降維近似中，由于具有相同變量的分量函數(shù)通過累乘組合在一起，其多維積分并不能轉(zhuǎn)化為低維積分，因此，除單變量乘積降維近似模型外，鮮有乘積降維近似模型用于統(tǒng)計矩估計；然而，隨著隨機系統(tǒng)維數(shù)的增加，單變量乘積降維近似模型并不能保證合理的精度。相比較而言，加和降維近似模型中，分量函數(shù)通過加和的方式相組合，其多維積分必定可轉(zhuǎn)化為一系列的低維積分，因此，加和降維近似模型在統(tǒng)計矩估計中應(yīng)用更為廣泛。Zhao等[12]首先將多變量函數(shù)的單變量降維近似模型用于統(tǒng)計矩估計；隨后Rahman 等[13]亦提出了基于單變量降維近似模型的統(tǒng)計矩估計方法，并對單變量降維近似模型進行了詳細的理論分析；李洪雙等[14]采用了在文獻[13]的基礎(chǔ)上引入了Nataf 變換，拓展了其適用范圍。由于僅涉及到多個單變量函數(shù)的概率積分，基于單變量降維近似模型的點估計法具有極高的計算效率，然而對于交叉項作用較復(fù)雜隨機系統(tǒng)，其精度亦不夠理想。為改善精度，Xu 等[15]進一步提出了針對矩函數(shù)進行降維多變量加和降維近似的統(tǒng)計矩估計方法；Huang 等[16]和Yu 等[17]將矩函數(shù)的雙變量降維近似分別與Rosenblatt 變換和Nataf 變換相結(jié)合，發(fā)展了較為實用的統(tǒng)計矩估計方法。相比于矩函數(shù)，原函數(shù)更為簡單，F(xiàn)an 等[18]提出了針對原函數(shù)采用多變量加和降維近似的直接積分法，在不增加計算量的基礎(chǔ)上，進一步提高了統(tǒng)計矩估計的精度。不難發(fā)現(xiàn)，隨著降維近似模型截斷維數(shù)的增加，模型的精度會逐步改善，但分量函數(shù)會逐步增多，進而導(dǎo)致計算量激增。因此，將截斷維數(shù)取為2是兼顧精度和效率的實用選擇[16?18]。然而，即使對于雙變量加和降維近似模型，一旦隨機系統(tǒng)的維數(shù)較高，雙變量分量函數(shù)將急劇增加，其計算效率并不總是很理想。事實上，雙變量分量函數(shù)的積分中，各積分點的權(quán)重系數(shù)是不同的，若允許權(quán)重較小的計算點的響應(yīng)存在一定的誤差，則可針對分量函數(shù)引入較為成熟的代理模型以改善計算效率。

為此，本文基于函數(shù)近似和數(shù)值積分中求積節(jié)點的特性，提出了Kriging 近似模型訓(xùn)練點集確定的“米”字型策略，并基于此建立了雙變量分量函數(shù)的Kriging 近似，將其與現(xiàn)有的基于雙變量降維近似模型的統(tǒng)計矩估計方法相結(jié)合，建立了更為高效的統(tǒng)計矩點估計法。

1 基于雙變量降維近似模型的統(tǒng)計矩點估計法

通常，基于雙變量降維近似模型的統(tǒng)計矩計算包括三個步驟：1)將分布類型不同的隨機變量X向獨立標(biāo)準(zhǔn)正態(tài)空間轉(zhuǎn)換，得到新的隨機向量U=(U1,U2,···,Un)T和響應(yīng)函數(shù)h(U)；2)利用雙變量降維近似模型將多維響應(yīng)函數(shù)分解為二維及以下分量函數(shù)的線性組合；3)基于雙變量降維近似模型進行響應(yīng)函數(shù)的統(tǒng)計矩計算。

1.1 變量的獨立化正態(tài)變換

隨機向量X中的變量概率信息(如分布類型和相關(guān)信息等)常常不同，在隨機分析和可靠度分析中常將隨機向量X向獨立標(biāo)準(zhǔn)正態(tài)空間轉(zhuǎn)換，得到新的隨機向量U=(U1,U2,···,Un)T。于是，響應(yīng)函數(shù)Z=g(X)可改寫為：

式中：X R為X中已知聯(lián)合概率分布的變量組成的向量；XN為剩余變量組成的向量，既包括獨立變量，亦包括已知邊緣密度函數(shù)和相關(guān)信息的變量；U={UN,UR}。N?1(·)、R?1(·)分別表示Nataf 變換[19]與Rosenblatt 變換[20]的逆變換。

1.2 基于雙變量降維近似模型的統(tǒng)計矩估計

由概率論可知，隨機響應(yīng)函數(shù)Z的b階統(tǒng)計矩矩MZ,b的表達式如下：

式中，當(dāng)b=1時，a=0，MZ,1表示響應(yīng)函數(shù)的均值；當(dāng)b>1時，a=MZ,1，MZ,b為響應(yīng)函數(shù)的b階中心矩。

根據(jù)降維近似的對象不同，基于雙變量降維近似模型的統(tǒng)計矩估計可以分為兩類：1)基于矩函數(shù)[h(U)?a]b的雙變量降維近似模型的統(tǒng)計矩估計；2)基于原函數(shù)h(U)的雙變量降維近似模型的統(tǒng)計矩估計。

1.3 基于矩函數(shù)雙變量降維近似模型的統(tǒng)計矩估計

若引入[h(U)?a]b的雙變量降維近似模型，并進行期望運算，可給出MZ,b的近似計算表達式如下[16]：

1.4 基于原函數(shù)雙變量降維近似模型的統(tǒng)計矩估計

若先引入h(U)的雙變量降維近似模型，然后再對[h(U)?a]b進行期望運算，可給出MZ,b的另一種近似計算表達式[18]：

為簡便，該方法記為N-2。

2 改進的雙變量降維近似統(tǒng)計矩估計

在雙變量降維近似統(tǒng)計矩估計中，無論是D-2方法還是N-2方法，其計算量主要體現(xiàn)于雙變量分量函數(shù)在積分節(jié)點處的函數(shù)值計算，對于隱式的g(X)，其函數(shù)值需通過結(jié)構(gòu)分析獲得。實用中，通常取r=7，則任一雙變量分量函數(shù)均涉及7×7=49次的結(jié)構(gòu)分析。若能在保持計算精度相當(dāng)?shù)那疤嵯聹p少原函數(shù)的調(diào)用次數(shù)，則可以有效改善計算效率，較好地兼顧精度和效率的平衡。

事實上，在雙變量分量函數(shù)的49個求積節(jié)點中，權(quán)系數(shù)的差異較為明顯，權(quán)系數(shù)較小節(jié)點的函數(shù)值即使存在著一定的誤差，總體計算結(jié)果的誤差仍是可以接受的。為此，本文將任一雙變量分量函數(shù)的49個求積節(jié)點分為計算節(jié)點和近似節(jié)點，其中計算節(jié)點的函數(shù)值通過函數(shù)調(diào)用或結(jié)構(gòu)分析獲得，近似節(jié)點的函數(shù)值則由基于計算節(jié)點確定的雙變量分量函數(shù)的近似模型計算得到，此處近似模型采用Kriging 模型。由于近似節(jié)點的函數(shù)值不再涉及原函數(shù)調(diào)用或結(jié)構(gòu)分析，將有效改善計算效率。

2.1 Kriging 近似模型

2.2 改進的雙變量降維近似統(tǒng)計矩估計

2.2.1雙變量分量函數(shù)的Kriging 模型

1)計算節(jié)點的選取策略

圖1(a)給出了以標(biāo)準(zhǔn)正態(tài)密度函數(shù)為權(quán)函數(shù)的二維高斯求積分公式的49個積分節(jié)點，記為集合DOE。理論上，定義域內(nèi)任意點皆可以作為Kriging 模型的訓(xùn)練點集，即計算節(jié)點。然而，統(tǒng)計矩估計中僅需要雙變量分量函數(shù)在49個積分節(jié)點的函數(shù)值，因此，計算節(jié)點可在此49個節(jié)點中選取，且應(yīng)遵循如下準(zhǔn)則：

圖1 選點方案示意圖Fig.1 Schematic diagram of point selection

①具有較大權(quán)系數(shù)的計算節(jié)點函數(shù)值誤差對整體積分的影響往往較大，因此計算節(jié)點優(yōu)先選擇權(quán)系數(shù)較大的積分節(jié)點，換言之，計算節(jié)點應(yīng)盡量靠近中心區(qū)域；

②函數(shù)近似中，內(nèi)插精度往往高于外插精度，因此計算節(jié)點應(yīng)包含一定的外圍節(jié)點；

③當(dāng)涉及對稱型分布的隨機變量，其參考點坐標(biāo)分量為0，二維積分節(jié)點中的軸點和一維積分節(jié)點重合，為充分利用這些節(jié)點，計算節(jié)點宜優(yōu)先選擇軸點。

2)雙變量分量函數(shù)的Kriging 近似

首先，針對集合DOEY中所有的計算節(jié)點u0，通過調(diào)用函數(shù)或結(jié)構(gòu)分析確定雙變量函數(shù)h(Ui,Uj,u ij,c)

2.2.2改進的雙變量降維近似統(tǒng)計矩估計

相應(yīng)于常規(guī)的基于雙變量降維近似模型的統(tǒng)計矩估計，結(jié)合雙變量分量函數(shù)Kriging 近似的改進統(tǒng)計矩估計亦分為兩類：基于矩函數(shù)雙變量降維和Kriging 近似的統(tǒng)計矩估計、基于原函數(shù)雙變量降維和Kriging 近似的統(tǒng)計矩估計。

基于上述考慮，本文給出的計算節(jié)點選取方案如圖1(b)所示。為簡便，可將計算節(jié)點的集合記為DOEY，其形狀呈“米”字形，故稱為“米”字形策略。

1)基于矩函數(shù)雙變量降維和Kriging 近似的統(tǒng)計矩估計

將式(14)得到的雙變量分量函數(shù)的Kriging 近似代入式(3)，可得到改進的基于矩函數(shù)雙變量降維近似模型的統(tǒng)計矩估計：

為簡便，該方法記為米D-2。

2)基于原函數(shù)雙變量降維和Kriging 近似的統(tǒng)計矩估計

將式(14)得到的雙變量分量函數(shù)的Kriging 近似函數(shù)代入式(4)中，可得到改進的基于原函數(shù)雙變量降維近似模型的統(tǒng)計矩估計，其表達式如下：

為簡便，該方法記為米N-2。

不難發(fā)現(xiàn)，與D-2、N-2方法相比，米D-2和米N-2方法亦充分考慮了交叉項的影響，只不過是真實雙變量分量函數(shù)被合適的Kriging 函數(shù)近似代替。

2.3 算法步驟

綜合上述過程，改進的雙變量降維近似統(tǒng)計矩估計實現(xiàn)過程如下：

1)根據(jù)已知概率分布信息，由式(1)將隨機向量X轉(zhuǎn)化為標(biāo)準(zhǔn)正態(tài)向量U。

2)確定各雙變量分量函數(shù)在計算節(jié)點處的函數(shù)值，并基于此由式(14)確定各雙變量分量函數(shù)的Kriging 近似。

3)確定各單變量分量函數(shù)在求積節(jié)點處的函數(shù)值，由各雙變量分量函數(shù)的Kriging 近似確定其在二維積分節(jié)點的函數(shù)值，將上述函數(shù)值以及參考點處的函數(shù)值代入式(15)或式(16)即可得到響應(yīng)Z的b階統(tǒng)計矩。

3 算例分析

式中：X1為對數(shù)正態(tài)分布變量，均值為22、標(biāo)準(zhǔn)差為2；X2為正態(tài)分布變量，均值為10、標(biāo)準(zhǔn)差為0.9；X3為極值I型變量，均值為2、標(biāo)準(zhǔn)差為0.6。

對于上述三變量函數(shù)，參考點為uc=T?1(μX)=

圖2 h(U1,U2,12,c)的選點方案和Kriging 近似模型的精度Fig.2 Scheme for selecting points and accuracy of Kriging approximation model of h(U1,U2,12,c)

表1 統(tǒng)計矩計算結(jié)果對比Table 1 Comparison of statistical moments

算例2.工程算例

考察圖3所示承受軸向壓縮載荷F的環(huán)形截面柱，其功能函數(shù)為[22]：

式中：E、F、L、D、T分別是材料彈性模量、軸荷載、柱高、內(nèi)徑和厚度；E、F為對數(shù)正態(tài)變量，均值為2.1×1011Pa 和6000 N，標(biāo)準(zhǔn)差為0.4×1011Pa和400 N；D、T、L是獨立的正態(tài)變量，均值分別為0.03、0.006和2.5，標(biāo)準(zhǔn)差分別為0.001、0.0004和0.09。

圖3 環(huán)形柱模型Fig.3 Model of an annular section column

分別采用D-2、N-2和米D-2、米N-2四種方法計算式(20)所示功能函數(shù)的前四階矩，標(biāo)準(zhǔn)解為多維Gauss-Hermite數(shù)值積分迭代收斂的結(jié)果，結(jié)果如表2所示。從表2可以發(fā)現(xiàn)：1)高階矩計算結(jié)果中，D-2較N-2精度偏低，米D-2較米N-2精度偏低，米D-2和米N-2分別較D-2、N-2略高，其原因在于建議兩方法中降維近似誤差和Kriging 近似誤差的相互抵消；2)函數(shù)調(diào)用次數(shù)上，D-2和N-2相同，米D-2和米N-2相同，且由于引入了Kriging 近似后兩者的計算效率提升1倍多，效果顯著。

算例3.工程算例

式中：MUF、UTS、δ、N、R、R0分別是材料利用率、極限拉伸強度、密度、轉(zhuǎn)子轉(zhuǎn)速、外半徑和內(nèi)半徑；上述各變量均為獨立正態(tài)分布變量，其均值分別0.8358、22 000、0.29、21 000、24和8，標(biāo)準(zhǔn)差分別為0.2、5000、0.005、1000、0.5和0.3。

分別采用D-2、N-2和米D-2、米N-2四種方法計算式(21)所示功能函數(shù)的前四階矩，標(biāo)準(zhǔn)解為3×107個樣本點的蒙特卡羅法的計算結(jié)果，結(jié)果如表3所示。從表3可以發(fā)現(xiàn)：1)D-2和N-2兩者精度相當(dāng)，米D-2較米N-2略低，且前兩者在二、三階矩精度較后兩者略好，但在四階矩精度略差，因此，建議方法中的降維近似誤差和Kriging近似誤差的疊加效應(yīng)和抵消效應(yīng)對于各階矩并不總是保持一致；2)建議方法的函數(shù)調(diào)用次數(shù)從現(xiàn)有方法的577次降為217，計算效率顯著提高。

表2 統(tǒng)計矩計算結(jié)果對比Table 2 Comparison of statistical moments

表3 統(tǒng)計矩矩計算結(jié)果對比Table 3 Comparison of statistical moments

進一步對比上述3個算例，隨機系統(tǒng)維數(shù)從3上升至6時，現(xiàn)有方法的函數(shù)調(diào)用次數(shù)從154上升至577，而本文建議方法則從82上升至217。事實上，由算例1的效率分析可以發(fā)現(xiàn)基于雙變量降維近似模型的統(tǒng)計矩估計方法的函數(shù)調(diào)用次數(shù)與參考點坐標(biāo)有著較大的關(guān)系?？疾靚變量函數(shù)，對于參考點坐標(biāo)均為0和參考點坐標(biāo)均不為0兩種極端情況，N-2(或D-2)與米N-2(或米D-2)的函數(shù)調(diào)用次數(shù)表達式為：

式中：N1、N2分別為參考點坐標(biāo)均不為0時N-2與米N-2的函數(shù)調(diào)用次數(shù)；N3、N4分別為參考點坐標(biāo)均為0時N-2與米N-2的函數(shù)調(diào)用次數(shù)。圖4給出了2 種情況下各方法函數(shù)調(diào)用次數(shù)隨維數(shù)的變化規(guī)律，雖然米N-2的函數(shù)調(diào)用次數(shù)與N-2具有類似的增長規(guī)律，但總的分析次數(shù)顯著減少。若定義：

顯然，γ1和γ2分別反映了兩種極端情況下米N-2 與N-2的函數(shù)調(diào)用次數(shù)的比值。圖5給出了γ1和γ2隨維數(shù)的變化規(guī)律，不難發(fā)現(xiàn)γ1約為52%，γ2約為35%。

圖4 功能函數(shù)調(diào)用次數(shù)對比圖Fig.4 Comparison of structural analysis times

易知，對于參考點坐標(biāo)不全為0的一般情況，米N-2與N-2的函數(shù)調(diào)用次數(shù)的比值應(yīng)在35%～52%。換言之，相比于已有方法，文中建議方法可提升1倍的計算效率，能更好地兼顧精度和效率。

圖5 效率相對比值圖Fig.5 Comparison of efficiency ratios

4 結(jié)論

為了進一步提高基于雙變量降維近似模型的統(tǒng)計矩估計方法的效率，本文發(fā)展了基于“米”字形選點策略的雙變量函數(shù)的Kriging 近似模型，并結(jié)合現(xiàn)有的矩函數(shù)或原函數(shù)雙變量降維近似模型，提出了更好地兼顧精度和效率的改進統(tǒng)計矩估計方法。文中分析結(jié)果表明：

(1)本文建議方法(米N-2和米D-2)的計算效率相比于現(xiàn)有方法(N-2和D-2)顯著提高，函數(shù)調(diào)用次數(shù)約為現(xiàn)有方法的50%左右；

(2)建議方法的誤差主要包括降維近似誤差和Kriging 近似誤差，兩種誤差有時疊加、有時抵消；總體而言，建議方法與現(xiàn)有方法的精度大體相當(dāng)；

(3)建議的米N-2和米D-2中，米N-2具有相對更好的精度，可優(yōu)先采用。