蓄滯洪區洪水演進模擬及堤防潰決損失評估方法

蔣水華， 黃中發， 黃勁松， 張 穎， 江先河

(1.南昌大學 建筑工程學院， 江西 南昌 330031； 2.江西水利職業學院， 江西 南昌 330013)

1 研究背景

蓄滯洪區災害模擬、風險分析、預警，以及洪水演進分析、評價、決策一直都是水文水資源領域和防災減災部門關注的重大課題[1]。蓄滯洪區在洪水來臨期間降低下游洪災損失的同時，也要考慮降低蓄滯洪區內進洪過程中的損失，而且隨著時間的推移，蓄滯洪區內產業經濟已得到發展，人口也有所增長，分洪時將造成較大損失，故有必要模擬蓄滯洪區洪水演進過程，估計對應的生命、經濟和生態環境損失，從而指導防洪搶險。

洪水演進數值模擬是洪水風險分析及其損失評估的關鍵步驟，是進行洪水風險管理的中心環節[2]。童漢毅等[3]、Tucciarelli等[4]、李大鳴等[5]分別采用有限差分法、有限元法、有限體積法模擬蓄滯洪區洪水演進過程，取得了較好的數值結果，但是采用這些方法建立的模型推導和計算過程復雜，計算量較大。隨著計算機和地理信息系統技術的發展，采用一維或二維水力學方法模擬洪水演進過程發展迅速，水文學方法得到較為廣泛的應用[6]。魏凱等[6]、穆聰[7]等、王崇浩等[8]、郭鳳清等[9]、朱世云等[10]、袁雄燕[11]等、Morales-Hernández等[12]和Bladé等[13]利用一維或二維水力學方法建立MIKE21 FM模型，分別模擬了蓄滯洪區、大壩潰壩、河湖和堤防潰決的洪水演進過程，結果表明二維水力學模擬方法可在短時間內直觀地顯示洪水淹沒范圍的變化，提供較詳細的淹沒區信息，可較好地根據數值模擬結果制定相關決策指導防洪搶險工作。黃琳煜等[14]基于MIKE FLOOD建立MIKE URBAN與MIAKE 21耦合模型，并構建了浦東川沙地區耦合模型來分析出現澇點和積水的原因，取得了較好的研究成果。王揚等[15]歸納總結了MIKE內部各個模型的基本原理、特點及其適用范圍，并且指出了MIKE模型存在的局限性。

針對蓄滯洪區堤防的規模較大，而現有的防洪標準較低、堤身堤型不達標、區域內人口密度大、水保工程措施不到位，且區內產業嚴重失衡、種植業較多以及蓄洪排澇成本較大等因素引起的洪水演進過程中損失估算的難題。本文發展了堤防潰決洪水演進數值模擬方法，建立了基于蓄滯洪區洪水演進淹沒數據的損失評估方法。依托鄱陽湖康山蓄滯洪區，采用MIKE21建立蓄滯洪區洪水演進數值模型，進而根據洪水演進過程模擬結果(淹沒范圍、水深和流速)估計堤防潰決造成的生命、經濟和生態環境損失。

2 洪水演進數值模擬方法

丹麥MIKE21模型[8-10,16]是常用的洪水演進數值模擬軟件，分為海岸水文學和海洋學、環境水文學、泥沙傳輸過程和波浪4個模塊，其原理是采用二維水動力學數學模型模擬水流運動，直觀得出水位、流量和流速等重要的水力指標隨時間的變化關系。相比于其他模型，MIKE21模型的功能較強大，不僅可以進行軟啟動，還可以設置多種控制性水工結構，進行干、濕節點及單元設置等[11]，在我國乃至世界范圍內得到了廣泛使用，比如在我國南水北調工程和長江綜合治理等工程中均得到了成功應用。

2.1 控制方程

MIKE21模型將河道水流看成不可壓縮的牛頓液體，采用Navier-Stokes方程組來描述河道水流運動規律[17]。為了考慮河床底部摩阻作用和紊動影響，引進渦黏系數，建立平面二維水動力數值模型。模型控制方程包括平面二維流連續方程，計算表達式為:

(1)

式中：t為計算時間，s；x、y、z為右手Catesian坐標系；h為靜止水深,m；u、v為流速在x、y方向上的分量，m/s；S為源匯項。x方向上的二維水流動量方程為:

