CN107765347A - 一种高斯过程回归和粒子滤波的短期风速预测方法 - Google Patents

Abstract

本发明公开了一种高斯过程回归和粒子滤波的短期风速预测方法，实现对异常值的在线动态检测与修正并提高了风速预测精度。首先，采用偏自相关函数确定与待预测时刻风速值相关性最大的输入变量集合，确定状态向量并构建合适的训练样本集合，在训练样本集合中建立高斯过程回归短期风速预测模型，并给出训练过程拟合残差；然后结合状态向量和高斯过程回归模型，建立粒子滤波状态空间方程，采用粒子滤波算法对当前量测值进行状态估计；最终，对粒子滤波的估计值和量测值残差进行分析，并根据“3σ”原则判断并修正异常值。本发明提供的方法能够有效地对异常值进行检测与修正，提高了短期风速预测精度，能够较好地解决电力***风速预测问题。

Description

一种高斯过程回归和粒子滤波的短期风速预测方法

技术领域

本发明涉及一种电力***短期风速预测方法，对电力***风速进行预测，属于电力***技术领域。

背景技术

目前用于风速预测主要有数值天气预报和统计模型两类方法。数值天气预报需要建立物理模型，通过微观气象学理论和计算流体力学得到风速、风向、温湿度等信息。统计模型方法主要采用数学统计的思想，通过挖掘数据间存在的内在规律进行预测。此类方法主要有时间序列、神经网络、支持向量机、卡尔曼滤波等。由于风速具有典型非线性、强波动性和强随机性特征，基于时间序列分析的线性模型难以描述风速变化特征。基于BP神经网络的短期风速预测模型具有较好的预测精度，但过程基于“黑箱”原理，难于建立显性的数学表达。支持向量机(support vector machines,SVM)模型采用结构风险最小化代替神经网络的经验风险最小化，从而具有更好的泛化能力，适应于处理小样本、高维、非线性回归问题，但模型超参数训练过程耗时长，在一定程度上限制了其广泛应用。通过时间序列分析建立的卡尔曼风速预测模型，取得了理想的预测精度。但卡尔曼滤波方法适用于线性数学模型，对于非线性过程处理能力弱。高斯过程回归(Gaussian process regression,GPR)输出具有概率分布特点，与SVM和BP神经网络模型相比，GPR具有更好预测精度。因此，本发明建立基于GPR的短期风速预测模型。

同时，短期风速预测往往存在两个重要问题：1)历史风速序列混杂的噪声影响了预测模型精度。风速序列在采集、传输、存储等过程中不可避免的会受到各种噪声因素的影响，如量测设备采集错误、数据传输过程中数据丢失等。在训练预测模型时，这些异常值将导致预测值偏离真实值，数据缺失严重情况下甚至难以给出预测结果。同时，模型参数的不准确估计也降低了预测精度。统计学方法指出：在训练过程中，拟合值远远偏离真实值的样本点即为异常值点。由此，产生了基于偏差的异常值检测方法。该方法主要过程为：首先利用已知数据建立数学模型，根据拟合数据和真实值间的残差来判断数据是否异常。其假设条件是拟合残差服从均值为零的高斯分布。2)由于统计模型方法机理自身存在的滞后性，导致预测值变化趋势落后于风速真实值，尤其在风速发生突变时刻，预测模型输出值往往难以反映真实值，需要对预测值进行修正。

针对风速预测过程中存在的第一个问题，本发明采用基于偏差的异常值检测方法。在利用GPR拟合真实值过程中，由于风速的强随机性和不确定性、预测模型自身的滞后性，使得拟合值和真实值有较大偏差，从而“淹没”了异常值处的偏差，导致异常值检测困难或造成误检、漏检。为消除正常值处较大的偏差，本发明引入了非线性非高斯滤波方法-粒子滤波(particle filter,PF)，从而提出基于高斯过程回归和粒子滤波(GPR-PF)相结合的短期风速预测模型，实现了异常值在线动态检测与修正。通过算例验证了本发明方法能够有效的对异常值进行检测并修正。

发明内容

