做风光资源评估的人，十有八九都干过这么一件事：拿一年测风数据，点开Matlab，histogram一下，再拿wblfit拟合一条Weibull曲线叠上去。光伏那边也类似，把归一化后的辐照度数据用betafit拟合一条Beta分布曲线。这两个分布，一个管风电，一个管光电，是新能源概率建模里用得最多的一对搭档。但多数教程只告诉你要用这个函数、那个工具箱，不讲清楚为什么偏偏是这两种分布，参数背后对应什么物理过程，更不会讲怎么把Weibull和Beta模型组合起来，系统性地研究风光互补特性。这篇文章就把这条完整链路捋一遍：先从数学和物理背景讲清楚选型逻辑，再给出可直接复制的Matlab代码，包括参数估计、拟合优度检验、蒙特卡洛组合抽样，最后放上我实际项目里踩过的坑和排查方法。不管你是做风电场测风数据分析、光伏电站辐照度建模，还是微电网容量规划，这套流程都能拿来直接用。

1. 为什么偏偏是Weibull和Beta

1.1 Weibull分布不是巧合，而是风速统计的必然

很多人第一次接触Weibull分布是在风资源评估报告里，往往只看结论：平均风速多少、Weibull形状参数k和尺度参数c是多少。但真正要理解的是，为什么风速数据会呈现出典型的右偏分布。

风速的物理来源是气压梯度力，但地表的摩擦拖拽、地形抬升、湍流交换会把这股规则的气流搅得支离破碎。实测风速从来不是对称的钟形曲线，而是拖着一条长长的右尾，也就是高风速时段虽然少，但确实存在。同时风速不可能为负，这就把分布严格限制在半轴上。Weibull密度函数恰好具备这两个特征：定义域在[0,∞)，形状参数k能灵活调节偏度。k=1时退化为指数分布，k=2时就是瑞利分布，k在3左右已经非常接近正态的形态。这个可塑性是正态分布给不了的，因为正态分布允许负值、形态固定对称，和风速的真实统计行为差得太远。

数学上，风速v的Weibull概率密度写成：

f(v) = (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)

其中c是无量纲化的尺度参数，跟平均风速正相关；k决定分布形状，工程上常说k越大风速越稳定。更有用的是期望风速可以解析算出来：

E[v] = c * Γ(1 + 1/k)

这里Γ是伽马函数。这意味着只要从测风数据里拟合出k和c，不用翻原始数据就能估算年平均风速，进而估算理论发电量。我实际做风电场预可研时，经常用这个公式交叉验证测风塔数据的合理性：如果拟合出的期望风速和实测平均风速差超过5%，基本可以判断拟合参数设置有问题，或者原始数据里有严重异常段。

1.2 Beta分布天生为有界辐照度服务

太阳辐照度和风速完全不同。辐照度的物理上限非常明确，大气上界的太阳常数是确定的，再经过大气衰减、云层遮挡，落到水平面上的值始终被约束在一个区间内。归一化之后，辐照度就是[0,1]之间的有界变量。Beta分布恰好就是定义在(0,1)区间上的分布族，概率密度为：

f(x) = x^(α-1) * (1-x)^(β-1) / B(α, β)

其中B(α,β)是Beta函数。α和β两个参数控制分布形态，自由度相当高：α<1、β<1时密度呈U形，对应那种晴雨边界特别明显的天气；α和β都大于1时，密度中间高两边低，对应天气相对稳定的情况；α=β=1时退化为均匀分布。做光伏发电量概率评估时，一个站点在不同季节的辐照度数据往往呈现出完全不同的形态特征，Beta分布用两个参数就能把这些形态差异都表达出来。

但这里有个实操中特别容易踩的坑：Beta分布的支撑集是开区间(0,1)，而实测辐照度数据经常包含极接近0的阴天数据和接近1的晴空极值。强行把包含0和1的数据喂给betafit，轻则参数估计不稳，重则直接算出NaN。合理做法是只对白天且辐照度大于某个阈值的样本建模，或者对数据进行微小边界偏移，把0和1映射到1e-6和1-1e-6这种级别。具体处理方式后面代码部分会详细说。

1.3 两款分布选型的对照逻辑

用一个表格把这套选型逻辑收拢一下，方便后续理解代码设计：

选型不当的后果，我在项目里见过不少。有人硬用正态分布拟合风速，拟合出来的分位数在低风速段直接变成负值，后续做储能容量配置时算出荒谬结果；也有人把夜间辐照度一起塞进Beta拟合，α、β参数被大量零值带偏，绘制出的分布曲线和白天真实光伏出力完全对不上。所以选定分布类型这步，不是"工具箱里有啥用啥"，而是要先想清楚变量的物理边界和统计形态。

