隨機參數(shù)系統(tǒng)不確定性量化的泛函觀點與整體靈敏度分析1)

萬志強 陳建兵,2) Michael Beer

?(同濟大學土木工程防災國家重點實驗室,上海 200092)

+(同濟大學土木工程學院,上海 200092)

??(萊布尼茨漢諾威大學風險與可靠性研究所,德國漢諾威30167)

引言

現(xiàn)代化土木工程建設正逐步走向整體化、精細化及全概率化的全新分析與設計范式.與之相應的是,數(shù)值仿真計算模型已然成為工程結構系統(tǒng)設計與分析中不可或缺的一部分.然而,設計與分析過程中不可避免地會遇到諸多不確定性因素的影響,如仿真計算模型的邊界約束可能具有不確定性、初始條件可能難以精確給定、系統(tǒng)的某些力學性能參數(shù)可能具有不可忽略的隨機性等[1-2].一般而言,這些不確定性因素均可表征為某種參數(shù)化的形式,如混凝土材料的抗壓強度與彈性模量[3]、巖土材料的有效黏聚力與有效內摩擦角[4]、隨機風場模型的地面粗糙度與剪切速度[5]等.

由于系統(tǒng)輸入?yún)?shù)存在不確定性,必然導致系統(tǒng)輸出亦存在一定程度的不確定性,因而仿真預測的結果不可避免地與實際情況存在一定偏離.研究表明,不同系統(tǒng)輸入?yún)?shù)對系統(tǒng)輸出不確定性的“貢獻”一般是不同的.因此,量化系統(tǒng)輸入?yún)?shù)對系統(tǒng)輸出不確定性的影響(或稱量化系統(tǒng)輸出關于系統(tǒng)輸入?yún)?shù)的靈敏度),即靈敏度分析,對工程設計、優(yōu)化與決策具有重要意義[6].

靈敏度分析的目標是量化系統(tǒng)輸入?yún)?shù)對系統(tǒng)輸出不確定性的影響,即選擇并計算靈敏度指標.例如,基于差商的基本效應指標[7](亦稱為Morris 指標)及基于導數(shù)的靈敏度指標[8],主要適用于確定性系統(tǒng),在本質上屬于局部靈敏度的范疇.另一方面,針對隨機參數(shù)系統(tǒng),人們發(fā)展了基于方差分解的Sobol’指標[9]以及基于條件概率密度函數(shù)的矩獨立指標[10]等整體靈敏度指標.這些指標側重于對靈敏度“大小”的計算,因而往往其值恒為非負.這容易使人產生輸入變量的不確定性總是使得輸出變量的不確定性增大的錯覺[11].事實上,系統(tǒng)輸入?yún)?shù)對輸出不確定性的影響是有“方向”的,而這一“方向”對結構分析、特別是考慮不確定性影響下的結構優(yōu)化設計至關重要[12].例如,基于失效概率的靈敏度指標便是一類能同時反映靈敏度“大小”與“方向”的指標[13-14],被廣泛用于基于可靠度的結構優(yōu)化與設計之中[12-15].

在上述研究工作的基礎上,Chen 等[16]提出了一類全新的基于Fr′echet 導數(shù)的整體靈敏度指標,并給出了相應的高效算法.其中,隨機系統(tǒng)中不確定性傳播的泛函觀點雖已初具雛形,但尚未得以深入論述.同時,如何針對具體的工程問題采用基于Fr′echet 導數(shù)的整體靈敏度指標,特別是對參數(shù)進行靈敏度排序等相關問題的研究尚不夠充分.

基于此,本文從隨機系統(tǒng)中不確定性傳播的泛函表達角度,進一步探討了該整體靈敏度指標的數(shù)學定義、物理意義和基本性質.進而,在ε?等效分布的定義下,討論了該整體靈敏度指標的參數(shù)化表達形式,給出了其對應的重要性測度和方向靈敏度的定義及其物理意義.建立了該整體靈敏度指標與矩獨立指標之間的內在聯(lián)系.結合概率密度演化?概率測度變換方法[17],提出了該整體靈敏度指標及其對應的重要性測度與方向靈敏度的數(shù)值計算方法和流程.最后,通過4 則算例分析詳細討論了該靈敏度指標的工程實用性.

1 隨機系統(tǒng)中不確定性量化與傳播的泛函觀點

1.1 隨機變量變換的泛函空間觀點

在概率論框架下,隨機系統(tǒng)中輸入變量的不確定性可以隨機變量(向量)的數(shù)學形式進行表征[18-19].記概率空間為(?,F,),其中? 為樣本空間、F 為σ-代數(shù)(? 的事件域)、為概率測度,? ∈? 代表樣本空間中的一個點(即基本事件).不失一般性,考慮一隨機系統(tǒng)

其中,g(·):Rn→Rm為一連續(xù)映射,其分量為(g1,g2,···,gm)T;Θ(?)=(Θ1(?),Θ2(?),···,Θn(?))T為n維輸入隨機向量,其聯(lián)合概率密度函數(shù)pΘ(θ)已知,其中θ=(θ1,θ2,···,θn)T為隨機向量Θ的可實現(xiàn)值;X(?)=(X1(?),X2(?),···,Xm(?))T為m維系統(tǒng)輸出向量,記其聯(lián)合概率密度函數(shù)為pX(x),其中x=(x1,x2,···,xm)T為隨機向量X的可實現(xiàn)值.

由概率論可知,若已知pΘ(θ),則根據(jù)上式可確定輸出隨機向量的概率密度函數(shù)pX(x).事實上,根據(jù)概率守恒原理[20-21],可得

