摘 要：

最大熵法被认为是可以准确推断岩土参数概率分布的一种有效方法，但其还存在直方图分组和积分界限选择问题。岩土参数测试样本的直方图分组是最大熵法推断其概率分布的前提和基础，既有直方图分组没有统一的方法，一般为研究人员主观取值，为此提出借助经验法中的Sturges规则与平均移位直方图的思想，通过计算平均积分平方误差最小值获得最优直方图分组的方法。既有的c32型区间取值方法在样本偏度系数绝对值大于0.5时失效，通过对实际测试样本的反复统计计算，考虑积分上、下界限的同时变化，提出在样本偏度系数的绝对值大于0.5时，采用c′31积分界限取值方法，从而弥补了c32区间取值法在样本偏度系数大于0.5时失效的不足。研究成果完善了最大熵法在推断岩土参数概率分布中的应用。

关键词:

岩土参数; 最大熵法; 概率分布; 最优直方图分组; 积分界限；

作者简介：

*孔令奇(1980—),男，讲师，硕士，从事岩土工程可靠度和数据挖掘技术研究。E-mail:5929985@qq.com;

基金：

国家自然科学基金资助项目(51308466);

引用：

孔令奇, 李翠娟. 岩土参数概率分布最大熵法直方图分组和积分界限选择问题研究[ J] . 水利水电技术(中英文), 2021, 52 (10): 172- 181.

KONG Lingqi, LI Cuijuan. Study on histogram grouping and integration boundary selection for maximum entropy method of probability distribution of geotechnical parameters[J]. Water Resources and Hydropower Engineering, 2021, 52(10): 172- 181.

0 引 言

岩土参数概率分布的研究一直是一项基础性研究工作，岩土参数的概率分布模型直接影响岩土工程可靠度指标的计算结果，其概率分布的研究具有极其重要的价值和意义。国内外研究人员已做了大量工作,得出了一些指导工程实际的结论：①研究指出正态分布、对数正态分布、贝塔分布等基于参数估计的方法在推断岩土参数概率分布时误差较大，如正态分布存在概率分布的定义区间与实际岩土参数分布区间不匹配的问题;②参数法拟合岩土参数概率分布时存在较大误差,故而提出非参数拟合方法，如正交多项式拟合法、正态信息扩散原理拟合法和最大熵法等，但基于非参数估计的核密度拟合方法存在过于依赖样本数据，可能出现拟合过度的问题。文献[3]指出最大熵原理避免了一些非参数法过分依赖样本数据的缺点，引入的外在主观信息量少，不依赖于研究人员的工程经验，分布精度和稳健性更高，被认为是可以准确推断岩土参数概率分布的一种有效方法。

最大熵法已广泛应用于岩石抗压强度的概率密度函数推断、基桩竖向承载力可靠指标计算及边坡可靠度分析计算中。但对既有研究工作总结分析发现，最大熵法在实际应用过程中，还存在如下问题：①目前采用最大熵法推断岩土参数概率分布过程中，主要将直方图作为评价概率分布拟合优劣的参考基准。但众所周知，直方图分组数的不同，概率密度拟合曲线形状会受到很大的影响。较小的分组会导致曲线不够光滑，较大的分组又会使曲线过度光滑而失去细节。实际岩土参数概率分布推断过程中，同一组岩土参数测试样本，不同研究人员给出的直方图分组数不同,目前岩土参数样本频率直方图分组个数和区间的划分主要是研究人员凭主观决定。②文献[15]采用最大熵法推断岩石力学参数概率分布模型时，指出在利用最大熵原理时，积分过程中必须选择一个合适的积分区间，如果选择不当，可能使最大熵分布的尾部出现一个不应有的升高现象。在实际研究过程中，样本界限的选择一直是个难题，如果上界选择过大，概率密度拟合曲线的尾部就会出现波动，因为尾部样本点分散，取值概率很小；反之，如果上界选择过小，拟合曲线的尾部会上翘，原因在这种界限附近，样本点尚未明显减少，使得这个点以外的点上的取值概率得以累加。如果样本下界选得过大，拟合曲线在首端的取值概率较大，会使曲线变得平坦。文献[11,13]探讨了最大熵法积分区间的确定方法，文中对比分析了3σ、4σ、c33、c31和c32型区间取值方法的优劣，最终确定c32型为最合理的积分区间。但深入研究、分析发现，c32型区间取值方法只考虑分布呈偏态时对样本某一界限的影响，而样本偏度系数较大时，忽略了样本另一界限的改变对分布曲线形状的影响，就会出现左端截尾点为负的情况，这与实际岩土参数的分布不符。因此，基于最大熵法推断岩土参数概率分布模型理论中积分界限的确定仍值得深入研究。

