基于GIS與水動力模型的陶樂防洪保護區漫溢洪水風險分析

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

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

0 引 言

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

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

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

1 模型原理

1.1 河道一維水動力學模型

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

(1)

(2)

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

1.2 離散求解

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

(3)

(4)

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

1.3 漫溢洪水風險評估方法

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

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

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

(5)

式中，h(x)為待插點x處的值；αi為自然鄰點xi的權重；h(xi)為自然鄰點xi處的值。采用自然鄰點插值法對計算網格進行插值，具有計算穩定、精度高、速度快等優點，能夠快速準確地反映洪水在保護區內最高水位的分布情況，適用于洪水風險應急分析。

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

h=hw-hg(hw>hg)

(6)

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

2 應用實例

2.1 研究區域概況

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

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

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

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

2.2 模型構建

采用水力學方法對陶樂防洪保護區進行洪水分析計算，建模流程如圖2所示。

圖2 建模流程

2.2.1河道一維水動力學模型構建

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

圖3 陶樂保護區斷面位置示意

根據資料統計，近40年來，寧夏出現具有代表性的洪水年份為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，洪水持續71天。2012年第3號洪水主要是受黃河上游持續降雨及水庫下泄等影響，加之近年來黃河寧夏段大洪水極少，河道萎縮淤積嚴重，導致此次洪水漫灘幾率增加，灘地淤積，主槽沖刷。綜合考慮河道變化及上游水庫調度的影響，本研究選擇2012年第3號洪水用于率定模型。

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

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

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

圖5 石嘴山水文站水位—流量關系曲線

圖6 斷面率定結果

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

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

2.2.2防洪保護區水深內插

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

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

圖8 陶樂計算區網格劃分示意

2.3 計算結果分析

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

圖9 陶樂防洪保護區100年一遇洪水淹沒水深分布

3 結 論

基于水動力學理論構建了寧夏黃河青銅峽站至石嘴山站河道一維非恒定流水動力模型，采用所建模型計算得到陶樂防洪保護區黃河河道最高水面線，利用GIS平臺按洪水淹沒模擬分析結果確定陶樂計算區100年一遇灘地漫溢洪水的淹沒水深及淹沒范圍。基于洪水分析計算結果，按照國家相關技術文件規定，繪制了標準、統一、規范的洪水淹沒水深圖。

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