射月溝水庫潰壩洪水模擬及潰壩原因分析

趙懸濤，劉昌軍，王文川，顧斌杰，楊 昆，張 淼，馬建明

（1.中國水利水電科學研究院防洪抗旱減災中心，北京 100038；2.華北水利水電大學水資源學院，鄭州450046；

3.水利部水旱災害防御司，北京 100053）

0 引 言

我國中小型水庫多且分布廣，水庫大壩一旦發生潰決，會造成嚴重的人員傷亡和經濟財產損失，已成為國家防災減災中一個亟待解決的問題［1］。潰壩洪水流量極大，流速快，沖擊力強，嚴重威脅了水庫下游人民的生命財產安全。因此，對水庫潰壩原因及過程進行系統分析和總結，為國家防災減災工作提供案例支撐是非常有必要的。對于潰壩機理、潰壩過程及潰壩洪水演進等，國內外多學者開展了大量研究和案例分析工作。蔣先剛［2］、牛志攀［3］、李云［4］、劉青泉［5］等學者對潰口的拓展過程作了詳細探討，王笑［6］、劉慧玲［7］、Hou［8］、張大偉［9］等學者在潰壩水流的模擬方面作了深入研究。Wu［10］提出的潰壩模型可用于模擬漫頂和管涌引起的潰壩，并采用物理實驗與實際案例的50 組數據進行了對比，在潰壩峰值流量、潰口寬度上都與觀測值較為接近；Wang 等［11］提出一種多線截面圖法，研究了在不規則斷面上潰壩波的解析解。ISSAKHOV 等［12］采用流體體積法（VOF，volume of fluid）模擬了壩體破裂時水面的運動，模擬結果與實驗數據基本吻合。馬利平等［13］采用源項法將潰口演變模型和二維水動力模型進行耦合，并通過對一個土石壩和兩個堰塞壩潰決算例的模擬，驗證了耦合模型的模擬精度。

然而，針對無資料小流域暴雨洪水、水庫潰壩過程及災害等進行全面調查和系統分析的研究較少。根據射月溝水庫“7·31”潰壩洪災現場的調查資料，采用劉昌軍等［14］提出的時空變源水文模型對“7·31”暴雨和入庫洪水過程進行計算分析。并結合現場實測和視頻資料，分析了射月溝水庫溢洪道泄洪、漫頂溢洪、潰壩洪水過程及下游淹沒和災害損失情況，計算得到“7·31”暴雨洪水過程中射月溝水庫入庫和出庫的洪水流量過程線。最后基于中國水科院自主研發的洪水分析平臺IFMS［15］構建二維水動力學模型對射月溝水庫潰壩水流的演進進行數值模擬，系統分析了潰壩洪水過程、下游淹沒情況和災害形成原因。

1 資料與方法

1.1 研究區概況

射月溝流域位于東天山末端哈爾里克山南坡，西鄰大天生圈流域，南抵哈密盆地，流域分水嶺北側是伊吾縣鹽池牧場。流域面積為406 km2，出山口以上山區面積為187 km2，流域平均坡降為35‰。徑流模數為0.94，年均降雨量150 mm 左右，多年平均徑流量1 581 萬m3。流域內有6 條支溝，河谷類型為寬淺的“U”型谷，區內地貌大部分為侵蝕構造的中低山，山頂山坡基巖裸露、基本無植被，流域內氣象站點、小流域及水系分布如圖1所示。

圖1 研究區域概況Fig.1 Overview of the research area

射月溝水庫位于哈密市伊州區沁城鄉東南約8 km，水庫壩址處地理位置為東經94°37′，北緯42°46′，是一座以工業供水、灌溉為主，兼顧防洪、生態等綜合利用的水利樞紐工程。該水庫為攔河式水利工程，工程規模為小（Ⅰ）型，大壩為瀝青混凝土心墻砂礫石壩［16］，壩體剖面為三角形，壩長403 m，壩頂高程為1 496.65 m，水庫總庫容677.9 萬m3，水庫下游沿河村落主要為二宮村，總人口816人，常住人口201人。