2. Matlab实现：从数据清洗到参数估计

2.1 数据预处理决定拟合质量的上限

很多人在这一步翻车。测风塔给的数据通常是10分钟或1小时平均风速，光伏电站给的是小时级水平面总辐照度GHI。直接拿原始数据去拟合之前，至少要做三件事。

第一，剔除物理上不可能出现的值。风速小于0或超过120m/s基本是传感器故障；辐照度出现负值、或超过该站点历史晴空上限太多，一般是辐射表进入阴影或信号漂移。第二，处理缺测。简单删除记录会导致时间序列不完整，而线性插值在长缺测段又会产生虚假平台。我常用的方案是：缺测少于连续3个点用线性插值，超过3个点就丢弃该段，避免把插值出来的假数据当成真实统计样本。第三，辐照度归一化的上限要选对。用历史绝对最大值会把分布右边界拉得很远，让大多数样本集中在小值区域，Beta拟合失真。更好的选择是取全年小时辐照度数据的95%分位数作为归一化基准，这样既保留极端晴空信息，又不会让个别传感器毛刺主导整个分布。

归一化这一步的处理逻辑，我单独说明一下。设原始辐照度为G，选定的上限为G_ref，则x = G / G_ref。理论上x可以等于1，但Beta拟合需要严格小于1，所以工程上再做一层保护：

x_safe = min(max(x, 1e-6), 1 - 1e-6)

这个偏移不影响统计结论，但能保证后续betafit稳定收敛。

2.2 Weibull参数估计的两条路线

Matlab统计工具箱提供了wblfit函数，可以直接对风速样本做极大似然估计，返回[k_hat, c_hat]以及置信区间。一行代码搞定：

[k_hat, c_hat] = wblfit(v);

但我还是建议手动把MLE推导一遍，因为理解机理能让你更快调试异常结果。风速样本v_1,...,v_n的对数似然函数对k求导后，可以得到k满足的方程：

1/k = sum(v_i^k * log(v_i)) / sum(v_i^k) - sum(log(v_i)) / n

这个方程没有解析解，但左侧是单调递减函数，右侧是单调递增函数，用fzero几下就能收敛。代码如下：

% 手动求解Weibull MLE v(v <= 0) = []; % 删除非正风速 n = length(v); k = fzero(@(kk) 1/kk - sum(v.^kk .* log(v)) / sum(v.^kk) + sum(log(v)) / n, 2); c = (sum(v.^k) / n)^(1/k);

fzero初始值取2很安全，因为绝大多数实际风速数据k都落在1到4之间。如果你发现初始值对结果敏感，说明数据本身有问题，比如大量静风零值混在里面，那种情况应该用零截断Weibull或者混合模型，而不是硬套标准Weibull。

拿到k和c后，立刻做两个检查：计算期望风速c*gamma(1+1/k)跟实测算术平均风速比一比；再画出经验CDF和理论CDF，看中段有没有明显系统偏离。这些检查是wblfit不会替你做的。

2.3 Beta参数估计的稳定做法

beta分布跟Matlab工具箱里的betafit也能直接对接：

% 输入x必须是(0,1)范围内的样本 alpha_hat = betafit(x);

betafit返回的是两个输出，实际上返回的是[alpha, beta]，别漏了。它内部使用的也是极大似然估计，通过数值优化迭代求解，不需要我们手动推导，但有一个前提：样本必须严格落在(0,1)开区间。这也是我在2.1节强调边界偏移的原因。

如果想增加对结果的掌控感，可以用矩估计结果作为优化初值。Beta分布的矩估计公式是：

alpha0 = mu * (mu * (1-mu) / s^2 - 1) beta0 = (1-mu) * (mu * (1-mu) / s^2 - 1)

其中mu是样本均值，s^2是样本方差。用这个初值在手写负对数似然函数的fminsearch里跑一遍，通常两三次迭代就收敛。这个手写过程在调试时很有价值，因为betafit如果报错，你看不到内部迭代信息，而手写版本能逐步追踪参数轨迹，定位是数据边界问题还是初值问题。

2.4 完整可复制的拟合与可视化脚本

把上面片段拼装成一个完整脚本，输入两个CSV列，输出拟合参数、直方图对比、拟合优度检验结果。这是我从项目里截出来的简化版，可直接套用：