本文在既有研究基础上主要完成如下工作：通过岩土参数实测样本的对比分析，指出既有岩土参数概率密度函数推断过程中直方图分组方法的不足，提出建立合理的直方图分组是岩土参数概率分布推断的前提，本文将经验法中的Sturges规则与平均移位直方图的思想相结合，通过计算平均积分平方误差的最小值获得岩土参数最优直方图分组方法；以实测的岩土参数测试样本为例，探讨既有最大熵法积分界限确定方法存在的不足，指出样本偏度系数较大时，应同时考虑样本另一界限的改变对分布曲线形状的影响，提出样本偏度系数的绝对值大于0.5时积分区间的确定方法，完善最大熵法在岩土参数概率密度拟合上的应用。

1 最大熵推断方法

最大熵法的概率分布估计是一个优化问题

i阶原点矩；f(z)是随机变量z的概率密度函数。

由式(1)可知，通过改变f(z)可使熵S[f(z)]达到最大值。最大熵分布概率密度函数的解析表达式为

式中，λ0,λ1,…,λN为拉格朗日乘子，详细推导过程参见文献[13]和文献[21]。

特别注意，为了避免出现计算zi时发生溢出的可能性，可将z值域变换到[0,1]之间来计算，为此，后文对岩土参数测试样本进行了极差归一化处理,即

通常，采用最大熵法拟合时矩的阶数要取大一点，但是必须注意，如果样本容量比较小，高阶矩的统计值会因误差较大而失去意义。因此一般矩的阶数取4～6时，就可以获得较为满意的拟合结果。对于岩土参数测试样本来说，由于受经济技术条件限制，岩土体参数的试验数据比较有限，样本个数通常不足50个。文献[13]的统计结果也表明，矩的阶数取8阶以上时，拟合误差均较大，矩的阶数为7时的拟合误差与矩的阶数为6时的拟合误差相差不大。故针对岩土参数测试样本，本文将矩的阶数取为5～6。

2 直方图分组方法

2.1 既有直方图分组

直方图是使拟合偏差最小的一种方法，所以目前对岩土参数概率分布的推断通常是基于岩土参数，绘制样本的频率直方图，然后对直方图进行拟合并评价拟合的优劣。

在概率论与数理统计中，对于正态分布总体的随机变量，其直方图子区间的划分与样本数量有最佳关系，取分组数m=1.87(n-1)2/5。但大量的研究工作已表明，岩土参数测试样本离散性严重，呈偏态分布。在既有岩土参数概率分布推断中，直方图分组主要是研究人员凭主观决定。

以文献[13,17,18]中提供的容重和抗拉强度参数测试样本为例，文献[13]将容重参数测试样本直方图分组取为9,文献[17]中容重参数测试样本直方图分组取为11。文献[13]中抗拉强度参数测试样本直方图分组取为9,文献[18]中抗拉强度参数测试样本直方图分组取为8。绘制两组岩土参数测试样本的直方图如图1和图3所示。

由图1和图3可见，直方图分组的不同，表征了岩土参数不同的概率分布形态。采用最大熵法拟合重度和抗拉强度参数样本直方图，拟合曲线如图1—图4所示。

