CN111596290B

CN111596290B - 一种基于最大相关熵扩展卡尔曼滤波的雷达目标跟踪方法

Info

Publication number: CN111596290B
Application number: CN202010483616.4A
Authority: CN
Inventors: 李春霞; 李明星; 王恒; 张德; 王翾
Original assignee: CETC Information Science Research Institute
Current assignee: CETC Information Science Research Institute
Priority date: 2020-06-01
Filing date: 2020-06-01
Publication date: 2022-07-29
Anticipated expiration: 2040-06-01
Also published as: CN111596290A

Abstract

本发明提供一种基于最大相关熵扩展卡尔曼滤波的雷达目标跟踪方法，包括：构建雷达状态方程和雷达量测方程、雷达状态先验更新方程和雷达量测先验更新方程；根据雷达状态方程、雷达状态先验更新方程和雷达量测方程构建雷达系统非线性模型；根据k+1时刻雷达目标点迹的量测预测值和k+1时刻雷达目标点迹的量测值得到k+1时刻核函数对角阵奇异参数；若k+1时刻核函数对角阵奇异参数大于预设的第一阈值，则将k+1时刻雷达目标点迹的状态预测值作为k+1时刻的雷达滤波值。本发明提出的方法对于非线性、非高斯滤波问题，可以获得误差二阶项信息，还可以捕捉滤波误差的高阶统计量，使系统跟踪性能得到极大改善，降低了系统跟踪性能在重尾非高斯噪声下严重恶化的影响。

Description

一种基于最大相关熵扩展卡尔曼滤波的雷达目标跟踪方法

技术领域

本发明涉及雷达跟踪技术，属于雷达信号处理技术领域，具体涉及一种基于最大相关熵扩展卡尔曼滤波的雷达目标跟踪方法。该方法是基于最大相关熵(MCC)准则实现的，该方法既可以处理雷达目标跟踪中非线性问题，同时，由于采用了最大相关熵而非传统的最小均方误差准则(MMSE)作为优化准则，也可以解决系统噪声为非高斯噪声时目标跟踪的问题。

背景技术

在雷达跟踪系统中，由于目标位置的量测通常是距离、方位角、俯仰角等信息，而这些信息一般是在极坐标系得到的，这就使得雷达跟踪系统必然是非线性的。对于非线性滤波问题，目前最常用的滤波算法有扩展卡尔曼滤波(EKF)、不敏卡尔曼滤波(UKF)以及粒子滤波(PF)等。而在雷达量测系统中，往往又存在重尾非高斯噪声，而目前处理非高斯噪声的滤波算法则有学生t滤波、Huber’s M估计滤波等方法。然而，常用的非线性滤波方法在高斯假设下可以取得很好的滤波效果，但在非高斯情况下则出现滤波恶化等问题。而处理非高斯噪声的滤波方法则又会出现参数选择困难、计算复杂等问题。

对于非线性、非高斯滤波问题，本发明将最大相关熵扩展卡尔曼滤波(MCEKF)应用到雷达跟踪系统中，该滤波方法采用最大相关熵准则(MCC)而非传统的最小均方误差准则(MMSE)作为优化准则，相比于传统准则只能得到误差二阶项信息，该准则不仅可以获得误差二阶项信息，还可以捕捉滤波误差的高阶统计量，从而使得系统跟踪性能得到极大改善，降低了系统跟踪性能在重尾非高斯噪声下严重恶化的影响。

为了方便描述，下面首先介绍雷达跟踪系统中的滤波系统模型。

雷达跟踪系统中的滤波模型包括离散状态方程和离散量测方程。对于雷达来说，目标的状态信息通常是在直角坐标系下获得的，而量测信息则是在空间极坐标系下获得的，所以将雷达系统中的直角坐标系记为雷达直角坐标系OXYZ，极坐标系记为雷达极坐标系，其对应关系图1中所示。

雷达直角坐标系下离散时间系统的状态方程为，

X(k+1)＝f(X(k))+V(k) (1)

其中，X(k+1)为k+1时刻目标运动的状态向量，f(·)为目标运动的状态转移函数，V(k)为k时刻目标运动的过程噪声，且假定过程噪声服从零均值的高斯分布，其方差为，

E[V(k)V^T(j)]＝Q(k)δ_kj (2)

其中，

Q(k)为k时刻的过程噪声的协方差矩阵。

