基于擴展型共軛無跡變換的隨機不確定性傳播分析方法

王 濤，李正良,2，范文亮,2

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

由于不可控制性，無論是結構的輸入(如地震、風、波浪荷載等)還是結構參數(shù)(如尺寸參數(shù)和材料參數(shù)等)均呈現(xiàn)隨機特性[1]。相對于將上述參數(shù)取固定值的確定性描述方式，引入隨機變量、隨機過程和隨機場的隨機系統(tǒng)描述更能合理地表征結構的真實狀態(tài)。由于輸入和參數(shù)的隨機性，結構響應亦需從隨機的觀點加以考察。量化輸入和參數(shù)的不確定性對結構響應的影響是典型不確定性傳播分析。理論上，隨機過程和隨機場均可由隨機變量描述，因此僅包含隨機變量系統(tǒng)的不確定性傳播分析是隨機系統(tǒng)不確定性分析的基礎。

眾所周知，隨機變量存在兩類描述方式：概率密度描述和矩描述。相對于概率密度函數(shù)信息獲取的困難，隨機變量的矩描述更易實現(xiàn)，因而獲得了廣泛的關注。理論上，隨機變量的精確描述需獲得其無窮階矩，但其不現(xiàn)實亦不必要。目前大多采用隨機變量的低階矩描述，且主要集中在前幾階矩。

隨機系統(tǒng)響應函數(shù)的統(tǒng)計矩估計實際上是求解響應函數(shù)的概率積分。由于響應函數(shù)通常比較復雜甚至無顯式表達式，獲取其解析解異常困難，相對而言，獲得其數(shù)值解是實用且有效的途徑，其主要包括Monte Carlo模擬法、降維積分法和選點策略方法等。

Monte Carlo模擬法原理簡單、易于實現(xiàn)，且可以給出比較好的結果，但由于計算費用昂貴，往往僅用于其他方法的校核。

降維積分法通過高維分解模型將原函數(shù)分解為多個低維函數(shù)的加和或者乘積形式，主要包括基于單變量加和或乘積降維近似的統(tǒng)計矩估計方法[2? 3]與基于廣義加和降維近似的統(tǒng)計矩估計方法[4?5]。基于單變量降維近似的統(tǒng)計矩估計方法效率高，但加和模型對于交叉項貢獻較為顯著的響應函數(shù)的精度較差，而乘積模型則不適用于以加和形式為主的響應函數(shù)?；趶V義降維近似的統(tǒng)計矩估計方法對于各類響應函數(shù)均具有較好的適用性，但對于高維隨機系統(tǒng)，即使降維的維數(shù)d取2或3[4?7]，由于雙變量分量函數(shù)或三變量分量函數(shù)的數(shù)量較多，積分計算所涉及的結構分析次數(shù)仍較多，計算效率有待于進一步提高。

選點策略方法著重于選取合適的點集進行統(tǒng)計矩的高效計算，本質上是數(shù)值積分方法在統(tǒng)計矩估計中的應用與拓展。該類方法主要包括基于高斯求積公式的張量積分方法[8?10]、稀疏網(wǎng)格積分方法[11? 14]以及無跡變換積分方法[15?16]。張量積分方法具有很高的計算精度，但維數(shù)災難問題不可避免[10]。稀疏網(wǎng)格積分方法利用一維配置點的特殊張量積操作進行線性組合來構建多維求積公式，較張量積分效率更高，可以一定程度上緩解維數(shù)災難；但稀疏網(wǎng)格積分對交叉項影響顯著的情形精度有所減弱[13]；此外，此方法可能會出現(xiàn)負權系數(shù)[13]，導致較大的誤差[14]。無跡變換積分方法是緩解維數(shù)災難的另一類途徑[15?16]，其思路是通過單維或多維單項式精確積分的約束方程獲取求積節(jié)點和權系數(shù)，具有較高的效率，但亦存在負權系數(shù)的問題[17]。

值得指出地是，研究者提出的共軛無跡變換方法(CUT)可以避免負權系數(shù)的產生，具有更好的穩(wěn)健性。然而，目前該方法僅應用于非線性濾波領域的期望計算[17?18]，且僅適用于包含少量高斯分布與均勻分布變量的低維問題，對于可能涉及任意隨機變量類型和高維問題的隨機系統(tǒng)分析并不適用。為此，本文擬將正態(tài)-非正態(tài)變換與共軛無跡變換方法相結合，發(fā)展可適用于任意隨機變量類型的第Ⅰ類擴展型共軛無跡變換方法；在此基礎上，引入高維分解模型，進一步將擴展型共軛無跡變換方法拓展至高維隨機系統(tǒng)分析。

