MATLAB长时间序列降水MK趋势分析：从数据预处理到突变点检测的完整指南

简介这份资源面向气候、水文及地理等专业的学生与科研初学者提供基于MATLAB实现长时间序列Mann-Kendall非参数趋势检验与突变检验的完整实验资料可用于干旱趋势分析等场景。压缩包内共1个doc文档约457KB系统梳理了MK趋势分析法的秩计算、符号判定、Z-score构造与显著性判断流程以及突变检验中秩序列、UF与UB统计量的计算和突变点识别方法并配有相应公式推导与程序说明。文档以1961—2018年三江源区40个站点的SPI3数据为实验对象展示了Z值矩阵输出与可视化结果得出该区域SPI3、SPI12整体呈上升趋势、约80%站点通过0.1显著性检验的结论说明干旱化趋势有所下降。目前已有6383人学习下载适合希望快速上手MK检验、理解趋势与突变分析思路的读者参考借鉴。1. 长时间序列降水做 MK 趋势分析为什么有人跑出趋势有人跑出噪声手里攒了三十年以上的月降水或年降水序列想判断它到底有没有显著上升或下降趋势Mann-KendallMK检验几乎是绕不开的一步。它不要求数据服从正态分布对缺失值和异常值也没那么敏感这是它比线性回归斜率更受欢迎的根本原因。但我在实际帮人看结果时发现一个很普遍的现象同一份降水数据有人跑出 p0.05 的显著上升有人跑出 p0.3 的“无趋势”差别往往不在数据本身而在预处理和参数选择上。这篇笔记就围绕“用 MATLAB 对长时间序列降水做 MK 趋势分析”这件事把从数据整理、检验计算、突变点识别到结果解读的完整链路拆开讲。适合手上有气象站、格点降水或再分析降水数据需要做趋势显著性判断的从业者和研究生。读完你应该能自己写出一套可复现的脚本并且知道哪些地方最容易翻车。2. MK 检验的数学底子与 MATLAB 落地前的三个决策MK 检验的核心是统计量 S对序列中每一对 (i, j)ij比较 x_j 与 x_i 的大小大于记 1小于记 -1相等记 0全部累加得到 S。S 的期望为 0方差在无 tied 值时是 n(n-1)(2n5)/18。标准化后得到 Z双侧检验下 |Z|1.96 对应 p0.05。原理不复杂但落到代码前有三个决策必须先定否则后面全是返工。2.1 用原始 MK 还是去趋势的 TFPW-MK原始 MK 有一个被反复讨论的假设序列是随机独立的。降水序列往往存在自相关尤其是月尺度数据今年 7 月和去年 7 月、前年 7 月之间可能有相关性。自相关会让方差被低估Z 值被高估本来不显著的趋势被误判为显著。常见的补救是预白化pre-whitening其中 TFPWTrend-Free Pre-Whitening用得最多先估斜率去掉趋势对残差做 AR(1) 预白化再把趋势加回去最后做 MK。我的建议是年降水序列自相关通常较弱可以直接用原始 MK月降水序列强烈建议跑一遍 TFPW-MK 做对照。如果两者结论一致可以放心报原始 MK如果原始显著、TFPW 不显著就要警惕自相关导致的假显著。2.2 显著性水平与双侧/单侧怎么选气象里默认双侧检验α0.05对应 |Z|1.96。但如果你有明确的物理预期比如“某区域暖季降水应该增加”可以用单侧阈值变成 1.645。问题在于很多人写论文时默认双侧却在解读时按单侧说“接近显著”这是不严谨的。我一般会在脚本里同时输出 Z、p 和置信区间让读者自己判断。2.3 缺失值、平局值和序列长度MK 对缺失值没有原生处理必须先把缺失值补上或剔除。剔除会破坏等间隔假设影响突变点检测插补又会引入人为趋势。常见做法是缺失比例低于 5% 用线性插值高于 5% 就考虑换站点或换时段。平局值tied values在降水里很常见比如很多天降水为 0年尺度上少见月尺度上如果某月常年为 0 就会出现。MATLAB 的 ktaub 或自己写的 MK 要正确处理 tied 校正项否则方差偏大、Z 偏小。序列长度方面MK 在 n10 时检验力很低n20 才比较稳。三十年逐年数据 n30 是够的但如果做季节尺度每季只有 30 个点也要注意。3. 从原始降水表到 MK 结果一套可复现的 MATLAB 脚本这一章给出一套我常用的脚本骨架输入是一列年份加一列降水值输出是 S、Z、p、Sen 斜率以及突变点。你可以直接改成自己的数据格式。3.1 数据读入与预处理假设数据存为 precip.csv第一列年份第二列年降水。先读入、检查缺失、做基本清洗。% 读入降水数据假设两列年份、降水 data readmatrix(precip.csv); year data(:,1); precip data(:,2); % 检查缺失值 if any(isnan(precip)) fprintf(存在缺失值数量%d\n, sum(isnan(precip))); % 线性插值补齐仅建议缺失比例低时使用 precip fillmissing(precip, linear); end % 检查平局值比例 tie_ratio 1 - numel(unique(precip)) / numel(precip); fprintf(平局值比例%.2f%%\n, tie_ratio * 100); % 基本统计 fprintf(序列长度%d均值%.2f标准差%.2f\n, ... numel(precip), mean(precip), std(precip));这段代码做了三件事缺失值插补、平局值比例检查、基本统计输出。平局值比例超过 20% 就要在 MK 里做 tied 校正否则 Z 会偏小。插补只在缺失少时用缺失多的话建议直接换数据。3.2 原始 MK 检验的实现MATLAB 没有内置 MK需要自己写。下面这个函数返回 S、Z、p。function [S, Z, p] mk_test(x) n numel(x); S 0; for i 1:n-1 for j i1:n S S sign(x(j) - x(i)); end end % 计算方差含 tied 校正 [~, ~, ic] unique(x); tie_counts accumarray(ic, 1); tie_term sum(tie_counts .* (tie_counts - 1) .* (2 * tie_counts 5)); varS (n * (n - 1) * (2 * n 5) - tie_term) / 18; % 标准化 if S 0 Z (S - 1) / sqrt(varS); elseif S 0 Z (S 1) / sqrt(varS); else Z 0; end % 双侧 p 值 p 2 * (1 - normcdf(abs(Z))); end关键点在方差计算tie_term 是平局校正项没有它方差会偏大。Z 的连续性校正S0 减 1S0 加 1在小样本时影响明显。p 用正态近似n10 时够用。3.3 Sen 斜率与 TFPW 预白化MK 只给显著性不给趋势大小。Sen 斜率是中位数斜率抗异常值。TFPW 则是在 MK 前做预白化。% Sen 斜率 function slope sen_slope(x) n numel(x); slopes []; for i 1:n-1 for j i1:n slopes(end1) (x(j) - x(i)) / (j - i); end end slope median(slopes); end % TFPW-MK function [Z, p, slope] tfpw_mk(x) n numel(x); t (1:n); % 估斜率 b sen_slope(x); % 去趋势 detrended x - b * t; % AR(1) 预白化 r corr(detrended(1:end-1), detrended(2:end)); y detrended(2:end) - r * detrended(1:end-1); % 加回趋势 y y b * t(2:end); [~, Z, p] mk_test(y); slope b; endSen 斜率直接给出“每年变化多少毫米”比回归斜率更稳。TFPW 里 AR(1) 系数 r 用相邻点相关估计去趋势后再白化最后加回趋势。注意 TFPW 后序列长度少 1这是正常的。3.4 突变点检测Mann-Kendall 突变检验趋势之外很多人还关心“哪一年发生了突变”。常用的是顺序 MK 突变检验构造正序和逆序统计量 UF、UB交点即突变点。function [UF, UB] mk_mutation(x) n numel(x); UF zeros(n,1); UB zeros(n,1); % 正序 for k 2:n s 0; for i 1:k-1 s s sign(x(k) - x(i)); end E k * (k - 1) / 4; Var k * (k - 1) * (2 * k 5) / 72; if s 0 UF(k) (s - 1) / sqrt(Var); elseif s 0 UF(k) (s 1) / sqrt(Var); end end % 逆序 xr flipud(x); for k 2:n s 0; for i 1:k-1 s s sign(xr(k) - xr(i)); end E k * (k - 1) / 4; Var k * (k - 1) * (2 * k 5) / 72; if s 0 UB(n - k 1) -(s - 1) / sqrt(Var); elseif s 0 UB(n - k 1) -(s 1) / sqrt(Var); end end endUF 和 UB 的交点如果在置信区间内就是显著突变点。注意这里没有做 tied 校正年降水平局少影响不大月降水要补上。3.5 结果汇总与出图把上面拼起来输出一张表加一张图。% 主流程 [S, Z, p] mk_test(precip); slope sen_slope(precip); [Z_tfpw, p_tfpw, ~] tfpw_mk(precip); [UF, UB] mk_mutation(precip); fprintf(原始 MK: S%d, Z%.3f, p%.4f\n, S, Z, p); fprintf(TFPW-MK: Z%.3f, p%.4f\n, Z_tfpw, p_tfpw); fprintf(Sen 斜率: %.3f mm/year\n, slope); figure; subplot(2,1,1); plot(year, precip, b-o); hold on; plot(year, polyval(polyfit(year, precip, 1), year), r--); ylabel(降水 (mm)); title(降水序列与线性趋势); subplot(2,1,2); plot(year, UF, b-); hold on; plot(year, UB, r-); yline(1.96, k--); yline(-1.96, k--); ylabel(统计量); title(MK 突变检验);出图时注意UF 和 UB 的交点要肉眼确认脚本只能给候选。线性趋势线只是辅助真正报的是 Sen 斜率。4. 参数怎么设、结果怎么读五个容易翻车的细节这一章集中讲踩坑。每一条都是我在实际数据上遇到过的按“现象→原因→解决”写。4.1 现象Z 值很大但 p 不显著原因通常是序列太短。n15 时即使 Z2.0p 也可能到 0.06。MK 的检验力随 n 增长很慢n20 时不要过度解读 p 值。解决要么延长序列要么用 bootstrap 做经验 p 值要么明确说明样本量限制。4.2 现象月降水 MK 显著年降水不显著月尺度自相关强Z 被高估。解决跑 TFPW-MK 对照或者先把月降水聚合到年再检验。如果研究问题就是月尺度那要报 TFPW 结果并说明自相关处理。4.3 现象突变点每年都变原因是没有做逆序检验或者 UF、UB 交点落在置信区间外。解决确认 UB 的计算是否正确交点必须在 ±1.96 之间才算显著突变。另外突变点对序列长度敏感前 5 年和后 5 年不建议作为突变点。4.4 现象Sen 斜率和线性回归斜率符号相反原因通常是异常值。Sen 是中位数斜率抗异常值线性回归受极端值影响大。解决检查数据里的极端年份确认是否记录错误。如果数据没问题以 Sen 为准并在文中说明。4.5 现象平局值多导致 Z 偏小原因是没有做 tied 校正。解决在方差计算里加上 tie_term或者用 ktaub 这类已经处理好的函数。平局值比例超过 20% 时校正前后 Z 可能差 0.3 以上足以改变显著性结论。5. 把 MK 用稳的几个进阶习惯最后一章说几个我自己的习惯都是被数据教出来的。第一永远先画图再看统计量。降水序列的图能告诉你很多 MK 不会说的东西有没有台阶、有没有极端年、有没有明显的非线性。我见过一份数据 MK 不显著但图上 1990 年后明显抬升后来发现是站点迁移导致的不是气候信号。第二MK 只回答“有没有单调趋势”不回答“趋势是不是线性”。如果序列先降后升MK 可能给出不显著但实际有转折。这时候要配合突变检验和分段趋势。第三报结果时给全S、Z、p、Sen 斜率、置信区间、序列长度、缺失处理方式。只报一个 p 值是不负责任的。我一般会在表里列这几项审稿人一看就知道你做过什么。第四脚本要可复现。把数据读入、预处理、检验、出图写成一条流水线换一份数据只改文件名。我习惯把参数放在脚本开头比如 alpha0.05、use_tfpwtrue这样别人能一眼看到你做了什么选择。第五别迷信 p0.05。在长时间序列里p0.049 和 p0.051 没有本质区别。趋势的物理意义、量级和持续性比一个阈值更重要。我现在的做法是同时报 p 和 Sen 斜率让读者判断这个趋势值不值得关注。希望帮到你。本文还有配套的精品资源点击获取