CN105022035A - 一种基于模型修正的弹道目标发射点估计装置及其方法 - Google Patents

Abstract

本发明公开了基于模型修正的弹道目标发射点估计装置及其方法，本发明通过对雷达测量的点迹数据进行反向跟踪滤波，并采用龙格库塔外推方法计算弹道导弹的飞行轨迹和发射点位置，利用模型偏差对估计的发射点位置进行修正。根据本发明公开的方法，采用基于弹道导弹运动方程的不敏滤波算法，实现了高精度的非线性滤波，且计算量适中，具有较高的费效比；通过4阶龙格库塔方法进行轨道外推计算，具有较高的数值精度；通过模型偏差修正的方法，去除了主动段所引起的发射点位置估计偏差，大大提高了弹道导弹发射点估计的精度。

Description

一种基于模型修正的弹道目标发射点估计装置及其方法

技术领域

本发明涉及一种弹道目标发射点估计装置、弹道目标发射点估计方法，尤其涉及一种基于模型修正的弹道目标发射点估计装置、基于模型修正的弹道目标发射点估计装置方法。

背景技术

弹道导弹具有射程远、威力大、精度高、机动性强等优越性，成为现代战争中的“撒手锏”。二战后世界各国开始竞相研制各种导弹武器，弹道导弹将成为战争初期或关键时刻的主要作战方式。现阶段，弹道导弹的各种突防技术已日趋成熟，仅在我国周边就有超过十余个国家和地区，已经掌握和即将拥有中程以上的弹道导弹，这对我国国家安全构成严重威胁。

雷达作为导弹防御系统中的核心探测器，在导弹防御的预警探测、跟踪、识别、制导以及杀伤效果评估等各关键环节均发挥着不可替代的关键作用。雷达性能的优劣对于整个导弹防御系统的性能都有很大的影响，而通过雷达的探测，对弹道导弹发射点估计的问题是导弹防御体系中至关重要的环节，该问题一直是理论和实际研究的重点和热点。

下面结合图1详细介绍现有的弹道目标发射点估计方法通常的处理流程。

步骤101：对雷达测量的点迹数据进行预处理，完成野值剔除，坐标转换功能。

步骤102：对预处理后的点迹数据进行平滑/滤波处理，常用的方法有最小二乘平滑(LS)和扩展卡尔曼滤波(EKF)方法。

步骤103：通过平滑/滤波后的数据对弹道目标采用二体问题公式进行定轨，计算目标的轨道根数。

步骤104：通过弹道目标的轨道根数，由于目标的弹道面是过地球地心的平面且该弹道上的发射点和落点与地球表面相交，满足地球曲面方程；可利用解析法直接按公式解算出发射点位置。

现有技术中的弹道目标发射点估计方法，都存在一些缺点，概括如下：

(1)现有技术通常采用的最小二乘平滑方法或扩展卡尔曼滤波方法；其中最小二乘平滑方法未考虑弹道目标运动模型，滤波的精度不高，且不适用于实时处理，多作为事后批处理分析。扩展卡尔曼滤波在非线性滤波问题上易于发散，且滤波的精度不高。

(2)基于二体问题进行的椭圆定轨方法，解算精度不够，同时该方程忽略了地球自转，地表高程等因素影响，而这对于中程及以上的弹道导弹发射点估计精度的影响是不可忽略的。

(3)现有技术发射点估计时忽略了弹道导弹主动段的影响，采用二体问题处理，将弹道目标飞行的全程建模为自由段，这忽略了主动段所引起的发射点估计偏差，导致估计结果偏差较大。

发明内容

本发明的目的在于提供一种基于模型修正的弹道目标发射点估计装置、基于模型修正的弹道目标发射点估计装置方法，综上所述，现有弹道目标发射点估计方法在数据平滑/滤波精度，弹道目标定轨和外推，主动段引起的估计偏差等方面还要进一步改进。因此，本发明提出的一种基于模型修正的弹道目标发射点估计装置及其方法具有及其重要的意义。

本发明的解决方案是：一种基于模型修正的弹道目标发射点估计装置，其包括：

点迹数据预处理模块，所述点迹数据预处理模块用于对点迹数据进行预处理，雷达量测点迹数据为Z＝{z₁,…,z_k}；

时间逆序处理模块，其用于对雷达量测点迹数据Z＝{z₁,…,z_k}进行时间逆序处理，转换为Z'＝{z_k,z_k-1,…,z₁}，并去除时间重复的点迹数据；

反向滤波处理模块，其用于对雷达量测点迹的时间逆序数据Z'＝{z_k,z_k-1,…,z₁}，利用基于弹道目标运动方程进行反向滤波处理；

反向轨道外推计算模块，其用于利用反向滤波的估计值，采用基于弹道目标运动方程的龙格库塔4阶方法进行反向轨道外推计算，得到反向轨道外推结果X_launch＝<X₁,Y₁,Z₁>_ecef；

发射点位置计算模块，其用于计算弹道目标发射点位置，利用X_launch＝<X₁,Y₁,Z₁>_ecef从地心地固坐标系下转换到WGS84地心坐标系<L₁,B₁,H₁>；

正向滤波处理模块，其用于对雷达量测点迹数据Z＝{z₁,…,z_k}，利用基于弹道目标运动方程进行正向滤波处理；

正向轨道外推计算模块，其用于利用正向滤波的估计值，采用基于弹道目标运动方程的4阶龙格库塔方法进行正向轨道外推，得到正向轨道外推结果X_impact＝<X₂,Y₂,Z₂>_ecef；

最大飞行高度和射程计算模块，其用于计算弹道目标最大飞行高度和射程，通过X_launch＝<X₁,Y₁,Z₁>_ecef和X_impact＝<X₂,Y₂,Z₂>_ecef统计得出弹道目标的最大飞行高度H_max，计算弹道目标落点估计位置，将从地心地固坐标系下转换到WGS84地心坐标系<L₂,B₂,H₂>，通过发射点估计位置X_launch＝<X₁,Y₁,Z₁>_ecef和落点估计位置X_impact＝<X₂,Y₂,Z₂>_ecef，计算出弹道目标的射程D；

模型偏差查询模块，其用于通过弹道目标最大飞行高度H_max和射程D到弹道目标发射点模型偏差数据库中查询模型偏差ΔX；

模型偏差修正模块，其用于利用基于模型偏差修正发射点方法计算修正后的弹道目标发射点位置<L_launch,B_launch,H_launch>。

本发明还提供一种基于模型修正的弹道目标发射点估计方法，其包括以下步骤：

步骤201，雷达量测点迹数据为Z＝{z₁,…,z_k}，对点迹数据进行预处理；

步骤202，对雷达量测点迹数据Z＝{z₁,…,z_k}进行时间逆序处理，转换为Z'＝{z_k,z_k-1,…,z₁}，并去除时间重复的点迹数据；

步骤203，对雷达量测点迹的时间逆序数据Z'＝{z_k,z_k-1,…,z₁}，利用基于弹道目标运动方程进行反向滤波处理；

步骤204，利用反向滤波的估计值，采用基于弹道目标运动方程的龙格库塔4阶方法进行反向轨道外推计算，得到反向轨道外推结果X_launch＝<X₁,Y₁,Z₁>_ecef；

步骤205，计算弹道目标发射点位置，利用X_launch＝<X₁,Y₁,Z₁>_ecef从地心地固坐标系下转换到WGS84地心坐标系<L₁,B₁,H₁>；

步骤206，对雷达量测点迹数据Z＝{z₁,…,z_k}，利用基于弹道目标运动方程进行正向滤波处理；

步骤207，利用正向滤波的估计值，采用基于弹道目标运动方程的4阶龙格库塔方法进行正向轨道外推，得到正向轨道外推结果X_impact＝<X₂,Y₂,Z₂>_ecef；

