★ FEATURED ARTICLE 气象塔风电行业里叫测风塔出来的历史风力数据看起来就是一堆时间和风速的表格但真正能把它变成风能资源评估结论的人并不算多。我刚接触这个领域的时候也犯过“拿数据就求平均”的错结果被前辈按在工位上从头教了一遍测风数据从传感器到报告之间隔着一整条质量控制、统计建模、发电量估算的链路。这篇博文就把这个过程完整走一遍用Matlab从导入、清洗气象塔历史风力数据开始一直到平均风速、威布尔分布、风功率密度、风玫瑰图和年发电量估算适合刚入行的风电开发工程师、新能源方向的研究生以及所有想弄明白风能资源评估数据怎么变成决策依据的人。1. 测风塔数据风电项目选址的“底牌”1.1 为什么风资源评估不能靠“拍脑袋”做风电项目前期开发设备的选型、机位的排布、投资收益的测算几乎全部压在一个东西上场址的风资源条件。一个风电场的投资规模动辄上亿如果资源评估偏差10%整个项目的收益率曲线都会明显偏移。而风资源评估最硬的数据来源就是立在现场的气象塔实测数据——它比卫星反演、数值再分析网格更接近未来风机轮毂高度处的真实气流。气象塔之所以叫“塔”是因为它通常要立到计划风机轮毂高度附近比如70米、90米、100米塔身上在不同高度布置风速计、风向标、温度、气压传感器。数据采集器会连续记录风速、风向等要素通常是10分钟一条记录。这个10分钟平均数据就是后续所有分析的原料。和测风塔数据相比其他数据源可以作为参考或长期订正但无法替代现场实测的地位。这篇文章想讲清楚的核心就是这批10分钟原始数据从“Excel表格”变成“资源评估结论”的完整过程。不管你是只写过几行Matlab的初学者还是已经做过几个项目但想理顺流程的工程师都可以照着下面的步骤走。1.2 测风塔数据的构成与“一塔一年”行业惯例拿到气象塔数据先别急着写代码。先看清楚手里是什么。典型的测风塔数据文件列结构大致是这样的字段含义典型单位/范围time记录时间yyyy-mm-dd HH:MM:SSws70 / ws5070米、50米高度平均风速m/s0~60wd70 / wd5070米、50米高度风向°0~360temp空气温度°Cpres大气压力hPahum相对湿度%行业里做风资源评估默认的“起步数据”是一整年的完整测风记录通常要求跨越四个季节数据完整率不低于90%。为什么是一年因为风资源有明显季节性只有覆盖完整年周期的数据才能代表场址的真实水平。如果只有半年数据评估结果会有明显的季节偏差。数据本身也不是拿来就能用的。10分钟数据里常见的问题包括传感器结冰导致的连续重复值、雷击或线路故障导致的野值、设备停机导致的缺测、风向标被塔影遮挡导致的错误风向。这些都在“数据处理”这一步解决。2. 一次完整的风资源评估工作流从测风塔到Matlab2.1 工作流总览与每个环节的意义我习惯把风能资源评估数据研究拆成八个环节数据获取、数据导入、质量控制、数据订正与插补、统计特征分析、垂直外推、发电量估算、报告输出。其中数据获取靠现场设备和采集器完成我们拿到手的是第二环节开始。质量控制是所有环节里最枯燥但又最重要的一步直接决定后面所有统计数字的可信度。统计特征分析包括平均风速、风向频率、威布尔分布参数、风功率密度、湍流强度等这些指标共同构成风资源评估报告的核心表。发电量估算则是在统计特征基础上结合风机功率曲线计算机组理论发电量再扣除尾流、利用率等损失得出风电场层面的发电量预期。从数据分析角度讲Matlab做这批时序数据处理很顺手。内置的timetable、retime、wblfit这些函数几行代码就能完成数据聚合和概率分布拟合画出来的图也足够进报告。当然用Python也能做同样的事只是工具不同分析逻辑完全一致本文以Matlab为例。2.2 环境准备与目录约定工欲善其事必先利其器。做数据研究之前先把项目目录搭好。我自己的习惯是下面这套结构project_root/ ├── data/ │ ├── raw_data/ # 原始测风数据只读不修改 │ └── processed_data/ # 清洗后的数据 ├── scripts/ # 分析脚本 ├── results/ # 统计结果表 └── figures/ # 输出图片原始数据目录只放从现场拷贝回来的文件保持“原样”处理后的数据单独放这样做的好处是数据可追溯。比如甲方问“这个异常值是不是删多了”你可以随时回到原始文件重新检查。Matlab版本建议R2019b以上需要统计和机器学习工具箱主要是为了用wblfit做威布尔分布拟合。没有这个工具箱也不慌后面会给出一个备用的手动拟合方法。2.3 检查文件与数据规模的第一个脚本拿到数据目录后先用一个最简单的脚本确认文件数量和格式clear; clc; rawDir fullfile(pwd, data, raw_data); files dir(fullfile(rawDir, *.csv)); fprintf(发现测风数据文件 %d 个\n, numel(files));这一步看着简单但很值得做。实际项目中有时数据会按月拆成12个文件有时一个文件包含全年有时还混着不同塔的数据。先看清楚文件规模再写批量导入逻辑能避免很多不必要的返工。3. 数据导入与清洗的实操细节坏数据会让整个评估前功尽弃3.1 批量导入CSV与时间戳解析原始CSV文件的格式千差万别有的列名是中文有的时间格式是“2023-01-01 00:00”有的是“2023/1/1 0:00”。我用统一的方式处理先读入再统一列名最后把时间字符串解析成datetime类型。opts detectImportOptions(filename); opts.VariableNames {time,ws70,wd70,ws50,wd50,temp,pres,hum}; raw readtable(filename, opts); % 统一时间格式 raw.time datetime(raw.time, InputFormat, yyyy-MM-dd HH:mm:ss); raw sortrows(raw, time);如果遇到解析报错多半是时间格式不一致可以先检查一下时间字符串的规律% 查看前20个时间字符串 head(raw.time, 20)时间解析这块经验是“宁慢勿错”。有些CSV文件的时间列带时区后缀或者包含“0:00”这种不补零的格式统一处理之前必须先看清楚实际字符串格式。解析完成后用diff检查时间间隔是否连续dt diff(raw.time); fprintf(最小采样间隔%s\n, char(min(dt))); fprintf(最大采样间隔%s\n, char(max(dt)));10分钟采样频率的数据最小采样间隔应该是10分钟。如果最大间隔远远大于10分钟说明存在缺测段需要在后面重点处理。3.2 质量控制逻辑过滤野值、标记卡滞、记录完整率质量控制是整个处理环节的“底线”。我的处理规则很简单但每条都有明确依据检查项判断条件处理方式负值/超量程风速0或60 m/s置为NaN风向超范围风向0或360置为NaN温度异常温度-40℃或60℃置为NaN传感器卡滞连续12个以上记录完全相同置为NaN相邻跳变相邻时刻风速变化20 m/s且不回落按野值标记代码实现方式并不复杂% 基础上下限过滤 valid1 raw.ws70 0 raw.ws70 60; valid2 raw.wd70 0 raw.wd70 360; valid3 raw.temp -40 raw.temp 60; % 连续重复值判断传感器卡滞 dup [false; diff(raw.ws70) 0]; % 统计连续重复段长度 runLen 0; dupRun false(size(raw.ws70)); for i 2:height(raw) if diff(raw.ws70(max(1,i-1):i)) 0 runLen runLen 1; else runLen 0; end if runLen 12 dupRun(i) true; end end % 总体有效标记 valid valid1 valid2 valid3; valid(dupRun) false; % 无效值置NaN raw.ws70(~valid) NaN; raw.wd70(~valid) NaN; % 统计完整率 validRate mean(valid) * 100; fprintf(数据完整率%.2f%%\n, validRate);这段逻辑里有个容易忽略的点风速为0在物理上是可能出现的无风天气所以下限判断用“0”不要写成“0”否则会把真实静风数据误删。风向判断也一样0°代表北风是有效值。这里说个实操体会质控处理不只是“删数据”删掉的每一段都应该能解释原因。比如冬天塔上传感器结冰风速数据会出现连续几小时完全相同的恒定值这段数据被标为无效后最好在报告里注明“冬季结冰导致传感器卡滞”。甲方、评审专家看报告时非常在意这些细节。3.3 从10分钟原始数据到小时平均序列原始数据是10分钟一条做统计分析和后续发电量估算时很多场景会用到小时平均风速。聚合方法用timetable的retime函数一行代码的事。但这里有个陷阱风向不能直接做算术平均。0°到360°的边界是连在一起的如果某个小时内有几条记录是350°、355°、0°、5°算术平均值会算出个错误方向比如用普通mean会得出约177°。正确做法是把风向转成u、v分量再平均TT table2timetable(raw(:, {time,ws70,wd70,temp,pres,hum})); TT rmmissing(TT); % 先删除含无效值的整行 % 将风向转为u、v分量 TT.u TT.ws70 .* cosd(TT.wd70); TT.v TT.ws70 .* sind(TT.wd70); % 聚合为小时平均 TT_hour retime(TT, hourly, mean); % 从u、v分量合成小时平均风向 TT_hour.wd70_avg atan2d(TT_hour.v, TT_hour.u); TT_hour.wd70_avg(TT_hour.wd70_avg 0) TT_hour.wd70_avg(TT_hour.wd70_avg 0) 360;聚合这块还有一个容易被忽略的业务规则一个小时内如果有效10分钟记录太少算出来的小时平均值是不可信的。行业里常见要求是每小时至少要有4条有效记录才保留该小时平均否则标记为NaN。实现上可以用retime加“聚合函数”也可以简单用movsum统计窗口内有效值数量再筛选。这一步保证了后面每个小时风速都有“足够的观测支撑”。提示风向平均值算出来后如果要切回360°制表示务必处理负值否则会出现风向-5°这种看起来就不严谨的数值。4. 风能分析的三个硬指标平均风速、威布尔分布与风功率密度4.1 平均风速与垂直外推风速统计最基础的就是平均风速。直接对小时序列求均值即可meanSpeed mean(TT_hour.ws70, omitnan); maxSpeed max(TT_hour.ws70);但实际评估中测风塔测风高度不一定是未来风机轮毂高度比如塔高70米风机轮毂高度100米。这时候需要做垂直外推。工程上最常用的是幂指数公式[ v_2 v_1 \left( \frac{h_2}{h_1} \right)^{\alpha} ]其中α是风切变指数。如果你在现场70米和50米两个高度都测了风速可以用实测数据估算α% 由50米和70米高度风速估算风切变指数 alpha log(meanSpeed70 / meanSpeed50) / log(70 / 50); % 外推到100米 meanSpeed100 meanSpeed70 * (100 / 70) ^ alpha;如果手里只有单层测风数据就用行业经验值α1/7约0.143。但注意这个经验值属于“海上或平坦地形近似值”复杂地形下用实测两层风速算own alpha会更靠谱。切变指数也直接反映大气稳定性夜间稳定层结下切变很强这些都是估算轮毂高度风速时的细节但很影响最终结果。4.2 威布尔分布拟合比“平均值”更有信息量风资源的概率分布特征行业标准做法是用两参数威布尔分布来描述。风速概率密度函数[ f(v) \frac{k}{A} \left( \frac{v}{A} \right)^{k-1} e^{-(v/A)^k} ]A是尺度参数和平均风速正相关k是形状参数决定分布形态。一般在陆上风电场k值落在1.5~3之间k越大说明风速变化越集中风资源越稳定。Matlab里拟合很简单speedValid TT_hour.ws70(~isnan(TT_hour.ws70)); params wblfit(speedValid); A_fit params(1); k_fit params(2); % 校验用A和k反推理论平均风速与直接均值对比 meanTheory A_fit * gamma(1 1 / k_fit); fprintf(拟合均值%.2f m/s实测均值%.2f m/s\n, meanTheory, mean(speedValid, omitnan));如果理论均值与实测均值差得比较多说明样本分布并不完全符合威布尔模型可能是数据清洗不彻底或者存在极端工况有必要回去检查质控步骤。为什么说威布尔参数比平均风速“更有信息量”因为平均风速相同不代表风速分布相同。A6.5、k2的场址与A6.5、k3的场址年均风速一样但高风速段的累计频率差异明显风机在高风速段发的电更多同样的平均风资源后者等效发电小时可能高出不少。只看平均风速完全看不出这种差异。没有统计工具箱时也可以手动拟合威布尔参数较简单的方法是矩估计[ k \approx \left( \frac{\sigma}{\bar{v}} \right)^{-1.086} ]然后用均值与k反推出A。这个估计精度对工程初判够用但正式报告还是建议用最大似然法。4.3 风功率密度衡量资源“质量”的关键参数风功率密度WPD是判断风资源等级的核心指标公式[ WPD \frac{1}{2} \rho v^3 ]注意正确的算法是用逐时风速三次方的平均值而不是“平均风速的三次方”两者差距很大。因为风速的三次方对高风速段非常敏感用平均风速的三次方会系统性低估风功率密度。空气密度ρ也要用实测数据计算简单方式是用理想气体公式% 用温度、气压估算空气密度kg/m³ rho TT_hour.pres ./ (287.05 .* (TT_hour.temp 273.15)); % 风功率密度 WPD 0.5 * mean(rho .* TT_hour.ws70 .^ 3, omitnan);打个比方平均值是“账面上的资源”风功率密度才是“真正能装进口袋的能量”。我经手的一个模拟项目里两个场址平均风速只差0.3 m/s但风功率密度差距达到80 W/m²直接决定了机型的IEC等级选择。评估报告里WPD是一个绕不开的硬数据。提示WPD对极端风速很敏感质控如果没做干净一个被保留下来的野值就足以让全年WPD明显偏高。算WPD之前务必确保所有风速超过合理阈值的记录都已经被标为NaN。5. 风向与风玫瑰别让“平均值”掩盖方向信息5.1 风向数据的“0°/360°边界陷阱”风向是圆周数据它的统计处理与风速很不一样。前面聚合小时数据的时候已经提醒过风向平均的问题。做扇区统计时同样要小心最典型的是0°和360°其实是一个方向但数值上差了360。处理方式简单粗暴先把360°统一改成0°wd TT_hour.wd70_avg; wd(wd 360) 0;然后做16扇区划分。16扇区就是每22.5°一个区间。用histcounts时要注意边缘误差我的习惯是让0°严格落在第一个扇区中心而不是边界edges -11.25:22.5:348.75; [~, ~, binIdx] histcounts(wd, edges); freq accumarray(binIdx, 1) / numel(wd) * 100;5.2 用Matlab绘制风玫瑰与扇区频率表风玫瑰图是表达风向和风速频率的标准方式。Matlab可以用极坐标直方图快速画figure; polarhistogram(deg2rad(wd), 16, FaceColor, [0.2 0.6 0.9]); title(风向频率玫瑰图);但更严谨的报告用图我一般会自己画扇区柱状图这样速度、频率、平均风速可以按扇区统计清楚做出来的扇区频率表长这样扇区风向区间频率占比扇区平均风速N348.75°~11.25°22.3%6.8 m/sNNE11.25°~33.75°8.7%6.2 m/s............从这张表能直观看出主风向。实际项目中主风向落在哪几个扇区直接决定风机排布方式和自动偏航控制策略。5.3 主导风向对风机排布与尾流评估的影响风玫瑰图不只是“画出来好看”它直接影响微观选址。在主导风向上风机之间需要拉开距离否则前排风机产生的尾流会明显降低下风向风机的来流风速。行业经验是主风向间距通常要留到8~10倍叶轮直径垂直主风方向可以稍密一些。评估阶段从测风塔数据分析出主导风向是风机排布的第一手依据。同时测风塔自身也会产生塔影遮挡。塔身对下风向的传感器有尾流干扰如果发现某个扇区的风速系统性偏低且该扇区和塔身位置吻合这部分数据要在风资源评估中标记为“塔影影响扇区”必要时剔除或修正。这点在现场踏勘和数据分析里都容易被忽略。6. 从数据到决策功率曲线插值、容量系数与年发电量估算6.1 用功率曲线把风速换算成出力资源分析做到这一步终于要和“电”挂钩了。方法是把小时平均风速序列代入某机型标准功率曲线得到每小时的输出功率。以某2MW主流陆上机型为例切入风速3 m/s额定风速约11 m/s切出风速25 m/s。功率曲线文件通常包含风速列和功率列pc readtable(power_curve.csv); % 假设功率曲线列为 wind_speed 和 power_kw powerEst interp1(pc.wind_speed, pc.power_kw, TT_hour.ws70, linear, 0); % 空气密度修正 rho_actual mean(rho, omitnan); powerEst powerEst .* (rho_actual / 1.225);interp1的最后一个参数0表示风速低于切入风速或高于切出风速时输出功率按0处理。这一步比手动if判断干净很多也避免数组维数不匹配。小时功率序列求和就是理论发电量annualEnergy sum(powerEst) / 1000; % kWh - MWh fprintf(年理论发电量%.0f MWh\n, annualEnergy);需要注意的是用小时平均风速估算发电量相对用10分钟数据存在一定偏差。风速的瞬时波动在高风速段带来的发电增量会被平均掉工程初判阶段可以用小时数据但正式评估报告通常建议用10分钟数据或分钟级数据。这是方法精度问题分析流程本身完全一致。6.2 容量系数外行人也能看懂的选址指标容量系数CF定义是年发电量除以机组额定功率×全年小时数[ CF \frac{E_{annual}}{P_{rated} \times 8760} ]如果某2MW机组理论年发电量为500万kWh那么ratedPower 2e3; % kW CF annualEnergy * 1000 / (ratedPower * 8760); fprintf(容量系数%.3f\n, CF);容量系数0.25~0.35是陆上风电场常见的区间。低于0.2的项目投资收益压力会非常大高于0.35的项目在资源评估报告里就是相当优质了。我习惯在报告里同时给出等效满发小时数比如容量系数0.27约等于等效满发2365小时这个说法对投资决策者更直观。6.3 尾流与利用率修正的实际经验需要提醒的是前面算的是单机理论发电量不等于风电场实际发电量。整场发电量还要做一系列折减常见折减项如下折减项典型范围备注尾流损失5%~15%机群排布越密尾流越大场用电与输配电损耗2%~4%集电线路、变压器损耗叶片污染与功率曲线偏差3%~8%污垢、结冰、制造公差停机维护与非计划损失3%~8%检修、电网限电、故障假设尾流8%、损耗3%、维护5%合计折减约15%则整场年发电量约为单机估算值乘机位数再乘0.85。做评估时如果把单机理论发电量直接乘台数会被评审专家一眼看穿这是报告里最常见的低级错误。我在经手的模拟项目X里就亲自把“乐观版”和“稳健版”两套发电量都算过一遍实际运维数据证明稳健版更接近真实值。7. 实操中绕不开的几个坑切塔、补缺与季节偏差7.1 检测时间断点数据完整性问题数据处理完成后一定要检查时间序列中有没有断点。设备停机、采集器重启、卡恢复都可能造成几小时甚至几天的数据缺口。一张全年时间间隔检查表就能看出问题dtHour minutes(diff(TT_hour.time)); gapIdx find(dtHour 90); % 大于1.5小时视为缺口 for i 1:length(gapIdx) fprintf(缺口%s 到 %s历时 %.1f 小时\n, ... TT_hour.time(gapIdx(i)), TT_hour.time(gapIdx(i)1), dtHour(gapIdx(i))); end现实中出现缺口不可怕可怕的是不对缺口做说明。还有些项目会遇到“切塔”情况测风塔因为台风倒了重新立塔后继续测中间缺了两个月。这种情况缺口前后的数据原则上是两个数据集是否拼接使用要很慎重。7.2 缺测补插能补的补不能补的硬着头皮也不能补很多人拿到缺测数据第一反应是用线性插值填满。我的原则是短缺口2小时以内可以补长缺口不能补。原因是风速序列的持续性和随机性没法用线性模型真实还原强行补出一段“看起来连续”的数据实际上是在给评估结论注入虚假信心。Matlab里可以用fillmissing配合MaxGap控制补插范围TT_hour.ws70 fillmissing(TT_hour.ws70, linear, MaxGap, hours(2));这一步把短缺口补上长缺口保持NaN在报告里明确写出“数据完整率92%其中缺失部分未做人工插补”。完整率是底补插是权宜。如果完整率都不到80%建议直接考虑延长测风周期而不是靠插值硬凑一年。7.3 季节代表性与长期订正最后一个大坑是季节代表性。假设项目方只提供了3月到10月的测风数据这段时期平均风速通常低于全年平均因为少了冬春大风季直接把8个月数据说成年平均结果会系统性偏低。反之只有冬季数据结果会偏高。解决思路是结合长期再分析数据把短期实测数据和长期参考数据建立相关关系将短期统计修正到长期平均水平这就是行业里说的“长期订正”或“代表性订正”。具体实现涉及相关分析、线性回归甚至概率映射是一套独立的分析流程这篇文章不展开了但作为风能资源评估数据分析的一部分它必须出现在报告的方法论说明里。提示任何风资源评估报告都应该包含“数据时段与代表性”说明清晰写明实测数据覆盖的时间范围、完整率、缺测原因及订正方法这是专业性的直接体现。我自己现在拿到一份风资源评估报告第一眼看的不是那行平均风速而是“数据完整率”和“质控过程描述”这两栏。只要这两栏写得潦草后面算得再精细我都会打个问号。测风塔数据研究这件事90%的功夫其实花在把数据处理干净、把逻辑讲清楚上Matlab代码只是把分析过程落地的工具。哪怕只是刚接触数据分析的读者从这份流程入手也能避开我最开始踩过的那些坑。
阅读完成 · 觉得有帮助？