其中?Θ為輸入隨機向量值域空間的任意子域,?X為由式(1)中的映射決定的在輸出隨機向量值域空間中與?Θ對應的子域.由此可得

為簡單計,先考察n=m=1 的情況.此時,式(1)成為標量形式X(?)=g(Θ(?)),而式(3)則相應地退化為∫?X pX(x)dx=∫?ΘpΘ(θ)dθ.注意此時積分區(qū)域為?X=g(?Θ).進一步,考慮下列情況.(1)函數(shù)g(·)為單調函數(shù).此時,有

其中λ(·)表示子域的Lebesgue 測度,對一維情況即為其子域的長度;為Jacobi 行列式,對一維情況即為其導數(shù)之絕對值,g?1(·)表示g(·)的反函數(shù).

(2)函數(shù)g(·)為非單調函數(shù).此時,總可以將其定義域剖分為若干個不相交子域,并使之在每個子域中為單調函數(shù)[21].記其相應的反函數(shù)為(·)(j=1,2,···,N),N為子域劃分段數(shù),此時有

對其進一步化簡可得

另一方面,也可以通過Dirac 函數(shù)δD[·]的篩選性質體現(xiàn)的映射關系獲得輸出變量的概率密度函數(shù).事實上,對于上述情形,有

其中pX|Θ(x|Θ=θ)為條件概率密度函數(shù).由于存在映射關系X(?)=g(Θ(?)),因而有pX|Θ(x|Θ=θ)=δD[x?g(θ)],由此可得式(8)中的第二個等式.

顯然,可認為式(8)是輸出概率密度函數(shù)pX(x)關于輸入概率密度函數(shù)pΘ(θ)的隱式表達,而式(4)、式(5)和式(7)則是對式(8)進行積分獲得的顯式解析表達.上述各式表明,輸出概率密度函數(shù)pX(x)由輸入概率密度函數(shù)pΘ(θ)所決定,因而可以進一步寫成如下算子形式的表達

其中下標g用來表示算子ψg決定于函數(shù)g(·).在實際工程中,函數(shù)g(·)描述了問題的物理機制,因此,式(9)表明:不確定性在系統(tǒng)中的傳播,即從輸入概率信息到輸出概率信息的變換,是由系統(tǒng)物理機制所決定的.從函數(shù)空間來看,這一變換的實現(xiàn)是通過由物理機制決定的算子作用來實現(xiàn)的.

對于m≥2 的情況,上述基本思想是完全類似的.限于篇幅,茲不贅述.

1.2 隨機動力系統(tǒng)中不確定性傳播的泛函空間觀點

這一泛函空間觀點,對于隨機動力系統(tǒng)中的不確定性傳播也是適用的.為了明確起見,考慮如下隨機動力系統(tǒng)[21]

其中X=(X1,X2,···,Xm)T為狀態(tài)向量,˙X=?X/?t;G=(G1,G2,···,Gm)T為確定性函數(shù)向量;Θ=(Θ1,Θ2,···,Θn)T為刻畫隨機源信息的隨機向量,其概率密度函數(shù)pΘ(θ)已知,這里θ=(θ1,θ2,···,θn)T為其對應的可實現(xiàn)值.式(10)的初始條件為X(0)=x0,其中x0為確定性初始值.

根據(jù)概率守恒原理[20],記(X(t),Θ)的聯(lián)合概率密度函數(shù)為pXΘ(x,θ,t),則其滿足如下廣義概率密度演化方程[21]

其初始條件為

求解上述方程即可得X(t)的概率密度函數(shù)

由此可見,若已知初始隨機源概率信息(輸入概率信息)pΘ(θ),即可通過求解式(11)、式(12)和式(13)獲得輸出概率信息pX(x,t).因此,上述3 個方程決定了從輸入概率密度函數(shù)到輸出概率密度函數(shù)的變換.從函數(shù)空間來看,即存在如下算子方程

這里用下標G表示算子ψG由方程(10)所決定,而方程(10)是描述系統(tǒng)演化內在物理機制的方程.因此,從這一角度來看,概率密度演化理論[20-21]正是算子方程(14)的具體實現(xiàn)過程,而不確定性在系統(tǒng)中從輸入概率信息到輸出概率信息的傳播,是由系統(tǒng)物理機制(通過確立傳播算子)來決定的.從這一意義上說,ψG可稱為不確定性傳播算子.

事實上,式(9)和式(14)中的算子ψg和ψG,正是Frobenius-Perron 算子[21-22].值得注意的是,上述算子是線性算子,僅依賴于物理方程(例如以上的g(·)和G(·)),而與輸入概率信息(例如以上pΘ(θ))的形式是無關的.換言之,物理系統(tǒng)本身不因隨機源概率信息的變化而變化.這一隨機系統(tǒng)客觀不變性原理,正是上述不確定性傳播算子具有不變性的內在原因[23].

2 基于Fr′echet 導數(shù)的整體靈敏度

2.1 基于Fr′echet 導數(shù)的整體靈敏度指標定義

基于上述泛函空間理解,Chen 等[16]提出了基于Fr′echet 導數(shù)的整體靈敏度指標.為清晰起見,下面以隨機系統(tǒng)(10)與相應的不確定性傳播算子方程(14)為基礎簡要闡述其基本思想.

考慮如下問題:若輸入概率密度函數(shù)發(fā)生了變化,例如從pΘ(θ)變?yōu)?pΘ(θ),則系統(tǒng)響應(輸出)的概率密度函數(shù)也將相應地發(fā)生改變,且由算子方程(14)可知

