推移质泥沙颗粒跃移运动参数贝叶斯估计

肖建峰王浩

推移质泥沙颗粒跃移运动参数贝叶斯估计

肖建峰1，2王浩2

一、前言

在研究单颗粒推移质泥沙运动规律方面，很多学者做了大量的工作。韩其为研究了单颗粒泥沙运动力学及统计规律，认为泥沙运动过程近似为马尔科夫（Markov）过程。孙志林指出泥沙交换的统计规律基本上建立在时间离散的Markov链基础之上，泥沙的运动和交换是时间连续的过程。徐俊峰指出床面层内的泥沙状态转移属于非周期不可约的Markov-MapkoB链。胡春宏利用高速摄影机对推移质泥沙颗粒运动进行了大量试验工作。侯晖昌依据水流底速符合正态分布，推导出颗粒运动参数符合卡方分布。Lee通过试验获得了颗粒跃移运动的无量纲参数并分析了颗粒跃移参数与摩阻流速之间的相关关系。Nino通过试验获得了颗粒无量纲参数跃长L/d、跃高H/d与泥沙剪切应力比值τb/τcr之间的关系。Sommerfeld et al.、Oesterle考虑了泥沙颗粒间碰撞运动，Sommerfeld等模拟了不同粒径大小间的碰撞运动，并分析了泥沙间颗粒发生碰撞的频率，导出了颗粒间碰撞发生的概率。Bialik讨论了泥沙颗粒随机运动数学模型，认为颗粒无量纲运动参数与摩阻流速、床面剪切应力比值之间成正相关的关系。

本文拟依据贝叶斯理论，对单颗粒运动参数进行贝叶斯估计和检验，进一步讨论泥沙颗粒运动规律。

二、颗粒运动参数的统计规律

窦国仁推导出颗粒单步跳跃长度近似Γ分布，并根据Γ分布的可加性推导出颗粒跳跃n步运动参数的分布密度。并认为指数分布是γ=1时Γ分布，单步跳跃长度是指数分布，则单次跳跃距离分布也是指数分布。韩其为指出泥沙单步距离分布近似服从指数分布，即：

i=2，3，4，t为时间，即为μ=1/L单步距离的倒数，Ui为单步运动的平均速度。

因而可得出运动颗粒单步距离的密度函数为：

胡春宏根据大量实验数据，依泥沙跃长和跃高的相对值，即L/L、H/H、Ui作为统计量，认为河床底部和加糙时均对上述统计分布参数无影响。并分析得出均符合Γ分布。

式中，α=7.51、β、α0为Γ分布的三个参数。

三、跃移参数分布的贝叶斯估计

1.样本选取

贝叶斯理论认为，任一个未知量θ都可看作一个随机变量，应用一个概率分布描述对的未知状况。p（x|θ）表示随机变量θ给定某一个值时，总体X的条件分布。一般的从先验的分布π（θ）中产生一个样本θ'，然后从总体分布p（x|θ'）中产生一个样本x=（x1，…xn），其联合概率密度为。无论是经典统计还是贝叶斯理论，均承认似然函数L（θ'）是存在的，对于样本观察值x=（x1，…xn），总体和样本关于参数θ的信息都包含在似然函数中，而再利用其做统计推断时会产生较大差别。

推移质泥沙颗粒运动试验在水槽里进行，水槽尺寸为100m×1m×0.8m，颗粒粒径范围0.9～1.5mm。采用PIV系统采集泥沙颗粒运动轨迹，颗粒运动参数提取基于图像识别的推移质输移检测系统提取。颗粒泥沙运动轨迹和样本散点分布见图1、图2。

图1 颗粒典型跃移运动轨迹图

图2 颗粒运动参数L/d、H/d散点及误差线图

2.先验分布的确定

对于先验分布确定，大致可以分为利用先验信息确定先验分布及无信息先验分布，本研究拟利用先验信息确定。若推移质泥沙颗粒单步运动参数样本服从指数分布E（θ），其密度函数为：

表1 模拟参数统计表

图3 迭代后参数误差图

图4 观测值及随机样本间误差条形图

图5 无量纲参数L/d、H/d密度分布图

式中：λ为无量纲数，D为泥沙颗粒粒径值。

上式的似然函数为：

则后验分布形式为：

可看出，该似然函数形式类似于伽玛分布概率密度函数。因此，泥沙颗粒单步运动参数的先验分布即为参数为Ga（α，β）的Γ分布，其密度函数为

式中，α、β为超产数。

因此，推移质泥沙颗粒运动参数的先验分布为指数分布的共轭分布，即伽玛分布。

3.抽样及参数估计

一般情况下，无法直接获得参数独立的后验分布，且难以对联合分布进行分解。可借助马尔科夫链（MCMC）方法进行模拟。采用Gibbs抽样算法对参数进行估计时，给定初始值θ0，然后从总体分布p（x|θ'）中产生新的随机数，获得马尔科夫链。采取Gibbs抽样进行迭代，经过一段时间n次迭代后，可得到θn=（θ1n，……，θin），最终得到θ1，θ2，θ3，…。此时各时刻θn的边际分布为平稳分布，模型收敛。

为便于模型的收敛，模拟计算迭代次数为10万次。整个迭代计算过程是基本稳定的，达到收敛，模型收敛后，即可后验分布的参数信息。收敛后样本均值、方差，误差范围等参数统计见表1。参数误差分布见图3。

4.后验分布确定

后验分布确定后，采用随机模拟（MonteCarlo）方法，产生符合后验分布的随机样本，进行统计推断。对观测数据及随机产生样本无量纲参数L/d、H/d进行误差对比分析，误差范围均在合理的区间，误差条形图见图4。进一步采取随机样本对颗粒运动参数分布进行模拟，两无量纲参数L/d、H/d分布密度见图5。

四、结论

模拟过程表明，估计偏差较小、稳定性较好，模型是收敛的，估计和检测方法是适用的。

对单颗粒运动参数进行贝叶斯估计和检验表明，推移质跃移颗粒运动无量纲参数符合指数分布。

影响颗粒运动的因素较多，且指数分布仅仅是伽马分布的特例，因此需进一步讨论更广义的颗粒运动参数分布■

（作者单位：1.水利部淮河水利委员会2330012.河海大学210098）

国家自然科学基金面上项目（编号：51279046），江苏省普通高校研究生科研创新计划项目（KYZZ_0146），中央高校基本科研基金（No.2014B02714）