異重流水卷吸經(jīng)驗(yàn)式不確定性對(duì)層平均數(shù)學(xué)模型的影響

（浙江大學(xué) 港口海岸與近海工程研究所，浙江舟山316000）

泥沙異重流常發(fā)生在水庫(kù)和高含沙河流等區(qū)域，如長(zhǎng)江深水航道區(qū)浮泥沿河床運(yùn)動(dòng)形成的異重流現(xiàn)象[1].異重流在人類(lèi)的生產(chǎn)生活中有著廣泛的應(yīng)用，合理利用異重流排沙可以減小水庫(kù)淤積和航道淤塞.國(guó)內(nèi)外學(xué)者通過(guò)野外觀測(cè)、實(shí)驗(yàn)和數(shù)值模擬對(duì)異重流取得了豐碩的研究成果[2－5].特別是異重流數(shù)值模擬，隨著計(jì)算機(jī)科學(xué)的進(jìn)步，取得了飛速的發(fā)展.層平均模型在異重流的數(shù)值模擬研究中得到了廣泛的應(yīng)用[6－11].層平均模型不可避免地要引用泥沙侵蝕、水卷吸和泥沙沉速等經(jīng)驗(yàn)關(guān)系式.其中，水卷吸經(jīng)驗(yàn)關(guān)系式的精確度對(duì)數(shù)學(xué)模型的計(jì)算結(jié)果影響很大.這是因?yàn)?，異重流的?qū)動(dòng)力來(lái)自其與環(huán)境流體之間的密度差，而水卷吸經(jīng)驗(yàn)式的計(jì)算結(jié)果代表環(huán)境水體進(jìn)入異重流的水量，若卷吸水較多，則使得含沙量較低，反之則含沙量高，并因此影響密度差和驅(qū)動(dòng)力.

早在1959年，Ellison等[12]建立了異重流卷吸水量與理查德森數(shù)的經(jīng)驗(yàn)式.在此基礎(chǔ)上，范家驊[2]認(rèn)為急流時(shí)泥沙異重流的水卷吸同鹽水異重流類(lèi)似.Parker等[13]結(jié)合少量的持續(xù)入流式異重流水卷吸實(shí)驗(yàn)數(shù)據(jù)，對(duì)Ellison經(jīng)驗(yàn)式進(jìn)一步修正，修正公式在異重流層平均模擬中得到了廣泛應(yīng)用.然而，異重流發(fā)生在水下，不易觀測(cè)到，且其與背景流體的交界面環(huán)境非常復(fù)雜（湍流、開(kāi)爾文－赫姆霍茨（K－H）不穩(wěn)定性），要量化這個(gè)區(qū)域是個(gè)挑戰(zhàn)，需要獲取大量高精度野外及實(shí)驗(yàn)數(shù)據(jù)來(lái)進(jìn)一步率定水卷吸經(jīng)驗(yàn)式.當(dāng)實(shí)測(cè)數(shù)據(jù)增加或存在誤差時(shí)，采用傳統(tǒng)最小二乘法擬合得到的經(jīng)驗(yàn)系數(shù)值極可能發(fā)生變化，這使得異重流水卷吸經(jīng)驗(yàn)式具有較大的不確定性.比如，Ottolenghi等[14]指出，Parker的水卷吸經(jīng)驗(yàn)式[15]可能會(huì)低估開(kāi)閘式異重流水卷吸量一個(gè)量級(jí).

