金沙江白格堰塞湖潰決洪水預(yù)報誤差與改進(jìn)

王敏　盧金友　姚仕明　朱勇輝　周建銀

摘要：金沙江“11.3”白格堰塞體潰決造成下游河道數(shù)百千米內(nèi)的超常規(guī)洪水過程，對其進(jìn)行準(zhǔn)確模擬是洪水預(yù)報和制定應(yīng)對措施的基礎(chǔ)。采用MIKE11和筆者開發(fā)的一維潰壩洪水模型（DBFM）對壩下洪水演進(jìn)進(jìn)行了對比模擬，兩個模型均為求解圣維南方程組，其中MIKE11采用Abbott六點(diǎn)隱式差分方法，DBFM模型采用基于有限體積的HLL+MUSCL界面插值方法顯式求解。結(jié)果表明：與實(shí)測數(shù)據(jù)相比，兩個模型模擬結(jié)果在洪水傳播時間和洪峰衰減上表現(xiàn)出顯著差異。隨后分析了河道基流、動量方程慣性、河道糙率、河道地形等因素對洪水模擬誤差的影響，論證了采用滾動預(yù)報方法來提高洪水預(yù)報精度的可行性。最后對山區(qū)河流堰塞湖潰決洪水模擬的難點(diǎn)和要求提出了展望。

關(guān)鍵詞：堰塞湖;洪水演進(jìn);數(shù)值模擬;滾動預(yù)報;金沙江

中圖法分類號：P33

文獻(xiàn)標(biāo)志碼：A

DOI：l0.16232/j.cnki.1001-4179.2019.03.007

文章編號：1001-4179（2019）03-0034-06

2018年10月11日和11月3日，西藏江達(dá)縣波羅鄉(xiāng)白格村川藏交界的金沙江上游河道右岸同一位置兩次發(fā)生山體滑坡形成堰塞湖。其中“11·3”堰塞壩是在“10·11”堰塞湖潰決殘余壩體基礎(chǔ)上形成的，其壩高較“10·11”堰塞壩還高約30m。為降低蓄洪量、減小洪水風(fēng)險，采取了人工挖槽泄流的措施。11月12日，泄流槽進(jìn)水，11月13日13：45，壩前水位達(dá)到最高點(diǎn)2956.40m。

為提前預(yù)估堰塞湖下游金沙江沿岸洪水風(fēng)險，給居民轉(zhuǎn)移、工程處置水庫調(diào)度等應(yīng)急決策提供科學(xué)依據(jù)，堰塞湖形成后，在水利部長江水利委員會組織下，長江科學(xué)院相關(guān)技術(shù)人員采用成熟的商業(yè)軟件MIKE11開展了洪水演進(jìn)計算。堰塞湖潰決后，又根據(jù)沿程各站實(shí)測流量過程，開展?jié)L動預(yù)報計算。應(yīng)急決策過程中，因時間倉促，資料不足，模擬結(jié)果難免存在一定的誤差。

潰壩洪水的模擬比常規(guī)水流模擬困難，其原因在于：①山區(qū)河道地形比降大、起伏劇烈;②流量在時空上變化極為劇烈，存在干濕變化，間斷特征顯著。此外，山區(qū)河流地形水文資料缺少、強(qiáng)非恒定非均勻流的阻力仍未明晰等原因，均可能造成潰壩洪水模擬的誤差。由于堰塞湖不太常發(fā)生且實(shí)測資料較少，模擬實(shí)際堰塞湖洪水的成果不多。2008年唐家山堰塞湖發(fā)生后，眾多學(xué)者進(jìn)行了堰塞湖潰決洪水演進(jìn)的模擬。劉帆采用具有人工黏性的修正MacCormack數(shù)值格式計算了潰壩洪水的一維演進(jìn)張細(xì)兵、朱勇輝、黃明海也計算了唐家山堰塞湖的潰決洪水及其演進(jìn)，采用的是MIKE11的Abbott六點(diǎn)隱格式和傳統(tǒng)的四點(diǎn)偏心格式求解圣維南方程組。

