简介这是一份面向气象、水文及环境数据分析人员的Python实战脚本聚焦季节尺度下的Mann-Kendall突变检测。M-K检验不依赖数据分布适合非正态时间序列的趋势与突变点识别而季节尺度分析则需先分离季节成分与趋势成分常用于干旱监测、气候变化与水资源管理研究。资源包共2个文件包含1个py脚本与1个xlsx数据表压缩包约11KB体量轻便便于直接运行与二次修改。脚本以SPEI标准化降水蒸散发指数为示例数据流程覆盖数据读取与预处理、STL季节分解、趋势成分的M-K检验、Z-score与P-value计算以及突变点可视化帮助读者快速掌握从原始序列到突变结论的完整分析链路。目前已有195人学习下载适合具备Python基础、希望将统计检验落地到实际气候或干旱数据中的学习者参考。1. 季节尺度 M-K 突变检测一份能直接跑通的 Python 脚本做气象、水文或者农业数据分析的人大概率都遇到过这样的场景手头有几十年的月平均气温、月降水量或者径流量数据导师或甲方突然问一句「这个序列有没有明显的突变点是哪个季节先变的」。这时候如果只会做全年尺度的 Mann-Kendall 趋势检验是交不了差的——因为很多变量的突变信号恰恰藏在季节里全年平均反而把它抹平了。这份「基于 Python 实现季节尺度 M-K 突变检测」的脚本解决的就是这个问题把原始月尺度序列按季节拆开对每个季节单独跑 M-K 突变检测输出突变年份和置信区间。它适合有基础 Python 能力、手头有长时间序列数据、需要做突变诊断的从业者和研究生。脚本本身不依赖冷门库numpy、pandas、matplotlib 加 scipy 就能跑拿到手改一下读取路径就能出图。2. 先搞懂 M-K 突变检测在算什么从全年到季节的拆解逻辑2.1 Mann-Kendall 突变检测的核心是两个统计量M-K 突变检测不是玄学它的数学骨架其实很清晰。对一条时间序列它构造两个序列顺序统计量 UF 和逆序统计量 UB。UF 是从前往后算的标准化秩和UB 是从后往前算再取负。两条曲线画在同一张图上如果它们在置信区间内出现交点且交点之后 UF 突破置信线这个交点对应的时刻就是候选突变点。具体来说对序列 $x_1, x_2, ..., x_n$先算秩序列 $s_k$$$s_k \sum_{i1}^{k} r_i, \quad r_i \begin{cases} 1 x_i x_j \ 0 x_i \leq x_j \end{cases} \quad (j1,...,i)$$然后对 $s_k$ 做标准化得到 UF。UB 就是把序列倒过来做同样的操作再取负。这个过程的直觉是UF 反映「从过去到现在」的变化趋势UB 反映「从未来回看」的趋势两者交叉的地方就是趋势发生翻转的位置。提示M-K 检验对异常值不敏感也不要求数据服从正态分布这是它比线性回归突变检测更常用的原因。但它要求序列基本独立如果数据有强自相关需要先做预白化处理。2.2 为什么要按季节拆开做而不是直接做全年全年尺度的 M-K 检测有一个致命问题季节信号会互相抵消。比如某地夏季降水在 1985 年之后显著增加但冬季降水同期在减少全年总降水可能看起来毫无变化UF 和 UB 曲线根本不相交。你如果只跑全年就会得出「无突变」的错误结论。季节尺度拆分的逻辑是把 12 个月的数据按气象季节分组——春季 3-5 月、夏季 6-8 月、秋季 9-11 月、冬季 12-次年 2 月。每个季节单独构造时间序列然后各自跑一遍 M-K。这样每个季节的突变信号不会被其他季节稀释。常见做法是先把月数据按年份和季节聚合得到每个季节每年的均值或总量再对这个季节序列做检测。这里有个细节冬季跨年1980 年 12 月、1981 年 1 月、1981 年 2 月应该归到 1980 年冬季还是 1981 年冬季我一般归到后一年即 1981 年冬季这样时间轴连续不会出现首尾缺项。2.3 脚本的整体结构和技术选型这份脚本的流程分四步读取原始月尺度数据、按季节聚合、对每个季节跑 M-K、输出突变点和置信区间图。技术选型上numpy 负责数组运算和秩统计pandas 负责时间索引和分组聚合scipy 的 norm 用于计算标准化统计量的 p 值matplotlib 负责出图。为什么不用现成的 pyMannKendall 库因为它主要做趋势检验突变检测的 UF/UB 双曲线输出不够直接而且季节拆分需要自己写循环。自己实现一遍反而更可控参数想怎么调就怎么调。脚本里置信水平默认 0.05对应置信线 ±1.96这个值在代码里是一个变量改一处就能切换 0.01 或 0.1。3. 动手跑通脚本数据准备、季节聚合与 M-K 核心实现3.1 数据格式要求和读取方式脚本对输入数据的要求很宽松一个 CSV 文件至少两列一列是日期格式YYYY-MM或YYYY-MM-DD一列是数值。列名不固定在代码开头的配置区改一下就行。我一般会把配置区放在脚本最前面方便改路径和列名。import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.stats import norm # 配置区 FILE_PATH data/monthly_data.csv # 数据文件路径 DATE_COL date # 日期列名 VALUE_COL value # 数值列名 ALPHA 0.05 # 显著性水平 SEASON_MAP { # 月份到季节的映射 3: 春季, 4: 春季, 5: 春季, 6: 夏季, 7: 夏季, 8: 夏季, 9: 秋季, 10: 秋季, 11: 秋季, 12: 冬季, 1: 冬季, 2: 冬季 } # df pd.read_csv(FILE_PATH, parse_dates[DATE_COL]) df[year] df[DATE_COL].dt.year df[month] df[DATE_COL].dt.month df[season] df[month].map(SEASON_MAP) # 冬季跨年处理12月归到当年1-2月归到前一年冬季 df.loc[df[month].isin([1, 2]), year] df[year] - 1这段代码做了三件事读取 CSV 并解析日期、提取年份和月份、映射季节。最后一行是冬季跨年的关键处理——1 月和 2 月的数据年份减 1这样 1981 年 1 月就归到 1980 年冬季和 1980 年 12 月对齐。如果不做这一步冬季序列每年只有 12 月有数据1-2 月会跑到下一年去序列直接错位。3.2 按季节聚合成年序列聚合方式有两种均值和总量。气温、水位这类用均值降水、径流这类用总量。脚本里默认用均值改一个参数就能切换。AGG_METHOD mean # 可选 mean 或 sum seasonal_series {} for season in [春季, 夏季, 秋季, 冬季]: sub df[df[season] season] if AGG_METHOD mean: ts sub.groupby(year)[VALUE_COL].mean() else: ts sub.groupby(year)[VALUE_COL].sum() # 去掉缺测年份保证序列连续 ts ts.dropna() seasonal_series[season] ts print(f{season}: {len(ts)} 年, {ts.index.min()}-{ts.index.max()})聚合后每个季节得到一个以年份为索引的 Series。这里有个容易翻车的地方如果某年某季节有月份缺测均值会偏总量会偏低。脚本里用dropna()直接丢掉缺测年份但更稳妥的做法是设一个阈值比如某季节缺测月份超过 1 个就整年剔除。这个逻辑可以根据数据质量自己加。3.3 M-K 突变检测的核心函数实现这是整个脚本最核心的部分。UF 和 UB 的计算要严格按公式来尤其是秩统计的边界条件。def mk_mutation_test(series, alpha0.05): 对一条时间序列做 Mann-Kendall 突变检测 返回: years, UF, UB, 突变点年份列表 x series.values.astype(float) n len(x) years series.index.values # ---- 计算顺序统计量 UF ---- s np.zeros(n) for i in range(1, n): # 统计 x[i] 之前比它小的个数 r np.sum(x[:i] x[i]) s[i] s[i-1] r # 标准化 E n * (n - 1) / 4.0 Var n * (n - 1) * (2 * n 5) / 72.0 UF np.zeros(n) for i in range(1, n): if Var 0: UF[i] (s[i] - E) / np.sqrt(Var) # ---- 计算逆序统计量 UB ---- x_rev x[::-1] s_rev np.zeros(n) for i in range(1, n): r np.sum(x_rev[:i] x_rev[i]) s_rev[i] s_rev[i-1] r UB_rev np.zeros(n) for i in range(1, n): if Var 0: UB_rev[i] (s_rev[i] - E) / np.sqrt(Var) UB -UB_rev[::-1] # 翻转并取负 # ---- 找突变点UF 和 UB 交点且在置信区间内 ---- z norm.ppf(1 - alpha / 2) mutations [] for i in range(1, n): if (UF[i] - UB[i]) * (UF[i-1] - UB[i-1]) 0: if abs(UF[i]) z: mutations.append(years[i]) return years, UF, UB, mutations, z这段代码有几个参数需要说清楚。E是秩序列的期望Var是方差标准公式里 $E n(n-1)/4$$Var n(n-1)(2n5)/72$。z是置信线alpha0.05 时 z≈1.96。突变点的判定条件是相邻两点 UF-UB 的差乘积小于等于 0说明发生了交叉且交点处 UF 的绝对值小于 z说明在置信区间内。两个条件缺一不可只判交叉不判置信区间会把噪声交点也算进去。注意当序列存在多个等值时标准公式的方差需要做结校正。如果数据里重复值很多比如降水为 0 的月份建议加上校正项否则 UF/UB 会偏大假突变点增多。3.4 出图和结果输出出图部分把四个季节的 UF/UB 曲线画成 2×2 子图置信线用虚线标出突变点用红色竖线标注。fig, axes plt.subplots(2, 2, figsize(14, 10)) axes axes.flatten() for idx, season in enumerate([春季, 夏季, 秋季, 冬季]): ts seasonal_series[season] years, UF, UB, mutations, z mk_mutation_test(ts, ALPHA) ax axes[idx] ax.plot(years, UF, b-, labelUF) ax.plot(years, UB, g-, labelUB) ax.axhline(z, colorr, linestyle--, linewidth0.8) ax.axhline(-z, colorr, linestyle--, linewidth0.8) ax.axhline(0, colorgray, linewidth0.5) for m in mutations: ax.axvline(m, colorred, linestyle:, alpha0.7) ax.text(m, z*0.9, str(m), fontsize8, colorred) ax.set_title(f{season} M-K 突变检测) ax.set_xlabel(年份) ax.set_ylabel(统计量) ax.legend(locupper left, fontsize8) plt.tight_layout() plt.savefig(mk_seasonal_result.png, dpi150) plt.show() # 打印突变点汇总 for season in [春季, 夏季, 秋季, 冬季]: ts seasonal_series[season] _, _, _, mutations, _ mk_mutation_test(ts, ALPHA) print(f{season} 突变点: {mutations if mutations else 未检测到显著突变})出图后重点看两件事一是 UF 和 UB 是否在置信线之间有明确交点二是交点之后 UF 是否持续上升或下降并突破置信线。如果只是碰了一下又缩回去那大概率不是真突变是序列波动造成的假交点。4. 避坑与排查季节 M-K 检测里最容易翻车的五个地方4.1 现象四个季节都检测出突变且年份几乎一样原因原始月数据存在全局性的趋势或跳变比如仪器更换、站点迁移导致整条序列在某个时间点整体偏移。这种情况下每个季节都会检测到同一个突变点看起来「结果很好」实际上是数据一致性问题不是气候突变。解决先对原始月序列做整体一致性检验画一下原始数据的时间曲线看是否有明显的台阶。如果有需要做均一化校正或者分段处理不能直接拿来做季节突变分析。4.2 现象冬季序列长度比其他季节少好几年原因冬季跨年处理的边界问题。如果数据从某年 1 月开始那么第一年冬季只有 1-2 月没有前一年的 12 月如果数据到某年 12 月结束最后一年冬季只有 12 月。聚合后首尾年份的冬季值不可靠。解决在聚合后检查每个季节的年份范围把首尾不完整的年份剔除。脚本里可以在dropna()之后加一步如果某季节第一年或最后一年只有 1 个月的数据直接删掉。更简单的做法是要求输入数据至少覆盖完整的年份边界。4.3 现象UF 和 UB 曲线在置信线外相交但脚本没报突变点原因突变点判定条件里加了abs(UF[i]) z的限制。如果交点出现在置信线外说明这个交叉不显著脚本正确地忽略了它。但有些人会误以为脚本漏检。解决这是预期行为不是 bug。如果确实想看到所有交点可以把判定条件里的置信区间限制去掉但那样输出的突变点会包含大量噪声。我一般会保留这个限制同时在图上把置信线画出来让读者自己判断。4.4 现象数据里有大量重复值UF/UB 曲线异常偏高原因M-K 检验的标准方差公式假设序列无结无重复值。当数据里有很多相同的值比如降水为 0 的月份、气温取整后重复方差被低估标准化统计量偏大假突变点增多。解决加入结校正。校正后的方差公式为$$Var \frac{n(n-1)(2n5) - \sum_{t} t(t-1)(2t5)}{72}$$其中 $t$ 是每个重复值组的长度。在代码里加一个计算结校正的步骤把Var替换掉即可。这一步对降水数据尤其重要。4.5 现象脚本跑完不出图或者图是空的原因最常见的是日期列解析失败。如果 CSV 里日期格式是198001这种纯数字pd.read_csv不会自动识别为日期parse_dates会报错或产生 NaT后续分组全部为空。解决先检查df[DATE_COL].dtype如果是 object 或 int需要手动转换。对于YYYYMM格式用pd.to_datetime(df[DATE_COL].astype(str), format%Y%m)。对于YYYY-MM格式parse_dates一般能处理但如果前面有空格或特殊字符先做str.strip()。5. 进阶技巧滑动窗口突变检测与结果验证5.1 用滑动窗口找突变点的置信区间标准 M-K 只给出一个突变年份但实际上突变往往不是一年完成的而是一个 3-5 年的过渡期。我一般会加一个滑动窗口版本对每个季节序列用 10 年窗口滑动计算 UF看 UF 突破置信线的起始年份作为突变开始时间的估计。def sliding_window_mk(series, window10, alpha0.05): 滑动窗口 M-K返回每个窗口起点的 UF 值 x series.values.astype(float) n len(x) years series.index.values z norm.ppf(1 - alpha / 2) results [] for start in range(0, n - window 1): seg x[start:startwindow] m len(seg) s np.zeros(m) for i in range(1, m): s[i] s[i-1] np.sum(seg[:i] seg[i]) E m * (m - 1) / 4.0 Var m * (m - 1) * (2 * m 5) / 72.0 uf (s[-1] - E) / np.sqrt(Var) if Var 0 else 0 results.append((years[start], uf)) return results这个函数的输出是一系列(窗口起始年, UF值)把 UF 首次超过 z 的年份作为突变起始年比单点检测更稳健。窗口长度一般取 10 年序列短的话可以降到 8 年但不要低于 6 年否则统计量不稳定。5.2 用 Pettitt 检验交叉验证突变点M-K 突变检测有时会给出多个候选点这时候可以用 Pettitt 检验做交叉验证。Pettitt 检验是基于秩的非参数突变检验和 M-K 互补。如果两种方法给出的突变年份接近相差 1-2 年那这个突变点的可信度就很高。def pettitt_test(series): Pettitt 突变检验返回最可能的突变年份 x series.values.astype(float) n len(x) years series.index.values U np.zeros(n) for t in range(1, n): U[t] U[t-1] np.sum(np.sign(x[t] - x[:t])) K np.max(np.abs(U)) # 显著性判断 p 2 * np.exp(-6 * K**2 / (n**3 n**2)) idx np.argmax(np.abs(U)) return years[idx], pPettitt 的 p 值小于 0.05 时认为突变显著。把它的结果和 M-K 的结果放在一起看如果一致基本可以下结论如果不一致需要回到原始数据检查是否有缺测、异常值或者处理方式差异。5.3 结果验证的三个习惯第一个习惯画原始季节序列的折线图肉眼确认突变年份附近是否有明显的均值跳变。统计检验是辅助数据本身的形态才是最终依据。第二个习惯把突变年份和已知的气候事件或人类活动对照。比如某流域 1990 年代的水利工程建设、某地区 2000 年后的城市化加速这些都能为突变点提供物理解释。没有物理解释的突变点需要谨慎对待。第三个习惯换一种聚合方式重跑。如果均值聚合检测到突变但总量聚合没有说明突变主要体现在年内分配上而不是总量变化。这个信息本身就有价值。提示脚本里的所有参数——置信水平、聚合方式、季节划分——都放在配置区改完直接跑。建议每次改参数后把结果图另存一份方便对比。从那以后我每次做季节突变检测都会先把原始月序列画出来看一眼再跑 M-K最后用 Pettitt 交叉验证。三步走完结论才敢往报告里写。希望帮到你。本文还有配套的精品资源点击获取