CN106154274A - 一种水下小目标多基地融合跟踪方法 - Google Patents

Abstract

本发明提供了一种水下小目标多基地融合跟踪方法，所述方法包括：步骤1)对t_k时刻多基地系统的目标量测集合Z_k进行重采样，得到重采样随机粒子集合，通过对所述重采样随机粒子的分布进行量化和融合处理，得到t_k时刻的等效量测z_k；步骤2)通过混合粒子滤波对t_k时刻小目标状态进行估计；所述混合粒子包括预测粒子和目标量测采样粒子；所述预测粒子为用t_k-1时刻的随机粒子预测得到的t_k时刻的随机粒子；所述目标量测采样粒子根据所述等效量测z_k及先验概率分布对等效量测进行随机采样得到；并对混合粒子进行重采样；步骤3)令k＝k+1，转入步骤1)；直至目标跟踪结束。本发明的方法提高了多基地水下小目标的跟踪精度。

Description

一种水下小目标多基地融合跟踪方法

技术领域

本发明涉及水下小目标的跟踪技术领域，特别涉及一种水下小目标多基地融合跟踪方法。

背景技术

冷战时期，各国处在战争警戒状态，海上探测和防御的对象主要是敌对国家的舰船和潜艇等大型目标。随着冷战的结束，特别是前苏联解体后，水下武器装备小型化迅速发展，蛙人、水下运载器及水下机器人等技术装备日益成熟，这类具有隐蔽性好，破坏力强等明显“非对称”优势的攻击方式颇受恐怖主义分子的青睐，已成为恐怖主义分子进行恐怖活动的重要方式。目前，大范围的水域监控，如岛礁防护、港池锚泊区防护等已成为水下安保趋势的主流，传统的重点水域的“点防护”方式已经难以适应水下安保的需求。为实现大范围水域的“面防护”，往往需要多基地组网的工作方式实现。

近年来，国内外研究人员对蛙人等水下小目标的探测及识别进行了相关研究，但针对水下小目标的多基地跟踪技术还没有进行系统地研究。Anders等(参考文献[1]：ANDERSYAAKOV BAR-SHALOM,ODDVARHALLINGSTAD JOHN GLATTETRE.Multitarget Multisensor Tracking in the Presenceof Wakes[J].JOURNAL OF ADVANCES IN INFORMATION FUSION,2009,4(2):117-145)研究了多传感器下开式蛙人呼吸形成气泡背景下的目标跟踪问题，改进了传统的数据关联算法，有效抑制了尾流背景干扰量测对数据关联精度的影响，并将序贯融合用于水下小目标多基地的融合估计，目标状态估计采用卡尔曼滤波算法。文献2(Dietmar Stiller and Bernd Nützel.Detection and tracking of divers[C].UDTEurope,2009)在单基地主动声纳CW信号和LFM信号联合探测后融合数据处理的基础上，通过多假设跟踪(MHT)进行蛙人目标的状态估计，有效抑制了虚假航迹，提高了目标航迹的连续性。文献3(K.Clare Xu,Y.Jay Guo,Eryk Dutkiewicz.SwimmerTracking with Underwater Acoustic Networks[C].International Symposium onCommunications and Information Technologies(ISCIT),2010:1045-1049)利用水下传感器网络在定位精度的优势，将球形插值算法用于水面游泳者的定位。