1 問題描述

對于一個典型的不確定性傳播分析問題，隨機系統(tǒng)的響應函數(shù)Z可表示為：

式中：Θ={Θ1,Θ2,···,ΘN}為N維隨機向量；g(?)為表示響應量Z與Θ 映射關系的函數(shù)。

通常，響應函數(shù)的統(tǒng)計矩和概率密度函數(shù)是不確定性傳播分析所關注的問題。由概率論可知，隨機響應函數(shù)Z的k階統(tǒng)計矩MZ,k的表達式如下：

式中：pΘ(θ)為Θ 的聯(lián)合概率密度函數(shù)；MZ,1為其均值；MZ,k為其k階中心矩。

從式(2)看出，響應函數(shù)統(tǒng)計矩的估計實質是一個N維積分求解的過程。雖然研究者發(fā)展了多種統(tǒng)計矩估計方法，但現(xiàn)有方法存在要么計算成本昂貴，要么計算精度不足等問題[10?16]。故而，兼顧精度與效率的統(tǒng)計矩估計方法有待進一步發(fā)展。

2 共軛無跡變換方法

2.1 共軛無跡變換方法

共軛無跡變換方法本質上屬于數(shù)值積分方法，與現(xiàn)有各類數(shù)值積分方法類似，確定求積節(jié)點及其權系數(shù)亦是其核心環(huán)節(jié)，但確定規(guī)則存在顯著差異。共軛無跡變換方法中，確定求積節(jié)點及其權系數(shù)主要基于如下兩個準則：1)定性準則：求積節(jié)點需位于主軸和共軛軸上，其中主軸為坐標軸，共軛軸為象限內的軸；2)定量準則：一定階次的單維和多維單項式積分可由節(jié)點及權系數(shù)精確計算。顯然階次越高，維數(shù)越高，確定節(jié)點和權系數(shù)所需的方程數(shù)越多，求解越困難；若積分的權函數(shù)為非對稱函數(shù)，則既存在奇數(shù)階方程，亦存在偶數(shù)階方程，而當權函數(shù)對稱性時則僅需偶數(shù)階方程；此外，對于非對稱的權函數(shù)，共軛軸亦不具有對稱性，確定過程更為復雜，因此目前的共軛無跡變換方法只發(fā)展了適用于對稱權函數(shù)的4階、6階與8階方法，為簡便，可分別記為CUT-4、CUT-6及CUT-8。

下面將簡單介紹共軛無跡變換方法確定主軸、共軛軸以及求積節(jié)點和權系數(shù)的原理。

2.1.1求積節(jié)點的定性確定準則

對于高斯變量或均勻變量，利用其概率分布的對稱性，構造其N維空間的主軸、共軛軸和縮放共軛軸以及其對應的基準點集如下。

1)主軸：稱以原點為中心的正交坐標軸為主軸，記為σ。根據(jù)定義可知，N維空間存在N個主軸，則主軸上對應的基準點集σi可表示為：

稱該點集與原點連線所構成的軸為Mth共軛軸，記為CM。

稱該點集與原點連線所構成的軸為Nth縮放共軛軸，記為SN(h)。

顯然，主軸、共軛軸以及縮放共軛軸上各自的基準點集均分布在不同的圓上，共軛無跡變換方法中求積節(jié)點確定的定性準則為：求積節(jié)點由上述各基準點集分別縮放后確定的點集組合而成。易知，縮放后的基準點集仍然位于主軸、共軛軸以及縮放共軛軸上。為了便于理解，以如圖1所示的二維空間為例，圖中的σ、C2及S2分別為構造出的主軸、共軛軸和縮放共軛軸，各軸上對應的點即為CUT 的求積節(jié)點。

圖1 二維空間主軸σ，共軛軸C2，縮放共軛軸S2 及CUT點Fig.1 Principal axisσ,conjugateaxis C2,scaled conjugate axes S2 and CUT points of 2D space

2.1.2求積節(jié)點和權系數(shù)的定量確定準則

根據(jù)求積節(jié)點確定的定性準則可知：各基準點集縮放后確定的求積節(jié)點集僅取決于距離尺度變量ri，且根據(jù)變量概率分布的對稱性可知上述節(jié)點集中各節(jié)點具有相同的權系數(shù)，記為wi。因此，合理確定ri與wi是共軛無跡變換方法的核心環(huán)節(jié)。目前，常用的求積節(jié)點與權系數(shù)確定的定量規(guī)則是通過Isserlis定理[19]導出ri和wi滿足的矩約束方程。

