河堤決口預測及潰堤洪水的礦井風險評價模型*

胡蘇安 何盈利 姜寶元 余子文 高 岳

(中國礦業大學資源與地球科學學院)

我國東部地區許多煤礦毗鄰水庫、河道等地表水系，開采使得諸如地表塌陷、地裂縫等一系列地質災害逐漸產生，并成為地表水下滲的良好通道，尤其是處于河道附近的礦區，在暴雨等極端天氣下，極易發生淹井、涌水、突泥等事故。近年來，針對河道潰堤及潰堤水流的研究主要采用數值模擬方法，劉麗玲等[1]進行了河道與漸潰堤壩耦聯的水力數值模擬；鄭國棟等[2]對Abbott六點中心格式的應用進行了研究，確定了其適用性；陸靈威等[3]通過在大型室內水槽中進行物理模型試驗，模擬分析了潰堤后洪水在潰口外洪泛區內的演進和落水波在河道內的傳播過程；國際上，美國氣象局的潰壩洪水預測模型(DAMBREAK)和簡化潰壩洪水預測模型(SMPDBK)，美國土保局的簡化潰壩演進模型(TR66)，以及荷蘭DELFT大學的洪水系統(DELFTFLS)等相繼問世，并在工程實踐中得到了廣泛應用[4]，為河堤決口預測及潰堤洪水演進分析提供了理論依據。在洪水風險評估方面，陳華麗等[5]綜合考慮了影響洪水災害危險性和社會經濟易損性的因素，對湖北省的洪水災害風險進行了評價；Masatoshi[6]對日本109個河流水系的洪水風險性進行了評價，并進行了不同水系的洪水風險分區研究。本研究在上述成果的基礎上，以河堤決口處的流量、水位等水力參數建立河堤決口預測模型，再引入礦井井口到河堤各區段路線的阻礙因子總值得到潰堤洪水淹沒分區模型，從而對堤壩決口的洪水淹井風險進行安全評價。

1 河堤決口預測模型

1.1 一維單一河道水力模型

1.1.1 一維河道水力模型的基本原理

描述一維水流運動的方程組為Saint-Venant方程組[7]，該方程組建立在質量和動量守恒的基礎上，以水位和流量為研究對象，表達式為

(1)

式中，Q為流量;h為斷面水位;A為過水斷面面積;g為重力加速度;q為河流單位旁側入流量;x為空間坐標;t為時間坐標，C為謝才系數;α為垂向速度分布系數；R為水力半徑，m。

1.1.2 方程組離散

選用Abbott六點中心格式[2]進行方程組求解，河道斷面按照水位(h)—流量(Q)—水位(h)的順序交替布置(其中，Q和h不在同一個斷面上(圖1)，Q總是布置于相鄰的h之間，距離可以不相同)，而后在每個時間段內，利用隱格式的差分法交替計算Q和h[2]。

圖1 河道斷點布置示意

在連續方程中，Q僅對x求偏導，所以方程容易寫成以h為中心的形式(圖2)，動量方程僅以Q為中心(圖3)。

圖2 連續方程Abbott六點中心格式[2]

圖3 動量方程Abbott六點中心格式[2]

在n+1/2時刻，連續方程的離散形式可以表示為

(2)

動量方程的離散形式可以表示為

(3)

式(2)、式(3)中的α、β、γ及δ為方程的計算系數。

1.1.3 離散方程求解

假設1條河流有n個節點，若河道首、末水位點都是h點，n為奇數，于是n個線性方程有n+2個未知數，多出的2個未知數來源于首、末水位點，Zj-1、Zj+1分別變成河流上下游的水位節點，即：h1=Hus，此時系數α1=-1，β1=1，γ1=0，δ1=0[2]。同理，末斷面上hn=Hds，此時系數αn=0，βn=1，γn=-1，δn=0。只要知道上下游水位Hus和Hds，如此n個方程左邊有n個未知數，可以使用標準消元法進行求解。此時，任意節點的變量Z(水位或流量)可以表示為上下游水位節點的函數：

(4)

式中，c、a、b分別為計算系數。

1.2 河堤決口過程模擬

DAMBRK模型[8]將河堤決口概化成梯形，本研究建模時對決口演變過程進行了一定程度簡化處理，不對泥沙輸運進行計算，并且假定決口尺寸按指數形式擴大。決口處的流量采用Freed堰流公式進行計算。對比于其他模型，該模型較簡單，求解所需參數較少，易于進行工程應用。

1.2.1 河堤決口模擬

本研究針對所研究的土堤將決口模擬為梯形(圖4)[9-13]，并按照河堤逐漸潰決進行計算。假設整個潰決過程的總歷時為T，河堤從頂部某一點開始潰決，在潰決過程中，潰口形狀按照線性比率擴大，到達T時刻時，決口寬度達到最大，決口底部高程為最終高程hbm。

圖4 河堤決口示意

潰決模擬過程中，決口的底部高程和底部寬度為時間的函數，隨著時間的變化而變化。決口底部高程計算公式為

(5)

