基于Copula理论的粗粒土渗透破坏临界水力比降估值

黄 达 ，曾 彬顾东明

（1.重庆大学 土木工程学院，重庆 400045；2.重庆大学 山地城镇建设与新技术教育部重点实验室，重庆 400045）

1 引 言

渗透稳定性是土力学的重要研究内容之一，而渗透破坏临界水力比降的确定是进行土体渗透稳定性分析和渗流控制的基础[1]。砂砾类无黏性粗粒土作为常用的坝体填筑材料和良好的土石坝地基土，其渗透破坏的临界水力比降理论计算方法大致可归纳为3类：①基于渗透力概念，按照其与土体（或土颗粒）浮重间的极限平衡推导出的临界水力比降公式，具有代表性的有Terzaghi公式[2]、刘杰公式[3]、沙金煊公式[4]、毛昶熙公式[5]等，其中前3种计算方法被有关规范[6]推荐使用；②基于毛管模型的土颗粒渗透破坏临界水力比降公式，以Kovacs[7]提出的毛管模型为基础，Indraratna等[8－9]推导了砂砾类土颗粒的临界水力比降公式，刘忠玉等[10]推导了考虑细颗粒流失量的临界水力比降公式；③基于试验数据的概率统计方法，如葛祖立[11]采用多元回归分析方法统计了418个试验结果，建立了计算临界水力比降的两个经验公式。由于粗粒土的渗透破坏机制复杂，临界水力比降离散程度大，不同计算公式得到的临界水力比降差异较大[12]，因此，试验方法仍是当前确定粗粒土临界水力比降最直接、最可靠的方法，也是绝大多数理论或经验取值方法的重要基础。充分利用有限的试验数据，采用概率统计等数学方法探求土体组构与临界水力比降统计相关性，进而估计类似土体的临界水力比降，这也是一种科学实用的途径。

目前，在概率统计分析方面，一般采用多维正态分布方法，即变量的边缘分布均为正态分布，但正态分布并不一定为小样本条件下土体物理力学参数分布的最优估计，甚至有时不能通过检验[13]。Sklar[14]提出的Copula理论为构建多变量联合分布提供了一种很好途径。该理论指出，可以将任何一有限维的联合分布分解为它的边缘分布和一个表示结构关系的Copula函数，其中Copula函数描述变量间的相关性和一致性[14]。采用Copula理论建立相关变量间的相互关系，优点在于可将边缘分布和结构关系分开研究，这样就可构造边缘分布为任意分布的多变量联合分布函数，使问题大大简化。Copula理论已经在国内外金融、保险以及水文等领域[15]广泛应用，如Shiau[16]利用Copula函数建立了干旱历时和干旱烈度的联合分布，Song等[17]基于Plackett Copula函数建立了干旱历时、干旱烈度以及干旱间隔时间之间的联合分布，并计算了相应的条件重现期。Copula理论在岩土工程中应用的报道并不多见，唐小松等[18]利用Copula函数对基桩荷载-位移双曲线进行了概率分析。

本文以朱崇辉[19]开展的砂砾类粗粒土室内渗透变形试验数据为基础，采用Copula理论构造其临界水力比降与孔隙比及级配特征参数（不均匀系数和曲率系数）的四维多变量联合分布函数，为建立粗粒土临界水力比降与孔隙比及级配特征的多变量概率相关关系和临界水力比降的估值提供了一种新的途径。

2 渗透变形试验及变量间相关性度量

2.1 渗透变形试验

以陕西杨凌渭河砂砾类粗粒土为研究对象，朱崇辉[19]开展了粗粒土渗透破坏临界水力比降 Jcr与土体孔隙比e、级配不均匀系数 Cu和曲率系数Cc的单因素渗透变形试验研究。试验装置包括渗透仪、供水设备和量测设备，渗透变形试验示意图如图1所示。渗透仪主要由直径约为100 mm的有机玻璃圆形试验筒、透水板及底座组成，试样高度为80 mm。试验样本颗粒粒径范围为0.1～20.0 mm。通过控制e、Cu和Cc中任两个变量为常数，改变另一参数进行颗粒人工级配（含水率控制在约5%～8%）。3组31次试验土体级配及孔隙比与相应 Jcr如表1所示。