模擬一維水流運(yùn)動的圣維南方程組是一個雙曲型方程組。為準(zhǔn)確模擬間斷，差分格式必須滿足守恒性。其中，有限體積方法天然滿足守恒性，適合用來構(gòu)建潰壩洪水模擬的模型。大部分潰壩模型隱含矩形河道假設(shè)，不適應(yīng)于山區(qū)河道。張大偉建立了采用Z和Q為求解變量的一維有限體積潰壩模型，采用HLL格式的近似黎曼解求解界面通量并采用MUSCL方法對界面兩側(cè)物理量進(jìn)行了重構(gòu)，宋利祥以此為基礎(chǔ)建立了河網(wǎng)模型，鄭川東進(jìn)一步將水面梯度項，放入對流項中，以更好地滿足動量守恒。這些模型理論上適合應(yīng)用于實(shí)際山區(qū)河道的潰壩洪水，但是尚待實(shí)踐檢驗(yàn)。向小華采用特征線方法求解了以A和Q為變量的非守恒形式的圣維南方程組，構(gòu)建了一維河網(wǎng)模型。以守恒變量A、Q為基本變量，求解守恒型方程，更符合方程守恒性的本質(zhì)，也是準(zhǔn)確捕捉間斷的要求。

本文采用MIKE11模型和筆者建立的基于有限體積法的一維潰壩洪水模型（DBFM），對金沙江“11·3”白格堰塞湖潰決洪水演進(jìn)過程進(jìn)行對比計算;通過模擬不同的參數(shù)和工況，分析模型、地形、水流狀況等因素對洪水模擬預(yù)報誤差的影響;論證了滾動預(yù)報對提高模擬預(yù)報精度的效果。

1 模型介紹

1.1 MIKE11模型

MIKE11模型求解一維圣維南方程組，其形式如下：

公式

式中，t，x為時間和一維坐標(biāo);A，z，，R，Q分別為過水?dāng)嗝婷娣e、水位、水力半徑以及斷面流量;Qi為單位長度內(nèi)的旁側(cè)入流量;β為動量方程慣性項修正因子。MIKE11模型默認(rèn)慣性項修正因子取值為

公式

MIKE11模型采用Abbott六點(diǎn)隱格式進(jìn)行離散，離散所得的三對角方程采用追趕法進(jìn)行求解。

1.2 一維潰壩洪水模型

為模擬山區(qū)復(fù)雜條件下可能存在間斷的潰壩洪水，筆者構(gòu)建了一維潰壩洪水模型（DBFM）。模型采用有限體積法離散圣維南方程組。采用HLL格式的近似Riemann解對界面進(jìn)行插值。界面通量的計算為

公式

式中，V為求解變量，S為傳播流速（此處為流速），F(xiàn)（V）為界面通量，下標(biāo)L、R代表界面的左、右。為提高界面插值的精度，對界面兩側(cè)的變量，采用MUSCL方法對界面兩側(cè)變量進(jìn)行數(shù)值重構(gòu)。以第i個控制體的右界面（i+1/2）為例，界面兩側(cè)的通量為

公式

式中，公式是限制器。本次研究所用模型采用VanLeer限制器：

公式

參考Ying對水位梯度源項采用順風(fēng)差分以保證穩(wěn)定性。有限體積模型在出入口各需兩個條件。Zhang在模型出入口，除了，上游給流量、下游給水位之外，還各補(bǔ)充了一個水位流量關(guān)系。該模型考慮到出口水位無法給定，除給定入口流量和出口水位流量關(guān)系之外，入口補(bǔ)充條件為根據(jù)連續(xù)方程離散得到的新時刻水位和流量的關(guān)系，出口補(bǔ)充條件為流量沿程梯度為零。

2 實(shí)測洪水過程模擬對比

