基于贝叶斯滤波算法的警戒相控阵雷达目标跟踪时间资源优化分配算法

李纪三， 侯娇， 刘溶， 任渊

(南京船舶雷达研究所， 江苏 南京 210003)

0 引言

相控阵雷达通过改变信号单元发射的电磁波相位来实现不同方向的波束指向，在计算机控制下，在方位和仰角上可达到微秒量级的波束捷变，具有扫描方式灵活、工作波形丰富、控制参数可变、工作模式可选等优点[1-2]，通过对相控阵雷达的优化控制，包括任务合理规划、资源有效调度以及工作参数优化选择等，可实现多功能多任务并能释放雷达作战潜能[3-5]。

雷达的主要任务是在各种电磁环境下探测、跟踪和识别空中及海面目标，其中跟踪任务的时间资源分配是资源调度的重要方面，跟踪任务消耗的时间资源与跟踪数据率呈正比，跟踪数据率是指单位时间内对目标的采样次数，是采样时间间隔的倒数。跟踪的数据率影响跟踪精度和稳定性，数据率越高，跟踪的稳定性和跟踪精度越高，但是消耗的时间资源也越多。跟踪任务的时间资源分配是指在保证一定精度和稳定性的前提下，使得跟踪花费的时间资源最少。根据当前跟踪的状况，自适应地确定目标的跟踪数据率是优化分配时间资源的有效途径。文献[6]介绍了基于跟踪残差自适应数据率跟踪算法，该算法通过残差与量测噪声方差间的比值对下一时刻采样周期进行调整，当目标机动时该算法可自适应地减小采样周期，但是当目标无机动性时采样周期却会不断增大。van Keuk等[7]研究了采样周期与信噪比、跟踪稳定度以及资源消耗的关系，提出了基于跟踪误差、相对跟踪精度和波束宽度的采样间隔计算公式。上面两个文献奠定了后续研究基础。Benoudnine等[8]基于交互多模型滤波算法实现了自适应采样间隔计算，下个采样间隔是各模型计算采样间隔的加权平均。文献[9]引入了目标威胁等级，作为调整数据率的参考值。文献[10]将van Keuk等[7]的自适应准则扩展到扩展卡尔曼滤波算法、不敏卡尔曼滤波算法、交互多模型滤波算法和粒子滤波算法。

以上研究把跟踪精度作为调整数据率的依据和准则，随着相控阵雷达体制的发展，和差测角波束因其具有实现简单和精度高的特点已得到广泛应用，从滤波算法角度，量测误差相对于传统机扫雷达大大减小，跟踪的主要矛盾变成跟踪稳定性，即跟踪波束能否罩住目标。

基于残差的自适应跟踪算法具有运算量小、工程实现简单等优点，该算法为了实时响应目标的机动，通常利用当前采样点的残差计算下个周期的采样间隔。在导弹攻击或者飞机逃离过程中会采取多种样式的机动方式，在蛇形机动反向机动阶段残差会下降到很小的值，此时会误判为机动结束而增加跟踪数据率，进而造成目标逃出雷达的跟踪波门，造成失跟。

相对于传统基于残差的变数据率跟踪算法把机动引起的残差看作一个确定变量，本文基于贝叶斯滤波算法，把残差看成随机变量，并用相应的概率密度加以描述。在目标匀速直线运动阶段，残差服从正态分布，通过对残差进行拉依达准则判决来抑制量测噪声起伏产生的数据率抖动问题。在机动阶段，残差服从修正的瑞利分布，对数据率采用指数平滑法进行平滑，通过在数据率增加和减小时段采用不同的平滑因子，在数据率下降时采用较大的平滑因子，即在计算下个采样间隔时本周期的残差中所占比重较大，降低噪声引起的机动误判和快速响应目标的机动，在数据率上升时采用较小的平滑因子，即历史残差所占比重大，从而减少误判机动结束所造成的目标失跟。

