CN105825058A - 超稀疏雷达数据摄动补偿初轨计算方法 - Google Patents

Abstract

本发明提供了一种超稀疏雷达数据摄动补偿初轨计算方法，根据测量数据信息，建立时间转移方程，融合轨道选优思想，在二体意义下计算多圈转移情况下的方程解，在此基础上，考虑一定精度的摄动补偿，对时间转移方程的解进行修正迭代，直到得到满足解析摄动条件的初始轨道根数。本发明有效解决了一升一降两点数据时间间隔较长、传统的初轨计算方法不适用的问题，且满足轨道改进对初轨精度的要求，填补了我国在超稀疏雷达数据大时间跨度情况下的高精度初轨确定方面的空白。

Description

超稀疏雷达数据摄动补偿初轨计算方法

技术领域

本发明属于航天测量与控制领域，涉及一种初轨确定方法。

背景技术

超稀疏雷达数据是指测量设备对单个空间目标每圈次获取极少采样点数的测量数据，比如，美国的NAVSPASUR空间篱笆系统，空间目标每次穿过时只有极短的测量弧段(1至2个点)。当目标连续两次不同方向穿过观测设备波束时(一升一降或一降一升)，留下的两个位置矢量的时间间隔一般相差若干小时。超稀疏数据特点是测量点数极少(对于低轨目标，每天可能就4至6个点数据)，数据时间跨度大(数小时以上)。由于测量时间跨度大且点数稀疏，传统意义上的初轨计算方法已不适用。

发明内容

为了克服现有技术的不足，本发明提供一种适应超稀疏雷达数据融合轨道选优思想的摄动补偿定初轨方法，能够提高基于超稀疏雷达测量数据定初轨的精度及实用性。

本发明解决其技术问题所采用的技术方案包括以下步骤：

步骤一、给定起始位置r₁和目标位置r₂，r₁对应的测量时间为t₁，r₂对应的测量时间为t₂；

步骤二、建立时间转移方程n(t₂-t₁)＝2πN+(ε-sinε)-(δ-sinδ)，其中，μ为地球引力常数，a为轨道半长轴，N为转移圈数，若N＝0为不足一圈的转移，N≥1为多圈转移，ε,δ为Lagrange参数，且 σ＝|r₁-r₂|；

令F(a)＝n(t₂-t₁)-(ε-sinε)+(δ-sinδ)，利用牛顿迭代法求解，其中，迭代收敛后的a_k+1即为轨道半长轴a；

通过和确定轨道参数e，然后得到轨道参数M＝n(t₂-t₁)+E₁-esinE₁；

通过和解出参数e_p、e_Q，e_P＝(e_Px,e_Py,e_Pz)^T，e_Q＝(e_Qx,e_Qy,e_Qz)^T，通过以下方程组确定轨道参数Ω,i,ω：

e_n＝(e_nx,e_ny,e_nz)^T＝e_P×e_Q

e_nx＝sinΩsini

e_ny＝-cosΩsini

e_nz＝cosi

e_Pz＝sinωsini

e_Qz＝cosωsini；

对于多圈转移情况，ε、δ有不同的象限组合，因此，解算得到多个t₁时刻的轨道根数σ_m，m＝1,4，其中σ_m＝(a_me_mi_mΩ_mω_mM_m)^T；

步骤三、令以此为基准，对步骤二计算出的多个轨道进行选优，计算度量值Δa_m＝|a_m-a₀|，选取度量值最小的Δa_i对应σ_i为最佳轨道根数σ^j；

步骤四、利用步骤三计算的t₁时刻的最优根数σ^j，考虑摄动补偿，采用拟平均根数法进行轨道外推，计算出t₂时刻的理论瞬时根数拟平均根数法进行轨道外推的计算方法如下：

$\mathrm{\σ}\left(t\right)=\stackrel{\‾}{\mathrm{\σ}}\left(t_0\right)+(\mathrm{\δ}n+\mathrm{\σ}1+\mathrm{\σ}2)(t-t_0)+{\mathrm{\σ}}_l\left(t\right)+{\mathrm{\σ}}_s\left(t\right)$

