成像声呐可以在昏暗或浑浊水域中生成几乎等同于光学影像质量的高清晰度图像，因此被广泛用于各类水下作业，例如渔业管理、结构检测、管道泄漏、水底探测、水下搜寻、水下安检等^[1-3]. 通过成像声呐获取水下鱼群的声学影像，从影像中提取个体目标进行多目标跟踪研究是现在渔业的重要研究方向，既可以实现目标计数，又可以观测鱼类行为，这将为渔业资源评估提供技术支撑^[4]. 在水下多目标跟踪中，经常面临观测量不确定、探测声呐分辨力有限、目标密集及低门限检测等问题，导致跟踪效果不佳. Handegard等^[5]通过成像声呐的拖网作业对水下目标进行观测，利用 $\alpha - \beta $ 滤波对个体目标进行跟踪研究，该方法对于稀疏目标状态下的个体目标具有良好的跟踪效果. Han等^[6]利用成像声呐对网箱中转移的鱼类进行个体计数与体长估计，采用卡尔曼滤波（Kalman filtering，KF）算法对水下目标进行状态估计，该研究未对多目标跟踪进行深入研究. Shahrestani等^[7]利用自适应分辨率成像声呐定点采集水下信息，提出自动图像处理机制，可以高效、准确地进行目标数量估计与数据分析，该研究没有涉及到多目标关联的问题. Huang等^[8]利用双频识别声呐采集水下鱼类信息，通过声学图像处理方法提取个体目标，采用扩展卡尔曼滤波（extended Kalman filtering，EKF）算法对非线性状态下的运动目标进行跟踪研究，重点分析鱼类运动的各种参数. 徐盼麟等^[9]利用光学摄像头对水下鱼类进行三维跟踪，该研究采用交互式多模型，结合联合概率数据关联算法（interacting multiple model joint probabilistic data association，IMMJPDA）对多目标进行有效跟踪. 该算法的计算量与目标数的阶乘成正比，因此在密集回波环境下，计算量呈指数规律增长.

针对该问题，本文提出基于成像声呐的水下多目标跟踪方法. 通过定点采集水下鱼群信息，建立基于信号强度的目标模型，提取个体目标；采用基于序贯蒙特卡罗的概率密度假设滤波，结合Auction航迹识别算法（sequential Monte Carlo probability hypothesis density combined with Auction track recognition，SMCPHD-Auction）对水下多目标进行跟踪，通过仿真分析该算法的可行性，开展现场试验.

该研究采用的成像声呐是利用声学透镜发射独立波束的多波束系统——双频识别声呐（dual-frequency identification sonar，DIDSON）. 相对于传统的多波束声呐系统采用时延阵列或者数字波束形成技术进行信号的收发，双频识别声呐通过声学透镜实现波束形成，既降低系统功耗、减小体积，又可以在水下生成高清晰度声学图像，达到识别目标的程度^[10-11].

在每个工作周期，DIDSON发射96个波束，声波遇到水下障碍物返回，阵元接收回波信号组成一帧声呐图像，如图1所示. 如图1（a）所示为96个阵元同时发射波束的示意图；如图1（b）所示为声呐发射波束的覆盖范围，其中水平视场角为 $\alpha = {29^ \circ }$，垂直视场角为 $\beta = {14^ \circ }$；如图1（c）所示为波束探测范围内的物体呈现在声学影像图上的示意图.

声呐图像中每个像素点代表回波信号的强度，对每个像素点建立信号强度模型：

(1) $I = \bar I + \sigma \sin \left( {\omega t} \right) + k\zeta .$

式中： $I$、 $\bar I$、 $\sigma $、 $\omega $、 $\zeta $ 分别为图像中某个像素点的强度、该像素点在连续若干帧图像中的平均强度、强度变化幅度、强度变化角频率及噪声等级； $k$ 为噪声等级系数，通常取为1，因此当图像中像素点的强度符合式（2）时，即为背景：

(2) $\bar I - \sigma - \zeta \leqslant I \leqslant \bar I + \sigma + \zeta .$

由于背景的强度低于目标的强度，目标的强度可以判决为

(3) $I > \bar I + \sigma + \zeta .$

$\bar I$ 和 $\sigma $ 通过下式进行更新：

(4) $\bar I' = {{( {n - 1}) }}\bar I /{n}+ I/{n},$

