简介面向气象水文领域的研究人员、学生及相关数据分析者这份以最大风速为例的均一化订正资源系统演示了如何消除因仪器更换、站点迁移或测量方法改变导致的系统性偏差使不同时期和不同站点的风速记录具备可比性从而为气候诊断、天气预报和环境评估提供更可靠的基础数据。压缩包共包含三个文件内有可直接读取的站点最大风速原始数据、实现完整订正流程的程序脚本以及配套的方法说明文档其中数据文件覆盖一九七一至二〇一五年程序脚本内含数据清洗、缺失值处理、订正算法、结果验证与可视化等模块整包约三百七十七KB轻量便捷目前已有五百九十人学习。通过运行该程序读者既能快速复现最大风速均一化订正的全过程也能参照说明文档理解每个步骤的统计原理和参数选择进而将同样思路迁移到气温、降水等其他气象要素的订正处理中切实提升气象数据质量控制与研究分析能力。1. 最大风速序列里的“假变异”均一化订正到底在修正什么如果手上有一份1971—2015年的最大风速序列画出来发现1989年前后突然掉了1.2 m/s第一反应往往是“气候变了”。但翻一翻观测簿那一年换了测风仪站址还向东挪了200米。仪器更换、站址迁移、周边建筑变化都会让风速序列产生系统性跳变这就是气象水文数据处理里常说的非均一性问题。如果不处理这种“假变异”会被误当成真实气候信号后续的趋势分析、重现期极值风速计算都会跟着偏。均一化订正要做的事就是结合统计检验与元数据记录识别断点并施加订正因子让不同时段的记录回到同一观测基准上。这篇博客以最大风速为例用压缩包里的 YJ_vmax_1971-2015.csv 和 homogenization.py 走通读取、断点检测、订正、验证的完整链路。从事气候数据、水资源分析、风荷载评估的同行都可以直接照着改。2. CSV读入与异常值清洗先把序列变成可检验的对象2.1 先判断你的表格是年度极值还是逐月极值打开 CSV 后第一件事不是跑代码而是确认数据粒度。最常见的两种结构年份最大风速或者年份月份最大风速。这直接决定后续滑动窗口参数的取值。年度极值一年一条窗口取11年会吞掉太多样本逐月数据则可以用24个月窗口做参考序列。如果表格是从 Excel 导出的还要留意编码问题。年份列正常风速列是数值但月份列可能带着“月”字后缀读进来变成字符串排序和索引都会出错。import pandas as pd import numpy as np df pd.read_csv(YJ_vmax_1971-2015.csv, encodingutf-8) print(df.head()) print(df.info())如果打印时提示 UnicodeDecodeError把 encoding 改成 gbk 重试。这是 Windows 下 Excel 导出 CSV 最常见的编码坑没有之一。列名确认没问题后再做时间索引if month in df.columns: df[dt] pd.to_datetime( dict(yeardf[year], monthdf[month], day1) ) df df.set_index(dt).sort_index() else: df[year] df[year].astype(int) df df.set_index( pd.to_datetime(df[year], format%Y) ).sort_index()有 month 字段就按月度索引否则按年度索引。统一成时间索引后滚动平均、断点检测、绘图可以直接复用 pandas 自带的时间序列接口。把这份代码存成 homogenization.py在命令行执行python homogenization.py就能从文本文档变成可运行脚本不用依赖 IDE。2.2 缺失值与粗差什么该插补什么该剔除风速记录里的缺失值不算多但处理策略要提前定。年度最大风速缺失一年那一年就不能参与极值分析月度缺失可以用邻月插值补前提是连续缺失不超过6个月。超过这个限度插值出来的“最大风速”已经没有任何物理意义不如直接标成 NaN 并在后续分析里跳过。处理对象常用判据处理动作说明缺失值连续缺失不超过6个月线性插值超限置 NaN 并打标记粗差abs(x - mean) 3 * std剔除后插值只处理孤立异常点系统性偏差SNHT 断点检测均一化订正清洗解决不了交给第3章miss df[vmax].isna().sum() print(f缺失记录数: {miss}) m, s df[vmax].mean(), df[vmax].std() out_mask (df[vmax] m - 3 * s) | (df[vmax] m 3 * s) df.loc[out_mask, vmax] np.nan print(f标记粗差: {out_mask.sum()} 条) df[vmax] df[vmax].interpolate(limit6)3-sigma 判据对近似正态的变量有效但最大风速更接近 Gumbel 分布右尾偏长。如果序列偏态明显我会先对风速取对数再做 3-sigma 检查避免把真正的台风过程误判成粗差。插值用线性即可风速年际自相关弱高阶插值反而引入虚假波动。2.3 参考序列怎么构造单站数据也有办法SNHT 这类均一化方法通常需要参考序列。有邻近站时优先选相关系数大于 0.6、且本身已经过均一化订正的站点按相关系数加权合成参考序列。没有邻近站时常见做法是对目标站自身做中心滑动平均生成“拟参考序列”。ref df[vmax].rolling( window11, centerTrue, min_periods5 ).mean()窗口取 11 是经验值。序列长度只有 20 年时窗口降到 7否则边缘年份会产生太多 NaN。单站拟参考序列存在自相关问题断点检测结果只能作为初筛不能当作确证最终结论要靠与元数据互证。数据不一致在这里很常见——统计上显著但观测簿上那年既没换仪器也没迁站我倾向于保留原始记录在说明文件里备注“疑似非均一”而不是强行订正。3. SNHT断点检测与订正因子计算锁定跳变年份3.1 为什么要检测目标站与参考站的差异序列直接对原始风速做均值变化检验会把逐年自然波动当作断点。SNHT 的核心思路是构造目标站与参考站的差异序列再检验差异序列是否存在均值突变。两个站经历相似的大尺度天气过程差异序列里的自然变率被抵消掉大部分剩下的主要是站点自身的系统性变化。对最大风速取对数后再做差相当于把乘性关系变成加性关系真实风速乘上站点影响因子取对数后站点因子变成加性常数均值突变检验就能识别出这个常数的改变。这个变换是 SNHT 区别于普通 t 检验的关键。3.2 SNHT 统计量的计算逻辑def snht_test(v, ref, min_seg10): 返回候选断点索引与T统计量 diff np.log(np.maximum(v, 0.1)) - np.log( np.maximum(ref, 0.1) ) z (diff - np.nanmean(diff)) / np.nanstd(diff) n len(z) best_t, best_k -1.0, 0 for i in range(min_seg, n - min_seg): z1, z2 np.nanmean(z[:i]), np.nanmean(z[i:]) t i * z1 ** 2 (n - i) * z2 ** 2 if t best_t: best_t, best_k t, i return best_k, best_t对每个可能断点位置 i计算前后两段标准化均值 z1、z2构造统计量 T i * z1^2 (n - i) * z2^2遍历后取最大值。T 越大说明 i 点前后均值差异越显著。min_seg10 保证断点两边至少各占 10 个样本否则订正因子的估计误差会大得没法用。T 的临界值依赖序列长度和显著性水平n45 时 T 超过 9.2 可作 95% 显著的经验阈值。要精确值可以查 Alexandersson 的临界值表我在代码里一般先用 9.3 做初筛再把可疑年份和观测簿逐条对。3.3 加性订正还是乘性订正温度偏差与绝对值关系不大用加性订正减去一个常数即可。最大风速的量级与测站地形、仪器性能直接挂钩换仪器造成的通常是整体缩放乘性订正更合理。降水这类含大量零值的变量先做 log1p 变换再乘性订正最大风速很少为 0直接乘因子即可。订正方式公式适用变量风险加性x x delta温度、气压对量级敏感的序列残留比例误差乘性x x * factor风速、降水、辐射因子依赖分段均值样本少时偏大加性订正和乘性订正的选择不是拍脑袋。如果把最大风速从 25 m/s 和 15 m/s 两段直接加一个常数低风速段的相对偏差会被放大乘性因子同时压缩整段分布更贴合测风仪灵敏度变化带来的系统效应。4. 把订正流程写进 homogenization.py迭代检测与序列重建4.1 多断点迭代订正一次还不够序列可能有两个以上断点换过两次仪器就至少两个。如果检测到一个断点就订正并退出第二次跳变会被第一次订正扭曲。常见做法是循环检测订正后重新运行 SNHT直到找不出显著断点才算结束。def detect_breakpoints(vmax, ref_builder, t_thresh9.3, min_seg10, max_iter5): 迭代检测断点返回断点列表与订正后序列 series vmax.copy() breaks [] for _ in range(max_iter): ref ref_builder(series) k, t snht_test(series, ref, min_segmin_seg) if t t_thresh: break breaks.append((series.index[k], t)) before series.iloc[:k].mean() after series.iloc[k:].mean() factor after / before if before 0 else 1.0 series.iloc[:k] series.iloc[:k] * factor return breaks, series这里有个关键选型基准段用了断点之后的“新时段”。假设1995年换仪器且观测簿记录新仪器更可靠就把1971—1994年的风速统一乘以“1995—2015均值 / 1971—1994均值”。反过来把后半段订正到前半段也能成立但序列会带着老仪器的偏差走完整条趋势线论文里容易被审稿人追问。所以我的默认选择是以后段为基准。max_iter5 是保险丝防止异常序列让循环跑不完。4.2 订正因子估计的可靠性两段各20年的因子精度比两段各10年的因子高很多。因子标准误差可以用 bootstrap 估计从前后两段分别有放回抽样各抽500次计算 factor 分布取 2.5% 到 97.5% 分位作为置信区间。如果区间宽度超过 0.1说明因子在 0.90 到 1.10 之间还晃来晃去断点的证据不足。这时候我倾向降低 T 阈值重新检测或者把该断点标记为“弱断点”在数据说明里注明。很多人只看 T 统计量过没过线就订正结果把 ENSO 引起的年代际波动也当成换仪器造成的偏差校正掉了。订正因子估计的不确定性直接决定了订正后序列的可信度这个步骤不能省。4.3 输出订正后的序列与订正记录corrected, breaks detect_breakpoints( df[vmax], ref_builderlambda s: s.rolling( 11, centerTrue, min_periods5 ).mean(), ) df[vmax_homogenized] corrected print(断点记录:, breaks) df.to_csv(YJ_vmax_homogenized.csv, encodingutf-8-sig)输出时保留原始列和订正列断点年份、T 值、订正因子单独存一个说明文件。千万别只覆盖原列——后续复现图表时原始序列和订正序列都必须保留。编码用 utf-8-sig 是为了 Excel 打开时不乱码这也是交付数据文件时的通用做法。5. 订正效果验证别只看图要过统计检验5.1 断点前后差异的显著性复检订正是基于 SNHT 检测结果做的但订正完成后还要独立验证订正后的断点位置前后均值差异应该不再显著。from scipy import stats bp breaks[0][0] if breaks else None if bp is not None: before df.loc[ df.index bp, vmax_homogenized ] after df.loc[ df.index bp, vmax_homogenized ] t_stat, p_val stats.ttest_ind(before, after) print(f断点 {bp.date()} 前后 t 检验: p{p_val:.3f})p 值大于 0.05 表示订正后断点两侧均值差异不显著。45 年序列分两段每段约 20 年t 检验效能足够。再看 Levene 方差检验换仪器有时只改变波动幅度、不改变均值方差检验能补上 t 检验看不到的问题。5.2 分布形态对比QQ 图与直方图订正后的序列应该更接近理论极值分布。最大风速一般用 Gumbel 分布拟合如果原始序列混入两套不同水平的记录QQ 图尾部会出现明显弯折订正后弯折减轻说明分布一致性提升。import matplotlib.pyplot as plt fig, axes plt.subplots(1, 2, figsize(10, 4)) for ax, col, title in zip( axes, [vmax, vmax_homogenized], [原始序列, 订正序列], ): stats.probplot( df[col].dropna(), diststats.gumbel_r, plotax, ) ax.set_title(title) plt.tight_layout() plt.show()如果 QQ 图尾部弯折还在多半是断点没检全。回头把 max_iter 调大或者把 min_seg 降到 8 再跑一轮。还有一种情况断点位置本身没找对订正因子施加到错误分段上QQ 图弯折会更严重这时候重新核对 SNHT 返回的索引和观测簿记录。5.3 CUSUM 图核查CUSUM 是把序列减去均值后累加断点在图上表现为斜率突变。SNHT 本质上就是 CUSUM 思想的统计化版本看 CUSUM 图能直观确认断点位置是否合理。diff df[vmax_homogenized] - df[vmax_homogenized].mean() cusum diff.cumsum() plt.plot(df.index, cusum, label订正后) plt.axvline( pd.Timestamp(breaks[0][0].date()), colorgray, linestyle--, label断点 ) plt.xlabel(年份) plt.ylabel(累积距平) plt.legend() plt.show()CUSUM 对线性趋势也敏感。如果序列本身存在真实气候趋势累积和会持续上升不要一看到单调上升就认定是断点。判断断点要看斜率变化点而不是累积量本身的大小。提示SNHT 检测结果只能说明“这里存在统计显著跳变”不能证明“一定是非均一性”。先查元数据再下结论。6. 从最大风速延伸到降水、温度参数调整的四个关键位置同一条流程换到其他气象要素只需要调整四个位置其余代码可以原样复用。第一处是变换函数。最大风速直接取对数做差降水含大量零值先 log1p 再进入 SNHT避免零值取对数产生负无穷温度用原始值做差订正时也改为加性因子不做乘法。第二处是参考序列的平滑窗口。年度数据用 11 年窗口对 45 年序列够用月度数据窗口要短我用 25 个月既能抑制季节波动又不会把年际信号抹平。第三处是 SNHT 阈值。阈值本质上是显著性水平的映射序列越短 T 的自然波动越大。我习惯用蒙特卡洛模拟生成白噪声序列跑 1000 次 SNHT取 95% 分位作当前序列长度的阈值不做硬编码。第四处是订正基准段的选择。有元数据记录时以仪器更换后的段为基准没有元数据时以最近一个完整 10 年段为基准并在输出文件里注明“基准段依据统计推断待元数据复核”。还有一个容易踩的坑不要过度订正。SNHT 检测出断点时先对照观测簿确认那一年有没有换仪、迁站记录。完全对不上且断点前后均值差小于 0.5 m/s 时我一般标记为“疑似非均一”而不是直接订正这样后续做重现期分析时还有选择余地。订正完成后把断点年份、T 值、订正因子、是否与元数据吻合整理成一张表连同订正前后序列一并写入 CSV。这张表就是论文方法部分最有力的附件审稿人追问“怎么确定断点真实存在”时直接对表逐条说明即可。本文还有配套的精品资源点击获取 最大风速序列均一化订正：SNHT断点检测与Python实现
简介面向气象水文领域的研究人员、学生及相关数据分析者这份以最大风速为例的均一化订正资源系统演示了如何消除因仪器更换、站点迁移或测量方法改变导致的系统性偏差使不同时期和不同站点的风速记录具备可比性从而为气候诊断、天气预报和环境评估提供更可靠的基础数据。压缩包共包含三个文件内有可直接读取的站点最大风速原始数据、实现完整订正流程的程序脚本以及配套的方法说明文档其中数据文件覆盖一九七一至二〇一五年程序脚本内含数据清洗、缺失值处理、订正算法、结果验证与可视化等模块整包约三百七十七KB轻量便捷目前已有五百九十人学习。通过运行该程序读者既能快速复现最大风速均一化订正的全过程也能参照说明文档理解每个步骤的统计原理和参数选择进而将同样思路迁移到气温、降水等其他气象要素的订正处理中切实提升气象数据质量控制与研究分析能力。1. 最大风速序列里的“假变异”均一化订正到底在修正什么如果手上有一份1971—2015年的最大风速序列画出来发现1989年前后突然掉了1.2 m/s第一反应往往是“气候变了”。但翻一翻观测簿那一年换了测风仪站址还向东挪了200米。仪器更换、站址迁移、周边建筑变化都会让风速序列产生系统性跳变这就是气象水文数据处理里常说的非均一性问题。如果不处理这种“假变异”会被误当成真实气候信号后续的趋势分析、重现期极值风速计算都会跟着偏。均一化订正要做的事就是结合统计检验与元数据记录识别断点并施加订正因子让不同时段的记录回到同一观测基准上。这篇博客以最大风速为例用压缩包里的 YJ_vmax_1971-2015.csv 和 homogenization.py 走通读取、断点检测、订正、验证的完整链路。从事气候数据、水资源分析、风荷载评估的同行都可以直接照着改。2. CSV读入与异常值清洗先把序列变成可检验的对象2.1 先判断你的表格是年度极值还是逐月极值打开 CSV 后第一件事不是跑代码而是确认数据粒度。最常见的两种结构年份最大风速或者年份月份最大风速。这直接决定后续滑动窗口参数的取值。年度极值一年一条窗口取11年会吞掉太多样本逐月数据则可以用24个月窗口做参考序列。如果表格是从 Excel 导出的还要留意编码问题。年份列正常风速列是数值但月份列可能带着“月”字后缀读进来变成字符串排序和索引都会出错。import pandas as pd import numpy as np df pd.read_csv(YJ_vmax_1971-2015.csv, encodingutf-8) print(df.head()) print(df.info())如果打印时提示 UnicodeDecodeError把 encoding 改成 gbk 重试。这是 Windows 下 Excel 导出 CSV 最常见的编码坑没有之一。列名确认没问题后再做时间索引if month in df.columns: df[dt] pd.to_datetime( dict(yeardf[year], monthdf[month], day1) ) df df.set_index(dt).sort_index() else: df[year] df[year].astype(int) df df.set_index( pd.to_datetime(df[year], format%Y) ).sort_index()有 month 字段就按月度索引否则按年度索引。统一成时间索引后滚动平均、断点检测、绘图可以直接复用 pandas 自带的时间序列接口。把这份代码存成 homogenization.py在命令行执行python homogenization.py就能从文本文档变成可运行脚本不用依赖 IDE。2.2 缺失值与粗差什么该插补什么该剔除风速记录里的缺失值不算多但处理策略要提前定。年度最大风速缺失一年那一年就不能参与极值分析月度缺失可以用邻月插值补前提是连续缺失不超过6个月。超过这个限度插值出来的“最大风速”已经没有任何物理意义不如直接标成 NaN 并在后续分析里跳过。处理对象常用判据处理动作说明缺失值连续缺失不超过6个月线性插值超限置 NaN 并打标记粗差abs(x - mean) 3 * std剔除后插值只处理孤立异常点系统性偏差SNHT 断点检测均一化订正清洗解决不了交给第3章miss df[vmax].isna().sum() print(f缺失记录数: {miss}) m, s df[vmax].mean(), df[vmax].std() out_mask (df[vmax] m - 3 * s) | (df[vmax] m 3 * s) df.loc[out_mask, vmax] np.nan print(f标记粗差: {out_mask.sum()} 条) df[vmax] df[vmax].interpolate(limit6)3-sigma 判据对近似正态的变量有效但最大风速更接近 Gumbel 分布右尾偏长。如果序列偏态明显我会先对风速取对数再做 3-sigma 检查避免把真正的台风过程误判成粗差。插值用线性即可风速年际自相关弱高阶插值反而引入虚假波动。2.3 参考序列怎么构造单站数据也有办法SNHT 这类均一化方法通常需要参考序列。有邻近站时优先选相关系数大于 0.6、且本身已经过均一化订正的站点按相关系数加权合成参考序列。没有邻近站时常见做法是对目标站自身做中心滑动平均生成“拟参考序列”。ref df[vmax].rolling( window11, centerTrue, min_periods5 ).mean()窗口取 11 是经验值。序列长度只有 20 年时窗口降到 7否则边缘年份会产生太多 NaN。单站拟参考序列存在自相关问题断点检测结果只能作为初筛不能当作确证最终结论要靠与元数据互证。数据不一致在这里很常见——统计上显著但观测簿上那年既没换仪器也没迁站我倾向于保留原始记录在说明文件里备注“疑似非均一”而不是强行订正。3. SNHT断点检测与订正因子计算锁定跳变年份3.1 为什么要检测目标站与参考站的差异序列直接对原始风速做均值变化检验会把逐年自然波动当作断点。SNHT 的核心思路是构造目标站与参考站的差异序列再检验差异序列是否存在均值突变。两个站经历相似的大尺度天气过程差异序列里的自然变率被抵消掉大部分剩下的主要是站点自身的系统性变化。对最大风速取对数后再做差相当于把乘性关系变成加性关系真实风速乘上站点影响因子取对数后站点因子变成加性常数均值突变检验就能识别出这个常数的改变。这个变换是 SNHT 区别于普通 t 检验的关键。3.2 SNHT 统计量的计算逻辑def snht_test(v, ref, min_seg10): 返回候选断点索引与T统计量 diff np.log(np.maximum(v, 0.1)) - np.log( np.maximum(ref, 0.1) ) z (diff - np.nanmean(diff)) / np.nanstd(diff) n len(z) best_t, best_k -1.0, 0 for i in range(min_seg, n - min_seg): z1, z2 np.nanmean(z[:i]), np.nanmean(z[i:]) t i * z1 ** 2 (n - i) * z2 ** 2 if t best_t: best_t, best_k t, i return best_k, best_t对每个可能断点位置 i计算前后两段标准化均值 z1、z2构造统计量 T i * z1^2 (n - i) * z2^2遍历后取最大值。T 越大说明 i 点前后均值差异越显著。min_seg10 保证断点两边至少各占 10 个样本否则订正因子的估计误差会大得没法用。T 的临界值依赖序列长度和显著性水平n45 时 T 超过 9.2 可作 95% 显著的经验阈值。要精确值可以查 Alexandersson 的临界值表我在代码里一般先用 9.3 做初筛再把可疑年份和观测簿逐条对。3.3 加性订正还是乘性订正温度偏差与绝对值关系不大用加性订正减去一个常数即可。最大风速的量级与测站地形、仪器性能直接挂钩换仪器造成的通常是整体缩放乘性订正更合理。降水这类含大量零值的变量先做 log1p 变换再乘性订正最大风速很少为 0直接乘因子即可。订正方式公式适用变量风险加性x x delta温度、气压对量级敏感的序列残留比例误差乘性x x * factor风速、降水、辐射因子依赖分段均值样本少时偏大加性订正和乘性订正的选择不是拍脑袋。如果把最大风速从 25 m/s 和 15 m/s 两段直接加一个常数低风速段的相对偏差会被放大乘性因子同时压缩整段分布更贴合测风仪灵敏度变化带来的系统效应。4. 把订正流程写进 homogenization.py迭代检测与序列重建4.1 多断点迭代订正一次还不够序列可能有两个以上断点换过两次仪器就至少两个。如果检测到一个断点就订正并退出第二次跳变会被第一次订正扭曲。常见做法是循环检测订正后重新运行 SNHT直到找不出显著断点才算结束。def detect_breakpoints(vmax, ref_builder, t_thresh9.3, min_seg10, max_iter5): 迭代检测断点返回断点列表与订正后序列 series vmax.copy() breaks [] for _ in range(max_iter): ref ref_builder(series) k, t snht_test(series, ref, min_segmin_seg) if t t_thresh: break breaks.append((series.index[k], t)) before series.iloc[:k].mean() after series.iloc[k:].mean() factor after / before if before 0 else 1.0 series.iloc[:k] series.iloc[:k] * factor return breaks, series这里有个关键选型基准段用了断点之后的“新时段”。假设1995年换仪器且观测簿记录新仪器更可靠就把1971—1994年的风速统一乘以“1995—2015均值 / 1971—1994均值”。反过来把后半段订正到前半段也能成立但序列会带着老仪器的偏差走完整条趋势线论文里容易被审稿人追问。所以我的默认选择是以后段为基准。max_iter5 是保险丝防止异常序列让循环跑不完。4.2 订正因子估计的可靠性两段各20年的因子精度比两段各10年的因子高很多。因子标准误差可以用 bootstrap 估计从前后两段分别有放回抽样各抽500次计算 factor 分布取 2.5% 到 97.5% 分位作为置信区间。如果区间宽度超过 0.1说明因子在 0.90 到 1.10 之间还晃来晃去断点的证据不足。这时候我倾向降低 T 阈值重新检测或者把该断点标记为“弱断点”在数据说明里注明。很多人只看 T 统计量过没过线就订正结果把 ENSO 引起的年代际波动也当成换仪器造成的偏差校正掉了。订正因子估计的不确定性直接决定了订正后序列的可信度这个步骤不能省。4.3 输出订正后的序列与订正记录corrected, breaks detect_breakpoints( df[vmax], ref_builderlambda s: s.rolling( 11, centerTrue, min_periods5 ).mean(), ) df[vmax_homogenized] corrected print(断点记录:, breaks) df.to_csv(YJ_vmax_homogenized.csv, encodingutf-8-sig)输出时保留原始列和订正列断点年份、T 值、订正因子单独存一个说明文件。千万别只覆盖原列——后续复现图表时原始序列和订正序列都必须保留。编码用 utf-8-sig 是为了 Excel 打开时不乱码这也是交付数据文件时的通用做法。5. 订正效果验证别只看图要过统计检验5.1 断点前后差异的显著性复检订正是基于 SNHT 检测结果做的但订正完成后还要独立验证订正后的断点位置前后均值差异应该不再显著。from scipy import stats bp breaks[0][0] if breaks else None if bp is not None: before df.loc[ df.index bp, vmax_homogenized ] after df.loc[ df.index bp, vmax_homogenized ] t_stat, p_val stats.ttest_ind(before, after) print(f断点 {bp.date()} 前后 t 检验: p{p_val:.3f})p 值大于 0.05 表示订正后断点两侧均值差异不显著。45 年序列分两段每段约 20 年t 检验效能足够。再看 Levene 方差检验换仪器有时只改变波动幅度、不改变均值方差检验能补上 t 检验看不到的问题。5.2 分布形态对比QQ 图与直方图订正后的序列应该更接近理论极值分布。最大风速一般用 Gumbel 分布拟合如果原始序列混入两套不同水平的记录QQ 图尾部会出现明显弯折订正后弯折减轻说明分布一致性提升。import matplotlib.pyplot as plt fig, axes plt.subplots(1, 2, figsize(10, 4)) for ax, col, title in zip( axes, [vmax, vmax_homogenized], [原始序列, 订正序列], ): stats.probplot( df[col].dropna(), diststats.gumbel_r, plotax, ) ax.set_title(title) plt.tight_layout() plt.show()如果 QQ 图尾部弯折还在多半是断点没检全。回头把 max_iter 调大或者把 min_seg 降到 8 再跑一轮。还有一种情况断点位置本身没找对订正因子施加到错误分段上QQ 图弯折会更严重这时候重新核对 SNHT 返回的索引和观测簿记录。5.3 CUSUM 图核查CUSUM 是把序列减去均值后累加断点在图上表现为斜率突变。SNHT 本质上就是 CUSUM 思想的统计化版本看 CUSUM 图能直观确认断点位置是否合理。diff df[vmax_homogenized] - df[vmax_homogenized].mean() cusum diff.cumsum() plt.plot(df.index, cusum, label订正后) plt.axvline( pd.Timestamp(breaks[0][0].date()), colorgray, linestyle--, label断点 ) plt.xlabel(年份) plt.ylabel(累积距平) plt.legend() plt.show()CUSUM 对线性趋势也敏感。如果序列本身存在真实气候趋势累积和会持续上升不要一看到单调上升就认定是断点。判断断点要看斜率变化点而不是累积量本身的大小。提示SNHT 检测结果只能说明“这里存在统计显著跳变”不能证明“一定是非均一性”。先查元数据再下结论。6. 从最大风速延伸到降水、温度参数调整的四个关键位置同一条流程换到其他气象要素只需要调整四个位置其余代码可以原样复用。第一处是变换函数。最大风速直接取对数做差降水含大量零值先 log1p 再进入 SNHT避免零值取对数产生负无穷温度用原始值做差订正时也改为加性因子不做乘法。第二处是参考序列的平滑窗口。年度数据用 11 年窗口对 45 年序列够用月度数据窗口要短我用 25 个月既能抑制季节波动又不会把年际信号抹平。第三处是 SNHT 阈值。阈值本质上是显著性水平的映射序列越短 T 的自然波动越大。我习惯用蒙特卡洛模拟生成白噪声序列跑 1000 次 SNHT取 95% 分位作当前序列长度的阈值不做硬编码。第四处是订正基准段的选择。有元数据记录时以仪器更换后的段为基准没有元数据时以最近一个完整 10 年段为基准并在输出文件里注明“基准段依据统计推断待元数据复核”。还有一个容易踩的坑不要过度订正。SNHT 检测出断点时先对照观测簿确认那一年有没有换仪、迁站记录。完全对不上且断点前后均值差小于 0.5 m/s 时我一般标记为“疑似非均一”而不是直接订正这样后续做重现期分析时还有选择余地。订正完成后把断点年份、T 值、订正因子、是否与元数据吻合整理成一张表连同订正前后序列一并写入 CSV。这张表就是论文方法部分最有力的附件审稿人追问“怎么确定断点真实存在”时直接对表逐条说明即可。本文还有配套的精品资源点击获取