基于貝葉斯定理的異重流水卷吸經(jīng)驗(yàn)式不確定性研究

張維凱

浙江大學(xué)海洋學(xué)院,浙江 舟山 316000

泥沙異重流是水庫(kù)和河流入?？诘貐^(qū)常發(fā)生的現(xiàn)象。當(dāng)密度較大的挾沙水流遇到密度較小的水體時(shí)下潛并沿底部運(yùn)動(dòng)形成。科學(xué)利用異重流可以給我們的生產(chǎn)生活帶來(lái)很大的方便，如合理利用水庫(kù)異重流排沙可以減小水庫(kù)淤積；異重流在海洋形成的沉積區(qū)常常富含豐富的資源。國(guó)內(nèi)外學(xué)者通過(guò)野外觀測(cè)、實(shí)驗(yàn)和數(shù)值模擬對(duì)異重流取得了豐碩的研究成果[1-5]。特別地，異重流數(shù)值模擬由于具有很好的經(jīng)濟(jì)性和時(shí)效性，越來(lái)越得到廣泛的應(yīng)用。但是，異重流數(shù)學(xué)模型也存在一定的局限性，即一些物理量的計(jì)算要使用經(jīng)驗(yàn)關(guān)系式，如異重流層平均模型中需要引用泥沙侵蝕經(jīng)驗(yàn)式和水卷吸經(jīng)驗(yàn)式[6]。而水卷吸經(jīng)驗(yàn)式對(duì)異重流數(shù)學(xué)模型的計(jì)算結(jié)果影響很大，這是因?yàn)楫愔亓髋c環(huán)境流體的密度差是其驅(qū)動(dòng)力，卷吸水多，則異重流泥沙濃度低，反之則高，由此影響異重流與環(huán)境流體的密度差。經(jīng)驗(yàn)式是基于實(shí)驗(yàn)數(shù)據(jù)，通過(guò)數(shù)學(xué)方法擬合得到的，存在很大的不確定性。這是因?yàn)椋愔亓靼l(fā)生在水下，難以觀測(cè)到，實(shí)驗(yàn)室條件下的異重流具有較強(qiáng)的非恒定性，使得獲取大量高精度的野外及實(shí)測(cè)數(shù)據(jù)來(lái)率定經(jīng)驗(yàn)式存在很多困難[7]。

Wang和Cao提出一套基于貝葉斯定理的Metropolis-Hastings采樣法（以下簡(jiǎn)稱概率法），能夠利用少量標(biāo)準(zhǔn)貫入實(shí)驗(yàn)的數(shù)據(jù)得到大量不排水條件下楊氏模量的等效樣本[8,9]。采用這種方法的好處是可以獲得經(jīng)驗(yàn)系數(shù)的大量樣本，供我們分析其統(tǒng)計(jì)特征，對(duì)于傳統(tǒng)的經(jīng)驗(yàn)系數(shù)率定方法，得到的固定的經(jīng)驗(yàn)系數(shù)取值，可以預(yù)見(jiàn)的是，當(dāng)實(shí)測(cè)數(shù)據(jù)增加或存在誤差時(shí)，經(jīng)驗(yàn)系數(shù)的取值將會(huì)發(fā)生變化。本文將這種方法應(yīng)用于水卷吸經(jīng)驗(yàn)式（式1）的不確定性分析中，并將采樣系數(shù)由一個(gè)擴(kuò)展為兩個(gè)，首先使用概率法獲得經(jīng)驗(yàn)系數(shù)的大量樣本（第一節(jié)），然后分析經(jīng)驗(yàn)系數(shù)樣本的統(tǒng)計(jì)特征。

1 水卷吸經(jīng)驗(yàn)式不確定性研究方法

1.1 水卷吸經(jīng)驗(yàn)式的介紹

Parker等提出了水卷吸經(jīng)驗(yàn)式[10]：

其中，ew為水卷吸系數(shù)，A=0.00153、B＝0.0204為Parker率定的經(jīng)驗(yàn)系數(shù)取值，Ri=Rgch/u2為理查德森數(shù)，R為有效重力，g為重力加速度，c為層平均泥沙濃度，h為厚度，u為層平均速度。

1.2 基于貝葉斯定理的Metropolis-Hastings采樣法