1.2 數據資料

用于構建水文和水動力學模型的數字化資料主要包括地形髙程數據、小流域、水系、節點、土壤質地及土地利用等基礎數據。其中高程信息數據由現場三維激光掃描得到激光點云數據和從地理空間數據云網站（http：//www.gscloud.cn/home#page1/4）下載的30 m 分辨率數字髙程數據融合而成。流域邊界、水系、節點、土壤質地及土地利用數據來自全國山洪災害調查評價成果數據，水庫設計基本情況和潰壩過程通過搜集水庫設計資料和現場調查獲取，降雨數據采用氣象部門共享的衛星、氣象站、雷達融合數據。表1 記錄了數據的詳細信息，圖2為庫區三維激光點云散點圖。

圖2 庫區三維激光點云Fig.2 3D laser point cloud of reservoir area

表1 基礎數據Tab.1 Basic data

1.3 研究方法

（1）模塊化小流域洪水模擬系統。模塊化小流域洪水模擬系統（FFMS，Flash Flood Module Simulation System）［17］主要特點包括：①軟件支持多個國內外水文模型自動建模和手動拖拽式建模，包括HEC、PRMS、新安江及時空變源分布式水文模型。②軟件集成了不同流域下墊面產匯流和河道參數庫，大大提高了建模效率。③自動提取地貌響應單元，自動對應不同產匯流模型，實現智能化建模。④支持不同模型要素、ARCGIS 工作底圖和計算結果可視化。⑤具體小流域劃分、參數提取、面雨量插值和自動率定等功能。FFMS軟件界面見圖3。

圖3 模塊化小流域洪水分析系統主界面Fig.3 The main interface of the modular small basin flood analysis system

入庫洪水過程推算采用的是時空變源分布式水文模型，如何高效穩定的求解非飽和土壤下滲過程是該模型要解決的核心問題。此模型的處理方式是基于Talbot 和Ogden［18］下滲計算方法，利用Green-Ampt 模型的下滲速率，結合Brooks- Corey［19］模型非飽和導水率曲線及Van-Genuchten 模型［20］土壤水分特征曲線計算方法，得到離散區間的含水量及濕潤鋒的計算公式，推導過程參見文獻［14］。

（2）洪水分析系統。洪水分析系統（IFMS，Integrated Flood Modeling System）由中國水科院聯合南京水科院、河海大學、山東大學及大連理工大學等單位開發，主要特點包括：①先進的內核計算引擎：能用于大型河網、山區陡坡河道、管網水流模擬，激波捕捉高性能二維格式，強大的一二維耦合功能用于漫堤/潰堤洪水模擬以及城市地表與地下管網水流交互。②面向對象的建模思路和數據管理模式，極大提高了基礎數據的復用率和建模的便捷性。③強大的輔助建模工具：提供網格剖分、管網拓撲處理、土地利用提取等功能強大的輔助工具。軟件界面如圖4所示。

圖4 洪水分析系統界面Fig.4 Flood Analysis System Interface

本文基于IFMS 平臺建立的二維水動力學模型采用Godunov算法［9］進行數值計算。其中黎曼問題采用Roe格式［21］的近似黎曼解進行計算，底坡源項采用特征分級離散以保證模型的守恒性，阻力源項采用隱式離散用于提高模型的穩定性。

（3）水力學公式。基于庫區激光點云數據（圖2）和現場潰壩視頻資料（圖7），計算得到不同時刻庫區的庫水位，采用薄壁堰公式對溢洪道泄水和大壩漫頂溢流過程進行計算。

式中：m0為流量系數，取值為0.424 767；b為堰頂寬度，大壩為395 m，溢洪道為30 m；g為重力加速度，m/s2；H為堰前水深，m。

潰口洪峰流量的推算采用鐵道部科學研究院的經驗公式。

式中：Qmax為潰壩最大流量，m3/s；L為庫區長度，取為1 100 m；當L/B＞5 時，取L/B=5；b為潰口平均寬度，取為48.71 m；h0為壩前水深，取為18.63 m；h為潰口處殘留壩體的平均高度，m；k為經驗系數，k和h取為0。

