基于聚类分析的模糊马尔科夫链在降雨量预测中的应用

宋 帆，杨晓华，武翡翡，刘 童

(北京师范大学环境学院 水环境模拟国家重点实验室，北京 100875)

降雨量是流域或地区水资源循环系统的重要输入项目，是最主要的水资源补给途径，在农业生产、工业建设、水利防洪等方面起着重要的作用。由于气象条件具有多变性以及复杂性，导致降雨过程存在较大的随机性[1]，如何较为准确地对降雨量进行预测也成了一个热点研究问题。目前较为常用的降雨量预测方法有灰色预测[2]、马尔科夫链[3]、人工神经网络[4]等。

马尔可夫过程是研究事物状态及状态转移规律的理论。它通过确定系统不同状态的初始概率及状态之间的转移概率关系, 来确定状态的变化趋势, 从而进行预测[5]。常用的马尔科夫链有基于绝对分布的马尔科夫链、叠加马尔科夫链、加权马尔科夫链等，其中加权马尔科夫链是效果最好、应用最广泛的[6]。本次研究在加权马尔科夫链的基础上进行了改进：首先利用欧氏距离聚类代替常用的均值分类；其次引入隶属度函数来确定初始状态向量；最后对预测结果进行分类改进。从而建立了基于聚类分析的模糊马尔科夫链降雨量预测模型，并利用中国气象站记录的1951-2013年降雨量数据对模型进行验证。

1 模型原理及步骤

1.1 马尔科夫预测原理

马尔科夫预测是根据马尔科夫理论与方法来研究未来发展趋势的一种动态分析方法，既适用于时间序列，又适用于空间序列。其基本特征是无后效性，即系统在未来某一时刻的状态只与现在的状态有关，而与过去的状态无关。因此马尔科夫链在模拟预测随机波动性较大的问题上有较好的效果[7]。

1.2 聚类-模糊马尔可夫链模型流程

本次研究所采用的聚类-模糊马尔可夫链模型如图1所示，首先通过聚类分析将降雨量序列进行分类，随后计算转移概率矩阵，并进行马氏性检验。为了确立初始状态向量，引入隶属度的概念，对参照样本初始状态进行测算。随后通过自相关系数的计算来确定各阶权重。三者结合即可得到待测样本降雨量的初步预测值，再经过结果改进即可得到最终预测值。

图1 模型流程Fig.1 The process of model

1.3 模型预测步骤

(1)历史数据聚类分析。用聚类分析代替均值分类法对历史降雨量数据进行聚类分析，将历史降雨量数据分为5～6组。具体操作在SPSS软件中进行，使用欧氏距离公式测算各年数据间的相似性，运用组间联结法进行聚类[8]。

欧式距离的计算公式为:

(1)

(2)求转移概率矩阵。根据聚类结果可将降雨量序列分为N个状态。即状态空间E={1,2,…,N}。若本次研究的降雨量序列{Sn}=S1,S2,…,Sn为非负参数，且第n+1个参数Sn+1的状态与其他参数无关，只是Sn的一个概率函数P，那么可得到此时的一步转移概率矩阵：

P(Sn+1=j|Sn=i)=Pij(1) (i,j∈E)

(2)

同理可得步长为k时的转移概率矩阵如下[9]：

(3)

其中Pij=fij/fi，fij表示由i状态到j状态对应的频数，fi表示样本总数。在步长的选择上，如果选择的步长过大结果会明显失真，为了保证预测结果的准确性，本次探究步长k最大取5,即k=1,2,3,4,5[10]。

(3)“马氏性”检验。检验历史降雨量序列是否具有马氏性, 是能否应用马尔可夫链模型进行预测的前提，通常用χ2统计量来检验。假设将一步转移频数矩阵的第j列之和与各行各列总和的比值称为“边际概率”，即：

(4)

统计量：

(5)

(4)求各参照样本状态向量。模糊马尔科夫链模型与普通马尔科夫链模型相比，不仅仅需要转移概率矩阵，并且需要各参照样本的初始状态向量。在此引入隶属度的概念，利用隶属度公式测算各参照样本分属于不同状态的隶属度，从而确定一个最准确的初始状态[6]。假设目标年份降雨量为A，聚类分析的各等级降雨量下临界值分别为A0,A1,A2,…,Ax,那么各参照样本降雨量对各状态的隶属度计算如下：

(6)

(5)求待测样本各阶状态向量。将参照样本年份的状态向量与各步转移概率矩阵相乘即可得到待测样本年份的各阶状态向量。例如，本次研究中2011年降雨量数据为待测数据，将2010年状态向量与一步转移概率相乘即可得到2011年一阶状态向量，同理由2009年状态向量以及二步转移概率矩阵可得到2011年二阶状态向量。

(6)计算各阶自相关系数以及权重。为了正确反映步长对马尔科夫链预测值的重要性，采用各阶自相关系数来计算权重，具体公式如下：

(7)

(8)

式中：m为步长k的最大值，即5。

