简介本资源为风能资源评估方向的Matlab代码工具包面向计算机、电子信息工程、数学等专业的大学生及科研入门者可用于课程设计、期末大作业与毕业设计。包内共25个文件以19个m脚本为核心辅以mat、txt、accdb与kml等数据与地理文件压缩包约5.99MB覆盖风速分布分析、风功率密度计算、风能玫瑰图绘制、数据质量控制与长期风数据抓取等模块。代码采用参数化编程关键参数可灵活调整注释清晰并附赠案例数据可直接运行便于快速搭建评估模型并验证结果。已有66人学习适合需要完整实现思路与可复用脚本的读者参考借鉴。1. 风能资源评估的 MATLAB 实现从测风数据到可投资结论拿到一份测风塔的原始数据直接算个平均风速就敢写资源评估报告这事我干过后来被现实狠狠教育。风能资源评估远不止“平均风速”四个字它要回答的是这个点位一年能发多少电、风机的选型是否匹配、湍流和切变会不会让机组提前报废、50 年一遇的极端风速能不能扛住。MATLAB 在这套流程里的价值是把 Weibull 拟合、风功率密度、风频玫瑰图、湍流强度、极大风速推算这些环节串成一条可复现的流水线而不是每次换一个测风塔就重写一遍脚本。这篇文章面向的是手里有测风数据、需要出评估结论的风资源工程师和新能源方向的研究生我会把每个环节的代码、参数含义和踩过的坑都摊开讲让你拿到数据就能跑通。2. 测风数据的预处理与质量控制脏数据不清理后面全是玄学风资源评估的第一步永远不是建模而是把测风塔导出的原始数据洗干净。测风塔通常按 10 分钟或 1 小时记录字段包括各高度层的风速、风向、温度、气压还有一堆状态码。原始数据里常见的脏东西有传感器结冰导致的持续零风速、通讯中断造成的整段缺失、风向在 0/360 度附近的跳变、以及超出物理量程的野值。这些不处理后面 Weibull 拟合出来的 k 和 c 参数会偏得离谱年发电量估算能差出 20% 以上。2.1 读入与字段对齐我一般先把数据统一成 MATLAB 的 timetable时间列做严格校验因为不同测风塔导出的时间戳格式五花八门有的是本地时间有的是 UTC混用会让风频玫瑰图整体旋转。% 读入测风塔导出的 CSV假设列顺序为时间, 10m风速, 30m风速, 50m风速, 风向, 温度 opts detectImportOptions(tower_raw.csv); opts.VariableNamingRule preserve; T readtimetable(tower_raw.csv, opts); % 统一时间列为 datetime 并设为行时间 T.Properties.DimensionNames{1} Time; T sortrows(T, Time); % 按时间排序防止乱序 T rmmissing(T, DataVariables, ... {WindSpeed_10m,WindSpeed_30m,WindSpeed_50m}); % 先删关键列缺失 % 物理量程过滤风速合理范围 0~40 m/s超出视为野值 validIdx T.WindSpeed_50m 0 T.WindSpeed_50m 40; T T(validIdx, :);这段代码的关键点有三个。readtimetable保证时间列被识别为 datetime 而不是字符串后续做时间重采样才不会出错。sortrows是血泪经验有些数据采集软件在断电恢复后会把缓存数据追加到文件末尾时间顺序是乱的。物理量程过滤用 40 m/s 作为上限是因为绝大多数陆上测风塔的杯式风速仪量程就在这个范围超过的基本是电磁干扰产生的野值。2.2 缺失值插补与结冰标记缺失值不能简单删掉否则会破坏时间序列的连续性影响后续湍流强度的计算。常见做法是连续缺失小于 3 个点的用线性插值大于 3 个点的整段标记为无效不参与统计。结冰的识别靠温度和风速的联合判断——温度低于 2 摄氏度且风速持续为 0大概率是传感器结冰而不是真的无风。% 线性插值填补短缺失最多 3 个连续点 T.WindSpeed_50m fillmissing(T.WindSpeed_50m, linear, MaxGap, 3); % 结冰标记温度 2 且风速 0标记为无效 iceFlag T.Temperature 2 T.WindSpeed_50m 0; T.WindSpeed_50m(iceFlag) NaN; % 统计有效数据完整率 validRate sum(~isnan(T.WindSpeed_50m)) / height(T); fprintf(50m 风速有效数据完整率%.2f%%\n, validRate * 100);MaxGap参数设成 3 是经验值对应 10 分钟数据里的半小时缺口。完整率低于 90% 的测风塔我一般会建议延长观测周期或者换塔因为 IEC 61400-12 对数据完整率有明确要求低于这个阈值做出来的年发电量估算置信度不够。结冰标记这一步很多人会漏结果冬季月份的平均风速被拉低Weibull 分布的尺度参数 c 偏小最后算出来的发电量偏保守投资决策就会误判。2.3 风向数据的环形处理风向是环形数据0 度和 360 度是同一个方向。做风频玫瑰图时如果直接按数值分箱0 度附近的数据会被拆到首尾两个箱子里玫瑰图会出现一个假的缺口。正确做法是用circ_系列函数或者手动做角度折叠。% 风向折叠到 [0, 360) T.WindDirection mod(T.WindDirection, 360); % 按 16 方位分箱每个方位 22.5 度 edges 0:22.5:360; binIdx discretize(T.WindDirection, edges, IncludedEdge, left); binIdx(isnan(binIdx)) 16; % 360 度归入第 16 个方位 % 统计各方位风频 freq histcounts(binIdx, 1:17) / sum(~isnan(binIdx));discretize的IncludedEdge参数设成left是为了让每个区间的左边界闭合避免 0 度被重复计数。最后那个binIdx(isnan(binIdx)) 16是处理 360 度整点的情况因为discretize对等于最后一个边界的值会返回 NaN。这些细节看起来琐碎但风频玫瑰图是给业主看的第一张图画错了会直接影响对主导风向的判断进而影响风机排布方案。3. Weibull 分布拟合与风功率密度计算参数选错发电量差一个档风资源评估的核心是把风速的统计分布拟合出来最常用的是两参数 Weibull 分布。形状参数 k 反映风速的集中程度尺度参数 c 反映风速的整体大小。这两个参数决定了风功率密度而风功率密度直接决定了这个点位值不值得开发。很多人直接用wblfit一把梭但测风数据里有效样本量、分箱方式、拟合方法的选择都会影响结果下面把几种常见做法和适用场景讲清楚。3.1 三种拟合方法的对比与选择MATLAB 里拟合 Weibull 分布至少有三种路径最大似然估计wblfit、矩估计、以及最小二乘法拟合累积分布。最大似然估计在样本量大于 1000 时最稳定矩估计计算快但对野值敏感最小二乘法适合数据已经分箱成风频直方图的情况。% 方法一最大似然估计推荐样本量充足时 pd fitdist(T.WindSpeed_50m(~isnan(T.WindSpeed_50m)), Weibull); k_mle pd.A; % 形状参数 c_mle pd.B; % 尺度参数 % 方法二矩估计快速估算用于初步筛查 v_mean mean(T.WindSpeed_50m, omitnan); v_std std(T.WindSpeed_50m, omitnan); k_mom (v_std / v_mean)^(-1.086); c_mom v_mean / gamma(1 1/k_mom); % 方法三最小二乘拟合累积分布数据已分箱时 [sortedV, idx] sort(T.WindSpeed_50m(~isnan(T.WindSpeed_50m))); empCDF (1:length(sortedV)) / length(sortedV); X log(sortedV); Y log(-log(1 - empCDF)); p polyfit(X, Y, 1); k_lsq p(1); c_lsq exp(-p(2) / p(1));三种方法的结果通常不会完全一致k 值差 0.1、c 值差 0.3 是常有的事。我的习惯是先用矩估计快速看一眼再用最大似然估计出正式结果最后用最小二乘法交叉验证。如果三种方法算出来的风功率密度差异超过 5%说明数据本身有问题得回头查预处理环节。fitdist返回的pd.A是形状参数pd.B是尺度参数这个命名和 Weibull 分布的标准写法是反的第一次用很容易搞混。3.2 风功率密度与可开发等级判定风功率密度是评估点位好坏的核心指标计算公式里用到了 gamma 函数。根据国家标准风功率密度等级从 1 级到 7 级3 级以上才具备较好的开发价值。% 空气密度修正按实测温度和气压 rho mean(T.Pressure ./ (287.05 * (T.Temperature 273.15)), omitnan); % 风功率密度计算 WPD 0.5 * rho * c_mle^3 * gamma(1 3/k_mle); % 年发电量初步估算假设风机扫风面积 A功率系数 Cp A 12000; % 扫风面积单位 m^2对应约 4MW 机组 Cp 0.45; % 功率系数现代机组典型值 hours 8760; E_annual WPD * A * Cp * hours / 1000; % 单位 kWh fprintf(空气密度%.3f kg/m^3\n, rho); fprintf(风功率密度%.1f W/m^2\n, WPD); fprintf(年发电量估算%.0f kWh\n, E_annual);空气密度这一项很多人直接取 1.225 kg/m³但在高原或者高温地区实际空气密度可能只有 1.0 甚至更低风功率密度会差出 20%。gamma(1 3/k_mle)这一项是 Weibull 分布三阶矩的解析解不能省。年发电量估算里的Cp取 0.45 是当前主流机组的典型值但实际应该查具体机型的功率曲线这里只是初步筛查用。如果算出来的风功率密度低于 150 W/m²我一般会直接建议放弃这个点位除非有特殊的上网电价政策。3.3 风频玫瑰图与主导风向分析风频玫瑰图不只是画着好看它决定了风机排布的主轴方向。主导风向对应的风能方向机组应该尽量垂直于这个方向排列减少尾流影响。% 计算各风向方位的风频和平均风速 dirEdges 0:22.5:360; dirIdx discretize(T.WindDirection, dirEdges, IncludedEdge, left); dirIdx(isnan(dirIdx)) 16; freqByDir zeros(1, 16); meanSpdByDir zeros(1, 16); for i 1:16 mask dirIdx i; freqByDir(i) sum(mask) / sum(~isnan(dirIdx)); meanSpdByDir(i) mean(T.WindSpeed_50m(mask), omitnan); end % 绘制风频玫瑰图 theta deg2rad(dirEdges(1:16) 11.25); figure; polarplot([theta, theta(1)], [freqByDir, freqByDir(1)], b-o, LineWidth, 1.5); title(50m 高度风频玫瑰图);polarplot画玫瑰图时角度需要取每个方位区间的中心值所以加了 11.25 度的偏移。风频和平均风速要分开统计因为主导风向不一定等于最大风能方向——有时候某个方向风频不高但风速很大风能贡献反而更大。这个区别在排布方案里很关键尾流计算要用风能玫瑰图而不是风频玫瑰图。4. 湍流强度与极端风速推算决定机组选型和塔筒寿命风资源评估报告里如果只有平均风速和风功率密度那是不完整的。湍流强度和极端风速直接关系到机组的安全等级和塔筒的疲劳寿命。湍流强度超标机组会频繁触发振动保护发电量损失不说关键部件寿命也会大幅缩短。极端风速算错塔筒设计载荷不够后果更严重。4.1 湍流强度的计算与分级湍流强度定义为风速标准差与平均风速的比值通常按 10 分钟窗口计算再取 90% 分位数作为代表值。IEC 61400-1 把湍流强度分为 A、B、C 三个等级对应不同的设计标准。% 按 10 分钟窗口计算湍流强度 T10 retime(T, regular, mean, TimeStep, minutes(10)); T10.TI T10.WindSpeed_50m_std ./ T10.WindSpeed_50m; % 取 90% 分位数作为代表湍流强度 TI_rep prctile(T10.TI(~isnan(T10.TI)), 90); % 按风速分箱统计湍流强度IEC 要求 15 m/s 处的值 binEdges 0:2:30; spdBin discretize(T10.WindSpeed_50m, binEdges); TI_by_bin zeros(1, length(binEdges)-1); for i 1:length(binEdges)-1 mask spdBin i; if sum(mask) 10 TI_by_bin(i) prctile(T10.TI(mask), 90); end endretime做 10 分钟重采样时标准差需要单独算不能直接用mean。我一般会写一个自定义函数传给retime同时输出均值和标准差。prctile取 90% 分位数是 IEC 规定的做法不能用平均值代替因为湍流强度的分布是右偏的平均值会低估极端情况。如果 15 m/s 处的湍流强度超过 0.16这个点位就需要选 IEC A 类机组机型选择范围会窄很多造价也更高。4.2 极大风速的极值推算极端风速推算用极值理论常用 Gumbel 分布拟合年最大风速序列。测风塔通常只有 1 到 2 年的数据直接拟合年最大风速样本量不够常见做法是用月度最大值或者用风速的日最大值序列做拟合。% 提取各月最大风速 T_monthly retime(T, monthly, (x) max(x, [], omitnan)); monthlyMax T_monthly.WindSpeed_50m; % Gumbel 分布拟合 pdGumbel fitdist(monthlyMax(~isnan(monthlyMax)), ExtremeValue); mu_g pdGumbel.mu; sigma_g pdGumbel.sigma; % 推算 50 年一遇极大风速 T_return 50; V_50 mu_g - sigma_g * log(-log(1 - 1/T_return)); fprintf(50 年一遇极大风速50m%.1f m/s\n, V_50);Gumbel 分布的mu是位置参数sigma是尺度参数。log(-log(1 - 1/T_return))是 Gumbel 分布的分位数函数T_return是重现期。这里用月度最大值而不是年最大值是因为测风塔观测周期通常不够长月度最大值能把样本量扩大 12 倍。但要注意月度最大值序列不满足极值理论的独立性假设算出来的结果会偏保守实际工程中会再乘一个 0.9 左右的折减系数。这个折减系数没有统一标准各家咨询公司用的值不一样报告里要写清楚。4.3 风切变与轮毂高度风速外推测风塔通常只测到 50m 或 70m而现代机组的轮毂高度在 100m 以上需要做风切变外推。常用的是幂律模型风切变指数 alpha 由两个高度的风速比取对数得到。% 计算风切变指数用 30m 和 50m 数据 mask ~isnan(T.WindSpeed_30m) ~isnan(T.WindSpeed_50m) ... T.WindSpeed_30m 2 T.WindSpeed_50m 2; alpha log(mean(T.WindSpeed_50m(mask)) / mean(T.WindSpeed_30m(mask))) / ... log(50 / 30); % 外推到 100m 轮毂高度 V_hub mean(T.WindSpeed_50m, omitnan) * (100 / 50)^alpha; fprintf(风切变指数%.3f\n, alpha); fprintf(100m 轮毂高度平均风速%.2f m/s\n, V_hub);风切变指数 alpha 的典型值在 0.1 到 0.3 之间小于 0.1 说明大气层结不稳定大于 0.3 说明可能有逆温层。mask里限制风速大于 2 m/s 是为了排除静风时段因为低风速下的风向和风速测量误差都很大会污染 alpha 的计算。外推到 100m 时如果 alpha 是用 30m 和 50m 算的外推距离比较远误差会放大最好用 50m 和 70m 的数据算。如果测风塔只有一层数据那风切变只能取经验值报告里要注明这个不确定性。5. 避坑与排查那些让评估报告翻车的细节风资源评估这活理论都不难难的是数据里藏着的各种坑。下面这几条是我这些年踩过的、也见过别人踩的每条都按现象、原因、解决来写你对着排查能省不少时间。5.1 现象Weibull 拟合的 k 值大于 3风功率密度异常偏高原因通常是数据里混入了大量低风速的无效记录比如传感器故障时的持续零风速或者静风时段的重复计数。这些低风速值会把分布的形状参数拉大c 值也跟着偏大。解决方法是回到预处理环节检查零风速的占比如果超过 20%要逐条排查是真实静风还是传感器问题。另外用fitdist之前先画一下风速直方图看看分布形状是不是双峰双峰说明数据里混了不同来源的记录。5.2 现象风频玫瑰图的主导风向和附近测风塔完全相反原因多半是风向传感器的安装方向错了或者数据导出时做了镜像翻转。有些测风塔的风向传感器在安装时没有对准正北而是按塔筒的某个标记对齐导致整体偏移一个固定角度。解决方法是找同一区域其他测风塔的数据做交叉验证如果偏移是固定的可以整体旋转修正。如果偏移不固定那可能是传感器故障这段数据不能用。5.3 现象年发电量估算比实际运行数据高 30% 以上原因通常是尾流模型没考虑或者空气密度用了标准值而实际点位海拔较高。尾流损失在大型风电场里能到 10% 到 20%单台机组的估算不能直接乘以机组数量。解决方法是先用 Jensen 模型或者更精细的 CFD 做尾流计算再乘一个综合折减系数包括可利用率、电气损耗、叶片污染等。空气密度必须用实测温度和气压算不能偷懒取 1.225。5.4 现象极端风速推算结果每次跑都不一样原因是极值拟合对样本量很敏感月度最大值序列里如果有一两个异常大的值拟合结果就会跳。解决方法是先做异常值检验用 Grubbs 检验或者箱线图法把离群点挑出来单独看。如果离群点是真实的极端天气事件要保留如果是传感器野值要剔除。另外重现期推算的置信区间也要算出来只给一个点估计值是不严谨的。5.5 现象MATLAB 脚本在换了一台电脑后报错原因通常是路径分隔符或者编码问题。Windows 和 Linux 的路径分隔符不一样中文路径在某些 MATLAB 版本里也会出问题。解决方法是所有路径用fullfile拼接不要手写斜杠。数据文件尽量用英文命名CSV 的编码用 UTF-8 无 BOM。如果用了自定义函数确保函数文件在 MATLAB 的搜索路径里或者用addpath显式添加。6. 从评估到报告把 MATLAB 结果变成可交付的结论跑完上面所有计算你手里会有一堆数字k、c、风功率密度、湍流强度、极大风速、年发电量。但这些数字本身不是报告报告要回答的是“这个点位能不能投、投多少、风险在哪”。我一般会把 MATLAB 的输出整理成三张核心表风资源参数汇总表、风频风能分布表、以及不同机型方案的年发电量对比表。下面这段代码把关键结果导出成结构化的表格方便后续写报告时直接引用。% 汇总关键参数 results table( ... {形状参数 k; 尺度参数 c; 平均风速 (m/s); 风功率密度 (W/m^2); ... 湍流强度 (90%); 50年一遇极大风速 (m/s); 风切变指数; 年发电量估算 (kWh)}, ... [k_mle; c_mle; mean(T.WindSpeed_50m,omitnan); WPD; TI_rep; V_50; alpha; E_annual], ... VariableNames, {Parameter, Value}); % 导出为 CSV供报告引用 writetable(results, wind_resource_summary.csv); % 不同机型方案对比假设三种轮毂高度和功率曲线 hubHeights [90, 100, 110]; ratedPower [3.6, 4.0, 4.5]; % MW for i 1:length(hubHeights) V_h mean(T.WindSpeed_50m,omitnan) * (hubHeights(i)/50)^alpha; % 简化估算实际应查功率曲线插值 CF 0.35; % 容量系数需根据功率曲线和风速分布计算 AEP ratedPower(i) * 1000 * CF * 8760; fprintf(轮毂高度 %dm额定功率 %.1fMW估算年发电量 %.0f kWh\n, ... hubHeights(i), ratedPower(i), AEP); endresults表里的每一项都要在报告里给出计算方法和不确定度。比如风功率密度的不确定度主要来自风速测量误差和空气密度估算误差一般给 ±10%。年发电量的不确定度更大能到 ±15%。这些不确定度不是摆设投资决策里的内部收益率计算要用到。机型对比那段代码里的CF是简化处理实际应该用功率曲线和 Weibull 分布做卷积积分但那个计算量大初步筛查用容量系数就够了。最后说一个我自己的习惯每次跑完评估我都会把原始数据、预处理后的数据、以及所有中间结果存成不同的文件命名里带上日期和版本号。风资源评估经常要反复迭代业主会提各种新要求比如换个机型、改个轮毂高度、加个尾流方案。如果中间结果没存每次都要从头跑一遍浪费时间不说还容易在重复计算中引入不一致。MATLAB 的save命令配合-v7.3参数能存大文件timetable 也能直接存恢复的时候用load就行。这个习惯看起来笨但关键时刻能救命。希望帮到你。本文还有配套的精品资源点击获取
