季节尺度M-K突变检测：Seasonal Kendall与UF/UB实现

简介：面向气象、水文与环境科学领域的学习者与科研人员，这份材料给出一套用Python实现季节尺度Mann-Kendall突变检测的完整脚本，用于在具有周期性波动的长序列中判断趋势走向并定位突变发生的时点，读者需具备pandas、numpy的基本操作能力。压缩包共2个文件，含1个py脚本与1个xlsx示例数据表，整体约11KB，轻量易读，解压即可查看与运行。脚本以SPEI标准化降水蒸散发指数为样例序列，串联数据读取与缺失值核查、季节成分与趋势成分分解、对趋势项开展M-K检验并输出Z值与P值、突变点定位及结果绘图等环节，形成从原始数据到判读结论的完整链路。目前已有200人学习。读者可将示例数据替换为自身观测序列，直接复用其中的分解流程、检验代码与绘图模板，省去从零搭建季节性突变检测框架的重复工作，也便于对照输出结果理解方法细节。

1. 季节尺度 M-K 突变检测：一份 .py 能替你省掉多少手工分组

一份 30 年的月尺度监测序列，直接丢进普通 M-K，Z 值常常大得吓人，p 值小到小数点后一串零；可把数据按季节拆开，冬季不显著、夏季反而下降。这不是数据错了，是 12 个月的周期把方差撑大了。季节尺度 M-K 突变检测用 Python 把“按月/季分组—季节 Kendall 趋势—UF/UB 突变点—Sen 斜率”串成一条可复现流水线，解决的就是“整体显著、分季翻车”这类问题。它适合手里有 20 年以上月/季序列、需要判断突变年份和季节差异的人，也适合被普通 M-K 的假显著坑过、想补一套交叉验证的人。边界后面讲，先把它跑通。

2. 为什么按月分组的 M-K 不能直接套：季节 Kendall、去季节化与突变检验的边界

2.1 普通 M-K 的假设与季节序列的冲突：方差被低估之后

Mann-Kendall 检验是非参数方法，核心统计量 S 由所有“后值减前值”的符号和构成，不要求数据正态，对异常值也不算敏感，所以在水文、气象、环境监测里被用得很泛。它的原假设是“序列无单调趋势”，备择假设是“存在单调上升或下降趋势”。公式不复杂：S 等于每一对样本符号函数之和，方差里再扣掉相同值（ties）的修正项，最后标准化成 Z，查正态表得 p。问题出在“样本独立”这个前提上。月尺度序列天然带 12 个月周期，lag-12 自相关很高，相邻年份同月的数据也往往相关。样本不独立时，S 的方差会被低估，Z 被放大，p 值偏小，于是本该谨慎的结论变成了“极显著”。

更隐蔽的是，普通 M-K 只回答“有没有单调趋势”，不回答“哪一年发生了突变”。很多人把 UF/UB 突变检验和趋势检验混在一起，看到 UF 超过 1.96 就宣布突变，这是两码事。趋势检验看整段序列的单调性，突变检验看序列内部均值或分布是否在某个时间点发生改变。季节尺度上，这两件事还要再拆一层：你是关心“每个季节内部有没有突变”，还是关心“去掉季节周期后的整体序列有没有突变”？路线不同，代码和解释都不同。

Seasonal Kendall 检验的做法是按季节分组，分别算每个季节的 S_g 和 Var(S_g)，再求和得到总 S 和总方差，最后做一次 Z 检验。它不要求季节之间独立，只要求每个季节子序列内部近似独立，实践上比直接对月序列做 M-K 稳得多。代价是每个季节的样本量不能太少，通常至少 10 个点，最好 20 个点以上。另一个路线是去季节化：对每个观测值减去对应季节的均值或中位数，把周期成分拿掉，再对残差做普通 M-K 或 UF/UB。两条路线没有绝对优劣：逐季节适合回答“春季突变年是不是比秋季早”，去季节化适合回答“整个系统在哪个年份发生了整体跃变”。我一般会两条都跑，结果方向一致才敢写进报告。

还有一点容易被忽略：Seasonal Kendall 是趋势检验，不是突变检验。它给你一个总 Z 和一个 p，告诉你整体趋势显不显著，但不会告诉你突变点。很多库里的seasonal_kendall返回的就是趋势结果，拿它当突变检测用，相当于把体温计当 CT 机。突变点还得靠 UF/UB、Pettitt 或滑动 t 去定位。把这条边界先划清，后面代码才不会用错函数。

2.2 UF/UB 突变检验的构造：顺序统计量、逆序与交叉点