(5) $\sigma ' = {{ ({n - 1}) }}\sigma/{n} + {\sqrt {2{{\left( {I - \bar I} \right)}^2}}}\Bigg/{n} .$

式中： $\bar I'$ 为 $\bar I$ 更新后的值； $\sigma '$ 为 $\sigma $ 更新后的值； $n$ 为正整数，根据图像背景的复杂程度取值.

基于数据关联的多目标跟踪算法，例如联合概率数据关联（joint probabilistic data association，JPDA）、多假设跟踪（multiple hypothesis tracking，MHT）等，该算法的计算量随着目标数、量测数呈指数增长，因此当水下目标密集、杂波较多时，多目标跟踪的计算量呈爆炸式增加^[12-14]. 为了提高目标跟踪的效率、降低计算量，该研究采用基于随机集理论的目标跟踪算法，即利用概率假设密度（probability hypothesis density，PHD）滤波算法实现目标状态估计与目标数量估计，从而避免数据关联所需的庞大计算，其中PHD滤波器利用粒子滤波实现^[15].

首先建立多目标跟踪系统的动态模型：

(6) $\begin{aligned} {{{x}}_k} = f\left( {{{{x}}_{k - 1}}, {{{n}}_{k - 1}}} \right), \; {{{z}}_k} = h\left( {{{{x}}_k}, {{{v}}_k}} \right). \end{aligned} $

式中： $f\left( \cdot \right)$ 为状态转移函数； $h\left( \cdot \right)$ 为观测函数； ${{{x}}_k}$ 和 ${{{z}}_k}$ 分别为目标状态量和观测量，对应的系统过程噪声为 ${{{n}}_{k - 1}}$，观测噪声为 ${{{v}}_k}$.

假设 $k$ 时刻系统的后验概率密度函数可以用一组加权粒子 ${{{X}}_k} = \left\{ {{{x}}_k^i, w_k^i} \right\}_{i = 1}^{{L_k}}$ 表征：

(7) ${f_k}\left( {{{{x}}_k}|{{{Z}}_{1:k}}} \right) \approx \sum\nolimits_{i = 1}^{{L_k}}\; {w_k^i} \delta ({{{x}}_k} - {{x}}_k^i).$

式中： $w_k^i$ 为权值， ${L_k}$ 为粒子个数， $\delta ( \cdot )$ 为狄拉克函数， ${{{Z}}_{1:k}}$ 为 ${{{z}}_0}$ 至 ${{{z}}_k}$ 的集合. SMCPHD滤波算法的实现步骤如下.

1）初始化. 当 $k = 0$ 时，先验分布和每个目标用粒子 ${{{X}}_0} = \left\{ {{{x}}_0^i, w_0^i} \right\}_{i = 1}^{{L_0}}$ 表示，假定估计目标数目为 ${\hat N_0}$，则 $w_0^i = {{{{\hat N}_0}} /{{L_0}}}$， $x_0^i$ 通过初始PHD函数 ${D_0}$ 获得：

(8) ${D_0}({{{x}}_0}) = \sum\nolimits_{i = 1}^{{L_0}} \;{w_0^i} \delta ({{{x}}_0}).$

2）预测. 当 $k \geqslant 1$ 时，选定 ${p_k}( \cdot |{{{z}}_k})$ 和 ${q_k}( \cdot |{{{x}}_{k - 1}}, {{{z}}_k})$ 这2个概率密度函数作为重要性函数，其中 ${p_k}( \cdot |{{{z}}_k})$ 用来产生新生目标强度函数 ${\gamma _k}({{{x}}_k})$ 的粒子， ${q_k}( \cdot |{{{x}}_{k - 1}}, {{{z}}_k})$ 用来产生衍生目标的强度函数和继续存活的目标强度函数的粒子. 通常情况下，选取状态转移概率密度函数为 ${q_k}( \cdot |{{{x}}_{k - 1}}, {{{z}}_k})$，新生目标概率密度函数为 ${p_k}( \cdot |{{{z}}_k})$，因此有