雷达极坐标系下离散时间系统的量测方程为，

其中，Z(k+1)＝[r(k+1) θ(k+1) φ(k+1)]^T，x(k+1)、y(k+1)和z(k+1)分别为k+1时刻目标在直角坐标系中x轴、y轴以及z轴的位置，r(k+1)、θ(k+1)和φ(k+1)分别为目标在k+1时刻的径向距离、方位角以及俯仰角的雷达量测信息，W(k+1)为量测噪声，其方差为，

E(W(k+1)W^T(j+1))＝R(k+1)δ_kj (5)

其中，R(k+1)为时刻k+1的量测噪声协方差矩阵，δ_kj的定义同式(3)。

其次介绍最大相关熵准则(MCC)。

最大相关熵是度量两个随机变量

之间相似性的一个重要的信息论量。假设它们的联合分布函数为F_XY(x,y)，它们之间的相关熵定义为，

V(X,Y)＝E[κ(X,Y)]＝∫κ(x,y)dF_XY(x,y) (6)

其中，E代表期望算子，κ(·，·)代表平移不变的Mercer核。在本发明中，所有的核函数是均是由高斯核函数给出的。其表达式为，

其中e＝x-y，σ＞0代表核带宽。

然而，在大多数雷达跟踪系统中，只有有限个数据可用，所以变量之间的联合概率密度F_XY通常情况下是未知的。在这种情况下，可以使用样本均值对相关熵进行估计。具体形式是，

这里e(i)＝x(i)-y(i)，

是从联合概率密度函数F_XY中抽取的N个样本。

对高斯核进行泰勒级数展开，可以得到，

可以看出，相关熵是X-Y所有偶阶矩的加权和，从中能够提取数据的高阶统计量。需要注意，当核带宽很大时，二阶矩将对相关熵起主要作用。

给定误差数据序列,基于MCC准则下的代价函数为：

发明内容

本发明旨在至少解决现有技术中存在的技术问题之一，提供一种基于最大相关熵扩展卡尔曼滤波的雷达目标跟踪方法。

本发明的一个方面提供一种基于最大相关熵扩展卡尔曼滤波的雷达目标跟踪方法，包括：

构建雷达状态方程和雷达量测方程；

使用预设的重尾非高斯噪声雷达量测噪声初始化雷达量测方程，使用预设的目标运动状态初始值初始化雷达状态方程；

根据k时刻雷达目标点迹的状态值和所述雷达状态方程得到k+1时刻雷达目标点迹的状态值；

根据所述k+1时刻的雷达目标点迹的状态值得到k+1时刻雷达目标点迹的量测值；

构建雷达状态先验更新方程和雷达量测先验更新方程；

根据所述雷达状态先验更新方程和k时刻雷达滤波值得到k+1时刻雷达目标点迹的状态预测值；

根据所述k+1时刻的雷达目标点迹的状态预测值和量测先验更新方程得到k+1时刻雷达目标点迹的量测预测值；

根据所述雷达状态方程、雷达状态先验更新方程和雷达量测方程构建雷达系统非线性模型，得到构建后的系统非线性模型；

根据所述k+1时刻雷达目标点迹的量测预测值和所述k+1时刻雷达目标点迹的量测值得到k+1时刻核函数对角阵奇异参数；

若所述k+1时刻核函数对角阵奇异参数大于预设的第一阈值，则将所述k+1时刻雷达目标点迹的状态预测值作为k+1时刻的雷达滤波值；若所述k+1时刻核函数对角阵奇异参数不大于预设的第一阈值，则使用所述构建后的雷达系统非线性模型进行定点迭代，将迭代结果作为k+1时刻的雷达滤波值；

若k+1＝N，其中，N为预设的雷达跟踪时长，则停止雷达跟踪；否则，循环执行上述步骤。

可选的，所述构建雷达状态方程、雷达状态先验更新方程、雷达量测方程，包括：

雷达状态方程：X(k+1)＝f(X(k))+V(k)

雷达量测方程：

Z(k+1)＝h(X(k+1))+W(k+1)

其中，

X(k+1)为k+1时刻雷达目标点迹的状态值；

X(k)为k时刻雷达目标点迹的状态值；

V(k)为k时刻雷达目标运动的过程噪声；

f(·)为雷达目标点迹的状态转移函数；

Z(k+1)为在k+1时刻雷达目标点迹的量测值；

