小样本条件下砂砾石坝料级配特征参数的贝叶斯估计方法

刘 彪，赵宇飞，陈祖煜，王 毅，王文博

（1．中国水利水电科学研究院 流域水循环模拟与调控国家重点实验室，北京 100038；2．中国水利水电第八工程局有限公司，湖南 长沙 410007；3．中国水利水电第六工程局有限公司，辽宁 沈阳 110179）

1 研究背景

砂砾石料作为一种无黏性的粗粒土材料，因其压实后具有较高的强度和变形模量，且对不同的地基条件具有良好的适应性等优良工程特性而在水利工程中得到了广泛的应用。新疆南疆地区地域辽阔，河床覆盖层厚，天然砂砾石料广泛分布于河床和戈壁沙漠，储量丰富。当地砂砾石坝因筑坝材料就地取材、结构简单、变形适应性好、施工工序少等优点，成为该地区水利水电开发过程中普遍选用的坝型［1］。在砂砾石坝的填筑施工过程中，准确把握施工过程中碾压参数和筑坝材料的级配特征参数是保证大坝填筑施工质量的关键。近年来随着大坝智能化建设水平的提升，大坝碾压智能监控系统实现了对碾压施工参数的全过程、精细化控制，为我们准确把控施工过程中的碾压参数提供了便利［2-3］。而对于筑坝材料的级配特征参数而言，由于坝料的级配特性受母岩特性、产生原因、空间分布等不确定因素的影响，当前采用试验方法仍然是确定其参数最直接、最可靠的方法，也是大多数理论和经验估值的重要依据［4］。已有研究表明，筑坝材料的压实质量（用干密度作为控制指标）具有强烈的级配相关性，受P5含量、最大粒径、曲率系数等级配参数的影响尤为显著［5］。准确把握填筑单元坝料的级配特征参数不仅有助于研究筑坝材料物理力学特性（渗透系数、强度参数及变形特性参数）［6］，而且有助于确定筑坝填筑标准及施工质量控制。

在实际工程中，筑坝材料的级配特征参数常采用筛分法进行确定，这种试验方法存在采样率低，操作过程繁琐，检测过程与结果受人为因素影响大等缺点以致检测结果代表性差，对于某一填筑单元而言，挖坑检测数目极少，仅3到5组，这使得现有的大坝碾压智能监控系统中所采用的压实质量评估模型不能有效地考虑级配特征参数对压实质量的影响。目前很多学者将连续压实指标与多种影响压实质量的碾压参数相耦合，采用回归方法或人工神经网络方法构建压实质量评估模型，但是评估模型中并没有充分考虑含水量、级配等筑坝材料属性的变化，这意味着在不同的筑坝材料物理力学特性下，相同的碾压施工指标值并不一定代表相同的密实度［7-8］。当坝料级配发生变化时，已有模型的评估精度会受到一定的影响，从而造成结果的误判。因此，如何根据已有的小样本数据去推断构建坝料级配的总体分布，从而准确把握坝料的物理力学特性，实现级配特征参数与海量施工数据深度融合，搭建施工过程大坝填筑压实质量评估模型来合理评估坝料碾压密实度，是当前大坝填筑碾压监控系统建设中需着重解决的难点问题。

现阶段根据小样本数据（样本容量n≤30）［9］去推断总体分布的方法主要有贝叶斯方法、基于计算机仿真的Bootstrap 方法、蒙特卡洛方法以及历史数据融合等方法［10］，在这些方法中贝叶斯方法充分利用历史试验数据作为先验信息，并结合最新的现场试验数据对先验信息进行动态修正以做出较准确的总体分布预测［11-13］，近年来在岩土工程领域得到了广泛的发展与应用。徐军等［14］将模糊综合评判法和贝叶斯理论相结合，探讨了由小样本数据确定岩土参数的概率分布；吴越等［15］提出了岩土强度参数的随机分布特征参数均值和方差的联合先验分布为正态—逆伽马分布，并根据贝叶斯公式推导了相应的共轭后验分布及最大后验估计量的计算公式；王俊杰等［16］利用卡方检验法对样本概率分布进行拟合，采用贝叶斯方法对分布参数进行优化，为岩土工程参数优化提供了新的途径；赵宇飞等［17］假定岩土强度参数服从二维正态分布的基础上，利用Bayes方法计算得到了后验分布密度函数中各参数的计算公式，并通过锦屏一级水电站中硬性结构面抗剪强度参数优化的实例验证了方法的可行性。从上述研究可以发现贝叶斯方法借助先验信息，有效降低了对评估样本数据量的需求，在小样本数据推断分析中具有明显的优势，现阶段贝叶斯理论在处理正态分布以及指数分布等简单分布的应用中已经比较成熟，然而对于复杂的分布如威布尔分布，由于没有共轭先验分布，先验分布的获取极为困难，贝叶斯推导繁杂且大概率无解［18］。

