最近在处理某风电场前期选址的测风数据分析工作突然觉得这类活儿挺值得拿出来讲讲。风能资源评估听着是个大词落到具体执行层面其实就是一件事把气象塔行业内常叫测风塔测到的历史风力数据导入Matlab完成清洗、预处理、统计分析最终算出平均风速、威布尔分布参数、风功率密度、主导风向这些核心指标。这篇文章就围绕这个完整链路展开适合刚接触风资源评估的工程师、相关专业的学生以及对新能源数据分析感兴趣的人参考。我会把整个流程拆开讲连踩过的坑一起交代。1. 风能资源评估到底在算什么测风塔的本质是一个架在目标场址上的“长期气象观测站”。塔身上不同高度会安装风速计、风向标、温度探头和气压计一般层高选10m、50m、70m、80m、100m这些典型风机轮毂高度附近。数据采集器通常每10分钟记录一次平均值然后存成时间序列。这个10分钟的平均风速就是后续所有分析的原始素材。为什么必须依赖实测数据而不是气象站数据因为风能资源评估要回答的是“本地风况到底适不适合发电”而不是“这个地区大气候如何”。这两个问题的尺度完全不同。周边气象站的数据虽然年份长、容易拿但海拔、地形粗糙度、山谷走向不一样风况差异非常大。这也是测风塔数据不可替代的原因——它是场址真实风况的采样精度和代表性都是模型再分析数据比不了的。用Matlab来做这件事有它非常实在的好处读入CSV、Excel毫不需要额外写解析模块统计工具箱自带威布尔拟合函数绘图接口能出一整套风资源图表。对比Excel处理几十万行时间序列时的卡顿Python虽然也能做但需要另配一套环境Matlab在这个领域依然是很多设计院和咨询机构的首选工具。整套评估链条有四个核心产出长期平均风速和风切变指数、风速频率分布及威布尔拟合参数、风功率密度、主导风向与湍流强度。这四个指标直接决定了某个风机机型是否适配场址、排布间距怎么拉、理论发电量能到多少。所以后面所有代码都是围绕这四个目标来组织的。2. 气象塔数据的导入与格式处理2.1 先弄清采集器导出的数据长什么样测风塔的数据采集器多数支持直接导出CSV或TXT文件。以我常碰到的某项目X为例导出文件的表头大致长这样时间戳风速70m风向70m风速50m风向50m温度气压2019-01-01 00:007.822486.912526.2923.12019-01-01 00:108.152557.232496.3923.42019-01-01 00:207.572406.442386.1924.0注意到几个细节时间戳一般是字符串格式风速单位是m/s风向单位是度0°表示正北按顺时针增大270°为正西气压单位是hPa或kPa。有些采集器还会附加瞬时风速、最大风速或仪器状态码比如9代表冻结、8代表校验失败这些字段在清洗阶段非常有用千万别一上来就删掉。还有一个容易踩的坑导出文件的编码。很多采集器软件在Windows下导出的是GBK编码Matlab默认用UTF-8读会乱码。处理办法是提前用记事本之类工具另存为UTF-8或者在读取时指定编码参数。2.2 Matlab导入数据的三种方式我最常用的方式是readtable它能把CSV和TXT自动识别成表格变量列名保留在变量名里时间戳会被读成字符串或datetime后续操作非常顺手。示例代码% 读取测风塔CSV文件 opts detectImportOptions(wind_tower_data.csv); opts.VariableNamingRule preserve; % 保留原始列名 data readtable(wind_tower_data.csv, opts);如果你手里的文件没有表头只有纯数字矩阵那就用readmatrix读取主体数据同时单独把时间戳用readcell读进来拼接。这种方式适合那种导出了两年数据但前几行是采集器配置信息的老式格式。第三种方式适合几十MB以上的大文件用datastore创建数据流按行迭代读取避免一次性把全部数据压进内存导致卡顿。不过测风塔数据一般不会太大——10分钟间隔、全年数据也就52560条记录Excel级别都已经能处理。真正要小心的是那种每秒钟存一次原始样本的采集器数据量直奔百万行才需要考虑分块读取。2.3 时间戳对齐和时区的坑时间戳是整个数据集的锚点。如果时间轴没对齐后面计算平均风速、日变化都会出错。常见的坑有三个第一个坑是时区。采集器默认记录UTC协调世界时但业务分析需要当地时间。中国地区的场址把UTC加8小时即可但要注意很多老式采集器在转换时会把“UTC”写进时间戳导致datetime解析时自带时区偏移最好是读进来以后统一用datetime的TimeZone属性处理再转成无时区的本地时间。第二个坑是夏令时和冬令时。国内不存在这个问题但如果你处理的是境外项目务必用datetime的自动转换特性避免自己写加减小时逻辑否则到了切换日必错。第三个坑是缺失时间戳。采集器偶尔会因为停电、故障错过某个10分钟窗口导致时间轴不连续。导入后第一件事就是检查时间间隔是否恒定time_diff diff(data.timestamp); idx_gap find(time_diff ~ minutes(10));如果有不连续的位置要把它们显式标出来否则后面做傅里叶滤波、日均曲线时会造成假信号。我的习惯是直接在时间轴上生成完整的10分钟序列用retime将原数据对齐到新的时间栅格上缺失位置自然变成NaN处理起来干净又明确。3. 数据预处理清洗决定了分析质量3.1 有效数据完整率是第一道质量门槛做风资源评估第一件事不是急着算平均风速而是先统计这一年数据到底“有效”了多少。行业里有个常用标准月有效数据完整率不低于90%年有效数据完整率不低于90%否则这段时间的统计结果不具备代表性需要补充测风。有效数据指的是“未受仪器故障、覆冰、停机等因素影响”的数据。完整率计算公式很简单valid_ratio sum(~isnan(data.ws_70m)) / height(data) * 100;如果完整率低于90%就需要排查是不是某几个月份仪器坏了。我处理过某项目X的数据发现冬季两个月风速通道大量缺失后来一查是风速计覆冰冻结转不动了这种情况不能简单插值补上因为补出来的数据不代表真实风况只能对缺失月份做标记并在报告中说明。3.2 缺失值处理插值还是剔除缺失值的去留取决于分析目的和缺失模式。如果缺失是零散的、孤立的可以使用线性插值填上对统计结果的影响很小。如果缺失是连续整段比如连续3天没有数据那就直接剔除这段不做任何填充。用Matlab补缺失值有两种推荐方式% 方式一线性插值填充适合零散缺失 data.ws_70m fillmissing(data.ws_70m, linear); % 方式二用前后12小时的中位数平滑填充适合小时级短段缺失 data.ws_70m fillmissing(data.ws_70m, movmedian, hours(24));强调一下千万别对长期缺失段用插值。我曾经见过有人把连续一个月缺失段用线性插值补出来结果算出来的年平均风速比真实情况偏低了将近10%整个发电量估算都被带偏。连续缺失段必须通过完整率扣减的方式处理或者在统计分析里直接屏蔽该类时段。3.3 异常值识别剔除有物理依据的坏点异常值检测不能只靠肉眼扫散点图要有物理依据。我一般按五个规则逐层过滤第一风速不能为负当风速等于0但风向正常时大概率是风速计停转或冻住要剔除。第二风速超出物理上限。虽然理论上强风可以超过60m/s但测风塔风速计的量程一般也就50m/s超出量程的读数必然是失真值直接剔除。第三风向必须在0~360°范围内超出范围的要么是传感器故障要么是记录错误。第四垂直层之间做一致性检查。比如70m风速和50m风速在同一时刻的比值正常情况下在1.0~1.3之间浮动风切变导致高处风速更大如果出现70m风速突然比50m低一大截多半是仪器受到扰动要检查原始记录。第五相邻时间段变化率检查。10分钟内风速突变超过20m/s基本不现实超过的属于异常波动需要标记。这些规则写成Matlab脚本逐条应用最后生成一个质量标记列完整保留原始数据但又能在统计时快速过滤。我一直强调不要在原表上直接删行这样会丢失溯源能力后期复查时很难交代。3.4 风向和风速的相关性校验有一个容易被忽略的细节风向传感器和风速传感器虽然在同一个支架上但它们各自独立可能其中一个坏了另一个还正常。如果发现风向全部指向同一方向或风速正常但风向乱跳先用圆形直方图看一下风向分布是否合理。如果主导风向不在该地区的典型风向上就要怀疑传感器是否装反或者数据记录通道接反。还可以做风速风向的联合相关性检查比如把风速按风向分扇区统计各扇区的平均风速。正常情况下相邻扇区平均风速应该平滑过渡如果某个扇区的平均风速突然离谱地高于周围扇区大概率是局部流场受到塔影效应或者传感器支架干扰这类数据在计算风玫瑰图时需要剔除对应风向区间或者用其他高度层的数据交叉验证后决定是否保留。4. 核心指标计算与统计建模4.1 平均风速、风速日变化与风切变指数平均风速是风资源评估最基本的指标。但直接对所有原始数据求算术平均并不严谨应该以10分钟平均风速为样本先按月按月求平均再按年加权平均这样能消除采样月份不均的影响。风速日变化和季节变化同样重要。白天和夜间的风速可能存在显著差异尤其是在复杂地形区域。用Matlab对多年数据做逐时刻平均可以得到一条“平均日风速曲线”这条曲线对风功率预测很有参考价值。风切变指数描述风速随高度变化的规律。多数情况下风速随高度近似遵循幂律% 计算风切变指数 alpha h1 50; h2 70; v1 nanmean(data.ws_50m); v2 nanmean(data.ws_70m); alpha log(v2/v1) / log(h2/h1);alpha通常在0.1到0.3之间。如果某个高度层的alpha异常偏大或偏小就要对该高度层的数据质量做进一步检查。这个指数在风机轮毂高度外推时非常关键——测风塔不一定刚好建到风机轮毂高度但评估时要利用已有高度层的风速推算出轮毂高度的风资源。4.2 威布尔分布拟合风速频率分布的核心建模风资源评估中最经典的分布模型是双参数威布尔分布。概率密度函数为f(v) (k/c)·(v/c)^(k-1)·exp(-(v/c)^k)其中k是形状参数c是尺度参数。k值越大说明风速分布越集中接近常数风c值近似对应平均风速的水平。拟合威布尔分布有两种常用方法极大似然估计MLE和最小二乘法LS。Matlab里直接调用wblfit函数就能完成MLE拟合% 拟合威布尔分布参数MLE法 ws data.ws_70m(~isnan(data.ws_70m)); % 剔除缺失值 [parmHat, parmCI] wblfit(ws); k parmHat(1); % 形状参数 c parmHat(2); % 尺度参数拟合后要画频率直方图叠加拟合曲线直观评估拟合效果。同时可以算拟合优度比如将风速区间分成若干桶比较每个桶的实际频率和理论频率计算决定系数R²。通常情况下R²大于0.95就算是比较好的拟合。如果拟合效果不理想常见原因是数据里混入了异常时段或风速双峰分布。双峰分布常见于存在明显海陆风交替的地区这时候需要做混合威布尔拟合或者分风向扇区分别拟合不能拿单一分布硬套。4.3 风功率密度真正决定资源等级的量平均风速高不代表资源一定好因为风功率与风速立方成正比风速的波动对能量贡献极大。风功率密度的公式是W/m² 0.5·ρ·v³其中ρ是空气密度取1.225 kg/m³。用实测风速序列逐条计算再平均得到的是实测风功率密度用威布尔参数也能算理论风功率密度W/m² 0.5·ρ·c³·Γ(13/k)其中Γ是伽马函数。对比实测值和理论值可以校验拟合好坏。在Matlab里计算实测风功率密度很直接rho 1.225; % 海平面标准空气密度 power_density 0.5 * rho * data.ws_70m.^3; mean_power_density nanmean(power_density);按国内风资源等级划分50m高度年平均风功率密度大于400W/m²基本具备较好开发条件300~400W/m²属于可开发但经济性一般低于200W/m²基本不适合大规模风电开发。4.4 风向玫瑰图与扇区统计分析风向玫瑰图是风资源报告中最直观的一张图。它把360°划分成16个或12个扇区统计各扇区的风频和平均风速。Matlab里可以用polarhistogram直接绘制% 风向玫瑰图频率分布 figure; polarhistogram(deg2rad(data.wd_70m), 16, FaceColor, blue, FaceAlpha, 0.5);更有价值的是“风向风速联合玫瑰图”每个扇区用不同颜色表示平均风速区间这样能直观看出哪个方向的风既多又强。扇区统计分析还要计算每个扇区的威布尔参数、湍流强度特别是主导风向扇区这些数据直接用于风机排布主导风方向上风机间距要拉大减少尾流影响。4.5 湍流强度与极端风况湍流强度TI定义为风速标准差除以平均风速是风机载荷计算的重要输入% 按扇区计算湍流强度 ws_std movstd(ws, 3); % 3个采样点窗口 ws_mean movmean(ws, 3); TI ws_std ./ ws_mean;通常TI低于0.1说明风况平稳0.1~0.25属于中等湍流超过0.25就属于高湍流区域对风机的疲劳载荷分析要格外谨慎。复杂山地和切变强烈的近海区域TI往往偏高选址时不能只看平均风速。极端风况分析包括年最大风速、极大风速和风向的最大变化率。这些数据用于确定风机的极限载荷工况对塔筒和叶片的结构设计至关重要。从实测数据中提取年最大3s阵风风速和50年一遇极大风速虽然Matlab做不了外推统计但至少可以把实测极值完整提取出来供结构工程师使用。5. 核心Matlab代码实现与逐段解读5.1 完整导入和清洗流程这节给出一个可以直接套用的模板。假设数据文件叫wind_data.csv包含时间戳、70m风速/风向、50m风速/风向、温度、气压六列。% 1. 导入 opts detectImportOptions(wind_data.csv); opts.VariableNamingRule preserve; data readtable(wind_data.csv, opts); % 2. 时间戳解析 data.timestamp datetime(data.timestamp, InputFormat, yyyy-MM-dd HH:mm:ss); % 3. 生成完整时间基准 fullTime (data.timestamp(1):minutes(10):data.timestamp(end)); data retime(table2timetable(data), fullTime, linear);这里用retime把数据对齐到完整时间栅格上缺失位置自动生成NaN。table2timetable是Matlab中推荐处理时间序列的方式后续很多时间相关操作都比普通table方便。5.2 异常值过滤函数把前面说的物理过滤规则封装成一个函数返回逻辑索引方便后续统计function valid filter_anomalies(ws, wd) valid true(size(ws)); valid(ws 0 | ws 60) false; % 风速范围 valid(wd 0 | wd 360) false; % 风向范围 d abs([NaN; diff(ws)]); valid(d 20) false; % 突变检查 end调用时逐高度层过滤生成质量标记列。这里要注意过滤之后的有效数据量如果太少统计结果会显著失真所以第一步先看完整率再决定是否继续。5.3 威布尔MLE拟合及可视化function [k, c, R2] fit_weibull(ws) ws ws(~isnan(ws)); [parmHat, ~] wblfit(ws); k parmHat(1); c parmHat(2); % 拟合优度计算 edges 0:1:ceil(max(ws)); obs histcounts(ws, edges, Normalization, probability); mid edges(1:end-1) 0.5; theo wblpdf(mid, k, c) * 1; % 积分宽度1m/s R2 1 - sum((obs - theo).^2) / sum((obs - mean(obs)).^2); end这里必须注意直方图桶宽的设置。如果桶宽取0.5m/s理论概率也要按对应宽度缩放否则R²算出来会偏小。5.4 风功率密度和扇区统计把扇区统计做成一个表输出每行对应一个扇区列出平均风速、风功率密度、湍流强度、频率占比function stats sector_stats(ws, wd) sectors 0:22.5:360; % 16扇区 stats table(); for i 1:length(sectors)-1 idx wd sectors(i) wd sectors(i1); s_ws ws(idx); stats.mean_ws(i) mean(s_ws, omitnan); stats.power_density(i) 0.5 * 1.225 * mean(s_ws.^3, omitnan); stats.freq(i) sum(idx) / length(ws) * 100; stats.ti(i) std(s_ws, omitnan) / mean(s_ws, omitnan); end end注意扇区边界处理风向正好等于360°的情况要归入第一扇区否则会漏掉正北方向的数据。代码里可以用取模运算处理。5.5 生成项目报告图表输出图一般至少四张风速频率直方图威布尔拟合曲线、风向玫瑰图、风功率密度日变化曲线、风切变指数垂直剖面图。前三张我们上面已经覆盖风切变剖面图实现如下figure; heights [10, 50, 70, 100]; mean_ws_by_height [4.8, 6.9, 7.6, 8.2]; plot(mean_ws_by_height, heights, -o, LineWidth, 1.5); xlabel(平均风速 (m/s)); ylabel(高度 (m)); title(风廓线); % 根据实际需要可改英文 grid on;这张图能直观看出风速随高度的变化趋势配合alpha值一起放到报告里能让审查工程师快速理解场址的风切变特征。6. 常见问题排查与实操心得6.1 时间戳解析报错datetime解析失败最常见的原因是格式不匹配。很多采集器导出的时间戳中带有“T”字符比如2019-01-01T00:00:00这时需要指定InputFormat为yyyy-MM-ddTHH:mm:ss。还有一种情况是导出的时间秒数被省略了yyyy-MM-dd HH:mm也能解析但要小心不同行格式不一致导致的丢数据。排查技巧先unique一下时间戳列查看实际字符格式再写解析规则不要瞎猜。6.2 威布尔拟合不收敛或结果怪异wblfit偶尔对双峰分布数据拟合效果很差表现为k值极小或c值异常偏大。这时候先做风速频率直方图看数据是否真的接近单峰分布。如果明显双峰建议分时段或分扇区拟合如果数据含有大量零风速静风期拟合也会被拉偏可以先剔除风速小于0.5m/s的数据再拟合并在报告里注明剔除比例。6.3 风向玫瑰图方向偏移如果绘制出来的玫瑰图主导风向和已知的区域气候资料明显不符优先检查风向传感器安装是否正对正北。很多塔在安装时由于磁偏角未修正导致整体风向偏差几十度直接用会造成扇区统计全部错位。处理办法是如果已知当地磁偏角对风向序列统一减去该角度后再分析。6.4 完整率超90%但结果不合理这种情况最常见的原因是多个高度层数据来源不同步。比如70m通道是采集器A记录的50m通道是采集器B记录后人工合并的两条时间序列没对齐。务必检查合并后每一行是否同一时间点必要时做插值重采样统一到同一时间基准。6.5 关于实操工具链的一点个人经验坦率说Matlab在风资源分析里最大的优势是“开箱即用”统计工具箱和绘图功能覆盖了90%的日常需求不需要像Python那样安装一堆第三方库。但要注意版本差异——老版本没有timetable和retime接口需要自己写循环处理时间序列。建议尽量用较新版本能省很多工作量。另外建议每个分析过程都保留原始数据副本所有清洗步骤用脚本串联而非手动在变量区操作。我自己的习惯是在脚本开头设置随机种子、保留所有中间变量为结构体这样无论何时回头看都能完整复现当初的分析过程。风资源评估报告的审查非常严格数据溯源能力比计算本身更重要。这块内容其实还可以继续往外延伸比如把威布尔参数用于理论发电量估算、结合SCADA数据做机群尾流分析甚至搭一个简单的风速预测模型。但核心基本功就是这个数据进得来、洗得净、算得准、画得出。这个闭环练熟了再上什么复杂算法都不慌。 
长河君 · 2026/10/19 15:05:41 · 阅读 1.2k 基于Matlab的测风塔数据处理与风能资源评估实战指南
最近在处理某风电场前期选址的测风数据分析工作突然觉得这类活儿挺值得拿出来讲讲。风能资源评估听着是个大词落到具体执行层面其实就是一件事把气象塔行业内常叫测风塔测到的历史风力数据导入Matlab完成清洗、预处理、统计分析最终算出平均风速、威布尔分布参数、风功率密度、主导风向这些核心指标。这篇文章就围绕这个完整链路展开适合刚接触风资源评估的工程师、相关专业的学生以及对新能源数据分析感兴趣的人参考。我会把整个流程拆开讲连踩过的坑一起交代。1. 风能资源评估到底在算什么测风塔的本质是一个架在目标场址上的“长期气象观测站”。塔身上不同高度会安装风速计、风向标、温度探头和气压计一般层高选10m、50m、70m、80m、100m这些典型风机轮毂高度附近。数据采集器通常每10分钟记录一次平均值然后存成时间序列。这个10分钟的平均风速就是后续所有分析的原始素材。为什么必须依赖实测数据而不是气象站数据因为风能资源评估要回答的是“本地风况到底适不适合发电”而不是“这个地区大气候如何”。这两个问题的尺度完全不同。周边气象站的数据虽然年份长、容易拿但海拔、地形粗糙度、山谷走向不一样风况差异非常大。这也是测风塔数据不可替代的原因——它是场址真实风况的采样精度和代表性都是模型再分析数据比不了的。用Matlab来做这件事有它非常实在的好处读入CSV、Excel毫不需要额外写解析模块统计工具箱自带威布尔拟合函数绘图接口能出一整套风资源图表。对比Excel处理几十万行时间序列时的卡顿Python虽然也能做但需要另配一套环境Matlab在这个领域依然是很多设计院和咨询机构的首选工具。整套评估链条有四个核心产出长期平均风速和风切变指数、风速频率分布及威布尔拟合参数、风功率密度、主导风向与湍流强度。这四个指标直接决定了某个风机机型是否适配场址、排布间距怎么拉、理论发电量能到多少。所以后面所有代码都是围绕这四个目标来组织的。2. 气象塔数据的导入与格式处理2.1 先弄清采集器导出的数据长什么样测风塔的数据采集器多数支持直接导出CSV或TXT文件。以我常碰到的某项目X为例导出文件的表头大致长这样时间戳风速70m风向70m风速50m风向50m温度气压2019-01-01 00:007.822486.912526.2923.12019-01-01 00:108.152557.232496.3923.42019-01-01 00:207.572406.442386.1924.0注意到几个细节时间戳一般是字符串格式风速单位是m/s风向单位是度0°表示正北按顺时针增大270°为正西气压单位是hPa或kPa。有些采集器还会附加瞬时风速、最大风速或仪器状态码比如9代表冻结、8代表校验失败这些字段在清洗阶段非常有用千万别一上来就删掉。还有一个容易踩的坑导出文件的编码。很多采集器软件在Windows下导出的是GBK编码Matlab默认用UTF-8读会乱码。处理办法是提前用记事本之类工具另存为UTF-8或者在读取时指定编码参数。2.2 Matlab导入数据的三种方式我最常用的方式是readtable它能把CSV和TXT自动识别成表格变量列名保留在变量名里时间戳会被读成字符串或datetime后续操作非常顺手。示例代码% 读取测风塔CSV文件 opts detectImportOptions(wind_tower_data.csv); opts.VariableNamingRule preserve; % 保留原始列名 data readtable(wind_tower_data.csv, opts);如果你手里的文件没有表头只有纯数字矩阵那就用readmatrix读取主体数据同时单独把时间戳用readcell读进来拼接。这种方式适合那种导出了两年数据但前几行是采集器配置信息的老式格式。第三种方式适合几十MB以上的大文件用datastore创建数据流按行迭代读取避免一次性把全部数据压进内存导致卡顿。不过测风塔数据一般不会太大——10分钟间隔、全年数据也就52560条记录Excel级别都已经能处理。真正要小心的是那种每秒钟存一次原始样本的采集器数据量直奔百万行才需要考虑分块读取。2.3 时间戳对齐和时区的坑时间戳是整个数据集的锚点。如果时间轴没对齐后面计算平均风速、日变化都会出错。常见的坑有三个第一个坑是时区。采集器默认记录UTC协调世界时但业务分析需要当地时间。中国地区的场址把UTC加8小时即可但要注意很多老式采集器在转换时会把“UTC”写进时间戳导致datetime解析时自带时区偏移最好是读进来以后统一用datetime的TimeZone属性处理再转成无时区的本地时间。第二个坑是夏令时和冬令时。国内不存在这个问题但如果你处理的是境外项目务必用datetime的自动转换特性避免自己写加减小时逻辑否则到了切换日必错。第三个坑是缺失时间戳。采集器偶尔会因为停电、故障错过某个10分钟窗口导致时间轴不连续。导入后第一件事就是检查时间间隔是否恒定time_diff diff(data.timestamp); idx_gap find(time_diff ~ minutes(10));如果有不连续的位置要把它们显式标出来否则后面做傅里叶滤波、日均曲线时会造成假信号。我的习惯是直接在时间轴上生成完整的10分钟序列用retime将原数据对齐到新的时间栅格上缺失位置自然变成NaN处理起来干净又明确。3. 数据预处理清洗决定了分析质量3.1 有效数据完整率是第一道质量门槛做风资源评估第一件事不是急着算平均风速而是先统计这一年数据到底“有效”了多少。行业里有个常用标准月有效数据完整率不低于90%年有效数据完整率不低于90%否则这段时间的统计结果不具备代表性需要补充测风。有效数据指的是“未受仪器故障、覆冰、停机等因素影响”的数据。完整率计算公式很简单valid_ratio sum(~isnan(data.ws_70m)) / height(data) * 100;如果完整率低于90%就需要排查是不是某几个月份仪器坏了。我处理过某项目X的数据发现冬季两个月风速通道大量缺失后来一查是风速计覆冰冻结转不动了这种情况不能简单插值补上因为补出来的数据不代表真实风况只能对缺失月份做标记并在报告中说明。3.2 缺失值处理插值还是剔除缺失值的去留取决于分析目的和缺失模式。如果缺失是零散的、孤立的可以使用线性插值填上对统计结果的影响很小。如果缺失是连续整段比如连续3天没有数据那就直接剔除这段不做任何填充。用Matlab补缺失值有两种推荐方式% 方式一线性插值填充适合零散缺失 data.ws_70m fillmissing(data.ws_70m, linear); % 方式二用前后12小时的中位数平滑填充适合小时级短段缺失 data.ws_70m fillmissing(data.ws_70m, movmedian, hours(24));强调一下千万别对长期缺失段用插值。我曾经见过有人把连续一个月缺失段用线性插值补出来结果算出来的年平均风速比真实情况偏低了将近10%整个发电量估算都被带偏。连续缺失段必须通过完整率扣减的方式处理或者在统计分析里直接屏蔽该类时段。3.3 异常值识别剔除有物理依据的坏点异常值检测不能只靠肉眼扫散点图要有物理依据。我一般按五个规则逐层过滤第一风速不能为负当风速等于0但风向正常时大概率是风速计停转或冻住要剔除。第二风速超出物理上限。虽然理论上强风可以超过60m/s但测风塔风速计的量程一般也就50m/s超出量程的读数必然是失真值直接剔除。第三风向必须在0~360°范围内超出范围的要么是传感器故障要么是记录错误。第四垂直层之间做一致性检查。比如70m风速和50m风速在同一时刻的比值正常情况下在1.0~1.3之间浮动风切变导致高处风速更大如果出现70m风速突然比50m低一大截多半是仪器受到扰动要检查原始记录。第五相邻时间段变化率检查。10分钟内风速突变超过20m/s基本不现实超过的属于异常波动需要标记。这些规则写成Matlab脚本逐条应用最后生成一个质量标记列完整保留原始数据但又能在统计时快速过滤。我一直强调不要在原表上直接删行这样会丢失溯源能力后期复查时很难交代。3.4 风向和风速的相关性校验有一个容易被忽略的细节风向传感器和风速传感器虽然在同一个支架上但它们各自独立可能其中一个坏了另一个还正常。如果发现风向全部指向同一方向或风速正常但风向乱跳先用圆形直方图看一下风向分布是否合理。如果主导风向不在该地区的典型风向上就要怀疑传感器是否装反或者数据记录通道接反。还可以做风速风向的联合相关性检查比如把风速按风向分扇区统计各扇区的平均风速。正常情况下相邻扇区平均风速应该平滑过渡如果某个扇区的平均风速突然离谱地高于周围扇区大概率是局部流场受到塔影效应或者传感器支架干扰这类数据在计算风玫瑰图时需要剔除对应风向区间或者用其他高度层的数据交叉验证后决定是否保留。4. 核心指标计算与统计建模4.1 平均风速、风速日变化与风切变指数平均风速是风资源评估最基本的指标。但直接对所有原始数据求算术平均并不严谨应该以10分钟平均风速为样本先按月按月求平均再按年加权平均这样能消除采样月份不均的影响。风速日变化和季节变化同样重要。白天和夜间的风速可能存在显著差异尤其是在复杂地形区域。用Matlab对多年数据做逐时刻平均可以得到一条“平均日风速曲线”这条曲线对风功率预测很有参考价值。风切变指数描述风速随高度变化的规律。多数情况下风速随高度近似遵循幂律% 计算风切变指数 alpha h1 50; h2 70; v1 nanmean(data.ws_50m); v2 nanmean(data.ws_70m); alpha log(v2/v1) / log(h2/h1);alpha通常在0.1到0.3之间。如果某个高度层的alpha异常偏大或偏小就要对该高度层的数据质量做进一步检查。这个指数在风机轮毂高度外推时非常关键——测风塔不一定刚好建到风机轮毂高度但评估时要利用已有高度层的风速推算出轮毂高度的风资源。4.2 威布尔分布拟合风速频率分布的核心建模风资源评估中最经典的分布模型是双参数威布尔分布。概率密度函数为f(v) (k/c)·(v/c)^(k-1)·exp(-(v/c)^k)其中k是形状参数c是尺度参数。k值越大说明风速分布越集中接近常数风c值近似对应平均风速的水平。拟合威布尔分布有两种常用方法极大似然估计MLE和最小二乘法LS。Matlab里直接调用wblfit函数就能完成MLE拟合% 拟合威布尔分布参数MLE法 ws data.ws_70m(~isnan(data.ws_70m)); % 剔除缺失值 [parmHat, parmCI] wblfit(ws); k parmHat(1); % 形状参数 c parmHat(2); % 尺度参数拟合后要画频率直方图叠加拟合曲线直观评估拟合效果。同时可以算拟合优度比如将风速区间分成若干桶比较每个桶的实际频率和理论频率计算决定系数R²。通常情况下R²大于0.95就算是比较好的拟合。如果拟合效果不理想常见原因是数据里混入了异常时段或风速双峰分布。双峰分布常见于存在明显海陆风交替的地区这时候需要做混合威布尔拟合或者分风向扇区分别拟合不能拿单一分布硬套。4.3 风功率密度真正决定资源等级的量平均风速高不代表资源一定好因为风功率与风速立方成正比风速的波动对能量贡献极大。风功率密度的公式是W/m² 0.5·ρ·v³其中ρ是空气密度取1.225 kg/m³。用实测风速序列逐条计算再平均得到的是实测风功率密度用威布尔参数也能算理论风功率密度W/m² 0.5·ρ·c³·Γ(13/k)其中Γ是伽马函数。对比实测值和理论值可以校验拟合好坏。在Matlab里计算实测风功率密度很直接rho 1.225; % 海平面标准空气密度 power_density 0.5 * rho * data.ws_70m.^3; mean_power_density nanmean(power_density);按国内风资源等级划分50m高度年平均风功率密度大于400W/m²基本具备较好开发条件300~400W/m²属于可开发但经济性一般低于200W/m²基本不适合大规模风电开发。4.4 风向玫瑰图与扇区统计分析风向玫瑰图是风资源报告中最直观的一张图。它把360°划分成16个或12个扇区统计各扇区的风频和平均风速。Matlab里可以用polarhistogram直接绘制% 风向玫瑰图频率分布 figure; polarhistogram(deg2rad(data.wd_70m), 16, FaceColor, blue, FaceAlpha, 0.5);更有价值的是“风向风速联合玫瑰图”每个扇区用不同颜色表示平均风速区间这样能直观看出哪个方向的风既多又强。扇区统计分析还要计算每个扇区的威布尔参数、湍流强度特别是主导风向扇区这些数据直接用于风机排布主导风方向上风机间距要拉大减少尾流影响。4.5 湍流强度与极端风况湍流强度TI定义为风速标准差除以平均风速是风机载荷计算的重要输入% 按扇区计算湍流强度 ws_std movstd(ws, 3); % 3个采样点窗口 ws_mean movmean(ws, 3); TI ws_std ./ ws_mean;通常TI低于0.1说明风况平稳0.1~0.25属于中等湍流超过0.25就属于高湍流区域对风机的疲劳载荷分析要格外谨慎。复杂山地和切变强烈的近海区域TI往往偏高选址时不能只看平均风速。极端风况分析包括年最大风速、极大风速和风向的最大变化率。这些数据用于确定风机的极限载荷工况对塔筒和叶片的结构设计至关重要。从实测数据中提取年最大3s阵风风速和50年一遇极大风速虽然Matlab做不了外推统计但至少可以把实测极值完整提取出来供结构工程师使用。5. 核心Matlab代码实现与逐段解读5.1 完整导入和清洗流程这节给出一个可以直接套用的模板。假设数据文件叫wind_data.csv包含时间戳、70m风速/风向、50m风速/风向、温度、气压六列。% 1. 导入 opts detectImportOptions(wind_data.csv); opts.VariableNamingRule preserve; data readtable(wind_data.csv, opts); % 2. 时间戳解析 data.timestamp datetime(data.timestamp, InputFormat, yyyy-MM-dd HH:mm:ss); % 3. 生成完整时间基准 fullTime (data.timestamp(1):minutes(10):data.timestamp(end)); data retime(table2timetable(data), fullTime, linear);这里用retime把数据对齐到完整时间栅格上缺失位置自动生成NaN。table2timetable是Matlab中推荐处理时间序列的方式后续很多时间相关操作都比普通table方便。5.2 异常值过滤函数把前面说的物理过滤规则封装成一个函数返回逻辑索引方便后续统计function valid filter_anomalies(ws, wd) valid true(size(ws)); valid(ws 0 | ws 60) false; % 风速范围 valid(wd 0 | wd 360) false; % 风向范围 d abs([NaN; diff(ws)]); valid(d 20) false; % 突变检查 end调用时逐高度层过滤生成质量标记列。这里要注意过滤之后的有效数据量如果太少统计结果会显著失真所以第一步先看完整率再决定是否继续。5.3 威布尔MLE拟合及可视化function [k, c, R2] fit_weibull(ws) ws ws(~isnan(ws)); [parmHat, ~] wblfit(ws); k parmHat(1); c parmHat(2); % 拟合优度计算 edges 0:1:ceil(max(ws)); obs histcounts(ws, edges, Normalization, probability); mid edges(1:end-1) 0.5; theo wblpdf(mid, k, c) * 1; % 积分宽度1m/s R2 1 - sum((obs - theo).^2) / sum((obs - mean(obs)).^2); end这里必须注意直方图桶宽的设置。如果桶宽取0.5m/s理论概率也要按对应宽度缩放否则R²算出来会偏小。5.4 风功率密度和扇区统计把扇区统计做成一个表输出每行对应一个扇区列出平均风速、风功率密度、湍流强度、频率占比function stats sector_stats(ws, wd) sectors 0:22.5:360; % 16扇区 stats table(); for i 1:length(sectors)-1 idx wd sectors(i) wd sectors(i1); s_ws ws(idx); stats.mean_ws(i) mean(s_ws, omitnan); stats.power_density(i) 0.5 * 1.225 * mean(s_ws.^3, omitnan); stats.freq(i) sum(idx) / length(ws) * 100; stats.ti(i) std(s_ws, omitnan) / mean(s_ws, omitnan); end end注意扇区边界处理风向正好等于360°的情况要归入第一扇区否则会漏掉正北方向的数据。代码里可以用取模运算处理。5.5 生成项目报告图表输出图一般至少四张风速频率直方图威布尔拟合曲线、风向玫瑰图、风功率密度日变化曲线、风切变指数垂直剖面图。前三张我们上面已经覆盖风切变剖面图实现如下figure; heights [10, 50, 70, 100]; mean_ws_by_height [4.8, 6.9, 7.6, 8.2]; plot(mean_ws_by_height, heights, -o, LineWidth, 1.5); xlabel(平均风速 (m/s)); ylabel(高度 (m)); title(风廓线); % 根据实际需要可改英文 grid on;这张图能直观看出风速随高度的变化趋势配合alpha值一起放到报告里能让审查工程师快速理解场址的风切变特征。6. 常见问题排查与实操心得6.1 时间戳解析报错datetime解析失败最常见的原因是格式不匹配。很多采集器导出的时间戳中带有“T”字符比如2019-01-01T00:00:00这时需要指定InputFormat为yyyy-MM-ddTHH:mm:ss。还有一种情况是导出的时间秒数被省略了yyyy-MM-dd HH:mm也能解析但要小心不同行格式不一致导致的丢数据。排查技巧先unique一下时间戳列查看实际字符格式再写解析规则不要瞎猜。6.2 威布尔拟合不收敛或结果怪异wblfit偶尔对双峰分布数据拟合效果很差表现为k值极小或c值异常偏大。这时候先做风速频率直方图看数据是否真的接近单峰分布。如果明显双峰建议分时段或分扇区拟合如果数据含有大量零风速静风期拟合也会被拉偏可以先剔除风速小于0.5m/s的数据再拟合并在报告里注明剔除比例。6.3 风向玫瑰图方向偏移如果绘制出来的玫瑰图主导风向和已知的区域气候资料明显不符优先检查风向传感器安装是否正对正北。很多塔在安装时由于磁偏角未修正导致整体风向偏差几十度直接用会造成扇区统计全部错位。处理办法是如果已知当地磁偏角对风向序列统一减去该角度后再分析。6.4 完整率超90%但结果不合理这种情况最常见的原因是多个高度层数据来源不同步。比如70m通道是采集器A记录的50m通道是采集器B记录后人工合并的两条时间序列没对齐。务必检查合并后每一行是否同一时间点必要时做插值重采样统一到同一时间基准。6.5 关于实操工具链的一点个人经验坦率说Matlab在风资源分析里最大的优势是“开箱即用”统计工具箱和绘图功能覆盖了90%的日常需求不需要像Python那样安装一堆第三方库。但要注意版本差异——老版本没有timetable和retime接口需要自己写循环处理时间序列。建议尽量用较新版本能省很多工作量。另外建议每个分析过程都保留原始数据副本所有清洗步骤用脚本串联而非手动在变量区操作。我自己的习惯是在脚本开头设置随机种子、保留所有中间变量为结构体这样无论何时回头看都能完整复现当初的分析过程。风资源评估报告的审查非常严格数据溯源能力比计算本身更重要。这块内容其实还可以继续往外延伸比如把威布尔参数用于理论发电量估算、结合SCADA数据做机群尾流分析甚至搭一个简单的风速预测模型。但核心基本功就是这个数据进得来、洗得净、算得准、画得出。这个闭环练熟了再上什么复杂算法都不慌。