多阵列水下多目标跟踪的分布式算法研究

徐 文 吴雨桑 张 婷

（1.浙江大学信息与电子工程学院，浙江杭州 310027；2.浙江大学海洋学院，浙江省海洋观测-成像试验区重点实验室，浙江舟山 316021）

1 引言

水下多目标跟踪算法可实时获取多个目标的状态信息，如位置、速度和加速度等，从而分辨各个目标的身份和意图，在军事和民用领域都具有重要意义［1］。目前，水下多目标跟踪算法主要分为基于数据关联的多目标跟踪算法和基于随机有限集（Random Finite Set，RFS）的多目标跟踪算法两类。前者的核心是数据关联，通常要求监测区域中目标数目已知且不随时间变化。当监测区域中目标数目较多时，该类算法计算复杂度较高。基于随机有限集的多目标跟踪算法将目标的状态和量测建模为随机有限集，将多目标跟踪问题建模为集值滤波问题，避免了复杂的数据关联过程，大大降低了算法复杂度，并适用于目标和干扰数目未知且时变的复杂水下环境。2014年，Melo等人将势概率假设密度（Cardinalized Probability Hypothesis Density，CPHD）滤波算法应用于多个自主式水下航行器（Autonomous Underwater Vehicle，AUV）的跟踪问题［2］。2015 年，Saucan 等人将水听器阵列接收的阵元域数据作为量测，改进了叠加近似概率假设密度（Probability Hypothesis Density，PHD）滤波算法，并应用于水下扩展目标的跟踪问题［3］。2016 年，Saucan 等人又将叠加近似CPHD 滤波算法应用于相似场景［4］。由于PHD 和CPHD 滤波算法在进行多目标状态提取时存在聚类过程，导致目标数目和目标状态的估计出现偏差，因此，近年来基于多伯努利随机有限集的水下多目标跟踪算法逐渐得到应用。2017 年，吴秋韵等人将多假设跟踪（Multiple Hypothesis Tracking，MHT）滤波算法和标签多伯努利（Labelled Multi-Bernoulli，LMB）滤波算法应用于水面强目标干扰下的水下单目标跟踪问题［5-6］。2020年，郑策等人将水听器阵列接收的阵元域数据作为量测，改进了多伯努利（Multi-Bernoulli，MB）滤波算法，并应用于水下多目标的波达方向（Direction of Arrival，DOA）跟踪问题［7］。同年，吴孙勇等人将MB滤波算法应用于水下脉冲噪声环境单矢量水听器（Acoustic Vector Sensor，AVS）多目标DOA 跟踪问题［8］。

上述成果都基于单水听器或单水听器阵列等单节点实现。由于声波能量在水下随着距离的增加而逐渐衰减，单节点的接收信噪比逐渐降低。对于较为安静的目标，单节点探测距离受限，只能跟踪一定范围内的运动目标，而距离节点较远的目标则难以准确定位与跟踪。分布式声学传感器网络在监测区域离散布放多个节点，每个节点独立工作，对各自的监测区域进行目标定位与跟踪。不同节点之间通过通信链路交换数据，利用信息融合获得更高的目标定位与跟踪性能［9］。相较于单水听器，使用具有一定孔径的水听器阵列作为分布式声学传感器网络的节点，可以获取相干信号处理增益，使得目标探测距离更远，定位分辨率更高。2018 年，孙文等人将多传感器多伯努利（Multi-Sensor Multi-Bernoulli，MS-MB）滤波算法应用于以到达时间（Time of Arrival，TOA）为量测的分布式多目标跟踪场景［10］，跟踪性能受各节点的序贯更新顺序影响较大。同年，陈中悦等人将伯努利（Bernoulli）滤波算法应用于分布式水听器阵列下的单目标声源跟踪［11］。2020年，郑策将δ-广义标签伯努利（δ-Generalized Labelled Multi-Bernoulli，δ-GLMB）滤波算法应用于多声呐阵列对多目标进行被动跟踪［7］。同年，王冠群提出了一种改进联合多站多伯努利（Improved Multi-Sensor Multi-Bernoulli，IMSMB）滤波算法［12］，得到了较为稳定的目标航迹。上述成果中的分布式声学传感器网络使用的都是集中式融合结构，需要设置一个融合中心处理来自所有节点的信息，并且节点间交换的是存在大量噪声干扰的量测信息，算法的计算复杂度和通信负担较高。目前，在水下多目标跟踪领域应用随机有限集理论的研究中，关于分布式融合结构以及如何选择有效量测集的研究成果还非常少［1］。

