CN108062595A - 基于wrf/cfd/sahde-rvm耦合的复杂地貌区域短时风能预测方法 - Google Patents

Abstract

发明提供基于WRF/CFD/SAHDE‑RVM耦合的复杂地貌区域短时风能预测方法。该预测方法包括WRF风向预测、SAHDE‑RVM短时风速预测、初步制作APG、完善APG和短时风能预测等步骤。该方法能够通过风速计的观测数据更加精确地匹配到对应的边界条件的信息。能够为复杂地貌区域的风电塔提供未来短时间内的风向和风速信息，便于区域内各个风电塔进行实时地调整，提高咬风能力。具有广阔的应用前景。

Description

基于WRF/CFD/SAHDE-RVM耦合的复杂地貌区域短时风能预测 方法

技术领域

本发明涉及风电领域，特别涉及一种风能预测方法。

背景技术

为了对一个复杂地貌区域进行短时风能预测,现行的方法中有一种基于天气研究与预报模式(The Weather Research and Forecasting Model，WRF)提供的气象数据在WAsP软件中进行区域的风场短时预测，WAsP模型的基本思想是将控制气流运动的动量方程线性化，这样可以显著减少求解非线性偏微分方程(the nonlinear partialdifferential equations)的难度，同时也节约了大量的时间成本。另外有一种基于WRF提供的气象数据进行实时模拟预测，随后将模拟结果输入CFD软件进行计算，获取整片区域预测风场信息的方法。此外还有一种测量与预测联系的方法(Measure-Correlation-prediction，简称MCP)。

发明人发现现有技术至少存在以下缺陷：

一是采用WAsP模型,是通过风速计观测资料计算地转风，假定一定范围内的地转风速不会变化，再根据周围地区的地表特征，由地转风推算任一点的近地层风速和风能参数。由于是从一个观测点外推，误差会随着与测风点距离的增加而加大，因此直接将WAsP用于较大区域的风能推算不可行；另外，由于WAsP本身采用线性模型进行计算，会造成随着流经地形的复杂流体的计算结果的不确定性，有其一定的局限性；同时现有的研究发现表明采用线性的模型会高估山顶处气流的加速度而且无法解决山体背风面的气体分流问题。

二是采用CFD软件实时进行数值模拟的方法，需要及时从WRF中提取短时预测气象数据，同时运行所耗费的计算资源巨大，需要一定的计算时间，效率难以保证，不便于风电塔根据短时风能预测结果实时进行应对调整。

三是采用MCP方法需要两个必要条件：(1)需要获得复杂地貌区域全局和局部一段时间内的观测数据，(2)该全局数据与局部数据要有好的相关性。但由于复杂地貌区域下的流体特性的非线性，目前很难去找到一个特定的方程去描述它们的相关性；同时，进行全局的实地数据测量需要设立多个风速计对整个地区进行长时间的风速及风向的观测，这种方法首先建设及维护成本高昂，其次个别测量地点位置的不当选择可能会影响复杂地形区域风场信息的表达。

因此，提供基于WRF/CFD/SAHDE-RVM耦合的复杂地貌区域短时风能预测方法,对复杂地形区域的风能评估有重要的意义。

发明内容

本发明的目的是提供一种可根据实际状况计算合理限速值的桥上行车限速控制***，以解决现有技术中存在的问题。

为实现本发明目的而采用的技术方案是这样的，基于WRF/CFD/SAHDE-RVM耦合的复杂地貌区域短时风能预测方法，包括以下步骤：

1)获取目标复杂地貌区域地形资料，基于WRF提供的气象数据进行短时风向预测。得出目标复杂地貌区域的短时主流风向。

2)SAHDE-RVM短时风速预测。

2.1)使用RVM进行回归预测，建立RVM回归预测模型。

2.2)收集目标复杂地貌区域风速计处历史风速风向数据作为SAHDE算法的训练样本。通过对样本数据集的训练，自适应获取步骤2.1)所述RVM回归预测模型的最优参数。

2.3)根据风速计实时数据，返回风速计处的短时风速预测结果。

