基于贝叶斯时变模型和期望年损失的沿海城市洪水风险分析方法

摘要

本发明公开了一种基于贝叶斯时变模型和期望年损失的沿海城市洪水风险分析方法，包括以下步骤：(1)融合时变参数模型和贝叶斯模型平均方法得到降雨和潮位的非平稳分布；(2)运用Copula函数分析非平稳状态下联合和同现概率的变化；(3)通过城市洪水淹没模型模拟不同工况生成不同降雨和潮位组合事件的洪水损失网格；(4)基于洪水发生概率以及洪水产生的损失，采用期望年损失方法，量化变化环境下的洪水风险变化以及不确定性。本发明的基于贝叶斯时变模型和期望年损失的沿海城市洪水风险分析方法，可以合理地生成降雨和潮位的非平稳分布，并综合考虑洪水风险可能性和损失，分析变化环境下洪水风险变化及不确定性。

说明书

技术领域

本发明属于洪水风险分析技术领域，具体涉及一种基于贝叶斯时变模型和期望年损失的沿海城市洪水风险分析方法。

背景技术

沿海城市由于独特的地理位置，面临多种洪水的驱动因素，如强降雨、高潮位等。当不同的驱动因素共同作用时，城市洪水的危害程度显著增加，会造成严重的社会经济损失。根据第六次政府间气候变化评估报告(IPCC AR6)可知，在全球气温上升的背景下，极端降雨和高潮位呈现增加趋势，导致沿海城市洪水风险进一步加重。此外，沿海城市人口聚集和社会经济发达，一旦发生洪水灾害损失严重。因而，研究变化环境下降雨和潮位联合洪水风险的量化评估方法对沿海城市洪水管理至关重要。

时变分布模型可以很好地描述水文序列的非平稳性，该方法已被广泛应用于变化环境下水文事件分析。在构建时变分布模型时，选择随时间变化的分布参数是至关重要的一步。以往研究多是直接假定某一分布参数随时间变化，忽视了在变化参数选择方面造成的不确定性。当分布函数有多个参数时(如GEV分布中的位置、尺度和形状参数)，选择哪一个参数以及如何描述其变化趋势(如线性变化还是指数变化)成为了重要问题。目前，鲜有研究聚焦此类问题，而不合理的参数选择会增加非平稳状态下水文频率分析的不确定性。

贝叶斯模型平均(BMA)方法在多模型融合时发挥着重要作用，该方法最大化了每个模型结构的优点，进而提供更准确的模拟值，评估了多模型融合模拟的不确定性。而目前的时变分析和贝叶斯模型平均的研究主要是针对单变量的非平稳性分析，而已知沿海洪水受降雨和潮位双变量的共同影响。

Copula函数具有不要求单变量要有相同的边缘分布的优势，被广泛运用于双变量联合分布求解问题。水动力模型在近些年也被利用于沿海城市洪水风险评估中。然而，Copula函数和水动力模型在洪涝风险评估方面仍存在不足。Copula函数只能反映洪水发生的可能性，不会考虑洪水的淹没损失。水动力模型又只能直观表示洪水产生的后果，忽视洪水发生的概率。目前为止，未来洪涝风险及不确定性，特别是同时考虑洪水概率和洪水损失的情况，在沿海城市仍然很少见。因此，本申请提供一种洪涝风险度量标准，即期望年损失，量化降雨和潮位联合作用产生的洪涝风险以及评估变化环境对于洪涝风险不确定性的影响。

发明内容

为了解决上述问题，本发明提出一种基于贝叶斯时变模型和期望年损失的沿海城市洪水风险分析方法。本发明通过综合时变分布模型和贝叶斯模型平均方法，分析降雨和潮位的非平稳分布，降低了不确定性的影响；采用一种洪水风险度量标准，即期望年损失，量化降雨潮位联合作用产生的洪水风险以及评估变化环境对于洪水风险不确定性的影响。

本发明的基于贝叶斯时变模型和期望年损失的沿海城市洪水风险分析方法，包括以下步骤：

(1)融合时变参数模型和贝叶斯模型平均方法得到降雨和潮位的非平稳分布和不确定性区间；

(2)利用Copula函数分析非平稳状态下联合和同现概率的变化，以表征变化环境对降雨和潮位组合事件的影响；

(3)通过城市洪水淹没模型模拟不同工况生成不同降雨和潮位组合事件的洪水损失网格；

(4)基于洪水发生概率以及洪水产生的损失，采用期望年损失方法，量化变化环境下的洪水风险的变化以及不确定性。