x(k+1)、y(k+1)和z(k+1)分别为k+1时刻雷达目标在直角坐标系中x轴、y轴以及z轴的位置；

W(k+1)为k+1时刻雷达目标点迹的量测噪声；

h(·)为雷达目标量测的非线性函数。

可选的，根据所述雷达状态方程、雷达状态先验更新方程和雷达量测方程构建雷达系统非线性模型，得到构建后的雷达系统非线性模型，包括：

重组所述雷达状态方程、雷达状态先验更新方程和雷达量测方程，得到非线性模型组成参量；

根据所述非线性模型组成参量的协方差矩阵得到Cholesky分解参量；

根据所述Cholesky分解参量得到构建后的雷达系统非线性模型。

可选的，所述重组所述雷达状态方程、雷达状态先验更新方程和雷达量测方程，得到非线性模型组成参量，包括；

所述非线性模型组成参量为：

可选的，所述根据所述非线性模型组成参量的协方差矩阵得到Cholesky分解参量，包括：

其中，M(k+1)为Cholesky分解参量，R(k+1)为k时刻雷达目标点迹的量测噪声协方差矩阵

可选的，所述根据所述Cholesky分解参量得到构建后的雷达系统非线性模型：

D(k+1)＝B(X(k+1))+e(k+1)

其中，非线性模型的第一参量为：

非线性模型的第二参量为：

非线性模型的第三参量为：e(k+1)＝M^-1(k+1)α(k+1)

可选的，根据所述k+1时刻雷达目标点迹的状态预测值和所述k+1时刻雷达目标点迹的量测值得到k+1时刻核函数对角阵奇异参数，包括：

β(k+1)＝ψ^T(k+1)η^-1(k+1)ψ(k+1)

其中，

为k时刻雷达目标点迹的量测预测值；

Z(k+1)为k时刻雷达目标点迹的量测值；

为

的协方差；

h_X(k)为k时刻雷达目标量测的非线性函数h(·)的雅可比矩阵；

β(k+1)为k+1时刻核函数对角阵奇异参数。

可选的，所述使用所述构建后的雷达系统非线性模型进行定点迭代，包括：

根据k+1时刻雷达目标点迹的状态预测值得到k+1时刻迭代初始值；

根据k+1时刻雷达目标量测的非线性函数的雅可比矩阵得到k+1时刻迭代量测矩阵；

使用所述构建后的雷达系统非线性模型和预设的高斯核带宽，构建k+1时刻高斯核函数，根据所述k+1时刻高斯核函数得到对应的k+1时刻高斯核函数对角阵；

根据所述迭代初始值和迭代量测矩阵对所述k+1时刻高斯核函数对角阵进行至少一次迭代，得到每次迭代的结果；

判断所述每次迭代的结果是否符合预设的第二阈值，若是，则停止迭代，得到对应的迭代结果。

可选的，使用所述构建后的雷达系统非线性模型和预设的高斯核带宽，构建k+1时刻高斯核函数，包括：

其中，σ为高斯核带宽，

为k+1时刻t-1次迭代的迭代参数。

可选的，所述根据所述k+1时刻高斯核函数得到对应的k+1时刻高斯核函数对角阵，包括：

其中，x代表前n个高斯核函数，y代表后m个高斯核函数。

本发明实施例的一种基于最大相关熵扩展卡尔曼滤波的雷达目标跟踪方法中，对于非线性、非高斯滤波问题，本发明将最大相关熵扩展卡尔曼滤波(MCEKF)应用到雷达跟踪系统中，该滤波方法采用最大相关熵准则(MCC)而非传统的最小均方误差准则(MMSE)作为优化准则，相比于传统准则只能得到误差二阶项信息，该准则不仅可以获得误差二阶项信息，还可以捕捉滤波误差的高阶统计量，从而使得系统跟踪性能得到极大改善，降低了系统跟踪性能在重尾非高斯噪声下严重恶化的影响。在雷达跟踪中，使用最大相关熵扩展卡尔曼滤波的雷达目标跟踪方法既可以实现非线性跟踪系统的滤波，又可以有效降低非高斯重尾噪声使雷达跟踪性能严重恶化的影响，相比于传统的非线性滤波方法，其跟踪性能及其鲁棒性得到极大改善。

附图说明

图1为本发明一实施例的基于最大相关熵扩展卡尔曼滤波的雷达目标跟踪方法的流程图；

