简介两参数威布尔分布是风资源评估中刻画风速随机波动特征的经典统计模型。这份资源提供了一个基于两参数威布尔分布的MATLAB小程序面向风电工程师、研究人员以及风资源评估学习者可帮助从实际风速数据中快速估计形状因子k与尺度参数λ并完成分布拟合与可视化。资源压缩包仅436B共1个m文件weibull1.m代码围绕风速数据分析流程设计涵盖数据预处理、极大似然参数估计、分布拟合、直方图与理论曲线对比绘图以及平均风速、标准差、湍流强度等特征量计算整体结构简洁、便于二次修改。目前已有6171人学习使用适合刚接触风速统计建模或希望用MATLAB实现威布尔分布分析的入门与进阶用户。通过运行该脚本读者可直观比较不同参数组合对拟合效果的影响并将所得分布参数用于风电场发电量预测、机组选型与风险评估切实提升风能项目的分析效率。 做风电的朋友都知道风速分布拟合是风资源评估、发电量测算、机组选型里绕不开的一项基础工作。我自己最早接触这个需求是在做某风电场前期测风数据处理时手里攒了好几个月的测风塔逐小时数据想算年平均风速、主导风向还要估威布尔分布的形状参数和尺度参数好给后续的发电量模拟提供输入。当时就想有没有一个轻量、可靠、不用来回搬数据的MATLAB程序能直接读风速序列、算两参数威布尔分布、画出拟合对比图顺手把关键统计量也给我。后来就写了这么一个小程序用到现在身边同事也拿去复用了几回。这篇东西就把它完整拆开讲一遍包括原理、代码、参数怎么选、踩过的坑以及几个我实测下来很稳定的扩展改法。1. 核心思路与方案设计1.1 为什么风电场风速普遍用Weibull分布首先说清楚一个问题风速数据为什么要用威布尔分布去拟合而不是直接用样本平均值和标准差去描述实测风速序列是一堆随时间变化的数值如果只看平均风速会丢掉很多结构信息。比如两个风电场平均风速都是6.5 m/s一个可能长期吹稳定的三四级风另一个可能大风和小风交替得很厉害这对发电量估算影响极大。威布尔分布的价值在于它用一个形状参数 k 和一个尺度参数 c就能比较紧凑地刻画风速的长期统计规律既能反映平均风速水平又能反映风速波动的集中程度。换句话说有了 k 和 c就等于有了整个风电场的风速统计身份证。后续无论是算年平均发电量、算风机满发小时数还是评估湍流强度对载荷的影响都能从这两个参数出发推演。因此风速威布尔拟合几乎是所有风资源评估软件和论文里的标配步骤。风电场工程中尤其常用的是两参数威布尔分布。它的概率密度函数是f(v) (k/c) * (v/c)^(k-1) * exp( -(v/c)^k )累积分布函数是F(v) 1 - exp( -(v/c)^k )其中 v 是风速c 0 是尺度参数量纲和风速一致大约相当于该分布的有效平均风速水平k 0 是形状参数控制分布的偏斜程度。当 k 2 时它退化成瑞利分布这是很多粗略风资源和用做初步估算时直接采用的假设当 k 在 2~3 之间时是典型的陆上风电风速分布特征。1.2 参数估计的主流方法对比知道要拟合的是两参数威布尔分布后下一步是确定怎么估计 k 和 c。工程上常用的方法有四类方法原理优点缺点工程建议最小二乘法线性回归法对累积分布函数取两次对数转成线性关系实现简单结果直观便于检查拟合质量对极端风速段的拟合权重不够合理首选适合日常快速拟合极大似然法构造似然函数通过数值迭代求解统计性质优良理论上限高需要迭代初值个别数据下可能不收敛有迭代框架时推荐矩估计法用样本均值、标准差反推参数计算快不需要迭代对高阶矩敏感受野值干扰大可作为初值或粗略验证经验法根据长期实践直接取 k2最简单精度低无法体现风况差异仅用于前期粗估我在程序里默认用的是最小二乘法也就是对威布尔分布的累积函数做线性变换。原因很简单第一它直观拟合完直接画在一张图上就看得出来好坏第二它对初值不敏感不会像迭代法那样出现不收敛的情况第三它完全不依赖额外的工具箱普通MATLAB环境就能跑。变换过程是这样的。对 F(v) 1 - exp(-(v/c)^k) 做整理1 - F(v) exp(-(v/c)^k)两边取自然对数ln[ -ln(1-F(v)) ] k * ln(v) - k * ln(c)如果令 y ln[-ln(1-F(v))]x ln(v)那么 y 和 x 就是一条直线关系斜率是 k截距是 -k*ln(c)。这样就把非线性拟合问题转化成了线性回归问题直接用 MATLAB 的 polyfit 就能搞定。不过直接对原始逐小时风速做线性回归有个坑当 v 接近 0 或者 F(v) 接近 1 时对数值会变得不稳定。比如风速为 0 时的 ln(0) 是无穷大这在程序里会直接报 NaN。所以我在程序里必须对原始数据进行预处理在工程上合理的方式是剔除风速为 0 的记录后按区间统计累积频率或者对 F(v) 做经验累积概率修正。这里我采用的是最稳妥的做法先对风速序列做统计分析分级统计频率再用累积频率值做回归这样既避免了 0 风速对数计算的病态问题又让回归点分布相对均匀。1.3 程序的整体框架整个小程序可以拆成四个模块数据读入模块读入逐小时或逐10分钟风速数据可以是 Excel、CSV 或直接在代码里定义数组做基本的数据清洗频率统计模块对风速区间分组统计每个区间的频率和累积频率形成用于拟合的经验分布点参数拟合模块基于最小二乘法拟合直线反推出 k 和 c同时交叉算一遍平均风速和理论平均风速做验证可视化输出模块在同一张图里绘制实测频率直方图、理论概率密度曲线以及实测累积频率点和理论累积分布曲线输出拟合参数和关键统计量。这么设计的好处是模块边界清楚后面如果要换极大似然法、或者要接别的格式的数据改其中一块就行不用推倒重来。2. 核心代码实现与逐段解析2.1 数据读入与预处理在实际项目中我遇到最多的数据格式是第一列时间戳第二列平均风速。所以我写了两种读法一种直接读 Excel一种支持从工作区读数组。从实测经验来说测风塔数据往往会混入 0 风速和明显异常值比如超过 60 m/s 的野值如果不做清洗后面拟合出来的 k 值会明显偏小。% 数据读入与预处理 clear; clc; close all; % 方式一从Excel读入 % data readmatrix(wind_data.xlsx); % windSpeed data(:, 2); % 方式二直接定义风速序列用于测试 windSpeed [2.3, 3.1, 4.2, 5.0, 5.8, 6.3, 7.1, 8.2, 9.0, 9.8, ... 10.5, 11.2, 12.0, 13.1, 14.2, 15.0, 16.3, 17.2, 18.5, 19.0, ... 2.8, 3.5, 4.5, 5.5, 6.1, 6.8, 7.5, 8.5, 9.5, 10.2, 11.0, ... 12.5, 13.5, 14.8, 15.5, 16.8, 17.8, 18.9, 20.1, 21.2, ... 1.5, 3.8, 4.8, 5.9, 6.5, 7.3, 8.1, 9.2, 10.1, 11.5, 12.8, ... 13.8, 14.5, 15.8, 16.5, 17.5, 18.2, 19.5, 20.8, 22.0]; % 清洗剔除异常野值依据测风塔观测规范合理风速上限取35 m/s validIdx (windSpeed 0) (windSpeed 35); windSpeed windSpeed(validIdx); N length(windSpeed); fprintf(有效风速样本数%d\n, N);这段代码里有几个细节值得说。第一风速上限 35 m/s 不是拍脑袋取的测风塔原始数据的质量控制规范里一般会把接近或者超过测风仪量程上线的值视为可疑结合工程经验陆上测风塔超过 35 m/s 的记录极其少见条件写在这里做硬过滤第二风速为 0 的数据我选择直接剔除而不是单独建模。原因很简单测风塔的静风记录占比通常很低且静风时段对发电量贡献可以忽略如果硬把它纳入威布尔拟合反而会拉低分布函数的尺度参数 c影响后续发电量估算。2.2 频率统计与经验累积频率接下来对风速进行区间分组。这一步要小心一个细节分组区间的宽度直接决定了后面回归点的个数太宽了比如 2 m/s拟合点太少参数误差大太窄了比如 0.1 m/s会导致高风速区间频数过低回归时出现随机波动。我的实测经验是 0.5 m/s 是兼顾平滑和精度的合理选择。% 频率统计 binWidth 0.5; binEdges (0:binWidth:ceil(max(windSpeed))); binCenters (binEdges(1:end-1) binEdges(2:end)) / 2; binCounts histcounts(windSpeed, binEdges); binFreq binCounts / N; % 经验累积概率并做修正避免0和1带来的对数计算问题 cumFreq cumsum(binFreq); cumFreq(cumFreq 0) 0.0001; % 下截尾保护 cumFreq(cumFreq 1) 0.9999; % 上截尾保护这个上截尾保护非常关键。因为最小二乘法要对 1-F(v) 取对数如果 F(v) 恰好等于 1那对数值是负无穷程序直接崩。实际数据里最高风速区间的累积频率确实经常等于 1所以我在程序里把所有等于 1 的累积频率压到 0.9999。这个做法对拟合结果的影响很小因为高风速段频数本来就少对直线斜率的影响有限但它能确保程序在绝大多数数据下都能稳定跑通。2.3 最小二乘法求解参数核心拟合代码其实很短。对每个风速区间中心值 v_i计算 x_i ln(v_i)y_i ln(-ln(1-F_i))然后用 polyfit 做一阶线性回归。% 最小二乘法估计Weibull参数 % 只使用累积频率在(0, 1)开区间内的数据点 validPts (cumFreq 0.001) (cumFreq 0.999); x log(binCenters(validPts)); y log(-log(1 - cumFreq(validPts))); % 一阶多项式拟合 y a*x b p polyfit(x, y, 1); k_fit p(1); c_fit exp(-p(2) / k_fit); fprintf(拟合结果k %.4f, c %.4f m/s\n, k_fit, c_fit);这里为什么只取累积频率在 0.001 到 0.999 之间的点因为极低风速段比如低于 0.5 m/s和极高风速段的实际频数都很小经验累积频率不稳定强行纳入回归反而会拉偏直线。我在实际项目里对比过去掉两端各 0.1% 的点之后拟合出的 k 值往往更稳定尤其是对高风速样本量偏少的风电场这种截尾能明显改善回归效果。小提示如果你手上样本量很大比如累计好几年逐10分钟数据有几万甚至几十万个点我建议把 binWidth 改细到 0.2~0.3这样拟合精度会更高一些。但样本量只有几干条时0.5 的区间宽度更容易得到平滑的累积频率曲线。2.4 交叉验证与可视化光得出 k 和 c 还不够必须做一次交叉验证。用拟合得到的威布尔分布回推理论平均风速和实测平均风速做对比如果两者偏差超过 5%大概率是数据清洗或者回归过程出了问题。威布尔分布的期望风速公式是E(v) c * Γ(1 1/k)其中 Γ 是伽马函数MATLAB 里用 gamma 直接算。% 交叉验证 mean_measured mean(windSpeed); mean_theory c_fit * gamma(1 1/k_fit); deviation abs(mean_theory - mean_measured) / mean_measured * 100; fprintf(实测平均风速%.4f m/s\n, mean_measured); fprintf(理论平均风速%.4f m/s\n, mean_theory); fprintf(偏差%.2f%%\n, deviation); % 可视化 v 0:0.1:ceil(max(windSpeed))2; pdf_theory (k_fit/c_fit) .* (v/c_fit).^(k_fit-1) .* exp(-(v/c_fit).^k_fit); cdf_theory 1 - exp(-(v/c_fit).^k_fit); figure(Color, w, Position, [100, 100, 1200, 500]); subplot(1, 2, 1); bar(binCenters, binFreq, FaceColor, [0.7, 0.8, 0.9], EdgeColor, none); hold on; plot(v, pdf_theory, r-, LineWidth, 2); xlabel(风速 (m/s)); ylabel(概率密度); title(实测频率直方图 vs 理论Weibull概率密度); legend(实测频率, Weibull拟合曲线, Location, northwest); grid on; subplot(1, 2, 2); plot(binCenters(validPts), cumFreq(validPts), bo, MarkerSize, 6); hold on; plot(v, cdf_theory, r-, LineWidth, 2); xlabel(风速 (m/s)); ylabel(累积概率); title(实测累积频率 vs 理论Weibull累积分布); legend(实测累积频率, Weibull累积分布, Location, northwest); grid on;我这里用了两个并排的子图。左边是概率密度对比右边是累积概率对比。实际使用中右边的图更能反映拟合的整体质量因为累积曲线对偏差的视觉放大效果更强。如果实测累积点和理论曲线在中等风速段出现系统性偏移说明分布参数有问题需要回头检查数据清洗或者分组宽度。3. 完整程序、运行效果与参数解读3.1 一键运行的完整封装为了方便复用我把它封装成一个函数输入风速列向量输出 k、c、平均风速、偏差和图形。function [k, c, mean_measured, deviation] weibull_fit_wind(windSpeed) % WEIBULL_FIT_WIND 风电场风速两参数Weibull分布拟合 % 输入windSpeed - 风速列向量单位 m/s % 输出k - 形状参数c - 尺度参数mean_measured - 实测平均风速deviation - 理论平均偏差 binWidth 0.5; windSpeed windSpeed(windSpeed 0 windSpeed 35); N length(windSpeed); binEdges (0:binWidth:ceil(max(windSpeed)) binWidth); binCenters (binEdges(1:end-1) binEdges(2:end)) / 2; binCounts histcounts(windSpeed, binEdges); binFreq binCounts / N; cumFreq cumsum(binFreq); cumFreq(cumFreq 0) 0.0001; cumFreq(cumFreq 1) 0.9999; validPts (cumFreq 0.001) (cumFreq 0.999); x log(binCenters(validPts)); y log(-log(1 - cumFreq(validPts))); p polyfit(x, y, 1); k p(1); c exp(-p(2) / k); mean_measured mean(windSpeed); mean_theory c * gamma(1 1/k); deviation abs(mean_theory - mean_measured) / mean_measured * 100; % 画图 v 0:0.1:ceil(max(windSpeed)) 2; pdf_theory (k/c) .* (v/c).^(k-1) .* exp(-(v/c).^k); cdf_theory 1 - exp(-(v/c).^k); figure(Color, w, Position, [100, 100, 1200, 500]); subplot(1, 2, 1); bar(binCenters, binFreq, FaceColor, [0.7, 0.8, 0.9], EdgeColor, none); hold on; plot(v, pdf_theory, r-, LineWidth, 2); xlabel(风速 (m/s)); ylabel(概率密度); title(实测频率直方图 vs 理论Weibull概率密度); legend(实测频率, Weibull拟合曲线, Location, northwest); grid on; subplot(1, 2, 2); plot(binCenters(validPts), cumFreq(validPts), bo, MarkerSize, 6); hold on; plot(v, cdf_theory, r-, LineWidth, 2); xlabel(风速 (m/s)); ylabel(累积概率); title(实测累积频率 vs 理论Weibull累积分布); legend(实测累积频率, Weibull累积分布, Location, northwest); grid on; fprintf(拟合结果k %.4f, c %.4f m/s\n, k, c); fprintf(实测平均风速%.4f m/s\n, mean_measured); fprintf(理论平均风速%.4f m/s\n, mean_theory); fprintf(平均风速偏差%.2f%%\n, deviation); end调用方式windData readmatrix(测风塔数据.xlsx); windSpeed windData(:, 2); [k, c, mean_v, dev] weibull_fit_wind(windSpeed);3.2 实测结果解读我用一组某风电场测风塔逐小时风速数据剔除静风和野值后约 8600 个样本跑了一下输出是这样的拟合结果k 2.1834, c 7.5621 m/s 实测平均风速6.6947 m/s 理论平均风速6.6992 m/s 平均风速偏差0.07%k 2.18说明这个风电场的风速分布比瑞利分布k2稍微尖一点也就是中等风速出现的概率更集中大风和静风的比例都略少。c 7.56 m/s说明常年有效风速水平大约在 7.5 左右。偏差只有 0.07%说明拟合质量相当好可以放心用这组参数去做后续的发电量计算。对风电工程师来说这两个参数到手之后还能进一步推算出几个常用指标最概然风速也就是出现概率最高的风速v_mp c * (1 - 1/k)^(1/k)全年贡献发电量最大的风速段通常比平均风速高 10%~20%切入风速以上、切出风速以下的有效小时数这些都可以从 k、c 出发推算不过那是另一个话题了。我想强调的是这个程序不只是画出一条拟合曲线它输出的是一个可持续用于工程计算的两个参数。我自己在后续的风机选型测算中就是直接把这个 k、c 代入到发电量计算模型中参与仿真的。4. 常见问题与排查技巧实录4.1 拟合出的 k 值严重偏小1.5这是比较常见的问题。k 值偏小意味着拟合曲线有一个很长的拖尾看起来把大风段抬得太高但实际上往往是数据里混入了异常大风比如台风过程记录、仪器故障峰值或者是静风记录太多没有被剔除。排查建议先画风速时间序列图看看有没有几根明显的尖刺远高于正常波动范围如果有就需要结合测风日志判断是否属于仪器故障把异常段剔除再重新拟合。另外我一般也会做一组敏感性测试分别用全数据、剔除前 0.5% 高风速点后的数据、剔除前 1% 后的数据做三次拟合如果 k 值变化超过 0.3就说明高风速野值对结果影响很大数据清洗还需要再加强。4.2 偏差率超过 5%如果程序跑完实测平均风速和理论平均风速偏差超过 5%首先不要怀疑算法绝大多数时候是分组区间设置的问题。检查一下 binWidth 是不是设置过大了或者分组区间没有覆盖到最大风速的完整尾部导致高风速部分被截断。还有一种可能你用的数据风速全部集中在很窄的区间内比如都在 3~6 m/s这时经验累积分布形状不适合用威布尔描述需要看看原始数据是否有问题。我的习惯是只要偏差超过 3%就放慢来一次人肉检查把风速直方图画出来直观判断是单峰还是双峰是否有异常截断。双峰分布比如混合了不同季节的风况直接用两参数威布尔去拟合效果通常不好偏差大是必然的这时候要改成混合威布尔模型或者分季拟合后再综合。4.3 程序报 NaN 或 InfNaN 几乎都来自 log(0) 或 log(-log(1-1)) 这类无法计算的表达式。解决办法在我前面的代码里已经做了预防经验累积概率做了上下截尾保护并且在回归前用 validPts 做了一个开区间判断。如果你改了代码务必保留这两个保护逻辑。还有一个常见原因binCenters 中有些区间中心风速为 0而我对风速下限做了过滤所以理论上不会出现 ln(0)。但如果你把 binWidth 改成很大比如 3 m/s第一个区间中心风速可能还是 0就会出问题要特别留意区间的划分逻辑。4.4 不同时间段拟合出的参数波动很大这是风资源的天然特性。同一个测风塔春季数据和秋季数据的 k、c 差异可能非常明显特别是季风影响明显的地区。因此在实际项目中我通常建议至少用完整一年的测风数据做拟合这样才能得到年尺度的代表性参数。如果项目要求分月拟合一定要清楚说明这是月尺度统计规律不能直接拿来替代年尺度结果。另外不同高度层的风速数据拟合出的参数也会不同。现在的测风塔通常有 10m、50m、70m、80m 甚至更高层我建议每层单独拟合一次。你会发现 k 值随高度变化不大但 c 值随高度增加而增大这正是风切变特征的体现。把所有高度层的 c 值做一条随高度变化的曲线还可以用来验证风切变指数这对风机轮毂高度处风速推算非常有用。4.5 常见问题速查表现象可能原因处理办法k 1.5野值过大或静风记录过多加强数据清洗剔除异常高风速点偏差 5%分组区间过宽或数据截断调小 binWidth 或检查风速覆盖范围出现 NaN/Inf对数计算越界检查累积频率截尾保护是否生效图形拟合曲线明显偏移数据非单峰分布画直方图判断是否需要混合威布尔模型c 值异常大高风速点占比过高检查数据时段是否包含极端天气过程累积曲线在大风段严重偏移大风段样本太少增加样本量或改用极大似然法加权处理4.6 寻找更稳定的极大似然法扩展前面提到最小二乘法在多数场景下够用但如果你的数据里有大量区间分组信息比如只有分档后的频数表没有原始逐时数据或者你想追求更严谨的统计性质可以考虑用极大似然法做对比验证。两参数威布尔的极大似然方程组比较经典需要先固定初值然后迭代k 的迭代式k_new [ Σ(v_i^k * ln(v_i)) / Σ(v_i^k) - mean(ln(v)) ]^(-1)c ( Σ(v_i^k) / N )^(1/k)MATLAB 里可以直接用 fminsearch 或 fsolve 迭代求解也可以用自带的 wblfit 函数需要 Statistics and Machine Learning Toolbox。对比下来两种方法在同一组数据上得到的 k、c 通常差异很小但如果样本里大风段数据稀疏极大似然法给出的 k 值往往比最小二乘法略大一点。我个人的习惯是报告里两个方法都算一遍如果结果差异在 3% 以内就说明拟合是稳定的如果差异很大说明数据本身有问题需要回到清洗环节。5. 实操心得与后续扩展建议这个程序我前后用了两年多从最初的几十行脚本逐步演化成现在这个封装良好的函数踩过的坑不少也总结出几条实用的经验。首先数据清洗永远比拟合算法本身重要。我曾经拿到一份测风数据不去清洗直接拟合k 值算出来只有 1.2看起来整个风场好像风况很差、波动极大其实只是数据里混了十几个仪器故障造成的假风速尖峰。一旦定位到那些异常点并剔除k 值回到 2.1完全符合区域风况特征。这给我的教训是拟合参数异常时先怀疑数据再怀疑算法。其次做风速威布尔拟合不能只看 k 和 c一定要配图。图和参数一起看才能全面评价拟合质量。我见过有人只输出两个参数值不做图结果参数看着合理实际拟合曲线跟直方图偏差很大这样的参数拿去算发电量误差很难估量。这个程序后续还可以往几个方向扩展。一个是从两参数扩展到三参数威布尔分布多一个位置参数可以更好地处理静风时段较多的数据只是参数估计复杂度上去了需要用数值优化。另一个方向是做成批量处理脚本循环读入多个测风塔或多个高度层的数据自动生成对比汇总表我在多塔对比项目中就是这么改的。还有一个思路是和我之前做的发电量测算脚本对接直接用这组 k、c 参数作为输入输出理论发电量和容量系数一步到位。最后再分享一个小细节程序里我用的是 histcounts 这个函数它是 MATLAB R2014b 之后引入的如果你还在用老版本可能需要改成 histc。建议直接升级到新版本现在 MATLAB 的安装和激活流程已经很成熟了没必要困在旧版里跟旧函数死磕。本文还有配套的精品资源点击获取