CN110146850B - 用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法 - Google Patents

Abstract

本发明公开一种用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法，属于多基地雷达系统数据融合技术领域。本发明序贯地利用失序量测更新失序量测产生时刻和最新顺序滤波时刻间的所有状态，并利用粒子滤波实现目标跟踪，克服了由于雷达站点数据处理时间不同和通信链路延迟等原因引起的雷达量测失序，从而导致目标融合跟踪性能恶化的问题。本发明判断当前量测是否是顺序量测，若是进行顺序量测更新，否则进行失序量测融合，有效解决了多基地雷达系统中因各雷达站点数据预处理时间不同和通信链路延迟等原因导致的多个任意时序的失序量测的问题，针对多个任意时序的失序量测问题提出了一种普适的解决方法，相比于直接忽略失序量测的融合精度更高。

Description

用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法

技术领域

本发明属于多基地雷达系统数据融合技术领域，特别涉及一种用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法。

背景技术

随着现代环境的越发复杂，单基地雷达的结构空间和功率等可用资源有限，探测能力受到了较大的限制，多基地雷达具有反隐身、抗反辐射导弹、抗干扰和反低空突防的特性，可以很好的适应现代战场环境。近年来，多基地雷达系统已经广泛的应用于海、陆、空等多个领域，在检测、识别、跟踪和成像等多个方面都受到了极大的关注。尤其地，对位于不同位置的雷达接收机进行数据融合，可以有效地提高对目标的定位和跟踪性能。多基地雷达系统数据融合有集中式和分布式两种，集中式融合将各接收机的局部量测送至融合中心进行统一处理，以获得全面一致的全局估计结果，具有最优的融合精度。

然而，多基地雷达系统中各雷达站点相距较远且独立工作，在集中式融合处理中由于各接收机数据预处理时间不同和通信链路延时等原因，雷达站点的量测到达融合中心时失序现象严重，这些失序量测频繁地以任意时序到达融合中心，即多个失序量测或连续到达，或多个失序量测交叉到达，或失序量测与顺序量测交叉到达。由于多基地雷达系统中往往会出现大量的失序量测，若直接忽略这些失序量测，则信息损失严重，导致目标融合跟踪性能严重恶化，因此需要寻找一种特定的方法来处理多基地雷达系统中的失序量测问题。一种理想的方法是每当一个失序量测到达融合中心，便对所有量测进行重新排序重新滤波，然而该方法要求缓存历史量测，这在存储能力有限的多基地雷达系统中是难以满足的。现有技术中，基于线性最小均方误差估计方法提出了一种最优更新方法，有效地解决了任意时序的多个失序量测更新融合问题，但是该方法基于卡尔曼滤波算法，仅适用于线性高斯系统，在多基地雷达系统中该方法跟踪性能退化严重；另外，还有针对多个失序量测提出的一种存储有效的高斯近似融合方法，有效地解决了非线性系统下多个失序量测的更新融合问题，但是该方法仅利用单个均值和协方差矩阵近似表征目标状态的概率密度函数，当多基地雷达系统面临的目标状态是多模分布时，该方法跟踪精度差。

发明内容

本发明的目的在于解决上述技术问题，提出了一种用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法，从确定的贝叶斯理论出发，推导了任意时序的多个失序量测的更新融合方法，并提出了粒子滤波实现方法，有效地解决了多基地雷达系统中由于各接收机数据预处理时间不同和通信链路延时等原因导致的多个任意时序的失序量测问题。

一种用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法，应用于多基地雷达系统，所述多基地雷达系统包括数个雷达，融合中心接收雷达量测，所述方法包括以下步骤：

S1、融合中心初始化，初始化粒子样本及其相应权值，初始化更新时刻t_k＝0，其中，时间索引k＝0；

S2、所述融合中心获取量测，判断当前量测的时戳t_z是否小于或等于所述多基地雷达系统的观测总时间t_total，若是则进入步骤S3；否则结束；

S3、判断当前量测的时戳t_z是否大于或等于上一更新时刻t_k，若是则进入步骤S4；否则进入步骤S5；