(9) $\tilde {{x}}_k^i \sim \left\{ \begin{array}{l} {q_k}( \cdot |{{x}}_{k - 1}^i, {{{z}}_k}), \;i = 1, \cdots , {L_{k - 1}};\\ {p_k}( \cdot |{{{z}}_k}),\;i = {L_{k - 1}} + 1, \cdots , {L_{k - 1}} + {J_k}. \end{array} \right.$

对应粒子的权值为

(10) $\tilde w_{k|k - 1}^i = \left\{ \begin{gathered} \frac{{{\phi _{k|k - 1}}({{x}}_k^i|{{x}}_{k - 1}^i)w_{k - 1}^i}}{{{q_k}({{x}}_k^i|{{x}}_{k - 1}^i, {{{z}}_k})}}, \;i = 1, \cdots , {L_{k - 1}}; \\ \frac{{{\gamma _k}({{x}}_k^i)}}{{{J_k}{p_k}({{x}}_k^i|{{{z}}_k})}},\; i = {L_{k - 1}} + 1, \cdots , {L_{k - 1}} + {J_k}. \\ \end{gathered} \right.$

式中：

$ {\phi _{k|k - 1}}({{{x}}_k}|{{{x}}_{k - 1}}) = {\beta _{k|k - 1}}({{{x}}_k}|{{{x}}_{k - 1}}) + {e_{k|k - 1}}({{{x}}_{k - 1}})$ $p({{{x}}_k}|{{{x}}_{k - 1}}),$

其中 ${\beta _{k|k - 1}}({{{x}}_k}|{{{x}}_{k - 1}})$ 为衍生目标的PHD函数， ${e_{k|k - 1}}$ $({{{x}}_{k - 1}})$为目标存活概率，因此当前时刻的PHD函数为

(11) ${D_{k|k - 1}}({{{x}}_k}|{{{Z}}_{1:k - 1}}) \approx {\hat D_{k|k - 1}}({{{x}}_k}|{{{Z}}_{1:k - 1}}) = \sum\limits_{i = 1}^{{L_{k - 1}} + {J_k}} {\tilde w_{k|k - 1}^i\delta ({{{x}}_k})} .$

3）更新. 当系统获得最新观测集合后，依据预测步骤得到粒子集合 $\{ w_{k|k - 1}^i, {{x}}_k^i\} _{i = 1}^{{L_{k - 1}} + {J_k}}$ 近似的PHD函数 ${D_{k|k - 1}}({{{x}}_k}|{{{Z}}_{1:k - 1}})$；接着进行粒子更新，每个粒子的更新权重为

(12) $\tilde w_k^i = \left[ {1 - {p_{D, k}}({{x}}_k^i) + \sum\limits_{z \in {Z_k}} {\frac{{{\psi _{k, z}}({{x}}_k^i)}}{{{\kappa _k}({{{z}}_k}) + {C_k}({{{z}}_k})}}} } \right]\tilde w_{k|k - 1}^i.$

式中： ${\psi _{k, z}}({{x}}_k^i) = {p_{D, k}}({{x}}_k^i)p({{{z}}_k}|{{x}}_k^i)$，

${C_k}({{{z}}_k}) = \displaystyle\sum\limits_{i = 1}^{{L_{k - 1}} + {J_k}} $ $ {{\psi _{k, z}}({{x}}_k^i)} \tilde w_{k|k - 1}^i.$

4）重采样. 计算粒子集合的权值累加和，即目标个数：

(13) ${\hat N_k} = \sum\nolimits_{i = 1}^{{L_{k - 1}} + {J_k}} {\tilde w_k^i} .$

判断粒子集合 $\left\{ {{{\tilde w_k^i} /{{{\hat N}_k}}}, \tilde {{x}}_k^i} \right\}_{i = 1}^{{L_{k - 1}} + {J_k}}$ 是否需要重采样，定义有效样本容量参数 ${N_{\rm eff}} = {1 {\left/\displaystyle\sum{_{i = 1}^N}\right. {{{\left( {\tilde w_k^i} \right)}^2}} }}$ 以及阈值 ${N_{\rm th}}$. 若 ${N_{\rm eff}} < {N_{\rm th}}$，则进行重采样，得到新的粒子集合 $\left\{ {{{w_k^i} /{{{\hat N}_k}}}, {{x}}_k^i} \right\}_{i = 1}^{{L_k}}$. 接着进行尺度变换，乘以 ${\hat N_k}$，得到重采样后的粒子集 $\left\{ {w_k^i, {{x}}_k^i} \right\}_{i = 1}^{{L_k}}$.

5）状态估计输出.

通过加权准则确定PHD函数：

(14) ${D_k}({{{x}}_k}|{{{Z}}_{1:k}}) = \sum\nolimits_{i = 1}^{{L_k}} {w_k^i\delta ({{{x}}_k})} .$

对所有粒子权值求和，获得目标数目估计值：

(15) ${\hat N_k} = \sum\nolimits_{i = 1}^{{L_k}} {w_k^i} .$

对于 ${\hat N_k}$ 个目标的状态量，采用 $k $-means 分类技术进行提取，步骤如下. 1）从当前时刻目标状态有限集中，任意选取 ${\hat N_k}$ 个数据作为初始聚类中心；2）计算剩余数据样本到各聚类中心的距离，依据计算结果将它们归类到最近的聚类中；3）计算每个聚类中所有数据的平均值，将该均值作为新聚类的中心；4）不断重复该过程直到相邻2次聚类中心没有任何变化，即样本已经收敛，选中的聚类中心为当前时刻的目标状态量^[16].

由于SMCPHD滤波获得的目标状态不具有身份标志，即已有的目标轨迹与当前时刻目标状态没有关联. 为了对水下目标进行良好的跟踪，设计基于Auction的航迹识别算法^[17]. 与二维分配算法相似，将当前 $k$ 时刻的目标状态 ${{{X}}_k}$ 分配给前一时刻已确认的目标航迹 ${{{T}}_{k - 1}}$，形成关联矩阵；依据最优分配准则计算 ${{{X}}_k}$ 与 ${{{T}}_{k - 1}}$ 中目标状态的关联代价，总的关联代价和最小时表示所有项目的匹配程度最高，可以将当前时刻的目标状态与已确认的每条航迹关联，航迹识别算法原理如图2所示.

假设 $k - 1$ 时刻的航迹 ${{{T}}_{k - 1}}$ 包含 ${L_{k - 1}}$ 个航迹 $\left\{ {{{t}}_{k - 1}^m} \right\}_{m = 1}^{{L_{k - 1}}}$， $k$ 时刻的目标状态 ${{{X}}_k}$ 包含 ${N_k}$ 个状态估计 $\left\{ {{{x}}_k^n} \right\}_{n = 1}^{{N_k}}$， ${{{T}}_{k - 1}}$ 与 ${{{X}}_k}$ 的关联目标函数定义为 $C =$ $\min \displaystyle\sum{_{m = 1}^{{L_{k - 1}}}} \sum\nolimits_{n = 1}^{{N_k}} {{c_{m,n}}{a_{m,n}}} $，其中 ${c_{m, n}}$ 为代价系数， ${a_{m, n}}$ 为二值变量，且满足： $\displaystyle\sum{_{n = 1}^{{N_k}}} {{a_{m, n}}} = 1 (m = 1, 2, \cdots , {L_{k - 1}})$ 与 $\displaystyle\sum{_{m = 1}^{{L_{k - 1}}}} {{a_{m, n}}} = 1(n = 1, 2, \cdots , {N_k})$. 代价函数采用巴氏距离衡量，距离越小，表示2个物理量越相似、关联代价越小，若 $k$ 时刻的一个目标状态量为 ${{x}}_k^j = \left[ {x_k^j, \dot x_k^j, y_k^j, \dot y_k^j} \right],\;k - 1$ 时刻的一条航迹为 ${{t}}_{k - 1}^i =$ $\left[ {x_{k - 1}^i, \dot x_{k - 1}^i, y_{k - 1}^i, \dot y_{k - 1}^i} \right]$，巴氏距离定义为

(16) $d\left( {{{t}}_{k - 1}^i, {{x}}_k^j} \right) = \sqrt {1 - \rho \left( {{{t}}_{k - 1}^i, {{x}}_k^j} \right)} .$

式中： $\rho \left( {{{t}}_{k - 1}^i, {{x}}_k^j} \right)$ 为巴氏系数，

$ \left.\rho \left( {{{t}}_{k - 1}^i, {{x}}_k^j} \right) =\displaystyle\sum\nolimits_{s = 1}^4 {{A_s}{B_s}}\right/{\left( {A_s} + {B_s}\right)}^2$，

其中 ${A_s}$ 与 ${B_s}$ 的值分别取自航迹状态量和目标状态量的元素值，因此关联代价系数为 $ {c_{i, j}} = $ $ \exp\,\left[ { - d\left( {{{t}}_{k - 1}^i, {{x}}_k^j} \right)} \right]$.

航迹识别算法的具体过程如图3所示，计算过程如下. 1）计算关联矩阵 $\left[ {{c_{i, j}}} \right]$，其中每个元素代表航迹 ${{{T}}_{k - 1}}$ 与目标状态 ${{{X}}_k}$ 的关联代价；2）将每条航迹的初始代价清零，即 ${P_n} = 0$；3）找到 ${{{X}}_k}$ 中的目标状态 ${{x}}_k^j$ 的最佳匹配航迹 ${{t}}_{k - 1}^i$，匹配原则：

${c_{i, j}} - {P_j} = \mathop {\max }\limits_{j = 1, 2, \cdots , {L_{k - 1}}} \left\{ {{c_{i, j}} - {P_j}} \right\}$；

4）判断 ${{{X}}_k}$ 是否全部匹配成功；5）将3）中已匹配的目标状态剔除；6）更新每条航迹的代价： ${P_j} =$ $ {P_j} + {d_n} + \xi $，其中 ${d_n}$ 为目标状态估计的最优与次最优分配代价的差值， $\xi $ 为代价常数；7）输出当前时刻的航迹 ${{{T}}_k}$.

如果 ${N_k} < {L_{k - 1}}$，执行完步骤1）~7）后，若有航迹剩余，则表示该条航迹对应的目标消失；如果 ${N_k} \geqslant {L_{k - 1}}$，执行完步骤1）~7）后，若有目标状态剩余，则表示新目标产生.

为了方便算法仿真，假设所有目标都在 $\left[ { - 100, 100} \right]{\rm m} \times \left[ { - 100, 100} \right]$ m的平面上运动，每个目标的运动都服从如下的线性高斯模型：

(17) $ {{{x}}_k} = \left[ {\begin{array}{*{20}{c}} \!\! 1&T&0&0 \!\!\\ \!\!0&1&0&0 \!\!\\ \!\!0&0&1&T \!\!\\ \!\! 0&0&0&1 \!\! \end{array}} \right]{{{x}}_{k - 1}} + \left[ \begin{array}{*{20}{c}} \!\! {{T^2} /2}\!\! & \!\!T \!\!\\ \!\!T\!\!&\!\!0\!\!\\ \!\!0\!\! &\!\! {{T^2} /2} \!\!\\ \!\! 0\!\!&\!\!T\!\! \end{array} \right]\left[ {\begin{array}{*{20}{c}} \!\!\!{{v_{1, k}}}\!\!\! \\ \!\!\!{{v_{2, k}}} \!\!\! \end{array}} \right].\!\!\!\!\!$

式中： ${{{x}}_k} = {\left[ {x, \dot x, y, \dot y} \right]^{\rm T}}$，采样周期 $T = 1$ s；过程噪声 $\left\{ {{v_{1, k}}} \right\}$、 $\left\{ {{v_{2, k}}} \right\}$ 是彼此独立的零均值高斯白噪声，各自的标准差为 ${\sigma _{v1}} = 1$， ${\sigma _{v2}} = 0.1$；每个目标的生存概率为 ${e_{k|k - 1}} = 0.95$，假设没有衍生目标；杂波采用泊松分布模型产生，每次扫描到的杂波平均数 ${\lambda _k} = 5$，杂波均匀地分布在 $\{ {\left( {x, y} \right)|x \in \left[ { - 100, 300} \right]}, $ $ {y \in \left[ { - 100, 100} \right]} \}$ 的平面内；目标的检测概率 ${p_{\rm D}} = 0.99$；新目标的产生服从泊松过程，强度函数为 ${\gamma _k} = 0.2N\left( { \cdot ;\bar x, {{Q}}} \right)$，其中 $N\left( { \cdot ;\bar x, {{Q}}} \right)$ 表示均值为 $\bar x$、方差为 ${{Q}}$ 的正态分布， $\overline x $ 为每个时刻该目标的状态量， ${{Q}}$ 都相同，取值为

${{Q}} = \left[ {\begin{array}{*{20}{c}} {10}&0&0&0 \\ 0&1&0&0 \\ 0&0&{10}&0 \\ 0&0&0&1 \end{array}} \right]$，

目标的量测方程为

(18) ${{ y}_k} = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 1&0&0&0 \end{array}} \\ {\begin{array}{*{20}{c}} 0&0&1&0 \end{array}} \end{array}} \right]{{{x}}_k} + \left[ {\begin{array}{*{20}{c}} {{w_{1, k}}} \\ {{w_{2, k}}} \end{array}} \right].$

