基于leg-by-leg机动的两级采样被动跟踪方法

奚 畅 蔡志明 袁 骏

(海军工程大学电子工程学院 武汉 430000)

1 引言

目标被动跟踪又称目标运动分析(Target Motion Analysis,TMA)，指利用被动观测的信息对目标运动状态在时间上进行连续估计的过程。线谱信号的方位和频率是被动声呐观测的重要参数，利用低频分析与记录(Low frequency analysis recording,Lofar)谱的方位-频率TMA[1,2]对于水下态势感知具有重要意义。检测前跟踪(Track-Before-Detecting,TBD)技术不对单帧数据做门限处理，利用多帧数据积累能量进行检测判决，可提高目标的发现概率。基于粒子滤波的TBD方法[3,4]通过引入目标运动模型和传感器观测模型，完整地体现了跟踪的思想，是当前弱目标TBD的研究热点。

但是，粒子滤波TBD算法在方位-频率观测情况下的应用尚存在障碍。同时观测方位和频率时目标状态向量维数较高，而粒子滤波TBD算法在高维状态空间的采样效率较低[5–7]，为满足一定的检测和估计要求，算法所需的粒子数量随采样维数呈指数级增长[8]，带来的计算量和存储量都是灾难性的。

文献[9]根据状态向量中各分量对量测有无直接影响，把粒子的高维状态采样转换为两个低维采样，缺陷是在第1级仅对位置和强度采样，目标位置发生变化时无法沿目标轨迹累积能量，未能体现TBD算法“时间换信噪比”的思想。文献[10]认为对量测有直接影响的状态向量分量的估计会较早接近真值，根据状态向量是否可测分别从后验状态分布和先验分布中采样，提出局部搜索采样的方法以提高粒子采样效率。Rao-Blackwellized粒子滤波[11]将粒子状态变量划分为线性状态变量和非线性状态变量，通过Kalman滤波方法估计线性状态变量，通过粒子滤波估计非线性状态变量，实现粒子滤波状态向量降维。标准的Rao-Blackwellized粒子滤波实现中要求对每个粒子运行一个Kalman滤波，虽然状态向量降维减少了所需的粒子数，但计算量的减少并不显著[12]。文献[13]在高维状态空间假设下比较了粒子滤波、马尔可夫蒙特卡罗粒子滤波、混合粒子滤波3种算法的性能，但没有给出提升算法性能的方法。一些学者尝试从粒子滤波原理入手，通过数值逼近等方法解决高维采样问题[14,15]，理论较为复杂且无法直接应用于目标跟踪的场景。

leg-by-leg机动[16,17]是一种易于实施的机动模式，适用于船艇等机动性欠佳的观测载体。本文针对观测站leg-by-leg机动模式下利用方位-频率信息的粒子滤波检测前跟踪算法，一方面利用leg-byleg机动可观测性特点，在前、后直行段分别建立极坐标系和直角坐标系下的目标状态模型，提出将极坐标系下的目标状态向量映射至直角坐标系的方法，通过两级采样将一个高维采样问题变为两个低维采样问题，从而改善粒子的采样效率；另一方面根据粒子的空间分布特征，自适应地调整过程噪声协方差矩阵，从而改善滤波的收敛性，避免陷入局部最优。仿真结果表明，所提方法可以有效地增大滤波收敛率、减小目标距离估计误差、缩短收敛时间，海试数据处理结果进一步验证了所提方法的可行性和有效性。

2 系统模型

2.1 目标状态模型

假设目标具备线谱特征，目标与观测站在同一水平面内做匀速直线运动，在直角坐标系下建立k时刻的目标状态向量为

其中，xk和yk分别为目标在x轴方向和y轴方向的坐标，分别为目标在x轴方向和y轴方向的速度，fk为目标线谱固有频率，snrk为目标线谱信噪比。

目标状态转移方程为

其中，F与Q分别为目标匀速直线运动模型的状态转移矩阵以及过程噪声协方差矩阵

T为观测时间间隔，q1,q2,q3分别为目标运动状态、线谱频率、线谱信噪比的过程噪声级。

2.2 量测模型

功率P(xk)与 状态向量中信噪比项snrk的关系为

