一种高分辨雷达空间目标跟踪方法研究

刘政玮, 陈 映, 鲁耀兵

(北京无线电测量研究所, 北京 100854)

0 引 言

高分辨雷达通过发射宽带信号获取目标的高分辨一维距离像,从而实现目标的分类与识别[1-2]。由于其高距离分辨率,高分辨雷达在弹道导弹防御和空间目标监视等领域受到了越来越多的关注[3-4]。但是由于宽带检测、跟踪等关键问题尚未完全解决,现役的高分辨雷达只能采用窄带-宽带信号交替发射的工作方式,即主要通过窄带信号实现目标跟踪,而宽带信号则用于目标特征提取[5]。这种工作方式不仅增加了雷达系统设计的复杂性,还会造成雷达资源的浪费。因此,实现高分辨雷达目标跟踪对空间目标跟踪有着极其重要的现实意义。

空间目标运动在笛卡尔坐标系中更易于描述,在该坐标系描述的目标运动特征也更符合目标的本质运动规律,因此目标的运动方程通常建立在笛卡尔坐标系中,而雷达通常将目标量测描述为球坐标系下的距离、方位、俯仰[6]。由于目标状态方程和雷达方程通常建立在不同的坐标系下,目标的跟踪方法基本可以划分为3类[7]:混合坐标系的目标跟踪方法、笛卡尔坐标系下的目标跟踪方法以及球坐标系下的目标跟踪方法。混合坐标系下的目标跟踪方法是目前被广泛应用的滤波方法。大多数用于空间目标跟踪的非线性滤波技术,如扩展Kalman滤波器(extended Kalman filter, EKF)[8-10]、无迹Kalman滤波器(unscented Kalman filter, UKF)[11-12]、容积Kalman滤波器(cubature Kalman filter, CKF)[13]等均属于混合坐标系下的目标跟踪方法。这类方法的缺点在于严重依赖坐标转换及非线性函数的线性化方式。笛卡尔坐标系的目标跟踪方法将球坐标下的测量转换到笛卡尔坐标系,使得转换后的量测及观测方程表现为一种伪线性形式[14],在笛卡尔坐标系完成目标状态预测和状态更新。由于转化后的量测是有偏差的,会使得跟踪精度下降,因此必须使用各种补偿方法[15]。此时状态更新中非线性处理体现于伪量测对应的误差协方差矩阵的计算。球坐标系下的目标跟踪方法则是将目标的运动方程和观测方程在球坐标系下进行建模。显然,这类建模的优点在于测量模型是线性、非耦合的;量测噪声符合高斯分布假设;量测中的多普勒信息可被直接利用[16]等。缺点是球坐标系下的目标运动方程推导困难,推导后的运动方程是高度非线性的,各状态间也是高度耦合的。复杂运动的弹道目标、空间目标,球坐标系下的目标运动学方程的非线性程度远远高于笛卡尔坐标系下的非线性程度,以至于量测方程的非线性优势被抵消,导致滤波精度下降[17]。

对于同一空间目标,高分辨雷达带宽越大、波束宽度越宽,即雷达的距离分辨率越高、角度分辨率越低,量测的非线性就越高。随着量测非线性程度的增大,会出现一类特殊的非线性问题[18-19]:量测空间内呈现高斯状的噪声分布在更新过程中被转换到状态空间内时,其不确定性区域会变为强烈的非高斯分布,这种现象在高分辨雷达跟踪弹道目标、空间目标时经常出现。当高分辨雷达在跟踪空间目标时出现上述特殊的非线性问题时,EKF、UKF等滤波方法会出现跟踪精度下降,甚至滤波器发散的问题。因此,本文首先针对这一类特殊非线性问题进行了分析,并提出一种基于中间过渡状态的Kalman滤波方法。改进的Kalman滤波方法将量测空间内的目标位置量(距离、方位、俯仰)及其速度量(径向速度、方位角速率、俯仰角速率)作为中间过渡状态,则Kalman滤波更新过程被转换到量测空间中进行,实现了量测方程的线性化。仿真实验结果验证了改进滤波方法可有效提高高分辨雷达跟踪空间目标的精度。

1 空间目标跟踪方法

在空间目标跟踪中,通常使用状态量的状态方程描述目标运动特征,用量测方程描述雷达探测目标的位置信息,两者通常建立在不同的坐标系中。同时,由于空间目标的动力学特性复杂,因此雷达对空间目标的跟踪涉及坐标系统、运动建模和滤波算法等方面的内容。

