岩土参数概率分布的系统推断研究*

孔令奇 李翠娟

(1.成都工业学院材料与环境工程学院， 成都 610031； 2.西南交通大学土木工程学院， 成都 610031)

岩土参数的概率分布模型是研究边坡可靠度的前提[1]。近年来，研究人员一直致力于岩土参数概率分布的推断方法研究。对于岩土参数这样的工程随机变量，其样本概率分布总体的估计无外乎分为参数估计法和非参数估计法[2]。基于参数估计法推断岩土参数概率密度函数的研究已经取得了一些重要结论——岩土参数不严格服从正态分布[3-5]。利用正态分布和对数正态分布表示岩土参数概率分布存在经典分布的定义区间和岩土参数实际分布区间不匹配的问题[6-8]。采用正态、对数正态等无界分布会使结果失去实际工程意义[4]。beta 分布可有效地反映岩土参数的分布特点,是经典分布有限范围内的最优分布[4-5]。研究人员发现某些情况下参数法拟合岩土参数概率分布时存在较大误差[9]，故而提出非参数拟合方法。文献[10]系统考察了非参数拟合的方法如：最大熵法[1,11-12]、正态信息扩散法(NIS)[13]和正交多项式法[15]在推断岩土参数概率分布时的优劣，最终给出正态信息扩散法为相对最优推断方法的结论。文献[14]进一步证明正态信息扩散法在大样本岩土参数概率分布推断中的有效性，但该方法在应用的过程中存在窗宽选择的难题[16]。

迄今，研究人员提出了众多的岩土参数概率分布推断方法，对于工程人员来说，面对一组新测试样本时仍存在着选择哪种理论分布推断样本概率分布的难题。

为此，对于大量的岩土参数测试样本，将借助数据挖掘技术中的模糊C-均值聚类方法，探寻岩土参数概率分布的统计规律；对于一组新测试的岩土参数样本，基于其数字特征值和判别分析方法，将其分类，获得概率分布模型；解析岩土参数分布特征对正态信息扩散法最优窗宽的影响，提出正态信息扩散法最优窗宽的确定方法。为岩土参数概率分布推断提供简单、有效的方法。

1 岩土参数统计特征值

对任意的随机变量，除了可用其概率密度函数来表示外，还可以用各阶矩来刻画该随机变量[17]。如果能确定岩土参数的概率密度函数和相应的数字特征，就可全面地描述岩土参数的概率特征[2,17]。大量的研究工作主要关注岩土参数概率密度函数的获得，忽视了其主要数值特征分布规律的研究。工程中可以用岩土参数的某些关键量或主要数值特征量去近似地描述它的概率特征，这些特征值需能表征岩土参数测试样本的波动性、离散型和分布的不对称性[2,17]。为此，定义了有效系数，为2阶原点矩与平均值之比，用来表征岩土参数测试样本的波动程度；变异系数表征岩土参数样本的离散性，为样本的方差与均值之比，一般当变异系数cv>0.2时，说明岩土参数样本的离散性很严重[2]；偏度系数表征岩土参数分布的不对称性。当样本数量较大时，可以采用更高阶矩来描述其数字特征，由于岩土参数测试样本数量通常较小，属于小样本，高阶矩将失去意义[2]。

为了从统计中寻找规律，收集了19组岩土参数测试样本，来说明统计的思路和得出的关键结论。19组岩土参数测试样本数据如表1所示。

工程中，当各个随机变量的量纲和数量级不一致时，为了探寻海量数据内在的规律性，需对样本进行变换处理，以消除各样本量纲和数量级的限制。首先对岩土参数测试样本进行了极差归一化变换处理[21]，即：

(1)

表1内收集得到的岩土参数测试样本采样于不同的地质单元，且采用的测试设备及测量控制方法等也不尽相同[22]。对表1内的样本不难获得其统计特征值(如表3所示)，但很难从统计特征值中直观地判断采用何种理论分布推断其概率密度函数。

表1 岩土参数测试样本

2 岩土参数的聚类分析

2.1 模糊C-均值聚类

要从众多的数据中抽取岩土参数固有的特征，获得简洁的系统行为，模糊C-均值聚类分析是解决上述问题的有效手段之一。

传统的系统聚类和K均值聚类法存在把每个测试样本严格地划分到某个类中的局限性，Bezdek借助模糊理论，提出模糊C-均值聚类方法[23]。该方法通过迭代计算不断修正隶属度矩阵，极小化所有测试样本到聚类中心的距离和隶属度的加权，最终通过优化目标获得每个测试样本对所有类中心的隶属度，实现对大量测试样本进行分类的目的。

设矩阵X由n个测试样本(x1x2…xn)的p个特征值构成。

(2)

