向枝葉， 叢沛桐

(華南農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院， 廣州 510642)

受氣候變化和人類活動(dòng)影響，極端降雨和大洪水發(fā)生的概率顯著增加[1-2]，給大中型水庫帶來了較大的超蓄淹沒風(fēng)險(xiǎn)，開展庫區(qū)超蓄淹沒風(fēng)險(xiǎn)動(dòng)態(tài)模擬和三維展示研究，對(duì)于防洪減災(zāi)工作的快速落實(shí)具有重要意義。

近年來，二維水動(dòng)力模型已成為研究洪水演進(jìn)規(guī)律、評(píng)估洪水風(fēng)險(xiǎn)的重要方法和手段。李大鳴等[3-4]提出了單元水量出流修正法，建立了基于改進(jìn)水量平衡模式的二維洪水演進(jìn)模型，并用于蓄泄洪區(qū)和水庫下游山區(qū)的洪水模擬。劉鵬飛等[5]采用基于二維淺水方程的CJK3D-WEM軟件模擬了長(zhǎng)江下游河段的潮流影響。張大偉等[6]為模擬存在挖沙坑、橋墩、丁壩等的復(fù)雜河段水流，建立了基于Godunov格式的精細(xì)化二維水動(dòng)力模型，可用于復(fù)雜研究區(qū)域的防洪評(píng)價(jià)工作。在洪水風(fēng)險(xiǎn)可視化研究方面，許小華等[7]基于洪水分析模型實(shí)現(xiàn)了洪水風(fēng)險(xiǎn)實(shí)時(shí)分析系統(tǒng)，結(jié)合地圖服務(wù)可進(jìn)行洪水風(fēng)險(xiǎn)二維平面展示。牛帥等[8]考慮了降雨、潮位等動(dòng)態(tài)因子輸入，建立了基于一、二維水動(dòng)力模型的易澇區(qū)動(dòng)態(tài)洪水風(fēng)險(xiǎn)分析模型。王峰等[9]利用傾斜攝影技術(shù)，整合了三維場(chǎng)景構(gòu)建、洪澇數(shù)據(jù)及數(shù)字高程模型(digital elevation model，DEM)數(shù)據(jù)，開發(fā)了Flood-3DGIS系統(tǒng)，可三維展示高精度的城市洪澇淹沒過程。

目前，在洪水風(fēng)險(xiǎn)分析及可視化展示方面的研究多集中在水庫下游、蓄泄洪區(qū)等區(qū)域，而關(guān)于水庫不同來水、不同調(diào)度方案下的庫區(qū)動(dòng)態(tài)超蓄淹沒實(shí)景三維可視化展示的研究較少。現(xiàn)以韓江高陂水利樞紐庫區(qū)為研究對(duì)象，利用激光雷達(dá)(light detection and ranging，LIDAR)數(shù)據(jù)獲取高分辨率地形，為更加準(zhǔn)確地概化水庫，反映真實(shí)庫容情況和淹沒情形，模擬堤防對(duì)洪水的阻擋作用，建立基于無結(jié)構(gòu)網(wǎng)格有限體積法和Villemonte公式的動(dòng)態(tài)超蓄淹沒風(fēng)險(xiǎn)分析模型，模型可輸入任意來水、調(diào)度方案進(jìn)行動(dòng)態(tài)洪水風(fēng)險(xiǎn)模擬，并通過系統(tǒng)實(shí)現(xiàn)實(shí)景地物三維淹沒情景展示，以期對(duì)風(fēng)險(xiǎn)預(yù)警預(yù)報(bào)、水庫調(diào)度運(yùn)行、避洪轉(zhuǎn)移提供重要的指導(dǎo)作用。

1 二維超蓄淹沒模型

1.1 基本控制方程

均勻流體流動(dòng)滿足Boussinesq理論且服從靜水壓力假設(shè)即為淺水流動(dòng)。模型采用淺水方程的守恒形式描述洪水在庫區(qū)的演進(jìn)過程，忽略水面風(fēng)應(yīng)力和柯氏力的二維淺水方程組的矢量形式[10-11]為

(1)

(2)

(3)

(4)

(5)

(6)

式中：U為守恒向量；F、G分別為x、y方向的通量向量；S為源項(xiàng)向量；h為水深，m；u、v分別為x、y方向的流速分量，m/s；zb為河床底高程，m；qL為單位面積旁側(cè)入流源項(xiàng)，m/s；Sfx、Sfy為河床底摩阻項(xiàng)；n為糙率。

淺水方程的空間離散方法為基于無結(jié)構(gòu)三角形網(wǎng)格單元中心格式的有限體積法，界面法向數(shù)值通量的計(jì)算采用Roe格式的近似Riemann解，時(shí)間積分采用Euler向前格式。

