做气候或水文序列突变检测的十有八九开口就是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)) / (n1n2-2))这个统计量服从自由度为n1n2-2的t分布。由于滑动t检验在大多数应用中都会取对称窗口也就是n1n2n所以自由度就是2n-2。给定显著性水平alpha之后临界值可以从t分布表查也可以用Matlab的tinv函数一行算出来。需要注意一个容易搞反的符号问题t统计量的分子是mean(seg1)-mean(seg2)所以t0说明前段均值更高意味着序列在i点附近发生了均值下降突变t0说明后段均值更高发生了均值上升突变。很多人跑完结果只看绝对值忽略了符号方向导致对突变方向的描述写反这是写论文时很常见的低级错误。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选太小比如n2两个子序列样本太少检验功效很低而且临界值因为自由度低而变得很大很难检出突变n选太大比如n20前后两段各取20个点窗口内部的均值信息被“抹平”了突变前后的差异被稀释突变点定位也变得模糊。文献里最常见的取法是n5或n10。n5对短时突变更敏感但容易受局部波动干扰n10更稳健但会漏掉一些持续时间短的突变。更稳妥的做法不是纠结选哪个而是把多个n值的结果放在一起比较分别取n3、5、7、10跑一遍如果不同窗口宽度下检出的突变年份基本一致说明这个突变是真实存在的如果结果随n变化漂移很大那多半只是局部的随机波动不能当作真正的突变。这个经验我在多个数据集上验证过非常管用。顺带提醒一下n的最小合法值是2因为n1时自由度2*1-20t分布根本没有定义。代码里最好加个判断n小于2直接报错。2.3 显著性水平和临界值的确定显著性水平alpha一般取0.05或0.01。由于前后子序列等长自由度恒为2n-2临界值只由n和alpha决定不随位置变化画图时直接画一条水平线就行。下面这张表是几个常用取法的双侧临界值方便没有统计工具箱的读者手动对照。滑窗长度n自由度df2n-2alpha0.05临界值alpha0.01临界值342.7764.604582.3063.3557122.1793.05510182.1012.878这里特别想提醒一点很多人习惯直接用1.96这个来自正态分布的临界值这在n较大的时候误差还能接受但n5时真实临界值是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 n1 : N-n1 seg1 data(i-n : i-1); % 突变点之前的n个样本 seg2 data(i : in-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) / (n1n2-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统计量公式。循环范围从n1到N-n1这样既能保证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的边界处是NaNplot画图时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我在实际项目中看结果有一套自己的习惯n3的结果参考价值有限因为自由度过低、临界值过大通常只有非常剧烈的突变才能超过n5和n7的结果最常用n10的结果可以作为稳健性佐证。如果某个突变年份在n5和n7下同时出现而且在n10下虽然不完全一致但落在邻近两三年内那基本可以认定这是个真实的突变。如果只有n3检出来其他窗口都没有那基本是短时噪声在作祟可以直接忽略。这套多窗口验证的思路不仅适用于滑动t检验其他突变检测方法也可以借鉴。突变检测本质上是一个从噪声里找信号的过程单一参数设定下检出的信号可能是偶然多个参数设定下稳定出现的信号才更可信。4.3 多重比较问题与显著性水平收紧滑动t检验在N个位置上执行了约N-2n次独立的t检验这天然存在多重比较问题。即使序列完全不包含任何突变按alpha0.05的判据每20个检验里就有一个会碰巧超过临界值。如果序列长度80、n5有效检验位置大约70个随机假阳性的期望就有大约3.5个。这就是为什么很多序列跑出来总是能“检到”几个突变点的原因之一。处理多重比较主流思路有两种。一种是把显著性水平收紧到alpha0.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统计量公式自由度也相应改为Neff1Neff2-2。Matlab里一阶自相关可以先用autocorr或直接把序列和它自身滞后一期的版本做相关系数。需要说明的是这个修正只是工程化的近似处理学术上更严格的做法是采用块状Bootstrap等重采样方法或者使用专门针对自相关时间序列的突变检测方法。不过作为日常数据分析的快速筛查有效样本量修正已经能明显减少假突变性价比很高。5.2 渐变型突变带来的“平台期”问题真实数据里的突变不一定都是阶梯式的更常见的是在一个较短的时间段内连续变化形成所谓的渐变型突变。这种情况下t统计量不是单个尖峰而是在连续几年内形成一个“平台”整个平台期都超过临界线。定位突变年份时要取平台内|t|最大的点而不是平台起点或终点。这个点对应两段均值差异最大的分割位置在物理意义上最接近突变发生的中心时刻。如果平台期特别宽说明这个“突变”更接近一个缓慢的均值转移过程。此时要谨慎使用“突变”这个词在论文里描述为“均值在X-X年期间发生显著变化”可能更严谨。我见过不少初学者把整个显著区段都标成突变年份画在图上就是密密麻麻几根竖线审稿人看到基本都会质疑务必注意。5.3 边界信息丢失与对齐问题滑动t检验有一个天生短板序列的前n个点和后n个点无法计算t统计量。一段80年的序列n5时前后各有5年没有t值有效覆盖范围只有70年。这不是bug而是方法本身的样本需求决定的。处理方式就是保持t_stat两端为NaN绘图时断开不要用插值去硬补。硬补会让边界处的“显著”假象看起来像真实信号。另一个对齐方面的坑是索引漂移。t_stat(i)对应的是data中第i个位置但很多人画图时直接把t_stat和年份向量plot在一起如果年份向量不是从1开始或者中间有跳年就会错位。建议永远保持数据、年份、t_stat三者长度一致并且用一个统一的索引号来对齐。我自己的习惯是代码里不写裸的数字所有位置都用year对应的逻辑索引来取很少出错。5.4 子序列方差为零及其他边界异常前面代码里已经处理了sp0的情况但在真实项目中类似这种“非典型数据”造成的异常远比想象中多。比如某段数据全是整数四舍五入后相同值大量重复子序列方差可能小到接近机器精度又比如序列中存在极端离群值某一段均值被一个极大值拉高造成t值爆表。这些异常会直接破坏检验结果甚至让后续的多窗口对比完全失效。我的建议是在跑滑动t检验之前先对数据做一个基本体检。用summary或直接画时间序列图检查有没有明显的离群值、有没有长时间恒定不变的区段、有没有整体趋势。有离群值先决定是保留、剔除还是做稳健变换有趋势先考虑是否要做去趋势有明显突变但不知道位置的先目测一个大致范围。数据体检看起来是笨功夫实际能省下后面排查异常结果的大量时间。还有一个小细节很多人会在代码里直接用ttest2函数来做滑动t检验。ttest2默认使用的是Welch检验不假定方差齐性而传统滑动t检验用的是合并方差的Student版本两者的统计量和自由度不同不能直接互换。如果确实想用ttest2需要设置Vartype,equal来指定方差齐性。我给出的函数是手写版好处是你完全清楚每一步在算什么后续想改成Welch版本、加入自相关修正都是几分钟的事。从实际操作来看滑动t检验最让人放心的地方就是可解释性强。跑完结果哪年突变、方向如何、强度多大全部能从t统计量序列里读出来。配合多窗口交叉验证和显著性水平收紧它在突变检测工具箱里是性价比很高的一个成员。最后分享一个我自己的小习惯每次输出结果时除了候选突变年份我会把该年份前后两段的均值差也一并打出来。这样报告里写“该站点年均气温在1997年前后发生显著上升突变增幅约0.8度”时每一个数字都有据可查而不是只有一张显著性图。