基于Bootstrap方法的巖土體參數(shù)聯(lián)合分布模型識別

唐小松 ，李典慶 ，周創(chuàng)兵 ，方國光

（1.武漢大學(xué) 水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072；2.武漢大學(xué) 水工巖石力學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072）

1 引 言

巖土工程可靠度分析中經(jīng)常包含相關(guān)非正態(tài)變量，如巖土體的抗剪強(qiáng)度參數(shù)黏聚力和內(nèi)摩擦角間存在統(tǒng)計(jì)負(fù)相關(guān)性[1－2]、基樁荷載-位移雙曲線參數(shù)具有負(fù)相關(guān)性[3－5]，非飽和土的土-水特征曲線參數(shù)亦具有明顯的負(fù)相關(guān)關(guān)系[6]，而且這些參數(shù)一般都服從非正態(tài)分布。眾所周知，巖土工程可靠度分析通常需要已知相關(guān)非正態(tài)巖土體參數(shù)的聯(lián)合概率分布函數(shù)。近年來，Copula理論[7]為相關(guān)非正態(tài)巖土體參數(shù)聯(lián)合概率分布函數(shù)的建立提供了一種新的途徑，它在巖土工程領(lǐng)域中逐漸得到應(yīng)用[8－10]。究其原因，在Copula理論框架下構(gòu)造巖土體參數(shù)聯(lián)合概率分布函數(shù)具有以下2方面不可替代的優(yōu)點(diǎn)[7]：其一是該方法將巖土體參數(shù)聯(lián)合概率分布函數(shù)的構(gòu)造簡化為巖土體參數(shù)邊緣分布函數(shù)估計(jì)和Copula函數(shù)選擇問題，且邊緣分布函數(shù)估計(jì)和Copula函數(shù)選擇分開獨(dú)立進(jìn)行；其二是該方法能夠?qū)⒉煌愋偷倪吘壏植己瘮?shù)和Copula函數(shù)結(jié)合在一起，從而構(gòu)造出具有任意邊緣分布函數(shù)及任意相關(guān)結(jié)構(gòu)巖土體參數(shù)的聯(lián)合概率分布函數(shù)。鑒于Copula理論的上述優(yōu)點(diǎn)，本文主要研究基于Copula理論的相關(guān)非正態(tài)巖土體參數(shù)聯(lián)合概率分布函數(shù)構(gòu)造方法。

基于Copula理論的巖土體參數(shù)聯(lián)合概率分布函數(shù)構(gòu)造方法中，巖土體參數(shù)的最優(yōu)邊緣分布函數(shù)和最優(yōu)Copula函數(shù)都是基于實(shí)測數(shù)據(jù)采用AIC（Akaike information criterion）準(zhǔn)則[11]優(yōu)化比較選擇確定的，即具有最小AIC 值的邊緣分布函數(shù)和Copula函數(shù)被認(rèn)為是擬合原始觀測數(shù)據(jù)概率分布特性和相關(guān)結(jié)構(gòu)最優(yōu)的概率分布模型。AIC 準(zhǔn)則中的AIC 值是巖土體參數(shù)實(shí)測數(shù)據(jù)的函數(shù)，因此，是樣本的函數(shù)。而樣本是一個(gè)N 維隨機(jī)向量，N為樣本容量，從而樣本的函數(shù)亦是一個(gè)隨機(jī)變量，這個(gè)隨機(jī)變量稱為統(tǒng)計(jì)量。統(tǒng)計(jì)量的變異性大小由樣本容量決定，樣本容量越小，變異性越大。一般來說，樣本數(shù)目小于100時(shí)屬于小樣本容量。巖土體參數(shù)具有小樣本容量是巖土工程中一個(gè)客觀存在的事實(shí)[1－6]，基于有限巖土體參數(shù)試驗(yàn)數(shù)據(jù)計(jì)算得到的AIC 值通常具有較大的變異性。傳統(tǒng)方法沒有考慮AIC 值的變異性，而僅僅基于巖土體參數(shù)實(shí)測數(shù)據(jù)求得不同備選邊緣分布函數(shù)和Copula函數(shù)的單一AIC值識別最優(yōu)邊緣分布函數(shù)和最優(yōu)Copula函數(shù)。眾所周知，隨機(jī)變量的單一取值并不能代表這個(gè)變量的實(shí)際分布，因而基于單一AIC 值的最優(yōu)邊緣分布函數(shù)和最優(yōu)Copula函數(shù)識別方法是不合理的。為了得到合理的識別結(jié)果，應(yīng)該充分考慮小樣本容量巖土體參數(shù)引起的統(tǒng)計(jì)量AIC 值的變異性，其關(guān)鍵難點(diǎn)主要體現(xiàn)在以下兩方面：①已知巖土體參數(shù)實(shí)測數(shù)據(jù)如何模擬AIC 值的變異性；②考慮AIC 值變異性時(shí)如何表述識別結(jié)果。對于上述問題，迄今為止還未見系統(tǒng)的研究。

