1 引言

随着诱饵技术和电磁控制技术等的不断发展，基于传统特征量的弹道目标识别技术已难以满足未来高科技战争的需求，而微动特征作为弹道目标的固有属性，难以被模仿，且弹头和诱饵、碎片等运动形式存在明显差异，因而被用来区分识别真弹头，近年来受到国内外研究学者的广泛关注^[1–3]。由于目标微动会对雷达回波产生微多普勒调制，无论是宽带雷达还是窄带雷达，均可通过对回波信号参数进行估计来提取目标的微动特征，进而进行目标的分类识别。然而研究表明^[4–6]，单基雷达获取的目标微多普勒信息均具有较强的姿态敏感性，不同视角获得的微动特征不同，因而难以准确反映目标的真实空间结构。考虑到组网雷达能获得目标在不同视角上的观测信息，于是有望利用多部雷达联合观测对目标进行3维微动特征提取，从而克服单一视角的局限性，提高目标识别精度。

考虑宽带雷达能得到目标的高分辨距离像，获取目标更细微的形状结构特征，目前基于宽带雷达组网的目标微动特征提取研究较多。文献[7]对MIMO雷达中的旋转目标微多普勒效应展开分析，并基于时间-距离像对目标的3维微动特征进行了提取；文献[8]利用分布式组网雷达对有翼锥体目标进行联合观测，并基于不同视角微多普勒特征的相关性，借助几何分析的方法，实现了目标的3维进动特征提取；文献[9]通过三站1维距离像融合实现了锥体目标的3维重构。由此可见，基于宽带雷达组网的目标微动特征提取技术研究已较为成熟，然而，就目前现有雷达实际装备情况来看，由于宽带高分辨雷达价格昂贵，一时间仍难以实现对窄带雷达的全面升级和替换，雷达网仍以窄带体制为主，因此研究利用窄带雷达网对弹道目标进行特征提取具有更大的实际意义，能有效解决目前理论研究与实际脱节的矛盾；文献[10]在窄带组网体制下推导了不同视角锥体散射中心瞬时频率变化关系，利用频谱熵实现了散射中心的匹配关联，并进一步提出了基于散射中心瞬时频率相关性的目标参数提取算法，但没能获取目标的3维微动特征，也无法实现空间目标的3维重构。

本文在文献[10]的基础上，进一步对窄带雷达网中的锥体目标3维进动特征提取展开了研究。在详细分析了锥体进动引发的微多普勒频率调制特性的基础上，利用锥顶微多普勒频率调制系数比，实现了不同视角下散射中心匹配关联，并获取了目标的3维锥旋矢量，进而利用锥顶和底面边缘散射中心微多普勒频率相关性，结合频率补偿的方法对锥体特征参数进行了提取，在此基础上解算出每一时刻锥顶坐标，从而实现了目标空间位置的3维重构。最后，仿真分析了本文方法的有效性。

2 进动模型

以无翼锥形弹头为例，建立进动模型如图 1所示，由于其存在旋转对称性，因此仅考虑其做锥旋运动。假定目标以角速度ω_c绕锥旋轴旋转，且锥旋轴与目标对称轴相交于点o，两者夹角为θ，以交点为坐标原点建立参考坐标系oxyz如下，目标对称轴初始方位为φ₀, LOS表示雷达视线方向，其在参考坐标系中的方位角和俯仰角分别为 (ε, χ)，与锥旋轴的夹角为α，与对称轴夹角为β，定义雷达视线方向与对称轴构成的平面为底面圆环的电磁波入射平面，该平面与圆环交于p、q两点。锥体高度为h，底面半径为r，锥顶与进动中心的距离为h₁，底面中心与进动中心的距离为h₂，且目标满足远场条件，雷达与进动中心的距离为R₀。

根据散射中心理论，对于旋转对称目标，其高频散射特性主要由锥顶D及底面边缘两个散射中心p, q确定^{[1, 5]}。由几何关系分析可得雷达视线与对称轴的夹角β满足：