针对复杂背景干扰下的水下目标跟踪，EvanHanusa(参考文献[4]E.Hanusa,D.W.Krout.Posterior distribution preprocessing with theJPDA algorithm:PACsim data set[C].15th International Conference on Information Fusion(FUSION),2012:2593-2599.、参考文献[5]E.Hanusa,D.Krout,M.R.Gupta.Clutter rejection by clustering likelihood-basedsimilarities[C].14th International Conference on Information Fusion(FUSION),2011:1-6.和参考文献[6]E.Hanusa,D.Krout,M.R.Gupta.Contact Clustering and Fusion forPreprocessing Multistatic Active Sonar Data[C]//16th International Conference onInformation Fusion,2013:522-529.)给出了利用目标量测的信噪比、幅度等特征信息，通过多传感器量测聚类融合得到等效目标观测。在假设多传感器系统各声纳节点目标观测数据比杂波、噪声等背景干扰量测一致性更优的前提下，RAMONAGEORGESCU等(参考文献[7]R.Georgescu,P.Willett,S.Marano,V.Matta.Predetection fusion in large sensor networks with unknown target locations[J].JOURNALOF ADVANCES IN INFORMATION FUSION,2012,7(1):61-77.、参考文献[8]R.Georgescu,P.Willett.Random finite set Markov chain Monte Carlo predetectionfusion[C]//14th International Conference on Information Fusion(FUSION),2011:1-8.和参考文献[9]R.Georgescu,P.Willett.Predetection fusion with Doppler measure-mentsand amplitude information[J].IEEE Journal of Oceanic Engineering 2012,37(1):56-65.)提出了多基地系统的预检测量测融合方法，该方法可较好地抑制虚假量测对目标量测融合的干扰。在水下小目标多基地探测系统中，具有混合干扰量测在各声纳节点的稳定性弱于目标观测特点，上述方法没有有效利用这一特点，导致多基地系统的量测融合方法过程复杂，计算量大。

粒子滤波算法在处理非线性、非高斯问题的状态估计方面有独特优势，利用带有一定权值的随机采样粒子的加权组合完成目标状态估计，可有效地完成非线性、非高斯系统下的目标状态估计。但粒子滤波的重采样过程引起“粒子贫化”，影响估计性能甚至导致发散(参考文献[10]Jaechan Lim.Performance Degradation Due toParticle Impoverishment in Particle Filtering[J].Journal of Electrical Engineering&Technology.2014,9(6):2107-2113.)。为了解决重采样带来的粒子贫化问题，Gordon等人提出对每个样本点增加高斯扰动，将马尔可夫链蒙特卡罗MCMC(Markov ChainMonte Carlo)引入粒子滤波算法中，该方法在改善粒子多样性的同时，也极大地增加了算法的运行时间(参考文献[11]GORDON N J,SALMOND D J,SMITH A F M.Novel approach to nonlinear/non-Gaussian Bayesian state estimation[J].IEEEProceedings of Radar and Signal Processing,1993,140(2):107-113.)。Tiancheng Li等(参考文献[12]Tiancheng Li,Tariq Pervez Sattar,Shudong Sun.Deterministic resampling:Unbiased sampling to avoid sample impoverishment in particle filters[J].SignalProcessing,2012,92(7):1637–1645)通过确定重采样的方法，避免了低权重随机粒子的剔除，缓解了传统粒子滤波过程中的“粒子贫化”现象。在粒子聚集的一定邻域内抽取产生新的粒子，也可降低多样性损失带来的粒子贫化现象。传统粒子滤波在水下小目标状态估计过程中，随机粒子存在可能偏离目标真实位置的问题，导致目标状态估计的精度降低。

发明内容

本发明的目的在于克服多基地测量中存在的虚假量测对目标量测融合的干扰问题，克服多基地传统粒子滤波在水下小目标状态估计过程中，存在的随机粒子可能偏离目标真实位置的问题，提出了一种水下小目标多基地融合跟踪方法，该方法首先对t_k时刻多基地系统的目标量测集合Z_k进行处理得到等效量测；然后用t_k-1时刻的随机粒子预测t_k时刻的随机粒子的基础上，进一步对等效量测值进行随机采样，得到采样粒子，并将预测粒子和采样粒子进行融合，得到目标的状态估计值，由此提高了了多基地水下小目标的跟踪精度。

为了实现上述目的，本发明提供了一种水下小目标多基地融合跟踪方法，所述方法包括：

步骤1)对t_k时刻多基地系统的目标量测集合Z_k进行重采样，得到重采样随机粒子集合通过对所述重采样随机粒子的分布进行量化和融合处理，得到t_k时刻的等效量测z_k；

