一种联合粒子滤波和跟踪置前检测的水下弱目标跟踪方法

摘要

本发明涉及一种联合粒子滤波和跟踪置前检测的水下弱目标跟踪方法。针对水下目标辐射声信号微弱的问题，提出一种序贯粒子滤波和跟踪置前检测相结合的方法。粒子滤波用以解决序贯目标跟踪中状态空间模型的非线性和非高斯问题，跟踪置前检测方法用以解决微弱目标的检测跟踪难问题。该方法的步骤包括1)系统建模；2)产生初始化粒子；3)利用跟踪置前检测思想对粒子进行更新预测；4)粒子重采样，并循环步骤3处理新的数据。仿真实验表明，对于水下低信噪比目标，联合粒子滤波和跟踪置前检测方法能够实现有效跟踪，并优于常规的粒子滤波目标跟踪方法。

说明书

技术领域

本发明属于水下目标探测跟踪领域的弱目标探测跟踪技术，具体的是涉及一 种联合粒子滤波和跟踪置前检测水下弱目标跟踪方法。

背景技术

海洋覆盖了地球70％以上的面积，蕴藏了巨大的物质资源。由于水的特殊性， 海水对无线电波和可见光吸收太快，衰减迅速，从而不能实现能量的远距离传输， 制约了这些技术在海洋探测和水下目标跟踪方面的应用。声波是唯一可以在水下 进行远程传播的载体，不同的传播距离有不同的应用场景，如传播距离为米到百 米级可以利用声波进行高分辨率声成像，传播距离为公里到百公里级可以利用声 波进行探测和定位，传播距离为百公里到千公里级可以利用声波进行信息传输。 因此，声探测跟踪技术在声纳探潜、海洋生物跟踪等军事和民用领域具有广泛的 应用。

使用水听器阵列探测跟踪水下目标是当前最主要的手段。但海洋作为声传播 的物理空间，其构成要素包含海水、水中物体和海底，因此声传播与海洋环境的 变化过程有紧密联系，即海洋环境中的声传播受到温度变化、密度变化和盐度变 化等多种物理过程引起的声速剖面扰动影响很大，导致观测的传递函数是空时变 化的，给水下目标探测跟踪带来了挑战。同时随着科技水平的提高，水下军事目 标辐射强度越来越小，我们关注的目标特征越来越细微，需要更强有力的技术手 段对低信噪比目标进行检测跟踪。探测设备输出端的信噪比对目标探测和跟踪的 性能有很大的影响，当输出信噪比较小时，会对传统的先检测后跟踪方法带来挑 战，甚至使得目标检测失败，导致目标无法跟踪。特别是当目标距离远，目标辐 射噪声小时，被动声呐接收到的信噪比微弱，导致目标检测和跟踪失败。

发明内容

本发明针对水下弱目标的跟踪检测问题，提出一种联合粒子滤波和跟踪置前 检测方法，既能适用于非线性、非高斯的声传播模型，又能改善低信噪比目标的 检测跟踪性能，实现对水下弱目标的有效检测和跟踪，在弱目标跟踪方面具有较 好的实用性。

本发明是通过以下技术方案实现的，其具体步骤是：

步骤1：系统建模，包括目标运动建模和目标观测建模。系统建模的精确性 直接影响到跟踪方法的性能，通过分析可以确定较为合理的物理模型。

11)目标运动建模。受限于水流体的阻力作用，水下目标的运动速度一般都 较小，因此在较短的时间间隔内，其距离、深度等状态参量可近似为一阶线性变 化附加二阶的扰动(噪声)，故其状态方程为线性方程。这里假定声源点目标在 距离深度平面内作定深匀速运动，设声源状态为s，包含目标的深度z_s、距离r_s和速度v_s等信息，则有s＝[z_s,r_s,v_s]^T，于是声源的状态变化方程可用如下运动方 程表示：

s_t＝Fs_t-1+w_t    (1)

状态转移矩阵如下所示：

      $F=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& T\\ 0& 0& 1\end{array}\right)---\left(2\right)$      

其中，T为两次测量的时间间隔，w_t是符合高斯分布的过程噪声，其均值为零， 协方差矩阵Q，表示如下：

      $Q=\left(\begin{array}{c}{v}_{z_s}\\ T^2{v}_{a_s}/2\\ v_{a_s}\end{array}\right)---\left(3\right)$      

