基于空气动力学模型含雷达散射截面观测的联合目标跟踪识别

郭鹏，白亮，武梦洁，蒋宏

(1.北京航空航天大学飞行器控制一体化技术重点实验室，北京100191;2.洛阳电光设备研究所 火力控制技术国防科技重点实验室，河南 洛阳471009)

0 引言

联合目标跟踪识别(JTC)技术是信息融合领域新兴的一个研究方向［1］，大致可分为以文献［2 -4］为代表的基于质点运动模型的JTC 技术和以文献［5 -8］为代表的基于刚体运动模型的JTC 技术两类。基于质点运动模型的JTC 技术主要是依赖运动学信息来进行目标分类的，通常只能实现广泛意义上的分类，如“战斗机/客机”;而要识别目标飞机型号，如F-15 战斗机/VFY-218 战斗机，则需要采用基于刚体运动模型的JTC 技术，并在观测中至少包含一个目标特征，该特征有可能区分特定的目标型号。例如，文献［5］建议使用高分辨距离成像特征，文献［6 -7］选择了雷达散射截面(RCS)特征，因为这些分类特征不仅依赖运动学信息，还包含目标的属性信息，能够使分类结果更准确，同时，目标飞机型号的正确识别也有助于提高目标跟踪精度。

为了正确识别出目标的型号，向指挥和制导系统提供更精确的信息。本文以地基被动警戒雷达为背景，利用地面上多个商业广播所发射的无线电调频(FM)信号，提出一种有效的基于空气动力学模型含RCS 观测的联合目标跟踪识别技术。

1 建立目标飞机的状态模型

在使用真实的空气动力学模型来描述目标的平动和转动时，需要考虑特定飞机的具体参数，因此目标状态和目标类型之间是完全耦合的。

1.1 对飞机目标的平动建模

有3 种类型的力作用于飞机上:重力，推力和气动力。重力表示在惯性坐标系下，向下为正，W =mg，其中m 为飞机质量。推力的大小由飞行员控制，对于一个给定类型的飞机，推力有上界和下界，上界由引擎最大的输出决定，下界由维持足够升力的最低飞行速度决定。推力用T 表示，并假定它沿着机体坐标系的x 轴。气动力表示在气流坐标系下:

式中:D、C、L 分别为阻力、侧力、升力;CD、CC、CL分别为无量纲的阻力系数、侧力系数、升力系数;ρ 为飞机飞行高度上的空气密度;v 为飞机飞行速度;S机翼面积。只要侧滑角β 足够小，可假定CC≈0.对亚音速和超音速飞行，在一定的α(迎角)范围内，可以用以下近似:

式中:CLα为升力线斜率;CDmin为零升阻力系数;KD为展弦比，它们主要是飞机形状和速度的函数。在有限的速度变化范围内，可忽略这些参数对速度的依赖。考虑飞机所受的重力、推力和气动力，将飞机平动的总力F 表示在惯性坐标系下:

式中:Oib为机体坐标系到惯性坐标系变换矩阵;Obw为气流坐标系到机体坐标系变换矩阵。由于有CD、CL，(4)式不仅是一个刚体运动模型，还是一个设计明确的飞行。由(4)式可知，要计算F 需要两套参数。第一套是目标类型参数Ac={CLα，CDmin，KD，S，m}，下标c 表示这些参数所代表的目标类型。目标身份类型一旦被确定，这套参数就完全被确定了;第二套是目标运动学参数{pz，v，ψ，θ，φ，T}，其中ψ、θ、φ 为3 个欧拉角，z 方向的高度pz和空气密度ρ 相关。

1.2 对飞机目标的转动建模

