基于贝叶斯定理的异重流水卷吸经验式不确定性研究

张维凯

浙江大学海洋学院,浙江 舟山 316000

泥沙异重流是水库和河流入海口地区常发生的现象。当密度较大的挟沙水流遇到密度较小的水体时下潜并沿底部运动形成。科学利用异重流可以给我们的生产生活带来很大的方便，如合理利用水库异重流排沙可以减小水库淤积；异重流在海洋形成的沉积区常常富含丰富的资源。国内外学者通过野外观测、实验和数值模拟对异重流取得了丰硕的研究成果[1-5]。特别地，异重流数值模拟由于具有很好的经济性和时效性，越来越得到广泛的应用。但是，异重流数学模型也存在一定的局限性，即一些物理量的计算要使用经验关系式，如异重流层平均模型中需要引用泥沙侵蚀经验式和水卷吸经验式[6]。而水卷吸经验式对异重流数学模型的计算结果影响很大，这是因为异重流与环境流体的密度差是其驱动力，卷吸水多，则异重流泥沙浓度低，反之则高，由此影响异重流与环境流体的密度差。经验式是基于实验数据，通过数学方法拟合得到的，存在很大的不确定性。这是因为，异重流发生在水下，难以观测到，实验室条件下的异重流具有较强的非恒定性，使得获取大量高精度的野外及实测数据来率定经验式存在很多困难[7]。

Wang和Cao提出一套基于贝叶斯定理的Metropolis-Hastings采样法（以下简称概率法），能够利用少量标准贯入实验的数据得到大量不排水条件下杨氏模量的等效样本[8,9]。采用这种方法的好处是可以获得经验系数的大量样本，供我们分析其统计特征，对于传统的经验系数率定方法，得到的固定的经验系数取值，可以预见的是，当实测数据增加或存在误差时，经验系数的取值将会发生变化。本文将这种方法应用于水卷吸经验式（式1）的不确定性分析中，并将采样系数由一个扩展为两个，首先使用概率法获得经验系数的大量样本（第一节），然后分析经验系数样本的统计特征。

1 水卷吸经验式不确定性研究方法

1.1 水卷吸经验式的介绍

Parker等提出了水卷吸经验式[10]：

其中，ew为水卷吸系数，A=0.00153、B＝0.0204为Parker率定的经验系数取值，Ri=Rgch/u2为理查德森数，R为有效重力，g为重力加速度，c为层平均泥沙浓度，h为厚度，u为层平均速度。

1.2 基于贝叶斯定理的Metropolis-Hastings采样法

本文采用Metropolis-Hastings采样法，对经验系数N1=B/A和N2=1/A进行采样，若假定经验系数N1服从均值为µ1，标准差为σ1的正态分布，经验系数N2服从均值为µ2，标准差为σ2的正态分布，那么等式2左边服从均值为µ1+µ2Ri，标准差为的正态分布。对采样过程叙述如下。

用Data指代实测值（即从实验中获得的异重流相关的Ri，ew数据），用m指代经验系数组合，m=N1,N2），用P(m/Data)表征给定实测数据时，经验系数某组取值的概率。根据贝叶斯定理，有[11-13]：

其中，下标j表示经验系数样本的序号，n为总样本数（n=8万），P(m)为先验分布（即无实测数据条件下，经验系数某组取值的概率），P(Data/m)为似然函数（即，对于给定的经验系数取值，公式1计算值与实测值之间的拟合程度），（3）式分母为一个归一化常量，在采样过程中不变化，无需计算，（4）、（5）、（6）式相乘的结果即为（3）式的分子。ξi为实测的1/ew值，nS为实测值组数（55组）。µ1min、µ1max、σ1min、σ1max分别为N1先验概率均值µ1的最小值、最大值，标准差σ1的最小值、最大值；µ2min、µ2max、σ2min、σ2max分别为N2先验概率均值µ2的最小值、最大值，标准差σ2的最小值、最大值。

首先任意选取一组经验系数值m0～(N1,0,N2,0)；分别以N1,0、N2,0为均值，以N1,0×σ1、N2,0×σ2为标准差，随机取得一组新的经验系数值样本，记为应用式（3）-（7）式计算即的接受概率。当ra大于一个[0,1]之间的随机数时，接受为新的起点重复上述过程，直至所得样本取值的概率分布趋于稳定[8]。

2 经验系数的样本分布