步骤2)通过混合粒子滤波对t_k时刻小目标状态进行估计；所述混合粒子包括预测粒子和目标量测采样粒子；所述预测粒子为用t_k-1时刻的随机粒子预测得到的t_k时刻的随机粒子；所述目标量测采样粒子根据所述等效量测z_k及先验概率分布对等效量测进行随机采样得到；并对混合粒子进行重采样；

步骤3)令k＝k+1，转入步骤1)；直至目标跟踪结束。

上述技术方案中，所述步骤1)具体包括：

步骤1-1)在t_k时刻，对所有小目标探测声纳的接收数据进行预处理获取多基地系统的目标量测集合Z_k；

所述预处理包括：带通采样、降采样和时空配准；得到t_k时刻多基地系统的混响干扰量测和目标量测集合为：

$Z_k=\{z_{\mathrm{1,1}},z_{\mathrm{1,2}},\·\·\·,z_{1,N_1},\·\·\·,z_{M,N_M}\};$

其中，M为多基地系统中声纳节点数目；N_i为第i个声纳节点得到的量测数目；

步骤1-2)根据多基地系统各声纳节点在量测位置处的观测噪声协方差，对目标量测集合Z_k进行重采样，得到重采样随机粒子的集合其中

步骤1-3)对所述重采样随机粒子的分布进行量化；

首先对所述重采样随机粒子所在的平面进行了P×Q栅格剖分；将所述重采样随机粒子所在的平面划分为边长h的栅格，则序号(i,j)的栅格π(i,j)左下角的位置(x_i,y_j)(x_i,y_j)为：

$\left\{\begin{array}{c}{x}_i=x_0+i*h\\ y_j=y_0+j*h\end{array}\right.---\left(3\right)$

式中，(x₀,y₀)为栅格化平面坐标起点；

计算落入每个栅格内随机粒子的数目；若随机粒子的二维坐标为(x,y)，则该随机样本所在栅格π(p,q)的栅格序号(p,q)为：

式中，符号是小数向下取整运算；

步骤1-4)对随机粒子进行融合处理，得到t_k时刻的等效量测z_k；

所述融合处理的过程为：统计每个栅格内随机粒子的数目，设定一个阈值，删除随机粒子个数小于阈值的栅格内的所有随机粒子，用于抑制概率密度函数拖尾区域对应的小概率随机粒子对多基地量测融合性能的影响；对剩下的随机粒子进行加权融合处理，得到k时刻的等效量测z_k的两个分量：

$\left\{\begin{array}{c}{z}_k\left(x\right)=\frac{1}{N^{\′}}\underset{i}{\mathrm{\Σ}}P({\mathrm{\π}}_{\mathrm{ij}},x)N\left({\mathrm{\π}}_{\mathrm{ij}}\right)\\ z_k\left(y\right)=\frac{1}{N^{\′}}\underset{j}{\mathrm{\Σ}}P({\mathrm{\π}}_{\mathrm{ij}},y)N\left({\mathrm{\π}}_{\mathrm{ij}}\right)\end{array}\right.---\left(5\right)$

式中，N(π_ij)表示经小概率采样粒子抑制后参与融合处理的序号为(i,j)的栅格内的随机粒子数目，P(π_ij,x)表示栅格π(i,j)的x坐标，P(π_ij,y)表示栅格π(i,j)的y坐标；N′为经小概率采样粒子抑制后所有栅格内粒子的总数；z_k＝(z_k(x),z_k(y))。

上述技术方案中，所述步骤2)具体包括：

步骤2-1)根据目标运动模型f，用t_k-1时刻的随机粒子预测t_k时刻的随机粒子；

