简介面向气象数据分析与科研学习者的一套MATLAB小波分析示例代码聚焦非平稳气象信号的时间-频率局部化处理可帮助理解小波基选择、小波系数计算、小波方差与模平方等核心概念并应用于实际降水序列。资源包为RAR压缩格式共2个文件含1个M脚本和1个MAT数据文件脚本实现完整分析流程数据文件提供如暴雨量等实测样本便于直接运行与结果对照整体仅2KB轻量易用。截至目前已有5382人学习下载适合需要入门或快速复现气象小波分析流程的学生、研究人员。通过运行这套代码可掌握从数据载入到小波模、方差解读的完整链路同时提升MATLAB编程能力与对气象数据多尺度变化特征的洞察力。 做气象数据诊断的时候我经常要回答一个问题某段降水序列的周期在哪个时段表现出来过是2到4年还是8到16年这类问题用傅里叶变换也能答但答得不完整——傅里叶只告诉你整段序列里有哪些周期成分却无法告诉你这些周期分别在什么时候出现。小波分析正好补上这个缺口它把时间序列映射到“时间-频率”二维平面能同时看清楚“什么周期”和“哪个时段”。MATLAB的小波工具箱让这件事落地容易了很多但网上很多代码只给一个画图脚本参数设置、显著性检验、结果解读全都没说透。这篇整理了一份完整的实操流程从数据预处理到显著性检验再到结果判读适合气象、水文、气候等方向刚开始接触小波分析的同学。1. 为什么读气象序列不能只靠傅里叶变换气象要素序列几乎都是非平稳的。降水的季节变化、年际变化、年代际变化叠在一起而且各个周期成分的强度会随时间改变。比如同一个站点的夏季降水可能1980年代以前以准3年周期为主1990年代以后变成了准6年周期。傅里叶变换做的是全局频谱分解它把整段序列当成平稳信号处理最后给出的能量分布是所有时段叠加后的平均结果。遇到这种“周期随时间漂移”的情况傅里叶只能看到两个周期的平均态看不到这段历史演变。小波分析本质上是给傅里叶变换加了一个“可移动的窗口”同时通过伸缩和平移母小波来控制窗口的尺度和位置。对于一维时间序列 (x(t))连续小波变换的定义是[ W(a,b)\int x(t)\psi^*{a,b}(t)dt,\quad \psi{a,b}(t)\frac{1}{\sqrt{a}}\psi\left(\frac{t-b}{a}\right) ]其中 (a) 是尺度参数控制小波的伸缩对应周期大小(b) 是平移参数控制小波在时间轴上的位置。(|W(a,b)|^2) 就是小波功率谱表示某个尺度在某个时刻的能量强度。气象里最常用的是复Morlet小波因为它是复数小波能同时提供振幅和相位信息适合检测振荡周期。MATLAB里的amor就是解析Morlet小波。它的时频局部化特性比较均衡时间分辨率和频率分辨率不会像某种窗函数那样顾此失彼又不会让谱图平滑到失去细节。实际经验是如果主要关心周期结构和显著性Morlet基本够用换成Morse小波结果差别不大但参数调整会多一点初学阶段没必要额外折腾。所以小波分析不是要替代傅里叶而是解决“周期是否随时间变化”这个问题。读完整幅小波功率谱你能看到的不只是一个平均周期而是一段完整的时频演变史这对研究气候模态的阶段性变化特别重要。2. 数据准备缺失值、季节循环与标准化这些细节小波分析对输入序列的质量其实挺敏感数据没处理好后面画图全是坑。最常见的数据来源有站点观测、格点再分析资料、模式输出格式通常是 NetCDF、CSV、Excel。第一步是确定你要分析哪个变量比如月平均气温、月降水量然后提取成一条等间隔的时间序列。等间隔是硬要求如果原始观测有缺测需要用fillmissing或interp1做插值。% 假设 data 是包含时间 t 和变量值 x 的表格 t data.time; x data.value; % 如果时间间隔不均匀先插值到等间隔月序列 dt_year 1/12; t_daily_start t(1); t_daily_end t(end); time_regular (t(1):dt_year:t(end)); x_regular interp1(t, x, time_regular, linear); x_regular fillmissing(x_regular, linear);对于月数据建议先去掉季节循环得到距平序列。如果直接用原始月值去做小波分析季节尺度12个月周期的功率会非常强把年际尺度的细节压在下面导致我们真正关心的2到8年信号被淹没。去掉季节循环的办法很简单对每个月份做多年平均得到该月的“气候态”再用原始值减去它。N length(x_regular); months repmat((1:12), N/12, 1); clim accumarray(months, x_regular, [], nanmean); anom x_regular - clim(months);有的研究还建议做标准化anom (anom - mean(anom)) / std(anom);标准化不会改变周期结构但会改变功率谱的绝对数值而且后续红噪声显著性检验需要用到方差所以做标准化不会影响是否显著反而让不同变量之间的谱图可以横向比较。另外一个容易忽略的步骤是去趋势。长序列如果存在明显的气候趋势会在小波功率谱的低频端产生一个很大的能量区看起来像存在超长周期但其实只是线性趋势带来的伪信号。最简单的处理是anom_detrend detrend(anom);。如果研究目标本身就是低频趋势那就不要去掉这要看你要回答什么科学问题。总的来说预处理的原则是只滤掉会干扰目标尺度的成分不要过度处理否则低频信号被滤掉反而看不到年代际周期。3. 核心实现cwt函数、功率谱和红噪声显著性检验3.1 小波基函数选择用MATLAB做连续小波变换核心函数是cwt属于Wavelet Toolbox。老代码里经常看到cwtft那是旧接口建议新代码直接用cwt参数更简洁返回结果也更直观。cwt支持多个解析小波气象分析我固定选amor也就是复Morlet。原因前面说过它能把时频分辨率调到比较均衡的位置而且给出的小波系数是复数后续可以算相位也能很自然地转成功率谱。3.2 cwt函数的基本调用假设预处理后的序列叫ts时间步长是dt调用一行就搞定dt 1/12; % 月数据单位为年若日数据则 dt1/365.25 [wt, f, coi] cwt(ts, amor, dt, VoicesPerOctave, 10); power abs(wt).^2; period 1./f;这里的第四行参数VoicesPerOctave表示每个倍频程内部划分的尺度数默认是10数值越大谱图在频率方向越平滑计算量也略增。对月尺度到年尺度的气象数据10到12比较合适再高也不会有本质变化。返回的wt是一个复数矩阵行对应频率列对应时间f是物理频率向量单位由dt决定coi是影响锥边界绘图时要配合使用。有个容易踩的点cwt返回的f是频率不是周期。画小波图时大家习惯纵轴显示周期所以要做一次period 1./f。频率单位跟着dt走dt用的是年f就是 cycles/yearperiod就是年。如果dt设成了1个月周期单位就变成月坐标轴直接差12倍。3.3 显著性检验与红噪声背景谱小波功率谱图单看颜色深浅不够还需要知道哪些区域是统计显著的。气象序列通常可以用一阶自回归过程AR(1)也叫红噪声作为零假设。做法是先估计序列的滞后1自相关系数 (\rho)然后对每个频率构造红噪声背景谱再按卡方分布求出置信水平阈值。对于复Morlet小波小波功率的统计分布近似于自由度为2的卡方分布所以95%置信阈值是[ \hat{p}{95}(f) \frac{1}{2} \chi^2{0.95}(2) \cdot P_{red}(f) ]对应MATLAB实现alpha 0.95; rho corr(ts(1:end-1), ts(2:end)); sig_level zeros(size(f)); for k 1:length(f) red_noise (1 - rho^2) ./ ... (1 - 2*rho*cos(2*pi*f(k)*dt) rho^2); sig_level(k) std(ts)^2 * red_noise * chi2inv(alpha, 2) / 2; end % 生成二维显著性掩膜 sig95 repmat(sig_level(:), 1, length(ts)); mask power sig95;注意std(ts)^2对应方差如果提前做了标准化这里就是1前面的理论谱计算会自动缩放。rho用的是整条序列的自相关系数这是很多文献的标准做法虽然不是最精细的局部背景估计但胜在稳健、可复现。绘图时可以把显著性区域叠加成黑色等值线time (0:length(ts)-1)*dt; figure(Color,w); contourf(time, period, power, 40, LineStyle, none); set(gca, YScale, log, YDir, reverse); % 大周期在上 ylim([0.5, 20]); % 根据研究目标调整单位年 xlabel(Time (year)); ylabel(Period (year)); colormap(parula); colorbar; hold on; plot(time, 1./coi, w--, LineWidth, 1.5); contour(time, period, mask, [1 1], k-, LineWidth, 0.8);这里纵轴用了对数坐标并且YDir设为reverse让长周期显示在图上方符合气象文献习惯。contourf的功率值我倾向于显示log2(power)或者power.^(0.5)因为原始功率谱的动态范围太大低频端可能把整个色标拉偏。但显著性掩膜mask必须在原始power上计算别用变换后数值去比。4. 结果图判读影响锥和显著性区域到底怎么看一幅完整的小波功率谱图横轴是时间纵轴是周期对数坐标颜色代表能量强度黑色等值线圈出显著区域白色虚线是影响锥。拿到图以后第一步不是找颜色最红的位置而是先找到影响锥边界。coi表示因数据截断在两端产生的边界效应区域。对于小波变换序列开头和结尾附近的小波系数是不可靠的因为母小波在边缘处没有完整覆盖数据周期越长的信号受边界影响的范围越大所以影响锥在低频端会迅速收窄看起来像个漏斗。图中白色虚线的绘制我用的是1./coi因为纵轴是周期。如果直接把coi画进去位置会完全错乱这是新手最容易犯的错之一。影响锥以内的功率颜色再深也不能当作真实周期信号只有在锥内且通过显著性检验的区域才值得解读。换句话说你要找的是黑线等值线圈住、且位于白色虚线“安全区”内的连续深色区域。读图顺序我一般是这样的先看“有哪些时间段出现了显著的周期能量”。如果某段时间内有一块显著区域横跨某个周期范围说明该时段内存在准周期振荡。再看这个区域的周期中心大致在哪范围是窄还是宽。窄说明信号很规律宽说明是频带更宽的准周期过程这本身也有物理含义。最后看这个显著时段的开始和结束年份与已知气候事件或指数变化对比。比如降水小波图上1965到1985年出现2到3年显著周期可能会联系到某种遥相关模态的年代际转变。还有一点要记住统计显著不等于物理真实。红噪声检验只说明“这个能量超过随机背景的可能性达到95%”但样本长度有限多重比较会增加假阳性概率。尤其是周期接近序列长度一半以上的低频信号即使检验显示显著也容易受趋势和边界效应影响稳妥的做法是同时做敏感性分析比如换用不同预处理方案看显著区域是否稳定。5. 那些让我重画了三遍图的小波分析细节5.1 时间步长单位搞反我最早用cwt分析月降水数据dt直接填了1结果周期坐标轴单位变成“月”。当时图上所有峰值周期都在12、24、36我心想这不就是季节周期吗但数据明明已经去掉季节循环了。后来查了一下输出频率单位才意识到dt是采样周期而不是采样频率。月数据想得到以年为单位的周期dt必须写1/12。这个错一次就能记住但很多人会在这里栽跟头。5.2 coi画法的坑前面提过coi是频率轴上的边界不是周期轴上的边界。第一版绘图脚本我直接plot(time, coi)结果白色虚线的位置压在高频区域连图像形状都不对。翻文档发现coi的单位和f相同绘图时做1./coi才和period对齐。还有个细节是coi的长度应该和时间点数一致如果发现维度对不上检查一下是不是cwt版本差异。5.3 数据里有NaN导致整片空白站点数据经常有缺测如果缺测被填充成了NaNcwt会直接返回 NaN 区域画出来的图一块白一块红。这不算报错但非常容易当成小波分析结果不稳定。解决方法是预处理阶段用fillmissing做插值或者用前向填充。如果缺口太多插值会产生假信号建议放弃该站点或者用更稳健的插值方法。我在处理某区域降水格点时就因为一个格点缺测太多没处理结果整个区域的小波图低频段全空浪费了半天时间。5.4 色标范围被低频功率带偏小波功率谱的绝对值往往集中在低频端或趋势区域直接用原始功率画图色标会被一两个大值拉满其他时频区域全部发蓝细节根本看不清。我现在的习惯是用log2(power)或者开根号画填色图能有效拉开中小功率的对比度。显著性掩膜仍然在原始功率上计算所以统计判断不受影响。5.5 频率向量里的极值周期cwt自动计算的频率范围很宽最高频对应极小尺度低频甚至可能接近序列长度。如果ylim从min(period)设到max(period)高频噪声和超低频边界效应会占据整个图面肉眼只能看到一片巨大的“八”字形色块。建议根据研究问题把周期范围限制在合理区间比如分析月降水时限制在0.5到30年分析日数据气温时限制在2到64天。配合对数坐标图面会干净得多。6. 从单变量到多变量交叉小波与相干小波扩展单变量小波功率谱回答的是“这个变量自己有哪些周期模态以及何时显著”但气象研究中很多问题是双变量关系降水与某个环流指数的耦合周期是多少两类指数之间的相关性在什么时候较强这就要从单变量小波谱走向交叉小波谱和小波相干。交叉小波谱类似于傅里叶交叉谱的时频版本它用两个变量的小波系数乘积来定义反映两个序列在某个时频域上共同能量更高的区域。MATLAB 中可以直接用wcoherence计算小波相干它能给出0到1的相干系数适合分析两个序列在特定周期上的相关强度随时间的变化。基本调用[Wxy, f, coi] wcoherence(x_anom, y_anom, dt, VoicesPerOctave, 10);其中x_anom和y_anom是等长、已做距平和标准化的序列dt与cwt用法一致。返回的相干系数可以直接画填色图同样需要注意影响锥和显著性。wcoherence内部有平滑操作结果比裸计算交叉谱更稳健但绘图前最好阅读文档确认coi的输出形式。如果想进一步分析两个序列之间的滞后关系可以看小波相位差。用cwt分别得到两个复小波系数矩阵wt_x和wt_y相位差就是angle(wt_x .* conj(wt_y))。相位差的箭头图在交叉小波分析中很常见能直观看出哪个序列在某个时频区域领先另一个领先多少周期。不过这个分析需要先掌握单变量cwt的系数结构否则很容易把二维矩阵的维度搞混。再往外扩展还可以对多站点或多格点批量做小波分析把每个格点的显著周期统计出来画成空间分布图。比如提取每个格点小波功率谱上最显著周期对应的时间段再插值成空间图能看出某区域周期模态的空间不均匀性。这些扩展都建立在前面单变量流程跑通的基础上。先把cwt、power、sig95、coi这几个要素的读写和绘图逻辑盘清楚后面学交叉小波和批量分析会快非常多。本文还有配套的精品资源点击获取 MATLAB小波分析实战：从傅里叶到时频域，完整流程与显著性检验
简介面向气象数据分析与科研学习者的一套MATLAB小波分析示例代码聚焦非平稳气象信号的时间-频率局部化处理可帮助理解小波基选择、小波系数计算、小波方差与模平方等核心概念并应用于实际降水序列。资源包为RAR压缩格式共2个文件含1个M脚本和1个MAT数据文件脚本实现完整分析流程数据文件提供如暴雨量等实测样本便于直接运行与结果对照整体仅2KB轻量易用。截至目前已有5382人学习下载适合需要入门或快速复现气象小波分析流程的学生、研究人员。通过运行这套代码可掌握从数据载入到小波模、方差解读的完整链路同时提升MATLAB编程能力与对气象数据多尺度变化特征的洞察力。 做气象数据诊断的时候我经常要回答一个问题某段降水序列的周期在哪个时段表现出来过是2到4年还是8到16年这类问题用傅里叶变换也能答但答得不完整——傅里叶只告诉你整段序列里有哪些周期成分却无法告诉你这些周期分别在什么时候出现。小波分析正好补上这个缺口它把时间序列映射到“时间-频率”二维平面能同时看清楚“什么周期”和“哪个时段”。MATLAB的小波工具箱让这件事落地容易了很多但网上很多代码只给一个画图脚本参数设置、显著性检验、结果解读全都没说透。这篇整理了一份完整的实操流程从数据预处理到显著性检验再到结果判读适合气象、水文、气候等方向刚开始接触小波分析的同学。1. 为什么读气象序列不能只靠傅里叶变换气象要素序列几乎都是非平稳的。降水的季节变化、年际变化、年代际变化叠在一起而且各个周期成分的强度会随时间改变。比如同一个站点的夏季降水可能1980年代以前以准3年周期为主1990年代以后变成了准6年周期。傅里叶变换做的是全局频谱分解它把整段序列当成平稳信号处理最后给出的能量分布是所有时段叠加后的平均结果。遇到这种“周期随时间漂移”的情况傅里叶只能看到两个周期的平均态看不到这段历史演变。小波分析本质上是给傅里叶变换加了一个“可移动的窗口”同时通过伸缩和平移母小波来控制窗口的尺度和位置。对于一维时间序列 (x(t))连续小波变换的定义是[ W(a,b)\int x(t)\psi^*{a,b}(t)dt,\quad \psi{a,b}(t)\frac{1}{\sqrt{a}}\psi\left(\frac{t-b}{a}\right) ]其中 (a) 是尺度参数控制小波的伸缩对应周期大小(b) 是平移参数控制小波在时间轴上的位置。(|W(a,b)|^2) 就是小波功率谱表示某个尺度在某个时刻的能量强度。气象里最常用的是复Morlet小波因为它是复数小波能同时提供振幅和相位信息适合检测振荡周期。MATLAB里的amor就是解析Morlet小波。它的时频局部化特性比较均衡时间分辨率和频率分辨率不会像某种窗函数那样顾此失彼又不会让谱图平滑到失去细节。实际经验是如果主要关心周期结构和显著性Morlet基本够用换成Morse小波结果差别不大但参数调整会多一点初学阶段没必要额外折腾。所以小波分析不是要替代傅里叶而是解决“周期是否随时间变化”这个问题。读完整幅小波功率谱你能看到的不只是一个平均周期而是一段完整的时频演变史这对研究气候模态的阶段性变化特别重要。2. 数据准备缺失值、季节循环与标准化这些细节小波分析对输入序列的质量其实挺敏感数据没处理好后面画图全是坑。最常见的数据来源有站点观测、格点再分析资料、模式输出格式通常是 NetCDF、CSV、Excel。第一步是确定你要分析哪个变量比如月平均气温、月降水量然后提取成一条等间隔的时间序列。等间隔是硬要求如果原始观测有缺测需要用fillmissing或interp1做插值。% 假设 data 是包含时间 t 和变量值 x 的表格 t data.time; x data.value; % 如果时间间隔不均匀先插值到等间隔月序列 dt_year 1/12; t_daily_start t(1); t_daily_end t(end); time_regular (t(1):dt_year:t(end)); x_regular interp1(t, x, time_regular, linear); x_regular fillmissing(x_regular, linear);对于月数据建议先去掉季节循环得到距平序列。如果直接用原始月值去做小波分析季节尺度12个月周期的功率会非常强把年际尺度的细节压在下面导致我们真正关心的2到8年信号被淹没。去掉季节循环的办法很简单对每个月份做多年平均得到该月的“气候态”再用原始值减去它。N length(x_regular); months repmat((1:12), N/12, 1); clim accumarray(months, x_regular, [], nanmean); anom x_regular - clim(months);有的研究还建议做标准化anom (anom - mean(anom)) / std(anom);标准化不会改变周期结构但会改变功率谱的绝对数值而且后续红噪声显著性检验需要用到方差所以做标准化不会影响是否显著反而让不同变量之间的谱图可以横向比较。另外一个容易忽略的步骤是去趋势。长序列如果存在明显的气候趋势会在小波功率谱的低频端产生一个很大的能量区看起来像存在超长周期但其实只是线性趋势带来的伪信号。最简单的处理是anom_detrend detrend(anom);。如果研究目标本身就是低频趋势那就不要去掉这要看你要回答什么科学问题。总的来说预处理的原则是只滤掉会干扰目标尺度的成分不要过度处理否则低频信号被滤掉反而看不到年代际周期。3. 核心实现cwt函数、功率谱和红噪声显著性检验3.1 小波基函数选择用MATLAB做连续小波变换核心函数是cwt属于Wavelet Toolbox。老代码里经常看到cwtft那是旧接口建议新代码直接用cwt参数更简洁返回结果也更直观。cwt支持多个解析小波气象分析我固定选amor也就是复Morlet。原因前面说过它能把时频分辨率调到比较均衡的位置而且给出的小波系数是复数后续可以算相位也能很自然地转成功率谱。3.2 cwt函数的基本调用假设预处理后的序列叫ts时间步长是dt调用一行就搞定dt 1/12; % 月数据单位为年若日数据则 dt1/365.25 [wt, f, coi] cwt(ts, amor, dt, VoicesPerOctave, 10); power abs(wt).^2; period 1./f;这里的第四行参数VoicesPerOctave表示每个倍频程内部划分的尺度数默认是10数值越大谱图在频率方向越平滑计算量也略增。对月尺度到年尺度的气象数据10到12比较合适再高也不会有本质变化。返回的wt是一个复数矩阵行对应频率列对应时间f是物理频率向量单位由dt决定coi是影响锥边界绘图时要配合使用。有个容易踩的点cwt返回的f是频率不是周期。画小波图时大家习惯纵轴显示周期所以要做一次period 1./f。频率单位跟着dt走dt用的是年f就是 cycles/yearperiod就是年。如果dt设成了1个月周期单位就变成月坐标轴直接差12倍。3.3 显著性检验与红噪声背景谱小波功率谱图单看颜色深浅不够还需要知道哪些区域是统计显著的。气象序列通常可以用一阶自回归过程AR(1)也叫红噪声作为零假设。做法是先估计序列的滞后1自相关系数 (\rho)然后对每个频率构造红噪声背景谱再按卡方分布求出置信水平阈值。对于复Morlet小波小波功率的统计分布近似于自由度为2的卡方分布所以95%置信阈值是[ \hat{p}{95}(f) \frac{1}{2} \chi^2{0.95}(2) \cdot P_{red}(f) ]对应MATLAB实现alpha 0.95; rho corr(ts(1:end-1), ts(2:end)); sig_level zeros(size(f)); for k 1:length(f) red_noise (1 - rho^2) ./ ... (1 - 2*rho*cos(2*pi*f(k)*dt) rho^2); sig_level(k) std(ts)^2 * red_noise * chi2inv(alpha, 2) / 2; end % 生成二维显著性掩膜 sig95 repmat(sig_level(:), 1, length(ts)); mask power sig95;注意std(ts)^2对应方差如果提前做了标准化这里就是1前面的理论谱计算会自动缩放。rho用的是整条序列的自相关系数这是很多文献的标准做法虽然不是最精细的局部背景估计但胜在稳健、可复现。绘图时可以把显著性区域叠加成黑色等值线time (0:length(ts)-1)*dt; figure(Color,w); contourf(time, period, power, 40, LineStyle, none); set(gca, YScale, log, YDir, reverse); % 大周期在上 ylim([0.5, 20]); % 根据研究目标调整单位年 xlabel(Time (year)); ylabel(Period (year)); colormap(parula); colorbar; hold on; plot(time, 1./coi, w--, LineWidth, 1.5); contour(time, period, mask, [1 1], k-, LineWidth, 0.8);这里纵轴用了对数坐标并且YDir设为reverse让长周期显示在图上方符合气象文献习惯。contourf的功率值我倾向于显示log2(power)或者power.^(0.5)因为原始功率谱的动态范围太大低频端可能把整个色标拉偏。但显著性掩膜mask必须在原始power上计算别用变换后数值去比。4. 结果图判读影响锥和显著性区域到底怎么看一幅完整的小波功率谱图横轴是时间纵轴是周期对数坐标颜色代表能量强度黑色等值线圈出显著区域白色虚线是影响锥。拿到图以后第一步不是找颜色最红的位置而是先找到影响锥边界。coi表示因数据截断在两端产生的边界效应区域。对于小波变换序列开头和结尾附近的小波系数是不可靠的因为母小波在边缘处没有完整覆盖数据周期越长的信号受边界影响的范围越大所以影响锥在低频端会迅速收窄看起来像个漏斗。图中白色虚线的绘制我用的是1./coi因为纵轴是周期。如果直接把coi画进去位置会完全错乱这是新手最容易犯的错之一。影响锥以内的功率颜色再深也不能当作真实周期信号只有在锥内且通过显著性检验的区域才值得解读。换句话说你要找的是黑线等值线圈住、且位于白色虚线“安全区”内的连续深色区域。读图顺序我一般是这样的先看“有哪些时间段出现了显著的周期能量”。如果某段时间内有一块显著区域横跨某个周期范围说明该时段内存在准周期振荡。再看这个区域的周期中心大致在哪范围是窄还是宽。窄说明信号很规律宽说明是频带更宽的准周期过程这本身也有物理含义。最后看这个显著时段的开始和结束年份与已知气候事件或指数变化对比。比如降水小波图上1965到1985年出现2到3年显著周期可能会联系到某种遥相关模态的年代际转变。还有一点要记住统计显著不等于物理真实。红噪声检验只说明“这个能量超过随机背景的可能性达到95%”但样本长度有限多重比较会增加假阳性概率。尤其是周期接近序列长度一半以上的低频信号即使检验显示显著也容易受趋势和边界效应影响稳妥的做法是同时做敏感性分析比如换用不同预处理方案看显著区域是否稳定。5. 那些让我重画了三遍图的小波分析细节5.1 时间步长单位搞反我最早用cwt分析月降水数据dt直接填了1结果周期坐标轴单位变成“月”。当时图上所有峰值周期都在12、24、36我心想这不就是季节周期吗但数据明明已经去掉季节循环了。后来查了一下输出频率单位才意识到dt是采样周期而不是采样频率。月数据想得到以年为单位的周期dt必须写1/12。这个错一次就能记住但很多人会在这里栽跟头。5.2 coi画法的坑前面提过coi是频率轴上的边界不是周期轴上的边界。第一版绘图脚本我直接plot(time, coi)结果白色虚线的位置压在高频区域连图像形状都不对。翻文档发现coi的单位和f相同绘图时做1./coi才和period对齐。还有个细节是coi的长度应该和时间点数一致如果发现维度对不上检查一下是不是cwt版本差异。5.3 数据里有NaN导致整片空白站点数据经常有缺测如果缺测被填充成了NaNcwt会直接返回 NaN 区域画出来的图一块白一块红。这不算报错但非常容易当成小波分析结果不稳定。解决方法是预处理阶段用fillmissing做插值或者用前向填充。如果缺口太多插值会产生假信号建议放弃该站点或者用更稳健的插值方法。我在处理某区域降水格点时就因为一个格点缺测太多没处理结果整个区域的小波图低频段全空浪费了半天时间。5.4 色标范围被低频功率带偏小波功率谱的绝对值往往集中在低频端或趋势区域直接用原始功率画图色标会被一两个大值拉满其他时频区域全部发蓝细节根本看不清。我现在的习惯是用log2(power)或者开根号画填色图能有效拉开中小功率的对比度。显著性掩膜仍然在原始功率上计算所以统计判断不受影响。5.5 频率向量里的极值周期cwt自动计算的频率范围很宽最高频对应极小尺度低频甚至可能接近序列长度。如果ylim从min(period)设到max(period)高频噪声和超低频边界效应会占据整个图面肉眼只能看到一片巨大的“八”字形色块。建议根据研究问题把周期范围限制在合理区间比如分析月降水时限制在0.5到30年分析日数据气温时限制在2到64天。配合对数坐标图面会干净得多。6. 从单变量到多变量交叉小波与相干小波扩展单变量小波功率谱回答的是“这个变量自己有哪些周期模态以及何时显著”但气象研究中很多问题是双变量关系降水与某个环流指数的耦合周期是多少两类指数之间的相关性在什么时候较强这就要从单变量小波谱走向交叉小波谱和小波相干。交叉小波谱类似于傅里叶交叉谱的时频版本它用两个变量的小波系数乘积来定义反映两个序列在某个时频域上共同能量更高的区域。MATLAB 中可以直接用wcoherence计算小波相干它能给出0到1的相干系数适合分析两个序列在特定周期上的相关强度随时间的变化。基本调用[Wxy, f, coi] wcoherence(x_anom, y_anom, dt, VoicesPerOctave, 10);其中x_anom和y_anom是等长、已做距平和标准化的序列dt与cwt用法一致。返回的相干系数可以直接画填色图同样需要注意影响锥和显著性。wcoherence内部有平滑操作结果比裸计算交叉谱更稳健但绘图前最好阅读文档确认coi的输出形式。如果想进一步分析两个序列之间的滞后关系可以看小波相位差。用cwt分别得到两个复小波系数矩阵wt_x和wt_y相位差就是angle(wt_x .* conj(wt_y))。相位差的箭头图在交叉小波分析中很常见能直观看出哪个序列在某个时频区域领先另一个领先多少周期。不过这个分析需要先掌握单变量cwt的系数结构否则很容易把二维矩阵的维度搞混。再往外扩展还可以对多站点或多格点批量做小波分析把每个格点的显著周期统计出来画成空间分布图。比如提取每个格点小波功率谱上最显著周期对应的时间段再插值成空间图能看出某区域周期模态的空间不均匀性。这些扩展都建立在前面单变量流程跑通的基础上。先把cwt、power、sig95、coi这几个要素的读写和绘图逻辑盘清楚后面学交叉小波和批量分析会快非常多。本文还有配套的精品资源点击获取