Python与大气科学气象时序学习

Python气象数据处理与绘图(16):去趋势与滤波

2020-03-14  本文已影响0人  摸鱼咯

在最开始,我首先给出一个原始序列,序列并没有任何意义,蓝线为原始序列,黄线为其线性回归,可以看到序列有一个明显的线性上升的趋势,并且序列有很强的年际变化和年代际变化特征。


原始序列及其趋势

去趋势

研究长时间的气候变化会发现很多序列是带有趋势的,有部分观点认为这是认为因素造成全球变暖导致的,研究本身的大气变化需要去掉序列的趋势,简单来说就是序列的每一个点都减去序列趋势的回归斜率乘以一个时间步长,当然前提是这个序列的线性趋势是显著的。还有一个要注意的是去趋势序列(若干年长度的序列)并不等于该序列的年际变化,获得年际变化分量需要去用滤波或者其他方法单独提取,我目前看到的大部分文献中的去趋势都是为了剔除人为因素造成的全球变暖。
scipy直接提供了一个去趋势函数,十分的方便。

scipy.signal.detrend(data, axis=-1, type='linear', bp=0, overwrite_data=False)
#data:输入数据,可以为任意维度
#axis,指定对哪一维度去趋势
#type,可设置为'linear'即为去线性趋势,设置为'constant',则为去平均值,即为求距平
#bp,断点,若设置,则为断点两侧分别去趋势,即将序列分成两个子序列各自计算
#overwrite_data,是否覆盖原数据

对上边的序列去线性趋势(红色为原始序列,蓝色为去趋势后序列):

time_series_dtrend = signal.detrend(time_series, axis=0, type='linear', 
                                    bp=0, overwrite_data=False)
plt.plot(x,time_series_dtrend,c = "b")
plt.plot(x,time_series,c = "r")
去趋势对比图

滤波

滤波有很多滤波器,我这里只给出其中一种,仍是scipy.signal提供的Butterworth(巴特沃斯)滤波器。

scipy.signal.butter(N, Wn, btype='low')
#N: 滤波器阶数
#Wn:临界频率
#btype:滤波器类型,{‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’}分别为高通,低通,带通,带阻。
#比如50年序列,提取年代际用低通,提取年际用高通

那比如说我们要对上边的原始序列分别提取年代际分量,我们的采样频率是1年,1个年代际周期是10年,也就是频率为0.1,那么截至频率Wn=2*0.1/1=0.2,比如我们构造一个8阶Butterworth低通滤波器,那么

b, a = signal.butter(8, 0.2, 'lowpass')  
#b,a分别为分子和分母系数,用法请继续往下看

滤波器构造好后就可以对原序列滤波了,这时候使用到滤波函数

scipy.signal.filtfilt(b, a, x, axis=-1)
#这时我们就用到了a和b
#X为要滤波的序列
#axis指定对哪一维滤波

完整的滤波:

b, a = signal.butter(8, 0.2, 'lowpass')  
time_series_filt = signal.filtfilt(b, a, time_series,axis = 0)   
plt.plot(x,time_series_filt,c = "r")
滤波

虽然看起来很奇怪,但是确实是这样,我给出九年滑动平均来对比(蓝色为原始序列,红色为低通滤波,黑色为九年滑动平均):


效果对比

所以这就不难发现为什么很多人用九年滑动平均代替年代际分量,因为确实结果差不多,实际上从原理来讲我认为滑动平均也可以算是一种滤波了。至于为什么滤波的结果这么平滑,这是滤波器的阶数造成的,降低滤波器的阶数就可以了。

滑动平均

最后部分给出滑动平均的方法,前阵子我在一个群里看到有人问,实际上python的滑动平均是被卷积函数代替的。通过构造一个卷积核来实现权重滑动平均,在这一点上,它的用法是比NCL中的smooth函数更广泛的。我们以等权重九年滑动平均举例(卷积的概念就不细讲了,感兴趣的可以百度,大致理解为滑动也是可以的。):

#首先构造一个等权重的卷积核,这个卷积核规定了滑动长度为9,每个点的权重是等权重的1/9
a = np.repeat(1/9, 9)
#卷积运算(滑动平均)
time_series_9 = np.convolve(time_series, a, mode='same')
#这里的mode是设置滑动后序列两端值的 参数有{‘full’, ‘valid’, ‘same’}
#‘full’ 默认值,返回每一个卷积值,长度是N+M-1,在卷积的边缘处,信号不重叠,存在边际效应,实际上是扩充了两端数值,通常不选这个选项。
#‘same’ 返回的数组长度为max(M, N),边际效应依旧存在。得到序列与原始序列长度一样。
#‘valid’  两端缺失,不自动补全。
上一篇下一篇

猜你喜欢

热点阅读