基于G‑L混合噪声特性核岭回归技术的风速预报方法及装置与流程

本发明涉及短期风速预报技术领域，具体涉及基于G-L混合噪声特性核岭回归技术的短期风速预报方法及装置。

背景技术：

对于线性系统而言,从Gauss时代起，就利用最小二乘技术把平面上的点拟合成直线，把高维空间的点拟合成超平面。经历了200多年的发展，经典最小二乘技术已经成为许多领域数据处理的最广泛使用的技术。但是，对于线性回归或非线性回归中的不适定问题，基于最小二乘回归技术的性能会变得很坏，针对这种情况，众多学者研究了最小二乘回归的改进模型，提出了许多新的回归算法。岭回归(Ridge regression，简记为RR)就是其中之一，岭回归模型是一种专门用于共线性数据分析的有偏估计回归方法，实质上是一种改良的最小二乘估计法，通过放弃最小二乘技术的无偏性，以损失部分信息、降低精度为代价获得回归系数更为符合实际、更可靠的回归方法，对病态数据的耐受性远远强于最小二乘技术。岭回归算法自A.E.Hoerl和R.W.Kennard于1962年提出以来，就得到了广泛的关注，成功应用于工程技术和社会科学等各个方面。设给定风速数据：

D_l＝{(x₁,y₁),(x₂,y₂),…,(x_l,y_l)} (1)

其中x_i∈Rⁿ,y_i∈R,i＝1,2,…,l，多元线性回归模型为f(x)＝ω^T·x+b，其中x＝(x₁,x₂,…,x_l)^T，参数向量ω∈Rⁿ决定最小二乘回归和岭回归模型，其中x_i∈X＝Rⁿ，Rⁿ表示n维欧式空间，R表示实数集，l表示样本个数，上标T表示转置运算。通过最小化目标函数：

其中C＞0是惩罚参数，目标函数(2)中假设样本噪声服从Gauss分布，其均值为0，方差为σ²，即y_i＝f(x_i)+ξ_i,i＝1,…,l，ξ_i～N(0,σ²)。式(2)中的g_GN-RR表示岭回归模型的目标函数。

最小二乘回归机和岭回归机在解决线性回归问题中取得了良好的效果。在实际应用中，样本D_l中x_i与y_i一般不满足线性关系，这样线性GN-RR在解决非线性回归问题时往往不能取得预期的效果。2000年，Suykens等提出的基于Gauss噪声模型的核岭回归机，通过核技术构造合适的核变换Φ:Rⁿ→H(H为Hilbert空间)，利用Hilbert空间中的內积运算引进非线性核函数K(x_i,x_j)＝(Φ(x_i)·Φ(x_j))，从而把非线性问题转化为线性问题来解决。一般地，基于Gauss噪声特性核岭回归机通过最小化目标函数来求解：

经典的回归模型(线性最小二乘回归、核岭回归、最小二乘支持向量回归等模型)一般都假设数据中的噪声分布服从高斯同方差分布。当噪声服从Gauss分布时，基于Gauss噪声特性的核岭回归机(Kernel ridge regression based on Gauss-noise，简记为GN-KRR)能够取得预期的效果。研究表明在许多实际应用领域中，噪声分布不一定服从高斯同方差分布，而服从Beta分布、拉普拉斯分布、高斯异方差分布、威布尔分布，或其他类型噪声特性的分布。尤其在风速/风功率预报中，噪声分布随着季节和地域而变化，需要由某种混合分布来联合表示，如利用Gauss-Laplace(简记为G-L)混合噪声分布来拟合不确定数据中的未知噪声特性。此时应用ν-SVR、GN-KRR等经典回归技术进行预测，则预报结果不能满足实际领域中对风速预报精度的要求。

技术实现要素：

本发明提供了一种基于G-L混合噪声特性核岭回归技术的风速预报方法及装置，以解决现有的单一的噪声特性的核岭回归技术不能满足实际应用中对风速预报精度要求的问题。

为解决上述技术问题，本发明的基于G-L混合噪声特性核岭回归技术的风速预报方法包括如下步骤：

