中国科技核心期刊      中国指挥与控制学会会刊     军事装备类重点期刊
Command Control and Simulation

Command Control and Simulation >

2019 , Vol. 41 >Issue 5: 47 - 53

DOI: https://doi.org/10.3969/j.issn.1673-3819.2019.05.011

Theory & Research

Prediction and Error Analysis for Radar Handover Time Based on Encke Method

  • WANG Biao 1, 2 ,
  • WU Nan 1 ,
  • MENG Fan-kun 1 ,
  • WANG Jian-tao 1
Expand
  • 1. Information Engineering University, Zhengzhou 450000
  • 2. Unit 95129 of PLA, Kaifeng 475000, China

Received date: 2019-04-23

  Revised date: 2019-05-15

  Online published: 2022-05-09

Copyright

Copyright reserved © 2019
Fold

Abstract

Trajectory prediction is an important step in radar handover time prediction in missile early warning. In order to improve the prediction accuracy of radar handover time, Enke method is applied to the ballistic missile prediction. Enke method establishes the reference track of the missile movement by using the two-body mechanical model, solves the differential equation of the deviation between the reference track and the actual track to obtain the deviation value, combines the reference track and the deviation value to get the actual orbit, then gets the radar handover predicted time by this trajectory prediction. In the error propagation analysis, firstly the state covariance matrix at the handover time is calculated by using the method of unscented transform, then the error ellipsoid at the handover time is calculated by the state covariance matrix. The simulation results show that, compared with the traditional trajectory prediction method, Enke method considers the influence of the perturbation factor in missile motion, can integrate over a long time, and take less time to obtain higher accuracy. It is of high application value for radar handover time prediction.

Key words： trajectory prediction; radar handover time; Encke method; error analysis; unscented transform

Cite this article

WANG Biao , WU Nan , MENG Fan-kun , WANG Jian-tao . Prediction and Error Analysis for Radar Handover Time Based on Encke Method[J]. Command Control and Simulation, 2019 , 41(5) : 47 -53 . DOI: 10.3969/j.issn.1673-3819.2019.05.011

弹道导弹射程远、威力大,已经成为大国战略威慑的重要武器。为应对这种威胁,近年来,用于导弹防御的导弹预警系统得到了很大发展。由于弹道导弹在自由段的飞行时间和距离占全弹道的80%以上[1],因此,自由段成为导弹预警研究的主要过程。导弹在自由段速度快、距离远,而单部雷达探测范围有限,往往需要多部雷达进行预警跟踪。在多雷达预警中,当前一部雷达与后一部雷达对导弹的探测范围不重叠时,需要进行雷达交接时刻的预报[2,3],来提高整体预警能力,使用有效的弹道预报方法可以更好地提高雷达交接时刻的预报性能[4,5]。
传统的弹道预报方法包括以求解轨道根数为核心的解析几何法和以求解目标动力学微分方程为核心的数值积分法[6]。文献[6]分析了两种方法的预报特性。解析几何法将导弹视为二体运动,未考虑运动中的摄动因素,计算简单,预报快捷,但是预报精度不高。数值积分法考虑了导弹运动中的摄动因素,建立了复杂的动力学微分方程,预报精度高,但计算量大。
恩克法由德国科学家恩克(Johann F. Encke)提出,曾应用在计算短周期彗星和小行星轨道[7]。本文改进恩克法,并将其应用在雷达交接时刻预报中,首先,介绍了恩克法弹道预报算法的具体步骤,其次,对雷达交接时刻计算流程进行了详细分析,并基于无迹变换的方法进行了误差分析[8]。最后,对算法进行了仿真实验,并与传统方法进行了对比,验证了算法的有效性。

1 雷达交接时刻预报算法

1.1 雷达探测模型