近年來在統(tǒng)計(jì)學(xué)領(lǐng)域應(yīng)用較多的Bootstrap方法[12]為小樣本條件下統(tǒng)計(jì)量分布的直接模擬提供了一種新的工具。Bootstrap方法最早由Efron[12]于1979年提出，近年來在巖土工程領(lǐng)域中得到了應(yīng)用[13－14]；如謝桂華等[13]采用Bootstrap方法模擬了地基沉降計(jì)算經(jīng)驗(yàn)系數(shù)均值的分布及置信區(qū)間；Luo等[14]利用Bootstrap方法研究了小樣本容量引起的土性參數(shù)統(tǒng)計(jì)量變異性對支護(hù)開挖失效概率的影響。Bootstrap方法的基本思想是通過對原始觀測樣本數(shù)據(jù)進(jìn)行有放回地隨機(jī)抽樣得到大量與原始觀測樣本相同樣本容量的Bootstrap子樣本，然后基于Bootstrap子樣本計(jì)算統(tǒng)計(jì)量的估計(jì)值，最終獲得統(tǒng)計(jì)量的變異系數(shù)及其概率分布。該方法只依賴于給定的原始觀測樣本數(shù)據(jù)，且能夠充分挖掘原始觀測數(shù)據(jù)所攜帶的總體信息，不需要對巖土體參數(shù)的實(shí)際分布作任何假設(shè)以及增加新的數(shù)據(jù)觀測，因此，該方法是一種非參數(shù)統(tǒng)計(jì)方法。從本質(zhì)上說，Bootstrap方法是一種將小樣本問題轉(zhuǎn)化為大樣本問題的分析方法。鑒于Bootstrap方法的上述優(yōu)勢，本文將Bootstrap方法應(yīng)用于小樣本容量巖土體參數(shù)AIC 值分布的模擬問題，提出了基于Bootstrap方法的巖土體參數(shù)最優(yōu)邊緣分布函數(shù)和最優(yōu)Copula函數(shù)識別方法。最后以基樁荷載-位移雙曲線參數(shù)試驗(yàn)數(shù)據(jù)為例證明了所提方法的有效性。

2 巖土體參數(shù)聯(lián)合分布函數(shù)構(gòu)造的Copula函數(shù)方法

本文采用Copula函數(shù)構(gòu)造巖土體參數(shù)的聯(lián)合概率分布函數(shù)，下面簡要介紹Copula函數(shù)基本定義以及基于Copula函數(shù)的巖土體參數(shù)聯(lián)合概率分布函數(shù)構(gòu)造方法。Copula函數(shù)是將變量的聯(lián)合分布與其邊緣分布聯(lián)結(jié)起來的函數(shù)，它惟一地描述了變量間相關(guān)性信息，包括變量間相關(guān)系數(shù)大小以及變量間相關(guān)結(jié)構(gòu)類型。對于二維情況，Copula函數(shù)定義為[0，1]2空間中邊緣分布為[0，1]區(qū)間內(nèi)均勻分布的二維聯(lián)合概率分布函數(shù)[7]：

式中：F(x1，x2)為變量X1和X2的聯(lián)合概率分布函數(shù)；u1=F1(x1)和u2=F2(x2)分別為變量X1和X2的邊緣分布函數(shù)；C為Copula函數(shù)；θ為Copula函數(shù)的相關(guān)參數(shù)。相應(yīng)地，若已知變量X1和X2的概率密度函數(shù)f1(x1)和f2(x2)，則通過對式（1）兩邊求導(dǎo)可得變量X1和X2的聯(lián)合概率密度函數(shù)f (x1，x2)為

式中：c（F1（x1），F(xiàn)2（x2）;θ）=?2C（u1，u2;θ）/?u1?u2，為Copula密度函數(shù)。因此，若已知變量X1和X2的邊緣分布函數(shù)和Copula函數(shù)，利用式（1）和式（2）就可以建立變量X1和X2的二維聯(lián)合概率分布模型。Copula函數(shù)的相關(guān)參數(shù)θ 可由變量X1和X2間的Kendall秩相關(guān)系數(shù)τ 由下式直接求出[7]：

因此，當(dāng)變量X1和X2間的Kendall秩相關(guān)系數(shù)τ 已知時(shí)，通過求解式（3）所示的非線性積分方程就可得出參數(shù)θ 。從式（1）和式（2）可以看出，相關(guān)非正態(tài)巖土體參數(shù)聯(lián)合概率分布模型的建立包括2個(gè)步驟：第1步是建立巖土體參數(shù)的邊緣分布函數(shù)；第2步是選擇最優(yōu)的Copula函數(shù)描述巖土體參數(shù)間的相關(guān)結(jié)構(gòu)。如前所述，上述2步可以分開獨(dú)立進(jìn)行，下面分別給出巖土體參數(shù)最優(yōu)邊緣分布函數(shù)和最優(yōu)Copula函數(shù)識別方法。

本文采用工程上常用的AIC準(zhǔn)則[11]識別最優(yōu)的邊緣分布函數(shù)，即具有最小AIC 值的邊緣分布函數(shù)被認(rèn)為是擬合原始觀測數(shù)據(jù)概率分布特性最優(yōu)的邊緣分布函數(shù)。AIC 值定義為變量原始觀測數(shù)據(jù)點(diǎn)處概率密度函數(shù)值對數(shù)和的-2倍與2倍概率密度函數(shù)分布參數(shù)數(shù)目之和，表達(dá)式為

式中：ε為邊緣分布函數(shù)的分布參數(shù)向量；k1為邊緣分布函數(shù)中分布參數(shù)的數(shù)目；N為原始觀測數(shù)據(jù)的樣本容量。因此，已知巖土體參數(shù)的一組觀測數(shù)據(jù)，就可以采用式（4）求出多種備選邊緣分布函數(shù)的AIC 值，從而識別出擬合原始觀測數(shù)據(jù)概率分布特性最優(yōu)的邊緣分布函數(shù)。本文選取巖土工程中常用的正態(tài)分布、對數(shù)正態(tài)分布、極值I型分布和威布爾分布作為備選邊緣分布函數(shù)擬合巖土體參數(shù)的概率分布特性[15－18]。為了避免出現(xiàn)負(fù)值，正態(tài)分布和極值I型分布在0處進(jìn)行左截尾。上述4種邊緣分布函數(shù)都包含2個(gè)分布參數(shù)，因此k1=2，這些分布參數(shù)ε 可以利用巖土體參數(shù)觀測數(shù)據(jù)的均值和標(biāo)準(zhǔn)差求出。