其中，fi为第i个 频率观测单元对应的频率值，βj为第j个方位观测单元对应的方位值。观测站与目标间的相对运动会导致线谱发生多普勒频移，fr,k为目标线谱固有频率与多普勒频移之和，βk为目标方位。

3 leg-by-leg机动可观测性分析

leg-by-leg机动由两个不同运动方向的匀速直行段及连接它们的圆弧段组成，又称转向机动、折线机动、两段直航式机动，是一种易于实施的机动模式，适用于船艇等机动性欠佳的观测载体。按照时间先后顺序，两个匀速直行段分别称为leg1段和leg2段，示意图如图1所示。

图1 leg-by-leg机动模式示意图

被动量测信息是目标运动参数的不完全描述，且状态空间模型的非线性程度较高，因此需分析系统可观测性[18]。目标运动参数可观测的必要条件是测量信息维度大于目标运动方程阶数[19]，在目标匀速直线运动(运动方程阶数为1)的假设下，为保证目标运动参数可观测，测量信息的维度应大于1。增大观测信息维度的方法包括：增加观测量数目，如同时观测目标的方位、频率、相位差变化率等信息；增加观测站数目，如采用双基地、多基地观测；增大观测站运动方程阶数，如观测站进行圆周机动、抛物线机动。

方位-频率观测情况下，当观测站处于leg-byleg机动模式中转向前的匀速直行段(leg1段)时，测量信息维度为2，在观测站和目标存在相对运动的情况下对于匀速直线运动的目标具备可观测性[20]。但是，在目标运动分析问题中，除了讨论是否可观测，还应讨论可观测性的强弱，这关系到目标状态估计的收敛率、收敛速度、估计精度等性能[21]。目标运动状态的估计依靠观测量的变化，变化幅度越大越有利于估计。在水下目标探测场景下，由于观测站和目标航速较低且相距较远，因此方位变化幅度较小；由于水下目标线谱固有频率较低，因此多普勒效应不显著，频率变化幅度较小。特别是在低信噪比情况下，目标方位变化和多普勒频移的检测更加困难。因此，leg1段虽然具备可观测性，但可观测性较弱。

当观测站处于leg-by-leg机动模式中转向后的匀速直行段(leg2段)时，观测站转向前后等效于两个不同位置的观测站同时对目标进行观测[22]，通过增加观测站数量的方式进一步增大了观测信息维度，构成观测维数冗余。转向后目标方位变化率会发生较大幅度的变化，且相对运动速度的变化会导致多普勒频移变化，观测量较大幅度的变化有利于目标运动状态估计，因此leg2段具有较强的可观测性。

由上述分析可知，在方位-频率观测情况下，对于匀速直线运动的目标，leg1段的可观测性较弱，很难准确估计目标的运动状态，leg2段的可观测性较强，在目标运动分析时具有较好的收敛率、收敛速度和估计精度。

4 两级采样方法

两级采样方法即在leg1段对极坐标系下的目标状态信息进行采样，在转向前将极坐标下的粒子状态向量映射至直角坐标下，在leg2段对完整的目标状态向量进行采样，从而将一个高维采样问题变为两个低维采样问题。

leg1段可观测性较弱，很难准确估计完整的目标运动状态，Lofar谱中的方位信息是极坐标系下目标运动参数的完全描述，因此可在leg1段构建极坐标系下的目标状态模型。对于目标与观测站在同一水平面内做匀速直线运动的情况，由于被动声呐探测场景下观测站和目标航速较低、相距较远，可近似认为目标位变率恒定；由于目标方位变化幅度较小且水下目标线谱频率较低，可近似认为相对运动导致的线谱多普勒频移恒定。设为目标的方位变化率，建立极坐标系下的目标状态向量为

设q4为目标位变率的过程噪声级，状态转移矩阵及过程噪声协方差矩阵分别定义为

图2 相对速度示意图

式(27)说明，对于一定的目标方位和相对运动切向速度，不论如何采样，粒子始终分布在一条与目标方位线平行的直线上。为避免陷入局部最优，的采样原则是使粒子在与目标方位线平行的直线上始终保持均匀分布，即粒子相对运动的法向速度在有效范围内服从均匀分布，法向速度表示为

