炮弹起始扰动姿轨耦合建模与微多普勒特征分析

宫志华 李开明 段鹏伟 陈春江

①(中国人民解放军63861部队 白城 137001)

②(空军工程大学信息与导航学院 西安 710077)

1 引言

依据弹道学理论，炮弹目标运动包括质心运动和绕心转动，据此构建的炮弹运动模型就是炮弹质心运动和绕心转动的动力学方程[1-4]。炮弹绕心运动，亦称为角运动，包括自旋、进动、章动、翻滚和摆动等，在武器试验鉴定领域，这些角运动表征参数的获取对武器系统的战术性能研判、故障判定、射表编拟和目标辨识等具有重要意义。

2000年，美国海军研究院的Chen[5]将目标质心平动以外的旋转、振动等微小运动称为微动(micromotion)，并将微动对雷达回波产生的附加频率调制称为微多普勒效应(micro-Doppler effect)。在高频区，目标回波可被等效为一系列散射中心电磁散射的合成，通过各散射中心在空间的变化可近似反映目标的运动情况[6]。炮弹角运动使弹体上某些强散射中心发生空间变化，对雷达回波信号形成的频率调制，也属于典型的微多普勒效应，其产生的微多普勒特征是炮弹目标状态判定的重要依据。因此，炮弹角运动的动力学模型是炮弹目标角运动参数与测量雷达回波信号映射模型的理论基础，是描述炮弹角运动特征量与雷达信号之间关系的依据，可以反映炮弹角运动规律与目标回波之间的内在联系，是炮弹目标特征参数辨识的重要参考。

近年来，雷达目标微多普勒效应研究是国内外雷达领域学术界和工程界的研究热点之一[7-12]。文献[10]提出一种基于物理光学法和等效电磁流法的扩展性微动目标回波模拟方法，并分析了圆锥与带翼弹头的进动特性；文献[11]针对相位污染的缺损回波提取了空间碎片目标的微动特征；文献[12]基于正交频分复用信号外辐射源雷达，对带旋翼低空目标进行回波建模与微动参数提取；文献[13]给出了飞行动目标及其微动电磁散射特征分析的系统方法，借助软件建模平台建立了目标多自由度运动模型，对目标的雷达回波信号进行了仿真计算分析。实际应用场景中，对微多普勒效应进行建模分析时结合目标环境与运动特征能够获得更好的解释性，从而对目标特性分析和识别更有价值[14]。同时，现有的对弹道目标的微动回波建模大多针对外太空飞行的弹道导弹目标展开[15]，弹道导弹目标的弹道中段通常在平流层，受大气影响可以忽略，而对于炮弹、火箭弹等战术类导弹、弹丸目标，其飞行基本在大气层内，受空气扰动影响明显。因此，其动力学模型与弹道导弹存在明显差异，常规的弹道导弹微动回波建模不再适用[1-4,6]。特别是，炮弹目标在弹道初始直线段受到起始扰动的影响[16]，其特殊的角运动对目标回波的调制作用与弹道导弹目标微动回波存在明显差异。

基于以上分析，本文针对炮弹类目标，重点分析炮弹目标角运动机理，在此基础上结合目标运动学方程建立目标角运动回波信号模型，从炮弹目标运动的物理原理出发，揭示目标微动与回波信号之间的内在关系，为有效分析炮弹目标的雷达测量信息、准确辨识炮弹目标运动参数提供参考依据。

2 炮弹起始扰动分析

在炮弹飞行过程中，根据弹丸在空中飞行的实际弹道轨迹，将弹道划分为直线弹道段和曲线弹道段。相较于飞行全弹道，刚出炮口的一小段初始弹道，如图1中红色粗线段所示，弹丸的初速度较大，弹体结构偏心、质量偏心以及风等强迫扰动项对弹丸的作用还不显著，即弹道倾角速度≈0、弹道倾角加速度≈0，一般而言弹道倾角下降量小于 30′，此时的弹道可近似为直线[1-4]，通常把这段离炮口不远的一段弹道称为“直线段”[2]。

