基于Copula与MBLS的光伏功率概率预测：从边缘分布到时空场景生成

这两年做光伏功率预测，接触了不少做电网调度和新能源场站的朋友。大家基本都有一个共识：光伏功率预测最难的不是把明天中午的出力值估准，而是把“可能出现的出力区间”估准。我见过很多点预测做得很漂亮的模型，一到实际调度就露馅——因为调度员真正需要的不是一个数，而是一个概率分布。这也是我今天要展开的Copula光伏功率预测思路的核心价值：前端用单调广义学习系统（MBLS）把条件均值回归做扎实，后端用Copula理论去描述时间和空间维度上的相依结构，最终输出的不是一条预测曲线，而是一批在时间上和空间上都互相协调的功率场景。这篇文章主要面向正在做光伏功率概率预测的研究生、算法工程师，以及需要把预测结果落到调度业务里的工程人员；后半部分的Matlab代码实现，是我在模拟项目中整理出来的可运行版本，照着自己的数据改一改，就能把整套流程搬走。

1. 光伏功率预测为什么绕不开“概率”这件事

1.1 一个具体的调度场景，点预测为什么会不够用

先讲一个我实际工作中遇到的场景。某天上午，调度员看到第二天的光伏功率点预测曲线，显示中午12点出力是180兆瓦。单看这个数字，火电、储能、负荷侧的安排都可以按这个值去匹配。结果第二天上午飘来一片云，12点实际出力只有90兆瓦，直接出现近百兆瓦的偏差，备用容量不得不在短时间内介入。

问题出在哪里？点预测把整个预测结果压成了一个期望值，但光伏出力的随机性远比“平均情况”复杂。晴天、多云、阴天、降水前兆，每种天气形态下，同一个时刻的出力可能落在完全不同的区间里。调度真正需要的，是“180兆瓦”背后的完整信息：这个值有多可信？上下浮动范围是多少？极端情况可能低到多少？换句话说，需要的是预测分布，而不是预测均值。

这就是概率预测的出发点。它不回答“明天中午出力是多少”，而是回答“明天中午每个出力值出现的概率有多大”。有了这个分布，调度员可以计算置信区间、评估风险、做鲁棒优化，甚至可以把预测分布直接馈入到随机规划框架里去。这比单纯追求点预测的RMSE下降，对实际业务的价值要大得多。

1.2 概率预测输出什么，时空预测又比单站时序多做了什么

概率预测常见输出有三种形式，在实际项目中按需求选：

• 分位数曲线：给出10%、50%、90%等分位点，直接画出预测区间，简单直观；
• 预测区间：两个分位点之间的区间，比如90%置信区间，业务侧最容易理解；
• 场景集合：一次性生成多条完整的时间序列轨迹，每条轨迹都满足时间和空间上的相关性，适合做随机优化、储能调度和风险评估。

单站时序概率预测通常只处理一个站点的时间不确定性，却忽略了空间上的联动。实际上，光伏电站之间因为云团移动、天气系统过境，出力曲线之间普遍存在时间滞后的相关性——A站先被云遮住半小时后，B站可能才被遮住。如果两个站点的预测结果各管各的，场景在空间上就是“散”的，整体可信度上不去。

时空概率预测要做的，就是同时把时间和空间两个维度的不确定性建模出来。而要把这两件事统一到同一个数学框架里，我的选择是用Copula理论：每个站点、每个时段的边缘分布单独建模，站点之间、时段之间的相依结构用Copula统一描述。这样既不过分简化问题，工程上又能拆成清晰的两步走。

2. 单调广义学习系统（MBLS）的原理与选取理由

2.1 先理解宽度学习系统：不加深反而加宽的回归底座

MBLS的前身是广义学习系统，它的思路和深度学习正好相反。深度学习依赖“深”，通过不断增加网络层数来提取抽象特征；广义学习系统依赖“宽”，把输入经过随机映射生成足够多的特征节点和增强节点，直接在宽度上把特征空间铺开，然后用岭回归一步求出输出权重。