UF/UB 的构造思路是“正序算一遍，逆序再算一遍”。对序列前 m 个点，计算正序统计量 S_m，减去期望 E(S_m)=m(m-1)/4，再除以标准差 sqrt(Var)，得到 UF_m。Var 的公式是 m(m-1)(2m+5)/72，有相同值时要做 ties 修正。UF 大于 0 表示上升趋势，小于 0 表示下降趋势。把序列反过来，用同样的方法算一遍，再把结果取负、按原顺序对齐，得到 UB。UB 可以理解成“从末尾往前看”的趋势。把 UF 和 UB 画在同一张图上，如果两条线在置信线之间相交，且相交后持续分离，这个交点就是突变候选点。

置信线通常取 ±1.96，对应双侧 0.05 的显著性水平。注意三个条件要同时满足：第一，交点落在 ±1.96 之间；第二，交点之后 UF 和 UB 要持续分离，不能来回穿；第三，交点不能太靠近序列两端，前后至少各留 5 个点，否则端点效应会把结果带偏。只满足“相交”而不看置信线，是常见的翻车点。交点落在置信线外，只能说明趋势显著，不能说明突变。UF 在交点前已经超过 1.96，交点后继续上升，说明趋势在交点前就开始了，突变年可能比交点更早，这时候要结合 Pettitt 或滑动 t 看。

季节尺度上，UF/UB 有两种跑法。第一种是逐季节：把春季、夏季、秋季、冬季各自拉出来，每个子序列单独算 UF/UB，横轴是年份。这样能看出“哪个季节先突变”。第二种是去季节化后跑整体：先减掉季节均值，再对残差算 UF/UB，横轴还是年份。两种跑法的交叉点含义不同，报告里要写清楚。逐季节跑的时候，每个季节的样本量是年数，30 年就只有 30 个点，UF/UB 在两端波动很大，交点年份的精度大概只能到“前后 1 到 2 年”，别写成精确到某一年。血泪经验是：UF/UB 图好看，但它是累积统计量，对渐变响应有滞后，不能当黑匣子直接吐年份。

2.3 前置检查清单：表结构、季节映射与最小样本量

在写代码之前，先把输入表压成三列：时间、值、站点（可选）。时间列要能解析成日期，至少能拆出年和月；值列必须是数值，缺失值用空值表示，不要用 -999 之类的哨兵值混进去。季节映射按北半球常用规则：3 到 5 月为春，6 到 8 月为夏，9 到 11 月为秋，12 到次年 2 月为冬。南半球要平移 6 个月，这个别照抄。如果数据本身是季度值，就直接用季度列，不要再按月拆。站点合并之前先做齐性检验，不同站点的突变年可能本来就不一样，合并会抹掉差异。

提示：min_n 不是拍脑袋定的，先把每个季节的样本量打印出来再定。样本量不够的季节，UF/UB 的交叉点会像玄学一样漂。

缺失值处理要克制。季节内缺一个月，不要用全年均值去填，那会把季节信号抹平；更不要用前后年份的平均值直接补，容易制造假突变。常见做法是保留缺失，让每个季节的子序列长度不一致，M-K 本身能处理不等长样本，只要样本量够。如果缺测集中在某几年，建议把序列切成两段分别跑，再比较两段的突变点。数据分辨率也要统一，月数据不要和年数据混在一个表里跑，季节映射会错位。

3. 把 .py 跑起来：输入表结构、命令行参数与四个核心函数

3.1 环境、文件结构与入口参数

这份 .py 是单文件脚本，拿到后先建虚拟环境，把依赖装齐。依赖不复杂：numpy 做数组和符号运算，pandas 读表和分组，scipy 提供正态分布和统计检验，matplotlib 出图，pymannkendall 可以用来交叉核对趋势结果，但突变检测部分建议用脚本里的手写实现，方便输出中间量。下面用mk_seasonal.py代指包里的脚本文件，实际文件名以项目包为准。

python -m venv .venv source .venv/bin/activate # Windows 用 .venv\Scripts\activate pip install numpy pandas scipy matplotlib pymannkendall # 最小可用命令：输入月尺度 CSV，按季节逐段跑 python mk_seasonal.py \ --input data/monthly.csv \ --date-col date \ --value-col value \ --station-col station \ --alpha 0.05 \ --min-n 10 \ --outdir output

