基于噪声模型支持向量回归技术的短期风速预报方法与流程

本发明涉及天气预报技术领域，具体讲，涉及基于噪声模型支持向量回归技术的短期风速预报方法。技术背景对于线性系统而言,从Gauss时代起,就利用最小二乘法把平面上的点拟合成直线,把高维空间的点拟合成超平面。经历了100多年的发展,经典最小二乘法已经成为许多领域数据处理的最广泛使用的方法。但是,对于线性回归中的不适定问题或非线性回归中的问题,基于最小二乘法的线性回归的性能可能变得很坏,针对这种情况,众多学者研究了最小二乘回归的改进问题,提出了许多新的回归算法。支持向量回归(Supportvectorregression,简记为SVR)就是其中之一。支持向量机方法是建立在统计学习理论和结构风险最小化原理基础上的，根据有限的样本信息在模型复杂性(即对特定训练样本的学习精度)和学习能力(即无错误地识别任意样本的能力)之间寻求最佳折衷，以期获得最好的推广能力。支持向量机方法的主要优点有：1.它是专门针对有限样本情况的，其目标是得到现有信息下的最优解而不仅仅是样本数趋于无穷大时的最优值；2.算法最终将转化成为一个二次型寻优问题，从理论上说，得到的将是全局最优点，解决了在神经网络方法中无法避免的局部极值问题；3.算法将实际问题通过非线性变换转换到高维的特征空间，在高维空间中构造线性决策函数来实现原空间中的非线性决策函数，特殊性质能保证机器有较好的推广能力，同时它巧妙地解决了维数问题，其算法复杂度与样本维数无关。支持向量回归方法自Vapnik于1995年提出以来,就得到了广泛的关注,它成功应用于科学技术和社会科学等各个领域。设给定数据:Dl＝{(x1,y1),...,(xi,yi),...,(xl,yl)}(1)其中xi∈Rn,yi∈R,i＝1,2,…,l,多元线性回归模型为f(xi)＝ωT·xi+b,xi＝(xi1,xi2,…,xin)T，参数向量ω∈Rn，其中xi∈Rn,Rn表示n维欧式空间,R表示实数集，l表示样本个数，上标T表示转置，在ν-SVR中，选取ε-不敏感损失函数线性ν-SVR的原问题为:其中C>0是惩罚因子，ξi,为松弛因子，(2)中目标函数假设样本不受噪声影响。(2)中的决策函数为线性回归模型。在实际应用中，样本Dl中xi与yi一般不满足线性关系，这样线性ν-SVR在解决非线性回归问题时往往不能取得预期的效果。通过核技巧构造合适的核变换Φ:Rn→H(H为Hilbert空间),利用Hilbert空间中的內积运算引进非线性核函数K(xi,xj)＝(Φ(xi)·Φ(xj)),把非线性问题转化为线性问题来解决。从而得到非线性ν-SVR的原问题为:不可避免地,样本都要受到噪声影响。一般地，假设样本受到高斯噪声影响。2010年Wu.提出基于高斯噪声影响的支持向量机，并利用遗传算法和粒子群优化算法进行了求解。其原问题为:当噪声服从Gauss分布时,基于Gauss噪声的支持向量回归(supportvectorregressionbasedontheGauss-noise,简记为GN-SVR)能够取得预期的效果。然而对于像风速预报的噪声、相关电磁波的冲击方向估计的噪声却不服从Gauss分布,它们分别服从Beta分布、Laplace分布，或其他分布。此时应用GN-SVR进行预测,则预报结果不能满足实际要求。