$S_k^{\left(i\right)}=f(S_{k-1}^{\left(i\right)},{\mathrm{\ω}}_{k-1}^{\left(i\right)},w_{k-1})---\left(6\right)$

其中，i＝1…N，N表示随机粒子的数目；f(·)表示目标运动模型，w_k-1表示t_k-1时刻系统过程噪声，为t_k-1时刻的第i个随机粒子；为t_k-1时刻的第i个随机粒子的权值；为t_k时刻的第i个预测粒子；

特别的，在t₀时刻，对探测结果的先验概率分布p(X₀)进行随机采样，生成N个随机粒子并将随机权值设为i＝1…N；

步骤2-2)通过步骤1)的等效量测z_k得到的似然函数进一步得到的归一化权重；

${\mathrm{\ω}}_k^{\left(i\right)}={\mathrm{\ω}}_{k-1}^{\left(i\right)}L\left(z_k\right|S_k^{\left(i\right)})---\left(7\right)$

${\mathrm{\ω}}_k^{\left(i\right)}=\frac{{\mathrm{\ω}}_k^{\left(i\right)}}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^{\left(i\right)}}---\left(8\right)$

步骤2-3)根据所述等效量测z_k及先验概率分布对等效量测进行随机采样，生成M个目标量测采样粒子其中j＝1,2,…,M；

步骤2-4)计算所述M个目标量测采样粒子的权重；

根据公式(6)结合目标运动规律得到t_k时刻预测粒子的位置分布，从而计算目标量测采样粒子的似然函数其中，Z_k|k-1是根据t_k-1时刻随机粒子经运动模型f转移到t_k时刻后得到的目标量测的预测；从而进一步求得目标量测采样粒子的权重为：

${\mathrm{\ω}}_{\mathrm{mk}}^{\left(j\right)}=\frac{L\left(S_{\mathrm{mk}}^{\left(j\right)}\right|Z_{k|k-1})}{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}L\left(S_{\mathrm{mk}}^{\left(j\right)}\right|Z_{k|k-1})}---\left(9\right)$

步骤2-5)将N个预测粒子和M个目标量测采样粒子进行融合，得到混合粒子滤波的目标状态估计X_k：

$X_k=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^{\left(i\right)}{S}_k^{\left(i\right)}+\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ω}}_{\mathrm{mk}}^{\left(j\right)}{S}_{\mathrm{mk}}^{\left(j\right)}}{2}---\left(10\right)$

步骤2-6)对N个预测粒子和M个目标量测采样粒子进行重采样处理得到N个重采样粒子，并将N个重采样粒子替换N个预测粒子。

本发明的优点在于：

1、本发明的方法利用了水下小目标多基地主动探测的混合干扰量测在各声纳节点的稳定性弱于目标观测的特点，通过量测重采样的随机粒子融合，实现了水下小目标的多基地量测融合；

2、本发明的方法针对水下小目标运动模型与实际运动难以较好吻合，粒子滤波中“粒子贫化”造成的目标跟踪性能下降的问题，将两种粒子进行融合处理，保留了传统粒子滤波过程在处理非线性、非高斯问题方面的优势；提高了水下小目标的状态估计精度。

附图说明

图1为本发明的水下小目标多基地融合跟踪方法的流程图；

图2为本发明的方法与传统粒子滤波方法的蛙人状态估计结果的对比图；

图3为本发明的方法与传统粒子滤波方法的蛙人状态估计的均方根误差对比图。

具体实施方式

本发明的方法利用两级粒子滤波实现水下小目标的多基地集中式融合跟踪，主要包括量测融合和状态估计两方面。根据多基地系统各声纳节点得到的目标观测和混响干扰观测稳定性不同，通过多基地量测的预检测随机粒子完成多基地系统下水下小目标的量测融合；为改善水下小目标机动时粒子滤波算法的状态估计性能，通过量测重采样处理得到的随机粒子和传统粒子滤波估计仅由初始化经状态转移的粒子组成的混合粒子集的加权求和，完成多基地系统下水下小目标的状态估计，降低了目标模型与目标实际运动不能准确吻合时粒子偏离目标真实位置对目标状态估计的影响，得到水下小目标更高精度、更稳健的多基地状态估计，改善水下小目标融合跟踪性能。

