CN108062722B - 基于风速变异系数的山地风电场模型风机机械功率计算方法 - Google Patents

Abstract

本发明公开了一种基于风速变异系数的山地风电场模型风机机械功率计算方法，包括以下步骤：步骤1、根据风能数据，构建山地风电场风速时空特征参数表；步骤2、选择某一工况S，获取风速非一致性指标（风速变异系数）σ _s 及其置信水平CL，进而计算变异系数上分位点σ _u 和下分位点σ _l ；步骤3、结合山地风电场风速时空特征参数表该工况的风速空间分布类型WSD，分别计算该σs置信区间对应的各机组风速；步骤4、选择风速时间分布类型WTD，并确定相关参数；步骤5、分别计算该σ _s 置信区间对应的机组机械输出功率。本发明法解决了现有技术中反映山地风电场风速非一致性指标缺失，无法量化评价和使用风速非一致性空间分布特征等技术问题。

Description

基于风速变异系数的山地风电场模型风机机械功率计算方法

技术领域

本发明属于风力发电技术领域，具体涉及一种基于风速变异系数的山地风电场模型风机机械功率计算方法。

背景技术

随着我国环境问题的凸显，风力发电作为解决环境问题的有效手段之一，近十年来不断增长，并已成为我国继煤电、水电之后的第三大电源。截至2016年底，我国风电装机已达10万余台，累计装机容量达到1.69亿千瓦。根据我国风电规划目标，到2020年底，风电累计年发电量确保达到4200亿千瓦时，约占全国总发电量的6％，可以说风电正在成为我国能源结构的重要组成部分。但风电出力具有随机性和间歇性以及反调峰特性，给电力系统安全稳定运行提出了一系列新的问题，引起了国内外电力与能源领域工业界专家和学术界学者的广泛关注和深入研究。据统计，我国可利用的低风速资源面积约占全国风能资源区的68％。近来，随着高塔筒、大叶轮等技术创新与实践，我国内陆风电重点开发区域由限电严重的高风速“三北地区”逐渐转向接近电网负荷受端地区，不存在远距离送电问题的低风速“华南、西南、华东地区”，但此类地区海拔高、地形和气象条件复杂，山地风电场运行特性具有较为明显的差异性特征。其中，贵州特有的高原山区气象性能决定了风能所处环境有着海拔高、湿度大、地形复杂等特点，山地风电场明显区别于我国其他区域的平原、海上风电场。以贵州来看，预计2017年底，风电装机规模将达到4490兆瓦，占全网装机比例达到9.8％，提升2％；全年发电量67亿千瓦时，占全网发电量4.2％。可以预见随着风电装机比例进一步增长，将对贵州电网为代表的交直流送端电力系统安全稳定运行带来新的挑战。因此，亟需面向云贵高原地理和气候环境特点，开展山地风电场机电暂态模型及其建模方法研究，更好地掌握高比例风能接入电力系统的动态特性，从而为类似环境地区电网规划建设和运行调度关键技术进一步奠定基础。

现有技术中，风电场原动机模型侧风速模拟上一般采用在一个时段内单一机组、测风塔甚至整个风电场所有机组的平均风速，导致应用于山地风电场时，模型输出功率特性准确性不高，难以适应含复杂地形山地风电场的电力系统动态特性分析的需要；此外，近期有研究提出考虑弱一致性风速分布山地风电场的机电暂态模型构建，但却没有提出关于该弱一致性或非一致性的量化评价指标，故无法量化评价和使用山地风电场的风速非一致性空间分布特征。

发明内容

本发明的目的在于克服上述缺点，提出一种能量化评价和使用山地风电场的风速非一致性空间分布特征，有益于风速空间分布的量化模拟，提高对于弱一致性分布山地风电场接入系统动态分析的准确性的基于风速变异系数的山地风电场模型风机机械功率计算方法。

本发明的一种基于风速变异系数的山地风电场模型风机机械功率计算方法，包括以下步骤：

步骤1、根据风能历史数据或评估数据，构建山地风电场风速时空特征参数表，参数表的参数包括工况S、机组数N、机组编号序列W_N、平均风速V_av、风速非一致性指标(风速变异系数)σ_s及其置信区间CL，风速空间分布类型WSD，风速时间分布类型WTD；

步骤2、由步骤1的山地风电场风速时空特征参数表，选择某一工况S，获取风速非一致性指标(风速变异系数)σ_s及其置信水平CL，进而计算变异系数上分位点σ_u和下分位点σ_l；

步骤3、针对所选择工况S，依据步骤2得到的上分位点σ_u、下分位点σ_l，结合山地风电场风速时空特征参数表中该工况的风速空间分布类型WSD，分别计算该σ_s置信区间对应的各机组2组风速V_i，i＝1,2,…,N，V_j，j＝1,2,…,N；

步骤4、针对所选择工况S，由山地风电场风速时空特征参数表，选择风速时间分布类型WTD，并确定参数；

步骤5、针对所选择工况S，计算该山地风电场风能利用系数C_p、折减因子α_TM-H、以及第3步骤所得各机组风速V_i和V_j，分别计算该σ_s置信区间对应的2组各机组机械输出功率P_{mech_i}，i＝1,2,…,N，P_{mech_j}，j＝1,2,…,N。

上述基于风速变异系数的山地风电场模型风机机械功率计算方法，其中：所述步骤1的山地风电场包括现有山地风电场接入系统分析场景和规划山地风电场接入系统分析场景。

上述基于风速变异系数的山地风电场模型风机机械功率计算方法，其中：所述现有山地风电场接入系统的运行方式分析场景，根据风能历史数据，借助数理统计和聚类方法，获得主要工况S集合，以及所对应的平均风速V_av、风速非一致性指标(风速变异系数)σ_s及其置信区间CL，风速空间分布类型WSD，风速时间分布类型WTD，其中：

(1)给定机组数N下，风电场风速平均值V_av的计算公式为

