本文共 10119 字,大约阅读时间需要 33 分钟。
在实际的工程应用中,经常会遇到初始结果噪声太多的问题,比如信号强度抖动的太厉害,比如视频流中的bbox抖动的太厉害,比如光谱信号抖动的太厉害等等,这时候就需要一些简单的滑动平均算法。滑动平均其实是一个很朴素的方法,但是要与实际结合,构造出合适的平滑方式,是需要一些思考的。下面我将分别介绍滑动平均法(Mean Average)、指数滑动平均法(Exponential Mean Average)、SG滤波法(Savitzky Golay Filter)。
简单来说,滑动平均法把前后时刻的一共2n+1个观测值做平均,得到当前时刻的滤波结果。这是一个比较符合直觉的平滑方法,在生活中、工作中很经常会用到,但是很少去思考这么做的依据是什么,下面我就来仔细分析一下其中的原理。
对于一个观测序列,我们有这样的假设:每一次的观测值是带有噪声的,而我们期望噪声的均值为0,方差为 σ 2 \sigma^{2} σ2,观测值和真实值之间的关系如下: g t = x t + ε t ( 1 ) g_{t}=x_{t}+\varepsilon_{t}\space\space\space\space\space\space(1) gt=xt+εt (1)其中, x t x_{t} xt为观测值, g t g_{t} gt为真实值, ε t \varepsilon_{t} εt为噪声。为了降低噪声的影响,我们把相邻时刻的观测值相加后平均,公式如下: p t = ∑ i = 1 n ( x t − i + x t + i ) + x t 2 n + 1 ( 2 ) p_{t}=\frac{\sum_{i=1}^{n}{(x_{t-i}+x_{t+i})+x_{t}}}{2n+1}\space\space\space\space\space\space(2) pt=2n+1∑i=1n(xt−i+xt+i)+xt (2) p t p_{t} pt表示 t 时刻的滤波结果, x t − 1 x_{t-1} xt−1表示 t-1 时刻的观测值, n 代表滑动窗口半径。将公式(1)代入公式(2),可以得到 p t = ∑ i = 1 n ( g t − i + g t + i ) + g t − ∑ i = 1 n ( ε t − i + ε t + i ) − ε t 2 n + 1 ( 3 ) p_{t}=\frac{\sum_{i=1}^{n}{(g_{t-i}+g_{t+i})}+g_{t}-\sum_{i=1}^{n}{(\varepsilon_{t-i}+\varepsilon_{t+i})}-\varepsilon_{t}}{2n+1}\space\space\space\space\space\space(3) pt=2n+1∑i=1n(gt−i+gt+i)+gt−∑i=1n(εt−i+εt+i)−εt (3)前面说到了,我们假设噪声的均值为0,所以 ∑ i = 1 n ( ε t − i + ε t + i ) + ε t \sum_{i=1}^{n}{(\varepsilon_{t-i}+\varepsilon_{t+i})+\varepsilon_{t}} ∑i=1n(εt−i+εt+i)+εt为0,那么我们得到的结果就是: p t = ∑ i = 1 n ( g t − i + g t + i ) + g t 2 n + 1 ( 4 ) p_{t}=\frac{\sum_{i=1}^{n}{(g_{t-i}+g_{t+i})}+g_{t}}{2n+1}\space\space\space\space\space\space (4) pt=2n+1∑i=1n(gt−i+gt+i)+gt (4)当观测数据的真实值变化较小时,或者变化为线性时,可以近似认为: p t = ∑ i = 1 n ( g t − i + g t + i ) + g t 2 n + 1 = g t ( 5 ) p_{t}=\frac{\sum_{i=1}^{n}{(g_{t-i}+g_{t+i})}+g_{t}}{2n+1}=g_{t}\space\space\space\space\space\space (5) pt=2n+1∑i=1n(gt−i+gt+i)+gt=gt (5)从上面的分析过程我们可以看到,当滑动窗口内的真实数据变化不大的时候,我们可以抑制掉很大一部分噪声,滤波结果近似真实值;当滑动窗口内的真实值变化较大时,这种滤波方式就会损失一部分精确度,滤波结果接近真实值的平均期望。所以,窗口的大小会对滤波结果有很大影响。窗口越大,滤波结果越平滑,但会一定程度上偏离真实值;窗口越小,滤波结果越接近观测值,但噪声偏大。滑动平均法还有一个升级版本,也就是加权滑动平均法。实际场景中,每个观测值的重要程度不同,忽略每个观测值的置信度直接平均不能得到精确的结果,所以就需要给观测值加权。加权滑动平均法的公式如下:
p t = ∑ i = 1 n ( x t − i ∗ w t − i + x t + i ∗ w t + i ) + x t ∗ w t 2 n + 1 ( 6 ) p_{t}=\frac{\sum_{i=1}^{n}{(x_{t-i}*w_{t-i}+x_{t+i}*w_{t+i})+x_{t}*w_{t}}}{2n+1} \space\space\space\space\space\space (6) pt=2n+1∑i=1n(xt−i∗wt−i+xt+i∗wt+i)+xt∗wt (6) w t w_{t} wt为 t 时刻的权重。(6)式表示的是把每个观测值乘以权重后再平均。这种方法适用于观测值本身带有置信度的情况。注意,这里有一个小问题,如果置信度的取值范围是0到1之间,那么加权之后计算得到的观测值往往小于真实值,我来解释一下为什么。首先,我们假设观测值和真实值的均值是相等的,也就是 X ˉ = ∑ i = 1 2 n + 1 x t 2 n + 1 = ∑ i = 1 2 n + 1 g t 2 n + 1 = G ˉ ( 7 ) \bar{X} = \frac{\sum_{i=1}^{2n+1}{x_{t}}}{2n+1}=\frac{\sum_{i=1}^{2n+1}{g_{t}}}{2n+1} = \bar{G} \space\space\space\space\space\space (7) Xˉ=2n+1∑i=12n+1xt=2n+1∑i=12n+1gt=Gˉ (7)当我们把观测值乘以权重了之后,观测值和真实值的均值就不相等了,因为真实值的权重均值为1,而观测值的权重均值为 w t ˉ = ∑ i = 1 n ( w t − i + w t + i ) + w t 2 n + 1 \bar{w_{t}}=\sum_{i=1}^{n}\frac{(w_{t-i}+w_{t+i})+w_{t}}{2n+1} wtˉ=∑i=1n2n+1(wt−i+wt+i)+wt,是小于等于1的,最终的预测值也是小于等于真实值的,而且大概率是小于。所以我们需要对(6)式增加一个修正: p t ~ = p t w t ˉ ( 8 ) \tilde{p_{t}} = \frac{p_{t}}{ \bar{ w_{t}} } \space\space\space\space\space\space (8) pt~=wtˉpt (8)这样,得到的预测值就会更加合理了。小结:滑动平均法使用的前提是,噪声的均值为0,真实值变化不大或线性变化的场景。如果真实值有较高频率的非线性突变的话,滑动平均法的效果就不够好了。同时,滑动平均法的窗口选取很重要,需要根据具体数据来选择。如果需要使用在线版本的滑动平均,那么就要把窗口前移,也就是把当前时刻的前n个观测值进行平均,但这样得到的结果会明显滞后于当前观测值,窗口越大,滞后的现象越严重。
class MovAvg(object): def __init__(self, window_size=7): self.window_size = window_size self.data_queue = [] def update(self, data): if len(self.data_queue) == self.window_size: del self.data_queue[0] self.data_queue.append(data) return sum(self.data_queue)/len(self.data_queue)
指数滑动平均法相当于加权滑动平均法的变体,主要区别在于,指数滑动平均法的权重是固定的、随时间推移呈指数衰减。指数滑动平均法的公式如下:
p t = w ∗ x t + ( 1 − w ) ∗ p t − 1 ( 9 ) p_{t}=w*x_{t}+(1-w)*p_{t-1} \space\space\space\space\space\space (9) pt=w∗xt+(1−w)∗pt−1 (9) p t p_{t} pt表示预测值, w w w表示衰减权重,通常我们设为固定值0.9, x t x_{t} xt表示观测值,这是一个递推公式。前面说了,指数滑动平均法的权重是随时间推移呈指数衰减的,那么上面的这个递推公式的指数体现在哪里呢?我们把(9)式进行延伸: p t − 1 = w ∗ x t − 1 + ( 1 − w ) ∗ p t − 2 ( 10 ) p_{t-1}=w*x_{t-1}+(1-w)*p_{t-2} \space\space\space\space\space\space (10) pt−1=w∗xt−1+(1−w)∗pt−2 (10)将(9)和(10)两式子联立,可得 p t = w ∗ x t + ( 1 − w ) ∗ ( w ∗ x t − 1 + ( 1 − w ) ∗ p t − 2 ) ( 11 ) p_{t}=w*x_{t}+(1-w)*(w*x_{t-1}+(1-w)*p_{t-2}) \space\space\space\space\space\space (11) pt=w∗xt+(1−w)∗(w∗xt−1+(1−w)∗pt−2) (11)发现没有,在(11)式中 p t p_{t} pt与 p t − 2 p_{t-2} pt−2的关系是 ( 1 − w ) 2 (1-w)^{2} (1−w)2倍,而在(9)式中 p t p_{t} pt与 p t − 1 p_{t-1} pt−1的关系是 1 − w 1-w 1−w倍,呈指数衰减关系。同时,在初始时刻有如下关系: p 0 = x 0 ( 12 ) p_{0}=x_{0} \space\space\space\space\space\space (12) p0=x0 (12)根据这一关系和上述的递推公式,我们就能够得到整个算法的公式了。由于这种指数衰减的特性,指数滑动平均法会比滑动平均法的实时性更强,更加接近当前时刻的观测值。在实际场景下,如果目标的波动较大时,指数滑动平均法会比滑动平均更加接近当前的真实值。那么是不是就说明,指数滑动平均法在任意场景下都比滑动平均法更好呢?不一定。我们来分析一下指数衰减法的误差项,这里为了简便表示,设定 t=2 ,同时,将(1)式和(12)式代入公式(11),可得到误差项:
ε = w ∗ ε 2 + ( 1 − w ) ∗ ( w ∗ ε 1 + ( 1 − w ) ∗ ε 0 ) ( 13 ) \varepsilon=w*\varepsilon_{2}+(1-w)*(w*\varepsilon_{1}+(1-w)*\varepsilon_{0}) \space\space\space\space\space\space (13) ε=w∗ε2+(1−w)∗(w∗ε1+(1−w)∗ε0) (13)所以误差项也是呈指数衰减的,越接近当前时刻的误差项权重越大。假如在当前的工程场景中,误差是固定的分布,不受目标的观测值大小影响的话,那么指数滑动平均法会更接近真实值;假如误差会受目标观测值影响,比如我们观测的是一个连续运动的目标,中间突然出现了一个偏离很远的观测点,那么这个点为误检的概率相当大,也就是该观测值的误差比之前其他点的误差要大得多,此时指数加权平均法的结果就会波动较大,结果就不如滑动平均了。小结:当误差不受观测值大小影响的话,指数滑动平均比滑动平均好;当误差随观测值大小变化时,滑动平均比指数滑动平均更好。
class ExpMovAvg(object): def __init__(self, decay=0.9): self.shadow = 0 self.decay = decay self.first_time = True def update(self, data): if self.first_time: self.shadow = data self.first_time = False return data else: self.shadow = self.decay*self.shadow+(1-self.decay)*data return self.shadow
SG滤波法(Savitzky Golay Filter)的核心思想也是对窗口内的数据进行加权滤波,但是它的加权权重是对给定的高阶多项式进行最小二乘拟合得到。它的优点在于,在滤波平滑的同时,能够更有效地保留信号的变化信息,下面我来介绍一下其原理。
我们同样对当前时刻的前后一共2n+1个观测值进行滤波,用k-1阶多项式对其进行拟合。对于当前时刻的观测值,我们用下面的公式进行拟合: x t = a 0 + a 1 ∗ t + a 2 ∗ t 2 + . . . + a k − 1 ∗ t k − 1 ( 14 ) x_{t}=a_{0}+a_{1}*t+a_{2}*t^{2}+...+a_{k-1}*t^{k-1} \space\space\space\space\space\space (14) xt=a0+a1∗t+a2∗t2+...+ak−1∗tk−1 (14)同样,对于前后时刻(如t-1、t+1、t-2、t+2等时刻)的预测值,我们同样可以用(14)式来计算,这样一共得到2n+1个式子,构成一个矩阵(似乎发不了矩阵,我放个图片吧):class SavGol(object): def __init__(self, window_size=11, rank=2): assert window_size % 2 == 1 self.window_size = window_size self.rank = rank self.size = int((self.window_size - 1) / 2) self.mm = self.create_matrix(self.size) self.data_seq = [] def create_matrix(self, size): line_seq = np.linspace(-size, size, 2*size+1) rank_seqs = [line_seq**j for j in range(self.rank)] rank_seqs = np.mat(rank_seqs) kernel = (rank_seqs.T * (rank_seqs * rank_seqs.T).I) * rank_seqs mm = kernel[self.size].T return mm def update(self, data): self.data_seq.append(data) if len(self.data_seq) > self.window_size: del self.data_seq[0] padded_data = self.data_seq.copy() if len(padded_data) < self.window_size: left = int((self.window_size-len(padded_data))/2) right = self.window_size-len(padded_data)-left for i in range(left): padded_data.insert(0, padded_data[0]) for i in range(right): padded_data.insert( len(padded_data), padded_data[len(padded_data)-1]) return (np.mat(padded_data)*self.mm).item()
本文图片制作的相关代码。
import cv2import numpy as npimport matplotlib.pyplot as pltimport imageio# 一维数据滤波ma, ema, sg = MovAvg(), ExpMovAvg(), SavGol()data_list, data_ma, data_ema, data_sg = [], [], [], []for i in range(200): data = i+np.random.randint(-50, 50) data_list.append(data) data_ma.append(ma.update(data)) data_ema.append(ema.update(data)) data_sg.append(sg.update(data))plt.plot(data_list, label='raw')plt.plot(data_ma, label='ma')plt.plot(data_ema, label='ema')plt.plot(data_sg, label='sg')plt.legend()plt.show()# 鼠标轨迹滤波ma_x, ma_y = MovAvg(), MovAvg()ema_x, ema_y = ExpMovAvg(), ExpMovAvg()sg_x, sg_y = SavGol(), SavGol()def draw_circle(event, x, y, flags, param): if event == cv2.EVENT_MOUSEMOVE: sx = np.random.randint(-50, 51) sy = np.random.randint(-50, 51) cv2.circle(show, (x+sx, y+sy), 5, (255, 255, 255), -1) x, y = ma_x.update(x+sx), ma_y.update(y+sy) cv2.circle(show, (int(x), int(y)), 5, (0, 0, 255), -1) x, y = ema_x.update(x+sx), ema_y.update(y+sy) cv2.circle(show, (int(x), int(y)), 5, (0, 255, 0), -1) x, y = sg_x.update(x+sx), sg_y.update(y+sy) cv2.circle(show, (int(x), int(y)), 5, (255, 0, 0), -1)show = np.zeros((1024, 1024, 3), np.uint8)cv2.namedWindow('image')buff = []while True: cv2.setMouseCallback('image', draw_circle) cv2.imshow('image', show) save = cv2.resize(show, (512, 512)) save = cv2.cvtColor(save, cv2.COLOR_BGR2RGB) buff.append(save) if cv2.waitKey(100) == ord('q'): breakcv2.destroyAllWindows()imageio.mimwrite('test.gif', buff, 'GIF', duration=0.1)
转载地址:http://xjqqd.baihongyu.com/