CN111948601B - 一种非高斯噪声条件下的单站纯角度目标定位与跟踪方法 - Google Patents

Abstract

本发明公开了一种非高斯噪声条件下的单站纯角度目标定位与跟踪方法，基于高斯混合容积卡尔曼滤波器进行非高斯噪声条件下的单站纯角度目标定位与跟踪算法设计，通过将非高斯噪声建模为两个不同方差的高斯噪声，并对状态量进行高斯分量建模，针对每一个高斯分量利用容积卡尔曼滤波器进行状态估计，并对滤波器输出进行高斯分量缩减合并，并迭代计算，最终得出对目标状态估计值。本发明的有益效果是，基于高斯混合容积卡尔曼滤波的单站纯角度目标定位与跟踪算法，将非高斯观测噪声建模为闪烁噪声，采用高斯分量分割与融合策略，利用容积卡尔曼对高斯分量进行状态估计，实现对目标的定位与跟踪，具有定位精度高，计算复杂度低等特点。

Description

一种非高斯噪声条件下的单站纯角度目标定位与跟踪方法

技术领域

本发明涉及纯角度目标定位与跟踪技术领域，尤其涉及的是，一种非高斯噪声条件下的单站纯角度目标定位与跟踪方法。

背景技术

无源定位系统本身不向外辐射电磁波，而是通过被动接收辐射信号来确定目标的位置和运动状态，具有隐蔽性强、作用距离远等优点。单站无源定位系统仅采用一个观测平台，避免了多站之间的协同工作，机动性强，不依赖站间通信，也降低了自身被探测与干扰的可能性。

现阶段对纯角度目标定位与跟踪的研究主要是基于高斯噪声假设，采用不同的滤波算法，提高对目标状态的估计精度，例如在李爱玲申请公开的发明专利《运动单站只测角位置信息无源定位方法》(申请公开号：CN 109959895 A)中，公开了一种运动单站只测角位置信息无源定位方法，可实现多目标轻型下的快速定位；在韩一娜等申请公开的发明专利《一种方位测量融合与多目标定位方法》(申请公开号：CN 110309599 A)中，公开一种可应用于多传感器网络中的多目标纯方位被动定位的空间概率模型，利用启发式搜索进行似然估计，上述方法都假设系统及观测噪声复合高斯分布，在复杂的实际环境中噪声特性并不确定，因此可能导致结果发散。在解决非高斯噪声方面，有将粒子滤波等随机采样方法应用于目标定位与跟踪问题，如鄢社锋等人申请公开的发明专利《一种基于简化容积粒子滤波的目标定位方法》(申请公开号：CN 109460539 A)中，公开了一种基于简化容积粒子滤波的目标定位方法，利用容积卡尔曼滤波器产生建议密度分布，指导粒子滤波器的采样，是目前精度最高的算法之一，但其计算量非常大，不适合用于实时目标定位与跟踪系统。

现有技术方案中存在如下不足：1、现有的单站纯角度目标定位与跟踪方法主要基于高斯假设的滤波算法，该类算法在高斯噪声条件下具有良好的跟踪性能，但是在非高斯条件下估计精度降低，甚至发散。2、可用于非高斯条件下单站纯角度目标定位与跟踪如粒子滤波算法，具有较好的估计精度，但是计算量巨大，不适用于实时目标定位与跟踪。

发明内容

本发明的目的在于提供一种基于高斯混合容积卡尔曼滤波的纯角度目标定位与跟踪方法，适用于非高斯噪声假设，可有效提高纯角度目标定位与跟踪精度，同时保持毫秒级的运算时间，可满足实时目标定位与跟踪的需要。所要解决的技术问题包括：

本发明的技术方案如下：一种非高斯噪声条件下的单站纯角度目标定位与跟踪方法，包括以下步骤：

步骤1：状态方程、观测方程建模；

单站纯角度目标定位与跟踪算法中状态方程和观测方程模型为公式1：

公式1：

公式1中，x_k表示目标状态，根据目标类型包含不同的状态量表示为代表目标的位置和速度；F为状态转移矩阵表示为/>T表示仿真步长；w_k-1表示过程噪声，建模为高斯分布；

测量方程中的测量值为对目标的角度测量可以表示为公式2：

公式2：

其中x_k,y_k表示目标位置，x_o,y_o表示观测站位置，本发明中表示无人机位置。将观测噪声v_k非高斯噪声建模为闪烁噪声，将其建模为两个不同方差的高斯噪声加和，表达式为公式3：

公式3：p(v)＝(1-ε)N(v；u₁,R₁)+εN(v；u₂,R₂)