%% 读取数据 data = readtable('site_data.csv'); v = data.WindSpeed; % 风速，m/s G = data.GHI; % 水平面总辐照度，W/m2 %% 预处理 % 风速 v(v <= 0 | v > 60) = NaN; v = rmmissing(v); % 辐照度：只保留白天且物理合理范围 G(G <= 5 | G > 1400) = NaN; G = rmmissing(G); % 归一化 G_ref = quantile(G, 0.95); x = G / G_ref; x = min(max(x, 1e-6), 1 - 1e-6); %% 拟合 [k_wbl, c_wbl] = wblfit(v); [alpha_beta, beta_beta] = betafit(x); %% 直方图和拟合曲线 figure('Position', [100 100 800 300]); subplot(1, 2, 1); histogram(v, 'Normalization', 'pdf', 'NumBins', 40); hold on; v_grid = linspace(0, max(v), 200); plot(v_grid, wblpdf(v_grid, k_wbl, c_wbl), 'r-', 'LineWidth', 1.5); xlabel('风速 (m/s)'); ylabel('概率密度'); title(sprintf('Weibull: k=%.2f, c=%.2f', k_wbl, c_wbl)); subplot(1, 2, 2); histogram(x, 'Normalization', 'pdf', 'NumBins', 40); hold on; x_grid = linspace(1e-6, 1-1e-6, 200); plot(x_grid, betapdf(x_grid, alpha_beta, beta_beta), 'b-', 'LineWidth', 1.5); xlabel('归一化辐照度'); ylabel('概率密度'); title(sprintf('Beta: alpha=%.2f, beta=%.2f', alpha_beta, beta_beta)); %% 拟合优度检验 [h_wbl, p_wbl] = kstest(v, 'CDF', makedist('Weibull', 'a', c_wbl, 'b', k_wbl)); [h_beta, p_beta] = kstest(x, 'CDF', makedist('Beta', 'a', alpha_beta, 'b', beta_beta));

单独说明一个细节：Matlab的makedist里Weibull参数顺序是(a, b)，其中a对应尺度参数c，b对应形状参数k。这个对应关系非常容易搞反，我见过多个项目因为参数顺序写错导致KS检验结果一团糟。wblfit的返回值顺序是先形状后尺度，但makedist的输入顺序是先尺度后形状，两个方向不一样，这个坑必须记住。

3. 风光组合建模与互补性量化

3.1 联合分布框架：独立假设什么时候成立

把Weibull和Beta放在同一个框架里，最自然的思路是把风速v和归一化辐照度x看成二维随机变量，联合密度表达成两个边缘密度相乘：

f(v, x) = f_Weibull(v; k, c) * f_Beta(x; α, β)

这个相乘隐含了一个重要前提：风速和辐照度相互独立。但气象过程往往让它们存在相关性。最典型的是锋面过境时，大风和阴雨常常同时出现；夏季晴天午后，辐照度很高但大气相对稳定，风速反而偏低。如果无视这种相关性，组合模拟出来的风光联合出力可能过于乐观，因为模型会把高风和高辐照度同时发生的概率算大了。

我的工程处理方式是分层近似：先把数据按季节拆分，再按典型天气类型拆分，在每一个子集内做独立性近似。这样做的原因是，季节和天气类型这两个变量已经吸收了大部分风速-辐照度相关性，剩余残差相关性对结果的影响在工程误差范围内。如果你要更严谨，可以引入Copula连接函数去建模尾部相关性，但多数容量配置场景下，分层近似已经足够。

3.2 蒙特卡洛抽样与互补性指标定义

得到拟合参数后，组合研究最常用的是蒙特卡洛方法。反变换抽样的思路非常直接：rand生成[0,1]均匀随机数，再分别通过wblinv和betainv的逆CDF函数转成风速和归一化辐照度。代码只需要几行：

rng(2024); % 固定种子，保证可复现 N = 1e5; u_wbl = rand(N, 1); u_beta = rand(N, 1); v_sim = wblinv(u_wbl, c_wbl, k_wbl); x_sim = betainv(u_beta, alpha_beta, beta_beta);

抽样之后要转成功率。风电功率曲线一般给三段式或分段线性近似：低于切入风速、达到切出风速出力为零，中间基本按线性或二次曲线爬升。光伏出力则近似为辐照度的线性函数，乘上温度修正系数。这里给一个简化的功率转换：