1)获取某地区具有噪声影响的风速数据集D_l，利用Bayesian原理，得到基于G-L混合噪声特性的损失函数c(ξ)；

2)利用统计学习理论和最优化理论，结合步骤1)中得到的基于G-L混合噪声特性的损失函数c(ξ)，建立基于G-L混合噪声特性核岭回归模型的原问题，推导并求解出基于G-L混合噪声特性核岭回归模型的对偶问题；

3)利用十折交叉验证技术确定基于G-L混合噪声特性核岭回归模型的对偶问题的惩罚参数及权重参数，选取核函数K(·,·)；构造基于G-L混合噪声特性核岭回归模型的决策函数f(x)；

4)构造基于G-L混合噪声特性核岭回归模型的风速预报模式：输入向量为其中i,j为相关联的两个时刻；输出值为x_i+step，其中step为预测间隔时间，利用该预报模式预报分析某一时刻i以后step间隔时刻的风速值。

所述基于G-L混合噪声特性的损失函数为其中，权重参数λ₁≥0,λ₂≥0且λ₁+λ₂＝1，ξ为已知的独立同分布的随机变量。

所述基于G-L混合噪声特性的核岭回归模型的原问题为：

s.t.ξ_i＝y_i-ω^T·Φ(x_i)-b

ξ_i≥0,i＝1,…,l

其中，C＞0是惩罚参数，权重参数λ₁,λ₂≥0且λ₁+λ₂＝1，l为样本个数，y_i(i＝1,…,l)为测量值，ξ_i(i＝1,…,l)为已知的独立同分布的随机变量，Φ:Rⁿ→H为核变换，H为Hilbert空间，ω为参数向量，T为转置，b∈R；式中的P_GLM-KRR表示基于G-L混合噪声特性核岭回归模型的原问题，表示基于G-L混合噪声特性核岭回归模型的原问题的目标函数。

步骤2)中，构造Lagrange泛函：

应用最优化理论，得到基于G-L混合噪声特性核岭回归模型原问题的对偶问题；其中，α_i(i＝1,…,l)为拉格朗日乘子。

所述基于G-L混合噪声特性的核岭回归模型对偶问题为：

其中，l为样本个数，惩罚参数C＞0，权重参数λ₁,λ₂≥0且λ₁+λ₂＝1；s.t.为subject to的缩写，表示基于G-L混合噪声特性核岭回归模型对偶问题的目标函数，α_i,α_j(i,j＝1,…,l)为拉格朗日乘子，y_i(i＝1,…,l)为测量值，K(x_i,x_j)为核函数。

基于G-L混合噪声特性核岭回归模型的决策函数为：

其中，Φ:Rⁿ→H为核变换，H为Hilbert空间，K(x_i,x_j)＝(Φ(x_i)·Φ(x_j))，ω∈Rⁿ为参数向量，(Φ(x_i)·Φ(x_j))表示H空间中的内积。

本发明基于G-L混合噪声特性核岭回归技术的风速预报装置包括损失函数获取模块、对偶问题求解模块、决策函数构造模块及风速预报模块；

所述损失函数获取模块，用于获取某地区具有噪声影响的风速数据集D_l，利用Bayesian原理，得到基于G-L混合噪声特性的损失函数c(ξ)；

所述对偶问题求解模块，用于利用统计学习理论和最优化理论，结合基于G-L混合噪声特性的损失函数c(ξ)，建立基于G-L混合噪声特性核岭回归模型的原问题，推导并求解出基于G-L混合噪声特性核岭回归模型对偶问题；

所述决策函数构造模块，用于利用十折交叉验证技术确定基于G-L混合噪声特性核岭回归模型对偶问题的惩罚参数及权重参数，选取核函数K(·,·)；构造基于G-L混合噪声特性核岭回归模型的决策函数f(x)；

所述风速预报模块，用于构造基于G-L混合噪声特性核岭回归模型的风速预报模式：输入向量为其中i,j为相关联的两个时刻；输出值为x_i+step，其中step为预测间隔时间，利用该预报模式预报分析某一时刻i以后step间隔时刻的风速值。

