CN103617462A - 一种基于地理统计学的风电场风速时空数据建模方法 - Google Patents

Abstract

本发明公开了一种基于地理统计学的风电场风速时空数据建模方法，包括以下步骤：1：通过风机地理位置信息、风电场风速的空间协方差函数和风电场风速的变异函数构造出风电场的空间结构矩阵，来表征各风机输入风速间的空间相关性；2：应用泛克里格法和贝叶斯算法对各风机输入风速进行分层建模，并采用Gibbs抽样对模型参数进行估计；3：对各台风机风速的向前P步进行预测，得到各台风机风速的P步向前预测分布的模拟样本，然后再通过抽样取平均求得台风机风速的最优预测结果，其中P≥1。本发明从风的物理特性出发，综合分析不同风机处的风速数据在时间和空间上存在的相关性，进而对整个风电场建立更精确的预测模型，所以其预测效果比传统方法好。

Description

一种基于地理统计学的风电场风速时空数据建模方法

技术领域

本发明涉及风力发电技术领域，具体涉及一种基于地理统计学的风电场风速时空数据建模方法。

背景技术

随着风电装机容量在电力系统中所占比例的增加，风能的随机性、间歇性、低能量密度等特点给电力系统的安全稳定和经济运行带来一系列严重的影响。尤其是风电穿透功率超过一定值后，风电的接入可能会严重影响电能质量和电力系统运行，并可能会危及常规发电方式。鉴于此，人们对风电场风速预测与风力发电功率预测进行了大量的研究。

传统的风电场风速预测方法是采用统计方法对平均风速的历史数据进行统计预测。虽然这些方法在某些条件下能取得很好的预测结果，但是对其预测机理和预测结果不能给出充分的解释。此外，这些模型都简单的将风电场当作一台风机来处理，而忽略了同一风电场不同风机处的风速间存在的相关关系和相互影响。

发明内容

针对风电场风速数据具有时空特性的这一特点，本发明提出了一种基于地理统计学的风电场风速时空数据建模方法，从而实现对风电场多风机的风速进行预测。

本发明所采用的技术方案是：一种基于地理统计学的风电场风速时空数据建模方法，其特征在于，包括以下步骤：

步骤1：通过风机地理位置信息、风电场风速的空间协方差函数和风电场风速的变异函数构造出风电场的空间结构矩阵，来表征各风机输入风速间的空间相关性；

步骤2：应用泛克里格法和贝叶斯算法对各风机输入风速进行分层建模，并采用Gibbs抽样对模型参数进行估计；

步骤3：对各台风机风速的向前P步进行预测，得到各台风机风速的P步向前预测分布的模拟样本，然后再通过抽样取平均求得台风机风速的最优预测结果，其中P≥1。

作为优选，步骤1中所述的风电场风速的空间协方差函数，其具体构建过程为：假设W(s,t)是定义在D×T的时空随机过程，其中D代表空间区域，T代表时间，h为两观测点(s_i,s_j)间的空间距离，在时刻t（t=1,2,…,T），定义任意两观测点间的协方差函数为：

C_t(h,t)=Cov(W(s_i,t),W(s_j,t))    （式壹）

采用参数化的协方差函数模型来表示：

C_t(h,t)=σ²κ(s_i,s_j;Φ)    （式贰）

式中σ²为随着时间变化的变异参数，κ(s_i,s_j;Φ)为核函数，其选用Matern族相关函数形式的核函数：

$\begin{array}{c}\mathrm{\κ}(s_i,s_j;\mathrm{\Φ};\mathrm{\υ})=\frac{1}{2^{\mathrm{\υ}-1}\mathrm{\Γ}\left(\mathrm{\υ}\right)}{\left(2\sqrt{\mathrm{\υ}}{d}_{\mathrm{ij}}\mathrm{\φ}\right)}^{\mathrm{\υ}}{K}_{\mathrm{\υ}}\left(2\sqrt{\mathrm{\υ}}{d}_{\mathrm{ij}}\mathrm{\φ}\right)\\ \mathrm{\φ}>0,\mathrm{\υ}>0\end{array}$     （式叁）

式中d_ij为风机站点s_i与s_j之间的测地线距离；Κ_υ(.)是阶次为υ的第二类修正Bessel函数，取υ=1；φ为空间衰减参数，取φ=0.02；在时刻t时，任意两点间的时空协方差是关于这两点距离的函数。

作为优选，步骤1中所述的风电场风速的变异函数，其具体构建过程为：假设W(s,t)是定义在D×T的时空随机过程，其中D代表空间区域，T代表时间，h为两观测点(s_i,s_j)间的空间距离，在时刻t（t=1,2,…,T），对于一个时空平稳过程W(s_i,t),假设N为距离等于h的点对数量，定义距离为h的变异函数γ(h)为：

