CN109802656B - 基于幅度信息的卡尔曼滤波方法 - Google Patents

Abstract

本发明公开了一种基于幅度信息的卡尔曼滤波方法，主要解决当观测噪声变化时，现有技术对目标状态估计精度不足的问题，其实现方案是：1)在不同扫描周期，获取雷达信号处理得到的量测信息、幅度信息、状态噪声协方差；2)由雷达系统参数获得初始状态误差协方差和观测噪声协方差；3)将1)和2)获得的参数设置为滤波器的初始值；4)利用获得的幅度信息计算观测噪声的协方差比例系数，并计算观测噪声协方差；5)根据获得的观测噪声协方差进行状态估计；6)将状态估计得到的信息输入到下一扫描周期。本发明减少了观测噪声与观测噪声统计信息失配的情况，提高了目标状态估计的准确率，可用于对雷达目标进行跟踪和参数估计。

Description

基于幅度信息的卡尔曼滤波方法

技术领域

本发明属于雷达技术领域，特别涉及一种卡尔曼滤波方法，可用于对雷达目标进行跟踪和参数估计。

背景技术

卡尔曼滤波算法是目前被广泛应用的一种雷达目标跟踪算法。在系统特性已知、系统噪声以及观测噪声的统计特性都先验已知的情况下，卡尔曼滤波算法能够实现最优估计。但在一般情况下，系统的状态并不都是先验已知的，状态噪声和观测噪声的统计特性并不是固定的，尤其是观测噪声的统计信息会随着目标RCS的起伏发生变化，在这种情况下用传统的卡尔曼滤波算法进行雷达目标跟踪很难实现最优性，容易出现滤波发散、跟踪精度下降以至目标跟踪丢失的问题。

为了提高卡尔曼滤波算法的性能，解决噪声不平稳的问题，实现最优估计,目前采用的方法主要是通过自适应地设计和调整滤波系统的状态噪声协方差来进行，如自适应卡尔曼滤波算法AKF等，这些方法虽然能够根据不同环境自适应的获取恰当的滤波参数，自适应的改变状态噪声的协方差矩阵，但是没有解决观测噪声统计信息变化的问题。一些调整观测噪声的协方差矩阵方法，其中包括两次Kalman滤波的观测噪声自适应调整方法，由于只利用了位置信息没有有效地利用量测除位置以外的其他信息量，因而需要进行多步存储，大大增加了运算量，难以进行工程实现。

在实际情况下，按照雷达目标回波模型，在不同的扫描周期下，目标的RCS是起伏的，这可等效为卡尔曼滤波时观测噪声统计信息在不断变化。而目前的卡尔曼滤波方法均没有涉及目标RCS起伏对观测噪声统计信息的影响，即没有与信号处理的实际情况进行有效结合，难以对观测噪声的统计信息进行有效估计，造成观测噪声与其统计信息失配的情况，影响状态估计以及跟踪精度。

发明内容

本发明目的在于针对上述已有技术的缺陷，提出一种基于幅度信息的卡尔曼滤波方法以提高卡尔曼滤波的状态估计精度，提高雷达目标跟踪的性能。

为实现上述目的，本发明的技术方案包括如下：

1.基于幅度信息的卡尔曼滤波方法，其特征在于，包括如下：

(1)在每次扫描时，获取雷达信号处理得到的量测信息Z、目标幅度信息A和过程噪声协方差信息Q；

(2)通过雷达系统参数获取初始观测噪声协方差信息R(1)和初始状态误差协方差信息P(1|1)；

(3)在第1个扫描周期，用(1)和(2)中获得的信息，设置滤波器的初始值，即第1个扫描周期的状态误差协方差P(1|1)，过程噪声协方差Q(1),观测噪声协方差信息R(1)、幅度信息A(1)及由量测信息Z(1)得到的状态估计信息

(4)在第k+1个扫描周期获取雷达信号处理报送的幅度信息A(k+1)，并利用上一个扫描周期的幅度信息A(k)，计算得到观测噪声的协方差比例系数L(k+1)为：

