简介：这套MATLAB程序包面向水文、气象等领域的科研人员与在校学生，实现Mann-Kendall趋势检验与突变点分析。趋势检验程序用于显著性检验，突变点分析程序负责突变点识别，秩统计量计算程序为其提供支撑；两个模块配合可完成从秩计算、趋势检验到突变识别的完整流程；程序未考虑序列中相等值（结）的情况，当数据存在重复值时，检验结果可能受影响，使用前建议先做相应处理。资源包共5个文件，包含3个M脚本、1个txt说明文档和1份上海气温示例数据，压缩包仅3KB，结构紧凑；说明文档及每个程序头部的参数说明都为快速调用和修改提供了便利；目前已有2488人学习下载。借助自带示例数据，读者能完整运行趋势检验与突变点分析流程，直观理解非参数方法在气候、水文序列诊断中的应用，也可替换为自己的时间序列数据开展科研或课程实验，是一份简洁实用的Matlab参考代码。

1. Mann-Kendall 趋势检验与突变点分析：为什么做水文气候序列分析的人绕不开它

做水文、气象、环境监测数据的长期趋势分析，Mann-Kendall 趋势检验和突变点分析几乎是默认第一选择。对一组年降水、月径流或气温序列，你需要回答两个问题：整体有没有显著变化趋势？如果有，趋势是从第几年开始变向的？Mann-Kendall 检验不要求数据服从正态分布，对缺失值和异常值的耐受度也比参数检验高很多，用它处理多年观测序列比直接拟合一条直线可靠得多。这套 MATLAB 程序把趋势显著性检验和序贯 MK 突变点识别做成了完整可运行的实现，适合水库调度、区域气候评估、土壤墒情这类长序列分析场景。有你要处理的数据，照着下文的步骤可以完整复现。

2. 先讲清原理：S 统计量、方差修正与显著性判断

MK 检验的原理不复杂，每一步都有对应的 MATLAB 代码，但细节直接决定结果可不可信。先理清统计量的构造方式，再谈参数怎么选，后面写代码时才不会把公式抄错。

2.1 S 统计量与方差修正：为什么并列值必须算

MK 检验的核心是统计量 S，它比较序列中每一对样本的大小关系。对长度为 n 的序列 x，S 等于所有满足 i<j 的样本点对做符号比较后累加：

S = ΣΣ sgn(xj − xi)，其中 i<j

sgn 是符号函数：后一个值大于前一个记 +1，小于记 −1，相等记 0。如果 S 显著大于 0，说明序列整体上升；显著小于 0，说明整体下降。举一个直观例子：假设序列是 [3, 5, 4, 6, 8]，从第一对 (3,5) 开始，5>3 记 +1；(3,4) 再记 +1；(3,6) 记 +1；(3,8) 记 +1；接着比较 (5,4) 记 −1；(5,6) 记 +1，依此类推。把所有点对符号加总，正值占多数时序列自然呈上升态势。如果数据是随机波动的，正负号大致抵消，S 会落在 0 附近。

S 本身不能直接用于假设检验，因为它的期望和方差只由样本长度和并列结构决定，所以要标准化成 Z 统计量。这里最容易被忽略的是并列值（tied groups）。如果数据里连续出现重复值——降水量按 0.1 mm 取整、气温按 0.1℃ 记录，或者缺测位置被填了 0——S 的方差会被高估，而方差高估会让 Z 偏小，可能导致“本来有趋势却被判定为不显著”。反过来，如果处理不当，某些情况下也会放大 Z。修正方式是从总方差里减掉并列组贡献：

Var(S) = [n(n−1)(2n+5) − Σ tp(tp−1)(2tp+5)] / 18

其中 tp 是第 p 组并列数据的个数。比如一组数据里出现三个相等的 12，tp 就是 3。MATLAB 里用 unique 配合 accumarray 可以快速统计分组计数，后面的代码会直接用。

为什么强调这个修正？因为 MK 检验本质上是排列比较，样本中所有点对的符号信息已经全部包含在 S 里。如果忽略了并列修正，对连续型观测数据影响不大，但一旦出现取整、缺测填 0 或浮点精度截断，Z 值偏差立刻显现。很多现成脚本只写了 S 和 Var(S) 的简单公式，对干净数据能跑，对有大量重复值的监测数据就会翻车。

2.2 从 S 到 Z：显著性水平、临界值和假设检验口径

标准化后的 Z 统计量写法如下：