图2为本发明一实施例的雷达极坐标与雷达直角坐标对应关系示意图；

图3为本发明一实施例的基于最大相关熵扩展卡尔曼滤波的雷达目标跟踪滤波位置均方根误差仿真结果图；

图4为本发明一实施例的基于最大相关熵扩展卡尔曼滤波的雷达目标跟踪滤波速度均方根误差仿真结果图；

图5为本发明一实施例的基于最大相关熵扩展卡尔曼滤波的雷达目标跟踪滤波位置均方根误差仿真结果图；

图6为本发明一实施例的基于最大相关熵扩展卡尔曼滤波的雷达目标跟踪滤波速度均方根误差仿真结果图。

具体实施方式

为使本领域技术人员更好地理解本发明的技术方案，下面结合附图和具体实施方式对本发明作进一步详细描述。

如图1所示，本实施例提出一种基于最大相关熵扩展卡尔曼滤波的雷达目标跟踪方法S100。使用该方法能够解决雷达目标跟踪的非线性、非高斯问题，降低重尾非高斯噪声对系统跟踪性能的影响，从而提高雷达目标跟踪精度。

下面介绍一下该发明的具体步骤。

步骤S110：构建雷达状态方程和雷达量测方程。

具体地，在本步骤中，通过雷达状态方程和雷达量测方程组成滤波模型，使用滤波模型对雷达目标点迹进行滤波，并将滤波结果作为雷达跟踪的结果。

构建k时刻到k+1滤波模型的雷达状态方程以及雷达量测方程，雷达状态方程满足式(11),雷达量测方程满足式(12)。

X(k+1)＝f(X(k))+V(k) (11)

Z(k+1)＝h(X(k+1))+W(k+1) (12)

其中，k≥0且k为整数，代表雷达目标跟踪的时刻，当k＝0时，代表初始时刻。

其中，X(k)为k时刻雷达目标点迹的状态值，通常用向量表示，且X(k)＝[x(k),vx(k),y(k),vy(k),z(k),vz(k)]^T，x(k)、y(k)和z(k)分别为k时刻目标在直角坐标系中x轴、y轴以及z轴的位置，vx(k)、vy(k)和vz(k)分别为k时刻目标在直角坐标系中x轴、y轴以及z轴的速度。

f(·)为目标运动的状态转移函数。

X(k+1)为k+1时刻雷达目标点迹的状态值，通常用向量表示，且维度及含义与X(k)一致。

V(k)为k时刻目标运动的过程噪声。

Z(k+1)为k+1时刻雷达目标点迹的量测值，且Z(k+1)＝[r(k+1) θ(k+1) φ(k+1)]^T，r(k+1)、θ(k+1)和φ(k+1)分别为目标在k+1时刻雷达目标点迹的径向距离、方位角以及俯仰角的雷达量测信息，h(·)为目标量测的非线性函数，W(k+1)为k+1时刻目标的量测噪声。

步骤S120：使用预设的重尾非高斯噪声雷达量测噪声初始化雷达量测方程，使用预设的目标运动状态初始值初始化雷达状态方程。具体的，本步骤还包括：

步骤S121：将雷达量测方程中量测噪声设置为重尾非高斯噪声，本发明使用的重尾非高斯噪声为混合高斯噪声且服从混合高斯分布，如下式所示，

W(k+1)～(1-μ)(0,R(k+1))+μ(0,qR(k+1)) (13)

其中，μ代表误差模型污染的程度，q表示混合高斯噪声方差因子，具体实施过程中q可以选择不同值，μ通常设置为0.1。观察不同混合高斯噪声方差下的滤波效果。

步骤S122：设置目标运动状态初始值，包括目标运动的初始位置和初始速度，并使用该目标运动状态初始值初始化雷达目标点迹的状态值X(0)。设置雷达跟踪时间长度N，用于后期判断雷达是否可结束目标跟踪。并使用雷达量测值对滤波进行起始，即根据X(0)和(12)式得到Z(0)，再通过式(11)和式(12)得到X(1)和Z(1),使用两点差分法得到滤波初始值。

步骤S130：根据k时刻雷达目标点迹的状态值和所述雷达状态方程得到k+1时刻雷达目标点迹的状态值。根据所述k+1时刻的雷达目标点迹的状态值得到k+1时刻雷达目标点迹的量测值。

