简介一份面向风电领域科研与工程人员的论文复现资源围绕风电功率曲线异常数据清洗与风能资源精细化评估提供完整 Python 代码及逐段中文解释。资源覆盖 k-means、DBSCAN、Thompson tau、Copula 理论等主流清洗方法并实现贝叶斯变点检测与四分位法组合清洗策略用于识别和剔除风电机组实测功率中的异常记录。在此基础上进一步构建考虑风速风向联合分布的风能评估模型采用混合 Weibull 分布与 von Mises 分布描述风况特征并通过遗传算法优化 Logistic 函数拟合高精度功率曲线帮助读者形成数据清洗—风能评估—功率预测的完整技术链路。压缩包为 1 个 PDF 文件总大小 863KB内容包含原理讲解、可运行代码及运行结果分析适合具备一定数据分析基础的研究生、科研人员和风电场工程师。目前已有 69 人学习下载可用于实测数据质量控制、风电场选址及机组性能评估等场景。1. 风电功率曲线异常数据清洗为什么散点图上的离群点不能一刀切风电场SCADA数据拿回来第一件事不是算发电量而是看风速-功率散点图。理想情况下的散点图是一条厚度均匀的带子但实际数据里总散落着各种不听话的点低风速段功率飙到四五百、额定风速附近功率掉到几十、甚至整段数据在某个区间内变成另一条平行线。这些异常数据不清理功率特性曲线测试的准确性就会被打折扣IEC标准里那些bin均值也会被拖偏。《风电功率曲线异常数据清洗及考虑风速风向的风能资源评估》这篇论文核心就是把这堆异常数据按成因拆开、分门别类地清洗再做考虑风向的风能资源评估。这篇复现把论文里提到的k-means、DBSCAN、Thompson tau、Copula四种方法以及一个贝叶斯变点-四分位组合算法的完整Python代码都补齐了拿来就能跑。适合正在做功率曲线测试、机组性能评估或者写风电方向论文的工程师和研究生尤其是那种读过论文但作者没开源代码、只能自己从头抠的场景。2. 四种经典清洗方法的选型与实现从k-means到Copula2.1 k-means先标准化再谈聚类k-means是异常清洗里最直白的一个思路把风速-功率数据聚成两类样本量大的那一类当作正常数据小簇当作异常。这个假设的前提是异常数据在散点图上会聚成独立的云团比如传感器持续漂移、通讯故障导致的功率值整体偏移这类异常确实会形成一个小簇。from sklearn.cluster import KMeans from sklearn.preprocessing import StandardScaler def kmeans_clean(data, n_clusters2, random_state42): # 风速和功率量纲差异太大风速10左右功率可能上千 # 不标准化的话欧氏距离会被功率完全主导聚类结果约等于按功率切一刀 scaler StandardScaler() scaled scaler.fit_transform(data[[wind_speed, power]]) kmeans KMeans(n_clustersn_clusters, random_staterandom_state) data data.copy() data[cluster] kmeans.fit_predict(scaled) # 假设正常数据集中在样本量最大的簇中 counts data[cluster].value_counts() normal_cluster counts.idxmax() cleaned data[data[cluster] normal_cluster].drop(cluster, axis1) return cleaned kmeans_cleaned kmeans_clean(data)这里做了两步关键处理。第一步是标准化直接用StandardScaler把风速和功率都压到均值0、方差1的空间里这是k-means这类基于欧氏距离的算法最基本的前置条件。第二步是保留最大簇背后的逻辑是正常数据在功率曲线附近密度最高异常不管聚成几个小簇样本总量通常远小于正常数据。如果实际数据里正常簇本身就有多个形态比如不同控制策略下的功率曲线差异明显可以把n_clusters从2调大到3或4然后按簇内平均功率曲线是否符合风机出厂功率曲线来选正常簇而不是只看样本数。2.2 DBSCAN密度聚类对散点型异常更友好k-means的短板是必须提前告诉它聚几类而且对簇形状敏感。DBSCAN不需要指定簇数它根据密度把数据分成若干簇密度不够的点直接标成噪声。功率曲线散点图上正常数据沿曲线密集分布异常点则以低密度形式散落在曲线外的空白区这种结构天生适合DBSCAN。from sklearn.cluster import DBSCAN from sklearn.preprocessing import StandardScaler def dbscan_clean(data, eps0.4, min_samples8): # 同样先标准化注意 eps 是在标准化空间里的距离阈值 scaler StandardScaler() scaled scaler.fit_transform(data[[wind_speed, power]]) db DBSCAN(epseps, min_samplesmin_samples) data data.copy() data[cluster] db.fit_predict(scaled) # DBSCAN 中标签为 -1 的点就是噪声点直接视为异常 cleaned data[data[cluster] ! -1].drop(cluster, axis1) return cleaned dbscan_cleaned dbscan_clean(data)eps是核心参数表示两个点被认定为一类的最远距离。标准化之后的数据空间里eps取0.3到0.8是比较常见的范围我习惯先取0.4试跑再看清洗率。min_samples表示一个核心点周围至少要有多少个点才算密集区取值5到10较为合理数据量大的话可以适当加大。这里最需要注意的是如果把原始坐标系下的eps值比如0.5直接套到标准化后的数据上结果往往会全变成噪声或者一个点都没清洗掉这个坑后面专门讲。2.3 Thompson tau单变量统计的快速筛查Thompson tau法的思路和前面两种完全不同它不关心风速只盯着功率这个单变量用均值加减若干倍标准差圈出一个正常区间落在区间外的就是异常。这个方法的好处是快、可解释性强适合在对数据还没什么先验认识的时候先跑一遍粗筛。def thompson_tau_clean(data, columnpower, tau1.5): values data[column] mean np.mean(values) std np.std(values) # 功率偏离均值超过 tau 倍标准差即视为异常 outliers np.abs(values - mean) tau * std cleaned data[~outliers] return cleaned thompson_cleaned thompson_tau_clean(data)原文里tau取的是1.0这个值偏激进会误伤很多正常数据。实际使用我一般从1.5起步正常数据量大、异常占比低的时候可以放到2.0。这个方法有个天然缺陷功率的方差随风速增大而增大靠一个全局的标准差来圈正常区间低风速段容易被误杀高风速段又可能漏掉真正的异常。所以我通常只把它当作预清洗手段把明显离谱的点先剔掉而不是当作最终方案。更稳的做法是先把数据按风速分箱在每个箱内单独做Thompson tau这样至少避开了方差随风速变化的问题。2.4 Copula把联合分布写进清洗逻辑Copula理论是这里建模层面最扎实的一种方法。风速和功率的关系不是简单的线性相关用二维高斯去拟合散点图会很勉强Copula的思路是把单变量的边缘分布和变量之间的依赖结构拆开分别建模然后重新组合成联合分布。每个数据点都能算出它在联合分布下的概率密度密度越低越可能是异常。from copulas.multivariate import GaussianCopula def copula_clean(data, threshold_percentile5): copula GaussianCopula() copula.fit(data[[wind_speed, power]]) # 计算每个点的联合概率密度低密度点视为异常 probs copula.pdf(data[[wind_speed, power]]) threshold np.percentile(probs, threshold_percentile) cleaned data[probs threshold] return cleaned copula_cleaned copula_clean(data)这个方法的门槛在于copulas库的安装和API版本兼容不同版本里GaussianCopula这个类名可能不一样有的版本叫GaussianMultivariate。另外pdf函数返回的是密度值而不是概率值密度值可以大于1所以不能用固定的绝对阈值而是用分位数来定默认取5%分位数意味着清洗掉联合概率密度最低的5%的点。这个比例需要根据实际异常率调整异常明显偏多的时候可以提高到10%反之降到2%。Copula方法对样本量也有要求样本少于500个点时拟合出来的联合分布不稳定PDF经常出现数值异常这一点在后面的避坑章节会展开。四种方法放在一起看选型逻辑其实很清晰方法异常形态假设关键参数适用场景主要短板k-means异常聚成少量小簇n_clusters传感器漂移、系统性偏移簇数难定量纲敏感DBSCAN异常是低密度孤立点eps, min_samples散点型离群数据eps需要调密度敏感Thompson tau异常分布在单变量尾部tau快速粗筛、单维筛查无法利用风速信息Copula低联合概率密度分位数阈值非线性关联、复杂分布模型是黑匣子样本要求高3. 贝叶斯变点-四分位组合算法处理工况切换型异常的实战方案3.1 为什么需要组合算法工况切换是聚类方法的盲区前面四种方法处理的是散点型异常和偏移型异常但风电数据里还有一类更难缠的功率曲线发生整体位移。风机限电、变桨策略切换、部分机组停机再启动都会让风速-功率关系在一段时间内变成另一条平行的曲线。这条另一条曲线上的点本身完全正常不落在任何低密度区域k-means和DBSCAN拿它没办法Copula也只能把它当成另一个正常模式。贝叶斯变点-四分位组合算法的核心思路就是先找到功率时间序列里的均值突变位置再把突变点附近的一段数据连同四分位法检出的离群点一起删掉。3.2 贝叶斯变点-四分位组合算法实现import numpy as np import pandas as pd class BayesianChangePointIQR: def __init__(self, window_size50, threshold0.95, iqr_multiplier1.5): self.window_size window_size # 变点检测的滑动窗口大小 self.threshold threshold # 贝叶斯因子的判定阈值 self.iqr_multiplier iqr_multiplier # 四分位法的倍数 def calculate_bayes_factor(self, data1, data2, mean1, mean2, std1, std2): 简化版贝叶斯因子结合效应量和样本量的大小 n1, n2 len(data1), len(data2) pooled_std np.sqrt(((n1 - 1) * std1 ** 2 (n2 - 1) * std2 ** 2) / (n1 n2 - 2)) if pooled_std 0: return 0 # 效应量均值差除以合并标准差乘以样本量调整因子 effect_size abs(mean1 - mean2) / pooled_std return effect_size * np.sqrt(n1 * n2 / (n1 n2)) def bayesian_change_point_detection(self, data): 在功率时间序列上滑动检测变点 n len(data) change_points [] for i in range(self.window_size, n - self.window_size): before data[i - self.window_size:i] # 变点前窗口 after data[i:i self.window_size] # 变点后窗口 bf self.calculate_bayes_factor( before, after, np.mean(before), np.mean(after), np.std(before), np.std(after) ) if bf self.threshold: change_points.append(i) return change_points def iqr_outlier_detection(self, data): 四分位法基于 IQR 圈定正常区间 q1, q3 np.percentile(data, [25, 75]) iqr q3 - q1 lower q1 - self.iqr_multiplier * iqr upper q3 self.iqr_multiplier * iqr return (data lower) | (data upper) def combined_cleaning(self, wind_speed, power): 组合清洗四分位法负责全局离群变点检测负责区间位移 change_points self.bayesian_change_point_detection(power) outliers self.iqr_outlier_detection(power).copy() # 把每个变点前后各 10 个点全部标记为异常 # 因为功率切换是渐变过程不能只删变点本身那一个点 for cp in change_points: start max(0, cp - 10) end min(len(power), cp 10) outliers[start:end] True clean_mask ~outliers cleaned pd.DataFrame({ wind_speed: wind_speed[clean_mask], power: power[clean_mask] }) return cleaned, outliers # 使用示例 bcp_iqr BayesianChangePointIQR() cleaned_data, outlier_flags bcp_iqr.combined_cleaning(wind_speed, power)变点检测的核心在calculate_bayes_factor。这个函数把前后两个窗口的均值差、合并标准差和样本量组合成一个标量值越大说明这两个窗口越可能来自不同的分布。严格来说这是贝叶斯因子的工程近似不是完整的边际似然比但作为变点的判据足够有效。窗口大小window_size50意味着拿前后各50个点做比较对应10分钟级采样就是约10分钟的数据窗口。如果采样频率更高比如1秒一个点窗口要相应调大到几百。3.3 参数调整与适用边界这个算法的三个参数各有各的调整逻辑。threshold控制变点检测的灵敏度取0.95意味着只有前后窗口差异足够大才认定为变点调低到0.8会检出更多变点但也更容易把正常波动误判成切换。iqr_multiplier和常用箱线图里的1.5倍IQR一致想让清洗更保守就调到2.0或2.5。变点前后各10个点的扩展范围也不是固定的如果是缓慢的功率爬坡过程比如限电指令下发后功率花了半小时才降下来这个范围要扩到30甚至更大。这个算法最适合的数据形态是带时间顺序的SCADA时序数据因为变点检测依赖时间序列的前后关系。如果手里的数据已经被随机打乱或者只有风速-功率两列而没有时间戳这个算法的变点检测部分就发挥不了作用只能退化成纯四分位清洗。另外要注意算法里的四分位法是直接作用在整个功率序列上的没有按风速分箱所以遇到风速分布极不均匀的数据还是要配合按风速段的预清洗使用。4. 风速-风向联合风能评估Weibull参数估计与扇区功率密度计算4.1 Weibull参数估计图解法与极大似然法风能评估的第一步是把风速的概率分布拟合出来。双参数Weibull分布是风电行业的事实标准形状参数k决定分布形态尺度参数c决定平均风速大小。论文里给了两种估计方法图解法适合快速估计和提供初值极大似然法更精确。import numpy as np from scipy import stats from scipy.optimize import minimize def graphical_estimation(wind_speed): 图解法对经验CDF做双对数变换用线性回归拟合 Weibull 参数 sorted_v np.sort(wind_speed) n len(sorted_v) # 经验累积分布函数用 (i)/(n1) 避免端点出现 ln(0) emp_cdf np.arange(1, n 1) / (n 1) # Weibull 分布的 CDF 经过变换后是线性关系 # ln(-ln(1-F(v))) k * ln(v) - k * ln(c) x np.log(sorted_v) y np.log(-np.log(1 - emp_cdf)) slope, intercept, _, _, _ stats.linregress(x, y) k slope c np.exp(-intercept / k) return k, c def weibull_mle(wind_speed): 极大似然估计最小化负对数似然 def neg_log_likelihood(params): k, c params if k 0 or c 0: return np.inf # 形状和尺度参数必须为正 return -np.sum(stats.weibull_min.logpdf(wind_speed, k, scalec)) # 初始值形状参数从2开始尺度参数用平均风速 init [2.0, np.mean(wind_speed)] result minimize(neg_log_likelihood, init, bounds[(0.1, 10), (0.1, 30)], methodL-BFGS-B) return result.x[0], result.x[1]图解法的数学基础是Weibull分布的CDF变换ln(-ln(1-F(v))) k ln(v) - k ln(c)所以对经验CDF做双对数变换后用线性回归回归斜率就是k截距换算得到c。这个方法的优点是计算简单、不需要迭代缺点是经验CDF两端的数据点波动大会影响回归精度。极大似然法通过scipy.optimize.minimize对负对数似然最小化这里用了stats.weibull_min.logpdf其中k是形状参数scalec是尺度参数和手写Weibull PDF是等价的。实际使用中我的习惯是先跑图解法拿到初值再喂给MLE做精细迭代这样能显著降低MLE陷入局部最优的概率。4.2 风向分箱与联合分布建模风速分布拟合好之后风向的处理是这篇论文复现里比较有工程味道的部分。风向是圆形变量0度和360度本质上是同一个方向不能拿来做简单的线性统计。论文摘要提到用von Mises分布建模风向但在工程复现里更常用也更容易验证的做法是等宽分箱加逐扇区拟合。def joint_wind_distribution(wind_speed, wind_direction, n_sectors12): 把风向分成 n_sectors 个扇区每个扇区独立拟合 Weibull 分布 sector_width 360 / n_sectors result [] for i in range(n_sectors): lower i * sector_width upper (i 1) * sector_width mask (wind_direction lower) (wind_direction upper) sector_speeds wind_speed[mask] # 样本量太少的扇区跳过避免 Weibull 拟合过拟合 if len(sector_speeds) 10: k, c weibull_mle(sector_speeds) result.append({ sector: i, frequency: len(sector_speeds) / len(wind_speed), k: k, c: c, mean_speed: np.mean(sector_speeds) }) return result扇区宽度选30度12扇区是权衡的结果扇区越多风向分辨率越高但每个扇区的样本量会变少。如果一个扇区的样本少于10个MLE拟合出来的Weibull参数波动会很大k值可能从1.5跳到3.5这种情况下把相邻扇区合并是更稳的做法。frequency这一项在后面的能量评估里会作为权重因为每个扇区的出现频率不同对总风能资源的贡献也不同。4.3 平均功率密度与风能资源量化风能资源评估的最终输出是这个场址每平方米扫风面积上平均能拿到多少功率。功率密度和风速成立方关系所以风速分布尾部哪怕只差一点点能量评估结果都会差很多。def wind_power_density(v, air_density1.225): 单位面积风功率密度W/m2 return 0.5 * air_density * v ** 3 def assess_wind_energy(joint_params, rotor_diameter80): 基于各扇区 Weibull 分布数值积分计算全场平均风功率 air_density 1.225 rotor_area np.pi * (rotor_diameter / 2) ** 2 # 风速积分范围取 0 到 25 m/s覆盖绝大多数风况 v_range np.linspace(0, 25, 1000) total_power 0 sector_results [] for p in joint_params: # 该扇区的风速概率密度 pdf_values stats.weibull_min.pdf(v_range, p[k], scalep[c]) # 该扇区的期望功率密度对 pdf * 功率密度做数值积分 expected_density np.trapz(pdf_values * wind_power_density(v_range), v_range) # 乘上扫风面积得到期望功率再乘扇区频率 sector_power expected_density * rotor_area * p[frequency] total_power sector_power sector_results.append({ sector: p[sector], frequency: p[frequency], expected_power_kw: sector_power / 1000, k: p[k], c: p[c] }) return total_power / 1000, sector_resultsnp.trapz是梯形数值积分1000个积分点在0到25m/s的范围内足够平滑。这里的计算逻辑是单个扇区的Weibull PDF描述了这个方向上风速的概率分布把PDF和功率密度函数相乘再积分得到这个扇区的期望功率密度乘上扫风面积和扇区出现频率后就是该扇区对全场平均功率的贡献。把12个扇区累加起来就得到整个场址的平均功率。需要注意如果传入的数据只覆盖机组运行时段评估结果会偏高因为停机时段的零功率没有计入平均。更严谨的做法是先统计机组可利用率把评估结果乘以可利用率系数。5. 实操避坑五个最常翻车的参数与数据陷阱5.1 参数与量纲陷阱坑一不标准化直接跑k-means或DBSCAN。现象是清洗出来的结果看起来完全合理但仔细对比会发现被保留的数据几乎完全按功率大小排列风速信息被忽略了。原因是风速的量纲是10左右功率的量纲是几百甚至上千欧氏距离的计算被功率完全主导聚类中心的位置基本由功率决定。解决方法是聚类前强制做标准化我用StandardScaler已经写进前面的代码里了自己实现的时候最容易漏这一步。坑二DBSCAN的eps在原始坐标系里瞎调。现象是eps取0.5聚类结果要么全部是噪声点要么全部归为一类两个极端之间几乎没有过渡。原因是在原始坐标系下风速和功率的尺度差一个数量级0.5这个距离阈值在功率维度上微不足道但在风速维度上又太大。解决方法是先把数据标准化到同一尺度再去调eps。我一般用k-距离图辅助选值把每个点到第min_samples近邻的距离排序画出来找曲线拐点对应的距离那就是eps的合理起点。坑三Thompson tau的tau值拍脑袋乱取。现象是清洗后低风速段的正常点被成片误删。原因是功率方差随风速增大而增大全局一个标准差低风速段的正常波动相对于全局std显得异常高风速段的真实异常反而被淹没。解决方法是按风速分箱在每个箱内单独计算均值和标准差或者用中位数和MAD绝对中位差代替均值和标准差MAD对异常值的鲁棒性要好得多。5.2 环境与验证陷阱坑四copulas库的API版本和数值稳定性问题。现象是from copulas.multivariate import GaussianCopula报导入错误或者pdf函数返回一堆NaN。原因是不同版本的copulas库类名不一致新版本改成了GaussianMultivariate另外数据里有缺失值或inf值时拟合过程会产生NaN。解决方法是先检查数据有没有NaN和inf用np.isfinite过滤一遍如果导入失败查一下已安装版本改用GaussianMultivariate样本量小于500时直接用scipy.stats.multivariate_normal代替效果接近且少一个依赖。坑五风向扇区分得太细导致Weibull拟合翻车。现象是某个扇区的k值突然跳到4以上或者c值明显偏离相邻扇区画出PDF后发现曲线形状完全不合理。原因是样本量不足MLE迭代到边界或者过拟合了少量样本。解决方法是扇区宽度至少保持30度并且保证每个扇区样本量不少于20个样本不够就把相邻扇区合并。我在写代码时把样本量阈值设在了10实际使用建议提高到20宁可牺牲一点风向分辨率也要保证参数估计的稳定性。坑六清洗完没有量化验证纯靠肉眼调参。现象是同一份数据两个人调出来的清洗结果差异巨大一个删了3%另一个删了15%但散点图看起来都挺干净。原因是清洗参数没有统一标准个人主观判断占了主导。解决方法是固定三个输出指标清洗率删除点数占比、异常样本的功率分布范围、清洗前后风速分箱均值变化每次都记录这三个值调参时对比指标变化而不是对比图片。6. 进阶用法用遗传算法拟合Logistic功率曲线并做清洗质量验证6.1 遗传算法拟合Logistic功率曲线清洗完数据之后下一步通常是要出一条平滑的标准功率曲线。IEC标准里的bin平均法对bin宽度敏感bin取太宽曲线太平取太窄又有很多空bin。论文摘要里提到用遗传算法优化Logistic函数来做拟合工程上我常用scipy.optimize.differential_evolution来替代手写遗传算法效果接近且不需要额外装DEAP库。from scipy.optimize import differential_evolution def logistic_power(v, a, b, c, d): 四参数 Logistic 功率曲线 d 是切入功率下限a 是额定功率附近的上限 c 是曲线中点拐点风速b 控制曲线陡峭程度 return d a / (1 np.exp(-b * (v - c))) def fit_logistic_curve(wind_speed, power): def loss(params): a, b, c, d params pred logistic_power(wind_speed, a, b, c, d) return np.sqrt(np.mean((pred - power) ** 2)) # RMSE # 参数边界a 和 d 限制功率范围b 限制陡峭度c 限制拐点风速 bounds [(100, 2000), (0.1, 2.0), (5, 15), (-50, 50)] result differential_evolution(loss, bounds, seed42, tol1e-6) return result.x a, b, c, d fit_logistic_curve( cleaned_data[wind_speed].values, cleaned_data[power].values )选择遗传算法而不是梯度下降是有原因的。这个四参数Logistic的损失函数是非凸的scipy.optimize.curve_fit这类基于梯度的优化器对初值非常敏感初值选得不好很容易陷进局部最优。differential_evolution在参数边界内做全局搜索不需要给定好的初值只要边界范围合理基本能稳定收敛。边界设置上a的上限2000要覆盖机组的额定功率c的范围5到15要根据机型的切入和切出风速来定不同机型的差异会很大。6.2 用拟合残差反查清洗质量拟合本身不是目的用拟合结果反过来验证清洗质量才是这套流程里最有价值的一步。拟合完成后把每个风速段的残差实际功率减拟合功率按风速分组做统计然后看下面的体检表检查项判断标准说明整体R²大于0.95低于0.9说明拟合形式或数据质量有问题分段残差均值每个风速段都接近0某段系统性偏正/偏负说明那里还藏着未清洗的异常分段残差标准差随风速增大而增大但平滑出现突跳说明该风速段有局部异常簇高风速段残差不能系统性为负负值代表清洗过度把正常满发数据误删了如果某个风速段的残差均值偏离0超过额定功率的3%我就会回到清洗环节把该风速段的数据单独拉出来看。常见的情况是某段风速区间内存在一批功率整体偏低的点它们看起来像一条平行线四分位法检不出来但残差分析一跑就露馅。这个时候针对这个风速段单独做一次DBSCAN或Copula清洗比全局重新调参有效得多。从那以后我每次做功率曲线清洗都会强制走一遍清洗拟合残差诊断的闭环不再只看散点图顺不顺眼。参数改了、样本换了残差分布能直接告诉你这套清洗参数站不站得住脚省掉了大量靠肉眼反复对比图表的重复劳动。希望这个流程对你有帮助。本文还有配套的精品资源点击获取