使飞机产生滚动，俯仰和偏航的力矩都是在飞行员的直接控制下。因此，可以通过对飞行员的意图建模来代替对力矩建模。飞行员的意图对转动的影响可以定义为一组转动模态，每个转动模态对应一个不同类型的飞行，如转弯或爬升。这样，飞机的转动只不过是飞行员命令的一系列未知机动{γ1，γ2，…，γk}的结果，其中γk为飞机从kΔt 到(k +1)Δt时刻的操作模态。

对于状态模型，假设目标的转动可以用两个纵向模态和一个横向模态描述。这3 个模态为:模态1(CV)，常速固定高度飞行;模态2(PT)，纵向爬升或下降;模态3(CT)，恒速率协调转弯。

状态模型需要随时间变化的欧拉角，因此可以绕过力矩建模问题，直接对每个转动模态进行更有意义的欧拉角建模。

表1总结了每个转动模态下的欧拉角模型。其中，带“* ”变量表示驱动飞机转动的未知飞行员输入，ω*表示转弯角速率，Δθ*表示俯仰角的变化，αss表示以规定速度维持水平飞行所需的稳态迎角。最后，假设模态序列{γk}按照模型转移概率矩阵为P的一阶马尔可夫链演变，矩阵P 中的元素被定义为P(i，j)=Pr (γk=j|γk-1=i).

表1 欧拉角模型Tab.1 Euler angle model

1.3 完整的飞行模型

图1给出一个由xk-1到xk的伪码描述。每次扫描，首先判断该模型是否发生了模态改变:如果没有改变，上一时刻的驾驶员输入被添加少量噪声而作为当前时刻的驾驶员输入;如果有改变，需要获取新的飞行员输入。图中Tss为以规定速度维持水平飞行所需的稳态推力，分别为纵向爬升的俯仰角变化、纵向爬升的推力以及恒速率协调转弯的转弯角速率，都服从均匀分布，且新的飞行员输入从这些均匀分布中取值。通过让这些分布依赖于目标类型，每个飞机实际的运动学限制可以被使用，从而有利于目标识别。

图1中伪码的最后一行调用一个数值积分子程序，它利用k-1 时刻的状态和当前飞行员的输入，完成了飞机的平动和转动，并将状态从k-1 时刻推进到k 时刻。图2给出了这个数值积分子程序的伪码。

图1 飞机状态模型伪码Fig.1 Pseudo-code of state model for flight

到此为止，建立了一个基于飞机实际的平移和转动的状态模型。

2 建立被动雷达的观测模型

被动雷达可以观测到时延、多普勒频移及RCS.在三维空间要唯一地定位目标，至少需要3 个源于不同发射器的延迟测量;同样，要唯一地确定目标的三维速度，也至少需要3 个源于不同发射器的多普勒频移测量。因此，为简单起见，采用3 个发射器，单被动雷达接收器。这样，k 时刻被动雷达的观测为Zk={(τm，k，dm，k，σm，k)}，m =1，2，3，其中(τm，k，dm，k，σm，k)分别为k 时刻源于第m 个广播电台的时延、多普勒及RCS 观测量。

图2 状态数值积分子程序伪码Fig.2 Pseudo-code of state integration subroutine

2.1 对时延和多普勒频移建模

式中:τm，k为时延;pk为k 时刻目标的位置;pr为被动雷达接收站的位置;c 为光速。

式中:dm，k为多普勒频移;vk为k 时刻目标的速度矢量;nm为位置的单位向量;λm为电磁波的波长。

式中:ptm为第m 个广播电台的位置;pr为被动雷达接收站的位置，如果让惯性坐标系的原点与被动雷达接收器重合，则pr=0.时延和多普勒频移的测量噪声可假设为零均值高斯白噪声，可得时延、多普勒观测量的似然函数为

式中:h(τ)m，k、h(d)m，k由(6)式和(7)式定义。

2.2 对RCS 建模