整个过程可以这么理解：原始特征经过若干组随机权重变换，映射成特征节点组，再对特征节点做非线性变换得到增强节点，拼成一个大的特征矩阵，求解一个带正则项的线性回归。因为不需要反向传播、不需要迭代训练，广义学习系统在保证拟合能力的同时，训练速度和增量更新能力都非常突出。这在光伏预测场景里很关键：站点数据持续过来，模型需要频繁增量更新，每来一批新样本都从头训练一遍深度学习模型，计算成本实在太高。

2.2 单调性从哪里来：物理规律对模型的约束

光伏功率预测有一个非常明确的物理规律：在同一个站点、同一时刻的天气条件下，辐照度越强，出力越高，二者应当是单调非减的关系。普通机器学习模型不管这个规律，拟合出来的预测曲线偶尔会出现“辐照度变大了，出力反而变小”的反常情况。这类违背物理常识的预测，调度员一看就很难信任。

MBLS的核心改进，就是把这层单调性约束注入到广义学习系统中。实际实现上通常有两种思路：

• 限制权重符号：对于连接辐照度特征的那部分权值，约束为非负值，保证特征响应方向始终正向；
• 加入单调性惩罚：在训练损失里加一个正则项，对违反单调性的样本对进行惩罚。

我在模拟项目中用的是第一种思路的简化版本：对辐照相关特征通道做非负映射处理，其余特征通道保持自由。这种做法的好处是不改变模型的基本求解框架，依然可以闭式求解，只是对特征矩阵做符号约束，训练过程非常稳定。如果不做这个约束，预测曲线在云量突变时段经常出现局部抖动，物理上很难解释。

2.3 为什么MBLS适合做概率预测的前端模型

概率预测分两步：前端预测条件均值，后端对残差或预测分布建模。前端模型的选择会影响整个后端的可靠性。MBLS在这个位置上有几个明显的优势：

第一，训练稳定。岭回归求解输出权重，不存在梯度消失、学习率调参等问题，换数据集也不需要太多手动调参。第二，训练快，适合做滚动预测和增量更新。第三，拟合函数的平滑度可控，通过调节正则化系数，可以让预测曲线不会出现太多突兀的波动，这对残差分布的一致性很有帮助——残差分布越稳定，后端Copula建模越轻松。

但也要说清楚，MBLS不是万能的。它的特征节点和增强节点本质上是从随机映射之后的特征空间做线性组合，对极强非线性关系的表达能力不如深度网络。实际应用中，如果数据量巨大、天气形态极其复杂，可以把它作为基线模型，再考虑升级方案。

3. Copula理论怎么把“几个预测”变成“一组场景”

3.1 Copula的核心思想：把边缘分布和相依结构分开

Copula的数学定义并不复杂，核心是这样一个事实：任意一组随机变量的联合分布，都可以拆成两部分——每个变量各自的边缘分布，以及一个描述变量之间依赖关系的“连接函数”，这个连接函数就叫Copula。之所以要在光伏预测里引入它，是因为它可以让我们不直接对复杂的多维联合分布建模，而是先分别把每个预测时点和站点的不确定性搞准，再用Copula把相关性串起来。

打个比方，如果把联合分布比作一个家庭，边缘分布就是每个家庭成员的性格，Copula就是家庭成员之间的互动关系。我们可以先分别了解每个人的性格（边缘分布），再用一个关系模型去描述他们如何相互影响（Copula）。这样既降低了建模难度，又不会牺牲相关性信息。

3.2 边缘分布先行的建模流程

实际操作中，Copula光伏功率预测不是一步到位的，而是分阶段走：

• 第一步，用MBLS对每个站点、每个预测时段的归一化功率做条件均值预测，得到预测值；
• 第二步，分析预测误差的分布特征。误差可能是异方差的——晴天误差小，阴天误差大。所以不能直接套一个全局正态分布，我会根据天气类型或辐照度水平分组建模，或者用核密度估计拟合经验分布；
• 第三步，把每个站点的预测功率映射到[0,1]区间上的均匀分布。这个映射用的就是各站点自己的边缘累积分布函数；
• 第四步，再对这些均匀分布数据做Copula拟合，把站点之间、时段之间的相关性捕捉出来；
• 第五步，从拟合好的Copula中随机采样，得到一组[0,1]区间上的多维样本，再进行逆变换，还原成功率值。这组样本就是最终的空间相关场景集合。

