利用HRRP序列估计弹道中段目标进动频率

韦楠楠 张兴敢

（南京大学电子科学与工程学院，江苏南京 210023）

1 引言

弹道导弹的飞行过程可分为助推段、中段和再入段，中段是弹道导弹飞行时间最长的一段，被认为是防御的关键阶段。构成弹道中段目标的物体有弹头、诱饵、弹体碎片、金属箔条、干扰机等，这些目标的形状较为简单，弹头和诱饵一般是锥、柱、锥台、球冠的组合体，碎片一般是不规则薄片；这些目标的运动形式也较为特定，是平动与微动（进动、章动、翻滚等）的复合运动。弹头及诱饵等目标由于受到质量及质量分布等因素影响，它们的进动频率、章动角等运动特性存在明显差异，因此常被用于弹道目标识别的研究中。

随着雷达技术的发展，弹道目标的运动细节能够更精细地刻画出来，由于微动能够对电磁波进行额外的调制，导致多普勒频移现象的产生，其表现在时频图上，则是频率周期性的变化，许多学者［1-3］对此问题进行了广泛研究。国内外对弹道目标的微动特性研究主要有两种思路：一是研究学者们［3-6］利用窄带回波信号获得的RCS 序列或时频图提取出目标微动参数，但是目标多散射点的回波信号在一个波束中会相互交叠，从而产生微多普勒折叠的现象。另一种研究思路则是通过分析宽带雷达的一维距离像（high range resolution profile，HRRP）序列的调制特征来提取微动参数，文献［7］利用相邻HRRP之间的差分值形成一维差分序列，再估计出目标进动频率，文献［8］利用逆Radon变换和伪Zernike矩方法提取出HRRP序列的微动特征，文献［9］通过正弦曲线拟合的方法提取出目标散射中心的微动参数。宽带雷达的分辨率较高，可以在距离维上分辨出目标对应的多个散射点，同时还能解决多分量微多普勒信号混叠问题。无论是利用窄带信号或是宽带信号提取微动特征，由于诱饵与弹头的相似度越来越高，并且反导作战对系统实时性要求越来越高，而微动特征提取与目标成像往往依赖于长时间观测，因此如何在数据率、积累时长、带宽等条件的限制下提取出稳定的、可分性高的微动特征是需要考虑的问题。

本文的研究则是通过宽带雷达获得的HRRP序列提取目标进动频率。本文以锥柱状目标作为分析对象，首先通过电磁仿真获得锥柱状目标的转台特性，分析其多个散射点在半空间范围内的可见范围；其次，通过构建锥柱状目标的进动模型分析HRRP 序列周期性变化的过程；建立锥柱体弹头滑动型散射中心的位置模型，在此基础上推导滑动型散射中心的微多普勒的数学表达式。基于以上分析，本文提出了自适应阈值法获取目标进动频率的算法，该算法提取HRRP的强散射中心单元，对该单元进行时频处理获得时频图，通过自适应阈值法提取出时频图多个散射点的时频脊线，据此估计目标微动频率。本文实验采用了卫星工具箱（STK）和电磁仿真工具（FEKO）联合仿真，生成弹道导弹宽带雷达回波数据，通过蒙特卡洛实验和对比实验证明了算法具有较好的稳定性和抗噪性。

2 弹道目标进动模型及散射特性

2.1 锥柱体目标的散射特性

旋转体目标的散射场，主要包括尖顶散射、锥面行波、曲面不连续处散射、柱面行波和角顶散射，这五类散射中，行波是较弱的一类散射场。雷达获得的回波信号可以描述为各散射中心的回波合成，当雷达发射电磁波照射在目标上，目标旋转一定角度后，如果电磁波照亮部分的形状、介质等影响回波因素未发生变化，那么目标的散射场也基本不变。

建立锥柱状弹头模型如图1 所示，雷达发射带宽为8.5 GHz～9.5 GHz 的步进频信号，通过FEKO计算得到的目标一维距离像在0°～180°范围内的变化情况如图2 所示。五个散射中心的可见范围如下：散射中心A为0～(π -γ)，散射中心B、C在(0～π)范围内均可见，散射中心D 在范围内可见，在散射中心E的可见范围为(0～γ)。通常防御雷达位于弹头的前下方，雷达视线和目标主轴的夹角为50°～80°。

2.2 锥柱体目标的进动模型分析