S4、当前量测为顺序量测，将时戳t_z记为t_k+1，采用标准粒子滤波方法获得t_k+1时刻的顺序滤波后验状态，将表征t_k+1时刻顺序滤波后验状态的粒子样本进行存储，令k＝k+1，返回所述步骤S2；

S5、当前量测为失序量测，将时戳t_z记为τ，通过判断t_k-l＜τ＜t_k-l+1确定失序量测的延迟步数l；

S6、采用标准粒子滤波方法获得τ时刻的顺序滤波后验状态，将表征τ时刻顺序滤波后验状态的粒子样本进行存储；对于时刻t_k-l+1和时刻t_k之间的所有中间时刻κ_j，初始化j，进入步骤S7；其中，j＝k-l+1,k-l+2,…,k，κ_j＝t_k-l+1,t_k-l+2,…t_k；

S7、采用平滑方法获得τ时刻的平滑概率密度函数；

S8、根据所述平滑概率密度函数和当前失序量测的似然函数，得到失序量测的异步似然函数；

S9、根据所述失序量测的异步似然函数更新κ_j时刻的顺序滤波后验状态，获得κ_j时刻融合了当前失序量测的融合后验状态，令j＝j+1；

S10、判断κ_j是否大于t_k，即j是否大于k，若是则进入步骤S11；否则返回步骤S7；

S11、将κ_j时刻的所述融合后验状态代替κ_j时刻的顺序滤波状态后验状态，将表征κ_j时刻融合后验状态的粒子样本代替表征κ_j时刻顺序滤波后验状态的粒子样本，将所述融合中心获得的后验状态按照对应的时戳顺序进行排列，即后验状态集合对应的时戳集合为[t₁,t₁,…,t_k-l,τ,t_k-l+1,…,t_k]，返回所述步骤S2。

进一步地，所述步骤S1中初始化粒子样本及其相应权值的步骤包括：

产生Q个初始粒子样本

服从

对应的权值为

其中，i₀＝1,2,…,Q表示粒子样本标号，Q表示粒子数目。

进一步地，步骤S4包括：

已知粒子样本

对于所求粒子样本

有

其中，～表示采样，即从高斯分布

中采样得到

p(x_k+1|x_k)表示状态转移函数，表示t_k时刻的目标状态x_k转移到t_k+1时刻的目标状态x_k+1的条件概率函数；z_k+1表示t_k+1时刻达到融合中心的量测，p(z_k+1|x_k+1)表示似然函数，

表示粒子

对应的权值；

目标顺序滤波后验状态为

进一步地，所述步骤S5包括：

当前量测为失序量测，时戳满足t_k-l＜τ＜t_k-l+1，其中，l表示当前失序量测的延迟步数，为正整数。

进一步地，所述步骤S6包括：

对于失序量测产生时刻τ，已知粒子样本

所求粒子样本为

有

则τ时刻目标顺序滤波后验状态为

对于时刻t_k-l+1和时刻t_k之间的所有中间时刻κ_j＝t_k-l+1,t_k-l+2,…t_k，其中，j＝k-l+1,k-l+2,…,k，初始化j＝k-l+1，进入步骤S7。

进一步地，所述步骤S7包括：

k时刻的顺序滤波先验概率密度函数和顺序滤波后验概率密度函数的粒子表达式为

其中，δ(·)表示标准狄利赫利函数；

当κ_j＝t_k-l+1时，τ时刻的平滑概率密度函数为

其中，x_j表示κ_j时刻的目标状态，z₁:_j表示从t₀时刻直到κ_j时刻的量测集合；将τ时刻和t_k-l+1时刻相应的顺序滤波先验概率密度函数的粒子表达式代入上式τ时刻的平滑概率密度函数，得到τ时刻的平滑概率密度函数的粒子表达式

其中，

表示在量测z₁:_j和目标状态

的条件下粒子

对应的平滑权值；当κ_j＝t_k-l+1时，τ时刻平滑概率密度函数p(x(τ)|x_j,z_1:j)表示为粒子样本

其中

其中，∝表示正比于；

当κ_j＞t_k-l+1时，τ时刻的平滑概率密度函数为

p(x(τ)|x_j,z_1:j)＝∫p(x(τ),x_j-1|x_j,z_1:j)dx_j-1