(5)利用由(4)得到的比例系数L(k+1)和第k个周期的观测噪声协方差R(k)计算第k+1个扫描周期观测噪声协方差R(k+1)：

R(k+1)＝L(k+1)R(k)；

(6)根据(5)得到的观测噪声协方差R(k+1)，通过卡尔曼滤波算法计算第k+1个扫描周期的状态估计信息

(7)将第k+1个扫描周期的状态估计信息

状态误差协方差P(k+1|k+1)，观测噪声协方差R(k+1)、过程噪声协方差Q(k+1)、目标幅度信息A(k+1)作为下一个周期滤波器的输入，进入下一个周期的状态估计。

本发明与现有技术相比，具有如下优点：

(1)本发明利用了雷达信号处理得到的幅度信息，通过计算协方差比例系数，得到观测噪声协方差，实现了对不同扫描周期观测噪声协方差的正确估计，减少了观测噪声与观测噪声协方差失配的情况，解决了随扫描周期变化观测噪声统计信息难以实时计算的难题，使得目标的状态估计更加准确，改善目标跟踪的效率和质量。

(2)本发明由于根据幅度信息获得的观测噪声协方差R(k+1)，通过卡尔曼滤波算法计算第k+1个扫描周期的状态估计信息

减小了计算量，便于工程实现。

附图说明

图1为本发明的实现流程图；

图2为在雷达扫描周期总个数为100时，用本发明方法与传统卡尔曼方法以及真实轨迹、量测轨迹对比图；

图3为在雷达扫描周期总个数为100时，用本发明方法与传统卡尔曼滤波方法仿真在每个扫描周期x方向上的测量误差图；

图4为在雷达扫描周期总个数为100时，用本发明方法与传统卡尔曼滤波方法仿真在每个扫描周期y方向上的测量误差图。

具体实施方式

以下结合附图对本发明进行详细说明：

本发明的基于幅度信息的卡尔曼滤波方法，适用于目标运动状态在离散时间系统中由线性方程表示的场景，该场景下的目标运动状态方程为：

X(k+1)＝F(k)X(k)+G(k)u(k)+V(k)，k＝1，2，…，

其中，k为扫描周期，F(k)为状态转移矩阵，X(k)为离散动态系统在第k个扫描周期的状态；G(k)为输入控制矩阵；u(k)为已知输入或控制信号；V(k)为过程噪声，其协方差为Q(k)。

在该场景下，根据目标运动状态方程X(k+1)，将滤波器接收到的量测信息Z(k+1)表示为：

Z(k+1)＝H(k+1)·X(k+1)+W(k+1)

其中，H(k+1)为量测矩阵，W(k+1)为具有协方差R(k+1)的观测噪声序列，观测噪声W(k)和过程噪声V(k)分别为均值为零，互不相关的白噪声，并满足如下关系：

E(W(k))＝0

Cov(W(k),W(j))＝E(W(k)W^T(j))＝Q_kδ_kj

E(V(k))＝0

Cov(V(k),V(j))＝E(V(k)V^T(j))＝R_kδ_kj

Cov(W(k),V(j))＝E(W(k)V^T(j))＝0

其中，E(W(k))为W(k)求数学期望，E(V(k))为V(k)的数学期望，Cov(W(k),W(j))为W(k)和W(j)的协方差，Cov(V(k),V(j))为V(k)和V(j)的协方差，Cov(W(k),V(j))为W(k)和V(j)的协方差，V^T(j)是V(j)的共轭转置，W^T(j)是W(j)的共轭转置，δ_kj为Kronecher-δ函数，定义为：

参照图1，本发明的实现步骤如下：

步骤1,获取雷达信号处理得到的参数。

滤波器会随扫描周期的更新不断从雷达信号处理的结果中获得目标的各种信息，包括目标的量测信息Z、目标幅度信息A和过程噪声协方差信息Q，其中，量测信息Z包含速度分量s和位置分量b。

步骤2,获得初始状态误差协方差P(1|1)和初始观测噪声协方差R(1)。

从雷达系统参数中获取初始观测噪声协方差信息R(1)和初始状态误差协方差信息P(1|1)，其中：

