临近空间升力式再入飞行器跳跃空域预测

喻晨龙， 谭贤四， 曲智国， 谢 非， 李陆军

(1. 空军预警学院研究生管理大队，湖北 武汉 430019； 2. 空军预警学院防空预警装备系， 湖北 武汉430019；3. 长江勘测规划设计研究院， 湖北 武汉 430019)

临近空间升力式再入飞行器的飞行高度高、飞行速度快、突防能力强、打击距离远，具有升力式气动布局，其由火箭助推至高空后，依靠地球重力再入大气层，可长距离、无动力滑翔飞行，典型的有HTV-2和AHW两种[1-3]。这类再入式飞行器在气动外形、运行环境、飞行受力、控制方式和机动模式等方面与弹道导弹存在明显差异，给预警探测系统带来了极大的挑战[4-8]。针对这一杀手锏武器，目前可行的探测手段是以反导预警系统为依托，以米波雷达、凝视雷达等新型装备为补充进行协同探测，从而尽早尽快发现目标并维持较高的跟踪精度，为后端的拦截系统提供连续稳定的信息支撑。

通常，远程预警相控阵雷达靠近边境线前沿部署，探测距离可达5 000 km，然而受地球曲率的限制，在探测临近空间目标时探测距离最远只有1 300 km。临近空间升力式再入飞行器在纵向上存在跳跃机动，在横向上存在侧向偏移，其中：只有当目标处于雷达探测范围内时，才能保持连续稳定跟踪；当目标飞出探测范围，即向下机动至地平线以下时会使情报中断。若能利用目标在雷达水平面波束层中飞行时获得的有限观测信息，对目标的运动包络和可能的跳跃空域进行预测，则可辅助雷达在远程预警相控阵雷达丢失目标后补充照射，从而保证情报的连续。

升力式再入飞行器的弹道包括平衡滑翔和跳跃滑翔[9]。研究者们针对目标的运动参数、攻角和倾侧角的变化规律进行了深入分析，如：李广华等[10]认为，在平衡滑翔状态下目标的速度与高度、速度与速度倾角之间存在解析关系，且随着倾侧角的变化而变化；文献[11-12]的作者通过分析气动力加速度对目标质心运动和绕质心运动的影响，说明气动力改变是目标跳跃机动的主要原因；李广华[13]认为目标的运动状态参数反映着目标的气动力，只有当参数满足特定数学关系时，才能维持平衡滑翔轨迹，当偏离平衡点时就会产生跳跃轨迹；王路等[14]通过目标的三自由度弹道方程计算了目标的攻角和最大升阻比。然而，上述研究未利用运动参数、攻角和倾侧角的变化规律进行弹道预测。鉴于此，笔者在分析目标的运动特性的基础上，估算目标的运动参数、攻角和倾侧角，合理假设攻角和倾侧角的变化规律，最后通过数值仿真方法预测目标飞离雷达探测范围后可能的跳跃空域。

1 模型建立

1.1 飞行器的三自由度运动模型

为了描述目标的运动过程和雷达的探测过程，定义地心(Earth Canter，EC)坐标系、飞行器东北天(Vehicle East-North-Up，VENU)坐标系和雷达站东北天(Radar East-North-Up, RENU)坐标系3个相关的参考坐标系[14]。这3个坐标系各自包括1个直角坐标系和1个球坐标系。

假设地球为均匀圆球，不考虑自转影响，则飞行器在无动力滑翔过程中满足EC球坐标系和VENU直角坐标系的三自由度运动方程[8，14],即

(1)

(2)

式中：r(t)、χ(t)、φ(t)分别为飞行器在EC球坐标系下的地心距、经度和纬度；v(t)、γ(t)、ψ(t)分别为飞行器在VENU直角坐标系下的速度、航迹倾角和航向角；m为飞行器的质量；g为重力加速度；σ(t)为倾侧角；D(t)、L(t)分别为飞行过程中所受阻力和升力,

(3)

(4)

其中CD(t)为阻力系数，CL(t)为升力系数，ρ(t)为大气密度，S为飞行器气动参考面积。

文献[15-18]作者认为：当飞行器的速度>8 Ma时，升力系数可简化为攻角α(t)的线性函数，阻力系数可简化为α(t)的二次函数。对气动参数插值拟合后可得