步骤208，计算弹道目标最大飞行高度和射程，通过X_launch＝<X₁,Y₁,Z₁>_ecef和X_impact＝<X₂,Y₂,Z₂>_ecef统计得出弹道目标的最大飞行高度H_max，计算弹道目标落点估计位置，将从地心地固坐标系下转换到WGS84地心坐标系<L₂,B₂,H₂>，通过发射点估计位置X_launch＝<X₁,Y₁,Z₁>_ecef和落点估计位置X_impact＝<X₂,Y₂,Z₂>_ecef，计算出弹道目标的射程D；

步骤209，通过弹道目标最大飞行高度H_max和射程D到弹道目标发射点模型偏差数据库中查询模型偏差ΔX；

步骤210，利用基于模型偏差修正发射点方法计算修正后的弹道目标发射点位置<L_launch,B_launch,H_launch>。

作为上述方案的进一步改进，反向滤波处理与正向滤波处理均采用不敏滤波算法。

优选地，步骤203与步骤206中基于弹道目标运动方程的不敏滤波方法包括以下步骤：

一、假定目标的状态向量为弹道导弹的运动模型为：

$\left[\begin{array}{c}\stackrel{\&CenterDot;}{x}\\ \stackrel{\&CenterDot;}{y}\\ \stackrel{\&CenterDot;}{z}\\ \stackrel{\&CenterDot;\&CenterDot;}{x}\\ \stackrel{\&CenterDot;\&CenterDot;}{y}\\ \stackrel{\&CenterDot;\&CenterDot;}{z}\end{array}\right]=\left[\begin{array}{c}\stackrel{\&CenterDot;}{x}\\ \stackrel{\&CenterDot;}{y}\\ \stackrel{\&CenterDot;}{z}\\ -\frac{\mathrm{\&rho;}\left(h\right)}{2\mathrm{\&beta;}}vx-\frac{\mathrm{\&mu;}G}{{\left(r\right)}^3}x+(2{\mathrm{\&omega;v}}_y+{\mathrm{\&omega;}}^{2}x)\\ -\frac{\mathrm{\&rho;}\left(h\right)}{2\mathrm{\&beta;}}vy-\frac{\mathrm{\&mu;}G}{{\left(r\right)}^3}y+(-2{\mathrm{\&omega;v}}_x+{\mathrm{\&omega;}}^{2}y)\\ -\frac{\mathrm{\&rho;}\left(h\right)}{2\mathrm{\&beta;}}{\mathrm{vv}}_z-\frac{\mathrm{\&mu;}G}{{\left(r\right)}^3}z\end{array}\right]$     公式9，

其中ρ(h)＝ρ₀e^-kh为地球大气密度函数，其中ρ₀＝1.22kg/m³，k＝0.14141×10^-3m^-1，h为目标的海拔高度；β为弹道系数；μ_G＝3.986005×10¹⁴m³/s²为地球万有引力常量；ω＝7.292116×10^-5rad/s为地球自转速度；

二、雷达测量模型h为：

$h\left(X\right)=\left[\begin{array}{c}\sqrt{x^2+y^2+z^2}\\ \arctan (x/y)\\ arc\tan (z/\sqrt{x^2+y^2})\end{array}\right]$       公式10，

假定雷达量测数据的噪声δ为不相关的零均值高斯白噪声，其中距离噪声的方差为σ_r，方位角噪声的方差为σ_α，仰角噪声的方差为σ_β；

三、目标的状态向量为β为弹道系数，则状态向量维数L＝7，不敏滤波的处理过程为：首先利用不敏变换生成采样点，采用常用的标准对称采样方法，在目标的运动模型和雷达量测模型是非线性函数，且目标在k-1时刻的估计均值为和P_k-1/k-1时，选择2L+1＝15个采样点按式11计算：

