基于相位分析与散射点关联的速度估计方法

陈学斌, 叶春茂, 张 彦, 胡庆荣

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

0 引 言

高精度的目标径向速度估计对于目标跟踪和识别而言具有重要意义[1]。针对雷达目标的速度估计,当前常用的方法主要是将目标物视作一个点,并通过对该点目标的相参相位回波信息进行分析,从而获得对该目标物的高精度的速度估计值[2-4]。

然而,当面对宽带条件下的扩展目标时,若直接将其视为多个散射点目标的集合,再用窄带多普勒测速方法获得各散射点的速度估计值,并取估计值的平均数作为扩展目标整体平动速度的估计,则该获得的平动估计值准确度并不理想。一方面,由于扩展目标会受到姿态变化、部件微动、噪声等因素的影响[5-7],各散射点的运动并不只含有平动的成分,简单的均值处理会降低平动速度估计的精度和稳健性;另一方面,同一距离单元内可能存在着与目标整体速度差异较大的部件散射点成分(如直升机旋翼、飞机螺旋桨等),也会使该单元慢时间回波的解模糊处理复杂。因此，直接使用窄带多普勒测速方法来估计扩展目标整体平动速度时其估计的准确性会受到限制,必须辅之以更加精细化的信号处理方式。

针对宽带条件下的扩展目标,另一类方法则是根据扩展目标整体的运动趋势来进行运动参数估计。其中,根据目标回波的包络走动情况进行平动速度估计的方法最为常用。该类方法通过对相邻复(实)包络作相关处理并取相关结果峰值的方式[8],或通过对一段时间内目标回波包络所移动的距离使用设计的多项式距离模型进行拟合的方式[9],以达到运动参数估计的目的。然而,该方法会受到纵向距离分辨能力的限制,测速精度有限。

为了更准确地获取扩展目标的平动速度信息,本文立足于脉冲类的相参雷达信号,提出了一种基于相位分析和散射点关联的估计方法。该方法仍采用分析回波相位中的多普勒信息进行速度参数提取的方式,以避开纵向距离分辨能力对估计精度的影响。在此基础上,使用预测-观测相关联的方式,将运动状态明显不属于平动方式的散射点运动分量进行剔除,从而通过精细处理的方式提高扩展目标物平动速度的估计精度。与此同时,本文所设计的方法是一种递进式的平动速度估计算法,在系统计算能力不断增强的今天,其在对宽带雷达信号的即时精确处理中也具备积极的意义。

本文的内容安排如下:第1节介绍本文方法所基于的回波信号的相位模型与用于相位分析的自适应Chirplet分解方法;第2节介绍用于区分平动/非平动成分的预测-观测匹配关联方法;第3节则介绍本文所提估计方法的处理过程,并用实验结果检验了该方法的有效性。

1 回波信号的相位模型与分析

1.1 扩展目标的运动分量

在窄带条件下,目标物可视为一个点,其运动可以用该点在雷达射线维度上相对雷达所在地的距离变化进行描述;而在宽带条件下,由于距离分辨能力的提升,扩展目标相对于雷达的运动应分为平动分量、转动分量以及微动分量来分别进行更加细致的描述[10-12]。

如图1所示,假设某飞机在一段时间内从位置①处飞到了位置②处。

图1 目标相对于雷达的运动Fig.1 Motion of target relative to radar

若在窄带条件下,飞机被视作一个点目标,其运动描述只需要获得该段时间内点目标所处位置到雷达处的雷达射线维度上的距离变化即可;然而,若在宽带条件下,飞机应视作一个扩展目标,飞机上的各处散射中心(如:尾翼、机翼上的散射中心)除了都具备飞机整体相对于雷达射线的运动外,还由于飞机相对于雷达射线的姿态变化而具备在雷达射线方向上各不相同的走动。由于姿态变化可视为飞机相对于雷达做了旋转运动,因此将这种不同的走动称为目标运动的转动分量;同时,称飞机整体相对于雷达射线的运动为平动分量[11]。此外,机翼的振动、螺旋桨的旋动等部件散射中心除平动、转动分量外还具备自身运动引起的微动分量,这就使对扩展目标上各散射中心的运动描述需要更加精细化。

