用Python复现气象顶刊图表：手把手教你做小波分析（附Torrence-Compo代码）

Python气象科研实战：从Torrence-Compo代码到顶刊级小波分析图表

第一次看到《地理科学》那篇长江中下游降水研究的功率谱图时，我被那些优美的等高线圈和醒目的显著性区域震撼到了——这不正是我论文需要的关键证据吗？但当我打开MATLAB准备复现时，面对一堆看不懂的cwt函数参数和神秘的"影响锥"概念，整整两周毫无进展。直到我发现科罗拉多大学Torrence和Compo教授那篇被引7000+次的经典论文，以及Predybaylo博士移植的Python版本，才真正打开了小波分析的大门。

1. 环境配置与数据准备

工欲善其事，必先利其器。我们先搭建一个专为气象数据分析优化的Python环境：

conda create -n wavelet python=3.8 conda install -c conda-forge numpy scipy matplotlib jupyter

Torrence-Compo官方代码库提供了两个关键文件：

• waveletFunctions.py：核心算法模块
• waveletAnalysis.py：示例分析脚本

建议直接克隆原始仓库：

import urllib.request urllib.request.urlretrieve( "https://raw.githubusercontent.com/chris-torrence/wavelets/master/python/waveletFunctions.py", "waveletFunctions.py" )

准备你的时间序列数据时要注意：

• 数据应该是一维数组，如某站点的年降水量序列
• 缺失值需提前处理（线性插值或剔除）
• 时间间隔dt要准确设置（月数据为1/12，年数据为1）

import numpy as np # 示例：南京站1901-2020年夏季降水量(mm) nanjing_rain = np.array([...]) dt = 1 # 年数据

2. 小波变换核心参数解析

运行小波变换前，这些参数决定了分析的质量：

# 关键参数设置 pad = 1 # 推荐填充零值到2的幂次长度 dj = 0.25 # 尺度间隔（每octave分4个子尺度） s0 = 2*dt # 最小尺度（2倍时间间隔） J1 = 7/dj # 尺度数量 mother = 'MORLET' # 小波母函数

Morlet小波的波形控制参数k0尤为重要：

• k0=6（默认）：平衡时频分辨率
• k0>6：时间分辨率↑，频率分辨率↓
• k0<6：频率分辨率↑，时间分辨率↓

气象数据通常使用k0=6，生物医学信号可能用k0=2

通过调整这些参数，我曾发现某期刊论文中声称的"20年周期"其实是参数设置不当导致的伪信号——把dj从0.5改为0.25后，那个显著峰就消失了。

3. 结果可视化与学术规范

得到小波功率谱后，如何制作符合顶刊要求的图表？看这段改进版的绘图代码：

def plot_wavelet(time, data, wave, period, coi, sig95, global_ws, global_sig): plt.figure(figsize=(10, 8)) # 时间序列子图 plt.subplot(3, 1, 1) plt.plot(time, data, 'k', linewidth=1.5) plt.ylabel('Precipitation (mm)') plt.title('a) Nanjing Summer Rainfall Anomalies') # 功率谱子图 plt.subplot(3, 1, 2) levels = np.linspace(0, np.max(power), 20) plt.contourf(time, period, power, levels, cmap='jet') plt.colorbar(label='Power (°C²)') plt.contour(time, period, sig95, [1], colors='k', linewidths=2) plt.plot(time, coi, 'k--') plt.yscale('log') plt.ylabel('Period (years)') plt.title('b) Wavelet Power Spectrum') # 全局谱子图 plt.subplot(3, 1, 3) plt.plot(global_ws, period, 'k', label='Power') plt.plot(global_sig, period, 'r--', label='95% significance') plt.yscale('log') plt.xlabel('Power (°C²)') plt.title('c) Global Wavelet Spectrum')

几个让审稿人眼前一亮的细节：

1. 使用jet色带增强冷暖对比
2. 显著性检验线加粗到2pt
3. 子图标题用a)、b)编号
4. 坐标标签带单位
5. 影响锥(COI)用虚线表示

4. 高级分析与论文写作技巧

小波分析结果如何转化为论文中的科学发现？这里有个真实案例：

发现1：交叉小波揭示遥相关

# 计算两个序列的交叉小波谱 cross_power = np.abs(wave1 * wave2.conj())

发现2：尺度平均方差突显年代际变化

# 计算2-8年尺度平均 scale_avg = (scale >= 2) & (scale < 8) avg_power = np.mean(power[scale_avg, :], axis=0)

在论文方法部分要注明： "小波分析采用Torrence和Compo(1998)的方法，显著性检验基于红噪声背景谱(滞后1自相关=0.72)"

讨论部分可以这样写： "3.6年的显著周期(图3b)与ENSO的典型周期一致，这与Wang等(2017)的发现相吻合。但值得注意的是，在1980年后该周期信号明显减弱，可能反映了太平洋年代际振荡(PDO)相位转变的影响。"

5. 常见陷阱与解决方案

问题1：边界效应失真

• 对策：始终关注COI曲线外的结果不可靠
• 代码检查：

coi_mask = np.tile(period, (len(time),1)).T > np.tile(coi, (len(period),1)) power[coi_mask] = np.nan

问题2：虚假显著周期

• 对策：尝试不同的lag1自相关值
• 诊断代码：

for lag1 in [0.6, 0.7, 0.8]: sig95 = wave_signif(variance, dt, scale, lag1=lag1) # 比较不同lag1下的显著区域变化

问题3：尺度分辨率不足

• 现象：周期轴出现阶梯状伪影
• 改进：减小dj到0.125，增加J1到56

那次我帮学弟检查论文，发现他漏掉了COI区域的数据筛选，差点把边界效应导致的伪信号当作重大发现。现在他的致谢里还留着我的名字——科学研究的严谨性就体现在这些细节里。