本文提出一种基于匹配场定位量测的分布式水下多目标跟踪方法。首先在本地节点进行阵列匹配场定位，在生成的模糊函数上选取多个较高峰坐标作为量测，然后使用势平衡多伯努利（Cardinality Balanced Multi-Bernoulli，CBMB）滤波算法滤除此量测集中的噪声干扰，得到本地节点的多目标后验概率密度，以降低常规匹配场定位中低信噪比导致的定位误差；最后，将各个本地节点与其邻近节点的多目标后验概率密度通过广义协方差交集（Generalized Covariance Intersection，GCI）融合法则进行序贯融合，以充分利用不同节点的量测信息，提高定位与跟踪精度。由于通过数据链路交换的是经滤波的量测集后验概率密度，而不是量测集本身，减少了噪声干扰带来的通信负担和滤波计算量。相较于单节点处理方法，本文提出的方法可以提升多目标联合定位与跟踪性能，且便于在分布式声学传感器网络中实现。

2 问题建模

2.1 目标状态模型

考虑水下固定深度平面运动目标建模，k时刻的目标状态向量为xk=[xk，vx，k，yk，vy，k]T，其中xk和yk对应水平面上x和y两个维度的目标位置，vx，k和vy，k对应两个维度的目标速度。水下目标可使用匀速直线模型建模目标运动状态［11］。目标在k时刻与k+1时刻之间的状态转移方程可以描述为：

2.2 量测模型

k时刻量测向量为zk=[zx，k，zy，k]T，对应匹配场定位模糊函数上高于设定阈值的峰的坐标，zx，k和zy，k分别为x和y位置分量。量测向量与目标状态向量之间的关系建模为：

3 算法原理

3.1 GCI融合法则

常规匹配场定位［13］通常取模糊函数与目标数目相同的峰对应的坐标作为定位输出［5］。实际应用中，当存在环境失配、信噪比低等问题时，定位结果与目标声源的实际位置偏差较大，需要选用多个定位输出结果生成量测集。k时刻阈值设置为，将归一化模糊函数中所有高于的峰对应的坐标选入量测集，数目记为ηk，量测集表示为：

使用多目标跟踪算法滤除量测集中的噪声干扰，得到目标的连续轨迹。

Vo 等人证明MB 滤波算法在跟踪过程中会出现目标数目估计值大于真实值的情况［14］，即产生明显的目标数目（即势）估计偏差。CBMB 滤波算法在MB 滤波算法的基础上加入了一种势修正策略，实现了滤波过程中目标数目的无偏估计，从而更加准确地跟踪目标［14］。

假设分布式水声传感器网络包含n个水听器阵列节点，节点集合记为S={1，2，…，n}。在时刻k，本地节点s∈S使用CBMB 滤波算法得到后验概率密度。根据GCI融合法则［15］，节点s将与邻近节点的后验概率密度进行序贯融合，得到融合后验概率密度。使用CBMB滤波算法处理融合后验概率密度，保留存在概率大于设定阈值的伯努利分量，并使用此分量的概率密度计算状态均值。由本地节点s得到当前时刻目标状态信息的最终估计结果，从而提高跟踪精度。假设节点s的邻近节点集合为Ns⊆S，节点s与其邻近节点的并集为Ss={s}∪Ns，则融合后验概率密度表示为：

