基于洪水演進模擬的潰堤洪災(zāi)損失評估

李文歡,支歡樂,蔣水華*,雷 聲,黃 河

(1.南昌大學(xué)工程建設(shè)學(xué)院,江西 南昌 330031;2.中鐵水利水電規(guī)劃設(shè)計集團有限公司,江西 南昌 330029;3.江西省水利科學(xué)院,江西 南昌 330029)

洪水作為世界上最常見且最具威脅的自然災(zāi)害之一,嚴(yán)重影響了人們生產(chǎn)生活,制約了中國經(jīng)濟社會的發(fā)展。中國水旱災(zāi)害公報資料表明[1],2019年,中國因洪澇災(zāi)害共造成4 766.6萬人受災(zāi),658人死亡失蹤,直接經(jīng)濟損失1 922.7億元。因此,開展洪水災(zāi)害風(fēng)險評估,進而建立洪水災(zāi)害預(yù)警系統(tǒng),制定有效的防災(zāi)減災(zāi)措施具有重要的工程應(yīng)用價值。

目前,眾多學(xué)者在洪水演進分析方面開展了大量有益的研究[2-4],并發(fā)現(xiàn)相較于Delft3D[2]、HEC-RAS[3]和Delft3D-Flow[4]等軟件,MIKE 21軟件不僅可以考慮風(fēng)速和入滲等多種因素的影響,而且模擬結(jié)果更準(zhǔn)確,計算效率更高。例如,劉衛(wèi)林等[5]針對水文資料缺乏的中小河流通過小流域水文計算和MIKE 11水動力模型進行研究,提出了臨界水位預(yù)警指標(biāo)確定方法。在水動力學(xué)模型方面,詹明強等[6]研究了潰壩對下游淹沒范圍和程度的影響,建立了一維非恒定流水動力學(xué)模型。

洪澇災(zāi)害損失評估主要分為生命損失、經(jīng)濟損失和生態(tài)環(huán)境損失這3類[7]。對于生命損失評估,目前主要參考了國外經(jīng)驗并結(jié)合國內(nèi)情況估算洪災(zāi)人口死亡率和生命損失[8]。郭磊等[9]通過夜間燈光的遙感數(shù)據(jù)和土地利用類型推求人口空間展布數(shù)據(jù),再疊加歷史洪水模擬得到受災(zāi)人口數(shù)量。經(jīng)濟損失評估常分為直接和間接損失評估,如Pariartha等[10]采用MIKE FLOOD計算洪水淹沒水深和淹沒面積,疊加分析住宅地理信息,再根據(jù)住房損失曲線和平均房價來評估洪水造成的經(jīng)濟損失。目前關(guān)于洪水造成的生態(tài)環(huán)境損失研究較為欠缺,這方面的成果不多,國外則主要以定性方法評估洪水對于生態(tài)環(huán)境的影響。近年來發(fā)展迅猛的地理信息技術(shù)受到了學(xué)者們的青睞,與洪水災(zāi)害損失評估結(jié)合進行洪水淹沒空間分布數(shù)據(jù)分析和承災(zāi)體空間信息災(zāi)害評估[8-10]。

采用MIKE 21軟件模擬潰堤洪水演進過程,根據(jù)土地利用類型不同合理確定下墊面糙率值,并基于Python對ArcGIS進行二次開發(fā),在此基礎(chǔ)上建立具有明確物理意義的損失定量評估模型,進而量化經(jīng)濟財產(chǎn)價值和生態(tài)系統(tǒng)服務(wù)價值,結(jié)合潰堤洪水時空分布數(shù)據(jù),建立洪水災(zāi)害損失評估模型。最后,以鄱陽湖區(qū)某重點實際工程為例驗證方法的有效性,并評估潰堤洪水造成的生命、經(jīng)濟和生態(tài)環(huán)境損失。

1 潰堤洪水演進模擬方法