建立锥柱体弹头目标的进动模型如图3 所示，假设以目标质心O 点为坐标原点建立参考坐标系OXYZ。目标围绕自身对称轴OA 旋转的运动称为自旋，自旋角速度设为ωZ；目标绕着某一旋转轴（非几何对称轴）旋转的运动称为锥旋，锥旋轴设为OZ，目标锥旋角速度设为ωp；OZ 与OA 的夹角θ，称为进动角。假设雷达视线（radar light of sight，RLOS）方向位于YOZ 平面内，与OZ 的夹角α称为雷达视线角；RLOS 与OA 轴的夹角ρ称为姿态角。假设在t=0时刻，自旋轴在YOZ平面的初始方位角为φ0，则姿态角ρ随时间的变化公式为

式（1）中雷达视线角α会随着弹道目标飞行时间的推移而改变，弹道导弹中段的飞行时间一般为十几分钟，而进动周期为秒级，α在一个进动周期内的变化几乎可以忽略，所以在求进动参数时α可看作固定值。进动角θ对ρ(t)的影响与α等效。式（1）表明了目标进动的周期性会导致姿态角ρ周期性变化。对于旋转对称体而言，自旋不会导致散射场变化，但锥旋和雷达视线的变化会影响散射场，所以目标周期性进动会导致一维距离像上的散射中心位置和幅度发生周期性变化。

2.3 锥柱体目标散射中心的微多普勒模型

假设锥柱体目标平动经过补偿，只考虑进动带来的频率调制。进动目标的散射场变化的影响因素主要是锥旋，因此仅考虑目标锥旋运动产生的微多普勒特性。目标模型如图3 所示，假设目标质心点O 到锥柱结合的平面的距离为d1，与柱体底面的距离为d2，柱体底面半径为r。目标对称轴OA 在XOY 平面内的投影为OV，OV 与OX 的夹角为φ，在弹体上取一个固连其上的坐标系OVWZ，OV、OW、OZ 构成右手系。记参考坐标系OXYZ 的单位向量基为[eX，eY，eZ]，连体坐标系OVWZ 的单位向量基为[eV，eW，eZ]，则目标体锥旋角速度向量可定义为：ω=φ⋅eZ=ωp⋅eZ，其中φ=φ0+ωpt，φ0为OA 的初始位置在XOY 平面内的投影与OX 轴的夹角，由于φ0仅改变微多普勒的初相，不影响微多普勒的周期及幅度，所以令φ0=0，则φ=ωpt。

锥柱体目标的散射中心可分为锥顶局部散射中心（图1中的点A）、边缘（棱线）型散射中心（图1中的B、C、D、E）。根据电磁散射理论，弹头顶部A通常会形成一个较强的散射中心，由于该散射中心位于弹体对称轴上，其运动规律和弹头运动规律相同，即以角速度ωp绕OZ做圆周运动。根据文献［10］的分析，锥顶散射点A 在OVWZ 坐标系的矢量rA=(sinθ⋅eV+cosθ⋅eW)⋅LA，LA为OA长度，则点A的速度可表示为

雷达视线方向为RLOS=[0，sinα，cosα]T，雷达发射电磁波的波长为λ，那么锥顶散射中心A 在雷达视线方向上引起的微多普勒可表示为：

由式（3）可以看出，锥顶散射中心A 的微多普勒表达式是正弦形式。

边缘（棱线）型散射中心的微多普勒，除了受弹体微运动的影响外，还和雷达视线与锥轴所构成的平面有关，当旋转对称弹头在自由空间进动时，散射中心B、C、D、E 的微运动规律和弹头本身的进动规律并不完全一致，因此，该处散射中心的微多普勒将不再是简单的正弦形式。以散射中心B、E 为例，根据文献［11-12］的分析，散射中心B、E 在雷达视线方向上引起的微多普勒可表示为：

其中：F(t) ≜sinθsinαsinωpt+cosθcosα，i=B，E。从式（4）可以看出，滑动型散射中心的微动引起的微多普勒在余弦形式之外，增加了的调制项。对式（4）进行泰勒展开，并忽略较小的高阶项，可以得到散射中心B、E的微多普勒近似表达式为

将式（5）化简为式（6）形式

底部边缘形成的散射中心C、D 的微多普勒形式也如上式（4）所给出，只是d1替换为坐标原点O 到底部所在平面的距离d2，同时只取加号即可。从式（5）、（6）来看，滑动型散射模型的微多普勒呈现显著的非单频性，但其微多普勒曲线表现出周期性，这种周期性与目标的锥旋频率相对应。滑动型散射中心的微多普勒在锥旋频率的基础上增加了倍频分量，但并未改变微多普勒的周期性。