根据建立的基于G-L混合噪声特性的核岭回归模型的原问题，构造Lagrange泛函，利用最优学习理论，得到基于G-L混合噪声特性核岭回归模型原问题的对偶问题。

所述基于G-L混合噪声特性的损失函数其中，权重参数λ₁≥0,λ₂≥0且λ₁+λ₂＝1，ξ为已知的独立同分布的随机变量。

基于G-L混合噪声特性的核岭回归模型的决策函数为：

其中，Φ:Rⁿ→H为核变换，H为Hilbert空间，K(x_i,x_j)＝(Φ(x_i)·Φ(x_j))，ω∈Rⁿ为参数向量，(Φ(x_i)·Φ(x_j))表示H空间中的内积，l为样本个数，α_i(i＝1,…,l)为拉格朗日乘子。

本发明的有益效果是：本发明提出一种应用基于G-L混合噪声特性核岭回归技术进行风速预报分析的新的统一理论模型。通过Bayesian原理推导出基于G-L混合噪声特性的损失函数，在此基础上利用统计学习理论、最优化理论构造基于G-L混合噪声特性核岭回归模型，该核岭回归模型具有较高稳定性和鲁棒性；最后利用基于G-L混合噪声特性核岭回归技术得到风速预报；本发明的方法及装置能够满足实际应用中，如风力发电、农业生产等，对风速预报精度的要求。

在风速预报中，需要根据噪声特性的实际分布确定预报技术，利用基于G-L混合噪声特性的核岭回归技术、基于考虑Gauss噪声特性的核岭回归技术和ν-支持向量回归技术进行风速预报，表现出了较好的预报效果。

附图说明

图1 Gauss概率密度函数(pdf)、Laplace概率密度函数(pdf)、Gauss损失函数和Laplace损失函数曲线图；

图2不同参数值下的Gauss-Laplace损失函数；

图3 10分钟以后风速预报结果ν-SVR(C＝181,step＝1)；

图4 10分钟以后风速预报结果GN-KRR(C＝181,step＝1)；

图5 10分钟以后风速预报结果GLM-KRR(C＝181,λ₁＝0.5,λ₂＝0.5,step＝1)；

图6 30分钟以后风速预报结果ν-SVR(C＝181,step＝3)；

图7 30分钟以后风速预报结果GN-KRR(C＝181,step＝3)；

图8 30分钟以后风速预报结果GLM-KRR(C＝181,λ₁＝0.5,λ₂＝0.5,step＝3)；

图9 60分钟以后风速预报结果ν-SVR(C＝181,step＝6)；

图10 60分钟以后风速预报结果GN-KRR(C＝181,step＝6)；

图11 60分钟以后风速预报结果GLM-KRR(C＝181,λ₁＝0.5,λ₂＝0.5,step＝6)；

图12基于G-L混合噪声特性核岭回归技术的风速预报方法总体流程图。

具体实施方式

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

本发明的基于G-L混合噪声特性核岭回归技术的风速预报方法实施例

该方法包括以下步骤：

1)获取某地区具有噪声影响的风速数据集D_l，利用Bayesian原理，得到Gauss-Laplace(简记为G-L)混合噪声特性的损失函数c(ξ)；

2)利用统计学习理论和最优化理论，结合步骤1)中得到的基于G-L混合噪声特性的损失函数，建立基于G-L混合噪声特性核岭回归模型的原问题，推导并求解出基于G-L混合噪声特性核岭回归模型对偶问题；

3)利用十折交叉验证技术确定基于G-L混合噪声特性核岭回归模型对偶问题的最优参数C、λ₁及λ₂，选取核函数K(·,·)；构造基于G-L混合噪声特性核岭回归模型的决策函数f(x)；

4)构造基于G-L混合噪声特性核岭回归模型的风速预报模式：输入向量为其中i,j为相关联的两个时刻，输出值为x_i+step，其中step为预测间隔时间；利用该预报模式预报分析某一时刻i以后step间隔时刻的风速值。

