張 蕾， 唐小松， 李典慶

( 武漢大學 a.水資源與水電工程科學國家重點實驗室；b.水工巖石力學教育部重點實驗室， 湖北 武漢 430072)

土體抗剪強度參數(shù)粘聚力和內(nèi)摩擦角是巖土結構物變形與穩(wěn)定性分析的重要指標，如邊坡穩(wěn)定性分析、地基承載力分析和擋土墻土壓力計算問題等。粘聚力和內(nèi)摩擦角的取值在很大程度上關系到巖土工程設計的安全性和經(jīng)濟性。眾所周知，土體抗剪強度參數(shù)粘聚力和內(nèi)摩擦角間存在明顯的統(tǒng)計負相關性，而且粘聚力和內(nèi)摩擦角大多服從非正態(tài)分布[1～3]。可靠度方法能夠有效地考慮抗剪強度參數(shù)間相關性對巖土工程安全分析的影響。巖土工程可靠度分析時通常需要建立抗剪強度參數(shù)的聯(lián)合分布函數(shù)，相關非正態(tài)抗剪強度參數(shù)聯(lián)合分布函數(shù)的建立需要大量的試驗數(shù)據(jù)，然而實際工程受經(jīng)濟技術條件限制只能獲得有限的試驗數(shù)據(jù)，在大多數(shù)情況下，基于這些數(shù)據(jù)只能獲得抗剪強度參數(shù)的邊緣分布函數(shù)和相關系數(shù)。因此，為了計算簡便，大多數(shù)研究忽略抗剪強度參數(shù)間相關性，即使有少部分研究考慮抗剪強度參數(shù)間相關性也是采用二維正態(tài)分布模型[4]和基于Nataf變換[5]的等效二維正態(tài)分布模型來建立抗剪強度參數(shù)二維分布。雖然二維正態(tài)分布模型可以考慮抗剪強度參數(shù)間相關性，但是它隱含著參數(shù)的邊緣分布也是正態(tài)分布。此外，它隱含著采用Gaussian copula函數(shù)表征抗剪強度參數(shù)間相關性。上述兩點在很大程度上限制了二維正態(tài)分布模型在土體抗剪強度參數(shù)分析中的應用。同時，Gaussian copula是否能夠精確地描述抗剪強度參數(shù)間相關性也有待進一步研究。為此，本文主要研究Copula函數(shù)[6]框架下土體抗剪強度參數(shù)二維分布模型建立方法。

Copula函數(shù)最早由Sklar[7]于1959年提出，他指出，任意一個多元聯(lián)合分布函數(shù)都可以分解為相應的邊緣分布和一個Copula函數(shù)，這個Copula函數(shù)描述了變量間的相關性。Copula函數(shù)不僅能有效考慮變量間的相關性，還能構造具有任意邊緣分布變量的聯(lián)合分布函數(shù)。Copula函數(shù)為相關非正態(tài)變量聯(lián)合分布函數(shù)的構造提供了一種新的途徑。目前，Copula函數(shù)在金融行業(yè)[8, 9]和水文領域[10～12]已得到廣泛的應用。然而Copula函數(shù)在巖土工程中應用還較少，唐小松等[13]采用Copula函數(shù)建立了基樁荷載—位移雙曲線參數(shù)的聯(lián)合分布函數(shù)。Uzielli和Mayne[14]也將Copula函數(shù)用于研究基樁荷載—位移雙曲線參數(shù)的聯(lián)合分布函數(shù)。然而目前還未見有將Copula函數(shù)用于研究土體抗剪強度參數(shù)二維分布模型問題，為此，本文嘗試將Copula函數(shù)用于巖土體抗剪強度參數(shù)二維分布模型分析問題。收集了小浪底水利樞紐工程三組土體抗剪強度參數(shù)試驗數(shù)據(jù)，分別采用A-D檢驗法和最小平方歐氏距離法識別出抗剪強度參數(shù)的最優(yōu)邊緣分布和最優(yōu)Copula函數(shù)，并與抗剪強度參數(shù)二維正態(tài)分布模型進行比較，驗證了基于Copula函數(shù)的抗剪強度參數(shù)二維分布模型的優(yōu)越性，得出了幾點有益的結論。