為方便起見,記輸入概率密度函數(shù)的變化為δpΘ=(θ)?pΘ(θ),輸出概率密度函數(shù)的變化為δpX=(x,t)?pX(x,t).由式(15)兩端分別減去式(14)可得由此可見,若輸入概率信息具有變化,則輸出概率信息的變化亦由不確定性傳播之線性算子所決定.

事實上,上述概率信息“變化”的傳播算子,就是不確定性傳播算子的Fr′echet 導數(shù),其嚴格的數(shù)學定義為[24]:

令V和W為Banach 空間,定義映射ψ:V→W.若存在一有界線性算子Fψ∈L(V,W),滿足

則稱Fψ為在u0∈U?V處的算子ψ 的Fr′echet 導數(shù),記為Fψ=ψ′[u0].式中,∥·∥V和∥·∥W分別對應Banach空間V和W上的范數(shù).

在本文中,Banach 空間V上的范數(shù)采用基于總變差距離的范數(shù)

其中,L1(?Θ)為定義在?Θ上的?1范數(shù),即

而空間W上的范數(shù)同樣可定義為

由于不確定性傳播算子是線性算子,因此概率信息“變化”的傳播算子就是不確定性傳播算子本身.

在實際計算中,可以利用Fr′echet 導數(shù)與G?ateaux導數(shù)的等價性.G?ateaux 導數(shù)的定義為:令V和W為Banach 空間,定義映射ψ :V→W.若存在一有界線性算子Gψ∈L(V,W),滿足

則稱Gψ為在u0∈U?V處的算子ψ 的G?ateaux 導數(shù),記為Gψ=ψ′[u0].

顯然,Fr′echet 導數(shù)即為G?ateaux 導數(shù).反之,若G?ateaux 導數(shù)在點u0處連續(xù)或式(21)對h一致成立且∥h∥V=1,則點u0處的G?ateaux 導數(shù)亦為Fr′echet導數(shù)[24].

無論是確定性系統(tǒng)還是隨機系統(tǒng),系統(tǒng)靈敏度的定義都可以闡述為[25]“系統(tǒng)靈敏度是輸入量的擾動導致輸出量的擾動這一效應的度量”.基于這一觀點,在經(jīng)典的確定性系統(tǒng)分析中往往采用函數(shù)的偏導數(shù)作為系統(tǒng)靈敏度的度量.顯然,在隨機系統(tǒng)中采用算子的Fr′echet 導數(shù)作為系統(tǒng)的靈敏度度量,乃是一種自然的推廣.但是,以函數(shù)的偏導數(shù)作為確定性系統(tǒng)的靈敏度,僅能反映系統(tǒng)在給定輸入值一點處的變化情況,因而是一種局部的度量;而對隨機系統(tǒng)而言,Fr′echet 導數(shù)給出了在系統(tǒng)輸入全部可能值(通過概率密度函數(shù)刻畫)處的變化導致在系統(tǒng)輸出全部可能值上的變化情況,因而其本質是一種整體靈敏度指標(global sensitivity index,GSI)[16].

基于Fr′echet 導數(shù)的整體靈敏度指標有4 個重要性質:

(1)無論物理系統(tǒng)本身是線性還是非線性的,因為Fr′echet 導數(shù)是線性算子,因此基于Fr′echet 導數(shù)的整體靈敏度指標與pΘ的具體形式是無關的(類似于線性函數(shù)x=kθ 中的局部靈敏度dx/dθ=k是與θ 取值無關的).這意味著,這一整體靈敏度指標不僅適用于輸入概率密度信息具有小擾動的情況,也適用于輸入概率密度信息具有大擾動、甚至若干輸入變量從隨機變量變?yōu)榇_定性變量的情況(例如概率密度函數(shù)坍塌為Dirac 函數(shù)).

(2)基于性質(1),若考慮輸入概率模型具有某種程度的認知不確定性[1-2,17-18],例如,輸入概率密度函數(shù)可能為,則這種認知不確定性傳播的量化可以直接由該整體靈敏度指標給出,即Fψ?(.

(3)對于隨機動力系統(tǒng),系統(tǒng)響應是時變的,因此基于Fr′echet 導數(shù)的整體靈敏度指標為一時變量.在此情況下,可根據(jù)問題的需要考慮系統(tǒng)響應的某代表性值,如系統(tǒng)響應時程的極大值,即Zext=max{X(t),t∈[0,T]},簡記為Z.

(4)注意到Fr′echet 導數(shù)是線性有界算子,在式(17)、式(18)和式(20)定義下的范數(shù)為1(證明見第2.3 節(jié)).由于Fr′echet 導數(shù)就是不確定性傳播算子本身,因此,不確定性傳播算子的范數(shù)為1.這在本質上是概率守恒原理[20]在泛函空間觀點下的數(shù)學表述之一.

2.2 整體靈敏度指標的參數(shù)化表達

在工程中,特別具有實際意義的是參數(shù)化表達的整體靈敏度指標.對絕大多數(shù)實際工程問題而言,系統(tǒng)輸入?yún)?shù)Θ一般為連續(xù)型隨機向量,記其概率密度函數(shù)為pΘ(θ;ζ),其中ζ=(ζ1,ζ2,···,ζs)T為其對應的分布參數(shù)向量.此時,pΘ(θ;ζ)的變分為

在文獻[16]中,考慮各分布參數(shù)之間相互獨立變化,即當僅考慮δζj≠0 而其他為0 時,給出了基于Fr′echet 導數(shù)的整體靈敏度指標的參數(shù)化表達

引入式(14),式(23)可以進一步改寫為