逻辑说明：先隔离环境，避免和系统里的旧版 numpy 打架；然后一次装齐依赖；最后用最小参数跑通。--date-col和--value-col是必填，--station-col可选，不填就当成单站点；--alpha控制置信线，默认 0.05 对应 ±1.96；--min-n是每个季节的最小子样本量，低于这个值的季节直接跳过；--outdir用来放图和数据表。

脚本内部一般分四层：读表和清洗、季节分组与 Seasonal Kendall、UF/UB 突变检测、Sen 斜率与结果汇总。读表层只做类型转换和缺失标记，不做插值；分组层负责把月份映射到季节；检测层对每个季节子序列算 UF/UB 和交叉点；汇总层把趋势、突变年、前后均值差拼成一张长表。这样拆的好处是，哪一层出问题都能单独打印中间量，不用整段重跑。

3.2 核心函数一：季节分组与 Seasonal Kendall 趋势检验

import numpy as np import pandas as pd from scipy.stats import norm def seasonal_kendall(df, value_col="value", season_col="season", alpha=0.05, min_n=10): """按季节分组计算 Seasonal Kendall 趋势检验。 返回总 Z、p 和每个季节的 S/Var 明细。""" s_total = 0.0 var_total = 0.0 detail = [] for season, g in df.groupby(season_col): x = g[value_col].dropna().to_numpy(dtype=float) n = len(x) if n < min_n: continue s = 0.0 for i in range(n - 1): for j in range(i + 1, n): s += np.sign(x[j] - x[i]) # ties 修正：相同值分组后扣减方差 _, counts = np.unique(x, return_counts=True) tie_term = np.sum(counts * (counts - 1) * (2 * counts + 5)) var = (n * (n - 1) * (2 * n + 5) - tie_term) / 18.0 s_total += s var_total += var detail.append({"season": season, "n": n, "S": s, "Var": var}) if var_total <= 0: return None, pd.DataFrame(detail) # 连续性修正：S 减去 sign(S) z = (s_total - np.sign(s_total)) / np.sqrt(var_total) p = 2 * (1 - norm.cdf(abs(z))) return {"Z": z, "p": p, "significant": p < alpha}, pd.DataFrame(detail)

逻辑说明：groupby(season_col)把春、夏、秋、冬拆开，每个季节单独算 S 和 Var，再求和。tie_term处理相同值，记录精度越低、相同值越多，这一项越不能省。最后用总 S 和总方差算 Z，p 取双侧。参数说明：value_col是数值列，season_col是季节列，alpha只影响显著性标记，不改 Z；min_n低于阈值的季节会被跳过，明细表里也不出现。这个函数只回答趋势，不回答突变，别拿它的 p 去标突变年。

3.3 核心函数二：UF/UB 突变检测与 Sen 斜率

def mk_mutation(x, alpha=0.05): """M-K 突变检验：正序 UF、逆序 UB，返回交叉点索引。 交叉点条件：落在置信线内，且相邻符号发生变化。""" x = np.asarray(x, dtype=float) n = len(x) zcrit = norm.ppf(1 - alpha / 2) def _uf_sequence(seq): m = len(seq) out = np.full(m, np.nan) for k in range(1, m): # k 是当前末尾索引 m_points = k + 1 # 当前参与计算的点数 s = 0.0 for i in range(k): for j in range(i + 1, k + 1): s += np.sign(seq[j] - seq[i]) e = m_points * (m_points - 1) / 4.0 var = m_points * (m_points - 1) * (2 * m_points + 5) / 72.0 out[k] = (s - e) / np.sqrt(var) if var > 0 else 0.0 return out uf = _uf_sequence(x) ub_rev = _uf_sequence(x[::-1]) ub = np.full(n, np.nan) for k in range(n): if not np.isnan(ub_rev[k]): ub[n - 1 - k] = -ub_rev[k] cross = [] for i in range(1, n): if np.isnan(uf[i]) or np.isnan(ub[i]): continue if (uf[i - 1] - ub[i - 1]) * (uf[i] - ub[i]) <= 0 \ and abs(uf[i]) <= zcrit and abs(ub[i]) <= zcrit: cross.append(i) return uf, ub, cross def sen_slope(x, t=None): """Sen 斜率：所有点对斜率的中位数，单位是 值/年。""" x = np.asarray(x, dtype=float) n = len(x) slopes = [] for i in range(n - 1): for j in range(i + 1, n): dt = (t[j] - t[i]) if t is not None else (j - i) if dt != 0: slopes.append((x[j] - x[i]) / dt) return np.median(slopes) if slopes else np.nan