因此,要对扩展目标的运动进行确切描述,就需要分别对各散射中心的平动、转动、微动分量进行相应的估计。而对于转动、微动分量的估计又建立在平动分量得到精确补偿的基础之上。故研究如何精确地获得扩展目标的平动速度估计并将之用于目标平动补偿,有利于精细描述扩展目标的运动情况,进而服务于宽带条件下的目标成像与识别需求。

1.2 回波信号的相位模型

对于一个较大型的雷达目标物体,当运动持续时间段足够短时,可以认为该物体在该段时间内所受到的外力合力为一个定值,即可认为该段时间内物体平动的加速度保持恒值,物体整体沿着某一直线做匀加速运动。现假设该物体上存在某一个散射点,若将该散射点的运动轨迹投影到雷达射线的径向方向上,则其所处的径向距离随时间的变化可表示为

(1)

式中:R0为初始时刻散射点所处的径向距离;v0表示初始时刻散射点的径向速度;a则为散射点的径向加速度。

现假设雷达向目标发射脉冲串信号,每一脉冲的发射信号形式表示为

st(t)=p(t)exp(j2πfct)

(2)

式中：fc表示载频。则处于目标物上的该散射点Pi的雷达回波经下变频后其可表达为

(3)

将式(1)代入式(3)中,只考虑相位项而忽略包络项的影响,则该散射点的回波可表示为

(4)

再令

(5)

则,忽略常数项,式(4)可化为

(6)

若视ω0为初始频率,γ为调频率,则式(6)表现为一个线性调频信号。

由上述可知,只要能够获取到某散射点在足够小段时间内的回波相位历程,并从中估计出其所含的ω0和γ的值,便能通过式(5)反算出该段时间内该散射点的初速度与加速度的值。而对ω0和γ的估计又建立在对相位信息的分析上,这就使得时频分布类的方法[13-16]对速度、加速度估计问题得以适用。

鉴于宽带条件下的扩展目标可视作多个散射点目标的集合,而单个散射点目标的回波相位模型又具有线性调频特性,本文采用Chirplet分解方法[17-20],来完成对目标上各散射点回波相位的时频分析。

1.3 自适应Chirplet分解方法的使用

对于宽带条件下的扩展目标,其脉冲串回波经脉冲压缩、包络对齐[21-23]操作后,其上的散射点因初始时刻所处的雷达射线径向位置不同而分布到不同的纵向距离单元中。而位于同一距离单元内的散射点又因受到姿态变化、部件微动、噪声等因素的影响,具备不同的运动参数。这就使得采用自适应Chirplet分解方法逐纵向距离单元地提取运动参数不同的散射点、获取其运动参数估计值,并最后综合地估计整体目标的平动参数成为可能。

自适应Chirplet分解是一种将待处理慢时间信号s(tm)近似表示成多个Chirplet小波基的线性叠加的方法[24],其可表示为

(7)

式中：≈表示信号s(tm)中可能还含有噪声等成分;Cn表示信号s(tm)在Chirplet小波基gn(tm)上的投影值;而Chirplet小波基则可表示为

(8)

要逐个纵向距离单元对信号进行Chirplet分解,就需要对每一单元内的慢时间信号所含的每一Chirplet分量进行参数估计,文献[24]提供了一种解决方法,该方法可在不发生速度模糊的条件下,通过Radon-Wigner变换[11]、解线性调频、窄带滤波、一维搜索等方式,从回波相位中获得高精度的ω0和γ的估计值。

在单一Chirplet分量参数估计的基础上,依照CLEAN准则[25]逐次地对回波中含有的每个Chirplet分量进行参数估计并提取,从而完成对单个纵向距离单元内信号的分解。此外,还可引入RELAX技术[26],以迭代方式使已估计出的运动参数组能以更新的方式进一步提升估计精度直至算法收敛。

需要注意的是,当回波信号的脉冲重复频率(pulse repeat frequency,PRF)较低时,频谱折叠、速度模糊的情况将会出现[27]。在这种情况下,便需要使用“折叠因子”这个先验信息,对在上述过程中获得的初速度估计值进行解折叠处理,具体计算可表示为