本文采用Metropolis-Hastings采樣法，對(duì)經(jīng)驗(yàn)系數(shù)N1=B/A和N2=1/A進(jìn)行采樣，若假定經(jīng)驗(yàn)系數(shù)N1服從均值為μ1，標(biāo)準(zhǔn)差為σ1的正態(tài)分布，經(jīng)驗(yàn)系數(shù)N2服從均值為μ2，標(biāo)準(zhǔn)差為σ2的正態(tài)分布，那么等式2左邊服從均值為μ1+μ2Ri，標(biāo)準(zhǔn)差為的正態(tài)分布。對(duì)采樣過(guò)程敘述如下。

用Data指代實(shí)測(cè)值（即從實(shí)驗(yàn)中獲得的異重流相關(guān)的Ri，ew數(shù)據(jù)），用m指代經(jīng)驗(yàn)系數(shù)組合，m=N1,N2），用P(m/Data)表征給定實(shí)測(cè)數(shù)據(jù)時(shí)，經(jīng)驗(yàn)系數(shù)某組取值的概率。根據(jù)貝葉斯定理，有[11-13]：

其中，下標(biāo)j表示經(jīng)驗(yàn)系數(shù)樣本的序號(hào)，n為總樣本數(shù)（n=8萬(wàn)），P(m)為先驗(yàn)分布（即無(wú)實(shí)測(cè)數(shù)據(jù)條件下，經(jīng)驗(yàn)系數(shù)某組取值的概率），P(Data/m)為似然函數(shù)（即，對(duì)于給定的經(jīng)驗(yàn)系數(shù)取值，公式1計(jì)算值與實(shí)測(cè)值之間的擬合程度），（3）式分母為一個(gè)歸一化常量，在采樣過(guò)程中不變化，無(wú)需計(jì)算，（4）、（5）、（6）式相乘的結(jié)果即為（3）式的分子。ξi為實(shí)測(cè)的1/ew值，nS為實(shí)測(cè)值組數(shù)（55組）。μ1min、μ1max、σ1min、σ1max分別為N1先驗(yàn)概率均值μ1的最小值、最大值，標(biāo)準(zhǔn)差σ1的最小值、最大值；μ2min、μ2max、σ2min、σ2max分別為N2先驗(yàn)概率均值μ2的最小值、最大值，標(biāo)準(zhǔn)差σ2的最小值、最大值。

首先任意選取一組經(jīng)驗(yàn)系數(shù)值m0～(N1,0,N2,0)；分別以N1,0、N2,0為均值，以N1,0×σ1、N2,0×σ2為標(biāo)準(zhǔn)差，隨機(jī)取得一組新的經(jīng)驗(yàn)系數(shù)值樣本，記為應(yīng)用式（3）-（7）式計(jì)算即的接受概率。當(dāng)ra大于一個(gè)[0,1]之間的隨機(jī)數(shù)時(shí)，接受為新的起點(diǎn)重復(fù)上述過(guò)程，直至所得樣本取值的概率分布趨于穩(wěn)定[8]。

2 經(jīng)驗(yàn)系數(shù)的樣本分布

由于水卷吸系數(shù)應(yīng)小于0.075[10]，故采樣時(shí)，給定限制條件N1≥13.3，在得到N1～N2的樣本后，通過(guò)簡(jiǎn)單的數(shù)學(xué)計(jì)算，轉(zhuǎn)換成經(jīng)驗(yàn)系數(shù)A-B?？紤]到A和B相差一個(gè)量級(jí)，在統(tǒng)計(jì)樣本的頻數(shù)分布時(shí)，首先將樣本A乘以1000，樣本B乘以100，然后以每個(gè)樣本點(diǎn)為圓心，0.2為半徑作圓，統(tǒng)計(jì)圓內(nèi)點(diǎn)的頻數(shù)，如圖1所示。然后統(tǒng)計(jì)依次頻數(shù)大于2500，2490，2480……的樣本數(shù)，最終得到頻數(shù)大于22的樣本總數(shù)為76403，得到總樣本數(shù)95%的A-B樣本所在的區(qū)間范圍，使用同樣方法統(tǒng)計(jì)得到25%、50%、75%樣本，如圖2，其中黑色十字為概率最大樣本[A,B]＝[0.0016,0.0249]。顏色由深至淺，依次為25%、50%、75%、95%樣本的區(qū)間范圍。占總樣本95%的A-B取值組合中，A-B的區(qū)間范圍分別是概率最大A-B值的[51%,630%]和[63%,1200%]，具體范圍見(jiàn)表1。

表1 不同比例經(jīng)驗(yàn)系數(shù)樣本所在區(qū)間范圍Table 1 Interval range of different proportions of empirical coefficients

圖1 經(jīng)驗(yàn)系數(shù)A和B樣本分布Fig.1 Sample distribution of empirical coefficients A and B