发明目的：针对历史风速序列存在的噪声或数据缺失影响风速预测精度的问题，本发明提供一种基于高斯过程回归和粒子滤波的短期风速预测方法，实现对异常值的在线动态检测与修正，从而准确估计模型参数并提高风速预测精度。为有效选取出与待预测时刻风速值相关性较大的输入变量集合，本发明采用偏自相关函数度量两变量间的相关性并确定训练样本集合，进一步的确定状态向量。在训练样本集中通过高斯过程回归建立状态空间方程，采用粒子滤波算法对当前量测值进行状态估计，对估计值和量测值的残差进行分析，并根据“3σ”原则判断异常值。然后，修正异常值，并对清洗后的风速序列重新建立高斯过程回归预测模型。在进行提前15分钟风速预测时，同样采用粒子滤波算法对最新量测值进行状态估计，实现了异常值在线检测并修正。最终，算例分析结果表明，粒子滤波算法能够有效检测出异常风速值，降低了风速预测误差。

技术方案：一种基于高斯过程回归和粒子滤波的短期风速预测方法，包括以下步骤：

1)获取电力***短期风速预测所需的基本数据，并对原始数据进行零均值化处理；

(2)采用偏自相关函数确定与待预测时刻风速值相关性最大的输入变量集合，确定状态向量并构建合适的训练样本集合；

(3)在训练样本集合中建立高斯过程回归短期风速预测模型，并给出训练过程拟合残差；

(4)结合状态向量和高斯过程回归模型，建立粒子滤波状态空间方程，采用粒子滤波算法对当前量测值进行状态估计；

(5)对粒子滤波的估计值和量测值残差进行分析，并根据“3σ”原则判断并修正异常值；

(6)对清洗后的风速序列重新建立高斯过程回归预测模型。在进行提前15分钟风速预测时，同样采用粒子滤波算法对最新量测值进行状态估计，实现异常值在线检测并修正。

进一步地，步骤(1)中对原始风速时间序列进行零均值化处理，所述零均值化公式为：

式中：x(t)为原始风速时间序列，为序列x(t)的平均值。

进一步地，步骤(2)采用偏自相关函数确定与待预测时刻风速值相关性最大的输入变量集合，并确定状态向量与构建合适的训练样本集合，所述偏自相关函数计算方法为：

3.1假设x_i是输出变量，当滞后阶数为k时，偏自相关系数值在95％置信区间内，则x_i-k可以作为其中一个输入向量，如果所有偏自相关系数值在95％置信区间内，则认为x_i-1是输入变量；

3.2对于时间序列{x₁,x₂,L,x_n}，滞后阶数为k时的协方差定义为γ_k，(k＝0时，γ为方差)，计算公式如下：

式中：k＝0,1,2,L,M，是时间序列的均值；M＝n/4是最大滞后阶数；

3.3滞后阶数为k的自相关函数(autocorrelation function，ACF)定义为ρ_k：

式中：为滞后阶数k＝0时的协方差。

则滞后阶数为k时的PACF定义为α_kk：

式中：k＝1,2,L,M。

进一步地，步骤(3)在训练样本集合中建立高斯过程回归短期风速预测模型，并给出训练过程拟合残差，所述高斯过程回归模型预测过程为：

4.1假设训练样本集合为D＝{(x_i,y_i)|i＝1,2,3,···,n}＝(X,y)，其中：x_i∈R^m为m维输入向量，m×n维输入矩阵则可表示为X＝[x₁,x₂,···,x_n]，n表示训练样本点数量，y_i∈R为对应于x_i的输出标量；

4.2定义函数空间f(x)＝Φ(x)^Tω，f(x⁽¹⁾)、f(x⁽²⁾)、…、f(x⁽ⁿ⁾)构成随机变量的一个集合，且服从联合高斯分布，高斯过程模型就可以表示为：

式中：独立高斯白噪声服从均值为0，方差为σ²的高斯分布，记做ε:N(0,σ²)；δ_ij为Kronecker delta函数，当i＝j时，函数δ_ij＝1；m(x)为有限维分布族的均值函数，描述风速均值输出结果；k(x,x′)为协方差函数，刻画风速方差大小；

