-
风资源评估是风电开发的重要环节[1]，风功率密度又是风能资源评估的重要参数之一[2]，是决定风电场是否具有开发潜能的重要评判标准。在年平均风速和空气密度相同的情况下，不同的风频分布将计算得出不同的风功率密度，因而准确分析风资源的风频分布具有重要意义。
目前，对风频分布的研究多采用威布尔分布，它能比较好地拟合风资源的实际情况。刘鹏等[3]采用了矩估计法模拟风频威布尔分布，此方法简单方便，缺点是不能完全利用风数据样本信息，模拟精度不高。Justus和Lysen提出的经验计算方法[4]，被2002版《全国风能评价技术规定》采纳，因简单、实用、耗时短，在风电场建设初期做初步评估时被广泛采用。Akdag 和Dinler提出了用能量因子法来模拟威布尔分布参数[5]，此方法目前被广泛应用于风资源流体力学评估模型软件中，如Meteodyn Universe。Stevens MJM 将最大似然法运用于风资源风频分布的模拟计算，并以此估算风电场发电量[6]。最小二乘法由法国勒让德创立，是处理实测值的纯粹代数方法，此方法被广泛应用于各行业，这是一种寻找最佳解，揭示最接近真实情形的状态[7]。
估算风功率密度的方法较多，准确地计算风功率密度有赖于风频威布尔分布拟合的准确性[8]。但是，目前还缺少基于威布尔分布拟合准确性计算风功率密度方面的研究。因此，文章对目前国际常用的5种风频威布尔拟合方法进行对比研究，采用决定系数对威布尔拟合结果准确度进行评价，并在此基础上利用威布尔分布计算风功率密度，并与实测风数据计算所得风功率密度相比较绝对误差和相对误差，从而确定哪些方法更可靠，更准确。文章选取我国不同经纬度和地形地貌的4座测风塔实测数据进行威布尔分布拟合和决定系数分析，并计算相应的风功率密度，用案例进行对比验证得出较优方法，可为工程师在实际的风资源评估工作中提供参考和指导。
-
威布尔分布是由瑞典科学家威布尔从材料强度的统计理论中推导出来的，是日常生活中一种常用的分布，在风电领域的风频分布中也被广泛采纳应用。
威布尔的形状参数$k$值反映了风速分布的范围宽窄，一般情况下，$k$值越小，风速分布越宽。当$0 < k < 1$时，风频随着风速的增大而减小，在现实情况中出现较少；当$k = 1$时，风频为指数分布；当$k = 2$时，为瑞利分布；当$k = 3.5$时，接近正态分布。威布尔分布累计概率密度函数如式（1）所示：
$$ {F_{(v)}} = 1 - {\rm{EXP}}\left[ { - {\left( {\dfrac{v}{c}} \right)^k}} \right] $$ (1) 式中：
${F_{(v)}}$——威布尔风频累计概率密度（%）；
$v$ ——风速（m/s）；
$c$ ——威布尔尺度参数（m/s）；
$k$ ——为威布尔形状参数。
对威布尔累计概率密度函数（1）求导数，得到威布尔风频分布概率密度函数，如式（2）所示：
$$ {f_{(v)}} = \dfrac{k}{c}{\left( {\dfrac{v}{c}} \right)^{(k - 1)}}{\rm{EXP}}\left[ { - {\left( {\frac{v}{c}} \right)^k}} \right] $$ (2) 式中：
${f_{(v)}}$——威布尔风频分布概率密度（%）。
风功率密度是指空气流垂直通过单位截面积的风能，是表征1个地方风能资源多少的1个指标，在风资源评估中一般指年均风功率密度，某地的风功率密度越大，说明该地区风能资源越好[9]。风功率密度计算方法有两种，一种是基于时间序列实测风数据[10]，风功率密度与风速的3次方成正比，其计算公式如式（3）所示：
$$ \overline P = \dfrac{1}{{2n}}\rho \displaystyle \sum\limits_{i = 1}^n {{v_i}^3} = \dfrac{1}{2}\rho \overline {{v_i}^3} $$ (3) 式中：
$\overline P $——年均风功率密度（W/m2）；
$\rho $——空气密度（kg/m3）；
${v_i}$——风速时间序列值（m/s)；
$n$——时间序列。
另一种是选择合适风速分布的数学模型来模拟风速变化，常见的有双参数威布尔分布、瑞利分布、伽马分布、对数正态分布，大多学者认为双参数威布尔分布可以较好的模拟风速变化[11-13]。基于威布尔分布的风功率密度计算公式如式（4）所示：
$$ \overline P = \dfrac{1}{2}\rho \int_0^n {{v^3}{f_{(v)}}} {\rm{d}}v $$ (4) -
该经验方法由Justus提出[5]，我国颁布的《全国风能评价技术规定》中也采用此种方法，此方法计算简单。其$k$值和$c$值计算公式如式（5）和（6）[14]所示：
$$ k = {\left( {\dfrac{\sigma }{{\overline {{v_i}} }}} \right)^{ - 1.086}} $$ (5) $$ c = \dfrac{{\overline {{v_i}} }}{{\Gamma (1 + 1/k)}} $$ (6) 式中：
$\sigma $ ——风速标准差（m/s）；
$ \Gamma （1+1/k） $——伽马函数；
$\overline {{v_i}} $ ——风速序列平均值（m/s)。
-
此方法由Lysen提出[15]，在Justus提出的方法上有所改进，其$k$值计算由式（5）给出，$c$值计算公式如式（7）所示：
$$ c={\overline{{v}_i}}(0.568+0.433 / k)^{-\tfrac{1}{k}} $$ (7) -
基于能量守恒规律，使其威布尔拟合计算得到的风功率密度与实测数据计算所得的风功率密度差值尽可能小的情况下求$k$和$c$值。此方法首先求出能量因子${E_{{\rm{pf}}}}$，其计算公式如式（8）所示，再计算$k$值[16]，$k$值计算公式如式（9）所示，$c$值计算公式同式（7）：
$$ {E_{{\rm{pf}}}} = \dfrac{{\dfrac{1}{n}\displaystyle \sum\limits_{i = 1}^n {v_i^3} }}{{{{\left( {\dfrac{1}{n}\displaystyle \sum \limits_{i = 1}^n {v_i}} \right)}^3}}} = \dfrac{{\overline {({v_i}^3)} }}{{{{(\overline {{v_i}} )}^3}}} = \dfrac{{\Gamma (1 + 3/k)}}{{{\Gamma ^3}(1 + 1/k)}} $$ (8) $$ k = 1 + \dfrac{{3.69}}{{{{({E_{{\rm{pf}}}})}^2}}} $$ (9) 式中：
${E_{{\rm{pf}}}}$——风功率密度能量因子（W/m2）。
-
1821年由德国数学家高斯提出最大似然法，此方法首次用于计算发电量可以追溯到Stevens MJM [6]，用这种模型来估算风频分布参数$k$和$c$值，其中蕴含的规则是：在求得的$k$和$c$这个特定的威布尔分布函数中，能得到现在的实测风数据序列样本概率最大。
假设某一风速时间序列样本为：
$$ V = ({v_1},{v_2},{v_3}...{v_n}) $$ 威布尔概率密度似然函数$f$在某一风速${v_i}$下的公式如式（10）所示：
$$ {f_{({v_i}\left| {k,c} \right.)}} = \dfrac{k}{c}{\left( {\dfrac{{{v_i}}}{c}} \right)^{(k - 1)}}{\rm{EXP}}\left[ { - {\left( {\dfrac{{{v_i}}}{c}} \right)^k}} \right] $$ (10) 则风速时间序列样本$V$产生的联合密度函数为：
$$ \begin{split} &f({v}_{i}|k,c)=f({v}_{1},{v}_{2},{v}_{3},\cdots,{v}_{n}|k,c) =\\&f({v}_{1}|k,c)\times f({v}_{2}|k,c) \times f({v}_{3}|k,c)\times \cdots\times\\& f({v}_{n}|k,c) ={\displaystyle \prod _{i=1}^{n}f({v}_{i}}|k,c）\end{split} $$ (11) 为了计算方便，等式两边取对数：
$$ \begin{split} &\ln [f({v_i}\left| {k,c)] = \ln } \right.[\prod\limits_{i = 1}^n {f({v_i}\left| {k,c)} \right.} ] =\\& \ln [f({v_1}\left| {k,c)]} \right. + \ln [f({v_2}\left| {k,c} \right.)] +\\& \ln [f({v_3}\left| {k,c} \right.)] + ...... + \ln [f({v_n}\left| {k,c} \right.)] \\ \end{split} $$ (12) 上述方程分别对$k$，$c$求导，并令其分别为0，见式（13）和式（15），即可求得威布尔分布的参数$k$和$c$[17]，具体见式（14）和式（16）：
$$ \dfrac{{\rm{d}}}{{\rm{d}}k}\mathrm{ln}[f(v|k,c)]=0 $$ (13) $$ k = {\left[ {\dfrac{{\displaystyle \sum\nolimits_{i = 1}^n {v_i^k\ln ({v_i})} }}{{\sum\limits_{i = 1}^n {v_i^k} }} - \dfrac{{\displaystyle \sum\limits_{i = 1}^n {\ln ({v_i})} }}{n}} \right]^{ - 1}} $$ (14) $$ \dfrac{{\rm{d}}}{{{\rm{d}}c}}\ln [f(v\left| {k,c} \right.)] = 0 $$ (15) $$ c = {\left[ {\dfrac{1}{n}\displaystyle \sum\limits_{i = 1}^n {{v_i}^k} } \right]^{1/k}} $$ (16) -
风速时间序列某一风速的累计概率密度分布函数公式如式（17）所示：
$$ F({v_i}) = 1 - {\rm{EXP}}\left[ { - {\left( {\dfrac{{{v_i}}}{c}} \right)^k}} \right] \Rightarrow \frac{1}{1-F\left( {{v}_{i}} \right)}={\rm{EXP}}{\left( {\dfrac{{v}_{i}}{c}} \right)}^{k} $$ (17) 为计算方便，等式两边取对数如式（18）所示：
$$ \ln \left[ {\frac{1}{{1 - F({v_i})}}} \right] = {\left( {\frac{{{v_i}}}{c}} \right)^k} \ln [ - \ln (1 - F({v_i})] = k\ln ({v_i}) - k\ln (c) \\ $$ (18) 累计概率密度函数经过上述的变化推导就变成常见的线性$y = mx + b$形式。最小二乘法的理论基石是实测值与模拟求得的结果值之间的差的平方总和最小。在用于威布尔拟合中，其公式如式（19）所示：
$$ {\rm{S}}{{\rm{S}}_{{\rm{err}}}} = {\displaystyle \sum\limits_{i = 1}^n {\left( {(k{\text{ln}}({v_i}) - k\ln (c)) - {\text{ln}}\left[ {\dfrac{1}{{1 - F({v_i})}}} \right]} \right)} ^2} $$ (19) 式中：
${\rm{S}}{{\rm{S}}_{{\rm{err}}}}$——为实测值与模拟结果值差的平方总和。
式（19）对$k$和$c$分别求导，并分别令其为0,得出$k$值和$c$值[18]，如式（20）和式（21）所示：
$$ k = \dfrac{{\displaystyle \sum\limits_{i = 1}^n {[(\ln ({v_i}) - \overline {\ln ({v_i})} ]\left[ {\ln \left( {\dfrac{1}{{1 - F({v_i})}}} \right) - \overline {\ln \left( {\dfrac{1}{{1 - F({v_i})}}} \right)} } \right]} }}{{\displaystyle \sum\limits_{i = 1}^n {[{{(\ln ({v_i}) - \overline {\ln ({v_i})} ]}^2}} }} $$ (20) $$ \ln c = \dfrac{{\left( {\overline {\ln \left( {\dfrac{1}{{1 - F({v_i})}}} \right)}} \right) - k\overline {\ln ({v_i})} }}{{ - k}} $$ (21) 上述最大似然法和最小二乘法所涉及的公式属于理论推导，运算较复杂，在涉及数据繁多的测风数据处理时，常常借助Excel中的规划求解器Slover来解决，Slover规划求解器是一种多次运算组合模拟参数使得目标函数取得最大值或者最小值的解决办法[19],此方法方便快捷，批量处理测风数据能提高计算的效率和准确性。具体运用于此文章时将多次组合模拟$k$和$c$值，使得联合密度函数$ {\displaystyle \prod _{i=1}^{n}f({v}_{i}}|k,c） $取得最大值，最小二乘法函数即式（19）中的${{\rm{SS}}_{{\rm{rr}}}}$取得最小值，此时相应的$k$和$c$即为最优解。
-
在实际运用中，针对风场实测数据的风频分布情况一般采用上述5种方法进行拟合，文章采用决定系数R2来衡量其拟合结果的精确度[20]，其计算公式如式（22）所示：
$$ {R^2} = 1 - \dfrac{{\displaystyle \sum\nolimits_{i = 1}^n {{{[{f_{\rm{w}}}({v_i}) - {f_{\rm{m}}}({v_i})]}^2}} }}{{\displaystyle \sum\nolimits_{i = 1}^n {{{[{f_{\rm{w}}}({v_i}) - \overline {{f_{\rm{m}}}({v_i})} ]}^2}} }} $$ (22) $$ \overline {{f_{\rm{m}}}({v_i})} = \dfrac{1}{{n + 1}}\displaystyle \sum\limits_{i = 1}^n {{f_{\rm{m}}}} ({v_i}) $$ (23) 式中：
${f_{\rm{w}}}({v_i})$——威布尔分布参数计算风速为${v_i}$时的风频（%）；
${f_{\rm{m}}}({v_i})$——实测风速序列为${v_i}$时的风频（%）；
$\overline {{f_{\rm{m}}}({v_i})}$——实测风速序列风频的平均值（%）；
$ 0.9 < {R^2} \leqslant 1 $，拟合精度好；
$ 0.8 < {R^2} \leqslant 0.9 $，拟合精度较好；
$ 0.7 < {R^2} \leqslant 0.8 $，拟合精度一般；
$ {R^2} \leqslant 0.7 $，拟合精度差。
-
文章选取我国不同经纬度，不同地形地貌及风资源条件各不相同的4座测风塔，1个完整年统计值为10 min的风速样本52560个数据进行分析[21]。测风塔1位于内蒙古茂明的平坦地形，植被稀少地貌非常简单；测风塔2位于云南昭通东部海拔2100 m左右的多山地带，植被较多，地形复杂；测风塔3位于新疆小草湖的半沙漠干旱地带；测风塔4位于广东低纬度甲东地区的沿海滩涂附近。其测风塔所在地理区域位置如图1所示。
Figure 1. The geographic position diagram of 4 masts
测风塔的经纬度，所在地空气密度，年平均风速和实测风功率密度等基本参数如表1所示。
测风塔编号 经度（E）/（°） 纬度（N）/（°） 空气密度/(kg·m−3) 年均风速/（m·s−1) 风速标准差/（m·s−1) 最大值/（m·s−1) 风功率密度/（W·m−2） 1 109.77 41.58 1.12 7.4 3.7 25.3 467.50 2 103.61 27.35 0.95 7.8 3.1 27.1 426.50 3 88.69 43.07 1.10 9.1 4.6 29.5 611.90 4 116.10 22.81 1.21 7.0 4.1 24.1 521.70 Table 1. Wind characteristic parameters
将测风数据以0.5 m/s作为一个风速区间范围进行频率统计，采用上文介绍的5种威布尔拟合方法分别对收集到的4座测风塔数据进行模拟计算得$k$和$c$值，用威布尔分布函数可以求出各风速段相对应的频率，在此基础上利用式（4）计算威布尔分布风功率密度和式（22）计算决定系数，各测风塔结果分别如表2~表5所示。
采用方法 经验法EMJ 经验法EML 能量因子法EPF 最大似然法MLE 最小二乘法LLSA 相应参数 威布尔参数K 2.12 2.12 2.10 2.08 2.03 威布尔参数C 8.39 8.40 8.40 8.40 8.40 风功率密度P/（W·m−2） 454.20 454.90 457.80 459.50 481.60 决定系数R2 0.96 0.96 0.97 0.97 0.96 Table 2. Simulation results of Weibull distribution and power density of mast 1
采用方法 经验法EMJ 经验法EML 能量因子法EPF 最大似然法MLE 最小二乘法LLSA 相应参数 威布尔参数K 1.80 1.80 2.20 2.10 1.78 威布尔参数C 8.80 8.81 8.90 8.80 8.70 风功率密度P/（W·m−2） 473.46 474.94 405.82 408.64 463.94 决定系数R2 0.85 0.85 0.92 0.91 0.77 Table 3. Simulation results of Weibull distribution and power density of mast 2
采用方法 经验法
EMJ经验法
EML能量因子法
EPF最大似然法
MLE最小二乘法
LLSA相应参数 威布尔参数K 2.10 2.10 2.30 2.10 1.77 威布尔参数C 10.27 10.28 10.00 9.70 9.50 风功率密度
P/（W·m−2）737.84 739.83 640.71 628.43 680.36 决定系数R2 0.69 0.69 0.83 0.85 0.75 Table 4. Simulation results of Weibull distribution and power density of mast 3
采用方法 经验法EMJ 经验法EML 能量因子法EPF 最大似然法MLE 最小二乘法LLSA 相应参数 威布尔参数K 1.65 1.65 1.54 1.56 1.45 威布尔参数C 7.80 7.81 7.80 7.90 8.02 风功率密度P/（W·m−2） 475.58 477.26 515.81 525.36 588.58 决定系数R2 0.92 0.92 0.96 0.98 0.95 Table 5. Simulation results of Weibull distribution and power density of mast 4
从表2~表5模拟结果可知，4座测风塔用经验法EMJ和EML计算的参数值K分别相等，C值均相差0.01，非常接近。
测风塔1和测风塔4用5种威布尔方法模拟风频结果均好，决定系数均高于0.90，能量因子法EPF和最大似然法MLE的决定系数相对更高，均高于0.96。相应威布尔计算所得的风功率密度与实测数据计算所得风功率密度绝对误差和相对误差也小，其中测风塔1使用能量因子法EPF和最大似然法MLE与实测风数据计算所得风功率467.50 W/m2绝对误差分别为9.70 W/m2和8.00 W/m2，相对误差值分别为2.1%和1.7%。测风塔4用5种威布尔方法计算所得风功率密度和实测风功率密度绝对误差分别为46.12 W/m2，44.44 W/m2，5.89 W/m2，3.66 W/m2和66.88 W/m2。
测风塔2和测风塔3威布尔模拟的决定系数横向比较也有相似的趋势，能量因子法EPF和最大似然法MLE的决定系数高于其他3种方法，其计算所得的风功率密度与实测数据所得风功率密度更加接近，绝对值误差和相对误差均小于其他3种方法。尽管测风塔3模拟决定系数整体偏低，其分布范围为0.69~0.85之间，但能量因子法EPF和最大似然法MLE结果拟合精度较好分别为0.83和0.85，其相对较高的趋势仍未变化。
测风塔3用能量因子法EPF和最大似然法MLE计算所得的风功率密度与实测数据计算所得风功率密度绝对误差为28.81 W/m2和16.53 W/m2，相对误差分别为4.7%和2.7%；经验法EMJ和EML绝对误差较大，分别为125.94 W/m2和127.93 W/m2，最小二乘法LLSA为68.46 W/m2，相对误差分别为20.6%，17.3%和11.2%。
测风塔3用经验法EMJ和EML及最小二乘法LLSA威布尔模拟结果一般，决定系数分别为0.69，0.69和0.75，在此基础上利用风频分布计算所得的风功率密度误差更大，这也更好地验证了风频分布拟合准确与否对风功率密度计算准确度影响至关重要，风速与风功率密度是3次方关系。
4座测风塔实测风频分布图和威布尔分布模拟分布情况如图2~图5所示。
Figure 2. Wind speed frequency of the measurement data and Weibull distribution simulation of mast 1
Figure 3. Wind speed frequency of the measurement data and Weibull distribution simulation of mast 2
Figure 4. Wind speed frequency of the measurement data and Weibull distribution simulation of mast 3
Figure 5. Wind speed frequency of the measurement data and Weibull distribution simulation of mast 4
-
文章针对风频分布拟合方法进行比较研究和风功率密度计算结果验证，引用决定系数对风频分布拟合度进行优劣评价，确定威布尔拟合常用的5种方法准确性高低。为了使结论具有普适性和代表性，文章采用我国不同经纬度，地形地貌各异和风资源禀赋程度各不相同的4座实测测风塔1年统计数据进行验证。
最终得出结论：威布尔拟合常用的5种计算方法中，能量因子法EPF和最大似然法MLE能更好地拟合实测风数据的风频，采用威布尔计算所得的风功率密度和实测风数据计算的风功率密度更为接近，误差较小，决定系数值较大，其风频的拟合度优于经验法和最小二乘法；经验法和最小二乘法拟合度相对差一些,相应计算的风功率密度误差也大一些。尽管某些测风塔用5种方法拟合所得的决定系数相对都较小，拟合度一般，但横向比较能量因子法EPF和最大似然法MLE的拟合度仍较经验法和最小二乘法有优势，相应计算所得的风功率密度误差仍小于其他3种方法。
这一结果与目前我国较流行的流体力学风资源模拟软件Meteodyn Universe的运算规则相吻合[22]，这款软件在做各高层绘图计算时用的是能量因子法EPF，在各风机机位点参数计算$k$值，$c$值时运用的是最大似然法MLE。
因此，在实际项目风资源评估中，建议采用能量因子法EPF和最大似然法MLE来模拟风数据的风频分布，计算得出更加接近实际的风功率密度，可为更好地了解掌握风场的风资源提供可靠依据。
Comparison of Wind Power Density Calculation Methods Based on Weibull Distribution
doi: 10.16516/j.ceec.2024.1.04
- Received Date: 2023-04-04
- Rev Recd Date: 2023-06-09
- Available Online: 2023-09-07
- Publish Date: 2024-01-10
-
Key words:
- Weibull distribution /
- wind power density /
- energy pattern factor (EPF) method /
- maximum likelihood estimation (MLE) method /
- coefficient of determination
Abstract:
| Citation: | LI Hua. Comparison of wind power density calculation methods based on Weibull distribution [J]. Southern energy construction, 2024, 11(1): 33-41 doi: 10.16516/j.ceec.2024.1.04
|
DownLoad:
