CN111860974B - 一种基于状态空间和联合分布的干旱多级预测方法 - Google Patents

Abstract

本发明提供了一种基于状态空间和联合分布的干旱多级预测方法，包括：首先提取干旱指数SPI时间序列，利用一阶状态空间模型解译干旱指数SPI时间序列的状态转移矩阵，并计算待预测时段的干旱状态概率矩阵；然后提取待预测时段前一时段的相关因子，构建每种干旱状态下的相关因子和干旱指数SPI时间序列的联合概率分布，推求相关因子给定情况下的干旱指数的条件概率模型；最后，计算待预测区域在待预测时段T的干旱指数的全状态分布；基于所述干旱指数的全状态分布，计算得到干旱指数的均值和90％概率区间。本发明的有益效果是：状态空间方法能反映干旱不确定性本质；通过多重预测和分类别引入影响因子能提高干旱预测的精度；过程考虑更加全面和精细化。

Description

一种基于状态空间和联合分布的干旱多级预测方法

技术领域

本发明涉及干旱预测技术领域，尤其涉及一种基于状态空间和联合分布的干旱多级预测方法。

背景技术

近几年我国云南、四川、重庆、湖南、湖北、安徽等省市接连遭受了五十年一遇或百年一遇的大旱，造成难以估量的经济损失。研究一种可靠的干旱预测模型，提前预估干旱事件的发生时间和影响程度将对社会经济可持续发展有重要意义。当前，最常用的干旱指数是标准化降雨指数(SPI)，对SPI指数预报的方法包括时间序列方法、多元统计回归方法、机器学习方法等。这些模型和理论简单可行，但是对干旱状态的延续性和概率性本质认识不足，预测精度也因此较为有限。

发明内容

为了解决现有技术存在的干旱指数的预测精度低以及提供可用信息较少的技术问题，本发明提供了一种基于状态空间和联合分布的干旱多级预测方法，主要包括以下步骤：

S101：提取待预测区域的历史月降雨数据，并计算历史月降雨数据对应的标准化降雨指数SPI，形成干旱指数SPI时间序列；

S102：利用一阶状态空间模型解译干旱指数SPI时间序列，得到状态转移矩阵；

S103：根据所述状态转移矩阵，计算待预测区域在待预测时段处于各干旱状态的概率，得到待预测时段的干旱状态概率矩阵；

S104：对每一种干旱状态对应的SPI值，提取待预测时段前一时段的相关因子，得到各干旱状态下的相关因子；

S105：构建每种干旱状态下的相关因子和干旱指数SPI时间序列的联合概率分布，推求相关因子给定情况下的干旱指数的条件概率模型；

S106：根据所述干旱状态概率矩阵和所述条件概率模型，计算待预测区域在待预测时段的干旱指数的全状态分布；

S107：基于所述干旱指数的全状态分布，计算得到干旱指数的均值和90％概率区间。

进一步地，步骤S101中，历史月降雨数据包括待预测区域的多个连续时段的月平均降雨量，且以一个月为一个时段。

进一步地，步骤S102中，利用一阶状态空间模型解译干旱指数SPI时间序列，得到状态转移矩阵；具体包括：

S201：采用τ分布概率密度函数拟合历史月降雨数据，得到分布概率曲线；τ分布概率密度函数f(x)如公式(1)所示：

上式中，为x历史月降雨数据；β和γ分别为待求的尺度参数和形状参数，将步骤S101中获取的历史月降雨数据代入公式(1)，采用最大似然法或者矩法求解得到β和γ；

S202：对所述分布概率曲线进行正态标准化处理，得到历史月降雨数据对应的标准化降雨指数SPI；所述历史月降雨数据中的所有时段的月平均降雨量对应的SPI值组成干旱指数SPI时间序列；

S204：利用一阶状态空间模型模拟干旱指数SPI时间序列，获得干旱指数SPI时间序列的状态转移矩阵P(S_t＝j|S_t-1＝i)；S_i为隐状态过程，S_i＝z表示第i个时段处于第z个干旱状态；i＝1,2,…,n，z＝1,2,…,m，n为历史月降雨数据中的时段总数，m为干旱状态总数，z表示第z个干旱状态，m>0，n>0；状态转移矩阵为P(S_t＝j|S_t-1＝i)，表示从t-1时段状态i转移到t时段状态j的概率。

进一步地，步骤S103中，根据所述状态转移矩阵，计算待预测区域在待预测时段处于各干旱状态的概率，得到T时段向T+1时段的干旱状态概率矩阵ξ_T，T时段为待预测时段的前一时段，T+1时段为待预测时段；包括：