4.3GPR预测模型在n维训练集D内建立先验分布，在n_*维测试集D_*＝{(x_i,y_i)|i＝n+1,L,n+n_*}下转变为后验分布，则训练样本观测值y和测试数据的输出向量f_*之间构成联合高斯分布

其中，K(X,X)＝K_n表示N×N核矩阵，其元素K_ij＝k(x_i,x_j)；K(X,X_*)＝K(X_*,X)^T为测试数据X_*与训练集的输入X之间的协方差矩阵；K(X_*,X_*)为X_*自身的协方差；

4.4由此得出预测值f_*后验分布为

其中

均值向量为GPR模型风速预测均值，对应于点预测输出，为对应于的方差，由此可获得具有概率分布意义的风速区间预测结果。

进一步地，步骤(4)中结合状态向量和高斯过程回归模型，建立粒子滤波状态空间方程，采用粒子滤波算法对当前量测值进行状态估计。所述的粒子滤波状态空间方程为：

式中：H＝[1 0 0 0]；非线性状态转移函数GP(g)反映由k时刻状态X_k得到k+1时刻风速预测值，该函数由训练样本集合通过训练GPR模型参数得到；W_k和V_k+1分别为***过程噪声和观测噪声；X(k)＝[X₁(k)X₂(k)X₃(k)X₄(k)]^T。

进一步地，步骤(5)中对粒子滤波的估计值和量测值残差进行分析，并根据“3σ”原则判断并修正异常值，所述“3σ”原则具体过程为：

6.1对于测试数据r₁,r₂,L,r_n，取其算术平均值

及剩余误差值由此求出其均方根偏差为

6.2异常值判别依据为：若该值为异常数据；若则r_i为正常数据。

有益效果：本发明建立了基于高斯过程回归与粒子滤波的短期风速预测方法，利用粒子滤波非线性、非高斯滤波能力实现对原始风速序列中存在的异常值进行检测并修正，同时采用偏自相关函数确定最优输入变量集合与状态向量，避免了人工经验选取输入变量的不足。算例分析结果表明，本发明的风速预测方法能够实现异常值的在线动态检测与修正，对清洗后的时间序列建立短期风速预测模型，进一步增强了模型的预测性能。本发明的方法提供了一种风速时间序列异常值在线动态检测与修正方法，提高了短期风速预测精度，具有一定的工程应用意义。

附图说明

图1为自相关和偏自相关函数图；

图2为风速拟合结果；

图3为风速拟合结果局部放大图；

图4为不同模型拟合残差；

图5为提前15min风速预测结果；

图6为本发明预测方法的流程图。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

本发明的思路是将粒子滤波的非线性、非高斯滤波特征用于对原始历史风速序列中存在的异常值进行在线动态检测与修正，从而对清洗后的数据建立风速预测模型，提高预测精度。首先，为确定输入变量集合及状态向量，采用偏自相关函数度量变量间的相关性。其次，在训练样本集中通过高斯过程回归建立状态空间方程，采用粒子滤波算法对当前量测值进行状态估计，对估计值和量测值的残差进行分析，并根据“3σ”原则判断异常值。然后，修正异常值，并对清洗后的风速序列重新建立高斯过程回归预测模型。在进行提前15分钟风速预测时，同样采用粒子滤波算法对最新量测值进行状态估计，实现了异常值在线检测并修正。最终，通过算例分析结果表明本发明方法的有效性。

本发明采用GPR建立短期风速预测模型，GPR以贝叶斯理论和统计学习理论为基础，在处理高维数、小样本和非线性等复杂回归问题时具有易编程实现、超参数自适应获取以及输出概率分布等优点，从而在时间序列分析、动态***模型辨识、***控制等多领域获得了广泛应用。

应用GPR进行短期风速预测建模时，假设训练样本集合为D＝{(x_i,y_i)|i＝1,2,3,···,n}＝(X,y)，其中：x_i∈R^m为m维输入向量，m×n维输入矩阵则可表示为X＝[x₁,x₂,···,x_n]，n表示训练样本点数量，y_i∈R为对应于x_i的输出标量。用数学语言描述风速预测过程为：

