简介面向气象、水文与环境科学领域的学习者与科研人员这份材料给出一套用Python实现季节尺度Mann-Kendall突变检测的完整脚本用于在具有周期性波动的长序列中判断趋势走向并定位突变发生的时点读者需具备pandas、numpy的基本操作能力。压缩包共2个文件含1个py脚本与1个xlsx示例数据表整体约11KB轻量易读解压即可查看与运行。脚本以SPEI标准化降水蒸散发指数为样例序列串联数据读取与缺失值核查、季节成分与趋势成分分解、对趋势项开展M-K检验并输出Z值与P值、突变点定位及结果绘图等环节形成从原始数据到判读结论的完整链路。目前已有200人学习。读者可将示例数据替换为自身观测序列直接复用其中的分解流程、检验代码与绘图模板省去从零搭建季节性突变检测框架的重复工作也便于对照输出结果理解方法细节。1. 季节尺度 M-K 突变检测一份 .py 能替你省掉多少手工分组一份 30 年的月尺度监测序列直接丢进普通 M-KZ 值常常大得吓人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)(2m5)/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 个月这个别照抄。如果数据本身是季度值就直接用季度列不要再按月拆。站点合并之前先做齐性检验不同站点的突变年可能本来就不一样合并会抹掉差异。检查项经验阈值不满足时怎么办每个季节样本量至少 10最好 20少于 10 的季节只做趋势描述不下突变结论连续缺测不超过 3 年缺测段前后分开跑不要硬插值lag-1 自相关绝对值小于 0.2考虑预白化或 Hamed-Rao 方差修正lag-12 自相关明显高于 lag-1优先逐季节或去季节化别直接普通 M-K相同值比例低于 30%检查记录精度ties 修正要写进代码序列长度20 年以上少于 20 年只做探索性分析提示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用来放图和数据表。参数含义常用值--input输入 CSV 路径月尺度长表--date-col日期列名date--value-col数值列名value--station-col站点列名可选station--alpha显著性水平0.05--min-n季节内最小样本量10 或 15--outdir输出目录output--modeseason 或 deseasonseason脚本内部一般分四层读表和清洗、季节分组与 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_colvalue, season_colseason, alpha0.05, min_n10): 按季节分组计算 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(dtypefloat) 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_countsTrue) 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 和总方差算 Zp 取双侧。参数说明value_col是数值列season_col是季节列alpha只影响显著性标记不改 Zmin_n低于阈值的季节会被跳过明细表里也不出现。这个函数只回答趋势不回答突变别拿它的 p 去标突变年。3.3 核心函数二UF/UB 突变检测与 Sen 斜率def mk_mutation(x, alpha0.05): M-K 突变检验正序 UF、逆序 UB返回交叉点索引。 交叉点条件落在置信线内且相邻符号发生变化。 x np.asarray(x, dtypefloat) 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, tNone): Sen 斜率所有点对斜率的中位数单位是 值/年。 x np.asarray(x, dtypefloat) 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决定zcrit0.05 对应 1.96x是单个季节子序列长度不足时不要调用sen_slope的t传年份数组单位就是“每年变化量”不传就按索引差算。注意 UB 对齐容易写错跑完先用一条已知上升序列验证 UF 是否为正、UB 是否从负转正。3.4 主流程逐季节跑与去季节化跑def run_seasonal(df, value_colvalue, season_colseason, year_colyear, alpha0.05, min_n10): rows [] for season, g in df.groupby(season_col): g g.sort_values(year_col) x g[value_col].to_numpy(dtypefloat) years g[year_col].to_numpy() if len(x) min_n: continue uf, ub, cross mk_mutation(x, alphaalpha) 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 年序列可以设 15alpha与前面的趋势检验保持一致不要一个 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 个靠近端点的交点直接降级为“疑似”。这三条同时满足才能把交点年份列为突变候选。只满足第一条可能是噪声只满足第二条可能是趋势显著但没有突变三条都满足再去找第二意见。现象可能解释处理置信线内交叉一次并持续分离候选突变记录年份做 Pettitt 或滑动 t 验证多次交叉、来回穿周期残留或序列非单调检查 lag-12 自相关改逐季节或去季节化交点在序列两端端点效应不报精确年份只报方向UF 超过 1.96 但没有交叉显著趋势不是突变只写趋势不写突变年UF 和 UB 都在 0 附近无明显趋势检查数据长度和量纲注意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, dtypefloat) 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, labelUF, colorC0) ax.plot(years, ub, labelUB, colorC1) ax.axhline(1.96, ls--, lw0.8, colorgray) ax.axhline(-1.96, ls--, lw0.8, colorgray) ax.axhline(0, lw0.6, colorblack) ax.set_xlabel(Year) ax.set_ylabel(Standardized statistic) ax.legend() fig.tight_layout() fig.savefig(ufub.png, dpi300)逻辑说明UF 和 UB 画在同一坐标系置信线用虚线零线用细实线。参数说明dpi300够报告用图里用英文标签避免服务器缺中文字体导致方块。matplotlib 中文字体在无桌面环境翻车是常事能避开就避开。结果表建议至少八个字段季节、样本量、趋势 Z、趋势 p、Sen 斜率、突变年、交叉次数、突变前后均值差。逐季节跑完拼成一张长表to_csv输出后面写报告直接筛。注意趋势 p 和突变年是两列不要合并成一个“显著性”字段否则读者会以为 p 小就是突变显著。5. 避坑与排查五个让突变点漂移的常见问题5.1 现象p 值小得离谱UF/UB 满屏交叉这是最常见的翻车现场。30 年月数据直接做普通 M-Kn 有 360Z 值动不动上 5p 值小到 1e-8UF/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 算出 Z2.8B 算出 Z2.3先别怀疑数据。原因可能是 ties 修正有没有加、方差公式用的是哪种、有没有预白化、季节分组方式不同、缺失值处理不同、p 值是单侧还是双侧。解决固定一套实现输出中间量 S、Var、ties 数用一条已知序列做单元测试比如 1 到 20 的严格上升序列Z 应该为正且很大不要本文还有配套的精品资源点击获取 
长河君 · 2026/11/2 6:46:11 · 阅读 1.2k 季节尺度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-KZ 值常常大得吓人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)(2m5)/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 个月这个别照抄。如果数据本身是季度值就直接用季度列不要再按月拆。站点合并之前先做齐性检验不同站点的突变年可能本来就不一样合并会抹掉差异。检查项经验阈值不满足时怎么办每个季节样本量至少 10最好 20少于 10 的季节只做趋势描述不下突变结论连续缺测不超过 3 年缺测段前后分开跑不要硬插值lag-1 自相关绝对值小于 0.2考虑预白化或 Hamed-Rao 方差修正lag-12 自相关明显高于 lag-1优先逐季节或去季节化别直接普通 M-K相同值比例低于 30%检查记录精度ties 修正要写进代码序列长度20 年以上少于 20 年只做探索性分析提示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用来放图和数据表。参数含义常用值--input输入 CSV 路径月尺度长表--date-col日期列名date--value-col数值列名value--station-col站点列名可选station--alpha显著性水平0.05--min-n季节内最小样本量10 或 15--outdir输出目录output--modeseason 或 deseasonseason脚本内部一般分四层读表和清洗、季节分组与 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_colvalue, season_colseason, alpha0.05, min_n10): 按季节分组计算 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(dtypefloat) 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_countsTrue) 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 和总方差算 Zp 取双侧。参数说明value_col是数值列season_col是季节列alpha只影响显著性标记不改 Zmin_n低于阈值的季节会被跳过明细表里也不出现。这个函数只回答趋势不回答突变别拿它的 p 去标突变年。3.3 核心函数二UF/UB 突变检测与 Sen 斜率def mk_mutation(x, alpha0.05): M-K 突变检验正序 UF、逆序 UB返回交叉点索引。 交叉点条件落在置信线内且相邻符号发生变化。 x np.asarray(x, dtypefloat) 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, tNone): Sen 斜率所有点对斜率的中位数单位是 值/年。 x np.asarray(x, dtypefloat) 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决定zcrit0.05 对应 1.96x是单个季节子序列长度不足时不要调用sen_slope的t传年份数组单位就是“每年变化量”不传就按索引差算。注意 UB 对齐容易写错跑完先用一条已知上升序列验证 UF 是否为正、UB 是否从负转正。3.4 主流程逐季节跑与去季节化跑def run_seasonal(df, value_colvalue, season_colseason, year_colyear, alpha0.05, min_n10): rows [] for season, g in df.groupby(season_col): g g.sort_values(year_col) x g[value_col].to_numpy(dtypefloat) years g[year_col].to_numpy() if len(x) min_n: continue uf, ub, cross mk_mutation(x, alphaalpha) 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 年序列可以设 15alpha与前面的趋势检验保持一致不要一个 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 个靠近端点的交点直接降级为“疑似”。这三条同时满足才能把交点年份列为突变候选。只满足第一条可能是噪声只满足第二条可能是趋势显著但没有突变三条都满足再去找第二意见。现象可能解释处理置信线内交叉一次并持续分离候选突变记录年份做 Pettitt 或滑动 t 验证多次交叉、来回穿周期残留或序列非单调检查 lag-12 自相关改逐季节或去季节化交点在序列两端端点效应不报精确年份只报方向UF 超过 1.96 但没有交叉显著趋势不是突变只写趋势不写突变年UF 和 UB 都在 0 附近无明显趋势检查数据长度和量纲注意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, dtypefloat) 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, labelUF, colorC0) ax.plot(years, ub, labelUB, colorC1) ax.axhline(1.96, ls--, lw0.8, colorgray) ax.axhline(-1.96, ls--, lw0.8, colorgray) ax.axhline(0, lw0.6, colorblack) ax.set_xlabel(Year) ax.set_ylabel(Standardized statistic) ax.legend() fig.tight_layout() fig.savefig(ufub.png, dpi300)逻辑说明UF 和 UB 画在同一坐标系置信线用虚线零线用细实线。参数说明dpi300够报告用图里用英文标签避免服务器缺中文字体导致方块。matplotlib 中文字体在无桌面环境翻车是常事能避开就避开。结果表建议至少八个字段季节、样本量、趋势 Z、趋势 p、Sen 斜率、突变年、交叉次数、突变前后均值差。逐季节跑完拼成一张长表to_csv输出后面写报告直接筛。注意趋势 p 和突变年是两列不要合并成一个“显著性”字段否则读者会以为 p 小就是突变显著。5. 避坑与排查五个让突变点漂移的常见问题5.1 现象p 值小得离谱UF/UB 满屏交叉这是最常见的翻车现场。30 年月数据直接做普通 M-Kn 有 360Z 值动不动上 5p 值小到 1e-8UF/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 算出 Z2.8B 算出 Z2.3先别怀疑数据。原因可能是 ties 修正有没有加、方差公式用的是哪种、有没有预白化、季节分组方式不同、缺失值处理不同、p 值是单侧还是双侧。解决固定一套实现输出中间量 S、Var、ties 数用一条已知序列做单元测试比如 1 到 20 的严格上升序列Z 应该为正且很大不要本文还有配套的精品资源点击获取