如果假设被动雷达系统是相干的(通过使用相当有限的调频无线电信号带宽)，理论上可以还原k时刻源自第m 个发射台的复数散射系数sm，k.然而，用于前端滤波器组的积分时间不能用散射系数的相位来分类，故使用|sm，k|，m =1，…，M 来分类，因为RCS 被定义为σm，k= |sm，k|2，这相当于使用RCS.RCS 是一个非负量，假设其测量噪声为高斯分布不恰当，因此可以对(相干)匹配滤波器的同步和正交分量的热噪声建模。对RCS 测量噪声，采用的具体模型为

式中:上标代表sm，k的实部和虚部;w1，w2为独立同分布、零均值、方差为的高斯随机变量。假定热噪声的方差独立于收到的FM 信号，可去掉的下标m.最终生成的似然函数可以表示为

式中:I0是第一种零阶修正贝塞尔函数;σm，k是k 时刻对源自第m 个发射台的RCS 观测量;σ2s 是影响散射系数sm，k测量的热噪声方差，m=1，…，M.

2.3 完整的观测模型

假设各观测量τm，k、dm，k、σm，k的噪声是独立的，粒子滤波所需要的观测似然函数将是各观测量似然函数的乘积，即

式中:p(σm，k|xk)由(12)式决定。

3 Matlab 实时调用FEKO 软件

FEKO 是一种用于RCS 理论计算的电磁仿真软件。可以由目标的实时位置、姿态角、电磁波发射站和雷达接收站的位置得到FEKO 软件的输入信息，即入射电磁波的方位角、高低角，散射电磁波的方位角、高低角。然后经Matlab 调用FEKO，可以得到不同型号飞机的实时RCS.调用过程如下:

1)运行PREFEKO，对飞机模型进行预处理，飞机模型如图3所示。

2)运行EDITFEKO，对“.pre”文件中的控制量、激励和远场等信息进行配置。配置好后的参数如图4所示。

3)运行FEKO 软件，对配置好的模型求解RCS值。

图4 EDITFEKO 界面Fig.4 The interface of EDITFEKO

4)运行Matlab 的m 函数实现调用，获取RCS.

4 粒子滤波执行

由于RCS 观测量，目标姿态和目标位置等状态量之间没有解析关系，并且飞机实际平动和转动的状态模型没有简单的显式表达式，因此不能使用扩展卡尔曼滤波，然而粒子滤波提供了解决方案。

正则化粒子滤波［9］(RPF)是一种基于重采样的改进算法，它与传统粒子滤波的区别在于:传统粒子滤波从离散近似分布中重采样，而正则化粒子滤波从连续近似分布中重采样。正则化粒子滤波可以有效缓解重采样造成的粒子多样性丧失问题，能够获得较高的精度。其算法如图5所示，S 表示统计协方差阵，K 表示Epanechnikov Kernel，一般情况下取近似值，用高斯核代替，hopt表示最优带宽。

目标有3 种飞行模态，本文将采用交互式多模型正则化粒子滤波算法来实现，并用交互式多模型粒子滤波算法与其比较。

5 仿真实验

5.1 仿真场景及参数设置

在仿真过程中，对战斗机VFY-218 进行跟踪和识别，无人机UAV 充当目标库中的参考目标。采用Monte Carlo 仿真来验证算法的有效性，仿真的次数为NMC=50，每次仿真的采样次数为Nk=75.仿真软件为Matlab 和FEKO，计算机的硬件为AMD Athlon(tm)II X2 250 Processor 3.01 GHz，1.75 GHz 内存。

图5 正则化粒子滤波Fig.5 Regularized particle filter

