滑动t检验的Matlab实现：气候水文突变检测与判读指南

做气候或水文序列突变检测的，十有八九开口就是Mann-Kendall检验。MK确实好用，但真到要判定具体哪一年发生突变的时候，它的UF/UB曲线经常给你画出一大片交叉区，反而让人犯难。相比之下，滑动t检验的思路朴素得多——把序列在每个时刻劈成两段，检验前后两段的均值差异是否显著。这个办法的结果极其直观：每个年份对应一个t统计量，超过临界值的位置就是突变可能发生的年份，还能顺带看出突变的方向和强度。本文就把这套东西的Matlab实现完整拆一遍，从原理到代码，再到判读方法和若干个实际使用中才碰得到的坑，争取让你看完就能自己跑起来。

1. 滑动t检验到底在做什么：原理与适用边界

1.1 一句话讲透原理

滑动t检验（Sliding t-test）的思想可以用一句话概括：对长度为N的序列，把每个时刻i当作潜在突变点，取该点之前n个样本作为第一组、之后n个样本作为第二组，对这两组做一次经典的两独立样本t检验。检验的原假设是“突变点前后两段均值没有显著差异”，如果t统计量落入拒绝域，就认为该点前后的均值发生了显著变化，也就是发生了均值型突变。

t统计量的形式非常标准：

t = (mean(seg1) - mean(seg2)) / (sp * sqrt(1/n1 + 1/n2))

其中sp是两段合并标准差，计算公式是：

sp = sqrt(((n1-1)*var(seg1) + (n2-1)*var(seg2)) / (n1+n2-2))

这个统计量服从自由度为n1+n2-2的t分布。由于滑动t检验在大多数应用中都会取对称窗口，也就是n1=n2=n，所以自由度就是2n-2。给定显著性水平alpha之后，临界值可以从t分布表查，也可以用Matlab的tinv函数一行算出来。

需要注意一个容易搞反的符号问题：t统计量的分子是mean(seg1)-mean(seg2)，所以t>0说明前段均值更高，意味着序列在i点附近发生了均值下降突变；t<0说明后段均值更高，发生了均值上升突变。很多人跑完结果只看绝对值，忽略了符号方向，导致对突变方向的描述写反，这是写论文时很常见的低级错误。

1.2 为什么不用MK检验替代它

经常有人问我：MK检验那么多论文在用，为什么还要单独实现滑动t检验？我的回答是：两种方法解决的问题不完全一样。

MK检验从本质上说是一种基于秩的非参数趋势检验，它的UF/UB曲线交叉定位突变点的方法虽然流行，但实际效果并不稳定。我处理过不少站点数据，UF和UB曲线常常在一段连续几年内反复交叉，根本无法给出唯一的突变年份。而滑动t检验把每个时间点当作候选突变点逐一检验，输出的是一整条t统计量序列，哪里显著、哪里不显著一目了然，定位能力比MK的曲线交叉法直观得多。

另外，MK检验对序列中的多个突变点不够敏感。如果序列先后经历了两次均值突变，MK的UF/UB曲线往往只会在第一次突变附近交叉，第二次突变容易被整体趋势掩盖。滑动t检验则不同，只要某个时间点前后两段均值差异足够大，就会被标记出来，更容易发现多重突变。

两者也各有前提。MK是非参数方法，不要求数据正态，对异常值稳健；滑动t检验是参数方法，要求数据近似正态、方差齐性、样本基本独立。所以严谨的做法其实是互补使用：先用MK判断序列有没有显著趋势变化，再用滑动t检验做突变定位，最后用Pettitt检验或B-G分割算法交叉验证突变年份。单独依赖任何一种方法都容易出问题。

1.3 适用场景和限制

滑动t检验最适合检测的对象是均值突变，也就是序列从一个平均水平跳变到另一个平均水平的情况。典型场景包括年平均气温序列的增暖突变、降水序列的阶段性变化、径流序列受人类活动影响后的均值跃变、植被指数序列的突变等。这类应用中，前后两段的方差通常相对稳定，滑动t检验的效果很好。