c个类的聚类中心记为V={v1,v2,…,vc}，其中vi=(vi1,vi2,…,vic)，i=1,2,…,c，2≤c≤n。

(3)

其中dik=‖xk-vi‖

式中：U=(uik)c×n为隶属度矩阵。

模糊C-均值聚类法的具体步骤如下[17,21]：

2)计算迭代至第l步的聚类中心V(l)：

(4)

3)修正U(l)，计算J(l)。

(5a)

i=1,2,…,c;k=1,2,…,c

(5b)

4)判断收敛的条件:

(6)

式中：εu为终止容限。

如满足收敛条件，终止迭代，否则l=l+1，然后转至第2步。

2.2 岩土参数聚类

针对表1的19组数据进行模糊C-均值聚类分析。分类的个数取决于聚类数形图和不一致系数[17]，将19组测试样本分成4类，迭代收敛的条件取εu=10-6。19组岩土参数的隶属度矩阵如表2所示。

表2 隶属度矩阵

由表2可见：如1号测试样本压缩系数属于第4类的隶属度最大为0.926 9，所以将其归为第4类，其他测试样本分类方法相同。表3给出了归聚为一类的岩土参数测试样本的统计特征值。

由表3明显可见，归聚为一类的岩土参数测试样本数字特征值相近，表明其概率分布相近。由概率论和数理统计可知，直方图以很高的概率接近随机变量的概率密度函数[2]，归聚为一类的岩土参数测试样本的直方图如图1～4所示。

a—样本序号2; b—样本序号7; c—样本序号9; d—样本序号11; e—样本序号13。——beta分布； ----正态信息扩散法。

表3 聚类之后的样本的统计特征值

3 岩土参数概率分布选择

3.1 推断岩土参数的概率密度函数

大量岩土参数测试样本进行分类后，目标是获得各类岩土参数测试样本的概率分布形式。在既有研究成果的基础上，采用的岩土参数概率密度拟合方法为beta分布和正态信息扩散法。

3.1.1标准beta分布

标准beta分布为[17]

(7)

式中：α1、α2为B分布参数。

文献[17]提出根据标准beta分布的均值、方差与分布参数α1、α2的关系，同时考虑概率性质和岩土参数的数字特征，只要获得岩土参数测试样本的平均值和有效值，即可获得beta分布的参数α1和α2。

3.1.2正态信息扩散法

正态信息扩散法的概率密度函数为：

(8)

式中：K()为核函数；h为窗宽；n为样本个数；X1,X2,…,Xn为取自总体X的样本。

详细的原理过程见文献[6,24-25]，窗宽h采用文献[14]给出的窗宽计算式：

(9)

式中：δ为样本观测值的标准差。

3.2 拟合效果评价指标

K-S检验法的检验量是整个取值范围内的最大偏差值，且与显著性水平有关，不能提供某一分布拟合是良好的绝对信息[2]。为了从全局角度评估beta分布和正态信息扩散法对岩土参数概率密度函数的拟合质量，采用均方差指标ηMSE对拟合曲线进行效果校验[17]：

a—样本序号4; b—样本序号6; c—样本序号8; d—样本序号16; e—样本序号17; f—样本序号18; g—样本序号19。——beta分布； ----正态信息扩散法。

a—样本序号3; b—样本序号5; c—样本序号14; d—样本序号15。——beta分布； ----正态信息扩散法。

a—样本序号1; b—样本序号10; c—样本序号12。——beta分布； ----正态信息扩散法。

(10)

式中：fi为理论频率；S*(xi)为频率直方图中每个子区间的相对频率；m为区间个数。

3.3 岩土参数测试样本的概率密度函数推断

3.3.1各类样本概率密度拟合效果分析

图1～4分别给出了第1～4类岩土参数测试样本beta分布和正态信息扩散法拟合效果。

由图1～4可见，归聚为一类的岩土参数样本概率分布形状基本相同。表4给出了岩土参数测试样本的beta分布参数和正态信息扩散法的窗宽。表5给出了两种拟合方法的误差分析。其中表4中的h指按式(9)计算所得值。

表4 拟合参数求解结果

表5 拟合误差分析

由表3～5可见：1)第1类测试样本有效系数和偏度系数均较小；第2类测试样本对比第1类测试样本，有效系数相当，偏度系数较大；这两类测试样本反映的概率分布特征为波动性较小，这两类岩土参数测试样本采用beta分布和正态信息扩散法进行拟合均有效，误差均在0.1以内[17]，且正态信息扩散法拟合最优窗宽均为h。2)第3类测试样本的有效系数和偏度系数明显变大，测试样本反映的概率分布特征为波动性较大，不对称性明显。第3类测试样本采用beta分布拟合时误差均大于0.1，正态信息扩散法拟合误差小于0.1，所以，对于第3类测试样本推荐采用正态信息扩散法进行拟合其概率密度函数，最优窗宽仍为h。3)对于第4类测试样本，显著特征是有效系数更大，变异系数四类中最大，测试样本反映的概率分布特征为波动性、离散性严重。beta分布拟合误差更大，当采用正态信息扩散法拟合其概率密度函数时，式(9)给出的窗宽h显然已不能体现样本的较大波动性。