为解决上述问题，本文首先采用参数化Bootstrap 方法扩大数据样本，其次利用非参数核密度估计法直接从抽样结果的样本特性出发拟合概率密度函数，确定小样本数据下待估参数的先验分布；然后利用贝叶斯理论结合现场试验数据对先验分布加以修正得出参数所服从的后验分布；最后采用混合Gibbs抽样方法对后验分布进行模拟仿真求解，给出了基于贝叶斯理论的两参数后验威布尔分布估计结果。本文方法简化了两参数威布尔分布直接运用贝叶斯公式时存在解析解难以计算的问题，为小样本数据威布尔分布的分析评估提供了新的思路。

本文研究中，依托新疆大石门水利枢纽工程开展了基于贝叶斯理论的砂砾石料级配特征参数概率分布估计的系统研究，所涉及的砂砾石料颗粒浑圆坚硬，与另外一种常见堆石坝料——爆破料相比，在碾压施工过程中砂砾石颗粒不易破碎，压实后具有较高的强度和变形模量，碾压施工后颗粒级配曲线与料场原始的级配曲线较为接近［19］，因此通过本文提出的小样本条件下坝料级配特征参数的贝叶斯估计方法可以动态的获取某碾压单元中坝料级配特征参数的总体分布，为大坝填筑施工过程中实时准确评估砂砾石坝料的压实特性提供了重要的数据支撑。

2 研究框架及方法介绍

本文的研究框架如图1所示，具体的求解步骤如下：

图1 研究框架

（1）以大石门水利枢纽工程砂砾石料碾压质量挖坑检测得到的级配参数小样本数据为研究对象，通过拟合优度检验确定威布尔分布作为总体分布模型，求解待估参数θ；

（2）采用参数化Bootstrap 法与非参数核密度估计法来对待估参数进行概率描述获得贝叶斯先验分布g（θ）；

（3）结合先验分布和现场试验数据，根据贝叶斯公式求得后验分布g（θ|x），为求解方便采用混合Gibbs抽样的方法对后验分布进行模拟仿真求解；

（4）通过后验分布的特征值确定在一定可靠度条件下各级配参数的估计值，此外根据贝叶斯方法更新后的威布尔分布作为总体分布，产生随机数作为压实质量评估模型中坝料级配特征参数的输入指标。

2.1 贝叶斯理论 贝叶斯理论充分利用了先验信息和现场试验信息，通过现场试验信息对先验信息进行修正得到更为准确的后验信息。贝叶斯公式一般表示为：

式中：θ为总体分布中的未知参数，即所要估计的威布尔分布中的形状参数ξ和尺度参数η；g（θ）为参数θ的先验分布密度函数；为给定θ时xi的条件概率分布，常称为似然分布；g（θ|x）为给定x时θ的条件分布即后验分布；不依赖于θ，在计算θ的后验分布中仅起到一个正则化的作用［20］，则公式（1）可以简化为：

式中：∝表示两端仅差一个常数因子；f（x|θ）g（θ）是后验分布g（θ|x）的核。

2.2 先验分布的计算 根据已有的小样本数据去推断总体分布时，一般需要将数据拟合成具体的概率分布模型，并估计其相应的分布参数，在数理统计分析中威布尔分布对于各类型试验数据有极强的适用性，可以得到小样本条件下级配特征参数分布较为准确的估计，是目前常用的概率统计方法。在采用贝叶斯理论进行小样本数据分析时，选取准确合适的先验分布是贝叶斯方法应用的关键，不同形式的先验分布产生的贝叶斯后验评估结果差异较大。目前常用的先验分布确定方法主要有无信息先验分布和共轭先验分布。威布尔分布作为一种较为复杂的经验分布函数，直接应用于贝叶斯分析时缺少相应的共轭先验分布，这使得先验分布的确定较为困难［10］。鉴于此本文采用最小二乘法初步拟合两参数威布尔分布并进行相应的拟合优度检验，而后采用参数化的Bootstrap 方法增加样本信息量，最后采用非参数核密度估计法从抽样数据特性出发拟合出对应的概率密度函数，计算出待估参数的先验分布。先验分布的具体求解过程如下：