图1 渗透变形试验示意图Fig.1 Schematic diagram of seepage-induced deformation test

表1 单因素渗透变形试验结果[19]Table 1 Results of seepage-induced deformation tests with single factor[19]

2.2 变量间相关性度量

在用Copula函数描述变量间的相关性结构之前，须进行变量间相关性分析，以考察各变量之间的相关程度。相关系数是度量随机变量之间相关性的常用指标，有Pearson线性相关系数γn、Spearman秩相关系数ρn和Kendall秩相关系数τn等[15]，其计算公式分别为

表2为2.1节试验（表1）的相关性系数计算结果，可看出：粗粒土的 Jcr与e、Cu、Cc均存在较好的负相关性。

第二届中国京菜节已经落下帷幕，其为振兴发展中国京菜奠定了更加坚实的基础。今后，在北京市商务局的指导下、在有关方面的支持下，北京烹饪协会将进一步实施振兴发展中国京菜行动计划，立足当下、着眼未来，整体规划、突出重点。在发展京菜的同时，还要高度重视京味餐饮企业和全国地方菜餐饮企业在京的发展。让我们继续努力，提升餐饮品质稳增长，为促进北京餐饮行业的健康有序发展而努力奋斗。

表2 孔隙比及级配特征参数与临界水力比降间的相关系数Table 2 Correlation coefficients of critical hydraulic gradient with void ratio and characteristic parameters of grain gradation

3 Copula联合分布函数构造

3.1 Copula理论简介

Copula函数是定义在[0，1]区间均匀分布的联合分布函数。设 X1，…，Xn为具有边缘分布 F1，…，Fn的随机变量，其联合分布函数为F (x1，x2，…，xn)，u1=F1(x1)，…，un=Fn(xn)，则存在一个n元Copula函数C，使得任意 xi∈ Rn，有[15]

式中：x1，x2，…，xn分别为随机变量 X1，X2…，Xn的观测值；u1，u2…， un分别为 X1，X2…，Xn的边缘分布函数。

特别地，当 X1，…，Xn相互独立或它们之间的相依性可以忽略不计时，则有

所以，对于相互独立的多个变量，其联合分布等于它们各自边缘分布的乘积，而无需利用Copula函数构建它们的联合分布。由此可看出，Copula函数对于构建多个非独立变量间的联合分布函数，特别是变量间的相关性不易明确地用数学表达式表示时，更能显现出其优势。显然，e、Cu、Cc三者之间并非独立。例如，Cu、Cc均是有效粒径 d10和限定粒径 d60的函数，e 与其颗粒组成也有密切关系[20]。采用传统方法很难将4个变量间的整体关系表示出来，而采用Copula理论则可以实现。

Copula理论之所以受到大量研究者的青睐，还在于它可以将边缘分布和变量间相关结构分开来研究，这样便大大减小了多变量概率模型建模和分析的难度。利用Copula函数构造变量的联合分布函数，第1步是估计各变量的边缘分布函数，该函数可以是任意形式的；第2步是选择最优的Copula函数来描述变量间的相关关系。

3.2 小样本单变量边缘分布推断及检验

正态信息扩散法基于信息扩散原理，强调充分利用样本数据信息，而非预先假定样品符合某经典的概率分布曲线（拟合检验），因此，数学和物理意义更为严密。正态信息扩散法是估计小样本条件下岩土物理力学参数概率密度函数较好的方法[13]。该方法提出的随机变量X 概率密度函数 f (x) 的正态信息扩散估计为

式中：h为标准正态扩散函数的窗宽。

对于随机变量X 的观测值 xi，其最大值为 xmax，最小值为 xmin，根据正态信息扩散的择近原则，可求得h=β(xmax-xmin)/(n-1)[21]。β 是与n 相关的常数，可查表获得，具体参见文献[21]。