从雷达站观测导弹飞行的观测数据中获得导弹的状态值,首先,需要建立雷达的观测模型。雷达站的位置通常使用经度、纬度、高度(L,B,H)来描述,而导弹目标的位置[X,Y,Z]T通常在地固系(ECEF)下进行描述,雷达对导弹目标位置的探测通常使用斜距、方位角、仰角(R,A,E)来确定,因此,需要进行坐标系间的转换,其转换过程如下。
1)将雷达站位置(L,B,H)转换为地固系坐标
(rR)ECEF= X Y Z= ( N + H ) · cos B · cos L ( N + H ) · cos B · sin L [ N ( 1 - e 2 ) + H ] · sin B
式中,N= a 1 - e 2 si n 2 B ,a为地球赤道半径,e为地球偏心率,L为雷达站的地理经度,B为雷达站的地理纬度,H为雷达站的大地高程。
2)将地固系下导弹对雷达的相对位置转换为雷达直角坐标系ENU下的坐标
x y z= T ECEF ENU X Y Z - ( r R ) ECEF= T ECEF ENU X - ( N + H ) · cos B · cos L Y - ( N + H ) · cos B · sin L Z - [ N ( 1 - e 2 ) + H ] · sin B
x · y · z ·= d dt T ECEF ENU X Y Z - ( r R ) ECEF= T ECEF ENU X · Y · Z ·
式中, T ECEF ENU为地固系与雷达直角坐标系的转换矩阵。
3)将导弹在雷达直角坐标系的坐标转化为雷达探测数据(R,A,E)
R A E= x 2 + y 2 + z 2 arctan x y arctan z x 2 + y 2

1.2 恩克法弹道预报方法

雷达交接时刻预报的重要环节是如何对导弹进行弹道预报。本文使用的恩克法弹道预报主要分为三个方面:利用二体力学模型建立导弹运动的基准轨道,求解基准轨道与实际轨道的偏差值,导弹状态在不同坐标系间的转换,具体流程如下。
1)求解基准轨道根数与t时刻的位置ρ
基准轨道是理想的二体轨道,可由6个独立的轨道根数(半长轴a、偏心率e、轨道倾角i、升交点赤经Ω、近地点幅角ω、真近点角f)来确定,轨道根数的求解则可由初始时刻的r0、v0算得[7]。
轨道确定后,轨道的偏近点角E,平近点角M,平均角速度n及过近地点时刻τ的计算如下
E = 2 arctan 1 - e 1 + e tan f 2 M = E - e sin E n = μ a - 3 τ = nt - M n
式中,μ为地球引力常数。由公式(5),可推导出预报时刻t的真近点角ft,进而,由ft计算t时刻下基准轨道位置ρ、速度 ρ ·
ρ = h 2 μ 1 1 + e cos f cos f t sin f t 0 ρ · = μ h - sin f t e + cos f t 0
式中,h= μa ( 1 - e 2 ) ,为动量矩常数值。
2)求解t时刻真实轨道与基准轨道的偏差量δr
本文分析导弹的受力情况,可得实际轨道和基准轨道的动力学微分方程分别为
r ¨+ μ r 3r=ap
ρ ¨+ μ ρ 3ρ=0
式中,ap为摄动力加速度表达式。导弹在自由段受到的主要摄动力为地球非球形引力,该摄动力中J2摄动影响最大,为 10-3量级,而J4项摄动与其他摄动力均为 10-6量级以下,因此通常只考虑J2摄动。
由基准轨道的定义可知δr=r-ρ,可得
δ r ¨= r ¨- ρ ¨
将式(7)、式(8)代入式(9)得到δr的微分方程δ r ¨=ap+ μ ρ 3 1 - ρ 3 r 3 r - δr。具体的表达式为
d dt δ r x δ r y δ r z δ v x δ v y δ v z= d dt δ r x δ r y δ r z δ r · x δ r · y δ r · z= δ r · x δ r · y δ r · z a px + μ ρ 3 1 - ρ 3 r 3 r x - δ r x a py + μ ρ 3 1 - ρ 3 r 3 r y - δ r y a pz + μ ρ 3 1 - ρ 3 r 3 r z - δ r z
因式(10)为非线性微分方程,无法获得解析解,故采用工程上常用的四阶龙格-库塔(R-K)方法计算其数值解。
3)计算t时刻导弹的状态量r
在获得t时刻导弹基准轨道的位置、速度及其偏差量后,t时刻真实轨道的位置、速度可由两者之和得到,具体表达式分别为:
r(t)= r x r y r z= ρ x ρ y ρ z+ δ r x δ r y δ r z
v(t)= r · x r · y r · z= ρ · x ρ · y ρ · z+ δ r · x δ r · y δ r · z
由于r(t)、v(t)是在地惯系下表示的,需要将其再统一到地固系下
X Y Z= T EC I ECEF r x r y r z
X · Y · Z ·= d dt T ECI ECEF r x r y r z= T ECI ECEF r · x r · y r · z
式中, T ECI ECEF为地惯系到地固系的坐标转换矩阵。
以上是恩克法预报的三个步骤,在步骤2)中,恩克法与传统的数值积分法一样,需要进行微分方程的求解,不同之处为:数值积分法积分求解的对象是导弹状态矢量,需要进行小时长积分多次迭代计算才能有较高精度;而恩克法积分求解的对象是实际轨道与基准轨道的摄动偏差量,是一个微小的量,在大时长积分的情况下,对整体位置的精度影响不大,从而可以进行单步积分,减小了整体的运算量节约了运算时间。恩克法弹道预报流程如图1所示。
图1 恩克法弹道预报流程图

