简介这份资源提供基于Python实现的季节尺度M-K突变检测脚本面向从事气候、水文、生态及环境数据分析的科研人员与高年级学生用于识别具有季节性模式的时间序列中的趋势突变点。Mann-Kendall检验不依赖数据分布适合非正态序列脚本围绕SPEI干旱指数数据展开涵盖数据读取、缺失与异常值处理、STL季节分解、趋势成分的M-K检验以及突变点可视化等环节帮助读者掌握从原始季节数据到突变诊断的完整分析链路。压缩包共2个文件包含1个py脚本与1个xlsx数据表整体约11KB体量轻便便于直接运行与二次修改。目前已有195人学习下载适合希望快速上手季节尺度突变检测、并借助现成脚本验证自身数据的读者参考使用。1. 季节尺度 M-K 突变检测为什么月数据直接跑会翻车拿一份 30 年的逐月气温或降水序列直接丢进经典的 Mann-Kendall 突变检测里十有八九会得到一个让你怀疑人生的结果突变点飘忽不定换个起止年份结论就变甚至一年里能给你标出三四个突变。这不是代码写错了而是季节尺度数据本身带着强烈的周期性——月均温的年内波动幅度往往比几十年趋势的累积变化还大。M-K 检验的前提是序列独立同分布月数据显然不满足季节循环会把趋势信号淹没UF/UB 曲线交叉点自然就成了玄学。这篇要讲清楚的就是面对气象、水文、农业里最常见的月尺度或季尺度序列怎么用 Python 把 M-K 突变检测做对。核心思路是先做季节分离或去季节化再对每个季节或去季节后的残差单独跑 M-K最后把结果合起来看。适合手上有长序列观测数据、需要判断某年某季开始发生显著转变的从业者也适合刚学 Python 想找个真实场景练手的人。下面从原理、数据准备、代码实现到避坑一步步拆开。2. 季节尺度 M-K 的原理与数据预处理2.1 经典 M-K 检验在季节数据上为什么会失效Mann-Kendall 检验是一种非参数趋势检验它统计序列中所有后值大于前值的配对数量构造统计量 S再标准化成 Z。原假设是序列无趋势、随机独立。突变检测则用滑动方式正向序列算一遍 UF反向序列算一遍 UB两条曲线在置信区间内的交点就是潜在突变点。问题出在随机独立这个前提。月均温序列里1 月到 7 月单调上升7 月到次年 1 月单调下降这种确定性周期让 S 统计量被季节循环主导。你检测到的突变很可能只是某年某月异常偏高而不是气候态发生了转变。同理降水序列虽然年内周期没气温那么规整但汛期/非汛期的量级差异同样会干扰。常见做法有两种一是按季节拆开把每年春季3-5 月的值单独组成一个序列30 年就有 30 个点对这个序列跑 M-K二是先算每个月的气候平均值用原始值减去对应月份的气候值得到距平序列再对距平跑 M-K。前者物理意义清晰后者保留更多样本量。我一般两个都跑互相验证。2.2 数据准备从 CSV 到季节序列的转换假设你手上是一份station_data.csv两列date和value逐月记录。先做几件事解析日期、检查缺失、按季节分组。季节划分按气象惯例春 3-5 月夏 6-8 月秋 9-11 月冬 12-次年 2 月。冬季跨年要特别处理把 12 月归到当年还是次年会影响序列对齐我习惯把当年 12 月和次年 1、2 月算作同一冬季标记为起始年份。import pandas as pd import numpy as np # 读取数据parse_dates 自动解析日期列 df pd.read_csv(station_data.csv, parse_dates[date]) df df.sort_values(date).reset_index(dropTrue) # 检查缺失值M-K 对缺失敏感最好先插补或剔除 print(df.isnull().sum()) df df.dropna(subset[value]) # 提取年和月 df[year] df[date].dt.year df[month] df[date].dt.month # 定义季节映射冬季特殊处理 def get_season(row): m row[month] if m in [3, 4, 5]: return MAM, row[year] elif m in [6, 7, 8]: return JJA, row[year] elif m in [9, 10, 11]: return SON, row[year] else: # 12 月归入当年冬季1、2 月归入上一年冬季 if m 12: return DJF, row[year] else: return DJF, row[year] - 1 df[[season, season_year]] df.apply( lambda r: pd.Series(get_season(r)), axis1 ) # 按季节聚合取平均也可取总和看变量性质 seasonal df.groupby([season, season_year])[value].mean().reset_index() seasonal seasonal.sort_values([season, season_year])这段代码的关键在get_season函数12 月归当年1、2 月归上一年保证冬季序列连续。聚合方式用mean适合气温降水应该用sum。聚合后每个季节得到一个独立序列长度等于年数。如果某年某季缺一个月mean会自动忽略 NaN但样本代表性下降建议在聚合前先检查每月是否齐全。提示季节序列长度最好不少于 25 年M-K 在小样本下检验功效低交点位置不稳定。2.3 去季节化距平序列的构造如果你不想损失样本量可以构造距平序列。做法是对每个月份计算它在整个研究期内的气候平均值然后用原始值减去这个平均值。这样得到的距平序列去掉了年内周期可以整体跑 M-K。# 计算各月气候平均值 monthly_clim df.groupby(month)[value].mean().to_dict() # 构造距平 df[anomaly] df.apply( lambda r: r[value] - monthly_clim[r[month]], axis1 ) # 距平序列按时间排序后可直接用于 M-K anomaly_series df.sort_values(date)[anomaly].values距平法的好处是样本量是月数而非年数M-K 统计量更稳健坏处是如果趋势本身有季节差异比如春季升温快、秋季升温慢距平序列里仍残留季节信号。所以距平法适合先做初步判断最终结论还是建议分季节验证。3. 用 Python 实现季节 M-K 突变检测3.1 核心算法UF/UB 曲线的计算M-K 突变检测的核心是构造正序统计量 UF 和逆序统计量 UB。对长度为 n 的序列 x正序统计量按时间顺序计算对每个 i计算前面有多少个 j i 满足 x[j] x[i]记为 s_i累积和 S_i sum(s_1...s_i)均值 E(S_i) i(i-1)/4方差 Var(S_i) i(i-1)(2i5)/72标准化 UF_i (S_i - E(S_i)) / sqrt(Var(S_i))逆序统计量把序列反转后同样计算再反号对应回原时间轴。两条曲线在置信区间通常 ±1.96对应 0.05 显著性内的交点即为突变点。def mk_mutation_test(x, alpha0.05): 返回 UF、UB 两条曲线和突变点索引 x: 一维数组按时间顺序 alpha: 显著性水平默认 0.05 n len(x) # 正序 s np.zeros(n) for i in range(1, n): s[i] np.sum(x[:i] x[i]) S np.cumsum(s) E np.array([i * (i - 1) / 4 for i in range(1, n 1)]) Var np.array([i * (i - 1) * (2 * i 5) / 72 for i in range(1, n 1)]) # 避免除零 Var[Var 0] 1e-10 UF (S - E) / np.sqrt(Var) # 逆序 x_rev x[::-1] s_rev np.zeros(n) for i in range(1, n): s_rev[i] np.sum(x_rev[:i] x_rev[i]) S_rev np.cumsum(s_rev) UB_rev (S_rev - E) / np.sqrt(Var) UB -UB_rev[::-1] # 反号并反转回原时间轴 # 找交点UF 和 UB 差值变号处 diff UF - UB crossings [] for i in range(1, n): if diff[i - 1] * diff[i] 0: crossings.append(i) # 置信区间 from scipy.stats import norm z norm.ppf(1 - alpha / 2) return UF, UB, crossings, z参数说明x是输入序列长度至少 10alpha控制置信水平0.05 对应 z≈1.960.01 对应 z≈2.58。crossings返回的是索引位置对应到年份需要加上起始年。注意Var在 i1 时为 0代码里做了保护。交点判定用差值变号比直接找 UFUB 更稳健。3.2 分季节批量跑 M-K 并汇总结果有了核心函数就可以对每个季节序列分别调用把结果整理成表。results {} for season in [MAM, JJA, SON, DJF]: sub seasonal[seasonal[season] season].sort_values(season_year) x sub[value].values years sub[season_year].values if len(x) 10: print(f{season} 样本不足跳过) continue UF, UB, crossings, z mk_mutation_test(x) # 交点对应年份 cross_years [years[i] for i in crossings] results[season] { years: years, UF: UF, UB: UB, crossings: cross_years, z: z } print(f{season}: 突变点年份 {cross_years}, 置信水平 z{z:.2f})这段代码遍历四个季节对每个季节提取序列、调用检测函数、记录交点年份。crossings可能为空无显著突变或有多个需结合物理意义判断主突变点。实际项目里我一般还会输出 UF 和 UB 曲线图肉眼确认交点是否落在置信区间内。3.3 结果可视化UF/UB 曲线与突变点标注图是给别人看的也是给自己复核的。用 matplotlib 画 UF、UB 和 ±1.96 参考线交点用竖线标出。import matplotlib.pyplot as plt fig, axes plt.subplots(2, 2, figsize(12, 8)) for ax, season in zip(axes.flat, [MAM, JJA, SON, DJF]): if season not in results: continue r results[season] ax.plot(r[years], r[UF], b-, labelUF) ax.plot(r[years], r[UB], r--, labelUB) ax.axhline(r[z], colorgray, linestyle:) ax.axhline(-r[z], colorgray, linestyle:) for cy in r[crossings]: ax.axvline(cy, colorgreen, alpha0.5) ax.set_title(f{season} M-K) ax.legend() plt.tight_layout() plt.savefig(mk_seasonal.png, dpi150)axvline标出交点年份绿色竖线。如果交点落在 ±z 之外说明该点统计显著。图里如果 UF 和 UB 在置信区间内多次交叉不要急着下结论可能是序列自相关导致需要做预白化处理。4. 避坑与常见问题排查4.1 突变点一年一个样序列自相关在作怪现象同一份数据换不同起止年份突变点从 1995 年跳到 2003 年。原因M-K 假设序列独立但气温、降水序列普遍存在一阶自相关导致方差被低估交点位置不稳定。解决先做预白化pre-whitening用 AR(1) 模型去掉自相关再跑 M-K。Python 里可以用statsmodels估计自相关系数构造去相关序列。from statsmodels.tsa.ar_model import AutoReg def prewhiten(x): 一阶自相关预白化 model AutoReg(x, lags1, old_namesFalse).fit() rho model.params[1] x_pw np.zeros_like(x) x_pw[0] x[0] for i in range(1, len(x)): x_pw[i] x[i] - rho * x[i - 1] return x_pwrho是估计的一阶自相关系数x_pw是去相关后的序列。注意预白化会改变序列量级但 M-K 是非参数检验只关心秩所以不影响结论。4.2 冬季序列跨年错位12 月归属搞反现象冬季突变点比实际偏晚一年。原因12 月归到次年冬季导致冬季序列整体后移。解决按 2.2 节的get_season逻辑12 月归当年1、2 月归上一年。检查方法打印冬季序列的年份和对应月份确认 1990 年冬季包含 1989 年 12 月、1990 年 1 月和 2 月。4.3 缺失值插补引入假突变现象某年某季缺测插补后 M-K 在该年附近标出突变。原因插补值如均值填充人为制造了序列的跳跃或平坦段。解决缺失超过 10% 的年份直接剔除不要插补缺失少量时用相邻年份同季均值填充并在结果里标注。更稳妥的做法是只保留完整年份牺牲样本量换可靠性。4.4 置信区间选 0.05 还是 0.01现象0.05 下突变点很多0.01 下几乎没有。原因置信水平控制的是假阳性率气象序列噪声大0.05 容易把波动当突变。解决先用 0.05 初筛再用 0.01 复核。如果 0.01 下交点消失说明该突变不稳健不宜作为主要结论。我一般报告两个水平的结果让读者自己判断。4.5 多个交点怎么选主突变点现象UF 和 UB 在置信区间内交叉三次。原因序列存在多尺度波动M-K 对中期波动敏感。解决结合滑动 t 检验或 Pettitt 检验交叉验证优先选物理意义明确的年份如政策实施年、仪器更换年如果交点密集考虑对序列做 5 年滑动平均后再检测突出长期趋势。5. 进阶技巧滑动窗口 M-K 与多方法交叉验证单次 M-K 只能给一个突变点但气候序列的转变往往是渐进的。我习惯用滑动窗口 M-K设定窗口长度如 20 年从序列起点滑到终点每个窗口内跑一次 M-K记录突变点年份最后看突变点随窗口移动的变化。如果突变点在某段时期集中出现说明转变是持续的如果跳来跳去说明不显著。def sliding_mk(x, years, window20, step1): 滑动窗口 M-K返回每个窗口的突变点 n len(x) records [] for start in range(0, n - window 1, step): end start window seg x[start:end] seg_years years[start:end] UF, UB, crossings, z mk_mutation_test(seg) for c in crossings: records.append({ window_start: seg_years[0], window_end: seg_years[-1], mutation_year: seg_years[c] }) return pd.DataFrame(records)window一般取 20-30 年太短统计功效不足太长突变点被平滑掉。step取 1 年输出 DataFrame 后可以画散点图横轴窗口起始年纵轴突变年份看是否有一条水平带。交叉验证方面Pettitt 检验和贝叶斯突变检测都可以作为补充。Pettitt 检验在scipy里没有现成函数但实现简单核心是找使 Mann-Whitney 统计量最大的分割点。贝叶斯方法可以用PyMC或ruptures库后者专门做变点检测pip install ruptures即可。import ruptures as rpt # Pelt 算法检测变点 algo rpt.Pelt(modelrbf).fit(x) change_points algo.predict(pen10) print(变点索引:, change_points)pen是惩罚项越大变点越少需要根据序列长度和噪声调整。ruptures的结果和 M-K 对照如果两者一致结论可信度大幅提升。最后说个我自己的习惯每次跑完 M-K先把 UF/UB 图存下来隔一天再看一遍。当天看容易先入为主隔天看经常发现交点其实落在置信区间边缘或者某个季节的序列长度不够。这个后悔药能帮你少发很多有问题的结论。希望帮到你。本文还有配套的精品资源点击获取