CN102542167B - 一种风电场风速时间序列预测方法 - Google Patents

Abstract

本发明是一种风电场风速时间序列预测方法，其特点是，包括的步骤有：利用风速采集仪器每隔一小时记录一次同一地区的风速数据，整理采集的原始风速数据形成风速时间序列用于分析预测；运用快速独立分量分析算法对风速时间序列进行多尺度分解，分解得到多个独立分量；通过对各个独立分量的延迟时间和嵌入维数的计算，运用相空间重构理论对各个独立分量进行相空间重构；利用最小二乘支持向量机回归模型对相空间重构后的各个独立分量进行建模预测；将预测结果进行叠加得到最终的风速时间序列的预测结果。具有风速时间序列预测科学合理，准确、可靠和适用性强等优点。

Description

一种风电场风速时间序列预测方法

技术领域

本发明涉及一种风电场风速时间序列预测方法，尤其涉及一种基于局域均值分解和多核最小二乘支持向量机的短期风速预测方法。 

背景技术

能源与环境是当今人类生存和发展所急需解决的问题。常规能源以煤、石油、天然气为主，它不仅资源有限，而且还造成了严重的大气污染。因此，对可再生能源的利用，尤其是对风能的开发利用，已受到各个国家的高度重视。随着风能利用的加速发展，越来越多的大型风电场将纳入电网的统调，风电在电网的比重越来越大，但是由于系统的最大负荷受限于风电场穿透功率极限，所以当负荷超过一定值，就会严重影响电网稳定的运行。而对风速的准确预测可以减少风电场的旋转设备和运行成本，提高风电穿透功率极限，可以帮助调度部门及时调整计划，从而减轻风能对电网的冲击。 

目前，国内外用于风速预测的方法主要有持续预测法、卡尔曼滤波法、时间序列分析法和神经网络方法。以上预测方法大多都针对原始风速序列，存在着风速预测精度低，平均相对预测误差通常达到20％。 

发明内容

本发明的目的在于克服现有风速预测技术上的缺陷，提出一种准确、可靠，应用价值高的基于局域均值分解和多核最小二乘支持向量机的短期风速预测方法。 

本发明的目的是由以下技术方案来实现的：一种风电场风速时间序列预测方法，其特征在于，它包括下述步骤： 

(1)利用风速采集仪器每隔一小时记录一次同一地区的风速数据，整理采集的原始风速数据，形成风速时间序列用于分析预测； 

(2)运用局域均值分解算法对风速时间序列进行多尺度分解，分解得到多个PF分量， 

局域均值分解算法步骤如下： 

(a)找出原始信号x(t)所有的局部极值点n_i，求出所有相邻的局部极值点的平均值： 

$m_i=\frac{n_i+n_{i+1}}{2}---\left(1\right)$

将所有相邻的平均值点m_i用直线连接起来，然后用滑动平均法进行平滑处理，得到局 部均值函数m₁₁(t)； 

(b)求出包络估计值 

$a_i=\frac{|n_i-n_{i+1}|}{2}---\left(2\right)$

将所有相邻的平均值点a_i用直线连接起来，然后用滑动平均法进行平滑处理，得到包络估计函数a₁₁(t)； 

(c)将局部均值函数m₁₁(t)从原始信号x(t)中分离出来，得到 

h₁₁(t)＝x(t)-m₁₁(t)    (3) 

(d)用h₁₁(t)除以包络估计函数a₁₁(t)以对h₁₁(t)进行解调，得到 

s₁₁(t)＝h₁₁(t)/a₁₁(t)    (4) 

对s₁₁(t)重复上述步骤便能得到s₁₁(t)的包络估计函数a₁₂(t)，假如a₁₂(t)不等于1，说明s₁₁(t)不是一个纯调频信号，需要重复上述迭代过程n次，直至s_1n(t)为一个纯调频信号，也就是s_1n(t)的包络估计函数a_1(n+1)(t)＝1，所以有 