当 S > 0 时：Z = (S − 1) / √Var(S)
当 S = 0 时：Z = 0
当 S < 0 时：Z = (S + 1) / √Var(S)

±1 的连续校正项是为了让离散的 S 分布向连续正态分布过渡更自然。显著性水平 α 一般取 0.05，比较 |Z| 是否超过 1.96。常用水平对应的临界值如下表：

我的选择习惯是：常规年序列数据用 α=0.05；如果样本量偏小，序列只有二十几年，用 α=0.10 很容易出假显著，这时候我会再切到 α=0.01 做一次交叉验证。水文气候领域很多报告要求同时给出 0.05 和 0.01 两档结果，所以程序里最好把 α 做成可配置参数而不是硬编码。

还要说清楚一个容易混淆的点：|Z| 超过临界值时，我们拒绝“序列无趋势”的原假设，说序列存在显著单调趋势。Z 只回答“有没有”，不回答“趋势强度多大”。趋势强度要用 Sen 斜率或线性倾向值来评估，那是另一组计算。实际写报告时不要把 Z 值直接当成变化幅度，这是新手最常见的一个误读。

3. 落到 MATLAB：趋势检验和序贯 MK 的完整实现

这一章给出可以直接运行的代码。先做数据预处理，再写趋势检验主函数，最后实现序贯 MK 的 UF/UB 曲线计算。代码都按 MATLAB 脚本组织，复制到本地就能跑。

3.1 数据预处理：读取、缺测过滤与列向量化

拿到原始数据后先做两件事：读入并统一成列向量，然后过滤 NaN。MK 检验的符号比较在遇到 NaN 时会传播，任何一步出现 NaN，累计的 S 都可能变成损坏值。常见做法是用 readmatrix 或 load 读数据，然后对存在缺测的序列做两种处理：缺测占比低于 5% 且是连续时间序列，用线性插值补上；占比高或集中在某一段，直接删除这段再分析。直接填 0 是很危险的操作，因为你把缺测变成了并列值，会污染方差修正项。

一个稳健入口写法如下，函数内部先做非 NaN 过滤：

function [x_clean, t_clean] = prepare_series(x, t) % 输入: x 原始观测序列, t 对应时间轴 % 输出: 剔除NaN后的序列和时间轴 valid = ~isnan(x); x_clean = x(valid); t_clean = t(valid); % 若剔除后样本长度过短, 抛错提示 if length(x_clean) < 10 error('清洗后样本不足10个, MK检验结果不可信'); end end

逻辑说明：把观测值和时间轴按同一掩码过滤，避免后续绘图或定位突变点时时间轴对不上。10 个样本的门槛来自 MK 检验最早的应用约定——样本少于 10 时 S 的精确分布仍远离正态，直接做正态近似容易出错，需要退回到查表或精确概率计算。

参数说明：这里的 10 是保守下限。如果是逐月数据做 40 年序列，有效样本接近 480 个，完全没问题；如果缺测清洗后只剩下 8 个点，程序直接拒绝执行比硬出结果更负责。

3.2 趋势检验主函数：完整可运行代码

下面这段函数实现了完整 MK 检验，包含并列值修正、Z 计算、p 值估计和结论判定。调用方式是：

x = load('series.txt'); [S, VarS, Z, p, trend] = mk_trend_test(x, 0.05);

function [S, VarS, Z, p, trend] = mk_trend_test(x, alpha) % Mann-Kendall 趋势检验主函数 % 输入: % x - 列向量观测序列(N x 1) % alpha - 显著性水平, 默认0.05 % 输出: % S - 原始统计量 % VarS - 方差(含并列修正) % Z - 标准化统计量 % p - 双侧p值 % trend - 趋势结论字符串 if nargin < 2, alpha = 0.05; end x = x(:); % 强制列向量 x = x(~isnan(x)); % 剔除NaN n = length(x); if n < 10, error('样本量不足10, 需改用精确检验'); end S = 0; for i = 1:n-1 for j = i+1:n diff_val = x(j) - x(i); if diff_val > 0 S = S + 1; elseif diff_val < 0 S = S - 1; end % diff_val == 0 时 S 不加不减 end end % 并列值分组统计: 用unique统计每组出现次数 [~, ~, ic] = unique(x); counts = accumarray(ic, 1); tie_adj = sum(counts .* (counts - 1) .* (2 .* counts + 5)); VarS = (n * (n - 1) * (2 * n + 5) - tie_adj) / 18; 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))); if p > 1, p = 1; end if abs(Z) > norminv(1 - alpha / 2) if Z > 0 trend = '显著上升'; else trend = '显著下降'; end else trend = '无显著趋势'; end end

