CN109344999A - 一种径流概率预报方法 - Google Patents

Abstract

本发明公开一种径流概率预报方法，其中，该方法主要包括：采用基于K‑medoids的聚类方法对训练集进行聚类，得到隐含马尔科夫模型HMM的初始化参数；利用Baum–Welch算法对HMM进行学习，得到HMM的状态转移概率矩阵以及观测模型的概率分布；根据贝叶斯信息准则BIC进行模型选择，选择适宜该训练集的HMM状态个数；最终根据给出的预报因子，结合高斯混合回归GMR推理得到条件概率分布函数作为径流概率预报。本发明的概率预报方法引入了径流隐含状态的概念，可利用水文、地形、气象等诸多因素训练得到隐含状态转移概率矩阵，得到有效可靠的未来径流概率预报分布，为水库优化调度决策提供科学依据。

Description

一种径流概率预报方法

技术领域

本发明涉及水文水资源领域，更具体地，涉及一种径流概率预报方法。

背景技术

径流预报能够在防洪、供水、发电等水库调度运行管理的各个方面发挥重要作用。高精度并且可靠的预报能为流域梯级水库优化调度决策提供科学依据，对流域水资源合理利用具有重大意义。然而，降雨径流的形成过程受到水文、地形、气象等诸多自然因素的影响，呈现高度的非线性、随机性和不确定性特征，给径流预测带来了困难，难以采用传统的方法开展径流预测。因此，如何构建考虑诸多影响径流过程的隐含因素的预测模型是亟需解决的理论和实际工程问题。

隐含马尔科夫模型(hidden Markov model，HMM)由隐含状态序列以及其对应的观测模型构成，其中隐含状态序列是一个典型的马尔科夫链，具有无后效性，对应的观测模型满足独立输出假设：观测模型跟当前隐含状态相关且仅跟当前隐含状态相关。可将隐含马尔科夫模型运用于水文的径流预报当中，视前期径流、土壤、气象等因素为隐含状态序列，对应的径流过程为观测概率模型，能够充分模拟具有非线性、随机性和不确定性的历史径流过程。如何通过训练后的隐含马尔科夫模型来预测未来的径流序列是目前所需解决的问题。

发明内容

针对现有技术的缺陷，本发明的目的在于解决现有降雨径流的形成过程受到水文、地形、气象等诸多自然因素的影响，呈现高度的非线性、随机性和不确定性特征，给径流预测带来了困难，难以采用传统的方法开展径流预测，以及如何通过训练后的隐含马尔科夫模型来预测未来的径流序列的技术问题。

为实现上述目的，本发明提供一种径流概率预报方法，包括如下步骤：

(1)构建由预报径流和其预报因子组成的数据集，对数据进行归一化处理得到训练集；

(2)基于K-medoids的聚类方法对训练集进行聚类，根据聚类后得到的各类的均值向量、协方差矩阵以及径流状态转移概率矩阵作为隐含马尔科夫模型的初始化参数，所述隐含马尔科夫模型包括隐含状态序列及其对应的观测模型；

(3)利用Baum–Welch算法对隐含马尔科夫模型进行学习，得到训练后的隐含马尔科夫模型径流状态转移概率矩阵以及观测模型的概率分布；

(4)设置不同的隐含状态个数K＝1～15；重复步骤(2)-(3)，并计算相应的贝叶斯信息准则值；

(5)选择贝叶斯信息准则值最小值所对应的K值，作为隐含马尔科夫模型隐含个数；

(6)根据训练得到的隐含马尔科夫模型，结合高斯混合回归(Gaussian MixtureRegression，GMR)方法，利用观测模型的概率分布以及测试集所给的预报因子得到各径流隐含状态下径流条件概率分布，利用径流状态转移概率矩阵求解测试集中每个时间节点落在K个径流隐含状态下的概率，综合推导得到测试集的条件概率分布函数，作为对测试集的径流概率预报结果。

