企业官网建设流程全解析

简介：面向气象数据分析与科研学习者的一套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; % 月数据，单位为年，若日数据则 dt=1/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/year，period就是年。如果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导致整片空白

站点数据经常有缺测，如果缺测被填充成了NaN，cwt会直接返回 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这几个要素的读写和绘图逻辑盘清楚，后面学交叉小波和批量分析会快非常多。

本文还有配套的精品资源，点击获取