李淑慧，邓志红，冯肖雪，潘 峰，2

（1.北京理工大学自动化学院，北京 100081；2.昆明北理工产业技术研究院有限公司，云南昆明 650000）

1 引言

机载雷达作为一种对地观测的重要传感器，具有全天候、全天时、高分辨的特点，能实现目标的探测、识别和估计等功能，因此受到军事领域的广泛关注.但是，机载雷达常常工作于复杂多变的恶劣环境，使得目标跟踪面临的不确定性问题异常突出：（1）机载雷达一般处于平视或下视工作状态，必然受到比地基雷达更为严重的地（海）杂波［1］.雷达杂波的脉冲干扰会使机载雷达系统的量测噪声呈现复杂、未知的非高斯长拖尾特性［2，3］.（2）实际作战环境中，机载雷达跟踪的敌方目标（如坦克、战车、无人机等）往往具有较强的机动性.敌方目标的强机动将会使目标跟踪模型中的状态量发生突变，进而诱导未知的非高斯长拖尾系统噪声［4］.（3）载机的运动导致其杂波分布范围广、强度大，尤其在城市和山区地带，副瓣杂波通常会淹没落入杂波区的目标，使机载雷达无法检测到目标，出现随机的目标航迹丢失（即量测丢失）现象［5］.总之，由于复杂的作战环境以及目标机动特性带来的不确定性因素给机载雷达目标跟踪带来了巨大的挑战.因此，研究含量测丢失和未知长拖尾非高斯噪声的机载雷达目标跟踪问题具有十分重要的意义.

目前，已有大量的文献研究含未知噪声的状态估计问题［6］.例如，Ardeshiri 等［6］针对过程噪声协方差和量测噪声协方差均未知的状态估计问题，采用变分贝叶斯方法设计了近似贝叶斯平滑器.然而，该算法只适用于高斯噪声场景.为此，一些学者也针对含未知长拖尾非高斯噪声的非线性系统状态估计问题开展了研究［7～10］.赵俊波等［8］针对含未知过程噪声协方差和量测噪声协方差的非线系统，设计了鲁棒的无迹卡尔曼估计器.黄玉龙等［10］针对含长拖尾过程和测量噪声的非线性系统，提出了鲁棒的状态平滑器.该算法凭借变分推理逼近真实后验PDF（Probability Density Function）的优势，联合估计出长拖尾非高斯噪声的自由度、尺度矩阵和系统状态.上述算法均不适用量测丢失的场景.

近年来，研究人员也针对含量测丢失的非线性系统状态估计问题进行研究［11～14］.马柯茂等［11］针对含高斯噪声的非线性系统，采用量测丢失补偿的方法设计了用于量测丢失的状态估计方法.然而，该算法难以应用于长拖尾非高斯噪声场景.江顺等［12］针对含延时和随机量测丢失的网络化非线性系统，设计了故障估计器.虽然该算法适用于非高斯噪声，但它要求噪声协方差有界且已知，量测丢失概率恒定已知，且状态估计的精度较低.上述算法只适用于噪声已知的场景，尚不能处理量测丢失概率和噪声均未知的情况.

由于机动目标跟踪场景复杂，未知的长拖尾非高斯噪声和量测丢失不可避免地同时存在于真实系统中.现有的算法均无法解决同时含量测丢失和未知长拖尾非高斯噪声的状态估计问题.为此，本文研究强杂波背景下含量测丢失的机动目标跟踪问题.

2 问题描述

考虑到复杂的作战环境、强杂波以及目标强机动等因素的影响，本文研究了如下的机载雷达目标跟踪系统［4，15，16］：

其中，xk∈Rn是k时刻目标的状态，yk∈Rm是笛卡尔坐标系下k时刻雷达接收到的测量，fk-1(·)和hk(·)分别为已知的函数.表示直接利用全球定位系统（Global Positioning System，GPS）得到机载平台的位置.随机变量ℓk服从概率未知的伯努利随机分布，表示量测数据的丢失情况.假设初始状态向量x0服从均值为协方差为P0|0的高斯分布，且与系统噪声和量测噪声互不相关.

