基于GIS與水動力模型的陶樂防洪保護(hù)區(qū)漫溢洪水風(fēng)險分析

周躍華，岳志春，，潘汀超，Rizwan Qadir

(1.寧夏回族自治區(qū)水旱災(zāi)害防御中心，寧夏 銀川 750001；2.天津大學(xué)水利工程仿真與安全國家重點實驗室，天津 300072)

0 引 言

黃河寧夏青銅峽到石嘴山河段洪水災(zāi)害易發(fā)。20世紀(jì)青銅峽水文站有記載的洪峰流量大于5 000 m3/s的洪水共5次，洪災(zāi)嚴(yán)重威脅了沿河兩岸廣大人民群眾的生命財產(chǎn)安全。

洪水演進(jìn)數(shù)值模擬是開展防洪保護(hù)區(qū)洪水風(fēng)險分析的主要方法。王揚等[1]基于一維、二維水力學(xué)耦合模型，建立了韓江南北堤防洪保護(hù)區(qū)潰壩洪水演進(jìn)數(shù)學(xué)模型。陳平等[2]構(gòu)建二維非恒定流數(shù)學(xué)模型模擬大陸澤及寧晉泊蓄滯洪區(qū)50年一遇洪水演進(jìn)過程。袁水龍等[3]耦合分布式水文模型MIKE SHE和一維水動力模型MIKE 11研究淤地壩建設(shè)對黃土高原小流域暴雨洪水過程的影響。苑希民等[4]針對山洪溝防洪標(biāo)準(zhǔn)偏低、遭遇暴雨洪水時潰漫堤風(fēng)險較大且較難準(zhǔn)確預(yù)測的問題，系統(tǒng)構(gòu)建潰漫堤洪水多維耦聯(lián)數(shù)值仿真模型。李大鳴等[5]以二維非恒定流控制方程為理論基礎(chǔ)，采用有限體積法，將改進(jìn)的水量平衡模式應(yīng)用于山區(qū)洪水演進(jìn)計算，建立了洪水演進(jìn)數(shù)學(xué)模型。郭立兵等[6]構(gòu)建山洪溝道防沖建筑物壅水風(fēng)險評估模型，能夠較好地模擬寧夏北武當(dāng)溝山洪演進(jìn)過程和防沖建筑物壅水風(fēng)險。田福昌等[7]建立山洪溝道潰堤洪水演進(jìn)一、二維水動力耦合數(shù)值模擬模型，分析評估潰堤山洪淹沒風(fēng)險。針對防洪保護(hù)區(qū)河道洪水的漫灘、潰堤及其與內(nèi)澇疊加耦合預(yù)測的問題，張靖雨等[8]構(gòu)建區(qū)域內(nèi)澇模型疊加聯(lián)用的一、二維水動力耦合模型，模擬中運河南片防洪保護(hù)區(qū)內(nèi)洪水演進(jìn)。苑希民等[9]建立了模擬河道和灌區(qū)洪水演進(jìn)過程的漫潰堤洪水全二維水動力模型。戚藍(lán)等[10]基于高精度DEM數(shù)據(jù)建立三亞市主城區(qū)暴雨、洪水與風(fēng)暴潮多元耦合精細(xì)化洪澇分析模型。

針對狹長河段計算區(qū)域采用二維水動力模型模擬計算效率低的問題，將河道及防洪保護(hù)區(qū)的洪水演進(jìn)概化為一維流動，既可以反映漫溢洪水的特征，又可以減少網(wǎng)格數(shù)量提高模型效率。本文在前人研究成果的基礎(chǔ)上，采用水力學(xué)方法對黃河青石段建立一維非均勻流水動力模型，以陶樂防洪保護(hù)區(qū)為研究對象，對陶樂計算區(qū)段河道進(jìn)行網(wǎng)格加密處理，計算各斷面最高水位，連成水面線；同時結(jié)合保護(hù)區(qū)地形計算保護(hù)區(qū)淹沒水深，為漫溢洪水演進(jìn)仿真模型構(gòu)建、防洪保護(hù)區(qū)洪水風(fēng)險分析提供一定的參考。

1 模型原理

1.1 河道一維水動力學(xué)模型

河道水流流動可采用基于Saint-Venant方程[11-13]構(gòu)建的一維水力學(xué)模型模擬。Saint-Venant方程組為

(1)

(2)

式中，Q為流量；A為斷面面積；q為旁側(cè)入流；t為時間步長；x為空間步長；g為重力加速度；β為流速分布系數(shù)；Z為水面高程；Sf為摩阻損失。

