基于MIKE模型的中小河流洪水風(fēng)險(xiǎn)分析

孫 依，王 潔，丁 曼，李鴻雁

(1.吉林大學(xué) 新能源與環(huán)境學(xué)院，吉林 長春130021；2.沈陽農(nóng)業(yè)大學(xué)，遼寧 沈陽 110866； 3.吉林省水利水電勘測設(shè)計(jì)研究院，吉林 長春 130021)

0 引 言

近年來，中小河流水害頻繁，且與大江大河相比防洪措施較為薄弱[1,2]。洪澇災(zāi)害作為世界上最為頻發(fā)的自然災(zāi)害之一，也日益成為經(jīng)濟(jì)社會(huì)發(fā)展中面臨的嚴(yán)重問題。2006年末，水利部啟動(dòng)了關(guān)于流域洪水預(yù)警與風(fēng)險(xiǎn)管理技術(shù)研究的項(xiàng)目，以風(fēng)險(xiǎn)管理為基礎(chǔ)，提出了洪水風(fēng)險(xiǎn)圖制作，并將洪水風(fēng)險(xiǎn)圖的編制工作積極推廣至全國[3]。目前，利用MIKE模型模擬洪水演進(jìn)已有大量研究成果。Mark[4]采用隱式差分法求解圣維南方程組，構(gòu)建一維水動(dòng)力學(xué)模型，對(duì)非恒定流過程進(jìn)行了模擬，結(jié)果表明模型適用性良好。田俊玲[5]在一二維模型中引入GIS，建立精細(xì)化仿真模型，結(jié)果可視化更好。Liang D等[6]將動(dòng)態(tài)鏈接庫技術(shù)加入模型當(dāng)中，建立耦合模型模擬潰堤條件下的洪水演進(jìn)情況。此外，由于中小河流日漸成為洪澇災(zāi)害損失的主體[7]，其洪水風(fēng)險(xiǎn)分析也成為眾多學(xué)者研究的重點(diǎn)。劉衛(wèi)林等[8]基于MIKEFLOOD模型對(duì)以羅塘河為代表的中小河流進(jìn)行洪水風(fēng)險(xiǎn)分析，并提出抵御洪水能力的防范措施。李龍[9]采用一二維模型對(duì)牡丹江上游敦化段的洪水特點(diǎn)進(jìn)行風(fēng)險(xiǎn)分析，為敦化市防洪規(guī)劃提供理論支撐。高鵬等[10]采用MIKE模型對(duì)輝發(fā)河干流進(jìn)行淹沒風(fēng)險(xiǎn)分析，并制定合理的避險(xiǎn)轉(zhuǎn)移方案。

本文對(duì)于中小河流的研究主要以飲馬河(石頭口門水庫壩址以下)為研究對(duì)象，通過采用MIKE軟件建立了一、二維水動(dòng)力模型及耦合，模擬多種計(jì)算方案下的洪水演進(jìn)過程，并進(jìn)行災(zāi)害分析，為日后編制中小河流洪水風(fēng)險(xiǎn)圖提供了重要依據(jù)。

1 研究區(qū)概況

飲馬河發(fā)源于伊通縣地局子鄉(xiāng)老爺嶺東南側(cè)[11]，是松花江(三岔河口以上)下游左岸最大的一條支流。流域面積18 247 km2，河道長度387.5 km，河道平均比降約為0.03%。其流域平均高程300 m左右，干流煙筒山水文站以上流域地形均為山區(qū)，生長著茂密的次生林，煙筒山水文站以下流域地形均為低山丘陵區(qū)，且部分土地類型已轉(zhuǎn)化為耕地。石頭口門水庫以下至飲馬河中、下游河段兩側(cè)為平原，地表波狀起伏支流較多[12]。