可选地，步骤(2)中K-medoids算法使用核函数代替欧几里得距离来度量两个数据点之间的相似性；

从训练集包括的数据点中随机选择K个中心点，当更新中心点时，每个数据点被分类为与中心点最相似的类别，并且计算每个数据点与其相同类别所有数据点的相似度，选择相似度最高的数据点作为此类别的新的中心点。

可选地，步骤(3)包括：

E步骤：根据现有的隐含马尔科夫模型，计算各个观测数据输入到模型中的计算结果，确定完全数据的对数似然函数：

其中，Q(θ,θ^old)表示对数似然函数，θ表示隐含马尔科夫模型的所有参数集θ＝{π,A,φ}，θ^old表示更新前的参数值，t表示数据集的索引，1≤t≤T，T表示数据集个数，π_k表示隐含马尔科夫模型在k状态的初始概率，1≤k≤K，K表示隐含状态个数，A_ij表示从状态i转移到状态j的概率，φ_k表示第k个观测概率模型的参数，p(x_t|φ_k)表示第t个观测数据点在第k个观测概率模型下的条件概率值；

γ_t(k)表示观测数据在时刻t处于第k个状态变量的概率，其概率计算公式如下：

其中，X表示所有观测数据，z_t表示t时刻的状态；

ξ_t(i,j)表示观测数据在时刻t-1处于第i个状态且在时刻t处于第j个状态的概率，其概率计算公式如下：

M步骤：根据以上E步骤得到的参状态参数γ_t(k)和ξ_t(i,j)，对隐含马尔科夫模型的参数θ＝{π,A,φ}进行重估计，假设隐含马尔科夫模型中观测概率模型为高斯分布，则φ_k＝{μ_k,∑_k}，新的参数估计值由下式给出：

可选地，步骤(4)中贝叶斯信息准则值通过以下公式计算：

其中，表示模型的极大似然值，n_p表示模型参数个数，隐含状态个数不同则对应的模型极大似然数和模型参数个数均不同。

可选地，步骤(6)包括：

在预报的过程中，对于每个隐含状态，观测概率模型的均值向量μ_k将按照预报因子、预报变量进行拆分：

其中，为表示预报因子均值的子向量，为表示预报变量均值的子向量，与此同时，协方差矩阵也被拆分：

其中，表示预报因子之间的协方差矩阵，表示预报因子与预报变量的协方差矩阵，表示预报变量与预报因子的协方差矩阵，表示预报变量之间的协方差矩阵；

根据联合正态分布的性质推理，给定预报因子y_t后，预报值的条件概率服从联合正态分布，其均值向量和协方差矩阵如下：

在原始的高斯混合回归框架中，不同高斯观测模型的权重h_k表示为每个观测值属于第k个高斯观测模型的比重；

通过递归计算来估计此权重值，从而不仅考虑预报因子，而且还考虑封装在隐含马尔科夫模型中的时间序列信息：

其中，h_k(y_t)表示隐含马尔科夫模型前向变量，表示高斯分布的条件概率值；

对于隐含马尔科夫模型，具有K个隐含状态对应的高斯观测模型，因此，给定预报因子y_t的情况下，预报变量z_t的条件概率分布函数p(z_t|y_t)为：

其中，表示第k个隐含状态下z_t的条件概率分布函数。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，具有以下有益效果：

本发明提供的一种基于隐含马尔科夫模型结合高斯混合回归的径流预报方法，将隐含马尔科夫模型运用与水文的径流预报当中，视前期径流、土壤、气象等因素为隐含状态序列，对应的径流过程为观测概率模型，能够充分的模拟具有非线性、随机性和不确定性的历史径流过程。混合高斯回归是一种概率密度估计的方法，其主要作用是通过联合概率密度函数衍生出该联合概率密度函数的条件概率密度函数作为预报概率。研究将此权重值推广到隐含马尔科夫模型中，通过递归计算来估计权重值，不仅考虑预报因子，而且还考虑封装在HMM中的时间序列信息。