α在有效范围内服从均匀分布，其概率密度函数为

由连续型随机变量函数分布定理[23]可得的概率密度函数

5 自适应过程噪声协方差矩阵

对于含有角度测量的目标运动分析，一般采用极坐标系、对数极坐标系、修正极坐标系[24,25]下的目标状态向量建模方法，以获得较高的滤波稳定性。本文leg2段在直角坐标系下对目标状态向量建模，估计结果可以直观地反映目标运动情况，但是稳定性不够，容易导致滤波发散。导致此问题的原因是目标状态更新时采用固定的过程噪声协方差矩阵，重采样后的粒子在2维空间内的分布近似为矩形(如图3(a)所示)。粒子相对于观测站方位线的切向偏移一定时，距离较远的粒子方位误差较小，可以获得更大的权重，导致粒子逐渐向远离观测站的方向聚集，使系统陷入局部最优，最终造成滤波发散。

针对此问题，以各粒子与观测站间距离作为其空间分布特征，利用此特征自适应地调整过程噪声协方差矩阵，令式(4)中的目标运动状态过程噪声级q1与粒子距离成正比，使重采样后的粒子在2维空间内的分布近似为以目标方位线为轴的扇形(如图3(b)所示)。由此可保证在滤波收敛前相距观测站不同距离带内的粒子数近似一致，提高滤波稳定性。

图3 粒子分布示意图

设目标距离上限rmax对应的目标运动状态过程噪声级为qmax，对于距离为的粒子，其过程噪声级为

6 算法实施流程

两级采样被动跟踪方法的实施步骤如下。步骤1，在极坐标系下，基于被动声纳Lofar谱，利用粒子滤波TBD方法跟踪目标。步骤2，将极坐标系下的粒子集映射至直角坐标系下。步骤3，在直角坐标系下，基于被动声纳Lofar谱，利用粒子滤波TBD方法跟踪目标。

针对观测站leg-by-leg机动模式的情况，须在观测站转向前完成步骤2，即步骤1、步骤2须在第1个直行段完成。

7 仿真试验与结果分析

7.1 仿真场景设置

为了验证前述假设的合理性和所提算法的有效性，并定量地分析算法性能，设置一种典型的被动声呐水下目标跟踪场景，采用仿真方法进行试验研究。下面说明仿真场景的设置流程。

第1步假设相对运动态势。设观测站在搜索目标时保持8 kn 航速、0°航向的匀速直线运动，以发现目标时刻为初始时刻，观测站转入leg-by-leg机动模式，在原工况下继续直行 600 s后，向左以300 m 转向半径进行转向，将航向调整至2 70°后直行 6 00 s。设具备线谱特征的目标保持匀速直线运动，目标航速 4 kn，目标航向4 5°，目标线谱固有频率 175 Hz，线谱信噪比1 2 dB，初始时刻目标与观测站距离1 0 km，目标方位1 20°。

上述假设的观测站、目标机动模式与实际情况相符，且设置的观测站航速、发现目标距离、目标航速、目标线谱固有频率及信噪比均为常规值，因此认为这是一种较为典型的水下目标跟踪相对运动态势，观测站与目标的运动轨迹如图4所示。

图4 相对运动态势

第2步基于假设的相对运动态势，计算目标方位的变化情况，以及受多普勒效应影响后到达观测站位置的目标线谱频率变化情况，结果如图5所示。

由图5可知，对于leg1段(0～600 s)，目标方位由 1 20°逐渐增大至1 24.6°，目标方位变化速度较为稳定，目标频率变化幅度小于0.02 Hz。由此可证，第4节构建低维目标状态向量时假设目标位变率恒定且线谱多普勒频移恒定是合理的。对于leg2段(720～1320 s)，观测站的转向造成相对运动态势变化，进而导致目标的方位及线谱频率均发生了较大程度的变化，此时系统具有较强的可观测性，利于目标运动参数的估计，与第3节可观测性分析的结论一致。

图5 目标方位及频率变化情况