1 土體抗剪強度參數(shù)二維分布模型

1.1 試驗數(shù)據(jù)

為了便于讀者理解本文方法，采用三組土體抗剪強度參數(shù)的試驗數(shù)據(jù)為例來說明具體的步驟。數(shù)據(jù)來自小浪底水利樞紐工程主壩防滲體1區(qū)挖坑取樣及現(xiàn)場檢測成果[15, 16]，小浪底1區(qū)土料即小浪底大壩壤土斜心墻土料，大壩斜心墻土料來自寺院坡料場，料場開采區(qū)土層厚度10～30 m，上部0～10 m以中粉質(zhì)壤土為主，下部20 m以重粉質(zhì)壤土為主，這類土廣泛分布于黃河中上游地區(qū)，關中地區(qū)以及陜北地區(qū)等。分析該類土的抗剪強度統(tǒng)計特性，對分析小浪底壩體抗?jié)B性和穩(wěn)定性方面具有一定的指導意義。此外，由于填筑過程中土料已被攪拌均勻，故近似地認為數(shù)據(jù)樣本來自同一母體。由于試驗誤差等因素影響，個別原始試驗數(shù)據(jù)存在異常值，本文采用工程上常用的3σ法將這些異常值進行了剔除，首先在坐標面上標出所有的數(shù)據(jù)點，可以看出有的試驗數(shù)據(jù)點遠離回歸直線，當某一數(shù)據(jù)點偏離回歸直線的距離達到一個特定值時，就判定為異常值，處理后的三組試驗數(shù)據(jù)分別為：三軸固結排水剪切試驗數(shù)據(jù)(CD)(63)、三軸固結不排水剪切試驗數(shù)據(jù)(CU)(64)、三軸不排水剪切試驗數(shù)據(jù)(UU)(61)，括號中數(shù)字代表抗剪強度參數(shù)的樣本數(shù)目，具體數(shù)值見表1。

表1 土料三軸試驗抗剪強度參數(shù)

續(xù)表1

1.2 抗剪強度參數(shù)邊緣分布

在建立抗剪強度參數(shù)二維分布模型之前需要確定其邊緣分布函數(shù)。鑒于A-D檢驗法對于小樣本數(shù)據(jù)識別能力強、精度高、更注重尾部權重的優(yōu)點，本文采用該檢驗法確定抗剪強度參數(shù)的邊緣分布函數(shù)。在文獻[1～3]的基礎上，假設抗剪強度參數(shù)邊緣分布服從正態(tài)分布、對數(shù)正態(tài)分布和極值Ⅰ型分布，然后采用A-D檢驗法識別出最優(yōu)的邊緣分布函數(shù)，選取擬合度最大的分布函數(shù)作為抗剪強度參數(shù)的最優(yōu)邊緣分布。

表2列出了上述三組數(shù)據(jù)粘聚力c的邊緣分布的A-D檢驗結果，同時圖1也給出了粘聚力c的經(jīng)驗頻率以及假設分布擬合結果?？梢钥闯觯龖B(tài)分布是CD、CU和UU數(shù)據(jù)粘聚力的最優(yōu)分布。同理，表3列出了上述三組數(shù)據(jù)內(nèi)摩擦角φ的邊緣分布的A-D檢驗結果，圖2給出了內(nèi)摩擦角的經(jīng)驗頻率以及假設分布擬合結果。可以看出，對于CD數(shù)據(jù)而言，對數(shù)正態(tài)分布是擬合內(nèi)摩擦角最優(yōu)的分布函數(shù)；對于CU和UU數(shù)據(jù)來說，相應的最優(yōu)邊緣分布函數(shù)分別為正態(tài)分布和極值Ⅰ型分布?？梢姵Ｓ玫恼龖B(tài)分布并不一定是擬合內(nèi)摩擦角的最優(yōu)分布函數(shù)。

表2 粘聚力c的邊緣分布檢驗結果(α=0.05)

