交叉小波变换和小波相干在地球物理时间序列中的应用

摘要

许多科学家利用小波方法来分析时间序列，通常使用流行的免费软件。 然而，目前还没有类似的易于使用的小波包来一起分析两个时间序列。 我们讨论交叉小波变换和小波相干性来检查两个时间序列之间的时频空间关系。 我们演示了如何使用相位角统计来获得对因果关系的信心并测试时间序列之间物理关系的机械模型。 作为此类分析被证明有用的典型数据示例，我们将这些方法应用于北极涛动指数和波罗的海最大海冰范围记录。 蒙特卡罗方法用于评估红噪声背景下的统计显着性。 已经开发了一个软件包，允许用户执行交叉小波变换和小波相干（http://www.pol.ac.uk/home/research/waveletcoherence/）。

1 Introduction

        地球物理时间序列通常是由我们知之甚少的复杂系统生成的。 因此，此类系统中的可预测行为（例如趋势和周期性）引起了人们的极大兴趣。 大多数检查频域周期性的传统数学方法（例如傅里叶分析）都隐含地假设基础过程在时间上是平稳的。 然而，小波变换将时间序列扩展到时频空间，因此可以找到局部间歇周期性。 小波变换有两类： 连续小波变换 (CWT) 及其离散小波变换 (DWT)。 DWT 是数据的紧凑表示，对于降噪和数据压缩特别有用，而 CWT 更适合特征提取。

        由于我们对提取时间序列中的低信噪比信号感兴趣，因此我们在本文中仅讨论 CWT。 虽然 CWT 是分析时间序列中局部间歇性振荡的常用工具，但通常需要同时检查两个可能以某种方式关联的时间序列。 特别是，检查时频空间中具有较大公共功率的区域是否具有一致的相位关系，从而暗示时间序列之间的因果关系。 许多地球物理时间序列不是正态分布的，我们建议将 CWT 应用于此类时间序列的方法。 我们从两个 CWT 构建交叉小波变换 (XWT)，它将揭示它们在时频空间中的共同功率和相对相位。 我们将进一步定义两个 CWT 之间的小波相干性 (WTC) 度量，即使公共功率较低，它也可以找到显着的相干性，并展示如何计算针对红噪声背景的置信水平。
在讨论 XWT 和 WTC 之前，我们将介绍基本的 CWT 理论。 新的发展，例如量化相位关系和计算 WTC 显着性水平，将得到更全面的处理。 在时间序列上使用这些方法时，重要的是要有坚实的机械基础来建立所发现的任何关系，并且我们警告不要在“散布枪”方法中使用这些方法（特别是如果时间序列概率密度函数是 已修改）。 为了说明如何使用各种方法，我们将它们应用于气象学和冰川学的两个数据集。 最后，我们将提供 MatLab 软件包的链接。

2 Data

我们期望通过考虑气候系统而联系起来的两个物理效应的一个例子是北极大气的平均冬季状态和冰条件反映的冬季严重程度。

北极涛动（AO）是北半球气候变化的一个重要方面。 AO 被定义为北半球 20°N 极地海平面压力异常的主要经验正交函数（EOF）（Thompson 和 Wallace，1998），其特征是北极和中纬度之间大气质量的交换 。 波罗的海是北大西洋地区和欧亚大陆地区之间的过渡带，导致冰情年际变化较大。 波罗的海每年冬季部分被冰覆盖，年最大冰面积占海域面积的10-100%，冰季长度为4-7个月，年最大冰厚度为50-120 cm（Jevrejeva，2001；Seiná 等人，2001）。 最近发表的结果表明，由 AO 或有点类似的北大西洋涛动遥相关描述的北极和北大西洋的大规模大气环流模式显着控制着波罗的海的冰况（Loewe 和 Koslowski，1998；Omstedt 和 Chen） ，2001 年；耶夫热耶娃和摩尔，2001 年；耶夫热耶娃，2002 年）。