式中，hb為瞬時決口的底部高程,m；hd為堤防總高,m；hbm為決口穩定時的底部高程,m；tb為開始發生潰決至計算時刻經歷的時間；T為潰決到達穩定的總歷時；ρ0為決口處的非線性程度因子，本研究認為河堤決口尺寸按線性變化，故ρ0=1。

決口底部寬度的計算公式為

(6)

式中，bi為瞬時決口的底部寬度，m;b為決口穩定時的底部寬度,m。

1.2.2 河堤決口處流量計算

針對線性梯形決口流量，本研究采用寬頂堰流公式[9-11]進行計算

Qb=cvks[3.1bi(h-hb)1.5]+2.45z(h-hb)2.5，

(7)

式中，cv為流量校正系數；ks為考慮尾水影響出流的淹沒校正系數，由于決口下游為干底，故不考慮決口下游尾水影響，取ks=1.0；h為決口上游水位，m；z為決口邊坡系數。

2 潰堤洪水淹沒分區模型

2.1 層次結構模型構建

首先將整個系統問題進行有條理有層次的分析，從而構建一個具有層次結構的模型，并在層次模型下，將復雜的問題分為系統元素的組成部分。綜合分析該類系統元素，按照其性和關系構成若干個層次，上層元素作為支配下層元素的準則，系統元素可以分為3個層次(最高層、中間層、最底層)。對決口洪水行進過程進行4個方面的影響因素指標綜合分析以及對各種主要評價因素進行劃分，可以構建決口洪水行進過程的影響指標體系，其中最高層為潰堤洪水行進過程影響因素，中間層有地形條件、地表狀況、地表狀況、地質特征和監督管理4個方面的因素，最底層中地形條件包含有坡向、高度和連通性，地表狀況包含有道路覆蓋、房屋覆蓋、植被覆蓋和裸地覆蓋，地質特征包含有砂土、黏土和基巖，監督管理包含有管理制度、日常維護運行、緊急備案維護和其他人為因素。決口洪水行進過程影響指標體系構建完畢后，選取一種科學可靠的分析方法來評價決口洪水行進過程影響指標等級是獲得切合實際的評價結果的關鍵所在。本研究選用層次分析法[14]進行評估。

2.2 建立判斷矩陣

兩兩比較判斷矩陣構建的目的是衡量各個準則指標之間的相對重要性，衡量方法一般采用Satty提出的1～9 值法(表1)，從而得到兩兩比較重要性標度構成的判斷矩陣[15-17]。

表1 判斷矩陣標度定義

2.3 判斷矩陣一致性檢驗

層次分析法中需要對判斷矩陣的一致性進行檢驗，本研究引入了檢查判斷思維的一致性指標CI(Consistency Index)：

CI=(λmax-n)/(n-1) ，

(8)

式中，λmax為判斷矩陣的最大特征值；n為判斷矩陣的階數。

在上述分析的基礎上，查找一次性指標RI(表2)[18]，而后計算一致性比率CR(Consistency Ratio)：

CR=CI/RI .

(9)

通過層次分析法并進行一致性檢驗后，當CR<0.10時，認為判斷矩陣的一致性檢驗是可以接受的，否則，應對判斷矩陣作適當修正。

表2 平均隨機一致性指標[18]

2.4 最終權重計算結果

經過層次分析法計算，得到了如表3所示的各中間層和最底層的權重。

表3 潰堤洪水行進過程影響因素評價體系

3 實例分析

3.1 計算實例

本研究以山東新泰煤礦潰水事故為例進行分析，2007年8月15日夜間開始，山東新汶突降暴雨，山洪暴發，導致柴汶河東都河堤被沖垮，致使河堤發生決口。河堤長度500 m，河堤頂部高程為170 m，計算時間為2007年8月15日0:00:00至2007年8月22日0:00:00，從2007年8月15日20:00:00 開始河堤發生決口。潰決歷時35 h，決口邊坡系數為1.0，決口變化過程呈線性，流量校正系數取1.0。河道上游水位為169 m，河道底部糙率為0.03，底坡坡降為0.000 5，河道總長116 km，發生潰決時河道上游恒定流量為1 800 m3/s，下游恒定流量為900 m3/s，潰決時從側向潰。

3.2 計算結果分析

根據一維河道水力模型和DAMBRK模型進行耦合計算，得到了如圖5～圖7所示的決口流量與決口水位、決口寬度及決口水流流速的過程曲線。

由圖5可知：2007年8月15日夜間開始潰決，隨著時間的增長流量急劇增大，大約35 h后達到峰值流量900 m3/s，隨后流量逐漸減小，大約85 h后流量起伏不大，并與河道中上下游流量相平衡；在此過程中水位變化初期與流量變化趨勢相似，在8月15日夜間潰決開始后流量隨著時間急劇增大達到最大水位169 m，之后水位隨時間降低，且降低速率較快，大約 72 h 后水位達到穩定。

圖6 潰決流量-決口寬度過程曲線

河堤潰決期間，本研究假設決口變化呈梯形線性增大。由圖6可見：在潰決初期，河堤決口隨時間急劇增大，此后決口寬度擴大速度有所減慢，最終達到最大值65 m，之后決口寬度基本固定，變化很小。

