去年做风光互补电站的资源评估咨询时甲方提了个很朴素的需求按手头一年多的风资源和光照实测数据算一个最合理的光伏/风电配比。我一开始也走了很多工程报告里的老路——把年均风速、年均辐照度代进容量公式出一个看起来没毛病的方案。结果等拿到真正的小时级时序数据重新用概率分布做拟合之后才发现问题远没有均值算得那么简单。同样年均风速的两个风场一个可能常年刮稳定的小风另一个是台风过境式的大起大落年发电量能差出百分之二三十。这件事让我把风速用Weibull分布、光强用Beta分布这套组合建模方法彻底用到了实际项目上而且直接在Matlab里实现、检验、做置信评估。本文就把这套流程完整拆开从原理到代码到踩坑一次说清楚。适合正在做风光资源评估、容量配置或并网调度建模的工程师也适合科研党作为分布拟合的参考模版。1. 为什么要用双分布描述风光资源单一均值方案的局限1.1 平均风速的陷阱一次由均值算容量引发的误判先说一个我真实踩过的例子。两个拟建风场A场年均风速6.2 m/sB场年均风速6.1 m/s按很多教材里的简化公式估算年利用小时数两者几乎一模一样。可到了实测数据分析阶段我发现两个场的风速概率密度形态差异非常大A场的风速集中在5~8 m/s区间高风速段衰减平缓B场的风速却是大量2~4 m/s的弱风日加上偶尔几阵十几米每秒的大风。如果按同一台风机功率曲线去积分年发电量B场因为大量时间运行在切入风速以下出力表现明显差于A场。问题就出在均值不携带分布形态信息。风资源的随机性太强工程评估必须同时刻画平均大小和波动幅度、偏斜方向。Weibull分布正好提供两个参数形状参数k描述曲线形态尺度参数c描述整体大小两个参数就能把小风稳定型和大风间歇型区分开。这是用均值做单点评估永远做不到的。1.2 光强为什么落在Beta的边界里光伏资源评估面临类似问题。光照辐照度G的波动比风速更守规矩它天然被限制在0到某个上限之间白天从零升到峰值再降回零而且受云层遮挡影响会出现大量低辐照度样本导致日辐照度分布明显偏斜、有一侧长尾。如果我们硬套正态分布问题很明显正态分布定义域是整个实数轴而辐照度不会小于0也不会超过物理上限更重要的是晴天和阴天的数据混在一起时分布往往呈现高峰值正偏斜的形态正态分布的对称性根本拟合不了。Beta分布的定义域是[0,1]形状完全由两个正参数α和β控制既能拟合偏左偏右、高矮胖瘦的各种形态又不需要人为规定中心位置。所以把光强先归一化到[0,1]再用Beta拟合天然契合物理边界这比正态或对数正态都合理。1.3 双分布组合研究的意义用Weibull拟合风速、用Beta拟合光强并不是为了好看而是给后续的风光互补评估打地基。新能源出力评估的真正问题是给定全年的风速序列和辐照度序列我们能不能回答在90%置信水平下这个风光电站至少能出力多少这个问题的答案必须建立在两套概率分布之上再通过蒙特卡洛抽样或解析卷积去求组合出力的分位点。没有分布拟合就只能用确定性场景拍脑袋精度和可信度都无从谈起。2. Weibull分布拟合风速参数含义、估算方法与Matlab实现2.1 形状参数k与尺度参数c在风速场景下的物理含义Weibull分布的概率密度函数和累积分布函数分别是f(v) (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)F(v) 1 - exp(-(v/c)^k)参数c的量纲和风速一致m/s它的取值大致对应平均风速所在的量级更严谨地说它对应累积概率63.2%处的风速值。参数k无量纲直接决定曲线的形态k1时Weibull退化为指数分布风速集中在很低的值高风速段概率衰减很快k2时近似Rayleigh分布这是很多风资源报告默认采用的简化情形k在2.5~3.5时曲线接近钟形但仍有正偏比较接近典型内陆风场的实测风速分布k越大风速围绕某个中心值的集中程度越高。实际拟合中不推荐默认为k2因为沿海和内陆、山地和平原的差异非常大应该让数据自己说话。有一个工程经验值参考大多数陆地风场的k在1.5~2.8之间强季风地区可能到3以上。如果拟合出来k1.2要检查数据里是否混入了大量停机风速或静风时段。2.2 三种参数估算方法对比图形法、矩估计、极大似然方法原理优点缺点工程适用性图形法对CDF取双对数线性化拟合直线求k和c直观、手算可行人为分段影响大精度低只适合快速预览矩估计用样本均值和方差反推参数计算简单、不需要迭代高阶矩不稳定易受离群值影响可用于初值极大似然估计MLE最大化样本的对数似然函数统计性质好Matlab直接支持需要数值迭代数据量大时求和耗时最推荐缺省首选用MLE求Weibull参数表面上是调包但背后逻辑值得理解对n个独立风速样本v1, v2, ..., vn对数似然函数为logL n*log(k/c) (k-1)*Σlog(vi/c) - Σ(vi/c)^k对k和c求偏导并令其为零得到的方程组没有解析解需要迭代。Matlab的wblfit封装的就是这个迭代过程。如果只想要一个粗糙的手算初值矩估计公式也很实用c ≈ 均值 / Γ(1 1/k)其中Γ是伽马函数k初值可以先用变异系数近似。2.3 直接可跑的Matlab拟合代码下面这段代码先用真实参数生成一组模拟风速数据作为演示再分别用wblfit和手写MLE两种方式估算参数最后绘制拟合对比图。实际项目里把wblrnd那行换成读取实测风速数据即可。% 模拟一年8760小时风速数据(单位 m/s) rng(42); N 8760; true_k 2.2; true_c 7.5; windSpeed wblrnd(true_c, true_k, N, 1); % 方法一Matlab自带极大似然估计 phat wblfit(windSpeed); c_est phat(1); % 尺度参数 k_est phat(2); % 形状参数 fprintf(wblfit 估计: c %.3f, k %.3f\n, c_est, k_est); % 方法二手写MLE迭代用矩估计提供初值 meanWind mean(windSpeed); stdWind std(windSpeed); k_init (stdWind / meanWind)^(-1.086); c_init meanWind / gamma(1 1/k_init); k k_init; c c_init; tol 1e-8; maxIter 1000; for it 1:maxIter % 当前参数下的预测值 term1 mean(log(windSpeed) .* windSpeed.^k) / mean(windSpeed.^k); term2 mean(log(windSpeed)); k_new 1 / (term1 - term2); c_new (mean(windSpeed.^k_new))^(1/k_new); if abs(k_new - k) tol abs(c_new - c) tol k k_new; c c_new; break; end k k_new; c c_new; end fprintf(手写MLE 估计: c %.3f, k %.3f\n, c, k); % 绘图对比 edges 0:0.5:25; histogram(windSpeed, edges, Normalization, pdf, FaceAlpha, 0.3); hold on; v 0:0.1:25; plot(v, wblpdf(v, c_est, k_est), r-, LineWidth, 1.8); legend({风速直方图, Weibull拟合曲线}); xlabel(风速 (m/s)); ylabel(概率密度); title(风速Weibull拟合结果); grid on;这段代码里有两个细节值得提第一wblfit返回的参数顺序是[尺度, 形状]也就是[c, k]和很多教材里写(x; k, c)的顺序相反初学者常在这里颠倒导致画出的曲线完全变形我刚开始也栽过这个跟头。第二手写MLE的核心是交替迭代先固定k更新c再固定c更新k循环直到收敛。这个迭代方法在Matlab里运行几千个样本不费力但我实际项目中仍然以wblfit结果为准手写版更多用来对照验证。3. Beta分布拟合光照强度归一化、边界处理与手写MLE3.1 归一化到[0,1]区间的工程细节光资源数据通常来自气象站或卫星遥感反演单位是W/m²量级从0到1000甚至1100左右。用Beta拟合前必须归一化x G / G_max其中G_max的选取有两种做法一是取当地历史观测最大值二是取晴天理论峰值。用历史最大值的好处是简单坏处是如果观测期内没有出现真正的晴空峰值G_max偏小会把晴天样本推到接近1的位置导致α参数失真。我一般优先用当地纬度、日期对应的理论晴空辐照度上限比如1000~1100 W/m²从物理上保证归一化后的数据不会顶到1。还有一点很多人忽略Beta分布的坚实边界是开区间(0,1)而归一化后的辐照度很容易出现0这种值夜间、阴天极端情况。直接丢给betafit会报错或得到不合理的估计。常规做法是进行分类处理而不是简单粗暴地把0改成1e-6。我的处理策略是先把夜间和零辐照度样本剔除单独建立晴空/非晴空的0-1状态模型对白天有光照的样本做归一化Beta拟合如果非要保留零点可以给x设定一个下限例如x max(x, 1e-3)但这会轻微改变右偏分布的形状务必在报告中注明。3.2 为什么不用正态分布去拟合光强一次直观对比为了说明问题我拿一组模拟的白天辐照度样本分别做正态和Beta拟合。Beta拟合后的对数似然明显更高Q-Q图尾部也更贴合。根源在于辐照度分布有硬边界正态分布在左尾会预测出负值概率这在物理上不成立即使你忽略这个理论瑕疵实际拟合时左尾的偏差也会拖动均值和方差估计影响后续出力置信区间的精度。Beta的两个参数α、β可以分别控制左右尾的肥瘦程度把一个有界、偏斜的分布描述得特别准。3.3 从数据清洗到betafit、手写MLE的完整代码% 模拟白天光照数据(单位 W/m^2)峰值约1050 rng(123); N_day 5000; G_theoretical_max 1050; irradiance_raw betarnd(2.5, 4.0, N_day, 1) * G_theoretical_max; % 剔除“极端不可能”的0值保留正值 irradiance irradiance_raw(irradiance_raw 0); % 归一化(用物理上限而非样本最大值) x irradiance / G_theoretical_max; % 处理边界为避免恰好等于0或1造成betafit报错做轻微收缩 x(x 0) 1e-4; x(x 1) 1 - 1e-4; % 方法一Matlab自带betafit ab_est betafit(x); alpha_est ab_est(1); beta_est ab_est(2); fprintf(betafit 估计: alpha %.3f, beta %.3f\n, alpha_est, beta_est); % 方法二利用digamma函数手写MLE牛顿迭代(供验证) % 目标方程: psi(alpha)-psi(alphabeta)mean(log(x)) % psi(beta)-psi(alphabeta)mean(log(1-x)) mean_logx mean(log(x)); mean_log1x mean(log(1 - x)); alpha 1.0; beta 1.0; for it 1:500 % 计算当前残差 g1 psi(alpha) - psi(alpha beta) - mean_logx; g2 psi(beta) - psi(alpha beta) - mean_log1x; % 简易梯度下降更新(实际可用牛顿法改进) alpha_new alpha 0.01 * g1 * alpha; beta_new beta 0.01 * g2 * beta; if abs(alpha_new - alpha) 1e-8 abs(beta_new - beta) 1e-8 break; end alpha alpha_new; beta beta_new; end fprintf(手写MLE 估计: alpha %.3f, beta %.3f\n, alpha, beta); % 绘图 t linspace(0.001, 0.999, 200); histogram(x, 40, Normalization, pdf, FaceAlpha, 0.3); hold on; plot(t, betapdf(t, alpha_est, beta_est), b-, LineWidth, 1.8); xlabel(归一化辐照度 x G / G_{max}); ylabel(概率密度); title(光照Beta分布拟合结果); grid on;需要提醒的是手写MLE里我用了梯度下降的简化形式步长是经验取法实际做研究建议用更稳健的牛顿-拉夫森更新或者直接用betafit的返回值。我在项目里很少手写Beta的MLE因为betafit速度足够快而且对边界有内置处理。手写这版的主要价值是帮助理解参数估计的迭代本质后续做变体模型比如混合Beta分布时你就知道该在哪里扩展。4. 拟合优度检验KS检验、Q-Q图与看起来好但检验不过的原因4.1 KS检验的逻辑最大竖直线段差距拟合完参数不能拿图看一眼就宣布效果不错。图形直观但主观必须有一个量的判断。Kolmogorov-Smirnov检验KS检验是比较经验分布函数和理论分布函数的经典工具。它的统计量是D sup|F_n(x) - F(x)|其中F_n(x)是样本经验CDFF(x)是用拟合参数算出的理论CDF。D代表两条累积曲线之间最大的垂直距离D越大说明理论分布和实际数据偏离越严重对应的p值就越小。在常规显著性水平下比如α0.05如果p0.05就要拒绝数据服从该分布的原假设。有个容易误解的点KS检验的p值并不等于拟合是对的的概率。它只是在说在当前样本量下我们没有足够证据认为数据偏离该分布。样本量越大检验越挑剔哪怕偏差很小也会给出极小的p值。这就引出了下面要说的坑。4.2 Matlab实现kstest与关键参数设置Matlab的kstest默认检验标准正态分布要做Weibull或Beta的拟合优度检验必须手动传入CDF函数句柄或直接在kstest中用CDF参数指定。代码写法如下% 对风速做Weibull拟合优度检验 [cHat, kHat] wblfit(windSpeed); [hWind, pWind, ksStatWind] kstest(windSpeed, CDF, (v) wblcdf(v, cHat, kHat)); fprintf(风速KS检验: h%d, p%.4f, D%.4f\n, hWind, pWind, ksStatWind); % 对光强做Beta拟合优度检验 [aHat, bHat] betafit(x); [hSolar, pSolar, ksStatSolar] kstest(x, CDF, (g) betacdf(g, aHat, bHat)); fprintf(光强KS检验: h%d, p%.4f, D%.4f\n, hSolar, pSolar, ksStatSolar);这里要注意参数是从同一批数据里估计出来的再用这批数据做KS检验实际上会让检验偏保守p值偏大因为拟合已经吸收了一部分偏差。严格统计意义上应该做Lilliefors修正或bootstrap但工程场景下直接报告这一版结果也说得过去只需在方法描述中注明参数由MLE估计KS检验未做参数估计修正。4.3 看着好但检验不过的三种典型原因这种情况我在项目中遇到不止一次概率密度曲线叠在直方图上几乎完美重合但kstest给出的p值小于0.001。原因通常有三个。一是样本量太大。小时级数据一年有8760个点KS统计量的鉴别力随样本量明显提升常规检验几乎必然拒绝任何理论分布。这时候我会降低对p值的执念转而看D的量级只要D在0.02以下我就认为拟合在实际工程意义上是可接受的。二是数据不是纯单峰分布。比如风速数据混入了静风时段或者光照数据混合了晴天和阴天两种工况用单一个Weibull或Beta去拟合混合分布曲线插值后看起来平滑但CDF的坡度对不上。处理办法是先按工况分段建模白天/夜间、晴/阴不要试图用一个分布覆盖全部情形。三是尾部数据污染。异常值、传感器误码、停机检修记录混进样本哪怕只有千分之一的脏数据也会在经验CDF的尾部拉出一条台阶从而放大D值。把这些异常点清洗干净后检验往往就过了。5. 风光互补组合评估从两组拟合分布到容量置信区间5.1 组合研究到底在算什么拿到风速的(ksest, c)和光强的(aHat, bHat)工作远没结束这只是资源侧的概率建模。接下来要真正组合研究回答工程问题风、光同时波动时一个包含风机P_w和光伏阵列P_pv的混合电站整体出力在某置信水平下能保证多少这个问题不能简单地把各自均值加起来回答因为风、光并不是同时满发。常规做法是建立归一化出力模型。风机出力可以简化成三段/四段式曲线切入风速以下出力为0切入到额定风速之间按三次方或二次方升功率额定风速到切出风速之间满发切出后停机。光伏出力更简单可以近似为P_pv P_rated * x其中x就是归一化辐照度如果考虑温度效应还可以加一个温度修正项。把风速样本经风机出力曲线映射为风电出力把光强样本映射为光电出力再按装机配比加权就得到组合出力的样本序列。5.2 独立假设什么时候成立、什么时候不成立做组合评估前必须交代一个前提风速样本和光强样本的独立性。我通常先算风速与辐照度的相关系数。如果只有0.1以下独立抽样足够如果实测数据里明显存在白天光强高时风速偏低、夜间风速偏大的日间负相关仍硬着头皮假设独立会高估互补性导致置信区间过宽或过窄这在调度方案里会出大问题。严格的组合建模应该用Copula把两套边缘分布连接起来高斯Copula或t-Copula在风光联合建模里都很常见。但工程起步阶段先用独立假设做一个粗略评估再通过后面的敏感性分析看偏差有多大是一条务实的路径。5.3 蒙特卡洛抽样与置信区间估计的Matlab实现下面的代码演示了从拟合好的分布中生成组合出力的完整流程算出了不同置信水平下的可保证出力。% 假设风机额定功率1000kW光伏额定功率500kW风光配比2:1 P_wind_rated 1000; P_pv_rated 500; v_in 3; % 切入风速 m/s v_r 12; % 额定风速 m/s v_out 25; % 切出风速 m/s % 蒙特卡洛抽样次数 M 50000; % 从拟合Weibull抽风速 v_samples wblrnd(c_est, k_est, M, 1); % 从拟合Beta抽归一化光强(注意收缩到(0,1)避免极端尾) x_samples betarnd(alpha_est, beta_est, M, 1); x_samples min(max(x_samples, 1e-4), 1 - 1e-4); % 风机出力映射 P_wind zeros(M, 1); idx_power (v_samples v_in) (v_samples v_out); idx_full (v_samples v_r) (v_samples v_out); idx_ramp (v_samples v_in) (v_samples v_r); P_wind(idx_full) P_wind_rated; P_wind(idx_ramp) P_wind_rated * (v_samples(idx_ramp).^3 / v_r^3); % 光伏出力映射(忽略温度修正) P_pv P_pv_rated * x_samples; % 组合出力 P_total P_wind P_pv; % 统计指标 mean_total mean(P_total); std_total std(P_total); q [0.1, 0.05, 0.01]; conf_out prctile(P_total, q * 100); fprintf(组合出力均值: %.1f kW\n, mean_total); fprintf(90%%置信水平下最小出力: %.1f kW\n, conf_out(1)); fprintf(95%%置信水平下最小出力: %.1f kW\n, conf_out(2)); fprintf(99%%置信水平下最小出力: %.1f kW\n, conf_out(3)); % 绘制组合出力分布 histogram(P_total, 80, Normalization, pdf, FaceAlpha, 0.4); hold on; xline(conf_out(1), r--, LineWidth, 1.5); xline(conf_out(2), g--, LineWidth, 1.5); xline(conf_out(3), b--, LineWidth, 1.5); legend({组合出力概率密度, 90%置信水平, 95%置信水平, 99%置信水平}); xlabel(组合出力 (kW)); ylabel(概率密度); title(风光互补组合出力分布); grid on;这段代码跑完之后你会得到一组可以直接用于容量配置评估的数值比如总装机1500 kW在90%置信水平下可保证出力可能是320 kW而两个单一电站的90%置信出力之和却远小于这个数——这正是互补性带来的收益量化。值得注意的是如果要把这个结果用于并网调度或储能容量规划我通常还会把蒙特卡洛抽样次数提高到20万以上并把随机种子固定保证报告里的数字可复现。另外出力映射里的额定风速v_r选取要跟具体风机型号对应不要拿不同机型混用。6. 实操中的坑与经验总结6.1 数据清洗的优先级问题做分布拟合的第一道工序永远是数据清洗排在拟合之前。风速数据里的负值、平台静止记录长时间完全相同的数值、超过物理上限的异常点都会严重扭曲参数估计。光照数据里的夜间零值、日出日落过渡段的半阴影值也需要先剔除或单独建模。我的习惯是先把数据画成时序图用肉眼看一遍再用上下百分位数做初步过滤最后才进入拟合流程。顺序不能颠倒否则后面所有检验和置信区间都会被脏数据带跑。6.2 betafit和wblfit的收敛与参数初值问题Matlab的wblfit一般很稳但偶发遇到数据集中在极窄区间时迭代可能不收敛或给出极端值。这时我会手动加一个参数范围限制或者改用fminsearch直接最小化负对数似然。betafit的问题更多出现在数据边界大量样本贴近0或1时α、β估计会急剧膨胀甚至警告矩阵奇异。遇到这种情况先检查归一化时是否用错了上限再把极端边界样本做收缩处理就像前面代码里min/max那两行基本能解决。6.3 图表输出与报告呈现拟合结果最终要落在报告里我的输出标配是三张图概率密度拟合对比图、CDF对比图、Q-Q图。Q-Q图在工程报告中特别直观因为哪怕你写再多统计量非技术背景的评审还是习惯看散点是否落在对角线上我一般还会在图下注一句如果散点贴近yx说明实测分位点和理论分位点吻合良好。另外把拟合出的k、c、α、β参数连同95%置信区间一并写入报告比只放图更有说服力也方便后续复现。最后再分享一个实际体会分布拟合只是整个风光评估链条的地基它可以让你从拍脑袋配比升级到用置信水平说话但真正决定项目成败的往往是在数据清洗和工况界定上的细心程度。同一套Matlab代码给干净的数据和给脏数据跑出来的置信区间能差出三倍。所以如果你在自己的风光项目里也遇到拟合很漂亮但决策总出问题的情况建议先回头检查数据清洗流程再检查独立性假设是否成立。把这两关把好了剩下的统计计算交给Matlab它足够可靠。