$\stackrel{\‾}{\mathrm{\σ}}\left(t_0\right)=\mathrm{\σ}\left(t_0\right)-\[{\mathrm{\σ}}_l\left(t_0\right)+{\mathrm{\σ}}_s\left(t_0\right)\]$

其中，σ＝(aeiΩωM)^T为卫星的六个轨道根数，σ(t₀)和σ(t)分别为t₀和t时刻的瞬时轨道根数，为t₀时刻的平根数，δ＝(000001)^T，σ1和σ2分别为摄动的一阶和二阶长期项摄动，σ_l(t₀)和σ_s(t₀)分别是t₀时刻摄动的长周期项和短周期项摄动，σ_l(t)和σ_s(t)分别是t时刻摄动的长周期项和短周期项摄动，μ为地球引力常数；

将t₂时刻的理论瞬时根数转化为相对应的理论位置

步骤五、计算t₂时刻空间目标理论位置与目标实际位置r₂的位置偏差

步骤六、若|Δr|＜ε，ε为收敛门限，ε的取值范围(10^-6,10^-3)，则迭代结束，步骤三计算的σ^j即为满足解析摄动条件的初始轨道根数；否则，以作为目标位置，转入步骤一。

本发明的有益效果是：有效解决了一升一降两点数据时间间隔较长、传统的初轨计算方法不适用的问题，且满足轨道改进对初轨精度的要求，填补了我国在超稀疏雷达数据大时间跨度情况下的高精度初轨确定方面的空白。

具体实施方式

下面结合实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。

本发明根据测量数据信息，建立时间转移方程，融合轨道选优思想，在二体意义下计算多圈转移(两个测量点的时间差大于1个轨道周期)情况下的方程解，在此基础上，考虑一定精度的摄动补偿，对时间转移方程的解进行修正迭代，直到得到满足解析摄动条件的初始轨道根数。

本发明具体包括如下步骤：

步骤一：初始化状态量。

给定起始位置r₁，目标位置r₂，r₁对应的测量时间为t₁，r₂对应的测量时间为t₂。

步骤二：求解时间转移方程，得出二体意义下的轨道根数。

时间转移方程为：

n(t₂-t₁)＝2πN+(ε-sinε)-(δ-sinδ)(1)

其中，

μ为地球引力常数，a为轨道半长轴，

N为转移圈数，若N＝0为不足一圈的转移，N≥1为多圈转移

ε,δ为Lagrange参数，且

$sin\frac{\mathrm{\ϵ}}{2}=\±\sqrt{\frac{r_1+r_2+\mathrm{\σ}}{4a}}---\left(2\right)$

$sin\frac{\mathrm{\δ}}{2}=\±\sqrt{\frac{r_1+r_2-\mathrm{\σ}}{4a}}---\left(3\right)$

σ＝|r₁-r₂|

由时间转移方程计算轨道半长轴a，再依次计算其它轨道参数。

1)计算轨道半长轴a，令

F(a)＝n(t₂-t₁)-(ε-sinε)+(δ-sinδ)(4)

利用

$a_{k+1}=a_k-\frac{F\left(a_k\right)}{F^{\′}\left(a_{k+1}\right)},k=0,1,2...---\left(5\right)$

牛顿迭代法求解，其中，

$F^{\′}\left(a\right)=-\frac{3}{2}\frac{n}{a}(t_2-t_1)+\frac{r_1+r_2+\mathrm{\σ}}{2a^2}tan\frac{\mathrm{\ϵ}}{2}-\frac{r_1+r_2-\mathrm{\σ}}{2a^2}tan\frac{\mathrm{\δ}}{2}$

迭代收敛后的a_k+1即为方程的解。

2)确定轨道参数e,M的公式如下：

$e\sin E_1=\[(1-\frac{r_1}{a})cos(\mathrm{\ϵ}-\mathrm{\δ})-(1-\frac{r_2}{a})/\sin (\mathrm{\ϵ}-\mathrm{\δ})---\left(6\right)$

$e\cos E_1=1-\frac{r_1}{a}---\left(7\right)$

其中，ε,δ由公式(2)、(3)求得，根据公式(6)、(7)解得e,E₁。

M＝n(t₂-t₁)+E₁-esinE₁(8)

3)确定轨道参数Ω,i,ω的公式如下：

$r_1=a(\cos E_1-e)e_p+\sqrt{1-e^2}\sin E_1{e}_Q---\left(9\right)$