1.3 雷达交接时刻预报流程

根据建立的雷达探测模型,使用弹道预报方法,就可以进行雷达交接时刻的预报。假设前置雷达稳定探测到导弹某一时刻的状态矢量(位置和速度),已知后置雷达的位置及威力空间γmax=(γmax, Amax,Emax),则雷达交接时刻的预报流程如下:
1)确定弹道预报的时间范围[t0,tf],其中t0为初始时刻,tf为导弹落点预报时刻,落点预报时刻可通过求解导弹的二体运动方程与地球理想球面方程计算获得;
2)将时间范围[t0,tf]进行2n(初始时令n=1)等分,得到时刻点 t1,t2,···, t 2 n + 1,使用弹道预报方法分别获得对应时间点的位置r1,r2,···, r 2 n + 1,进而,可计算出对应时刻点下,导弹与后置雷达的相对位置(γ1γ2··· γ 2 n + 1);
3)设t时刻雷达对导弹的观测量γ=(γ,A,E),规定:当且仅当γ<γmax,A<Amax,E<Emax都成立时,才有γ<γmax,否则为γ>γmax;若存在ti时刻满足γi<γmax,且ti-1时刻满足γi-1>γmax,说明预报交接点在区间[ti-1,ti]内(记为[tA,tB]),转入步骤4);若不存在ti时刻满足上述条件,令n=n+1,并转到步骤2);
4)使用弹道预报方法获得(tA+tB)/2时刻导弹的位置r,并转换为雷达与导弹的斜距γ,如果γ<γmax,则说明预报交接点在区间 t i - 1 , ( t i - 1 + t i ) / 2,并更新tB=(tA+tB)/2;如果γ>γmax,说明预报交接点在区间[(ti-1+ti)/2,ti],更新tA=(tA+tB)/2;
5)设ε为交接时刻计算的精确度,当|tA-tB|≥ε说明精确度还不够,需转到步骤4),当|tA-tB|<ε时,即满足所要求的精确度时,计算终止,得出最后的预报交接时刻为(tA+tB)/2,再次使用预报方法,获得交接时刻导弹的位置和速度。
雷达交接时刻的计算流程如图2所示。
图2 雷达交接时刻的预报流程图

1.4 基于无迹变换的交接时刻误差传播分析

在对预报时刻的误差传播情况进行分析时,需要计算预报时刻的状态协方差矩阵,现有的主要方法包括处理线性系统的求解雅克比矩阵方法、将非线性系统拟线性化的协方差分析描述函数法、基于无迹变换(Unscented Transform,UT)的协方差矩阵传播分析法。本文采用的是基于UT的协方差传播分析方法,该方法不需要求解偏导数矩阵,且不要求系统具有一阶可微性,适用于各类线性及复杂的非线性系统,与恩克法弹道预报融合后,其协方差传播的流程如图3所示。
图3 基于无迹变换的雷达交接时刻预报分析
设初始t0时刻导弹的状态矢量均值和协方差分别为 X –和PX,采用无迹变换计算交接时刻的均值和协方差过程如下。
1)按照如下样点选取规则,构造2n+1个σ点和相应的权重
σ ( i ) = X – , W i = κ n + κ , i = 0 . σ ( i ) = X – + ( ( n + κ ) P X ) i , 　　 W i = 1 2 ( n + κ ) , i = 1,2 , ··· , n σ ( n + i ) = X – - ( ( n + κ ) P X ) i , 　　 W n + i = 1 2 ( n + κ ) , i = 1,2 , ··· , n
式中,$\left(\sqrt{\boldsymbol{P}_{X}}\right)_{i}$为对协方差矩阵进行Cholesky分解获得的上三角阵的第i行;n为状态矢量维数,κ为标量参数,n+κ=3,Wi为对应的权重值。
2)由恩克法弹道预报方法计算交接时刻下对应样点的预报状态矢量ϕ
ϕ(i)=f(σ(i)),i=0,1,···,2n
式中,函数f(x)表示使用预报方法求交接时刻对应的状态矢量。
3)根据对应样点的预报值和相应的权重,计算交接时刻的状态均值 X – t和协方差 P X t
X – t = ∑ i = 0 2 n W i ϕ ( i ) P X t = ∑ i = 0 2 n W i [ ϕ ( i ) - X – t ] [ ϕ ( i ) - X – t ] T
式中,Wi为权重值,与公式(15)相同。
由雷达交接时刻的均值和协方差矩阵,可以在三维坐标系中以误差椭球和二维坐标系中以误差椭圆的形式,对均值和误差情况进行可视化直观展示。误差椭球(圆)的中心为预报均值,其长半轴、短半轴的大小由各个方向上的标准差和误差描述倍数决定,其中心轴与坐标轴的夹角可由协方差矩阵求得,具体求解公式可参阅文献[9]。

