风能作为极具开发潜力的可再生能源，为减少温室气体排放提供重要途径，在电力行业的重要性日渐突显^[1]. 近年来，随着风电场建设规模的不断扩大，风能的利用与开发已从规模化集中开发转向规模化开发与分散式开发并举的新模式^[2]. 风资源丰富地区的风电场的选址、大规模风电场的运行特性及其对大气边界层影响等的评估将成为风能开发与利用过程的重要任务.

传统的风资源评估主要依赖于当地测风塔的观测数据，即基于若干测点的数据进行外推获得区域的风资源分布. 然而，现场实测数据往往较稀缺，且无法覆盖整个场区；利用若干点进行的外推估算也存在较大的误差^[3]. 由此可见，仅仅利用传统测风塔进行的风资源评估存在一定的局限性. 对未建风电场区域风资源的预测，和对已建成风电场周围气候特性的评估，在时空尺度上均属于β天气预报范畴. 因此，采用数值模拟与实测研究相结合的研究方法，发展高精度的大规模风电场区域的风资源评估的天气预报数值模式可为风电产业高效可持续发展提供有力的技术保障.

近些年，基于同化资料开展的数值天气预报（numerical weather prediction, NWP）模式成为风资源评估的热点^[4]. 相比于传统的观测方法，天气预报数值模式可以更灵活地获得任意时空尺度的风资源数据，从而实现对特定区域风资源的有效评估^[5]. 因此，关于利用数值模式开展风资源评估的研究屡见不鲜. Jiménez等^[6]利用天气预报气象模型MM5和风资源评估软件WAsP对北海德国海域的风资源进行评估，发现两者在近海的误差基本相同，MM5气象模式的优势在于无须提供基础的测量数据，但须进一步开发；Lu等^[7]利用MM5模型对我国东部沿海地区的风资源特性进行评估；Lazic等^[8]在对该模型的精确性分析的基础上，对瑞典Gotland风电场的功率输出进行预报. 同时，天气预报模式（weather research and forecasting，WRF）在天气预报风资源评估数值模拟中得到了广泛应用. Zhao等^[9]基于WRF模型开发风电场短期功率预报系统. Wang等^[10-11]分析WRF模式对风电场近地风速模拟的误差特征，指出模式误差与大气稳定性密切相关，同时验证了同化常规资料对低层大气初始场的质量和风的预报性能都有所改善^[12-13]. 此外，国外诸多学者已将WRF模型成功应用于对巴基斯坦^[14]、丹麦^[15]、葡萄牙^[16]等国家陆上及海上风能资源的评估和预测. Fitch等^[17]提出风电场参数化（wind farm parameterization, WFP）模型，利用该模型分析风电场对大气环境的影响规律. 然而，现有关于风电场风资源特性的研究大多忽略了已建风电场与大气边界层之间的相互作用. 因此，在大规模风电场风资源评估的天气预报数值研究中，须进一步考虑大规模风电场与大气边界层之间的相互作用.

以我国张北县百万千瓦风能基地为研究对象，建立风电场参数化模型，开展天气预报（中尺度）数值模拟. 利用风电场附近测风塔的观测数据，分别对考虑风电场参数化模型及未考虑风电场参数化模型的天气预报模式的模拟结果进行验证分析，从而探究耦合风电场参数化模型的天气预报数值模式的有效性，可为风能富集区的风电场选址、风电场运行特性及其对气候影响的分析提供重要的技术支撑.

大规模风电场位于我国张北县百万千瓦风能基地，该风能基地位于河北省西北部、内蒙古高原南缘的坝上地区，地处华北内地连接内蒙古的咽喉地段，地理坐标为40°57'N~41°34'N、114°10'E~115°27'E. 东南部多复杂高地，最高海拔为2 128 m；北、中部地势平坦，向西北渐低，最低点海拔约为1 300 m. 该基地是风电场较佳的选址区域，截止2016年底该基地风电总装机容量约为2.82×10⁶ kW. 此外，该地属于温带半干旱大陆季风气候，具有显著的季度性，大风常出现在冬季和春季.

