降雨数据这东西看起来就是一串数字但里面藏着的规律远比表面上多得多。我第一次拿到60年逐月降雨记录时最想知道的问题只有一个这些年的旱涝变化到底有没有周期性哪个年代干旱主导哪个时期湿润信号突然增强起初我用傅里叶变换跑了一圈全局周期倒是能看出来可那些周期具体是什么时候出现的、中间又消失多久完全看不出来。后来换成连续小波变换CWT在Matlab里把降雨时间序列拆成“时间-尺度-能量”三个维度问题一下就清楚了。这篇文章完整记录我用小波变换分析降雨数据的全过程包括原理、代码、绘图和周期显著性检验的实操要点适合正在做水文数据分析、气候统计或者相关课程设计的同学直接参考。1. 为什么用连续小波变换分析降雨序列1.1 傅里叶变换的天然局限先说傅里叶变换。FFT做出来的功率谱确实能告诉你整个序列里有哪些频率成分比如在周期3年和周期12年处出现明显峰值。但问题是这个结果是“全时段平均”的。也就是说如果某个2-4年的准周期只在1980年代出现、到2000年以后彻底消失FFT仍然会给一个“平均意义上存在2-4年周期”的结论至于这个周期是什么时候出现的、何时中断的谱图上一律不体现。这就像你听一整场音乐会FFT能告诉你曲子里有钢琴、有小提琴、有鼓但完全说不出钢琴声是在第几分钟出现的。降雨序列恰恰是典型的非平稳信号旱涝模式的持续时间会随年代变化气候系统内部的阶段性特征非常明显。所以做降雨周期分析如果只用FFT等于把时间维度的信息全部丢掉了。小波变换解决的就是这个问题。它在时间和频率两个维度上同时展开信号输出是一张二维谱图横轴是时间纵轴是尺度或周期颜色深浅代表能量强弱。这样你既能看出所有周期又能知道这些周期在哪些时段显著。1.2 小波变换的工作原理连续小波变换的数学形式是W(a,b) 1 / sqrt(a) ∫ x(t) ψ*((t-b)/a) dt其中a是尺度参数对应周期大小b是平移参数对应时间位置ψ是小波基函数。核心思想很简单拿一个小波基函数在信号上平移改变b同时伸缩改变a不断计算信号与基函数的相似程度。尺度a小的时候基函数被压缩能捕捉高频短周期的细节尺度a大的时候基函数被拉伸能捕捉低频长周期的整体趋势。这个过程相当于一个焦距自动变化的显微镜。低频成分用宽窗看分辨率低但看得长久高频成分用窄窗看分辨率高但看得局促。相比之下FFT的窗是固定的所以没法对每个频率都做到最优的时间定位。在降雨和气象数据周期分析里最常用的基函数是复Morlet小波。Matlab中名为amor通常取中心频率ω0≈6。它是复值小波能同时提取振幅和相位而且时频分辨率组合比较平衡适合用来探测振荡型信号。如果你用实值小波比如mexh相位信息就没了做周期显著性检验也略麻烦。所以一般降雨分析直接用amor最省事。1.3 为什么我坚定选Matlab做小波分析的工具不少Python有PyWaveletsR也有相关包。但我自己最顺手的还是Matlab尤其是从R2016b开始cwt函数的封装已经非常成熟一行就能完成连续小波变换并返回频率向量和锥形影响区COI。不需要手动构造尺度序列也不需要在各个库之间调试版本冲突。Matlab还有一个优势是绘图生态。小波功率谱这种东西需要tiledlayout控制多子图、对Y轴做对数刻度、叠加colorbar、画显著性等值线Matlab里几个命令就能组合得很干净。另外如果你要做双变量分析自带的wcoherence、xwt函数可以直接生成交叉小波谱和相干谱省掉大量自写代码的时间。下面这个表是我给新手解释小波和FFT区别时经常贴出来的对比先存个印象后面实操会反复用到对比项FFTCWT时间定位无有频率分辨率固定高频低、低频高适合信号平稳信号非平稳信号输出形式频率-功率曲线时间-周期-功率二维图边界问题隐含周期假设有COI需注意降雨适用场景初步查看整体周期定位各时段周期变化2. 数据准备与预处理2.1 拿到降雨数据后的第一件事理论上小波变换可以直接套在一串原始数字上但实操中如果你不好好处理数据后面谱图上出现鬼影一样的假周期排查起来特别头疼。所以我每次拿到降雨数据第一步都是做数据体检。先确认时间分辨率和长度。逐月数据和逐日数据能分析的周期范围完全不同。我下面以逐月数据为例做演示采样间隔是1个月。如果数据是逐日的你就要把dt改成1天分析与历法无关的周期。再做缺失值检查。很多气象站点数据或多或少有缺测尤其连续缺测超过几个月的情况。直接用fillmissing线性插值是最快的方案但有一个前提缺测段不能太长。缺1-2个月线性插值问题不大缺了一整年插值出来的序列在那个时段没有任何真实波动会在小波谱里产生低能量伪特征。最后是单位与量纲。降雨量单位尽量统一成毫米不要一会儿英寸一会儿毫米。这虽然不影响小波分析的最终结论但提示数据源的一致性尤其当你拿多站数据做对比的时候量纲不齐就会让功率谱的可比性变差。% 读取示例降雨数据 rainData readtable(rainfall_monthly.xlsx); t rainData.date; % datetime 类型 rain rainData.precip; % 单位mm % 缺失值线性插值 rain fillmissing(rain, linear);2.2 标准化和去趋势要慎重预处理里最容易忽略的是标准化。小波功率谱的绝对值与信号方差成正比如果你拿两个不同站点对比方差大的站点谱图颜色整体都会偏深这不代表它周期更强只代表它波动更剧烈。为了可比性我习惯先把序列标准化rain (rain - mean(rain)) / std(rain);标准化之后序列均值是0方差是1小波功率谱的能量值就有了相对统一的量纲可以跨站点看相对强度。去趋势就没那么标准化了。降雨序列在很多地区并没有明显的长期单调趋势但有的站点会呈现缓慢增加或减少。这个趋势会污染小波谱的低频端让长周期出现大块高能量区而这其实是趋势导致的假象不是真实周期。我的做法是先做线性拟合看斜率是否显著。显著就去掉线性趋势不显著就保留。用detrend(rain)把线性趋势抽掉但注意这只会去掉线性成分如果是非线性缓慢变化可以用低通滤波大致分离后再提取周期特征不过这样容易过头把真正的年代际信号也滤掉。宁可只去线性趋势不去做过度滤波。2.3 采样间隔决定周期下限采样定理在小波里同样成立。如果数据间隔是1个月那么你可以分析的最短周期是2个月小于这个基本没有意义。一般实际分析会建议最小周期为2-3个采样间隔也就是6到8个月。如果你做逐月数据却看到3个月的显著能量要么是原始数据本身就含这种高频信号要么就是插值或者噪声搞出来的伪谱需要多留一个心眼。最大周期则受序列长度限制。经验上说最长可分析的周期大约为序列长度的1/4到1/2。60年逐月数据全长720个月理论上可以分析到180-360个月的周期也就是15到30年。但这样的长周期会在数据两端严重受到边界效应影响真正可信的中间区域很小。所以大部分降水研究会把周期分析范围控制在2到20年既有气候物理意义又避免边界带来的大范围污染。3. Matlab小波分析实操全流程3.1 核心函数选型Matlab的小波工具箱里有多个相关函数初学者容易搞混。cwt是连续小波变换的主入口现在的版本支持直接指定小波和采样周期返回小波系数、频率和COI。官方文档推荐用这个做连续小波时频分析。cwtft是旧接口用来做傅里叶变换基础上的CWT灵活性更高但参数复杂一些很多老教程还在用它。除非你要精细控制尺度序列否则没必要碰。cwtfilterbank用于构造滤波器组一般用来做信号分解、重构或者自定义时频图不是做降雨周期分析的首选。dwt和wavedec是离散小波变换做去噪、压缩、特征提取才用别混进来。我看到有人拿小波去噪的细节层直接说成“周期”这是概念错位后面在常见问题里我会专门提。所以我这篇博文用的就是最普通的cwt(rain, amor, dt)一行出结果。3.2 参数设置小波基、频率范围和尺度数小波基上面说了用amor即复Morlet。关于voices参数它控制每倍频程内的尺度数。默认10对多数情况够了提高到16可以让谱图更平滑代价是计算时间增加和低频端边界变宽。除非你某个特定周期带的边缘特别毛糙不然默认就足够。频率范围方面cwt默认会根据信号长度自动选择最大尺度但你真正关心的可能是某个区间。我不建议把计算范围人为缩小直接全频算绘图时再用ylim限定显示范围就好。这样后面想再检查更高或更低频段不用重新跑一遍变换。采样周期参数dt设置成1代表单位是月。此时返回值f的单位是cycles/month周期等于1./f单位是月。如果我处理的是日数据dt就设为1天周期就是天。千万别把dt理解成“序号间隔”它的单位直接影响整个频率轴的解释。3.3 完整代码从读取到绘制小波功率谱下面是一段可以直接运行的完整流程以60年月降雨为例%% 读取与预处理 rainData readtable(rainfall_monthly.xlsx); t rainData.date; rain rainData.precip; rain fillmissing(rain, linear); rain (rain - mean(rain)) / std(rain); N length(rain); dt 1; % 月 fs 1 / dt; % 每月1个样本 %% 连续小波变换 [wt, f, coi] cwt(rain, amor, fs); power abs(wt).^2; period 1 ./ f; % 周期单位月 %% 绘图 t_years year(t) (month(t) - 1) / 12; figure; tiledlayout(2, 1, TileSpacing, compact); nexttile; plot(t_years, rain, LineWidth, 1); xlabel(年份); ylabel(标准化降雨); title(月降雨时间序列); grid on; nexttile; imagesc(t_years, period, power); set(gca, YDir, normal, YScale, log, YLim, [2 120]); xlabel(年份); ylabel(周期月); title(小波功率谱); colorbar;这段代码的关键点在于imagesc的Y轴设置。period是从高频到低频还是低频到高频排列取决于cwt输出的频率向量顺序。用YDir控制方向再用log让周期轴在低周期段拉开区分度不然2个月和12个月会挤在底部看不清。如果你要对不同子图统一比较还可以用clim([0, 3])或者设置统一的颜色范围。后面常见问题里我会讲怎么调颜色范围避免谱图“一片蓝”。3.4 全局小波功率谱与显著性检验小波功率谱已经告诉你每个时段的周期能量了但如果要判断哪些周期是统计显著的而不是噪声偶尔引发的波动就得做显著性检验。有两个常用路线我分别说一下。第一是理论分布法。Torrence和Compo在1998年那篇经典论文里指出在红噪声背景下小波功率谱的局部值服从χ²分布对于Morlet小波来说自由度近似为2。因此可以用AR(1)红噪声模型估计背景谱再乘上对应置信度的χ²临界值就能得到95%显著性线。这个做法在数学上很漂亮但初学时容易把公式的系数弄混。我更推荐蒙特卡洛方法思路直接也不容易错生成大量与原始序列一阶自相关特征相同的AR(1)随机序列对每个序列都做CWT计算全局平均功率谱取95%分位数作为显著性阈值。真实数据的全局功率谱如果超出这条线就认为这个周期上的能量强度显著高于随机背景。% 估计AR(1)参数 alpha corr(rain(1:end-1), rain(2:end)); sigma std(rain); % 生成单个AR(1)序列 ar1 (n) filter(1, [1 -alpha], sigma * randn(n, 1)); rng(2024); ns 300; surrPower zeros(length(f), ns); for i 1:ns r ar1(N); r (r - mean(r)) / std(r); [wt_s, ~, ~] cwt(r, amor, fs); surrPower(:, i) mean(abs(wt_s).^2, 2); end sig95 quantile(surrPower, 0.95, 2); % 全局小波功率谱 globalPower mean(power, 2); figure; semilogx(period, globalPower, LineWidth, 1.5); hold on; semilogx(period, sig95, --, LineWidth, 1); set(gca, XDir, reverse); % 让长周期在左还是右自己按习惯设定 xlabel(周期月); ylabel(全局功率); legend(全局小波谱, 95%显著性阈值, Location, best);蒙特卡洛的优点是稳健适合所有数据形态。缺点是慢300次AR(1)生成和CWT在普通电脑上大概要跑几十秒可以接受。如果你对频率分量的数量敏感可以适当减少到100次但一般建议不少于200次。注意AR(1)序列要同样做标准化否则背景功率尺度对不上。4. 结果解读与应用场景4.1 一张小波功率谱怎么读拿到谱图以后很多人盯着五彩斑斓的图发愁不知道哪些区域能信。我有一套固定的读图顺序。先看横轴时间、纵轴周期、颜色能量。颜色越浅能量越高。一般用parula或jet配色高能量区域在浅色调位置。接着找能量集中带。稳定周期往往表现为一条连续的高能亮带比如12个月处全年都很亮说明年周期稳定存在。如果亮带在某个时间段中断说明那个时段该周期变弱。再看周期带的宽度。能量带越宽代表周期越不精确可能是多个相近尺度混合。这是正常现象降雨序列很少像精密仪器一样输出一个绝对精确的周期。最后一定要看COI。cwt返回的coi参数标出了数据边界影响区。边界附近的计算因为小波部分落在信号之外能量被低估所以处于COI之内的能量带不能解读成真实周期。这是所有小波分析新手最容易犯的错看到左右两边一片亮就觉得有超长周期其实那只是边界效应。4.2 识别降雨的多尺度周期结构我拿自己处理过的一个站点数据举例。60年逐月降雨序列做完CWT后谱图上大概能分出三层结构。年周期亮带贯穿始终位置稳定这是因为站点受季节降水控制每年雨季都有稳定表现。其次在3-7年周期带上出现了一个分段明显的亮块集中在1980到1998年之间前后时段几乎看不到高能特征。这个就说明该站点在那些年份存在准周期性的湿润增强信号可能与当时降水带的年代际摆动有关。而这种时段性特征FFT根本表达不出来。再往长周期看15年以上区域在左右两端几乎全是COI覆盖中间只有一小块可信区域所以我不会对它下太强的结论。做小波分析最重要的一课就是学会在没必要下结论的地方守口如瓶。4.3 扩展到多站点交叉小波分析如果你手头有两个站点的降雨序列想看看它们之间的周期关联或者想对比降雨与某个气候指数之间的关系单变量CWT就不够用了。这时候用交叉小波变换XWT和交叉小波相干WTC。Matlab里wcoherence(x, y, fs)可以直接返回相干谱能量高的区域表示两个序列在对应时频区域有较强共同振荡。另外相位箭头的方向还能告诉你两个序列之间的滞后关系箭头水平朝右同相竖直朝上领先90度。这个在地学里特别常用比如看降雨对气候响应的滞后时间。我做多站对比的时候一般先把每站标准化然后计算两两之间的wcoherence。如果某个周期带上相干谱出现大面积高亮说明这些站点共享同一个气候驱动因素对该区域水资源管理有直接参考价值。5. 常见问题与排查技巧实录5.1 边界效应COI是底线不是装饰边界效应是CWT的固有缺陷。小波窗口的一部分掉出信号边缘时数据量不够计算出的系数不可靠。coi就是用来标记这个不可信区域的。如果你的序列只有30年却想分析20年周期那整个谱图几乎有一半都在COI里结论基本没法用。遇到这种情况我的经验有两个方向一是换数据拉长约40年以上二是把分析目标压缩到短周期段不去硬碰长周期。千万不要用“两端补零再分析”来假装延长数据。补零会让边界处多出一堆人为的低频能量谱图看起来更“饱满”但那全是假的。补数据不如不补。5.2 谱图一片蓝怎么办很多新手第一次跑出来的小波功率谱图整体颜色都很暗高能区域只有一小撮。这通常不是计算错误而是颜色标尺范围被极端值拉大了。降雨序列中偶尔会出现某个月份功率极高比如一次异常极端的降水事件它的能量会把colorbar的最大值顶得很高其他正常信号相比之下都显得很暗。解决办法是限制颜色显示范围。用clim(prctile(power(:), [0 95]))之类的方式把95%分位数作为上限丢弃那一片孤立的极端能量让主要能量区的颜色层次拉开。我几乎每次出图都会调整这一步默认的clim很少能直接达到出版级效果。不然就对功率取对数再画。imagesc显示log10(power)也是一种常用做法它的好处是让强弱能量之间的差异变小适合动态范围大的数据。需要说明的是如果做显著性检验必须用原始功率值只在显示环节变换。5.3 小波基选不同结果完全对不上有人拿db4做离散小波去噪有人拿mexh做CWT两种结果放在一起比较周期特征很容易得出互相矛盾的东西。这不是软件坏了而是两种小波分析根本不是一个层面的工具。db4这类离散小波适合做多分辨率分解输出是细节系数和近似系数用于去噪、压缩效果好。但它的频带划分固定不像Morlet小波那样能灵活地连续扫描频率。如果你想研究降雨的连续周期变化请坚持用连续小波特别是复Morlet。同样在连续小波里mexh是实值二阶高斯小波对单纯峰值位置敏感适合定位“事件发生时刻”比如洪水过程峰值但对周期振荡的检测不如复Morlet。做降雨周期工作首选amor基本不会翻车。5.4 频率与周期换算别搞反Matlab的cwt返回的f的单位是cycles/sample或者乘以采样率后的Hz很多人直接拿f当周期用结果图的横轴完全乱掉。我已经强调过周期 1 / f。月数据下f的单位是cycles/month周期单位是月。日数据下就是天。如果period序列方向是从高频到低频画图时ylim就会比较奇怪。我的建议是绘图前先disp(period(1:5))看一眼确认方向再用set(gca, YDir, normal)调整显示方向。这个细节能省去大半绘制时的懊恼。5.5 问题排查速查表常见现象可能原因排查与处理谱图两边出现高亮大块区域边界效应查看COI忽略锥内区域全图颜色过暗colorbar范围被极端值拉爆使用clim限制到95%分位12个月亮带但中心偏到13个月数据时间段或插值引入偏移检查原始数据是否有闰月、缺测处理结果出现小于2个月的周期采样频率不足或插值噪声压低显示周期下限确认数据质量周期带一片模糊看不出结构voices数太低调高voices到16重算与文献周期结论对不上小波基或参数不一致确认是否都用复Morlet和类似频率范围COI信息不会画对坐标转换不熟在imagesc绘图后手动叠加coi曲线提示做小波分析之前先花10分钟把数据的基本统计量、趋势和缺测情况看一遍。数据质量决定了谱图质量这一步省不掉。最后分享一下我的个人体会前前后后跑了几十组降雨数据之后我的感觉是小波变换是最适合用来“看”降雨多尺度结构的手段但它不是什么灵丹妙药。它不会自动告诉你这个周期是物理机制还是随机波动显著性检验只能排除一部分随机性真正让结论站稳脚跟的还是你在实际气候背景下的解释能力。还有个细节Matlab不同版本之间cwt的语法有过调整早期用cwt(x, scales, wavename)现在用cwt(x, wavename, fs)。如果你把网上的老代码直接贴进新版本大概率会报错。遇到这种情况不用慌看一眼当前版本的帮助文档把参数顺序调整过来即可。如果后续想把这个项目扩展可以往两个方向走一是把所有站点的数据整理成一个大表格批量计算功率谱后聚类出几个典型周期区域二是把CWT结果作为特征和季节预测模型结合起来看能不能从中提取对干旱预测有用的信号。但无论哪个方向小波分析这一步打下的底子都会让你在处理非平稳降雨序列时心里有个清清楚楚的谱图。