企业官网建设流程全解析

做空气质量数据分析那会儿，我最常用的是把两个变量的相关系数算一算，比如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)红噪声的蒙特卡洛显著性检验。

我做了个简单的对比：

如果只是自己分析数据，优先用内置函数；如果要写论文且希望图和审稿人预期一致，建议用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函数把这一步封装起来了，你看不到细节，但原理是类似的。

这里有一个常见误区：显著性检验的蒙特卡洛次数太少，结果不稳定。比如我用MonteCarloCount=100跑出来的显著区域和MonteCarloCount=500跑出来的图有明显差异，后来把次数提升到300以上才稳定下来。所以如果你最后得到的图看起来“这里红一块那里红一块”很零碎，先别急着解读，试试增大蒙特卡洛次数、看看显著性区域是否变得平滑。零碎的显著区域大概率是检验次数不够或者数据本身存在短促的突发事件带来的。

4.2 边界效应与COI：别把边缘的假结果当宝

小波变换需要对有限长度的序列做卷积，在序列开头和结尾，小波基函数会伸出数据范围之外，缺少真实数据的部分，只能靠补零或其他方式填充，结果就是边界附近的谱值被严重低估或扭曲。为了标记这个不可靠区域，绘图时会在左右两侧画出锥形影响区（COI），通常用半透明阴影或细线围起来。

我看到太多人忽略了COI：分析一年数据，结果图上最显眼的一大片红色出现在年初和年末的边界处，还得出了“1月和12月关系特别强”的结论。实际上那很可能只是边界效应在起作用。正确做法是，凡是在COI范围内的结果一律不作为核心证据，最多只能作为“趋势暗示”，下结论之前必须回到原始时间序列去核实。

如果低频段整体都落在COI里，比如小时数据想要分析年尺度周期，那数据长度根本不够，延长观测时间才是合理的解决方式，换任何小波工具都救不了。我在数据处理阶段就会估算最长可用周期，一般取数据长度的1/3到1/2作为可靠上限，避免后面做无用功。

4.3 代码报错与结果异常的快速排查

用Matlab做小波相干时，我整理过一张速查表，出现对应问题直接按顺序检查：

我遇到过最折磨人的问题是一次跑固定代码，换了台电脑之后wcoherence报错，查了半天发现是那台机器上Matlab版本太老，wcoherence函数还不存在。后来处理方式很直接：第一优先更新工具箱，第二才去考虑代码兼容性。如果项目需要长期延续，建议在README里写清楚依赖的Matlab版本号，省得换环境时抓瞎。

4.4 我自己踩过的坑和几个实操心得

第一个坑是只看相干图，不做单序列CWT检查。两个序列如果本身都没有周期性，小波相干偶尔出现的红色区域经常是假象。我现在的习惯是先分别画出PM2.5和风速的连续小波功率谱，确认各自的显著周期带，再做相干分析。这样一来，后面对显著区域的解释会可靠很多。

第二个坑是相位箭头。刚开始读图时，我把箭头方向完全当作“因果方向”，后来认真看文档才明白，箭头揭示的是相位关系，不等同于因果关系。比如PM2.5与风速在某段时间内呈现反相，只说明风速大时PM2.5低、风速小时PM2.5高，并不能直接断言风速导致PM2.5下降。因果性需要结合物理过程、滞后分析和外部干预证据去论证。

第三个心得是结果解释一定要回到原始序列做局部验证。图上发现6月前后有一个显著的4天周期相干区域，我就把那段时间的原始序列截出来，手工计算移动窗口相关，确认关系确实存在且方向一致。这种验证动作耗时不多，但能让你在汇报和写作时更有底气。如果原始序列上根本看不出对应关系，那很可能整个分析过程某个环节出了问题，这时候回头看数据清洗反而是最有效的手段。

给自己的一个建议

小波相干分析这个方法，上手门槛并不高，Matlab内置函数已经把大部头工作都封装好了；但它真正值钱的地方在于你怎么解读那张时频图。空气质量和气象要素之间的关系往往具有很强的时变性和尺度依赖性，小波相干能帮你把这些动态关系可视化出来，不过它不能替你判断因果。我个人的实操体会是，从一对变量开始跑通流程，理解每一个输出矩阵的含义，再批量扩展到多个要素交叉分析，这样的节奏最稳。遇到奇怪的图不要急着改参数，先回到原始数据检查，很多时候问题不是出在算法上，而是数据本身没洗干净。希望这份经验对你手头的项目有帮助。