表3 內(nèi)摩擦角φ的邊緣分布檢驗結果(α=0.05)

現(xiàn)有文獻也對巖土體抗剪強度參數(shù)分布進行了分析，如文獻[1]指出黏土質(zhì)圍巖c和φ的最優(yōu)分布概型分別為正態(tài)分布和對數(shù)正態(tài)分布;文獻[2]中也曾對小浪底心墻料c和φ的邊緣分布做了分析，指出對數(shù)正態(tài)分布更適合作為摩擦角的最優(yōu)分布形式，粘聚力的最優(yōu)分布概型為正態(tài)分布；文獻[3]對我國17座水庫土壩290個土樣的固結快剪試驗進行了分析，得到粘聚力為正態(tài)分布占總數(shù)42.3%，對數(shù)正態(tài)的占38.5%，內(nèi)摩擦角為正態(tài)分布占總數(shù)61.5%，對數(shù)正態(tài)的占26.9%，其他的符合極值Ⅰ型分布。然而由于數(shù)據(jù)量有限，各研究成果之間還存在一定的差異，但總體規(guī)律基本是一致的，本文得出的抗剪強度參數(shù)概率分布與文獻[1～3]中結論基本一致。

圖1 參數(shù)c的經(jīng)驗分布及不同分布擬合曲線

圖2 參數(shù)φ的經(jīng)驗分布及不同分布擬合曲線

1.3 基于Copula函數(shù)的抗剪強度參數(shù)二維分布模型

傳統(tǒng)的抗剪強度參數(shù)二維分布模型大多數(shù)采用二維正態(tài)或二維對數(shù)正態(tài)分布，這些模型雖然可以考慮抗剪強度參數(shù)間相關性，但是它們的邊緣分布只能是正態(tài)分布或對數(shù)正態(tài)分布。為了能夠構造出具有任意邊緣分布函數(shù)及任意相關性的變量聯(lián)合分布函數(shù)，本文采用Copula函數(shù)建立土體抗剪強度參數(shù)二維分布模型。Copula函數(shù)是將變量的聯(lián)合分布與其邊緣分布聯(lián)結起來的函數(shù)，本質(zhì)上它也是一種聯(lián)合分布函數(shù)?；贑opula函數(shù)可得土體抗剪強度參數(shù)c和φ的二維聯(lián)合分布函數(shù)為

F(c,φ)=C(F1(c),F2(φ);θ)=C(u1,u2;θ)

(1)

式中：u1=F1(c)和u2=F2(φ)分別為c和φ的邊緣分布函數(shù)；C為Copula函數(shù)；θ為Copula函數(shù)的參數(shù)。當抗剪強度參數(shù)c和φ的邊緣分布函數(shù)及Copula函數(shù)已知時，采用式(1)就可以得到抗剪強度參數(shù)c和φ的二維累積分布函數(shù)。若變量c和φ的概率密度函數(shù)f1(c)和f2(φ)存在且連續(xù)，對式(1)兩邊求導可得變量c和φ的二維聯(lián)合概率密度函數(shù)f(c，φ)：

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

(2)

式中：c(F1(c),F2(φ);θ)=c(u1,u2;θ)=?2C(u1,u2;θ)/?u1?u2為Copula函數(shù)的概率密度函數(shù)。在已知變量c和φ的邊緣分布函數(shù)和Copula函數(shù)的前提下，利用式(1)和式(2)就可以分別構造出變量c和φ的聯(lián)合累積分布函數(shù)和聯(lián)合概率密度函數(shù)。由式(2)以及條件分布的定義，可得變量φ已知的條件下變量c的條件累積分布函數(shù)F(c|φ2≤φ≤φ1)，如式(3)所示：

F(c|φ2≤φ≤φ1)=

(3)

式中：φ1和φ2分別為變量φ取值的上限和下限。同理，可得變量c已知的條件下變量φ的條件累積分布函數(shù)。