最優(yōu)的Copula函數(shù)也可以采用AIC準(zhǔn)則[11]識別，即具有最小AIC 值的Copula函數(shù)被認(rèn)為是擬合原始觀測數(shù)據(jù)相關(guān)結(jié)構(gòu)最優(yōu)的Copula函數(shù)。AIC 值定義為變量原始觀測數(shù)據(jù)點(diǎn)處Copula密度函數(shù)值對數(shù)和的-2倍與2倍Copula函數(shù)相關(guān)參數(shù)數(shù)目之和，表達(dá)式為

式中：k2為Copula函數(shù)中相關(guān)參數(shù)的數(shù)目，對于二維Copula函數(shù)來說k2=1；(u1i，u2i)，為原始觀測數(shù)據(jù)(x1i，x2i)的經(jīng)驗(yàn)分布值，可由下式計(jì)算：

式中：rank(x1i)和rank(x2i)分別為按升序排列時(shí)實(shí)測值x1i和x2i在整列觀測數(shù)據(jù)x1={x11，x12，···，x1N}和x2={x21，x22，···，x2N}中的秩次。因此，已知巖土體參數(shù)的一組觀測數(shù)據(jù)，就可以采用式（5）求出多種備選Copula函數(shù)的AIC 值，從而識別出擬合原始觀測數(shù)據(jù)相關(guān)結(jié)構(gòu)最優(yōu)的Copula函數(shù)。本文選取巖土工程中常用的Gaussian、Plackett、Frank和No.16 Copula函數(shù)作為備選Copula函數(shù)擬合巖土體參數(shù)間的相關(guān)結(jié)構(gòu)[8－10]。上述4種Copula函數(shù)都能描述變量間的負(fù)相關(guān)性，且相關(guān)系數(shù)都能達(dá)到-1，它們非常適合描述負(fù)相關(guān)性較強(qiáng)的巖土體參數(shù)間相關(guān)關(guān)系[8－10]。這些Copula函數(shù)的相關(guān)參數(shù)θ 可以利用巖土體參數(shù)觀測數(shù)據(jù)的Kendall秩相關(guān)系數(shù)由式（3）求出。

從式（4）和式（5）可以看出各備選邊緣分布函數(shù)和Copula函數(shù)的AIC值都是巖土體參數(shù)試驗(yàn)數(shù)據(jù)的函數(shù)，因此，都是樣本的函數(shù)。由統(tǒng)計(jì)學(xué)理論可知，樣本是一個(gè)N 維隨機(jī)向量，從而AIC 值這個(gè)樣本的函數(shù)亦是一個(gè)隨機(jī)變量，它稱為統(tǒng)計(jì)量。統(tǒng)計(jì)量的變異性大小由樣本容量決定，樣本容量越小變異性越大。實(shí)際工程中巖土體參數(shù)試驗(yàn)數(shù)據(jù)通常具有較小的樣本容量，基于這些有限試驗(yàn)數(shù)據(jù)計(jì)算得到的AIC 值具有較大的變異性。傳統(tǒng)方法沒有考慮AIC 值變異性，而僅僅基于巖土體參數(shù)實(shí)測數(shù)據(jù)求得不同備選邊緣分布函數(shù)和Copula函數(shù)的單一AIC 值識別最優(yōu)邊緣分布函數(shù)和Copula函數(shù)。眾所周知，隨機(jī)變量的單一取值并不能代表這個(gè)變量的實(shí)際分布，因而基于單一AIC 值的最優(yōu)邊緣分布函數(shù)和Copula函數(shù)識別方法是不合理的。為了得到合理結(jié)果，應(yīng)該充分考慮小樣本容量巖土體參數(shù)引起的統(tǒng)計(jì)量AIC 值的變異性。下面介紹統(tǒng)計(jì)量AIC 值變異性模擬的Bootstrap方法以及基于Bootstrap的最優(yōu)邊緣分布函數(shù)和Copula函數(shù)識別方法。

3 巖土體參數(shù)聯(lián)合分布函數(shù)識別的Bootstrap方法

本文采用Bootstrap方法模擬統(tǒng)計(jì)量AIC值的變異性，下面簡要介紹Bootstrap方法基本原理以及基于Bootstrap的巖土體參數(shù)最優(yōu)邊緣分布函數(shù)和最優(yōu)Copula函數(shù)識別方法。Bootstrap方法最早由Efron[12]于1979年提出。它通過對原始觀測樣本數(shù)據(jù)進(jìn)行有放回地隨機(jī)抽樣獲得大量Bootstrap子樣本，然后基于這些子樣本計(jì)算統(tǒng)計(jì)量的估計(jì)值，最終獲得統(tǒng)計(jì)量的變異系數(shù)及其概率分布。該方法只依賴于給定的巖土體參數(shù)的原始觀測樣本，且能夠充分挖掘原始觀測樣本所攜帶的巖土體參數(shù)總體信息，不需要對巖土體參數(shù)的實(shí)際分布作任何假設(shè)和增加新的數(shù)據(jù)觀測。盡管簡單，但該方法的理論依據(jù)及其良好的收斂性早已被統(tǒng)計(jì)學(xué)家所證明[12]。Bootstrap方法不僅被證明對于大部分統(tǒng)計(jì)量滿足大樣本的相合性，而且對于小樣本分析更具有不可替代的優(yōu)越性?；贐ootstrap的統(tǒng)計(jì)量AIC 值的變異性模擬主要包括以下3步：

（1）令巖土體參數(shù)的原始觀測樣本數(shù)據(jù)為X0={(x1i，x2i)，i=1，2，···，N}，從中有放回地隨機(jī)抽樣N次，每次抽取巖土體參數(shù)的一次觀測，從而得到第1個(gè)與X0相同樣本容量的Bootstrap子樣本