(2)风电场风速标准差V_δ的计算公式为

(3)提出以变异系数σ_s作为衡量风速非一致性的指标，其计算公式为

其中，V_δ为风电场风速标准差，V_av为风电场风速平均值；

(4)现有山地风电场的变异系数置信区间CL，可取为，CL＝0.99或0.95；

(5)风速空间分布类型WSD，涉及泊松分布、均匀分布、指数分布、正态分布、瑞利分布、韦布尔分布；

(6)风速时间分布类型WTD，包括斜坡风、阵风、“墨西哥草帽风”，以及自定义分段线性函数风速。

上述基于风速变异系数的山地风电场模型风机机械功率计算方法，其中：所述规划山地风电场接入系统的预想事件分析场景，根据风能评估数据，借助数理统计和聚类方法，获得主要工况S集合，以及所对应的平均风速V_av、风速非一致性指标(风速变异系数)σ_s及其置信区间CL，风速空间分布类型WSD，风速时间分布类型WTD，其中：

(1)给定机组数N下，风电场风速平均值V_av的计算公式为

(2)风电场风速标准差V_δ的计算公式为

(3)提出以变异系数σ_s作为衡量风速非一致性的指标，其计算公式为

其中，V_δ为风电场风速标准差，V_av为风电场风速平均值；

(4)规划山地风电场的变异系数置信区间CL，可取为，CL＝0.95或0.9；

(5)风速空间分布类型WSD，涉及泊松分布、均匀分布、指数分布、正态分布、瑞利分布、韦布尔分布；

(6)风速时间分布类型WTD，包括斜坡风、阵风、“墨西哥草帽风”，以及自定义分段线性函数风速。

上述基于风速变异系数的山地风电场模型风机机械功率计算方法，其中：所述步骤2中由步骤1的山地风电场风速时空特征参数表，选择某一工况S，获取风速非一致性指标(风速变异系数)σ_s及其置信水平CL，进而计算变异系数上分位点σ_u和下分位点σ_l，其具体步骤包括：

步骤2.1、依据查表法，由Z检验表可确定置信水平下的Z_cl系数；

步骤2.2、计算相应置信区间中变异系数下分位点σ_l：