圖7 潰決流量-決口流速過程曲線

由圖7可知：潰決流量與流速的變化趨勢一致，在潰決流量增大的同時決口處水流流速也急劇增大。并且，潰決流量達到最大的同時，決口處水流速度也達到最大值，峰值流速之后有所減小并逐漸趨于穩定。

通過潰決流量與決口水位、決口寬度及決口水流流速的過程曲線可以看出，上述模型模擬的河堤潰決過程符合河堤潰決的實際情況。

4 基于ArcGIS的河堤潰決風險分級

4.1 數據來源

為保證最終計算結果的精確性與真實性，必須盡可能采集研究區內的地形、地貌、土地利用情況等資料。對于地形、地貌情況等一般選用DEM數據來表現，土地利用類型等數據可從區域遙感影像上獲取，本研究選用精度為30 m×30 m的DEM數據，DEM數據和遙感影像均來自地理空間數據云。數據在使用前應進行填洼處理等校對處理，剔除“偽地形”等所帶來的影響。

4.2 分析原理

傳統意義上的自然災害風險表達式為[19]

風險=危險性×易損性 .

(10)

一般而言，災害的危險性是研究的前提，承災體的易損性是災害的具體表現，風險則是研究結果。洪水災害所造成的損失和危害在很大程度上取決于承災體的承受能力，即經濟效益的損失程度。受研究區范圍限定，承災體的易損性無法準確估算，且各區段差異基本可忽略，因此該方法的可行性受到了限制。

假設研究區內所有承災體的易損性相同，轉換思維研究洪水行進過程中由于下墊面所產生的水流摩阻損失，將洪水視為“承災體”來求算洪水的“易損性”，用洪水行進路程及河堤距原承災體的距離來估計下墊面對洪水的“危險性”，從而得出洪水的流損情況，進而對堤壩進行風險評估及分級。

4.3 ArcGIS區劃分析

應用ArcGIS軟件強大的空間分析功能，可以對研究區的洪水災害危險性、易損性以及綜合風險進行分析制圖，能有效提高分析結果的精確性，減少不必要的人力、物力和時間上的消耗[20]。首先根據衛星圖與遙感圖進行匹配，選定研究區及附近范圍，進行相應調整，再經ArcGIS軟件操作得出新泰市土地分類圖，如圖8所示。由圖8可知：人類居住區及植被覆蓋區主要分布于河谷地帶，建筑群和植被對洪水行進有較大的阻礙作用；裸巖、裸土區主要分布于近山體側，研究區分布范圍較少，且對洪水行進過程阻礙作用小。

圖8 新泰市土地分類結果

根據校正處理過后的DEM數據，在ArcGIS軟件中構建研究區的基本框架，從而確定區內的基本地形地貌情況，地形越高則表明對洪水行進過程的阻礙作用越大。其中地形標準差主要反映的是地形坡度變化，該值變化越大，對洪水的阻礙能力越弱。

土地分類圖中包含有土地利用類型信息及地質信息，按表3中得到的權重對各影響因素進行賦值；根據研究區DEM分布情況，對各個不同區間的DEM進行賦值；由于監督管理因素在本模型中影響較小，故可忽略該因素的影響。為分析不同潰點對礦井的影響，首先假設洪水行進路線呈直線，計算洪水行進過程阻礙因子，則可得到不同潰點到礦井的阻礙因子總值。若阻礙因子總值越大，說明該潰點礦井危險越小；反之，則說明該潰點礦井危險越大。為計算方便，將研究區的堤壩分為7個區間段，根據得到的各影響因子權重、土地分類和DEM數據，通過ArcGIS軟件進行計算，得到了如表4所示的各區間段路線阻礙因子總值。

表4 各區間段路線阻礙因子總值 ×10-3

根據各區段路線阻礙因子總值，將堤壩等級分為高風險區、較高風險區、中風險區、較低風險區、低風險區，最終得到研究區(圖9)的堤壩風險區劃圖，如圖10所示。

圖9 研究區衛星遙感圖

5 結 語

根據河堤決口處的流量、水位等水力參數建立了河堤決口預測模型，模擬分析了河堤潰決過程。根據ArcGIS軟件建模計算得到了研究區堤壩風險區劃分結果，可知第4區段礦井風險最大，而2007年山東新泰煤礦發生的礦井淹井事故堤壩決口點也在第4區段，說明該模型符合實際，對堤壩潰決預防工作有一定的指導意義。本研究的不足在于：①采用的數據精度不高，在結果處理上可能存在一些系統誤差，與真實情況存在一定的偏差；②洪水災害“危險性”評價指標過于單一，就理論而言強度指標應包括淹沒水深、洪峰流速、洪峰流量、淹沒歷時等，但考慮到研究區范圍受限等各種因素的影響，洪水準確的流速、持續時間等數據難于獲取，本研究僅根據洪水溢流路程進行研究區洪水災害“危險性”評價；③對洪水溢流過程中摩阻損耗的模擬分析采用的是靜態模型，而實際情況可能復雜多樣。

圖10 研究區堤壩風險區劃分結果(單位：m)