(2)

2.2 邊界條件

邊界條件主要包含上、下邊界條件、特殊邊界條件以及動邊界條件，具體如下：

(1)上邊界條件。蓄滯洪區上邊界條件一般選取邊界處的流量或水位過程，根據研究目的不同，可選取歷史最高水位或流量過程、設計水位或流量過程、實測時段水位或流量過程等作為對應的上邊界條件。

(2)下邊界條件。通常滯洪區在運行時為封閉區域，所以模型不設置下邊界條件；若不是封閉區域，下邊界條件一般為邊界處的水位或流量過程。

(3)特殊邊界條件。對起到排水作用的構筑物(涵管、涵洞等)和阻水作用的建筑物(道路、堤防等)進行概化，使得模擬更加符合客觀實際。

(4)動邊界條件。動邊界是指計算區域中隨水位變化的邊界，為了保證模型計算的連續性，采用MIKE21模型中的干、濕邊界來設置動邊界條件[9]。一般來說，干濕水深分別采用系統默認值0.005 m與0.1 m可以滿足計算精度，具體含義為當計算區域的水深小于0.005 m時，該邊界條件則為干邊界，不參與計算；當水深大于0.1 m時，該邊界條件則為濕邊界，重新參與計算。

2.3 參數選取

MIKE21模型參數包括數值參數和物理參數。數值參數一般取默認值，物理參數主要包括河床糙率系數、動邊界計算參數和渦黏系數等。MIKE21水動力數學模型的主要參數為糙率系數，該參數是一個反映水流阻力的綜合系數。研究區下墊面糙率系數可依據《洪水風險圖編制導則》[18]給出的糙率變化范圍取值或考慮研究區域河床現狀進行率定取值。

2.4 潰口設置

蓄滯洪區分洪口門包括自然分洪、閘門分洪和爆破分洪這3種分洪方式，都須預留分洪口，根據國家防汛抗旱總指揮部《關于長江洪水調度方案的批復》(國訊[2011]22號)[19]，當湖口水位達到分洪水位時，采取相關工程措施將堤防預留的分洪口扒開，蓄滯洪區開始進洪。參照分洪口門的設置參數，進口形態接近流線型，故可采用流線型寬頂堰公式計算潰口流量[20]：

(3)

式中：Qb為潰口處出流流量，m3/s；Z為潰口處河道水位,m；m和σ分別為流量系數和淹沒系數，根據文獻[20-21]計算確定，沒有考慮淹沒系數與流量系數存在的不確定性；Zb為潰口頂部高程，m；B為潰口寬度，m。

3 堤防潰決損失評估方法

當發生超標準洪水時，為防洪搶險，將預留的堤防分洪口扒開，進而通過堤防分洪口的洪水會對下游蓄滯洪區造成一定的經濟、生命、環境損失。生命損失是指堤防下游地區受潰堤洪水影響造成的死亡或受災人口數量。經濟損失是指潰堤洪水帶來的直接經濟損失和間接經濟損失，可用貨幣量表示。生態環境損失指的是堤防潰決引起的資源破壞和環境污染對生態環境造成的影響。

關于生命損失，結合宋敬衖等[22]的研究成果，根據人口聚居特點結合行政區劃將淹沒區劃分為若干子區域，使用風險人群積分算法得到實用的生命損失LOL計算公式如下：

(4)

式中：PARij為第i個區域中第j組風險人口數量；a為風險區域總數；c為風險人口分組數；IRij為第i個區域中第j組風險人群個體生命損失率，計算公式為:

(5)

式中：Pi為潰堤概率；fi為發生第i等級洪水個體死亡率，可根據洪水演進模擬的淹沒范圍、水深以及預警時間等因素查表得到；l為洪水等級數目[22]。此處fi為傳遞參量，需要根據淹沒區受災程度以及預警時間來確定。

經濟損失主要包括提防的破壞、下游淹沒區城鎮的損害、工程收益的損失等。根據損失特征采用損失率方法估算直接經濟損失，該方法適用于估算各類流動資產與固定資產的損失[22]，間接經濟損失一般取直接經濟損失的0.63倍[23-24]。

(6)