其中融合权重ωi满足0 ≤ωi≤1 且，使用Metropolis 权重计算［16］。其中，节点i与节点j(i，j∈S)间的权重为：

di表示网络拓扑结构中节点i的度，即与节点i相连的边的数目。

式（4）的分子中存在复杂的积分运算，实际应用中通常采用序贯蒙特卡洛（Sequential Monte Carlo，SMC）方法进行近似，关于CBMB 滤波器GCI融合的SMC具体实现方法可参见文献［17］。

3.2 集中式、分布式融合方法的性能对比分析

本小节从通信负担、跟踪精度和计算效率三方面比较分布式融合MFP-GCI-CBMB 方法与集中式融合MFP-CEN-CBMB（Matched Field Processing Centralized Cardinality Balanced Multi-Bernoulli）方法的性能。

在集中式融合结构中，其他节点将处理后的数据通过通信链路传输到中心节点，传输数据有两种类型：1.直接传输MFP 处理得到的模糊函数，记为MFP-CEN-CBMB-1 方法；2.传输模糊函数中的所有峰值，记为MFP-CEN-CBMB-2 方法；在分布式融合结构中，存在本地节点与其邻近节点两类节点，本地节点将其邻近节点传输而来的后验概率密度进行融合。两类融合方法的具体实现信息见表1。

1）通信负担

对于集中式融合，MFP-CEN-CBMB-1 方法需要传输模糊函数中每个离散网格点的模糊度值，因此每帧单条通信链路需传输Nmesh个浮点数，其中Nmesh为监测区域的网格点数；MFP-CEN-CBMB-2 方法需要传输模糊函数中的所有峰的坐标和该坐标对应的模糊度值，因此每帧单条通信链路需要传输4×个浮点数，其中为节点i在第k帧中产生的峰的数目。

对于分布式融合方法MFP-GCI-CBMB，在第k帧，节点i表征多目标后验概率密度的参数集为，其中为节点i在第k帧的伯努利分量总数，满足，实现中要求跟踪目标数目≤伯努利分量最大数目Mmax为伯努利分量m单目标后验概率密度的近似粒子总数，满足为伯努利分量m的存在概率，为一标量；ω(m，l)为伯努利分量m第l号粒子的权重，为一标量；x(m，l)为伯努利分量m第l号粒子的状态向量，维度为6×1。节点i需要发送此参数集中每个伯努利分量的存在概率及权重较大的粒子的权重及状态向量。选取粒子时，先将它们按照权重从大到小排序选取，直到它们的权重和≥ω，此时所选取的粒子数为。节点i将参数集传输至所有邻近节点j∈Ni，则节点i第k帧发送的数据量为：

此外，节点i也需要接收来自邻近节点j∈Ni的表征多目标后验概率密度的参数集，j∈Ni。网络中的总通信次数为拓扑结构中所有节点的连通度之和，即。故网络第k帧需传输的数据总量Dk满足：

2）跟踪精度

本文使用最优子模式分配（Optimal Subpattern Assignment，OSPA）距离［18］表征同一网络中的单节点、集中式和分布式融合方法的跟踪精度。

OSPA 距离可以同时衡量复杂多目标场景下目标数目的估计精度和不同目标状态的估计精度。假设某一帧中各个目标的真实状态集为X={x1，x2，…，xNt}，估计状态集为其中Nt和Ne分别代表当前帧的真实目标数目和估计目标数目，p为距离敏感性参数，c为势惩罚参数，OSPA距离可表示为：

3）计算负担

集中式融合中，中心节点需要承担大部分的计算任务，而其他节点只需承担MFP 步骤的计算。MFP-CEN-CBMB-1 方法通过数据链路传输模糊函数，并由中心节点计算峰值点坐标，网络中除中心节点外，其他节点计算负担较轻；而MFP-CENCBMB-2 方法需要其他节点计算模糊函数中的峰值点，因此计算负担较方法1大。分布式融合方法中，本地节点仅需要融合其邻近节点的后验概率密度，每个本地节点的计算负担与其所需融合的邻近节点数目成正比，不同本地节点的计算负担相当。本文4.3节将用平均单帧数据处理时间来衡量同一网络中的集中式、分布式融合方法的计算负担。

