简介MK检验是气候诊断与预测中常用的非参数统计方法不依赖样本分布假设适用于判断气候序列是否存在突变及其发生时间也可用于降水、干旱频次等要素的趋势显著性检验。这套Matlab源码包面向气候、水文、环境等专业研究者与学生包含三个脚本文件分别实现主流程、Z统计量计算与改进型M-K检验代码结构清晰注释完整运行后可直接输出UF和UB曲线帮助快速定位突变时刻并支持调整置信水平。压缩包仅3KB轻量易用适合二次开发或嵌入现有研究流程。资源已有3405人学习无论是处理长序列数据还是撰写科研论文都可作为直观的算法参考与可复用模板尤其适合希望弄懂UF与UB含义的Matlab使用者。通过观察UF、UB曲线交点位置也可辅助判断气候突变的起始年份。1. 在气候诊断里我们经常遇到这种局面一段 30 多年的降水序列肉眼看上去像有变化但当场说不出它是缓慢上升还是某一年突然跳变。Mann-KendallM-K检验的价值就在这里它基于符号差分不要求样本服从正态分布对异常值也不敏感所以被广泛用于降水、径流、干旱指数序列的趋势检测与突变定位。源码包里通常有三个文件MK_do.m、MannKendallZ.m 和 MMK.m。UF 表示按原始序列正序计算出的标准化统计量UB 表示按逆序计算后翻转到原时间轴的对应统计量UF 与 UB 在置信区间内的交点就是突变发生的时间位置。下面我把这条链路的原理、MATLAB 实现和参数调整完整拆开。2. UF与UB统计量的构造正序累计、逆序回代与交点判据2.1 符号差分与 S 统计量Mann-Kendall 检验的第一步不是直接看均值或方差而是对序列做两两比较任意 ij 时计算 sign(xj−xi)xjxi 取 1xjxi 取 −1相等取 0全部累加得到一个 S 统计量。S 的绝对值越大说明序列整体朝一个方向推进的确定性越强S 为正意味着上升趋势为负意味着下降趋势。标准化时最常用的方差公式是Var(S) n(n-1)(2n5)/18这个公式默认数据中不存在大量并列值。如果序列里有很多相同观测值比如降水被取整到 0.1 mm或者气温只保留一位小数直接套公式会低估方差。标准修正是先把每类并列值的个数 tp 找出来计算 Σ tp(tp−1)(2tp5)/18再从理论方差中减去。这一步对显著性判定影响很大处理不当会把假趋势判成真趋势。标准化后的统计量为S0 时Z(S−1)/√Var(S)S0 时Z0S0 时Z(S1)/√Var(S)。相比直接拿 S 除以标准差这条式子多了连续性修正在小样本下更稳定。对 α0.05 的双侧检验|Z|1.96 就认为趋势显著。与线性回归的趋势检验不同M-K 检验不看斜率大小只统计“升的样本对”和“降的样本对”谁更多。回归会被极端年份拉高斜率M-K 则让每个样本只出一票因此它天然适合非正态、厚尾的降水和径流资料。2.2 UF正序列的累计形式UF 不是对整个序列只算一次 Z而是从第二个数据点开始逐步推进。对每个时刻 k取前 k 个样本计算内部所有成对符号差之和 S_k并按当前子样本容量 k 的方差做标准化UF_k S_k / sqrt(k(k-1)(2k5)/18)这样得到的曲线每个点表示“如果把序列截断到该时刻前段子列是否已经出现统计显著的趋势”。UF 曲线突破 ±1.96 临界线说明截至该时刻上升或下降趋势已经不能被随机波动解释。实际画图时UF 通常是一条随 k 缓慢爬升或下降的折线在长序列中偶尔还会出现多次穿越临界线的情况说明趋势强弱在不同年代有明显差异。2.3 UB反向序列的构造与对齐UB 的计算方法和 UF 完全一致只是输入序列要先反转。把 xn, xn−1, …, x1 当新序列逐点累计最后将结果翻转回原时间轴。这一翻转是为了让 UB 曲线上第 k 个点仍然对应原始第 k 年含义则是“从末端回看截至该年的趋势累积量”。UF 与 UB 交叉时说明正序趋势和逆序趋势发生了方向性的反转。但要特别说明交叉只是必要条件不是充分条件。只有交点落在两条临界线之间时突变点才可信。如果交点出现在临界线之外或两条线在整段区间内反复缠绕属于典型的边界效应不应强行解读成气候突变。2.4 曲线形态速查表UF 与 UB 形态常见结论UF 穿越 ±1.96UB 在临界线内两者在临界线内相交存在显著突变交点附近即突变时刻UF 穿越 1.96UB 未穿始终无交点持续上升趋势无明显突变点UF 与 UB 都在临界线内反复交叉序列以波动为主不建议判突变UF 和 UB 同时集中在临界线同侧整段趋势极强突变定位不可靠这张表在批量处理中非常实用。先用形态过滤一遍再对少量特殊曲线做人工诊断比逐个看图要高效得多。3. MATLAB源码拆解MannKendallZ与MK_do的实现细节3.1 MannKendallZ.m含并列值修正的趋势统计量MannKendallZ.m 负责输出标准化 Z 值用它可以快速判断序列趋势方向、强度和显著性。一个可运行的版本如下function [Z, pval] MannKendallZ(x) % 输入x - 按时间顺序排列的序列向量 % 输出Z - 标准化MK统计量Z0上升Z0下降 % pval - 双侧检验p值 x x(:); x x(~isnan(x)); % 剔除NaN n length(x); if n 8 error(样本量不足MK检验至少需要8个有效数据); end % 计算 S 统计量 S 0; for i 1:n-1 for j i1:n S S sign(x(j) - x(i)); end end % 并列值修正每类取值 tp ul unique(x); ties 0; for k 1:length(ul) tp sum(x ul(k)); if tp 1 ties ties tp * (tp - 1) * (2 * tp 5); end end % 方差与 Z 值 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 pval 2 * (1 - normcdf(abs(Z))); end逻辑说明代码先用双重循环完成所有成对比较再计算并列值的修正项。if S0 分支里减 1S0 分支里加 1目的是对离散 S 分布做连续性修正避免小样本时 Z 值被高估。pval 用 normcdf 近似适用于 n≥8 的常规序列。参数说明x 必须是严格按时间顺序排列的等间隔序列。如果数据是反序输入的Z 的符号会反转整个结论完全变样。NaN 的简单处理是直接剔除但会改变样本量缺失比例超过 10% 时我一般先用 fillmissing 做线性插值或按年均值做距平再进入检验。这段代码在 MATLAB R2019b 及之后的版本都能直接运行不需要 Statistics Toolbox 之外的额外模块。3.2 MK_do.mUF/UB 突变检测的实现MK_do.m 承担突变定位核心是分别计算正向和反向的累计标准化统计量再找交点。核心代码如下function [UF, UB, tc] MK_do(x) % 输入x - 时间序列向量 % 输出UF - 正序列标准化统计量 % UB - 反序列标准化统计量 % tc - 临界线内交点位置数组 x x(:); n length(x); UF zeros(n, 1); S 0; for k 2:n for j 1:k-1 S S sign(x(k) - x(j)); end varS k * (k - 1) * (2 * k 5) / 18; UF(k) S / sqrt(varS); end y flipud(x); % 反向序列 UB zeros(n, 1); S 0; for k 2:n for j 1:k-1 S S sign(y(k) - y(j)); end varS k * (k - 1) * (2 * k 5) / 18; UB(k) S / sqrt(varS); end UB flipud(UB); % 翻转回原时间坐标 % 在临界线内寻找交点 tc []; for k 2:n-1 if (UF(k) - UB(k)) * (UF(k1) - UB(k1)) 0 if abs(UF(k1)) 1.96 abs(UB(k1)) 1.96 tc [tc, k1]; end end end end逻辑说明外层循环 k 决定累计长度内层循环把新增样本 x(k) 与前面 k−1 个样本逐一比较符号差累进 S。这样 UF(k) 对应“前 k 个样本的内部净秩序增量”。反向序列部分对翻转后的序列执行相同过程再把 UB 翻转回来保证横坐标与原始时间对齐。tc 记录的是满足“UF 与 UB 异号”且统计量在临界线内的位置直接对应突变时间点。参数说明如果要求更严格可以把临界线判定从 abs(UF(k1))1.96 改成同时检查 k 和 k1 两个点如果只是初筛可以只保留第一条异号判断交点再人工复核。批量处理时用宽松版论文定稿时用严格版是我比较常用的做法。3.3 绘图与临界线叠加得到 UF 和 UB 之后标准绘图代码如下figure(Color, w); plot(1:n, UF, r-, LineWidth, 1.4); hold on; plot(1:n, UB, b--, LineWidth, 1.4); plot([1 n], [1.96 1.96], k:, LineWidth, 1.0); plot([1 n], [-1.96 -1.96], k:, LineWidth, 1.0); xlabel(时间序号); ylabel(标准化统计量); legend({UF, UB, ±1.96临界线}, Location, best); grid on; if ~isempty(tc) plot(tc, UF(tc), ko, MarkerFaceColor, k); end临界线默认按 α0.05 取 ±1.96α0.01 时替换为 ±2.576α0.10 时替换为 ±1.645。图中黑色圆点标注的是自动找到的交点位置复查时重点看这个点附近两条曲线是否干净利落地交叉还是贴在一起拖了一段距离。拖尾型交叉通常说明趋势是渐变不适合写“突变”。4. 实战降水序列的MK趋势检验与突变点定位4.1 读取数据并运行假设你拿到的数据是源码包里的 house6x6 或 lastd77 这类表格文件第一列是年份后面是站点或变量列。读取并处理一列的流程如下T readtable(lastd77.csv); x T.precip; % 降水量列 yrs T.year; % 年份列 fprintf(NaN count: %d\n, sum(isnan(x))); if any(isnan(x)) x fillmissing(x, linear); end [UF, UB, tc] MK_do(x); [Z, pval] MannKendallZ(x); fprintf(MK Z%.3f, p%.4f, 突变点位置%s\n, ... Z, pval, mat2str(tc));运行后控制台会输出趋势 Z 值和交点位置。若 tc 为空说明在 α0.05 下没有检测到可信突变点。第一次跑通时建议顺手画一下 3.3 节的图确认曲线形态不要只依赖自动输出的数字。4.2 结果解读假设打印结果是MK Z2.31, p0.021, 突变点位置[37]。这说明三件事Z 值为 2.31超过 1.96序列存在显著上升趋势突变发生在第 37 个时间点对应年份可以直接用yrs(37)取出来p 值为 0.021在 0.05 水平下显著。拿到输出后还要再看一眼曲线。交点处两条线如果呈尖锐交叉说明趋势方向在短时间里完成反转突变特征清晰如果两条线在很长区间内贴近后再分开说明转折过程持续了数年写作时更适合用“趋势转折期”而非“突变点”。4.3 常见误用与参数调整检查项处理方式判读建议缺失值fillmissing 插值缺失占比高时先做月距平并列值过多加 10^-4 级微噪或修正方差并列比例超 5% 时谨慎时间间隔不规则resample 重采样不连续序列直接检验会失真季节性周期按月分组或年度聚合混入季节周期会得到伪趋势并列值的问题在降水数据中最常见尤其是取整后的日降水大量 0 值会让 sign 函数大量返回 0UF/UB 曲线变得迟钝。对这类数据我通常先做月或年总量聚合再执行 MK 检验而不是直接对日值跑。样本量也是经常被忽略的一个参数。n 小于 15 时UF/UB 的交点置信度很低至少需要 20 个样本才能保证交点位置稳定经验中 30 个以上样本的结论才适合写进报告或论文。5. 进阶MMK修正检验与多站点批量计算5.1 什么情况下需要 MMK标准 Mann-Kendall 检验要求样本相互独立。降水、径流序列往往存在一阶正自相关今年的值受去年影响有效样本量低于实际长度方差公式算出来的 Var(S) 偏小Z 值偏大原本不显著的序列被误判为显著。MMKModified Mann-Kendall的思路是先估计序列的一阶秩自相关如果显著就按有效样本量放大方差再做标准化从而得到一个更保守的结论。5.2 MMK.m 的实现逻辑MMK 的实现分为三步秩变换、有效样本量估计、修正方差后的 Z 值计算function [Z, var_adj] MMK(x, alpha) % alpha - 显著性水平默认0.05 if nargin 2, alpha 0.05; end x x(:); n length(x); % 1. 计算秩序列的一阶自相关 r tiedrank(x); dr r - mean(r); rho1 sum(dr(1:end-1) .* dr(2:end)) / sum(dr.^2); % 2. 自相关显著性判断近似使用正态临界值 crit norminv(1 - alpha/2) / sqrt(n); if abs(rho1) crit n_eff n / (1 2 * (n - 1) / n * rho1); else n_eff n; end % 3. 计算 S 统计量 S 0; for i 1:n-1 for j i1:n S S sign(x(j) - x(i)); end end varS n * (n - 1) * (2 * n 5) / 18; var_adj varS * n / n_eff; if S 0 Z (S - 1) / sqrt(var_adj); elseif S 0 Z (S 1) / sqrt(var_adj); else Z 0; end end逻辑说明先用 tiedrank 把原始数据转成秩再用秩序列计算一阶自相关 rho1这样能避免原始尺度对相关性估计的影响。当 |rho1| 超过临界值时按 n_eff 调整方差rho1 为正时 n_eff 小于 n修正后方差变大Z 值缩小结论更保守rho1 为负时则相反。修改后的 var_adj 可以直接用于显著性判断。参数说明crit 的取值随 alpha 变化alpha0.05 时n30 的临界值约为 0.358n60 时约为 0.253。实际使用中建议同时输出修正前后的 Z 值便于判断序列自相关对趋势结论的影响程度。5.3 批量计算与筛选站点数量多时建议把 MK_do 和 MannKendallZ 封装成函数在一个循环里统一处理files dir(stations/*.csv); out table(); for i 1:length(files) T readtable(fullfile(files(i).folder, files(i).name)); x T.Value; [UF, UB, tc] MK_do(x); [Z, pval] MannKendallZ(x); if ~isempty(tc) mutation T.Year(tc(1)); % 取第一个交点 else mutation NaN; end out [out; {files(i).name, Z, pval, mutation}]; end writetable(out, MK_summary.xlsx);汇总表生成后先筛趋势显著的站点再对有突变的站点回看曲线。对 MK_do 的 tc 判定还有一个调整技巧想要初筛更宽就保留所有 UF(k) 与 UB(k) 异号的交点想要结论更稳就让交点前后至少 3 个点的 UF/UB 方向保持一致。放松适合批量扫数据收紧适合论文级结论具体阈值需要根据研究区域的序列长度和波动特性来调整。本文还有配套的精品资源点击获取 Mann-Kendall检验原理与MATLAB实现：气候序列趋势与突变检测
简介MK检验是气候诊断与预测中常用的非参数统计方法不依赖样本分布假设适用于判断气候序列是否存在突变及其发生时间也可用于降水、干旱频次等要素的趋势显著性检验。这套Matlab源码包面向气候、水文、环境等专业研究者与学生包含三个脚本文件分别实现主流程、Z统计量计算与改进型M-K检验代码结构清晰注释完整运行后可直接输出UF和UB曲线帮助快速定位突变时刻并支持调整置信水平。压缩包仅3KB轻量易用适合二次开发或嵌入现有研究流程。资源已有3405人学习无论是处理长序列数据还是撰写科研论文都可作为直观的算法参考与可复用模板尤其适合希望弄懂UF与UB含义的Matlab使用者。通过观察UF、UB曲线交点位置也可辅助判断气候突变的起始年份。1. 在气候诊断里我们经常遇到这种局面一段 30 多年的降水序列肉眼看上去像有变化但当场说不出它是缓慢上升还是某一年突然跳变。Mann-KendallM-K检验的价值就在这里它基于符号差分不要求样本服从正态分布对异常值也不敏感所以被广泛用于降水、径流、干旱指数序列的趋势检测与突变定位。源码包里通常有三个文件MK_do.m、MannKendallZ.m 和 MMK.m。UF 表示按原始序列正序计算出的标准化统计量UB 表示按逆序计算后翻转到原时间轴的对应统计量UF 与 UB 在置信区间内的交点就是突变发生的时间位置。下面我把这条链路的原理、MATLAB 实现和参数调整完整拆开。2. UF与UB统计量的构造正序累计、逆序回代与交点判据2.1 符号差分与 S 统计量Mann-Kendall 检验的第一步不是直接看均值或方差而是对序列做两两比较任意 ij 时计算 sign(xj−xi)xjxi 取 1xjxi 取 −1相等取 0全部累加得到一个 S 统计量。S 的绝对值越大说明序列整体朝一个方向推进的确定性越强S 为正意味着上升趋势为负意味着下降趋势。标准化时最常用的方差公式是Var(S) n(n-1)(2n5)/18这个公式默认数据中不存在大量并列值。如果序列里有很多相同观测值比如降水被取整到 0.1 mm或者气温只保留一位小数直接套公式会低估方差。标准修正是先把每类并列值的个数 tp 找出来计算 Σ tp(tp−1)(2tp5)/18再从理论方差中减去。这一步对显著性判定影响很大处理不当会把假趋势判成真趋势。标准化后的统计量为S0 时Z(S−1)/√Var(S)S0 时Z0S0 时Z(S1)/√Var(S)。相比直接拿 S 除以标准差这条式子多了连续性修正在小样本下更稳定。对 α0.05 的双侧检验|Z|1.96 就认为趋势显著。与线性回归的趋势检验不同M-K 检验不看斜率大小只统计“升的样本对”和“降的样本对”谁更多。回归会被极端年份拉高斜率M-K 则让每个样本只出一票因此它天然适合非正态、厚尾的降水和径流资料。2.2 UF正序列的累计形式UF 不是对整个序列只算一次 Z而是从第二个数据点开始逐步推进。对每个时刻 k取前 k 个样本计算内部所有成对符号差之和 S_k并按当前子样本容量 k 的方差做标准化UF_k S_k / sqrt(k(k-1)(2k5)/18)这样得到的曲线每个点表示“如果把序列截断到该时刻前段子列是否已经出现统计显著的趋势”。UF 曲线突破 ±1.96 临界线说明截至该时刻上升或下降趋势已经不能被随机波动解释。实际画图时UF 通常是一条随 k 缓慢爬升或下降的折线在长序列中偶尔还会出现多次穿越临界线的情况说明趋势强弱在不同年代有明显差异。2.3 UB反向序列的构造与对齐UB 的计算方法和 UF 完全一致只是输入序列要先反转。把 xn, xn−1, …, x1 当新序列逐点累计最后将结果翻转回原时间轴。这一翻转是为了让 UB 曲线上第 k 个点仍然对应原始第 k 年含义则是“从末端回看截至该年的趋势累积量”。UF 与 UB 交叉时说明正序趋势和逆序趋势发生了方向性的反转。但要特别说明交叉只是必要条件不是充分条件。只有交点落在两条临界线之间时突变点才可信。如果交点出现在临界线之外或两条线在整段区间内反复缠绕属于典型的边界效应不应强行解读成气候突变。2.4 曲线形态速查表UF 与 UB 形态常见结论UF 穿越 ±1.96UB 在临界线内两者在临界线内相交存在显著突变交点附近即突变时刻UF 穿越 1.96UB 未穿始终无交点持续上升趋势无明显突变点UF 与 UB 都在临界线内反复交叉序列以波动为主不建议判突变UF 和 UB 同时集中在临界线同侧整段趋势极强突变定位不可靠这张表在批量处理中非常实用。先用形态过滤一遍再对少量特殊曲线做人工诊断比逐个看图要高效得多。3. MATLAB源码拆解MannKendallZ与MK_do的实现细节3.1 MannKendallZ.m含并列值修正的趋势统计量MannKendallZ.m 负责输出标准化 Z 值用它可以快速判断序列趋势方向、强度和显著性。一个可运行的版本如下function [Z, pval] MannKendallZ(x) % 输入x - 按时间顺序排列的序列向量 % 输出Z - 标准化MK统计量Z0上升Z0下降 % pval - 双侧检验p值 x x(:); x x(~isnan(x)); % 剔除NaN n length(x); if n 8 error(样本量不足MK检验至少需要8个有效数据); end % 计算 S 统计量 S 0; for i 1:n-1 for j i1:n S S sign(x(j) - x(i)); end end % 并列值修正每类取值 tp ul unique(x); ties 0; for k 1:length(ul) tp sum(x ul(k)); if tp 1 ties ties tp * (tp - 1) * (2 * tp 5); end end % 方差与 Z 值 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 pval 2 * (1 - normcdf(abs(Z))); end逻辑说明代码先用双重循环完成所有成对比较再计算并列值的修正项。if S0 分支里减 1S0 分支里加 1目的是对离散 S 分布做连续性修正避免小样本时 Z 值被高估。pval 用 normcdf 近似适用于 n≥8 的常规序列。参数说明x 必须是严格按时间顺序排列的等间隔序列。如果数据是反序输入的Z 的符号会反转整个结论完全变样。NaN 的简单处理是直接剔除但会改变样本量缺失比例超过 10% 时我一般先用 fillmissing 做线性插值或按年均值做距平再进入检验。这段代码在 MATLAB R2019b 及之后的版本都能直接运行不需要 Statistics Toolbox 之外的额外模块。3.2 MK_do.mUF/UB 突变检测的实现MK_do.m 承担突变定位核心是分别计算正向和反向的累计标准化统计量再找交点。核心代码如下function [UF, UB, tc] MK_do(x) % 输入x - 时间序列向量 % 输出UF - 正序列标准化统计量 % UB - 反序列标准化统计量 % tc - 临界线内交点位置数组 x x(:); n length(x); UF zeros(n, 1); S 0; for k 2:n for j 1:k-1 S S sign(x(k) - x(j)); end varS k * (k - 1) * (2 * k 5) / 18; UF(k) S / sqrt(varS); end y flipud(x); % 反向序列 UB zeros(n, 1); S 0; for k 2:n for j 1:k-1 S S sign(y(k) - y(j)); end varS k * (k - 1) * (2 * k 5) / 18; UB(k) S / sqrt(varS); end UB flipud(UB); % 翻转回原时间坐标 % 在临界线内寻找交点 tc []; for k 2:n-1 if (UF(k) - UB(k)) * (UF(k1) - UB(k1)) 0 if abs(UF(k1)) 1.96 abs(UB(k1)) 1.96 tc [tc, k1]; end end end end逻辑说明外层循环 k 决定累计长度内层循环把新增样本 x(k) 与前面 k−1 个样本逐一比较符号差累进 S。这样 UF(k) 对应“前 k 个样本的内部净秩序增量”。反向序列部分对翻转后的序列执行相同过程再把 UB 翻转回来保证横坐标与原始时间对齐。tc 记录的是满足“UF 与 UB 异号”且统计量在临界线内的位置直接对应突变时间点。参数说明如果要求更严格可以把临界线判定从 abs(UF(k1))1.96 改成同时检查 k 和 k1 两个点如果只是初筛可以只保留第一条异号判断交点再人工复核。批量处理时用宽松版论文定稿时用严格版是我比较常用的做法。3.3 绘图与临界线叠加得到 UF 和 UB 之后标准绘图代码如下figure(Color, w); plot(1:n, UF, r-, LineWidth, 1.4); hold on; plot(1:n, UB, b--, LineWidth, 1.4); plot([1 n], [1.96 1.96], k:, LineWidth, 1.0); plot([1 n], [-1.96 -1.96], k:, LineWidth, 1.0); xlabel(时间序号); ylabel(标准化统计量); legend({UF, UB, ±1.96临界线}, Location, best); grid on; if ~isempty(tc) plot(tc, UF(tc), ko, MarkerFaceColor, k); end临界线默认按 α0.05 取 ±1.96α0.01 时替换为 ±2.576α0.10 时替换为 ±1.645。图中黑色圆点标注的是自动找到的交点位置复查时重点看这个点附近两条曲线是否干净利落地交叉还是贴在一起拖了一段距离。拖尾型交叉通常说明趋势是渐变不适合写“突变”。4. 实战降水序列的MK趋势检验与突变点定位4.1 读取数据并运行假设你拿到的数据是源码包里的 house6x6 或 lastd77 这类表格文件第一列是年份后面是站点或变量列。读取并处理一列的流程如下T readtable(lastd77.csv); x T.precip; % 降水量列 yrs T.year; % 年份列 fprintf(NaN count: %d\n, sum(isnan(x))); if any(isnan(x)) x fillmissing(x, linear); end [UF, UB, tc] MK_do(x); [Z, pval] MannKendallZ(x); fprintf(MK Z%.3f, p%.4f, 突变点位置%s\n, ... Z, pval, mat2str(tc));运行后控制台会输出趋势 Z 值和交点位置。若 tc 为空说明在 α0.05 下没有检测到可信突变点。第一次跑通时建议顺手画一下 3.3 节的图确认曲线形态不要只依赖自动输出的数字。4.2 结果解读假设打印结果是MK Z2.31, p0.021, 突变点位置[37]。这说明三件事Z 值为 2.31超过 1.96序列存在显著上升趋势突变发生在第 37 个时间点对应年份可以直接用yrs(37)取出来p 值为 0.021在 0.05 水平下显著。拿到输出后还要再看一眼曲线。交点处两条线如果呈尖锐交叉说明趋势方向在短时间里完成反转突变特征清晰如果两条线在很长区间内贴近后再分开说明转折过程持续了数年写作时更适合用“趋势转折期”而非“突变点”。4.3 常见误用与参数调整检查项处理方式判读建议缺失值fillmissing 插值缺失占比高时先做月距平并列值过多加 10^-4 级微噪或修正方差并列比例超 5% 时谨慎时间间隔不规则resample 重采样不连续序列直接检验会失真季节性周期按月分组或年度聚合混入季节周期会得到伪趋势并列值的问题在降水数据中最常见尤其是取整后的日降水大量 0 值会让 sign 函数大量返回 0UF/UB 曲线变得迟钝。对这类数据我通常先做月或年总量聚合再执行 MK 检验而不是直接对日值跑。样本量也是经常被忽略的一个参数。n 小于 15 时UF/UB 的交点置信度很低至少需要 20 个样本才能保证交点位置稳定经验中 30 个以上样本的结论才适合写进报告或论文。5. 进阶MMK修正检验与多站点批量计算5.1 什么情况下需要 MMK标准 Mann-Kendall 检验要求样本相互独立。降水、径流序列往往存在一阶正自相关今年的值受去年影响有效样本量低于实际长度方差公式算出来的 Var(S) 偏小Z 值偏大原本不显著的序列被误判为显著。MMKModified Mann-Kendall的思路是先估计序列的一阶秩自相关如果显著就按有效样本量放大方差再做标准化从而得到一个更保守的结论。5.2 MMK.m 的实现逻辑MMK 的实现分为三步秩变换、有效样本量估计、修正方差后的 Z 值计算function [Z, var_adj] MMK(x, alpha) % alpha - 显著性水平默认0.05 if nargin 2, alpha 0.05; end x x(:); n length(x); % 1. 计算秩序列的一阶自相关 r tiedrank(x); dr r - mean(r); rho1 sum(dr(1:end-1) .* dr(2:end)) / sum(dr.^2); % 2. 自相关显著性判断近似使用正态临界值 crit norminv(1 - alpha/2) / sqrt(n); if abs(rho1) crit n_eff n / (1 2 * (n - 1) / n * rho1); else n_eff n; end % 3. 计算 S 统计量 S 0; for i 1:n-1 for j i1:n S S sign(x(j) - x(i)); end end varS n * (n - 1) * (2 * n 5) / 18; var_adj varS * n / n_eff; if S 0 Z (S - 1) / sqrt(var_adj); elseif S 0 Z (S 1) / sqrt(var_adj); else Z 0; end end逻辑说明先用 tiedrank 把原始数据转成秩再用秩序列计算一阶自相关 rho1这样能避免原始尺度对相关性估计的影响。当 |rho1| 超过临界值时按 n_eff 调整方差rho1 为正时 n_eff 小于 n修正后方差变大Z 值缩小结论更保守rho1 为负时则相反。修改后的 var_adj 可以直接用于显著性判断。参数说明crit 的取值随 alpha 变化alpha0.05 时n30 的临界值约为 0.358n60 时约为 0.253。实际使用中建议同时输出修正前后的 Z 值便于判断序列自相关对趋势结论的影响程度。5.3 批量计算与筛选站点数量多时建议把 MK_do 和 MannKendallZ 封装成函数在一个循环里统一处理files dir(stations/*.csv); out table(); for i 1:length(files) T readtable(fullfile(files(i).folder, files(i).name)); x T.Value; [UF, UB, tc] MK_do(x); [Z, pval] MannKendallZ(x); if ~isempty(tc) mutation T.Year(tc(1)); % 取第一个交点 else mutation NaN; end out [out; {files(i).name, Z, pval, mutation}]; end writetable(out, MK_summary.xlsx);汇总表生成后先筛趋势显著的站点再对有突变的站点回看曲线。对 MK_do 的 tc 判定还有一个调整技巧想要初筛更宽就保留所有 UF(k) 与 UB(k) 异号的交点想要结论更稳就让交点前后至少 3 个点的 UF/UB 方向保持一致。放松适合批量扫数据收紧适合论文级结论具体阈值需要根据研究区域的序列长度和波动特性来调整。本文还有配套的精品资源点击获取