简介面向气候与环境数据分析人员及Python初学者这份基于季节尺度的Mann-Kendall突变检测脚本解决了SPEI等干旱指数中趋势突变点难以自动识别的问题。资源包共2个文件包含一个Python程序与一个SPEI3.xlsx示例数据表压缩包仅11KB整体轻量易用。脚本覆盖数据读取、缺失值检查、季节性分解STL、MK统计量计算与结果可视化等完整流程运行后可直接输出突变点年份及显著性判断。已有195人学习下载适合需要快速掌握MK检验思路、理解季节数据分解方法的读者。脚本可替换为其他时间序列数据使用也可作为论文或科研项目中干旱趋势分析的基础工具。1. 季节尺度 M-K 突变检测一条曲线定位气候突变年份很多人手里攒着一大串日尺度的气象或水文数据想找“从哪一年开始变了”上来就拿原始序列跑 Mann-Kendall 突变检验结果 UF 和 UB 两条线缠成一团交点多到没法看。真正的问题是逐日序列噪声太大MK 突变检验对高频扰动极其敏感趋势突变信号被淹没在天气尺度的波动里。把数据按季节聚合成季节均值序列再去做 M-K 突变检测交点数会大幅收敛突变年份一目了然这也是水文气候分析里更稳的做法。这份基于 Python 实现的季节尺度 M-K 突变检测脚本核心就干两件事按季节聚合序列、计算并绘制 UF/UB 统计量及置信区间帮你快速定位突变年份。适合做气候变化诊断、水文年径流分析、植被物候突变研究的从业者和研究生以及想把 MK 检验用对、用稳的 Python 用户。2. M-K 突变检验的原理与季节序列构造先搞懂 UF/UB 再动手2.1 Mann-Kendall 统计量S、方差和标准化 Z 是怎么来的Mann-Kendall 突变检验以下简称 MK其实是从趋势检验衍生出来的。先对长度为 n 的时间序列 x1, x2, ..., xn计算所有配对比较的符号和得到 S 统计量def calculate_s(seq): s 0 n len(seq) for i in range(n - 1): for j in range(i 1, n): if seq[j] seq[i]: s 1 elif seq[j] seq[i]: s - 1 return s逻辑说明S 等于正序对个数减去逆序对个数。如果序列整体上升正序对占主导S 为正的大数整体下降则 S 为负。逐对比较的时间复杂度是 O(n²)对于季节序列这种长度动辄几十年的数据完全能接受没必要优化成 O(n log n)。接下来是方差。当 n 较大通常 n10时S 近似服从正态分布方差计算公式里有一个修正项用于处理序列中存在相等值平局的情况def calculate_var_s(seq): n len(seq) # 统计每个重复值出现的次数用于平局修正 ties {} for v in seq: ties[v] ties.get(v, 0) 1 tie_term 0 for v, cnt in ties.items(): if cnt 1: tie_term cnt * (cnt - 1) * (2 * cnt 5) var_s (n * (n - 1) * (2 * n 5) - tie_term) / 18 return var_s参数说明tie_term是平局修正项数据里有重复值比如连续几年季节均值相同时必须算否则方差会被高估或低估直接影响后面标准化统计量的置信区间判断。如果序列里完全没有重复值这项为 0。有了 S 和方差就能算出标准化统计量 Z并查标准正态分布表得到显著性 p 值。但 M-K 突变检验不只是算一个 Z它需要把正向序列和反向序列的 Z 值演化过程分别画出来也就是 UF 和 UB。2.2 季节尺度的含义为什么不能拿逐日数据直接算我见过有同事把 30 年的逐日气温直接丢进 MK 突变检验画出来 UF 线像锯齿一样来回穿越置信线根本没法判断突变年份。原因在于 M-K 突变检验对序列的单调性变化敏感但逐日数据的自相关性强天气系统的短期波动会产生大量伪交点。季节尺度的思路很简单先把日数据按季节聚合成春夏秋冬四个季节的平均值或总量取决于研究变量每个季节得到一条长度为年份数的序列再分别对四条序列做 MK 突变检验。这样做有三个实际好处降噪季节平均消除了高频天气扰动保留的是年际到年代际的信号。对齐气象学定义很多突变是季节性的比如春季升温突变逐日序列根本看不出边界但春季均值序列能清楚看到转折点。满足 MK 对数据独立性的近似要求逐日数据强自相关MK 检验的前提会破季节均值相关性弱得多结果更可信。如果你研究的是日尺度事件如极端降水发生频次那就不该聚合到季节而是按年统计极端事件次数再做 MK这是另一个话题。总之别让 MK 检验的输入数据包含太多高频噪声。2.3 数据准备pandas 按季节聚合与缺失值处理写代码前先把数据格式统一。最常见的是 Pandas DataFrame列至少包含日期和数值。聚合的第一步是把日期设为索引然后映射季节import pandas as pd import numpy as np def monthly_to_seasonal(df, value_colvalue): # 要求 df 的索引是 DatetimeIndex df df.copy() df[year] df.index.year df[month] df.index.month # 气象季节划分12-2月为冬3-5月为春6-8月为夏9-11月为秋 season_map {12: DJF, 1: DJF, 2: DJF, 3: MAM, 4: MAM, 5: MAM, 6: JJA, 7: JJA, 8: JJA, 9: SON, 10: SON, 11: SON} df[season] df[month].map(season_map) # 注意冬季跨年12月应归属到下一个年份的冬季这里先不处理后面避坑章节再讲 seasonal df.groupby([year, season])[value_col].mean().unstack() return seasonal逻辑说明groupby按年和季节分组mean()取季节平均unstack()把季节变成列行是年份。这样得到的就是一个“年份 × 季节”的表格每一列就是一条待检测的季节序列。参数说明value_col指定数值列名season_map是季节映射表这里用英文缩写 DJF/MAM/JJA/SON 是为了和专业文献保持一致。需要注意冬季跨年问题1 月和 2 月属于日历年还是冬季年的问题我一般把 12 月、1 月、2 月统一划到 12 月所在年份的冬季具体取舍后面避坑章节单独说。缺失值处理上如果某个季节缺一个月直接取均值会偏。我一般要求先检查seasonal.isnull().sum()缺测月份占比超过 20% 的季节直接置为 NaN并在后续 MK 计算前剔除不要用插值硬填否则可能把突变点填出来。3. 代码实现季节尺度 M-K 突变检测的核心函数与绘图3.1 构建 UF/UB 统计量的核心函数UF正序列统计量和 UB反序列统计量的算法是一致的只是正序列是从头到尾逐步增加子序列反序列是从尾到头逐步增加子序列最后把反序列的结果取负再倒序。核心代码from scipy.stats import norm def mk_uf_ub(seq): n len(seq) uf np.zeros(n) # 正向统计量 ub np.zeros(n) # 反向统计量 # 正向计算 UF for k in range(2, n 1): sub seq[:k] s calculate_s(sub) var_s calculate_var_s(sub) if var_s 0: continue if s 0: z (s - 1) / np.sqrt(var_s) elif s 0: z (s 1) / np.sqrt(var_s) else: z 0 uf[k - 1] z # 反向计算 UB先倒序重复正向计算再取负并倒序 rev seq[::-1] ub_rev np.zeros(n) for k in range(2, n 1): sub rev[:k] s calculate_s(sub) var_s calculate_var_s(sub) if var_s 0: continue if s 0: z (s - 1) / np.sqrt(var_s) elif s 0: z (s 1) / np.sqrt(var_s) else: z 0 ub_rev[k - 1] z ub -ub_rev[::-1] uf[0] 0 ub[0] 0 return uf, ub逻辑说明calculate_s和calculate_var_s复用上一节的函数。s-1或s1是连续性修正把离散的 S 分布近似到连续正态分布。反向序列计算完成后取负是因为倒序序列的上升趋势对应原序列的下降趋势必须在符号上纠正回来。参数说明seq是一维 numpy 数组或列表建议传入 float 类型的季节均值序列长度至少 10。代码里research if var_s 0直接跳过这是处理那些子序列全相等的情况比如某几年季节均值完全相同此时 Z 无定义保持为 0。这里有个细节经常被忽略UF 和 UB 的第一位都强制置 0。因为在 k1 时子序列只有 1 个数S 统计量不存在但图上连线时需要一个起点置 0 表示从原点出发。3.2 季节序列构造与调用主流程把数据读进来、聚合、然后对每个季节跑一遍import matplotlib.pyplot as plt def run_mk_on_seasons(seasonal_df, alpha0.05): st norm.ppf(1 - alpha / 2) # 置信区间临界值 results {} for season in seasonal_df.columns: seq seasonal_df[season].dropna().values if len(seq) 10: print(f{season} 序列长度不足10跳过) continue uf, ub mk_uf_ub(seq) results[season] { uf: uf, ub: ub, seq: seq, years: seasonal_df.index.values[:len(seq)] } return results, st逻辑说明norm.ppf(1 - alpha/2)是标准正态分布的双侧分位数α0.05 时约等于 1.96。UF 和 UB 曲线在这个临界值内的交点才被认为是显著的突变点在临界值外交叉或交叉不明显的只能当作疑似突变。返回的results字典里存了每条季节序列的 UF、UB、原始序列和对应的年份数组方便绘图。参数说明alpha是显著性水平水文气候领域常用 0.05 或 0.1。样本量小时建议放宽到 0.1样本量大时用 0.01 也行。判断突变点是否显著看的是交点处的 UF 是否越过置信线而不是仅仅看两条线是否交叉。3.3 突变点可视化置信区间与交点标注绘图时把置信区间画成两条水平线UF 和 UB 分别用不同颜色交点自动标注出年份def plot_mk(results, st, save_pathNone): fig, axes plt.subplots(2, 2, figsize(14, 10)) seasons list(results.keys()) for idx, season in enumerate(seasons): ax axes[idx // 2][idx % 2] r results[season] years r[years] ax.plot(years, r[uf], labelUF, color#d62728) ax.plot(years, r[ub], labelUB, color#1f77b4) ax.axhline(st, colorgray, linestyle--, alpha0.7) ax.axhline(-st, colorgray, linestyle--, alpha0.7) # 寻找交叉点符号发生变化 cross_idx [] for i in range(len(years) - 1): if (r[uf][i] - r[ub][i]) * (r[uf][i1] - r[ub][i1]) 0: cross_idx.append(i) for ci in cross_idx: ax.axvline(years[ci], colorgreen, linestyle:, alpha0.8, linewidth1.2) ax.set_title(f{season} M-K Mutation Test) ax.legend() ax.grid(alpha0.3) ax.set_xlabel(Year) ax.set_ylabel(Z Value) fig.tight_layout() if save_path: fig.savefig(save_path, dpi300, bbox_inchestight)逻辑说明交叉点通过判断相邻两点UF - UB的乘积小于 0 来定位即两条曲线在该区间内发生了穿越。绿色虚线把可能突变年份标出来最终的人为判断还得结合置信区间只有交叉点落在 ±1.96 两条灰线之间且 UF 在交叉后继续穿越置信线才能认定为显著突变。参数说明save_path为 None 时用plt.show()显示否则保存为高分辨率 PNG。dpi300适合论文投稿bbox_inchestight防止坐标轴被截断。这里的交叉点检测是线性近似的如果序列数据点稀疏实际突变年份可能落在两个采样点之间输出结果建议写成“约 1998 年”。这套流程跑完每个季节会得到一张四象限图突变年份直接标在图上。实际项目里我会再往results里写一份 CSV把交叉点年份和是否显著存下来方便后续统计。4. 避坑与常见问题排查季节尺度 MK 最容易翻车的 5 个地方4.1 季节序列构造期的坑坑 1冬季跨年归属错位导致序列出现伪突变现象冬季DJF序列跑出来的突变年份恰好是 1990 年左右但手动看原始数据发现 1989 年 12 月异常偏暖1990 年 1、2 月正常突变纯粹是 12 月的暖冬拉高了冬季均值。原因把 12 月归到当前日历年导致冬季由 1、2 月当年和 12 月当年年底组成但气象学冬季应当跨年12 月应归入冬季年的年末即 1989 年 12 月应该和 1990 年 1、2 月组成同一个冬季。解决聚合时单独做跨年处理def assign_winter_year(df): df df.copy() # 12月归到下一年 df[winter_year] np.where(df[month] 12, df[year] 1, df[year]) # 冬季只取12、1、2月按 winter_year 聚合 winter df[df[month].isin([12, 1, 2])].groupby(winter_year)[value].mean() return winter注意这里的winter_year是冬季的年份标签这样 1989 年 12 月和 1990 年 1/2 月会一起算进 1990 年的冬季。跑完再看突变年份可能会偏移 1 年这是两种季节定义造成的正常差异。坑 2某个季节缺测严重dropna 后序列长度缩水但没提示现象夏季序列明明有 50 年跑 MK 时说序列长度不足 10 年被跳过。排查发现是因为某个夏季连续缺了两个月seasonal_df[JJA]里全是 NaNdropna()之后只剩 8 个有效年份。原因run_mk_on_seasons里seq seasonal_df[season].dropna().values把 NaN 全删了但没有记录原始长度和缺失比例。解决在聚合函数里加一个可用年份计数valid_cnt seasonal_df[season].notna().sum() if valid_cnt 10: print(f{season} 仅有 {valid_cnt} 年有效不满足最小长度要求) continue同时建议把有效年份大于 10 但小于序列总长度 60% 的季节打上“低置信”标记这类结果只能参考不宜直接下突变结论。4.2 统计量与绘图期的坑坑 3UF 和 UB 交叉了但交点不在置信区间内被误读为显著突变现象画出来的图上绿色交叉点标在 2002 年但两条线交叉那条竖线的位置UF 值只有 0.8远没到 ±1.96结果文章里被写成了“2002 年发生显著突变”审稿人直接质疑。原因判断标准不严。MK 突变检验的判定条件是UF 和 UB 在置信区间内交叉且交叉后 UF 突破置信线。只交叉不穿越只能说明序列在那个点附近出现转折不足以认定统计显著。解决代码里加显著的过滤条件def find_significant_crossing(uf, ub, years, st): crossing_years [] for i in range(len(years) - 1): if (uf[i] - ub[i]) * (uf[i1] - ub[i1]) 0: # 检查交叉点处的UF是否在置信区间内 mid_uf (uf[i] uf[i1]) / 2 if abs(mid_uf) st: crossing_years.append(years[i 1]) return crossing_years只有在abs(mid_uf) st条件下的交叉点才输出。如果交叉处 UF 已经跑出置信区间说明趋势早已确立交点不是突变起点而是趋势持续点。坑 4序列带强自相关MK 检验的 z 值系统偏大现象把季节均值的原始序列交给 MK 跑结果几乎所有季节都检测出突变突变年份分布在全时间段看起来很假。查了发现原始序列的 lag-1 自相关系数高达 0.7。原因MK 检验假设数据独立但气候变量存在明显的年际持续性。正自相关会让 S 统计量的方差被低估导致 z 值偏大虚报突变。解决做预白化处理常见做法是先用 AR(1) 模型滤出自相关部分def prewhiten(seq): n len(seq) r1 np.corrcoef(seq[:-1], seq[1:])[0, 1] if abs(r1) 0.1: return seq # 生成预白化序列y_i x_i - r1 * x_{i-1} y np.zeros(n) y[0] seq[0] for i in range(1, n): y[i] seq[i] - r1 * seq[i-1] return y预白化会损失第一个数据点但我这里保留首值当作初始值。处理后再跑 MK突变点数量通常会回到合理范围。注意预白化后得到的突变年份是“去持续性后的突变”和原始序列的突变年份可能有 1~2 年偏差报告里要写清楚处理流程。坑 5序列长度不一致导致绘图错位现象某季节等序列因缺测被 dropna 后years数组和uf数组长度不匹配matplotlib直接报错“x and y must have same first dimension”。原因我在run_mk_on_seasons里seq是 dropna 后的而years seasonal_df.index.values[:len(seq)]用的是切片如果 NaN 在序列中间而不是末尾切片会错位。解决构造季节序列时先剔除无效年份保持年份和数值一一对应sub seasonal_df[season].dropna() seq sub.values years sub.index.values.astype(int)这类问题属于数据处理顺序没理清写脚本时要把“对齐”时刻放在心上否则画图画到一半才暴露最费时间。5. 进阶用滑动窗口与显著性检验验证突变点并批量跑多个站点5.1 三条验证手段交点在置信区间内、正反统计量交叉、滑动窗口稳定性单跑一遍 UF/UB 不能急着下结论我一般会加一个滑动窗口验证。思路是把季节序列按一个固定窗长比如 10 年切割成多个子段分别做 MK 突变检验看突变年份是否集中在某一段时间。如果 10 个窗口里 8 个都说突变在 2005 年前后这个突变才可信如果突变年份在 1995 到 2010 之间乱跳基本可以判定是噪声。def sliding_window_mk(seq_years, seq, window10): candidates [] n len(seq) for start in range(0, n - window 1, 3): # 每次滑动3年 sub_years seq_years[start:startwindow] sub seq[start:startwindow] uf, ub mk_uf_ub(sub) # 找交叉点 for i in range(len(sub_years) - 1): if (uf[i] - ub[i]) * (uf[i1] - ub[i1]) 0: candidates.append((sub_years[i1], uf[i1])) return candidates参数说明window10是滑动窗口长度小于 10 年的话 MK 的正态近似太勉强step3是滑动步长步长越小越平滑但计算量增大。返回的candidates是 (年份, UF 值) 列表统计年份频率分布就能看出突变年份是否稳定。另一个验证手段是看 UF 曲线是否在突变点之后持续穿越置信区间。真正突变的话UF 从某个时刻起会持续上升或下降并保持在置信区间外如果 UF 只在突变点附近探一下头又缩回来这种突变不值得深究。我常用的做法是算突变点之后 10 年的 UF 平均值绝对值大于 1.96 才算稳定突变。5.2 批量处理多站点与结果导出实际项目中很少只用一份数据可能一个流域几十个站点或者一套再分析资料里几十个格点。批量跑的时候我会把整条流程包成一个函数def batch_mk(file_list, save_dir): summary [] for f in file_list: df pd.read_csv(f, parse_dates[date], index_coldate) seasonal monthly_to_seasonal(df, value) results, st run_mk_on_seasons(seasonal, alpha0.05) for season, r in results.items(): for yr in find_crossings(r, st): summary.append([f, season, yr, significant if abs(r[uf]) st else candidate]) pd.DataFrame(summary, columns[file, season, cross_year, flag]).to_csv( f{save_dir}/mk_summary.csv, indexFalse) return summary批量导出的 CSV 里每一行是一个突变候选点后续可以在地图上打点或者和 ENSO、降水异常做关联分析。我习惯把“是否显著”直接写进表里避免事后翻图。跑完批量之后还有个习惯动作把突变前后两个时段的季节均值做差值算出突变幅度。比如春季序列 2005 年突变突变前 1980–2005 均值是 12.3°C突变后 2006–2023 均值是 13.8°C突变幅度 1.5°C这个数字比突变年份本身更能说明实际影响。从那以后我每次跑季节尺度 MK 都会强制走一遍“聚合—预白化—UF/UB—滑动窗口验证—输出幅度”这个流程少了哪一步都觉得不踏实。做突变检测最怕的就是拿着一条曲线就宣告“发现了气候突变”数据预处理和交叉验证省不了。希望帮到你。本文还有配套的精品资源点击获取