在完成对雷达状态方程、雷达量测方程的初始化后，即得到了k＝0时刻的雷达目标点迹的状态值和雷达目标点迹的量测值，即可据此通过雷达状态方程和雷达量测方程进一步得到后续时刻(k＝1、2、3……)的雷达目标点迹的状态值和雷达目标点迹的量测值。

S140：构建雷达状态先验更新方程和雷达量测先验更新方程：

雷达状态先验更新方程：

雷达量测先验更新方程：

其中，

为k时刻的雷达滤波值，在初始时刻，使用Z(0)和Z(1)根据极坐标系和直角坐标系对应关系得到滤波初始值

为k+1时刻的雷达目标点迹的状态预测值；

为k+1时刻的雷达目标点迹的量测预测值；

f(·)为雷达目标点迹的状态转移函数。

S150：根据所述雷达状态先验更新方程和k时刻雷达滤波值得到k+1时刻雷达目标点迹的状态预测值。根据所述k+1时刻的雷达目标点迹的状态预测值和量测先验更新方程得到k+1时刻雷达目标点迹的量测预测值。

在完成对雷达状态先验更新方程的初始化后，即得到了k＝2时刻的雷达目标点迹的状态预测值

即可据此通过雷达状态先验更新方程和雷达量测先验更新方程进一步得到后续时刻(k＝3，4,5……)的雷达目标点迹的状态预测值和雷达目标点迹的量测预测值。

S160：根据所述雷达状态方程、雷达状态先验更新方程和雷达量测方程构建雷达系统非线性模型，得到构建后的系统非线性模型。具体的，在本步骤中，还包括：

步骤S161：重组所述雷达状态方程、雷达状态先验更新方程和雷达量测方程，得到非线性模型组成参量。

结合雷达系统量测方程Z(k+1)＝h(X(k+1))+W(k+1)以及

重组非线性模型如下形式，

其中，非线性模型组成参量为

步骤S162：根据所述非线性模型组成参量的协方差矩阵得到Cholesky分解参量。

α(k+1)的协方差矩阵可以表示为，

其中，M(k+1)可以通过对上式进行Cholesky分解得到，为Cholesky分解参量。

步骤S163：根据所述Cholesky分解参量得到构建后的雷达系统非线性模型。

将式(14)两边分别左乘M^-1(k+1)，可以得到如下的非线性模型，

D(k+1)＝B(X(k+1))+e(k+1) (17)

其中，D(k+1)为构建非线性模型的第一参量，B(k+1)为构建非线性模型的第二参量，e(k+1)为构建非线性模型的第三参量，且分别满足

e(k+1)＝M^-1(k+1)α(k+1)。

完成上述步骤，下面进一步实现雷达目标跟踪过程，主要为实现目标状态由k时刻递推得到k+1时刻的目标跟踪实现步骤。

根据雷达状态先验更新方程，得到k+1时刻的雷达目标点迹的状态预测值，

然后得到其对应的一步预测协方差，

其中，f_X(k)为非线性函数f(·)的雅克比矩阵。

而后得到当前k+1时刻的雷达目标点迹的状态预测值所对应的雷达目标点迹的量测预测值，并求出雷达目标点迹的目标量测的非线性函数h(·)的雅可比矩阵h_X(k+1)。

由步骤S162中得到的M(k+1)分解获得M_P(k+1)与M_W(k+1)。

S170：根据所述k+1时刻雷达目标点迹的量测预测值和所述k+1时刻雷达目标点迹的量测值得到k+1时刻核函数对角阵奇异参数。

定点迭代过程开始前，对量测异常可能导致的核函数对角阵奇异问题进行处理。

β(k+1)＝ψ^T(k+1)η^-1(k+1)ψ(k+1) (23)

其中，β(k+1)为k+1时刻核函数对角阵奇异参数

S180：若所述k+1时刻核函数对角阵奇异参数大于预设的第一阈值，则将所述k+1时刻雷达目标点迹的状态预测值作为k+1时刻的雷达滤波值；若所述k+1时刻核函数对角阵奇异参数不大于预设的第一阈值，则使用所述构建后的雷达系统非线性模型进行定点迭代，将迭代结果作为k+1时刻的雷达滤波值；所述第一阈值可以为60、70、100。

如果|β(k+1)|＞θ且θ＞0(θ为预设的正阈值)，则执行预测步骤，将k+1时刻雷达目标点迹的状态预测值作为k+1时刻滤波结果，即

