一种考虑历史极端洪水特征的多站洪水过程随机模拟方法与流程

1.本发明涉及流域防洪减灾技术领域，具体涉及一种考虑历史极端洪水特征的多站洪水过程随机模拟方法。


背景技术：

2.在气候变化与人类活动综合影响下，全球下垫面条件与流域水文循环要素时空演变规律发生了深刻变化，流域极端气候和水文事件特性出现变异，水文节律非稳态性加剧，不确定性增强。21世纪以来，我国洪涝灾害进入了多发、群发时期，1991至2017年间我国洪水造成的直接经济损失以平均每年64％的速度增长，洪涝灾害给社会经济可持续发展带来了巨大威胁，尤其是，当遭遇类似1870年超过流域防洪体系防御标准的洪水时，会给流域带来毁灭性的影响。
3.现阶段，我国在流域超标准洪水应对方面关注较少，且受限于实测资料序列长度有限等因素，现有技术手段难以准确模拟极端气候条件下流域超标准洪水的水文特征，进一步加大了流域超标准洪水防御预案编制工作的难度。在现有数据资料背景下，大多采用“历史实测序列法”或者“典型年法”，选取某一个大洪水典型或者通过“同倍比”、“同频率”法放大得到设计洪水过程，以此来对极端洪水条件下水工程运用方案进行规划设计。然而，随着流域水工程防洪体系建成，各区防洪标准各异，区域整体防洪安全难以协调，采用传统某一频率或重现期洪水或典型大洪水进行防洪规划设计难以统筹考虑局部与整体、上下游、左右岸、干支流各区防洪安全。因此，随机模拟法作为一种合理有效的数据增扩的方法，被广泛用于解决洪水序列样本多样性不足的问题。该方法能根据历史实测资料的统计特性，无需处理复杂的洪水组合遭遇问题，模拟出大量可能的来水情况可以满足水库群常规运行、超保准运用等多种应用需求。
4.然而，现有洪水随机模拟模型受模型复杂度的限制，大多仅能保证序列低阶的统计特征，难以保证模拟的多站点水文时间序列的高阶自相关、互相关的特征。同时，流域内各区(各站点)防洪标准不一，现有多站随机模拟方法大多仅根据各站历史连续性实测洪水资料进行建模，建模过程较少考虑防洪区内各站点出现的历史超标准水文极值。因此，模拟出的极端条件下洪水序列的洪峰、洪量等模拟特征值与实际设计洪水特征值偏差较大，且模拟序列中各站点洪水特征值间的相关关系难以保持，模拟出的流域超标准洪水序列可靠性较低，这将会给超标准洪水防御方案编制工作带来数据偏差。
5.因此，如何利用有限的、不连续的水文数据资料，综合考虑流域多站点高阶的自互相关特性以及流域历史极端水文特征情况对现有随机模拟技术进行提档升级，进而模拟出合理可靠的极端气候情景下流域超标准洪水过程，是一项亟待解决的科学技术难题。


技术实现要素：