堰塞湖模擬計算時，以實(shí)際潰口流量過程為入口條件。潰口流量過程由壩前水位過程與堰塞湖水位-庫容曲線推算。河道糙率取0.05，模擬范圍從潰口至其下游河道600km。模擬所用河道地形由DEM地形數(shù)據(jù)提取。沿程主要站點(diǎn)的斷面形態(tài)如圖1所示。可見，DEM無法獲取水下地形高程，故計算采用的斷面地形有一定的概化，可能與實(shí)際河槽有一定的區(qū)別，這將不可避免地帶來模擬誤差。

圖2是各站實(shí)測流量過程。需要指出的是，由于洪水陡漲陡落，現(xiàn)場水位的測量較為困難，有時不得不借助臨時測站，而本次奔子欄以上各站的洪峰流量超過了有記錄歷史最大洪水，這使得按照水位流量關(guān)系推算的流量也存在一定的誤差，加上非恒定流的水位流量關(guān)系本身也存在爭議，因此實(shí)測流量和水位本身也可能存在一定的誤差。

實(shí)際洪水并非發(fā)生在干河道上，河道本身有一定的基流。根據(jù)各站洪水起漲前的流量資料，可以得到沿程的河道基流情況。其中有些來自支流入?yún)R，因缺少支流流量資料，本次統(tǒng)一按照區(qū)間均勻入?yún)R考慮。各站基流流量如表1所示。

計算和實(shí)測各站洪峰過程過程對比如圖3所示。由圖可見：與實(shí)測相比，蘇洼龍以。上，DBFM模型比MIKE11模型符合性更好，而蘇洼龍以下則相反。DBFM模型洪峰衰減較快，洪峰傳播時間符合較好，而MIKE11模型計算洪峰衰減較慢，但是洪峰傳播速度較快。預(yù)報時若能結(jié)合兩者的模擬結(jié)果，可以提高預(yù)報的精度。從蘇洼龍壩上和壩下的實(shí)測流量過程差異可見，蘇洼龍電站攔蓄了部分洪水，顯著改變了洪水過程。由于模型并未考慮蘇洼龍的攔蓄作用，因此，蘇洼龍以下模擬結(jié)果與實(shí)測的符合是一種“巧合”。

3 誤差分析

3.1 河道基流的影響

為分析河道基流的影響，模型模擬的條件除了潰壩下游初始河床為干河床，其他與上一節(jié)一致。圖4是潰口下游2個測站的實(shí)測流量過程與計算流量過程的對比。與圖3對比可見，干河床條件下洪水傳播速度變慢，洪峰流量減小。

3.2 動量方程慣性項的作用

模型計算表明，此次潰壩洪水整個過程均為緩流，且Fr數(shù)絕大多數(shù)在0.5以下。但是，潰壩洪水本身可能出現(xiàn)急流，因此，此處仍對慣性項的作用予以分析討論。MIKE11模型采用的是Abbott六點(diǎn)隱格式，具有很強(qiáng)的適應(yīng)能力，但是，對于急流仍有局限。為了克服急流模擬的困難，引入了LPI（Local Particial Intia）因子。該因子在緩流時不起作用。DBFM模型可以正常模擬紊流情況，不需要對動量方程慣性項加以限制。

圖5對比了慣性項正常計算和歸零時各站的洪水過程。

由圖5可見，慣性項歸零后，洪水過程更為“尖瘦”，但作用有限。由于動量方程慣性項對流量過程實(shí)際起著擴(kuò)散坦化的作用，在潰壩洪水模擬中，若出現(xiàn)急流，由于慣性項被抑制，MIKE11模型可能銳化洪水過程、增大洪峰流量。

3.3 河道糙率的影響

河道糙率對洪水模擬的精度也有一定的影響。糙率越大，阻力越大，同流量下水位越高槽蓄量越大、流速越小。圖6對比了糙率取0.05和0.06兩種情況下葉巴灘和巴塘站洪水過程。可見，糙率增大后，洪水過程被推遲、洪峰流量減小。