＝∫p(x(τ)|x_j-1,z_1:j-1)p(x_j-1|x_j,z_1:j-1)dx_j-1

其中，p(x(τ)|x_j-1,z_1:j-1)表示在量测z_1:j-1和目标状态x_j-1的条件下τ时刻的平滑概率密度函数，p(x_j-1|x_j,z_1:j-1)表示为

将κ_j-1时刻的顺序滤波后验概率密度函数的粒子表达式代入上式中p(x_j-1|x_j,z_1:j-1)，得到其粒子表达式

将在量测z_1:j-1和目标状态条件x_j-1下τ时刻的平滑概率密度函数的粒子表达式和上式

代入κ_j＞t_k-l+1时τ时刻的平滑概率密度函数，得到τ时刻平滑概率密度函数的粒子表达式

当κ_j＞t_k-l+1时，τ时刻平滑概率密度函数p(x(τ)|x_j,z_1:j)表示为粒子样本

其中

进一步地，所述步骤S8包括：

根据所述步骤S7获得的平滑概率密度函数p(x(τ)|x_j,z_1:j)和当前失序量测的似然函数p(z(τ)|x(τ))，得到失序量测的异步似然函数p(z(τ)|x_j,z_1:j)，即

p(z(τ)|x_j,z_1:j)＝∫p(z(τ)|x(τ))p(x(τ)|x_j,z_1:j)dx(τ)

将τ时刻平滑概率密度函数p(x(τ)|x_j,z_1:j)的粒子表达式代入上式异步似然函数中，得到异步似然函数p(z(τ)|x_j,z_1:j)的粒子表达式

其中，

为表征平滑概率密度函数的粒子样本对应的平滑权值，

表示失序量测的似然函数。

进一步地，所述步骤S9包括：

采用所述异步似然函数对κ_j时刻的顺序滤波后验概率密度函数进行更新，得到融合了失序量测的融合后验概率密度函数

将所述异步似然函数的粒子表达式和κ_j时刻的顺序滤波后验概率密度函数的粒子表达式代入所述融合后验概率密度函数，得到融合后验概率密度函数的粒子表达式

由所述融合后验概率密度函数的粒子表达式，得到κ_j时刻的粒子的融合后验权值表达式

经过重采样后，获得融合了该失序量测后的κ_j时刻的融合后验状态

进一步地，所述步骤S1之前，还包括初始化系统参数，包括：多基地雷达系统监测平面大小、具有发射器的单基地雷达总个数、具有接收器的单基地雷达总个数、观测总时间及目标初始状态。

本发明的有益效果：本发明提供了一种用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法，本发明的方法序贯地利用失序量测更新失序量测产生时刻和最新顺序滤波时刻之间的所有状态，并利用粒子滤波实现目标跟踪。首先判断当前量测是否是顺序量测，若是则进行顺序量测更新，即利用标准粒子滤波算法获得顺序滤波后验状态，否则进行失序量测融合，即对于失序量测产生时刻，利用标准粒子滤波算法获得失序量测时戳处的顺序滤波后验状态，对于失序量测产生时刻和最新顺序滤波时刻之间的所有中间时刻，先利用基于贝叶斯框架的平滑方法求解失序量测产生时刻处的平滑概率密度函数，然后联合平滑概率密度函数和失序量测的似然函数，获得失序量测的异步似然函数，最后利用失序量测的异步似然函数对这些中间时刻的顺序滤波后验状态进行更新，获得相应的融合后验状态，有效地解决了多基地雷达系统中由于各接收机数据预处理时间不同和通信链路延迟等原因导致的多个任意时序的失序量测问题；并且具备以下优点：

(1)与现有的非线性系统下多个任意时序的失序量测数据融合方法相比，本发明的方法从确定的贝叶斯理论出发，融合精度高；

(2)本发明的方法利用粒子样本近似目标状态的概率密度函数，可以很好的表征具有多模分布的目标，因此可以很好的适应多基地雷达系统面临的复杂环境，且当粒子数目趋于无穷大时，跟踪精度非常接近于理想顺序重滤波方法的跟踪精度；