這表明基于Fr′echet 導數(shù)的整體靈敏度指標參數(shù)化表達是不確定性傳播算子作用于單位化的輸入函數(shù)變化上的結果.因此,上述整體靈敏度指標參數(shù)化表達的幾何意義為單位化的輸入函數(shù)變化作為抽象空間中的點(向量)經(jīng)系統(tǒng)作用后的“方向”及“長度”變化.

更具體地,在實際工程問題中,有時通過隨機變量Θ 的統(tǒng)計矩參數(shù)來近似構造其概率密度函數(shù).例如,四階矩正態(tài)變換方法假定隨機變量Θ 可表達為標準正態(tài)隨機變量U的多項式函數(shù)[26]

其中,參數(shù)a0,a1,a2和a3由Θ 的前四階矩即其均值μΘ、標準差σΘ、偏度αΘ和峰度βΘ計算得到

由隨機變量的函數(shù)變換法則[19-21],隨機變量Θ的概率密度函數(shù)可表示為

其中,?(u)為標準正態(tài)分布的概率密度函數(shù),反函數(shù)u=u?1(θ)的顯式表達詳見文獻[27].

可見,此時Θ 的概率密度函數(shù)由其前四階矩確定.不妨考慮隨機向量Θ=(Θ1,Θ2,···,Θn)T相互獨立,而對于非獨立情況可通過Nataf 變換[28]、Copula函數(shù)[29]或隨機函數(shù)模型[3]轉換為獨立情況.此時,該隨機向量的概率密度函數(shù)可表示為

式中,pΘ?為第? 個隨機變量Θ?的概率密度函數(shù),分別對應Θ?的均值、標準差、偏度和峰度.由式(23)和式(28)可得隨機系統(tǒng)(1)基于Fr′echet 導數(shù)的整體靈敏度指標為

矩表達形式中前四階矩的物理意義是十分明確的,其對應的全局靈敏度指標反映了各個輸入變量的前四階矩信息變化對輸出概率密度函數(shù)的影響.

更一般地,若Θ的概率密度函數(shù)是已知的但是非參數(shù)化的,如pΘ(θ)是通過核密度估計[1]或概率密度演化理論[21]得到的數(shù)值解,則可通過ε-等效分布的方法將其參數(shù)化.為此,記隨機變量Θ的概率密度函數(shù)為pΘ(θ).若Θ可以由有限個參數(shù)進行表達,則稱(θ;ζ)為pΘ(θ)的ε-等效分布,滿足

其中,ζ為有限維參數(shù)向量,可由優(yōu)化算法[12]求解式(30)給出;ε 為給定的容許誤差,如0.01.注意,這里(θ;ζ)可以是任意形式的.例如,有確定形式的概率密度函數(shù)[2-3]、由矩表達形式給出的概率密度函數(shù)[26-27]或由級數(shù)展開方法[30]計算得到的概率密度函數(shù)等.在上述基礎上將pΘ(θ;ζ)替換為(θ;ζ),則隨機系統(tǒng)(1)基于Fr′echet 導數(shù)的整體靈敏度指標的參數(shù)化表達可統(tǒng)一為式(23).

2.3 重要性測度

如式(23)所示的整體靈敏度指標描述了系統(tǒng)響應X的值域范圍內每一點x處概率密度的靈敏度.進一步,對整體靈敏度指標取∥·∥W范數(shù)可得其重要性測度

現(xiàn)分別考察式(31)的幾何意義和物理意義.

2.3.1 幾何意義

考慮Fr′echet 導數(shù)的算子范數(shù),其定義為[24]

式中各符號同式(17)中定義.任意取Θ的概率密度函數(shù)為p2和p1,令p=p2?p1,由均值定理有[24]

此處利用了不確定性傳播算子是線性算子的性質.

注意tp2+(1 ?t)p1為概率密度函數(shù).因此,由概率相容性和不確定性傳播算子方程(9)有

因此,重要性測度總是小于1.結合前述式(24)的幾何意義表明,重要性測度的幾何意義是單位化(即長度為1)的輸入函數(shù)變化作為抽象空間中的點(向量)經(jīng)系統(tǒng)作用后的長度,這一長度不可能超過1,且該長度越大,則其重要性越大.

2.3.2 物理意義

其中,ej為s維單位向量,且僅其第j位上元素為1而其它為0;dTV(p,q)為概率密度函數(shù)p和q的總變差距離(total variation distance),即[16]

類比于連續(xù)介質力學中的變形梯度并觀察式(36)不難發(fā)現(xiàn),由式(31)所定義的整體靈敏度對應的重要性測度反映了某種“體積變化率”(如圖1 所示):當?shù)趈個分布參數(shù)ζj發(fā)生擾動,使得輸入概率空間產生“體積變化”后,隨機系統(tǒng)對應的輸出概率空間亦相應地產生了“體積變化”,即.由式(36)定義可知,這兩者的比值的物理意義便是由不確定性傳播導致的“體積變化率”.需要注意的是,由式(36)定義的“體積變化率”恒為正值,因此可視其為“體積變化率”的絕對值或將其稱為“幅值”.如此,由式(31)給出的整體靈敏度指標對應的重要性測度的物理意義是十分明確的.

圖1 基于Fr′echet 導數(shù)的整體靈敏度指標對應的重要性測度的物理意義Fig.1 Physical meanings of the importance measure with respect to the Fr′echet derivative based GSI

2.4 方向靈敏度

另一方面,考慮如(23)所示的整體靈敏度指標在系統(tǒng)響應X值域范圍內的方向靈敏度

式中,sgn(A)為符號函數(shù):當A≥0 時取+1(代表“正”或“+”)、否則取?1(代表“負”或“?”).計算整體靈敏度指標在?X上的積分,易得

事實上,由式(31)可知