1.2 離散求解

采用 Abbott 六點隱式差分格式來求解Saint-Venant方程組，即

(3)

(4)

式中，h為水位點；Q為流量點；j為計算斷面號；n為時間步。Abbott六點隱式差分格式具體求解過程見參考文獻(xiàn)[14]，該格式具有計算穩(wěn)定、精度高、可靠性強等優(yōu)點。

1.3 漫溢洪水風(fēng)險評估方法

對于灘地漫溢洪水淹沒的淹沒范圍和水深分布，可利用數(shù)值模型計算出的洪水高程作為洪水水位開展淹沒分析[15]。具體來說就是基于河道一維水動力模型計算洪水演進(jìn)過程，然后依據(jù)一維水動力模型模擬獲得的各斷面最高水位開展漫溢洪水影響評估。

首先，利用地理信息技術(shù)，在河道實測斷面的基礎(chǔ)上加密，從DEM數(shù)據(jù)提取河岸高程，將河道斷面線向保護(hù)區(qū)外延，形成加密的保護(hù)區(qū)斷面線，保護(hù)區(qū)斷面線與河道斷面線一一對應(yīng)。繼而從保護(hù)區(qū)DEM數(shù)據(jù)提取河道和灘區(qū)離散點高程，整理為河道一維水動力模型計算需要的斷面數(shù)據(jù)，同一個斷面線上的離散點最高水位通過一維水動力模型求解。

其次，對計算區(qū)域進(jìn)行二維網(wǎng)格剖分，由各離散點的最高水位數(shù)據(jù)內(nèi)插得到各網(wǎng)格水深數(shù)據(jù)。河道一維水動力模型的計算結(jié)果是各個斷面上離散點的最高水位，難以直觀體現(xiàn)洪水的淹沒范圍和水深分布。需要通過空間插值方法，將各個斷面上離散點的最高水位插值生成連續(xù)的洪水水位高程面。自然鄰點插值方法是一種基于空間自相關(guān)性的面積權(quán)重線性內(nèi)插法[16]。其原理是對于一組泰森多邊形，當(dāng)在數(shù)據(jù)集中加入一個新的數(shù)據(jù)點時，就會修改這些泰森多邊形。鄰點的權(quán)重將決定待插點的權(quán)重，待插點的權(quán)重和目標(biāo)泰森多邊形成比例，即

(5)

式中，h(x)為待插點x處的值；αi為自然鄰點xi的權(quán)重；h(xi)為自然鄰點xi處的值。采用自然鄰點插值法對計算網(wǎng)格進(jìn)行插值，具有計算穩(wěn)定、精度高、速度快等優(yōu)點，能夠快速準(zhǔn)確地反映洪水在保護(hù)區(qū)內(nèi)最高水位的分布情況，適用于洪水風(fēng)險應(yīng)急分析。

最后，將插值得到的各網(wǎng)格洪水水面高程與研究區(qū)的地面高程相減得到洪水淹沒水深分布，即

h=hw-hg(hw>hg)

(6)

式中，h為淹沒水深；hw為水面高程插值結(jié)果；hg表示地面高程。利用GIS平臺生成保護(hù)區(qū)淹沒水深圖。

2 應(yīng)用實例

2.1 研究區(qū)域概況

黃河是世界著名的天然多沙河流，全長約5 464 km，流域面積約789.6 km2，黃河寧夏段自中衛(wèi)市南長灘入境至石嘴山市麻黃溝出境，河段長度為397 km，區(qū)間流域面積為520 km2，穿越11個市、縣(區(qū))，地理位置如圖1所示。

圖1 黃河寧夏段地理位置示意

隨著黃河寧夏段上游水庫的建成使用和氣候條件變化，該河段水沙關(guān)系發(fā)生變異并表現(xiàn)為非協(xié)調(diào)性變化。20世紀(jì)80年代以前，黃河寧夏段多年沖淤基本平衡，主河槽能夠保持一定的泄洪輸沙能力，隨后在人類活動與自然氣候變化的雙重影響下，黃河寧夏段上游來水量、汛期輸沙量和造床洪峰流量均大幅度減少，中常洪水過程加長使泥沙主要淤積在主河槽內(nèi)，使得主河槽淤積萎縮日趨嚴(yán)重，平灘流量逐年減少，主河槽輸水輸沙能力嚴(yán)重下降，灘地漫溢淹沒現(xiàn)象頻發(fā)，威脅兩岸人民生命財產(chǎn)安全。陶樂防洪保護(hù)區(qū)位于黃河寧夏段下游，面積約182km2，保護(hù)區(qū)西靠黃河，東臨鄂爾多斯臺地，地形狹長，類似于峽谷地區(qū)河流防洪保護(hù)區(qū)。