粒子滤波状态估计方法以贝叶斯估计理论为基础，通过t_k时刻目标状态X_k后验概率分布p(X_k|Y_k)，按照一定的估计规则完成目标状态X_k的估计，其中，Y_k＝{y₀,y₁,…y_k}表示从起始t₀时刻到t_k时刻声纳得到目标量测集合。

粒子滤波估计方法借助蒙特卡罗思想利用满足一定分布规律的带权随机粒子集表征目标状态的后验概率。设从后验概率密度p(x_k|Y_k)中抽取N个独立同分布的随机样本z^(s)，s＝1,…,N，则后验概率分布p(X_k|Y_k)可近似逼近，如式(1)所示：

$p\left(X_k\right|Z_k)\≈\underset{s=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^s\mathrm{\δ}(z_k-z_k^{\left(s\right)})---\left(1\right)$

其中，表示t_k时刻粒子s对应的权重，N表示粒子总数，

$\mathrm{\δ}(z_k-z_k^{\left(s\right)})=\left\{\begin{array}{c}1,z_k=z_k^{\left(s\right)}\\ 0,z_k\≠z_k^{\left(s\right)}\end{array}\right.---\left(2\right)$

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

如图1所示，一种水下小目标多基地融合跟踪方法，所述方法包括：

步骤1)对t_k时刻多基地系统的目标量测集合Z_k进行重采样，得到重采样随机粒子集合通过对所述重采样随机粒子的分布进行量化和融合处理，得到t_k时刻的等效量测z_k；具体包括：

步骤1-1)在t_k时刻，对所有小目标探测声纳的接收数据进行预处理，获取多基地系统的目标量测集合Z_k；

所述预处理包括：带通采样、降采样和时空配准；得到t_k时刻多基地系统的目标量测集合为：

$Z_k=\{z_{\mathrm{1,1}},z_{\mathrm{1,2}},\·\·\·,z_{1,N_1},\·\·\·,z_{M,N_M}\};$

其中，M为多基地系统中声纳节点数目；N_i为第i个声纳节点得到的量测数目；

步骤1-2)根据多基地系统各声纳节点在量测位置处的观测噪声协方差，对目标量测集合Z_k进行重采样，得到重采样随机粒子的集合其中

步骤1-3)对所述重采样随机粒子的分布进行量化；

首先对所述重采样随机粒子所在的平面进行了P×Q栅格剖分；将所述重采样随机粒子所在的平面划分为边长h的栅格，则序号(i,j)的栅格π(i,j)左下角的位置(x_i,y_j)为：

$\left\{\begin{array}{c}{x}_i=x_0+i*h\\ y_j=y_0+j*h\end{array}\right.---\left(3\right)$

式中，(x₀,y₀)为栅格化平面坐标起点；

计算落入每个栅格内随机粒子的数目；若随机粒子的二维坐标为(x,y)，则该随机样本所在栅格π(p,q)的栅格序号(p,q)为：

式中，符号是小数向下取整运算。

步骤1-4)对随机粒子进行融合处理，得到t_k时刻的等效量测z_k；

所述融合处理的过程为：统计每个栅格内随机粒子的数目，设定一个阈值，删除随机粒子个数小于阈值的栅格内的所有随机粒子，用于抑制概率密度函数拖尾区域对应的小概率随机粒子对多基地量测融合性能的影响。对剩下的随机粒子进行加权融合处理，得到k时刻的等效量测z_k的两个分量：