式中，φ为初始相位角，且由文献[6]可知是$\varphi ={{\phi }_{0}}-\varepsilon $。进一步对各散射中心在雷达视线上的投影关系分析可知，锥体3个散射中心到雷达的距离分别为：

考虑到锥体目标在实际运动中各部分之间存在相互遮挡，目标上各散射中心不能始终保持同时可见，使得式 (2) 的使用范围受到限制。但锥顶D和近散射点p在大部分情况下都能被观测到^[3]，并能够获得二者的稳定连续观测信息，因此本文主要利用D, p的微动信息展开后继研究。假设雷达波长为λ，由式 (1) 和式 (2) 可得，进动引发的D, p两点的微多普勒调制为^[11]：

由式 (3) 可以看出，D点的微多普勒频率变化服从正弦规律，而p点的微多普勒频率由两部分之和组成，不再服从简单的正弦调制规律，且两点微多普勒频率均与目标的进动和结构特征有关，共包含w_c, h₁, h₂, r, θ, α, φ 7个未知参数，其中ω_c, φ可通过提取正弦曲线特征得到，而θ, α两者之间存在耦合，仅通过单部雷达，仍无法实现对目标进动角及尺寸大小的求解。考虑到多视角观测能获得更加丰富的目标信息，具有较好的解耦合性能，因此本文将采用雷达组网方式对目标特征进行提取，并进一步实现3维重构。

3 散射中心匹配关联

首先建立窄带雷达网系统观测模型如图 2所示，图中OXYZ为全局坐标系，与参考坐标系oxyz平行，假定系统中共有N部窄带雷达同时进行观测，并都已满足时空同步要求，各雷达视线在OXYZ坐标系中的方位角和俯仰角为 (ε_i, χ_i), n_i (i=1, 2, ···, N) 为雷达视线方向，满足：

由第2节分析可知，当采用多部雷达同时进行观测时，由于各雷达观测视角不同，同一时刻目标各散射中心在雷达视线上的投影位置排列顺序将存在差异，相对应地，同一时刻各散射中心的微多普勒频率也会不同。因此，在利用组网雷达进行特征提取之前，首先得实现不同视角散射中心的匹配关联。

由式 (3) 可以看出，对于同一观测目标而言，锥顶D的微多普勒频率调制系数A仅与雷达观测视角有关，任取雷达网中两部雷达，其调制系数比满足：

而p点调制规律更为复杂，不具备上述比例关系，因此，通过比较观察不同雷达间的调制系数比即可实现散射中心的匹配关联。文献[10]采用频谱分析的方法，通过计算不同散射中心的频谱熵来实现散射中心的匹配关联，然而在两个散射中心回波信号无法分离的情况下，散射中心无法与各自频谱一一对应起来，因此该方法存在较大的局限性。相比而言，本文方法则更加简单实用。

4 3维重构 4.1 3维锥旋矢量提取

为更好地实现对锥体目标的3维重构，首先对锥旋轴方向进行估计。采用Viterbi算法提取锥顶微多普勒曲线振幅得到：

Viterbi算法^[12–14]作为信号隐状态估计的有效手段之一，能够依据各信号成分强度对信号进行逐次分离，因此常被用来对多目标信号瞬时频率进行估计。与逆Randon变换、Hough变换等曲线参数提取方法相比，Viterbi算法对各信号分量形式依赖性不高，即使是非正弦信号，同样能实现瞬时频率的准确估计。

此外，雷达观测视角α_i还满足：

式 (7) 中，$\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\omega }\text{ = }{{\left( {{\omega }_{x}}, {{\omega }_{y}}, {{\omega }_{z}} \right)}^{\text{T}}}$，联立式 (4)、式 (6)、式 (7)，令B=h₁sinθ，此时方程组中共包含B, ω_x, ω_y, ω_z 4个未知参数，因此，至少需要3部雷达同时进行观测才能实现对上述参数的求解。进一步将求得的参数回代到方程组，还可以确定sina_i的大小。文献[10]采用基于视线角方差最小准则的频率搜索补偿方法对雷达视线角进行估计，需要经过多次循环迭代才能得到准确估计值，计算复杂，且易受噪声影响，而本文所提方法用到的只是锥顶微多普勒频率曲线的振幅和周期，且这两个曲线参数均可由Viterbi算法准确提取得到，在求得锥旋矢量的同时也能估计出各雷达视线角的大小，计算更为简单，算法稳定性更好。