1.1 雷达坐标系

雷达探测到的位置信息是对空间目标进行跟踪的基础,现有的单脉冲雷达多提供在雷达站球坐标系下描述的径向距离、方位和俯仰信息。在雷达站球坐标系下,空间目标的运动方程比较复杂,为简化跟踪过程,通常在直角坐标系下描述运动方程。根据实际应用情况,本文在雷达北天东直角坐标系中分析空间目标运动情况,在雷达球坐标系中探测目标径向距离、角度信息,两种坐标系的相互关系如图1所示。

图1 雷达站北天东坐标系与雷达站球坐标系Fig.1 Radar station north sky east coordinate system and radar station spherical coordinate system

雷达站北天东坐标系是一种站心坐标系,是以雷达中心为坐标原点Or,OrXr和OrZr轴是地球参考椭球的切线,分别指向北和东,而OrYr轴垂直于当地水平面向上。

雷达站球坐标系原点为雷达站,R为斜距,方位角A为正北轴OrXr顺时针旋转到目标与雷达站在水平面上的投影角度,俯仰角E定义为目标与雷达站连线与水平面的夹角。

1.2 空间目标运动模型及观测模型

目前,空间目标跟踪常用的3种地球模型分别是平地球模型、球地球模型和椭球地球模型[20],在实际应用中,可以根据空间目标的飞行距离相对于地球半径的大小选择合适的模型。对于长时间的空间运动目标跟踪需要考虑由地球形状的不规则引起的摄动,因此采用椭球地球模型可以获得更准确的重力加速度模型ag:

(1)

在雷达北天东直角坐标系中分析空间目标运动受力情况更加直接,过程噪声也更易于描述,因此大多数文献在直角坐标系中建立目标的运动方程。在雷达北天东坐标系中,设空间目标的状态变量为x=[x,vx,y,vy,z,vz],则空间目标运动状态微分方程可描述为[21]

(2)

(3)

式中:ae、ac为雷达北天东坐标系中由地球自转引起的牵连加速度矢量和柯氏加速度矢量,ad为气动阻力加速度矢量。3项加速度矢量的具体表达式可参考文献[20]。

通常在雷达测量坐标系下获得目标量测Z=[R,A,E],因此目标的量测方程为

(4)

式中:vR、vA、vE为相互独立的零均值高斯分布。

1.3 传统的空间目标滤波方法

目前空间目标跟踪方法中应用最为广泛的是混合坐标系下的目标跟踪方法,其因为直观的运动方程及量测方程建模方法,受到众多学者的关注[22-23]。其中,EKF通过泰勒级数展开的方法,获得非线性系统的线性近似,进而采用卡尔曼滤波处理非线性系统的滤波问题。UKF、CKF则是通过确定性采样点的变换求出后验概率密度函数的均值和协方差,避免了对非线性函数的近似。

另外一种将非线性过程线性化的方式是将目标的量测值转换到直角坐标系,在直角坐标系中进行目标状态的预测更新,这种方法称为转换量测卡尔曼滤波器(converted measurement Kalman filter, CMKF)。传统的CMKF存在偏差,可以通过将其乘以一个偏移因子实现去偏,但该方法存在兼容性问题[24],文献[25]提出的改进方法对其计算结果进行了修正。这些算法在计算转换误差统计特性时均以量测值为条件,会造成增益矩阵与量测误差相关,最终导致估计结果有偏。因此,有学者提出以目标预测信息为条件计算转换误差特性的无偏CMKF,从而克服滤波参数与量测误差之间的相关性[26-27]。但是,无论是以量测值为条件还是以预测值为条件的转换方法都会造成转换后的量测噪声是非高斯的,且在各坐标方向上互耦。

球坐标系下的目标跟踪方法在量测空间内选取目标运动状态,同样可以实现量测方程的线性化。这类方法需要直接在球坐标系下描述目标运动方程,但是对于空间运动目标,其在运动过程中受到各种力的作用,在球坐标系中难以直接进行建模。在实际应用中,当目标运动的假设条件发生变化时,需要重新对其进行建模分析。另一方面,在目标的运动方程转换到球坐标系时,通常不考虑过程噪声的转换,这时极易出现模型失配。

1.4 高分辨雷达跟踪空间目标的非线性问题

