一种基于卡尔曼滤波的QAR数据噪声处理方法

本发明涉及qar数据预处理方法领域，具体涉及一种基于卡尔曼滤波的qar数据噪声处理方法。

背景技术：

利用qar数据可以还原真实的飞行轨迹，但是由于gps误差以及其他各种各样不可预知的影响因素的存在，qar数据中总是存在部分噪声，如果不对轨迹数据加以处理，则可能会出现轨迹“打结”、轨迹“回跳”等现象，严重影响了轨迹的展示效果。而为了尽可能真实地还原飞机的飞行轨迹，提升轨迹展示效果，需要对qar数据进行噪声处理。

卡尔曼滤波是一种利用观测值对预测值进行修正，来得到最优状态估计的算法，其实质是寻找满足均方误差最小条件下的状态估计值。卡尔曼滤波算法在实际应用中仅需要保存上一时刻的状态估计值，而无需保存所有的历史数据，借助系统的状态转移方程，利用上一时刻的状态估计值就可以对当前时刻的状态观测值进行修正，从而得到当前时刻的状态估计值，并以此为一个周期进行迭代，继续对下一时刻的状态值进行估计。卡尔曼递推算法计算的时间复杂度较低，所需要的内存空间较小，非常适用于实时连续系统。

技术实现要素：

本发明的目的在于提供一种基于卡尔曼滤波的qar数据噪声计算方法，该方法能够利用qar数据中记录的丰富信息对飞行轨迹点进行有效的纠正，使纠正后的轨迹点更接近于真实的飞行轨迹点，具体步骤如下：

步骤1，定义t时刻状态向量、t时刻状态转移矩阵、t时刻输入传递矩阵；

步骤2，将t时刻飞机状态向量和t时刻实际状态转移协方差矩阵结合状态转移矩阵、输入传递矩阵进行计算，得到t+1时刻预测状态向量、t+1时刻飞机状态预测误差协方差矩阵、t+1时刻转移状态误差协方差矩阵。

步骤3，将t+1时刻状态向量观测值与t+1时刻观测值向量，根据t+1时刻转移状态误差协方差矩阵得到t+1时刻卡尔曼系数；

步骤4，将于步骤2得到的t+1时刻预测状态向量与观测状态通过t+1时刻卡尔曼系数进行融合得到滤波后的t+1时刻状态，并根据状态误差协方差矩阵得到t+1时刻实际状态转移协方差矩阵。

步骤5，重复步骤2到步骤4，如此循环往复直到所有时刻的状态完成处理，得到完整的航迹数据。

作为优选，步骤1所述t时刻状态向量为：

其中，statet为t时刻的状态向量，t用来指代数据的序号，k代表最大序号，与采样次数即数据长度相等，lngt为t时刻飞机所在的经度，latt为时刻t飞机所在的纬度，vt为时刻t飞机相对于地面的速度，当t＝0时，statet为初始状态向量；

步骤1所述t时刻状态转移矩阵为：

ft表示t时刻的状态转移矩阵，当t＝0时，ft为初始状态转移矩阵。其中，ωx^t表示在t时刻把飞行距离转化成经度需要乘上的系数，ωy表示把飞行距离转化成纬度需要乘上的系数，θt为t时刻飞机的朝向角，δt为qar数据的采样间隔，上述ωx^t，ωy为：

其中，latt表示t时刻飞机所处的纬度，earth-radius_long表示地球长半轴半径，earth_radius_short表示地球短半轴半径。

步骤1所述t时刻输入传递矩阵为：

其中，bt为t时刻的输入传递矩阵，ωx^t表示在t时刻把飞行距离转化成经度需要乘上的系数，ωy表示把飞行距离转化成纬度需要乘上的系数，由于在对轨迹进行采样时δt足够小，根据牛顿运动学定则，此时在做近似匀加速直线运动的状态转移方程为：

xt+1＝xt+(vt*δt+0.5*at*δt²)*sin(θt)*ωx

yt+1＝yt+(vt*δt+0.5*at*δt²)*cos(θt)*ωy

vt+1＝vt+at*δt

其中，at表示t时刻飞机的水平加速度，vt*δt+0.5*at*δt²表示δt内飞机的水平位移，vt为t时刻飞机的水平速度，δt为离散时间的采样间隔；

作为优选，步骤2所述t+1时刻预测状态向量为：

statet+1^-＝ft*statet+bt*at+xt，0≤t≤k-1

式中，ft为t时刻的状态转移矩阵，bt为t时刻的输入传递矩阵，at为t时刻飞机的水平加速度，ft*statet表示飞机在两个时刻间的状态转移过程，bt*at表示水平加速度对飞机状态的影响；xt为t时刻系统状态向量误差，该误差使得飞机动力装置对飞机的控制与记录数据之间存在偏差，符合正态分布xt～n(0，qt)，qt为t时刻状态向量误差的协方差矩阵；

