CN111412919A

CN111412919A - 一种空间目标的初始轨道误差计算方法及装置

Info

Publication number: CN111412919A
Application number: CN202010263472.1A
Authority: CN
Inventors: 李大卫; 杨旭; 程昊文; 金勇�
Original assignee: National Astronomical Observatories of CAS
Current assignee: National Astronomical Observatories of CAS
Priority date: 2020-04-07
Filing date: 2020-04-07
Publication date: 2020-07-14
Anticipated expiration: 2040-04-07
Also published as: CN111412919B

Abstract

本公开提供了一种空间目标的初始轨道误差计算方法及装置，包括：基于观测数据，计算空间目标在J2000坐标系下的初始轨道以及初始轨道协方差；根据初始轨道与初始轨道协方差，计算空间目标在J2000坐标系下指定时刻的预报位置、预报速度以及在J2000坐标系下指定时刻的初始轨道协方差；根据J2000坐标系下指定时刻的初始轨道协方差，确定在测量坐标系下指定时刻的初始轨道协方差，以得到指定时刻赤经与赤纬理论预报值的误差，测量坐标系以观测站为原点；根据J2000坐标系下指定时刻的初始轨道协方差，确定在UNW坐标系下指定时刻的初始轨道协方差，以得到指定时刻在U轴、N轴及W轴方向的预报误差。本公开可以提高初始轨道误差估计的准确度。

Description

一种空间目标的初始轨道误差计算方法及装置

技术领域

本公开涉及空间碎片监测领域，特别涉及一种空间目标的初始轨道误差计算方法及装置。

背景技术

通常情况下，对于新发现的空间碎片，需要先确定该空间碎片的初始轨道，再利用初始轨道引导其他观测设备开展后续观测，因此，需要在初始确定时提供轨道的误差和误差随时间的发散量，从而为后续的观测设备给出准确的搜索范围。

目前，地基监测对初始轨道误差的计算基于广义拉普拉斯法的最小二乘估计，但是该方法对监测数据误差引起的误差模型的估计不完整，误差的估计偏小。因此，亟需一种空间目标的初始轨道误差计算方法，以提高误差估计的准确度。

发明内容

本公开的目的是提供一种空间目标的初始轨道误差计算方法及装置，以提高初始轨道误差估计的准确度。

为实现上述目的，本公开实施例提供一种空间目标的初始轨道误差计算方法，所述方法包括：

基于观测数据，计算空间目标在J2000坐标系下的初始轨道以及初始轨道协方差；

根据所述初始轨道与初始轨道协方差，计算所述空间目标在所述J2000坐标系下指定时刻的预报位置、预报速度以及在所述J2000坐标系下指定时刻的初始轨道协方差；

根据所述J2000坐标系下指定时刻的初始轨道协方差，确定在测量坐标系下指定时刻的初始轨道协方差，以得到所述指定时刻赤经与赤纬理论预报值的误差，所述测量坐标系以观测站为原点；

根据所述J2000坐标系下指定时刻的初始轨道协方差，确定在UNW坐标系下指定时刻的初始轨道协方差，以得到所述指定时刻在U轴、N轴及W轴方向的预报误差。

进一步的，本公开实施例还提供一种空间目标的初始轨道误差计算装置，所述装置包括：

第一计算模块，用于基于观测数据，计算空间目标在J2000坐标系下的初始轨道以及初始轨道协方差；

第二计算模块，用于根据所述初始轨道与初始轨道协方差，计算所述空间目标在所述J2000坐标系下指定时刻的预报位置、预报速度以及在所述J2000坐标系下指定时刻的初始轨道协方差；

第三计算模块，用于根据所述J2000坐标系下指定时刻的初始轨道协方差，确定在测量坐标系下指定时刻的初始轨道协方差，以得到所述指定时刻赤经与赤纬理论预报值的误差，所述测量坐标系以观测站为原点；