(9)

(1) 对回波先进行包络对齐操作,并记录下包络对齐时各脉冲包络所平移的距离;

(10)

其中,

(11)

则当Δ取整数时,折叠因子k的估计值为

(12)

而当Δ不取整数时,k的估计值为

(13)

其中,ceil(Δ)表示取大于Δ的最小整数。

应注意到,式与式的取值是建立于速度不模糊范围为(-|vmax|,|vmax|]的基础上,即相对的多普勒不模糊范围为[-fr/2,fr/2)。

2 预测-观测的散射点关联

在使用上述方法获得扩展目标上各纵向距离单元上的各散射点的解折叠初速度值后,便可经式(5)转换成初始频率,从而进一步得到扩展目标上各散射点在初始时刻t0于距离-多普勒域上的分布。

(14)

(15)

而用式(14)和式(15)进行逐散射点的预测,便能获得目标在时刻t1的距离-多普勒分布预测值。

在此基础上,若在时刻t1时,雷达获取到了目标实际的、解频域折叠后的散射点距离-多普勒分布观测值,便能将预测值与观测值关联起来,实现散射点的跟踪,并根据跟踪轨迹区分出非平动运动分量,以提高平动速度估计的精度。

2.1 同量纲特征向量的构建

由于噪声、散射点闪烁、参数估计误差等影响,各散射点在距离-多普勒维度上的预测值和观测值并不完全重合。此时,需要建立一个特征向量,来量化地表征各个值之间的匹配程度。其中,最为直观的特征向量构建方式便为

(16)

在预测和观测进行关联时,每一个预测值或观测值都对应一个特征向量。通过这样的映射,预测和观测两个值之间的差异便转化为二者所对应的特征向量之间的差异。而向量间的差异则引入向量间距r的概念进行表示:

(17)

式中：X和Y分别为预测、观测值的特征向量。

同时,应注意到,式(17)所示的向量距离的计算中存在一个加号,而加号两侧项的量纲是不同的,这不符合加法的意义。为此,应将向量内的元素置为同一量纲,将式(16)化为

(18)

式中:R表示纵向距离位置;f表示多普勒位置;T表示预测的间隔时间长度。此时,向量内的两个元素便具有了同种量纲(m/s),而向量间距也转化为

(19)

式中：ΔR和Δf分别表示预测和观测两种值之间纵向距离位置和多普勒位置的差值。

2.2 正确关联的概率表示

由式(19)可以看出,每一特征向量间距都对应了一条预测点与观测点间的连接。而某一连接为正确匹配的概率则与该连接所对应的向量间距大小呈反相关关系,即表示为

pip{第i个连接为正确匹配}

(20)

式中：ri表示第i个连接所对应的向量间距。

为构造正确关联的概率表达式,现先考虑从单一预测点向所有观测点映射的情况。图2(a)展示了这种映射,并标明了各个连接所对应的向量间距ri, j和概率pi, j,其中i表示第i个预测点,j表示第j个观测点。现假设对于预测点分布中的任意一个点,在观测点分布中一定存在一个点与之匹配,则图2(a)所标出的所有连接就构成了对第2预测点的正确关联概率的一个完备空间,即满足:

(21)

而由式(20)所示的关系,可做出以下定义:

(22)

(23)

图2 预测-观测的两种映射情况Fig.2 Two kinds of mappings for prediction-observation

需要注意的是,图2所示的两种映射情况下的概率定义是建立在不同的概率空间上的。因此,需要用一个方法将这两种概率综合起来,使每一个预测-观测连接都只对应一个特定的、可区分的概率值。

针对不同的概率空间,本文使用权数相加的方式,构造一个综合的概率模型:

(24)

式中：C表示权重因子。

ri, j≡p

(25)

此时,各连接综合概率的差异与特征向量间距无关,只取决于预测点与观测点的点数差异。若令M表示观测分布中的所有点数,N表示预测分布中的所有点数。则由式(22)与式(23)有:

(26)

(27)

此时,式(24)则化为

(28)

C/M=(1-C)/N

(29)

经变化,可得到因子C的表达式:

C=M/(M+N)

