简介这份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 的时间序列逐对比较所有观测值的大小关系。如果序列整体呈上升趋势那么靠后的值普遍大于靠前的值正向比较的个数会显著多于负向比较反之亦然。这种比较只关心谁大谁小不关心具体数值差距所以它天然免疫异常值的影响。数学表达上定义统计量 SS Σ(i1 to n-1) Σ(ji1 to n) sign(x_j - x_i)其中 sign() 是符号函数。当 n 较大时S 近似服从正态分布均值为 0方差由下式给出Var(S) [n(n-1)(2n5) - Σ t_p(t_p-1)(2t_p5)] / 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 i1: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_p5); end end VarS (n*(n-1)*(2*n5) - 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 和 UBn 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-k1))); UB(k) -Z; % 注意符号取反 endUB 计算时的符号取反经常被忽略。原因是逆序序列的趋势方向与原始序列相反标准化后的 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]; %#okAGROW 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)]; %#okAGROW 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% 置信水平alpha0.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.96UB 在中段与之交叉输出突变年份 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-k1))); 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最终值})]; %#okAGROW 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 k1: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)); endPettitt 检验用双重循环计算时间复杂度 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-k1))); 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 会退化为普通循环不报错但不会提速。本文还有配套的精品资源点击获取