季节尺度M-K突变检测的Python实现与代码拆解

简介：基于Python实现季节尺度曼肯德尔突变检测的完整脚本，面向气候与环境数据分析者，也适用于水文、生态等时间序列突变研究。脚本利用曼肯德尔非参数检验法识别时间序列中的显著趋势突变点，该方法不依赖数据分布，对非正态数据依然稳健；同时结合基于Loess的季节趋势分解处理周期性波动，可对标准化降水蒸散发指数逐季节进行突变检测，并自动标记突变发生年份与显著性水平。资源共包含两个文件，一个为可执行的Python主程序，一个为配套的Excel示例数据表，压缩包大小仅为11KB，非常轻量。程序完整覆盖数据读取、缺失值处理、季节分解、突变统计量计算及结果可视化，可直接运行验证，帮助理解季节尺度突变检测的原理与编码实现。已有195人学习使用，适合需要快速开展气候突变分析的研究者与学生。

1. 季节尺度 M-K 突变检测：Python 实现思路与代码资源拆解

做气象水文序列分析的人，基本都绕不开 Mann-Kendall 检验。全网搜 python M-K 突变检测，教程不少，但绝大多数只教你处理全年序列，一说到季节尺度，要么含糊带过，要么直接把四季数据拼在一起跑——那结果就是突变点分布成了玄学，UF 和 UB 曲线交叉得像一团乱麻，没法跟实际气候背景对上。这份基于 Python 实现季节尺度 M-K 突变检测的代码资源，解决的核心问题就是：把序列按季节拆分，分别计算统计量，再叠合展示，让每个季节的突变位置独立且可解释。它面向的是我这类需要跟降水、气温、径流序列打交道的气象水文从业者，也适合写论文需要汇报突变年份的学生。结论先行：季节性拆分不是简单分组，正序反序序列的定义窗口、交点判读规则、季节样本量不足时的处理策略，才是这个代码真正值钱的地方。

2. 为什么季节尺度要单独处理：M-K 方法在分组序列上的边界

2.1 M-K 检验的数学前提与季节分组的冲突

先回顾经典 Mann-Kendall 检验的原理。它本质上是基于秩的非参数检验，通过比较序列中所有元素对的大小关系来构造 Kendall 秩相关系数。对时间序列 xt（t = 1, 2, ..., n），统计量 S 定义为：

S = Σ(i<j) sign(xj - xi)

其中 sign 函数只在 xj > xi 时取 1，两者相等取 0，xj < xi 时取 -1。S 的期望值为 0，方差在序列不存在趋势的假设下可以用 n(n-1)(2n+5)/18 近似。

这地方有个容易被忽略的前提：M-K 检验假设数据是独立同分布的，并且序列中不存在季节周期性。如果直接把 30 年逐月数据混在一起跑，周期性会污染秩比较，导致 S 统计量被季节性波动主导，趋势信号被淹没。我做某流域降水序列分析时犯过这个错，全序列 M-K 检测显示显著上升趋势，按季节拆开看，夏季降水根本无显著趋势，趋势完全来自春季那一支。

季节尺度的做法是把序列按自然季节（通常 3 个月一组）切成四个子序列，每个子序列分别做 M-K 检验。这个拆分的意义不只是统计口径，更关键是：突变检验使用正序和逆序两个方向累积计算，当子序列长度只有原序列 1/4 时，正序、逆序序列的构建方式直接影响交点数量和位置。

2.2 季节尺度的输入数据格式约定

这份代码对输入数据有明确的格式预期，我在拆解时整理了一份最小样例：

注意，代码不做缺失年份补齐。如果你的数据有断年，比如 1995 年缺一整年，最简单的处理是先把该年份置为 NaN，再调用代码的预处理函数做线性插值；如果缺的年份太多（超过总年份的 10%），建议放弃该站点的季节突变分析。

我一般习惯把数据整理成 CSV 格式，三列无表头或带表头均可，代码里用 pandas read_csv 读取后统一把列名映射为 year、month、value。读取后第一件事是检查月份范围：