式中： $\left\{ {{w_{1, k}}} \right\}$ 和 $\left\{ {{w_{2, k}}} \right\}$ 是彼此独立的零均值高斯白噪声，各自的标准差分别为 ${\sigma _{w1}} = 2$， ${\sigma _{w2}} = 2$；测量噪声和过程噪声是相互独立的. 对于每个可能的目标，固定分配 $\rho $ =1 000个粒子，即每个采样间隔加入1 000个新生粒子，以表征可能产生的新目标，新粒子按照 $p \sim 0.2N\left( { \cdot ;\bar x, {{Q}}} \right)$ 的正态分布生产.

假设共有3个目标，各目标的初始位置与出现、消失情况如表1所示，其中目标1在初始时刻出现，目标2和目标3是新生目标.

经过SMCPHD滤波，绘出目标沿着 $x$、 $y$ 轴跟踪的效果图，如图4所示. 如图4（a）、（b）所示为不同时刻下观测量分别沿着 $x$、 $y$ 轴的分布，包括3个目标的有效观测值和杂波，即图中的小圆圈；如图4（c）、（d）所示为经过滤波后，目标的状态估计值与真实轨迹的对照. 可以看出，目标的估计量与真实轨迹基本吻合，证明利用SMCPHD滤波可以有效地降低虚警率，实现多目标的准确跟踪.