飲馬河流域水系良好且支流較多，其中伊通河為飲馬河流域最大的支流，集水面積8 440 km2，占飲馬河總面積的46.3%，河道長度343.5 km，河道平均比降0.03%。干流設(shè)有石頭口門、九臺(tái)、德惠水文站，其中石頭口門站多年平均徑流量為7.61 億m3，德惠站多年平均徑流量達(dá)9.54 億m3。此外，流域內(nèi)均勻分布九臺(tái)、雙陽、德惠、農(nóng)安等氣象站，觀測資料精度較高，全年四季分明，降雨多發(fā)生在6-9月，易于7-8月造成流域大暴雨。各河流特征值見表1，水系圖如圖1。

表1 飲馬河主要支流特征值Tab.1 Characteristic value of main tributaries of Yinma River

圖1 飲馬河流域水系示意圖Fig.1 Water system diagram of Yinma River Basin

2 研究方法及模型原理

2.1 洪水分析方法

飲馬河(石頭口門水庫壩址以下)洪水主要來源為石頭口門水庫下泄洪水以及汛期飲馬河干流與支流洪水的匯入。針對(duì)這一特征，本次研究選取一、二維水動(dòng)力學(xué)模型模擬河道和潰堤淹沒區(qū)的洪水演進(jìn)過程，即通過一維水動(dòng)力學(xué)模型建立飲馬河從石頭口門水庫壩址至河口的洪水演算模型，當(dāng)洪水潰堤進(jìn)入淹沒區(qū)后再采用二維水力學(xué)模型計(jì)算。

由于洪水分析計(jì)算分為河道和淹沒區(qū)的洪水演進(jìn)過程的計(jì)算，操作十分復(fù)雜。為了在提高效率的同時(shí)保證計(jì)算結(jié)果的穩(wěn)定可靠，本次洪水分析采用DHI Water and Environment研發(fā)的MIKE 模型系列軟件，包括MIKE 11、MIKE 21和MIKE Flood模塊，可用于一維非恒定流洪水與二維非恒定流洪水的分析模擬。MIKE 11軟件的原理是在圣維南方程和能量方程的基礎(chǔ)上對(duì)河道水動(dòng)力模擬和結(jié)構(gòu)物進(jìn)行計(jì)算[13]，通過隱式有限差分格式對(duì)模型求解。MIKE 21軟件是在二維淺水波方程的基礎(chǔ)上，通過單元中心的顯式有限體積法對(duì)模型進(jìn)行求解[13]。最后，應(yīng)用MIKE Flood模塊將一、二維模型進(jìn)行耦合計(jì)算。

2.2 模型的構(gòu)建范圍

本次洪水風(fēng)險(xiǎn)圖編制范圍為飲馬河石頭口門水庫到飲馬河河口及沿岸的防洪保護(hù)區(qū)。其中，一維河道模型的構(gòu)建范圍為石頭口門水庫到飲馬河河口，范圍內(nèi)的主干河道長度為110 km；二維水流模型的構(gòu)建范圍為飲馬河沿岸的防洪保護(hù)區(qū)，包括A、B、C、D、E五個(gè)分區(qū)，模型構(gòu)建范圍如圖2所示。

圖2 模型構(gòu)建范圍圖Fig.2 Map of model building scope

2.3 網(wǎng)格的剖分

基于已經(jīng)收集到的防洪保護(hù)區(qū)范圍1∶10 000比例尺數(shù)字高程(DEM)、國家2000大地坐標(biāo)系和1985年國家高程基準(zhǔn)，將飲馬河防洪保護(hù)區(qū)按主要水系、地物分割、地形地貌以及堤防情況等劃分成若干區(qū)域。飲馬河左岸霧開河、伊通河及右岸小南河將干流堤防分隔成上、下游相對(duì)獨(dú)立的堤段，初步將左岸分區(qū)為A、B、E三個(gè)計(jì)算分區(qū)，右岸分區(qū)為C、D兩個(gè)分區(qū)。由于無結(jié)構(gòu)網(wǎng)格在復(fù)雜性區(qū)域具有良好的適應(yīng)性，且在局部加密方面具有一定的靈活性，在模擬自然邊界和條件相對(duì)復(fù)雜的水下地形時(shí)，能夠提高邊界精度，使模擬結(jié)果更好。因此本次采用無結(jié)構(gòu)三角形對(duì)防洪保護(hù)區(qū)進(jìn)行剖分，進(jìn)行洪水分析。分析區(qū)域總面積約為1568 km2，控制最大網(wǎng)格面積不大于0.02 km2，網(wǎng)格總數(shù)128 898個(gè)，最大網(wǎng)格為0.019 km2。網(wǎng)格分布如圖3～4所示。