附图说明

图1为本发明提供的隐含马尔科夫模型结合高斯回归预报方法整体流程图；

图2为本发明提供的小湾水库来水预报值、90％预报区间与实际径流对比示意图；

图3为本发明提供的小湾水库来水预报概率积分转换(PIT)统计图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

混合高斯回归是一种概率密度估计的方法，其主要作用是通过联合概率密度函数衍生出该联合概率密度函数的条件概率密度函数作为预报概率，一般运用于混合模型的预测中，而隐含马尔科夫模型可以视为一种条件混合模型，因此可将混合高斯回归拓展到隐含马尔科夫模型作为预测方法。

本发明提供一种基于隐含马尔科夫模型结合高斯混合回归的径流预报方法，可以有效的考虑影响径流过程的各种因素，具有较高的预报精度以及预报可靠性。

图1所示为本发明基于隐含马尔科夫模型结合高斯混合回归的径流预报方法总体流程图，具体包括以下步骤：

(1)构建由预报径流y和其预报因子z组成的数据集x＝[y,z]，对数据进行归一化处理。

(2)采用基于K-medoids的聚类方法对训练集进行聚类，与传统的K-means方法不同，K-medoids算法使用核函数k(x_i,x_i′)代替欧几里得距离||x_i-x′_i||²用于度量两个数据点之间的相似性，本发明采用的是径向基(RBF)核函数：

其中，x_i表示第i个数据，x′_i表示第i'个数据，d(i,i′)表示第i个数据与第i'个数据的相似性也即核函数值，σ表示RBF的参数。

首先，从原始数据点中随机选择K个中心点m_1:K from{1,…,N}，当更新中心点时，每个点被分类为与中心点最相似的类别，并且计算它与相同类别所有点的相似度，选择相似度最高(与其他所有点的核函数值之和最大)的点作为此类别的新的中心点m_k：

其中，z′_i表示第i'数据的类别。

根据聚类后的各类的均值向量μ_k、协方差矩阵∑_k以及状态转移概率矩阵A得到隐含马尔科夫模型的初始化参数。

(3)利用Baum–Welch算法对隐含马尔科夫模型进行学习，得到隐含马尔科夫模型的径流状态转移概率矩阵A以及观测模型的概率分布p(x_t|s_t＝k,φ_k)；

E过程：根据现有的模型，计算各个观测数据输入到模型中的计算结果，首先确定完全数据的对数似然函数：

γ_t(k)表示给定模型和观测数据，在时刻t处于第k个状态变量的概率，其概率计算公式如下：

ξ_t(i,j)表示给定模型和观测数据，在时刻t-1处于第i个状态且在时刻t处于第j个状态的概率，其概率计算公式如下：

M过程：M过程是重新估计参数。

对于观测概率模型为高斯分布的隐含马尔科夫模型来说，为了得到每个状态的均值μ_k和Σ_k项的M步，新的参数估计值由下式给出：

(4)设置不同的隐含状态个数K＝1～15；重复步骤(2)-(3)，并计算相应的贝叶斯信息准则值BIC：

(5)选择贝叶斯信息准则值最小值所对应的K值，作为径流预测模型的隐含个数；

(6)根据训练得到的隐含马尔科夫模型，结合高斯混合回归方法，推导得到测试集的条件概率分布函数，作为概率预报结果。

隐含马尔科夫模型也称为非独立混合模型，它可以解释为一个混合模型的扩展，其中每个观测的混合分布概率不是独立的，而是取决于先前观察的潜在状态变量。在预报的过程中，隐含马尔科夫模型的所有变量被区分为预报因子(前期径流和气象指标)和预报变量(未来径流)x_t＝[y_t,z_t]。相应的对于每个隐含状态，观测概率模型的均值向量将按照预报因子、预报变量进行拆分：

