蒙西地区不同边界层参数化方案的近地层风场预报效果评估*

2018-08-08李玉鹏王东海尹金方

中山大学学报(自然科学版)(中英文) 2018年4期

李玉鹏，王东海，尹金方

(1.中国气象科学研究院灾害天气国家重点实验室，北京100081；2.中山大学大气科学学院，广东广州 510275；3.吉林省气象科学研究所，吉林长春 130062)

可再生能源的使用为社会和经济的发展供应能源，并且减少化石燃料对环境和人类健康的负面影响[1]。在众多可再生资源中，风能是一种清洁且储量丰富的可再生能源。据估计全球的风能约为2.74×109MW，其中可利用的风能为2×107MW，具有很高的开发利用价值[2]。现阶段，全世界的风电装机容量以每年30%的速度增长[3]。根据国家能源局发布的国家可再生能源发展目标，至2050年左右，风电在国家总体电源结构中所占的比例将达到1/5左右，成为我国电力供应中的主力电源之一[4]。

风是大气中最基本的自然现象，可以表征大气的流动性和能量特征。风电产业不断发展，再加上极端天气事件频发，风场的模拟和预报变得越来越重要。早期针对风的预测方法很多是基于概率和数理统计的，包括时间序列法(ARMA)、模糊逻辑法(Fuzzy Logic)、人工神经网络法(ANN)、卡尔曼滤波法(Kalman Filters)等[5-8]，但风受地形、天气系统演变、辐射、地面加热以及海陆热力差异等多种因素作用，而以统计为基础的预测模型缺少影响风速变化的动力学信息，因而预报时效较短,预报准确率低，只有通过求解大气运动物理方程组为基础的数值天气预报(NWP)才能较好地提高风电预测的准确率[9]。近年来WRF(Weather Research Forecast)等中尺度天气数值预报模式趋于成熟，现已较为广泛地运用于风能资源评估领域[10]，中尺度数值模式在风电预报中的应用也得到国内众多学者的关注[9, 11-15]。与此同时WRF等模式针对大气运动中的各种物理过程已经发展出一系列参数化方案，包括微物理参数化方案，辐射方案(长波方案和短波方案)，积云过程方案以及与边界层有关的方案：表面层(SL),行星边界层(PBL)和陆面过程模式(LSM)等等。各种物理过程参数化方案都是建立在一定的假设之上，不仅彼此之间是非线性的关系，而且与模式的动力核心之间也是非线性的。这参数化方案的选择会给预报结果带来重要的影响[16]，特别是边界层参数化方案。

风机轮毂高度一般在70～100 m，该高度范围内的近地层风场演变是边界层内的天气现象，近地层风场的预报效果很大程度上依赖于模式边界层参数方案的合理性[9]。在植被稀疏、地形复杂的环境下，模式边界层的参数化是决定模式能否较好模拟出该地区近地层风场特征的关键[17]。大气边界层物理过程与下垫面物理性质具有紧密关系，蒙西地区独特的地理位置和地形地貌会对该地区的边界层特性产生影响。

本文研究的蒙西地区处于37-46°N，97-120°E，以蒙古高原为主体，同时具有较为复杂多样的形态。平均海拔在1 000 m以上，北部地势较为平坦，中部自西向东分布着阴山山脉(狼山、色尔腾山、大青山、灰腾梁)，西南部还有贺兰山山脉，这些山脉共同构成内蒙古高原地貌的脊梁。该地区属典型的中温带季风气候，具有降水量少而不匀、寒暑变化剧烈的显著特点[18-19]。除此之外，大气边界层物理过程也受到天气尺度背景场的制约与影响[20]，蒙西地区常年受到蒙古高压的作用，这也影响着该地区的边界层特性。下垫面和天气背景的共同作用使得该地区的边界层特性有独特的特点，边界层内的风速明显高于其他地区，其年平均风速在3 m/s以上，蕴藏丰富的风能资源[21]，合理地选择边界层参数化方案对于蒙西地区近地层风场的预报至关重要。

虽然国内外关于边界层参数化方案对近地形风场影响的研究并不少见[22-27]，但大多是基于10 m风场的天气学分析，与风电预报结合较少，而且这些研究很可能得不到“普适”的结论，因为在某一局地模式中与本地结合最好的参数化方案很有可能在其他区域完全不适用[28-29]。如果不对特定的地区进行边界层参数化方案的试验，只靠以往的研究很难选出最优的方案，因此有必要对蒙西地区进行不同的边界层参数化方案数值试验，从而利用该地区丰富的风塔观测资料对结果进行评估检验，以便为该地区的风资源预报服务提供参考。