（1）威布尔分布参数求解及拟合优度检验。两参数威布尔概率分布函数的表达式可表示为：

其密度函数为：

式中：F（x）为分布函数；f（x）为密度函数；x为随机变量；ξ为形状参数；η为尺度参数。

本文首先采用最小二乘法［21-22］对威布尔分布参数进行初步求解，求得待估参数后，利用相关系数公式（5）评估最小二乘法拟合直线的效果，该值的绝对值越接近于1说明拟合效果越好。然后为了使初步得到的威布尔分布模型在一定程度上具有鲁棒性，本文采用K-S检验［23］对其进行拟合优度检验。

（2）参数化Bootstrap 法增加样本信息量。采用最小二乘法初步拟合得到小样本威布尔分布后，为了提高统计分析的精度，常需要设法增加样本信息量，目前工程中常用的方法是美国斯坦福大学B．Efron教授提出的Bootstrap 法［24］，其基本思想是依据小样本信息来模拟未知分布，通过不断获取再生样本从而实现小样本转换为大样本。传统的非参数Bootstrap 方法在建立先验分布时存在抽样误差大的问题，为了减少此误差，使经验分布函数的泛化能力更强，本文以威布尔分布作为先验样本的经验分布Fn，然后再采用Bootstrap 法从Fn中抽取再生样本，重复抽取n次，对每组样本分别使用最小二乘法进行威布尔分布拟合得到n组形状参数ξ和尺度参数η。

（3）非参数核密度估计确定先验分布。为了减少抽样带来的误差，提高经验分布函数的泛化能力，采用非参数核密度估计法从抽样数据特性出发拟合出对应的概率密度函数，计算出待估参数的先验分布。核密度估计公式为［25］：

式中：K（•）为核函数，当样本量比较大时，核函数的选择对估计结果的影响不大，本文选取最为常用的高斯核函数作为核函数；h为窗宽，最佳窗宽的估计值为其中σ为样本标准差，n为样本数。

将Bootstrap 方法抽样得到的n组形状参数ξ和尺度参数η分别代入式（6）即可得到待估参数θ（形状参数ξ或尺度参数η）的先验分布公式为：

本文先验分布确定方法有效降低了对评估样本数据量的需求，同时也避免了因威布尔分布参数多而导致的计算复杂的问题。

2.3 后验分布的计算 在获得先验分布后，结合现场试验数据X=（x1，x2，…，xn），可得到后验分布公式为：

根据以上后验分布函数来计算待估参数的估计值时需要进行数值积分来逼近待估参数的期望和方差，本文后验分布是非标准形式的密度函数，采用数值积分方法计算时会有高维积分出现，使得逼近误差随着维数的增加而增加，同时也增加了计算难度［26］。马尔科夫链蒙特卡罗（Markov Chain Monte-Carlo，MC-MC）［27-29］方法是解决此类问题的一种行之有效的方法，基本思想是通过建立平稳分布为g（θ|x）的Markov 链来得到后验分布g（θ|x）的样本，随数值仿真模拟的变化而实时抽取随机样本，来动态模拟求取积分，基于这些样本可做各种统计推断。在贝叶斯方法中常见的构造马尔科夫链的方法是Gibbs 抽样，其常用来处理高维、非标准形式的后验分布。在抽样过程中后验分布g（ξ，η|x）若给定η，则g（ξ，η|x）仅为ξ的函数，此时称g（ξ|x，η）为参数ξ的满条件分布。在Gibbs 抽样过程中，从没有显式形式的满条件分布中抽样比较困难，故本文引入Metropolis 算法抽取随机数，将Gibbs 抽样和Metropolis 算法［30］结合起来采用混合Gibbs 抽样的方法来构造马尔科夫链。其中Gibbs 抽样方法［31］研究已较为成熟，本文不再赘述，Metropolis算法［32-33］的实现方法如下：

