MATLAB实现Mann-Kendall突变检验：原理、代码与实战

简介：这份PDF教程面向环境科学、气候学、地质学等领域需要分析时间序列突变点的研究人员与学生，以MATLAB为工具，系统讲解Mann-Kendall非参数检验的实现方法。手册从方法起源、统计原理讲起，逐步拆解UFk、UBk统计量的计算流程，完整覆盖秩序列构造、均值方差计算、显著性交叉点判断等关键步骤；随附的MATLAB脚本包含数据读取、秩次计算、绘图和显著性检验等功能，还给出了批量处理Excel数据并输出统计量表的扩展代码，便于读者对照验证和修改复用。资源为单个PDF文件，大小仅409KB，轻量便携，适合在科研和课程设计中随时查阅。该资源已有305人在线学习，具备一定参考价值。

1. 从一条水文序列里读出"哪一年真的变了"

做 Mann-Kendall 突变检验的人，手里通常不止一条曲线，而是一条十几年甚至几十年的观测序列：年径流量、降水距平、NDVI、气温……你早就从折线图里看出趋势有转折，但写论文或者出报告时不能只说"看起来变了"，需要一个不依赖分布假设、能指出具体突变年份的统计量。Mann-Kendall 检验正好补这个位置，它本质上是基于秩的非参数方法，对正态性和线性关系都不挑剔，所以水文、气象、环境领域几乎把它当默认的突变检测工具。

MATLAB 做这件事，常见路径有三条：自己写循环统计量、用 Statistics and Machine Learning Toolbox 的kruskalwallis凑合、或者直接调 File Exchange 上的现成函数。我的建议是自己写核心部分，因为 M-K 检验本身没有复杂矩阵运算，几十行代码就能跑完，而自己写的好处是可以完全控制置信区间、突变点定位逻辑和绘图样式。这篇博文按"理论 → 代码 → 参数调优 → 陷阱规避"的顺序来讲，最后给出一个可以直接改数据就跑的完整脚本框架。

2. M-K 检验的统计量构造与突变点定位原理

2.1 为什么基于秩构造统计量能绕开分布假设

Mann-Kendall 检验的核心思想非常简单：对长度为 n 的时间序列，逐对比较所有观测值的大小关系。如果序列整体呈上升趋势，那么靠后的值普遍大于靠前的值，正向比较的个数会显著多于负向比较；反之亦然。这种比较只关心"谁大谁小"，不关心具体数值差距，所以它天然免疫异常值的影响。

数学表达上，定义统计量 S：

S = Σ(i=1 to n-1) Σ(j=i+1 to n) sign(x_j - x_i)

其中 sign() 是符号函数。当 n 较大时，S 近似服从正态分布，均值为 0，方差由下式给出：

Var(S) = [n(n-1)(2n+5) - Σ t_p(t_p-1)(2t_p+5)] / 18

这个公式里第二项是处理并列值（ties）的修正项，t_p 是第 p 组并列值的个数。没有并列值时这一项为 0。标准化的检验统计量 Z 就很好算了：

Z = (S - 1) / sqrt(Var(S)) (S > 0) Z = 0 (S = 0) Z = (S + 1) / sqrt(Var(S)) (S < 0)

Z 的正负告诉你是上升还是下降趋势，绝对值超过 1.96（对应 95% 置信水平）就说明趋势显著。到这里为止，算的还只是趋势显著性，突变检验在这个基础上多做一步：对序列中每个时间点 k，把序列分成前后两段，分别计算两段的 M-K 统计量，然后看两个统计量的差异是否在某一点达到最大。这个"差异最大点"就是候选突变点。

2.2 顺序统计量 UF 与逆序统计量 UB 的几何含义

实际用得最多的突变定位方法是"顺序/逆序曲线交叉法"，也叫 UF-UB 曲线法。它的做法是：对原始序列 x_1, x_2, ..., x_n，从第一个数据点开始，依次计算到第 k 个点为止的子序列的 M-K 统计量，标准化后记为 UF_k。这个值刻画的是"截至 k 时刻，序列相对起始时刻的变化趋势是否显著"。

