CN114462684B - 耦合数值天气预报与实测数据的风速多点同步预测方法 - Google Patents

Abstract

本发明公开了一种耦合数值天气预报与实测数据的风速多点同步预测方法：首先，将不同测点的历史实测阵风风速数据及数值天气预报数据作为模型输入，测点处的当前时刻实测数据作为输出，采用特征选择方法比较各输入变量的重要性程度，并选取程度较高的变量作为最终输入；然后，采用移动窗口法建立基于贝叶斯动态回归的阵风风速预测模型，并对预测结果进行评估，获取高精度风速预测模型；最后，利用该预测模型对多测点处的阵风风速进行长短期预测。本发明可保证阵风风速预测结果准确可靠的基础上，提升计算效率并实现多测点阵风风速的同步预测。

Description

耦合数值天气预报与实测数据的风速多点同步预测方法

技术领域

本发明涉及风速预测领域，特别是涉及一种耦合数值天气预报与实测数据的风速多点同步预测方法。

背景技术

风速预测已广泛应用于多种领域中，比如：风力发电调度、列车车速调节以及保障户外运动安全等方面。准确预测风速对于提高风能利用效率、降低运营成本以及保证人员安全等均具有重要意义。风速监测与数值模拟技术的快速发展为风速预测提供了重要的研究基础和手段。然而，风速的高随机性、强波动性等特点极大地影响了其预测精度和计算效率。因此，高精度风速实时预测方法研究日益受到关注。

现有的主要风速预测方法之一为基于数值天气预报的风速预测方法，该方法主要根据所测地区的风速监测数据，利用计算机进行数值计算，得到风速预测结果。基于数值天气预报的预测结果整体与实际情况相符，但该方法对于某一具体地点的预测精度仍有待提高，且其对数据和硬件要求苛刻，计算量巨大，难以推广使用。另一种为基于数据驱动的统计预测方法，该方法需采集大量风速数据，经过处理分析后对风速进行预测。相较于数值天气预报方法，统计方法更易实现，但由于风速非平稳性特性，导致该方法的预测精度较低。近年来，随着计算机技术性能的不断提升，机器学习算法在解决非线性问题中表现出独特优势，因此被广泛应用于风速预测领域。然而，当处理较多的预测回归变量时，机器学习方法的计算效率较低。贝叶斯线性回归具有计算简单高效的优点常被用于解决预测问题。相较于多元线性回归，该方法可定量估计回归过程中的不确定性。事实上，贝叶斯线性回归属于静态回归方法，即回归系数为常数且不随时间发生变化。然而，静态线性模型仅适合描述回归变量与输出之间关系的局部特征，而非全局变化。因此，亟需兼具高精度和计算效率的风速长短期预测方法。

发明内容

发明目的：针对目前风速预测方法预测精度低、计算费时、难以推广使用等问题，本发明提供了一种耦合数值天气预报与实测数据的阵风风速多点同步预测方法，以有效提升阵风风速不同步长的预测精度，并实现多测点的同步预测。

技术方案：为了实现上述目标，一种耦合数值天气预报与实测数据的风速多点同步预测方法采用如下的技术方案：

第一步：获取目标区域内多测点的实测阵风风速数据及数值天气预报数据；

第二步：以目标区域内t时刻数值天气预报数据及t时刻前n个时刻测点k的实测阵风风速数据构成的向量作为输入，将t时刻多测点的实测阵风风速数据作为输出数据，构建样本集；其中/>表示t时刻数值天气预报数据的第p个变量，x_k,t-i表示第t-i时刻测点k的实测阵风风速数据，i＝1,2,…,n；

第三步：采用特征选择方法比较第二步中样本集各输入变量的重要性程度，选取程度较高的若干变量作为最终输入，得到训练样本集；

第四步：采用移动窗口法，选取不同的预测步长，基于第三步中的训练样本集构建短期风速预测训练样本集和长期风速预测训练样本集；

第五步：基于第四步中的短期和长期风速预测训练样本集，分别对贝叶斯线性动态回归模型和贝叶斯二次动态回归模型进行训练，构建短期风速预测模型和长期风速预测模型；

第六步：基于第五步中的两种风速预测模型，实现多测点阵风风速的长短期预测。

可选的，第五步中贝叶斯线性动态回归模型表示为：

Y_t＝Φ(X_t)w_t+v_t,v_t～N(0,V)

其中，

Y_t＝[Y_1,t Y_2,t … Y_k,t]^T

X_t＝blockdiag(X_1,t X_2,t … … … X_k,t)

w_t＝[w_1,t w_2,t … … … w_k,t]^T

w_k,t＝[α_k,t β_k,t]^T＝[α_k,t β_k1,t β_k2,t … β_kd,t]^T