定义函数空间f(x)＝Φ(x)^Tω，f(x⁽¹⁾)、f(x⁽²⁾)、…、f(x⁽ⁿ⁾)构成随机变量的一个集合，且服从联合高斯分布，高斯过程模型就可以表示为

式中：独立高斯白噪声服从均值为0，方差为σ²的高斯分布，记做ε:N(0,σ²)；δ_ij为Kroneckerdelta函数，当i＝j时，函数δ_ij＝1；m(x)为有限维分布族的均值函数，描述风速均值输出结果；k(x,x′)为协方差函数，刻画风速方差大小。

为简化推导，风速均值m(x)进行数据预处理使之为0。GPR预测模型在n维训练集D内建立先验分布，在n_*维测试集D_*＝{(x_i,y_i)|i＝n+1,L,n+n_*}下转变为后验分布，则训练样本观测值y和测试数据的输出向量f_*之间构成联合高斯分布

其中，K(X,X)＝K_n表示N×N的核矩阵，其元素K_ij＝k(x_i,x_j)；K(X,X_*)＝K(X_*,X)^T为测试数据X_*与训练集的输入X之间的协方差矩阵；K(X_*,X_*)为X_*自身的协方差。

由此得出预测值f_*后验分布为

其中

均值向量为GPR模型风速预测均值，对应于点预测输出，为对应于的方差，由此可获得具有概率分布意义的风速区间预测结果。

PF算法具有良好的非线性非高斯***状态滤波能力，并且随机量无需满足高斯分布的制约条件，因而在信号处理、通信、人工智能等领域获得应用。本发明选择粒子滤波算法对风速序列存在的异常值进行处理，提出基于GPR-PF的短期风速预测模型。

PF基本思想是用一组粒子近似表示***的后验概率分布，然后用这一近似的表示估计非线性***的状态。用非线性***动态状态空间模型描述粒子滤波过程为

式中：x_k和z_k分别为k时刻***状态向量和量测值；f(g)和h(g)分别为***状态转移函数和测量模型函数；w_k-1和ν_k分别为***过程噪声和观测噪声；从重要性密度函数q(x_k|x_1:k-1,z_k)抽样得到N个样本，并将这些样本表示为：状态k时刻的预测值为即

式中：是***已知噪声分布中采样得到的对应于粒子i的独立样本。完成预测阶段的所有粒子就构成k时刻的先验概率样本，表示为也就是得到了先验概率密度p(x_k|Z_k-1)。在获得新的观测量z_k后，各个粒子根据权值公式进行更新

对权值进行归一化处理

利用重采样方法克服粒子退化现象。根据每个粒子归一化权值对样本重新采样，复制权值较大的粒子，删除权值较小的粒子，得到等权值粒子集，即可得到后验分布密度函数

式中：δ(g)是狄拉克Kronecker delta函数，一个循环迭代过程结束。

最小均方误差准则意义下的状态最优估计为

自相关函数和偏自相关函数在识别模型类型及估计阶数过程中具有重要的意义。本发明根据自相关函数和偏自相关函数衡量X_k与X_k-τ间的相关关系，从而有效分析出时间延迟，并确定输入变量集合及状态向量。其中，τ为延迟时间。

将江苏省某风电场实测风速(数据采样时间间隔为15min)用于GPR-PF短期风速预测建模，图1为时间序列分析的自相关函数图和偏自相关函数图。从图中可以看出，自相关函数具有拖尾特征，而偏自相关函数截尾，从而风速序列满足AR模型。结合偏自相关函数图，本文选取4个输入变量，即预测k+1时刻的风速X_k+1时，将k,k-1,k-2,k-3时刻的风速X_k,X_k-1,X_k-2,X_k-3作为输入变量。令X₁(k)＝X(k)，X₂(k)＝X(k-1)，X₃(k)＝X(k-2)，X₄(k)＝X(k-3)，则k时刻状态向量为X(k)＝[X₁(k)X₂(k)X₃(k)X₄(k)]^T。

