基于多重現(xiàn)期洪水情境模擬和綜合風險度方法的入海河口地區(qū)洪水風險區(qū)劃研究

王秀杰，王艷鵬，苑希民，楊 航

(1.天津大學水利工程仿真與安全國家重點實驗室，天津 300350； 2.中國水利水電第五工程局，四川 成都 610066)

近年來，受全球變暖和氣候變化影響，我國暴雨、洪水及風暴潮成災事件頻發(fā)，防汛形勢日趨嚴峻。我國三角洲地區(qū)河網(wǎng)縱橫，地勢平坦，遭遇洪澇災害時洪水難以外排，且下游風暴潮沿河道上溯、頂托，導致上游洪水下泄不暢，容易形成較大的潰漫堤風險。因此，開展三角洲地區(qū)的多源洪水風險區(qū)劃對于指導流域防洪規(guī)劃和防洪減災非工程措施建設具有重要意義。

國內外學者針對洪水風險區(qū)劃進行了大量研究，美國學者[1]根據(jù)河流歷史洪水發(fā)生情況，計算了5年、10年和20年一遇洪水條件下河流兩側不同洪水風險緩沖區(qū)，Plate等[2]選用洪水頻率指標確定了受山洪威脅區(qū)域的禁止開發(fā)區(qū)，王艷艷等[3]從洪災概率、風險指標和防洪標準等方面整合了中國和荷蘭現(xiàn)有風險評估方法，并以安慶市為例分析計算了淹沒區(qū)的洪水風險，此類洪水區(qū)劃方法依據(jù)特定頻率洪水淹沒范圍進行洪水風險劃分，但對于同一頻率洪水淹沒區(qū)域，洪水淹沒特征不同，災害影響也不同，因此，對于洪水致災特征的風險區(qū)劃同樣重要。英國洪災研究中心(FHRC)Penning-Roswell等[4]研究了基于確定性框架的建筑物脆弱性與洪水淹沒水深及淹沒歷時的關系，劉樹坤教授[5]在二維洪水淹沒分析的基礎上，根據(jù)水深和流速對生命財產(chǎn)損失程度分為極危險區(qū)、重災區(qū)、輕災區(qū)和安全區(qū)，澳大利亞和尼泊爾洪泛區(qū)管理委員會[6]考慮洪水動量因素引起的洪水致災風險，基于水深流速之積劃分洪水危害級別，苑希民等[7-8]從致災因子、孕災環(huán)境和承災體三個方面分析遴選評價指標，基于模糊層次分析法建立了京津冀地區(qū)洪災風險綜合評價模型并實現(xiàn)了風險等級的區(qū)劃，羅賢等[9]利用二維Gumbel模型，開展了針對秦淮河流域暴雨與洪水水位兩變量的洪水風險分析，蔣衛(wèi)國、方建等[10, 11]基于系統(tǒng)理論從指標選取原則、指標量化、風險等級等方面構建了洪災風險評估指標體系，實現(xiàn)了區(qū)域洪水的災害風險評估，郭偉萍等[12]將數(shù)學期望概念引入洪災損失評估，提出基于潰堤洪水數(shù)值模擬的洪災損失期望值計算方法，確定洪水風險附加費率。以上可歸類為洪水致災特征區(qū)劃，此類方法通?？紤]單一洪水頻率，針對單一風險特征指標或者綜合風險指標進行洪水風險區(qū)劃，考慮的洪水風險因子較少，不能全面地反映多重現(xiàn)期的真實洪水風險。

基于此，本文考慮最大淹沒水深，最大行進流速，最大淹沒歷時洪水風險要素指標，提出了融合不同洪水頻率的各洪水風險特征的區(qū)劃方法，基于凸風險度量理論構建洪水頻率和洪水風險特征曲線，利用期望理論計算多重現(xiàn)期洪水的綜合風險度，以江新聯(lián)圍防洪保護區(qū)為例開展洪水風險區(qū)劃研究。

1 洪水風險分析模型與區(qū)劃方法

1.1 洪水風險分析模型