(8)年降雨量预测值。根据状态概率可得到当年降雨量预测值Xn+1。具体公式如下：

Sn+1=HD/i

(9)

(10)

式中：H为状态特征值；Pi为各分类状态的预测概率；D为各预测状态对应的下临界值。

(9)结果改进。上述步骤在应用时结果并不十分理想，因此本次研究在实际计算的基础上对结果进行了分类改进。具体公式如下：

(11)

式中：T是马尔科夫预测的待测样本所处的状态等级；H为状态特征值。

2 数据来源

本次研究所用数据全部来自中国气象数据网中的《中国地面气候年值数据集》。选用了16个气象数据站的1951-2013年降雨量数据，将1951-2010年数据作为参照样本，2011-2013年为待测样本。采用基于聚类分析的模糊马尔科夫链来进行预测。

3 实例研究

3.1 研究步骤

以江苏省为例进行模型验证如下：

(1)用SPSS软件将降雨量数据进行聚类分析，共分为枯水年、偏枯年、正常年、偏丰年和丰水年5类，结果如表1所示(1951-2010年数据按顺序标为1～60)。

(2)根据公式(3)对分类结果进行统计计算，得到步长1～5时的转移概率矩阵如下，它决定了降雨量指标状态转移过程的概率法则。

表1 降雨量序列分类Tab.1 Classification of rainfall sequence

(4)根据表1中分类以及公式(6)可得2006-2010年降雨量状态向量，再分别与相应的概率转移矩阵相乘可得2011年归一化状态向量如下：

(5)根据公式(7)、(8)可求得各阶自相关系数以及相应的权重，如表2所示。

表2 各阶自相关系数以及权重Tab.2 Autocorrelation coefficients and weights

(7)根据改进后的公式(11)对2011年降雨量数值进行预测，同理重复步骤(2)～(6)对2012和2013年降雨量进行预测，结果见表3。

表3 南京市降雨量预测结果Tab.3 Rainfall prediction results of Nanjing

3.2 结果讨论

3.2.1 精度分析

使用基于结果改进的聚类—模糊马尔科夫链对 2011-2013年3年降水量进行预测，结果显示3年降雨量的相对误差全部都在10%之内，平均误差为8.37%。说明建立的马尔科夫链模型可用于降雨量的预测，并且结果改进后具有较高精度。

3.2.2 遍历性与平稳分布

由于降雨量的5个状态枯水年、偏枯年、正常年、偏丰年和丰水年是互通的，并且都没有周期性，因此全部状态构成的空间是一个闭集，那么此马尔科夫链是不可约的。因此本次研究的马尔科夫链是一个正常返链，具有遍历性，当其达到平稳状态时的分布为平稳分布，也就是其极限分布。此时极限分布{πj,j∈E}以及各状态出现的周期Ti公式如下：

(12)

(13)

(14)

本次研究中步长为3时的模型拟合性最好，因此采用三步转移概率矩阵求平稳分布，结果如表4所示。

表4 平稳分布与状态周期Tab.4 Stationary distribution and period of state

由表4可知，根据南京市1951-2013年历史降雨量数据推算，偏枯年降雨状态出现的几率最大，平均2.47年出现一次，而丰水年出现的几率最小，平均每隔20.62年才出现一次。这个数据与各状态实际占比几乎相同，说明该周期预测合理，可为城市用水规划提供依据。

3.3 拓展预测

为了验证模型的正确性以及适用范围，在数据可获得且完整的基础上，随机取了江苏省另外两个城市降雨测站以及其他省13个降雨量测站数据进行进一步验证。各城市2011-2013年降雨量预测结果以及相对误差见表5。

表5 拓展预测结果Tab.5 The results of expanded prediction

由表5以及南京市预测数据可知，16个城市共48个预测数据中，绝对误差超过15%的年份有15个，占31.25%；绝对误差低于10%的年份有27个，占56.25%。48个预测数据总体相对误差为12.4%，总的来说预测结果相对较好。但由于气象条件的复杂性和偶然性，有一些年份会出现较大误差，这是不可避免的。从各城市平均误差来看，超过15%的城市有5个，占31.25%，与高误差年份占比一致。但误差低于10%的城市只有7个，占43.75%，可见个别极端年份对整体误差的影响较大，而这可作为今后降雨量预测研究的重要方向。

4 结 语

本文以降雨量预测为研究目标，建立了基于聚类分析的模糊马尔科夫链降雨量预测模型，并以南京市降雨量数据为实例进行分析，结果显示2011-2013年降雨量预测平均误差为8.37%，预测精度较高。随后将结果改进后的模型推广到全国范围内的其他15个城市，结果显示误差低于10%的城市占43.75%，而误差低于10%的年份占56.25%。表明模型在预测降雨量方面效果较好，可为水资源的规划提供一定依据。但模型在异常值的预测上能力较差，整体预测精度也因此下滑，如何进一步改善极端值年份的影响将是今后的研究重点。