由于闪烁噪声服从长拖尾的非高斯分布［15］，同时考虑到目标的闪烁、杂波、强电磁干扰、以及目标强机动等因素的影响，本文采用闪烁噪声来模拟未知的长拖尾过程噪声∈Rn和量测噪声∈Rm.此外，由于学生t分布可以视为无穷多个高斯分布与长拖尾Gamma分布的混合，具有长拖尾非高斯特性，常将闪烁噪声建模为学生t分布［17］.于是分别写成如下等式：

由于长拖尾的非高斯噪声未知，常常无法获得准确的尺度矩阵和自由度参数.假设尺度矩阵Q，R的先验分布分别服从逆Wishart分布［10］，即

其中，t0，u0分别表示自由度参数，T0，U0分别为逆尺度参数.未知的自由度参数ω，ν分别服从Gamma先验分布，即

此外，x0，Q，R，ω，v的联合先验概率分布服从如下等式：

此外，由于地（海）杂波以及载机运动的影响，副瓣杂波通常会淹没落入杂波区的目标，使机载雷达无法检测的目标，出现随机的量测丢失现象.当ℓk=1 时，表示k时刻雷达检测到目标.当ℓk=0 时，表示k时刻雷达丢失目标.随机变量ℓk满足如下PDF：

其中，e0∈[0，1]和e1∈[0，1]为Beta 分布的先验形状参数.

本文的主要目的在于针对含量测丢失和未知长拖尾非高斯噪声的机载雷达目标跟踪系统式（1），利用变分贝叶斯方法设计出强杂波背景下含量测丢失的机载雷达机动目标跟踪算法，进而估计出目标的状态.

3 基于变分贝叶斯推理的机载雷达机动目标跟踪算法

3.1 近分层概率图模型

本节的目的是设计分层概率图模型，为后续变分贝叶斯推理做准备.当量测丢失概率已知时，似然概率密度p(yk|xk)为

由于式（12）中似然PDFp(yk|xk)为加和的形式，很难直接利用变分贝叶斯方法估计出未知的状态.于是，将伯努利随机变量ℓk的概率转化为乘积形式的概率质量函数，

于是，似然PDFp(yk|xk)可写为

根据式（1）和式（2）得

根据式（14）和式（15）得

于是，式（4）～式（7）、式（11）、式（13）、式（16）和式（17）共同构成了如图1所示的分层概率图模型.

图1 分层概率图模型

3.2 基于变分贝叶斯方法的机载雷达机动目标跟踪算法设计

因为式（18）表达形式复杂，所以联合PDFp(Θ1：K，y1：K)不存在解析解.本文利用平均场理论［18］和变分贝叶斯方法，通过寻找一个解析、易分解的近似PDF来逼近真实的平滑后验PDF，即

令θk为Θ1：K中的任意元素，即θk∈Θ1：K.式（19）中的q(θk)是θk的近似后验PDF.通过最小化近似平滑后验PDFq(Θ1：K) 与真实后验 PDFp(Θ1：K|y1：K) 之间的Kullback-Leibler散度就可以获得q(θk)的表达形式［10］，即

其中，E(·) 表示期望，cθ是与θk无关的常数.由于q(λ1：K)的变分参数相互耦合，不能直接求解式（20），而需要利用固定点迭代法求解式（20）.固定点迭代法在迭代更新中的任意PDF 时，需要使用其他PDF 上一次的迭代值［18］.

3.3 近似后验PDF的计算

为了计算近似的平滑后验PDF，对式（18）取对数，即log(p(Θ1：K，y1：K)).利用式（20）可依次计算各个近似后验PDFq(θk)的变分参数.

（1）当Θ1：K中的元素为x0：K时，利用式（18）和式（20），q(i+1)(x0：K)计算如下：

根据式（22）可知，q(i+1)(x0：K)等同于具有修正似然PDFN(yk；hk(xk)，)的非线性后验平滑PDF.于是，将含有量测丢失的非线性系统就转化为不含量测丢失等价的非线性高斯系统.基于高斯平滑器设计原理［19］，状态的更新过程主要包括前向传播和后向传播.因此，可得如下等式：