其中和分别表示在深度方向上和径向上的加速度扰动。

本发明在目标的状态变量中增加随机变量表征样本粒子存在状 态，其中e表示目标存在，表示目标不存在，该变量具有两种状态的马尔科夫 链，其马尔科夫转移矩阵为：

      $\Pi =\left(\begin{array}{cc}1-P_d& P_b\\ P_d& 1-P_b\end{array}\right)---\left(4\right)$      

其中P_b表示新生粒子概率，P_d表示死亡粒子概率。

12)目标观测建模。用y_t表示t时刻水听器阵的量测声压，设水听器阵有N 个阵元，则有y_t＝[p₀,p₁,…,p_N]^T，则量测方程为：

y_t＝h(s_t)+ω_t    (5)

其中，h(·)为声传播模型，ω_t为服从高斯分布的量测误差，其均值为零，协方差 矩阵为R。

为了适应粒子滤波和跟踪置前检测方法的应用场景，利用匹配场处理思想对 声压测量数据进行自相关处理。将声场按单位距离进行网格化，令z_t＝y_t′y_t表示 声场不同网格点的强度值，经过转换后的测量值类似于红外观测中的灰度图，其 表达式如下：

      $z_t=\left(\begin{array}{cc}{h}^{\prime }\left(s_t\right)h\left(s_t\right)+{w_t}^{\prime }{w}_t& E_t=e\\ {w_t}^{\prime }{w}_t& E_t=\bar{e}\end{array}\right)---\left(6\right)$      

步骤2：产生初始化粒子。我们定义初始化粒子为新生粒子即在上一时刻粒子存在状态为不存在，而在下一时刻为存在。状态向量从先验概 率密度函数q_b(s_k|z_k)中采样得到，建议分布q_b(s_k|z_k)服从观测域中z_k>γ条件区 域的均匀分布。水下运动目标的深度、距离和速度均服从均匀分布。包含：

(21)设定总的粒子数，粒子出生概率和死亡概率。

(22)根据处理得到的匹配场相关值，设置对应的门限值γ。当信躁比较大时， 门限相对可以严格，当信噪比较低时，门限相对可以宽松。

(23)获取所有相关值大于门限的状态量，并作为采样样本。

(24)按照均匀分布规律对所有的状态量进行采样，并作为初始化粒子。

步骤3：利用跟踪置前检测思想对粒子进行更新预测。对水听器阵列获取的 连续声压数据进行匹配场处理，得到目标状态的新息数据；针对水下弱目标检测 困难的特点，不对目标的存在性进行判断，而是按照一定的出生概率和死亡概率 对上一时刻的粒子状态进行预测，得到当前时刻的粒子状态；然后利用全部的新 息数据对不同类型的粒子进行更新，计算各个粒子的权重并输出当前状态量。包 含以下几个步骤：

(31)按照表1预测每个粒子下一时刻的存在变量

表1粒子存在变量更新流程

            

(32)对存活粒子状态进行预测更新，包含新生粒子和存活粒子两种，其中新生 粒子是上个时刻为死亡状态，而当前时刻为存活状态；存活粒子则是上个时刻和 当前时刻均为存活状态。

32a)新生粒子状态向量从q_b(s_k|z_k)中采样得到，此时先验概 率密度函数与新息相关。建议分布q_b(s_k|z_k)服从观测域中z_k>γ条件区域的均匀 分布。水下运动目标的深度、距离和速度均服从均匀分布。

32b)存活粒子状态向量通过状态转移方程计算更新。

(33)计算粒子的权重，进行归一化处理，输出当前估计状态值。包含两个步骤：

33a)计算粒子的重要性权值

      ${\tilde{w}}_k^n=\left(\begin{array}{cc}\underset{i\in C_i\left(s_k^n\right)}{\Pi }\underset{j\in C_j\left(s_k^n\right)}{\Pi }l\left(z_k^{(i,j)}\right|s_k)& ,E_k^n=1\\ 0& ,E_k^n=0\end{array}\right)---\left(7\right)$      