σ_l＝σ_s-Z_cl×σ_d`

以及相应置信区间中变异系数上分位点σ_u：

σ_u＝σ_s+Z_cl×σ_d

其中，变异系数标准差σ_d可取为0.05，或根据实际数据样本计算得到。

上述基于风速变异系数的山地风电场模型风机机械功率计算方法，其中：所述步骤3中针对所选择工况S，依据步骤2得到的上分位点σ_u、下分位点σ_l，结合山地风电场风速时空特征参数表中该工况的风速空间分布类型WSD，分别计算该变异系数置信区间对应的各机组2组风速V_i，i＝1,2,…,N，V_j，j＝1,2,…,N,具体步骤如下：

步骤3.1、根据步骤2选择的工况S，结合步骤1山地风电场风速时空特征参数表中该工况下可选的风速空间分布类型，选择一个风速空间分布类型WSD，选择依据为：现有风电场或类似风电场典型工况或极端工况下的风速空间分布与上述可选风速空间分布类型WSD最为接近，即可选择；

步骤3.2、依据步骤2所得风速非一致性指标(风速变异系数)上分位点σ_u、下分位点σ_l，以及待分析工况的风速平均值V_av，分别可得该工况下变异系数上分位点和下分位点对应的风速标准差V_δu和V_δl。

步骤3.3、将上述参数代入相应概率分布函数中，得到对应于变异系数上下分位点的2组N个服从该风速空间分布类型WSD(正态分布)的伪随机数据，即该变异系数置信区间对应的2组N台机组的风速V_i,i＝1,2…,N，V_j，j＝1,2,…,N。

上述基于风速变异系数的山地风电场模型风机机械功率计算方法，其中：所述步骤4中针对所选择工况S，由山地风电场风速时空特征参数表，选择风速时间分布类型WTD，并确定参数，相关WTD模型及其参数如下：

(1)WTD＝1，即斜坡风模型，其数学模型如下：

式中，斜坡风起始时间t_sr，斜坡风结束时间t_er，风速增加幅值A_r，D_r＝t_er–t_sr；

(2)WTD＝2，即阵风模型，其数学模型如下：

式中，阵风起始时间t_sg，阵风结束时间t_eg，最大风速V_max，A_g＝(V_max–V₀)/2，D_g＝t_eg–t_sg；

(3)WTD＝3，即“墨西哥草帽风”模型，其数学模型如下：

式中，V₀为初始风速、V_max为最大风速、V_min为最小风速，草帽风起始时间t₀，草帽风结束时间t₃，草帽上升沿时间t₁，草帽下降沿时间t₂；

(4)WTD＝4，即一种自定义分段线性函数风速模型，其数学模型如下：

式中，V₀为初始时刻t₀风速，V₁、V₂、V₃、V₄分别为中间过程各时间点t₁、t₂、t₃、t₄的风速，V₅为结束时刻t₅的风速，其中，V₁＝K₁V₀，V₂＝K₂V₁，V₃＝K₃V₂，V₄＝K₄V₃，V₅＝K₅V₄，K₁，K₂，K₃，K₄，K₅为各比例系数；

上述基于风速变异系数的山地风电场模型风机机械功率计算方法，其中：所述步骤5中针对所选择工况S，计算该山地风电场风能利用系数C_p、折减因子α_TM-H、以及第3步骤所得各机组风速V_i和V_j，分别计算该σ_s置信区间对应的2组各机组机械输出功率P_{mech_i}，i＝1,2,…,N，P_{mech_j}，j＝1,2,…,N，具体步骤如下：

步骤5.1、求取该山地风电场风能利用系数Cp，

C_p＝0.5(r-0.022β²-5.6)e^-0.17r

式中：β为桨距角；r满足公式r＝2.237Vw/ω；V_w为风速，m/s；ω为风机转子角速度，rad/s；

步骤5.2、计算山地风电场空气密度折减因子a_TM-H，

式中：ρ_H为海拔高度为H时的空气密度，g/m³；ρ₀为常温、标准大气压力标准状态下空气密度，海平面、15℃条件下空气的密度是1.225g/m³；H为海拔高度，单位m；T₀为绝对温度，取273℃；α为空气温度梯度，取0.0065℃/m；α_H为海拔高度折减因子；

密度与温度、相对湿度、大气压的关系为

其中，t为气温，℃；P为大气压，hP_a；

为相对湿度，％；a_TM为该温度、相对湿度下的空气密度折减因子；

最终得到空气密度与海拔高度、温度、相对湿度、大气压的关系为：

ρ＝α_TMα_Hρ₀＝α_TM-Hρ₀

步骤5.3、求取风力发电机机械功率P_mech

根据标准空气动力模型，计算该风力发电机的机械功率，

式中S为风轮扫过面积(S＝πR²＝3770m²)；R为风轮叶片半径，m；ρ为空气密度,g/m³；V_w为风速,m/s。

总体来看，本发明与现有技术相比，具有以下明显的有益效果：

本发明的建模方法通过利用现有或规划风电场各台风电机组的发电量、平均风速、最大风速、最小风速等历史数据或风能评估数据，计算或估算相应山地风电场的风速变异系数，并以此作为衡量该风电场的风速非一致性指标。从而解决了现有建模技术往往采用一个时段内单一机组、测风塔，甚至整个风电场所有机组的平均风速，不能反映山地区别于平原和海上风电场的风力时空分布弱一致性特征，缺失关于该弱一致性或非一致性特征的量化评价指标，导致难以准确模拟和分析山地风电场接入系统动态特性等技术问题。该指标相较于标准差、离差平方和等，更直观、灵敏、有效地反映了风速的非一致性，解决了现有技术中反映山地风电场风速非一致性指标缺失的问题，有益于风速空间分布的量化模拟，对于弱一致性分布山地风电场接入系统动态分析的准确性提高，有明显的有益效果。引入置信水平，将风速空间分布的点估计问题转变为区间估计问题，有益于预想场景下风电场出力随机性和间歇性的模拟。

总之，本发明通过引入风速非一致性指标(风速变异系数)σ_s，进行山地风电场模型风机机械功率计算方法，相较于现有技术更为准确地模拟山地风电场的时域动态特性。

附图说明

图1为山地风电场风速非一致性指标备选项计算结果的比较图；

图2为本发明的山地风电场模型风机机械功率计算流程图；

图3为实施例中的贵州某山地风电场风机编号及分布示意图；

图4为实施例中的贵州某山地风电场接入系统示意图；

图5为实施例中的山地风电场平均风速略高于切入风速(V_av＝3.5m/s)的情形下各机组风速空间分布(起始时刻)及时间分布图；

图6为实施例中的山地风电场平均风速约为典型低风速(V_av＝7m/s)的情形下各机组风速空间分布(起始时刻)及时间分布图；

图7为实施例中的V_av＝3.5m/s情形下传统等值建模法与本发明所提方法获得的风电场送出线路有功功率、无功功率以及升压站高压侧母线电压时域仿真结果对比图；

图8为实施例中的V_av＝7m/s情形下传统等值建模法与本发明所提方法获得的风电场送出线路有功功率、无功功率以及升压站高压侧母线电压时域仿真结果对比图；

图9为实施例中的不同弱一致性风速指标情形下本发明所提机械功率计算方法所得有功功率仿真结果偏差比例。

具体实施方式

以下结合附图和实施例详细描述本发明的具体实施方式，但本发明不受所述具体实施例所限。

通过以下3个方面，即特性剖析、指标提出、置信水平的分析、讨论与思考，本发明提出了一种基于风速变异系数的山地风电场模型风机机械功率计算方法。

1、特性剖析。

山地风电场风速空间分布非一致性特性及诱因分析如下：(1)高程因素。山地风电场地域起伏，与平原或沿海滩涂区域相比，其风能分布除受到粗糙度、风机尾流、障碍物的影响，还会受到地形和高程变化的影响。其中由于风速随高度切变的原理，而风能与风速呈三次方关系，高程变化是山地风电场风能或者说风速分布变化的最主要影响因素。(2)地形因素。复杂地形变化形成了山脊、山谷、山凹、陡壁、盆地等地貌形式，可能产生迎风面、背风面、喇叭口等情况，也会造成风电场内各处风速与风向变化大、紊流强度不一、风切变、极端风况等不同情况。(3)风机布置因素。有别于平原和海上风电场，主导风向稳定的山地风电场中风机布置间距并非风机发电量差异的决定性因素，而往往与风机所处地及其周围的地形地貌，特别是风机所在山脊的地势与走向更为相关。

2、指标提出。

山地风电场风速非一致性指标的选择与确定。为了选择和确定山地风电场风速不一致性量化指标，以贵州某山区两期风电场的2016年风速实测值为例，进行了数理统计分析与指标比选。采用标准差、变异系数、离差平方和以及最大最小值之差与最大值之比4个备选指标，对该两期风电场逐月风速进行了上述指标值的计算，并给出了相应月平均风速，结果如图1所示。其中图1a、1b、1c以及图1d分别为标准差、变异系数、离差平方和以及最大最小值之差与最大值之比作为山区风电场风速分布非一致性指标时的全年逐月变化情况。对比观察上述4张图可见，该风电场两期工程平均风速有一定差异的情况下，最具全面性、差别性、形象地反映它们风速分布非一致性的指标就是变异系数，例如：对比图1a、1b、1c,3月，5月，6月以及12月计算结果显示两期风电场的风速标准差以及离差平方和差异很小甚至相同，但变异系数确能给出清晰的差异指示。此外，对比图1b和图1d,变异系数和最大最小值之差与最大值之比指标有较好的一致性和差异性指示效果，但8月和9月的计算结果显示后者的数值变化不太稳定，过于激烈，其中8月后者两期工程数值差异接近25％，而前者两期工程差异仅8％左右；9月后者两期工程数值几乎无差异，但前者两期工程差异仍有8％左右。因此，可以认为变异系数作为风速非一致性指标，最能更直观、灵敏、有效的反映山地风电场内风速非一致性水平。

3、置信水平。

在抽样对总体参数作出估计时，由于样本的随机性，其结论总是不确定的。因此，采用一种概率的陈述方法，也就是数理统计中的区间估计法，即估计值与总体参数在一定允许的误差范围以内，其相应的概率有多大，这个相应的概率称作置信水平，可用CL表示。而置信区间是指在某一置信水平下，样本统计值与总体参数值间误差范围。显然，山地风电场风速变化随机性也将导致非一致性指标值变得不确定，因此，很有必要在指标中引入数理统计理论中的置信水平概念，将风速空间分布的点估计问题转变为区间估计问题。其中，对于现有风电场运行方式分析来说，其风速非一致性指标置信水平结合历史数据的质量，可适当取高一些，如0.95，甚至0.98等。而对于规划风电场规划方案分析来说，结合测风塔数据、风能评估数据和风机微观选址情况进行取值，一般来看，一个10～20平方公里的风电场场区会安装1～2个测风塔，但是复杂地形条件各异，测风塔的代表性大大降低，同时计算流体动力学模型也无法保证计算的精确性。因此，其风速非一致性指标置信水平取值相较于现有风电场，可适当调低，如：0.9等。

本发明的技术方案：

参见图2，一种基于风速变异系数的山地风电场模型风机机械功率计算方法，它包括：

步骤1、根据风能历史数据或评估数据，构建山地风电场风速时空特征参数表，涉及工况S、机组数N、机组编号序列W_N、平均风速V_av、风速非一致性指标(风速变异系数)σ_s及其置信区间CL，风速空间分布类型WSD，风速时间分布类型WTD。

步骤2、由步骤1的山地风电场风速时空特征参数表，选择某一工况S，获取风速非一致性指标(风速变异系数)σ_s及其置信水平CL，进而计算变异系数上分位点σ_u和下分位点σ_l。

步骤3、针对所选择工况S，依据步骤2得到的上分位点σ_u、下分位点σ_l，结合山地风电场风速时空特征参数表中该工况的风速空间分布类型WSD，分别计算该σ_s置信区间对应的各机组2组风速V_i，i＝1,2,…,N，V_j，j＝1,2,…,N。

步骤4、针对所选择工况S，由山地风电场风速时空特征参数表，选择风速时间分布类型WTD，并确定参数。

步骤5、针对所选择工况S，计算该山地风电场风能利用系数C_p、折减因子α_TM-H、以及第3步骤所得各机组风速V_i和V_j，分别计算该σ_s置信区间对应的2组各机组机械输出功率P_{mech_i}，i＝1,2,…,N，P_{mech_j}，j＝1,2,…,N。

其中步骤1包括现有山地风电场接入系统分析场景和规划山地风电场接入系统分析场景两类，具体步骤包括：

步骤1.1、对于现有山地风电场接入系统的运行方式分析场景，根据风能历史数据，借助数理统计和聚类方法，获得主要工况S集合，以及所对应的平均风速V_av、风速非一致性指标(风速变异系数)σ_s及其置信区间CL，风速空间分布类型WSD，风速时间分布类型WTD，其中，

(1)给定机组数N下，风电场风速平均值V_av的计算公式为

(2)风电场风速标准差V_δ的计算公式为

(3)提出以变异系数σ_s作为衡量风速非一致性的指标，其计算公式为

其中，V_δ为风电场风速标准差，V_av为风电场风速平均值。

(4)现有山地风电场的变异系数置信区间CL，一般来说可取为，CL＝0.99或0.95。

(5)风速空间分布类型WSD，一般来说涉及，泊松分布、均匀分布、指数分布、正态分布、瑞利分布、韦布尔分布。

(6)风速时间分布类型WTD，一般来说包括，斜坡风、阵风、“墨西哥草帽风”，以及1种自定义分段线性函数风速。

步骤1.2、对于规划山地风电场接入系统的预想事件分析场景，根据风能评估数据，采用类似于步骤1.1的方法，获得主要工况S集合，平均风速V_av、风速非一致性指标(风速变异系数)σ_s及其置信区间CL，风速空间分布类型WSD，风速时间分布类型WTD。但对于规划山地风电场的变异系数置信区间CL，一般来说可取为，CL＝0.95或0.9。

步骤2由步骤1的山地风电场风速时空特征参数表，选择某一工况S，获取风速非一致性指标(风速变异系数)σ_s及其置信水平CL，进而计算变异系数上分位点σ_u和下分位点σ_l。具体步骤包括：

步骤2.1、依据查表法，由Z检验表可确定置信水平下的Z_cl系数，例如：当置信水平CL＝90％时，Z_cl＝1.645；当置信水平CL＝95％时，Z_cl＝1.96；当置信水平CL＝99％时，Z_cl＝2.576。

步骤2.2、计算相应置信区间中变异系数下分位点σ_l：

σ_l＝σ_s-Z_cl×σ_d` (4)