式中，Y_t表示目标区域内t时刻各测点的阵风风速预测结果构成的矩阵，Y_k,t表示目标区域内t时刻测点k的阵风风速预测值，blockdiag(·)表示块对角矩阵，X_t表示t时刻的输入数据，w_t表示t时刻的动态回归系数矩阵，w_k,t表示t时刻对应测点k的动态回归系数，α_k,t和β_k,t分别表示t时刻对应测点k的截距和回归系数，β_kd,t表示t时刻测点k的第d个回归系数，v_t表示服从0均值、方差V正态分布的回归误差；

动态回归系数的变化为：

w_t＝w_t-1+q_t,q_t～N(0,Q)

其中，w_t-1表示t-1时刻的动态回归系数矩阵，q_t表示服从0均值、方差Q正态分布的状态矩阵误差；

Q＝blockdiag(Q_1,t Q_2,t … … …Q_k,t)

Q_k,t＝blockdiag(Q_kα Q_kβ)

式中，Q_k,t表示t时刻测点k的状态误差，Q_kα和Q_kβ为t时刻对应测点k的截距和回归系数的状态误差。

可选的，第五步中，贝叶斯二次动态回归模型为：

Y_t＝[Y_1,t Y_2,t … Y_k,t]^T

α_t＝[α_1,t α_2,t … α_K,t]^T

式中，Y_t表示目标区域内t时刻各测点的阵风风速预测结果构成的矩阵，Y_k,t表示目标区域内t时刻测点k的阵风风速预测值，α_k,t和β_k,t分别表示t时刻对应测点k的截距和回归系数，X_l,t表示t时刻测点l的输入数据，β_kk,t表示t时刻测点k的平方系数，β_kl,t表示t时刻测点k和l的交叉系数，v_t表示服从0均值、方差V正态分布的回归误差，K表示目标区域内的测点数。

可选的，第五步中，采用卡尔曼滤波与期望最大化EM算法估计贝叶斯线性动态回归模型和贝叶斯二次动态回归模型参数。

可选的，第一步中，数值天气预报数据为欧洲中期天气预报中心数据ECMWF、或第五代大气研究中心中尺度模式MM5、或天气研究与预报模型WRF。

可选的，第三步中，选取数值天气预报数据中10m高度阵风风速、10m高度平均风速、温度以及目标测点对应大气压处的温度和风速。

可选的，第一步中的数值天气预报数据为目标测点所在位置处经纬度0.25°×0.25°范围内九个网格节点的数据。

有益效果：与现有技术相比，本发明具有以下优点：

(1)与数值天气预报预测平均风速不同，贝叶斯动态回归模型采用概率预测方法，通过结合数值天气预报和实测数据，有效提高了长短期时间尺度上的风速预测精度；

(2)以多个网格节点的天气预报数据和现场实测数据作为贝叶斯动态回归模型的输入，考虑了输入变量的空间相关性和阵风的局部变化特性，可实现多个测点的同步预测；

(3)阵风风速预测模型采用移动窗口法，提升了计算效率并实现了阵风风速的多步预测，可广泛用于风速预测领域。

附图说明

图1是耦合数值天气预报与实测数据的风速多点同步预测方法技术流程图；

图2是数值天气预测网格节点和测点位置示意图；

图3是基于移动窗口法的风速多步预测示意图。

具体实施方式

以下结合附图对本发明的技术方案作进一步的详细说明。

本发明实现方案的主要流程具体如下(参见图1)：

第一步：获取目标区域测点的实测阵风风速数据及数值天气预报数据。其中，数值天气预报数据为欧洲中期天气预报中心数据(ECMWF)、第五代大气研究中心中尺度模式(MM5)或天气研究与预报模型(WRF)，且数值天气预报数据为目标测点所在位置的经纬度0.25°×0.25°范围内九个网格点的数据(参见图2)。

第二步：以目标区域内t时刻数值天气预报数据及t时刻前n个时刻第k个测点的实测阵风风速数据作为输入，将t时刻多测点的实测阵风风速数据作为输出数据，构建样本集；其中/>表示t时刻数值天气预报数据的第p个变量，x_k,t-i表示第t-i时刻测点k的实测阵风风速数据，i＝1,2,…,n。

第三步：采用特征选择方法比较第二步中样本集各输入变量的重要性程度，选取程度较高的若干变量作为最终输入，得到训练样本集，其中，最终输入为分别为数值天气预报数据包括10m处阵风风速、10m处平均风速、温度、与目标测点相近大气压处的温度及风速。

第四步：采用移动窗口法(参见图3)，选取不同的预测步长，基于第三步中的训练样本集构建短期风速预测训练样本集和长期风速预测训练样本集。