Copula函數(shù)的種類很多，且不同Copula函數(shù)描述變量間的相關性是不一樣的。鑒于抗剪強度參數(shù)具有較強的統(tǒng)計負相關性，本文選取Gaussian、Frank、Plackett和No.16 copula函數(shù)描述c和φ間相關關系。上述4種Copula函數(shù)的分布函數(shù)及參數(shù)取值范圍見文獻[6]，限于篇幅，這里不再列出。基于Copula函數(shù)的抗剪強度參數(shù)二維分布模型構造過程中重要的一步是Copula函數(shù)相關參數(shù)θ的估計。目前估計參數(shù)θ的方法很多，本文采用與抗剪強度參數(shù)邊緣分布函數(shù)無關的非參數(shù)法估計θ。這種方法基于Kendall秩相關系數(shù)τ和Copula函數(shù)相關參數(shù)θ間的對應關系求得Copula函數(shù)的相關參數(shù)θ。Kendall秩相關系數(shù)τ與二維Copula函數(shù)C(u1,u2;θ)之間有如下關系[6]：

(4)

特別地，對于Gaussian Copula有更簡潔的關系式：

(5)

因此，當抗剪強度參數(shù)間Kendall秩相關系數(shù)τ已知時，通過求解式(4)所示的積分方程就可得出參數(shù)θ?？辜魪姸葏?shù)的Kendall秩相關系數(shù)τ可由下式計算[6]：

(6)

表4 Copula函數(shù)相關參數(shù)計算結果

2 土體抗剪強度參數(shù)二維分布模型識別

2.1 最優(yōu)Copula函數(shù)的識別

本節(jié)重點是如何在多種Copula函數(shù)中識別出建立抗剪強度參數(shù)二維分布模型的最優(yōu)Copula函數(shù)。識別最優(yōu)Copula函數(shù)有多種方法，如最小平方歐氏距離法、AIC信息準則法、BIC信息準則法及Nash-Sutcliffe模型指數(shù)法等。由于最小平方歐氏距離法計算簡便而且識別精確高，本文采用該方法進行最優(yōu)Copula函數(shù)識別。平方歐氏距離d2為原始觀測數(shù)據(jù)點處Copula分布函數(shù)值與經(jīng)驗累積分布函數(shù)值之差的平方和，計算公式為

(7)

式中：pi和pei分別為樣本點處的理論累積頻率和經(jīng)驗累積頻率；N為原始觀測數(shù)據(jù)數(shù)目。對應于具有最小平方歐氏距離的Copula函數(shù)即認為是擬合抗剪強度參數(shù)實測數(shù)據(jù)最優(yōu)的Copula函數(shù)。

表5 4種Copula函數(shù)的平方歐氏距離計算結果

表5給出了四種Copula函數(shù)的平方歐氏距離結果?？梢钥闯?，對于UU數(shù)據(jù)來說，Plackett copula函數(shù)具有最小的平方歐氏距離，它是擬合UU數(shù)據(jù)相關性最優(yōu)的Copula函數(shù)；而對于CD和CU數(shù)據(jù)來說，最優(yōu)的Copula函數(shù)都是Frank copula函數(shù)。由上述結果可以看出，對于本文研究的三組數(shù)據(jù)來說，Gaussian copula函數(shù)都不是最優(yōu)的Copula函數(shù)，這說明常用的Gaussian copula函數(shù)并不一定是擬合抗剪強度參數(shù)間相關性最優(yōu)的Copula函數(shù)。因此，實際工程中應該基于抗剪強度參數(shù)的實測數(shù)據(jù)盡可能選取表征抗剪強度參數(shù)間相關性的最優(yōu)Copula函數(shù)建立其二維分布模型，盡量避免為了計算簡便而盲目采用Gaussian copula函數(shù)建立抗剪強度參數(shù)二維分布模型。

圖3 不同Copula函數(shù)粘聚力c條件累積分布函數(shù)的比較

