简介面向从事气候变化、水文气象等领域的数据分析人员这份Matlab源码包实现了经典的Mann-Kendall(MK)非参数检验方法可用于判断降水、气温、干旱频次等气候序列中是否存在显著趋势并借助UF、UB统计量曲线识别突变发生的时间。资源共包含3个m文件涵盖主执行脚本与核心计算函数代码结构简洁便于读者根据自身数据调整输入格式与检验参数适合初学者对照学习MK检验原理也为科研工作者提供可直接调用的工具脚本。整个压缩包仅3KB轻量便携已有3405人学习下载。通过运行这些源码可以直观理解MK检验中UF和UB曲线的含义掌握突变点判定方法并快速迁移到自己的研究数据中。1. MK检验在MATLAB里要先看懂UF和UB两条曲线的含义拿到一段逐年径流数据或NDVI时序先问趋势是否显著再问趋势从哪一年开始变。Mann-Kendall检验MK检验的Z值只能回答前者要回答后者得看MATLAB画出来的UF和UB两条曲线。UF是正序序列逐点累积得到的MK统计量UB是逆序序列回算后取负得到的对应曲线两条线的交点落在置信区间内时一般把这个位置视为突变点候选。很多人照着教程跑通了图却说不清每条线每个点是什么调参数时只能靠试。下面把定义、代码和解读一步步捋清楚适合做水文、气候、遥感趋势分析想自己控制计算过程的工程师。2. UF与UB在MK检验里的定义和计算顺序2.1 MK统计量S和标准化Z值UF曲线的每一帧都在算什么MK检验的核心不复杂就是对序列里所有数据对 (x_i, x_j)ij比较大小S Σ sign(x_j - x_i)sign(z) 在 z0 时取 1z0 时取 -1z0 时取 0。S 的符号代表整体趋势方向S 的绝对值代表趋势强度。在无趋势的零假设下S 的期望为 0方差为 Var(S) n(n-1)(2n5)/18其中 n 是序列长度。当 n≥8 时S 近似服从正态分布于是可以构造标准化统计量S0 时Z (S-1) / √Var(S)S0 时Z (S1) / √Var(S)S0 时Z 0加减 1 是连续性修正目的是让离散的 S 分布更贴近正态分布。实际计算时别漏掉这一步否则临界值对比会偏松。UF曲线就是把这个过程按“逐帧”来做对前 k 个样本单独算一次 S 和 Z得到 UF_kk 从 2 取到 n。所以 UF 曲线上的第 k 个点表示“只看从第 1 年到第 k 年的数据趋势方向和显著性是什么样”。它是一根不断把新数据点并入计算窗口的累积趋势曲线。下面是向量化示意x randn(100, 1); % 模拟无趋势序列 n length(x); UF zeros(n, 1); for k 2:n S sum(sign(x(k) - x(1:k-1))); % 第k个点与前k-1个点逐对比较 VarS k*(k-1)*(2*k5)/18; % 用当前窗口长度k而不是总长n if S 0 UF(k) (S - 1) / sqrt(VarS); elseif S 0 UF(k) (S 1) / sqrt(VarS); end end这段代码里最容易被忽略的是 VarS 用 k 计算而不是用总长 n。很多网上流传的旧脚本在这里写死 n序列较长时 UF 曲线会被人为压扁突变点位置也会偏移。用当前窗口长度 k 才是标准做法。2.2 逆序序列与UB为什么要倒过来再算一遍UF 能看出趋势从开始到当前时刻的累积变化但看不出“从当前时刻往结尾看”是什么状态。为了找到趋势方向发生变化的位置需要另一条从尾部回看的曲线。做法是把原始序列倒过来得到逆序列 y_i x_{n-i1}对 y 重复逐帧计算得到 UF再倒序映射回原始时间轴并取负UB_k -UF_{n-k1}取负的原因很简单逆序列里 y_i y_j 表示原始序列中靠后的值比靠前的小也就是原序列在下降。如果不取负UB 的符号会和 UF 完全相反两条曲线在图上很难形成有意义的交叉。取负之后UB 和 UF 处于同一个“上升为正、下降为负”的坐标系里。对逆序列做计算时MATLAB 里用 flipud 或 fliplr 即可。逆序列的起点对应原始序列的终点所以 UB 曲线在右侧末端对应逆序列的前几个点画图时要把统计量序列再倒回来才能和 UF 对齐到同一时间轴上。这一步的先后顺序错了交叉点位置就对不上。需要提醒的是单调上升序列也会让 UF 和 UB 在中段产生交点。UF 从 0 向上走UB 从正值方向往 0 回落两条线必然在中间相交。因此交叉点这个信号本身只说明“趋势特征从某个时刻起前后不一致”不能立刻当成突变结论。还要看交点的显著性、原始序列的形态和实际物理背景。2.3 显著性水平和置信区间怎么读MK检验的 UF 和 UB 曲线通常和一条置信区间带一起画。常用的显著性水平与临界值如下alpha临界值 ±z适用场景0.05±1.9695%置信区间最常用0.10±1.6490%置信区间短序列或探索性分析0.01±2.5899%置信区间长序列且要求高可靠判断规则分两步。第一步看 UF 是否超出上下临界线UF 超过 1.96 说明存在显著上升趋势低于 -1.96 说明存在显著下降趋势。第二步看 UF 和 UB 的交点位置交点落在两条临界线之间才说明突变发生在统计上可辨识的时段里交点在区间外时突变点附近的数据变异性太大不能下明确结论。画图时不需要额外工具箱直接画两条水平参考线即可plot(t, 1.96*ones(size(t)), k--)再画一条对称的 -1.96。临界值按 alpha 参数换用即可alpha0.1 时换成 1.64alpha0.01 时换成 2.58。3. 在MATLAB中实现MK检验UF与UB曲线一份可直接调用的函数3.1 主函数mk_uf_ub完整代码把逐帧计算封装成函数便于重复调用。下面是一个不依赖额外工具箱核心功能的版本function [UF, UB, t] mk_uf_ub(x, alpha) % MK检验UF与UB曲线计算 % 输入: % x - 一维时间序列列向量或行向量 % alpha - 显著性水平默认0.05 % 输出: % UF - 正序序列逐段标准化MK统计量 % UB - 逆序序列回算取负后的统计量 % t - 与UF/UB对应的序号从1开始 if nargin 2, alpha 0.05; end x x(:); n length(x); if n 8 error(序列长度至少为8否则MK统计量不满足正态近似); end UF zeros(n, 1); UB zeros(n, 1); for k 2:n S sum(sign(x(k) - x(1:k-1))); % 正序逐对比较 VarS k * (k - 1) * (2*k 5) / 18; if S 0 UF(k) (S - 1) / sqrt(VarS); elseif S 0 UF(k) (S 1) / sqrt(VarS); end end xr flipud(x); % 逆序序列 UR zeros(n, 1); for k 2:n S sum(sign(xr(k) - xr(1:k-1))); VarS k * (k - 1) * (2*k 5) / 18; if S 0 UR(k) (S - 1) / sqrt(VarS); elseif S 0 UR(k) (S 1) / sqrt(VarS); end end UB -flipud(UR); % 倒序映射并取负 t (1:n); if nargout 0 z norminv(1 - alpha/2); % 需统计工具箱可换成查表值1.96/1.64/2.58 plot(t, UF, b-, t, UB, r-, ... t, z*ones(n,1), k--, t, -z*ones(n,1), k--); legend(UF, UB, 上界, 下界, Location, best); grid on; xlabel(时间序号); ylabel(MK统计量); end end这个函数把正序和逆序两段分开写逻辑清晰。如果不想依赖统计工具箱可以把 nargout0 分支里的 norminv 改成查表值alpha0.05 用 1.96alpha0.1 用 1.64alpha0.01 用 2.58把对应关系直接写在注释里。3.2 函数参数表和调用方式下面这张表列出了输入输出参数及其含义方便复制到自己的工程里对照使用。参数类型说明xdouble vector待检验的时间序列函数内部会转为列向量alphadouble scalar显著性水平决定置信区间临界值默认0.05UFdouble vector正序MK统计量与x等长UF(1)0UBdouble vector逆序MK统计量与x等长UB(1)0tdouble vector时间序号对应绘图横轴调用方式x randn(80, 1) linspace(0, 2, 80); [UF, UB, t] mk_uf_ub(x);如果只想要图直接调用mk_uf_ub(x)就会走 nargout0 分支自动绘图。需要注意该分支依赖统计工具箱的 norminv在没有工具箱的 MATLAB 环境里建议在函数外面单独画图或者把 norminv 换成固定临界值。提示函数内部 nargout0 分支只做示意实际工程里我更推荐把绘图放在函数外这样可以把 UF/UB 的计算结果保存下来方便后续做阈值标注或多图对比。3.3 边界条件NaN、平局值和序列长度实际数据很少是理想时间序列喂给函数之前要处理三个问题。NaN 处理直接用 sign(x(k)-x(1:k-1)) 时NaN 会进入 sum统计量直接变成 NaN。建议先插值fillmissing(x, linear)可以处理中间缺测点缺测集中在两端时直接截断更稳妥。平局值当 x(k)x(j) 时 sign0不贡献 S 统计量。MK检验对少量平局是稳健的但平局占比超过 30% 时要谨慎解读因为连续相等的值会压缩方差项的可靠性。序列长度n8 时方差公式误差大函数直接报错。对于长度 8-20 的短序列建议用 alpha0.1 的临界值否则很容易漏检。短序列的 Z 值离散度高紧贴 ±1.96 的结果本质上不可靠。4. 用带突变点的仿真数据验证UF和UB的交叉位置4.1 构造一组“趋势跳跃噪声”的测试数据要验证算法先得知道自己埋进去的突变点在哪。构造长度 100 的时间序列前 50 个点线性上升第 51 个点整体抬高之后继续线性上升再加上白噪声。rng(2024); % 固定随机种子保证可复现 t (1:100); x 0.05*t 2.0*randn(100,1); x(51:100) x(51:100) 3; % 第51个点起整体抬升构造突变0.05*t 是基础趋势随机噪声标准差为 2突变幅度为 3。相对于噪声这个跳变很容易被 MK检验捕获。如果把噪声标准差改成 5突变幅度 3 就会被淹没UF 曲线的上穿也不稳定这正好说明 MK检验检测不了信噪比过低的突变。4.2 绘制UF和UB曲线并解读关键交叉调用函数并画图[UF, UB, tt] mk_uf_ub(x); figure(Color, w); plot(tt, UF, b-, LineWidth, 1.2); hold on; plot(tt, UB, r-, LineWidth, 1.2); plot(tt, 1.96*ones(size(tt)), k--); plot(tt, -1.96*ones(size(tt)), k--); plot([51 51], [-4 4], g:, LineWidth, 1); legend(UF, UB, 95%上限, 95%下限, 埋入突变点, Location, northwest); xlabel(时间序号); ylabel(MK统计量); grid on; hold off;从图上能看到三个关键特征。UF 曲线从第 50 多点开始向上抬升随后突破上置信限说明前段上升趋势显著UB 曲线从右侧回看在中段与 UF 形成交叉交叉位置基本落在第 51 号附近交叉之后两条线反向张开说明趋势特征发生明显改变。绿色虚线标出的是实际埋入突变点与 UF/UB 交叉位置吻合。plot([51 51], [-4 4])是画一条竖直参考线新版MATLAB可以直接用xline(51)代替。如果交叉点比真实突变点滞后 3-5 个点属于正常现象因为 MK统计量是累积量需要窗口包含足够多的突变后样本后才会产生明显转向。4.3 三个容易误判的情形第一是单调序列的假交叉。纯上升序列的 UF 和 UB 在中间位置也会相交但这种交叉对应的是趋势过程的几何对称不是突变。判断方法是看交叉点两侧是否出现趋势方向的实质性改变或者直接看原始序列的一阶差分是否有符号翻转。第二是交点落在置信区间外。两条线相交了但交点在 -1.96 以下或 1.96 以上说明背景变异性太大这个时刻的位置不可信。此时不应下突变结论而是把原始序列按不同长度分段做 MK检验对比各段趋势是否一致。第三是出现多个交点。当序列含周期成分或多次跳跃时UF 与 UB 会在多个位置交叉。先看年内季节性比如 NDVI 数据做 MK检验前应先去季节循环如果仍有多个交点优先选择连续 2-3 个点都靠近同一位置的交叉群孤立的单一交叉点一般不作为突变依据。5. 实战进阶自相关修正、Sen斜率估计与季节MK检验5.1 用Sen斜率估计趋势幅度MK只给方向和显著性不给趋势幅度。实际写报告时经常同时输出 Sen 斜率Theil-Sen估计它是所有点对斜率的全组合中位数function beta sen_slope(x) n length(x); slopes []; for i 1:n-1 for j i1:n slopes [slopes; (x(j) - x(i)) / (j - i)]; end end beta median(slopes); endSen 斜率的单位是“每时间步的变量变化量”。时间轴不均匀时j-i 要替换成实际时间差。它的优点是抗离群点一个极端值最多影响一条点对斜率不会像最小二乘那样直接拉动整条回归线。数据量在几千个点以内时双重循环可接受更大数据量用 nchoosek 生成组合下标但注意内存占用。5.2 序列存在自相关时先做TFPW预白化再跑MK水位、气温这类序列普遍存在滞后一阶自相关会让 MK检验的有效样本量虚高本不显著的趋势被判定为显著。常见做法是趋势无关预白化TFPW先用 Sen 斜率去掉趋势对残差做 AR(1) 白化再把趋势加回去b sen_slope(x); n length(x); detrended x - b * (1:n); % 去掉线性趋势 rho corr(detrended(1:end-1), detrended(2:end)); if abs(rho) 0.05 x_new detrended; % 自相关太小不需要白化 else res detrended(2:end) - rho * detrended(1:end-1); x_new res b * (2:n); % 把趋势加回白化后的残差 end白化后序列长度少 1加回趋势时时间轴从 2 开始否则和原始序列长不对齐。阈值 0.05 是经验值也可以改成显著性检验来决定。TFPW 后的序列再传给 mk_uf_ub输出的 UF/UB 会更保守交叉点也更稳定。5.3 季节数据用季节MK检验避免伪趋势含完整季节周期的数据直接跑 MK检验会把季节差异当成趋势生成大量假交叉。标准做法是季节 MK检验每个季节或月份单独计算 S 统计量和方差再把各季节的 S 加总方差按季节独立累加% s为季节序号S_season为各季节内部MK统计量 % Var_season为各季节对应方差 S_total sum(S_season); Var_total sum(Var_season); if S_total 0 Z (S_total - 1) / sqrt(Var_total); elseif S_total 0 Z (S_total 1) / sqrt(Var_total); end如果不同季节间存在显著的跨季节相关Hirsch-Slack 方法会在总方差里加入顺序成对样本的协方差项。这个修正项让 Var_total 变大Z 值更保守。做完季节 MK检验后再回到第 4 章的绘图流程UF 和 UB 交叉点才会有实际分析意义seasonal Mann-Kendall 的方差修正把连续月份间的协方差项考虑进去后原来的临界值判断才可靠。本文还有配套的精品资源点击获取