如果|β(k+1)|≤θ，则进行定点迭代过程。

具体地，在本步骤中，还包括：

步骤S181：根据k+1时刻雷达目标点迹的状态预测值得到k+1时刻迭代初始值。根据k+1时刻雷达目标量测的非线性函数的雅可比矩阵得到k+1时刻迭代量测矩阵。

开始迭代过程。其中t表示迭代索引，将k+1时刻雷达目标点迹的状态预测值、也就是状态一步预测值作为迭代初始值

t＝1开始且量测矩阵为

步骤S182：使用所述构建后的雷达系统非线性模型和预设的高斯核带宽，构建k+1时刻高斯核函数，根据所述k+1时刻高斯核函数得到对应的k+1时刻高斯核函数对角阵。

确定滤波过程中的核带宽σ，例如σ＝6、8、10，以及使定点迭代过程终止的ε(预设为很小的正阈值，ε＝10^-9)。

步骤S183：根据式(7)和式(17)构造k+1时刻高斯核函数G_σ，具体方式为，

根据核函数构造对角阵，具体为

其中，x代表前n个高斯核函数，y代表后m个高斯核函数。

步骤S184：根据所述迭代初始值和迭代量测矩阵对所述k+1时刻高斯核函数对角阵进行至少一次迭代，得到每次迭代的结果。根据下列式求解第t次迭代的迭代结果

及其相应的协方差

步骤S185：判断所述每次迭代的结果是否符合预设的第二阈值ε，若是，则停止迭代，得到对应的迭代结果。

比较当前迭代步的估计以及上一步的迭代估计，满足式(32)时结束迭代过程。如果不满足，重复步骤S183至S184，直至满足式(32)时结束迭代过程。

满足式(32)时的滤波结果及其相应的协方差则为k+1时刻时的雷达滤波值及对应的协方差值。

S190：若k+1＝N，其中，N为预设的雷达跟踪时长，则停止雷达跟踪；否则，循环执行上述步骤。

本发明的有益效果在于：在雷达跟踪中，使用最大相关熵扩展卡尔曼滤波的雷达目标跟踪方法既可以实现非线性跟踪系统的滤波，又可以有效降低非高斯重尾噪声使雷达跟踪性能严重恶化的影响，相比于传统的非线性滤波方法，其跟踪性能及其鲁棒性得到极大改善。

下面结合附图2至6的仿真结果，对本发明技术方案做进一步解释。

附图2是本发明雷达极坐标与雷达直角坐标对应关系示意图。参照附图2，并假设目标在三维平面内做匀速直线运动，那么，具体实施步骤描述如下。

步骤S00：构建k时刻到k+1滤波模型的状态方程以及量测方程，状态方程满足式(33),量测方程满足式(34)。

X(k+1)＝f(X(k))+V(k) (33)

Z(k+1)＝h(X(k+1))+W(k+1) (34)

其中，X(k)为k时刻目标运动的状态向量且X(k)＝[x(k),vx(k),y(k),vy(k),z(k),vz(k)]^T，x(k)、y(k)和z(k)分别为k时刻目标在直角坐标系中x轴、y轴以及z轴的位置，vx(k)、vy(k)和vz(k)分别为k时刻目标在直角坐标系中x轴、y轴以及z轴的速度，f(·)为目标运动的状态转移函数，X(k+1)为k+1时刻目标运动的状态向量且维度及含义与X(k)一致，V(k)为k时刻目标运动的过程噪声。Z(k+1)为k+1时刻雷达量测值，且Z(k+1)＝[r(k+1) θ(k+1) φ(k+1)]^T，r(k+1)、θ(k+1)和φ(k+1)分别为目标在k+1时刻的径向距离、方位角以及俯仰角的雷达量测信息，h(·)为目标量测的非线性函数，W(k+1)为k+1时刻目标的量测噪声。

步骤S01：将量测方程中量测噪声设置为重尾非高斯噪声，本发明使用的重尾非高斯噪声为混合高斯噪声且服从混合高斯分布，如下式所示，

W(k)～(1-μ)(0,R(k))+μ(0,qR(k)) (35)

步骤S02：设置目标运动状态，包括目标在雷达直角坐标系中的初始位置(100，500,300)Km和速度分别为(-200，200，200)m/s，及雷达跟踪时间长度为100s。并由雷达量测值对滤波进行起始值设置。