(3)本发明的方法适用于以任意时序到达融合中心的多个失序量测融合问题，具体地，来自各雷达接收机的多个失序量测的时序可以是连续的、相互交叉的、和顺序量测相互交叉的；

(4)本发明的方法操作简单，只需要对多个失序量测进行序贯处理，利用每个失序量测的异步似然函数对中间时刻的顺序滤波后验状态进行更新。

附图说明

图1为本发明实施例的流程图。

图2为本发明实施例的多基地雷达系统集中式融合示意图。

图3为本发明实施例的多基地雷达系统场景与目标运动轨迹示意图。

图4为本发明实施例的各雷达接收机的量测到达融合中心的时序示意图。

图5为本发明实施例的多个任意时序失序量测的时序示意图。

图6为本发明方法与理想顺序重滤波方法、忽略失序量测方法、高斯近似融合方法的跟踪精度对比图。

具体实施方式

下面结合附图对本发明的实施例做进一步的说明。

请参阅图1，本发明提供了一种用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法，应用于多基地雷达系统，多基地雷达系统包括数个雷达，融合中心接收雷达量测，如图2所示。

本实施例中，本发明首先初始化系统参数，包括：多基地雷达系统观测平面大小；具有发射器的单基地雷达总个数M；具有接收器的单基地雷达总个数N；雷达系统观测总时间t_total；目标的初始状态

其中，(x(0),y(0))表示目标的初始位置，

表示目标的初始速度；目标初始状态偏差服从高斯分布

如图3所示，多基地雷达系统由M＝2部发射机和N＝4部接收机组成，雷达系统对一个观测平面大小为200km×200km的二维平面区域一运动目标进行监视，目标的初始状态x(0)＝(50,1.5,30,1.7)′；即目标的初始位置是(50,30)，并以(1.5,1.7)的速度运动，图3给出了雷达系统观测总时间t_total＝80s的目标运动轨迹示意图。目标初始状态偏差服从高斯分布

并且C₀＝diag(0.5,0,0.5,0)。

本实施例中，接收机2、3、4由于数据预处理时间不同和通信链路延时等原因，所有量测均延迟某一周期数；本实施例设置接收机2的所有量测均延迟四步，接收机3的所有量测均延迟一步，接收机4的所有量测均延迟一步，如图4所示(图4给出了前20s到达融合中心的量测时序示意图)。

本发明方法通过以下步骤实现：

S1、融合中心初始化，初始化粒子样本及其相应权值，初始化更新时刻t_k＝0，其中，时间索引k＝0。

本实施例中，产生Q个初始粒子样本

服从

对应的权值为

其中，i₀＝1,2,…,Q表示粒子样本标号，Q表示粒子数目。

S2、融合中心获取量测，判断当前量测的时戳t_z是否小于或等于多基地雷达系统的观测总时间t_total，若是则进入步骤S3；否则结束。

本实施例中，步骤S2通过以下子步骤实现：

S21、融合中心获取量测，时戳为t_z；

S22、判断当前量测的时戳t_z是否小于或等于多基地雷达系统的观测总时间t_total，若是则进入步骤S3；否则结束。

S3、判断当前量测的时戳t_z是否大于或等于上一更新时刻t_k，若是则进入步骤S4；否则进入步骤S5。

本实施例中，若当前量测的时戳t_z大于或等于上一个更新时刻，说明当前量测是顺序量测，则执行步骤S4，即利用标准粒子滤波算法进行状态更新；否则，当前量测的时戳t_z小于上一个更新时刻，说明当前量测是失序量测，则执行步骤S5，即进行失序量测融合。

S4、当前量测为顺序量测，将时戳t_z记为t_k+1，采用标准粒子滤波方法获得t_k+1时刻的顺序滤波后验状态，将表征t_k+1时刻顺序滤波后验状态的粒子样本进行存储，令k＝k+1，返回步骤S2。

本实施例中，已知粒子样本

对于所求粒子样本

有

其中，～表示采样，即从高斯分布

中采样得到

p(x_k+1|x_k)表示状态转移函数，表示t_k时刻的目标状态x_k转移到t_k+1时刻的目标状态x_k+1的条件概率函数；z_k+1表示t_k+1时刻达到融合中心的量测，p(z_k+1|x_k+1)表示似然函数，

