做风光资源评估的人十有八九都干过这么一件事拿一年测风数据点开Matlabhistogram一下再拿wblfit拟合一条Weibull曲线叠上去。光伏那边也类似把归一化后的辐照度数据用betafit拟合一条Beta分布曲线。这两个分布一个管风电一个管光电是新能源概率建模里用得最多的一对搭档。但多数教程只告诉你要用这个函数、那个工具箱不讲清楚为什么偏偏是这两种分布参数背后对应什么物理过程更不会讲怎么把Weibull和Beta模型组合起来系统性地研究风光互补特性。这篇文章就把这条完整链路捋一遍先从数学和物理背景讲清楚选型逻辑再给出可直接复制的Matlab代码包括参数估计、拟合优度检验、蒙特卡洛组合抽样最后放上我实际项目里踩过的坑和排查方法。不管你是做风电场测风数据分析、光伏电站辐照度建模还是微电网容量规划这套流程都能拿来直接用。1. 为什么偏偏是Weibull和Beta1.1 Weibull分布不是巧合而是风速统计的必然很多人第一次接触Weibull分布是在风资源评估报告里往往只看结论平均风速多少、Weibull形状参数k和尺度参数c是多少。但真正要理解的是为什么风速数据会呈现出典型的右偏分布。风速的物理来源是气压梯度力但地表的摩擦拖拽、地形抬升、湍流交换会把这股规则的气流搅得支离破碎。实测风速从来不是对称的钟形曲线而是拖着一条长长的右尾也就是高风速时段虽然少但确实存在。同时风速不可能为负这就把分布严格限制在半轴上。Weibull密度函数恰好具备这两个特征定义域在[0,∞)形状参数k能灵活调节偏度。k1时退化为指数分布k2时就是瑞利分布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 两款分布选型的对照逻辑用一个表格把这套选型逻辑收拢一下方便后续理解代码设计对比维度Weibull分布Beta分布适用变量风速归一化辐照度定义域[0, ∞)(0, 1)典型形态右偏、可近似对称U形、钟形、偏态均可核心参数k形状、c尺度α形状、β形状物理含义k反映风速稳定性c和平均风速正相关α、β共同刻画晴空比例和波动强度估计复杂度极大似然有迭代解析格式极大似然需数值优化常见失误零风速样本过多导致k偏低数据含边界值导致拟合失败选型不当的后果我在项目里见过不少。有人硬用正态分布拟合风速拟合出来的分位数在低风速段直接变成负值后续做储能容量配置时算出荒谬结果也有人把夜间辐照度一起塞进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(11/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对应尺度参数cb对应形状参数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 六个高频问题速查表做风光联合建模这几年我把遇到过的典型问题整理成了下面这个速查表每个问题都对应一个明确的修改建议现象可能原因处理办法betafit返回NaN或负数参数数据包含0或1Beta开区间被破坏做边界偏移xmin(max(x,1e-6),1-1e-6)只保留白天样本Weibull的k值明显偏小数据里静风零值太多标准Weibull被拖低用零截断Weibull或对v0子集建模并记录零风概率kstest在图形拟合良好时仍然拒绝原假设样本量大微小偏差也会被检出综合看QQ图和K-S距离不要只盯着p值全年拟合效果差但分季节拟合效果好季节异质性被混入同一组参数按春/夏/秋/冬或按月分组分别拟合归一化辐照度后Beta曲线头部严重失真G_ref选成了历史最大值样本挤压在低值区改用95%分位数作为G_ref蒙特卡洛模拟两次结果完全不同没有固定随机数种子在脚本开头设置rng(可复现的固定数)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、α、β这四个数把按季度分组的拟合参数和滚动窗口趋势图一并附上。决策者真正关心的不是分布公式有多优美而是资源波动在什么时段最剧烈、系统该在什么时候预留多大的裕度。把这些讲清楚了你的建模工作才算真正闭环。