% 简化的功率转换模型 v_cut_in = 3; v_r = 12; v_cut_out = 25; % m/s P_wind = zeros(N, 1); idx = v_sim >= v_cut_in & v_sim <= v_r; P_wind(idx) = (v_sim(idx) - v_cut_in) / (v_r - v_cut_in); idx = v_sim > v_r & v_sim <= v_cut_out; P_wind(idx) = 1; P_pv = x_sim; % 在简化模型里光伏出力近似正比于归一化辐照度

然后就可以算互补性指标。我常用的是一个简洁的归一化波动指标。设σ_wind和σ_pv分别为单独出力序列的标准差，σ_joint为风光联合出力序列的标准差，定义互补系数：

η = 1 - σ_joint / (σ_wind + σ_pv)

如果η接近0，说明联合出力和单独出力标准差之和不差多少，互补性弱；η越大，说明联合出力波动被显著压平，互补性强。这个指标虽然简单，但在方案比选时非常直观。做容量配置时，我会遍历不同风电/光电装机比，画出η曲线，找出让系统波动最小的比例区间，再结合成本模型选定最终方案。

3.3 概率约束下的容量配置逻辑

分布组合模型的真正价值，在于把"资源不确定性"变成"系统可靠性指标"。平均值方法算出来的电量往往很乐观，但极端天气段才是系统瓶颈。用蒙特卡洛样本，可以构造可靠性指标。

最常用的是LPSP，即电力不足概率。对每一组模拟的风光出力样本，判断是否满足负荷需求，统计满足不了的比例：

load_profile = 0.3; % 单位化负荷，MW battery_capacity = 0; % 简化模型里先不考虑储能 deficit = max(load_profile - (P_wind * cap_wind + P_pv * cap_pv), 0); LPSP = mean(deficit > 0);

把实际发电功率乘上对应的装机容量cap_wind和cap_pv，再和负荷比较，只要LPSP低于某一设计阈值比如5%，就认为该装机配比可接受。这样做比单纯看年发电量大数要靠谱得多，因为它显式地惩罚了资源差的那段时段。如果再叠加储能，还可以把LPSP继续压低，代价是储能容量成本上升。这个框架本质上是一个容量配置优化问题，而Weibull和Beta的拟合质量会直接影响所有下游结论，所以前面两步的扎实程度至关重要。

4. 实战踩坑记录与常见问题排查

4.1 六个高频问题速查表

做风光联合建模这几年，我把遇到过的典型问题整理成了下面这个速查表，每个问题都对应一个明确的修改建议：

4.2 参数敏感性经验

这类建模最容易忽视的是参数估计对拟合结果的敏感性。实测下来，Weibull的形状参数k每变化0.1，年发电量估算可能偏差3%到5%；Beta的α和β如果因为边界处理不当变化20%，联合出力的95%分位区间可能被整体拉偏。所以我的习惯是每次拟合都做一次参数稳健性检查：对原始样本做bootstrap重抽样，重复拟合一两百次，看看k、c、α、β的分布是否稳定。在Matlab里用bootci一行就能实现，这个验证成本很低，但能让下游结论可靠很多。

4.3 一张直方图之外的认知升级

最后提一个经验层面的建议。很多人做完拟合就收手了，只留下一张叠加了理论分布曲线的直方图。但我更推荐把参数做成滚动窗口曲线：按月滑动计算k、c、α、β，画在时间轴上。这样能直观看到风速稳定性在不同季节的迁移、辐照度Beta参数从夏天的高α低β模式切换到冬天的低α高β模式。这种图对项目汇报特别有说服力，也比单张直方图信息量大得多。

我在实际项目里发现，把滚动参数曲线放入评估报告后，新能源电站的投资方和电网调度人员都能很快理解资源特性，而不需要去读复杂的概率密度函数公式。这个附加产出的价值经常被低估，但它恰恰是建模工作从"技术结果"变成"决策依据"的关键一环。

结尾

做风光资源概率建模这几年，最深的体会是：Weibull和Beta这两个分布并不是工具箱里冷冰冰的数学函数，而是把复杂气象不确定性压缩成几个有物理含义参数的桥梁。k值记录着风速"乖不乖"，α和β则刻画出日照"稳不稳"。真正理解这层含义之后，你看到的不再是一条拟合曲线，而是这个场站一整年的天气故事。

最后再分享一个小技巧：项目交付时不要只给k、c、α、β这四个数，把按季度分组的拟合参数和滚动窗口趋势图一并附上。决策者真正关心的不是分布公式有多优美，而是资源波动在什么时段最剧烈、系统该在什么时候预留多大的裕度。把这些讲清楚了，你的建模工作才算真正闭环。