简介这份资料面向气象与水文领域的研究人员及数据分析学习者聚焦最大风速的均一化订正问题帮助消除因仪器更换、测量方法调整或站点迁移带来的系统性偏差使不同站点与时期的风速记录具备可比性。压缩包共3个文件包含1个Python脚本、1个CSV数据文件和1份PDF说明文档整体约377KB分别对应订正算法实现、最大风速原始数据与代码使用说明。内容覆盖数据读取与预处理、缺失值和异常值处理、订正因子计算与应用、结果验证及可视化等环节读者可借助脚本与示例数据完整跑通均一化流程理解Thorne-Wyatt、Renfrew、HOM等方法的实现思路并将相关经验迁移到其他气象参数的均一化处理中。目前已有591人学习下载适合希望提升气象数据分析实操能力、需要可靠风速序列开展气候研究的读者参考。1. 风速均一化订正为什么最大风速序列不能直接拿来做趋势分析拿到一份气象站 30 年最大风速年极值序列直接画趋势线、跑 Mann-Kendall 检验大概率会得到一个“风速显著下降”的结论。但这个结论很可能是假的——不是因为风速没变而是因为测风仪器换了、站址迁了、周围盖楼了、观测规范改了。这些非气候因素造成的断点会让一条本来平稳的序列看起来像在“下降”或“上升”。风速均一化订正要解决的就是这个问题把台站历史最大风速序列里由非气候因素引起的系统性偏差识别出来并校正掉让订正后的序列真正反映气候信号。这套方法在气象水文领域属于基础但绕不开的环节尤其在做风资源评估、极值重现期推算、风灾风险区划时不订正的数据基本不能用。适合有基本 Python 或 R 能力、手头有台站历史风速数据、需要做长序列趋势分析的气象水文从业者和相关方向的研究生。2. 均一化订正的技术路线从 SNHT 到分位数映射2.1 为什么最大风速的订正比气温降水更难气温和降水的均一化订正方法已经比较成熟但最大风速有几个特殊之处导致不能直接套用。第一最大风速是极值统计量不是均值统计量。气温的均一化订正通常针对月均值或年均值样本量大、分布接近正态SNHT标准正态均一性检验和 Pettitt 检验都能较好地识别均值突变。但最大风速年极值一年只有一个值30 年也就 30 个点样本量极小检验功效天然不足。第二最大风速的物理分布是偏态的。年最大风速通常服从 Gumbel 分布、Weibull 分布或广义极值分布GEV不是正态分布。SNHT 的前提是序列近似正态直接用在最大风速上会出问题。常见做法是先做概率变换把最大风速转成接近正态的变量再做检验。第三最大风速的断点往往不是“均值突变”而是“方差突变”或“分布形态突变”。比如仪器从风杯换成超声风速仪可能均值变化不大但方差明显缩小——超声风速仪对小尺度湍流的响应不同。这时候只检验均值的 SNHT 就漏掉了。第四参考站的选择比气温降水更苛刻。气温和降水的空间相关性好几十公里外的参考站仍然可用。但最大风速受局地地形和地表粗糙度影响极大参考站必须非常近通常要求直线距离小于 50 km高差小于 200 m而且下垫面条件要相似。这就导致很多台站根本找不到合适的参考站只能做单站订正。2.2 方法选型SNHT、Pettitt、贝叶斯方法怎么选实际业务和科研中最大风速均一化订正的主流方法有以下几类方法适用场景优点局限SNHT标准正态均一性检验有参考站序列样本量≥20能定位断点位置可做多断点要求近似正态对极值序列需先变换Pettitt 检验单站序列无参考站非参数不要求分布只能检出一个断点对尾部变化不敏感贝叶斯突变点检测样本量小需要不确定性量化能给出断点概率分布计算量大先验选择影响结果分位数映射QM断点已识别需要校正能校正整个分布不只是均值需要断点前后都有足够样本多元线性回归残差检验有多个参考站能同时考虑多个因子参考站质量差时引入新偏差我一般会先用 Pettitt 检验做快速筛查再用 SNHT 配合参考站做确认。如果两种方法识别的断点位置一致基本可以确定。如果不一致就需要看元数据——台站有没有迁站记录、仪器更换记录、观测规范变更记录。元数据是均一化订正的“后悔药”没有元数据的时候统计方法的结论要打折扣。分位数映射是校正阶段的核心方法。它的思路是假设断点前后两个时段的最大风速分布之间存在一个变换关系用断点后的分布去拟合断点前的分布。具体做法有经验分位数映射empirical QM和参数化分位数映射parametric QM。经验 QM 直接用经验累积分布函数做映射不假设分布形式适合样本量较大的情况。参数化 QM 先拟合 GEV 或 Weibull 分布再用拟合的分布做映射适合样本量小的情况——最大风速年极值序列通常只有 30-50 个点参数化 QM 更稳妥。2.3 参考站选取的硬性条件和软性条件参考站选得好不好直接决定订正结果可不可信。硬性条件包括直线距离一般要求 50 km地形复杂地区 30 km高差 200 m山区可放宽到 300 m序列长度至少覆盖待订正站的全时段且缺测率 5%参考站自身经过均一化检验确认无断点软性条件包括下垫面相似都是平坦草地、都是城市站、都是山地站风向一致性最大风速的主导风向要相近相关性检验待订正站和参考站的年最大风速序列相关系数应 0.6否则参考站信息量不足如果找不到满足条件的参考站就退化为单站订正。单站订正只能依赖 Pettitt 检验和元数据可靠性会下降但总比不订正好。3. 用 Python 跑通最大风速均一化订正的最小流程3.1 数据准备与预处理假设手头有一个 CSV 文件包含台站号、年份、年最大风速m/s、最大风速出现日期。先做基本预处理。import pandas as pd import numpy as np from scipy import stats import matplotlib.pyplot as plt # 读取数据假设列名为 station_id, year, max_wind_speed, date df pd.read_csv(station_max_wind.csv, parse_dates[date]) # 筛选目标台站 target_station 54511 df_target df[df[station_id] target_station].copy() df_target df_target.sort_values(year).reset_index(dropTrue) # 检查缺测年份 full_years pd.DataFrame({year: range(df_target[year].min(), df_target[year].max() 1)}) df_target full_years.merge(df_target, onyear, howleft) # 标记缺测 df_target[is_missing] df_target[max_wind_speed].isna() print(f总年数: {len(df_target)}, 缺测年数: {df_target[is_missing].sum()}) # 简单插补线性插值仅适用于连续缺测不超过3年的情况 df_target[max_wind_speed_filled] df_target[max_wind_speed].interpolate( methodlinear, limit3, limit_directionboth) # 如果缺测太多考虑剔除该站或使用更复杂的插补方法 if df_target[is_missing].sum() len(df_target) * 0.1: print(警告缺测率超过10%订正结果可靠性下降)这段代码做了三件事读取数据、检查缺测、线性插补。参数说明limit3表示最多连续插补 3 年超过 3 年的缺测不插补因为线性插值在长缺测段会引入虚假趋势。limit_directionboth表示序列首尾的缺测也尝试插补。如果缺测率超过 10%建议换站或使用更复杂的时空插补方法。3.2 Pettitt 检验识别突变点Pettitt 检验是一种非参数突变点检测方法不要求数据服从特定分布适合最大风速这种偏态序列。def pettitt_test(x, alpha0.05): Pettitt 突变点检验 x: 输入序列一维数组 alpha: 显著性水平 返回: (突变点位置索引, p值, 是否显著) n len(x) U np.zeros(n) for t in range(1, n): U[t] U[t-1] np.sign(x[t] - np.median(x[:t])) * (t) # 简化版 # 标准 Pettitt 统计量 U_full np.zeros(n) for t in range(n): for i in range(t1): for j in range(t1, n): U_full[t] np.sign(x[i] - x[j]) K np.max(np.abs(U_full)) # 近似 p 值 p_value 2 * np.exp(-6 * K**2 / (n**3 n**2)) # 突变点位置 change_point np.argmax(np.abs(U_full)) return change_point, p_value, p_value alpha # 对最大风速序列做 Pettitt 检验 x df_target[max_wind_speed_filled].dropna().values cp, p, sig pettitt_test(x) print(f突变点年份: {df_target[year].iloc[cp]}, p值: {p:.4f}, 显著: {sig})Pettitt 检验的核心思想是如果序列存在突变点那么突变点前后的子序列之间的秩和差异会很大。U_full[t]统计的是前 t1 个点和后 n-t-1 个点之间的符号差异累积。K取最大绝对值p_value用近似公式计算。参数alpha0.05是显著性水平p 值小于 0.05 认为突变显著。注意Pettitt 检验只能检出一个突变点如果序列有多个断点需要分段检验或改用其他方法。3.3 SNHT 检验与参考站对比SNHT 需要参考站序列。假设已经选好参考站数据在同一 CSV 中。def snht_test(target, reference, alpha0.05): SNHT 均一性检验 target: 待检序列 reference: 参考序列 返回: (最大T值位置, p值, 是否显著) n len(target) # 计算比值序列或差值序列取决于变量类型 # 最大风速用差值更合适因为比值在风速接近0时不稳定 diff target - reference # 标准化 diff_std (diff - np.mean(diff)) / np.std(diff) T np.zeros(n) for t in range(1, n): # 前段均值和后段均值 mean1 np.mean(diff_std[:t]) mean2 np.mean(diff_std[t:]) T[t] t * mean1**2 (n - t) * mean2**2 T_max np.max(T) # 临界值近似n20时 # 更精确的临界值需要查表或蒙特卡洛模拟 T_critical 8.0 # 对应 alpha0.05 的近似值 change_point np.argmax(T) return change_point, T_max, T_max T_critical # 假设参考站数据在同一 DataFrame 中 ref_station 54512 df_ref df[df[station_id] ref_station].sort_values(year) # 对齐年份 merged df_target.merge(df_ref[[year, max_wind_speed]], onyear, suffixes(_target, _ref)) merged merged.dropna(subset[max_wind_speed_target, max_wind_speed_ref]) cp_snht, T_max, sig_snht snht_test( merged[max_wind_speed_target].values, merged[max_wind_speed_ref].values ) print(fSNHT 突变点年份: {merged[year].iloc[cp_snht]}, T{T_max:.2f}, 显著: {sig_snht})SNHT 的核心是构造一个统计量 T衡量突变点前后两段均值的差异。diff是待检站和参考站的差值序列标准化后消除量纲影响。T[t]在突变点处达到最大。临界值T_critical8.0是 n20 时的近似值更精确的做法是用蒙特卡洛模拟生成临界值表。注意SNHT 对序列两端的突变点检测能力较弱如果突变发生在序列开头或结尾 5 年内结果要谨慎。3.4 分位数映射校正识别出突变点后用分位数映射做校正。假设突变点在 1995 年1995 年之前是“旧”时段之后是“新”时段需要把旧时段校正到新时段的分布。def quantile_mapping(data_old, data_new, data_to_correct): 分位数映射校正 data_old: 旧时段数据参考分布 data_new: 新时段数据目标分布 data_to_correct: 需要校正的数据 返回: 校正后的数据 # 经验累积分布 old_sorted np.sort(data_old) new_sorted np.sort(data_new) # 对每个待校正值找到它在旧分布中的分位数 # 然后映射到新分布的对应分位数 corrected np.zeros_like(data_to_correct) for i, val in enumerate(data_to_correct): # 在旧分布中的分位数 p np.searchsorted(old_sorted, val) / len(old_sorted) p np.clip(p, 0.01, 0.99) # 避免极端分位数 # 在新分布中的对应值 idx int(p * len(new_sorted)) idx min(idx, len(new_sorted) - 1) corrected[i] new_sorted[idx] return corrected # 分段 break_year 1995 old_data merged[merged[year] break_year][max_wind_speed_target].values new_data merged[merged[year] break_year][max_wind_speed_target].values # 校正旧时段数据 corrected_old quantile_mapping(old_data, new_data, old_data) # 合并校正后的序列 corrected_series np.concatenate([corrected_old, new_data]) years np.concatenate([ merged[merged[year] break_year][year].values, merged[merged[year] break_year][year].values ]) # 可视化对比 fig, axes plt.subplots(1, 2, figsize(12, 4)) axes[0].plot(years, np.concatenate([old_data, new_data]), b-, label原始) axes[0].plot(years, corrected_series, r-, label校正后) axes[0].axvline(break_year, colorgray, linestyle--, label突变点) axes[0].set_xlabel(年份) axes[0].set_ylabel(最大风速 (m/s)) axes[0].legend() axes[0].set_title(序列对比) axes[1].hist(old_data, bins10, alpha0.5, label旧时段, densityTrue) axes[1].hist(new_data, bins10, alpha0.5, label新时段, densityTrue) axes[1].hist(corrected_old, bins10, alpha0.5, label校正后旧时段, densityTrue) axes[1].set_xlabel(最大风速 (m/s)) axes[1].legend() axes[1].set_title(分布对比) plt.tight_layout() plt.savefig(homogenization_result.png, dpi150) plt.show()分位数映射的逻辑是对旧时段的每个值找到它在旧分布中的分位数然后取新分布中同一分位数的值作为校正值。np.searchsorted用于快速定位分位数np.clip把分位数限制在 0.01-0.99 之间避免极端分位数导致校正值失真。参数说明break_year是突变点年份需要根据前面的检验结果设定。校正后的序列在突变点前后分布一致可以用于后续趋势分析。4. 避坑指南最大风速均一化订正中最容易翻车的 5 个地方4.1 坑一参考站选得太远订正后引入新偏差现象订正后的序列趋势和参考站趋势高度一致但和待订正站周边的其他站趋势不一致。原因参考站距离太远最大风速的空间相关性已经很低参考站的信息主要反映的是它自己的局地气候不是待订正站的气候信号。强行用远距离参考站做 SNHT会把参考站的局地特征“传染”给待订正站。解决参考站距离严格控制在 50 km 以内山区控制在 30 km 以内。如果找不到宁可做单站订正。单站订正虽然可靠性下降但不会引入虚假的空间信号。可以用多个参考站做交叉验证如果不同参考站给出的断点位置差异很大说明参考站信息不可靠。4.2 坑二忽略元数据统计断点和实际断点对不上现象Pettitt 检验和 SNHT 检验都显示 1995 年有突变但台站元数据里 1995 年没有任何仪器更换或迁站记录。反而 2003 年有明确的仪器更换记录但统计检验在 2003 年没有检出显著突变。原因统计检验检出的是“数据分布变化”不一定是“仪器变化”。1995 年的突变可能来自观测规范变更比如最大风速的统计时段从 10 分钟改为 2 分钟这种变更在元数据里可能没有详细记录。而 2003 年的仪器更换可能恰好没有改变最大风速的统计特性比如两种仪器在最大风速量级上响应一致。解决统计检验和元数据必须交叉验证。元数据是“金标准”统计检验是“辅助工具”。如果统计检验检出的断点在元数据中有对应记录优先采信。如果统计检验检出但元数据没有记录需要进一步排查——可能是元数据不完整也可能是统计检验的假阳性。可以用贝叶斯方法给出断点的概率分布而不是硬性判定“有”或“没有”。4.3 坑三分位数映射在样本量小时外推过度现象校正后的旧时段最大风速出现了不合理的极值比如校正后出现了 50 m/s 的风速但原始序列最大值只有 35 m/s。原因经验分位数映射在样本量小时极端分位数比如 0.95 以上的估计非常不稳定。旧时段可能只有 20 个点0.95 分位数对应的是第 19 个点稍微换个样本这个值就变了。如果新时段的 0.95 分位数恰好很大校正后就会产生虚假极值。解决样本量小于 30 时优先用参数化分位数映射。先拟合 GEV 分布再用拟合的分布做映射。GEV 分布的尾部有参数控制不会像经验分位数那样剧烈波动。如果一定要用经验分位数映射把分位数范围限制在 0.05-0.95 之间超出范围的值用参数化方法外推。4.4 坑四多断点序列只做单断点订正现象序列在 1985 年和 2005 年各有一个断点但只做了 1985 年的订正2005 年之后的序列仍然有偏差。原因Pettitt 检验和 SNHT 检验默认只检出一个断点。如果序列有多个断点单断点方法只能找到最显著的那个其他断点被忽略。解决先用单断点方法找到最显著的断点把序列分成两段然后在每段内再做单断点检验。重复这个过程直到没有显著断点。这就是“分段检验”的思路。更严谨的做法是用多断点方法比如贝叶斯多断点模型或动态规划方法。但多断点方法计算量大而且断点越多不确定性越大。实际业务中如果断点超过 3 个建议直接剔除该站因为订正后的序列可靠性已经很低了。4.5 坑五订正后不做独立验证现象订正后的序列趋势合理但用来做极值重现期推算时得到的 50 年一遇风速和周边站差异很大。原因订正只保证了序列内部的均一性没有保证序列和周边站的时空一致性。如果订正过程中引入了偏差序列内部看起来均一但和周边站对比就会暴露问题。解决订正后必须做独立验证。常用方法有① 和周边未订正站做空间一致性检验订正后的序列应该和周边站的相关性更高② 用订正后的序列做极值推算和用原始序列、周边站序列的结果对比差异应该在合理范围内③ 如果有可能用独立时段的数据做交叉验证——比如用 1980-2000 年数据建模用 2001-2020 年数据验证。5. 订正后的序列怎么用趋势检验与极值推算的衔接订正完序列下一步通常是做趋势分析或极值推算。这两件事对订正质量的要求不同衔接时需要注意几个技巧。趋势分析对订正质量最敏感。如果订正不彻底残留的断点会直接污染趋势估计。我一般会做“双保险”先用 Mann-Kendall 检验做非参数趋势检验再用 Sens slope 估计趋势幅度。Mann-Kendall 检验不要求正态分布对最大风速这种偏态序列比较稳健。Sens slope 比线性回归斜率更抗 outliers。如果两种方法给出的趋势方向一致基本可信。如果不一致说明序列里还有未识别的断点或 outliers。import pymannkendall as mk # 对订正后的序列做 Mann-Kendall 趋势检验 result mk.original_test(corrected_series) print(f趋势: {result.trend}, p值: {result.p:.4f}, Sens slope: {result.slope:.4f}) # 如果 p 0.05趋势显著 # Sens slope 的单位是 m/s per year极值推算对订正质量的要求更高因为极值推算依赖分布的尾部。订正后的序列如果尾部被扭曲重现期估计会严重偏差。我一般会先用订正后的序列拟合 GEV 分布然后用轮廓似然法估计重现期。轮廓似然法比矩估计法更稳健尤其在样本量小时。from scipy.stats import genextreme as gev # 拟合 GEV 分布 params gev.fit(corrected_series) shape, loc, scale params # 推算 50 年一遇最大风速 return_period 50 p 1 - 1/return_period wind_50 gev.ppf(p, shape, loc, scale) print(f50年一遇最大风速: {wind_50:.2f} m/s) # 轮廓似然法估计置信区间简化版 # 实际应用中建议用专门的极值分析包如 pyextremes这里有个血泪经验订正后的序列在做极值推算时不要直接用原始序列的极值。订正可能会改变极值的大小尤其是分位数映射校正后旧时段的极值可能被“拉”到新时段的分布上。如果旧时段的极值被拉得过高重现期估计会偏大。我一般会对比订正前后极值推算结果如果差异超过 20%就需要回头检查订正过程。最后一个技巧订正后的序列建议保留“订正标记”。在数据表里加一列is_corrected标记哪些年份的数据被校正过。这样后续分析时可以区分“原始观测”和“校正值”在做敏感性分析时能快速切换。这个习惯帮我省了很多后悔药——有一次合作方质疑订正结果我直接拉出标记列对比了订正前后的趋势问题当场定位。希望帮到你。本文还有配套的精品资源点击获取 最大风速均一化订正：从SNHT到分位数映射的Python实战
简介这份资料面向气象与水文领域的研究人员及数据分析学习者聚焦最大风速的均一化订正问题帮助消除因仪器更换、测量方法调整或站点迁移带来的系统性偏差使不同站点与时期的风速记录具备可比性。压缩包共3个文件包含1个Python脚本、1个CSV数据文件和1份PDF说明文档整体约377KB分别对应订正算法实现、最大风速原始数据与代码使用说明。内容覆盖数据读取与预处理、缺失值和异常值处理、订正因子计算与应用、结果验证及可视化等环节读者可借助脚本与示例数据完整跑通均一化流程理解Thorne-Wyatt、Renfrew、HOM等方法的实现思路并将相关经验迁移到其他气象参数的均一化处理中。目前已有591人学习下载适合希望提升气象数据分析实操能力、需要可靠风速序列开展气候研究的读者参考。1. 风速均一化订正为什么最大风速序列不能直接拿来做趋势分析拿到一份气象站 30 年最大风速年极值序列直接画趋势线、跑 Mann-Kendall 检验大概率会得到一个“风速显著下降”的结论。但这个结论很可能是假的——不是因为风速没变而是因为测风仪器换了、站址迁了、周围盖楼了、观测规范改了。这些非气候因素造成的断点会让一条本来平稳的序列看起来像在“下降”或“上升”。风速均一化订正要解决的就是这个问题把台站历史最大风速序列里由非气候因素引起的系统性偏差识别出来并校正掉让订正后的序列真正反映气候信号。这套方法在气象水文领域属于基础但绕不开的环节尤其在做风资源评估、极值重现期推算、风灾风险区划时不订正的数据基本不能用。适合有基本 Python 或 R 能力、手头有台站历史风速数据、需要做长序列趋势分析的气象水文从业者和相关方向的研究生。2. 均一化订正的技术路线从 SNHT 到分位数映射2.1 为什么最大风速的订正比气温降水更难气温和降水的均一化订正方法已经比较成熟但最大风速有几个特殊之处导致不能直接套用。第一最大风速是极值统计量不是均值统计量。气温的均一化订正通常针对月均值或年均值样本量大、分布接近正态SNHT标准正态均一性检验和 Pettitt 检验都能较好地识别均值突变。但最大风速年极值一年只有一个值30 年也就 30 个点样本量极小检验功效天然不足。第二最大风速的物理分布是偏态的。年最大风速通常服从 Gumbel 分布、Weibull 分布或广义极值分布GEV不是正态分布。SNHT 的前提是序列近似正态直接用在最大风速上会出问题。常见做法是先做概率变换把最大风速转成接近正态的变量再做检验。第三最大风速的断点往往不是“均值突变”而是“方差突变”或“分布形态突变”。比如仪器从风杯换成超声风速仪可能均值变化不大但方差明显缩小——超声风速仪对小尺度湍流的响应不同。这时候只检验均值的 SNHT 就漏掉了。第四参考站的选择比气温降水更苛刻。气温和降水的空间相关性好几十公里外的参考站仍然可用。但最大风速受局地地形和地表粗糙度影响极大参考站必须非常近通常要求直线距离小于 50 km高差小于 200 m而且下垫面条件要相似。这就导致很多台站根本找不到合适的参考站只能做单站订正。2.2 方法选型SNHT、Pettitt、贝叶斯方法怎么选实际业务和科研中最大风速均一化订正的主流方法有以下几类方法适用场景优点局限SNHT标准正态均一性检验有参考站序列样本量≥20能定位断点位置可做多断点要求近似正态对极值序列需先变换Pettitt 检验单站序列无参考站非参数不要求分布只能检出一个断点对尾部变化不敏感贝叶斯突变点检测样本量小需要不确定性量化能给出断点概率分布计算量大先验选择影响结果分位数映射QM断点已识别需要校正能校正整个分布不只是均值需要断点前后都有足够样本多元线性回归残差检验有多个参考站能同时考虑多个因子参考站质量差时引入新偏差我一般会先用 Pettitt 检验做快速筛查再用 SNHT 配合参考站做确认。如果两种方法识别的断点位置一致基本可以确定。如果不一致就需要看元数据——台站有没有迁站记录、仪器更换记录、观测规范变更记录。元数据是均一化订正的“后悔药”没有元数据的时候统计方法的结论要打折扣。分位数映射是校正阶段的核心方法。它的思路是假设断点前后两个时段的最大风速分布之间存在一个变换关系用断点后的分布去拟合断点前的分布。具体做法有经验分位数映射empirical QM和参数化分位数映射parametric QM。经验 QM 直接用经验累积分布函数做映射不假设分布形式适合样本量较大的情况。参数化 QM 先拟合 GEV 或 Weibull 分布再用拟合的分布做映射适合样本量小的情况——最大风速年极值序列通常只有 30-50 个点参数化 QM 更稳妥。2.3 参考站选取的硬性条件和软性条件参考站选得好不好直接决定订正结果可不可信。硬性条件包括直线距离一般要求 50 km地形复杂地区 30 km高差 200 m山区可放宽到 300 m序列长度至少覆盖待订正站的全时段且缺测率 5%参考站自身经过均一化检验确认无断点软性条件包括下垫面相似都是平坦草地、都是城市站、都是山地站风向一致性最大风速的主导风向要相近相关性检验待订正站和参考站的年最大风速序列相关系数应 0.6否则参考站信息量不足如果找不到满足条件的参考站就退化为单站订正。单站订正只能依赖 Pettitt 检验和元数据可靠性会下降但总比不订正好。3. 用 Python 跑通最大风速均一化订正的最小流程3.1 数据准备与预处理假设手头有一个 CSV 文件包含台站号、年份、年最大风速m/s、最大风速出现日期。先做基本预处理。import pandas as pd import numpy as np from scipy import stats import matplotlib.pyplot as plt # 读取数据假设列名为 station_id, year, max_wind_speed, date df pd.read_csv(station_max_wind.csv, parse_dates[date]) # 筛选目标台站 target_station 54511 df_target df[df[station_id] target_station].copy() df_target df_target.sort_values(year).reset_index(dropTrue) # 检查缺测年份 full_years pd.DataFrame({year: range(df_target[year].min(), df_target[year].max() 1)}) df_target full_years.merge(df_target, onyear, howleft) # 标记缺测 df_target[is_missing] df_target[max_wind_speed].isna() print(f总年数: {len(df_target)}, 缺测年数: {df_target[is_missing].sum()}) # 简单插补线性插值仅适用于连续缺测不超过3年的情况 df_target[max_wind_speed_filled] df_target[max_wind_speed].interpolate( methodlinear, limit3, limit_directionboth) # 如果缺测太多考虑剔除该站或使用更复杂的插补方法 if df_target[is_missing].sum() len(df_target) * 0.1: print(警告缺测率超过10%订正结果可靠性下降)这段代码做了三件事读取数据、检查缺测、线性插补。参数说明limit3表示最多连续插补 3 年超过 3 年的缺测不插补因为线性插值在长缺测段会引入虚假趋势。limit_directionboth表示序列首尾的缺测也尝试插补。如果缺测率超过 10%建议换站或使用更复杂的时空插补方法。3.2 Pettitt 检验识别突变点Pettitt 检验是一种非参数突变点检测方法不要求数据服从特定分布适合最大风速这种偏态序列。def pettitt_test(x, alpha0.05): Pettitt 突变点检验 x: 输入序列一维数组 alpha: 显著性水平 返回: (突变点位置索引, p值, 是否显著) n len(x) U np.zeros(n) for t in range(1, n): U[t] U[t-1] np.sign(x[t] - np.median(x[:t])) * (t) # 简化版 # 标准 Pettitt 统计量 U_full np.zeros(n) for t in range(n): for i in range(t1): for j in range(t1, n): U_full[t] np.sign(x[i] - x[j]) K np.max(np.abs(U_full)) # 近似 p 值 p_value 2 * np.exp(-6 * K**2 / (n**3 n**2)) # 突变点位置 change_point np.argmax(np.abs(U_full)) return change_point, p_value, p_value alpha # 对最大风速序列做 Pettitt 检验 x df_target[max_wind_speed_filled].dropna().values cp, p, sig pettitt_test(x) print(f突变点年份: {df_target[year].iloc[cp]}, p值: {p:.4f}, 显著: {sig})Pettitt 检验的核心思想是如果序列存在突变点那么突变点前后的子序列之间的秩和差异会很大。U_full[t]统计的是前 t1 个点和后 n-t-1 个点之间的符号差异累积。K取最大绝对值p_value用近似公式计算。参数alpha0.05是显著性水平p 值小于 0.05 认为突变显著。注意Pettitt 检验只能检出一个突变点如果序列有多个断点需要分段检验或改用其他方法。3.3 SNHT 检验与参考站对比SNHT 需要参考站序列。假设已经选好参考站数据在同一 CSV 中。def snht_test(target, reference, alpha0.05): SNHT 均一性检验 target: 待检序列 reference: 参考序列 返回: (最大T值位置, p值, 是否显著) n len(target) # 计算比值序列或差值序列取决于变量类型 # 最大风速用差值更合适因为比值在风速接近0时不稳定 diff target - reference # 标准化 diff_std (diff - np.mean(diff)) / np.std(diff) T np.zeros(n) for t in range(1, n): # 前段均值和后段均值 mean1 np.mean(diff_std[:t]) mean2 np.mean(diff_std[t:]) T[t] t * mean1**2 (n - t) * mean2**2 T_max np.max(T) # 临界值近似n20时 # 更精确的临界值需要查表或蒙特卡洛模拟 T_critical 8.0 # 对应 alpha0.05 的近似值 change_point np.argmax(T) return change_point, T_max, T_max T_critical # 假设参考站数据在同一 DataFrame 中 ref_station 54512 df_ref df[df[station_id] ref_station].sort_values(year) # 对齐年份 merged df_target.merge(df_ref[[year, max_wind_speed]], onyear, suffixes(_target, _ref)) merged merged.dropna(subset[max_wind_speed_target, max_wind_speed_ref]) cp_snht, T_max, sig_snht snht_test( merged[max_wind_speed_target].values, merged[max_wind_speed_ref].values ) print(fSNHT 突变点年份: {merged[year].iloc[cp_snht]}, T{T_max:.2f}, 显著: {sig_snht})SNHT 的核心是构造一个统计量 T衡量突变点前后两段均值的差异。diff是待检站和参考站的差值序列标准化后消除量纲影响。T[t]在突变点处达到最大。临界值T_critical8.0是 n20 时的近似值更精确的做法是用蒙特卡洛模拟生成临界值表。注意SNHT 对序列两端的突变点检测能力较弱如果突变发生在序列开头或结尾 5 年内结果要谨慎。3.4 分位数映射校正识别出突变点后用分位数映射做校正。假设突变点在 1995 年1995 年之前是“旧”时段之后是“新”时段需要把旧时段校正到新时段的分布。def quantile_mapping(data_old, data_new, data_to_correct): 分位数映射校正 data_old: 旧时段数据参考分布 data_new: 新时段数据目标分布 data_to_correct: 需要校正的数据 返回: 校正后的数据 # 经验累积分布 old_sorted np.sort(data_old) new_sorted np.sort(data_new) # 对每个待校正值找到它在旧分布中的分位数 # 然后映射到新分布的对应分位数 corrected np.zeros_like(data_to_correct) for i, val in enumerate(data_to_correct): # 在旧分布中的分位数 p np.searchsorted(old_sorted, val) / len(old_sorted) p np.clip(p, 0.01, 0.99) # 避免极端分位数 # 在新分布中的对应值 idx int(p * len(new_sorted)) idx min(idx, len(new_sorted) - 1) corrected[i] new_sorted[idx] return corrected # 分段 break_year 1995 old_data merged[merged[year] break_year][max_wind_speed_target].values new_data merged[merged[year] break_year][max_wind_speed_target].values # 校正旧时段数据 corrected_old quantile_mapping(old_data, new_data, old_data) # 合并校正后的序列 corrected_series np.concatenate([corrected_old, new_data]) years np.concatenate([ merged[merged[year] break_year][year].values, merged[merged[year] break_year][year].values ]) # 可视化对比 fig, axes plt.subplots(1, 2, figsize(12, 4)) axes[0].plot(years, np.concatenate([old_data, new_data]), b-, label原始) axes[0].plot(years, corrected_series, r-, label校正后) axes[0].axvline(break_year, colorgray, linestyle--, label突变点) axes[0].set_xlabel(年份) axes[0].set_ylabel(最大风速 (m/s)) axes[0].legend() axes[0].set_title(序列对比) axes[1].hist(old_data, bins10, alpha0.5, label旧时段, densityTrue) axes[1].hist(new_data, bins10, alpha0.5, label新时段, densityTrue) axes[1].hist(corrected_old, bins10, alpha0.5, label校正后旧时段, densityTrue) axes[1].set_xlabel(最大风速 (m/s)) axes[1].legend() axes[1].set_title(分布对比) plt.tight_layout() plt.savefig(homogenization_result.png, dpi150) plt.show()分位数映射的逻辑是对旧时段的每个值找到它在旧分布中的分位数然后取新分布中同一分位数的值作为校正值。np.searchsorted用于快速定位分位数np.clip把分位数限制在 0.01-0.99 之间避免极端分位数导致校正值失真。参数说明break_year是突变点年份需要根据前面的检验结果设定。校正后的序列在突变点前后分布一致可以用于后续趋势分析。4. 避坑指南最大风速均一化订正中最容易翻车的 5 个地方4.1 坑一参考站选得太远订正后引入新偏差现象订正后的序列趋势和参考站趋势高度一致但和待订正站周边的其他站趋势不一致。原因参考站距离太远最大风速的空间相关性已经很低参考站的信息主要反映的是它自己的局地气候不是待订正站的气候信号。强行用远距离参考站做 SNHT会把参考站的局地特征“传染”给待订正站。解决参考站距离严格控制在 50 km 以内山区控制在 30 km 以内。如果找不到宁可做单站订正。单站订正虽然可靠性下降但不会引入虚假的空间信号。可以用多个参考站做交叉验证如果不同参考站给出的断点位置差异很大说明参考站信息不可靠。4.2 坑二忽略元数据统计断点和实际断点对不上现象Pettitt 检验和 SNHT 检验都显示 1995 年有突变但台站元数据里 1995 年没有任何仪器更换或迁站记录。反而 2003 年有明确的仪器更换记录但统计检验在 2003 年没有检出显著突变。原因统计检验检出的是“数据分布变化”不一定是“仪器变化”。1995 年的突变可能来自观测规范变更比如最大风速的统计时段从 10 分钟改为 2 分钟这种变更在元数据里可能没有详细记录。而 2003 年的仪器更换可能恰好没有改变最大风速的统计特性比如两种仪器在最大风速量级上响应一致。解决统计检验和元数据必须交叉验证。元数据是“金标准”统计检验是“辅助工具”。如果统计检验检出的断点在元数据中有对应记录优先采信。如果统计检验检出但元数据没有记录需要进一步排查——可能是元数据不完整也可能是统计检验的假阳性。可以用贝叶斯方法给出断点的概率分布而不是硬性判定“有”或“没有”。4.3 坑三分位数映射在样本量小时外推过度现象校正后的旧时段最大风速出现了不合理的极值比如校正后出现了 50 m/s 的风速但原始序列最大值只有 35 m/s。原因经验分位数映射在样本量小时极端分位数比如 0.95 以上的估计非常不稳定。旧时段可能只有 20 个点0.95 分位数对应的是第 19 个点稍微换个样本这个值就变了。如果新时段的 0.95 分位数恰好很大校正后就会产生虚假极值。解决样本量小于 30 时优先用参数化分位数映射。先拟合 GEV 分布再用拟合的分布做映射。GEV 分布的尾部有参数控制不会像经验分位数那样剧烈波动。如果一定要用经验分位数映射把分位数范围限制在 0.05-0.95 之间超出范围的值用参数化方法外推。4.4 坑四多断点序列只做单断点订正现象序列在 1985 年和 2005 年各有一个断点但只做了 1985 年的订正2005 年之后的序列仍然有偏差。原因Pettitt 检验和 SNHT 检验默认只检出一个断点。如果序列有多个断点单断点方法只能找到最显著的那个其他断点被忽略。解决先用单断点方法找到最显著的断点把序列分成两段然后在每段内再做单断点检验。重复这个过程直到没有显著断点。这就是“分段检验”的思路。更严谨的做法是用多断点方法比如贝叶斯多断点模型或动态规划方法。但多断点方法计算量大而且断点越多不确定性越大。实际业务中如果断点超过 3 个建议直接剔除该站因为订正后的序列可靠性已经很低了。4.5 坑五订正后不做独立验证现象订正后的序列趋势合理但用来做极值重现期推算时得到的 50 年一遇风速和周边站差异很大。原因订正只保证了序列内部的均一性没有保证序列和周边站的时空一致性。如果订正过程中引入了偏差序列内部看起来均一但和周边站对比就会暴露问题。解决订正后必须做独立验证。常用方法有① 和周边未订正站做空间一致性检验订正后的序列应该和周边站的相关性更高② 用订正后的序列做极值推算和用原始序列、周边站序列的结果对比差异应该在合理范围内③ 如果有可能用独立时段的数据做交叉验证——比如用 1980-2000 年数据建模用 2001-2020 年数据验证。5. 订正后的序列怎么用趋势检验与极值推算的衔接订正完序列下一步通常是做趋势分析或极值推算。这两件事对订正质量的要求不同衔接时需要注意几个技巧。趋势分析对订正质量最敏感。如果订正不彻底残留的断点会直接污染趋势估计。我一般会做“双保险”先用 Mann-Kendall 检验做非参数趋势检验再用 Sens slope 估计趋势幅度。Mann-Kendall 检验不要求正态分布对最大风速这种偏态序列比较稳健。Sens slope 比线性回归斜率更抗 outliers。如果两种方法给出的趋势方向一致基本可信。如果不一致说明序列里还有未识别的断点或 outliers。import pymannkendall as mk # 对订正后的序列做 Mann-Kendall 趋势检验 result mk.original_test(corrected_series) print(f趋势: {result.trend}, p值: {result.p:.4f}, Sens slope: {result.slope:.4f}) # 如果 p 0.05趋势显著 # Sens slope 的单位是 m/s per year极值推算对订正质量的要求更高因为极值推算依赖分布的尾部。订正后的序列如果尾部被扭曲重现期估计会严重偏差。我一般会先用订正后的序列拟合 GEV 分布然后用轮廓似然法估计重现期。轮廓似然法比矩估计法更稳健尤其在样本量小时。from scipy.stats import genextreme as gev # 拟合 GEV 分布 params gev.fit(corrected_series) shape, loc, scale params # 推算 50 年一遇最大风速 return_period 50 p 1 - 1/return_period wind_50 gev.ppf(p, shape, loc, scale) print(f50年一遇最大风速: {wind_50:.2f} m/s) # 轮廓似然法估计置信区间简化版 # 实际应用中建议用专门的极值分析包如 pyextremes这里有个血泪经验订正后的序列在做极值推算时不要直接用原始序列的极值。订正可能会改变极值的大小尤其是分位数映射校正后旧时段的极值可能被“拉”到新时段的分布上。如果旧时段的极值被拉得过高重现期估计会偏大。我一般会对比订正前后极值推算结果如果差异超过 20%就需要回头检查订正过程。最后一个技巧订正后的序列建议保留“订正标记”。在数据表里加一列is_corrected标记哪些年份的数据被校正过。这样后续分析时可以区分“原始观测”和“校正值”在做敏感性分析时能快速切换。这个习惯帮我省了很多后悔药——有一次合作方质疑订正结果我直接拉出标记列对比了订正前后的趋势问题当场定位。希望帮到你。本文还有配套的精品资源点击获取