高轨预警卫星预警信息生成方法研究*

周童，张雅声，周海俊，汉京滨

(航天工程大学，北京 101400)

0 引言

弹道导弹的早期预警设备包括高轨预警卫星和远程预警雷达，高轨预警卫星由于位置较高将率先发现处于主动段的弹道导弹，并利用导弹尾焰的红外辐射对导弹进行跟踪。远程预警雷达可以根据高轨预警卫星的指引在导弹关机后对导弹进行接力跟踪[1-2]。对远程弹道导弹来说，高轨预警卫星探测到导弹关机后，远程预警雷达由于探测距离限制并不能及时捕获到导弹目标，在这种情况下就需要对高轨预警卫星的探测数据进行处理，生成预警信息，对远程预警雷达的搜索与捕获过程进行指引[3-5]。

高轨预警卫星与远程预警雷达的交接问题主要涉及2个过程：一是预警卫星探测到导弹关机后的弹道及误差预报问题；二是2种设备交接班利用的预警信息的表示方法。在弹道预报方面，当前的弹道预报方法主要有解析几何法和数值积分法[6]。解析几何法的依据是椭圆弹道理论，将导弹自由段的飞行轨迹近似为椭圆弹道，利用解析几何法实现弹道的预报[7-9]。数值积分法的依据是导弹的动力学方程，利用数值逼近的方法计算导弹弹道[10-11]。解析几何法计算过程简单，速度较快，但计算精度较低；数值积分法计算过程复杂，速度较慢，但计算精度较高。考虑到计算精度要求，数值积分法在实际工程中应用更为广泛。在预警信息表示方面，相关研究还较少。文献[12]提出了一种导弹预警卫星引导能力评价方法，在预警信息含义方面没有详细表述。文献[13]给出了粗略的预警信息生成模型，但缺乏理论依据且空域的划分比较粗略。

针对以上问题，本文首先给出了基于数值积分法的弹道及误差预报过程，然后将弹道预报的结果与概率统计理论中的高斯概率误差分布相结合提出了4种预警信息的表示方法，最后通过仿真验证了预警信息生成方法的正确性。

1 基于数值积分法的弹道预报

对于远程弹道导弹来说，发动机关机后导弹将在大气层外无动力飞行，此时可以利用考虑多种摄动因素的二体运动模型对导弹运动进行描述，并以此建立导弹的运动学方程。由于导弹的动力学方程为变系数非线性常微分方程，因此可以采用数值积分的方法对方程进行求解。数值积分法计算弹道的一般过程为首先建立动力学方程，而后选取适当的积分方法按照一定的积分步长进行计算，在给定初始状态的情况下，可以得到任意时刻导弹的运动状态。

1.1 建立目标运动学方程

(1)

式(1)为时变的非线性微分方程，无法获得解析解，通常采用数值积分的方法获得其数值解。数值积分的过程可以描述为：已知目标在k-1时刻的状态xk-1，那么可以通过积分获得目标在下一时刻k的运动状态为

(2)

利用式(2)，对导弹动力学方程采用合适的数值积分方法，就可以计算得到导弹弹道[14]。

1.2 弹道误差传播分析

在利用数值积分法计算弹道的过程中，通常利用协方差分析描述函数法(CADET)对弹道误差的传播进行分析，CADET的一般过程为首先对系统进行拟线性化处理，然后利用线性化系统协方差分析原理计算随机状态变量的协方差。

(3)

式中：m为位置矢量r的均值。

因此拟线性系统的动态矩阵描述函数N可以表示为

(4)

式中：I为3×3的单位矩阵；

利用CADET方法可得导弹的状态均值和协方差的传播方程：

(5)

2 基于统计理论的状态预测误差概率分布

弹道导弹在发动机关机后进入中段飞行，由于中段飞行环境的复杂性、弹道目标和卫星运动的非线性以及预报算法模型的建模精度等因素的影响，对目标运动状态预测误差分布进行精确建模分析非常困难。通常的做法是将这些误差源带来的随机误差统一处理并建模为服从零均值的高斯分布，并将传感器视线观测噪声建模为零均值高斯噪声。根据概率统计理论中心极限定理，可以将目标运动状态预测误差建模为高斯分布。

基于统计理论的高斯分布概率模型给出了状态估计量在高斯假设条件下的概率分布计算[15]，下面对其进行分析。

设ξ为状态量m的估计，ξ为随机变量且满足高斯分布，则ξ的概率分布函数为

(6)

式中:N为状态矢量ξ的维数；m为状态矢量真实值；P为状态估计量ξ的状态误差协方差矩阵，即

