简介：基于Python的季节尺度M-K突变检测脚本，面向气象、水文、环境科学等领域关注季节性时间序列突变点的研究者与分析人员。脚本以SPEI3.xlsx中的标准化降水蒸散发指数为例，完整演示了从数据读取、预处理、季节分解到Mann-Kendall突变检验的流程，并通过可视化定位潜在突变时刻，可直接替换为自有数据运行。压缩包共2个文件，含1个.py脚本与1个.xlsx数据文件，整体仅11KB，结构精简、便于阅读修改。脚本内置STL等分解思路，将季节性成分与趋势成分分离后再执行M-K检验，适合处理周期性明显的月度或季节数据；配套Excel文件提供示例数据，帮助快速验证方法效果。目前已有196人学习下载，适合正在学习非参数趋势检验、或需要复现干旱指数突变分析的Python使用者作为参考模板。通过动手实践，可掌握季节性数据突变检测的关键步骤与代码实现，节省从零搭建算法框架的时间。
1. 季节尺度 M-K 突变检测：为什么年平均会把突变点“藏起来”
做水文和气候序列分析的人，十有八九会遇到这种情况：年均值的 Mann-Kendall 趋势检验平平淡淡，看不出任何突变，但把数据按季节拆开一看，夏季序列在第 N 年附近有一个非常明显的跳变，只是被冬季和春季的平稳段“平均”掉了。这不是数据错，而是尺度选错了。年尺度的 M-K 突变检测，本质上在做“一整年的平均过程是否发生了均值结构变化”，而季节尺度的 M-K 是分别对春、夏、秋、冬四条序列各跑一遍序贯检验，得到各自的 UF/UB 曲线，用交叉点定位突变发生的时间窗口。这份基于 Python 实现的季节尺度 M-K 突变检测脚本，解决的就是“年尺度看不见、季节尺度才暴露”的那一类问题，适合手上拿着逐月或逐季数据、正在做径流、降水、气温或 NDVI 序列突变定位的从业者。
2. 序贯 M-K 的 UF/UB 怎么算：统计量、方差与分钟级最小实现
2.1 先搞懂序贯统计量：秩序列 S(k) 到底是什么
季节尺度 M-K 突变检测用的是“序贯 M-K”版本，也叫前进-后退 M-K。它和普通趋势 M-K 的最大区别是：趋势检验只给你一个终值 Z，而突变检测要给你一条曲线。曲线怎么来的？把原有时间序列从第 1 个样本逐步增长到第 n 个样本，每增加一个样本就计算一次 M-K 统计量，这样就能得到一个随样本长度变化的序列。
具体公式是这样的：设序列 x_1, x_2, …, x_n，定义秩序列 S(k) = sum_{i=1}^{k} sum_{j=1}^{i-1} sign(x_i - x_j)。sign 函数取值只有三种情况：x_i 大于 x_j 记 1，小于记 -1，相等记 0。换句话说，S(k) 统计的是“前 k 个样本里，后面出现的样本比前面样本更大或更小的次数差”。如果序列整体上升，越到后面 S(k) 越大；如果整体下降，S(k) 会持续偏向负值。S(k) 本身不是正态的，但可以标准化。
在零假设“序列没有趋势”下，S(k) 的期望是 k(k-1)/4，方差是 k(k-1)(2k+5)/72。标准化之后得到 UF(k) = [S(k) - k(k-1)/4] / sqrt( k(k-1)(2k+5)/72 )。UF 曲线就是这条标准化序列。反向序列从 x_n 倒推到 x_1，按相同方法计算，再把结果取负并镜像回原时间轴，得到 UB 曲线。UF 和 UB 在置信区间内发生交叉，交叉点位置就是突变发生时刻。置信区间一般取 α=0.05 对应的 ±1.96。
用大白话说：UF 是从头往尾看的累计趋势强度，UB 是从尾往头看的累计趋势强度。如果序列真的在第 m 年突变，从前往后看，突变前趋势不明显，突变后趋势强烈，UF 会在 m 附近剧烈抬升；从后往前看，UB 会在同一位置反向剧烈变动。两条曲线必然“打架”，“打架”的交点就是突变点。这个逻辑在季节尺度下完全复用，只是把“全部年份”替换成“某个季节的年份子序列”。
2.2 最小可用的核心实现：不带第三方统计库也能算
我写项目内部代码时，习惯把核心检验函数和控制绘图的部分分开，这样以后换数据集不用改算法层。下面的实现完全基于 NumPy，不依赖 SciPy 的统计函数，自己实现 tie 校正和顺序统计量。
import numpy as np def _sign_arr(d): # 向量化 sign：避免在双循环里反复判断 return (d > 0).astype(np.float64) - (d < 0).astype(np.float64) def _tie_correction(k, x): # 前 k 个样本中，等值组带来的方差修正项 # 公式：sum t_j*(t_j-1)*(2*t_j+5) / 72 seg = x[:k] _, counts = np.unique(seg, return_counts=True) t = counts[counts > 1] if t.size == 0: return 0.0 return float(np.sum(t * (t - 1) * (2 * t + 5))) / 72.0 def mann_kendall_forward(x): # 返回正向序列标准化后的 UF 曲线 x = np.asarray(x, dtype=np.float64) n = len(x) if n < 4: raise ValueError("样本量太短，序贯 M-K 没有意义") uf = np.zeros(n) for k in range(2, n + 1): # 计算前 k 个样本的 S(k) total = 0.0 for i in range(1, k): # 第 i 个样本与它之前所有样本比较 total += np.sum(_sign_arr(x[i] - x[:i])) s_k = total mean_k = k * (k - 1) / 4.0 var_k = k * (k - 1) * (2 * k + 5) / 72.0 var_k -= _tie_correction(k, x) uf[k - 1] = (s_k - mean_k) / np.sqrt(var_k) return uf def mann_kendall_sequential_full(x): # 正向 UF 与反向 UB 一起算 uf = mann_kendall_forward(x) # 反向序列：将 x 倒序后重复正向计算，再取负并反演 ub_tmp = mann_kendall_forward(x[::-1]) ub = -ub_tmp[::-1] return uf, ub这段代码的逻辑，我拆开说明一下。mann_kendall_forward里的外层循环 k 从 2 跑到 n，每次扩展一个样本进入“当前视野”。内层x[i] - x[:i]是拿当前样本和它之前全部样本比较，_sign_arr生成一组 -1、0、1 的向量，求和就是这一轮对 S(k) 的新增量。这样写虽然时间复杂度是 O(n^2)，但实际工程里季节子序列长度通常只有 30～60 个样本，跑起来没有任何压力。
_tie_correction是容易被忽略的。当序列里有大量等值样本，比如连续多个月降水量为 0 时，不修正方差会把 Var(S) 算小，UF 绝对值被放大，趋势显著性被高估。修正项就是前面公式里的 sum t(t-1)(2t+5)/72，这里的 t 是每个等值组的大小。
反向 UB 的构造有一个容易写错的地方：先倒序计算uf_tmp，得到的是“从尾巴开始往前累积”的标准曲线，然后要取负并[::-1]翻转回原始时间轴。顺序不能反，符号也不能丢。为什么取负？因为倒着看序列时，原来上升的趋势在反向视角里会变成下降，标准化之后符号相反，必须乘以 -1 才能让两条曲线在同一坐标系下有意义。
2.3 绘制和交叉点查找：NaN 开头与 ±1.96 边界
得到 UF/UB 之后，前两个点通常不参与展示，因为 k=1 和 k=2 时 S(k) 太小，标准化结果恒为 0，画出来会引出诡异的“零段”。我一般会把 uf[0]、uf[1] 直接置为 NaN，让 matplotlib 自动断开这两段。
找交叉点也不能只看“两条线有没有交点”。真正可靠的突变点需要满足两个条件：第一，UF 与 UB 在 ±1.96 置信线之间或附近发生符号切换；第二，交点之后 UF 或 UB 至少有一方继续朝同方向延伸，而不是反复震荡。下面这段是实际脚本里用的交叉点判断函数。
def find_crossing(uf, ub, crit=1.96, min_k=3): d = uf - ub crossings = [] for k in range(min_k, len(d) - 1): # 符号变化或直接相交 if d[k] * (d[k + 1]) < 0 or d[k] == 0.0: # 交叉点落在置信区间附近时才认为有效 if abs(uf[k]) <= crit or abs(ub[k]) <= crit: crossings.append(k) return crossings这个函数的参数说明：min_k=3跳过开头两个无效点，你也可以把它设成 5 或 10，取决于序列总长度。crit默认 1.96，对应 α=0.05。判断时我用“或”而不是“且”，因为实际数据中 UF 和 UB 往往不会恰好同时越过临界线，只要有一条在临界线附近，交叉位置就有参考价值。如果你追求更严格，可以把or改成and，代价是会漏掉一部分边缘突变。
3. 季节拆分与批量绘图：把逐月数据变成四条独立检验曲线
3.1 数据准备：按日历季还是按水文年拆分
季节尺度 M-K 的第一步，是把原始时间序列拆成四季的子序列。拆法不同，结果含义完全不同。最常见的是自然历法季度：3-5 月为春季，6-8 月为夏季，9-11 月为秋季，12-2 月为冬季。按年份将每年同一季节的均值或总量串起来，得到四条长度相等的子序列。
第二种是水文意义上的季节，比如把汛期（5-9 月）和非汛期（10-4 月）拆开，或者把水文年从某个月份起算。拆分时最关键的是“如何把 12 个月映射到季节标签”。我通常直接在 DataFrame 里加一列 season 标签，映射规则用一个字典一次定死。
import pandas as pd def assign_season(df, month_col='month', style='calendar'): """ style='calendar': 自然历法季度 style='flood': 汛期/非汛期，汛期为5-9月 """ if style == 'calendar': df['season'] = ((df[month_col] - 1) // 3 + 1).astype(int) elif style == 'flood': df['season'] = np.where(df[month_col].isin([5, 6, 7, 8, 9]), 1, 2) return df这里(month - 1) // 3 + 1把月份压缩成 1、2、3、4 四个数字，分别对应春夏秋冬。对于汛期/非汛期，我直接判断 5-9 月。注意 12-2 月跨年，单纯按自然年分组会把 2023 年 12 月和 2024 年 1 月分到两个不同年份，而这在气候统计里其实是同一个冬季。遇到跨年季节，建议在水文年数据里额外加一列hydro_year，规则是：如果 month 大于等于起始月，则该季节归属于当年，否则归属于上一年。
拆分前还要检查数据的连续性。逐月数据里有缺失值很正常，但季节均值的统计口径一旦变化，M-K 检验会跟着错。我一般先做 pivot，把年份、季节、数值整理成宽表，看每个季节有多少个非空年份；季节子序列少于 20 个样本时，后续统计量方差会很不稳定，后面避坑章节再细说。
3.2 批量循环四季：单脚本输出四条 UF/UB 曲线
核心算法没问题后，批量循环四季就很简单了。下面这段函数把四季的 UF/UB、交叉点和绘图都封装在一起，输出一个字典结果和一张 2×2 对比图。
from scipy import stats import matplotlib.pyplot as plt def seasonal_mk_analysis(df, alpha=0.05, style='calendar'): crit = stats.norm.ppf(1 - alpha / 2) results = {} # 按季节分组 for season in sorted(df['season'].unique()): sub = df[df['season'] == season].sort_values('year') x = sub['value'].to_numpy(dtype=np.float64) years = sub['year'].to_numpy() uf, ub = mann_kendall_sequential_full(x) # 前两个点无效，置为 NaN uf[:2] = np.nan ub[:2] = np.nan crossings = find_crossing(uf, ub, crit=crit, min_k=3) results[season] = { 'years': years, 'uf': uf, 'ub': ub, 'crossings': crossings } return results, crit这个函数的关键参数：alpha=0.05决定临界线宽度；style直接传给前面拆分函数，用来控制季节划分。循环内部，每个季节独立排序、独立跑检验、独立找交叉点，互不干扰。这是季节尺度 M-K 和年尺度 M-K 的最大差别——四条序列完全并行。
拿到results之后可以按你自己的绘图风格画图。我给出一个最少参数的绘图版本，方便你把脚本改成可执行文件时直接复用。
def plot_seasonal_mk(results, crit, output_path='seasonal_mk.png'): fig, axes = plt.subplots(2, 2, figsize=(12, 8)) season_names = sorted(results.keys()) for ax, s in zip(axes.ravel(), season_names): res = results[s] ax.plot(res['years'], res['uf'], label='UF', lw=1.6) ax.plot(res['years'], res['ub'], label='UB', lw=1.6) ax.axhline(crit, color='gray', linestyle='--', lw=1) ax.axhline(-crit, color='gray', linestyle='--', lw=1) for c in res['crossings']: ax.axvline(res['years'][c], color='red', linestyle=':', lw=1) ax.set_title(f'season {s}') ax.legend() fig.tight_layout() fig.savefig(output_path, dpi=300)绘图时我习惯把交叉点用竖直虚线标出来，这样能直观看出四季突变是否同步。之前跑某模拟项目 X 时，我印象最深的是：夏季交叉点落在第 29 年，秋季交叉点在同年，冬季没有交叉——这说明突变不是“全年一致变化”，而是季节性的，年尺度 M-K 根本捕捉不到这种结构性差异。
多季节同时检验隐藏着一个多重比较问题：4 条序列同时跑，纯随机情况下至少有一条出现假交叉的概率会高于 5%。如果数据量够，可以把单个季节的 α 下调为 α/4，也就是约 0.0125，对应临界线 2.5。对一般业务分析，我建议先按 0.05 看整体，再单独关注那些交叉点前后趋势方向一致的季节，不必一上来就做 Bonferroni 校正。
4. 季节 M-K 避坑记录：四类经典翻车与修复办法
以下每一条都是我在实际处理水文序列时碰到过的状况，按“现象 → 原因 → 解决”写，方便你排查自己脚本输出的异常结果。
4.1 交叉点总是出现在第三第四个样本附近，且每次换数据都变
现象：季节 UF/UB 曲线在序列开头约第 3～5 个点位置出现明显交叉，全年四个季节无一例外，交叉年份落在前五年内。但把数据起始年份改一下，交叉点又跟着跑到新的“开头”。
原因：序贯 M-K 的 S(k) 在样本量很小时方差极小，标准化统计量对前几个样本的相对大小异常敏感。如果早期恰好有两年极端高值或低值，UF 会被瞬间推向临界线之上，与偏态反向的 UB 构成一个“虚拟交叉”。这个交叉不代表真实突变，只是短序列的边界效应。
解决：把find_crossing的min_k提高到总长度的 10% 左右，比如 40 年序列至少从第 5 年开始找。更稳妥的做法是：要求交叉点出现后，UF 或 UB 至少在后续连续 5 个样本保持同方向，不满足就直接标记为“边缘交叉”。在脚本里加一行过滤条件len(uf[k:k+5]) == 5且符号方向一致即可。
4.2 降水序列里大量零值，UF 曲线冲上天花板
现象：干旱区冬季月降水常年为 0，季节均值序列里有超过三分之一的值等值。运行 M-K 后 UF 绝对值经常超过 5，甚至 10，远高于 1.96，图上曲线直接冲顶。
原因：等值组破坏了 S(k) 的方差结构。公式中 Var(S) 的 t(t-1)(2t+5) 修正项没实现，或实现有 bug，导致方差被低估。零值越多，低估越严重，UF 被过度放大。
解决：确保每个 k 都调用_tie_correction，且修正项参与减法的位置正确。我踩过坑之后养成的习惯是写一个自检：构造长度为 30 的全零序列，UF 应该为 0 或接近 0；再构造一半为 0 的序列，观察 UF 不会超过 3。如果全零序列 UF 异常，先查 tie 修正，再查np.unique在浮点精度下的等值判断是否失效。
4.3 季节子序列太短，只有 15～20 个样本
现象：拆开季节后，每年只有 3 个月的数据参与计算，40 年数据变成 40 个点还好；如果数据只覆盖 2005-2020 年，每个季节只剩 16 个样本，交叉点极不稳定，不同 α 下结论完全相反。
原因：序贯 M-K 是渐近正态统计，n 低于 20 时标准化公式的近似精度明显下降，尤其当序列存在一两个极值时，单个样本就能扭转整条 UF 曲线方向。
解决：一是合并相邻季节，把“冬季+早春”作为滑动窗口单元；二是给季节子序列做滑动段合成，把每 3 年合并成一段，序列长度可以拉伸到 40 以上，代价是突变时间分辨率下降到 3 年。前者保突变精度，后者保统计稳健性，按你业务对时间精度的要求选。脚本里我默认在seasonal_mk_analysis加了一个min_len=20的硬性提示，小于该值直接抛警告。
4.4 UF 和 UB 交叉很多次，不知道信哪一个
现象：图上三条虚线纵横交错，UF 和 UB 在临界线内外反复穿越，每次交叉看起来都“有效”，导致突变年份无法判断。
原因：序列本身可能由多个波动周期叠加而成，或者趋势不是单次突变，而是渐变过程。单次突变检测模型对这类数据天然失效。
解决：不要死盯单一交叉点。先看 UF 曲线是否在临界线以上持续超过 5 个样本，如果 UF 反复回落临界线，说明不存在清晰的均值突变，更应该考虑“趋势分段”而不是突变点分析。另一个技巧是输出交叉点两侧的均值差：计算交叉点前后各 5 年样本均值差，如果小于一个标准差，即使统计量交叉，实际突变幅度也可忽略。这个“幅度门槛”往往比统计显著更符合业务判断。
5. 验证一套检测结果对不对：置换检验与一个固定习惯
拿到交叉点之后，先别急着写报告。任何非参数统计量都有可能被偶然的序列结构骗到，我强烈建议在脚本里加一个置换检验，检验目的是：把季节序列随机打乱，得到的交叉数量会不会也和原始数据一样多。如果打乱后同样频繁出现交叉，说明这套检验结果更多是噪音而非突变。
rng = np.random.default_rng(2024) n_perm = 200 def permutation_check(x, crit, n_perm=200): uf0, ub0 = mann_kendall_sequential_full(x) obs_cross = len(find_crossing(uf0, ub0, crit=crit)) cnt = 0 for _ in range(n_perm): xp = rng.permutation(x) uf_p, ub_p = mann_kendall_sequential_full(xp) cross_p = len(find_crossing(uf_p, ub_p, crit=crit)) if cross_p >= obs_cross: cnt += 1 p_value = cnt / n_perm return obs_cross, p_value这段代码做了一件很朴素的事：原序列交叉点数量是“观测值”，随机置换序列里交叉点数量是“零分布”，p 值表示纯随机情况下出现同样数量交叉点的概率。p 值大于 0.05 时，我一般直接放弃该季节的突变解读；p 值很小，才把交叉点写进结论。
置换检验跑完之后，最后一个技巧属于使用习惯层面。我从那以后每次做季节尺度 M-K，都强制走一遍这个固定流程：先看原始逐月数据的季节均值，确认有没有缺失和异常峰；再跑核心检验，拿到 UF/UB；最后必须跑置换检验，并保留随机种子。为什么要固定这套流程？因为温度和降水数据往往带有明显的年际振荡，单靠一条 UF 曲线容易把周期波动误判为突变，置换检验恰恰能把“偶发性结构”和“系统性突变”区分开来。
这份基于 Python 的季节尺度 M-K 脚本拿到手后，你只需要改三处：数据文件路径、季节划分规则（日历季还是汛期/非汛期）、置信水平 alpha。改完直接运行就能得到四张季节突变图和对应的交叉年份表。希望帮到你。
本文还有配套的精品资源，点击获取
