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

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

1 研究背景

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

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

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

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

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

2 研究框架及方法介紹

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

圖1 研究框架

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

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

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

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

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

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

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

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

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

其密度函數(shù)為：

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

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

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

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

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

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

本文先驗(yàn)分布確定方法有效降低了對(duì)評(píng)估樣本數(shù)據(jù)量的需求，同時(shí)也避免了因威布爾分布參數(shù)多而導(dǎo)致的計(jì)算復(fù)雜的問(wèn)題。

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

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

（1）初始化：t=0，選擇一個(gè)對(duì)稱提議分布q（x），滿足q（x|y）=q（y|x），給定一個(gè)起始樣本點(diǎn)xt，迭代終止值為T；

（2）令t=t+1，從q（x|xt-1）中生成候選樣本x′；

（3）計(jì)算接收概率α，α=min{1，x′/xt-1}；

（4）從均勻分布中抽取隨機(jī)αt，若αt≤α，則接受候選樣本，xt=x′；否則，拒絕候選樣本，并令xt=xt-1；

（5）重復(fù)（2）—（4）步直到迭代終止。

采用混合Gibbs抽樣同時(shí)產(chǎn)生多條Markov鏈，若這幾條鏈穩(wěn)定下來(lái)則認(rèn)為混合Gibbs抽樣收斂。

2.4 可靠度評(píng)價(jià) 在實(shí)際工程中，往往需要給出在一定可靠度條件下參數(shù)的估計(jì)值，設(shè)參數(shù)估計(jì)的可靠度為R，分布函數(shù)與可靠度的關(guān)系為：

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

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

3 算例驗(yàn)證

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

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

表1 后驗(yàn)分布參數(shù)的模擬估計(jì)值

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

表2 不同方法對(duì)比分析

4 工程實(shí)例驗(yàn)證與應(yīng)用

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

表3 試驗(yàn)數(shù)據(jù)

圖3 先驗(yàn)小樣本顆粒級(jí)配曲線圖

圖4 現(xiàn)場(chǎng)試驗(yàn)顆粒級(jí)配曲線圖

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

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

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

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

圖6 核密度估計(jì)擬合

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

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

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

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

圖7 P5含量模擬結(jié)果

表6 后驗(yàn)分布參數(shù)的模擬估計(jì)值

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

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

表9 K-S檢驗(yàn)計(jì)算結(jié)果

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

5 結(jié)論與展望

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

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

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

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

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