简介基于Python的季节尺度M-K突变检测脚本资源面向气候、水文与环境科学领域的分析人员解决带有季节性周期的时间序列中趋势变化与突变点识别问题。Mann-Kendall方法不依赖数据分布结合STL等季节分解手段可有效分离季节成分与趋势成分适用于SPEI等干旱指数的变异诊断。包内含1个Python脚本和1个Excel数据文件.xlsx共2个文件压缩包约11KB数据与代码分离便于替换输入文件后直接复用脚本逻辑。脚本覆盖数据读取、缺失检查、季节分解、M-K检验、突变点定位与结果可视化等完整流程并基于实际SPEI3数据给出可复现的参考案例。已有196人浏览学习适合入门者快速建立季节尺度M-K突变检测的代码框架也可为干旱演变特征分析提供可直接改造的工具基础。1. 为什么季节尺度的M-K检验比全年跑一遍更靠谱一个水文序列实战场景A同学在整理某水文观测站30年的月径流数据时发现一个反直觉的现象直接对全年序列做Mann-Kendall趋势检验Z值只有0.8完全不显著但把数据按春夏秋冬拆开夏季序列的Z值到了2.3还检测出一个明显的突变年份。问题就出在季节性波动上——冬天的剧烈起伏把夏天的趋势信号给“稀释”了。这个基于Python的季节尺度M-K突变检测脚本就是专门为这类场景设计的它先把时间序列按季节切片做分季节的趋势检验再对去季节化后的序列计算UF和UB曲线定位突变发生的大致年份。适合水文、气象、环境质量等领域的研究者以及所有被“全年检验不显著但分季明显”这类问题困扰的人。2. Mann-Kendall原理与季节尺度设计从统计假设到代码结构2.1 M-K检验的核心公式与两种实现方式Mann-Kendall检验是一个非参数秩检验它不要求数据服从正态分布也不受少数异常值影响所以水文气象里用得非常多。核心思想是看序列里“后一个值比前一个值大”的配对数量设序列长度为n统计量S等于所有(ij)中sign(xj - xi)的和。当n大于8时S近似服从正态分布方差Var(S) n(n-1)(2n5)/18于是可以用Z值来判断趋势是否显著。这里的关键是标准M-K检验只给一个整体趋势方向不告诉你趋势什么时候发生了变化所以需要另一个基于顺序和逆序迭代的UF/UB方法来做突变检测。在季节尺度上常见做法有两种这个脚本把两种都实现了。第一种是分组检验按3-5月为春、6-8月为夏、9-11月为秋、12-2月为冬的方法把数据切块每块单独算Z值和p值这样能回答“哪个季节在变”。第二种是去季节化后再做全年突变检测先用多年同季节均值算出季节内平均原值减掉这个平均得到“季节调整残差”再对残差序列做UF/UB突变检测。这么做的好处是把季节性的周期波动先剥掉突变点就不会被虚假的周期起伏干扰。下面这一段是S统计量计算的Python核心项目里简化为一个可复用的函数import numpy as np def mk_s_statistic(x): n len(x) s 0 for i in range(n): for j in range(i 1, n): s np.sign(x[j] - x[i]) return s说明双重循环是O(n^2)n很大时会偏慢但一般水文序列都在几百个点以内实际跑起来没压力。np.sign返回1、-1或0分别代表上升、下降和相等。若想加速可以先对序列做排序用逆序数法把复杂度降到O(n log n)但代码可读性会差一些我只在长序列里才换成那个版本。参数不需要额外设置直接传入一维数值数组即可。注意M-K检验要求数据独立如果序列有明显自相关S的方差会被低估导致虚假显著。这个脚本没有内置去自相关修正后面避坑章节会专门讲。2.2 项目目录结构与主流程整个项目是一个单文件脚本加一个示例数据目录结构如下seasonal_mk/ ├── seasonal_mk.py # 主脚本包含全部函数 ├── sample_data/ │ └── discharge_30y.csv # 示例数据date, value └── output/ ├── season_trend.csv # 分季节趋势检验结果 ├── mutation_points.csv # 突变点列表 └── plot.png # UF/UB曲线图主脚本的逻辑很直白从main()进去后依次调用五个函数load_data()读CSVpreprocess()做插值和季节标记seasonal_trend_test()算每个季节的M-K统计量mutation_test()计算UF/UB序列最后plot_results()把图画出来。没有用面向对象封装因为数据流是单向的用函数更符合“拿到数据跑一遍出结果”的场景。这样的目录设计有一个好处输出目录和输入目录分开二次开发时可以直接替换sample_data不需要改代码。seasonal_mk.py里所有的函数都只依赖pandas、numpy、matplotlib和scipy环境安装成本很低。下面是一个简化的入口调用片段if __name__ __main__: df load_data(sample_data/discharge_30y.csv) df preprocess(df) season_result seasonal_trend_test(df) mutation_result mutation_test(df) plot_results(df, mutation_result)说明这里seasonal_trend_test返回的是DataFrame每一行是一个季节列包含season、n、S、Var、Z、p、trend。mutation_test返回UF和UB两个数组以及它们在置信区间内的交叉点。整个流程控制在5个步骤内方便调试如果结果不对可以从preprocess开始逐步打印中间表。如果想替换自己的数据只需要保证CSV包含两列第一列是日期可以解析成时间第二列是观测值。日期格式不限于YYYY-MM-DD只要是pandas能识别的格式就行例如2023/1/31、2023-01-31。如果数据里还有站点编号、经纬度这类附加字段脚本当前会直接忽略不影响主流程运行。再看主流程中几个函数的职责划分load_data只负责读文件preprocess只负责清洗和分组seasonal_trend_test和mutation_test只做统计计算plot_results只做可视化。这样拆分的最大好处是当你觉得结果不对时可以单独在某一步后面加print看中间值不需要运行完整脚本。对于新手我建议先跑通示例数据再替换自己的文件这样能更快定位是数据问题还是代码问题。3. 核心代码走读数据清洗、季节分组与突变点计算3.1 数据加载与缺失值处理打开CSV后第一件事是解析日期。示例数据只有两列但真实数据里常有“年、月、日分开”的情况所以load_data里用pd.to_datetime统一解析并生成一个season列。缺失值直接用线性插值补不做前后填充因为对于趋势检验来说前后填充会人为制造长段相同值让S统计量偏低。import pandas as pd def preprocess(df): df[date] pd.to_datetime(df[date]) df[year] df[date].dt.year df[month] df[date].dt.month df[season] df[month].map( {12: 1, 1: 1, 2: 1, # 冬12、1、2月 3: 2, 4: 2, 5: 2, # 春3、4、5月 6: 3, 7: 3, 8: 3, # 夏6、7、8月 9: 4, 10: 4, 11: 4} # 秋9、10、11月 ) df[value] df[value].interpolate(methodlinear) return df.dropna(subset[value])说明季节划分我用的气象常规定义注意把12月归到冬、1月和2月连在一起这样冬季是一个连续跨年的集合分组检验时不会把同一个冬天的数据切成两半。interpolate是pandas自带方法默认按索引顺序线性插值如果缺失值出现在首尾插值会补不上所以最后用dropna去掉这些头尾残缺行。参数上df[value]只要保证是数值类型否则插值会报错建议读入时用pd.to_numeric(df[value], errorscoerce)。这里有一个细节值得展开插值方法的选择会影响结果。线性插值适合缓变的水文数据如果换成多项式插值很容易在缺失段附近产生振荡制造出虚假的上升/下降趋势。我一般只在连续缺失不超过5个点时才用线性插值超过5个点的直接删掉那一整段或者用季节平均填充。因为趋势检验对局部大块缺失非常敏感硬补的结果会让S统计量出现大量0差配对。加载完数据后最好先看一眼季节分布。如果某个季节只有两三年数据M-K检验会因为没有足够的配对而无法计算方差。这时脚本会在seasonal_trend_test里返回Z0, p1表示“不显著”但实际原因是数据量不足而非真的没有趋势。为了避免误导我通常会在预处理后加一条打印语句print(df.groupby(season).size())检查每个季节的样本量。3.2 季节分组与M-K统计量计算有了season列后分组检验就很简单。按季节分组每组重新编序调用前面写的mk_s_statistic然后用正态近似算出Z值。这里有一个容易忽略的点如果某个季节有并列值方差公式要减掉结(tie)修正项否则Var会被高估。from scipy import stats as st def seasonal_trend_test(df): results [] for season, group in df.groupby(season): x group[value].values n len(x) s mk_s_statistic(x) unique, counts np.unique(x, return_countsTrue) ties sum(c * (c - 1) * (2 * c 5) for c in counts if c 1) / 18 var_s n * (n - 1) * (2 * n 5) / 18 - ties if var_s 0: results.append({season: season, n: n, S: s, Z: 0.0, p: 1.0}) continue z (s - 1) / np.sqrt(var_s) if s 0 else ((s 1) / np.sqrt(var_s) if s 0 else 0) p 2 * (1 - st.norm.cdf(abs(z))) results.append({season: season, n: n, S: s, Z: z, p: p}) return pd.DataFrame(results)逻辑说明unique统计每个值出现的次数ties是标准结修正项注意这里修正项是乘以c(c-1)(2c5)/18而不是直接用c因为每个结内部会产生重复比较S的方差需要扣掉这些重复信息。当var_s小于等于0时说明序列太短或全是相同值直接把Z设0、p设1是稳妥的兜底做法。st.norm.cdf是标准正态累积分布函数双侧p值等于2倍右尾概率。参数说明season是整数1-4对应冬春夏秋n是季节内样本数比如月数据一个季节有30年*3个月90个点S是原始秩统计量Z是标准化后的值正为上升趋势负为下降趋势。判断显著性的阈值通常取1.9695%置信水平。如果使用月数据每个季节有3个月季节内部其实是连续的月度值所以样本量足够大检验结果比较稳定。实际运行时我通常会把results排序后打印成表格。季节趋势的意义往往比全年结果更直观比如“冬季无显著变化”、“夏季显著上升”这样一句话比单一Z值容易理解得多。代码输出的trend列会在后面填上increasing、decreasing或no trend这个字符串是根据Z值和p值组合判断的。若是p小于alpha且Z为正则标记increasing反之同理。3.3 突变点定位UF与UB序列的交叉判定突变检测用的是Forward和Backward两个序列UF从序列头开始逐步计算M-K统计量UB从序列尾开始反着算。如果UF和UB两条曲线在95%置信区间内出现交叉交叉点对应的位置就是突变发生的时段。实现的时候注意UF[t]需要基于“前t个样本”重新计算S和Var而不是直接复用整条序列的结果。def compute_uf(x): n len(x) uf np.zeros(n) for t in range(1, n): s mk_s_statistic(x[:t 1]) unique, counts np.unique(x[:t 1], return_countsTrue) ties sum(c * (c - 1) * (2 * c 5) for c in counts if c 1) / 18 var_s t * (t 1) * (2 * t 6) / 18 - ties if var_s 0: uf[t] 0 else: uf[t] (s - 1) / np.sqrt(var_s) if s 0 else ((s 1) / np.sqrt(var_s) if s 0 else 0) return uf def mutation_test(df): x df[value].values uf compute_uf(x) ub_rev compute_uf(x[::-1]) ub ub_rev[::-1] critical 1.96 cross_idx [] for i in range(1, len(uf)): if (uf[i-1] - ub[i-1]) * (uf[i] - ub[i]) 0: if abs(uf[i]) critical and abs(ub[i]) critical: cross_idx.append(df[date].iloc[i]) return {uf: uf, ub: ub, cross_points: cross_idx}逻辑说明compute_uf内部对每个t都调用mk_s_statistic复杂度是O(n^3)在数据点超过500后会很慢项目里默认数据规模不大如果未来接入逐日数据需要改用增量更新S的方式后面二次开发章会提到。cross_idx判断的是UF和UB符号相反且差值变号的位置并且要求两条曲线都在临界值内避免在曲线都跑出置信区间时还误判。交叉点对应的是df[date].iloc[i]也就是突变发生的具体日期或年份。参数说明critical1.96对应95%置信水平如果你想更严格可以改成2.5899%但一般水文文献用95%较多。cross_idx里存的是所有满足条件的交叉时刻实际突变通常取第一个或最后一个需要结合原始序列图人工确认。如果交叉点有多个我一般用第一个交叉点作为突变开始年份然后结合趋势图看是否跟肉眼判断一致。一个常见的困惑是UF和UB的定义方向。UF是正序列迭代结果UB是逆序列迭代结果。当序列存在上升趋势时UF会从零开始一路上升超过临界值UB则从反向来看通常一开始就在临界值附近最后回到零。两条线的交叉点往往出现在趋势方向改变的年份。如果你画的图里两条线完全平行说明序列没有突变只有单调趋势如果交叉点在图形边缘则可能只是数据开头或结尾的边界效应不应解读为真实突变。4. 运行与结果解读跑通脚本并看懂输出4.1 参数配置与命令行入口脚本用argparse接收参数不传时会走默认的示例数据。常用参数有三个--input指定CSV路径--alpha设置显著性水平--output指定输出目录。跑起来的命令很简单python seasonal_mk.py --input sample_data/discharge_30y.csv --alpha 0.05 --output output说明--alpha是判断趋势显著性的阈值脚本内部会转成Z值的双尾临界值默认0.05就是1.96。--output目录不存在时脚本会自动创建避免手动mkdir。如果你想把图画得更大可以额外加--fig-size参数项目里默认宽度是12英寸高度是6英寸。运行前需要确认Python环境里安装了pandas、numpy、scipy和matplotlib。这一步常被新手忽略报错时往往是一个ModuleNotFoundError。我一般用下面这条命令安装pip install pandas numpy scipy matplotlib安装完成后把命令行切到seasonal_mk目录下再运行脚本。如果CSV路径写错了脚本会在load_data阶段崩溃提示FileNotFoundError。解决方法是先看当前工作目录再写相对路径最稳妥的是直接使用绝对路径。运行后终端会打印两段信息第一段是每个季节的Z值和p值第二段是突变点年份。如果某个季节的样本数少于3脚本会跳过并给出warning因为M-K检验最少需要3个数据点才能算出S统计量。--alpha传0.01时临界值变成2.58传0.1时变成1.64这个参数只影响突变检测的置信区间宽度不影响分季节趋势检验的p值计算。4.2 结果输出与显著性判断结果文件里最需要注意的是season_trend.csv中的Z和p两列。我的习惯是先看p值p0.05再看Z的正负。比如示例数据跑出来夏季Z2.3p0.021趋势标记为increasing冬季Z-1.2p0.23标记为no trend。注意p0.05不完全是硬边界如果样本量很小p值本身就不可靠这时我会看S的正负作为辅助。season_trend.csv会保存成下面这种结构season,n,S,Var,Z,p,trend 1,90,-45,1234,-1.2,0.23,no trend 2,90,12,1300,0.8,0.42,no trend 3,90,220,1280,2.3,0.021,increasing 4,90,-30,1250,-0.9,0.36,no trend说明这个文件可以直接用Excel打开如果你不喜欢终端里对齐的输出。trend列根据Z的正负和p值生成但要注意p值只是显著性判断不直接代表趋势大小。两个季节可能p值相近但Z值一个0.2一个2.0说明样本量或方差不同不能只对比p值。突变点输出在mutation_points.csv每条记录包含date和uf_ub_cross两个字段。这里的日期是原始序列中交叉点所在的那一天不是突变开始的准确时刻通常可以把这一年当作趋势突变的分界年。举个例子如果交叉点在2008年那说明2008年前后序列的统计特征发生了显著变化。需要特别说明的是UF/UB方法只能定位“大概从哪年开始变化”不能告诉你变化幅度幅度要用后续章节提到的Theil-Sen斜率来补充。看输出图时有一种比较容易误判的情况UF和UB在置信区间外交叉比如两条线都大于1.96。这时候交叉点不能算突变点因为序列的统计特征已经发生了极端变化交叉位置不稳定。脚本里特意把这种情况过滤掉了所以图里只画区间内的交叉标记点。如果你看到图中UF和UB交叉但脚本没有报告突变年份多半是交叉点跑到了置信区间外这是正常现象不是bug。5. 常见问题与避坑指南季节M-K检验中的五个典型坑5.1 数据层面的三个坑第一个坑是季节划分不连续。有人把冬季定义成“1、2、3月”结果3月的序列和后面春季重叠分组时同一个日期出现在两个季节里。解决方法是先把月份映射成季节码再按季节码分组不要用两个分隔条件分别筛选。示例代码里的map字典就是标准写死的方式改动季节定义时只改字典即可。第二个坑是闰年带来的2月29日。如果你的数据是逐日序列闰年会多出一天导致某一年的冬季样本数比其他年份多1。这在突变检测里影响不大但在按季节汇总均值时会引入微小偏差。解决办法是在预处理时主动过滤掉2月29日或者把日期先归一化到年份-月格式。我常用的是在preprocess里加一行df df[~(df[date].dt.month 2) | ~(df[date].dt.day 29)]简单直接。第三个坑是缺失值补过头。线性插值适合中间缺失但如果你用ffill前向填充遇到连续缺失十天以上的片段会产生一段平台期S统计量里会出现大量0差配对直接压低了突变信号的灵敏度。遇到连续缺失超过总长度5%时我一般会放弃插值直接删掉这段连续缺失的前后数据而不是硬补。5.2 算法层面的两个坑第四个坑是忽略结修正。实际数据里经常出现相同的观测值比如流量数据取整后有很多1、2这样的重复值。如果你直接用没有结修正的方差公式S的方差会偏大Z值偏小趋势会被掩盖。项目里已经加了unique/counts修正但如果你复制代码到别处很容易把这十几行删掉。判断有没有修正算一下Var(S)是否总等于n(n-1)(2n5)/18如果总等说明没做结修正。第五个坑是自相关导致的高估显著性。M-K检验假设样本独立但水文序列基本都有明显的1阶自相关比如今天的径流受昨天影响。如果序列有自相关而你不处理UF/UB曲线的交叉点会变得特别多看起来像有很多突变实际上都是虚假信号。轻量做法是在做M-K检验前先对序列做一次滞后1阶回归取残差作为新序列严格做法是计算有效样本量n*做方差修正。这个脚本版本为了保持简洁没有内置自相关处理二次开发时我在下一章给出一个简单的补救函数。这里再补充一个不是所有人都遇到过的问题时间轴排序错误。CSV里的日期如果不是升序排列mk_s_statistic会把排序错误的配对也当成真实比较结果毫无意义。解决方法是读入后立刻执行df.sort_values(date, inplaceTrue)。这个操作成本极低但能避免一整晚的困惑。还有如果你的数据是多年逐日序列但中间跳过了某些整段年份M-K检验仍然会把这些年份之间的间隔当成连续时间导致突变点偏移。遇到这种情况我会在缺失年份处插入断点序列或者按年份分组后只计算年际变异而不是直接跑逐日突变。6. 二次开发把这个脚本扩展成通用的季节检验工具拿到这个脚本后最常见的改造需求是把它从“只能跑季节”扩展成“能跑月、旬、候等多种时间尺度”。我一般会在preprocess中增加一个--scale参数值为season、month、quarter或dekad然后动态生成分组列。下面是一个最小改法def set_scale(df, scale): if scale season: df[group] df[month].map({12:1, 1:1, 2:1, 3:2, 4:2, 5:2, 6:3, 7:3, 8:3, 9:4, 10:4, 11:4}) elif scale month: df[group] df[month] elif scale quarter: df[group] df[month].map({1:1, 2:1, 3:1, 4:2, 5:2, 6:2, 7:3, 8:3, 9:3, 10:4, 11:4, 12:4}) return df这样一来seasonal_trend_test里的groupby(season)改成groupby(group)就能复用。我还建议把mk_s_statistic改成增量更新版本避免O(n^3)在逐日数据上跑不动。具体做法是维护一个前缀数组每次新来一个样本时用当前前缀里的较小值个数更新S复杂度降到O(n^2)整体跑得快很多。另一个有价值的扩展是加一个Theil-Sen斜率估计。M-K只能告诉你趋势是否显著斜率大小才是实际关心的。对每个季节的序列计算所有配对点的斜率中位数这个值不受异常值影响配合M-K的p值可以写进报告“夏季显著上升斜率上升0.12个单位/年”。实现代码只有十几行但能让结果的专业度提升一个档次。验证改造是否成功的方法很简单拿同一份数据分别跑原版和改造版对比p值和交叉点年份差异应小于0.01。如果你还想要一个自相关修正函数可以写成这样def prewhiten(x): # 滞后1阶回归取残差替代原序列 x np.asarray(x, dtypefloat) lag1 x[:-1] y x[1:] beta np.polyfit(lag1, y, 1)[0] residual y - beta * lag1 return np.r_[x[0], residual]这个函数不是万能药但对于1阶自相关较强的序列能让UF/UB曲线变得平滑不少交叉点数量也会回归正常水平。用了预白化之后原本跑出来五六个突变点的序列通常能减少到一两个真实突变点。从那以后我每次拿到新序列都会先跑一遍分季节的趋势表再看全年突变点两个结果结合起来定结论改代码时也强制自己走一遍“原版对照”这个流程。希望帮到你。本文还有配套的精品资源点击获取