上式中，为表示预报因子均值的子向量，为表示预报变量均值的子向量，与此同时，协方差矩阵也被拆分：

根据联合正态分布的性质推理，给定预报因子y_t后，预报值的条件概率服从联合正态分布，其均值向量和协方差矩阵如下：

在原始的高斯混合回归框架中，不同高斯观测模型的权重h_k表示，为每个观测值属于第k个高斯观测模型的比重，此权重值并未考虑观测序列的顺序信息。研究将此权重值推广到隐含马尔科夫模型中，通过递归计算来估计此权重值，从而不仅考虑预报因子，而且还考虑封装在HMM中的时间序列信息：

对于隐含马尔科夫模型，具有K个隐含状态对应的高斯观测模型。因此，给定预报因子y_t的情况下，预报变量z_t的条件概率分布函数p(z_t|y_t)为：

本发明以澜沧江下游小湾电的月平均径流量时间序列为对象，数据采用1953年到2010年月径流数据。以前1月径流、前12月径流以及南方涛动指数和亚洲经向环流指数为预报因子，按照图1所示的基于隐含马尔科夫模型结合高斯混合回归的径流预报方法进行预报。

表1给出小湾电站入库径流概率预报的评分表，其中，表1中SS_CRPS表示连续分级概率评分，SS_MSE表示均方误差评分，SS_RMSE表示均方根误差评分，SS_MAE表示平均绝对误差评分。

表1概率预报评分表

所有评分SS均与历史重采样预报对比，若预报精度高于历史重采样预报则为正值，若预报精度低于历史重采样预报则为负值，评分范围在-1～1之间，评分越高代表预报精度越好。其中SS_CRPS为连续分级概率评分，其同属考虑了预报的偏差以及概率预报的范围，是一种概率预报常用的评价指标；SS_MSE为均方误差评分，其评分公式与水文预报常用的评价指标“纳什效率系数”相同；SS_RMSE为均方根误差评分即Brier评分，再概率预报中广泛使用；SS_MAE为平均绝对误差评分。从表中可以看出，隐含马尔科夫模型的预报评分均大于0，代表研究所提出的概率预报精度均高于历史重采样预报。从表中看出：其中12月～4月的评分普遍比5月～11月的评分要高。这一结果表明，月预报精度与降雨的季节性密切相关：枯水期径流主要取决于流域“干湿”情况，以前期径流作为预报因子，能够较好的进行预报；相比而言，汛期径流主要取决于降水情况，由于降水可预报性很低，汛期径流预报精度也不高。

图2展示了本预报方法在测试期间的预报值、90％预报区间与实际月径流的对比图，图中可以看出本发明的预报方法具有较高的精度。

在预报精度评价的基础上，进一步对预报可靠性进行分析。分析观测值相对于预报概率积分转换(PIT)值，如图3所示。可以看到，概率积分转换整体上呈现均匀分布的特征。这一结果说明：观测值基本上可以看作来自对应概率预报的随机采样。这一结果表明，此预报具有较高的可靠性。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种径流概率预报方法，其特征在于，包括如下步骤：

(1)构建由预报径流和其预报因子组成的数据集，对数据进行归一化处理得到训练集；

(2)基于K-medoids的聚类方法对训练集进行聚类，根据聚类后得到的各类的均值向量、协方差矩阵以及径流状态转移概率矩阵作为隐含马尔科夫模型的初始化参数，所述隐含马尔科夫模型包括隐含状态序列及其对应的观测模型；

(3)利用Baum–Welch算法对隐含马尔科夫模型进行学习，得到训练后的隐含马尔科夫模型径流状态转移概率矩阵以及观测模型的概率分布；

(4)设置不同的隐含状态个数K＝1～15；重复步骤(2)-(3)，并计算相应的贝叶斯信息准则值；

(5)选择贝叶斯信息准则值最小值所对应的K值，作为隐含马尔科夫模型隐含个数；

