基于Copula函数的土体抗剪强度参数二维分布模型

张 蕾， 唐小松， 李典庆

( 武汉大学 a.水资源与水电工程科学国家重点实验室；b.水工岩石力学教育部重点实验室， 湖北 武汉 430072)

土体抗剪强度参数粘聚力和内摩擦角是岩土结构物变形与稳定性分析的重要指标，如边坡稳定性分析、地基承载力分析和挡土墙土压力计算问题等。粘聚力和内摩擦角的取值在很大程度上关系到岩土工程设计的安全性和经济性。众所周知，土体抗剪强度参数粘聚力和内摩擦角间存在明显的统计负相关性，而且粘聚力和内摩擦角大多服从非正态分布[1～3]。可靠度方法能够有效地考虑抗剪强度参数间相关性对岩土工程安全分析的影响。岩土工程可靠度分析时通常需要建立抗剪强度参数的联合分布函数，相关非正态抗剪强度参数联合分布函数的建立需要大量的试验数据，然而实际工程受经济技术条件限制只能获得有限的试验数据，在大多数情况下，基于这些数据只能获得抗剪强度参数的边缘分布函数和相关系数。因此，为了计算简便，大多数研究忽略抗剪强度参数间相关性，即使有少部分研究考虑抗剪强度参数间相关性也是采用二维正态分布模型[4]和基于Nataf变换[5]的等效二维正态分布模型来建立抗剪强度参数二维分布。虽然二维正态分布模型可以考虑抗剪强度参数间相关性，但是它隐含着参数的边缘分布也是正态分布。此外，它隐含着采用Gaussian copula函数表征抗剪强度参数间相关性。上述两点在很大程度上限制了二维正态分布模型在土体抗剪强度参数分析中的应用。同时，Gaussian copula是否能够精确地描述抗剪强度参数间相关性也有待进一步研究。为此，本文主要研究Copula函数[6]框架下土体抗剪强度参数二维分布模型建立方法。

Copula函数最早由Sklar[7]于1959年提出，他指出，任意一个多元联合分布函数都可以分解为相应的边缘分布和一个Copula函数，这个Copula函数描述了变量间的相关性。Copula函数不仅能有效考虑变量间的相关性，还能构造具有任意边缘分布变量的联合分布函数。Copula函数为相关非正态变量联合分布函数的构造提供了一种新的途径。目前，Copula函数在金融行业[8, 9]和水文领域[10～12]已得到广泛的应用。然而Copula函数在岩土工程中应用还较少，唐小松等[13]采用Copula函数建立了基桩荷载—位移双曲线参数的联合分布函数。Uzielli和Mayne[14]也将Copula函数用于研究基桩荷载—位移双曲线参数的联合分布函数。然而目前还未见有将Copula函数用于研究土体抗剪强度参数二维分布模型问题，为此，本文尝试将Copula函数用于岩土体抗剪强度参数二维分布模型分析问题。收集了小浪底水利枢纽工程三组土体抗剪强度参数试验数据，分别采用A-D检验法和最小平方欧氏距离法识别出抗剪强度参数的最优边缘分布和最优Copula函数，并与抗剪强度参数二维正态分布模型进行比较，验证了基于Copula函数的抗剪强度参数二维分布模型的优越性，得出了几点有益的结论。

1 土体抗剪强度参数二维分布模型

1.1 试验数据

为了便于读者理解本文方法，采用三组土体抗剪强度参数的试验数据为例来说明具体的步骤。数据来自小浪底水利枢纽工程主坝防渗体1区挖坑取样及现场检测成果[15, 16]，小浪底1区土料即小浪底大坝壤土斜心墙土料，大坝斜心墙土料来自寺院坡料场，料场开采区土层厚度10～30 m，上部0～10 m以中粉质壤土为主，下部20 m以重粉质壤土为主，这类土广泛分布于黄河中上游地区，关中地区以及陕北地区等。分析该类土的抗剪强度统计特性，对分析小浪底坝体抗渗性和稳定性方面具有一定的指导意义。此外，由于填筑过程中土料已被搅拌均匀，故近似地认为数据样本来自同一母体。由于试验误差等因素影响，个别原始试验数据存在异常值，本文采用工程上常用的3σ法将这些异常值进行了剔除，首先在坐标面上标出所有的数据点，可以看出有的试验数据点远离回归直线，当某一数据点偏离回归直线的距离达到一个特定值时，就判定为异常值，处理后的三组试验数据分别为：三轴固结排水剪切试验数据(CD)(63)、三轴固结不排水剪切试验数据(CU)(64)、三轴不排水剪切试验数据(UU)(61)，括号中数字代表抗剪强度参数的样本数目，具体数值见表1。

