基于多重现期洪水情境模拟和综合风险度方法的入海河口地区洪水风险区划研究

王秀杰，王艳鹏，苑希民，杨 航

(1.天津大学水利工程仿真与安全国家重点实验室，天津 300350； 2.中国水利水电第五工程局，四川 成都 610066)

近年来，受全球变暖和气候变化影响，我国暴雨、洪水及风暴潮成灾事件频发，防汛形势日趋严峻。我国三角洲地区河网纵横，地势平坦，遭遇洪涝灾害时洪水难以外排，且下游风暴潮沿河道上溯、顶托，导致上游洪水下泄不畅，容易形成较大的溃漫堤风险。因此，开展三角洲地区的多源洪水风险区划对于指导流域防洪规划和防洪减灾非工程措施建设具有重要意义。

国内外学者针对洪水风险区划进行了大量研究，美国学者[1]根据河流历史洪水发生情况，计算了5年、10年和20年一遇洪水条件下河流两侧不同洪水风险缓冲区，Plate等[2]选用洪水频率指标确定了受山洪威胁区域的禁止开发区，王艳艳等[3]从洪灾概率、风险指标和防洪标准等方面整合了中国和荷兰现有风险评估方法，并以安庆市为例分析计算了淹没区的洪水风险，此类洪水区划方法依据特定频率洪水淹没范围进行洪水风险划分，但对于同一频率洪水淹没区域，洪水淹没特征不同，灾害影响也不同，因此，对于洪水致灾特征的风险区划同样重要。英国洪灾研究中心(FHRC)Penning-Roswell等[4]研究了基于确定性框架的建筑物脆弱性与洪水淹没水深及淹没历时的关系，刘树坤教授[5]在二维洪水淹没分析的基础上，根据水深和流速对生命财产损失程度分为极危险区、重灾区、轻灾区和安全区，澳大利亚和尼泊尔洪泛区管理委员会[6]考虑洪水动量因素引起的洪水致灾风险，基于水深流速之积划分洪水危害级别，苑希民等[7-8]从致灾因子、孕灾环境和承灾体三个方面分析遴选评价指标，基于模糊层次分析法建立了京津冀地区洪灾风险综合评价模型并实现了风险等级的区划，罗贤等[9]利用二维Gumbel模型，开展了针对秦淮河流域暴雨与洪水水位两变量的洪水风险分析，蒋卫国、方建等[10, 11]基于系统理论从指标选取原则、指标量化、风险等级等方面构建了洪灾风险评估指标体系，实现了区域洪水的灾害风险评估，郭伟萍等[12]将数学期望概念引入洪灾损失评估，提出基于溃堤洪水数值模拟的洪灾损失期望值计算方法，确定洪水风险附加费率。以上可归类为洪水致灾特征区划，此类方法通常考虑单一洪水频率，针对单一风险特征指标或者综合风险指标进行洪水风险区划，考虑的洪水风险因子较少，不能全面地反映多重现期的真实洪水风险。

基于此，本文考虑最大淹没水深，最大行进流速，最大淹没历时洪水风险要素指标，提出了融合不同洪水频率的各洪水风险特征的区划方法，基于凸风险度量理论构建洪水频率和洪水风险特征曲线，利用期望理论计算多重现期洪水的综合风险度，以江新联围防洪保护区为例开展洪水风险区划研究。

1 洪水风险分析模型与区划方法

1.1 洪水风险分析模型

本文基于水力学理论进行洪水风险分析，主要包括一维水动力学方法、二维水动力学方法及一、二维耦合水动力学方法。理论方法介绍如下：

(1)一维水动力模型

一维水动力学模型采用的基本方程如下：

(1)

(2)

式中：Q为流量；A为断面面积；x为河道沿程坐标；t为时间；q为旁侧来水量；α为修正系数；g为重力加速度；h为水位；C为谢才系数；R为水力半径。