以二維系統(tǒng)的CUT-4為例，節(jié)點和權系數(shù)確定的定量規(guī)則為：不高于4階的單項式期望均可精確計算。于是可建立如下約束方程：

式中：U1、U2表示隨機變量；u1,i、u2,i表示求積節(jié)點的坐標；κ和ν為非負的整數(shù)且1≤κ+ν≤4；n為除零點外求積節(jié)點的數(shù)量。略去式(7)中的奇數(shù)階方程且將u1,i、u2,i用相應的ri替換，則式(7)可改寫為：

類似地，確定N維系統(tǒng)CUT-4的節(jié)點和權系數(shù)的約束方程為：

據(jù)此，CUT-4的積分節(jié)點u i={u1,i,u2,i,···,uN,i}和權系數(shù)αi可由表1給出。

表1 CUT-4方法的積分節(jié)點和權系數(shù)Table 1 Integral nodes and weightsof CUT-4 method

類似地，亦可確定CUT-6和CUT-8的積分節(jié)點與權系數(shù)，推導過程及結果可參考文獻[17]。

值得指出地是，隨著維度增加，非線性方程不易求解，且可能出現(xiàn)負權系數(shù)?，F(xiàn)有的CUT方法中，對于CUT-6當N≤9或CUT-8當N≤6時，可確保權系數(shù)為正值，因此當前僅發(fā)展了N≤9的CUT-6及N≤6的CUT-8[17]。

2.2 存在的問題

不難發(fā)現(xiàn)，現(xiàn)有CUT盡管較好地兼顧了精度和效率，但其適用范圍非常有限。若欲將其拓展至一般的隨機系統(tǒng)分析，以下問題亟待解決：

1)CUT利用了積分權函數(shù)的對稱性，目前僅在包含高斯分布或均勻分布隨機系統(tǒng)的期望計算中獲得了應用[17]，并不適用于涉及任意分布隨機變量的隨機系統(tǒng)分析。

2)現(xiàn)有CUT-6與CUT-8僅適用于低維隨機問題，對于高維問題不再適用。

3 擴展型共軛無跡變換方法

針對CUT的上述兩個問題，本文提出了兩類擴展型共軛無跡變換方法。

3.1 基于正態(tài)-非正態(tài)變換的第Ⅰ類擴展型共軛無跡變換方法

為將CUT 拓展至可適用于包含任何隨機變量的系統(tǒng)，可引入正態(tài)-非正態(tài)變換，即對于任意相互獨立的隨機向量Θ={Θ1,Θ2,···,ΘN}與相互獨立的標準正態(tài)向量U={U1,U2,···,UN}存在如下關系：

相應地，mZ,k可表示為：

式中：φ(?)表示標準正態(tài)分布的概率密度函數(shù)；u i=(u1,i,u2,i, ···,uN,i)，u i和αi分別為CUT(如CUT-4、CUT-6或CUT-8)的求積節(jié)點和權系數(shù)。

通過式(12)的變換，可將CUT拓展至包含任意隨機變量的隨機系統(tǒng)響應的統(tǒng)計矩估計，適用范圍顯著拓寬。本文稱此方法為第Ⅰ類擴展型共軛無跡變換方法，為簡便記為ECUTI。顯然，當所有隨機變量均為高斯分布時，ECUTI 退化為傳統(tǒng)的CUT。

3.2 基于高維分解模型的第Ⅱ類擴展型共軛無跡變換方法

盡管ECUTI 可適用于任意隨機變量，但其中ECUTI-6、ECUTI-8與CUT-6、CUT-8類似，僅分別適用于N≤9及N≤6的情形。然而，工程實際中，高維隨機系統(tǒng)并不鮮見，為此，可將高維分解模型與ECUTI 相結合，發(fā)展具有更廣泛適用性的方法。