表示粒子

对应的权值。

目标顺序滤波后验状态为

S5、当前量测为失序量测，将时戳t_z记为τ，通过判断t_k-l＜τ＜t_k-l+1确定失序量测的延迟步数l。

本实施例中，当前量测为失序量测，时戳满足t_k-l＜τ＜t_k-l+1，其中，l表示当前失序量测的延迟步数，为正整数。

S6、采用标准粒子滤波方法获得τ时刻的顺序滤波后验状态，将表征τ时刻顺序滤波后验状态的粒子样本进行存储；对于时刻t_k-l+1和时刻t_k之间的所有中间时刻κ_j，初始化j，进入步骤S7；其中，j＝k-l+1,k-l+2,…,k，κ_j＝t_k-l+1,t_k-l+2,…t_k。

本实施例中，对于失序量测产生时刻τ，已知粒子样本

所求粒子样本为

有

则该失序量测产生时刻τ时刻的目标顺序滤波后验状态为

对于时刻t_k-l+1和时刻t_k之间的所有中间时刻κ_j＝t_k-l+1,t_k-l+2,…t_k，其中，j＝k-l+1,k-l+2,…,k，初始化j＝k-l+1，进入步骤S7，从时刻t_k-l+1开始直到时刻t_k，循环执行步骤S7-S9。如图5所示为多个任意时序失序量测的时序示意图。

S7、采用平滑方法获得τ时刻的平滑概率密度函数。

本实施例中，k时刻的顺序滤波先验概率密度函数和顺序滤波后验概率密度函数的粒子表达式为

其中，δ(·)表示标准狄利赫利函数。

采用平滑方法获得τ时刻的平滑概率密度函数，分为κ_j＝t_k-l+1(即初始情况)和κ_j＞t_k-l+1两种情况。

(1)当κ_j＝t_k-l+1时，τ时刻的平滑概率密度函数为

其中，x_j表示κ_j时刻的目标状态，z₁:_j表示从t₀时刻直到κ_j时刻的量测集合；将τ时刻和t_k-l+1时刻相应的顺序滤波先验概率密度函数的粒子表达式代入式(10)，得到τ时刻的平滑概率密度函数(10)的粒子表达式

其中，

表示在量测z₁:_j和目标状态

的条件下粒子

对应的平滑权值。

因此，当κ_j＝t_k-l+1时，τ时刻平滑概率密度函数p(x(τ)|x_j,z_1:j)表示为粒子样本

其中，

由步骤S6得到，对应的平滑权值由式(11)可得，即

其中，符号∝表示正比于。

(2)当κ_j＞t_k-l+1时，τ时刻的平滑概率密度函数为

其中，p(x(τ)|x_j-1,z_1:j-1)表示在量测z_1:j-1和目标状态x_j-1的条件下τ时刻的平滑概率密度函数，由上一次迭代获得；p(x_j-1|x_j,z_1:j-1)表示为

将κ_j-1时刻的顺序滤波后验概率密度函数的粒子表达式代入式(14)，得到式(14)的粒子表达式

将在量测z_1:j-1和目标状态条件x_j-1下τ时刻的平滑概率密度函数的粒子表达式和式(15)代入式(13)，得到τ时刻平滑概率密度函数的粒子表达式

因此，κ_j＞t_k-l+1时，τ时刻平滑概率密度函数p(x(τ)|x_j,z_1:j)表示为粒子样本

其中，

由步骤S6得到，其对应的平滑权值由式(16)可得，即

通过将κ_j分为κ_j＝t_k-l+1和κ_j＞t_k-l+1两种情况，求得所有τ时刻的平滑概率密度函数。

S8、根据平滑概率密度函数和当前失序量测的似然函数，得到失序量测的异步似然函数。

本实施例中，根据步骤S7获得的平滑概率密度函数p(x(τ)|x_j,z_1:j)和失序量测的似然函数p(z(τ)|x(τ))，得到该失序量测的异步似然函数p(z(τ)|x_j,z_1:j)，即

p(z(τ)|x_j,z_1:j)＝∫p(z(τ)|x(τ))p(x(τ)|x_j,z_1:j)dx(τ) (18)