圖2 經(jīng)驗(yàn)系數(shù)A和B不同比例樣本所在的區(qū)間范圍Fig.2 Interval of different proportions of empirical coefficients A and B

圖3 顯示了概率最大系數(shù)組合對(duì)應(yīng)的經(jīng)驗(yàn)式與原經(jīng)驗(yàn)式的對(duì)比情況，圖3a為代入概率最大系數(shù)組合的公式1與實(shí)測(cè)數(shù)據(jù)的擬合圖像，圖3b為原經(jīng)驗(yàn)式與實(shí)測(cè)數(shù)據(jù)的擬合圖像?？梢钥闯?，概率最大系數(shù)組合與原公式擬合結(jié)果相近。

圖3 概率法擬合的經(jīng)驗(yàn)式與原經(jīng)驗(yàn)式對(duì)比Fig.3 Comparison between the empirical formula fitted by probability method and the original formula

3 結(jié)論

本文應(yīng)用基于貝葉斯定理Metropolis-Hastings采樣法對(duì)異重流水卷吸經(jīng)驗(yàn)式進(jìn)行了不確定性分析，得到了其中經(jīng)驗(yàn)系數(shù)組合A-B的大量樣本，并對(duì)其進(jìn)行了分析，結(jié)論如下：

（1）通過(guò)概率法擬合的最大概率經(jīng)驗(yàn)系數(shù)組合[A,B]＝[0.0016,0.0249]與原經(jīng)驗(yàn)式擬合結(jié)果相近，這也說(shuō)明概率法在擬合經(jīng)驗(yàn)式在保證傳統(tǒng)擬合方法的精確度基礎(chǔ)上，也具有更大的優(yōu)勢(shì)（能夠獲得經(jīng)驗(yàn)系數(shù)其它取值的樣本分布，繼而得到概率密度）。

（2）95%經(jīng)驗(yàn)系數(shù)A-B的樣本區(qū)間范圍比25%樣本區(qū)間范圍擴(kuò)大明顯（A擴(kuò)大790%，B擴(kuò)大832%），說(shuō)明經(jīng)驗(yàn)系數(shù)的不確定性很大，有可能會(huì)對(duì)異重流數(shù)值模擬的結(jié)果產(chǎn)生較大影響。

（3）將相同比例經(jīng)驗(yàn)系數(shù)組合相對(duì)于最大概率系數(shù)組合的比例范圍一起比較，可以看出，經(jīng)驗(yàn)系數(shù)B區(qū)間范圍遠(yuǎn)遠(yuǎn)大于A的范圍，說(shuō)明經(jīng)驗(yàn)系數(shù)B的不確定性更大，在率定時(shí)要尤其注意。

值得注意的是，這種基于貝葉斯定理的Metropolis-Hastings采樣法還可以應(yīng)用于其它的公式，分析其中經(jīng)驗(yàn)系數(shù)的不確定性，且通過(guò)公式與公式之間的對(duì)比，可以判斷由于實(shí)測(cè)數(shù)據(jù)的誤差或稀少造成的經(jīng)驗(yàn)式不確定性的大小，當(dāng)我們模擬異重流需要選取經(jīng)驗(yàn)公式時(shí)，可以提供一定的參考。

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

[1]嚴(yán)忠鑾,安瑞冬,李 嘉,等.濁度型清渾水交界面識(shí)別方法及其在水庫(kù)異重流觀測(cè)中的應(yīng)用[J].水利水電科技進(jìn)展,2013,33(6)：71-75

[2]范家驊.渾水異重流水量摻混系數(shù)的研究[J].水利學(xué)報(bào),2011,42(1)：19-26

[3]賀治國(guó),林挺,趙亮,等.異重流在層結(jié)與非層結(jié)水體中沿斜坡運(yùn)動(dòng)的實(shí)驗(yàn)研究[J].中國(guó)科學(xué)：技術(shù)科學(xué),2016,46(6)：570

[4]趙 琴,李 嘉,安瑞冬.水庫(kù)渾水異重流的兩相流模型適用性研究[J].水動(dòng)力學(xué)研究與進(jìn)展,2010,25(1)：76-84

[5]胡 鵬,胡元園,賀治國(guó),等.泥沙異重流與環(huán)境物質(zhì)交換經(jīng)驗(yàn)式對(duì)比[J].水科學(xué)進(jìn)展,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]曹子君,趙騰遠(yuǎn),王宇,等.基于貝葉斯等效樣本的土體楊氏模量的統(tǒng)計(jì)特征確定方法[J].防災(zāi)減災(zāi)工程學(xué)報(bào),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