孫 騫，馮曉波

(1. 武漢大學(xué) 水資源與水電工程科學(xué)國家重點實驗室，武漢 430072；2. 武漢大學(xué) 水工巖石力學(xué)教育部重點實驗室，武漢 430072)

0 引 言

很多工程可靠度分析涉及土體抗剪強(qiáng)度參數(shù)黏聚力c和內(nèi)摩擦角φ，確定土體抗剪強(qiáng)度參數(shù)的概率分布模型是土體工程可靠度分析的重要前提，其概率分布模型的準(zhǔn)確性直接影響著工程設(shè)計的安全性和經(jīng)濟(jì)性[1-3]。錯誤地使用某種分布當(dāng)做真實分布模型可能會導(dǎo)致較大的設(shè)計偏差，造成嚴(yán)重的后果[4]。

近年來，Copula 函數(shù)為土體抗剪強(qiáng)度參數(shù)模型表征提供了新的方法[4]。文獻(xiàn)中有許多表征相關(guān)結(jié)構(gòu)的Copula函數(shù)，例如Gaussian、t、Plackett、Frank、Clayton和Gumbel Copula[4]。每一種Copula函數(shù)都有獨立的結(jié)構(gòu)。Copula理論認(rèn)為任意聯(lián)合概率分布可以分解為參數(shù)邊緣分布函數(shù)和表示參數(shù)間相關(guān)結(jié)構(gòu)的Copula 函數(shù)。該方法實現(xiàn)了土體體參數(shù)間相關(guān)結(jié)構(gòu)的優(yōu)化，有效減小了傳統(tǒng)方法認(rèn)為相關(guān)結(jié)構(gòu)均服從Gaussian相關(guān)結(jié)構(gòu)而帶來的模型誤差[5-11]。

一般情況下常用AIC(Akaike Information Criterion)準(zhǔn)則來識別最優(yōu)邊緣分布和最優(yōu)Copula函數(shù)，該方法認(rèn)為在所有備選邊緣分布和Copula函數(shù)中，具有最小AIC值的邊緣分布和Copula函數(shù)即為最優(yōu)邊緣分布和最優(yōu)Copula函數(shù)。值得注意的是，無論是AIC值，還是用于估計邊緣分布參數(shù)和Copula函數(shù)參數(shù)的樣本均值、標(biāo)準(zhǔn)差、相關(guān)系數(shù)都源于有限的試驗數(shù)據(jù)，因此其準(zhǔn)確性和可靠性取決于試驗數(shù)據(jù)的樣本大小[12]。試驗表明，要想得到可靠的樣本均值、標(biāo)準(zhǔn)差、相關(guān)系數(shù)，樣本容量不得小于30[13]；要想得到可靠的邊緣分布和Copula函數(shù)識別結(jié)果，樣本容量需要達(dá)到100[14]。眾所周知，土體工程的樣本量屬于小樣本，且通常小于30，從小樣本中得出的二維分布模型顯然是不可靠的。因此，基于小樣本條件下的二維分布模型識別是一個具有挑戰(zhàn)性的土體工程難題。

此外，在有限的樣本外還可以得到一些先驗信息，包括文獻(xiàn)、工程師的判斷、當(dāng)?shù)氐墓こ探?jīng)驗等，這些先驗信息有助于得到最優(yōu)的二維分布模型。隨著貝葉斯理論的發(fā)展，越來越多的研究者證實了貝葉斯理論能夠充分地利用工程數(shù)據(jù)和先驗信息，將二者有機(jī)結(jié)合來描述土體工程參數(shù)模型的不確定性，在土體工程小樣本問題上是一種高效而實用的方法[21-25]。貝葉斯理論已被用于土體工程模型比較，例如利用貝葉斯理論識別土體中土層數(shù)目、土層厚度[15]，利用貝葉斯理論識別出了表征土體參數(shù)空間變異性的自相關(guān)函數(shù)[16]。