$\left\{\begin{array}{c}{\mathrm{\&chi;}}_0={\hat{X}}_{k-1/k-1},W_0^{\left(m\right)}=\frac{\mathrm{\&lambda;}}{2(L+\mathrm{\&lambda;})},W_0^{\left(c\right)}=\frac{\mathrm{\&lambda;}}{2(L+\mathrm{\&lambda;})}+1-{\mathrm{\&alpha;}}^2+\mathrm{\&gamma;};\\ {\mathrm{\&chi;}}_i={\hat{X}}_{k-/k-1}+{\left(\sqrt{(L+\mathrm{\&lambda;})P_{k-1/k-1}}\right)}_i,W_i^{\left(m\right)}=W_i^{\left(c\right)}=\frac{1}{2(L+\mathrm{\&lambda;})},i=1,\mathrm{...},L;\\ {\mathrm{\&chi;}}_{i+L}={\hat{X}}_{k-1/k-1}-{\left(\sqrt{(L+\mathrm{\&lambda;})P_{k-1/k-1}}\right)}_i,W_{i+L}^{\left(m\right)}=W_{i+L}^{\left(c\right)}=\frac{1}{2(L+\mathrm{\&lambda;})},i=1,\mathrm{...},L;\end{array}\right.$     公式11，

其中λ＝α²(L+κ)为比例参数，作为控制采样点到均值的距离；一般取α＝0.5；κ＝3-L；γ＝2；式中为(L+λ)P_k-1/k-1均方根矩阵的第i行；

对时刻k状态的一部提前预测，计算目标状态预测：

    公式12，χ_k/k-1＝f(χ_k/k-1,k-1)  公式13；

计算目标预测协方差矩阵：

$P_{k|k-1}=\underset{i=0}{\overset{2n}{\&Sigma;}}(W_i^{\left(c\right)}\&lsqb;{\mathrm{\&chi;}}_{i,k/k-1}-{\hat{X}}_{k/k-1}\&rsqb;{\&lsqb;{\mathrm{\&chi;}}_{i,k/k-1}-{\hat{X}}_{k/k-1}\&rsqb;}^T)+Q_{k-1}$    公式14；

计算时刻k量测的预测：

${\hat{Z}}_{k/k-1}=\underset{i=0}{\overset{2n}{\&Sigma;}}{W}_i{\mathrm{\&zeta;}}_{i,k/k-1}$      公式15，

ζ_k/k-1＝h(χ_k/k-1,k)    公式16，

$P_{{\tilde{z}}_k}=\underset{i=0}{\overset{2n}{\&Sigma;}}{W}_i^{\left(m\right)}({\mathrm{\&zeta;}}_{k/k-1}^i-{\hat{z}}_{k|k-1}){({\mathrm{\&zeta;}}_{k|k-1}^i-{\hat{z}}_{k|k-1})}^T+R_k$      公式17，

$P_{{\tilde{x}}_k{\tilde{z}}_k}=\underset{i=0}{\overset{2n}{\&Sigma;}}{W}_i^{\left(m\right)}({\mathrm{\&zeta;}}_{k|k-1}^i-{\hat{x}}_{k|k-1}){({\mathrm{\&zeta;}}_{k/k-1}^i-{\hat{x}}_{k|k-1})}^T$      公式18；

计算新息：

$v_k=Z_k-{\hat{Z}}_{k/k-1}$       公式19；

计算增益矩阵：

$K_k=P_{{\tilde{x}}_k{\tilde{z}}_k}{P}_{{\tilde{z}}_k}^{-1}$      公式20；

计算目标状态更新：

${\hat{X}}_{k|k}={\hat{X}}_{k|k-1}+K_k{v}_k$        公式21；

计算目标状态误差协方差的更新：

${\hat{P}}_{k|k}={\hat{P}}_{k|k-1}-K_k{P}_{{\tilde{z}}_k}{K}_k^T$      公式22。

作为上述方案的进一步改进，X_launch＝<X₁,Y₁,Z₁>_ecef从地心地固坐标系下转换到WGS84地心坐标系<L₁,B₁,H₁>、以及X_impact＝<X₂,Y₂,Z₂>_ecef从地心地固坐标系下转换到WGS84地心坐标系<L₁,B₁,H₁>均通过公式6得到：

$r=\sqrt{x^2+y^2}$

E²＝a²-b²

F＝54b²z²

G＝r²+(1-e²)z²-e²E²

$C=\frac{e^4{\mathrm{Fr}}^2}{G^3}$

$S=\sqrt[3]{1+C+\sqrt{C^2+2C}}$

$P=\frac{F}{3{(S+\frac{1}{S}+1)}^2{G}^2}$

$Q=\sqrt{1+2e^{4}P}$

$r_0=\frac{-\left({\mathrm{Pe}}^{2}r\right)}{1+Q}+\sqrt{\frac{1}{2}{a}^2(1+\frac{1}{Q})-\frac{P(1-e^2)Z^2}{Q(1+Q)}-\frac{1}{2}{\mathrm{Pr}}^2}$

$U=\sqrt{{(r-e^2{r}_0)}^2+Z^2}$

$V=\sqrt{{(r-e^2{r}_0)}^2+(1-e^2)Z^2}$

$Z_0=\frac{b^{2}Z}{aV}$       公式6；

$H_1=U(1-\frac{b^2}{aV})$

$B_1=arctan\&lsqb;\frac{Z+e^2{Z}_0}{r}\&rsqb;$

L₁＝arctan2(Y,X)

其中，x，y，z分别表示地心地固坐标系中的三个坐标轴；

L，B，H分别表示雷达站心的经度、纬度和高度；

a＝6378137m为地球长半轴，b＝6356752.3142m为地球短半轴，e＝0.006694380为地球椭球第一偏心率；

E，F，G，C，S，P，Q，r₀，U，V，Z₀分别为临时计算参数。

作为上述方案的进一步改进，所述最大飞行高度和射程计算模块通过式7计算出弹道目标的射程D：

$f_{dis\tan ce}(\left[\begin{array}{c}{X}_1\\ Y_1\\ Z_1\end{array}\right],\left[\begin{array}{c}{X}_2\\ Y_2\\ Z_2\end{array}\right])=R_e*\mathrm{\&beta;}$

$r_1=\sqrt{{X_1}^2+{Y_1}^2+{Z_1}^2}$       公式7，

$r_2=\sqrt{{X_2}^2+{Y_2}^2+{Z_2}^2}$

$\mathrm{\&beta;}=arccos\left(\frac{(X_1{X}_2+Y_1{Y}_2+Z_1{Z}_2)}{r_1{r}_2}\right)$

其中，<X₁,Y₁,Z₁>_ecef为地心地固坐标系下弹道目标发射点的位置参数，<X₂,Y₂,Z₂>_ecef为地心地固坐标系下弹道目标落地的位置参数，R_e为地球半径，β为方位，r₁，r₂分别为临时计算参数。

作为上述方案的进一步改进，步骤204与步骤207中基于弹道目标运动方程的4阶龙格库塔轨道外推方法包括以下步骤：

采用4阶龙格库塔方法进行外推，所用的计算增函数公式的形式为：

Δy_n＝y_n+1-y_n＝αk₁+βk₂+γk₃+δk₄      公式23，

其中

$\left\{\begin{array}{c}{k}_1=y^{\&prime;}(t_n,y_n)h\\ k_2=y^{\&prime;}(t_n+h/2,y_n+k_1/2)h\\ k_3=y^{\&prime;}(t_n+h/2,y_n+k_2/2)h\\ k_4=y^{\&prime;}(t_n+h/2,y_n+k_3/2)h\end{array}\right.$       公式24；

在满足各对应项系数相等的条件下选取未定系数α＝h/6，β＝h/3，γ＝h/3，δ＝h/6的适当值，y′(t_n,y_n)为弹道导弹运动方程见公式9，然后将这些值代入上式可得

$y_{n+1}-y_n=\frac{h}{6}(k_1+2k_2+2k_3+k_4)$      公式25。

作为上述方案的进一步改进，步骤209中构建弹道目标发射点模型偏差数据库方法采用基于仿真计算构建模型偏差库：

依据射程和最大飞行高度选择典型的弹道导弹主动段模型进行仿真：

a)假设发点位置O₁，仿真计算出关机点的位置和速度[r,v]；

b)通过关机点的位置和速度[r,v]利用自由段模型反推得到步骤205计算的结果即发射点估计值O₂；

c)依据发点O₁和O₂，计算模型偏差ΔX＝f_distance(O₁,O₂)，其中f_distance()为计算两点的距离函数为公式7。

在模型偏差数据库中记录该仿真的条件：主动段模型类型，射程最大飞行高度和计算的结果模型偏差ΔX。

作为上述方案的进一步改进，步骤209中构建弹道目标发射点模型偏差数据库方法采用基于实验数据构建模型偏差库：

通过已知真实的发射点位置O₁，计算的发射点估计O₂，计算模型偏差ΔX＝f_distance(O₁,O₂)，在模型偏差数据库中记录该仿真的条件：导弹型号，射程最大飞行高度和计算的结果模型偏差ΔX。

作为上述方案的进一步改进，基于模型偏差修正发射点方法包括以下步骤：

目标是估计发点的真实位置需要对发射点估计结果修正模型偏差ΔX，弹道目标的外推轨迹在地表的投影，即为发射平面与地球球面的交线，沿该投影线修正模型偏差ΔX可得到

${\hat{X}}_O={\hat{X}}_{O^{\&prime;}}-\mathrm{\&Delta;}X$       公式26，

其中模型偏差ΔX依赖于弹道导弹的模型，通过弹道目标发射点模型偏差数据库中查询获得。

与现有技术相比，本发明的有益效果如下：

(1)本发明针对弹道目标的运动方程进行建模，采用不敏滤波算法进行滤波处理。弹道目标的跟踪滤波是典型的非线性滤波问题，现有技术中采用的最小二乘/扩展卡尔曼滤波方法滤波精度不高，本发明中采用不敏滤波算法的滤波精度可达到二阶扩展卡尔曼滤波的精度，这对提高弹道目标的定轨精度和发射点位置估计精度都是非常重要的前提。

(2)本发明在进行发射点估计中采用基于地表高度的4阶龙格库塔数值积分方法，相对于传统方法采用了精度更高的数值解算方法，同时也考虑了不同地形具有不同地表的高度，进一步提高了发射点估计的精度。

(3)本发明利用对发射点估计位置运用模型偏差来消除弹道导弹主动段所引起的估计偏差，避免了现有技术直接忽略主动段所引起的误差，修正后的发射点位置大大提高了估计的精度。

附图说明

图1是现有技术中弹道目标发射点估计的流程图；

图2是本发明的弹道目标发射点估计的流程图；

图3是本发明中典型的弹道导弹飞行轨迹示意图；

图4是本发明实施例的弹道目标发射点估计原理示意图。

具体实施方式

以下结合附图以及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不限定本发明。

本发明的基于模型修正的弹道目标发射点估计装置，其包括点迹数据预处理模块、时间逆序处理模块、反向滤波处理模块、反向轨道外推计算模块、发射点位置计算模块、正向滤波处理模块、正向轨道外推计算模块、最大飞行高度和射程计算模块、模型偏差查询模块、模型偏差修正模块。