将目标跟踪结果显示在笛卡尔坐标系下，如图5所示为3个目标的真实运动轨迹以及经过SMCPHD滤波后的位置. 如图6所示为不同时刻下目标个数的估计值与真实个数的对比. 图中，n为目标个数. 可以看出，目标跟踪效果良好，目标数量估计准确.

基于SMCPHD滤波算法获得的目标状态估计值与真实目标状态量的偏差可以通过Wasserstein距离进行定量评估. 假设真实目标状态集为 ${{X}}$，估计值为 $\hat {{X}}$，且均非空，则Wasserstein距离定义为

(19) ${d_p}(\hat {{X}}, {{X}}) = \mathop {\min }\limits_c\left( {{\sum\limits_{i = 1}^{\left| {\hat {{X}}} \right|} {\sum\limits_{j = 1}^{\left| {{X}} \right|} {{C_{i, j}}{{\left\| {{{\hat x}_i} - {x_j}} \right\|}^p}} } }}\right)^{1/p}.$

式中： $\mathop {\min }\limits_c \left( \cdot \right)$ 表示在所有的传输矩阵C构成的集合中取最小的集合，元素 ${C_{i, j}}$ 满足条件： ${C_{i, j}} \geqslant 0, $ $ \displaystyle\sum\nolimits_{j = 1}^{\left| {{X}} \right|} {{C_{i, j}}} = {1}/{{\left| {\hat {{X}}} \right|}}, \displaystyle\sum\nolimits_{i = 1}^{\left| {\hat {{X}}} \right|} {{C_{i, j}}} = {1}/{{\left| {{X}} \right|}}$；