3)在CFD软件中，均匀进行若干个风向角下目标复杂地貌的数值模拟计算。根据CFD结果制作初始APG。

4)根据风速计的长期实测数据，通过算法完善APG。

5)将步骤1)所述的WRF的风向预测结果与步骤4)所述APG进行匹配，提取出对应的边界速度信息和流场信息。获得整片区域风能的短时预测结果。

进一步，步骤2.1)中所述RVM的回归预测模型为：

式中，w_i为加权系数，K(x,x_i)为核函数，N为样本数量。

进一步，所述RVM回归预测模型的核函数选用高斯核函数与二项式核函数组合得到的组合核函数。其中，所述组合核函数为：

K(x_i,x_j)＝λG(x_i,x_j)+(1-λ)P(x_i,x_j) (2)

式中，G(x_i,x_j)为高斯核函数，P(x_i,x_j)为二项式核函数，λ为权重。其中，0≤λ≤1。

进一步，所述步骤3)具体包括：

3.1)CFD软件获得目标复杂地貌区域的详细场地信息，运用前处理软件进行对应的模型建立，建立计算域并设定好边界条件。

3.2)将输入的边界速度无量纲化后分为x向及y向。其中，x向的速度用u表示，y向速度用v表示。边界处风速在不同风向角θ下记为u_B＝Vcosθ，v_B＝Vsinθ，B表示边界。计算无量纲化后的流场信息，用f(i,j,k)表示。

3.3)记录对应风速计侧风点在其观测高度下的风速数据(u_A,v_A)，以及对应全场地表面的风速信息。

3.4)将步骤3.2)中所述的无量纲化的边界速度和流场信息与步骤3.3)所述的无量纲条件下的风速计所在点(u_A,v_A)相关联。

3.5)根据不同入口风速下数值模拟获得的多个(u_A,v_A)结果，将其单位化得记为(U_A,V_A)。

3.6)将(U_A,V_A)点及其包含的边界速度信息和流场信息***APG中，绘制初始APG。其中，APG仅显示风速计实测的数据，但其中每个点都包含了与其对应的边界风速信息和流场信息。

本发明的技术效果是毋庸置疑的：

A.实现了所需计算资源的最小化；

B.能够通过风速计的观测数据更加精确地匹配到对应的边界条件的信息；

C.能够为复杂地貌区域的风电塔提供未来短时间内的风向和风速信息，便于区域内各个风电塔进行实时地调整，提高咬风能力；同时还能通过风速与发电功率的关系计算出短期内的发电量。

附图说明

图1为预测方法流程图；

图2为WRF短时风向预测示意图；

图3为APG的制作、完善与应用示意图；

图4为SAHDE-RVM短时风速预测示意图；

图5为根据CFD结果制作初始APG示意图；

图6为根据算法完善APG示意图；

图7为SAHDE-RVM短时风速预测示意图；

图8为WRF短时主导风向预测结果及观测数据对比图；

图9为SHADE-RVM短期风速预测结果与实测数据对比图；

图10为CFD软件中对复杂地貌区域模型的网格划分图；

图11为初始APG图；

图12为复杂地貌区域风能密度图。

具体实施方式

下面结合实施例对本发明作进一步说明，但不应该理解为本发明上述主题范围仅限于下述实施例。在不脱离本发明上述技术思想的情况下，根据本领域普通技术知识和惯用手段，做出各种替换和变更，均应包括在本发明的保护范围内。

实施例1：

本实施例公开一种基于WRF/CFD/SAHDE-RVM耦合的复杂地貌区域短时风能预测方法，参见图1，包括以下步骤：

1)参见图2，获取目标复杂地貌区域地形资料，并对其划分多重嵌套网格。设置边界层参数方案。基于WRF提供的气象数据进行短时风向预测。模拟未来一段时间内水平面内两个相互垂直方向的风速，每隔一段时间输出一次数据，用临近的九个网格点的均值作为风速大小的输出结果，参见图8，通过两个垂直方向速度的合成即可得到未来一段时间的该复杂地貌区域的短时主流风向。其中，所述气象数据包括同化全球环流模式再分析资料、常规地面气象资料和探空气象资料。

