简介这份MATLAB代码资源面向气象、水文及环境科学领域的学生与研究人员用于对非平稳气象信号开展小波分析。代码围绕小波基选择、小波系数生成、小波方差计算、小波模及模平方求解等环节展开并以暴雨量数据为例演示时间-频率局部化分析流程帮助使用者识别不同时间尺度上的极端降水事件与能量分布特征。资源包共2个文件包含1个m脚本与1个mat数据文件压缩包约2KB脚本负责核心分析逻辑数据文件提供可复现的暴雨量样本便于直接运行与调试。目前已有5384人学习下载适合希望将小波理论落地到气象数据处理、提升MATLAB编程与结果解读能力的读者参考。1. 气象数据小波分析到底在分析什么做气象数据处理的人迟早会撞上这样一个场景手里攥着一段气温、风速或降水的长序列想做趋势项和周期项的分离FFT 一上去就发现频谱糊成一团——因为气象序列几乎都是非平稳的季节波动、年际振荡、突变过程叠在一起傅里叶变换那套「全局频率」的假设直接失效。这时候气象数据小波分析就成了绕不开的工具而 MATLAB 因为自带 Wavelet Toolbox、矩阵操作顺手、绘图即出图成了很多人落地这套方法的第一选择。这篇笔记面向两类人一类是刚拿到气象站观测数据、想快速跑出小波功率谱和显著性检验的新手另一类是已经用过cwt但被边界效应、红噪声检验、尺度参数折腾过的熟手。我会把从数据预处理到出图的完整链路拆开讲重点放在参数怎么设、坑在哪、结果怎么读而不是复述小波变换的数学定义。整套流程用 MATLAB 原生函数加少量自写脚本就能跑通不需要额外工具箱之外的依赖。2. 从原始气象序列到可分析的小波输入2.1 为什么气象序列必须先做这几步预处理小波分析对输入序列的「干净程度」比 FFT 更敏感原因是它同时在时间轴和尺度轴上做局部化任何一个缺测值或异常跳变都会在多个尺度上产生虚假能量。常见的气象数据来源无非是站点观测、再分析资料、卫星反演三类它们共同的毛病是缺测NaN 或哨兵值如 -999、时间轴不均匀尤其早期人工观测、量纲不统一气温 ℃ 与 K、风速 m/s 与 knot。我一般按这个顺序处理先统一时间轴到等间隔月尺度或日尺度再处理缺测最后去趋势。去趋势这一步很多人跳过结果小波功率谱在最大尺度上翘起来误以为是长周期信号其实是趋势项泄漏。去趋势用线性拟合或滑动平均都行关键是要在论文或报告里写清楚用了哪种因为它直接影响低频段的解释。% 假设原始数据为两列datenum 和 value raw readmatrix(station_data.csv); t raw(:,1); % MATLAB datenum x raw(:,2); % 1) 缺测值替换哨兵值 -999 转 NaN x(x -999) NaN; % 2) 时间轴等间隔化重采样到月尺度 t0 datenum(1980-01-01); t1 datenum(2020-12-01); tq (t0:31:t1); % 近似月步长实际项目建议用 calmonths xq interp1(t, x, tq, linear, NaN); % 3) 缺测插补线性插值 首尾外推 nanIdx isnan(xq); xq(nanIdx) interp1(tq(~nanIdx), xq(~nanIdx), tq(nanIdx), linear, extrap); % 4) 去趋势减去最小二乘线性趋势 p polyfit(tq, xq, 1); xdetr xq - polyval(p, tq); % 5) 标准化零均值单位方差便于跨变量比较 xn (xdetr - mean(xdetr)) / std(xdetr);这段代码里几个参数值得说清楚。interp1的linear是最保守的选择气象序列不建议用样条插值因为样条会在缺测段制造出虚假的平滑振荡小波一分析就露馅。polyfit的阶数固定为 1如果你研究的是降水这种本身没有明显线性趋势的量可以跳过第 4 步但气温、海温这类有增暖背景的变量必须去趋势。标准化不是必须的但如果你要对比不同站点的功率谱标准化能让纵轴量纲统一省去后面反复解释。2.2 采样间隔和序列长度怎么定才不出错小波分析的尺度范围直接受序列长度和采样间隔约束。一个经验规则最大可分析尺度不要超过序列长度的三分之一否则边界效应会吃掉大半结果。比如 40 年月度数据480 个点最大尺度控制在 160 个月以内比较稳妥再往上功率谱两端的可信度就明显下降。采样间隔决定了你能分辨的最小尺度。月度数据理论上能分辨 2 个月但实际中 2 个月尺度的能量往往被噪声淹没我一般从 3 个月起画。日尺度数据可以下探到 2 天但要注意日气象序列有强烈的 1 天周期这个周期会在小波谱上形成一个贯穿所有尺度的竖条分析前最好先做日平均或带阻滤波去掉。还有一个容易被忽略的点序列长度最好是 2 的整数次幂附近不是硬性要求但用 FFT 加速的cwt在长度接近 2^n 时速度明显更快。480 个点不算快如果你有 1000 个站点要批量跑把序列截断到 512 或补零到 512 能省不少时间。补零要谨慎它会在两端引入虚假能量截断比补零安全。3. 用 MATLAB 跑通连续小波变换的最小命令集3.1 cwt 和 wcoherence 该选哪个MATLAB 里做小波分析主要有两条路cwt做单序列的连续小波变换wcoherence做两个序列的交叉小波和相干性分析。选哪个取决于你的科学问题。如果只是想知道单个气象要素有哪些周期用cwt如果想看两个要素比如气温和降水、海温和季风指数在哪些时段哪些尺度上同步用wcoherence。cwt的默认小波是 Morse 小波比老版本的 Morlet 更灵活但如果你要跟文献对比很多经典气象小波论文用的是 Morlet这时候需要显式指定。指定方法是在cwt里传Wavelet参数。Morse 小波的默认参数(3, 20)对大多数气象序列够用但如果你关注的是窄带信号比如 ENSO 的 3-7 年Morlet 的时间分辨率更好代价是频率分辨率下降。% 连续小波变换单序列 Fs 1; % 采样频率月度数据取 1单位1/月 wname morl; % 或 morse跟文献对比时用 morl cwt(xn, Fs, wname); % 指定尺度范围单位月 scales 3:1:160; cwt(xn, scales, wname); % 交叉小波相干两个序列 % y 需与 x 等长、同采样 wcoherence(xn, yn, Fs);cwt直接调用会出图但做批量分析时我一般用输出参数拿数据自己画因为默认图的配色和标注在报告里往往要改。cwt的输出包括系数矩阵wt、频率向量f、以及用于显著性检验的coi影响锥和sig9595% 置信水平。这三个东西是后面出图的核心缺一个图就不完整。wcoherence的输出是相干谱和相位谱相位谱用箭头表示两个序列的领先滞后关系读图时箭头向右表示同相向左表示反相向上表示 y 领先 x 四分之一周期。这个相位读法是交叉小波最容易搞错的地方很多人把箭头方向理解反了结论就全拧了。3.2 尺度序列怎么设线性还是对数尺度序列的设置直接影响图的纵轴分布。线性尺度比如 3:1:160在纵轴上等间距低频段被压缩高频段被拉伸对数尺度logspace让纵轴按对数分布低频段展开高频段压缩。气象文献里两种都有但周期分析通常用对数尺度因为周期本身是倍频关系对数轴更能体现 2 年、4 年、8 年这种倍频结构。% 对数尺度从 3 个月到 160 个月共 60 个尺度 scales logspace(log10(3), log10(160), 60); % 对应的实际周期月 periods scales; % 对 Morlet周期约等于尺度 % 画图时纵轴用 log2 刻度更直观 figure; wt cwt(xn, scales, morl); imagesc(tq, log2(scales), abs(wt)); set(gca, YDir, normal, YTick, log2([3 6 12 24 48 96 160]), ... YTickLabel, {3,6,12,24,48,96,160}); ylabel(周期 (月)); xlabel(时间); colorbar;这里有个细节Morlet 小波的「尺度」和「周期」在数值上近似相等但严格来说有一个转换系数约 1.03。做定性分析可以忽略做定量周期估计时最好乘上这个系数或者直接用cwt输出的频率向量f取倒数得到周期避免手算误差。imagesc画图时纵轴用log2(scales)而不是scales这样倍频周期在图上等间距读图时 2 年、4 年、8 年的位置一目了然。4. 显著性检验和影响锥小波分析最容易翻车的地方4.1 红噪声背景谱怎么算小波功率谱的显著性检验和 FFT 一样需要一个背景谱作为零假设。气象序列的默认零假设是红噪声AR(1) 过程因为大多数气象要素有自相关性白噪声假设太强会把大量真实信号误判为显著。红噪声背景谱由序列的 lag-1 自相关系数决定这个系数必须从去趋势后的序列算用原始序列算会偏大。% 计算 lag-1 自相关系数 xlag xn(1:end-1); xcur xn(2:end); r corrcoef(xlag, xcur); r1 r(1,2); % 红噪声背景谱简化版Morlet % 完整实现需按 Torrence Compo (1998) 的公式 n length(xn); dt 1; freq 1 ./ scales; alpha (2 r1) / (1 r1); % 红噪声谱指数 Pk (1 - r1^2) / (1 - 2*r1*cos(2*pi*freq*dt) r1^2);这段代码是简化版完整的红噪声检验需要按 Torrence Compo 的公式算每个尺度上的 95% 置信阈值再和实际功率谱比较。MATLAB 的cwt输出里已经带了sig95直接用就行不用自己重算。但你要知道sig95是怎么来的否则图上看不懂为什么有些区域被判定显著、有些没有。r1的计算有个坑如果序列本身有强周期比如明显的年循环lag-1 自相关会被周期信号污染算出来的红噪声背景偏高显著性检验偏保守。这时候可以先做带阻滤波去掉年循环再算r1或者改用更复杂的背景谱模型。我一般会同时画出去年循环和不去年循环两个版本对比看结论是否稳健。4.2 影响锥为什么必须画在图上影响锥Cone of Influence, COI是小波分析区别于 FFT 的核心概念之一。小波变换在序列两端需要补零或镜像延拓导致两端的大尺度结果不可信。COI 就是标出这个不可信区域的边界线通常画成一条从两端向中间收拢的弧线弧线以外的区域功率谱被边界效应污染不能用于物理解释。很多新手出的图不画 COI审稿人一眼就能看出来。更隐蔽的错误是画了 COI 但没把 COI 以外的区域遮掉或打阴影读者容易误读。我习惯用半透明阴影把 COI 以外区域盖住或者直接把那部分功率谱设为 NaN 不显示。% cwt 输出自带 coi [wt, f, coi] cwt(xn, scales, morl); % 画图时叠加 COI figure; pcolor(tq, log2(scales), abs(wt)); shading interp; hold on; plot(tq, log2(coi), w--, LineWidth, 1.5); % 把 COI 以外区域打阴影 for i 1:length(tq) if ~isnan(coi(i)) fill([tq(i) tq(i) tq(min(i1,end)) tq(min(i1,end))], ... [log2(coi(i)) max(log2(scales)) max(log2(scales)) log2(coi(i))], ... [0.5 0.5 0.5], FaceAlpha, 0.3, EdgeColor, none); end endCOI 的形状取决于小波类型和尺度序列。Morlet 的 COI 在两端收拢得比较快Morse 稍慢。如果你的结论严重依赖 COI 边缘附近的信号那这个结论本身就不够稳健应该缩短分析时段或换更长的序列重做。5. 避坑与排查气象小波分析里最常见的五个翻车点5.1 现象功率谱在最大尺度上翘起误判为长周期信号原因序列没去趋势或者去趋势不彻底。线性趋势在小波谱上表现为最大尺度处的能量堆积看起来像一个超长周期。另一个可能是序列两端有突变比如仪器更换导致的跳变小波在两端补零时把这个跳变放大。解决先画原始序列确认有没有趋势或跳变有趋势就做线性或多项式去趋势有跳变就分段处理或做均一化校正。去趋势后重跑如果最大尺度能量明显下降说明之前是趋势泄漏。5.2 现象显著性区域连成一片几乎所有尺度都显著原因红噪声背景谱的 lag-1 自相关系数算错了或者序列本身有强自相关但没做预白化。还有一种可能是序列长度太短自由度不足显著性检验本身不可靠。解决检查r1是否在合理范围气象序列一般 0.3-0.8如果接近 1 说明序列接近随机游走红噪声假设不适用。序列长度少于 100 个点时显著性检验结果只能作参考不要下强结论。5.3 现象交叉小波的相位箭头方向读反领先滞后结论颠倒原因wcoherence的相位定义是 y 相对 x 的相位差箭头向右表示同相向上表示 y 领先 x 四分之一周期。很多人凭直觉把箭头方向理解成时间流向结论就反了。解决拿两个已知相位关系的合成序列比如 y sin(t pi/2)跑一遍wcoherence看箭头方向建立正确的直觉。实际分析时相位结论最好和互相关函数或滞后回归交叉验证不要只靠小波相位下结论。5.4 现象批量跑几百个站点时内存爆掉或速度极慢原因cwt对每个序列都重新计算小波基和尺度序列批量循环时重复开销大。另外cwt默认输出完整系数矩阵长序列加多尺度时矩阵很大。解决把尺度序列预先算好循环里复用只保留需要的输出功率谱和 COI不要保留完整系数矩阵序列长度截断到 2 的整数次幂附近加速 FFT。如果站点超过 500 个考虑用parfor并行但要注意每个 worker 的内存占用。5.5 现象换一个小波类型或尺度范围结论就变了原因小波分析的结果对参数敏感这是方法本身的特性不是 bug。Morlet 和 Morse 的时间-频率分辨率不同线性尺度和对数尺度在图上呈现不同结论自然会有差异。解决做敏感性测试至少用两种小波类型和两种尺度范围各跑一遍只报告在所有配置下都稳健的信号。如果某个周期只在特定参数下出现那它很可能是方法 artifact不是真实信号。这一步很多论文省略但恰恰是审稿人最关注的。6. 把结果做成能发表的图配色、标注和批量出图技巧出图这件事做过的人都知道它是玄学重灾区。同样的数据配色差一点、COI 画法差一点审稿人的观感能差出两个档次。我踩过的坑包括用jet配色被批「不专业」功率谱用线性色标导致弱信号全看不见COI 阴影盖住了显著区域。配色上功率谱推荐用感知均匀的色标MATLAB 自带的parula或第三方cmocean里的thermal、balance都行。避免jet和hsv这两个在功率谱上会产生虚假的边界感。色标范围用caxis手动设一般取功率谱的 5% 到 95% 分位数避免极端值把整体对比度压掉。% 批量出图每个站点一张统一色标 stations {S001,S002,S003}; % 虚构站点编号 clim [0 1.5]; % 统一色标范围需先统计所有站点 for k 1:length(stations) fn sprintf(%s_monthly.mat, stations{k}); S load(fn); xn S.xn; [wt, f, coi] cwt(xn, scales, morl); P abs(wt).^2; figure(Visible,off,Position,[100 100 900 400]); pcolor(tq, log2(scales), P); shading interp; caxis(clim); colormap(cmocean(thermal)); hold on; plot(tq, log2(coi), k--, LineWidth, 1.2); set(gca, YDir, normal, YTick, log2([3 6 12 24 48 96 160]), ... YTickLabel, {3,6,12,24,48,96,160}, FontSize, 11); xlabel(年份); ylabel(周期 (月)); title(sprintf(站点 %s 小波功率谱, stations{k})); colorbar; print(gcf, sprintf(%s_cwt.png, stations{k}), -dpng, -r300); close(gcf); end批量出图的关键是统一色标和统一尺度范围否则站点之间没法对比。clim要先跑一遍所有站点统计全局分位数再定不能每个站点各设各的。Visible设为off避免弹窗print用-r300保证分辨率够投稿。如果图要进 LaTeX 或 Word建议同时输出 PDF 矢量图print的-dpdf选项即可。最后一个习惯每张图旁边存一份参数记录写清楚小波类型、尺度范围、去趋势方法、红噪声参数。半年后你回头看图没有这份记录根本想不起来当时怎么设的。这个习惯帮我省过好几次返工希望帮到你。本文还有配套的精品资源点击获取