$\left\{\begin{array}{c}{z}_k\left(x\right)=\frac{1}{N^{\′}}\underset{i}{\mathrm{\Σ}}P({\mathrm{\π}}_{\mathrm{ij}},x)N\left({\mathrm{\π}}_{\mathrm{ij}}\right)\\ z_k\left(y\right)=\frac{1}{N^{\′}}\underset{j}{\mathrm{\Σ}}P({\mathrm{\π}}_{\mathrm{ij}},y)N\left({\mathrm{\π}}_{\mathrm{ij}}\right)\end{array}\right.---\left(5\right)$

式中，N(π_ij)表示经小概率采样粒子抑制后参与融合处理的序号为(i,j)的栅格内的随机粒子数目，P(π_ij,x)表示栅格π(i,j)的x坐标，P(π_ij,y)表示栅格π(i,j)的y坐标；N′为经小概率采样粒子抑制后所有栅格内粒子的总数；z_k＝(z_k(x),z_k(y))；

步骤2)通过混合粒子滤波对t_k时刻小目标状态进行估计；并对混合粒子进行重采样；具体包括：

步骤2-1)根据目标运动模型f，用t_k-1时刻的随机粒子预测t_k时刻的随机粒子；

$S_k^{\left(i\right)}=f(S_{k-1}^{\left(i\right)},{\mathrm{\ω}}_{k-1}^{\left(i\right)},w_{k-1})---\left(6\right)$

其中，i＝1…N，N表示随机粒子的数目；f(·)表示目标运动模型，w_k-1表示t_k-1时刻系统过程噪声，为t_k-1时刻的第i个随机粒子；为t_k-1时刻的第i个随机粒子的权值；为t_k时刻的第i个预测粒子；

特别的，在t₀时刻，对探测结果的先验概率分布p(X₀)进行随机采样，生成N个随机粒子并将随机权值设为i＝1…N。

步骤2-2)通过步骤1)的等效量测z_k得到的似然函数进一步得到的归一化权重；

${\mathrm{\ω}}_k^{\left(i\right)}={\mathrm{\ω}}_{k-1}^{\left(i\right)}L\left(z_k\right|S_k^{\left(i\right)})---\left(7\right)$

${\mathrm{\ω}}_k^{\left(i\right)}=\frac{{\mathrm{\ω}}_k^{\left(i\right)}}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^{\left(i\right)}}---\left(8\right)$

步骤2-3)根据所述等效量测z_k及先验概率分布对等效量测进行随机采样，生成M个目标量测采样粒子其中j＝1,2,…,M；

步骤2-4)计算所述M个目标量测采样粒子的权重；

根据公式(6)结合目标运动规律得到t_k时刻预测粒子的位置分布，从而计算目标量测采样粒子的似然函数其中，Z_k|k-1是根据t_k-1时刻随机粒子经运动模型f转移到t_k时刻后得到的目标量测的预测；从而进一步求得目标量测采样粒子的权重为：

${\mathrm{\ω}}_{\mathrm{mk}}^{\left(j\right)}=\frac{L\left(S_{\mathrm{mk}}^{\left(j\right)}\right|Z_{k|k-1})}{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}L\left(S_{\mathrm{mk}}^{\left(j\right)}\right|Z_{k|k-1})}---\left(9\right)$

步骤2-5)将N个预测粒子和M个目标量测采样粒子进行融合，得到混合粒子滤波的目标状态估计X_k：

$X_k=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^{\left(i\right)}{S}_k^{\left(i\right)}+\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ω}}_{\mathrm{mk}}^{\left(j\right)}{S}_{\mathrm{mk}}^{\left(j\right)}}{2}---\left(10\right)$

步骤2-6)对N个预测粒子和M个目标量测采样粒子进行重采样处理得到N个重采样粒子，并将N个重采样粒子替换N个预测粒子；

步骤3)令k＝k+1，转入步骤1)；直至目标跟踪结束。

下面结合具体实例对本发明作进一步详细说明。

多基地探测系统包括发射端节点S1和接收端节点R1、R2，声纳信号采用中心频率70kHz，带宽8kHz，脉冲长度为6ms的线性调频信号LFM信号，对水下蛙人目标进行探测。