在采用粒子滤波进行风速序列异常值检测与修正时，首先采用历史风速值预测当前时刻风速状态，然后根据当前时刻风速量测值对预测值进行修正，从而获得此时刻量测值的最优估计，并获得风速估计值与量测值间的残差。

结合偏自相关函数确定的状态向量，本发明的GPR-PF短期风速预测状态空间模型为：

式中：H＝[1 0 0 0]；非线性状态转移函数GP(g)反映由k时刻状态X_k得到k+1时刻风速预测值，该函数由训练样本集合通过训练GPR模型参数得到。

粒子滤波通过对k时刻的状态来估计得到k+1时刻风速估计值，并根据k+1时刻量测值修正估计值，从而获得了滤波后风速最优估计值。

本发明采用基于偏差的数据异常值检测方法检测并修正原始风速序列中存在的异常值，通过分析GPR-PF风速估计值和量测值之间存在的残差r，进而根据“3σ”准则来判断异常值。

对于测试数据r₁,r₂,L,r_n，取其算术平均值

及剩余误差值由此求出其均方根偏差为

则异常值判别依据为：若该值为异常数据；若则r_i为正常数据。

为量化预测值接近真实值的程度，本文选择平均绝对百分比误差(mean absolutepercentage error,MAPE)和均方根误差(root mean square error,RMSE)作为模型预测效果评价指标，计算公式分别为：

式中：T为预测点个数，y_i为第i个预测点风速真实值，为第i个预测点模型预测值。

本发明采用某风电场2008年5月14日12:00时至5月25日23:45时共1104个实测风速值作为训练样本序列，数据采样时间间隔为15min，建立GPR-PF风速预测模型，对5月26日的96个风速值做提前一步(即提前15min)预测。

首先，由训练样本集建立GPR短期风速预测模型并求解模型超参数，建立状态空间方程，分析模型拟合值与量测值间的残差，从而采用基于偏差的方法进行异常值检测。图2为训练过程中采用GPR与GPR-PF两种方法得到的预测结果。从图中可以看出，采用GPR模型进行短期风速预测时，由于风速强波动性与随机性，使得GPR模型难以很好的跟踪风速变化趋势，存在一定的滞后性，从而产生较大的偏差。而GPR-PF模型采用当前时刻量测值对风速估计值进行动态更新，能够很好的估计出风速的真实状态。从图3预测结果局部放大图可以更细致比较两种模型的拟合效果。对风速预测值与实际量测值进行残差分析，图4分别为采用GPR和GPR-PF模型进行拟合时的残差结果。单独采用GPR模型时，可以看出残差分布较分散，不利于异常值的分析与检测。而采用GPR-PF混合方法时，残差分布较集中，可较容易的采用“3σ”准则判断异常值。

通过对残差分布分析及“3σ”准则，对异常值检测与修正结果如表1所示。

表1异常值检测与修正结果

对清洗后的风速序列重新建立GPR预测模型，并进行提前15min的短期风速预测。图5为采用BP神经网络、SVM、GPR和GPR-PF四种模型的风速预测结果。预测模型定量评价指标结果如表2所示。可以看出，相对于BP神经网络、SVM两种预测模型，GPR可以给出较好的预测精度。通过采用粒子滤波算法对异常值进行修正后，本发明的GPR-PF风速预测模型削弱了噪声影响，从而获得最优预测结果。

表2提前15min风速预测误差

综上所述，本发明基于高斯过程回归与粒子滤波的短期风速预测方法具有如下优势：1)采用偏自相关函数选取输入变量集合并确定状态向量，从而避免了人工经验选取输入变量的不足。2)建立基于粒子滤波的状态空间方程，实现对历史风速序列中存在的异常值进行在线动态检测与修正，对清洗后的风速时间序列建立高斯过程回归短期风速预测模型，从而进一步提高了预测精度。3)相对于支持向量机和BP神经网络预测方法，本发明采用高斯过程回归建立短期风速预测模型，具有更好的预测性能表现，其模型超参数可以自适应获取。本发明方法可以对原始时间序列中存在的异常值进行在线动态检测与修正，从而提高了短期风速预测精度，对电力***安排风功率出力计划及保证电网安全稳定运行具有一定的参考价值。

