简介面向气象科研与教学场景这份Matlab代码演示了基于小波变换的气象时序数据分析方法适合需要处理非平稳、非线性信号的研究人员或学生快速上手。资源压缩包仅2KB共2个文件一个M脚本负责完整小波分析流程包括小波基选择、小波系数与方差计算、小波模及模平方求解代码结构清晰、注释简明便于按步骤理解小波分析的数学实现另一个MAT数据文件提供了实际暴雨量样本可直接加载运行并对照结果进行气象解释。通过简单修改脚本参数可迁移至温度、湿度、风速等变量的多尺度变化研究帮助识别不同时间尺度上的极端事件与能量分布。目前已有5382人学习使用对于想掌握Matlab小波工具箱在气象领域应用的人来说是一份精炼实用的入门范例。 干气象数据分析这些年小波分析可以说是绕不开的一个工具。不管是做降水周期识别、气温突变检测还是研究海温异常这类大尺度过程的年代际变化小波分析都能帮我们同时回答“什么周期在变强”和“这个周期在什么时候最强”这两个问题。MATLAB凭借成熟的Wavelet Toolbox和相对简洁的语法是我做这块最顺手的实现环境。这篇文章不绕弯子直接结合我的实际项目讲讲怎么用MATLAB对气象数据进行小波分析从数据预处理、小波变换实现、显著性检验到绘图再到实际踩过的几个大坑尽量一次说清楚。1. 气象数据为什么要用小波分析1.1 傅里叶变换的先天短板做气象序列分析大家最先接触的往往是傅里叶变换。它能把一条时间序列拆成不同频率的正弦波叠加告诉你信号里有哪些主要周期。比如对一条40年的月降水量序列做FFT可以算出“年周期最强、3个月周期其次”之类的结果。问题是傅里叶变换丢掉了时间信息。它默认信号的频率成分在整个时间段内是恒定不变的可气象序列恰恰不是这样。以长江中下游夏季降水为例它的年际变化可能与El Niño事件有关但这种关系在某些年代很强某些年代又几乎消失。如果只用傅里叶变换你只能得到“存在一个2-7年的周期”却完全看不出这个周期在1970年代强、1980年代弱、90年代又恢复——这对于气候变化研究几乎是致命的。还有一个实际问题很多气象序列是非平稳的统计特征随时间漂移。直接对整条序列做傅里叶变换等于是把所有时段的信息平均在一起。哪怕中间出现过一次极端事件它在频谱里也只是一根不起眼的小峰很容易被淹没在背景噪声里。1.2 小波分析能回答气象领域的哪些问题小波分析的核心思路是拿一个有限长度的波形小波母函数去和信号做内积然后把这个波形在时间轴上平移、在尺度轴上伸缩得到一张“时间-频率”二维图谱。这张图上每个位置的数值就代表“在某个时刻某个频率成分的强度”。做气象分析时小波能直接回答这样几类问题周期随时间如何变化。比如ENSO的2-7年周期在不同年代有没有强弱变化是否在某些时段断掉。突变发生的大致时点和尺度。比如某地气温在1990年代前后是否发生了显著性跃变这种跃变在哪个周期尺度上最明显。两个气象要素之间的关系强度随时间怎么变。比如降水与海温的小波相干性哪些时间段两者强相关哪些时间段脱钩。常见的气象小波分析场景包括降水和气温序列的周期识别、干湿变化阶段的划分、海气耦合信号的年代际演变、甚至空气质量浓度序列的周期提取。这些场景下小波分析往往比传统滤波或滑动平均更能保留多尺度信息。2. MATLAB小波工具箱选型与代码框架2.1 新版cwt函数与老工具箱的区别MATLAB的Wavelet Toolbox经过几次更新小波函数接口变化比较大。老版本R2015b左右以前常用cwt(x, scales, morl)这种写法返回一个复数系数矩阵如果调用了cwtft参数又不一样。新版本R2016b之后推荐直接用cwt(x, fs, amor)用法直观很多% 新版cwt写法 fs 12; % 采样频率月数据就是每年12个点 [wt, f] cwt(x, fs, amor); % wt为小波系数f为对应的频率 period 1 ./ f; % 频率转周期这里amor是解析Morlet小波气象领域最常用。如果不指定小波名称MATLAB默认用morse效果也可以但物理意义和论文里常见表述不够统一。我建议做气象周期分析时显式指定amor这样代码可读性好结果也容易和文献对照。老接口也不是不能用但我在实际工作中踩过一个坑旧代码在R2021a上跑cwt默认参数变了输出的系数矩阵含义不一样画出来的图“看着正常但尺度轴对不上”。所以如果你在复习旧脚本务必先确认doc cwt里的说明而不是直接复制老代码跑。2.2 连续小波变换还是离散小波变换气象分析里周期识别和时频特征分析优先选连续小波变换CWT去噪和信号重构选离散小波变换DWT。两者选型依据可以参考下面这张表对比项连续小波变换(CWT)离散小波变换(DWT)核心用途时频功率谱、周期识别、显著性检验多分辨率分解、去噪、重构输出结果时间-频率二维系数图不同尺度分解分量(AD)频率分辨率高适合精细周期分析二进网格相对粗糙常用函数cwt、wcoherencewavedec、wrcoef、wdenoise气象典型应用降水周期、ENSO年代际变化趋势提取、信号去噪、突变检测预处理我见过有人用DWT去替代CWT做周期显著性检验结果因为DWT的频带是二进划分的2-8年的周期被硬切到好几个尺度里导致谱峰变宽、显著性判断失真。反过来用CWT做去噪也不合适它的基函数之间存在冗余重构过程麻烦且容易引入伪震荡。这个选型问题看起来简单实际项目里却经常成为后续分析误差的根源。3. 气象数据小波分析完整实操3.1 数据预处理宁可多花三十分钟也不要直接跑小波气象数据预处理是小波分析里最容易被跳过、但影响最大的环节。我见过不少代码读取数据后直接对原始序列做小波变换出来的谱图“看着热闹”实际上全是季节循环和趋势项的影子。我的标准流程分三步第一步缺失值处理。气象站数据经常有缺口尤其是日降水、风速这类序列。小波变换理论上不允许序列中有NaN你需要先做插补。常用的方法有线性插值、三次样条插值、或者用临近站点回归插补。对月尺度气候分析我个人偏好三次样条插值它不会像线性插值那样产生明显折角。第二步去趋势和去季节循环。如果序列里有长期趋势比如全球变暖导致的升温趋势这个趋势会在小波功率谱的低频端形成巨大能量泄漏把真正的年代际周期全部盖住。如果数据有明显的季节循环则会在1年周期处形成一个超高峰相关尺度附近的信号都看不清。处理办法是对温度数据先做“月距平”每月减去该月的多年平均值对降水数据可考虑做标准化距平再对处理后的序列做去趋势。第三步方差稳定化或标准化。小波功率谱对序列的绝对数值很敏感为了便于比较不同站点、不同要素我通常会做Z-score标准化% 去季节循环示例月均温度数据 temp_anom temp - repmat(mean(reshape(temp, 12, []), 2), 1, length(temp)/12); temp_anom temp_anom(:); % 去趋势 x detrend(temp_anom); % 标准化 x (x - mean(x)) / std(x);需要说明的是如果你做的是降水数据且想保留“量级”信息比如对比不同区域降水变率可以不去标准化但必须要去季节循环。具体取舍看你研究的科学问题。3.2 连续小波变换代码与参数选择预处理完成后核心代码其实很短% 参数设置 fs 12; % 月数据 wname amor; % 解析Morlet小波 % 连续小波变换 [wt, f] cwt(x, fs, wname); period 1 ./ f; % 周期单位月 power abs(wt).^2; % 小波功率 % 绘制小波功率谱 figure(Color, w); imagesc(t, period, power); set(gca, YScale, log, YDir, reverse); ylim([2 120]); ylabel(周期月); xlabel(时间年); title(小波功率谱); h colorbar; ylabel(h, 功率); colormap(jet);这里几个参数的考虑采样频率fs月数据是12日数据是365年数据是1。cwt内部会用fs把频率轴映射到物理频率不要想当然填1。周期范围我不建议直接套用默认范围。先看序列长度40年月数据能分析的可靠周期上限大约是N/(2倍)左右也就是10-20年以内。超出这个范围的结果基本落在影响锥COI外不可信。小波函数实验对比下来amor最稳。bump的频域局部化更好但时间局部化差一些适合频率分离要求高的时候morse作为通用默认也够用但参数解释起来比较绕。做气候诊断我基本只用amor。3.3 显著性检验设置与结果解读小波功率谱上的亮点不代表就通过了显著性检验。气象序列大多不是白噪声尤其是温度和降水存在明显的自相关也就是红噪声特征。直接拿白噪声背景做检验会把一堆虚假周期判成显著。标准的做法是Torrence和Compo1998提出的方法用AR(1)红噪声过程生成背景谱然后在每个尺度上做卡方检验。MATLAB官方工具箱没有直接给出这个函数但有两条路可以走路径1自己写简化版。对于经验丰富的用户可以直接调用lag1计算滞后一阶自相关系数再用解析公式构造红噪声背景谱最后用卡方分布计算95%置信线。这个过程不复杂但公式细节多容易出错。路径2使用开源小波包。我一般用现有实现再做二次封装比如经典的wavelet_coherence开源代码里面包含了Torrence和Compo的显著性检验逻辑MATLAB和Python版本都有。接入自己的数据后修改数据格式即可。% 伪代码调用开源小波包做显著性检验 % 假设已获得 wt、period、x % [signif, fft_theor] wavelet_signif(x, dt, scales, lag1, 0.95, morlet); % sig95 (signif) * ones(1, length(t)); % 扩展为等时间维度的置信线 % 然后把 sig95 叠加到功率谱图上只显示通过检验的区域说实话显著性检验这部分是我反复试错最多的地方。早先我直接拿MATLAB的cwt跑完功率谱就发论文审稿人一眼看出问题“95%置信区域呢”后来补上了红噪声检验图件质量完全不是一个档次。强烈建议读者不要跳过这一步。3.4 小波全谱与主周期提取小波功率谱图能直观看到亮色区域但要定量说“哪个周期显著”还得看小波全谱Global Wavelet Spectrum。它就是把每个频率上的小波功率沿时间方向取平均相当于把二维谱压成一维识别整体最强的周期% 计算小波全谱 global_power mean(power, 2); % 绘制全谱 figure(Color, w); semilogx(period, global_power, k-, LineWidth, 1.5); xlabel(周期月); ylabel(平均小波功率); grid on;这里有一点务必要注意小波全谱的“平均”是把边界附近不可信数据也算进去的。如果某个时间段边缘效应区域占比很大全谱峰值就会失真。严谨做法是先对每列做COI裁剪只在COI范围内求平均。MATLAB的cwt不直接返回COI掩膜你可以通过计算傅里叶周期对应的衰减来近似或者用开源包返回的COI数组。主周期提取后还能进一步看特定周期段的时变特征。比如你发现2-8年的周期显著可以在小波系数里把周期位于2-8年的分量取出来求平均功率得到一条随时间变化的曲线用来分析“该周期段的强度在哪些年代增强”。这种尺度带平均功率曲线是很多气候诊断论文的标准图件。4. 常见问题与排查心得4.1 边缘效应到底影响多大小波变换本质上是在每个时间点上用有限长度的小波窗口去匹配信号所以在序列两端窗口有一部分落在已知数据之外结果不可信。这个区域叫影响锥COI。在TORRENCE论文的图里COI通常画成一条带状阴影落在里面的功率谱基本不能解释。我自己的经验是数据长度至少要有目标最大周期的4-6倍COI的影响才算可控。比如你想分析年代际周期约8-20年那序列至少要有80-120年短了的话宁可把分析周期上限降低。你如果只有30年月数据硬着头皮看10年以上周期结论基本站不住脚审稿人一眼就会发现。MATLAB中画COI需要根据小波函数的尺度自动计算。官方cwt没有直接输出COI我一般用工具函数算好再叠加。如果不想深究开源包里通常有现成的coi输出直接hold on; plot(t, coi, w--)就可以。4.2 小波功率谱“糊成一片”怎么办这是新手最常见的视觉问题整张图都是红红黄黄的一大片看不出哪里有明显的局部峰。造成这个问题的原因通常有三个第一个原因是数据没去趋势。线性趋势会在低频端产生巨大能量导致整个谱图的颜色映射被低频端主导中高频细节全部被压缩成同一颜色。解决办法就是先用detrend或者对趋势项做去除。第二个原因是颜色映射范围太大。我建议用caxis或新版clim把颜色范围设到功率值的合适分位数范围比如取95%分位数作为上限这样主体结构才能显出来。第三个原因是小波函数选得不对。如果用莫尔斯小波且时间带宽积设置不当时频分辨率的平衡就会变差。回到amor这种默认参数往往能缓解。4.3 Morlet小波参数omega0怎么调老版本cwt(x, scales, morl, omega0)可以设置Morlet小波的带宽参数omega0。这个参数决定了小波在频率域和时间域的平衡默认值是6此时频率分辨率较好时间分辨率适中是气象领域公认的合理选择。如果你把omega0调大比如8或者10频率分辨率更高能更好地区分相近周期但时间定位会变差突变的起始时间会模糊调小则反之。实际项目中我基本不调这个参数一直用6。只有在需要区分非常接近的周期比如8年和10年时才会考虑在对比实验里调整一下。要注意的是改了omega0之后显著性检验的背景谱也要重新对应计算否则前后不一致结果没法解释。4.4 两个序列的相关性/相干性分析怎么实现有时你需要分析两个气象要素的关系随时间变化比如降水与海温、气温与气压。此时要用小波相干性Wavelet Coherence而不是分别画两张小波功率谱。MATLAB里对应函数是wcoherence% 计算两个序列的小波相干 [wcoh, wcs, f] wcoherence(x, y, fs, amor);这个函数输出三个变量相干系数wcoh、交叉小波谱wcs和频率轴f。绘制相干谱图非常方便figure(Color, w); imagesc(t, period, wcoh); set(gca, YScale, log, YDir, reverse); ylim([2 120]); colorbar;小波相干解读要特别小心相干系数高不等于因果关系成立它只说明两个序列在该时段该频段存在一致的相位行为且窗口越宽相干值越容易被高估。做结论前一定要结合相位箭头和物理过程一起看不能拿一张相干热力图直接说“A导致B”。我在这上面吃过亏后来结论都改得很谨慎。5. 从单要素分析到业务化应用的扩展思路说了这么多其实小波分析的价值不只是画一张漂亮的谱图。我在实际项目中逐步形成了一套标准扩展流程分享出来供参考第一步周期诊断。用CWT和显著性检验确认主要周期形成“该要素有哪几个显著周期带、每个周期带的时段演变”的基本判断。第二步滤波重构。在确认周期带后可以用DWT或带通滤波把特定周期分量提取出来分析它在时间轴上的高低相位变化。比如提取出准两年震荡分量后再去和同期的其他要素对比。第三步多要素联合。用wcoherence做两两相干分析找到关系最稳定的频段和时间段再结合物理机制去解释。第四步业务化展示。小波谱图不太直观面向决策者时我会把它简化成“某周期强度时间序列”叠加到普通时间序列图上用简单的折线图说明“2015年后2-4年周期增强”比直接甩一张彩色热力图有效得多。这套流程在气候变化诊断、极端事件归因、空气质量周期识别等方向都跑得通。说起来小波分析只是工具核心还是你对气象物理过程的理解——工具给出证据物理判断才是最终落点。希望这篇实操总结能让你少走几步弯路把时间多留给真正需要思考的科学问题。本文还有配套的精品资源点击获取