本文基于水力學理論進行洪水風險分析，主要包括一維水動力學方法、二維水動力學方法及一、二維耦合水動力學方法。理論方法介紹如下：

(1)一維水動力模型

一維水動力學模型采用的基本方程如下：

(1)

(2)

式中：Q為流量；A為斷面面積；x為河道沿程坐標；t為時間；q為旁側來水量；α為修正系數(shù)；g為重力加速度；h為水位；C為謝才系數(shù)；R為水力半徑。

(2)二維水動力模型

二維水動力計算模型基本原理如下：

連續(xù)方程：

(3)

動量方程：

(4)

(5)

式中：H為水深；Z為水位；q為連續(xù)方程中的源匯項；M與N分別為x和y方向的垂向平均單寬流量；u和v分別為垂向平均流速在x與y方向的分量；n為曼寧糙率系數(shù)；g為重力加速度。

(3)一、二維水動力耦合模型

堤防發(fā)生漫溢和潰決后，河道洪水由堤防和潰口進入保護區(qū)，為準確描述一維河道和二維保護區(qū)的水流相互作用，實現(xiàn)一、二維模型的時空動態(tài)耦合，本文構建了河道-保護區(qū)的一、二維耦合模型，耦合處采用寬頂堰公式進行側向連接，以實現(xiàn)耦合處水流信息的實時交互。寬頂堰流公式表述如下：

Q=μb(h1-Zc)(h1-h2)1/2.

(6)

式中：Q為過堰流量；b為堰頂寬度；h1和h2分別為上下游水位；Zc為堰頂高程；μ為流量系數(shù)。

一維模型為二維模型提供水位及流量邊界，二維模型則以旁側入流的方式將邊界流量傳遞給一維模型，河道和洪泛區(qū)之間水量不斷交換，達到動態(tài)平衡。其中，針對一維非恒定流方程組，采用六點隱式Abbott差分格式對模型進行離散求解[13]，對于二維水動力學基本方程組，采用單元中心的顯式有限體積法離散求解[14]，從而確保水流滿足質量守恒定律和動量守恒定律[15]。

1.2 洪水風險區(qū)劃方法

本文基于凸風險度量理論構建洪水頻率和洪水風險特征曲線，利用期望理論計算多個量級洪水的綜合風險度，理論方法介紹如下：

(1)凸風險度量理論

2002年，F(xiàn)ollmer和Schied[16-17]提出凸風險度量公理，定義如下：根據(jù)給定概率空間(Ω,F,P)和此概率空間上的一個實值隨機變量X，它代表某個資產(chǎn)或投資組合未來的不確定性收益或現(xiàn)金流，在洪水風險區(qū)劃中它代表可能會發(fā)生的某一量級的洪水頻率或洪水風險程度，令U為隨機變量X的所有情景集的集合，一個風險度量實質上是從U到實數(shù)集R的一個映射關系ρ:U→R， 凸風險度量是指其對應的映射ρ滿足如下的三條性質[18-22]：

(A)轉移不變性：ρ(X+a)=ρ(X)-a,?X∈U,a∈R；

(B)凸性：ρ(λX+(1-λ)Y)≤λρ(X)+(1-λ)ρ(Y), ?X,Y∈U,λ∈[0,1]；

(C)單調性：X≤Y?ρ(X)≥ρ(Y), ?X,Y∈U。

(2)基于凸風險度量和期望理論的綜合風險度法

針對入?？诘暮闈吵睆碗s工況，選取最大淹沒水深、最大行進流速、淹沒歷時指標作為評價此地區(qū)的洪水風險因子是合適的，因此考慮指標融合的洪水風險指標H；同時洪水風險指標H又隨著洪水頻率凸性遞減，因此借鑒凸風險度量表示性定理，構建洪水頻率和洪水風險特征曲線，如圖1所示，在洪水頻率和洪水風險特征曲線圖中，單調性和轉移不變性是指洪水頻率越低，相應的洪水風險越大；凸性是指隨著洪水頻率的增加，對應的洪水風險并不是線性變化，而呈凸性遞減。洪水頻率和洪水風險特征滿足凸風險度量公理。為反映多個量級洪水綜合淹沒情況下洪水風險的空間分布特征以及區(qū)域間洪水風險程度的差異性，本文結合期望理論提出了基于多重現(xiàn)期與風險特征的綜合風險度期望公式如公式(7)所示，并基于綜合風險度劃分洪水風險等級，如表1所示。