（a）前向传播

其中，tr(·)表示迹操作.

令θ=ℓ1：K，利用式（18）得

根据式（38）和式（40），q(i+1)(ℓ1：K)可以更新为伯努利分布.随机变量ℓk取0和1时的概率分别为

其中，

（5）当Θ1：K中的元素为ν时，利用式（18）得

其中，

接下来，计算求解近似后验PDF 时所需的期望.根据Gamma分布、伯努利分布和Beta分布的性质可得

3.4 基于变分贝叶斯推理的机动目标跟踪算法

本文所提用于量测丢失的机动目标跟踪算法总结如算法1所示.

4 仿真实验验证

以机载雷达跟踪空中目标为例，验证本文所提RVBSD 算法的有效性.参考文献［15，16，20］设置了本文的仿真场景.假设空中目标做恒速转弯模型运动，即

4.1 算法估计性能验证

本小节设置如下形式的目标检测概率：

其中，K=200表示总的仿真步长.为了验证所提RVBSD 算法的有效性，本文分别针对上述5 种情况，采用如下的算法进行对比：标准的容积卡尔曼平滑器（Cubature Kalman Smoother，CKS）［21］，基于量测丢失的容积卡尔曼平滑器（Cubature Kalman Smoother for measurement Drop，CKSD）［11］，基于未知过程和量测噪声协方差的变分贝叶斯容积卡尔曼平滑器（Variational Bayesian based Cubature Kalman Smoother with unknown parameters QR，VBCKS_QR）［6］，基于未知自由度参数的变分贝叶斯容积卡尔曼平滑器（Variational Bayesian based Cubature Kalman Smoother with unknown parameters wv，VBCKS_wv）［10］，基于噪声参数未知的变分贝叶斯平滑器（Variational Bayesian based Cubature Kalman Smoother with unknown parameters QRwv，VBCKS_QRwv）［10］.CKS 和CKSD 方法的噪声协方差分别设为VBCKS_QR 方法的初始噪声参数为t0=6，T0=u0=4，U0=VBCKS_QRwv 方法的初始参数为s0=r0=g0=3 且h0=3.本文所提RVBSD 算法初始参数值分别为r0=h0=g0=s0=3，T0=U0=t0=6，u0=4，e1=e0=0.5.根据高斯分布N(x0，P0|0)随机选择初始的状态估计蒙特卡洛迭代次数和变分迭代次数N分别为500和10.为了评价上述算法的性能，选用均方根误差（Root Mean Square Errors，RMSEs）、平均均方根误差（Average Root Mean Square Errors，ARMSEs）、绝对误差（Absolute Value Biases，AVBs）以及平均绝对误差（Averaged Absolute Value Biases，AAVBs）作为性能评价指标［22］.上述所有的状态估计算法都在Matlab2020a 仿真平台上运行，选用的计算机型号为：Intel Core i7-9750H CPU@2.60 GHz.

图2 展示了采用上述6 种算法所得的目标轨迹跟踪示意图.图3 分别描述了上述6 种算法的位置RMSEs、速度RMSEs 以及加速度RMSEs 结果图.图4 是与之相对应的位置AVBs、速度AVBs 以及加速度AVBs 结果图.表1 和表2 分别是基于500 次蒙特卡洛仿真得出的位置ARMSEs、速度ARMSEs、加速度ARMSEs、位置AAVBs、速度AAVBs 以及加速度AAVBs.此外，表1 还给出了上述6种算法的计算代价.

图2 雷达目标跟踪轨迹