(30)

将式(30)代入式(24),可得到综合概率表达式:

(31)

由式(31),便可定义并计算得出所有预测-观测连接所对应的关联正确概率值。

2.3 基于概率的关联选择

由式(31)给出的定义,可以计算出所有预测-观测连接的正确概率值。根据这些计算出的概率值,便能从大到小地依次选出所需要的关联匹配连接,直到剩余的连接对应的概率值全都小于某一阈值为止。关联选择的步骤如图3所示,图中将已选出点移出搜索范围的操作则能避免多对一的映射情况的发生。

图3 关联选择流程Fig.3 Flow chart of choosing association

通过完成预测-观测之间的连接,相邻两次观测所得的各散射点也得以关联,从而实现目标上散射点的跟踪,便于平动参数的估计。

3 扩展目标的平动速度值估计

3.1 扩展目标的平动速度值估计步骤

针对扩展目标,本文所提出的基于相位分析与散射点关联的速度估计算法步骤如下。

步骤1获得以t0为初始时刻的一段经脉冲压缩后、时间跨度为T的脉冲串回波信号。

步骤2对该段信号做包络对齐操作,并对平移所得的数据做二次三项式拟合处理,根据式(11)～式(13)计算得到折叠因子的估计值,作为先验信息用于初速度解模糊处理。

步骤3在包络对齐后,逐纵向距离单元地进行慢时间域的自适应Chirplet分解,获取扩展目标上各散射点在初始时刻的初速度、加速度的估计值。并根据式(5)获得初始时刻散射点的距离-多普勒二维分布。

步骤4根据估计的初速度、加速度值,以匀加速直线运动模型,预测出经过一段时间T后,t0+T时刻的各散射点的距离-多普勒二维分布。

步骤5获得以t0+T为初始时刻的一段时间跨度为T的脉压后脉冲串回波信号,重复步骤2与步骤3,获得观测所得的t0+T时刻散射点的距离-多普勒二维分布。

步骤6将预测、观测所得的t0+T时刻散射点的距离-多普勒二维分布图,用上文所述的基于概率的预测-观测关联方法进行两图中的散射点关联。

步骤7以预测的分布图为中介,将观测所得的t0时刻、t0+T时刻的散射点分布图进行关联,确定各被关联的散射点在时段T内的运动轨迹,并去除斜率差异大的轨迹。

(32)

计算出该段时间内各个时刻的平动速度估计。

值得一提的是,步骤7中去除了斜率差异大的轨迹,其原因在于:目标在做平动运动时,在平面波照射的近似条件下,目标上各散射点相对于雷达的距离变化量相同,距离像是不变的[11]。而运动轨迹中斜率差异大的轨迹会使得目标距离像发生较大变化,因此不能近似视作平动分量,应予以剔除。

上述步骤描述了如何获得某段时间内目标平动速度的估计值,图4展示了如何使用本文所提方法对扩展目标进行多时段平动速度估计的算法整体流程图。从图4所示流程图可以看出,本文所提方法是一种递进式的平动速度估计算法,通过不断输入时间跨度为T的脉压后脉冲串数据,算法可不断输出t0+mT到t0+(m+1)T时段中的各时刻目标所具备的平动速度的估计值,其中m=0,1,…,M-1,而M则表示需要输出速度估计值的、跨度为T的时间段数目。因此,当系统的计算能力足够强时,该递进式的平动速度估计算法能即时地对不断接收、采样、获得的脉冲串信号进行处理,并估计得到各时刻目标物整体的平动速度值,以服务于扩展目标物的即时成像与识别的需求。

图4 平动速度估计方法流程图Fig.4 Flow chart of translational velocity estimation method

3.2 扩展目标的平动速度值估计实例

现有一段展示民航飞机转向飞行时,其径向距离随时间变化的L波段雷达实测数据。该回波的波形成分为线性调频脉冲信号,系统使用先直接采样、再匹配滤波的方式进行脉冲压缩,以保证各脉冲间相参的相位信息。图5展示了该数据中的一段包络分布情况。

图5 飞机包络走动情况Fig.5 Change of airplane’s envelope