4.2 锥体弹头参数提取

在上述分析的基础上，若要提取锥体弹头参数，还需p对点的微多普勒频率进行充分利用，观察式 (3) 可知，f_d–p由正弦部分和非正弦部分之和组成，且正弦部分${{f}_{k}}=-2{{\omega }_{\text{c}}}{{h}_{2}}\sin \theta \sin \alpha \cos ({{\omega }_{\text{c}}}t+\varphi )\text{ }/\text{ }\!\!\lambda\!\!\text{ }$满足f_k=(-h₂/h₁)f_d-D，而非正弦部分此时仅包含 (r, θ) 两个未知参数。若能将正弦部分完全补偿，便可利用多视角观测对 (r, θ) 联立求解。考虑到锥体目标尺寸信息仍然未知，先假设补偿系数为η，且η₀=h₂/h₁，当η=η₀时，便可实现完全补偿，于是p点补偿后的微多普勒频率满足：

对η进行遍历，利用补偿后的频率两两联立方程可求得：

对每个η取值所对应求得的所有结果$\left( {{{\hat{r}}}_{\xi }}\left| \eta, {{{\hat{\theta }}}_{\xi }}\left| \eta \right. \right. \right)$做进一步处理，并定义归一化标准差σ为：

式中，$\Delta r=\left[{{{\hat{r}}}_{1}}-\bar{r}\ \ {{{\hat{r}}}_{2}}-\bar{r}\cdots {{{\hat{r}}}_{C_{N}^{2}}}-\bar{r} \right], \Delta \theta =\left[{{{\hat{\theta }}}_{1}}-\bar{\theta }\ \ {{{\hat{\theta }}}_{2}}-\bar{\theta }\cdots {{{\hat{\theta }}}_{C_{N}^{2}}}-\bar{\theta } \right], \bar{r}, \bar{\theta }$为平均值，按照上述归一化标准差定义，对于每一个η取值均能得到对应的σ。若η=η₁时，σ取得最小值，则说明此时f_k被补偿得最完全，补偿系数η₁也越接近η₀，由此可求得：

结合4.1节分析，将${\hat{\theta }} $代入B=h₁sinθ中，于是求得${{h}_{1}}=B\sin \hat{\theta }, {{h}_{2}}={{h}_{1}}\cdot {{\eta }_{1}}$。

4.3 锥顶空间位置估计

在求得锥体目标结构参数及旋转轴方向的基础上，为实现对目标空间位置的3维重构，还需确定各散射中心的相对位置，由于底面边缘两个散射中心会随雷达视线方向改变产生滑动，位置坐标不易确定，因而本文从锥顶散射中心入手，在锥体目标结构参数已知的条件下，只要能够求得每一时刻锥顶坐标，同样能实现对目标空间位置的3维重构。由于窄带雷达距离分辨力较低，难以直接从目标回波中获得各散射中心的径向距离变化规律，因此本文考虑在已知各散射中心运动形式和参数基础上，通过微多普勒频率反推每一时刻各散射中心相对应的径向距离变化。

令d=R_D–R₀，由式 (2) 可得：

式中，${{\varphi }_{i}}={{\phi }_{0}}-{{\chi }_{i}}$，由于每一时刻D点的微多普勒频率f_d–D均已获得，且f_d–D与R_D满足导数关系，因而d_i在每一时刻的值也能求解得到。此外，结合图 1可知，d_i为oD在第i部雷达视线上的投影，同时还应满足：

若令oD=(D_x, D_y, D_z)^T，通过3部雷达同时进行观测可以解算出oD为：

