2013, 24 (1): 109-116
2. 中国科学院研究生院，北京 100049;
3. 中国气象局风能太阳能资源评估中心，北京 100081
2. Graduate University of Chinese Academy of Sciences, Beijing 100049;
3. Center for Wind and Solar Energy Resources Assessment, CMA, Beijing 100081
风力发电厂的选址和功率输出推测需要一段较长时间的风记录数据，一般而言，测风塔只用作短期试验，满足不了评估所需的风记录信息。因此，将附近台站的数据与测风塔短期同步的风记录建立起良好关系，利用台站长期的风记录推测选址地区的风速、风向非常实用且有效。这种方法叫做测量相关推测法 (measure-correlate-predict，简称MCP方法)，是反映规划风电场风能长期平均水平的主要技术手段，也称为数据订正法。
MCP方法起源于20世纪90年代，最早由Derrick[1-2]将两站点间的风速建立起一阶线性回归关系，在研究中Derrick利用8个月的数据建立关系，推测出目标站两年的风速序列，随后他又提出幂指数的关系式。Nielsen等[3]提出一元二阶回归、二元回归的关系式。Riedel等[4]在二阶回归式中加入参数来降低二次方系数的量级。Draper等[5]利用同步数据的平均值和方差比建立关系 (方差比MCP方法)。Woods等[6]提出对样本数据进行矩阵分组计算，Mortimer[7]则提出一种将参照站和目标站的风记录按照风速大小、风向角度分组计算的方法，他们在如何分组上进行了大量研究。
以上这些方法都只适用于同一海拔高度上的两个站点。因此已有研究都没有考虑由于海拔高度差异带来的误差。然而现实中，同高度站点很少存在，这是由于地势上存在海拔高度差异；大多数参照站自身高度只有10 m[8]，而风电场风机真实高度至少有70 m[9]；目标站短期测风塔高度达不到所需的高度；需要得到目标站在高度上加密的长期风速分布情况。这些情况下直接使用MCP方法推测不同海拔高度站点风速将引起偏差。
本文选取内蒙古锡林郭勒地区风能测风塔数据进行研究，这是因为内蒙古草原地区地势平坦、下垫面均匀，能够排除地形、下垫面等因素带来的干扰。利用塔上6层不同高度的风杯数据，在方差比MCP方法基础上，利用风速的Weibull双参数分布规律建立起一种能纠正海拔高度所带来影响的修正的MCP方法，它能够利用一个与目标站不同海拔高度的参照站风记录推测出目标站的长期风速情况。
1 数据本文数据来源于锡林浩特草原区域内架设的两座风能资源观测塔 (表 1)，塔间距离约为7 km，该塔属于国家气候中心主持的风能资源详查和评价项目，共有约1年 (2009年5月21日—2010年4月30日) 的风速观测资料。原始数据分别为安装在塔上4 m，10 m，30 m，50 m，70 m，100 m高度处支架上、采样频率为1 Hz的风杯数据。根据Derrick的建议[1-2]，本文选用前8个月 (2009年5月21日—2010年1月21日) 数据进行10 min平均和6点滑动 (1 h) 后再用于推测以减小误差，余下数据 (2010年1月22日—2010年4月30日) 则进行结果推测并对比检验。
|
表 1 两座测风塔信息 Table 1 The geographic information of two anemometer towers |
2 MCP方法简介
MCP方法基于风场的空间相关原理，当目标站和气候站之间的风记录 (风速、风向) 具有比较好的相关关系时 (通过显著性检验)，就可以直接利用参照站数据对目标站长期风速、风向估算，并运用气候学和统计学处理方法进行确定[1-2]。气候学处理方法建立在一个基本假设之上，即风速具有一定程度的空间相关，并且依赖于参照站的观测。图 1为MCP方法的整体流程示意图。将参照站与目标站的同步风向风速数据建立MCP方法的模式，再依照参照站长期数据通过建立好的模型推测风速。
|
|
| 图 1. MCP方法流程图 Fig 1. Flowchart of MCP method | |
MCP方法旨在利用某地区某一风能测量塔来评估附近几个目标点年平均风能蕴藏量，从而确定哪一个点更适合建立风力发电厂；利用参照站的长期数据来订正目标站由于仪器受损、太阳能电池供电不足等因素造成的较长时间的数据丢失。这样可以节省用于风能选址大量的观测时间和试验费用，提高准确度和效率[10]。
20世纪90年代开始，大部分模型都基于线性回归模型进行修正[1-6]：
|
(1) |
式 (1) 中，ŷ为待推测的目标站长期风速；x为参照站的长期风记录；a，b均为待定系数。
本文采用的是MCP方法中具有代表性的、计算简洁、结果准确的方差比MCP方法。该方法是由Draper [5]将风记录按照风速分组 (＜3 m/s，3～4 m/s，4～5 m/s，……，19～20 m/s，＞20 m/s)，每组数据都直接利用参照站和目标站风速的平均值和标准差创建模型，即
|
(2) |
式 (2) 中，μx，μy分别是参照站和目标站风速的平均值，单位：m·s-1；σx，σy则是参照站和目标站风速的标准差，单位：m·s-1。
3 建立修正的MCP方法修正的MCP方法是在方差比MCP方法基础之上，利用风速的Weibull双参数分布建立起的一种方法，它能够通过一个与目标站有着不同经纬度、不同海拔高度的参照站风记录，推测出目标站的长期风速情况，从而消除海拔高度带来的误差。在建立修正的MCP方法之前，先要了解风速的Weibull双参数分布规律。
3.1 风速的Weibull双参数分布规律在过去30年，许多人探讨了风速分布的统计模型，包括Weibull双参数分布模型[11-14]，Weibull双参数分布是拟合实际风速最好的一种模型，允许通过调整参数来适应不同时间长度的风速序列。
Weibull双参数分布仅由k，c双参数来决定，k是形状参数，c是尺度参数。Weibull双参数分布概率密度函数为
|
(3) |
分布函数为
|
(4) |
其中，v为瞬时风速，单位：m·s-1。
形状参数k决定了密度函数曲线的形状，当k=1时，Weibull双参数分布呈指数分布；当0＜k＜1时，图像为以y轴和x轴为渐进的曲线；当k＞1时，峰值随k值增大而增大，图像越来越窄。尺度参数c决定着分布的尺度大小，即当k固定时，c增加图像峰值将降低，曲线变得平缓。Weibull双参数分布由这两个参数决定，起重要作用。只要确定了k，c双参数的大小就能得到风速的时间序列。
因此，只需要抓住双参数随高度分布的特点，得到反映高低空双参数规律相关的方程，就能由低 (高) 海拔高度的双参数去推测高 (低) 海拔高度的双参数，最后将推测出的Weibull双参数分布还原出待推测的风速序列[15]。
本文选取前8个月塔2的6层不同高度风记录，利用最大似然法来作估计，将数据按照风速大小每1 m·s-1为单位分类 (小于1 m·s-1的风速忽略不计，每组不少于1个数据)[16-17]。计算Weibull分布的双参数k和c值，其统计情况如表 2所示。
|
|
表 2 2009年5—12月Weibull双参数k，c统计表 Table 2 The estimated value of c and k for wind speed of six layers from May to November in 2009 |
图 2表示Weibull双参数分布形状参数k随高度的变化情况。从图 2可以看到，在40 m以下，k值随着高度的增加而增大；在50 m高度左右曲线开始出现拐点，k值达到最大值；60 m以上随着高度的增加k值反而减小。参数k与高度z之间可以拟合出这样的关系式：
|
(5) |
|
|
| 图 2. Weibull双参数分布参数k随高度变化曲线 Fig 2. Curve of parameter k of Weibull distribution varying with height | |
式 (5) 中，a为0.0116，b为0.1313，d为2.2186。
由图 3可知，尺度参数c随着高度变化也呈现出特别的规律。拟合关系式为
|
(6) |
|
|
| 图 3. Weibull双参数分布参数c随高度变化曲线 Fig 3. Curve of parameter c of Weibull distribution varying with height | |
式 (6) 中，z为高度，单位：m；α为待估计的指数函数参数。利用最小二乘法来拟合α的值，拟合结果如表 3所示。全年α平均值为0.135，因此将0.135作为指数函数参数的值。这与Tony等[18]的结果一致，说明这样的拟合函数具有普适性，在各种地形地区都适用。
|
|
表 3 α拟合值统计表 Table 3 The value of α |
得到了Weibull双参数分布中双参数随高度变化的规律后，还需要将推测到的Weibull双参数分布还原为待推测的风速时间序列。
在理想条件下，边界层内风随高度呈单调递增。这样不同高度上的风速大小虽然不同，但应该有着相同的等级排列顺序。设Y为高海拔高度风速，X为低海拔高度的风速。将两个高度N个同步数据按由小到大的顺序排列，最小的风速记为r0，最大的风速记为rN-1。若高低空风速呈单调关系，则对应的两个站点同一时刻风速各自的r应相等，即rXn=rYn，0＜n＜N。这样可以认为高低空风速存在着相等的累积分布函数，即rX/(N-1)=rY/(N-1)，这种等式关系可以作为建立风速序列还原的理论基础。
有了这些假设基础，风速的累积分布函数就可以用来代替表达风速的等级排列顺序，这种近似处理合情合理[18]。而风速累积分布遵循着Weibull双参数分布，Weibull双参数分布由k，c双参数决定，只要k，c两个参数已知，就能得到风速时间序列。
3.2 利用Weibull双参数分布修正方差比MCP方法得到Weibull双参数分布中双参数随高度的变化规律后，即可将这种规律代入方差比MCP方法得到一种修正的MCP方法。为了更直观地说明如何建立修正的MCP方法并更好地分析该方法的结果，这里设计了A和B两个对比试验方案 (图 4)。方案A (方差比MCP方法)：直接利用MCP方法中的方差比方法[5-6]，由塔1的10 m高风记录推测塔2的100 m高风速。方案B (修正的MCP方法)：利用方差比MCP方法由塔1的10 m高风杯数据推测塔2的10 m高风速数据；再由塔2的6层风杯数据的风速时间序列拟合出Weibull分布双参数随高度变化的统计公式；最后由Weibull双参数分布建立塔2的10 m高风速与塔2的100 m风速序列的非线性关系，结合方差比MCP方法推测出100 m风速，即一种修正的MCP方法。
|
|
| 图 4. 对比试验方案模拟图 Fig 4. Comparative flowchart of two methods | |
修正的MCP方法不同于方差比MCP方法，关键在于它将Weibull双参数分布的这种规律引入，从而能够修正海拔高度带来的误差。为了在方案A和方案B的对比试验中更方便地说明，这里统一符号[19]。将两个方案参照站的低高度 (参照站10 m处) 风速记作x；最终待推测目标站 (高海拔) 风速都记作ŷ。不同之处，方案B中，将高海拔高度 (目标站100 m处) 的Weibull双参数记为
根据Weibull双参数分布规律，高低空之间的风速关系就可以依靠目标站的10 m和100 m的Weibull双参数分布的4个参数建立关系：
|
(7) |
可得
|
(8) |
将式 (3) 和式 (4) 带入式 (8)，并且结合式 (2)，可以得到修正的MCP方法表达式：
|
(9) |
其中，k10，c10是参照站 (10 m处) 风速Weibull双参数分布的形状参数和尺度参数。由于k，c参数随高度变化特点在很大范围内具有普适性，因此只要得出某地区的式 (9) 就能长期广泛地适用于该地区。
4 推测结果首先对方案A，B中参与建模的2009年5月21日—2010年1月21日数据进行显著性检验，只有通过检验的数据才能进行MCP方法推测[1-7]。
由图 5可知，参照站塔1的10 m风速与目标站塔2的100 m风速的相关系数为0.6174，参照站塔1的10 m风速与目标站塔2的10 m风速的相关系数为0.8821。根据α=0.05的t分布显著性检验可知，39600个样本数据 (8个月) 的相关系数为0.0098，可以认为方案A与B的两个站点相关均显著。但很明显，从图 5中可以看出，方案A中参照站塔1的10 m风速与目标站塔2的100 m风速相关性比方案B中两个站点风速相关性小很多，因此误差要偏大很多。为便于观察，本文从8个月的风速时间序列图截取了前15天 (2009年5月21日—6月4日，含2160个 (15 d×24 h×6)10 min平均数据) 的风速时间序列图。
|
|
| 图 5. 参照站塔1的10 m风速和目标站塔2的100 m, 10 m风速散点图和风速时间序列图 Fig 5. Scatter plotand time series of 10 m wind speed ofreference tower and that of 100 m, 10 m of target tower | |
据前文可知，用来建模的前8个月的两组数据均通过了显著性检验，则认为其在统计学上显著相关[19]，可以用作风速推测。方案A：由塔1的2010年1月21日—4月30日的10 m风速利用方差比MCP方法直接推测出后相同时段塔2的100 m风速。方案B：先由塔1的2010年1月21日—4月30日的10 m风速利用方差比MCP方法推测出相同时段塔2的10 m风速，再由塔2的10 m推测风速数据求出Weibull双参数k10=2.2741，c10=6.9792。将k10，c10，塔2 10 m风速，z=100 m代入式 (9)，就能推算出塔2的100 m风速。图 6则是按照方案A和方案B推测步骤得到的结果与真实测量风速时间序列。
|
|
| 图 6. 方案A和方案B推测结果与实测值比较 Fig 6. Plots of wind speed time series for Method A and Method B with observations | |
由图 6可以看出，方案B明显比方案A推测结果更接近于实测值，方案A的推测结果与实测值相差很大，方案B的推测结果与实测值之间不仅未出现滞留和超前的现象，峰值推测也比较准确。为了更具说服力，本文采用4种检验因子来分别检验推测结果：① 推测值和实测值的相关系数R；② 长期平均风速检验因子m1；③ 卡方分布拟合优度检验因子m2；④ 年度发电总量检验因子m3(年度发电量的计算选取具有代表性、唯一国产化的风力发电机，型号为FD-24-200，其涡轮机启动风速 (启动涡轮机所需的最小风速) 为4.5 m·s-1，额定风速 (达到额定功率的最小风速) 为14 m·s-1，截止风速 (为保护涡轮机而停止转动的最小风速) 为28 m·s-1，额定功率 (平稳发电量) 为200 kW[17]。将数据分组，每一组600个数据，分别计算m1，m2，m3并求出相应的平均值和标准差[20]，统计结果见表 4。
|
|
表 4 两种方案预报结果检验因子统计 Table 4 Summary of results for Method A and Method B |
从表 4可以看到，方案B中4种检验因子的检验结果平均值和标准差均较方案A理想，m1和m3的平均值较方案A更接近于1，标准差比方案A更接近于0。m2平均值和标准差都比方案A接近于0。这表明修正的MCP方法较方差比MCP方法直接推测更接近于实测值，它在高度上起到了一定的修正作用，降低了MCP方法由于海拔高度的影响造成较大误差的局限性。
5 结论与讨论本文利用内蒙古锡林郭勒草原地区两座风能测风塔近1年6层高度上的风杯资料，得出风速的Weibull双参数分布k，c随高度变化的拟合关系，将变化规律加入到方差比MCP方法中，得到了一种能够修正由于高度引起推测偏差的修正的MCP方法，并选取了4种检验因子与方差比MCP方法推测结果进行了比较。得到以下结论：
1) 这种考虑了海拔高度带来的误差的修正的MCP方法推测结果比方差比MCP方法推测结果更接近于实测值，无论在相关系数、长期风速平均值检验因子，还是卡方拟合优度、年度发电量的检验因子检验结果中都得到了较好的结果。
2) 修正的MCP方法不但可以使用自动观测站的观测数据，还可以利用较高的气象塔或风廓线雷达、探空气球等更先进的手段进行加密观测。得到Weibull双参数分布k，c参数随高度分布更精确的变化规律，这样能运用到难于测量的地区，服务于风能评估。
本方法建立在风速序列服从Weibull双参数分布的基础之上，但研究表明约有低于5%的风速并不服从Weibull双参数分布，因此这种修正的MCP方法具有一定的适用范围。目前仍未找到风速不符合Weibull双参数分布的原因和规律，所以很难罗列出本方法精确的适用范围。因此需要在使用之前，对当地的风速进行Weibull分布检验，若不符合，本文方法可能失效。由于数据有限，未能考虑到季节变化的影响。另外，如上述风能测量塔只有6层数据，最高只延伸到100 m，拟合时会有误差出现并存在局限性，若希望得到精密的拟合结果则需要加密观测。
| [1] | Derrick A. Development of the Measure-correlate-predict Strategy for Site Assessment. Proc BWEA, 1992: 259–265. |
| [2] | Derrick A. Development of the Measure-correlate-predict Strategy for Site Assessment. Proc EWEC, 1993: 311–315. |
| [3] | Nielsen M, Landberg L, Mortensen N, et al. Application of Measure-correlate-predict Approach for Wind Resource Measurement. Proc EWEA, 2001: 12–16. |
| [4] | Riedel V, Strack M, Waldl H P. Robust Approximation of Functional Relationships Between Meteorological Data: Alternative Measure-correlate-predict Algorithms. Proc EWEA, 2001: 17–18. |
| [5] | Draper N R, Smith H. Applied Regression Analysis. New York: John Wiley and Sons Inc, 1966. |
| [6] | Woods J, Watson S. A new matrix method of predicting long-term wind roses with MCP. Journal of Wind Engineering and Industrial Aerodynamics, 1997, 66, (2): 85–94. DOI:10.1016/S0167-6105(97)00009-3 |
| [7] | Mortimer A A. A New Correlation Prediction Method for Potential Windfarm Sites. Proc BWEC, 1994: 20–27. |
| [8] | 杨萍, 刘伟东, 仲跻芹, 等. 北京地区自动气象站气温观测资料的质量评估. 应用气象学报, 2011, 22, (6): 706–715. DOI:10.11898/1001-7313.20110608 |
| [9] | 牟聿强, 王秀丽, 别朝红, 等. 风电场风速随机性及容量系数分析. 电力系统保护与控制, 2009, 37, (1): 65–70. DOI:10.7667/j.issn.1674-3415.2009.01.013 |
| [10] | Justus G. Nationwide assessment of potential output from wind-powered generators. J Applied Meteor, 1976, 15, (7): 673–678. DOI:10.1175/1520-0450(1976)015<0673:NAOPOF>2.0.CO;2 |
| [11] | Dorvlo A. Estimating wind speed distribution. Energy Conversion and Management, 2002, 43, (17): 2311–2318. DOI:10.1016/S0196-8904(01)00182-0 |
| [12] | Lun I Y F, Lam J C. A study of Weibull parameters using long-term wind observations. Renewable Energy, 2000, 20, (2): 145–153. DOI:10.1016/S0960-1481(99)00103-2 |
| [13] | Seguro V J, Lambert T W. Modern estimation of the parameters of the Weibull wind speed distribution for wind energy analysis. Journal of Wind Engineering and Industrial Aerodynamics, 2000, 85, (1): 75–84. DOI:10.1016/S0167-6105(99)00122-1 |
| [14] | Genc A. Estimation of wind power potential using Weibull distribution. Energy Sources, 2005, 27, (9): 809–822. DOI:10.1080/00908310490450647 |
| [15] | 徐宝清, 田德, 王海宽, 等. 风速的Weibull分布参数确定方法研究. 农业工程学报, 2007, 23, (10): 31–33. DOI:10.3321/j.issn:1002-6819.2007.10.006 |
| [16] | 杨维军, 王斌. 二参数Weibull分布函数对近地层风速的拟合及应用. 应用气象学报, 1999, 10, (1): 119–123. |
| [17] | 徐大海, 朱蓉. 大气平流扩散的箱格预报模型与污染潜势指数预报. 应用气象学报, 2000, 1, (4): 1–7. |
| [18] | Tony Burton. 风能技术. 武鑫, 谷海涛, 李海东, 等译. 北京: 科学出版社, 2007. |
| [19] | 王国复, 徐枫, 吴增祥. 气象元数据标准与信息发布技术研究. 应用气象学报, 2005, 16, (1): 114–121. DOI:10.11898/1001-7313.20050115 |
| [20] | 许小峰, 胡欣, 王卫丹, 等. 国内大气科学发展状况及优先领域分析. 应用气象学报, 2006, 17, (6): 657–664. DOI:10.11898/1001-7313.20060603 |