K-S检验法[15]（即Kolmogorov-Smirnov法）是一种适用于样本数较少情况下的概率分布类型的检验方法，当计算的随机变量Dn小于具有显著水平α 的上限值Dn，α时，便认为通过检验，具体检验过程参见文献[15]。利用K-S检验法分别对由正态信息扩散法得到的各变量的概率分布进行检验，结果如表4所示。各变量Dn均小于具有显著水平0.05的上限值Dn，α，因此，由正态信息扩散法得到的各变量的概率密度函数都通过拟合良好性检验，能够真实地反映随机变量的概率分布情况。

表3 系数及1/(2h2)值Table3 Valuesof and1/(2h2)forcoefficients

表3 系数及1/(2h2)值Table3 Valuesof and1/(2h2)forcoefficients

表4 K-S检验法结果Table 4 Results of Kolmogorov-Smirnov test procedure

3.3 最优Copula函数的构造与检验

关于最优Copula函数构造常用的方法是选用不同的Copula函数拟合数据，然后对拟合结果进行拟合优度检验[12]。具体的步骤如下：

（1）根据变量的边缘分布和相关结构等性质选取符合条件的Copula函数。考虑到粗粒土 Jcr与e、Cu和 Cc的相关关系，本文选用单参数的四维对称Archimedean Copula函数来描述变量间的相关结构，其函数可定义为

式中：φ (·)为Archimedean Copula函数的生成函数；φ-1为 φ (·) 的逆函数；u1、u2、u3、u4分别对应e、Cu、Cc和 Jcr的边缘分布函数。

Archimedean Copula函数很多，Nelsen[23]系统地介绍了Archimedean Copula函数的构造方法，Matteis[24]在Nelsen[23]的基础上总结了20种对称Archimedean Copula函数的生成函数。本文即选取Matteis给出的20种生成函数推导相应的四维对称Archimedean Copula函数。

（2）Copula函数相关参数的估计，常用的参数估计方法有非参数法和极大似然法等[24]，对于四维Copula函数，采用非参数法进行参数估计不再适用[24]，因此，选用一般性的极大似然法对参数进行估计。构造Copula函数的对数似然函数为

式中：样本容量n=31；x1i、x2i、x3i、x4i和 F1(x1i)、F2(x2i)、F3(x3i)、F4(x4i)分别为e、Cu、Cc、Jcr的样本观测值和对应的边缘分布函数值；θ为Copula函数的待估参数；c为C (u1，u2，u3，u4)的密度函数，按下式计算：

式中：L为式（6）确定的对数似然函数；Lmax为L函数所取得的最大值。

根据式（6）～（8）分别计算 20种对称Archimedean Copula四维函数参数θ 的估计值，结果见表5。

（3）拟合优度检验，一般是将各个观测点数据的经验联合分布和Copula函数理论联合分布进行比较，计算均方根误差RMSE（root mean square errors）、AIC（Akaike information criteria）[25]等指标，指标值最小的即为拟合原始观测数据的最优Copula函数，工程上一般采用AIC 值进行判断[26]，其计算公式如下：

式中：k为Copula函数中相关参数的个数；C (ui，vi，wi，ti)为随机变量(xi，yi，zi，ri)的理论概率，其中ui、 vi、 wi和 ti分别为对应随机变量 xi、 yi、 zi和ri的边缘分布函数，即ui=Fx(xi)、vi=Fy(yi)、wi=Fz(zi)、ti=Fr(ri)；Femp(xi，yi，zi，ri)为随机变量(xi，yi，zi，ri)的经验概率，计算式为

表5 四维对称Archimedean Copula函数及参数θ、AIC 值计算结果Table 5 Four dimensional symmetrical Archimedean Copula functions and calculation results of parameters θ，AIC

式中：nl，m，k，s为同时满足X ≤ xi，Y ≤ yi，Z ≤ zi，R ≤ ri时联合观测值的个数。

根据式（9）、（10）计算各函数的拟合优度评价指标AIC 值，结果列于表5，可看出Nelsen No 6函数建立的联合分布模型具有最小的AIC 值。结合式（6）～（8）估计的θ 值，确定描述粗粒土 Jcr与e、Cu和 Cc的最优Copula函数为

图2为Nelsen No 6四维最优Copula函数式（11）计算得到各试验点的理论概率值 C (ui，vi，wi，ti)与经验概率值 Femp(xi，yi，zi，ri)的对比图。可以看出，其基本分布在45°对角线附近，表明所建立的联合分布函数拟合效果良好。