图1 容重参数样本最大熵拟合曲线(矩的阶数为5)

图2 容重参数样本最大熵拟合曲线(矩的阶数为6)

图3 抗拉强度参数样本最大熵拟合曲线(矩的阶数为5)

图4 抗拉强度参数样本最大熵拟合曲线(矩的阶数为6)

为了从全局角度评价拟合方法的优劣，在此引入一个无量纲的判定系数，其克服了K-S检验法仅通过最大值进行局部判断的不足。判定系数定义如下

式中，yi为实测值；y¯为实测值的均值；Yi为拟合估计值。

由式(3)可知，若实测值yi与拟合估计值Yi之间误差越小，则判定系数越接近1,拟合效果也就越好。容重和抗拉强度参数测试样本拟合检验值如表1和表2所列。

由表1和表2可见：①最大熵法矩的阶数确定，直方图分组不同，拟合误差不同。如容重参数样本矩的阶数为5时，两种不同直方图分组下拟合误差差别较大；抗拉强度参数样本在矩的阶数为6时，误差差别也较大。②直方图分组不同，与之匹配的最大熵法矩的阶数不同。对于容重参数样本直方图分组为11时，显然最大熵法矩的阶数为5时与之匹配的误差最小；而直方图分组为9时，最大熵法矩的阶数为6时与之匹配的误差最小。抗拉强度参数样本直方图分组为9时，最大熵法矩的阶数为5时与之匹配的误差最小；而直方图分组为8时，最大熵法矩的阶数为6时与之匹配的误差最小。

由此可见，最终用于后续可靠度分析的岩土参数最大熵法的概率密度函数，取决于直方图的分组。不同的直方图分组，与之匹配的最大熵法矩的阶数是不同的，概率密度函数形式不同，而哪一直方图分组取值更科学、合理，值得深入研究。

2.2 最优直方图分组

直方图分组常用的方法有经验法和基于渐进估计的方法,这些方法的主要缺点是存在一定的局限性，如经验法不针对任何的最优属性。将经验法中的Sturges规则与平均移位直方图的思想相结合，通过计算平均积分平方误差MISE的最小值可获得最优的直方图分组数。其主要思想是将估计最优直方图分组数的大小转化为估计MISE评价指标的最小值,MISE的评价指标为

式中，f⌢h(x)为样本直方图的概率密度估计；f(x)为样本真实的概率密度函数；E为期望。

ΔMISE(f⌢h)无法直接通过计算获得，只能通过样本数据估计，这样将估计直方图最佳分组的大小转化为估计ΔMISE(f⌢h)的最小值。

设R为样本的极差，m为直方图分组数，h=R/m为组距，k为落入[0,h]中样本的数量，样本在[0,h]的直方柱的预期高度θ为

θ的无偏估计为θ⌢=k/nh,n为样本的个数。则式(4)可表示为

进一步分解，有

设f(x)的平均值为θ¯,则式(7)可进一步分解为

采用无偏估计量E(θ⌢)=θ的分解规则，有

n个随机样本中有k个样本落入区间[0,h]的概率近似服从Poisson分布，则有

通过编制程序，输入实际岩土参数测试样本，根据式(11)求相邻Cn(h)的最小值，可得到最优的直方图组距h,从而获得最优的直方图分组数。

仍以容重和抗拉强度参数测试样本为例，进行直方图分组寻优之后，容重参数测试样本的最优直方图分组为6,抗拉强度参数测试样本的最优直方图分组亦为6。采用最大熵法绘制最优直方图分组下的概率密度曲线如图5和图6所示，检验结果如表1和表2所列。