然后对原始序列做逆序处理：从最后一个点倒着往前，同样逐点计算 M-K 统计量，标准化后记为 UB_k。几何上，UF 曲线从前往后看趋势累积，UB 曲线从后往前看趋势累积。如果序列在某个时刻发生了真正的突变，UF 和 UB 会在该时刻附近交叉，且交叉点落在置信区间边界（通常取 ±1.96）之内。交叉点对应的横坐标就是突变发生年份。

这个方法的直观解释是：突变前，从前往后看的趋势和从后往前看的趋势在突变点处取值必然交错；如果序列完全平稳，两条曲线会在零附近小幅波动，不会形成明显交叉。需要说明的是，UF 和 UB 的曲线交叉只是突变存在的必要条件，不是充分条件——两条线可能在置信带外交叉，那种情况通常说明序列变化是渐变的而非突变的。

2.3 突变点与趋势检验的边界：什么时候该信交叉点

用 UF-UB 曲线判断突变点，有一个容易忽略的前提：序列里只能有一个主导性的突变信号。如果序列存在两个以上显著突变点，比如 1990 年一次跃升、2005 年一次跃降，UF 和 UB 的交叉可能出现在两个位置，或者只显示其中一个更明显的信号，另一个被淹没。这种情况下需要配合滑动 t 检验或 Pettitt 检验做交叉验证。

另一个边界条件是序列长度。n < 8 时，S 的分布离散性太强，正态近似完全不成立，UF-UB 曲线会剧烈抖动，交叉点不可信。n 在 8~15 之间时，建议用连续校正公式（上面给出的那个），同时考虑用蒙特卡洛模拟校准 p 值。n > 30 时正态近似效果很好，直接用就行。

3. MATLAB 手工实现 M-K 突变检验的完整代码

3.1 数据准备：从 Excel 读入并预处理时间序列

写代码之前先说数据格式。M-K 检验要求输入是一维等间距时间序列，缺失值要提前处理。实测中常见的问题是 Excel 里日期列格式不统一，读进来后序列错位。我一般用readtable读入，明确指定列名，然后用ismissing检查缺失值。数据处理遵循"宁缺毋滥"的原则：缺失比例低于 5% 用线性插值补，高于 5% 就直接删除对应年份，因为插值会人为改变秩次关系，影响检验结果。

% 读入数据：假设 Excel 有两列，Year 和 Value data = readtable('hydrological_data.xlsx', 'VariableNamingRule', 'preserve'); year = data.Year; value = data.Value; % 缺失值处理：少于5%缺失用线性插值，否则剔除整行 missing_idx = ismissing(value); if mean(missing_idx) < 0.05 value(missing_idx) = interp1(year(~missing_idx), value(~missing_idx), year(missing_idx), 'linear'); fprintf('缺失值已通过线性插值填补，共 %d 个\n', sum(missing_idx)); else valid_idx = ~missing_idx; year = year(valid_idx); value = value(valid_idx); warning('缺失值超过5%%，已剔除对应年份'); end

这段代码里VariableNamingRule设为'preserve'是为了防止 MATLAB 把中文或特殊字符列名自动转换，实测中很关键。插值用interp1的线性方法是最保守的选择，不要用spline或pchip，它们会引入不存在的波动模式。

数据长度检查逻辑也很重要：numel(value)小于 10 就提示用户数据太短，因为 M-K 检验在短序列上的功效很低，就算算出结果也不适合下结论。

3.2 核心计算：UF 序列、UB 序列与置信带

计算 UF/UB 序列的核心是写一个子函数，输入一段子序列，输出标准化后的 M-K 统计量。按照前面公式，需要处理并列值修正、连续校正、以及 S 和 Var(S) 的计算。以下是完整实现：