高分辨雷达在跟踪空间目标时,受到雷达带宽、波束宽度、信噪比等因素影响,高分辨雷达的测距精度高,测角精度较低。当笛卡尔坐标系中的状态使用来自球坐标系的非线性量测更新时,量测的不确定性区域会呈现出非常明显的非高斯分布。这种分布呈现出一种曲线形状(在二维空间内类似于香蕉形状,在三维空间内类似于隐形眼镜形状),因此将这一分布称为“隐形眼镜(contact lens, CL)”分布,并将这一类非线性问题简称为CL问题[28]。本文以极坐标系到二维直角坐标系转换为例,研究状态空间内的量测分布特性及高分辨雷达量测非线性的产生机理,雷达站球坐标系到雷达北天东坐标系的转换可由此类推。

令(x0,y0)为雷达直角坐标系下目标的真实位置坐标,(r0,θ0)为雷达极坐标系下目标的真实位置,则目标在雷达极坐标中的距离和方位量测值为

(5)

(6)

将其在目标位置真值(x0,y0)处进行泰勒展开得

(7)

由式(7)可以看出,在雷达直角坐标系下,目标的量测值不再符合高斯分布,且其真实分布复杂,难以获得。同时,式(7)表明直角坐标系下目标量测分布与目标的真实位置及测距和测角噪声有关。

图2从几何角度描述了这种量测的不确定性区域的形成原因。对于远程空间目标,在雷达测距精度不变的情况下,若雷达探测角度σθ误差增加,则跨角度误差范围rσθ相应增大,由于其远远大于测距误差σr,最终导致量测的不确定性区域逐渐呈现为曲线状。

图2 量测非线性CL问题Fig.2 CL problem of measurement nonlinearity

Lerro等[29]提出一种非线性度量指标来描述量测在坐标转换中偏离高斯分布的程度:

(8)

假设高分辨雷达对空间目标的作用范围为2 000 km,图3描述了不同β值对量测非线性的影响程度。图中蓝色点迹为104个蒙特卡罗仿真点迹,绿色实线为根据直角坐标系下量测值泰勒展开得到的3σ误差椭圆理论值,红色实线为根据蒙特卡罗仿真点迹得到的3σ误差椭圆真实值。由图3(a)可以看出,β值越大,高分辨雷达对空间目标量测的距离噪声误差越小。测角噪声增大时,由量测高斯假设得到的3σ误差椭圆并不能体现直角坐标系下量测分布的真实特性,会导致经典非线性滤波器的跟踪性能下降。从图3(b)可以看出,3σ误差椭圆可以很好地描述量测的真实分布,各非线性滤波器仍适用,但此时雷达带宽不符合高分辨雷达带宽。

图3 不同β值对量测非线性的影响Fig.3 Influence of different β value’s on measurement nonlinearity

当空间目标跟踪中出现强非线性CL问题时,传统的非线性滤波器的跟踪精度会急剧下降,甚至出现滤波发散。在图4中,黑色区域表示量测不确定性区域和预测状态不确定性区域的交集区。在滤波的初始阶段,目标状态估计不够准确,此时交集区域明显是弯曲且非高斯的,如图4(a)所示。基于高斯假设的EKF滤波器会出现一致性问题,EKF滤波器极有可能出现发散。当目标状态估计准确时,交集区可以近似地视作符合高斯分布,如图4(b)所示,此时EKF可以较好地跟踪目标,但由于其在线性化时只保留一阶偏导项,跟踪精度较差。在这一过程中,由于UKF、CKF通过采样点逼近量测不确定区域,因此其可以较好地实现目标跟踪。与EKF、UKF、CKF不同,CMKF是通过增加量测不确定性区域的“厚度”,即通过牺牲一定的测距精度来保证滤波器的一致性,使得交集区可以近似为高斯分布,如图4(c)所示,因此CMKF可以应对CL问题,但其测距误差会急剧增大[30]。

图4 滤波器性能几何图解Fig.4 A geometric illustration diagram of filter performance

2 改进的Kalman滤波器

2.1 矢量分析法求解不同坐标系下的目标状态转换

(9)

球坐标系中,目标的距离、方位、俯仰为

(10)

令ex、ey、ez分别为直角坐标系中OrXr、OrYr、OrZr轴的单位矢量,在球坐标系中取eR为沿径向方向的单位矢量,eA为垂直于径向方向并指向方位角A增大一方的单位矢量,eE定义与eA同理。因此,目标径向矢量可表示为