本研究需要收集的數(shù)據(jù)主要包括研究區(qū)域內(nèi)的基礎(chǔ)地理資料、歷史洪水災(zāi)害資料、河道斷面資料、1∶10 000地形圖(DLG)、1∶10 000數(shù)字高程模型(DEM)、土地利用圖或近期遙感影像資料、水利普查數(shù)據(jù)、黃河大斷面數(shù)據(jù)、溝渠資料等。對地形圖進(jìn)行配準(zhǔn)后，疊加保護(hù)區(qū)的行政邊界、水系、防洪工程等圖層得到洪水風(fēng)險圖的底圖。

2.2 模型構(gòu)建

采用水力學(xué)方法對陶樂防洪保護(hù)區(qū)進(jìn)行洪水分析計算，建模流程如圖2所示。

圖2 建模流程

2.2.1河道一維水動力學(xué)模型構(gòu)建

為了保證計算精度，一維水動力學(xué)模型計算斷面在2012年實測斷面的基礎(chǔ)上按間距200 m進(jìn)行加密。將河道斷面線向保護(hù)區(qū)外延，形成間距為200 m左右的保護(hù)區(qū)斷面線，保護(hù)區(qū)斷面線與河道斷面線一一對應(yīng)，通過建立一維水動力學(xué)模型，計算河道各斷面最高水位，從而求得保護(hù)區(qū)各斷面最高水位。陶樂保護(hù)區(qū)斷面線的位置如圖3所示。

圖3 陶樂保護(hù)區(qū)斷面位置示意

根據(jù)資料統(tǒng)計，近40年來，寧夏出現(xiàn)具有代表性的洪水年份為1981年和2012年。1981年洪水，下河沿站測得最大洪峰流量為5 880 m3/s，青銅峽站測得最大洪峰流量為6 040 m3/s，石嘴山站測得最大洪峰流量為5 660 m3/s；2012年第3號洪水，下河沿站最大洪峰流量3 520 m3/s，石嘴山站最大洪峰流量3 400 m3/s，青銅峽最大洪峰流量為3 070 m3/s，洪水持續(xù)71天。2012年第3號洪水主要是受黃河上游持續(xù)降雨及水庫下泄等影響，加之近年來黃河寧夏段大洪水極少，河道萎縮淤積嚴(yán)重，導(dǎo)致此次洪水漫灘幾率增加，灘地淤積，主槽沖刷。綜合考慮河道變化及上游水庫調(diào)度的影響，本研究選擇2012年第3號洪水用于率定模型。

模型計算時間步長為5 s，各斷面水深初始值為起算時刻各斷面水深，設(shè)定有利于計算穩(wěn)定和結(jié)果精度。結(jié)果輸出時間步長為1 h，方便結(jié)果展示。

對寧夏黃河2012年第3號洪水過程(2012年7月22日～2012年9月20日)進(jìn)行模擬，提取相關(guān)斷面水位流量數(shù)據(jù)，根據(jù)2012年大斷面實測最高水位進(jìn)行模型驗證。上游入流邊界采用青銅峽水文站2012年第3號洪水過程進(jìn)行模擬水位過程曲線(見圖4)；下游出流控制邊界采用石嘴山水文站水位-流量關(guān)系曲線(見圖5)。利用2012年遙感影像實測最高水位與模擬45個大斷面最高水位對比，模型率定結(jié)果見圖6，其中44個斷面其中誤差在0.2 m以內(nèi)，滿足洪水風(fēng)險分析精度要求，計算精度保證率高達(dá)97.8%。

圖4 青銅峽水文站2012年第3號洪水過程曲線

圖5 石嘴山水文站水位—流量關(guān)系曲線

圖6 斷面率定結(jié)果

用模型對青銅峽水文站100年一遇洪水進(jìn)行分析。上游入流邊界為青銅峽水文站100年一遇流量過程曲線(見圖7)；下游邊界條件為石嘴山水位—流量關(guān)系曲線(見圖4)。陶樂防洪保護(hù)區(qū)灘區(qū)初始視為無積水，按照河道斷面水面線向外延伸確定淹沒范圍和淹沒水深。

圖7 青銅峽水文站100年一遇流量過程曲線