观测资料来自中部平坦地形风电场及东南复杂山地地形区域风电场附近2台测风塔M1（114.821 5ºE，41.238 6ºN）、M2（115.148 1ºE，41.185 7ºN）的2016年全年的测风数据. 选取1月、4月、7月和10月4个月份的测风塔70 m高度处的10 min平均风速及风向. 该观测数据连续稳定且有代表性，通过了有效值检验和一致性检验等质量控制^[18]，可为模式准确性验证提供可靠的对比数据.

采用新一代数值天气预报模式和数据同化系统相结合的Advanced Research WRF（WRF-ARW, V-3.7.1）开展相关试验研究. 该模式是完全可压的非静力模式^[19]，采用地形追随非静力气压垂直坐标、Arakawa C交错网格、4阶Runge-Kutta显式离散差分格式，满足质量、动量、干熵守恒. 利用σ垂直坐标对可压缩Euler方程进行转化及动力学求解，通量形式为

(1) $\left[ {{{U}},\varOmega ,\varTheta } \right] = {\mu _{\rm d}}\left[ {\left( {u,v,w} \right),\eta ,\theta } \right].$

式中：μ_d为干空气质量，(u, v, w)为水平和垂直协变速度，η为垂直坐标变量，θ为位温.

为了便于构造守恒的差分格式以及讨论大气中的能量转换关系，通常采用湿空气Euler通量方程：

(2) ${\partial _t}U + (\nabla \cdot {{U}}u) + {\mu _{\rm d}}\alpha {\partial _x}p + (\alpha /{\alpha _{\rm d}}){\partial _\eta }p{\partial _x}\varphi = {F_U},$

(3) ${\partial _t}V + (\nabla \cdot {{U}}v) + {\mu _{\rm d}}\alpha {\partial _y}p + (\alpha /{\alpha _{\rm d}}){\partial _\eta }p{\partial _y}\varphi = {F_V},$

(4) ${\partial _t}W + (\nabla \cdot {{U}}w) - g[(\alpha /{\alpha _{\rm d}}){\partial _\eta }p - {\mu _{\rm d}}] = {F_W},$

(5) ${\partial _t}\varTheta + (\nabla \cdot {{U}}\theta ) = {F_\varTheta },$

(6) ${\partial _t}\varphi + \mu _{\rm d}^{ - 1}[({{U}} \cdot \nabla \varphi ) - gW] = 0,$

(7) ${\partial _t}{Q_{ m}} + (\nabla \cdot {{U}}{q_{ m}}) = {F_{{Q_{ m}}}}.$

式中：α=1/ρ=α_d(1+∑q_i)⁻¹为比体积，α_d为干空气的比体积，∑q_i为水汽、云量、雨量等相对干空气的质量混合比总和；F_U、F_V、F_W、F_Θ分别为由物理量、湍流掺混、球面投影和地球自转引起的作用力项；φ=gz为重力势；p为压力；U、V、W分别为水平（经纬）方向速度通量和垂直方向速度通量；Q_m = μ_dq _m，q_m为通用湿空气混合比，即水蒸气、云或雨等相对于干空气的质量混合比； ${F_{{Q_{ m}}}} $为通用湿空气变量所引起的作用力项.

干空气比体积的诊断方程和总压的关系为

(8) ${\partial _\eta }\varphi = - {\alpha _{\rm d}}{\mu _{\rm d}},$

(9) $p = {p_0}{[{R_{\rm d}}{\theta _{\rm m}}/({p_0}{\alpha _{\rm d}})]^\gamma },$

(10) ${\partial _t}{\mu _{\rm d}} + (\nabla \cdot {{U}}) = 0.$

式中：γ=c_p/c_V=1.4为干空气定压比热容和定容比热容之比；θ_m=θ(1+1.16q_v)为湿空气位温，q_v为水汽质量混合比；R_d为干空气气体常数，p₀为标准大气压力. 该模式的优点在于能够为求解大规模风电场周围流动特征及其对周边大气边界层的作用提供真实的天气背景.

