基于證據理論和代理模型的QMU分析

謝朝陽，李貴杰，彭忠明，黃洪鐘

(1. 中國工程物理研究院總體工程研究所 四川 綿陽 621000；2. 電子科技大學機械電子工程學院 成都 611731)

在結構可靠性以及系統性能分析中，需要通過不確定性傳播研究輸入參數的不確定性對系統響應的影響，獲得系統響應的不確定性。工程實際中隨機和認知不確定性同時存在，設結構系統的響應函數抽象為：

不確定性普遍存在于工程結構的分析和設計中。根據不確定性的不同來源和人們的認知程度，在工程實際研究中將不確定性分為隨機不確定性(aleatory uncertainty)和認知不確定性(epistemic uncertainty)[1-2]。傳統的結構可靠性分析主要考慮隨機不確定性，一般基于概率方法，需要大量實驗樣本構造參數的精確概率分布[3-4]。對于信息難以獲取或掌握信息較少的情況，國內外學者提出了基于模糊數學、區間分析、證據理論等非概率模型來描述處理結構可靠性分析中的認知不確定性，并發展了相應的非概率可靠性指標和分析方法[5-8]。近年來，美國圣地亞等國家實驗室提出了裕量與不確定度量化(quantification of margins and uncertainty, QMU)的概念用于具有隨機和認知不確定性的禁試武庫可靠性、安全性評估認證以及風險決策[9-10]。國內外學者圍繞隨機和認知不確定性分別研究了QMU理論及其應用于系統性能評估和分析時的具體實施流程。文獻[9]描述了如何采用QMU處理復雜系統的可靠性以及代替傳統可靠性估計的置信度，提出了QMU方法的關鍵要素即性能通道、裕量、不確定性。文獻[11]提出的QMU度量指標考慮了系統操作區域和性能要求的不確定性。文獻[12]應用QMU方法研究了環狀結構的可靠性。針對混合不確定性，文獻[13]集成概率區間分析和證據理論研究了可應用于QMU的計算工具。文獻[14]基于Bayes網絡充分考慮復雜分層系統各級組件的隨機和認知不確定性，建立了QMU方法的Bayes網絡框架。

從上述文獻看，不確定性量化與傳播是QMU方法應用于工程實際的關鍵環節，而實際工程中系統性能響應常常由隱式模型表示，需進行復雜的數值模擬計算，因此需要建立代理模型進行不確定性傳播分析從而降低計算成本。本文在QMU概念的基礎上，考慮隨機和認知混合不確定性，提出了基于證據理論和Kriging代理模型的結構性能裕量與不確定性量化計算方法，以提高不確定性分析精度與效率。

1 QMU的基本概念與分析流程

QMU是一套考慮不確定性并以性能閾值和裕量為關鍵點的復雜系統風險決策支持方法，如圖1所示，其裕量(M)的定義為系統名義響應值和性能通道邊界閾值之間的距離[9-10]。系統設計或服役結構的性能是否滿足要求通過置信因子CF(confidence factor)進行定義度量。置信因子CF定義為性能裕量M與該性能的不確定度U之比：

如果對系統評估過程中的不確定性認識量化足夠充分，當置信因子CF＞1時，就判定為可靠；當置信因子CF≤1時，就判定系統存在失效風險。

對于復雜系統或大型結構，運用QMU方法評估其性能狀態的分析流程涉及3個環節(如圖1所示)：1)建立觀測清單；2)建立性能通道；3)性能裕量預測及其不確定性量化。對于一般的結構可靠性分析，如果已知結構性能通道特征方程，即性能響應的功能函數，則采用QMU進行評估的關鍵是通過不確定性傳播計算結構響應的不確定性分布。

2 基于證據理論的混合不確定性分析

2.1 證據理論的基本原理

D-S證據理論源于Dempster和Shafer提出的一套基于“證據”和“組合”來處理不確定性推理問題的數學方法。由于證據理論能夠對各種不完全信息、不確定信息、不可靠信息甚至沖突信息進行合理的描述和處理，因而被廣泛應用于各種領域的不確定性研究。其基本原理主要包括基本可信度分配函數、信任與似然函數[15]。

圖1 QMU方法的關鍵要素流程示意圖

設Θ為識別框架，即對于某一問題或事物所能認識到的所有可能結果的集合。函數m為從Θ的冪集2Θ到[0,1]的映射且滿足以下3條性質：則函數m為識別框架Θ上的基本概率分配(BPA)。m(A)稱為集合A的基本概率分配值，表示對A的信任程度，其中m(A)＞0的集合稱為該分配函數的焦元。若B表示識別框架的任一子集，則B的真實程度由信任函數和似然函數表示：