但上述貝葉斯理論僅限于單參數(shù)。為了表征兩個土體參數(shù)間的關(guān)系，Wang和Aladejare[17]采用貝葉斯方法導(dǎo)出單軸抗壓強(qiáng)度(UCS)和楊氏模量(E)的聯(lián)合概率分布。然而，他們直接使用二維正態(tài)分布模型來表征單軸抗壓強(qiáng)度(UCS)和楊氏模量(E)的分布，假定單軸抗壓強(qiáng)度(UCS)和楊氏模量(E)的分布均遵循單變量的一維正態(tài)分布，參數(shù)間相關(guān)性均遵循Gaussian Copula，完全沒有考慮二維分布模型選擇的不確定性。眾所周知，土體參數(shù)不一定遵循單變量正態(tài)分布，他們可能遵循單變量對數(shù)正態(tài)分布、Gumbel分布等等。同樣，土體參數(shù)間相關(guān)結(jié)構(gòu)也不一定遵循Gaussian Copula，他們可能遵循Plackett Copula、Frank Copula、No.16 Copula等等。更重要的是，土體工程二維分布模型中邊緣分布和Copula函數(shù)選擇對于土體工程可靠度有重要的影響。過去的研究[18]表明，選擇不同的邊緣分布和Copula函數(shù)所產(chǎn)生的失效概率的差異為幾個數(shù)量級。因此，在土體抗剪強(qiáng)度參數(shù)識別中，采用貝葉斯方法識別出最優(yōu)二維分布模型具有重要的意義。

本文的目的在于提出基于貝葉斯理論的土體抗剪強(qiáng)度參數(shù)最優(yōu)二維分布模型識別方法。簡要闡述了在擁有一定先驗信息和較少試驗數(shù)據(jù)的情形下，貝葉斯理論識別最優(yōu)二維分布模型的原理。采用蒙特卡洛模擬方法驗證了貝葉斯理論識別最優(yōu)二維分布模型的有效性。對比了貝葉斯獨立識別、貝葉斯非獨立識別、AIC一步識別，AIC兩步識別4種識別方法的識別能力，并分析了影響貝葉斯理論識別精度的主要因素。最后，搜集了29組實際工程土體抗剪強(qiáng)度參數(shù)試驗數(shù)據(jù)，研究了貝葉斯獨立識別和非獨立識別在實際工程土體抗剪強(qiáng)度參數(shù)最優(yōu)二維分布模型識別中的應(yīng)用。

1 表征土體抗剪強(qiáng)度參數(shù)相關(guān)結(jié)構(gòu)的Copula方法

根據(jù) Copula 理論，任意一個多元聯(lián)合分布都可以由相應(yīng)的邊緣分布和一個 Copula 函數(shù)組合而成， Copula 函數(shù)明確了變量間相關(guān)系數(shù)的大小和相關(guān)結(jié)構(gòu)的類型。對于土體抗剪強(qiáng)度參數(shù)c和φ來說，它們的聯(lián)合概率分布函數(shù)F(c,φ)和聯(lián)合概率密度函數(shù)f(c,φ)可以分別表示為[4]：