图1 弹道轨迹与直线弹道段示意图Fig.1 Schematic diagram of ballistic and the straight line segment

在直线弹道段，炮弹发射后效期结束瞬间，仅受起始扰动的影响，炮弹的动力学方程可用复攻角运动的2阶变系数齐次微分方程表示[17]：

起始扰动[4]，即当初始s=0时，由发射瞬间随机产生的复攻角∆0和复攻角速度引起的弹道攻角变化。采用“系数冻结法”求解式(1)，则复攻角∆的解可写为

由上述分析可知，由起始扰动引起的炮弹攻角变化是由两个角频率分别为ω1和ω2的快、慢圆运动复合成的，快圆运动是陀螺力矩作用的结果，慢圆运动是外力矩作用的结果。对于动态稳定良好的炮弹目标来说，必须满足λ1<0,λ2<0，则攻角的两圆运动的半径不断缩小，即模态振幅K1,K2的值不断变小。此时，每个圆运动都成为半径递减的螺线，随着攻角的不断缩小，两个圆运动逐渐合并，趋于单圆运动，再逐渐衰减而消失，复攻角的矢端将在复平面上画出不断缩小的外摆线或椭圆曲线。

3 雷达信号建模

3.1 坐标系解析

通过上述对炮弹在直线段仅受初始扰动影响的绕心运动特性分析，如图2所示，可将炮弹角运动分解为3个运动的叠加：即绕弹轴矢量OA的自旋、弹轴绕陀螺动量矩矢量G(即章动轴)旋转的章动(快圆运动)和陀螺动量矩矢量G绕速度矢量V1(即进动轴)旋转的进动(慢圆运动)。

图2 炮弹直线段角运动空间几何关系Fig.2 Spatial geometric relation of projectile angular motion in the straight line segment

进一步，结合弹道学理论[2]和测量雷达相关坐标系，建立炮弹目标直线段角运动空间几何关系，以雷达对炮弹目标的尾追式跟踪为例，设炮弹目标自旋角速度为，章动角为ϕ，章动角速度为ω1，进动角为ψ，进动角速度为ω2。

图2中引入了以下4种坐标系，包括全局坐标系(即雷达坐标系)Q-UV W、炮口坐标系O1-xyz、参考坐标系O-XY Z和弹轴坐标系O-ξηζ。炮口坐标系以炮口中心为原点O1，水平轴O1x为射击面与炮口水平面的交线，顺时针射向为正，铅直轴O1y在射击面内并与水平轴O1z相垂直；雷达坐标系以站址中心为原点Q，水平轴QU在水平面内指向正北方向，QV为铅垂线方向并与水平轴QW相垂直；参考坐标系以炮弹质心为原点O，与雷达坐标系始终平行；弹轴坐标系原点在炮弹质心上，Oξ轴沿弹轴向前为正，Oη轴垂直于弹轴向上为正，Oζ轴由右手法则确定。

设发射初始时刻，自旋轴OA、章动轴G和进动轴V1重合，参考坐标系原点与炮口坐标系原点重合，雷达坐标系Q-UV W、炮口坐标系O-xyz、参考坐标系O-XY Z和弹轴坐标系O-ξηζ之间的几何关系如图3所示。

图3 初始时刻4种坐标系空间几何关系Fig.3 Spatial geometric relation of the four coordinate systems as the projectile is stationary

此时，炮口坐标系Oz轴和弹轴坐标系Oζ轴重合，Oξη和Oxy共面，并绕Oζ轴逆时针旋转θ0角，θ0称为射角，即弹丸飞离炮口瞬间的弹轴矢量线与炮口水平面间的夹角(位于垂直面内)；炮口坐标系Oy轴和参考坐标系OY轴重合，Oxz和OXZ共面，并绕OY轴逆时针旋转α0角，α0称为射向角，即发射前炮口坐标系和参考坐标系水平面重合时，起始方向存在的夹角(位于水平面内)。则弹轴坐标系到参考坐标系的坐标变换矩阵可表示为