第四计算模块，用于根据所述J2000坐标系下指定时刻的初始轨道协方差，确定在UNW坐标系下指定时刻的初始轨道协方差，以得到所述指定时刻在U轴、N轴及W轴方向的预报误差。

本公开提供的空间目标的初始轨道误差计算方法，可以利用光学望远镜的观测数据、观测数据的测量误差，计算初始轨道和轨道误差，基于简化的轨道动力学模型，预报任意时刻空间目标的位置和位置误差，根据位置误差制定后续观测站观测该空间目标的观测策略。

附图说明

图1是本公开实施例提供的一种空间目标的初始轨道误差计算方法流程图；

图2是赤经误差、赤纬误差、斜距误差随时间变化的示意图；

图3是UNW坐标系中U轴、N轴及W轴方向的误差随时间变化的示意图。

具体实施方式

本公开实施例提供一种空间目标的初始轨道误差计算方法及装置。

为了使本技术领域的人员更好地理解本公开中的技术方案，下面将结合本公开实施方式中的附图，对本公开实施方式中的技术方案进行清楚、完整地描述，显然，所描述的实施方式仅仅是本公开一部分实施方式，而不是全部的实施方式。基于本公开中的实施方式，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施方式，都应当属于本公开保护的范围。

结合图1，本公开实施方式中提供的一种空间目标的初始轨道误差计算方法，可以包括：

S1：基于观测数据，计算空间目标在J2000坐标系下的初始轨道以及初始轨道协方差。

S2：根据所述初始轨道与初始轨道协方差，计算所述空间目标在所述J2000坐标系下指定时刻的预报位置、预报速度以及在所述J2000坐标系下指定时刻的初始轨道协方差。

S3：根据所述J2000坐标系下指定时刻的初始轨道协方差，确定在测量坐标系下指定时刻的初始轨道协方差，以得到所述指定时刻赤经与赤纬理论预报值的误差，所述测量坐标系以观测站为原点。

S4：根据所述J2000坐标系下指定时刻的初始轨道协方差，确定在UNW坐标系下指定时刻的初始轨道协方差，以得到所述指定时刻在U轴、N轴及W轴方向的预报误差。

在本实施方式中，步骤S1可以包括如下步骤S11～S13，具体的：

S11：在J2000坐标系下，基于广义拉普拉斯初始轨道确定算法，建立所述观测数据与所述初始轨道之间的线性模型，以计算得到所述初始轨道，所述线性模型包含误差项。

S12：计算所述误差项的协方差。

S13：根据最小二乘理论以及所述误差项的协方差，计算得到所述初始轨道协方差。

在本实施方式中，步骤S11还可以包括如下步骤S111～S112，具体的：

S111：基于广义拉普拉斯初始轨道确定算法，建立所述观测数据与所述初始轨道的线性方程。

S112：根据所述线性方程与所述误差项，建立所述观测数据与所述初始轨道之间的线性模型。

具体的，线性方程如下：

其中，

λ＝cosδcosα

μ＝cosδsinα

ν＝sinδ

α为赤经；δ为赤纬；λ、μ、ν是关于赤经α和赤纬δ的函数；x₀、y₀、z₀为初始时刻的位置；

为初始时刻的速度；F、F_z、G、G_z是关于时间的函数；X_e、Y_e、Z_e为所述观测站的坐标。

具体的，线性模型如下：

y＝Xβ+ε

ε＝ε_y+ε_X

ξ＝[Δα Δδ]^T～σ²I

其中，

y＝[y₁ y₂ y₃]^T

y₁＝(νX_e-λZ_e)

y₂＝(νY_e-μZ_e)

y₃＝(μX_e-λY_e)

α为赤经；δ为赤纬；β为初始轨道；y为线性模型的观测量；X为设计矩阵；ε为误差项；ε_X为来自设计矩阵X的误差；ε_y为来自观测量y的误差；λ、μ、ν是关于赤经α和赤纬δ的函数；x₀、y₀、z₀为初始时刻的位置；

