★ FEATURED ARTICLE 简介基于Python实现季节尺度曼肯德尔突变检测的完整脚本面向气候与环境数据分析者也适用于水文、生态等时间序列突变研究。脚本利用曼肯德尔非参数检验法识别时间序列中的显著趋势突变点该方法不依赖数据分布对非正态数据依然稳健同时结合基于Loess的季节趋势分解处理周期性波动可对标准化降水蒸散发指数逐季节进行突变检测并自动标记突变发生年份与显著性水平。资源共包含两个文件一个为可执行的Python主程序一个为配套的Excel示例数据表压缩包大小仅为11KB非常轻量。程序完整覆盖数据读取、缺失值处理、季节分解、突变统计量计算及结果可视化可直接运行验证帮助理解季节尺度突变检测的原理与编码实现。已有195人学习使用适合需要快速开展气候突变分析的研究者与学生。1. 季节尺度 M-K 突变检测Python 实现思路与代码资源拆解做气象水文序列分析的人基本都绕不开 Mann-Kendall 检验。全网搜 python M-K 突变检测教程不少但绝大多数只教你处理全年序列一说到季节尺度要么含糊带过要么直接把四季数据拼在一起跑——那结果就是突变点分布成了玄学UF 和 UB 曲线交叉得像一团乱麻没法跟实际气候背景对上。这份基于 Python 实现季节尺度 M-K 突变检测的代码资源解决的核心问题就是把序列按季节拆分分别计算统计量再叠合展示让每个季节的突变位置独立且可解释。它面向的是我这类需要跟降水、气温、径流序列打交道的气象水文从业者也适合写论文需要汇报突变年份的学生。结论先行季节性拆分不是简单分组正序反序序列的定义窗口、交点判读规则、季节样本量不足时的处理策略才是这个代码真正值钱的地方。2. 为什么季节尺度要单独处理M-K 方法在分组序列上的边界2.1 M-K 检验的数学前提与季节分组的冲突先回顾经典 Mann-Kendall 检验的原理。它本质上是基于秩的非参数检验通过比较序列中所有元素对的大小关系来构造 Kendall 秩相关系数。对时间序列 xtt 1, 2, ..., n统计量 S 定义为S Σ(ij) sign(xj - xi)其中 sign 函数只在 xj xi 时取 1两者相等取 0xj xi 时取 -1。S 的期望值为 0方差在序列不存在趋势的假设下可以用 n(n-1)(2n5)/18 近似。这地方有个容易被忽略的前提M-K 检验假设数据是独立同分布的并且序列中不存在季节周期性。如果直接把 30 年逐月数据混在一起跑周期性会污染秩比较导致 S 统计量被季节性波动主导趋势信号被淹没。我做某流域降水序列分析时犯过这个错全序列 M-K 检测显示显著上升趋势按季节拆开看夏季降水根本无显著趋势趋势完全来自春季那一支。季节尺度的做法是把序列按自然季节通常 3 个月一组切成四个子序列每个子序列分别做 M-K 检验。这个拆分的意义不只是统计口径更关键是突变检验使用正序和逆序两个方向累积计算当子序列长度只有原序列 1/4 时正序、逆序序列的构建方式直接影响交点数量和位置。2.2 季节尺度的输入数据格式约定这份代码对输入数据有明确的格式预期我在拆解时整理了一份最小样例字段类型说明示例yearint年份必须连续或不含重复1991monthint1 到 12 的整数月份6valuefloat观测值降水/气温/流量128.4注意代码不做缺失年份补齐。如果你的数据有断年比如 1995 年缺一整年最简单的处理是先把该年份置为 NaN再调用代码的预处理函数做线性插值如果缺的年份太多超过总年份的 10%建议放弃该站点的季节突变分析。我一般习惯把数据整理成 CSV 格式三列无表头或带表头均可代码里用 pandas read_csv 读取后统一把列名映射为 year、month、value。读取后第一件事是检查月份范围# 检查月份合法性 if not ((raw_data[month] 1) (raw_data[month] 12)).all(): raise ValueError(月份必须为 1 到 12 之间的整数)这段预校验很有必要我见过一份数据里月份列混入了 13 和 0结果按季节分组时1 月的数据归到冬季13 月的数据被丢进 NaN 组直接影响后续统计量的秩计算。2.3 季节划分边界固定季节还是滑动窗口代码默认按固定自然季节分组春季MAM3 月到 5 月夏季JJA6 月到 8 月秋季SON9 月到 11 月冬季DJF12 月到次年 2 月需要注意12 月属于冬季但 12 月来自上一年。代码在实现时把 year 和 month 拼接成时间索引然后通过dt.quarter直接映射季度但这里有个隐性约定冬季是从当年 12 月开始的。比如 2020 年 12 月、2021 年 1 月和 2 月会被归为同一个冬季组跨了年份。代码内部处理方式是把 12 月的年份减 1统一到冬季结束年份这样每个季节组内的年份是连续且对齐的。如果你要的是水文学意义上的水文年季节划分比如从 10 月到次年 9 月这个代码不直接支持。我改过一版把季节映射函数改成可传入偏移参数的模式def assign_season(month, offset0): # offset0 为自然年季节offset1 表示从 10 月开始划分水文年季节 adjusted_month month - offset if adjusted_month 0: adjusted_month 12 ...对降水序列来说水文年划分往往比自然年更合理但这个改动会牵涉到后面 UF 序列的累积窗口构建不仅仅是分组函数的问题建议只在数据量充足每季 30 年以上时尝试。3. 核心算法实现UF、UB 序列如何正确计算3.1 正序累积统计量 UF 的逐季算法M-K 突变检测的核心输出是两个序列UF正序累积统计量和 UB逆序累积统计量。两者的交叉点就是突变发生的候选位置。这个代码的实现思路是在每个季节的子序列内逐样本递增地计算累积统计量。UF 的递推公式是这样走的对季节内子序列 x1, x2, ..., xm在第 k 步计算S_k Σ(i1..k-1) Σ(ji1..k) sign(xj - xi)这个双层累加在代码里不是直接暴力循环而是利用 Python 的itertools.combinations把前 k 个元素的所有组合生成出来再统计时间复杂度 O(k^2)。对于 30 年序列一个季节只有 30 个数据点组合数 C(30,2) 435暴力计算完全可接受。我用 70 年序列试过单季运行时间也就毫秒级不需要优化成 Fenwick 树。关键代码片段如下def compute_uf_subseries(sub_series): 计算正序累积统计量 UF 序列 n len(sub_series) uf [] s_k 0.0 var_s_k 0.0 for k in range(2, n 1): # 所有前 k 个元素的组合 comb list(itertools.combinations(range(k), 2)) sign_sum 0 for i, j in comb: diff sub_series[j] - sub_series[i] if diff 0: sign_sum 1 elif diff 0: sign_sum - 1 s_k sign_sum var_s_k k * (k - 1) * (2 * k 5) / 18.0 if var_s_k 0: uf.append(0.0) else: if s_k 0: uf.append((s_k - 1) / math.sqrt(var_s_k)) elif s_k 0: uf.append((s_k 1) / math.sqrt(var_s_k)) else: uf.append(0.0) return uf逻辑说明这里循环从 k 2 开始因为前 1 个元素无法构成元素对。s_k是第 k 步的累积 S 值var_s_k是第 k 步的方差。分子加 1 或减 1 是连续性修正这是 M-K 突变检测的标准做法不加修正会导致 UF 值在临界区附近整体偏小交点位置容易滞后。uf列表的索引 0 对应第 2 个数据点的统计量在后续画图中要注意这个偏移。参数说明sub_series是无缺失值的季节子序列。代码要求每个季节子序列的缺失值已在前置处理中插补否则会出现组合索引越界。我自己在实测中发现如果子序列里有 NaNdiff 0和diff 0的判断都会失效比较结果恒为 False最终 sign_sum 是 0UF 序列会平白多出一段 0 值干扰交点位置判读。3.2 逆序 UB 序列与序列反转的陷阱UB 序列的计算是这份代码最容易踩坑的部分。教科书上的标准做法是把序列反转对反转序列重新做一遍正序 M-K 累积再对结果的 UF 序列做反向翻转和平移得到 UB。这个描述说三句话就完了但实现细节差一毫米结果差千里。代码中的实现方式是def compute_ub_from_reversed(uf_reversed): 从反转序列的 UF 结果还原 UB 序列 n len(uf_reversed) ub [0.0] * n for i in range(n): ub[i] -uf_reversed[n - 1 - i] return ub这里的关键坑在于反转序列的 UF 计算的起始端是原序列的末端。举例说原序列是 1991 到 2020 年 30 个值反转后第一个值是 2020 年那么 UF 的第一个值反映的是 2020 年与之前年份的比较而不是 2020 年与之后年份的比较——因为反转后没有“之后”了。正确的解读是UB 在第 k 个点的值代表从末端反向累积到第 k 个点时的统计显著性。我的实践体会是调试阶段最有效的验证方式是用一个已知突变的人造序列来检查 UB 曲线的形态。比如构造一段前 15 年均值 50、后 15 年均值 80 的序列理想的 UF 曲线在突变点前缓慢攀升或波动在突变点后快速上升突破临界线UB 曲线在突变点后向下突破。如果你的 UB 曲线形态完全反过来说明反转时把-uf_reversed[i]的符号或索引弄错了。这段逻辑放在独立函数里是有意的代码作者在注释里专门标注了这一点可见这也是他踩过坑之后的总结。3.3 季节子序列的样本量与 UF 起点偏移季节尺度处理的一个特有问题是样本量。总序列 30 年时每个季节只有 30 个样本UF 序列从第 2 个点开始计算。如果数据从 1991 年到现在只有 20 年那么春季子序列只有 20 个点UF 序列前段前 5 年左右因为统计量累计不足方差小UF 值波动特别剧烈经常会在没有突变的位置打出虚假交点。我建议在做季节 M-K 突变检测前先对每个季节的子序列长度做检查if len(sub_series) 15: print(f警告{season}季样本量 {len(sub_series)} 不足 15突变点结果仅供参考)15 是经验阈值。教科书上对 M-K 检验的样本量要求一般建议 n ≥ 10考虑到季节子序列的独立样本更少我一般把红线设在 12低于 12 直接放弃该季节的突变判断只报告趋势方向。这不是原代码的默认行为但我强烈建议在调用时加这个判断——尤其在处理高原站点、极地区域站点时很多数据序列本身就只有十几年。还有一个容易忽略的起点偏移问题。UF 序列的第一个有效值对应原始序列的第二个数据点。在绘制 UF/UB 曲线时横坐标应该是该数据点对应的年份。比如春季序列包含 1991 到 2020 年的 3 个值第一个 UF 值对应 1992 年。代码里在处理横坐标时用season_years[1:]来对齐如果你直接用 0 到 n-1 的整数索引画出来的突变点年份会比实际早一年。4. 突变点判读与可视化交点到结论的完整链路4.1 交点判读逻辑UF 与 UB 交叉的三种情况完成了 UF 和 UB 序列的计算真正的挑战才开始怎么判断突变点这套代码内置了三条判断规则我会逐一拆开说。第一条规则是基础过滤只保留 UF 与 UB 曲线在置信区间内的交点。置信区间取的是显著性水平 α 0.05 对应的 ±1.96。两条曲线可能在区间外频繁交叉但这种交叉没有统计意义。代码里先用布尔索引筛出 UF 和 UB 同时处于 [-1.96, 1.96] 之内的区间再在这个区间段内搜索交点。第二条规则处理同一季节内多个交点的情况。气象序列现实里很常见UF 和 UB 在一个季节里交叉四五次。代码的默认策略是取连续交叉段中最靠前的交点作为突变候选点但因为 UF 序列前段方差小、波动大最靠前的交点往往是个假信号。我一般把代码改成只考察 UF 与 UB 在临界线内侧交叉、且之后 UF 持续上升或下降超过 1.96 的那一段取该段的起始交叉位置。这样过滤后多个交点通常能收敛到 1 到 2 个有效候选。第三条规则是方向约束突变点必须与趋势方向一致。假设春季序列整体是升温趋势UF 曲线最终在临界线上方那么突变点在 UF 从负转正并突破 0 线的位置附近而不是 UF 在负值区域交叉的位置。代码里会在候选点列表中检查如果 UF 终点值大于 1.96则只保留候选点中 UF 从负变正的那个索引反之亦然。规则落实到代码是这样的def locate_mutation_point(uf, ub, years, conf_level1.96): candidates [] in_conf (np.abs(uf) conf_level) (np.abs(ub) conf_level) for i in range(1, len(uf)): # 检测 UF 与 UB 的符号交叉 if np.sign(uf[i] - ub[i]) ! np.sign(uf[i-1] - ub[i-1]): if in_conf[i] and in_conf[i-1]: candidates.append(i) if len(candidates) 1: return candidates[0], years[candidates[0]] if len(candidates) 1: # 取与最终趋势方向一致的候选 trend_dir 1 if uf[-1] 0 else -1 for idx in candidates: if trend_dir 1 and uf[idx] 0 and uf[idx-1] 0: return idx, years[idx] if trend_dir -1 and uf[idx] 0 and uf[idx-1] 0: return idx, years[idx] return candidates[0], years[candidates[0]] return None, None逻辑说明np.sign的比较本质上是在检查两条曲线的相对位置是否发生翻转。in_conf确保交叉点两端的统计量都在临界线内过滤掉置信区间外的交叉。trend_dir由 UF 终点方向决定这与你的气候假设有关——我在实际使用中遇到过春季降水序列 UF 终点稍微大于 0 但整条曲线来回穿越此时直接取第一个 UF 变号位置可信度不高建议配合趋势显著性p 值来判断。参数说明conf_level1.96是硬编码的 95% 置信度。如果样本量很小小于 20我建议把置信水平放宽到 90% 即 1.645这样能避免因为临界线过宽而把所有交点都滤掉。4.2 交叉点的可视化呈现与年份标注这套代码的绘图部分用 matplotlib 实现核心是双曲线加置信线以及突变点的散点标注。我改过一版更适合论文使用的画法效果很直观def plot_mk_seasonal(results, season_name, output_path): fig, ax plt.subplots(figsize(8, 5)) years results[years] uf results[uf] ub results[ub] ax.plot(years, uf, labelUF, color#1f77b4, linewidth1.8) ax.plot(years, ub, labelUB, color#ff7f0e, linewidth1.8) ax.axhline(1.96, linestyle--, colorgray, linewidth0.8) ax.axhline(-1.96, linestyle--, colorgray, linewidth0.8) ax.axhline(0, linestyle:, colorblack, linewidth0.6) # 标注突变点年份 if results[mutation_year] is not None: ax.axvline(results[mutation_year], colorred, linestyle-., linewidth0.9) ax.annotate(f突变年份: {results[mutation_year]}, xy(results[mutation_year], ax.get_ylim()[0]), xytext(results[mutation_year], ax.get_ylim()[0] 0.3), rotation90, fontsize9, colorred) ax.set_title(f{season_name}季 M-K 突变检测, fontsize13) ax.legend(locbest) fig.tight_layout() fig.savefig(output_path, dpi150) plt.close(fig)这个绘图函数有几个细节值得说明。一是红色虚竖线标注的突变年份应该从results[mutation_year]取而不是从交点索引直接取因为年份可能不连续如数据有缺失年份时二是画图前必须确保 UF 和 UB 的年份轴对齐否则曲线错位会误导判读三是保存时用dpi150期刊插图一般要求 300 dpi投稿前记得改。我实际用的版本里还加了一个功能把四个季节的图拼成 2x2 的大图共享一个图例。这个做法在汇报场景下很实用审阅人一眼就能看到四个季节突变年份的差异。拼图时注意每个子图单独设置标题颜色统一按季节分配——我习惯春季用绿色、夏季用红色、秋季用橙色、冬季用蓝色跟季节属性匹配方便非专业读者快速关联。4.3 显著性判定与结果输出格式代码最终输出的结果除了图还有一张结构化表格字段包括季节、突变年份、UF 终值、是否显著、趋势方向。这个输出我改造过一版增加了置信区间的上下界方便直接写进论文的方法部分季节突变年份UF 终值0.05 显著性趋势方向春季19982.84显著上升夏季20061.32不显著上升秋季无0.76不显著无明显趋势冬季2003-2.21显著下降这里有一个方法论上的坑我当时在写某高校的项目报告时被导师指出过M-K 突变检验的结果本质上是筛选与定位不具有假设检验的严格 p 值基础。UF/UB 的显著性判断依赖于每一步累积统计量的独立检验序列内部存在自相关时方差公式的假设会被打破。因此季节尺度的 M-K 结果适合作为突变年份的初筛证据但要在论文里拿到硬结论还需要配合其他方法验证这点我在第 6 章专门展开。5. 季节尺度 M-K 的高频翻车点排查5.1 现象UB 曲线走势与 UF 完全相同有次我拿一个 1961-2020 年降水序列跑代码输出的春季 UB 曲线跟 UF 曲线几乎重叠季节突变年份一下子冒出四五个。排查后发现原因在数据组织读取 CSV 时月份列被 pandas 读成了字符串排序时字符串字典序导致 1 月在 10 月前面季节分组完全错位序对比较落在了错误的元素对之间。解决方法是强制转换类型并检验月份唯一性raw_data raw_data.astype({month: int, year: int})从那以后我每次读取数据都会跑一遍data.dtypes检查不再信任 CSV 里看似正常的数字列。5.2 现象所有季节突变点都集中在序列首段冬季序列的 UF 曲线在第一二个点就冲到 2.5 以上交点全部被过滤输出“无突变”。原因很清楚冬季样本量少且极端值集中前几个点的秩比较受一个极端值主导方差还没累积起来UF 值虚高。解决策略是增加一个最小计算窗口。我把compute_uf_subseries里的循环起始值从 2 改成 6即前 5 个点只做累积不参与突变判读这样首段虚高的 UF 值不会生成候选点。注意这只是工程回避如果极端值是真实物理过程要保留在序列里不要为了平滑而删除。5.3 现象交点年份与已知气候突变年份不一致某流域夏季降水序列代码给出的突变年份是 1994但公开发表文献一致指向 1992 年区域环流调整。检查发现问题出在插补策略——我把 1992 年上半年连续的缺失值用线性插值填充插值造成序列中间段的秩比较信息失真。M-K 检验对秩敏感插值序列与真实观测的秩关系不可靠。解决方式是改用季节均值填充。对降水这种强季节序列用该季节多年气候平均值替换缺失值对秩结构的影响远小于线性插值。我改成了sub_series sub_series.fillna(sub_series.mean())这方法不完美但至少不会引入趋势性偏差。5.4 现象代码运行无报错但输出结果明显偏离理论预期冬季 M-K 结果显著上升而我人工绘制距平曲线看到的却是明显下降。此事一度让我怀疑算法实现有问题后来发现是年份排序方向问题。CSV 中数据按年份降序排列代码内部做sort_values(year)时默认升序但冬季跨年份的归属处理在排序前执行导致 2020 年 12 月与 2020 年 1 月被归入同一季节从而把最末端的值错误地带入了序列开头。解决方法是严格规定数据格式原始文件必须按年升序保存若使用水文年季节划分则先处理跨年归属再排序。这个坑我花了整整一个下午定位最后打了一行print(season_years[:10])才发现排序错乱。5.5 现象不同季节的突变年份彼此矛盾无法解释四个季节最终输出春季 1994、夏季 2003、秋季 无、冬季 1998。这种情况在真实数据中非常常见不代表代码有 bug而是气候系统的多尺度响应本身存在季节差异。我处理方式是不强行统一而是把四个季节的突变年份连同相应季节的趋势显著性写成对比矩阵。若论文要求给出一个总突变年份可以采用季节面积加权或最主要季节降水贡献率加权但必须在方法部分说明加权依据。6. 交叉验证技巧用 Pettitt 检验与滑动 t 检验增强突变结论可信度季节尺度 M-K 突变点往往不止一个候选位置特别是当 UF 与 UB 在临界线附近反复缠绕时算法选出的交点年份需要外部验证。我常用的验证手法是 Pettitt 检验它本质上是非参数的变点检测对单一变点的定位能力比 M-K 的目视判读更严格适合作为 M-K 结果的交叉印证。Pettitt 检验代码简短我在 M-K 跑出候选突变年份后会再跑一遍这个def pettitt_test(data): n len(data) ranks pd.Series(data).rank().values U 0.0 max_U 0.0 k_star 0 for k in range(n - 1): # 计算候选变点前后累积等级差 U np.sum(ranks[:k1]) - k * (k 1) / 2.0 if abs(U) abs(max_U): max_U U k_star k # 显著性近似公式 p_value 2.0 * math.exp(-6.0 * max_U**2 / (n**3 n**2)) return k_star, p_value逻辑说明这个实现利用了一种递推关系将原 Pettitt 的双重循环化简为单次遍历。U 值绝对值最大点对应的 k_star 就是变点位置。p 值的公式是渐近近似当 n 超过 50 时可用性较好n 在 20 到 50 之间时我通常用蒙特卡洛模拟获取更可靠 p 值。参数说明data是某一季节的年值序列。注意输入顺序必须与 M-K 分析时的季节子序列一致否则交叉验证就失效了。对应关系上Pettitt 给出的 k_star 是基于 0 的索引对应原文中第 k_star1 个样本年份与 M-K 输出的年份相差一年以内都算吻合。我个人的工作习惯是对于 M-K 疑难的、UF 与 UB 多次交叉的季节序列用滑动 t 检验做第二轮验证。滑动 t 检验的做法是设定窗口长度常见取 5 或 7 年从序列第三个窗口开始逐步平移计算窗口内前后两段均值的 t 统计量。t 值超过临界值的年份如果与 M-K 交点差在两年内我就把该交点视为可靠的突变年份并在报告中注明。滑动 t 检验适合前置判断如果滑动 t 检验全程无显著 t 值即使 M-K 给出了清晰交点也建议对该突变结论保持保守表述。反过来如果滑动 t 检验显著但 M-K 交点找不到往往说明趋势是渐变而非突变应该考虑修改分析框架。最后说一个血泪经验M-K 突变检测的结果尤其依赖序列端点。如果你要对近几年的变化趋势做判断注意序列末端数据往往处于累积窗口尚未成熟的阶段末端交点看可以但不能作为严格论文结论。我在处理某地区降水序列时M-K 检测显示 2021 年有一个明显上升突变加长序列到 2023 年后发现这个突变消失了——纯属末端数据不足导致的假象。从那以后我每次跑完分析都会强制做一遍尾端截断验证把序列截去最后 3 年再跑如果突变点位置发生超过 2 年的漂移那这个突变点我宁可不上报也不要冒险写进结论。希望这个验证习惯也能帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