逻辑说明：UF 按正序逐点累积，UB 把序列反转后算同样的统计量，再取负并按原索引对齐。交叉点判定用了相邻两点差值乘积小于等于 0，同时要求当前点落在 ±1.96 内。参数说明：alpha决定zcrit，0.05 对应 1.96；x是单个季节子序列，长度不足时不要调用；sen_slope的t传年份数组，单位就是“每年变化量”，不传就按索引差算。注意 UB 对齐容易写错，跑完先用一条已知上升序列验证 UF 是否为正、UB 是否从负转正。

3.4 主流程：逐季节跑与去季节化跑

def run_seasonal(df, value_col="value", season_col="season", year_col="year", alpha=0.05, min_n=10): rows = [] for season, g in df.groupby(season_col): g = g.sort_values(year_col) x = g[value_col].to_numpy(dtype=float) years = g[year_col].to_numpy() if len(x) < min_n: continue uf, ub, cross = mk_mutation(x, alpha=alpha) slope = sen_slope(x, years) mut_year = years[cross[0]] if cross else None rows.append({ "season": season, "n": len(x), "sen_slope": slope, "mutation_year": mut_year, "cross_count": len(cross), }) return pd.DataFrame(rows)

逻辑说明：按季节分组后先按年份排序，保证 UF/UB 的时间方向一致；样本量不足的季节直接跳过；交叉点取第一个落在置信线内的，输出突变年；cross_count记录交叉次数，用来判断结果稳不稳。参数说明：year_col必须是数值年份，不能是字符串；min_n建议 10 起步，30 年序列可以设 15；alpha与前面的趋势检验保持一致，不要一个 0.05 一个 0.01。

去季节化模式走另一条路：先按季节减均值，得到残差列，再对残差整体跑mk_mutation。代码上只是多一步df["resid"] = df[value_col] - df.groupby(season_col)[value_col].transform("mean")，然后把resid当输入。两种模式的结果不要混在一张表里比较，逐季节的突变年是“季节内突变”，去季节化的是“整体突变”，含义不同。

4. UF/UB 曲线怎么读：显著线、交叉点与 Sen 斜率的联合判读

4.1 置信线与交叉点：三个必须同时满足的条件

拿到 UF/UB 图，先看置信线，再看交叉。第一个条件是交点落在 ±1.96 之间；第二个条件是交点之后两条线持续分离，不是刚交完又穿回去；第三个条件是交点前后各有足够多的点，至少各 5 个，靠近端点的交点直接降级为“疑似”。这三条同时满足，才能把交点年份列为突变候选。只满足第一条，可能是噪声；只满足第二条，可能是趋势显著但没有突变；三条都满足，再去找第二意见。

注意：UF/UB 交叉点只是候选突变点，不是精确年份。序列越短，交点越容易漂。

多个交叉点的时候，不要随便取最后一个。常见做法是取最早出现且后续稳定的那个；如果两个交点相隔很近，说明序列在那段时间反复调整，突变年的不确定性很大，报告里要写区间而不是单点。另一个办法是用 Pettitt 检验取最小 p 对应的突变点，和 UF/UB 交叉点对照，两个方法指向同一年，结论才硬。滑动 t 也能用，但窗口长度是人为选的，窗口一变结果就变，适合做辅助，不适合单独定案。

4.2 Sen 斜率与突变前后均值差：把统计量翻译成业务量

UF/UB 告诉你“哪一年可能变了”，Sen 斜率告诉你“每年变多少”，突变前后均值差告诉你“变了多少”。三个一起看，才能把统计结论翻译成业务语言。Sen 斜率是所有点对斜率的中位数，单位是“值/年”，对异常值不敏感；趋势检验的 p 值告诉你斜率是否显著；突变前后均值差告诉你跃变的幅度。斜率方向为正、突变后均值也上升，结论一致；斜率方向为正、突变后均值反而下降，先查异常值和缺失值，别急着写“先升后降”。

def mutation_stats(x, years, mut_year): """突变年前后均值、绝对差和相对变化。""" x = np.asarray(x, dtype=float) years = np.asarray(years) pre = x[years < mut_year] post = x[years >= mut_year] pre_mean = np.nanmean(pre) post_mean = np.nanmean(post) delta = post_mean - pre_mean return { "pre_n": len(pre), "post_n": len(post), "pre_mean": pre_mean, "post_mean": post_mean, "delta": delta, "relative_change": delta / abs(pre_mean) if pre_mean != 0 else np.nan, }