Claims (6)

1.一种高斯过程回归和粒子滤波的短期风速预测方法，其特征在于：包括以下步骤：

(1)获取电力***短期风速预测所需的基本数据，并对原始风速时间序列进行零均值化处理；

(2)采用偏自相关函数确定与待预测时刻风速值相关性最大的输入变量集合，确定状态向量并构建合适的训练样本集合；

(3)在训练样本集合中建立高斯过程回归短期风速预测模型，并给出训练过程拟合残差；

(4)结合状态向量和高斯过程回归模型，建立粒子滤波状态空间方程，采用粒子滤波算法对当前量测值进行状态估计；

(5)对粒子滤波的估计值和量测值残差进行分析，并根据“3σ”原则判断并修正异常值；

(6)对清洗后的风速序列重新建立高斯过程回归预测模型，在进行提前15分钟风速预测时，同样采用粒子滤波算法对最新量测值进行状态估计，实现异常值在线检测并修正。

2.如权利要求1所述的基于高斯过程回归和粒子滤波的短期风速预测方法，其特征在于：步骤(1)中对原始风速时间序列进行零均值化处理，所述零均值化公式为：

<mrow> <msup> <mi>x</mi> <mo>&amp;prime;</mo> </msup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <mover> <mi>x</mi> <mo>&amp;OverBar;</mo> </mover> </mrow>

式中：x(t)为原始风速时间序列，为序列x(t)的平均值。

3.如权利要求1所述的基于高斯过程回归和粒子滤波的短期风速预测方法，其特征在于：步骤(2)采用偏自相关函数确定与待预测时刻风速值相关性最大的输入变量集合，并确定状态向量与构建合适的训练样本集合，所述偏自相关函数计算方法为：

3.1假设x_i是输出变量，当滞后阶数为k时，偏自相关系数值在95％置信区间内，则x_i-k可以作为其中一个输入向量，如果所有偏自相关系数值在95％置信区间内，则认为x_i-1是输入变量；

3.2对于时间序列{x₁,x₂,L,x_n}，滞后阶数为k时的协方差定义为γ_k，(k＝0时，γ为方差)，计算公式如下：

<mrow> <msub> <mover> <mi>&amp;gamma;</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>n</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mi>k</mi> </mrow> </munderover> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <mover> <mi>x</mi> <mo>&amp;OverBar;</mo> </mover> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mrow> <mi>i</mi> <mo>-</mo> <mi>k</mi> </mrow> </msub> <mo>-</mo> <mover> <mi>x</mi> <mo>&amp;OverBar;</mo> </mover> <mo>)</mo> </mrow> </mrow>

式中：k＝0,1,2,L,M，是时间序列的均值；M＝n/4是最大滞后阶数；

3.3滞后阶数为k的自相关函数定义为ρ_k：

<mrow> <msub> <mover> <mi>&amp;rho;</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> <mo>=</mo> <mfrac> <msub> <mover> <mi>&amp;gamma;</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> <msub> <mover> <mi>&amp;gamma;</mi> <mo>^</mo> </mover> <mn>0</mn> </msub> </mfrac> </mrow>

式中：为滞后阶数k＝0时的协方差；

则滞后阶数为k时的PACF定义为α_kk：