$r_2=a(\cos E_2-e)e_p+a\sqrt{1-e^2}\sin E_2{e}_Q---\left(10\right)$

由(9)、(10)解出e_p，e_Q

e_P＝(e_Px,e_Py,e_Pz)^Te_Q＝(e_Qx,e_Q,_ye^T

e_n＝(e_nx,e_ny,e_nz)^T＝e_P×e_Q(11)

e_nx＝sinΩsini(12)

e_ny＝-cosΩsini(13)

e_nz＝cosi(14)

e_Pz＝sinωsini(15)

e_Qz＝cosωsini(16)

对于多圈转移情况，ε、δ有不同的象限组合，因此，解算得到多个t₁时刻的轨道根数σ_i(i＝1,4)，其中σ_i＝(a_ie_ii_iΩ_iω_iM_i)^T。

步骤三：轨道选优。

令以此为基准，对步骤二计算出的多组轨道进行选优，计算度量值Δa_i＝|a_i-a₀|(i＝1,4)，选取度量值最小的Δa_i对应σ_i为最佳轨道根数，即σ^j＝σ_i。

步骤四：考虑摄动补偿，计算理论转移位置。

1)利用步骤三计算的t₁时刻的最优根数σ^j，考虑摄动补偿，采用拟平均根数法进行轨道外推，计算出t₂时刻的理论瞬时根数为满足精度要求，在轨道外推时需考虑地球非球形引力及大气阻力摄动的全部一阶项和部分二阶项。

拟平均根数法进行轨道外推的计算方法如下：

$\mathrm{\σ}\left(t\right)=\stackrel{\‾}{\mathrm{\σ}}\left(t_0\right)+(\mathrm{\δ}n+\mathrm{\σ}1+\mathrm{\σ}2)(t-t_0)+{\mathrm{\σ}}_l\left(t\right)+{\mathrm{\σ}}_s\left(t\right)---\left(17\right)$

$\stackrel{\‾}{\mathrm{\σ}}\left(t_0\right)=\mathrm{\σ}\left(t_0\right)-\[{\mathrm{\σ}}_l\left(t_0\right)+{\mathrm{\σ}}_s\left(t_0\right)\]---\left(18\right)$

其中

σ＝(aeiΩωM)^T为卫星的六个轨道根数

σ(t₀)和σ(t)分别为t₀和t时刻的瞬时轨道根数

为t₀时刻的平根数

δ＝(000001)^T

σ1和σ2分别为摄动的一阶和二阶长期项摄动

σ_l(t₀)和σ_s(t₀)分别是t₀时刻摄动的长周期项和短周期项摄动

σ_l(t)和σ_s(t)分别是t时刻摄动的长周期项和短周期项摄动

$n=\sqrt{u}{a}^{-3/2}$

μ为地球引力常数

2)将t₂时刻的理论瞬时根数转化为相对应的理论位置

步骤五：摄动补偿。

由于摄动影响的存在，t₂时刻空间目标理论位置与目标实际位置r₂存在着偏差，计算位置偏差

步骤六：对转移方程的解进行迭代修正。

1)若|Δr|＜ε(ε为收敛门限，ε的取值范围(10^-6,10^-3))，迭代结束，步骤三计算的σ^j即为满足解析摄动条件的初始轨道根数；

2)否则，以作为目标位置，转入步骤一。

以某观测站数据为例：

设定某一观测站的位置，给定空间目标的瞬时轨道数据，利用精密轨道外推模型进行外推2天，找到在该观测站观测弧段内的一升一降两点数据，并转换成测站坐标系下的测量数据(测距、方位角、俯仰角)，加入一定的随机差和系统差(测距：随机差50m，系统差20m；方位角和俯仰角：随机差0.01°，系统差0.01°)。

观测站位置如表2：

表2测站的大地坐标系位置

大地经度(度)	大地纬度(度)	大地高度(米)
			107.5	25.0	500.0

目标的理论J2000瞬时根数如表3：

表3目标的理论J2000瞬时根数

历元时间(北京时)	a(米)	e	i(°)	Ω(°)	ω(°)	M(°)
							2010-12-610:8:2.78188	7009411.313	0.0015060	97.804	331.332	95.134	290.088

