基于雙變量降維模型和Kriging近似的統(tǒng)計(jì)矩點(diǎn)估計(jì)法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1 Kriging 近似模型

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

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

1)計(jì)算節(jié)點(diǎn)的選取策略

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

圖1 選點(diǎn)方案示意圖Fig.1 Schematic diagram of point selection

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 算法步驟

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

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

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

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

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ù)，參考點(diǎn)為uc=T?1(μX)=

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

表1 統(tǒng)計(jì)矩計(jì)算結(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是獨(dú)立的正態(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四種方法計(jì)算式(20)所示功能函數(shù)的前四階矩，標(biāo)準(zhǔn)解為多維Gauss-Hermite數(shù)值積分迭代收斂的結(jié)果，結(jié)果如表2所示。從表2可以發(fā)現(xiàn)：1)高階矩計(jì)算結(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 近似后兩者的計(jì)算效率提升1倍多，效果顯著。

算例3.工程算例

式中：MUF、UTS、δ、N、R、R0分別是材料利用率、極限拉伸強(qiáng)度、密度、轉(zhuǎn)子轉(zhuǎn)速、外半徑和內(nèi)半徑；上述各變量均為獨(dú)立正態(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四種方法計(jì)算式(21)所示功能函數(shù)的前四階矩，標(biāo)準(zhǔn)解為3×107個(gè)樣本點(diǎn)的蒙特卡羅法的計(jì)算結(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，計(jì)算效率顯著提高。

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

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

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

式中：N1、N2分別為參考點(diǎn)坐標(biāo)均不為0時(shí)N-2與米N-2的函數(shù)調(diào)用次數(shù)；N3、N4分別為參考點(diǎn)坐標(biāo)均為0時(shí)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

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

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

4 結(jié)論

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

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

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

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