简介：一套面向气象与水文研究者的风速均一化订正实操资料包，以最大风速数据为例，系统解决因仪器更换、测量方式改变或站点迁移带来的系统性偏差，使长序列风速数据能够跨站点、跨时期可靠比较。包内共3个文件，包括一个Python订正脚本、一份多年风速CSV原始数据表、一份PDF代码说明文档，整体仅377KB，结构清晰。目前已有592人学习，特别适合气象数据分析初学者以及需要复核历史资料的研究人员。通过实操可完整走通风速数据读取、缺失值插补、异常值识别、订正算法实现、结果验证与可视化绘图全流程；说明文档还对方法原理和参数调整做了补充，并且这套处理思路可延展到气温、降水等其他气象要素的均一化订正，为气候统计与研究提供扎实支撑。

1. 风速均一化订正：为什么最大风速序列最容易藏“虚假趋势”

做气象水文数据分析的人，最怕的不是缺数据，而是数据“看着连续、其实断过”。我拆过不少站点的最大风速序列，很多年份曲线乍看很平滑，但一旦知道某年迁站、某年换了测风仪，再看那条曲线就觉得处处可疑。以最大风速为例做均一化订正，是因为年极值对观测环境的系统性变化最敏感——迁站后风速突然整体变小、换仪器后风向代表性变了，这些都会让序列均值产生阶跃，直接用原始序列做趋势分析，结论很容易翻车。这份资源把最小可用的订正流程做成了一套 Python 脚本加数据，homogenization.py负责从读取数据到输出订正序列，YJ_vmax_1971-2015.csv是真实的 45 年单站年最大风速序列，配合代码说明.pdf可以完整跑通。适合气象水文从业者、数据工程师和学生，至少能帮你把“非均一序列”变成“可对比序列”，而不是拿到数据就画趋势线。

2. 数据读取与预处理：把 YJ_vmax_1971-2015.csv 变成可订正的序列

2.1 先摸清 CSV 的底：列结构、缺失值、异常值

压缩包里的 CSV 不是随便就能送进订正算法的，第一步永远是探明数据结构。我拿到文件的第一件事，是解压后用 Pandas 快速读一遍，看看列名、索引和统计描述，别急着跑订正。这一步能暴露很多问题，比如年份列被读成了字符串、风速列里混着NaN、甚至有的行是备注文字。

import pandas as pd df = pd.read_csv('YJ_vmax_1971-2015.csv', encoding='utf-8') print(df.info()) print(df.head(10)) print(df.describe())

读数据的逻辑很简单：read_csv读取原始 CSV，info()能看出每列的非空数量和数据类型，head(10)让你直观看到前几行的格式，describe()给出风速列的均值、标准差、最小最大值。如果df.info()里显示year为 object 类型，说明年份列带特殊字符或空格，需要先清洗；如果vmax列的非空数量小于总行数，说明存在缺失值。参数方面，编码建议先试utf-8，如果出现乱码再换gbk，老的气象站点数据经常会用到 GBK 编码。这个 CSV 里只有三个字段：year、vmax、station，其中station全是 YJ，单站序列。

2.2 预处理三件事：时间索引、缺失填补、异常截断

预处理的目的不是让数据变漂亮，而是让后续断点检测不被“脏数据”带偏。常见做法是先把年份设成索引，然后检查缺失值和异常值。缺失值我会分两级处理：连续缺 1 年的用前后两年线性插值；连续缺 3 年及以上的，插值没有任何意义，直接给这一整段打标记，宁可不参与断点检测，也不要硬补。异常值则要结合风速的物理上限来判断，单站年最大风速超过 80 m/s 的情况极少，如果出现基本可以怀疑是单位换算错误或人工录入问题。

df['year'] = pd.to_numeric(df['year'], errors='coerce') df['vmax'] = pd.to_numeric(df['vmax'], errors='coerce') df = df.set_index('year').sort_index() # 缺失值线性插值，连续缺失3年以上则保留 NaN df['vmax_filled'] = df['vmax'].interpolate(limit=2, limit_direction='both') df['gap_flag'] = df['vmax'].isna().astype(int) # 异常值筛查：超出物理范围直接置空 phys_upper, phys_lower = 80.0, 0.0 df.loc[(df['vmax_filled'] < phys_lower) | (df['vmax_filled'] > phys_upper), 'vmax_filled'] = None print(df[df['gap_flag'] == 1])

