简介面向气象、水文、环境等领域研究人员及MATLAB初级用户这套MKMann-Kendall趋势检验源码可用于时间序列的趋势与突变点检测。资源共3个MATLAB脚本文件压缩包体积仅2KB内容聚焦非参数统计检验流程覆盖S统计量与Z值计算、显著性水平判断及突变分析辅助功能可直接在MATLAB中调用。已有7774人学习下载参考价值得到验证。三个脚本分工明确分别承担主检验、趋势统计与突变点可视化支持适合初学者对照描述中的步骤逐一理解实现原理也便于有经验的开发者在此基础上扩展。特别提醒MK适用于线性趋势检测非线性场景需谨慎使用这套源码能帮助快速跑通标准流程是理解该方法的轻量入门工具。1. 用MK趋势检验分析时间序列S和Z值才是关键做气象和水文监测的人应该都有体会拿着十年的逐月数据最想回答的问题不是“增幅有多大”而是“这个上升到底是真的还是随机波动”。线性回归能给出斜率和R方但数据一旦带着几个极端年份回归结果就会跟着异常值跑偏。MK趋势检验Mann-Kendall test恰恰绕开了这个坑它只比较数据两两之间的高低关系不要求序列服从正态分布遇到野值也更抗干扰。这套MATLAB源码包里的核心文件分工比较明确mk.m负责最底层的统计量计算trendMK.m做趋势检验主流程而MannKendall.m更像是突变检验版本输出的UF/UB曲线用于定位突变点。下面从S统计量的构造开始拆再用MATLAB代码把趋势、显著性和突变点三段流程全部跑通。2. 秩次、S统计量与方差修正MK检验的非参数底座2.1 两两比较S统计量到底在数什么MK检验的基础操作是把时间序列里所有数据对拿出来比较晚出现的值比早出现的值更大还是更小。对于一个长度为n的序列一共有n(n-1)/2个数据对对每个数据对赋值sign(x_j - x_i) 1 当 x_j x_i sign(x_j - x_i) -1 当 x_j x_i sign(x_j - x_i) 0 当 x_j x_i其中i j。把所有符号加起来就得到S统计量。如果上升趋势存在晚期的值普遍大于早期的值加出来的S是正数下降趋势则反过来。S接近于0不代表没有波动而是说明上涨和下跌的配对数量大致抵消没有系统性方向。这里要注意i和j的顺序不能写反。早期数据放前面晚期数据放后面一旦反了S的正负号就会翻转趋势方向判断直接出错。在MATLAB里如果写双重循环一般习惯是for i 1:n-1 for j i1:n S S sign(x(j) - x(i)); end end这段代码的边界条件是外层循环到n-1内层从i1开始恰好覆盖所有i j的配对。sign函数返回1、-1或0不需要手动判断分支。2.2 方差公式与平局修正为什么不能用S直接比大小S本身的大小受样本量影响很大n从20变到50S的数值范围会成倍扩大所以直接用S判断显著性没有意义。把S标准化成Z值前需要知道在“没有趋势”这个零假设下S的方差。理论推导给出Var(S) n(n-1)(2n5) / 18但如果数据里存在相等数值也就是平局tie方差会被高估。常见修正是引入平局修正项设第t组相同值有t_i个修正后的方差为Var(S) (n(n-1)(2n5) - Σ t_i(t_i-1)(2t_i5)) / 18从工程实现角度我一般先统计平局分布再决定用哪个公式。下面的表格总结了三种情况的处理方式数据情况方差计算常见场景无平局n(n-1)(2n5)/18连续观测值如气温、水位存在平局修正公式降水量取整、等级评分数据样本量小于8建议用精确分布或查表短序列用MATLAB实现修正方差时可以用unique和histcounts统计相同值的个数function varS mk_var(x) n numel(x); [vals, ~, ic] unique(x); counts accumarray(ic, 1); varS n*(n-1)*(2*n5)/18; % 对平局做修正 t counts(counts 1); if ~isempty(t) varS varS - sum(t.*(t-1).*(2*t5))/18; end endunique取出所有不同的观测值accumarray统计每个值出现次数。最后一行减去平局带来的方差虚增部分。这个修正项在降水量这类整数型数据里尤其重要一个值重复出现五六次很常见。2.3 为什么不是线性回归秩次带来的抗差异线性回归对异常值敏感因为平方误差会给极端值更高权重。而MK检验每个数据对只贡献1、-1或0一个异常值最多影响n-1个配对影响幅度有限。换句话说把最大观测值放大10倍MK的S值不会变但回归斜率会被明显拉偏。这就是环境监测里MK比回归更常用的原因。不过需要说清楚MK检验检测的是“单调趋势方向”不是“趋势幅度”。如果你需要量化每十年上升多少得配合Sen斜率估计trendMK.m源码包里有对应实现时可以直接用。3. 用MATLAB实现S、Z和p值trendMK.m从零到可用3.1 数据预处理排序、NaN与缺失值MK检验要求数据按时间顺序排列但很多Excel导出的数据会混入缺失值或倒序行。在调用trendMK.m之前先做两步清理% 删除NaN并确保按时间递增排序 t t(~isnan(x)); x x(~isnan(x)); [t, idx] sort(t); x x(idx);isnan找出缺失位置sort按时间重排。如果缺失值超过总长度的20%不建议直接补插值因为插值会改变数据对的秩次关系让S统计量失真。宁可把缺测段截掉也不要让补出来的假数据影响显著性判断。3.2 S、Z和双侧p值核心计算代码trendMK.m的完整逻辑可以拆成三层计算S计算方差再做标准化。这里给出一个不依赖任何额外工具箱的MATLAB实现function [Z, p, S, varS] trendMK(x) n numel(x); if n 8 error(样本量小于8时正态近似不成立请查阅临界值表); end S 0; for i 1:n-1 for j i1:n S S sign(x(j) - x(i)); end end varS mk_var(x); if S 0 Z (S - 1) / sqrt(varS); elseif S 0 Z (S 1) / sqrt(varS); else Z 0; end p 2 * (1 - normcdf(abs(Z))); end代码里最容易被忽略的是S标准化时的连续性修正。S大于0时减1S小于0时加1这个操作把离散的S值映射到连续的正态分布上能让p值更接近精确结果。mk_var就是上一章的方差修正函数。p 2 * (1 - normcdf(abs(Z)))得到双侧p值对应“上升或下降趋势是否显著”的判断。如果你的MATLAB没有Statistics and Machine Learning Toolbox可以用erfc替代p erfc(abs(Z) / sqrt(2));erfc是MATLAB基础函数不需要额外工具箱。这一点在只装了MATLAB本体、还没来得及装工具箱的机器上很管用。3.3 显著性判定p值、临界Z值怎么用Z值的正态近似给出一个简洁的判定规则在95%置信水平下|Z| 1.96则认为趋势显著。工程上更推荐直接用p值因为p值可以告诉你显著到什么程度。置信水平临界Z值结论90%1.645较弱趋势证据95%1.960常用显著性标准99%2.576强显著证据比如返回Z 2.31p 0.021。p小于0.05方向为正结论就是“该时间序列呈显著上升趋势”。如果Z 1.80p 0.072只能说有上升倾向不能下显著结论。很多论文里写“趋势不显著”时指的就是p 0.05。3.4 调用示例与输出解读假设有一份2013到2022年的年均水位数据存放在water_level.csv里data readtable(water_level.csv); x data.level; [Z, p, S, varS] trendMK(x); fprintf(S%d, Z%.3f, p%.4f\n, S, Z, p); if p 0.05 Z 0 disp(显著上升); elseif p 0.05 Z 0 disp(显著下降); else disp(无明显趋势); end这里的输出会同时给出统计量和业务结论。S值是原始秩次统计量Z值和p值才是论文里需要报告的。mk.m在源码包里的角色就是把S和方差的计算独立封装trendMK.m负责组装流程这种分层对二次开发很友好换成突变检验时只需要替换标准化方式。4. UF和UB曲线用MannKendall.m定位突变点4.1 趋势检验与突变检验的差异趋势MK检验回答“整体有没有单调方向”突变MK检验回答“在哪个时间点附近序列的统计特征发生跳变”。两者的数学底座相似但处理对象不同。突变检验用的是顺序统计量UF和逆序统计量UB它们分别从序列正方向和反方向计算累积统计量。UF序列的具体做法是对每个前缀子序列x(1:k)计算一次类似S的统计量然后标准化。UB序列则把整段数据倒过来再做一遍最后取负并反向排列。两条曲线在置信区间边界内出现交点对应位置就是可能的突变点。4.2 UF和UB的MATLAB实现下面这个函数直接对每个前缀子序列计算S并转换为UF值逻辑清晰但运算量稍大适合样本量在几百量级的序列function [UF, UB] mk_mutation(x, alpha) n numel(x); uf zeros(1, n); ub zeros(1, n); % 顺序计算UF for k 2:n S 0; for i 1:k-1 for j i1:k S S sign(x(j) - x(i)); end end E k * (k - 1) / 4; V k * (k - 1) * (2 * k 5) / 72; uf(k) (S - E) / sqrt(V); end % 逆序计算UB最后取负反转 rev fliplr(x); for k 2:n S 0; for i 1:k-1 for j i1:k S S sign(rev(j) - rev(i)); end end E k * (k - 1) / 4; V k * (k - 1) * (2 * k 5) / 72; ub(n - k 1) -(S - E) / sqrt(V); end UF uf; UB ub; end注意UF和UB的计算里期望E用的是k(k-1)/4方差V用的是k(k-1)(2k5)/72这和整段序列的方差公式不同因为这里针对的是前缀长度k的累积统计量。UB的取负和反向是突变检验的关键两条线如果不做这个处理图形会完全错位。4.3 绘制突变点分析图拿到UF和UB后还要画两条临界线作为显著性边界。95%置信水平对应的临界Z值是1.96因此alpha 0.05; crit norminv(1 - alpha/2); figure; plot(1:n, UF, b-, LineWidth, 1.2); hold on; plot(1:n, UB, r--, LineWidth, 1.2); plot(1:n, ones(1,n)*crit, k:); plot(1:n, ones(1,n)*(-crit), k:); grid on; legend(UF, UB, 1.96, -1.96, Location, best); xlabel(时间序号); ylabel(统计量);代码中norminv(0.975)返回1.96用作上下边界。水平的两条黑虚线构成置信区间。判定规则是UF和UB在临界线之间出现交叉该交叉点就是突变点如果交点在临界线之外说明突变幅度过大原始数据的跳变甚至超出了正常波动范围。绘制完成后把交叉点所在的时间序号对应回具体时间再结合业务数据判断。比如某流域径流序列在2003年前后UF和UB交叉那一年通常有大坝蓄水、土地利用变化或重大气候事件可对应。这里要强调的是MK突变检验只能给“候选时间窗”不能替代机理分析。图形特征含义UF持续上升并超过上界序列存在显著上升趋势UF下降并跌破下界序列存在显著下降趋势UF与UB交叉于上下界内交叉点为潜在突变点UF与UB交叉于上下界外突变剧烈需要复核数据5. 显著性检验的可复现边界随机重排、样本量和自相关排查5.1 样本量小时Z的近似会失准MK检验的正态近似在n小于8时基本不可靠。如果你的序列只有六七个年份不要直接用normcdf算p值。常见做法是查MK检验精确临界值表或者用蒙特卡洛置换检验。置换检验的思路是把原始数据的位置顺序随机打乱一万次每次计算Z得到一个无趋势条件下的Z分布再看原始Z落在哪个分位。这个做法不依赖正态近似样本小一点也能给出经验p值。5.2 自相关会让MK太容易显著水文气象序列常有一阶自相关也就是今年的值受去年影响。存在正自相关时MK检验的方差公式会低估真实波动导致Z值偏大容易把随机波动误判为显著趋势。预处理时可以先拟合AR(1)模型取残差做MK检验也可以用去趋势预白化先做线性回归去掉趋势对残差估计自相关系数再用修正后的等效样本量调整方差。后者的缺点是把线性趋势带回检验过程和MK的非参数初衷有一点冲突但在应用研究中很常用。5.3 用随机重排拿到经验p值最后分享一个验证显著性的通用技巧用randperm打乱数据顺序重复计算Z得到经验分布。nBoot 1000; Zboot zeros(nBoot, 1); for b 1:nBoot xb x(randperm(n)); [Zb, ~, ~, ~] trendMK(xb); Zboot(b) Zb; end p_emp mean(abs(Zboot) abs(Z_obs));randperm(n)每次生成新的排列相当于把时间顺序彻底破坏原本的趋势信号被抹掉。mean(abs(Zboot) abs(Z_obs))计算的是打乱后出现比原始Z更极端值的比例这个比例就是经验p值。比较经验p值和trendMK返回的正态近似p值如果差异很大说明原始数据的自相关或异常结构在干扰显著性判断优先回去检查数据预处理而不是改置信水平。本文还有配套的精品资源点击获取