图2 理论与经验概率的比较Fig.2 Empirical probabilities versus theoretical ones

4 临界水力比降估值

4.1 估值的保证率

粗粒土 Jcr、e、Cu和 Cc间的相关关系通过构造的Nelsen No 6四维最优Copula函数式（11）来描述。定义其 Jcr估值为 Jcr0条件下，实际值小于或等于 Jcr0的概率为保证率 α(0≤ α≤ 1)，即

显然，当保证率α 过大时，估值偏大；α 较小时，估值偏低。在已知e=e0、Cu=Cu0和 Cc=Cc0情况下求Nelsen No 6四维最优Copula函数的条件概率可得到保证率α ：

式（13）所表示的意义为：在e、Cu和 Cc已知条件下，Jcr估值为 Jcr0的保证率(α)，即利用式（13），一方面可以计算其他临界水力比降估值方法的保证率α ；另一方面，则可计算在一定保证率条件下 Jcr的估值。

4.2 最优拟合保证率

根据式（13），在保证率α 一定的情况下可求得临界水力比降的估值 Jcr0。不同的α 值将得到不同的 Jcr0值，显然存在一个 α0值，使得临界水力比降估值与试验值拟合程度最优。统计学中可用残差平方和来衡量估计值与试验值的拟合度，越小，估值拟合程度越高。其计算公式为

基于最优Copula函数式（11），利用式（13）、（14），计算不同α 值条件下表1中各试验点和。在α 取值范围内（0≤ α≤ 1），取得最小值时所对应的α 值即为Copula函数最优拟合的保证率 α0。计算表明，当α=0.459时，取得最小值为2.993 4，即式（11）所构造的Copula函数的最优拟合保证率为0.459。

4.3 估值及比较

为验证本文所构造的最优Copula函数用于类似粗粒土临界水力比降估值的合理性，补充了5个砂砾类粗粒土室内渗透变形试验。试验样品的颗粒级配曲线如图3所示。e、Cu、Cc如表6所示。

利用构造的Nelsen No 6四维最优Copula函数在最优拟合保证率α0=0.459的条件下对5组试样进行了临界水力比降的估值。为了对比分析，参考《水利水电工程地质勘查规范》[5]，采用Terzaghi公式[2]计算流土破坏的临界水力比降，采用刘杰公式[3]计算管涌或过渡型破坏的临界水力比降，并与试验值比较，计算结果见表6。

图3 渗透变形试验试样颗粒级配曲线Fig.3 Curves of grain gradation for test specimens

表6 试样级配特征参数及试验临界水力比降与理论估值Table 6 Characteristic parameters of grain gradation，calculating and lab’ vaules of critical hydraulic gradients

从表可以看出：对于流土破坏的粗粒土，Copula方法与Terzaghi估值总体上偏小；而对于管涌或过渡型破坏的粗粒土，其估值与试验值接近，估值效果优于刘杰公式。由此可见，Copula理论方法用于粗粒土临界水力比降估值总体上效果良好。

5 结 论

（1）具有单参数的四维对称Archimedean Copula函数的Nelsen No 6是拟合粗粒土临界水力比降与孔隙比、级配不均匀系数和曲率系数间相关关系的最优Copula函数。

（2）利用构造的最优Copula函数求条件概率，便可得到粗粒土临界水力比降估值方法的保证率，或者计算在一定保证率条件下临界水力比降的估值。

（3）通过室内渗透变形试验及其他理论计算结果的对比分析，利用构造的最优Copula函数对粗粒土临界水力比降的估值合理。在已知粗粒土孔隙比及级配参数条件下，为估算临界水力比降提供了一种新途径。

[1]毛昶熙，段祥宝，毛佩郁，等.提防渗流与防冲[M].北京:中国水利水电出版社，2003.

[2]TERZAGHI K.Theoretical soil mechanics[M].New York:Wiley，1943.

[3]刘杰.土的渗透稳定与渗流控制[M].北京:水利电力出版社，1992.

[4]沙金煊.多孔介质中的管涌研究[J].水利水运科学研究，1981，(3):89－93.