$p$ 为Wasserstein系数，本文取值为2. 图7显示了在整个时间序列中，Wasserstein距离d不超过4，相对于目标的运动范围 $\left[ {0, 140} \right] \times \left[ { - 60, 30} \right]$ m，目标跟踪误差很小.

对提取的目标状态集，采用Auction航迹识别算法进行航迹识别，分别标示出轨迹1、2和3，结果如图8所示.

为了定量比较提出的多目标跟踪算法与基于数据关联的目标跟踪算法的效率，开展数据仿真. 设定有效目标数为5，蒙特卡罗仿真次数为50步，每个时刻的量测平均数为20，SMCPHD-Auction算法中粒子数为1 000，MHT算法中假设生成树深度为3，利用Matlab进行仿真计算. 仿真环境如下：CPU inter（R） Core（TM） i5-4590 CPU 3.3 GHz，内存为16 GB，MATLAB-2015（b）. JPDA和MHT相对于SMCPHD-Auction的仿真运行时间比分别为2.1和8.3.

选定一个鱼塘作为试验场地，将DIDSON固定在试验场一侧的水下，利用定点采集模式进行水下数据采集. DIDSON布置示意图如图9所示，将声呐布置在距离水面约0.5 m深的水下，并以10°的倾角向下，使得声呐探测范围尽可能大地覆盖住水域的深度. 如图9（a）所示为DIDSON布置的俯视图，如图9（b）所示为DIDSON布置的截面图，即为垂直视场图. 调节DIDSON的探测距离为10 m，使得声呐探测范围基本覆盖水域的宽度，并将DIDSON接入岸边的计算机进行水下数据的记录与存储. DIDSON的具体参数设置如表2所示.