在发射瞬间，由于受到初始扰动的影响，弹轴出现章动和进动，章动角为ϕ，进动角为ψ，如图2所示。不考虑章动角和进动角初始相位，则弹轴坐标系向参考坐标系转换的坐标变换矩阵可表示为

可见，∠AOV1值 即是攻角值δ。式(5)用于下文验证回波建模的正确性。

3.2 回波信号建模

其中，根据Euler-Rodrigues绕矢量轴旋转公式，各旋转矩阵可展开为[5,18]

设雷达发射单频连续波信号，其表达式为

其中，f0为载频。则目标回波可表示为

将发射信号作为参考信号，与目标回波做相参处理，得到回波信号为

3.3 微多普勒分析

对式(10)中的相位项Φm(t)关于时间t求导，可得到回波的多普勒频率为

由式(12)可以看出，由弹轴在空间的角运动引起的微多普勒表达式较为复杂[19]。进一步分析可得如下结论[20]：即雷达微多普勒频率随时间的变化具有周期性，弹顶散射点运动周期Ttop为章动周期Tu和进动周期Tc的最小公倍数；尾翼散射点运动周期Tbom为自旋周期Ts、进动周期Tc和章动周期Tu的最小公倍数，即有如下关系式成立：

其中，k1,k2,k3∈N，N为自然数集合。

4 仿真分析

为验证第3节炮弹角运动对测量雷达回波信号的映射关系模型构建是否准确，分别选取典型高速自旋炮弹和低速自旋尾翼炮弹进行雷达微动信号建模和时频分析。采用单频连续波雷达测量炮弹目标，设载频f0=10 GHz。在炮口坐标系下，射向角α0=5°，射角θ=53°。考虑到实际中炮弹的回波将以顶点和底面圆环边缘点的回波为主，尾翼炮弹的回波将以顶点和尾翼散射点的回波为主，仿真中分别设置3个强散射点，分别为弹顶散射点1个、弹底边缘散射点或尾翼散射点2个。

4.1 高速旋转炮弹分析

如图2所示，在弹轴坐标系O-ξηζ下，设弹体有3个强散射点：弹顶散射点A和弹底边缘散射点B1,B2，仿真参数设置如表1所示。

表1 仿真参数设置(高速旋转炮弹)Tab.1 Simulation parameters setting (high-speed rotating projectile)

根据式(12)可推算得，弹顶散射点运动周期理论值Ttop为0.50 s，弹底边缘散射点运动周期理论值Tbom为0.50 s。考虑章动角和进动角衰减情况，仿真结果分别如图4-图6所示。

图4 快圆和慢圆不衰减情况(高速旋转炮弹)Fig.4 In the case of neither the fast nor the slow circular motion decays (high-speed rotating projectile)

图4(a)和图4(b)为运动补偿前的目标微多普勒理论值和Gabor变换时频图像，图4(d)和图4(e)为进行运动补偿后，得到的目标微多普勒理论值和Gabor变换时频图像。由图4(d)和图4(e)可以明显看出，运动补偿后弹底边缘散射点由弹体自旋、章动和进动的共同作用产生微多普勒频率调制，且包络幅度主要受章动和进动的影响，弹顶散射点的微多普勒仅由章动和进动产生。由于炮弹自旋频率高，时频图像出现自旋边带，且由于章动的调制，各散射点均出现锯齿波动。由图4(c)可见，当炮弹角运动的快圆和慢圆角度幅值均不衰减时，由弹顶散射点体现的攻角变化曲线呈现出明显周期性，且幅度较为平稳，体现了快圆和慢圆运动在空间上的耦合，符合炮弹运动攻角的变化规律。对图4(e)的时频图做自相关处理，其结果如图4(f)所示，可见，自相关处理后得到的时频图像中的微多普勒信号周期为0.50 s，与理论值一致；进一步应用骨架提取和Hough变换，也可以从时频图像中提取微多普勒信号的周期，经过5,10,20次蒙特卡罗试验获得的弹顶和弹底边缘微动周期均值如表2所示。