水動力學(xué)模型常被用于模擬洪水過程演進過程。MIKE 21軟件是目前國際上流行的水動力學(xué)模擬軟件[11],常用于模擬潰堤/潰壩洪水演進過程。Symonds等[2]分析對比Delft3D和MIKE 21的模擬效果和計算效率,表明兩者的模擬結(jié)果都與實測數(shù)據(jù)吻合良好,但是MIKE 21在多核運行時使用非結(jié)構(gòu)化模型,計算效率更高。Shrestha等[3]分別采用HEC-RAS和MIKE 21對比分析了布倫特伍德地區(qū)Deer Creek河流的洪水演進過程,表明兩者都可以精確模擬洪水演進過程,但是相對于HEC-RAS,MIKE 21還可以考慮風(fēng)速、入滲等因素,更真實地模擬實際洪水過程。Parsapour等[4]分別采用Delft3D-Flow和MIKE 21進行水動力建模,在同等網(wǎng)格分辨率的情況下,MIKE 21會獲得更準(zhǔn)確的模擬結(jié)果。MIKE 21采用的二維水流模型基本方程組包括二維流水流連續(xù)方程與水流動量方程,其中二維流水流連續(xù)方程見式(1):

(1)

式中ξ——水面到基準(zhǔn)面的高度,m;t——計算時間,s;u——流速在x方向上的分量,m/s;v——流速在y方向上的分量,m/s。

x和y方向上的二維水流動量方程見式(2)、(3):

(2)

(3)

式中h——靜止水深,m;g——重力加速度,m2/s;V——風(fēng)速;Vx——x方向的風(fēng)速矢量,m/s;Vy——y方向的風(fēng)速矢量,m/s;C——阻力系數(shù);f——風(fēng)摩擦因素,f=γa2ρa;γa——風(fēng)應(yīng)力系數(shù);ρa——空氣密度;Ω——科氏力系數(shù),與緯度相關(guān)。

2 潰堤洪水災(zāi)害損失評估方法

基于研究區(qū)域洪災(zāi)承災(zāi)體類別、價值和空間分布情況,構(gòu)建空間人群分布、經(jīng)濟價值量化、生態(tài)系統(tǒng)服務(wù)價值分布數(shù)據(jù)空間化拓?fù)鋽?shù)據(jù)庫,發(fā)展可融合洪水時空分布數(shù)據(jù)的災(zāi)害損失評估方法,并結(jié)合洪水淹沒的時空分布特征,動態(tài)評估潰堤洪水所造成的生命、經(jīng)濟和生態(tài)環(huán)境損失。對于生命損失,一方面需要考慮人群的空間分布數(shù)據(jù),另一方面還需要考慮洪災(zāi)的生命損失率,采用的生命損失LOL計算見式(4)[12]:

(4)

式中 PAR——風(fēng)險人口數(shù)量,即研究區(qū)域內(nèi)居民受洪水影響的總數(shù)量;a——人群分布數(shù)據(jù)的像元網(wǎng)格總數(shù)量;IRi——第i個像元網(wǎng)格上的風(fēng)險個體生命損失率,計算見式(5)、(6)[12]:

IR=f0αβ

(5)

α=qm1+(1-q)m2

(6)

式中f0——中國風(fēng)險人口死亡率;α——洪災(zāi)嚴(yán)重程度系數(shù);β——修正系數(shù),一般取值為1.4;m1——災(zāi)難程度主要影響因子;m2——災(zāi)難程度次要影響因子;q——權(quán)重系數(shù),取值為0.8。

其中,m1可視為各個主要影響因素對風(fēng)險人口死亡率的影響程度與其權(quán)重系數(shù)的乘積,見式(7)[13]:

(7)

式中si——第i個主要影響因素對風(fēng)險人口死亡率的影響程度,主要包含風(fēng)險人口數(shù)量、潰堤洪水嚴(yán)重性SD、預(yù)警時間WT以及對洪水嚴(yán)重性的理解程度UD;θi——權(quán)重系數(shù),分別取值為θ1=θ2=0.2,θ3=θ4=0.3。

m2為各個次要影響因素對風(fēng)險人口死亡率的影響程度與其權(quán)重系數(shù)的乘積,見式(8)[13]:

(8)

式中ti——第i個次要影響因素對風(fēng)險人口死亡率的影響程度;ni——ti對應(yīng)的權(quán)重系數(shù)值。

