简介一套面向气象与水文研究者的风速均一化订正实操资料包以最大风速数据为例系统解决因仪器更换、测量方式改变或站点迁移带来的系统性偏差使长序列风速数据能够跨站点、跨时期可靠比较。包内共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, encodingutf-8) print(df.info()) print(df.head(10)) print(df.describe())读数据的逻辑很简单read_csv读取原始 CSVinfo()能看出每列的非空数量和数据类型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], errorscoerce) df[vmax] pd.to_numeric(df[vmax], errorscoerce) df df.set_index(year).sort_index() # 缺失值线性插值连续缺失3年以上则保留 NaN df[vmax_filled] df[vmax].interpolate(limit2, limit_directionboth) 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(limit2)表示最多连续插补两个缺失值超过两个的仍然保留为NaN后面断点检测时要跳过。limit_directionboth允许序列开头的缺失值也能用后面的值回推但年份序列开头缺数据其实很常见比如 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[k1:] - ranks[k]) u[k1] 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[k1]的递推式保证了计算复杂度在 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_star1] seg2 v[k_star1:] delta np.nanmean(seg2) - np.nanmean(seg1) corrected v.copy() corrected[:k_star1] corrected[:k_star1] 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 年的所有风速都要加上delta1999 年及之后不动。如果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, encodingutf-8-sig) print(output.head(10)) # 输出订正参数日志 with open(correction_log.txt, w, encodingutf-8) as f: f.write(fstation: YJ\n) f.write(fbreak_year: {df.index[k_star]}\n) f.write(fdelta: {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, window5): x np.asarray(x) t_stats [] for i in range(window, len(x) - window): seg1 x[i-window:i] seg2 x[i:iwindow] mean1, mean2 np.mean(seg1), np.mean(seg2) s1, s2 np.var(seg1, ddof1), np.var(seg2, ddof1) 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], labelRaw, alpha0.7, linewidth1.5) plt.plot(df.index, df[vmax_corr], labelCorrected, linewidth1.5) plt.axvline(df.index[k_star], colorred, linestyle--, labelBreak point) plt.legend() plt.title(YJ Maximum Wind Speed: Raw vs Corrected) plt.ylabel(Wind speed (m/s)) plt.grid(alpha0.3) plt.show()sliding_t里的window5表示左右各取 5 年做均值比较窗口太小容易受单年极值影响太大又会把断点附近的差异稀释。这个函数跑出来的 t 统计量如果围绕 0 附近波动说明序列平稳如果在某一年前后出现明显的尖峰则说明仍有残余的非均一性。可视化是必须做的一步我见过不少订正后的序列统计检验通过了但图上明显能看到整段断崖式抬升这种往往是订正量计算时混进了缺失值。绘图参数里alpha0.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, window5) years df.index[5:-5] print(最大 |t| 值对应的年份:, years[np.argmax(np.abs(t_stats))])从那以后我每次拿到一个新站点的风速序列都强制自己先把原始序列图、迁站记录、数据规格说明三样东西摆齐再跑订正脚本。订正不是数学游戏它是为了还原观测系统变更前的真实气候信息。希望帮到你。本文还有配套的精品资源点击获取