式中：S為按損失率計算得直接經濟損失；Si為第i類財產損失；βijk、Wijk分別為第i類第j種財產在第k種淹沒程度下的損失率和財產價值；e為財產分類總數；h為第i類財產分類總數；l為淹沒水深等級。此處βijk、Wijk為傳遞參量，Wijk需要根據受災程度選擇相應研究范圍，再到經濟庫行政面積中去查找當地經濟數據確定；βijk需要根據淹沒程度確定。

國際上對堤防潰決造成的環境損失研究較少，并且由于考慮的因素多且復雜，沒有統一的風險標準，一般較難直接統計計算。為此，本文采取條件價值評估方法評估堤防潰決造成的環境損失。條件價值評估方法是評估環境非使用價值的主要方法[25]，利用該方法評價環境損失的本質就是通過問卷調查方式來調查風險區內人口愿意為改善該環境功能的支付意愿，或者是放棄該環境功能的補償意愿，以此揭示環境對于被調查者的價值，即該區域環境資源的經濟價值[25]，其中人均支付意愿E(元/a)和環境損失量F的計算公式為:

(7)

F=E·N·f

(8)

式中：Ai為每人愿意支付的金額；Di為個體選擇該金額的概率；u為金額的具體數量；N為年限；f為人口數。此處E和F為傳遞參量，需要根據淹沒歷時來計算環境損失隨時間的變化關系。

4 工程應用

4.1 工程概況

康山大堤興建于1966年，屬于四級堤防，全長36.25 km。設計洪水標準：保證水位22.50 m(吳淞高程，湖口站水位)，堤頂高24.55 m，內外邊坡1∶3，防汛警界水位19.50 m。穿堤建筑物有梅溪咀閘(中型水閘，設計排澇流量為103 m3/s)、鑼鼓山電排站(大(2)型泵站，設計排澇流量73.8 m3/s)、大湖口閘(中型水閘，設計排澇流量為160.5m3/s)、落腳湖電排站(中型泵站，設計排澇流量15.5 m3/s)、里溪閘(小型水閘，設計排澇流量為7.1 m3/s)[26]。

康山大堤內垸1985年被國務院批準列為國家蓄滯洪區，有效蓄洪量16.58×108m3，是江西省最大的蓄滯洪區，也是國家重點蓄滯洪區，擔負著15×108m3分蓄洪任務。總集雨面積450.31 km2，蓄洪面積312.37 km2。區內有石口鎮、三塘鄉、大塘鄉、瑞洪鎮、康墾總場、康山鄉等6個鄉鎮場和1個現代農業示范區，常住人口2萬余戶、10余萬人，耕地面積117.21 km2，養殖水面面積139.01 km2。經統計，康山大堤保護面積343.4 km2，保護耕地面積94.2 km2，保護人口103 437人、22 560戶，康山蓄滯洪區居民土地財產登記情況有：承包土地262.49 km2，包括農作物117.21 km2，水產類139.01 km2，經濟林6.27km2；專業養殖家畜類11 919頭，家禽類198 384只；住房19 432間，合計843 353 m2；農業機械7 244臺，役畜類1 395頭，居民主要耐用消費品17 250臺[26]。

4.2 研究區地形概化及網格剖分

康山大堤潰堤洪水演進過程模擬的計算區域基礎地理信息數據比例尺為1∶10000，行政區劃數據比例尺為1∶50000。所收集的電子地圖需按要求用國家大地2000坐標系(CGCS2000)，高程基準統一采用1985年國家高程基準。其中1∶10000的數字高程模型如圖1所示。為了分析潰堤洪水演進過程，獲得潰口流量變化過程和堤防保護區淹沒信息，采用公式(1)和(2)建立康山蓄滯洪區水流運動平面二維非恒定數學模型。

模型計算范圍為整個康山蓄滯洪區，以北面康山大堤的分洪口為進洪邊界與退洪邊界，以北面康山大堤及南面陸地為物理邊界。由于計算區域地形復雜，為了合理布置網格，采用MIKE21模型中的FM非結構化網格(不規則的三角形網格)類型，利用有限體積法進行數值模擬。最大網格面積不超過0.01 km2，對居民區、建筑物密集區和種植區的計算網格適當進行加密處理。研究區域共剖分了58 095個網格節點和63 880個計算單元，網格中心點之間的最大間距為300 m，最小間距為26 m，蓄滯洪區的網格剖分如圖2所示。所涉及的物理參數取值如表1和2所示，具體包括邊界條件、初始條件、渦黏系數和經率定后的糙率等。需要說明的是，表1中采用1954年4月15日0時至1954年4月18日0時的洪水數據進行洪水過程模擬，這是因為1954年長江中下游發生了近100年最大洪水，洪水數據非常具有代表性。