（1）初始化：t=0，选择一个对称提议分布q（x），满足q（x|y）=q（y|x），给定一个起始样本点xt，迭代终止值为T；

（2）令t=t+1，从q（x|xt-1）中生成候选样本x′；

（3）计算接收概率α，α=min{1，x′/xt-1}；

（4）从均匀分布中抽取随机αt，若αt≤α，则接受候选样本，xt=x′；否则，拒绝候选样本，并令xt=xt-1；

（5）重复（2）—（4）步直到迭代终止。

采用混合Gibbs抽样同时产生多条Markov链，若这几条链稳定下来则认为混合Gibbs抽样收敛。

2.4 可靠度评价 在实际工程中，往往需要给出在一定可靠度条件下参数的估计值，设参数估计的可靠度为R，分布函数与可靠度的关系为：

对上式取对数可得：

式中xR为两参数威布尔分布在可靠度R条件下的估计值。将基于贝叶斯理论更新得到的两参数后验威布尔分布估计结果代入式（10），即可计算出在可靠度R条件下的砂砾料级配特征参数的估计值xR。此外根据贝叶斯方法更新后的威布尔分布作为总体分布，产生的随机数可作为压实质量评估模型中坝料级配特征参数的输入指标。

3 算例验证

假设砂砾石料P5 含量服从形状参数ξ=3、尺度参数η=2 的威布尔分布，从中产生25 个随机数作为先验样本数据集，5 个随机数作为现场试验数据集，采用本文方法编制的程序对上述小样本数据进行威布尔分布的贝叶斯估计求解。在对后验分布进行模拟仿真求解时，本文构建了迭代次数为10 000、20 000、30 000 的三条Markov 链，仿真过程中监控形状参数ξ和尺度参数η输出的模拟估计值，经过统计分析用于判断迭代过程是否收敛。三条马尔科夫链输出结果的统计分析结果如表1 所示。从表1 和图2 可以看出，随着迭代次数的增加，MC 误差在逐步的减小，且三条链在2．5%、50%和97．5%这三个分位数对应的估值基本相同，这表明在迭代过程中Markov 链逐渐收敛并趋于稳定。当Markov 链达到平稳分布时，将迭代计算结果的均值作为后验形状参数和尺度参数的估计值，后验形状参数ξ=3．1950 和尺度参数η=2．0380 与真实值的相对误差分别为6．5%和1．9%，在误差允许范围内，由此验证了本文小样本条件下贝叶斯估计方法和编制程序的准确性。

图2 待估参数的后验分位数

表1 后验分布参数的模拟估计值

上述通过一次蒙特卡罗模拟验证了本文方法的有效性，为了消除随机性的影响，先后进行了100次蒙特卡罗模拟，并分别采用了本文方法与传统Bootstrap 方法进行贝叶斯估计对比分析。后验分布时通过构建迭代次数为10 000 的Markov 链进行模拟仿真，计算得到的贝叶斯估计均值、与真实值的相对误差和均方误差如表2所示。从表2中可以看出，采用本文方法计算得到的贝叶斯估计均值、相对误差和均方误差都比采用传统Bootstrap 方法计算得到的值小，说明本文方法相比传统Bootstrap 方法具有较高的精度，这在一定程度上也说明了本文算法的有效性及稳定性。

表2 不同方法对比分析

4 工程实例验证与应用

4.1 项目背景及试验数据 大石门水利枢纽工程位于新疆维吾尔自治区巴音郭楞蒙古族自治州且末县车尔臣河干流之上，坝址位于车尔臣河出山口、与支流托其里萨依河汇合口下游300 m处。工程是车尔臣河流域规划中确定的近期重点开发控制性枢纽工程，属国家节水供水172 项重大水利工程之一，是一项承担防洪、发电和灌溉等任务的综合性水利枢纽工程。为了有效指导施工，保证大坝砂砾料回填质量，相关单位对新疆车尔臣河大石门水利枢纽工程大坝填筑的砂砾料进行现场原型级配碾压试验，共获得23组小样本级配数据（级配曲线如图3），此外还有4组为现场某一单元工程确定的挖坑试验数据（级配曲线如图4），根据级配曲线计算得到的P5 含量、最大粒径以及曲率系数数据详见表3。