CD(t)=0.044+0.08α(t)+2.2α2(t),

(5)

CL(t)=-0.151 7+3.151α(t)。

(6)

大气密度满足指数模型，即

ρ(t)=ρ0exp[-(r(t)-Re)/hs]，

(7)

式中：ρ0=1.752 kg/m3，为参考大气密度；hs=6.7 km，为参考高度；Re=6 371 km，为地球半径。

1.2 雷达的观测模型

飞行器穿越水平波束屏的过程中，雷达每间隔一段时间Ts(即搜索间隔时间)照射一次目标，得到了N组观测值{z1,z2,…,zN}，这些观测值是在RENU球坐标系下获取的，记为

zi=[ϑr(ti)θr(ti)φr(ti)]T，

(8)

式中：ϑr(ti)、θr(ti)、φr(ti) 分别为第i次照射时的目标径向距离、方位角、俯仰角，其中ti为第i次照射时对应的时刻。将观测值在RENU直角坐标系中表示，有

(9)

通过式(9)位置差分可得到速度信息Vr(ti)=[vrx(ti)vry(ti)vrz(ti)]T。

2 参数估计

防御方预警系统可直接获取的信息仅包括目标在穿越水平波束屏过程中的观测信息，即RENU直角坐标系下的位置和速度。而目标的面质比、控制量，以及在EC球坐标系下的地理位置和VENU直角坐标系下的速度是未知的。

由式(1)-(7)可知：目标的弹道轨迹取决于运动状态初值和攻角/倾侧角的变化规律。为了实现轨迹预测，首先把观测信息由RENU直角坐标系转换到EC球坐标系和VENU直角坐标系，得到目标飞离水平波束屏的运动状态，然后结合式(1)-(7)估算目标的攻角/倾侧角和面质比。

2.1 运动参数初值估计

1) 将位置观测信息从RENU直角坐标系转换到EC球坐标系。首先从RENU直角坐标系转换到EC直角坐标系，然后从EC直角坐标系转换到EC球坐标系。

假设雷达站OR所在地的经度为J，纬度为W，海拔高度为H。要从RENU直角坐标系变换到EC直角坐标系，根据欧拉角关系，首先把RENU直角坐标系原点先沿地心矢移动Re+H到地心O，然后绕x轴顺时针旋转π/2-W，最后绕z轴顺时针旋转π/2+J。即

(10)

式中：Xr(ti)=[xr(ti)yr(ti)zr(ti)]T，为RNEU直角坐标系下的位置；Xe(ti)=[xe(ti)ye(ti)ze(ti)]T，为EC直角坐标系下的位置；

(11)

为从RENU直角坐标系到EC直角坐标系的坐标变换矩阵，其中Tx、Tz分别为绕x、z轴的基元旋转矩阵。

将位置观测信息从EC直角坐标系转换到EC球坐标系，即

(12)

假设在第k次照射后飞行器飞离水平波束屏，雷达测量得到目标的质心位置为[ϑr(tk)θr(tk)φr(tk)]T(其中tk为第k次照射时对应的时刻)，这是一个RENU球坐标，通过坐标变换得到此时质心在EC球坐标系下位置为[re(tk)χe(tk)φe(tk)]T。

2) 将速度信息从RENU直角坐标系转换到VENU直角坐标系。先考虑位置坐标的变换，将位置坐标从RENU直角坐标系转换到EC直角坐标系，再从EC直角坐标系转换到VENU直角坐标系。由于从EC直角坐标系转换到VENU直角坐标系，与从RENU直角坐标系转换到EC直角坐标系的变换关系刚好相反，因此可类比式(10)中的方法转换。特别需注意的是，此时坐标系的原点变为了飞行器的质心。

要从EC直角坐标系转换到VENU直角坐标系，根据欧拉角关系，首先把EC直角坐标系先绕z轴逆时针旋转π/2+χe(tk)，然后绕x轴逆时针旋转π/2-φe(tk)，最后坐标原点由地心O沿地心矢移动re(tk)到飞行器的质心。即

(13)