3.4 其他影響因素分析

前面分析表明，河道基流、糙率、動量方程慣性項的處理均會對潰壩洪水模擬有一定的影響。因缺少足夠的資料，蘇洼龍圍堰攔蓄洪水過程未能在模擬中考慮，這是造成蘇洼龍及其下游洪水模擬誤差的重要原因。從基流的作用可見，支流入?yún)R對洪水演進(jìn)也有影響，特別是支流來流量較大的情況。

本次洪水模擬誤差最大的影響因素是地形誤差。由于缺少水下地形資料，河槽變得比實(shí)際寬淺得多。一方面，同樣水深情況下，過水面積要大得多，這將導(dǎo)致洪水下泄時流量過程的坦化;另一方面，缺少水下地形相當(dāng)于抬高了河底，由于沿程并非均勻抬高，則可能增大河床縱比降，加快洪水運(yùn)動。其綜合效果，與具體河道情況有關(guān)。

4 滾動預(yù)報

由于本次潰壩洪水下游傳播距離長達(dá)600km以上，長距離計算必然導(dǎo)致誤差的累積。潰口的流量過程本身也存在不小的誤差，尤其是在峰值流量上。本次潰壩洪水模擬預(yù)報時，由于蘇洼龍圍堰的攔蓄未能在模擬中加以考慮，蘇洼龍下游洪水過程的模擬預(yù)測更難以實(shí)現(xiàn)。在此情況下，利用潰決過程中沿程多個測站的實(shí)測洪水資料，不斷校正模型參數(shù)，進(jìn)行滾動模擬，就可以成為減小模擬誤差、提高預(yù)報精度的有效手段。此處仍采用DBFM模型，以奔子欄水文站為入口起算下游洪水演進(jìn)過程。圖7是石鼓斷面計算與實(shí)測流量過程的對比。可見，計算結(jié)果在峰值和相位上，均與實(shí)測符合較好，比圖3從潰口開始推算的石鼓流量過程大為改善。

5 結(jié)語

本文采用兩種基于不同求解方法的數(shù)值模型進(jìn)行對比計算，分析了河道基流、糙率、動量方程慣性項、河道地形等因素對模擬誤差的影響，論證了采用滾動預(yù)報的方法可以提高模擬精度。

（1） 山區(qū)河流河床大比降、高起伏，潰壩洪水過程強(qiáng)非恒定，且可能出現(xiàn)急流或干濕間斷，這些復(fù)雜條件對數(shù)值模型的模擬能力提出了挑戰(zhàn)。

（2） 河道基流、糙率、動量方程慣性項、河道地形均會對潰壩洪水模擬有一定的影響，其中地形的影響相對較大。基于潰決過程中實(shí)時監(jiān)測數(shù)據(jù)的滾動預(yù)報方法可提高洪水模擬和預(yù)報的精度。

（3） 實(shí)際預(yù)報中，綜合采用多個模型對比預(yù)測，能夠提高洪水預(yù)報的精度和把握，更好地為應(yīng)急決策提供支撐。

（4） 山區(qū)河流資料缺乏，使得潰壩洪水模擬預(yù)報存在一定程度的不確定性。快速、精確地獲得山區(qū)河流河道地形數(shù)據(jù)，是山區(qū)河流潰壩洪水模擬和預(yù)報迫切需要解決的問題。

參考文獻(xiàn)：

[1]劉帆.一維潰壩模擬及四川堰塞湖應(yīng)急計算[D].北京：清華大學(xué)，2008.

[2]王光謙，付旭東，李鐵建，等.汶川地震災(zāi)區(qū)堰塞湖應(yīng)急處置中的計算分析[C]//中國科協(xié)防災(zāi)減災(zāi)論壇，2008.

