风电项目的第一步往往不是选风机而是搞清楚场址的风到底怎么样。我在好几个前期评估项目里拿到过从测风塔导出的历史风力数据Matlab里一读缺失值、野值、塔影效应、采集器时钟漂移各种问题扎堆出现。这篇内容我会从“拿到测风塔原始数据”开始讲怎么用 Matlab 完成从导入、清洗、计算风功率密度到威布尔拟合、发电量估算的完整流程把里面值得注意的坑和经验一并写清楚。如果你正在做风能资源评估、可再生能源前期选址或者只是拿到了一个测风塔数据集不知道如何下手这篇内容就是按我实际工作的思路整理的你可以直接照着跑。1. 测风塔数据从采集器到 Matlab先搞明白你手里这堆数是什么1.1 典型的测风塔配置和输出结构测风塔也叫气象塔是风资源评估最重要的数据来源。一个规范配置的测风塔通常会在多个高度安装风速计和风向标比如 10m、30m、50m、70m甚至更高。除了风速风向绝大多数塔还会配温度传感器、气压传感器和湿度传感器这些数据直接关系到空气密度修正后续计算风功率密度时根本绕不开。数据采集器一般按 10 分钟的平均间隔记录一组数据。每 10 分钟一组一天就是 144 组一年下来一个通道大约 52560 个数据点。原始记录里通常不只有平均风速还有最大风速、最小风速、标准方差甚至每个 10 分钟窗口里的逐秒采样瞬时值。不同采集器的存储格式差别很大有的是直接 CSV有的是 NRG、Campbell 等采集器特有的格式还有的会导出为 Excel 表格。首次拿到数据时我建议先建一个“数据地图”高度层风速通道风向通道是否含最大值/最小值是否含标准差70m轮毂高度附近WS70_avgWD70_avg有有50mWS50_avgWD50_avg有有30mWS30_avgWD30_avg有有10mWS10_avgWD10_avg有有做完这张表你就知道自己手里有哪些通道哪些高度可用于风切变分析哪些通道能用来算湍流强度。没有这张地图后面每一步都会很被动。1.2 从原创采集文件到可分析表的常规转换不同采集器导出的文件第一行往往是变量名中间还可能穿插着注释行或设备日志。直接 readtable 经常会把时间列读成文本风速列读成 cell气压列带个“NaN”字符串总之第一步不会太干净。我的做法是分两步走先用人工过一遍文件头把设备日志、双标题行清理掉只保留表头和数据行。然后用 readtable 统一读取再用 detectImportOptions 把各列类型显式指定好。% 读取测风塔CSV数据 opts detectImportOptions(tower_data.csv); % 显式指定各列类型 opts setvartype(opts, {time,WS70_avg,WD70_avg,WS50_avg,WD50_avg,Temp,Press}, ... {datetime,double,double,double,double,double,double}); T readtable(tower_data.csv, opts); % 检查读取结果 head(T)有一个非常容易掉的坑时间列在 Excel 里显示为“2023-05-01 00:10”但读取进来之后格式可能是“01/05/2023 00:10”如果直接用小时索引处理顺序完全错乱。因此读取后立刻统一时间列格式并给时间列排序T sortrows(T, time); T.time dateshift(T.time, start, minute); % 对齐到整分钟很多数据采集器还会把没有风速记录的时间留成空格或 -9999这一步顺手统一换成 NaN方便后续清洗。做完这些数据才算真正进到 Matlab 的“操作台”上。1.3 时间戳连续性检查采集器时钟漂移不得不防测风塔在野外常年运行太阳能供电、冬季低温采集器的实时时钟经常出现漂移。有些塔几个月后会慢几分钟甚至慢到半小时。这种漂移对于长期平均风速影响不大但对逐时变化、逐日变化和极端风事件分析会直接造成“时间错位”。我从实际项目里总结了一个快速检查法用 diff(T.time) 检查相邻两行时间差是否为 10 分钟把所有非 10 分钟的间隔全列出来。timeGap minutes(diff(T.time)); badGap find(abs(timeGap - 10) 0.01); if isempty(badGap) disp(时间连续性正常); else disp([发现 , num2str(length(badGap)), 处间隔异常]); disp(T.time(badGap(1:min(end,5)))); end如果异常间隔的区域是小范围且只有一分钟级别的偏移一般可以在后续清洗时做重采样对齐。如果漂移严重就必须联系测风塔运维单位调取现场维护日志确认哪些时间段不可靠。这部分工作是数据代表性的底线不查清楚后面所有分析结果都会存疑。2. 数据清洗实操野值、缺测、塔影效应逐个干掉2.1 风速负值、超量程和逻辑错误怎么判定测风数据的清洗不能只靠“看着不对劲就删”。我一般先做三道硬过滤再做一道逻辑过滤。硬过滤风速小于 0直接判定为野值。理论上风速不可能为负出现负值往往是传感器信号短路或采集器故障。风速大于传感器量程。大多数测风杯的极限量程是 60m/s 或 75m/s超过量程上限的读数直接剔除。风向小于 0 或大于 360同样剔除。逻辑过滤需要结合通道关系。比如 70m 高度风速低于 10m 高度风速且相差特别悬殊不一定是错但如果是稳定低于 20% 以上很可能是低层传感器结冰或高层传感器轴承磨损这类数据要单独标记不能盲目删除。下面这一段清洗代码可以直接套用% 硬过滤 T.WS70_avg(T.WS70_avg 0 | T.WS70_avg 60) NaN; T.WS50_avg(T.WS50_avg 0 | T.WS50_avg 60) NaN; T.WD70_avg(T.WD70_avg 0 | T.WD70_avg 360) NaN; % 逻辑检查70m 风速应为塔内最高若低于 10m 风速 30% 以上标记为可疑 idx_logic (T.WS70_avg 0.7 .* T.WS10_avg) ~isnan(T.WS10_avg); T.WS70_avg(idx_logic) NaN;不过要注意逻辑过滤不能做得太激进。强逆温层结、低空急流、局地山谷风环流都可能让上层风速低于下层这种情况虽然少见但不是错误。我在项目里一般只标记不删除后续通过时间切片确认是否为天气系统造成的正常现象。2.2 卡死值和跳变值的识别技巧测风杯轴承卡死是野外塔最常见的机械故障之一。卡死时风速通道会持续输出同一个值持续几小时甚至几天乍一看曲线很正常但标准差接近 0。判断卡死值的方法很简单计算连续相同值的持续时间超过 30 分钟3 个连续采样点就该怀疑了。% 找出连续卡死 3 个点以上的区间 runLen 1; badRun false(height(T), 1); for i 2:height(T) if T.WS70_avg(i) T.WS70_avg(i-1) runLen runLen 1; else if runLen 3 badRun(i-runLen:i-1) true; end runLen 1; end end % 卡死区间置为 NaN T.WS70_avg(badRun) NaN;跳变值则是相邻两个点风速突然变化超过一定阈值。10 分钟平均风速一般不会瞬间变化超过 5m/s 到 8m/s除非是阵风锋面或者雷暴出流这才可能发生剧烈突变。所以跳变值也不能一删了之我通常是先把突变点标记出来再结合气压、温度变化判断是否属于真实天气过程。一个比较实用的替代思路是先画全年的风速时间序列散点图把肉眼可见的“尖刺”区间放大看确认是孤立点还是连续变化。孤立点直接剔除连续变化则保留并记录说明。2.3 塔影效应和表头朝向剔除风向数据里藏着的系统性偏差测风塔的塔身、拉线、横臂都会对气流产生扰动。当风向吹向塔体的某一侧时塔身下游的尾流会让测风杯的读数系统性偏低这就是“塔影效应”。如果直接把所有风向扇区的风速都纳入统计平均风速会被轻微拉低风功率密度会被低估。一条最常规的做法是根据测风塔安装日志里的传感器朝向把风向在塔身方位 ±30° 范围内的数据剔除。比如测风塔的横臂朝向为 90°正东当风向在 60°–120° 区间时风速数据不可用。% 假设传感器朝向为 90° 和 270° 两侧塔身两侧 60°–120° 与 240°–300° 为受扰扇区 towerShadowSec [60 120; 240 300]; isShadow (T.WD70_avg towerShadowSec(1,1) T.WD70_avg towerShadowSec(1,2)) | ... (T.WD70_avg towerShadowSec(2,1) T.WD70_avg towerShadowSec(2,2)); T.WS70_avg(isShadow) NaN;前面的逻辑过滤会引发一个问题大量数据被置成 NaN 后剩余数据量可能不够。行业里对有效数据率的要求一般是 90% 以上低于这个水平年度代表性就要打问号。所以塔影剔除法务必要记录剔除比例如果超过 25%说明塔的布局本身有问题不能硬着头皮继续分析。2.4 缺测数据的插补短期线性插值、长期需要相关回归插补策略取决于缺失数据的时长。缺 1–2 个 10 分钟点直接用前后线性插值完全可以误差很小。缺 2–3 小时线性插值还行但要谨慎。缺一整天甚至更久唯一的可靠办法是用邻近高度的风速相关性做回归或者用测风塔附近长期气象站数据做拟合。% 短时缺测线性插值 T.WS70_avg fillmissing(T.WS70_avg, linear, SamplePoints, T.time); % 长时间缺测如果 70m 缺得太多用 50m 回归到 70m 的幂律关系 validIdx ~isnan(T.WS70_avg) ~isnan(T.WS50_avg); p polyfit(log(T.WS50_avg(validIdx)), log(T.WS70_avg(validIdx)), 1); T.WS70_avg(isnan(T.WS70_avg)) exp(polyval(p, log(T.WS50_avg(isnan(T.WS70_avg)))));需要特别留意插补数据不代表真实测量值在最终报告中必须单独标注插补比例。能作为依据的气象标准里插补数据比例过高时风功率密度的置信度会明显下降这一点在后面的不确定性分析里会继续使用。3. 平均风速、风功率密度与湍流强度三个核心指标一次算清3.1 逐时平均风速与有效样本数量风速的统计基础是 10 分钟平均序列但做资源评估时通常还会按小时、月、年三个尺度统计。小时平均可以滤掉阵风脉动的短时干扰月平均可以看到季节变化年平均则是场址长期资源水平的基本盘。% 按小时平均 T.hourly dateshift(T.time, start, hour); hourlyWS groupsummary(T, hourly, mean, WS70_avg); % 年平均 annualWS groupsummary(T, year, mean, WS70_avg);这里有个非常关键的原理平均风速必须在“小时风速”和“10分钟风速算术平均”之间选择正确的口径。规范做法是用 10 分钟平均风速序列直接参与后续风功率密度计算因为风功率密度与风速立方成正比而风速立方平均不等于平均风速的立方。要是直接用年平均风速去估算发电量结果会严重偏低。3.2 风功率密度公式与空气密度修正风功率密度表示单位面积上气流的功率公式是[ P \frac{1}{2} \rho V^3 ]实际计算时用 10 分钟平均风速序列逐个点求风速三次方再取平均最后乘上 0.5 和空气密度。Matlab 写法% 计算空气密度根据实测气压和温度 Rd 287.05; % 干空气气体常数 J/(kg·K) T.K T.Temp 273.15; T.rho T.Press ./ (Rd .* T.K); % Press 单位 Pa rho_mean mean(T.rho, omitnan); % 风功率密度 W/m2 P_wd 0.5 * rho_mean * mean(T.WS70_avg.^3, omitnan);空气密度往往是新手会忽略的地方。标准空气密度 1.225 kg/m³ 只适用于海平面、15°C 的标准情况。高原风电场空气密度可能只有 1.05 甚至更低直接套 1.225风功率密度会被高估 10% 以上。如果你手上的测风塔没有气压传感器可以用海拔粗略估算[ \rho \approx 1.225 \times e^{-0.0001 \times z} ]其中 z 是海拔单位米。这个近似公式在 3000m 以下误差可以控制在 2%–3% 左右前期预评估够用施工图阶段必须要实测气压。3.3 湍流强度计算直接决定风机等级选型湍流强度反映风速脉动的剧烈程度定义为标准差除以平均风速[ I \frac{\sigma V}{\overline{V}} ]计算时需要 10 分钟窗口内的高频采样记录或者采集器输出的风速标准差通道。如果没有标准差通道只能用相邻 10 分钟平均值的变差来近似但这会低估真实湍流强度不太推荐。% 用采集器输出的标准差通道计算湍流强度 T.I70 T.WS70_std ./ T.WS70_avg; T.I70(T.I70 0.6) NaN; % 明显错误值剔除 I70_mean mean(T.I70, omitnan);为什么要专门算湍流强度因为风机选型时IEC 等级直接和湍流强度挂钩比如 IEC 1A 类风机要求 15m/s 风速下的湍流强度不超过 0.16。如果场址湍流强度太高风机载荷安全性不达标要么降低等级选型要么就得避开高湍流区域排布机位。这个参数对风机排布方案的影响比很多人想象中大得多。4. 威布尔分布拟合与风玫瑰图判断风资源质量的两把尺子4.1 威布尔分布参数估计最小二乘与极大似然风速频率分布是描述风资源最重要的统计特征。威布尔分布有两个参数形状参数 k 和尺度参数 c。k 值决定了风速分布的“尖陡”程度k 越小风速分布越矮胖说明风速变化范围大k 越大风速越集中比如 2–3 左右。c 值近似代表年平均风速换算后的水平。实际拟合思路把所有有效风速从小到大排序用中位秩公式计算经验累积频率。对威布尔累积分布函数做线性化变换用最小二乘拟合直线斜率和截距反推 k 和 c。v sort(T.WS70_avg(~isnan(T.WS70_avg))); n length(v); F (1:n) / (n 1); % 中位秩 x log(v); y log(-log(1 - F)); p polyfit(x, y, 1); k_fit p(1); c_fit exp(-p(2) / p(1));如果想省事Matlab 自带的 wblfit 可以直接做极大似然估计params wblfit(v); c_mle params(1); k_mle params(2);两种方法各有侧重最小二乘比较直观可控极大似然统计效率更高。我建议两个都算一遍对比如果差值超过 5%说明数据里还有未清理的异常段或者双峰分布比如季风区冬夏风向差异极大这时候单纯拟合单一威布尔分布就不够精确要考虑分季节或分风向扇区拟合。4.2 风玫瑰图主导风向比平均风速更能说明问题风玫瑰图的核心是展示不同方向的频率和风速贡献。选址阶段主要关注两点主风向是否集中主风向上的平均风速是否够高。如果主风向是东南风但强风速全在西北方向机位排布和道路设计就得有不同的侧重。Matlab 里用 polarhistogram 和 polaraxes 就能画dir_rad deg2rad(T.WD70_avg(~isnan(T.WD70_avg))); figure; polaraxes; polarhistogram(dir_rad, 16, Normalization, probability); title(70m 高度风玫瑰图);实际项目里我一般画两张图一张是风向频率玫瑰一张是风速贡献玫瑰。频率玫瑰告诉你哪边风多风速贡献玫瑰告诉你哪边风带来的能量多后者对排布的影响更直接。很多风电场排布方案失败就是只看了频率玫瑰没看能量玫瑰。4.3 风切变指数与轮毂高度风速推算测风塔各高度层的平均风速可以用来计算风切变指数 α[ \alpha \frac{\ln(V_2 / V_1)}{\ln(z_2 / z_1)} ]alpha log(mean(T.WS70_avg,omitnan) / mean(T.WS30_avg,omitnan)) / log(70 / 30);α 的值直接决定从已知高度推算轮毂高度风速的可信度。平坦开阔地形 α 大约在 0.10–0.20 之间复杂山地可达 0.3 甚至更高。如果测风塔最高层是 70m而风机轮毂高度是 100m就需要用 α 外推。hub_h 100; ref_h 70; V_hub mean(T.WS70_avg, omitnan) * (hub_h / ref_h)^alpha;需要警示一点α 不是恒定常数它会随风速区间变化夜间稳定层结下 α 可能异常偏大。所以在做外推时最好分风速段统计 α再取轮毂高度对应风速段的值而不是用全年平均 α 一把梭。否则 100m 高度风速容易被高估。5. 发电量估算与机位排布资源评估最终要回答投资方的问题5.1 功率曲线法的 AEP 理论计算资源评估做完了最终要落到“这个场址一年能发多少电”上。最常用的方法是功率曲线法把风速频率分布威布尔分布与风机功率曲线逐风速点相乘再乘以全年小时数。% 假设功率曲线以表格格式存储列1为风速列2为功率kW powerCurve readtable(turbine_power_curve.csv); v_pc powerCurve.V; P_pc powerCurve.kW; % 用威布尔概率密度函数生成风速分布 v_bin 0:0.5:30; f_v wblpdf(v_bin, c_fit, k_fit); % 确保功率曲线覆盖 v_bin 范围 P_interp interp1(v_pc, P_pc, v_bin, linear, 0); % 计算年发电量 (单位 kWh) AEP_theoretical 8760 * sum(f_v .* P_interp);注意两个细节。第一功率曲线的风速一般是标准空气密度下的值如果场址空气密度低于 1.225需要做密度修正。工程上常见的做法是按实际密度与标准密度的比值等比修正功率曲线密度越低同等风速下功率越低。第二风机的切入切出风速区间会直接影响有效发电小时数如果风速分布的高风速占比大切入切出后的截断效应更明显。5.2 综合折减系数把理论值拉回现实理论 AEP 是“风机全年无故障、无停机、无尾流损失、无控制损失、电网永远消纳”的理想值。现实里必须乘一个折减系数。我常用的一组参考值是折减项典型范围说明可利用率和维护停机93%–97%含故障检修、定期维护尾流损失5%–12%与机位间距、排布有关电气与输变电损失2%–5%含线损、变压器损耗控制与偏航损失1%–3%含偏航误差、功率限制叶片污染损耗1%–3%沙尘、盐雾、结冰影响场址不确定性5%–10%测风年代表性、长期订正误差综合下来年净发电量大约要乘 0.75–0.85。如果某个项目算出来的综合折减超过 90%要么你的可利用率数据异常优秀要么你漏了某个大项要回去重新核对。netFactor 0.95 * 0.90 * 0.97 * 0.98 * 0.98 * 0.90; AEP_net AEP_theoretical * netFactor;5.3 从资源评估到机位排布两个容易被忽视的图机位排布图纸上通常看两样东西一个是风速分布色块图一个是湍流强度色块图。风速高的区域优先布置风机这是直觉湍流强度高的区域尽量避开这是安全。如果高湍流区和高风速区重叠就需要逐机位做载荷安全评估不能拍脑袋决定。我处理过的项目中出现过两个典型情况。第一个场址西南角年平均风速高出东部 15%但西南角紧邻断崖湍流强度超过 0.25最终那排机位全部后移了 200m发电量下降但载荷安全性达标。第二个某个平坦湖岸场址主风向与岸线平行排布时按正方形网格布置结果尾流损失达到 14%改成顺主风向错列排布后降到 8%年发电量提升明显。这些都说明资源评估的最终交付物不只是报告还包括真正能指导排布的分析图和数据表。行业里的成熟做法是用 CFD 或尾流模型做多轮迭代优化但在没有商业软件的情况下用 Matlab 基于威布尔分布、风速分布和功率曲线做一版初步的静态估算已经能给项目可行性判断提供足够的方向性参考。我在实际项目里最大的体会是测风塔数据的处理流程尽量标准化清洗规则、插补方法、塔影扇区、折减系数全都用脚本固化下来。这样每次拿到新项目的数据只需改文件路径和塔结构参数就能快速产出结果。风电项目的周期往往很短前期评估要快、要准、要可追溯一套可复用的 Matlab 处理流程比临时写脚本省下不止一半时间。如果你手头也有一套测风塔数据建议先按第 1 节的方法建好数据地图再按清洗、插补、指标计算、分布的流程走一遍。遇到什么奇怪的数据现象欢迎来交流。