其中，$C_i\left(s_k^n\right)=\{i_0-p,...i_0-1,i_0,i_0+1,...i_0+p\},$i₀表示最接近状态向量位置的整 数，p表示目标对邻近分辨单元的影响程度。其中表示相关处理的似 然函数。定义互谱密度函数C＝E[yy′]＝hh′+R，则似然函数可表示为：

      $l\left(s\right)=\exp (-\frac{{\phi }_j\left(s\right)}{R})---\left(8\right)$      

      ${\phi }_j\left(s\right)=\mathrm{trC}-\frac{h^{\prime }\mathrm{Ch}}{h^{\prime }h}---\left(9\right)$      

其中trC是对矩阵C进行求迹运算。

33b)对粒子重要性权值进行归一化处理：

      $w_k^n=\frac{{\tilde{w}}_k^n}{\underset{n=1}{\overset{N}{\Sigma }}{\tilde{w}}_k^n}---\left(10\right)$      

33c)输出目标的状态估计向量

步骤4：粒子重采样，并循环步骤3处理新的数据。对于常规的粒子滤波而 言，粒子退化是其主要的缺陷，即随着迭代次数的增加，粒子的权重会集中在少 数几个粒子上，而其他大多数粒子的权值很小，结果导致粒子集无法表示真实的 后验概率密度。

为了解决粒子退化的问题，引入序贯重采样技术。通过引入有效粒子数N_eff表示粒子的退化程度，其计算公式如下：

      $N_{\mathrm{eff}}=\frac{1}{\underset{i=1}{\overset{N}{\Sigma }}{\left({\omega }_k^i\right)}^2}---\left(11\right)$      

同时，引入一个阈值N_threshold对有效粒子数进行判断，当有效粒子数小于阈 值时，则进行重采样，如此可以降低算法的计算复杂度。

本发明的优点是：可以有效的降低检测失败导致跟踪失败的概率。该方法是 一种非相参积累方法，强调跟踪处理的结果，也就是航迹提取。为了更好的提取 目标的观测信息，算法直接使用处理数据，而将检测判断作为处理环节的最后一 步，可以不对目标的出现与消失进行判断，而是通过概率的形式进行描述；同时， 由于没有对原始数据进行门限处理，不存在点迹概念，也就无需进行航迹关联， 有效的保证了跟踪的持续性。

粒子滤波本身能够处理非线性、非高斯系统问题，适应声传播波导环境具有 复杂多变的空时变化特点，因此基于粒子滤波的跟踪置前检测算法在水下弱目标 跟踪中具有广泛的应用前景。

附图说明

图1为本发明的系统模型示意图。

图2为信噪比等于-2dB时，初始粒子的分布情况图。

图3为信噪比等于-2dB时，跟踪稳定后粒子收敛情况图。

图4为信噪比等于-2dB时，本发明的跟踪方法的深度跟踪误差图。

图5为信噪比等于-2dB时，本发明的跟踪方法的距离跟踪误差图。

图6为信噪比等于-2dB时，常规粒子滤波处理的深度跟踪误差图。

图7为信噪比等于-2dB时，常规粒子滤波处理的距离跟踪误差图。

图8为信噪比等于-2dB时，本发明的跟踪方法的跟踪结果状态存在概率图。

具体实施方式

结合附图说明本发明的具体实施方式和所获得的性能。

本发明所述的联合粒子滤波和跟踪置前检测方法，包括：

步骤1：系统建模，包括目标的状态变化方程和观测方程两个方面。系统建 模的精确性直接影响到跟踪方法的性能，通过分析可以确定较为合理的物理模 型。

11)目标运动建模。受限于水流体的阻力作用，水下目标的运动速度一般都 较小，因此在较短的时间间隔内，其距离、深度等状态参量可近似为一阶线性变 化附加二阶的扰动(噪声)，故其状态方程为线性方程。这里假定声源点目标在 距离深度平面内作定深匀速运动，设声源状态为s，包含目标的深度z_s、距离r_s和速度v_s等信息，则有s＝[z_s,r_s,v_s]^T，于是声源的状态变化方程可用如下运动方 程表示：

s_t＝Fs_t-1+w_t    (1)

状态转移矩阵如下所示：

      $F=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& T\\ 0& 0& 1\end{array}\right)---\left(2\right)$      

