简介这套Matlab小波周期分析实操包面向水文、气象、环境及地理等领域的科研与教学用户围绕Morlet小波时频分析全流程解决时间序列周期性检测、多尺度振荡识别与主周期提取等常见需求。压缩包共20个文件、4.33MB以m运行脚本、xls/xlsx/txt数据样本、jpg可视化图谱、avi操作录像及doc文档为主内含年径流量鲁台子、唐乃海与年降温数据等真实时间序列便于直接替换数据验证方法。包内提供两个主运行脚本适配Matlab 2021a及以上版本并配套清晰AVI操作录像手把手演示数据导入、参数设置、系数提取、图形绘制与结果解读doc文档补充方法原理与扩展接口提示帮助理解小波系数模平方、方差图与主周期判读。已有12人学习下载适合需要快速上手小波周期分析的水文气象初学者及相关课程实验。 Matlab做周期分析,我最早用的是快速傅里叶变换,后来和小波周期分析一比,才明白为什么大家都在写“小波”的论文。原因很简单:FFT告诉你“这个序列有显著的2年周期”,小波分析却能告诉你“这个2年周期只在1965到1975年显著,1980年之后就不明显了”。对径流、气温这类受气候波动影响的序列来说,这种“周期随时间变化”的信息,恰恰比单纯一个周期值重要得多。这篇内容是我整理Matlab小波周期分析实操包时的完整笔记,包含径流、气温两类时间序列的完整案例、可视化图谱的生成与判读,以及我在反复调试中遇到的各种坑。适合水文、气象、地理、生态方向的硕博生和科研搬砖党参考。如果你正在处理降水、径流、NDVI、气温这类时间序列,想提取主周期却不想被一堆公式劝退,这篇文章应该能帮你少走不少弯路。1. 为什么做周期分析要从小波变换入手1.1 傅里叶变换的局限:周期是全局的,但水文序列是非平稳的先聊一个接触过周期分析的人都会遇到的问题。FFT给出的频谱是整段时间的“平均结果”,它默认信号在整个时间范围内周期性是恒定不变的。但径流序列根本不是这么回事——长江宜昌站年径流序列里,丰枯变化受副高、季风、太阳活动、厄尔尼诺共同影响,某个周期可能只在特定年代显著,之后被其他信号淹没。FFT把30年数据压成一条频谱曲线,看似干净,实际把最关键的信息丢了。小波变换的核心思想来自“加窗”:把一个窗口在时间轴上滑动,在每个窗口内做傅里叶变换,窗口大小还可以随频率自动伸缩。低频用宽窗口,频率分辨率高;高频用窄窗口,时间分辨率高。这就是“自适应”的时频分析,数学上叫多分辨率分析。对径流、气温这类非平稳序列,小波变换能同时告诉你“什么周期”和“何时出现”,这是FFT做不到的。1.2 小波函数怎么选:Morlet是我的默认选项选择小波母函数,是我被问得最多的问题。Matlab小波工具箱里有一堆:db系列、sym系列、mexh、morl、gaus等。做地质统计和周期分析,绝大多数文献用的是复值Morlet小波,原因有三个:Morlet是复值小波,可以同时提取振幅和相位,便于做周期显著性分析和重建信号。它的时频分辨率乘积接近不确定性原理的下限,也就是说时频分辨率搭配比较均衡。水文气象领域形成了一套约定俗成的显著性检验标准(Torrence Compo 1998那套),这套方法默认伴随Morlet小波使用,审稿人也认。要注意的是,Matlab自带的cwt函数用的是小波工具箱内置工具,默认分析中可以直接指定cwt(ts, amor)使用复Morlet小波。如果不想自己推导公式,直接从这条路径上手最稳。1.3 连续小波变换(CWT)和离散小波变换(DWT)别混淆实操中经常有人把这两者搞混。DWT(离散小波变换)一般用于信号压缩、去噪、特征提取,输出的是不同尺度下的系数序列;而周期分析用的是CWT(连续小波变换),它把尺度参数看成连续值,输出的是一个二维的时频功率谱。简单说,DWT回答“信号分解成哪些成分”,CWT回答“哪个周期在什么时间显著”。小波周期分析里核心的“小波功率谱”和“全局小波谱”,都是从CWT结果出发的。2. 实操前必须理清的几个核心概念2.1 尺度(scale)和周期(period)的换算关系这里有个新手必踩的坑:小波分析里的“尺度”不等于“周期”。对Morlet小波,当中心频率约为0.8125 Hz时,尺度与傅里叶周期近似相等(周期约等于1.033倍的尺度)。但具体换算关系取决于小波函数的中心频率设定,不同代码里出来的结果可能略有差异。实操包里的做法是直接指定周期性单位:% 设置角频率和尺度范围 s0 2; % 最小尺度,对应约2个时间单位 ds 0.25; % 尺度步长,越小精度越高,计算越慢 j1 32; % 总尺度层数 scale s0 * 2.^((0:j1)*ds); % 对数分布的尺度序列这里选择s0 2很讲究:对逐月径流序列,最小周期是2个月,低于这个值没有实际物理意义;对逐年年径流序列,最小尺度往往设到2年。尺度序列按2的指数递增,比如2、2.38、2.83、3.36……这样在低频段覆盖更多层次,符合水文数据的周期特征(常见周期如2-3年、5-8年、10-16年分布范围很宽)。2.2 小波功率谱、全局小波谱和COI区域小波功率谱(CWT功率谱):一个二维矩阵,行是时间,列是尺度/周期,颜色代表功率值,反映某个周期在某个时间的波动强度。全局小波谱(Global Wavelet Spectrum):把每个尺度上的功率对时间取平均,得到一条“平均功率谱曲线”,峰值对应的周期就是序列的主周期。COI(锥形影响区):数据两端的边界区域,因为计算时补零导致功率值失真,这个区域内的峰值不可信。Matlab绘图时用半透明阴影表示,很多初学者没意识到COI的意义,把边界上的高功率当成了重大发现,这是很常见的错误。2.3 显著性检验:别把噪声当信号小波功率高不一定是“真周期”。判断显著性,主流做法是采用红噪声(或白噪声)背景谱作为原假设,计算置信水平(通常95%)对应的功率阈值。Torrence Compo提供了标准的检验方法,具体思路是:假设时间序列由一个一阶自回归过程(AR(1))生成,计算这个背景过程的小波功率谱,如果某个时间-周期点的观测功率超过这个阈值的95%分位数,就被认定为显著周期成分。Matlab实操包里的显著性计算部分,我建议使用经大量期刊论文验证过的wt函数库(可以搜“Wavelet_coherence”等开源代码,或者参考Torrence那套代码改写的Matlab版本),而不是完全用Matlab小波工具箱内置的显著性和置信水平函数。原因是两者置信区间算法略有不同,制图时等值线画法也有差异,审稿阶段选用文献支持更广泛的方案更稳妥。3. 径流序列小波周期分析的完整实操步骤3.1 数据准备:去趋势和标准化这里用月径流数据做例子。原始数据往往带明显的年际变化和季节性,如果直接丢给CWT,低频部分的超强信号会压制其他周期的显示。实操中建议的预处理顺序是:读入序列后,先检查缺失值和异常值,缺失的月份用线性插值补齐;做季节标准化(消除月均值差异),或对序列做标准化处理(z-score);对序列做去趋势处理(去除线性趋势),否则趋势信号会污染最低频部分的功率谱。我用的预处理示例:data xlsread(monthly_runoff.xlsx); ts data(:,2); ts fillmissing(ts, linear); ts_standard (ts - mean(ts)) / std(ts); % 去趋势 t (1:length(ts)); p polyfit(t, ts_standard, 1); ts_detrend ts_standard - polyval(p, t);有人会问:去趋势会不会把长周期信号也削弱?会,所以去趋势要“适度”。如果序列本身只有年际变化而没有长期趋势,只做标准化就够了。判断方式很简单:先画一下序列图,如果有明显整体上升/下降沿,才需要去趋势。3.2 计算小波功率谱的核心代码fs 12; % 采样频率:每月一个点 dt 1/fs; % 时间步长 period 2:0.1:80; % 要展示的周期范围(月) % 用cwt或wt函数计算 [wt, f, coi] cwt(ts_detrend, fs, amor); % 从频率f转成周期 periods 1 ./ f; power abs(wt).^2;需要说明的是,Matlab新版cwt返回的频率是Hz,转换成周期时是取倒数。如果你用Torrence那套经典代码,它还自带coi输出,可以直接绘图时叠加。第一版实操包我用的是经典wt函数,因为它在显著性检验部分返回了signif阈值,方便后续画“显著区域等高线”。3.3 绘制小波功率谱图谱小波功率谱图谱是文章的核心图之一。绘制时我习惯了以下设置:figure(Color,w,Position,[100 100 1000 700]); contourf(t_years, log2(periods), power, 30) colormap(flipud(jet)) shading interp hold on % 画COI plot(t_years, log2(coi), k--, LineWidth, 2) % 画95%置信等值线 contour(t_years, log2(periods), significance_mask, [1 1], k) xlabel(时间(年)); ylabel(周期(年));这里有两个制图细节值得注意:第一,纵轴用log2刻度而不是线性刻度,因为小波尺度是对数分布的,线性刻度会把低周期部分挤成一团;第二,置信度等值线一定要叠加在功率谱上,否则读者无法区分哪些功率区域有意义。许多论文图谱里等高线就是显著性检验结果的体现。4. 结果判读:怎么看图、怎么提取主周期4.1 从图谱中读出的典型特征以长江某站月径流序列为例,图谱显示:1年周期附近有明显的功率高值带,贯穿整个时间段——这是季节变化,基本每个水文站都有,属于可预期的结果;在5~8年周期附近有一个功率峰,且在1980年前后最显著;用等高线检验后,可见该区域跨过95%显著性边界,说明这个准周期确实存在;20~30年的长周期区域功率也较高,但大部分在COI内或未通过显著性检验,只能谨慎描述为“趋势成分的体现”,不能直接判定为显著周期。4.2 用全局小波谱确定主周期全局小波谱可以避免在图谱上“凭感觉”挑周期。代码:global_wavelet mean(power, 2); % 对所有时间求平均 plot(log2(periods), global_wavelet)这时你会得到一条曲线,峰值对应的周期就是序列的“平均主要周期”。实际操作中,主峰往往对应1年;排掉1年季节峰后,找一个次峰作为“准周期”。如果次峰出现在5.3年,那么论文里可以写“该站径流序列存在显著的5~6年准周期”。记得在图上画一条显著性水平的横线,给出客观判断依据。4.3 分频段重构波形有些情况下,单看功率谱还不够。比如你发现8~16年周期区域有信号,你想要这个周期带的时间形态——什么时候偏丰,什么时候偏枯。这时要利用小波逆变换,把目标周期带对应的系数重构出来。Matlab中可以通过icwt实现:% 重建8~16年周期分量的信号 period_band [8 16]; ts_recon icwt(wt, [], period_band);重建后的曲线叠加在原序列上,可以直观看到该周期带的相位变化,比如“1950-1970年间处于偏丰期,1970-1990年转为偏枯期”。这种做法对分析径流丰枯阶段性非常有价值。5. 气温序列与径流序列的分析差异5.1 气温序列的显著性特征气温序列和径流序列的小波分析看着一样,实际差异很大。气温序列通常平滑、噪声少、长周期信号更强;径流序列受降水、下垫面、人类活动影响,噪声大,短周期成分复杂。因此:气温序列往往更容易出现显著的“全球变暖长期趋势”信号,图谱中最低频部分功率很高,而且覆盖COI外的大片区域;径流序列在1~2年的高频分量明显更活跃,特别是受ENSO影响的地区,2~7年周期带显著;气温序列的显著性检验通常可用白噪声背景谱,径流序列优先用红噪声背景谱,因为径流本身具有较强的一阶自相关特性。5.2 季节标准化对气温序列的影响对气温序列做季节标准化时要非常小心。气温的月均值本身就随季节大幅波动,如果粗暴地做z-score,会保留季节分量;如果去掉季节循环(比如减去多年月均值),那分析的焦点就集中到了“异常气温”的准周期上,这种做法的结果更适合讨论气候异常事件。两种方式对应不同的科学问题,先问自己“要分析的是原始气温的周期性,还是气温异常的周期性”,再做预处理。这个决定直接影响结果能说明什么问题,不可不慎重。5.3 多站点对比时的统一参数做完单站后,如果要做流域内多站点的周期对比,最常犯的错误是每站用不同的尺度范围或显著性水平。实操包中做多站对比时,建议所有站点统一使用相同的尺度序列、相同的背景谱假设、相同的显著性阈值(如95%),否则站点间的差异可能只是参数不同带来的假象。我遇到过有人用5个站点分别跑不同的尺度范围,最后图上各站主周期差异很大,实际上只是尺度采样粗细不同导致的视觉偏差。6. 制图排查与常见报错处理这里把实操包中反复出现的几个问题和排查思路列出来,这些都是文档里不太会写的东西。6.1 图谱颜色一片红或一片蓝导致这个现象的原因通常是功率值跨度太大,色标被少数极端值拉爆。对策是先对功率取对数:power_log log10(power),再用contourf绘制。这样中低功率区域的差异才能体现出来,不然整张图就会变成几块纯色。6.2 cwt函数提示Data lengthening错误部分旧版本Matlab对边界进行处理时要求序列长度满足特定条件,而cwt不支持NaN,如果数据里有NaN,会报错。解决方案是提前用fillmissing处理缺失值,或者用rmmissing删掉不连续的段(但删除会改变时间轴,不建议)。更稳妥的是保留原始时间轴,用插值补齐。6.3 COI画不出来使用Matlab内置cwt时,输出参数没有直接给出COI,需要通过cwtfreqbounds或coi变量处理。如果发现coi返回的是频率单位,记得转换成周期再绘图。我用过几个版本的Matlab,coi的返回格式确实不太一致,打印出来先看一眼单位再画图,能省掉很多麻烦。6.4 纵轴周期刻度标注混乱绘图时如果直接指定YTick和YTickLabel,容易因log2坐标而产生错位。推荐用如下写法,干净且不易出错:yticks [1 2 4 8 16 32 64]; set(gca, YTick, log2(yticks), YTickLabel, yticks);配合这套刻度,图上显示的周期值一目了然,不会出现坐标轴标注和实际位置对不上的情况。6.5 多图排版输出做论文时,小波功率谱、全局小波谱和重建序列通常要拼在一起。用subplot时注意调整图幅大小,建议输出宽度在17cm以上(双栏期刊的半栏宽度约8.5cm,整页约17cm),并用exportgraphics导出tiff格式,300dpi以上,这样审稿阶段不会因为清晰度被挑刺。7. 配套视频里的分步操作逻辑实操包里附带的分步操作视频,总时长约40分钟,我刻意按照“数据准备→参数设置→图谱绘制→显著性检验→结果判读”的顺序录制,因为这是我日常带学生做项目时最常见的完整流程。视频里用同一套径流数据从头跑到尾,保证操作和本文代码完全对应。视频中多次强调的一个检查动作:每算完一步,都先把变量工作区的维度检查一遍。power矩阵的行数应该等于periods的长度,列数等于时间序列长度。很多学员代码报错,原因都在于矩阵维度对不上——比如尺度序列是从0开始编号,而功率矩阵从1开始编号,最终绘图时错位。这个检查习惯了之后,能避开至少一半的bug。另外视频里也演示了怎么把小波谱图和传统FFT谱图放到一起做对比。这种对比展示在论文里非常有用:FFT给出的是“平均峰”,小波谱给出的是“峰何时显著”,两者配合,审稿人会觉得你对该方法有完整理解,而不是只会套软件。8. 经验总结与后续扩展建议我用Matlab做小波周期分析几年下来,最大的体会有三点。第一,预处理比分析本身更容易影响结果,一个没做去趋势、没处理缺失值的序列,跑出的功率谱可能完全是误导性的,所以不要在数据清洗上省时间。第二,显著性检验是周期分析的“底线”,要么用Torrence的经典检验方案,要么选用成熟的工具箱函数,不建议自己写一个简化的随机检验然后宣称“显著”。第三,小波功率谱图的配色和标注虽然看起来只是细枝末节,但期刊审稿阶段,杂乱无章的图谱往往会被质疑整个分析的严谨性,值得花时间调好。后续想继续扩展的话,可以从两个方向入手。一是把单序列的CWT扩展到交叉小波变换和相干小波分析,用于分析两个序列(如降水与径流、Nino3.4指数与径流)在时频域的共振关系,这也是当前水文气候研究的热门内容;二是把Matlab的流程迁移到Python,用PyWavelets或cwt结合SciPy实现类似功能,方便后续做批量化和自动化分析。不过对大多数只需要周期性结论的科研场景,Matlab这一套流程已经足够完善。最后再分享一个小技巧:保存图谱时,记得同时保存一份fig格式的源文件。审稿人提出“把色标改成别的配色”或“把周期范围放宽”这类意见时,有源文件几分钟就能改出来,没源文件就得重跑一遍画图代码,那种返工的感觉确实不太好受。本文还有配套的精品资源点击获取