P=E[(ξ-m)(ξ-m)T].

(7)

若P为对称正定矩阵，则存在逆阵P-1，且P-1为对称正定矩阵。令

κ=(ξ-m)P-1(ξ-m)T.

(8)

(9)

(10)

可知，式(9)和式(10)不但描述了误差概率密度随误差绝对值的变化，而且描述了误差概率密度随误差方向的变化规律。

对位置误差分布概率展开研究,令

(11)

(12)

根据矩阵分析理论，存在且唯一存在正交矩阵A，有AT=A(AT为A的转置)，且使得

(13)

式中:λ1,λ2,λ3为协方差矩阵Pr的特征值。

令ζ=AT,γ，并将式(13)代入式(12)中，整理得

(14)

(15)

(16)

式中:Γ()为gamma函数;erf()为误差函数,表达形式为

(17)

3 预警信息生成模型

不同的交接班任务对预警信息的要求是不一样的。较宽视场的接班设备覆盖空域大，只需预警信息提供目标的概略空域分布；接班设备视场再小一些，就需要缩小接班传感器搜索的空域，此时需要预警信息提供目标较为精确的空域分布；对于采用窄视场搜索的接班设备来说，为保证搜索数据率，快速发现目标，则需要预警信息提供尽可能精确的目标空域分布，同时还能指示目标在空域中出现的概率。本节将按照上述思路，由精细到粗略分别研究建立概率椭球模型、概率球模型、均匀椭球模型和均匀球模型4种粒度的预警信息模型。

3.1 概率椭球模型

概率椭球模型是以目标预测状态的分布概率密度为基础，用误差椭球精细描述目标状态的分布概率。根据式(9)和式(10)可知，第k时刻目标位置与速度状态预测值的概率分布密度函数分别为

(18)

式中:Pr,k和Pv,k都是3×3的非负正定Hermitian矩阵，因此存在3×3的正定矩阵A和B，有AAT=I，BBT=I，且使得

(19)

(20)

式中：λr,i，λv,i(i=1,2,3)分别为协方差矩阵Pr,k和Pv,k的特征值，且有如下关系：

(21)

(22)

由式(22)可知，对于给定的概率Pr和Pv，相应的κ的取值与状态变量无关，可以直接求解κr和κv的值。下面不加区分的给出概率P取不同值时，利用Matlab符号工具箱求解得到的对应的κ值，即

(23)

概率椭球模型是以预测值为原点、等概率面为椭球面的三维联合高斯分布模型。它不仅描述了目标状态的分布概率密度随偏离预测值距离的变化规律和变化方向，而且描述了目标状态位于椭球体内的实际分布概率大小。

3.2 概率球模型

概率球模型以目标预测状态的分布概率密度为基础，用误差椭球最大误差半径所在轴方向的概率分布密度来描述状态的概率分布密度，用误差球来描述目标状态的分布概率。概率球模型的概率分布密度仍是高斯型概率分布密度函数的形式，但是，概率分布的大小以球体来描述。

假设已由式(19)和式(20)得到第k时刻目标位置与速度状态估计误差协方差矩阵的特征值，并且特征值满足如式(21)所示，则第k时刻目标位置与速度状态预测值的概率分布密度分别表示为

(24)

概率球在概率椭球的基础上，将原来概率椭球各轴向不同的概率分布密度统一以同一个概率分布密度来描述目标预测状态的概率分布密度。其概率积分本质上还是对高斯概率分布密度积分，利用前文推导的结果，直接给出以下结论。

概率球模型是以预测值为原点、等概率面为球面的三维联合高斯分布模型。它描述了目标状态的分布概率密度随偏离预测值距离的变化规律，但没有描述概率密度变化方向，即认为在概率球内目标状态分布概率密度是各向同性分布。概率球模型也描述了目标状态位于概率球内的概率分布大小，但是，由于其弱化了目标状态的实际概率分布密度，所以目标在概率球内的实际概率分布大于理论得到的概率分布。

3.3 均匀椭球模型

均匀椭球模型是在概率椭球模型的基础上，认为在椭球内部目标状态的概率分布密度是均匀的。因此，采用均匀椭球模型的预警信息首先需要基于概率椭球模型，然后根据给定的概率分布分别计算椭球体三轴半径。

假设已由式(19)和式(20)得到第k时刻目标位置与速度状态估计误差协方差矩阵的特征值，并且特征值满足如式(21)所示，则第k时刻目标位置与速度状态预测值的概率分布密度分别表示为

(25)

