CN102183802B - 一种基于卡尔曼滤波及演化建模的短期气候预测方法 - Google Patents

Abstract

本发明公开了一种基于卡尔曼滤波与演化建模结合的气候预测方法，本方法首先通过卡尔曼滤波建立关于预测因子的线性模型，在此基础上利用非线性常微分方程数学模型去模拟逼近卡尔曼滤波中的误差序列，并进行误差预测。演化算法是用计算机模拟大自然的演化过程，特别是模拟生物演化过程来求解复杂问题的一类计算方法，具有自适应、自组织、自学习以及内在并行性等智能特性。将两种算法相结合，可以比单纯线性模型更好的模拟气候的自然特性，提高气候预测精度。通过本发明方法，可以对短期的日照时数、温度、降水量进行预测，从而提供短期气候的先验知识。

Description

一种基于卡尔曼滤波及演化建模的短期气候预测方法

技术领域

本发明属于一种短期气候预测的方法，尤其涉及一种基于卡尔曼滤波及演化建模的短期气候预测方法。该方法可用于解决如日照时数、温度等随时间波动较大的短期气候的预测，也同样适用于长期气候的预测。

背景技术

随着我国国民经济的迅速发展，国家越来越需要更加准确的短期气候预测，尤其是在国家计划、农业、水利和防灾减灾等部门。短期气候预测是国际大气科学和地球科学领域的前沿课题，也是极其困难的跨学科难题。目前短期气候预测主要面临三大方面难题，首先是气候的变化存在很多不确定因素，短期气候数据时间序列波动性太大，无法确定其准确模型。第二是对短期气候变化的预测方程组极其复杂，而且从数学和物理上如何正确表征各因素之间的相互作用，还没有完全解决。因而面临着预测理论的困难，所以短期气候预测不能沿用长期气候预测的原理和方法，必须面对整个复杂的气候***及其变化研制新的理论和方法，其技术水平也受到其它相关学科发展的制约和限制。第三资料和数据的不足，为了揭示短期气候预测的宏观和微观规律，需要大量这方面的历史数据。

气候要素观测数据时间序列可以被视作是一个精确完美的气候数值预测模式在完备条件下的解，是各气候因子以及自身相互间非线性作用的综合演化表现形式。随着气象观测方法的多元化和大型计算机的应用，特别是气象学家对气候问题的拓展和深入，短期气候数值预测取得了一定进展，但大部分对于短期气候的预测仍然是对观测数据时间序列进行线性建模分析，而一些气候预测的非线性模型又必须要求模型非线性函数已知，由人为主观估计。因此有必要找到一种数学方法客观地模拟具有非线性特征的气象学要素时间序列。

基于卡尔曼滤波及演化建模的短期气候预测方法是在卡尔曼滤波建立的线性模型与演化建模的非线性模型相互作用的结果。卡尔曼滤波根据不断更新的数据与预测值之间的误差调整回归系数，达到实时跟踪预测的目的。而演化建模方法是以遗传算法为基础，模拟生物演化过程来求解复杂问题的一类计算方法，其具有自适应、自组织、自学习以及内在的并行性等智能特性。因此本发明将线性模型与非线性模型相叠加，在目前不能确定气候数值模型即非线性偏微分方程组中气候因子非线性相互作用物理机制的显示数学表达式的情况下，用非线性常微分方程数学模型区逼近气象要素历史观测时间序列，并进而作短期气候预测。同时他们都具有利用实时数据不断调整模型参数或表达式的能力，使预测精度更高。目前，基于卡尔曼滤波及演化建模的短期气候预测未见有文献公开报道。

发明内容

本发明的目的在于克服上述现有技术不足，提供一种基于卡尔曼滤波及演化建模的短期气候预测方法，利用卡尔曼滤波与演化建模相结合的方式，用卡尔曼滤波递推***描述预测因子之间的线性关系，根据其预测误差不断调整回归系数；同时在线性模型的基础上叠加由演化建模描述的针对卡尔曼滤波误差所建立的非线性模型，从而优化预测模型，提高预测精度。