1 贝叶斯滤波原理

贝叶斯原理是利用已知信息来构造系统状态变量的后验概率密度，即用系统状态转移方程来计算状态的先验概率密度，再利用观测值对先验概率密度函数进行修正，从而得到状态变量的后验概率密度[11]。

设随机、离散系统S的状态空间模型为

(1)

式中：Xk为k时刻(k≥1)的系统状态向量；F(·)为系统状态转移模型；Uk为系统随机噪声；Zk为k时刻的系统观测向量；H(·)为系统观测模型；Vk为随机观测噪声。

求解状态后验分布P(Xk|Zk)是贝叶斯滤波的关键问题，可由如下贝叶斯滤波方程递推求解[12-14]：

1)预测方程(Chapman-Kolmogorov方程)

(2)

2)更新方程

(3)

若系统S是线性的，系统噪声Uk和观测噪声Vk为零均值高斯白噪声(即服从高斯分布)，其协方差矩阵记为Var(U)=Q，Var(V)=R，则基于贝叶斯滤波算法的卡尔曼滤波算法[15]可表述为

步骤1一步预测：

(4)

Pk|k-1=APk-1|k-1AT+Q，

(5)

步骤2观测更新：

Kk=Pk|k-1CT(CPk|k-1CT+R)-1，

(6)

(7)

Pk|k=Pk|k-1-KkCPk|k-1,

(8)

经考察得知，Pk|k-1的变化趋势是非负递减的，当k趋于无穷时，系统达到稳态存在极限：

(9)

式中：P为协方差矩阵；K为增益矩阵；M为预测方差矩阵。

当目标运动方程采用常速度模型时，得到稳态下的常增益滤波器为

(10)

(11)

2 常增益滤波器的残差计算

1)一个加速度为u的匀加速模型

X(k)=AX(k-1)+Gu，

(12)

2)测量方程

Y(k)=CX(k)+ω(k),

(13)

式中：Y(k)为k时刻的量测值；ω(k)为量测噪声，并满足高斯分布。

3)滤波方程

(14)

4)状态预测方程

(15)

5)量测预测方程

(16)

将滤波的误差定义为真值减去估计值，有

(17)

于是，当系统稳定之后，k时刻的误差应该与k-1时刻的误差相同，即

(18)

(19)

对于匀加速直线模型，不考虑量测噪声ω(k)下，由纯加速度u引起的位置误差为

(20)

同理可推导出由u引起的位置残差[16]为

e(k)=T2u/β.

(21)

3 自适应数据率算法

3.1 基于残差的采样周期计算

如果误差在加速度增加时要自适应地保持恒定，则由(20)式就得出采样周期应按与加速度的平方根呈反比的方式减小，文献[6]给出了周期递推的计算公式：

(22)

式中：D(k)为跟踪的数据率，是采样周期T(k)的倒数；e0(·)为归一化残差。(22)式可展开为

(23)

跟踪残差包含系统的过程噪声(也称为模型噪声)和量测噪声，过程噪声是由选择的目标运动动态方程与实际目标运动过程不符产生的，例如目标当前在做圆周运动，而用直线运动方程来表述，则会产生误差。

在目标做匀速直线运动、没有机动情况时过程噪声为0，残差主要是由量测噪声贡献的，由于量测噪声服从正态分布，e0(1)，…，e0(k)相互独立且服从正态分布N(μ,σ2)，密度函数[17]为

(24)

式中：x为跟踪残差；σ为标准差；μ为均值。

图1所示为在目标做匀速运动时按照(22)式直接计算出的数据率。从图1中可观察到，由于量测噪声的正态分布特性，噪声的起伏引起虚警(目标没有机动，误判为目标机动，从而增加数据率)。

图1 匀速直线运动的数据率Fig.1 Data rate of uniform linear motion

针对图1由观察噪声引起的数据率起伏问题，本文引入统计决策中2σ准则[18]:

P(μ-2σ

(25)

式中：P(·)为概率密度函数。由(25)式可知，如果当前归一化残差小于2σ，则认为残差是由量测噪声起伏引起的，数据率保持不变；如果当前归一化残差大于2σ，则认为残差是由目标的机动引起的，需要根据残差调整数据率。

3.2 目标机动模型

目标的机动可由以下模型进行描述：匀速模型、匀加速模型、Jerk模型、Singer模型和当前统计模型等。在这些模型中，当前统计模型较为典型，该模型假定加速度服从修正的瑞利分布，其均值是当前加速度，其均值与方差之间的关系可以用来建立机动加速度的均值和方差自适应算法。修正的瑞利密度函数[19]可描述如下。

1)当目标的当前加速度为正时，概率密度函数为

(26)

式中：umax为已知目标加速度的上限，umax>0；u的均值和方差为

(27)

(28)

2)当目标的当前加速度为负时，概率密度函数为

(29)

式中：u-max为已知目标加速度的负下限，u-max<0；u的均值和方差为

(30)

(31)

3)当目标在机动时，加速度u为修正的瑞利分布，分布的起伏特性决定了加速度引起的残差T2u/β(见(21)式)，也服从瑞利分布。如果通过(22)式直接计算数据率，则会引起数据率的剧烈振荡，不利于稳定跟踪，特别是在机动阶段的反向机动时刻，即由正负加速度相互切换时会出现残差为0的情况，此时按照公式要增加采样间隔，会引起目标的失跟，为此利用指数平滑法对数据率进行平滑。

一次指数平滑法[20]的公式为

Ds,k=ρDk+(1-ρ)Ds,k-1,

(32)

式中：ρ为平滑系数；Ds,k为时刻k的平滑值；Dk为时刻k的实际值；Ds,k-1为时刻k-1的平滑值。(22)式代入(32)式，有

(33)

(34)

从(34)式中可以看出，当前的数据率与第1次的数据率和残差序列以及平滑系数ρ有关。k时刻计算数据率时，ρ越大则当前的残差在计算中占比越大，ρ越小则以前的残差占比越大。因此可以通过实时调整ρ的取值来调节当前残差占的比重：在机动开始时刻(即残差突然变大)采用较大的值，使得系统响应迅速；在残差变小时采用较小的值，防止机动未结束、误判为结束而采用较小数据率，造成目标的失跟。

3.3 变数据率跟踪算法流程

本文算法是在直角坐标系中进行残差计算，雷达的量测在球坐标系中进行，转换到直角坐标系时存在非线性问题。若利用过程噪声对残差进行归一化，则过程噪声计算误差较大，而文献[6]介绍的算法利用0.1 n mile的常量对残差进行归一化会引起较大的误差。本文算法在直角坐标系下滤波，在量测坐标系下计算残差，利用雷达本身的量测噪声对残差进行归一化，具体实施步骤[21]如下：

步骤1对k时刻状态进行预测:

(35)

式中：xp(k)、yp(k)、zp(k)分别为x轴、y轴、z轴的位置预测值；xs(k)、ys(k)、zs(k)分别为x轴、y轴、z轴的位置平滑值；vs，x(k)、vs，y(k)、vs，z(k)分别为x轴、y轴、z轴的速度平滑值；T(k)为时刻t(k)与t(k+1)间的采样周期。

步骤2直角坐标系下的预测值转换为球坐标系：

(36)

式中：rp(k)、βp(k)、θp(k)分别为k时刻预测值在极坐标下的距离、仰角和方位。

步骤3雷达在(βp(k),θp(k))照射后得到量测值(rm(k),βm(k),θm(k))，将量测值转为直角坐标系[21]：

(37)

式中：xm(k)、ym(k)、zm(k)分别为x轴、y轴、z轴在k时刻的量测值。

步骤4在直角坐标系下对k时刻目标的位置进行滤波：

(38)

