简介面向风电领域科研人员、工程技术人员及高校研究生的论文复现资料包主题为风电功率曲线异常数据清洗及考虑风速风向的风能资源评估。内容系统分析异常数据成因及其对功率特性建模的影响实现k-means、DBSCAN、Thompson tau法与Copula理论等多种清洗方法并给出贝叶斯变点检测与四分位法相结合的组合清洗算法。在此基础上构建考虑风速风向联合分布的风能资源评估模型采用混合Weibull分布与von Mises分布建模风速风向通过遗传算法优化Logistic函数实现高精度功率曲线拟合形成从数据清洗、风能评估到功率预测的完整技术链条。压缩包共1个PDF文件大小863KB内含完整可运行的Python代码及详细中文解释已有69人学习下载。读者可重点研读组合清洗策略、风速风向联合建模与遗传算法曲线拟合的实现细节通过仿真与实测数据验证各项技术有效性。1. 风电功率曲线异常数据清洗先过数据关再谈风资源评估拿到一座风电场一整年的SCADA数据直接在坐标轴上画风速-功率散点图你大概率看到的不是一条教科书式的平滑曲线而是一朵“蒲公英”风速10m/s时功率可能只有几十千瓦风速4m/s时功率却冲到额定值。这些异常点如果不做清洗就进入评估后面算风速频率、风能密度都会严重失真反过来说如果清洗过度把高风速段的样本误删评估结果又会系统性偏低。这篇按论文复现的通用思路把风电功率曲线异常数据清洗和考虑风速风向的风能资源评估串成一条可执行的链路其中清洗部分代码可以直接抄走改数据。2. 异常数据从哪来五类噪声的物理来源与清洗选型思路2.1 功率曲线的四个区间与五类典型噪声先讲功率曲线的分段。常见风机的切入风速在3m/s左右额定风速在11m/s左右切出风速在25m/s左右功率曲线大致分成四段切入前区间风速低于3m/s理论上功率为0机组处于待机或自耗电状态切入-额定区间3~11m/s功率随风速近似按三次方关系爬升这是风资源转化的主区间额定区间11~25m/s功率被控制在额定值附近叶片变桨距限功率切出区间接近或超过25m/s风机逐步收桨停机功率快速回0理解了区间再看异常点的形态就清楚多了。我平时看散点图习惯先找五类典型噪声负功率点。出现在风速低、机组待机或自耗电时段SCADA把“从电网取电”记成了负值。超额定功率点。功率冲到额定值的120%甚至更高多发生在切入瞬间的功率尖峰属于控制瞬态或传感器采样失真。切入前高功率点。风速只有3~4m/s功率却跑到几百千瓦最常见原因是风向标结冰导致风速计读数偏低或者上一时段数据被错误对齐到当前时间戳。限功率点。风速在额定区间功率却被调度指令压低形成一条平行于风速轴的水平带。这类点在物理上真实存在不是传感器错误论文复现时怎么处理它直接决定高风速段的样本量。卡死数据点。风速长时间恒定在一个值附近功率却在波动对应风速计轴承磨损或结冰后的假读数。这五类噪声里第4类最容易被没经验的人误删。一条原则先区分“脏数据”和“真实但偏离理想曲线”的数据前者必须删后者要视评估目标决定是否保留。做风资源评估时我倾向于不保留限功率点做功率曲线标定时则要看论文的具体定义。提示限功率点是真实物理数据它的去留要按评估目标决定不要一刀切。2.2 清洗基线怎么选分箱IQR为主DBSCAN做补充论文复现里最常用的清洗方法有两类基于统计的分箱四分位距bin-IQR和基于密度的DBSCAN聚类。我的建议是用bin-IQR做主力用DBSCAN处理稀疏边界箱不要一上来就上机器学习。bin-IQR的逻辑是把风速轴切成若干个等宽窗口在每个窗口内单独看功率的分布再用四分位距识别离群点。它最大的优点是可解释性强每个点为什么被删、删了多少都能对着散点图讲清楚。复现论文时审稿人和读者最看重的就是这种可追溯的清洗逻辑。DBSCAN的优点是它不需要预先指定分箱宽度能自动发现任意形状的密集区域。但它的两个参数eps和min_samples对结果非常敏感而且对量纲极其敏感——风速是0~30m/s功率是0~2000kW两个维度直接算欧氏距离的话功率维度会完全主导聚类结果结果基本是玄学。所以我在实际复现时采用“物理筛查 → bin-IQR → 稀疏箱DBSCAN补充”的组合。另外注意别用K-means做这个任务。K-means假设簇是凸的、各向同性的而功率曲线的正常点是一条从原点向右上延伸的长条形带用K-means会把长条带中段切碎清洗出的结果很难看。表bin-IQR与DBSCAN的选型对比对比维度分箱IQRDBSCAN参数数量箱宽、k值eps、min_samples可解释性高逐箱可追溯中需结合散点图判断稀疏数据处理容易误删相对稳健计算开销极小随样本量增大明显变慢2.3 动手之前画一张不糊的功率散点图数据没进清洗流程前先画一张全量散点图这一步能帮你后面判断清洗参数是否合理。直接plt.scatter画出来几十万点叠在一起就是一团黑看不到任何结构。我一般把点调小、把透明度压低import pandas as pd import matplotlib.pyplot as plt df pd.read_csv(scada_10min.csv, parse_dates[time]) df df[[time, wind_speed, wind_direction, active_power]].dropna() fig, ax plt.subplots(figsize(12, 6)) ax.scatter(df[wind_speed], df[active_power], s1, alpha0.3, colorsteelblue) ax.set_xlabel(风速 (m/s)) ax.set_ylabel(有功功率 (kW)) ax.set_title(原始功率散点图先找结构再定清洗参数) ax.grid(alpha0.3) plt.tight_layout() plt.show()s1是把点的面积压到极小alpha0.3让重叠区域自动变深色。这样低密度区域的异常点会以稀疏浅色出现正常点聚成深色带。看这张图时重点确认三件事有没有明显的上下水平线异常带、切入风速附近有没有高功率集团、高风速段是不是已经很稀疏。散点图本身不参与清洗它是用来帮你定参数的。如果图上高风速段本来就稀疏那么第3章的稀疏箱处理策略就必须格外小心这一步比清洗算法本身更影响结果。3. 一个能跑的清洗流程物理筛查、分箱IQR与稀疏箱边界处理3.1 第一关物理规则筛查清洗的第一步永远不是花哨算法而是把不可能存在的点先删掉。物理规则筛查的作用是拦截传感器级错误为后续统计清洗减轻压力。我一般用四个边界风速必须大于0且小于35m/s超过35m/s的风速计读数基本是仪器故障功率必须大于等于0且小于等于额定功率的1.05倍SCADA偶尔会记录到短暂超发风向在0~360度之间风速明显高于切入风速但功率为0的点大多是数据对齐错误SCADA数据本身也经常有字段缺失有时没有风向字段有时功率字段有跳变。先做dropna()去掉时间戳、风速、功率为空的记录风向字段如果后面要用必须在物理筛查阶段就一起清洗不要等到算扇区时才想起来了。def physical_filter(df, rated_power2000.0, cut_in3.0): 物理规则筛查按SCADA字段的物理边界剔除不可能数据 df df.copy() # 基础字段非空 df df.dropna(subset[wind_speed, active_power]) # 风速边界0~35 m/s df df[(df[wind_speed] 0.0) (df[wind_speed] 35.0)] # 功率边界0 ~ 额定功率*1.05 df df[(df[active_power] 0.0) (df[active_power] rated_power * 1.05)] # 风向边界0~360需要风向字段时用 if wind_direction in df.columns: df df[(df[wind_direction] 0.0) (df[wind_direction] 360.0)] # 可疑的“零功率高风速”点风速高于切入风速2且功率0多为对齐错误 df df[~((df[active_power] 0.0) (df[wind_speed] cut_in 2.0))] return df这个函数里rated_power、cut_in要按手里的风机型号改上面是2MW机型的默认值。最后一条零功率高风速的筛选用的是“风速大于切入风速2”为什么加2因为切入风速本身就有2~3m/s的过渡带风速4m/s时功率为0可能是正常的欠速状态留出余量是为了不过早误删。这一条是可选规则论文里如果没有相应清洗步骤可以去掉。3.2 第二关分箱四分位距物理筛查之后剩下的异常大多是限功率、湍流、偏航误差造成的偏离型数据。分箱四分位距的做法是把风速以固定宽度分箱在每个箱内统计功率的Q1和Q3将落在[Q1 - k*IQR, Q3 k*IQR]之外的样本标记为离群。通常k取1.5这是箱线图的标准值但功率曲线数据往往呈偏态分布一些论文会用2.0或更大的k来减少误删。实际复现时建议先按论文参数跑一遍再对比清洗前后散点图看低风速段有没有被过度清洗。def bin_iqr_filter(df, bin_width0.5, k1.5, min_samples10): 按风速分箱逐箱做四分位距离群剔除 df df.copy() # 分箱截断到bin_width的整数倍 df[bin] (df[wind_speed] // bin_width) * bin_width keep_list [] for bin_edge, group in df.groupby(bin): # 单箱样本量太少时统计量不可靠直接保留该箱 if len(group) min_samples: keep_list.append(group) continue q1 group[active_power].quantile(0.25) q3 group[active_power].quantile(0.75) iqr q3 - q1 # 物理下界不能小于0上界不设硬限交给IQR判断 lower max(q1 - k * iqr, 0.0) upper q3 k * iqr mask (group[active_power] lower) (group[active_power] upper) keep_list.append(group[mask]) cleaned pd.concat(keep_list).drop(columnsbin) return cleaned三个参数各有讲究。bin_width常见取0.5m/s因为SCADA风速精度通常在0.1m/s级别0.5m/s一个箱能保证每箱有足够样本又能分辨曲线局部的斜率变化如果整体数据量不足放宽到1.0m/s。min_samples是箱内样本阈值一台风机一年10分钟平均数据约5万条0.5m/s分箱后大多数箱有几百条样本个别高风速箱可能只有几十条设10条能避免统计噪声。k是离群容忍度k越大保留越多论文里从1.5到3都有做敏感性分析时一般就围绕k展开。提示论文给的k值不一定适合你的机组调参前先记录原始散点图的形态。3.3 第三关边界箱与稀疏箱怎么处理分箱四分位距在数据密集区表现很好但只要接近切出风速样本量就会断崖式下降。风速20m/s以上一年可能只有几十条记录这时IQR会被少量正常点拉宽离群点反而被保留或者反过来正常点被当成离群点误删。稀疏箱的常见处理策略有三种策略做法适用场景直接保留样本量小于min_samples时整箱不清理高风速数据少且无明显异常放宽k稀疏箱把k放大到3~4高风速段存在少量噪声但不想误删DBSCAN补充对稀疏箱内点做密度聚类只删标签为-1的点箱内样本形状不规则时最稳我一般用“直接保留 DBSCAN补充”的组合。直接保留保证不删错DBSCAN只针对那些聚成小团、粘在正常带边缘的点。因为稀疏箱里点本来就少DBSCAN反而比IQR更稳。from sklearn.cluster import DBSCAN from sklearn.preprocessing import RobustScaler def sparse_bin_dbscan(group_df, eps0.2, min_samples3): 稀疏箱专用清洗归一化后做密度聚类删除噪声点 features group_df[[wind_speed, active_power]].values # RobustScaler用中位数和四分位距归一化离群点不会拉偏尺度 scaler RobustScaler() scaled scaler.fit_transform(features) labels DBSCAN(epseps, min_samplesmin_samples).fit_predict(scaled) return group_df[labels ! -1]RobustScaler对离群点不敏感比StandardScaler更适合这种场景。eps0.2是在归一化空间里的邻域半径min_samples3是针对稀疏箱的特殊值。这段代码只接收单箱数据如果你的稀疏箱有几十个点min_samples可以调大到5。3.4 清洗效果怎么量化清洗不是“删了就完”。论文复现里你要向读者证明清洗确实有效三个指标够用保留样本比例、分箱平均功率曲线、清洗前后功率RMSE。def cleaning_summary(raw_df, cleaned_df, bin_width0.5): 输出清洗前后的样本量与分箱平均功率对比 raw_df raw_df.copy() cleaned_df cleaned_df.copy() for d in (raw_df, cleaned_df): d[bin] (d[wind_speed] // bin_width) * bin_width raw_curve raw_df.groupby(bin)[active_power].mean() clean_curve cleaned_df.groupby(bin)[active_power].mean() # RMSE按公共箱计算 common_bins raw_curve.index.intersection(clean_curve.index) rmse np.sqrt(((raw_curve[common_bins] - clean_curve[common_bins]) ** 2).mean()) print(f原始样本数: {len(raw_df)}, 清洗后样本数: {len(cleaned_df)}, f保留率: {len(cleaned_df) / len(raw_df):.1%}, RMSE: {rmse:.1f} kW) return rmse保留率一般落在75%~95%区间。如果保留率低于70%大概率是清洗参数过严回头看k和bin_width如果高于98%说明数据本身质量好IQR可能没删掉该删的点。4. 考虑风速风向的风能资源评估扇区统计、空气密度修正与结果对比4.1 风向数据的坑机舱风向标、死区与基准方向SCADA风向来自机舱顶部的风向标它测的不是气象学里的来流风向而是“机舱相对风”。风机在偏航过程中机舱轴线与来流之间总有一个偏航误差角小风时段能达到正负10~20度。所以拿到原始风向不要直接画风玫瑰先做两件事检查风向频率分布有没有明显的缺角。如果某个角度区间样本数几乎为零可能不是真实风况而是风向标死区或数据记录问题。检查风向是否长时间固定在某几个值。这种卡死数据如果不处理扇区统计会失真。把风向做成直方图看一眼最直接fig, ax plt.subplots(figsize(10, 5)) df[wind_direction].hist(bins72, axax, edgecolorwhite) ax.set_xlabel(风向 (deg)) ax.set_ylabel(样本数) ax.set_title(风向频率直方图检查死区与卡死) plt.show()如果直方图有角度断层再决定是删除对应区间的数据还是回去查SCADA日志。还有一个复现时特别容易踩的坑论文里用的是气象学风向0°正北、顺时针递增还是机舱相对风向0°机舱朝向。两者的零度基准不同直接照抄论文代码最容易在这里翻车。4.2 扇区划分与风玫瑰统计风向统计口径常见的是16扇区每个22.5°或12扇区每个30°。年数据量大用16扇区数据量小用12扇区更稳。不要用循环逐行数数用pandas的groupby就够def wind_rose_stats(df, sectors16): 按扇区统计风向频率、平均风速、风速立方均值 df df.copy() dir_step 360.0 / sectors # 把风向映射到扇区编号 df[sector] (df[wind_direction] // dir_step).astype(int) # 扇区中心角用于画图标注 df[sector_center] (df[sector] * dir_step dir_step / 2) % 360 stats [] total len(df) for sector, group in df.groupby(sector): n len(group) stats.append({ sector: sector, sector_center: group[sector_center].iloc[0], freq: n / total, # 风向频率 v_mean: group[wind_speed].mean(), # 扇区平均风速 v_cube_mean: (group[wind_speed] ** 3).mean(), # 扇区风速立方均值 }) return pd.DataFrame(stats) stats_df wind_rose_stats(df, sectors16) stats_df.head()sectors决定扇区宽度freq是风向频率v_cube_mean是后面算风能密度的关键输入。这里特别提醒别对平均风速先取三次方再算那会系统性低估因为三次方是个凸函数风速波动越激烈E(v³)和(E(v))³的差距越大。4.3 有效风能密度公式、空气密度修正与代码风能资源评估的核心指标是风能密度WPD公式是 WPD 0.5 × ρ × E(v³)其中ρ是空气密度E(v³)是风速立方的期望。这里的E(v³)必须按每个样本的风速立方取平均而不是“平均风速的三次方”差别在强湍流下能到15%以上。空气密度ρ也别直接拿1.225用。高原风电场的实际密度可能只有1.05~1.10直接套标准密度会高估5%~15%的风能密度。如果SCADA里有机舱温度和气压优先实测计算ρ P / (R × T)P为气压PaT为开尔文温度R287.05 J/(kg·K)。def air_density(temp_cNone, pressure_paNone, altitudeNone): 计算空气密度优先用SCADA实测温度气压其次用海拔公式估算 if temp_c is not None and pressure_pa is not None: T temp_c 273.15 return pressure_pa / (287.05 * T) if altitude is not None: # 标准大气压海拔修正无实测时用 p0 101325.0 T0 288.15 g 9.80665 M 0.0289644 R_univ 8.31447 p p0 * (1 - 0.0065 * altitude / T0) ** (g * M / (R_univ * 0.0065)) return p / (R_univ * T0) return 1.225 # 海平面标准密度 def wpd(speed_series, rho): 风能密度 0.5 * rho * mean(v^3) return 0.5 * rho * (speed_series ** 3).mean() rho air_density(temp_c15.0, pressure_pa96000.0) # 或者: rho air_density(altitude500) wpd_all wpd(df[wind_speed], rho) print(f不考虑风向的全场风能密度: {wpd_all:.1f} W/m²)如果数据里有逐小时的温度和气压记录空气密度应该逐条计算再乘到每条样本上而不是先对温度和气压取平均再算一个密度。逐样本计算的结果更稳代码里为了简洁只演示了单值版本。此外论文里常会附带各扇区的威布尔拟合参数用scipy.stats的weibull_min做极大似然估计注意固定location0否则拟合出的形状参数会和论文对不上。4.4 考虑风向和不考虑风向差异到底在哪很多读者会问WPD公式里根本没有风向考虑了风向结果会变化吗数学上不会变——全体样本的风速立方均值恰好等于各扇区风速立方均值的频率加权和。所以“考虑风速风向”改变的不是全场风能密度这个数字而是它的空间分布和场址适用性。考虑风向之后你能看到三个论文真正关心的结论主导风向扇区的风能占比。好的场址通常60%以上的风能集中在两三个相邻扇区这决定了风机排布方案里主风向方向的行间距。风向频率和风速频率是否一致。有些扇区风向频率高但风速低风能贡献反而小只看风向频率会被带偏。同一风速下不同风向的功率曲线是否一致。偏航误差较大的扇区功率曲线会整体右移这直接影响发电量预测的修正系数。复现时做一张扇区WPD对比表就够了def sector_wpd_table(df, rho, sectors16): 输出每个扇区的风向频率与风能密度贡献 rose wind_rose_stats(df, sectors) rose[wpd_sector] 0.5 * rho * rose[v_cube_mean] rose[wpd_contribution] rose[freq] * rose[wpd_sector] # 全场WPD等于各扇区贡献之和 wpd_total rose[wpd_contribution].sum() rose[energy_share] rose[wpd_contribution] / wpd_total * 100.0 return rose.sort_values(energy_share, ascendingFalse) table sector_wpd_table(df, rho, sectors16) print(table.head(5))energy_share就是每个扇区的风能占比主导风向扇区一目了然。论文里的风向玫瑰图和风能玫瑰图本质上就是这张表的不同画法。5. 论文复现中的五个常见坑现象、原因与解决办法5.1 清洗后曲线仍有大量离群点现象跑完bin-IQR散点图上依然有很多明显偏离主带的点。原因最常见是分箱太宽。bin_width取2m/s时一个箱里混入了好几个风速区间的数据IQR被拉得很宽离群点自然溜进去。另一种情况是稀疏箱没有单独处理高风速段直接套IQR离群点被保留下来。解决先把bin_width调回0.5m/s再检查每个箱的样本量给样本量小于min_samples的箱单独走DBSCAN分支。5.2 风向扇区评估和论文差异巨大现象算出来的主导风向扇区角度和论文里完全对不上。原因风向基准方向不一致最常翻车。SCADA里风向可能是机舱相对方位角0°机舱朝向而论文用的是气象风0°正北另外风向字段可能被记录成了度/10或其他编码开工前先确认量纲。解决先看风向直方图是否有规律再对照时间戳和机舱朝向推导基准角。如果始终对不上可以把风向整体旋转一个固定角度用实测风速与功率的相关性来校准偏航角。5.3 清洗前后平均风速变化很小但功率评估变化很大现象清洗前平均风速6.8m/s清洗后6.7m/s看着几乎没变可是WPD从380W/m²掉到了310W/m²。原因这是功率与风速三次方关系的必然结果。被清洗的异常点大多集中在低风速段和高风速段两个尾部这些点对一次方平均风速影响很小对三次方平均风速影响巨大。这不是bug是正常现象。解决报告里分别列出原始和清洗后的平均风速、风速立方均值、WPD三个指标不要只报一个数。5.4 DBSCAN把正常低风速点删了现象DBSCAN清洗后切入风速附近的正常点大量消失。原因特征没归一化或者归一化方式不对。风速范围0~30功率范围0~2000直接用原始值算距离功率维度会完全主导邻域判定正常点被误判为噪声。解决改RobustScaler做归一化eps在归一化空间里通常设为0.1~0.3。如果低风速段数据太密可以单独缩小eps再试。5.5 复现结果图“不像论文”现象论文里的功率曲线是一条很细的带子自己画出来却是一大团粗散的云。原因论文里的曲线大多是清洗后按风速分箱取平均得到的汇总曲线或加上了置信区间带你只是在散点图上叠加了原始数据点没有做二次汇总。解决清洗后的数据按0.5m/s分箱算每一箱的P50和90%区间画出一条带置信带的曲线。画完这张图和论文对比清洗质量一眼就见分晓。6. 进阶把你的清洗与评估流程做成可复用的一行命令6.1 把整条链路封装成流水线把前面的步骤打包成一个函数以后换一个风电场的数据只需要改参数。这个函数接收原始DataFrame返回清洗后的数据和扇区评估表。def clean_and_assess(df, rated_power2000.0, cut_in3.0, bin_width0.5, iqr_k1.5, sectors16): 一站式物理筛查 - bin-IQR - 扇区风资源评估 df_clean physical_filter(df, rated_power, cut_in) df_clean bin_iqr_filter(df_clean, bin_widthbin_width, kiqr_k) rho air_density(altitude0.0) # 按场址实际海拔修改 table sector_wpd_table(df_clean, rho, sectorssectors) return df_clean, table6.2 用拟合优度验证清洗质量判断清洗流程是否真的可用可以拿清洗后的数据重新拟合功率曲线看拟合误差是否比清洗前小。常见做法是用分段多项式或非线性最小二乘拟合分箱平均点再计算RMSE。如果拟合RMSE没有下降说明清洗参数过松或过紧需要回头检查散点图。我自己的习惯是每调一组清洗参数就记录“保留率、分箱宽度、RMSE”三个数排成一张小表。换数据后先跑一遍默认参数再用散点图检查每个清洗步骤的中间结果。清洗参数没有绝对正确的答案论文里的取值只适用于他们的机组和场址你的风机叶片、地形、尾流都不一样照抄参数一定会出偏差。希望这篇能帮你把数据清洗和风资源评估这条链路走通少走些弯路。本文还有配套的精品资源点击获取