Python实现季节尺度M-K突变检测：UF/UB曲线定位突变年份

简介基于Python实现季节尺度M-K突变检测.py是一份面向气候与环境科学研究者、数据分析人员的实用脚本资源聚焦SPEI指数等季节性时间序列的突变点识别解决通用M-K检验难以直接处理周期波动数据的问题。资源共2个文件包含1个Python脚本与1个SPEI3.xlsx示例数据文件压缩包整体仅11KB结构精简便于快速启动项目。脚本覆盖数据读取、缺失值检查、STL季节分解、Mann-Kendall突变检测及结果可视化完整流程可直接替换数据运行相比通用M-K检验该脚本特意先分解出季节成分与趋势成分再对趋势项计算Z-score和P-value能更准确地捕捉干旱等气候指标的突变年份并输出突变点位图供深入研判。目前已有195人学习下载适合需要系统掌握季节性M-K突变检测实现路径的入门及进阶用户尤其对干旱监测、水资源管理方向的研究者极具参考价值既可作为教学示例也能作为方法模板开展二次开发。1. 季节尺度M-K突变检测一个被年序列掩盖的信号做水文或气象趋势分析时我经常遇到一类现象某站点近60年降水年序列的Mann-Kendall趋势检验“不显著”但把数据按季节拆开夏季降水量从某个年份开始明显减少冬季反而小幅增多。年序列把方向相反的季节信号互相抵消了突变点也就淹没在均值里。基于Python实现季节尺度M-K突变检测.py 这个脚本要解决的就是把“按季节拆开、逐季节做M-K突变检验”这件事固化下来输入月值数据按季节聚合为四组逐年序列再分别计算UF/UB统计量并定位突变年份。适合水文、气象、环境监测、生态遥感领域需要批量处理站点序列的分析人员也适合学生毕业论文里需要做突变检验但不想每次复制粘贴代码的读者。下面这套实现不依赖GIS类重型库用pandas、numpy、scipy和matplotlib就能跑通。2. M-K突变检测原理与季节尺度的必要性UF/UB曲线为何能定位突变2.1 Mann-Kendall突变检验的统计原理Mann-Kendall突变检验建立在秩序列之上不要求样本服从正态分布也不受少数异常值影响因此在水文气候领域被广泛当作“看突变年份”的佐证工具。它的核心是构造两个序列UF正向秩序列统计量和UB反向秩序列统计量。设时间序列 x1, x2, …, xn对每个时刻 i统计 xi 大于其之前所有时刻值的个数记为 ri。定义秩序列s_k Σ_{i1}^{k} ri当 k 从 2 到 n 时s_k 的期望和方差有明确的解析表达式E(s_k) k(k-1)/4Var(s_k) k(k-1)(2k5)/72于是 UF_k (s_k - E(s_k)) / sqrt(Var(s_k))k2,…,n。UF 是一条随年份推进变化的标准化统计量曲线超过 ±1.96即0.05显著性水平表示该时段存在显著上升或下降趋势。UB 则是把序列倒过来重算一遍再取负号并翻转时间轴相当于从序列末端向起点回溯得到的相同统计量。如果序列在某一时刻发生真实突变UF 和 UB 会在这个时刻前后出现交叉——交叉点对应的年份就是突变候选年份。注意交叉点本身只是一个统计现象还需要结合交叉点附近两条曲线是否都处于置信区间内以及交叉后曲线是否继续分离并越界综合判断突变是否真实显著。2.2 季节尺度与年尺度的差异为什么必须按季节拆开年尺度序列把一年内12个月的值累加或平均这相当于对季节扰动做了一次低通滤波。举个我实际处理过的例子某站夏季降水逐年减少、冬季降水逐年增加年总量基本持平。对该站的年降水序列做M-K检验UF曲线在置信区间内小幅波动根本找不到突变点但把夏季降水单独提出来做M-K发现UF在1997年前后迅速从正的平稳区跌穿-1.96UB从另一端交叉过来突变特征非常清晰。下表是同一组模拟数据在年尺度和季节尺度下的检测结果对比序列突变年份趋势方向是否显著年总量未检测到基本平稳否春季1-3月2003上升是夏季4-6月1997下降是秋季7-9月无无趋势否冬季10-12月2010上升是这就是“季节尺度”的价值它把一条年序列拆成四条平行序列每一条都保留了该季节的多年演变信息。做季节尺度的M-K突变检测不是对整条月值序列做一次检验而是先把每年同一季节的数据聚合为一个值形成“年份—季节值”序列再对每个季节的序列分别做M-K检验。换句话说春季序列的长度等于观测年数不是月数。2.3 环境依赖与核心库选型这个脚本只需要四个库numpy负责向量化秩次计算scipy.stats提供标准正态分位数pandas完成季节聚合和数据清洗matplotlib用于绘制UF/UB曲线。安装命令简洁pip install numpy pandas scipy matplotlib不推荐用statsmodels里的趋势检验代替因为M-K突变检验需要同时输出正反向统计量statsmodels只提供趋势显著性检验没有UF/UB曲线交叉这一套。另外要注意scipy版本对norm.ppf结果精度的影响很小但不同版本对数组索引和csr_矩阵等无关功能没有影响真正需要关注的是pandas版本在groupby之后时序索引的行为我会在对应章节标注。3. 数据准备与季节划分从月度观测到四组年际序列3.1 构造标准的月值时间序列表假设你手头数据是站点观测的月值最省心的格式是三列year、month、value。无论原始数据来自Excel还是TXT先统一转成这种长表。下面代码用pandas读取并做基础检查import pandas as pd # 原始csv至少包含 year, month, value 三列 df pd.read_csv(monthly_values.csv) # 检查缺失和类型 print(df.head()) print(df.isnull().sum()) # 确保数值列不是字符串 df[value] pd.to_numeric(df[value], errorscoerce) df df.dropna(subset[value, year, month]) # 构造时间索引便于后续处理 df[time] pd.to_datetime(df[year].astype(int).astype(str) - df[month].astype(int).astype(str) -01)说明pd.to_numeric(errorscoerce)会把无法转换的脏数据变成NaN随后用dropna剔除。这一步必做否则后面groupby聚合时字符串类型会导致求和直接拼接或报错。如果你的数据本身就带时间列也可以用pd.to_datetime统一解析再提取dt.year和dt.month效果一样。3.2 季节映射与跨年问题季节划分方式直接决定结果含义。我这里采用自然年四等分1-3月春季、4-6月夏季、7-9月秋季、10-12月冬季用一个简单函数映射def season_of_month(month): 自然年四等分季节返回字符串代号。 if month 3: return spring elif month 6: return summer elif month 9: return autumn else: return winter df[season] df[month].apply(season_of_month)为什么不用气候学标准的DJF12-2月、MAM3-5月、JJA6-8月、SON9-11月因为12月所属年份存在跨年归属问题2023年1月和2月应该与2022年12月组成同一个冬季但按自然年分组时2022年12月会被分到2022年冬季2023年1-2月分到2023年冬季导致冬季序列被硬生生截断。很多初学者在这里翻车。本脚本先采用自然年四等分好处是年份边界干净、每季都有完整三个月如果你必须用气候季我建议在聚合时额外处理将12月数据年份加1再按新年份分组这样12月就归到了下一个冬季的年份段里。3.3 按季节聚合sum与mean的选择季节尺度突变检测的关键步骤是把月值压缩为“每年某季节的单值”。降水和径流用总量时选sum气温和面指数用平均值时选mean。聚合代码seasonal df.groupby([season, year])[value].agg(sum).reset_index() print(seasonal[seasonal[season] summer].head())输出示例season year value 0 summer 1961 312.5 1 summer 1962 287.3 2 summer 1963 350.1这里有一个关键参数选择用sum降水的季节总量对缺失月份极度敏感。如果某年夏季只有6月和8月有数据7月缺测那么该年夏季总量会凭空少三分之一聚合后变成一个异常低值M-K突变检测可能把它误判为一次“突变”。因此聚合前应该补做月度完整性检查# 统计每个季节年份下有效数据月数少于3个月的标记为NaN month_count df.groupby([season, year])[month].count().reset_index() valid month_count[month_count[month] 3] seasonal seasonal.merge(valid, on[season, year], howinner)这段代码先把每个月记录数算出来只保留完整包含3个月的季节年份。如果你的序列本身允许季节平均基于2个月比如站点冬季经常缺测12月可以用2个月阈值但必须注明这会降低结果可靠性。4. 用Python计算季节尺度的M-K突变检测UF/UB核心函数4.1 UF/UB核心算法实现这一段是整个脚本的心脏。我按魏凤英《现代气候统计诊断与预测技术》中的经典公式实现用numpy向量化代替三重循环几十年的序列瞬间出结果import numpy as np def calc_uf(seq): 计算M-K正向秩序列统计量UF。 返回数组长度为n索引i对应原始序列第i1个时刻的UF值。 n len(seq) s np.zeros(n) # 第i个时刻统计其前面所有值中比它小的个数 for i in range(1, n): r_i np.sum(seq[:i] seq[i]) s[i] s[i-1] r_i uf np.zeros(n) for k in range(2, n 1): # s[k-1] 就是标准公式中的 s_k e k * (k - 1) / 4.0 var k * (k - 1) * (2 * k 5) / 72.0 uf[k - 1] (s[k - 1] - e) / np.sqrt(var) # uf[0] 0对应UF_1 return uf def calc_ub(seq): 计算反向秩序列统计量UB并与原始时间轴对齐。 rev seq[::-1] uf_rev calc_uf(rev) # 逆序翻转并对齐取负号后反向 ub -uf_rev[::-1] return ub逻辑说明calc_uf中s[i]是累加的秩序列s[k-1]对应标准的s_k。k从2开始是因为k1时期望和方差均为0数学上未定义。np.sum(seq[:i] seq[i])统计的是“前面的值小于当前值”的数量等价于标准式xi xj的个数。对于逆序序列uf_rev计算后ub -uf_rev[::-1]完成两次操作取负号是把反方向统计量转回正向坐标系索引翻转则是把逆序计算的时间映射回原始时间让UB曲线与UF曲线画在同一根时间轴上。4.2 对每个季节批量检测输出突变年份与显著性拿到UF和UB后需要找交点。交点是UF和UB符号从同侧变到异侧的位置程序实现如下def find_change_points(uf, ub, years): 找UF与UB曲线的交点返回(年份, 方向)列表。 方向1表示上升突变-1表示下降突变。 n len(uf) pts [] for i in range(1, n): prev_diff uf[i-1] - ub[i-1] curr_diff uf[i] - ub[i] if prev_diff * curr_diff 0: # 判定突变方向看UF在交点附近的斜率趋势 direction 1 if uf[i] 0 else -1 pts.append((years[i], direction)) return pts说明window里没有额外参数方向用uf[i]的符号近似判断。如果UF从负转正说明突变后序列上升反之下降。严格一点可以比较交点前后两个时刻UF均值符号direction 1 if (uf[i] uf[i-1]) * 0.5 0 else -1我再补充一个限制只有交点处统计量绝对值小于1.96才有资格称为“显著突变候选点”因为超出置信区间意味着趋势已经显著即使交叉也可能只是两条越界曲线的偶然汇合。所以找交点后加一道过滤conf 1.96 filtered_pts [] for yr, direction in pts: idx years.index(yr) if abs((uf[idx] ub[idx]) * 0.5) conf: filtered_pts.append((yr, direction))4.3 输出结果表每个季节的突变年份、方向、显著性把上述函数串起来对四个季节循环seasons [spring, summer, autumn, winter] result_rows [] for season in seasons: sub seasonal[seasonal[season] season].sort_values(year) if len(sub) 10: print(season, 样本量不足10年跳过) continue years sub[year].values seq sub[value].values uf calc_uf(seq) ub calc_ub(seq) pts find_change_points(uf, ub, years) for yr, direction in pts: result_rows.append({ season: season, year: yr, direction: 上升 if direction 0 else 下降 }) result_df pd.DataFrame(result_rows) print(result_df)运行后你会得到一张类似这样的表seasonyeardirectionspring2003上升summer1997下降summer2001下降winter2010上升先别急着把每个交点都当突变点。尤其“夏季”出现1997和2001两个临近交点时需要结合绘图判断哪一个是主突变。两条曲线可能在真实突变附近来回交叉形成一个“交叉带”这时候要选第一个穿透置信区间后长期分离的点。5. 季节尺度M-K突变检测避坑与常见问题排查5.1 序列太短少于10年就出图会怎样现象只有6年数据的季节序列UF和UB曲线像心电图一样上下穿插交点有七八个根本没法看。原因M-K统计量中的方差公式基于渐近正态近似样本量太小时方差估计误差极大统计量分布严重偏离正态。季节尺度白噪声在短序列中也会被误判为趋势。解决每个季节至少保留10年数据最好能达到30年即气候标准期长度。我在代码里显式写了if len(sub) 10: continue如果数据不够建议改用置换检验或Bootstrap确定突变点置信区间而不是直接套M-K。5.2 交点很多但都不在置信区间内现象UF和UB交叉了但交叉点处两条曲线都在±1.96之内看起来像在“游泳”。原因序列整体没有显著趋势突变交点只是随机波动形成的假交叉。这种情况下没有统计学意义上的突变点。解决不要硬找交点。先看UF曲线是否至少存在一段连续越界过程比如连续3年以上UF 1.96再去看UB是否在对应时段交叉。如果UF从未越界结论应写“未检测到显著突变”而不是“存在多次突变可能性”。5.3 季节聚合时跨年冬季被割裂现象采用气候学冬季12-2月定义时某年12月被分到前一年的冬季但反查数据发现那一年只有12月和1月2月缺失冬季序列出现一个异常低值。原因自然年分组把12月归到当前年份与冬季定义冲突。如果不做特殊处理冬季序列会在每年12月与1-2月之间被切开。解决如果你的业务必须用DJF冬季请给12月单独处理在季节映射前把月份为12的记录年份加1然后再按新年份分组。这样2022年12月会进入2023年冬季序列。本脚本默认自然年四等分的spring1-3月不涉及这个问题读者按需修改。5.4 数据缺失导致季节总量虚低现象某年夏季只有一个月有记录聚合后夏季总量只有平时的一半M-K结果在该年份附近出现一个虚假突变点。原因groupby([season,year]).sum()不感知缺测月它只对现有值求和。解决聚合前做完整性计数。我在第3.3节已给出month_count的过滤代码用3个月完整记录约束。对于气温数据如果用的是mean缺失一个月对季节均值影响相对小但同样建议至少保留2个月并在论文或报告中注明。5.5 多个突变点并存时UF/UB交叉多次现象曲线在置信区间内外交叉了三四次结果表里每个季节都有多个年份。原因序列可能存在多阶段突变也可能存在二次趋势叠加导致UF/UB反复穿越。解决优先选择“交叉后UF持续越界超过3年”的交点。具体做法对每个候选交点向后看后续UF曲线是否持续处于越界状态若连续越界段长度最长即为主突变点。其余次要交点单独列出但不要全部打包为突变结论。数据分析报告里可以这样写“主突变为1997年另在2001、2005年出现次级波动。”6. 把季节尺度M-K突变检测做成可复用的脚本参数化、验证与导出6.1 脚本参数化从硬编码到命令行参数当你手里有几十个站点目录时硬编码文件路径会把人逼疯。我把脚本主体包装成可传参的命令行程序核心结构如下import argparse parser argparse.ArgumentParser(description季节尺度M-K突变检测) parser.add_argument(--input, requiredTrue, help输入CSV文件) parser.add_argument(--value-col, defaultvalue, help数值列名) parser.add_argument(--agg, defaultsum, choices[sum, mean], help季节聚合方式) parser.add_argument(--min-years, typeint, default10, help季节序列最小年数) parser.add_argument(--out-prefix, defaultmk_result, help输出文件前缀) args parser.parse_args()我把主流程放在if __name__ __main__:里先读数据再做季节聚合然后循环季节检测最后把结果写成CSV并绘制曲线。这样每个站点只需要换--input参数即可。参数--min-years控制前面说的样本量最短门槛--agg控制聚合方式非常直观。6.2 用模拟序列验证算法正确性我每次换环境或换数据源都会先用一个已知突变点在40年处的模拟序列做自检这一步是防止算法写错还硬套真实数据。验证代码rng np.random.default_rng(42) mock np.r_[rng.normal(10, 1, 40), rng.normal(13, 1, 30)] # 前40年均值10后30年均值13 years_mock np.arange(1970, 2040) uf calc_uf(mock) ub calc_ub(mock) pts find_change_points(uf, ub, years_mock) print(pts) # 预期输出年份在2009-2011之间正常情况会输出(2010, 1)或(2010, -1)方向取决于UF交汇时的正负号组合。如果输出年份离40年基座即2009年偏差超过2年多半是秩序列里相等值处理或者索引映射出了问题。6.3 结果导出把突变年份和曲线保存为CSV与PNG结果落盘是写报告的刚需。我用to_csv保存结果表用matplotlib把四条季节曲线合成一张四子图方便贴进分析报告import matplotlib.pyplot as plt fig, axes plt.subplots(2, 2, figsize(12, 8)) season_list [spring, summer, autumn, winter] for ax, season in zip(axes.flatten(), season_list): sub seasonal[seasonal[season] season].sort_values(year) if len(sub) args.min_years: continue years sub[year].values uf calc_uf(sub[value].values) ub calc_ub(sub[value].values) ax.plot(years, uf, labelUF) ax.plot(years, ub, labelUB) ax.axhline(1.96, colorgray, linestyle--) ax.axhline(-1.96, colorgray, linestyle--) ax.set_title(season) ax.legend() plt.tight_layout() plt.savefig(f{args.out_prefix}_mk_curves.png, dpi300)导出前要确认seasonal是全局变量或者把它作为参数传入绘图函数。这里说明一个常见坑如果季节序列不是按年份排序的UF/UB会乱跳画图前必须sort_values(year)。我吃过大亏有一次没排序曲线像锯齿一样还以为是算法写错了。这套脚本我留着当模板用每次遇到新的站点数据第一件事是跑一遍模拟序列确认代码在当前numpy/scipy版本下正常再套真实数据。判断突变点时也保留“是否越界”的硬约束宁缺毋滥。希望这个思路能帮你在处理四季信号时少走弯路。本文还有配套的精品资源点击获取