下面对上述步骤作进一步详细介绍：

步骤1)中，获取具有噪声特性的数据集D_l＝{(x₁,y₁),(x₂,y₂),…,(x_l,y_l)}，其中x_i∈Rⁿ,y_i∈R,i＝1,2,…,l，Rⁿ表示n维欧式空间，R表示实数集，l表示样本个数；利用Bayesian原理，推导基于G-L混合噪声特性的最优损失函数，具体过程为：

给定具有噪声特性的数据集D_l，回归函数f(x)是未知的。一般地，最小化目标函数：

其中c(ξ_i)＝c(y_i-f(x_i))表示在样本点(x_i,y_i)∈D_l进行预测时所得到预测值f(x_i)与测量值y_i(i＝1,…,l)比较所产生的损失值，λ是正数。假设噪声特性是加性的，即y_i＝f_i(x_i)+ξ_i(i＝1,…,l)，且ξ_i(i＝1,…,l)是独立同分布(i.i.d.)的随机变量，ξ_i(i＝1,…,l)的方差为σ²，均值为μ。利用数据g∈D_l估计函数f(x)，根据Bayesian原理可得噪声特性的最优损失函数为：

c(x,y,f(x))＝-logp(y-f(x)) (5)

其中p(y-f(x))＝p(ξ)表示误差ξ的概率密度函数，c(x_i,y_i,f(x_i))＝c(ξ_i)(i＝1,…,l)表示在样本点(x_i,y_i)进行预测时所得到预测值f(x_i)与y_i比较所产生的损失值，c(ξ)表示损失函数。

如图1所示，Gauss噪声特性的损失函数为：

如图1所示，Laplace噪声特性的损失函数为：c(ξ)＝c(y-f(x))＝|ξ| (7)

如图2所示，Gauss-Laplace混合噪声特性的损失函数为：

步骤2)中，构造并求解最优化问题进一步具体为：

基于G-L混合噪声特性核岭回归模型的原问题为：

其中ξ_i＝y_i-ω^T·Φ(x_i)-b(i＝1,2,…,l)，C＞0是惩罚参数，权重参数λ₁,λ₂≥0且λ₁+λ₂＝1。式(9)中的P_GLM-KRR表示基于G-L混合噪声特性核岭回归模型的原问题，表示基于G-L混合噪声特性核岭回归模型原问题的目标函数。

构造Lagrange泛函L(ω,b,α,ξ)：

根据最优化理论，为求L(ω,b,α,ξ)的极小值，分别对ω,b,ξ求偏导数.由KKT(Karush-Kusn-Tucke)条件▽_ωL＝0,▽_bL＝0,▽_ξL＝0，得:

把上述极值条件代入L(ω,b,α,ξ)，并对α求极大值，可得到基于G-L混合噪声特性核岭回归模型原问题(9)的对偶问题(简记为GLM-KRR)为：

其中C＞0是惩罚参数，权重参数λ₁,λ₂≥0，且λ₁+λ₂＝1，且有

其中的D_GLM-KRR表示基于G-L混合噪声特性核岭回归模型的对偶问题，表示基于G-L混合噪声特性核岭回归模型对偶问题的目标函数。

拉格朗日乘子法是Powel和Hestenes于1969年针对等式约束优化问题同时独立提出的优化算法，也称为PH算法；其基本思想是：从原问题的拉格朗日函数出发，再加上适当的罚函数，从而将原问题转化为求解一系列的无约束优化子问题；增广拉格朗日乘子法(Augmented Lagrange Multiplier method，简记为ALM)是1973年Rockfellar将PH算法推广到求解不等式约束优化问题，也称为PHR算法；ALM法是解决同时带有等式和不等式约束问题的一类优化方法。其基本思想是：把解等式约束优化问题的乘子法推广到不等式约束优化问题，即先引进辅助变量把不等式约束转化为等式约束，然后再利用最有效条件消去辅助变量。