2)SAHDE-RVM短时风速预测。

2.1)使用RVM进行回归预测，建立RVM回归预测模型。其中，给定训练样本的输入集X＝{x₁,x₂,x₃,…,x_n,}与相应的输出集T＝{t₁,t₂,t₃,…,t_n}，其中n为样本个数，设t_i为目标值且有t_i∈R。输出值t_i的函数模型可表示为：

t_i＝y(x_i,w)+ε_i (1)

式中，ε_i表示高斯白色噪声，且ε_i服从分布ε_i～N(0,σ²)，则p(t_i|x_i)＝N(t_i|y(x_i,w),σ²)。RVM模型的输出可表示为非线核函数的组合，核函数无需满足mercer条件。RVM的回归预测模型为：

式中，w_i为加权系数，K(x,x_i)为核函数，N为样本数量。

对于独立分布的输出集t_i的似然估计为：

式中，Φ＝(Φ₁,Φ₂,…,Φ_N)，Φ_i＝(1,K(x_i,x₁),...,K(x_i,x_n))^T(i＝1,…,N)。

由稀疏贝叶斯原理定义的权值参数为零的高斯先验分布为：

式中，α_i为先验高斯分布的超参数，α＝(α₀,…,α_N)^T。每个独立的超参数α_i控制着权参数w_i的先验分布，使相关向量机模型具有稀疏性。

根据贝叶斯原理，由式(3)和(4)计算权值矢量w的后验分布：

式中，

由式(5)可知，若要确定权值矢量w需对超参数α、σ²进行确定。

采用贝叶斯框架计算超参数的似然分布：

p(t|α,σ²)＝∫p(^t|w,σ²)p(w|α)dW＝N(0,C)(6)

式中，C为协方差且C＝σ²I+ΦA^-1Φ^T。

对式(6)中的两个超参数α、σ²分别求偏导数且令偏导数等于零。

得到求解最优超参数过程中的更新迭代算法：

式中，u_j是第j个后验平均权值，γ_j＝1-α_jM_jj(γ_j∈[0,1])，M_jj为权值协方差矩阵Σ的对角线元素。

RVM学习过程中，首先初始化式(7)、(8)中的两个参数，并通过更新迭代公式不断跟新计算两个超参数，当模型中的参数均达到最大训练次数时，更新计算停止，此时得到的α与σ²为最优值。

若给定***一个输入值X_*，则输出的概率分布为：

y_*＝u^TΦ(X_*)＝u_MPΦ(X_*) (10)

式中，当给定一个输入值X_*，则有相应的预测输出值y_*。为预测方差，表示不确定性。

2.2)参见图7，收集目标复杂地貌区域风速计处历史风速风向数据作为自适应混合差分进化算法(simulated annealing-based hybrid differential evolution，SAHDE)的训练样本。通过对样本数据集的训练，自适应获取步骤2.1)所述RVM回归预测模型的最优参数。其中，在本实施例中，SAHDE操作步骤如下：

a)变异操作：DE算法由多种变异操作方法，选用如下方法进行变异操作。

x_m＝x_best+F[(x₁-x₂)+(x₃-x₄)] (12)

随机挑选出四个父代个体，分别记为x₁,x₂,x₃,x₄，x_best是父代中的最优个体，x_m为变异产生的变异个体，F为变异率(F∈[0,1.2])。

b)交叉操作：选取两个个体x_i和x_m进行交叉操作，交叉操作后生成的新个体为x_c，具体的操作方法如下式：

其中，rand()表示[0,1]之间的随机函数，randr(i)∈{1,2,…,D}为随机产生的整数，D为优化变量的维数，C_R为交叉率(C_R∈(0,1))。

c)选择操作：DE的子代个体由贪婪选择法产生，经过上述交叉与变异的操作后，产生了个体x_c和目标个体二者进行竞争，较优者被选为子代个体，选择操作法如式11所示：

式中，x_c为交叉操作与变异操作后产生的个体，为子代中第i个个体。

