中国的风能资源大多分布在丘陵和山区。因此，在风电场选址时普遍采用计算流体力学(CFD)软件进行复杂地形风能资源评估，为风机的合理排布提供依据。由于风能资源评估需要获得气候平均态的平均风速和平均风功率密度分布，CFD技术采用基于Navier-Stokes大气湍流动量方程的有限元数值求解，计算过程复杂，对计算资源耗费很大，不可能根据拟建风电场1 a的观测资料逐时进行CFD数值模拟。因此，目前应用CFD软件模拟复杂地形风场的基本方法是：首先模拟中性大气层结不同风向条件下的风速分布，称为定向计算；然后用测风塔的实测风速对定向计算的风速分布结果进行调整。由于不同的垂直风速廓线输入导致的CFD计算结果差异很大，如果采用不同大气稳定度条件下的垂直风速廓线进行定向计算，会显著提高复杂地形风场CFD数值模拟的准确度。怎样确定不同大气稳定度条件下的垂直风廓线及其相应的湍流通量参数，是进行风电场风能资源CFD数值模拟的关键问题。

边界层大气科学实际上是一门实验科学，大气边界层湍流理论是建立在外场观测实验基础上的。1945—1951年苏联地球物理研究所前后开展了4次大气边界层湍流观测试验，在距离地面15 m的高度范围内，设置了6层湍流风速、温度等观测。Monin等(1954)在此实验的基础上共同创立了适用于50 m高度以下的近地层大气的莫宁-奥布霍夫相似理论，导出了近地层大气的平均风、温廓线表达式。1968年美国空军在堪萨斯西南部的小麦田中进行了大气边界层观测实验，在32 m的观测塔上共设置了7层湍流观测。Wyngaard等(1971)通过该实验将相似理论推广到了自由对流层，得出了气象参数平均廓线特定的幂律形式，构成了通用于对流边界层的对流相似原理。Businger等(1971)同样在此实验基础上研究了各个湍流量之间的关系，并得出相应的经验公式。 1977年9月至1979年2月，Nieuwstadt等(1984)在荷兰Cabauw进行了一系列的夜间观测实验，实验选择晴朗的夜晚进行，在日落之后约2—3 h开始观测直到日出，使其有一个稳定的湍流边界层。风速观测高度分别为10、20、40、80、120、160和200 m。在已有研究成果的基础上，首次提出了稳定条件下的局地相似理论，并通过观测数据验证，成功地将其应用于夜间稳定边界层。在1987和1988年的3月，澳大利亚弗林德斯大气和海洋科学研究所在斯宾塞湾进行了两次飞机探测实验，飞机配备了一组全面的工具，可以测量高分辨率的空气温度、湿度、表面温度和三维风要素，以及垂直湍流通量等，用于研究整个沿海地区大气湍流特性。Shao等(1990)利用该实验观测资料将局地相似理论应用于存在平流的非均匀下垫面大气边界层，结果表明局地相似理论同样适用于对流层和混合层。Gryning等(2007)采用丹麦116 m观测塔、160 m灯塔和250 m电视塔上的大气湍流观测资料，将近地层垂直风速廓线关系延伸到全边界层。其研究结果表明，在近地层，垂直风廓线满足莫宁-奥布霍夫相似理论，相同稳定度条件下，莫宁-奥布霍夫长度随高度线性增加；近地层以上长度特征尺度与稳定度有关，但与高度无关；大气边界层顶的长度特征尺度可以忽略。最终Gryning导出了全边界层的风廓线表达形式。

中国同时开展具有梯度和通量观测的大气边界层观测实验比较少，有些研究工作采用了国际上的实验数据。例如，胡隐樵等(1993)使用日本气象厅213 m观测塔的大气湍流观测资料验证了局地相似理论；张强等(1994)使用瑞典在农田中24 m气象观测塔上的大气湍流观测资料，验证了近地面层大气的局地相似性，并指出在非均匀下垫面近地面层大气中使用局地相似理论比使用莫宁-奥布霍夫相似理论更为合适。中国在大气环境影响评价和风工程应用领域中，主要采用指数函数表达风速垂直廓线，廓线指数是经验常数，与地表粗糙度、大气稳定度有关(李宗恺等，1985； GB/T 13201-91；GB 50009-2001；GB/T 18710-2002)。这些经验方法虽然能够描述近地层风速廓线分布，但是不能给出CFD计算所需的莫宁-奥布霍夫长度等近地层湍流参数。