为了达到以上目的，本发明采取如下技术方案予以实现：

包括下列步骤：

A)收集气候预测模型中预测因子的历史数据，包括近年来的日照时数、温度、相对湿度和降雨量；

B)确定预测因子，整理与预测因子相关的历史数据，将历史数据分段处理，将每段数据的综合作用结果作为一个数据项目，同时将处理后的数据分为初值计算样本和检测样本两部分；

C)根据上述样本计算卡尔曼滤波初值：

Y_t＝X_tβ_t+e_t

上式为卡尔曼滤波测量方程，首先利用多元线性回归计算回归系数β₀，

再计算出动态噪声方差阵W和测量噪声方差阵V；

D)根据上述样本建立卡尔曼递推模型：利用步骤C)中计算的初值及卡尔曼递推式对检测样本的数据进行预测，并与检测样本进行比较，得到预测误差时间序列；

E)对上述预测误差进行预处理：

$x\left(t\right)=\stackrel{\‾}{x}\left(t\right)+\tilde{x}\left(t\right)$

$\tilde{x}\left(t\right)=x\left(t\right)-\stackrel{\‾}{x}\left(t\right)$

将预测误差分为光滑部分和粗糙部分，并假设光滑部分是由宏观气候因素控制，粗糙部分是由微观气候因素控制，针对分解后的数据分别采用以下步骤F)和步骤G)的数学模型来描述；

F)针对步骤E)中的光滑部分时间序列项，用染色体来表示问题的可行解，包括可行解的形式和参数，然后利用遗传操作在解空间中进行搜索，最后根据拟合表达式与实际值之间的偏差平方和进行优胜劣汰，求得最终的解，即拟合光滑部分时间序列最好的非线性表达式；

G)针对步骤E)中的粗糙部分时间序列项，用自然基小波对其进行分析，从这些时间序列中找出其受气候因子影响下的周期波动规律，实现粗糙部分时间序列的拟合和预测功能；

H)将步骤D)中实现的卡尔曼滤波以及步骤F)和步骤G)中分别针对卡尔曼滤波误差而进行的演化建模和自然基小波分别对各自的气候数据时间序列模型做下一步预测，并将预测结果相叠加，从而实现卡尔曼滤波及演化建模的短期气候预测。

将步骤D)中，卡尔曼滤波递推模型对检测样本进行预测的同时，其卡尔曼滤波回归系数根据预测误差做出实时调整。

步骤E)中，可根据具体情况调整光滑参数l的值，l可用于调节时间序列

的光滑程度，l越大，

越光滑。

步骤F)中，染色体表示的微分方程模型为一个树表达式，其微分方程解由Runge-Kutta方法求出。

步骤G)中，针对时间序列的粗糙部分建模，采用的是多尺度自然分形模型。

步骤H)中，短期气候预测模型为线性预测模型与非线性预测模型相结合的方式。

本发明提出基于线性卡尔曼滤波与非线性常微分方程相结合的气候预测模型：用卡尔曼滤波建立气候预测中各预测因子之间的线性关系，同时用非线性微分方程拟合预测因子的非线性表现，其中非线性微分方程采用演化建模方法求得，卡尔曼滤波参数和微分方程都随时间移动做即时调整，能够更好反映预测量随时间的变化规律，提高预测精度。

附图说明

图1是本发明的建模流程图。

图2是常微分方程树结构图。

图3为本发明的演化建模流程图。

图4是本发明的基于自然基小波算法流程图。

下面结合附图对本发明的技术方案作进一步详细说明。

具体实施方式

参照图1、图2、图3、图4所示，一种基于卡尔曼滤波及演化建模的短期气候预测方法，包括下列步骤：

A)收集气候预测模型中预测因子的历史数据，主要包括近年来的日照时数、温度、相对湿度和降雨量等；