(2)二维水动力模型

二维水动力计算模型基本原理如下：

连续方程：

(3)

动量方程：

(4)

(5)

式中：H为水深；Z为水位；q为连续方程中的源汇项；M与N分别为x和y方向的垂向平均单宽流量；u和v分别为垂向平均流速在x与y方向的分量；n为曼宁糙率系数；g为重力加速度。

(3)一、二维水动力耦合模型

堤防发生漫溢和溃决后，河道洪水由堤防和溃口进入保护区，为准确描述一维河道和二维保护区的水流相互作用，实现一、二维模型的时空动态耦合，本文构建了河道-保护区的一、二维耦合模型，耦合处采用宽顶堰公式进行侧向连接，以实现耦合处水流信息的实时交互。宽顶堰流公式表述如下：

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

(6)

式中：Q为过堰流量；b为堰顶宽度；h1和h2分别为上下游水位；Zc为堰顶高程；μ为流量系数。

一维模型为二维模型提供水位及流量边界，二维模型则以旁侧入流的方式将边界流量传递给一维模型，河道和洪泛区之间水量不断交换，达到动态平衡。其中，针对一维非恒定流方程组，采用六点隐式Abbott差分格式对模型进行离散求解[13]，对于二维水动力学基本方程组，采用单元中心的显式有限体积法离散求解[14]，从而确保水流满足质量守恒定律和动量守恒定律[15]。

1.2 洪水风险区划方法

本文基于凸风险度量理论构建洪水频率和洪水风险特征曲线，利用期望理论计算多个量级洪水的综合风险度，理论方法介绍如下：

(1)凸风险度量理论

2002年，Follmer和Schied[16-17]提出凸风险度量公理，定义如下：根据给定概率空间(Ω,F,P)和此概率空间上的一个实值随机变量X，它代表某个资产或投资组合未来的不确定性收益或现金流，在洪水风险区划中它代表可能会发生的某一量级的洪水频率或洪水风险程度，令U为随机变量X的所有情景集的集合，一个风险度量实质上是从U到实数集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)基于凸风险度量和期望理论的综合风险度法

针对入海口的洪涝潮复杂工况，选取最大淹没水深、最大行进流速、淹没历时指标作为评价此地区的洪水风险因子是合适的，因此考虑指标融合的洪水风险指标H；同时洪水风险指标H又随着洪水频率凸性递减，因此借鉴凸风险度量表示性定理，构建洪水频率和洪水风险特征曲线，如图1所示，在洪水频率和洪水风险特征曲线图中，单调性和转移不变性是指洪水频率越低，相应的洪水风险越大；凸性是指随着洪水频率的增加，对应的洪水风险并不是线性变化，而呈凸性递减。洪水频率和洪水风险特征满足凸风险度量公理。为反映多个量级洪水综合淹没情况下洪水风险的空间分布特征以及区域间洪水风险程度的差异性，本文结合期望理论提出了基于多重现期与风险特征的综合风险度期望公式如公式(7)所示，并基于综合风险度划分洪水风险等级，如表1所示。

图1 综合风险度计算示意图(阴影部分面积即为R)Fig.1 Schematic diagram of comprehensive risk calculation(the shaded area is R)