1 资料与研究方法

1.1 观测资料简介

本文用于检验模式预报效果的观测数据为蒙西地区70个国家级地面气象站的10 m风观测资料，和蒙西地区18座测风塔的10 m和70 m风向风速资料。地面站和测风塔的分布如图1所示，可以看出，蒙西的中部站点分布较密，东部站点较稀疏。

1.2 边界层参数化方案简介及试验设置

现阶段WRF中的行星边界层参数化方案有10多种，主要有MRF方案[30]、YSU方案[31]、ACM2方案[32]、TEMF方案[33]、MYJ方案[34]、QNSE方案[35]、BouLac方案[36]、MYNN2方案[37]、MYNN3方案[38]以及UW方案[39]，其中前4种边界层方案为非局地方案，其它为局地边界层方案，详见表 1。本研究将选取其中7种常用的方案进行试验。

表1 WRF模式中常用的边界层参数化方案Table 1 The common boundary layer parameterization schemes used in WRF model

本研究采用WRF-ARW 3.7.1模式对蒙西地区2015年5月的低层风场进行回报，预报每天世界时12时冷启，一次预报36 h。试验中的初边条件来自于美国环境预报中心的全球预报系统提供的 0.5°×0.5°的全球预报场，每6 h更新一次边界值。模式采用12 km和4 km的双层嵌套，两重模拟区域的水平网格数分别为448×373和601×403，中心经纬度为102.6°E和41°N，模拟区域如图1所示。试验中还利用GSI系统同化地面观测、探空、飞机报等资料。模式的物理过程参数化方案选择如下：短波辐射采用Dudhia方案[40]；长波辐射采用RRTM方案[41]；微物理过程采用Lin方案[42]；陆面过程采用Noah方案[43]；积云对流采用KF方案[44]，在此基础上采用7种不同的边界方案进行回报试验，7种方案分别为YSU、MYJ、MYNN2、QNSE、ACM2、UW以及TEMF，其中YSU、ACM2和TEMF为非局地闭合方案，其余为局地闭合方案。近地层方案与边界层相匹配的，所对应近地层过程参数化方案见表 2。

表2 本文试验选取的7种边界层试验方案及其近地层方案Table 2 Physical parameterization configuration of the simulation

图1 模式区域设置及测风站点分布Fig.1 Simulation domain and spatial distributions of wind observation stations in western Inner Mongolia

1.3 研究方法

本研究主要通过一些统计量来检验不同边界层对风场的预报情况，主要的统计量如下：

风速误差

ERR=xf-x0

风向误差为公式(2)，因为风向是0°～360°的循环的标量，其误差的绝对值不会超过180°，所以与风速误差的公式有所区别。

绝对误差为，

AE=ERR

相对误差为，

有了上面的基础定义，可以得到平均误差(ME)、平均绝对误差(MAE)、平均相对误差(MRE)、均方根误差(RMSE)和相关系数(CORR)等一系列统计量。

这些公式中，N为观测数量，x0为观测值，xf为预报值。

Weibull分布是用于描述风速分布的概率函数，风速v的Weibull分布概率密度函数表达式为，

其中k为形状参数，是一个无量纲量，λ为尺度因子，单位为m/s，两个参数可用如下公式来计算

其中，σ和μ分别表示风速的标准差和平均值。

为了研究风场预报效果与风切变之间的关系，定义风切变指数[45]为

其中，v(z1)和v(z2)分别表示z1和z2高度处的风速，本文z1和z2分别取为10 m和70 m。风切变指数α为无量纲数，一般分为5个等级，分别为强稳定、稳定、中性、对流及强对流，相应的α取值为分别为大于0.3、0.2～0.3、0.1～0.2、0～0.1及小于0。

2 结果与分析

2.1 试验的整体效果分析

