CN115932907B - 基于卡尔曼滤波器和粒子群优化的单目标直接跟踪方法 - Google Patents

Abstract

本发明公开了基于卡尔曼滤波器和粒子群优化的单目标直接跟踪方法，首先得到接收信号；其次，跟踪目标每个时刻的信号子空间，并且利用子空间理论，获得每个时刻的噪声子空间；接着，联合所有基站接收到的信号进行数据融合，依据MUSIC谱构建和当前时刻目标位置有关的代价函数；然后，利用粒子群算法求解代价函数的最小值，得到每个时刻目标的位置；最后通过卡尔曼滤波器，对每个时刻目标的位置进行处理，提高了定位精度。本发明能够极大的降低直接定位算法的复杂度，在实时跟踪的场景下，可以满足快速跟踪的要求。相对于传统的两步定位算法、子空间数据融合算法，在高信噪比时有着相当的定位精度；在低信噪比时，有着更高的定位精度。

Description

基于卡尔曼滤波器和粒子群优化的单目标直接跟踪方法

技术领域

本发明属于无线定位技术领域，尤其涉及基于卡尔曼滤波器和粒子群优化的单目标直接跟踪方法。

背景技术

传统的无源定位技术，大都采用两步定位算法，接收机收集到信号后，需要进行中间参数的估计，之后利用几何关系确定目标位置。在此过程中，信息存在损失，因此两步定位算法在需要高精度定位的场景下，难以完成任务。直接定位算法解决了这个问题，采用直接构建和目标位置有关的代价函数，直接求解代价函数获得目标位置的方法，没有信息损失，因而定位精度高，特别是在低信噪比时表现远远优于两步定位。

但是所构建的代价函数如何求解，是一个重要的问题。现有的方法，最常用的是网格搜索法，将所有的可能目标位置划分成网格，逐一计算代价函数的值，从而选取最优值。因此，该方法的计算复杂度很高，特别是要求高定位精度情况时，需要非常精细的网格，带来巨大的运算量。

实际情况下，目标位置是移动的情况更多，因此我们需要每个时刻不停地对目标位置进行跟踪，因此若仍然采用网格搜索法，其巨大的运算量将会导致难以进行实时地跟踪目标位置。为了解决上述问题，本发明提出了一种基于卡尔曼滤波器和粒子群优化的单目标直接跟踪方法。

发明内容

发明目的：为了解决现有技术存在的问题，本发明提供了一种基于卡尔曼滤波器和粒子群优化的单目标直接跟踪方法，在保证估计性能的同时显著降低了计算复杂度，易于实时处理。

技术方案：本发明所述的一种基于卡尔曼滤波器和粒子群优化的单目标直接跟踪方法，包括以下步骤：

基于卡尔曼滤波器和粒子群优化的单目标直接跟踪方法，包括以下步骤：

1)构造多阵列单目标直接跟踪模型，获取时刻t的接收信号r_l(t)；

2)跟踪目标当前时刻的信号子空间W_l，并且利用空间谱理论，获得当前时刻每个阵列的噪声子空间；

3)联合所有观测基站进行数据融合，构建直接定位的代价函数P_MUSIC；

4)利用粒子群算法解算代价函数P_MUSIC，得到目标位置；

5)利用卡尔曼滤波处理得到的目标位置，得到更准确的结果。

优选的，步骤1)中，所述多阵列单目标直接跟踪模型包括若干个观测基站和一个目标即辐射源，每个所述观测基站均配备有阵列，将得到的接收信号r_l(t)定义为：

r_l(t)＝a_l(p)s_l(t)+n_l(t)

其中，为导向矢量，d表示阵列中阵元的间距，λ表示波长，M表示阵元数量；u_l＝[x_l,y_l]^T表示第l个观测基站的位置，辐射源的位置为p＝[x,y]^T，u_l(1)和p(1)表示u_l(1)和p(1)的第一个元素，θ_l表示第l个观测基站接收到的辐射源的波达角度；s_l(t)表示第l个观测基站接收到辐射源的回波信号；n_l(t)表示高斯白噪声矢量。

优选的，步骤2)中，获取噪声子空间的实现过程为：

2.1)适当选择自相关矩阵逆P(j)的初值P(0)和信号子空间W_l的初值W_l(0)；

2.2)计算中间参数y(j)＝W_l ^H(j-1)x(j)和中间参数h(j)＝P(j-1)y(j)；

2.3)计算增益向量g(j)＝h(j)/[β+y^H(j)h(j)]；