在本文中，我们研究了冬季AO与波罗的海冰范围之间的联系，特别是根据预期的因果关系探讨了该系列之间的相位关系。 冰况由 1720 年至 2000 年期间波罗的海最大年冰范围（BMI）的时间序列表示（Seiná 等人，2001 年）。 我们使用 Thompson 和 Wallace (1998) 的冬季 AO 指数（1851 年 12 月至 1997 年 2 月）。
许多统计测试假设概率密度函数 (pdf) 接近正态。 我们对地球物理时间序列 CWT 的经验表明，远离正态分布的序列会产生相当不可靠且不太显着的结果。 因此，有时转换时间序列的 pdf 可能是个好主意。 但是，我们警告不要贸然更改 pdf。 BMI 指数呈双峰分布，最大概率约为 70 000 km2 和 420 000 km2。 这两种模式之间的简单周期性振荡几乎具有方波的形状。 方波的功率泄漏到基本周期之外的频带中。 因此，我们将 BMI 指数转换为百分位数记录（就其累积分布函数而言），从而迫使 pdf 成为矩形。 这具有缩小间歇振荡的带宽的效果。 我们对两个时间序列进行标准化（零均值，单位标准差），并将标准化版本简单地称为 AO 和 BMI。 AO和BMI时间序列的时间序列如图1和2所示。 1 和 2。

3 Methods

3.1 The Continuous Wavelet Transform (CWT)

小波是均值为零且在频率和时间上都局部化的函数。 我们可以通过小波在时间 (Δt) 和频率（Δω 或带宽）上的局部性来表征它。 海森堡不确定性原理的经典版本告诉我们，时间局部化和频率局部化之间总是存在权衡。 如果没有正确定义 Δt 和 Δω，我们会注意到不确定性乘积 Δt·Δω 的大小是有限的。 一种特殊的小波 Morlet 定义为

其中 ω0 是无量纲频率，η 是无量纲时间。 当使用小波进行特征提取时，Morlet 小波（ω0=6）是一个不错的选择，因为它在时间和频率局部化之间提供了良好的平衡。 因此，我们将进一步的处理限制在这个小波上，尽管我们提出的方法通常是适用的（参见例如Foufoula-Georgiou，1995）。

CWT 背后的想法是将小波作为带通滤波器应用于时间序列。 通过改变小波的尺度 (s) 来及时拉伸小波，使得 η = s·t，并将其归一化为具有单位能量。

对于 Morlet 小波（ω0=6），傅立叶周期 (λwt ) 几乎等于尺度 (λwt =1.03 s)。 具有统一时间步 δt 的时间序列 (xn, n=1,...,N) 的 CWT 被定义为 xn 与缩放和归一化小波的卷积。 我们写道，在实践中，在傅里叶空间中实现卷积速度更快（详细信息请参阅 Torrence 和 Compo，1998）。 我们将小波功率定义为 |WnX(s)|2。 WnX(s) 的复数参数可以解释为局部相位。

CWT 具有边缘伪影，因为小波没有完全及时定位。 因此，引入影响锥（COI）是有用的，其中边缘效应不能被忽略。 这里我们将 COI 视为由边缘不连续性引起的小波功率下降到边缘值的 e^−2 的区域。
小波功率的统计显着性可以相对于原假设进行评估，即信号是由具有给定背景功率谱 (Pk) 的平稳过程生成的。 许多地球物理时间序列具有独特的红噪声特征，可以通过一阶自回归 (AR1) 过程很好地建模。 具有滞后 1 自相关 α 的 AR1 过程的傅立叶功率谱（根据观察到的时间序列估计，例如 Allen 和 Smith，1996）由下式给出

其中 k 是傅里叶频率指数。
小波变换可以被认为是应用于时间序列的一系列连续的带通滤波器，其中小波尺度与滤波器的特征周期（λwt）线性相关。 因此，对于具有功率谱 Pk 的平稳过程，通过调用傅立叶卷积定理，给定小波尺度的方差就是 Pk 相应频带中的方差。 如果 Pk 足够平滑，那么我们可以使用转换 k−1=λwt 简单地用 Pk 来近似给定尺度下的方差。 Torrence 和 Compo (1998) 使用蒙特卡罗方法表明这种近似对于 AR1 谱来说非常好。 然后他们表明，给定功率谱 (Pk) 的过程的小波功率大于 p 的概率为

其中，对于实小波，v 等于 1；对于复小波，v 等于 2。
AO 和 BMI 的 CWT 如图 1 和 2 所示。 分别为1和2。 两个时间序列的小波功效有明显的共同特征，例如 1940 年左右的 5 年波段有显着的峰值。尽管这两个系列在 1860 年至 1900 年期间的 2 年至 7 年波段也有高功效。 对于 AO，功效不高于 5% 显着性水平。 然而，这一时期所描绘的图案之间的相似度很低，因此很难判断这是否仅仅是巧合。 交叉小波变换在这方面有所帮助。

3.2 The cross wavelet transform