目标初始位置p0=［-95 000，50 000，4 000］T，初始速度为v0=［110，160，0］T，初始俯仰角、滚转角和偏航角分别为θ=4°，φ=0°和ψ=55.49°.3 个发射器和被动雷达的位置分别为pt1=［48 000，-18 000，0］T，pt2=［- 23 000，38 000，0］T，pt3=［42 000，23 000，0］T，pr=［0，0，0］T.采样周期Δt=0.4 s，积分步长Δi =0.2 s，3 个发射器交替向被动雷达发送信息以获得观测值，每轮处理3 个发射器信息的时间为1.2 s.粒子滤波的粒子数Np=600，模型转移概率矩阵P =［0.98，0.01，0.01;0.01，0.98，0.01;0.01;0.01，0.98］，初 始 模 型 概 率μcv(0)=0.8，μct(0)= 0.1，μpt(0)=0.1;传输时延的观测噪声σ2τ=0.000 1 ms2，多普勒频移的观测噪声σ2d=1 Hz2，RCS 的观测噪声σ2s=64 m4.

目标的运动模态如表2所示，其真实轨迹和跟踪轨迹如图6所示。

5.2 跟踪结果及分析

图7、图8和图9分别给出了姿态角跟踪图、速度跟踪图和位置均方根误差(RMSE)，图中的竖直虚线表示该时刻目标运动模态发生变化。

表2 目标的运动模态及参数Tab.2 The motion model and parameters of target

图6 目标运动轨迹Fig.6 Target motion trajectory

1 ～20 采样区间在固定高度做匀速直线运动，飞机受力平衡，速度和姿态角基本不变，此时交互式多模型粒子滤波(IMMPF)和交互式多模型正则化粒子滤波(IMMRPF)在x、y 和z 方向的RMSE 相当，都较小。

21 ～35 采样区间在竖直平面内做爬升运动，飞行员给飞机施加的推力T 大于空气阻力与重力的分力，俯仰角变化Δθ*=6° ＞0°，故飞机向上爬升。由于飞机在纵向发生机动，故此时的俯仰角变化和速度变化存在跟踪滞后现象。与匀速运动相比，爬升运动的IMMPF 和IMMRPF 在x、y 和z 方向的均方根误差都增大，且增大的幅度相当。

36 ～50 采样区间在固定高度做匀速直线运动，由于受到爬升运动俯仰角和速度滞后的影响，z 方向速度不能减小到0.在x 方向上，IMMRPF 均方根误差要小于IMMPF，而在y 和z 方向上均方根误差相当。与爬升运动相比，匀速运动的IMMPF 和IMMRPF 在x、y 和z 方向的均方根都减小，表明目标在非机动情况下的跟踪性能要优于机动情况下。

51 ～75 采样区间做恒速率协调转弯运动，在水平方向飞行员施加的推力和空气阻力平衡，飞机机身倾斜使升力在竖直方向的分力平衡重力，在水平方向的分力提供向心力。此时x 方向的速度不断增大，y 方向的速度不断减小，滚转角从0°变到-75°.IMMRPF 的偏航角在协调转弯过程中都要优于IMMPF，滚转角在51 ～60 采样区间也优于IMMPF.在x 和y 方向，IMMRPF 的速度相对于IMMPF 更接近真实值，因此IMMRPF 的均方根误差明显小于IMMPF.

图7 姿态角跟踪图Fig.7 Euler angle tracking

由图7和图8可以看出，目标发生机动之后，对于目标状态的跟踪存在滞后现象。这种滞后反映在模态概率上为机动检测时刻总是滞后于机动发生时刻。如图10所示，IMMPF 的模型概率在51 ～60 采样区间要略优于IMMPF，其他时候相当。

在第n 次Monte Carlo 仿真中，第k 次采样时刻的位置误差定义为

图8 速度跟踪图Fig.8 Velocity tracking

式中:NL为跟踪失败的次数。

表3 跟踪性能Tab.3 Tracking performance

从表3可以看出，IMMRPF 的综合跟踪性能要优于IMMPF.

图9 位置均方根误差Fig.9 The position RMSEs

图10 3 个模态的概率Fig.10 The probabilities of three models

5.3 分类结果及分析