圖1 綜合風險度計算示意圖(陰影部分面積即為R)Fig.1 Schematic diagram of comprehensive risk calculation(the shaded area is R)

表1 風險區(qū)化等級劃定范圍Table 1 Grade delineation range of risk zoning

(7)

式中：pi為某一洪水淹沒頻率，Hi為該計算單元對應pi的洪水風險指標H值。由于利用上述公式計算期望值時，計算單元的洪水淹沒指標值Hi在起淹洪水頻率處存在跳躍，故假定在計算時p0始終為起淹洪水頻率的下一級洪水頻率，且對應的H0=0；而p1，pn則分別為該計算單元的起淹洪水頻率和最高洪水計算頻率。

其中，洪水風險指標(H)表征計算單元在某一量級洪水頻率下的風險程度大小，指標選取應以能全面反映洪水淹沒特征為要，因此指標計算以“最大淹沒水深”為主要因子，綜合考慮“最大行進流速”、“最大淹沒歷時”風險要素的影響，公式如下：

H=α1α2h.

(8)

式中：α1為“最大行進流速”修正系數(shù)，α2為“最大淹沒歷時”修正系數(shù)[23]。其中，當v<1.5 m/s時，α1=1.0，1.5 m/s≤v<3.0 m/s時，α1=1.2,v≥3.0 m/s時，α1=1.5；當t<3d(天)時，α2=1.0，3 d≤t<7 d時，α2=1.2，t≥7d時，α2=1.5。

2 研究區(qū)概況與資料

2.1 研究區(qū)概況

廣東江新聯(lián)圍防洪保護區(qū)屬于珠江流域，地處珠江三角洲下游，涉及江門市蓬江區(qū)、江海區(qū)、新會區(qū)等行政區(qū)，區(qū)內河網(wǎng)縱橫，水系復雜，河道眾多，周邊水系涉及西江干流、磨刀門水道、潭江、虎跳門水道等50條河流，河道總長度為1 016.3 km，江新聯(lián)圍堤線北起天河頂，沿西江自上而下到達潭江左岸的梅林沖，干堤全長91.76 km，防洪保護區(qū)內水情復雜，洪源眾多，不僅受暴雨、洪水相互影響，而且汛期受下游潮位頂托作用，區(qū)域內暴雨積水難以外排，導致洪澇潮災害頻發(fā)。

2.2 保護區(qū)內洪水來源

根據(jù)珠江三角洲洪潮水文特性和洪潮影響程度，本文在洪水風險分析的基礎上確定江新聯(lián)圍防洪保護區(qū)內的洪水來源主要有：

(1)河道洪水

保護區(qū)內洪水來源主要為西江洪水和潭江洪水，洪水具有峰高、量大、歷時長、漲落較緩慢等特點，歷史洪水災害損失嚴重。本文在模型中以1998年6月大洪水作為典型，采用同頻率放大法推求西江高要站設計洪水過程線，并以思賢滘設計洪峰流量作為控制條件，考慮石角站洪峰流量與高要站洪峰流量在思賢滘遭遇，對石角站1998年典型洪水最大7天過程進行同倍比縮放并平移，使得石角站洪峰流量與高要站洪峰流量相碰。圖2為高要站和石角站相應設計洪水過程線。

(2)風暴潮

保護區(qū)入海河口支流較多，水流分散，洪潮摻雜，受潮汐作用和臺風登陸影響，歷史上發(fā)生多次風暴潮災害。保護區(qū)內下游潮位站有老鴨崗、石咀、西炮臺、官沖、黃金、萬頃沙西、燈籠山、橫門、焦門和三沙口，潮位過程以2008年“0814”號“黑格比”臺風風暴潮最為典型，本文采用同倍比放大法放大各站點典型洪水過程，計算得到各頻率的潮位過程。各潮位站潮位邊界條件如圖2所示。