两个时间序列 xn 和 yn 的交叉小波变换 (XWT) 定义为 WXY =WXWY *，其中 * 表示复共轭。 我们进一步定义交叉小波幂为|WXY |。 复数参数 arg(W xy ) 可以解释为时频空间中 xn 和 yn 之间的局部相对相位。 Torrence 和 Compo (1998) 给出了具有背景功率谱 P X 的两个时间序列的交叉小波功率的理论分布

其中 Z (p) 是与由两个 Χ2 分布乘积的平方根定义的 pdf 的概率 p 相关的置信水平。 例如，图 3 中标记的 5% 显着性水平是使用 Z2(95%)=3.999 计算的。
AO 和 BMI 的 XWT 如图 4 所示。在这里，我们注意到，我们通过肉眼从各个小波变换中发现的共同特征在 5% 的水平上非常显着。 我们还注意到，1940 年至 1980 年的 10 至 16 年区间也存在显着的共同权力。

为了在时间序列中记录的现象之间存在简单的因果关系，我们期望振荡是锁相的。 因此，令人欣慰的是，XWT 显示 AO 和 BMI 在所有具有显着共同权力的部门中处于反相状态。 由于 AO 和 BMI 在所有尺度上都是反相的，因此我们得出结论，BMI 在很大程度上只是反映了 AO。 在具有显着功率的区域之外，相位关系也主要是反相位的。 因此，我们推测 AO 和 BMI 之间的联系比交叉小波幂所暗示的联系更强。

3.3 Cross wavelet phase angle

由于我们对两个时间序列分量之间的相位差感兴趣，因此我们需要估计相位差的平均值和置信区间。 我们使用 COI 之外统计显着性高于 5% 的区域的相位循环平均值来量化相位关系。 这是计算平均相位的有用且通用的方法。 一组角度 (ai , i=1...n) 的圆平均值定义为（例如 Zar，1999）

由于相位角不是独立的，因此很难可靠地计算平均角度的置信区间。 只需增加刻度分辨率即可将计算中使用的角度数量设置为任意高。 然而，了解平均值周围角度的分散情况很有趣。 为此，我们将循环标准差定义为

其中 R=p(X2+Y 2)。 圆形标准差与线性标准差类似，它的变化范围是从零到无穷大。 当角度紧密分布在平均角度周围时，它会给出与线性标准差类似的结果。 在某些情况下，计算每个尺度的平均相位角可能是有原因的，然后相位角可以量化为几年。
5% 有效区域内和 COI 之外的 XWT 相位角的平均相位为 –176±12（其中±表示圆标准差）。 这基本上证实了AO和BMI反相的结论。 请注意，由于四月份冰面积最大，时间序列已经有 3 个月的滞后。 3 个月的滞后与对流层平流层强迫的机制一致（Baldwin 等，2001；Jevrejeva 等，2003）。 相位角在所有尺度上都是恒定的观察结果表明，由于信号从 AO 传播到冰范围的物理机制，存在恒定的时滞。 完全反相的偏差表明 AO 稍微领先 BMI，但是，循环标准偏差太大，无法得出任何确定的结论。

3.4 Wavelet coherence

交叉小波功率揭示了具有高公共功率的区域。 另一个有用的度量是交叉小波变换在时频空间中的相干程度。 遵循 Torrence 和 Webster (1998)，我们将两个时间序列的小波相干性定义为

其中 S 是平滑算子。 请注意，这个定义与传统相关系数的定义非常相似，并且将小波相干性视为时频空间中的局部相关系数是有用的。 我们将平滑算子 S 写为

其中Sscale表示沿小波尺度轴的平滑，Stime表示时间上的平滑。 设计平滑算子使其具有与所使用的小波相似的足迹是很自然的。 对于 Morlet 小波，Torrence 和 Webster (1998) 给出了合适的平滑算子

其中 c1 和 c2 是归一化常数，5 是矩形函数。 因子 0.6 是根据经验确定的 Morlet 小波的尺度去相关长度（Tor-rence 和 Compo，1998）。 实际上，两个卷积都是

