具体实施方式

为使本发明的目的、内容和优点更加清楚，下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。

本发明涉及环境风场分析领域，具体是通过技术研究，实现测风塔多层风向风速概率的计算方法，为风能资源评估、工程结构设计、浅层风风场预报提供支撑。

有鉴于此，本发明通过用gamma分布拟合风速，混合von mises分布拟合风向，得到风向、风速概率密度函数。接着通过计算风向、风速相关系数，风向、风速累积概率密度，得到风向风速联合概率密度。在风向风速概率密度数据基础上，通过matplotlib绘制了多层风场的概率密度图形，得到多层风场风向风速的概率分布模型。该方法适用于风能资源评估、工程结构设计、浅层风风场预报。

为了实现上述目的本发明采用的技术方案如下：测风塔多层风向风速的概率分布实现方法，包括以下步骤：

S1、逐层风向、风速概率分布：

数据处理，舍弃静风场。

拟合风速概率分布。

拟合风向概率分布。

S2、风向风速联合概率分布：

计算风向、风速累积概率。

计算风向、风速相关系数。

计算风向风速联合概率。

S3、多层风向风速概率分布模型：

各层风向风速联合概率数据处理。

多层风向风速概率分布。

多层风向风速概率距平分布。

本发明的方法具体包括如下步骤：

S1、逐层风向、风速概率分布：

S11、数据处理

将分钟风场数据按照逐小时进行聚合，分别得到10米、50米、90米逐小时最大风速及对应风向。小于0.5m/s的风规定为静风场，对静风场数据进行舍弃。

S12、用gamma分布来拟合各层风速数据

gamma分布的概率密度函数为：

其中v为风速，η为形状参数，γ为尺度参数。

对各层风速数据，根据极大似然估计方法计算gamma分布的形状参数、位置参数、尺度参数，得到相对应的风速分布概率密度函数。具体由python的科学计算库scipy实现。

S13、用混合von mises分布来拟合各层风向数据

混合von mises分布的概率密度函数如下：

其中，θ为风向，i为von mises分布混合组数，w_i为各组分布所占权重；μ_i为vonmises分布的位置参数，k_i为尺度参数；I₀(k_i)为修正的零阶第一类Bessel 函数。

混合组数设置为8。将0-360°等间隔划成8方位，各方位阈值分别为 (45°,90°,135°,180°,225°,270°,315°,360°)。获得风向数据在各方位的概率累积值，即权重w。

采用极大似然估计方法计算各组von mises分布的参数。计算各组对应的平均正弦值和余弦值进而得到各组分布的位置参数μ。

在平均正弦值和余弦值基础上，根据经验函数得到尺度参数k。

如果各组权重w有接近或等于0的情况，去除这一组，剩下的组数再进行 w、μ、k的迭代计算。

S2、逐层风向风速联合概率分布：

S21、计算风向、风速累积分布概率

在概率密度函数基础上，积分计算风向、风速的累积概率密度F(θ)和F(υ)。风向累积概率密度F(θ)的积分步长为1°，风速累积概率密度F(υ)的积分步长为0.5m/s。

S22、计算风向、风速相关系数

对风向数据进行正弦值和余弦值的计算。风向风速相关系数计算公式如下：

其中r_vc为各风速与所对应风向余弦值的相关系数，r_vs为各风速与所对应风向正弦值的相关系数，r_cs为风向余弦值和正弦值之间的相关系数。风向风速相关系数为r。

S23、风向风速联合概率分布函数

风向风速联合概率分布函数如下：

f(ν,θ)＝f(ν)f(θ)[1+δ(1-2F(υ)(1-2F(θ))]

ν为风速，θ为风向；f(ν)为风速概率密度函数，f(θ)为风向概率密度函数；r 为风向风速相关系数；F(υ)为风速累积密度函数，F(θ)为风向累积密度函数。将相应要素和函数带入公式，就得到风向风速联合概率分布的数据集。

S3、多层风向风速概率分布模型：

S31、各层风向风速联合概率数据处理

获取10米、50米、90米各层的风向风速联合概率f(ν,θ)的最大值，将这个最大值作为多层建模Z轴方向的间隔Δh。

S32、多层风向风速概率分布

以10米、50米、90米各层的风向风速联合概率f(ν,θ)为输入，利用python 的matplotlib进行三维绘图，视角中的仰角设置为10°、方位角设置为210°，绘图颜色设置为“rainbow”，利用plot_surface函数进行三维曲面绘图。10 米绘制在第一层，然后间隔Δh绘制50米为第二层，再间隔Δh绘制第三层即90 米。得到多层风向风速联合概率。

S33、多层风向风速概率距平分布

对10米、50米、90米风向风速联合概率f(ν,θ)取平均，得到联合概率平均值。然后再分别用10米、50米、90米联合概率减去概率平均值，得到各层联合概率距平。重复1)2)步骤，绘图颜色设置为“seismic”(负距平为冷色调蓝色，正距平为暖色调红色)，得到多层风场的风向风速联合概率的距平分布。

