萬志強，陳建兵

(同濟大學(xué)土木工程防災(zāi)國家重點實驗室，土木工程學(xué)院，上海 200092)

土木工程中存在諸多不確定性因素。例如，混凝土抗壓強度一般存在15%~20%的變異性[1]；巖土體抗壓強度參數(shù)，如黏聚力和內(nèi)摩擦角，可存在20%甚至更大的不確定性[2-3]；鋼筋開始銹蝕的時間，也存在顯著的隨機性[4]。合理表征與量化上述不確定性，是進(jìn)行工程設(shè)計決策的基礎(chǔ)。

通常，工程中的不確定性可分為固有不確定性與認(rèn)知不確定性兩大類。前者是由事物的客觀屬性決定的，這種不確定性無法降低；而后者則是由知識不完備導(dǎo)致的，這種不確定性將隨著人們對問題認(rèn)識的深入和獲得信息的增加而降低[5-6]。對于固有不確定性，目前普遍采用概率論方法，即用確定性概率密度或概率分布函數(shù)描述該不確定性量[7]。而對于認(rèn)知不確定性，同樣可用概率的方法對其進(jìn)行描述，即考慮該概率密度或分布函數(shù)是不確定的，例如其分布類型或分布參數(shù)是不定的[8-10]。Bootstrap抽樣技術(shù)[2,11]和Bayes更新方法[6]為這種分布類型及分布參數(shù)的不確定性提供了可行途徑。

與此相應(yīng)，在工程設(shè)計中需要反映兩類不確定性在系統(tǒng)中的傳播及其對可靠性的影響。數(shù)十年來，針對固有不確定性的傳播與分析問題，人們發(fā)展了Monte Carlo方法[12]、隨機攝動方法[13]、正交多項式展開理論[14]和概率密度演化理論[15]等多種方法。近年來，為了進(jìn)一步提升效率，各種代理模型方法的研究得到了高度重視[16]。工程實際問題中，固有不確定性與認(rèn)知不確定性往往是耦合的。上述方法原則上可以用于同時存在兩類不確定性的問題。例如，唐小松等[2]基于Bootstrap方法、采用Monte Carlo模擬對由于數(shù)據(jù)稀缺導(dǎo)致存在認(rèn)知不確定性時的邊坡可靠度進(jìn)行了分析。但在上述方法中，通常需要兩重迭代，因而計算工作量極大[2,8-9]。最近，Chen和Wan[17]結(jié)合概率密度演化理論與概率測度變換(PDEM-COM)，提出了統(tǒng)一處理兩類不確定性的邏輯框架，對于具有認(rèn)知不確定性和固有不確定性耦合的不確定性量化與傳播問題發(fā)展了高效分析策略。

在上述研究工作的基礎(chǔ)上，本文基于工程實測數(shù)據(jù)，首先分別引入Bootstrap方法和Bayes更新方法進(jìn)行認(rèn)知不確定性表征。在此基礎(chǔ)上，結(jié)合PDEM-COM方法，給出了基于數(shù)據(jù)進(jìn)行工程系統(tǒng)不確定性量化、傳播與可靠性分析的基本途徑與數(shù)值方法。實例分析表明，本文建議方法具有很高的精度和效率，便于實際工程應(yīng)用。

1 兩類不確定性的表征方法

不失一般性，記系統(tǒng)響應(yīng)為：

式中：Θ為模型的不確定性參數(shù)；Λ為確定性參數(shù)；f(·)為線性或非線性算子。進(jìn)一步，記具有固有不確定性的參數(shù)為ΘAU，具有認(rèn)知不確定性的參數(shù)為ΘEU，則有：Θ=(ΘAU,ΘEU)T，或

即有NAU個固有不確定性參數(shù)、NEU個認(rèn)知不確定性參數(shù)。