它的限制也很明确。第一，对渐变型趋势不友好：如果序列是持续上升或持续下降的趋势，任意前后两段的均值都会有差异，滑动t检验几乎处处显著，产生大量虚假突变点。第二，对短序列效果差：子序列长度n本身就占掉了2n个样本，序列总长度最好在30年以上，否则统计功效非常有限。第三，t检验的独立性假设在气候水文序列中经常被违反，需要做自相关处理，这一点后面会专门讲。

2. 动手前先想清楚：数据要求与滑窗长度

2.1 对输入数据的基本要求

写代码前先把数据准备这件事说清楚。滑动t检验对数据的基本要求是等间隔、连续、完整。所谓等间隔，指的是时间步长一致，比如逐年值、逐月值或者逐日值，不能中间有跳变。逐月数据直接拿来做突变检测会引入明显的季节周期影响，通常要先聚合成年值，或者去掉季节循环后再分析。

缺失值也是个容易踩的坑。Matlab的mean和var函数遇到NaN会直接返回NaN，如果序列中间有缺测，滑动t检验会在缺测位置及其附近全部产生NaN，画图时出现一大段空缺。所以数据进来第一件事就是检查缺测。常用的做法是用fillmissing做线性插值，或者用前后多年的气候平均值插补。对突变检测来说，插补会略微平滑序列，但只要缺测比例不太高，影响可以接受。

正态性问题也需要心里有数。t检验对正态偏离有一定的稳健性，尤其是样本量较大时。但如果你处理的是极端降水这类偏态很强的变量，最好先做对数变换或开方变换，让数据更接近正态，再跑滑动t检验。否则检验的显著性水平会有偏差。

2.2 滑窗长度n怎么选

滑窗长度n是整个方法里最关键的参数，没有之一。n选太小，比如n=2，两个子序列样本太少，检验功效很低，而且临界值因为自由度低而变得很大，很难检出突变；n选太大，比如n=20，前后两段各取20个点，窗口内部的均值信息被“抹平”了，突变前后的差异被稀释，突变点定位也变得模糊。

文献里最常见的取法是n=5或n=10。n=5对短时突变更敏感，但容易受局部波动干扰；n=10更稳健，但会漏掉一些持续时间短的突变。更稳妥的做法不是纠结选哪个，而是把多个n值的结果放在一起比较：分别取n=3、5、7、10跑一遍，如果不同窗口宽度下检出的突变年份基本一致，说明这个突变是真实存在的；如果结果随n变化漂移很大，那多半只是局部的随机波动，不能当作真正的突变。这个经验我在多个数据集上验证过，非常管用。

顺带提醒一下，n的最小合法值是2，因为n=1时自由度2*1-2=0，t分布根本没有定义。代码里最好加个判断，n小于2直接报错。

2.3 显著性水平和临界值的确定

显著性水平alpha一般取0.05或0.01。由于前后子序列等长，自由度恒为2n-2，临界值只由n和alpha决定，不随位置变化，画图时直接画一条水平线就行。下面这张表是几个常用取法的双侧临界值，方便没有统计工具箱的读者手动对照。

这里特别想提醒一点：很多人习惯直接用1.96这个来自正态分布的临界值，这在n较大的时候误差还能接受，但n=5时真实临界值是2.306，用1.96会凭空多出一堆“显著”的假突变点。滑动t检验是小样本t检验，临界值必须查t分布或调用tinv，不能拿正态近似偷懒。

3. 核心代码实现与逐行拆解

3.1 完整函数代码

下面这个函数是核心实现。我用一个独立的m文件来写，输入是原始序列、滑窗长度n和显著性水平alpha，输出是与原始序列等长的t统计量序列、临界值和自由度。所有能够提前确定的信息都放在函数里，调用侧只需关心业务逻辑。

