山区风能资源的模拟研究[1]

刘永强　朱瑞兆

（国家气象局气象科学研究院）

风能指大气运动所具有的动能，通常用风能密度w表示。在山区，由于不均匀下垫面的增幅作用，常常形成局地强风能区。

然而，由于山区测站稀疏及测值代表性差等原因，不能像全国及区域范围风能的研究那样，直接由实测资料计算风能密度［1］。作为气象学研究的一个基本问题，山区气流分布的理论及场地试验方面的研究较多，不过尚未有较成熟的定量计算方法。近年来，山地气流的中尺度数值模拟研究发展较快，有些模拟结果可直接用于风能的估价［2］，但通常使用的模式包含的物理过程复杂，要求的网格较细，因而计算量较大，故目前多用于理论研究。

为了提供一种既能定量计算，又较为实用的估价山区风能资源的方法，本文根据Sherman研究大气扩散背景流场的基本思想［3］和地表气流基本平行于地表面的事实，采用地形坐标系中的水平无辐散关系作为约束条件，在变分的意义下对地表面风场进行调整。这样，可避免文献［3］中由于地表作为固定下边界而致使地表流场模拟的精度不高之不足（这对风能估价的精度有很大的影响）。此外，为了增加模式的实用性，模拟计算中还采用了自然正交展开的方法，可大大减少无辐散调整的次数。

一、模拟计算区的地形特征和资料

本文选择内蒙古自治区乌兰察布盟南部山丘地带作为模拟计算的区域（以下简称“计算区”）。该区域位于大青山的东部，由几个相对高度为数百米的小山组成，可近似地看作孤立的山地（见图1）。

用于模拟计算的资料主要是地面网格点高度和地面风。前者从内蒙古自治区分区地图（比例尺1∶20万，等高距40m）读数获得。地面风资料包括计算区内部及周围六个地面气象台站［见图1（a）］1980年5月的实测风速。观测间隔为1h，故每站资料容量为744。

二、模式及计算方案

（一）无辐散调整过程

在测站稀疏的山区，客观分析得到的流场往往是不太可靠的。通常所采用的补救方法之一是对客观分析的结果［以下称“初始场”，记为V0=（u0，v0）］在某种物理约束下调整到一个新的流场［记为V=（u，v）］，这样似乎能更客观地反映实际情况。

观测表明［4］，气流过山时，可分为绕流和越山气流两支；在风速不是太强时，除在背风面及迎风面山脚处有垂直涡旋产生外，越山气流通常沿地面呈水平运动。这种特征对于孤立山地更为明显。在这种情况下，可近似认为气流满足地形坐标系中水平两维无辐散关系（这相当于连续方程取零级近似）。这一关系常用于山地气流理论的研究［5］。本文将以此作为风场调整的约束关系。

引入地形坐标

这里，H（x，y）和h（x，y）分别为行星边界层高度和地面高度；z为局地直角坐标系中垂直方向的坐标。在地形起伏不是很大的情况下，可近似地取地面的空气密度为常数。于是，水平无辐散关系可写为

设，这里Hd=H-h。

式（2）即为风场调整的约束条件。除此以外，还要求调整幅度在变分的意义下达到最小。相应的数学问题为：寻找一个新的流场V*，使

达到最小。式中σ为调整的区域，λ=λ（x，y）为拉格朗日乘数。

根据变分原理，有

在边界上，取

式（3a）和（3b）即为无辐散调整的基本方程。

（二）自然正交函数的应用

以上调整过程是对某一时刻的风场而言的。在风能评价中，确定某地风能潜力所需的资料序列长度通常至少为一年。设序列容量为M，则对于3h间隔的观测，M至少为103的量级。这样，上述调整方程就须重复求解数千次；显然，从计算时间的角度来说，是相当花费的。然而，根据调整方程（3a）的性质，将风场进行自然正交展开，则可大大减少调整的次数［6］。

方程（3a）可写为

这里，F为相应方程的右端项。这是一个线性偏微分方程，其解具有线性可加性。即若

则对任意常数Ci，有