表1 土料三轴试验抗剪强度参数

续表1

1.2 抗剪强度参数边缘分布

在建立抗剪强度参数二维分布模型之前需要确定其边缘分布函数。鉴于A-D检验法对于小样本数据识别能力强、精度高、更注重尾部权重的优点，本文采用该检验法确定抗剪强度参数的边缘分布函数。在文献[1～3]的基础上，假设抗剪强度参数边缘分布服从正态分布、对数正态分布和极值Ⅰ型分布，然后采用A-D检验法识别出最优的边缘分布函数，选取拟合度最大的分布函数作为抗剪强度参数的最优边缘分布。

表2列出了上述三组数据粘聚力c的边缘分布的A-D检验结果，同时图1也给出了粘聚力c的经验频率以及假设分布拟合结果。可以看出，正态分布是CD、CU和UU数据粘聚力的最优分布。同理，表3列出了上述三组数据内摩擦角φ的边缘分布的A-D检验结果，图2给出了内摩擦角的经验频率以及假设分布拟合结果。可以看出，对于CD数据而言，对数正态分布是拟合内摩擦角最优的分布函数；对于CU和UU数据来说，相应的最优边缘分布函数分别为正态分布和极值Ⅰ型分布。可见常用的正态分布并不一定是拟合内摩擦角的最优分布函数。

表2 粘聚力c的边缘分布检验结果(α=0.05)

表3 内摩擦角φ的边缘分布检验结果(α=0.05)

现有文献也对岩土体抗剪强度参数分布进行了分析，如文献[1]指出黏土质围岩c和φ的最优分布概型分别为正态分布和对数正态分布;文献[2]中也曾对小浪底心墙料c和φ的边缘分布做了分析，指出对数正态分布更适合作为摩擦角的最优分布形式，粘聚力的最优分布概型为正态分布；文献[3]对我国17座水库土坝290个土样的固结快剪试验进行了分析，得到粘聚力为正态分布占总数42.3%，对数正态的占38.5%，内摩擦角为正态分布占总数61.5%，对数正态的占26.9%，其他的符合极值Ⅰ型分布。然而由于数据量有限，各研究成果之间还存在一定的差异，但总体规律基本是一致的，本文得出的抗剪强度参数概率分布与文献[1～3]中结论基本一致。

图1 参数c的经验分布及不同分布拟合曲线

图2 参数φ的经验分布及不同分布拟合曲线

1.3 基于Copula函数的抗剪强度参数二维分布模型

传统的抗剪强度参数二维分布模型大多数采用二维正态或二维对数正态分布，这些模型虽然可以考虑抗剪强度参数间相关性，但是它们的边缘分布只能是正态分布或对数正态分布。为了能够构造出具有任意边缘分布函数及任意相关性的变量联合分布函数，本文采用Copula函数建立土体抗剪强度参数二维分布模型。Copula函数是将变量的联合分布与其边缘分布联结起来的函数，本质上它也是一种联合分布函数。基于Copula函数可得土体抗剪强度参数c和φ的二维联合分布函数为

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

(1)

式中：u1=F1(c)和u2=F2(φ)分别为c和φ的边缘分布函数；C为Copula函数；θ为Copula函数的参数。当抗剪强度参数c和φ的边缘分布函数及Copula函数已知时，采用式(1)就可以得到抗剪强度参数c和φ的二维累积分布函数。若变量c和φ的概率密度函数f1(c)和f2(φ)存在且连续，对式(1)两边求导可得变量c和φ的二维联合概率密度函数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函数的概率密度函数。在已知变量c和φ的边缘分布函数和Copula函数的前提下，利用式(1)和式(2)就可以分别构造出变量c和φ的联合累积分布函数和联合概率密度函数。由式(2)以及条件分布的定义，可得变量φ已知的条件下变量c的条件累积分布函数F(c|φ2≤φ≤φ1)，如式(3)所示：

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