SAHDE算法在初始时进行全局搜索而在后期进行局部搜索，C_R的大小由当前进化代数g_now与最大进化代数g_max自适应调整，具体方法如下所示：

C_R0为交叉算子C_R的初值，C_R的值根据上式自适应调整，初始值C_R0较小，而后其取值逐步增大。

使用模拟退火算法对SAHDE产生的当前最优个体进行二次搜索。在SAHDE当中，选定当前最优个体为初始个体，既y₀＝x_best，初始温度选定为T₀,产生新个体的方式如下所示：

y_r+1,j＝y_r,j+ηε(x_jmax-x_jmin),j＝1,2,...,D (16)

式中，r为模拟退火算法的迭代次数。y_r为r次迭代后产生的新个体。η为控制扰动幅度。ε为服从均值或正态分布的随机变量。x_jmax、x_jmin分别为第j维优化变量的取值范围。

SAHDE产生的最优个体再经过模拟退火算法二次搜索后，适应度变化为ΔF，ΔF＝f(y_r+1)-f(y_r)。若ΔF＜0则接受新个体并用新个体替换原来的最优个体。若e^(-ΔF/T)＞rand()，此时也接受新产生的个体，并用其替换种群中的一个非最优个体。否则拒绝。若选择了接受新的个体，按T_r+1＝aT_r(0＜a＜1)降温。否则不降温。

为了加快样本的训练速度并提高模型的预测精度，本实施例中使用归一化的方法对风速实测样本数据进行预处理：

式中，为数据归一化的值。x_max为风速样本数据的最大值。x_min为风速样本数据的最小值。

选用高斯核函数与二项式核函数组合得到的组合核函数作为步骤2.1)所述RVM模型的核函数，取长补短，发挥各自的优势。组合后的核函数如下：

K(x_i,x_j)＝λG(x_i,x_j)+(1-λ)P(x_i,x_j) (18)

其中，G(x_i,x_j)为高斯核函数，其表达式为：G(x_i,x_j)＝exp(-||x_i-x_j||²/σ²)。P(x_i,x_j)为二项式核函数。λ为权重，0≤λ≤1。

收集风速计处过去一年的风速风向数据作为上述SAHDE-RVM算法的训练样本，通过对样本数据集的训练，自适应获取模型的最优参数α、σ²。

2.3)参见图4，根据风速计实时数据，返回风速计处的短时风速预测结果。此后每对其输入一个风速计处的瞬时风速数据都能返回一个风速计处的短时(指定时间，如15min)风速预测结果。附图9为100个测点处SHADE-RVM短期风速预测结果与实测数据的对比图。

3)参见图5，在CFD软件中，均匀进行多个风向角下目标复杂地貌的数值模拟计算。根据CFD结果制作初始风速计相位图(Anemometer phase graph，APG)。其中，CFD软件中对复杂地貌区域模型的网格划分图如图10所示。初始APG图如图11所示。

3.1)CFD软件获得目标复杂地貌区域的详细场地信息，运用前处理软件(如ICEM)进行对应的模型建立，同时建立与之相匹配的计算域，在划分网格后导入FLUENT，设定好边界条件各个边界如何设置后进行数值模拟计算。

3.2)入口风速V考虑参考高度z_ref＝10m处以10m/s进行一次该场地均匀划分的360个风向角下的数值模拟。记录下对应风速计侧风点记为A点在其观测高度下的风速数据u_A,v_A，及对应全场地表面的风速信息。将输入的边界速度无量纲化后分为x向及y向。其中，x向的速度用u表示，y向速度用v表示。边界处风速在不同风向角θ下记为u_B＝Vcosθ，v_B＝Vsinθ，B表示边界。计算无量纲化后的流场信息，用f(i,j,k)表示。

3.3)记录对应风速计侧风点在其观测高度下的风速数据(u_A,v_A)，以及对应全场地表面的风速信息。

3.4)将步骤3.2)中所述的无量纲化的边界速度和流场信息与步骤3.3)所述的无量纲条件下的风速计所在点(u_A,v_A)相关联。

3.5)根据不同入口风速下数值模拟获得的360个(u_A,v_A)结果，将其单位化得记为(U_A,V_A)。

