Mann-Kendall检验全解析从原理到11种变体及Python实现做水文趋势分析、气候变化评估、环境监测数据研究的朋友大概率都跟Mann-Kendall检验打过交道。这个非参数趋势检验方法几乎成了论文里的默认选项一句采用Mann-Kendall检验对序列进行了趋势分析就能交代完方法部分。可真到自己动手处理数据时问题就接踵而至序列有明显自相关时MK结果还可靠吗月径流这种强季节性数据怎么处理多个站点的趋势怎么整体判断想找突变点又该用什么这些问题仅靠最基础的MK检验是无法回答的——好在MK家族已经发展出了一系列变体针对不同数据特征各有解法。这篇文章我打算把MK检验这件事讲透从最底层的数学原理到11种常见变体各自解决什么问题、适用什么场景、Python代码怎么落地最后给出一套完整的数据分析流程。内容偏实战代码可直接复制运行适合刚接触趋势分析的研究生也适合已经用了很久MK但想系统梳理一遍的从业者。1. 为什么MK检验能成为趋势分析的默认选项1.1 非参数方法的核心优势先聊一个基本问题做趋势分析的方法很多线性回归是最直觉的——把时间当成自变量数据当成因变量拟合一条直线看斜率显著不显著就够了。那为什么还要用MK检验关键在于非参数三个字。线性回归依赖数据服从正态分布、残差独立同分布这些假设而水文气象数据恰恰经常不满足这些条件降雨量偏态分布、极端事件造成离群值、数据经过对数变换后分布依然奇怪。MK检验完全不关心数据服从什么分布它只看数据点之间的相对大小关系。就算你对整个序列做了对数变换、标准化甚至单调函数映射MK检验的结果都不会变——因为秩次关系被完整保留了。这一点在实际工作中的价值非常大。比如处理径流数据时丰水年和枯水年的数值可以差出几十倍如果直接做线性回归少数极端年份可能把整条拟合线带偏而MK检验对异常值天然稳健。再加上MK检验能比较自然地处理缺失值——两两比较时只要跳过缺失数据即可不需要插补。这对于动辄几十年的水文观测序列来说非常友好早期记录不全的情况太常见了。1.2 MK检验的数学内核MK检验的思路可以用一句话概括如果序列存在单调趋势那么先出现的观测值和后出现的观测值之间应当系统性地偏大或偏小。具体来说对于长度为n的序列x₁, x₂, ..., xₙ定义S统计量[ S \sum_{i1}^{n-1} \sum_{ji1}^{n} \text{sign}(x_j - x_i) ]其中sign是符号函数括号内为正取1为负取-1相等取0。这个公式的含义是把所有数据点两两配对看后出现的点相对先出现的点是变大还是变小最后求和。举一个具体例子。假设有5个数据3, 5, 2, 8, 6。S统计量的计算过程如下配对比较符号(3, 5)5 31(3, 2)2 3-1(3, 8)8 31(3, 6)6 31(5, 2)2 5-1(5, 8)8 51(5, 6)6 51(2, 8)8 21(2, 6)6 21(8, 6)6 8-1S 6。正值说明整体上后面的比前面的大暗示存在上升趋势负值则对应下降趋势。这里的核心问题变成S多大才算统计显著在原假设序列没有趋势的前提下S的期望值为0其方差有理论公式无并列值时 [ \text{Var}(S) \frac{n(n-1)(2n5)}{18} ]有并列值时需要在方差中减去tied组带来的修正项 [ \text{Var}(S) \frac{n(n-1)(2n5) - \sum_{p} t_p(t_p-1)(2t_p5)}{18} ]其中t_p是第p个并列组的观测个数。之后用正态近似将S标准化为Z统计量当S 0Z (S - 1) / √Var(S)当S 0Z 0当S 0Z (S 1) / √Var(S)Z服从标准正态分布查表或用误差函数就能得到p值。通常取显著性水平α0.05|Z| 1.96时判定趋势显著。这套逻辑有个隐含前提S的分布可以用正态分布近似。这个近似在样本量较大时一般要求n ≥ 8~10效果良好但样本量很小时需要使用精确分布表。这也是为什么MK检验不太适合用于特别短的序列。S统计量本身还可以进一步归一化为Kendall tau相关系数[ \tau \frac{S}{n(n-1)/2} ]tau的取值范围在[-1, 1]之间表示趋势的强度方便不同序列之间进行比较。后面讲Theil-Sen斜率时还会再用到S逻辑。2. 先手写一遍核心逻辑再上现成库2.1 用Python手写MK检验理解原理最好的方式是自己写一遍代码。下面这段实现涵盖了MK检验的核心流程——S统计量、方差修正、Z值以及p值。import numpy as np from scipy.stats import norm def mann_kendall(x): 手写Mann-Kendall检验 参数: x: 一维数组, 时间序列数据 返回: dict: S统计量, Kendall tau, Z值, p值 x np.asarray(x, dtypefloat) n len(x) # 计算S统计量, 同时记录并列组 s 0 ties {} for i in range(n - 1): for j in range(i 1, n): diff x[j] - x[i] if diff 0: s 1 elif diff 0: s - 1 # 统计当前元素的并列次数 count np.sum(x x[i]) if count 1: ties[i] count # 方差修正 unique, counts np.unique(x, return_countsTrue) tie_terms counts[counts 1] var_s n * (n - 1) * (2 * n 5) / 18 if tie_terms.size 0: var_s - np.sum(tie_terms * (tie_terms - 1) * (2 * tie_terms 5)) / 18 # 标准化Z值 if s 0: z (s - 1) / np.sqrt(var_s) elif s 0: z (s 1) / np.sqrt(var_s) else: z 0 # 双侧检验p值 p 2 * (1 - norm.cdf(abs(z))) # Kendall tau tau s / (n * (n - 1) / 2) return { S: s, tau: tau, z: z, p_value: p } # 测试 data np.array([3, 5, 2, 8, 6, 9, 7, 11, 10, 13]) result mann_kendall(data) print(result)这段代码的时间复杂度是O(n²)对动辄几百上千个数据点的时间序列足够用了。如果数据量大到几万以上可以考虑用基于归并排序的方法把复杂度降到O(n log n)但实际MK检验很少会遇到这种规模不必过度优化。2.2 pymannkendall一句话调用所有变体手写一遍是为了理解原理但实际分析时我推荐直接使用pymannkendall库。这个库把MK家族大部分变体都封装好了接口统一结果结构清晰。pip install pymannkendall基础用法非常简洁import pymannkendall as mk import numpy as np # 生成一组带上升趋势的示例数据 np.random.seed(42) trend np.linspace(0, 2, 50) noise np.random.normal(0, 1, 50) data trend noise # 原始MK检验 result mk.original_test(data) # result是一个NamedTuple主要字段 print(result.trend) # increasing / decreasing / no trend print(result.h) # True/False是否拒绝原假设 print(result.p) # p值 print(result.z) # Z统计量 print(result.tau) # Kendall tau print(result.s) # S统计量 print(result.slope) # Theil-Sen斜率这里有个小细节值得注意original_test返回的结果里包含了slope字段它直接给出了趋势的速率估计。也就是说仅靠这一个函数你就能获得趋势方向显著性变化速率三件套这也是我在实际项目中输出结论时的标准配置。pymannkendall库支持的方法列表比较完整后面各变体的代码示例我都会基于这个库来实现。个别库没有直接封装的方法比如区域MK和部分MK我会给出自己实现的代码。3. 11种变体逐一拆解各自解决什么问题MK检验这么多变体核心原因只有一个实际数据太不守规矩了。有的序列有季节性波动有的有明显的自相关有的需要同时考虑多个站点或变量有的不仅要判断趋势方向还想定位趋势从哪年开始——每一种不守规矩都催生了对应的变体。下面这11种是我在实际工作中真正用过、也认为值得记在工具箱里的。3.1 原始MKOriginal Mann-Kendall这是最基础的版本上一节已经完整介绍过。它适用于数据序列基本独立、无明显周期性、无显著自相关的场景比如年径流量、年降水量这类已经按年聚合、时间跨度不太长的序列。使用时有一个建议先画一下序列的ACF图自相关函数图确认没有显著滞后相关性再放心使用原始版本。这个习惯能帮你避开后面几个大坑。3.2 季节MKSeasonal Mann-Kendall河流的月径流、城市的月降水量、近海的海表温度——这些数据天然带有强烈的季节性周期。如果直接对整条序列做MK检验季节性的高低起伏会严重干扰S统计量。比如夏季径流普遍高于冬季这种系统性的季节差异会掩盖真正的长期趋势甚至造成虚假趋势。季节MK的思路是把数据按季节分组1月的数据只和1月的比2月的只和2月的比。每个季节组内计算一个S_k然后所有季节的S相加得到总S。方差部分除了各组方差之和还要加上组间协方差如果季节之间存在相关性的话。import pymannkendall as mk # 假设monthly_data是12个月的逐月数据按月份排列成二维数组 # 每行代表一年每列代表一个月 # 或者直接传入一维序列并指定月份从1到12循环 result mk.seasonal_test(monthly_data) print(result.trend, result.p, result.Tau)实际处理时我习惯先把数据整理成年份×季节的矩阵结构这样可以直观地看出每年哪个季节贡献了主要趋势。特别注意季节MK处理的是周期性季节数据如果你的数据本身就是年尺度的不需要用这个方法。3.3 趋势预白化MKTrend-Free Pre-whitening MK自相关问题可以说是MK检验在实践中最大的隐形杀手。当序列存在正自相关即t时刻的值与t-1时刻的值正相关S统计量的方差会被低估导致检验更容易虚假显著——明明没有趋势却判定为有趋势。预白化Pre-whitening的思路很直接先把序列中的趋势成分去掉用残差估计自回归系数再用这个系数对残差做白化处理最后把趋势加回来重新检验。具体流程用Theil-Sen斜率β估计序列的线性趋势计算去趋势序列 Y_t Y_t - βt对Y_t拟合一阶自回归模型AR(1)估计系数ρ₁如果ρ₁显著通常检验|ρ₁|是否大于某个阈值计算白化序列 Y_t Yt - ρ₁·Y{t-1}将趋势加回Y_t Y_t βt对Y做标准的MK检验result mk.trend_free_pre_whitening_test(data) print(result.trend, result.p, result.slope)这个方法在pymannkendall里就是一行调用。但要注意一个问题预白化会损失一点检验功效——也就是说如果趋势真的存在预白化之后可能更难检测出来。这是控制第一类错误虚假报警与牺牲第二类错误的权衡。3.4 方差修正MKVariance Correction MK预白化MK的思路是改造数据方差修正MK的思路则是修正检验统计量本身。以Hamed和Rao1998提出的方法为代表它不改变S统计量的计算方式而是通过有效样本量来放大方差[ \text{Var}^(S) \text{Var}(S) \times \frac{n}{n^} ]其中n是实际样本量n*是考虑了自相关结构后的有效样本量。有效样本量的估计基于序列秩次的样本自相关函数实际计算时会截断到显著滞后期数。# Hamed Rao修正 result mk.hamed_rao_modification_test(data) # Yue Wang修正另一种方差修正方式 result2 mk.yue_wang_modification_test(data)方差修正的好处是不会改变原始序列的信息量——预白化会在白化过程中抹掉一部分信号而方差修正保留了完整序列。如果序列的自相关结构比较复杂或者你不想因为预白化损失太多功效方差修正通常是更好的选择。3.5 自举MKBootstrap MK理论上MK检验的方差公式是在一定假设下推导的如果数据结构过于复杂不满足这些假设理论方差可能不准确。这时可以用重采样方法绕开理论分布。自举MK的做法是对原始序列进行有放回的块状重采样Block Bootstrap生成大量与原始序列长度相同的重采样序列对每个重采样序列计算S统计量从而得到S的经验分布最终用经验分布计算p值。result mk.bootstrap_test(data, nsim2000) print(result.trend, result.p, result.slope)我把这个方法归为不守规矩数据的兜底方案。当序列存在复杂的周期结构、非线性模式或者你用了很多修正方法依然不放心时自举MK能给出更稳健的推断。代价是计算量大但现代计算机跑2000次重采样通常也就是几秒钟的事。3.6 多变量MKMultivariate MK如果你同时监测了水温、溶解氧、pH值等多个水质指标想判断整个水体的环境状态是否存在显著变化趋势单独对每个指标做MK检验会面临多重比较问题——指标一显著、指标二不显著结论怎么下多变量MK把多个变量联合建模利用变量之间的协方差结构构造一个综合检验统计量。原假设是所有变量都没有趋势备择假设是至少有一个变量存在趋势。这样能给出一个全局性的判断。# 生成3列数据代表3个变量 data_multi np.column_stack([var1, var2, var3]) result mk.multivariate_test(data_multi) print(result.trend, result.p)实际使用中需要注意多变量MK要求各变量的数据长度一致且站点或变量之间需要有一定的物理关联。如果把几个完全无关的指标强行放一起做多变量检验结果的含义会很模糊。3.7 区域MKRegional MK区域MK解决的问题和多变量MK类似但不完全一样。多变量MK检验的是多个变量在同一地点的趋势区域MK检验的是多个站点在同一区域的趋势。比如你要评估整个长江上游流域的降水趋势把所有站点分别做MK然后数一数几个显著这种做法很常见但不够严谨——站点之间通常存在空间相关不应该当作独立样本看待。区域MK的思路是将各站点的S统计量求和同时用各站点之间的协方差修正总方差然后计算区域S的标准化Z值。def regional_mk(site_data_list): 简化版区域MK site_data_list: 列表每个元素是一个站点的时间序列 n_sites len(site_data_list) s_values [] for site_data in site_data_list: result mk.original_test(site_data) s_values.append(result.s) s_total np.sum(s_values) # 简化为假设站点独立实际应考虑空间协方差 var_set [] for site_data in site_data_list: n len(site_data) var_set.append(n * (n - 1) * (2 * n 5) / 18) var_total np.sum(var_set) if s_total 0: z (s_total - 1) / np.sqrt(var_total) elif s_total 0: z (s_total 1) / np.sqrt(var_total) else: z 0 from scipy.stats import norm p 2 * (1 - norm.cdf(abs(z))) return {regional_S: s_total, z: z, p: p}简化版假设站点间相互独立但实际区域站点往往存在空间相关性此时方差会被低估。更严谨的做法是用交叉站点数据的秩次计算协方差矩阵或者用空间自举方法。如果只是想快速给一个区域整体趋势的判断这个简化版也够用了。3.8 部分MKPartial MK有时候趋势会受到第三方变量的干扰。比如你想评估某条河流的水质改善趋势但同期降水量明显增加可能稀释了污染物浓度——你需要判断的是在控制了降水这个协变量之后污染物浓度是否仍然有显著下降趋势。部分MKPartial Mann-Kendall test就是在控制一个或多个协变量的情况下检验目标变量的趋势。它的原理类似于偏相关系数先对目标变量和协变量分别做秩变换再计算偏秩相关最后检验部分秩相关是否显著。# var是目标变量, covar是协变量 result mk.partial_test(var, covar) print(result.trend, result.p)这个变体在环境监测、流行病学这类经常需要调整混杂因素的研究中特别实用。但注意它只能控制一个协变量多个协变量时可能需要基于秩回归的更复杂方法。3.9 序贯MKSequential MK / Mann-Kendall-Sneyers前面所有变体都是在回答整段序列有没有趋势序贯MK回答的是一个不同的问题如果趋势存在它是什么时候开始的这个方法由Sneyers提出核心在于计算两条曲线从序列开头逐步向后计算的UF统计量序列和从序列末尾逐步向前计算的UB统计量序列。当UF和UB在显著性水平线通常±1.96之间出现交叉交叉点就指示了趋势开始突变的位置。def sequential_mk(x): 序贯MK统计量UF和UB计算 n len(x) UF np.zeros(n) UB np.zeros(n) # 正向计算UF s 0 for k in range(1, n): for j in range(k): if x[k] x[j]: s 1 var_s k * (k - 1) * (2 * k 5) / 18 if s 0: UF[k] (s - 1) / np.sqrt(var_s) elif s 0: UF[k] (s 1) / np.sqrt(var_s) else: UF[k] 0 # 反向计算UB x_rev x[::-1] s 0 for k in range(1, n): for j in range(k): if x_rev[k] x_rev[j]: s 1 var_s k * (k - 1) * (2 * k 5) / 18 if s 0: UB[n - 1 - k] -((s - 1) / np.sqrt(var_s)) elif s 0: UB[n - 1 - k] -((s 1) / np.sqrt(var_s)) else: UB[n - 1 - k] 0 return UF, UB # 绘制突变点分析图 import matplotlib.pyplot as plt UF, UB sequential_mk(data) years np.arange(len(data)) plt.figure(figsize(10, 5)) plt.plot(years, UF, labelUF) plt.plot(years, UB, labelUB) plt.axhline(1.96, colorgray, linestyle--, linewidth0.8) plt.axhline(-1.96, colorgray, linestyle--, linewidth0.8) plt.legend() plt.show()UF和UB交叉点的位置就是趋势突变的大致年份。这个方法在气候突变研究中用得非常多——比如某流域径流在某个年代出现系统性减少序贯MK能帮助定位突变发生的时间窗口。3.10 分段MKSegmented MK序贯MK定位到了突变点接下来自然的问题是突变前后各段序列的趋势如何此时需要把序列在突变点处切开分别对前后两段做MK检验。分段MK一般配合变点检测方法使用。先用Pettitt检验或Buishand U检验找出最可能的变点位置然后分段做趋势检验。这种方法能揭示一种容易被忽视的现象整段序列可能显示无显著趋势但分段后其实前半段显著上升、后半段显著下降——整体相互抵消了。from scipy.stats import mannwhitneyu def pettitt_test(x): Pettitt变点检测返回最可能的变点位置 n len(x) U np.zeros(n) for t in range(n): sign_sum 0 for i in range(t 1): for j in range(t 1, n): sign_sum np.sign(x[i] - x[j]) U[t] abs(sign_sum) k np.argmax(U) 1 # 变点位置从1开始 return k # 找到变点后分两段做MK change_point pettitt_test(data) seg1_result mk.original_test(data[:change_point]) seg2_result mk.original_test(data[change_point:]) print(f变点位置: {change_point}) print(f第一段趋势: {seg1_result.trend}, p{seg1_result.p}) print(f第二段趋势: {seg2_result.trend}, p{seg2_result.p})这个组合拳变点检测分段MK是长序列分析中价值最高的工作流。它能帮你回答趋势是否发生了转折比笼统地报告56年来径流呈显著下降趋势有意义得多。3.11 Theil-Sen斜率估计趋势大小的标配严格来说Theil-Sen斜率不是MK检验的变体而是MK检验的最佳搭档——MK告诉你趋势是否存在Theil-Sen告诉你趋势有多大。计算方法也很符合MK的非参数气质对序列中所有ij的配对计算斜率 (x_j - x_i) / (j - i)最后取所有斜率的中位数。中位数对异常值不敏感比最小二乘回归的斜率稳健得多。from scipy.stats import theilslopes slope, intercept, low_slope, high_slope theilslopes(data, np.arange(len(data))) print(f趋势速率: {slope:.3f} 单位/年) print(f95%置信区间: [{low_slope:.3f}, {high_slope:.3f}])实际报告结果时我看到的规范格式是径流序列呈显著下降趋势p 0.05平均下降速率为XX万立方米/年Theil-Sen斜率。趋势显著性和趋势速率一个都不能少——p值能说明规律性斜率能说明重要性两者搭配才能完整描述趋势特征。变体总览编号变体名称解决的核心问题典型应用1原始MK无趋势检验基线年尺度序列基础分析2季节MK季节性周期干扰月径流、月降水分析3TFPW-MK自相关导致虚假显著有AR(1)结构的连续序列4方差修正MK自相关但不想牺牲信号长序列趋势再验证5自举MK复杂结构难以理论推断不规律数据兜底方案6多变量MK多指标联合趋势判断水质多指标监测7区域MK多站点空间相关流域或区域评估8部分MK控制协变量影响环境措施效果评估9序贯MK定位突变点年份气候突变分析10分段MK变点前后分段趋势长序列转折分析11Theil-Sen斜率估计趋势大小变化速率定量报告4. 面对数据怎么选型一个决策参考学了这么多变体真正面对一份数据时怎么选我的经验是按两步走先看数据本身的结构再看你要回答的研究问题。4.1 按数据特征选型第一步永远是画图和做诊断。拿到一条时间序列我会先做三件事画时间序列图看有没有明显的季节性波动每年的高点和低点是否有规律地重复有没有异常跳变画ACF图看自相关在滞后1阶或几阶是否显著超出置信区间做Pettitt或Buishand检验确认是否存在明显变点根据诊断结果参考下面的表格数据特征推荐方法年尺度、无自相关原始MK月尺度、有明显季节周期季节MK连续序列、ACF显示AR(1)显著TFPW-MK或方差修正MK有多站点/多变量区域MK或多变量MK存在协变量干扰部分MK疑似存在突变点序贯MK分段MK4.2 按研究问题选型数据特征之外还要看你想回答什么问题。同样是径流数据如果只想写一句近50年径流趋势显著那原始MK就够了如果研究报告要求说明趋势的速率必须配上Theil-Sen斜率如果想为水资源管理政策提供依据最好做变点检测因为哪一年开始减少比整体减少更有管理价值如果是流域尺度评估单独对每个站做MK然后报有3个站显著、2个站不显著这很容易被审稿人质疑应该用区域MK给出一个空间整体的结论有一个容易被忽视的点分析之前先想清楚检验的零假设是什么。MK的原假设是序列不存在趋势不同变体只是在不同数据假设下检验同一个原假设。把变体当成某种数据情况下的MK来理解选型就会自然很多。4.3 一个务实的组合套路我自己在项目里常用的套路是这样的先用原始MK和Theil-Sen斜率拿到基础结果——趋势方向、p值、速率做ACF检查自相关如果显著用TFPW-MK或方差修正MK跑一遍作为对照如果顺便想看突变点做序贯MK画出UF/UB曲线有变点的话用Pettitt定位然后分段做MK多站点数据就直接用区域MK或者把所有站点的MK结果列成表但用这几种方法对单一变量做结果更清晰的区域表现最终报告结果时我会并排展示原始MK的z值和p值、修正后的z值和p值、Theil-Sen斜率、突变点位置如果有。这样做的好处是透明——审稿人或决策者能看到数据在不同假设下的表现而不是得到一个被优化过的单一结论。5. 完整实战56年水文序列的趋势分析讲完理论和方法下面走一个实际案例。我用一组模拟数据模拟某水文站1961-2016年共56年的年径流序列单位亿立方米数据带有明显的下降趋势但被噪声覆盖。整个分析流程我会尽量还原真实项目中会做的事情。5.1 数据准备与自相关检查import numpy as np import pandas as pd import matplotlib.pyplot as plt import pymannkendall as mk from scipy.stats import theilslopes from statsmodels.graphics.tsaplots import plot_acf np.random.seed(2024) years np.arange(1961, 2017) n len(years) # 模拟年径流基准值 下降趋势 噪声 base 350 decline -1.2 # 每年下降1.2亿立方米 noise_std 30 runoff base decline * np.arange(n) np.random.normal(0, noise_std, n) # 绘制时序图和自相关图 fig, axes plt.subplots(1, 2, figsize(14, 4)) axes[0].plot(years, runoff, linewidth1.5) axes[0].set_title(Annual Runoff Series) axes[0].set_xlabel(Year) axes[0].set_ylabel(Runoff (10^8 m³)) plot_acf(runoff, lags15, axaxes[1]) plt.tight_layout() plt.show()ACF图如果显示滞后1阶的自相关明显超出蓝色置信区间说明序列存在显著自相关直接用原始MK需要谨慎。这个模拟数据是我用线性趋势加入噪声生成的真实数据的结构通常更复杂——可能存在非线性趋势、周期波动、均值突变等。无论如何第一步都是把图画出来用眼睛看一遍总不会错。5.2 三路并跑原始MK、TFPW-MK、方差修正MK# 原始MK result_raw mk.original_test(runoff) print( 原始MK ) print(f趋势: {result_raw.trend}, p值: {result_raw.p:.4f}, Z: {result_raw.z:.4f}) # TFPW-MK result_tfpw mk.trend_free_pre_whitening_test(runoff) print( TFPW-MK ) print(f趋势: {result_tfpw.trend}, p值: {result_tfpw.p:.4f}, Z: {result_tfpw.z:.4f}) # Hamed Rao方差修正 result_hr mk.hamed_rao_modification_test(runoff) print( Hamed Rao修正MK ) print(f趋势: {result_hr.trend}, p值: {result_hr.p:.4f}, Z: {result_hr.z:.4f}) # 自举MK result_bt mk.bootstrap_test(runoff) print( Bootstrap MK ) print(f趋势: {result_bt.trend}, p值: {result_bt.p:.4f}, Z: {result_bt.z:.4f}) # Theil-Sen斜率 slope, intercept, low, high theilslopes(runoff, np.arange(n)) print(f\nTheil-Sen斜率: {slope:.3f} 亿立方米/年) print(f95%置信区间: [{low:.3f}, {high:.3f}])运行结果会给出各方法下趋势方向、p值和Z值的对比。如果数据没有明显的自相关几种方法的结果应该高度一致如果自相关较强你可能会看到原始MK的p值明显比修正方法小——这就是自相关造成的虚假显著在起作用。我在实际报告里通常把上述结果汇总成一个表方法Z值p值趋势判断原始MK-x.xx0.00x显著下降TFPW-MK-x.xx0.0xx显著下降方差修正MK-x.xx0.0xx显著下降Bootstrap MK-x.xx0.0xx显著下降5.3 突变点分析与分段检验趋势分析不能只看整体定位趋势变化的时间节点同样重要。下面结合Pettitt检验和序贯MK判断径流序列是否存在突变。# Pettitt变点检测 def pettitt_test(x): n len(x) U np.zeros(n) for t in range(n): sign_sum 0 for i in range(t 1): for j in range(t 1, n): sign_sum np.sign(x[i] - x[j]) U[t] abs(sign_sum) k np.argmax(U) 1 return k change_point pettitt_test(runoff) print(fPettitt检测变点位置: {years[change_point - 1]}年) # 分段MK seg1 mk.original_test(runoff[:change_point]) seg2 mk.original_test(runoff[change_point:]) print(f变点之前: 趋势{seg1.trend}, p值{seg1.p:.4f}) print(f变点之后: 趋势{seg2.trend}, p值{seg2.p:.4f}) # 序贯MK可视化 UF, UB sequential_mk(runoff) # sequential_mk函数见3.9节 plt.figure(figsize(10, 5)) plt.plot(years, UF, labelUF, linewidth1.5) plt.plot(years, UB, labelUB, linewidth1.5) plt.axhline(1.96, colorgray, linestyle--, linewidth0.8, label0.05 significance) plt.axhline(-1.96, colorgray, linestyle--, linewidth0.8) plt.axvline(xyears[change_point - 1], colorred, linestyle-., alpha0.7, labelChange point) plt.xlabel(Year) plt.ylabel(Statistic) plt.legend() plt.show()这种变点检测分段MK的组合能揭示很多整体检验看不出来的结构。比如整体检验可能显示显著下降但分段后你会发现真正的下降集中在某一个时段而不是均匀地贯穿整个记录期——这对归因分析和水资源管理决策的意义完全不同。分析任务完成后完整的结论应该包括几个层次的信息整段趋势方向与显著性MK部分、趋势速率Theil-Sen斜率部分、趋势是否均匀或存在突变Pettitt/序贯MK部分。我在项目报告里通常就是这样组织结果。6. 从实际项目里踩过的坑6.1 自相关不检查分析等于白做这是我在早期项目里付出过代价的经验。当时处理一组逐月水质数据直接用原始MK跑出来p值0.01非常显著当时很兴奋。后来一位前辈提醒我检查自相关才发现滞后1阶的自相关系数高达0.7经过TFPW-MK修正后p值变成了0.08完全不显著。差别如此之大就是因为自相关导致方差被严重低估。现在我的流程固定为先做ACF诊断再选方法宁可多一步也不要让结论建立在错误的方法上。尤其是在发布论文或者给决策部门出报告的场景下虚假显著的结论比不显著的危害大得多。6.2 样本量太小正态近似不靠谱MK检验的Z统计量依赖正态近似这个近似在小样本下会失真。我用模拟实验验证过当n10时原始MK的p值经常低估真实值。实际工作中如果遇到短序列比如某些站点只有七八年的监测数据我会直接改用精确检验或使用自举MK而不是硬套正态近似。6.3 显著性 ≠ 趋势大小p值只有0.001不代表趋势量级很大只代表趋势存在的统计证据很强。反过来p值0.06也不代表趋势不存在可能只是数据噪声太大、样本量不足。我在报告中习惯把p值和Theil-Sen斜率放在一起讨论。一个统计学上显著的下降趋势如果斜率为-0.01单位/年实际意义可能非常有限一个边缘显著的趋势如果斜率达到-2单位/年管理层面反而需要认真对待。统计显著性和实际显著性永远是两回事。6.4 多重比较别数几个站点显著有些论文的做法是把20个站点的MK结果列一个表然后用其中12个站点呈显著下降趋势来支撑区域结论。这种做法在统计上是有问题的——即使所有站点都没有真实趋势按α0.05的显著性水平平均也会有1个站点呈现虚假显著。要给出区域整体的趋势判断应该优先用区域MK或多变量MK这类联合检验方法而不是靠数显著性个数。6.5 数据清洗时别做多余的处理MK检验只关心秩次数据变换对数、标准化等不会改变检验结果但插补方法的选择会影响结果——尤其当缺失值比例较高时。我的建议是缺失值较少时直接跳过即可MK检验天然能处理缺失比例高时先分析缺失是否随机这个前提再决定用线性插值还是别的填充方法。千万不要为了好看而用某种平均值填充这会强行降低序列的方差影响检验的可靠性。我自己现在做长时间序列的趋势分析已经形成了一套固定流程先画时序图和ACF图把数据结构和自相关摸清楚再根据数据结构在原始MK、季节MK、TFPW-MK和方差修正MK里选合适的方法然后配上Theil-Sen斜率报告趋势速率最后用Pettitt和序贯MK检查是否存在突变点。这套流程跑完基本能把一条序列的趋势特征描述完整了。如果你刚开始用MK检验我的建议是先把原始MK的原理彻底吃透再逐步接触各种变体——毕竟所有变体都是在原始逻辑上打补丁底层的秩序比较思想是贯通的。 Mann-Kendall检验全解析：从原理、11种变体到Python实现
Mann-Kendall检验全解析从原理到11种变体及Python实现做水文趋势分析、气候变化评估、环境监测数据研究的朋友大概率都跟Mann-Kendall检验打过交道。这个非参数趋势检验方法几乎成了论文里的默认选项一句采用Mann-Kendall检验对序列进行了趋势分析就能交代完方法部分。可真到自己动手处理数据时问题就接踵而至序列有明显自相关时MK结果还可靠吗月径流这种强季节性数据怎么处理多个站点的趋势怎么整体判断想找突变点又该用什么这些问题仅靠最基础的MK检验是无法回答的——好在MK家族已经发展出了一系列变体针对不同数据特征各有解法。这篇文章我打算把MK检验这件事讲透从最底层的数学原理到11种常见变体各自解决什么问题、适用什么场景、Python代码怎么落地最后给出一套完整的数据分析流程。内容偏实战代码可直接复制运行适合刚接触趋势分析的研究生也适合已经用了很久MK但想系统梳理一遍的从业者。1. 为什么MK检验能成为趋势分析的默认选项1.1 非参数方法的核心优势先聊一个基本问题做趋势分析的方法很多线性回归是最直觉的——把时间当成自变量数据当成因变量拟合一条直线看斜率显著不显著就够了。那为什么还要用MK检验关键在于非参数三个字。线性回归依赖数据服从正态分布、残差独立同分布这些假设而水文气象数据恰恰经常不满足这些条件降雨量偏态分布、极端事件造成离群值、数据经过对数变换后分布依然奇怪。MK检验完全不关心数据服从什么分布它只看数据点之间的相对大小关系。就算你对整个序列做了对数变换、标准化甚至单调函数映射MK检验的结果都不会变——因为秩次关系被完整保留了。这一点在实际工作中的价值非常大。比如处理径流数据时丰水年和枯水年的数值可以差出几十倍如果直接做线性回归少数极端年份可能把整条拟合线带偏而MK检验对异常值天然稳健。再加上MK检验能比较自然地处理缺失值——两两比较时只要跳过缺失数据即可不需要插补。这对于动辄几十年的水文观测序列来说非常友好早期记录不全的情况太常见了。1.2 MK检验的数学内核MK检验的思路可以用一句话概括如果序列存在单调趋势那么先出现的观测值和后出现的观测值之间应当系统性地偏大或偏小。具体来说对于长度为n的序列x₁, x₂, ..., xₙ定义S统计量[ S \sum_{i1}^{n-1} \sum_{ji1}^{n} \text{sign}(x_j - x_i) ]其中sign是符号函数括号内为正取1为负取-1相等取0。这个公式的含义是把所有数据点两两配对看后出现的点相对先出现的点是变大还是变小最后求和。举一个具体例子。假设有5个数据3, 5, 2, 8, 6。S统计量的计算过程如下配对比较符号(3, 5)5 31(3, 2)2 3-1(3, 8)8 31(3, 6)6 31(5, 2)2 5-1(5, 8)8 51(5, 6)6 51(2, 8)8 21(2, 6)6 21(8, 6)6 8-1S 6。正值说明整体上后面的比前面的大暗示存在上升趋势负值则对应下降趋势。这里的核心问题变成S多大才算统计显著在原假设序列没有趋势的前提下S的期望值为0其方差有理论公式无并列值时 [ \text{Var}(S) \frac{n(n-1)(2n5)}{18} ]有并列值时需要在方差中减去tied组带来的修正项 [ \text{Var}(S) \frac{n(n-1)(2n5) - \sum_{p} t_p(t_p-1)(2t_p5)}{18} ]其中t_p是第p个并列组的观测个数。之后用正态近似将S标准化为Z统计量当S 0Z (S - 1) / √Var(S)当S 0Z 0当S 0Z (S 1) / √Var(S)Z服从标准正态分布查表或用误差函数就能得到p值。通常取显著性水平α0.05|Z| 1.96时判定趋势显著。这套逻辑有个隐含前提S的分布可以用正态分布近似。这个近似在样本量较大时一般要求n ≥ 8~10效果良好但样本量很小时需要使用精确分布表。这也是为什么MK检验不太适合用于特别短的序列。S统计量本身还可以进一步归一化为Kendall tau相关系数[ \tau \frac{S}{n(n-1)/2} ]tau的取值范围在[-1, 1]之间表示趋势的强度方便不同序列之间进行比较。后面讲Theil-Sen斜率时还会再用到S逻辑。2. 先手写一遍核心逻辑再上现成库2.1 用Python手写MK检验理解原理最好的方式是自己写一遍代码。下面这段实现涵盖了MK检验的核心流程——S统计量、方差修正、Z值以及p值。import numpy as np from scipy.stats import norm def mann_kendall(x): 手写Mann-Kendall检验 参数: x: 一维数组, 时间序列数据 返回: dict: S统计量, Kendall tau, Z值, p值 x np.asarray(x, dtypefloat) n len(x) # 计算S统计量, 同时记录并列组 s 0 ties {} for i in range(n - 1): for j in range(i 1, n): diff x[j] - x[i] if diff 0: s 1 elif diff 0: s - 1 # 统计当前元素的并列次数 count np.sum(x x[i]) if count 1: ties[i] count # 方差修正 unique, counts np.unique(x, return_countsTrue) tie_terms counts[counts 1] var_s n * (n - 1) * (2 * n 5) / 18 if tie_terms.size 0: var_s - np.sum(tie_terms * (tie_terms - 1) * (2 * tie_terms 5)) / 18 # 标准化Z值 if s 0: z (s - 1) / np.sqrt(var_s) elif s 0: z (s 1) / np.sqrt(var_s) else: z 0 # 双侧检验p值 p 2 * (1 - norm.cdf(abs(z))) # Kendall tau tau s / (n * (n - 1) / 2) return { S: s, tau: tau, z: z, p_value: p } # 测试 data np.array([3, 5, 2, 8, 6, 9, 7, 11, 10, 13]) result mann_kendall(data) print(result)这段代码的时间复杂度是O(n²)对动辄几百上千个数据点的时间序列足够用了。如果数据量大到几万以上可以考虑用基于归并排序的方法把复杂度降到O(n log n)但实际MK检验很少会遇到这种规模不必过度优化。2.2 pymannkendall一句话调用所有变体手写一遍是为了理解原理但实际分析时我推荐直接使用pymannkendall库。这个库把MK家族大部分变体都封装好了接口统一结果结构清晰。pip install pymannkendall基础用法非常简洁import pymannkendall as mk import numpy as np # 生成一组带上升趋势的示例数据 np.random.seed(42) trend np.linspace(0, 2, 50) noise np.random.normal(0, 1, 50) data trend noise # 原始MK检验 result mk.original_test(data) # result是一个NamedTuple主要字段 print(result.trend) # increasing / decreasing / no trend print(result.h) # True/False是否拒绝原假设 print(result.p) # p值 print(result.z) # Z统计量 print(result.tau) # Kendall tau print(result.s) # S统计量 print(result.slope) # Theil-Sen斜率这里有个小细节值得注意original_test返回的结果里包含了slope字段它直接给出了趋势的速率估计。也就是说仅靠这一个函数你就能获得趋势方向显著性变化速率三件套这也是我在实际项目中输出结论时的标准配置。pymannkendall库支持的方法列表比较完整后面各变体的代码示例我都会基于这个库来实现。个别库没有直接封装的方法比如区域MK和部分MK我会给出自己实现的代码。3. 11种变体逐一拆解各自解决什么问题MK检验这么多变体核心原因只有一个实际数据太不守规矩了。有的序列有季节性波动有的有明显的自相关有的需要同时考虑多个站点或变量有的不仅要判断趋势方向还想定位趋势从哪年开始——每一种不守规矩都催生了对应的变体。下面这11种是我在实际工作中真正用过、也认为值得记在工具箱里的。3.1 原始MKOriginal Mann-Kendall这是最基础的版本上一节已经完整介绍过。它适用于数据序列基本独立、无明显周期性、无显著自相关的场景比如年径流量、年降水量这类已经按年聚合、时间跨度不太长的序列。使用时有一个建议先画一下序列的ACF图自相关函数图确认没有显著滞后相关性再放心使用原始版本。这个习惯能帮你避开后面几个大坑。3.2 季节MKSeasonal Mann-Kendall河流的月径流、城市的月降水量、近海的海表温度——这些数据天然带有强烈的季节性周期。如果直接对整条序列做MK检验季节性的高低起伏会严重干扰S统计量。比如夏季径流普遍高于冬季这种系统性的季节差异会掩盖真正的长期趋势甚至造成虚假趋势。季节MK的思路是把数据按季节分组1月的数据只和1月的比2月的只和2月的比。每个季节组内计算一个S_k然后所有季节的S相加得到总S。方差部分除了各组方差之和还要加上组间协方差如果季节之间存在相关性的话。import pymannkendall as mk # 假设monthly_data是12个月的逐月数据按月份排列成二维数组 # 每行代表一年每列代表一个月 # 或者直接传入一维序列并指定月份从1到12循环 result mk.seasonal_test(monthly_data) print(result.trend, result.p, result.Tau)实际处理时我习惯先把数据整理成年份×季节的矩阵结构这样可以直观地看出每年哪个季节贡献了主要趋势。特别注意季节MK处理的是周期性季节数据如果你的数据本身就是年尺度的不需要用这个方法。3.3 趋势预白化MKTrend-Free Pre-whitening MK自相关问题可以说是MK检验在实践中最大的隐形杀手。当序列存在正自相关即t时刻的值与t-1时刻的值正相关S统计量的方差会被低估导致检验更容易虚假显著——明明没有趋势却判定为有趋势。预白化Pre-whitening的思路很直接先把序列中的趋势成分去掉用残差估计自回归系数再用这个系数对残差做白化处理最后把趋势加回来重新检验。具体流程用Theil-Sen斜率β估计序列的线性趋势计算去趋势序列 Y_t Y_t - βt对Y_t拟合一阶自回归模型AR(1)估计系数ρ₁如果ρ₁显著通常检验|ρ₁|是否大于某个阈值计算白化序列 Y_t Yt - ρ₁·Y{t-1}将趋势加回Y_t Y_t βt对Y做标准的MK检验result mk.trend_free_pre_whitening_test(data) print(result.trend, result.p, result.slope)这个方法在pymannkendall里就是一行调用。但要注意一个问题预白化会损失一点检验功效——也就是说如果趋势真的存在预白化之后可能更难检测出来。这是控制第一类错误虚假报警与牺牲第二类错误的权衡。3.4 方差修正MKVariance Correction MK预白化MK的思路是改造数据方差修正MK的思路则是修正检验统计量本身。以Hamed和Rao1998提出的方法为代表它不改变S统计量的计算方式而是通过有效样本量来放大方差[ \text{Var}^(S) \text{Var}(S) \times \frac{n}{n^} ]其中n是实际样本量n*是考虑了自相关结构后的有效样本量。有效样本量的估计基于序列秩次的样本自相关函数实际计算时会截断到显著滞后期数。# Hamed Rao修正 result mk.hamed_rao_modification_test(data) # Yue Wang修正另一种方差修正方式 result2 mk.yue_wang_modification_test(data)方差修正的好处是不会改变原始序列的信息量——预白化会在白化过程中抹掉一部分信号而方差修正保留了完整序列。如果序列的自相关结构比较复杂或者你不想因为预白化损失太多功效方差修正通常是更好的选择。3.5 自举MKBootstrap MK理论上MK检验的方差公式是在一定假设下推导的如果数据结构过于复杂不满足这些假设理论方差可能不准确。这时可以用重采样方法绕开理论分布。自举MK的做法是对原始序列进行有放回的块状重采样Block Bootstrap生成大量与原始序列长度相同的重采样序列对每个重采样序列计算S统计量从而得到S的经验分布最终用经验分布计算p值。result mk.bootstrap_test(data, nsim2000) print(result.trend, result.p, result.slope)我把这个方法归为不守规矩数据的兜底方案。当序列存在复杂的周期结构、非线性模式或者你用了很多修正方法依然不放心时自举MK能给出更稳健的推断。代价是计算量大但现代计算机跑2000次重采样通常也就是几秒钟的事。3.6 多变量MKMultivariate MK如果你同时监测了水温、溶解氧、pH值等多个水质指标想判断整个水体的环境状态是否存在显著变化趋势单独对每个指标做MK检验会面临多重比较问题——指标一显著、指标二不显著结论怎么下多变量MK把多个变量联合建模利用变量之间的协方差结构构造一个综合检验统计量。原假设是所有变量都没有趋势备择假设是至少有一个变量存在趋势。这样能给出一个全局性的判断。# 生成3列数据代表3个变量 data_multi np.column_stack([var1, var2, var3]) result mk.multivariate_test(data_multi) print(result.trend, result.p)实际使用中需要注意多变量MK要求各变量的数据长度一致且站点或变量之间需要有一定的物理关联。如果把几个完全无关的指标强行放一起做多变量检验结果的含义会很模糊。3.7 区域MKRegional MK区域MK解决的问题和多变量MK类似但不完全一样。多变量MK检验的是多个变量在同一地点的趋势区域MK检验的是多个站点在同一区域的趋势。比如你要评估整个长江上游流域的降水趋势把所有站点分别做MK然后数一数几个显著这种做法很常见但不够严谨——站点之间通常存在空间相关不应该当作独立样本看待。区域MK的思路是将各站点的S统计量求和同时用各站点之间的协方差修正总方差然后计算区域S的标准化Z值。def regional_mk(site_data_list): 简化版区域MK site_data_list: 列表每个元素是一个站点的时间序列 n_sites len(site_data_list) s_values [] for site_data in site_data_list: result mk.original_test(site_data) s_values.append(result.s) s_total np.sum(s_values) # 简化为假设站点独立实际应考虑空间协方差 var_set [] for site_data in site_data_list: n len(site_data) var_set.append(n * (n - 1) * (2 * n 5) / 18) var_total np.sum(var_set) if s_total 0: z (s_total - 1) / np.sqrt(var_total) elif s_total 0: z (s_total 1) / np.sqrt(var_total) else: z 0 from scipy.stats import norm p 2 * (1 - norm.cdf(abs(z))) return {regional_S: s_total, z: z, p: p}简化版假设站点间相互独立但实际区域站点往往存在空间相关性此时方差会被低估。更严谨的做法是用交叉站点数据的秩次计算协方差矩阵或者用空间自举方法。如果只是想快速给一个区域整体趋势的判断这个简化版也够用了。3.8 部分MKPartial MK有时候趋势会受到第三方变量的干扰。比如你想评估某条河流的水质改善趋势但同期降水量明显增加可能稀释了污染物浓度——你需要判断的是在控制了降水这个协变量之后污染物浓度是否仍然有显著下降趋势。部分MKPartial Mann-Kendall test就是在控制一个或多个协变量的情况下检验目标变量的趋势。它的原理类似于偏相关系数先对目标变量和协变量分别做秩变换再计算偏秩相关最后检验部分秩相关是否显著。# var是目标变量, covar是协变量 result mk.partial_test(var, covar) print(result.trend, result.p)这个变体在环境监测、流行病学这类经常需要调整混杂因素的研究中特别实用。但注意它只能控制一个协变量多个协变量时可能需要基于秩回归的更复杂方法。3.9 序贯MKSequential MK / Mann-Kendall-Sneyers前面所有变体都是在回答整段序列有没有趋势序贯MK回答的是一个不同的问题如果趋势存在它是什么时候开始的这个方法由Sneyers提出核心在于计算两条曲线从序列开头逐步向后计算的UF统计量序列和从序列末尾逐步向前计算的UB统计量序列。当UF和UB在显著性水平线通常±1.96之间出现交叉交叉点就指示了趋势开始突变的位置。def sequential_mk(x): 序贯MK统计量UF和UB计算 n len(x) UF np.zeros(n) UB np.zeros(n) # 正向计算UF s 0 for k in range(1, n): for j in range(k): if x[k] x[j]: s 1 var_s k * (k - 1) * (2 * k 5) / 18 if s 0: UF[k] (s - 1) / np.sqrt(var_s) elif s 0: UF[k] (s 1) / np.sqrt(var_s) else: UF[k] 0 # 反向计算UB x_rev x[::-1] s 0 for k in range(1, n): for j in range(k): if x_rev[k] x_rev[j]: s 1 var_s k * (k - 1) * (2 * k 5) / 18 if s 0: UB[n - 1 - k] -((s - 1) / np.sqrt(var_s)) elif s 0: UB[n - 1 - k] -((s 1) / np.sqrt(var_s)) else: UB[n - 1 - k] 0 return UF, UB # 绘制突变点分析图 import matplotlib.pyplot as plt UF, UB sequential_mk(data) years np.arange(len(data)) plt.figure(figsize(10, 5)) plt.plot(years, UF, labelUF) plt.plot(years, UB, labelUB) plt.axhline(1.96, colorgray, linestyle--, linewidth0.8) plt.axhline(-1.96, colorgray, linestyle--, linewidth0.8) plt.legend() plt.show()UF和UB交叉点的位置就是趋势突变的大致年份。这个方法在气候突变研究中用得非常多——比如某流域径流在某个年代出现系统性减少序贯MK能帮助定位突变发生的时间窗口。3.10 分段MKSegmented MK序贯MK定位到了突变点接下来自然的问题是突变前后各段序列的趋势如何此时需要把序列在突变点处切开分别对前后两段做MK检验。分段MK一般配合变点检测方法使用。先用Pettitt检验或Buishand U检验找出最可能的变点位置然后分段做趋势检验。这种方法能揭示一种容易被忽视的现象整段序列可能显示无显著趋势但分段后其实前半段显著上升、后半段显著下降——整体相互抵消了。from scipy.stats import mannwhitneyu def pettitt_test(x): Pettitt变点检测返回最可能的变点位置 n len(x) U np.zeros(n) for t in range(n): sign_sum 0 for i in range(t 1): for j in range(t 1, n): sign_sum np.sign(x[i] - x[j]) U[t] abs(sign_sum) k np.argmax(U) 1 # 变点位置从1开始 return k # 找到变点后分两段做MK change_point pettitt_test(data) seg1_result mk.original_test(data[:change_point]) seg2_result mk.original_test(data[change_point:]) print(f变点位置: {change_point}) print(f第一段趋势: {seg1_result.trend}, p{seg1_result.p}) print(f第二段趋势: {seg2_result.trend}, p{seg2_result.p})这个组合拳变点检测分段MK是长序列分析中价值最高的工作流。它能帮你回答趋势是否发生了转折比笼统地报告56年来径流呈显著下降趋势有意义得多。3.11 Theil-Sen斜率估计趋势大小的标配严格来说Theil-Sen斜率不是MK检验的变体而是MK检验的最佳搭档——MK告诉你趋势是否存在Theil-Sen告诉你趋势有多大。计算方法也很符合MK的非参数气质对序列中所有ij的配对计算斜率 (x_j - x_i) / (j - i)最后取所有斜率的中位数。中位数对异常值不敏感比最小二乘回归的斜率稳健得多。from scipy.stats import theilslopes slope, intercept, low_slope, high_slope theilslopes(data, np.arange(len(data))) print(f趋势速率: {slope:.3f} 单位/年) print(f95%置信区间: [{low_slope:.3f}, {high_slope:.3f}])实际报告结果时我看到的规范格式是径流序列呈显著下降趋势p 0.05平均下降速率为XX万立方米/年Theil-Sen斜率。趋势显著性和趋势速率一个都不能少——p值能说明规律性斜率能说明重要性两者搭配才能完整描述趋势特征。变体总览编号变体名称解决的核心问题典型应用1原始MK无趋势检验基线年尺度序列基础分析2季节MK季节性周期干扰月径流、月降水分析3TFPW-MK自相关导致虚假显著有AR(1)结构的连续序列4方差修正MK自相关但不想牺牲信号长序列趋势再验证5自举MK复杂结构难以理论推断不规律数据兜底方案6多变量MK多指标联合趋势判断水质多指标监测7区域MK多站点空间相关流域或区域评估8部分MK控制协变量影响环境措施效果评估9序贯MK定位突变点年份气候突变分析10分段MK变点前后分段趋势长序列转折分析11Theil-Sen斜率估计趋势大小变化速率定量报告4. 面对数据怎么选型一个决策参考学了这么多变体真正面对一份数据时怎么选我的经验是按两步走先看数据本身的结构再看你要回答的研究问题。4.1 按数据特征选型第一步永远是画图和做诊断。拿到一条时间序列我会先做三件事画时间序列图看有没有明显的季节性波动每年的高点和低点是否有规律地重复有没有异常跳变画ACF图看自相关在滞后1阶或几阶是否显著超出置信区间做Pettitt或Buishand检验确认是否存在明显变点根据诊断结果参考下面的表格数据特征推荐方法年尺度、无自相关原始MK月尺度、有明显季节周期季节MK连续序列、ACF显示AR(1)显著TFPW-MK或方差修正MK有多站点/多变量区域MK或多变量MK存在协变量干扰部分MK疑似存在突变点序贯MK分段MK4.2 按研究问题选型数据特征之外还要看你想回答什么问题。同样是径流数据如果只想写一句近50年径流趋势显著那原始MK就够了如果研究报告要求说明趋势的速率必须配上Theil-Sen斜率如果想为水资源管理政策提供依据最好做变点检测因为哪一年开始减少比整体减少更有管理价值如果是流域尺度评估单独对每个站做MK然后报有3个站显著、2个站不显著这很容易被审稿人质疑应该用区域MK给出一个空间整体的结论有一个容易被忽视的点分析之前先想清楚检验的零假设是什么。MK的原假设是序列不存在趋势不同变体只是在不同数据假设下检验同一个原假设。把变体当成某种数据情况下的MK来理解选型就会自然很多。4.3 一个务实的组合套路我自己在项目里常用的套路是这样的先用原始MK和Theil-Sen斜率拿到基础结果——趋势方向、p值、速率做ACF检查自相关如果显著用TFPW-MK或方差修正MK跑一遍作为对照如果顺便想看突变点做序贯MK画出UF/UB曲线有变点的话用Pettitt定位然后分段做MK多站点数据就直接用区域MK或者把所有站点的MK结果列成表但用这几种方法对单一变量做结果更清晰的区域表现最终报告结果时我会并排展示原始MK的z值和p值、修正后的z值和p值、Theil-Sen斜率、突变点位置如果有。这样做的好处是透明——审稿人或决策者能看到数据在不同假设下的表现而不是得到一个被优化过的单一结论。5. 完整实战56年水文序列的趋势分析讲完理论和方法下面走一个实际案例。我用一组模拟数据模拟某水文站1961-2016年共56年的年径流序列单位亿立方米数据带有明显的下降趋势但被噪声覆盖。整个分析流程我会尽量还原真实项目中会做的事情。5.1 数据准备与自相关检查import numpy as np import pandas as pd import matplotlib.pyplot as plt import pymannkendall as mk from scipy.stats import theilslopes from statsmodels.graphics.tsaplots import plot_acf np.random.seed(2024) years np.arange(1961, 2017) n len(years) # 模拟年径流基准值 下降趋势 噪声 base 350 decline -1.2 # 每年下降1.2亿立方米 noise_std 30 runoff base decline * np.arange(n) np.random.normal(0, noise_std, n) # 绘制时序图和自相关图 fig, axes plt.subplots(1, 2, figsize(14, 4)) axes[0].plot(years, runoff, linewidth1.5) axes[0].set_title(Annual Runoff Series) axes[0].set_xlabel(Year) axes[0].set_ylabel(Runoff (10^8 m³)) plot_acf(runoff, lags15, axaxes[1]) plt.tight_layout() plt.show()ACF图如果显示滞后1阶的自相关明显超出蓝色置信区间说明序列存在显著自相关直接用原始MK需要谨慎。这个模拟数据是我用线性趋势加入噪声生成的真实数据的结构通常更复杂——可能存在非线性趋势、周期波动、均值突变等。无论如何第一步都是把图画出来用眼睛看一遍总不会错。5.2 三路并跑原始MK、TFPW-MK、方差修正MK# 原始MK result_raw mk.original_test(runoff) print( 原始MK ) print(f趋势: {result_raw.trend}, p值: {result_raw.p:.4f}, Z: {result_raw.z:.4f}) # TFPW-MK result_tfpw mk.trend_free_pre_whitening_test(runoff) print( TFPW-MK ) print(f趋势: {result_tfpw.trend}, p值: {result_tfpw.p:.4f}, Z: {result_tfpw.z:.4f}) # Hamed Rao方差修正 result_hr mk.hamed_rao_modification_test(runoff) print( Hamed Rao修正MK ) print(f趋势: {result_hr.trend}, p值: {result_hr.p:.4f}, Z: {result_hr.z:.4f}) # 自举MK result_bt mk.bootstrap_test(runoff) print( Bootstrap MK ) print(f趋势: {result_bt.trend}, p值: {result_bt.p:.4f}, Z: {result_bt.z:.4f}) # Theil-Sen斜率 slope, intercept, low, high theilslopes(runoff, np.arange(n)) print(f\nTheil-Sen斜率: {slope:.3f} 亿立方米/年) print(f95%置信区间: [{low:.3f}, {high:.3f}])运行结果会给出各方法下趋势方向、p值和Z值的对比。如果数据没有明显的自相关几种方法的结果应该高度一致如果自相关较强你可能会看到原始MK的p值明显比修正方法小——这就是自相关造成的虚假显著在起作用。我在实际报告里通常把上述结果汇总成一个表方法Z值p值趋势判断原始MK-x.xx0.00x显著下降TFPW-MK-x.xx0.0xx显著下降方差修正MK-x.xx0.0xx显著下降Bootstrap MK-x.xx0.0xx显著下降5.3 突变点分析与分段检验趋势分析不能只看整体定位趋势变化的时间节点同样重要。下面结合Pettitt检验和序贯MK判断径流序列是否存在突变。# Pettitt变点检测 def pettitt_test(x): n len(x) U np.zeros(n) for t in range(n): sign_sum 0 for i in range(t 1): for j in range(t 1, n): sign_sum np.sign(x[i] - x[j]) U[t] abs(sign_sum) k np.argmax(U) 1 return k change_point pettitt_test(runoff) print(fPettitt检测变点位置: {years[change_point - 1]}年) # 分段MK seg1 mk.original_test(runoff[:change_point]) seg2 mk.original_test(runoff[change_point:]) print(f变点之前: 趋势{seg1.trend}, p值{seg1.p:.4f}) print(f变点之后: 趋势{seg2.trend}, p值{seg2.p:.4f}) # 序贯MK可视化 UF, UB sequential_mk(runoff) # sequential_mk函数见3.9节 plt.figure(figsize(10, 5)) plt.plot(years, UF, labelUF, linewidth1.5) plt.plot(years, UB, labelUB, linewidth1.5) plt.axhline(1.96, colorgray, linestyle--, linewidth0.8, label0.05 significance) plt.axhline(-1.96, colorgray, linestyle--, linewidth0.8) plt.axvline(xyears[change_point - 1], colorred, linestyle-., alpha0.7, labelChange point) plt.xlabel(Year) plt.ylabel(Statistic) plt.legend() plt.show()这种变点检测分段MK的组合能揭示很多整体检验看不出来的结构。比如整体检验可能显示显著下降但分段后你会发现真正的下降集中在某一个时段而不是均匀地贯穿整个记录期——这对归因分析和水资源管理决策的意义完全不同。分析任务完成后完整的结论应该包括几个层次的信息整段趋势方向与显著性MK部分、趋势速率Theil-Sen斜率部分、趋势是否均匀或存在突变Pettitt/序贯MK部分。我在项目报告里通常就是这样组织结果。6. 从实际项目里踩过的坑6.1 自相关不检查分析等于白做这是我在早期项目里付出过代价的经验。当时处理一组逐月水质数据直接用原始MK跑出来p值0.01非常显著当时很兴奋。后来一位前辈提醒我检查自相关才发现滞后1阶的自相关系数高达0.7经过TFPW-MK修正后p值变成了0.08完全不显著。差别如此之大就是因为自相关导致方差被严重低估。现在我的流程固定为先做ACF诊断再选方法宁可多一步也不要让结论建立在错误的方法上。尤其是在发布论文或者给决策部门出报告的场景下虚假显著的结论比不显著的危害大得多。6.2 样本量太小正态近似不靠谱MK检验的Z统计量依赖正态近似这个近似在小样本下会失真。我用模拟实验验证过当n10时原始MK的p值经常低估真实值。实际工作中如果遇到短序列比如某些站点只有七八年的监测数据我会直接改用精确检验或使用自举MK而不是硬套正态近似。6.3 显著性 ≠ 趋势大小p值只有0.001不代表趋势量级很大只代表趋势存在的统计证据很强。反过来p值0.06也不代表趋势不存在可能只是数据噪声太大、样本量不足。我在报告中习惯把p值和Theil-Sen斜率放在一起讨论。一个统计学上显著的下降趋势如果斜率为-0.01单位/年实际意义可能非常有限一个边缘显著的趋势如果斜率达到-2单位/年管理层面反而需要认真对待。统计显著性和实际显著性永远是两回事。6.4 多重比较别数几个站点显著有些论文的做法是把20个站点的MK结果列一个表然后用其中12个站点呈显著下降趋势来支撑区域结论。这种做法在统计上是有问题的——即使所有站点都没有真实趋势按α0.05的显著性水平平均也会有1个站点呈现虚假显著。要给出区域整体的趋势判断应该优先用区域MK或多变量MK这类联合检验方法而不是靠数显著性个数。6.5 数据清洗时别做多余的处理MK检验只关心秩次数据变换对数、标准化等不会改变检验结果但插补方法的选择会影响结果——尤其当缺失值比例较高时。我的建议是缺失值较少时直接跳过即可MK检验天然能处理缺失比例高时先分析缺失是否随机这个前提再决定用线性插值还是别的填充方法。千万不要为了好看而用某种平均值填充这会强行降低序列的方差影响检验的可靠性。我自己现在做长时间序列的趋势分析已经形成了一套固定流程先画时序图和ACF图把数据结构和自相关摸清楚再根据数据结构在原始MK、季节MK、TFPW-MK和方差修正MK里选合适的方法然后配上Theil-Sen斜率报告趋势速率最后用Pettitt和序贯MK检查是否存在突变点。这套流程跑完基本能把一条序列的趋势特征描述完整了。如果你刚开始用MK检验我的建议是先把原始MK的原理彻底吃透再逐步接触各种变体——毕竟所有变体都是在原始逻辑上打补丁底层的秩序比较思想是贯通的。