简介：这份资源提供基于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(drop=True) # 检查缺失值，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)), axis=1 ) # 按季节聚合，取平均（也可取总和，看变量性质） 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']], axis=1 ) # 距平序列按时间排序后可直接用于 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)(2i+5)/72
- 标准化 UF_i = (S_i - E(S_i)) / sqrt(Var(S_i))
逆序统计量把序列反转后同样计算，再反号对应回原时间轴。两条曲线在置信区间（通常 ±1.96，对应 0.05 显著性）内的交点即为突变点。
def mk_mutation_test(x, alpha=0.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是输入序列，长度至少 10；alpha控制置信水平，0.05 对应 z≈1.96，0.01 对应 z≈2.58。crossings返回的是索引位置，对应到年份需要加上起始年。注意Var在 i=1 时为 0，代码里做了保护。交点判定用差值变号，比直接找 UF==UB 更稳健。
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-', label='UF') ax.plot(r['years'], r['UB'], 'r--', label='UB') ax.axhline(r['z'], color='gray', linestyle=':') ax.axhline(-r['z'], color='gray', linestyle=':') for cy in r['crossings']: ax.axvline(cy, color='green', alpha=0.5) ax.set_title(f'{season} M-K') ax.legend() plt.tight_layout() plt.savefig('mk_seasonal.png', dpi=150)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, lags=1, old_names=False).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, window=20, step=1): """滑动窗口 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(model='rbf').fit(x) change_points = algo.predict(pen=10) print('变点索引:', change_points)pen是惩罚项，越大变点越少，需要根据序列长度和噪声调整。ruptures的结果和 M-K 对照，如果两者一致，结论可信度大幅提升。
最后说个我自己的习惯：每次跑完 M-K，先把 UF/UB 图存下来，隔一天再看一遍。当天看容易先入为主，隔天看经常发现交点其实落在置信区间边缘，或者某个季节的序列长度不够。这个后悔药，能帮你少发很多有问题的结论。希望帮到你。
本文还有配套的精品资源，点击获取