在本文中，考慮一類重要的認(rèn)知不確定性來源，即固有不確定性參數(shù)ΘAU的概率分布具有不確定性[1,8]。例如，在巖土工程中，由于試驗數(shù)據(jù)十分有限，導(dǎo)致邊坡可靠度分析所需的工程抗剪強度參數(shù)的分布參數(shù)變異性極大，這是分布參數(shù)不確定性的體現(xiàn)[2-3]；在結(jié)構(gòu)工程中，混凝土抗壓強度的試驗數(shù)據(jù)可能同時通過正態(tài)分布和Weibull分布的假設(shè)檢驗[18]，這是分布類型不確定性的表現(xiàn)。上述概率分布的不確定性，可通過Bootstrap方法和Bayes更新方法進(jìn)行表征和量化，詳述如下。

1.1 Bootstrap方法

Bootstrap方法[10]是對原始試驗數(shù)據(jù)進(jìn)行有放回隨機抽樣、以獲取分布參數(shù)不確定性信息的方法。當(dāng)實測數(shù)據(jù)較為稀缺時，該方法是量化認(rèn)知不確定性的有效途徑[2]。

例如，假定原始試驗數(shù)據(jù)為D=(θ1,θ2,…,θn)，進(jìn)行K次Bootstrap抽樣后可得：

若僅根據(jù)原始數(shù)據(jù)D進(jìn)行統(tǒng)計分析，可以得到前二階矩估計值，即均值μΘ和標(biāo)準(zhǔn)差σΘ。通過Bootstrap抽樣技術(shù)，可進(jìn)一步獲得一系列矩估計值當(dāng)隨機變量Θ的分布類型確定時，Bootstrap提供了表征分布參數(shù)認(rèn)知不確定性的途徑。此時，隨機變量Θ的概率描述不再是唯一的概率密度函數(shù)pΘ(θ,μΘ,σΘ)，而變?yōu)楦怕拭芏群瘮?shù)族相應(yīng)地，工程結(jié)構(gòu)或系統(tǒng)分析模型的輸出量X不再具有唯一確定的概率密度函數(shù)，而是對應(yīng)的概率密度函數(shù)族在第2節(jié)中，將詳細(xì)闡述如何在不顯著增加計算成本的前提下，進(jìn)行上述不確定性的量化。

1.2 Bayes更新方法

另一種可能的工程情況是：在結(jié)構(gòu)設(shè)計階段，通過現(xiàn)場實測得到一批試驗數(shù)據(jù)，需要對設(shè)計參數(shù)進(jìn)行更新[19]；或在工程監(jiān)測時獲得一批新的實測數(shù)據(jù)，需要對設(shè)計階段的分析結(jié)果進(jìn)行驗算。這兩種情況可統(tǒng)一為如下表述：

記隨機變量Θ的原始觀測數(shù)據(jù)集為D1、新增觀測數(shù)據(jù)集為D2。同時，記由數(shù)據(jù)D1通過統(tǒng)計推斷得到的概率密度函數(shù)為其中ζ1為由數(shù)據(jù)D1統(tǒng)計推斷得到的分布參數(shù)。若考慮Θ的推斷過程中存在認(rèn)知不確定性，即ζ1本身是不確定性的。在一般情況下，ζ1的先驗分布pζ(ζ)是人為假定的。事實上，根據(jù)1.1節(jié)的討論，該先驗分布可基于Bootstrap方法給出。當(dāng)引入新增數(shù)據(jù)D2后，通過Bayes更新方法，可得參數(shù)的后驗分布：

其中，L(·)為極大似然函數(shù)，有：

式中，card(·)表示集合中的元素個數(shù)。由式(4)和式(5)可得：

由此，隨機變量Θ的概率密度函數(shù)由pΘ(θ|ζ1)更新為pΘ(θ|ζ2)。此時，工程模型輸出對應(yīng)的概率分布亦應(yīng)由更新為這一更新過程，可通過下節(jié)介紹的高效方法完成。

2 概率密度演化與概率測度變換

本節(jié)將簡述概率密度演化理論(PDEM)與概率測度變換(COM)的基本原理。在此基礎(chǔ)上，介紹固有不確定性與認(rèn)知不確定性同時存在的情況下不確定性量化與傳播的PDEM-COM統(tǒng)一理論框架。

2.1 概率密度演化理論(PDEM)

不失一般性，考慮隨機微分方程：

