肖建峰王浩

推移質(zhì)泥沙顆粒躍移運動參數(shù)貝葉斯估計

肖建峰1，2王浩2

一、前言

在研究單顆粒推移質(zhì)泥沙運動規(guī)律方面，很多學者做了大量的工作。韓其為研究了單顆粒泥沙運動力學及統(tǒng)計規(guī)律，認為泥沙運動過程近似為馬爾科夫（Markov）過程。孫志林指出泥沙交換的統(tǒng)計規(guī)律基本上建立在時間離散的Markov鏈基礎之上，泥沙的運動和交換是時間連續(xù)的過程。徐俊峰指出床面層內(nèi)的泥沙狀態(tài)轉(zhuǎn)移屬于非周期不可約的Markov-MapkoB鏈。胡春宏利用高速攝影機對推移質(zhì)泥沙顆粒運動進行了大量試驗工作。侯暉昌依據(jù)水流底速符合正態(tài)分布，推導出顆粒運動參數(shù)符合卡方分布。Lee通過試驗獲得了顆粒躍移運動的無量綱參數(shù)并分析了顆粒躍移參數(shù)與摩阻流速之間的相關(guān)關(guān)系。Nino通過試驗獲得了顆粒無量綱參數(shù)躍長L/d、躍高H/d與泥沙剪切應力比值τb/τcr之間的關(guān)系。Sommerfeld et al.、Oesterle考慮了泥沙顆粒間碰撞運動，Sommerfeld等模擬了不同粒徑大小間的碰撞運動，并分析了泥沙間顆粒發(fā)生碰撞的頻率，導出了顆粒間碰撞發(fā)生的概率。Bialik討論了泥沙顆粒隨機運動數(shù)學模型，認為顆粒無量綱運動參數(shù)與摩阻流速、床面剪切應力比值之間成正相關(guān)的關(guān)系。

本文擬依據(jù)貝葉斯理論，對單顆粒運動參數(shù)進行貝葉斯估計和檢驗，進一步討論泥沙顆粒運動規(guī)律。

二、顆粒運動參數(shù)的統(tǒng)計規(guī)律

竇國仁推導出顆粒單步跳躍長度近似Γ分布，并根據(jù)Γ分布的可加性推導出顆粒跳躍n步運動參數(shù)的分布密度。并認為指數(shù)分布是γ=1時Γ分布，單步跳躍長度是指數(shù)分布，則單次跳躍距離分布也是指數(shù)分布。韓其為指出泥沙單步距離分布近似服從指數(shù)分布，即：

i=2，3，4，t為時間，即為μ=1/L單步距離的倒數(shù)，Ui為單步運動的平均速度。

因而可得出運動顆粒單步距離的密度函數(shù)為：

胡春宏根據(jù)大量實驗數(shù)據(jù)，依泥沙躍長和躍高的相對值，即L/L、H/H、Ui作為統(tǒng)計量，認為河床底部和加糙時均對上述統(tǒng)計分布參數(shù)無影響。并分析得出均符合Γ分布。

式中，α=7.51、β、α0為Γ分布的三個參數(shù)。

三、躍移參數(shù)分布的貝葉斯估計

1.樣本選取

貝葉斯理論認為，任一個未知量θ都可看作一個隨機變量，應用一個概率分布描述對的未知狀況。p（x|θ）表示隨機變量θ給定某一個值時，總體X的條件分布。一般的從先驗的分布π（θ）中產(chǎn)生一個樣本θ'，然后從總體分布p（x|θ'）中產(chǎn)生一個樣本x=（x1，…xn），其聯(lián)合概率密度為。無論是經(jīng)典統(tǒng)計還是貝葉斯理論，均承認似然函數(shù)L（θ'）是存在的，對于樣本觀察值x=（x1，…xn），總體和樣本關(guān)于參數(shù)θ的信息都包含在似然函數(shù)中，而再利用其做統(tǒng)計推斷時會產(chǎn)生較大差別。