1.2 堰流公式

一般在二維水動(dòng)力模型中，線狀阻水構(gòu)筑物可通過加密網(wǎng)格抬高地形或使用堰流公式兩種方法添加至模型中[12-14]。采用Villemonte堰流公式精細(xì)化模擬線性構(gòu)筑物[15]，可避免因局部加密網(wǎng)格，網(wǎng)格過渡不均勻造成模型失穩(wěn)問題的發(fā)生，表達(dá)式為

(7)

式(7)中：Qs為過堤流量，m3/s；W為堤頂寬度，m；Hw為堤頂高程，m；C為流量系數(shù)，取值為1.838；Hus為上游水位，m；Hds為下游水位，m。

2 模型構(gòu)建

2.1 項(xiàng)目區(qū)概況

韓江高陂水利樞紐(以下簡(jiǎn)稱樞紐)位于廣東省梅州市大埔縣高陂鎮(zhèn)上游約5 km處的韓江干流上，如圖1所示，是韓江流域防洪供水體系中的控制性水利工程，壩址以上集雨面積26 590 km2，正常蓄水位38 m，設(shè)計(jì)洪水位、校核洪水位47.44 m，總庫容3.656×108m3，防洪庫容達(dá)2.673×108m3。項(xiàng)目主體工程已于2021年1月下閘蓄水，與上游永定棉花灘水庫及下游堤防共同構(gòu)成“堤庫結(jié)合”的防洪體系，保護(hù)人口624.29萬人，為粵東地區(qū)經(jīng)濟(jì)社會(huì)的可持續(xù)發(fā)展提供安全保障。

圖1 韓江高陂水利樞紐位置示意圖Fig.1 Schematic diagram of the location of Gaobei hydro-junction

工程移民征地范圍由停發(fā)流量(6 700 m3/s)回水線確定，樞紐防洪庫容全部置于正常蓄水位以上，下閘蓄洪庫區(qū)將產(chǎn)生超蓄淹沒，居民遷移線至設(shè)計(jì)洪水位回水線之間仍存有大量的人口、房屋及文物單位等，一旦庫區(qū)超蓄將造成難以估量的損失，甚至?xí)?duì)社會(huì)穩(wěn)定形成沖擊[16]。對(duì)庫區(qū)開展超蓄淹沒風(fēng)險(xiǎn)研究是十分必要的。

2.2 地形數(shù)據(jù)處理及網(wǎng)格剖分

模型地形數(shù)據(jù)采用近期的LIDAR實(shí)測(cè)數(shù)據(jù)，LIDAR數(shù)據(jù)經(jīng)解算、分類處理成點(diǎn)數(shù)據(jù)，高程點(diǎn)經(jīng)檢查后生成Delaunay不規(guī)則三角網(wǎng)并輸出為匹配模型網(wǎng)格大小的地形，得到分辨率為2 m的數(shù)字高程模型(digital elevation model，DEM)，詳盡的地形數(shù)據(jù)為精細(xì)化建模提供了基礎(chǔ)。

二維建模范圍由試算的100年一遇洪水淹沒范圍外擴(kuò)確定，總面積約為188 km2，覆蓋庫區(qū)正常運(yùn)用淹沒面積27.78 km2，網(wǎng)格剖分邊長(zhǎng)沿河道邊界向河道兩岸外側(cè)由最小10 m到最大60 m漸變，網(wǎng)格角度控制在30°～120°，網(wǎng)格質(zhì)量滿足精度和計(jì)算速度的要求。地形插值方法為自然鄰點(diǎn)插值法，建模區(qū)地形如圖2所示。

圖2 建模區(qū)地形Fig.2 Topography of modeling area

2.3 定解條件

邊界條件和初始條件合稱定解條件，是求解非恒定流問題的必要條件。為減少模擬累積誤差、提高模型穩(wěn)定性，上游邊界給定上游來水過程，下游控制邊界給定調(diào)度模型輸出的庫水位過程；閉合邊界選擇無滑移條件。

以壩址遭遇30、50、100年設(shè)計(jì)頻率洪水的情況的計(jì)算工況為例，進(jìn)行模型分析，各重現(xiàn)期入庫洪峰流量、庫水位如表1所示。

干濕邊界采用單元分區(qū)法處理，設(shè)置干水深hdry=0.005 m，濕水深hwet=0.05 m，由此為依據(jù)將研究區(qū)網(wǎng)格單元分為干區(qū)、半干半濕區(qū)及濕區(qū)，干區(qū)不參與計(jì)算，半干半濕區(qū)僅進(jìn)行水量通量計(jì)算。