(3)

式中：φ1和φ2分别为变量φ取值的上限和下限。同理，可得变量c已知的条件下变量φ的条件累积分布函数。

Copula函数的种类很多，且不同Copula函数描述变量间的相关性是不一样的。鉴于抗剪强度参数具有较强的统计负相关性，本文选取Gaussian、Frank、Plackett和No.16 copula函数描述c和φ间相关关系。上述4种Copula函数的分布函数及参数取值范围见文献[6]，限于篇幅，这里不再列出。基于Copula函数的抗剪强度参数二维分布模型构造过程中重要的一步是Copula函数相关参数θ的估计。目前估计参数θ的方法很多，本文采用与抗剪强度参数边缘分布函数无关的非参数法估计θ。这种方法基于Kendall秩相关系数τ和Copula函数相关参数θ间的对应关系求得Copula函数的相关参数θ。Kendall秩相关系数τ与二维Copula函数C(u1,u2;θ)之间有如下关系[6]：

(4)

特别地，对于Gaussian Copula有更简洁的关系式：

(5)

因此，当抗剪强度参数间Kendall秩相关系数τ已知时，通过求解式(4)所示的积分方程就可得出参数θ。抗剪强度参数的Kendall秩相关系数τ可由下式计算[6]：

(6)

表4 Copula函数相关参数计算结果

2 土体抗剪强度参数二维分布模型识别

2.1 最优Copula函数的识别

本节重点是如何在多种Copula函数中识别出建立抗剪强度参数二维分布模型的最优Copula函数。识别最优Copula函数有多种方法，如最小平方欧氏距离法、AIC信息准则法、BIC信息准则法及Nash-Sutcliffe模型指数法等。由于最小平方欧氏距离法计算简便而且识别精确高，本文采用该方法进行最优Copula函数识别。平方欧氏距离d2为原始观测数据点处Copula分布函数值与经验累积分布函数值之差的平方和，计算公式为

(7)

式中：pi和pei分别为样本点处的理论累积频率和经验累积频率；N为原始观测数据数目。对应于具有最小平方欧氏距离的Copula函数即认为是拟合抗剪强度参数实测数据最优的Copula函数。

表5 4种Copula函数的平方欧氏距离计算结果

表5给出了四种Copula函数的平方欧氏距离结果。可以看出，对于UU数据来说，Plackett copula函数具有最小的平方欧氏距离，它是拟合UU数据相关性最优的Copula函数；而对于CD和CU数据来说，最优的Copula函数都是Frank copula函数。由上述结果可以看出，对于本文研究的三组数据来说，Gaussian copula函数都不是最优的Copula函数，这说明常用的Gaussian copula函数并不一定是拟合抗剪强度参数间相关性最优的Copula函数。因此，实际工程中应该基于抗剪强度参数的实测数据尽可能选取表征抗剪强度参数间相关性的最优Copula函数建立其二维分布模型，尽量避免为了计算简便而盲目采用Gaussian copula函数建立抗剪强度参数二维分布模型。

图3 不同Copula函数粘聚力c条件累积分布函数的比较

为了进一步比较不同Copula函数构造的抗剪强度参数模型的差异，图3给出了采用式(3)计算的CD组实测数据在变量φ取均值附近(μφ-0.25σφ≤φ≤μφ+0.25σφ)以及变量φ取较小值(μφ-2σφ≤φ≤μφ-1.5σφ)和较大值(μφ+1.5σφ≤φ≤μφ+2σφ)时变量c的条件累积分布函数。可以看出，在变量φ的均值附近不同Copula函数构造的变量c的条件累积分布函数差别较小。相反，当变量φ取较小值或较大值时不同Copula函数构造的c的条件累积分布函数差异显著。为了从定量的角度比较不同Copula函数得出的条件累积分布函数差异，表6给出了不同Copula函数的粘聚力c条件累积分布函数分位数值。同样可以看出当变量φ取较小值时，不同Copula函数得出的变量c的条件累积分布函数分位数差别较大，这点尤其应该引起重视，因为水利水电工程中确定抗剪强度参数设计值都是采用小值平均法。

2.2 基于Copula函数的抗剪强度参数二维分布模型与二维正态分布模型的对比

