锐单电子商城 , 一站式电子元器件采购平台!
  • 电话:400-990-0325

两个传感器卡尔曼滤波python优化实现

时间:2023-02-28 14:30:00 传感器kinw8486a传感器t3419p传感器

两个传感器卡尔曼滤波器python优化实现

  • 两个传感器滤波器
  • 优化:提前计算Kk
  • Python实现

两个传感器滤波器

一辆车做变速直线运动,车上有GPS定位仪和速度表。两个传感器的数据通过卡尔曼滤波器优化了汽车的位置。 如果观察到的数据是向量z=[z1 z2].T。 z1为GPS定位位置,z二是速度表的速度读数。z协方差矩阵R~N(0,R)。 车辆状态方程的噪声是wk~N(0,Q)。 系统状态矩阵A=[1]预测值为上一位,最优估计为预测值,加测值与预测值差乘以比例。 卡尔曼滤波计算公式如下

在这里插入图片描述

优化:提前计算Kk

关键要计算k,就是测量差的分配比例,通过(1.11) (1.10) (1.13)可以计算三种类型k,k只和系统噪声Q
与测量误差协方差矩阵R有关。因此,可以提前迭代几次k,然而,估计大量数据的计算位置。只需计算(1.9)式和(1.12)式。k提前计算。将卡尔曼五步简化为两步,大大降低矩阵计算量。

Python实现

如果车辆直线位置的真正方程是100 *(1 - np.exp(-t/10)), 传感器z1的方差为2,速度表z2的方差为0.2,
R=[[2 0] ,系统噪声方差Q=[1]。 系统状态A=[1],观察矩阵H=[1 1].T
[0 0.2]]

函数pre_kalman,通过参数A H Q R ,卡尔曼滤波参数可以提前计算K_K。迭代几次收敛。
函数obv()按方程100 *(1 - np.exp(-t/10)生成汽车的真实位置,加入高斯分布噪声,返回传感器观察
的数据z。z0为GPS位置数据,z一是速度。需要正确z1获得积分的位置。然后根据z和k_k估计卡尔曼滤波器。

一般GPS干扰位置方差大,速度表面差小,但积分后累积误差大。从图中可以看出,卡尔曼滤波器的估计值非常好。如图所示

附程序

import numpy as np import matplotlib import matplotlib.pyplot as plt from matplotlib.font_manager import FontProperties  # 预设字体 font = {'family': 'SimHei', "size": 14} matplotlib.rc('font', **font)  ''' 卡尔曼滤波,2路观察 '''  ''' 真实 x_k = A*xk_1   w_k  w_k方差Q 观察 z_k = H*x_k   v_k   v_k方差R  预测 x_k_predict = A*x_k_est p_k_predict = A * p_k_est * A_T   Q k_k = p_k_predict * H_T *( H*p_k_predict*H_T   R)^-1  最优估计 x_k_est = x_k_predict   k_k *( z_k - H*x_k_predict) p_k_est = (I - k_k * H) * p_k_predict  '''  # 系统状态A A = np.mat([1])  # 两个观察仪器观察矩阵 H = np.mat([[1], [1]])  #观察仪器方差sigma1 sigma2 sigma1 = 2 sigma2 = 0.2  #系统噪声wk方差Q Q = np.mat(1)  '''由Q R算kk''' def pre_kalman(A,H):     R=np.mat([[sigma1, 0],[0, sigma2]])     I=np.mat([1])     p_est = np.mat([0]) #估计协方差矩阵     for _ in np.arange(0,10,1):         p_predict = A * p_est * A.T   Q #预测协方差矩阵         kk = p_predict * H.T * (H * p_predict * H.T   R).I         p_est = (I - kk * H) * p_predict     return kk #返回k_k  ''' 两个测量仪器测量两组信号,z位置测量值为0,z1.速度测量值 z0的方差为sigma1,z1方差为sigma2 真实信号为1-exp(t/10) ''' def obv():     t = np.linspace(1,100,100)     real = 100 *(1 -  np.exp(-t/10))     #real = t**2      #求速度xdot     dt = t[1:] - t[:-1]     dx = real[1:] - real[:-1]     xdot = dx / dt     #补第0个对齐     xdot = np.hstack([real[0],xdot])      #真实信号真实噪声,模拟测量仪器的误差噪声     t = t[:-1]     real=real[:-1]     xdot=xdot[:-1]     z1= real   np.random.normal(0,sigma1,size=t.shape[0])     z2= xdot   np.random.normal(0,sigma2,size=t.shape[0])     z = np.mat([(a,b) for (a,b) in zip(z1,z2)])     return [t, dt, real,z]   ''' 计算卡尔曼 z是观察到的数据 kk是卡尔的慢比例系数 ''' def calc_kalman(A,H,z_k,kk):     global  x_est     x_pred = A * x_est     x_est = x_pred   kk * (z_k - H * x_pred)     return x_est   global x_est #最佳估计值 if __name__ == '__main__':      #算k值     kk = pre_kalman(A,H)      #生成观察仪器信息Z,返回时间,真实值,两个观察仪器的噪声值     [t, dt, xreal, z] = obv()      x_est=np.mat([0])     x_est_list=[]      z1_integral=0     for k in range(z.shape[0]):         #积分速度         z1_integral  = z[k][0,1] * dt[k]         z[k][0,1] = z1_integral          #计算卡尔曼,结果是x_est         calc_kalman(A, H, z[k].T, kk)         x_est_list.append(x_est[0,0])      x_est_list=np.array(x_est_list)       plt.title(卡尔曼滤波器)     plt.plot(t, xreal, label='真实',color='g')     plt.plot(t, z[:,0], label='位置测量',color='b')     plt.plot(t, z[:,1], label=速度积分,color='y')     plt.plot(t, x_est_list, label='估计',color='r')     plt.legend(loc="lower right")     plt.show()   
锐单商城拥有海量元器件数据手册IC替代型号,打造电子元器件IC百科大全!

相关文章