对于一个求解函数最小值的优化问题(求函数最大值也类似)，一般可以描述为下列数学规划模型：

式中x为决策变量，f(x)为目标函数，式为约束条件，U是基本空间，R是U的子集；满足约束条件的解X称为可行解，集合R表示所有满足约束条件的解所组成的集合，称为可行解集合；式(5)中的x与式(12)中x的含义相同，x＝(x₁,x₂,…,x_l)^T,y＝(y₁,y₂,…,y_l)^T,(x_i,y_i)∈D_l,i＝1,2,…,l，上标T表示转置。p(y-f(x))＝p(ξ)表示误差ξ的概率密度函数；c(x_i,y_i,f(x_i))＝c(ξ_i)表示在样本点(x_i,y_i)进行预测时所得到预测值f(x_i)与y_i比较所产生的损失值，c(ξ)表示损失函数。

步骤2)中所述的统计学习理论和最优化理论为现有技术，请参考文献《数据挖掘中的新方法：支持向量机》(作者：邓乃杨,田英杰；出版社：科学出版社,2004.6)、《The Nature of Statistical Learning Theory》(Vapnik V.New York:Springer-Verlag,1995)及《Pattern recognition and machine learning》(C.M.Bishop.Springer,New York,2006)这里不再详细说明。

步骤3)中，利用十折交叉验证技术确定最优参数C、λ₁、λ₂具体为：利用ALM法求解基于G-L混合噪声特性核岭回归模型对偶问题(11)及利用十折交叉验证技术确定最优参数C、λ₁、λ₂。提出的基于G-L混合噪声特性核岭回归机应用Matlab 7.8程序语言实现，模型GLM-KRR的参数C∈[1,201]，λ₁,λ₂≥0且λ₁+λ₂＝1。

选取合适的核函数K(·,·)是利用核技术构造核函数K(·,·)，把基于G-L混合噪声特性岭回归模型推广为基于G-L混合噪声特性核岭回归模型；其中K(x_i,x_j)＝(Φ(x_i)·Φ(x_j))，Φ:Rⁿ→H，H为Hilbert空间，(Φ(x_i)·Φ(x_j))为H空间中的内积。常见核函数有：

(1)多项式核函数：K(x_i,x_j)＝((x_i·x_j)+1)^d；

(2)Gauss径向基核函数：K(x_i,x_j)＝exp(-||x_i-x_j||²/σ²)；

其中d是正数，取d＝2或3；σ是正数，取σ＝0.2。

可得到基于G-L混合噪声特性核岭回归模型的决策函数为

其中Φ:Rⁿ→H(H为Hilbert空间)为核变换，K(x_i,x_j)＝(Φ(x_i)·Φ(x_j))，ω∈Rⁿ为参数向量，(Φ(x_i)·Φ(x_j))表示H空间中的内积。

步骤3)中所述的十折交叉验证技术为现有技术，请参考文献《数据挖掘中的新方法：支持向量机》(作者：邓乃杨,田英杰；出版社：科学出版社,2004.6)及《The Nature of Statistical Learning Theory》(Vapnik V.New York:Springer-Verlag,1995)，这里不再详细说明。

步骤4)中，将基于G-L混合噪声特性核岭回归技术应用于风速预报中，构造预报模式为：输入向量为其中i,j为相关联的两个时刻；输出值为x_i+step，其中step为预测间隔时间。利用这种预报模式预报分析某一时刻i以后step间隔时刻的风速值。如取step＝3表示预测某一时刻i以后30分钟的风速值；取step＝12表示预测某一时刻i以后120分钟的风速值。

本发明的基于G-L混合噪声特性核岭回归技术的风速预报装置实施例

该装置包括损失函数获取模块、对偶问题求解模块、决策函数构造模块及风速预报模块；

其中，损失函数获取模块，用于获取某地区具有噪声影响的风速数据集D_l，利用Bayesian原理，得到G-L混合噪声特性的损失函数c(ξ)；