4 计算机仿真

4.1 定位仿真实验

监测区域范围设置为［0，5］km×［0，5］km×［5，100］m，假设有2个目标声源先后出现在监测区域内做匀速直线运动，具体运动状态见表2，真实运动轨迹见图1(a)，仿真中的波导环境参数见图1(b)。

图1 左：目标真实轨迹与阵列监测区域；右：仿真波导环境参数图Fig.1 Left：True trajectory of targets and detection area of arrays； Right：Environmental parameters of simulation

表2 目标运动状态Tab.2 Motion state of two targets

2 个目标声源都在水下50 m 处定深航行，信号频率为200 Hz、300 Hz、400 Hz，声源级为130 dB，使用简正波模型计算水平-深度维的传播损失。仿真中使用3 个水平线阵（Horizontal Line Array，HLA）接收声压信号，3 个HLA 均平行于x轴，阵元数为24，总长度为128 m，深度为100 m。3 个HLA 的阵元中心分别位于（1，1） km、（2.5，2） km 和（4，1.5） km 处，具体位置见图1。通过设置噪声方差使HLA2 在各频率处的信噪比等于10 dB。HLA1 和HLA3 的噪声方差与HLA2 相同，计算得到的3 个HLA 的信噪比见表3。仿真中使用滑动窗口生成每一帧处理的频率快拍，快拍时间长度为T=12 s。一帧信号包含11 个快拍，两帧信号之间有5 个快拍重叠，监测历时100帧。

表3 各HLA的信噪比（单位：dB）Tab.3 Array-based signal-to-noise ratio of each horizontal line array （unit：dB）

设置监测区域的x、y维度拷贝场网格间距为100 m，z（深度）维度网格间距为5 m。截取深度为50 m 的模糊函数x-y截面分析MFP 算法的定位性能。图2展示了HLA2 MFP结果（k=21和k=59两帧），其中伪彩图表示归一化的模糊函数x-y截面。白色圆圈表示2 个目标的真实位置，白色十字表示估计位置，白色横线表示HLA 的位置。截面上还存在多个明显的亮色位置，这些位置将影响定位结果，导致偏差。

图2 归一化模糊函数x-y截面（伪彩图）：目标真实位置（白色圆圈）、估计位置（白色十字）和HLA的位置（白色横线）Fig.2 x-y intersection of the ambiguity function （pseudo-color）：true positions of two targets （white circles），Estimated positions of two targets （white crosses），positions of three arrays （white horizontal lines）

4.2 阈值对跟踪性能的影响

本文采用的跟踪算法基于SMC 方法实现。不失一般性，将每个伯努利分量的采样粒子数设置为最小值Lmin=300，最大值Lmax=1000，目标存活概率PS=0.99，目标检测概率PD=0.98。不同阈值下MFP-CBMB方法的单次跟踪情况如图4所示。当阈值较大时，杂乱的噪声干扰会使跟踪结果相对于真实位置出现较大偏移；当阈值逐渐减小时，算法对噪声干扰的滤除效果变好，但过小的阈值又会使目标轨迹的跟踪不连续。综合考虑上述问题，选取α=0.5作为分布式融合实验的阈值参数。

图4 不同阈值下MFP-CBMB方法的单次跟踪结果Fig.4 Single-time tracking results of MFP-CBMB by setting different thresholds

从跟踪精度角度出发，监测区域内目标声源的数目无具体限制，只要在同一快拍中不同目标声源的位置处于MFP 中的不同网格点处，在分辨力允许的情况下可以分辨即可。分辨力需要根据监测环境设置，且不同监测环境中目标声源的运动范围和运动速度不同，因此无法对不同的监测环境给出通用的可分辨数量限制，需要根据具体环境分析。从计算效率和通信负担角度出发，每帧中需要传输的后验概率密度中的伯努利分量总数有最大限制Mmax，要求目标声源总数≤Mmax。