步骤2所述t+1时刻飞机状态预测误差协方差矩阵为：

cov(statet+1)＝cov(ft*statet)＝ftptft^t，0≤t≤k-1

步骤2中所述t+1时刻转移状态误差协方差矩阵为：

pt+1^-＝ftptft^t+qt，0≤t≤k-1

式中，pt+1^-为t+1时刻的转移状态误差协方差矩阵，pt为t时刻的状态误差协方差矩阵，ft为t时刻的状态转移矩阵，qt为t时刻状态向量误差的协方差矩阵；

作为优选，步骤3所述t+1时刻状态向量观测值为：

zt+1＝h*statet+1^-+r，0≤t≤k-1

其中，zt+1为t+1时刻状态向量观测值，h为变换矩阵，r为观测误差协方差矩阵；

步骤3所述t+1时刻观测值向量为：

其中，lngt+1为t+1时刻观测到的飞机所在经度，latt+1为t+1时刻观测到的飞机所在纬度。

步骤3所述t+1时刻卡尔曼系数为：

kt+1＝pt+1^-*h^t*(h*pt+1-*h^t+r)^-1，0≤t≤k-1

其中，pt+1^-为t+1时刻转移状态误差协方差矩阵，h为变换矩阵，r为观测误差协方差矩阵。卡尔曼系数的主要作用是衡量根据状态协方差矩阵pt+1^-和观测状态协方差矩阵r的大小，决定卡尔曼滤波得到的估计值是更偏向于预测值还是更偏向观测值；

作为优选，步骤4所述滤波后的t+1时刻状态为：

statet+1＝statet+1^-+kt+1*(zt+1-h*statet+1^-)，0≤t≤k-1

式中，statet+1^-为t+1时刻的预测状态，zt+1为t+1时刻的观测状态，h为变化矩阵，kt+1为t+1时刻的卡尔曼系数。zt+1-h*statet+1^-为实际的观测与预测值之间的残差，其与kt+1的乘积用来修正飞机位置的预测值。

步骤4所述t+1时刻实际状态转移协方差矩阵为：

pt+1＝(i-kt*h)*pt+1^-，0≤t≤k-1

其中，i为单位矩阵，kt为t时刻的卡尔曼系数，pt+1^-为t+1时刻转移状态误差协方差矩阵，h为变换矩阵；随着迭代的进行，最佳状态估计值的噪声分布情况需要不断地进行更新；状态的不确定性随着迭代的进行总体上会逐渐变小，但是由于传递噪声的存在，状态的不确定性又会呈现出逐渐变大的趋势；卡尔曼滤波算法就是通过不断地调整从而在这一状态中达到一个平衡；

本发明可以纠正由于测量误差导致的qar数据中的空间位置误差，过程中结合了预测值与观测值，经过滤波后的空间位置更具有连贯性，符合实际情况。

附图说明

图1：是本发明的基本方法流程示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

下面结合图1介绍本发明的具体实施方式如下。

本发明实施例为一种基于卡尔曼滤波的qar数据噪声处理方法，具体如下：

步骤1，定义t时刻状态向量、t时刻状态转移矩阵、t时刻输入传递矩阵；

步骤1所述t时刻状态向量为：

其中，statet为t时刻的状态向量，t用来指代数据的序号，k代表最大序号，与采样次数即数据长度相等，lngt为t时刻飞机所在的经度，latt为时刻t飞机所在的纬度，vt为时刻t飞机相对于地面的速度，当t＝0时，statet为初始状态向量，k＝7289；

步骤1所述t时刻状态转移矩阵为：

ft表示t时刻的状态转移矩阵，当t＝0时，ft为初始状态转移矩阵。其中，ωx^t表示在t时刻把飞行距离转化成经度需要乘上的系数，ωy表示把飞行距离转化成纬度需要乘上的系数，θt为t时刻飞机的朝向角，δt为qar数据的采样间隔，k＝7289，上述ωx^t，ωy为：

其中，latt表示t时刻飞机所处的纬度，earth_radius-long表示地球长半轴半径，单位为米，earth_radius_short表示地球短半轴半径，单位为米。

步骤1所述t时刻输入传递矩阵为：

其中，bt为t时刻的输入传递矩阵，ωx^t表示在t时刻把飞行距离转化成经度需要乘上的系数，k＝7289，ωy表示把飞行距离转化成纬度需要乘上的系数，由于在对轨迹进行采样时δt足够小，根据牛顿运动学定则，此时在做近似匀加速直线运动的状态转移方程为：