综上所述，很多近地层观测实验结果证明湍流相似理论适用于各类大气稳定度的近地层，但是如何表达各类大气稳定度条件下地面至100 m的风速垂直分布，为复杂地形的风场CFD数值模拟提供近地层环境入口风速廓线，是一个需要研究的问题。本研究拟采用中国国家气候中心2009—2010年在内蒙古自治区锡林浩特开展的100 m铁塔大气湍流观测实验数据，分析不同稳定度条件下地面至100 m的垂直风廓线变化规律，并给出各类层结稳定条件下垂直风廓线的数学表达式。

2009—2010年中国国家气候中心在内蒙古自治区锡林浩特市东北部开展了为期1 a的近地层大气湍流观测实验。图 1为实验场区周边10 km范围内的地形图，采用3°分带的高斯-克吕格投影的西安80大地坐标系(单位：m)，其中五角星表示测风塔所处位置，塔高100 m。实验场区10 km范围内地形高度变化平缓，海拔高度整体在1100 m以上，从东南到西北、宽约3—6 km的带状区域海拔高度在1150 m以上，东南角约5 km×5 km区域地势最高，海拔高度达到1200 m。测风塔周边约1 km×1 km范围内地势平坦、开阔，位于海拔高度1150 m的平台上。观测实验区为草原下垫面，夏季草高30—40 cm，覆盖率70%—80%，冬季有高度5 cm左右的枯草。

图 1 测风塔的位置及其周边地形(单位：m) Fig. 1 Location of the observation tower and the terrain nearby

图选项

使用的资料时段为2009年8月15—31日、2010年 4月1—24日。观测分为两种：平均量梯度观测和湍流通量梯度观测。风速、温度和湿度的平均量观测仪器采用江苏省无线电科学研究所生产的ZQZ-TF风速、风向传感器和气温传感器，以及芬兰维萨拉公司生产的HMP45D气温、湿度传感器，观测高度分别为4、10、30、50、70和100 m(表 1)。湍流通量观测采用美国Campbell公司生产的CSAT3三维超声风温仪和美国Li-COR公司生产的CO₂/H₂O分析仪(Li7500)，观测高度分别为10、30、50和70 m(表 2)。

考虑到受仪器故障、检修或天气因素以及供电等的影响会出现数据缺测和异常的情况，首先对观测资料进行了如下质量控制：

(1)剔除不符合物理常识的数据和不符合一般气候统计特征的数据。

(2)把大于5倍方差的数据视为野点剔除。

(3)剔除由于传感器上的水汽凝结等原因导致的数据接收系统和数据传输系统产生的一些随机脉冲。将|X_i+2-X_i|>5σ的值视为随机脉冲进行剔除，其中X为时间序列，σ为方差。

(4)剔除风速小于1.5 m/s和摩擦速度小于0.1 m/s的数据。

锋面等天气过程过境时，也会造成大气边界层湍流运动加剧，使近地层大气不满足定常假设条件。通过查询与观测同期的地面天气观测资料，剔除了2009年8月有两次明显降水过程的资料，即2009年8月19日02—18时(北京时，下同)和8月27日05—20时。2010年4月1—24日没有发现明显的天气过程影响，因此没有进行相应数据处理。

为了分析计算方便，将坐标系的x轴旋转到主导风向，使y方向的平均速度为0，按照右手螺旋法则建立坐标系，旋转时用湍流梯度观测得到的瞬时u，v，w三个方向的风速乘以旋转矩阵

点击浏览公式

为了消除仪器安装不够水平导致的误差，可以通过坐标旋转使平均垂直速度为0。用式(1)中得到的u₁，v₁，w₁三个方向的风速乘以旋转矩阵

点击浏览公式

3 观测资料分析 3.1 大气稳定度

Monin等(1954)在近地层莫宁-奥布霍夫相似理论中，风速和温度随高度的变化为

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

分别采用以下两种方法计算莫宁-奥布霍夫长度(L)：(1)通量法，根据风速和温度脉动观测数据，采用涡度相关法计算各高度的摩擦速度(u_*)、感热通量(w′T′)和莫宁-奥布霍夫长度(L)；(2)梯度法，根据两个高度层上风速和温度观测，采用迭代法计算莫宁-奥布霍夫长度(L)，Businger等(1971)得到的不同大气稳定度条件下的近地层风速梯度和温度梯度表达式为