其中I{A}為示性函數(shù),當A為真時取1,否則取0.

由式(38)所定義的方向靈敏度的意義在于,系統(tǒng)響應的值域?X可以通過式(42)劃分為正、負兩個子值域(可能是多連通域).具體而言:(1)輸入分布參數(shù)的微小正增長,如?ζj>0 將導致正值域內的概率密度函數(shù)增大、負值域內的概率密度函數(shù)減小;反之,輸入分布函數(shù)的微小負增長將導致正值域內概率密度函數(shù)減小,而負值域內的概率密度函數(shù)增大.因此,與分別表示了由于參數(shù)變化導致的概率凈減和概率凈增值域;(2)若按式(38)定義可得有限個正、負值域,則正或負值域的不連續(xù)子值域個數(shù)(如圖2)可反映物理系統(tǒng)的復雜程度.一般而言,系統(tǒng)非線性程度越復雜,其不交叉子值域個數(shù)越多.在第4 節(jié)將結合具體實例對其進行詳細討論.

2.5 與矩獨立靈敏度指標的聯(lián)系

考慮隨機系統(tǒng)(1),基于矩獨立的全局靈敏度指標定義為[10]

而此時pΘi與pΘ的總變差距離為

上式意味著當?shù)趇個隨機變量Θi坍塌為確定性變量時,輸入概率空間產生了單位體積的變化.此時,基于矩獨立的全局靈敏度指標等價于

在第2.3 節(jié)中,式(36)解釋了基于Fr′echet 導數(shù)的整體靈敏度指標的重要性測度的物理意義:它是隨機系統(tǒng)在點pΘ處沿其分布參數(shù)ζj方向的(絕對)體積變化率.類比于此,式(47)可理解為第i個隨機變量Θi沿其所有可實現(xiàn)值ˉθi方向(絕對)體積變化的期望值.

在此理解的基礎上,記第i個隨機變量Θi向其均值μi坍塌前Θ的聯(lián)合概率密度函數(shù)為pΘ(θ;μi,σi)、坍塌后為pΘi(θ;μi,0),則由均值定理有[24]

這里Fψ,i為基于Fr′echet 導數(shù)的整體靈敏度指標的參數(shù)化表達.注意到Fψ,i為一線性算子,式(48)可進一步寫為

式中不等式右邊兩項分別對應點pΘi和pΘ處的關于標準差參數(shù)的重要性測度.在點pΘ處有

其中Si(σi)為由式(31)給出的第i個隨機變量Θi關于其標準差參數(shù)σi基于Fr′echet 導數(shù)的整體靈敏度指標的重要性測度.

3 數(shù)值方法

基于Fr′echet 導數(shù)的整體靈敏度計算分為兩部分:(1)系統(tǒng)響應的概率密度函數(shù)關于輸入分布參數(shù)的導數(shù);(2)輸入概率密度函數(shù)關于分布參數(shù)的導數(shù)的范數(shù).對于第(2)部分,由于?pΘ(θ;ζ)是已知的,因此計算該范數(shù)項是容易的;而對于第(1)部分則可以通過中心差分法對其進行數(shù)值近似,即

對一些簡單隨機系統(tǒng)而言,可通過顯式解析表達的Frobenius-Perron 算子[31]直接求解式(14),從而獲得系統(tǒng)響應的概率密度函數(shù).然而,對絕大多數(shù)實際工程問題而言,Frobenius-Perron 算子難以直接進行理論求解.因此,本文將基于概率密度演化理論(probability density evolution method,PDEM)[21]進行求解,從而獲得pX(x),詳見第3.1 節(jié).

一般而言,對j=1,2,···,s分別進行概率密度演化分析即可獲得基于Fr′echet 導數(shù)的整體靈敏度指標.這總計需要2Ns次確定性分析.具體而言,不妨考慮Θ為10 維獨立的隨機向量、且各隨機變量均為兩參數(shù)分布.一般取概率密度演化分析中的概率空間剖分數(shù)為500,則總共需要求解20000 次確定性物理系統(tǒng)！這毫無疑問地給實際工程問題分析帶來了極大的挑戰(zhàn).最近,Chen 和Wan[17]從概率測度變換(change of probability measure,COM)的角度提出了降低確定性分析次數(shù)的方法,詳見第3.2 節(jié).

3.1 隨機系統(tǒng)不確定性傳播的概率密度演化理論

基于概率密度演化理論,對系統(tǒng)(10)的數(shù)值求解算法的基本步驟包括[21]:

(1)概率空間剖分與點集優(yōu)選:記Θ的樣本空間為?Θ,將其剖分為N個子域?1,?2,···,?N,并使其滿足且對任意1≤p≠q≤N有?p∩?q=?.對q=1,2,···,N,在每個子域中選取一個代表點,記為θq∈?q,并計算其對應的賦得概率.賦得概率定義為子域?q相對于全域?Θ的Voronoi體積,即

此步驟所涉及的點集優(yōu)選策略及其點集誤差分析可進一步見文獻[32-33].

(2)對q=1,2,···,N,給定代表點θ=θq并求解系統(tǒng)(10)得到˙X(θq,t)=G(X(θq,t),θq,t).

(3)對q=1,2,···,N,求解半離散化的廣義概率密度演化方程

需要說明的是,若考慮系統(tǒng)(10)的極值分布pZ(z),其中Z=max{X(t),t∈[0,T]}.此情況下,基于等價極值事件的基本思想、通過構造虛擬隨機過程即可獲得該極值分布pZ(z),詳見文獻[34].

3.2 隨機系統(tǒng)不確定性傳播的概率測度變換