将τ时刻平滑概率密度函数p(x(τ)|x_j,z_1:j)的粒子表达式代入式(18)，得到异步似然函数p(z(τ)|x_j,z_1:j)的粒子表达式

其中，

为表征平滑概率密度函数的粒子样本对应的平滑权值，由步骤S7可得；

表示失序量测的似然函数，由量测方程得到。

S9、根据失序量测的异步似然函数更新κ_j时刻的顺序滤波后验状态，获得κ_j时刻融合了当前失序量测的融合后验状态，令j＝j+1。

本实施例中，采用异步似然函数对κ_j时刻的顺序滤波后验概率密度函数进行更新，得到融合了失序量测的融合后验概率密度函数

将式(19)和κ_j时刻的顺序滤波后验概率密度函数的粒子表达式代入式(20)，得到融合后验概率密度函数的粒子表达式

由式(21)可得，κ_j时刻的粒子的融合后验权值表达式

经过重采样后，获得融合了该失序量测后的κ_j时刻的融合后验状态

完成该κ_j时刻的S7-S9的失序量测融合过程，令j＝j+1，即进入下一κ_j时刻的失序量测融合，循环S7-S9，直到κ_j＝t_k。

S10、判断κ_j是否大于t_k，即j是否大于k，若是则进入步骤S11；否则返回步骤S7。

本实施例中，判断κ_j是否大于t_k，j是否大于k，即是否到达了迭代终止条件。

S11、将κ_j时刻的融合后验状态代替κ_j时刻的顺序滤波状态后验状态，将表征κ_j时刻融合后验状态的粒子样本代替表征κ_j时刻顺序滤波后验状态的粒子样本，将融合中心获得的后验状态按照对应的时戳顺序进行排列，即后验状态集合对应的时戳集合为[t₁,t₁,…,t_k-l,τ,t_k-l+1,…,t_k]，返回步骤S2。

本实施例中，用κ_j时刻的融合后验状态

代替κ_j时刻的顺序滤波后验状态

得到新的状态估计集合

其对应的时戳集合为

[t₁,t₂,…,t_k-l,τ,t_k-l+1,…,t_k] (25)

如图6所示，为多基地雷达系统中以任意时序到达融合中心的多个失序量测场景下，忽略失序量测方法、理想顺序重滤波方法、高斯近似融合方法和本发明方法(粒子近似融合方法)的跟踪精度对比图。其中，跟踪精度通过跟踪轨迹和真实轨迹之间的均方根误差进行描述。可以得出，由于忽略失序量测方法忽略了大量的失序量测，导致信息损失严重，所以它有最差的跟踪性能；理想顺序重滤波方法假设所有量测都是顺序量测，是一种理想的时序处理方法，其有着最优的跟踪性能，在本实施例中设置为性能下界；而粒子近似融合方法，相比于忽略失序量测方法有着更好的跟踪结果，这说明本发明方法可以有效的处理多基地雷达系统中多个任意时序的失序量测问题；同时，粒子近似融合方法相较于高斯近似融合方法有着更小的跟踪误差，这说明本发明方法在跟踪精度方面优于现存的可用于多基地雷达系统中的多个失序量测融合算法；另外，粒子近似融合方法的均方根误差曲线非常接近于理想顺序重滤波方法，进一步说明了本方法的有效性。

综上，本发明通过综合考虑目标跟踪精度和存储需求等问题，提出的一种用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法，克服了由于雷达站点数据预处理时间不同和通信链路延迟等原因引起的雷达量测失序，从而导致目标融合跟踪性能恶化的问题。本发明的特点是序贯地利用失序量测更新失序量测产生时刻和最新顺序滤波时刻之间的所有状态，并利用粒子滤波实现目标跟踪。首先判断当前量测是否是顺序量测，若是则进行顺序量测更新，即利用标准粒子滤波算法获得顺序滤波后验状态，否则进行失序量测融合，即对于失序量测产生时刻，利用标准粒子滤波算法获得失序量测时戳处的顺序滤波后验状态，对于失序量测产生时刻和最新顺序滤波时刻之间的所有中间时刻，先利用基于贝叶斯框架的平滑方法求解失序量测产生时刻处的平滑概率密度函数，然后联合平滑概率密度函数和失序量测的似然函数，获得失序量测的异步似然函数，最后利用失序量测的异步似然函数对这些中间时刻的顺序滤波后验状态进行更新，获得相应的融合后验状态。有效地解决了多基地雷达系统中由于各雷达站点数据预处理时间不同和通信链路延迟等原因导致的多个任意时序的失序量测的问题，本发明针对多基地雷达系统中多个任意时序的失序量测问题提出了一种普适的解决方法，相比于忽略失序量测方法和高斯近似融合方法存在优势。