$\begin{array}{c}\mathrm{\γ}\left(h\right)=\frac{1}{2}\mathrm{var}[W(s_i,t)-W(s_j,t)]\\ =\frac{1}{2(N-1)}\underset{t=1}{\overset{T-1}{\mathrm{\Σ}}}{\{W(s_i,t)-W(s_j,t)\}}^2\end{array}$     （式肆）

s_i,s_j表示所有距离为h的点对，γ(h)是以任意两点之间距离h为变量的函数。

作为优选，步骤1中所述的构造出风电场的空间结构矩阵，其具体实现过程为：空间结构矩阵X_t的行数n由风机数量决定，列数为p，由两部分组成，其中第一部分q列定义为矩阵F，F为n×q维矩阵，表示样本点的空间坐标信息，取q=3，即矩阵F的第1列元素为1，第2、3两列为相对应的风机空间坐标经度和纬度值，其中第二部分p-q列定义为与空间协方差结构相关的空间方向场，假设协方差矩阵∑_v及

均为非奇异矩阵，引入混合能量矩阵B：

$B={\mathrm{\Σ}}_v^{-1}-{\mathrm{\Σ}}_v^{-1}F{\left(F^{\′}{\mathrm{\Σ}}_v^{-1}F\right)}^{-1}{F}^{\′}{\mathrm{\Σ}}_v^{-1}$

对它进行谱分解，可得：

B=UEU^-1

Bu_i=e_iu_i

其中，U=(u₁,...,u_n)，E=diag(e₁,...,e_n)，u₁,…,u_n是矩阵U的特征向量，矩阵E是一个对角矩阵，对角线上的元素为特征值，按照特征值非降的方向将特征值进行重新排序，则有：e₁=...=e_q=0<e_q+1≤...≤e_n，规定空间结构矩阵X_t的后p-q列等于e_i∑_vu_i，其中，i=q+1,...,p，最终，空间结构矩阵X_t可由式伍的形式进行构造：

X_t=(F,e_q+1∑_vu_q+1,...,e_m∑_vu_p)    （式伍）

其中，矩阵X_t的列向量组两两正交。

作为优选，步骤2中所述的应用泛克里格法和贝叶斯算法对各风机输入风速进行分层建模，并采用Gibbs抽样对模型参数进行估计，其具体实现过程为：假设W(s,t)是定义在D×T的时空随机过程，其中D代表空间区域，T代表时间，h为两观测点(s_i,s_j)间的空间距离，在时刻t（t=1,2,…,T）；

首先，对于非平稳的时空随机过程W(s,t)，假设数据中存在主导趋势和残差分量两部分，即：

W(s,t)=μ(s,t)+v(s,t)    （式陆）

其中μ(s,t)为在点s_i处的数学期望，表示W(s,t)的规则而连续的变化，且该趋势可以用一个确定的函数或多项式来拟合，v(s,t)称为残差分量，是围绕主导趋势μ(s,t)摆动的小范围区域内的空间变异值，其数学期望为零；

假设任意时刻t残差向量v(s)服从均值为零，协方差矩阵为∑_v的多元正态分布，即：v(s)～N(0,∑_v)，协方差矩阵∑_v由以下元素组成：

$\begin{array}{c}{\mathrm{\&Sigma;}}_{\mathrm{ij}}=\mathrm{Cov}(Z\left(s_i\right),Z\left(s_j\right))\\ ={\mathrm{\&sigma;}}_v^2\mathrm{\&kappa;}(s_i,s_j;\mathrm{\&phi;})\end{array}$     （式柒）

其中

为变异参数，s_i，s_j为风机站点，i,j=1,2,…,n，Κ_υ(.)是阶次为υ的第二类修正Bessel函数，φ为空间衰减参数，假设

服从超参数为与

的Gamma先验分布，即：

$p\left({\mathrm{\&tau;}}_v\right)~\mathrm{Gamma}(a_{{\mathrm{\&tau;}}_v},b_{{\mathrm{\&tau;}}_v})$     （式捌）；

其次，对于时空随机过程Z_t=(Z(s₁,t),…,Z(s_n,t))'为对应观测点s_i在时刻t的观测值，基于泛克里格的贝叶斯动态模型如式玖和拾所示：

Z_t=X_tθ_t+v_t,v_t～N[0,∑_v]    （式玖）

θ_t=G_tθ_t-1+ω_t,ω_t～N[0,∑_ω]    （式拾）

其中Z_t是n×1维观测向量，X_t是n×p维表征观测点空间相关性的空间结构矩阵，量化了多时空尺度数据模型的空间构成，θ_t是p×1维表征观测点随着时间变化的时间状态变量，反映了各个风机风速的随着时间变化而变化的动态特性，G_t是p×p维状态转移矩阵，v_t是n维空间变异向量，ω_t是p维观测噪声向量，v_t和ω_t相互独立，v_t的协方差矩阵为∑_v，对于ω_t的协方差矩阵∑_ω，设定一个与∑_ω有关的矩阵W，其中