2.4)更新自相关矩阵的逆矩阵

2.5)计算误差向量e(j)＝x(j)-W_l(j-1)y(j)；

2.6)更新信号子空间W_l(j)＝W_l(j-1)+e(j)g^H(j)；

2.7)重复2.2)-2.6)，利用每个时刻的快拍来不断更新信号子空间得到信号子空间W_l；

2.8)根据信号子空间W_l获取噪声子空间E_l；

其中，j为快拍数的循环值，x(j)为每个快拍的信号，β为遗忘因子，Tri{O}＝triu(O)+triu(O)^H，O表示任意矩阵，triu()表示取矩阵的上三角元素。

优选的，步骤3)中：

根据空间谱理论，信号子空间和噪声子空间具有正交性，当阵列的导向矢量a_l(p)由辐射源的位置p构成时，导向矢量到噪声子空间的投影为零，即满足：

其中，是噪声子空间；

用模来表示：

接着改写为：

联合所有阵列，构建代价函数：

其中，L为阵列的个数。

优选的，步骤4)中利用粒子群算法求解直接获得目标该时刻位置的测量值。

优选的，步骤5)中：首先建立卡尔曼滤波器的状态方程和测量方程：

y(t+1)＝Fy(t)+Bh+c(t)

z(t)＝Hy(t)+e(t)

其中，y(t)＝[x₁(t),x₂(t),v₁(t),v₂(t)]表示t时刻的状态向量，x₁(t)、x₂(t)、v₁(t)、v₂(t)分别表示t时刻目标辐射源的X轴方向坐标、Y轴方向坐标、沿X轴方向的速度、沿Y轴方向的速度；h表示控制量，且h＝[0,0]^T；F表示状态转移矩阵、B表示控制量转移矩阵、H表示测量矩阵；c(t)为过程噪声，服从零均值高斯分布，其协方差矩阵是Q(t)；e(t)是测量噪声，也服从零均值高斯分布，其方差为R(t)，卡尔曼滤波器的初值设置为：

其中T_s是跟踪的时间间隔，假设已知t-1时刻的状态向量和状态协方差矩阵/>通过卡尔曼方程得到t时刻的状态向量和状态协方差矩阵的先验估计：

其中，中的前两个元素为t-1时刻目标的X轴和Y轴坐标；利用目标的位置坐标，获得目标位置的测量值z(t)；

修正先验估计，获得位置的后验估计和协方差矩阵的后验估计：

I表示单位矩阵，K(t)表示卡尔曼增益：

每个时刻都进行卡尔曼滤波，获得更准确的估计结果。

有益效果：与现有技术相比，本发明的有益效果：本发明所提算法的估计精度优于传统两步定位算法、使用网格搜索的子空间数据融合算法、使用粒子群算法的子空间数据融合算法；与传统两步定位算法、使用网格搜索的子空间数据融合算法相比，本发明在低信噪比时，有着明显更好的定位精度；同时在高信噪比时可以在保证估计性能的情况下显著降低计算复杂度。

附图说明

图1为本发明的流程图；

图2为单目标直接跟踪场景图；

图3为本发明与传统定位方法在不同观测站数目下的计算复杂度示意图；

图4为本发明与传统定位方法在不同信噪比下的均方根误差性能示意图；

图5为本发明与传统定位方法在不同快拍数下的均方根误差性能示意图。

图6为本发明与传统定位方法在不同阵元数下的均方根误差性能示意图。

具体实施方式

下面结合附图对本发明做进一步详细说明。

本发明提供一种基于卡尔曼滤波器和粒子群优化的单目标直接跟踪方法，如图1所示，具体包括以下步骤：

步骤1：构造如图2所示的多阵列单目标直接跟踪模型；得到每个时刻的接收信号信息r_l(t)。

假设对信源进行了T秒的跟踪，跟踪间隔为T_s秒，则满足T＝nT_s，其中n为正整数。在短时间[(n-1)T_s,nT_s],n＝1,2,…N_s内，我们可以近似认为信源的位置不变，此时对数据进行J次采样可以用来估计此时信源的位置。

每个采样时刻，第l个观测基站的接收信号为其中为导向矢量，d表示阵元间距。u_l＝[x_l,y_l]^T表示第l个观测站的位置，辐射源的位置为p＝[x,y]^T，θ_l表示第l个观测站接收到的辐射源的波达角度。s_l(t)表示第l个观测基站接收到辐射源的回波信号；n_l(t)表示高斯白噪声矢量。