离散地进行，因此标准化系数是通过数字确定的。
小波相干性的统计显着性水平是使用蒙特卡罗方法估计的。 我们生成了一个大型的代理数据集对集合，其 AR1 系数与输入数据集相同。 对于每一对，我们计算小波相干性。 然后，我们仅使用 COI 之外的值来估计每个量表的显着性水平。 经验检验表明，AR1 系数对显着性水平影响不大。 然而，平滑算子的细节有很大的影响。 例如，计算尺度平滑时选择的分辨率对显着性水平有重大影响（见图 4）。 显着性水平的蒙特卡罗估计需要 1000 个替代数据集对。 每个八度音阶的尺度数量应该足够高，以捕获尺度平滑算子的矩形形状，同时最大限度地减少计算时间。 根据经验，我们发现每个八度 10 个音阶就足够了。
AO 和 BMI 的 WTC 平方如图 5 所示。与 XWT 相比，更大的部分显得很重要，所有这些区域都显示 AO 和 BMI 之间的反相位关系。 高于 5% 显着性水平的时频图面积并不是因果关系的可靠指示。 即使对尺度进行适当权衡以求平均值，两个系列也有可能在一个特定尺度上完全相关，而显着相关的面积远小于 5%。 然而，图 5 的重要区域是如此广泛，以至于这不太可能只是偶然的。 AO 的振荡在 BMI 中表现出 2-20 年不等的波长，这表明 BMI 被动地反映了 AO。 低相干性区域与 AO 中的低小波功率一致，因此是预期的。 可能是因为 AO 代表的 EOF 并没有真正捕捉到小冰河时期行动中心的实际位置。 与 XWT 一样，我们可以计算重要区域的平均相位角。 小波相干性显着的区域和 COI 之外的平均相位角为 174±15°。

4 Summary

CWT 将时间序列扩展到时频空间，可以以高度直观的方式看到振荡。 当使用小波进行特征提取时，Morlet 小波（！0=6）是一个不错的选择，因为它在时间和频率上都有合理的局部化。 从两个时间序列的 CWT 可以构建 XWT。 XWT 暴露了具有高公共功率的区域，并进一步揭示了有关相位关系的信息。 如果这两个系列在物理上相关，我们预计会出现一致或缓慢变化的相位滞后，可以根据物理过程的机械模型进行测试。 WTC 可以被认为是两个 CWT 之间的局部相关性。 通过这种方式，局部锁相行为就被发现了。 WTC 更理想的功能是以时频空间局部性稍差为代价的。 WTC 的显着性水平必须使用蒙特卡罗方法来确定。

4.1 Practical tips

交叉小波分析和小波相干性是测试两个时间序列之间所提出的联系的强大方法。

–    检查时间序列的直方图，确保它们与正态分布相差不远。 如果时间序列的 pdf 远离高斯分布，请考虑转换时间序列。 选择变换时，最好选择解析变换，例如如果数据呈对数正态分布，则取对数。 在其他情况下，我们用于 BMI 的简单“百分位”变换可能有用。 使用该特定转换的优点是它没有任何异常值。
–   考虑给定提议的链接机制对分析结果的期望是什么。 我们警告不要盲目地将这些方法应用于随机选择的数据集。 与其他统计测试一样，某些数据集会偶然显示出具有高度统计意义的链接。
–    选择小波后，将计算两个时间序列的 CWT。 我们建议音阶分辨率为每八度 10 个音阶，并使用 Morlet 小波，除非有充分理由不这样做。 对于地球物理时间序列，AR1 红噪声假设通常是十个合适的，并且等式： (3)和(4)可以用来计算小波功率的显着性水平。 请记住要特别小心，不要误解 COI 内的结果。

–    根据两个 CWT 计算 XWT。 XWT 暴露了具有高公共功率的区域，并进一步揭示了有关相位关系的信息。 如果这两个系列在物理上相关，我们预计会出现一致或缓慢变化的相位滞后，可以根据物理过程的机械模型进行测试。 相位角的圆形平均值可用于量化相位关系。
–    此外，可以从两个 CWT 计算出 WTC，它可以被视为时频空间中时间序列之间的局部相关性。 XWT 揭示了高公共功率，而 WTC 则发现了局部锁相行为。 WTC 更理想的功能是以时频空间局部性稍差为代价的。 WTC 的显着性水平必须使用蒙特卡罗方法来确定。

作者用于执行 XWT 和 WTC 的 MatLab 软件包可以在 http://www.pol.ac.uk/home/research/waveletcoherence/ 找到。
致谢。 图勒研究所提供了财政援助。 我们的一些软件包含最初由 C. Torrence 和 G. Compo 编写的代码，可在以下网址获取：http://paos.colorado.edu/research/wavelets/ 以及由阿拉斯加大学的 E. Breiten-berger 编写的代码，并进行了改编 来自免费软件 SSA-MTM 工具包：http://www.atmos.ucla.edu/tcd/ssa。
编辑：M.蒂尔
评审人： 两名裁判员