均匀椭球模型弱化了目标状态的概率密度分布特性，将内部每一点都视为均匀分布，但其描述了目标状态的概率分布的变化方向。因此，采用均匀椭球模型生成的预警信息能够用椭球三轴半径的大小来描述目标状态概率分布的变化方向，但是不能描述目标状态概率分布密度的变化趋势。

3.4 均匀球模型

均匀球模型是在概率球的基础上，认为误差球内部目标状态的概率分布密度是均匀的。

假设已由式(19)和式(20)得到第k时刻目标位置与速度状态估计误差协方差矩阵的特征值，则第k时刻目标位置与速度状态预测值的概率分布密度分别表示为

(26)

均匀球模型弱化了目标状态的概率密度分布特性，认为目标状态概率分布各向同性，且内部每一点都为均匀分布。因此，采用均匀球模型生成的预警信息粒度是最粗的，在4种模型中能够提供的指示信息最少，精确最差，能力最弱。但是，该模型形式简单且易于描述，在接班设备视场覆盖能力较强、覆盖空域较大的情况下，能够达到较好的引导效率。

4 仿真校验

本节将通过仿真，对上文提出的高轨预警卫星预警信息生成过程进行模拟分析，并利用Monte Carlo仿真验证预警信息生成方法的正确性。

4.1 弹道误差传播与预警信息生成

设预警卫星估计得到导弹在关机点时刻t0的状态均值m0和协方差P0分别为m0=(2 155 145 m,6 002 997 m，-1 069 m，-426 m/s，5 867 m/s,1 843 m/s)，

利用数值积分法对导弹状态均值和协方差进行预报计算，得到t=600 s时的导弹状态均值m和协方差P。数值积分采用四阶龙格-库塔方法，积分步长设为0.1 s。图1给出了导弹状态变量均值随时间变化的仿真结果。

利用计算得到的t=600 s时的导弹状态均值和协方差可以根据本文提出的预警信息生成模型得到导弹的预警信息。由前文分析可知，基于概率椭球模型生成的预警信息最具代表性，因此本文以概率椭球模型预警信息生成为例开展仿真。图2为根据t=600s时的导弹状态生成的地心惯性坐标系下导弹位置和速度的预警信息。

图中椭球体所包含的区域为基于本文预警信息生成方法计算得到的理论误差分布，外圈椭球体以99.5%的概率涵盖导弹位置和速度分布区域，内圈椭球体以95%的概率涵盖导弹位置和速度分布区域。值得注意的是误差椭球的三轴方向并不一定平行于地心惯性系的3个坐标轴，椭球的三轴方向为误差变化量最大的3个方向，可由协方差矩阵的特征向量确定。

4.2 预警信息生成方法正确性验证

本节采用的验证方法如下：根据t0时刻的导弹关机点状态估计值及误差协方差矩阵作为状态初值，根据数值积分法预报t=600 s时的导弹运动状态，并基于本文提出的预警信息生成方法计算t=600 s时的导弹状态的理论误差分布。同时，取t0时刻导弹真实弹道数据作为真值，以t0时刻导弹运动状态误差为标准差，预报导弹在t=600 s时的运动状态，以500次Monte Carlo仿真结果作为导弹状态预测的实际分布。最后，对比考察理论误差分布和实际分布在统计意义上是否符合一致性。

图3为导弹位置预警信息生成方法正确性验证结果，三维图中*为500次Monte Carlo生成的导弹实际位置，外圈椭球体以99.5%的概率涵盖导弹位置分布区域，内圈椭球体以95%的概率涵盖导弹位置分布区域。为直观比较，将导弹理论状态与实际状态分别投影到地心惯性坐标系的3个主平面，椭圆中心为导弹状态均值。其中‘—’代表误差分布概率为99.5%，‘--’代表误差分布概率为95%。

图4为导弹速度预警信息生成方法正确性验证结果，图中坐标轴与各标识含义与图3中相同。

5 结束语

本文以导弹防御为背景，研究了导弹在发动机关机后自由飞行过程中运动误差的传播特性及误差的描述方法。文章首先给出了基于数值积分法的弹道误差传播方法，而后基于数值积分法获得的导弹状态由细到粗提出了4种导弹预警信息模型，并通过仿真验证了所提出模型的正确性。文章的研究内容可以推广到我国导弹预警系统的建设中，但在实际运用中预警信息并非越精细越好，预警信息越详细，交班设备生成预警信息的时间就越多，整个交接班过程效率反而可能达不到最优，根据具体任务选择合适的预警信息生成方法将会大大提高导弹预警系统的防御能力。