采用四维变分同化处理的NCEP-FNL资料作为初始条件，利用MODIS遥感土地覆盖分级的土地用途类别和修改Noah土壤表面模型来替代美国地质勘探局的数据，静态地形数据分辨率为“modis_30s+30s”，采用Lambert投影方式. 模式在水平方向采用三重网格嵌套，其中大规模风电场位于最内层区域，如图1所示. 图中，z为高程，d03中M1、M2为2台测风塔位置、黑色小点为风力机位置. 各嵌套层的水平网格数分别为160×160、106×106、130×130，对应水平分辨率分别为9、3、1 km. 由于场区内存在复杂的山地地形区域，相比以往文献^[20-22]，对垂直方向的网格进行了非均匀加密处理，垂直方向共81层（顶高20 km），将近地层200 m以下高度区域划分为23层网格，风力机风轮旋转区域划分为13层网格. 最内层积分时间步长为3.2 s，每10 min输出一次结果；为了获得稳定的求解，每次求解起始的24 h作为缓冲时长（spin-up time），积分时长为744 h.

为了增加垂直风廓线的精度及减少模型对时间步长的敏感性，采用WRF Single-Moment 6-class（WSM6）方案的微物理过程^[23]；采用RRTM长波辐射方案^[24]及Goddard短波辐射方案^[25]；采用Grell-Freitas积云参数化方案^[26]. 真实大气环境下风电场的运行涉及动量输运、热量输运、水汽输运、剪切效应和地形强迫等物理过程，决定下垫面与大气之间的动量、热量及水汽交换，MYNN-2.5边界层参数化方案能够准确地描述风电场对边界层过程的影响，能保证风电场参数化模型的有效^[27-28].

为了在WRF模式下研究风电场与大气之间的相互作用，Fitch等^[17]提出在平均流场上施加动量源的风电场参数化（WFP）模型. 该模型假设风力机从环境中获得的能量一部分转化为电能，其余转化为湍流动能（turbulence kinetic energy, TKE）. 其中，获取的总能量用风力机的推力系数C_T表征，转化为电能的部分用功率系数C_P表征. 根据致动盘理论^[29]，风力机对大气边界层的作用力表达式为

(11) ${{ F}_{{\rm{drag}}}} = {C_{T}}\rho \left| {{U}} \right|{{U}}A/2 .$

式中：U = [u, v]为来自风轮上游的水平速度矢量；ρ为空气密度，取1.23 kg/m³；A = (π/4)D²为风力机的风轮旋转面积，其中D为风轮直径；C_T取决于来流风速及风力机类型，根据来流风速及风力机推力系数曲线获得.

由于模式最内层的水平网格分辨率为1 km，有可能多台风力机位于同一个网格中. 假设风力机的风轮旋转面与来流风向正交，且忽略同一网格内风力机之间的尾流干扰效应^[17]. 那么，运行过程中的风力机将对当地大气环境产生动量亏损. 在笛卡尔坐标系下，网格(i, j, k)内由风力机的扰动而引发的动能损失的变化率为

(12) ${{\partial {\rm{KE}}_{{\rm{drag}}}^{ijk}}}/{{\partial t}} = - N_{\rm T}^{}\Delta x\Delta y{C_{ T}}{\rho _{ijk}}\left| {{U}} \right|_{ijk}^3{A_{ijk}}/2.$

式中：N_T为单位网格内风力机的台数；Δx、Δy为网格在经纬度方向的长度；|U|_ijk为网格(i, j, k)处速度的模；ρ_ijk为网格（i, j, k）处空气的密度；A_ijk为网格(i, j)处垂直方向k与k+1层之间的风轮旋转面截面面积. 相比较小的垂直速度，水平方向的速度分量为主导因素，故网格内流体的动能损失的变化率为

(13) ${{\partial {\rm{KE}}_{{\rm{cell}}}^{ijk}}}/{{\partial t}} = {\rho _{ijk}}{\left| {{U}} \right|_{ijk}}\left( {{z_{k + 1}} - {z_k}} \right)\Delta x\Delta y {{\partial {{\left| {{U}} \right|}_{ijk}}}}/{{\partial t}}.$

式中：z_k、z_k+1分别为模型垂直方向k、k+1层的高度. 根据能量守恒原理，某网格内流场能量的变化率和由风力机旋转造成的能量亏损速率相等，即式（12）与式（13）相等，则可得动量源项为

(14) $\frac{{\partial {{\left| {{U}} \right|}_{ijk}}}}{{\partial t}} = - \frac{{N_{\rm T}^{}{C_{ T}}\left| {{U}} \right|_{ijk}^2{A_{ijk}}}}{2{\left( {{z_{k + 1}} - {z_k}} \right)}}.$