式中，Bel(B)表示完全支持命題B的證據的基本概率函數值之和；Pl(B)是完全或者部分支持命題B的基本概率函數之和。Bel(B)和Pl(B)分別構成對B信任度的下限和上限，記為[Bel(B),Pl(B)]。

2.2 基于證據理論的不確定性傳播模型

根據證據理論，設不確定性變量Y的描述用證據空間(Y,Y,m)表示。其中Y表示不確定性變量Y的所有可能值的集合，構成Y的樣本空間。令U表示集合Y的子集，Y表示子集U的可數集合，即Y的焦元的集合。m為子集U的基本概率分配函數，滿足當U∈Y時，m(U)≥0；當U?Y時，m(U)=0且對所有U∈Y，

在結構可靠性以及系統性能分析中，需要通過不確定性傳播研究輸入參數的不確定性對系統響應的影響，獲得系統響應的不確定性。工程實際中隨機和認知不確定性同時存在，設結構系統的響應函數抽象為：

式中，X為客觀不確定性向量，采用概率分布來描述其不確定性；Y為n個證據空間(Y1,Y1,m1),(Y2,Y2,m2),…,(Yn,Yn,mn)描述的不確定性變量，其聯合證據空間可表示為(Y,Y,mY)。其中Y=Y1×Y2×…×YnY，U=U1×U2×…×UnY，U∈Y，Uk∈Yk，(k=1,2,…,nY)，mY(U)為聯合BPA，可以通過下式計算：

根據全概率理論，G小于某個閾值c的概率P可寫為：

式中，乘積空間Y中的焦元Ui的概率等于Ui的聯合BPA，即Pr{Yi∈Ui}=mY(Ui)。根據認知不確定性變量的性質，響應G小于某個閾值c的概率將會是一個區間，其邊界用信任函數和似然函數度量：

改變不同的閾值c，即可得到系統響應G的不確定性描述。

3 混合不確定性傳播求解方法

3.1 Kriging代理模型

式中，fT(x)=[f1(x),f2(x),…,fb(x)]為回歸模型的基函數，是一個確定性部分提供對模擬全局的近似，一般用x的多項式表示；α為回歸模型對應基函數的系數；S(x)是一個高斯隨機過程，提供對模擬局部偏差的近似，其均值為0，方差為σ2。通常情況下多項式f(x)可以取為平均響應值μ，而不影響模擬的精度，式(10)可改寫為：

對于空間某點x′，系統性能函數的響應可通過Kriging代理模型的預測值和預測均方差分別表示為：

式中，r為點x′和樣本點之間的相關向量；A是n×1的單位向量。

3.2 基于最大置信水平期望提高加點方法

如何確定Kriging代理模型訓練樣本的位置和數量是保證代理模型計算精度和效率的關鍵。文獻[18]計算結構可靠性時采用自適應抽樣方法，采用最大置信水平期望提高準則建立Kriging代理模型訓練樣本的加點方法。

基于已構建的Kriging代理模型，可通過MCS計算結構的響應不確定性分布。采用MCS方法通過代理模型M計算時，假設通過隨機抽樣得到N個輸入向量Xm=[xm,1,xm,2,…,xm,N]，結構響應函數在點xm,i的預測值G(xm,i,)為服從分布的隨機變量。對于樣本空間中的第i個樣本向量xm,,的結構響應預測值可按其是否小于或大于閾值c進行分類。針對該點的預測置信水平CL(xm,i)可根據正確分類的概率表示如下：

式中，Φ(?)表示標準正態分布函數；|?|為求絕對值符號。式(14)表示在響應預測值的條件下，實際真實值G(xm,i)＜c的概率。對于由MCS產生的N個樣本點，根據閾值c正確分類概率獲得的置信水平期望為：

ECL表示代理模型替代原響應函數的近似精度。如果Kriging代理模型M計算響應預測值的近似精度低于設定的目標值，則需添加新的訓練樣本點構建新的代理模型M*。設響應函數在點xm,i處的真實響應值為G(xm,i)，增加(xm,i，G(xm,i))為訓練數據后構建的代理模型M*，其近似精度的置信水平期望提高可表示為：

由于尚未獲得樣本點xm,i處的真實響應值G(xm,i)，式(16)無法直接計算獲得，因此通過置信水平定義EI計算表達式為[18]：

式中，CL(xm,i)是樣本xm,i的響應值正確分類的概率；fX(xm,i)是樣本xm,i的概率密度函數值；為響應預測值的均方差。式中第一項表示xm,i成為訓練樣本后代理模型針對該點預測置信水平提高的潛在值，前兩項相乘表示增加樣本點xm,i后獲得置信水平提高潛在值的可能性，第三項表示響應預測值的均方差對置信水平期望提高的影響。通過式(17)計算MCS樣本xm中每個點xm,i的置信水平的期望提高，其中使得該值最大的點作為下一個加入Kriging代理模型的訓練樣本點，即：