3.3 基于Fre′chet 導數(shù)的整體靈敏度指標的概率密度演化概率測度變換算法

基于Fr′echet 導數(shù)的整體靈敏度指標的數(shù)值計算流程如下:

4 算例分析

4.1 算例一:正態(tài)隨機變量的線性組合函數(shù)

考慮一隨機系統(tǒng)

4.1.1 基于Fr′echet 導數(shù)的整體靈敏度指標

記相對于分布參數(shù)μ1,μ2和μ3的整體靈敏度指標分別為F (μ1),F (μ2)和F (μ3);相對于分布參數(shù)σ1、σ2和σ3的整體靈敏度指標分別為F (σ1),F (σ2)和F (σ3).本例中已知Θ1,Θ2和Θ3為正態(tài)隨機變量,因此其ε-等效分布即為正態(tài)分布.基于Fr′echet 導數(shù)的整體靈敏度指標繪于圖3 和圖4 之中,對應的重要性測度如表1 所示,按其大小排序有:μ3>μ2=μ1以及σ3>σ2=σ1.從表中可見,所有的重要性測度均滿足式(35),即其值不超過1.

圖3 相對于均值參數(shù)的基于Fr′echet 導數(shù)的整體靈敏度指標(算例一)Fig.3 Fr′echet derivative based GSI in terms of the mean values(Example 1)

圖4 相對于標準差參數(shù)的基于Fr′echet 導數(shù)的整體靈敏度指標(算例一)Fig.4 Fr′echet derivative based GSI in terms of the standard deviations(Example 1)

表1 基于Fr′echet 導數(shù)的整體靈敏度指標的重要性測度(算例一)Table 1 Importance measure of the Fr′echet derivative based GSI(Example 1)

在式(57)所示的隨機系統(tǒng)中,各輸入隨機變量關于系統(tǒng)輸出是線性且獨立同分布的關系.因此,線性系數(shù)反映了各個輸入隨機變量與系統(tǒng)響應的物理關系.具體而言,分別就均值參數(shù)和標準差參數(shù)考察Θ1,Θ2和Θ3的重要性測度(表1)如下.

除了通過表1 計算所得的重要性測度外,直接觀察圖3 和圖4 中整體靈敏度指標的峰值點亦能得到上述相關結論,即Θ1,Θ2和Θ3的重要性測度關于均值參數(shù)有比例關系、而關于標準差參數(shù)有比例關系1:1:2.

由式(38)所定義的方向靈敏度繪于圖5 和圖6 中.從圖5 可見,F (μ1)與F (μ3)的方向是相同的、與F (μ2)的方向則是相反的.而在圖6 中,F (σ1)、F (σ2)和F (σ3)的方向一致.主要原因如下:

(1)均值參數(shù)反映了隨機變量的總體趨勢.這意味著,若μ1或μ3有一個正增量?μ,則在x<0 范圍內的概率密度函數(shù)會相應地減小,而在x≥0 范圍內的概率密度函數(shù)則會相應地增大,最終導致系統(tǒng)響應的概率密度函數(shù)整體向右移動;反之則反.這一結論,正是由物理方程(57)中Θ1,Θ2和Θ3前的符號“+”、“?”和“+”所決定的.

圖5 相對于均值參數(shù)的基于Fr′echet 導數(shù)的整體靈敏度指標的方向靈敏度(算例一)Fig.5 Direction sensitivity of the Fr′echet derivative based GSI in terms of the mean values(Example 1)

圖6 相對于標準差參數(shù)的基于Fr′echet 導數(shù)的整體靈敏度指標的方向靈敏度(算例一)Fig.6 Direction sensitivity of the Fr′echet derivative based GSI in terms of the standard deviations(Example 1)

(2)對于標準差參數(shù)而言,由圖6 可知,增大σ1,σ2或σ3均會使得x∈(?2,+2)范圍內的概率密度函數(shù)相應地減小,而x∈(?∞,?2]∪[+2,+∞)范圍內的概率密度函數(shù)會相應地增大,最終導致系統(tǒng)響應的概率密度函數(shù)變得更為扁平.換言之,本例中系統(tǒng)響應標準差將隨著σ1,σ2或σ3的增大而增大.顯然,由可知這一結論是正確的.

4.1.2 與矩獨立靈敏度指標的關系

考慮隨機系統(tǒng)(57)的矩獨立靈敏度指標[10],對i=1,2,3,si(Θi=μi)分別為

需注意的是,雖然本例中分布參數(shù)μ1和μ2的重要性測度相同,但它們的方向靈敏度是相反的.可見,單憑重要性測度提供的靈敏度信息是不夠充分的,這側面反映了基于Fr′echet 導數(shù)的整體靈敏度指標的優(yōu)越性.

4.2 算例二:隧道錨桿支護系統(tǒng)

考慮一隧道錨桿支護系統(tǒng)[27,36],如圖7 所示,假定系統(tǒng)是軸對稱與各向同性的.在巖體壓力Prock、噴漿混凝土襯砌承壓力Pshot與錨桿承壓力Pbolt的共同作用下,該支護系統(tǒng)的功能函數(shù)可表達為[36]

圖7 隧道錨桿支護系統(tǒng)Fig.7 Rock bolting system of the tunnel

式中,巖體壓力Prock可隱式表達為[36]

其中P0為地應力,為修正后巖體有效黏聚力,有

噴漿混凝土襯砌承壓力Pshot為徑向變形的函數(shù)

其中,Ks為支護剛度模量(考慮d?r0)

式(62)中u(·)為徑向變形場,它是離隧道中心距離r的反比例函數(shù),即

錨桿承壓力Pbolt按下式計算