圖2 各水文站設計洪水和潮位過程Fig.2 Design flood and tide level process of each hydrological station

(3)內澇洪水

江新聯(lián)圍內暴雨洪水主要來自上游的天沙河、會城河，暴雨常伴隨河道洪水或風暴潮產(chǎn)生，由于外江水位高，圍內漬水不能自排流，容易形成內澇災害。河道附近的雨量氣象站有天河、江門、新會、大敖站，其中新會站資料系列最長最完整，位于江新聯(lián)圍中部，能較好地反映圍內洪水情況，4個雨量站中新會站設計暴雨值最大，從最大洪水風險角度考慮，本次采用新會氣象站作為雨量代表站。新會站72 h設計暴雨成果如圖3所示。

圖3 新會站50年一遇設計暴雨過程Fig.3 Design rainstorm process of Xinhui station once in fifty years

2.3 洪水計算情境擬定

江新聯(lián)圍防洪保護區(qū)處于洪水、潮水綜合影響區(qū)域，洪水來源為西江洪水、潭江洪水、外海風暴潮、暴雨內澇。保護區(qū)堤防屬2級堤防，防洪標準為50年一遇，經(jīng)過堤防整體防洪能力復核后確定防洪標準為20年一遇。通過西江與潭江段河道洪水與外海風暴潮遭遇相關性分析，考慮可能遭遇的最大風險，并根據(jù)現(xiàn)場調查和歷史資料分析，保護區(qū)遭遇外海風暴潮時常伴有降雨發(fā)生，因此考慮洪水、風暴潮與降雨的組合；同時由于以往研究大多針對單一洪水頻率工況，無法體現(xiàn)多重現(xiàn)期洪水風險程度的聯(lián)系與差異性，因此本文考慮并制定多重現(xiàn)期洪水組合情境方案，具體情況如表2所示。

表2 廣東江新聯(lián)圍防洪保護區(qū)洪水區(qū)劃計算方案Table 2 Flood zoning calculation scheme of Jiangxinlianwei flood protection area in Guangdong Province

潰口設置方面考慮現(xiàn)狀設防標準高一等級及以上頻率洪水的堤防漫潰，對于漫溢段堤防經(jīng)論證不易形成漫溢潰口的，只設漫溢段，不設潰口。漫溢潰口位置及潰決時間點根據(jù)洪水水位變化與堤防高程的相互關系進行試算確定。

3 結果與討論

3.1 洪水風險因素分布

為研究江新聯(lián)圍防洪保護區(qū)洪水風險因素分布特征，需根據(jù)洪水分析情景計算方案建立多源洪水一、二維水動力耦合模型。模型構建步驟如下：(1)采用1:5萬河道一維矢量數(shù)據(jù)和珠江三角洲2008年實測斷面資料建立一維河網(wǎng)，建模范圍為西江高要站以下、北江石角站以下、潭江石咀站以下及珠江三角洲部分，共包括50條河流，涉及河道總長度為1 016.3 km；(2)二維模型計算范圍在保護區(qū)原有基礎上向西適當延展，面積達437.6 km2，為充分反映保護區(qū)內堤防、公路、鐵路、河道等線狀地物的位置等特征，考慮線狀地物的阻水作用，應用內邊界對線狀地物進行概化處理，并根據(jù)保護區(qū)內土地利用情況，對保護區(qū)內的村莊、道路、耕地、河流等地物設置了糙率分區(qū)，以反映保護區(qū)不同下墊面對洪水演進的影響，較為真實地還原洪水演進過程和水流交換效果，線狀地物和糙率分區(qū)分別如圖4(b)(c)所示；(3)基于河道一維與保護區(qū)二維水動力模型，疊加河道上下游洪潮邊界及區(qū)域內降雨過程，通過潰口與漫溢段實現(xiàn)河道與保護區(qū)的耦合連接及水流信息的實時交互，以此分析洪、澇、潮多源洪水在防洪保護區(qū)的淹沒情況。一、二維耦合模型如圖4(a)所示。

