一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法

1.本发明涉及水利工程领域中的中长期水文预报技术，特别涉及一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法。

背景技术：

2.中长期径流预测是指预见期超过月、季等期长的水文径流预测，它是水利工程调度计划编制的基本依据，是水资源优化调度的基础信息支撑。不同于短历时的洪水预报，长历时的径流成因复杂，难以准确预测，导致中长期径流预测误差大、不确定性高。目前，中长期径流单值预测方案仍难以为水资源规划与管理提供准确可靠的信息。由于水资源系统可通过水库径流调节对径流预测误差实现缓冲和容错，在难以获取准确的单值预测结果的条件下，通过辨识单值预报方案误差的变化特征，融合单值预报和误差模拟构建概率预测方案，开展径流区间及其可能性的定量预测，可以提高预测信息对调度的支撑作用。
3.决定概率预测方案质量的指标主要包括精准度、可靠性和稳定性三类指标。精准度主要评估单值预测方案与实测值的偏差；可靠性主要评估正常条件下概率预测方案覆盖实测值的可能性；稳定性主要评测在非常(极端)条件下预测区间无法有效覆盖实测值时的可能偏差。追求精准、可靠和稳定的概率预测方案是径流预测的最终目标。影响这三项指标的预测模型因素主要包括预测因子(自变量)、预测方程(模型结构)和预测误差：预测因子和预测方程是径流形成的确定性物理过程的反映，前者代表影响径流形成的水文循环关键要素，后者通过构建数学模型描述这些要素与径流之间的响应关系；而预测误差是径流预测及水文循环中不确定性的综合反映，一般采用预测径流与实际径流的偏差来衡量。选取不同的预测因子、方程既影响单值预测的结果，同时也决定了误差的分布特性，最终将产生不同的预测效果。因此，通过遴选预测因子、改变模型结构、优化模型参数可以改变预测整体效果。
4.现有研究针对上述三项工作展开了一系列探索，在遴选因子方面，对气象、水文因子进行了分类筛选；在模型结构方面，研发了时间序列、统计回归、神经网络等一系列数据驱动类模型方法；在参数设定方面，通常假定模型参数具有平稳性特征，采用离差平方和最小的方式进行参数率定。然而，静态平稳性结构的预测模型方法越来越难以客观反映径流过程，尤其是随机项的非平稳及时变规律：例如径流过程的响应受气象、水文要素在产汇流过程中时、空延续效应和复杂性影响，具有时空关联、尖峰厚尾等特征；传统平稳时不变的预测模型方法面临重大挑战。因此，针对中长期径流预测中误差及不确定性的时变特点，如何构建精准、可靠且稳定的径流时变预测模型是亟需解决的科学问题。

技术实现要素：