2 仿真分析

2.1 仿真参数设置

为验证算法,本文使用了Matlab进行仿真,仿真平台如下:计算机CPU为Intel Core i7-7700HQ(2.80 GHz),内存为8 GB,显卡为NVIDIA GeForce GTX 1050,仿真软件为Matlab2017b。
已知某弹道导弹飞行状态的实际数据,假设前置雷达跟踪探测到了导弹在t0=300 s的状态数据,此时时间为2013年5月4日8时5分0秒(协调世界时),其状态矢量(地固系位置和速度矢量,单位分别为m和m/s)的估计均值 x ˙ 0=(-4 768 663.5-2 824 156.1 3 969 707.2 2 982.790 8-4 584.489 9 3 644.722 7)T,对应的协方差矩阵为
P x 0= 10000 0 0 0 0 0 00 10000 0 0 0 0 0 0 10000 0 0 0 0 0 0 4 0 0 0 0 0 0 4 0 0 0 0 0 0 4
后置雷达的位置(纬度40°、西经80°、高度500 m),最大探测距离为1 600 km,仰角、方位角的范围不做限制,利用以上已知条件,对雷达交接时刻进行预报。
在与传统弹道预报方法计算的雷达交接时刻进行对比时,由于数值法与本文算法都使用了龙格-库塔积分方法,两者的计算时间、计算精度都与积分时长及积分次数有关,在预报时间相同的条件下,随着积分步数的增多,计算精度越高,但计算时间也会越长。为有效检验本文算法的效能,分析不同方法计算的精度与运算耗时及误差传播情况,综合考虑了积分步数和预报精度,恩克法采用一步积分,数值法采用三步积分。同时由于导弹运动速度高,将交接时刻计算的精确度设置到0.001 s。
根据真实数据,导弹进入后置雷达探测范围时刻(雷达交接时刻)的真实值为1 333.831 s,此时导弹的状态矢量xh=(-20 170.2-5 679 330.5 5 109 216.7 5 310.1-755.9-1 356.6)T。

2.2 仿真结果及分析