采用公式(2)计算待预测区域在待预测时段处于各干旱状态的概率，即待预测时段处于各干旱状态的概率组成的干旱状态概率矩阵ξ_T：

ξ_T＝α_TΓ/α_T1′_m (2)

上式中，ξ_T为m维的向量，1′_m是单位向量；α_T是前向概率向量。

进一步地，步骤S104中，对每一类干旱状态对应的多个SPI值组成的集合，采用相关性分析方法，提取待预测时段前一时段与其具有相关性的影响因子，作为其相关因子R；相关因子R包括：前期SPI、气温、太阳辐射、植被状态和遥相关因子。

进一步地，步骤S105中，采用高斯过程回归模型构建第z类干旱状态下的相关因子R和干旱指数SPI的条件概率模型f_z(SPI_*|R_*,R,SPI)。

进一步地，步骤S106中，待预测区域在待预测时段T的干旱指数的全状态分布P(O_T+1＝SPI_*|O_1:T)的计算公式如公式(3)：

上式中，T+1表示待预测时段，T表示待预测时段的前一时段；ξ_Tz表示干旱状态概率矩阵ξ_T中的第z个元素，即时段T处于第z类干旱状态的概率值。

本发明提供的技术方案带来的有益效果是：本发明所提出的技术方案与现有技术相比，能够取得下列有益效果：

(1)采用状态空间模型分析干旱指数(SPI)的隐状态(Z,状态)并预测未来干旱处于各个状态的概率。状态空间模型揭示了SPI指数序列下的重旱，轻旱等的转移机制和概率。与确定性的干旱分类相比，状态空间方法合理，有效，能反映干旱不确定性本质；

(2)利用联合分布对各个状态分别分析相关因子并建立条件概率预测模型，并将其作为干旱预测的子单元，相比SPI的直接预测方法，这种方法通过多重预测和分类别引入影响因子能提高干旱预测的精度。

(3)通过干旱状态权重结合多个SPI条件概率分布，构造了SPI的全概率公式，与其他方法相比，此过程考虑更加全面和精细化，提供的SPI概率分布结果能提供均值，概率区间，分位值等更加有用的信息，支持防灾减灾决策。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明实施例中一种基于状态空间和联合分布的干旱多级预测方法的流程图；

图2是本发明实施例中基于状态空间和联合分布的干旱多级预测方法的详细流程图；

图3(a)和图3(b)是本发明实施例中SPI与影响因子的联合分布示意图；

图4是本发明实施例中提供的预测结果图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

本发明的实施例提供了一种基于状态空间和联合分布的干旱多级预测方法。

请参考图1，图1是本发明实施例中一种基于状态空间和联合分布的干旱多级预测方法的流程图，具体包括如下步骤：

S101：提取待预测区域的历史月降雨数据，并计算历史月降雨数据对应的标准化降雨指数(SPI)，形成干旱指数SPI时间序列；本发明实施例中，采用标准化降雨指数(SPI)作为干旱指数；

S102：利用一阶状态空间模型解译干旱指数SPI时间序列，得到状态转移矩阵；

S103：根据所述状态转移矩阵，计算待预测区域在待预测时段处于各干旱状态的概率，得到待预测时段的干旱状态概率矩阵；

S104：对每一种干旱状态对应的SPI值(每一个干旱状态对应着多个历史时段，每个时段的月平均降雨量均对应着一个SPI值，因此组成一个SPI值集合)，提取待预测时段前一时段的相关因子，得到各干旱状态下的相关因子；

S105：构建每种干旱状态下的相关因子和干旱指数SPI时间序列的联合概率分布，推求相关因子给定情况下的干旱指数的条件概率模型；

S106：根据所述干旱状态概率矩阵和所述条件概率模型，计算待预测区域在待预测时段的干旱指数的全状态分布(SPI的概率分布函数)；

S107：基于所述干旱指数的全状态分布，计算得到干旱指数的均值和90％概率区间。