為了進一步比較不同Copula函數(shù)構造的抗剪強度參數(shù)模型的差異，圖3給出了采用式(3)計算的CD組實測數(shù)據(jù)在變量φ取均值附近(μφ-0.25σφ≤φ≤μφ+0.25σφ)以及變量φ取較小值(μφ-2σφ≤φ≤μφ-1.5σφ)和較大值(μφ+1.5σφ≤φ≤μφ+2σφ)時變量c的條件累積分布函數(shù)?？梢钥闯?，在變量φ的均值附近不同Copula函數(shù)構造的變量c的條件累積分布函數(shù)差別較小。相反，當變量φ取較小值或較大值時不同Copula函數(shù)構造的c的條件累積分布函數(shù)差異顯著。為了從定量的角度比較不同Copula函數(shù)得出的條件累積分布函數(shù)差異，表6給出了不同Copula函數(shù)的粘聚力c條件累積分布函數(shù)分位數(shù)值。同樣可以看出當變量φ取較小值時，不同Copula函數(shù)得出的變量c的條件累積分布函數(shù)分位數(shù)差別較大，這點尤其應該引起重視，因為水利水電工程中確定抗剪強度參數(shù)設計值都是采用小值平均法。

2.2 基于Copula函數(shù)的抗剪強度參數(shù)二維分布模型與二維正態(tài)分布模型的對比

對于抗剪強度參數(shù)c和φ的二維分布模型，目前最常用的是二維正態(tài)分布[4]。二維正態(tài)分布是中心對稱分布，它假設變量的邊緣分布函數(shù)相同且都為正態(tài)分布?；贑opula函數(shù)的二維分布模型則不受上述條件的限值，它可以建立任意邊緣分布函數(shù)及任意相關性抗剪強度參數(shù)的聯(lián)合分布函數(shù)。為了比較兩種模型優(yōu)缺點，基于三組數(shù)據(jù)計算的平方歐氏距離如表7所示?？梢娀贑opula函數(shù)建立的抗剪強度參數(shù)二維分布模型具有更小的平方歐氏距離，與二維正態(tài)分布模型相比，基于最優(yōu)Copula函數(shù)建立的二維分布模型能更好地擬合原始觀測數(shù)據(jù)。

2.3 邊緣分布函數(shù)對二維分布模型的影響

為了研究抗剪強度參數(shù)邊緣分布對基于Copula函數(shù)建立的二維分布模型的影響，以CD組試驗數(shù)據(jù)為例，圖4給出了抗剪強度參數(shù)邊緣分布如下三種情況時二維分布模型等概率密度線：(1)c和φ均服從正態(tài)分布；(2)c和φ均服從對數(shù)正態(tài)分布；(3)c和φ均服從極值Ⅰ型分布。圖中Kendall秩相關系數(shù)均為-0.384。篇幅所限，這里僅給出了采用Frank copula和No.16 copula建立的二維分布函數(shù)?？梢钥闯?，盡管相關系數(shù)和Copula函數(shù)相同，不同邊緣分布函數(shù)構造的聯(lián)合概率密度函數(shù)存在明顯的差別。因此，在采用Copula函數(shù)建立抗剪強度參數(shù)二維分布模型時，要盡量避免因為簡單或者從經(jīng)驗角度出發(fā)而選取正態(tài)分布，應根據(jù)實測數(shù)據(jù)確定最優(yōu)邊緣分布，進而建立基于Copula函數(shù)的二維分布模型。

表6 不同Copula函數(shù)的粘聚力c條件累積分布函數(shù)分位數(shù)值

圖4 不同邊緣分布函數(shù)構造的抗剪強度參數(shù)聯(lián)合概率密度函數(shù)的等概率密度線

分布類型CDCUUU最優(yōu)Copula0.023(Frank)0.026(Frank)0.016(Placket)二維正態(tài)分布0.0590.0560.083

3 結 論

(1)小浪底水利樞紐工程1區(qū)土料三軸剪切試驗抗剪強度參數(shù)三組實測數(shù)據(jù)的邊緣分布擬合結果表明該區(qū)粘聚力的最優(yōu)分布均為正態(tài)分布，內(nèi)摩擦角的最優(yōu)分布不都是正態(tài)分布?？辜魪姸葏?shù)邊緣分布函數(shù)的選取對抗剪強度參數(shù)二維分布模型具有明顯的影響。