其中，T为两次测量的时间间隔，w_t是符合高斯分布的过程噪声，其均值为零， 协方差矩阵Q，表示如下：

      $Q=\left(\begin{array}{c}{v}_{z_s}\\ T^2{v}_{a_s}/2\\ v_{a_s}\end{array}\right)---\left(3\right)$             

其中和分别表示在深度方向上和径向上的加速度扰动。

本发明在目标的状态变量中增加随机变量表征样本粒子存在状 态，该变量具有两种状态的马尔科夫链，其马尔科夫转移矩阵为：

      $\Pi =\left(\begin{array}{cc}1-P_d& P_b\\ P_d& 1-P_b\end{array}\right)---\left(4\right)$      

其中P_b表示新生粒子概率，P_d表示死亡粒子概率。

12)目标观测建模。用y_t表示t时刻水听器阵的量测声压，设水听器阵有N 个阵元，则有y_t＝[p₀,p₁,…,p_N]^T，则量测方程为：

y_t＝h(s_t)+ω_t    (5)

其中，h(·)为声传播模型，ω_t为服从高斯分布的量测误差，其均值为零，协方差 矩阵为R。

为了适应粒子滤波和跟踪置前检测方法的应用场景，利用匹配场处理思想对 声压测量数据进行自相关处理。将声场按单位距离进行网格化，令z_t＝y_t′y_t表示 声场不同网格点的强度值，经过转换后的测量值类似于红外观测中的灰度图，其 表达式如下：

      $z_t=\left(\begin{array}{cc}{h}^{\prime }\left(s_t\right)h\left(s_t\right)+{w_t}^{\prime }{w}_t& E_t=e\\ {w_t}^{\prime }{w}_t& E_t=\bar{e}\end{array}\right)---\left(6\right)$      

步骤2：产生初始化粒子。我们定义初始化粒子为新生粒子即在上一时刻粒子存在状态为不存在，而在下一时刻为存在。状态向量从先验概 率密度函数q_b(s_k|z_k)中采样得到，建议分布q_b(s_k|z_k)服从观测域中z_k>γ条件区 域的均匀分布。水下运动目标的深度、距离和速度均服从均匀分布。包含：

(21)设定总的粒子数，粒子出生概率和死亡概率。

(22)根据处理得到的匹配场相关值，设置对应的门限值γ。当信躁比较大时， 门限相对可以严格，当信噪比较低时，门限相对可以宽松。

(23)获取所有相关值大于门限的状态量，并作为采样样本。

(24)按照均匀分布规律对所有的状态量进行采样，并作为初始化粒子。

步骤3：利用跟踪置前检测思想对粒子进行更新预测。对水听器阵列获取的 连续声压数据进行匹配场处理，得到目标状态的新息数据；针对水下弱目标检测 困难的特点，不对目标的存在性进行判断，而是按照一定的出生概率和死亡概率 对上一时刻的粒子状态进行预测，得到当前时刻的粒子状态；然后利用全部的新 息数据对不同类型的粒子进行更新，计算各个粒子的权重并输出当前状态量。包 含以下几个步骤：

(31)按照表1预测每个粒子下一时刻的存在变量

表1粒子存在变量更新流程

            

(32)对存活粒子状态进行预测更新，包含新生粒子和存活粒子两种，其中新生 粒子是上个时刻为死亡状态，而当前时刻为存活状态；存活粒子则是上个时刻和 当前时刻均为存活状态。

32a)新生粒子状态向量从q_b(s_k|z_k)中采样得到，此时先验概 率密度函数与新息相关。建议分布q_b(s_k|z_k)服从观测域中z_k>γ条件区域的均匀 分布。水下运动目标的深度、距离和速度均服从均匀分布。

32b)存活粒子状态向量通过状态转移方程计算更新。

(33)计算粒子的权重，进行归一化处理，输出当前估计状态值。包含两个步骤：

33a)计算粒子的重要性权值

      ${\tilde{w}}_k^n=\left(\begin{array}{cc}\underset{i\in C_i\left(s_k^n\right)}{\Pi }\underset{j\in C_j\left(s_k^n\right)}{\Pi }l\left(z_k^{(i,j)}\right|s_k)& ,E_k^n=1\\ 0& ,E_k^n=0\end{array}\right)---\left(7\right)$             