2.2.2防洪保護(hù)區(qū)水深內(nèi)插

陶樂防洪保護(hù)區(qū)洪水風(fēng)險分析依據(jù)青石段一維水動力模型計算獲得的各斷面最高水位開展。計算區(qū)二維淹沒統(tǒng)計時采用的邊界數(shù)據(jù)為一維模型河道水面線計算結(jié)果。根據(jù)水動力模型計算結(jié)果，利用自然鄰點法插值計算得到各網(wǎng)格最高水位，將其與保護(hù)區(qū)地面高程疊加分析，得到洪水的淹沒范圍與淹沒水深分布。

采用離散求解的方法求解陶樂計算區(qū)二維水深。將200 m左右間距的計算區(qū)斷面線離散成50 m間距的離散點，同一個斷面線上的離散點最高水位通過一維水動力模型求解。對陶樂計算區(qū)進(jìn)行二維網(wǎng)格剖分，最大網(wǎng)格面積按1 000 m2進(jìn)行控制，各網(wǎng)格水深數(shù)據(jù)由各離散點的最大水深數(shù)據(jù)內(nèi)插得到。保護(hù)區(qū)網(wǎng)格如圖8所示。各網(wǎng)格地面高程通過保護(hù)區(qū)地形數(shù)據(jù)提取得到，根據(jù)各網(wǎng)格最高水位和地面高程做差求得防洪保護(hù)區(qū)最大淹沒水深分布，結(jié)果大于0的網(wǎng)格即被洪水淹沒的區(qū)域，得到的每個網(wǎng)格的值即為淹沒水深，得到淹沒水深分布圖層。將其與工作底圖疊加和渲染后，繪制得到陶樂防洪保護(hù)區(qū)100年一遇洪水淹沒水深分布圖。

圖8 陶樂計算區(qū)網(wǎng)格劃分示意

2.3 計算結(jié)果分析

陶樂防洪保護(hù)區(qū)靠黃河側(cè)無堤保護(hù)，為狹長的帶狀區(qū)域，地勢較低的區(qū)段主要有紅崖子鄉(xiāng)河段上游部分、五堆子鄉(xiāng)附近河段靠下游部分、陶樂鎮(zhèn)附近河段靠下游部分、高仁鎮(zhèn)河段上游部分，洪水分析計算結(jié)果顯示，地勢較低的部分淹沒范圍和淹沒水深較大，符合實際情況。保護(hù)區(qū)內(nèi)線狀建筑物主要為省道S203，沿程高程均高于河道水位，擋水作用明顯。洪水分析計算結(jié)果顯示道路起到了明顯擋水作用，如圖9所示，計算結(jié)果符合實際情況，滿足合理性要求。當(dāng)陶樂計算區(qū)遭遇100年一遇洪水時，洪水淹沒總?cè)丝跀?shù)量約0.19萬，主要影響月牙湖鄉(xiāng)、高仁鄉(xiāng)、陶樂鎮(zhèn)、紅崖子鄉(xiāng)，總淹沒面積約37.5 km2。

圖9 陶樂防洪保護(hù)區(qū)100年一遇洪水淹沒水深分布

3 結(jié) 論

基于水動力學(xué)理論構(gòu)建了寧夏黃河青銅峽站至石嘴山站河道一維非恒定流水動力模型，采用所建模型計算得到陶樂防洪保護(hù)區(qū)黃河河道最高水面線，利用GIS平臺按洪水淹沒模擬分析結(jié)果確定陶樂計算區(qū)100年一遇灘地漫溢洪水的淹沒水深及淹沒范圍?；诤樗治鲇嬎憬Y(jié)果，按照國家相關(guān)技術(shù)文件規(guī)定，繪制了標(biāo)準(zhǔn)、統(tǒng)一、規(guī)范的洪水淹沒水深圖。

研究成果表明，在寧夏黃河河床淤積抬高的情況下，遭遇洪水水位抬升造成漫灘現(xiàn)象風(fēng)險較大，陶樂防洪保護(hù)區(qū)受災(zāi)較為嚴(yán)重的區(qū)域為蘆葦、王家溝、小紅柳灘、中灘村、下八頃村和黃泥崗等村莊。本文所建模型精度較高，分析結(jié)果可為黃河寧夏段整治工程規(guī)劃設(shè)計及防汛指揮決策提供支持，文中采用的GIS與水動力模型融合方法可為其他區(qū)域淹沒風(fēng)險分析提供技術(shù)參考與應(yīng)用借鑒。