本发明的基于贝叶斯时变模型和期望年损失的沿海城市洪水风险分析方法，可以合理地生成降雨和潮位的非平稳分布，并综合考虑洪水风险可能性和损失，分析变化环境下洪水风险变化及不确定性。研究方法和结果可以帮助决策者评估沿海城市未来的洪水风险，并为适应气候变化的可持续洪水管理提供支持。

所述步骤(1)中通过对设置的不同时变分布模型进行加权之后，得到最终的降雨和潮位频率分布，即非平稳边缘分布。

所述步骤(1)中贝叶斯模型平均的权重由期望最大值算法求解，采用蒙特卡洛抽样方法在每个水文值处生成模拟值，进而推导出模拟序列的不确定性区间。

所述不确定性区间为5％到95％分位数之间的范围。

所述步骤(2)中Copula函数采用两种Archimedean族的Copula函数和两种椭圆型Copula函数构建二元联合分布模型，并根据AIC、OLS和KS检验，选择最优的Copula函数。

所述两种Archimedean族的Copula函数为Gumbel Copula和Frank Copula，所述两种椭圆型Copula函数为Gaussian Copula和Student’s t-Copula。

所述步骤(3)中城市洪水淹没模型采用PCSWMM模型，其中一二维水动力模型的耦合方式为孔口连接方式。

所述步骤(4)中的期望年损失方法是在步骤(2)中非平稳状态下的洪水事件概率和步骤(3)中对应的洪水损失的基础上，量化得到相应地区的洪水风险变化。

附图说明

图1是本发明的基于贝叶斯时变模型和期望年损失的沿海城市洪水风险分析方法的流程图。

图2是本发明时变分布模型示意图。

图3是本发明实施例的降雨(左)和潮位(右)经验分布和候选边缘分布拟合比较。

图4是本发明实施例的城市洪水淹没模型中一二维模型耦合方式。

图5是本发明实施例的观测点的实测淹没水深和模拟水深比较。

图6是本发明实施例的降雨和潮位边缘分布参数随时间变化趋势。图中，黑色散点表示分布参数的时间序列，灰色虚线表示拟合的参数变化模型，灰色条带表示5％到95％的预测区间，(a)降雨位置参数(b)潮位位置参数(c)降雨尺度参数(d)潮位尺度参数(e)降雨形状参数(f)潮位形状参数。

图7是本发明实施例的降雨和潮位的非平稳分布和5％到95％预测区间。

图8是本发明实施例的降雨和潮位的超越概率分布和增长率分布。

图9是本发明实施例优选的Frank Copula函数理论分布和经验分布相关图。

图10是本发明实施例的降雨和潮位组合事件的非平稳状态联合概率变化，图中，(a)增加值(b)增加率(ⅰ)代表中值(ⅱ)和(ⅲ)表示5和95％的预测边界。

图11是本发明实施例的降雨和潮位组合事件的非平稳状态同现概率变化，图中，(a)增加值(b)增加率(ⅰ)代表中值(ⅱ)和(ⅲ)表示5和95％的预测边界。

图12是本发明实施例的不同降雨和潮位组合的洪水损失，图中，(a)洪水损失网格(b)典型潮位下的洪水损失随降雨变化(c)典型降雨下的洪水损失随潮位变化。

图13是本发明实施例的模拟降雨和潮位组合事件。图中，(a)模拟降雨和潮位组合分布和洪水损失等高线，线条为洪水损失等高线，单位：百万元；(b)洪水损失和相对概率乘积结果，所有小方形求和为非平稳状态下的EAD。

图14是本发明不同设计年期望年损失的特征值及增加率。图中，(a)不同设计年的EAD箱线图。箱内黑线表示中值；(b)不同设计年EAD各分位数的增加率。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

如图1所示，一种基于贝叶斯时变模型和期望年损失的沿海城市洪水风险分析方法，包括以下步骤：

(1)融合时变参数模型和贝叶斯模型平均方法得到降雨和潮位的非平稳分布和不确定性区间；

(2)利用Copula函数分析非平稳状态下联合和同现概率的变化，以表征变化环境对降雨和潮位组合事件的影响；

(3)通过城市洪水淹没模型模拟不同工况生成不同降雨和潮位组合事件的洪水损失网格；

(4)基于洪水发生概率以及洪水产生的损失，采用期望年损失方法，量化变化环境下的洪水风险的变化以及不确定性。