从图2可以看出，本文所提RVBSD算法估计的轨迹更贴近真实的目标运动轨迹.此外，从图3、图4可以看出，本文所提的RVBSD算法具有较小的位置RMSEs、速度RMSEs、加速度RMSEs、位置AVBs、速度AVBs以及加速度AVBs.需要特别指出的是，当采用标准的CKS算法估计含量测丢失和长拖尾非高斯噪声的系统状态时，不可避免地会发生滤波发散.滤波发散主要的原因在于标准的CKS 算法只适用于高斯系统且对丢失的量测未做任何的处理.虽然基于量测丢失补偿的CKSD算法采用了量测预测补偿的策略，但是由于将系统噪声等价为高斯噪声，进行状态估计时并不能很好地刻画非高斯噪声的长拖尾特性.与CKS 算法和CKSD 算法相比，剩余的三种对比算法均是将系统噪声建模为学生t分布.在某种程度上，它们能刻画非高斯噪声的长拖尾特性.但是VBCKS_QR 算法和VBCKS_wv算法分别针对未知学生t分布的尺度参数和自由度参数进行估计，不适用于二者均未知的场景.虽然VBCKS_QRwv 算法可以联合估计未知学生t 分布的尺度和自由度参数，但与VBCKS_QR算法和VBCKS_wv 算法类似，未考虑量测丢失的情况.而本文所提的RVBSD 算法不仅能联合估计未知学生t分布中的自由度参数和尺度参数，还考虑了量测丢失对状态估计的影响.综上所述，本文所提的RVBSD算法能更好地跟踪目标的轨迹.

图3 状态估计RMSEs对比

图4 状态估计AVBs对比

此外，如表1 和表2 所示，本文所提RVBSD 算法有较小的位置ARMSEs、速度ARMSEs、加速度ARMSEs、位置AAVBs、速度AAVBs 以及加速度AAVBs.这也验证了本文所提算法的有效性.虽然本文所提的RVBSD算法与其他5种算法相比耗时相对较长，但是为了获得较高的估计精度，牺牲一定的时间成本也是允许的.总体而言，本文所提RVBSD算法是有效的.

表1 状态估计的RMSEs对比

表2 状态估计的AAVBs对比

4.2 不同检测概率下算法估计性能验证

为了验证不同检测概率下本文算法的估计性能，除情形1外，还采用如下几种形式的目标检测概率进行对比分析.

鉴于4.1 节已经验证本文算法的优势，为节省篇幅，本小节省略了其他算法的结果分析.除了目标检测概率，其他参数的设置均同4.1节.表3和表4分别给出了不同目标检测情形下位置ARMSEs、速度ARMSEs、加速度ARMSEs、位置AAVBs、速度AAVBs 以及加速度AAVBs的统计结果.

表3 不同目标检测情形下状态估计的RMSEs对比

表4 不同目标检测情形下状态估计的AAVBs对比

首先，值得肯定的是，本文所提的算法在情形1～3下取得了较好的跟踪结果.这进一步验证了本文所提算法的有效性.此外，从目标检测概率πk在K/3 ≤k＜2K/3区间取不同概率值来看，随着目标检测概率πk的减少，状态估计的性能大致呈下降趋势.从相同的目标检测概率πk在不同区间的取值来看，越长的采样区间遭受强杂波（即较低的目标检测概率），状态估计的性能总体呈下降的趋势.需要特别注意的是，当目标检测概率分别取情形4和情形6时，都出现了滤波发散问题.这主要是因为随着目标检测概率的降低，收集的量测中包含有用的状态估计信息缺失严重.同样，当目标检测概率受强杂波影响的区间越多，强杂波淹没目标的有用信息越多，进而使雷达传感器难以接收到有效的量测对状态进行估计，出现了滤波发散问题.因此，为了有效的进行目标状态的估计，目标的检测概率不能太低，遭受强杂波的时间（即区间）也不能太长，否则无法收集到目标的有用信息对状态进行估计.总体上来说，本文提出的算法能鲁棒地、有效地应对强杂波的干扰，对目标的状态进行估计.

5 结论

本文设计了强杂波背景下含量测丢失的目标跟踪算法.该算法采用学生t分布来模拟非高斯噪声的长拖尾特性.通过引入伯努利随机变量，将求和形式的后验概率密度函数转换成乘积形式的概率质量函数，并构建了分层状态空间模型.在此基础上，设计了用于量测丢失的鲁棒变分贝叶斯平滑器.以机载雷达跟踪空中目标为例验证了本文算法的有效性.