$\left\{\begin{array}{c}{h}_{11}\left(t\right)=x\left(t\right)-m_{11}\left(t\right)\\ h_{12}\left(t\right)=s_{11}\left(t\right)-m_{12}\left(t\right)\\ M\\ h_{1n}\left(t\right)=s_{1(n-1)}\left(t\right)-m_{1n}\left(t\right)\end{array}\right.---\left(5\right)$

式中 

$\left\{\begin{array}{c}{s}_{11}\left(t\right)=h_{11}\left(t\right)/a_{11}\left(t\right)\\ s_{12}\left(t\right)=h_{12}\left(t\right)/a_{12}\left(t\right)\\ M\\ s_{1n}\left(t\right)=h_{1n}\left(t\right)/a_{1n}\left(t\right)\end{array}\right.---\left(6\right)$

迭代终止条件为 

$\underset{n\→\∞}{\lim }{a}_{1n}\left(t\right)=1---\left(7\right)$

为了减少迭代次数，降低运算时间，设置一个变量Δ＝10^-4，使得当满足1-Δ≤a_1n(t)≤1+Δ时，终止迭代； 

(e)将迭代过程中所产生的全部包络估计函数相乘，得到包络信号瞬时幅值函数： 

$a_1\left(t\right)=a_{11}\left(t\right)a_{12}\left(t\right){\mathrm{\Λa}}_{1n}\left(t\right)=\underset{q=1}{\overset{n}{\mathrm{\Π}}}{a}_{1q}\left(t\right)---\left(8\right)$

(f)将包络信号a₁(t)和纯调频信号s_1n(t)相乘得： 

PF₁(t)＝a₁(t)s_1n(t)    (9) 

它为信号x(t)的第一个PF分量，包含了原始信号的最高频率成分，是一个单分量的调幅、调频信号，瞬时频率f₁(t)则由调频信号s_1n(t)求出： 

$f_1\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{d\left[\arccos \left(s_{1n}\left(t\right)\right)\right]}{\mathrm{dt}}---\left(10\right)$

(g)从原始信号x(t)中将第一个PF分量PF₁(t)分离出来，得到剩余信号u₁(t)，由于剩余信号u₁(t)中还包含有较多的频率成分，因此将u₁(t)作为原始数据重复以上步骤对其进行分解，得到第二个PF分量，重复这个过程直到u_k为一个单调函数为止，得到一定数量的PF分量： 

$\left\{\begin{array}{c}{u}_1\left(t\right)=x\left(t\right)-{\mathrm{PF}}_1\left(t\right)\\ u_2\left(t\right)=u_1\left(t\right)-{\mathrm{PF}}_2\left(t\right)\\ M\\ u_k\left(t\right)=u_{k-1}\left(t\right)-{\mathrm{PF}}_k\left(t\right)\end{array}\right.---\left(11\right)$

最终信号x(t)可表示为k个PF分量和余量之和： 

$x\left(t\right)=\underset{p=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{PF}}_p\left(t\right)+u_k\left(t\right)---\left(12\right)$

(3)通过对各个PF分量的延迟时间τ和嵌入维数m的计算，运用相空间重构理论对各个PF分量进行相空间重构； 

y_PF(t)＝{PF_p(t)，PF_p(t+τ)，Λ，PF_p[t+(m-1)τ]}    (13) 

y_u(t)＝{u_k(t)，u_k(t+τ)，Λ，u_k[t+(m-1)τ]}        (14) 

t＝1，2，Λ.，p＝1，2，Λ，k 

(4)利用多核最小二乘支持向量机回归模型对相空间重构后的各个PF分量进行建模预测； 

$y_{\mathrm{PFl}}(t+\mathrm{m\τ})=\underset{i=t}{\overset{t+m-1}{\mathrm{\Σ}}}{\mathrm{\α}}_{i}K(y_{\mathrm{PFi}},x)+b,l=\mathrm{1,2},\mathrm{\Λ},k---\left(15\right)$

$y_{\mathrm{uk}}(t+\mathrm{m\τ})=\underset{i=t}{\overset{t+m-1}{\mathrm{\Σ}}}{\mathrm{\α}}_{i}K(y_{\mathrm{ui}},x)+b---\left(16\right)$

(5)将预测结果进行叠加得到最终的风速时间序列的预测结果 

$\hat{y}(t+\mathrm{\τ})=\underset{l=1}{\overset{k}{\mathrm{\Σ}}}{y}_{\mathrm{PFl}}(t+\mathrm{\τ})+y_{\mathrm{uk}}(t+\mathrm{\τ})---\left(17\right)$

所述原始风速数据中相邻两个风速数据点之间的采集间隔为15～60min。 

所述步骤(3)中，所述各个PF分量(p＝1，2，Λ，k)的延迟时间使用自相关法获得，所述各个PF分量的嵌入维数使用假近邻法获得，自相关法计算延迟时间τ的方法为：按照式(18)计算各PF_p分量跨度为jτ的自相关函数R(jτ) 

$R\left(\mathrm{j\τ}\right)=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}[{\mathrm{PF}}_p\left(i\right)\·{\mathrm{PF}}_p(i+\mathrm{j\τ})]---\left(18\right)$

式中PF_pi(i＝1，2，Λ，N)为PF分量，固定j，作出自相关函数R(jτ)关于时间τ(即取τ＝1，2，Λ)的函数图像，则R(jτ)第一次达到零点时所对应的τ即为延迟时间； 

假近邻法确定嵌入维数方法具体为：在m维空间中，每个相点矢量为X(i)＝{PF_p(i)，PF_p(i+τ)，Λ，PF_p(i+(m-1)τ)}，都有一个某距离内的最近邻点X^NN(i)，其距离为R_m(i)，则 

$R_m^2\left(i\right)={\left|\right|X\left(i\right)-X^{\mathrm{NN}}\left(i\right)\left|\right|}^2=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\left\{{\mathrm{PF}}_p\right[(i+(j-1)\mathrm{\τ}]-{{\mathrm{PF}}_p}^{\mathrm{NN}}[i+(j-1)\mathrm{\τ}\left]\right\}}^2---\left(19\right)$

当相空间的维数增加到m+1维时，这两个相点的距离就会发生变化，两者的距离成为R_m+1(i)，且有 

$R_{m+1}^2\left(i\right)=\underset{j=1}{\overset{m+1}{\mathrm{\Σ}}}{\left\{{\mathrm{PF}}_p\right[(i+(j-1)\mathrm{\τ}]-{{\mathrm{PF}}_p}^{\mathrm{NN}}[i+(j-1)\mathrm{\τ}\left]\right\}}^2=R_m^2\left(i\right)+{|{\mathrm{PF}}_p(i+\mathrm{m\τ})-{{\mathrm{PF}}_p}^{\mathrm{NN}}(i+\mathrm{m\τ})|}^2---\left(20\right)$

由m增加1而引起的两近邻点间距离的变化时 

${[R_{m+1}^2\left(i\right)-R_m^2\left(i\right)]}^{1/2}=|{\mathrm{PF}}_p(i+\mathrm{m\τ})-{{\mathrm{PF}}_p}^{\mathrm{NN}}(i+\mathrm{m\τ})|---\left(21\right)$

定义新距离E(m) 

$E\left(m\right)=\frac{1}{N-\mathrm{m\τ}}\underset{i=1}{\overset{N-\mathrm{m\τ}}{\mathrm{\Σ}}}|{\mathrm{PF}}_p(i+\mathrm{m\τ})-{{\mathrm{PF}}_p}^{\mathrm{NN}}(i+\mathrm{m\τ})|---\left(22\right)$

E₁(m)＝E(m+1)/E(m)    (23) 

式中：E₁(m)为两个新距离的比值，当E₁(m)在m＞m₀时，E₁(m)＝1不再变化，则m₀即为确定的嵌入维数。 

所述步骤(4)中，所述的多核最小二乘支持向量机回归模型中的多核核函数，表示为 

K＝λ₁K_poly+λ₂K_RBF+λ₃K_exp

式中，K为多核函数，K_ploy为多项式核函数，K_poly＝(x_ix_j)+c)^d；K_RBF为径向基函数，K_RBF＝exp(-(x_i-x_j)²/(2σ²))；K_exp为指数核函数，K_exp＝exp(-||x_i-x_j||/(2σ′²))；λ₁，λ₂，λ₃分别为核函数的权值， 且0≤λ_k≤1。 

本发明的一种风电场风速时间序列预测方法的优点体现在： 

1.由于运用局域均值分解算法，对原始风速时间序列进行分解，得到多个不同尺度的PF分量，从而简化了不同特征信息间的藕合，降低了预测建模的难度； 

2.由于采用最小二乘支持向量机回归模型进行预测，通过对多个不同类型核函数的线性加权组合构造新的等价核函数，降低单一核函数及其参数的选择对建模精度的影响； 

3.由于将局域均值分解和多核最小二乘支持向量机二者相结合来建立预测模型，充分发挥两种算法的优点，从而提高了预测的科学性，准确、可靠性和适用性。 

附图说明

图1是本发明的一种风电场风速时间序列预测方法流程图。 

图2是本发明具体实施方式中东北某风电场2010年12月份的实测风速时间数列。 

图3是本发明中具体实施方式中利用局域均值分解将原始风速时间序列进行多尺度分解的分解结果。 

图4是本发明具体实施方式中的风速的预测效果图。 

具体实施方式

参照图1，本发明的一种风电场风速时间序列预测方法是基于局域均值分解Local Mean Decomposition，LMD和多核最小二乘支持向量机Multiple Kernel Least Square Support Vector Machine，MK-LSSVM的短期风速预测方法，包括下述步骤： 

(1)利用风速采集仪器每隔一小时记录一次同一地区的风速数据，整理采集的原始风速数据，形成风速时间序列用于分析预测； 

(2)运用LMD算法对风速时间序列进行多尺度分解，分解得到多个PF分量。 

LMD算法步骤如下： 

(a)找出原始信号x(t)所有的局部极值点n_i，求出所有相邻的局部极值点的平均值： 

$m_i=\frac{n_i+n_{i+1}}{2}---\left(1\right)$

将所有相邻的平均值点m_i用直线连接起来，然后用滑动平均法进行平滑处理，得到局部均值函数m₁₁(t)。 

(b)求出包络估计值 

$a_i=\frac{|n_i-n_{i+1}|}{2}---\left(2\right)$

将所有相邻的平均值点a_i用直线连接起来，然后用滑动平均法进行平滑处理，得到包络估计函数a₁₁(t)。 

(c)将局部均值函数m₁₁(t)从原始信号x(t)中分离出来，得到 

h₁₁(t)＝x(t)-m₁₁(t)    (3) 

(d)用h₁₁(t)除以包络估计函数a₁₁(t)以对h₁₁(t)进行解调，得到 

s₁₁(t)＝h₁₁(t)/a₁₁(t)    (4) 

对s₁₁(t)重复上述步骤便能得到s₁₁(t)的包络估计函数a₁₂(t)，假如a₁₂(t)不等于1，说明s₁₁(t)不是一个纯调频信号，需要重复上述迭代过程n次，直至s_1n(t)为一个纯调频信号，也就是s_1n(t)的包络估计函数a_1(n+1)(t)＝1，所以有 

$\left\{\begin{array}{c}{h}_{11}\left(t\right)=x\left(t\right)-m_{11}\left(t\right)\\ h_{12}\left(t\right)=s_{11}\left(t\right)-m_{12}\left(t\right)\\ M\\ h_{1n}\left(t\right)=s_{1(n-1)}\left(t\right)-m_{1n}\left(t\right)\end{array}\right.---\left(5\right)$

式中 

$\left\{\begin{array}{c}{s}_{11}\left(t\right)=h_{11}\left(t\right)/a_{11}\left(t\right)\\ s_{12}\left(t\right)=h_{12}\left(t\right)/a_{12}\left(t\right)\\ M\\ s_{1n}\left(t\right)=h_{1n}\left(t\right)/a_{1n}\left(t\right)\end{array}\right.---\left(6\right)$

迭代终止条件为 

$\underset{n\→\∞}{\lim }{a}_{1n}\left(t\right)=1---\left(7\right)$

为了减少迭代次数，降低运算时间，可以设置一个变量Δ＝10^-4，使得当满足 1-Δ≤a_1n(t)≤1+Δ时，终止迭代。 

(e)将迭代过程中所产生的全部包络估计函数相乘，得到包络信号(瞬时幅值函数)： 

$a_1\left(t\right)=a_{11}\left(t\right)a_{12}\left(t\right){\mathrm{\Λa}}_{1n}\left(t\right)=\underset{q=1}{\overset{n}{\mathrm{\Π}}}{a}_{1q}\left(t\right)---\left(8\right)$

(f)将包络信号a₁(t)和纯调频信号s_1n(t)相乘得： 

PF₁(t)＝a₁(t)s_1n(t)    (9) 

它为信号x(t)的第一个PF分量，包含了原始信号的最高频率成分，是一个单分量的调幅-调频信号，瞬时频率f₁(t)则可由调频信号s_1n(t)求出： 

$f_1\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{d\left[\arccos \left(s_{1n}\left(t\right)\right)\right]}{\mathrm{dt}}---\left(10\right)$

(g)从原始信号x(t)中将第一个PF分量PF₁(t)分离出来，得到剩余信号u₁(t)，由于剩余信号u₁(t)中还包含有较多的频率成分，因此将u₁(t)作为原始数据重复以上步骤对其进行分解，得到第二个PF分量，重复这个过程直到u_k为一个单调函数为止，得到一定数量的PF分量： 

$\left\{\begin{array}{c}{u}_1\left(t\right)=x\left(t\right)-{\mathrm{PF}}_1\left(t\right)\\ u_2\left(t\right)=u_1\left(t\right)-{\mathrm{PF}}_2\left(t\right)\\ M\\ u_k\left(t\right)=u_{k-1}\left(t\right)-{\mathrm{PF}}_k\left(t\right)\end{array}\right.---\left(11\right)$

最终信号x(t)可表示为k个PF分量和余量之和： 

$x\left(t\right)=\underset{p=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{PF}}_p\left(t\right)+u_k\left(t\right)---\left(12\right)$

(3)通过对各个PF分量的延迟时间τ和嵌入维数m的计算，运用相空间重构理论对各个PF分量进行相空间重构； 

y_PF(t)＝{PF_p(t)，PF_p(t+τ)，Λ，PF_p[t+(m-1)τ]}    (13) 

y_u(t)＝{u_k(t)，u_k(t+τ)，Λ，u_k[t+(m-1)τ]}        (14) 

t＝1，2，Λ.，p＝1，2，Λ，k 

(4)利用多核最小二乘支持向量机回归模型对相空间重构后的各个PF分量进行建模预测； 

$y_{\mathrm{PFl}}(t+\mathrm{m\τ})=\underset{i=t}{\overset{t+m-1}{\mathrm{\Σ}}}{\mathrm{\α}}_{i}K(y_{\mathrm{PFi}},x)+b,l=\mathrm{1,2},\mathrm{\Λ},k---\left(15\right)$

$y_{\mathrm{uk}}(t+\mathrm{m\τ})=\underset{i=t}{\overset{t+m-1}{\mathrm{\Σ}}}{\mathrm{\α}}_{i}K(y_{\mathrm{ui}},x)+b---\left(16\right)$

(5)将预测结果进行叠加得到最终的风速时间序列的预测结果 

$\hat{y}(t+\mathrm{\τ})=\underset{l=1}{\overset{k}{\mathrm{\Σ}}}{y}_{\mathrm{PFl}}(t+\mathrm{\τ})+y_{\mathrm{uk}}(t+\mathrm{\τ})---\left(17\right)$

所述原始风速数据中相邻两个风速数据点之间的采集间隔为15～60min。 

所述步骤(3)中，所述各个PF分量(p＝1，2，Λ，k)的延迟时间使用自相关法获得，所述各个PF分量的嵌入维数使用假近邻法获得，自相关法计算延迟时间τ的方法为：按照式(18)计算各PF_p分量跨度为jτ的自相关函数R(jτ) 

$R\left(\mathrm{j\τ}\right)=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}({\mathrm{PF}}_p\left(i\right)\·{\mathrm{PF}}_p(i+\mathrm{j\τ}))---\left(18\right)$

式中PF_pi(i＝1，2，Λ，N)为PF分量。固定j，作出自相关函数R(jτ)关于时间τ(即取τ＝1，2，Λ)的函数图像，则R(jτ)第一次达到零点时所对应的τ即为延迟时间τ。 

假近邻法确定嵌入维数方法具体为：在m维空间中，每个相点矢量为X(i)＝{PF_p(i)，PF_p(i+τ)，Λ，PF_p(i+(m-1)τ)}，都有一个某距离内的最近邻点X^NN(i)，其距离为R_m(i)，则 

$R_m^2\left(i\right)={\left|\right|X\left(i\right)-X^{\mathrm{NN}}\left(i\right)\left|\right|}^2=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\left\{{\mathrm{PF}}_p\right[(i+(j-1)\mathrm{\τ}]-{{\mathrm{PF}}_p}^{\mathrm{NN}}[i+(j-1)\mathrm{\τ}\left]\right\}}^2---\left(19\right)$

但相空间的维数增加到m+1维时，这两个相点的距离就会发生变化，两者的距离成为R_m+1(i)，且有 

$R_{m+1}^2\left(i\right)=\underset{j=1}{\overset{m+1}{\mathrm{\Σ}}}{\left\{{\mathrm{PF}}_p\right[(i+(j-1)\mathrm{\τ}]-{{\mathrm{PF}}_p}^{\mathrm{NN}}[i+(j-1)\mathrm{\τ}\left]\right\}}^2=R_m^2\left(i\right)+{|{\mathrm{PF}}_p(i+\mathrm{m\τ})-{{\mathrm{PF}}_p}^{\mathrm{NN}}(i+\mathrm{m\τ})|}^2---\left(20\right)$

由m增加1而引起的两近邻点间距离的变化是 

${[R_{m+1}^2\left(i\right)-R_m^2\left(i\right)]}^{1/2}=|{\mathrm{PF}}_p(i+\mathrm{m\τ})-{{\mathrm{PF}}_p}^{\mathrm{NN}}(i+\mathrm{m\τ})|---\left(21\right)$

定义新的距离E(m) 

$E\left(m\right)=\frac{1}{N-\mathrm{m\τ}}\underset{i=1}{\overset{N-\mathrm{m\τ}}{\mathrm{\Σ}}}|{\mathrm{PF}}_p(i+\mathrm{m\τ})-{{\mathrm{PF}}_p}^{\mathrm{NN}}(i+\mathrm{m\τ})|---\left(22\right)$

E₁(m)＝E(m+1)/E(m)    (23) 

式中：E₁(m)为两个新距离的比值，当E₁(m)在m＞m₀时，E₁(m)＝1不再变化，则m₀即为 确定的嵌入维数。 

所述步骤(4)中，所述的多核最小二乘支持向量机回归模型中的多核核函数，表示为 

K＝λ₁K_poly+λ₂K_RBF+λ₃K_exp

式中，K为多核函数，K_ploy为多项式核函数，K_poly＝(x_ix_j)+c)^d；K_RBF为径向基函数，K_RBF＝exp(-(x_i-x_j)²/(2σ²))；K_exp为指数核函数，K_exp＝exp(-||x_i-x_j||/(2σ′²))；λ₁，λ₂，λ₃分别为核函数的权值， 且0≤λ_k≤1。 

参照图2，以东北某风电场2010年12月份的实测风速时间数列作为实验样本，对本发明的预测模型进行了验证。图2所示为该月的实测原始风速时间序列，每小时作为一个采样点，共720点。可见，原始风速序列波动较剧烈，且无明显变化规律。将后10天的风速时间序列，采样点481-720，共240点作为测试样本，按照相空间重构算法，依次利用最近480点预测下一个时刻的风速。以第481点为例，预测过程如下： 

(1)将前480个风速样本点进行局域均值分解分解，得到4个PF分量和剩余分量如图3所示。可以看出，除PF₁和PF₂外，其余分量均具有明显变换规律且比较平稳，故一定程度上降低了建模难度。 

(2)对PF分量建立多核最小二乘支持向量机预测模型。根据各自的预测模型对每一个分量进行预测，叠加得到第481点的风速值。 

为了提高预测精度，降低等价核函数K对单个核函数参数的依赖性，权值λ_k(k＝1，2，3)以每个核函数预测的均方根误差σ_RMS作为选择标准，σ_RMS较小的给较大的权值。λ_k选择方法如下： 

${\mathrm{\λ}}_k=\frac{\underset{k=1}{\overset{3}{\mathrm{\Σ}}}{\mathrm{\σ}}_{{\mathrm{RMS}}_k}-{\mathrm{\σ}}_{{\mathrm{RMS}}_k}}{2\underset{k=1}{\overset{3}{\mathrm{\Σ}}}{\mathrm{\σ}}_{{\mathrm{EMS}}_k}}$

通过计算得到核函数的权值分别为：λ₁＝0.43，λ₂＝0.31，λ₃＝0.26。对采样点481-720逐个进行预测，预测结果和测试样本对比见图4。从图4可见，预测值和实测值在大多时刻吻合较好，说明所建立的预测模型符合该风速的变化规律。 

选择合理的误差指标，对评定预测效果的好坏有着重要的作用。本发明运用用下式对 预测效果进行检验： 

${\mathrm{\σ}}_{\mathrm{MAE}}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}|\stackrel{\‾}{y}\left(i\right)-y\left(i\right)|$

${\mathrm{\σ}}_{\mathrm{RMS}}=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{[\stackrel{\‾}{y}\left(i\right)-y\left(i\right)]}^2}$

式中：σ_MAE为平均绝对误差，σ_RMS为均方根误差；N为预测样本数； 为序列的真实值；y(i)为预测结果。 

为了进一步验证本发明的有效性，分别对前480个风速样本点进行经验模式分解Empirical Mode Decomposition，EMD分解，得到的预测结果如表1所示。 

从表中可以看出，当采用均方根误差衡量预测结果时，较经验模式分解法降低了约16.5％，说明本发明的模型具有良好的预测精度。 

表1 

具体实施方式给出本发明的一种风电场风速时间序列预测方法并非穷举，本领域技术人员不经过创造性劳动的简单复制和改进，仍属于本发明权利要求保护的范围。 

Claims (3)

1.一种风电场风速时间序列预测方法，其特征在于，它包括下述步骤：

(1)利用风速采集仪器每隔一小时记录一次同一地区的风速数据，整理采集的原始风速数据，形成风速时间序列用于分析预测；

(2)运用局域均值分解算法对风速时间序列进行多尺度分解，分解得到多个PF分量，局域均值分解算法步骤如下：

(a)找出原始信号x(t)所有的局部极值点n_i，求出所有相邻的局部极值点的平均值：

$m_i=\frac{n_i+n_{i+1}}{2}---\left(1\right)$

将所有相邻的平均值点m_i用直线连接起来，然后用滑动平均法进行平滑处理，得到局部均值函数m₁₁(t)；

(b)求出包络估计值

$a_i=\frac{{|n}_i-n_{i+1}|}{2}---\left(2\right)$

将所有相邻的平均值点a_i用直线连接起来，然后用滑动平均法进行平滑处理，得到包络估计函数a₁₁(t)；

(c)将局部均值函数m₁₁(t)从原始信号x(t)中分离出来，得到

h ₁₁ (t)＝x(t)-m ₁₁ (t)         (3)

(d)用h₁₁(t)除以包络估计函数a₁₁(t)以对h₁₁(t)进行解调，得到

s₁₁(t)＝h₁₁(t)/a₁₁(t)          (4)

对s₁₁(t)重复上述局域均值分解算法的步骤便能得到s₁₁(t)的包络估计函数a₁₂(t)，假如a₁₂(t)不等于1，说明s₁₁(t)不是一个纯调频信号，需要重复步骤(d)中迭代过程n次，直至s_1n(t)为一个纯调频信号，也就是s_1n(t)的包络估计函数a_1(n+1)(t)＝1，所以有

$\left\{\begin{array}{c}{h}_{11}\left(t\right)=x\left(t\right)-m_{11}\left(t\right)\\ h_{12}\left(t\right)=s_{11}\left(t\right)-m_{12}\left(t\right)\\ .\\ .\\ .\\ h_{1n}\left(t\right)=s_{1(n-1)}\left(t\right)-m_{1n}\left(t\right)\end{array}---\left(5\right)\right.$

式中

$\left\{\begin{array}{c}{s}_{11}\left(t\right)=h_{11}\left(t\right)/a_{11}\left(t\right)\\ s_{12}\left(t\right)=h_{12}\left(t\right)/a_{12}\left(t\right)\\ .\\ .\\ .\\ s_{1n}\left(t\right)=h_{1n}\left(t\right)/a_{1n}\left(t\right)\end{array}---\left(6\right)\right.$

迭代终止条件为

$\underset{n\→\∞}{\lim }{a}_{1n}\left(t\right)=1---\left(7\right)$

为了减少迭代次数，降低运算时间，设置一个变量Δ＝10^-4，使得当满足1-Δ≤a_1n(t)≤1+Δ时，终止迭代；

(e)将迭代过程中所产生的全部包络估计函数相乘，得到包络信号瞬时幅值函数：

$a_1\left(t\right)=a_{11}\left(t\right)a_{12}\left(t\right)\·\·\·a_{1n}\left(t\right)=\underset{q=1}{\overset{n}{\mathrm{\Π}}}{a}_{1q}\left(t\right)---\left(8\right)$

(f)将包络信号a₁(t)和纯调频信号s_1n(t)相乘得：

PF₁(t)＝a₁(t)s_1n(t)        (9)

它为信号x(t)的第一个PF分量，包含了原始信号的最高频率成分，是一个单分量的调幅、调频信号，瞬时频率f₁(t)则由调频信号s_1n(t)求出：

$f_1\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{d\left[\arccos \right]\left(s_{1n}\left(t\right)\right)}{\mathrm{dt}}---\left(10\right)$

(g)从原始信号x(t)中将第一个PF分量PF₁(t)分离出来，得到剩余信号u₁(t)，由于剩余信号u₁(t)中还包含有较多的频率成分，因此将u₁(t)作为原始数据重复以上步骤对其进行分解，得到第二个PF分量，重复这个过程直到u_k为一个单调函数为止，得到一定数量的PF分量：

$\left\{\begin{array}{c}{u}_1\left(t\right)=x\left(t\right)-{\mathrm{PF}}_1\left(t\right)\\ u_2\left(t\right)=u_1\left(t\right)-{\mathrm{PF}}_2\left(t\right)\\ .\\ .\\ .\\ u_k\left(t\right)=u_{k-1}\left(t\right)-{\mathrm{PF}}_k\left(t\right)\end{array}\right.---\left(11\right)$

最终信号x(t)可表示为k个PF分量和余量之和：

$\left(t\right)=\underset{p=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{PF}}_p\left(t\right)-u_k\left(t\right)---\left(12\right)$

(3)通过对各个PF分量的延迟时间τ和嵌入维数m的计算，运用相空间重构理论对各个PF分量进行相空间重构；

y_PF(t)＝{PF_p(t),PF_p(t+τ),…,PF_p[t+(m-1)τ]}      (13)

y_u(t)＝{u_k(t),u_k(t+τ),…,u_k[t+(m-1)τ]}       (14)

t＝1,2,….,p＝1,2,…,k

(4)利用多核最小二乘支持向量机回归模型对相空间重构后的各个PF分量进行建模预测；

$y_{\mathrm{PFl}}(t+\mathrm{m\τ})=\underset{i=t}{\overset{t+m-1}{\mathrm{\Σ}}}{\mathrm{\α}}_{i}K(y_{\mathrm{PFi}},x)+b,l=\mathrm{1,2}\·\·\·,k---\left(15\right)$

$y_{\mathrm{uk}}(t+\mathrm{m\τ})=\underset{i=t}{\overset{t+m-1}{\mathrm{\Σ}}}{\mathrm{\α}}_{i}K(y_{\mathrm{ui}},x)+b---\left(16\right)$

(5)将预测结果进行叠加得到最终的风速时间序列的预测结果

$\hat{y}(t+\mathrm{\τ})=\underset{l=1}{\overset{k}{\mathrm{\Σ}}}{y}_{\mathrm{PFl}}(t+\mathrm{\τ})+y_{\mathrm{uk}}(t+\mathrm{\τ})---\left(17\right)$

所述各个PF分量(p＝1,2,…,k)的延迟时间使用自相关法获得，所述各个PF分量的嵌入维数使用假近邻法获得，自相关法计算延迟时间τ的方法为：按照式(18)计算各PF_p分量跨度为jτ的自相关函数R(jτ)

$R\left(\mathrm{j\τ}\right)=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}[{\mathrm{PF}}_p\left(i\right)\·{\mathrm{PF}}_p(i+\mathrm{j\τ})]---\left(18\right)$

式中PF_pi(i＝1,2,…,N)为PF分量，固定j，作出自相关函数R(jτ)关于时间τ(即取τ＝1,2,…)的函数图像，则R(jτ)第一次达到零点时所对应的τ即为延迟时间；

假近邻法确定嵌入维数方法具体为：在m维空间中，每个相点矢量为X(i)＝{PF_p(i),PF_p(i+τ),…,PF_p(i+(m-1)τ)}，都有一个某距离内的最近邻点X^NN(i)，其距离为R_m(i)，则

$R_m^2\left(i\right)={\left|\right|X\left(i\right)-X^{\mathrm{NN}}\left(i\right)\left|\right|}^2=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\left\{{\mathrm{PF}}_p\right[(i+(j-1)\mathrm{\τ}]-{{\mathrm{PF}}_p}^{\mathrm{NN}}[i+(j-1)\mathrm{\τ}\left]\right\}}^2---\left(19\right)$

当相空间的维数增加到m+1维时，这两个相点的距离就会发生变化，两者的距离成为R_m+1(i)，且有

$R_{m+1}^2\left(i\right)=\underset{j=1}{\overset{m+1}{\mathrm{\Σ}}}{\left\{{\mathrm{PF}}_p\right[(i+(j+1)\mathrm{\τ}]-{{\mathrm{PF}}_p}^{\mathrm{NN}}[i+(j-1)\mathrm{\τ}\left]\right\}}^2=R_m^2\left(i\right)+{|{\mathrm{PF}}_p(i+\mathrm{m\τ})-{{\mathrm{PF}}_p}^{\mathrm{NN}}(i+\mathrm{m\τ})|}^2---\left(20\right)$

由m增加1而引起的两近邻点间距离的变化时

${[R_{m+1}^2\left(i\right)-R_m^2\left(i\right)]}^{1/2}=|{\mathrm{Pf}}_p(i-\mathrm{m\τ})-{{\mathrm{PF}}_p}^{\mathrm{NN}}(i+\mathrm{m\τ})|---\left(21\right)$

定义新距离E(m)

$E\left(m\right)=\frac{1}{N-\mathrm{m\τ}}\underset{i=1}{\overset{N-\mathrm{m\τ}}{\mathrm{\Σ}}}|{\mathrm{PF}}_p(i+\mathrm{m\τ})-{{\mathrm{PF}}_p}^{\mathrm{NN}}(i+\mathrm{m\τ})|---\left(22\right)$

E₁(m)＝E(m+1)/E(m)         (23)

式中：E₁(m)为两个新距离的比值，当E₁(m)在m>m₀时，E₁(m)＝1不再变化，则m₀即为确定的嵌入维数。

2.根据权利要求1所述的一种风电场风速时间序列预测方法，其特征在于：所述原始风速数据中相邻两个风速数据点之间的采集间隔为15～60min。

3.根据权利要求1所述的一种风电场风速时间序列预测方法，其特征在于：所述步骤(4)中，所述的多核最小二乘支持向量机回归模型中的多核核函数,表示为

K＝λ₁K_poly+λ₂K_RBF+λ₃K_exp

式中，K为多核函数，K_ploy为多项式核函数，K_poly＝(x_ix_j)+c)^d；K_RBF为径向基函数，K_RBF＝_exp(-(x_i-x_j)²/(2σ²))；K_exp为指数核函数，K_exp＝exp(-||x_i-x_j||/(2σ'²))；λ₁,λ₂,λ₃分别为核函数的权值，且0≤λ_k≤1。