这里interpolate(limit=2)表示最多连续插补两个缺失值，超过两个的仍然保留为NaN，后面断点检测时要跳过。limit_direction='both'允许序列开头的缺失值也能用后面的值回推，但年份序列开头缺数据其实很常见，比如 1971 年没观测。异常值筛查的阈值设成 0 和 80 是保守方案，如果你处理的是台风多发区，80 m/s 以上也可能是真实值，这种情况下不做硬截断，改成打印警告手动复核。预处理完成后，vmax_filled列才是进入订正算法的输入，原始列保留不动，方便后面回退对比。

3. 订正算法核心：从断点检测到订正量计算

3.1 为什么选 Pettitt 检验：非均一性本质是均值突变

风速均一化订正的核心假设是：一条受非均一因素影响的风速序列，可以看成若干段均值不同的“平稳段”拼接而成。迁站会让整体风速普遍降低，换仪器会让记录风速普遍偏高，这些影响不会改变年与年之间的自然波动，但会整体抬高或压低某一时期的均值。所以订正不是去平滑曲线，而是找到这些均值发生阶跃的时间点，然后把各段对齐到同一基准上。

断点检测的方法有 SNHT、惩罚最大 t 检验、Pettitt 检验等。我会优先选 Pettitt，因为它是非参数方法，不要求风速序列服从正态分布，而年最大风速本身是极值序列，分布很偏，用参数假设反而容易出错。Pettitt 检验的做法是把序列遍历一遍，对每个可能的分割点计算前后两部分的秩和差异，取差异最大且统计显著的位置作为断点。它天然适合在资源里的这类 45 年单站序列上运行，计算量小，结果也直观。

import numpy as np from scipy.stats import rankdata def pettitt_test(x): x = np.asarray(x) n = len(x) ranks = rankdata(x) u = np.zeros(n) for k in range(n - 1): sign = np.sign(ranks[k+1:] - ranks[k]) u[k+1] = u[k] + np.sum(sign) k_star = np.argmax(np.abs(u[:-1])) kappa = np.abs(u[k_star]) p_value = 2.0 * np.exp(-6.0 * kappa**2 / (n**3 + n**2)) return k_star, p_value, u

这段代码实现了 Pettitt 检验的核心流程。rankdata把原始风速转成秩次，避免极值对距离计算的干扰；u[k+1]的递推式保证了计算复杂度在 O(n) 量级，45 年数据跑起来是毫秒级。输出的k_star是断点对应的索引位置，p_value用来判断这个断点是否显著，常见做法是取 0.05 作为阈值。需要注意，这个公式里的显著性水平有近似性，当p_value小于 0.05 时才认为存在可订正的非均一断点，否则序列视为均一，直接跳过订正。对于年最大风速这种长度通常不超过 100 的序列，这个近似足够用了。

3.2 实现断点检测与订正量计算

找到断点后，订正量本质上是断点前后两段均值的差值。在绝大多数场景里，我们默认最新的观测规范最可靠，因此以最后一段为基准，把前面各段整体加上一个订正量，让早期数据在均值上与近期对齐。这里有一个关键设定：订正量是一个常数，不是按比例缩放。风速均一化订正通常不做乘法，因为乘性订正会改变年际变率，而迁站和换仪器的影响主要是系统性的位置偏移，常数加法更稳。

def compute_correction(vmax, k_star): v = vmax.copy() seg1 = v[:k_star+1] seg2 = v[k_star+1:] delta = np.nanmean(seg2) - np.nanmean(seg1) corrected = v.copy() corrected[:k_star+1] = corrected[:k_star+1] + delta return corrected, delta k_star, p_value, u = pettitt_test(df['vmax_filled']) print(f"断点位置: {df.index[k_star]}, p值: {p_value:.4f}") if p_value < 0.05: df['vmax_corr'], delta = compute_correction(df['vmax_filled'].values, k_star) print(f"订正量: {delta:.3f} m/s") else: df['vmax_corr'] = df['vmax_filled']

这里seg1是断点之前的序列，seg2是断点之后的序列，delta是两段均值的差。要注意k_star是索引下标，对应的是df.index[k_star]那一年。比如断点定位在 1998 年，那么 1971-1998 年的所有风速都要加上delta，1999 年及之后不动。如果p_value不小于 0.05，说明没有显著断点，序列均一，不需要订正。这时强行订正反而会引入虚假趋势，这是初学者最容易犯的错。计算订正量时必须跳过NaN，也就是缺失补插后仍未覆盖的位置，否则均值会被空值拉小。