表3 试验数据

图3 先验小样本颗粒级配曲线图

图4 现场试验颗粒级配曲线图

4.2 实例应用 首先采用最小二乘法对样本容量为23的先验小样本级配参数数据进行两参数威布尔分布拟合，参数的拟合结果如图5所示，对应拟合分布的统计参数详见表4。从表中可以看出各级配参数的相关系数均大于0．9，拟合结果较为可观。在显著性水平α=0．05条件下，计算得到的K-S统计量均小于临界值D（23，0．05）=0．624，说明本研究所选用的威布尔分布模型能够较好的适用于砂砾石坝料级配参数的描述及参数估计。

图5 砂砾石料级配参数威布尔拟合

表4 拟合分布的统计参数

接着利用Bootstrap抽样从上述级配参数对应的Fn（x）中抽取1000组再生子样本，每组样本的数量为23个，对每组样本分别使用最小二乘法进行威布尔分布拟合得到1000组形状参数ξ和尺度参数η。在对得到的抽样结果进行统计分析时，为了克服结果对分布类型过分依赖而造成主观性过大的问题，本文采用了非参数核密度估计法直接从抽样结果的样本特性出发拟合概率密度函数，计算得到1000组形状参数和尺度参数的样本标准差和最优窗宽，进而得到待估参数的先验密度函数，计算结果见表5。形状参数和尺度参数对应的概率直方图和密度函数图如图6所示。

图6 核密度估计拟合

表5 核密度估计参数及密度函数

在求得先验分布后，结合4组现场试验数据对先验分布进行修正，代入式（8）可以得到贝叶斯后验分布表达式为：

针对上述非标准形式的Bayes后验分布密度函数，为求解方便本文依据混合Gibbs抽样原理编制了相应的Python程序对后验分布进行抽样，产生具有初始值的Markov链，经过一段时间的迭代，当Markov链收敛到平稳分布时即可求得待估参数的后验分布样本。为了判断迭代过程是否收敛，分别构建了迭代次数为10 000、20 000、30 000的三条Markov链，通过监控后验形状参数ξ和尺度参数η的自相关函数图和抽样迭代轨迹进行收敛判断分析。限于篇幅，图7仅列出P5含量的模拟结果。图7（a）中自相关函数图表示当前迭代结果与之前迭代结果的相关性，用以判断收敛速度的快慢，从图中可以发现在Lag=10左右样本的自相关系数已接近于零，表明收敛速度很快。为提高模拟效率，本文每隔10个数取一个样本进行统计分析，待估参数的统计模拟估计值见表6。从表中可以看出，随着迭代次数的增加形状参数ξ和尺度参数η的均值和标准差均在小范围内波动，MC误差在逐步的减小，这表明在迭代过程中Markov链逐渐收敛并趋于稳定。图7（b）为形状参数ξ和尺度参数η在迭代30 000次后的动态轨迹，从迭代轨迹图中同样可以发现Markov链基本收敛在一定区域内，这时的抽样结果可以认为是后验分布的样本，其对应的概率密度函数见图7（c）。上述当Markov链达到平稳分布时，将迭代计算结果的均值作为后验形状参数ξ和尺度参数η的估计值，可以得到砂砾石坝料级配参数经贝叶斯方法更新后的威布尔分布概率分布函数以及在可靠度R条件下的坝料级配特征参数的估计值xR，详见表7。

表7 贝叶斯方法更新后的威布尔分布概率分布函数

图7 P5含量模拟结果

表6 后验分布参数的模拟估计值

4.3 有效性评价 为了验证本文贝叶斯方法更新后的威布尔分布作为总体分布的有效性，本文随机抽取了大石门水利工程不同填筑单元的十组现场挖坑检测级配数据X（如表8所示），运用K-S检验方法验证现场数据X是否服从本文贝叶斯方法更新后的威布尔分布。首先建立原假设H0：F0（x）=Fn（x），其中Fn（x）为本文贝叶斯方法更新后的威布尔分布函数，F0（x）为样本观测值的累积分布函数：F0（xi）=i/n，i=1，2，…，n。然后计算K-S 检验统计量：Dn=max{|Fn（x）-F0（x）|}，计算结果如表9 所示，从表中可以看出在在显著性水平α=0．05条件下，坝料各级配特征参数计算得到的K-S统计量均小于临界值D（10，0．05）=0．410，说明随机抽取的大石门工程不同填筑单元的级配数据服从本文贝叶斯方法更新后的威布尔分布函数，验证了本文贝叶斯方法更新后的威布尔分布作为总体分布的有效性。