第五步：基于第四步中的短期和长期风速预测训练样本集，分别对贝叶斯线性动态回归模型和贝叶斯二次动态回归模型进行训练，构建短期风速预测模型和长期风速预测模型。采用期望最大化(EM)算法估计贝叶斯线性动态回归模型和贝叶斯二次动态回归模型中的参数。

贝叶斯线性动态回归模型表示为：

Y_t＝Φ(X_t)w_t+v_t,v_t～N(0,V)

其中，

Y_t＝[Y_1,t Y_2,t … Y_k,t]^T

X_t＝blockdiag(X_1,t X_2,t … … … X_k,t)

w_t＝[w_1,t w_2,t … … … w_k,t]^T

w_k,t＝[α_k,t β_k,t]^T＝[α_k,t β_k1,t β_k2,t … β_kd,t]^T

式中，Y_t表示目标区域内t时刻各测点阵风风速预测结果，Y_k,t表示目标区域内t时刻测点k的阵风风速预测值，Φ(·)为未知的基本函数，blockdiag(·)表示块对角矩阵，X_t表示第t时刻的实测阵风风速数据，X_k,t表示t时刻第k个测点的实测阵风风速数据，w_t为动态回归系数矩阵，w_k,t为t时刻第k个动态回归系数，v_t为服从0均值、方差V正态分布的回归误差，α_k,t和β_k,t分别表示截距和回归系数，β_kd,t为t时刻第k个测点的第d个回归系数。

动态回归系数的变化为：

w_t＝w_t-1+q_t,q_t～N(0,Q)

其中，

Q＝blockdiag(Q_1,t Q_2,t … … … Q_k,t)

Q_k,t＝blockdiag(Q_kα Q_kβ)

式中，q_t为服从0均值、方差Q正态分布的状态矩阵误差，Q_k,t表示t时刻第k个测点的状态误差，Q_kα和Q_kβ为截距和回归系数的状态误差。

矩阵形式的贝叶斯动态模型可表示为：

Y_t＝F_tθ_t+v_t,v_t～N(0,V)

θ_t＝Gθ_t-1+q_t,q_t～N(0,Q)

其中，

F_t＝Φ(X_t)

θ_t＝w_t

G＝I

θ_t～N(f_t,P_t)

式中，F_t表示观测矩阵：G为状态矩阵，I为单位矩阵，θ_t表示均值为f_t、方差为P_t的隐藏状态变量。

可选的，贝叶斯二次动态回归模型为：

Y_t＝[Y_1,t Y_{2, t}… Y_k,t]^T

α_t＝[α_1,t α_2,t … α_K,t]^T

式中，Y_t表示目标区域内t时刻各测点的阵风风速预测结果构成的矩阵，Y_k,t表示目标区域内t时刻测点k的阵风风速预测值，α_k,t和β_k,t分别表示t时刻对应测点k的截距和回归系数，X_l,t表示t时刻测点l的输入数据，β_kk,t表示t时刻测点k的平方系数，β_kl,t表示t时刻测点k和l的交叉系数，v_t表示服从0均值、方差V正态分布的回归误差，K表示目标区域内的测点数。

可选的，采用卡尔曼滤波与期望最大化算法估计贝叶斯动态回归模型的参数Θ＝{V,Q}，通过迭代计算直到期望最大化算法的相邻时间步对数似然函数趋于收敛：

状态变量在t+1时刻的先验分布为p(θ_t+1|Y_1:t)～N(f_t+1|t,P_t+1|t)。卡尔曼滤波可以表示为：

f_t+1|t＝f_t|t

P_t+1|t＝P_t|t+Q

f_t+1＝f_t+1|t+K_t+1(Y_t+1-F_t+1f_t+1|t)

P_t+1＝(I-K_t+1F_t+1)P_t+1|t

其中，(·)_t+1|t表示已知t时刻时t+1时刻的值，K为增益矩阵。

卡尔曼滤波为向前递推算法，其向后递推为初始值f_T|T,P_T|T的卡尔曼平滑，表达式如下

f_t|T＝f_t|t+L_t(f_t+1|T-f_t+1|t)

T为训练集的大小，θ_t的后验分布为N(f_t|T,P_t|T)。

预测结果表示为：

其中，与/>分别表示阵风风速的预测值与方差。

期望最大化算法的E步计算观测值的对数似然函数为：

M步令对数似然函数关于参数的偏导数为零，从而得到新的参数如下：

根据上式，θ_t、及/>的期望函数为

E(θ_t|Y)＝f_t|T

然后，贝叶斯动态回归的参数通过以上步骤，E步与M步不断迭代，直到相邻时间步的对数似然函数趋于收敛。

根据最优的贝叶斯动态回归模型的参数Θ＝{V,Q}，建立基于贝叶斯动态回归算法的阵风风速预测模型。

第六步：基于第五步中的两种风速预测模型，实现多测点阵风风速的长短期预测。

