气象数据小波分析实战：MATLAB实现CWT、XWT与WTC周期诊断

简介这份气象数据小波分析MATLAB代码面向气象、水文及大气科学方向的学生与科研人员用于处理非平稳、非线性气象信号帮助理解小波变换在时间-频率局部化分析中的实际应用。资源包共2个文件包含1个m脚本与1个mat数据文件压缩包约2KB脚本负责实现小波分析主流程数据文件提供暴雨量等实测或模拟序列便于直接运行与验证。代码覆盖小波基选择、小波系数生成、小波方差计算、小波模及模平方求解等关键环节可结合暴雨量数据研究极端降水事件的时空演变与能量分布特征。目前已有5384人学习下载适合希望提升MATLAB编程能力、加深对气象数据多尺度特性理解并需要可复现代码范例的读者参考使用。1. 气象数据小波分析从非平稳序列里挖出周期信号的那把刀做气象数据分析的人迟早会撞上一个尴尬傅里叶变换告诉你这段气温序列里有 365 天的周期但你盯着图看发现这个周期在最近十年明显变弱了而傅里叶变换对此一无所知。原因很简单傅里叶变换把整个时间轴上的信息平均掉了它假设信号的频率成分不随时间变化。可气象数据——气温、降水、风速、气压——恰恰是最典型的非平稳序列周期和趋势都在漂移。小波分析就是冲着这个痛点来的它把时间序列拆成不同尺度的成分同时保留时间信息和频率信息让你能看到“哪个周期在哪个时间段活跃”。这份 MATLAB 代码包围绕气象数据的连续小波变换CWT、小波相干XWT和小波功率谱展开适合有基本 MATLAB 功底、正在做气候诊断或水文周期分析的人。下面我从数据准备一路拆到出图把参数和坑都摆出来。2. 小波变换选型为什么气象领域偏爱 Morlet而不是 Daubechies2.1 连续小波与离散小波的分工小波变换分两大阵营连续小波变换CWT和离散小波变换DWT。很多人第一次接触小波是从 DWT 开始的因为它有快速算法做去噪和压缩很顺手。但气象周期分析几乎清一色用 CWT原因在于 CWT 的输出是一张二维的时频图——横轴时间、纵轴尺度对应周期颜色表示能量强弱。你要找的是“1990 年到 2005 年之间有一个 8 年左右的周期”这种问题只有 CWT 能直接回答。DWT 把信号分解成若干层细节系数和近似系数每一层对应一个频带但时间分辨率随层数下降而且分界频率是固定的不适合追踪周期漂移。代码包里 CWT 是主线DWT 只在预处理阶段出现——比如用wdenoise做序列的初步去噪把观测噪声压掉再送进 CWT。这个顺序不要搞反先 CWT 再去噪会把噪声的能量也摊进时频图里出来的图全是碎斑。2.2 母小波的物理意义决定你能看到什么母小波的选择直接决定时频分辨率的平衡。气象领域最常用的是 Morlet 小波因为它是复小波能同时给出振幅和相位振幅谱看周期强弱相位谱看两个序列的滞后关系。Morlet 的波形是一个高斯包络调制的复指数形状像一段衰减的正弦波和气象序列里那种“持续几个周期的振荡”匹配得很好。对比一下其他常见选择Daubechies 系列是正交小波适合 DWT 做去噪和压缩但它是实小波没有相位信息做 CWT 出来的图不好解读Mexican hat 小波适合检测突变点但频率分辨率差周期定位不准。代码包里默认wname morl如果你做的是降水序列Morlet 同样适用因为降水虽然有零值截断但年际振荡的形态还是接近正弦。注意Morlet 小波不满足容许条件均值为零MATLAB 的cwt函数内部做了修正但如果你自己写卷积实现必须减去均值否则低频部分会发散。2.3 尺度与周期的换算关系CWT 输出的纵轴是尺度a不是周期T。对 Morlet 小波换算公式是T (4 * pi * a) / (omega0 sqrt(2 omega0^2))其中omega0是 Morlet 小波的中心频率MATLAB 默认omega0 6。代入后近似为T ≈ 1.03 * a也就是说尺度和周期几乎一一对应。但如果你改了omega0这个系数会变。代码包里有一个scale2period.m函数专门做这个换算出图时把纵轴从尺度换成周期读图的人才能直接说“这是 5 年周期”而不是“这是尺度 48”。2.4 数据预处理的三个必要步骤气象数据进 CWT 之前有三件事必须做少一件图就废。第一步去趋势。气温序列有全球变暖的线性趋势降水序列可能有站点迁移导致的阶跃。趋势的能量集中在最低频会污染大尺度小波的功率谱让整个时频图底部发红。用detrend做线性去趋势或者用polyfit拟合高阶多项式再减掉。第二步标准化。不同气象要素的量纲不同气温是摄氏度降水是毫米风速是米每秒。CWT 的功率谱是振幅平方量纲不统一没法比较。代码包里统一做zscore标准化变成零均值单位方差。第三步补缺测。气象观测序列常有缺测CWT 对缺测敏感一个 NaN 会让整条小波系数链断掉。常见做法是线性插值补缺如果缺测段超过连续 5 个点插值不可靠考虑用同站其他要素回归补或者直接截断序列。% 气象序列预处理去趋势、标准化、补缺测 raw readmatrix(temp_series.csv); % 读入原始气温序列 t raw(:,1); x raw(:,2); % 补缺测线性插值 nanIdx isnan(x); x(nanIdx) interp1(t(~nanIdx), x(~nanIdx), t(nanIdx), linear, extrap); % 去线性趋势 x_detrend detrend(x, linear); % 标准化 x_norm zscore(x_detrend); % 检查 fprintf(序列长度 %d均值 %.4f标准差 %.4f\n, length(x_norm), mean(x_norm), std(x_norm));这段代码的逻辑是先插值补缺再detrend去掉线性趋势最后zscore标准化。参数说明interp1的linear是插值方法extrap允许端点外推detrend的linear表示去线性趋势如果要同时去季节循环改成先减去月气候态再detrend。输出前打印均值和标准差确认标准化后均值接近 0、标准差接近 1否则检查是否有异常值没处理。3. 跑通 CWT 主流程从 cwt 函数到 cone of influence3.1 用 MATLAB 内置 cwt 还是自己写卷积MATLAB 从 R2016b 开始提供了cwt函数底层用的是 FFT 加速的卷积比自己写循环快一到两个数量级。代码包里主流程用内置cwt但保留了一个cwt_conv.m作为教学对照方便理解小波系数是怎么算出来的。自己写卷积的版本长这样function W cwt_conv(x, scales, omega0) % 简化版 Morlet 连续小波变换仅用于理解原理 n length(x); W zeros(length(scales), n); for i 1:length(scales) a scales(i); % 构造 Morlet 小波母函数 t_wav -4:0.01:4; psi pi^(-1/4) * exp(1i*omega0*t_wav) .* exp(-t_wav.^2/2); % 尺度伸缩 psi_a conj(psi) / sqrt(a); % 卷积 W(i,:) conv(x, psi_a, same); end end逻辑说明对每个尺度a把母小波伸缩后与信号做卷积得到该尺度下的小波系数。参数说明omega0默认取 6scales是尺度向量通常取 2 的幂次序列。这个版本没有做 FFT 加速长序列会慢但每一步都看得见适合调试。实际跑数据用内置cwt% 主流程CWT 计算与出图 fs 1; % 采样间隔年数据为 1 scales 2:0.25:128; % 尺度范围对应周期约 2 到 128 年 wname morl; % Morlet 小波 % 计算 CWT [cfs, f] cwt(x_norm, scales, wname, fs); % 计算小波功率谱 power abs(cfs).^2; % 出图 figure; contourf(t, 1./f, power, 20, linecolor, none); colormap(jet); colorbar; xlabel(年份); ylabel(周期 (年)); title(气温序列小波功率谱);参数说明scales的范围决定你能看到的最小和最大周期下限受采样间隔限制上限一般不超过序列长度的三分之一wname选morlfs是采样频率年数据为 1月数据为 12。cwt返回的f是频率向量1./f换算成周期。出图用contourf填色等值线比imagesc更平滑。3.2 cone of influence边缘效应不是玄学CWT 在序列两端有边缘效应因为卷积核在端点处没有足够的数据支撑小波系数会被人为放大。MATLAB 的cwt函数自动计算了 cone of influenceCOI在图上用一条 U 形曲线标出受影响区域。COI 以外的区域功率谱不可信读图时直接忽略。COI 的边界公式是sqrt(2) * a也就是说尺度越大受影响的边缘越宽。如果你做的是 50 年气温序列128 年尺度的 COI 会覆盖整个序列那个尺度的结果完全没有意义。所以尺度上限一般取序列长度的四分之一到三分之一50 年序列最多看到 15 到 20 年周期。代码包里plot_coi.m函数负责在图上叠加 COI 曲线function plot_coi(t, period, coi) % 在时频图上叠加 cone of influence hold on; plot(t, coi, k--, linewidth, 1.5); hold off; end调用时把cwt返回的coi向量传进去黑色虚线以下区域就是可信区。这个细节很多论文里不标审稿人一眼就能看出来。3.3 显著性检验红噪声背景谱小波功率谱出来以后哪些峰是真实周期哪些是随机涨落需要做显著性检验。气象序列的零假设通常用红噪声AR(1) 过程作为背景谱因为气象数据普遍存在一阶自相关。代码包里wave_signif.m实现了标准的红噪声检验% 红噪声显著性检验 lag1 0.72; % 序列的一阶自相关系数需从数据估计 alpha 0.05; % 显著性水平 [signif, fft_theor] wave_signif(1, 1, scales, lag1, alpha); % 在功率谱图上叠加显著性等值线 hold on; contour(t, 1./f, power, signif, k, linewidth, 1.2); hold off;参数说明lag1是一阶自相关系数用corr(x(1:end-1), x(2:end))估计alpha取 0.05 或 0.1signif是每个尺度对应的显著性阈值。图上黑色等值线包围的区域就是通过 95% 显著性检验的周期。注意显著性检验只说明“这个峰不太可能是红噪声产生的”不代表物理机制别过度解读。4. 交叉小波与相干分析两个序列的周期怎么对上4.1 交叉小波变换的物理含义单序列 CWT 告诉你“气温有 5 年周期”但你想知道“这个 5 年周期和太平洋海温的 5 年周期是不是同步的”。这时候需要交叉小波变换XWT。XWT 的定义是两个序列 CWT 的乘积再取共轭W_xy W_x * conj(W_y)|W_xy|表示两个序列在某个时间和尺度上共同拥有的能量相位角angle(W_xy)表示两个序列的滞后关系。如果相位差为 0说明同相如果为 pi/2说明一个领先另一个四分之一周期。代码包里xwt.m函数封装了这个计算function [Wxy, period, coi] xwt(x, y, fs) % 交叉小波变换 scales 2:0.25:128; [Wx, f] cwt(x, scales, morl, fs); [Wy, ~] cwt(y, scales, morl, fs); Wxy Wx .* conj(Wy); period 1./f; coi sqrt(2) * scales / fs; end逻辑说明分别对两个序列做 CWT然后逐点相乘取共轭。参数说明x和y必须等长且已经过预处理fs采样频率输出Wxy是复数矩阵abs(Wxy)是交叉功率谱angle(Wxy)是相位谱。4.2 小波相干归一化后的耦合强度XWT 的缺点是功率谱受两个序列自身振幅影响如果气温振幅大、海温振幅小XWT 会偏向气温。小波相干WTC做了归一化R_xy |S(W_xy)| / sqrt(S(|W_x|^2) * S(|W_y|^2))S是平滑算子在时间和尺度两个方向做滑动平均。R_xy取值在 0 到 1 之间表示两个序列在时频面上的局部相关系数。WTC 比 XWT 更适合判断“耦合强度”因为它是无量纲的。代码包里wtc.m实现了平滑和归一化function [R, period] wtc(x, y, fs) % 小波相干 [Wxy, period, ~] xwt(x, y, fs); [Wx, ~] cwt(x, 2:0.25:128, morl, fs); [Wy, ~] cwt(y, 2:0.25:128, morl, fs); % 时间和尺度方向平滑 S_Wxy smooth2d(Wxy, 3, 1); S_Wx2 smooth2d(abs(Wx).^2, 3, 1); S_Wy2 smooth2d(abs(Wy).^2, 3, 1); R abs(S_Wxy) ./ sqrt(S_Wx2 .* S_Wy2); end参数说明smooth2d是自定义的二维滑动平均函数第一个参数是平滑窗口在尺度方向的大小第二个是时间方向的大小。窗口越大相干图越平滑但细节丢失越多。一般尺度方向取 3 到 5时间方向取 1 到 2。4.3 相位箭头怎么读WTC 图上通常叠加相位箭头箭头方向表示两个序列的相位差。箭头向右表示同相向左表示反相向上表示 y 领先 x 四分之一周期向下表示 x 领先 y。读图时注意箭头只在相干值高的区域有意义相干值低于 0.5 的区域箭头是噪声。代码包里phase_arrows.m在相干图上每隔几个点画一个箭头function phase_arrows(t, period, R, phase, skip) % 在相干图上叠加相位箭头 [T, P] meshgrid(t, period); idx_t 1:skip:length(t); idx_p 1:skip:length(period); quiver(T(idx_p, idx_t), P(idx_p, idx_t), ... cos(phase(idx_p, idx_t)), sin(phase(idx_p, idx_t)), ... 0.5, k, linewidth, 0.8); end参数说明skip控制箭头密度一般取 4 到 60.5是箭头缩放因子。箭头太密会糊成一片太疏看不出趋势。5. 避坑与排查那些让时频图变成废图的细节5.1 现象功率谱底部一片红看不出任何周期原因序列没去趋势低频能量淹没了所有周期信号。气温序列的线性增暖趋势在 CWT 里表现为最大尺度上的高功率整个图底部发红把 5 到 10 年周期的信号盖住了。解决进 CWT 之前强制detrend。如果去线性趋势后底部还是红检查是否有阶跃或异常值用findchangepts找突变点分段去均值。5.2 现象小波相干图全图相干值都很高接近 1原因两个序列都有强自相关平滑窗口又开得太大导致归一化后的相干值虚高。这是小波相干分析里最常见的翻车点。解决减小平滑窗口尺度方向从 5 降到 3时间方向从 2 降到 1。同时做蒙特卡洛检验生成 1000 对同长度的 AR(1) 序列算它们的 WTC取 95 分位数作为显著性阈值只有超过这个阈值的相干值才可信。5.3 现象COI 曲线以下区域出现强信号以为是新发现原因COI 以下区域受边缘效应影响小波系数被人为放大功率谱不可信。很多人不看 COI 直接读图把边缘效应当成周期信号。解决出图时强制叠加 COI 曲线读图只看 COI 以上区域。如果目标周期落在 COI 以下说明序列长度不够要么延长序列要么降低尺度上限。5.4 现象换了omega0参数后周期位置对不上原因尺度和周期的换算系数依赖omega0。MATLAB 默认omega0 6换算系数约 1.03如果改成omega0 8系数变成约 1.3同样的尺度对应的周期变了。解决固定omega0 6或者每次改omega0后重新用scale2period换算纵轴。不要在不同omega0的图之间直接比较周期位置。5.5 现象月数据直接当年数据跑周期全错原因采样频率fs没设对。年数据fs 1月数据fs 12日数据fs 365。如果月数据用fs 1CWT 算出来的“周期”单位是天不是月整个纵轴偏移。解决进 CWT 之前确认fs与数据时间单位一致。代码包里check_fs.m函数根据时间列自动推断采样频率推断逻辑是取时间差的中位数再取倒数。6. 进阶技巧把 CWT 结果做成可复现的周期诊断报告跑通单次 CWT 只是起点实际项目里你要对比多个站点、多个要素、多个时间段。我一般会把整个流程封装成一个函数输入站点列表和要素名输出每个站点的周期诊断表。核心思路是用cwt的功率谱做峰值检测自动提取显著周期。function report period_diagnosis(data_cell, station_names, fs) % 批量周期诊断输入多个站点的序列输出显著周期表 n_station length(data_cell); report table(); for i 1:n_station x data_cell{i}; x detrend(x, linear); x zscore(x); [cfs, f] cwt(x, 2:0.25:128, morl, fs); power abs(cfs).^2; period 1./f; % 对时间维平均得到全局功率谱 global_power mean(power, 2); % 找峰值功率大于两侧相邻点的局部极大值 [pks, locs] findpeaks(global_power, MinPeakProminence, 0.1); % 取前三个显著周期 [~, sortIdx] sort(pks, descend); top_n min(3, length(sortIdx)); for j 1:top_n idx locs(sortIdx(j)); report [report; {station_names{i}, period(idx), pks(sortIdx(j))}]; end end report.Properties.VariableNames {Station, Period_Year, Power}; end逻辑说明对每个站点序列做预处理后跑 CWT沿时间维平均得到全局功率谱用findpeaks找局部极大值按功率排序取前三个周期。参数说明MinPeakProminence控制峰值的最小突出度设太小会检出噪声峰设太大会漏掉弱周期一般取全局功率谱标准差的 0.5 倍。输出是一张表每个站点三行分别是前三个显著周期及其功率。这个报告的好处是可复现同样的输入数据跑出来的周期表完全一致不会因为手动读图产生主观偏差。我一般还会在报告里加一列 COI 标记如果某个周期落在 COI 以下标注为“不可信”避免后续分析踩坑。从那以后我每次做周期诊断都强制走一遍“去趋势 → 标准化 → CWT → 红噪声检验 → COI 过滤 → 峰值提取”这条链中间任何一步跳过后面的结论都不敢用。希望帮到你。本文还有配套的精品资源点击获取