基于测风塔数据的地表粗糙度计算方法研究

项 雯，彭秀芳，施 晨，冯 浩

(中国能源建设集团江苏省电力设计院有限公司，江苏 南京 211102)

0 引言

风能资源评估作为建设风电场可行性的基础，是发展风力发电的重要前提条件。采用数值模拟手段可以获取开发利用高度处风能资源量的区域分布状况，而地表粗糙度对近地面风速的数值模拟结果影响较大，因此，在风能资源评估中合理确定地表粗糙度值是一个非常重要的问题。

理论上，地表粗糙度z0是平均风速随高度减小到零时的高度[1]。赵晓松[2]等根据牛顿迭代法计算零平面位移，进而得到地表粗糙度；Bruin[3]等根据质量守恒定律计算地表粗糙度。Martano[4]根据Monin Obukhov相似理论，利用单层超声风温资料，将计算零平面位移和空气动力学粗糙度问题简化为一个可由最小二乘法求单变量的过程，单独求出z0。这些方法基本上需要利用多层风温资料或者一层三维超声数据，而在实际工程项目中，一般收集到的测风塔数据通常只有一层温度数据，基本没有三维超声数据，这就给运用上述类似方法估算地表粗糙度带来极大阻碍。

基于此，本文对平坦地形地表粗糙度计算方法展开研究，提出一种仅以测风塔实测数据为基础计算地表粗糙度值的方法，并通过实例进行验证。本文方法所需资料简单易获取，为今后设计人员估算地表粗糙度值、准确评估风能资源提供帮助。

1 粗糙度计算方法

1.1 风廓线法

在近地面层(即产生地点为贴地层至其上方50～100 m高度处[5])风速随高度是否遵从对数规律分布主要由当时的大气层结决定。大气层结可分为不稳定、稳定和中性三种，在静力中性条件下，风速随高度呈现对数变化规律。

在描述近地面层结构时最常用的就是Monin Obukhov相似理论，即近地面层的无量纲风速[6]：

式中：z为距离地面高度，m；u为z高度上平均水平风速，m/s；κ为von Karman常数，通常取值0.4；u*为摩擦风速，m/s；L为Monin-Obukhov长度，m。

在中性层结时，L→∞，φ(Z/L)→1，对式(1)进行高度积分，得到各层高度风速：

式中：z0为地表粗糙度。zi为测风塔第i个风速仪安装高度(以下简称“第i层高度”)，m；为 第i层高度在中性大气层结下的风速平均值(以下简称“第i层高度中性风速平均值”)，m/s。考虑到近地面层高度通常不超过100m，因此本文规定zi≤100 m。

将测风塔第一层高度中性风速平均值作为参照风速，可以由式(2)得：

相邻高度风速之间两两相减，得：

由式(3)可以看出若粗糙度z01、第一层高度z1和与之对应的中性风速平均值确定，则与lnzi呈线性相关关系。采用最小二乘法拟合散点(lnzi,)得到斜率a，即：

考虑到参照风速也可选为其他层高度风速，因此分别计算第i层高度风速作为参照风速时与之对应的粗糙度z0i，再对所有z0i求平均值得到最终所求粗糙度z0，即：

1.2 中性风速选取方式

根据测风塔实测数据，可得到不同高度处风速与湍流强度相关曲线，中性温度层结对湍流发生发展基本没有影响[5]，即当湍流强度基本保持不变时，认为热效应的影响消失，大气稳定度可考虑为中性，此时把湍流强度开始变稳定的风速定义为中性大气强风阈值[7]，把湍流强度基本稳定时对应的风速区间本文在此称为中性风速组。

运用式(2)～(7)求取粗糙度的关键点在于各层高度中性风速的选取。测风塔各层高度(其中，zmax≤100 m；zmin建议取20～30 m，若测风塔周边地表环境非常简单，对10 m高度测风数据影响程度很小，则zmin可取10 m)的中性风速确定方式为：

1)选取一个或多个完整年时段的10 min实测数据并剔除其中无效数据得到测风塔有效数据，各高度有效数据完整率均需超过90%。

2)根据测风塔有效数据，绘制不同高度处风速与湍流强度曲线，查找相邻风速段平均湍流强度变化率最先满足|ε|≤5%时对应的风速值，该值即为中性大气强风阈值Ur，同时参考我国法规中推荐的修订帕斯奎尔分类法[8]，规定各层高度的Ur不小于6 m/s。为保证中性风速组内对应的湍流强度基本稳定，需在u≥Ur范围内进一步剔除湍流强度变化率|ε|超过10%对应的风速段，此外，考虑到较高风速段(尤其是12～15 m/s以上)因数据量过少可能会造成湍流强度计算结果存在一定不确定性，由此引起的湍流强度变化率过大往往不具有代表性，因此，本文对于风速出现频率小于0.5%的风速段不管湍流强度变化率是否超过10%均予以保留。综上，得到的风速区间即为某层高度的初选中性风速，对应的测风时段即为初选中性风速时段。