Traer等[11]采用基于貝葉斯定理的 Monte Carlo方法（以下簡(jiǎn)稱(chēng)概率法），首次嘗試了對(duì)異重流水卷吸經(jīng)驗(yàn)式進(jìn)行不確定性分析，并探討了其不確定性對(duì)異重流經(jīng)典四方程模型的影響.Traer所采用的概率法相較于以最小二乘法為代表的傳統(tǒng)擬合方法有很大優(yōu)勢(shì)，因?yàn)橥ㄟ^(guò)經(jīng)驗(yàn)系數(shù)樣本與實(shí)測(cè)數(shù)據(jù)擬合優(yōu)劣的比較篩選過(guò)程，可以得到經(jīng)驗(yàn)系數(shù)的大量樣本，用以分析其均值和標(biāo)準(zhǔn)差等統(tǒng)計(jì)特性，而傳統(tǒng)擬合方法得到的是固定的經(jīng)驗(yàn)系數(shù)取值.然而，Traer所采用的異重流四方程模型具有恒定流和忽略地形沖淤過(guò)程的假設(shè).一方面，恒定流的假設(shè)顯然不能全面地描述異重流隨時(shí)間的變化；另一方面，異重流運(yùn)動(dòng)過(guò)程中造成的地形沖淤，以及地形沖淤對(duì)異重流演化的反作用不能輕易忽略[9－10].Traer等采用四方程模型的可能原因是：異重流三方程模型（包括水沙混合體質(zhì)量守恒、動(dòng)量守恒和泥沙質(zhì)量守恒）模擬自加速異重流時(shí)高估了泥沙侵蝕量，消耗過(guò)多能量，產(chǎn)能小于耗能從而不滿(mǎn)足能量平衡，與異重流自加速本質(zhì)相反.然而，Hu等[16]通過(guò)理論推導(dǎo)和數(shù)值實(shí)驗(yàn)發(fā)現(xiàn)，三方程模型并沒(méi)有破壞能量平衡.實(shí)際上，泥沙侵蝕量越多，異重流傳播時(shí)有越多勢(shì)能轉(zhuǎn)化為動(dòng)能，增加產(chǎn)能，滿(mǎn)足能量平衡.

因此，本文考慮異重流非恒定特性、地形沖淤和地形沖淤對(duì)異重流演化的影響，研究異重流水卷吸經(jīng)驗(yàn)式對(duì)三方程模型框架下異重流層平均水沙耦合數(shù)學(xué)模型計(jì)算結(jié)果的影響.具體而言，本文結(jié)合現(xiàn)有異重流水卷吸實(shí)驗(yàn)數(shù)據(jù)，應(yīng)用基于貝葉斯定理的Gibbs隨機(jī)采樣方法，采樣獲取了異重流水卷吸經(jīng)驗(yàn)式中經(jīng)驗(yàn)系數(shù)組合的大量樣本；然后，將經(jīng)驗(yàn)系數(shù)樣本值代入異重流層平均水沙耦合異重流數(shù)學(xué)模型（以下簡(jiǎn)稱(chēng)耦合模型），模擬陡坡上異重流演化過(guò)程，從而討論經(jīng)驗(yàn)系數(shù)不確定性對(duì)模型計(jì)算結(jié)果的影響.

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

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

Parker等[13]基于異重流實(shí)驗(yàn)數(shù)據(jù)，得到了水卷吸經(jīng)驗(yàn)式：

式中：ew為水卷吸系數(shù)；E1＝2.4，E2＝0.5，為Parker率定的經(jīng)驗(yàn)系數(shù)取值；Ri＝Rgφh／u2為理查德森數(shù)，R為有效重力，g為重力加速度，φ為層平均泥沙的體積分?jǐn)?shù)，h為層平均厚度，u為層平均速度.

1.2 基于貝葉斯定理的Gibbs采樣法

用Data指代實(shí)測(cè)值（即Parker等[13]在可侵蝕底床上所做持續(xù)入流式異重流實(shí)驗(yàn)中獲得的相關(guān)Ri和ew的數(shù)據(jù)，共53組），用m指代經(jīng)驗(yàn)系數(shù)組合，m＝ （E1，E2），用P（m｜Data）表征給定實(shí)測(cè)數(shù)據(jù)時(shí)，經(jīng)驗(yàn)系數(shù)某組取值的概率.根據(jù)貝葉斯定理[17]有