# 检查月份合法性 if not ((raw_data['month'] >= 1) & (raw_data['month'] <= 12)).all(): raise ValueError('月份必须为 1 到 12 之间的整数')

这段预校验很有必要，我见过一份数据里月份列混入了 13 和 0，结果按季节分组时，1 月的数据归到冬季，13 月的数据被丢进 NaN 组，直接影响后续统计量的秩计算。

2.3 季节划分边界：固定季节还是滑动窗口

代码默认按固定自然季节分组：

• 春季（MAM）：3 月到 5 月
• 夏季（JJA）：6 月到 8 月
• 秋季（SON）：9 月到 11 月
• 冬季（DJF）：12 月到次年 2 月

需要注意，12 月属于冬季，但 12 月来自上一年。代码在实现时把 year 和 month 拼接成时间索引，然后通过dt.quarter直接映射季度，但这里有个隐性约定：冬季是从当年 12 月开始的。比如 2020 年 12 月、2021 年 1 月和 2 月会被归为同一个冬季组，跨了年份。代码内部处理方式是把 12 月的年份减 1，统一到冬季结束年份，这样每个季节组内的年份是连续且对齐的。

如果你要的是水文学意义上的水文年季节划分（比如从 10 月到次年 9 月），这个代码不直接支持。我改过一版，把季节映射函数改成可传入偏移参数的模式：

def assign_season(month, offset=0): # offset=0 为自然年季节，offset=1 表示从 10 月开始划分水文年季节 adjusted_month = month - offset if adjusted_month <= 0: adjusted_month += 12 ...

对降水序列来说，水文年划分往往比自然年更合理，但这个改动会牵涉到后面 UF 序列的累积窗口构建，不仅仅是分组函数的问题，建议只在数据量充足（每季 30 年以上）时尝试。

3. 核心算法实现：UF、UB 序列如何正确计算

3.1 正序累积统计量 UF 的逐季算法

M-K 突变检测的核心输出是两个序列：UF（正序累积统计量）和 UB（逆序累积统计量）。两者的交叉点就是突变发生的候选位置。这个代码的实现思路是：在每个季节的子序列内，逐样本递增地计算累积统计量。

UF 的递推公式是这样走的：对季节内子序列 x1, x2, ..., xm，在第 k 步计算：

S_k = Σ(i=1..k-1) Σ(j=i+1..k) sign(xj - xi)

这个双层累加在代码里不是直接暴力循环，而是利用 Python 的itertools.combinations把前 k 个元素的所有组合生成出来再统计，时间复杂度 O(k^2)。对于 30 年序列，一个季节只有 30 个数据点，组合数 C(30,2) = 435，暴力计算完全可接受。我用 70 年序列试过，单季运行时间也就毫秒级，不需要优化成 Fenwick 树。

关键代码片段如下：

def compute_uf_subseries(sub_series): """计算正序累积统计量 UF 序列""" n = len(sub_series) uf = [] s_k = 0.0 var_s_k = 0.0 for k in range(2, n + 1): # 所有前 k 个元素的组合 comb = list(itertools.combinations(range(k), 2)) sign_sum = 0 for i, j in comb: diff = sub_series[j] - sub_series[i] if diff > 0: sign_sum += 1 elif diff < 0: sign_sum -= 1 s_k = sign_sum var_s_k = k * (k - 1) * (2 * k + 5) / 18.0 if var_s_k == 0: uf.append(0.0) else: if s_k > 0: uf.append((s_k - 1) / math.sqrt(var_s_k)) elif s_k < 0: uf.append((s_k + 1) / math.sqrt(var_s_k)) else: uf.append(0.0) return uf

逻辑说明：这里循环从 k = 2 开始，因为前 1 个元素无法构成元素对。s_k是第 k 步的累积 S 值，var_s_k是第 k 步的方差。分子加 1 或减 1 是连续性修正，这是 M-K 突变检测的标准做法，不加修正会导致 UF 值在临界区附近整体偏小，交点位置容易滞后。uf列表的索引 0 对应第 2 个数据点的统计量，在后续画图中要注意这个偏移。

