简介这份资源提供基于Python实现季节尺度M-K突变检测的完整脚本面向从事气候、水文、生态等时间序列分析的环境科研人员与数据分析学习者。M-K检验不依赖数据分布适合非正态序列的趋势与突变点识别脚本进一步将方法拓展到季节尺度可对具有周期性特征的数据做季节性分解后再检测突变弥补常规方法忽略季节波动的不足。压缩包共2个文件含1个py脚本与1个xlsx数据表整体约11KB脚本负责读取数据、预处理、季节分解、M-K检验与结果可视化数据表则提供SPEI干旱指数样例便于直接运行验证。已有195人学习下载。通过该脚本读者可掌握从数据读取、缺失值处理、STL分解到Z-score与P值判读的完整流程并借助可视化定位突变年份为干旱监测、水资源评估等场景提供可复用的分析工具与排错参考。1. 季节尺度 M-K 突变检测为什么月尺度跑不出结论换成季节就翻盘很多人第一次做突变检测都是拿一整条月序列直接丢进mannwhitneyu或者自己写的 M-K 里跑出来 p 值忽大忽小突变点飘来飘去最后只能归结为“玄学”。问题往往不在算法而在尺度月序列里混着强烈的季节循环12 个月的周期信号会把趋势和突变淹没M-K 检验的秩统计量被周期性拉偏突变点自然不稳。季节尺度 M-K 突变检测要做的就是先把月数据按季节聚合春 3-5 月、夏 6-8 月、秋 9-11 月、冬 12-次年 2 月再对每个季节的年度序列单独做 Mann-Kendall 趋势检验和突变点识别让季节内部的信号浮出来。这套方法适合做气象、水文、农业物候、能源负荷这类有明显季节周期的长序列分析Python 生态里pymannkendall、scipy、pandas足够撑起整条链路。下面按“数据怎么整 → M-K 怎么算 → 突变点怎么定 → 坑在哪”推一遍能直接抄去改自己的数据。2. 季节聚合与 M-K 检验的落地链路从月表到季节年度序列2.1 为什么必须先聚合再检验而不是逐月跑M-K 检验的原假设是序列独立同分布月尺度序列相邻月份高度自相关冬季和夏季的均值差了几十个百分点这种结构性差异会让方差估计偏大趋势检验的灵敏度被吃掉。更麻烦的是突变检测M-K 突变检验常用的是 Pettitt 检验或 Buishand 检验或者基于 M-K 统计量的滑动窗口依赖序列的平稳性假设月序列的季节循环会让统计量在每年固定月份出现伪突变。把同一季节的月份抽出来组成年度序列等于把季节循环这个干扰项消掉剩下的年际变化才是我们关心的信号。常见做法是先按季节分组求平均或求和看变量性质得到每个季节一条长度为年份数的序列再对每条序列独立做 M-K。2.2 用 pandas 做季节聚合的最小代码import pandas as pd import numpy as np # 假设 df 有两列datedatetime和 value月值 df pd.read_csv(monthly_data.csv, parse_dates[date]) df df.sort_values(date).reset_index(dropTrue) # 定义季节映射12、1、2 归为冬季注意跨年 def season_of(month): if month in (3, 4, 5): return MAM # 春季 elif month in (6, 7, 8): return JJA # 夏季 elif month in (9, 10, 11): return SON # 秋季 else: return DJF # 冬季 df[season] df[date].dt.month.map(season_of) # 关键冬季跨年年份要按“水文年”或“季节年”重算 # 这里用简单规则12 月归入下一年1、2 月归入当年 df[season_year] df[date].dt.year mask_dec df[date].dt.month 12 df.loc[mask_dec, season_year] df.loc[mask_dec, date].dt.year 1 # 按季节 季节年聚合得到每个季节每年的均值 seasonal ( df.groupby([season, season_year])[value] .mean() .reset_index() .rename(columns{value: season_mean}) ) # 检查每个季节的年份是否连续缺年要补 for s in [MAM, JJA, SON, DJF]: sub seasonal[seasonal[season] s].sort_values(season_year) full_years range(sub[season_year].min(), sub[season_year].max() 1) missing set(full_years) - set(sub[season_year]) if missing: print(f{s} 缺失年份: {sorted(missing)})这段代码的核心逻辑是三步季节映射、跨年处理、分组聚合。参数上最容易被忽略的是冬季的年份归属——12 月如果不归到下一年冬季序列会被拆成两段突变检测直接失效。season_year的规则要和你的研究惯例一致水文领域常用水文年10 月起算气象领域常用日历年选一种并在全文统一。聚合方式用mean还是sum取决于变量温度、浓度用均值降水量、径流量用求和。聚合完一定要检查年份连续性缺年会让 M-K 的秩统计量错位后面突变点定位会偏。2.3 用 pymannkendall 做趋势检验与突变点初筛import pymannkendall as mk results {} for s in [MAM, JJA, SON, DJF]: sub seasonal[seasonal[season] s].sort_values(season_year) x sub[season_mean].values # 原始 M-K 趋势检验 trend mk.original_test(x) # Pettitt 突变检验返回突变点位置 pettitt mk.pettitt_test(x, alpha0.05) results[s] { trend: trend.trend, p_value: trend.p, tau: trend.Tau, slope: trend.slope, change_point_index: pettitt.cp, change_point_year: sub[season_year].iloc[pettitt.cp] if pettitt.cp is not None else None, pettitt_p: pettitt.p } for s, r in results.items(): print(s, r)original_test返回趋势方向、p 值、Kendall 的 Tau 和 Sen 斜率pettitt_test返回突变点索引cp和对应 p 值。参数alpha控制显著性阈值默认 0.05做多季节比较时建议用 0.05 再配合 FDR 校正否则四个季节里至少一个假阳性的概率接近 18%。cp是数组索引映射回年份要用聚合后的season_year不要用原始月表的年份。如果pettitt.cp为 None说明没有检测到显著突变点这时候不要硬报一个年份直接标注“未检出”。2.4 突变点精定位滑动 t 检验与 M-K 互验Pettitt 给的是单一突变点实际序列可能有多个突变或者突变是渐进的。常见做法是用滑动 t 检验做互验在候选突变点前后各取 n 年窗口比较两段均值差异。n 一般取 5 到 10太短方差估计不稳太长会平滑掉突变。from scipy import stats def sliding_t(x, years, window8): scores [] for i in range(window, len(x) - window): left x[i - window:i] right x[i:i window] t, p stats.ttest_ind(left, right, equal_varFalse) scores.append((years[i], abs(t), p)) return pd.DataFrame(scores, columns[year, t_abs, p]) for s in [MAM, JJA, SON, DJF]: sub seasonal[seasonal[season] s].sort_values(season_year) st sliding_t(sub[season_mean].values, sub[season_year].values, window8) peak st.loc[st[t_abs].idxmax()] print(s, 滑动t峰值年:, peak[year], t, round(peak[t_abs], 2), p, round(peak[p], 4))滑动 t 的峰值年如果和 Pettitt 的突变年接近相差 1-2 年可信度就高如果差得远说明突变不显著或者序列有多个变点需要进一步用 Buishand 范围检验或贝叶斯变点模型。窗口window是敏感参数建议做 5、8、10 三档对比结论一致才下判断。3. 参数怎么设显著性、窗口、聚合方式与多重比较3.1 显著性水平与多重比较校正四个季节各做一次检验族错误率FWER会膨胀。控制 FWER 用 Bonferroni把 alpha 除以检验次数四个季节就是 0.0125控制 FDR 用 Benjamini-Hochberg适合探索性分析。实操里我一般两套都跑Bonferroni 通过的结论写进正文BH 通过的作为参考。pymannkendall不直接返回校正后 p 值需要自己收集 p 值后用statsmodels的multipletests处理。from statsmodels.stats.multitest import multipletests pvals [results[s][p_value] for s in [MAM, JJA, SON, DJF]] reject_bh, p_bh, _, _ multipletests(pvals, alpha0.05, methodfdr_bh) reject_bonf, p_bonf, _, _ multipletests(pvals, alpha0.05, methodbonferroni) print(BH:, p_bh, reject_bh) print(Bonferroni:, p_bonf, reject_bonf)3.2 滑动窗口长度的选择依据滑动 t 的窗口长度没有万能值取决于序列长度和突变尺度。经验规则窗口取序列长度的 1/10 到 1/5且不少于 5。比如 60 年序列窗口 6-12 年合理。窗口太短t 统计量对单年异常值敏感窗口太长突变点前后各 15 年1990 年的突变可能被平移到 1985。建议先跑 Pettitt 拿到候选年再在候选年附近用不同窗口验证看峰值年是否稳定。3.3 聚合方式对结果的影响同一份月数据用均值聚合和用中位数聚合突变点可能差几年。均值对极端月份敏感中位数稳健但会损失信息。温度类变量用均值降水类变量如果月内分布偏态严重可以考虑中位数或者先做对数变换再均值。另外季节划分方式也有影响气象季节3-5、6-8、9-11、12-2和天文季节春分到夏至等差几天长序列下影响不大短序列少于 30 年建议用气象季节并说明。参数常用取值影响alpha0.05校正后 0.0125控制假阳性滑动窗口序列长度 1/10~1/5≥5窗口短敏感窗口长平滑聚合方式均值 / 中位数 / 求和均值敏感中位数稳健季节划分气象季节 / 水文季节跨年处理必须统一突变检验Pettitt 滑动 t 互验单方法易漏检多突变4. 避坑与排查季节 M-K 突变检测里最容易翻车的 5 个点4.1 现象冬季突变点总是落在 12 月或 1 月明显不合理原因冬季跨年处理错误。12 月被归到当年1、2 月也归到当年导致一个冬季被拆成两段序列在年份边界出现伪突变。解决统一season_year规则12 月归下一年1、2 月归当年或者全部按水文年重算。聚合后打印每个季节的年份列表确认连续且无重复。4.2 现象pymannkendall 报错或返回的 cp 是 None原因序列里有 NaN 或者长度太短。pymannkendall对缺失值不自动处理长度少于 8 的序列 Pettitt 检验不稳定。解决聚合后先dropna()再检查每个季节的样本量少于 10 年的季节标注“样本不足结果仅供参考”。如果 cp 为 None不要强行取argmax直接报未检出。4.3 现象四个季节里三个显著一个不显著结论怎么写原因多重比较未校正或者该季节本身变异大。解决先做 BH 校正再看该季节的 Sen 斜率方向和量级。如果校正后仍不显著但斜率方向和其他季节一致可以在讨论里提“趋势一致但未达显著”不要硬说显著。另外检查该季节是否有缺年或异常值拉偏方差。4.4 现象滑动 t 峰值年和 Pettitt 突变年差 5 年以上原因序列存在多个突变点或者突变是渐进的单一突变点模型不适用。解决用 Buishand 范围检验或者贝叶斯变点检测如ruptures库做多变点识别。ruptures的 Pelt 算法可以自动定变点数量代价是参数调优麻烦。实操里先用ruptures粗筛再用滑动 t 在候选点附近精定位。import ruptures as rpt for s in [MAM, JJA, SON, DJF]: sub seasonal[seasonal[season] s].sort_values(season_year) x sub[season_mean].values.reshape(-1, 1) algo rpt.Pelt(modelrbf).fit(x) bkps algo.predict(pen3) # pen 越大变点越少 years sub[season_year].values cp_years [years[i-1] for i in bkps if i len(years)] print(s, 多变点候选年:, cp_years)pen是惩罚项控制变点数量建议从 1 到 10 试结合业务知识定。modelrbf适合均值突变如果方差也变用modell2或normal。4.5 现象换一份数据或换一个季节划分结论完全反转原因序列太短或者突变信号本身弱。解决做敏感性分析把季节划分换一种比如滑动季节、聚合方式换一种、窗口换一档看结论是否稳健。如果三种设置下结论一致可信如果反转说明数据不支持强结论如实写“突变点对参数敏感需更长序列验证”。这是血泪经验审稿人最爱问的就是敏感性提前做掉比被问倒好。5. 进阶把季节 M-K 突变检测做成可复用的批量脚本5.1 批量处理多站点或多变量的函数封装单站点跑通后下一步是批量。把聚合、检验、突变定位、校正、输出封装成一个函数输入是月表 DataFrame 和变量名输出是每个季节的结果表。这样几十个站点可以一次跑完避免复制粘贴出错。def seasonal_mk(df, value_col, date_coldate, alpha0.05, window8): df df[[date_col, value_col]].dropna().sort_values(date_col).copy() df[season] df[date_col].dt.month.map(season_of) df[season_year] df[date_col].dt.year df.loc[df[date_col].dt.month 12, season_year] 1 out [] for s in [MAM, JJA, SON, DJF]: sub (df[df[season] s] .groupby(season_year)[value_col] .mean() .reset_index() .sort_values(season_year)) if len(sub) 10: out.append({season: s, note: 样本不足}) continue x sub[value_col].values years sub[season_year].values trend mk.original_test(x) pettitt mk.pettitt_test(x, alphaalpha) st sliding_t(x, years, windowwindow) peak st.loc[st[t_abs].idxmax()] out.append({ season: s, n_years: len(sub), trend: trend.trend, p: trend.p, tau: trend.Tau, sen_slope: trend.slope, pettitt_cp_year: years[pettitt.cp] if pettitt.cp is not None else None, pettitt_p: pettitt.p, sliding_t_peak_year: peak[year], sliding_t_p: peak[p] }) res pd.DataFrame(out) # 对 p 值做 BH 校正 valid res[p].notna() if valid.sum() 1: _, p_bh, _, _ multipletests(res.loc[valid, p], alphaalpha, methodfdr_bh) res.loc[valid, p_bh] p_bh return res调用时seasonal_mk(df, temperature)就能拿到四个季节的汇总表。注意season_of函数要定义在外部或者作为参数传入批量处理不同区域时季节定义可能不同。5.2 结果验证用合成序列做已知突变点回测脚本写完不能直接信要用合成数据回测。构造一条已知突变点的序列比如前 30 年均值 10后 30 年均值 12加高斯噪声看脚本能不能把突变点定在 30 附近。np.random.seed(42) n 60 x np.concatenate([np.random.normal(10, 1, 30), np.random.normal(12, 1, 30)]) years np.arange(1961, 1961 n) syn pd.DataFrame({date: pd.to_datetime([f{y}-01-01 for y in years]), value: x}) # 注意合成数据是年序列直接调 seasonal_mk 会因月份缺失报错 # 回测时直接调 mk.pettitt_test(x) 验证突变点 print(mk.pettitt_test(x, alpha0.05).cp) # 应接近 30回测通过后再上真实数据。如果合成序列都定不准真实数据更没戏。这一步很多人跳过后面结论被质疑时拿不出证据后悔药没处买。5.3 输出与可视化突变点标注的横坐标处理画图时横坐标年份太密是常见问题matplotlib里用MaxNLocator或者手动设xticks间隔。突变点用竖虚线标注Pettitt 和滑动 t 的结果用不同颜色区分。import matplotlib.pyplot as plt from matplotlib.ticker import MaxNLocator fig, axes plt.subplots(2, 2, figsize(12, 8)) for ax, s in zip(axes.flat, [MAM, JJA, SON, DJF]): sub seasonal[seasonal[season] s].sort_values(season_year) ax.plot(sub[season_year], sub[season_mean], markero, ms3) cp_year results[s][change_point_year] if cp_year: ax.axvline(cp_year, colorred, ls--, labelfPettitt {cp_year}) ax.xaxis.set_major_locator(MaxNLocator(nbins8)) ax.set_title(s) ax.legend() plt.tight_layout() plt.savefig(seasonal_mk.png, dpi150)MaxNLocator(nbins8)控制横坐标最多 8 个刻度避免标签重叠。保存用dpi150以上投稿够用。图里同时标 Pettitt 和滑动 t 的突变年读者一眼能看出两种方法是否一致。这套流程我从单站点跑到多站点批量最大的教训是季节聚合的跨年规则和缺年检查必须在第一步做掉后面所有检验都建立在这条序列上第一步错了后面全白算。另外别迷信单一检验的 p 值Pettitt、滑动 t、Buishand 三种方法互验结论一致才写进报告。希望帮到你。本文还有配套的精品资源点击获取