因此，对包含运行工况风力机的网格（i, j, k），风轮截面A_ijk对应部分从大气中获得能量并转化为有用的输出功率P及耗散在大气中的湍流动能TKE的表达式分别为

(15) $\frac{{\partial {P_{ijk}}}}{{\partial t}} = - \frac{{N_{T}^{}{C_{ P}}\left| {{U}} \right|_{ijk}^3{A_{ijk}}}}{2{\left( {{z_{k + 1}} - {z_k}} \right)}},$

(16) $\frac{{\partial {\rm{TK}}{{\rm{E}}_{ijk}}}}{{\partial t}} = - \frac{{N_{\rm T}^{}{C_{\rm {TKE}}}\left| {{U}} \right|_{ijk}^3{A_{ijk}}}}{2{\left( {{z_{k + 1}} - {z_k}} \right)}}.$

式中：C_P取决于来流风速及风力机的类型，可依据来流风速及风力机功率输出曲线求得；模型假设耗散在大气中湍流动能的表征系数C_TKE=C_T−C_P. 基于此，WFP模型可通过在平均流场上施加动量源及湍流动能源来体现风力机对大气的扰动；由于扰动增加的风剪切效应将进一步促进大气在垂直方向的掺混，该效应将由WRF模式中的边界层方案考虑^[17].

根据目标风电场资料，包括大规模风电场内1 880台风力机的地理坐标及总装机容量，选定额定功率为1.5 MW的主流风电机组进行参数化建模. 其中风力机的风轮直径为89 m，额定风速、切入和切出风速分别为10、3、25 m/s，风力机的基本参数、推力系数曲线及功率输出曲线，如表1和图2所示.

大规模风电场会对其内部及周围大气边界层产生影响，从而影响局地或区域风资源特征^[30]. 为了检验WFP模型对大规模风电场风资源特性的评估能力，选取和观测资料一致的模拟周期，设计并开展2组数值试验：1）仅采用WRF模式、未耦合WFP模型，记为WRF-NWFP；2）耦合WFP模型的WRF模式，记为WRF-WFP. 利用评估指标对风速、风向的模拟效果进行定量统计和对比分析.

为了定量地评价模式对场区风资源评估的准确性，分别计算模拟结果与观测资料的平均偏差（mean bias error, MBE）、均方根误差（root mean square error, RMSE）、相对均方根误差（relative root mean square error, RRMSE）和一致度（index of agreement, IA）等评价指标.

1）平均偏差表征模拟结果与观测结果的总体偏差情况，表达式为

(17) ${\rm MBE} = \frac{1}{n}\sum\limits_{i=1}^n {\left( {{S_i} - {O_i}} \right)} .$

式中：S_i、O_i分别为第i个模拟和观测的统计样本.

2）均方根误差表征模型的精确度，通常为正值，表达式为

(18) ${\rm RMSE} = {\left({\frac{1}{n}\sum\limits_{i=1}^n {{{\left( {{S_i} - {O_i}} \right)}^2}} } \right)}^{1/2}.$

3）相对均方根误差由均方根误差除以观测结果的平均值得到：

(19) ${\rm RRMSE} = { {\left({ {\frac{1}{n}\sum\limits_{i=1}^n {{{\left( {{S_i} - {O_i}} \right)}^2}} } } \right)^{1/2}} \Bigg/ \left({ \sum\limits_{i=1}^n {{O_i}} }/n \right)}. $

式中：RRMSE < 0.1表示模型精度为优，0.1 < RRMSE < 0.2表示模型精度为良，0.2 < RRMSE < 0.3表示模型精度为中，RRMSE > 0.3表示模型精度较差 ^[9].

4）一致度通常表征模拟结果对观测结果预测的精确程度：

(20) ${\rm IA} = 1 - \frac{{\displaystyle\sum\nolimits_{i = 1}^n {\left| {{S_i} - {O_i}} \right|} }}{{\displaystyle\sum\nolimits_{i = 1}^n {\left( {\left| {{S_i} - {S_{\rm avg}}} \right| + \left| {{O_i} - {O_{\rm avg}}} \right|} \right)} }}.$

式中：S_avg、O_avg分别为模拟和观测样本的平均值. 本试验分别统计了模拟结果和观测结果的10 min平均样本，样本总数n = 4 320.