function [Z, S, VarS] = mk_stat(x) % 输入：x 为一段序列向量 % 输出：Z 标准化统计量, S 原始统计量, VarS 方差 n = length(x); S = 0; for i = 1:n-1 for j = i+1:n S = S + sign(x(j) - x(i)); end end % 计算并列值修正项 unique_vals = unique(x); ties = 0; for k = 1:length(unique_vals) t_p = sum(x == unique_vals(k)); if t_p > 1 ties = ties + t_p*(t_p-1)*(2*t_p+5); end end VarS = (n*(n-1)*(2*n+5) - ties) / 18; % 连续校正并标准化 if S > 0 Z = (S - 1) / sqrt(VarS); elseif S < 0 Z = (S + 1) / sqrt(VarS); else Z = 0; end end

这个函数里有个细节值得注意：计算并列值修正时，用了unique找全部不重复值，然后逐个统计出现次数。如果数据是连续的浮点数，并列值几乎不会出现，但水文数据里常出现"径流量恰好相等"或"降水天数取整后相同"的情况，这个修正项就不能省。此外，sign函数是 MATLAB 内置的，返回 -1、0、1，不需要自己写判断逻辑。

接着用这个子函数计算 UF 和 UB：

n = length(value); % 计算 UF 序列（正向） UF = zeros(n, 1); for k = 2:n [Z, ~, ~] = mk_stat(value(1:k)); UF(k) = Z; end % 计算 UB 序列（逆向） UB = zeros(n, 1); for k = 2:n % 对逆序序列取反，统一用 mk_stat 计算 [Z, ~, ~] = mk_stat(value(end:-1:(end-k+1))); UB(k) = -Z; % 注意符号取反 end

UB 计算时的符号取反经常被忽略。原因是逆序序列的趋势方向与原始序列相反，标准化后的 Z 值需要取反才能和 UF 在同一坐标轴下对比。比如原始序列从第 20 年开始上升，倒着看就是从第 20 年开始下降，反映在逆序统计量上是负值，取反后变成正值，才能正向和 UF 的上升趋势对应。

3.3 突变年份提取与 95% 置信区间判定

算出 UF 和 UB 后，突变点定位逻辑是：找两条曲线在置信带内的交叉点。严谨的 MATLAB 实现如下：

alpha = 0.05; Z_alpha = norminv(1 - alpha/2); % 1.96 % 找交叉点索引 cross_idx = []; for k = 2:n if sign(UF(k-1) - UB(k-1)) ~= sign(UF(k) - UB(k)) cross_idx = [cross_idx, k]; %#ok<AGROW> end end % 筛选在置信带内的交叉点 mutate_year = []; for k = cross_idx % 交叉点前后各取一个点，判断是否在置信带内 if abs(UF(k)) < Z_alpha || abs(UB(k)) < Z_alpha mutate_year = [mutate_year, year(k)]; %#ok<AGROW> end end if isempty(mutate_year) fprintf('未检测到置信水平 %.0f%% 下的显著突变点\n', (1-alpha)*100); else fprintf('检测到突变年份: '); fprintf('%d ', mutate_year); fprintf('\n'); end

这里有一个判定细节：交叉点对应的 UF 或 UB 值必须在置信带内（绝对值小于 1.96）。实际做的时候经常遇到两条曲线在置信带外反复交叉的情况，这通常说明序列波动剧烈，突变信号不干净，不建议当作突变年份输出。此外，如果交叉点出现在序列前 10% 或后 10% 的位置，可信度也要打折扣，因为端点附近的 UF/UB 计算基于极短子序列，方差大、统计量波动剧烈。

3.4 绘制 UF-UB 曲线图：论文级排版与标注

绘图时我用下面的模板，包含置信区间上下界虚线、UF/UB 曲线、突变年份参考线：

figure('Position', [100, 100, 800, 450]); hold on; grid on; % 置信带 fill([year(1), year(end), year(end), year(1)], ... [Z_alpha, Z_alpha, -Z_alpha, -Z_alpha], ... [0.9, 0.9, 0.9], 'EdgeColor', 'none', 'FaceAlpha', 0.3); plot(year, UF, 'b-', 'LineWidth', 1.5, 'DisplayName', 'UF'); plot(year, UB, 'r-', 'LineWidth', 1.5, 'DisplayName', 'UB'); % 置信区间线 yline(Z_alpha, 'k--', 'LineWidth', 0.8); yline(-Z_alpha, 'k--', 'LineWidth', 0.8); % 标记突变年份 if ~isempty(mutate_year) xline(mutate_year(1), 'g--', 'LineWidth', 1.2, ... 'Label', sprintf('%d', mutate_year(1)), 'LabelOrientation', 'top'); end xlabel('年份'); ylabel('M-K 统计量'); legend('Location', 'northwest'); set(gca, 'FontSize', 11, 'LineWidth', 0.6); hold off;

