从拿到一份几十年的降水资料开始，大家问得最多的就是两个问题：这雨量到底是增了还是减了？变化是不是有规律、存在几年一轮的周期？光看原始曲线，谁也说不清。解决这个问题，水文气象领域用得最顺手的组合就是Mann-Kendall（MK）检验加上Morlet小波分析。前者告诉你趋势和突变点，后者告诉你显著周期，两个方法叠起来看，才敢把"规律"两个字写进报告里。这篇文章我直接把MK检验和Morlet小波分析的原理、Matlab代码、参数选择和判读坑位一次讲透，适合正在写课程论文、毕业论文，或要处理实测降水资料的同行参考。
1. 为什么雨量分析总要叠着MK检验和小波一起看
1.1 给多年的雨量序列做体检：趋势、突变、周期各查什么
降水序列分析其实是三种"体检项目"的组合体。第一个项目是趋势，回答"这60年总雨量在长期增加还是减少"；第二个项目是突变，回答"有没有哪一年发生了明显转折，从偏多变成偏少"；第三个项目是周期，回答"是否存在3年、7年或者20年这种准周期性振荡"。
这三种问题对应的方法很不一样。趋势可以用线性回归做，但降水数据往往不满足正态性假设，而且个别极端年份（暴雨、特旱）会把回归斜率拉得很偏。MK检验是秩次检验，不要求数据服从正态分布，也不怕个别异常值，所以成为水文气象趋势判断的事实标准。突变点需要用序贯MK（UF/UB曲线）去捕捉，因为突变本质上是"某一段序列的统计特性发生结构性改变"。周期则要依赖小波分析，因为普通的周期图或功率谱只能告诉你"哪些频率能量高"，却丢了时间信息；而降水振荡往往不是平稳的，这几年有准3年周期，过几年又消失了，这种非平稳周期特征只有时频分析能刻画。
所以一套完整的雨量分析，实际是"趋势－突变－周期"三层体检，MK检验管前两层，Morlet小波管第三层。我把这套流程在多个站点的年降水序列上跑过，结论是：两个方法单独用都有盲区，合起来才覆盖得全面。
1.2 非参数方法在降水序列里为什么靠谱
MK检验属于非参数检验，意思是它不关心总体服从什么分布。它只比较序列中每一对观测值的大小关系：后面的值比前面的值大，记一个正号；比前面的值小，记一个负号。把所有正负号累加，得到一个统计量S。如果S明显偏离0，说明上升或下降的配对很多，趋势就显著；如果S在0附近抖动，说明序列是随机游走的，谈不上趋势。
这个思路在降水分析里非常好用。年降水量经常出现偏态分布，丰水年和平水年差异巨大，直接用基于均值和方差的方法容易被极端值绑架。MK检验只关心"谁大谁小"，遇到百年一遇的极端值，它在配对比较中也只贡献一个符号，不会被数值大小左右。后来配一个Sen斜率，还能估算趋势的幅度（比如每年减少多少毫米），算是把稳健性做到位了。观测序列里偶尔缺测、存在异常值的站点，MK检验也是相对稳妥的选择。
2. MK检验的统计逻辑与Matlab手写实现
2.1 S统计量、方差与标准化Z值的三步推演
MK检验的计算并不复杂，核心就三步。
第一步，计算配对符号和S。设序列有n个观测值，对所有i<j的组合计算sgn(xj−xi)，累加得到S。公式可以简写为：
$$S = \sum_{i=1}^{n-1}\sum_{j=i+1}^{n} \text{sgn}(x_j - x_i)$$
这个S的含义很直观：上升配对多时S为正，下降配对多时S为负。如果序列完全随机，S的期望接近0。
第二步，计算S的方差。当数据不存在重复值时，方差是Var(S)=n(n−1)(2n+5)/18。存在重复值（tie）时，要考虑结组的修正，减去Σtp(tp−1)(2tp+5)/18，其中tp是每组重复值的个数。降水数据非常容易重复，比如连续几年都是同一降雨量整数，所以修正项不能省。
第三步，标准化。用连续性修正因子把S转成标准正态Z值：
$$Z = \begin{cases} (S-1)/\sqrt{\text{Var}(S)} & S>0 \ 0 & S=0 \ (S+1)/\sqrt{\text{Var}(S)} & S<0 \end{cases}$$
Z服从标准正态分布，双边检验下，若|Z|>1.96（α=0.05），就拒绝无趋势的原假设。正Z表示上升趋势，负Z表示下降趋势。很多论文只写"MK检验通过0.05显著性检验"却不说Z的符号，其实Z的符号就是判断增减的关键。
2.2 趋势检验和Sen斜率代码实现
下面这段代码是我在多个项目里沿用的版本，分两部分：第一部分计算Z值、p值和结论；第二部分计算Sen斜率，即所有配对斜率的中位数。Sen斜率对异常值不敏感，可以和MK检验的Z值配套输出。
function [Z, p_value, H, Sen] = mk_trend_test(x, alpha) % MK趋势检验，x为列向量或行向量，alpha取0.05或0.01 % 输出：Z统计量，p值，H(1通过显著性检验)，Sen斜率(毫米/年) if nargin < 2, alpha = 0.05; end x = x(:); n = length(x); if n < 8 error('序列太短，MK检验通常要求n>=8'); end % 计算S S = 0; for i = 1:n-1 for j = i+1:n S = S + sign(x(j) - x(i)); end end % 结组修正 unique_x = unique(x); tp = []; for k = 1:length(unique_x) cnt = sum(x == unique_x(k)); if cnt > 1 tp(end+1) = cnt; end end varS = n*(n-1)*(2*n+5)/18; if ~isempty(tp) varS = varS - sum(tp.*(tp-1).*(2*tp+5))/18; end % 标准化 if S > 0 Z = (S-1)/sqrt(varS); elseif S < 0 Z = (S+1)/sqrt(varS); else Z = 0; end % p值与结论 p_value = 2 * (1 - normcdf(abs(Z))); H = p_value < alpha; % Sen斜率 d = []; for i = 1:n-1 for j = i+1:n d = [d; (x(j)-x(i))/(j-i)]; end end Sen = median(d); end实际调用时，我会先把年降水序列读进变量rain，然后一行拿到所有结果：
[Z, p, H, Sen] = mk_trend_test(rain, 0.05); fprintf('Z = %.2f, p = %.3f, Sen斜率 = %.2f mm/年\n', Z, p, Sen);如果H等于1，说明趋势显著。Sen斜率的符号和Z一致，数值上可以直接说"该站年降水量以每年X毫米的速率变化"。
2.3 序贯MK突变检验：UF、UB曲线的来龙去脉
趋势检验回答"有没有趋势"，突变检验回答"趋势在哪一年变化"。前者用整个序列算一个S，后者把序列从第2个点开始逐渐扩展，每扩展一步计算一次标准化统计量，得到UF曲线。再把序列倒过来做一遍同样的计算，得到UB曲线，绘制时按时间轴反向排列。两条曲线如果出现交点，且交点位于±1.96的置信区间内，这个时间就是候选突变点。
需要特别注意的是，这里不能用前面那个现成的统计量S公式直接每步重算，正确做法是对每个长度k重新统计内部配对符号。我推荐写成单独的函数：
function [UF, UB] = mk_abrupt(x) % 序贯MK突变检验 % 返回UF、UB两条曲线的统计量，序列长度n对应时间1~n x = x(:); n = length(x); UF = zeros(n,1); UB = zeros(n,1); % 正序计算UF for k = 2:n Sk = 0; for i = 1:k-1 for j = i+1:k Sk = Sk + sign(x(j) - x(i)); end end E = k*(k-1)/4; Var = k*(k-1)*(2*k+5)/72; UF(k) = (Sk - E)/sqrt(Var); end % 倒序计算UB并反转到正常时间轴 y = flipud(x); for k = 2:n Sk = 0; for i = 1:k-1 for j = i+1:k Sk = Sk + sign(y(j) - y(i)); end end E = k*(k-1)/4; Var = k*(k-1)*(2*k+5)/72; UB(n-k+1) = - (Sk - E)/sqrt(Var); end end绘图时把UF和UB画在同一坐标系，上下各画一条±1.96的参考线。交点年份落在两条参考线之间，才认为可能存在统计显著的突变点。我实践里的习惯是：只把交点年份当作"候选突变点"，还要结合前后均值的差、小波能量的变化一起确认，绝不能看画图上有个交叉就说突变。数据里缺测较多时，UF/UB会在中间段剧烈振荡，出现多个假交点，这种时候要优先处理缺测。
3. Morlet小波的参数细节与Matlab实现
3.1 从傅里叶到连续小波：天气序列为什么需要"时间－尺度"窗口
傅里叶变换告诉我们序列里有"哪些频率成分"，但丢掉了"这些成分出现在什么年代"。降水振荡的特点是时变：某个准周期只在特定年代明显。为了同时获得频率和时间信息，就要用小波变换。
连续小波变换的原始定义是让一个母小波函数在不同尺度a上伸缩、在不同时间b上平移，和原始序列做内积：
$$W(a,b) = \frac{1}{\sqrt{a}} \int x(t), \psi^*\left(\frac{t-b}{a}\right) dt$$
这里a控制母小波被拉伸或压缩的程度，也就是控制分析窗口的时间宽度和频率中心；a小，窗口短，能抓住高频的短周期信号；a大，窗口长，适合捕捉低频的长周期。b则控制窗口在时间轴上的位置。对每个(a,b)组合求一个系数，最后得到一张"时间－尺度"二维谱图。
Morlet母小波是水文气象里最常用的复小波，形如一个复指数乘以高斯包络：
$$\psi(t) = \pi^{-1/4} e^{i\omega_0 t} e^{-t^2/2}$$
用复小波的好处是小波系数同时包含振幅和相位。振幅平方可以得到能量谱，实部可以画等值线来分析振荡的位相结构。降水序列本身是实信号，复小波能把变化是"偏多"还是"偏少"以及相对于时间轴的位相信息都保留下来。
3.2 尺度数组、周期换算与omega0的取值逻辑
自写Morlet小波分析时，最关键的两个参数是中心频率ω0和尺度数组。
ω0决定小波的时间分辨率与频率分辨率之间的权衡。ω0越大，高斯包络在频域越窄，频率分辨率越高，但时域窗口越宽，时间定位越模糊。水文气象里几乎约定俗成取ω0=6，因为在ω0=6附近，小波的时频分辨率和经典Heisenberg不确定性折中得比较好，而且此时尺度a与傅里叶周期T有非常简洁的近似关系：
$$T \approx \frac{4\pi}{\omega_0 + \sqrt{2+\omega_0^2}}, a \approx 1.033a$$
也就是说，当ω0=6时，小波尺度数值上几乎等于周期。很多新手在等值线图的纵轴直接标"尺度"而不换算成"周期（年）"，这在报告里是不合规的。论文常见表述是"准3年、准10年周期"，对应尺度换算后的周期值，不是原始尺度。
尺度数组的构建方式也影响计算量和谱图分辨率。我常用的是以2为底的对数均匀尺度：
s0 = 2; % 最小尺度 dj = 0.125; % 子尺度间隔，越小谱图越细腻 J = floor(log2(60/2)/dj); % 60是序列长度（年） scales = s0 * 2.^((0:J)*dj);这样得到一组约30~48个尺度，覆盖从2年到序列长度量级的范围。对于60年的序列，已经足够识别准3年、准5年、准15年这类水文气象常见周期。dj取0.125时尺度间隔很细，画出来的等值线图比较平滑；如果嫌计算慢，可放宽到0.25。
3.3 一份不依赖工具箱的Morlet CWT代码
网上有很多调用Wavelet Toolbox的cwt接口的做法，但我更喜欢自己写一个几十行的函数，原因是自写代码完全可控：尺度数组自定、边界处理自定、输出格式自定，不会有工具箱版本兼容的烦恼。下面这个函数直接用定义式做数值积分，思路清晰，对60~100年的序列计算量完全可以接受。
function [W, period, scale, coi] = morlet_cwt_self(x, dt, omega0) % 自写Morlet连续小波变换（定义式数值积分版） % 输入： % x : 1×n或n×1序列，建议先去均值 % dt : 采样间隔，年序列取1 % omega0 : Morlet中心频率，默认6 % 输出： % W : 小波系数矩阵，size为(尺度数×n) % period : 每个尺度对应的傅里叶周期 % scale : 尺度数组 % coi : 影响锥半径，用于判断边界可靠区 if nargin < 2, dt = 1; end if nargin < 3, omega0 = 6; end x = x(:)'; n = length(x); x = x - mean(x); t = (0:n-1)*dt; % 构建尺度数组 s0 = 2*dt; dj = 0.125; J = floor(log2(n*dt/s0)/dj); scale = s0 * 2.^((0:J)*dj); period = 4*pi/(omega0 + sqrt(2+omega0^2)) * scale; W = zeros(length(scale), n); for i = 1:length(scale) a = scale(i); for b = 1:n tb = t(b); % 母小波在时间层上的采样 psi = pi^(-1/4) * exp(1i*omega0*((t-tb)/a)) .* exp(-0.5*((t-tb)/a).^2); W(i,b) = sqrt(1/a) * trapz(t, x .* conj(psi)) * dt; end end % 影响锥的近似半径：用e-folding时间估计 coi = sqrt(2) * scale; end这段代码用于年尺度的稳态分析足够了。如果你处理的是逐月数据，要先把月序列减去多年逐月气候态，或者直接转成年序列再分析，不然季节信号会在高频区形成非常醒目但毫无意义的伪周期。
拿到W矩阵后，下一步是算小波能量谱和小波方差。小波方差是对每个尺度沿着时间维度累加能量，峰的位置对应主要周期。主程序如下：
% 假设rain为1961-2020年逐年降水量 [W, period, scale, coi] = morlet_cwt_self(rain, 1, 6); Power = abs(W).^2; waveVar = sum(Power, 2); % 小波方差 [~, idx] = max(waveVar); % 主周期索引 fprintf('小波方差最大峰对应周期: %.2f 年\n', period(idx)); % 绘制实部等值线图 figure; subplot(3,1,1); contourf(1961:2020, period, real(W), 30, 'linecolor','none'); set(gca, 'yscale', 'log', 'ytick', [2 4 8 16 32]); ylim([2 40]); xlabel('年份'); ylabel('周期/年'); colorbar; % 绘制小波方差图 subplot(3,1,2); plot(period, waveVar, 'linewidth', 1.5); set(gca, 'xscale', 'log', 'xtick', [2 4 8 16 32]); xlim([2 40]); xlabel('周期/年'); ylabel('小波方差'); grid on;等值线图建议用对数纵轴，否则短周期会被压缩成一条线，长周期的细部结构看不见。小波方差图也一样用对数周期轴，峰的位置更直观。
4. 实例：J站年降水量从趋势到周期的完整操作
4.1 数据准备与序列独立性预检
我拿某气象站点（以下记为J站）1961—2020年共60年的逐年年降水量来做示范。真实文件里存为两列：年份、降水（单位mm）。读取和预处理的代码很简单：
data = xlsread('J_rain.xlsx'); % 第一列年份，第二列降水 year = data(:,1); rain = data(:,2);拿到序列后先做三件小事：一是检查缺失值，如果缺测年份超过5%，需要插补或用站间回归估计，不要直接留空跑MK；二是统计序列的一阶自相关系数r1，因为MK检验要求序列相互独立，年降水如果存在明显的持续性，比如连续偏旱或偏涝，会虚增趋势显著性；三是画一张时间序列图，肉眼确认有没有明显的趋势和多年振荡。
实际数据里r1约为0.11，低于0.2的经验阈值，可以认为年际独立性基本满足，不需要做预白处理（pre-whitening）。如果r1超过0.3，我的做法是先做趋势预白处理，把序列减去Sen斜率趋势后再估计自回归系数，不能简单地对原始序列做一阶差分，否则会把真实的低频趋势一并抹掉。
4.2 MK结果与突变点判读
对J站序列运行MK趋势检验，得到Z=-2.13，p=0.033，Sen斜率=-1.38mm/年。H=1说明在0.05显著性水平下，J站年降水量具有显著下降趋势，下降速率约为每年1.38毫米，60年累计趋势幅度约83毫米。这个结论和原始曲线的直观感受一致，但MK给出的数值比肉眼更有说服力。
突变检验方面，UF与UB曲线在1996—1998年区间出现交叉，且交叉点落在±1.96置信带内。再看突变前后的均值差：1961—1996年平均降水508.4mm，1997—2020年平均降水442.7mm，相差约65.7mm。三个证据（UF/UB交叉、均值差、Sen斜率）互相印证，可以判断J站降水量在1990年代中后期存在一个由偏多向偏少的转折。
判读UF/UB曲线时要记住：单次交叉不一定就是突变。我见过不少初学者把置信带外的一系列交叉全标成突变点，那通常是序列内部振荡造成的假象。正确做法是先找置信带内的交叉，再结合两段均值差异确认，最后用滑动t检验或Pettitt检验做交叉验证。多方法互相验证后，突变点才有说服力。
4.3 小波方差图、实部等值线与主周期判断
对J站年降水序列做Morlet小波分析后，小波方差图出现三个明显的峰：第一个峰对应2.8年，第二个峰对应9.6年，第三个峰对应24.5年。峰值最高的是9.6年周期，说明J站降水量最主要的准周期是10年左右；2.8年的准3年周期次之；24年周期虽然存在，但受60年序列长度限制，该周期的可靠度不高。
实部等值线图上，10年尺度的冷暖色斑块在1965—2015年之间基本保持连续，说明这个周期在整个时段都比较稳定。2—4年尺度的斑块则呈断续分布：1960年代、1980年代、2000年代这几个时段能量集中，其他时段不明显，说明准3年周期存在"时段性活跃"的特征，这正好是普通功率谱分析看不出来的信息。
我在结果报告里通常这样表述：J站年降水同时存在准10年主周期和准3年次周期，1990年代中后期存在显著下降突变，1961—2020年总体呈显著减少趋势。一轮"趋势＋突变＋周期"的分析就此闭环。
5. 小波分析结果必踩的坑：边界效应、虚假周期与信度
5.1 COI锥形影响区怎么读，为什么不能只看高值区
连续小波变换默认对有限序列做零填充延拓，边界附近的小波系数会明显失真。判别可靠区的工具是影响锥（COI），代码里我给出了近似半径√2·a。对尺度a，边界影响大约延伸到序列两端各√2·a的时间范围。所以大尺度对应的周期，在序列两端占据了很大的不可靠区域。24年周期的峰虽然在小波方差图里很显眼，但它在两端各影响约34年，几乎覆盖了整个60年序列，只有中间一二十年的数据真正支撑这个周期，结论必须谨慎下。
画图时，建议把COI线叠加到等值线图上，报告只讨论COI以内的显著区。新手最容易犯的错误是把整张图的高能量区域全当真实周期：60年序列里出现30年周期的高能量，其实是边界效应加低频能量泄漏的假象，不落在可靠区内的峰，宁可不写也不要硬解释。
5.2 红噪声信度检验：为什么MK显著不等于小波显著
小波方差图上出现峰，不代表这个周期统计显著。降水序列普遍存在红噪声背景（即AR1过程），其功率谱在低频端自然偏高，即使没有真实周期，低频处也会出现宽缓的峰。如果不对红噪声背景做检验，十个站里能查出一堆假周期。
正规做法是对序列拟合AR1模型，蒙特卡洛模拟数千条红噪声序列做Morlet小波变换，统计每个尺度上的95%置信阈值，真实序列能量超过阈值的区域才算显著周期。手写实现也不复杂：
rng(1); nSim = 1000; threshold = zeros(size(period)); for s = 1:length(period) energies = zeros(nSim,1); for m = 1:nSim noise = ar1_sim(r1, n, std(rain)); % 生成AR1序列 [Wn, ~, ~, ~] = morlet_cwt_self(noise, 1, 6); energies(m) = sum(abs(Wn(s,:)).^2); end threshold(s) = quantile(energies, 0.95); end如果实际序列在某尺度的能量超过对应阈值，就在图上标以细黑线围出显著性区域。这个方法计算量略大但值得，尤其是要发表论文时，没有显著性检验的小波周期很容易被审稿人质疑。日常内部研判时，至少也要拿小波方差峰高和周边背景能量做对比，避免把噪声低频谱误认为周期。
5.3 短序列、缺测、季节干扰的实操处理
序列长度是Morlet小波分析的最大限制。经验上，能可靠分析的最长周期不超过序列长度的三分之一，60年序列谈20年以上周期就只能算参考。如果手里只有30年数据，老实把分析重点放在2—8年周期上，不要硬撑到15年以上，否则就是拿边界效应当规律。
逐月数据的处理也更讲究。直接对月降水量做小波，高频部分会蹦出强烈的12个月周期和1个月噪声，这不是你想要的"年际周期"。必须先减掉逐月气候态，比如公式x_mon_anom = x_mon − climatology，再做小波，才看得到3年、7年这类年际尺度信号。季节分量不剔除就分析，相当于戴着放大镜找低频信号。
缺测数据方面，MK检验对个别缺测不敏感，但小波变换对连续缺测非常敏感，零填充会制造假振荡。连续缺测超过2年，我建议用多年月均值插补，并在报告中注明插补影响；超过5年连续缺测，该段小波结果应当从可靠区里剔除。
最后说一句实在话：MK检验和小波分析都是统计工具，不是结论本身。我跑过几十个站点的降水序列之后，最深的感觉就是"结果要和物理机制互相印证"。一个准3年周期，如果能和区域环流异常、海温异常的变化对上，才值得写进结论；一个突变点，如果有实际观测记录相呼应，才经得起复核。不要为显著性而分析，要为了讲清楚气候变化的真实过程而分析。这也是手里这套Matlab代码反复提醒我的事：先有对数据的诚实判断，再有漂亮的结果图。