逻辑说明：主流程是双层循环做全部样本点的符号比较，计算复杂度 O(n²)。水文气候序列长度一般几百到几千，跑起来很快；如果到十万级，可以改用归并排序思想统计顺序对数量来加速，但常规用不到。S 算完后用 unique 和 accumarray 统计并列组，这一步最容易踩坑。

参数说明：diff_val > 0记 +1、<0记 −1、等于 0 不记，这是 MK 检验对“相同观测值”的处理方式；tie_adj就是前面方差公式里的 Σ 项；输出 Z 的同时给了 p 值，比单独看 Z 更方便判断显著程度。norminv是临界值的编程写法，和直接写 1.96 效果一样，好处是 alpha 可配置。注意p = 2 * (1 - normcdf(abs(Z)))在 Z 接近 0 时可能略大于 1，所以加了上限保护。

3.3 序贯 MK：UF/UB 曲线计算与突变点定位

趋势检验回答“有没有趋势”，突变点分析回答“什么时候变向”。实现方式是序贯 MK：对序列的每个前缀计算一次 MK 统计量，得到正序列 UFk；把原序列倒序后同样计算，再取负并倒序，得到反序列 UBk。两条曲线在置信区间内交叉的位置，就是候选突变点。完整实现：

function [UFk, UBk] = sequential_mk(x) % 序贯 Mann-Kendall 突变检验 % 输入: x 列向量 % 输出: UFk, UBk 统计量曲线(长度与x相同) x = x(:); n = length(x); UFk = zeros(n, 1); for k = 2:n Sk = 0; for i = 1:k-1 if x(k) > x(i) Sk = Sk + 1; elseif x(k) < x(i) Sk = Sk - 1; end end E = k * (k - 1) / 4; V = k * (k - 1) * (2 * k + 5) / 72; UFk(k) = (Sk - E) / sqrt(V); end % 倒序序列重算 x_rev = flipud(x); UB_rev = zeros(n, 1); for k = 2:n Sk = 0; for i = 1:k-1 if x_rev(k) > x_rev(i) Sk = Sk + 1; elseif x_rev(k) < x_rev(i) Sk = Sk - 1; end end E = k * (k - 1) / 4; V = k * (k - 1) * (2 * k + 5) / 72; UB_rev(k) = (Sk - E) / sqrt(V); end UBk = flipud(-UB_rev); end

逻辑说明：UFk 在时间点 k 处，比较的是第 k 个观测与其之前所有观测的大小关系。序列持续上升时，UFk 倾向为正且不断上穿正值区域；持续下降时倾向为负。UBk 等于倒序结果取负并翻转回来，和 UFk 的起点不同但尺度相同，所以当趋势真正发生转向时，两条曲线会在转向时刻附近交错。

参数说明：这里的方差公式是简化版，因为它针对每个前缀内的顺序统计量，不含并列修正。如果数据取整严重，建议把并列修正加进 UF/UB 计算，或者至少意识到边界线附近的交点可信度要打折扣。绘图时再叠加 ±1.96 的临界线。绘图部分：

figure; plot(1:n, UFk, 'b-', 'LineWidth', 1.2); hold on; plot(1:n, UBk, 'r-', 'LineWidth', 1.2); plot([1 n], [1.96 1.96], 'k--'); plot([1 n], [-1.96 -1.96], 'k--'); xlabel('时间序列索引'); ylabel('统计量 UFk / UBk'); legend({'UFk', 'UBk', '临界线'}, 'Location', 'best'); grid on;

图形看的就是红蓝两线在临界带内外的交叉情况。但肉眼看图只适合初筛，批量分析时需要一个能自动提取交点的辅助判断，这是下一章的重点。

4. 突变点解读与自动识别：别把“相交”直接当成“突变”

拿到 UF/UB 图之后，第一反应往往是“找交叉点”。但交叉点本身不等于突变点，必须结合显著性水平和曲线形态一起判断，否则很容易把随机波动当成趋势转折。

4.1 UF/UB 双曲线：交叉点、越线与置信区间

画完图后，我第一眼先看 UFk 有没有越线。如果整条 UFk 都在 ±1.96 以内，说明序列整体上连趋势显著性都没通过，这种情况下即便 UFk 和 UBk 在某些位置相交，最多只是随机波动，谈不上突变。值得分析的形态是：UFk 持续上升并越过 +1.96 线一段时间，UBk 也在负向区域快速下穿，两条线在临界区间内交叉。这个时候交叉点对应的横坐标才是趋势方向转变的时刻。