4.3 分布式融合实验及与集中式融合方法的性能对比

本节通过分布式融合实验验证论文所提方法的有效性。监测区域内3 个HLAs 组成一个分布式声学传感器网络，HLA 都固定放置在海底，并采取水声通信方式进行信息交互，每个HLA 进行独立的MFP 和量测筛选处理，且仅与邻近的HLA 进行数据传输。网络中2条链路的直线距离分别约为1.8 km、1.6 km。分布式融合过程中，HLA1和HLA3分别对两个后验概率密度（自身和HLA2）进行融合，HLA2对三个后验概率密度进行融合（自身、HLA1 和HLA3）。HLA 的融合权重使用 Metropolis 权重计算获得，如表4。

表4 各个HLA的融合权重Tab.4 Fusion weights of each horizontal line array

仿真设置OSPA 距离参数c=1000，p=1。图5给出了100 次蒙特卡洛实验中，单节点MFP-CBMB方法和多节点MFP-GCI-CBMB 方法的平均OSPA 距离随时间变化的情况。图5（左）展示了HLA1 单节点及2 个节点分布式融合处理的平均OSPA 距离，图5（右）展示了HLA2 单节点及3 个节点分布式融合处理的平均OSPA 距离。在目标出现和消失的邻近帧中，分布式融合能够有效缓解目标数目改变导致的跟踪精度突增问题。在目标1 出现的第10 帧，相较于单节点处理，分布式融合方法的平均OSPA距离由1 km 下降到0.75 km；在目标2 出现的第20 帧，由于跟踪误差的累积，单节点处理的平均OSPA距离在持续7～8帧内均保持在较高水平，而分布式融合方法可以使突增的跟踪误差较快回落，且回落时间不超过2～3 帧。这表明，分布式融合MFPGCI-CBMB 方法可以有效降低跟踪误差的突升水平，提高跟踪误差突升的回落速度。

图5 单节点处理（local）与分布式融合方法（GCI）的平均OSPA距离.左：HLA1；右：HLA2Fig.5 Average OSPA distance between local processing （local） and distributed fusion method （GCI）.Left：HLA1； Right：HLA2

分布式融合的平均OSPA 距离呈现总体下降状态。表5 展示了监测时间内HLA1 与HLA2 的平均单帧OSPA距离。当本地节点为HLA1时，分布式融合平均OSPA 距离从150.3 m 下降到101.1 m；当本地节点为HLA2 时，分布式融合平均OSPA 距离从119.6 m 下降到90.1 m。统一参照HLA1 进行对比，融合2 个HLA 的平均OSPA 距离下降了33%，融合3 个HLA 的平均OSPA 距离下降了41%。相比于单节点处理，MFP-GCI-CBMB 提高了监测时间内的跟踪精度，验证了在分布式融合处理中的有效性。

表5 平均单帧OSPA距离（单位：m）Tab.5 Average OSPA distance per frame （unit：m）