3.6)将(U_A,V_A)点及其包含的边界速度信息和流场信息***APG中，绘制APG。其中，APG仅显示风速计实测的数据，但其中每个点都包含了与其对应的边界风速信息和流场信息。

4)根据风速计的长期实测数据，通过算法完善APG。

参见图6，根据风速计的长期实测数据将APG完善。假设u_*,v_*为风速计测得的某一个数据，单位化后得U_*,V_*，定义其在APG中的点为点P_*。如果在APG中，点P_*足够接近某一点T，则用点T作为P*的近似值，开始进行下一个风速计数据的流程。如果不接近则进行以下流程：将点P_*投影到离其最近的APG曲线中，投影点设为P，定义该曲线段的两个端点为T₁和T₂，计算P_*在该段的相对位置：根据r可对点P_*包含在APG图中的边界风速信息进行计算：u_B’＝(1-r)U_B1+rU_B2，v_B’＝(1-r)V_B1+rV_B2，U_B1,V_B1U_B2,V_B2分别为点T₁、T₂处对应的边界速度信息。随后通过u_B’,v_B’计算出流场的信息f’i,j,k，并将其***到风速计数据u_A’,v_A’中。u_A’,v_A’的模设为|u_A|，将各个上述求得的量单位化即可得 将点U_A’,V_A’及其相关联的边界速度信息U_B’,V_B’及流场信息F’(i,j,k)***到APG中。

随后再次判断点P_*是否足够接近某一点T，重复上述流程。

随着越来越多的信息被风速计收集及录入，APG中会变得越来越完善，能够提供更全面的风速计数据及边界速度之间的相关信息。

5)参见图3，将步骤1)所述的WRF的风向预测结果与步骤4)所述APG进行匹配，提取出对应的边界速度信息和流场信息。获得整片区域风能的短时预测结果。由预测结果得到的复杂地貌区域风能密度图如图12所示。将WRF的风向预测结果风速计处数据与绘制出的APG进行匹配，提取出对应的无量纲化的边界速度信息以及流场信息，根据同一个风向角不同参考风速情况下流场的形状是线性相关这一点，利用SAHDE-RVM的风速大小预测结果，与无量纲化的边界速度信息及流场信息相乘即可获得整片区域风能的短时预测结果。如SAHDE-RVM预测短时风速为u_C,v_C，根据WRF提供的风向预测结果对应的APG中的边界速度信息U_BC,V_BC及流场信息F_C(i,j,k)，计算相应的边界速度 此时即可得到整片区域任一点的短时风能预测结果。

值得说明的是，本实施例展现了中尺度数值天气预报模式与小尺度数值模式全新的结合思路，通过CFD制作初始APG之后用算法完善APG的方法，实现了所需计算资源的最小化，不需再通过导入WRF预测数据后实时利用CFD对整片区域进行数值模拟计算。同时由于中尺度模式中的网点计算值是对于一个网格大小内的区域一段时间内的平均值，无法具体代表某一确定地点的风速情况，为此我们引入了SAHDE-RVM。通过使用区域中的单个风速计的历史数据对SAHDE-RVM进行训练，实现了相对中尺度WRF模式更短时间内的风速预测。

与线性WAsP模型是一个根据以往数据作出的回归结果再外推不同，本实施例包含有基于计算流体动力学的物理模型进行的数值模拟计算，能够更广阔的在复杂地貌区域应用。同时本实例不仅不需要进行短时风能预测地点的历史风速、风向数据，只需获得整片复杂地貌区域的高程数据，结合WRF和CFD技术，就能通过在该区域设立一个风速计获得的观测数据预测出整个复杂地貌区域一段时间后的风速风向状态。在短时风速预测过程中，对比传统DE算法，本方法引入了模拟退火算法，对SAHDE产生的最优个体进行二次搜索，摆脱了传统DE算法的早熟收敛的缺点。此外，当上述方法真正应用到了风能评估中去时，APG将会随着观测时间的增长逐渐变成一个大数据集，使其与传统MCP方法提供一个固定的相关方程相比，能够通过风速计的观测数据更加精确地匹配到对应的边界条件的信息。