B)确定预测因子，整理与预测因子相关的历史数据。将历史数据分段处理，将每段数据的综合作用结果作为一个数据项目。同时将处理后的数据分为初值计算样本和检测样本两部分；

C)根据上述样本计算卡尔曼滤波初值：

Y_t＝X_tβ_t+e_t

上式为卡尔曼滤波测量方程，首先利用多元线性回归计算回归系数β₀，

再计算出动态噪声方差阵W和测量噪声方差阵V。

D)根据上述样本建立卡尔曼递推模型：利用步骤C)中计算的初值及卡尔曼递推式对检测样本的数据进行预测，并与检测样本进行比较，得到预测误差时间序列；

E)对上述预测误差进行预处理：

$x\left(t\right)=\stackrel{\‾}{x}\left(t\right)+\tilde{x}\left(t\right)$

$\tilde{x}\left(t\right)=x\left(t\right)-\stackrel{\‾}{x}\left(t\right)$

将预测误差分为光滑部分和粗糙部分，并假设光滑部分是由宏观因素控制，粗糙部分是由微观因素控制，针对分解后的数据分别采用以下步骤F)和步骤G)的数学模型来描述；

F)针对步骤E)中的光滑部分时间序列项，用染色体来表示问题的可行解，包括可行解的形式和参数，然后利用遗传操作在解空间中进行搜索，最后根据拟合表达式与实际值之间的偏差平方和进行优胜劣汰，求得最终的解，即拟合光滑部分时间序列最好的非线性表达式；

G)针对步骤E)中的粗糙部分时间序列项，用自然基小波对其进行分析，从这些时间序列中找出其微观规律，实现粗糙部分时间序列的拟合和预测功能；

H)将步骤D)中实现的卡尔曼滤波以及步骤F)和步骤G)中分别针对卡尔曼滤波误差而进行的演化建模和自然基小波分别对各自的气候数据时间序列模型做下一步预测，并将预测结果相叠加，从而实现卡尔曼滤波及演化建模的短期气候预测。

上述方案中，步骤D)所述的针对气候预测的卡尔曼滤波递推模型为：

${\hat{Y}}_t=X_t{\widehat{\mathrm{\β}}}_{t-1}---\left(1\right)$

R_t＝C_t-1+WY                                    (2)

${\mathrm{\σ}}_t=X_t{R}_t{X}_t^T+V---\left(3\right)$

$A_t=R_t{X}_t^T{\mathrm{\σ}}_t^{-1}---\left(4\right)$

${\widehat{\mathrm{\β}}}_t={\widehat{\mathrm{\β}}}_{t-1}+A(Y_t-{\hat{Y}}_t)---\left(5\right)$

$C_t=R_t-A_t{\mathrm{\σ}}_t{A}_t^T---\left(6\right)$

中(1)式为预测方程，

为预报值，X_t为预测因子，

为回归系数，是β_t-1的估计值。式(2)中的C_t-1为

误差方差阵，W是ε_t的方差矩阵。式(3)中的σ_t是预报误差方差阵，R_t为递推值

的误差方差阵，

是X_t的转置矩阵，V是e_t的方差矩阵。式(4)中A_t是增益矩阵，

是σ_t的逆矩阵。式(5)中Y_t是预报值的实际测量值。式(5)和式(6)根据预报误差递推下一步的

和C_t。

上述方案中，步骤F)所述的针对光滑部分时间序列的常微分方程组串行演化建模算法采用的是两级演化建模思想，即在以GP(genetic programming)技术优化模型结构的同时，嵌套以GA(genetic algorithm)优化模型参数的过程，其基本步骤如下：

1)初始化模型种群

2)采用GP技术优化模型结构

3)对模型进行简洁化和规范化处理

4)采用GA算法对模型的参数进行优化

5)检查是否满足终止条件，如果没有，则跳转到步骤2)继续执行，最后输出种群中匹配度最好的模型。

上述方案中，步骤G)所述的粗糙部分时间序列的自然基小波分析，其表达的主要思想为：针对粗糙部分时间序列X＝(x(t₁)，x(t₂)，...，x(t_i)，...，x(t_n))^T，按每组k个数据依次排序，构成二维数组x₂(1：m_k，1：k)，则波函数