表2 Hough变换提取的时频图像周期(高速旋转炮弹)Tab.2 Extracted period of time-frequency image with Hough Transform (high-speed rotating projectile)

从表2的统计结果可以看出，经过20次蒙特卡罗试验得到时频图像中两个微动信号的周期都约为0.49 s，同时可以发现，随着试验次数的增加，微动参数提取结果的均值逐渐逼近真实值。此时，相对应的弹顶散射点运动周期估计值约为0.49 s，弹底散射点运动周期估计值约为0.49 s，与理论值及自相关法提取的微动周期十分接近，且对比可以看出自相关法提取的微动周期更加精确，后续将主要采用自相关法进行高速自旋炮弹的微动周期提取。

进一步，考虑快圆角度幅值快速衰减的情况，设衰减规律为ϕ′=ϕ·exp(-3.5t)，由图5(a)和图5(b)可见，运动补偿后，当章动幅值衰减时，弹底边缘散射点自旋微多普勒仅被进动所调制，此时弹顶散射点的微多普勒仅由进动产生，各散射点的锯齿波动均消失。同样，对图5(b)的时频图采用自相关法可得：弹顶散射点运动周期Ttop为0.50 s，弹底散射点运动周期Tbom为0.50 s，与理论值一致；由图5(c)可见，由弹顶散射点体现的攻角幅值曲线经章动角幅值衰减后，攻角幅值整体呈现衰减趋势，但均值整体接近进动角，符合炮弹运动攻角变化规律。

图5 快圆衰减而慢圆不衰减情况(高速旋转炮弹)Fig.5 In the case of the fast circular motion decays and the slow circular motion remains (high-speed rotating projectile)

当同时考虑快圆和慢圆角度幅值衰减的情况下，快圆衰减规律同上，设慢圆衰减规律为ψ′=ψ·exp(-0.05t)，由图6(a)和图6(b)可见，当进动角和章动角同时衰减时，由于进动角幅值相对较大，弹底边缘散射点和弹顶散射点的微多普勒主要由自旋和进动产生，各散射点锯齿波动均消失，进动衰减使弹顶和弹底散射点的微多普勒包络幅值在不断减小。对图6(b)的时频图采用自相关法可得，弹顶散射点运动周期Ttop为0.50 s，弹底散射点运动周期Tbom为0.50 s，与理论值一致；由图6(c)可见，由弹顶散射点体现的攻角幅值曲线经章动角幅值衰减后，同理攻角幅值整体呈现衰减趋势，但上下包络的中心表现为以进动角为初值的指数衰减形式，符合炮弹运动攻角变化规律。

图6 快圆和慢圆均在衰减情况(高速旋转炮弹)Fig.6 In the case of both the fast and the slow circular motion decays (high-speed rotating projectile)

4.2 旋转尾翼炮弹分析

如图2所示，在弹轴坐标系O-ξηζ下，设弹体有3个强散射点：弹顶散射点A和尾翼散射点B1,B2，仿真参数设置如表3所示。

表3 仿真参数设置(旋转尾翼炮弹)Tab.3 Simulation parameters setting (rotating tail projectile)

根据式(12)可计算得到弹顶散射点运动周期的理论值Ttop为2.00 s，尾翼散射点运动周期的理论值Tbom为2.00 s。考虑章动角和进动角衰减情况，仿真结果分别如图7-图9所示。

图7 快圆和慢圆不衰减情况(旋转尾翼炮弹)Fig.7 In the case of neither the fast nor the slow circular motion decays (rotating tail projectile)