下面讨论自然正交函数在减少调整次数方面的作用。设第k个测站l时刻的风速为Ukl，Vkl（k=1，2，…，n；l=1，2，…，m。n，m分别为测站数和测值序列容量），该站平均风速为。Ukl，Vkl可表示为自然正交函数ukj，vkj（j=1，2，…，N；N=2n）的线性组合：

这里，当k为奇数时，第一式中以取代。

图2描述了自然正交函数应用于无辐散调整的过程。其中，“模式输入i模式输出i”（i=Ⅰ，Ⅱ，…，N+1）表示对第i个特征向量（即自然正交函数）或平均风向量进行无辐散调整，而下标s表示计算区格点标号。整个过程包括分解、调整和合成三部分。由此获得的调整风场应与直接对实测风进行调整获得的结果一致。

自然正交展开实现了风场的时、空分离，而无辐散调整只需对仅随空间变化的N个自然正交函数和平均风向量进行就行了。这样，大大减少了调整次数，因而减少了计算量，使得计算模式更为实用。此外，自然正交函数由整个资料序列的协方差矩阵决定，因而包含了历史资料的信息。从某种意义上来说，这部分地弥补了山区测站稀疏之不足。

在对各正交函数和平均风进行了无辐散调整以后，某一时刻调整风场的获得就可简单地归结为：用该时刻的实测风由式（6）第二式计算展开系数aj，再由式（6）第一式便可得到所需结果。

（三）计算方案

1.计算区网格和初始场。

计算区分为16×16的正方形网格，网格距为16km。

初始风场由各站的实测风加权平均获得

式中：r——格点S与第k个测站之间的距离；

Wk（r）——该站在S处的权重系数，取为［7］

其中，a为大于0的参数。这里忽略了文献［7］中考虑的另一个因子——风向的影响。

2.边界层高度。

由于探空站很少且探测间隔较大，加上计算区的地形复杂，因而很难得到边界层顶的时空结构。本文采用经验公式［8］：

式中：hs——计算区地形最高处的高度；

Ha——边界层平均厚度；

K——边界层顶斜率。

Ha和K的取值见表1。这种取法只能反映边界层变化的一般特征，即夜晚低平，而白天与之相反。

3.调整方程的求解。

采用差分近似将调整方程（3）变为一组线性代数方程，并以超张弛迭代法求解。张弛系数取为1.8。

三、计算结果

（一）自然正交展开的一些特征

表2给出了自然正交展开的一些计算结果。若以前N0（＜N）个正交函数（以相应的特征值λj的数值从大至小依次排列）逼近原风场，则相对误差R2可由下式估计。

其中，表示λj对应的特征向量（即正交函数）在风场中所占的比重。N0=3时，σ2=78.3%或R2=21.7%，可见收敛速度是较快的。表中还给出这三个特征值对应的特征向量（正交函数），由于空间测站太少，故不易反映具有明显意义的流场分布特征。

图3为展开系数序列a1～a4在19—22日的变化。一个明显的特征是每一序列（尤其是a1和a2正负相间）呈准周期变化，周期近似为24h，这与实测风的变化特征（图略）是一致的。图3的另一个特征是每个序列的变化依一定的秩序进行。a1达到极大时（第19日02时），a2开始增加；a2达到极大时（约19日10时），a3又开始增加；……正负中心依次传递。前三个序列之间的关系似更明显一点。这一特征表明了不同序列之间的某种联系，它对序列的预报有一定的意义［9］。

（二）无辐散调整风场

我们以平均风（表2最右端一列向量）为例，讨论客观分析和无辐散调整的计算结果，并用中旗站的实测资料进行验证。计算中取Ha=800m，K=0.5，hs=2000m，a=0.1km-2。由式（7）求初始场时，先后用5站（不包括中旗站）和4站（再去掉卓资站）的风测值进行内插和调整，其结果见图4。

由图可见，无辐散调整后，中部出现风速负变区，其两侧为正变区。中旗站位于正变区内。与图1比较，可以发现正（负）变值区基本上对应于地形较高（低）处。调整结果与参数的取值有一定的关系，Ha或K增大时，调整幅度减小。