具体实施例如下：

测风塔多层风向风速概率分布实现方法主要包括：拟合得到风向、风速概率分布；风向风速联合概率分布；多层风向风速概率分布模型建立。

S1、逐层风向、风速概率分布：

S11、数据处理

测风塔高100米，其中在10米、50米、90米安装的有风速、风向传感器，通过测量可以得到逐分钟风场数据。将分钟风场数据按照逐小时进行聚合，得到10米、50米、90米逐小时最大风速及对应风向。小于0.5m/s的风规定为静风场，对静风场数据进行舍弃。

S12、用gamma分布来拟合各层风速数据

gamma分布的概率密度函数为：

其中v为风速，η为形状参数，γ为尺度参数。

对各层风速数据，根据极大似然估计方法计算gamma分布的形状参数、尺度参数，得到相对应的风速分布概率密度函数。具体实现由python的科学计算库scipy实现(scipy.stats.gamma.fit())。

S13、用混合von mises分布来拟合各层风向数据

混合von mises分布的概率密度函数为：

其中，θ为风向，i为von mises分布混合组数，w_i为各组分布所占权重；μ_i为vonmises分布的位置参数，k_i为尺度参数；I₀(k_i)为修正的零阶第一类Bessel 函数。

采用迭代方式计算混合von mises分布的概率密度函数。

混合组数可自行定义，考虑到风向实际情况，一般不宜过大。根据风向常用方位，这里将混合组数设置为8。即将0-360°等间隔划成8方位，各方位阈值分别为(45°,90°,135°,180°,225°,270°,315°,360°)。计算各方位风向数据样本占总风向样本的比例，视为风向数据在各方位的概率累积值，即权重w。

采用极大似然估计方法计算各组von mises分布的参数。计算各组风向数据对应的平均正弦值和余弦值进而得到各组分布的位置参数μ。

在平均正弦值和余弦值基础上，根据经验函数得到尺度参数k。

如果各组权重w有接近或等于0的情况，去除这一组，剩下的组数再进行 w、μ、k的迭代计算。

S2、逐层风向风速联合概率分布：

S21、计算风向、风速累积分布概率

在概率密度函数基础上，积分计算风向、风速的累积概率密度F(θ)和F(υ)。风向累积概率密度F(θ)的积分步长为1°，风速累积概率密度F(υ)的积分步长为0.5m/s。

S22、计算风向、风速相关系数

对所有的风向数据一一进行正弦值和余弦值的计算，得到各风向数据的正弦和余弦值。然后计算风向风速相关系数。

风向风速相关系数计算公式如下：

其中r_vc为各风速与所对应风向余弦值的相关系数，r_vs为各风速与所对应风向正弦值的相关系数，r_cs为风向余弦值和正弦值之间的相关系数。风向风速相关系数为r。

S23、风向风速联合概率分布函数

风向风速联合概率分布函数如下：

f(ν,θ)＝f(ν)f(θ)[1+r(1-2F(υ)(1-2F(θ))]

ν为风速，θ为风向；f(ν)为风速概率密度函数，f(θ)为风向概率密度函数；r 为风向风速相关系数；F(υ)为风速累积密度函数，F(θ)为风向累积密度函数。将相应要素和函数带入公式，就得到风向风速联合概率分布的数据集。

S3、多层风向风速概率分布模型：

S31、各层风向风速联合概率数据处理

获取10米、50米、90米各层的风向风速联合概率f(ν,θ)的最大值，将这个最大值作为多层风向风速概率分布建模中Z轴方向的各层显示间隔Δh，便于概率分布模型的展示。

S32、多层风向风速概率分布

以10米、50米、90米各层的风向风速联合概率f(ν,θ)为输入，利用python 的matplotlib进行三维绘图，视角中的仰角设置为10°、方位角设置为210°，绘图颜色设置为“rainbow”，利用plot_surface函数进行三维曲面绘图。10 米绘制在第一层，然后间隔Δh绘制50米为第二层，再间隔Δh绘制第三层即90 米。得到多层风向风速联合概率。图示可见附图1。

S33、多层风向风速概率距平分布

对10米、50米、90米三层的风向风速联合概率f(ν,θ)取平均，得到联合概率平均值然后再分别用10米、50米、90米联合概率减去概率平均值，得到各层联合概率距平。重复S32步骤，利用python的matplotlib进行三维绘图，利用plot_surface函数进行三维曲面绘图,10米绘制在第一层，然后间隔Δh绘制50米为第二层，再间隔Δh绘制第三层即90米,绘图颜色设置为“seismic”(负距平为冷色调蓝色，正距平为暖色调红色)，得到多层风场的风向风速联合概率的距平分布。图示可见附图2。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。