<mfenced open = "" close = "}"> <mtable> <mtr> <mtd> <mrow> <msub> <mover> <mi>&amp;alpha;</mi> <mo>^</mo> </mover> <mn>11</mn> </msub> <mo>=</mo> <msub> <mover> <mi>&amp;rho;</mi> <mo>^</mo> </mover> <mn>1</mn> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>&amp;alpha;</mi> <mo>^</mo> </mover> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mover> <mi>&amp;rho;</mi> <mo>^</mo> </mover> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>k</mi> </msubsup> <msub> <mover> <mi>&amp;rho;</mi> <mo>^</mo> </mover> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>-</mo> <mi>j</mi> </mrow> </msub> <msub> <mover> <mi>&amp;alpha;</mi> <mo>^</mo> </mover> <mrow> <mi>k</mi> <mi>j</mi> </mrow> </msub> </mrow> <mrow> <mn>1</mn> <mo>-</mo> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>k</mi> </msubsup> <msub> <mover> <mi>&amp;rho;</mi> <mo>^</mo> </mover> <mi>j</mi> </msub> <msub> <mover> <mi>&amp;alpha;</mi> <mo>^</mo> </mover> <mrow> <mi>k</mi> <mi>j</mi> </mrow> </msub> </mrow> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mover> <mi>&amp;alpha;</mi> <mo>^</mo> </mover> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>=</mo> <msub> <mover> <mi>&amp;alpha;</mi> <mo>^</mo> </mover> <mrow> <mi>k</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>-</mo> <msub> <mover> <mi>&amp;alpha;</mi> <mo>^</mo> </mover> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>&amp;CenterDot;</mo> <msub> <mover> <mi>&amp;alpha;</mi> <mo>^</mo> </mover> <mrow> <mi>k</mi> <mo>,</mo> <mi>k</mi> <mo>-</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>,</mo> <mrow> <mo>(</mo> <mi>j</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced>

式中：k＝1,2,L,M。

4.如权利要求1所述的基于高斯过程回归和粒子滤波的短期风速预测方法，其特征在于：步骤(3)在训练样本集合中建立高斯过程回归短期风速预测模型，并给出训练过程拟合残差，所述高斯过程回归模型预测过程为：

4.1假设训练样本集合为D＝{(x_i,y_i)i＝1,2,3,…,n}＝(X,y)，其中：x_i∈R^m为m维输入向量，m×n维输入矩阵则可表示为X＝[x₁,x₂,···,x_n]，n表示训练样本点数量，y_i∈R为对应于x_i的输出标量；

4.2定义函数空间f(x)＝Φ(x)^Tω，f(x⁽¹⁾)、f(x⁽²⁾)、…、f(x⁽ⁿ⁾)构成随机变量的一个集合，且服从联合高斯分布，高斯过程模型就可以表示为：

<mrow> <mi>y</mi> <mo>~</mo> <mi>GP</mi> <mrow> <mo>(</mo> <mi>m</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>,</mo> <mi>k</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <msup> <mi>x</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>n</mi> <mn>2</mn> </msubsup> <msub> <mi>&amp;delta;</mi> <mi>ij</mi> </msub> <mo>)</mo> </mrow> </mrow>

式中：独立高斯白噪声服从均值为0，方差为σ²的高斯分布，记做ε:N(0,σ²)；δ_ij为Kroneckerdelta函数，当i＝j时，函数δ_ij＝1；m(x)为有限维分布族的均值函数，描述风速均值输出结果；k(x,x′)为协方差函数，刻画风速方差大小；

4.3GPR预测模型在n维训练集D内建立先验分布，在n_*维测试集D_*＝{(x_i,y_i)|i＝n+1,L,n+n_*}下转变为后验分布，则训练样本观测值y和测试数据的输出向量f_*之间构成联合高斯分布

<mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mi>y</mi> </mtd> </mtr> <mtr> <mtd> <msub> <mi>f</mi> <mo>*</mo> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>~</mo> <mi>N</mi> <mrow> <mo>(</mo> <mn>0</mn> <mo>,</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mi>K</mi> <mrow> <mo>(</mo> <mi>X</mi> <mo>,</mo> <mi>X</mi> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>n</mi> <mn>2</mn> </msubsup> <mi>I</mi> </mrow> </mtd> <mtd> <mrow> <mi>K</mi> <mrow> <mo>(</mo> <mi>X</mi> <mo>,</mo> <msub> <mi>X</mi> <mo>*</mo> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>K</mi> <mrow> <mo>(</mo> <msub> <mi>X</mi> <mo>*</mo> </msub> <mo>,</mo> <mi>X</mi> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>K</mi> <mrow> <mo>(</mo> <msub> <mi>X</mi> <mo>*</mo> </msub> <mo>,</mo> <msub> <mi>X</mi> <mo>*</mo> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>)</mo> </mrow> </mrow>

其中，K(X,X)＝K_n表示N×N核矩阵，其元素K_ij＝k(x_i,x_j)；K(X,X*)＝K(X*,X)^T为测试数据X*与训练集的输入X之间的协方差矩阵；K(X*,X*)为X*自身的协方差；

4.4由此得出预测值f*后验分布为

<mrow> <msub> <mi>f</mi> <mo>*</mo> </msub> <mo>|</mo> <mi>X</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <msub> <mi>x</mi> <mo>*</mo> </msub> <mo>~</mo> <mi>N</mi> <mrow> <mo>(</mo> <mover> <msub> <mi>f</mi> <mo>*</mo> </msub> <mo>&amp;OverBar;</mo> </mover> <mo>,</mo> <mi>cov</mi> <mo>(</mo> <msub> <mi>f</mi> <mo>*</mo> </msub> <mo>)</mo> <mo>)</mo> </mrow> </mrow>

其中

<mrow> <mi>cov</mi> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mo>*</mo> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mi>k</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mo>*</mo> </msub> <mo>,</mo> <msub> <mi>x</mi> <mo>*</mo> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mi>K</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mo>*</mo> </msub> <mo>,</mo> <mi>X</mi> <mo>)</mo> </mrow> <mo>&amp;times;</mo> <msup> <mrow> <mo>&amp;lsqb;</mo> <mi>K</mi> <mrow> <mo>(</mo> <mi>X</mi> <mo>,</mo> <mi>X</mi> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>n</mi> <mn>2</mn> </msubsup> <msub> <mi>I</mi> <mi>n</mi> </msub> <mo>&amp;rsqb;</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mi>K</mi> <mrow> <mo>(</mo> <mi>X</mi> <mo>,</mo> <msub> <mi>x</mi> <mo>*</mo> </msub> <mo>)</mo> </mrow> </mrow>

均值向量为GPR模型风速预测均值，对应于点预测输出，为对应于的方差，由此可获得具有概率分布意义的风速区间预测结果。

5.如权利要求1所述的基于高斯过程回归和粒子滤波的短期风速预测方法，其特征在于：步骤(4)中结合状态向量和高斯过程回归模型，建立粒子滤波状态空间方程，采用粒子滤波算法对当前量测值进行状态估计。所述的粒子滤波状态空间方程为：

<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <msub> <mi>X</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <mi>G</mi> <mi>P</mi> <mo>(</mo> <msub> <mi>X</mi> <mi>k</mi> </msub> <mo>)</mo> <mo>+</mo> <msub> <mi>W</mi> <mi>k</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>y</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <msub> <mi>HX</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>+</mo> <msub> <mi>V</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced>

式中：H＝[1 0 0 0]；非线性状态转移函数GP(g)反映由k时刻状态X_k得到k+1时刻风速预测值，该函数由训练样本集合通过训练GPR模型参数得到；W_k和V_k+1分别为***过程噪声和观测噪声；X(k)＝[X₁(k) X₂(k) X₃(k) X₄(k)]^T。

6.如权利要求1所述的基于高斯过程回归和粒子滤波的短期风速预测方法，其特征在于：步骤(5)中对粒子滤波的估计值和量测值残差进行分析，并根据“3σ”原则判断并修正异常值，所述“3σ”原则具体过程为：

6.1对于测试数据r₁,r₂,L,r_n，取其算术平均值

<mrow> <mover> <mi>r</mi> <mo>&amp;OverBar;</mo> </mover> <mo>=</mo> <mfrac> <mn>1</mn> <mi>n</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msub> <mi>r</mi> <mi>i</mi> </msub> </mrow>

及剩余误差值由此求出其均方根偏差为

<mrow> <mi>&amp;sigma;</mi> <mo>=</mo> <msqrt> <mrow> <mfrac> <mn>1</mn> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msup> <msub> <mi>v</mi> <mi>i</mi> </msub> <mn>2</mn> </msup> </mrow> </msqrt> </mrow>

6.2异常值判别依据为：若该值为异常数据；若则ri为正常数据。