将此方法运用到风电相关领域中，从长期看，能够对整片复杂地貌区域各个地点常吹风向及风速大小进行汇总，综合各方数据后分析能够得出该地区最适宜设立风力发电塔的地点及发电塔的朝向。从短期来看，能够为复杂地貌区域的风电塔提供未来短时间内的风向和风速信息，便于区域内各个风电塔进行实时地调整，提高咬风能力。同时还能通过风速与发电功率的关系计算出短期内的发电量。

Claims (4)

1.基于WRF/CFD/SAHDE-RVM耦合的复杂地貌区域短时风能预测方法，其特征在于，包括以下步骤：

1)获取目标复杂地貌区域地形资料，基于WRF提供的气象数据进行短时风向预测；得出目标复杂地貌区域的短时主流风向；

2)SAHDE-RVM短时风速预测；

2.1)使用RVM进行回归预测，建立RVM回归预测模型；

2.2)收集目标复杂地貌区域风速计处历史风速风向数据作为SAHDE算法的训练样本；通过对样本数据集的训练，自适应获取步骤2.1)所述RVM回归预测模型的最优参数；

2.3)根据风速计实时数据，返回风速计处的短时风速预测结果；

3)在CFD软件中，均匀进行若干个风向角下目标复杂地貌的数值模拟计算；根据CFD结果制作初始APG；

4)根据风速计的长期实测数据，通过算法完善APG；

5)将步骤1)所述的WRF的风向预测结果与步骤4)所述APG进行匹配，提取出对应的边界速度信息和流场信息；获得整片区域风能的短时预测结果。

2.根据权利要求1所述的基于WRF/CFD/SAHDE-RVM耦合的复杂地貌区域短时风能预测方法，其特征在于：步骤2.1)中所述RVM的回归预测模型为：

<mrow> <mi>y</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>w</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msub> <mi>w</mi> <mi>i</mi> </msub> <mi>K</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>w</mi> <mn>0</mn> </msub> <mo>+</mo> <msub> <mi>&amp;epsiv;</mi> <mi>i</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>

式中，w_i为加权系数，K(x,x_i)为核函数，N为样本数量。

3.根据权利要求2所述的基于WRF/CFD/SAHDE-RVM耦合的复杂地貌区域短时风能预测方法，其特征在于：所述RVM回归预测模型的核函数选用高斯核函数与二项式核函数组合得到的组合核函数；其中，所述组合核函数为：

K(x_i,x_j)＝λG(x_i,x_j)+(1-λ)P(x_i,x_j) (2)

式中，G(x_i,x_j)为高斯核函数，P(x_i,x_j)为二项式核函数，λ为权重；其中，0≤λ≤1。

4.根据权利要求1所述的基于WRF/CFD/SAHDE-RVM耦合的复杂地貌区域短时风能预测方法，其特征在于，所述步骤3)具体包括：

3.1)CFD软件获得目标复杂地貌区域的详细场地信息，运用前处理软件进行对应的模型建立，建立计算域并设定好边界条件；

3.2)将输入的边界速度无量纲化后分为x向及y向；其中，x向的速度用u表示，y向速度用v表示；边界处风速在不同风向角θ下记为u_B＝Vcosθ，v_B＝Vsinθ，B表示边界；计算无量纲化后的流场信息，用f(i,j,k)表示；

3.3)记录对应风速计侧风点在其观测高度下的风速数据(u_A,v_A)，以及对应全场地表面的风速信息；

3.4)将步骤3.2)中所述的无量纲化的边界速度和流场信息与步骤3.3)所述的无量纲条件下的风速计所在点(u_A,v_A)相关联；

3.5)根据不同入口风速下数值模拟获得的多个(u_A,v_A)结果，将其单位化得记为(U_A,V_A)；

3.6)将(U_A,V_A)点及其包含的边界速度信息和流场信息***APG中，绘制初始APG。其中，APG仅显示风速计实测的数据，但其中每个点都包含了与其对应的边界风速信息和流场信息。