时变分布模型在描述变化环境下的水文特征起到了很好的效果，分布参数随时间变化，进而反映出气候变化等因素对于变量的影响(如图2)。时变GEV分布的计算式如下：

式中，μ

不同的变化参数选择对应着不同的时变分布模型。贝叶斯模型平均(BMA)是通过对不同模型的计算值进行加权得到一个综合值的方法。贝叶斯模型平均能够在模型加权时最大程度地保证最优地模型拥有最大地权重，进而降低多模型融合时的不确定性。本申请通过对设置的不同时变模型进行加权之后，得到综合的频率分布，即非平稳边缘分布。

假设y为多模型的模拟变量，D＝[R,T]为实测资料，f＝[f

式中：p(y丨D)是非平稳分布。p(f

期望最大算法(EM)可用于推求贝叶斯平均中的权重值ω

Copula函数的优点在于构建过程中不要求变量拥有相同的边缘分布，按照二元Sklar理论，H(r,t)为具有边缘分布F

H(r,t)＝C(F

本申请采用两种Archimedean族的Copula函数(Gumbel Copula，Frank Copula)和两种椭圆型Copula函数(Gaussian Copula，Student’s t-Copula)构建二元联合分布模型。根据AIC、OLS和KS检验，选择最优的Copula函数。其二元Copula函数形式如下：

表1候选Copula函数形式

注：u和v是边缘分布；

本申请选择联合概率和同现概率表征降雨和潮位联合事件的概率。降雨R和潮位T有一个变量超过某一量级的概率称为联合概率，记作P∪(r,t)。

P∪(r,t)＝P((R>r)∪(T>t))＝1-H(r,t)#(4)

降雨R和潮位T均超过某一量级的概率称为同现概率，记作P∩(r,t)。

P∩(r,t)＝P((R>r)∩(T>t))＝1-F

PCSWMM是基于SWMM开发的水文水动力模型。在本申请中，选择孔口连接方式构建一维河道、管网和二维地表洪水耦合模拟模型。动力波方程被用于水动力计算，其计算式如下：

式中：A为过水断面面积，m

期望年损失(EAD)方法能同时利用Copula函数和水动力模型，综合反映洪水发生的概率和洪水损失，更好地量化洪水风险变化。通常的计算公式为：

EAD是期望年损失，D(f)是发生概率f的洪水事件造成的洪水损失。在本申请中是评估降雨和潮位联合作用产生的复合洪水风险，因而采用二元积分方法估计洪水损失(见公式9)。为了简化积分运算，我们选择随机模拟足够数量的降雨潮位联合事件。将每次事件的概率密度函数乘以对应的洪水损失，最后求和得到期望年损失。

式中：D(x,y)为特定降雨和潮位造成的损失；f(x,y)为特定降雨和潮位对应的概率密度函数。

下面参考图3-图14，以海甸岛为例，利用本发明的基于贝叶斯时变模型和期望年损失的沿海城市洪水风险分析方法来分析海甸岛的具体情况。

本实施例涉及的数据有三个部分。第一部分是构建二元联合分布的历史降雨和潮位数据。第二部分是构建城市洪水淹没模型的基础数据，包括数字高程数据、排水管网数据和河道数据。第三部分是历史降雨潮位过程和实测的淹没数据。上述数据的来源及主要用途如下：1974年到2012年历史日降雨和潮位数据由中国海南省海口市水务局提供。年最大降雨量和当日对应的最高潮位被选择用于建立边缘分布和联合分布。数字高程数据用于构建城市洪水淹没模型中的二维模型，该数据由中国科学院资源环境科学于数据中心提供。排水管网和河道数据用于搭建城市洪水淹没模型中的一维模型，由中国海南省海口市水务局提供。历史淹没数据是城市洪水淹没模型中参数率定的基础数据，包括2014年7月威马逊台风的降雨潮位过程和观测点的实测淹没深度。

边缘分布选择：选择四种常用的函数(即GEV、Lognorm、Norm和Gamma)用于确定降雨和潮位的边缘分布。采用极大似然估计拟合分布函数中的参数。选择Kolmogorov-Simirnov(KS)检验法进行拟合检验，采用纳什效率系数(NSE)、相对误差(RE)、离差平方和(OLS)最小原则和赤池信息准则(AIC)对分布函数进行拟合优度评价。图3比较了候选边缘分布函数和经验分布函数的拟合效果，表2给出了候选边缘分布函数的检验与评价结果。

选取KS检验显著性水平α＝0.05，与样本数n＝39对应的临界值近似为0.2178，由表2可以看出所有的候选边缘分布均通过了KS检验。对于拟合优度评价而言，GEV分布的NSE指数均大于其他三种分布函数，同时RE、OLS和AIC的计算值均小于其他三种分布函数，说明广义极值(GEV)分布是降雨和潮位的最优拟合分布。

表2候选边缘分布函数的检验与评价

PCSWMM模型构建包括一二维模型的构建和模型参数率定，最终构建的模型包括2071个一维检查井，2038个管渠，23580个二维网格，13个子汇水区。一维模型通过ArcGIS进行处理构建，将研究数据中的管网、河道和数字高程数据导入到PCSWMM模型中，进行连接处理，生成一维模型。当一维模型建立后，基于设置的边界，建立二维网格，并通过孔口的方式将一二维模型进行耦合(如图4)。在建立二维模型时，采用六角形的网格，分辨率为25m，粗糙度为0.085，采样因子取3。在进行不同降雨潮位组合的模拟计算时，选择霍顿模型作为下渗模型，通过动力波的方法洪水模拟计算，其模拟演算的时间步长选取0.5s。

选取2014年7月威马逊台风的实测数据进行城市洪水淹没模型参数率定，实际的降雨和潮位过程作为模型输入的边界条件。观测点实测以及模拟深度如图5所示。由于数据的精度，如地形等，导致模拟深度和实测深度完全一致。但是，总体上实测深度和模拟深度的误差值较小，相对误差均小于15％，且NSE值达到了0.73，则认为该模型能够合理地进行城市洪水模拟。

图6显示了通过累加降雨和潮位数据的方法得到分布参数变化趋势，再运用Levenberg-Marquardt算法拟合生成不同分布参数的时变模型。从图6d和图6e看出，降雨分布的形状参数和潮位分布的尺度参数变化趋势均不明显，且相关系数小于0.4，因而这两个参数不纳入考虑时变分布模型的构建。对于其他四个分布参数，参数值与拟合模型的相关系数都超过了0.7，表明了这些拟合模型是合理的。线性和指数拟合模型如图6所示。降雨分布的位置和尺度参数呈现上升趋势，潮位的位置参数呈现上升趋势，而形状参数呈现出下降趋势。同时5％到95％的预测区间能完全覆盖参数的时间序列，反映了参数变化的所有可能，因而能够合理地表征时变分布模型地不确定性。

基于选择的四种分布参数对降雨和潮位设置不同的时变分布模型。对于降雨，模型一：位置参数时变，尺度、形状参数不变；模型二：位置、尺度参数时变，形状参数不变。对于潮位，模型一：位置参数时变，尺度和形状参数不变；模型二：位置、形状参数时变，尺度参数不变。

表3给出了各个时变分布模型在贝叶斯模型平均之后的权重值。可以明显地发现不论是降雨还是潮位，模型二的权重都高于模型一，这表明了模型二对于最后的非平稳分布有着更高的贡献度，综合考虑多个时变参数的分布模型具有更好的精度。贝叶斯模型平均的方法确保了更好的时变分布模型有更高的权重，令加权得到的非平稳分布模型更加的合理。

表3BMA模拟中不同模型的权重

为验证基于时变模型的贝叶斯模型平均方法的有效性，我们采用1974年到2003年的历史降雨和潮位数据建立非平稳模型，将预测结果推导至2012年。图7显示了降雨和潮位的非平稳分布以及5％和95％的预测区间。整体来看，降雨和潮位的非平稳分布与经验分布拟合较优，这表明了该方法是合理的。同时，预测区间能够描述贝叶斯模型平均的不确定性，5％和95％的预测区间将所有的分布模型包含，可以反映贝叶斯模型平均的所有结果。由图7也可以看出，越大的降雨和潮位有着更大的不确定性。

图8给出了基于时变分布模型的贝叶斯模型平均方法的2030年平稳和非平稳状态下的降雨和潮位的超越概率以及增加率变化。整体上，相较于平稳状态，降雨和潮位的超越概率在非平稳态显著提升。其中，降雨概率增长率随着降雨的增加而增长，平均增长率为61.09％。而潮位的概率增长率弱于降雨，平均增长率为15.08％。最大的增长率为31.11％，对应于2.82m潮位。因此表明变化环境将对强降雨和高潮位产生显著影响，进而提升沿海城市洪水发生的概率。

四种Copula函数被采纳构建降雨和潮位的二元联合分布。如表4所见，全部的候选Copula函数都通过了显著性水平0.05的KS检验。Frank Copula的AIC和OLS指数均是最小值，因此Frank Copula被选择认为是最优的函数。另外，结合图9，Frank Copula与分布点的相关系数达到了0.979，表明了该函数是表征非平稳状态下降雨和潮位相关性的最优函数。

表4候选Copula函数拟合的检验与评价

图10a和图10b显示了非平稳状态下的降雨和潮位组合下的联合概率的增加值和增加率。同时考虑了非平稳分布的不确定性，在图10的子图ii和iii给出了概率增加值和增加率的5％和95％的边界值。总体来看，由于强降雨和高潮位自身的低概率发生情况，使得概率的增加值随着降雨潮位的增加而减小，但增加率反而逐渐增大(如图10b所示)。概率的平均增加率为33.22％，最大增加率为74.72％。图10b表征了降雨潮位非平稳分布带来的不确定性，5％和95％对应的平均增加率分别为18.97％和48.59％。

同现概率的增加值随着降雨和潮位的增加而逐渐降低(如图11a所示)，这和联合概率的变化特征是一致的。概率的增加率在强降雨和高潮位时的增加情况更加明显，平均增加率为64.82％，最大增加率为153.05％(如图11b所示)。另外，由于降雨潮位边缘分布的不确定性，在5％和95％的增加率分别为27.17％和114.53％(见图11b子图ii和iii)。因而，我们注意到在非平稳状态时，强降雨和高潮位的发生概率都大幅增加，意味着变化环境对于极端降雨潮位事件有着显著影响。

图12a显示不同降雨潮位组合的洪水损失，由洪水损失成本和淹没深度计算得到。基于我们过去的研究，该研究区域的单位洪水损失成本由实地调查2010年10月和2014年7月两次洪水事件设置为258元。淹没深度通过城市洪水淹没模型模拟获得。由图12a可以可知，洪水损失随着降雨潮位的增加而增大。图12b显示了典型潮位下的洪水损失随降雨的变化。洪水损失随降雨的变化幅度很微弱，不同的降雨情境下的洪水损失差距很小。图12c显示了典型降雨下的洪水损失随潮位的变化。不同于图12b，洪水损失随着潮位的增加显著增加。例如，当降雨为400mm，潮位从3.0m增加到4.0m后，洪水损失从462.89百万元增长到1479.7百万元，增长了223.69％。因而，潮位在该研究区域的洪水损失方面发挥了主导作用。同时，我们发现潮位一旦超过3.0m后，洪水损失就会骤然增加，需要对预测超过3.0m的潮位设置更多的防洪措施。

为了同时反映洪水概率以及淹没损失，期望年损失方法用于评估洪水风险。随机模拟10000组降雨和潮位组合事件，图13a显示了模拟的降雨和潮位组合事件分布和洪水损失等高线。这些洪水损失等高线均来自于洪水损失网格(见图12a)。

图13b表示了模拟的降雨和潮位组合事件洪水损失与其相对概率的乘积结果分布。其中不同颜色的小方形代表不同降雨和潮位组合事件对期望年损失的贡献程度。在图13b中，小方形代表的值越大表明该组合事件对期望年损失的贡献度越大。我们发现贡献度较大的组合事件主要集中在低降雨和低潮位，涉及到重现期4到10年的降雨和重现期1到5年的潮位。因而，尽管一次强降雨和高潮位带来的洪水损失远高于一次低降雨和低潮位的组合事件，但由于低降雨和低潮位的组合事件发生的次数频繁，它们对洪水的期望年损失的贡献要高于强降雨和高潮位事件。

图14a给出了不同设计年期望年损失的特征值。灰箱中的黑线表示期望年损失的中值，箱子的边缘以及须边表示特定的百分位数，箱子边缘对应25和75分位数，须边对应5和95分位数。相较于平稳状态，在变化环境下期望年损失值都显著增加。对于损失中值，增加率从2030年的20.56％上升到2060年的69.84％(如图14b所示)，表明洪水风险在变化环境下将会进一步加重。另外，期望年损失的百分位数跨度随着设计年份的增加而逐步增大。如2030年期望年损失的五分位数与九十五分位数之间的差值为26.88百万元，而2060年这一数值达到了91.6百万元。这就表明了未来变化环境下会使洪水风险具有更大的不确定性。

本发明通过综合时变分布模型和贝叶斯模型平均方法，分析降雨和潮位的非平稳分布，降低了不确定性的影响；采用期望年损失，量化降雨和潮位联合作用产生的洪水风险以及评估变化环境对于洪水风险不确定性的影响。