由于水卷吸系数应小于0.075[10]，故采样时，给定限制条件N1≥13.3，在得到N1～N2的样本后，通过简单的数学计算，转换成经验系数A-B。考虑到A和B相差一个量级，在统计样本的频数分布时，首先将样本A乘以1000，样本B乘以100，然后以每个样本点为圆心，0.2为半径作圆，统计圆内点的频数，如图1所示。然后统计依次频数大于2500，2490，2480……的样本数，最终得到频数大于22的样本总数为76403，得到总样本数95%的A-B样本所在的区间范围，使用同样方法统计得到25%、50%、75%样本，如图2，其中黑色十字为概率最大样本[A,B]＝[0.0016,0.0249]。颜色由深至浅，依次为25%、50%、75%、95%样本的区间范围。占总样本95%的A-B取值组合中，A-B的区间范围分别是概率最大A-B值的[51%,630%]和[63%,1200%]，具体范围见表1。

表1 不同比例经验系数样本所在区间范围Table 1 Interval range of different proportions of empirical coefficients

图1 经验系数A和B样本分布Fig.1 Sample distribution of empirical coefficients A and B

图2 经验系数A和B不同比例样本所在的区间范围Fig.2 Interval of different proportions of empirical coefficients A and B

图3 显示了概率最大系数组合对应的经验式与原经验式的对比情况，图3a为代入概率最大系数组合的公式1与实测数据的拟合图像，图3b为原经验式与实测数据的拟合图像。可以看出，概率最大系数组合与原公式拟合结果相近。

图3 概率法拟合的经验式与原经验式对比Fig.3 Comparison between the empirical formula fitted by probability method and the original formula

3 结论

本文应用基于贝叶斯定理Metropolis-Hastings采样法对异重流水卷吸经验式进行了不确定性分析，得到了其中经验系数组合A-B的大量样本，并对其进行了分析，结论如下：

（1）通过概率法拟合的最大概率经验系数组合[A,B]＝[0.0016,0.0249]与原经验式拟合结果相近，这也说明概率法在拟合经验式在保证传统拟合方法的精确度基础上，也具有更大的优势（能够获得经验系数其它取值的样本分布，继而得到概率密度）。

（2）95%经验系数A-B的样本区间范围比25%样本区间范围扩大明显（A扩大790%，B扩大832%），说明经验系数的不确定性很大，有可能会对异重流数值模拟的结果产生较大影响。

（3）将相同比例经验系数组合相对于最大概率系数组合的比例范围一起比较，可以看出，经验系数B区间范围远远大于A的范围，说明经验系数B的不确定性更大，在率定时要尤其注意。

值得注意的是，这种基于贝叶斯定理的Metropolis-Hastings采样法还可以应用于其它的公式，分析其中经验系数的不确定性，且通过公式与公式之间的对比，可以判断由于实测数据的误差或稀少造成的经验式不确定性的大小，当我们模拟异重流需要选取经验公式时，可以提供一定的参考。

[1]Meiburg E,Kneller B.Turbidity currents and their deposits[J].Annual Review of Fluid Mechanics,2010,42(1)：135-156

[1]严忠銮,安瑞冬,李 嘉,等.浊度型清浑水交界面识别方法及其在水库异重流观测中的应用[J].水利水电科技进展,2013,33(6)：71-75

[2]范家骅.浑水异重流水量掺混系数的研究[J].水利学报,2011,42(1)：19-26

[3]贺治国,林挺,赵亮,等.异重流在层结与非层结水体中沿斜坡运动的实验研究[J].中国科学：技术科学,2016,46(6)：570

[4]赵 琴,李 嘉,安瑞冬.水库浑水异重流的两相流模型适用性研究[J].水动力学研究与进展,2010,25(1)：76-84

[5]胡 鹏,胡元园,贺治国,等.泥沙异重流与环境物质交换经验式对比[J].水科学进展,2017,28(2)：257-264

[6]Wells M,Nadarajah P.The intrusion depth of density currents flowing into stratified water bodies[J].Journal of Physical Oceanography,2009,39(39)：1935-1947

[7]Wang Y,Cao Z.Probabilistic characterization of Young's modulus of soil using equivalent samples[J].Engineering Geology,2013,159：106-118

[8]曹子君,赵腾远,王宇,等.基于贝叶斯等效样本的土体杨氏模量的统计特征确定方法[J].防灾减灾工程学报,2015,35(5)：581-585

[9]Parker G,Fukushima Y,Pantin HM.Self-acceleratingturbiditycurrents[J].Journalof Fluid Mechanics,1986,171：145-181

[10]Traer MM,Hilley GE,Fildani A,etal.The sensitivity of turbidity currents to mass and momentum exchanges between these under flows and their surroundings[J].Journal of Geophysical Research Atmospheres,2012,117(F1)：1009T

[11]Hilley GE,Mynatt I,Pollard DD.Structural geometry of Raplee Ridge mono cline and thrust fault imaged using inverse Boundary Element Modeling and ALSM data[J].Journal of Structural Geology,2010,32(1)：45-58

[12]Ang AHS,Tang WH.Probability concepts in engineering：emphasis on applications to civil and environmental Engineering[M].Beijing：China Architecture&Building Press,2017