python 小波分析和小波相干
文章目录
- python 小波分析和小波相干
- 写在前面
- 小波分析
- 小波相干
- 小结
写在前面
这个东西单纯为了记录一下用Python做小波分析和小波相干的过程，以及怎么理解得的的结果。MATLAB中有很成熟的小波工具包，python中也有好几个可以做小波分析的包。我自己使用了Pyleoclim，这个包主要是用来做古气候数据处理和分析的，主要的功能如下图。这个文章只涉及小波分析部分，代码在链接的教程中都要。Pyleoclim包目前支持Python 3.9以上版本。再次强调，代码都来自于链接教程！！！。
小波分析
例子中用到了两个古气候数据，一个是同位素氘，另一个是。
首先加载需要用到的包
import matplotlib.pyplot as plt
import pandas as pd
import pyleoclim as pyleo
import numpy as np然后读取数据，并将数据存入pyleoclim包中的 Series类型，主要给定序列的名称和时间。
dDdf = pd.read_csv('../data/edc3deuttemp2007.csv')
dDdf.head()
dDts = pyleo.Series(time=dDdf['Age']/1000,value=dDdf['Deuterium'],time_name='Age',
time_unit='kyr BP',value_name = r'$\deltaD$',
value_unit=u'\u2030',label=r'EPICA Dome C $\delta D$')
dDts.plot()
需要注意的是，这个数据的采样点在时间上并不平均，即相邻两个数据点的时间间隔并不一致。
这不满足一些分析的基本要求，或者会增加一些分析的难度和不确定性。
fig, ax = plt.subplots()
plt.hist(np.diff(dDts.time))
plt.xlabel(r'$\Delta t (kyr)$')
plt.ylabel('# occurrences')
plt.title('Distribution of age increments',weight='bold')
因此，需要在时间上以500年的时间分辨率进行插值，并对序列进行标准化。
dD05 = dDts.interp(step=0.5).standardize()小波分析 直接通过调用 wavelet 方法来实现
scal = dD05.wavelet()
scal.plot()
pyleoclim给了两种计算显著性的方式，分别是参数化逼近方法 ar1asym，另一种是模拟方法 ar1sim
ar1asym方法
scal_sig = scal.signif_test(method='ar1asym')
scal_sig.plot()
ar1sim方法
scal_sig2 = scal.signif_test(method='ar1sim', number = 1000)
scal_sig2.plot()
这个结果显示该时间序列在300-600 kyr BP存在显著的100kyr周期。
小波相干
同位素氘可以作为指示温度变化的指标， 接下来读入序列，并使用小波相干分析两者关系。
url = 'ftp://ftp.ncdc.noaa.gov/pub/data/paleo/icecore/antarctica/antarctica2015co2composite.txt'
co2df = pd.read_csv(url,skiprows=137,sep='\t')
co2df.head()
co2ts = pyleo.Series(time=co2df['age_gas_calBP']/1000,value= co2df['co2_ppm'],
time_name='Age',time_unit='kyr BP',value_name = r'$CO_2$',
value_unit='ppm',label='EPICA Dome C CO2')
co2ts.plot(color='C1')
从变化图中可以明显看到和同位素氘序列类似的冰期间冰期循环，以及工业革命后的快速爬升。
将两者画在一张图上
ms = pyleo.MultipleSeries([dDts,co2ts])
ms.plot()
对两者都进行标准化
ms.standardize().plot()
调用wavelet_coherence方法进行小波相干分析
coh = co2ts.wavelet_coherence(dDts,method='wwz')
fig, ax = coh.plot()
小波相干-wavelet transform coherency (WTC) 的结果可以理解为两个变量在时频域的相关性的平方，0表示不相关，1表示相关； 箭头表示两个变量在某个时频上的相位关系，根据三角函数转换，右 = 0度, 上 = 90度, 左 = 180度, 下 = -90 or + 270度。
图中基本向右的箭头表示 和 氘 的变化在相应的时间和频域上的相位角接近零。
可以计算一定频率区间内的平均相位角。
coh.phase_stats(scales=[80,100])结果为：
Results(mean_angle=0.2855352895679047, kappa=25.30163136494105, sigma=0.20083246106901098, kappa_hi=486.4211381097546, sigma_lo=0.045364651843896514)这个结果说明，在80-100kyr的轨道周期上，这两个序列的变化基本是同相位的，两者之间不存在显著的提前或者滞后。
同样，也可以计算小波想干的显著性。
coh_sig = coh.signif_test(number=100)
fig, ax = coh_sig.plot()
借用matab小波工具箱文档 来解释箭头方向的意义。
下面图片是matlab文档中 nino3指数（热带太平洋某区域的海温，ENSO）与印度降水的小波相干图。
图中表示nino3指数和印度降水在2-7年周期上（ENSO的主要周期）具有较强的相关，箭头的方向说明两者之前存在3/8到1/2（箭头向左） 周期的滞后（即1-3.5年）。这张图就说明热带太平洋海温（ENSO）与印度的降水存在相关，但是这影响存在1-3.5年的滞后。
小结
本文主要是记录一下学习小波分析和小波相干分析的笔记。相关原始内容在文中的链接中都可以找到，当然我个人的理解可能存在一定问题，希望大家有自己的认识。再次强调，虽然本人已经使用pyleoclim包做了自己的分析，但是为了方便，本文中使用的代码和图片均来自于官方文档。
