做风电前期的人一定都有同感一座测风塔真正值钱的不是它站在那儿的姿态而是它在场址上空吹了一年风之后攒下的那几十万个10分钟平均记录。这些数据可以用来回答最核心的问题——这个场址到底有多少风、风从哪来、能不能赚钱。但数据是原始的夹杂噪声、缺测、仪器故障痕迹必须经过一整套导入、清洗、统计、拟合的过程才能变成平均风速、风功率密度、威布尔分布参数这些能写进评估报告的硬指标。这篇文章我就围绕这套完整流程展开用Matlab把从气象塔测量的历史风力数据从excel/csv文件一直处理到图表输出。其中会涉及数据质量控制的几条实用规则、威布尔分布拟合的参数解读、风功率密度的修正算法以及我怎么识别和处理那些让我们头疼的低完整性数据、结冰失真、塔影效应。代码都是可以照着跑的细节里也藏了不少实际项目里才能积累的避坑经验。1. 风能资源评估与气象塔数据分析的整体思路1.1 气象塔测风数据为什么是前期选址的“硬通货”风能资源评估的意义往小了说是判断一个地方值不值得装风机往大了说是给整个风电项目投资决策提供基础依据。而所有评估工作的地基就是测风塔数据。气象塔和气象站不一样气象站测的是“区域气候平均值”数据可能来自城市边缘、机场周边或者某个山头离你实际要装风机的位置可能几十公里远。但测风塔是实实在在立在项目场址中间、甚至立在风机未来要排布的位置附近测得的就是将来风机轮毂高度处实际的风况。行业标准里测风塔一般会在多个高度同时布设风速仪和风向标常见的高度组合是10m、30m、50m、70m有时也会在80m、100m甚至更高处加密测量。高度选择的逻辑很简单风机的轮毂高度大概在什么位置测风塔就必须在那个高度上下都有测量这样才有足够的资料做风切变外推估算轮毂高度处的真实风速。另外一个关键点是时间代表性。测风周期要求至少完整一年这是为了覆盖冬、夏不同季节、不同大气环流形势下的风况变化。如果只测了风季的几个月评估出来的平均风速会比全年真实值明显偏高后面发电量测算就全偏了。所以我在看一套测风数据时第一件事永远是看数据的时间长度和完整性而不是先跑统计。在实际项目里测风塔数据还是后续所有数字化工作的基础输入风机选型需要知道平均风速和湍流强度机位排布优化需要风向玫瑰图发电量估算需要威布尔分布参数和空气密度修正甚至银行评估项目贷款时也要看第三方机构复核过的测风数据结论。测风塔数据的质量真的就是整个项目前期的生命线。1.2 从原始记录到风资源结论的分析链条拿到一套完整的测风塔历史数据我一般会按照下面的链路去处理数据导入把数据采集器导出的csv/txt文件读入Matlab统一时间格式。数据质量控制对风速、风向、温度、气压做范围检查、相关性检查、突变检查把异常值标成缺失。数据统计计算各高度年平均风速、主导风向、风功率密度、湍流强度、风切变指数。分布拟合用威布尔分布拟合风速频率分布得到形状参数和尺度参数。长期订正可选但评估报告里常做将测风期数据与附近长期气象站数据做相关分析把短期测风数据订正到多年平均风况水平。结果输出生成图表和指标表格写入风资源评估报告。这篇文章的重点放在第1到第4步因为它们决定了后续所有结论的底数。第5步长期订正我自己经常觉得是个很微妙的环节——气象塔数据本身准确性高但如果测风期只有8个月那就需要找一个相关性良好的长期参考站来做订正。相关性的好坏直接决定了订正结果的可靠度用Matlab做线性回归的时候我一般会要求相关系数至少到0.7以上才敢去做订正。链条上最容易出问题的其实是第2步数据质量控制。因为数据采集器和传感器在野外工作一整年低温、结冰、沙尘、动物啃咬线缆这些情况都可能让记录出问题。如果一开始没洗干净后面拟合出来的威布尔参数、算出来的风功率密度就都带着污染。2. 数据导入与预处理Matlab面前的准备2.1 你手里的气象塔数据到底是什么格式接触过不同项目后你会发现不同数据采集器导出的格式比想象中更多变。但核心内容其实都差不多逐条记录的时间戳、各高度风速平均值、各高度风向、空气温度、气压有时候还包括风速标准差、最大风速、最小风速。这些记录的时间间隔基本都是10分钟一条对应行业通用的“10分钟平均”标准。举个例子某个项目的数据可能长这样TimestampWS_70mWS_50mWS_30mWD_70mTemp_2mPressure2024-01-01 00:008.327.155.96245.3-8.41012.52024-01-01 00:108.026.955.81243.7-8.11012.2风速单位一般是m/s风向单位是度0~3600为正北温度摄氏度气压hPa或kPa。但强烈建议你拿到数据先确认一遍单位不要默认因为有些采集器输出的风速是m/s有些可能因为标定问题带着偏移量。我就见过一个项目数据表头写的是m/s实际导出的是mph差了一截后面对比才发现不对劲。数据文件大多按月份拆分一个测风塔一年的数据就是12个csv文件。处理时最便捷的方式是把所有月份文件读进来然后拼接成一个长表再用统一的清洗逻辑处理。Matlab的readtable函数在这个场景下非常好用可以自动识别csv表头、数据类型遇到空值会填充为NaN省掉了很多手工解析的麻烦。这里要特别提一个细节时间戳必须转成Matlab的datetime格式后再合并否则直接把12个csv纵向拼接会出现时间顺序错乱。尤其跨年、跨月时那些写成01/01/2024 00:00的文本字符串排序规则和真实时间不一定一致。2.2 用Matlab把原始数据“搬进来”用Matlab导入这类测风塔历史数据最稳的做法是先用readtable逐月读取再统一合并。列名处理上如果csv文件本身有表头且列名规范直接使用。如果表头是多语言或者带单位我建议设置VariableNamingRule,preserve保留原始列名读取后再手动重命名避免中文表头变成V1/V2导致后面代码可读性变差。时间方面如果采集器导出的是“2024-01-01 00:00:00”这种标准格式直接rawData.Timestamp datetime(rawData.Timestamp, InputFormat, yyyy-MM-dd HH:mm:ss);如果遇到“2024/1/1 0:00”这种格式就得给datetime指定对应的InputFormat。最麻烦的是有的采集器会把时间拆成两列一列日期一列时间合起来也不是标准的ISO格式这时候用datetime的InputFormat, yyyy-MM-ddHH:mm:ss这种组合方式处理或者干脆拆开再拼接处理完之后用sortrows按时间排序。读进来的数据我会先看一眼摘要summary(rawData);这一步能快速发现明显的类型问题比如某列被读成了cell数组、时间列还停留文本状态、风速列缺测太多等等。然后再进入清洗。2.3 数据清洗的几条硬原则数据清洗的“硬原则”来自多年的经验教训就是不能看到异常值就随手删掉。我建议按顺序执行下面三类检查第一范围检查。这是最基础的。风速合理范围一般定在0~60m/s之间如果某个测量值小于0或者超过60直接置为缺失。风向正常范围是0~360度温度范围根据项目所在地气候设定比如北方项目可以放宽到-50~50摄氏度气压范围比如950~1080hPa。范围检查的区间不要定得太紧否则会把真实的极端天气数据给误删了。第二突变检查。相邻两条10分钟平均风速之间正常变化不会特别极端。如果突然从3m/s跳到18m/s绝大多数情况下是传感器瞬时故障或者接线松动。我一般用相邻差的阈值做检测风速时变超过10m/s就判定为异常点。这个阈值也可以根据场址是山地还是平原调整——山地环境阵风变化会剧烈一些。第三相关性检查。测风塔不同高度之间的风速序列理论上应该有很强的相关性。如果70m高度风速是8m/s50m高度的风速无论如何不会只有0.2m/s。所以做清洗时可以把不同高度风速序列放在同一个图上看或者直接计算两个高度风速的相关系数。某个时段相关性突然变差往往说明其中一个传感器出了问题。清洗之后的处理原则是能标缺失就标缺失不要随便用相邻值填。缺失值在后续统计里可以用omitnan方式跳过但填进去的是假数据算法再精细也没用。只有在缺失时段非常短比如只有一两个点且后面要进行连续谱分析时才用线性插值补一下插值比例也要在报告里明确写清楚。3. 核心指标计算风速统计、风功率密度与威布尔拟合3.1 平均风速与风功率密度的计算逻辑风资源评估里最常用但最容易误用的指标一个是平均风速一个是风功率密度。平均风速就是所有有效风速记录的算术平均值这个很好理解。但风功率密度不一样它正比于风速的立方。为什么是立方因为空气的动能在单位时间内流过单位面积时功率是0.5ρv³。这带来了一个非常重要的实际影响风速大的那部分记录对平均风功率密度的贡献远远超过风速较小的部分。举个例子假设风速是均匀的每小时都在8m/s平均风速8m/s平均风功率密度大约0.5×1.225×512313W/m²。但如果一天里一半时间风速6m/s一半时间风速10m/s平均风速也是8m/s但平均风功率密度却是0.5×1.225×(2161000)/2372W/m²比均匀风况高了不少。这就是为什么风资源好坏不能用年平均风速单一指标衡量风功率密度才是最贴近实际发电能力的指标。计算风功率密度的时候风速单位要用m/s空气密度单位kg/m³结果单位才是W/m²。常温海平面空气密度大约取1.225kg/m³但实际项目中如果场址海拔高、气温低空气密度会变化。高海拔地区的空气更稀薄同样风速下风机出力会变小所以用实际空气密度修正是评估报告中不可跳过的一步。3.2 威布尔分布拟合Matlab实现与参数解读威布尔分布是描述风速频率分布最常用的统计模型。原因很直观风速的分布通常是右偏的低风速出现频率高高风速出现频率低而威布尔分布正好能拟合这种形状。威布尔分布的概率密度函数里有俩关键参数形状参数k无纲量和尺度参数A单位与风速相同。k决定了这个分布的“形状”k比较小比如1.5~2分布右偏明显说明风速分布离散、常有强风k比较大比如3~4分布更对称、更接近正态说明风速稳定。尺度参数A大致对应平均风速量级但注意它不等于平均风速。平均风速和威布尔参数的关系是平均风速 A × Γ(1 1/k)其中Γ是伽马函数。所以在Matlab中用fitdist(v, wbl)拟合出A和k之后估算平均风速可以直接套这个公式和实际记录的平均风速做交叉验证。如果两者偏差太大说明数据没有较好地被威布尔分布描述这时候就要回头检查数据的质量或分布特征。在风资源评估报告里威布尔分布参数不仅用于描述风速分布还会被进一步用来估算不同风速区间的小时数、计算理论年发电量。风电投资测算时风机的功率曲线会和风速频率分布做卷积卷积结果的好坏很大程度上就取决于这个威布尔拟合的参数是否可靠。所以拟合之前的数据清洗一定要做扎实否则后面都是小数点前的故事。3.3 风向、湍流强度与风切变指数平均风速和风功率密度解决的是“有多少风”的问题风向、湍流强度和风切变解决的是“风机能不能高效带走这些风”的问题。风向分析的核心产物是风向玫瑰图。用Matlab的polarhistogram画36个扇区每个10度的概率分布看得非常清楚。主导风向集中的场址风机排布可以设计成顺着主风向的行列机位之间尾流影响会小一些如果风向很分散排布设计就要更加均匀。风向数据每年、每季都可能有差异所以最好按季度分开画风季和冬季的主导风向可能会差出几十度这对机位排布影响很大。湍流强度TI等于风速标准差除以平均风速。它衡量的是风速在10分钟内的脉动程度。湍流强度高的地方风机叶片承受的交变载荷大疲劳寿命受影响选型时就需要选择高湍流等级的机型或者适当降低功率曲线控制策略来保护机组。行业里一般把TI按10分钟序列计算而且通常只统计风速大于3m/s以上的数据因为低风速段的湍流强度没太大工程意义。风切变指数是用不同高度风速推导出来的。基本的指数公式为风速比 高度比的α次方所以α ln(v2/v1) / ln(h2/h1)。切变指数大说明地面上风速随高度增长得快这对高塔筒、大叶轮风机更有利切变指数小甚至出现负值说明近地面风速混合强低层风能相对丰富但高层的增益不明显。切变指数也会随昼夜和季节变化白天地面受热对流强切变小夜间大气层结稳定切变可能变大。这些细节在选塔筒高度和做发电量估算时都必须考虑进去。4. 完整流程代码从导入到出图的落地实现4.1 数据导入与清洗代码这里给出一套可以直接落地跑的Matlab实现。假设你的测风塔数据已经按月导出成csv文件名类似met_tower_202401.csv每个文件的列结构一致字段包括时间戳Timestamp、70m风速WS_70m、50m风速WS_50m、30m风速WS_30m、70m风向WD_70m、2m温度Temp_2m、气压Pressure。%% 气象塔测风数据导入与清洗 clc; clear; close all; % 读取一个月的实测数据文件可按需循环读入多个月 rawData readtable(met_tower_202401.csv, VariableNamingRule, preserve); % 时间列转为datetime格式 rawData.Timestamp datetime(rawData.Timestamp, InputFormat, yyyy-MM-dd HH:mm:ss); % 提取常用列统一为单独变量方便后期调用 ws70 rawData.WS_70m; ws50 rawData.WS_50m; ws30 rawData.WS_30m; wd70 rawData.WD_70m; temp rawData.Temp_2m; pres rawData.Pressure; %% 范围检查 % 风速合理范围0~60 m/s ws70(ws70 0 | ws70 60) NaN; ws50(ws50 0 | ws50 60) NaN; ws30(ws30 0 | ws30 60) NaN; % 风向范围0~360° wd70(wd70 0 | wd70 360) NaN; % 温度范围-40~50℃根据项目地调整 temp(temp -40 | temp 50) NaN; %% 突变检查 % 相邻10分钟风速差超过10 m/s判定为异常 delta abs(diff(ws70)); idxSpike [false; delta 10]; ws70(idxSpike) NaN; %% 完整率统计 comp70 sum(~isnan(ws70)) / length(ws70) * 100; fprintf(70m风速数据完整率: %.2f%%\n, comp70);这段代码里有一个细节值得留意突变检查只针对70m高度做了示例实际项目中我会对每个高度都跑一遍同时还会检查不同高度之间的风速序列相关性偏差比如某个时间段内70m和50m风速差超过一定阈值那就说明至少一个传感器数据不可信需要进一步人工核验。4.2 威布尔拟合与风功率密度计算代码导入并完成数据清洗后就可以计算风资源评估的核心指标了。下面的代码先做威布尔分布拟合然后计算风功率密度并加入了空气密度修正。%% 风速分布分析与指标计算 % 提取有效风速剔除缺失和零风速保证威布尔拟合要求v0 vPos ws70(~isnan(ws70) ws70 0); % 威布尔分布拟合 pd fitdist(vPos, wbl); A pd.A; % 尺度参数 k pd.B; % 形状参数 fprintf(威布尔参数: 尺度A%.3f, 形状k%.3f\n, A, k); % 由威布尔参数估算平均风速与实际平均风速交叉验证 meanVFit A * gamma(1 1/k); meanVActual mean(vPos); fprintf(威布尔估算平均风速: %.2f m/s, 实际平均风速: %.2f m/s\n, meanVFit, meanVActual); %% 风功率密度计算标准空气密度 rhoStd 1.225; % kg/m3 wpdStd 0.5 * rhoStd * mean(vPos.^3); fprintf(标准空气密度下平均风功率密度: %.2f W/m2\n, wpdStd); %% 空气密度修正 % 依据实测平均温度和气压修正空气密度 meanTempK mean(temp, omitnan) 273.15; meanPresPa mean(pres, omitnan) * 100; % 假设气压单位hPa rhoActual meanPresPa / (287.05 * meanTempK); wpdActual 0.5 * rhoActual * mean(vPos.^3); fprintf(实际空气密度: %.4f kg/m3\n, rhoActual); fprintf(修正后平均风功率密度: %.2f W/m2\n, wpdActual);使用fitdist(vPos, wbl)拟合时有一个新手经常踩的坑fitdist要求数据全部大于0风速等于0的点如果不剔除拟合会直接报错或出现异常结果。所以我对数据做了一步ws70 0筛选剔除零风速。这样处理在工程上是合理的因为威布尔分布本就是描述正风速分布的模型。另一个常见的争论点是年平均风速要不要包含零风速时段从工程统计角度我们更关心有风时段的风能资源水平零风速静风时段在机制模型里信息量有限剔除后计算的指标对风机选型更有指导意义。4.3 风向玫瑰图、切变指数与湍流强度代码第三个代码块把风向、切变指数和湍流强度一次算明白同时输出图表。%% 风向分析、切变指数与湍流强度 % 风向玫瑰图 dirValid wd70(~isnan(wd70)); figure(Color, w); polarhistogram(deg2rad(dirValid), 36, Normalization, probability, FaceColor, [0.1 0.4 0.8]); title(70m高度风向玫瑰图); set(gca, ThetaZeroLocation, top); % 正北朝上 % 求得主风向扇区出现频率最高的风向区间 [N, edges] histcounts(dirValid(:), 0:10:360); [~, idx] max(N); dominantDir (edges(idx) edges(idx1)) / 2; fprintf(主导风向扇区中心角度: %.1f°\n, dominantDir); %% 风切变指数计算 % 70m与50m高度之间的切变指数 alpha log(ws70 ./ ws50) ./ log(70 / 50); alphaValid alpha(isfinite(alpha) ~isnan(alpha)); alphaMed median(alphaValid); fprintf(70m-50m中位切变指数: %.3f\n, alphaMed); % 如果数据表中有10分钟风速标准差列则计算湍流强度 if ismember(Std_70m, rawData.Properties.VariableNames) std70 rawData.Std_70m; ti std70 ./ ws70; % 通常只统计风速大于3m/s时的湍流强度更有工程意义 idxTI ~isnan(ti) isfinite(ti) ws70 3; medianTI median(ti(idxTI)); fprintf(70m高度中位湍流强度: %.3f\n, medianTI); end这段代码里polarhistogram是我最喜欢用的一个函数36个扇区基本能把风向分布特征清楚地展示出来。有个细节值得说明默认情况下polarhistogram的0度是在右侧但对风资源分析来说正北朝上是约定俗成的习惯所以我加了set(gca, ThetaZeroLocation, top)设置。切变指数的计算里我只取了有限且不为NaN的比值。两个高度风速都非常小的时候比值会不稳定计算出的切变指数可能偏离物理实际。所以遇到比值异常大或异常小比如超出-0.5~1.5区间的值在统计中位值时应先剔除。5. 常见问题与排查技巧实录5.1 数据完整率低于80%怎么办这是实际项目中最常见也最棘手的问题。测风塔运行一年因为故障、维护、存储异常数据完整率低于90%很正常但如果跌到80%以下评估结论的可靠性就要打问号了。如果完整率不足第一步先看缺失的时段集中在哪个季节。如果缺的是冬季风最大的那两个月那年平均风速会明显偏低如果缺的是夏季风小的时段平均风速又会偏高。缺失时段的风况特征直接决定了偏差方向。第二步要判断缺失是随机零散缺失还是整段长时间缺失。零散缺失相对好办用保留的数据做统计并在报告中说明缺失时段比例与代表性问题整段长时间缺失比如连续一个月没有数据就比较麻烦需要把这一个月用临近年份同期的气象塔记录或者附近长期气象站的相关关系来补足。第三步是补数据的操作手法。我不建议直接拿平均值去填缺失段那样等于人为降低方差。更合理的做法是用附近参考站的同期数据通过回归方程把目标时段的风速估算出来再以插值替换缺失。但这里有个原则插值得到的数据只作为辅助参考不能参与到威布尔拟合和风功率密度统计里否则会平滑掉真实的高风速分布让评估结果偏保守。5.2 结冰导致的“零风速长尾巴”怎么识别在北方、高海拔场址冬季测风塔结冰是常态。风杯结冰后转不起来记录的风速会连续保持在接近零的数值甚至连续几天都是0.01、0.02这样的微小值。这类数据如果混入统计会让冬春季平均风速被严重拉低。识别这种问题的关键线索有两个一是气温如果当时温度低于0℃且风速长时间只有0~0.5m/s很大程度上就是结冰导致的虚假静风二是同时刻不同高度的风速值如果50m高度正常、70m高度却非常接近零基本可以说70m传感器已经冻住了。处理上首先是搞清楚测风塔的传感器是否带加热功能。带加热的传感器结冰概率低但也不能完全避免。冬季低温时段的数据最好结合温度序列做辅助判断当温度低于某个阈值比如-5℃且风速持续低于正常下限时就该把这段记录整体标记为无效。另一个实践经验是对比结冰时段前后一两个星期同一时刻的风速如果存在系统性偏低就说明这段时间数据不可用。有的从业者会用“零风速比例超过多少就判定传感器结冰”的规则其实不太稳妥。真实静风也是存在的尤其在山谷、盆地地形里。最有效的还是温度结合多高度相关性一起判断宁可漏标也不要错杀因为错杀会丢掉真实的低风速分布信息。5.3 测风塔影子效应与传感器故障判断测风塔本身是结构件大风从塔体背面吹过时会在塔后形成一个低速尾流区。安装在同一高度的风速仪如果位于这个尾流区里测到的风速会明显偏低这就是所谓的“塔影效应”。判断塔影效应最直接的办法是把风速按风向扇区做统计。比如塔体朝向某个方向那风向在该方向附近时风速比其他风向扇区显著偏低尤其是偏低的扇区范围稳定出现在一个方向区间时基本可以断定存在塔影影响。按行业处理惯例这些受影响扇区的数据会在统计时剔除或者在报告中注明修正方法。传感器故障又是另一类问题。风速计的轴承磨损、信号线接触不良都会导致数据出现大幅跳变或持续异常。我的排查习惯是每次处理新数据集都用plot把风速序列整体画出来用眼睛扫一遍全局。图上那些突兀的“毛刺”或者“台阶”就是需要重点标记的时段。这一步看图的功夫省不掉自动化清洗只能处理常规逻辑真正的异常往往需要人眼判断。风速仪方位偏移也是高频问题。风向标被强风打偏后风向数据会整体偏移一个固定角度。判断方法是看风向玫瑰图是否出现明显不合理的单峰或者对比两个高度风向传感器之间差值分布是否长期偏离某一固定值。一旦发现偏移就要在数据处理阶段做方位校正否则主导风向分析就会出错直接影响机位排布。这套流程走下来我的体会是Matlab的代码实现并不是最复杂的地方真正花时间的是把数据从“原始记录”变成“可信样本”的过程。威布尔参数也好风功率密度也好都是建立在高质量数据之上的上层建筑。数据质量差再漂亮的拟合曲线也只是纸上谈兵。所以我每次处理测风塔数据都会留出专门的时间做可视化检查和质量记录把每一步清洗操作的逻辑和剔除量都写清楚——这不仅是对项目负责也是对自己后面写报告、解释异常时省事。如果你手头正握着一套当地气象塔的历史测风数据可以先从最简单的范围检查和完整率统计开始一点一点把分析流水线搭起来。等你跑通一次完整流程后面再做类似项目就会快很多。最后分享一个个人习惯处理完一批测风数据后我会把关键指标和当时的天气过程记录放在一起保存下来下次再遇到异常数据时翻一翻这些备注往往比重新分析更快找到原因。