做空气质量数据分析那会儿我最常用的是把两个变量的相关系数算一算比如PM2.5和风速的负相关、PM2.5和相对湿度的正相关然后写进报告里。后来发现一个尴尬的问题同一组数据用全年数据算相关和用冬季数据算相关结论经常对不上甚至符号都能反过来。原因倒不难理解气象和污染物序列是非平稳的相关性本身就在随时间尺度和时间段变化。这时候“基于Matlab的二维小波相干分析”就派上了用场——它能在时间和频率两个维度上同时展示两个序列的关系强度告诉你这两个变量在哪个时间段、哪个周期尺度上显著相关、相位是正还是反。这篇博文我会用空气质量数据当例子把Matlab里小波相干分析的原理、工具选型、完整实现流程和常见坑一次讲透适合刚接触小波分析、或者已经会用相关分析但对时频方法不熟的朋友参考。1. 方法原理与适用场景为什么单算一个相关系数不够1.1 普通相关性分析在气象环境数据里的痛点先从一个简单的例子说起。假设手上有一整年的PM2.5小时浓度序列和同时段风速序列用Pearson相关去算大概率会得到显著负相关。这个结论看似合理风吹散了污染物。但如果你把数据切成不同季节看冬天静稳天气下风速较小的时段往往也是污染物堆积最快的时段风速和PM2.5反而可能呈现弱正相关夏天对流强暴雨前后风速和污染物浓度的关系又完全是另一套模式。普通相关分析输出的只是一个全局数值它假设两个变量的关系是平稳、线性的。但空气质量数据从来不是这样排放源随时间波动气象条件随季节和天气过程变化污染物浓度既有日尺度上的早晚高峰规律也有周尺度上的工作日-周末差异还有月、季尺度上的冷暖交替影响。关系是分层的、动态的。如果只用一条相关系数去概括必然丢掉大量信息。另一个痛点是非平稳性。很多时间序列方法比如简单回归、傅里叶变换都要求序列平稳或者至少均值方差稳定。可污染物浓度和气象要素几乎都带趋势和明显的季节循环直接套用传统频域分析很容易失真。小波分析的优势就在这里它不要求整段序列平稳而是把“整体相关”拆解成“局部时频相关”相当于给分析加了一个可以滑动的时间窗口同时又保留了频率信息。1.2 从交叉小波到小波相干的基本逻辑要理解小波相干得先知道连续小波变换CWT在干什么。CWT把一维的时间信号用一个小波基函数最常用的是Morlet小波在不同尺度上做卷积输出一个二维的系数矩阵横轴是时间纵轴是尺度通常换算成周期颜色代表信号在该时频位置的功率强度。这个矩阵可以理解为把信号的“能量”铺在了时间-频率平面上。对两个序列分别做CWT后就可以计算交叉小波谱cross wavelet spectrum。交叉小波谱衡量的是两个序列在某个时间和周期上共同出现高能量区域的强度。但有个问题如果一个序列本身功率就特别大交叉小波谱会被它的能量牵着走看起来像两者强相关其实只是单边强。所以实际分析中更常用的是小波相干wavelet coherence。它的计算思路是先对两个CWT谱做二维平滑再对交叉谱做同样的平滑然后用平滑后的交叉谱振幅平方除以两个平滑后功率谱的乘积。数学上看有点类似于频域里的决定系数结果被归一化到0到1之间。值越高说明在某个时频位置上两个序列的局域相关越强。比起交叉小波谱小波相干更不容易被单个变量的强能量干扰对环境数据这种“方差随时变”的序列更友好。这里需要说明一下“二维”的含义。在Matlab和大多数文献语境里小波相干分析的输出本身就是一个二维矩阵横轴时间、纵轴周期颜色代表相干强度。有些场景还会把两个空间维度都做小波变换比如用二维小波多分辨分析处理污染物浓度空间场那是另一个方向。本文聚焦的还是最常用的“两条一维时间序列在时频二维平面上的相干分析”这也是目前空气质量与气象因素关系研究中用得最多的方法。1.3 结果图的读法和关注要点小波相干图几乎成了这类分析的标配输出。横轴是时间纵轴建议直接用周期比如小时、天避免用频率给人造成直觉上的别扭。颜色深浅代表相干系数从蓝到红表示相关性从低到高。图上通常还有黑色等值线表示通过显著性检验的区域比如95%置信水平。相位信息也很关键通常用箭头表示。箭头指向右侧表示两个序列在该时频点同相即同时增大或同时减小指向左侧表示反相一个增大另一个减小箭头略微向上或向下则表示两个序列之间存在相位超前或滞后关系。比如PM2.5与风速在小尺度小时级别上常表现为反相箭头多偏向左侧而PM2.5与相对湿度在某些周期上可能有滞后箭头会明显向上或向下倾斜。刚看到小波相干图时最容易犯的错误是只看颜色深浅忽略显著性区域。颜色再红如果没有通过显著性检验多半是噪声。我的习惯是先锁定有黑色等值线的区域再看其中的周期范围和箭头方向最后回到原始时间序列上做局部验证这样才能形成一个完整的证据链。边界附近的锥形区域COI则要格外小心那里受边界效应影响大尽量不要下结论。2. 数据准备与Matlab工作环境搭建2.1 空气质量数据该怎么准备在做小波相干之前数据清洗往往决定结果质量。很多人拿到的原始数据是监测站点提供的逐小时平均值常见问题是时间戳对不齐。比如PM2.5用了某省站点的数据风速用的却是机场气象站的数据两个序列的采样时间不一致或者中间缺了几个小时这时候直接塞进函数里轻则报错重则得到一条离谱的相干谱。我通常按四步处理。第一步把时间统一成datetime类型并按小时对齐第二步删除或插补缺失值如果缺失比例低于5%用线性插值或者前一天同时刻平均值补上超过10%就建议该时段直接截断第三步剔除明显异常的极值比如PM2.5出现负值、风速出现几十米每秒的野值等最好结合原始记录的质控标记处理第四步决定是否标准化。关于标准化有一点容易困惑。小波相干对两个序列的绝对量级不敏感因为最后会做归一化处理所以你就算一个序列是100微克每立方米的PM2.5、另一个是0到10米每秒的风速也不影响相干数值。但标准化会影响相位箭头的解释吗不会相位反映的是变化方向。所以是否标准化更多是个人选择我喜欢保留原始量纲去做分析解读结果时更贴近物理量有时为了画图方便也会用z-score标准化后放在同一张图上对比这时结果中的相干数值不变。2.2 Matlab工具箱选型内置函数和经典工具包怎么选这里要给新手提个醒别急着下载各种来路不明的工具箱先看清楚Matlab自带的Wavelet Toolbox能做什么。从R2016b之后MathWorks在Wavelet Toolbox里加入了wcoherence函数可以直接计算两个信号的小波相干并返回谱矩阵和频率向量。对绝大多数处理空气质量数据的需求这个函数完全够用。如果你做研究需要复现论文里那种风格的小波相干图或者要调整显著性检验的蒙特卡洛次数、想深度定制绘图样式那Grinsted等人发布的开源工具包俗称wtc工具箱就更有用。这套工具包早期用于地球物理研究思路清晰很多环境、气候、经济领域的论文都引用了它。它提供了wtc函数计算小波相干和wtcplot函数绘制带显著性检验和相位箭头的图还自带基于AR(1)红噪声的蒙特卡洛显著性检验。我做了个简单的对比对比项Matlab内置wcoherenceGrinsted工具包wtc上手难度低参数少直接调用即可中需要理解几个尺度参数输出完整性返回相干谱、交叉谱、频率返回完整结构体含显著性结果绘图风格简洁原生风格经典论文风格可定制性强显著性检验内置蒙特卡洛内置AR(1)红噪声蒙特卡洛批次处理遍历调用方便同样方便但参数要多写几行适合场景快速分析、内部报告期刊论文、复现经典方法如果只是自己分析数据优先用内置函数如果要写论文且希望图和审稿人预期一致建议用Grinsted工具包。不过无论选哪种都要记住核心步骤是一样的算CWT、算平滑交叉谱、算相干谱、做显著性检验。2.3 文件组织与运行准备我习惯把这类分析项目组织成固定结构。根目录下放数据文件夹data、脚本文件夹scripts、结果输出文件夹output再加一个README备忘。脚本方面我一般写四个脚本01_load_data.m负责读取和清洗数据02_preprocess.m统一时间轴和缺失值03_wavelet_coherence.m做核心计算04_plot_results.m负责出图和导出这样每个环节出了问题都好排查。如果你要使用开源的wtc工具箱下载后把整个文件夹放到项目目录下然后在脚本开头用addpath(genpath(path/to/wtc_toolbox))把文件加进路径。注意不要在Matlab默认路径下乱丢文件否则升级Matlab之后容易出兼容问题。这一步省不了因为调用函数时Matlab是按路径搜索的路径没配好就会报“Undefined function or variable”。3. 核心实现用Matlab跑一遍二维小波相干分析3.1 用内置wcoherence函数快速上手我先讲最省力的方式。假设数据已经清洗好存在一个表格里至少有三列时间t、pm25、wind_speed时间列是datetime格式逐小时数据。用内置函数做小波相干核心代码非常短data readtable(air_quality_hourly.csv); t data.time; x data.pm25; y data.wind_speed; % 剔除缺失/无效数据简单做法任一方为NaN的行直接去掉 valid ~(isnan(x) | isnan(y)); t t(valid); x x(valid); y y(valid); % 指定采样间隔为1小时 dt hours(1); % 调用wcoherence [wcoh, wcs, freq] wcoherence(x, y, dt); % 画图 figure imagesc(t, freq, wcoh); axis xy; set(gca, YScale, log); colorbar; xlabel(时间); ylabel(周期/频率); title(PM2.5与风速的小波相干);wcoherence函数会自己选择合适的Morlet小波参数返回三个输出wcoh是二维相干矩阵wcs是交叉小波谱freq是频率向量。需要注意freq是频率而不是周期画图时如果你想让纵轴展示“周期”要对freq取倒数并调整坐标方向period 1 ./ freq; imagesc(t, period, wcoh); axis xy; set(gca, YScale, log);这一步我经常忘所以单独拿出来说。刚接触时容易直接用freq画图结果纵轴是“频率”看图的人还得换算成周期增加沟通成本。空气质量研究里大家更关心“几天一个周期”所以纵轴用周期更直观。wcoherence还支持一些可选参数比如FrequencyLimits和NumOctaves控制频段范围和倍频程数。默认情况下结果覆盖整个可用频段但数据长度较短时低频段几乎都被COI覆盖画出来没意义。我的做法是提前估算数据长度支持的最大周期只展示有意义的那一段避免图被无效区域占满。3.2 用Grinsted工具包实现自定义风格如果你想走论文风格那就用wtc工具包。以同样数据为例代码大概是这样的data readtable(air_quality_hourly.csv); t data.time; x data.pm25; y data.wind_speed; valid ~(isnan(x) | isnan(y)); t t(valid); x x(valid); y y(valid); % 预估AR(1)系数用于置信检验 Lag1 ar1(x); % 设置小波尺度参数 S0 1/4; % 最小尺度这里取0.25小时也可设成2个采样间隔 Dj 1/12; % 尺度序列间隔值越小计算越慢 J1 2*round(log2(length(x)/S0)/Dj); % 最大尺度指数 pad 1; % 补零提升计算效率 MonteCarloCount 200; % 显著性检验重复次数 % 计算小波相干 [WCO, Coi, WT, Coherency] wtc(x, y, Pad, pad, Dj, Dj, S0, S0, J1, J1, ... Lag1, Lag1, MonteCarloCount, MonteCarloCount); % 绘图 figure; wtcplot([WCO, Coi, WT, Coherency], xdate, t); title(PM2.5与风速的小波相干wtc工具箱);这里有几个参数值得解释。S0是最小尺度通常取两个采样间隔内的值比如小时数据可以设成1/4小时也可以设成2小时的倒数形式。Dj控制尺度轴的粒度常用1/12数值越小计算越慢但尺度分辨率越细。J1是尺度级数的上限工具包作者给出的经验公式是2*round(log2(N/S0)/Dj)它会根据数据长度自动计算最大尺度。MonteCarloCount是显著性检验的重采样次数默认通常是100到300分析论文级别的数据时可以设成300耗时增加但结果更稳定。wtc函数返回的变量顺序和含义不同版本略有差异。我上面写的四输出形式[WCO, Coi, WT, Coherency]里WCO是包含相干值、显著性水平等信息的结构体Coi是边界锥形区域的坐标WT是交叉小波谱相关的结构体Coherency是小波相干矩阵。直接用wtcplot函数可以一键出图图上会自动画黑色显著性等值线和箭头方向非常省事。3.3 批量处理多组序列从单对扩展到多要素组合实际项目里不可能只分析一对变量比如要比较PM2.5与风速、PM2.5与相对湿度、PM2.5与气压、O3与温度这时候手写一堆重复代码很笨也容易抄错变量名。我用循环来批量处理把结果和图片都保存到独立目录。vars {pm25, wind_speed, humidity, pressure, o3, temperature}; targetVar pm25; mkdir(output/coherence_png); for k 1:length(vars) if strcmp(vars{k}, targetVar) continue; end v1 data.(targetVar); v2 data.(vars{k}); % 跳过缺失过多的情况 valid ~(isnan(v1) | isnan(v2)); if sum(valid) length(v1) * 0.8 warning(变量 %s 的有效数据不足80%%跳过。, vars{k}); continue; end [wcoh, ~, freq] wcoherence(v1(valid), v2(valid), dt); figure(Visible, off); imagesc(t(valid), 1./freq, wcoh); axis xy; set(gca, YScale, log); colorbar; title(sprintf(%s 与 %s 的小波相干, targetVar, vars{k})); saveas(gcf, fullfile(output/coherence_png, sprintf(%s_%s.png, targetVar, vars{k}))); close(gcf); end批量处理时我最担心的是个别变量里有异常数据导致某一次循环计算特别慢甚至报错所以我习惯在循环里加try-catch块捕获错误后打印变量名继续下一组。宁可一次跑完看日志也不要跑一半中断然后整个人对着屏幕发懵。另外输出文件名要规范最好带上时间范围或者周期范围避免之后想找某个图却得一个个点开确认。3.4 绘图与导出让图表达到投稿标准Matlab默认的colormap是parula画出来颜色层次够用但不算出彩。我自己的经验是谈空气质量的小波相干图用蓝色到红色的渐变色会更直观中间用白或黄过渡。Matlab内置的turbo或者hot风格也都可以但注意不要让大面积的红色在没有显著性的区域出现那样会误导读者。我一般先把未通过显著性检验的区域透明度调低或者叠加斜线阴影这个过程在wtcplot里已经内置了在纯内置函数绘图时需要自己处理一下。导图也是个细节。如果是日常汇报直接保存png就够了如果打算放进论文最好输出矢量图。Matlab里可以这样写exportgraphics(gcf, coherence_pm25_ws.pdf, ContentType, vector);这样保存的是PDF矢量图放到LaTeX里或者Adobe Illustrator里再编辑都不会糊。早年我都是print -dpng后来换了exportgraphics才彻底解决图片模糊问题。还要注意图里字号期刊一般要求文字可读尺寸至少在7到8磅以上坐标轴标签用稍微大一点的字体也能提高审稿人印象分。4. 常见问题与排查技巧实录4.1 显著性检验图上的黑色等值线到底怎么来的小波相干图中那些黑色粗线不是凭空画出来的。无论是Matlab内置函数还是Grinsted工具包显著性检验的基本思路都是蒙特卡洛模拟假设两个序列之间不存在真实相关但各自保留与原始序列相似的随机结构然后反复生成多对替代数据计算每对替代数据的小波相干值统计出该时频位置上的95%分位数。如果原始数据的相干值超过这个阈值就认为在该时频点上两者的关系显著。替代数据的生成通常假设序列服从AR(1)过程也就是一阶自回归过程系数由原始序列估计得到。这也是为什么Grinsted工具包里要专门调用ar1函数去估算Lag1。内置wcoherence函数把这一步封装起来了你看不到细节但原理是类似的。这里有一个常见误区显著性检验的蒙特卡洛次数太少结果不稳定。比如我用MonteCarloCount100跑出来的显著区域和MonteCarloCount500跑出来的图有明显差异后来把次数提升到300以上才稳定下来。所以如果你最后得到的图看起来“这里红一块那里红一块”很零碎先别急着解读试试增大蒙特卡洛次数、看看显著性区域是否变得平滑。零碎的显著区域大概率是检验次数不够或者数据本身存在短促的突发事件带来的。4.2 边界效应与COI别把边缘的假结果当宝小波变换需要对有限长度的序列做卷积在序列开头和结尾小波基函数会伸出数据范围之外缺少真实数据的部分只能靠补零或其他方式填充结果就是边界附近的谱值被严重低估或扭曲。为了标记这个不可靠区域绘图时会在左右两侧画出锥形影响区COI通常用半透明阴影或细线围起来。我看到太多人忽略了COI分析一年数据结果图上最显眼的一大片红色出现在年初和年末的边界处还得出了“1月和12月关系特别强”的结论。实际上那很可能只是边界效应在起作用。正确做法是凡是在COI范围内的结果一律不作为核心证据最多只能作为“趋势暗示”下结论之前必须回到原始时间序列去核实。如果低频段整体都落在COI里比如小时数据想要分析年尺度周期那数据长度根本不够延长观测时间才是合理的解决方式换任何小波工具都救不了。我在数据处理阶段就会估算最长可用周期一般取数据长度的1/3到1/2作为可靠上限避免后面做无用功。4.3 代码报错与结果异常的快速排查用Matlab做小波相干时我整理过一张速查表出现对应问题直接按顺序检查现象可能原因排查/解决方式Undefined function or variable工具箱未安装或路径未添加检查ver命令是否显示Wavelet Toolbox用addpath(genpath(...))添加工具箱路径Matrix dimensions must agree两个输入序列长度不一致打印length(x)和length(y)对齐时间轴后重试NaN/Inf出现在结果矩阵原始数据含NaN或Inf清洗数据剔除缺失计算前用isfinite过滤图片全是深蓝色两个序列几乎没有共性先分别画单变量CWT确认特征检查数据是否已经标准化或去趋势显著区域太碎蒙特卡洛次数不足或数据中有强奇异值提高MonteCarloCount到300以上检查异常值是否被剔除干净纵轴频率和周期乱套使用freq还是1./freq混淆明确到底要展示频率还是周期画图前统一换算低频段全被COI覆盖数据长度不足以支撑该周期截断只看COI以外的区域或延长数据记录时间我遇到过最折磨人的问题是一次跑固定代码换了台电脑之后wcoherence报错查了半天发现是那台机器上Matlab版本太老wcoherence函数还不存在。后来处理方式很直接第一优先更新工具箱第二才去考虑代码兼容性。如果项目需要长期延续建议在README里写清楚依赖的Matlab版本号省得换环境时抓瞎。4.4 我自己踩过的坑和几个实操心得第一个坑是只看相干图不做单序列CWT检查。两个序列如果本身都没有周期性小波相干偶尔出现的红色区域经常是假象。我现在的习惯是先分别画出PM2.5和风速的连续小波功率谱确认各自的显著周期带再做相干分析。这样一来后面对显著区域的解释会可靠很多。第二个坑是相位箭头。刚开始读图时我把箭头方向完全当作“因果方向”后来认真看文档才明白箭头揭示的是相位关系不等同于因果关系。比如PM2.5与风速在某段时间内呈现反相只说明风速大时PM2.5低、风速小时PM2.5高并不能直接断言风速导致PM2.5下降。因果性需要结合物理过程、滞后分析和外部干预证据去论证。第三个心得是结果解释一定要回到原始序列做局部验证。图上发现6月前后有一个显著的4天周期相干区域我就把那段时间的原始序列截出来手工计算移动窗口相关确认关系确实存在且方向一致。这种验证动作耗时不多但能让你在汇报和写作时更有底气。如果原始序列上根本看不出对应关系那很可能整个分析过程某个环节出了问题这时候回头看数据清洗反而是最有效的手段。给自己的一个建议小波相干分析这个方法上手门槛并不高Matlab内置函数已经把大部头工作都封装好了但它真正值钱的地方在于你怎么解读那张时频图。空气质量和气象要素之间的关系往往具有很强的时变性和尺度依赖性小波相干能帮你把这些动态关系可视化出来不过它不能替你判断因果。我个人的实操体会是从一对变量开始跑通流程理解每一个输出矩阵的含义再批量扩展到多个要素交叉分析这样的节奏最稳。遇到奇怪的图不要急着改参数先回到原始数据检查很多时候问题不是出在算法上而是数据本身没洗干净。希望这份经验对你手头的项目有帮助。