以及相应置信区间中变异系数上分位点σ_u：

σ_u＝σ_s+Z_cl×σ_d (5)

其中，变异系数标准差σ_d一般可取为0.05，或根据实际数据样本计算得到。

步骤3所述针对所选择工况S，依据步骤2得到的上分位点σ_u、下分位点σ_l，结合山地风电场风速时空特征参数表中该工况的风速空间分布类型WSD，分别计算该变异系数置信区间对应的各机组2组风速V_i，i＝1,2,…,N，V_j，j＝1,2,…,N,具体步骤如下：

步骤3.1、根据步骤2选择的工况S，结合步骤1山地风电场风速时空特征参数表中该工况下可选的风速空间分布类型，如：泊松分布、均匀分布、指数分布、正态分布、瑞利分布、韦布尔分布，选择一个风速空间分布类型WSD，选择依据为：现有风电场或类似风电场典型工况或极端工况下的风速空间分布与上述可选风速空间分布类型WSD最为接近，即可选择。

步骤3.2、依据步骤2所得风速非一致性指标(风速变异系数)上分位点σ_u、下分位点σ_l，以及待分析工况的风速平均值V_av，分别可得该工况下变异系数上分位点和下分位点对应的风速标准差V_δu和V_δl。

步骤3.3、将上述参数代入Matlab软件中相应概率分布函数中。以风速空间分布类型WSD为正态分布为例，该函数命令格式为：normrnd(μ1，σ1，m1，n1)，其中，平均值μ1＝V_av，标准差σ1＝V_δu或V_δl，所得结果矩阵行数m1＝N，n1为所得结果矩阵列数，即n1＝1。这样，可得到对应于变异系数上下分位点的2组N个服从该风速空间分布类型WSD(正态分布)的伪随机数据，即该变异系数置信区间对应的2组N台机组的风速V_i,i＝1,2…,N，V_j，j＝1,2,…,N。