[3]張細(xì)兵，盧金友，范北林，等.唐家山堰塞湖潰壩洪水演進(jìn)及下泄過程復(fù)演[J].人民長江，2008，39（22）：76-78.

[4]朱勇輝，范北林，盧金友，等.唐家山堰塞湖潰壩洪水分析及泄流沖刷模擬[J].人民長江，2008，39（22）：79-82.

[5]黃明海，金峰，楊文俊.不同河道狀況下的堰塞壩潰壩洪水特性研究[J].人民長江，2008，39（22）：66-68，82.

[6]Michel-Dansac V，Berthon C，Clain S，et al. A well-balanced

scheme for the shallow-water equations with topography or Manning friction[J]. Journal of Computational Physics，2017（335）：115-154.

[7]Valiani A，Caleffi V. Momentum balance in the shallow water equations on bottom discontinuities[J]. Advances in Water Resources，2017（100）：1-13.

[8]LeFloch P C，Thanh M D. A Godunov-type method for the shallow

water equations with discontinuous topography in the resonant regime[J].Jourmal of Computational Physics，2011，230（20）：7631-7660.

[9]張大偉，程曉陶，黃金池，等.復(fù)雜明渠水流運(yùn)動的高適用性數(shù)學(xué)模型[J].水利學(xué)報，2010（5）：531-536.

[10]宋利祥，李清清，胡曉張，等.基于有限體積法的河網(wǎng)水動力并行計算模型研究[J].長江科學(xué)院院報，2018（14）：1-7.

[11]鄭川東，白鳳朋，楊中華.求解守恒形式的圣維南方程中處理不規(guī)則斷面的一種改進(jìn)方法[J].水電能源科學(xué)，2017（12）：95-99，103.

[12]向小華，吳曉玲，牛帥，等.通量差分裂方法在一維河網(wǎng)模型中的應(yīng)用[J].水科學(xué)進(jìn)展，2013（6）：894-900.

[13]DHI. MIKE 1D Reference manual]EB/OL].2017.

[14]Ying X，Khan A A， WangS S Y.Upwind Conservative Scheme for the Saint Venant Equations[J]. Journal of Hydraulic Engineering，2004，130（10）：977-987.

[15]Zhang s，Duan J G. 1D finite volume model of unsteady flow over mobile bed[J]. Jourmal of Hydrology，2011，405（1）：57-68.

引用本文：王敏，盧金友，姚仕明，朱勇輝，周建銀.金沙江白格堰塞湖潰決洪水預(yù)報誤差與改進(jìn)[J].人民長江，2019，50（3）：34-39.

Study on dam break flood forecast error and improvement for Baige barrier lake on Jinsha River

WANG Min，LU Jinyou，YAO Shiming，ZHU Yongwei，ZHOU Jianyin

（Changjiang River Scientific Research Institute，Changjiang Water Resource Commission，Wuhan 430010，China）

Abstract：The broken of "11.3" Baige barrier lake on Jinsha River resulted in the ultra-normal flood in hundreds of kilometers of the downstream. Accurate simulation of the flood is the basis of the flood forecast and the development of the response measures. The flood routing law in the downstream of the dam was simulated by MIKE11 and a one-dimensional dam break flood model developed by the author. The two models both solve the Saint-Venant equations，while the MIKE1 1 uses the Abbott six-point implicit difference method and the DBFM model uses a finite volume of the HLL + MUSCL interface interpolation method. The simulation results showed that the two models show differences in the flood propagation and the peak attenuation comparing with the measured data. Then，the influence of the factors such as the channel base flow，the momentum equation inertia，the channel roughness and the channel terrain on the simulation error of the flood were further analyzed，and the feasibility of using rolling forecast method to improve the flood forecast accuracy was demonstrated. In the end，the difficulties and requirements of the dam break flood simulation of the barrier lake in the mountainous area were put forward.

Key words：barrier lake;flood routing;numerical simulation;rolling forecast;Jinsha River