$ \rm{\pmb{OD}} = {\left[\begin{gathered} \mathit{\boldsymbol{n}}_1^{\rm{T}} \hfill \\ \mathit{\boldsymbol{n}}_2^{\rm{T}} \hfill \\ \mathit{\boldsymbol{n}}_3^{\rm{T}} \hfill \\ \end{gathered} \right]^{ - 1}} \cdot \left[\begin{gathered} {d_1} \hfill \\ {d_2} \hfill \\ {d_3} \hfill \\ \end{gathered} \right] $

(14)

综上所述，基于窄带雷达组网的弹道目标3维微动特征提取及重构步骤为：

步骤1建立弹道目标进动模型，分析各散射中心微多普勒调制规律；

步骤2对目标回波进行时频分析，采用Viterbi算法提取各散射中心微多普勒曲线；

步骤3基于锥顶微多普勒频率调制系数比，实现不同视角下散射中心匹配关联；

步骤4提取3维锥旋矢量和目标结构参数，在此基础上解算出每一时刻锥顶坐标，从而实现锥体目标空间3维重构。

5 仿真分析 5.1 目标参数估计

在下述仿真中设定目标为锥体，目标参数设置为：h₁=2.0 m, h₂=0.5 m, r=0.5 m, h=2.5 m, θ=13°，目标对称轴初始方位角φ₀=60°，目标的锥旋频率为f_c=4 Hz，锥旋矢量为$\left( 2\sqrt{3}\text{ }\!\!\pi\!\!\text{, 4}\sqrt{3}\text{ }\!\!\pi\!\!\text{ }, 2\text{ }\!\!\pi\!\!\text{ } \right) $。雷达参数设置为：载频f=8×10⁹ Hz，信号带宽为5 MHz，雷达脉冲重复频率为2000 Hz，积累时间1 s，信噪比为10 dB。在全局坐标系中3部雷达M₁, M₂, M₃测得的目标方位角和俯仰角 (ε_i, χ_i) 分别为 (40°, 84°), (45°, 48°), (30°, 17°)。图 3分别为该3部雷达获得的同一时间段内目标回波Cohen类时频分布重排结果，可以看出重排后的谱图不仅具有更好的时频聚集性，同时还有效抑制了各分量之间的交叉项^{[15, 16]}，有利于提高各散射中心瞬时频率的估计精度。

在上述谱图重排的基础上，进一步采用Viterbi算法对各散射中心瞬时频率进行提取，并通过拟合更好地削减了交叉项所带来的不利影响，得到各雷达散射中心瞬时频率估计结果如图 4所示。然后依据第3节匹配关联准则，对各散射中心曲线幅度做进一步处理，区分出锥顶和锥底边缘散射中心。此时提取到各雷达锥顶正弦曲线的频率为$\left\| {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\omega }} \right\|$ =25.14 rad/s，振幅分别为567.6 Hz, 361.8 Hz, 307.8 Hz，而由式 (6) 计算得到的振幅理论值分别为569.2636 Hz, 362.0643 Hz, 309.5390 Hz，两者相当接近，代入方程组式 (7) 可以求解得到3维锥旋矢量$\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\omega }={{\left( 10.8598, 21.7937, 6.2545 \right)}^{\text{T}}}$，与理论值基本吻合，同时可以求得各雷达视线角分别为${{{\hat{\alpha }}}_{1}}=70.7025, {{{\hat{\alpha }}}_{2}}=36.9851, {{{\hat{\alpha }}}_{3}}=30.7848$。