类似地，亦可通过Matlab软件中函数poissrnd(lambda，m2，n2)，unifrnd(a，b，m3，n3)、exprnd(MU，m4，n4)分别产生符合风速空间分布类型WSD(泊松分布、均匀分布、指数分布)的伪随机风速数据。其中，lambda，a和b，MU分别为上述分布函数的参数；m2，m3，m4分别为上述函数所得结果矩阵行数；n2，n3，n4分别为上述函数所得结果矩阵列数。

步骤4针对所选择工况S，由山地风电场风速时空特征参数表，选择风速时间分布类型WTD，并确定参数，相关WTD模型及其参数如下：

(1)WTD＝1，即斜坡风模型，其数学模型如下：

式中，斜坡风起始时间t_sr，斜坡风结束时间t_er，风速增加幅值A_r，D_r＝t_er–t_sr。

(2)WTD＝2，即阵风模型，其数学模型如下：

式中，阵风起始时间t_sg，阵风结束时间t_eg，最大风速V_max，A_g＝(V_max–V₀)/2，D_g＝t_eg–t_sg。

(3)WTD＝3，即“墨西哥草帽风”模型，其数学模型如下：

式中，V₀为初始风速、V_max为最大风速、V_min为最小风速，草帽风起始时间t₀，草帽风结束时间t₃，草帽上升沿时间t₁，草帽下降沿时间t₂。

(4)WTD＝4，即一种自定义分段线性函数风速模型，其数学模型如下：

式中，V₀为初始时刻t₀风速，V₁、V₂、V₃、V₄、为中间过程各时间点t₁、t₂、t₃、t₄的风速，V₅为结束时刻t₅的风速，其中，V₁＝K₁V₀，V₂＝K₂V₁，V₃＝K₃V₂，V₄＝K₄V₃，V₅＝K₅V₄，K₁，K₂，K₃，K₄，K₅为各比例系数。

步骤5针对所选择工况S，计算该山地风电场风能利用系数C_p、折减因子α_TM-H、以及第3步骤所得各机组风速V_i和V_j，分别计算该σ_s置信区间对应的2组各机组机械输出功率P_{mech_i}，i＝1,2,…,N，P_{mech_j}，j＝1,2,…,N。具体步骤如下：

步骤5.1、求取该山地风电场风能利用系数Cp，

C_p＝0.5(r-0.022β²-5.6)e^-0.17r (11)

式中：β为桨距角；r满足公式r＝2.237Vw/ω；V_w为风速，m/s；ω为风机转子角速度，rad/s；

步骤5.2、计算山地风电场空气密度折减因子a_TM-H，

式中：ρ_H为海拔高度为H时的空气密度，g/m³；ρ₀为常温、标准大气压力标准状态下空气密度，海平面、15℃条件下空气的密度是1.225g/m³；H为海拔高度，单位m；T₀为绝对温度，取273℃；α为空气温度梯度，取0.0065℃/m；α_H为海拔高度折减因子；

密度与温度、相对湿度、大气压的关系为

其中，t为气温，℃；P为大气压，hP_a；

为相对湿度，％；a_TM为该温度、相对湿度下的空气密度折减因子；

最终得到空气密度与海拔高度、温度、相对湿度、大气压的关系为：

ρ＝α_TMα_Hρ₀＝α_TM-Hρ₀ (14)

步骤5.3、求取风力发电机机械功率P_mech

根据标准空气动力模型，计算该风力发电机的机械功率，

式中S为风轮扫过面积(S＝πR²＝3770m²)；R为风轮叶片半径，m；ρ为空气密度,g/m³；V_w为风速,m/s。

通过山地风电场的风速特性剖析，提出了一种风速非一致性指标(风速变异系数)，同时，继而提出了一种基于风速变异系数的山地风电场模型风机机械功率计算方法，相较于现有技术更为准确地模拟山地风电场的时域动态特性。

具体实施例：

下面以贵州某山地风电场为例进一步对本发明技术方案进行说明，本案例所列举的两种场景，为某山地风电场所有可能出现的场景中的两种，不代表全部。该山地区域属于低风速范畴，风速易受地形的影响，不同位置的风速差异较大。机组切入风速3m/s，平均风速6.5m/s，额定风速9.5m/s，切出风速20m/s。风电场风机分布如图3所示。

场景一平均风速为3.5m/s

1、由步骤1，本场景选取的各参数如表1所示，机组数N为50，机组编号W_N(如图4所示)，平均风速V_av选取稍高于切入风速的3.5m/s，选取变异系数σ_s为15％，置信区间CL为0.95，风速空间分布类型(WSD)选取斜坡风，选取正态分布作为风速时间分布类型(WTD)。