为初始时刻的速度；F、F_z、G、G_z是关于时间的函数；X_e、Y_e、Z_e为所述观测站的坐标；I为单位矩阵；ξ为观测数据的测量误差构成的矩阵，包含赤经误差和赤纬误差；σ为观测数据的误差值。

具体的，误差项的协方差通过如下方法计算：

ε＝ε_y+ε_X＝(U-A)ξ

ξ＝[Δα Δδ]^T

cov(ξ)＝σ²I

ε_yU＝[U₁ U₂]＝M_LXM_ZL

z＝(α δ)^T

L＝[λμν]^T

计算得到：

其中，

为初始时刻的速度；F、F_z、G、G_z是关于时间的函数；Σ为误差项的协方差；z为观测数据；σ为观测数据的误差值；I为单位矩阵。

具体的，初始轨道协方差通过如下方式计算得到：

计算所述初始轨道的最小二乘估计：

β＝(X^TΣ^-X)^-1Σ^-y

计算所述初始轨道协方差：

P＝(X^TWX)^-1＝(X^TΣ^-X)^-1

β为初始轨道；X为设计矩阵；Σ为误差项的协方差；P为在J2000坐标系下的初始轨道协方差；W为权矩阵；y为线性模型的观测量。

具体的，J2000坐标系下指定时刻的初始轨道协方差通过如下方法计算得到：

其中，

为指定时刻的预报位置；

为指定时刻的预报速度；Φ为状态转移矩阵；P为在J2000坐标系下的初始轨道协方差；P_t为在J2000坐标系下指定时刻的初始轨道协方差；t为任意指定历元；t₀为初始轨道的历元；

为t时刻的轨道参数，包括t时刻的位置和速度矢量；

为t₀时刻的轨道参数，包括t₀时刻的位置和速度矢量。

具体的，指定时刻赤经和赤纬理论预报值的误差通过如下方法计算得到：

P_t(α,δ)＝MPM^T

计算得到：

[σ(α),σ(δ)]＝diag(P_t(α,δ))

H为观测数据

与轨道参数

间的测量矩阵；Φ为状态转移矩阵；P_t(α,δ)为在测量坐标系下t时刻的初始轨道协方差；σ(α),σ(δ)分别为在赤经和赤纬的理论预报值的误差；

为t时刻第k组观测数据；P为在J2000坐标系下的初始轨道协方差；α_kt为t时刻的赤经；δ_kt为t时刻的赤纬；

为n组观测数据的全体；

为t₀时刻的轨道参数；观测数据共有n组，k表示第k组观测数据。

具体的，指定时刻在U轴、N轴及W轴方向的预报误差通过如下方法计算得到：

P_t＝ΦP₀Φ′

P_UNW＝MP_tM′

计算得到：

σ_UNW＝diag(P_UNW)

r为位置矢量；v为速度矢量；M_r为J2000坐标系下的位置转到UNW坐标系的转换矩阵；M_v为J2000坐标系下的速度转到UNW坐标系的转换矩阵；Φ为J2000坐标系下零时刻的位置速度到t时刻位置速度的状态转移矩阵；P_t为J2000坐标系下指定时刻的初始轨道协方差；P₀为J2000坐标系下初始轨道历元时刻的初始轨道协方差；P_UNW为在UNW坐标系下任意指定时刻的初始轨道协方差；σ_UNW为指定时刻在U轴、N轴及W轴方向的预报误差。

在一些实施方式中，所述方法还包括：

S5：根据所述UNW坐标系下指定时刻的初始轨道协方差以及所述空间目标在指定时刻的预报位置、预报速度，确定所述空间目标穿越误差椭球的时间差。

S6：计算在所述时间差内，观测设备的指向。

具体的，空间目标穿越误差椭球的时间差通过如下方法计算得到：

其中，Δt为时间差；

为所述指定时刻的预报速度；σ(r_u)为U轴方向的预报位置误差。

具体的，所述观测设备的指向通过如下方法计算得到：

其中，