参数说明：sub_series是无缺失值的季节子序列。代码要求每个季节子序列的缺失值已在前置处理中插补，否则会出现组合索引越界。我自己在实测中发现，如果子序列里有 NaN，diff < 0和diff > 0的判断都会失效，比较结果恒为 False，最终 sign_sum 是 0，UF 序列会平白多出一段 0 值，干扰交点位置判读。

3.2 逆序 UB 序列与序列反转的陷阱

UB 序列的计算是这份代码最容易踩坑的部分。教科书上的标准做法是：把序列反转，对反转序列重新做一遍正序 M-K 累积，再对结果的 UF 序列做反向翻转和平移，得到 UB。这个描述说三句话就完了，但实现细节差一毫米结果差千里。

代码中的实现方式是：

def compute_ub_from_reversed(uf_reversed): """从反转序列的 UF 结果还原 UB 序列""" n = len(uf_reversed) ub = [0.0] * n for i in range(n): ub[i] = -uf_reversed[n - 1 - i] return ub

这里的关键坑在于：反转序列的 UF 计算的起始端是原序列的末端。举例说，原序列是 1991 到 2020 年 30 个值，反转后第一个值是 2020 年，那么 UF 的第一个值反映的是 2020 年与之前年份的比较，而不是 2020 年与之后年份的比较——因为反转后没有“之后”了。正确的解读是：UB 在第 k 个点的值，代表从末端反向累积到第 k 个点时的统计显著性。

我的实践体会是，调试阶段最有效的验证方式是用一个已知突变的人造序列来检查 UB 曲线的形态。比如构造一段前 15 年均值 50、后 15 年均值 80 的序列，理想的 UF 曲线在突变点前缓慢攀升或波动，在突变点后快速上升突破临界线；UB 曲线在突变点后向下突破。如果你的 UB 曲线形态完全反过来，说明反转时把-uf_reversed[i]的符号或索引弄错了。

这段逻辑放在独立函数里是有意的，代码作者在注释里专门标注了这一点，可见这也是他踩过坑之后的总结。

3.3 季节子序列的样本量与 UF 起点偏移

季节尺度处理的一个特有问题是样本量。总序列 30 年时，每个季节只有 30 个样本，UF 序列从第 2 个点开始计算。如果数据从 1991 年到现在只有 20 年，那么春季子序列只有 20 个点，UF 序列前段（前 5 年左右）因为统计量累计不足，方差小，UF 值波动特别剧烈，经常会在没有突变的位置打出虚假交点。

我建议在做季节 M-K 突变检测前，先对每个季节的子序列长度做检查：

if len(sub_series) < 15: print(f'警告：{season}季样本量 {len(sub_series)} 不足 15，突变点结果仅供参考')

15 是经验阈值。教科书上对 M-K 检验的样本量要求一般建议 n ≥ 10，考虑到季节子序列的独立样本更少，我一般把红线设在 12，低于 12 直接放弃该季节的突变判断，只报告趋势方向。这不是原代码的默认行为，但我强烈建议在调用时加这个判断——尤其在处理高原站点、极地区域站点时，很多数据序列本身就只有十几年。

还有一个容易忽略的起点偏移问题。UF 序列的第一个有效值对应原始序列的第二个数据点。在绘制 UF/UB 曲线时，横坐标应该是该数据点对应的年份。比如春季序列包含 1991 到 2020 年的 3 个值，第一个 UF 值对应 1992 年。代码里在处理横坐标时用season_years[1:]来对齐，如果你直接用 0 到 n-1 的整数索引，画出来的突变点年份会比实际早一年。

4. 突变点判读与可视化：交点到结论的完整链路

4.1 交点判读逻辑：UF 与 UB 交叉的三种情况

完成了 UF 和 UB 序列的计算，真正的挑战才开始：怎么判断突变点？这套代码内置了三条判断规则，我会逐一拆开说。

第一条规则是基础过滤：只保留 UF 与 UB 曲线在置信区间内的交点。置信区间取的是显著性水平 α = 0.05 对应的 ±1.96。两条曲线可能在区间外频繁交叉，但这种交叉没有统计意义。代码里先用布尔索引筛出 UF 和 UB 同时处于 [-1.96, 1.96] 之内的区间，再在这个区间段内搜索交点。