步骤5根据k时刻的量测值和预测值计算k时刻的残差，残差e(k)定义为t(k)时刻测量值与预测值之差，即

(39)

式中:er(k)、eβ(k)、eθ(k)分别为距离残差、方位残差和仰角残差。

步骤6残差归一化：

e0(k)=max(eβ(k)/σβ,eθ(k)/σθ)，

(40)

式中:σβ、σθ分别为仰角量测标准差和方位量测标准差。

步骤7根据k时刻的残差计算k+1时刻的采样周期，计算下个周期的数据率

(41)

步骤8对计算的数据率进行平滑，

D(k)=ρD(k)+(1-ρ)D(k-1)，

(42)

式中：

(43)

(43)式表明：当e0(k)>1时，要残差变大需要降数据率，此时平滑因子取值较大，能够较快地响应目标的机动，防止目标失跟；当e0(k)=1时，保持数据率不变；当e0(k)<1时，要增加数据率，此时选取的平滑因子较小。为了防止误判，即机动没有结束而增大数据率造成失跟，此时采用这种策略会在某种程度上牺牲一些时间资源，保证了目标的稳定跟踪。

4 仿真实验

为验证本文所提算法的有效性，与常规天线周期作为数据率的定数据率Cohen跟踪算法作比较，利用300次蒙特卡洛算法计算得到的均方根误差(RMSE)作为滤波性能的衡量标准，对比分析两种算法在目标作匀速直线运动和机动情况下的跟踪性能。

4.1 实验1

仿真的雷达参数为：雷达位于坐标原点(0 km，0 km)处，量测精度为距离60 m，方位0.25°，仰角0.25°.

仿真场景如图2所示，飞机从方位45°、距离120 km、高度15 km，以速度830 km/h沿直线飞向雷达，匀速直线飞行40 s后开始作蛇形机动，转弯向心加速度10g、转弯半径6.9 km，先向右转弯90°，然后向左转弯180°，再向右转弯90°，最后沿着直线以速度830 km/h匀速直线飞行40 s结束。

图2 目标飞行轨迹Fig.2 Target flight trajectory

跟踪的残差为预测值和量测值的差，残差反映了模型预测的准确性和对目标的跟踪效果，用来作为调整数据率的依据。对于相控阵雷达跟踪和搜索(TAS)，采用了单波束进行跟踪，为了提高测量精度，波束宽度设计得通常不宽。当跟踪的残差大于波束宽度的一半后就会造成失跟，虽然可以用更宽的波束进行重新捕获，但是捕获花费的时间资源是普通跟踪的数十倍，失去了资源优化调度的意义。图3和图4为两种算法进行50次仿真的分布。从图3和图4中可以看出：在0～40 s和100～140 s阶段是匀速直线运动阶段，两种算法的残差分布相当；在50～100 s阶段蛇形机动阶段，本文算法的跟踪残差要小。

图3 Cohen算法的跟踪残差Fig.3 Tracking residual of Cohen algorithm

图4 本文算法的跟踪残差Fig.4 Tracking residual of the proposed algorithm

为了更明显地对比两种算法的性能，计算了残差的方差，结果如图5所示。从图5中可以看出本文算法优于传统算法，本文算法的跟踪稳定性高。

图5 跟踪残差对比Fig.5 Comparison between tracking residuals of traditional algorithm and the proposed algorithm

由于目标主要在方位上进行机动，比较了方位上的跟踪精度，跟踪的误差分布如图6和图7所示，可见误差分布和残差的分布相似，如果采用常系数阿尔法- 贝塔滤波器，则误差和残差有固定关系，误差是残差的一半。从图6和图7中也可以看出这种趋势。

图6 Cohen算法的跟踪误差Fig.6 Tracking error of Cohen algorithm

图7 本文算法的跟踪误差Fig.7 Tracking error of the proposed algorithm