所述点迹数据预处理模块用于对点迹数据进行预处理，雷达量测点迹数据为Z＝{z₁,…,z_k}Z＝{z₁,…,z_k}。时间逆序处理模块用于对雷达量测点迹数据Z＝{z₁,…,z_k}进行时间逆序处理，转换为Z'＝{z_k,z_k-1,…,z₁}，并去除时间重复的点迹数据。反向滤波处理模块用于对雷达量测点迹的时间逆序数据Z'＝{z_k,z_k-1,…,z₁}，利用基于弹道目标运动方程进行反向滤波处理。反向轨道外推计算模块用于利用反向滤波的估计值，采用基于弹道目标运动方程的龙格库塔4阶方法进行反向轨道外推计算，得到反向轨道外推结果X_launch＝<X₁,Y₁,Z₁>_ecef。发射点位置计算模块，其用于计算弹道目标发射点位置，利用X_launch＝<X₁,Y₁,Z₁>_ecef从地心地固坐标系下转换到WGS84地心坐标系<L₁,B₁,H₁>。正向滤波处理模块用于对雷达量测点迹数据Z＝{z₁,…,z_k}，利用基于弹道目标运动方程进行正向滤波处理。正向轨道外推计算模块用于利用正向滤波的估计值，采用基于弹道目标运动方程的4阶龙格库塔方法进行正向轨道外推，得到正向轨道外推结果X_impact＝<X₂,Y₂,Z₂>_ecef。最大飞行高度和射程计算模块用于计算弹道目标最大飞行高度和射程，通过X_launch＝<X₁,Y₁,Z₁>_ecef和X_impact＝<X₂,Y₂,Z₂>_ecef统计得出弹道目标的最大飞行高度H_max，计算弹道目标落点估计位置，将从地心地固坐标系下转换到WGS84地心坐标系<L₂,B₂,H₂>，通过发射点估计位置X_launch＝<X₁,Y₁,Z₁>_ecef和落点估计位置X_impact＝<X₂,Y₂,Z₂>_ecef，计算出弹道目标的射程D。模型偏差查询模块用于通过弹道目标最大飞行高度H_max和射程D到弹道目标发射点模型偏差数据库中查询模型偏差ΔX。模型偏差修正模块用于利用基于模型偏差修正发射点方法计算修正后的弹道目标发射点位置<L_launch,B_launch,H_launch>。

图2是本发明实施例的弹道目标发射点估计的流程图，如图2所示，本发明实施例的弹道目标发射点估计方法的工作流程包括：

步骤201：雷达量测点迹数据为Z＝{z₁,…,z_k}，对点迹数据进行预处理，包括：

(1)对雷达测量点迹数据进行野值剔除；对点迹数据的距离、方位、俯仰、径向速度进行逐个检测，采用3-Sigma的方法判别：其中x_k为k时刻的测量值，为该时刻的预测值，σ_x为测量误差，如差值超过3倍的测量误差则判别为野值。

(2)将雷达测量点迹从雷达站心的球坐标系转换到地心地固坐标系。

雷达量测数据在以雷达站为中心的极坐标系下获得的量测Z＝<r,β,α>，其中距离r，方位β，仰角α，雷达的测量误差分别为σ_r、σ_θ和σ_α(本实施例中σ_r＝50m，σ_θ＝0.1deg，σ_η＝0.1deg)。

第一步将点迹从雷达站心的球坐标系Z＝<r,β,α>转换到雷达站心东北天坐标系X_enu＝<e,n,u>：

e＝rcos(α)sin(β)

n＝rcos(α)cos(β)           (1)

u＝rsin(α)

第二步将点迹从雷达站心东北天坐标系X_enu＝<e,n,u>转换到地心地固坐标系X_ecef＝<x,y,z>，假设雷达站心的经度、纬度和高度为<L,B,H>，则雷达站心东北天坐标系到地心地固坐标系的转换矩阵如下：

$M=\left[\begin{array}{ccc}-\sin L& -\sin B\cos L& \cos B\cos L\\ \cos L& -\sin B\sin L& \cos B\sin L\\ 0& \cos B& \sin B\end{array}\right]---\left(2\right)$

计算雷达站心在地心地固坐标系的位置X₀＝<x,y,z>：

$\begin{array}{c}{x}_0=(N+H)\cos B\cos L\\ y_0=(N+H)\cos B\sin L\\ z_0=\left(N\right(1-e_1^2)+H)\sin B\end{array}---\left(3\right)$

其中

$N=\frac{a}{\sqrt{1-e_1^2{\sin }^{2}B}}---\left(4\right)$

且a＝6378137m为地球长半轴，b＝6356752.3142m为地球短半轴，为地球第一偏心率。

将点迹转换到地心地固坐标系X_ecef＝<x,y,z>：

$\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=M\left(\begin{array}{c}e\\ n\\ u\end{array}\right)+\left(\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right)---\left(5\right)$

步骤202：对雷达量测点迹数据Z＝{z₁,…,z_k}进行时间逆序处理，转换为Z'＝{z_k,z_k-1,…,z₁}，并去除时间重复的点迹数据。

步骤203：对雷达量测点迹的时间逆序数据Z'＝{z_k,z_k-1,…,z₁}，利用基于弹道目标运动方程进行反向滤波处理，滤波方法采用不敏滤波算法。

步骤204：利用滤波的估计值，采用基于弹道目标运动方程的龙格库塔4阶方法进行反向轨道外推计算，考虑地表高度。

步骤205：计算弹道目标发射点位置，利用外推结果X_launch＝<X₁,Y₁,Z₁>_ecef转换到WGS84地心坐标系<L₁,B₁,H₁>，这里给出一种确定解的方法：

$\begin{array}{c}r=\sqrt{x^2+y^2}\\ E^2=a^2-b^2\\ F=54b^2{z}^2\\ G=r^2+(1-e^2)z^2-e^2{e}^2\\ C=\frac{e^4{\mathrm{Fr}}^2}{G^3}\\ S=\sqrt[3]{1+C+\sqrt{C^{2}2C}}\\ P=\frac{F}{3{(S+\frac{1}{S}+1)}^2{G}^2}\\ Q=\sqrt{1+2e^{4}P}\\ r_0=\frac{-\left({\mathrm{Pe}}^{2}r\right)}{1+Q}+\sqrt{\frac{1}{2}{a}^2(1+\frac{1}{Q})-\frac{P(1-e^2)Z^2}{Q(1+Q)}-\frac{1}{2}{\mathrm{Pr}}^2}\\ U=\sqrt{{(r-e^2{r}_0)}^2+Z^2}\\ V=\sqrt{{(r-e^2{r}_0)}^2+(1-e^2)Z^2}\\ Z_0=\frac{b^{2}Z}{aV}\\ H_1=U(1-\frac{b^2}{aV})\\ B_1=\arctan \&lsqb;\frac{Z+e^2{Z}_0}{r}\&rsqb;\\ L_1=\arctan 2(Y,X)\end{array}---\left(6\right)$

步骤206：对雷达量测点迹数据Z＝{z₁,…,z_k}，利用基于弹道目标运动方程进行正向滤波处理，滤波方法采用不敏滤波算法。

步骤207：利用滤波的估计值，采用基于弹道目标运动方程的4阶龙格库塔方法进行正向轨道外推，考虑地表高度。

步骤208：计算弹道目标最大飞行高度和射程；通过反向外推轨道和正向外推轨道可统计得出弹道目标的最大飞行高度H_max。计算弹道目标落点估计位置，利用公式(6)将外推结果X_impact＝<X₂,Y₂,Z₂>_ecef从地心地固坐标系下转换到WGS84地心坐标系<L₂,B₂,H₂>。通过发射点估计位置X_launch＝<X₁,Y₁,Z₁>_ecef和落点估计位置X_impact＝<X₂,Y₂,Z₂>_ecef，通过式7可计算出弹道目标的射程D：

$f_{dis\tan ce}(\left[\begin{array}{c}{X}_1\\ Y_1\\ Z_1\end{array}\right],\left[\begin{array}{c}{X}_2\\ Y_2\\ Z_2\end{array}\right])=R_e*\mathrm{\&beta;}---\left(7\right)$

$\begin{array}{c}{r}_1=\sqrt{{X_1}^2+{Y_1}^2+{Z_1}^2}\\ r_2=\sqrt{{X_2}^2+{Y_2}^2+{Z_2}^2}\\ \mathrm{\&beta;}=\arccos \left(\frac{(X_1{X}_2+Y_1{Y}_2+Z_1{Z}_2)}{r_1{r}_2}\right)\end{array}---\left(8\right)$