为了验证本文方法的有效性，选择文献[1]中提供的黏聚力和内摩擦角测试样本(内摩擦角1#样本),共81 组，取直方图分组均为10,采用最优直方图分组寻优后，两组测试样本的最优直方图分组均为7。选择文献[14]中提供的两组内摩擦角测试样本，样本容量分别为42(内摩擦角2#样本)和63(内摩擦角3#样本),直方图分组分别取为13和16, 采用最优直方图分组寻优后，内摩擦角2#样本的最优直方图分组为6,内摩擦角3#样本的最优直方图分组为7。针对直方图不同分组，采用最大熵法拟合样本的概率分布，拟合效果检验值如表3—表6所列。限于篇幅，本文仅给出每组测试样本在矩的阶数为6,既有文献样本直方图分组取值及最优直方图分组下的拟合曲线如图7—图10所示。

图5 容重参数测试样本最大熵法拟合曲线

图6 抗拉强度参数测试样本最大熵法拟合曲线

图7 黏聚力参数测试样本拟合曲线

图8 内摩擦角参数测试1#样本拟合曲线

图9 内摩擦角参数测试2#样本拟合曲线

图10 内摩擦角参数测试3#样本拟合曲线

由图3—图6及表1—表6可知，无论是通过拟合曲线的直观对比，还是通过检验值的定量检验，均证明采用最优直方图分组后的最大熵法拟合精度均大于既有文献中直方图分组下的拟合精度。最优直方图分组方法具有严格的数学基础，不依赖外在的主观信息，是一种科学、合理的直方图分组方法。

3 最大熵法取值区间

3.1 取值区间的确定

积分界限R的确定很关键。通常取R∈[zmin,zmax],或者下界值取稍小于zmin,上界值取稍大于zmax。尽管这种方法简单，但对小容量样本是不合适的。原因是有限容量的样本与总体之间总会存在一定的差异，这样利用样本边界作为总体的界限亦必然有差异。样本的三阶矩可用来表征分布曲线的偏态特性，能充分考虑分布的不对称性。

文献[11,13]探讨了最大熵法积分区间的确定方法，文中对比分析了3σ、4σ、c33、c31、和c32型区间取值方法的优劣，最终确定c32型为最合理的积分区间。c32型区间取值方法为：以[μ-3σ,μ+3σ]为基础，参考偏度系数c进行调整，当c<0,左端边界取a=μ-(3-c)σ,减小下限值，右端b=μ+3σ;当c>0,右端边界取b=μ+(3+c)σ,增大上限值，左端a=μ-3σ。其中，μ为均值，σ为标准差，偏度系数c=m3/δ3,m3为第3阶中心矩。上述区间取值方法在实际中得到了有效应用，但当样本偏度系数较大，绝对值大于0.5时，上述方法并不有效。如3.2节中，标准风压参数测试样本，采用c32区间取值方法时出现左端截尾点为负的情况，这与实际岩土参数的分布不符；黏土液性指数参数样本出现右端波动，不应有的升高现象，说明上界选择过大。深入分析发现，c32型区间取值方法只考虑分布呈偏态时对样本某一界限的影响，而样本偏度系数较大时，忽略了样本另一界限的改变对分布曲线形状的影响。

在既有研究基础上，通过对实际测试样本的反复、统计计算，在c32型区间取值基础上加以补充，命名为c′31区间取值法，即以[μ-3σ,μ+3σ]为基础，参考偏度系数c进行调整：

当-0.5<c<0,左端边界取a=μ-(3-c)σ,减小下限值，右端b=μ+3σ。

当0<c<0.5,右端边界取b=μ+(3+c)σ,增大上限值，左端a=μ-3σ。

若c<-0.5,说明分布左偏，应减小下限，同时也适当减小上限，左端边界取a=μ-(3-c)σ,右端边界取b=μ+(3+c)σ。

若c>0.5,说明分布右偏，应增大上限，同时也适当增大下限，左端边界取a=μ-(3-c)σ,右端边界取b=μ+(3+c)σ。

3.2 两种方法拟合效果对比

以文献[13]中的标准风压参数样本为例，样本最大、最小值区间为[798,4 367],样本偏度系数为1.143 6。按照c32区间取值方法只增大上限，为[-623.091 7,5 642.286 7],左端取值为-623.091 7,显然这与岩土参数实际分布不符，文中做了截断取0处理。采用c31区间取值方法，在增大上限的同时，增大下限，为[405.76, 5 642.286 7],对比可见，c′31方法取值区间更接近样本区间。c′31区间取值拟合曲线如图11所示，c32区间取值拟合曲线如图12所示。

再以文献[1]中提供的内摩擦角测试样本为例，样本偏度系数为0.534 7,样本最大、最小值区间为[9.6,25.9],按照c32区间取值方法只增大上限，为[4.997,29.477 2]。按照c′31区间取值方法，在增大上限的同时，增大下限，为[6.995,29.477 2]。c′31区间取值拟合曲线如图13所示，c32区间取值拟合曲线如图14所示。由图14可见，按照c32区间取值方法，在区间[4.997,8.277]之间概率为0,过小的下界取值无意义。

以文献[11]提供的黏土容重参数样本为例，样本的偏度系数0.898 1,样本最大、最小值区间为[15.8,19.1],按照c32区间取值为[14.88,19.81]。按照c′31的方法区间取值为[15.31,19.81]。c′31区间取值拟合曲线如图15所示，c32区间取值拟合曲线如图16所示。由图16可见，按照c32区间取值方法，在区间[14.88,15.4]之间概率为0,过小的下界取值无意义。

再以文献[11]提供的黏土液性指数参数样本为例，样本的偏度系数-0.513 2,样本最大、最小值区间为[0.59,1.15],按照c32区间取值方法，减小下限值，为[0.513 5,1.370 1],按照c′31的区间取值方法，在减小下限的同时，减小上限，为[0.513 5,1.296 5]。c′31区间取值拟合曲线如图17所示，c32区间取值拟合曲线如图18所示。由图18可见，按照c32区间取值方法，在曲线的尾部出现了不应有的升高现象。

图11 标准风压样本c′31区间取值拟合曲线

图12 标准风压样本c32区间取值拟合曲线

图13 内摩擦角样本c′31区间取值拟合曲线

图14 内摩擦角样本c32区间取值拟合曲线

图15 黏土容重参数c′31区间取值拟合曲线

图16 黏土容重参数c32区间取值拟合曲线

图17 黏土液性指数c′31区间取值拟合曲线

图18 黏土液性指数c32区间取值拟合曲线

通过对测试样本的统计结果及图11—图18可见，c′31的区间取值方法，在样本偏度系数较大时，考虑了上、下界限的同时变化，弥补了c32区间取值法在样本偏度系数大于0.5时出现左端为负和右端波动的不足，最大熵分布不仅包络样本区间，且与样本的实际分布区间更接近。

4 结 论

本文在既有研究基础上，完善了最大熵法在岩土参数概率密度拟合上的应用，主要结论如下：

(1)直方图分组是最大熵法推断岩土参数概率分布的前提，基于平均积分平方误差的最小值获得最优直方图方法，具有严格的数学基础，不依赖外在的主观信息，是一种科学、合理的直方图分组方法。在最优的直方图分组下，最大熵法拟合岩土参数测试样本的概率分布误差最小。

(2)c′31的区间取值方法在样本偏度系数较大时，能同时考虑上、下界的变化，弥补了c32区间取值法在样本偏度系数绝对值大于0.5时出现左端为负和右端波动的不足，完善了最大熵法在岩土参数概率密度拟合上的应用。

[1] 朱唤珍, 李夕兵, 宫凤强. 大样本岩土参数概率分布的正态信息扩散推断[J]. 岩土力学, 2015, 36(11): 3275-3282.

ZHU Yuanzhen, LI Xibing, GONG Fengqiang. Inference of probability distribution of rock and soil parameters by normal information diffusion [J]. Rock and Soil Mechanics, 2015, 36(11): 3275-3282.

[3] 陈旺旺, 李典庆, 唐小松, 等. 岩土参数概率分布的最大熵估计及边坡可靠度分析 [ J]. 岩土力学, 2018, 39 ( 4): 1469- 1478.

CHEN Wangwang, LI Dianqing, TANG Xiaosong, et al. Probability distribution of shear strength parameters using maximum[ J]. Rock and Soil Mechanics, 2018, 39(4): 1469-1478.

[11] 黄天朗, 宫凤强, 伍婷玉. 小样本岩土参数概率分布的正交多项式推断方法[J]. 工程力学, 2017, 34(8): 85-95.

HUANG Tianlang, GONG Fengqiang, WU Tingyu. Orthogonal polynomial inference method of the probability distribution function for small sample of geotechnical parameters[J]. Engineering Mechanics, 2017, 34(8): 85-95.

[13] 宫凤强, 黄天朗, 李夕兵. 岩土参数最优概率分布推断方法及判别准则的研究[ J]. 岩石力学与工程学报, 2016, 35 ( 12): 2452-2460.

GONG Fengqiang, HUANG Tianlang, LI Xibing. Study on inference method and criterion of optimal probability distribution of geotechnical parameters [ J]. Chinese Journal of Rock Mechanics And Engineering, 2016, 35(12): 2452-2460.

[14] 宫凤强, 侯尚骞, 李夕兵. 岩土参数截尾分布的正态信息扩散推断方法[ J]. 武汉大学学报(工学版), 2016, 49( 5): 661- 667.

GONG Fengqiang, HOU Shangqian, LI Xibing. Truncated distribution deduction method for geotechnical parameters based on noramal information diffusion method [ J ]. Engineering Journal of Wuhan university, 2016, 49(5): 661-667.

[15] 邓建, 李夕兵, 古德生. 岩石力学参数概率分布的信息熵推断 [J]. 岩石力学与工程学报, 2004, 23(13): 2177-2181.

DENG Jian, LI Xibing, GU Desheng. Probability distribution of rock mechanics parameters by using maximum entropy method [J]. Chines Journal of Rock Mechanics and Engineering, 2004, 23(13): 2177- 2181.

[17] 李红英, 谭跃虎, 赵辉. 某滑坡体岩土参数概率分布统计分析方法研究[J]. 地下空间与工程学报, 2012, 8(3): 659-665.

LI Hongying, TAN Yuehu, ZHAO Hui. Research on the statistical analysis method of probability distribution of geotechnical parameters of a landslide[ J]. Journal of Underground Space and Engineering, 2012, 8(3): 659-665.

[18] 姚多喜, 鲁海峰, 邵亚红. 煤系岩体参数概率分布及Bayes 优 化[J]. 安徽理工大学学报(自然科学版), 2014, 34( 4): 15- 18.

YAO Duoxi, LU Haifeng, SHAO Yahong. Probability distribution of coal rock parameters and Bayes optimization [ J]. Journal of Anhui University of Technology (Natural Science Edition), 2014, 34(4): 15-18.

[21] 陈立周, 何晓峰, 翁海珊. 工程变量随即优化方法设计: 原理与应用[M]. 北京: 科学出版社, 1997.

CHEN Lizhou, HE Xiaofeng, WENG Haishan. Design of engineering variable optimization method-principle and application [M]. Beijing: Science Press, 1997.

水利水电技术（中英文）

水利部《水利水电技术（中英文）》杂志是中国水利水电行业的综合性技术期刊（月刊），为全国中文核心期刊，面向国内外公开发行。本刊以介绍我国水资源的开发、利用、治理、配置、节约和保护，以及水利水电工程的勘测、设计、施工、运行管理和科学研究等方面的技术经验为主，同时也报道国外的先进技术。期刊主要栏目有:水文水资源、水工建筑、工程施工、工程基础、水力学、机电技术、泥沙研究、水环境与水生态、运行管理、试验研究、工程地质、金属结构、水利经济、水利规划、防汛抗旱、建设管理、新能源、城市水利、农村水利、水土保持、水库移民、水利现代化、国际水利等。