3.3.2正态信息扩散法最优窗宽确定

正态信息扩散法的窗宽选择一直是个难题，为此，提出对波动性较大的岩土参数测试样本窗宽选择方法，具体步骤如下：

1)按式(9)计算窗宽h的初始值；2)依次改变窗宽的大小为h/2、h/3、h/4、…，进行粗选窗宽，设为h′，粗选窗宽的确定依据为均方差最小；3)确定粗选窗宽h′后，以更小的步长微调粗选窗宽h′(文中以0.000 5为步长，增大、减小窗宽h′)，确定最优窗宽h″，最优窗宽的确定依据为均方差最小。

仿真分析发现，对于波动性较大的测试样本，工程误差允许范围内，可选择粗选窗宽h′为最优窗宽。

对于第4类测试样本，窗宽取为h和寻优之后的窗宽h″的对比曲线如图5所示。

a—样本序号1; b—样本序号10; c—样本序号12。

当窗宽取为h时，3组测试样本的均方差分别为0.239 8、0.466 4、0.293 1。当窗宽取为h″时，均方误差分别为0.093 8、0.036 3、0.043 4。证明本文方法的有效性。

4 工程实例分析

测试样本取自南昌地铁二号线翠苑路站[18]，该地铁二号线粉质黏土压缩系数参数30个测试样本：27.3，18.6，22.9，23.7，21.9，20.7，22.7，21.7，23.1，22.0，23.0，21.5，21.5，22.8，21.3，23.7，25.2，23.0，23.1，24.1，21.5，23.7，22.2，25.6，20.4，23.6，22.3，20.7，22.5，19.0 MP-1。通过计算，该测试样本经过极差归一化处理后的有效系数为1.1，变异系数为0.21，偏度系数为0.25。根据判别分析方法，判断该组测试样本属于第1类。采用beta分布和正态信息扩散法都可推断其概率密度函数。采用beta分布推断其概率分布，测试样本的平均值为0.449 4，有效值为0.493 2，求得beta分布的参数α1=2.560 6和α2=2.471 1。获得beta分布的概率密度函数为：

(11)

绘制其频率直方图和beta分布拟合曲线如图6所示。均方差为0.017 1。

图6 beta分布拟合曲线

为验证方法的有效性，提取文献[18]提供的23个黏聚力测试样本：14.515 5，14.645 0，14.671，14.877 9，15.018 7，15.232 9，15.257 5，15.269 8，15.287，15.377 5，15.395 3，15.402 7，15.495 5，15.587 6，15.635，15.654 1，15.710 8，16.003，16.029 2，16.121，16.141 5，16.878 7，17.142 3 kPa。通过计算，该测试样本经过极差归一化处理后的有效系数为1.18，变异系数为0.25，偏度系数为0.78。采用判别分析方法，判断该组测试样本属于第3类。采用正态信息扩散法推断其概率密度函数，最优窗宽为h=0.139 8。采用正态信息扩散法拟合其概率分布如图7所示，为了对比，亦给出beta分布拟合曲线。

——beta分布； —·—正态信息扩散法。

正态信息扩散分布法拟合曲线的均方差指标为0.087 2，beta分布拟合曲线的均方差指标为0.317 9。

以上工程实例表明，对于一组新的测试样本，只要先计算出其数值特征值，即可简单、有效地确定采用何种理论模型推断其概率分布。

5 结束语

通过模糊C-均值聚类方法分析，对19组岩土参数测试样本进行分类统计分析，统计结果能有效地展现岩土参数测试样本概率分布选择的规律性，结论如下：

1)基于参数估计的beta分布只适用于推断波动性较小的岩土参数测试样本的概率分布(如第1类和第2类)。

2)正态信息扩散法当窗宽取为h时，仅适用于推断波动性和离散性不严重的岩土参数测试样本的概率分布(如第1、2、3类)；对于波动性、离散性严重的岩土参数测试样本，正态信息扩散法在推断其概率分布时最优窗宽需根据具体样本寻优选择，采用较小的窗宽才能反映出样本的波动性。

限于篇幅，仅以参数估计法中的beta分布和非参数估计的正态信息扩散法两种典型概率模型为例，其他概率分布的判断过程与此相同。采用聚类分析后，岩土参数测试样本的概率分布推断不再是一个复杂的、反复试错的过程，整个工作将变得系统化、程序化。