其中R_e为地球半径。

步骤209：通过弹道目标最大飞行高度H_max和射程D到弹道目标发射点模型偏差数据库中查询模型偏差ΔX；

步骤210：利用基于模型偏差修正发射点方法计算修正后的弹道目标发射点位置<L_launch,B_launch,H_launch>。

下面详细介绍上述步骤203和步骤206中基于弹道目标运动方程的不敏滤波方法，步骤204和步骤207中基于弹道目标运动方程的4阶龙格库塔轨道外推方法，步骤209中构建弹道目标发射点模型偏差数据库方法和步骤210中基于模型偏差修正发射点方法。

(一)基于弹道目标运动方程的不敏滤波方法：

假定目标的状态向量为弹道导弹的运动模型为

$\left[\begin{array}{c}\stackrel{\&CenterDot;}{x}\\ \stackrel{\&CenterDot;}{y}\\ \stackrel{\&CenterDot;}{z}\\ \stackrel{\&CenterDot;\&CenterDot;}{x}\\ \stackrel{\&CenterDot;\&CenterDot;}{y}\\ \stackrel{\&CenterDot;\&CenterDot;}{z}\end{array}\right]=\left[\begin{array}{c}\stackrel{\&CenterDot;}{x}\\ \stackrel{\&CenterDot;}{y}\\ \stackrel{\&CenterDot;}{z}\\ -\frac{\mathrm{\&rho;}\left(h\right)}{2\mathrm{\&beta;}}vx-\frac{\mathrm{\&mu;}G}{{\left(r\right)}^3}x+(2{\mathrm{\&omega;v}}_y+{\mathrm{\&omega;}}^{2}x)\\ -\frac{\mathrm{\&rho;}\left(h\right)}{2\mathrm{\&beta;}}vy-\frac{\mathrm{\&mu;}G}{{\left(r\right)}^3}y+(-2{\mathrm{\&omega;v}}_x+{\mathrm{\&omega;}}^{2}y)\\ -\frac{\mathrm{\&rho;}\left(h\right)}{2\mathrm{\&beta;}}{\mathrm{vv}}_z-\frac{\mathrm{\&mu;}G}{{\left(r\right)}^3}z\end{array}\right]---\left(9\right)$

其中ρ(h)＝ρ₀e^-kh为地球大气密度函数，其中ρ₀＝1.22kg/m³，k＝0.14141×10^-3m^-1，h为目标的海拔高度；β为弹道系数；μ_G＝3.986005×10¹⁴m³/s²为地球万有引力常量；ω＝7.292116×10^-5rad/s为地球自转速度。

雷达测量模型h为

$h\left(X\right)=\left[\begin{array}{c}\sqrt{x^2+y^2+z^2}\\ \arctan (x/y)\\ \arctan (z/\sqrt{x^2+y^2})\end{array}\right]---\left(10\right)$

假定雷达量测数据的噪声δ为不相关的零均值高斯白噪声，其中距离噪声的方差为σ_r，方位角噪声的方差为σ_α，仰角噪声的方差为σ_β。

目标的状态向量为β为弹道系数，则状态向量维数L＝7。不敏滤波的处理过程为：首先利用不敏变换生成采样点，采用常用的标准对称采样方法。在目标的运动模型(见式9)和雷达量测模型(见式10)是非线性函数，且目标在k-1时刻的估计均值为和P_k-1/k-1时，选择2L+1＝15个采样点按式11计算：

$\left\{\begin{array}{c}{\mathrm{\&chi;}}_0={\hat{X}}_{k-1/k-1},W_0^{\left(m\right)}=\frac{\mathrm{\&lambda;}}{2(L+\mathrm{\&lambda;})},W_0^{\left(c\right)}=\frac{\mathrm{\&lambda;}}{2(L+\mathrm{\&lambda;})}+1-{\mathrm{\&alpha;}}^2+\mathrm{\&gamma;};\\ {\mathrm{\&chi;}}_i={\hat{X}}_{k-/k-1}+{\left(\sqrt{(L+\mathrm{\&lambda;})P_{k-1/k-1}}\right)}_i,W_i^{\left(m\right)}=W_i^{\left(c\right)}=\frac{1}{2(L+\mathrm{\&lambda;})},i=1,\mathrm{...},L;\\ {\mathrm{\&chi;}}_{i+L}={\hat{X}}_{k-1/k-1}-{\left(\sqrt{(L+\mathrm{\&lambda;})P_{k-1/k-1}}\right)}_i,W_{i+L}^{\left(m\right)}=W_{i+L}^{\left(c\right)}=\frac{1}{2(L+\mathrm{\&lambda;})},i=1,\mathrm{...},L;\end{array}\right.---\left(11\right)$

其中λ＝α²(L+κ)为比例参数，作为控制采样点到均值的距离；一般取α＝0.5；κ＝3-L；γ＝2；式中为(L+λ)P_k-1/k-1均方根矩阵的第i行。

对时刻k状态的一部提前预测，计算目标状态预测：

${\hat{X}}_{k/k-1}=\underset{i=0}{\overset{2n}{\&Sigma;}}{W}_i^{\left(m\right)}{\mathrm{\&chi;}}_{i,k/k-1}---\left(12\right)$

χ_k/k-1＝f(χ_k/k-1,k-1)            (13)

计算目标预测协方差矩阵：

$P_{k|k-1}=\underset{i=0}{\overset{2n}{\&Sigma;}}(W_i^{\left(c\right)}\&lsqb;{\mathrm{\&chi;}}_{i,k/k-1}-{\hat{X}}_{k/k-1}\&rsqb;{\&lsqb;{\mathrm{\&chi;}}_{i,k/k-1}-{\hat{X}}_{k/k-1}\&rsqb;}^T)+Q_{k-1}---\left(14\right)$

计算时刻k量测的预测：

${\hat{Z}}_{k/k-1}=\underset{i=0}{\overset{2n}{\&Sigma;}}{W}_i{\mathrm{\&zeta;}}_{i,k/k-1}---\left(15\right)$

ζ_k/k-1＝h(χ_k/k-1,k)         (16)

$P_{{\tilde{z}}_k}=\underset{i=0}{\overset{2n}{\&Sigma;}}{W}_i^{\left(m\right)}({\mathrm{\&zeta;}}_{k/k-1}^i-{\hat{z}}_{k|k-1}){({\mathrm{\&zeta;}}_{k|k-1}^i-{\hat{z}}_{k|k-1})}^T+R_k---\left(17\right)$

$P_{{\tilde{x}}_k{\tilde{z}}_k}=\underset{i=0}{\overset{2n}{\&Sigma;}}{W}_i^{\left(m\right)}({\mathrm{\&zeta;}}_{k|k-1}^i-{\hat{x}}_{k|k-1}){({\mathrm{\&zeta;}}_{k/k-1}^i-{\hat{x}}_{k|k-1})}^T---\left(18\right)$

计算新息：

$v_k=Z_k-{\hat{Z}}_{k/k-1}---\left(19\right)$

计算增益矩阵：

$K_k=P_{{\tilde{x}}_k{\tilde{z}}_k}{P}_{{\tilde{z}}_k}^{-1}---\left(20\right)$

计算目标状态更新：

${\hat{X}}_{k|k}={\hat{X}}_{k|k-1}+K_k{v}_k---\left(21\right)$

计算目标状态误差协方差的更新：

${\hat{P}}_{k|k}={\hat{P}}_{k|k-1}-K_k{P}_{{\tilde{z}}_k}{K}_k^T---\left(22\right)$

(二)基于弹道目标运动方程的4阶龙格库塔轨道外推方法：

采用4阶龙格库塔方法进行外推，所用的计算增函数公式的形式为：

Δy_n＝y_n+1-y_n＝αk₁+βk₂+γk₃+δk₄     (23)

其中

$\left\{\begin{array}{c}{k}_1=y^{\&prime;}(t_n,y_n)h\\ k_2=y^{\&prime;}(t_n+h/2,y_n+k_1/2)h\\ k_3=y^{\&prime;}(t_n+h/2,y_n+k_2/2)h\\ k_4=y^{\&prime;}(t_n+h/2,y_n+k_3/2)h\end{array}\right.---\left(24\right)$

在满足各对应项系数相等的条件下选取未定系数α＝h/6，β＝h/3，γ＝h/3，δ＝h/6的适当值，y′(t_n,y_n)为弹道导弹运动方程见式9，然后将这些值代入上式可得

$y_{n+1}-y_n=\frac{h}{6}(k_1+2k_2+2k_3+k_4)---\left(25\right)$

只要所取间隔h大小合适，用龙格库塔法解微分方程可确保弹道外推具有足够的准确性，本实例中取h＝0.1秒。

外推结束条件是判别弹道目标的海拔高度是否小于等于地表高度，其中地表高度通过目标的经度和纬度在地理信息数据中查询获得。

(三)构建弹道目标发射点模型偏差数据库方法：

弹道导弹的主动段的运动模型复杂且不同弹种并不相同，传统的发射点估计中直接采用弹道目标的自由段模型进行外推，它与真实的发点位置存在一定偏差。

理论上在关机点前，应该使用主动段模型而不是自由段模型反推发点，因此产生了模型偏差的ΔX。下面分别给出基于仿真计算和实验数据构建模型偏差库的方法。

(1)基于仿真计算构建模型偏差库

依据射程和最大飞行高度选择典型的弹道导弹主动段模型进行仿真：

a)假设发点位置O₁，仿真计算出关机点的位置和速度[r,v]；

b)通过关机点的位置和速度[r,v]利用自由段模型反推得到发射点估计值O₂(步骤205计算的结果)；