步骤2：利用PAST算法跟踪目标当前时刻的信号子空间W_l，并且利用空间谱理论，获得当前时刻每个阵列的噪声子空间。

根据PAST算法原理，我们可以在每个时刻利用收集到的J次快拍不断迭代更新信号子空间。下面给出PAST算法在单目标场景下的流程：

(1)适当选择自相关矩阵逆的初值P(0)和信号子空间的初值的W_l(0)；

(2)计算中间参数y(j)＝W_l ^H(j-1)x(j)和h(j)＝P(j-1)y(j)；

(3)计算增益向量g(j)＝h(j)/[β+y^H(j)h(j)]；

(4)更新自相关矩阵的逆矩阵

(5)计算误差向量e(j)＝x(j)-W_l(j-1)y(j)；

(6)更新信号子空间W_l(j)＝W_l(j-1)+e(j)g^H(j)；

(7)重复(2)-(6)，利用每个时刻的快拍来不断更新信号子空间。

其中，Tri{P(j-1)-g(j)h^H(j)}＝triu(P(j-1)-g(j)h^H(j))+triu(P(j-1)-g(j)h^H(j))^H，triu()表示取矩阵的上三角元素。简单的初值设置可以是将P(0)设置为单位矩阵，此处即为数字1，而W_l(0)设置为M×M单位矩阵的前导单位向量。算法最后获得的W_l即为估计到的信号子空间，对于不同阵列来说，获得了不同的W_l。然后，我们利用空间谱理论，可以计算出信号的噪声子空间，

步骤3：联合所有基站进行数据融合，构建直接定位的代价函数P_MUSIC。

根据空间谱理论，信号子空间和噪声子空间具有正交性，当阵列的导向矢量a_l(p)由真实的辐射源位置p构成时，导向矢量到噪声子空间的投影为零。即满足

将上式用模来表示

又可以改写为

因此我们联合L个阵列，构建以下代价函数

步骤4：利用粒子群算法解算代价函数P_MUSIC，得到目标位置。

本发明讨论的目标跟踪问题为一个D＝2维的优化问题。粒子群算法作为一种智能优化算法，可以很好的解决该问题。粒子群算法相较于传统的网格搜索法，有着不错的精度，同时算法的复杂性很低，因此在实时跟踪场景下，运用更加广泛。我们利用粒子群算法可以由上式直接获得目标该时刻位置的测量值。

步骤5：利用卡尔曼滤波处理目标位置，得到更准确的结果。

卡尔曼滤波器可以根据预测模型，对测量值给出修正，从而减少噪声的影响，目标跟踪的性能更好。首先建立卡尔曼滤波器的状态方程和测量方程

y(t+1)＝Fy(t)+Bu+c(t)

z(t)＝Hy(t)+e(t)

其中，y＝[x₁(t),x₂(t),v₁(t),v₂(t)]表示状态向量，x₁(t)、x₂(t)、v₁(t)、v₂(t)分别表示目标信源的X轴方向坐标、Y轴方向坐标、沿X轴方向的速度、沿Y轴方向的速度。u表示控制量，在本发明中为两个方向的加速度，u＝[0,0]^T。F表示状态转移矩阵、B表示控制量转移矩阵、H表示测量矩阵。c(t)为过程噪声，服从零均值高斯分布，协方差矩阵是Q(t)。e(t)是测量噪声，也服从零均值高斯分布，方差为R(t)。相应量的初值设置为

我们假设已知t-1时刻的状态向量和协方差矩阵/>通过卡尔曼方程可以得到t时刻的状态向量和协方差矩阵的先验估计

其中，中的前两个元素就是t-1时刻目标的X轴和Y轴坐标。利用目标的位置坐标，我们可以利用上文步骤提及的PAST算法来获得信号子空间、噪声子空间，然后利用MUSIC谱构建代价函数，使用粒子群算法解算代价函数，获得位置的测量值z(t)。

接下来，通过下式来修正先验估计，获得位置的后验估计和协方差矩阵的后验估计

I表示单位矩阵，此处为4×4单位矩阵。K(t)表示卡尔曼增益，由下式计算

至此，我们介绍完了卡尔曼滤波的五个方程，每个时刻都进行卡尔曼滤波，我们就可以获得更准确的估计结果。