以张北县百万千瓦风能基地大规模风电场为研究对象，选取2016年4月、7月、10月、1月作为春季（SPR）、夏季（SUM）、秋季（AUT）、冬季（WIN）4个季度的代表月份，分别采用WRF-NWFP模式和WRF-WFP模式进行2组天气预报数值试验. 分析不同季度下大规模风电场平坦地形区域及复杂山地区域风速的时域分布及其威布尔概率密度函数分布和风向特征，与观测结果进行对比分析，从而评估WRF-WFP模式对大规模风场区域不同季度及地形条件下风资源特性的模拟能力.

如图3所示为由WRF-WFP和WRF-NWFP模式获得的各季度风电场内平坦和山地区域2台测风塔（平坦地形为M1，山地地形为M2）70 m高度处的模拟风速及观测风速的时间序列分布曲线（灰色区域为观测风速的偏差范围）以及相应时间序列下WRF-NWFP与WRF-WFP模式模拟风速差值ΔWS的柱状图. 由图3可知，不论在平坦地形还是在复杂地形，模拟结果均与观测结果较吻合，呈现出显著的昼夜交替规律；各季度风速的模拟值偏高. 这主要由于模式采用的陆面数据与实际地形存在偏差，低估了亚格子尺度地形引发的阻力效应^[29]. 对比各季度ΔWS可以看出，虽然WRF-NWFP和WRF-WFP模式对风速均存在高估现象，但后者的高估幅度较低，尤其在平坦地形下的夏、秋季节较明显. 可见，考虑风电场与大气边界层相互作用的WRF-WFP模式对风速的模拟精度更高.

如表2所示为WRF-WFP和WRF-NWFP模式各季度平坦地形和复杂地形下模拟与观测风速对比的准确性评价指标. 可以看出，模拟结果与观测结果的平均偏差为1.55~2.88 m/s，均方根误差为1.99~3.65 m/s；夏季以外的其他季度的相对均方根误差均小于0.2；一致度均接近或大于0.7. 模式对夏季风速的模拟效果略弱，主要原因在于大规模场区的气候特征属于温带大陆性季风气候，在夏季易受到来自东南方向的海洋热带气旋的影响，从而增加了模式对风速预测的不确定性；夏季风速较低也是预测精度较低的原因之一. 进一步对比WRF-NWFP和WRF-WFP模拟风速与测风塔观测风速，可以看出，两者均具有较好的模拟精度，其中后者表现出了一定的优势：相对WRF-NWFP模式，在平坦地形区域，WRF-WFP模式的IA提高了3.3%（SPR）、3.6%（SUM）、4.6%（AUT）、2.8%（WIN）；在复杂地形区域，提高了3.2%（SPR）、3.1%（SUM）、4.3%（AUT）、2.6%（WIN）. 总体而言，在平坦地形下WRF-WFP模式对风速的模拟能力提高了约3.6%，而在复杂地形下提高了约3.3%.

此外，该模式对风速的模拟能力表现出季节差异性. 例如，在夏、秋两季场区的主风向为南，M1位于区域风电场的下游，易受到风电场尾流效应的影响，利用耦合风电场参数化模型的WRF-WFP模式对区域风资源评估十分必要. 然而，该模型对风速的模拟结果仍偏高，除了上文所述的模式对地形解析存在偏差或物理参数化方案带来的系统误差等主要原因^[31]，还由于WFP模型假设忽略了同一网格内部风力机之间的尾流干涉效应^[17]，低估了风电场对周边大气边界层的尾流干涉作用.

除了对风速的时间序列进行评估外，利用风速的概率密度函数判断模式在各风速区间的准确性. 采用威布尔概率密度分布函数，对2个站点70 m高处度的模拟及观测风速进行统计分析. 经典威布尔（Weibull）概率密度函数（probability density function，PDF）表达式为

(21) $f\left( v \right) = {k}/{c}{\left( {{v}/{c}} \right)^{k - 1}}\exp \;\left[ { - {{\left( {{v}/{c}} \right)}^k}} \right].$

式中：v为风速，k为形状参数，c为尺度参数. k越大，风速频率分布越向平均风速集中，即较小或较大风速发生的概率较低. 如图4所示为由WRF-WFP和WRF-NWFP模式获得的各季度M1和M2站点的风速直方图及其概率密度函数曲线. 其中，f为频率；k_O、k_S1、k_S2和c_O、c_S1、c_S2分别为观测结果、WRF-WFP模式和WRF-NWFP模式的形状参数和尺度参数；图4(a)~(d)为平坦地形、图4(e)~(f)为复杂地形.