c)依据发点O₁和O₂，计算模型偏差ΔX＝f_distance(O₁,O₂)，其中f_distance()为计算两点的距离函数，见式7。

在模型偏差数据库中记录该仿真的条件：主动段模型类型，射程最大飞行高度和计算的结果模型偏差ΔX。

(2)基于实验数据构建模型偏差库

通过已知真实的发射点位置O₁，步骤205计算的发射点估计O₂，计算模型偏差ΔX＝f_distance(O₁,O₂)，在模型偏差数据库中记录该仿真的条件：导弹型号，射程最大飞行高度和计算的结果模型偏差ΔX。

(四)基于模型偏差修正发射点方法：

如图4所示，目标是估计发点的真实位置需要对步骤205的发射点估计结果修正模型偏差ΔX。弹道目标的外推轨迹在地表的投影，即为发射平面与地球球面的交线，沿该投影线修正模型偏差ΔX可得到

${\hat{X}}_O={\hat{X}}_{O^{\&prime;}}-\mathrm{\&Delta;}X---\left(26\right)$

其中模型偏差ΔX依赖于弹道导弹的模型，通过步骤209查询获得。

综上所述，本发明的弹道目标发射点估计方法，该方法包括：

(1)对雷达测量的点迹数据进行预处理，进行反向跟踪滤波，采用4阶龙格库塔方法计算发射点位置；

(2)对雷达测量的点迹数据进行正向跟踪滤波，采用龙格库塔方法计算落点位置，并计算弹道目标的最大飞行高度和射程；

(3)利用模型偏差进行发射点位置的修正。

对雷达测量的点迹数据进行预处理的步骤包括：对雷达测量点迹数据进行野值剔除，并进行坐标转换。对雷达测量点迹数据进行野值剔除的步骤包括：对点迹数据的距离、方位、俯仰、径向速度进行逐项检测，采用3-Sigma的方法判别该点迹是否是野值，如是则删除该点迹数据。3-Sigma的野值判别方法是测量数据值与预测数据值的差值是否超过3倍的测量误差，如超过则判别为野值。坐标转换是将雷达测量点迹从雷达站心的球坐标系转换到地心地固(Earth-Centered Earth-Fixed，ECEF)坐标系。反向跟踪滤波的步骤包括：将雷达测量点迹按时间逆序后，进行迭代滤波，滤波方法采用基于弹道目标运动方程的不敏滤波(UKF)方法。

采用4阶龙格库塔方法计算发射点位置的步骤包括：通过反向跟踪滤波后得到弹道目标估计的位置和速度，采用基于弹道目标运动方程的4阶龙格库塔方法进行反向轨道迭代外推，外推结束条件是弹道目标的海拔高度小于等于地表高度，其中地表高度通过地理信息数据查询获得；通过坐标转换计算出发射点的位置(经度、纬度和高度)。

对雷达测量的点迹数据进行正向跟踪滤波的步骤包括：将雷达测量点迹按时间顺序，进行迭代滤波，滤波方法采用基于弹道目标运动模型的不敏滤波方法。

采用龙格库塔方法计算落点位置的步骤包括：通过正向滤波后得到弹道目标估计的位置和速度，采用基于弹道目标运动方程的4阶龙格库塔方法进行正向轨道迭代外推，外推结束条件是弹道目标的海拔高度小于等于地表高度，其中地表高度通过地理信息数据查询获得；通过坐标转换计算出落点的位置(经度、纬度和高度)。

计算弹道目标的最大飞行高度和射程的步骤包括：通过反向外推轨道和正向外推轨道可统计得出弹道目标的最大飞行高度，通过发射点估计位置和落点估计位置可计算出弹道目标的射程。

利用模型偏差进行发射点位置的修正的步骤包括：通过弹道目标最大飞行高度和射程可在弹道目标发射点模型偏差数据库中查询到模型偏差；将通过外推的弹道轨道和估计得发射点位置进行模型偏差修正，可最终计算出修正后的发射点估计位置。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于模型修正的弹道目标发射点估计装置，其包括点迹数据预处理模块，所述点迹数据预处理模块用于对点迹数据进行预处理，雷达量测点迹数据为Z＝{z₁,...,z_k}Z＝{z₁,...,z_k}；

其特征在于：所述装置还包括：

时间逆序处理模块，其用于对雷达量测点迹数据Z＝{z₁,...,z_k}进行时间逆序处理，转换为Z'＝{z_k,z_k-1,...,z₁}，并去除时间重复的点迹数据；

反向滤波处理模块，其用于对雷达量测点迹的时间逆序数据Z'＝{z_k,z_k-1,...,z₁}，利用基于弹道目标运动方程进行反向滤波处理；

反向轨道外推计算模块，其用于利用反向滤波的估计值，采用基于弹道目标运动方程的龙格库塔4阶方法进行反向轨道外推计算，得到反向轨道外推结果X_launch＝<X₁,Y₁,Z₁>_ecef；

发射点位置计算模块，其用于计算弹道目标发射点位置，利用X_launch＝<X₁,Y₁,Z₁>_ecef从地心地固坐标系下转换到WGS84地心坐标系<L₁,B₁,H₁>；

正向滤波处理模块，其用于对雷达量测点迹数据Z＝{z₁,...,z_k}，利用基于弹道目标运动方程进行正向滤波处理；

正向轨道外推计算模块，其用于利用正向滤波的估计值，采用基于弹道目标运动方程的4阶龙格库塔方法进行正向轨道外推，得到正向轨道外推结果X_impact＝<X₂,Y₂,Z₂>_ecef；

最大飞行高度和射程计算模块，其用于计算弹道目标最大飞行高度和射程，通过X_launch＝<X₁,Y₁,Z₁>_ecef和X_impact＝<X₂,Y₂,Z₂>_ecef统计得出弹道目标的最大飞行高度H_max，计算弹道目标落点估计位置，将从地心地固坐标系下转换到WGS84地心坐标系<L₂,B₂,H₂>，通过发射点估计位置X_launch＝<X₁,Y₁,Z₁>_ecef和落点估计位置X_impact＝<X₂,Y₂,Z₂>_ecef，计算出弹道目标的射程D；