圖3 飲馬河(石頭口門水庫壩址以下)洪水風(fēng)險(xiǎn)圖二維模型計(jì)算網(wǎng)格Fig.3 Two dimensional model computing grid of flood risk map for Yinma River

圖4 飲馬河洪水分析局部網(wǎng)格剖分情況Fig.4 Local grid division for flood analysis of Yinma River

2.4 設(shè)計(jì)洪水

飲馬河(石頭口門水庫壩址以下)洪水風(fēng)險(xiǎn)圖編制區(qū)域的洪水來源為石頭口門水庫下泄洪水以及汛期飲馬河干流與各支流洪水，支流主要包括小南河、霧開河、三道溝與伊通河。洪水主要由暴雨形成，造成暴雨的主要天氣系統(tǒng)是臺(tái)風(fēng)，受暴雨特性的制約，一次暴雨歷時(shí)在1～3 d，但雨量主要集中在24 h以內(nèi)。

飲馬河現(xiàn)有堤防已達(dá)到10～50年一遇標(biāo)準(zhǔn)，設(shè)計(jì)洪水標(biāo)準(zhǔn)為20年、50年和100年一遇。洪水分析方案主要針對(duì)設(shè)計(jì)低標(biāo)準(zhǔn)洪水、標(biāo)準(zhǔn)內(nèi)洪水與超標(biāo)準(zhǔn)洪水3級(jí)[14,15]，本次飲馬河(石頭口門水庫壩址以下)洪水風(fēng)險(xiǎn)圖分析標(biāo)準(zhǔn)內(nèi)洪水采用經(jīng)復(fù)核后的現(xiàn)狀標(biāo)準(zhǔn)。超標(biāo)準(zhǔn)洪水按高一個(gè)等級(jí)設(shè)計(jì)標(biāo)準(zhǔn)的洪水考慮，洪水量級(jí)選取20年一遇洪水，50年一遇洪水和100年一遇洪水。洪水頻率共分以下幾個(gè)方案： ①飲馬河20年一遇洪水與支流小南河、霧開河、伊通河同頻洪水組合；②飲馬河50年一遇洪水與支流小南河、霧開河、伊通河同頻洪水組合；③飲馬河100年一遇洪水與支流小南河、霧開河、伊通河同頻洪水組合。

2.5 潰口設(shè)置方案

潰口的設(shè)置主要根據(jù)“最可能”、“最不利”、“代表性”三個(gè)原則[16,17]，結(jié)合現(xiàn)場調(diào)查歷史災(zāi)害、河道險(xiǎn)工險(xiǎn)段、歷史潰口統(tǒng)計(jì)及洪災(zāi)損失等確定潰口位置，依據(jù)堤防現(xiàn)狀和歷史實(shí)際潰口確定潰口尺寸。最終確定11種方案，其中農(nóng)村堤段7處，城區(qū)段4處。

(1)A區(qū)擬在雙豐村設(shè)置一個(gè)潰口，位于飲馬河左岸九臺(tái)農(nóng)村段堤防，河道凹岸的迎流頂沖位置，且堤線與主流方向接近垂直，歷史上出現(xiàn)潰壩決口。