绘图时三个参数值得调：FaceAlpha控制置信带透明度，设 0.3 比较合适，太大遮住曲线，太小看不出区间；LabelOrientation设为'top'避免年份标签压住坐标轴；Position里的宽高比按 800x450 设置，放到 Word 或 LaTeX 里不需要二次裁剪。

4. 参数敏感性与多场景实战验证

4.1 置信水平选择：95% vs 90% 对突变点判定差异

默认用 95% 置信水平（alpha=0.05）是惯例，但实际数据里经常出现"UF 和 UB 在置信带边缘交叉"的情况——交叉点对应的统计量绝对值在 1.8 到 2.0 之间，这时候 95% 水平下不显著，90% 水平下（1.645）就显著。我在多个数据集上跑过，降水类数据尤其常见这种边界情况。

处理原则是：如果 95% 下检测不到，但 UF 的累积趋势线确实存在明显转向，就把置信水平放宽到 90% 再跑一次，并在报告里注明"在 90% 置信水平下显著"。不建议直接修改数据或挑对自己有利的置信水平，而是把两个水平的结果都列出来，交给读者判断。

% 同时计算两个置信水平的临界值 Z_90 = norminv(1 - 0.10/2); % 1.645 Z_95 = norminv(1 - 0.05/2); % 1.960 fprintf('90%% 置信区间边界: %.3f\n', Z_90); fprintf('95%% 置信区间边界: %.3f\n', Z_95);

norminv是标准正态分布的逆累积分布函数，参数里1 - alpha/2对应双侧检验。注意不要写成norminv(alpha/2)，那样取到的是负边界。

4.2 趋势成分干扰：先做预白化还是一阶差分

M-K 检验有个被论文里反复讨论的弱点：当序列存在显著的长期趋势（比如持续升温或持续城市化导致的径流递减）时，UF/UB 曲线会整体向一个方向漂移，导致交叉点定位偏移。我遇到过一个典型案例：某流域年径流量因为上游取水持续下降，UF 曲线一路向下穿越 -1.96，UB 在中段与之交叉，输出"突变年份 2008"，但实际上 2008 年没有发生任何突变性事件，纯粹是趋势的附属产物。

处理这种问题的常见做法是趋势预白化（Trend-Free Pre-Whitening）。步骤是：先用线性回归拟合整个序列的斜率，从原始序列中减去拟合趋势，对残差做 M-K 检验。MATLAB 实现如下：

% 线性去趋势 t = (1:n)'; p = polyfit(t, value, 1); trend_line = polyval(p, t); detrended = value - trend_line; % 对去趋势后的序列做 M-K 检验 [Z_detrended, S_detrended, VarS_detrended] = mk_stat(detrended); fprintf('去趋势后 Z = %.3f, p = %.4f\n', Z_detrended, 2*(1-normcdf(abs(Z_detrended))));

去趋势后如果 Z 不再显著，说明原序列的变化主要是线性趋势而非突变；如果依然显著且 UF/UB 出现新的交叉点，才适合下突变结论。这个步骤在降水、径流、气温数据上都适用，代价是可能掩盖真实的突变信号——如果突变本身是阶梯式的，线性去趋势会削弱它的强度。所以去趋势后的结果只作参考，不替代原序列的检验。

4.3 多站点数据批处理：一个脚本跑完所有结果

实际项目里很少只分析一个站点。水文站网动辄几十个站点，每个站点都要做突变检验并输出结果表。写一个循环批处理脚本，配合站点元数据表，可以大幅减少手工操作。代码结构如下：