經(jīng)濟損失評估主要包括直接經(jīng)濟損失和間接經(jīng)濟損失兩方面。其中,直接經(jīng)濟損失采取損失率方法計算,通過建立各類財產(chǎn)洪水損失率,結(jié)合洪水淹沒特征、經(jīng)濟損失統(tǒng)計指標(biāo)以及淹沒水深等級,進行疊加計算見式(9)[14]:

(9)

式中D——總直接經(jīng)濟損失;Sp——第p類財產(chǎn)的直接經(jīng)濟損失;e——研究區(qū)域財產(chǎn)分布在二維空間上的像元網(wǎng)格數(shù)量;h——第j個像元網(wǎng)格上的財產(chǎn)種類;l——不同淹沒水深等級;βijk——第i個像元網(wǎng)格上第j類財產(chǎn)在第k種淹沒水深對應(yīng)的損失率;Wijk——第i個像元網(wǎng)格上第j類財產(chǎn)對應(yīng)的財產(chǎn)價值。

目前間接經(jīng)濟損失R主要通過洪災(zāi)間接系數(shù)進行確定,根據(jù)直接經(jīng)濟損失與間接經(jīng)濟損失的關(guān)系計算見式(10)[14]:

(10)

式中Ki——第p類財產(chǎn)對應(yīng)的間接系數(shù),農(nóng)業(yè)取值為15%～28%,住宅區(qū)取值為15%,公路取值為25%。

對于生態(tài)環(huán)境損失計算,改進了文獻(xiàn)[15]方法,見式(11):

E=c×ESV

(11)

式中c——經(jīng)驗損失率,取c為0.1;ESV——被洪水淹沒區(qū)域的生態(tài)系統(tǒng)服務(wù)價值;E——按經(jīng)驗系數(shù)法計算的生態(tài)環(huán)境損失。

3 洪災(zāi)損失評估系統(tǒng)實現(xiàn)

為定量評估潰堤洪水災(zāi)害損失,借助面向?qū)ο蟮挠嬎銠C程序語言——Python,調(diào)用Arcpy站點包實現(xiàn)內(nèi)嵌于ArcGIS的內(nèi)置插件條,基于堤防保護區(qū)域人群空間分布數(shù)據(jù)、經(jīng)濟財產(chǎn)價值量化空間分布數(shù)據(jù)和生態(tài)系統(tǒng)服務(wù)價值空間分布數(shù)據(jù),有機結(jié)合洪水時空分布特征,建立洪水災(zāi)害損失動態(tài)評估模型,實現(xiàn)對生命、經(jīng)濟以及生態(tài)環(huán)境損失的動態(tài)定量評估。

在損失評估模型的實現(xiàn)過程中,洪水分布數(shù)據(jù)和承載體空間分布數(shù)據(jù)的類型及精度是不一致的。在后續(xù)的損失計算過程中,要求數(shù)據(jù)可以相互疊加計算,需要保證洪水?dāng)?shù)據(jù)和承載體數(shù)據(jù)的投影坐標(biāo)系、柵格像元大小以及數(shù)據(jù)類型等一致,所以需要提前進行數(shù)據(jù)柵格標(biāo)準(zhǔn)化處理。具體步驟如下:①根據(jù)潰堤洪水演進模型,提取未被標(biāo)準(zhǔn)化的洪水?dāng)?shù)據(jù);②將洪水?dāng)?shù)據(jù)進行面轉(zhuǎn)柵格處理,然后將柵格化的洪水?dāng)?shù)據(jù)與承載體數(shù)據(jù)空間分布數(shù)據(jù)進一步進行標(biāo)準(zhǔn)化處理。