初始状态误差协方差P(1|1)与常数r有关，r是从雷达系统参数直接获取的常数，根据r得到初始状态误差协方差P(1|1)为：

T为扫描周期的时间间隔；

初始观测噪声协方差R(1)由第1个扫描周期的测距误差

测角误差

组成，

和

直接从雷达系统参数中获取,表示如下：

步骤3,在第1个扫描周期，设置滤波器的初始值。

(3a)由第1个扫描周期获取的量测信息得到的状态估计信息

将第1个扫描周期的量测信息Z(1)表示为：

其中，b(1)是Z(1)的位置分量，s(1)是Z(1)的速度分量；

根据b(1)和s(1)得到第1个周期的状态估计信息

其中，

和

分别表示状态估计信息

的位置分量和速度分量；

(3b)将第1个扫描周期的状态误差协方差P(1|1)，过程噪声协方差Q(1),观测噪声协方差信息R(1)、幅度信息A(1)及由量测信息Z(1)得到的状态估计信息

设置为滤波器的初始值。

步骤4,利用幅度信息，计算观测噪声协方差比例系数L(k+1)。

在第k+1个扫描周期滤波器获取雷达信号处理端报送的各种信息后，根据本周期的幅度信息A(k+1)，上一个扫描周期的幅度信息A(k)，计算得到观测噪声的协方差比例系数L(k+1)：

步骤5,通过比例系数L(k+1)计算观测噪声协方差R(k+1)。

(5a)在第k个扫描周期，将观测噪声协方差R(k)由信号处理阶段的测角误差σ_θ(k)和测距误差计σ_ρ(k)表示如下:

(5b)将第k个扫描周期的测角误差σ_θ(k)表示为：

其中，U₁为常数，其与测角方法有关；Q_BW为波束宽度，SNR_k为第k个扫描周期目标点迹的检测信噪比；

将第k个扫描周期的测距误差σ_ρ(k)表示为：

其中，U₂为常数，其与雷达系统参数及性能有关；D为雷达的距离分辨率，其与信号的有效带宽有关，SNR_k为第k个扫描周期目标的检测信噪比；

(5c)计算目标的检测信噪比SNR_k：

其中，A(k)为第k个扫描周期的目标点迹的幅度值，n²为接收机噪声功率，在一个雷达系统中，由于接收机噪声功率n²为定值，所以检测信噪比SNR_k只与目标幅度信息A有关；

(5d)将信噪比SNR_k带入到(5b)公式中，得第k个扫描周期中观测噪声协方差R(k)为：

(5e)根据(5d)的结果，得到第k+1个扫描周期的观测噪声的协方差矩阵R(k+1)为：

步骤6,通过观测噪声协方差R(k+1)，进行第k+1个扫描周期的状态估计。

在现有技术中，进行目标状态估计的方法主要有最小二乘法算法、α-β滤波算法、粒子滤波算法以及卡尔曼滤波算法等，本实例采用卡尔曼滤波算法，进行目标的状态估计。

目标状态估计是一个不断迭代的过程，需要利用上一个周期的状态估计信息

状态误差协方差P、过程噪声协方差Q进行本周期的状态估计。在完成第k个扫描周期的目标状态估计后，第k个扫描周期的目标状态信息

状态误差协方差P(k|k)，过程噪声协方差Q(k)成为已知量，这些信息成为第k+1个扫描周期滤波器的输入，用于第k+1个扫描周期目标的状态估计。

本步骤的具体实现如下:

(6a)在第k+1个扫描周期进行状态估计前，滤波器获取第k个扫描周期状态估计信息

状态误差协方差P(k|k)，过程噪声协方差Q(k)；

(6b)利用第k个扫描周期的状态估计信息

进行k+1时刻的状态一步预测：

其中，

为由第k+1个扫描周期的状态预测值，F(k)为状态转移矩阵，G(k)为输入控制矩阵，u(k)为已知输入或控制信号；

(6c)由第k+1个扫描周期的状态预测值

得到量测预测值

其中H(k+1)为量测矩阵；

(6d)由第k+1个扫描周期量测信息Z(k+1)和量测预测值

计算得到新息v(k+1)：