其中，状态转移矩阵为，

其中，T为采样间隔，T＝1s。并假定直角坐标系下各个坐标轴的过程噪声为高斯白噪声，即V(k)为，

其中，n_x、n_y和n_z分别为直角坐标系中x维、y维和z维的过程噪声，它们的均值均为0，方差为n_x＝n_y＝n_z＝1m/s²。

步骤S03：结合雷达系统量测方程Z(k+1)＝h(X(k+1))+W(k+1)以及

重组非线性模型如下形式，

其中，h(·)为目标量测的非线性函数，W(k+1)为k+1时刻目标的量测噪声，X(k+1)为k+1时刻目标的状态信息，通常为位置、速度、等信息，Z(k+1)为k+1时刻目标的量测信息，通常为径向距离、方位角、俯仰角等信息，且

α(k+1)的协方差矩阵可以表示为，

其中，M(k+1)可以通过对上式进行Cholesky分解得到。将式(38)两边分别左乘M^-1(k+1)，可以得到如下的非线性模型，

D(k+1)＝B(X(k+1))+e(k+1) (41)

e(k+1)＝M^-1(k+1)α(k+1)。

步骤S04：确定滤波过程中的核带宽，σ可以取2,4,6，以及使定点迭代过程终止的ε(预设为很小的正阈值，可以取ε＝10^-6、10^-7)。

步骤S05：目标状态由时刻k递推得到时刻k+1的目标跟踪实现步骤。

I、由于在直角坐标系中，目标运动通常是匀速的，所以状态方程通常是线性的，且状态先验更新值为，

然后得到其对应的一步预测协方差，

而后得到当前先验更新值所对应的量测预测值，并求出向量h的雅可比矩阵h_X(k+1)。

其中，三维极坐标下量测矩阵h_X(k+1)的具体计算方法为，

其中，

II、由步骤S03得到的M(k+1)分解获得M_P(k+1)与M_W(k+1)。

III、定点迭代过程开始前，对量测异常可能导致的核函数对角阵奇异问题进行处理。如果|β(k+1)|＞θ且θ＞0(为预设的正阈值，可以取θ＝100)，则执行预测步骤即

将先验更新值作为最终滤波结果。如果|β(k+1)|≤θ，则进行定点迭代过程。

IV：开始迭代过程。其中t表示迭代索引，将状态一步预测值作为迭代初始值

t＝1开始且初始量测矩阵为

(1)根据式(7)和式(17)构造高斯核函数，根据核函数构造对角阵。具体方式为，

(2)根据下列式求解第t次迭代的滤波值

及其相应的协方差

(3)比较当前迭代步的估计以及上一步的迭代估计，满足式(53)时结束迭代过程。如果不满足，重复步骤(1)和(2)直至满足式(53)时结束迭代过程，需要注意的是，在每次迭代步后，需要重新更新量测值的一步预测及其对应的量测矩阵，以及重新构造式(17)中的B(X(k+1))，以便在新的迭代步中使用。

(4)满足式(53)时的滤波结果及其相应的协方差则为k+1时刻时的滤波值及协方差值。

V：重复步骤S05直至滤波结束。

附图3到附图6为应用上述实施例中一种基于最大相关熵扩展卡尔曼滤波的雷达目标跟踪方法。的仿真结果图。在实施例中选择的滤波算法除上述最大相关熵扩展卡尔曼滤波算法(MCEKF)外，还选择了扩展卡尔曼滤波算法(EKF)进行仿真，并在量测噪声为高斯噪声和重尾非高斯噪声下进行对比。设置仿真参数如下表所示。

表1重尾非高斯噪声仿真参数

图3到图6表示了两种滤波算法的位置均方根误差及速度均方根误差。其中，滤波的位置均方根误差、速度均方根误差越小越好。

由仿真结果图3及图4可知，当量测噪声为高斯噪声时，传统扩展卡尔曼滤波(EKF)算法的滤波效果与本发明最大相关熵扩展卡尔曼滤波(MCEKF)算法的滤波效果近似相同，都能实现很好的跟踪，且传统EKF的滤波效果要略优于MCEKF的滤波效果。但是当雷达量测受到重尾非高斯噪声干扰时，由图5以及图6可以看到，EKF的跟踪性能出现明显恶化，目标跟踪性能严重下降。从滤波结果来看，MCEKF的滤波效果要优于EKF的滤波效果。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。