第3步仿真各观测时刻的Lofar谱数据。设被动声呐Lofar谱观测的方位范围为[0,360]°，方位单元间隔为0.2°，频率范围为[150,200] Hz，频率单元间隔为0.1 Hz，观测区域包括1 800×500个方位-频率观测单元，观测时间间隔1 0 s。根据上述设置的Lofar谱显示特性，结合第2步计算的目标方位及线谱频率变化情况，利用2.2节量测模型，计算各观测时刻的Lofar谱数据。初始时刻Lofar谱仿真值如图6所示，图6(a)中目标方位、频率对应的观测单元位于红色圆中，将目标所在的局部区域放大显示如图6(b)。

图6 初始时刻Lofar谱仿真结果

7.2 可观测性验证

对于7.1节设置的仿真场景，计算机动过程中的目标距离估计误差的Cramer-Rao下界(Cramer-Rao Lower Bound,CRLB)，结果如图7所示。

图7 目标距离估计误差CRLB

由图7分析可知，随着观测时间增长，累计观测量增多，目标距离估计误差CRLB逐渐下降。在leg1段(0～600 s)CRLB并非无穷大，认为对目标具备可观测性，但leg1段最终时刻的CRLB依然大于73.6 km，即可观测性很弱，无法实现对目标运动状态的估计。第 6 00 s开始转向后，CRLB迅速减小，第 7 50 s CRLB减小至1.0 km，第9 60 s CRLB减小至0.1 km，说明转向对于系统可观测性能的提升是明显的，目标在leg2具备较强的可观测性，利于目标运动状态估计。上述仿真结果与第3节可观测性分析的结论一致。

7.3 目标距离估计性能验证

对于7.1节设置的仿真场景，分别用传统的高维采样粒子滤波TBD方法以及本文所提两级采样粒子滤波TBD方法对数据进行处理，对于粒子数为1000～10000间隔1000变化的情况，各自进行50次蒙特卡罗仿真。对于两级采样方法，将极坐标系下的粒子集映射至直角坐标系时，式(18)中初始分布方差设置为极坐标系下粒子集状态的方差乘以0.5。

最终时刻的距离估计误差小于真实距离的10%时判定滤波收敛，可计算各粒子数情况下的滤波收敛率。在滤波收敛的情况下，将各次蒙特卡罗仿真的最终时刻距离估计误差进行平均，得到各粒子数情况下的平均距离估计误差。在滤波收敛的情况下，将目标距离估计误差随时间变化趋势作为滤波收敛情况，对各粒子数情况下各次蒙特卡罗仿真的滤波收敛情况进行平均，可得高维采样方法和两级采样方法的平均滤波收敛情况。

需要说明的是，高维采样方法在各观测时刻均对式(1)所示的6维目标状态向量进行采样，因此在初始时刻即可估计目标位置进而计算距离估计误差；两级采样方法在第1级采样时并不包含目标位置信息，在第2级采样时才可估计目标位置，因此距离估计误差在第2级采样后才可计算。各粒子数情况的收敛率如图8所示，各粒子数情况的平均距离估计误差如图9所示，传统方法和所提方法的距离估计误差变化情况如图10所示。

由图8分析可知，两级采样方法的滤波收敛率接近100%，高维采样方法的滤波收敛率随着粒子数的增加逐渐增大，所提方法相较于传统方法的滤波收敛率增大均值约47.6%。

图8 各粒子数情况的滤波收敛率

由图9分析可知，随着粒子数的增加，两级采样方法的距离估计误差基本呈单调下降趋势，所提方法相较于传统方法的距离估计误差有显著减小，在不同粒子数情况下的误差减小均值约为329 m。

图9 各粒子数情况的距离估计误差

由图10分析可知，高维采样方法和两级采样方法首次计算的距离估计误差较为接近，这是由于粒子初始化时依据同样的目标距离分布特性，初次采样时刻的距离估计值近似为目标距离期望值。在leg1段(0～600 s)，高维采样方法的距离估计误差逐渐增大，原因是目标状态更新时采用固定的过程噪声协方差矩阵，导致粒子逐渐向远离观测站的方向聚集；两级采样方法由于采用了第5节所示的自适应过程噪声协方差矩阵，距离估计误差并未增大，具备一定的滤波稳定性。第6 00 s开始转向后，两级采样方法以较短时间实现滤波收敛，高维采样方法在转向时刻的距离估计误差较大且收敛速度较慢导致收敛时间较长，所提方法相较于传统方法的滤波收敛时间缩短均值约450 s。