6.本发明的目的就是针对现有技术的缺陷，提供一种考虑历史极端洪水特征的多站洪水过程随机模拟方法，其能够解决采用现有随机模拟技术在对多站点洪水模拟过程中由
于模型复杂度等原因的限制，大多仅仅考虑各站点洪水过程低阶的相关特性，难以刻画流域各站点超标准洪水过程间高阶自相关、互相关性结构等问题；并解决采用现有随机模拟技术解决现有随机模拟方法难以刻画流域各站点洪水过程间高阶自相关、互相关性结构等问题，并且由于无法考虑站点历史极端洪水特征，导致模拟出的极端洪水序列可靠性较低，无法精准刻画洪水在防洪保护区内超标准的不确定性的问题。
7.本发明一种考虑历史极端洪水特征的多站洪水过程随机模拟方法，其技术方案为：
8.构建考虑高阶自互相关的多站洪水模拟模型，生成多站洪水模拟样本库；
9.结合历史极端洪水特征，生成考虑历史极端洪水的洪水特征序列集；
10.根据主站洪水特征序列，在所述多站洪水模拟样本库中匹配搜索，得到考虑历史洪水特征的多站模拟洪水序列。
11.较为优选的，所述构建考虑高阶自互相关的多站洪水模拟模型，生成多站洪水模拟样本库包括：
12.将多站点实测洪水序列协整成多站实测洪水序列矩阵，并分析多站实测洪水序列矩阵不同列变量的分布特征以及列变量之间的相关性；
13.基于多站实测洪水序列矩阵列变量间的相关性构建多变量的联合分布；
14.对所述联合分布进行随机采样，生成多站模拟洪水样本库。
15.较为优选的，结合历史极端洪水特征，生成考虑历史极端洪水的洪水特征序列集包括：
16.根据流域防洪特点，选取水库群防洪体系中的主要控制站点，并结合历史洪水资料完整程度，选取需要模拟的主站洪水特征；
17.对不连续的主站洪水特征的边缘分布以及联合分布进行拟合；
18.对所述主站洪水特征的联合分布进行相关性抽样，结合主站洪水特征边缘分布的反函数，模拟生成主站洪水特征序列集。
19.较为优选的，所述根据主站洪水特征序列，在所述多站洪水模拟样本库中匹配搜索，得到考虑历史洪水特征的多站模拟洪水序列包括：
20.以所述多站模拟洪水样本库为基础，计算所述多站模拟样本库中需要模拟的主站洪水特征的特征值，生成多站模拟洪水样本库的洪水特征序列；
21.根据模拟生成的主站洪水特征序列集，按照设定的抽样目标，在所述多站模拟洪水样本库中进行匹配搜索，得到考虑历史洪水特征的多站模拟洪水序列。
22.较为优选的，所述设定的抽样目标为：
23.洪水特征抽样误差最小、模拟洪水序列年与年之间的连续性特征与实测序列一致。
24.较为优选的，所述生成多站模拟洪水样本库后，还包括检验采样生成的多站模拟洪水样本库的合理性，所述检验过程包括：
25.根据多站模拟洪水序列计算多站模拟洪水序列矩阵的相关系数矩阵ρ
sim
；
26.基于ρ
sim
和实测洪水序列相关系数矩阵ρ
obs
计算矩阵误差指标eor；
27.若矩阵误差指标eor小于预设值，则采样生成的多站模拟洪水样本库合理。
28.较为优选的，所述根据流域防洪特点，选取水库群防洪体系中的主要控制站点，并
结合历史洪水资料完整程度，选取需要模拟的主站洪水特征包括：
29.根据研究区域各站点地理位置，选取各子流域下游站点作为主站，根据主站的洪水特性，选取需要模拟的主站洪水特征；
30.所述需要模拟的主站洪水特征包括年洪峰流量、年均流量、年最大一日洪量、年最大三日洪量、年最大七日洪量。
31.较为优选的，对不连续的主站洪水特征的边缘分布以及联合分布进行拟合包括：
32.根据选取的主站实测洪水资料，考虑主站历史大洪水的洪水特征值，对主站各洪水特征分布进行拟合；
33.计算不同主站洪水特征之间的相关系数，建立多个洪水特征之间的联合分布。
34.较为优选的，以所述多站模拟洪水样本库为基础，计算所述多站模拟样本库中需要模拟的主站洪水特征的特征值，生成多站模拟洪水样本库的洪水特征序列包括：
35.统计所述多站洪水模拟样本库中每行的洪水特征；
36.根据需要模拟的主站洪水特征进行筛选，生成多站模拟洪水样本库的洪水特征序列。
37.较为优选的，根据模拟生成的主站洪水特征序列集，按照设定的抽样目标，在所述多站模拟洪水样本库中进行匹配搜索，得到考虑历史洪水特征的多站模拟洪水序列：
38.根据模拟生成的主站洪水特征序列集cm，按照设定的抽样目标进行n1次不放回采样，在多站模拟洪水样本库的洪水特征序列cs中找到与之匹配的n1年洪水特征，并记录抽取的洪水特征在矩阵cs中所处的矩阵行号；
39.根据n1个矩阵行号在多站实测洪水序列矩阵s
obs
中取出对应行的各站洪水序列，得到n1年多站模拟洪水序列；
40.其中，在按照设定的抽样目标进行n1次不放回采样时，第一次抽样仅需满足抽样目标中的部分条件，剩余的n
1-1次抽样需满足抽样目标中的全部条件。
41.本发明的有益效果为：
42.1、解决了采用现有随机模拟技术在对多站点洪水模拟过程中由于模型复杂度等原因的限制，大多仅仅考虑各站点洪水过程低阶的相关特性，难以刻画流域各站点超标准洪水过程间高阶自相关、互相关性结构等问题；解决了采用现有随机模拟技术解决现有随机模拟方法难以刻画流域各站点洪水过程间高阶自相关、互相关性结构等问题，并且由于无法考虑站点历史极端洪水特征，导致模拟出的极端洪水序列可靠性较低，无法精准刻画洪水在防洪保护区内超标准的不确定性的问题。
43.2、本方法能模拟年、月、日不同时间尺度的多站洪水过程，模拟的洪水过程在满足各截口均值、偏态系数、变差系数等低阶统计特征的同时，还能刻画各站点洪水过程高阶自、互相关特征。
44.3、本方法克服了考虑历史洪水的不连续洪水序列随机模拟建模的难题，模拟序列能有效保持下游关键控制站点不同洪水特征的边缘分布，使模拟的极端洪水情况更加可靠。
45.4、最终抽样模拟序列能在保证各站点洪水年内相关结构的基础上，进一步保持洪水年际间的相关关系下游主要控制站点的模拟洪水的不同洪水特征。研究成果可根据实际工程需求生成考虑极端水文特征的多站洪水过程，可进一步为水工程群防洪调度方案规划
设计、风险评估提供数据支持。
附图说明
46.图1为本发明流程示意图；
47.图2～5为本发明宜昌站不同洪水特征的分布拟合结果示意图。
具体实施方式
48.为了使本技术所要解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本技术进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本技术，并不用于限定本技术。
49.需要说明的是，当元件被称为“固定于”或“设置于”另一个元件，它可以直接在另一个元件上或者间接在该另一个元件上。当一个元件被称为是“连接于”另一个元件，它可以是直接连接到另一个元件或间接连接至该另一个元件上。
50.需要理解的是，术语“长度”、“宽度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本技术和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本技术的限制。
51.此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。在本技术的描述中，“多个”的含义是两个或两个以上，除非另有明确具体的限定。
52.如图1所示，一种考虑历史极端洪水特征的多站洪水过程随机模拟方法实施步骤如下：
53.步骤1：构建考虑高阶自互相关的多站洪水模拟模型，生成多站洪水模拟样本库。该步骤主要通过计算各站点相关性，分析各站点不同截口流量过程之间的相关性，对构造各站点洪水流量联合分布，并进行联合采样，获取的模拟洪水样本库不仅能满足各个洪水变量的边缘分布，而且能有效刻画历史多站实测洪水序列的高阶自互相关结构。
54.步骤2：结合历史极端洪水特征，生成考虑历史极端洪水的洪水特征序列集。该步骤通过选取主要控制站点以及关键洪水特征，结合历史大洪水特征资料对选取的洪水特征进行分布拟合，并采取步骤1所述的相关性抽样方法对洪水特征进行联合抽样，生成洪水特征序列集。
55.步骤3：根据主站洪水特征序列，在所述多站洪水模拟样本库中匹配搜索，得到考虑历史洪水特征的多站模拟洪水序列。该步骤用于根据洪水特征序列cm，在步骤1生成的洪水过程样本库s
sim
中匹配搜索，得到考虑历史洪水特征的多站超标准洪水序列。
56.步骤1包括步骤101～103，步骤101～103具体如下：
57.步骤101，将多站点实测洪水序列协整成多站实测洪水序列矩阵，并分析多站实测洪水序列矩阵不同列变量的分布特征以及列变量之间的相关性；
58.步骤102，基于多站实测洪水序列矩阵列变量间的相关性构建多变量的联合分布；
59.步骤103，对所述联合分布进行随机采样，生成多站模拟洪水样本库。
60.其中，步骤101实施过程如下：
61.①
：将k个站点的实测洪水序列协整成洪水序列矩阵，如果模拟时间尺度为天，则实测洪水序列矩阵s
obs
如下所示：
[0062][0063]
式中n为模拟变量总数，若模拟站点总数为k，则n＝k
×
365；每一列代表模拟时间尺度。
[0064]
②
：选取合适的分布函数对每列洪水流量进行拟合，并计算矩阵s
obs
的列与列之间的皮尔逊相关系数ρ
obs
。
[0065]
其中，步骤102实施过程如下：
[0066]
③
：选取合适的copula函数，基于多站洪水序列变量间的相关性构建联合分布，联合分布函数参数可用皮尔逊相关系数ρ
obs
进行估计：
[0067]
f(q1,q2,...,qn)＝c(f1(x),f2(x),...,fn(x))
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0068]
式中，qn为第n个截口的洪水向量；c(x)为不同截口洪水的联合分布函数；f(x)为n个洪水变量的联合分布函数，fn(x)为第n个洪水变量的边缘分布。
[0069]
其中，步骤103实施过程如下：
[0070]
④
：随机生成服从均匀分布的向量(ε1,ε2,...,εn)，通过边缘分布的反函数可求得站点t时刻的流量值；然后令u1＝ε1，由于u1和ε2已知，通过条件分布c(u2|u1)＝ε2可计算得到u2，以此类推，通过条件分布c(ui|u1,...,u
i-1
)＝εi可得唯一值ui，最后得到模拟向量(u1,u2,
…
,un)；重复m次即可得到m年多站洪水过程。
[0071]
⑤
：通过每个抽样变量的边缘分布的反函数i＝1,...,n，可求得模拟洪水序列矩阵s
sim
。
[0072][0073]
式中，s
sim
为模拟序列矩阵；其结构与s
obs
类似，不同的是行数为模拟总年数m。
[0074]
⑥
：计算模拟序列矩阵的相关系数矩阵ρ
sim
，结合实测洪水序列相关系数矩阵ρ
obs
通过以下的矩阵误差指标eor来判断模拟序列的相关结构保持程度。
[0075][0076]
式中，eor代表实测序列与模拟序列相关系数矩阵的rmse，其值越小，代表模拟序
列的相关结构保持的越好。
[0077]
步骤2包括步骤201～203，步骤201～203具体如下：
[0078]
步骤201，根据流域防洪特点，选取水库群防洪体系中的主要控制站点，并结合历史洪水资料完整程度，选取需要模拟的主站洪水特征；
[0079]
步骤202，对不连续的主站洪水特征的边缘分布以及联合分布进行拟合；
[0080]
步骤203，对所述主站洪水特征的联合分布进行相关性抽样，结合主站洪水特征边缘分布的反函数，模拟生成主站洪水特征序列集，检验采样生成的多站模拟洪水样本库的合理性。
[0081]
其中，步骤201实施过程如下：
[0082]
①
：根据研究区域各站点地理位置，选取各子流域下游站点作为主站，根据主站的洪水特性，选取需要模拟的主站洪水特征。主站洪水特征具体可包括：年洪峰流量、年均流量、年最大一日洪量、年最大三日洪量、年最大七日洪量等。
[0083]
其中，步骤202实施过程如下：
[0084]
②
：根据选取的主站实测洪水资料，考虑主站历史大洪水的洪水特征值，对主站各洪水特征分布进行拟合。对于中国大部分地区水文站点，多采用皮尔逊三型分布(p-iii)对水文特征进行分布拟合，并结合历史大洪水资料，采用不连续矩法估计分布的参数。
[0085]
③
：计算不同主站洪水特征之间的皮尔逊相关系数，采用copula函数建立n个洪水特征之间的联合分布：
[0086]
f'(ch1,ch2,...,chn)＝c'(f1(x),f'2(x),...,f'n(x))
[0087]
式中，chn为第n个主站实测洪水特征向量，由所选主站实测连续性洪水资料统计得到；c'(x)为不同洪水特征间的联合分布函数；f'(x)为n个洪水特征的联合分布函数，f'n(x)为第n个洪水特征变量的边缘分布，其参数由步骤(2)估算得到。
[0088]
其中，步骤203实施过程如下：
[0089]
根据多站模拟洪水序列计算多站模拟洪水序列矩阵的相关系数矩阵ρ
sim
；
[0090]
基于ρ
sim
和实测洪水序列相关系数矩阵ρ
obs
计算矩阵误差指标eor；
[0091]
若矩阵误差指标eor小于预设值，则采样生成的多站模拟洪水样本库合理。
[0092]
其具体可采用以下方式进行验证。
[0093]
④
：基于采用步骤1所述的(4)～(5)对n个洪水特征的联合分布f'(x)进行抽样，其中将步骤1第(4)步中的径流变量间的条件分布c(u2|u1)替换成洪水特征间的条件分布c'(u2|u1)，将步骤1第(5)步的抽样变量的边缘分布的反函数替换成洪水特征边缘分布的反函数对主站洪水特征进行相关性抽样，生成主站洪水特征模拟序列矩阵cm。
[0094][0095]
式中，n1为模拟洪水特征序列的总年数；矩阵中每一行代表每一个模拟年中各主站考虑洪水特征，总数为l，若选取k个主站，r个洪水特征，则l＝k
×
r；cm
i,j
是第i个模拟年份中第j个模拟洪水特征，i＝1,2
…
n1，j＝1,2,
…
l。矩阵cm可作为洪水特征序列，其中各洪水特征的边缘分布计算过程中考虑了历史极端洪水特征值，与此同时，样本库内各洪水特
征间的相关性结构也与实际洪水特征间的相关性结构保持一致。
[0096]
步骤3包括步骤301～302，步骤301～302具体如下：
[0097]
步骤301，以所述多站模拟洪水样本库为基础，计算所述多站模拟样本库中需要模拟的主站洪水特征的特征值，生成多站模拟洪水样本库的洪水特征序列。
[0098]
步骤302，根据模拟生成的主站洪水特征序列集，按照设定的抽样目标，在所述多站模拟洪水样本库中进行匹配搜索，得到考虑历史洪水特征的多站模拟洪水序列。所述设定的抽样目标为：洪水特征抽样误差最小、模拟洪水序列年与年之间的连续性特征与实测序列一致。
[0099]
其中，步骤301实施过程如下：
[0100]
①
：统计步骤1生成的洪水特征序列s
sim
中每年(行)的洪水特征，并根据步骤2选取的各站点洪水特征进行筛选，生成样本库的主站洪水特征矩阵cs。
[0101][0102]
式中，fs有m行r列，csi,j代表模拟序列中主站第i年第j个洪水特征，其中i＝1,2
…
m，j＝1,2,
…
l。
[0103]
其中，步骤302实施过程如下：
[0104]
②
：基于步骤2生成的洪水特征矩阵cm，在cs中找到与之匹配的n1年洪水特征，记录抽取的洪水特征在矩阵cs中所处的矩阵行号。
[0105]
其抽样目标分为两部分，第一部分为最小化洪水特征抽样误差，采用cs与cm之间各洪水特征的相对误差值之和来量化表征。第二部分则采用相关系数差值来代表各站点模拟洪水序列年与年之间的连续性特征与实测序列一致性。
[0106][0107]
式中，cmk与csk分别为特征矩阵cm与cs对应行数的第k个洪水特征值；ns是站点总个数；na是一年内的总的截口数；和分别是模拟序列与实测序列第j个站点第p个截口和第q个截口洪水的相关系数，其中截口p与截口q的时滞为l。
[0108]
③
：根据抽样目标obj，进行n1次不放回采样，对于第一次抽样，仅需满足模拟不同洪水特征相对误差值之和最小的目标，对于剩余n1-1次不放回抽样，则需严格按照目标obj进行抽样。
[0109]
④
：根据n1个行号在步骤1模拟得到的洪水序列矩阵s
obs
中取出对应行的各站洪水序列，进而模拟出n1年多站模拟洪水序列。
[0110]
根据此目标搜索得到的洪水过程能使模拟序列洪水特征的边缘分布与考虑历史大洪水情况下的边缘分布吻合，模拟出的极端情况更加可靠，并且模拟洪水特征具有时间序列特有的连续性，且与实测序列的连续性保持一致。
[0111]
实施例一
[0112]
本实施例以金沙江下游屏山站以及长江干流寸滩站、宜昌站为研究对象，进行实例仿真，以验证本发明的效果。具体实施步骤如下：
[0113]
步骤1.首先将三站洪水序列协整成实测洪水序列矩阵，选取特定分布对各站点日洪水过程进行分布拟合，并分析计算实测洪水序列矩阵相关结构。在此基础上，采用gaussian copula函数拟合各站点不同日洪水的联合分布，并进行相关性抽样，结合每日的边缘分布模拟得到能有效保证高阶自互相关特征的n1年三站日洪水过程。
[0114]
步骤2.选取下游主要控制站点宜昌站作为主要控制站，选取年最大洪峰流量、年最大7日洪量、年最大15日洪量、年最大30日洪量四个洪水特征作为极端洪水特征量化指标。采用步骤1所述的方法对主站四个洪水特征进行模拟，得到宜昌站n2年模拟洪水特征序列集。
[0115]
步骤3.根据步骤1得到的n1年宜昌站洪水序列，统计其宜昌站每一年的年最大洪峰流量、年最大7日洪量、年最大15日洪量、年最大30日洪量值，并根据步骤2得到宜昌站n2年模拟洪水特征序列集，以最小化宜昌站各洪水特征间的误差与保证各站点抽样洪水的年际连续性为目标，在步骤1得到的n1年三站模拟洪水序列中进行采样，得到满足条件的n2年三站同期洪水过程。
[0116]
根据步骤1模拟得到10000年三站模拟日洪水过程，并将模拟结果与多站季节性自回归模型模拟结果进行对比，根据式(4)计算采用所提方法得到的洪水序列相关系数矩阵误差eor值为0.0146，而采用msar法得到的eor值为0.149，远大于所提方法计算得到的eor值。由此可得，所提方法模拟的洪水序列能精确刻画实测洪水序列的高阶自相关、互相关特性。
[0117]
同时根据宜昌站历史设计资料以及实测洪水资料，模拟出1000年洪水特征序列集，以满足主站洪水特征误差最小与洪水序列连续性特征为目标，在10000年三站模拟序列中进行采样，得到的1000年三站同期模拟洪水序列。其中宜昌站不同洪水特征边缘分布如图2～4所示，从图中可得，模拟序列中洪峰、年最大7日洪量、年最大15日洪量以及年最大30日洪量拟合得到的分布与考虑极端大洪水特征得到的洪水特征理论分布差别较小，所提方法得到的模拟洪水序列能够刻画洪水特征整体分布特征。
[0118]
表1展示了宜昌站不同设计频率下模拟值与实测值，从表中可以看出两者相差很小，证明了通过该方法模拟出的洪水特征设计值与实际设计值偏差不大，所提方法模拟的极端超标准洪水情景更加可靠。
[0119]
表1宜昌站设计洪水值对比结果
[0120][0121]
此外，统计了各站点实测序列与模拟序列每年第一日与第二日两个截口与前一年最后两日的相关系数矩阵，并统计了实测相关系数矩阵与模拟相关系数矩阵的绝对误差值，如表2所示。由表可知，年际相关系数绝对误差值最大仅为0.0289，几乎可忽略不计，总体而言，该模型生成的模拟洪水序列能较好保持实测洪水的年际相关特征。
[0122]
表2各站点洪水年际相关系数矩阵绝对误差
[0123][0124]
以上所述实施例仅用以说明本技术的技术方案，而非对其限制；尽管参照前述实施例对本技术进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本技术各实施例技术方案的精神和范围，均应包含在本技术的保护范围之内。