做风资源评估的人手里最重要的是什么很多刚从学校出来的朋友会觉得是WAsP或者WindPRO这类高级软件但真正落到“我手里只有一座70米测风塔一年的10分钟平均数据要怎么判断这片区域值不值得继续往下做”的时候最踏实的第一步反而是用Matlab把数据老老实实读进来清洗干净算出一批基础参数画几张能说服自己的图。这篇文章就用一个典型的气象塔历史风力数据案例从导入、处理到计算与出图完整走一遍Matlab流程。无论你是刚入行的评估工程师、准备做毕业设计的研究生还是临时接手测风数据项目的开发者都可以照着这套思路把数据摸透。为了避免光贴代码我会把每个步骤背后的原理也讲清楚。你不需要事先掌握很多统计学知识跟着操作就能理解为什么平均风向不能用普通mean为什么风功率密度要用风速的三次方为什么Weibull拟合对估算年发电量这么重要。1. 风能资源评估的核心指标与气象塔数据基础1.1 气象塔测风数据里到底有什么测风塔不是简单的立一根杆不同高度会安装不同传感器。常见测风高度是10m、30m、50m、70m再高的塔还会测到90m、100m或更高具体看项目需求。每个高度通常有一台风速计和一个风向标塔身上还会配置温度、气压和湿度传感器。更专业一点的测风塔会加装超声风速仪用来记录高频湍流数据。对常规风资源评估来说气象塔输出的数据一般是10分钟平均记录其中包含该时段平均风速、最大风速、最小风速以及标准差。如果是1Hz采样再聚合成10min均值原始文件里可能还有每秒原始数据但这类文件体积很大通常只保留统计值。你看一个典型数据表列名大概长这样TIMESTAMP、WD10、WS10、WD30、WS30、TEMP70、PRESS70。时间列是采样时刻WD和WS分别代表风向和风速后面的数字是测量高度。但真正打开这种文件之后坑特别多。时间戳格式不统一、个别行列错位、风向标半夜冻住、风速传感器被扬尘堵住这些问题我在项目里都遇到过。所以数据处理的第一步不是急着算数而是先搞清楚传感器类型、采样原理、数据约定再动手清洗。否则后面所有计算结果都建立在一个不牢靠的基础上。1.2 风能资源评估的关键指标风资源评估看起来指标很多真正核心的可以压缩成下面这六项。第一是平均风速。通常指某个高度上10min平均风速的算术平均是最常用的资源描述。但平均风速不能单独使用因为它没有反映风速分布。第二是风功率密度。这是风能资源丰歉的核心指标公式是0.5乘以空气密度再乘以风速三次方的平均。因为风机捕获的功率与风速三次方成正比所以风功率密度比平均风速更能说明问题。很多初学者会用平均风速直接立方这是错的正确做法是对每个样本的风速三次方取平均再乘以0.5和空气密度。第三是Weibull分布参数。风速频率分布通常用双参数Weibull分布描述形状参数k和尺度参数c。k决定了分布的形态和风速的集中程度c则反映整体风速量级。机组发电量估算和极端风计算都离不开这两个参数。第四是湍流强度等于风速标准差除以平均风速公式是TIσ/v。湍流强度影响机组疲劳载荷和实际出力在复杂地形尤其重要。第五是主导风向。包括风向频率分布和风能密度方向分布这决定了机位怎么排布、机组之间的尾流影响有多大。第六是风切变指数。风速随高度变化通常用幂律公式v2/v1(h2/h1)^α来推算不同高度的风速α就是风切变指数。它影响从测风高度向轮毂高度折算风速的可靠性。1.3 数据的时间分辨率与质量控制前提评估一个场址测风数据至少要覆盖一个完整年因为风资源的季节性非常强。夏天和冬天的平均风速可能相差一倍以上。如果项目只给了你三个月的数据顶多能做初步判断不能直接拍板年发电量。时间分辨率上10min平均是行业最常用的基线。1h平均会把阵风信息平滑掉算出来的湍流强度会偏低1s原始数据则用于湍流和极端风分析文件量巨大。常规风资源评估我建议至少保留10min平均数据再做小时和日聚合分析。质量控制必须放在所有统计计算之前不能跳过。数据有效率低于90%不代表不能使用关键要看缺失分布。如果一年里只有某一个月集中缺测哪怕总有效率95%也不能简单插补后当作完整年计算。野点剔除也需要规则比如风速超出0~60m/s、风向超出0~360°、一个小时内风速突变超过阈值的点都要打标签。气象塔数据永远不可能完美但只要你把清洗规则记录下来评估结论就可追溯可复现。2. Matlab数据导入与预处理实操2.1 不同格式数据的导入方法比较说回Matlab。我常用的导入方式是根据文件大小和格式来选。小规模Excel文件readtable最省事几十万行的CSVreadmatrix或textscan更快如果文件自带可解析的时间列readtimetable能一步到位。先用readtable作为通用入门方案代码是这样的% 读取CSV保留原始列名 raw readtable(met_tower.csv, PreserveVariableNames, true); % 把列名统一成后续方便使用的名称 % 假定原始列顺序时间、风向、风速、温度、气压 raw.Properties.VariableNames {Time,WindDir,WindSpeed,Temp,Press}; % 解析时间戳常见格式 yyyy-MM-dd HH:mm:ss raw.Time datetime(raw.Time, InputFormat, yyyy-MM-dd HH:mm:ss); % 转成timetable后面很多函数都围绕timetable设计 tt table2timetable(raw); % 如果你有多个高度比如WS10和WS50用cell数组统一列名也可以 % raw.Properties.VariableNames {Time,WindDir,WindSpeed,Temp,Press,WS10,WS50};如果文件时间列本身能被Matlab自动识别直接用readtimetable更简单tt readtimetable(met_tower.csv);但readtimetable对时间格式比较敏感时间戳如果带有时区后缀“08:00”或者日期直接用“20230101”这种整数表示它不一定能解析。此时还是要走readtable加datetime转换的老路。Excel文件也不难处理多个sheet分别对应不同高度时最好用readtable分别读取再按时间戳做outerjoin合并。我见过一个项目把70m风速放在一个sheet温度放在另一个sheet日期还不对齐最后只能先分别读入用时间戳join把所有列拼在一起。这种时候千万不要凭行号对齐时间戳才是唯一基准。2.2 缺失值处理插补还是剔除导入数据后第一件事是看缺失情况。用summary可以快速知道每列有多少missing值也可以自己统计% 统计每列缺失数量 missCount sum(ismissing(tt)); disp(missCount);缺失处理的原则是“零散缺失可插补连续成片要剔除”。温度和气压存在明显的日周期零散缺失用线性插补影响不大。但某几个连续小时的10min风速缺口如果直接用线性插补填上相当于人为制造了一段并不存在的持续风会污染湍流强度计算。这种情况应该保留NaN后续计算时用omitnan跳过。我自己的经验阈值是连续缺失超过6个10min点也就是1小时就不插补了。代码可以这么做% 先对零散缺失做线性插补 tt.WindSpeed fillmissing(tt.WindSpeed, linear); % 找出连续缺失超过6个点的缺口重新置为NaN missVec double(ismissing(tt.WindSpeed)); d diff([0; missVec(:); 0]); starts find(d 1); ends find(d -1) - 1; runLengths ends - starts 1; % 只把长缺失段设为NaN longRuns runLengths 6; for i find(longRuns) tt.WindSpeed(starts(i):ends(i)) NaN; end注意fillmissing默认是线性插补它会完全用前后相邻点连一条直线。如果缺口只有两三个点这种处理几乎不影响统计量。但如果你做湍流和极端阵风分析连1个点的插补都要谨慎最好保留原始标记。2.3 异常值与野点剔除异常值不能只靠“看着不对”来剔除要有可解释的规则。第一道规则是物理合理性判断% 风速合理范围单位m/s validWS tt.WindSpeed 0 tt.WindSpeed 60; % 风向合理范围单位° validWD tt.WindDir 0 tt.WindDir 360; % 温度合理范围单位°C validTemp tt.Temp -50 tt.Temp 60; % 组合成综合有效标记 valid validWS validWD validTemp; tt.WindSpeed(~valid) NaN;第二道规则是针对突变点。常用的两种方法滑动中位数和变化率阈值。滑动中位数适合捉孤立的尖峰变化率阈值适合抓传感器跳变。% 用1小时滑动中位数判断离群点 out isoutlier(tt.WindSpeed, movmedian, hours(1)); tt.WindSpeed(out) NaN; % 10min平均风速前后变化超过10m/s属于异常跳变 dV abs([NaN; diff(tt.WindSpeed)]); tt.WindSpeed(dV 10) NaN;这里有个很容易犯的错不要直接把大于“mean3*std”的值全删掉。风速本身是右偏分布强风过程完全属于正常现象3σ会把正常的极端大风误删。物理范围和变化率才是更安全的过滤器。风向传感器还有一个典型故障冻结。温度低于0°C时风向标可能被冰卡住输出长时间保持同一个角度。处理方法是检查每个风向值连续重复的次数超过一定数量就标记为异常。风速和气温过低可以作为辅助判断条件。这类问题只靠自动算法不一定抓得准最好再叠加人工抽查时间序列。3. 风能资源参数计算与Matlab代码实现3.1 平均风速与风向的统计预处理完成后先做最基础的统计。Matlab里的retime函数可以把10min数据聚合成小时平均或日平均在分析季节性变化时非常方便% 小时平均风速 ws_hourly retime(tt(:, WindSpeed), hourly, mean); % 日平均风速 ws_daily retime(tt(:, WindSpeed), daily, mean); % 年平均风速 ws_yearly retime(tt(:, WindSpeed), yearly, mean);风向就不能这么算了。风向是角度0°和360°本身是同一个方向。如果你直接对350和10取平均得到的是180°完全错误。正确做法是把风向当成圆上的向量对sin和cos分别取平均再通过atan2还原角度% 风向圆平均 u mean(sin(deg2rad(tt.WindDir)), omitnan); v mean(cos(deg2rad(tt.WindDir)), omitnan); dirMean rad2deg(atan2(u, v)); if dirMean 0 dirMean dirMean 360; end这样做350和10的平均就会接近0°而不是180°。除了平均方向还可以用合成向量长度R来量化风向稳定性R接近1说明风向非常集中R接近0说明风向散乱。这个R值在报告里作为风向稳定性指标也很有用。3.2 风功率密度与空气密度修正风功率密度是整个评估里最不能算错的一个数。空气密度ρ并不是常数它随气温、气压和湿度变化。标准海平面条件下大约是1.225kg/m³但在海拔3000m的高原密度可能只有0.9kg/m³直接影响风功率密度和发电量。常规做法是用温度和气压反算密度% 假设温度单位是°C气压单位是kPa Tk tt.Temp 273.15; Ppa tt.Press * 1000; rho Ppa ./ (287.05 .* Tk); % 风功率密度逐个样本计算后取平均单位W/m² wpd mean(0.5 .* rho .* tt.WindSpeed.^3, omitnan);注意这里用的是逐样本的rho乘以v³再取平均不是用平均rho乘以平均风速的三次方。这两个结果差异很大因为风速立方放大了大风时段的贡献。我在项目里见过有人把平均风速代入v³算出来的风功率密度低估了20%以上。风功率密度对高风速段极其敏感数据里只要有几个强风过程就能明显抬升整体值如果用平均风速替代恰恰把最关键的贡献丢掉了。3.3 Weibull分布参数估计与拟合质量判断风速频率分布通常用双参数Weibull分布描述。Matlab里的fitdist可以直接做极大似然估计% 取出有效风速忽略NaN ws_valid rmmissing(tt.WindSpeed); % Weibull分布拟合 pd fitdist(ws_valid, Weibull); c pd.A; % 尺度参数 k pd.B; % 形状参数如果用的Matlab版本比较老可以用wblfit返回参数顺序也是尺度参数在前形状参数在后别搞反了。另一个我常用的快速估计方法是矩估计mu mean(ws_valid); sigma std(ws_valid); k0 (sigma/mu)^(-1.086); c0 mu / gamma(1 1/k0);这个公式在绝大多数情况下给出的初值都很稳也可以用来核对fitdist的结果。k通常在1.5到3.5之间c大约等于平均风速的1.05到1.15倍。如果你拟合出来的k只有0.8或者高达6就要警惕数据分布是不是被异常值污染了。参数算出来不是终点。要把频率直方图和拟合曲线叠加重点看大风速段是否拟合合理。如果拟合曲线在高风速端明显低于柱状图说明尾段被低估年发电量可能会被低估反之高估则要检查野点。做这种判断时最好配合年发电量计算一起看单纯看曲线贴合程度是不够全面的。3.4 风玫瑰图与主导风向识别风玫瑰图是判断主导风向最直观的工具。Matlab里的polarhistogram可以快速画出风向频率图dirRad deg2rad(tt.WindDir); figure; polaraxes(ThetaZeroLocation, top, ThetaDir, clockwise); polarhistogram(dirRad, 16, Normalization, probability); title(风向频率玫瑰图);这段代码把北放在正上方角度按顺时针方向递增符合国内气象报告的习惯。16个扇区每个跨22.5°是风资源报告里最常见的扇区划分方式。风能玫瑰比频率玫瑰更实用因为它反映了哪个方向传来的风真正携带更多的能量。polarhistogram不好直接支持加权概率我通常在Matlab里算扇区统计表再决定画图方式edges 0:22.5:360; dirBin discretize(tt.WindDir, edges); ok ~isnan(dirBin) ~isnan(tt.WindSpeed); % 每个扇区风向频率 freq accumarray(dirBin(ok), 1, [16 1]) / sum(ok) * 100; % 每个扇区平均风速 meanWS accumarray(dirBin(ok), tt.WindSpeed(ok), [16 1], nanmean, NaN); % 每个扇区风能密度 E 0.5 .* rho .* tt.WindSpeed.^3; meanE accumarray(dirBin(ok), E(ok), [16 1], nanmean, NaN); pctE 100 .* meanE ./ nansum(meanE); T table(edges(1:end-1) 11.25, freq, meanWS, pctE, ... VariableNames, {扇区中心方向, 频率百分比, 平均风速, 风能占比}); disp(T);这张统计表本身就可以写进报告。要注意的是扇区边界不要重复0°和360°其实同一个方向discretize如果直接用0:22.5:360作为边界360°会被分到最后一个扇区容易多算。稳妥的做法是把所有等于360°的风向先改成0再参与分扇区。3.5 湍流强度与风切变湍流强度是风速标准差除以平均风速。如果你的气象塔数据文件已经包含每10min的风速标准差列直接相除即可% 假设有标准偏差列 WS_std TI tt.WS_std ./ tt.WindSpeed; % 风速越低TI越不稳定一般只统计风速大于4m/s时的湍流强度 sel tt.WindSpeed 4; TI_mean mean(TI(sel), omitnan);风速接近0时TI会极端偏大这种情况下统计出来的平均值没有工程意义所以要加一个风速下限。如果你的数据没有标准差列只能用原始高频数据算或者用小时标准差代替但结果会偏小报告里必须注明。风切变指数用两高度风速计算。假设有10m和50m风速按下式求αv10 mean(tt.WS10, omitnan); v50 mean(tt.WS50, omitnan); alpha log(v50 / v10) / log(50 / 10);切变指数通常在0.1到0.5之间。白天太阳加热、大气混合充分α接近0.1夜间地面冷却、大气层结稳定α可能超过0.3。如果项目关心轮毂高度风速最好分季节和昼夜统计α用同一个年均值推算可能产生明显偏差。切变指数还和地形、地表粗糙度密切相关复杂地形区域更要用实测数据统计不能随便套公式。4. 结果可视化和评估报告输出4.1 时间序列图与高度对比评估报告里图表比表格更能直接说明问题。最基础的时间序列图要解决两件事整体趋势和细节突变。拿一整年10min数据直接画一条曲线会密成一片黑什么也看不清。所以我通常先选前30天看细节再用月平均或日平均图看长期趋势。% 前30天大约4320个10min点 figure; plot(tt.Time(1:4320), tt.WindSpeed(1:4320), LineWidth, 0.8); hold on; plot(tt.Time(1:4320), tt.WS50(1:4320), LineWidth, 0.8); legend(70m风速, 50m风速, Location, best); xlabel(时间); ylabel(风速m/s); grid on;全年数据改成月平均柱状图更清楚一眼能看出哪几个月风大、哪几个月风小。箱线图也能看每个月的风速分布范围% 提取月份 tt.Month month(tt.Time); % 绘制月风速箱线图 figure; boxchart(tt.Month, tt.WindSpeed); xlabel(月份); ylabel(风速m/s);如果你的Matlab版本没有boxchart可以用传统的boxplot代替效果类似。重点在于报告里除了“年平均风速”还要有季节变化的描述。4.2 风速频率直方图与Weibull叠加频率直方图和Weibull拟合曲线叠加是风资源报告里最常出现的图也是最能看出拟合问题的图。前面已经算出了c和k画图这一步很直接ws_valid rmmissing(tt.WindSpeed); pd fitdist(ws_valid, Weibull); c pd.A; k pd.B; figure; histogram(ws_valid, Normalization, pdf, NumBins, 40, ... FaceColor, [0.7 0.7 0.7]); hold on; vline linspace(0.5, max(ws_valid), 100); plot(vline, wblpdf(vline, c, k), r-, LineWidth, 2); legend(实测风速概率密度, Weibull拟合); xlabel(10min平均风速m/s); ylabel(概率密度); grid on;观察曲线时重点看高风速尾段。如果实测柱状最右边还有不少超过20m/s或30m/s的点而拟合曲线已经衰减到接近0说明Weibull分布没有完全抓住强风过程。这时候不要急着手动改参数先检查这些大风是不是异常值确认无误后可以在报告中单独指出强风过程对评估的影响必要时通过极值分布模型补充描述。4.3 风玫瑰图的美化与导出风玫瑰图画完后导出图片时要注意分辨率和尺寸。Matlab里exportgraphics是一个很稳定的选择figure; polaraxes(ThetaZeroLocation, top, ThetaDir, clockwise); polarhistogram(deg2rad(tt.WindDir), 16, Normalization, probability); title(风向频率玫瑰图); % 导出300dpi的PNG exportgraphics(gcf, wind_rose.png, Resolution, 300);风能玫瑰的绘制稍微麻烦一点。polarhistogram不支持权重我通常把扇区统计表算好后在Matlab里用polaraxes加自定义扇形填充或者直接导出表后用其他绘图工具画。如果你愿意折腾polaraxes底层的cartesian坐标转换能实现扇形填充但代码量不小。更建议的做法是先输出上文中的T表格再根据表格画图这样统计与绘图解耦数据核对也更容易。4.4 一键生成评估表格最后把关键参数汇总成一个表方便写报告也方便项目交接。不要在分析结束后临时拼凑从导入那天起就建立汇总变量每算完一个指标就写入一行后面引用会非常方便。result table; result.meanWS mean(tt.WindSpeed, omitnan); result.maxWS max(tt.WindSpeed, [], omitnan); result.wpd wpd; result.TI_mean TI_mean; result.k k; result.c c; result.dirMean dirMean; writetable(result, wind_resource_summary.csv);如果项目有多个测风高度把高度信息也加到表里比如多一行“测风高度70m”。这样汇总表本身就是一个小型数据库。导出图片时统一用exportgraphics设置300dpi保存为PNG插入Word或PDF报告都很清晰。5. 实操中的常见问题与避坑经验5.1 塔影效应和传感器故障识别塔影是测风塔支架和传感器之间相互作用产生的测量偏差。风从塔支架下游方向吹来时风速会被扭曲典型特征是某个固定扇区的平均风速明显低于相邻扇区。这个问题不画风玫瑰很难被发现。处理方法是识别受扰扇区并剔除该扇区数据。常见的受扰扇区是塔体正对来流的下风向区域范围大约30°到60°。具体做法是绘制不同扇区的平均风速图找持续低值扇区再对照现场照片或塔结构图确认受扰方向最后在统计计算时把该扇区数据标记为无效。有些评估软件内部也有扇区有效率概念但如果你在用Matlab处理必须自己完成这个过程。5.2 Weibull拟合的收敛性问题风速数据如果大量集中在0附近fitdist可能出现Warning甚至拟合出非常不合理的k值。这是因为Weibull分布的对数似然函数在风速等于0时会出现极端值数据里的静风点越多影响越大。一个稳的办法是自定义负对数似然函数用mle优化并设置参数边界和初值ws_valid rmmissing(tt.WindSpeed); ws_valid ws_valid(ws_valid 0.1); % 过滤静风点但记录数量变化 nlogLL (params) -sum(log(wblpdf(ws_valid, params(1), params(2)))); param0 [c0, k0]; opts statset(MaxIter, 5000); phat mle(ws_valid, pdf, wblpdf, start, param0, ... Options, opts, LowerBound, [0.01, 0.2], ... UpperBound, [100, 10]); c phat(1); k phat(2);过滤掉静风点会影响发电量估算所以要在报告中注明过滤比例特别是静风占比高的地区。更稳妥的做法是把静风时段单独统计不参与Weibull拟合但后续发电量计算时再把静风小时数加回来。5.3 风向圆统计的边界陷阱前面提过风向平均不能用普通mean。这里再强调一个容易踩的坑当定义扇区边界时0°和360°是同一个方向但程序会把它们当成两个数。如果你用0:22.5:360作为边界discretize会把0分到第一个扇区360分到最后一个扇区相当于最后多算了一次北极向。解决办法很简单先把所有风向修正到[0, 360)区间wd mod(tt.WindDir, 360); wd(wd 0) 360; % 如果希望以360代表正北也可以但要保证边界统一实际处理时要固定一种约定并在脚本注释里写清楚否则下一个人接手时很难看懂。这种细节最容易在多人协作时引入错误所以代码注释不要省。5.4 数据代表性不足怎么办如果测风塔数据只有八个月不能直接把评估结果当作长期资源结论。常规做法是找附近长年代气象站或再分析数据的同期风速序列和测风塔数据做线性回归再把测风塔数据订正到长周期平均风速。这个流程叫MCP本质上是把短期测量“挂靠”到长期参考序列上。MCP听起来高级实际做起来有几个前提。参考站距离不能太远地形特征不能差太多参考序列质量也要检查。再分析数据用来补空间覆盖可以但直接替代现场测站风险很大因为再分析数据的分辨率不一定能反映局地地形影响。所以拿到一份历史测风数据时先不要急着写代码先问清楚数据的时间范围、原始采样频率、维护记录和故障标记。这些背景信息比任何算法都重要。数据文件本身不会告诉你传感器坏过只有维护日志会。5.5 报告里的“平均风速”陷阱很多评估报告只给年平均风速这远远不够。平均风速相同的两个场址风功率密度可能差30%因为风功率密度对高风速段极其敏感。正确做法是同时给出年平均风速、风功率密度、Weibull参数、湍流强度、主风向扇区、静风小时占比。另外季节变化一定要展示。单一个年均值会掩盖“冬季大风、夏季小风”或“春秋稳定、夏冬波动”这些重要特征。报告评审时一张月平均风速表和一张风向频率表往往比一段文字更有说服力。最后再分享一个我自己早期的教训有一次我在算湍流强度时直接用全套数据求标准差和均值结果12月晚上几个夜间强风过程把TI平均值拉得很高。后来把数据按月、按昼夜分组发现真正异常的是那几天传感器结冰。从那之后我养成了“先分组观察、再整体统计”的习惯也建议你们都保留一个清洗日志。气象塔数据不可能完美但只要你把每一步清洗规则和剔除量记录下来评估结论就是可追溯、可复现的。这比任何高级算法都重要。