式(59)～式(65)中各確定性模型參數(shù)取值及其物理意義詳見表2.

此外,考慮式(59)～式(65)中的地應力P0(MPa)由四階矩正態(tài)變換方法進行表征[27],其前四階矩分別為μP0=9.975,σP0=2.300,αP0=?0.249 1和 βP0=2.874 8;巖體彈性模量E(GPa)服從對數(shù)正態(tài)分布,其均值為μE=4.370,標準差為σE=0.524 4;初始徑向位移u0(cm)服從對數(shù)正態(tài)分布,其均值為μu0=3.200,標準差為σu0=0.400.通常,對數(shù)正態(tài)分布的概率密度函數(shù)為

表2 隧道錨桿支護系統(tǒng)確定性參數(shù)(算例二)Table 2 Deterministic parameters in the rock bolting system of the tunnel(Example 2)[36]

因此,實際分析中對數(shù)正態(tài)分布函數(shù)的參數(shù)a和b按下式進行換算[18]

其中μ為均值、σ 為標準差.記E的分布參數(shù)為aE和bE,u0的分布參數(shù)為au0和bu0,由式(67)計算得到.本例靈敏度分析結果繪于圖8 和圖9 之中.

圖8 相對于參數(shù)P0前四階矩的基于Fr′echet 導數(shù)的整體靈敏度指標(算例二)Fig.8 Fr′echet derivative based GSI in terms of the first fourth-moments of the parameter P0(Example 2)

圖9 相對于參數(shù)E 和u0的分布參數(shù)a 和b 的基于Fr′echet導數(shù)的整體靈敏度指標(算例二)Fig.9 Fr′echet derivative based GSI in terms of the distribution parameters a and b of the parameters E and u0(Example 2)

本例中的8 個分布參數(shù)的重要性測度排序為

對應的重要性測度分別為0.910 5,0.793 8,0.713 6,0.635 3,0.275 9,0.161 1,0.122 5 和0.086 2,且均滿足不超過1 的性質.可見,地壓力P0對功能函數(shù)的影響是最大的,次之為巖體彈性模量E和初始徑向位移u0.另一方面,從圖8 和圖9 可以進一步觀察到:

(1)相較于低階矩而言,高階矩參數(shù)對響應概率密度函數(shù)的影響更為復雜,即其方向靈敏度有更多的不相交正值域或負值域.

(2)若僅考慮均值分布參數(shù),可以看到參數(shù)P0將使得g的概率密度函數(shù)向正向移動,而參數(shù)E則使g的概率密度函數(shù)向負向移動.這一點可以在錨桿承壓力計算式(65)與徑向變形場式(64)中得以理解:當P0增大或E減小時,徑向變形增大,從而使得錨桿承壓力Pbolt增大,最終導致g增大,因而其概率密度函數(shù)向正向移動.

4.3 算例三:大壩下穩(wěn)態(tài)受限滲流模型

考慮一個大壩下穩(wěn)態(tài)受限滲流分析模型[14,37],如圖10 所示.在分析中,假定滲流是各向同性和均勻的、且滲流僅發(fā)生在AB段和CD段.記粉質碎石層的透水率為Tg,x(x方向)和Tg,y(y方向);粉質砂土層的透水率為Ts,x(x方向)和Ts,y(y方向).記AB段的水頭差為hW=hD+20,則AB段到CD段的滲流控制方程為[14]

圖10 大壩下穩(wěn)態(tài)受限滲流模型的原理圖[14]Fig.10 Schematic representation of the steady-state confined seepage below the dam [14]

本文用有限元法(共計3413 個節(jié)點,1628 個單元)求解式(69)所示控制方程[14]得到分析區(qū)域的水頭差分布函數(shù)hW(x,y),最后計算滲流量q

其中q的單位為L/(h·m).式(69)各參數(shù)取值列于表3 中,其中U(7,10)為區(qū)間[7,10]上均勻分布,LN(μ,σ)為均值為μ、標準差為σ 的對數(shù)正態(tài)分布.

式中B(ah,bh)=Γ(ah)Γ(bh)/Γ(ah+bh),Γ(·)為Gamma 函數(shù).

表3 大壩下穩(wěn)態(tài)受限滲流模型參數(shù)(算例三)Table 3 Parameters in the dam model of the steady-state confined seepage(Example 3)[37]

靈敏度分析結果繪于圖11 與圖12 中.由表4 可知,按參數(shù)a進行重要性測度排序有

而按參數(shù)b進行重要性測度排序有

這意味著:(1)粉質砂土層透水率是非常重要的參數(shù)、而粉質碎石層透水率則是相對而言重要性較低;(2)除在hD中二者接近外,參數(shù)a相較于參數(shù)b更為重要.但同樣所有變量的重要性測度也均不超過1.

圖11 相對于參數(shù)a 的基于Fr′echet 導數(shù)的整體靈敏度指標(算例三)Fig.11 Fr′echet derivative based GSI in terms of the parameter a(Example 3)

圖12 相對于參數(shù)b 的基于Fr′echet 導數(shù)的整體靈敏度指標(算例三)Fig.12 Fr′echet derivative based GSI in terms of the parameter b(Example 3)

表4 基于Fr′echet 導數(shù)的整體靈敏度指標的重要性測度(算例三)Table 4 Importance measure of the Fr′echet derivative based GSI(Example 3)

4.4 算例四:鋼筋混凝土平面框架結構