(2)B區(qū)擬設(shè)置2個(gè)潰口，一個(gè)擬建在望河堡，位于德惠市城區(qū)段堤防，飲馬河干流霧開河匯入口以下，河道凹岸處，主河道臨近堤防，偏于不利考慮，選擇該處潰口。另一個(gè)擬建在陳家崴子，位于德惠市農(nóng)村段堤防，凹岸的迎流頂沖位置，堤防為砂基，堤內(nèi)地勢低洼，積水浸泡堤防，經(jīng)常出現(xiàn)險(xiǎn)情，偏于不利考慮，選擇該處潰口。

(3)C區(qū)擬建2個(gè)潰口，一個(gè)擬建在趙家屯，另一個(gè)擬建在下洼子，二者均位于飲馬河右岸九臺(tái)空港新城區(qū)段河道右側(cè)凹岸，此處有一小溪流匯入，受洪水沖刷易形成塌岸，威脅堤防安全。

(4)D區(qū)擬在飲馬河右岸設(shè)置4個(gè)潰口，分別為蓮花村潰口、明月樓潰口、五家子潰口、合發(fā)屯潰口。其中蓮花村潰口位于老城區(qū)堤防，小南河匯入凹岸的迎流頂沖位置，且堤線與主流方向接近垂直，堤防為砂基，歷史上出現(xiàn)過險(xiǎn)工險(xiǎn)情。明月樓潰口位于九臺(tái)農(nóng)村段堤防，河道凹岸處，經(jīng)歷洪水對(duì)堤岸的沖刷，水流對(duì)堤岸易產(chǎn)生淘刷，有潰堤危險(xiǎn)。五家子潰口位于德惠市農(nóng)村段堤防，凹岸的迎流頂沖位置，由于其附近有一拆除的五家子涵洞，且堤防為砂基，在堤內(nèi)產(chǎn)生沖坑及管涌，易產(chǎn)生潰堤危險(xiǎn)故選擇該處潰口。合發(fā)屯潰口位于德惠市農(nóng)村段堤防，飲馬河干流伊通河匯入口以下河道凹岸處，考慮支流洪水匯入后易產(chǎn)生沖刷潰堤危險(xiǎn)，故選擇此處潰口。

(5)E區(qū)擬在農(nóng)村段堤防設(shè)置2個(gè)潰口，分別為臥牛石丁字堤潰口和東排木涵洞潰口，位于飲馬河干流伊通河匯入口以下河道凹岸處，由于其歷史上出現(xiàn)潰壩決口，易發(fā)生潰堤危險(xiǎn)，因此選擇此處潰口。

2.6 計(jì)算方案及邊界條件

2.6.1 計(jì)算方案

本次潰口寬度參考1985年堤防潰口的情況確定。根據(jù)調(diào)查了解，飲馬河1985年洪峰流量為921 m3/s，相當(dāng)于5年一遇洪水標(biāo)準(zhǔn)。干流堤防潰口寬度從40～200 m不等，一般潰口寬度在100 m以下，最大的潰口寬度達(dá)到200 m左右。本次潰口寬度的設(shè)置根據(jù)洪水量級(jí)情況，分別設(shè)置3個(gè)等級(jí)：①當(dāng)發(fā)生20年一遇洪水時(shí)，各潰口寬度取50 m；②當(dāng)發(fā)生50年一遇洪水時(shí)，各潰口寬度取100 m；③當(dāng)發(fā)生100年一遇洪水時(shí)，各潰口寬度取200 m。根據(jù)計(jì)算分區(qū)、洪水量級(jí)和潰口分布，共設(shè)置計(jì)算方案 20 個(gè)。

2.6.2 邊界條件

(1)一維水動(dòng)力模型邊界條件。由河網(wǎng)分布情況可知，以石頭口門水庫的設(shè)計(jì)放流過程作為模型的上邊界條件，以飲馬河河口水位-流量關(guān)系作為模型下邊界條件，各支流匯入點(diǎn)的設(shè)計(jì)流量過程作為各區(qū)間來水的邊界條件。邊界入流及各區(qū)間匯流均按1953年典型洪水過程推算而得。上、下邊界及區(qū)間來水邊界情況見圖5。

