风速均一化订正实战：从断点检测到分位数匹配的Python全流程

简介这份资料包面向气象与水文领域的研究人员及数据分析学习者聚焦最大风速序列的均一化订正问题帮助消除因仪器更换、测量方法调整或站点迁移带来的系统性偏差使不同站点、不同时期的风速记录具备可比性。包内共3个文件包含1个Python脚本、1个CSV数据文件和1份PDF说明文档压缩包约377KB分别对应订正算法实现、示例风速数据与代码使用说明便于直接运行与对照理解。目前已有591人学习下载。通过脚本与数据读者可完整走通数据读取、缺失值与异常值预处理、订正因子计算、均一化序列生成及结果验证等环节并借助可视化对比订正前后的风速变化趋势掌握Thorne-Wyatt、Renfrew、HOM等方法的实现思路为气候研究中的非均一数据处理提供可复用的实践参考。1. 风速均一化订正为什么最大风速序列不能直接拿来做趋势分析如果你手上有某个气象站 1960 年至今的日最大风速序列直接画一条趋势线大概率会得到一个“风速在显著下降”的结论。但这个结论很可能是假的。原因不在天气而在观测系统本身测风仪器的型号换过、塔高变过、周边建筑物长起来了、观测时次从 4 次变成 24 次、甚至站址搬过。这些非气候因素造成的跳变会淹没真实的气候信号。风速均一化订正要解决的就是把这些“人为断点”识别出来并订正掉让订正后的序列只保留气候意义上的变化。这篇笔记以日最大风速为例把均一化订正的完整链路走一遍从数据准备、断点检测、到订正量计算和效果验证每一步都给可复现的代码和参数说明。适合手里有长序列风速数据、准备做趋势分析或极值统计的气象水文从业者也适合刚接触均一化、想找一个完整案例上手的人。需要先说清楚一个边界均一化订正不是“把数据修好看”它是有假设的。核心假设是——参考序列邻近站或再分析与被检站经历相同的气候变化但参考序列本身是均一的。如果参考站也换过仪器那订正结果就是错的。所以选参考站这件事比后面跑什么检测算法都重要。我一般会花一半时间在选站和元数据核对上剩下的一半才交给算法。2. 均一化订正的技术路线从元数据到统计检验怎么选2.1 先搞清楚断点从哪来元数据优先统计检验兜底风速序列的非均一性来源按影响从大到小排通常是这几类仪器更换尤其从风杯换成超声启动风速阈值变了、观测高度变化10m 换到 2m 或反之、站址迁移、周边环境变化新建高楼、树木生长、观测时次和统计方法变化日最大风速的统计窗口变了。这些信息最可靠的来源是台站元数据——历史沿革表、仪器更换记录、站址变更文件。但现实是很多老站的元数据残缺甚至根本没有。这时候才轮到统计检验上场。统计检验的逻辑是如果序列在某点前后相对于参考序列的差值发生了显著变化那这个点就可能是断点。注意是“相对于参考序列”不是看序列本身的均值跳变——因为气候变化本身也会造成均值变化只有“差值”的变化才指向非均一性。常见做法是元数据和统计检验结合先用元数据圈定可疑年份再用统计检验确认或者反过来统计检验找出候选断点再回元数据里找解释。两者对上了断点可信度就高。2.2 四种主流方法的选择SNHT、PMT、MASH、RHtest 各适合什么场景风速均一化订正常用的方法有这么几类选哪个取决于你的数据条件和目标方法全称/来源适合场景局限SNHT标准正态均一性检验单断点、序列较长、参考序列质量好多断点场景容易漏检PMTPairwise Multiple Test多断点、有多个参考站计算量大参考站需均一MASHMultiple Analysis of Series for Homogenization多站联合检测适合台站网实现复杂对参考站数量有要求RHtestR 语言 rhTest 包单站、可结合元数据、支持多种检验需要 R 环境参数较多我一般会这样选如果只有一个站、参考序列是再分析资料用 RHtest 的 PMT 或 SNHT 模式如果是区域台站网、有多个邻近站用 MASH 或 PMT 做多站联合检测。风速这个变量比气温噪声大单站检测的漏检率不低所以有条件尽量用多站方法。2.3 参考序列怎么选邻近站、再分析、还是区域平均参考序列的质量直接决定订正成败。三种常见选择邻近站选距离近、海拔相近、气候背景一致的站。要求参考站本身经过均一化检验或者至少有完整元数据证明没换过仪器。距离一般控制在 50km 以内山区要更近。再分析资料比如 ERA5 的 10m 风速。优点是时空连续、没有断点缺点是分辨率粗对局地风速的刻画能力有限尤其是复杂地形。用再分析做参考时通常要先做尺度匹配——把再分析风速插值到站点或者用回归建立两者关系。区域平均把多个邻近站平均成一个区域序列。平均能削弱单站噪声但如果区域内有站也非均一会污染参考序列。我的习惯是优先用经过检验的邻近站没有就用再分析两者都有就做交叉验证——用邻近站订正一遍用再分析订正一遍看结果差多少。差得大说明参考序列本身有问题得回头查。3. 用 Python 跑通最大风速均一化订正的最小流程3.1 数据准备日最大风速序列的读取与质控假设你拿到的是 CSV 格式的日最大风速数据列包括日期、风速值、站点号。第一步不是直接跑检测而是质控。风速数据的常见问题缺测、异常大值比如台风天记录到 50m/s 但实际是仪器故障、连续相同值仪器卡死、单位不统一m/s 和 km/h 混用。import pandas as pd import numpy as np # 读取日最大风速数据 df pd.read_csv(daily_max_wind.csv, parse_dates[date]) df df.sort_values(date).reset_index(dropTrue) # 基本质控 # 1. 风速为负或超过合理上限的标记为缺测 df.loc[(df[wind_speed] 0) | (df[wind_speed] 60), wind_speed] np.nan # 2. 连续相同值超过5天的标记为可疑仪器卡死 same_count (df[wind_speed] df[wind_speed].shift()).astype(int) same_group (same_count 0).cumsum() group_sizes df.groupby(same_group)[wind_speed].transform(size) df.loc[group_sizes 5, wind_speed] np.nan # 3. 统计缺测率 missing_rate df[wind_speed].isna().mean() print(f缺测率: {missing_rate:.2%}) # 4. 按月统计检查是否有整月缺测 monthly_count df.set_index(date)[wind_speed].resample(M).count() print(monthly_count[monthly_count 20])这段代码做了四件事把物理上不可能的值置为缺测、识别仪器卡死造成的连续相同值、统计整体缺测率、检查是否有整月缺测。参数说明风速上限 60m/s 是通用阈值沿海台风影响区可以放宽到 70连续相同值阈值 5 天是经验值干燥地区可以放宽到 7 天。缺测率超过 20% 的年份后续订正要谨慎因为断点检测对缺测敏感。3.2 断点检测用 RHtest 的 SNHT 模式找候选断点RHtest 是加拿大环境部开发的均一化检验工具有 R 包和命令行版本。这里用 Python 调用 R 的方式演示也可以直接用 R 跑。核心是rhTest函数支持 SNHT、PMT 等多种检验。import subprocess import os # 准备 RHtest 输入文件两列第一列参考序列第二列待检序列 # 这里假设参考序列是邻近站已经过质控 ref pd.read_csv(ref_station.csv, parse_dates[date]).set_index(date)[wind_speed] target df.set_index(date)[wind_speed] # 对齐日期 combined pd.concat([ref, target], axis1, joininner).dropna() combined.columns [ref, target] combined.to_csv(rhtest_input.txt, sep , indexFalse, headerFalse) # 调用 RHtest假设已安装 R 和 rhTest 包 r_script library(rhTest) data - read.table(rhtest_input.txt) # SNHT 检验返回候选断点 result - rhTest(data$V2, data$V1, testSNHT) print(result$breaks) write.csv(result$breaks, breaks.csv, row.namesFALSE) with open(run_rhtest.R, w) as f: f.write(r_script) subprocess.run([Rscript, run_rhtest.R], checkTrue) breaks pd.read_csv(breaks.csv) print(breaks)逻辑说明RHtest 需要两列输入——参考序列和待检序列按日期对齐。rhTest函数返回候选断点及其显著性。参数说明testSNHT指定用标准正态均一性检验适合单断点如果怀疑多断点改成testPMT。显著性水平默认 0.05可以通过p.value参数调整。注意RHtest 对缺测敏感输入前要确保两列都没有缺测或者用插值补齐——但插值本身会引入误差缺测多的年份建议单独处理。3.3 订正量计算用分位数匹配做逐日订正找到断点后下一步是计算订正量。风速的订正不能简单加减均值差因为风速分布是偏态的不同分位数的偏差不一样。常用做法是分位数匹配Quantile Matching在断点前后各取一段窗口分别计算参考序列和待检序列的分位数关系然后把这个关系外推到整个时段。def quantile_matching_correction(target, ref, break_date, window5): 用分位数匹配计算订正量 target: 待检序列Series索引为日期 ref: 参考序列Series索引为日期 break_date: 断点日期 window: 断点前后各取多少年做拟合 break_year pd.Timestamp(break_date).year # 断点前窗口 pre_mask (target.index.year break_year - window) (target.index.year break_year) # 断点后窗口 post_mask (target.index.year break_year) (target.index.year break_year window) # 计算断点前后的分位数关系 quantiles np.arange(0.05, 1.0, 0.05) pre_target_q target[pre_mask].quantile(quantiles) pre_ref_q ref[pre_mask].quantile(quantiles) post_target_q target[post_mask].quantile(quantiles) post_ref_q ref[post_mask].quantile(quantiles) # 断点前的比值关系target/ref pre_ratio pre_target_q / pre_ref_q post_ratio post_target_q / post_ref_q # 订正系数把断点后的关系调整到断点前 correction_factor pre_ratio / post_ratio # 对断点后的数据做逐日订正 corrected target.copy() post_all target.index pd.Timestamp(break_date) # 用插值把分位数订正系数映射到每个值 for i, val in target[post_all].items(): # 找到 val 在断点后分位数中的位置 q_rank (post_target_q val).mean() # 插值得到对应的订正系数 factor np.interp(q_rank, quantiles, correction_factor) corrected.loc[i] val * factor return corrected # 应用订正 corrected_series quantile_matching_correction(target, ref, 1985-01-01, window5)逻辑说明这段代码的核心思想是——如果断点前后 target 和 ref 的比值关系发生了变化那这个变化就是非均一性造成的需要把断点后的比值调整回断点前的水平。参数说明window5表示用断点前后各 5 年做拟合窗口太短拟合不稳太长会混入气候变化信号一般 5 到 10 年。quantiles从 0.05 到 0.95步长 0.05覆盖了风速的主要分布范围。注意分位数匹配假设断点前后的偏差是乘性的如果偏差是加性的应该用差值而不是比值。风速一般用乘性更合理因为高风速的绝对偏差通常更大。3.4 效果验证订正前后序列对比与趋势检验订正完不能直接信要做验证。验证分两步一是看订正后断点是否消除二是看订正对趋势的影响。import matplotlib.pyplot as plt from scipy import stats # 1. 订正前后对比图 fig, axes plt.subplots(2, 1, figsize(12, 8)) axes[0].plot(target.index, target.values, label原始序列, alpha0.7) axes[0].plot(corrected_series.index, corrected_series.values, label订正后, alpha0.7) axes[0].axvline(pd.Timestamp(1985-01-01), colorr, linestyle--, label断点) axes[0].legend() axes[0].set_ylabel(日最大风速 (m/s)) # 2. 订正前后趋势对比 def trend_test(series): 用 Mann-Kendall 检验计算趋势 n len(series) s 0 for i in range(n-1): for j in range(i1, n): s np.sign(series.iloc[j] - series.iloc[i]) var_s n*(n-1)*(2*n5)/18 z (s - np.sign(s)) / np.sqrt(var_s) if s ! 0 else 0 p 2 * (1 - stats.norm.cdf(abs(z))) # Sens slope slopes [] for i in range(n-1): for j in range(i1, n): slopes.append((series.iloc[j] - series.iloc[i]) / (j - i)) slope np.median(slopes) return slope, p slope_orig, p_orig trend_test(target.dropna()) slope_corr, p_corr trend_test(corrected_series.dropna()) print(f原始序列趋势: {slope_orig:.4f} m/s/年, p{p_orig:.4f}) print(f订正后趋势: {slope_corr:.4f} m/s/年, p{p_corr:.4f})逻辑说明Mann-Kendall 是非参数趋势检验不要求数据正态分布适合风速这种偏态变量。Sens slope 给出趋势的稳健估计。参数说明p 值小于 0.05 认为趋势显著。如果订正前后趋势方向或显著性发生根本变化说明订正起了作用——但也要警惕订正过度。我一般会对比订正前后的趋势如果订正后趋势从显著变不显著或者斜率变化超过 50%就要回头检查断点是否找多了。4. 风速均一化订正的避坑清单五个血泪教训4.1 坑一参考站自己就不均一订正结果全歪现象订正后序列在断点处确实平滑了但整体趋势变得很奇怪和邻近区域其他站对不上。原因参考站本身在同期也换过仪器或迁过站但没做均一化检验。用非均一的参考去订正待检站相当于用一把不准的尺子去量另一把尺子。解决参考站必须先做均一化检验。如果参考站也有断点要么换站要么先订正参考站。我一般会要求参考站的元数据完整或者至少用 RHtest 跑一遍确认没有显著断点。4.2 坑二断点检测对缺测太敏感缺测多的年份误报现象检测出一堆断点但元数据里这些年份没有任何仪器或站址变化。原因缺测导致序列的统计特性在缺测前后发生变化被算法误判为断点。风速数据在早期1960-1980缺测率普遍偏高。解决检测前先统计逐月缺测率缺测率超过 30% 的年份标记出来检测时排除或单独处理。RHtest 有处理缺测的选项但效果有限。我的做法是缺测率高的时段不参与断点检测但订正时仍然处理——用邻近时段的关系外推。4.3 坑三订正窗口选太长把气候信号也订正掉了现象订正后序列的趋势比原始序列还大或者趋势方向反转。原因分位数匹配的窗口如果取 10 年以上窗口内本身包含气候变化趋势导致订正系数里混入了气候信号订正时把真实趋势也放大了。解决窗口一般取 5 年最多不超过 8 年。如果断点密集窗口还要缩短。另外订正前先对序列做去趋势处理订正后再把趋势加回去——这个做法更稳妥但实现复杂一些。4.4 坑四多断点场景用单断点方法漏检严重现象序列明显有多个跳变但 SNHT 只检测出一个断点。原因SNHT 是单断点检验序列有多个断点时它只能找到最显著的那个其他断点被掩盖。解决用 PMT 或 MASH 做多断点检测。如果只能用 SNHT就分段检测——先找最显著的断点把序列分成两段再在每段里继续找直到没有显著断点。但分段检测会累积误差断点多了要谨慎。4.5 坑五订正后不做独立验证自说自话现象订正报告里只有订正前后的对比图没有独立验证。原因对比图只能说明订正起了作用不能说明订正对了。订正可能过度也可能不足。解决至少做两种独立验证。一是用另一套参考序列比如再分析重复订正看结果是否一致二是留出部分时段不参与订正用订正后的关系去预测看预测误差。我一般会用 ERA5 做交叉验证——如果邻近站订正和 ERA5 订正的结果差在 10% 以内就认为订正可信。5. 进阶技巧用元数据约束断点检测把误报率压下来前面讲的流程是纯统计的但实际业务里元数据才是最强的约束。我现在的做法是先用元数据圈定“可疑年份”——仪器更换、站址迁移、观测时次变化的年份然后在这些年份附近做统计检验而不是全序列盲扫。这样误报率能降一大半。具体操作把元数据整理成一张表列包括年份、变更类型、变更描述。然后用 Python 把这张表和断点检测结果做匹配。# 元数据表 metadata pd.DataFrame({ year: [1975, 1985, 1998, 2005], change_type: [仪器更换, 站址迁移, 观测时次变化, 仪器更换], description: [风杯换型, 搬迁至新址, 4次改24次, 换超声风速仪] }) # 断点检测结果 detected_breaks pd.DataFrame({ break_date: [1975-06-01, 1985-03-01, 2005-08-01], p_value: [0.01, 0.03, 0.02] }) detected_breaks[year] pd.to_datetime(detected_breaks[break_date]).dt.year # 匹配断点年份与元数据变更年份相差不超过1年 matched detected_breaks.merge(metadata, onyear, howleft) matched[confirmed] matched[change_type].notna() print(matched)逻辑说明这段代码把统计检测出的断点和元数据变更记录做匹配。如果断点年份附近有元数据变更记录这个断点就是“有解释的”可信度高如果没有就是“无解释的”需要进一步排查——可能是参考序列的问题也可能是气候变化造成的假断点。参数说明匹配窗口设为 ±1 年因为元数据记录的年份和实际变更时间可能有偏差。如果元数据精确到月窗口可以缩小到 ±3 个月。还有一个技巧对“无解释的断点”不要急着订正先检查参考序列在同期有没有异常。如果参考序列在同期也有跳变那这个断点很可能是参考序列的问题不是待检站的。这时候应该换参考序列而不是硬订正。最后说一个我自己的习惯每次订正完我都会把订正前后的序列、断点位置、元数据变更记录画在一张图上人工过一遍。算法再先进也不如人眼扫一遍来得可靠。有些断点算法觉得显著但你看一眼就知道是台风年造成的假信号——这种就得手动排除。均一化订正这件事算法是工具判断在人。希望帮到你。本文还有配套的精品资源点击获取