图6 （左）展示了平均OSPA 距离随信噪比的变化情况，当信噪比较低时，跟踪误差较大，随着信噪比的上升，跟踪误差稳定到0.1 km 左右。需要注意的是，MFP 处理使用离散网格点对监测区域进行划分，即便信噪比较高，也会存在固有的定位误差。由于仿真设置帧内时间长度为72 s，且存在2 个目标声源，因此该误差与目标声源的运动速度相比可以接受。实际应用中信噪比阈值需要根据具体场景对跟踪精度的宽容程度确定，与监测区域大小、帧内时间长度、目标声源数目以及运动速度等影响因素有关。例如，在仿真实验中，为了获得较好的跟踪效果，如平均OSPA距离控制在0.2 km内，信噪比要求≥-5 dB。集中式融合结构将HLA2作为中心节点，将HLA1和HLA3作为其他节点。根据3.2节，MFP-CEN-CBMB-1方法中，HLA1和HLA3将模糊函数传输至HLA2，MFP-CEN-CBMB-2 方法中，主要传输模糊函数中的所有峰值。分布式融合结构以HLA2 作为本地节点，融合其邻近节点（HLA1、HLA3）的后验概率密度，取MFP-GCI-CBMB 方法的粒子权重和阈值ω=0.99。对于集中式融合方法，中心节点HLA2 接收其他节点（HLA1、HLA3）的数据，对数据做归一化，选取大于阈值α=0.5 的峰值点的坐标作为量测，送入CBMB方法中进行滤波。经阈值筛选后，两种集中式方法的量测相同，滤波结果也相同。如表5 所示，在监测时间内，集中式、分布式融合方法的平均单帧OSPA 距离分别为61.1 m 和90.1 m。集中式融合方法中心节点获取的量测信息更为原始完整，因此其跟踪精度更高，而分布式融合方法中本地节点获取的是量测滤除干扰后的后验概率密度，在通信负担大幅降低的情况下，可达到与集中式融合方法相当的跟踪精度。

图6 （左） MFP-GCI-CBMB方法在不同信噪比下的平均跟踪性能；（右） 集中式、分布式融合方法平均每条数据链路的通信量Fig.6 （Left） Average tracking performance of MFP-GCI-CBMB method under different signal-to-noise ratio； （Right） Average traffic load of centralized and distributed fusion methods per data link

三种融合方法的单条链路平均每帧通信量如图6（右）所示，分别约为203.2 KBytes、50.6 KBytes以及1.9 KBytes，均较平稳。两帧时间间隔为72 s，对应通信速率分别约为22.6 kbps、5.6 kbps 以及0.2 kbps。相较于两种集中式融合方法，分布式融合方法传输的通信量大幅下降，减轻了系统的通信负担。表6展示了三种融合方法在每个节点上的平均单帧处理时间（包括量测获取、滤波和融合时间），实验采用MATLAB 2016b 仿真软件，处理器为Intel i7-4790、3.60 GHz，内存为8GB。分布式融合方法的优点在于，三个节点的处理时间较为均衡，计算负担最重的本地节点HLA2的处理时间并不明显高于其邻近节点HLA1、HLA3；集中式融合方法的优点在于，中心节点HLA2 承担了大部分的计算负担，其处理时间约为其他节点的2～3倍，而其他节点的计算负担较轻，其中MFP-CEN-CBMB-1 方法中其他节点的计算负担更轻一些。

表6 集中式、分布式融合方法的平均单帧数据处理时间（单位：s）Tab.6 Average processing time of centralized and distributed fusion methods per frame （unit：s）

综合考虑计算效率，集中式、分布式融合方法各有优势。若网络中存在硬件计算较强的节点，采用集中式融合方法，将此节点作为中心节点，分担其他普通节点的计算负担。但需要考虑中心节点失效造成的网络瘫痪；若网络中节点的硬件计算相差不大，更适合采用分布式融合方法，使节点的计算负担较为均衡。

5 结论

本文提出了一种基于匹配场定位量测的分布式势平衡多伯努利联合定位与跟踪方法。针对分布式处理中单节点可获取信息受限问题，建立了一种基于水下环境匹配场定位的适合水声多目标跟踪的量测模型，选取合适的阈值筛选模糊度图上的峰值，结合CBMB滤波算法，降低了噪声干扰引起的定位误差。基于GCI融合法则序贯融合本地节点与其邻近节点的后验概率密度，充分利用了分布式水声传感器网络中各节点的量测信息，有效提高了跟踪精度，并可推广至具有更多节点的网络。仿真结果表明，相比于单节点处理方法，分布式融合方法的跟踪精度有所提升；相比于集中式融合方法，网络的通信负担更低，计算效率更高。