（2）重復(fù)上述步驟B 次，即可獲得B 個(gè)Bootstrap子樣本i=1，2，···，N;b=1，2，···，B}；

（3）基于每個(gè)Bootstrap子樣本估計(jì)巖土體參數(shù)的樣本均值、標(biāo)準(zhǔn)差和Kendall秩相關(guān)系數(shù)，然后采用式（4）和式（5）分別計(jì)算4種備選邊緣分布函數(shù)和Copula函數(shù)的AIC 值，從而得到AIC 值的變異系數(shù)及其概率分布。

一般來說，為了達(dá)到良好的收斂效果，Bootstrap子樣本數(shù)目通常取值較大。如文獻(xiàn)[14]研究表明，B=104時(shí)能夠保證Bootstrap模擬的土性參數(shù)的樣本均值和標(biāo)準(zhǔn)差收斂于總體均值和標(biāo)準(zhǔn)差。因此，本文亦采用B=104模擬AIC 值的變異性?；诿總€(gè)Bootstrap子樣本計(jì)算的4種備選邊緣分布函數(shù)和Copula函數(shù)的AIC 值即可識別出該子樣本的最優(yōu)邊緣分布函數(shù)和最優(yōu)Copula函數(shù)。統(tǒng)計(jì)所有B個(gè)Bootstrap子樣本的4種備選邊緣分布函數(shù)和Copula函數(shù)被識別為最優(yōu)邊緣分布和最優(yōu)Copula的次數(shù)[19－20]，從而得到4種備選邊緣分布函數(shù)和Copula函數(shù)為最優(yōu)邊緣分布和最優(yōu)Copula的權(quán)重系數(shù)。與傳統(tǒng)方法不同，本文在考慮AIC 值變異性基礎(chǔ)上，將識別結(jié)果表示為4種備選邊緣分布函數(shù)，Copula函數(shù)被識別為最優(yōu)邊緣分布和最優(yōu)Copula的權(quán)重系數(shù)集合。這種表述方法充分考慮了小樣本容量巖土體參數(shù)引起的AIC 值的變異性，從而更加合理地描述了巖土體參數(shù)最優(yōu)邊緣分布函數(shù)和最優(yōu)Copula函數(shù)識別結(jié)果。與傳統(tǒng)方法一樣，本文方法得出的識別結(jié)果亦可直接代入巖土工程可靠度分析中，從而反映小樣本容量巖土體參數(shù)引起的AIC 值變異性對巖土結(jié)構(gòu)物可靠度結(jié)果的影響。

4 算 例

4.1 基樁荷載-位移雙曲線參數(shù)試驗(yàn)數(shù)據(jù)

本文以基樁荷載-位移雙曲線參數(shù)試驗(yàn)數(shù)據(jù)為例，研究Bootstrap方法在小樣本巖土體參數(shù)最優(yōu)邊緣分布函數(shù)和最優(yōu)Copula函數(shù)識別中的應(yīng)用。對于基樁荷載-位移關(guān)系的預(yù)測問題，國內(nèi)外研究人員提出了多種基樁荷載-位移關(guān)系曲線模型，如雙曲線模型、冪函數(shù)模型、指數(shù)函數(shù)模型、GM（1，1）模型等，其中雙曲線模型在基樁中應(yīng)用最為廣泛。原因在于雙曲線模型中參數(shù)較少，而且物理意義明確。為了降低實(shí)測基樁荷載-位移雙曲線數(shù)據(jù)的離散性，文獻(xiàn)[3－5]提出了采用標(biāo)準(zhǔn)化的荷載-位移雙曲線模型去表征實(shí)測荷載-位移關(guān)系曲線，該標(biāo)準(zhǔn)化的荷載-位移雙曲線模型為

式中：Q為軸向荷載；QSTC為基樁實(shí)測極限承載力；y為樁端位移；a 和b 是雙曲線兩個(gè)參數(shù)，它們具有明確的物理意義，其中a為雙曲線初始斜率的倒數(shù)，b為雙曲線極限值的倒數(shù)，其物理意義見圖1。

圖1 基樁標(biāo)準(zhǔn)化荷載-位移雙曲線模型Fig.1 Definition of the normalized hyperbolic load-settlement curve of foundation pile

考慮文獻(xiàn)[5]中4類基樁荷載-位移雙曲線參數(shù)試驗(yàn)數(shù)據(jù)，這4類數(shù)據(jù)分別是無黏性土中打入樁（D-NC）、無黏性土中灌注樁（B-NC）、黏性土中打入樁（D-C）和黏性土中灌注樁（B-C），它們的樣本容量分別為28、30、59和53，可見這些數(shù)據(jù)都屬于小樣本容量數(shù)據(jù)。圖2給出了上述4類雙曲線參數(shù)數(shù)據(jù)a 和b 的散點(diǎn)圖。由這些數(shù)據(jù)計(jì)算的雙曲線參數(shù)的均值、標(biāo)準(zhǔn)差和相關(guān)系數(shù)見表1?？梢钥闯?，雙曲線參數(shù)具有明顯的負(fù)相關(guān)關(guān)系，它們的Kendall秩相關(guān)系數(shù)分別為-0.597、-0.750、-0.740和-0.755?；谶@些數(shù)據(jù)采用本文方法就可以模擬AIC 值的變異性并識別出擬合數(shù)據(jù)概率分布特性和相關(guān)結(jié)構(gòu)最優(yōu)的邊緣分布函數(shù)和Copula函數(shù)，從而建立雙曲線參數(shù)a 和b 的聯(lián)合概率分布模型，為基樁正常使用極限狀態(tài)可靠度分析奠定基礎(chǔ)[8]。