模型偏差查询模块，其用于通过弹道目标最大飞行高度H_max和射程D到弹道目标发射点模型偏差数据库中查询模型偏差ΔX；

模型偏差修正模块，其用于利用基于模型偏差修正发射点方法计算修正后的弹道目标发射点位置<L_launch,B_launch,H_launch>。

2.一种基于模型修正的弹道目标发射点估计方法，其包括以下步骤：

步骤201，雷达量测点迹数据为Z＝{z₁,...,z_k}，对点迹数据进行预处理；

其特征在于：所述估计方法还包括以下步骤：

步骤202，对雷达量测点迹数据Z＝{z₁,...,z_k}进行时间逆序处理，转换为Z'＝{z_k,z_k-1,...,z₁}，并去除时间重复的点迹数据；

步骤203，对雷达量测点迹的时间逆序数据Z'＝{z_k,z_k-1,...,z₁}，利用基于弹道目标运动方程进行反向滤波处理；

步骤204，利用反向滤波的估计值，采用基于弹道目标运动方程的龙格库塔4阶方法进行反向轨道外推计算，得到反向轨道外推结果X_launch＝<X₁,Y₁,Z₁>_ecef；

步骤205，计算弹道目标发射点位置，利用X_launch＝<X₁,Y₁,Z₁>_ecef从地心地固坐标系下转换到WGS84地心坐标系<L₁,B₁,H₁>；

步骤206，对雷达量测点迹数据Z＝{z₁,...,z_k}，利用基于弹道目标运动方程进行正向滤波处理；

步骤207，利用正向滤波的估计值，采用基于弹道目标运动方程的4阶龙格库塔方法进行正向轨道外推，得到正向轨道外推结果X_impact＝<X₂,Y₂,Z₂>_ecef；

步骤208，计算弹道目标最大飞行高度和射程，通过X_launch＝<X₁,Y₁,Z₁>_ecef和X_impact＝<X₂,Y₂,Z₂>_ecef统计得出弹道目标的最大飞行高度H_max，计算弹道目标落点估计位置，将从地心地固坐标系下转换到WGS84地心坐标系<L₂,B₂,H₂>，通过发射点估计位置X_launch＝<X₁,Y₁,Z₁>_ecef和落点估计位置X_impact＝<X₂,Y₂,Z₂>_ecef，计算出弹道目标的射程D；

步骤209，通过弹道目标最大飞行高度H_max和射程D到弹道目标发射点模型偏差数据库中查询模型偏差ΔX；

步骤210，利用基于模型偏差修正发射点方法计算修正后的弹道目标发射点位置<L_launch,B_launch,H_launch>。

3.如权利要求2所述的基于模型修正的弹道目标发射点估计方法，其特征在于：反向滤波处理与正向滤波处理均采用不敏滤波算法。

4.如权利要求3所述的基于模型修正的弹道目标发射点估计方法，其特征在于：步骤203与步骤206中基于弹道目标运动方程的不敏滤波方法包括以下步骤：

一、假定目标的状态向量为弹道导弹的运动模型为：

$\left[\begin{array}{c}\stackrel{\&CenterDot;}{x}\\ \stackrel{\&CenterDot;}{y}\\ \stackrel{\&CenterDot;}{z}\\ \stackrel{\&CenterDot;\&CenterDot;}{x}\\ \stackrel{\&CenterDot;\&CenterDot;}{y}\\ \stackrel{\&CenterDot;\&CenterDot;}{z}\end{array}\right]=\left[\begin{array}{c}\stackrel{\&CenterDot;}{x}\\ \stackrel{\&CenterDot;}{y}\\ \stackrel{\&CenterDot;}{z}\\ -\frac{\mathrm{\&rho;}\left(h\right)}{2\mathrm{\&beta;}}vx-\frac{\mathrm{\&mu;}G}{{\left(r\right)}^3}x+(2{\mathrm{\&omega;v}}_y+{\mathrm{\&omega;}}^{2}x)\\ -\frac{\mathrm{\&rho;}\left(h\right)}{2\mathrm{\&beta;}}vy-\frac{\mathrm{\&mu;}G}{{\left(r\right)}^3}y+(-2{\mathrm{\&omega;v}}_x+{\mathrm{\&omega;}}^{2}y)\\ -\frac{\mathrm{\&rho;}\left(h\right)}{2\mathrm{\&beta;}}{\mathrm{vv}}_z-\frac{\mathrm{\&mu;}G}{{\left(r\right)}^3}z\end{array}\right]$    公式9，

其中ρ(h)＝ρ₀e^-kh为地球大气密度函数，其中ρ₀＝1.22kg/m³，k＝0.14141×10^-3m^-1，h为目标的海拔高度；β为弹道系数；μ_G＝3.986005×10¹⁴m³/s²为地球万有引力常量；ω＝7.292116×10^-5rad/s为地球自转速度；

二、雷达测量模型h为：

$h\left(X\right)=\left[\begin{array}{c}\sqrt{x^2+y^2+z^2}\\ \arctan (x/y)\\ \arctan (z/\sqrt{x^2+y^2})\end{array}\right]$    公式10，

假定雷达量测数据的噪声δ为不相关的零均值高斯白噪声，其中距离噪声的方差为σ_r，方位角噪声的方差为σ_α，仰角噪声的方差为σ_β；

三、目标的状态向量为β为弹道系数，则状态向量维数L＝7，不敏滤波的处理过程为：首先利用不敏变换生成采样点，采用常用的标准对称采样方法，在目标的运动模型和雷达量测模型是非线性函数，且目标在k-1时刻的估计均值为和P_k-1/k-1时，选择2L+1＝15个采样点按式11计算：