考慮一榀鋼筋混凝土平面框架結構,見圖13.結構底層柱截面為600 mm×800 mm、二層柱截面為500 mm×700 mm,其他層柱截面為500 mm×600 mm;中跨梁截面統(tǒng)一設計為300 mm×450 mm、邊跨梁截面均采用300 mm×600 mm.截面配筋信息詳見文獻[17].采用OpenSEES 軟件對其進行有限元建模,具體而言:(1)單元類型設置為纖維梁單元;(2)采用與混凝土結構設計規(guī)范[38]對應的混凝土單軸損傷本構模型(“ConcreteD”命令)模擬混凝土材料的應力?應變關系;(3)考慮各向同性強化的Giuffe-Menegotto-Pinto模型[17]模擬鋼筋材料的應力?應變關系(“Steel02”命令).

分析中,考慮第1 層,第2 層,第3 ～6 層及第7 ～10 層的混凝土材料的抗壓強度fc1,fc2,fc3和fc4為獨立同分布的隨機變量,其概率信息詳見表5.注意到混凝土彈性模量Ec(104MPa)與其抗壓強度fc(MPa)是相關的.本例中,采用混凝土結構設計規(guī)范中建議的經(jīng)驗關系[38],即

圖13 鋼筋混凝土平面框架結構[17]Fig.13 Reinforced concrete frame structure[17]

表5 混凝土抗壓強度參數(shù)(算例四)Table 5 Parameters in the reinforced concrete frame structure(Example 4)

考察該結構在南北向El Centro 地震波作用下的隨機動力響應,地震動幅值調整為中震作用下對應的0.2 g.記感興趣量為頂層位移響應dtop(t)的極值,即其中為∞范數(shù).

表中,N(μ,σ)代表正態(tài)分布,μ 為均值,σ 為標準差,并對i=1,2,3,4,記fci的均值參數(shù)為μi,標準差參數(shù)為σi.

本例分析結果如圖14 和圖15 所示,所得整體靈敏度指標的重要性測度列于表6 之中.同樣,計算所得的各重要性測度均不超過1.可見,與fc2和fc3相比,fc1和fc4對頂層位移響應極值的影響更大,且尤以fc4影響最大.有意思的是,對底層位移做靈敏度分析,即取dext=∥dbottom(t)∥∞,此時相對于均值與標準差參數(shù)的重要性測度如表7 所示.此時,底層位移響應極值幾乎完全由fc1所影響.表8 給出了在0.1g峰值加速度下底層位移與頂層位移的重要性測度(此時結構響應基本處于線性狀態(tài)).值得注意的是,此時除fc1外,fc3和fc4均對底層位移響應極值的重要性測度亦有不可忽略的影響.

圖14 相對于均值的基于Fr′echet 導數(shù)的整體靈敏度指標(算例四)Fig.14 Fr′echet derivative based GSI in terms of the mean value(Example 4)

圖15 相對于標準差的基于Fr′echet 導數(shù)的整體靈敏度指標(算例四)Fig.15 Fr′echet derivative based GSI in terms of the standard deviation(Example 4)

表6 基于Fr′echet 導數(shù)的整體靈敏度指標的重要性測度(算例四/頂部位移/0.2g)Table 6 Importance measure of the Fr′echet derivative based GSI(Example 4/top displacement/0.2g)

表7 基于Fr′echet 導數(shù)的整體靈敏度指標的重要性測度(算例四/底部位移/0.2g)Table 7 Importance measure of the Fr′echet derivative based GSI(Example 4/bottom displacement/0.2g)

表8 基于Fr′echet 導數(shù)的整體靈敏度指標的重要性測度(算例四/0.1g)Table 8 Importance measure of the Fr′echet derivative based GSI(Example 4/0.1g)

從定性上看,這可能是由于結構在中震作用下進入了一定程度的非線性[11],且非線性主要表現(xiàn)在混凝土材料中,如圖16 所示.這也再次表明隨機系統(tǒng)的整體靈敏度與物理系統(tǒng)的復雜程度密切相關,且在非線性與隨機性耦合情況下,靈敏度指標排序可能出現(xiàn)與線性情況下不同的結果.

圖16 某底層柱混凝土材料的應力?應變曲線(算例四)Fig.16 Stress-strain curves of concrete materials of one bottom column(Example 4)

5 結論

基于隨機系統(tǒng)中不確定性傳播的泛函觀點及其泛函表達,引入了基于Fr′echet 導數(shù)的整體靈敏度指標.進一步,通過定義ε-等效分布,討論了該整體靈敏度的參數(shù)化表達形式,定義了其對應的重要性測度與方向靈敏度,并論述了其物理意義.結合概率密度演化理論(PDEM)和概率測度變換(COM),具體給出了該整體靈敏度對應的數(shù)值計算方法.通過4 個算例對該整體靈敏度的工程實用性進行了分析與論證.本文主要結論如下:

(1)與傳統(tǒng)整體靈敏度不同的是,本文所提整體靈敏度可看作是靈敏度分析從確定性系統(tǒng)到隨機性系統(tǒng)的一個自然延拓,因而更具有物理指導意義.

(2)特別地,本文所提整體靈敏度不僅能反映系統(tǒng)響應空間中任一點處的靈敏度,亦能通過其重要性測度給出宏觀的重要性排序、或通過其方向靈敏度標定特定的設計區(qū)域.而且,通過理論證明重要性測度不超過1.因此,該整體靈敏度所包含的不確定性量化信息十分豐富.

(3)概率密度演化理論與概率測度變換具有良好的計算精度與分析效率,可用于解決實際工程的靈敏度分析問題.

此外,如何根據(jù)實測數(shù)據(jù)進行靈敏度分析、如何確保高維概率空間下PDEM-COM 方法的計算精度、以及如何考慮非精確概率模型下的靈敏度分析等一系列問題,還有待進一步深入研究.