对平均风和各正交函数进行调整后，以展开系数为权重作线性组合，便可得调整风场。设V′为调整风V与实测风V0的差值，它们的均值分别为和，并定义：

分别表示相对误差和相对方差。表3给出了中旗站的计算结果。表中与同义，只是仅对Vi0≥8m/s的风速（共158个）进行计算。表中给出三次计算结果，其中Ⅰ、Ⅲ分别为使用4站和5站实测风进行内插和调整所得的结果，Ⅱ除未进行调整外，其余与Ⅰ类似。对于Ⅰ，还计算了月平均风速日变化（见图5）。

由表3和图5可见：①与Ⅱ相比，Ⅰ中Ra较小，这表明无辐散调整对风速估计的精度有一定改进。②Ⅰ的结果较Ⅲ好。Ⅲ中多使用了地处计算区谷地的卓资站的资料。由此可见，由于局地地形的影响，该站风资料的代表性较差，似不能较好地表示计算区其他区域流场的特征。③计算风速相对于实测风有一负的系统性误差，其值介于-1.0～-2.0。此外，比较而言，白天（风速较大）误差较小。

四、讨论

以上叙述了应用无辐散调整进行山区风能资源模拟计算的过程和实例。尽管计算模式包含的物理过程简单，计算网格较粗，且计算结果依赖于一些经验参数，但由于具有计算量小、易于修改、且根据通常可以获得的地面资料就可计算等优点，因而对于实际计算有一定的意义。尤其是引入自然正交展开后，调整次数大大减少，从而增加了模式的应用性。计算结果表明，应用这种方法可以得到一定精度的山区风能估计值。

最后，我们对下面有关问题作一些说明和讨论。

1.由于资料缺乏，本文只对位于中旗站的网格点处的计算结果进行了验证，因而不能完全说明整个计算区（尤其是调整幅度最大区）模拟结果的可靠程度。

2.当过山气流很强时，越山气流易在背风面出现波动及涡旋，同时绕流也由于较强的侧向切变而导致涡旋的产生。这时若用本文采用的水平无辐散关系作为风场调整的约束关系，显然会引起较大误差。

3.在自然正交展开中，正交函数只是空间的函数这一性质仅仅是对确定的资料序列而言的；随着序列的改变，正交函数也会改变。为使正交函数具有一般性（或相对稳定性），可适当增加资料序列的容量。

4.在调整方程的右端顶中，；由于本文规定Hd只取几种状态，因而可认为只有V0随时间变化，对其展开，便可减少调整次数。在实际大气中，Hd也随时间而变，这时若仅对V0展开就不再能够达到减少调整次数的目的。在这种情况下，可对Hd（或直接对）进行自然正交展开。

参考文献

［1］朱瑞兆，薛桁.风能的计算和我国风能的分布［J］.气象，1981（8）：26-28.

［2］Segal M，Mahrer Y，Pielke R A.Numerical model study of wind energy characteristics over heterogeneous terrain-Central Israel case study［J］.Boundary-Layer Meteor.，1982，22（3）：373-392.

［3］Sherman C A.A mass-consistent model for wind fields over complex terrain［J］.J.Appl.Meteor，1978，17（3）：312-319.

［4］傅抱璞.山地气候［M］.北京：科学出版社，1983.

［5］杨大升.动力气象学（修订本）［M］.北京：气象出版社，1983.

［6］Ludwig F L，Byrd G.An efficient method for deriving mass-consistent flow fields from wind observations in rough terrain［J］.Atmos.Enviro.，1980，14（5）：585-587.

［7］Endlich R M，Mancuso R L.Objective analysis of environmental conditions associated with severe thunderstorms and tornado［J］.Mon.Wea.Rev，1968，96（6）：342-350.

［8］Endlich R M.Wind energy estimates by use of a diagnostic model［J］.Boundary-Layer Meteor.，1984，30（1-4）：375-386.

［9］Paegle J N，Haslam R B.Empirical orthogonal function estimates of local predictability［J］.J.Appl.Meteor.，1984，21（2）：117-126.

---

[1]本文发表在《气象学报》，1988年第1期，收录在《风能、太阳能资源研究论文集》，气象出版社2008年版。