公式3中，量测噪声的协方差为u₁＝u₂＝0，R₁＝R，R₂＝100R₁，R＝diag(10,10)闪烁系数ε＝0.1；

步骤2：状态初始化；

对系统状态协方差等进行初始化，P₀＝cov(x₀)，公式4：

公式4：

步骤3：时间更新；

对k-1时刻的状态误差协方差矩阵进行Cholesky分解为公式5：

公式5：

公式5中，r＝(i-1)I+j，I表示系统状态分解为高速分量的个数，J表示系统过程噪声分解成高斯分量个数；j取值从1到J；计算容积点，公式6：

公式6：

利用状态方程计算容积点预测估计值，公式7：

公式7：ξ_c,k|k-1(i)＝f(ξ_c,k-1|k-1(i))

k时刻系统状态的预测估计值，公式8：

公式8：

k时刻系统状态误差协方差矩阵预测估计值，公式9：

公式9：

k时刻高斯分量权值计算如下，公式10：

公式10：

公式10中m表示容积点的个数，表示时间更新后各高斯分量的权值；/>表示过程噪声的第j个高斯分量的权值；/>表示k-1时刻系统状态第i个高斯分量的权值；/>表示经过状态方程传递的高斯分量的权值；

步骤4：量测更新；

对公式9的状态误差协方差矩阵预测值进行Cholesky分解，公式11：

公式11：

计算容积点，公式12：

公式12：

利用量测方程对容积点进行传递，公式13：

公式13：ε_c,k|k-1(r)＝h(ξ_c,k|k-1(r))

计算k时刻的量测预测值，公式14：

公式14：

计算k时刻量测预测自相关协方差矩阵，公式15：

公式15：

计算k时刻互相关协方差矩阵，公式16：

公式16：

计算k时刻卡尔曼增益，公式17：

公式17：

计算k时刻的状态更新值，公式18：

公式18：

计算k时刻状态误差协方差矩阵，公式19：

公式19：

公式19中，n＝(r-1)L+l，L表示系统观测噪声分解成的高斯分量个数，l

取值从1到L；

步骤5：计算滤波器输出结果；

k时刻的系统状态和协方差的估计值分别为，公式20及公式21：

公式20：

公式21：

公式21中，ω_k(n)表示第n个高斯分量的权值，在每次迭代时按照公式22计算：

公式22：

公式22中，表示量测噪声第l个高斯分量的权值，p(z_k|x_k,n)表示第n个高斯分量的量测似然分布，其计算方式为公式23：

公式23：步骤6：高斯分量缩减合并；

缩减后的高斯分量权值为公式24：

公式24：

其均值和协方差分别为公式25及公式26：

公式25：

公式26：

上述公式25及公式26中，表示被缩减的高斯分量归一化权值；

步骤7：判断K+1是否大于N，其中K为当前时刻，N为总迭代时间，如果是则迭代计算结束，输出目标估计结果；如果否，则迭代K＝K+1后，返回步骤3继续进行时间更新步骤。

本发明的有益效果是：1、基于高斯混合容积卡尔曼滤波的单站纯角度目标定位与跟踪算法，将非高斯观测噪声建模为闪烁噪声，采用高斯分量分割与融合策略，利用容积卡尔曼对高斯分量进行状态估计，实现对目标的定位与跟踪，具有定位精度高，计算复杂度低等特点。2、将状态量分解为多个高斯分量，分别利用容积卡尔曼滤波进行目标状态估计，提高非高斯条件下的估计精度；3、将非高斯观测噪声建模为闪烁噪声，利用多个不同方差的高斯分布表示，便于仿真实现；4、采用高斯分量分割与融合策略，降低算法复杂度，提高实时性。

附图说明

图1为本发明方法流程图。

图2为本发明仿真状态估计曲线图。

图3为本发明仿真位置均方根误差曲线图。

具体实施方式

为了便于理解本发明，下面结合附图和具体实施例，本说明书所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本说明书中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是用于限制本发明。本说明书所使用的术语“和/或”包括一个或多个相关的所列项目的任意的和所有的组合。

本发明的一个实施例是，基于高斯混合容积卡尔曼滤波器进行非高斯噪声条件下的单站纯角度目标定位与跟踪算法设计，通过将非高斯噪声建模为两个不同方差的高斯噪声，并对状态量进行高斯分量建模，针对每一个高斯分量利用容积卡尔曼滤波器进行状态估计，并对滤波器输出进行高斯分量缩减合并，并迭代计算，最终得出对目标状态估计值。本发明所述方法流程图如图1所示。

本发明具体实现过程如下：

步骤1：状态方程、观测方程建模；

单站纯角度目标定位与跟踪算法中状态方程和观测方程模型为公式1：