在k个时间间隔的值为：

自然基小波从波长

开始，依次搜索

的自然基小波。

$F=\left(\frac{D^2}{l-1}\right)/\left(\frac{S^2}{n-1}\right)$

其中： $D^2=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{m}_k{({{\tilde{x}}_i}^{\left(w\right)}-\stackrel{\‾}{x})}^2,$ $S^2=\underset{j=1}{\overset{m_k}{\mathrm{\Σ}}}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{(x_2(j,i)-{{\tilde{x}}_i}^{\left(w\right)})}^2$

满足自由度为(l-1，n-l)的F分布，若F≥F_a(α取0.05)，就以第k-1次的残差序列扣除本次自然基小波在时间域的周期展开，得到一个第k次残差时间序列，然后以波长递减直到k＝2终止。

为了验证本发明基于卡尔曼滤波及演化建模的短期气候预测方法的优越性，以下结合实施例对本发明作进一步的详细描述。

本发明是基于卡尔曼滤波及演化建模的短期气候预测方法，本发明建模的建模原理是，首先多元线性回归利用气候历史数据确定卡尔曼滤波初始值，然后结合气候历史数据和初始值建立卡尔曼滤波递推***模型。再利用演化建模和自然基小波对卡尔曼滤波预测值与实际测量值之间的误差进行优化，以克服单纯基于卡尔曼滤波建立的线性模型的缺陷，从而实现卡尔曼滤波和演化建模混合建模方法。以提高短期气候预测精度。

针对上述具体问题，本发明设计的基于卡尔曼滤波及演化建模的短期气候预测方法具体描述如下：

1)数据采集：针对气候预测的大量历史气候数据，主要从气象站获得，气候包括往年每一天的日照时数、平均温度、相对湿度等。

2)确定卡尔曼滤波初值：为了构建卡尔曼滤波预测模型，必须首先确定卡尔曼滤波递推***中的回归系数误差方差阵C₀、动态噪声方差阵W、测量噪声方差阵V的初值。首先根据预测量确定预测因子，根据已知气候数据建立多元线性回归方程，通过线性回归方程解：

Y＝Xβ+ε

β＝A^-1B＝(X^TX)^-1(X^TY)

求出回归系数初值。

C₀的确定：C₀是

的误差方差阵，这里直接取为：C₀＝[0]_m×m

W的确定：W是动态噪声ε_t的方差矩阵，这里假设其为白噪声，W可表示为：

$W=\left[\begin{array}{cccc}{\mathrm{\ω}}_1& \·\·\·& \·\·\·& 0\\ \·\·\·& {\mathrm{\ω}}_2& \·\·\·& \·\·\·\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ 0& \·\·\·& \·\·\·& {\mathrm{\ω}}_m\end{array}\right]$

这里按照多元线性回归计值的办法求得，这里我们选取的是2004年的气候测量数据，这样

V的确定：V是e_t的方差阵，同样假设其为白噪声，这V可表示为：

$V=\left[\begin{array}{cccc}{v}_1& 0& \·\·\·& 0\\ 0& v_2& \·\·\·& 0\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ 0& 0& \·\·\·& v_m\end{array}\right]$

利用预报量Y的n个分量建立回归方程，求出n个残差(q₁，q₂，...，q_n)，于是：

$V=\left[\begin{array}{cccc}\frac{q_1^2}{k-m-1}& 0& \·\·\·& 0\\ 0& \frac{q_2^2}{k-m-1}& \·\·\·& 0\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ 0& 0& \·\·\·& \frac{q_n^2}{k-m-1}\end{array}\right]$

其中k是样本容量，k＞m+1。

3)建立卡尔曼滤波预测模型：卡尔曼滤波预测模型建立在如下递推***的基础之上：

${\hat{Y}}_t=X_t{\widehat{\mathrm{\β}}}_{t-1}---\left(1\right)$