function [t_stat, crit, df] = sliding_t_test(data, n, alpha) % SLIDING_T_TEST 滑动t检验突变检测 % 输入: % data - 一维数值序列，行向量或列向量均可 % n - 前后子序列长度，默认5 % alpha - 显著性水平，默认0.05 % 输出: % t_stat - 与data等长的t统计量序列，边界位置为NaN % crit - 双侧检验临界值（正值） % df - 自由度 if nargin < 2 || isempty(n) n = 5; end if nargin < 3 || isempty(alpha) alpha = 0.05; end if n < 2 error('滑窗长度n必须不小于2'); end data = data(:)'; % 统一成行向量 N = numel(data); t_stat = nan(1, N); % 边界处保留NaN，绘图时自动断开 for i = n+1 : N-n+1 seg1 = data(i-n : i-1); % 突变点之前的n个样本 seg2 = data(i : i+n-1); % 突变点之后的n个样本（含当前点i） n1 = numel(seg1); n2 = numel(seg2); mean1 = mean(seg1); mean2 = mean(seg2); var1 = var(seg1); var2 = var(seg2); sp2 = ((n1-1)*var1 + (n2-1)*var2) / (n1+n2-2); sp = sqrt(sp2); if sp == 0 % 两段都无波动，单独处理，避免除零 if mean1 == mean2 t_stat(i) = 0; else t_stat(i) = sign(mean1 - mean2) * Inf; end continue; end t_stat(i) = (mean1 - mean2) / (sp * sqrt(1/n1 + 1/n2)); end df = 2*n - 2; crit = tinv(1 - alpha/2, df); end

这个函数的主体逻辑其实只有不到二十行。每次循环计算两个子序列的均值、方差，然后套用标准t统计量公式。循环范围从n+1到N-n+1，这样既能保证seg1和seg2都恰好有n个样本，又不会越界。t值赋在位置i上，i就是突变点的候选位置，语义和后续的判读逻辑直接对应。

关于除零的处理，我单独写了一段。实际数据里确实会遇到这种情况，比如某段子序列恰好全是同一个数值（干旱区降水序列常年为0，或者数据被四舍五入后出现很多重复值），此时合并方差为0，直接的除法会得到NaN或者Inf。我这里的处理原则是：如果两段均值也相同，说明没有突变，t赋0；如果均值不同，那差异就是无穷显著，赋正负Inf。这样不至于程序崩溃，不过Inf在画图时会被Matlab当作异常值处理，后面主脚本里需要把这些非有限值替换掉，我下面会展示。

3.2 主脚本调用与绘图

函数写好后，主脚本非常清爽。下面这段代码做了三件事：调用滑动t检验函数、把非有限值处理掉、绘制t统计量序列和临界值参考线。

% 加载你的数据，这里用随机数据做演示 rng(42); N = 80; data = randn(1, N) + 5; data(41:end) = data(41:end) + 1.2; % 模拟第40~41年处的一次均值突变 year = 1961:2040; n = 5; alpha = 0.05; [t_stat, crit, df] = sliding_t_test(data, n, alpha); t_stat(~isfinite(t_stat)) = NaN; % 将Inf替换为NaN，避免画图拉伸 figure('Color', 'w', 'Position', [100 100 900 500]); plot(year, t_stat, 'b-', 'LineWidth', 1.5); hold on; yline(crit, 'r--', 'LineWidth', 1.2); yline(-crit, 'r--', 'LineWidth', 1.2); xlabel('年份'); ylabel('t统计量'); title(sprintf('滑动t检验突变检测 (n=%d, \\alpha=%.2f)', n, alpha)); legend('t统计量', sprintf('临界值 ±%.3f', crit), 'Location', 'best'); grid on; set(gca, 'FontSize', 12, 'Box', 'on');

这里用了yline画水平参考线，Matlab R2018b及以上版本都支持。如果你的版本比较老，可以退回到plot([year(1) year(end)],[crit crit],'r--')的写法。随机模拟数据在分段点处设置了一个1.2的均值跳变，跑出来的t统计量应该会在第40年附近形成一个明显的尖峰并超过临界线。

有几点提醒。第一，t_stat的边界处是NaN，plot画图时NaN会导致线段断开，这正是我们想要的视觉效果，不会在序列两端画出误导性的曲线。第二，Inf替换成NaN是为了防止整张图的纵坐标被拉伸到离谱的范围，这在真实数据中出现子序列方差为0时必须处理。第三，t_stat和year是同长度的，所以直接plot(year, t_stat)就能保证横坐标对齐，不需要额外裁剪。

4. 突变点判读要诀：检验统计量结果怎么解读

4.1 定位突变年份的两个判据

拿到t统计量序列后，怎么判定哪一年是突变年？我的经验是两个判据配合使用。

第一个判据是单点显著：某个年份i的|t_stat(i)| > crit，说明以该点为分界的前后两段均值差异显著，i是一个候选突变点。这是最基础的判据，但单独使用容易误判，因为随机波动造成个别点超过临界值的情况并不少见。