[5]毛昶熙，段祥宝，吴良骥.砂砾土各级颗粒的管涌临界坡降研究[J].岩土力学，2009，30(12):3705－3709.MAO Chang-xi，DUAN Xiang-bao，WU Liang-ji.Study of critical gradient of piping for various grain sizes in sandy gravels[J].Rock and Soil Mechanics，2009，30(12):3705－3709.

[6]中华人民共和国水利部.GB50487－2008水利水电工程地质勘察规范[S].北京:中国计划出版社，2009.

[7]KOVACS G.Seepage hydraulics[M].Amsterdam:Elsevier Scientific Publishing Company，1981.

[8]INDRARATNA B，VAFAI F.Analytical model for particle migration within base soil-filter system[J].Journal of Geotechnical and Geoenvironmental Engineering，1997，123(2):100－109.

[9]INDRARATNA B，ASCE M，RADAMPOLA S.Analysis of critical hydraulic gradient for particle movement in filtration[J].Journal of Geotechnical and Geoenvironmental Engineering，2002，128(4):347－350.

[10]刘忠玉，乐金朝，苗天德.无黏性土中管涌的毛管模型及其应用[J].岩石力学与工程学报，2004，23(22):3871－3876.LIU Zhong-yu，YUE Jin-chao，MIAO Tian-de.Capillarytube model for piping in noncohesive soils and its application[J].Chinese Journal of Rock Mechanics and Engineering，2004，23(22):3871－3876.

[11]葛祖立.砂砾管涌类型的判别和临界坡降计算[J].水利水电技术，1987，(5):36－41.

[12]刘黎.粗粒料渗透特性及渗透规律试验研究[D].成都:四川大学，2006.

[13]宫凤强，李夕兵，邓建.小样本岩土参数概率分布的正态信息扩散法推断[J].岩石力学与工程学报，2006，25(12):2559－2564.GONG Feng-qiang，LI Xi-bing，DENG Jian.Probability distribution of small samples of geotechnical parameters using normal information spread method[J].Chinese Journal of Rock Mechanics and Engineering，2006，25(12):2559－2564.

[14]SKLAR A.Fonctions de répartition àn dimensions et leurs magres[J].Publications de l’Institut de Statistique de l’Université de Paris，1959，8:229－231.

[15]宋松柏，蔡焕杰，金菊良，等.Copula函数及其在水文中的运用[M].北京:科学出版社，2012.

[16]SHIAU J T.Fitting drought duration and severity with two-dimensional Copulas[J].Management，2006，20(5):795－815.

[17]SONG SONGBAI，VIJAY P S.Frequency analysis of droughts using the Plackett Copula and parameter estimation by genetic algorithm[J].Stochastic Environmental Research and Risk Assessment，2010，24(5):783－805.

[18]唐小松，李典庆，周创兵，等.基于Copula函数的基桩-荷载位移双曲线概率分析[J].岩土力学，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.

[19]朱崇辉.粗粒土的渗透特性研究[D].杨凌:西北农林科技大学，2006.

[20]VUKOVIC M，SORO A.Determination of hydraulic conductivity of porous media from grain-size composition[M].Littleton:Water Resources Publications，1992.

[21]王新洲.基于信息扩散原理的估计理论、方法及其抗差性[J].武汉测绘科技大学学报，1999，24(3):240－244.WANG Xin-zhou.The theory，method and robustness of the parameter estimation based on the principle of information spread[J].Journal of Wuhan Technical University of Surveying and Mapping，1999，24(3):240－244.

[22]樊妮，赫孝良，赵谦.基于Bayesian思想的最优Copula函数选择[J].工程数学学报，2012，29(4):516－522.FAN Ni，HE Xiao-liang，ZHAO Qian.Optimal Copula function selection based on Bayesian methodology[J].Chinese Journal of Engineering Mathematics，2012，29(4):516－522.

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

[24]DE MATTEIS R.Fitting Copulas to data[D].Zurich:Diploma Thesis，Institute of Mathematics of the University Zurich，2001.

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

[26]唐小松，李典庆，周创兵，等.基于Copula函数的抗剪强度参数间相关性模拟及边坡可靠度分析[J].岩土工程学报，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.