F(c,φ=C[F1(c),F2(φ);θ]=C(u,v;θ)

(1)

f(c,φ)=f1(c)f2(φ)D(F1,c),F2(φ);θ

(2)

式中：u=F1(c)、v=F2(φ)分別為c、φ的邊緣累積分布函數(shù)；f1(c)、f2(φ)分別為c、φ的概率密度函數(shù)；C(u,v;θ)為Copula函數(shù)；φ為Copula函數(shù)的相關(guān)參數(shù)；D[F1(c),F2(φ);θ]=D(u,v;θ)=?2C(u,v,θ)/?u?v為Copula密度函數(shù)。

Copula 函數(shù)構(gòu)造聯(lián)合分布模型的關(guān)鍵是確定Copula 函數(shù)相關(guān)參數(shù)和識別最優(yōu)Copula 函數(shù)。在確定Copula 函數(shù)相關(guān)參數(shù)時常采用秩相關(guān)系數(shù)法。Kendall 秩相關(guān)系數(shù)τ與Copula 函數(shù)C(u,v;θ)的關(guān)系式為[4]：

(3)

因此給定Kendall 秩相關(guān)系數(shù) 便可通過公式(3)求出參數(shù)θ。

在識別最優(yōu)Copula 函數(shù)方面，本文選取Gaussian、Plackett、Frank 和 No.16 Copula 函數(shù)來描述土體抗剪強(qiáng)度參數(shù)的統(tǒng)計負(fù)相關(guān)性[5-11]。它們的 Copula函數(shù)和參數(shù)取值范圍見表1。

2 土體抗剪強(qiáng)度參數(shù)最優(yōu)二維分布模型識別的AIC準(zhǔn)則

2.1 AIC一步識別

根據(jù)AIC準(zhǔn)則定義，具有最小AIC值的二維分布模型被認(rèn)為是最優(yōu)二維分布模型。AIC一步識別是指同時識別出土體抗剪強(qiáng)度參數(shù)二維分布模型中最優(yōu)邊緣分布和最優(yōu) Copula 函數(shù)。由于AIC 值定義為試驗數(shù)據(jù)點處二維分布模型密度函數(shù)值對數(shù)和的負(fù) 2倍與 2 倍二維分布模型參數(shù)數(shù)目之和，故AIC一步識別的計算表達(dá)式為[19]：

表1 4種Copula函數(shù)和參數(shù)取值范圍Tab.1 Summary of adopted 4 copula functions and

(4)

式中：Nd為試驗數(shù)據(jù)的樣本數(shù)目；m為二維分布模型中參數(shù)的數(shù)目，模型中共包含μc、σc、μφ、σφ、θ等參數(shù)，所以m=5；(ui,vi)為土體抗剪強(qiáng)度參數(shù)試驗數(shù)據(jù)(ci,φi)的經(jīng)驗分布值，將土體抗剪強(qiáng)度參數(shù)試驗數(shù)據(jù)D={(ci,φi),i=1,2,…,Nd} 轉(zhuǎn)化為標(biāo)準(zhǔn)均勻分布變量U={(ui,vi),i=1,2,…,Nd}的目的是為了方便計算，具體的計算公式為：

(5)

Rank為排序函數(shù)，排列順序為升序。

此外，AIC準(zhǔn)則中涉及的均值、標(biāo)準(zhǔn)差和相關(guān)系數(shù)是基于試驗數(shù)據(jù)采用極大似然估計方法得到的。

2.2 AIC兩步識別

AIC的兩步識別顧名思義是先識別出邊緣分布模型，再識別出Copula 函數(shù)，具有最小AIC值的邊緣分布模型和最小AIC值的Copula 函數(shù)組合后即為最優(yōu)的二維分布模型。識別邊緣分布模型的公式為[19]：

(6)

(7)

式中：Nd為試驗數(shù)據(jù)的樣本數(shù)目；m為邊緣分布模型中參數(shù)的數(shù)目，本文所選邊緣分布模型中的參數(shù)包括μ、σ，所以m=2。

識別Copula 函數(shù)的公式為：

(8)

式中：Nd為試驗數(shù)據(jù)的樣本數(shù)目；m為Copula 函數(shù)中參數(shù)的數(shù)目，由于本文 4 種備選 Copula 函數(shù)都是單參數(shù) Copula，因此m=1。

3 土體抗剪強(qiáng)度參數(shù)最優(yōu)二維分布模型識別的貝葉斯方法

基于Copula 理論構(gòu)造的二維分布模型可以拆分為兩個邊緣分布和一個Copula函數(shù)，通常根據(jù)土體抗剪強(qiáng)度參數(shù)二維分布模型中最優(yōu)邊緣分布和最優(yōu) Copula 函數(shù)是否被同時識別出將識別方法分為兩種：一步識別法和兩步識別法。貝葉斯一步識別需要五維積分，單次識別耗時長，計算效率低。兩步識別是指首先識別出表征參數(shù)c和φ的邊緣分布，再識別出表征參數(shù)間相關(guān)性的 Copula 函數(shù)。兩步識別法在識別過程中避免高維積分的運(yùn)算，大大提高了計算效率。對于工程師而言，計算效率的高低是決定是否采用該種識別方法的關(guān)鍵所在，能在可接受的精度范圍內(nèi)快速得出結(jié)論才是最重要的，所以本文只研究貝葉斯兩步識別。根據(jù)貝葉斯兩步識別中第二步識別過程與第一步識別過程是否有關(guān)，分為貝葉斯兩步獨立識別和貝葉斯兩步非獨立識別。顯然，獨立識別指前后兩步識別沒有關(guān)系，非獨立識別指第二次識別過程基于第一步的識別結(jié)果。為了簡化表達(dá)，文中用貝葉斯獨立識別和貝葉斯非獨立識別分別代替貝葉斯兩步獨立識別和貝葉斯兩步非獨立識別。

3.1 邊緣分布的識別

根據(jù)前人的研究成果，本文的備選邊緣分布類型采用正態(tài)分布、對數(shù)正態(tài)分布、極值Ⅰ型分布、伽馬分布，其對應(yīng)的分布函數(shù)和概率密度函數(shù)如表2所示。

表2 4種邊緣分布的分布函數(shù)和概率密度函數(shù)Tab.2 Summary of distribution and their probability density functions of 4 marginal distributions

下面以黏聚力c為例詳細(xì)闡述貝葉斯識別方法。

(9)

當(dāng)邊緣分布模型Mck中參數(shù)μc、σc未知時，可利用全概率公式將Pr(Dc|Mck)表示為：

(10)

式中：Pr(Dc|μc,σc,Mck)為給定Mck及其參數(shù)μc、σc條件下Dc發(fā)生的概率；Pr(μc,σc|Mck)為參數(shù)μc、σc的概率密度函數(shù)。假設(shè)Dc為黏聚力參數(shù)的Nd次獨立同分布觀測，那么Pr(Dc|μc,σc,Mck)可以表示為試驗數(shù)據(jù)點處備選邊緣概率函數(shù)的乘積：

(11)

當(dāng)先驗信息較少時，根據(jù)最大熵原理，通常假定μc、σc服從二維均勻分布，則概率密度函數(shù)Pr(μc,σc|Mck)的計算公式為：

(12)

3.2 Copula函數(shù)的識別

Copula 函數(shù)識別前需要確定的要素有：參數(shù)u和v數(shù)據(jù)，先驗信息以及備選 Copula 函數(shù)模型。參考前人研究成果，本文選取 Gaussian、Plackett、Frank 和No.16 Copula 函數(shù)為表征土體抗剪強(qiáng)度參數(shù)相關(guān)結(jié)構(gòu)的 Copula 函數(shù)，則C={Ck,k=1,2,3,4}={Gaussian, Plackett, Frank, No.16}。分別計算這4種備選 Copula函數(shù)Ck在給定土體抗剪強(qiáng)度參數(shù)試驗數(shù)據(jù)D={(ci,φi),i=1,2,…,Nd}條件下發(fā)生的概率Pr(Ck|D)。根據(jù)貝葉斯理論的基本原理，具有最大發(fā)生概率Pr(Ck|D)的Copula 函數(shù)即為最優(yōu)Copula 函數(shù)。

(13)

當(dāng)Copula 函數(shù)Ck中參數(shù)θ未知時，可利用全概率公式將Pr(D|Ck)表示為：

(14)

(15)

式中：Pr(D|τ,Ck)為給定Ck及τ在的條件下D出現(xiàn)的概率；Pr(τ|Ck)為τ的概率密度函數(shù)。同理，可假設(shè)τ在區(qū)間[τmin,τmax]內(nèi)服從均勻分布，則Pr(τ|Ck)=1/(τmax-τmin)。假設(shè)D為獨立同分布觀測，那么Pr(D|τ,Ck)可以表示為試驗數(shù)據(jù)點處備選Copula 密度函數(shù)的乘積：

(16)

顯然，在貝葉斯理論框架下識別Copula函數(shù)的首要任務(wù)是確定變量U={(ui,vi),i=1,2,…,Nd}。由Copula函數(shù)的定義可知， Copula 函數(shù)的邊緣分布在區(qū)間[0,1]上服從均勻分布。根據(jù)變量U確定方法的不同，貝葉斯方法分為獨立識別和非獨立識別。獨立識別指兩步識別沒有依賴關(guān)系，非獨立識別指第二步識別以第一步識別結(jié)果為基礎(chǔ)。

3.2.1 獨立識別

由于第一步識別和第二步識別沒有關(guān)系，所以變量U是基于原始數(shù)據(jù)的經(jīng)驗分布值確定的。由公式(5)將數(shù)據(jù)D={(ci,φi),i=1,2,…,Nd}轉(zhuǎn)化為標(biāo)準(zhǔn)均勻分布變量U={(ui,vi),i=1,2,…,Nd}。

3.2.2 非獨立識別

根據(jù)Copula理論可知，在變量U中u和v可分別視為參數(shù)c和φ的累積分布函數(shù)，即u=F1(c|μc,σc)和v=F2(φ|μφ,σφ)。第一步邊緣分布模型識別后可以確定參數(shù)c和φ的累積分布函數(shù)類型，在確定μc和σc以及μφ和σφ后便可以確定c和φ的累積分布函數(shù)F1(c|μc,σc)和F2(φ|μφ,σφ)。下面以參數(shù)c為例詳述如何在已知累積分布函數(shù)類型后得到μc和σc的最可能值，下文簡稱MPV值。依據(jù)貝葉斯原理可得出μc和σc的后驗分布為：

(17)

4 最優(yōu)二維分布模型識別方法的對比

本節(jié)采用蒙特卡洛模擬方法驗證貝葉斯理論識別最優(yōu)二維分布模型的有效性。為了驗證該方法的有效性，將貝葉斯識別結(jié)果與AIC 準(zhǔn)則識別結(jié)果進(jìn)行對比。試驗表明1 000 次模擬能夠得到非常穩(wěn)健的最優(yōu)二維分布模型識別結(jié)果，所以本節(jié)試驗方案設(shè)計為：在給定真實二維分布模型情況下，重復(fù)模擬1 000組樣本數(shù)目為Nd的服從給定真實二維分布模型的樣本，分別采用貝葉斯獨立識別、貝葉斯非獨立識別、AIC一步識別、AIC兩步識別共計4種方法給出每組樣本的最優(yōu)二維分布模型，最后統(tǒng)計并對比1 000次模擬樣本中真實二維分布模型被識別為最優(yōu)二維分布模型的次數(shù)。本文擬采用正態(tài)(Normal)分布、對數(shù)正態(tài)分布(Lognormal)、極值Ⅰ型(Gumbel)分布、伽馬(Gamma)分布作為模擬算例的邊緣分布類型，選取Gaussian、Plackett、Frank 和No.16 Copula 函數(shù)為表征模擬算例土體抗剪強(qiáng)度參數(shù)相關(guān)結(jié)構(gòu)的 Copula 函數(shù)，分析這4種方法的識別能力與不同真實二維分布模型、樣本數(shù)目以及參數(shù)相關(guān)性的關(guān)系。對于真實二維分布模型，假定μc=66，σc=22，μφ=28，σφ=3.5。選取樣本數(shù)目Nd分別等于 30、50和100以及相關(guān)系數(shù)τ分別等于-0.25、-0.50 和-0.75 ，因此c和φ各4種備選邊緣分布模型、4種Copula 函數(shù)構(gòu)成了4×4×4×3×3=576種備選二維分布模型，再加上3種樣本數(shù)目和3種相關(guān)系數(shù)，共組成 組模擬方案。表3 給出了 64 種備選二維分布模型。

表3 64種備選二維分布模型Tab.3 64 kinds of alternative bivariate distribution models

注：在c和φ中， “1”，“2”，“3”和“4”分別代表“正態(tài)(Normal)分布”，“對數(shù)正態(tài)(Lognormal)分布”，“極值Ⅰ型(Gumbel)分布”和“伽馬(Gamma)分布”；在Copula中，“1”，“2”，“3”和“4”分別代表“Gaussian”，“Plackett”，“Frank”和“No.16”。

這里假設(shè)64種備選二維分布模型的先驗概率相等，由于本文采用兩步識別，每步識別中的備選模型數(shù)量均為4，因此先驗概率皆為1/4。在每組方案模擬的1000組樣本中，當(dāng)真實的二維分布模型被識別為最優(yōu)二維分布模型的次數(shù)大于其余63種模型時，認(rèn)為識別成功；否則認(rèn)為識別失敗。表4為上述模擬試驗中4種方法的識別結(jié)果。

如表4所示，經(jīng)對比發(fā)現(xiàn)，相同條件下兩種貝葉斯方法成功識別的概率均高于兩種AIC方法成功識別的概率。

就單個模型在1 000次模擬試驗的識別精度而言，貝葉斯方法也優(yōu)于AIC準(zhǔn)則。由于模型種類多，且4種識別方法的識別精度規(guī)律保持相對穩(wěn)定，所以只選取了部分試驗數(shù)據(jù)作圖1來展示4種方法識別精度的高低。

表4 576組模擬試驗中4種方法成功識別真實二維分布模型的組數(shù)Tab.4 Number of successful identifications over 576 types using Bayesian method

圖1 4種識別方法成功識別次數(shù)比較Fig.1 Number of successful identifications of 4 methods

顯然，兩種貝葉斯識別方法的識別精度普遍優(yōu)于兩種AIC識別方法的識別精度，且貝葉斯非獨立識別在多數(shù)情形下優(yōu)于貝葉斯獨立識別，AIC一步識別和兩步識別在識別精度上并無明顯差異。造成貝葉斯獨立識別和非獨立識別精度差異的主要原因是Copula函數(shù)的變量不同。以小浪底大壩斜心墻三軸固結(jié)排水(CD)土體抗剪強(qiáng)度參數(shù)數(shù)據(jù)為例，通過兩種方法得到變量U，變量U的分布情況如圖2所示。由圖2可知，兩種方法使用的變量U在整體趨勢上一致，但是并不完全相同， Copula函數(shù)的變量不同導(dǎo)致了識別結(jié)果的差異。由于這4種方法在計算機(jī)上運(yùn)行一次的速度均在3秒以內(nèi)，所以在識別土體抗剪強(qiáng)度參數(shù)最優(yōu)二維分布模型時，建議優(yōu)先采用貝葉斯非獨立法進(jìn)行識別。

圖2 貝葉斯獨立法和非獨立法中Copula函數(shù)變量散點圖Fig.2 Scatter plots of Copula function in Bayesian theory bivariate distribution identification

5 影響識別精度的因素分析

前文驗證了貝葉斯非獨立識別方法對最優(yōu)二維分布模型識別的高效性，下面進(jìn)一步分析影響該方法識別精度的主要因素。比較后發(fā)現(xiàn)，影響最優(yōu)二維分布模型識別精度的主要因素是樣本數(shù)目、參數(shù)間相關(guān)性、備選二維分布模型集合。當(dāng)然，備選二維分布模型的先驗概率對真實二維分布模型的識別結(jié)果也有明顯的影響。由于篇幅有限，圖3僅展示了4種真實二維分布模型的識別結(jié)果，真實二維分布模型分別為[1,1,1]、[1,1,2]、[1,1,3]、[1,1,4]。

5.1 樣本數(shù)目

從圖3(a)～3(d)中可以看出，真實二維分布模型被識別為最優(yōu)二維分布模型的次數(shù)隨樣本數(shù)目的增加而增大。此外，當(dāng)樣本數(shù)目逐漸增大時，單位樣本數(shù)目的增加對于識別精度的提高效果逐漸減小，這說明識別精度在樣本數(shù)目較少時對樣本數(shù)目的變化更敏感。因此若能在實際工程中獲得更多的土體抗剪強(qiáng)度參數(shù)試驗數(shù)據(jù)，則可以得到更可靠的二維分布模型，提高工程設(shè)計的可靠度。

5.2 參數(shù)間相關(guān)性

從圖3(a)～3(d)中可以看出，參數(shù)間相關(guān)性對最優(yōu)二維分布模型的識別結(jié)果具有重要的影響。參數(shù)間相關(guān)性主要是通過影響Copula函數(shù)的識別從而影響二維分布模型的識別。通過對比發(fā)現(xiàn)，當(dāng)真實Copula函數(shù)為Gaussian、Plackett 和 Frank Copula 函數(shù)時，真實二維分布模型被識別為最優(yōu)二維分布模型的次數(shù)隨參數(shù)間負(fù)相關(guān)系數(shù)的增加而增大。而當(dāng)真實Copula函數(shù)為No.16 Copula 函數(shù)時，真實二維分布模型被識別為最優(yōu)二維分布模型的次數(shù)隨參數(shù)間負(fù)相關(guān)系數(shù)的增加而減小。這是因為大多數(shù)Copula函數(shù)在相關(guān)系數(shù)趨近于0時都收斂于獨立 Copula 函數(shù)，此時不同的Copula函數(shù)差異性很小，可以忽略，當(dāng)相關(guān)系數(shù)的絕對值逐漸增大時，不同的Copula函數(shù)差異性也隨之增大，因此更容易識別，本文中Gaussian、Plackett 和 Frank Copula 函數(shù)便屬于這一類；當(dāng)然，還存在少數(shù)的Copula函數(shù)，它們在相關(guān)系數(shù)趨近于0時不收斂于獨立 Copula 函數(shù)，且當(dāng)相關(guān)系數(shù)的絕對值逐漸增大時，它們與其他Copula 函數(shù)的差異性逐漸減小，因此識別難度增加，本文中No.16 Copula 函數(shù)便屬于這一類。因此，參數(shù)間相關(guān)性大小對于最優(yōu)二維分布模型的識別精度沒有完全統(tǒng)一的規(guī)律，換言之，大部分二維分布模型的識別精度隨相關(guān)性的增加而增大，少數(shù)二維分布模型的識別精度隨相關(guān)性的增加而降低。

圖3 4種真實二維分布模型識別能力對比圖Fig.3 Identification ability of 4 real bivariate distribution models

5.3 備選二維分布模型集合

前文分析了參數(shù)間相關(guān)性對二維分布模型識別精度的影響。由于真實二維分布模型只選取了一種邊緣分布類型，故只分析備選二維分布模型中Copula 函數(shù)的差異對識別能力的影響。對比圖3(a)～3(d)不難發(fā)現(xiàn)，圖3(a)和圖3(d)的成功識別次數(shù)明顯高于圖3(b)和圖3(c)，這是因為一般來說， Plackett 和 Frank Copula 函數(shù)在相關(guān)系數(shù)相同時具有相似的相關(guān)結(jié)構(gòu)，而 Gaussian 和 No.16 Copula 函數(shù)相關(guān)結(jié)構(gòu)與其余 3 種備選 Copula 函數(shù)存在較大差別。因此，當(dāng)真實Copula函數(shù)為Gaussian或No.16 Copula 函數(shù)時，真實二維分布模型被識別為最優(yōu)二維分布模型的次數(shù)遠(yuǎn)遠(yuǎn)大于真實Copula函數(shù)為 Plackett 或Frank Copula 函數(shù)的二維分布模型。因此，真實 Copula 函數(shù)與其余備選 Copula 函數(shù)之間存在的差異越大就越容易被識別成功，需要的樣本數(shù)目越少。

5.4 先驗信息

由于信息缺乏，前文的分析假定64 種備選二維分布模型的先驗概率都相等。當(dāng)有足夠的證據(jù)表明某種二維分布模型作為最優(yōu)二維分布模型的概率明顯大于其他二維分布模型作為最優(yōu)二維分布模型的概率時，應(yīng)代入不等的先驗概率進(jìn)行計算。就本文的模擬試驗而言，由于樣本均生成于真實的二維分布模型，故樣本來自于真實二維分布模型的概率明顯大于其他備選二維分布模型。

假定真實二維分布模型的概率為2/5，剩余63種備選模型均為1/5，相關(guān)系數(shù)τ=-0.5。計算結(jié)果如圖4所示：

圖4 先驗信息對識別精度的影響Fig.4 The prior information influences on identification

從圖4中可知，考慮不相等的先驗概率的識別結(jié)果明顯優(yōu)于相等的先驗概率的識別結(jié)果，真實二維分布模型被識別為最優(yōu)二維分布模型的概率顯著提高。

6 工程應(yīng)用

本節(jié)探討貝葉斯理論在實際工程土體抗剪強(qiáng)度參數(shù)最優(yōu)二維分布模型識別中的應(yīng)用。本文共搜集了全世界范圍內(nèi) 29 組土體抗剪強(qiáng)度參數(shù)的現(xiàn)場或者室內(nèi)試驗數(shù)據(jù)，大部分的試驗數(shù)據(jù)的樣本數(shù)目小于30，屬于小樣本。本文提出的基于貝葉斯理論的土體抗剪強(qiáng)度參數(shù)最優(yōu)二維分布模型識別方法相較AIC準(zhǔn)則的識別方法能在小樣本條件下更有效地識別出最優(yōu)二維分布模型。本節(jié)選擇小浪底大壩斜心墻三軸固結(jié)排水(CD)土體抗剪強(qiáng)度參數(shù)數(shù)據(jù)為例，闡述基于貝葉斯理論識別最優(yōu)二維分布模型的計算步驟：

(1)獲取土體抗剪強(qiáng)度參數(shù)試驗數(shù)據(jù)D={(ci,φi),i=1,2,…,63} 轉(zhuǎn)化為標(biāo)準(zhǔn)均勻分布變量U={(ui,vi),i=1,2,…,n}。

(2)土體抗剪強(qiáng)度參數(shù)存在較強(qiáng)的統(tǒng)計負(fù)相關(guān)性，Kendall 秩相關(guān)系數(shù)τ=-0.38。選取本文中的64種二維分布模型為備選模型，設(shè)定τ的積分區(qū)間為[-1,0]，此時τmin=-1，τmax=0。

(3)采用公式(9)計算出64種備選模型在給定試驗數(shù)據(jù)D條件下的后驗概率Pr(Mk|D)， 具有最大后驗概率的備選模型即為最優(yōu)二維分布模型。

對于黏聚力c而言，根據(jù)公式(9)分別計算正態(tài)分布、對數(shù)正態(tài)分布、極值Ⅰ型分布、伽馬分布發(fā)生的概率分別為： 96.88%、 0.12%、 2.54%、 0.46%，因而表征本例中黏聚力邊緣分布的最優(yōu)模型為正態(tài)分布；對于內(nèi)摩擦角φ而言，正態(tài)分布、對數(shù)正態(tài)分布、極值Ⅰ型分布、伽馬分布發(fā)生的概率分別為：30.69%、 18.42%、 2.08%、48.81%，因而表征本例中內(nèi)摩擦角邊緣分布的最優(yōu)模型為伽馬分布；同理，對于Copula函數(shù)而言，Gaussian、Plackett、Frank 和No.16發(fā)生的概率分別為：24.81%、 32.26%、 42.66%、 0.27%，因而表征本例中Copula函數(shù)的最優(yōu)模型為Frank Copula函數(shù)。可見，黏聚力邊緣分布為正態(tài)分布，內(nèi)摩擦角邊緣分布為伽馬分布，相關(guān)結(jié)構(gòu)為Frank Copula函數(shù)的二維分布模型為64種備選模型中能最優(yōu)地表征小浪底樞紐工程固結(jié)排水(CD)土體抗剪強(qiáng)度參數(shù)試驗數(shù)據(jù)二維分布的模型。表5給出了模擬樣本c和φ的均值和標(biāo)準(zhǔn)差，等效樣本的統(tǒng)計量均值和標(biāo)準(zhǔn)差與原始樣本數(shù)據(jù)相比，其相對誤差都在1%以內(nèi)，模擬數(shù)據(jù)的Kendall 秩相關(guān)系數(shù) ，與樣本數(shù)據(jù)相比相對誤差僅為2.63%，這些微小的相對誤差說明蒙特卡洛模擬樣本能夠比較準(zhǔn)確地還原原始樣本的數(shù)據(jù)特征。圖5給出了模擬數(shù)據(jù)和樣本數(shù)據(jù)的分布情況，從圖5中可以看出，模擬數(shù)據(jù)基本覆蓋了樣本數(shù)據(jù)的分布區(qū)域，兩者擬合度良好，能夠有效地反映試驗數(shù)據(jù)的二維分布情況。

表5 土體抗剪強(qiáng)度參數(shù)模擬樣本的統(tǒng)計特征值Tab.5 Statistical feature values of shear strength parameter simulation samples

圖5 土體抗剪強(qiáng)度參數(shù)等效樣本散點圖Fig.5 Scatter plots of shear strength parameter simulation samples

重復(fù)以上3步，即可得出29組數(shù)據(jù)對應(yīng)的最優(yōu)二維分布模型，計算結(jié)果見表6，該表統(tǒng)計了貝葉斯獨立和非獨立兩種方法的識別結(jié)果。經(jīng)統(tǒng)計，黏聚力的最優(yōu)分布模型依次為：正態(tài)分布，伽馬分布，對數(shù)正態(tài)分布和極值Ⅰ型分布；對于內(nèi)摩擦角的最優(yōu)分布模型依次為：正態(tài)分布，伽馬分布，極值Ⅰ型分布和對數(shù)正態(tài)分布；表征參數(shù)相關(guān)結(jié)構(gòu)的最優(yōu)Copula 函數(shù)依次為： Gaussian， Frank， Plackett， No.16?？梢娫趯嶋H工程中土體抗剪強(qiáng)度參數(shù)存在各種類型的分布情況，要具體工程具體分析，不能一概而論，否則在計算土體結(jié)構(gòu)物可靠度時會低估或高估實際失效概率，造成工程失事風(fēng)險或工程資源浪費，這進(jìn)一步證明了考慮土體抗剪強(qiáng)度參數(shù)二維分布模型的不確定性是非常有必要的。

表6 29組工程數(shù)據(jù)貝葉斯獨立識別和非獨立識別的結(jié)果統(tǒng)計Tab.6 Summary of 29 sets of shear strength parameters and the results of their best-fit models

7 結(jié) 論

(1)貝葉斯理論能夠有效地識別表征土體抗剪強(qiáng)度參數(shù)最優(yōu)的二維分布模型。該方法識別二維分布模型時，不需要估計二維分布模型的參數(shù)，并能與工程上現(xiàn)有的工程經(jīng)驗等先驗信息相結(jié)合，為降低模型選擇的不確定性提供了一條高效且穩(wěn)定的途徑。

(2)與常用的基于AIC準(zhǔn)則的識別方法相比，貝葉斯理論在識別能力和識別精度上表現(xiàn)更出色，特別是在樣本數(shù)目較小的情況下，貝葉斯理論的優(yōu)勢表現(xiàn)得更為明顯。其中貝葉斯非獨立識別方法可以認(rèn)為是兼顧識別能力和識別精度的最優(yōu)選擇。

(3)土體抗剪強(qiáng)度參數(shù)的樣本數(shù)目、參數(shù)間相關(guān)性、備選二維分布模型集合以及先驗信息都顯著影響貝葉斯理論的識別精度。通常，樣本數(shù)目越大、參數(shù)間相關(guān)性越強(qiáng)、備選二維分布模型集合中真實 Copula 函數(shù)與其余備選Copula 函數(shù)差異越大、真實二維分布模型具有越高的先驗概率，貝葉斯理論的識別精度就會越高。

□