表1 模型計算物理參數取值

表2 研究區域內不同下墊面糙率取值

圖1 康山蓄滯洪區數字高程模型圖(1∶10000)

圖2 蓄滯洪區網格剖分圖

4.3 洪水潰決過程及流量分析

康山蓄滯洪區在康山大堤預留有分洪口門。分洪口位于大堤樁號20+070～20+450 m之間，分洪口門寬380 m，在分洪口門兩端布置深層水泥攪拌樁以防止口門無限制擴大。

根據國家防汛抗旱總指揮部《關于長江洪水調度方案的批復》(國訊[2011]22號)[19]，當鄱陽湖湖口水位達到20.59 m(吳淞高程22.50 m)時，對應的康山站水位20.68 m(吳淞高程22.55 m)，采用人工爆破的方式將預留潰口扒開，口門底寬為12.5 m，堤防潰決，蓄滯洪區開始進洪。口門爆破并沖刷擴大的過程按2 h計，最終口門底寬取300 m。參考《江西省鄱陽湖蓄滯洪區安全建設工程可行性研究報告》[26]關于規劃口門底高程的選取，取口門底高程15.93 m(黃海高程)。潰決口門發展所需時間為口門初始形態發展到終止形態的時間，其中口門結構參數和口門發展時間需要結合當地或類似地區已發生過的實際潰決調查資料合理確定。不考慮區間降雨計算口門的發展過程，得到潰口流量變化過程曲線見圖3。由圖3可見，當防洪堤潰決后，潰口流量在2 h內很快達到最大值14 000 m3/s。24 h以后潰口流量開始下降，至48 h完成整個進洪過程，蓄滯洪區內水位達到20.68 m(黃海高程)，與外湖水位一致，潰口不再進流，流量為0。

圖3 康山大堤潰口流量過程線

4.4 洪水演進過程模擬

為了便于分析，在蓄滯洪區的康山鄉、瑞洪鎮、大塘鄉、三塘鄉、石口鎮古竹片區、康山墾殖場6個行政區的23個行政村(自然村)內布置了23個特征點，各特征點布置及潰口位置見圖4。堤防潰決后研究區各典型時段的淹沒深度見圖5。

圖4 康山蓄滯洪區特征點布置及潰口位置

通過模型計算得到堤防潰決后研究水域各典型時段的淹沒水深，見圖5。由圖5可知，堤防潰決1h后，水流到達康山墾殖場和大湖管理局，12 h后，洪水可以推進到瑞洪鎮片區的東源村，瑞洪鎮片區的東三村由于地勢較高，洪水需要在24 h左右才能達到淹沒深度。相比之下，地勢較高、遠離龍口的大塘安全區在36 h后洪水才到達。由于康山墾殖場落腳湖分場和康山鄉的臨湖區域在口門附近，這兩處的最大流速大于1.0 m/s，受洪水沖刷影響較大。瑞洪鎮江家村附近部分地形高程小于20 m，當蓄滯洪區分洪36 h后，西面安全區開始受影響。由圖5(c)可知，分洪結束時蓄滯洪區平均淹沒水深達6.95 m，康山墾殖場及石口鎮部分區域淹沒水深達到8 m。

圖5 堤防潰決后研究水域各典型時段的淹沒深度

4.5 潰堤洪水淹沒面積

康山大堤分洪口門潰決時，潰口寬380 m，口門底高程15.93 m。潰決前，康山蓄滯洪區內蓄洪底部水位為15.10 m，約69 km2(去除水域面積)被漬水淹沒。整個潰堤洪水演進過程歷時48 h，進洪總量為16.58×108m3，總的淹沒面積為288 km2。

按照淹沒嚴重程度對淹沒區域進行分區，并統計各分區的淹沒面積如表3所示。定義危險區淹沒深度大于2.0 m，重災區淹沒深度在1.0～2.0 m之間，中災區淹沒深度在0.5～1.0 m之間，輕災區淹沒深度小于0.5 m，安全區洪水未到達區域，該區域是人員轉移安置區。