% 假设站点数据存储在 sites_data.mat 中，结构为数组 load('sites_data.mat'); % 包含 site_names, years, values_matrix results = table(); for i = 1:length(site_names) value_i = values_matrix(:, i); valid = ~ismissing(value_i); if sum(valid) < 10 warning('站点 %s 有效数据不足，跳过', site_names{i}); continue; end % 截取有效区间 year_i = years(valid); value_i = value_i(valid); % 计算 UF/UB n_i = length(value_i); UF_i = zeros(n_i, 1); UB_i = zeros(n_i, 1); for k = 2:n_i [Z_uf, ~, ~] = mk_stat(value_i(1:k)); [Z_ub, ~, ~] = mk_stat(value_i(end:-1:(end-k+1))); UF_i(k) = Z_uf; UB_i(k) = -Z_ub; end % 找交叉点 diff_sign = sign(UF_i(2:end) - UB_i(2:end)) .* sign(UF_i(1:end-1) - UB_i(1:end-1)); cross_pos = find(diff_sign < 0) + 1; % 筛选置信带内突变点 valid_cross = cross_pos(abs(UF_i(cross_pos)) < Z_95 | abs(UB_i(cross_pos)) < Z_95); if ~isempty(valid_cross) mutate_yr = year_i(valid_cross(1)); else mutate_yr = NaN; end % 存储结果 results = [results; table(site_names{i}, mutate_yr, Z_uf, 'VariableNames', ... {'站点', '突变年份', 'UF最终值'})]; %#ok<AGROW> end % 导出结果表 writetable(results, 'mk_results.xlsx');

这段代码里用了diff_sign做交叉检测，比for循环加if判断快很多。需要注意find里的 +1 偏移，因为diff后的向量比原向量短一个位置。

5. 处理并列值、季节性与短序列的进阶策略

5.1 并列值修正后被忽略的边界情况

前文代码里包含了并列值修正公式，但实际使用中，当数据存在大量离散取值（如降水天数、径流等级）时，并列值的占比会很高。模拟数据测试表明，当并列值组数超过总长度的一半时，修正后的 Var(S) 可能趋近于零，标准化后的 Z 值会异常偏大或偏小。

遇到这种情况有两个出路：一是改用 Exact Test（精确检验）而非正态近似，MATLAB 里没有直接内置函数，但可以用置换检验模拟 p 值；二是检查数据是否有太多不合理的相等值——有时是数据记录精度问题，比如把径流量四舍五入到整数导致大量相等值，这种情况应该退回原始精度数据。

5.2 季节性数据：先按月拆分再合并结果

M-K 检验原始形式假设数据是独立同分布的，季节性数据（如月降水量）直接跑会导致强自相关，虚假的交叉点大幅增多。水文气象领域处理这类数据的标准做法是先按月拆分，对每个月份分别做 M-K 检验，再通过 Fisher 方法合并各月 p 值，或者只在特定季节（汛期、雨季）上单独分析。

% 按月拆分后分别计算 Z 值 years_unique = unique(year); nt = length(years_unique); monthly_Z = NaN(12, 1); for m = 1:12 idx = (month == m); if sum(idx) < 10 continue; end monthly_data = value(idx); [Z_m, ~, ~] = mk_stat(monthly_data); monthly_Z(m) = Z_m; end % 可视化各月趋势强度 bar(1:12, monthly_Z); xlabel('月份'); ylabel('M-K Z 值'); yline(1.96, 'r--'); yline(-1.96, 'r--');

这个"按月独立检验"的思路，比把所有月份混在一起跑更符合水文规律。混在一起跑的问题在于：汛期和非汛期的数据掺在一起，趋势方向可能互相抵消，或者因为季节差异产生额外的秩次混乱，导致 UF/UB 曲线震荡剧烈。

5.3 与 Pettitt 检验的交叉验证：避免单一方法误判

UF-UB 曲线法的局限在于它本质是"趋势累积过程的可视化诊断"，突变点定位有一定主观性（交叉不一定是干净的单点交叉）。我通常再用 Pettitt 检验做交叉验证，两者结果一致的突变年份可信度较高。