且W服从自由度为2a_ω，协方差矩阵为2b_ω的Wishart先验分布，即：

p(W)～Wishart(2a_ω,2b_ω)    （式十一）；

然后，模型预测；模型的对数似然函数表达式为：

$\begin{array}{c}\log \left\{f(z_1,...,z_T|\mathrm{\&Lambda;})\right\}\&Proportional;-\frac{T}{2}\log \left|{\mathrm{\&Sigma;}}_v\right|\\ -\frac{1}{2}\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{(z_t-{\mathrm{\&theta;}}_t)}^{\&prime;}{\mathrm{\&Sigma;}}_v^{-1}(z_t-{\mathrm{\&theta;}}_t)\end{array}$     （式十二）

其中Λ代表参数θ,τ,W,G的集合；则函数的联合后验密度函数为：

π(Λ|z₁,…,z_T)∝f(z₁,…,z_T|Λ)π(Λ)    （式十三）

其中π(Λ)为Λ中参数的先验分布；

最后，并采用Gibbs抽样对模型参数进行估计；假设第k-1次开始时参数估计值为θ^(k-1)，

W^(k-1)，G^(k-1)，则第k次抽样的迭代步骤如下：

a.由后验分布π(θ|Λ^(k-1);Z_t)抽取θ^(k)；

b.由后验分布π(τ_v|Λ^(k-1);Z_t)抽取

c.由后验分布π(W|Λ^(k-1);Z_t)抽取W^(k)；

d.由后验分布π(G|Λ^(k-1);Z_t)抽取G^(k)；

上述步骤完成一次Gibbs抽样，重复上述过程直至参数的后验条件分布为平稳分布为止；假设经过M次迭代，参数的后验条件分布变为平稳，则得到长度为M的Gibbs抽样序列，舍弃前m个值，则序列变为Λ^m+1,…,Λ^M，作为最终的参数估计数据；

在此基础上对风速预报进行贝叶斯分析，向前1步的后验预测分布为：

π(z_T+1|z₁,…,z_T)=∫π(z_T+1|Λ)π(Λ|z₁,…,z_T)dΛ    （式十四）

上式所示密度函数的期望E(z_T+1|z₁,…,z_T)即为向前1步预测的最优结果，本发明通过对密度函数π(z_T+1|Λ)抽样取平均获得E(z_T+1|z₁,…,z_T)。

作为优选，步骤5中所述的对各台风机风速的向前P步进行预测，其具体实现过程为：首先对密度函数π(z_T+1|Λ)进行抽样，得到

然后，对于j=2,3,…,P，从后验密度函数

抽取

从风速预测密度函数

中抽取由此可得到风速的P步向前预测分布的模拟样本，然后再通过抽样取平均求得最优预测结果。

作为优选，τ_v的先验分布满足式捌所示Gamma分布，其中

W的先验分布满足式十一所示Wishart分布，其中2a_ω=p,2b_ω=0.01I_p,I_p为p×p的单位矩阵；状态变量θ在t=0时刻的初值θ₀服从以下分布p(θ₀)～N[0,C]，其中：C=100I_p；状态转移矩阵G在t=0时刻的初值G₀：G₀=I_p。

本发明考虑风电场风速数据时空特性的风速预测方法就是从风的物理特性出发，综合分析不同风机处的风速数据在时间和空间上存在的相关性，进而对整个风电场建立更精确的预测模型，所以其预测效果比传统方法好。

附图说明

图1：为本发明的算法流程图。

图2：为本发明实施例的风电场风机位置分布图。

具体实施方式

以下结合附图和具体实施例对本发明做进一步的阐述。

请见图1、图2，本发明所采用的技术方案是：一种基于地理统计学的风电场风速时空数据建模方法，其特征在于，包括以下步骤：

步骤1：通过风机地理位置信息、风电场风速的空间协方差函数和风电场风速的变异函数构造出风电场的空间结构矩阵，来表征各风机输入风速间的空间相关性；

其中，风速数据预处理的过程为：由于风速实测数据均是非负值，且由于气候的原因，出现短时间的大风速的概率并不大，风速统计分布实际上是向一侧偏斜的，并非服从正态分布。为了更好的建立概率分层模型，需先对多尺度实测风速数据进行适当的预处理，使之经过适当的变换后能更适合模型的要求。根据风速数据可知该风电场风速主要集中在6-17m/s，将风速小于6m/s的实测记录值全设置为6m/s，大于17m/s的实测记录值设置为17m/s，利用Box-Cox变换对多时空尺度实测风速数据Z进行适当的预处理，处理过程如下：