3)为保证各高度中性风速测风时段同期，对步骤2)得到的各层高度初选中性风速的测风时段求交集得到共同的测风时段，该测风时段内各高度的风速即为符合中性大气层结条件下的风速组。

2 实例分析

2.1 粗糙度取值对风能资源评估的影响

一般在风电场宏观选址阶段，规划区域内基本上仅有少量已建测风塔，有的甚至没有或在附近区域才有已建测风塔，这就很难直接反映整个规划风电场区域的风能资源。对于平坦地形，WAsP软件是目前比较常用的风能资源评估软件，利用WAsP软件推算规划区域风能资源的主要过程为：①对场内或附近测风塔资料进行整理分析；②提取规划风电场及附近区域地形地貌资料；③根据地形地貌资料绘制粗糙度线，并给粗糙度线两侧的地表粗糙度赋值；④将整理分析后的测风数据文件和处理后的地形图文件导入软件中进行相关计算，即可得出规划区域风能资源分布。可见，地表粗糙度取值对风能资源的评估结果至关重要。为了说明此问题，本文选用某一平原区域来分析粗糙度取值对风能资源评估的影响。

2.1.1 研究区域概况

研究区域地面高程为3～5 m，地形平坦，基本地貌特征如图1所示，整个区域由北至南地表粗糙度逐渐增大，风速由沿海向内陆地区呈衰减趋势。

图1 研究区域地形地貌及测风塔分布图

区域①为大片鱼塘养殖区，区域②为以农田为主的开阔地貌，区域③为有一些房屋和树木的田野，区域④为有较多房屋和树木的田野，区域⑤为有很多房屋接近市郊的田野，其余区域⑥均为有密集建筑群的小城镇。现阶段在区域②～⑤这四个区域内均建有一座测风塔，各测风塔基本信息如表1所示。

表1 测风塔基本信息表

2.1.2 粗糙度取值影响分析

对于平坦地形，可依照文献[9]粗糙度划分方法，区域①粗糙度选择为0.02 m，区域②粗糙度选择在0.03～0.10 m，区域③粗糙度选择在0.05～0.20m，区域④粗糙度选择在0.10～0.30 m，区域⑤粗糙度选择在0.15～0.40 m，区域⑥粗糙度均选择为0.6 m。

在上述粗糙度值划分方法区间内，采用两种设置方案。方案一：选择区域②粗糙度值为0.08、区域③粗糙度值为0.1、区域④粗糙度值为0.12，区域⑤粗糙度值为0.15；方案二：区域②粗糙度值为0.03、区域③粗糙度值为0.15、区域④粗糙度值为0.2、区域⑤粗糙度值为0.4。以1#测风塔100 m高度实测资料作为输入条件，模拟出的整个区域风速分布规律如图2～图3所示。

图2 方案一的风速等值线分布

图3 方案二的风速等值线分布

由图2～图3可知，随着粗糙度取值的不同，整个区域的风能资源计算结果差别也较大。方案一情况下(即各区域的粗糙度相对变化值较小时)，区域③模拟风速变化范围在6.1～6.3 m/s，自北向南风速衰减了0.2 m/s；区域④模拟风速变化范围在6.0～6.15 m/s，自北向南风速衰减了0.15 m/s。方案二情况下(即各区域的粗糙度相对变化值较大时)，区域③模拟风速变化范围在5.7～6.3 m/s，自北向南风速衰减了0.6 m/s；区域④模拟风速变化范围在5.6～5.9 m/s，自北向南风速衰减了0.3 m/s。可见，随着区域间粗糙度相对变化值的增大，区域③和区域④的风速变化范围增大，整个区域自北向南风速衰减速率加快。

2.2 本文推荐方法计算结果分析

继续以上述研究区域为例，采用本文提出的新的地表粗糙度计算方法来推求各局部区域的粗糙度值，并通过模拟结果来验证本文推荐方法的可行性和可靠性。

2.2.1 中性风速组的选取

本文收集到的4个测风塔基本信息如表1所示，各测风塔数据均为10 min实测数据，其中1#和2#测风塔的测风年时段同期，3#和4#测风塔的测风年时段同期。采用1.2节提出的中性风速组选取方式来确定各测风塔的中性风速组。

以2#测风塔100 m和30 m两个高度为例，100 m和30 m高度中性风速组初选过程如图4所示。

图4 2#塔100 m 和30 m高度中性风速组初选过程

由图4(a)可知，100 m高度中性大气强风阈值为8 m/s；在8～14 m/s风速段内各级风速出现频率均大于0.5%，且该范围内12 m/s风速段湍流强度和13 m/s风速段湍流强度的变化率超过10%，因此剔除风速区间[12 m/s, 12 m/s]后，得到2#塔100 m高度初选中性风速组为[8 m/s,12 m/s)∪(13 m/s, 20.6 m/s]。