2、由步骤2，选取风速变异系数的标准差σ_d＝0.15，根据式(8)、(9)计算得到风速变异系数上分位点σ_u＝15.294％和下分位点σ_l＝14.706％，再通过式(10)得到各变异系数和平均风速对应的风速标准差，结果如表2所示。

表1山地风电场风速时空特征参数表

表2不同风速下各变异系数对应风速标准差

3、由步骤3，将表2中的参数代入Matlab软件中相应概率分布函数中。函数命令格式为：normrnd(μ1，σ1，m1，n1)。本场景中μ1＝3.5，σ1所取的值如表2所示，m1＝50，n1＝1，运行normrnd函数可得到两组对应的50台机组的风速V_i，i＝1,2…,N，V_j，j＝1,2,…,N，所得两组风机风速初始值结果如图5a所示。

4、由步骤4，选取风速时间分布类型WTD＝1，即斜坡风模型，数学模型如式7所示，V₀为上述3中所得到的某一台风机的风速V_n，起始时间t_sr＝5s，结束时间t_er＝10s，风速增加幅值A_r＝V₀。选取21号风机作为参考，其风速随时间变化曲线如图5b所示。

5、由步骤5，本场景中计算得到该山地风电场风能利用系数C_p、折减因子α_TM-H分别为0.47和0.92。

6、仿真时间为25s，得到V_av＝3.5m/s情形下传统等值建模法与本发明所提方法获得的风电场送出线路有功功率、无功功率以及升压站高压侧母线电压时域仿真结果对比图，如图7所示。

由图7可以清楚地看出，(1)图7a、图7b、图7c中变异系数为15％的置信区间所得到的结果和不考虑变异系数的传统平均风速建模法相比，风电场有出力偏差。图7a中在t＝17s左右有功功率输出偏差最大，接近20MW，输出偏差比例达到了28.6％。图7b中在t＝17s左右无功功率输出偏差最大，接近2Mvar，输出偏差比例达到了23.5％。图7c中在t＝17s左右升压站高压侧母线电压标幺值偏差最大，接近0.001，这主要受风电场SVC动作的影响，所以偏差较小。因此，说明山地风电场的建模需充分考虑风速空间分布非一致性。(2)本场景中变异系数的置信区间上下界随时间的变化呈现出相同的趋势。图7a中上下界最大值相差小于10MW，图7b中上下界最大值相差小于1Mvar，图7c中上下界最大值相差小于0.001。由此可以得到以下结论，置信水平取值合理，置信区间能够较为准确的反映变异系数为15％的情况下各参数的取值范围，对风电场的规划和运行有一定的价值。

场景二平均风速为7m/s

1、由步骤1，本场景选取的各参数如表1所示，机组数N为50，机组编号W_N如图4所示，平均风速V_av选取典型低风速的7m/s，选取变异系数σ_s为15％，置信区间CL为0.95，风速空间分布类型(WSD)选取一种自定义分段线性函数风速，选取正态分布作为风速时间分布类型(WTD)。

2、由步骤2，选取风速变异系数的标准差σ_d＝0.15，根据式(8)、(9)计算得到风速变异系数上分位点σ_u＝15.294％和下分位点σ_l＝14.706％，再通过式(10)得到各变异系数和平均风速对应的风速标准差，结果如表2所示。

3、由步骤3，将表2中的参数代入Matlab软件中相应概率分布函数中。函数命令格式为：normrnd(μ1，σ1，m1，n1)。本场景中μ1＝7，σ1所取的值如表2所示，m1＝50，n1＝1，运行normrnd函数可得到两种组对应的50台机组的风速V_i,i＝1,2…,N，V_j，j＝1,2,…,N，所得两组风机风速初始值结果如图6a所示。

4、由步骤4，选取风速时间分布类型WTD＝4，即一种自定义分段线性函数风速模型，数学模型如式10所示，V₀为上述3中所得到的某一台风机的风速V_n，各参数取值如下：K₁＝1.1，K₂＝1.2，K₃＝0.9，K₄＝0.95，K₅＝0.9，t₀＝3s，t₁＝4.5s，t₂＝5.5s，t₃＝6.5s，t₄＝7.3s，t₅＝8.1s，t₆＝9s。选取21号风机作为参考，其风速随时间变化曲线如图6b所示。

5、由步骤5，本场景中计算得到该山地风电场风能利用系数C_p、折减因子α_TM-H分别为0.47和0.92。

6、仿真时间为25s，得到V_av＝7m/s情形下传统等值建模法与本发明所提方法获得的风电场送出线路有功功率、无功功率以及升压站高压侧母线电压时域仿真结果对比图，如图8所示。

由图8可以清楚的看出，(1)图8a、图8b、图8c中变异系数为15％的置信区间所得到的结果跟不考虑变异系数的传统平均风速建模法相比，风电场有出力偏差。图8a中在t＝6.5s左右有功功率输出偏差最大，接近20MW，输出偏差比例达到了18.75％。图8b中在t＝9.5s左右无功功率输出偏差最大，超过1.5Mvar，输出偏差比例达到了12.5％。图8c中在t＝9.5s左右升压站高压侧母线电压标幺值偏差最大，超过0.001，这主要受风电场SVC动作的影响，所以偏差较小。因此，说明山地风电场的建模需充分考虑风速空间分布非一致性。(2)本场景中变异系数的置信区间上下界随时间的变化呈现出相同的趋势。图8a中上下界最大值相差小于5MW，图8b中上下界最大值相差小于0.5Mvar，图8c中上下界最大值相差小于0.001。由此可以得到以下结论，置信水平取值合理，置信区间能够较为准确的反映变异系数为15％的情况下各参数的取值范围，对风电场的规划和运行有一定的价值。

此外，进一步地针对不同平均风速下，对风电场各变异系数对应的有功功率输出进行时域仿真，结果如图9所示。