對于生命損失評估模塊,在數(shù)據(jù)標(biāo)準(zhǔn)化處理模塊的基礎(chǔ)上,實現(xiàn)動態(tài)生命損失評估的具體步驟如下:①輸入洪水的時間尺度,本文輸入時間是潰堤洪水演進過程中的模型起始時間和末尾時間,設(shè)定的間隔時間為2 h,生命損失模塊可根據(jù)間隔時間循環(huán)讀取淹沒水深數(shù)據(jù)和流速數(shù)據(jù);②設(shè)定模型參數(shù),計算個體風(fēng)險生命損失率IR分布情況,累加計算所有像元網(wǎng)格的個體風(fēng)險生命損失率總值,疊加計算人群數(shù)據(jù),將獲得的生命損失值統(tǒng)計情況存入表格中。

對于經(jīng)濟損失評估模塊,洪水淹沒是1個動態(tài)變化的過程,某一瞬時的洪水?dāng)?shù)據(jù)僅能分析這個時間點的經(jīng)濟損失空間分布情況。經(jīng)濟損失評估具體步驟如下:①為實現(xiàn)經(jīng)濟損失動態(tài)計算,輸入潰堤洪水演進過程中起始、末尾時間以及間隔時間,表示每過一次間隔時間計算一次經(jīng)濟損失;②輸入經(jīng)濟分布數(shù)據(jù)和淹沒水深數(shù)據(jù),而財產(chǎn)種類數(shù)據(jù)和損失率關(guān)系可根據(jù)蓄滯洪區(qū)或者防洪保護區(qū)的實際社會情況進行修改;③輸出各類財產(chǎn)損失的統(tǒng)計表以及經(jīng)濟損失的分布數(shù)據(jù)。

對于生態(tài)環(huán)境損失評估模塊,根據(jù)式(11)可知,需要輸入生態(tài)系統(tǒng)服務(wù)價值空間分布數(shù)據(jù)。生態(tài)環(huán)境損失評估具體步驟如下:①輸入始末以及間隔時間;②通過疊加洪水淹沒的面積篩選受影響的生態(tài)系統(tǒng)服務(wù)價值空間分布數(shù)據(jù),設(shè)定經(jīng)驗損失率,計算生態(tài)損失;③輸出生態(tài)損失統(tǒng)計表。

本文建立的洪災(zāi)損失動態(tài)評估模型有效利用了潰堤洪水演進模擬獲取的洪水空間分布數(shù)據(jù),其主要優(yōu)勢如下:①基于ArcGIS平臺構(gòu)建空間化拓?fù)鋽?shù)據(jù)庫,并且與國家數(shù)據(jù)統(tǒng)計年鑒內(nèi)容緊密結(jié)合,保證了數(shù)據(jù)有效性,同時解決了數(shù)據(jù)的及時更新問題;②考慮了洪水和承載體數(shù)據(jù)的空間分布差異特征,有效提升了洪水損失評估結(jié)果準(zhǔn)確性;③通過有限的輸入和輸出以及模型參數(shù)修改反映了洪災(zāi)損失動態(tài)變化情況;④可以根據(jù)研究對象的洪災(zāi)損失典型資料,調(diào)查分析和調(diào)整洪災(zāi)損失率關(guān)系,兼顧了模型的普適性;⑤實現(xiàn)了潰堤洪水演進模擬、空間化拓?fù)鋽?shù)據(jù)庫動態(tài)連接以及洪災(zāi)損失的實時評估。

4 工程應(yīng)用

鄱陽湖區(qū)某重點堤防地處江西省九江市永修縣境內(nèi),圩堤起于永修縣城郊的觀音巖,止于楊公腦,全長33.57 km,保護面積56.28 km2,保護耕地5.03×104畝(1畝約等于666.67 m2,下同),常住人口2.34萬人。據(jù)2005年現(xiàn)場勘察,三角聯(lián)圩各類險工險段多達(dá)6處,累計堤線長30多km,大多堤段的斷面較為單薄且內(nèi)外坡均較陡,汛期水流逐年沖刷,極易導(dǎo)致堤身滲透破壞,存在部分堤段填筑壓實不均勻且高程不足,脫坡險情突出。2020年7月12日19時40分,潰口寬度達(dá)20余m,潰堤后,圩內(nèi)平均水位在0.5 h內(nèi)上漲0.5 m多,由于洪水沖刷,潰口迅速擴大,截至7月13日10點,三角聯(lián)圩潰口寬度已經(jīng)擴大到200余m,并于7月14日14時啟動封堵作業(yè),7月16日21時43分,三角聯(lián)圩完成合攏。淹沒耕地面積5.54×104畝,水產(chǎn)業(yè)面積1.05×104畝,淹沒5 100余棟房屋,其中倒塌了40余棟,潰口經(jīng)過55 h的搶修后合攏。