圖2 雙曲線參數(shù)試驗(yàn)數(shù)據(jù)的散點(diǎn)圖Fig.2 Scattering plots of the measured parameters of the hyperbolic curve

表1 雙曲線參數(shù)試驗(yàn)數(shù)據(jù)的統(tǒng)計(jì)特性Table 1 Statistical properties of the measured parameters of the hyperbolic curve

4.2 雙曲線參數(shù)最優(yōu)邊緣分布函數(shù)的識別

在采用本文方法識別雙曲線參數(shù)a 和b 的最優(yōu)邊緣分布函數(shù)之前，首先給出基于傳統(tǒng)方法的最優(yōu)邊緣分布函數(shù)識別結(jié)果。為此，表2給出了基于雙曲線參數(shù)原始觀測樣本數(shù)據(jù)（見圖2）計(jì)算的4種備選邊緣分布函數(shù)的單一AIC 值。可以看出，除了B-C數(shù)據(jù)的參數(shù)a 和D-NC數(shù)據(jù)的參數(shù)b 最優(yōu)邊緣分布函數(shù)分別是極值I型分布和對數(shù)正態(tài)分布外，其余數(shù)據(jù)參數(shù)a 和b 的最優(yōu)邊緣分布都是威布爾分布?？梢?，采用威布爾分布能夠較好地描述雙曲線參數(shù)試驗(yàn)數(shù)據(jù)的概率分布特性。

下面采用本文方法識別雙曲線參數(shù)a 和b 的最優(yōu)邊緣分布函數(shù)。在給出AIC 值的概率分布之前，首先研究小樣本數(shù)據(jù)對樣本均值和標(biāo)準(zhǔn)差的影響。圖3給出了基于Bootstrap模擬的參數(shù)a 的樣本均值和標(biāo)準(zhǔn)差概率密度函數(shù)。為了比較，圖中還給出了基于參數(shù)a 的原始觀測數(shù)據(jù)計(jì)算的樣本均值和標(biāo)準(zhǔn)差（見豎直線）。可以看出，小樣本巖土體參數(shù)數(shù)據(jù)的樣本均值和標(biāo)準(zhǔn)差具有較大的離散性，如D-NC、B-NC、D-C和B-C數(shù)據(jù)的樣本均值變異系數(shù)分別為0.10、0.14、0.07和0.10，它們的樣本標(biāo)準(zhǔn)差變異系數(shù)分別為0.14、0.11、0.10和0.22?？偟膩碚f，小樣本巖土體參數(shù)數(shù)據(jù)對樣本標(biāo)準(zhǔn)差變異性的影響要大于對樣本均值的影響，這符合統(tǒng)計(jì)學(xué)理論的一般規(guī)律。此外，基于參數(shù)a 的原始觀測數(shù)據(jù)計(jì)算的樣本均值和標(biāo)準(zhǔn)差僅僅是樣本均值和標(biāo)準(zhǔn)差實(shí)際分布的單一取值，該取值能夠較好地被Bootstrap方法所模擬。因此，Bootstrap方法除了保留原始觀測數(shù)據(jù)的相關(guān)特性外，還能充分挖掘原始觀測數(shù)據(jù)所隱藏的總體信息，這也是Bootstrap方法被廣泛應(yīng)用于統(tǒng)計(jì)量變異性模擬的原因所在。同理，采用本文方法可以得出參數(shù)b 的樣本均值和標(biāo)準(zhǔn)差的概率密度函數(shù)，限于篇幅，這里不再列出。小樣本巖土體參數(shù)數(shù)據(jù)樣本均值和標(biāo)準(zhǔn)差的這種變異性將直接導(dǎo)致計(jì)算的邊緣分布函數(shù)AIC 值具有較大變異性，這是因?yàn)橛?jì)算AIC 值時(shí)各備選邊緣分布函數(shù)的分布參數(shù)都是基于樣本均值和標(biāo)準(zhǔn)差確定的。

圖3 參數(shù)a 的樣本均值和標(biāo)準(zhǔn)差概率密度函數(shù)Fig.3 PDFs of mean and standard deviation of a

為了研究小樣本巖土體參數(shù)試驗(yàn)數(shù)據(jù)對AIC 值變異性的影響，圖4以參數(shù)a為例給出了4種備選邊緣分布函數(shù)AIC 值的概率密度函數(shù)。小樣本巖土體參數(shù)試驗(yàn)數(shù)據(jù)對AIC 值變異性具有重要的影響，如B-C數(shù)據(jù)計(jì)算的正態(tài)分布、對數(shù)正態(tài)分布、極值I型分布和威布爾分布AIC 值的變異系數(shù)分別為0.10、0.14、0.08和0.08。不同邊緣分布函數(shù)AIC值的概率密度函數(shù)重疊區(qū)域較大，這進(jìn)一步說明識別最優(yōu)邊緣分布函數(shù)時(shí)4種備選邊緣分布都有可能被識別為最優(yōu)的邊緣分布函數(shù)。傳統(tǒng)方法基于單一AIC 值識別最優(yōu)的邊緣分布函數(shù)忽略了其余備選邊緣分布函數(shù)為最優(yōu)邊緣分布的概率，因此，是不合理的。此外，基于參數(shù)a 的原始觀測數(shù)據(jù)計(jì)算的邊緣分布函數(shù)AIC 值是實(shí)際分布的單一取值，該取值亦能較好地被Bootstrap方法所模擬。