方位跟踪精度如图8所示。由图8可见，本文算法的精度稍优，传统算法随着蛇形机动的方向改变，精度发生了起伏，主要原因是反向机动时残差变小，跟踪的数据率增大，但此时机动加速度大小没有变化，采样周期的变大造成了跟踪精度的下降。

图8 方位RMSE对比Fig.8 Comparison between bearing RMSEs of traditional algorithm and the proposed algorithm

图9和图10为两种算法采样周期的分布图。由图9和图10可以看出：在反向机动的过程中由于没有采用平滑因子，50 s、70 s和100 s附近是残差最大的时刻，采样周期最小，但是采样周期变大；采用本文算法在目标的蛇形机动过程中，采样周期仍采用较低的值；本文算法全程平均采样周期为1.28 s，传统算法平均采样周期为1.79 s，精度和稳定性增加的代价是时间资源增多。

图9 Cohen算法的跟踪采样周期Fig.9 Sampling period of Cohen algorithm

图10 本文算法的跟踪采样周期Fig.10 Sampling period of the proposed algorithm

4.2 实验2

目标的初始位置为：x轴25 km，y轴50 km，指向为90°，初始速度为0 m/s，加速度为0g. 雷达坐标在(0 km，0 km)位置。目标的运动情况为：经过30 s，以10 m/s2纵向加速度达到300 m/s的速度；从40 s开始持续40 s，以-20 m/s2的横向加速度进行转弯；从200 s开始持续15 s，以20 m/s2的横向加速度进行转弯，调整姿态，飞向雷达；从300 s开始持续30 s，以30 m/s2的横向加速度进行转180°，转弯半径约为2 km；从400 s开始持续10 s，以-90 m/s2的横向加速度进行180°转弯，转弯半径约为1 km；从450 s开始持续30 s，以10 m/s2的纵向加速度加速到600 m/s；从500 s开始持续22.5 s，以-90 m/s2的横向加速度进行180°转弯，转弯半径约为2 km；从620 s开始持续60 s，以-10 m/s2的纵向加速度，减速停止，目标运动轨迹如图11所示。

图11 目标飞行轨迹Fig.11 Targer flight trajectory

实验2相对于实验1机动更丰富，更能反映雷达的跟踪情况。图12所示为跟踪的残差对比。由图12可见，在弱机动阶段两种算法的残差水平基本相当，在突然机动阶段由于对目标的运动描述模型与目标的实际模型相差很大，因此按照理论模型对目标的预测位置与目标的实际位置相差很大，残差很大，如图12中的223 s处和600 s处。残差超过波束宽度的一半就容易造成目标失跟，图12中传统算法的残差方位最大值达到1.8°，本文算法最大值为1.2°.

图12 残差RMSE对比Fig.12 Comparison between residual RMSEs of traditional algorithm and the proposed algorithm

图13为两种算法方位的跟踪误差对比，可见两种算法跟踪误差曲线与残差一致，只是幅度不同，在机动时误差大，在弱机动时误差小。图14为两种算法平均数据率，可见本文算法采用平滑因子，在采样周期减少阶段采用较大的平滑因子，因此采用周期不会突然增大，保证了目标的稳定跟踪。

图13 方位RMSE对比Fig.13 Comparison between bearing RMSEs of traditional algorithm and the proposed algorithm

图14 跟踪数据率对比Fig.14 Comparison between tracking data rates of traditional algorithm and the proposed algorithm

5 结论

本文针对警戒相控阵雷达目标跟踪的时间资源优化分配问题，基于贝叶斯滤波算法提出一种基于残差的跟踪时间资源自适应分配算法，在直角坐标下进行滤波在球坐标系下计算残差，并利用量测精度对残差进行归一化，引进平滑因子并在采样间隔上升和下降时期采用不同取值。得到主要结论如下：

1)本文算法比传统算法计算归一化残差简单，避免了坐标转换的非线性问题。

2)平滑因子能够有效抑制目标在蛇形机动反向机动阶段因残差减小导致的采样间隔扩大引起的失跟问题。