圖5 邊界條件示意圖Fig.5 Schematic diagram of boundary conditions

(2)潰堤洪水模型邊界條件。潰堤模型需要處理的邊界條件為潰口形態(tài)的變化過程，采用預(yù)先制定的潰口形態(tài)變化過程結(jié)果。由潰堤洪水模型與一、二維水流模型進(jìn)行聯(lián)立求解，三者連接處的邊界條件不需設(shè)定。

(3)二維水動(dòng)力模型邊界條件。堤防潰口處與潰堤模型連接，不需人工設(shè)定，其他堤防及各交通干線都只考慮漫溢，不考慮潰決。

3 洪水演進(jìn)及結(jié)果分析

3.1 模型參數(shù)的選取及率定

一維河道模型參數(shù)的選取主要分為主槽和灘地，主槽糙率主要由河床組成、水流形態(tài)和岸壁特性三方面因素確定[18]。由于飲馬河干流主槽河底底坡較均勻，石頭口門水庫上游段河床多為礫質(zhì)粗砂、中砂，下游多為細(xì)砂、灰壤土，水流通暢且岸壁整齊，滿足天然河道糙率表中單式斷面第Ⅳ類，故主槽糙率值選用0.03。邊灘糙率值主要由斷面形態(tài)、床質(zhì)和植被因素確定[18]，由于研究區(qū)兩側(cè)灘面寬闊平坦、土壤肥沃，多數(shù)被墾為良田，局部筑有民堤，局部河段形成阻水，滿足灘地糙率表中第Ⅲ類特征，因此，邊灘糙率值選用0.05。二維水動(dòng)力模型參數(shù)主要為糙率、初始水深、干濕水深等[19]，淹沒區(qū)下墊面類型及糙率的選取直接影響著洪水分析的結(jié)果。故本次依據(jù)研究區(qū)內(nèi)土地類型對(duì)區(qū)域地物進(jìn)行概化，將淹沒區(qū)分為6種地表類型[10,20]，各類糙率值見表2。

表2 下墊面糙率表Tab.2 Roughness of underlying surface

在對(duì)模型進(jìn)行率定的過程中，主要選取2010年汛期7月21日-9月7日實(shí)測洪水過程，來率定一維河道模型參數(shù)。為了保證模型的穩(wěn)定性和在沿程水位精度計(jì)算上的準(zhǔn)確性，將模型中的計(jì)算時(shí)間步長設(shè)置為30 s，河道斷面最小間距設(shè)置為50 m，水位計(jì)算點(diǎn)的最小間距也設(shè)置為50 m。以石頭口門水庫實(shí)測放流過程和各支流實(shí)測流量過程作為飲馬河的入流過程，通過調(diào)整飲馬河的河床主槽及灘區(qū)的糙率參數(shù)，使德惠水文站的水位、流量與實(shí)測水位、實(shí)測流量在最大水平上達(dá)到吻合。模擬值與實(shí)測值對(duì)比成果見圖6。當(dāng)洪峰流量誤差小于10%，水位誤差小于20 cm時(shí)[21]，河道參數(shù)取值比較合理。將率定后的河道糙率分布對(duì)1973年洪水進(jìn)行驗(yàn)證，德惠站的計(jì)算與實(shí)測結(jié)果如圖7所示，由圖7可見，洪水的漲落趨勢與實(shí)測數(shù)據(jù)擬合較好，流量過程誤差小于10%。驗(yàn)證結(jié)果符合要求。

圖6 德惠站模擬流量、水位與實(shí)測流量、水位過程對(duì)比圖(2010年)Fig.6 Comparison between simulated flow (water level) and measured flow (water level) of Dehui station (2010)

圖7 德惠站模擬流量與實(shí)測流量過程對(duì)比圖(1973年)Fig.7 Comparison between simulated flow and measured flow of Dehui station (1973)

3.2 洪水演進(jìn)