(6e)根据第k个扫描周期的过程噪声协方差Q(k)和状态误差协方差P(k|k)进行第k+1个扫描周期状态误差协方差的一步预测P(k+1|k)：

P(k+1|k)＝F(k)P(k|k)F'(k)+Q(k)，

其中，F'(k)是状态转移矩阵F(k)的转置；

(6f)根据协方差一步预测P(k+1|k)和由(5)得到的观测噪声的协方差R(k+1)，得到第k+1个扫描周期的新息协方差S(k+1)：

S(k+1)＝H(k+1)P(k+1|k)H'(k+1)+R(k+1)，

其中，H'(k+1)是量测矩阵H(k+1)的转置；

(6g)由新息协方差S(k+1)计算得到第k+1个扫描周期的增益矩阵C(k+1)：

C(k+1)＝P(k+1|k)H'(k+1)S^-1(k+1)，

其中，S^-1(k+1)是新息协方差S(k+1)的逆；

(6h)根据增益矩阵C(k+1)和新息协方差S(k+1)，按如下公式进行第k+1个扫描周期的状态误差协方差P(k+1|k+1)的更新：

P(k+1|k+1)＝P(k+1|k)-C(k+1)S(k+1)C'(k+1)，

其中，C'(k+1)为增益矩阵C(k+1)的转置；

(6i)根据新息v(k+1)和增益矩阵C(k+1)，得到第k+1个扫描周期的状态估计信息

步骤7,将第k+1个扫描周期滤波器的信息输入到下一个扫描周期。

将第k+1个扫描周期的状态估计信息

状态误差协方差P(k+1|k+1)，观测噪声协方差R(k+1)、过程噪声协方差Q(k+1)、目标幅度信息A(k+1)作为下一个周期滤波器的输入，用于下一个周期的状态估计。

本发明的效果可以通过以下计算机仿真进一步说明：

一、仿真条件

设定雷达进行单目标跟踪，扫描周期总个数为100，扫描周期时间间隔为1s，目标起始坐标分别为x＝1000m,y＝2000m,x坐标和y坐标上的速度分别为200m/s、5m/s，前70个扫描周期为匀速直线运动，后30个扫描周期为x和y方位上加速度为10m/s²和-7m/s²的匀加速直线运动过程噪声方差为0.1m²，观测噪声是均值为0m，观测噪声径向距离的方差在[15²,150²]区间内的随机变化，观测噪声角度的方差在[0.01²,0.1²]区间内随机变化，幅度在[1,10]区间随机变化，各扫描周期内观测噪声相互独立，状态转移矩阵和观测矩阵分别为：F＝[1,T,0,0；0,1,0,0；0,0,1,T；0,0,0,1]、H1＝[1,0,0,0；0,0,1,0]；

二、仿真内容

仿真1：在上述仿真条件下，用本发明估计的目标运动轨迹，并与观测样本轨迹和真实轨迹对比，结果如图2所示，从图2中可以看出本发明估计的轨迹接近真实轨迹，减少了观测噪声的影响。

仿真2：在进行200次蒙特卡洛试验后，对比本发明实现的滤波算法与传统的卡尔曼滤波算法在直角坐标系下X方向和Y方向上的均方根误差，结果分别如图3、图4所示，从图3和图4中可以看出，在目标进行前70个扫描周期匀速运动和后30个扫描周期匀加速运动状态下，本发明方法进行目标状态估计得到的均方根误差均要小于传统算法，这说明本发明方法有更好的状态估计精度。

Claims (2)

1.基于幅度信息的卡尔曼滤波方法，其特征在于，包括如下：

(1)在每次扫描时，获取雷达信号处理得到的量测信息Z、目标幅度信息A和过程噪声协方差信息Q；

(2)通过雷达系统参数获取初始观测噪声协方差信息R(1)和初始状态误差协方差信息P(1|1)；

(3)在第1个扫描周期，用(1)和(2)中获得的信息，设置滤波器的初始值，即第1个扫描周期的状态误差协方差信息P(1|1)，过程噪声协方差信息Q(1),观测噪声协方差信息R(1)、目标幅度信息A(1)及由量测信息Z(1)得到的状态估计信息