References
Allen, M. R. and Smith, L. A.: Monte Carlo SSA: Detecting irreg-
ular oscillations in the presence of coloured noise, J. Clim., 9, 3373–3404, 1996.
Baldwin, M. P., Gray, L. J., Dunkerton, T. J., Hamilton, K., Haynes,P. H., Randel, W. J., Holton, J. R., Alexander, M. J., Hirota, I., Horinouchi, T., Jones, D. B. A., Kinnerslay, J. S., Marquardt, C., Sato, K., and Tarahashi, M.: The Quasi-Biennial Oscillation, Rev. Geophys., 39, 170–229, 2001.

Foufoula-Georgiou, E. and Kumar, P.: Wavelets in Geophysics,
Academic Press, 373, 1995.
Jevrejeva, S. and Moore, J. C.: Singular Spectrum Analysis of Baltic Sea ice conditions and large-scale atmospheric patterns since 1708, Geophys. Res. Lett., 28, 4503–4507, 2001.

Jevrejeva, S.: Association between the ice conditions in the Baltic Sea and the North Atlantic Oscillation, Nordic Hydrol., 33, 319–330, 2002.

Jevrejeva, S., Moore, J. C., and Grinsted, A.: Influence of the Arc-tic Oscillation and El Nin˜o-Southern Oscillation (ENSO) on ice conditions in the Baltic Sea: The wavelet approach, J. Geophys. Res., 108(D21), 4677, doi:10.1029/2003JD003417, 2003.

Loewe, P. and Koslowski, G. : The Western Baltic Sea ice season in terms of a mass-related severity index 1879–1992, II Spectral characteristics and associations with NAO, QBO and solar cycle, Tellus, 50A, 219–241, 1998.

Omstedt, A. and Chen, D.: Influence of atmospheric circulation on the maximum ice extent in the Baltic Sea, J. Geophys. Res., Vol. 106, 4493–4500, 2001.

Sein¨a, A., Gr¨onvall, H., Kalliosaari, S., and Vainio, J.: Ice seasons 1996–2000 in Finnish sea areas / J¨a¨atalvet 1996–2000 Suomen merialueilla, Meri, Report Series of the Finnish Institute of Ma-rine Research, 43, 2001.

Thompson, D. W. J. and Wallace, J. M.: The Arctic Oscillation sig-nature in the winter geopotential height and temperature fields, Geophys. Res. Lett., 25, 1297–1300, 1998.

Torrence, C. and Compo, G. P.: A practical guide to wavelet analy-
sis, Bull. Am. Meteorol. Soc., 79, 61–78, 1998.
Torrence, C. and Webster, P.: Interdecadal Changes in the ESNO-
Monsoon System, J. Clim., 12, 2679–2690, 1999.
Zar, J. H.: Biostatistical Analysis, Prentice hall, New Jersey, 1999.

Fig. 1. The standardized time series of winter (DJF) AO (bottom) and its continuous wavelet power spectrum (top). The thick black contour designates the 5% significance level against red noise and the cone of influence (COI) where edge effects might distort the picture is shown as a lighter shade. The standardized AO series has an AR1 coefficient of 0.02.

Fig. 2. The standardized BMI percentile time series (bottom) and its continuous wavelet power (top). The thick contour designates the 5% significance level against red noise and the cone of influence (COI) where edge effects might distort the picture is shown as a lighter shade. The standardized BMI percentile series has an AR1 coefficient of 0.08.

Fig. 3. Cross wavelet transform of the standardized AO and BMI time series. The 5% significance level against red noise is shown as a thick contour. The relative phase relationship is shown as arrows (with in-phase pointing right, anti-phase pointing left, and BMI leading AO by 90 pointing straight down).

Fig. 4. Wavelet coherence 5% significance level determined using Monte Carlo generated noise (with 10 000 surrogate data set pairs). The legend shows the 2 AR1 coefficients of the surrogate data sets and the numbers of scales per octave (s/o) used in calculating the scale smoothing. The color of the noise has little impact on the significance level, whereas the specifics of the smoothing have a large impact. The large values at either end of the spectrum are due to the scale smoothing operator reaching the scale boundaries.

Fig. 5. Squared wavelet coherence between the standardized AO and BMI time series. The 5% significance level against red noise is shown as a thick contour. All significant sections show anti-phase behavior.