R_t＝C_t-1+WY                                                        (2)

${\mathrm{\σ}}_t=X_t{R}_t{X}_t^T+V---\left(3\right)$

$A_t=R_t{X}_t^T{\mathrm{\σ}}_t^{-1}---\left(4\right)$

${\widehat{\mathrm{\β}}}_t={\widehat{\mathrm{\β}}}_{t-1}+A(Y_t-{\hat{Y}}_t)---\left(5\right)$

$C_t=R_t-A_t{\mathrm{\σ}}_t{A}_t^T---\left(6\right)$

中(1)式为预测方程，

为预报值，X_t为预测因子，

为回归系数，是β_t-1的估计值。式(2)中的C_t-1为

误差方差阵，W是ε_t的方差矩阵。式(3)中的σ_t是预报误差方差阵，R_t为递推值的误差方差阵，

是X_t的转置矩阵，V是e_t的方差矩阵。式(4)中A_t是增益矩阵，

是σ_t的逆矩阵。式(5)中Y_t是预报值的实际测量值。式(5)和式(6)根据预报误差递推下一步的和C_t。

通过卡尔曼滤波递推***进行短期气候预测，并用测量值减去预测值，从而得到关于预测差值的时间序列。

4)演化建模：为进一步提高卡尔曼滤波的预测精度，本发明采用在卡尔曼滤波线性模型的基础上叠加由演化建模建立的非线性模型。演化建模是用计算机模拟大自然的演化过程，用染色体来表示问题的可行解，然后执行遗传操作在解空间中进行搜索，选出优秀个体，这样一代一代演化，最终求得问题的最好可行解。

图3为本发明的演化建模流程图。

本发明采用遗传程序设计的演化算法，用树型结构来表示如下微分方程：

$\frac{{\mathrm{dy}}^{\left(n\right)}\left(t\right)}{\mathrm{dt}}=f(t,y^{\left(1\right)},y^{\left(2\right)},...,y^{\left(n\right)})$

例如将

y⁽⁴⁾(t)＝y⁽³⁾+4y⁽¹⁾*y⁽²⁾-t*e^y

表示成图2形式，对应于树的叶节点的内容区自然数集D＝{R，y₁，y₂，y₃，...，y_n}，树的非叶节点的内容取自运算集O＝{+，-，×，/，exp，sin，cos，ln}，于是模型空间由数集D、运算集O以及数的高度H确定。当树型表达式确定，即可用Runge-Kutta方法计算出解y^*(t)在相继时刻t₂，...，t_m的值矩阵Y^*，从而计算模型拟合误差：‖Y-Y^*‖以评估该模型的好坏，当最佳模型确定后，就可以用t_m的值做初始条件进行预测。

因此演化建模主要步骤为：

(1)确定演化建模的参数(包括种群数目P、每个种群规模n、最大演化代数Maxgeno、最大树深度D等)以及并行算法的一些控制参数(包括种群间的迁移率ρ，迁移代频g等)；

(2)各种群独立地进行演化建模的过程，包括种群的初始化、模型的结构优化和参数优化(执行选择-杂交-变异)和个体的适应度值评估。同时各种群在演化建模的过程中与其他种群进行信息交换：

①每隔g代将当前种群中的m个最好的个体发送到其他的δ个种群(在原种群中保留其副本)。其中m由迁移率ρ和子种群大小n_c决定。

②每隔g代接受由其他种群发送过来的m个个体，并替换当前种群中的m个最差个体。

(3)当种群运行到最大代数Maxgeno时，程序终止，选出最优个体作为最终预测模型。

5)自然分型模型：图4是本发明的基于自然基小波算法流程图，针对粗糙部分时间序列X＝(x(t₁)，x(t₂)，...，x(t_i)，...，x(t_n))^T，按每组k个数据依次排序，构成二维数组x₂(1：m_k，1：k)，则波函数

在k个时间间隔的值为：

自然基小波从波长

开始，依次搜索

的自然基小波，