α为赤经；δ为赤纬；Δt为时间差；x'、y'、z'为空间目标相对观测站的位置；R_e为观测站位置的地心距；w为地球自转角速度；m₀为观测站在初始时刻的赤经；α₀为空间目标在预报时刻赤经的理论观测值；

为空间目标的相对位置在赤道平面的投影；

为t时刻空间目标的相对位置。

在一些实施方式中，所述方法还包括：

S7：建立观测坐标系，所述观测坐标系为以观测设备为原点的直角坐标系，所述观测坐标系的z轴沿观测方向，x轴指北，y轴与所述x轴、z轴构成右手系。

S8：计算在所述观测坐标系下指定时刻的初始轨道协方差。

S9：根据所述观测坐标系下指定时刻的初始轨道协方差，确定垂直于所述观测方向的长轴的值。

S10：根据所述长轴的值，确定所述观测设备的视场。

具体的，观测坐标系下指定时刻的初始轨道协方差通过如下方法计算得到：

P′_t＝R_{J2000→观测坐标系}P_tR′_{J2000→观测坐标系}

其中，P′_t为观测坐标系下指定时刻的初始轨道协方差；P_t为J2000坐标系下指定时刻的初始轨道协方差；R_{J2000→观测坐标系}为从J2000坐标系向观测坐标系的旋转矩阵；R′_{J2000→观测坐标系}为R_{J2000→观测坐标系}的转置矩阵。

具体的，观测设备的视场通过如下方法计算得到：

其中，Fov为视场；2σ_l为长轴；

为空间目标相对于观测站的位置。

此外，本公开实施方式还提供一种空间目标的初始轨道误差计算装置，所述装置包括：

第一计算模块100，用于基于观测数据，计算空间目标在J2000坐标系下的初始轨道以及初始轨道协方差。

第二计算模块200，用于根据所述初始轨道与初始轨道协方差，计算所述空间目标在所述J2000坐标系下指定时刻的预报位置、预报速度以及在所述J2000坐标系下指定时刻的初始轨道协方差。

第三计算模块300，用于根据所述J2000坐标系下指定时刻的初始轨道协方差，确定在测量坐标系下指定时刻的初始轨道协方差，以得到所述指定时刻赤经与赤纬理论预报值的误差，所述测量坐标系以观测站为原点。

第四计算模块400，用于根据所述J2000坐标系下指定时刻的初始轨道协方差，确定在UNW坐标系下指定时刻的初始轨道协方差，以得到所述指定时刻在U轴、N轴及W轴方向的预报误差。

下面以一个具体的实施例，来进一步地说明本公开提供的空间目标的初始轨道误差计算方法。

首先，选取低轨空间碎片轨道，观测站的坐标为(40°N 120°E 1.0km)，仿真光学观测数据，弧段长度约6分钟，测角误差约5角秒，随机分布，采用本公开提供的方法得到初始轨道协方差的估计值。

(1)计算初始轨道和初始轨道协方差。如下面的表1为初始轨道(位置、速度)协方差；表2为初始轨道(开普勒根数)协方差。

X(re)

Y(re)

z(re)

Vx(re/tu)

Vy(re/tu)

Vz(re/tu)

X(re)

1.6330E-12

7.6088E-12

-7.9701E-12

-4.6360E-11

7.8632E-12

5.3165E-12

Y(re)

7.6088E-12

3.8199E-11

-3.9864E-11

-2.3178E-10

4.2646E-11

2.3272E-11

Z(re)

-7.9701E-12

-3.9864E-11

4.1797E-11

2.4231E-10

-4.4426E-11

-2.4714E-11

Vx(re/tu)

-4.6360E-11

-2.3178E-10

2.4231E-10

1.4136E-09

-2.5879E-10

-1.4223E-10

Vy(re/tu)

7.8632E-12

4.2646E-11

-4.4426E-11

-2.5879E-10

5.8752E-11

1.8645E-11

Vz(re/tu)

5.3165E-12

2.3272E-11

-2.4714E-11

-1.4223E-10

1.8645E-11

2.8908E-11

表1

a(km)

e