3.3 混合不確定性傳播計算流程

基于最大置信水平期望提高準則的自動加點Kriging代理模型和證據理論，以求解式(9)的似然函數為例，計算混合不確定性傳播的主要步驟如下：

1)根據證據理論，計算響應函數中主觀不確定性變量Y的聯合BPA，并確定每個焦元Ui的邊界。

2)根據隨機變量X的分布類型，采用拉丁超立方抽樣(Latin Hypercube sample)方法取其中n為隨機變量X的維數。

3)采用優化方法計算所抽取樣本點在焦元Ui內對應的響應函數的極值Gmin(xt)。

4)根據xt和Gmin(xt)建立Kriging代理模型M。

5)根據MCS抽樣和代理模型，采用MCEI準則搜尋點x*并計算Gmin(x*)，將該數據加入到原樣本中返回步驟4)，更新代理模型。

6)重復步驟4)、步驟5)，直至滿足ECL大于設定的限制。

7)根據MCS樣本計算響應函數小于閾值c的概率邊界。

8)重復步驟3)～步驟7)步計算出每個焦元Ui內響應函數的輸出不確定性小于某一閾值c的概率邊界。

9)根據聯合BPA和步驟8)步計算結果求解在隨機和認知混合不確定性輸入下響應函數的不確定性分布。

式(8)表示的信任函數則需將步驟3)～步驟5)中的求最小值Gmin(xt)替換為求最大值Gmax(xt)，即可獲得系統響應小于某一特定閾值c的概率下限。

4 基于證據理論框架的置信因子計算

從QMU的概念和流程可知，性能裕量M由系統響應和閾值的估計值描述，不確定性U也由系統響應和閾值的不確定性進行確定。在證據理論的框架下，結構性能響應的不確定性分布用似然函數和信任函數表示。引入安全系數γ，從支撐風險決策的角度保守估計，針對性能通道上邊界閾值要求定義系統響應估計值和不確定性取值見表1。下邊界閾值要求對應的系統響應估計值和不確定性取值見表2。

基于表1和表2的取值，上下邊界的性能裕量M和不確定性U均可由下式計算：

如果性能通道同時有上下邊界閾值要求，表征系統QMU度量的置信因子則為上下邊界置信因子的最小值。

表1 系統響應估計值和不確定性取值(上邊界閾值)

5 算例研究

圖2 曲柄滑塊機構示意圖

圖2為曲柄滑塊機構示意圖[19]。曲柄的長度為a，連桿的長度為b，其內徑為d1，外徑為d2，受到外力載荷為F，連桿的楊氏模量為E，連桿的屈服強度為S。它們均為隨機變量，其分布參數列于表3中。此外，滑塊與地面NN的摩擦力系數μ和偏移量e的精確分布不可獲知?；趯＜医涷灪陀邢薜臍v史數據，變量μ和e的區間和BPA可以獲得，它們的BPA及其聯合BPA分別列于表4和表5中。針對連桿的強度，傳統的結構可靠性分析定義極限狀態函數為材料的強度與最大應力之間的差，即：

文獻[19]基于證據理論采用UAA方法求解了上述極限狀態函數的非精確概率可靠性。本文根據QMU方法的框架，以連桿的最大應力作為系統的性能響應函數Y，S作為性能響應函數的上邊界閾值，采用置信因子評估結構的可靠性。

表3 曲柄滑塊機構中隨機變量X的分布參數

表4 曲柄滑塊機構算例認知不確定性變量

表5 曲柄滑塊機構認知變量的聯合BPA

由上述自動加點Kriging代理模型方法計算得到的系統響應函數的不確定性分布由信任函數和似然函數描述如圖3所示，其中虛線表示系統響應Y的信任函數和似然函數，實線表示閾值YT的累積概率分布函數。

根據表1和式(10)～式(20)提供的計算方法，取不同的安全系數γ得到的系統QMU分析結果如表6所示。根據QMU的概念，當置信因子CF＞1時，判定為可靠；當置信因子CF≤1時，判定系統存在失效風險。表6中當風險決策時的安全系數γ=0.95時，置信因子CF＜1，說明系統存在失效風險，當可接受的決策風險變大，即安全系數降低時，置信因子CF＞1，說明系統可靠為可接受狀態。算例表明安全系數的大小與QMU分析得到的置信因子密切相關，從一定程度上表示了風險決策者對系統決策風險的可接受程度。

圖3 系統響應函數和閾值的不確定性分布

表6 不同安全系數取值計算得到的CF

從結構可靠性分析的視角出發，根據蒙特卡洛方法式(21)定義的極限狀態方程，可得到響應失效概率的真實值介于信任函數與似然函數值之間，即4.68×10-4≤Pf≤1.164 8×10-3。說明滑塊存在較低的失效概率。γ=0.95時的置信因子CF能夠反映出系統的失效狀態，這與傳統可靠性分析的結果相一致。

