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

肖建峰王浩

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

肖建峰1，2王浩2

一、前言

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

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

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

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

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

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

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

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

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

1.樣本選取

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

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

圖1 顆粒典型躍移運(yùn)動軌跡圖

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

2.先驗(yàn)分布的確定

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

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

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

圖4 觀測值及隨機(jī)樣本間誤差條形圖

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

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

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

則后驗(yàn)分布形式為：

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

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

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

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

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

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

4.后驗(yàn)分布確定

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

四、結(jié)論

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

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

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

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

國家自然科學(xué)基金面上項(xiàng)目（編號：51279046），江蘇省普通高校研究生科研創(chuàng)新計(jì)劃項(xiàng)目（KYZZ_0146），中央高校基本科研基金（No.2014B02714）