式中：Xe(ti)=[xe(ti)ye(ti)ze(ti)]T和Xv(ti)=[xv(ti)yv(ti)zv(ti)]T分别表示第i次照射时目标在EC直角坐标系和VENU直角坐标系下的位置；

(14)

为从EC直角坐标系到VENU直角坐标系的坐标变换矩阵。

那么，从RENU直角坐标系到VENU直角坐标系的位置坐标变换关系为

(15)

因此，从RENU直角坐标系下的速度矢量[vrx(ti)vry(ti)vrz(ti)]T转换到VENU直角坐标系下的速度矢量[vvx(ti)vvy(ti)vvz(ti)]T，其变换关系为

(16)

由于VNEU直角坐标系下的速度、航迹倾角与航向角有如下关系:

(17)

因此，根据RENU直角坐标系下的一组量测信息差分得到飞行器飞离水平波束屏的速度[vrx(tk)vry(tk)vrz(tk)]T，然后根据坐标转换关系把各轴上的分速度变换到VENU直角坐标系下，最后得到飞行器的速度vv(tk)、航迹倾角γv(tk)和航向角ψv(tk)。

2.2 攻角和倾侧角估计

式(1)、(2)描述了目标的地理位置和速度的迭代关系，式(1)-(7)隐含了目标控制量的解算方法。因此，将式(2)改写为[20]

(18)

(19)

(20)

式中：ρ(ti)=ρ0exp[-(re(ti)-Re)/hs]。

(21)

(22)

(23)

式中：K(ti)为目标的升阻比。

将式(21)代入式(5)、(6)，可以解算出目标第i次照射时的攻角α(ti)、升力系数CL(ti)和阻力系数CD(ti)。

3 跳跃空域预判

当初值确定后，飞行器的弹道由攻角和倾侧角决定，攻角和倾侧角并不是随机变化的。其中：1) 攻角直接决定飞行过程中的气动特性，再入过程中若调整频繁且幅值较大，将严重影响飞行器的热防护系统，因此攻角的变化要尽量简单，一般选择分段线性函数，大小通常控制为0°～ 20°，此外攻角与升阻比紧密相关，最常用的制导律是以最大升阻比飞行;2) 倾侧角是主控制量，其确定升力在纵向和侧向上的分量，并决定飞行器横向机动程度，改变倾侧角是为了规避禁飞区和突防，倾侧角最常用的制导律是方位角误差走廊，其幅值不变、符号翻转，幅值通常控制在±40°以内。

因此，升力式再入飞行器在远端运动时，为了维持最小能耗和最大升阻比飞行，并且实现有效突防，对攻角和倾侧角的变化规律可做如下假设：

(24)

(25)

式中：α(tk)为飞行器飞出水平玻束屏时估计的攻角；αmaxK(tk)为最大升阻比对应的攻角；σ(tk)为飞行器飞出水平玻束屏时估计的倾侧角；t为飞行时间；T为攻角改变的时间；T1为倾侧角改变的时间。

将飞行器飞离水平波束屏时的运动参数、攻角和倾侧角的预测值代入式(1)-(7)进行数值仿真，随着迭代的进行，当式(2)中第1次出现γ(t)=0时停止。

4 仿真验证

假设某型远程预警相控阵雷达部署地点在东经110°、北纬10°、海拔高度0 km，波束宽度1.2°，沿着正东设置单层水平波束屏，量测噪声服从高斯分布，测量误差ϖϑ=100 m，ϖθ=ϖφ=0.001 5°，搜索间隔时间Ts=1 s。飞行器穿越水平波束屏的过程中雷达在RENU球坐标系内得到了14组观测数据，如表1所示。

表1 RENU球坐标系内观测数据

按照2.1节的估算方法对观测数据进行坐标转换，计算得到目标飞离水平波束屏时的地理位置和速度分别为：re(tk)=6 425 km，χe(tk)=117.472 8°，φe(tk)=9.889 5°；vv(tk)=19.850 9 Ma，γv(tk)=-1.868 6°，ψv(tk)=92.029 4°。进一步，根据2.2节的估算方法，计算得到目标的面质比为S/m=0.002 5，在飞离水平波束屏时α(tk)≈5°，σ(tk)≈40°。