逻辑说明：以突变年为界切前后两段，分别算均值和差值。参数说明：mut_year是 UF/UB 交叉点对应的年份，归到后段；pre_n和post_n要一起输出，前后样本量差太多时均值差不可靠；relative_change在均值接近 0 时没有意义，直接置空。Sen 斜率的单位和业务量纲一致，比如“毫克/升/年”，写报告时别只写“显著上升”，要写“每年上升多少”。

4.3 出图与结果表：一张图三行代码，一张表八个字段

import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(8, 4)) ax.plot(years, uf, label="UF", color="C0") ax.plot(years, ub, label="UB", color="C1") ax.axhline(1.96, ls="--", lw=0.8, color="gray") ax.axhline(-1.96, ls="--", lw=0.8, color="gray") ax.axhline(0, lw=0.6, color="black") ax.set_xlabel("Year") ax.set_ylabel("Standardized statistic") ax.legend() fig.tight_layout() fig.savefig("ufub.png", dpi=300)

逻辑说明：UF 和 UB 画在同一坐标系，置信线用虚线，零线用细实线。参数说明：dpi=300够报告用；图里用英文标签，避免服务器缺中文字体导致方块。matplotlib 中文字体在无桌面环境翻车是常事，能避开就避开。

结果表建议至少八个字段：季节、样本量、趋势 Z、趋势 p、Sen 斜率、突变年、交叉次数、突变前后均值差。逐季节跑完拼成一张长表，to_csv输出，后面写报告直接筛。注意趋势 p 和突变年是两列，不要合并成一个“显著性”字段，否则读者会以为 p 小就是突变显著。

5. 避坑与排查：五个让突变点漂移的常见问题

5.1 现象：p 值小得离谱，UF/UB 满屏交叉

这是最常见的翻车现场。30 年月数据直接做普通 M-K，n 有 360，Z 值动不动上 5，p 值小到 1e-8，UF/UB 图上一堆交叉点。原因不是趋势有多强，而是月序列的 lag-12 自相关把方差压低了，Z 被放大。解决顺序是：先算 lag-1 和 lag-12 自相关，lag-12 明显高就不要再跑普通 M-K；改跑 Seasonal Kendall 看整体趋势，或者逐季节跑 UF/UB；如果坚持用普通 M-K，至少做预白化或 Hamed-Rao 方差修正，并且把修正前后的 Z 都打出来。预白化会削弱趋势，别过度，去掉 lag-1 就够了，不要连 lag-12 一起硬去。

5.2 现象：同一个季节的突变年每年跑出来不一样

某一年跑出来突变在 2003，换一年数据变成 2007，再跑一次又回到 2005。原因通常是季节内样本量太小、缺测年份多、排序不稳定，或者交叉点判定把平台期也算进去了。解决：min_n提到 10 或 15；缺测不要用全年均值填，保留空值让子序列变短；季节映射固定死，不要这次按气象季节、下次按水文季节；输出cross_count，交叉次数大于 2 的季节标记为“不稳定”，人工复核；用 Pettitt 做第二意见，两个方法指向同一年再报。

5.3 现象：交叉点在置信线外，却被当成突变

只看 UF 和 UB 相交，不看 ±1.96，是新手最容易犯的错。两条线在 3.5 附近交叉，说明趋势很强，但交点不在置信区间内，不能叫突变。解决：交叉点必须落在置信线内；交叉后至少连续 3 个点 UF 或 UB 超过 1.96，或者两条线持续分离；如果交点在线外，只写“趋势显著”，不写“突变年”。血泪经验是：UF/UB 图上的交叉点越多，越要克制，别把每个交叉都解释成事件。

5.4 现象：突变年份比实际晚 1-2 年

UF/UB 是累积统计量，对渐变响应有滞后，突变年后移一两年很常见。另一个原因是预处理：有人先做 5 年滑动平均再跑突变检验，滑动窗口会把突变点抹平，交点自然往后漂。解决：突变检验用原始年值，不要先平滑；检查突变年前后 3 年的原始数据，看跃变到底发生在哪一年；报告里写“突变年约在 2005 年前后”，不要写死“2005 年”；用 Pettitt 或滑动 t 交叉验证，多个方法取交集。

5.5 现象：不同库、不同人算出的 Z 值不一样

同一份数据，A 算出 Z=2.8，B 算出 Z=2.3，先别怀疑数据。原因可能是 ties 修正有没有加、方差公式用的是哪种、有没有预白化、季节分组方式不同、缺失值处理不同、p 值是单侧还是双侧。解决：固定一套实现，输出中间量 S、Var、ties 数；用一条已知序列做单元测试，比如 1 到 20 的严格上升序列，Z 应该为正且很大；不要

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