为说明本发明的方法在进行水下小目标多基地集中式量测跟踪方面的优势，分别利用传统粒子滤波状态估计方法和本发明的方法对多基地观测数据进行处理，实现蛙人的状态估计。为进一步刻画不同方法的性能，这里将上述估计结果的均方根误差评估不同方法的估计性能。t_k时刻状态估计的均方根误差公式为：

$\mathrm{RMSE}=\sqrt{{({\hat{x}}_k-x_k)}^2+{({\hat{y}}_k-y_k)}^2}---\left(11\right)$

其中，表示t_k时刻的目标状态估计得到的位置坐标，(x_k,y_k)表示t_k时刻目标的真实位置。

如图2和图3所示，通过本发明的方法得到的目标状态估计结果是最优的，该方法在目标状态估计时首先利用多基地系统的目标量测值进行融合，得到等效量测，并通过增加采样粒子，丰富了传统粒子滤波估计随机粒子的多样性，得到多基地系统下水下小目标精度更高、稳健性更好的目标状态估计。

Claims (3)

1.一种水下小目标多基地融合跟踪方法，所述方法包括：

步骤1)对t_k时刻多基地系统的目标量测集合Z_k进行重采样，得到重采样随机粒子集合通过对所述重采样随机粒子的分布进行量化和融合处理，得到t_k时刻的等效量测z_k；

步骤2)通过混合粒子滤波对t_k时刻小目标状态进行估计；所述混合粒子包括预测粒子和目标量测采样粒子；所述预测粒子为用t_k-1时刻的随机粒子预测得到的t_k时刻的随机粒子；所述目标量测采样粒子根据所述等效量测z_k及先验概率分布对等效量测进行随机采样得到；并对混合粒子进行重采样；

步骤3)令k＝k+1，转入步骤1)；直至目标跟踪结束。

2.根据权利要求1所述的水下小目标多基地融合跟踪方法，其特征在于，所述步骤1)具体包括：

步骤1-1)在t_k时刻，对所有小目标探测声纳的接收数据进行预处理，获取多基地系统的目标量测集合Z_k；

所述预处理包括：带通采样、降采样和时空配准；得到t_k时刻多基地系统的目标量测集合为：

$Z_k=\{z_{\mathrm{1,1}},z_{\mathrm{1,2}},...,z_{1,N_1},...,z_{M,N_M}\};$

其中，M为多基地系统中声纳节点数目；N_i为第i个声纳节点得到的量测数目；

步骤1-2)根据多基地系统各声纳节点在量测位置处的观测噪声协方差，对目标量测集合Z_k进行重采样，得到重采样随机粒子的集合其中

步骤1-3)对所述重采样随机粒子的分布进行量化；

首先对所述重采样随机粒子所在的平面进行了P×Q栅格剖分；将所述重采样随机粒子所在的平面划分为边长h的栅格，则序号(i,j)的栅格π(i,j)左下角的位置(x_i,y_j)为：

$\left\{\begin{array}{c}{x}_i=x_0+i*h\\ y_j=y_0+j*h\end{array}\right.---\left(3\right)$

式中，(x₀,y₀)为栅格化平面坐标起点；

计算落入每个栅格内随机粒子的数目；若随机粒子的二维坐标为(x,y)，则该随机样本所在栅格π(p,q)的栅格序号(p,q)为：

式中，符号是小数向下取整运算；

步骤1-4)对随机粒子进行融合处理，得到t_k时刻的等效量测z_k；

所述融合处理的过程为：统计每个栅格内随机粒子的数目，设定一个阈值，删除随机粒子个数小于阈值的栅格内的所有随机粒子，用于抑制概率密度函数拖尾区域对应的小概率随机粒子对多基地量测融合性能的影响；对剩下的随机粒子进行加权融合处理，得到t_k时刻的等效量测z_k的两个分量：