在第k 次扫描时，目标类型的概率可以表示为式中:Nc表示类型的数目;c 表示从离散类型标签中获取信息的判断准则;z1:k表示1 到k 时刻的观测序列。而Pr(c(xk)=j|z1:k)的表达式可以利用粒子滤波的权值近似表示为

从表4可以看出，在联合目标跟踪识别中，IMMPF 和IMMRPF 正确识别目标的概率都较高，因而能够较准确识别出目标的类型，并且IMMRPF 的类型平均识别率优于IMMPF.IMMRPF 利用了正则化粒子滤波的优势，从后验密度的连续近似采样得到新的粒子，因此避免了大权值粒子被多次复制，保持了粒子的多样性，故类型平均概率要优于IMMPF.

表4 类型平均识别率Tab.4 The average identification probability

在正确识别出战斗机的类型后，分类器将战斗机的RCS 值和类型参数Ac={CLα，CDmin，KD，S，m}输入到跟踪器中，其中CLα= 2.75 rad-1，CDmin=0.013，KD=2.68，S=42.7 m2，m=13 397 kg.

6 结论

本文对飞机目标建立了平动和转动的空气动力学模型，利用各飞行模态的特点将全状态变量简化为关注目标类型、飞行模态和飞行员输入;通过Matlab 调用FEKO 电磁仿真软件，将跟踪器的输出与各传感器的位置信息综合换算得到FEKO 电磁仿真软件的输入，有效缩短了求解RCS 的时间，然后将分类器得到的实时RCS 值和类型参数输入到跟踪器当中，完成目标的联合跟踪识别。

仿真实验结果表明，粒子滤波能够较好地处理高度非线性的系统模型，IMMRPF 在位置均方根误差、跟踪成功的次数和类型平均概率等方面都要优于IMMPF，具有较好的综合跟踪性能。

References)

［1］ 单甘霖，梅卫，王春平.联合目标跟踪与分类技术的进展及存在问题［J］.兵工学报，2007，28(6):733 -738.SHAN Gan-lin，MEI Wei，WANG Chun-ping.The development and problem of joint target tracking and classification technology［J］.Acta Armamentarii，2007，28(6):733 -738.(in Chinese)

［2］ Yang W，Fu Y W，Long J Q，et al.Joint detection，tracking and classification of multiple targets in clutter using the PHD filter［J］.IEEE Transactions on Aerospace and Electronics Systems，2012，48(4):3594 -3609.

［3］ Melzi M，Ouldali A.Joint multiple target tracking and classification using the unscented Kalman particle PHD filter［C］∥2011 IEEE 9th International New Circuits and System Conference.Bordeaux，FR:IEEE，2011:534 -537.

［4］ MEI W，SHAN G L，LI X R.Simultaneous tracking and classification:a modularized scheme［J］.IEEE Transactions on Aerospace and Electronic Systems，2007，43(2):487 -506.

［5］ Miller M I，Srivastava A，Grenander U.Conditional-mean estimation via jump-diffusion process in multiple target tracking/recognition［J］.IEEE Transactions on Signal Processing，1995，43(11):2678 -2690.

［6］ Guo Y F，Peng D L，Chen H J，et al.A kernel particle filter algorithm for joint tracking and classification［C］∥2012 15th International Conference on Information Fusion.Singapore:IEEE，2012:2386 -2391.

［7］ Herman S，Moulin P.A particle filtering approach to FM-band passive radar tracking and automatic target recognition［C］∥2002 IEEE Aerospace Conference Proceedings.Big Sky，Montana:IEEE，2002:1789 -1808.

［8］ Herman S M.A particle filter approach to joint passive radar tracking and target classification［D］.US:Graduate College of the University of Illinois，2002.

［9］ 罗飞腾.目标跟踪的粒子滤波技术研究［D］.合肥:中国科学技术大学，2010:16 -17.LUO Fei-teng.Research on particle filtering in target tracking［D］.Hefei:University of Science and Technology of China，2010:16 -17.(in Chinese)