式中：下標(biāo)j表示經(jīng)驗(yàn)系數(shù)樣本的序號(hào)；nT＝2×105為總樣本數(shù)；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è)值之間的擬合程度.由于每個(gè)經(jīng)驗(yàn)系數(shù)的取值受到實(shí)驗(yàn)條件和測(cè)量數(shù)據(jù)準(zhǔn)確性等很多因素的影響，符合中心極限定理，所以本文假設(shè)經(jīng)驗(yàn)系數(shù)取值符合正態(tài)分布，采用 Gibbs采樣法[18－19]，分別對(duì)E1和E2進(jìn)行循環(huán)采樣.采樣過(guò)程敘述如下.

任意選取一組經(jīng)驗(yàn)系數(shù)值 （E1，0，E2，0）作為初始值；先固定E2，0不變，以E1，0為均值，以E1，0σE1為標(biāo)準(zhǔn)差（σE1為E1，0先驗(yàn)概率的標(biāo)準(zhǔn)差），隨機(jī)取得一個(gè)新的備選系數(shù)樣本值，將其記為E＊1，則E＊1的接受概率[20]為

式中：f（E1，0｜E＊1）為以E＊1為均值，E＊1σE1為標(biāo)準(zhǔn)差的正態(tài)分布在E1，0處的值.當(dāng)ra大于一個(gè)[0，1]之間的隨機(jī)數(shù)時(shí)，接受E＊1并以（E＊1，E2，0）為新的起點(diǎn).固定E＊1，使用同樣的方法對(duì)E2，0進(jìn)行采樣.如此循環(huán)，在采樣的過(guò)程中，不斷采取以下方法對(duì)樣本頻次進(jìn)行統(tǒng)計(jì)：以每個(gè)樣本點(diǎn)為圓心，0.2為半徑作圓，統(tǒng)計(jì)圓內(nèi)點(diǎn)的頻次.結(jié)果表明，當(dāng)采樣所得的（E1，E2）樣本數(shù)達(dá)到2×105后，繼續(xù)增加樣本數(shù)量對(duì)樣本集合的頻次分布不再產(chǎn)生影響，本文認(rèn)為采樣已經(jīng)收斂，停止采樣.

Gibbs法的關(guān)聯(lián)性采樣體現(xiàn)在兩個(gè)方面：① 由于每個(gè)樣本的入選僅與上一個(gè)樣本相關(guān)聯(lián)，所以在采樣過(guò)程中，如果待選（E1，E2）樣本能夠比上一個(gè)樣本更好地?cái)M合實(shí)測(cè)水卷吸數(shù)據(jù)（似然函數(shù)值較高），則該樣本能夠以更高的概率被接受；② 每組（E1，E2）樣本之間相關(guān)聯(lián)，如果一組（E1，E2）樣本使得水卷吸計(jì)算值和實(shí)驗(yàn)值差距較大（似然函數(shù)小），則該組樣本的接受概率就小.基于上述認(rèn)識(shí)可以預(yù)期的是，在通過(guò)足夠數(shù)量的采樣后，能夠使得水卷吸實(shí)驗(yàn)值和計(jì)算值擬合越好的（E1，E2）樣本頻次越高.

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

得到樣本點(diǎn)在（E1，E2）平面的分布后，進(jìn)一步得到25%nT、50%nT、75%nT和95%nT的（E1，E2）樣本所在的區(qū)間范圍，如圖1所示，其中黑色十字為概率最大樣本（E1，E2）＝（2.38，0.52）.顏色由深至淺，依次為25%nT、50%nT、75%nT、95%nT樣本的區(qū)間范圍.占總樣本95%的E1、E2取值組合中，E1和E2的區(qū)間范圍分別是[0.96Em1，2.72Em1]和[3×10－4Em2，6.59Em2]，Em1和Em2分別為概率最大的E1和E2值（即以該點(diǎn)為圓心，0.2為半徑的圓內(nèi)，樣本頻次最高的值）.