潰壩流量過程與潰壩最大流量、下游水位及潰壩可泄庫容有關，計算采用公式（3）。

式中：Qt為t時刻的流量，m3/s；Qm為潰壩最大流量，m3/s；T為潰壩歷時。

2 暴雨洪水分析

2.1 暴雨特征分析

2018年7月29日夜間至8月2日上午，新疆哈密市北部山區出現暖濕背景條件下的強降水天氣過程，尤其是7月30日夜間到31日。收集了射月溝水庫及周邊各氣象站、衛星降水融合數據、衛星遙感等相關資料，分析得到射月溝水庫“7·31”暴雨特征如下：

（1）持續時間長：此次暖區強降水對流云團始于7月31日凌晨2時，從哈密市東部平原區發展演變北上，強降水對流云團持續覆蓋在以射月溝流域水庫上游區為中心的區域，其上游區沁城鄉小堡自動氣象站降水時間長達10 h（圖5），于上午12 時減弱北上蒙古國出境。

（2）影響區域大：此次強降水對流云團影響覆蓋哈密市伊州區東北區域、伊吾縣大部，覆蓋最大面積可達近4 萬km2，約占哈密市近1/3 的面積，射月溝流域及其以東區域為強降水中心覆蓋區域，累計降雨量在55 mm 以上的面積占水庫上游區域面積的84%。

（3）降水強度強：在射月溝流域的沁城鄉小堡、射月溝水庫、沁城鄉三個自動氣象站（站點位置分布見圖1）分別觀測到115、100、79 mm 的強降水，均突破有記錄以來極值。沁城鄉水庫上游小堡站2018年7月31日凌晨5 時至11 時6 h 降雨量為97.9 mm，超過日降雨量歷史極值52.4 mm。氣象站點逐時降水量分布如圖5 所示。各時段降雨量均遠超有記錄以來歷史極值，為短歷時極端強降雨。

圖5 射月溝流域內氣象站降水情況Fig.5 Precipitation at weather station in Sheyuegou watershed

2.2 入庫洪水過程計算

根據國家氣象局衛星、雷達融合降雨數據和地面站實測降雨數據，可得到射月溝流域“7·31”暴雨全過程降雨數據。利用全國山洪災害調查評價成果，獲取射月溝水庫所屬流域的DEM、土壤質地、土地利用類型等模型輸入資料。采用FFMS對射月溝流域“7·31”暴雨洪水過程的模擬結果如圖6所示。

圖6 射月溝水庫7·31入庫洪水過程線Fig.6 The flood process line of Sheyuegou Reservoir on July 31

水文模型的產流模塊采用時空變源混合產流模型，坡面和河道匯流模塊采用運動波方法，降雨數據采用氣象局提供的融合降水產品。計算得到射月溝水庫的入庫洪峰流量為1 915 m3/s，時間為7月31日9點，新疆自治區水文局后期通過洪水調查推算入庫洪峰流量為1 848 m3/s，FFMS 計算結果與自治區水文局調查計算結果相比，洪峰流量計算誤差為3.5%，模擬結果較好。

2.3 出庫洪水過程計算

射月溝水庫最大庫容為678 萬m3，溢洪道最大流量為380 m3/s，壩頂高程為1 496.65 m，溢洪道堰頂高程為1 492.53 m。7月31日2∶00，射月溝流域開始降雨，此時水庫庫容為406 萬m3。5∶00，降雨強度逐漸增大。6∶25，庫水位達到了1 491.4 m，入庫洪峰流量為30 m3/s，溢洪道尚未溢洪。7∶12，庫水位到達1 492.5 m，入庫洪峰流量為476 m3/s，溢洪道開始溢洪。7∶50，溢洪道水深達到了2 m。9∶10，射月溝水庫開始漫頂溢流。10∶06，左壩壩肩開始形成潰口，最終導致潰壩。11∶00，水庫基本泄空。圖7為水庫漲水-漫壩-潰壩過程中幾個關鍵時間節點的視頻圖像資料。

