从拿到一份几十年的降水资料开始大家问得最多的就是两个问题这雨量到底是增了还是减了变化是不是有规律、存在几年一轮的周期光看原始曲线谁也说不清。解决这个问题水文气象领域用得最顺手的组合就是Mann-KendallMK检验加上Morlet小波分析。前者告诉你趋势和突变点后者告诉你显著周期两个方法叠起来看才敢把规律两个字写进报告里。这篇文章我直接把MK检验和Morlet小波分析的原理、Matlab代码、参数选择和判读坑位一次讲透适合正在写课程论文、毕业论文或要处理实测降水资料的同行参考。1. 为什么雨量分析总要叠着MK检验和小波一起看1.1 给多年的雨量序列做体检趋势、突变、周期各查什么降水序列分析其实是三种体检项目的组合体。第一个项目是趋势回答这60年总雨量在长期增加还是减少第二个项目是突变回答有没有哪一年发生了明显转折从偏多变成偏少第三个项目是周期回答是否存在3年、7年或者20年这种准周期性振荡。这三种问题对应的方法很不一样。趋势可以用线性回归做但降水数据往往不满足正态性假设而且个别极端年份暴雨、特旱会把回归斜率拉得很偏。MK检验是秩次检验不要求数据服从正态分布也不怕个别异常值所以成为水文气象趋势判断的事实标准。突变点需要用序贯MKUF/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个观测值对所有ij的组合计算sgn(xj−xi)累加得到S。公式可以简写为$$S \sum_{i1}^{n-1}\sum_{ji1}^{n} \text{sgn}(x_j - x_i)$$这个S的含义很直观上升配对多时S为正下降配对多时S为负。如果序列完全随机S的期望接近0。第二步计算S的方差。当数据不存在重复值时方差是Var(S)n(n−1)(2n5)/18。存在重复值tie时要考虑结组的修正减去Σtp(tp−1)(2tp5)/18其中tp是每组重复值的个数。降水数据非常容易重复比如连续几年都是同一降雨量整数所以修正项不能省。第三步标准化。用连续性修正因子把S转成标准正态Z值$$Z \begin{cases} (S-1)/\sqrt{\text{Var}(S)} S0 \ 0 S0 \ (S1)/\sqrt{\text{Var}(S)} S0 \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检验通常要求n8); end % 计算S S 0; for i 1:n-1 for j i1: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(end1) cnt; end end varS n*(n-1)*(2*n5)/18; if ~isempty(tp) varS varS - sum(tp.*(tp-1).*(2*tp5))/18; end % 标准化 if S 0 Z (S-1)/sqrt(varS); elseif S 0 Z (S1)/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 i1: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 i1:k Sk Sk sign(x(j) - x(i)); end end E k*(k-1)/4; Var k*(k-1)*(2*k5)/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 i1:k Sk Sk sign(y(j) - y(i)); end end E k*(k-1)/4; Var k*(k-1)*(2*k5)/72; UB(n-k1) - (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越大高斯包络在频域越窄频率分辨率越高但时域窗口越宽时间定位越模糊。水文气象里几乎约定俗成取ω06因为在ω06附近小波的时频分辨率和经典Heisenberg不确定性折中得比较好而且此时尺度a与傅里叶周期T有非常简洁的近似关系$$T \approx \frac{4\pi}{\omega_0 \sqrt{2\omega_0^2}}, a \approx 1.033a$$也就是说当ω06时小波尺度数值上几乎等于周期。很多新手在等值线图的纵轴直接标尺度而不换算成周期年这在报告里是不合规的。论文常见表述是准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(2omega0^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.13p0.033Sen斜率-1.38mm/年。H1说明在0.05显著性水平下J站年降水量具有显著下降趋势下降速率约为每年1.38毫米60年累计趋势幅度约83毫米。这个结论和原始曲线的直观感受一致但MK给出的数值比肉眼更有说服力。突变检验方面UF与UB曲线在1996—1998年区间出现交叉且交叉点落在±1.96置信带内。再看突变前后的均值差1961—1996年平均降水508.4mm1997—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代码反复提醒我的事先有对数据的诚实判断再有漂亮的结果图。
从拿到一份几十年的降水资料开始大家问得最多的就是两个问题这雨量到底是增了还是减了变化是不是有规律、存在几年一轮的周期光看原始曲线谁也说不清。解决这个问题水文气象领域用得最顺手的组合就是Mann-KendallMK检验加上Morlet小波分析。前者告诉你趋势和突变点后者告诉你显著周期两个方法叠起来看才敢把规律两个字写进报告里。这篇文章我直接把MK检验和Morlet小波分析的原理、Matlab代码、参数选择和判读坑位一次讲透适合正在写课程论文、毕业论文或要处理实测降水资料的同行参考。1. 为什么雨量分析总要叠着MK检验和小波一起看1.1 给多年的雨量序列做体检趋势、突变、周期各查什么降水序列分析其实是三种体检项目的组合体。第一个项目是趋势回答这60年总雨量在长期增加还是减少第二个项目是突变回答有没有哪一年发生了明显转折从偏多变成偏少第三个项目是周期回答是否存在3年、7年或者20年这种准周期性振荡。这三种问题对应的方法很不一样。趋势可以用线性回归做但降水数据往往不满足正态性假设而且个别极端年份暴雨、特旱会把回归斜率拉得很偏。MK检验是秩次检验不要求数据服从正态分布也不怕个别异常值所以成为水文气象趋势判断的事实标准。突变点需要用序贯MKUF/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个观测值对所有ij的组合计算sgn(xj−xi)累加得到S。公式可以简写为$$S \sum_{i1}^{n-1}\sum_{ji1}^{n} \text{sgn}(x_j - x_i)$$这个S的含义很直观上升配对多时S为正下降配对多时S为负。如果序列完全随机S的期望接近0。第二步计算S的方差。当数据不存在重复值时方差是Var(S)n(n−1)(2n5)/18。存在重复值tie时要考虑结组的修正减去Σtp(tp−1)(2tp5)/18其中tp是每组重复值的个数。降水数据非常容易重复比如连续几年都是同一降雨量整数所以修正项不能省。第三步标准化。用连续性修正因子把S转成标准正态Z值$$Z \begin{cases} (S-1)/\sqrt{\text{Var}(S)} S0 \ 0 S0 \ (S1)/\sqrt{\text{Var}(S)} S0 \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检验通常要求n8); end % 计算S S 0; for i 1:n-1 for j i1: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(end1) cnt; end end varS n*(n-1)*(2*n5)/18; if ~isempty(tp) varS varS - sum(tp.*(tp-1).*(2*tp5))/18; end % 标准化 if S 0 Z (S-1)/sqrt(varS); elseif S 0 Z (S1)/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 i1: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 i1:k Sk Sk sign(x(j) - x(i)); end end E k*(k-1)/4; Var k*(k-1)*(2*k5)/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 i1:k Sk Sk sign(y(j) - y(i)); end end E k*(k-1)/4; Var k*(k-1)*(2*k5)/72; UB(n-k1) - (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越大高斯包络在频域越窄频率分辨率越高但时域窗口越宽时间定位越模糊。水文气象里几乎约定俗成取ω06因为在ω06附近小波的时频分辨率和经典Heisenberg不确定性折中得比较好而且此时尺度a与傅里叶周期T有非常简洁的近似关系$$T \approx \frac{4\pi}{\omega_0 \sqrt{2\omega_0^2}}, a \approx 1.033a$$也就是说当ω06时小波尺度数值上几乎等于周期。很多新手在等值线图的纵轴直接标尺度而不换算成周期年这在报告里是不合规的。论文常见表述是准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(2omega0^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.13p0.033Sen斜率-1.38mm/年。H1说明在0.05显著性水平下J站年降水量具有显著下降趋势下降速率约为每年1.38毫米60年累计趋势幅度约83毫米。这个结论和原始曲线的直观感受一致但MK给出的数值比肉眼更有说服力。突变检验方面UF与UB曲线在1996—1998年区间出现交叉且交叉点落在±1.96置信带内。再看突变前后的均值差1961—1996年平均降水508.4mm1997—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代码反复提醒我的事先有对数据的诚实判断再有漂亮的结果图。