图7(a)和图7(b)为未进行运动补偿的自旋尾翼炮弹微多普勒理论值和Gabor变换时频图像，图7(d)和图7(e)为进行运动补偿后得到的微多普勒理论值和Gabor变换时频图像。由图7(d)和图7(e)可见，尾翼散射点的微多普勒同时被自旋、章动和进动调制，曲线包络清晰可见，由于章动幅度偏小，其微多普勒受自旋和进动的影响较大；弹顶散射点的微多普勒主要由章动和进动产生。在快圆和慢圆角度幅值均不衰减的情况下，由图7(c)可见，由弹顶散射点体现的攻角变化曲线体现了快圆和慢圆角运动在空间上的合成，符合炮弹运动攻角周期性的变化规律。同样，对图7(e)的时频图做自相关处理，可以得到时频图像中的微多普勒信号周期为2.00 s，如图7(f)所示；进一步，应用骨架提取和Hough变换从时频图像中提取时频图像周期，同样经过多次蒙特卡罗试验进行微动周期的提取，表4分别给出5,10,20次蒙特卡罗试验获得的弹顶和尾翼的微动周期估计的均值。

表4 Hough变换提取的时频图像周期(旋转尾翼炮弹)Tab.4 Extracted period of time-frequency image with Hough Transform (rotating tail projectile)

从表4的统计结果可以看出，经过20次蒙特卡罗试验得到时频图像中两个微动信号的周期估计值分别约为2.00 s和1.99 s，相对应的弹顶散射点运动周期Ttop约 为2.00 s，弹底散射点运动周期Tbom约为1.99 s，与理论值及自相关法提取的微动周期十分接近，且可以看出自相关法提取的微动周期更加精确，后续将主要采用自相关法进行自旋尾翼炮弹的微动周期提取。

在考虑快圆角度幅值快速衰减的情况下，衰减规律为ϕ′=ϕ·exp(-1.5t)，由图8(a)和图8(b)可见，当章动角幅值衰减后，尾翼散射点的微多普勒主要受自旋和进动的影响，弹顶散射点的微多普勒主要受进动的影响。对图8(b)的时频图同样采用自相关法可得，弹顶散射点运动周期Ttop为2.00 s，尾翼散射点运动周期Tbom为2.00 s，与理论值一致；由图8(c)可见，由弹顶散射点体现的攻角幅值曲线经章动角幅值衰减过后，攻角幅值整体呈现衰减趋势，但均值整体接近进动角，符合炮弹运动攻角变化规律。

图8 快圆衰减而慢圆不衰减情况(旋转尾翼炮弹)Fig.8 In the case of the fast circular motion decays and the slow circular motion remains (rotating tail projectile)

在同时考虑快圆和慢圆角度幅值均衰减的情况下，快圆衰减规律同上，慢圆衰减规律为ψ′=ψ·exp(-0.05t)，由图9(a)和图9(b)可见，当进动角和章动角幅值同时衰减后，由于进动角幅值相对较大，尾翼散射点的微多普勒包络主要被自旋和进动调制，弹顶散射点微多普勒主要由进动产生，进动幅值的衰减使得弹顶和尾翼散射点的微多普勒包络幅值在不断减小。对图9(b)的时频图采用自相关法可得，弹顶散射点运动周期为2.00 s，尾翼散射点运动周期为2.00 s，与理论值一致；由图9(c)可见，由弹顶散射点体现的攻角幅值经章动角幅值衰减过后，攻角幅值整体呈现衰减趋势，但上下包络的中心表现为以进动角为初值的指数衰减形式，符合炮弹运动攻角变化规律。

图9 快圆和慢圆均在衰减情况(旋转尾翼炮弹)Fig.9 In the case of both the fast and the slow circular motion decays (rotating tail projectile)

5 实例验证

首先给出某型炮弹和某型末制导炮弹(带4个尾翼)尾追式跟踪得到的直线段回波时频分析仿真结果，考虑到实际炮弹目标的瞬时速度变化，图10和图11分别给出瞬时径向速度随时间的变化，仿真采用连续波雷达，载频与实际装备完全相同。