飞行器飞离水平波束屏后，其攻角和倾侧角按照式(24)、(25)的规律变化，即攻角和倾侧角先是维持不变，再分别在时刻T和T1发生改变，攻角由飞离时的攻角α(tk) 变化为最大升阻比攻角αmaxK(tk) ，倾侧角幅值不变、符号翻转。这时，随着迭代的进行，轨迹发生较大变化。由于攻角和倾侧角改变的时刻是无法预知的，当式(2)中第1次出现γ(t)=0时迭代停止，此时的轨迹曲线就围成了一个空域范围Ω，在EC球坐标系下的目标跳跃空域范围如图1所示。

5 结论

临近空间武器的兴起将打破常规作战规则，在探测这类目标时，远程预警雷达因受到地球地平线遮挡而丢失目标，本文基于目标弹道特性，利用有限的观测数据对目标可能的跳跃空域进行预判，能提前指示预警系统调度雷达资源探测，具有较强的军事应用价值。

参考文献：

[1] PAUL L M,VINCENT L R,LUAT T N,et al.NASA hypersonic fight demonstration overview,status,and future plans[J].Acta astronautica,2004,55(7):619-630.

[2] 李益翔.美国高超声速飞行器发展历程研究[D].哈尔滨：哈尔滨工业大学,2016：18-32.

[3] STEVEN H W,JEFFREY S,DALE S,et al.The DAPPA/AF falcon program:the hypersonic technology vehicle #2(HTV-2) flight

demonstration phase[C]∥15th AIAA International Space Planes and Hypersonic Systems and Technologies Conference.28 April-1 2008,Dayton,Ohio:2008-2539.

[4] RANDALL T V,LAERENCE D H,CHARLES R,M,et al.X-43A Hypersonic vehicle technology development[J].Acta Astronautica,2006(59):181-191.

[5] STEVEN W,MING T,SUE M,et al.Falcon HTV-3X-A reusable hypersonic test bed[C]∥15th AIAA International Space Planes and Hypersonic Systems and Technologies Conference.28 April-1 May,2008,Dayton,Ohio：2008-2544.

[6] 秦雷，李君龙，周荻.国外临近空间高速滑翔弹头运动特征研究[J].战术导弹技术,2015(1):16-20.

[7] 李君龙，秦雷，谢晓瑛.临近空间非弹道式机动模式与跟踪滤波问题[J].现代防御技术,2017,45(1):99-106.

[8] 李慧峰.高超声速飞行器制导与控制技术[M].北京：中国宇航出版社,2012:10.

[9] 陈小庆，候中喜，刘建霞.高超声速滑翔飞行器弹道特性分析[J].导弹与航天运载技术,2011(2):5-9.

[10] 李广华，张洪波，汤建国.高超声速滑翔飞行器典型弹道特性分析[J].宇航学报，2015,36(4):397-403.

[11] OLIVIER D M,ROBERT H B.Tracking and identification of a maneuvering reentry vehicle[C]∥AIAA Guidance,Navigation,and Control Conference and Exhibit.11-14 August 2003,Austin,Texas,USA,2003:2003-5446.

[12] 张凯，熊家军，韩春耀，等.一种基于气动力模型的高超声速滑翔目标跟踪算法[J].宇航学报,2017,38(2):123-130.

[13] 李广华.高超声速滑翔飞行器运动特性分析及弹道跟踪预报方法研究[D].长沙:国防科学技术大学,2016：17-49.

[14] 王路，邢清华，毛艺帆.基于升阻比变化规律的再入高超声速滑翔飞行器轨迹预测算法[J].系统工程与电子技术,2015,37(10):2335-2339.

[15] PHILIPS T H.A Common Aero Vehicle(CAV) model,description,and employment guide[R].USA:Schafer Corporation,2003.

[16] 赵汉元.飞行器再入动力学与制导[M].长沙:国防科技大学出版社,1997：48-97.

[17] 黄长强，国海峰，丁达理.高超声速滑翔飞行器轨迹优化与制导综述[J].宇航学报,2014,35(4):369-379.

[18] 王璐璐，秦玉亮，王宏强，等.飞行器升阻比估计误差分析方法研究[J].雷达科学与技术,2012,10(2):175-179.