5.发明目的：本发明的目的是提供一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法，通过对预测因子、模型结构、模型参数三要素分别优化，提高
径流预测的精准度、可靠性、稳定性，为水资源高效利用提供有效、丰富的决策信息。
6.技术方案：本发明的一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法，包括以下步骤：
7.s1、作为“输入-结构-参数”全要素分层组合优选的第一层，对模型输入进行优选：收集研究区长系列中长期包括大气环流因子、降雨、蒸发、径流在内的气象、水文因子数据资料构建径流预测初始因子集，基于copula熵与径流成因分析筛选预测因子组合得到驱动因子集，采用长短期记忆神经网络lstm模型生成中长期径流单值预测方案集；
8.s2、作为“输入-结构-参数”全要素分层组合优选的第二层，对模型结构进行优选：以信息熵作为评价单值预测方案集信息价值的指标，采用具有内外层循环嵌套结构的逐步回归对单值预测方案集中的方案组合进行优选，对最优子集采用贝叶斯模型平均bma的方式生成同时包含单值过程与概率区间的径流集合预测方案；
9.s3、作为“输入-结构-参数”全要素分层组合优选的第三层，对模型参数进行优选：基于广义自回归条件异方差模型garch对径流预测误差进行精细表征，利用误差的相依关系进行实时误差校正，结合蒙特卡洛抽样生成校正后的误差情景序列，叠加至步骤s2中的单值过程上得到考虑误差校正的径流时变概率预测方案，依据径流预测“精准度-可靠性-稳定性”多维评价指标体系对得到的概率预测方案进行评价。进一步的，步骤s1包括以下步骤：
10.s11、径流预测因子集生成；
11.根据研究区域气象条件，确定控制研究区域降雨、径流的主要天气要素，选取相应的指标与前期降雨、蒸散发量以及径流量共同组成径流预测因子集；
12.s12、基于copula熵与径流成因分析的驱动因子识别；
13.设n维随机变量(x1,x2,...,xn)具有边缘分布ui＝fi(xi),i＝1,2,...,n，其中f为累积概率密度函数，copula熵定义为：
[0014][0015]
其中，ce为copula熵，c为copula函数，基于copula熵对所有与径流成因有关的预测因子进行关联度排序，选取排序靠前的若干因子作为驱动因子集；
[0016]
s13、lstm模型的径流单值预测方案集生成；
[0017]
采用lstm对月入湖径流进行预测，以不同因子组合作为因变量进行学习、训练、模拟并进行预测。
[0018]
更进一步的，步骤s13中lstm的计算方法为：
[0019]
(a)遗忘门：
[0020]ft
＝σ(wf·
[h
t-1
,x
t
]+bf)
[0021]
其中，f
t
为当前时段的遗忘门门控，表示控制部分历史输出信息转换为当前时段的输入，反映历史信息在当前时段的丢失程度，wf为遗忘门的权重矩阵，bf为遗忘门的偏置项，σ为激活函数；
[0022]
(b)输入门：
[0023]it
＝σ(wi·
[h
t-1
,x
t
]+bi)
[0024]
其中，i
t
为当前时段的输入门门控，表示控制部分当前时段的输入转化为系统状
态变量，反映当前时段的信息输入效率，wi为输入门的权重矩阵，bi为输入门的偏置项；
[0025]
(c)当前输入的长期状态c
t
'：
[0026]ct
'＝tanh(wc·
[h
t-1
,x
t
]+bc)
[0027]
其中，wc为当前输入的长期状态的权重矩阵，bc为wc的偏置项，tanh为激活函数；
[0028]
(d)当前时刻的长期状态c
t
：
[0029]ct
＝f
t
·ct-1
+i
t
·ct
'
[0030]
(e)输出门：
[0031]ot
＝σ(wo·
[h
t-1
,x
t
]+bo)
[0032]
其中，o
t
为当前时段的输出门门控，表示控制部分系统状态变量转化为当前时段的输出，反映当前时段的信息输出效率，wo为输出门的权重矩阵，bo为wo的偏置项；
[0033]
(f)lstm的输出值：
[0034]ht
＝o
t
·
tanh(c
t
)
[0035]
将计算时段划分为训练期与校验期，将训练期中经步骤s12中识别的驱动因子集w作为上述计算步骤中的输入项xt，将径流量作为其中的输出项ht，经神经网络的训练后率定得到参数wf、bf、wi、bi、wo、bo，并与校验期的输入项xt一同代入模型中计算得到对应的输出项ht，得到一组单值预测过程，构成单值预测方案集。
[0036]
进一步的，步骤s2包括以下步骤：
[0037]
s21、基于信息熵与逐步回归的方案集组合优选；
[0038]
信息熵是度量信息系统不确定性的统计量，其表达式为：
[0039][0040]
其中，h为信息熵，x为一个具有n个离散取值的随机变量，x＝{x1,x2,...,xn}，其概率分布列此外，对于一个二元随机变量(x,y)，其信息熵又称为联合熵，定义为：
[0041][0042]
其中，p(x,y)为随机变量x与y的联合概率密度函数；联合熵表示多个随机变量包含的信息总量，而其信息同量被定义为互信息，即：
[0043]
t(x,y)＝h(x)+h(y)-h(x,y)
[0044]
因此，用信息熵中的联合熵与互信息分别描述单值预测方案集中各候选方案间的信息总量与信息同量，并以此作为方案组合优选的依据；
[0045]
s22、基于bma的径流单值-概率集合预测方案生成；
[0046]
设y为预测变量，y
obs
为实测样本序列，f＝{f1,f2,...,fk}为候选模型空间，其中k表示步骤s21中优选的单值预测方案组合内的候选方案组数；根据贝叶斯理论，假定实测与预报序列均服从正态分布，进一步细化：
[0047]
[0048]
其中，p(y|y
obs
)表示预测变量的概率分布，ωi表示第i个模型的权重，表示均值为fi，方差为的正态分布；
[0049]
因此，采用期望最大化算法率定bma参数ωi和并采用monte-carlo抽样给出概率预测情景序列，获得包含单值过程与概率区间的径流集合预测方案。
[0050]
更进一步的，步骤s22中采用monte-carlo抽样给出概率预测情景序列，具体步骤如下：
[0051]
(a)采用box-cox变换将非正态分布的径流序列转化成正态分布，具体步骤为：
[0052][0053]
其中，λ为转换系数，其取值由极大似然估计给出；y'为经过转化后的服从正态分布的径流预测序列；
[0054]
(b)采用e-m算法率定bma参数ωi和
[0055]
(c)根据各候选方案的bma权重[ω1,ω2,...,ωk]，定义累积权重值
[0056][0057]
生成[0,1]均匀随机数u，如果则表示此次抽样选取的是第i个候选方案fi；
[0058]
(d)根据fi在时段j的概率分布随机生成服从正态分布的径流预测序列y'；
[0059]
(e)基于累积权重值公式对生成的服从正态分布的径流预测序列y'进行还原，计算公式为：
[0060][0061]
(f)重复步骤(c)-(e)m次，得到一组具有m个情景的预测情景集，即bma的径流概率预测方案，其均值过程即为bma的径流单值预测方案；二者共同构成径流单值-概率集合预测方案。
[0062]
进一步的，步骤s3包括以下步骤：
[0063]
s31、基于garch模型的径流预测误差精细表征；
[0064]
假定残差项ε
t
与条件方差h
t
服从：
[0065]
结构-参数”全要素分层组合优选的第三层。
[0082]
本发明的一种装置设备，包括存储器和处理器，其中：
[0083]
存储器，用于存储能够在处理器上运行的计算机程序；
[0084]
处理器，用于在运行所述计算机程序时，执行如上述一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法的步骤。
[0085]
本发明的一种存储介质，所述存储介质上存储有计算机程序，所述计算机程序被至少一个处理器执行时实现如上述一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法的步骤。
[0086]
有益效果：与现有技术相比，本发明的显著技术效果为：本发明通过对水文径流预测模型的预测因子、模型结构、模型参数三要素分别开展优化与优选，通过初筛并优选预测驱动因子的组合、剔除冗余的预测信息，实现对预测模型“输入-结构-参数”的分层组合优选，提高径流预测的精准度、可靠性、稳定性，能够为水资源高效利用提供有效、丰富的决策信息。
附图说明
[0087]
图1是本发明方法的流程图；
[0088]
图2是方案集组合优选时具有内外层双循环嵌套结构的逐步回归过程流程图；
[0089]
图3是实施案例中全要素优选以及部分要素优选的月径流概率预测对比示意图。
具体实施方式
[0090]
下面结合附图和具体的实施例对本发明进行详细说明。然而，对于本领域技术人员而言显而易见的是，本发明可以无需一个或多个这些细节而得以实施。在其他的例子中，为了避免与本发明发生混淆，对于本领域公知的一些技术特征未进行描述。本发明中的相关步骤的顺序并非是限定性，即本领域的技术人员可以调整，本发明中的顺序是一种案例性的撰写方式，而非限定性描述。
[0091]
本发明的一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法，包括收集研究区长系列中长期径流、降雨、蒸发资料以及大气环流因子数据，构成径流预测因子集，基于copula熵与径流成因分析对预测因子筛选并自由组合得到驱动因子集，采用基于长短期记忆神经网络lstm模型中生成中长期径流单值预测方案集；以信息熵作为评价单值预测方案集信息价值的指标，采用具有内外层循环嵌套结构的逐步回归对方案集中的方案组合进行优选，对最优子集采用贝叶斯模型平均bma的方式生成同时包含单值过程与概率区间的径流集合预测方案；基于广义自回归条件异方差模型garch对径流集合预测误差进行精细表征，结合蒙特卡洛抽样生成误差情景序列，叠加至bma单值过程上得到考虑误差校正的径流概率预测方案。本发明针对水文径流预测模型的预测因子、模型结构、模型参数三要素，采用“输入-结构-参数”分层组合优选的方式对预测模型进行升级调控，能够有效提升预测精度、可靠性与稳定性，为水资源规划与管理提供有效、丰富的预报决策信息。
[0092]
如图1所示，本发明的一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法，包括以下步骤：
[0093]
s1、输入优选：径流单值预测方案生成。
[0094]
作为“输入-结构-参数”全要素分层组合优选的第一层，对模型输入进行优选：收集研究区长系列中长期包括大气环流因子、降雨、蒸发、径流在内的气象、水文因子数据资料构建径流预测初始因子集，基于copula熵与径流成因分析筛选预测因子组合得到驱动因子集，采用lstm模型生成中长期径流单值预测方案集。
[0095]
具体包括以下步骤：
[0096]
s11、径流预测初始因子集生成；
[0097]
根据区域产汇流特征，选取对当前时段有直接或间接影响的水文要素作为预报因子，构建物理或数学模型描述预报因子(自变量)与径流量(因变量)的复杂非线性关系，即：
[0098]qt
＝f(r
t-1
,r
t-2
,
…
,r
t-τ
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0099]
式中，q
t
为t时刻的径流量，r为预测因子，其下标表示对应时刻，f为预测方程，τ为预测因子的预见期。
[0100]
预测因子的选取直接影响预测精度。一般地，影响中长期径流成因的水文要素包括前期径流、降雨和蒸散发，气象要素包括气温、气压以及环流指数，其中气象要素提供定性的趋势预测，而水文要素提供定量的数值预测。根据研究区域气象条件，确定控制研究区域降雨、径流的主要天气要素(如季风、海平面气压、副热带高压等)，选取相应的指标(如季风指数、南方涛动指数、副高面积指数等)，与前期降雨、蒸散发量以及径流量等共同组成径流预测初始因子集。
[0101]
s12、基于copula熵与径流成因分析的驱动因子识别；
[0102]
采用全因子作为模型输入不仅会显著增加计算耗时，而且可能导致模型泛化，因此需要对预报因子进行初筛，其关键在于寻与径流序列相关性强、成因解释得当的“驱动因子”。copula熵是一种具有全阶次、多变量的非线性相关性度量指标，可用于描述预测因子与径流序列的相关程度。设n维随机变量(x1,x2,...,xn)具有边缘分布ui＝fi(xi),i＝1,2,...,n，其中f为累积概率密度函数，copula熵可定义为：
[0103][0104]
式中，ce为copula熵，c为copula函数。因此，本发明基于copula熵对所有与径流成因有关的预测因子进行关联度排序，选取排序靠前的r项因子作为驱动因子集w＝{w1,w2,
…
wr}。
[0105]
s13、lstm模型的径流单值预测方案集生成；
[0106]
采用lstm对径流进行预测，以不同因子组合作为因变量进行学习、训练、模拟并进行预测。lstm是循环神经网络(rnn)的变体，能够解决一般rnn模型在长序列回归上的梯度爆炸与梯度消失问题。lstm由一个输入层、一个或多个隐含层和一个输出层组成，其中每个隐含层包含两个状态变量h和c，分别用于保存短期状态与长期状态。二者通过门控单元进行连接并构成了隐含层，门控单元由遗忘门、输入门和输出门构成，其中遗忘门决定上一时刻的输入有多少可以保留至当前时刻，而输入门决定当前时刻的输入有多少可以保存到长期状态c中，输出门决定长期状态c有多少可以作为lstm的输出值。lstm的计算方法为：
[0107]
(a)遗忘门：
[0108]ft
＝σ(wf·
[h
t-1
,x
t
]+bf)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0109]
式中，f
t
为当前时段的遗忘门门控，表示控制部分历史输出信息转换为当前时段的输入，反映历史信息在当前时段的丢失程度，wf为遗忘门的权重矩阵，bf为遗忘门的偏置项，σ为激活函数，一般用sigmoid函数；
[0110]
(b)输入门：
[0111]it
＝σ(wi·
[h
t-1
,x
t
]+bi)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0112]
式中，i
t
为当前时段的输入门门控，表示控制部分当前时段的输入转化为系统状态变量，反映当前时段的信息输入效率，wi为输入门的权重矩阵，bi为输入门的偏置项；
[0113]
(c)当前输入的长期状态c
t
'：
[0114]ct
'＝tanh(wc·
[h
t-1
,x
t
]+bc)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0115]
式中，wc为当前输入的长期状态的权重矩阵，bc为wc的偏置项，tanh为激活函数；
[0116]
(d)当前时刻的长期状态c
t
：
[0117]ct
＝f
t
·ct-1
+i
t
·ct
'
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0118]
(e)输出门：
[0119]ot
＝σ(wo·
[h
t-1
,x
t
]+bo)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0120]
式中，o
t
为当前时段的输出门门控，表示控制部分系统状态变量转化为当前时段的输出，反映当前时段的信息输出效率，wo为输出门的权重矩阵，bo为wo的偏置项。
[0121]
(f)lstm的输出值h
t
：
[0122]ht
＝o
t
·
tanh(c
t
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0123]
将计算时段划分为训练期与校验期，将训练期中经步骤s12中识别的驱动因子集w作为上述计算步骤中的输入项x
t
，将径流量作为其中的输出项h
t
，经神经网络的训练后率定得到参数wf、bf、wi、bi、wo、bo，并与校验期的输入项x
t
一同代入模型中计算得到对应的输出项h
t
，即可得到一组单值预测过程，构成单值预测方案集，如图3中(c)所示，考虑预测因子优选的预测方案能够尽可能覆盖实测真值，具有较高的可靠性，但其区间宽度过大，精准性有待提高。
[0124]
s2、结构优选：基于bma的径流集合预测方案生成；
[0125]
作为“输入-结构-参数”全要素分层组合优选的第二层，对模型结构进行优选：以信息熵作为评价单值预测方案集信息价值的指标，采用具有内外层循环嵌套结构的逐步回归对单值预测方案集中的方案组合进行优选，对最优子集采用贝叶斯模型平均bma的方式生成同时包含单值过程与概率区间的径流集合预测方案。
[0126]
具体包括以下步骤：
[0127]
s21、基于信息熵与逐步回归的方案集组合优选；
[0128]
由于水文模型的“同参异效性”以及驱动因子的一致性，单值预测方案集内各候选方案间存在一定的信息冗余，导致计算规模的扩大与计算效率的降低，过量的输入可能会给模型带来大量噪声信息，导致模型输出不确定性的放大，降低了模型的稳定性，因此须通过组合优选的方式去除冗余信息，保留关键信息。该过程是一个具有内外层循环嵌套结构的逐步回归过程，如图2所示，即：内层循环包括，从单值预测方案初始全集出发，按编号顺序从小到大依次删去某一方案，并计算该方案与剩余方案子集间的信息同量，选取信息冗余最多的方案，将其剩余方案子集更新为新的全集，并重复上述操作直至剩余方案子集内元素个数为1，循环终止；外层循环包括，对具有不同元素个数的剩余方案子集，计算其信息
总量，并将其与初始全集的信息总量相比，当该比值小于某一给定阈值时，取此时的全集为最优组合，循环终止。具体的，方案集组合优选时具有内外层双循环嵌套结构的逐步回归过程流程为：
[0129]
(1)令m＝f；m＝|f|；其中，f为单值预测方案初始全集，m表示单值预测方案更新全集；m表示更新全集的方案个数；
[0130]
(2)令i＝1，si＝m\{mi}，其中，mi表示第i个方案，si表示若删去第i个方案后的单值预测方案剩余子集；
[0131]
(4)计算obji＝cv(mi,si)，cv表示剩余子集si与更新全集mi的信息同量，此时完成一次内层循环；
[0132]
(5)对i从1取到m，即遍历更新全集中的每个方案，内层循环继续，重复上述步骤，即可得到若删去集合内任一方案后剩余子集的信息冗余情况，选取对应剩余子集冗余最高的方案确定删去，得到最优剩余子集s
min
，内层循环终止；
[0133]
(6)将最优剩余子集s
min
设置成新的更新全集m，分别计算其与初始全集f的信息总量h，此时完成一次外层循环；
[0134]
(7)当更新全集m与初始全集f的信息总量比(即相对信息总量)大于某一给定阈值β时(即h(m)/h(f)≤β不成立)，外层循环继续，继续计算具有不同元素个数的更新全集m的相对信息总量；当相对信息总量小于某一给定阈值或更新全集m元素个数减少为2时(即m＝2)，外层循环终止，此时的更新全集m即为确定的最优子集sf。
[0135]
本发明选取信息熵作为衡量预测方案信息价值的指标。信息熵是度量信息系统不确定性的统计量，它从随机变量的概率分布出发，能够刻画不同分布间的独立性、关联性以及冗余性，因此被广泛应用于预测、优化、决策等多个场景。信息熵的表达式为：
[0136][0137]
式中，h为信息熵，x为一个具有n个离散取值的随机变量，x＝{x1,x2,...,xn}，其概率分布列此外，对于一个二元随机变量(x,y)，其信息熵又称为联合熵，定义为：
[0138][0139]
式中，p(x,y)为随机变量x与y的联合概率密度函数。联合熵表示多个随机变量包含的信息总量，而其信息同量可被定义为互信息，即：
[0140]
t(x,y)＝h(x)+h(y)-h(x,y)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(11)
[0141]
因此，可用信息熵中的联合熵与互信息分别描述单值预测方案集中各候选方案间的信息总量与信息同量，并以此作为方案组合优选的依据。
[0142]
s22、基于bma的径流单值-区间集合预测方案生成；
[0143]
贝叶斯模型平均(bma)在集合概率预报与水文不确定性分析中有着广泛的应用，该方法基于贝叶斯理论，分析各候选模型能够成为最优模型的概率，并在最终的模型选择中倾向于选择概率更大的模型，因此能够有效消除模型不确定性的不利影响。
[0144]
设y为预报变量，y
obs
为实测样本序列，f＝{f1,f2,...,fk}为候选模型空间，其中k
表示步骤s21中优选的单值预测方案组合内的候选方案组数。根据贝叶斯理论，给定样本y
obs
下预报变量y的概率密度函数p(y|y
obs
)为：
[0145][0146]
式中，p(y|fi,y
obs
)为在给定样本y
obs
和模型fi的条件下预报变量y的概率密度函数，即“后验分布”；p(fi|y
obs
)为在给定样本y
obs
的情况下第i个预测模型为最优模型的概率，即bma权重ω。
[0147]
显然，变量y的预报能够以确定性单值结果或概率密度函数的形式给出，因此可以开展单值过程预测与概率区间预测。假定实测与预报序列均服从正态分布，可进一步细化：
[0148][0149]
式中，p(y|y
obs
)表示预测变量的概率分布，ωi表示第i个模型的权重，表示均值为fi，方差为的正态分布。
[0150]
因此，采用期望最大化(e-m)算法率定bma参数ωi和进行率定，并采用monte-carlo抽样能够给出概率预测情景序列，具体步骤如下：
[0151]
(a)采用box-cox变换将非正态分布的径流序列转化成正态分布，具体步骤为：
[0152][0153]
式中，λ为转换系数，其取值由极大似然估计给出；y'为经过转化后的服从正态分布的径流预测序列；
[0154]
(b)采用e-m算法率定bma参数ωi和令迭代次数iter＝0，假定初始权重为均匀权重，即计算初始方差其中t为计算时段数。此时初始似然函数值为当进行到第iter次迭代时，有隐藏变量：
[0155][0156]
此时，对权重进行更新，即重新计算方差与似然函数值l
(iter+1
)，计算公式分别为：
[0157][0158][0159]
重复上述操作直至l
(iter+1
)-l
(iter)
≤0.001，此时迭代停止，有参数率定完成；
[0160]
(c)假定步骤s21中优选的单值预测方案组合内有候选方案k组，根据各候选方案的bma权重[ω1,ω2,...,ωk]，定义累积权重值
[0161][0162]
生成[0,1]均匀随机数u，如果则表示此次抽样选取的是第i个候选方案fi；
[0163]
(d)根据fi在时段j的概率分布随机生成服从正态分布的径流预测序列y'；
[0164]
(e)基于公式(13)对生成的服从正态分布的径流预测序列y'进行还原，计算公式为：
[0165][0166]
(f)重复步骤(c)-(e)m次，即可得到一组具有m个情景的预测情景集，即bma的径流概率预测方案，其均值过程即为bma的径流单值预测方案。二者共同构成包含单值过程与概率区间的径流集合预测方案，如图3中(b)所示，考虑预测因子与模型结构优选的预测方案能够在尽可能覆盖实测真值的前提下降低区间宽度，具有较高的精准性，但其在预测失准情形下的预测偏差过大，可靠性有待提高。
[0167]
s3、参数优选：考虑误差校正的径流概率预测方案生成；
[0168]
作为“输入-结构-参数”全要素分层组合优选的第三层，对模型参数进行优选：基于广义自回归条件异方差模型(garch)构建径流集合预测误差的时变模拟模型，利用误差的相依关系进行实时误差校正，结合蒙特卡洛抽样生成校正后的误差情景序列，叠加至bma单值预报结果上得到考虑误差校正的径流概率预测方案，依据径流预测“精准度-可靠性-稳定性”多维评价指标体系对得到的概率预测方案进行评价。
[0169]
具体包括以下步骤：
[0170]
s31、基于garch模型的径流预测误差精细表征；
[0171]
广义自回归条件异方差模型(garch)是一种针对金融数据所量体订做的回归模型，除去和普通回归模型(如arma、arch等)相同之处外，garch对误差的方差进行了进一步的建模，特别适用于具有波动聚集性的时间序列的分析和预测，这样的分析对投资者的决
策能起到非常重要的指导性作用。而径流预测的方差往往也存在这波动聚集性，且具有“尖峰、厚尾”的统计分布特性，与金融序列的规律高度一致，因此可用于描述径流预测误差序列的时程变化情况。
[0172]
garch模型一般由条件均值和条件方差两个基本方程来描述，条件均值由步骤s22中的单值预测过程给出，故只需建立方差(即误差项)的模型，从而达到精细表征径流预测误差的目的。假定残差项ε
t
与条件方差h
t
服从：
[0173][0174]
式中，c是常数项，ai和bj分别是滞后i阶与j阶的条件方差h
t-i
和残差项ε
t-j
的系数，m与n分别为garch条件方差与残差项的滞后阶数，可记为ε
t
～garch(m,n)。采用极大似然估计对参数ai、bj进行了率定。
[0175]
s32、基于蒙特卡洛抽样的径流预测误差校正；
[0176]
基于步骤s31中得到的garch模型及其各参数，可结合蒙特卡洛抽样生成径流预测误差模拟情景序列。对于每一个时刻t，首先基于蒙特卡罗方法随机生成n组服从标准正态分布的随机数v
t
，代入公式(16)计算得到条件方差h
t
的模拟情景序列，将条件方差项开方后即可得到误差模拟情景序列，该序列即为径流预测误差的校正。
[0177]
s33、考虑误差校正的径流概率时变预测方案生成；
[0178]
一般地，预测误差与实测值的关系为：
[0179][0180]
式中，y
t
为时段t的预测值，为时段t的实测真值，e
t
为该时段对应的预测误差。若考虑误差校正，即当e
t
'＝z(e
t
)时，有：
[0181][0182][0183]
式中，e
t
'为校正后的径流预测误差，z为误差校正函数，反映历史预测误差的统计特性；为径流“理论值”，即考虑误差校正后对径流实测值的精细预估。结合步骤s32中生成的径流预测误差模拟情景序列，即e
t
'，将其叠加在步骤s22中bma的径流单值预测过程上，即可得到径流概率预测情景序列，即考虑误差校正的径流时变概率预测方案，如图3中(a)所示，考虑全要素优选后的预测方案能够在尽可能覆盖实测真值的前提下进一步降低区间宽度，同时对于预测失准的情形也有更小的预测偏差，具有较高的稳定性。
[0184]
s34、径流预测精度评价；
[0185]
一般地，概率预测结果的检验要从精度与可靠度两方面着手。精度是指概率预测结果分布的某一分位数(中位数q50或均值)与实测值的接近程度，一个精准的径流预测方案满足：(a)预测误差无偏；(b)预测误差方差尽可能小，可采用均方根误差(rmse)进行评价；可靠度是指概率预测结果的在一定置信水平下的置信区间覆盖实测值的能力，通常采
用覆盖率(cr)以及brier评分(bs)进行评价。此外，从风险管理与系统安全性的角度出发，好的预测方案应具有较强的抗风险能力，即稳定性，当概率预测失准(即给定的置信水平下概率预测的区间未覆盖实测真值)时，其区间阈值距离实测值越近，该预测方案为决策者提供的有效信息越多，抗风险能力越强。失准偏差(md)考虑预测失准情形下的最小预测偏差，可用于表征径流概率预测的稳定性。因此，本发明构建了径流预测“精准度-可靠性-稳定性”多维评价指标，选取rmse、bs、md分别为评价径流概率预测精准度、可靠性与稳定性的指标来衡量预测方案的优劣，对预测方案进行全面、系统的评价。具体的计算公式为：
[0186][0187][0188][0189]
式中，α为置信水平，t为时段数，p为分位数，probi与oi分别为第i次预测发生的概率及其实际情况(如果预测结果在许可误差范围内，则oi＝1，反之则oi＝0)，与分别为第t个时段径流预测的许可误差的上下限，一般取相对误差
±
30％，与分别为第t个时段径流概率预测区间的上下限。
[0190]
本发明的一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测系统，包括：
[0191]
模型输入优选模块，用于收集研究区长系列中长期气象、水文因子数据资料构建径流预测初始因子集，基于copula熵与径流成因分析筛选预测因子组合得到驱动因子集，采用基于深度学习的长短期记忆神经网络模型生成中长期径流单值预测方案集，作为“输入-结构-参数”全要素分层组合优选的第一层；
[0192]
模型结构优选模块，用于以信息熵作为评价单值预测方案集信息价值的指标，采用具有内外层循环嵌套结构的逐步回归对方案集中的方案组合进行优选，对最优子集采用贝叶斯模型平均的方式生成同时包含单值过程与概率区间的径流集合预测方案，作为“输入-结构-参数”全要素分层组合优选的第二层；
[0193]
模型参数优选模块，用于基于广义自回归条件异方差模型构建径流集合预测误差的时变模拟模型，利用误差的相依关系进行实时误差校正，结合蒙特卡洛抽样生成校正后的误差情景序列，叠加至贝叶斯模型平均单值预报结果上得到考虑误差校正的径流概率预测方案，依据径流预测“精准度-可靠性-稳定性”多维评价指标体系对得到的概率预测方案进行评价，作为“输入-结构-参数”全要素分层组合优选的第三层。
[0194]
本发明的一种装置设备，包括存储器和处理器，其中：
[0195]
存储器，用于存储能够在处理器上运行的计算机程序；
[0196]
处理器，用于在运行所述计算机程序时，执行如上述一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法的步骤，并达到如上述方法一致的技术效果。
[0197]
本发明的一种存储介质，所述存储介质上存储有计算机程序，所述计算机程序被至少一个处理器执行时实现如上述一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法的步骤，并达到如上述方法一致的技术效果。