R=ReR=RcosEcosAex+RsinEey+RcosEsinAez

(11)

由此可得,球坐标系中各单位向量为

(12)

所以

(13)

所以,球坐标系中速度量为

(14)

将式(12)代入得

(15)

综上所述,直角坐标系中的目标速度为

(16)

球坐标系中,目标径向速度、方位角速率、俯仰角速率分别为

(17)

2.2 改进的Kalman滤波器

高分辨雷达在跟踪空间目标时易出现非线性问题——CL问题,主要原因在于球坐标系中的量测转换到直角坐标系时,其不确定性区域呈现出非常强烈的非高斯分布。因此,改进的Kalman滤波器将滤波器的更新过程转换到量测空间内进行,选取量测量(距离、方位、俯仰)及其一阶导数(径向速度、方位角速率、俯仰角速率)作为目标在量测空间内的状态,对目标进行更新,这样在量测分布符合高斯假设的同时,量测方程也相应地变为线性方程。

为方便后续方法的比较,假设空间目标的非线性系统状态估计的一般描述为

(18)

则改进的Kalman滤波器计算步骤如下:

步骤 1预测。假设预测滤波器输入为k时刻滤波结果xk=[xk,vxk,yk,vyk,zk,vzk]和Pk,通过下式获得状态预测结果xk+1|k和Pk+1|k。

xk+1|k=f(xk)

(19)

(20)

(21)

式中:Qk为过程噪声wk的协方差。

容积准则能够将笛卡尔坐标系下的高斯分布的运动状态转换为某个多维几何体容积的计算,具有较高的计算效率,能够获得较高的数值精度。因此,改进的Kalman滤波器根据三阶球面径向准则将状态空间内的预测协方差Pk+1|k转换为量测空间内的协方差PMk+1|k。

通过下式对xk+1|k采样,得到2n个采样点:

(22)

式中:ei表示第i个元素为1的单位向量。

采样点经过式(10)和式(17)转换后的值为Ci(i=1,2,…,2n)。

则量测空间内中间状态的协方差PMk+1|k为

(23)

(24)

步骤 3更新。由状态xMk+1|k及其协方差PMk+1|k得到量测空间内的状态估计结果xMk+1PMk+1。

由于更新过程在量测空间进行,因此原本的非线性量测矩阵h(xk)变为线性矩阵Hz,即

(25)

量测空间内目标运动状态更新过程为

(26)

(27)

xMk+1=xMk+1|k+KMk+1(zk+1-HzxMk+1|k)

(28)

PMk+1=(I6-KMk+1Hz)PMk+1|k

(29)

步骤 4状态估计结果提取。根据式(9)、式(16)及三阶球面径向准则将量测空间内的状态估计结果xMk+1和PMk+1转换得到k时刻直角坐标系中的目标滤波估计结果xk+1和Pk+1,这一步可认为是中间过渡状态转换的逆过程。

改进的Kalman滤波器与CKF均采用了容积变换准则,但是二者有本质的不同。由于CKF在混合坐标系中进行目标的预测更新,失去了Kalman滤波器本身具有的一致性和无偏性。改进的Kalman滤波器则是在预测后将状态空间内的目标预测值通过容积准则映射为量测空间内的中间过渡状态,从而得到线性的Kalman更新过程,更新后再将中间过渡状态逆变换到直角坐标系内,得到目标状态的一致性估计。同时在转换过程中,不同于CKF仅仅将空间目标的六维运动状态转换为量测空间内的三维位置预测信息,改进的Kalman滤波器选取了量测空间内目标的六维运动信息,进一步保证了转换过程的一致性。另一方面,相比CKF需要对复杂空间目标的运动方程进行采样,改进的Kalman滤波器能够在保证滤波一致性的同时,减少通过容积变换求取目标预测值及协方差的步骤,目标的状态估计是通过线性Kalman滤波器进行的,因此计算量会减少。

3 仿真实验

3.1 仿真场景

仿真的弹道目标射程为3 700 km,在雷达北天东坐标系下,仿真生成的弹道轨迹如图5所示,其中,“∘”和“△”分别表示弹道目标发射点位置和落点位置。

图5 雷达北天东坐标系弹道目标轨迹图Fig.5 Radar north sky east coordinate system trajectory diagram of ballistic target

雷达的探照范围为2 000 km,数据率为10 Hz,雷达跟踪弹道目标时长为100 s。在雷达球坐标系下,导弹的径向距离、方位角和俯仰角随时间的变化关系分别如图6、图7、图8所示,其中蓝色段为弹道目标轨迹,红色段为雷达探测弧段。