第二个判据是连续显著区段加峰值定位：如果从i1到i2连续若干年都满足|t| > crit，说明突变发生在这个时段内，具体年份取这个时段内|t|最大的位置。这个判据更稳健。我以前处理一个站点的年平均气温序列时，突变区段覆盖了连续五六年，一开始直接看单点挑了第一年，后来发现取峰值点年份才能和实际观测记录对应上。从那以后我一直用峰值定位，显著改善。

下面这段代码实现了连续显著区段筛选和峰值定位：

over = abs(t_stat) > crit; i = 1; while i <= numel(over) if over(i) j = i; while j <= numel(over) && over(j) j = j + 1; end seg_idx = i : j-1; [~, imax] = max(abs(t_stat(seg_idx))); fprintf('显著时段 %d-%d, 候选突变年: %d\n', ... year(seg_idx(1)), year(seg_idx(end)), year(seg_idx(imax))); i = j; else i = i + 1; end end

注意t统计量赋给的年份i，含义是“以i为分界点，前n个样本与后n个样本的均值差异”。读结果时，如果候选突变年是1978，意思是序列在1978年前后发生了均值跳变，而不是说1978年本身是异常的一年。这两个表述在论文里差很多，写的时候要拎清。

4.2 多窗口交叉验证判断突变是否真实

前面提到过，n的取值会影响结果。把交叉验证落实到操作层面，方法很简单：写一个循环，对多个n逐一调用滑动t检验，输出每个n检出的显著突变年份集合，然后求公共部分。

for n_try = [3, 5, 7, 10] [t_stat_try, crit_try] = sliding_t_test(data, n_try, alpha); over_try = abs(t_stat_try) > crit_try; if any(over_try) fprintf('n=%2d: 突变候选年 ', n_try); fprintf('%d ', year(over_try)); fprintf('\n'); else fprintf('n=%2d: 无显著突变\n', n_try); end end

我在实际项目中看结果有一套自己的习惯：n=3的结果参考价值有限，因为自由度过低、临界值过大，通常只有非常剧烈的突变才能超过；n=5和n=7的结果最常用；n=10的结果可以作为稳健性佐证。如果某个突变年份在n=5和n=7下同时出现，而且在n=10下虽然不完全一致但落在邻近两三年内，那基本可以认定这是个真实的突变。如果只有n=3检出来，其他窗口都没有，那基本是短时噪声在作祟，可以直接忽略。

这套多窗口验证的思路不仅适用于滑动t检验，其他突变检测方法也可以借鉴。突变检测本质上是一个从噪声里找信号的过程，单一参数设定下检出的信号可能是偶然，多个参数设定下稳定出现的信号才更可信。

4.3 多重比较问题与显著性水平收紧

滑动t检验在N个位置上执行了约N-2n次独立的t检验，这天然存在多重比较问题。即使序列完全不包含任何突变，按alpha=0.05的判据，每20个检验里就有一个会碰巧超过临界值。如果序列长度80、n=5，有效检验位置大约70个，随机假阳性的期望就有大约3.5个。这就是为什么很多序列跑出来总是能“检到”几个突变点的原因之一。

处理多重比较，主流思路有两种。一种是把显著性水平收紧到alpha=0.01，这样单次检验的假阳性率大幅下降。另一种是使用多重比较校正，比如Bonferroni校正——将alpha除以有效检验次数N-2n，得到更严格的单次检验阈值。Bonferroni虽然保守，但作为快速筛查还是够用的。实际项目中我通常的做法是：先用0.05水平筛选候选突变点，再用0.01水平或连续显著区段判据二次筛选，把两个判据都满足的位置作为最终突变点。不要机械地套用任何一个单一判据，结合人工判断是必须的。

5. 踩过才知道的坑：几个高频问题与调试经验

5.1 数据自相关导致的假突变

这是我踩过最深的坑。气候和水文序列几乎都有正自相关，也就是当年的值受前一年影响。t检验的基本假设之一是样本独立，当这个假设被违反时，实际自由度远小于名义自由度，检验的标准误被低估，结果会过于“激进”，大量伪突变被检出来。

有一个简化的修正办法：用有效样本量替代实际样本量。有效样本量的近似公式是：

Neff = N * (1 - r1) / (1 + r1)