(4)在第k+1个扫描周期获取雷达信号处理报送的目标幅度信息A(k+1)，并利用上一个扫描周期的目标幅度信息A(k)，计算得到观测噪声的协方差比例系数L(k+1)为：

(5)利用由(4)得到观测噪声的协方差比例系数L(k+1)和第k个周期的观测噪声协方差信息R(k)计算第k+1个扫描周期观测噪声协方差信息R(k+1)：

(5a)在第k个扫描周期，将观测噪声协方差信息R(k)由信号处理阶段的测角误差σ_θ(k)和测距误差计σ_ρ(k)表示如下:

(5b)将第k个扫描周期的测角误差σ_θ(k)表示为：

其中，U₁为常数，其与测角方法有关；Q_BW为波束宽度，SNR_k为第k个扫描周期目标点迹的检测信噪比；

将第k个扫描周期的测距误差σ_ρ(k)表示为：

其中，U₂为常数，其与雷达系统参数及性能有关；D为雷达的距离分辨率，其与信号的有效带宽有关，SNR_k为第k个扫描周期目标的检测信噪比；

(5c)计算目标的检测信噪比SNR_k：

其中，A(k)为第k个扫描周期的目标点迹的目标幅度信息，n²为接收机噪声功率，在一个雷达系统中，由于接收机噪声功率n²为定值，所以检测信噪比SNR_k只与目标幅度信息A有关；

(5d)将信噪比SNR_k带入到(5b)公式中，得第k个扫描周期中观测噪声协方差信息R(k)为：

(5e)根据(5d)的结果，得到第k+1个扫描周期的观测噪声协方差信息R(k+1)的矩阵为：

(6)根据(5)得到的观测噪声协方差信息R(k+1)，通过卡尔曼滤波算法计算第k+1个扫描周期的状态估计信息

(7)将第k+1个扫描周期的状态估计信息

状态误差协方差信息P(k+1|k+1)，观测噪声协方差信息R(k+1)、过程噪声协方差信息Q(k+1)、目标幅度信息A(k+1)作为下一个周期滤波器的输入，进入下一个周期的状态估计。

2.根据权利要求1所述的方法,其中(6)中通过卡尔曼滤波算法计算第k+1个扫描周期的状态估计信息

其实现如下:

(6a)利用第k个扫描周期的状态估计信息

进行状态的一步预测：

其中，

为k+1时刻的状态预测值，F(k)为状态转移矩阵，G(k)为输入控制矩阵，u(k)为已知输入或控制信号；

(6b)由状态预测值

得到量测预测值

其中H(k+1)为量测矩阵；

(6c)由量测信息Z(k+1)和量测预测值

计算得到新息v(k+1)：

其中，Z(k+1)为第k+1次扫描获取的量测信息；

(6d)根据第k个扫描周期的过程噪声协方差信息Q(k)和状态误差协方差信息P(k|k)进行状态误差协方差信息的一步预测P(k+1|k)：

P(k+1|k)＝F(k)P(k|k)F'(k)+Q(k)，

其中，F'(k)是状态转移矩阵F(k)的转置；

(6e)根据一步预测P(k+1|k)和(5)得到的观测噪声协方差信息R(k+1)，得到新息协方差S(k+1)：

S(k+1)＝H(k+1)P(k+1|k)H'(k+1)+R(k+1)，

其中，H'(k+1)是量测矩阵H(k+1)的转置；

(6f)由新息协方差S(k+1)计算得到增益矩阵C(k+1)：

C(k+1)＝P(k+1|k)H'(k+1)S^-1(k+1)，

其中，S^-1(k+1)是新息协方差S(k+1)的逆；

(6g)根据增益矩阵C(k+1)和新息协方差S(k+1)按如下公式进行状态误差协方差信息P(k+1|k+1)的更新：

P(k+1|k+1)＝P(k+1|k)-C(k+1)S(k+1)C'(k+1)

其中，C'(k+1)为增益矩阵C(k+1)的转置；

(6h)根据新息v(k+1)和增益矩阵C(k+1)，得到状态估计信息

Images