两个观测点数据如表4：

表4测量数据

初轨计算方法如下：

1)初始化状态量。根据两点的测量信息，进行坐标转换，由测站大地坐标系转换到J2000惯性系，得到J2000惯性系下的r₁、r₂，以及r₁对应的测量时间为t₁，r₂对应的测量时间为t₂。

t₁＝2010年12月6日10时8分2.78188秒

t₂＝2010年12月6日22时4分25.53687秒

表5两点的J2000惯性系下的空间位置

空间位置	x(米)	y(米)	z(米)
				r1	5375202.47183465	-3398046.20955139	2940011.19463203
r2	-5779854.27886388	2636625.44128764	2951945.77024189

2)求解二体意义下的时间转移方程

根据步骤二，根据不同的ε、δ象限组合，得到4个迭代最终值：

a₁＝6859615.85255609m

a₂＝7009544.21694125m

a₃＝6857243.25848122m

a₄＝7006331.46480324m

根据a不同值求解其它的轨道参数，结果如表6所示：

表6ε、δ不同象限组合下的J2000惯性系轨道根数

组	历元时间(北京时)	a(米)	e	i(°)	Ω(°)	ω(°)	M(°)5 -->
								1	2010-12-6 10:8:2.78188	6859615.852	0.377714635	98.325	331.579	269.936	71.630
2	2010-12-6 10:8:2.78188	7009544.216	0.001336095	98.325	331.579	92.526	292.709
								3	2010-12-6 10:8:2.78188	6857243.258	0.0574032547	81.674	151.579	270.114	250.850
4	2010-12-6 10:8:2.78188	7006331.464	0.4252122915	81.674	151.579	270.062	294.370

3)轨道选优

对于低轨空间目标，偏心率一般很小(e≤0.05)，a若越接近a₀，则a对应的轨道越接近真实轨道。

a₀＝7005557.730，得出第2组轨道为最优轨道。

4)理论位置的计算

超稀疏雷达测量数据，主要是低轨目标，而低轨目标主要摄动源为地球非球形引力摄动及大气阻力摄动。对于轨道高度为500km的空间目标，地球非球形引力摄动影响量级为O(10^-3)，大气阻力摄动影响量级为O(10^-6)，其它摄动力摄动影响量级大于O(10^-6)，故结合实效性要求及初轨精度要求，摄动源只考虑地球非球形引力摄动及大气阻力摄动。

为了使目标轨道计算达到10-⁶精度(对于位置误差10m，速度误差1cm/s)，需考虑运动方程的全部一阶解(包括一阶、二阶长期摄动和一阶周期摄动)和部分二阶解。为避免小偏心率问题，采用第一类无奇点根数σ＝(a,i,Ω,ξ＝ecosω,η＝-esinω,λ＝M+ω)，同时，为了避免临界倾角问题，采用拟平均根数法。

将第2组轨道作为轨道外推初值，采用步骤四中拟平均根数的解析法外推，得出t₂时刻理论空间目标位置及理论轨道根数，空间目标位置如表7、表8所示：

表7t₂时刻理论位置

x(米)	y(米)	z(米)
			5375202.471	-3398046.209	2940011.194

表8理论轨道

历元时间(北京时)	a(米)	e	i(°)	Ω(°)	ω(°)	M(°)
							2010-12-610:8:2.78188	7009501.421	0.001337250	98.326	332.104	85.749	68.907

5)摄动补偿计算

根据步骤五，计算Δr结果如表9所示：

表9t₂时刻理论位置与外推起始值偏差

Δx(米)	Δy(米)	Δz(米)
			-24963.228	-52418.978	-343.565

计算位置偏差：

若|Δr|＜ε，则表8则为最终初轨，否则，令转入步骤一，进行再次迭代。

6)最终迭代结果即为满足解析摄动条件的初始轨道根数，结果如下：

表10t₂时刻理论位置

x(米)	y(米)	z(米)6 -->
			-5756756.663	2686376.991	2949836.949

表11初轨根数

历元时间(北京时)	a(米)	e	i(°)	Ω(°)	ω(°)	M(°)
							2010-12-6 10:8:2.78188	7009445.660	0.0014785	97.806	331.334	95.365	289.854