由表3可知，輕災區、中災區、重災區所占淹沒面積的百分比較小，分別為0.66%、0.82%、1.77%。相比之下，危險區淹沒面積為278.63 km2，占總淹沒面積的96.75%，說明在這種工況下，康山蓄滯洪區受災程度大且涉及范圍廣。將MIKE計算結果文件導入ArcGIS，計算得到洪水未到達區域的面積為31.3 km2，該區域亦為安全區，可作為人員轉移安置區。

通過數值模擬可以得出不同淹沒程度的面積分布情況，從而可為康山蓄滯洪區洪水風險分析、抗洪搶險以及人員轉移安置提供較為準確的參考依據。

表3 潰堤洪水淹沒面積統計表

4.6 特征點淹沒水深、流速

以康山鄉、瑞洪鎮、大塘鄉、三塘鄉、石口鎮古竹、康山墾殖場片區內23個特征點為代表，分析蓄滯洪區淹沒水深及流速變化情況，如圖6所示，從圖6中可獲得洪水到達特征點需要的時間、特征點的最大淹沒水深、流速峰值和峰現時間等。

由圖5和6可知，康山鄉雖離潰口較近，但地勢較高，直至堤防潰決后9 h 40 min洪水才達到府前村，最大淹沒水深和流速峰值均出現在府前村，分別為2.65 m和0.0133 m/s，峰現時間為潰堤后30 h 20 min，大山、金山村未被洪水淹及。

在潰堤9 h 10 min后洪水達到瑞洪鎮西崗村，瑞洪鎮最大淹沒水深出現在西崗村為4.55 m，最大流速出現在羅家村為0.0536 m/s，峰現時間為潰堤后24 h 30 min。其中該片區的把下村因地勢較高，未被洪水淹及。

圖6 不同片區特征點淹沒水深及流速過程線

大塘鄉在潰堤50 min后開始淹沒，最大淹沒水深和最大流速均出現在大湖村，分別為7.65 m和0.0731 m/s，峰現時間為2 h 40 min。其中該片區陳家塘村、南垅村因地勢較高，未被洪水淹及。

石口鎮古竹片區在潰堤50 min后開始淹沒，最大淹沒水深出現在古竹村為6.24 m；最大流速出現在西村為0.027 m/s，峰現時間為潰堤后28 h 50min。劉埠村和湖濱村由于地勢較高，直至潰堤后48 h進洪完成，未被洪水淹及。

三塘鄉因地處康山蓄滯洪區南部山區，離潰口較遠，且地勢較高，潰堤后23 h 40 min，洪水才達到三塘鄉馬山村，最大淹沒水深3.67 m，最大流速為0.2635 m/s，峰現時間為潰堤后29 h。魏家、余家橋村未被洪水淹及，該未被洪水淹及的區域可作為避險轉移的重點區域。

以上6個區域中康山墾殖場片區最先被淹沒，其中插旗分場在潰堤后30 min開始淹沒，里溪分場、甘泉洲分場在潰堤后40 min開始淹沒。康山墾殖場片區最大淹沒水深和最大流速均出現在里溪分場，最大淹沒水深為8.32 m，最大流速為0.074 m/s，峰現時間為潰堤后2 h 30 min。

此外由圖6可知，康山墾殖場片區最先被淹沒，幾乎都被淹沒，且淹沒水深平均超過2 m，故該片區人口物資需要在預警時間內迅速撤離至離該片區較近的石口鎮古竹片區安全樓或大塘鄉片區安全樓；古竹片區的安全樓可建立在劉埠村和湖濱村，這兩個村未被淹及；大塘鄉片區的安全樓可建立在陳家塘村和南垅村；瑞洪片區的把下村也未被淹沒，可將該區域的西崗村和羅家村的人口和物資遷至把下村避險；康山鄉片區的大山村和金山村從蓄洪開始到結束均未被淹沒，可在此處建立安全樓以便康山鄉片區人口物資避險；三塘片區最后被淹及，且僅有馬山村被淹沒，其余地區未遭到淹沒，故該片區人口物資相對安全。

4.7 淹沒區損失統計