图10 距离估计误差随时间变化情况

由上述分析可知，本文所提两级采样方法可增大滤波收敛率、减小距离估计误差、缩短滤波收敛时间，算法有效性得证。以滤波收敛率和目标距离估计误差为评价指标，两级采样方法在粒子数为1000情况下的性能优于高维采样方法在粒子数为10000情况下的性能，说明本文提出方法与传统方法相比可将粒子采样效率提高至少10倍。

8 海试数据运用

海试数据来源于某次综合性水声试验，探测装备为被动拖线阵声呐，目标为具备线谱特征的合作声源。数据时长 1 800 s，初始时刻目标相对观测站的方位为352°，目标与观测站间距离为9km。观测站以7 kn航速、265°航向直行620s后，以固定转向半径进行逆时针转向，转向段时长580s，完成转向后的航向角为90°，以5.3kn航速直行600 s。leg1为转向前的620s，leg2为转向后的600s。合作目标航速稳定在7.8 kn 左右，初始时刻至第3 00 s目标航向为250°，第3 00 ～8 00 s 目标航向约为1 83°，第8 00 s 至最终时刻目标航向为2 42°。观测站和目标航迹如图11，经、纬度坐标刻度的整数部分用字母替代显示。

图11 观测站和目标航迹

初始时刻Lofar谱如图12所示，图12(a)中目标方位、频率对应的观测单元位于红色圆中，将目标所在的局部区域放大显示如图12(b)。

由于被动拖线阵声呐存在左右舷模糊现象，因此图12(a)以观测站航向为轴左右对称。由图12(b)可知，目标方位、频率对应的观测单元及其相邻方位观测单元能量较强，原因是单帧处理利用的时间长度内目标方位发生变化。利用本文提出的两级采样粒子滤波TBD方法对数据进行处理，设粒子数为10000，进行100次蒙特卡罗试验，将各次计算最终时刻的目标距离估计误差从小到大排列，结果如图13所示。

图12 初始时刻Lofar谱实测结果

图13 目标距离估计误差

根据目标和观测站经纬度可计算最终时刻目标的真实距离为9.9 km。假设最终时刻的距离估计误差小于真实距离的10%时判定滤波收敛，由图13可知滤波收敛率为92%。在滤波收敛的情况下，将各次蒙特卡罗计算的最终时刻距离估计误差进行平均，可得平均距离估计绝对误差为425.6 m，平均相对误差为4.3%。计算结果可基本满足实际探测需求，所提算法可行性和有效性得证。

基于海试数据计算的滤波收敛率略小于仿真计算结果，一方面原因是系统模型假设目标信噪比恒定，而海洋环境的复杂性导致目标能量较大幅度变化；另一方面原因是系统模型假设目标保持匀速直线运动，而利用的这段数据中目标进行了两次转向机动。由此也得到一个初步印象，该方法对目标信噪比和目标状态模型具有一定的稳健性。

9 结论

本文针对被动声呐方位-频率观测情况下粒子滤波检测前跟踪算法中高维采样效率低的问题，利用leg-by-leg机动可观测性特点，在前、后直行段分别建立极坐标系和直角坐标系下的目标状态模型，通过两级采样将一个高维采样问题变为两个低维采样问题，从而改善粒子的采样效率。为避免粒子滤波陷入局部最优，提出根据各粒子的空间分布特征自适应地调整过程噪声协方差矩阵，从而改善滤波的收敛性。利用仿真试验和海试数据对算法的可行性和有效性进行了验证。仿真结果表明，对于典型的水下目标跟踪场景，所提方法与传统方法相比，可使滤波收敛率增大约47.6%，距离估计误差减小约329 m，滤波收敛时间缩短约450 s。

本文基于粒子滤波检测前跟踪算法进行讨论，所提方法同样适用于先检测再跟踪的情况。对于水下目标线谱频率极低导致转向前后多普勒频移难以被观测的情况，利用方位-频率测量的被动跟踪退化为基于线谱特征的纯方位特征辅助跟踪。