这样一步步拆下来，每一步都有清晰的数学工具和评估方法，每个环节出了问题也能单独定位。

3.3 Copula族选型：高斯、t和Clayton怎么选

Matlab的数据工具箱里提供了协方差估计和高斯Copula、t Copula、Clayton Copula等常用类型。我实际测试下来，选型经验如下：

实际项目里我一般先用高斯Copula做基线，再用t Copula做对比，看对数似然和场景评价指标再定。如果站点之间存在明显的厚尾联动，比如强对流天气导致出力同时骤降，t Copula往往有明显优势。

3.4 时空相关性如何被Copula编码进去

有人会问，Copula不是只描述变量之间的静态相关性吗，时间动态怎么处理？这个问题在实现上需要一点技巧。我会把预测时段的维度也展平成变量维度的一部分。比如有5个站点，预测未来4个小时、每15分钟一个点，那就是16个时段，展平后一共80个维度，对这80维的联合分布直接做Copula建模。这样时间相关性、空间相关性、时空交叉相关性全部包含在同一个相关矩阵或Copula结构里。

当然，维度升高会带来相关矩阵估计的方差问题。折中的办法是分两步：先做时间维度上的自相关建模，比较常用的做法是用多维高斯Copula或ARMA类结构近似；再做空间维度上的Copula拟合。或者干脆对展平矩阵直接做主成分压缩，降维之后再拟合Copula。简单场景下，直接展平拟合高斯Copula在小规模站点数下表现已经足够。

4. 时空概率预测整体框架：数据怎么来、特征怎么管

4.1 数据集结构与预处理要点

我做模拟项目时，用的是几座相邻光伏电站的历史数据，每条记录至少包含以下字段：

数据预处理阶段最容易踩的坑是时间对齐。各站点的数据可能因为采集设备差异出现时间戳错位，如果不先按照统一时间粒度重采样，后面算空间相关性时会出现虚假的滞后或超前。我通常先把所有时序数据统一重采样到15分钟网格，再检查缺失值比例，超过5%的时段直接剔除，少量缺失用线性插值补。

4.2 特征设计：时序、空间滞后和天气分型一个都不能少

特征设计直接决定MBLS的上限。我常用的特征组包括：

• 时序特征：当前时刻往前推4到8个时刻的出力、辐照度序列，让模型感知数据趋势；
• 空间滞后特征：相邻站点的前几个时刻出力，尤其是上游风向一侧的站点滞后特征，对下游站点预测很有帮助；
• 气象特征：预报辐照度、温度、湿度等，有预报数据一定要用，模型才能对未来的天气突变有反应；
• 天气类型编码：把天气类型做成分类变量，让模型针对不同天气形态学习不同的映射关系。

空间滞后特征是时空预测的关键变量。因为云团移动方向通常与风向一致，上游站点的出力变化往往领先于下游站点，把这种滞后关系学进MBLS，比单纯堆时间序列特征效果要好很多。

4.3 训练集和验证集怎么切

时间序列数据严禁随机打乱切分。光伏预测模型如果训练集和测试集时间重叠，会严重高估模型能力。我采用滚动切分方案：用过去30天数据训练，预测未来一天，然后滑动窗口不断向前推进。验证集必须是模型完全没有见过的时间段，最好是跨越不同天气季节的连续时间窗口，这样评价出来的指标才可靠。

5. Matlab代码实现：从训练MBLS到生成Copula场景

5.1 数据读取与特征构造代码骨架

Matlab实现的第一步是组织数据。下面这段代码是数据读取和特征构造的骨架，可以直接替换成自己的数据源。

% 读取多站点数据 % data: 表格，列依次为timestamp, site1_power, site2_power, ghi, temp, wind data = readtable('pv_data.csv'); data.timestamp = datetime(data.timestamp, 'InputFormat', 'yyyy-MM-dd HH:mm:ss'); % 站点功率列 power_cols = {'site1_power', 'site2_power'}; n_site = length(power_cols); % 归一化功率 for k = 1:n_site data.(power_cols{k}) = data.(power_cols{k}) / rated_power(k); end % 构造时序特征：前p个时刻的历史功率和辐照度 p = 6; feat = []; targets = []; for k = 1:n_site for lag = 1:p feat_name = sprintf('site%d_lag%d', k, lag); data.(feat_name) = [NaN; data.(power_cols{k})(1:end-1)]; end end % 剔除开头NaN行 data = rmmissing(data);