技术特征：

1.一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法，其特征在于，包括以下步骤：s1、作为“输入-结构-参数”全要素分层组合优选的第一层，对模型输入进行优选：收集研究区长系列中长期气象、水文因子数据资料构建径流预测初始因子集，基于copula熵与径流成因分析筛选预测因子组合得到驱动因子集，采用基于深度学习的长短期记忆神经网络lstm模型生成中长期径流单值预测方案集；s2、作为“输入-结构-参数”全要素分层组合优选的第二层，对模型结构进行优选：以信息熵作为评价单值预测方案集信息价值的指标，采用具有内外层循环嵌套结构的逐步回归对单值预测方案集中的方案组合进行优选，对最优子集采用贝叶斯模型平均bma的方式生成同时包含单值过程与概率区间的径流集合预测方案；s3、作为“输入-结构-参数”全要素分层组合优选的第三层，对模型参数进行优选：基于广义自回归条件异方差模型garch对径流预测误差进行精细表征，利用误差的相依关系进行实时误差校正，结合蒙特卡洛抽样生成校正后的误差情景序列，叠加至步骤s2中的单值过程上得到考虑误差校正的径流时变概率预测方案，依据径流预测“精准度-可靠性-稳定性”多维评价指标体系对得到的概率预测方案进行评价。2.根据权利要求1所述的一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法，其特征在于，步骤s1包括以下步骤：s11、径流预测初始因子集生成；根据研究区域气象条件，确定控制研究区域降雨、径流的主要天气要素，选取相应的指标与前期降雨、蒸散发量以及径流量共同组成径流预测初始因子集；s12、基于copula熵与径流成因分析的驱动因子识别；设n维随机变量(x1,x2,...,x
n
)具有边缘分布u
i
＝f
i
(x
i
),i＝1,2,...,n，其中f为累积概率密度函数，copula熵定义为：其中，ce为copula熵，c为copula函数，基于copula熵对所有与径流成因有关的预测因子进行关联度排序，选取排序靠前的若干因子作为驱动因子集；s13、基于lstm模型的径流单值预测方案集生成；采用lstm对各月的入湖径流进行预测，以不同因子组合作为因变量进行学习、训练、模拟并进行预测，获得单值预测方案集。3.根据权利要求2所述的一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法，其特征在于，步骤s13中lstm的计算方法为：(a)遗忘门：f
t
＝σ(w
f
·
[h
t-1
,x
t
]+b
f
)其中，f
t
为当前时段的遗忘门门控，表示控制部分历史输出信息转换为当前时段的输入，反映历史信息在当前时段的丢失程度，w
f
为遗忘门的权重矩阵，b
f
为w
f
的偏置项，σ为激活函数；(b)输入门：i
t
＝σ(w
i
·
[h
t-1
,x
t
]+b
i
)
其中，i
t
为当前时段的输入门门控，表示控制部分当前时段的输入转化为系统状态变量，反映当前时段的信息输入效率，w
i
为输入门的权重矩阵，b
i
为w
i
的偏置项；(c)当前输入的长期状态c
t
'：c
t
'＝tanh(w
c
·
[h
t-1
,x
t
]+b
c
)其中，w
c
为当前输入的长期状态的权重矩阵，b
c
为w
c
的偏置项，tanh为激活函数；(d)当前时刻的长期状态c
t
：c
t
＝f
t
·
c
t-1
+i
t
·
c
t
'(e)输出门：o
t
＝σ(w
o
·
[h
t-1
,x
t
]+b
o
)其中，o
t
为当前时段的输出门门控，表示控制部分系统状态变量转化为当前时段的输出，反映当前时段的信息输出效率，w
o
为输出门的权重矩阵，b
o
为w
o
的偏置项；(f)lstm的输出值h
t
：h
t
＝o
t
·
tanh(c
t
)将计算时段划分为训练期与校验期，将训练期中经步骤s12中识别的驱动因子集w作为上述计算步骤中的输入项x
t
，将径流量作为其中的输出项h
t
，经神经网络的训练后率定得到参数w
f
、b
f
、wi、b
i
、w
o
、b
o
，并与校验期的输入项x
t
一同代入模型中计算得到对应的输出项h
t
，得到一组单值预测过程，构成单值预测方案集。4.根据权利要求1所述的一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法，其特征在于，步骤s2包括以下步骤：s21、基于信息熵与逐步回归的方案集组合优选；信息熵的表达式为：其中，h为信息熵，x为一个具有n个离散取值的随机变量，x＝{x1,x2,...,x
n
}，其概率分布列此外，对于一个二元随机变量(x,y)，其信息熵又称为联合熵，定义为：其中，p(x,y)为随机变量x与y的联合概率密度函数；联合熵表示多个随机变量包含的信息总量，而其信息同量被定义为互信息，即：t(x,y)＝h(x)+h(y)-h(x,y)因此，用信息熵中的联合熵与互信息分别描述单值预测方案集中各候选方案间的信息总量与信息同量，并以此作为方案组合优选的依据；s22、基于bma的径流单值-概率集合预测方案生成；设y为预测变量，y
obs
为实测样本序列，f＝{f1,f2,...,f
k
}为候选模型空间，其中k表示步骤s21中优选的单值预测方案组合内的候选方案组数；根据贝叶斯理论，假定实测与预报序列均服从正态分布，进一步细化：
其中，p(y|y
obs
)表示预测变量的概率分布，ω
i
表示第i个模型的权重，表示均值为f
i
，方差为的正态分布；采用期望最大化算法率定bma参数ω
i
和并采用monte-carlo抽样给出概率预测情景序列，获得包含单值过程与概率区间的径流集合预测方案。5.根据权利要求4所述的一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法，其特征在于，步骤s22中采用monte-carlo抽样给出概率预测情景序列，具体步骤如下：(a)采用box-cox变换将非正态分布的径流序列转化成正态分布，具体为：其中，λ为转换系数，其取值由极大似然估计给出；y'为经过转化后的服从正态分布的径流预测序列；(b)采用e-m算法率定bma参数ω
i
和(c)根据各候选方案的bma权重[ω1,ω2,...,ω
k
]，定义累积权重值]，定义累积权重值生成[0,1]均匀随机数u，如果则表示此次抽样选取的是第i个候选方案f
i
；(d)根据f
i
在时段j的概率分布随机生成服从正态分布的径流预测序列y'；(e)基于累积权重值公式对生成服从正态分布的径流预测序列y'进行还原，计算公式为：(f)重复步骤(c)-(e)m次，得到一组具有m个情景的预测情景集，即bma的径流概率预测方案，其均值过程即为bma的径流单值预测方案；二者共同构成径流单值-概率集合预测方案。6.根据权利要求1所述的一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法，其特征在于，步骤s3包括以下步骤：s31、基于garch模型的径流预测误差精细表征；garch模型由条件均值和条件方差两个基本方程来描述，条件均值由步骤s2中的单值预测过程给出，条件方程的计算为：
假定残差项ε
t
与条件方差h
t
服从：其中，c是常数项，a
i
和b
j
分别是滞后i阶与j阶的条件方差h
t-i
和残差项ε
t-j
的系数，m与n分别为garch条件方差与残差项的滞后阶数，记为ε
t
～garch(m,n)；采用极大似然估计对参数a
i
、b
j
进行率定；s32、基于蒙特卡洛抽样的径流预测误差校正；基于步骤s31中得到的garch模型及其各参数，结合蒙特卡洛抽样生成径流预测误差模拟情景序列，对于每一个时刻t，首先基于蒙特卡罗方法随机生成n组服从标准正态分布的随机数v
t
，代入步骤s31中公式计算得到条件方差h
t
的模拟情景序列，将条件方差项开方后得到误差模拟情景序列，该序列即为径流时变预测误差的校正；s33、考虑误差校正的径流概率时变预测方案生成；预测误差与实测值的关系为：其中，y
t
为时段t的预测值，为时段t的实测真值，e
t
为该时段的预测误差；若考虑误差校正，即当e
t
'＝z(e
t
)时，有：)时，有：其中，e
t
'为校正后的径流预测误差，z为误差校正函数，反映历史预测误差的统计特性；为径流“理论值”，即考虑误差校正后对径流实测值的精细预估；结合步骤s32中生成的径流预测误差模拟情景序列，即e
t
'，将其叠加在步骤s22中bma的径流单值预测过程上，得到径流概率预测情景序列，即考虑误差校正的径流概率预测方案；s34、径流预测精度评价；构建径流预测“精准度-可靠性-稳定性”多维评价指标，选取均方根误差rmse、brier评分bs、失准偏差md分别为评价径流概率预测精准度、可靠性与稳定性的指标来衡量预测方案的优劣，对预测方案进行全面、系统的评价。7.一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测系统，其特征在于，包括：模型输入优选模块，用于收集研究区长系列中长期气象、水文因子数据资料构建径流预测初始因子集，基于copula熵与径流成因分析筛选预测因子组合得到驱动因子集，采用基于深度学习的长短期记忆神经网络模型生成中长期径流单值预测方案集，作为“输入-结构-参数”全要素分层组合优选的第一层；模型结构优选模块，用于以信息熵作为评价单值预测方案集信息价值的指标，采用具有内外层循环嵌套结构的逐步回归对方案集中的方案组合进行优选，对最优子集采用贝叶
斯模型平均的方式生成同时包含单值过程与概率区间的径流集合预测方案，作为“输入-结构-参数”全要素分层组合优选的第二层；模型参数优选模块，用于基于广义自回归条件异方差模型构建径流集合预测误差的时变模拟模型，利用误差的相依关系进行实时误差校正，结合蒙特卡洛抽样生成校正后的误差情景序列，叠加至贝叶斯模型平均单值预报结果上得到考虑误差校正的径流概率预测方案，依据径流预测“精准度-可靠性-稳定性”多维评价指标体系对得到的概率预测方案进行评价，作为“输入-结构-参数”全要素分层组合优选的第三层。8.一种装置设备，其特征在于，包括存储器和处理器，其中：存储器，用于存储能够在处理器上运行的计算机程序；处理器，用于在运行所述计算机程序时，执行如权利要求1-7任一项所述一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法的步骤。9.一种存储介质，其特征在于，所述存储介质上存储有计算机程序，所述计算机程序被至少一个处理器执行时实现如权利要求1-7任一项所述一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法的步骤。

技术总结

本发明公开了一种“输入-结构-参数”全要素分层组合优选的中长期径流时变概率预测方法，包括构建径流预测因子集，对预测因子筛选并自由组合得到驱动因子集，生成中长期径流单值预测方案集；以信息熵作为评价单值预测方案集信息价值的指标，采用具有内外层循环嵌套结构的逐步回归对方案集中的方案组合进行优选，对最优子集生成同时包含单值过程与概率区间的径流集合预测方案；对径流集合预测误差进行精细表征并生成误差情景序列，得到考虑误差校正的径流时变概率预测方案。本发明提出全要素分层优选的中长期径流预测方法，为水资源规划与管理提供精准、可靠、稳定的预报方案信息。稳定的预报方案信息。稳定的预报方案信息。