根据上述实验条件设置,本文分别使用恩克法和传统的解析法、数值法对雷达交接时刻及其相应的位置速度进行预报,并使用无迹变换的方法,计算交接时刻的状态协方差矩阵,三种算法的运算时间采用运算100次取平均值获得。
三种方法计算的雷达交接时刻值及与真实值的误差如表1所示。对比可知本文恩克法预报的交接时刻的预报精度最高,误差仅为0.115 s。解析法与数值法精度相当,误差分别为0.388 s、0.389 s。
表1 三种方法预报雷达交接时刻值及误差
解析法 数值法
(3步长)		 恩克法
(1步长)
交接时刻预报值/s 1 333.443 1 333.442 1 333.716
与实际值误差/s 0.388 0.389 0.115
三种方法计算交接时刻时的位置速度均值及真实值如表2所示。由表可知,Y轴方向为三个方向中偏离最大,其中解析法偏离大于3 km,数值法偏离大于1 km,恩克法偏离近1 km。说明解析法预报的状态均值误差较大,数值法与恩克法预报的位置精度较高且精度相当。
表2 三种方法预报的位置速度均值及与真实值对比
X/m Y/m Z/m Vx/(m·s-1) Vy/(m·s-1) Vz/(m·s-1)
解析法 -19 694.436 -5 676 148.030 5 112 253.22 5 313.699 -751.812 -1 351.383
数值法 -21 849.155 -5 678 021.154 5 108 882.084 5 311.109 -756.738 -1 355.566
恩克法 -18 255.850 -5 680 105.118 5 110 173.812 5 313.377 -757.804 -1 355.178
真实值 -20 170.260 -5 679 330.179 5 109 216.450 5 310.179 -755.961 -1 356.642
三种方法使用无迹变换计算交接时刻的位置速度标准差结果如表3所示。分析可知,三种方法在不同方向上计算的位置标准差都相差不大,其最大差值仅为4 m,说明三种方法的误差传播发散程度相当。
表3 三种方法在交接时刻下位置速度的标准差
X/m Y/m Z/m Vx/(m·s-1) Vy/(m·s-1) Vz/(m·s-1)
解析法 1 961.080 2 201.790 2 307.433 1.622 2.874 3.001
数值法 1 959.776 2 200.381 2 307.943 1.621 2.869 3.001
恩克法 1 962.787 2 196.821 2 311.710 1.623 2.865 3.005
三种方法的运算耗时如表4所示。对比可知,解析法耗时最短,恩克法次之,数值法耗时最长。原因在于解析法运算只进行了一次二体预报,其计算耗时短,恩克法较解析法多了一次积分运算,而数值法则需要进行三次积分计算其耗时最长。
表4 三种预报方法运算耗时
解析法 数值法(3步长) 恩克法(1步长)
计算耗时/s 0.037 0.138 0.082
为直观地展示三种方法在交接时刻的误差传播情况,本文以蒙特卡洛(Monte Carlo)打靶仿真1 000次获得的交接时刻散点分布为基准,分别与三种方法计算的交接时刻三维误差椭球及两个方向的二维误差椭圆进行对比,用来验证算法的准确性和有效性,计算结果比较如图4、图5所示。图中,圆圈表示1 000次蒙特卡洛打靶点,实线表示解析法的计算结果,虚线表示数值法的计算结果,点划线表示恩克法的计算结果。由图中可以看出,三种方法中,解析法误差椭球(圆)离仿真打靶点最远,位置精度最差,数值法、恩克法位置预报精度相当,与前文分析一致。
图4 仿真打靶(1 000次)交接时刻位置与误差椭球三维分布情况
图5 仿真打靶(1 000次)交接时刻位置与误差椭圆二维分布情况
综上分析,解析法虽然耗时最短,但是,其预报的交接时刻及状态均值的精度最低。数值法虽然预报的状态均值精度与恩克法相当,但交接时刻的时间预报精度低,存在时间与空间不完全匹配的问题,同时运算耗时最长。恩克法对比解析法考虑了摄动因素影响,对比数值法能够进行大步长积分,其运算耗时较短,交接时刻及状态均值的计算精度最高。

3 结束语

本文介绍了导弹预警中的雷达探测模型及恩克法弹道预报的具体步骤,对雷达交接时刻计算的流程进行了详细说明,使用无迹变换对预报交接时刻的误差情况进行了分析。实验表明,与传统的解析几何法和数值积分法相比,本文算法考虑了导弹运动中受到的摄动力影响,能够进行大步长的积分运算,在相对较低的运算时间内,提高了雷达交接时刻的预报精度。

References

Publishing order | Descend order by publishing year | Descend order by cited within
[1]
刘蛟龙, 尉建利, 于云峰. 弹道导弹自由段拦截窗口建模与仿真[J]. 指挥控制与仿真, 2012, 34(6):72-75.

[2]
毛艺帆, 张多林, 王路. 基于轨道根数交接的弹道导弹落点快速预报[J]. 空军工程大学学报(自然科学版), 2017, 18(1):33-38.

[3]
Lee D G, Cho K S, Shin J H. A Simple Prediction Method of Ballistic Missile Trajectory to Designate Search Direction and its Verification Using a Test Bench[C]. 10th Asian Control Conference (ASCC), Sabah, 2015.

[4]
蔺美青, 高玉良, 阮菲. 防空反导作战雷达交接方案优化设计方法研究[J]. 空军雷达学院学报, 2012, 26(2):119-123.

[5]
张荣涛. 多雷达跟踪弹道导弹交接预报技术研究[J]. 现代雷达, 2010, 32(8):31-32.

[6]
孙瑜, 孟凡坤, 吴楠, 等. 解析几何法和数值积分法的弹道预报性能分析[J]. 信息工程大学学报, 2016, 17(5):635-640.

[7]
张洪波. 航天器轨道力学理论与方法[M]. 北京: 国防工业出版社, 2015.

[8]
Angrisani L, D"Apuzzo M, Moriello R S L. Unscented Transform: a Powerful Tool for Measurement Uncertainty Evaluation[J]. IEEE Transactions on Instrumentation and Measurement, 2006, 55(3):737-743.

DOI

[9]
吴楠, 王锋, 孟凡坤. 无再入观测弹道导弹气动参数和落点联合预报[J]. 弹道学报, 2018, 30(3):18-24.

Options
Outlines

/

〈 〉