以图1目标模型进行仿真分析，利用式（1）生成不同雷达视线角下的目标姿态角数据，并结合目标电磁仿真数据，产生不同视线角下的目标HRRP 序列。设目标进动频率为1.5 Hz，进动角度为10°，雷达工作参数同2.1，脉冲重复频率为500 Hz，观察时间设为2 s。时频分析方法采用100 点的汉明窗和128点的短时傅里叶变化处理。

图4、图5给出了目标分别在雷达视线角为50°、80°下的微多普勒时频分析结果。图4（a）的HRRP序列显示50°视线角下目标散射点出现越距离单元走动，而图5（a）显示的80°视线角下的目标散射点出现重合，这一现象与图2所显示的半空间HRRP散射点变化一致，由于目标进动的周期性，一维距离像的幅度与能量强弱产生周期性变化。各自取不同距离单元的回波进行时频分析，当雷达视线角为50°时，图4 显示目标回波在第32、33 个距离单元的时频图可以观察到区分明显的时频脊线，而在第34、35、36个距离单元的时频图出现了不同程度的交叠与重合；图5 显示目标在80°的雷达视线角下，其回波第32个距离单元的时频图出现严重交叠现象。

旋转体目标的进动导致各散射点的微动变化曲线为正弦函数或正弦函数的叠加形式，而在某些视角下，可能会出现散射点无法分离、时频脊线交叠重合的现象，但我们发现时频脊线的振幅与能量强弱均具有周期性，该周期性与进动频率相对应。在本文的研究中，试图从图像处理技术的角度从时频图中提取出不同散射点的时频脊线，从而估计出目标的进动频率。

3 基于HRRP序列的进动频率估计方法

3.1 提取HRRP的强散射中心单元

假设雷达回波信号经过平动补偿，收到Np帧HRRP 序列，其中每幅距离像的距离单元数Ns。为进行进动频率的提取，将其构建为一维距离像数据矩阵：

3.2 微多普勒瞬时频率估计

进动频率提取和分离的前提，首先通过对H0r的Nst个强散射中心单元依次进行STFT处理，得到散射点的时频分布图，由于HRRP在经过平动补偿后，时频曲线会聚焦在零频附近，在时频曲线交点处，容易产生错误关联而无法分离各个曲线。在这里提出利用自适应阈值的方法从时频图中提取时频脊线。

STFT的定义为：

x(t)为待处理信号，g(t)为窗函数。假设时频图中存在p条时频脊线，通过遍历每一时刻的频谱提取时频脊线，设t时刻的频谱为：

Ei为第i个频点的幅值，ωi为第i个频点的频率值。首先对ft(i)进行归一化：操作步骤如下：

步骤1首先，设立能量幅度阈值e为0.15，判断满足ft(ωi)＞e条件的频点值个数，假设满足该条件的频点值有k个，即ωi1，ωi2，...，ωik；若k＞p，则令e=e+0.01，通过提高阈值，去掉旁瓣；若k=0，则令e=e-0.005；若k＜p，则令ωi[(k+1)：p]=ωik，直至k=p，则退出循环。此时，在t时刻可提取出p个频点值。

步骤2提取完所有时刻点的频谱值后，此时获得一个p×m维矩阵表示第k个散射点，在第tm时刻提取的频率值大小。由于目标短时间运动过程中，可认为回波信号是平稳的，其时频变化曲线的幅度大小相对平稳。利用此特点，对Wp×m矩阵进行重排，每一列按照频率幅值由小到大排序，，重排后的矩阵为。

步骤3利用卡尔曼滤波方法对的每一条时频脊线进行平滑处理。

需要说明的是，步骤1 中的阈值e的设定没有严格的要求，由于能量幅度已经过归一化处理，因此需要满足e∈(0，1)，该步骤的最终目的是提取出p个频点值，e的值越大，能提取到的频点值越少，e的值越小，则提取到的频点值越精确。

3.3 进动频率估计

在2.3 节的分析中，进动时目标的微多普勒频率曲线为周期函数，其频率分量一定是进动频率的倍数，因此这为进动频率的提取提供了途径。该曲线的振动频率等于进动频率。通过3.2 的三个步骤，提取了散射中心的时频线Xk(tm)

其中Xk(tm)为第k个散射点的微多普勒曲线。时频脊线去除直流分量的影响，利用快速傅里叶变换分别对p条时频脊线曲线进行频谱分析，

获得相应的频谱曲线Fk(fn)，fn=f1，...，fN，fN为采样频率，取能量较大的频率分量作为进动频率估计值

为第k个散射点进动频率估计值。最后取p条时频曲线的进动频率的平均值作为最后估计的进动频率为：