由图9可以清楚的看出，(1)随着变异系数的增大，考虑与不考虑变异系数的风电场出力偏差比例逐渐增大，特别是超过15％以后，对于额定或高平均风速工况，出力偏差比例逼近变异系数百分比，而对接近切入风速的低平均风速工况，受低速停机的影响，出力偏差比例已显著超过变异系数百分比，理论极端情况可超过100％。(2)风电场平均风速越低，其风速空间分布非一致性指标，变异系数变化导致，考虑与不考虑变异系数的风电场出力偏差比例越剧烈，或者说偏差值曲线的斜率越大。(3)总的来看，初步可以得出以下结论，含山地风电场的电力系统暂态或动态分析中，对于平均风速接近风机切入风速3m/s或2.5m/s的山地风电场，应考虑空间分布非一致性指标，即变异系数的影响；对于平均风速超过风机额定风速的山地风电场，若该工况下变异系数超过15％，应在风电场建模与风机风速模拟中考虑空间分布非一致性指标，即变异系数的影响。

综上所述，本发明提出了一种基于风速变异系数的山地风电场模型风机机械功率计算方法，上述实施案例计算结果表明：

1)本发明针对现有建模技术往往采用一个时段内单一机组、测风塔甚至整个风电场所有机组的平均风速，不能反映山地区别于平原和海上风电场的风力时空分布弱一致性特征，缺失关于该弱一致性或非一致性特征的量化评价指标，导致难以准确模拟和分析山地风电场接入系统动态特性等技术问题，该建模方法通过利用现有或规划风电场各台风电机组的发电量、平均风速、最大风速、最小风速等历史数据或风能评估数据，计算或估算相应山地风电场的风速变异系数，并以此作为衡量该风电场的风速非一致性指标。该指标相较于标准差、离差平方和等，更直观、灵敏、有效地反映了风速的非一致性，解决了现有技术中反映山地风电场风速非一致性指标缺失的问题，有益于风速空间分布的量化模拟，对于弱一致性分布山地风电场接入系统动态分析的准确性提高，有明显的有益效果。

2)此外，引入置信水平，将风速空间分布的点估计问题转变为区间估计问题，有益于预想场景下风电场出力随机性和间歇性的模拟。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，任何未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (8)

1.一种基于风速变异系数的山地风电场模型风机机械功率计算方法，包括以下步骤：

步骤1、根据风能数据，构建山地风电场风速时空特征参数表，参数表的参数包括工况S、机组数N、机组编号序列W_N、平均风速V_av、风速非一致性指标的风速变异系数σ_s及其置信区间CL，风速空间分布类型WSD，风速时间分布类型WTD，其中：

(1)给定机组数N下，风电场风速平均值V_av的计算公式为

(2)风电场风速标准差V_δ的计算公式为

(3)以变异系数σ_s作为衡量风速非一致性的指标，其计算公式为

其中，V_δ为风电场风速标准差，V_av为风电场风速平均值；

步骤2、由步骤1的山地风电场风速时空特征参数表，选择某一工况S，获取风速非一致性指标的风速变异系数σ_s及其置信水平CL，进而计算变异系数上分位点σ_u和下分位点σ_l；

步骤3、针对所选择工况S，依据步骤2得到的上分位点σ_u、下分位点σ_l，结合山地风电场风速时空特征参数表中该工况的风速空间分布类型WSD，分别计算该σ_s置信区间对应的各机组风速；

步骤4、针对所选择工况S，由山地风电场风速时空特征参数表，选择风速时间分布类型WTD，并确定相关参数；

步骤5、针对所选择工况S，计算该山地风电场风能利用系数C_p、折减因子α_TM-H、以及第3步骤所得各机组风速，分别计算该σ_s置信区间对应的机组机械输出功率。

2.如权利要求1所述的基于风速变异系数的山地风电场模型风机机械功率计算方法，其特征在于：所述步骤1的山地风电场包括现有山地风电场接入系统分析场景和规划山地风电场接入系统分析场景。

3.如权利要求2所述的基于风速变异系数的山地风电场模型风机机械功率计算方法，其特征在于：所述现有山地风电场接入系统的运行方式分析场景，根据风能历史数据，借助数理统计和聚类方法，获得主要工况S集合，以及所对应的平均风速V_av、风速非一致性指标的风速变异系数σ_s及其置信区间CL，风速空间分布类型WSD，风速时间分布类型WTD，其中：

(1)给定机组数N下，风电场风速平均值V_av的计算公式为

(2)风电场风速标准差V_δ的计算公式为

(3)以变异系数σ_s作为衡量风速非一致性的指标，其计算公式为

其中，V_δ为风电场风速标准差，V_av为风电场风速平均值；

(4)现有山地风电场的变异系数置信区间CL，可取为，CL＝0.99或0.95；

(5)风速空间分布类型WSD，涉及泊松分布、均匀分布、指数分布、正态分布、瑞利分布、韦布尔分布；

(6)风速时间分布类型WTD，包括斜坡风、阵风、“墨西哥草帽风”，以及自定义分段线性函数风速。

4.如权利要求2所述的基于风速变异系数的山地风电场模型风机机械功率计算方法，其特征在于：所述规划山地风电场接入系统的预想事件分析场景，根据风能评估数据，借助数理统计和聚类方法，获得主要工况S集合，以及所对应的平均风速V_av、风速非一致性指标的风速变异系数σ_s及其置信区间CL，风速空间分布类型WSD，风速时间分布类型WTD，其中：

(1)给定机组数N下，风电场风速平均值V_av的计算公式为

(2)风电场风速标准差V_δ的计算公式为

(3)以变异系数σ_s作为衡量风速非一致性的指标，其计算公式为

其中，V_δ为风电场风速标准差，V_av为风电场风速平均值；

(4)规划山地风电场的变异系数置信区间CL，可取为，CL＝0.95或0.9；