为研究不同边界层方案在蒙西风场预报中的表现，对2015年5月份蒙西地区风场进行回报试验，选取每次预报的第12～35 h共计24个时次的预报，然后采用反向距离方法将风场预报结果插值到观测站点以便与观测进行比较。首先分析近地层10 m高度(图 2)和70 m高度(图 3)的站点观测和模拟结果的平均风速和平均风向的逐小时变化曲线。10 m站点观测的平均风速为1～10 m/s，70 m的站点观测平均风速则在2～14 m/s。对风向而言，10 m跨度为50°～350°，70 m则为100°～350°，两者相差不大。就实况风速的变化而言，10 m风速相对于70 m的风速脉动性更强，风速值的峰谷交替频率更高，且峰谷的起伏更为明显，这种情况在1～20日时段内尤为突出，20～30日时段内波动强度变弱。各边界层方案对10 m风速都有一定的高估，且风速起伏越明显的地方，这种高估也就更为明显，这种系统性误差在前人的研究中也有所揭示[46]。模式对风向的把握整体较好，10 m和70 m都没有明显的偏差。总体而言虽然模式的预报结果抓住了实际风速风向的变化趋势和位相变化，但是不同的边界层方案之间存在明显差异，10 m风速的差异尤为突出。

图2 2015年5月试验时段内观测与模拟的10 m风速和风向随时间的变化Fig.2 Time series of the observed and various BL scheme-simulated speed & direction for 10 m wind at in May 2015

图3 2015年5月试验时段内观测与模拟的70 m风速和风向随时间的变化Fig.3 Time series of the observed and various BL scheme-simulated wind speed & direction at 70 m in May 2015

为了更直观地分析各方案对风向的模拟效果，分别对10 m(图 4)和70 m(图 5)的风场模拟结果绘制风玫瑰图与实况进行对比。由图 4中可知，在试验期间10 m风场频率较高的风向在第Ⅱ和第Ⅲ象限，频率多在7%以上，而第Ⅰ和第Ⅳ象限风向频率较低，都在4%左右。主导风向为W、WNW和NW，这三者出现的频率都在10%左右，差距不大，其中WNW方向的频率是三者中稍高的。各方案都抓住了实际风向的整体分布特征，但一些细节方面跟实况仍有一些差异。各试验方案模拟的主导风向跟观测相比有顺时针偏差，其中W风较观测顺转了一个方位且频率偏大，同时与WNW风的频率差增大，其中ACM2、MYJ、MYNN2和UW等方案表现最为明显。而在第Ⅰ、Ⅳ象限的非主导风向区内从NE沿顺时针到SSE这个区间内，各试验方案对其风向频率都有所低估。综合来看，YSU和TEMF两种方案对10 m风向的把握要优于其余5种边界层方案。图 5中70 m风场的主导风向和10 m的相似，也位于坐标系的第Ⅱ和第Ⅲ象限，非主导风向区为第Ⅰ和第Ⅳ象限。W、WNW和NW是出现频率最大的3个风向，频率都在11%左右，最主要的风向为NW风。各方案模拟的风向的总体分布形势与实况相符，但同样有使主导风向“加强”的趋势，即高估主导风向发生的频率，从而加大与非主导风向之间的频率差，使得模拟的风向更加趋于集中在某一风向范围内。

图4 2015年5月试验时段内10 m高度观测与各边界层方案的风玫瑰图Fig.4 Wind rose plots for the observed and various BL scheme-simulated results for 10 m wind in May 2015

图5 2015年5月试验时段内70 m高度观测与各边界层方案的风玫瑰图Fig.5 Wind rose plots for the observed and various BL scheme-simulated results for 70 m wind in May 2015

表3和表4分别为不同高度风速风向预报结果的统计特征，针对风速的统计量分别为相对误差(RE)、误方根误差(RMSE)、时序相关(time CORR)和站点相关(station CORR)，因为风向是一个循环型的变量，不存在“大小”的概念，所以针对风向不讨论相对误差，只计算其偏差(ME)，其余3个统计量与风速相同。由表 3不难看出，YSU方案对10 m风场预报的平均相对误差(19.34%)和平均均方根误差(1.95 m/s)都很明显地低于其他方案。就各方案的平均站点相关系数而言，MYJ的最高，为0.50，最低的TEMF为0.44，各方案站点相关系数差异不大，各方案对风速的空间分布形势的预报能力相近。MYNN2的平均时序相关系数最大为0.71，最低的是TEMF，为0.59，比MYNN2低0.12，表明TEMF方案对风速的突变的预报能力较弱。而70 m风速预报结果的统计量与10 m风相反，YSU的平均相对误差和平均均方根误差为各方案中最大的，分别为16.87%和3.19 m/s。70 m风速的站点相关系数和时序相关系数都是MYNN2和ACM2的两种方案较高，站点相关系数都为0.56，时序相关都为0.69，而TEMF方案两种相关系数最低，分别为0.50和0.62。如表 4所示，10 m风向的统计结果表明，QNSE方案预报的风向均方根误差最小，为4.66°，误差最大的方案为UW方案，为15.20°，MYJ的平均均方根误差最小，为53.71°，平均均方根误差最大的为TEMF方案，为57.11°。站点相关系数和时序相关系数都是MYNN2方案最大，分别为0.40和0.48，而TEMF方案两种相关系数各方案中都是最小，分别为0.37和0.44。QNSE方案的70 m风向误差最小，为6.12°，是，而UW方案的误差最大，为15.37°。均方根误差则是ACM2最小，为39.34°，最大的为TEMF方案，为43.68。站点相关系数最小的方案都是TEMF，两系数分别为0.45和0.49，而最大的方案分别为UW和ACM2，值分别为0.50和0.55。综合来看，各方案预报的10 m和70 m风向都有顺时针偏差，但各方案的平均风向偏差都在20°以内。总体来讲TEMF方案在风向的预报上表现较差，MYNN2表现较好。