$\left\{\begin{array}{c}{z}_k\left(x\right)=\frac{1}{N^{\′}}\underset{i}{\mathrm{\Σ}}P({\mathrm{\π}}_{\mathrm{ij}},x)N\left({\mathrm{\π}}_{\mathrm{ij}}\right)\\ z_k\left(y\right)=\frac{1}{N^{\′}}\underset{j}{\mathrm{\Σ}}P({\mathrm{\π}}_{\mathrm{ij}},y)N\left({\mathrm{\π}}_{\mathrm{ij}}\right)\end{array}\right.---\left(5\right)$

式中，N(π_ij)表示经小概率采样粒子抑制后参与融合处理的序号为(i,j)的栅格内的随机粒子数目，P(π_ij,x)表示栅格π(i,j)的x坐标，P(π_ij,y)表示栅格π(i,j)的y坐标；N′为经小概率采样粒子抑制后所有栅格内粒子的总数；z_k＝(z_k(x),z_k(y))。

3.根据权利要求1或2所述的水下小目标多基地融合跟踪方法，其特征在于，所述步骤2)具体包括：

步骤2-1)根据目标运动模型f，用t_k-1时刻的随机粒子预测t_k时刻的随机粒子；

$S_k^{\left(i\right)}=f(S_{k-1}^{\left(i\right)},{\mathrm{\ω}}_{k-1}^{\left(i\right)},w_{k-1})---\left(6\right)$

其中，i＝1…N，N表示随机粒子的数目；f(·)表示目标运动模型，w_k-1表示t_k-1时刻系统过程噪声，为t_k-1时刻的第i个随机粒子；为t_k-1时刻的第i个随机粒子的权值；为t_k时刻的第i个预测粒子；

特别的，在t₀时刻，对探测结果的先验概率分布p(X₀)进行随机采样，生成N个随机粒子并将随机权值设为

步骤2-2)通过步骤1)的等效量测z_k得到的似然函数进一步得到的归一化权重；

${\mathrm{\ω}}_k^{\left(i\right)}={\mathrm{\ω}}_{k-1}^{\left(i\right)}L\left(z_k\right|S_k^{\left(i\right)})---\left(7\right)$

${\mathrm{\ω}}_k^{\left(i\right)}=\frac{{\mathrm{\ω}}_k^{\left(i\right)}}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^{\left(i\right)}}---\left(8\right)$

步骤2-3)根据所述等效量测z_k及先验概率分布对等效量测进行随机采样，生成M个目标量测采样粒子其中 $S_{\mathrm{mk}}^{\left(j\right)}~p\left(z_k\right),j=\mathrm{1,2},...,M;$

步骤2-4)计算所述M个目标量测采样粒子的权重；

根据公式(6)结合目标运动规律得到t_k时刻预测粒子的位置分布，从而计算目标量测采样粒子的似然函数其中，Z_k|k-1是根据t_k-1时刻随机粒子经运动模型f转移到t_k时刻后得到的目标量测的预测；从而进一步求得目标量测采样粒子的权重为：

${\mathrm{\ω}}_{\mathrm{mk}}^{\left(j\right)}=\frac{L\left(S_{\mathrm{mk}}^{\left(j\right)}\right|Z_{k|k-1})}{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}L\left(S_{\mathrm{mk}}^{\left(j\right)}\right|Z_{k|k-1})}---\left(9\right)$

步骤2-5)将N个预测粒子和M个目标量测采样粒子进行融合，得到混合粒子滤波的目标状态估计X_k：

$X_k=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^{\left(i\right)}{S}_k^{\left(i\right)}+\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ω}}_{\mathrm{mk}}^{\left(j\right)}{S}_{\mathrm{mk}}^{\left(j\right)}}{2}---\left(10\right)$

步骤2-6)对N个预测粒子和M个目标量测采样粒子进行重采样处理得到N个重采样粒子，并将N个重采样粒子替换N个预测粒子。