4.1 模型計算范圍確定

基于中國大地坐標(biāo)系,建立堤防保護區(qū)域數(shù)字高程模型見圖1。采用三角形網(wǎng)格進行模型剖分見圖2,將東南部堤防潰口(紅色邊界)設(shè)置為入流邊界,北部三角聯(lián)圩和西部陸地(綠色邊界)設(shè)置為閉邊界。

圖1 堤防保護區(qū)域數(shù)字高程模型

圖2 堤防水動力學(xué)模型網(wǎng)格剖分

本數(shù)值模型共存在4個片區(qū),分別為堤防南部片區(qū)、北部片區(qū)、東部片區(qū)和西部片區(qū),在每個片區(qū)內(nèi)分別設(shè)置4個特征點(即自然村)。假設(shè)決口發(fā)生在三角聯(lián)圩南部片區(qū)下方見圖3,其具體地理位置為九江市永修縣三角鄉(xiāng)、永豐墾殖場和南昌市新建區(qū)大塘坪鄉(xiāng)交界,樁號27k+880～28k+010 m處。

圖3 三角堤防區(qū)域特征點和潰口位置

4.2 模型參數(shù)設(shè)置及有效性

提出的洪水演進模擬方法涉及的物理參數(shù)(包括邊界條件、網(wǎng)格形式、初始條件、渦黏系數(shù))見表1。

表1 洪水演進數(shù)值模型主要物理參數(shù)取值[16]

本模型中土地利用類型主要包括旱地、水田、草地、灘地、居民區(qū)和水面6種,其相應(yīng)糙率參數(shù)按照下墊面分布情況,并參照相關(guān)經(jīng)驗進行率定選取[19],見表2。為了驗證模型的有效性,圖4比較了采用所建模型模擬的潰口淹沒水位過程線和實測的潰口區(qū)域內(nèi)淹沒水位過程線。由圖4可知,在水位基本不變后,水位模擬值為22.65 m,水位實測值為22.58 m,兩條水位過程線較為吻合,表明基于MIKE 21軟件能夠合理模擬潰堤洪水在堤防保護區(qū)域內(nèi)的演進過程。此外,所建模型模擬的淹沒面積為56.15 km2,與堤防潰口合攏后的實測淹沒面積(56.28 km2)相差僅0.13 km2,進一步說明了本模型的有效性,可以滿足計算精度要求。

表2 糙率參數(shù)參考取值

圖4 潰口區(qū)域內(nèi)模擬水位與實測水位的比較

4.3 潰堤洪水演進過程分析

圖5給出了三角聯(lián)圩潰口流量隨潰決時間的變化情況。2020年7月12日19時為圖5中的起始時刻,在7月12日19時40分堤防潰決,隨后堤防潰口寬度逐漸擴大至最大值200 m,潰堤后的21 h內(nèi)流量從0急劇增加到最大值1 724 m3/s,隨后流量開始下降,于7月14日14時啟動封堵作業(yè),7月16日21時43分,三角聯(lián)圩潰口完成合攏,流量接近于零,整個堤防保護區(qū)域進洪過程在98 h內(nèi)完成。

圖5 三角堤防潰口流量過程線