(5)风速空间分布类型WSD，涉及泊松分布、均匀分布、指数分布、正态分布、瑞利分布、韦布尔分布；

(6)风速时间分布类型WTD，包括斜坡风、阵风、“墨西哥草帽风”，以及自定义分段线性函数风速。

5.如权利要求1、2、3或4任一项所述的基于风速变异系数的山地风电场模型风机机械功率计算方法，其特征在于：所述步骤2中由步骤1的山地风电场风速时空特征参数表，选择某一工况S，获取风速非一致性指标的风速变异系数σ_s及其置信水平CL，进而计算变异系数上分位点σ_u和下分位点σ_l，其具体步骤包括：

步骤2.1、依据查表法，由Z检验表可确定置信水平下的Z_cl系数；

步骤2.2、计算相应置信区间中变异系数下分位点σ_l：

σ_l＝σ_s-Z_cl×σ_d

以及相应置信区间中变异系数上分位点σ_u：

σ_u＝σ_s+Z_cl×σ_d

其中，变异系数标准差σ_d可取为0.05，或根据实际数据样本计算得到。

6.如权利要求1、2、3或4任一项所述的基于风速变异系数的山地风电场模型风机机械功率计算方法，其特征在于：所述步骤3中针对所选择工况S，依据步骤2得到的上分位点σ_u、下分位点σ_l，结合山地风电场风速时空特征参数表中该工况的风速空间分布类型WSD，分别计算该变异系数置信区间对应的各机组2组风速V_i，i＝1,2,…,N，V_j，j＝1,2,…,N，具体步骤如下：

步骤3.1、根据步骤2选择的工况S，结合步骤1山地风电场风速时空特征参数表中该工况下可选的风速空间分布类型，选择一个风速空间分布类型WSD，选择依据为：现有风电场或类似风电场典型工况或极端工况下的风速空间分布与上述可选风速空间分布类型WSD最为接近，即可选择；

步骤3.2、依据步骤2所得风速非一致性指标的风速变异系数的上分位点σ_u、下分位点σ_l，以及待分析工况的风速平均值V_av，分别可得该工况下变异系数上分位点和下分位点对应的风速标准差V_δu和V_δl

步骤3.3、将上述参数代入相应概率分布函数中，得到对应于变异系数上下分位点的2组N个服从该风速空间分布类型WSD的正态分布的伪随机数据，即该变异系数置信区间对应的2组N台机组的风速V_i,i＝1,2…,N，V_j，j＝1,2,…,N。

7.如权利要求1、2、3或4任一项所述的基于风速变异系数的山地风电场模型风机机械功率计算方法，其特征在于：所述步骤4中针对所选择工况S，由山地风电场风速时空特征参数表，选择风速时间分布类型WTD，并确定参数，相关WTD模型及其参数如下：

(1)WTD＝1，即斜坡风模型，其数学模型如下：

式中，V₀为初始风速，t_sr为斜坡风起始时间，t_er为斜坡风结束时间，A_r为风速增加幅值，D_r＝t_er–t_sr；

(2)WTD＝2，即阵风模型，其数学模型如下：

式中，V₀为初始风速，t_sg为阵风起始时间，t_eg为阵风结束时间，V_max为最大风速，A_g＝(V_max–V₀)/2，D_g＝t_eg–t_sg；

(3)WTD＝3，即“墨西哥草帽风”模型，其数学模型如下：

式中，V₀为初始风速、V_max为最大风速、V_min为最小风速，草帽风起始时间为t₀，草帽风结束时间为t₃，草帽上升沿时间为t₁，草帽下降沿时间为t₂；

(4)WTD＝4，即一种自定义分段线性函数风速模型，其数学模型如下：

式中，V₀为初始时刻t₀风速，V₁、V₂、V₃和V₄分别为中间过程各时间点t₁、t₂、t₃、t₄的风速，V₅为结束时刻t₅的风速，其中，V₁＝K₁V₀，V₂＝K₂V₁，V₃＝K₃V₂，V₄＝K₄V₃，V₅＝K₅V₄，K₁，K₂，K₃，K₄，K₅为各比例系数。

8.如权利要求1、2、3或4任一项所述的基于风速变异系数的山地风电场模型风机机械功率计算方法，其特征在于：所述步骤5中针对所选择工况S，计算该山地风电场风能利用系数C_p、折减因子α_TM-H、以及第3步骤所得各机组风速V_i和V_j，分别计算该σ_s置信区间对应的2组各机组机械输出功率P_{mech_i}，i＝1,2,…,N，P_{mech_j}，j＝1,2,…,N，具体步骤如下：

步骤5.1、求取该山地风电场风能利用系数Cp，

C_p＝0.5(r-0.022β²-5.6)e^-0.17r

式中：β为桨距角；r满足公式r＝2.237Vw/ω；V_w为风速，m/s；ω为风机转子角速度，rad/s；

步骤5.2、计算山地风电场空气密度折减因子a_TM-H，

式中：ρ_H为海拔高度为H时的空气密度，g/m³；ρ₀为常温、标准大气压力标准状态下空气密度，海平面、15℃条件下空气的密度是1.225g/m³；H为海拔高度，单位m；T₀为绝对温度，取273℃；α为空气温度梯度，取0.0065℃/m；α_H为海拔高度折减因子；

密度与温度、相对湿度、大气压的关系为

其中，t为气温，℃；P为大气压，hP_a；

为相对湿度，％；a_TM为该温度、相对湿度下的空气密度折减因子；

最终得到空气密度与海拔高度、温度、相对湿度、大气压的关系为：

ρ＝α_TMα_Hρ₀＝α_TM-Hρ₀

步骤5.3、求取风力发电机机械功率P_mech，

根据标准空气动力模型，计算该风力发电机的机械功率，

式中S为风轮扫过面积，S＝πR²＝3770m²；R为风轮叶片半径，m；ρ为空气密度,g/m³；V_w为风速,m/s。

Images