i(rad)

Omg(rad)

w(rad)

M(rad)

a(km)

3.1073E-01

1.6951E-05

-4.4663E-06

2.1932E-06

-3.0306E-02

3.0305E-02

E

1.6951E-05

9.3031E-10

-2.4290E-10

1.2661E-10

-1.6516E-06

1.6516E-06

i(rad)

-4.4663E-06

-2.4290E-10

6.5396E-11

-2.2670E-11

4.3579E-07

-4.3578E-07

Omg(rad)

2.1932E-06

1.2661E-10

-2.2670E-11

1.0386E-10

-2.1204E-07

2.1197E-07

w(rad)

-3.0306E-02

-1.6516E-06

4.3579E-07

-2.1204E-07

2.9563E-03

-2.9562E-03

M(rad)

3.0305E-02

1.6516E-06

-4.3578E-07

2.1197E-07

-2.9562E-03

2.9562E-03

表2

(2)计算初始轨道时刻后200秒时初始轨道和初始轨道协方差。如下面的表3为预报时间为200秒时的初始轨道(位置、速度)；表4为预报时间为200秒时的初始轨道(位置、速度)协方差。

	X(km)	Y(km)	z(km)	Vx(km/s)	Vy(km/s)	Vz(km/s)
							X(re)	-43.06592595	5862.678613	3178.552604	-7.733816259	-0.0485555	-0.033736865

表1

X(re)

Y(re)

z(re)

Vx(re/tu)

Vy(re/tu)

Vz(re/tu)

X(re)

6.2185E-11

-6.2281E-11

4.1392E-11

2.8884E-10

-5.7528E-11

-4.7554E-11

Y(re)

-6.2281E-11

6.3217E-11

-4.1894E-11

-2.8911E-10

6.0102E-11

4.6543E-11

Z(re)

4.1392E-11

-4.1894E-11

2.8390E-11

1.9195E-10

-3.9638E-11

-2.8885E-11

Vx(re/tu)

2.8884E-10

-2.8911E-10

1.9195E-10

1.3442E-09

-2.6474E-10

-2.2431E-10

Vy(re/tu)

-5.7528E-11

6.0102E-11

-3.9638E-11

-2.6474E-10

6.3921E-11

3.7459E-11

Vz(re/tu)

-4.7554E-11

4.6543E-11

-2.8885E-11

-2.2431E-10

3.7459E-11

5.1344E-11

表4

(3)计算赤经、赤纬的理论预报值的误差。

当预报时间为200秒时，赤经的误差为0.000170118度，赤纬的误差为0.000160759度，斜距的误差为0.178625549千米。每隔200秒，预报24小时，得到赤经误差、赤纬误差、斜距误差、UNW误差随时间变化图，即图2所示。

(4)计算UNW坐标系下U轴、N轴及W轴三个方向的预报误差。

当预报时间为200秒时，UNW三个方向的预报误差分别为0.05km、0.03km、0.05km。每隔200秒，预报24小时，得到UNW误差随时间变化图，即图3所示。

(5)计算观测设备的等待时间和指向

当预报时间为15400秒时，指向赤经为0.2756弧度，赤纬为-1.0497弧度，探测距离在赤道面的投影为1682km，观测站的赤经为0.3090弧度。当等待时间为提前和延后7秒，则赤纬保持不变，赤经的变化率为-2.1170e-04弧段/秒。

本说明书中的上述各个实施方式均采用递进的方式描述，各个实施方式之间相同相似部分相互参照即可，每个实施方式重点说明的都是与其他实施方式不同之处。

以上所述仅为本申请的几个实施方式，虽然本申请所揭露的实施方式如上，但所述内容只是为了便于理解本申请的技术方案而采用的实施方式，并非用于限定本申请。任何本申请所属技术领域的技术人员，在不脱离本申请所揭露的精神和范围的前提下，可以在实施方式的形式上及细节上作任何的修改与变化，但本申请的专利保护范围，仍须以所附权利要求书所界定的范围为准。