根據(jù)高維分解模型方法，Z=h(U)可近似為多個最高維數(shù)為d(d

式中：uq,i和ωi分別為Gauss-Hermite求積公式的求積節(jié)點和權系數(shù)；m為一維積分的節(jié)點數(shù)量，Np為二維CUT的節(jié)點數(shù)量，且m與Np的取值應保證一維積分與二維積分具有相同或相近的代數(shù)精度；u i={us,i,ut,i}與αi分別為二維CUT的節(jié)點與權系數(shù)，可由式(17)確定：

文中稱上述方法為第Ⅱ類擴展型共軛無跡變換方法，為簡便記為ECUTII。

表2 二維系統(tǒng)CUT方法對應的ri 和wiTable 2 ri and wi in CUT method for two-dimensional system

3.3 效率分析

若記CUT、ECUTI及ECUTⅡ的函數(shù)調用次數(shù)分別為Npo、NpI及NpII，顯然，Npo=NpI，且有：

式中：Npo4、Npo6及Npo8分別為CUT-4、CUT-6及CUT-8的函數(shù)調用次數(shù)；NpⅡ4、NpⅡ6及NpⅡ8分別為ECUTII-4、ECUTII-6及ECUTII-8的函數(shù)調用次數(shù)。

三類方法的函數(shù)調用次數(shù)Np隨著N的變化規(guī)律如圖2所示。不難發(fā)現(xiàn)，當N小于等于維度界限Nb時，ECUTⅠ(CUT)的計算效率高于ECUTⅡ；然而當N>Nb時，ECUTⅡ的計算效率將高于ECUTⅠ(CUT)。具體而言，對于CUT-4、CUT-6及CUT-8，Nb分別取8、6及3。

圖2 各類CUT 方法效率對比Fig.2 Efficiency comparison of various CUT methods

4 算例分析

為了研究提出的方法的精度、效率以及適用性，本節(jié)將對三個算例進行分析。算例1考慮了一個隨機變量均為高斯變量的問題，用于驗證提出的方法的計算效率以及計算精度；算例2和算例3分別考慮了具有非高斯變量的隨機問題與高維隨機問題，用于驗證提出的方法對于各類隨機問題的普適性。為了驗證方法的計算效率與精度，分別將提出的方法分別與Monte Carlo模擬方法(MCS)、7點張量數(shù)值積分方法(FTM)[17]、雙變量降維近似方法(BDRM)[5]及2階精度的稀疏網(wǎng)格方法(SGI-2)[12]進行了對比。本文將MCS計算結果視為精確解，其他各方法的相對誤差為：

式中：Value表示對應方法的計算結果；MCV表示Monte Carlo方法計算的結果。

4.1 算例1：隨機變量為高斯分布類型的非線性問題

考慮一個如圖3所示的屋架結構，屋架的上弦桿和其他壓桿采用鋼筋混凝桿，下弦桿和其他拉桿采用鋼桿。設屋架承受均布載荷q作用，將均布載荷q轉化成節(jié)點載荷P=0.25ql后,C點處延垂直面方向的位移對應的功能函數(shù)可表示為[22]：

圖3 屋架結構Fig.3 Roof structure

式中，共包含6個隨機變量，其物理意義與統(tǒng)計特性如表3所示。采用各類方法計算出功能函數(shù)的前四階矩如表4所示。

表3 算例1中隨機變量的統(tǒng)計特征Table 3 Statistical characteristics of the random variablesfor Example1

通過表4可以看，就精度而言，采用各類方法可以較為精確地估計前兩階矩，然而對于高階矩的估計，采用SGI 和ECUTⅡ-4的相對誤差較大，F(xiàn)TM的結果與精確解相同，而BDRM、ECUTⅡ-6及ECUTⅡ-8由于高維分解模型的截斷，存在一定的誤差，但其相對誤差低于7%，原始的CUT均具有很高的精度，和精確解的結果相同；從效率上看，F(xiàn)TM所需函數(shù)調用次數(shù)較多；而CUT-4和CUT-6具有明顯的優(yōu)勢，且優(yōu)于ECUTⅡ-4、ECUTⅡ-6及BDRM；ECUTⅡ-8的函數(shù)調用次數(shù)高于CUT-8；雖然SGI方法僅需要較少的函數(shù)調用次數(shù)，然其精度不足。整體而言，對于服從正態(tài)分布的低維隨機問題，本文推薦采用CUT；當容許較小誤差且注重計算效率時，相比于CUT-8，ECUTⅡ-8亦被推薦。

表4 算例1中極限狀態(tài)函數(shù)前四階統(tǒng)計矩的結果Table 4 Results of the first four order moments of the limit state function for Example 1

4.2 算例2：具有不同分布類型的隨機問題

考慮一個倒三角荷載作下的6層框架結構，如圖4所示?？蚣艿目缍萀=7.5 m，層高H=3 m，梁與柱的截面尺寸分別為300 mm×400 mm 和500 mm×500 mm，水平荷載F可表示為：

圖4 6層框架結構Fig.4 Six story frame structure

該問題中的隨機變量服從非正態(tài)分布類型，傳統(tǒng)的CUT 不再適用。首先，基于式(12)，將響應函數(shù)Z中各隨機變量變換為獨立標準正態(tài)向量，可得：

進而，分別采用MCS、FTM、SGI、BDRM以及提出的ECUTI和ECUTⅡ計算Z的前四階矩，其結果如表6所示。

表5 算例2中隨機變量的統(tǒng)計特征Table 5 Statistical characteristics of the random variables for Example 2

表6 算例2中極限狀態(tài)函數(shù)前四階統(tǒng)計矩的結果Table 6 Results of the first four order moments of the limit state function for Example 2

通過表6可以看出，對于非正態(tài)隨機變量類型，響應函數(shù)Z的前兩階矩均可以被精確估計。對于三階、四階矩的估計，SGI-2、ECUTI-4和ECUTⅡ-4估計的相對誤差較大；FTM 和ECUTI-8的估計結果和精確解相同，但FTM 所需的有限元分析次數(shù)較多，而ECUTI-8僅需要355次有限元分析；BDRM、ECUTI-6、ECUTⅡ-6及ECUTⅡ-8估計的三階、四階矩的相對誤差較低，最大相對誤差為8.25%，且提出的ECUTI和ECUTⅡ的有限元分析次數(shù)均低于BDRM，具有更好的效率。值得指出地是，表6中ECUTⅡ-6的計算精度高于ECUTⅡ-8，造成這一現(xiàn)象的原因可能是高維分解模型引入的誤差與CUT-6本身的誤差產生了抵消效應。整體而言，對于低維問題，本文推薦使用ECUTI-6和ECUTI-8；同時，當容許少許誤差且注重計算效率時，ECUTⅡ-6和ECUTⅡ-8亦可被推薦采用。

4.3 算例3：高維隨機問題

考慮一個高維隨機問題，其極限狀態(tài)函數(shù)的表達式為：

式中，Θi(i=1,2,···,d)是獨立同分布的隨機變量，均服從對數(shù)正態(tài)分布，其均值為1，標準差為0.2，維度d=20。

對于此類高維問題，CUT-6、CUT-8、ECUTI-6以及ECUTI-8均不再適用，ECUTI-4及ECUTⅡ可解決此類問題。采用各類方法估計高維隨機問題的前四階矩，其計算結果如表7所示。

表7 算例3中極限狀態(tài)函數(shù)前四階統(tǒng)計矩的結果Table 7 Resultsof the first four order momentsof the limit state function for Example 3

通過表7可以看出，對于前兩階矩，各類方法均具有較高的精度；對于三階、四階矩的估計，F(xiàn)TM 方法的結果接近精確解但函數(shù)調用次數(shù)高達1016量級；雖然采用精度水平為2的SGI-2方法最為高效，但方法精度不足，若需達到較高精度，SGI方法(SGI-3)需要的函數(shù)調用次數(shù)為120 321，其次數(shù)遠遠高于BDRM 和本文提出的方法；ECUTI-4 和ECUTⅡ-4估計的三階矩的相對誤差較大且ECUTI-4函數(shù)調用次數(shù)達百萬量級；本文提出的ECUTⅡ-6 與ECUTⅡ-8能夠兼顧精度與效率。整體而言，對于高維問題，本文推薦采用ECUTⅡ-6與ECUTⅡ-8，其中，ECUTⅡ-6函數(shù)調用次數(shù)較少而ECUTⅡ-8的精度更高。

5 結論

統(tǒng)計矩是描述隨機系統(tǒng)不確定性傳播分析的有效方式之一，實用時以低階矩描述為主。本文將共軛無跡變換方法引入到隨機系統(tǒng)不確定性傳播分析中，并針對現(xiàn)有方法的不足，結合正態(tài)-非正態(tài)變換以及高維分解模型，提出了兩類適用于不確定性傳播分析的擴展型共軛無跡變換方法。文中分析結果表明：

(1)相較于傳統(tǒng)的統(tǒng)計矩估計方法，本文建議的兩類方法均可兼顧精度與效率，對于隨機系統(tǒng)可以進行高效的響應統(tǒng)計矩估計。

(2)建議的兩類方法均可適用于任意隨機變量類型的隨機系統(tǒng)，擴展了傳統(tǒng)的共軛無跡變換方法的適用范圍。

(3)對于低維問題，建議直接采用ECUTI；對于高維問題，建議采用ECUTⅡ，其誤差主要包括CUT 本身的誤差和高維分解模型的截斷誤差，兩種誤差有時疊加、有時抵消。

需指出地是，盡管建議方法具有較高的計算效率，但是對于特別高維的問題，仍可能涉及成千上萬次的結構分析次數(shù)。若能夠在建議方法的基礎上進一步改善計算效率將是非常有益的，亦是值得進一步深入研究的方向。