$\left\{\begin{array}{c}{\mathrm{\&chi;}}_0={\hat{X}}_{k-1/k-1},W_0^{\left(m\right)}=\frac{\mathrm{\&lambda;}}{2(L+\mathrm{\&lambda;})},W_0^{\left(c\right)}=\frac{\mathrm{\&lambda;}}{2(L+\mathrm{\&lambda;})}+1-{\mathrm{\&alpha;}}^2+\mathrm{\&gamma;};\\ {\mathrm{\&chi;}}_i={\hat{X}}_{k-/k-1}+{\left(\sqrt{(L+\mathrm{\&lambda;})P_{k-1/k-1}}\right)}_i,W_i^{\left(m\right)}=W_i^{\left(c\right)}=\frac{1}{(L+\mathrm{\&lambda;})},i=1,\mathrm{...},L;\\ {\mathrm{\&chi;}}_{i+L}={\hat{X}}_{k-1/k-1}-{\left(\sqrt{(L+\mathrm{\&lambda;})P_{k-1/k-1}}\right)}_i,W_{i+L}^{\left(m\right)}=W_{i+L}^{\left(c\right)}=\frac{1}{2(L+\mathrm{\&lambda;})},i=1,\mathrm{...},L;\end{array}\right.$    公式11，

其中λ＝α²(L+κ)为比例参数，作为控制采样点到均值的距离；一般取α＝0.5；κ＝3-L；γ＝2；式中为(L+λ)P_k-1/k-1均方根矩阵的第i行；

对时刻k状态的一部提前预测，计算目标状态预测：

   公式12，χ_k/k-1＝f(χ_k/k-1,k-1)   公式13；

计算目标预测协方差矩阵：

$P_{k|k-1}=\underset{i=0}{\overset{2n}{\&Sigma;}}(W_i^{\left(c\right)}\&lsqb;{\mathrm{\&chi;}}_{i,k/k-1}-{\hat{X}}_{k/k-1}\&rsqb;{\&lsqb;{\mathrm{\&chi;}}_{i,k/k-1}-{\hat{X}}_{k/k-1}\&rsqb;}^T)+Q_{k-1}$    公式14；

计算时刻k量测的预测：

${\hat{Z}}_{k/k-1}=\underset{i=0}{\overset{2n}{\&Sigma;}}{W}_i{\mathrm{\&zeta;}}_{i,k/k-1}$    公式15，

ζ_k/k-1＝h(χ_k/k-1,k)   公式16，

$P_{{\tilde{z}}_k}=\underset{i=0}{\overset{2n}{\&Sigma;}}{W}_i^{\left(m\right)}({\mathrm{\&zeta;}}_{k/k-1}^i-{\hat{z}}_{k|k-1}){({\mathrm{\&zeta;}}_{k|k-1}^i-{\hat{z}}_{k|k-1})}^T+R_k$    公式17，

$P_{{\tilde{x}}_k{\tilde{z}}_k}=\underset{i=0}{\overset{2n}{\&Sigma;}}{W}_i^{\left(m\right)}({\mathrm{\&zeta;}}_{k|k-1}^i-{\hat{x}}_{k|k-1}){({\mathrm{\&zeta;}}_{k/k-1}^i-{\hat{x}}_{k|k-1})}^T$    公式18；

计算新息：

$v_k=Z_k-{\hat{Z}}_{k/k-1}$    公式19；

计算增益矩阵：

$K_k=P_{{\tilde{x}}_k{\tilde{z}}_k}{P}_{{\tilde{z}}_k}^{-1}$    公式20；

计算目标状态更新：

${\hat{X}}_{k|k}={\hat{X}}_{k|k-1}+K_k{v}_k$    公式21；

计算目标状态误差协方差的更新：

${\hat{P}}_{k|k}={\hat{P}}_{k|k-1}-K_k{P}_{{\tilde{z}}_k}{K}_k^T$    公式22。

5.如权利要求2所述的基于模型修正的弹道目标发射点估计方法，其特征在于：X_launch＝<X₁,Y₁,Z₁>_ecef从地心地固坐标系下转换到WGS84地心坐标系<L₁,B₁,H₁>、以及X_impact＝<X₂,Y₂,Z₂>_ecef从地心地固坐标系下转换到WGS84地心坐标系<L₁,B₁,H₁>均通过公式6得到：

$r=\sqrt{x^2+y^2}$

E²＝a²-b²

F＝54b²z²

G＝r²+(1-e²)z²-e²E²

$C=\frac{e^4{\mathrm{Fr}}^2}{G^3}$

$S=\sqrt[3]{1+C+\sqrt{C^2+2C}}$

$P=\frac{F}{3{(S+\frac{1}{S}+1)}^2{G}^2}$

$Q=\sqrt{1+2e^{4}P}$

$r_0=\frac{-\left({\mathrm{Pe}}^{2}r\right)}{1+Q}+\sqrt{\frac{1}{2}{a}^2(1+\frac{1}{Q})-\frac{P(1-e^2)Z^2}{Q(1+Q)}-\frac{1}{2}{\mathrm{Pr}}^2}$

$U=\sqrt{{(r-e^2{r}_0)}^2+Z^2}$

$V=\sqrt{{(r-e^2{r}_0)}^2+(1-e^2)Z^2}$

$Z_0=\frac{b^{2}Z}{aV}$

                                            公式(6)

$H_1=U(1-\frac{b^2}{aV})$

$B_1=arctan\&lsqb;\frac{Z+e^2{Z}_0}{r}\&rsqb;$

L₁＝arctan 2(Y,X)

其中，x，y，z分别表示地心地固坐标系中的三个坐标轴；

L，B，H分别表示雷达站心的经度、纬度和高度；

a＝6378137m为地球长半轴，b＝6356752.3142m为地球短半轴，e＝0.006694380为地球椭球第一偏心率；

E，F，G，C，S，P，Q，r₀，U，V，Z₀分别为临时计算参数。

6.如权利要求2所述的基于模型修正的弹道目标发射点估计方法，其特征在于：所述最大飞行高度和射程计算模块通过式7计算出弹道目标的射程D：

$f_{dis\tan ce}(\left[\begin{array}{c}{X}_1\\ Y_1\\ Z_1\end{array}\right],\left[\begin{array}{c}{X}_2\\ Y_2\\ Z_2\end{array}\right])=R_e*\mathrm{\&beta;}$

$r_1=\sqrt{{X_1}^2+{Y_2}^2+{Z_1}^2}$    公式7，

$r_2=\sqrt{{X_2}^2+{Y_2}^2+{Z_2}^2}$

$\mathrm{\&beta;}=arccos\left(\frac{(X_1{X}_2+Y_1{Y}_2+Z_1{Z}_2)}{r_1{r}_2}\right)$

其中，<X₁,Y₁,Z₁>_ecef为地心地固坐标系下弹道目标发射点的位置参数，<X₂,Y₂,Z₂>_ecef为地心地固坐标系下弹道目标落地的位置参数，R_e为地球半径，β为方位，r₁，r₂分别为临时计算参数。

7.如权利要求2所述的基于模型修正的弹道目标发射点估计方法，其特征在于：步骤204与步骤207中基于弹道目标运动方程的4阶龙格库塔轨道外推方法包括以下步骤：

采用4阶龙格库塔方法进行外推，所用的计算增函数公式的形式为：

Δy_n＝y_n+1-y_n＝αk₁+βk₂+γk₃+δk₄   公式23，

其中

$\left\{\begin{array}{c}{k}_1=y^{\&prime;}(t_n,y_n)h\\ k_2=y^{\&prime;}(t_n+h/2,y_n+k_1/2)h\\ k_3=y^{\&prime;}(t_n+h/2,y_n+k_2/2)h\\ k_4=y^{\&prime;}(t_n+h/2,y_n+k_3/2)h\end{array}\right.$    公式24；

在满足各对应项系数相等的条件下选取未定系数α＝h/6，β＝h/3，γ＝h/3，δ＝h/6的适当值，y′(t_n,y_n)为弹道导弹运动方程见公式9，然后将这些值代入上式可得

$y_{n+1}-y_n=\frac{h}{6}(k_1+2k_2+2k_3+k_4)$    公式25。

8.如权利要求6所述的基于模型修正的弹道目标发射点估计方法，其特征在于：步骤209中构建弹道目标发射点模型偏差数据库方法采用基于仿真计算构建模型偏差库，依据射程和最大飞行高度选择典型的弹道导弹主动段模型进行仿真：a)假设发点位置O₁，仿真计算出关机点的位置和速度[r,v]；b)通过关机点的位置和速度[r,v]利用自由段模型反推得到步骤205计算的结果即发射点估计值O₂；c)依据发点O₁和O₂，计算模型偏差ΔX＝f_distance(O₁,O₂)，其中f_distance()为计算两点的距离函数为公式7；在模型偏差数据库中记录该仿真的条件：主动段模型类型，射程最大飞行高度和计算的结果模型偏差ΔX。

9.如权利要求2所述的基于模型修正的弹道目标发射点估计方法，其特征在于：步骤209中构建弹道目标发射点模型偏差数据库方法采用基于实验数据构建模型偏差库：通过已知真实的发射点位置O₁，计算的发射点估计O₂，计算模型偏差ΔX＝f_distance(O₁,O₂)，在模型偏差数据库中记录该仿真的条件：导弹型号，射程最大飞行高度和计算的结果模型偏差ΔX。

10.如权利要求2所述的基于模型修正的弹道目标发射点估计方法，其特征在于：基于模型偏差修正发射点方法包括以下步骤：

目标是估计发点的真实位置需要对发射点估计结果修正模型偏差ΔX，弹道目标的外推轨迹在地表的投影，即为发射平面与地球球面的交线，沿该投影线修正模型偏差ΔX可得到

${\hat{X}}_O={\hat{X}}_{O^{\&prime;}}-\mathrm{\&Delta;}X$    公式26，

其中模型偏差ΔX依赖于弹道导弹的模型，通过弹道目标发射点模型偏差数据库中查询获得。