Pettitt 检验的核心思想是把序列分成两段，找到使两段经验分布差异最大的分割点。用 MATLAB 写主循环代码如下：

function [k, U_max, p_value] = pettitt_test(x) n = length(x); U = zeros(1, n-1); for k = 1:n-1 U(k) = 0; for i = 1:k for j = k+1:n U(k) = U(k) + sign(x(i) - x(j)); end end end [U_max, k_idx] = max(abs(U)); K = k_idx; % 近似 p 值（基于正态近似） p_value = 2 * exp(-6 * U_max^2 / (n^3 + n^2)); end

Pettitt 检验用双重循环计算，时间复杂度 O(n²)，对 50 年以内的序列完全够用。超过 100 年建议用向量化或 C 代码加速。如果 UF-UB 和 Pettitt 给出的突变年份相差超过 3 年，基本可以断定序列中有多个突变信号或趋势成分干扰，需要结合专业判断取舍。

5.4 一个容易踩的坑：unique排序语义

计算并列值修正时用unique(x)得到的是排序后的唯一值列表，但 MATLAB 的unique对浮点数默认按数值大小排序，对 NaN 的处理是排在最后。如果数据里有 NaN 没剔除干净，unique会把 NaN 当作一个值参与计数，sum(x == NaN)恒为 0，不会报错但会静默丢失修正项。所以代码里第一步必须是先过滤 NaN 再进 mk_stat，而不是在函数内部处理。

另外，数据量超过 10000 时，双重循环的 mk_stat 会成为性能瓶颈。用meshgrid向量化替代循环是常见的优化方向，但在典型的水文序列（30-100 年）里，这个优化收益不明显，代码可读性优先。

6. 用蒙特卡洛模拟验证突变检验的可靠性

如果你对自己的数据特别不放心，或者审稿人质疑突变年份的显著性，可以用蒙特卡洛模拟做一个"零假设检验"：对原始序列做随机重排，每次重排后跑一遍 UF-UB 交叉检测，统计随机情况下出现交叉点的概率。如果原始序列的突变强度超过 95% 的随机重排，说明突变点不是随机波动造成的。

rng(42); n_sim = 1000; cross_count = zeros(n_sim, 1); for s = 1:n_sim % 随机重排序列（无放回） perm_idx = randperm(n); perm_value = value(perm_idx); % 重排序列的 UF-UB UF_p = zeros(n, 1); UB_p = zeros(n, 1); for k = 2:n [Z_uf_p, ~, ~] = mk_stat(perm_value(1:k)); [Z_ub_p, ~, ~] = mk_stat(perm_value(end:-1:(end-k+1))); UF_p(k) = Z_uf_p; UB_p(k) = -Z_ub_p; end % 检测交叉点 diff_p = sign(UF_p(2:end) - UB_p(2:end)) .* ... sign(UF_p(1:end-1) - UB_p(1:end-1)); cross_p = find(diff_p < 0) + 1; % 是否有交叉点在置信带内 in_band = any(abs(UF_p(cross_p)) < Z_95 | abs(UB_p(cross_p)) < Z_95); cross_count(s) = in_band; end p_sim = mean(cross_count); fprintf('蒙特卡洛模拟 p 值: %.3f\n', p_sim);

模拟结果如果 p > 0.05，说明在纯随机序列里也经常出现类似的交叉点，此时应当谨慎解释突变检验结果。这个结论虽然略显"扫兴"，但在学术写作中反而是加分项——审稿人看到你做过负向验证，会认为你的方法使用是严谨的。

蒙特卡洛模拟的计算量主要在内层双重循环，1000 次重排对 50 年数据大概需要一分钟左右，可以接受。如果数据量翻倍，建议用parfor并行重排来提速：

parfor s = 1:n_sim % 循环体同上，cross_count(s) 赋值 end

使用parfor前先跑parpool检查并行池是否可用，否则 MATLAB 会退化为普通循环，不报错但不会提速。

本文还有配套的精品资源，点击获取