基于上述獲得的淹沒范圍、深度、流速、預警時間等信息采用公式(4)和(5)估計淹沒區生命損失，步驟如下：(1)通過淹沒信息劃分災區，統計計算各災區的風險人口。(2)通過洪水歷時、水深、淹沒范圍來查表確定死亡率，將死亡率與潰堤概率相乘得到損失率。(3)將步驟(1)和(2)計算結果相乘得到生命損失。獲得的淹沒區風險人口累積曲線和不同預警時間內的人口損失如圖7(a)～7(c)所示。

圖7 潰堤造成的研究區內生命、經濟和環境損失圖

經濟損失估算步驟如下：(1)在經濟數據庫中查詢康山蓄滯洪區的行政面積，求出模擬單元格中的受災單元格相對于行政面積的比例；各類財產貨幣化后的價值與比例的乘積即為該網格單元內各類型的財產價值。(2)由洪水演進模擬得到的淹沒信息(水深、流速和淹沒歷時)確定單元內各類財產的損失率。(3)將步驟(1)得到的財產價值與步驟(2)得到的損失率相乘，得到該單元格的財產損失，累積計算淹沒區域財產損失即為該蓄滯洪區的經濟損失。根據統計蓄滯洪區淹沒范圍內的直接經濟組成成分，采用公式(6)計算得到區內直接經濟損失和農林牧及家庭財產損失累積曲線如圖7(d)～7(e)所示。淹沒范圍內所有單元的損失累加，得到直接經濟損失為9.53×108元，間接損失取直接損失的63%，為6.00×108元，則康山大堤潰決總經濟損失為15.53×108元。

環境損失通過問卷調查的方式評估，采用公式(7)和(8)計算得到湖區重點堤防整治工程居民的平均支付意愿為每人776.8～780元/a。目前康山蓄滯洪區共有居民29 895人，按人均776.8元/a，按20 a計算，堤防環境整治工程總價值評估為4.64×108元，即堤防潰決導致的社會環境損失為4.64×108元。進而按照洪水淹沒進程得到環境損失累積曲線，如圖7(f)所示。

基于蓄滯洪區洪水演進過程數值模擬得到的洪水淹沒信息，采用第3節生命、經濟和環境損失評估方法可得到風險人口、淹沒-損失人口、經濟損失和環境損失，據此可根據不同預警時間下的淹沒損失制定相應的應急預案和具體的抗洪搶險措施。比如由圖7(a)可知，該蓄滯洪區內的風險人口隨著洪水歷時的一個變化過程，可結合淹沒范圍信息制定準確詳細的人口遷移方案，以減少淹沒損失和降低抗洪成本。圖7(b)～7(c)除了評估人口損失之外，還可根據不同預警時間編制搶險方案，制定預警級別和相應的響應措施。根據圖7(e)可判斷直接損失中損失比值最大的成分為農業經濟損失，從而為提前采取相關止損措施提供參考依據。

5 結 論

本文發展了堤防潰決洪水演進數值模擬方法，依托鄱陽湖康山蓄滯洪區，采用MIKE21建立洪水演算數值模型模擬蓄滯洪區洪水演進過程，獲得了最高水位條件下的淹沒范圍、洪水淹沒深度及流速過程線，進而估算了堤防潰決造成的蓄滯洪區內的生命、經濟和生態環境損失。主要結論如下：

(1)獲得了康山蓄滯洪區歷史最高水位條件下的潰堤淹沒信息(包括淹沒范圍、水深和流速等)，并根據洪水淹沒范圍、水深劃分了區內災害等級，最后根據潰決洪水到達特征點所需的時間、特征點最大淹沒水深、流速峰值和峰現時間確定了各片區安全樓的位置，以便制訂人口和物資遷移方案，這些研究結果可以有效指導蓄滯洪區的防洪搶險和人口轉移工作。

(2)建立了基于蓄滯洪區洪水演進淹沒數據的損失評估方法，基于潰堤淹沒數據得到了蓄滯洪區內各片區的三大損失累積曲線。并根據獲得的風險人口、淹沒區人口損失、經濟損失和環境損失判斷了風險嚴重程度，從而為堤防潰決風險評估、康山及類似蓄滯洪區防汛搶險決策方案的制定提供了重要的參考依據。