圖1 E1和E2不同比例樣本所在的區(qū)間范圍Fig.1 Interval of different proportions of E1and E2

圖2顯示了概率最大系數(shù)組合與最小二乘法擬合結(jié)果（E1，E2）＝（2.38，0.52）的對(duì)比情況，圖2（a）為代入概率最大系數(shù)組合的式（1）與實(shí)測(cè)數(shù)據(jù)的擬合圖像，圖2（b）為原經(jīng)驗(yàn)式與實(shí)測(cè)數(shù)據(jù)的擬合圖像.可以看出，概率最大系數(shù)組合與原公式擬合結(jié)果相近.

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

3 系數(shù)不確定性對(duì)數(shù)學(xué)模擬的影響

3.1 層平均水沙耦合異重流數(shù)學(xué)模型

在層平均假設(shè)下，非恒定的水沙耦合異重流數(shù)學(xué)模型可寫(xiě)成如下守恒形式[9－10]：

式中：U為守恒變量向量；F為通量向量；Sb為幾何源項(xiàng)向量；Sf為底床摩擦、底床變形、水體卷吸等源項(xiàng)向量；z為底床高程；E＝ωEs和D＝ωφb分別為泥沙上揚(yáng)通量和沉降通量，ω為泥沙沉速，Es＝8.9×10－9Z5m／（1＋8.9×10－9Z5m／0.3）為泥沙侵蝕系數(shù)[21]為泥沙雷諾數(shù)，R＝ （ρs－ρw）／ρ為有效重力，g為重力加速度，Ds為泥沙粒徑，ν為水的運(yùn)動(dòng)黏滯系數(shù)，φb＝r0φ為近底泥沙的體積分?jǐn)?shù)，r0為近底與水深平均泥沙體積分?jǐn)?shù)的比值；x為水平坐標(biāo)；t為時(shí)間；g′＝Rgφ為水下有效重力加速度為摩阻流速，cD為拖曳系數(shù)，rw為異重流上下界面阻力比值；ρ＝ρw（1－φ）＋ρsφ為水沙混合體密度，ρ0＝ρwp＋ρs（1－p），ρw為清水密度，ρs為泥沙密度；p為床沙孔隙率.控制方程使用有限體積法進(jìn)行離散[9－10]，并采用正方形網(wǎng)格，則有

式中：Δx為空間步長(zhǎng)；Δt為時(shí)間步長(zhǎng)；上標(biāo)n和 ＊代表時(shí)間節(jié)點(diǎn)號(hào)；i代表空間節(jié)點(diǎn)號(hào)；Fi＋1／2和Fi－1／2為體積單元間的通量向量；Sb，i為底坡源項(xiàng)在第i個(gè)空間單元的計(jì)算值.時(shí)間步長(zhǎng)需滿(mǎn)足CFL收斂條件.

本文將Parker推導(dǎo)的陡坡上異重流自加速臨界“點(diǎn)火條件”作為邊界條件，在這些初值條件下的異重流，從底床和環(huán)境水體獲得的能量高于其耗散，在沿程運(yùn)動(dòng)時(shí)會(huì)不斷加速，故稱(chēng)為“自加速”異重流，這種異重流在實(shí)際生活中應(yīng)用廣泛[8－9，11]：x＝0處，h＝2.0m，u＝0.801m／s，φ＝0.006 09.其他參數(shù)：g＝9.8m／s2，ω＝ 0.75cm／s，r0＝ 1.6[13]，cD＝0.004，Ds＝0.01cm，ν＝1×10－6m2／s，ρw＝1 000 kg／m3，ρs＝2 650kg／m3，p＝ 0.4，rw＝0.43.計(jì)算區(qū)域?yàn)? 000m，空間步長(zhǎng)1m，時(shí)間步長(zhǎng)1s，出口為自由出流.

3.2 （E1，E2）組合對(duì)異重流模型的影響分析