图10 某型炮弹直线段回波时频分析仿真结果Fig.10 Simulated time-frequency image of projectile in the straight line segment

图11 某型末制导炮弹时频分析仿真结果Fig.11 Time-frequency image of terminal guided projectile in the straight line segment

在上述理论和仿真分析的基础上，进一步结合实测数据进行验证。在试验鉴定中，采用某型连续波雷达尾追式跟踪测量，分别得到炮弹和末制导炮弹在初始直线弹道段的两幅时频图像，分别如图12和图13所示。

图13 某型末制导炮弹直线段回波真实时频图像Fig.13 Time-frequency image of terminal guided projectile in the straight line segment

从图12可以看出，炮弹目标真实回波的时频分析结果整体呈现衰减形式，这是由于目标远离雷达的速度引起的，且对比图10可以看出，实际目标回波受到遮挡效应的影响，时频图像主要表现为弹底散射点的多普勒，与图10中间亮线的变化趋势一致。弹底边缘散射点的微多普勒由自旋、章动和进动共同作用产生，其微多普勒包络主要受章动和进动的影响。由于炮弹的结构尺寸小、自旋频率高，进动的调制形成比较明显的包络，而章动的调制对自旋微多普勒边带形成更为微小的锯齿波动。随着攻角逐渐变小，进动调制包络也逐渐衰减。但整体而言，由于炮弹底面半径较小且弹底散射点微动幅度小，其微多普勒受章动和进动的调制效果不明显。

图13给出末制导炮弹直线段真实回波的时频图像，同样与图11的仿真结果变化趋势相同。但实际应用场景中，由于末制导炮弹自旋频率低，尾翼散射点微多普勒受自旋、章动和进动的复合调制，同时进动的调制形成明显的包络，而章动幅值相对较小，章动调制效果不明显。随着攻角变小，调制包络也逐渐衰减。由于带尾翼的炮弹弹底散射点微动幅度较大，受进动调制效果明显可见。

通过仿真分析与实际回波时频分析结果的进一步对比可以发现，两者变化的整体趋势是一致的，证明了本文的建模方法的正确性和可行性。

6 结论

本文通过分析炮弹目标在直线段运动过程中弹轴和速度矢量之间角运动的变化机理，构建了炮弹角运动对雷达测量信息的投影变换模型，通过选取典型炮弹目标进行仿真分析和实测数据验证得出如下结论：(1)炮弹目标受空气阻力的影响，在随弹道前进过程中，微动受到自旋、进动和章动的复合调制，且进动角和章动角的初值存在差异，随大气干扰而逐步衰减，与弹道导弹目标的微动建模存在明显不同，两者在弹道术语(如进动、章动等含义)方面的界定存在明显差异；(2)由于炮弹自旋频率较高，章动衰减快、幅度小，进动产生的微多普勒效果比较明显，体现在雷达时频图像中进动特征比较明显；末制导炮弹自旋频率较低，章动衰减快、幅度小，由自旋和进动产生的微多普勒效果比较明显，体现在雷达时频图像中自旋和进动特征明显。实测数据验证了所建立的基于炮弹目标角运动回波模型是正确有效的，符合炮弹运动规律。

以上研究结论为炮弹类目标的微动回波建模探索了新的途径，为寻求弹道目标与炮弹目标微动特征差异进行了崭新的尝试，为后续利用实际雷达测量信息辨识炮弹目标运动特性、气动特性、结构特性和参数匹配辨识提供了参考依据。

在实际工程应用中，由于受到雷达信号形式、雷达照射角度、目标散射特性、结构特性、信噪比和自然环境、气象条件等多种因素的共同影响，基于各型炮弹目标的角运动雷达回波信号建模更为复杂，下一步需要持续开展相关研究，为炮弹这一类特殊目标的微动特征提取与识别提供参考。