考虑到飞机所受合力在足够短时间内可视为恒定,即加速度为恒定值，若假设1 s的时间内飞机的加速度保持不变,且已知该段回波的脉冲重复周期约为0.004 s,则可将回波以256个脉冲作为一个时间跨度T进行分段,而在每一分段时间内都认为飞机在雷达射线的径向维度上做匀加速直线运动。

现对图5所示的扩展目标分段回波用滑窗的相邻相关法[21]做包络对齐的预处理,并逐段通过拟合获取到各段初始速度的折叠因子估计值。然后,将对齐后的每一个纵向距离单元进行慢时间域的自适应Chirplet分解与分析,获取各散射点在各段时间内具备的初速度、加速度的估计值。在获得了运动参数估计值后,便能依据匀加速直线运动模型预测出一段时间后各散射点所处的位置,并与实际观测值进行关联匹配。

针对图5所示回波,本文将对该回波的每一分段的相位分析都视作一次观测,并获取到散射点在每一分段初始时刻的距离-多普勒二维分布图,再逐分段地进行相邻两次观测之间的分布图关联。

以第4分段为例,图6展示了该分段时间内散射点的预测变化轨迹。图6中,三角形表示本次观测所得的初始时刻散射点的距离-初始多普勒频率分布,菱形则表示对这些点在下一次观测时散射点所处的初始位置的预测。

图6 第4次观测值向其预测值的映射Fig.6 Mapping from No.4 observation values to their prediction values

有了预测值,便可以与第5次的观测值进行匹配关联。本文使用上文的基于概率的预测-观测关联方法,先对每一连接建立同量纲特征向量,再使用式(31)所示的概率定义进行关联匹配的选择。最终得到的关联结果如图7所示,菱形表示第5次观测值,三角形表示第4次预测值。由图7可以看出,预测值与观测值都分布在同一区域内,两种点值的分布中有部分得到关联,其余则不被关联。

图7 第4次预测值向第5次观测值的映射Fig.7 Mapping from No.4 prediction values to No.5 observation values

根据图6和图7的映射,便可将预测值点作为中介,将第4和第5次观测值连接起来。连接出的散射点运动轨迹如图8所示,三角形表示第4次观测值,菱形表示第5次观测值。对比图6与图8,可以发现一些运动轨迹明显与整体目标不符的散射点被剔除在了关联之外。同时,两次观测中也有部分散射点被视为闪烁点而未被关联。这些未关联有效地降低了非平动因素的影响,提高了目标整体平动速度值的估计精度。

图8 第4次观测值与第5次观测值的关联Fig.8 Mapping between No.4 and No.5 observation values

另一方面,从图8中还可以看出，存在一些并不遵循整体运动规律的散射点。为降低其影响,一种办法是计算图8所示的各连接线的斜率值,将其中斜率明显偏差过大的连接剔除,只留下近似平行的轨迹线。而将这些被近似平行关联的散射点的运动参数估计值进行取均值操作,便能获得扩展目标整体在第4分段时间内的平动速度估计。

这种相位分析、散射点关联筛选、估计参数取均值的方式同样也能用于该回波其余各分段的速度估计。为作对比,现同时使用包络对齐并分段进行二次三项式拟合的方法(以下简称包络法)和本文所提的相位分析关联筛选的方法(以下简称相位法)分别对同种分段条件下的上述实测回波进行处理和速度估计,得到前7分段的速度估计结果如图9所示。

图9 飞机速度变化估计曲线Fig.9 Estimation curve of airplane’s velocity changing

从图9中各分段速度估计值的间断情况来看,相位法的估计曲线没有明显的断裂情况,而包络法的曲线则围绕相位法的曲线上下波动。由于目标是实际在连续运动的,其速度的变化曲线在理想的情况下,应该具备连续性。而相位法所估计的曲线较包络法具有更好的连续性,说明在针对该实测数据段的速度值估计上,相位法的效果更佳。

3.3 扩展目标平动速度的估计效果分析