為了進一步展示潰堤洪水演進的淹沒過程,圖6給出了6個典型進洪時段的淹沒水深分布。洪水在決堤后10 h內(nèi)到達(dá)堤防南部片區(qū),雖然紅旗村處于潰口附近,但是因紅旗村的地勢較高,在潰決40 h內(nèi)受災(zāi)情況較輕。由于堤防北部片區(qū)地勢較低,洪水在20 h內(nèi)向周家湖、三角圩等村莊挺進,潰堤40 h后,堤防北部片區(qū)基本被淹沒,洪水逐漸向東西部蔓延,堤防保護區(qū)域內(nèi)水深高達(dá)5.2 m。潰堤60 h后,堤防北部片區(qū)淹沒水深最高達(dá)6.0 m;潰堤80 h后,因堤防西部片區(qū)比東部片區(qū)地勢較高,故東部區(qū)域淹沒深度最深達(dá)到6.5 m,而西部區(qū)域仍舊有部分地區(qū)未被淹沒。進洪100 h后,堤防保護區(qū)域絕大部分被淹沒,此時需要盡快對圩內(nèi)居民實施避險轉(zhuǎn)移搶險救援等舉措。

a)10 h

4.4 潰堤洪水損失評估

根據(jù)收集的洪水影響對象統(tǒng)計數(shù)據(jù),對人群數(shù)據(jù)進行空間展布,并量化經(jīng)濟財產(chǎn)和生態(tài)系統(tǒng)服務(wù)價值。然后,搭建空間化拓?fù)鋽?shù)據(jù)庫來反映生命、經(jīng)濟、生態(tài)環(huán)境在空間上的分布差異,接著在潰堤洪水演進模型中提取洪水淹沒時空分布數(shù)據(jù),并將數(shù)據(jù)柵格進行標(biāo)準(zhǔn)化處理。最后,進行包括生命、經(jīng)濟和生態(tài)環(huán)境損失在內(nèi)的綜合損失動態(tài)評估。圖7給出了潰堤造成的生命損失隨潰口歷時的變化曲線。由圖7可知,隨著預(yù)警時間的延長,生命損失不斷減少。同時,居民對洪水嚴(yán)重性的理解程度也對生命損失具有重要的影響因素。例如,在預(yù)警時間不超過0.25 h且居民對洪水嚴(yán)重性理解程度為模糊的情況下,生命損失高達(dá)454人,見圖7a;相較之下,如果預(yù)警時間足夠長,且居民對洪水嚴(yán)重性的理解程度為明確的情況下,生命損失驟降至2人,見圖7b。

a)居民模糊理解洪水嚴(yán)重性情況下生命損失曲線

圖8給出了潰堤造成的經(jīng)濟損失隨潰口歷時的變化關(guān)系曲線。由圖8可知,直接經(jīng)濟損失隨著潰口歷時的增加而增大,其中居民住房直接經(jīng)濟損失最大,高達(dá)2.86億元,其次為家庭財產(chǎn)直接經(jīng)濟損失。同時,圖9給出了堤防保護區(qū)域內(nèi)總經(jīng)濟損失分布情況,不難看出居民用地?fù)p失情況最為嚴(yán)重。

a)耕地-家庭財產(chǎn)-住房直接經(jīng)濟損失

圖9 潰堤造成的總經(jīng)濟損失分布

5 結(jié)論

以江西省鄱陽湖區(qū)某重點堤防為例,建立了基于MIKE 21軟件的洪水演進數(shù)值模型,獲得了潰堤洪水下淹沒水深和流速等洪水空間分布數(shù)據(jù),豐富了潰堤洪水演進模擬方法。同時,考慮生命、經(jīng)濟、生態(tài)環(huán)境的空間分布特性,基于Python平臺建立了一種具有明確物理意義的,可考慮生命、經(jīng)濟和環(huán)境損失的定量動態(tài)評估模型。主要結(jié)論如下。

a)根據(jù)研究區(qū)下墊面土地利用情況設(shè)置不同的糙率值,以及與實測的潰口處淹沒水位過程線進行對比,驗證了洪水演進數(shù)值模型的有效性,表明所建模型具有較高的計算精度,可有效用于獲取淹沒水深、流速等洪水空間分布數(shù)據(jù)。

b)發(fā)展了潰堤洪水災(zāi)害損失評估方法,并采用Python語言編寫了基于ArcGIS的嵌入式耦合插件工具,建立了洪災(zāi)損失動態(tài)評估模型,可以有效評估不同工況下潰堤洪水造成的生命、經(jīng)濟以及生態(tài)環(huán)境損失,為洪水災(zāi)害防控提供了有效途徑。