本发明方法中，粒子群算法解算代价函数时，每次更新粒子位置和速度只要进行5次实数加法运算和5次实数乘法运算，因此相对于计算适应度函数可以忽略。计算适应度函数时，算法仍需依据MUSIC谱来建立，所以带来了一些运算量，复杂度为O[γL(αM(α+M))]，其中γ表示迭代次数，L表示基站个数，α表示粒子数，M表示阵元数。然而在传统的网格搜索方法中，所需要的复杂度为O(gL(M(M+1))，g表示进行二维搜索的次数，通常是一个很大的数。在单目标跟踪问题中，高信噪比时，网格通常取的非常小，这样就g的值将会非常大。图3为本发明所述方法与传统网格搜索法直接定位的计算复杂度随观测基站数目变化的示意图。仿真条件为：每个基站配备一个阵元数为10的均匀线阵，全局搜索范围为(0,1000)米，搜索布长为0.5米，因为算法都使用了PAST来获得子空间，因此我们只比较解算部分的复杂度。粒子群算法中，粒子数为100，迭代次数也为100次。从图3中可以看出来，使用粒子群算法，避免了网格搜索的高额计算量，计算复杂度显著降低。

本发明性能估计标准为均方根误差(root mean square error,RMSE)定义为：

其中，MON为蒙特卡罗实验次数，T表示跟踪时间，表示第mon次实验，t时刻目标位置的估计值，p_mon,t表示第mon次实验，t时刻目标位置的真实值。

图4为本发明所述方法与传统的两部定位方法、子空间数据融合算法、未使用卡尔曼滤波的粒子群算法的性能随信噪比变化的曲线。仿真条件为：目标位于(200，200)米，横纵坐标都以每秒10米的速度增加，跟踪总时间为40秒，每隔1秒进行一次跟踪，收集40个快拍，信噪比从-10dB到10dB，间隔为5dB。使用5个观测基站，位置为(0,0)米、(0,250)米、(0,500)米、(0,750)米、(0,1000)米，每个基站配备一个阵元数为10的均匀线阵，全局搜索范围为(0,1000)米，搜索布长为0.5米。粒子群算法中，粒子数为100，迭代次数也为100次。从图4中可以看出来，使用基于卡尔曼滤波和粒子群优化的算法，在低信噪比时，可以减少干扰，实现了更高的定位精度，在高信噪比时，也可以与网格搜索法保持相当的精度。

图5为发明所述方法与传统的两部定位方法、子空间数据融合算法、未使用卡尔曼滤波的粒子群算法的性能随快拍数变化的曲线。仿真条件为：目标位于(200，200)米，横纵坐标都以每秒10米的速度增加，跟踪总时间为40秒，每隔1秒进行一次跟踪，快拍数从20到60，间隔10，信噪比为-5dB。使用5个观测基站，位置为(0,0)米、(0,250)米、(0,500)米、(0,750)米、(0,1000)米，每个基站配备一个阵元数为10的均匀线阵，全局搜索范围为(0,1000)米，搜索布长为0.5米。粒子群算法中，粒子数为100，迭代次数也为100次。从图5中可以看出来，使用所提方法，快拍数增加时，PAST算法跟踪信号子空间更加准确，算法精度更高。

图6为发明所述方法与传统的两部定位方法、子空间数据融合算法、未使用卡尔曼滤波的粒子群算法的性能随阵元个数变化的曲线。仿真条件为：目标位于(200，200)米，横纵坐标都以每秒10米的速度增加，跟踪总时间为40秒，每隔1秒进行一次跟踪，快拍数为40，间隔10，信噪比为-5dB。使用5个观测基站，位置为(0,0)米、(0,250)米、(0,500)米、(0,750)米、(0,1000)米，每个基站配备一个均匀线阵，阵元数从6变化到10，间隔为1，全局搜索范围为(0,1000)米，搜索步长为0.5米。粒子群算法中，粒子数为100，迭代次数也为100次。从图6中可以看出来，使用所提方法，阵元数增加时，算法获得更高的定位精度。

综上所述，从仿真效果图的分析可知，本发明提出的一种基于卡尔曼滤波器和粒子群优化的单目标直接跟踪方法能够极大的降低计算复杂度，可以满足目标跟踪的实时需求，同时该算法在环境恶劣情况下，信噪比较低时，可以明显的提高定位精度，有着更好的性能。

本发明未详述之处，均为本领域技术人员的公知技术。

上面结合附图对本发明的实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。

Claims (4)

1.基于卡尔曼滤波器和粒子群优化的单目标直接跟踪方法，其特征在于，包括以下步骤：

1)构造多阵列单目标直接跟踪模型，获取时刻t的接收信号r_l(t)；