圖7 不同時刻潰壩現場照片Fig.7 Photos of the dam breach at different moments

水庫初始庫容406 萬m3，水庫實時進水量由入庫洪水過程線推算得到，由水位庫容曲線實時反推庫水位變化。溢洪道泄水和漫頂溢流采用式（1）計算，由視頻資料目估壩體潰決開始于10：08，潰壩洪峰流量由式（2）計算得到，潰壩歷時由式（3）估算，11∶00 水庫泄空，將幾個關鍵時間節點作為控制點，編制自動計算程序，綜合計算得到整個過程水庫出流過程線如圖8 所示，將其作為二維水動力模型的輸入。

圖8 射月溝水庫出流過程線Fig.8 Outflow process line of Sheyuegou Reservoir

3 潰壩洪水演進模擬

3.1 模型構建

（1）計算模型構建。研究區域面積為133.07 km2。利用三維激光掃描儀獲取庫區1 m 精度的高程數據，水庫下游地形采用精度為30 m 的DEM 數據，對精度不同的兩類數據進行融合處理，得到本次潰壩水流演進數值模擬的地形數據。

利用融合后的地形資料對研究區域進行非結構化網格剖分，對河道、村莊等關鍵位置進行了加密，地形信息儲存在網格節點上，網格尺寸從1～30 m 不等，其中壩址附近網格尺寸限制在2 m 以內，研究區域內的河道網格尺寸限制在10 m 以內，其他區域網格尺寸為30 m左右不等，共剖分網格836 154個。

（2）模型邊界條件和參數設置。模型的輸入資料為地形數據、入流數據等，其中入流數據為上文推算出的水庫出流流量過程。上游邊界為入流邊界，下游邊界為自由出流的開邊界，其余邊界定義為閉合邊界。

射月溝水庫屬于山丘區水庫且水庫下游地表形態單一，不失一般性考慮，模擬區域糙率值的選取共分為三類；村莊所在地糙率值取0.05，河道取值0.04，其他位置統一取值為0.035［22］。

3.2 模擬結果分析

模型模擬了溢洪道泄洪、大壩全斷面過水及壩體潰決至水庫泄空全過程，歷時226 min，模型計算耗時15 min。模擬過程涉及四個關鍵時間節點，包括7∶12 溢洪道開始泄水、9∶08 大壩全斷面過水、10∶08 潰口形成及11 點水庫基本泄空。模擬結果可得到任意時刻模型區域內的流場、流速、水深等水力要素的分布特征。

圖9 為不同時刻洪水淹沒示意圖，從中可以看出：8∶30 時左右，洪水已經演進到二宮村；在10時左右，整個二宮村基本完全處于淹沒狀態。

為直觀反映水庫下游區域淹沒水深隨時間變化情況，選取頭宮村、出山口、二宮村三個代表點（分別距大壩2.06、4.76、6.24 km，位置分布如圖9 所示），計算得到3 個代表點水深變化過程線如圖10 所示，3 個代表點的最大淹沒水深分別為3.54、12.11、3.49 m。據新疆水文局工作組調查結果，潰壩洪峰到達下游6 km（二宮村2 隊）處產生3 米多高的涌浪；據當地群眾反映，事發過程中，二宮村2 隊河道水位在8 點至10 點30 分左右出現過明顯漲-落-急漲的變化過程，水位下落時間有數十分鐘；兩種描述均和模擬結果較為一致。出山口地形較為狹窄，為典型的埡口，導致水位急速上涌，此位置也是模擬區域的最大淹沒水深點，模擬結果和現場調查結果較為一致。

圖9 研究區域不同時刻淹沒情況Fig.9 The submergence of the study area at different times

圖10 關鍵位置淹沒水深變化過程線Fig.10 Change process line of submerged water depth at key locations

根據模型計算結果，對下游洪水淹沒區不同水深面積進行提取統計，統計結果如表2 所示。從表2 中可以看出，洪水淹沒水深在2～6 m 之間的面積占了總面積的23.02%，受災最嚴重的二宮村二隊、三隊均在這個范圍內。

