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

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

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

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

1 河堤決口預測模型

1.1 一維單一河道水力模型

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

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

(1)

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

1.1.2 方程組離散

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

圖1 河道斷點布置示意

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

圖2 連續(xù)方程Abbott六點中心格式[2]

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

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

(2)

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

(3)

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

1.1.3 離散方程求解

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

(4)

式中，c、a、b分別為計算系數(shù)。

1.2 河堤決口過程模擬

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

1.2.1 河堤決口模擬

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

圖4 河堤決口示意

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

(5)

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

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

(6)

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

1.2.2 河堤決口處流量計算

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

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

(7)

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

2 潰堤洪水淹沒分區(qū)模型

2.1 層次結(jié)構(gòu)模型構(gòu)建

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

2.2 建立判斷矩陣

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

表1 判斷矩陣標度定義

2.3 判斷矩陣一致性檢驗

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

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

(8)

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

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

CR=CI/RI .

(9)

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

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

2.4 最終權(quán)重計算結(jié)果

經(jīng)過層次分析法計算，得到了如表3所示的各中間層和最底層的權(quán)重。

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

3 實例分析

3.1 計算實例

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

3.2 計算結(jié)果分析

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

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

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

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

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

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

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

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

4.1 數(shù)據(jù)來源

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

4.2 分析原理

傳統(tǒng)意義上的自然災害風險表達式為[19]

風險=危險性×易損性 .

(10)

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

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

4.3 ArcGIS區(qū)劃分析

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

圖8 新泰市土地分類結(jié)果

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

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

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

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

圖9 研究區(qū)衛(wèi)星遙感圖

5 結(jié) 語

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

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