点击浏览公式

点击浏览公式

表 3和4为观测实验期间分别采用通量法和梯度法计算得到的L在不同稳定度条件下的分布，其中稳定度分类参照Gryning(2007)和Pea(2009)。可以看出，无论是在4月或是在8月、无论在哪一高度层以及无论是用通量法还是梯度法计算，大气层结非常不稳定的出现频率是最高的(43%—67%)；中性、弱稳定和弱不稳定出现频率相对非常小。通量法计算得到的非常稳定大气层结出现频率高于稳定层结；梯度法计算表明稳定层结出现频率略高于非常稳定层结。总体而言，两种方法计算得到的大气稳定度出现频率分布基本相同。表 3和4中的结果与欧洲一些观测结果大不相同，图 2是Gryning等(2007)在丹麦Hovsore、德国汉堡和Holtslag(1984)在荷兰Cabauw观测得到的不同稳定度出现的频率。丹麦Hovsore地处西部沿海，稳定层结比不稳定层结出现概率高，非常不稳定和非常稳定出现频率相对较低；而地处德国汉堡的测风塔离市区5 km，不论风的来向是郊区还是市区，中性层结占大多数，不稳定层结比稳定层结出现概率高；荷兰Cabauw实验场地处于平坦乡村，非常不稳定和中性出现频率低，其他各类稳定度出现频率相差不大。考虑到上述文献中的测风塔都地处平原，而锡林浩特观测试验场区处于高原台地，统计了各类稳定度下的风向频率(图 3、4)。发现不稳定和非常不稳定出现时的主导风向为偏西风，测风塔处于局地地形的东侧，说明地形对近地层湍流交换起到了加强的作用。

图 2 部分文献中不同稳定度出现频率(a.Gryning，b.Holtslag) Fig. 2 Frequency of the different stability classes as shown in some literatures(a.Gryning，b.Holtslag)

图选项

图 3 2009年8月15—31日不同稳定度条件下的风向玫瑰图 Fig. 3 Wind rose under the different stability classes during 15-31 August 2009

图选项

图 4 2010年4月1—24日不同稳定度条件下的风向玫瑰图 Fig. 4 As in Fig. 3 but during 1-24 April 2010

图选项

根据表 3和4中各类大气稳定度的出现频次，以下将稳定度划分为稳定(0<L<200)、弱稳定(L>200)、不稳定(L<-200)和弱不稳定(-200<L<0)共4类。分别采用10、30、50、70、100 m高度的10 min平均风速除以10 m高度的平均摩擦速度，得到不同稳定度条件下的归一化风速随高度的变化。同时，根据Businger等(1971)提出的非中性层结的垂直风廓线表达式(14)—(17)，计算得到以上4类大气稳定度条件下平均归一化风速随高度的变化曲线。

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

选取中性层结条件下观测得到的两层风速计算出粗糙度长度(z₀)，Stull(1988)将1/L在-0.005—0.005看作近中性层结状态，根据10、30、50 m高度的1/L进行筛选，选出3个高度的1/L计算值都在-0.005—0.005范围内的同步观测数据，结果2009年8月共得到36组，2010年4月共得到55组(表 5)，然后选择4和10 m两个高度的观测风速计算Z₀

点击浏览公式

图 5a和b分别为2009年8月15—31日和2010年4月1—24日在稳定、弱稳定、弱不稳定和不稳定层结条件下，观测得到平均归一化风速垂直分布与根据Businger等(1971)经验公式计算结果的对比。从图 5可以看出，在不稳定、弱不稳定和弱稳定层结条件下，观测值与经验公式计算值比较相近，观测值均略有偏大；在稳定条件下，10和30 m高度的观测值与经验公式计算值吻合很好，之后随着高度的升高，观测值与经验公式计算值的偏差越来越大，归一化风速观测值随高度的增加远没有计算值增加得快。这是由于Businger经验公式仅适用于常通量层，大气层结越稳定，常通量层越低，因此图 5中稳定层结的偏差最大。不稳定条件下，大气运动比较剧烈，有利于上下层空气的混合，使得大气特性上下趋于均匀；而在稳定条件下，上下层交换比较少，大气特性在垂直方向差异较大，具有分层的特性。从2009年8月20日和2010年4月16日3 h间隔的小时平均风速廓线日变化(图 6)可以看出，风速廓线形状非常复杂，尤其是在夜间大气非常稳定时，经常出现30 m高度以上风速随高度加速很慢、甚至减小的情况，难以用单调函数整体表达，需要分层次逐一描述。事实上，在中国气象局风能资源专业观测网的内蒙古72座70和100 m测风塔的2009—2011年观测数据中，有大量稳定层结条件时风速随高度升高增大很慢的事实。