其中，pXΘ(x,θ,t)為(X,Θ)的聯(lián)合概率密度函數(shù)。相應(yīng)的初始條件為：

式中，δ(·)為Dirac函數(shù)。求解式(8)并對Θ積分可得：

求解方程組式(7)~式(9)的數(shù)值算法如下[20]：

1) 對概率空間ΩΘ進(jìn)行剖分，記剖分子域為：

由此，可獲得對應(yīng)的代表點Θ=θq,q=1,2,…,nsel及其對應(yīng)的賦得概率：

其中，nsel為代表點個數(shù)。本文采用使得GF-偏差最小化[21]的點集優(yōu)選策略獲得上述代表性點集。

2) 對每一個代表點求解物理方程式(7)，獲得速度狀態(tài)量

3) 對每一個代表點求解如下概率密度演化方程：

4) 累計求和：

由此，一旦給定一組系統(tǒng)輸入Θ的概率信息pΘ(θ)，通過上述概率密度演化理論可十分方便地獲得相應(yīng)的系統(tǒng)輸出概率信息pX(x)。

2.2 概率測度變換(COM)

考慮μ和ν是σ-有限測度，且定義于可測空間(ΩΘ,F)；若ν相對于μ絕對連續(xù)，則存在可測函數(shù)T :ΩΘ→[0,∞)，使其滿足[22]：

式中，函數(shù)T=dν/dμ稱為Radon-Nikodym導(dǎo)數(shù)。換言之，不同測度之間可以通過算子T進(jìn)行變換，而不需要引入額外計算工作量。

基于此，當(dāng)同時存在固有不確定性與認(rèn)知不確定性時，Chen和Wan[17]結(jié)合概率密度演化理論(PDEM)與概率測度變換(COM)，提出了一類不確定性量化與傳播的高效方法(簡記為PDEM-COM)，其算法如下：

1) 給定一組概率信息pΘ(θ)，進(jìn)行一輪PDEM分析。分析時，儲存代表性點集及速度狀態(tài)量和pX(x,t)。當(dāng)研究的是靜力分析問題時，可采用等價極值分布方法[23]。

4) 重復(fù)步驟2)~3)，直至所有需要重新量化的認(rèn)知不確定性完成不確定性傳播過程。

注意到，在PDEM-COM中，僅在第一輪概率密度演化分析中進(jìn)行確定性物理過程分析(一般為300次~500次左右)。在隨后的認(rèn)知不確定性量化分析過程中，將集中于測度變換的計算(新的賦得概率)。此時，與復(fù)雜物理系統(tǒng)確定性分析的計算成本相比，新增的計算工作量很小。而且，物理系統(tǒng)越復(fù)雜，這一優(yōu)勢越明顯。

3 工程算例

3.1 無限邊坡穩(wěn)定分析模型

考慮一個無限邊坡穩(wěn)定分析模型[2]，如圖1所示，其極限狀態(tài)方程為：

取邊坡幾何參數(shù)H=5m和α=30°，土體重度γW=17kN/m3?？辜魪姸葏?shù)c和φ分別代表土體有效黏聚力和有效內(nèi)摩擦角，分別服從對數(shù)正態(tài)分布和正態(tài)分布。

抗剪強度參數(shù)的試驗數(shù)據(jù)往往不易獲得，如文獻(xiàn)[2]僅給出了24組觀測值。唐小松等[2]根據(jù)這24組觀測值，采用Bootstrap方法考察了抗剪強度參數(shù)分布的兩類認(rèn)知不確定性，并給出了邊坡可靠度在這兩類不確定性下的可靠度指標(biāo)量化范圍。文獻(xiàn)[2]同時還指出，該分析需要與Bootstrap抽樣次數(shù)等量的可靠度計算次數(shù)，因此對于復(fù)雜的邊坡可靠度分析問題而言，這種需要進(jìn)行雙重循環(huán)的認(rèn)知不確定性分析方法計算工作量巨大。

根據(jù)第2節(jié)的基本理論，對認(rèn)知不確定性分析，可采用如下流程：

1) 令原始數(shù)據(jù)集為D=(ci,φi),i=1,2,…,n，本例取n=24；定義Bootstrap子樣本為Dj,j=1,2,…,B，其中B為Bootstrap抽樣次數(shù)，本例取B=104。