表8 大石门工程随机抽取的现场级配特征参数数据

表9 K-S检验计算结果

在实际工程中，大坝填筑施工之前都会进行室内外试验，确定坝料合理的最大最小干密度、坝料级配特征参数，并且进一步确定大坝填筑碾压施工参数与控制指标。但是在实际施工过程中常以有限的已碾压单元挖坑试验结果进行级配特征参数描述，无法获得整个工作仓级配特征参数的总体分布规律，也无法为施工过程中坝料碾压质量的评价提供数据支撑。通过本文提出的小样本条件下坝料级配特征参数的贝叶斯估计方法可以动态的获取某填筑单元在一定可靠度条件下级配特征参数的估计值，以此作为评价某单元工程级配特征参数的合理依据。此外根据本文贝叶斯方法更新后的威布尔分布作为总体分布，产生随机数作为压实质量评估模型中坝料级配特征参数的输入指标用于坝料压实质量特性研究。

5 结论与展望

针对砂砾石坝，坝料的级配特征参数是确定坝料压实特性、强度及变形参数所需要的最重要的指标之一，对于同一来源的砂砾石料，其颗粒组成具有明显的规律性，服从一定的概率分布规律。为实现通过已有的砂砾石坝挖坑检测小样本级配参数数据来推断构建总体分布，本文以大石门水利枢纽沥青混凝土心墙砂砾石坝的坝壳料碾压质量挖坑检测得到的小样本级配数据为研究对象，根据贝叶斯理论将先验小样本数据和现场试验数据进行融合，最终求解出了砂砾石料级配参数的后验分布函数。得出以下结论：

（1）由于双参数威布尔分布不存在共轭先验分布，本文采用了参数化Bootstrap 方法和非参数核密度估计法确定了小样本数据下威布尔分布参数的先验分布，有效地解决了威布尔分布下贝叶斯公式难以求解的问题；针对后验分布积分计算存在高维积分，计算量大的问题，本文采用了混合Gibbs抽样方法对后验分布进行模拟仿真求解，模拟精度较高，算例验证了本文方法的有效性与稳定性。

（2）通过算例对本文方法与传统Bootstrap 方法进行贝叶斯评估进行了比较，结果表明本文方法精度更高，计算结果更为稳定，说明了本文提出方法的优越性与准确性。

（3）本研究通过随机抽取大石门工程不同填筑单元的级配数据验证了本文贝叶斯方法更新后的威布尔分布作为坝料级配特征参数总体分布的有效性，通过本文提出的小样本条件下坝料级配特征参数的贝叶斯估计方法可以动态的获取某填筑单元工程在一定可靠度条件下级配特征参数的估计值，以此作为评价某单元工程级配特征参数的合理依据。此外，本文贝叶斯方法更新后的威布尔分布可作为坝料级配特征参数的总体分布，这对于研究大坝填筑碾压施工质量控制、压实后坝体材料力学参数估计及坝体变形分析与预测都具有重要的数据支撑作用。

土石坝筑坝材料的级配特征参数分布规律是一个具有很强实践性的问题，因此本文贝叶斯方法将不断把现场挖坑检测获得的真实级配数据融合到先验分布中，使获得的后验威布尔分布更具有实践性和真实性。需要注意的是，级配的表征其实是坝料级配特征参数（P5含量，曲率系数和最大粒径）之间相关联的联合分布，但是由于计算的难度以及为了满足实际工程切实的需求，本文将不同级配特征参数进行了单独的统计分析，这也为后续深入开展级配特征参数间的联合分布奠定了基础。此外本文仅根据大石门水利枢纽工程现场挖坑检测的砂砾石料级配数据进行了统计分析，后续将收集更多水利工程现场碾压级配数据以建立不同地质条件下坝料的先验分布，由此获取更准确的后验分布估计。