第二条规则处理同一季节内多个交点的情况。气象序列现实里很常见：UF 和 UB 在一个季节里交叉四五次。代码的默认策略是取连续交叉段中最靠前的交点作为突变候选点，但因为 UF 序列前段方差小、波动大，最靠前的交点往往是个假信号。我一般把代码改成只考察 UF 与 UB 在临界线内侧交叉、且之后 UF 持续上升或下降超过 1.96 的那一段，取该段的起始交叉位置。这样过滤后，多个交点通常能收敛到 1 到 2 个有效候选。

第三条规则是方向约束：突变点必须与趋势方向一致。假设春季序列整体是升温趋势，UF 曲线最终在临界线上方，那么突变点在 UF 从负转正并突破 0 线的位置附近，而不是 UF 在负值区域交叉的位置。代码里会在候选点列表中检查：如果 UF 终点值大于 1.96，则只保留候选点中 UF 从负变正的那个索引；反之亦然。

规则落实到代码是这样的：

def locate_mutation_point(uf, ub, years, conf_level=1.96): candidates = [] in_conf = (np.abs(uf) <= conf_level) & (np.abs(ub) <= conf_level) for i in range(1, len(uf)): # 检测 UF 与 UB 的符号交叉 if np.sign(uf[i] - ub[i]) != np.sign(uf[i-1] - ub[i-1]): if in_conf[i] and in_conf[i-1]: candidates.append(i) if len(candidates) == 1: return candidates[0], years[candidates[0]] if len(candidates) > 1: # 取与最终趋势方向一致的候选 trend_dir = 1 if uf[-1] > 0 else -1 for idx in candidates: if trend_dir == 1 and uf[idx] > 0 and uf[idx-1] <= 0: return idx, years[idx] if trend_dir == -1 and uf[idx] < 0 and uf[idx-1] >= 0: return idx, years[idx] return candidates[0], years[candidates[0]] return None, None

逻辑说明：np.sign的比较本质上是在检查两条曲线的相对位置是否发生翻转。in_conf确保交叉点两端的统计量都在临界线内，过滤掉置信区间外的交叉。trend_dir由 UF 终点方向决定，这与你的气候假设有关——我在实际使用中遇到过春季降水序列 UF 终点稍微大于 0 但整条曲线来回穿越，此时直接取第一个 UF 变号位置可信度不高，建议配合趋势显著性（p 值）来判断。

参数说明：conf_level=1.96是硬编码的 95% 置信度。如果样本量很小（小于 20），我建议把置信水平放宽到 90% 即 1.645，这样能避免因为临界线过宽而把所有交点都滤掉。

4.2 交叉点的可视化呈现与年份标注

这套代码的绘图部分用 matplotlib 实现，核心是双曲线加置信线，以及突变点的散点标注。我改过一版更适合论文使用的画法，效果很直观：

def plot_mk_seasonal(results, season_name, output_path): fig, ax = plt.subplots(figsize=(8, 5)) years = results['years'] uf = results['uf'] ub = results['ub'] ax.plot(years, uf, label='UF', color='#1f77b4', linewidth=1.8) ax.plot(years, ub, label='UB', color='#ff7f0e', linewidth=1.8) ax.axhline(1.96, linestyle='--', color='gray', linewidth=0.8) ax.axhline(-1.96, linestyle='--', color='gray', linewidth=0.8) ax.axhline(0, linestyle=':', color='black', linewidth=0.6) # 标注突变点年份 if results['mutation_year'] is not None: ax.axvline(results['mutation_year'], color='red', linestyle='-.', linewidth=0.9) ax.annotate(f'突变年份: {results["mutation_year"]}', xy=(results['mutation_year'], ax.get_ylim()[0]), xytext=(results['mutation_year'], ax.get_ylim()[0] + 0.3), rotation=90, fontsize=9, color='red') ax.set_title(f'{season_name}季 M-K 突变检测', fontsize=13) ax.legend(loc='best') fig.tight_layout() fig.savefig(output_path, dpi=150) plt.close(fig)

