简介这是一份面向风能行业工程师、科研人员及高校学生的MATLAB小程序用于对风电场风速数据进行两参数威布尔分布拟合与分析。资源包内仅含1个.m脚本文件压缩包体积约436B文件小巧但功能聚焦可直接在MATLAB环境中运行完成风速样本的读取、参数估计与分布拟合。脚本围绕形状因子k与尺度参数λ展开涵盖数据预处理、极大似然法或矩匹配法参数估计、概率密度函数构建、风速直方图与威布尔曲线对比绘制以及平均风速、标准差、湍流强度等特征值计算最终输出估计参数与统计指标为风电场功率预测、风险评估和设备选型提供量化依据。目前已有6233人学习下载说明该工具在风速统计建模场景中具有较高的实用参考价值适合希望快速上手威布尔分布拟合、理解风速概率特性的读者借鉴与二次开发。1. 风电场风速建模的起点两参数威布尔分布到底能算什么如果你手头有一整年的测风塔数据领导让你评估这片风场的风能资源你大概率绕不开一个东西——风速的概率分布。风速不是均匀出现的它有时静风、有时狂风直接拿平均风速去算发电量误差能大到让你怀疑人生。两参数威布尔分布就是干这个的用形状参数 k 和尺度参数 c 两个数字把一整年的风速概率结构压缩成一个可计算的数学表达。这个 MATLAB 小程序解决的就是从实测风速序列出发拟合出 k 和 c再画出概率密度曲线和累积分布曲线顺带算出平均风速、风功率密度这些风电场前期评估里绕不开的指标。适合谁用做风资源评估的工程师、学可再生能源方向的学生、以及任何需要快速验证风速分布拟合效果的开发者。它不依赖任何工具箱纯 MATLAB 基础函数就能跑拿到数据改个路径就能出图。2. 项目结构与核心算法从数据到曲线的完整链路2.1 目录结构与文件职责这个项目结构非常扁平没有复杂的包管理典型的小工具型 MATLAB 项目。拿到源码后你会看到类似这样的组织方式weibull_wind/ ├── main.m % 主入口串联数据读取、拟合、绘图 ├── data/ │ └── wind_speed.csv % 实测风速序列单列单位 m/s ├── funcs/ │ ├── weibull_fit.m % 两参数威布尔拟合核心函数 │ ├── weibull_pdf.m % 概率密度函数计算 │ └── weibull_cdf.m % 累积分布函数计算 └── output/ └── figures/ % 生成的图片存放目录main.m是唯一需要手动运行的脚本其余都是函数文件。data目录放原始数据格式要求很宽松单列 CSV 即可第一行可以是表头也可以不是代码里用readmatrix读取遇到非数值行会自动跳过。output/figures目录如果不存在main.m里会先判断再创建避免保存图片时报错。提示如果你的风速数据是 Excel 格式把readmatrix换成readtable再取.Variables即可不需要改其他逻辑。2.2 两参数威布尔拟合的实现拟合的核心思路是用最大似然估计MLE求解形状参数 k 和尺度参数 c。风速序列记作 $v_1, v_2, ..., v_n$两参数威布尔的概率密度函数为$$f(v) \frac{k}{c} \left(\frac{v}{c}\right)^{k-1} \exp\left[-\left(\frac{v}{c}\right)^k\right]$$对数似然函数对 k 求导后没有解析解所以代码里用fminsearch做数值优化。下面是weibull_fit.m的关键实现function [k, c] weibull_fit(v) % v: 风速序列列向量单位 m/s % k: 形状参数无量纲 % c: 尺度参数单位 m/s v v(:); % 强制列向量 v(v 0) []; % 剔除静风或负值威布尔要求 v 0 n length(v); % 用矩估计给初值避免 fminsearch 陷入局部最优 mu mean(v); sigma std(v); k0 (sigma / mu)^(-1.086); % 经验公式 c0 mu / gamma(1 1/k0); % gamma 函数来自 MATLAB 内置 % 负对数似然函数优化目标是最小化它 negloglik (params) -sum( log(params(1)/params(2)) ... (params(1)-1)*log(v/params(2)) - (v/params(2)).^params(1) ); options optimset(Display, off, TolX, 1e-6); params_hat fminsearch(negloglik, [k0, c0], options); k params_hat(1); c params_hat(2); end逻辑说明先剔除零值和负值因为威布尔分布的定义域是正实数。初值用矩估计法给出这一步很关键——如果直接给[1, 1]这种拍脑袋初值fminsearch有可能收敛到局部极值拟合出来的曲线和直方图对不上。负对数似然函数里每一项都做了向量化没有循环风速序列上万条也能秒出结果。TolX设为1e-6是精度和速度的折中再小意义不大再大曲线会明显偏离。参数说明输入v是实测风速单位必须是 m/s如果你拿到的数据是 km/h先除以 3.6。输出k一般在 1.5 到 3.0 之间风电场常见值是 2.0 左右c通常比平均风速略大比如平均风速 7 m/s 的风场c大概在 7.8 到 8.2 之间。2.3 概率密度与累积分布的计算拟合出 k 和 c 之后画图需要计算理论曲线。weibull_pdf.m和weibull_cdf.m分别对应概率密度和累积分布function y weibull_pdf(v, k, c) % 两参数威布尔概率密度函数 y (k/c) * (v/c).^(k-1) .* exp(-(v/c).^k); end function y weibull_cdf(v, k, c) % 两参数威布尔累积分布函数 y 1 - exp(-(v/c).^k); end这两个函数都做了向量化传入风速向量直接返回对应向量。weibull_cdf在算风功率密度的时候会用到因为风功率密度正比于风速立方需要用到 Gamma 函数代码里用gamma(1 3/k)直接算不需要数值积分。2.4 主流程与绘图main.m把上面这些串起来核心步骤就四步读数据、拟合、算曲线、画图。下面是关键片段%% 读取数据 data readmatrix(data/wind_speed.csv); v data(~isnan(data)); % 去掉缺失值 %% 拟合 [k, c] weibull_fit(v); %% 计算理论曲线 v_grid linspace(0, max(v)*1.2, 200); pdf_theory weibull_pdf(v_grid, k, c); cdf_theory weibull_cdf(v_grid, k, c); %% 绘图 figure(Position, [100, 100, 900, 400]); subplot(1,2,1); histogram(v, 30, Normalization, pdf, FaceColor, [0.7 0.8 0.9]); hold on; plot(v_grid, pdf_theory, r-, LineWidth, 2); xlabel(风速 (m/s)); ylabel(概率密度); legend(实测直方图, sprintf(Weibull k%.2f, c%.2f, k, c)); title(风速概率密度分布); grid on; subplot(1,2,2); [F, x] ecdf(v); plot(x, F, b-, LineWidth, 1.5); hold on; plot(v_grid, cdf_theory, r--, LineWidth, 2); xlabel(风速 (m/s)); ylabel(累积概率); legend(实测累积分布, Weibull 拟合); title(风速累积分布); grid on; %% 输出统计量 v_mean c * gamma(1 1/k); v_std c * sqrt(gamma(1 2/k) - gamma(1 1/k)^2); fprintf(形状参数 k %.4f\n, k); fprintf(尺度参数 c %.4f m/s\n, c); fprintf(拟合平均风速 %.4f m/s\n, v_mean); fprintf(实测平均风速 %.4f m/s\n, mean(v));逻辑说明左图把实测直方图和理论 PDF 叠在一起肉眼就能判断拟合好坏右图用ecdf画实测累积分布再叠理论 CDF如果两条线基本重合说明威布尔假设成立。最后输出拟合平均风速和实测平均风速做对比偏差超过 5% 就要回头检查数据里有没有异常值。参数说明histogram的Normalization设为pdf是为了和理论曲线量纲一致不设的话纵轴是频数曲线会压得很低。linspace的上限取max(v)*1.2是为了让曲线尾部完整显示取太小会截断取太大右边全是空白。3. 运行环境与实操步骤从零跑通第一张图3.1 MATLAB 版本与依赖检查这个项目对 MATLAB 版本要求不高R2018b 及以上都能跑。用到的函数里readmatrix是 R2019a 引入的如果你用的是更早版本换成csvread或xlsread就行。histogram是 R2014b 引入的再早的版本用hist代替但归一化参数要手动算。不需要任何工具箱fminsearch、gamma、ecdf都是基础 MATLAB 自带的。检查依赖的方法很简单在命令行敲which fminsearch、which readmatrix能返回路径就没问题。如果报未定义函数说明版本太老升级或者替换函数。3.2 数据准备与格式要求数据文件放在data/wind_speed.csv单列每行一个风速值单位 m/s。常见的数据来源是测风塔的十分钟平均风速一年大概五万多条。如果你的数据是多列的比如同时有风速和风向只取风速那一列在main.m里改成data(:, 2)这种索引。数据里如果有-999或者NaN这种缺失标记代码里v data(~isnan(data))只能去掉 NaN去不掉-999。这种情况要在读取后加一句v(v 0) []把负值一并剔除。我一般会在读数据之后先画个boxplot看看有没有离谱的异常值有的话手动处理掉再拟合。注意风速数据不要做平滑或滤波威布尔拟合用的是原始分布平滑会改变分布形状拟合出来的 k 和 c 会偏。3.3 运行与结果解读在 MATLAB 里把当前目录切到项目根目录直接运行main.m。命令行会输出 k、c 和平均风速同时弹出两张图。如果图片没有自动保存在main.m末尾加saveas(gcf, output/figures/weibull_fit.png)即可。结果解读看三点第一直方图和红色 PDF 曲线是否贴合如果直方图明显双峰而曲线是单峰说明风速分布不服从威布尔可能是地形复杂导致的风向分区第二累积分布两条线是否重合尾部偏差大说明极端风速的拟合不好第三拟合平均风速和实测平均风速的偏差一般在 2% 以内算正常超过 5% 要查数据质量。3.4 批量处理多个测风塔实际项目里往往有多个测风塔一个个跑太慢。可以把main.m改成函数外层套一个循环towers {tower1.csv, tower2.csv, tower3.csv}; results zeros(length(towers), 4); % k, c, 拟合均值, 实测均值 for i 1:length(towers) data readmatrix(fullfile(data, towers{i})); v data(~isnan(data)); [k, c] weibull_fit(v); results(i, :) [k, c, c*gamma(11/k), mean(v)]; end disp(array2table(results, VariableNames, ... {k, c, 拟合均值, 实测均值}, ... RowNames, towers));这样一次跑完所有塔结果输出成表格方便对比不同点位的风资源差异。如果某个塔的 k 明显偏离其他塔大概率是数据有问题单独拎出来检查。4. 避坑与常见问题拟合翻车的几种典型情况4.1 拟合曲线和直方图对不上现象运行完main.m红色 PDF 曲线明显偏离直方图要么峰值位置不对要么尾部太厚或太薄。原因最常见的是初值给得不好fminsearch收敛到了局部极值。其次是数据里混了异常值比如传感器故障导致的60 m/s这种离谱风速把分布尾部拉偏了。解决先检查数据用v(v 40) []把明显不合理的值剔掉。然后手动改初值把k0和c0的公式换成更稳的迭代法或者直接用mle函数需要 Statistics Toolbox做对比验证。如果还是对不上考虑风速分布本身不是威布尔可能是混合分布。4.2 报错「未定义函数 readmatrix」现象运行main.m直接报错提示readmatrix未定义。原因MATLAB 版本低于 R2019areadmatrix还没引入。解决把readmatrix(data/wind_speed.csv)换成csvread(data/wind_speed.csv)注意csvread要求文件里全是数值不能有表头。如果有表头用importdata或者readtable再取.Variables。4.3 拟合出的 k 值异常大或异常小现象k 输出 5 以上或者 1 以下和风电场常见值 1.5 到 3.0 差很远。原因k 大于 4 通常说明数据里静风比例太高或者风速变化范围很窄k 小于 1 说明风速分布极其分散可能是数据单位搞错了比如把 km/h 当成了 m/s。解决先确认单位把数据除以 3.6 再试。如果是静风比例高考虑用三参数威布尔加一个位置参数。但风电场评估里两参数够用k 异常大往往是数据质量问题不是模型问题。4.4 累积分布尾部偏差大现象右图里实测累积分布和理论 CDF 在风速大于 15 m/s 的区域分叉明显。原因威布尔分布对尾部拟合天生偏弱尤其是极端风速样本少的时候。另外如果数据里高风速段有截断比如传感器量程上限 25 m/s尾部会失真。解决尾部偏差在风资源评估里影响不大因为风功率密度主要看 3 到 12 m/s 这段。如果非要改善尾部可以用 Weibull 和 Gumbel 混合分布但代码复杂度会上升。我一般会在报告里注明尾部偏差不强行拟合。4.5 批量处理时结果表格行名报错现象用array2table输出结果时RowNames报错提示行名数量不匹配。原因towers是元胞数组array2table的RowNames要求是字符串数组或元胞数组但长度必须和行数一致。如果某个塔的数据读取失败results少了一行行名就对不上。解决在循环里加try-catch读取失败就填NaN保证results行数始终等于towers长度。或者干脆不用RowNames单独输出一列塔名。5. 二次开发与改造从单点拟合到风资源评估工具5.1 增加风功率密度计算风功率密度是风电场选址的核心指标公式是 $P/A \frac{1}{2}\rho c^3 \Gamma(1 3/k)$其中 $\rho$ 是空气密度一般取 1.225 kg/m³。在main.m里加几行rho 1.225; % 空气密度 kg/m^3 power_density 0.5 * rho * c^3 * gamma(1 3/k); fprintf(风功率密度 %.2f W/m^2\n, power_density);这个值算出来一般在 100 到 400 W/m² 之间低于 100 说明风资源较差高于 300 算不错。注意空气密度随海拔变化高原风场要按实际气压温度修正不能直接用 1.225。5.2 封装成函数供其他脚本调用把main.m改造成function [k, c, pd] wind_weibull(csv_path)输入 CSV 路径输出拟合参数和风功率密度。这样其他项目可以直接调用不用每次跑脚本。改造时注意把绘图部分用if nargout 0包起来不输出参数时才画图避免批量调用时弹出一堆窗口。5.3 加入拟合优度检验光看曲线重合不够严谨可以加 R² 或者 Kolmogorov-Smirnov 检验。R² 计算简单[F_emp, x_emp] ecdf(v); F_theory weibull_cdf(x_emp, k, c); R2 1 - sum((F_emp - F_theory).^2) / sum((F_emp - mean(F_emp)).^2); fprintf(拟合 R² %.4f\n, R2);R² 大于 0.98 算拟合良好低于 0.95 要警惕。KS 检验用kstest需要 Statistics Toolbox没有的话 R² 够用。5.4 改造时的验证习惯每次改完拟合逻辑我会拿一组已知 k 和 c 的合成数据做验证用wblrnd生成 10000 个样本跑一遍拟合看输出的 k 和 c 是否接近设定值。偏差在 5% 以内说明代码没问题偏差大就是拟合环节有 bug。这个习惯帮我省了很多次回头查数据的时间。提示合成数据验证时样本量要够大1000 个样本的拟合波动能到 10% 以上10000 个才稳定。5.5 一个容易忽略的细节风速数据的时间分辨率会影响拟合结果。十分钟平均风速和小时平均风速拟合出的 k 值不一样前者通常更大。如果你拿到的数据是小时平均而报告里要求用十分钟数据不能直接混用。我一般会在数据读取后加一个字段记录时间分辨率输出结果时一并打印避免后期扯皮。从那以后我每次做风资源评估都强制先跑一遍合成数据验证拟合代码再拿实测数据出结果。这个习惯看着笨但能挡住大部分玄学问题。希望帮到你。本文还有配套的精品资源点击获取