2) 由Bootstrap方法對24組觀測值進(jìn)行認(rèn)知不確定性度量，即給出抗剪強度參數(shù)c和φ的均值μc和μφ、標(biāo)準(zhǔn)差σc與σφ的Bootstrap分布及其95%置信區(qū)間(表1)，下橫線對應(yīng)區(qū)間下限值，上橫線對應(yīng)區(qū)間上限值，即：

3) 基于概率密度演化理論，以下列分布參數(shù)進(jìn)行一次可靠度分析[16,20-21,23]：

4) 針對每一個Bootstrap子樣本Dj計算統(tǒng)計量并計算對應(yīng)的測度變換算子

5) 由3)得到的概率密度演化分析結(jié)果，結(jié)合4)計算得到的測度變換算子，給出每個Bootstrap子樣本對應(yīng)的概率密度演化分析結(jié)果。這里的概率密度演化分析結(jié)果，可以是響應(yīng)的一階、二階矩信息，亦可是可靠度指標(biāo)信息，由問題本身確定。在本例中，考慮邊坡可靠度指標(biāo)，即：

式中，Φ(·)為標(biāo)準(zhǔn)正態(tài)分布的分布函數(shù)。

由此可得到存在認(rèn)知不確定性的情況下的邊坡可靠度指標(biāo)β={β(j),j=1,2,…,B}及其5%分位值和95%分位值，如圖2所示。為驗證本方法的計算精度，采用500萬次Monte Carlo模擬求解本例的可靠度分布，如圖3所示。從圖2與圖3的對比中可見，不僅可靠指標(biāo)分布直方圖總體形狀一致，而且相應(yīng)的均值與分位值均具有很高的計算精度。特別指出，在PDEM-COM中，僅需要進(jìn)行一輪完整的概率密度演化分析，因而全部分析中僅需要進(jìn)行500次確定性分析。

表1 由Bootstrap方法給出的95%置信區(qū)間(例3.1)Table 1 95% confidence interval by Bootstrap (E.3.1)

圖1 無限邊坡穩(wěn)定分析模型[2]Fig.1 Infinite slope model for stability analysis[2]

圖2 邊坡可靠度指標(biāo)的Bootstrap分布(PDEM-COM)Fig.2 Bootstrap distribution of reliability index of infinite slope model (PDEM-COM)

圖3 邊坡可靠度指標(biāo)的Bootstrap分布(5×106 MCS)Fig.3 Bootstrap distribution of reliability index of infinite slope model (5×106 MCS)

3.2 擋土墻穩(wěn)定性分析模型

考慮擋土墻穩(wěn)定性分析模型[3]，如圖4所示。在承載力滿足的情況下，該模型可能存在兩種失效模式，其對應(yīng)的功能函數(shù)分別為：

1) 滑移破壞模式

式中：cint=1/3·cb和φint=2/3·φb分別為接觸面處的粘聚力和內(nèi)摩擦角，cb和φb是對應(yīng)地基處的粘聚力與內(nèi)摩擦角；Pa為主動土壓力(忽略墻體表面裂縫影響)，有：

式中，Ka為主動土壓力系數(shù)，由Rankine理論給出：

圖4 擋土墻穩(wěn)定性分析模型[3]Fig.4 Stability analysis model of retaining wall[3]

式(20)~式(22)其他參數(shù)取值及定義詳見表2。

表2 擋土墻穩(wěn)定性分析模型參數(shù)(例3.2)Table 2 Model parameters of retaining wall (E.3.2)

2) 傾覆破壞模式

其中，L1、L2和La分別取為：

顯然，兩種失效模式均為有效粘聚力c與有效內(nèi)摩擦角φ的函數(shù)。文獻(xiàn)[3]提供了63組固結(jié)排水情況下的抗剪強度參數(shù)試驗數(shù)據(jù)，本節(jié)將利用該數(shù)據(jù)考察認(rèn)知不確定性對擋土墻穩(wěn)定性分析的影響。

3.2.1 情況1：由Bootstrap方法度量認(rèn)知不確定性