这段代码的思路是先把功率归一化到0-1区间，然后为每个站点构造滞后特征。滞后阶数p需要根据数据颗粒度和预测时效调整，15分钟数据一般取6到8个滞后点比较稳妥。

5.2 MBLS训练与单调性实现

MBLS的核心是特征节点、增强节点和岭回归求解。Matlab实现的关键如下：

% 训练集和测试集划分 X_train = train_features; % n x d 矩阵 Y_train = train_targets; % n x 1 向量 % 特征节点映射 n_feature = 50; W_f = randn(size(X_train, 2), n_feature) / sqrt(size(X_train, 2)); b_f = rand(1, n_feature) * 0.1; Z = tanh(X_train * W_f + b_f); % 增强节点映射 n_enhance = 100; W_e = randn(n_feature, n_enhance) / sqrt(n_feature); b_e = rand(1, n_enhance) * 0.1; H = tanh(Z * W_e + b_e); A = [Z, H]; % 岭回归求解输出权重 lambda = 0.1; W_out = (A' * A + lambda * eye(size(A, 2))) \ A' * Y_train; % 单调性处理：对辐照度特征通道做非负权重修正 % 这里是关键技巧，把连接辐照度特征的那组权重限制为非负 gamma_idx = find(contains(feature_names, 'ghi')); W_out(gamma_idx, :) = max(W_out(gamma_idx, :), 0);

窄看这段代码会发现，我在岭回归求解之后又做了一步“权重非负修正”，这是单调性约束的工程近似处理。对于小规模数据，效果比较直观；如果数据规模变大、非线性增强，可以考虑把非负约束放进优化目标内部，用带约束的最小二乘求解，Matlab的lsqlin函数可以直接支持。

5.3 Copula参数估计与时空场景采样

这一步是整套模型的精华。边缘分布+相依结构分开建模完了之后，用Copula采样生成时空场景。

% 假设pred_power是n_steps x n_site的MBLS预测功率 % 假设residual是残差矩阵，每列对应一个站点 % 第一步：对每个站点、每个时段拟合残差分布（这里用简单高斯核密度） n_steps = 16; n_site = 2; margin_cdf = cell(n_steps, n_site); for s = 1:n_site for t = 1:n_steps r = residual_history{t, s}; % 用核密度估计累积分布 pd = fitdist(r(:), 'Kernel', 'Kernel', 'normal'); margin_cdf{t, s} = pd; end end % 第二步：将历史残差转为均匀分布序列，用于Copula拟合 U = zeros(length(residual_history), n_steps * n_site); for i = 1:length(residual_history) k = 1; for s = 1:n_site for t = 1:n_steps U(i, k) = cdf(margin_cdf{t, s}, residual_history{i, t, s}); k = k + 1; end end end % 防止0和1的极端值影响估计 U = min(max(U, 1e-6), 1 - 1e-6); % 第三步：拟合高斯Copula [rho, nu] = copulafit('Gaussian', U); % nu在Gaussian下无效，统一接口需要 % 如果改用t Copula，使用： % [rho, nu] = copulafit('t', U, 'Method', 'ApproximateML'); % 第四步：生成M个随机场景 M = 500; U_sim = copularnd('Gaussian', rho, M); % 生成M行、n_steps*n_site列的[0,1]样本 % t Copula对应：U_sim = copularnd('t', rho, nu, M); % 第五步：逆变换回功率残差空间，叠加MBLS预测均值 scenarios = zeros(M, n_steps, n_site); for m = 1:M k = 1; for s = 1:n_site for t = 1:n_steps res_sample = icdf(margin_cdf{t, s}, U_sim(m, k)); scenarios(m, t, s) = pred_power(t, s) + res_sample; k = k + 1; end end end % 截断到[0,1]（功率已归一化） scenarios = min(max(scenarios, 0), 1);

这段代码的注释里已经包含了高斯Copula和t Copula两种切换方式。实际应用时，先用高斯Copula跑通流程，再换t Copula做对比即可。

5.4 概率性评价指标与可视化思路

场景集合生成之后，必须用标准指标量化评价。我常用的指标有两个：

• CRPS（连续排名概率得分）：衡量预测分布与真实值之间的差距，同时反映准确性和锐度，值越小越好；
• PIT（概率积分变换）直方图：检查概率预测是否校准。如果直方图近似均匀分布，说明预测分布与真实分布匹配良好；如果是U型或倒U型，说明分布过窄或过宽。

Matlab里CRPS没有现成函数，自己实现也不算麻烦。核心逻辑是把预测ECDF和真实观测的阶跃分布函数之差做积分，实际计算时可以转成对场景集合的直接计算：对每个时刻，统计场景集合中小于观测值的比例，结合场景分布宽度综合算出CRPS分数。

可视化方面，我会画三张图：第一张是某站点未来16个时段的场景集合图，用半透明色带显示规模区域；第二张是不同分位数（P10、P50、P90）与真实值的对照曲线；第三张是两个站点在同一天的场景散点图，直接展示时空相关性是否被Copula有效捕获。

6. 实测中的关键坑位与参数调优经验

6.1 边缘分布拟合不当是最大的精度杀手

Copula建模的前提是边缘分布准确。我在测试中发现，残差分布直接用全局高斯拟合时，场景集合很容易在极端出力处失真。光伏功率残差往往有厚尾，尤其是阴天和晴天交界时段。第一个坑是核密度估计的带宽需要按站点调节，带宽过小会让PIT直方图出现尖峰，带宽过大会让场景集合过于平坦。第二个坑是残差分布可能随时间变化，所以我会按天气类型分组拟合边缘分布，并且用滑动时间窗口定期更新，而不是一次性拟合所有历史数据。

6.2 空间相关矩阵的估计窗口不是越长越好

有些朋友一开始会拿一整年的数据去估计站点间的相关系数，结果模型在线运行后场景相关性反而失真。原因在于光伏站点的空间相关性受季节和天气系统影响显著，春夏和秋冬的云团形态、主导风向都不一样。我实测下来的经验是：用最近15到30天的数据滚动估计空间相关矩阵，并每3到5天更新一次，效果明显优于用全年固定矩阵。相关性矩阵在阴雨季节和晴朗季节之间切换时，站点相关性的变化幅度超出预期，固定矩阵完全吃不到这个信息。

6.3 单调性约束的强度必须交给验证集决定

单调性约束加得太强，模型会变得过于保守，预测曲线虽然物理上合理，但细节拟合不足，尤其在云量突变时预测值可能被“顶”住；约束太弱，又会出现局部非物理抖动。实操上，我把约束强度作为一个超参数，在验证集上扫描三个档位：

验证集上同时看CRPS和物理矛盾样本比例。选出CRPS最低且物理矛盾比例很低的档位，再进入正式测试集评估。

6.4 场景数量要够，但也不是越多越好

场景生成数量在500到1000个左右已经能覆盖大多数调度需求。场景太少，极端出力区间的统计会不稳定；场景太多，计算和存储成本增加，但指标改善非常有限。我常用500个场景做在线预测，2000个场景做离线评估，这算是一个比较稳妥的实践标准。

6.5 极端天气时段的代表性问题

最后分享一个容易被忽略的点：不管Copula拟合得再好，训练数据里极端天气事件始终是少数样本。对电网调度来说，恰恰是强对流、台风边缘、连续阴雨这类场景风险最大。我的建议是，为极端天气单独建立小样本故障场景库，在常规概率场景之外额外叠加少数极端场景，让调度员能看到“最坏情况”的边界。这个方法虽然技术含量不高，但在实际业务对接中评价非常高。

整套模型跑下来，最深的体会是：Copula光伏功率预测本质上不是某一个模型有多强，而是把“回归精度”和“相依结构”两个问题拆开、分别做到位。MBLS负责前者，Copula负责后者，各司其职，模型的可解释性也大大提升。如果你正在做多站点光伏概率预测，不妨按这个框架先跑通一遍，再逐步替换组件。