图 5 2009年8月15—31日(a)和2010年4月1—24日(b)平均归一化风速u/u*10随高度的变化 Fig. 5 Variation of the average normalized velocity u/u*10 with the height during 15-31 August 2009(a) and 1-24 April 2010(b)

图选项

图 6 相邻两高度层之间摩擦速度和感热通量变化不超过10%的出现次数统计 Fig. 6 Diurnal changes of the averaged wind profile in terms of every 3 h(a. 20 Augst 2009，b.16 April 2010)

图选项

一般认为，满足常通量假设的条件是近地层大气上下层的湍流通量变化不超过10%。通过对不同高度的动量通量和感热通量的统计分析发现，实验期间10—50 m或10—70 m高度层内，都没有出现过摩擦速度或感热通量的变化率不超过10%的情况，这意味着常通量层高度不超过30 m。表 6给出了相邻两高度层之间摩擦速度和感热通量变化不超过10%的出现频次，可以看出，摩擦速度满足常通量假设的占19%—23%，感热通量满足常通量假设的占12%—14%。在满足常通量假设的数据中，稳定层结占16%—34%。进一步表明了在稳定层结大气在垂直方向上由于湍流交换较弱，显示出分层的特性。

表 6 京津冀地区雷暴大风天气过程 Table 6 The number of occurrences of the event in which the variations of friction velocity and sensible heat flux between two adjacent heights is less than 10%

表选项

Nieuwstadt(1984)局地相似理论假定空间任一点都具有其局地湍流特征尺度，定义局地莫宁-奥布霍夫长度为

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

分别采用10、30、50和70 m高度的风速、温度脉动量观测，计算30 min平均的摩擦速度(u_*)和感热通量(w′T′)以及局地莫宁-奥布霍夫长度尺度(Λ)，再分别求取10与30、30与50、50与70 m之间的平均莫宁-奥布霍夫长度尺度，然后用各层平均高度除以平均莫宁-奥布霍夫长度尺度得到各层的大气稳定度(ζ)。此外，分别采用10、30、50和70 m高度的风速和气温30 min平均量观测数据，根据式(23)和(24)分别计算10—30、30—50和50—70 m层的无因次风速梯度。图 7和8分别为min平均的10—30、30—50和50—70 m层的无因次风速梯度随稳定度的变化，可以看出，随着高度的升高，无因次梯度与稳定度的函数关系越来越不清晰。

图 7 2009年8月15—31日的30 min平均的10—30 m(a)、30—50 m(b)和50—70 m层(c)的无因次风速梯度随稳定度的变化 Fig. 7 Dimensionless gradient of velocity(Φm)changes with ζ for 30 min average during 15-31 August 2009 for the layers of(a)10-30 m，(b)30-50 m and (c)50-70 m

图选项

图 8 2010年4月1—24日的30 min平均的10—30 m(a)、30—50 m(b)和50—70 m层(c)的无因次风速梯度随稳定度的变化 Fig. 8 As in Fig. 7 but during 1-24 April 2010

图选项

Businger等(1971)根据Kansas实验数据得到的无因次风速梯度与大气稳定度普适关系的BWIB模式为

点击浏览公式

随后的很多野外近地层观测(胡隐樵等，1993，张强等，1994；王丙兰等，2012)证明了BWIB模式的适用性。本研究采用Mathlab软件对2009年8月15—31日和2010年4月1—24日的30 min平均无因次风速梯度与稳定度的关系，按照10—30、30—50和50—70 m分层进行拟合，表 7中列出了各高度层无因次风速梯度与稳定度关系的拟合公式及其复相关系数。可以看出，每个高度层上的无因次风速梯度与稳定度关系的BWIB拟合复相关系数都在0.76以上，高度越低，复相关系数越大。在大气稳定层结条件下，无因次风速梯度与稳定度的线性相关(β_m)参数随高度的增加而减小，证明了层结稳定大气的分层特性。由于前文图 5中的不稳定大气层结条件下，观测值与式(14)计算值偏差较大。因此针对稳定大气，根据表 7中无因次梯度与稳定度的拟合结果，将式(14)改为