本领域的普通技术人员将会意识到，这里的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (9)

1.一种用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法，应用于多基地雷达系统，其特征在于，所述多基地雷达系统包括数个雷达，融合中心接收雷达量测，所述方法包括以下步骤：

S1、融合中心初始化，初始化粒子样本及其相应权值，初始化更新时刻t_k＝0，其中，时间索引k＝0；

S2、所述融合中心获取量测，判断当前量测的时戳t_z是否小于或等于所述多基地雷达系统的观测总时间t_total，若是则进入步骤S3；否则结束；

S3、判断当前量测的时戳t_z是否大于或等于上一更新时刻t_k，若是则进入步骤S4；否则进入步骤S5；

S4、当前量测为顺序量测，将时戳t_z记为t_k+1，采用标准粒子滤波方法获得t_k+1时刻的顺序滤波后验状态，将表征t_k+1时刻顺序滤波后验状态的粒子样本进行存储，令k＝k+1，返回所述步骤S2；

S5、当前量测为失序量测，将时戳t_z记为τ，通过判断t_k-l<τ<t_k-l+1确定失序量测的延迟步数l；

S6、采用标准粒子滤波方法获得τ时刻的顺序滤波后验状态，将表征τ时刻顺序滤波后验状态的粒子样本进行存储；对于时刻t_k-l+1和时刻t_k之间的所有中间时刻κ_j，初始化j，进入步骤S7；其中，j＝k-l+1,k-l+2,…,k，κ_j＝t_k-l+1,t_k-l+2,…t_k；

S7、采用平滑方法获得τ时刻的平滑概率密度函数；

S8、根据所述平滑概率密度函数和当前失序量测的似然函数，得到失序量测的异步似然函数；

S9、根据所述失序量测的异步似然函数更新κ_j时刻的顺序滤波后验状态，获得κ_j时刻融合了当前失序量测的融合后验状态，令j＝j+1；

S10、判断κ_j是否大于t_k，即j是否大于k，若是则进入步骤S11；否则返回步骤S7；

S11、将κ_j时刻的所述融合后验状态代替κ_j时刻的顺序滤波状态后验状态，将表征κ_j时刻融合后验状态的粒子样本代替表征κ_j时刻顺序滤波后验状态的粒子样本，将所述融合中心获得的后验状态按照对应的时戳顺序进行排列，即后验状态集合对应的时戳集合为[t₁,t₁,…,t_k-l,τ,t_k-l+1,…,t_k]，返回所述步骤S2。

2.如权利要求1所述的用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法，其特征在于，所述步骤S1中初始化粒子样本及其相应权值的步骤包括：

产生Q个初始粒子样本

服从

对应的权值为

其中，i₀＝1,2,…,Q表示粒子样本标号，Q表示粒子数目。

3.如权利要求2所述的用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法，其特征在于，步骤S4包括：

已知粒子样本

对于所求粒子样本

有

其中，～表示采样，即从高斯分布

中采样得到

p(x_k+1|x_k)表示状态转移函数，表示t_k时刻的目标状态x_k转移到t_k+1时刻的目标状态x_k+1的条件概率函数；z_k+1表示t_k+1时刻达到融合中心的量测，p(z_k+1|x_k+1)表示似然函数，

表示粒子

对应的权值；i_k表示粒子样本标号，取值范围为1,2,…,Q；

目标顺序滤波后验状态为

4.如权利要求3所述的用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法，其特征在于，所述步骤S5包括：

当前量测为失序量测，时戳满足t_k-l<τ<t_k-l+1，其中，l表示当前失序量测的延迟步数，为正整数。