由图4可以看出，不同季度下WRF模式（包括WRF-WFP和WRF-NWFP）与站点观测获得的风速概率密度分布之间呈现明显的季度差异性. 例如，春季、秋季和冬季模拟的概率密度分布与观测结果较吻合（如图4(a)、(c)、(d)、(e)、(g)、(h)），而夏季两者偏差较大（如图4(b)、(f)）. 同时，不同季度下均存在“高估高风速区间的风速概率密度、低估低风速区间的风速概率密度”的特点. 主要表现为：冬季的模拟效果最佳，M1站点模拟与观测结果的误差最小，且相应的PDF形状参数及尺度参数高度一致；M2站点的模拟结果略微高估了6~9 m/s区间的风速；春、秋两季的模拟结果总体上低估了10 m/s以下风速区间的风速而高估了10 m/s以上风速区间的风速；夏季模拟结果较大幅度地低估了低风速区而高估了高风速区. 类似结果在文献[32]中也有报道. 进一步比较各季度下WRF-NWFP和WRF-WFP的模拟风速概率密度函数可以发现，WRF-WFP获得的风速PDF分布特征更接近观测风速的PDF，特别是对M1站点秋季风速模拟的改善效果较明显. 主要原因是WRF-WFP模型充分考虑了风电场对大气边界层的影响，更加真实地反映了风电场周围的风资源特征.

因此，相比WRF-NWFP模式，WRF-WFP模式对中低风速区间概率密度的低估幅度以及对高风速区间概率密度的高估幅度均有所降低，能够较准确地预测风速的概率密度分布，从而较好地评估当地风资源的分布情况.

风向是风电场风资源评估过程须考虑的重要因素. 分别提取WRF模式（WRF-WFP和WRF-NWFP）模拟的风向与观测站点实测的风向数据，绘制风向玫瑰图，如图5所示.

可以看出，不论是在平坦地形区域还是在复杂山地地形区域，WRF-NWFP和WRF-WFP模式均能够准确地预测场区各季度的主导风向，即模拟得到的主导风向与观测到的主导风向保持高度的一致性. 此外，在大规模风电场内平坦地区和山地区域的季度风向玫瑰图分布基本一致，春、夏、秋、冬4个季节的主导风向分别为西北（NW）、南（S）、南（S）、北北西（NNW）. 如表3所示，定量地统计了各季度WRF-NWFP和WRF-WFP模式对平坦地形和复杂地形区域风向模拟的评价指标. 结果表明，两者的平均偏差为20.76°~29.83°，均方根误差为27.06°~43.34°，相对均方根误差均小于0.1，一致度高于0.7，表明两者对风向的评估均达到较好的精度. 比较WRF-NWFP和WRF-WFP模式可以发现，后者对风向的模拟的精度比前者高约3%.

（1）WRF-NWFP模式和WRF-WFP模式的模拟结果与观测结果较吻合. 后者有一定的优势：在平坦地形区域、复杂地形区域，其模拟精度增幅分别为4.6%、4.3%. 不论对风资源评估的准确性还是对风电场功率预测的精确性均有所提升.

（2）模拟结果存在高估高风速区间的概率密度和低估低风速区间的概率密度的特点，但WRF-WFP模式对中低风速区间概率密度的低估幅度以及对高风速区间概率密度的高估幅度均有所降低，能够较准确地预测风速的概率密度分布，从而较好地评估当地风资源的分布情况.

（3）WRF-WFP模式模拟的主导风向与观测的主风向具有高度的一致性，能够对不同季度大规模风电场内平坦地形和复杂地形区域的主导风向进行准确的预测.

（4）耦合风电场参数化模型的天气预报数值模式能够提高对大规模风电场风资源特性评估的精度.

（5）发展精度较高的天气预报数值模式对我国风能利用与开发技术创新具有重要意义，可为风资源丰富地区风电场的微观选址、大规模风电场的运行特性及其对大气边界层影响的评估提供重要的技术支撑，从而为我国风电产业高效可持续发展提供技术保障.