公式1：

公式1中，x_k表示目标状态，根据目标类型包含不同的状态量，在本发明中表示为代表目标的位置和速度；F为状态转移矩阵，本发明中表示为T表示仿真步长；w_k-1表示过程噪声，建模为高斯分布；测量方程中的测量值为对目标的角度测量可以表示为公式2：

公式2：

其中x_k,y_k表示目标位置，x_o,y_o表示观测站位置，本发明中表示无人机位置。

在本发明中将观测噪声v_k非高斯噪声建模为闪烁噪声，闪烁噪声一般可建模为高斯噪声和拉普拉斯噪声的混合形式，或多个不同方差的高斯分布混合。为了便于仿真实现，本发明中将其建模为两个不同方差的高斯噪声加和，表达式为公式3：

公式3：p(v)＝(1-ε)N(v；u₁,R₁)+εN(v；u₂,R₂)

公式3中，量测噪声的协方差为u₁＝u₂＝0，R₁＝R，R₂＝100R₁，R＝diag(10,10)闪烁系数ε＝0.1。

步骤2：状态初始化；

对系统状态协方差等进行初始化，P₀＝cov(x₀)，公式4：

公式4：

步骤3：时间更新；

对k-1时刻的状态误差协方差矩阵进行Cholesky分解为公式5：

公式5：

公式5中，r＝(i-1)I+j，I表示系统状态分解为高速分量的个数，J表示系统过程噪声分解成高斯分量个数；j取值从1到J。

计算容积点，公式6：

公式6：

利用状态方程计算容积点预测估计值，公式7：

公式7：ξ_c,k|k-1(i)＝f(ξ_c,k-1|k-1(i))

k时刻系统状态的预测估计值，公式8：

公式8：

k时刻系统状态误差协方差矩阵预测估计值，公式9：

公式9：

k时刻高斯分量权值计算如下，公式10：

公式10：

公式10中m表示容积点的个数，表示时间更新后各高斯分量的权值；/>表示过程噪声的第j个高斯分量的权值；/>表示k-1时刻系统状态第i个高斯分量的权值；/>表示经过状态方程传递的高斯分量的权值。

步骤4：量测更新；

对公式9的状态误差协方差矩阵预测值进行Cholesky分解，公式11：

公式11：

计算容积点，公式12：

公式12：

利用量测方程对容积点进行传递，公式13：

公式13：ε_c,k|k-1(r)＝h(ξ_c,k|k-1(r))

计算k时刻的量测预测值，公式14：

公式14：

计算k时刻量测预测自相关协方差矩阵，公式15：

公式15：

计算k时刻互相关协方差矩阵，公式16：

公式16：

计算k时刻卡尔曼增益，公式17：

公式17：

计算k时刻的状态更新值，公式18：

公式18：

计算k时刻状态误差协方差矩阵，公式19：

公式19：

公式19中，n＝(r-1)L+l，L表示系统观测噪声分解成的高斯分量个数，l取值从1到L。

步骤5：计算滤波器输出结果；

k时刻的系统状态和协方差的估计值分别为，公式20及公式21：

公式20：

公式21：

公式21中，ω_k(n)表示第n个高斯分量的权值，在每次迭代时按照公式22计算：

公式22：

公式22中，表示量测噪声第l个高斯分量的权值，p(z_k|x_k,n)表示第n个高斯分量的量测似然分布，其计算方式为公式23：

公式23：

步骤6：高斯分量缩减合并

缩减后的高斯分量权值为公式24：

公式24：

其均值和协方差分别为公式25及公式26：

公式25：

公式26：

公式25及公式26中，表示被缩减的高斯分量归一化权值。

步骤7：判断K+1是否大于N，其中K为当前时刻，N为总迭代时间，如果是则迭代计算结束，输出目标估计结果；如果否，则迭代K＝K+1后，返回步骤3继续进行时间更新步骤。

本发明上述方法的仿真验证，将本发明中仿真条件设置如下：运动单站假设为无人机，起始点为(0,0)，以15m/s的速度做螺旋轨迹运动，目标初始点为(3000,-500)，以9m/s速度沿45°航向做匀速直线运动，无人机对目标的测角误差为3°，仿真曲线如下图2-3所示，从图2-3中可以看出，本发明设计的单站纯角度目标定位与跟踪算法能够快速准确实现对运动目标的定位与跟踪，定位精度可收敛到40m以内，与基于粒子滤波的算法相比，定位精度差距不大，但收敛速度更快，同时本发明设计的算法属于毫秒级计算，能够满足实时单站纯角度目标定位与跟踪的应用。