依据4.2节中锥体弹头参数提取算法，由求得的${{{\hat{\alpha }}}_{i}}$和各散射中心瞬时频率对 (r, θ) 进行估计，可以得到遍历时η归一化标准差σ随之变化的结果如图 5所示，当η=0.253时，σ取得最小值，接近于理论分析值η₀=0.25，进一步将该η值代到方程式 (11)，从而求得参数r, θ, h₁, h₂估计值。定义相对误差=|理论值–估计值|/理论值，各参数估计结果如表 1所示，其中α₁, α₂, α₃估计精度相对较高，r, θ, h₁, h₂由于受到各散射中心瞬时频率提取误差的影响，估计精度则相对较低，通过进一步提高时频分辨率或者增加雷达网中雷达观测数量可实现上述参数估计精度的提升。但总的来说，各参数估计相对误差均小于5%，满足目标识别的精度要求，可用于下一步对锥顶坐标的求解。而在相同仿真条件下，文献[10]中的目标参数平均估计相对误差则接近于15%，明显高于本文误差，这也从侧面反映出本文参数提取算法的稳定性。

在上述目标特征参数提取的基础上，进一步按照4.3节所提算法对锥顶坐标进行求解。当时t=0.25 s，求得锥顶坐标$\left( {{{\hat{D}}}_{x}}, {{{\hat{D}}}_{y}}, {{{\hat{D}}}_{z}} \right)=\left( 1.1739, 1.1888, 0.5898 \right) $，与理论值$ \left( {{D}_{x}}, {{D}_{y}}, {{D}_{z}} \right)=\left( 1.2050, 1.2495, 0.6093 \right)$基本相符，并最终得到在观测时间内0~0.25 s锥顶散射中心的实际轨迹如图 6所示，与其理论轨迹近乎重合，从而更加充分地说明了本文重构方法的准确性和有效性。

5.2 算法鲁棒性分析

为了充分验证本文所提算法的鲁棒性，仿真分析了曲线参数估计误差对目标参数提取及重构精度的影响。由于在4.1节中3维锥旋矢量的求解精度主要受到锥顶微多普勒幅度和微动周期的影响，而微动周期通常都能被准确估计，因而在此主要分析微多普勒幅度的影响。在4.2节中，目标参数的准确提取关键在于式 (9) 的求解，且求解精度主要受瞬时频率提取误差制约，因此也有必要对瞬时频率提取误差带来的影响进行分析。为便于分析，定义归一化误差如下：

上式中$ {\hat{X}}$为估计值，X为真实值，进一步定义|η|为归一化绝对误差。假设锥顶微多普勒曲线幅度A_D提取值以及各散射中心瞬时频率的提取值归一化误差服从[–a a]上的均匀分布，采用蒙特卡洛方法进行分析，仿真100次，可以得到目标特征参数归一化绝对误差平均值的变化如图 7所示。从图 7(a)中可以看出，当a在区间[0, 0.1]变化时，锥旋矢量各方向分量估计误差均呈线性增加趋势，但总的来说，各方向分量估计精度仍然较高。而在图 7(b)中，当a在区间[0, 0.03]变化时，进动角θ变化较为平和，h₂/h₁则几乎不受影响，只有底面半径r值的估计精度变化最为敏感，随着各种高性能时频分析工具的发展，完全可以将瞬时频率提取精度进一步提高，从而满足目标参数的高精度提取。

综上所述可知，本文所提锥体目标参数提取及重构算法在一定程度上受到微多普勒曲线参数提取精度的影响，但由于文中所采用的基于Cohen类时频重排的Viterbi算法较好地实现瞬时频率曲线的提取，因此本文算法能保证目标参数提取及重构的可靠性，可用于目标识别。

6 结束语

本文对基于窄带雷达网的锥体目标3维进动特征提取问题展开了研究。依据目标的多视角微多普勒频率调制特性，利用3部雷达获取了目标的3维锥旋矢量及特征参数，并进一步解算出每一时刻锥顶坐标，实现了目标空间位置的3维重构。仿真结果表明，本文所提方法目标参数估计精度高，重构性能好，能够有效克服目标散射中心遮挡和姿态敏感性的不利影响，为基于窄带雷达的空间目标准确识别提供了解决方案。考虑到在实际的导弹防御体系当中，可能会出现窄带雷达和宽带雷达同时对目标进行观测识别的情况，后继工作将就不同体制雷达对目标3维特征参数的融合提取问题展开研究。