由于不同的（E1，E2）樣本值代表了不同的環(huán)境流體進(jìn)入異重流的清水量，進(jìn)而影響異重流與環(huán)境流體的密度差及其運(yùn)動(dòng)速度，為了使所有單個(gè)異重流算例均能流過(guò)最后一個(gè)統(tǒng)計(jì)斷面（1 000m），且盡可能提高計(jì)算效率，現(xiàn)隨機(jī)選擇5 000組樣本輸入模型并記錄其運(yùn)動(dòng)時(shí)間和運(yùn)動(dòng)距離.通過(guò)測(cè)試，當(dāng)計(jì)算時(shí)間為2 200s，計(jì)算區(qū)域?yàn)? 000m時(shí)，所有算例均符合要求.

由于篇幅所限，現(xiàn)選擇概率最大經(jīng)驗(yàn)系數(shù)值所對(duì)應(yīng)的異重流在4個(gè)典型時(shí)刻（100，200，300，400 s）的厚度h、速度u、泥沙體積分?jǐn)?shù)φ、底床形變?chǔ)b（初始時(shí)刻時(shí)為0），如圖3所示.其中，x表示與進(jìn)口邊界的距離.

圖3 不同時(shí)刻異重流水力參數(shù)Fig.3 Turbidity current hydraulic parameters at different times

從圖3中可以看出：隨著異重流運(yùn)動(dòng)時(shí)間的增加，厚度和速度增加明顯，且增速在加快；泥沙體積分?jǐn)?shù)先增加后減??；而底床被侵蝕得越來(lái)越多.

將所有（E1，E2）樣本值輸入上述異重流耦合模型，計(jì)算異重流過(guò)程，以x＝100m斷面為例，將計(jì)算所得水力參數(shù)范圍（0.53m≤h≤6.19m，0.98 m／s≤u≤2.30m／s，0.01≤φ≤0.11，－6.18m≤Δzb≤－2.91m）分為50等份，對(duì)每個(gè)區(qū)間作頻次（n）統(tǒng)計(jì)，結(jié)果如圖4所示.若計(jì)算水力參數(shù)的范圍越大，（E1，E2）不確定性對(duì)異重流計(jì)算結(jié)果的影響就越大.從圖4可以看出，厚度最小值頻次較高，速度和泥沙體積分?jǐn)?shù)的較大值頻次高.比如，厚度區(qū)間0.53～0.64m的頻次為69 981，占比約37%，隨著厚度的增大，頻次逐漸降低，計(jì)算速度和泥沙體積分?jǐn)?shù)有類(lèi)似但相反的規(guī)律.底床形變的最小值和中間值（－4m左右）頻次都較高.－6.19～－6.12m區(qū)間的頻次為14 159，占比7%；－4.22～4.15m區(qū)間的頻次為33 947，占比18%；其他區(qū)間占比普遍低于3%.

為解釋這一現(xiàn)象，并了解樣本空間內(nèi)不同樣本輸入模型計(jì)算所得的異重流水力參數(shù)分布情況，將E1樣本范圍[－3，7]以及E2樣本范圍[0，4]等分為100份，得到10 201組（E1，E2）.首先輸入耦合模型，得到10 201組計(jì)算結(jié)果，并取x＝100m斷面處計(jì)算結(jié)果作等值線，結(jié)果如圖5所示；然后輸入式（1）得到的10 201組ew計(jì)算值，作等值線，結(jié)果如圖6所示.