初始條件為匹配下游邊界條件設(shè)為正常蓄水位38 m。

表1 各重現(xiàn)期洪峰流量-水位成果表Table 1 Flood peak discharge and water level results

2.4 線性構(gòu)筑物概化

庫區(qū)內(nèi)需考慮堤防、公路等工程對(duì)洪水的阻擋作用，將堤防和公路概化為線性構(gòu)筑物添加至二維模型中。建模區(qū)內(nèi)堤防工程包括三河鎮(zhèn)匯城堤、大麻鎮(zhèn)的附麻堤、黨溪村的塘溪堤，其中附麻堤和塘溪堤現(xiàn)狀已完成堤防達(dá)標(biāo)建設(shè)，建設(shè)標(biāo)準(zhǔn)為20 a一遇洪水，堤頂高程如表2所示。區(qū)域內(nèi)的主要交通通道為省道(S222、S333)縣道、鄉(xiāng)道以及村道，主要研究省道和縣道的阻水影響，鄉(xiāng)道和村道因阻水作用不明顯而不作概化。

表2 堤頂高程表Table 2 Elevation of embankment top

2.5 參數(shù)設(shè)置

二維模型糙率取值采取分區(qū)設(shè)置，將下墊面分類為河道、村鎮(zhèn)、道路、水田、旱地、林地、空地共七類，其中韓江主槽糙率依據(jù)廣東省水利電力勘測(cè)設(shè)計(jì)研究院實(shí)測(cè)資料推求的斷面高水位糙率分段設(shè)置，在0.030～0.045；除河道以外的區(qū)域根據(jù)下墊面的土地利用類型，參照糙率經(jīng)驗(yàn)取值表賦初值，不同下墊面糙率取值如表3所示。

表3 不同下墊面糙率取值Table 3 Roughness values of different underlying surfaces

2.6 結(jié)果與分析

2.6.1 合理性分析

(1)水量平衡分析。由入流邊界、模型模擬結(jié)果分別統(tǒng)計(jì)出上游來水量、下游出流量以及樞紐時(shí)段庫容量，水面蒸發(fā)損失折算取值為0.39 m3/s，計(jì)算的相對(duì)誤差為1.3×10-7，小于1×10-6，滿足《洪水風(fēng)險(xiǎn)圖編制導(dǎo)則》(SL 483—2017)技術(shù)要求。

(2)如圖3所示，30、50、100年一遇洪水工況下，超蓄淹沒模型模擬的水面線結(jié)果與已批復(fù)的初步設(shè)計(jì)報(bào)告成果差值絕對(duì)值在0.00～0.13 m，說明模擬結(jié)果與初設(shè)計(jì)算成果基本相符。

(3)淹沒范圍合理性分析。結(jié)合DEM數(shù)據(jù)進(jìn)行對(duì)比分析，30、50、100年一遇方案淹沒形態(tài)及范圍大致相同，洪水演進(jìn)范圍主要集中在原河道，兩岸平原低洼地區(qū)被淹沒，結(jié)果與研究區(qū)地形走向相符，地形高區(qū)未發(fā)生淹沒，且本次模擬的淹沒范圍與已批復(fù)的高陂水利樞紐超蓄淹沒影響專題報(bào)告基本一致，說明淹沒范圍分布合理。

綜上，水量相對(duì)誤差、水面線、淹沒范圍等成果均與實(shí)際吻合，說明模型合理可靠。

圖3 韓江干流水面線成果Fig.3 Water surface profile results of Hanjiang main stream

2.6.2 超蓄淹沒分析

當(dāng)壩址遭遇30年一遇的洪水時(shí)，天然情況下河道兩岸已產(chǎn)生淹沒，樞紐下閘蓄洪將加大庫區(qū)上游淹沒風(fēng)險(xiǎn)。如表4所示，遭遇30、50、100年一遇洪水時(shí)，臨時(shí)淹沒面積(居民遷移線對(duì)應(yīng)范圍不納入統(tǒng)計(jì))分別為9.09、11.43、13.78 km2，其中水深大于4.0 m的淹沒面積分別占總面積的56.6%、62.9%、70.6%。庫區(qū)兩岸大多為陡峭的山體，洼地平原分布其中，隨著洪水量級(jí)的增大，兩岸山體的阻擋使洪水演進(jìn)范圍集中在主河道，河道水位迅速上漲，因此0.0～4.0 m水深范圍的淹沒面積增幅不大；同時(shí)洪水涌入洼地平原區(qū)，積水加深，水深大于4 m的淹沒面積增幅明顯。

表4 淹沒水深分布統(tǒng)計(jì)Table 4 Statistics of submerged depth distribution