由于飲馬河設(shè)置的潰口較多，且不同潰口在潰決后的淹沒過程具有較大差異，本次僅以飲馬河干流發(fā)生50年一遇洪水，各支流按相應(yīng)洪水組合為例，進(jìn)行洪水演進(jìn)，并分析其淹沒情況。

飲馬河農(nóng)村段堤防五家子潰口在發(fā)生潰決1 h后，保護(hù)區(qū)淹沒面積為0.63 km2；潰決發(fā)生3 h后，淹沒面積為3.24 km2；潰決發(fā)生8 h后，淹沒面積為8.97 km2；潰決發(fā)生1 d后，洪水受到高速公路G1和長圖線的限制，在高速公路和鐵路中間的地帶向下游演進(jìn)，淹沒面積為19.17 km2；潰決發(fā)生3 d后，洪水穿過高速公路G1向下游演進(jìn)，淹沒面積為47.45 km2；潰決發(fā)生8 d后，洪水穿過高速鐵路向下游演進(jìn)，淹沒面積為108.48 km2。整個(gè)潰決分洪過程中，受到鐵路、道路以及高地勢的阻水作用，洪水始終沿著鐵路和堤防中間的地帶向下游演進(jìn)，隨著淹沒水深的增大，直到淹沒鐵路及道路后向下游演進(jìn)。從不同時(shí)刻的淹沒水深及分布情況，可以看出模型模擬的結(jié)果符合由高向低的水流走勢，同時(shí)也能看出阻水建筑物對(duì)洪水演進(jìn)具有一定影響。

3.3 洪水影響分析

針對(duì)本次洪水風(fēng)險(xiǎn)圖20個(gè)計(jì)算方案及11個(gè)潰口的社會(huì)經(jīng)濟(jì)情況，對(duì)飲馬河(石頭口門水庫壩址以下)防洪保護(hù)區(qū)進(jìn)行災(zāi)情統(tǒng)計(jì)分析，主要包括淹沒面積、淹沒房屋面積、受影響公路、鐵路、人口以及GDP等方面。洪水影響圖見圖8，統(tǒng)計(jì)成果見表3。

圖8 飲馬河干流洪水影響圖Fig.8 Flood impact map of main stream of Yinma River

表3 飲馬河洪水影響分析成果表Tab.3 Analysis of flood impact of Yinma River

4 結(jié)論與建議

本次中小河流洪水風(fēng)險(xiǎn)分析能夠很好地基于MIKE模型演示洪水在時(shí)空上的變化過程，模擬洪水演進(jìn)過程。這不僅為日后中小河流洪水預(yù)報(bào)奠定了基礎(chǔ)，同時(shí)也為洪水風(fēng)險(xiǎn)圖的繪制提供了指導(dǎo)作用。就計(jì)算成果來看，飲馬河兩岸淹沒面積大，受災(zāi)損失較為嚴(yán)重，因此，在避險(xiǎn)轉(zhuǎn)移方面，應(yīng)著重做好防范措施，減少人員傷亡和財(cái)產(chǎn)損失。

此外。在研究過程中仍存在不足，飲馬河石頭口門水庫壩址以下到河口，有若干條支流匯入，本次計(jì)算只考慮相對(duì)較大支流匯入，實(shí)際水系水流情況復(fù)雜。后續(xù)研究應(yīng)增加實(shí)地調(diào)研，定期收集水文數(shù)據(jù)、水力特性數(shù)據(jù)和洪災(zāi)損失調(diào)查，積累大量基礎(chǔ)資料。在進(jìn)行二維水流計(jì)算時(shí)，網(wǎng)格大小的設(shè)定對(duì)計(jì)算結(jié)果及精度具有顯著影響，網(wǎng)格剖分越小，計(jì)算精度越高，但同時(shí)也會(huì)降低運(yùn)算軟件的運(yùn)行速度，所以在后續(xù)的研究中對(duì)網(wǎng)格的剖分應(yīng)做到更加細(xì)致，使結(jié)果更精確。

□