本发明的有益效果是：1、基于高斯混合容积卡尔曼滤波的单站纯角度目标定位与跟踪算法，将非高斯观测噪声建模为闪烁噪声，采用高斯分量分割与融合策略，利用容积卡尔曼对高斯分量进行状态估计，实现对目标的定位与跟踪，具有定位精度高，计算复杂度低等特点。2、将状态量分解为多个高斯分量，分别利用容积卡尔曼滤波进行目标状态估计，提高非高斯条件下的估计精度；3、将非高斯观测噪声建模为闪烁噪声，利用多个不同方差的高斯分布表示，便于仿真实现；4、采用高斯分量分割与融合策略，降低算法复杂度，提高实时性。

需要说明的是，上述各技术特征继续相互组合，形成未在上面列举的各种实施例，均视为本发明说明书记载的范围；并且，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (1)

1.一种非高斯噪声条件下的单站纯角度目标定位与跟踪方法，其特征在于，包括以下步骤：

步骤1：状态方程、观测方程建模；

单站纯角度目标定位与跟踪算法中状态方程和观测方程模型为公式1：

公式1：

公式1中，x_k表示目标状态，根据目标类型包含不同的状态量表示为代表目标的位置和速度；F为状态转移矩阵表示为/>T表示仿真步长；w_k-1表示过程噪声，建模为高斯分布；

测量方程中的测量值为对目标的角度测量可以表示为公式2：

公式2：

其中x_k,y_k表示目标位置，x_o,y_o表示观测站位置；将观测噪声v_k非高斯噪声建模为闪烁噪声，将其建模为两个不同方差的高斯噪声加和，表达式为公式3：

公式3：p(v)＝(1-ε)N(v；u₁,R₁)+εN(v；u₂,R₂)

公式3中，量测噪声的协方差为u₁＝u₂＝0，R₁＝R，R₂＝100R₁，R＝diag(10,10)闪烁系数ε＝0.1；

步骤2：状态初始化；

对系统状态协方差进行初始化，P₀＝cov(x₀)，公式4：

公式4：

步骤3：时间更新；

对k-1时刻的状态误差协方差矩阵进行Cholesky分解为公式5：

公式5：

公式5中，r＝(i-1)I+j，I表示系统状态分解为高速分量的个数，J表示系统过程噪声分解成高斯分量个数；j取值从1到J；计算容积点，公式6：

公式6：

利用状态方程计算容积点预测估计值，公式7：

公式7：ξ_c,k|k-1(i)＝f(ξ_c,k-1|k-1(i))

k时刻系统状态的预测估计值，公式8：

公式8：

k时刻系统状态误差协方差矩阵预测估计值，公式9：

公式9：

k时刻高斯分量权值计算如下，公式10：

公式10：

公式10中m表示容积点的个数，表示时间更新后各高斯分量的权值；/>表示过程噪声的第j个高斯分量的权值；/>表示k-1时刻系统状态第i个高斯分量的权值；表示经过状态方程传递的高斯分量的权值；

步骤4：量测更新；

对公式9的状态误差协方差矩阵预测值进行Cholesky分解，公式11：

公式11：

计算容积点，公式12：

公式12：

利用量测方程对容积点进行传递，公式13：

公式13：ε_c,k|k-1(r)＝h(ξ_c,k|k-1(r))

计算k时刻的量测预测值，公式14：

公式14：

计算k时刻量测预测自相关协方差矩阵，公式15：

公式15：

计算k时刻互相关协方差矩阵，公式16：

公式16：

计算k时刻卡尔曼增益，公式17：

公式17：

计算k时刻的状态更新值，公式18：

公式18：

计算k时刻状态误差协方差矩阵，公式19：

公式19：

公式19中，n＝(r-1)L+l，L表示系统观测噪声分解成的高斯分量个数，l取值从1到L；

步骤5：计算滤波器输出结果；

k时刻的系统状态和协方差的估计值分别为，公式20及公式21：

公式20：

公式21：

公式21中，ω_k(n)表示第n个高斯分量的权值，在每次迭代时按照公式22计算：

公式22：

公式22中，表示量测噪声第l个高斯分量的权值，p(z_k|x_k,n)表示第n个高斯分量的量测似然分布，其计算方式为公式23：

公式23：

步骤6：高斯分量缩减合并；

缩减后的高斯分量权值为公式24：

公式24：

其均值和协方差分别为公式25及公式26：

公式25：

公式26：

上述公式25及公式26中，表示被缩减的高斯分量归一化权值；

步骤7：判断K+1是否大于N，其中K为当前时刻，N为总迭代时间，如果是则迭代计算结束，输出目标估计结果；如果否，则迭代K＝K+1后，返回步骤3继续进行时间更新步骤。