由于本例中系統響應函數是認知不確定性變量e和μ的單調函數，函數的極值出現在兩個變量的端點處。因此不考慮優化求極值過程，采用基于Kriging方法調用系統響應函數的次數為2×9×50=900(每個焦元建立代理模型的訓練樣本數為50)，而采用MCS方法9個焦元計算信任與似然函數分布共需極限狀態函數次數為9×106。這說明本文建立的基于Kriging方法用于隨機和認知混合不確定性下的QMU分析具有很好的可行性。

6 結 束 語

本文建立的基于證據理論和代理模型能夠有效處理混合不確定性的傳播計算，從而支撐系統QMU分析過程中的不確定性量化這一關鍵要素。算例研究表明，基于最大置信水平期望提高加點的Kriging模型在處理混合不確定性分析方面具有較高的精度和效率。QMU是一套考慮不確定性并以性能閾值和裕量為關鍵點的復雜系統風險決策支持方法，算例分析表明QMU方法可作為傳統可靠性分析的補充，用于度量系統可接受風險時的可靠性分析。

[1]HELTON J C, JOHNSON J D, OBERKAMPF W L, et al.Representation of analysis results involving aleatory and epistemic uncertainty[J]. International Journal of General Systems, 2010, 39(6)： 605-646.

[2]HELTON J C. Alternative representations of epistemic uncertainty[J]. Reliability Engineering and System Safety,2004, 85(1-3)： 1-10.

[3]RACKWITZ R. Reliability analysis-a review and some perspectives[J]. Structural Safety, 2001, 23(4)： 365-395.

[4]AU S K, BECK J L. Estimation of small failure probabilities in high dimensions by subset simulation[J]. Probabilistic Engineering Mechanics, 2001, 16(4)： 263-277.

[5]DU X, VENIGELLA P K, LIU D. Robust mechanism synthesis with random and interval variables[J]. Mechanism& Machine Theory, 2009, 7(7)： 1321-1337.

[6]BEN-HAIM Y. A non-probabilistic concept of reliability[J].Structural Safety, 1994, 14(4)： 227-245.

[7]ELISHAKOFF I, ELISSEEFF P, GLEGG S A L.Nonprobabilistic, convex-theoretic modeling of scatter in material properties[J]. AIAA Journal, 1994, 32(4)： 843-849.

[8]HUANG H Z, WANG Z L, LI Y F, et al. A nonprobabilistic set model of structural reliability based on satisfaction degree of interval[J]. Mechanika, 2011, 8(1)： 85-92.

[9]EARDLEY D. Quantification of margins and uncertainties(QMU)[R/OL].[2016-03-05]. https：//fas.org/irp/agency/dod/jason/margins.pdf.

[10]NAS, NRC. Evaluation of quantification of margins and uncertainties for assessing and certifying the reliability of the nuclear stockpile[R]. Washington： National Academy Press, 2008.

[11]PEPIN J E, RYTHERFORD A C, HEMEZ F M. Define a practical QMU metric[C]//49th AIAA Structures, Structural Dynamics and Materials Conference. [S.l.]： [s.n.], 2008.

[12]LUCAS L J, OWHADI H, ORTIZ M. Rigorous verification, validation, uncertainty quantification and certification through concentration-of-measure inequalities[J].Computer Methods in Applied Mechanics and Engineering,2008, 197(51-52)： 4591-4609.

[13]SENTZ K, FERSON S. Probabilistic bounding analysis in the quantification of margins and uncertainties[J]. Reliability Engineering and System Safety, 2011, 96： 1126-1136.

[14]URBINA A, MAHADEVAN S, PAEZ T L. Quantification of margins and uncertainties of complex systems in the presence of aleatory and epistemic uncertainty[J].Reliability Engineering and System Safety, 2011, 96：1114-1125.

[15]BAE H R, GRANDHI R V, CANFIELD R A. An approximation approach for uncertainty quantification using evidence theory[J]. Reliability Engineering and System Safety, 2004, 86(3)： 215-225.

[16]MARTIN J D, SIMPSON T W. Use of kriging models to approximate deterministic computer models[J]. AIAA Journal, 2005, 43(4)： 853-863.

[17]BANG I K, HAN D S, HAN G J, et al. Structural optimization for a jaw using iterative Kriging metamodels[J].Journal of Mechanical Science and Technology, 2008,22(9)： 1651-1659.

[18]WANG Z, WANG P. A maximum confidence enhancement based sequential sampling scheme for simulation-based design[J]. Journal of Mechanical Design, 2014, 136(2)：1-10.

[19]DU X P. Unified uncertainty analysis by the first order reliability method[J]. Journal of Mechanical Design, 2008,130(9)： 1-10.