5.如权利要求3所述的用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法，其特征在于，所述步骤S6包括：

对于失序量测产生时刻τ，已知粒子样本

所求粒子样本为

有

则τ时刻目标顺序滤波后验状态为

对于时刻t_k-l+1和时刻t_k之间的所有中间时刻κ_j＝t_k-l+1,t_k-l+2,…t_k，其中，j＝k-l+1,k-l+2,…,k，初始化j＝k-l+1，进入步骤S7。

6.如权利要求5所述的用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法，其特征在于，所述步骤S7包括：

k时刻的顺序滤波先验概率密度函数和顺序滤波后验概率密度函数的粒子表达式为

其中，δ(·)表示标准狄利赫利函数；

当κ_j＝t_k-l+1时，τ时刻的平滑概率密度函数为

其中，x_j表示κ_j时刻的目标状态，z_1:j表示从t₀时刻直到κ_j时刻的量测集合；将τ时刻和t_k-l+1时刻相应的顺序滤波先验概率密度函数的粒子表达式代入上式τ时刻的平滑概率密度函数，得到τ时刻的平滑概率密度函数的粒子表达式

其中，

表示在量测z_1:j和目标状态

的条件下粒子

对应的平滑权值；当κ_j＝t_k-l+1时，τ时刻平滑概率密度函数p(x(τ)|x_j,z_1:j)表示为粒子样本

其中

其中，∝表示正比于；

当κ_j>t_k-l+1时，τ时刻的平滑概率密度函数为

其中，p(x(τ)|x_j-1,z_1:j-1)表示在量测z_1:j-1和目标状态x_j-1的条件下τ时刻的平滑概率密度函数，p(x_j-1|x_j,z_1:j-1)表示为

将κ_j-1时刻的顺序滤波后验概率密度函数的粒子表达式代入上式中p(x_j-1|x_j,z_1:j-1)，得到其粒子表达式

将在量测z_1:j-1和目标状态条件x_j-1下τ时刻的平滑概率密度函数的粒子表达式和上式

代入κ_j>t_k-l+1时τ时刻的平滑概率密度函数，得到τ时刻平滑概率密度函数的粒子表达式

当κ_j>t_k-l+1时，τ时刻平滑概率密度函数p(x(τ)|x_j,z_1:j)表示为粒子样本

其中

7.如权利要求6所述的用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法，其特征在于，所述步骤S8包括：

根据所述步骤S7获得的平滑概率密度函数p(x(τ)|x_j,z_1:j)和当前失序量测的似然函数p(z(τ)|x(τ))，得到失序量测的异步似然函数p(z(τ)|x_j,z_1:j)，即

p(z(τ)|x_j,z_1:j)＝∫p(z(τ)|x(τ))p(x(τ)|x_j,z_1:j)dx(τ)

将τ时刻平滑概率密度函数p(x(τ)|x_j,z_1:j)的粒子表达式代入上式异步似然函数中，得到异步似然函数p(z(τ)|x_j,z_1:j)的粒子表达式

其中，

为表征平滑概率密度函数的粒子样本对应的平滑权值，

表示失序量测的似然函数。

8.如权利要求7所述的用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法，其特征在于，所述步骤S9包括：

采用所述异步似然函数对κ_j时刻的顺序滤波后验概率密度函数进行更新，得到融合了失序量测的融合后验概率密度函数

将所述异步似然函数的粒子表达式和κ_j时刻的顺序滤波后验概率密度函数的粒子表达式代入所述融合后验概率密度函数，得到融合后验概率密度函数的粒子表达式

由所述融合后验概率密度函数的粒子表达式，得到κ_j时刻的粒子的融合后验权值表达式

经过重采样后，获得融合了该失序量测后的κ_j时刻的融合后验状态

9.如权利要求1-8任一项所述的用于多基地雷达失序量测融合的粒子滤波集中式跟踪方法，其特征在于，所述步骤S1之前，还包括初始化系统参数，包括：多基地雷达系统监测平面大小、具有发射器的单基地雷达总个数、具有接收器的单基地雷达总个数、观测总时间及目标初始状态。

Images