圖4 一、二維耦合模型構建Fig.4 Construction of one-dimensional and two-dimensional coupling models

為了保證模型計算精度要求及結果合理性，采用三水、馬口、北街、三江口等24個水文站點的實測潮位數(shù)據(jù)對模型進行率定和驗證，其中與二維保護區(qū)臨近的北街站和三江口站的率定結果與驗證結果如圖5所示。由圖可知，各水文站洪水過程計算結果與實測數(shù)據(jù)吻合；峰值、峰現(xiàn)時間接近，率定和驗證效果較好，模型能夠較為準確地模擬多源洪水條件下的水流變化特征，滿足計算精度要求。

圖5 各水文站率定驗證結果Fig.5 Calibration and verification results of various hydrological stations

利用率定與驗證后的模型參數(shù)建立一維河道與二維保護區(qū)的多源洪水耦合模型，對多重現(xiàn)期洪水情境方案依次進行模擬計算，以“西江、潭江100年一遇洪水、50年一遇內澇、50年一遇設計潮位”計算方案為例，提取研究區(qū)內各計算單元的風險要素指標進行洪水風險要素分析，包括最大淹沒水深、最大行進流速、最大淹沒歷時，如圖6所示。

由圖4和圖6對比可知，最大淹沒水深較多出現(xiàn)在地形較低洼的地區(qū)以及潰口附近，受線狀地物的阻水作用，道路和堤防兩側的淹沒水深略有差異，最大洪水流速較多出現(xiàn)在地形坡度較大的地區(qū)，針對區(qū)域暴雨內澇情況，保護區(qū)內除局部線狀地物淹沒歷時較小外，整個區(qū)域淹沒歷時均較大，參考以往研究中的洪水風險空間分布特征合理性分析與研究區(qū)域有關洪水風險分析成果[24-28]，綜合分析可知，本次基于多源洪水耦合模型獲得的洪水風險要素分布較為合理，結果能較為準確地反映保護區(qū)的淹沒特征以及不同重現(xiàn)期洪水風險程度的差異性，模型的精度和可靠性較高，可為后續(xù)洪水風險區(qū)劃提供詳實可靠的數(shù)據(jù)支撐。

圖6 洪水風險要素指標Fig.6 Flood risk factor index

3.2 區(qū)劃結果分析

以“綜合風險度(R)”為指標，以“最大淹沒水深”為主要影響因子，綜合考慮“最大行進流速”、“最大淹沒歷時”的影響，疊加計算多重現(xiàn)期洪水的3種計算方案的洪水風險要素特征，根據(jù)計算結果選取具有代表性的風險區(qū)劃計算單元，繪制基于凸風險度量和期望理論的風險特征曲線，如圖7所示；并根據(jù)綜合風險度法區(qū)劃等級劃定范圍對江新聯(lián)圍防洪保護區(qū)進行洪水風險等級區(qū)劃，區(qū)劃結果如圖8所示。

圖7 基于凸風險度量和期望理論的風險特征曲線Fig.7 Risk characteristic curves based on convex risk measure and expectation theory

為了檢驗綜合風險度方法的可靠性，本文采用澳大利亞和尼泊爾洪泛區(qū)管理委員會[5]提出的洪水風險區(qū)劃方法進行驗證，該方法考慮洪水動量因素引起的洪水致災風險，運用水深流速之積來劃分洪水危害級別，區(qū)劃標準參考《Designing Safer Subdivisions-Guidance On Subdivision Design In Flood Prone Areas》中風險等級劃分方法，如表3所示，區(qū)劃結果如圖9所示。

表3 最大水深×最大流速洪水風險等級劃定范圍Table 3 The flood risk level range of the maximum water depth multiplied by the maximum flood velocity method

圖8 綜合風險度法風險區(qū)劃結果 圖9 最大水深×最大流速法風險區(qū)劃結果Fig.7 Risk zoning results of comprehensive risk method Fig.8 Risk zoning results of the maximum water depth multiplied by the maximum flood velocity method