首先考察數(shù)據(jù)樣本大小帶來的認(rèn)知不確定性。記該63組數(shù)據(jù)的前30組數(shù)據(jù)為D1(數(shù)據(jù)集1)，全部數(shù)據(jù)為D2(數(shù)據(jù)集2)。同樣地，由Bootstrap方法可給出抗剪強度參數(shù)c和φ的均值μc和μφ及標(biāo)準(zhǔn)差σc和σφ的Bootstrap分布及其95%置信區(qū)間，見表3所示。將8組置信區(qū)間繪制于圖5可見，在樣本量較小時置信區(qū)間范圍也較大，擴大樣本量確能縮小該不確定性范圍；二階統(tǒng)計矩的縮小程度要比一階統(tǒng)計矩大。

表3 由Bootstrap方法給出的95%置信區(qū)間(例3.2.1)Table 3 95% confidence interval by Bootstrap method(E.3.2.1)

圖5 兩組數(shù)據(jù)集給出的Bootstrap分布的95%置信區(qū)間Fig.5 95% confidence intervals of two data sets by Bootstrap method

由數(shù)據(jù)集D1和D2給出的擋土墻可靠度的Bootstrap分布如圖6和圖7所示。可以看到，當(dāng)認(rèn)知不確定性較大的時候(如數(shù)據(jù)集D1相對于數(shù)據(jù)集D2)，其可靠度指標(biāo)的變化范圍也較大。需要注意的是，當(dāng)增加統(tǒng)計樣本數(shù)量以減少認(rèn)知不確定性時，可靠度指標(biāo)的不確定性范圍雖然隨之變窄，但其期望值卻下降了。這說明，若不考慮該認(rèn)知不確定性帶來的影響，如將可靠度指標(biāo)的期望作為最終分析評估的結(jié)果，將使得分析結(jié)果偏于不安全。因此，有必要根據(jù)實際工程案例進(jìn)行認(rèn)知不確定性量化工作。

3.2.2 情況2：由Bayes更新度量認(rèn)知不確定性

首先由數(shù)據(jù)集D1給出抗剪強度參數(shù)c和φ的分布參數(shù)(先驗概率密度函數(shù))，然后通過額外數(shù)據(jù)集(算子?代表取差集)對進(jìn)行更新，可得到(后驗概率密度函數(shù))。特別地，若考慮用所有已知數(shù)據(jù)D2進(jìn)行參數(shù)識別，亦可得到對應(yīng)的(全信息概率密度函數(shù))。圖8給出上述三種情況對應(yīng)的概率密度函數(shù)。在PDEM-COM分析中，以為分析工況進(jìn)行一次概率密度演化分析即可，而的情況用測度變換方法給出，由此給出的可靠指標(biāo)如表4所示，括號中數(shù)值為采用100萬次MCS直接模擬得到的結(jié)果。從表中可見，采用PDEM-COM可以得到精度很高的可靠指標(biāo)估計結(jié)果。

圖6 由數(shù)據(jù)集D1得到的兩種失效模式的可靠度分布Fig.6 Distributions of reliability index of two failure modes by data set D1

圖7 由數(shù)據(jù)集D2得到的兩種失效模式的可靠度分布Fig.7 Distributions of reliability index of two failure modes by data set D2

圖8 三種情況下抗剪強度參數(shù)的概率密度函數(shù)Fig.8 PDFs of shear parameters on three conditions

表4 兩種失效模式的可靠度指標(biāo)(例3.2.2)Table 4 Reliability indexes of two failure modes (E.3.2.2)

3.3 屋面桁架結(jié)構(gòu)

考慮一個屋面桁架結(jié)構(gòu)[24](如圖9所示)，其功能函數(shù)定義為：

式中：q為均布荷載，簡化為三個集中力P=ql/4；Ac和As分別為混凝土構(gòu)件與鋼構(gòu)件的截面積；Ec和Es分別為混凝土材料與鋼材的彈性模量值；ΔS為屋頂豎向位移；極限變形位移取ΔR=l/400。

圖9 屋面桁架結(jié)構(gòu)[24]Fig.9 Roof truss structure[24]