4. 运行订正与结果验证：一张图看出订正前后差异

4.1 执行订正：把订正量写回原始序列

订正量算出来了，剩下就是执行和输出。执行的过程很简单，但输出文件一定要保留原始值和订正值两列，方便后面核对。我会把结果写成一个新 CSV，包含year、vmax_raw、vmax_corr，同时把断点位置和订正量写进说明文件的附表里。这样别人拿到结果，可以清楚看到哪段被改动了，避免订正过程变成黑匣子。

df['vmax_raw'] = df['vmax'] df['vmax_corr'] = df['vmax'] # 先拷贝原始值 df.loc[df.index <= df.index[k_star], 'vmax_corr'] = df['vmax'] + delta output = df[['vmax_raw', 'vmax_corr']].dropna() output.to_csv('YJ_vmax_corrected.csv', encoding='utf-8-sig') print(output.head(10)) # 输出订正参数日志 with open('correction_log.txt', 'w', encoding='utf-8') as f: f.write(f"station: YJ\n") f.write(f"break_year: {df.index[k_star]}\n") f.write(f"delta: {delta:.3f} m/s\n")

utf-8-sig编码是给 Excel 用户准备的，直接双击 CSV 不会出现乱码。订正之后的vmax_corr列，断点之前的数值全部被平移，断点之后保持原样。这里要特别说明：df['vmax_corr'] = df['vmax'] + delta看起来很直接，但必须用df.loc限定行范围，df.index <= df.index[k_star]这个条件等价于“断点当年及之前”。如果断点正好在某一年中间，年序列数据没有月信息，只能把整年归到前半段。订正参数日志是我个人的习惯，哪怕结果不理想，也能知道是哪一步出了问题。

4.2 验证：滑动 t 检验与可视化的双重确认

订正完不能直接拿去画趋势线，要先验证订正是否合理。验证至少要做两件事：一是重新在订正后序列上跑一遍断点检测，看原来的断点是否消失；二是画图对比订正前后，看序列是否出现了不自然的阶梯感。我会用滑动 t 检验作为补充，它和 Pettitt 检验不同，能显示每个时间窗口内的均值差异，对于确认断点位置很直观。

import matplotlib.pyplot as plt def sliding_t(x, window=5): x = np.asarray(x) t_stats = [] for i in range(window, len(x) - window): seg1 = x[i-window:i] seg2 = x[i:i+window] mean1, mean2 = np.mean(seg1), np.mean(seg2) s1, s2 = np.var(seg1, ddof=1), np.var(seg2, ddof=1) se = np.sqrt(s1/window + s2/window) t_stats.append((mean1 - mean2) / se) return np.array(t_stats) plt.figure(figsize=(12, 5)) plt.plot(df.index, df['vmax_raw'], label='Raw', alpha=0.7, linewidth=1.5) plt.plot(df.index, df['vmax_corr'], label='Corrected', linewidth=1.5) plt.axvline(df.index[k_star], color='red', linestyle='--', label='Break point') plt.legend() plt.title('YJ Maximum Wind Speed: Raw vs Corrected') plt.ylabel('Wind speed (m/s)') plt.grid(alpha=0.3) plt.show()

sliding_t里的window=5表示左右各取 5 年做均值比较，窗口太小容易受单年极值影响，太大又会把断点附近的差异稀释。这个函数跑出来的 t 统计量如果围绕 0 附近波动，说明序列平稳；如果在某一年前后出现明显的尖峰，则说明仍有残余的非均一性。可视化是必须做的一步，我见过不少订正后的序列，统计检验通过了，但图上明显能看到整段断崖式抬升，这种往往是订正量计算时混进了缺失值。绘图参数里，alpha=0.7让原始序列显得淡一些，订正后的线用实线突出，红色虚线标注断点位置。跑完图之后，我会同时算一下订正前后标准差的变化，如果订正后标准差明显增大，不一定是坏事，但需要检查是不是订正量过大导致早期段被过度抬升。

5. 避坑：风速均一化订正的几个高频翻车点

5.1 把“极端年”当成断点：台风年带来的误判

现象：序列里某一年风速特别大，比如 12 级台风过境，Pettitt 检验直接把那一年定位成断点，前后两段均值差巨大。原因：年最大风速本身就是极值序列，单个极端台风年的值会让非参数检验的秩分布发生突变，检验算法分不清“自然极值”和“观测环境改变”。解决：先对序列做异常值筛查，把超过 3 倍标准差的值标记出来，然后对比去掉该年后的断点位置。如果去掉后断点消失，说明这个断点是由极端年伪造成的，不能采信。常见做法是对风速取对数再跑 Pettitt，对数变换能压制极端值的影响，通常能得到更稳的断点定位。