上述说明仅是本发明的优选实施方式，应当指出对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和等同替换，这些对本发明权利要求进行改进和等同替换后的技术方案，均落入本发明的保护范围。

Claims (5)

1.耦合数值天气预报与实测数据的风速多点同步预测方法，其特征在于，该方法包括以下步骤：

第一步：获取目标区域内多测点的实测阵风风速数据及数值天气预报数据；

第二步：以目标区域内t时刻数值天气预报数据及t时刻前n个时刻测点k的实测阵风风速数据构成的向量作为输入，将t时刻多测点的实测阵风风速数据作为输出数据，构建样本集；其中/>表示t时刻数值天气预报数据的第p个变量，x_k,t-i表示第t-i时刻测点k的实测阵风风速数据，i＝1,2,…,n；

第三步：采用特征选择方法比较第二步中样本集各输入变量的重要性程度，选取程度较高的若干变量作为最终输入，得到训练样本集；

第四步：采用移动窗口法，选取不同的预测步长，基于第三步中的训练样本集构建短期风速预测训练样本集和长期风速预测训练样本集；

第五步：基于第四步中的短期和长期风速预测训练样本集，分别对贝叶斯线性动态回归模型和贝叶斯二次动态回归模型进行训练，构建短期风速预测模型和长期风速预测模型；

第六步：基于第五步中的两种风速预测模型，实现多测点阵风风速的长短期预测；

第五步中贝叶斯线性动态回归模型表示为：

Y_t＝Φ(X_t)w_t+v_t,v_t～N(0,V)

其中，

Y_t＝[Y_1,t Y_2,t … Y_k,t]^T

X_t＝blockdiag(X_1,t X_2,t … … … X_k,t)

w_t＝[w_1,t w_2,t … … … w_k,t]^T

w_k,t＝[α_k,t β_k,t]^T＝[α_k,t β_k1,t β_k2,t … β_kd,t]^T

式中，Y_t表示目标区域内t时刻各测点的阵风风速预测结果构成的矩阵，Y_k,t表示目标区域内t时刻测点k的阵风风速预测值，blockdiag(·)表示块对角矩阵，X_t表示t时刻的输入数据，w_t表示t时刻的动态回归系数矩阵，w_k,t表示t时刻对应测点k的动态回归系数，α_k,t和β_k,t分别表示t时刻对应测点k的截距和回归系数，β_kd,t表示t时刻测点k的第d个回归系数，v_t表示服从0均值、方差V正态分布的回归误差；

动态回归系数的变化为：

w_t＝w_t-1+q_t,q_t～N(0,Q)

其中，w_t-1表示t-1时刻的动态回归系数矩阵，q_t表示服从0均值、方差Q正态分布的状态矩阵误差；

Q＝blockdiag(Q_1,t Q_2,t … … … Q_k,t)

Q_k,t＝blockdiag(Q_kα Q_kβ)

式中，Q_k,t表示t时刻测点k的状态误差，Q_kα和Q_kβ为t时刻对应测点k的截距和回归系数的状态误差；

第五步中，贝叶斯二次动态回归模型为：

Y_t＝[Y_1,t Y_2,t … Y_k,t]^T

α_t＝[α_1,t α_2,t … α_K,t]^T

式中，Y_t表示目标区域内t时刻各测点的阵风风速预测结果构成的矩阵，Y_k,t表示目标区域内t时刻测点k的阵风风速预测值，α_k,t和β_k,t分别表示t时刻对应测点k的截距和回归系数，X_l,t表示t时刻测点l的输入数据，β_kk,t表示t时刻测点k的平方系数，β_kl,t表示t时刻测点k和l的交叉系数，v_t表示服从0均值、方差V正态分布的回归误差，K表示目标区域内的测点数。

2.根据权利要求1所述的耦合数值天气预报与实测数据的风速多点同步预测方法，其特征在于，第五步中，采用卡尔曼滤波与期望最大化EM算法估计贝叶斯线性动态回归模型和贝叶斯二次动态回归模型参数。

3.根据权利要求1所述的耦合数值天气预报与实测数据的风速多点同步预测方法，其特征在于，第一步中，数值天气预报数据为欧洲中期天气预报中心数据ECMWF、或第五代大气研究中心中尺度模式MM5、或天气研究与预报模型WRF。

4.根据权利要求1所述的耦合数值天气预报与实测数据的风速多点同步预测方法，其特征在于，第三步中，选取数值天气预报数据中10m高度阵风风速、10m高度平均风速、温度以及目标测点对应大气压处的温度和风速。

5.根据权利要求1所述的耦合数值天气预报与实测数据的风速多点同步预测方法，其特征在于，第一步中的数值天气预报数据为目标测点所在位置处经纬度0.25°×0.25°范围内九个网格节点的数据。