以圖5（a）為例，根據(jù)0.9m厚度等值線，將E1－E2平面分為兩個(gè)區(qū)域，區(qū)域I中ew較小，區(qū)域II中ew較大（圖6）.當(dāng)（E1，E2）取區(qū)域I樣本值時(shí)（該區(qū)間樣本頻次為97 357，占比達(dá)51%），ew較小，計(jì)算的厚度相應(yīng)較?。划?dāng)（E1，E2）取區(qū)域II樣本值時(shí)（頻次93 214，占比達(dá)48%），ew較大，計(jì)算所得厚度較大，但是由于區(qū)域I的厚度范圍（0.53～0.9m）較區(qū)域II的厚度范圍（0.9～6.19m）小，所以頻率分布圖4（a）中出現(xiàn)了最小值頻次最高的現(xiàn)象.以異重流速度等值線圖5（b）為例，根據(jù)2.118m／s速度等值線，可將平面分為3個(gè)區(qū)域.當(dāng)（E1，E2）取區(qū)域I樣本值時(shí)（頻次22 658，占比12%），速度區(qū)間范圍為2.113～2.118m／s；當(dāng)（E1，E2）取區(qū)域II樣本值時(shí)（頻次80 235，占比42%），速度區(qū)間范圍為2.118～2.3m／s；當(dāng)（E1，E2）取區(qū)域 III 樣本值時(shí)（頻次87 079，占比45%），速度區(qū)間范圍為0.98～2.113 m／s，雖然區(qū)間III的樣本頻次較高，但是其速度的范圍較大，所以圖4（b）中出現(xiàn)了較大速度頻次高的現(xiàn)象，泥沙體積分?jǐn)?shù)和底床形變的頻率分布圖可以使用同樣的方法分析予以分析，不再贅述.

圖4 異重流水力參數(shù)的頻次統(tǒng)計(jì)（x＝100m）Fig.4 Frequency statistics of the hydraulic parameters（x＝100m）

圖5 E1－E2平面上計(jì)算異重流水力參數(shù)的等值線圖（x＝100m）Fig.5 Equivalent line of the hydraulic parameters on the E1－E2plane（x＝100m）

圖6 E1－E2 平面ew 的等值線與（E1，E2）樣本分布圖（Ri＝0.3）Fig.6 Equivalent line of ewon E1－E2plane and （E1，E2）sample distribution（Ri＝0.3）

將概率最大（E1，E2）取值輸入數(shù)學(xué)模型，計(jì)算所得各斷面水力參數(shù)分別記為hr、ur、φr和 Δzbr（以下簡(jiǎn)稱(chēng)標(biāo)準(zhǔn)值）.將（E1，E2）樣本值輸入數(shù)學(xué)模型計(jì)算所得水力參數(shù)（h、u、φ、Δzb）分別除以各斷面水力參數(shù)（hr、ur、φr、Δzbr），比較該倍數(shù)的變化情況，可以定量分析經(jīng)驗(yàn)系數(shù)不確定性對(duì)數(shù)學(xué)模型計(jì)算結(jié)果的影響.圖7給出了4個(gè)典型斷面（100，200，500，1 000m）上h／hr、u／ur、φ／φr、Δzb／Δzbr的范圍，并通過(guò)顏色的深淺區(qū)分了25%、50%、75%、95%樣本對(duì)應(yīng)的該比值的范圍.圖中黑色虛線代表標(biāo)準(zhǔn)值隨空間的變化.

圖7 不同比例（E1，E2）樣本計(jì)算的異重流水力參數(shù)在4個(gè)斷面的區(qū)間范圍Fig.7 Interval of turbidity current hydraulic parameters calculated by different proportions of（E1，E2）

