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

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


原始序列及其趋势

去趋势

研究长时间的气候变化会发现很多序列是带有趋势的,有部分观点认为这是认为因素造成全球变暖导致的,研究本身的大气变化需要去掉序列的趋势,简单来说就是序列的每一个点都减去序列趋势的回归斜率乘以一个时间步长,当然前提是这个序列的线性趋势是显著的。还有一个要注意的是去趋势序列(若干年长度的序列)并不等于该序列的年际变化,获得年际变化分量需要去用滤波或者其他方法单独提取,我目前看到的大部分文献中的去趋势都是为了剔除人为因素造成的全球变暖。
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’  两端缺失,不自动补全。
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 199,711评论 5 468
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 83,932评论 2 376
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 146,770评论 0 330
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 53,799评论 1 271
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,697评论 5 359
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,069评论 1 276
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,535评论 3 390
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,200评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,353评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,290评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,331评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,020评论 3 315
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,610评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,694评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,927评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,330评论 2 346
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 41,904评论 2 341

推荐阅读更多精彩内容