图 6是10 m和70 m观测及模拟的风速频率分布直方图及其对应的Weibull分布拟合曲线。图 6a 中实测的10 m风速众数偏左，即低风速发生的频率较大，高风速出现的频率较小，其形状参数为1.45，尺度因子λ大小为3.85 m/s。而各方案的模拟风速拟合的k值都要大于实际风速，最小值为1.88，最大值为2.09，λ值也都大于实况，从最小值5.37 m/s到最大值6.22 m/s。各方案对小于2.5 m/s的低风速模拟的不是很好，大于2.5 m/s的风速则都表现为一定的高估。YSU方案的拟合曲线形状跟观测最为接近，但是λ值较众值偏右。各试验方案的众值相对观测向右偏移，趋向于正态分布。图 6b 的70 m实际风速分布情况跟10 m的有所不同，其k值为2.02，λ为8.95 m/s，两个参数都比10 m的大，所以相较10 m实测风速整体右移，且形状更趋向于正态分布。各模拟方案的拟合参数相较10 m时与实况也更为接近。通过风速频率分布直方图和Weibull拟合曲线可以发现，各方案在低风速和高风速时表现为对实际风速的低估，在中等风速时则表现为一定程度的高估。

表3 2015年5月试验时段内10 m和70 m 风速预报结果的统计特征Table 3 Different metrics for comparing simulations of different BL schemes with observations at 10 m and 70 m for wind speed in May 2015

表4 2015年5月试验时段内10 m和70 m 风向预报结果的统计特征Table 4 Different metrics used for comparing simulations of different BL schemes with observations at 10 m and 70 m for wind direction in May 2015

图6 2015年5月试验时段内10 m和70 m风速频率分布及Weibull拟合Fig.6 Wind speed distribution and Weibull fitting of simulations of different BL schemes and observations at 10 m and 70 m in May 2015

2.2 风速误差的分段分析

太小的风速不能使风机叶片转动，无法将动能转化为电能，风速太大则有损坏风机的危险，因此风电领域中一般将风速分为4个等级[14]，分别为无效风速(0～3 m/s)、有效风速(3～15 m/s)、满发风速(15～25 m/s)和切出风速(大于25 m/s)，以此为基础分析不同风速区间的各边界层方案的预报效果，可以有针对性地选择边界层参数化方案，并且优化模式的预报效果。这里取前3个风速等级来分析。

图 7为10 m处风速的分段误差统计量结果，可以看到3个风速区间的平均相对误差和平均误差存在很大差异。各方案在无效风速(0～3 m/s)区间的平均相对误差最大，各方案对低风速的预报都有不同程度的高估，其中偏差最小为YSU方案，大小为81.80%，最大为MYJ，是136.87%。在有效风速(3～15)m/s区间内，整体模拟结果的平均相对误差是3个风速区间内最小的，7种方案平均相对误差在6.84%～32.79%，平均相对误差最小的仍然是YSU方案，最大的仍是MYJ方案。而在(15～25 m/s)风速区间内，各方案平均相对误差又有所增大，范围为33.10%～46.72%，这时MYJ表现最好，而YSU表现最差。而10 m风速的平均绝对误差方面，切入风速段和有效风速段(0～15 m/s)平均绝对误差小，最大为2.36 m/s，满发风速段(15～25 m/s)平均误差大，最大为7.30 m/s。图 8为70 m风速的分段检验的结果,图 8b显示切入风速区间各方案的平均相对误差最大，满发风速段次之，有效风速段最小。切入风速区间内各方案的误差在130.79%～156.43%，TEMF最高，QNSE最低。有效风速区间内各方案的平均相对误差在5.36%～14.29%，最低和最高的方案分别为ACM2和YSU。满发风速区间内平均相对误差仍然是YSU方案最高，为31.50，最低的方案是MYNN2，大小为22.44%。各风速段的平均绝对误差和10 m风速类似。