為了充分考慮不同備選邊緣分布為最優(yōu)邊緣分布的概率，表3給出了104個(gè)Bootstrap子樣本中4種備選邊緣分布被識別為最優(yōu)邊緣分布的次數(shù)及權(quán)重系數(shù)?？梢钥闯?，由于小樣本巖土體參數(shù)試驗(yàn)數(shù)據(jù)的影響，沒有一種備選邊緣分布能夠被100%地識別為最優(yōu)的邊緣分布函數(shù)。除了B-C數(shù)據(jù)的參數(shù)a 和D-NC數(shù)據(jù)的參數(shù)b 權(quán)重系數(shù)最大邊緣分布分別是極值I型分布和對數(shù)正態(tài)分布外，其余數(shù)據(jù)參數(shù)a 和b 的權(quán)重系數(shù)最大邊緣分布都是威布爾分布。上述結(jié)果與傳統(tǒng)方法識別結(jié)果保持了較好的一致性。然而，傳統(tǒng)方法沒有考慮其余備選邊緣分布函數(shù)為最優(yōu)邊緣分布的概率。如對于D-NC數(shù)據(jù)的參數(shù)b 來說，除去對數(shù)正態(tài)分布后其余3種備選邊緣分布函數(shù)為最優(yōu)邊緣分布的概率達(dá)到了64.40%。因此，本文提出的基于Bootstrap的巖土體參數(shù)最優(yōu)邊緣分布函數(shù)識別方法比傳統(tǒng)方法要優(yōu)越。

圖4 不同邊緣分布計(jì)算的參數(shù)a 的AIC 值概率密度函數(shù)Fig.4 PDFs of AIC values of various margins for a

表2 基于雙曲線參數(shù)試驗(yàn)數(shù)據(jù)計(jì)算的4種備選邊緣分布函數(shù)AIC 值Table 2 AIC values associated with various margin distribution functions using the measured data for hyperbolic curve-fitting parameters

表3 10 000個(gè)Bootstrap子樣本中不同邊緣分布被識別為最優(yōu)邊緣分布的次數(shù)及權(quán)重系數(shù)Table 3 Numbers of times and weight factors being the best-fit one for different margin distribution functions over 10 000 bootstrap samples

4.3 雙曲線參數(shù)最優(yōu)Copula函數(shù)的識別

在采用本文方法識別雙曲線參數(shù)a 和b 間最優(yōu)Copula函數(shù)前，首先給出基于傳統(tǒng)方法的最優(yōu)Copula函數(shù)識別結(jié)果。為此，表4給出了基于雙曲線參數(shù)原始觀測數(shù)據(jù)計(jì)算的4種備選Copula函數(shù)的單一AIC 值?？梢钥闯觯瑢τ贒-NC數(shù)據(jù)來說，Gaussian Copula函數(shù)計(jì)算的AIC 值最小，它是擬合D-NC數(shù)據(jù)間相關(guān)結(jié)構(gòu)最優(yōu)的Copula函數(shù)。而對于B-NC、D-C和B-C數(shù)據(jù)來說，最優(yōu)的Copula函數(shù)都是Plackett Copula函數(shù)。因此，大部分情況下采用Plackett Copula函數(shù)能夠較好地描述雙曲線參數(shù)試驗(yàn)數(shù)據(jù)的相關(guān)結(jié)構(gòu)。

表4 基于雙曲線參數(shù)試驗(yàn)數(shù)據(jù)計(jì)算的4種備選Copula函數(shù)AIC 值Table 4 AIC values associated with various Copula functions using measured data for hyperbolic curve-fitting parameters

下面采用本文方法識別雙曲線參數(shù)a 和b 間最優(yōu)Copula函數(shù)。在給出AIC 值的概率分布之前，首先研究小樣本巖土體參數(shù)試驗(yàn)數(shù)據(jù)對相關(guān)系數(shù)的影響。圖5給出了基于Bootstrap模擬的參數(shù)a 和b間Kendall秩相關(guān)系數(shù)的概率密度函數(shù)。圖中同時(shí)給出了基于參數(shù)a 和b 的原始觀測樣本數(shù)據(jù)計(jì)算的Kendall秩相關(guān)系數(shù)?？梢钥闯觯颖編r土體參數(shù)試驗(yàn)數(shù)據(jù)的相關(guān)系數(shù)亦具有較大的離散性，如D-NC、B-NC、D-C和B-C數(shù)據(jù)的相關(guān)系數(shù)變異系數(shù)分別為0.19、0.12、0.06和0.08。此外，基于參數(shù)a 和b 的原始觀測樣本數(shù)據(jù)計(jì)算的Kendall秩相關(guān)系數(shù)僅僅是相關(guān)系數(shù)實(shí)際分布的單一取值，該取值亦能夠較好地被Bootstrap方法所模擬。Bootstrap方法不僅能夠重復(fù)原始觀測樣本數(shù)據(jù)的Kendall秩相關(guān)系數(shù)，它還能模擬相關(guān)系數(shù)的變異性及其概率分布。小樣本巖土體參數(shù)試驗(yàn)數(shù)據(jù)相關(guān)系數(shù)的這種變異性將直接導(dǎo)致計(jì)算的Copula函數(shù)AIC值具有較大變異性，這是因?yàn)橛?jì)算AIC 值時(shí)各備選Copula函數(shù)的相關(guān)參數(shù)都是基于Kendall秩相關(guān)系數(shù)采用式（3）確定的。

圖5 雙曲線參數(shù)間相關(guān)系數(shù)的概率密度函數(shù)Fig.5 PDFs of correlation coefficient between a and b