计算得出的初始轨道根数与实际轨道根数比较如表12

表12轨道误差比较

Δa(米)	Δe	Δi(°)	ΔΩ(°)	Δω(°)	ΔM(°)
						37.347	0.0000275	0.002	0.002	0.231	0.234

为了进一步验证该算法，利用6000多个低轨空间目标进行初轨计算，结果表明，收敛率达到99％。利用一升一降两点所定初轨与实际轨道根数精度比较统计结果见表13。

表13精度统计结果

采用该方法计算的初轨完全能够满足后续轨道改进对初轨的精度要求。

Claims (1)

1.一种超稀疏雷达数据摄动补偿初轨计算方法，其特征在于包括下述步骤：

步骤一、给定起始位置r₁和目标位置r₂，r₁对应的测量时间为t₁，r₂对应的测量时间为t₂；

步骤二、建立时间转移方程n(t₂-t₁)＝2πN+(ε-sinε)-(δ-sinδ)，其中，μ为地球引力常数，a为轨道半长轴，N为转移圈数，若N＝0为不足一圈的转移，N≥1为多圈转移，ε,δ为Lagrange参数，且 σ＝|r₁-r₂|；

令F(a)＝n(t₂-t₁)-(ε-sinε)+(δ-sinδ)，利用牛顿迭代法求解，其中，迭代收敛后的a_k+1即为轨道半长轴a；

通过和确定轨道参数e，然后得到轨道参数M＝n(t₂-t₁)+E₁-esinE₁；

通过和解出参数e_p、e_Q，e_P＝(e_Px,e_Py,e_Pz)^T，e_Q＝(e_Qx,e_Qy,e_Qz)^T，通过以下方程组确定轨道参数Ω,i,ω：

e_n＝(e_nx,e_ny,e_nz)^T＝e_P×e_Q

e_nx＝sinΩsini

e_ny＝-cosΩsini

e_nz＝cosi

e_Pz＝sinωsini

e_Qz＝cosωsini；

对于多圈转移情况，ε、δ有不同的象限组合，因此，解算得到多个t₁时刻的轨道根数σ_m，m＝1,4，其中σ_m＝(a_me_mi_mΩ_mω_mM_m)^T；

步骤三、令以此为基准，对步骤二计算出的多个轨道进行选优，计算度量值Δa_m＝|a_m-a₀|，选取度量值最小的Δa_i对应σ_i为最佳轨道根数σ^j；

步骤四、利用步骤三计算的t₁时刻的最优根数σ^j，考虑摄动补偿，采用拟平均根数法进行轨道外推，计算出t₂时刻的理论瞬时根数拟平均根数法进行轨道外推的计算方法如下：

$\mathrm{\σ}\left(t\right)=\stackrel{\‾}{\mathrm{\σ}}\left(t_0\right)+(\mathrm{\δ}n+\mathrm{\σ}1+\mathrm{\σ}2)(t-t_0)+{\mathrm{\σ}}_l\left(t\right)+{\mathrm{\σ}}_s\left(t\right)$

$\stackrel{\‾}{\mathrm{\σ}}\left(t_0\right)=\mathrm{\σ}\left(t_0\right)-\[{\mathrm{\σ}}_l\left(t_0\right)+{\mathrm{\σ}}_s\left(t_0\right)\]$

其中，σ＝(aeiΩωM)^T为卫星的六个轨道根数，σ(t₀)和σ(t)分别为t₀和t时刻的瞬时轨道根数，为t₀时刻的平根数，δ＝(000001)^T，σ1和σ2分别为摄动的一阶和二阶长期项摄动，σ_l(t₀)和σ_s(t₀)分别是t₀时刻摄动的长周期项和短周期项摄动，σ_l(t)和σ_s(t)分别是t时刻摄动的长周期项和短周期项摄动，μ为地球引力常数；

将t₂时刻的理论瞬时根数转化为相对应的理论位置

步骤五、计算t₂时刻空间目标理论位置与目标实际位置r₂的位置偏差

步骤六、若|Δr|＜ε，ε为收敛门限，ε的取值范围(10^-6,10^-3)，则迭代结束，步骤三计算的σ^j即为满足解析摄动条件的初始轨道根数；否则，以作为目标位置，转入步骤一。