这个绘图函数有几个细节值得说明。一是红色虚竖线标注的突变年份应该从results['mutation_year']取，而不是从交点索引直接取，因为年份可能不连续（如数据有缺失年份时）；二是画图前必须确保 UF 和 UB 的年份轴对齐，否则曲线错位会误导判读；三是保存时用dpi=150，期刊插图一般要求 300 dpi，投稿前记得改。

我实际用的版本里还加了一个功能：把四个季节的图拼成 2x2 的大图，共享一个图例。这个做法在汇报场景下很实用，审阅人一眼就能看到四个季节突变年份的差异。拼图时注意每个子图单独设置标题，颜色统一按季节分配——我习惯春季用绿色、夏季用红色、秋季用橙色、冬季用蓝色，跟季节属性匹配，方便非专业读者快速关联。

4.3 显著性判定与结果输出格式

代码最终输出的结果除了图，还有一张结构化表格，字段包括季节、突变年份、UF 终值、是否显著、趋势方向。这个输出我改造过一版，增加了置信区间的上下界，方便直接写进论文的方法部分：

这里有一个方法论上的坑，我当时在写某高校的项目报告时被导师指出过：M-K 突变检验的结果本质上是筛选与定位，不具有假设检验的严格 p 值基础。UF/UB 的显著性判断依赖于每一步累积统计量的独立检验，序列内部存在自相关时，方差公式的假设会被打破。因此，季节尺度的 M-K 结果适合作为突变年份的初筛证据，但要在论文里拿到硬结论，还需要配合其他方法验证，这点我在第 6 章专门展开。

5. 季节尺度 M-K 的高频翻车点排查

5.1 现象：UB 曲线走势与 UF 完全相同

有次我拿一个 1961-2020 年降水序列跑代码，输出的春季 UB 曲线跟 UF 曲线几乎重叠，季节突变年份一下子冒出四五个。排查后发现原因在数据组织：读取 CSV 时月份列被 pandas 读成了字符串，排序时字符串字典序导致 1 月在 10 月前面，季节分组完全错位，序对比较落在了错误的元素对之间。解决方法是强制转换类型并检验月份唯一性：

raw_data = raw_data.astype({'month': int, 'year': int})

从那以后我每次读取数据都会跑一遍data.dtypes检查，不再信任 CSV 里看似正常的数字列。

5.2 现象：所有季节突变点都集中在序列首段

冬季序列的 UF 曲线在第一二个点就冲到 2.5 以上，交点全部被过滤，输出“无突变”。原因很清楚：冬季样本量少且极端值集中，前几个点的秩比较受一个极端值主导，方差还没累积起来，UF 值虚高。解决策略是增加一个最小计算窗口。我把compute_uf_subseries里的循环起始值从 2 改成 6，即前 5 个点只做累积不参与突变判读，这样首段虚高的 UF 值不会生成候选点。注意这只是工程回避，如果极端值是真实物理过程，要保留在序列里，不要为了平滑而删除。

5.3 现象：交点年份与已知气候突变年份不一致

某流域夏季降水序列，代码给出的突变年份是 1994，但公开发表文献一致指向 1992 年区域环流调整。检查发现，问题出在插补策略——我把 1992 年上半年连续的缺失值用线性插值填充，插值造成序列中间段的秩比较信息失真。M-K 检验对秩敏感，插值序列与真实观测的秩关系不可靠。解决方式是改用季节均值填充。对降水这种强季节序列，用该季节多年气候平均值替换缺失值，对秩结构的影响远小于线性插值。我改成了：

sub_series = sub_series.fillna(sub_series.mean())

这方法不完美，但至少不会引入趋势性偏差。

5.4 现象：代码运行无报错，但输出结果明显偏离理论预期