为什么强调“临界区间内”？因为如果交叉发生在临界线以上的区域，说明那是极端高值点造成的瞬时反转，不是稳定的突变。交叉点的位置也有讲究：交叉发生在索引中点，代表突变时刻。实际分析时，我一般会在报告里同时给出索引号和对应的日历年份，避免后续引用时出现错位。

如果 UFk 越线后一直保持在临界线外，但 UFk 与 UBk 没有交叉，典型情况是序列从分析起点就已经处在上升通道中，突变起点早于观测范围。这种结果不能支撑“序列在某年发生突变”的结论，只能说明序列相对历史均值发生了显著偏移。

4.2 交点识别策略：多个交点时如何筛选

气候水文序列经常出现多个交点。筛选原则我总结成以下几点：取第一个满足两条曲线在连续两个时间点之间符号相反的交叉；这个交叉必须落在 ± 临界线范围内；交叉之后两条线持续分离至少数个时间步，而不是立刻再次相交；如果序列存在明显周期性，先做季节差分或聚合到年尺度再判断。

这些规则用代码实现，核心就是找相邻点的符号变化并加置信区间条件：

function cp = find_mk_change_points(UFk, UBk, alpha) % 自动寻找UFk与UBk的显著交叉点 crit = norminv(1 - alpha / 2); cp = []; n = length(UFk); for i = 1:n-1 % 两线在相邻两点间位置关系发生翻转, 说明在i和i+1之间有交叉 d1 = UFk(i) - UBk(i); d2 = UFk(i+1) - UBk(i+1); if d1 * d2 < 0 % 交叉点须落在临界带附近, 至少一个端点在线带范围内 if (abs(UFk(i)) < crit) && (abs(UBk(i+1)) < crit) cp(end+1, 1) = i; elseif (abs(UFk(i+1)) < crit) && (abs(UBk(i)) < crit) cp(end+1, 1) = i; end end end end

逻辑说明：d1和d2是两线间距在相邻索引上的符号，一正一负说明中间存在交叉。临界区间约束用的是“交叉点左右两端点都不超出临界线”，比只要求交叉点坐标落在区间内更严格。

参数说明：这里的crit随 alpha 自动变化；cp返回的是交叉前那个时间索引。实际使用中还要加第二个约束——两线交叉后在后续 10 个点内不再反向相交，可以通过对cp做二次过滤实现。这个辅助函数的价值在于：分析 40 年以上的年序列时，肉眼找交叉点费劲还容易漏，批量处理几十个测站更是如此。

多个交叉点的后处理：做一次分组过滤，把间隔小于 5 个时间步的交叉点合并成最早的那个，被合并的往往只是小幅震荡。这些规则融入分析脚本后，突变定位基本稳定，不需要每次再人工指认。端部交叉要单独处理：如果交叉点落在序列前 10% 或后 10% 范围，说明趋势变化发生在观测起点之前或尚未完成，这种结果只可作为趋势方向的补充证据，不建议直接写进报告当成突变时刻。

每个交叉点判定为突变前，最好还做一次滑动 t 检验佐证。做法是取交叉点左右各 m 个样本构成两个窗口，计算 t 统计量，t 值显著才能放到结论里。m 常见取 5 到 10，样本太少时滑动 t 不稳定，样本太多时突变信息被抹平，两个方向都试一下再定。

5. 避坑记录：MK 检验和突变点分析最容易翻车的五处

5.1 三个高频运行与数值问题

现象 1：代码跑完，Z 显示 NaN，p 值也是 NaN，结论全是 NaN。原因：输入序列里有 NaN 没有过滤，双循环里的符号比较把 NaN 传播开了。S 统计量无论哪一步遇到 NaN，都会污染后续所有累加。解决：在函数入口统一过滤并检查剩余长度。缺测少用插值补，缺测多直接删除对应时段。批量跑十几组传感器数据时最容易遇到一条数据流从某个点开始全是 NaN，结果整批结果里出现一条 NaN，检查时先看清洗函数。