本節(jié)考察認(rèn)知不確定性與固有不確定性同時存在時對結(jié)構(gòu)可靠性分析的影響。若進(jìn)行結(jié)構(gòu)彈性位移下的可靠度分析(關(guān)于桿件內(nèi)力可靠度分析可詳見文獻(xiàn)[25])，可考慮混凝土彈性模量與截面積存在不確定性，其他設(shè)計參數(shù)均為確定性的。進(jìn)一步，由于既有試驗數(shù)據(jù)有限，考慮混凝土截面積的不確定性僅為固有不確定性，即服從給定的概率分布；而混凝土彈性模量則同時存在認(rèn)知不確定性和固有不確定性，即其概率分布是不確定性的，由既有試驗數(shù)據(jù)并結(jié)合Bootstrap抽樣技術(shù)進(jìn)行表征。本例取文獻(xiàn)[26]給出的100組商品混凝土彈性模量試驗值，將其分為數(shù)據(jù)集D1(前30組數(shù)據(jù))和數(shù)據(jù)集D2(所有數(shù)據(jù))兩部分，分別進(jìn)行Bootstrap分析，數(shù)據(jù)分析結(jié)果如圖10所示。

圖10 兩組數(shù)據(jù)集給出的Bootstrap分布的95%置信區(qū)間Fig.10 95% confidence intervals of two data sets by Bootstrap method

從圖中可見，增大樣本數(shù)量降低了不確定性的變化范圍，這與圖5中的特征是一致的。本例中的結(jié)構(gòu)分析參數(shù)詳見表5和表6。

表5 屋面桁架結(jié)構(gòu)確定性參數(shù)(例3.3)Table 5 Deterministic parameters of roof truss model (E.3.3)

表6 屋面桁架結(jié)構(gòu)不確定性參數(shù)(例3.3)Table 6 Uncertain parameters of roof truss model (E.3.3)

采用極限狀態(tài)函數(shù)式(25)，基于PDEM-COM方法得到可靠度分析結(jié)果如圖11所示。注意到，在初始時，由混凝土彈性模量數(shù)據(jù)進(jìn)行Bootstrap分析得到的95%置信區(qū)間(見圖10)范圍很窄，這說明該物理參數(shù)的認(rèn)知不確定性程度較低。因此，從圖11可見結(jié)構(gòu)可靠指標(biāo)的變化范圍亦較小。

圖11 屋面桁架結(jié)構(gòu)可靠度分布Fig.11 Reliability index distributions of roof truss structure via

4 結(jié)論

針對數(shù)據(jù)稀缺與信息更新導(dǎo)致的認(rèn)知不確定性，引入Bootstrap方法和Bayes更新的方法進(jìn)行表征。進(jìn)而，針對固有不確定性與認(rèn)知不確定性耦合問題，結(jié)合概率密度演化理論與概率測度變換(PDEM-COM)，提出了從數(shù)據(jù)出發(fā)到系統(tǒng)中的不確定性量化、傳播與可靠性分析的基本途徑和數(shù)值算法，實現(xiàn)了考慮認(rèn)知不確定性的工程可靠性高效分析。通過3個工程實例分析，驗證了所提方法的有效性。主要結(jié)論如下：

(1) 采用Bootstrap方法和Bayes更新方法進(jìn)行數(shù)據(jù)稀缺與信息更新時的認(rèn)知不確定性表征，可方便地結(jié)合概率密度演化-測度變換，形成從數(shù)據(jù)出發(fā)到系統(tǒng)不確定性量化的有效途徑；

(2) 由于認(rèn)知不確定性的存在，結(jié)構(gòu)可靠度不再是一個確定值，而是一個隨機變量，其置信區(qū)間與認(rèn)知不確定性的程度密切相關(guān)；

(3) 基于PDEM-COM的不確定性傳播方法具有較高的精度和效率，可適用于工程實際問題分析。

對于認(rèn)知不確定性導(dǎo)致分布信息發(fā)生大偏差、隨機輸入下復(fù)雜工程系統(tǒng)整體動力可靠度及其不確定性量化等問題，尚需今后進(jìn)一步深入研究。