简介面向气候与环境数据分析人员这份资源提供了一个可直接运行的Python脚本及其配套的SPEI示例数据用于对季节尺度的干旱指数序列执行Mann-Kendall突变检测。脚本覆盖Excel数据读取、缺失值检查、STL季节分解、趋势成分M-K检验以及突变点可视化输出的完整流程可帮助定位干旱状况发生显著变化的年份适合需要快速上手非平稳时间序列分析的气象、水文从业者或相关专业学生。资源包共2个文件以.py脚本为主、.xlsx数据为辅压缩包大小仅11KB轻量易用已有196人学习下载。借助脚本中的STL季节分解思路与M-K检验实现读者既能复现SPEI指数突变点分析也可将方法迁移至其他季节数据开展扩展研究甚至用于科研图表复刻与教学演示是理解统计检验原理与Python实践相结合的实用参考。1. 一个反直觉的结论年尺度的 M-K 检测会漏掉季节突变年初我在处理某气象站的降水数据时把整条年序列直接丢进 M-K 突变检验结果一个突变点都没有测出来。换成“基于 Python 实现季节尺度 M-K 突变检测.py”按春、夏、秋、冬分别跑冬季的 UF/UB 曲线在某个年份出现了清晰交叉交叉点还正好落在置信区间内。后来我反复确认才发现年尺度的平均效应会把季节突变信号磨平年均值看起来稳定但单个季节内部可能早就发生了均值跳跃。这个脚本就是解决这个问题的——它把 Mann-Kendall 突变检测从“整条时间序列”细化到“每个季节独立检测”自动完成季节切分、UF/UB 曲线计算、突变年份提取和显著性判定。适合手里有日尺度或月尺度长序列、想找数据变点的数据分析同学尤其是水文、气象、生态这类长时间观测场景不依赖原始数据服从正态分布拿到序列就能直接用。2. 项目定位与技术选型为什么 M-K 要按季节拆开检测2.1 非参数检验原理为什么长序列变点检测首选 M-KM-KMann-Kendall检验之所以在气象水文领域用了这么多年核心原因是它几乎不做分布假设。常规的 t 检验要求数据近似正态、方差齐性但降水、径流、气温这类观测序列经常带明显的偏态和离群值一条强降雨记录就能把均值拉偏。M-K 只比较两个点之间的大小关系大于记 1小于记 -1相等记 0然后累计这些符号。这样极端值最多影响几个符号比较不会像均值那样被整体带偏。突变检测用的不是单个统计量而是两条曲线UF 和 UB。UF 是从序列开头向结尾逐点累积计算出来的标准化统计量反映“从前往后看”的趋势变化UB 是把序列倒过来再算一遍标准化统计量取负号后按原序排列反映“从后往前看”的趋势变化。如果序列确实在某一年发生了均值跳跃UF 和 UB 通常会在那一年附近形成交叉点。而这个交叉点是否可信要看它是否落在给定显著性水平的临界线以内常用的双边 95% 临界值是正负 1.96。把这个逻辑从年度序列搬到季节序列技术选型上的好处很明显季节序列本质上还是“一个季节一个样本点”的年序列M-K 的方差近似公式、临界值查表逻辑完全不用改只需要在数据分组阶段把“年”替换成“季节年”。这也是我后来在二次开发时觉得这个脚本设计舒服的地方——检验核心和数据预处理是解耦的换数据尺度不动算法。2.2 季节尺度的拆分方式自然季节、水文季节与跨年冬季季节尺度拆分听起来简单实际上有细节。最常见的做法是自然季节3、4、5 月为春季6、7、8 月为夏季9、10、11 月为秋季12、1、2 月为冬季。春、夏、秋三季都在自然年内完整闭环冬季比较特殊——它跨了一个自然年真正的“1971 年冬季”应该是 1970 年 12 月加上 1971 年 1、2 月。如果不做跨年处理直接把每年 1、2、12 月拼一起那么 1、2 月来自 1971 年12 月来自 1971 年但这个冬季的完整时间跨度应该是 1970 年 12 月到 1971 年 2 月。拼错之后冬季序列的第一个样本点会缺少前一年的 12 月最后一个样本点会缺少次年的 1、2 月看上去只是首尾各少一个点但 M-K 对样本量很敏感季节序列通常只有二三十年长度少一个点就会改变整个 UF/UB 的走向。水文上还有一个常见做法是水文年拆分比如从每年 10 月到次年 9 月作为一个水文年再把水文年分成四个三个月的小段。这个脚本里采用自然季节加跨年冬季处理我后面看代码时发现它把 12 月单独拎出来往后错一年再归组逻辑是对的。如果你手里有明确的枯水期、丰水期需求可以在二次开发时把季节映射表改掉这个我会在最后一章展开。2.3 项目结构与数据约定单脚本如何划分职责这个资源的核心是一个 Python 脚本配合示例数据和输出目录整体结构大概是season_mk/ ├── mann_kendall_seasonal.py # 主脚本含检测、绘图、结果导出 ├── requirements.txt # numpy / pandas / matplotlib / scipy ├── data/ │ └── example_daily_prcp.csv # 示例日降水数据date prcp 两列 └── output/ ├── spring_ufub.png ├── summer_ufub.png ├── autumn_ufub.png ├── winter_ufub.png └── result_summary.csv数据约定很直接CSV 文件至少包含一列日期和一列数值。日期列会被pd.to_datetime解析所以你传1971-01-01或者1971/1/1都能识别数值列要求是浮点数缺测值要么提前处理要么脚本里做插值。输入如果是日尺度脚本会先算月均值再按季节聚合如果已经是月度数据可以直接跳过中间步骤。输出物方面每个季节一张 UF/UB 曲线图外加一个汇总 CSV。汇总 CSV 里每行对应一个季节的检测结果包含突变年份、UF/UB 交点位置、显著性水平和突变前后均值差。这个 CSV 是后续批量分析的关键接口后面讲运行参数和二次改造都会用到它。3. 核心模块拆解UF/UB 曲线从计算到突变判定的完整链路3.1 mk_uf 函数顺序统计量与标准化计算整个脚本里最核心的函数是正序 UF 统计量的计算。它接收一个季节的逐年均值序列返回一条长度相同的 UF 曲线import numpy as np def mk_uf(x): 计算正序序列的UF统计量。 x: 某个季节的逐年均值序列例如 [458.2, 462.1, ..., 612.7] 返回: 与x等长的UF数组UF[0] 固定为0 n len(x) uf np.zeros(n) s 0.0 for i in range(1, n): # 第i个样本与之前所有样本做符号比较并累加 s np.sum(np.sign(x[i] - x[:i])) # M-K方差近似公式此时参与比较的样本数为 i var_s i * (i 1) * (2 * i 5) / 18.0 uf[i] s / np.sqrt(var_s) if var_s 0 else 0.0 return uf这段代码的关键在np.sign(x[i] - x[:i])它一次性完成当前样本与全部历史样本的符号比较结果只有 -1、0、1 三种然后累加到s上。随着i增大s会反映越来越长时间尺度上的趋势特征。var_s用的是无平局情形的方差近似公式对降水这类连续观测数据基本够用如果数值列是离散的整数等级可能需要额外加上平局修正项否则方差会被低估。这里有一个容易看漏的细节uf[0]被初始化为 0i从 1 开始循环。也就是说第一个样本点本身不参与任何符号比较它只是曲线的起点。后面画图时UF 曲线是从第 2 个点开始才有实际意义的。3.2 逆序序列 UB 曲线镜像计算与交叉点提取UB 曲线的计算比 UF 多两个步骤逆序、取负、再反序。脚本里的实现长这样def mk_detect(x, alpha0.05): n len(x) uf mk_uf(x) # 逆序计算UF取负号后反序排列得到UB ub -mk_uf(x[::-1])[::-1] # 双边显著性水平对应的临界值 from scipy import stats z_crit stats.norm.ppf(1 - alpha / 2) # 找UF与UB的交叉点符号乘积由正变负或由负变正 cross_idx [ i for i in range(1, n) if (uf[i - 1] - ub[i - 1]) * (uf[i] - ub[i]) 0 ] return uf, ub, cross_idx, z_crit为什么要对逆序序列取负号原因在于方向一致性。对x[::-1]做正序 UF 计算得到的是“从末端往回首”的趋势累积量。如果原序列在某处突然上升那么逆序序列在对应位置看就是下降两者的符号天然相反。为了让 UB 和 UF 画在同一坐标系里比较必须把逆序结果取负再反序回来这样两条曲线统一了方向交叉点才对应突变候选年份。交叉点的提取用的是相邻两点符号乘积小于 0 的判别方式也就是两条曲线在某个区间内发生了穿越。这里我实际跑数据时发现如果序列存在连续多年的趋势波动UF 和 UB 可能交叉好几次交叉点未必全是真实突变。这个脚本的判定方式是看交叉点是否落在临界线以内同时看突变前后均值差是否显著我后面在避坑部分会再展开。3.3 季节聚合跨年冬季与逐季节均值序列的生成季节聚合这步直接决定后续 M-K 算的是什么东西。脚本里用映射表加跨年偏移来处理冬季import pandas as pd def to_seasonal(df, value_colvalue): df df.copy() df[date] pd.to_datetime(df[date]) df[month] df[date].dt.month df[year] df[date].dt.year # 自然季节映射12月归入冬季 season_map { 12: winter, 1: winter, 2: winter, 3: spring, 4: spring, 5: spring, 6: summer, 7: summer, 8: summer, 9: autumn, 10: autumn, 11: autumn } df[season] df[month].map(season_map) # 跨年处理12月并入下一个自然年的冬季 df[season_year] df[year] df.loc[df[month] 12, season_year] df[year] 1 # 按季节季节年分组求均值得到每季一个样本点 seasonal df.groupby([season, season_year])[value_col].mean() return seasonal.unstack(level0)这段代码对初学者有两个容易绕晕的点。第一个是season_year的偏移1971 年 12 月的数据被标记为season_year 1972于是它跟 1972 年 1、2 月分到同一个冬季组里这个冬季的标签是 1972。第二个是unstack(level0)的效果分组后索引是季节和季节年两层unstack把季节变成列每一列就是一条独立的季节均值序列正好可以喂给 M-K 函数。输出之后每个季节序列的长度等于有效年数而不是总月数这点在参数估计里很重要。M-K 的方差公式对样本量是很敏感的冬季序列如果因为跨年处理导致首尾缺数据长度会比其它季节少一两个点显著性判断可能直接翻车。这也是我建议拿到脚本先把seasonal.shape打印出来看一眼的原因。3.4 突变点判定规则与输出字段说明脚本最终不是简单输出所有交叉点而是按一套规则过滤。我根据代码逻辑整理成了这样一个判定链判定步骤具体条件说明候选点UF 与 UB 交叉符号乘积出现负值显著性过滤交叉点位于正负临界线之间超出临界线的交叉点多为强趋势变化不代表均值突变均值差验证交叉点前后各取 5 年均值计算差值差值过小视为振荡连续交叉约束交叉发生后的 5 年内趋势方向保持一致避免短期波动被误判为突变落在临界线外侧的交叉点在气象应用里通常被解释为序列整体趋势剧烈改变而不是局部突变。脚本默认采用这个解释框架所以它综合了交叉位置和均值差两个信息。输出的 CSV 里包含的字段大约有season、trend_z、突变年份、突变点 UF/UB 值、突变前均值、突变后均值、p_value 等。我自己在项目里会把 trend_z 单独拿出来做趋势显著性排序用它筛选掉那些交叉点勉强存在但整体趋势不强的季节。4. 运行指南数据准备、参数配置与结果解读4.1 环境准备依赖安装与快速验证脚本依赖的第三方库是 numpy、pandas、matplotlib、scipy。前三个做数据处理和画图scipy 用来查标准正态分布临界值。建议直接用虚拟环境装python -m venv venv source venv/bin/activate # Windows 下改为 venv\Scripts\activate pip install numpy pandas matplotlib scipy装完之后可以先跑一个最简单的自检构造一条长度 40、中间年份均值跳变 1.5 个标准差的合成序列喂给检测函数看脚本能不能在预设年份附近找到突变点。这一步只需要几秒钟但它能确认你本地环境和脚本逻辑是匹配的避免后面拿着真实数据跑出结果却不知道算得对不对。4.2 一条命令跑通全流程参数说明数据准备成 CSV 后运行命令大概是这样的python mann_kendall_seasonal.py \ --input data/example_daily_prcp.csv \ --date-col date \ --value-col prcp \ --alpha 0.05 \ --min-years 8 \ --outdir output核心参数含义如下参数默认值作用--input无输入 CSV 路径必填--date-coldate日期列名脚本用pd.to_datetime解析--value-colvalue数值列名传给季节聚合--alpha0.05显著性水平双边检验临界值随之变化--min-years8少于该年数的季节直接跳过避免样本量不足--outdiroutput输出目录不存在会自动创建--min-years这个参数容易被忽略但它很实用。M-K 的方差近似公式在小样本下会低估不确定性季节序列如果只有六七年数据算出来的 UF/UB 曲线波动会非常剧烈交叉点一抓一大把但基本没有统计意义。设置最小年数可以自动过滤掉这种季节我一般习惯设为 10。4.3 结果文件怎么读突变年份、显著性、均值差跑完之后output目录下会生成四个季节的 PNG 图和一份result_summary.csv。CSV 打开大概是这样的格式season 突变年份 UF UB Z统计量 p_value 突变前均值 突变后均值 是否显著 spring 1998 2.81 -2.42 2.81 0.0049 458.2 612.7 是 summer 2003 1.32 -1.05 1.32 0.1867 522.4 506.8 否 autumn 无 0.87 0.91 0.87 0.3843 477.2 480.1 否 winter 1994 2.47 -2.31 2.47 0.0135 120.4 89.6 是读这个表时我自己的习惯是先看是否显著列把显著的行筛出来再看突变年份结合原始数据确认这一年附近是否真的发生了极端事件或政策变化。比如上面这个示例里春季突变后均值从 458.2 升到 612.7差值超过 150 毫米这通常说明不是随机振荡值得进一步追踪原因。夏季的突变年份虽然存在但 p 值不显著均值差也只有 15 毫米这种我会直接放下。4.4 参数调整alpha、季节过滤与输入频率实际使用中你大概率需要调三个地方。第一个是--alpha如果只想保留最强的突变信号可以收紧到 0.01临界值会从 1.96 变成 2.58很多模棱两可的交叉点会被过滤掉如果反过来做探索性分析用 0.1 让更多候选点浮出来再逐个排查。第二个是季节过滤默认四个季节全算但如果你只关心冬季降水可以在代码里把to_seasonal返回结果的列名筛一下或者加一个--seasons winter参数省去画其它季节图的时间。第三个是输入频率。脚本默认从日数据开始聚合成月再聚合到季但如果你的数据源本身已经是月尺度建议直接复用to_seasonal函数不要先拆成日再重采样一遍。月尺度数据只需要在读取后直接按month分组少一层日到月的聚合也能避免月初月末的边界时间戳被解析出问题。5. 避坑记录季节尺度 M-K 检测的五个典型翻车点5.1 冬季序列首尾缺样本突变年份偏差一年现象跑完结果发现冬季序列的突变年份总是和春季差一年或者冬季有效样本量比其它季节少。原因数据从 1 月开始第一年的冬季只有 1、2 两个月缺上一年 12 月最后一年的冬季又缺次年 1、2 月。解决办法在to_seasonal里对season_year完成偏移后把首尾不完整的季节样本直接丢弃用dropna()清理。从那以后我每拿到新数据都会先打印每个季节的年份范围确认首尾年份一致才开始检测。5.2 UB 曲线方向画反交叉点全部错位现象UF 和 UB 两条曲线几乎平行看不出交点或者交叉点分布在完全不符合直觉的位置。原因计算 UB 时没有对逆序结果取负号画出来的曲线实际上反映的是“逆序趋势”而不是“原序逆推趋势”。解决办法严格按ub -mk_uf(x[::-1])[::-1]计算并且用已知有突变点的合成序列做验证。这个坑最容易在照着论文公式手搓代码时出现因为论文里经常只写“对逆序列重复上述过程”省略了符号处理和反序排列。5.3 缺测值没有预处理方差公式分母异常现象代码运行不报错但输出的 UF 曲线在某一段突然跳变成很大的值或者出现 NaN。原因pandas 在 groupby 聚合时默认跳过 NaN但如果某个季节内三个月两个有值一个缺测均值会被低估季节序列出现明显凹坑符号比较逻辑被打乱。解决办法在进入mk_uf之前对季节序列做线性插值或者用seasonal.interpolate(methodlinear)补齐。那种一个季节全部缺测的情况直接丢弃该年份不要硬补。5.4 多季节同时检测显著性泛滥现象四个季节跑完每季都能看到交叉点甚至两三处交叉真实原因却说不清楚。原因多重检验问题。对同一个站点做 4 次独立检验即使数据完全是随机噪声出现至少一个“显著交叉点”的概率也会明显高于单次检验的 0.05。解决办法一种方式是收紧 alpha 到 0.01另一种方式是要求交叉点前后连续 5 年的趋势方向一致才确认为突变。实际项目里我更喜欢后者因为它利用了时间连续性比单纯降 alpha 更接近气象意义上的突变。5.5 季节样本量不足正态近似失真现象某个季节只有 6~7 年数据UF/UB 曲线的波动幅度却非常大置信区间形同虚设。原因M-K 的方差近似公式在大样本下才稳定水文站常见的情况是某条支流的数据只完整记录了 5 个汛期硬跑检验没有统计意义。解决办法在脚本里强制--min-years 10把样本量不足的季节从结果表中排除。样本量不足时与其挑一个侥幸显著的交叉点不如承认数据还不够支撑突变判定。6. 二次开发把单文件脚本改造成批量检测工具6.1 抽出检测核心快速接入对照方法如果你手里的数据不止一个站点第一个重构动作应该是把mk_detect的输入从“单条序列”改成“可迭代的序列集合”。这个改造只要在函数外层加一个循环def batch_detect(all_series: dict, alpha0.05): rows [] for (site, season), series in all_series.items(): if len(series) 10: continue uf, ub, cross_idx, z_crit mk_detect(series.values, alpha) for idx in cross_idx: if -z_crit uf[idx] z_crit: rows.append({ site: site, season: season, year: series.index[idx], z: uf[idx] }) return pd.DataFrame(rows)同时可以再加一个 Pettitt 检验做交叉对照。Pettitt 是另一种非参数突变检验它的输出只有一个突变点和 M-K 的 UF/UB 交叉点互相印证。两个方法都指向同一年份时这个突变点的可信度会高很多。6.2 多站点、多季节结果汇总成一张宽表批量检测跑完后把结果pivot成站点行、季节列的宽表就能快速对比不同季节的突变年份分布。常见做法是热力图横轴是季节纵轴是站点颜色表示突变年份早晚。这个可视化能直接暴露区域性的突变同步性比逐张看图高效得多。我在某次多站点分析里发现周边五个站点的冬季突变年份集中在前后三年内只有最上游的站点要晚五年验证了那是一次区域性气候调整而非局部数据噪声。6.3 用合成序列验证每一次改动改动代码之后先用构造数据做回归测试是我的习惯。合成序列的核心逻辑是先生成随机噪声再从某个指定年份开始给均值加一个偏移量最后让脚本检测这个偏移发生在哪一年。如果脚本返回的突变年份和预设年份偏差超过两年说明改动引入了问题。把这条验证链路写成pytest测试函数之后后续无论怎么改季节映射、换临界值算法都能在几分钟内确认逻辑没有退化。这类时序突变检测工具最大的风险是“结果看起来合理但过程有隐蔽错误”。从那以后我每次跑季节尺度 M-K 之前都会强制走一遍合成数据自检确认当前环境下的 UF/UB 计算、季节聚合、突变点提取三个环节都正常再去分析真实序列。希望这个脚本能让你少踩一些我踩过的坑拿到数据后先把季节序列结构摸清楚再放心读突变结果。本文还有配套的精品资源点击获取