圖6給出了4種備選Copula函數(shù)計(jì)算的AIC值概率密度函數(shù)?？梢钥闯?，小樣本巖土體參數(shù)試驗(yàn)數(shù)據(jù)對AIC 值變異性具有重要的影響，如樣本容量最小的D-NC數(shù)據(jù)計(jì)算的Gaussian、Plackett、Frank和No.16 Copula函數(shù)AIC 值的變異系數(shù)分別為0.40、0.41、0.42和0.65?？梢姸SCopula函數(shù)AIC 值變異系數(shù)明顯大于一維邊緣分布函數(shù)AIC 值變異系數(shù)。與邊緣分布函數(shù)一樣，不同Copula函數(shù)AIC 值的概率密度函數(shù)重疊區(qū)域較大，識別最優(yōu)Copula函數(shù)時(shí)4種備選Copula函數(shù)都有可能被識別為最優(yōu)的Copula函數(shù)。傳統(tǒng)方法忽略了其余備選Copula函數(shù)為最優(yōu)Copula的概率，因此，是不合理的。同理，基于參數(shù)a 和b 的原始觀測數(shù)據(jù)計(jì)算的不同備選Copula函數(shù)AIC值是實(shí)際分布的單一取值，該取值能較好地被Bootstrap所模擬。此外，Bootstrap方法還能有效模擬AIC 值變異性及其概率分布，這是傳統(tǒng)方法所無法比擬的。

為了考慮不同備選Copula函數(shù)為最優(yōu)Copula的概率，表5給出了104個(gè)Bootstrap子樣本中4種備選Copula函數(shù)被識別為最優(yōu)Copula的次數(shù)及權(quán)重系數(shù)。由于小樣本巖土體參數(shù)試驗(yàn)數(shù)據(jù)的影響，4種備選Copula函數(shù)都不能被100%地識別為最優(yōu)的Copula函數(shù)。對于D-NC數(shù)據(jù)而言，權(quán)重系數(shù)最大的Copula函數(shù)是Gaussian Copula。對于D-C數(shù)據(jù)來說，權(quán)重系數(shù)最大的Copula函數(shù)是Frank Copula。而B-NC和B-C數(shù)據(jù)權(quán)重系數(shù)最大的Copula函數(shù)都是Plackett Copula。除了樣本容量最大的D-C數(shù)據(jù)外，上述結(jié)果與傳統(tǒng)方法識別結(jié)果亦保持了較好的一致性。由于傳統(tǒng)方法沒有考慮AIC 值變異性，由傳統(tǒng)方法識別的D-C數(shù)據(jù)最優(yōu)的Plackett Copula函數(shù)的權(quán)重系數(shù)只有17.87%，遠(yuǎn)低于非最優(yōu)的Frank Copula函數(shù)的66.14%。可見，傳統(tǒng)方法容易給出錯(cuò)誤的識別結(jié)果。此外，傳統(tǒng)方法亦沒有考慮其余備選Copula函數(shù)為最優(yōu)Copula的概率。如對于B-NC數(shù)據(jù)來說，除去Plackett Copula函數(shù)后其余3種備選Copula函數(shù)為最優(yōu)Copula的概率達(dá)到了54.01%。因此，本文提出的基于Bootstrap的巖土體參數(shù)最優(yōu)Copula函數(shù)識別方法遠(yuǎn)比傳統(tǒng)方法優(yōu)越。

圖6 不同Copula函數(shù)計(jì)算的AIC 值概率密度函數(shù)Fig.6 PDFs of AIC values of various Copulas functions

表5 10 000個(gè)Bootstrap子樣本中不同Copula函數(shù)被識別為最優(yōu)Copula的次數(shù)及權(quán)重系數(shù)Table 5 Numbers of times and weight factors being best-fit Copula function for different Copula functions over 10 000 bootstrap samples

5 結(jié) 論

（1）Bootstrap方法能夠有效地模擬統(tǒng)計(jì)量的變異性及其概率分布，該方法只依賴于巖土體參數(shù)的原始觀測樣本數(shù)據(jù)，且能夠充分挖掘原始觀測樣本所攜帶的巖土體參數(shù)總體信息，不需要對巖土體參數(shù)的實(shí)際分布作任何假設(shè)以及增加新的數(shù)據(jù)觀測，它為小樣本容量巖土體參數(shù)統(tǒng)計(jì)量的變異性模擬提供了一條有效的途徑。

（2）提出的基于Bootstrap的最優(yōu)邊緣分布函數(shù)和最優(yōu)Copula函數(shù)識別方法，不僅可以有效地考慮統(tǒng)計(jì)量AIC 值的變異性，而且能夠綜合地反映不同備選邊緣分布函數(shù)和Copula函數(shù)為最優(yōu)邊緣分布和最優(yōu)Copula的概率，相比傳統(tǒng)的基于單一AIC值的最優(yōu)邊緣分布函數(shù)和最優(yōu)Copula函數(shù)識別方法要優(yōu)越，建議優(yōu)先采用。

（3）基于小樣本容量巖土體參數(shù)試驗(yàn)數(shù)據(jù)估計(jì)的樣本均值、標(biāo)準(zhǔn)差和相關(guān)系數(shù)具有較大的離散性，這種離散性進(jìn)一步導(dǎo)致了識別最優(yōu)邊緣分布函數(shù)和最優(yōu)Copula函數(shù)時(shí)統(tǒng)計(jì)量AIC 值存在較大的變異性。為了提高巖土體參數(shù)最優(yōu)邊緣分布函數(shù)和最優(yōu)Copula函數(shù)的識別精度，建議盡可能多地收集巖土體參數(shù)試驗(yàn)數(shù)據(jù)。

（4）盡管Bootstrap方法只依賴于給定的巖土體參數(shù)原始觀測樣本數(shù)據(jù)，該原始觀測樣本的樣本容量以及數(shù)據(jù)質(zhì)量仍對最優(yōu)邊緣分布函數(shù)和最優(yōu)Copula函數(shù)的識別結(jié)果有重要的影響。為了使得識別結(jié)果收斂于總體情況，建議采取合理措施保證原始觀測樣本數(shù)據(jù)的代表性和準(zhǔn)確性。