步骤S101中，从国家气象科学数据中心(http://data.cma.cn/)提取待预测区域的历史月降雨数据；待预测区域的历史月降雨数据包括待预测区域的多个连续时段的月平均降雨量(以1个月为1个时段)，例如2018年1月至2020年1月的每个月的月平均降雨量。

步骤S102中，利用一阶状态空间模型解译干旱指数SPI时间序列，得到状态转移矩阵；具体包括：

S201：采用τ分布概率密度函数拟合历史月降雨数据，得到分布概率曲线；τ分布概率密度函数f(x)如公式(1)所示：

上式中，为x历史月降雨数据；β和γ分别为待求的尺度参数和形状参数，将步骤S101中获取的历史月降雨数据代入公式(1)，采用最大似然法或者矩法求解得到β和γ；

u为积分变量；

S202：对所述分布概率曲线进行正态标准化处理，得到历史月降雨数据对应的标准化降雨指数(SPI)；具体为：将所述历史月降雨数据中的各月平均降雨量代入公式(1)，计算得到对应的SPI值；所述历史月降雨数据中的所有时段的月平均降雨量对应的SPI值组成干旱指数SPI时间序列；标准化降雨指数(SPI)的具体计算过程可参考：气象干旱等级GB/T20481—2017[S].北京:中国标准出版社，2017。

S204：利用一阶状态空间模型模拟干旱指数SPI时间序列，获得干旱指数SPI时间序列的状态转移矩阵P(S_t＝j|S_t-1＝i)；S_i为隐状态过程，S_i＝z表示第i个时段处于第z个干旱状态；i＝1,2,…,n，z＝1,2,…,m，n为历史月降雨数据中的时段总数，m为干旱状态总数(根据气象干旱国家标准预先设定多个干旱状态)，z表示第z个干旱状态，m>0，n>0；状态转移矩阵为P(S_t＝j|S_t-1＝i)，表示从t-1时段状态i转移到t时段状态j的概率。

步骤S204中，利用一阶状态空间模型模拟干旱指数SPI时间序列，获得干旱指数SPI时间序列的状态转移矩阵P(S_t＝j|S_t-1＝i)；具体包括：

一阶状态空间模型具有属性P(S_t|S_1:t-1)＝P(S_t|S_t-1)，P(O_t|O_1:t-1,S_1:t)＝P(O_t|S_t)，表示序列当前状态只与上一状态有关，当前的观测值(SPI值)只与当前状态有关；

在本发明实施例中，O_t代表时段t(某一个月)的观测值(SPI值)(即第t个时段的SPI值)，S_t代表t第t个时段的干旱状态(根据SPI值(O_t)，结合气象干旱国家标准来确定)，如重旱、轻旱、无旱、轻涝、重涝；一阶状态空间模型L(S_1:T,O_1:T)表达式如公式(2)：

上式中，P(O_t|S_t)代表状态S_t＝z下的观测值激发概率，z＝{1,L,m}，t＝{1,L,T}，m为状态总数，T为历史月降雨数据中的时段总数；所述激发概率为高斯分布：

μ_z,

为状态S_t＝z下的高斯分布均值和方差；P(S_t＝j|S_t-1＝i)是从干旱状态i到干旱状态j的转移概率q_ij；q_ij(i＝1,2,…,m；j＝1,2,…,m)构成(m,m)维度的状态转移矩阵Γ；为了使Γ的元素之和为1，Γ中元素表示如公式(3)：

上式中，Γ_ij为状态转移矩阵Γ中的第i行第j列的元素，即为q_ij；

由上可知，一阶状态模型参数包括：状态总数m、矩阵β和各干旱状态下的高斯分布均值和方差μ_z,

将干旱指数SPI时间序列中的各SPI值代入上述公式(2)、(3)，采用最大似然法求解一阶状态模型参数，得到状态转移矩阵Γ。一阶状态空间模型的具体求解过程可参考：Pender D,Patidar S,Pender G,et al.Stochastic simulation of daily streamflowsequences using a hidden Markov model[J].Hydrology Research,2015:nh2015114.

步骤S103中，根据所述状态转移矩阵，计算待预测区域在待预测时段处于各干旱状态的概率，得到干旱状态概率矩阵ξ_T；包括：

采用公式(4)计算待预测区域在待预测时段处于各干旱状态的概率，即待预测时段处于各干旱状态的概率组成的干旱状态概率矩阵ξ_T：

ξ_T＝α_TΓ/α_T1′_m (4)

上式中，ξ_T为m维的向量，1′_m是单位向量；α_T是前向概率向量，求解参考一阶隐马尔科夫状态空间模型，如公式(5)：

α_T＝δP(O₁)ΓP(O₂)LΓP(O_T) (5)

上式中，δ是一阶状态空间模型初始分布，P(O_t)(t＝1,2,…,T)是概率矩阵(观测值激发概率)，元素为P(O_t＝o_t|S_t＝z)，O_t代表变量，o_t代表O_t的某个具体取值。

步骤S104中，采用相关性分析方法，对每一种干旱状态，继续挖掘其历史数据(包括前期SPI数据，气温数据，太阳辐射数据，植被状态数据和遥相关因子数据等)，分析与该类干旱状态的SPI值具有较高相关性的前期影响因子，作为其相关因子；对属于第z类干旱状态的SPI值，提取待预测时段前一时段的相关因子R；相关因子R有多个，包括：前期SPI、气温、太阳辐射、植被状态和遥相关因子等。

步骤S105中，构建每种干旱状态下的相关因子和干旱指数SPI时间序列的联合概率分布，推求相关因子给定情况下的干旱指数(SPI)的条件概率模型；

根据高斯过程回归模型，给定条件因子x，预测对象y得条件分布如公式(6)：

上式中

σ是高斯噪声标准差，

是协方差矩阵，p为协方差矩阵的维度，x为什么维度，p就为什么维度；x_*，y_*代表x,y某个具体值；

采用高斯过程回归模型构建干旱状态z(即第z类干旱状态)下的相关因子R和干旱指数(SPI)的条件概率模型f_z(SPI_*|R_*,R,SPI)；具体可参考：Zhu S,Luo X,Xu Z,etal.Seasonal streamflow forecasts using mixture-kernel GPR and advancedmethods of input variable selection[J].Nordic Hydrology,2019,50(1-2):200-214.

步骤S106中，待预测区域在待预测时段T+1的干旱指数的全状态分布P(O_T+1＝SPI_*|O_1:T)的计算公式如公式(7)：

上式中，T+1表示待预测时段，T表示待预测时段的前一时段；ξ_Tz表示干旱状态概率矩阵ξ_T中的第z个元素，即时段T处于第z类干旱状态的概率值；f_z(SPI_*|R_*,R,SPI)表示干旱状态z下的相关因子R和干旱指数(SPI)的条件概率模型。

本发明实施例以湖北省应城市干旱预测为例，该区域位于湖北省中部偏东，介于东经113°19′-113°45′，北纬30°43′-31°08′之间。市境内有大富水、汉北河、府河三条河流，水资源分布东南多西北少，年内分配不均。全市水资源多年人均日占有量为0.24立方米/人·日，亩均日占有量0.26立方米/亩·日，均低于湖北省水平。从水资源的时空和总量分布看，应城是一个水资源相对贫乏的地区，基本是三年一大灾，年年有小灾。分析表明应城出现春旱和秋旱的机率较大，有时还出现秋冬春、春夏秋三季连旱的情况。

采用本发明实施例提供的基于状态空间和联合分布的干旱多级预测方法的流程(如图2)，具体包括如下步骤：

(1)根据实例1961～2015年月尺度的降水，提取657个SPI 3个月尺度的干旱指数。

(2)将SPI数据分成训练期和验证期，1-527为训练期，528-657为验证期，对训练期SPI指数进行干旱状态空间分析，求解状态转移矩阵。分析结果表面SPI序列分为4个状态，均值分别为-1.593，-0.448，0.448和1.537，方差均为0.25。

(3)依据确定的状态矩阵，以527个位置SPI状态预测第528位置的干旱状态，以第528位置SPI状态预测529位置的干旱状态，滚动进行。得到各个状态的概率数据。

(4)选取每个状态下的前一个月SPI指数作为影响因子(相关因子)，构建SPI与影响因子的联合分布，如图3(a)和图3(b)所示。

(5)以预测528位置SPI指数为例，提取前一个月SPI数据,从四个状态的联合分布函数重推求给定前一个月SPI的当前SPI条件分布，因此SPI的条件分布有四组。

(6)以预测的528位置的状态概率作为权重，结合四个状态的条件分布，得到SPI指数的最后预测分布。

(7)通过以上步骤，可预测检验期所有的SPI指数分布。SPI预测情况如附图4.

预测结果不是单一的点预报，附加90％的置信区间能比点预报提供更可多信息用于决策。

将基于状态空间和联合分布的干旱多级预测方法和一些常用干旱指数预测方法(自回归滑动平均和人工神经网络)的预测性能做比较，如表1.本发明实施例提出的方法确定性系数最高(R²)，平均绝对误差(MAE)、均方根误差最小(RMSE)，基于均方根(SS_MSE)的预测能力得分最高，证明了本方法能有效提高干旱指数预测精度。

通过以上实例的验证可以认为，本方法能够精确预测干旱事件及其概率，可以推广应用。

本发明的有益效果是：本发明所提出的技术方案与现有技术相比，能够取得下列有益效果：

(1)采用状态空间模型分析干旱指数(SPI)的隐状态(Z,状态)并预测未来干旱处于各个状态的概率。状态空间模型揭示了SPI指数序列下的重旱，轻旱等的转移机制和概率。与确定性的干旱分类相比，状态空间方法合理，有效，能反映干旱不确定性本质；

(2)利用联合分布对各个状态分别分析相关因子并建立条件概率预测模型，并将其作为干旱预测的子单元，相比SPI的直接预测方法，这种方法通过多重预测和分类别引入影响因子能提高干旱预测的精度。

(3)通过干旱状态权重结合多个SPI条件概率分布，构造了SPI的全概率公式，与其他方法相比，此过程考虑更加全面和精细化，提供的SPI概率分布结果能提供均值，概率区间，分位值等更加有用的信息，支持防灾减灾决策。

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

Claims (3)

1.一种基于状态空间和联合分布的干旱多级预测方法，其特征在于：包括以下步骤：

S101：提取待预测区域的历史月降雨数据，并计算历史月降雨数据对应的标准化降雨指数SPI，采用标准化降雨指数SPI作为干旱指数，形成干旱指数SPI时间序列；

S102：利用一阶状态空间模型解译干旱指数SPI时间序列，得到状态转移矩阵；

步骤S102中，利用一阶状态空间模型解译干旱指数SPI时间序列，得到状态转移矩阵；具体包括：

S201：采用τ分布概率密度函数拟合历史月降雨数据，得到分布概率曲线；τ分布概率密度函数f(x)如公式(1)所示：

上式中，为x历史月降雨数据；β和γ分别为待求的尺度参数和形状参数，将步骤S101中获取的历史月降雨数据代入公式(1)，采用最大似然法或者矩法求解得到β和γ；

u为积分变量；

S202：对所述分布概率曲线进行正态标准化处理，得到历史月降雨数据对应的标准化降雨指数SPI；所述历史月降雨数据中的所有时段的月平均降雨量对应的SPI值组成干旱指数SPI时间序列；

S204：利用一阶状态空间模型模拟干旱指数SPI时间序列，获得干旱指数SPI时间序列的状态转移矩阵P(S_t＝j|S_t-1＝i)；S_t为隐状态过程，S_t＝z表示第t个时段处于第z个干旱状态；t＝1,2,…,n，z＝1,2,…,m，n为历史月降雨数据中的时段总数，m为干旱状态总数，z表示第z个干旱状态，m>0，n>0；状态转移矩阵为P(S_t＝j|S_t-1＝i)，表示从t-1时段状态i转移到t时段状态j的概率；

S103：根据所述状态转移矩阵，计算待预测区域在待预测时段处于各干旱状态的概率，得到待预测时段的干旱状态概率矩阵；

步骤S103中，根据所述状态转移矩阵，计算待预测区域在待预测时段处于各干旱状态的概率，得到T时段向T+1时段的干旱状态概率矩阵ξ_T，T时段为待预测时段的前一时段，T+1时段为待预测时段；包括：

采用公式(2)计算待预测区域在待预测时段处于各干旱状态的概率，即待预测时段处于各干旱状态的概率组成的干旱状态概率矩阵ξ_T：

ξ_T＝α_TΓ/α_T1′_m (2)

上式中，ξ_T为m维的向量，1′_m是单位向量；α_T是前向概率向量，Γ是状态转移矩阵；

S104：对每一种干旱状态对应的SPI值，提取待预测时段前一时段的相关因子，得到各干旱状态下的相关因子；

S105：构建每种干旱状态下的相关因子和干旱指数SPI时间序列的联合概率分布，推求相关因子给定情况下的干旱指数的条件概率模型；

步骤S105中，采用高斯过程回归模型构建第z类干旱状态下的相关因子R和干旱指数SPI的条件概率模型f_z(SPI_*|R_*,R,SPI)；

S106：根据所述干旱状态概率矩阵和所述条件概率模型，计算待预测区域在待预测时段的干旱指数的全状态分布；

步骤S106中，待预测区域在待预测时段的干旱指数的全状态分布P(O_T+1＝SPI_*|O_1:T)的计算公式如公式(3)：

上式中，T+1表示待预测时段，T表示待预测时段的前一时段；ξ_Tz表示干旱状态概率矩阵ξ_T中的第z个元素，即时段T处于第z类干旱状态的概率值；

S107：基于所述干旱指数的全状态分布，计算得到干旱指数的均值和90％概率区间。

2.如权利要求1所述的一种基于状态空间和联合分布的干旱多级预测方法，其特征在于：步骤S101中，历史月降雨数据包括待预测区域的多个连续时段的月平均降雨量，且以一个月为一个时段。

3.如权利要求1所述的一种基于状态空间和联合分布的干旱多级预测方法，其特征在于：步骤S104中，对每一类干旱状态对应的多个SPI值组成的集合，采用相关性分析方法，提取待预测时段前一时段与其具有相关性的影响因子，作为其相关因子R；相关因子R包括：前期SPI值、气温、太阳辐射、植被状态和遥相关因子。

Images