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

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

0 引 言

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

本文對于中小河流的研究主要以飲馬河(石頭口門水庫壩址以下)為研究對象，通過采用MIKE軟件建立了一、二維水動力模型及耦合，模擬多種計算方案下的洪水演進過程，并進行災(zāi)害分析，為日后編制中小河流洪水風(fēng)險圖提供了重要依據(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è)有石頭口門、九臺、德惠水文站，其中石頭口門站多年平均徑流量為7.61 億m3，德惠站多年平均徑流量達9.54 億m3。此外，流域內(nè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 洪水分析方法

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

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

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

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

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

2.3 網(wǎng)格的剖分

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

圖3 飲馬河(石頭口門水庫壩址以下)洪水風(fēng)險圖二維模型計算網(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è)計洪水

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

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

2.5 潰口設(shè)置方案

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

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

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

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

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

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

2.6 計算方案及邊界條件

2.6.1 計算方案

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

2.6.2 邊界條件

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

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

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

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

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

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

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

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

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

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

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

3.2 洪水演進

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

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

3.3 洪水影響分析

針對本次洪水風(fēng)險圖20個計算方案及11個潰口的社會經(jīng)濟情況，對飲馬河(石頭口門水庫壩址以下)防洪保護區(qū)進行災(zāi)情統(tǒng)計分析，主要包括淹沒面積、淹沒房屋面積、受影響公路、鐵路、人口以及GDP等方面。洪水影響圖見圖8，統(tǒng)計成果見表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)險分析能夠很好地基于MIKE模型演示洪水在時空上的變化過程，模擬洪水演進過程。這不僅為日后中小河流洪水預(yù)報奠定了基礎(chǔ)，同時也為洪水風(fēng)險圖的繪制提供了指導(dǎo)作用。就計算成果來看，飲馬河兩岸淹沒面積大，受災(zāi)損失較為嚴重，因此，在避險轉(zhuǎn)移方面，應(yīng)著重做好防范措施，減少人員傷亡和財產(chǎn)損失。

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

□