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

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

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

1 河堤決口預(yù)測(cè)模型

1.1 一維單一河道水力模型

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

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

(1)

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

1.1.2 方程組離散

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

圖1 河道斷點(diǎn)布置示意

在連續(xù)方程中，Q僅對(duì)x求偏導(dǎo)，所以方程容易寫(xiě)成以h為中心的形式(圖2)，動(dòng)量方程僅以Q為中心(圖3)。

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

圖3 動(dòng)量方程Abbott六點(diǎn)中心格式[2]

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

(2)

動(dòng)量方程的離散形式可以表示為

(3)

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

1.1.3 離散方程求解

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

(4)

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

1.2 河堤決口過(guò)程模擬

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

1.2.1 河堤決口模擬

本研究針對(duì)所研究的土堤將決口模擬為梯形(圖4)[9-13]，并按照河堤逐漸潰決進(jìn)行計(jì)算。假設(shè)整個(gè)潰決過(guò)程的總歷時(shí)為T(mén)，河堤從頂部某一點(diǎn)開(kāi)始潰決，在潰決過(guò)程中，潰口形狀按照線(xiàn)性比率擴(kuò)大，到達(dá)T時(shí)刻時(shí)，決口寬度達(dá)到最大，決口底部高程為最終高程hbm。

圖4 河堤決口示意

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

(5)

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

決口底部寬度的計(jì)算公式為

(6)

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

1.2.2 河堤決口處流量計(jì)算

針對(duì)線(xiàn)性梯形決口流量，本研究采用寬頂堰流公式[9-11]進(jìn)行計(jì)算

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

(7)

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

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

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

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

2.2 建立判斷矩陣

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

表1 判斷矩陣標(biāo)度定義

2.3 判斷矩陣一致性檢驗(yàn)

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

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

(8)

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

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

CR=CI/RI .

(9)

通過(guò)層次分析法并進(jìn)行一致性檢驗(yàn)后，當(dāng)CR<0.10時(shí)，認(rèn)為判斷矩陣的一致性檢驗(yàn)是可以接受的，否則，應(yīng)對(duì)判斷矩陣作適當(dāng)修正。

表2 平均隨機(jī)一致性指標(biāo)[18]

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

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

表3 潰堤洪水行進(jìn)過(guò)程影響因素評(píng)價(jià)體系

3 實(shí)例分析

3.1 計(jì)算實(shí)例

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

3.2 計(jì)算結(jié)果分析

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

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

圖6 潰決流量-決口寬度過(guò)程曲線(xiàn)

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

圖7 潰決流量-決口流速過(guò)程曲線(xiàn)

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

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

4 基于ArcGIS的河堤潰決風(fēng)險(xiǎn)分級(jí)

4.1 數(shù)據(jù)來(lái)源

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

4.2 分析原理

傳統(tǒng)意義上的自然災(zāi)害風(fēng)險(xiǎn)表達(dá)式為[19]

風(fēng)險(xiǎn)=危險(xiǎn)性×易損性 .

(10)

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

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

4.3 ArcGIS區(qū)劃分析

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

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

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

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

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

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

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

5 結(jié) 語(yǔ)

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

圖10 研究區(qū)堤壩風(fēng)險(xiǎn)區(qū)劃分結(jié)果(單位：m)