对于抗剪强度参数c和φ的二维分布模型，目前最常用的是二维正态分布[4]。二维正态分布是中心对称分布，它假设变量的边缘分布函数相同且都为正态分布。基于Copula函数的二维分布模型则不受上述条件的限值，它可以建立任意边缘分布函数及任意相关性抗剪强度参数的联合分布函数。为了比较两种模型优缺点，基于三组数据计算的平方欧氏距离如表7所示。可见基于Copula函数建立的抗剪强度参数二维分布模型具有更小的平方欧氏距离，与二维正态分布模型相比，基于最优Copula函数建立的二维分布模型能更好地拟合原始观测数据。

2.3 边缘分布函数对二维分布模型的影响

为了研究抗剪强度参数边缘分布对基于Copula函数建立的二维分布模型的影响，以CD组试验数据为例，图4给出了抗剪强度参数边缘分布如下三种情况时二维分布模型等概率密度线：(1)c和φ均服从正态分布；(2)c和φ均服从对数正态分布；(3)c和φ均服从极值Ⅰ型分布。图中Kendall秩相关系数均为-0.384。篇幅所限，这里仅给出了采用Frank copula和No.16 copula建立的二维分布函数。可以看出，尽管相关系数和Copula函数相同，不同边缘分布函数构造的联合概率密度函数存在明显的差别。因此，在采用Copula函数建立抗剪强度参数二维分布模型时，要尽量避免因为简单或者从经验角度出发而选取正态分布，应根据实测数据确定最优边缘分布，进而建立基于Copula函数的二维分布模型。

表6 不同Copula函数的粘聚力c条件累积分布函数分位数值

图4 不同边缘分布函数构造的抗剪强度参数联合概率密度函数的等概率密度线

分布类型CDCUUU最优Copula0.023(Frank)0.026(Frank)0.016(Placket)二维正态分布0.0590.0560.083

3 结 论

(1)小浪底水利枢纽工程1区土料三轴剪切试验抗剪强度参数三组实测数据的边缘分布拟合结果表明该区粘聚力的最优分布均为正态分布，内摩擦角的最优分布不都是正态分布。抗剪强度参数边缘分布函数的选取对抗剪强度参数二维分布模型具有明显的影响。

(2)基于Copula函数的二维分布模型分析方法可以建立具有任意边缘分布函数及任意相关性的抗剪强度参数的联合分布函数。与传统抗剪强度参数二维正态分布模型相比，基于最优Copula函数和最优边缘分布函数构造的抗剪强度参数二维分布模型具有更强的拟合原始观测数据的能力，而且适用范围广，它是建立相关非正态岩土体物理力学参数联合分布函数的一种有效的方法。

(3)常用的Gaussian copula函数并不一定是表征抗剪强度参数间相关性的最优的Copula函数。在实际工程中应该尽可能收集更多的抗剪强度参数试验数据，并基于这些试验数据识别表征抗剪强度参数间相关性的最优的Copula函数来建立抗剪强度参数二维分布模型，不能单单是为了计算简便而盲目采用Gaussian copula函数建立抗剪强度参数二维分布模型。

[1] 陈炜韬, 王玉锁, 王明年,等. 黏土质隧道围岩抗剪强度参数的概率分布及优化实例[J]. 岩石力学与工程学报, 2006, 25(S2): 3782-3787.

[2] 陈立宏, 陈祖煜, 刘金梅. 土体抗剪强度指标的概率分布类型研究[J]. 岩土力学, 2005, 26(1): 37-40.

[3] 范明桥. 粘性填筑土强度指标φ,c的概率特性[J]. 水利水运科学研究, 2000, (1): 49-53.

[4] 赵宇飞, 杨 钧. 岩土强度参数二维正态分布拟合研究[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] 张尧庭. 连接函数(Copula)技术与金融风险分析[J]. 统计研究, 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] 唐小松, 李典庆, 周创兵,等. 基于Copula函数的基桩荷载-位移双曲线概率分析[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] 于清波. 土石坝坝坡稳定可靠度研究及其工程应用[D]. 南京: 河海大学, 2006.

[16] 王奇峰, 王 碧. 小浪底高土石坝心墙土料应用探讨[J]. 西北水电, 2001, (1): 65-69.