5.2 订正基准选错：把早期数据全改成近期的“平均主义”

现象：订正后的序列变成一条水平线，长期趋势完全消失，和历史气候记录矛盾。原因：基准段选错了。有些脚本会把整个序列的均值作为基准，然后把每一段都向全序列均值对齐，结果就是所有段都被人为拉到同一条水平线上，真实的气候趋势被抹掉。解决：订正基准必须选择当前仍在执行的观测规范所对应的时段。比如 YJ 站 2003 年迁站，那么 2003 年以后的观测规范是现行的，订正量应该让 2003 年以前的数据向 2003 年以后的均值靠拢，而不是让所有年份向 45 年均值靠拢。判断基准的标准只有一个：哪个时间段的观测条件最接近现在，就以哪段为基准。

5.3 插值破坏序列尾部特性：缺失插补后断点偏移

现象：原始数据缺失 3 年，线性插补后跑断点检测，断点从 2003 年跑到了 2005 年，两个位置都有统计显著性。原因：线性插值本质上是在两个观测值之间拉一条直线，如果这段跨越了真实的断点，插补出来的值就带着“假趋势”，把断点的位置扯偏了。解决：连续缺失超过 2 年的位置不参与断点检测，检测时把缺失段直接留空，检测完成后如果需要完整序列，再用插补值填回去，但订正量计算只基于非插补值。另一个更稳的办法是先用没有缺失的年份做断点检测，然后再对缺失段单独插值，这样插补不会影响断点定位。

5.4 单位与规格混乱：最大风速与极大风速混用

现象：订正后序列在某个年份出现不自然的跳变，但那个年份并没有迁站记录。查元数据发现，1996 年之前记录的是 10 分钟平均最大风速，1996 年之后换成 3 秒瞬时极大风速，两者数值上差了一截。原因：风速观测有两个常见规范，一个是平均时距，一个是瞬时极大值，混在同一个序列里，任何统计检验都会把它当作断点。解决：在跑订正之前，先核对数据说明文档。代码说明.pdf里如果注明了各年份的风速规格，就要先把瞬时风速统一换算成平均风速，或者反过来。换算系数在不同高度和地形下差别很大，不能直接用简单的比例因子，最好找同站同时段的两套观测做回归。这一步没做好，后面所有订正都是白算。

6. 进阶：把订正推广到月最大风速与多站点互检

年最大风速序列只是入门，实际工程里更常见的是月最大风速序列，甚至日最大风速。把homogenization.py里的输入从年序列改成月序列，Pettitt 检测部分完全不用动，但要注意月序列存在明显的季节周期，直接跑 Pettitt 会把每年的季节性跳变误认为断点。常见做法是先按月份拆成 12 个子序列，分别做订正，最后再拼回去，或者对序列做年内的去季节化，比如用逐月距平代替原始值。多站点场景下，更可靠的是相对订正：选一个长序列参考站，用目标站与参考站的差值序列来检测断点，因为差值序列已经自然剔除了共同的气候波动，剩下的突变几乎都是站点自身变化引起的。这份资源里的单站流程可以看作这个思路的简化版，你把vmax换成vmax_target - vmax_ref，后面的订正步骤一模一样的。

验证订正效果时，我最近的习惯是跑一段滑动 t 检验作为“后悔药”检查：订正后的序列如果还有显著的 t 统计量尖峰，说明断点没订干净，要回头调整订正量或者检查基准段。更进一步，可以引入标准正态均一性检验的现成实现，像Climatol这类 R 包里的homogen()函数，能一次处理多站数据，但原理和我上面拆的 Pettitt 是相通的。写代码验证时注意滑窗大小的选择。

# 订正后序列的滑动 t 检验，窗口设为5年 t_stats = sliding_t(df['vmax_corr'].dropna().values, window=5) years = df.index[5:-5] print("最大 |t| 值对应的年份:", years[np.argmax(np.abs(t_stats))])

从那以后，我每次拿到一个新站点的风速序列，都强制自己先把原始序列图、迁站记录、数据规格说明三样东西摆齐，再跑订正脚本。订正不是数学游戏，它是为了还原观测系统变更前的真实气候信息。希望帮到你。

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