xt+1＝xt+(vt*δt+0.5*at*δt2)*sin(θt)*ωx

yt+1＝yt+(vt*δt+0.5*at*δt²)*cos(θt)*ωy

vt+1＝vt+at*δt

其中，at表示t时刻飞机的水平加速度，vt*δt+0.5*at*δt²表示δt内飞机的水平位移，vt为t时刻飞机的水平速度，δt为离散时间的采样间隔；

步骤2，将t时刻飞机状态向量和t时刻实际状态转移协方差矩阵结合状态转移矩阵、输入传递矩阵进行计算，得到t+1时刻预测状态向量、t+1时刻飞机状态预测误差协方差矩阵、t+1时刻转移状态误差协方差矩阵。

步骤2所述t+1时刻预测状态向量为：

statet+1^-＝ft*statet+bt*at+xt，0≤t≤k-1

式中，ft为t时刻的状态转移矩阵，bt为t时刻的输入传递矩阵，at为t时刻飞机的水平加速度，ft*statet表示飞机在两个时刻间的状态转移过程，bt*at表示水平加速度对飞机状态的影响；xt为t时刻系统状态向量误差，该误差使得飞机动力装置对飞机的控制与记录数据之间存在偏差，符合正态分布xt～n(0，qt)，qt为t时刻状态向量误差的协方差矩阵，k＝7289；

步骤2所述t+1时刻飞机状态预测误差协方差矩阵为：

cov(statet+1)＝cov(ft*statet)＝ftptft^t，0≤t≤k-1

步骤2中所述t+1时刻转移状态误差协方差矩阵为：

pt+1^-＝ftptft^t+qt，0≤t≤k-1

式中，pt+1^-为t+1时刻的转移状态误差协方差矩阵，pt为t时刻的状态误差协方差矩阵，ft为t时刻的状态转移矩阵，qt为t时刻状态向量误差的协方差矩阵，k＝7289；

步骤3，将t+1时刻状态向量观测值与t+1时刻观测值向量，根据t+1时刻转移状态误差协方差矩阵得到t+1时刻卡尔曼系数；

步骤3所述t+1时刻状态向量观测值为：

zt+1＝h*statet+1^-+r，0≤t≤k-1

其中，zt+1为t+1时刻状态向量观测值，h为变换矩阵，r为观测误差协方差矩阵，k＝7289；

步骤3所述t+1时刻观测值向量为：

其中，lngt+1为t+1时刻观测到的飞机所在经度，latt+1为t+1时刻观测到的飞机所在纬度，k＝7289。

步骤3所述t+1时刻卡尔曼系数为：

kt+1＝pt+1^-*h^t*(h*pt+1^-*h^t+r)^-1，0≤t≤k-1

其中，pt+1^-为t+1时刻转移状态误差协方差矩阵，h为变换矩阵，r为观测误差协方差矩阵，k＝7289。卡尔曼系数的主要作用是衡量根据状态协方差矩阵pt+1^-和观测状态协方差矩阵r的大小，决定卡尔曼滤波得到的估计值是更偏向于预测值还是更偏向观测值；

步骤4，将于步骤2得到的t+1时刻预测状态向量与观测状态通过t+1时刻卡尔曼系数进行融合得到滤波后的t+1时刻状态，并根据状态误差协方差矩阵得到t+1时刻实际状态转移协方差矩阵。

步骤4所述滤波后的t+1时刻状态为：

statet+1＝statet+1^-+kalt+1*(zt+1-h*statet+1^-)，0≤t≤k-1

式中，statet+1^-为t+1时刻的预测状态，zt+1为t+1时刻的观测状态，h为变化矩阵，kalt+1为t+1时刻的卡尔曼系数，zt+1-h*statet+1^-为实际的观测与预测值之间的残差，其与kalt+1的乘积用来修正飞机位置的预测值，k＝7289。

步骤4所述t+1时刻实际状态转移协方差矩阵为：

pt+1＝(i-kt*h)*pt+1^-，0≤t≤k-1

其中，i为单位矩阵，kalt为t时刻的卡尔曼系数，pt+1^-为t+1时刻转移状态误差协方差矩阵，h为变换矩阵，k＝7289；随着迭代的进行，最佳状态估计值的噪声分布情况需要不断地进行更新；状态的不确定性随着迭代的进行总体上会逐渐变小，但是由于传递噪声的存在，状态的不确定性又会呈现出逐渐变大的趋势；卡尔曼滤波算法就是通过不断地调整从而在这一状态中达到一个平衡；

步骤5，重复步骤2到步骤4，如此循环往复直到所有时刻的状态完成处理，得到完整的航迹数据。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。