一种基于目标方位与线谱的小平台水声被动探测信息融合方法，它属于水下信息融合技术领域。本发明解决了现有方法中线谱和方位的跟踪精度低、滤波过程中易发散难收敛以及易出现目标频率特征关联错误的问题。本发提出了基于多信息矢量关联的卡尔曼滤波算法，采用自适应波门，将目标信息矢量加权融合处理，在频率与方位两个维度对目标信息矢量进行估计，解决了滤波过程中易发散难收敛以及方位历程交叉、临近时易出现错跟或漏跟的问题，提高了目标方位、线谱的跟踪精度，保证了目标频率特征关联的正确性。本发明方法可以应用于水下信息融合。
一种基于目标方位与线谱的小平台水声被动探测信息融合
方法
技术领域
[0001] 本发明属于水下信息融合技术领域，具体涉及一种基于目标方位与线谱的小平台水声被动探测信息融合方法。
背景技术
[0002] 水下小平台被动声呐系统具有性价比高、隐蔽性好、使用灵活与多平台协同探测可提高探测性能等优势，受到了广泛的应用，如潜标、浮标、UUV等。相较于大平台的大孔径阵列，小平台受尺寸限制，其测向稳定性与空间增益的获取受到了很大的影响，这给后续目标定位跟踪带来了较大的困难。水下目标线谱信息较为丰富，目标方位是辐射源较为可靠的参数，故常将目标的线谱、方位等信息进行滤波、关联和融合处理，以获得每个目标的信息矢量估计(包括方位估计、线谱估计、能量估计、频率特征集)。
[0003] 现有方法对基于方位的多目标被动探测研究较多，且基于方位的多目标被动探测研究方法大致可分为两种：
[0004] 1、利用Hough变换或Radon变换(Magli,Enrico,Olmo,et al.On‑Board Selection of Relevant Images:An Application to Linear Feature Recognition.[J].IEEE Transactions on Image Processing,2001.以及XieYonghong,Ji Qiang.Effective line detection with error propagation[C]International Conference on Image Processing.IEEE,2001.)建立图像空间与状态参数之间的映射，根据统计特性在方位历程图中提取出方位历程的几何曲线，该方法已较为成熟。但由于目标的方位历程通常为不规则的几何曲线，采用该方法时空开销过大(郑援,胡成军,李启虎,一种多目标方位历程实时提取方法[J].声学学报,2005,(01):83‑8.)，难以实时处理。
[0005] 2、对方位历程进行峰值提取以获得目标的方位估计(孟彧仟,石超雄,朱宗斌.基于方位历程的多目标自动检测与跟踪[C]2016中国西部声学学术交流会,F,2016)，该方法也被称为最大值法。但该方法存在一定的缺陷：在跟踪时，若强干扰或多个稳定目标同时存在，易导致目标方位跟踪错误，稳定性较差。故学者们对其进行了改进：文献(方尔正,崔凯.基于矢量水听器的一种时间方位历程目标跟踪方法[J].应用声学,2005,(05):311‑6.)不仅提取了方位最大(峰)值，又提取了方位次最大值与第三最大值并采用了最小二乘拟合与加权平均法对数据进行处理，在一定程度上改善了上述缺陷；文献(丁超,高鑫,刘启军.卡尔曼滤波在时间方位历程图处理中的应用[C]2016中国西部声学学术交流会,中国四川成都,F,2016.以及金盛龙,李宇,黄海宁.水下多目标方位的联合检测与跟踪[J].声学学报,2019,44(04):503‑12.)则分别引入了卡尔曼滤波与粒子滤波对提取出的峰值数据进行处理，结果表明可以较准确得估计出多目标方位信息。
[0006] 现有方法对基于方位与线谱的多目标信息融合研究较少，且目前基于目标方位与线谱的小平台被动探测信息融合主要面临以下几个困难：
[0007] 一、受小平台尺寸限制与海洋声场的影响，目标线谱、方位量测误差较大，易出现线谱丢失、中断的情况，线谱和方位的跟踪精度会降低；
[0008] 二、目标的方位量测与状态参数之间呈非线性关系，滤波过程中易发散难收敛；
[0009] 三、多目标方位历程存在交叉、临近的情况下，易出现方位跟踪错误、跟踪中断与目标频率特征关联错误的问题。
发明内容
[0010] 本发明的目的是为了解决现有方法中线谱和方位的跟踪精度低、滤波过程中易发散难收敛以及易出现目标频率特征关联错误的问题，而提出的一种基于目标方位与线谱的小平台水声被动探测信息融合方法。
[0011] 本发明为解决上述技术问题所采取的技术方案是：一种基于目标方位与线谱的小平台水声被动探测信息融合方法，所述方法具体包括以下步骤：
[0012] 步骤一、根据每时刻检测到的多目标线谱频率与频率能量对线谱频率进行跟踪，获得目标线谱的频率信息矢量估计结果；
[0013] 步骤二、对步骤一中获得的频率信息矢量估计结果所对应的方位以及方位能量信息进行跟踪，获得线谱下目标的方位信息矢量估计结果；
[0014] 步骤三、对所有线谱下目标的方位信息矢量估计结果进行直方图统计，并对统计结果进行信息矢量跟踪和目标频率特征关联，获得每个目标的信息矢量估计结果。
[0015] 进一步地，所述步骤一的具体过程为：
[0016] 步骤一一、线谱频率轨迹起始
[0017] 对于k时刻输入的目标线谱频率与频率能量集合{Fd(k),Ad(k),d＝1,2,…,D}，分别将集合内的元素Fd(k)作为线谱频率轨迹的伪起始点；
[0018] 其中，D为检测到的线谱数，Fd(k)为k时刻第d个线谱的频率，Ad(k)为k时刻第d个线谱的频率能量；
[0019] 对于伪起始点Fj(k)，Fj(k)∈{Fd(k),Ad(k),d＝1,2,…,D}，判断在Fj(k)的自适应波门G内是否能找到检测到的k+1时刻频率Fj(k+1)，若能找到检测到的频率Fj(k+1)，则继续执行步骤1，否则继续执行步骤2；
[0020] 步骤1、判断伪起始点Fj(k)是否在已经判定为跟踪或起始线谱的自适应波门G内，即判断Fj(k)是否满足：Fj(k)∈{[Fjj(k)‑G,Fjj(k)+G],jj＝1,2,…,J′}，Fjj(k)代表已经判定为跟踪或起始线谱的第jj个线谱在k时刻的频率，J′代表已经判定为跟踪或起始线谱的个数；
[0021] 若Fj(k)∈{[Fjj(k)‑G,Fjj(k)+G],jj＝1,2,…,J′}，则不将Fj(k)作为k时刻的线谱频率轨迹起始点 ，若 则构建k时 刻目 标状态其中 T为采样时间，并将{Fj(k),Aj(k),Xj(k|
k)}作为k时刻的线谱频率轨迹起始点；
[0022] 步骤2、不将Fj(k)作为k时刻的线谱频率轨迹起始点；
[0023] 同理，遍历每一个伪起始点后，获得k时刻线谱频率轨迹起始点集{Fj(k),Aj(k),Xj(k|k),j＝1,2,…,J}，J代表k时刻线谱频率轨迹起始点的个数；
[0024] 步骤一二、线谱频率轨迹跟踪
[0025] 对于起始点集中的第j个线谱频率，目标状态的一步预测和量测的一步预测分别为：
[0026] Xj(k+1|k)＝F(k)Xj(k|k)
[0027] Zj(k+1|k)＝H(k+1)Xj(k+1|k)
[0028] 其中，F(k)代表k时刻的状态转移矩阵，Xj(k+1|k)代表预测出的k+1时刻目标状态，H(k+1)为k+1时刻的量测矩阵，Zj(k+1|k)代表预测出的k+1时刻线谱频率值；
[0029] 判断在Zj(k+1|k)的自适应波门G内是否检测到线谱；若在Zj(k+1|k)的自适应波门G内检测到线谱，则继续执行步骤a，否则继续执行步骤b；
[0030] 步骤a、将检测到的线谱频率与线谱频率对应的能量组成的备用量测集记为{Fs(k+1),As(k+1)，s＝1,2,…,S}，S表示在Zj(k+1|k)的自适应波门G内检测到的线谱数Fs(k+1)代表检测到的第s个线谱在k+1时刻的频率，As(k+1)代表检测到的第s个线谱在k+1时刻的频率能量；
[0031] 对备用量测集进行处理，得到第j个线谱在k+1时刻的频率量测Zj(k+1)以及第j个线谱在k+1时刻的频率能量Aj(k+1)：
[0032]
[0033]
[0034] 其中，var(·)代表求方差；
[0035] 计算协方差矩阵的一步预测P(k+1|k)：
[0036] P(k+1|k)＝F(k)P(k|k)FT(k)+Q(k)
[0037] 其中，上角标T代表转置，P(k|k)为k时刻的协方差矩阵，Q(k)为k时刻的状态噪声协方差矩阵；
[0038] 计算k+1时刻的新息协方差矩阵S(k+1)：
[0039] S(k+1)＝H(k+1)P(k+1|k)HT(k+1)+R(k+1)
[0040] 其中，R(k+1)是k+1时刻的量测噪声协方差矩阵；
[0041] 计算k+1时刻的增益K(k+1)：
[0042] K(k+1)＝P(k+1|k)HT(k+1)S‑1(k+1)
[0043] 其中，上角标‑1代表矩阵的逆；
[0044] Xj(k+1|k+1)＝Xj(k+1|k)+K(k+1)(Zj(k+1|k)‑Zj(k+1))
[0045] 其中，Xj(k+1|k+1)代表k+1时刻的目标状态；
[0046] P(k+1|k+1)＝P(k+1|k)‑K(k+1)S(k+1)KT(k+1)
[0047] 其中，P(k+1|k+1)为k+1时刻的协方差矩阵；
[0048] 步骤b、将Zj(k+1|k)作为k+1时刻的目标状态Xj(k+1|k+1)，将k+1时刻的频率能量Aj(k+1)取为Aj(k)；
[0049] 在第j个线谱的频率轨迹跟踪过程中，若在连续的Tp个跟踪时刻均执行步骤b，则终止第j个线谱的频率轨迹跟踪，否则继续对第j个线谱的频率轨迹进行跟踪；
[0050] 同理，对起始点集中的每个线谱频率轨迹进行跟踪；
[0051] 步骤一三、线谱频率轨迹更新
[0052] 对于k+1时刻输入的目标线谱频率与频率能量集合，采用步骤一一的方法对k+1时刻输入的目标线谱频率与频率能量集合进行处理，若k+1时刻存在线谱频率轨迹起始点，则采用步骤一二的方法对线谱轨迹起始点进行轨迹跟踪；
[0053] 同理，对后面其它时刻输入的目标线谱频率与频率能量集合进行处理。
[0054] 进一步地，所述自适应波门G的计算方法为：
[0055]
[0056]
[0057] 式中，自适应波门G的单位为°， 为输入信噪比，为目标方位俯仰角，远距离时取90°，B为带宽，T0为积分时间，对于单线谱而言，取BT0＝1。
[0058] 进一步地，所述步骤二的具体过程为：
[0059] 步骤二一、方位轨迹起始
[0060] 将k时刻某根线谱e的频率信息矢量表示为{Fe(k),Ae(k)}，e＝1,2,…,E，E表示在k时刻已经跟踪的线谱的个数，将线谱e在k时刻的方位和方位能量集合记为代表线谱e在k时刻的第c个方位， 代表线谱e在k时
刻的第c个方位能量；
[0061] 分别以 作为方位轨迹伪起始点，Ce代表方位轨迹伪起始点总数，对于伪起始点 判断在 的自适应波门G内是否能
找到检测到的k+1时刻方位 若能找到检测到的方位 则继续执行步骤1)，
否则继续执行步骤2)；
[0062] 步骤1)、判断伪起始点 是否在已经判定为跟踪或起始方位的自适应波门G内，即判断 是否满足： θnn(k)代表已经判定为跟踪或起始的第nn个方位，N代表已经判定为跟踪或起始方位的个数；
[0063] 若 则不将 作为k时刻的方位轨迹起始点，若 则构建k时刻目标状态
其中 T为采样时间，并将
作为k时刻的方位轨迹起始点；
[0064] 步骤2)、不将 作为k时刻的方位轨迹起始点；
[0065] 同理 ，遍历每 一个伪 起始点 后，获得 k时刻 方位轨迹 起始点 集N代表k时刻方位轨迹起始点的个数；
[0066] 步骤二二、方位轨迹跟踪
[0067] 对于起始点集中的第n个方位，目标状态的一步预测和量测的一步预测分别为：
[0068]
[0069]
[0070] 其中， 代表预测出的k+1时刻目标状态， 代表预测出的k+1时刻方位值；
[0071] 判断在 的自适应波门G内是否检测到方位，若在 的自适应波门G内检测到方位，则继续执行步骤a)，否则继续执行步骤b)；
[0072] 步骤a)、将检测到的方位与方位对应的方位能量组成的备用量测集记为S′表示在 的自适应波门G内检测到的方位
数， 代表检测到的第s′个方位， 代表检测到的第s′个方位能量；
[0073] 对备用量测集进行处理，得到k+1时刻的方位量测 以及k+1时刻的方位能量
[0074]
[0075]
[0076] 再采用步骤a中的方法计算k+1时刻目标状态
[0077] 步骤b)、将 作为k+1时刻的目标状态 将k+1时刻的方位能量 取为
[0078] 在第n个方位的轨迹跟踪过程中，若在连续的Tp个跟踪时刻均执行步骤b)，则终止第n个方位的轨迹跟踪，否则继续对第n个方位的轨迹进行跟踪；
[0079] 同理，对起始点集 中的其它各个方位以及k时刻的其它线谱的各个方位进行轨迹跟踪。
[0080] 更进一步地，所述步骤三的具体过程为：
[0081] 步骤三一、对所有线谱下目标的方位信息矢量估计结果进行直方图统计：
[0082]
[0083] 其中，{θ(k),A(k)}代表k时刻各线谱的方位、方位能量以及频率的融合结果；
[0084] 寻找{θ(k),A(k)}中检测到的全部方位轨迹起始点，并将检测到的全部方位轨迹起始点组成的集合表示为{θp(k),Ap(k),p＝1,2,…,P}，P为检测到的全部方位轨迹起始点的总数，即目标的个数，若θp(k)∈[θn(k)‑G,θn(k)+G]，θn(k)来自于步骤二跟踪获得的k时刻方位的集合，则将方位θn(k)所在的线谱对应的频率赋予{θp(k),Ap(k)}作为{θp(k),Ap(k)}的频率特征Fp(k)，则k时刻目标信息矢量起始集为{θp(k),Ap(k),Fp(k),p＝1,2,…,P}；
[0085] 步骤三二、采用步骤一二的方法分别对k时刻目标信息矢量起始集中的每个起始点进行跟踪，并采用步骤一三的方法对目标信息矢量起始点进行更新；
[0086] 且在目标信息矢量跟踪过程中，对于k+1时刻的方位θp(k+1)，若存在θp(k+1)∈[θn′(k+1)‑G,θn′(k+1)+G]，则将方位θn′(k+1)所在的线谱在k+1时刻的频率赋予{θp(k+1),Ap(k+1)}，即将方位θn′(k+1)所在的线谱在k+1时刻的频率与Fp(k)取并，将取并结果作为{θp(k+1),Ap(k+1)}的频率特征；否则，不存在θp(k+1)∈[θn′(k+1)‑G,θn′(k+1)+G]，则继续将Fp(k)作为{θp(k+1),Ap(k+1)}的频率特征；
[0087] 步骤三三、同理，直至获得所有时刻目标的信息矢量估计结果。
[0088] 本发明的有益效果是：
[0089] 本发明引入了目标线谱信息，联合方位信息进行多目标被动探测与信息融合，通过引入自适应波门，将波门内所有检测到的方位与能量提出，通过赋予目标频率特征，在不同维度下利用基于多信息矢量关联的卡尔曼滤波算法对多目标方位进行跟踪，解决了滤波过程中易发散难收敛以及方位历程交叉、临近时易出现错跟或漏跟的问题，提高了目标方位、线谱的跟踪精度，保证了目标频率特征关联的正确性。
附图说明
[0090] 图1为本发明仿真实验涉及的目标在x‑y平面坐标下的运动轨迹图；
[0091] 本发明中的x‑y平面是指平面直角坐标系的x‑y平面；
[0092] 图2a为本发明仿真实验涉及的频率历程图；
[0093] 图2b为本发明仿真实验涉及的方位历程图；
[0094] 图3为本发明仿真实验涉及的线谱跟踪结果图；
[0095] 图4为本发明仿真实验涉及的方位历程跟踪结果图；
[0096] 图5a为本发明仿真实验涉及的目标1频率特征图；
[0097] 图5b为本发明仿真实验涉及的目标2频率特征图；
[0098] 图5c为本发明仿真实验涉及的目标3频率特征图；
[0099] 图5d为本发明仿真实验涉及的目标4频率特征图。
具体实施方式
[0100] 本发明在频率与方位两个维度对目标信息矢量进行估计，避免多目标方位历程交叉、临近时易出现跟踪错误的情况。相较于只给出方位历程跟踪结果的现有方法，本发明方法能够输出每个目标的信息矢量估计(包括方位估计、线谱估计、能量估计、频率特征集)。而且，本发明定义了基于多信息矢量关联的卡尔曼滤波算法，将目标信息矢量(方位、线谱、能量)加权融合处理，相较于传统的最大值法，采用了自适应波门，在方位量测误差较大时，不会输入误差过大的量测信息，减少了方位估计轨迹发生畸变的情况。
[0101] 下面将本发明方法采取的主要步骤介绍如下：
[0102] 步骤一：输入每时刻检测到的多目标线谱与能量{Fd(k),Ad(k),d＝1,2,…,D}，其中F为线谱频率，A为对应能量，D为检测到的线谱数。对线谱{Fd,d＝1,2,…,D}进行跟踪(包括线谱轨迹起始、线谱轨迹跟踪与线谱轨迹更新)。输出目标的线谱信息矢量估计结果(线谱估计、能量估计等)。
[0103] 1.线谱轨迹起始：针对k时刻输入的目标信息矢量{Fd(k),Ad(k),d＝1,2,…,D}，将该集合内的元素作为线谱轨迹伪起始点，D代表线谱轨迹伪起始点总数。
[0104] 以伪起始点Fj(k),Fj(k)∈{Fd(k),Ad(k),d＝1,2,…,D}为例，给出具体步骤：
[0105] 判断k+1时刻在其自适应波门G内是否存在量测{Fj(k+1),Aj(k+1)}，即Fj(k+1)∈{[Fj(k)‑G,Fj(k)+G]}。若存在量测，进入步骤A，否则进入步骤B。
[0106] 步骤A：构建k时刻目标状态 其中 T为采样时间。输出k时刻线谱轨迹起始点集{Fj(k),Aj(k),Xj(k|k),j＝1,2,…,J}。
[0107] 步骤B：不进行线谱轨迹起始点并删除该伪起始点Fj(k)，判断下一个伪起始点Fj+1(k)。
[0108] 以此逻辑遍历所有的伪起始点，即可得到线谱轨迹起始点集{Fj(k),Aj(k),Xj(k|k),j＝1,2,…,J}。
[0109] 若找到的线谱轨迹起始点Fj(k)在已经判定为跟踪或起始的线谱Fjj自适应波门G内，即Fj(k)∈{[Fjj(k)‑G,Fjj(k)+G],jj＝1,2,…,J}，则不进行轨迹起始。
[0110] 波门G由自适应波门提供，下面给出自适应波门的定义：按照和信噪比有关的方位估计克拉美罗下界，首先估计局域信噪比，将方位估计克拉美罗下界的标准差乘以3倍，作为自适应波门。
[0111]
[0112]
[0113] 式中，波门G单位为°， 为输入信噪比， 为目标方位俯仰角，远距离时取90°。对于单线谱而言，取BT0＝1，B为带宽，T0为积分时间。
[0114] 2.线谱轨迹跟踪：输入k时刻线谱轨迹起始点集{Fj(k),Aj(k),Xj(k|k),j＝1,2,…,J}。以某线谱Fj为例，构建k时刻目标状态方程以及量测方程：
[0115] X(k+1|k)＝F(k)X(k|k)+Γ(k)v(k) (3)
[0116] Z(k+1|k)＝H(k+1)X(k+1|k)+w(k) (4)
[0117] 其中，X表示目标状态量，包含频率值以及频率变化率。Z表示量测量，包含频率量测信息。v(k),w(k)分别为过程噪声向量与量测噪声向量。
[0118] 利用卡尔曼滤波的状态一步预测与量测一步预测公式：
[0119] X(k+1|k)＝F(k)X(k|k) (5)
[0120] Z(k+1|k)＝H(k+1)X(k+1|k) (6)
[0121] 预测出k+1时刻的频率值Z(k+1|k)。若在预测值Z(k+1|k)的自适应波门G内存在量测，则进入步骤A，否则进入步骤B。
[0122] 步骤A：找出预测值Z(k+1|k)自适应波门G内检测到的线谱，记录频率与能量于备用量测集{Zs(k+1),As(k+1)，s＝1,2,…,S}中，S表示检测到的线谱数。利用公式(7)对备用量测集进行处理，即可输出k+1时刻的量测Z(k+1)与能量A(k+1)。
[0123]
[0124]
[0125] 进而利用多信息矢量关联的卡尔曼滤波算法公式(5‑10)，即可完成对k+1时刻频率值X(k+1|k+1)的估计。
[0126] P(k+1|k)＝F(k)P(k|k)FT(k)+Q(k) (8)
[0127] S(k+1)＝H(k+1)P(k+1|k)HT(k+1)+R(k+1) (9)
[0128] K(k+1)＝P(k+1|k)HT(k+1)S‑1(k+1) (10)
[0129] X(k+1|k+1)＝X(k+1|k)+K(k+1)(Z(k+1|k)‑Z(k+1)) (11)
[0130] P(k+1|k+1)＝P(k+1|k)‑K(k+1)S(k+1)KT(k+1) (12)
[0131] 步骤B：若在预测值Z(k+1|k)的自适应波门G内不存在量测，则将一步预测值Z(k+1|k)作为k+1时刻的滤波值，即X(k+1|k+1)＝Z(k+1|k)，k+1时刻的能量取A(k+1)＝A(k)。记录一步预测的步数，若连续预测超过Tp(Tp与信噪比有关)次，则终止该线谱轨迹的跟踪；否则，继续利用公式(5‑10)完成对线谱轨迹的跟踪。
[0132] 3.线谱轨迹更新：在进行线谱轨迹跟踪的同时，要进行线谱轨迹更新。利用步骤一中的线谱轨迹起始逻辑判断出k+1时刻的线谱轨迹起始点集合{Fm(k+1),Am(k+1),m＝1,2,…,M}，M代表k+1时刻的线谱轨迹起始点数目，若线谱轨迹起始点Fm(k+1)不在已经判定为跟踪或起始的线谱Fs自适应波门G内，即
S代表已经判定为跟踪或起始的线谱轨
迹数目，则更新一个新的线谱轨迹起始点；否则不对该线谱轨迹进行更新。
[0133] 以此逻辑遍历所有线谱轨迹起始点即可输出各个时刻目标的线谱信息矢量估计结果{Fs(k),As(k),s＝1,2,…,S}。
[0134] 步骤二：输入步骤一输出的k时刻线谱信息矢量估计结果{Fs(k),As(k),s＝1,2,…,S}与线谱对应的方位、能量信息 对所
有单线谱{Fs(k),s＝1,2,…,S}下的方位进行跟踪(包括方位轨迹起始、方位轨迹跟踪)。输出所有线谱下目标的方位信息矢量估计结果(方位估计、能量估计等)。
[0135] 1 .方位轨迹起始：针对某根线谱F在k时刻对应的目标信息矢量将跟踪出的方位与能量记为集合{θc(k),Ac
(k),c＝1,2,…,C}，将该集合内的元素作为方位轨迹伪起始点，C代表方位轨迹伪起始点总数。
[0136] 以伪起始点θn(k),θn(k)∈{θc(k),Ac(k),c＝1,2,…,C}为例，给出具体步骤：
[0137] 判断k+1时刻在其自适应波门G内是否存在量测{θn(k+1),An(k+1)}，即θn(k+1)∈{[θn(k)‑G,θn(k)+G]}。若存在量测，进入步骤A，否则进入步骤B。
[0138] 步骤A：构建k时刻目标状态 其中 T为采样时间。输出k时刻方位轨迹起始点集{θn(k),An(k),Xn(k|k),n＝1,2,…,N}。
[0139] 步骤B：不进行方位轨迹起始点并删除该伪起始点θn(k)，判断下一个伪起始点θn+1(k)。
[0140] 以此逻辑遍历所有的伪起始点，即可得到方位轨迹起始点集{θn(k),An(k),Xn(k|k),n＝1,2,…,N}。
[0141] 若找到的某个方位起始点θn(k)在已经判定为跟踪或起始的方位θnn(k)自适应波门G内，即θn(k)∈{[θnn(k)‑G,θnn(k)+G],nn＝1,2,…,N}，则不进行方位轨迹起始。
[0142] 2.方位轨迹跟踪：输入k时刻线谱F对应的方位轨迹起始点{θn(k),An(k),n＝1,2,…,N}，构建k时刻目标状态 利用公式(5‑6)预测出k+1时刻的方位值Z(k
+1|k)，此时量测Z仅包含方位信息。再利用步骤一中线谱轨迹跟踪的步骤A与步骤B，即可完成对方位轨迹的跟踪。
[0143] 以此逻辑将所有线谱对应的方位跟踪完毕后。即可初步得到k时刻目标的信息矢量估计结果
[0144] 步骤三：输入步骤二输出的k时刻所有线谱下目标的方位信息矢量估计结果对其进行直方图统计并对统计结果进行信息矢量跟踪与目标频率特征关联(包括目标信息矢量起始、目标信息矢量跟踪与目标信息矢量更新)。最后输出每个目标的信息矢量结果(包括方位估计、线谱估计、能量估计、频率特征集)。
[0145] 1.目标信息矢量起始：步骤二输出了所有线谱下目标的方位信息矢量估计结果，现要形成每个目标对应的方位信息矢量。故对步骤二中输出的结果融合处理，利用公式(13)对所有线谱k时刻的估计结果进行直方图统计，形成新的目标信息矢量{θ(k),A(k)}。
[0146]
[0147] 将目标的信息矢量估计{θF(k),AF(k),F(k)}转换为新的目标信息矢量{θ(k),A(k)}后，寻找k时刻集合{θ(k),A(k)}中检测到的方位点，记录方位与能量于集合{θp(k),Ap(k),p＝1,2,…,P}中，P代表方位轨迹起始点总数。判断集合中的元素θp(k)在哪根线谱对应方位的自适应波门G内，若θp(k)∈[θn(k)‑G,θn(k)+G]。其中θn(k)∈{θF(k),AF(k),F(k)}，则将频率特征信息F赋予给该目标信息矢量，即{Fp(k)}＝{F(k)}，其中{Fp(k)}为目标信息矢量中的频率特征集。
[0148] 至此即可输出k时刻目标信息矢量起始集{θp(k),Ap(k),Fp(k),p＝1,2,…,P}。
[0149] 2.目标信息矢量跟踪：输入k时刻目标信息矢量起始集{θp(k),Ap(k),Fp(k),p＝1,2,…,P}，利用步骤二中的方位轨迹起始逻辑，构建k时刻目标状态 再利
用公式(5‑6)预测出k+1时刻的方位值Z(k+1|k)。在该目标频率特征Fp对应的k+1时刻方位跟踪结果中，寻找在Z(k+1|k)自适应波门G内的方位值与其对应的能量，构成备选量测集再将备选量测集输入给公式(7)，
可以计算出k+1时刻的量测Z(k+1)与能量A(k+1)，利用步骤一中线谱轨迹跟踪的步骤A与步骤B，即可完成对目标信息矢量的跟踪。
[0150] 3.目标信息矢量更新：在进行目标信息矢量跟踪的同时，要进行信息矢量更新。对于k+1时刻的信息矢量{θ(k+1),A(k+1)}，利用步骤三中的目标矢量信息起始逻辑判断出k+1时刻的信息矢量起始点集合{θpp(k+1),App(k+1),pp＝1,2,…,PP}，若新的方位轨迹起始点θpp(k+1)不在已经判定为跟踪或起始的方位θp(k+1)自适应波门G内，即
则更新一个新的方位轨迹起始点；否
则不对该方位轨迹进行更新。
[0151] 跟踪时，判断k+1时刻估计的方位θp(k+1)在哪根线谱对应方位的自适应波门G内，即若θp(k+1)∈{θFo(k+1),AFo(k+1),Fo(k+1)}，则将频率特征信息Fo也赋予给该目标，目标信息矢量中的频率特征集更新为{Fp(k+1)}＝{Fp(k)}∪{Fo(k+1)}。
[0152] 至此即可得到k+1时刻目标的信息矢量{θp(k+1),Ap(k+1),Fp(k+1),p＝1,2,…,P}；
[0153] 以此逻辑跟踪完毕后。即可输出所有时刻目标的信息矢量估计结果{θp(k),Ap(k),Fp(k),p＝1,2,…,P,k＝1,2,…,T}，其中θ为目标的方位信息，A为目标的能量信息，F为目标的频率特征集，P为目标数，T为总采样时间。
[0154] 仿真实验
[0155] 仿真条件：假设共有4个做匀速直线运动的目标，每个目标的初始位置、初始速度、辐射的中心频率与声源级由表1给出，4个目标在x‑y坐标下的运动轨迹如图1所示。静止观测站坐标为(‑100，0)m。总共跟踪150帧数据，每帧间隔1s。环境噪声为75dB。方位估计的最小均方误差由与信噪比有关的方位估计克拉美罗下界给出。频率量测精度为0.01Hz，水下声速为1500m/s。
[0156] 表1各目标初始位置、初始速度、辐射中心频率与声源级
[0157]
[0158] 仿真总结：
[0159] 如图2a以及如图2b所示的频率历程图和方位历程图，本发明对四个匀速运动目标进行仿真并给出仿真结果，得到如图3所示的多目标线谱跟踪结果以及如图4所示的方位跟踪结果，可以看出在目标方位历程存在交叉时，没有出现错跟漏跟的现象。图5a至图5d给出了各目标的频率特征，经验证没有出现目标频率特征赋予错误的情况,该仿真结果验证了本发明的可行性。
[0160] 本发明的上述算例仅为详细地说明本发明的计算模型和计算流程，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。