其中，$C_i\left(s_k^n\right)=\{i_0-p,...i_0-1,i_0,i_0+1,...i_0+p\},$i₀表示最接近状态向量位置的整 数，p表示目标对邻近分辨单元的影响程度。其中表示相关处理的似 然函数。定义互谱密度函数C＝E[yy′]＝hh′+R，则似然函数可表示为：

      $l\left(s\right)=\exp (-\frac{{\phi }_j\left(s\right)}{R})---\left(8\right)$      

      ${\phi }_j\left(s\right)=\mathrm{trC}-\frac{h^{\prime }\mathrm{Ch}}{h^{\prime }h}---\left(9\right)$      

其中trC是对矩阵C进行求迹运算。

33b)对粒子重要性权值进行归一化处理：

      $w_k^n=\frac{{\tilde{w}}_k^n}{\underset{n=1}{\overset{N}{\Sigma }}{\tilde{w}}_k^n}---\left(10\right)$      

33c)输出目标的状态估计向量

步骤4：粒子重采样，并循环步骤3处理新的数据。对于常规的粒子滤波而 言，粒子退化是其主要的缺陷，即随着迭代次数的增加，粒子的权重会集中在少 数几个粒子上，而其他大多数粒子的权值很小，结果导致粒子集无法表示真实的 后验概率密度。

为了解决粒子退化的问题，引入序贯重采样技术。通过引入有效粒子数N_eff表示粒子的退化程度，其计算公式如下：

      $N_{\mathrm{eff}}=\frac{1}{\underset{i=1}{\overset{N}{\Sigma }}{\left({\omega }_k^i\right)}^2}---\left(11\right)$      

同时，引入一个阈值N_threshold对有效粒子数进行判断，当有效粒子数小于阈 值时，则进行重采样，如此可以降低算法的计算复杂度。

图1为水下弱目标跟踪系统模型，由一个垂直水听器阵列，一个声源和一个 海洋水声传播波导环境组成，其中垂直水听器阵列安装于起始位置，声源位于水 中。以阵接收平均信噪比等于-2dB为例，说明本发明方法的具体步骤：

步骤1：设定目标运动模型、波导环境参数和观测模型。波导区域平均水深 为106m，底的密度为ρ_bot＝1.7g/cm³，吸收系数为α_bot＝0.05dB/λ，声速为 c_bot＝1610m/s；接收阵为16元垂直均匀线阵，垂直放置于起始位置，其深度区 间位于15～75m，每隔5s接收一次声源的信号。声源辐射频率为400Hz，移动声 源起始距离为800m，初始深度45m，速度为v_s＝1.4m/s。

步骤2：设粒子数为5000，粒子出生概率为0.05，死亡概率为0.05，将距离 处于700-1400m，深度处于1-106m之间的网格水域作为目标检测区域，网格水 域的单位距离为1m。对第一帧获取的数据进行相关处理，并根据门限对相关值 进行采样，作为初始化粒子。

步骤3：垂直均匀线阵接收来自水下目标的声压信号，进行相关处理得到目 标状态新息。对粒子的存在状态进行预测，并利用新息对其进行更新处理，同时 计算每个粒子的归一化权重，输出估计状态量。

步骤4：当有效粒子数低于预定的门限值时，对粒子进行重要性重采样，并 转入步骤3处理下一时刻数据。

图2和图3分别表示粒子的分布情况，从图中可以看出，初始时刻粒子分布 在整个波导环境，随着跟踪的持续，粒子快速有效的收敛到目标的真实位置，表 明该方法的有效性和准确性。

图4到图7分别表示本发明方法和常规方法在深度和距离维上的跟踪精度， 从中可知，对于水下弱目标来说，本发明方法能够有效跟踪并保证足够的跟踪精 度，而常规方法则因为信噪比太低而导致跟踪发散，无法实现正常跟踪。因此本 发明方法比常规的跟踪方法具有更强的适用性。

图8表示跟踪结果目标状态存在概率，可以看出随着测量数据的增加，本发 明方法能够有效的对低信噪比目标进行稳定跟踪，具有较好的可靠性。

本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的 保护范围不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也及 于本领域技术人员根据本发明构思所能够想到的等同技术手段。