把风光的随机性建模成一组“有代表性、可计算”的离散场景，是电力系统随机优化里绕不开的一步。我最早接触这个课题时，第一反应是把风速和辐照度当作两个独立变量分别采样，结果做出来的联合场景严重失真：极端天气下的风光互补风险完全没被捕捉到。后来把分布之间的关联用 Copula 显式建模，再用 K-means 把几千个场景削减成几十个，才算真正把“基于 Copula 理论与 K-means 的风光场景生成与削减”这套流程跑通。这篇就围绕这个主题，把我在实际项目里验证过的建模细节、代码实现和踩坑记录整理出来。
适合谁看：准备做微电网容量优化、储能配置、日前调度随机规划的工程师和研究生，以及想用场景法但被大量文献术语劝退的入门者。我会尽量少用公式黑话，多用实际例子说话。
1. 场景生成与削减到底是什么
1.1 为什么不能直接用原始历史数据
先把背景说清楚。风光出力是典型的不确定变量，历史数据只是这条随机过程的一条现实轨迹。随机优化里需要把未来的不确定性描述成“多个等概率且不相容的场景”，这样优化模型才能看到最坏情况和平均情况。举个例子，做微网储能容量配置时，如果只拿一年的8760小时历史数据去跑优化，结果是“只对这一年负责”，第二年来了个极端天气就直接失守；而场景法则会把风速和辐照度采样出成千上万条可能曲线，然后从中挑选代表性子集，让规划结果对“统计意义上的未来”负责。
直接用原始数据还有三个麻烦。第一是数量太大，8760个小时样本做机组组合就已经吃力，更不要说在样本外继续扩充；第二是极端组合覆盖不全，历史里可能压根没有“高风速+强辐照同时出现”的样本，但 Copula 可以把这种联合分布的尾部外推出来；第三是历史曲线之间存在强自相关，直接把每小时样本打乱重排会破坏时间连续性。所以标准做法是分两步走：先生成足够大的人工场景集，再把它们削减到优化模型能处理的数量。
1.2 为什么选 Copula + K-means 这套组合
场景生成方法其实不少，例如独立采样、马尔可夫链、Copula 采样，还有近年很火的生成对抗网络。独立采样的问题很明显：风速和辐照度实际上不是独立的，阴天通常伴随大风，晴天往往风小，这种负相关性直接影响风光互补容量。马尔可夫链能抓住时序转移特性，但状态空间一大会出现组合爆炸，两个变量、每变量分成10个状态就要100维转移矩阵。深度学习生成模型效果好，但对数据量、训练时间和调参要求高，在工程上不容易落地。
Copula 的优势在于把“单变量的边际分布”和“变量之间的依赖结构”分开建模，这在统计上是经过严格证明的思路。它允许风速继续用 Weibull 分布、辐照度用 Beta 分布，再把两者的相关性用一个 Copula 函数连接起来，灵活且可解释。K-means 作为削减工具则胜在简单稳定、计算快、可复现，对成千上万个场景做聚类压缩，聚类中心就是典型场景，簇内样本个数就是概率权重，非常好用。这套组合的计算成本基本控制在分钟级，拿来做方案比选和参数敏感性分析非常顺手。
2. Copula 建模：为风光相关性建立数据基础
2.1 数据准备中的关键决策
数据是建模的地基。我用的数据来自某地区实测气象站的逐小时风速和辐照度记录，时间跨度为三年，共约 26280 个时刻点。拿到数据第一件事不是建模型，而是清洗，这部分花费的时间往往比想象中多。
第一步是处理缺失值。气象站偶尔会出现连续几个小时的设备故障，直接删掉会导致时间序列断裂，插值又会引入虚假相关性。我的处理方式是：缺失长度小于三小时的用前后线性插值，超过三小时的直接剔除当日数据，因为长时间缺失段严重时会造成伪趋势。第二步是处理异常值，比如凌晨两点辐照度突然出现 800 W/m²，这种显然不是数据记录错误就是传感器问题，需剔除。第三步是时间对齐，风速和辐照度必须精确到同一小时，否则相关性计算会失真。
预处理完后还有个容易被忽略的环节：数据分型。夜间辐照度是 0，白天才有值，如果把全天 24 小时混在一起拟合分布，会产生大量零值堆积。我的做法是把场景建模对象分成两段：白天段（比如 6 点到 18 点）用 Beta 分布拟合辐照度，夜间段单独处理或直接令出力为零；风速则不分昼夜统一建模。这个细节直接决定后面生成的场景是否合理。
2.2 边际分布拟合怎么选
边际分布是 Copula 建模的前置步骤。Copula 只处理变量之间的相关性结构，每个变量自身的概率形态要由边际分布负责。风速用 Weibull 分布是最经典的选择，因为它的双参数（形状参数 k、尺度参数 λ）能很好地刻画风速右偏分布特性，公式也直观：PDF 为 f(v) = (k/λ)(v/λ)^(k-1) exp(-(v/λ)^k)。拟合用极大似然估计，Python 里scipy.stats.weibull_min.fit直接能出参数。
辐照度稍微麻烦点。晴空条件下辐照度近似 Beta 分布，但全天数据里有大量零值（夜间和阴天），直接把零值和非零值混在一起拟合 Beta 会让参数严重偏掉。我在实操里用了一个很简单的“混合模型”思路：P(辐照度=0) 用一个伯努利变量刻画，非零部分用 Beta 分布拟合。生成场景时先以概率 p 抽出“是否为晴天段”，再对非零段从 Beta 分布采样。
拟合完成后必须做一次 QQ 图检验。把理论分布的分位数和实测数据的分位数画在一起，如果点近似落在 45 度直线上，说明边际拟合合格。这一步很多人跳过，直接进 Copula 拟合，结果相关性参数全是歪的。我记得某次检验发现风速的尾部偏离严重，后来给 Weibull 加了个零风速概率项（混合 Weibull），效果才达标。
2.3 Copula 函数的理解与选择
现在到了核心环节。简单说，Copula 就是把多个变量的联合分布函数拆成两部分：每个变量的边际分布函数，加上一个“连接函数”来描述变量间的关系。Sklar 定理是这一切的基石：如果 F 是联合 CDF，F1、F2 是边际 CDF，那么必然存在一个 Copula C，使得 F(x1, x2) = C(F1(x1), F2(x2))。反过来说，只要给定边际分布和 Copula，就能唯一确定联合分布。
风光场景建模中，我第一推荐的是 Gaussian Copula（高斯 Copula）。它的表达式为 C(u1, u2) = Φ_ρ(Φ⁻¹(u1), Φ⁻¹(u2))，其中 Φ 是标准正态分布 CDF，Φ⁻¹ 是逆 CDF，ρ 是标准正态空间里的相关系数。它相当于把两个均匀变量先映射到正态空间，再按相关系数 ρ 做二维正态采样。Gaussian Copula 胜在计算简单、参数只有 ρ 一个、拟合稳定，适合大多数风光联合建模场景。
但也有需要谨慎的场景。如果数据表现出明显的尾部相关性，比如极端大风天往往伴随极端低辐照，Gaussian Copula 就会低估这种尾部联动风险。这时可以尝试 t-Copula（t-Copula），它在尾部比 Gaussian Copula 更厚，能捕捉到极端的联合失效概率，代价是多了个自由度参数要拟合。更复杂的 Gumbel 和 Clayton Copula 属于阿基米德族，分别擅长刻画上尾相关和下尾相关，但参数解释性不如 Gaussian 直观，我在工程上很少用。
如何选？我的经验是先算数据在正态空间下的 Spearman 秩相关系数，用scipy.stats.spearmanr能得到一个稳健的 ρ 初值；再分别拟合 Gaussian 和 t-Copula，用 AIC 比较拟合优度；最后把两种 Copula 采样生成的相关性与经验分布做核对。绝大多数情况下 Gaussian Copula 的 AIC 已经足够低，t-Copula 只有在自由度小于 15 时才有实质性提升，否则没必要自找麻烦。
3. 场景生成实操：从 Copula 采样到完整场景
3.1 蒙特卡洛采样原理与步骤
Copula 建好后，生成场景的过程就是一组纯采样操作。以风速 v 和辐照度 s 两个变量为例，流程分五步：
第一步，取历史风速数据并拟合 Weibull 分布，得到 F_v(v)；取非零辐照度数据并拟合 Beta 分布，得到 F_s(s)，同时统计零辐照度占比 p0。第二步，利用历史观测值计算 Spearman 秩相关系数 ρ，构造 Gaussian Copula 的相关矩阵。第三步，从二维标准正态分布 N(0, Σ) 采样出一大堆 (z1, z2)，其中相关系数为 ρ；然后做变换 u1 = Φ(z1)、u2 = Φ(z2)，得到均匀分布空间里带相关性的样本对。第四步，对 u2 表示辐照度均匀样本做零值还原：以概率 p0 令辐照度为 0，否则采样 Beta 逆 CDF。第五步，用逆变换把 u1 映射回风速实际值，非零的 u2 映射回辐照度实际值，最终得到一组模拟的“风速-辐照度”联合样本。
这里的关键点是“逆变换”必须严格对应边际分布。如果边际拟合的是 Weibull，逆变换就是v = weibull_min.ppf(u, k, loc, scale)。很多初学者在第四、五步直接套用正态逆变换，生成出来数值范围完全对不上，这就是典型的没有理解 Copula 的边际分离思想。
3.2 时序相关性增强
上面描述的是静态采样，每次只生成一个风速和辐照度时刻。但实际场景是一个时间序列，哪怕静态相关是对的，时序上也不能一天从早到晚都是同一个值。最直接的方式是在 24 小时的维度上，用每个小时的边际分布分别采样，再逐小时用 Gaussian Copula 串联。这样一天的场景就是一串“风速-辐照度”对，每个小时内部满足相同的相关结构。
但这样还不够，因为相邻小时的风速往往有很强的自相关。某天风速 10 m/s，下一小时突然掉到 2 m/s 的概率并不高；纯独立采样会让时间曲线剧烈震荡，看起来像噪声而非真实天气过程。我通常会在 Copula 采样后加一个自回归平滑环节：用一阶自回归模型 AR(1) 对风速序列做滤波，比如 v_t = φ·v_(t-1) + (1-φ)·v̄ + ε_t，其中 φ 从历史数据自相关函数估计，通常在 0.8~0.9 之间。辐照度的时间相关性主要靠昼夜周期天然约束，夜间必然为零，白天可由云层遮蔽的块状过程产生，用高斯核做一个小时尺度的平滑即可，不用上复杂模型。
对于只想先跑通流程的读者，我的建议是：先把静态 Copula 采样和削减做出来，时序平滑作为二期优化。因为时序增强会引入参数，参数越多，验证工作量越大。
3.3 Python 实现代码示例
下面给出一段我实际在用的核心代码，基于scipy.stats、numpy和copulas库（Copulas 库由数据科学团队维护，接口直观），演示从边际拟合到 Copula 采样再到逆变换的完整链路：
import numpy as np from scipy import stats from copulas.multivariate import GaussianMultivariate # 1. 边际拟合 params_wind = stats.weibull_min.fit(wind_data, floc=0) cdf_wind = lambda x: stats.weibull_min.cdf(x, *params_wind) ppf_wind = lambda u: stats.weibull_min.ppf(u, *params_wind) # 非零辐照度 params_solar = stats.beta.fit(solar_data[solar_data > 0]) zero_prob = (solar_data == 0).mean() # 零占比 cdf_solar = lambda x: stats.beta.cdf(x, *params_solar) ppf_solar = lambda u: stats.beta.ppf(u, *params_solar) # 2. 把边际转换成均匀样本 u_wind = cdf_wind(wind_data) u_solar = cdf_solar(solar_data) # 只对非零部分计算 # 3. 拟合 Gaussian Copula copula_model = GaussianMultivariate() copula_model.fit(np.column_stack([u_wind, u_solar])) # 4. 从 Copula 采样生成 10000 个场景 samples = copula_model.sample(10000) u1_sample, u2_sample = samples[:, 0], samples[:, 1] # 5. 逆变换回物理量 wind_scenarios = ppf_wind(u1_sample) solar_scenarios = np.where(u2_sample < zero_prob, 0, ppf_solar(np.clip(u2_sample, 1e-6, 1-1e-6)))注意几个细节：weibull_min.fit的floc=0表示把位置参数固定为 0，因为风速不能为负；np.clip用来防止逆变换时出现概率值恰好等于 0 或 1 导致溢出；对辐照度零值的判断必须在逆变换之前完成。
3.4 场景标准化与格式组织
生成完原始场景后，我习惯统一组织成标准矩阵格式：每行是一个场景，每列是某个时刻点的风速或辐照度。假设生成 10000 个场景、每个场景 24 小时，那么矩阵形状就是 (10000, 48)，前 24 列风速、后 24 列辐照度。这个格式后面喂给 K-means 非常方便。
在保存场景之前还有一个标准化问题：风速的量纲是 m/s，辐照度的量纲是 W/m²，如果直接把原始值合并进特征矩阵，K-means 的距离计算会天然偏向幅值大的辐照度，风速的信息被稀释。我通常用StandardScaler先做 z-score 标准化（减均值除标准差），聚类完成后再把聚类中心还原到原始量纲。这一步在压缩场景时是必须的，很多人漏掉后聚类结果全是辐照度主导，风速偏差完全不体现。
4. K-means 场景削减：把 10000 个场景压缩到 30 个
4.1 削减的必要性与数学目标
生成几千上万个场景是为覆盖不确定性，可优化模型承受不了这么大的输入。以我们的微网规划模型为例，每个场景会引入一组连续性变量和约束，10000 个场景直接让混合整数线性规划的求解时间从几分钟变成几小时，甚至在内存上就过不去。所以必须削减。
削减不是随机抽几个场景，而是要求削减后的 K 个场景在某种概率距离下尽量逼近原始场景集。最常用的是 Wasserstein 距离：它度量两个概率分布之间的“最小搬运代价”，比直接比较均值方差更全面。但算真正的 Wasserstein 距离复杂度高，所以工程上通常用 K-means 这种先聚类、后取聚类中心的办法逼近，聚类中心作为代表性场景，簇内样本数量占比作为概率权重。
这样一来，优化模型里的一个不确定性集合就变成了：K 个典型日 + 每个典型日出现的概率。这正是随机规划标准输入格式。削减的数学目标可以理解为最小化“每个原始场景到其所属聚类中心的距离的加权平均”，K-means 的代价函数恰好就是它。
4.2 K-means 聚类削减具体操作
我的标准流程如下。
先把标准化后的场景矩阵输入sklearn.cluster.KMeans，设置n_clusters=K、random_state=42、n_init=20。为什么n_init=20？K-means 初始聚类中心是随机选的，不同的初始点会收敛到不同的局部最优，多跑几次取最好的一次能显著提高稳定性。第一次跑时我建议直接设n_init=50观察一下惯性数值是否稳定，若再增加也不变化，说明收敛良好。
聚类完成后，聚类中心本身需要反标准化回物理坐标，因为后续优化模型使用的是真实的风速和辐照度数值。权重计算用np.bincount(kmeans.labels_) / len(labels)得到每个簇的概率权重。比如 K=20 时，权重可能是 0.08、0.12、0.05 等，所有权重之和必须等于 1。
最后把中心矩阵拆分回风速列和辐照度列，按场景组织保存。一个直接可用的输出格式是 CSV：每行一个场景，第一列为概率权重，后面 48 列分别为 24 小时风速和 24 小时辐照度。
from sklearn.cluster import KMeans from sklearn.preprocessing import StandardScaler # 场景矩阵 X: (10000, 48)，前 24 列风速，后 24 列辐照度 scaler = StandardScaler() X_scaled = scaler.fit_transform(X) kmeans = KMeans(n_clusters=20, random_state=42, n_init=20).fit(X_scaled) centers = scaler.inverse_transform(kmeans.cluster_centers_) # 反标准化 weights = np.bincount(kmeans.labels_) / len(kmeans.labels_) np.save("scenario_centers.npy", centers) np.save("scenario_weights.npy", weights)4.3 如何选 K 值和评价削减质量
K 值没有绝对标准，取决于下游模型承受能力和精度要求。工程上一般从 K=5 试到 K=30，逐步观察两个指标：削减前后统计量偏差和模型目标函数稳定度。
我自己的选型流程是：先算削减前场景全集的均值、标准差、90% 分位数、风力-辐照度相关系数这四个统计量；然后对每个 K 做一次削减，重算后四者的绝对误差。误差收敛的拐点就是比较合适的 K。实际操作中 K=15 到 K=25 之间往往已经能把均值误差压到 1% 以内，方差和分位数误差在 5% 以内，再往上增加 K 的效果有限，对优化模型反而加大计算负担。
还要检查削减后的相关系数。因为 K-means 是按欧氏距离聚类的，聚类中心天然会让变量之间相关性略微弱化，极端情况下可能从 -0.35 变成 -0.2。如果相关性损失太大，说明 K 太小或标准化方式不合适。一个补丁手段是把风速和辐照度各自的方差做等权缩放，而不是用全数据方差去缩放，这样聚类时两个变量对距离贡献更均匀。
下边是一个我以前项目里实际得到的削减评估表，K 值取 20 时的结果：
| 统计量 | 削减前 10000 场景 | 削减后 20 场景 | 相对误差 |
|---|---|---|---|
| 风速均值 (m/s) | 6.84 | 6.88 | 0.6% |
| 辐照度均值 (W/m²) | 176.3 | 174.5 | 1.0% |
| 风速标准差 (m/s) | 3.12 | 3.03 | 2.9% |
| 风速 90% 分位数 (m/s) | 10.6 | 10.3 | 2.8% |
| 风速-辐照度相关系数 | -0.31 | -0.27 | 12.9% |
前四行都还行，但相关系数偏差 12.9% 明显偏大。我当时的处理是把风速和辐照度分别做 PCA 白化后聚类，即对两个变量先做无关性变换再缩放，聚类中心再映射回来，相关系数偏差降到了 4% 以内。不过这个方法对可解释性略有损伤，适合分析场景内部关系质量要求高、不太在意中心可读性的场合。
5. 验证与分析：场景到底靠不靠谱
5.1 生成场景的统计特性校验
场景生成完了，不是看一眼散点图像不像就完事，要做三项校验。第一是边际分布还原是否准确：取生成场景的风速边缘分布和原始历史数据做 QQ 图，如果点没有偏离 45 度线，说明逆变换正确。第二是相关性结构是否保留：计算生成样本对的 Spearman 秩相关，应该和原始数据相差不大。第三是极值覆盖是否合理：看生成场景里是否存在同时出现“强风+强辐照”的组合，以及它的概率和历史经验是否一致。
我遇到过最典型的校验失败案例是：生成场景的边际分布完全正确、相关性也对，但是极端组合概率严重偏高。排查后发现是 Copula 拟合时用了 Pearson 相关而不是 Spearman 秩相关，导致在正态映射后尾部关联被夸大。换成秩相关后问题立即消失。这说明在 Gaussian Copula 里，参数 ρ 的估计方式看着只差一个公式，实际影响非常大。
5.2 削减前后误差评估
削减效果的评估不能只看统计量，还要看它进入优化模型后的表现。我的做法是把削减后的 20 个场景放进一个微网储能规划模型，求解得到经济最优配置；再把这 20 个场景全量随机抽样成 500 个新场景，用同样的模型去回测，看目标函数差异。如果两者差距在 3% 以内，说明削减场景已经足够代表不确定性。
比如某次微网项目里，20 场景给出的储能容量配置是 2.4 MWh，500 场景回测后的最优解是 2.5 MWh，差距 4%。略微超过 3% 的阈值，我尝试把 K 提升到 25，配置变成 2.46 MWh，差距降到 1.6%，可接受。所以我的经验是：如果下游是规划模型，K 值宁多勿少，因为规划决策对尾部风险很敏感。
5.3 边界条件与工程建议
这套方法不是万能的。Copula 适合捕捉稳定的相关结构，但如果研究区域有明显的季节切换，比如冬夏两季的风光相关性方向相反，一定要分季节建 Copula，否则混合数据拟合出来的相关参数两头不靠。我的做法是分春夏秋冬四个子数据集分别建模，生成时再按月拼合成全年场景集。
另一个边界是数据量。Copula 拟合对样本量要求不高，几百个小时的数据就能估出参数，但边际分布的尾部需要足够数据支撑。如果只有一年的数据，极端 1% 分位数基本是拍脑袋，建议在生成场景时给边际分布的尾部参数加涨落区间，做敏感性分析，别把单次拟合结果当真理。
6. 常见问题与实操避坑记录
6.1 边际拟合和相关性前置的“先有鸡还是先有蛋”
很多初学者会先估 Copula 相关参数，再回头拟合边际，顺序反了直接导致 Sklar 定理不成立。严格顺序一定是先拟合边际，再用边际 CDF 转换数据到均匀空间，最后在均匀空间里拟合 Copula。任何一步反过来，后面的逆变换都会出错。
6.2 零值和边界值导致的 ppf 溢出
逆变换时概率值取到 0 或 1，stats.weibull_min.ppf(0)会返回 0，ppf(1)会返回无穷大。在采样高斯分布随机数时理论上不会出现极端均匀值，但浮点运算总会产生极少数 1e-16 这样的数，我用np.clip(u, 1e-6, 1-1e-6)做过截断，简单有效。
6.3 K-means 聚类中心数量上限
K-means 不是聚类数越多越好。当 K 超过 50 后，聚类中心的代表性提升很慢，但下游优化模型的计算负担线性增长。而且 K-means 的代价函数随 K 单调下降，靠肘点法判断时经常找不到明显拐点。我建议结合模型求解时间来定 K，而不是纯看统计指标。
6.4 初始聚类中心的随机性
K-means 对初始中心敏感，尤其在场景数量大、特征维度高时，不同随机种子可能给出完全不同的聚类结果。我的经验是固定一个随机种子作为基线，做对比实验时所有方案都用同一随机种子；如果发现两次运行结果差异过大，就增大n_init或者改用 K-means++ 初始化，能大幅降低对初始值的敏感度。
6.5 风力和辐照度单位不同导致的聚类失衡
这个问题值得再强调一次。风速只有个位数，辐照度是几百，直接聚类时距离计算几乎被辐照度主导。StandardScaler是基础，但还不够，因为标准化后风速方差可能远小于辐照度方差，仍需考虑是否对两变量做额外权重调整。我的技巧是先分别标准化，再把辐照度乘以 0.7 的权重系数，手动调整变量对聚类距离的贡献，效果往往比单纯缩放更好。
这套方法在实际项目里我已经用了三轮迭代，从最初的手工采样脚本到现在的完整 pipeline，效率提升非常明显。最后分享一个小技巧：场景生成后一定保留一份原始未削减的 10000 场景文件，不要只存削减后的 20 个。因为验证阶段你经常需要对某个统计量做重新回测，有原始场景文件就不用重新采样，能省不少时间。