$F=\left(\frac{D^2}{l-1}\right)/\left(\frac{S^2}{n-1}\right)$

其中： $D^2=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{m}_k{({{\tilde{x}}_i}^{\left(w\right)}-\stackrel{\‾}{x})}^2,$ $S^2=\underset{j=1}{\overset{m_k}{\mathrm{\Σ}}}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{(x_2(j,i)-{{\tilde{x}}_i}^{\left(w\right)})}^2$

满足自由度为(l-1，n-l)的F分布，若F≥F_a(α取0.05)，就以第k-1次的残差序列扣除本次自然基小波在时间域的周期展开，得到一个第k次残差时间序列，然后以波长递减直到k＝2终止。

6)用所建模型进行气候预测：将本发明中的卡尔曼滤波预测模型和针对卡尔曼滤波误差而进行的演化建模和自然基小波分别对各自的气候数据时间序列模型做下一步预测，并将预测结果相叠加，从而实现卡尔曼滤波及演化建模的短期气候预测。

Claims (3)

1.一种基于卡尔曼滤波及演化建模的短期气候预测方法，其特征在于，包括下列步骤：

A)收集气候预测模型中预测因子的历史数据，包括近年来的日照时数、温度、相对湿度和降雨量；

B)确定预测因子，整理与预测因子相关的历史数据，将历史数据分段处理，将每段数据的综合作用结果作为一个数据项目，同时将处理后的数据分为初值计算样本和检测样本两部分；

C)根据上述样本计算卡尔曼滤波初值：

Y_t=X_tβ_t+e_t

上式为卡尔曼滤波测量方程，首先利用多元线性回归计算回归系数β₀，再计算出动态噪声方差阵W和测量噪声方差阵V；

D)根据上述初值计算样本建立卡尔曼递推模型：利用步骤C)中计算的初值及卡尔曼递推式对检测样本的数据进行预测，并与检测样本进行比较，得到预测误差时间序列；其中，卡尔曼滤波递推模型对检测样本进行预测的同时，其卡尔曼滤波回归系数根据预测误差做出实时调整；

E)对上述预测误差进行预处理：

$x\left(t\right)=\stackrel{\‾}{x}\left(t\right)+\tilde{x}\left(t\right)$

其中，l为光滑参数；

将预测误差分为光滑部分和粗糙部分，并假设光滑部分是由宏观气候因素控制，粗糙部分是由微观气候因素控制，针对分解后的数据分别采用以下步骤F)和步骤G)的数学模型来描述；

F)针对步骤E)中的光滑部分时间序列项，用染色体来表示问题的可行解，包括可行解的形式和参数，其中，染色体表示的微分方程模型为一个树表达式，其微分方程解由Runge-Kutta方法求出；然后利用遗传操作在解空间中进行搜索，最后根据拟合表达式与实际值之间的偏差平方和进行优胜劣汰，求得最终的解，即拟合光滑部分时间序列最好的非线性表达式；

G)针对步骤E)中的粗糙部分时间序列项，用自然基小波对其进行分析，采用的是多尺度自然分形模型，从这些时间序列中找出其受气候因子影响下的周期波动规律，实现粗糙部分时间序列的拟合和预测功能；

H)将步骤D)中实现的卡尔曼滤波以及步骤F)和步骤G)中分别针对卡尔曼滤波误差而进行的演化建模和自然基小波分别对各自的气候数据时间序列模型做下一步预测，并将预测结果相叠加，从而实现卡尔曼滤波及演化建模的短期气候预测。

2.根据权利要求1所述的一种基于卡尔曼滤波及演化建模的短期气候预测方法，其特征在于步骤E)中，可根据具体情况调整光滑参数l的值，l可用于调节时间序列

的光滑程度，l越大，

越光滑。

3.根据权利要求1所述的一种基于卡尔曼滤波及演化建模的短期气候预测方法，其特征在于，步骤H)中，短期气候预测模型为线性预测模型与非线性预测模型相结合的方式。

Images