(2)基于Copula函數(shù)的二維分布模型分析方法可以建立具有任意邊緣分布函數(shù)及任意相關性的抗剪強度參數(shù)的聯(lián)合分布函數(shù)。與傳統(tǒng)抗剪強度參數(shù)二維正態(tài)分布模型相比，基于最優(yōu)Copula函數(shù)和最優(yōu)邊緣分布函數(shù)構造的抗剪強度參數(shù)二維分布模型具有更強的擬合原始觀測數(shù)據(jù)的能力，而且適用范圍廣，它是建立相關非正態(tài)巖土體物理力學參數(shù)聯(lián)合分布函數(shù)的一種有效的方法。

(3)常用的Gaussian copula函數(shù)并不一定是表征抗剪強度參數(shù)間相關性的最優(yōu)的Copula函數(shù)。在實際工程中應該盡可能收集更多的抗剪強度參數(shù)試驗數(shù)據(jù)，并基于這些試驗數(shù)據(jù)識別表征抗剪強度參數(shù)間相關性的最優(yōu)的Copula函數(shù)來建立抗剪強度參數(shù)二維分布模型，不能單單是為了計算簡便而盲目采用Gaussian copula函數(shù)建立抗剪強度參數(shù)二維分布模型。

[1] 陳煒韜, 王玉鎖, 王明年,等. 黏土質(zhì)隧道圍巖抗剪強度參數(shù)的概率分布及優(yōu)化實例[J]. 巖石力學與工程學報, 2006, 25(S2): 3782-3787.

[2] 陳立宏, 陳祖煜, 劉金梅. 土體抗剪強度指標的概率分布類型研究[J]. 巖土力學, 2005, 26(1): 37-40.

[3] 范明橋. 粘性填筑土強度指標φ,c的概率特性[J]. 水利水運科學研究, 2000, (1): 49-53.

[4] 趙宇飛, 楊 鈞. 巖土強度參數(shù)二維正態(tài)分布擬合研究[J]. 華北水利水電學院學報, 2006, 27(2): 79-81.

[5] Nataf A. Détermination des distributions de probabilité dont les marges sont données[J]. Comptes Rendus de l’Académie des Sciences，1962, 225: 42-43.

[6] Nelsen R B. An Introduction to Copulas[M]. New York: Springer Verlag, 2006.

[7] Sklar A. Fonctions de répartition à n dimensions et leurs marges[J]. Publications de l’Institut de Statistique de l’Université de Paris, 1959, 8: 229-231.

[8] Kole E, Koedijk K, Verbeek M. Selecting copulas for risk management[J]. Journal of Banking and Finance, 2007, 31(8): 2405-2423.

[9] 張堯庭. 連接函數(shù)(Copula)技術與金融風險分析[J]. 統(tǒng)計研究, 2002, (4): 48-51.

[10] Grimaldi S, Serinaldi F. Asymmetric copula in multivariate flood frequency analysis[J]. Advances in Water Resources, 2006, 29(8): 1155-1167.

[11] Zhang L, Singh V P. Bivariate rainfall frequency distributions using archimedean copulas[J]. Journal of Hydrology, 2007, 332(1-2): 93-109.

[12] Zhang Q, Chen Y D, Chen X, et al. Copula-based analysis of hydrological extremes and implications of hydrological behaviors in the Pearl River basin, China[J]. Journal of Hydrologic Engineering, 2011, 16(7): 598-607.

[13] 唐小松, 李典慶, 周創(chuàng)兵,等. 基于Copula函數(shù)的基樁荷載-位移雙曲線概率分析[J]. 巖土力學, 2012, 33(1): 171-178.

[14] Uzielli M，Mayne P W. Load-displacement uncertainty of vertically loaded shallow footings on sands and effects on probabilistic settlement[J]. Georisk, 2012, 6(1): 50-69.

[15] 于清波. 土石壩壩坡穩(wěn)定可靠度研究及其工程應用[D]. 南京: 河海大學, 2006.

[16] 王奇峰, 王 碧. 小浪底高土石壩心墻土料應用探討[J]. 西北水電, 2001, (1): 65-69.