[1]范明橋，盛金保.土強(qiáng)度指標(biāo)φ，c 的互相關(guān)性[J].巖土工程學(xué)報(bào)，1997，19(4):100－104.FAN Ming-qiao，SHENG Jin-bao.Cross correlation of soil strength indexes φ and c[J].Chinese Journal of Geotechnical Engineering，1997，19(4):100－104.

[2]TANG X S，LI D Q，CHEN Y F，et al.Improved knowledge-based clustered partitioning approach and its application to slope reliability analysis[J].Computers and Geotechnics，2012，45:34－43.

[3]PHOON K K，KULHAWY F H.Characterization of model uncertainties for laterally loaded rigid drilled shafts[J].Geotechnique，2005，55(1):45－54.

[4]DITHINDE M，PHOON K K，DE WET M，et al.Characterization of model uncertainty in the static pile design formula[J].Journal of Geotechnical and Geoenvironmental Engineering，2011，137(1):70－85.

[5]LI D Q，TANG X S，PHOON K K，et al.Bivariate simulation using Copula and its application to probabilistic pile settlement analysis[J].International Journal for Numerical and Analytical Methods in Geomechanics，2013，37(6):597－617.

[6]PHOON K K，SANTOSO A，QUEK S T.Probabilistic analysis of soil-water characteristic curves[J].Journal of Geotechnical and Geoenvironmental Engineering，2010，136(3):445－455.

[7]NELSEN R B.An introduction to Copulas[M].New York:Springer，2006.

[8]唐小松，李典慶，周創(chuàng)兵，等.基于Copula函數(shù)的基樁荷載-位移雙曲線概率分析[J].巖土力學(xué)，2012，33(1):171－178.TANG Xiao-song，LI Dian-qing，ZHOU Chuang-bing，et al.Probabilistic analysis of load-displacement hyperbolic curves of single pile using Copula[J].Rock and Soil Mechanics，2012，33(1):171－178.

[9]唐小松，李典慶，周創(chuàng)兵，等.基于Copula函數(shù)的抗剪強(qiáng)度參數(shù)間相關(guān)性模擬及邊坡可靠度分析[J].巖土工程學(xué)報(bào)，2012，34(12):2284－2291.TANG Xiao-song，LI Dian-qing，ZHOU Chuang-bing，et al.Modeling dependence between shear strength parameters using Copulas and its effect on slope reliability[J].Chinese Journal of Geotechnical Engineering，2012，34(12):2284－2291.

[10]唐小松，李典慶，周創(chuàng)兵，等.不完備概率信息條件下邊坡可靠度分析方法[J].巖土工程學(xué)報(bào)，2013，35(6):1027－1034.TANG Xiao-song，LI Dian-qing，ZHOU Chuang-bing，et al.Reliability analysis of slopes with incomplete probability information[J].Chinese Journal of Geotechnical Engineering，2013，35(6):1027－1034.

[11]AKAIKE H.A new look at the statistical model identification[J].IEEE Transactions on Automatic Control，1974，19(6):716－723.

[12]EFRON B.Bootstrap methods:another look at the jackknife[J].The Annals of Statistics，1979，7(1):1－26.

[13]謝桂華，張家生，劉榮桂.沉降計(jì)算經(jīng)驗(yàn)系數(shù)的Bootstrap法置信區(qū)間估計(jì)[J].中南大學(xué)學(xué)報(bào)(自然科學(xué)版)，2011，42(9):2843－2847.XIE Gui-hua，ZHANG Jia-sheng，LIU Rong-gui.Estimation of confidence interval of empirical coefficients in calculating foundation settlement by Bootstrap method[J].Journal of Central South University(Science and Technology)，2011，42(9):2843－2847.

[14]LUO Z，ATAMTURKTUR S，JUANG C H.Bootstrapping for characterizing the effect of uncertainty in sample statistics for braced excavations[J].Journal of Geotechnical and Geoenvironmental Engineering，2013，139(1):13－23.

[15]陳立宏，陳祖煜，劉金梅.土體抗剪強(qiáng)度指標(biāo)的概率分布類型研究[J].巖土力學(xué)，2005，26(1):37－40.CHEN Li-hong，CHEN Zu-yu，LIU Jin-mei.Probability distribution of soil strength[J].Rock and Soil Mechanics，2005，26(1):37－40.

[16]蘇永華，何滿潮，孫曉明.大子樣巖土隨機(jī)參數(shù)統(tǒng)計(jì)方法[J].巖土工程學(xué)報(bào)，2001，23(1):117－119.SU Yong-hua，HE Man-chao，SUN Xiao-ming.Approach on asymptotic approximations of polynomials for probability density function of geotechnics random parameters[J].Chinese Journal of Geotechnical Engineering，2001，23(1):117－119.

[17]蔣水華，馮曉波，李典慶，等.邊坡可靠度分析的非侵入式隨機(jī)有限元法[J].巖土力學(xué)，2013，34(8):2347－2354.JIANG Shui-hua，F(xiàn)ENG Xiao-bo，LI Dian-qing，et al.Reliability analysis of slope using non-intrusive stochastic finite element method[J].Rock and Soil Mechanics，2013，34(8):2347－2354.

[18]蔣水華，李典慶，周創(chuàng)兵，等.考慮參數(shù)空間變異性的非飽和土坡可靠度分析[J].巖土力學(xué)，2014，35(9):2569－2578.JIANG Shui-hua，LI Dian-qing，ZHOU Chuang-bing，et al.Reliability analysis of unsaturated slope considering spatial variability[J].Rock and Soil Mechanics，2014，35(9):2569－2578.

[19]SILVA R S，LOPES H F.Copula，marginal distributions and model selection:a Bayesian note[J].Statistics and Computing，2008，18(3):313－320.

[20]HUARD D，éVIN G，F(xiàn)AVRE A C.Bayesian copula selection[J].Computational Statistics &Data Analysis，2006，51(2):809－822.