简介基于Python实现的季节尺度M-K突变检测脚本面向气候、水文和生态环境研究中常见的SPEI、径流、气温等具有明显季节周期的时间序列数据旨在帮助科研人员客观识别序列中的趋势突变点为干旱演变、气候转型等分析提供统计依据。整个压缩包以zip格式提交共2个文件包含1个Python脚本和1个xlsx数据文件脚本覆盖数据读取、缺失值检查、STL季节分解、Mann-Kendall检验及突变点判读等完整流程压缩包整体仅11KB便于下载与复用。目前已有195人浏览/学习适合气候、环境及数据分析方向的学生与科研人员参考。脚本除了直接输出突变点位置和显著性指标外还会生成包含原始序列、季节成分、趋势成分与突变点位置的对比图便于直观理解分解与检验结果由于Mann-Kendall检验基于秩次、不依赖数据分布算法稳健性好替换数据后即可迁移用于其他季节尺度序列检测。 做水文气象分析的人应该都绕不开突变检测。无论是分析降水的年代际变化还是评估某个站点气温序列的稳定性M-KMann-Kendall检验几乎是最常用、也最容易被审稿人认可的非参数方法。我最早接触它是在处理长时间序列的降水数据当时直接用现成的R包跑完全年序列结果出来一个不温不火的突变点导师看了一眼就问“你按季节拆开看过吗”这一问问住了我。全年序列的M-K检验会把季节差异“平均”掉很多时候明明春季降水的突变点发生在1998年夏季在2005年合在一起反而什么都看不出来。后来我把这套逻辑整理成了基于Python的季节尺度M-K突变检测脚本今天把完整思路、代码和踩过的坑都写出来希望能帮你少走点弯路。1. 为什么要做季节尺度的M-K突变检测1.1 全年序列检测的局限先举个直观的例子。假设你有一份1960年到2020年的月降水数据站点在华北地区。如果直接对逐年降水总量做M-K检验会发现UF和UB曲线在置信区间内交缠在一起很难找到一个清晰的突变年份。这是因为华北的降水年际变率大春季和秋季的降水趋势可能完全相反全年总量序列把所有季节的信号揉在一起趋势相互抵消突变点自然就模糊了。另外一个问题是气候研究里“季节”本身就有明确的物理意义。夏季降水的突变往往和季风环流的调整有关冬季气温的突变又可能受到西风带位置变化的影响。把不同物理机制的信号放在同一个序列里检验统计上说得通但气候学解释上站不住脚。所以对季节尺度分别做突变检测不是“精细化”的加分项而是必要步骤。1.2 季节尺度检测的分析策略季节尺度的核心思想很简单把原始月序列按季节分组比如春季取3月、4月、5月的数据然后对每个季节的序列单独做M-K检验。但这里有一个关键细节——“季节序列”的构建方式有两种。一种是用季节累积量比如春季总降水量、夏季平均气温另一种是用季节代表值比如夏季平均气温用6、7、8三个月的均值。我建议做降水分析时用累积量做气温分析时用均值这符合气象业务里的习惯。实际写代码时我用的是pandas的groupby按月份提取子序列再对每个季节子序列分别计算M-K统计量序列最后把四个季节的结果放在一张图里对比。这样做的好处很明显每个季节的突变年份清晰独立而且可以对比不同季节突变时间的前后关系。比如我处理过的某个站点数据里春季降水突变点出现在1998年前后夏季降水突变点出现在2005年前后这种时间差本身就值得写进文章讨论里。2. M-K检验的核心原理与关键代码逻辑2.1 统计量UF和UB的计算过程M-K检验的突变检测部分核心是计算两个统计量序列UF和UB。UF是正序列的标准化统计量UB是反序列计算后倒序处理得到的。两条曲线在置信区间内的交点对应的时刻就是突变开始的时间。计算UF的过程分几步。对于长度为n的序列x对每个时刻i统计第i个样本之前所有样本中大于第i个样本的个数记为s_i。然后计算累计值S_k S_{k-1} s_k。在序列独立且随机的原假设下S_k的均值为0方差由公式Var(S_k) n(n-1)(2n5)/18给出。最后标准化UF_k (S_k - mean) / sqrt(Var)。UB的计算是把原序列倒序后用相同的过程算一遍然后对结果取负并倒序排列。这里有一个容易踩坑的地方很多资料里对UB的表述很含糊如果直接把反序列的统计量拿来用而不做取负和倒序画出来的UB曲线方向是反的交点位置全错。我最初就栽在这里结果把突变年份判错了三年。2.2 为什么直接用现成库还不够Python生态里其实有现成的M-K检验库比如pymannkendall功能很全能算趋势检验、季节检验甚至能处理结数据。但实际用下来我发现直接用库有两个不顺手的地方。第一pymannkendall主要输出的是趋势检验结果和突变点位置但画不了UF和UB曲线图。可实际业务里审稿人也好、导师也罢最想看的恰恰是这两条曲线的交叉图。第二季节尺度的处理逻辑并不复杂但用库时要对四个季节分别调一遍接口再手动整理结果代码量反而比直接实现核心函数要多。所以我最后选择自己实现核心的M-K统计量计算函数配合matplotlib画图。这样做还有一个好处你可以完全控制数据处理流程比如在进入统计量计算之前可以先对序列做预处理去除异常值或者做三阶平滑这在用现成库的时候很难插进去。3. Python完整实现从数据准备到突变点判定3.1 模拟数据生成方便复现先说数据准备。为了让你能直接跑通代码我构造了一份模拟数据从1960年到2020年每月一条降水数据春季降水在1998年之后显著增加夏季降水在2005年之后有趋势变化其他季节趋势不明显。import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy import stats # 生成模拟数据 np.random.seed(42) years np.arange(1960, 2021) months np.arange(1, 13) dates pd.date_range(1960-01-01, 2020-12-31, freqMS) n len(dates) # 基础降水季节正弦波动 加上 年际随机噪声 seasonal 30 20 * np.sin(2 * np.pi * (dates.month - 2) / 12) noise np.random.normal(0, 8, n) precip seasonal noise # 春季3-5月在1998年后增加10mm spring_mask (dates.month 3) (dates.month 5) (dates.year 1998) # 夏季6-8月在2005年后线性增加 summer_mask (dates.month 6) (dates.month 8) (dates.year 2005) precip[spring_mask] 10 trend_increase np.linspace(0, 12, (dates.year 2005).sum()) precip[summer_mask] trend_increase df pd.DataFrame({time: dates, year: dates.year, month: dates.month, precip: precip}) print(df.head())注意我用numpy的随机种子固定了随机数这样你本地跑出来的图和我的完全一致方便你检查代码逻辑哪里出了问题。3.2 季节分组与M-K统计量计算接下来是核心部分季节分组和M-K统计量计算。先把月份映射成季节我这里采用气象季节的划分3-5月春季、6-8月夏季、9-11月秋季、12-2月冬季跨年处理。def get_season(month): if month in [3, 4, 5]: return 春季 elif month in [6, 7, 8]: return 夏季 elif month in [9, 10, 11]: return 秋季 else: return 冬季 df[season] df[month].apply(get_season)然后实现M-K统计量计算函数。这个函数对输入序列计算UF和UB序列同时给出趋势检验的p值和tau值。def mk_statistic_series(data): 计算M-K检验的UF和UB序列 返回: uf, ub, tau, p_value n len(data) # 计算正序列的UF s np.zeros(n) for i in range(1, n): s[i] s[i-1] np.sum(data[:i] data[i]) uf np.zeros(n) for k in range(1, n): mean_s k * (k - 1) / 4 var_s k * (k - 1) * (2 * k 5) / 72 uf[k] (s[k] - mean_s) / np.sqrt(var_s) # 计算反序列的UB data_rev data[::-1] s_rev np.zeros(n) for i in range(1, n): s_rev[i] s_rev[i-1] np.sum(data_rev[:i] data_rev[i]) uf_rev np.zeros(n) for k in range(1, n): mean_s k * (k - 1) / 4 var_s k * (k - 1) * (2 * k 5) / 72 uf_rev[k] (s_rev[k] - mean_s) / np.sqrt(var_s) # UB 是反序列统计量的负值再倒序 ub -uf_rev[::-1] # 计算tau和p值用scipy做趋势检验 tau, p_value stats.kendalltau(data, np.arange(n)) return uf, ub, tau, p_value这里有个细节值得解释正态化方差用的分母和我在2.1节写的公式略有不同。原因是在代码实现里s_k的方差需要根据当前的k值动态计算所以用k(k-1)(2k5)/72而不是用整个序列长度n的公式。这是一个常见的实现细节很多初学者在这含糊过去但直接决定曲线的置信区间是否合理。3.3 累计距平辅助判断与可视化有了UF和UB序列还不够。我强烈建议同时计算累计距平曲线它的作用是辅助验证突变方向。累计距平的计算很直观先算序列的均值再用逐点累积值和均值的差。def cumulative_anomaly(data): mean_val np.mean(data) anomaly data - mean_val cumsum np.cumsum(anomaly) return cumsum累计距平曲线如果从下降转为上升说明序列从偏低阶段转入偏高阶段这和UF曲线穿越置信区间上界的方向应该是对应的。我见过很多文章只用UF/UB交叉点判断突变结果交叉点出现在置信区间外或者方向矛盾这种结果拿到审稿人那里站不住脚。累计距平曲线补充的是“趋势方向”的证据。最后把四个季节的UF/UB曲线、累计距平曲线放在一起画出来。画图的时候我习惯用2×2的子图每个子图一个季节上下两行上面放UF/UB曲线下面放累计距平曲线。seasons [春季, 夏季, 秋季, 冬季] fig, axes plt.subplots(4, 2, figsize(12, 16)) for idx, season in enumerate(seasons): season_data df[df[season] season].groupby(year)[precip].sum() years_arr season_data.index.values values season_data.values uf, ub, tau, p mk_statistic_series(values) cum cumulative_anomaly(values) ax1 axes[idx, 0] ax1.plot(years_arr, uf, labelUF, colorsteelblue) ax1.plot(years_arr, ub, labelUB, colorcoral) ax1.axhline(1.96, linestyle--, colorgray, linewidth0.8) ax1.axhline(-1.96, linestyle--, colorgray, linewidth0.8) ax1.set_title(f{season} M-K突变检验) ax1.legend() ax1.set_xlabel(年份) ax1.set_ylabel(统计量) ax2 axes[idx, 1] ax2.plot(years_arr, cum, colordarkgreen) ax2.set_title(f{season} 累计距平) ax2.set_xlabel(年份) ax2.set_ylabel(累计距平) plt.tight_layout() plt.show()运行这段代码你会看到春季的UF和UB曲线在1998年前后交叉夏季在2005年前后交叉秋季和冬季没有明显的交叉点。累计距平曲线也在对应的年份附近出现转折。这说明模拟数据里设置的突变点被正确识别出来了整个检测流程是有效的。4. 实测中的常见问题与避坑经验4.1 季节边界的定义要统一我最初犯过一个低级错误在数据分组时把12月、1月、2月当作冬季但没有跨年处理。也就是说冬季序列把1960年12月、1961年1月、1961年2月拼在一起与1961年12月、1962年1月、1962年2月混在一起边界不干净。正确的做法是如果采用12-2月作为冬季应当把年份标记重新对齐。具体来说把12月的数据归到下一年的年份标签里。我在代码里是这样处理的df[season_year] df[year] df.loc[df[month] 12, season_year] df.loc[df[month] 12, year] 1这个细节看起来小但影响很大。如果不处理冬季序列的第一个点和最后一个点可能不是完整的冬季年份M-K统计量在序列两端的表现会异常导致在边缘年份出现虚假的交叉点。4.2 样本量与结数据的处理季节尺度的M-K检验样本量等于年数不是总月数。比如60年的数据每个季节序列只有60个点。M-K检验对样本量没有严格要求样本量小的情况下检验功效会下降也就是说突变点可能检测不出来。如果样本量少于20年我建议直接用计算UF/UB曲线的方式做诊断但结论要谨慎不要只依赖p值说话。另外降水数据有个特殊性很多月份的降水量可能为0导致序列里出现大量结ties。M-K检验对结数据并不友好原假设要求序列连续计算方差时也需要考虑结的影响。pymannkendall库内置了修正方法但自己实现时如果发现UF曲线的走势特别“毛糙”先检查是不是零值太多。我习惯在进入M-K计算之前先统计一下序列中相同值的占比超过20%就要考虑使用修正方差公式或者改用季节M-K检验Seasonal Mann-Kendall那是另一个话题这里不多展开。4.3 交叉点不等于突变点这是新手最容易理解偏的地方。UF和UB曲线交叉不一定就是突变点。统计上的严格定义是交叉点必须落在±1.96的置信区间内才认为是显著的突变点。我实测过一个站点UF和UB在2008年交叉了但交叉的纵坐标是1.3小于1.96说明这个交叉在统计上不显著。如果只按交叉点判定突变年份就相当于把统计噪声当成了气候信号。碰到这种情况我会再看累计距平曲线有没有明显转折如果也没有就判定该序列无显著突变。另外有时UF曲线在置信区间外用了一段时间后UB曲线才跟上来交叉而交叉点已经出了置信区间。这种情况下突变发生的实际时间应该取UF曲线首次超出置信区间的时刻。这个判定规则我是从一篇文献里学到的后来实测确实比单纯取交叉点更准确。5. 结果解读与后续扩展5.1 判定突变的参考规则最后总结一下我实际使用时的判定规则你可以直接抄作业如果UF和UB曲线在置信区间内出现交叉点交叉点对应年份为突变年份如果交叉点在置信区间外录取UF首次超过临界线的年份作为突变年份辅助以累计距平曲线确认如果多条曲线多次交叉优先选择交叉后分离方向稳定、持续时间超过5年的节点如果春夏秋冬四个季节里只有一个季节检测出突变点写文章时仍然可以说“该区域XX季节降水在XX年发生突变”但不能推广为“该区域降水发生突变”。5.2 扩展多站点批量检测和显著性检验这份代码只处理了单站点数据但在实际研究里你可能要面对全省上百个站点的数据。扩展思路很简单写一个外层循环把每个站点的数据都传入mk_statistic_series函数然后把突变年份提取出来做出空间分布图。另外M-K检验结果里的tau值和p值也值得充分利用。tau值反映趋势强度p值反映显著性。我习惯在每个季节子图右上角加一行文字标注tau和p值让读者一眼看出趋势情况。ax1.text(0.02, 0.95, ftau{tau:.3f}, p{p:.3f}, transformax1.transAxes, vatop, fontsize10, bboxdict(facecolorwhite, alpha0.7))关于显著性检验还有一个细节如果做了多个季节的检验需要考虑多重比较问题。4个季节同时检验纯靠运气也可能出现一个季节显著。稳妥的做法是在文章里汇报原始p值同时注明“经过FDR校正后夏季的突变仍然显著”。这样一来审稿人挑不出毛病。我个人在实际操作中的体会是M-K突变检测的代码本身并不难写真正值钱的是对数据含义的理解以及对季节划分、边界处理、显著性判定这些细节的把控。你在跑通这套代码之后强烈建议用自己站点的数据做一遍完整的季节尺度检测再把结果和已知的物理机制对比——当春季降水的突变年份恰好和一次大尺度的环流调整对应上时那种感觉比调通代码本身还要有成就感。本文还有配套的精品资源点击获取