由圖7可知，隨著距進(jìn)口距離的增大，計(jì)算厚度和速度的標(biāo)準(zhǔn)值迅速增大、泥沙體積分?jǐn)?shù)標(biāo)準(zhǔn)值迅速減小、底床形變標(biāo)準(zhǔn)值變化不大.95%（E1，E2）樣本計(jì)算得到的比值（h／hr、u／ur、φ／φr、Δzb／Δzbr）范圍也隨著距離增大而增大.這說(shuō)明，（E1，E2）不確定性對(duì)計(jì)算水力參數(shù)的影響隨著與進(jìn)口邊界距離增大而增大.在x＝500m處，95%樣本計(jì)算的厚度是標(biāo)準(zhǔn)值的0.06～5.35倍、計(jì)算速度是標(biāo)準(zhǔn)值的0.49～1.24倍、計(jì)算泥沙體積分?jǐn)?shù)是標(biāo)準(zhǔn)值的0.37～4.35倍，計(jì)算底床形變是標(biāo)準(zhǔn)值的0.1～1.11倍；到x＝1 000m處，計(jì)算水力參數(shù)與標(biāo)準(zhǔn)值之間倍數(shù)擴(kuò)大至0.03～5.73（厚度）、0.46～1.23（速度）、0.39～5.31（泥沙體積分?jǐn)?shù)）和0.01～1.08（底床形變）.由圖7可見(jiàn)，隨著距離的增加，95%樣本計(jì)算所得水力參數(shù)的范圍遠(yuǎn)大于25%樣本計(jì)算所得水力參數(shù)范圍，這說(shuō)明異重流數(shù)學(xué)模型對(duì)水卷吸經(jīng)驗(yàn)式中經(jīng)驗(yàn)系數(shù)（E1，E2）的取值十分敏感.比如，在x＝500m處，95%樣本計(jì)算的厚度是標(biāo)準(zhǔn)值的0.1～3.16倍，比25%樣本的范圍擴(kuò)大了72%，計(jì)算速度、泥沙體積分?jǐn)?shù)和床面形變也有類(lèi)似現(xiàn)象.通過(guò)上述分析可以看出：經(jīng)驗(yàn)系數(shù)（E1，E2）對(duì)數(shù)學(xué)模型的計(jì)算結(jié)果影響很大，且隨著距離的增加，經(jīng)驗(yàn)系數(shù)不確定性對(duì)模型的影響在逐漸增大；概率最大（E1，E2）值代入數(shù)學(xué)模型可能低估異重流厚度和泥沙體積分?jǐn)?shù)，高估異重流速度和床面形變.

4 結(jié)論

（1）本文應(yīng)用基于貝葉斯定理的Gibbs隨機(jī)采樣法得到 概率最大 系數(shù)組 合 （E1，E2）＝ （2.38，0.52），與原公式的取值（E1，E2）＝（2.4，0.5）相近，該方法可以作為最小二乘法等傳統(tǒng)經(jīng)驗(yàn)式擬合方法的重要補(bǔ)充.

（2）比較（E1，E2）不同比例樣本計(jì)算結(jié)果在不同斷面的區(qū)間范圍，可以看出：相同比例樣本輸入模型計(jì)算得到異重流水力參數(shù)區(qū)間范圍隨運(yùn)動(dòng)距離的增加逐漸增大；經(jīng)驗(yàn)系數(shù)不確定性對(duì)數(shù)學(xué)模型的計(jì)算結(jié)果影響逐漸增大；采用概率最大經(jīng)驗(yàn)系數(shù)值可能低估異重流厚度和泥沙體積分?jǐn)?shù)，高估速度和底床形變.

（3）比較不同比例（E1，E2）樣本輸入數(shù)學(xué)模型得到的同一斷面計(jì)算水力參數(shù)區(qū)間范圍可以發(fā)現(xiàn)，在x＝1 000m 處，95% （E1，E2）樣本計(jì)算得到的厚度范圍比25%樣本計(jì)算所得厚度范圍擴(kuò)大了83%，說(shuō)明模型對(duì)（E1，E2）系數(shù)組合的取值十分敏感.

由于概率最大經(jīng)驗(yàn)系數(shù)取值輸入模型預(yù)測(cè)的異重流水力參數(shù)可能出現(xiàn)偏差，所以下一步將在兩個(gè)方面開(kāi)展工作：① 增加實(shí)驗(yàn)數(shù)據(jù)獲取，以便對(duì)經(jīng)驗(yàn)系數(shù)的取值進(jìn)行更加準(zhǔn)確的率定；② 開(kāi)展不同水卷吸經(jīng)驗(yàn)式不確定性對(duì)數(shù)學(xué)模型影響的比較，為經(jīng)驗(yàn)式選擇提供依據(jù).