现象 2：肉眼看着趋势平缓，结果 Z 高达 2.5，p 值小到离谱。原因：数据里存在大量重复值。最常见的是缺测位置被填 0，或观测本身按整数毫米取整，导致并列值分组计数巨大，方差修正项异常缩小，Z 被放大。解决：统计序列里重复值的占比。占比超过 10% 时优先怀疑两类问题：一是缺测填了固定值，二是取整精度过粗。填 0 缺测要改成 NaN 再过滤；取整问题考虑恢复更高精度，或先做月/年聚合再检验。

现象 3：样本量少于 10 个，程序照样输出漂亮的 Z 值。原因：程序没有做样本量门槛检查，直接用正态近似公式。小样本下 S 的真实分布接近离散分布，Z 的分布远不是标准正态。解决：头铁的做法是加样本量检查，n<10 直接报错；温和的做法是输出 S 符号用于方向判断，并注明这个结论不具备统计显著性。正规报告里少于 15 年的序列做 MK 检验基本都会被审稿人质疑，分析前注意样本量。

5.2 两个高频结果误读问题

现象 4：UFk 和 UBk 画出来交叉了七八处，每处看起来都像突变。原因：序列本身波动幅度大或存在明显季节周期，把逐月数据直接塞进序贯 MK，月际振荡制造了大量假交叉点。解决：先把序列聚合到年尺度再跑突变分析。年序列也不是绝对安全，如果每年内仍有脉冲式高值，可先做 3 年滑动平均预处理再做 MK，但预处理方式必须在报告里写清楚，否则审稿人第一个质疑就是数据修改。

现象 5：找到的突变点和文献里差了好几年。原因：常见原因是时间轴没对齐，或者两套数据用了不同采样口径。比如同一条水文资料，一份是汛期统计得到的年最大值序列，另一份是全年年平均值序列，两者突变年份天然不可能完全一致。解决：分析前统一时间口径，检查数据的开始年份和样本量是否一致。报告突变年份时，应写作“第 k 个时间点”并同时给出对应的日历年份，避免索引错位。交叉点落在两个索引之间时，默认取较早索引是一种惯例，但也要注明。

6. 最后一步验证：合成序列与假阳性测验，跑通前别下结论

6.1 合成序列验证

一套新到手的 MK 代码，我习惯先用已知答案的合成序列跑一遍。构造一个长度 100 的序列：前 50 个点线性上升叠加小幅随机噪声，第 51 点起整体下降并持续 50 个点。这个序列对应“某变量在中间年份从升转降”的典型突变模式。跑完序贯 MK 后，交叉点应落在第 50 到第 52 个索引之间；UFk 在前半段稳定超过 +1.96，后半段下穿负侧。

如果程序返回的交叉点不在这个区间，优先怀疑 UFk/UBk 的计算顺序和翻转有没有出错。UBk 的构造非常容易写反：先倒序计算、再取负、再翻转，三步顺序错一步，曲线形态就会完全变形。我的排查习惯是先打印 UFk 和 UBk 的前五个值，和手算结果对一遍再批量做分析。

趋势检验部分也可用相同序列验证：前半段单独做 MK 检验，Z 应显著为正；后半段 Z 应显著为负。这个验证跑通后，程序基本上就没有方向性错误了。

6.2 白噪声假阳性测验

第二个验证是假阳性控制。用 randn 生成 1000 个长度 100 的标准正态白噪声序列，每个都跑一遍 mk_trend_test。理论上假阳性率应接近设定 α=0.05，即 50 个左右序列被误判为显著趋势。实际跑出来在 30 到 70 之间都算正常，这能验证临界值和 p 值计算没有严重偏差。

rng(42); n = 100; trials = 1000; alpha = 0.05; hits = 0; for t = 1:trials x = randn(n, 1); [~, ~, Z, p, ~] = mk_trend_test(x, alpha); if p < alpha hits = hits + 1; end end fprintf('假阳性率: %.2f%%\n', hits / trials * 100);

逻辑说明：这里的随机数据满足 MK 检验的“无趋势”原假设，被错误拒绝的次数比例就近似 α。跑一次白噪声检验成本很低，但对代码可信度是很好的背书。

参数说明：randn 生成的都是独立同分布白噪声；n=100 对应常见序列长度；trials 取 1000 是平衡运行时间和统计精度的选择。如果这个比例和 5% 偏差过大，说明 Z 计算或并列修正写错。

从那以后，我每次拿到一批新的水文站数据或换一套复现代码，都会强制先跑这两个验证，再进正式分析。最典型的一次，某套现成脚本跑白噪声假阳性率高达 18%，查下来是并列值修正漏了取整数据的分组统计。希望这些记录能帮到你。

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