$Z^{\left(\mathrm{\&lambda;}\right)}=\left\{\begin{array}{cc}(Z^{\mathrm{\&lambda;}}-1)/\mathrm{\&lambda;},& \mathrm{\&lambda;}\&NotEqual;0\\ \ln Z,& \mathrm{\&lambda;}=0\end{array}\right.$

以上变换可以使模型满足线性、独立性、方差齐性、正态性的同时又不丢失实测数据中的信息。通过以上的变换，可以在一定程度上减小多时空尺度实测数据中不可观测到的误差和预测值之间的相关性。但是要适当选择参数λ的值，选取一段观测数据，通过使用极大似然估计得到参数λ的值。对于预测结果应先进行Box-Cox逆变换。

协方差函数与变异函数是用于考察时空数据相关性的两种度量方法。

其中，求取风电场空间协方差函数过程为：假设W(s,t)是定义在D×T的时空随机过程，其中D代表空间区域，T代表时间，h为两观测点(s_i,s_j)间的空间距离，在时刻t（t=1,2,…,T），定义任意两观测点间的协方差函数为：

C_t(h,t)=Cov(W(s_i,t),W(s_j,t))    （式壹）

采用参数化的协方差函数模型来表示：

C_t(h,t)=σ²κ(s_i,s_j;Φ)    （式贰）

式中σ²为随着时间变化的变异参数，κ(s_i,s_j;Φ)为核函数，其选用Matern族相关函数形式的核函数：

$\begin{array}{c}\mathrm{\&kappa;}(s_i,s_j;\mathrm{\&Phi;};\mathrm{\&upsi;})=\frac{1}{2^{\mathrm{\&upsi;}-1}\mathrm{\&Gamma;}\left(\mathrm{\&upsi;}\right)}{\left(2\sqrt{\mathrm{\&upsi;}}{d}_{\mathrm{ij}}\mathrm{\&phi;}\right)}^{\mathrm{\&upsi;}}{K}_{\mathrm{\&upsi;}}\left(2\sqrt{\mathrm{\&upsi;}}{d}_{\mathrm{ij}}\mathrm{\&phi;}\right)\\ \mathrm{\&phi;}>0,\mathrm{\&upsi;}>0\end{array}$     （式叁）

式中d_ij为风机站点s_i与s_j之间的测地线距离；Κ_υ(.)是阶次为υ的第二类修正Bessel函数，取υ=1；φ为空间衰减参数，取φ=0.02；在时刻t时，任意两点间的时空协方差是关于这两点距离的函数。

对于一个时空平稳过程W(s_i,t),假设N为距离等于h的点对数量，定义距离为h的变异函数γ(h)为：

$\begin{array}{c}\mathrm{\&gamma;}\left(h\right)=\frac{1}{2}\mathrm{var}[W(s_i,t)-W(s_j,t)]\\ =\frac{1}{2(N-1)}\underset{t=1}{\overset{T-1}{\mathrm{\&Sigma;}}}{\{W(s_i,t)-W(s_j,t)\}}^2\end{array}$     （式肆）

s_i,s_j表示所有距离为h的点对，γ(h)是以任意两点之间距离h为变量的函数。

不难看出变异函数与协方差函数间有如下关系：

γ(h)=C(0)-C(h)

若以h为横坐标，γ(h)为纵坐标作图，便得到变差函数散点图。运用数学公式拟合各点,就可得到变差函数曲线。变异函数γ(h)是一个单调递增函数，当h超过某一数值（变程a）后，γ(h)不再继续单调地增大，而稳定在一个极限值附近，该极限值称为基台值，即C(0)，基态值的大小反映变量变化幅度的大小。变程a反映了变量变化的空间范围。

步骤2：应用泛克里格法和贝叶斯算法对各风机输入风速进行分层建模，并采用Gibbs抽样对模型参数进行估计；其具体实现过程为：空间结构矩阵X_t的行数n由风机数量决定，列数为p，由两部分组成，其中第一部分q列定义为矩阵F，F为n×q维矩阵，表示样本点的空间坐标信息，取q=3，即矩阵F的第1列元素为1，第2、3两列为相对应的风机空间坐标经度和纬度值，其中第二部分p-q列定义为与空间协方差结构相关的空间方向场，假设协方差矩阵∑_v及均为非奇异矩阵，引入混合能量矩阵B：

$B={\mathrm{\&Sigma;}}_v^{-1}-{\mathrm{\&Sigma;}}_v^{-1}F{\left(F^{\&prime;}{\mathrm{\&Sigma;}}_v^{-1}F\right)}^{-1}{F}^{\&prime;}{\mathrm{\&Sigma;}}_v^{-1}$

对它进行谱分解，可得：

B=UEU^-1

Bu_i=e_iu_i

其中，U=(u₁,...,u_n)，E=diag(e₁,...,e_n)，u₁,…,u_n是矩阵U的特征向量，矩阵E是一个对角矩阵，对角线上的元素为特征值，按照特征值非降的方向将特征值进行重新排序，则有：e₁=...=e_q=0<e_q+1≤...≤e_n，规定空间结构矩阵X_t的后p-q列等于e_i∑_vu_i，其中，i=q+1,...,p，最终，空间结构矩阵X_t可由式伍的形式进行构造：

X_t=(F,e_q+1∑_vu_q+1,...,e_m∑_vu_p)    （式伍）

其中，矩阵X_t的列向量组两两正交。

5.根据权利要求1所述的基于地理统计学的风电场风速时空数据建模方法，其特征在于：

步骤2中所述的应用泛克里格法和贝叶斯算法对各风机输入风速进行分层建模，并采用Gibbs抽样对模型参数进行估计，其具体实现过程为：假设W(s,t)是定义在D×T的时空随机过程，其中D代表空间区域，T代表时间，h为两观测点(s_i,s_j)间的空间距离，在时刻t（t=1,2,…,T）；

首先，对于非平稳的时空随机过程W(s,t)，假设数据中存在主导趋势和残差分量两部分，即：

W(s,t)=μ(s,t)+v(s,t)    （式陆）

其中μ(s,t)为在点s_i处的数学期望，表示W(s,t)的规则而连续的变化，且该趋势可以用一个确定的函数或多项式来拟合，v(s,t)称为残差分量，是围绕主导趋势μ(s,t)摆动的小范围区域内的空间变异值，其数学期望为零；

假设任意时刻t残差向量v(s)服从均值为零，协方差矩阵为∑_v的多元正态分布，即：v(s)～N(0,∑_v)，协方差矩阵∑_v由以下元素组成：

$\begin{array}{c}{\mathrm{\&Sigma;}}_{\mathrm{ij}}=\mathrm{Cov}(Z\left(s_i\right),Z\left(s_j\right))\\ ={\mathrm{\&sigma;}}_v^2\mathrm{\&kappa;}(s_i,s_j;\mathrm{\&phi;})\end{array}$     （式柒）

其中

为变异参数，s_i，s_j为风机站点，i,j=1,2,…,n，Κ_υ(.)是阶次为υ的第二类修正Bessel函数，φ为空间衰减参数，假设

服从超参数为

与

的Gamma先验分布，即：

$p\left({\mathrm{\&tau;}}_v\right)~\mathrm{Gamma}(a_{{\mathrm{\&tau;}}_v},b_{{\mathrm{\&tau;}}_v})$     （式捌）；

其次，对于时空随机过程Z_t=(Z(s₁,t),…,Z(s_n,t))'为对应观测点s_i在时刻t的观测值，基于泛克里格的贝叶斯动态模型如式玖和拾所示：

Z_t=X_tθ_t+v_t,v_t～N[0,∑_v]    （式玖）

θ_t=G_tθ_t-1+ω_t,ω_t～N[0,∑_ω]    （式拾）

其中Z_t是n×1维观测向量，X_t是n×p维表征观测点空间相关性的空间结构矩阵，量化了多时空尺度数据模型的空间构成，θ_t是p×1维表征观测点随着时间变化的时间状态变量，反映了各个风机风速的随着时间变化而变化的动态特性，G_t是p×p维状态转移矩阵，v_t是n维空间变异向量，ω_t是p维观测噪声向量，v_t和ω_t相互独立，v_t的协方差矩阵为∑_v，对于ω_t的协方差矩阵∑_ω，设定一个与∑_ω有关的矩阵W，其中

且W服从自由度为2a_ω，协方差矩阵为2b_ω的Wishart先验分布，即：

p(W)～Wishart(2a_ω,2b_ω)    （式十一）；

然后，模型预测；模型的对数似然函数表达式为：

$\begin{array}{c}\log \left\{f(z_1,...,z_T|\mathrm{\&Lambda;})\right\}\&Proportional;-\frac{T}{2}\log \left|{\mathrm{\&Sigma;}}_v\right|\\ -\frac{1}{2}\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{(z_t-{\mathrm{\&theta;}}_t)}^{\&prime;}{\mathrm{\&Sigma;}}_v^{-1}(z_t-{\mathrm{\&theta;}}_t)\end{array}$     （式十二）

其中Λ代表参数θ,τ,W,G的集合；则函数的联合后验密度函数为：

π(Λ|z₁,…,z_T)∝f(z₁,…,z_T|Λ)π(Λ)    （式十三）

其中π(Λ)为Λ中参数的先验分布；

最后，并采用Gibbs抽样对模型参数进行估计；假设第k-1次开始时参数估计值为θ^(k-1)，W^(k-1)，G^(k-1)，则第k次抽样的迭代步骤如下：

a.由后验分布π(θ|Λ^(k-1);Z_t)抽取θ^(k)；

b.由后验分布π(τ_v|Λ^(k-1);Z_t)抽取

c.由后验分布π(W|Λ^(k-1);Z_t)抽取W^(k)；

d.由后验分布π(G|Λ^(k-1);Z_t)抽取G^(k)；

上述步骤完成一次Gibbs抽样，重复上述过程直至参数的后验条件分布为平稳分布为止；假设经过M次迭代，参数的后验条件分布变为平稳，则得到长度为M的Gibbs抽样序列，舍弃前m个值，则序列变为Λ^m+1,…,Λ^M，作为最终的参数估计数据；

在此基础上对风速预报进行贝叶斯分析，向前1步的后验预测分布为：

π(z_T+1|z₁,…,z_T)=∫π(z_T+1|Λ)π(Λ|z₁,…,z_T)dΛ    （式十四）

上式所示密度函数的期望E(z_T+1|z1,…,z_T)即为向前1步预测的最优结果，本发明通过对密度函数π(z_T+1|Λ)抽样取平均获得E(z_T+1|z₁,…,z_T)。

步骤3：对各台风机风速的向前P步进行预测，得到各台风机风速的P步向前预测分布的模拟样本，然后再通过抽样取平均求得台风机风速的最优预测结果，其中P≥1；其具体实现过程为：首先对密度函数π(z_T+1|Λ)进行抽样，得到

然后，对于j=2,3,…,P，从后验密度函数抽取从风速预测密度函数

中抽取

由此可得到风速的P步向前预测分布的模拟样本，然后再通过抽样取平均求得最优预测结果。

本实施例中，τ_v的先验分布满足式捌所示Gamma分布，其中

W的先验分布满足式十一所示Wishart分布，其中2a_ω=p,2b_ω=0.01I_p,I_p为p×p的单位矩阵；状态变量θ在t=0时刻的初值θ₀服从以下分布p(θ₀)～N[0,C]，其中：C=100I_p；状态转移矩阵G在t=0时刻的初值G₀：G₀=I_p。

本实施例采用常用的均方根误差（Root Mean Squared Error，RMSE）和平均绝对误差（Mean Absolute Error,MAE）作为性能评价指标，其定义为：

$\mathrm{RMSE}=\sqrt{\frac{1}{N^{\#}}\underset{i=1}{\overset{N^{\#}}{\mathrm{\&Sigma;}}}{(Z_{\left(i\right)}-{\hat{Z}}_{\left(i\right)})}^2}$

$\mathrm{MAE}=\frac{1}{N^{\#}}\underset{i=1}{\overset{N^{\#}}{\mathrm{\&Sigma;}}}|Z_{\left(i\right)}-{\hat{Z}}_{\left(i\right)}|$

其中，N^#为风电场群风速预测值中有效的数据个数，Z_(i)为风电场群风速实际值，

为风电场群风速预测值。

以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围，因此，凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于地理统计学的风电场风速时空数据建模方法，其特征在于，包括以下步骤：

步骤1：通过风机地理位置信息、风电场风速的空间协方差函数和风电场风速的变异函数构造出风电场的空间结构矩阵，来表征各风机输入风速间的空间相关性；

步骤2：应用泛克里格法和贝叶斯算法对各风机输入风速进行分层建模，并采用Gibbs抽样对模型参数进行估计；

步骤3：对各台风机风速的向前P步进行预测，得到各台风机风速的P步向前预测分布的模拟样本，然后再通过抽样取平均求得台风机风速的最优预测结果，其中P≥1。

2.根据权利要求1所述的基于地理统计学的风电场风速时空数据建模方法，其特征在于：

步骤1中所述的风电场风速的空间协方差函数，其具体构建过程为：假设W(s,t)是定义在D×T的时空随机过程，其中D代表空间区域，T代表时间，h为两观测点(s_i,s_j)间的空间距离，在时刻t（t=1,2,…,T），定义任意两观测点间的协方差函数为：

C_t(h,t)=Cov(W(s_i,t),W(s_j,t))    （式壹）

采用参数化的协方差函数模型来表示：

C_t(h,t)=σ²κ(s_i,s_j;Φ)    （式贰）

式中σ²为随着时间变化的变异参数，κ(s_i,s_j;Φ)为核函数，其选用Matern族相关函数形式的核函数：

$\begin{array}{c}\mathrm{\&kappa;}(s_i,s_j;\mathrm{\&Phi;};\mathrm{\&upsi;})=\frac{1}{2^{\mathrm{\&upsi;}-1}\mathrm{\&Gamma;}\left(\mathrm{\&upsi;}\right)}{\left(2\sqrt{\mathrm{\&upsi;}}{d}_{\mathrm{ij}}\mathrm{\&phi;}\right)}^{\mathrm{\&upsi;}}{K}_{\mathrm{\&upsi;}}\left(2\sqrt{\mathrm{\&upsi;}}{d}_{\mathrm{ij}}\mathrm{\&phi;}\right)\\ \mathrm{\&phi;}>0,\mathrm{\&upsi;}>0\end{array}$       （式叁）

式中d_ij为风机站点s_i与s_j之间的测地线距离；Κ_υ(.)是阶次为υ的第二类修正Bessel函数，取υ=1；φ为空间衰减参数，取φ=0.02；在时刻t时，任意两点间的时空协方差是关于这两点距离的函数。

3.根据权利要求1所述的基于地理统计学的风电场风速时空数据建模方法，其特征在于：

步骤1中所述的风电场风速的变异函数，其具体构建过程为：假设W(s,t)是定义在D×T的时空随机过程，其中D代表空间区域，T代表时间，h为两观测点(s_i,s_j)间的空间距离，在时刻t（t=1,2,...,T），对于一个时空平稳过程W(s_i,t),假设N为距离等于h的点对数量，定义距离为h的变异函数γ(h)为：

$\begin{array}{c}\mathrm{\&gamma;}\left(h\right)=\frac{1}{2}\mathrm{var}[W(s_i,t)-W(s_j,t)]\\ =\frac{1}{2(N-1)}\underset{t=1}{\overset{T-1}{\mathrm{\&Sigma;}}}{\{W(s_i,t)-W(s_j,t)\}}^2\end{array}$                （式肆）

s_i,s_j表示所有距离为h的点对，γ(h)是以任意两点之间距离h为变量的函数。

4.根据权利要求1所述的基于地理统计学的风电场风速时空数据建模方法，其特征在于：

步骤1中所述的构造出风电场的空间结构矩阵，其具体实现过程为：空间结构矩阵X_t的行数n由风机数量决定，列数为p，由两部分组成，其中第一部分q列定义为矩阵F，F为n×q维矩阵，表示样本点的空间坐标信息，取q=3，即矩阵F的第1列元素为1，第2、3两列为相对应的风机空间坐标经度和纬度值，其中第二部分p-q列定义为与空间协方差结构相关的空间方向场，假设协方差矩阵∑_v及均为非奇异矩阵，引入混合能量矩阵B：

$B={\mathrm{\&Sigma;}}_v^{-1}-{\mathrm{\&Sigma;}}_v^{-1}F{\left(F^{\&prime;}{\mathrm{\&Sigma;}}_v^{-1}F\right)}^{-1}{F}^{\&prime;}{\mathrm{\&Sigma;}}_v^{-1}$

对它进行谱分解，可得：

B=UEU^-1

Bu_i=e_iu_i

其中，U=(u₁,...,u_n)，E=diag(e₁,...,e_n)，u₁,…,u_n是矩阵U的特征向量，矩阵E是一个对角矩阵，对角线上的元素为特征值，按照特征值非降的方向将特征值进行重新排序，则有：e₁=...=e_q=0<e_q+1≤...≤e_n，规定空间结构矩阵X_t的后p-q列等于e_i∑_vu_i，其中，i=q+1,...,p，最终，空间结构矩阵X_t可由式伍的形式进行构造：

X_t=(F,e_q+1∑_vu_q+1,...,e_m∑_vu_p)              （式伍）

其中，矩阵X_t的列向量组两两正交。

5.根据权利要求1所述的基于地理统计学的风电场风速时空数据建模方法，其特征在于：

步骤2中所述的应用泛克里格法和贝叶斯算法对各风机输入风速进行分层建模，并采用Gibbs抽样对模型参数进行估计，其具体实现过程为：假设W(s,t)是定义在D×T的时空随机过程，其中D代表空间区域，T代表时间，h为两观测点(s_i,s_j)间的空间距离，在时刻t（t=1,2,...,T）；

首先，对于非平稳的时空随机过程W(s,t)，假设数据中存在主导趋势和残差分量两部分，即：

W(s,t)=μ(s,t)+v(s,t)                    （式陆）

其中μ(s,t)为在点s_i处的数学期望，表示W(s,t)的规则而连续的变化，且该趋势可以用一个确定的函数或多项式来拟合，v(s,t)称为残差分量，是围绕主导趋势μ(s,t)摆动的小范围区域内的空间变异值，其数学期望为零；

假设任意时刻t残差向量v(s)服从均值为零，协方差矩阵为∑_v的多元正态分布，即：v(s)～N(0,∑_v)，协方差矩阵∑_v由以下元素组成：

$\begin{array}{c}{\mathrm{\&Sigma;}}_{\mathrm{ij}}=\mathrm{Cov}(Z\left(s_i\right),Z\left(s_j\right))\\ ={\mathrm{\&sigma;}}_v^2\mathrm{\&kappa;}(s_i,s_j;\mathrm{\&phi;})\end{array}$                （式柒）

其中为变异参数，s_i，s_j为风机站点，i,j=1,2,...,n，Κ_υ(.)是阶次为υ的第二类修正Bessel函数，φ为空间衰减参数，假设服从超参数为

与

的Gamma先验分布，即：

$p\left({\mathrm{\&tau;}}_v\right)~\mathrm{Gamma}(a_{{\mathrm{\&tau;}}_v},b_{{\mathrm{\&tau;}}_v})$                （式捌）；

其次，对于时空随机过程Z_t=(Z(s₁,t),...,Z(s_n,t))'为对应观测点s_i在时刻t的观测值，基于泛克里格的贝叶斯动态模型如式玖和拾所示：

Z_t=X_tθ_t+v_t,v_t～N[0,∑_v]              （式玖）

θ_t=G_tθ_t-1+ω_t,ω_t～N[0,∑_ω]                （式拾）

其中Z_t是n×1维观测向量，X_t是n×p维表征观测点空间相关性的空间结构矩阵，量化了多时空尺度数据模型的空间构成，θ_t是p×1维表征观测点随着时间变化的时间状态变量，反映了各个风机风速的随着时间变化而变化的动态特性，G_t是p×p维状态转移矩阵，v_t是n维空间变异向量，ω_t是p维观测噪声向量，v_t和ω_t相互独立，v_t的协方差矩阵为∑_v，对于ω_t的协方差矩阵∑_ω，设定一个与∑_ω有关的矩阵W，其中

且W服从自由度为2a_ω，协方差矩阵为2b_ω的Wishart先验分布，即：

p(W)～Wishart(2a_ω,2b_ω)              （式十一）；

然后，模型预测；模型的对数似然函数表达式为：

$\begin{array}{c}\log \left\{f(z_1,...,z_T|\mathrm{\&Lambda;})\right\}\&Proportional;-\frac{T}{2}\log \left|{\mathrm{\&Sigma;}}_v\right|\\ -\frac{1}{2}\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{(z_t-{\mathrm{\&theta;}}_t)}^{\&prime;}{\mathrm{\&Sigma;}}_v^{-1}(z_t-{\mathrm{\&theta;}}_t)\end{array}$               （式十二）

其中Λ代表参数θ,τ,W,G的集合；则函数的联合后验密度函数为：

π(Λ|z₁,...,z_T)∝f(z₁,...,z_T|Λ)π(Λ)              （式十三）

其中π(Λ)为Λ中参数的先验分布；

最后，并采用Gibbs抽样对模型参数进行估计；假设第k-1次开始时参数估计值为θ^(k-1)，

W^(k-1)，G^(k-1)，则第k次抽样的迭代步骤如下：

a.由后验分布π(θ|Λ^(k-1);Z_t)抽取θ^(k)；

b.由后验分布π(τ_v|Λ^(k-1);Z_t)抽取

c.由后验分布π(W|Λ^(k-1);Z_t)抽取W^(k)；

d.由后验分布π(G|Λ^(k-1);Z_t)抽取G^(k)；

上述步骤完成一次Gibbs抽样，重复上述过程直至参数的后验条件分布为平稳分布为止；假设经过M次迭代，参数的后验条件分布变为平稳，则得到长度为M的Gibbs抽样序列，舍弃前m个值，则序列变为Λ^m+1,…,Λ^M，作为最终的参数估计数据；

在此基础上对风速预报进行贝叶斯分析，向前1步的后验预测分布为：

π(z_T+1|z₁,…,z_T)=∫π(z_T+1|Λ)π(Λ|z₁,…,z_T)dΛ    （式十四）

上式所示密度函数的期望E(z_T+1|z₁,…,z_T)即为向前1步预测的最优结果，本发明通过对密度函数π(z_T+1|Λ)抽样取平均获得E(z_T+1|z₁,…,z_T)。

6.根据权利要求5所述的基于地理统计学的风电场风速时空数据建模方法，其特征在于：

步骤5中所述的对各台风机风速的向前P步进行预测，其具体实现过程为：首先对密度函数π(z_T+1|Λ)进行抽样，得到然后，对于j=2,3,…,P，从后验密度函数抽取

从风速预测密度函数

中抽取

由此可得到风速的P步向前预测分布的模拟样本，然后再通过抽样取平均求得最优预测结果。

7.根据权利要求5所述的基于地理统计学的风电场风速时空数据建模方法，其特征在于：

τ_v的先验分布满足式捌所示Gamma分布，其中

W的先验分布满足式十一所示Wishart分布，其中2a_ω=p,2b_ω=0.01I_p,I_p为p×p的单位矩阵；状态变量θ在t=0时刻的初值θ₀服从以下分布p(θ₀)～N[0,C]，其中：C=100I_p；状态转移矩阵G在t=0时刻的初值G₀：G₀=I_p。

Images