2)跟踪目标当前时刻的信号子空间W_l，并且利用空间谱理论，获得当前时刻每个阵列的噪声子空间；

3)联合所有观测基站进行数据融合，构建直接定位的代价函数P_MUSIC；

4)利用粒子群算法解算代价函数P_MUSIC，得到目标位置；

5)利用卡尔曼滤波处理得到的目标位置，得到更准确的结果；

步骤2)中，获取噪声子空间的实现过程为：

2.1)适当选择自相关矩阵逆P(j)的初值P(0)和信号子空间W_l的初值W_l(0)；

2.2)计算中间参数y(j)＝W_l ^H(j-1)x(j)和中间参数h(j)＝P(j-1)y(j)；

2.3)计算增益向量g(j)＝h(j)/[β+y^H(j)h(j)]；

2.4)更新自相关矩阵的逆矩阵

2.5)计算误差向量e(j)＝x(j)-W_l(j-1)y(j)；

2.6)更新信号子空间W_l(j)＝W_l(j-1)+e(j)g^H(j)；

2.7)重复2.2)-2.6)，利用每个时刻的快拍来不断更新信号子空间得到信号子空间W_l；

2.8)根据信号子空间W_l获取噪声子空间E_l；

其中，j为快拍数的循环值，x(j)为每个快拍的信号，β为遗忘因子，Tri{O}＝triu(O)+triu(O)^H，O表示任意矩阵，triu()表示取矩阵的上三角元素；

步骤3)中：

根据空间谱理论，信号子空间和噪声子空间具有正交性，当阵列的导向矢量a_l(p)由辐射源的位置p构成时，导向矢量到噪声子空间的投影为零，即满足：

其中，是噪声子空间；

用模来表示：

接着改写为：

联合所有阵列，构建代价函数：

其中，L为阵列的个数。

2.如权利要求1所述的基于卡尔曼滤波器和粒子群优化的单目标直接跟踪方法，其特征在于，步骤1)中，所述多阵列单目标直接跟踪模型包括若干个观测基站和一个目标即辐射源，每个所述观测基站均配备有阵列，将得到的接收信号r_l(t)定义为：

r_l(t)＝a_l(p)s_l(t)+n_l(t)

其中，为导向矢量，d表示阵列中阵元的间距，λ表示波长，M表示阵元数量；u_l＝[x_l,y_l]^T表示第l个观测基站的位置，辐射源的位置为p＝[x,y]^T，u_l(1)和p(1)表示u_l(1)和p(1)的第一个元素，θ_l表示第l个观测基站接收到的辐射源的波达角度；s_l(t)表示第l个观测基站接收到辐射源的回波信号；n_l(t)表示高斯白噪声矢量。

3.如权利要求2所述的基于卡尔曼滤波器和粒子群优化的单目标直接跟踪方法，其特征在于，步骤4)中利用粒子群算法求解直接获得目标该时刻位置的测量值。

4.如权利要求3所述的基于卡尔曼滤波器和粒子群优化的单目标直接跟踪方法，其特征在于，步骤5)中：首先建立卡尔曼滤波器的状态方程和测量方程：

y(t+1)＝Fy(t)+Bh+c(t)

z(t)＝Hy(t)+e(t)

其中，y(t)＝[x₁(t),x₂(t),v₁(t),v₂(t)]表示t时刻的状态向量，x₁(t)、x₂(t)、v₁(t)、v₂(t)分别表示t时刻目标辐射源的X轴方向坐标、Y轴方向坐标、沿X轴方向的速度、沿Y轴方向的速度；h表示控制量，且h＝[0,0]^T；F表示状态转移矩阵、B表示控制量转移矩阵、H表示测量矩阵；c(t)为过程噪声，服从零均值高斯分布，其协方差矩阵是Q(t)；e(t)是测量噪声，也服从零均值高斯分布，其方差为R(t)，卡尔曼滤波器的初值设置为：

其中T_s是跟踪的时间间隔，假设已知t-1时刻的状态向量和状态协方差矩阵通过卡尔曼方程得到t时刻的状态向量和状态协方差矩阵的先验估计：

其中，中的前两个元素为t-1时刻目标的X轴和Y轴坐标；利用目标的位置坐标，获得目标位置的测量值z(t)；

修正先验估计，获得位置的后验估计和协方差矩阵的后验估计：

I表示单位矩阵，K(t)表示卡尔曼增益：

每个时刻都进行卡尔曼滤波，获得更准确的估计结果。