如圖4所示，遭遇30～100年一遇洪水時(shí)，韓江沿岸地區(qū)的淹沒水深基本在0～14 m，淹沒深度由上游向下游呈遞增趨勢(shì)。高陂鎮(zhèn)、銀江鎮(zhèn)的銀灘村、大麻鎮(zhèn)、三河鎮(zhèn)以及茶陽鎮(zhèn)均為易受洪災(zāi)區(qū)域。

2.6.3 淹沒損失評(píng)估

根據(jù)本次淹沒模擬結(jié)果和區(qū)域社會(huì)經(jīng)濟(jì)情況，對(duì)不同頻率的洪水淹沒損失進(jìn)行統(tǒng)計(jì)評(píng)估，包括受災(zāi)人口、淹沒房屋、淹沒耕地面積、主要專業(yè)項(xiàng)目等指標(biāo)，統(tǒng)計(jì)結(jié)果見表5。

3 動(dòng)態(tài)淹沒模擬三維可視化系統(tǒng)

因洪水現(xiàn)象的隨機(jī)性和調(diào)度方案的多樣性，為滿足實(shí)際調(diào)度會(huì)商需求，設(shè)計(jì)了動(dòng)態(tài)淹沒模擬三維可視化系統(tǒng)。系統(tǒng)采用B/S(瀏覽器/服務(wù)器)的模式，數(shù)據(jù)交換接口接入任意來水、調(diào)度方案等動(dòng)態(tài)因子，通過模型管控平臺(tái)調(diào)用水文模型、水庫調(diào)度模型，計(jì)算輸出來水過程和庫水位過程進(jìn)入超蓄淹沒模型進(jìn)行風(fēng)險(xiǎn)模擬；傾斜攝影原始數(shù)據(jù)經(jīng)后處理軟件處理，完成圖片和地形的融合，構(gòu)建三維實(shí)景模型，最后將超蓄淹沒結(jié)果與三維實(shí)景模型相融合，進(jìn)行淹沒風(fēng)險(xiǎn)(淹沒范圍、水深等)的三維可視化展示，系統(tǒng)技術(shù)流程如圖5所示，三維實(shí)景效果如圖6所示，100年一遇樞紐超蓄淹沒風(fēng)險(xiǎn)三維可視化展示如圖7所示。

表5 高陂水利樞紐庫區(qū)超蓄淹沒損失統(tǒng)計(jì)表Table 5 Statistics of flooding loss of over-impoundment of Gaobei hydro-junction

圖4 不同頻率淹沒水深圖Fig.4 Different frequency submerged depth map

圖6 三維實(shí)景效果圖Fig.6 3D terrain rendering

4 結(jié)論

為研究高陂水利樞紐庫區(qū)的超蓄淹沒風(fēng)險(xiǎn)，建立了動(dòng)態(tài)超蓄淹沒風(fēng)險(xiǎn)分析模型，模擬了30、50、100年一遇設(shè)計(jì)洪水淹沒情景，得出以下結(jié)論。

(1)從水量平衡、水面線誤差、淹沒范圍與地形走向相關(guān)性進(jìn)行了合理性分析，結(jié)果表明，超蓄淹沒模型結(jié)果合理可靠，符合實(shí)際，可用于高陂水利樞紐超蓄淹沒風(fēng)險(xiǎn)分析。

(2)當(dāng)遭遇30年一遇以上洪水時(shí)，高陂水利樞紐庫區(qū)超蓄淹沒影響的主要范圍為高陂鎮(zhèn)、銀江鎮(zhèn)的銀灘村、大麻鎮(zhèn)、三河鎮(zhèn)以及茶陽鎮(zhèn)等區(qū)域。

(3)動(dòng)態(tài)超蓄淹沒三維可視化系統(tǒng)以模型管控平臺(tái)為智能引擎，驅(qū)動(dòng)水文預(yù)報(bào)模型、水庫調(diào)度模型，實(shí)現(xiàn)來水和調(diào)度方案的動(dòng)態(tài)輸入，模擬庫區(qū)超蓄淹沒風(fēng)險(xiǎn)，并利用三維可視化技術(shù)，直觀展示超蓄淹沒范圍、水深、轉(zhuǎn)移路線等風(fēng)險(xiǎn)信息。這對(duì)于災(zāi)前預(yù)警預(yù)報(bào)、調(diào)度決策和快速撤離受災(zāi)群眾具有重要的指導(dǎo)意義，能為其他地區(qū)洪水風(fēng)險(xiǎn)模擬及三維可視化展示系統(tǒng)設(shè)計(jì)提供思路和參考。

圖7 100年一遇洪水淹沒風(fēng)險(xiǎn)三維展示Fig.7 3D display of 100 a flood inundation risk