简介基于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_{i1}^{k} sum_{j1}^{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)(2k5)/72。标准化之后得到 UF(k) [S(k) - k(k-1)/4] / sqrt( k(k-1)(2k5)/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_j5) / 72 seg x[:k] _, counts np.unique(seg, return_countsTrue) 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, dtypenp.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)但实际工程里季节子序列长度通常只有 3060 个样本跑起来没有任何压力。_tie_correction是容易被忽略的。当序列里有大量等值样本比如连续多个月降水量为 0 时不修正方差会把 Var(S) 算小UF 绝对值被放大趋势显著性被高估。修正项就是前面公式里的 sum t(t-1)(2t5)/72这里的 t 是每个等值组的大小。反向 UB 的构造有一个容易写错的地方先倒序计算uf_tmp得到的是“从尾巴开始往前累积”的标准曲线然后要取负并[::-1]翻转回原始时间轴。顺序不能反符号也不能丢。为什么取负因为倒着看序列时原来上升的趋势在反向视角里会变成下降标准化之后符号相反必须乘以 -1 才能让两条曲线在同一坐标系下有意义。2.3 绘制和交叉点查找NaN 开头与 ±1.96 边界得到 UF/UB 之后前两个点通常不参与展示因为 k1 和 k2 时 S(k) 太小标准化结果恒为 0画出来会引出诡异的“零段”。我一般会把 uf[0]、uf[1] 直接置为 NaN让 matplotlib 自动断开这两段。找交叉点也不能只看“两条线有没有交点”。真正可靠的突变点需要满足两个条件第一UF 与 UB 在 ±1.96 置信线之间或附近发生符号切换第二交点之后 UF 或 UB 至少有一方继续朝同方向延伸而不是反复震荡。下面这段是实际脚本里用的交叉点判断函数。def find_crossing(uf, ub, crit1.96, min_k3): 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_k3跳过开头两个无效点你也可以把它设成 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_colmonth, stylecalendar): stylecalendar: 自然历法季度 styleflood: 汛期/非汛期汛期为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, alpha0.05, stylecalendar): 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(dtypenp.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, critcrit, min_k3) results[season] { years: years, uf: uf, ub: ub, crossings: crossings } return results, crit这个函数的关键参数alpha0.05决定临界线宽度style直接传给前面拆分函数用来控制季节划分。循环内部每个季节独立排序、独立跑检验、独立找交叉点互不干扰。这是季节尺度 M-K 和年尺度 M-K 的最大差别——四条序列完全并行。拿到results之后可以按你自己的绘图风格画图。我给出一个最少参数的绘图版本方便你把脚本改成可执行文件时直接复用。def plot_seasonal_mk(results, crit, output_pathseasonal_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], labelUF, lw1.6) ax.plot(res[years], res[ub], labelUB, lw1.6) ax.axhline(crit, colorgray, linestyle--, lw1) ax.axhline(-crit, colorgray, linestyle--, lw1) for c in res[crossings]: ax.axvline(res[years][c], colorred, linestyle:, lw1) ax.set_title(fseason {s}) ax.legend() fig.tight_layout() fig.savefig(output_path, dpi300)绘图时我习惯把交叉点用竖直虚线标出来这样能直观看出四季突变是否同步。之前跑某模拟项目 X 时我印象最深的是夏季交叉点落在第 29 年秋季交叉点在同年冬季没有交叉——这说明突变不是“全年一致变化”而是季节性的年尺度 M-K 根本捕捉不到这种结构性差异。多季节同时检验隐藏着一个多重比较问题4 条序列同时跑纯随机情况下至少有一条出现假交叉的概率会高于 5%。如果数据量够可以把单个季节的 α 下调为 α/4也就是约 0.0125对应临界线 2.5。对一般业务分析我建议先按 0.05 看整体再单独关注那些交叉点前后趋势方向一致的季节不必一上来就做 Bonferroni 校正。4. 季节 M-K 避坑记录四类经典翻车与修复办法以下每一条都是我在实际处理水文序列时碰到过的状况按“现象 → 原因 → 解决”写方便你排查自己脚本输出的异常结果。4.1 交叉点总是出现在第三第四个样本附近且每次换数据都变现象季节 UF/UB 曲线在序列开头约第 35 个点位置出现明显交叉全年四个季节无一例外交叉年份落在前五年内。但把数据起始年份改一下交叉点又跟着跑到新的“开头”。原因序贯 M-K 的 S(k) 在样本量很小时方差极小标准化统计量对前几个样本的相对大小异常敏感。如果早期恰好有两年极端高值或低值UF 会被瞬间推向临界线之上与偏态反向的 UB 构成一个“虚拟交叉”。这个交叉不代表真实突变只是短序列的边界效应。解决把find_crossing的min_k提高到总长度的 10% 左右比如 40 年序列至少从第 5 年开始找。更稳妥的做法是要求交叉点出现后UF 或 UB 至少在后续连续 5 个样本保持同方向不满足就直接标记为“边缘交叉”。在脚本里加一行过滤条件len(uf[k:k5]) 5且符号方向一致即可。4.2 降水序列里大量零值UF 曲线冲上天花板现象干旱区冬季月降水常年为 0季节均值序列里有超过三分之一的值等值。运行 M-K 后 UF 绝对值经常超过 5甚至 10远高于 1.96图上曲线直接冲顶。原因等值组破坏了 S(k) 的方差结构。公式中 Var(S) 的 t(t-1)(2t5) 修正项没实现或实现有 bug导致方差被低估。零值越多低估越严重UF 被过度放大。解决确保每个 k 都调用_tie_correction且修正项参与减法的位置正确。我踩过坑之后养成的习惯是写一个自检构造长度为 30 的全零序列UF 应该为 0 或接近 0再构造一半为 0 的序列观察 UF 不会超过 3。如果全零序列 UF 异常先查 tie 修正再查np.unique在浮点精度下的等值判断是否失效。4.3 季节子序列太短只有 1520 个样本现象拆开季节后每年只有 3 个月的数据参与计算40 年数据变成 40 个点还好如果数据只覆盖 2005-2020 年每个季节只剩 16 个样本交叉点极不稳定不同 α 下结论完全相反。原因序贯 M-K 是渐近正态统计n 低于 20 时标准化公式的近似精度明显下降尤其当序列存在一两个极值时单个样本就能扭转整条 UF 曲线方向。解决一是合并相邻季节把“冬季早春”作为滑动窗口单元二是给季节子序列做滑动段合成把每 3 年合并成一段序列长度可以拉伸到 40 以上代价是突变时间分辨率下降到 3 年。前者保突变精度后者保统计稳健性按你业务对时间精度的要求选。脚本里我默认在seasonal_mk_analysis加了一个min_len20的硬性提示小于该值直接抛警告。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_perm200): uf0, ub0 mann_kendall_sequential_full(x) obs_cross len(find_crossing(uf0, ub0, critcrit)) 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, critcrit)) 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。改完直接运行就能得到四张季节突变图和对应的交叉年份表。希望帮到你。本文还有配套的精品资源点击获取 季节尺度Mann-Kendall突变检测：Python实现UF/UB与交叉点定位
简介基于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_{i1}^{k} sum_{j1}^{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)(2k5)/72。标准化之后得到 UF(k) [S(k) - k(k-1)/4] / sqrt( k(k-1)(2k5)/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_j5) / 72 seg x[:k] _, counts np.unique(seg, return_countsTrue) 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, dtypenp.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)但实际工程里季节子序列长度通常只有 3060 个样本跑起来没有任何压力。_tie_correction是容易被忽略的。当序列里有大量等值样本比如连续多个月降水量为 0 时不修正方差会把 Var(S) 算小UF 绝对值被放大趋势显著性被高估。修正项就是前面公式里的 sum t(t-1)(2t5)/72这里的 t 是每个等值组的大小。反向 UB 的构造有一个容易写错的地方先倒序计算uf_tmp得到的是“从尾巴开始往前累积”的标准曲线然后要取负并[::-1]翻转回原始时间轴。顺序不能反符号也不能丢。为什么取负因为倒着看序列时原来上升的趋势在反向视角里会变成下降标准化之后符号相反必须乘以 -1 才能让两条曲线在同一坐标系下有意义。2.3 绘制和交叉点查找NaN 开头与 ±1.96 边界得到 UF/UB 之后前两个点通常不参与展示因为 k1 和 k2 时 S(k) 太小标准化结果恒为 0画出来会引出诡异的“零段”。我一般会把 uf[0]、uf[1] 直接置为 NaN让 matplotlib 自动断开这两段。找交叉点也不能只看“两条线有没有交点”。真正可靠的突变点需要满足两个条件第一UF 与 UB 在 ±1.96 置信线之间或附近发生符号切换第二交点之后 UF 或 UB 至少有一方继续朝同方向延伸而不是反复震荡。下面这段是实际脚本里用的交叉点判断函数。def find_crossing(uf, ub, crit1.96, min_k3): 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_k3跳过开头两个无效点你也可以把它设成 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_colmonth, stylecalendar): stylecalendar: 自然历法季度 styleflood: 汛期/非汛期汛期为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, alpha0.05, stylecalendar): 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(dtypenp.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, critcrit, min_k3) results[season] { years: years, uf: uf, ub: ub, crossings: crossings } return results, crit这个函数的关键参数alpha0.05决定临界线宽度style直接传给前面拆分函数用来控制季节划分。循环内部每个季节独立排序、独立跑检验、独立找交叉点互不干扰。这是季节尺度 M-K 和年尺度 M-K 的最大差别——四条序列完全并行。拿到results之后可以按你自己的绘图风格画图。我给出一个最少参数的绘图版本方便你把脚本改成可执行文件时直接复用。def plot_seasonal_mk(results, crit, output_pathseasonal_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], labelUF, lw1.6) ax.plot(res[years], res[ub], labelUB, lw1.6) ax.axhline(crit, colorgray, linestyle--, lw1) ax.axhline(-crit, colorgray, linestyle--, lw1) for c in res[crossings]: ax.axvline(res[years][c], colorred, linestyle:, lw1) ax.set_title(fseason {s}) ax.legend() fig.tight_layout() fig.savefig(output_path, dpi300)绘图时我习惯把交叉点用竖直虚线标出来这样能直观看出四季突变是否同步。之前跑某模拟项目 X 时我印象最深的是夏季交叉点落在第 29 年秋季交叉点在同年冬季没有交叉——这说明突变不是“全年一致变化”而是季节性的年尺度 M-K 根本捕捉不到这种结构性差异。多季节同时检验隐藏着一个多重比较问题4 条序列同时跑纯随机情况下至少有一条出现假交叉的概率会高于 5%。如果数据量够可以把单个季节的 α 下调为 α/4也就是约 0.0125对应临界线 2.5。对一般业务分析我建议先按 0.05 看整体再单独关注那些交叉点前后趋势方向一致的季节不必一上来就做 Bonferroni 校正。4. 季节 M-K 避坑记录四类经典翻车与修复办法以下每一条都是我在实际处理水文序列时碰到过的状况按“现象 → 原因 → 解决”写方便你排查自己脚本输出的异常结果。4.1 交叉点总是出现在第三第四个样本附近且每次换数据都变现象季节 UF/UB 曲线在序列开头约第 35 个点位置出现明显交叉全年四个季节无一例外交叉年份落在前五年内。但把数据起始年份改一下交叉点又跟着跑到新的“开头”。原因序贯 M-K 的 S(k) 在样本量很小时方差极小标准化统计量对前几个样本的相对大小异常敏感。如果早期恰好有两年极端高值或低值UF 会被瞬间推向临界线之上与偏态反向的 UB 构成一个“虚拟交叉”。这个交叉不代表真实突变只是短序列的边界效应。解决把find_crossing的min_k提高到总长度的 10% 左右比如 40 年序列至少从第 5 年开始找。更稳妥的做法是要求交叉点出现后UF 或 UB 至少在后续连续 5 个样本保持同方向不满足就直接标记为“边缘交叉”。在脚本里加一行过滤条件len(uf[k:k5]) 5且符号方向一致即可。4.2 降水序列里大量零值UF 曲线冲上天花板现象干旱区冬季月降水常年为 0季节均值序列里有超过三分之一的值等值。运行 M-K 后 UF 绝对值经常超过 5甚至 10远高于 1.96图上曲线直接冲顶。原因等值组破坏了 S(k) 的方差结构。公式中 Var(S) 的 t(t-1)(2t5) 修正项没实现或实现有 bug导致方差被低估。零值越多低估越严重UF 被过度放大。解决确保每个 k 都调用_tie_correction且修正项参与减法的位置正确。我踩过坑之后养成的习惯是写一个自检构造长度为 30 的全零序列UF 应该为 0 或接近 0再构造一半为 0 的序列观察 UF 不会超过 3。如果全零序列 UF 异常先查 tie 修正再查np.unique在浮点精度下的等值判断是否失效。4.3 季节子序列太短只有 1520 个样本现象拆开季节后每年只有 3 个月的数据参与计算40 年数据变成 40 个点还好如果数据只覆盖 2005-2020 年每个季节只剩 16 个样本交叉点极不稳定不同 α 下结论完全相反。原因序贯 M-K 是渐近正态统计n 低于 20 时标准化公式的近似精度明显下降尤其当序列存在一两个极值时单个样本就能扭转整条 UF 曲线方向。解决一是合并相邻季节把“冬季早春”作为滑动窗口单元二是给季节子序列做滑动段合成把每 3 年合并成一段序列长度可以拉伸到 40 以上代价是突变时间分辨率下降到 3 年。前者保突变精度后者保统计稳健性按你业务对时间精度的要求选。脚本里我默认在seasonal_mk_analysis加了一个min_len20的硬性提示小于该值直接抛警告。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_perm200): uf0, ub0 mann_kendall_sequential_full(x) obs_cross len(find_crossing(uf0, ub0, critcrit)) 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, critcrit)) 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。改完直接运行就能得到四张季节突变图和对应的交叉年份表。希望帮到你。本文还有配套的精品资源点击获取