技术实现要素：
本发明旨在克服现有技术的不足，满足实际应用中(如风力发电、农业生产等)对短期风速预报的要求，为达到上述目的，本发明采取的技术方案是，应用Bayesian原理导出基于一般噪声模型的损失函数,给定具有噪声影响的数据集Dl＝{(x1,y1),...,(xi,yi),...,(xl,yl)},其中xi∈Rn,yi∈R,i＝1,...,l,Rn表示n维欧式空间,R表示实数集，l表示样本个数，求取最优损失函数；在此基础上构造基于噪声模型的支持向量回归机,最后利用基于噪声模型的支持向量回归技术得到短期风速预报；构造基于噪声模型的支持向量回归机具体为：(1)利用增广拉格朗日乘子法(AugmentedLagrangeMultipliermethod，简记为ALM)求解噪声支持向量回归模型，确定最优参数C、ν、m、n；选取合适的核函数K(·,·)；其中m、n为Beta噪声模型的损失函数中的参数，由噪声分布的期望μ和方差σ2确定,即m＝[(1-μ)·μ2/σ2]-μ,n＝[(1-μ)/μ]·m；(2)构造并求解最优化问题其中C>0是惩罚因子，0<ν<1是常数，s.t.为subjectto的缩写，DN-SVR表示基于噪声模型支持向量回归的对偶问题,表示基于噪声模型支持向量回归对偶问题的目标函数；得到最优解为拉格朗日乘子；(3)构造基于噪声模型支持向量回归的决策函数其中RSV为对应的样本,称为支持向量,Φ:Rn→H为核变换，H为Hilbert空间，K(xi,xj)＝(Φ(xi)·Φ(xj))，ω∈Rn为参数向量，(Φ(xi)·Φ(xj))表示H空间中的内积。求取最优损失函数：利用数据集Dl估计函数f(x),应用Bayesian原理的方法来得到最优损失函数c(x,y,f(x))＝-logp(y-f(x))，(7)其中p(y-f(x))＝p(ξ)表示噪声ξ的概率密度函数，c(xi,yi,f(xi))＝c(ξi)表示在样本点(xi,yi)进行预测时所得到预测值f(xi)与测量值yi比较所产生的损失，c(ξ)表示损失函数；得Gauss噪声模型的损失函数为:得Beta噪声模型的损失函数为:c(ξi)＝c(yi-f(xi))＝(1-m)logξi+(1-n)log(1-ξi)(9)。确定最优参数C、ν、m、n具体为：利用ALM法求解基于噪声模型的支持向量回归对偶问题式(5)及利用十折交叉验证方法确定最优参数C、ν、m、n，提出的基于噪声模型的支持向量回归机应用Matlab7.1程序语言实现,取N-SVR的参数集C∈[1,201],ν∈(0,1)，m,n∈(1,+∞)。利用核技巧构造合适的核函数K(·,·),把基于噪声模型的线性支持向量回归机推广为基于噪声模型的非线性支持向量回归机；其中K(xi,xj)＝(Φ(xi)·Φ(xj)),Φ:Rn→H,H为Hilbert空间,(Φ(xi)·Φ(xj))为H空间中的内积:(1)多项式核函数：K(xi,xj)＝((xi·xj)+1)d,(2)Gauss径向基核函数：K(xi,xj)＝exp(-||xi-xj||2/σ2),其中d是正整数,如取d＝2或3；σ是正数。构造并求解最优化问题进一步具体为：基于噪声模型支持向量回归的原问题：其中ξi＝yi-ωT·Φ(xi)-b,式中的PN-SVR表示基于噪声模型支持向量回归的原问题，表示基于噪声模型支持向量回归原问题的目标函数；通过构造Lagrange泛函可得到基于噪声模型支持向量回归原问题式(10)的对偶问题(简记为N-SVR)为:其中C>0为常数，是惩罚因子，0<ν<1是常数,包含ξi,包含αi,γ、ηi、为引进的辅助变量。在构造决策函数后，将基于噪声模型支持向量回归技术应用于短期风速预报中，构造模式为：输入向量为输出值为xi+6+step，其中step为预测间隔时间,利用这种模式预报分析某一时刻i+6以后step时刻的短期风速。本发明具备下列技术效果：本发明提出一种应用基于噪声模型的支持向量回归技术进行短期风速预报分析的新的统一理论模型及技术方法。利用最优化理论和一般噪声最优损失函数(7)构造的具有较高稳定性和鲁棒性的一类支持向量回归机,能够满足实际应用中(如风力发电、农业生产等)对短期风速预报的要求；在短期风速预报中,基于Beta噪声模型的支持向量回归技术的预报效果相比基于不考虑噪声模型的ν-SVR，和考虑Gauss噪声模型的GN-SVR和HN-SVR技术而言，表现出了更好的性能。附图说明图1Betapdf和Gausspdf。图2Lossfunctionofsixnoisemodels。图3预报结果ν-SVR(C＝81,ν＝0.5,step＝1)。图4预报结果GN-SVR(C＝81,ν＝0.5,step＝1)。图5预报结果HN-SVR(C＝81,ν＝0.5,step＝1)。图6预报结果BN-SVR(C＝81,ν＝0.5,m＝1.41,n＝1.71,step＝1)。图7预报结果ν-SVR(C＝81,ν＝0.5,step＝3)。图8预报结果GN-SVR(C＝81,ν＝0.5,step＝3)。图9预报结果HN-SVR(C＝81,ν＝0.5,step＝3)。图10预报结果BN-SVR(C＝81,ν＝0.5,m＝1.41,n＝1.71,step＝3)。图11预报结果ν-SVR(C＝81,ν＝0.5,step＝6)。图12预报结果GN-SVR(C＝81,ν＝0.5,step＝6)。图13预报结果HN-SVR(C＝81,ν＝0.5,step＝6)。图14预报结果BN-SVR(C＝81,ν＝0.5,m＝1.41,n＝1.71,step＝6)。图15本发明总体流程图。具体实施方式本发明基于噪声模型支持向量回归(Supportvectorregressionbasedonthenoise，简记为N-SVR)技术的短期风速预报方法，是以最优化理论和统计学习理论为基础,应用Bayesian原理的方法导出基于一般噪声模型的损失函数,在此基础上构造基于噪声影响的支持向量回归机的统一模型,最后利用基于噪声模型的支持向量回归技术得到短期风速预报。提出了基于噪声模型的支持向量回归技术，其原问题:其中ξi＝yi-ωT·Φ(xi)-b。可证得基于噪声模型支持向量回归原问题(12)关于ω的解存在且唯一。通过构造Lagrange泛函L(ω,b,α(*),ξ(*))可得到基于噪声模型支持向量回归原问题(12)的对偶问题(简记为N-SVR)为:其中C>0是惩罚因子,0<ν<1是常数,且有则可得到基于噪声模型支持向量回归的决策函数为:其中RSV为对应的样本，称为支持向量。Φ:Rn→H(H为Hilbert空间)为核变换，K(xi,xj)＝(Φ(xi)·Φ(xj))，ω∈Rn为参数向量，(Φ(xi)·Φ(xj))表示H空间中的内积。本发明目的是应用基于噪声模型的支持向量回归技术进行短期风速预报。由于各种原因,现实世界中的数据都受到噪声影响,因此基于噪声数据的机器学习模型具有重要的应用价值。当噪声服从Gauss分布时,基于Gauss噪声的支持向量回归能够取得预期的效果。然而对于像风速预报的噪声、气象预报(如雨量、温度等)的噪声却不服从Gauss分布。S.Bofinger等发现风功率预报的噪声并不服从Gauss分布,通过实验结果显示风功率预报的噪声服从Beta分布。A.Fabbri等研究了当预测值xp和测量值xm间的误差ε服从Beta分布时,ε的概率密度函数(probabilitydistributionfunction，简记为pdf)为f(x)＝xm-1·(1-x)n-1·h(图1),其中m,n是参数,由噪声分布的期望和方差确定,m＝(1-μ)·μ2/σ2-μ,n＝1-μ/μ·m,h＝Γ(m+n)/Γ(m)·Γ(n)是归一化因子。给定具有噪声影响的数据Dl,决策函数f(x)是未知的。一般地，最小化目标函数：其中λ是正数,c(ξi)＝c(yi-f(xi))是损失函数。假设噪声是加性的,即yi＝fi(xi)+ξi(i＝1,…,l)，且ξi是独立同分布(i.i.d.)的随机变量，ξi的方差为σ2，均值为μ。利用数据g∈Dl估计函数f(x),最大化后验概率P[f|g]来求得最优损失函数。由Bayesian原理得：P[f|g]∝P[g|f]·P[f]其中P[g|f]为条件概率，P[f]为先验概率，且假设P[f]∝exp(-λ·||f||2)，exp(-λ·||f||2)为光滑函数。由于噪声是加性的，且ξi是独立同分布的，故：由式(14)和(15)得最优损失函数为：c(x,y,f(x))＝-logp(y-f(x))(16)从而得Gauss噪声模型的损失函数为(图2):得Laplace噪声模型的损失函数为(图2):c(ξi)＝c(yi-f(xi))＝|ξi|(18)得Beta噪声模型的损失函数为(图2):c(ξi)＝c(yi-f(xi))＝(1-m)logξi+(1-n)log(1-ξi)(19)当噪声服从Beta分布或Laplace分布时,应用GN-SVR进行预测,则预报结果不能满足实际要求。为了解决上述问题,我们利用最优化理论和选取噪声损失函数(16),提出了基于一般噪声的支持向量回归模型(N-SVR),并成功应用于短期风速预报中,取得了良好的效果。本发明提出一种应用基于一般噪声的支持向量回归模型(N-SVR)和预报技术进行短期风速预报分析的一般理论。利用最优化理论和噪声损失函数构造的具有较高稳定性和鲁棒性的一种支持向量回归机,能够满足实际应用中(如风力发电、农业生产等)对短期风速预报的要求。应用基于噪声模型的支持向量回归技术进行短期风速预报主要通过以下几个步骤来实现：1设给定具有噪声影响的数据集Dl＝{(x1,y1),...,(xi,yi),...,(xl,yl)},其中xi∈Rn,yi∈R,i＝1,2,…,l，Rn表示n维欧式空间，R表示实数集，l表示样本个数；求取最优损失函数；2.利用ALM法求解噪声模型支持向量回归问题。确定最优参数C、ν、m、n；选取合适的核函数K(·,·)，其中m、n为Beta噪声模型的损失函数中的参数，由噪声分布的期望μ和方差σ2确定；3构造并求解最优化问题其中C>0是惩罚因子，0<ν<1是常数。得到最优解为拉格朗日乘子；4构造基于噪声模型支持向量回归的决策函数其中(RSV为对应的样本，称为支持向量)，Φ:Rn→H为核变换，H为Hilbert空间，K(xi,xj)＝(Φ(xi)·Φ(xj))，ω∈Rn为参数向量，(Φ(xi)·Φ(xj))表示H空间中的内积。在短期风速预报中,基于Beta噪声模型的支持向量回归技术的预报效果相比基于不考虑噪声模型的ν-SVR，和考虑Gauss噪声模型的GN-SVR和HN-SVR技术而言，表现出了更好的性能。风速预报性能的评价,一般用两个时间序列中基于预测值xp和测量值xm间的相似度来度量,即εi＝xp,i-xm,i,i＝1,…,l最常用的评价误差度量的指标是平均值绝对误差(themeanabsoluteerror，简记为MAE):相对平均值绝对误差(themeanabsolutepercentageerror，简记为MAPE):根平方值法(therootmeansquareerror，简记为RMSE)是预测误差中应用比较广泛的方法,RMSE在两个时间序列中基于预测值xp和测量值xm定义为:xp,i、xm,i分别表示时间序列中第i个预测值与测量值。用平均值绝对误差、相对平均值绝对误差、根平方值误差对ν-SVR、GN-SVR、HN-SVR、BN-SVR等四种回归技术模型进行了评价。在黑龙江省的风速数据集Dl中，其中的样本是每10分钟测量一次,样本容量为62466,Dl各列属性分别包括均值、方差、最小值、最大值等多个因子。我们取训练样本432个(从1至432,即3天的样本)，测试样本432个(从433至864,即3天的样本)进行了实验分析。输入向量为,输出值为xi+6+step,其中step＝1,3,6。即用向量分别预报某一时刻i+6以后10分钟、30分钟、60分钟的风速。1.时刻i+6以后10分钟的风速预报结果应用ν-SVR、GN-SVR、HN-SVR和BN-SVR进行预报某一时刻i+6以后10分钟的短期风速预报结果如图3-6所示。利用指标MAE、MAPE和RMSE评价四种模型的预报结果如表1.表1:四种模型的短期风速预报的误差统计(测试样本432)2.时刻i+6以后30分钟的风速预报结果应用ν-SVR、GN-SVR、HN-SVR和BN-SVR进行预报某一时刻i+6以后30分钟的短期风速预报结果如图7-10所示。利用指标MAE、MAPE和RMSE评价四种模型的预报结果如表2.表2:四种模型的短期风速预报的误差统计(测试样本432)3.时刻i+6以后60分钟的风速预报结果应用ν-SVR、GN-SVR、HN-SVR和BN-SVR进行预报某一时刻i+6以后60分钟的短期风速预报结果如图11-14所示。利用指标MAE、MAPE和RMSE评价四种模型的预报结果如表3.表3:四种模型的短期风速预报的误差统计(测试样本432)结论:应用ν-SVR、GN-SVR、HN-SVR和BN-SVR进行预报某一时刻i+6以后10分钟、30分钟、60分钟的短期风速预报的实验结果说明,BN-SVR的预报结果比ν-SVR、GN-SVR、HN-SVR的效果更好。下面结合附图和具体实施方式进一步详细说明本发明。1.计算噪声样本的损失函数设给定具有噪声影响的数据集Dl＝{(x1,y1),...,(xi,yi),...,(xl,yl)},其中xi∈Rn,yi∈R,i＝1,2,...，l.利用数据集Dl估计函数f(x),应用Bayesian原理的方法来得到最优损失函数:c(x,y,f(x))＝-logp(y-f(x))(25)2.利用增广拉格朗日乘子法(AugmentedLagrangeMultipliermethod，简记为ALM)求解噪声支持向量回归模型，确定最优参数C、ν、m、n。拉格朗日乘子法是Powel和Hestenes于1969年针对等式约束优化问题同时独立提出的优化算法，也称为PH算法。其基本思想是：从原问题的拉格朗日函数出发，再加上适当的罚函数，从而将原问题转化为求解一系列的无约束优化子问题。增广拉格朗日乘子法(AugmentedLagrangeMultipliermethod，简记为ALM)是1973年Rockfellar将PH算法推广到求解不等式约束优化问题，也称为PHR算法。ALM法是解决同时带有等式和不等式约束问题的一类优化方法。ALM法把解等式约束优化问题的乘子法推广到不等式约束优化问题，即先引进辅助变量把不等式约束转化为等式约束，然后再利用最有性条件消去辅助变量。对于一个求函数最小值的优化问题(求函数最大值也类似)，一般可以描述为下列数学规划模型:式中x为决策变量,f(x)为目标函数,式为约束条件,U是基本空间,R是U的子集。满足约束条件的解X称为可行解,集合R表示所有满足约束条件的解所组成的集合,称为可行解集合。式(16)、(25)c(x,y,f(x))＝-logp(y-f(x))中的x与式(26)中的x的含义相同，x＝(x1,x2,…,xl)T，y＝(y1,y2,…,yl)T，(xi,yi)∈Dl,i＝1,2,…,l，上标T表示转置。p(y-f(x))＝p(ξ)表示噪声ξ的概率密度函数。c(xi,yi,f(xi))＝c(ξi)表示在样本点(xi,yi)进行预测时所得到预测值f(xi)与yi比较所产生的损失，c(ξ)表示损失函数。利用十折交叉验证方法确定基于噪声模型支持向量回归机的最优参数C、ν、m、n。提出的基于噪声模型的支持向量回归机应用Matlab7.1程序语言实现,取N-SVR的参数集C∈[1,201],ν∈(0,1),m,n∈(1,+∞)。3.选取合适的核函数K(·,·)。利用核技巧构造核函数K(·,·),把基于噪声模型的线性支持向量回归机推广为噪声模型的非线性支持向量回归机。其中K(xi,xj)＝(Φ(xi)·Φ(xj)),Φ:Rn→H,H为Hilbert空间,(Φ(xi)·Φ(xj))为H空间中的内积。多项式核函数和Gauss径向基核函数在诸多实际应用中具有较强的鲁棒性和稳定的性能,且容易操作实施。(1)多项式核函数：K(xi,xj)＝((xi·xj)+1)d是正整数,(2)Gauss径向基核函数：K(xi,xj)＝exp(-||xi-xj||2/σ2),其中d是正数,取d＝2或3；σ是正数。4.构造并求解最优化问题。基于噪声模型支持向量回归的原问题为:其中ξi＝yi-ωT·Φ(xi)-b,i＝1，2,...,l，式中的PN-SVR表示基于噪声的支持向量回归的原问题，表示基于噪声模型支持向量回归原问题的目标函数；通过构造Lagrange泛函可得到基于噪声模型支持向量回归原问题(27)的对偶问题(简记为N-SVR)为:其中C>0是惩罚因子，0<ν<1是常数。式(5)、(11)、(13)、(20)、(28)等中的DN-SVR表示基于噪声模型支持向量回归的对偶问题,表示基于噪声模型支持向量回归对偶问题的目标函数。式(5)、(6)、(11)、(13)、(20)、(21)、(28)、(29)中，αi,为拉格朗日乘子，即为对偶问题(28)中的自变量。5.构造基于噪声模型支持向量回归的决策函数其中(RSV为对应的样本，称为支持向量)。6.将基于噪声模型的支持向量回归技术应用于短期风速预报中。构造预报模式为：输入向量为输出值为xi+6+step，其中step为预测间隔时间。利用这种预报模式预报分析某一时刻i+6以后step时刻的短期风速。如取step＝1表示预测某一时刻i+6以后10分钟的风速；取step＝6表示预测某一时刻i+6以后60分钟的风速。外文字符的中文注释：1.Beta：贝塔(希腊字母)。2.Gauss：高斯(著名数学家)。3.Supportvectorregression：支持向量回归。4.Vapnik，Wu,A.Fabbri，S.Bofinger，J.Holland：人名。5.Hilbert：希尔伯特(著名数学家)。6.Supportvectorregressionbasedonthenoise:基于噪声的支持向量回归,简记为N-SVR。7.SupportvectorregressionbasedontheGauss-noise:基于高斯-噪声的支持向量回归,简记为GN-SVR。8.SupportvectorregressionbasedontheBeta-noise:基于Beta-噪声的支持向量回归,简记为BN-SVR。9.C、ν、m、n:C、ν是ν-SVR、GN-SVR、HN-SVR、BN-SVR模型中的参数，C>0是惩罚因子，0<ν<1是常数；m、n是BN-SVR模型中的参数，m,n∈(1,+∞)。10.Lagrange:拉格朗日(著名数学家)。11.Powel，Hestenes，Rockfellar：人名。12.probabilitydistributionfunction(pdf):概率密度函数。13.themeanabsoluteerror(MAE)：平均值绝对误差。14.themeanabsolutepercentageerror(MAPE)：相对平均值绝对误差。15.therootmeansquareerror(RMSE)：根平方值法。16.输入向量，i表示某一时刻；xi+6：表示某一时刻i以后60分钟时刻的风速。17.xi+6+step，step＝1,3,6：表示某一时刻i+6以后10分钟、30分钟、60分钟时刻的风速。18.AugmentedLagrangeMultipliermethod:增广拉格朗日乘子法,简记为ALM。19.minf(x)：求函数f(x)的最小值；maxf(x)：求函数f(x)的最大值。20.PN-SVR:表示基于噪声模型支持向量回归的原问题；表示基于噪声模型支持向量回归原问题的目标函数。21.DN-SVR:表示基于噪声模型支持向量回归的对偶问题；表示基于噪声模型支持向量回归对偶问题的目标函数。22.Φ:表示核变换Φ:Rn→H,Rn为n维欧式空间,H为Hilbert空间；核函数K(xi,xj)＝(Φ(xi)·Φ(xj)),(Φ(xi)·Φ(xj))为H空间中的内积。23.ξ:表示噪声函数,ξ＝y-f(x),其中f(x)＝ωT·Φ(x)-b,ω∈Rn为参数向量。24.ξi:表示样本点(xi,yi)∈Dl的噪声,ξi＝yi-f(xi),i＝1,2,...,l′。25.p(ξ)表示噪声ξ的概率密度函数；c(ξi)表示在样本点(xi,yi)进行预测时所得到预测值f(xi)与测量值yi比较所产生的损失,c(ξ)表示噪声ξ的损失函数。26.为的函数,其中包含ξi,包含αi,27.c(x,y,f(x)):c(x,y,f(x))＝-logp(y-f(x)),(x,y)∈Dl为损失函数.当点(x,y)取Dl中某固定点(xi,yi),即x＝xi,y＝yi时,c(x,y,f(x))在点(xi,yi)的损失为c(xi,yi,f(xi))＝-logp(yi-f(xi))。28.L(ω,b,α(*),ξ(*)):为Lagrange泛函,其中γ、ηi、为引进的辅助变量。29.ωT,ω表示参数向量,T表示转置。