点击浏览公式

点击浏览公式

首先将采用通量法计算得到的局地莫宁-奥布霍夫长度尺度(Λ)代入式(26)—(27)，分别得到了稳定层结条件下2009年8月15—31日和2010年4月1—24日的归一化风速随高度的分布(图 9)。可以发现稳定条件下的归一化风速的观测值与计算值基本吻合。一方面证明了局地相似理论在内蒙古草原近地层大气的适用性，另一方面说明稳定层结条件下，分层描述大气特性可以取得较为准确的结果。

图 9 2009年8月15—31日(a)和2010年4月1—24日(b)通量法计算的平均归一化风速u/u*0随高度的变化 Fig. 9 Variation of the average normalized velocity u/u*0 with the height during 15-31 August 2009(a) and 1-24 April 2010(b)as calculated by the flux method

图选项

考虑到大气湍流观测设备比较昂贵，尚不能纳入常规气象观测业务，而70 m高度的测风塔观测在风电场选址阶段和风电场运营后的风电功率预报业务中属于常规观测项目。因此，只要测风塔的每层都有风速和温度的观测，即可采用梯度法计算得到各高度局地莫宁-奥布霍夫长度尺度Λ、u_*、T_*，然后采用表 7中的无因次风速梯度与稳定度关系的拟合关系系数，代入式(26)—(27)即可求出不同稳定度条件下的风速垂直廓线分布。图 10分别是稳定层结条件下2009年8月15—31日和2010年4月1—24日的归一化风速随高度的变化，可以发现，稳定条件下的归一化风速观测值与计算值同样基本吻合。

图 10 2009年8月15—31日(a)和2010年4月1—24日(b)梯度法计算的平均归一化风速u/u*0随高度的变化 Fig. 10 As in Fig. 9 but by the gradient method

图选项

随着中国风电产业的快速发展，越来越多的山地风能资源得到开发利用，风场的计算流体力学数值模拟计算(CFD)已成为风电场选址风能资源评估的必要技术手段。复杂地形风场微尺度数值模拟的CFD技术应用方法也不断地发展，由原来的采用中性大气条件进行初始定向计算，发展到对不同稳定度条件分别进行定向计算，显著提高了复杂地形风场微尺度数值模拟的准确率。本研究根据复杂地形微尺度数值模拟的需要，通过分析中国国家气候中心内蒙古草原大气湍流梯度观测实验数据，研究了不同稳定度条件下的近地层大气风速垂直廓线的表达方法和湍流通量计算，主要结论如下：

(1)内蒙古草原近地层大气中，不稳定层结和稳定层结条件发生的概率远高于中性和弱不稳定、弱稳定层结的发生概率。春季不稳定层结发生概率高达54%—56%，稳定层结发生概率29%—35%；夏季不稳定层结发生概率高达57%—60%，稳定层结发生概率28%—34%。

(2)在距地面100 m高度的近地层内，动量通量和感热通量与大气稳定度的关系满足局地相似原理。Businger等(1971)的BWIB模式可以较好地描述不稳定、中性和弱稳定大气层结条件下的无因次风速梯度与稳定度的关系，但在稳定层结条件下，随着高度的增加偏差显著增大。

(3)稳定层结时大气特性在垂直方向上变化较明显，需要分层考虑。根据内蒙古锡林浩特草原大气湍流观测资料得到，稳定层结时无因次风速梯度随稳定度(ζ)的变化关系可以表达为1+β_mζ，其中β_m的取值与距地面高度有关：30 m高度及以下，β_m取4.3—5.4；高度在(30 m，50 m］，β_m取2.9—3.4；50 m及以上高度，β_m取2.0—2.9。

建议复杂地形风场微尺度数值模拟方法为，首先在铁塔观测时，不仅需要设立风速观测，还需要温度观测；CFD初始定向计算按照不稳定、弱不稳定、中性、弱稳定和稳定分别进行，然后将测风塔观测资料按照不同稳定度分类，分别与相应的定向计算结果进行综合分析；稳定层结定向计算的初始风速廓线和湍流通量参数需在不同的高度层采用不同的拟合函数。

对于不同的地形和下垫面类型，如戈壁、海滩、长满灌木的丘陵等，或是在冬季和秋季，风速垂直廓线的分布规律是否仍然可以用局地莫宁-奥布霍夫相似原理分层表达，还需要进一步开展观测实验研究。