其中r1是序列的一阶自相关系数。对滑动t检验来说，可以分别对seg1和seg2计算一阶自相关，用修正后的Neff1和Neff2替代n1和n2放进t统计量公式，自由度也相应改为Neff1+Neff2-2。Matlab里一阶自相关可以先用autocorr或直接把序列和它自身滞后一期的版本做相关系数。

需要说明的是，这个修正只是工程化的近似处理，学术上更严格的做法是采用块状Bootstrap等重采样方法，或者使用专门针对自相关时间序列的突变检测方法。不过作为日常数据分析的快速筛查，有效样本量修正已经能明显减少假突变，性价比很高。

5.2 渐变型突变带来的“平台期”问题

真实数据里的突变不一定都是阶梯式的，更常见的是在一个较短的时间段内连续变化，形成所谓的渐变型突变。这种情况下，t统计量不是单个尖峰，而是在连续几年内形成一个“平台”，整个平台期都超过临界线。定位突变年份时要取平台内|t|最大的点，而不是平台起点或终点。这个点对应两段均值差异最大的分割位置，在物理意义上最接近突变发生的中心时刻。

如果平台期特别宽，说明这个“突变”更接近一个缓慢的均值转移过程。此时要谨慎使用“突变”这个词，在论文里描述为“均值在X-X年期间发生显著变化”可能更严谨。我见过不少初学者把整个显著区段都标成突变年份，画在图上就是密密麻麻几根竖线，审稿人看到基本都会质疑，务必注意。

5.3 边界信息丢失与对齐问题

滑动t检验有一个天生短板：序列的前n个点和后n个点无法计算t统计量。一段80年的序列，n=5时前后各有5年没有t值，有效覆盖范围只有70年。这不是bug，而是方法本身的样本需求决定的。处理方式就是保持t_stat两端为NaN，绘图时断开，不要用插值去硬补。硬补会让边界处的“显著”假象看起来像真实信号。

另一个对齐方面的坑是索引漂移。t_stat(i)对应的是data中第i个位置，但很多人画图时直接把t_stat和年份向量plot在一起，如果年份向量不是从1开始，或者中间有跳年，就会错位。建议永远保持数据、年份、t_stat三者长度一致，并且用一个统一的索引号来对齐。我自己的习惯是代码里不写裸的数字，所有位置都用year对应的逻辑索引来取，很少出错。

5.4 子序列方差为零及其他边界异常

前面代码里已经处理了sp==0的情况，但在真实项目中，类似这种“非典型数据”造成的异常远比想象中多。比如某段数据全是整数，四舍五入后相同值大量重复，子序列方差可能小到接近机器精度；又比如序列中存在极端离群值，某一段均值被一个极大值拉高，造成t值爆表。这些异常会直接破坏检验结果，甚至让后续的多窗口对比完全失效。

我的建议是：在跑滑动t检验之前，先对数据做一个基本体检。用summary或直接画时间序列图，检查有没有明显的离群值、有没有长时间恒定不变的区段、有没有整体趋势。有离群值先决定是保留、剔除还是做稳健变换；有趋势先考虑是否要做去趋势；有明显突变但不知道位置的先目测一个大致范围。数据体检看起来是笨功夫，实际能省下后面排查异常结果的大量时间。

还有一个小细节，很多人会在代码里直接用ttest2函数来做滑动t检验。ttest2默认使用的是Welch检验（不假定方差齐性），而传统滑动t检验用的是合并方差的Student版本，两者的统计量和自由度不同，不能直接互换。如果确实想用ttest2，需要设置'Vartype','equal'来指定方差齐性。我给出的函数是手写版，好处是你完全清楚每一步在算什么，后续想改成Welch版本、加入自相关修正，都是几分钟的事。

从实际操作来看，滑动t检验最让人放心的地方就是可解释性强。跑完结果，哪年突变、方向如何、强度多大，全部能从t统计量序列里读出来。配合多窗口交叉验证和显著性水平收紧，它在突变检测工具箱里是性价比很高的一个成员。最后分享一个我自己的小习惯：每次输出结果时，除了候选突变年份，我会把该年份前后两段的均值差也一并打出来。这样报告里写“该站点年均气温在1997年前后发生显著上升突变，增幅约0.8度”时，每一个数字都有据可查，而不是只有一张显著性图。