表1 风险区化等级划定范围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为“最大行进流速”修正系数，α2为“最大淹没历时”修正系数[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 研究区概况与资料

2.1 研究区概况

广东江新联围防洪保护区属于珠江流域，地处珠江三角洲下游，涉及江门市蓬江区、江海区、新会区等行政区，区内河网纵横，水系复杂，河道众多，周边水系涉及西江干流、磨刀门水道、潭江、虎跳门水道等50条河流，河道总长度为1 016.3 km，江新联围堤线北起天河顶，沿西江自上而下到达潭江左岸的梅林冲，干堤全长91.76 km，防洪保护区内水情复杂，洪源众多，不仅受暴雨、洪水相互影响，而且汛期受下游潮位顶托作用，区域内暴雨积水难以外排，导致洪涝潮灾害频发。

2.2 保护区内洪水来源

根据珠江三角洲洪潮水文特性和洪潮影响程度，本文在洪水风险分析的基础上确定江新联围防洪保护区内的洪水来源主要有：

(1)河道洪水

保护区内洪水来源主要为西江洪水和潭江洪水，洪水具有峰高、量大、历时长、涨落较缓慢等特点，历史洪水灾害损失严重。本文在模型中以1998年6月大洪水作为典型，采用同频率放大法推求西江高要站设计洪水过程线，并以思贤滘设计洪峰流量作为控制条件，考虑石角站洪峰流量与高要站洪峰流量在思贤滘遭遇，对石角站1998年典型洪水最大7天过程进行同倍比缩放并平移，使得石角站洪峰流量与高要站洪峰流量相碰。图2为高要站和石角站相应设计洪水过程线。

(2)风暴潮

保护区入海河口支流较多，水流分散，洪潮掺杂，受潮汐作用和台风登陆影响，历史上发生多次风暴潮灾害。保护区内下游潮位站有老鸭岗、石咀、西炮台、官冲、黄金、万顷沙西、灯笼山、横门、焦门和三沙口，潮位过程以2008年“0814”号“黑格比”台风风暴潮最为典型，本文采用同倍比放大法放大各站点典型洪水过程，计算得到各频率的潮位过程。各潮位站潮位边界条件如图2所示。

图2 各水文站设计洪水和潮位过程Fig.2 Design flood and tide level process of each hydrological station

(3)内涝洪水

江新联围内暴雨洪水主要来自上游的天沙河、会城河，暴雨常伴随河道洪水或风暴潮产生，由于外江水位高，围内渍水不能自排流，容易形成内涝灾害。河道附近的雨量气象站有天河、江门、新会、大敖站，其中新会站资料系列最长最完整，位于江新联围中部，能较好地反映围内洪水情况，4个雨量站中新会站设计暴雨值最大，从最大洪水风险角度考虑，本次采用新会气象站作为雨量代表站。新会站72 h设计暴雨成果如图3所示。

图3 新会站50年一遇设计暴雨过程Fig.3 Design rainstorm process of Xinhui station once in fifty years

2.3 洪水计算情境拟定

江新联围防洪保护区处于洪水、潮水综合影响区域，洪水来源为西江洪水、潭江洪水、外海风暴潮、暴雨内涝。保护区堤防属2级堤防，防洪标准为50年一遇，经过堤防整体防洪能力复核后确定防洪标准为20年一遇。通过西江与潭江段河道洪水与外海风暴潮遭遇相关性分析，考虑可能遭遇的最大风险，并根据现场调查和历史资料分析，保护区遭遇外海风暴潮时常伴有降雨发生，因此考虑洪水、风暴潮与降雨的组合；同时由于以往研究大多针对单一洪水频率工况，无法体现多重现期洪水风险程度的联系与差异性，因此本文考虑并制定多重现期洪水组合情境方案，具体情况如表2所示。

表2 广东江新联围防洪保护区洪水区划计算方案Table 2 Flood zoning calculation scheme of Jiangxinlianwei flood protection area in Guangdong Province

溃口设置方面考虑现状设防标准高一等级及以上频率洪水的堤防漫溃，对于漫溢段堤防经论证不易形成漫溢溃口的，只设漫溢段，不设溃口。漫溢溃口位置及溃决时间点根据洪水水位变化与堤防高程的相互关系进行试算确定。

3 结果与讨论

3.1 洪水风险因素分布

为研究江新联围防洪保护区洪水风险因素分布特征，需根据洪水分析情景计算方案建立多源洪水一、二维水动力耦合模型。模型构建步骤如下：(1)采用1:5万河道一维矢量数据和珠江三角洲2008年实测断面资料建立一维河网，建模范围为西江高要站以下、北江石角站以下、潭江石咀站以下及珠江三角洲部分，共包括50条河流，涉及河道总长度为1 016.3 km；(2)二维模型计算范围在保护区原有基础上向西适当延展，面积达437.6 km2，为充分反映保护区内堤防、公路、铁路、河道等线状地物的位置等特征，考虑线状地物的阻水作用，应用内边界对线状地物进行概化处理，并根据保护区内土地利用情况，对保护区内的村庄、道路、耕地、河流等地物设置了糙率分区，以反映保护区不同下垫面对洪水演进的影响，较为真实地还原洪水演进过程和水流交换效果，线状地物和糙率分区分别如图4(b)(c)所示；(3)基于河道一维与保护区二维水动力模型，叠加河道上下游洪潮边界及区域内降雨过程，通过溃口与漫溢段实现河道与保护区的耦合连接及水流信息的实时交互，以此分析洪、涝、潮多源洪水在防洪保护区的淹没情况。一、二维耦合模型如图4(a)所示。

图4 一、二维耦合模型构建Fig.4 Construction of one-dimensional and two-dimensional coupling models

为了保证模型计算精度要求及结果合理性，采用三水、马口、北街、三江口等24个水文站点的实测潮位数据对模型进行率定和验证，其中与二维保护区临近的北街站和三江口站的率定结果与验证结果如图5所示。由图可知，各水文站洪水过程计算结果与实测数据吻合；峰值、峰现时间接近，率定和验证效果较好，模型能够较为准确地模拟多源洪水条件下的水流变化特征，满足计算精度要求。

图5 各水文站率定验证结果Fig.5 Calibration and verification results of various hydrological stations

利用率定与验证后的模型参数建立一维河道与二维保护区的多源洪水耦合模型，对多重现期洪水情境方案依次进行模拟计算，以“西江、潭江100年一遇洪水、50年一遇内涝、50年一遇设计潮位”计算方案为例，提取研究区内各计算单元的风险要素指标进行洪水风险要素分析，包括最大淹没水深、最大行进流速、最大淹没历时，如图6所示。

由图4和图6对比可知，最大淹没水深较多出现在地形较低洼的地区以及溃口附近，受线状地物的阻水作用，道路和堤防两侧的淹没水深略有差异，最大洪水流速较多出现在地形坡度较大的地区，针对区域暴雨内涝情况，保护区内除局部线状地物淹没历时较小外，整个区域淹没历时均较大，参考以往研究中的洪水风险空间分布特征合理性分析与研究区域有关洪水风险分析成果[24-28]，综合分析可知，本次基于多源洪水耦合模型获得的洪水风险要素分布较为合理，结果能较为准确地反映保护区的淹没特征以及不同重现期洪水风险程度的差异性，模型的精度和可靠性较高，可为后续洪水风险区划提供详实可靠的数据支撑。

图6 洪水风险要素指标Fig.6 Flood risk factor index

3.2 区划结果分析

以“综合风险度(R)”为指标，以“最大淹没水深”为主要影响因子，综合考虑“最大行进流速”、“最大淹没历时”的影响，叠加计算多重现期洪水的3种计算方案的洪水风险要素特征，根据计算结果选取具有代表性的风险区划计算单元，绘制基于凸风险度量和期望理论的风险特征曲线，如图7所示；并根据综合风险度法区划等级划定范围对江新联围防洪保护区进行洪水风险等级区划，区划结果如图8所示。

图7 基于凸风险度量和期望理论的风险特征曲线Fig.7 Risk characteristic curves based on convex risk measure and expectation theory

为了检验综合风险度方法的可靠性，本文采用澳大利亚和尼泊尔洪泛区管理委员会[5]提出的洪水风险区划方法进行验证，该方法考虑洪水动量因素引起的洪水致灾风险，运用水深流速之积来划分洪水危害级别，区划标准参考《Designing Safer Subdivisions-Guidance On Subdivision Design In Flood Prone Areas》中风险等级划分方法，如表3所示，区划结果如图9所示。

表3 最大水深×最大流速洪水风险等级划定范围Table 3 The flood risk level range of the maximum water depth multiplied by the maximum flood velocity method

图8 综合风险度法风险区划结果 图9 最大水深×最大流速法风险区划结果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可知，根据区划计算单元绘制的风险特征曲线满足凸风险度量公理，对于不同的计算单元，随着洪水频率的增加，各洪水风险指标值呈现连续性、同质化递减趋势，这表明融合淹没水深、洪水流速、淹没历时的综合指标值能够较为准确地反映洪涝潮组合下的多源洪水风险，同时也验证了流速和历时指标区间范围划定的合理性。根据期望理论可以得知，图中阴影面积对应区划单元的综合风险度值，代表计算单元的洪水风险，洪水风险指标值增大时，基于凸风险度量和期望理论的特征曲线上移，相应的洪水风险增大。

分析图8和图9可知，洪水高风险和极高风险主要分布在棠下镇、北街街道、滘头街道、外海街道、礼乐街道以及溃口附近，统计综合风险度法和最大水深×最大流速法的各风险区域面积如表4所示，综合风险度法低风险区、中风险区、高风险区、极高风险区分别占比44.34%、28.30%、18.39%、8.97%，最大水深×最大流速法(以100洪+50潮+50涝为例)低风险区、中风险区、高风险区、较高风险区、极高风险区分别占比34.34%、21.55%、23.17%、14.71%、6.23%，最大水深×最大流速法的区划结果相较于综合风险度法风险值整体偏高，且随着洪水频率的降低风险逐渐增加。

表4 不同区划方法各风险等级区域面积统计表Table 4 Statistical table of regional area for different risk levels of different zoning methods 单位：m2

江新联围中部地区地势平坦，上游暴雨洪水在此泛滥成灾，且受下游潮位顶托作用，积水难以外排，对比分析图6中洪水风险要素指标，这些地区对应的淹没水深和淹没历时也较高，综合风险度法区划结果的合理性得以体现；江新联围北部棠下镇和南部三江镇、睦洲镇地形坡度较大，洪水流速较高，因此最大水深×最大流速法的区划结果风险较高，然而最大水深×最大流速法未考虑淹没历时的影响，强化了流速指标的作用，导致会城街道和礼乐街道的风险值偏高。综合分析淹没水深、洪水流速、淹没历时和地形高程分布情况，综合风险法洪水风险区划结果较为合理，综合风险度区划方法具有一定的合理性。

4 结论

本文借鉴金融风险管理中凸风险度量和期望理论提出基于综合风险度法的多源洪水耦合风险区划方法，在多源洪水和多重现期情境模拟分析的基础上综合考虑洪水频率、最大淹没水深、最大行进流速、最大淹没历时等风险指标，开展江新联围防洪保护区的洪水风险区划研究。研究结果表明，基于洪水风险分析模型获得的洪水风险要素分布较为合理，计算结果能较为准确地反映保护区内的淹没特征以及不同重现期洪水风险程度的联系和差异性，模型的精度和可靠性较高；综合风险度方法弥补了国内外针对单一频率洪水风险区划的不足与局限性，考虑的洪水风险因子能够较为准确地反映洪水淹没特性，区划结果反映了多个量级洪水综合淹没情况下洪水风险的空间分布特征，与最大水深乘最大流速法相比区划结果更具合理性。总之，本文提出的综合风险度区划方法更能全面直观地反映多重现期的洪水风险空间分布特征，适应于洪、涝、潮多源洪水条件下的入海河口地区的洪水风险区划，研究结果可为洪水风险社会化管理、防灾减灾规划和国土空间规划等提供决策依据和技术支撑。