图6 目标距离随时间的变化关系Fig.6 Change of target distance with time

图7 目标方位角随时间的变化关系Fig.7 Change of target azimuth with time

图8 目标俯仰角随时间的变化关系Fig.8 Change of target pitch angle with time

3.2 高分辨雷达跟踪弹道目标滤波仿真

基于第3.1节仿真给定的弹道目标运动轨迹,本文采用典型高分辨雷达,测距误差为0.1 m。针对不同的测角误差,本文通过50次蒙特卡罗求均方根误差(root mean square error, RMSE)平均值的方法,进行UKF、CKF、EKF、两种无偏CMKF分别为无偏CMKF(unbiased CMKF, UCMKF[31])、基于预测信息的CMKF(prediction position CMKF, PRE_CMKF),以及本文改进Kalman滤波器的性能对比。滤波器均采用两点初始化方法,其他仿真参数设置如表1所示。

表1 跟踪场景

在仿真场景1中,雷达跟踪弹道目标时的测距误差较大,测角精度高,量测的非线性程度低。在这种情形下,通过图9可以看出,UKF、CKF、EKF及本文改进的Kalman滤波器的性能相近,而UCMKF的测距误差较大,PRE_CMKF相较于UCMKF滤波性能有所提升。

图9 仿真场景1滤波结果Fig.9 Simulation scenario 1: filtering results

随着雷达性能提升,高分辨雷达的测距精度提高,图10为仿真场景2时目标的状态估计结果。从图10中可以看出,高分辨雷达测角精度较高时,由量测转换引起的非线性程度较低,本文提出的改进Kalman滤波器和UKF、CKF测距、测角精度相近。而EKF滤波误差急剧增大,这是由于EKF在完成系统模型的线性化近似时忽略了高阶项,所以导致EKF估计的状态后验分布产生较大误差。而此时两种CMKF虽然可以保证目标的稳定跟踪,但它们的测距误差较大。

图10 仿真场景2滤波结果Fig.10 Simulation scenario 2: filtering results

由于当前高分辨雷达测角技术发展的限制,高分辨雷达跟踪弹道目标时测量角度误差较大,因此仿真场景3中假设方位角误差为σA=0.5°。在高分辨雷达中,随着测角误差的增大,由量测值转换引起的CL问题的非线性程度增大。从图11中可以看出,EKF会出现滤波发散的情形。此时,两类CMKF滤波器通过牺牲一定程度的测距精度,保证了滤波器的一致性,其仍然能够跟踪目标。本文提出的Kalman滤波器测距精度和UKF、CKF基本一致,但由于测角误差增大,UKF、CKF的方位角精度急剧下降,在弹道目标防御过程中将会导致一系列严重的后果,而本文提出的改进Kalman滤波器的角度RMSE仍可收敛至较小的测角误差。因此,当雷达球坐标系量测转换到直角坐标系中出现CL分布时,本文提出的改进Kalman滤波器比UKF、CKF具有更优的滤波性能,弹道目标跟踪精度更高。

图11 仿真场景3滤波结果Fig.11 Simulation scenario 3: filtering results

为了比较各滤波器的计算复杂度,图12列出了不同滤波器单次状态估计的平均时间损耗。从图中可以看出,基于中间过渡状态的Kalman滤波器在保证滤波精度的同时,可以避免对复杂的运动方程进行采样,计算复杂度比UKF、CKF更小。

图12 平均时间消耗Fig.12 Average time consumption

4 结 论

本文研究了高分辨雷达跟踪空间目标时,量测值由于坐标转换在状态空间内不确定性区域呈现出严重非高斯状态的问题。此时,传统的非线性滤波器的状态估计精度会急剧下降。本文提出一种基于中间过渡状态的Kalman滤波器,在直角坐标系下进行目标状态预测,避免了球坐标系下对目标运动的复杂建模;预测后选取量测空间内目标位置和速度量及其协方差作为中间过渡状态,在量测坐标系下进行目标状态线性更新;更新后再将中间过渡状态映射到直角坐标系内,从而得到目标状态的一致性估计。仿真结果表明,相对于现有的非线性滤波器,改进算法可有效提高高分辨雷达对空间目标的状态估计精度,这将对雷达系统快速准确形成目标轨迹有着重要作用。