对采集的每帧数据进行图像处理，主要包括将 $96 \times 512$ 的二维数据还原成扇形图以及目标提取. 如图10（a）所示为经过图像还原后的扇形图，其中扇环两侧的数字代表测量范围，单位是m，图中的亮色区代表目标或者噪声，暗色区代表背景. 如图10（b）所示为经过背景去除即目标提取后得到的声呐图像.

利用SMCPHD-Auction算法对多目标进行跟踪，记录每个目标的运动轨迹. 表3记录了连续6帧图像中11个目标的位置，其中的位置 $\left( {\rho , \theta } \right)$ 对应图11所示的坐标系， $\rho $ 为目标的斜距， $\theta $ 为目标的方位角；6帧图像对应的时刻分别是 ${t_0}, {t_1}, \cdots , {t_5}$，且 ${t_{i + 1}} - {t_i} = \Delta T$，其中 $\Delta T$ 为连续两帧图像的时间间隔， $i = 0, 1, 2, 3, 4$. 根据极坐标与直角坐标系的转换关系，获得目标的直角坐标 $\left( {x, y} \right)$：

(20) $ \begin{gathered} x = \rho \cos \theta ,\; y = \rho \sin\theta . \end{gathered} $

由于每一个目标从出现到消亡会跨越若干帧图像，将不同帧图像中同一目标的位置依次用线段连接，并显示在同一幅图像中，是这个目标的运动轨迹. 为了显示所有目标的轨迹，需要保存目标在不同帧中出现的位置，并用不同颜色将不同目标轨迹进行标识. 为了使不同轨迹对应不同颜色，每条轨迹需要有一个颜色标记位，当新目标出现的时候，系统自动为该目标分配一种新颜色；当一个目标消亡时，系统回收该目标的颜色，这样可以保证同一帧图像中不同目标对应不同颜色，不同帧图像中，同一目标对应同样的颜色. 将表3中各目标的坐标依次连接，如图12（a）所示. 为了更好地展示不同时刻下多目标的跟踪轨迹，利用Visual Studio（VS）中的Microsoft Foundation Classes（MFC）模块结合Open Graphics Library（OpenGL），将连续多帧图像堆叠在同一空间中，不同的目标轨迹将呈现在三维空间中，如图12（b）所示. 此处不是真正的“三维空间”跟踪轨迹，仅仅是将原来堆叠在二维平面的6帧声呐图拉伸成三维，即各条轨迹不是由 $\left( {x, y, z} \right)$ 坐标构成，而是由 $\left( {x, y, t} \right)$ 坐标组成. 图12（b）中共有6帧连续图像按照时间顺序排列在空间中，时刻依次为 ${t_0}, {t_1}, \cdots , {t_5}$.

对于试验结果的验证，由于数据采集于人工渔场，声呐探测范围宽度约为10 m，深度约为4 m，而光学摄像头在水下很难捕捉到这么大范围的影像信息，因此可以在室内水池或水族馆等光线充足、背景简单的环境中利用光学摄像头进行对比试验.

提出基于声学信号强度模型的目标提取算法，对成像声呐采集的水下数据进行个体目标提取. 采用SMCPHD滤波器对目标状态与数量进行估计，利用Auction航迹识别算法关联当前时刻目标状态与已确定的航迹. 通过仿真发现，利用提出的水下多目标跟踪算法不仅可以准确获取目标状态与数量估计，而且大大降低了数据运算量，提高了跟踪效率. 采集试验场数据进行处理分析，利用OpenGL绘制出目标随时间变化的轨迹，真正实现了多目标跟踪.