表2 研究區域最大淹沒水深占比情況Fig.2 The proportion of the largest submerged water depth in the study area

4 潰壩原因分析和經驗總結

結合暴雨洪水特征、入庫洪水、潰壩過程及潰壩洪水演進和淹沒過程等多個方面對潰壩原因進行系統分析，發生漫頂潰壩和洪災主要原因如下。

（1）短歷時極端強降雨引發遠超水庫防洪能力的特大洪水是造成漫頂潰壩的決定性因素。本次暴雨洪水過程具有降雨量大、強度大和匯流快等特點，根據FFMS 軟件計算結果，射月溝水庫2 h以上就形成入庫洪峰，最大洪峰流量達到1 915 m3/s，遠超過歷史最大入庫洪峰流量170 m3/s。本次暴雨洪水入庫徑流總量達2 197 萬m3（年平均徑流總量1 819 萬m3），遠超水庫蓄洪（水庫總庫容670 萬m3）和泄洪能力（最大泄洪能力380 m3/s）。水庫水位在短時間內猛漲，7∶02 分水位1 491.8 m，尚未開始溢洪，9時左右水位就達到1 497.65 m（防浪墻頂高度），僅僅2個小時水庫水位就上漲了5.85 m，造成漫頂溢流。

（2）潰壩洪水是造成嚴重洪災的主要原因。由于潰壩流量極大，洪水流速快，沖擊力強，給水庫下游地區造成嚴重經濟損失和重大人員傷亡。根據現場調查結果和IFMS 軟件計算得到潰壩洪水演進和淹沒過程結果，潰壩洪峰經過壩下2 km處（二宮村1隊）時洪峰涌浪高達10 m，到達下游6 km（二宮村2隊）處仍有3 m高的涌浪，造成了二宮村1隊和2隊的重大洪水災害。

（3）壩下游出山口束窄加劇了潰壩洪水災害。根據IFMS軟件的計算結果，在潰壩下游4.76 km 處，由于出山口束窄，造成了水位抬高，出現了約12 m 高水位，給下游二宮村造成了嚴重人員傷亡和經濟損失。

（4）兩次洪峰過程是造成部分人員傷亡的間接原因。整個潰壩過程中水庫下游共出現兩次洪峰，第一次洪峰為溢洪道出流和壩體全斷面過水形成，第二次洪峰為潰壩洪峰。在潰壩洪峰出現前，由于入庫流量減少，在下游二宮村出現第一次洪峰后，水位開始下降，使當地百姓誤認為洪水已退去，已經轉移人員返回村里而造成人員傷亡。

（5）水庫管理站房選址不合理是造成了水庫管理人員傷亡的直接原因。水庫管理站建在壩后下游河灘地，不符合水利樞紐工程布置原則，關鍵時刻自身處于險地，不能正常發揮監測預警等作用。

5 結 論

基于洪災現場實測資料，采用水文、水動力學方法較好地解決了無資料小流域入庫洪水、水庫漲水、溢洪道泄水、漫頂溢水和壩體潰決全過程計算問題，構建二維水動力學模型對潰壩水流的演進進行了數值模擬，并系統分析和總結了射月溝水庫“7·31”暴雨洪水特征及災害形成原因。主要結論如下：

（1）本次射月溝流域各時段降雨量均遠超該地區有記錄以來的歷史極值，短歷時極端強降雨引發遠超水庫防洪能力的特大洪水是造成漫頂潰壩的決定性因素，

（2）采用FFMS 計算得到射月溝水庫入庫洪水過程線，與調查結果較為吻合，其中，模擬洪峰流量為1 915 m3/s，誤差為3.5%。

（3）聯合計算得到水庫出流洪水過程線，其中潰口洪峰流量為7 700 m3/s，此過程線作為二維水動力模型的邊界條件。對模擬結果統計分析，淹沒水深在2～6 m 之間的面積占了總面積的23.02%，受災最重的二宮村二隊、三隊多在這個范圍。

（4）下游河道出現兩次洪峰過程及大壩下游河道束窄是造成下游村莊人員傷亡和經濟財產損失的重要原因。