图7 2015年5月试验时段内10 m风速平均相对误差及平均绝对误差的分段统计结果Fig.7 Mean relative error and mean absolute error of wind speed for different wind categories at 10 m in May 2015

图8 2015年5月试验时段内70 m风速平均相对误差及平均绝对误差的分段统计结果Fig.8 Mean relative error and mean absolute error of wind speed for different wind categories at 70 m in May 2015

2.3 风场模拟误差的日变化特征分析

近地层风场变化很容易受到地表特性的影响，尤其是地表加热作用。在天气背景一定时，局地热力特性在一天中的改变必然会影响近地层风场，因此分析各试验方案对风场日变化的预报效果是十分重要的。根据世界时和地方时的关系，文本白天对应的是世界时00-10时，夜间指世界时11-23时。图 9a和b分别为风速和风向月平均得到日变化特征，可见实测的平均风速、平均风向存在很明显的日变化特点，当日22时-次日08时，风速逐渐增大，有大致2 m/s的增幅，08-21时风速逐渐减小。风向也大致呈现出同样的日变化特征，00-08时，平均风向不断顺转，这期间转动的幅度大约为40°，08-14时，风向又较快地逆转回到原本的方向，接着14-23时，风向较为稳定地维持在190°～200°之间，期间有小幅波动。结合图 9中的c和d发现，虽然各边界层方案都很好地抓住了这种风速风向的日变化特征，但各方案模拟得到的月平均风速要比观测偏高，其中最为明显的是MYJ方案，偏高最小的是YSU方案。大部分方案风速的平均相对误差在11时出现最低，且夜晚的误差大于白天的。对风向的模拟大部分方案都表现为正的偏差，即相对于观测顺偏，只有QNSE方案的平均风向偏差在05-12时出现了负的偏差。由图 10可以看到70 m的实测风速风向也存明显的日变化特点，观测的平均风速在00-07时不断增大，07-23时，风速不断减小。大部分的试验方案也都给出了类似的结果，只有YSU方案在00-11时其模拟的风速结果比其它边界层方案偏低，跟观测比约偏低1～1.5 m/s。70 m风向白天与10 m的风向变化相似，但是夜间14-23时平均风向一直为顺转，与10 m时有所差异。另外，结合图 9和图 10综合来看可以发现，在10 m高度上实测平均风速小于各方案模拟的风速，而在70 m高度上的实测平均风速大部分时间都大于各方案模拟的平均风速，这说明10 m和70 m之间的实际风速切变要大于大部分模拟的风速切变，这种情况在白天更为突出。

2.4 稳定度指数分析

大气层结对近地层风场有很大影响，而风切变指数可以反映近地层的大气稳定度。上文风速日变化的分析也发现了实际风切变大于模拟风切变的现象。为了更为细致地分析各方案对风切变指数在试验时段内的模拟情况，图 11给出了观测及各方案在一个月的试验段内逐小时的风切变指数变化特征。由图中可知风切变指数呈现出很明显的日变化特征，在每日白天大气以中性层结为主，有时会出现对流层结；而每日的夜间，稳定层结和强稳定层结占主导。7种边界层方案的模拟结果也都呈现出明显的日变化特征，但是跟实况相比都存在一定误差。在白天，除了TEMF方案外，其余各方案的预报相较观测都趋向于更不稳定层结，MYJ和QNSE表现最为明显，YSU和ACM2最不明显；而夜间各方案模拟的风切变指数大部分方案小于观测，只有MYNN2和YSU的风切变指数稍高于观测。综合来看，TEMF、YSU以及ACM2方案对风切变指数的预报较其他方案稍好，结合前面的分析可知，在稳定度预报较为准确的时刻，10 m风速的预报也比较准确，说明大气层结的稳定性对低层风场的影响更为突出。

图9 2015年5月试验时段内10 m各方案模拟风速、风向及其误差的日变化Fig.9 The diurnal variation of wind speed, wind direction and their errors for 10 m wind in May 2015

图10 2015年5月试验时段内70 m各方案模拟风速、风向及其误差的日变化Fig.10 The diurnal variation of wind speed, wind direction and their errors for 70 m wind in May 2015

图11 2015年5月试验时段内模拟和观测的风切变指数Fig.11 The wind shear index of simulations of different BL schemes and observations during the hindcasting period in May 2015