对偶问题求解模块，用于利用统计学习理论和最优化理论，结合步骤1)中得到的G-L混合噪声特性的损失函数，建立基于G-Le混合噪声特性核岭回归模型的原问题，推导并求解出基于G-L混合噪声特性核岭回归模型对偶问题；

决策函数构造模块，用于利用十折交叉验证技术确定基于G-L混合噪声特性核岭回归模型对偶问题的最优参数C、λ₁及λ₂，选取核函数K(·,·)；构造基于G-L混合噪声特性核岭回归模型的决策函数f(x)；

风速预报模块，用于构造基于G-L混合噪声特性核岭回归模型的风速预报模式：输入向量为其中i,j为相关联的两个时刻；输出值为x_i+step，其中step为预测间隔时间，利用该预报模式预报分析某一时刻i以后step间隔时刻的风速值。

该实施例中的G-L混合噪声特性的损失函数c(ξ)、基于G-L混合噪声特性的核岭回归模型的原问题、对偶问题、决策函数f(x)的具体公式及求解方法请参考方法实施例，这里不再一一详细阐述。

利用上述方法及装置进行风速预报时，表现了较好的性能，下面给出了风速预报性能评价的方法：

风速预报性能的评价一般用两个时间序列中基于预测值x_p和测量值x_m间的误差来度量，即ε_i＝x_p,i-x_m,i(i＝1,…,l)；

最常用的评价误差度量的指标是平均值绝对误差(the mean absolute error，简记为MAE)：

相对平均值绝对误差(the mean absolute percentage error，简记为MAPE)：

根平方值法(the root mean square error，简记为RMSE)、标准误差(the standard error of prediction,简记为SEP)是预测误差中应用比较广泛的方法，RMSE、SEP在两个时间序列中基于预测值x_p和测量值x_m定义为：

l为选择样本的尺度，x_p,i、x_m,i分别表示时间序列中第i个样本的预测值与测量值，为样本测量值的平均值。用平均值绝对误差、相对平均值绝对误差、根平方值误差、标准误差对ν-SVR、GN-KRR、GLM-KRR三种回归技术进行了评价。

在黑龙江省的风速数据集D_l中，其中的样本是每10分钟测量一次，共收集有62466个样本，D_l各列属性分别包括均值、方差、最小值、最大值等多个因子。我们取训练样本2880个(从1至2880，即20天的样本)，测试样本720个(从2161至2880，即5天的样本)进行了实验分析。输入向量为其中i,j为相关联的两个时刻；输出值为x_i+step，其中step＝1,3,6。即用向量分别预报某一时刻i以后10分钟、30分钟、60分钟的风速值，其中i＝1,…,720。

1、时刻i以后10分钟的风速预报结果

应用模型ν-SVR、GN-KRR、GLM-KRR预报某一时刻i以后10分钟以后的风速预报结果分别见图3～图5所示。

利用指标MAE、MAPE、RMSE和SEP评价三种模型10分钟以后的风速预报结果如表1所示。

表1：三种模型10分钟以后的风速预报的误差统计(测试样本720)

2、时刻i以后30分钟的风速预报结果

应用模型ν-SVR、GN-KRR、GLM-KRR预报某一时刻i以后30分钟以后的风速预报结果分别见图6～图8所示。

利用指标MAE、MAPE、RMSE和SEP评价三种模型30分钟以后的风速预报结果如表2所示。

表2：三种模型30分钟以后的风速预报的误差统计(测试样本720)

3、时刻i以后60分钟的风速预报结果

应用模型ν-SVR、GN-KRR、GLM-KRR预报某一时刻i以后60分钟以后的风速预报结果分别见图9～图11所示。

利用指标MAE、MAPE、RMSE和SEP评价三种模型60分钟以后的风速预报结果如表3所示。

表3：三种模型60分钟以后的风速预报的误差统计(测试样本720)

可见，应用模型ν-SVR、GN-KRR和GLM-KRR进行预报某一时刻i以后10分钟、30分钟和60分钟的风速预报的实验结果说明，模型GLM-KRR的预报结果比模型ν-SVR和GN-KRR的效果更好。