冬季 M-K 结果显著上升，而我人工绘制距平曲线看到的却是明显下降。此事一度让我怀疑算法实现有问题，后来发现是年份排序方向问题。CSV 中数据按年份降序排列，代码内部做sort_values('year')时默认升序，但冬季跨年份的归属处理在排序前执行，导致 2020 年 12 月与 2020 年 1 月被归入同一季节，从而把最末端的值错误地带入了序列开头。解决方法是严格规定数据格式：原始文件必须按年升序保存，若使用水文年季节划分，则先处理跨年归属再排序。这个坑我花了整整一个下午定位，最后打了一行print(season_years[:10])才发现排序错乱。

5.5 现象：不同季节的突变年份彼此矛盾，无法解释

四个季节最终输出：春季 1994、夏季 2003、秋季 无、冬季 1998。这种情况在真实数据中非常常见，不代表代码有 bug，而是气候系统的多尺度响应本身存在季节差异。我处理方式是不强行统一，而是把四个季节的突变年份连同相应季节的趋势显著性写成对比矩阵。若论文要求给出一个总突变年份，可以采用季节面积加权或最主要季节降水贡献率加权，但必须在方法部分说明加权依据。

6. 交叉验证技巧：用 Pettitt 检验与滑动 t 检验增强突变结论可信度

季节尺度 M-K 突变点往往不止一个候选位置，特别是当 UF 与 UB 在临界线附近反复缠绕时，算法选出的交点年份需要外部验证。我常用的验证手法是 Pettitt 检验，它本质上是非参数的变点检测，对单一变点的定位能力比 M-K 的目视判读更严格，适合作为 M-K 结果的交叉印证。

Pettitt 检验代码简短，我在 M-K 跑出候选突变年份后，会再跑一遍这个：

def pettitt_test(data): n = len(data) ranks = pd.Series(data).rank().values U = 0.0 max_U = 0.0 k_star = 0 for k in range(n - 1): # 计算候选变点前后累积等级差 U += np.sum(ranks[:k+1]) - k * (k + 1) / 2.0 if abs(U) > abs(max_U): max_U = U k_star = k # 显著性近似公式 p_value = 2.0 * math.exp(-6.0 * max_U**2 / (n**3 + n**2)) return k_star, p_value

逻辑说明：这个实现利用了一种递推关系，将原 Pettitt 的双重循环化简为单次遍历。U 值绝对值最大点对应的 k_star 就是变点位置。p 值的公式是渐近近似，当 n 超过 50 时可用性较好，n 在 20 到 50 之间时我通常用蒙特卡洛模拟获取更可靠 p 值。

参数说明：data是某一季节的年值序列。注意输入顺序必须与 M-K 分析时的季节子序列一致，否则交叉验证就失效了。对应关系上，Pettitt 给出的 k_star 是基于 0 的索引，对应原文中第 k_star+1 个样本年份，与 M-K 输出的年份相差一年以内都算吻合。

我个人的工作习惯是：对于 M-K 疑难的、UF 与 UB 多次交叉的季节序列，用滑动 t 检验做第二轮验证。滑动 t 检验的做法是设定窗口长度（常见取 5 或 7 年），从序列第三个窗口开始逐步平移，计算窗口内前后两段均值的 t 统计量。t 值超过临界值的年份，如果与 M-K 交点差在两年内，我就把该交点视为可靠的突变年份并在报告中注明。

滑动 t 检验适合前置判断：如果滑动 t 检验全程无显著 t 值，即使 M-K 给出了清晰交点，也建议对该突变结论保持保守表述。反过来，如果滑动 t 检验显著，但 M-K 交点找不到，往往说明趋势是渐变而非突变，应该考虑修改分析框架。

最后说一个血泪经验：M-K 突变检测的结果尤其依赖序列端点。如果你要对近几年的变化趋势做判断，注意序列末端数据往往处于累积窗口尚未成熟的阶段，末端交点看可以，但不能作为严格论文结论。我在处理某地区降水序列时，M-K 检测显示 2021 年有一个明显上升突变，加长序列到 2023 年后发现这个突变消失了——纯属末端数据不足导致的假象。从那以后，我每次跑完分析都会强制做一遍尾端截断验证：把序列截去最后 3 年再跑，如果突变点位置发生超过 2 年的漂移，那这个突变点我宁可不上报，也不要冒险写进结论。希望这个验证习惯也能帮到你。

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