由图4(b)可知，由于30 m高度湍流强度变化率最先满足|ε|≤5%时对应的风速值小于6 m/s，因此，30 m高度中性大气强风阈值直接取为6 m/s；在6～12 m/s风速段内各级风速出现频率均大于0.5%，且该范围内没有湍流强度变化率超过10%的风速段，因此，2#塔30 m高度初选中性风速组为u≥6 m/s。

以此类推，2#塔70 m高度初选中性风 速 组 为 [7 m/s, 11 m/s)∪(12 m/s, 19.3 m/s]；90 m高度初选中性风速组为[7 m/s, 12 m/s)∪(13 m/s, 19.8 m/s]。

在上述得到的各高度初选中性风速组的基础上进一步筛选出同期时段的数据即为2#塔最终所求中性风速组。

同理，其余1#、3#和4#三座测风塔的中性风速组分别为：

1#塔30 m高度初选中性风速组为u≥6 m/s，50 m～90 m高度初选中性风速组均为u≥7 m/s，100 m高度初选中性风速段组u≥8 m/s，以此为基础筛选出各高度同期时段数据即为1#塔终选中性风速组。

3#塔30～85 m高度初选中性风速组均为u≥6 m/s，100 m高度初选中性风速组为[7 m/s,11 m/s)∪(12 m/s, 19.9 m/s]，以此为基础筛选出各高度同期时段数据即为3#塔终选中性风速组。

4#塔30 m和50 m高度初选中性风速组均为u≥6 m/s，70m高度初选中性风速组 为 [6 m/s, 9 m/s)∪ (10 m/s, 16.9 m/s]，85 m高度初选中性风速组为[6 m/s, 10 m/s)∪(12 m/s, 17.4 m/s]，100m高度初选中性风速组为[7 m/s, 11 m/s)∪(12 m/s, 18.1 m/s]，以此为基础筛选出各高度同期时段数据即为4#塔终选中性风速组。

1#～4#四座塔各高度中性风速组的平均值如表2所示。

表2 各测风塔各高度中性风速平均值表

2.2.2 粗糙度计算结果

以表2中 4个测风塔中性风速平均值计算结果为基础，采用最小二乘法分别拟合4个塔的散点(lnzi,)，各测风塔拟合结果图如图5所示。

图5 各测风塔散点(lnzi,)拟合结果图

根据图5，1#塔拟合直线斜率a=1.304 3，结合式(7)和表2结果求得1#塔对应的区域②粗糙度值为0.040 m。以此类推，分别求得区域③粗糙度值为0.095 m，区域④粗糙度值为0.216 m，区域⑤粗糙度值为0.347 m。

以上各区域的粗糙度计算值均与文献[9]粗糙度划分方法结果基本吻合。

2.2.3 风资源模拟结果分析

考虑到1#和2#测风塔测风年时段同期，3#和4#测风塔测风年时段同期，因此，本文分别以1#和3#塔的100 m高度实测值作为输入条件，采用WAsP软件分别计算对比在方案一、方案二以及本文推荐方法这3种粗糙度取值方法下2#和4#塔位置处的计算结果。各方案下，1#塔推求2#塔对比结果如表3所示，3#塔推求4#塔对比结果如表4所示。

表3 1#塔推求2#塔的计算结果分析表

表4 3#塔推求4#塔的计算结果分析表

根据表3，当采用方案一时，2#塔计算风速与实测值相差0.11 m/s，风速误差为1.82%，风功率密度误差为5.86%；当采用方案二时，2#塔计算风速与实测值相差-0.13 m/s，风速误差达-2.15%，风功率密度误差达-8.11%；当采用本文推荐的方法时，2#塔计算风速与实测值仅相差-0.03 m/s，风速误差仅为-0.49%，风功率密度误差仅为-2.25%。

根据表4，当采用方案一时，4#塔计算风速与实测值相差0.24 m/s，风速误差达4.31%，风功率密度误差达14.04%；当采用方案二时，4#塔计算风速与实测值相差-0.07 m/s，风速误差为-1.26%，风功率密度误差为-4.49%；当采用本文推荐的方法时，4#塔计算风速与实测值仅相差-0.02 m/s，风速误差仅为-0.36%，风功率密度误差仅为-1.69%。

综上，采用本文推荐方法计算的粗糙度值准确性较高，风能资源模拟结果更加接近真实值。

3 结语

地表粗糙度是影响陆上风电场风能资源评估的一项重要因素。以往在实际工程中，一般常用文献[9]平坦地形粗糙度划分法选取粗糙度值，但该方法是根据定性描述的地形特点来划分粗糙度值，易受个人主观影响，粗糙度取值存在较大不确定性，影响风能资源评估准确性。鉴于此，本文提出一种仅以测风塔实测数据为基础计算地表粗糙度值的方法，通过实例验证，采用新方法计算的地表粗糙度值准确性高，以此模拟的风能资源更加接近真实值，较大程度地降低了粗糙度取值的不确定性，而且计算所需基础资料也简单易获取。