企业官网建设流程全解析

简介：这份资料面向气象与水文领域的研究人员及数据分析学习者，聚焦最大风速的均一化订正问题，帮助消除因仪器更换、测量方法调整或站点迁移带来的系统性偏差，使不同站点与时期的风速记录具备可比性。压缩包共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、贝叶斯方法怎么选

实际业务和科研中，最大风速均一化订正的主流方法有以下几类：

我一般会先用 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(drop=True) # 检查缺测年份 full_years = pd.DataFrame({'year': range(df_target['year'].min(), df_target['year'].max() + 1)}) df_target = full_years.merge(df_target, on='year', how='left') # 标记缺测 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( method='linear', limit=3, limit_direction='both') # 如果缺测太多，考虑剔除该站或使用更复杂的插补方法 if df_target['is_missing'].sum() > len(df_target) * 0.1: print("警告：缺测率超过10%，订正结果可靠性下降")

这段代码做了三件事：读取数据、检查缺测、线性插补。参数说明：limit=3表示最多连续插补 3 年，超过 3 年的缺测不插补，因为线性插值在长缺测段会引入虚假趋势。limit_direction='both'表示序列首尾的缺测也尝试插补。如果缺测率超过 10%，建议换站或使用更复杂的时空插补方法。

3.2 Pettitt 检验识别突变点

Pettitt 检验是一种非参数突变点检测方法，不要求数据服从特定分布，适合最大风速这种偏态序列。

def pettitt_test(x, alpha=0.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(t+1): for j in range(t+1, 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]统计的是前 t+1 个点和后 n-t-1 个点之间的符号差异累积。K取最大绝对值，p_value用近似公式计算。参数alpha=0.05是显著性水平，p 值小于 0.05 认为突变显著。注意：Pettitt 检验只能检出一个突变点，如果序列有多个断点，需要分段检验或改用其他方法。

3.3 SNHT 检验与参考站对比

SNHT 需要参考站序列。假设已经选好参考站，数据在同一 CSV 中。

def snht_test(target, reference, alpha=0.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) # 临界值（近似，n>20时） # 更精确的临界值需要查表或蒙特卡洛模拟 T_critical = 8.0 # 对应 alpha=0.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']], on='year', 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(f"SNHT 突变点年份: {merged['year'].iloc[cp_snht]}, T={T_max:.2f}, 显著: {sig_snht}")

SNHT 的核心是构造一个统计量 T，衡量突变点前后两段均值的差异。diff是待检站和参考站的差值序列，标准化后消除量纲影响。T[t]在突变点处达到最大。临界值T_critical=8.0是 n>20 时的近似值，更精确的做法是用蒙特卡洛模拟生成临界值表。注意：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, color='gray', linestyle='--', label='突变点') axes[0].set_xlabel('年份') axes[0].set_ylabel('最大风速 (m/s)') axes[0].legend() axes[0].set_title('序列对比') axes[1].hist(old_data, bins=10, alpha=0.5, label='旧时段', density=True) axes[1].hist(new_data, bins=10, alpha=0.5, label='新时段', density=True) axes[1].hist(corrected_old, bins=10, alpha=0.5, label='校正后旧时段', density=True) axes[1].set_xlabel('最大风速 (m/s)') axes[1].legend() axes[1].set_title('分布对比') plt.tight_layout() plt.savefig('homogenization_result.png', dpi=150) 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 检验做非参数趋势检验，再用 Sen's slope 估计趋势幅度。Mann-Kendall 检验不要求正态分布，对最大风速这种偏态序列比较稳健。Sen's slope 比线性回归斜率更抗 outliers。如果两种方法给出的趋势方向一致，基本可信。如果不一致，说明序列里还有未识别的断点或 outliers。

import pymannkendall as mk # 对订正后的序列做 Mann-Kendall 趋势检验 result = mk.original_test(corrected_series) print(f"趋势: {result.trend}, p值: {result.p:.4f}, Sen's slope: {result.slope:.4f}") # 如果 p < 0.05，趋势显著 # Sen's 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(f"50年一遇最大风速: {wind_50:.2f} m/s") # 轮廓似然法估计置信区间（简化版） # 实际应用中建议用专门的极值分析包，如 pyextremes

这里有个血泪经验：订正后的序列在做极值推算时，不要直接用原始序列的极值。订正可能会改变极值的大小，尤其是分位数映射校正后，旧时段的极值可能被“拉”到新时段的分布上。如果旧时段的极值被拉得过高，重现期估计会偏大。我一般会对比订正前后极值推算结果，如果差异超过 20%，就需要回头检查订正过程。

最后一个技巧：订正后的序列建议保留“订正标记”。在数据表里加一列is_corrected，标记哪些年份的数据被校正过。这样后续分析时，可以区分“原始观测”和“校正值”，在做敏感性分析时能快速切换。这个习惯帮我省了很多后悔药——有一次合作方质疑订正结果，我直接拉出标记列，对比了订正前后的趋势，问题当场定位。

希望帮到你。

本文还有配套的精品资源，点击获取