推移質(zhì)泥沙顆粒運動試驗在水槽里進行，水槽尺寸為100m×1m×0.8m，顆粒粒徑范圍0.9～1.5mm。采用PIV系統(tǒng)采集泥沙顆粒運動軌跡，顆粒運動參數(shù)提取基于圖像識別的推移質(zhì)輸移檢測系統(tǒng)提取。顆粒泥沙運動軌跡和樣本散點分布見圖1、圖2。

圖1 顆粒典型躍移運動軌跡圖

圖2 顆粒運動參數(shù)L/d、H/d散點及誤差線圖

2.先驗分布的確定

對于先驗分布確定，大致可以分為利用先驗信息確定先驗分布及無信息先驗分布，本研究擬利用先驗信息確定。若推移質(zhì)泥沙顆粒單步運動參數(shù)樣本服從指數(shù)分布E（θ），其密度函數(shù)為：

表1 模擬參數(shù)統(tǒng)計表

圖3 迭代后參數(shù)誤差圖

圖4 觀測值及隨機樣本間誤差條形圖

圖5 無量綱參數(shù)L/d、H/d密度分布圖

式中：λ為無量綱數(shù)，D為泥沙顆粒粒徑值。

上式的似然函數(shù)為：

則后驗分布形式為：

可看出，該似然函數(shù)形式類似于伽瑪分布概率密度函數(shù)。因此，泥沙顆粒單步運動參數(shù)的先驗分布即為參數(shù)為Ga（α，β）的Γ分布，其密度函數(shù)為

式中，α、β為超產(chǎn)數(shù)。

因此，推移質(zhì)泥沙顆粒運動參數(shù)的先驗分布為指數(shù)分布的共軛分布，即伽瑪分布。

3.抽樣及參數(shù)估計

一般情況下，無法直接獲得參數(shù)獨立的后驗分布，且難以對聯(lián)合分布進行分解?？山柚R爾科夫鏈（MCMC）方法進行模擬。采用Gibbs抽樣算法對參數(shù)進行估計時，給定初始值θ0，然后從總體分布p（x|θ'）中產(chǎn)生新的隨機數(shù)，獲得馬爾科夫鏈。采取Gibbs抽樣進行迭代，經(jīng)過一段時間n次迭代后，可得到θn=（θ1n，……，θin），最終得到θ1，θ2，θ3，…。此時各時刻θn的邊際分布為平穩(wěn)分布，模型收斂。

為便于模型的收斂，模擬計算迭代次數(shù)為10萬次。整個迭代計算過程是基本穩(wěn)定的，達到收斂，模型收斂后，即可后驗分布的參數(shù)信息。收斂后樣本均值、方差，誤差范圍等參數(shù)統(tǒng)計見表1。參數(shù)誤差分布見圖3。

4.后驗分布確定

后驗分布確定后，采用隨機模擬（MonteCarlo）方法，產(chǎn)生符合后驗分布的隨機樣本，進行統(tǒng)計推斷。對觀測數(shù)據(jù)及隨機產(chǎn)生樣本無量綱參數(shù)L/d、H/d進行誤差對比分析，誤差范圍均在合理的區(qū)間，誤差條形圖見圖4。進一步采取隨機樣本對顆粒運動參數(shù)分布進行模擬，兩無量綱參數(shù)L/d、H/d分布密度見圖5。

四、結(jié)論

模擬過程表明，估計偏差較小、穩(wěn)定性較好，模型是收斂的，估計和檢測方法是適用的。

對單顆粒運動參數(shù)進行貝葉斯估計和檢驗表明，推移質(zhì)躍移顆粒運動無量綱參數(shù)符合指數(shù)分布。

影響顆粒運動的因素較多，且指數(shù)分布僅僅是伽馬分布的特例，因此需進一步討論更廣義的顆粒運動參數(shù)分布■

（作者單位：1.水利部淮河水利委員會2330012.河海大學210098）

國家自然科學基金面上項目（編號：51279046），江蘇省普通高校研究生科研創(chuàng)新計劃項目（KYZZ_0146），中央高?；究蒲谢穑∟o.2014B02714）