分析圖7可知，根據(jù)區(qū)劃計算單元繪制的風險特征曲線滿足凸風險度量公理，對于不同的計算單元，隨著洪水頻率的增加，各洪水風險指標值呈現(xiàn)連續(xù)性、同質化遞減趨勢，這表明融合淹沒水深、洪水流速、淹沒歷時的綜合指標值能夠較為準確地反映洪澇潮組合下的多源洪水風險，同時也驗證了流速和歷時指標區(qū)間范圍劃定的合理性。根據(jù)期望理論可以得知，圖中陰影面積對應區(qū)劃單元的綜合風險度值，代表計算單元的洪水風險，洪水風險指標值增大時，基于凸風險度量和期望理論的特征曲線上移，相應的洪水風險增大。

分析圖8和圖9可知，洪水高風險和極高風險主要分布在棠下鎮(zhèn)、北街街道、滘頭街道、外海街道、禮樂街道以及潰口附近，統(tǒng)計綜合風險度法和最大水深×最大流速法的各風險區(qū)域面積如表4所示，綜合風險度法低風險區(qū)、中風險區(qū)、高風險區(qū)、極高風險區(qū)分別占比44.34%、28.30%、18.39%、8.97%，最大水深×最大流速法(以100洪+50潮+50澇為例)低風險區(qū)、中風險區(qū)、高風險區(qū)、較高風險區(qū)、極高風險區(qū)分別占比34.34%、21.55%、23.17%、14.71%、6.23%，最大水深×最大流速法的區(qū)劃結果相較于綜合風險度法風險值整體偏高，且隨著洪水頻率的降低風險逐漸增加。

表4 不同區(qū)劃方法各風險等級區(qū)域面積統(tǒng)計表Table 4 Statistical table of regional area for different risk levels of different zoning methods 單位：m2

江新聯(lián)圍中部地區(qū)地勢平坦，上游暴雨洪水在此泛濫成災，且受下游潮位頂托作用，積水難以外排，對比分析圖6中洪水風險要素指標，這些地區(qū)對應的淹沒水深和淹沒歷時也較高，綜合風險度法區(qū)劃結果的合理性得以體現(xiàn)；江新聯(lián)圍北部棠下鎮(zhèn)和南部三江鎮(zhèn)、睦洲鎮(zhèn)地形坡度較大，洪水流速較高，因此最大水深×最大流速法的區(qū)劃結果風險較高，然而最大水深×最大流速法未考慮淹沒歷時的影響，強化了流速指標的作用，導致會城街道和禮樂街道的風險值偏高。綜合分析淹沒水深、洪水流速、淹沒歷時和地形高程分布情況，綜合風險法洪水風險區(qū)劃結果較為合理，綜合風險度區(qū)劃方法具有一定的合理性。

4 結論

本文借鑒金融風險管理中凸風險度量和期望理論提出基于綜合風險度法的多源洪水耦合風險區(qū)劃方法，在多源洪水和多重現(xiàn)期情境模擬分析的基礎上綜合考慮洪水頻率、最大淹沒水深、最大行進流速、最大淹沒歷時等風險指標，開展江新聯(lián)圍防洪保護區(qū)的洪水風險區(qū)劃研究。研究結果表明，基于洪水風險分析模型獲得的洪水風險要素分布較為合理，計算結果能較為準確地反映保護區(qū)內的淹沒特征以及不同重現(xiàn)期洪水風險程度的聯(lián)系和差異性，模型的精度和可靠性較高；綜合風險度方法彌補了國內外針對單一頻率洪水風險區(qū)劃的不足與局限性，考慮的洪水風險因子能夠較為準確地反映洪水淹沒特性，區(qū)劃結果反映了多個量級洪水綜合淹沒情況下洪水風險的空間分布特征，與最大水深乘最大流速法相比區(qū)劃結果更具合理性?？傊疚奶岢龅木C合風險度區(qū)劃方法更能全面直觀地反映多重現(xiàn)期的洪水風險空間分布特征，適應于洪、澇、潮多源洪水條件下的入海河口地區(qū)的洪水風險區(qū)劃，研究結果可為洪水風險社會化管理、防災減災規(guī)劃和國土空間規(guī)劃等提供決策依據(jù)和技術支撐。