鉴于上述实测数据无法提供目标物的真实平动速度值,为评价速度估计效果,本文先对图5所示的整段回波以滑窗方式进行包络对齐操作,获得如图10(a)所示的对齐结果。再对包络对齐中各脉冲包络的迁移量在慢时间域用最小二乘法进行8次多项式拟合,得到包络走动曲线。最后,对该拟合曲线在慢时间域进行求导,获得如图10(b)所示的速度变化曲线。同时应注意到,该L波段的回波的脉冲重复频率约为250 Hz,故由图10(b)所示速度范围可知本文所用回波存在多普勒模糊情况。

现将图10(b)所示的曲线代替真实的速度变化曲线,与图9所示的两条曲线作对比与分析。先将图9中两条曲线同图10(b)的曲线作差值,得到图11所示曲线。由图11可见,包络法的差值曲线比相位法的差值曲线起伏更大,最大处的偏差值接近0.6 m/s。

图11 两种速度估计曲线与高阶多项式拟合曲线之差Fig.11 Difference between two velocity estimation curves and high-order polynomial fitting curve

为体现估计准确度情况,本文采用根均方误差(root mean square error,RMSE)这一指标来衡量估计值,RMSE的计算公式为

(33)

现将平动速度的估计值再次用于平动补偿中的包络对齐操作,即先用梯形法则将图9和图10(b)中的速度曲线积分为距离曲线,再根据距离曲线对原回波的包络进行移位,得到三幅对齐的包络图。由相位法、包络法、高阶多项式拟合法曲线所得的包络对齐效果如图12所示。

图12 包络对齐效果比较Fig.12 Comparison of envelope alignment effect

再对图12所示的3幅包络对齐效果图,沿慢时间轴进行叠加,获得积累的一维距离像。并根据文献[11]所提及的波形熵概念:

(34)

(35)

分别计算出3条一维距离像的波形熵:高阶多项式拟合曲线作对齐所对应波形熵为17.175 7;包络法曲线作对齐所对应波形熵为17.360 8;相位法曲线作对齐所对应波形熵为17.188 3。可以看出,相位法熵值更接近高阶(8次)多项式拟合熵值,包络法熵值则相对较大。由于波形熵表示一维像的能量集中程度,而包络对齐越精确,一维距离像积累的能量越集中,波形熵越小。因此,相位法的对齐精度接近高阶多项式拟合法,并优于包络法的对齐效果。

值得一提的是,本文使用高阶多项式拟合法所得速度变化曲线替代真实速度曲线进行比较的原因在于:该曲线具备良好的连续性和光滑性,比较符合目标物的实际运动是连续的这一物理属性。然而,当数据点较少时高阶多项式拟合效果并不理想。图13展示了对各个分段分别用8次多项式进行拟合所获得的速度曲线,并与图10(b)的曲线对比。其中,红线为对每一分段分别独立地进行8次多项式拟合所得结果,即每一次拟合都建立在256个数据点上;而黑线为图10(b)所示曲线,其对整段共2 048个数据点直接进行8次多项式拟合。

图13 不同方式下高阶多项式拟合结果对比Fig.13 Comparison of high-order polynomial fitting results in different ways

由图13可见,当数据点较少时高阶多项式拟合法在平动速度估计中会产生较大误差,对数据量的需求使得高阶多项式拟合法并不适用于对雷达信号的即时处理。而本文所提方法在获得一个分段数据时可对目标下一分段所处的距离-多普勒位置进行预测;再通过与下一分段的实测值进行关联处理,进一步提升对该分段速度估计的准确度。这种递进式的速度估计方式更适用于宽带信号的即时处理。

4 结 论

针对雷达扩展目标的平动速度估计,本文提出了一种基于相位分析与散射点关联相结合的方法。该方法通过分段建立匀加速直线运动模型,采用包络对齐与自适应Chirplet分解相结合的方式,能够从相参回波的相位中提取到目标的运动信息;再通过预测-观测匹配关联的方法,去除其中的非平动因素的影响,进一步提高目标速度的估计精度。实测数据的试验结果证实了该方法的有效性。

同时,本文所提方法为一种递进式的扩展目标平动速度估计方法,随着当今运算系统的计算能力不断提升,在未来该算法将更加适用于对宽带雷达信号的即时处理,并将对扩展目标的即时成像与识别起到积极的作用。