简介一款MATLAB小程序专为风电场风速数据分析而设计实现两参数威布尔分布的参数估计与曲线拟合适合风资源评估人员、风电工程师及能源方向研究生使用。压缩包为rar格式仅包含1个M脚本文件总体积约436B短小精悍便于直接运行和二次修改。已有6171人学习下载这一工具。代码围绕风速数据预处理、极大似然估计、分布拟合与可视化展开可自动生成风速直方图与威布尔分布对比曲线并输出形状参数k、尺度参数λ、平均风速、标准差及湍流强度等关键统计量为风电场功率预测、设备选型和风险分析提供统计依据。借助这一脚本读者既能快速掌握威布尔分布的MATLAB实现思路也可将整套流程迁移到自己的风速实测数据中完成从原始数据到概率密度曲线、再到工程指标输出的完整分析闭环。 做风电的人都知道风速数据分析是所有宏观选址、微观选址和发电量估算的起点。而在风速统计的所有模型里两参数威布尔Weibull分布基本是工程界默认的“标配”——精度够用、参数少、物理意义相对清楚。但问题是软件算的是一回事自己能不能用MATLAB把分布参数拟合出来并且拟合得让人放心又是另一回事。这篇就完整记录一下我用MATLAB写“风电场风速两参数威布尔分布计算小程序”的全过程从原理到代码到坑点一次讲清楚。1. 为什么偏偏选两参数威布尔分布做风速统计的时候很多人会纠结风速数据是一堆时序点怎么用一条曲线去描述它的概率规律这就涉及“风速概率密度分布”的概念。简单说风速不是随机乱变它的长期统计规律是可以用特定分布函数去逼近的。工程上最常用的就是两参数威布尔分布概率密度函数是f(v) (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)累积分布函数是F(v) 1 - exp(-(v/c)^k)表达式里只有两个参数——形状参数 k 和尺度参数 c。k 决定曲线的形状反映风速波动的剧烈程度c 和平均风速直接相关决定分布的中心位置。这两个参数一确定平均风速、极大风速概率、风功率密度全都能推导出来所以做风资源评估的人特别依赖这个模型。用别的分布行不行伽马分布、瑞利分布也能凑合但威布尔分布在大量风电场实测数据里的拟合优度普遍更高尤其是对“低风速段概率高、高风速段概率长尾”这种特征刻画得很准。而且参数少意味着估计稳定不会像三参数或四参数模型那样动不动就过拟合。这里必须强调一个应用上的事实切入风速以下的小风速和有功功率为0的时间段在实际采集数据里占比很高。如果你的样本没有经过预处理直接丢进拟合函数威布尔分布的参数会被这一堆“零风速”样本拉偏导致后续发电量估算偏高。这也是我要在程序里单独做数据清洗的原因。2. 参数估计的几种做法和选型两参数威布尔分布的参数估计工程上主流的有三种最小二乘法配合P-Q图线性回归、极大似然法、矩估计法。三种方法我都实际跑过先对比一把。最小二乘法的思路是把威布尔分布转化为线性关系。对累积分布函数两边取两次对数可以得到ln(-ln(1-F)) kln(v) - kln(c)令 Y ln(-ln(1-F))X ln(v)这就变成了一条直线斜率为 k截距为 -k*ln(c)。于是先用实测风速排序算累积频率再做线性回归就能从斜率截距里解出 k 和 c。优点是原理简单、代码好写、计算极快初值找得准缺点是累积频率的估算方式期望公式或中位秩公式会直接影响结果直线拟合对高低风速段的偏差也敏感。极大似然法的思路完全不同求一组参数让当前样本“出现”的概率最大。它的优点是统计性质好参数估计一致且渐近有效所以精度通常更高缺点是需要迭代求解而且对初值敏感——初值给不好迭代可能跑到很离谱的局部极值。矩估计法就是拿样本均值、样本方差联立反解参数但工程上用它的比较少因为矩估计对异常值更敏感在小样本下精度也一般。我的程序里最终采用了“最小二乘法找初值 极大似然法精修”的组合策略。先用线性回归直接给出一个可靠的初值再在初值附近做极大似然迭代既避开了初值难选的问题又保证最终的参数精度。后面完整代码会体现这个过程。3. MATLAB程序的结构和核心代码3.1 主程序框架程序分四层数据输入层、数据清洗层、参数估计层、结果输出层。数据输入考虑两种来源——Excel如测风塔的逐小时风速序列和MATLAB工作区变量代码里用 exist 判断再自动切换。主程序结构如下%% 风电场风速两参数威布尔分布拟合主程序 clc; clear; close all; % 步骤1读取风速数据 % data 是列向量单位m/s [filename, pathname] uigetfile({*.xlsx;*.xls, Excel文件; *.mat, MAT数据文件}, 选择风速数据文件); % ...详细读取分支见下文3.2 数据清洗拟合之前的必要一步实测风速数据基本上不会干净。测风塔故障、极端天气、传感器漂移都会产生各类离群点。我在程序里做了四级处理%% 数据清洗 v_raw data(:); % 强制转列向量 v_raw(v_raw 0 | isnan(v_raw)) []; % 剔除负值和NaN v_raw(v_raw 0) []; % 丢弃静风样本或单独统计 v_raw(v_raw 50) []; % 剔除超大离群点按项目实际调整 upper_bound mean(v_raw) 4*std(v_raw); % 基于均值和标准差的粗差剔除 v v_raw(v_raw upper_bound);这是程序里最容易忽略但最关键的环节。用实测风速做威布尔拟合跟用“数值模拟生成的理论样本”完全是两码事前者永远带着噪声、间断、异常尖峰。处理不好后面拟合出的 k 和 c 会明显偏移。比如我接过一个内陆山地风电场的逐小时风速数据原始记录里有几天的风速值恒为 0.00传感器冻结还有几条文件名里都标了“可疑”的数据直接拿全量去拟合形状参数 k 被拉到 1.2明显低于同类场址的 1.8~2.2。清洗之后就正常了k 稳定在 2.0 左右。所以数据清洗不是为了“显得严谨”而是直接决定参数是否可信。3.3 最小二乘法初值估计核心是用 P-Q 图做线性回归。先把风速排序算累积频率注意频率估计公式的选择。我用的是中位秩公式P_i (i - 0.3) / (n 0.4)这个公式在威布尔概率图检验里表现比较稳。然后用 polyfit 做一阶拟合注意这里要把数量级差异处理好%% 最小二乘初值估计 n length(v); vs sort(v); % 升序排序 P ((1:n) - 0.3) / (n 0.4); % 中位秩累积频率 X log(vs); Y log(-log(1 - P)); % 线性回归: Y a*X b, 对应 a k, b -k*ln(c) idx isfinite(X) isfinite(Y); p_lr polyfit(X(idx), Y(idx), 1); k0 p_lr(1); c0 exp(-p_lr(2) / k0);这一步算出的 k0、c0 能直接当最终结果用也可以作为后面极大似然法的迭代初值。实测下来只要样本量超过 500最小二乘法的初值就已经相当接近极大似然的终值。3.4 极大似然法精修极大似然的似然函数对两参数威布尔分布求偏导后化简得到形状参数 k 的迭代方程。常用的固定点迭代形式是k_new [ sum(v^k * ln(v)) / sum(v^k) - mean(ln(v)) ]^(-1)得到 k 后尺度参数 c 直接用解析式c ( (1/n) * sum(v^k) )^(1/k)MATLAB 里直接用 fzero 或 fsolve 都能解我为了少依赖工具箱用的是 fzero%% 极大似然精修 % 对数似然方程的残差函数 loglik_eq (k) n/k n*log(k) - n*log(c0) ... sum(log(v)) - sum((v/c0).^k .* log(v/c0)); k_ml fzero(loglik_eq, [0.3, 6]); % 在 [0.3,6] 内找根 c_ml (mean(v.^k_ml))^(1/k_ml);迭代区间 [0.3, 6] 覆盖了风电行业几乎所有的实际工况——海上风电 k 通常 1.5~2.5内陆平原 1.8~2.2复杂山地可能到 2.5 左右个别强阵风地区也会低到 1.2。把区间放宽一点永远没坏处fzero 在区间端点符号相反的条件下一定会返回一个根。3.5 完整计算函数封装实际项目中最好把拟合功能封装成函数方便复用于多个测风塔、多个高度层和批量处理function [k, c, stats] fit_weibull_2p(v) % 输入: v - 风速样本(m/s) % 输出: k - 形状参数, c - 尺度参数(m/s), stats - 拟合评价结构体 % 数据清洗 v v(:); v v(isfinite(v)); v(v 0) []; v(v 50) []; % 最小二乘初值 n length(v); vs sort(v); P ((1:n) - 0.3) / (n 0.4); X log(vs); Y log(-log(1 - P)); p_lr polyfit(X(isfinite(Y)), Y(isfinite(Y)), 1); k0 max(p_lr(1), 0.1); c0 exp(-p_lr(2) / k0); % 极大似然精修 if n 1 loglik_eq (k) n/k n*log(k) - n*log(c0) ... sum(log(v)) - sum((v/c0).^k .* log(v/c0)); k fzero(loglik_eq, [0.3, 6]); c (mean(v.^k))^(1/k); else k k0; c c0; end % 拟合优度评价 theoretical_cdf 1 - exp(-(vs/c).^k); [~, stats.r2] corr(P, theoretical_cdf); stats.k0 k0; stats.c0 c0; stats.mean_v mean(v); stats.std_v std(v); end这段代码在 i5 处理器上处理 8760 个逐小时风速数据点跑完只需要几十毫秒。加了批处理之后对三座测风塔、四个高度层共 12 组数据做循环拟合总耗时不到一秒完全满足“快速评估”的场景需求。4. 拟合效果的可视化与评估参数算完了不能直接收工。不做可视化、不做拟合优度检验就没法跟业主或审图方交代。我的程序里画两组图概率密度对比图和 P-Q 图。概率密度对比图的画法是把实测风速直方图的频率密度注意是密度不是频数跟威布尔理论密度曲线叠在一起。直方图必须归一化成密度否则量纲不对曲线会错位%% 概率密度对比图 figure(Color, w); histogram(v, 0:1:max(v), Normalization, pdf, FaceColor, [0.7 0.8 0.9]); hold on; v_plot linspace(0, max(v), 200); pdf_theory (k/c) .* (v_plot/c).^(k-1) .* exp(-(v_plot/c).^k); plot(v_plot, pdf_theory, r-, LineWidth, 2); xlabel(风速 (m/s)); ylabel(概率密度); legend(实测频率密度, 威布尔拟合, Location, northeast); title(sprintf(风速概率密度对比 (k%.2f, c%.2f m/s), k, c)); grid on;P-Q 图概率图则用来做“直线性检验”如果数据点基本落在一条直线上说明威布尔分布拟合良好。偏差集中在两端是正常现象但中间段如果呈 S 形弯曲就要怀疑数据是否混合了不同风况。%% P-Q 图 figure(Color, w); plot(X, Y, b., MarkerSize, 6); hold on; x_fit linspace(min(X), max(X), 50); y_fit p_lr(1) * x_fit p_lr(2); plot(x_fit, y_fit, r-, LineWidth, 1.5); xlabel(ln(v)); ylabel(ln(-ln(1-F))); title(威布尔概率图); legend(样本点, 线性拟合, Location, northwest); grid on;这里补充一个常被忽视的操作直方图的 bin 宽度会直接影响视觉判断。如果 bin 太宽峰值位置会被抹平太窄又会出现锯齿状噪声。我用 1 m/s 的 bin 宽度通常效果不错如果你的数据量很大或很小可以调整为 0.5 m/s 或 2 m/s。拟合优度的量化指标程序里输出了两组R²实测累积频率与理论累积频率的相关系数平方和 RMSE。R² 在 0.98 以上算拟合良好0.99 以上算优秀。但要注意R² 对“中间段”数据天然友好对两端偏差不敏感所以别只看 R² 就判定数据合格还是要眼神扫一遍 P-Q 图。5. 常见问题和排查技巧这个程序虽然不长但实际运行中常出问题。把我在不同类型风电场数据上踩过的坑整理一下。第一类数据读取失败或者格式不对。用 uigetfile 选择 Excel 时如果风速数据是用文本格式存的数字比如从 SCADA 直接导出的 CSV 里带单位readmatrix 会读出一堆 NaN。解决办法是读文件之前先 preview 一下或者强制用 detectImportOptions 设置 VariableNamingRule。程序内部我用了一个分支如果是 Excel尝试 xlsread 和 readmatrix 两套方案失败则给出明显报错提示。try data readmatrix(fullfile(pathname, filename)); catch data xlsread(fullfile(pathname, filename)); end第二类负风速和零风速的处理逻辑。负风速肯定是错误数据直接删。但零风速要斟酌如果时间序列里包含大量检修停机、传感器故障导致的零值直接删合理如果是真实静风气象无风天那么在评估“风机实际出力”时应该保留但在评估“该场址的潜在风资源”时可以单独统计静风比例。我在代码里默认丢弃零值但在结果输出里单独统计了 zero_ratio这样后续需要时可以快速切换策略。第三类数据量太少或风速范围偏窄。只有两三个月的短期测风数据或者数据全集中在 3~8 m/s 的狭窄区间最小二乘法的回归容易不稳定极大似然的收敛也会慢。低于 200 个有效风速样本时建议直接把极大似然迭代去掉只用最小二乘结果或者手动检查初值的合理性。第四类k 值极端异常小于 0.8 或大于 4。出现这种情况要先怀疑数据来源不是程序算错了。我在实际项目中遇到过一次 k 0.9 的情况后来查明数据是多个测风塔不同高度混合在一起的拼接文件风速序列有明显的多峰特征这种混合样本当然不符合单一威布尔分布假设。遇到这种情况程序里的拟合评价会把 R² 拉低到 0.9 以下P-Q 图也会显示明显的弯折——一眼就能看出来。解决办法是拆分成单塔单高度重新拟合或者改用混合威布尔分布模型。第五类fzero 无法找到根。少数情况下数据清洗后样本集为空、或者所有风速集中在一个常数值比如全是 6.0对数累积频率会变成无穷fzero 找不到合适的根。程序里加了一个保护判断如果样本标准差小于 0.01就直接用均值作为 c、k 设为 10并给出“数据退化为近似恒速风”的提示。这种情况虽然少见但在检修记录数据里确实可能发生。6. 结果输出的实际应用参数 k 和 c 拟合出来后风资源评估里最重要的几个量可以直接算平均风速估算、最大可能风速也就是若干年一遇的极值风速、风功率密度。平均风速的威布尔估计是V_mean c * Gamma(1 1/k)MATLAB 里 gamma 函数直接用V_mean_fit c * gamma(1 1/k);风功率密度是WPD 0.5 * rho * c^3 * Gamma(1 3/k)其中 rho 取 1.225 kg/m³标准空气密度。这两个量是发电量估算输入的关键。我在程序里把这两个量直接打印在输出区还加上一句提示如果实测平均风速和威布尔估计平均风速相差超过 0.5 m/s大概率是数据清洗环节出了问题。这个差值本身就是一个很好的内部校验指标。有一说一做风电项目的同事肯定也总结过类似经验。有一次我用实测全年 8760 小时的数据去拟合威布尔估计的平均风速 6.9 m/s而样本直接平均风速 7.0 m/s两个数几乎重合。这说明拟合效果非常好。反过来如果两个数差得远除非样本量特别小否则一定是数据有问题而不是模型有问题。这个“内外一致”的检查方法比单纯看 R² 更直观。7. 扩展方向和小结把一个基础的两参数威布尔拟合小程序做好之后很多延伸工作都变得很容易。我后续在这个程序基础上扩展了三个方向一是批量处理多测风塔数据输出各高度层的 k/c 参数表格二是做逐月威布尔参数滚动分析观察不同季节的分布形态变化三是把拟合结果直接输出成 WindPRO 和 WAsP 需要的格式减少人工转数据的环节。从我个人的实际体验来看这个程序最有价值的点在于把“数据清洗—初值估计—极大似然—图形诊断”的完整流程固化下来。特别是数据清洗和图形诊断两步是区分“会用 fitdist”和“能做风资源工程分析”的分水岭。你光调用 fitdist 函数得到两个参数不检查数据、不看 P-Q 图、不做平均风速交叉验证很难保证结果是可信的。最后一个小建议程序里的中位秩公式、迭代区间、离群值剔除阈值都不是拍脑袋定的而是根据风电行业测风数据的典型特征调的。如果你要把这套程序用到其他领域比如风速以外的气象要素、或者工业设备震动风速传感器数据这些超参需要重新校核不要照搬。代码本身不复杂真正的复杂度永远在数据里。本文还有配套的精品资源点击获取
基于MATLAB的风速威布尔分布拟合：从原理到代码实现
简介一款MATLAB小程序专为风电场风速数据分析而设计实现两参数威布尔分布的参数估计与曲线拟合适合风资源评估人员、风电工程师及能源方向研究生使用。压缩包为rar格式仅包含1个M脚本文件总体积约436B短小精悍便于直接运行和二次修改。已有6171人学习下载这一工具。代码围绕风速数据预处理、极大似然估计、分布拟合与可视化展开可自动生成风速直方图与威布尔分布对比曲线并输出形状参数k、尺度参数λ、平均风速、标准差及湍流强度等关键统计量为风电场功率预测、设备选型和风险分析提供统计依据。借助这一脚本读者既能快速掌握威布尔分布的MATLAB实现思路也可将整套流程迁移到自己的风速实测数据中完成从原始数据到概率密度曲线、再到工程指标输出的完整分析闭环。 做风电的人都知道风速数据分析是所有宏观选址、微观选址和发电量估算的起点。而在风速统计的所有模型里两参数威布尔Weibull分布基本是工程界默认的“标配”——精度够用、参数少、物理意义相对清楚。但问题是软件算的是一回事自己能不能用MATLAB把分布参数拟合出来并且拟合得让人放心又是另一回事。这篇就完整记录一下我用MATLAB写“风电场风速两参数威布尔分布计算小程序”的全过程从原理到代码到坑点一次讲清楚。1. 为什么偏偏选两参数威布尔分布做风速统计的时候很多人会纠结风速数据是一堆时序点怎么用一条曲线去描述它的概率规律这就涉及“风速概率密度分布”的概念。简单说风速不是随机乱变它的长期统计规律是可以用特定分布函数去逼近的。工程上最常用的就是两参数威布尔分布概率密度函数是f(v) (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)累积分布函数是F(v) 1 - exp(-(v/c)^k)表达式里只有两个参数——形状参数 k 和尺度参数 c。k 决定曲线的形状反映风速波动的剧烈程度c 和平均风速直接相关决定分布的中心位置。这两个参数一确定平均风速、极大风速概率、风功率密度全都能推导出来所以做风资源评估的人特别依赖这个模型。用别的分布行不行伽马分布、瑞利分布也能凑合但威布尔分布在大量风电场实测数据里的拟合优度普遍更高尤其是对“低风速段概率高、高风速段概率长尾”这种特征刻画得很准。而且参数少意味着估计稳定不会像三参数或四参数模型那样动不动就过拟合。这里必须强调一个应用上的事实切入风速以下的小风速和有功功率为0的时间段在实际采集数据里占比很高。如果你的样本没有经过预处理直接丢进拟合函数威布尔分布的参数会被这一堆“零风速”样本拉偏导致后续发电量估算偏高。这也是我要在程序里单独做数据清洗的原因。2. 参数估计的几种做法和选型两参数威布尔分布的参数估计工程上主流的有三种最小二乘法配合P-Q图线性回归、极大似然法、矩估计法。三种方法我都实际跑过先对比一把。最小二乘法的思路是把威布尔分布转化为线性关系。对累积分布函数两边取两次对数可以得到ln(-ln(1-F)) kln(v) - kln(c)令 Y ln(-ln(1-F))X ln(v)这就变成了一条直线斜率为 k截距为 -k*ln(c)。于是先用实测风速排序算累积频率再做线性回归就能从斜率截距里解出 k 和 c。优点是原理简单、代码好写、计算极快初值找得准缺点是累积频率的估算方式期望公式或中位秩公式会直接影响结果直线拟合对高低风速段的偏差也敏感。极大似然法的思路完全不同求一组参数让当前样本“出现”的概率最大。它的优点是统计性质好参数估计一致且渐近有效所以精度通常更高缺点是需要迭代求解而且对初值敏感——初值给不好迭代可能跑到很离谱的局部极值。矩估计法就是拿样本均值、样本方差联立反解参数但工程上用它的比较少因为矩估计对异常值更敏感在小样本下精度也一般。我的程序里最终采用了“最小二乘法找初值 极大似然法精修”的组合策略。先用线性回归直接给出一个可靠的初值再在初值附近做极大似然迭代既避开了初值难选的问题又保证最终的参数精度。后面完整代码会体现这个过程。3. MATLAB程序的结构和核心代码3.1 主程序框架程序分四层数据输入层、数据清洗层、参数估计层、结果输出层。数据输入考虑两种来源——Excel如测风塔的逐小时风速序列和MATLAB工作区变量代码里用 exist 判断再自动切换。主程序结构如下%% 风电场风速两参数威布尔分布拟合主程序 clc; clear; close all; % 步骤1读取风速数据 % data 是列向量单位m/s [filename, pathname] uigetfile({*.xlsx;*.xls, Excel文件; *.mat, MAT数据文件}, 选择风速数据文件); % ...详细读取分支见下文3.2 数据清洗拟合之前的必要一步实测风速数据基本上不会干净。测风塔故障、极端天气、传感器漂移都会产生各类离群点。我在程序里做了四级处理%% 数据清洗 v_raw data(:); % 强制转列向量 v_raw(v_raw 0 | isnan(v_raw)) []; % 剔除负值和NaN v_raw(v_raw 0) []; % 丢弃静风样本或单独统计 v_raw(v_raw 50) []; % 剔除超大离群点按项目实际调整 upper_bound mean(v_raw) 4*std(v_raw); % 基于均值和标准差的粗差剔除 v v_raw(v_raw upper_bound);这是程序里最容易忽略但最关键的环节。用实测风速做威布尔拟合跟用“数值模拟生成的理论样本”完全是两码事前者永远带着噪声、间断、异常尖峰。处理不好后面拟合出的 k 和 c 会明显偏移。比如我接过一个内陆山地风电场的逐小时风速数据原始记录里有几天的风速值恒为 0.00传感器冻结还有几条文件名里都标了“可疑”的数据直接拿全量去拟合形状参数 k 被拉到 1.2明显低于同类场址的 1.8~2.2。清洗之后就正常了k 稳定在 2.0 左右。所以数据清洗不是为了“显得严谨”而是直接决定参数是否可信。3.3 最小二乘法初值估计核心是用 P-Q 图做线性回归。先把风速排序算累积频率注意频率估计公式的选择。我用的是中位秩公式P_i (i - 0.3) / (n 0.4)这个公式在威布尔概率图检验里表现比较稳。然后用 polyfit 做一阶拟合注意这里要把数量级差异处理好%% 最小二乘初值估计 n length(v); vs sort(v); % 升序排序 P ((1:n) - 0.3) / (n 0.4); % 中位秩累积频率 X log(vs); Y log(-log(1 - P)); % 线性回归: Y a*X b, 对应 a k, b -k*ln(c) idx isfinite(X) isfinite(Y); p_lr polyfit(X(idx), Y(idx), 1); k0 p_lr(1); c0 exp(-p_lr(2) / k0);这一步算出的 k0、c0 能直接当最终结果用也可以作为后面极大似然法的迭代初值。实测下来只要样本量超过 500最小二乘法的初值就已经相当接近极大似然的终值。3.4 极大似然法精修极大似然的似然函数对两参数威布尔分布求偏导后化简得到形状参数 k 的迭代方程。常用的固定点迭代形式是k_new [ sum(v^k * ln(v)) / sum(v^k) - mean(ln(v)) ]^(-1)得到 k 后尺度参数 c 直接用解析式c ( (1/n) * sum(v^k) )^(1/k)MATLAB 里直接用 fzero 或 fsolve 都能解我为了少依赖工具箱用的是 fzero%% 极大似然精修 % 对数似然方程的残差函数 loglik_eq (k) n/k n*log(k) - n*log(c0) ... sum(log(v)) - sum((v/c0).^k .* log(v/c0)); k_ml fzero(loglik_eq, [0.3, 6]); % 在 [0.3,6] 内找根 c_ml (mean(v.^k_ml))^(1/k_ml);迭代区间 [0.3, 6] 覆盖了风电行业几乎所有的实际工况——海上风电 k 通常 1.5~2.5内陆平原 1.8~2.2复杂山地可能到 2.5 左右个别强阵风地区也会低到 1.2。把区间放宽一点永远没坏处fzero 在区间端点符号相反的条件下一定会返回一个根。3.5 完整计算函数封装实际项目中最好把拟合功能封装成函数方便复用于多个测风塔、多个高度层和批量处理function [k, c, stats] fit_weibull_2p(v) % 输入: v - 风速样本(m/s) % 输出: k - 形状参数, c - 尺度参数(m/s), stats - 拟合评价结构体 % 数据清洗 v v(:); v v(isfinite(v)); v(v 0) []; v(v 50) []; % 最小二乘初值 n length(v); vs sort(v); P ((1:n) - 0.3) / (n 0.4); X log(vs); Y log(-log(1 - P)); p_lr polyfit(X(isfinite(Y)), Y(isfinite(Y)), 1); k0 max(p_lr(1), 0.1); c0 exp(-p_lr(2) / k0); % 极大似然精修 if n 1 loglik_eq (k) n/k n*log(k) - n*log(c0) ... sum(log(v)) - sum((v/c0).^k .* log(v/c0)); k fzero(loglik_eq, [0.3, 6]); c (mean(v.^k))^(1/k); else k k0; c c0; end % 拟合优度评价 theoretical_cdf 1 - exp(-(vs/c).^k); [~, stats.r2] corr(P, theoretical_cdf); stats.k0 k0; stats.c0 c0; stats.mean_v mean(v); stats.std_v std(v); end这段代码在 i5 处理器上处理 8760 个逐小时风速数据点跑完只需要几十毫秒。加了批处理之后对三座测风塔、四个高度层共 12 组数据做循环拟合总耗时不到一秒完全满足“快速评估”的场景需求。4. 拟合效果的可视化与评估参数算完了不能直接收工。不做可视化、不做拟合优度检验就没法跟业主或审图方交代。我的程序里画两组图概率密度对比图和 P-Q 图。概率密度对比图的画法是把实测风速直方图的频率密度注意是密度不是频数跟威布尔理论密度曲线叠在一起。直方图必须归一化成密度否则量纲不对曲线会错位%% 概率密度对比图 figure(Color, w); histogram(v, 0:1:max(v), Normalization, pdf, FaceColor, [0.7 0.8 0.9]); hold on; v_plot linspace(0, max(v), 200); pdf_theory (k/c) .* (v_plot/c).^(k-1) .* exp(-(v_plot/c).^k); plot(v_plot, pdf_theory, r-, LineWidth, 2); xlabel(风速 (m/s)); ylabel(概率密度); legend(实测频率密度, 威布尔拟合, Location, northeast); title(sprintf(风速概率密度对比 (k%.2f, c%.2f m/s), k, c)); grid on;P-Q 图概率图则用来做“直线性检验”如果数据点基本落在一条直线上说明威布尔分布拟合良好。偏差集中在两端是正常现象但中间段如果呈 S 形弯曲就要怀疑数据是否混合了不同风况。%% P-Q 图 figure(Color, w); plot(X, Y, b., MarkerSize, 6); hold on; x_fit linspace(min(X), max(X), 50); y_fit p_lr(1) * x_fit p_lr(2); plot(x_fit, y_fit, r-, LineWidth, 1.5); xlabel(ln(v)); ylabel(ln(-ln(1-F))); title(威布尔概率图); legend(样本点, 线性拟合, Location, northwest); grid on;这里补充一个常被忽视的操作直方图的 bin 宽度会直接影响视觉判断。如果 bin 太宽峰值位置会被抹平太窄又会出现锯齿状噪声。我用 1 m/s 的 bin 宽度通常效果不错如果你的数据量很大或很小可以调整为 0.5 m/s 或 2 m/s。拟合优度的量化指标程序里输出了两组R²实测累积频率与理论累积频率的相关系数平方和 RMSE。R² 在 0.98 以上算拟合良好0.99 以上算优秀。但要注意R² 对“中间段”数据天然友好对两端偏差不敏感所以别只看 R² 就判定数据合格还是要眼神扫一遍 P-Q 图。5. 常见问题和排查技巧这个程序虽然不长但实际运行中常出问题。把我在不同类型风电场数据上踩过的坑整理一下。第一类数据读取失败或者格式不对。用 uigetfile 选择 Excel 时如果风速数据是用文本格式存的数字比如从 SCADA 直接导出的 CSV 里带单位readmatrix 会读出一堆 NaN。解决办法是读文件之前先 preview 一下或者强制用 detectImportOptions 设置 VariableNamingRule。程序内部我用了一个分支如果是 Excel尝试 xlsread 和 readmatrix 两套方案失败则给出明显报错提示。try data readmatrix(fullfile(pathname, filename)); catch data xlsread(fullfile(pathname, filename)); end第二类负风速和零风速的处理逻辑。负风速肯定是错误数据直接删。但零风速要斟酌如果时间序列里包含大量检修停机、传感器故障导致的零值直接删合理如果是真实静风气象无风天那么在评估“风机实际出力”时应该保留但在评估“该场址的潜在风资源”时可以单独统计静风比例。我在代码里默认丢弃零值但在结果输出里单独统计了 zero_ratio这样后续需要时可以快速切换策略。第三类数据量太少或风速范围偏窄。只有两三个月的短期测风数据或者数据全集中在 3~8 m/s 的狭窄区间最小二乘法的回归容易不稳定极大似然的收敛也会慢。低于 200 个有效风速样本时建议直接把极大似然迭代去掉只用最小二乘结果或者手动检查初值的合理性。第四类k 值极端异常小于 0.8 或大于 4。出现这种情况要先怀疑数据来源不是程序算错了。我在实际项目中遇到过一次 k 0.9 的情况后来查明数据是多个测风塔不同高度混合在一起的拼接文件风速序列有明显的多峰特征这种混合样本当然不符合单一威布尔分布假设。遇到这种情况程序里的拟合评价会把 R² 拉低到 0.9 以下P-Q 图也会显示明显的弯折——一眼就能看出来。解决办法是拆分成单塔单高度重新拟合或者改用混合威布尔分布模型。第五类fzero 无法找到根。少数情况下数据清洗后样本集为空、或者所有风速集中在一个常数值比如全是 6.0对数累积频率会变成无穷fzero 找不到合适的根。程序里加了一个保护判断如果样本标准差小于 0.01就直接用均值作为 c、k 设为 10并给出“数据退化为近似恒速风”的提示。这种情况虽然少见但在检修记录数据里确实可能发生。6. 结果输出的实际应用参数 k 和 c 拟合出来后风资源评估里最重要的几个量可以直接算平均风速估算、最大可能风速也就是若干年一遇的极值风速、风功率密度。平均风速的威布尔估计是V_mean c * Gamma(1 1/k)MATLAB 里 gamma 函数直接用V_mean_fit c * gamma(1 1/k);风功率密度是WPD 0.5 * rho * c^3 * Gamma(1 3/k)其中 rho 取 1.225 kg/m³标准空气密度。这两个量是发电量估算输入的关键。我在程序里把这两个量直接打印在输出区还加上一句提示如果实测平均风速和威布尔估计平均风速相差超过 0.5 m/s大概率是数据清洗环节出了问题。这个差值本身就是一个很好的内部校验指标。有一说一做风电项目的同事肯定也总结过类似经验。有一次我用实测全年 8760 小时的数据去拟合威布尔估计的平均风速 6.9 m/s而样本直接平均风速 7.0 m/s两个数几乎重合。这说明拟合效果非常好。反过来如果两个数差得远除非样本量特别小否则一定是数据有问题而不是模型有问题。这个“内外一致”的检查方法比单纯看 R² 更直观。7. 扩展方向和小结把一个基础的两参数威布尔拟合小程序做好之后很多延伸工作都变得很容易。我后续在这个程序基础上扩展了三个方向一是批量处理多测风塔数据输出各高度层的 k/c 参数表格二是做逐月威布尔参数滚动分析观察不同季节的分布形态变化三是把拟合结果直接输出成 WindPRO 和 WAsP 需要的格式减少人工转数据的环节。从我个人的实际体验来看这个程序最有价值的点在于把“数据清洗—初值估计—极大似然—图形诊断”的完整流程固化下来。特别是数据清洗和图形诊断两步是区分“会用 fitdist”和“能做风资源工程分析”的分水岭。你光调用 fitdist 函数得到两个参数不检查数据、不看 P-Q 图、不做平均风速交叉验证很难保证结果是可信的。最后一个小建议程序里的中位秩公式、迭代区间、离群值剔除阈值都不是拍脑袋定的而是根据风电行业测风数据的典型特征调的。如果你要把这套程序用到其他领域比如风速以外的气象要素、或者工业设备震动风速传感器数据这些超参需要重新校核不要照搬。代码本身不复杂真正的复杂度永远在数据里。本文还有配套的精品资源点击获取