在3.1 的预处理中，提取出了Nst个强散射中心单元，对每个强散射中心单元重复以上步骤，获得Nst个进动频率估计值，对其求均值，得到最终的频率估计值。

4 仿真数据的分析与处理

实验数据由仿真获得，仿真流程如图6所示，首先采用卫星工具箱（STK）生成弹道数据，再利用电磁散射仿真工具（FEKO）生成电磁仿真数据。

STK软件负责对战场环境、弹头平动微动、地球自转、引力等因素进行综合仿真，并生成任意时刻下雷达视线角、飞行距离、速度等数据。实验所用弹道发射地经纬度为(43.32，63.11)，落点经纬度为(116.40，39.89)，关机点位置(82.76，60.52)，取弹道中段的飞行时间进行实验，观测时间总长为297 s。在观测段中雷达视线与导弹运动轨迹切线方向夹角（即雷达视线角）变化如图7所示。

目标散射模型由FEKO 软件进行仿真。电磁仿真参数设置如下：雷达信号为X波段的步进频信号，信号载频起始频率为8.5 GHz，终止频率为9.5 GHz，工作频率步长为15.675 MHz，电磁计算的求解方式为PO算法。目标尺寸如图1所示，目标运动参数见表1，三个目标的目标尺寸相同，微动参数不同。

表1 目标微动参数Tab.1 Simulation parameters

通过FEKO 软件计算生成电磁数据包含在.ffe文件中，其包括了各入射点和反射点的散射电场E(n)。目标频域响应序列Sr(t，n)与E(n)的对应关系如式（17）所示为：

对序列Sr(t，n)进行逆傅里叶变换即获得HRRP。

4.1 进动频率估计

图8 是雷达前1000 次回波得到的弹头及诱饵的一维距离像序列。提取弹头和诱饵的一维距离像序列的强散射中心单元，并应用STFT 获得时频分布图如图9 所示，从图中可观察到时频脊线的幅度强弱存在周期性变化，并且也不是严格意义上的正弦曲线。

接下来对时频图提取时频脊线，效果如图10所示，从提取效果来看，本文所提的数据关联的方法，能够准确分离出时频脊线。

对获得的时频脊线去直流分量后分别进行FFT，如图11 所示，T1 目标的P1、P2散射点时频脊线的频谱最大值分别为1.4648 Hz、1.4648 Hz，T2 目标的P1、P2散射点时频脊线的FFT 频率最大值分别为1.9531 Hz、1.9531 Hz。

最后获得T1 的估计值为1.4726 Hz，T2 的估计值为2.0215 Hz，与真实值误差较小。

4.2 性能分析

为了验证本文算法在不同噪声影响下的稳定性，本文对仿真的信号加入高斯白噪声，信噪比从SNR=-5 以1 dB 的步进递增到20 dB，不同信噪比下进行100次蒙特卡洛实验获得进动参数的估计误差。取T2目标进行进动频率估计实验，雷达脉冲重复频率设为500 Hz，回波积累时长为2 s。结果如图12所示（图12（b）的均方误差结果利用20 log（）处理成dB 单位），从结果来看，当信噪比大于1 dB 时，进动频率的估计均方根误差趋于平稳，与文献［7］提出的序列差分方法进行比较，本文算法在信噪比较低的情况下，依然具有良好的提取效果。

4.3 弹道中段目标进动特征提取

为了验证本文算法在弹道目标整个观测段的提取效果，对弹道回波叠加信噪比为12 dB 的高斯白噪声，回波积累窗长为2 s，窗口步进1 s，共有297 个统计窗口。雷达工作方式不变，PRF 设为500 Hz，目标参数不变，对每个统计窗中的HRRP 序列提取目标进动频率，最后对297个估计值求平均，结果如表2 所示，三个目标的估计结果均靠近真实值。

表2 弹道中段目标进动频率估计结果Tab.2 Precession frequency estimation of ballistic missiles in midcourse phase

5 结论

本文通过分析弹道导弹的锥柱体目标的进动模型与散射点模型，利用散射点在HRRP 序列上的周期特性，运用时频分析方法获得时频图，提取目标的微多普勒频率，并估计目标进动参数。仿真结果表明，该方法利用图像处理方式能够在一定程度上消除噪声的影响，并且能解决不同视角下时频脊线重合的问题，有效地估计出进动参数，具有较好的稳定性和抗噪性。

中段弹头目标结构一般较为简单，散射中心位置要满足在观测期间有一个强散射中心的条件。文中实验只分析了两个散射中心的情况，实际中，目标可能存在两个或两个以上的散射中心，本文所提的方法依然适用，并且散射点对应位置关联的方法也相对简单。