(6)根据训练得到的隐含马尔科夫模型，结合高斯混合回归方法，利用观测模型的概率分布以及测试集所给的预报因子得到各径流隐含状态下径流条件概率分布，利用径流状态转移概率矩阵求解测试集中每个时间节点落在K个径流隐含状态下的概率，综合推导得到测试集的条件概率分布函数，作为对测试集的径流概率预报结果。

2.根据权利要求1所述的径流概率预报方法，其特征在于，所述步骤(2)中K-medoids算法使用核函数代替欧几里得距离来度量两个数据点之间的相似性；

从训练集包括的数据点中随机选择K个中心点，当更新中心点时，每个数据点被分类为与中心点最相似的类别，并且计算每个数据点与其相同类别所有数据点的相似度，选择相似度最高的数据点作为此类别的新的中心点。

3.根据权利要求1所述的径流概率预报方法，其特征在于，所述步骤(3)包括：

E步骤：根据现有的隐含马尔科夫模型，计算各个观测数据输入到模型中的计算结果，确定完全数据的对数似然函数：

其中，Q(θ,θ^old)表示对数似然函数，θ表示隐含马尔科夫模型的所有参数集θ＝{π,A,φ}，θ^old表示更新前的参数值，t表示数据集的索引，1≤t≤T，T表示数据集个数，π_k表示隐含马尔科夫模型在k状态的初始概率，1≤k≤K，K表示隐含状态个数，A_ij表示从状态i转移到状态j的概率，φ_k表示第k个观测概率模型的参数，p(x_t|φ_k)表示第t个观测数据点在第k个观测概率模型下的条件概率值；

γ_t(k)表示观测数据在时刻t处于第k个状态变量的概率，其概率计算公式如下：

其中，X表示所有观测数据，z_t表示t时刻的状态；

ξ_t(i,j)表示观测数据在时刻t-1处于第i个状态且在时刻t处于第j个状态的概率，其概率计算公式如下：

M步骤：根据以上E步骤得到的参状态参数γ_t(k)和ξ_t(i,j)，对隐含马尔科夫模型的参数θ＝{π,A,φ}进行重估计，假设隐含马尔科夫模型中观测概率模型为高斯分布，则φ_k＝{μ_k,∑_k}，新的参数估计值由下式给出：

4.根据权利要求3所述的径流概率预报方法，其特征在于，所述步骤(4)中贝叶斯信息准则值通过以下公式计算：

其中，表示模型的极大似然值，n_p表示模型参数个数，隐含状态个数不同则对应的模型极大似然数和模型参数个数均不同。

5.根据权利要求1所述的径流概率预报方法，其特征在于，所述步骤(6)包括：

在预报的过程中，对于每个隐含状态，观测概率模型的均值向量μ_k将按照预报因子、预报变量进行拆分：

其中，为表示预报因子均值的子向量，为表示预报变量均值的子向量，与此同时，协方差矩阵也被拆分：

其中，表示预报因子之间的协方差矩阵，表示预报因子与预报变量的协方差矩阵，表示预报变量与预报因子的协方差矩阵，表示预报变量之间的协方差矩阵；

根据联合正态分布的性质推理，给定预报因子y_t后，预报值的条件概率服从联合正态分布，其均值向量和协方差矩阵如下：

在原始的高斯混合回归框架中，不同高斯观测模型的权重h_k表示为每个观测值属于第k个高斯观测模型的比重；

通过递归计算来估计此权重值，从而不仅考虑预报因子，而且还考虑封装在隐含马尔科夫模型中的时间序列信息：

其中，h_k(y_t)表示隐含马尔科夫模型前向变量，表示高斯分布的条件概率值；

对于隐含马尔科夫模型，具有K个隐含状态对应的高斯观测模型，因此，给定预报因子y_t的情况下，预报变量z_t的条件概率分布函数p(z_t|y_t)为：

其中，表示第k个隐含状态下z_t的条件概率分布函数。