小樣本條件下砂礫石壩料級配特征參數(shù)的貝葉斯估計方法

劉 彪，趙宇飛，陳祖煜，王 毅，王文博

（1．中國水利水電科學研究院 流域水循環(huán)模擬與調(diào)控國家重點實驗室，北京 100038；2．中國水利水電第八工程局有限公司，湖南 長沙 410007；3．中國水利水電第六工程局有限公司，遼寧 沈陽 110179）

1 研究背景

砂礫石料作為一種無黏性的粗粒土材料，因其壓實后具有較高的強度和變形模量，且對不同的地基條件具有良好的適應性等優(yōu)良工程特性而在水利工程中得到了廣泛的應用。新疆南疆地區(qū)地域遼闊，河床覆蓋層厚，天然砂礫石料廣泛分布于河床和戈壁沙漠，儲量豐富。當?shù)厣暗[石壩因筑壩材料就地取材、結構簡單、變形適應性好、施工工序少等優(yōu)點，成為該地區(qū)水利水電開發(fā)過程中普遍選用的壩型［1］。在砂礫石壩的填筑施工過程中，準確把握施工過程中碾壓參數(shù)和筑壩材料的級配特征參數(shù)是保證大壩填筑施工質(zhì)量的關鍵。近年來隨著大壩智能化建設水平的提升，大壩碾壓智能監(jiān)控系統(tǒng)實現(xiàn)了對碾壓施工參數(shù)的全過程、精細化控制，為我們準確把控施工過程中的碾壓參數(shù)提供了便利［2-3］。而對于筑壩材料的級配特征參數(shù)而言，由于壩料的級配特性受母巖特性、產(chǎn)生原因、空間分布等不確定因素的影響，當前采用試驗方法仍然是確定其參數(shù)最直接、最可靠的方法，也是大多數(shù)理論和經(jīng)驗估值的重要依據(jù)［4］。已有研究表明，筑壩材料的壓實質(zhì)量（用干密度作為控制指標）具有強烈的級配相關性，受P5含量、最大粒徑、曲率系數(shù)等級配參數(shù)的影響尤為顯著［5］。準確把握填筑單元壩料的級配特征參數(shù)不僅有助于研究筑壩材料物理力學特性（滲透系數(shù)、強度參數(shù)及變形特性參數(shù)）［6］，而且有助于確定筑壩填筑標準及施工質(zhì)量控制。

在實際工程中，筑壩材料的級配特征參數(shù)常采用篩分法進行確定，這種試驗方法存在采樣率低，操作過程繁瑣，檢測過程與結果受人為因素影響大等缺點以致檢測結果代表性差，對于某一填筑單元而言，挖坑檢測數(shù)目極少，僅3到5組，這使得現(xiàn)有的大壩碾壓智能監(jiān)控系統(tǒng)中所采用的壓實質(zhì)量評估模型不能有效地考慮級配特征參數(shù)對壓實質(zhì)量的影響。目前很多學者將連續(xù)壓實指標與多種影響壓實質(zhì)量的碾壓參數(shù)相耦合，采用回歸方法或人工神經(jīng)網(wǎng)絡方法構建壓實質(zhì)量評估模型，但是評估模型中并沒有充分考慮含水量、級配等筑壩材料屬性的變化，這意味著在不同的筑壩材料物理力學特性下，相同的碾壓施工指標值并不一定代表相同的密實度［7-8］。當壩料級配發(fā)生變化時，已有模型的評估精度會受到一定的影響，從而造成結果的誤判。因此，如何根據(jù)已有的小樣本數(shù)據(jù)去推斷構建壩料級配的總體分布，從而準確把握壩料的物理力學特性，實現(xiàn)級配特征參數(shù)與海量施工數(shù)據(jù)深度融合，搭建施工過程大壩填筑壓實質(zhì)量評估模型來合理評估壩料碾壓密實度，是當前大壩填筑碾壓監(jiān)控系統(tǒng)建設中需著重解決的難點問題。

現(xiàn)階段根據(jù)小樣本數(shù)據(jù)（樣本容量n≤30）［9］去推斷總體分布的方法主要有貝葉斯方法、基于計算機仿真的Bootstrap 方法、蒙特卡洛方法以及歷史數(shù)據(jù)融合等方法［10］，在這些方法中貝葉斯方法充分利用歷史試驗數(shù)據(jù)作為先驗信息，并結合最新的現(xiàn)場試驗數(shù)據(jù)對先驗信息進行動態(tài)修正以做出較準確的總體分布預測［11-13］，近年來在巖土工程領域得到了廣泛的發(fā)展與應用。徐軍等［14］將模糊綜合評判法和貝葉斯理論相結合，探討了由小樣本數(shù)據(jù)確定巖土參數(shù)的概率分布；吳越等［15］提出了巖土強度參數(shù)的隨機分布特征參數(shù)均值和方差的聯(lián)合先驗分布為正態(tài)—逆伽馬分布，并根據(jù)貝葉斯公式推導了相應的共軛后驗分布及最大后驗估計量的計算公式；王俊杰等［16］利用卡方檢驗法對樣本概率分布進行擬合，采用貝葉斯方法對分布參數(shù)進行優(yōu)化，為巖土工程參數(shù)優(yōu)化提供了新的途徑；趙宇飛等［17］假定巖土強度參數(shù)服從二維正態(tài)分布的基礎上，利用Bayes方法計算得到了后驗分布密度函數(shù)中各參數(shù)的計算公式，并通過錦屏一級水電站中硬性結構面抗剪強度參數(shù)優(yōu)化的實例驗證了方法的可行性。從上述研究可以發(fā)現(xiàn)貝葉斯方法借助先驗信息，有效降低了對評估樣本數(shù)據(jù)量的需求，在小樣本數(shù)據(jù)推斷分析中具有明顯的優(yōu)勢，現(xiàn)階段貝葉斯理論在處理正態(tài)分布以及指數(shù)分布等簡單分布的應用中已經(jīng)比較成熟，然而對于復雜的分布如威布爾分布，由于沒有共軛先驗分布，先驗分布的獲取極為困難，貝葉斯推導繁雜且大概率無解［18］。

為解決上述問題，本文首先采用參數(shù)化Bootstrap 方法擴大數(shù)據(jù)樣本，其次利用非參數(shù)核密度估計法直接從抽樣結果的樣本特性出發(fā)擬合概率密度函數(shù)，確定小樣本數(shù)據(jù)下待估參數(shù)的先驗分布；然后利用貝葉斯理論結合現(xiàn)場試驗數(shù)據(jù)對先驗分布加以修正得出參數(shù)所服從的后驗分布；最后采用混合Gibbs抽樣方法對后驗分布進行模擬仿真求解，給出了基于貝葉斯理論的兩參數(shù)后驗威布爾分布估計結果。本文方法簡化了兩參數(shù)威布爾分布直接運用貝葉斯公式時存在解析解難以計算的問題，為小樣本數(shù)據(jù)威布爾分布的分析評估提供了新的思路。

本文研究中，依托新疆大石門水利樞紐工程開展了基于貝葉斯理論的砂礫石料級配特征參數(shù)概率分布估計的系統(tǒng)研究，所涉及的砂礫石料顆粒渾圓堅硬，與另外一種常見堆石壩料——爆破料相比，在碾壓施工過程中砂礫石顆粒不易破碎，壓實后具有較高的強度和變形模量，碾壓施工后顆粒級配曲線與料場原始的級配曲線較為接近［19］，因此通過本文提出的小樣本條件下壩料級配特征參數(shù)的貝葉斯估計方法可以動態(tài)的獲取某碾壓單元中壩料級配特征參數(shù)的總體分布，為大壩填筑施工過程中實時準確評估砂礫石壩料的壓實特性提供了重要的數(shù)據(jù)支撐。

2 研究框架及方法介紹

本文的研究框架如圖1所示，具體的求解步驟如下：

圖1 研究框架

（1）以大石門水利樞紐工程砂礫石料碾壓質(zhì)量挖坑檢測得到的級配參數(shù)小樣本數(shù)據(jù)為研究對象，通過擬合優(yōu)度檢驗確定威布爾分布作為總體分布模型，求解待估參數(shù)θ；

（2）采用參數(shù)化Bootstrap 法與非參數(shù)核密度估計法來對待估參數(shù)進行概率描述獲得貝葉斯先驗分布g（θ）；

（3）結合先驗分布和現(xiàn)場試驗數(shù)據(jù)，根據(jù)貝葉斯公式求得后驗分布g（θ|x），為求解方便采用混合Gibbs抽樣的方法對后驗分布進行模擬仿真求解；

（4）通過后驗分布的特征值確定在一定可靠度條件下各級配參數(shù)的估計值，此外根據(jù)貝葉斯方法更新后的威布爾分布作為總體分布，產(chǎn)生隨機數(shù)作為壓實質(zhì)量評估模型中壩料級配特征參數(shù)的輸入指標。

2.1 貝葉斯理論 貝葉斯理論充分利用了先驗信息和現(xiàn)場試驗信息，通過現(xiàn)場試驗信息對先驗信息進行修正得到更為準確的后驗信息。貝葉斯公式一般表示為：

式中：θ為總體分布中的未知參數(shù)，即所要估計的威布爾分布中的形狀參數(shù)ξ和尺度參數(shù)η；g（θ）為參數(shù)θ的先驗分布密度函數(shù)；為給定θ時xi的條件概率分布，常稱為似然分布；g（θ|x）為給定x時θ的條件分布即后驗分布；不依賴于θ，在計算θ的后驗分布中僅起到一個正則化的作用［20］，則公式（1）可以簡化為：

式中：∝表示兩端僅差一個常數(shù)因子；f（x|θ）g（θ）是后驗分布g（θ|x）的核。

2.2 先驗分布的計算 根據(jù)已有的小樣本數(shù)據(jù)去推斷總體分布時，一般需要將數(shù)據(jù)擬合成具體的概率分布模型，并估計其相應的分布參數(shù)，在數(shù)理統(tǒng)計分析中威布爾分布對于各類型試驗數(shù)據(jù)有極強的適用性，可以得到小樣本條件下級配特征參數(shù)分布較為準確的估計，是目前常用的概率統(tǒng)計方法。在采用貝葉斯理論進行小樣本數(shù)據(jù)分析時，選取準確合適的先驗分布是貝葉斯方法應用的關鍵，不同形式的先驗分布產(chǎn)生的貝葉斯后驗評估結果差異較大。目前常用的先驗分布確定方法主要有無信息先驗分布和共軛先驗分布。威布爾分布作為一種較為復雜的經(jīng)驗分布函數(shù)，直接應用于貝葉斯分析時缺少相應的共軛先驗分布，這使得先驗分布的確定較為困難［10］。鑒于此本文采用最小二乘法初步擬合兩參數(shù)威布爾分布并進行相應的擬合優(yōu)度檢驗，而后采用參數(shù)化的Bootstrap 方法增加樣本信息量，最后采用非參數(shù)核密度估計法從抽樣數(shù)據(jù)特性出發(fā)擬合出對應的概率密度函數(shù)，計算出待估參數(shù)的先驗分布。先驗分布的具體求解過程如下：

（1）威布爾分布參數(shù)求解及擬合優(yōu)度檢驗。兩參數(shù)威布爾概率分布函數(shù)的表達式可表示為：

其密度函數(shù)為：

式中：F（x）為分布函數(shù)；f（x）為密度函數(shù)；x為隨機變量；ξ為形狀參數(shù)；η為尺度參數(shù)。

本文首先采用最小二乘法［21-22］對威布爾分布參數(shù)進行初步求解，求得待估參數(shù)后，利用相關系數(shù)公式（5）評估最小二乘法擬合直線的效果，該值的絕對值越接近于1說明擬合效果越好。然后為了使初步得到的威布爾分布模型在一定程度上具有魯棒性，本文采用K-S檢驗［23］對其進行擬合優(yōu)度檢驗。

（2）參數(shù)化Bootstrap 法增加樣本信息量。采用最小二乘法初步擬合得到小樣本威布爾分布后，為了提高統(tǒng)計分析的精度，常需要設法增加樣本信息量，目前工程中常用的方法是美國斯坦福大學B．Efron教授提出的Bootstrap 法［24］，其基本思想是依據(jù)小樣本信息來模擬未知分布，通過不斷獲取再生樣本從而實現(xiàn)小樣本轉換為大樣本。傳統(tǒng)的非參數(shù)Bootstrap 方法在建立先驗分布時存在抽樣誤差大的問題，為了減少此誤差，使經(jīng)驗分布函數(shù)的泛化能力更強，本文以威布爾分布作為先驗樣本的經(jīng)驗分布Fn，然后再采用Bootstrap 法從Fn中抽取再生樣本，重復抽取n次，對每組樣本分別使用最小二乘法進行威布爾分布擬合得到n組形狀參數(shù)ξ和尺度參數(shù)η。

（3）非參數(shù)核密度估計確定先驗分布。為了減少抽樣帶來的誤差，提高經(jīng)驗分布函數(shù)的泛化能力，采用非參數(shù)核密度估計法從抽樣數(shù)據(jù)特性出發(fā)擬合出對應的概率密度函數(shù)，計算出待估參數(shù)的先驗分布。核密度估計公式為［25］：

式中：K（?）為核函數(shù)，當樣本量比較大時，核函數(shù)的選擇對估計結果的影響不大，本文選取最為常用的高斯核函數(shù)作為核函數(shù)；h為窗寬，最佳窗寬的估計值為其中σ為樣本標準差，n為樣本數(shù)。

將Bootstrap 方法抽樣得到的n組形狀參數(shù)ξ和尺度參數(shù)η分別代入式（6）即可得到待估參數(shù)θ（形狀參數(shù)ξ或尺度參數(shù)η）的先驗分布公式為：

本文先驗分布確定方法有效降低了對評估樣本數(shù)據(jù)量的需求，同時也避免了因威布爾分布參數(shù)多而導致的計算復雜的問題。

2.3 后驗分布的計算 在獲得先驗分布后，結合現(xiàn)場試驗數(shù)據(jù)X=（x1，x2，…，xn），可得到后驗分布公式為：

根據(jù)以上后驗分布函數(shù)來計算待估參數(shù)的估計值時需要進行數(shù)值積分來逼近待估參數(shù)的期望和方差，本文后驗分布是非標準形式的密度函數(shù)，采用數(shù)值積分方法計算時會有高維積分出現(xiàn)，使得逼近誤差隨著維數(shù)的增加而增加，同時也增加了計算難度［26］。馬爾科夫鏈蒙特卡羅（Markov Chain Monte-Carlo，MC-MC）［27-29］方法是解決此類問題的一種行之有效的方法，基本思想是通過建立平穩(wěn)分布為g（θ|x）的Markov 鏈來得到后驗分布g（θ|x）的樣本，隨數(shù)值仿真模擬的變化而實時抽取隨機樣本，來動態(tài)模擬求取積分，基于這些樣本可做各種統(tǒng)計推斷。在貝葉斯方法中常見的構造馬爾科夫鏈的方法是Gibbs 抽樣，其常用來處理高維、非標準形式的后驗分布。在抽樣過程中后驗分布g（ξ，η|x）若給定η，則g（ξ，η|x）僅為ξ的函數(shù)，此時稱g（ξ|x，η）為參數(shù)ξ的滿條件分布。在Gibbs 抽樣過程中，從沒有顯式形式的滿條件分布中抽樣比較困難，故本文引入Metropolis 算法抽取隨機數(shù)，將Gibbs 抽樣和Metropolis 算法［30］結合起來采用混合Gibbs 抽樣的方法來構造馬爾科夫鏈。其中Gibbs 抽樣方法［31］研究已較為成熟，本文不再贅述，Metropolis算法［32-33］的實現(xiàn)方法如下：

（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抽樣同時產(chǎn)生多條Markov鏈，若這幾條鏈穩(wěn)定下來則認為混合Gibbs抽樣收斂。

2.4 可靠度評價 在實際工程中，往往需要給出在一定可靠度條件下參數(shù)的估計值，設參數(shù)估計的可靠度為R，分布函數(shù)與可靠度的關系為：

對上式取對數(shù)可得：

式中xR為兩參數(shù)威布爾分布在可靠度R條件下的估計值。將基于貝葉斯理論更新得到的兩參數(shù)后驗威布爾分布估計結果代入式（10），即可計算出在可靠度R條件下的砂礫料級配特征參數(shù)的估計值xR。此外根據(jù)貝葉斯方法更新后的威布爾分布作為總體分布，產(chǎn)生的隨機數(shù)可作為壓實質(zhì)量評估模型中壩料級配特征參數(shù)的輸入指標。

3 算例驗證

假設砂礫石料P5 含量服從形狀參數(shù)ξ=3、尺度參數(shù)η=2 的威布爾分布，從中產(chǎn)生25 個隨機數(shù)作為先驗樣本數(shù)據(jù)集，5 個隨機數(shù)作為現(xiàn)場試驗數(shù)據(jù)集，采用本文方法編制的程序對上述小樣本數(shù)據(jù)進行威布爾分布的貝葉斯估計求解。在對后驗分布進行模擬仿真求解時，本文構建了迭代次數(shù)為10 000、20 000、30 000 的三條Markov 鏈，仿真過程中監(jiān)控形狀參數(shù)ξ和尺度參數(shù)η輸出的模擬估計值，經(jīng)過統(tǒng)計分析用于判斷迭代過程是否收斂。三條馬爾科夫鏈輸出結果的統(tǒng)計分析結果如表1 所示。從表1 和圖2 可以看出，隨著迭代次數(shù)的增加，MC 誤差在逐步的減小，且三條鏈在2．5%、50%和97．5%這三個分位數(shù)對應的估值基本相同，這表明在迭代過程中Markov 鏈逐漸收斂并趨于穩(wěn)定。當Markov 鏈達到平穩(wěn)分布時，將迭代計算結果的均值作為后驗形狀參數(shù)和尺度參數(shù)的估計值，后驗形狀參數(shù)ξ=3．1950 和尺度參數(shù)η=2．0380 與真實值的相對誤差分別為6．5%和1．9%，在誤差允許范圍內(nèi)，由此驗證了本文小樣本條件下貝葉斯估計方法和編制程序的準確性。

圖2 待估參數(shù)的后驗分位數(shù)

表1 后驗分布參數(shù)的模擬估計值

上述通過一次蒙特卡羅模擬驗證了本文方法的有效性，為了消除隨機性的影響，先后進行了100次蒙特卡羅模擬，并分別采用了本文方法與傳統(tǒng)Bootstrap 方法進行貝葉斯估計對比分析。后驗分布時通過構建迭代次數(shù)為10 000 的Markov 鏈進行模擬仿真，計算得到的貝葉斯估計均值、與真實值的相對誤差和均方誤差如表2所示。從表2中可以看出，采用本文方法計算得到的貝葉斯估計均值、相對誤差和均方誤差都比采用傳統(tǒng)Bootstrap 方法計算得到的值小，說明本文方法相比傳統(tǒng)Bootstrap 方法具有較高的精度，這在一定程度上也說明了本文算法的有效性及穩(wěn)定性。

表2 不同方法對比分析

4 工程實例驗證與應用

4.1 項目背景及試驗數(shù)據(jù) 大石門水利樞紐工程位于新疆維吾爾自治區(qū)巴音郭楞蒙古族自治州且末縣車爾臣河干流之上，壩址位于車爾臣河出山口、與支流托其里薩依河匯合口下游300 m處。工程是車爾臣河流域規(guī)劃中確定的近期重點開發(fā)控制性樞紐工程，屬國家節(jié)水供水172 項重大水利工程之一，是一項承擔防洪、發(fā)電和灌溉等任務的綜合性水利樞紐工程。為了有效指導施工，保證大壩砂礫料回填質(zhì)量，相關單位對新疆車爾臣河大石門水利樞紐工程大壩填筑的砂礫料進行現(xiàn)場原型級配碾壓試驗，共獲得23組小樣本級配數(shù)據(jù)（級配曲線如圖3），此外還有4組為現(xiàn)場某一單元工程確定的挖坑試驗數(shù)據(jù)（級配曲線如圖4），根據(jù)級配曲線計算得到的P5 含量、最大粒徑以及曲率系數(shù)數(shù)據(jù)詳見表3。

表3 試驗數(shù)據(jù)

圖3 先驗小樣本顆粒級配曲線圖

圖4 現(xiàn)場試驗顆粒級配曲線圖

4.2 實例應用 首先采用最小二乘法對樣本容量為23的先驗小樣本級配參數(shù)數(shù)據(jù)進行兩參數(shù)威布爾分布擬合，參數(shù)的擬合結果如圖5所示，對應擬合分布的統(tǒng)計參數(shù)詳見表4。從表中可以看出各級配參數(shù)的相關系數(shù)均大于0．9，擬合結果較為可觀。在顯著性水平α=0．05條件下，計算得到的K-S統(tǒng)計量均小于臨界值D（23，0．05）=0．624，說明本研究所選用的威布爾分布模型能夠較好的適用于砂礫石壩料級配參數(shù)的描述及參數(shù)估計。

圖5 砂礫石料級配參數(shù)威布爾擬合

表4 擬合分布的統(tǒng)計參數(shù)

接著利用Bootstrap抽樣從上述級配參數(shù)對應的Fn（x）中抽取1000組再生子樣本，每組樣本的數(shù)量為23個，對每組樣本分別使用最小二乘法進行威布爾分布擬合得到1000組形狀參數(shù)ξ和尺度參數(shù)η。在對得到的抽樣結果進行統(tǒng)計分析時，為了克服結果對分布類型過分依賴而造成主觀性過大的問題，本文采用了非參數(shù)核密度估計法直接從抽樣結果的樣本特性出發(fā)擬合概率密度函數(shù)，計算得到1000組形狀參數(shù)和尺度參數(shù)的樣本標準差和最優(yōu)窗寬，進而得到待估參數(shù)的先驗密度函數(shù)，計算結果見表5。形狀參數(shù)和尺度參數(shù)對應的概率直方圖和密度函數(shù)圖如圖6所示。

圖6 核密度估計擬合

表5 核密度估計參數(shù)及密度函數(shù)

在求得先驗分布后，結合4組現(xiàn)場試驗數(shù)據(jù)對先驗分布進行修正，代入式（8）可以得到貝葉斯后驗分布表達式為：

針對上述非標準形式的Bayes后驗分布密度函數(shù)，為求解方便本文依據(jù)混合Gibbs抽樣原理編制了相應的Python程序對后驗分布進行抽樣，產(chǎn)生具有初始值的Markov鏈，經(jīng)過一段時間的迭代，當Markov鏈收斂到平穩(wěn)分布時即可求得待估參數(shù)的后驗分布樣本。為了判斷迭代過程是否收斂，分別構建了迭代次數(shù)為10 000、20 000、30 000的三條Markov鏈，通過監(jiān)控后驗形狀參數(shù)ξ和尺度參數(shù)η的自相關函數(shù)圖和抽樣迭代軌跡進行收斂判斷分析。限于篇幅，圖7僅列出P5含量的模擬結果。圖7（a）中自相關函數(shù)圖表示當前迭代結果與之前迭代結果的相關性，用以判斷收斂速度的快慢，從圖中可以發(fā)現(xiàn)在Lag=10左右樣本的自相關系數(shù)已接近于零，表明收斂速度很快。為提高模擬效率，本文每隔10個數(shù)取一個樣本進行統(tǒng)計分析，待估參數(shù)的統(tǒng)計模擬估計值見表6。從表中可以看出，隨著迭代次數(shù)的增加形狀參數(shù)ξ和尺度參數(shù)η的均值和標準差均在小范圍內(nèi)波動，MC誤差在逐步的減小，這表明在迭代過程中Markov鏈逐漸收斂并趨于穩(wěn)定。圖7（b）為形狀參數(shù)ξ和尺度參數(shù)η在迭代30 000次后的動態(tài)軌跡，從迭代軌跡圖中同樣可以發(fā)現(xiàn)Markov鏈基本收斂在一定區(qū)域內(nèi)，這時的抽樣結果可以認為是后驗分布的樣本，其對應的概率密度函數(shù)見圖7（c）。上述當Markov鏈達到平穩(wěn)分布時，將迭代計算結果的均值作為后驗形狀參數(shù)ξ和尺度參數(shù)η的估計值，可以得到砂礫石壩料級配參數(shù)經(jīng)貝葉斯方法更新后的威布爾分布概率分布函數(shù)以及在可靠度R條件下的壩料級配特征參數(shù)的估計值xR，詳見表7。

表7 貝葉斯方法更新后的威布爾分布概率分布函數(shù)

圖7 P5含量模擬結果

表6 后驗分布參數(shù)的模擬估計值

4.3 有效性評價 為了驗證本文貝葉斯方法更新后的威布爾分布作為總體分布的有效性，本文隨機抽取了大石門水利工程不同填筑單元的十組現(xiàn)場挖坑檢測級配數(shù)據(jù)X（如表8所示），運用K-S檢驗方法驗證現(xiàn)場數(shù)據(jù)X是否服從本文貝葉斯方法更新后的威布爾分布。首先建立原假設H0：F0（x）=Fn（x），其中Fn（x）為本文貝葉斯方法更新后的威布爾分布函數(shù)，F(xiàn)0（x）為樣本觀測值的累積分布函數(shù)：F0（xi）=i/n，i=1，2，…，n。然后計算K-S 檢驗統(tǒng)計量：Dn=max{|Fn（x）-F0（x）|}，計算結果如表9 所示，從表中可以看出在在顯著性水平α=0．05條件下，壩料各級配特征參數(shù)計算得到的K-S統(tǒng)計量均小于臨界值D（10，0．05）=0．410，說明隨機抽取的大石門工程不同填筑單元的級配數(shù)據(jù)服從本文貝葉斯方法更新后的威布爾分布函數(shù)，驗證了本文貝葉斯方法更新后的威布爾分布作為總體分布的有效性。

表8 大石門工程隨機抽取的現(xiàn)場級配特征參數(shù)數(shù)據(jù)

表9 K-S檢驗計算結果

在實際工程中，大壩填筑施工之前都會進行室內(nèi)外試驗，確定壩料合理的最大最小干密度、壩料級配特征參數(shù)，并且進一步確定大壩填筑碾壓施工參數(shù)與控制指標。但是在實際施工過程中常以有限的已碾壓單元挖坑試驗結果進行級配特征參數(shù)描述，無法獲得整個工作倉級配特征參數(shù)的總體分布規(guī)律，也無法為施工過程中壩料碾壓質(zhì)量的評價提供數(shù)據(jù)支撐。通過本文提出的小樣本條件下壩料級配特征參數(shù)的貝葉斯估計方法可以動態(tài)的獲取某填筑單元在一定可靠度條件下級配特征參數(shù)的估計值，以此作為評價某單元工程級配特征參數(shù)的合理依據(jù)。此外根據(jù)本文貝葉斯方法更新后的威布爾分布作為總體分布，產(chǎn)生隨機數(shù)作為壓實質(zhì)量評估模型中壩料級配特征參數(shù)的輸入指標用于壩料壓實質(zhì)量特性研究。

5 結論與展望

針對砂礫石壩，壩料的級配特征參數(shù)是確定壩料壓實特性、強度及變形參數(shù)所需要的最重要的指標之一，對于同一來源的砂礫石料，其顆粒組成具有明顯的規(guī)律性，服從一定的概率分布規(guī)律。為實現(xiàn)通過已有的砂礫石壩挖坑檢測小樣本級配參數(shù)數(shù)據(jù)來推斷構建總體分布，本文以大石門水利樞紐瀝青混凝土心墻砂礫石壩的壩殼料碾壓質(zhì)量挖坑檢測得到的小樣本級配數(shù)據(jù)為研究對象，根據(jù)貝葉斯理論將先驗小樣本數(shù)據(jù)和現(xiàn)場試驗數(shù)據(jù)進行融合，最終求解出了砂礫石料級配參數(shù)的后驗分布函數(shù)。得出以下結論：

（1）由于雙參數(shù)威布爾分布不存在共軛先驗分布，本文采用了參數(shù)化Bootstrap 方法和非參數(shù)核密度估計法確定了小樣本數(shù)據(jù)下威布爾分布參數(shù)的先驗分布，有效地解決了威布爾分布下貝葉斯公式難以求解的問題；針對后驗分布積分計算存在高維積分，計算量大的問題，本文采用了混合Gibbs抽樣方法對后驗分布進行模擬仿真求解，模擬精度較高，算例驗證了本文方法的有效性與穩(wěn)定性。

（2）通過算例對本文方法與傳統(tǒng)Bootstrap 方法進行貝葉斯評估進行了比較，結果表明本文方法精度更高，計算結果更為穩(wěn)定，說明了本文提出方法的優(yōu)越性與準確性。

（3）本研究通過隨機抽取大石門工程不同填筑單元的級配數(shù)據(jù)驗證了本文貝葉斯方法更新后的威布爾分布作為壩料級配特征參數(shù)總體分布的有效性，通過本文提出的小樣本條件下壩料級配特征參數(shù)的貝葉斯估計方法可以動態(tài)的獲取某填筑單元工程在一定可靠度條件下級配特征參數(shù)的估計值，以此作為評價某單元工程級配特征參數(shù)的合理依據(jù)。此外，本文貝葉斯方法更新后的威布爾分布可作為壩料級配特征參數(shù)的總體分布，這對于研究大壩填筑碾壓施工質(zhì)量控制、壓實后壩體材料力學參數(shù)估計及壩體變形分析與預測都具有重要的數(shù)據(jù)支撐作用。

土石壩筑壩材料的級配特征參數(shù)分布規(guī)律是一個具有很強實踐性的問題，因此本文貝葉斯方法將不斷把現(xiàn)場挖坑檢測獲得的真實級配數(shù)據(jù)融合到先驗分布中，使獲得的后驗威布爾分布更具有實踐性和真實性。需要注意的是，級配的表征其實是壩料級配特征參數(shù)（P5含量，曲率系數(shù)和最大粒徑）之間相關聯(lián)的聯(lián)合分布，但是由于計算的難度以及為了滿足實際工程切實的需求，本文將不同級配特征參數(shù)進行了單獨的統(tǒng)計分析，這也為后續(xù)深入開展級配特征參數(shù)間的聯(lián)合分布奠定了基礎。此外本文僅根據(jù)大石門水利樞紐工程現(xiàn)場挖坑檢測的砂礫石料級配數(shù)據(jù)進行了統(tǒng)計分析，后續(xù)將收集更多水利工程現(xiàn)場碾壓級配數(shù)據(jù)以建立不同地質(zhì)條件下壩料的先驗分布，由此獲取更準確的后驗分布估計。