一种依据天基卫星测角资料确定低轨目标卫星初始轨道的方法

摘要

本发明公开了一种依据天基卫星测角资料确定低轨目标卫星初始轨道的方法，该方法利用高速运动低轨目标卫星相对于低轨观测卫星的测角资料，得到目标星相对观测星单位矢量的条件方程。在给条件方程加约束的前提下得到测距ρ0的条件方程和利用选择一定步长、给约束方程中的ρ0以及目标星的半长轴a0加一定权重的方法，得到迭代解作为下一步迭代的初值：位置速度进而迭代计算直到收敛得到低轨目标卫星初始轨道，由此屏蔽假解，得到真解，确定目标星的位置和速度。本发明通过对求解目标加约束条件的方式，解决了由于目标星相对观测星的高速动态性和距离观测量的缺失导致的假解问题，成功使计算程序收敛到真解附近，使定轨结果接近于理论最佳解C.R下界。

权利要求

1.一种依据天基卫星测角资料确定低轨目标卫星初始轨道的方法，其特征在于，所述方法采用天基测角资料确定低轨卫星初始轨道的运动方程，并利用该运动方程使得计算程序对所有低轨运动目标求解初始轨道时收敛到真解；所述方法首先利用高速运动低轨目标卫星相对于低轨观测卫星的测角资料，得到目标星相对观测星单位矢量 的条件方程；

接着，在给条件方程加约束的前提下得到测距ρ0加约束的条件方程 和最后，利用选择一定步长、给约束方程中的ρ0以及目标星的半长轴a0加一定权重的方法，得到迭代解作为下一步迭代的初值：位置 速度 进而迭代计算直到收敛得到低轨目标卫星初始轨道，由此屏蔽假解，得到真解，确定目标星的位置和速度；

具体的实现步骤如下：

步骤1，

依据地基定初轨的单位矢量法，得到目标相对观测星的单位矢量

式中 为目标星相对于观测星的单位矢量，ρk为某一时刻目标星相对于观测星的距离，观测量(αk,δk)为 在惯性坐标系下的赤经赤纬；

步骤2，

引进两个与 正交的单位矢量

单位矢量 由观测量赤经赤纬(αk,δk)计算，由此对每个观测量(αk,δk)产生下述两个加约束的条件方程：其中加权1/ρk*σ

式中σ为迭代过程中得到的中误差；

步骤3，

增加对ρ0的约束条件；

初始迭代只对测距ρ0加约束求解 对ρ0加约束条件方程为：加权1/wr,(wr＝100米)其解作为以下迭代计算的初值

第二次迭代时，改为求解 的改正量

此时对每个观测量(αk,δk)形成下述两个条件方程：

式中 由 计算，而 由 计算；

步骤4,

依据 再对测距ρ0加约束条件，得到测距ρ0加约束的条件方程：

ρ0的初值由圆轨道计算得到，将条件方程加上述约束条件的解作为下一步迭代的初值步骤5，在迭代计算后期对半长轴a0再加入一个约束条件，需要适当加权，以匹配ρ0，以提高定轨精度。

说明书

一种依据天基卫星测角资料确定低轨目标卫星初始轨道的

方法

技术领域

[0001] 本发明涉及人造卫星精密定轨中的卫星定初轨技术领域，具体涉及依据观测星测得的低轨目标卫星测角资料确定其初始轨道的方法。

背景技术

[0002] 依据天基卫星测角资料确定低轨目标卫星轨道，最初认为与地基观测站测角资料定初轨没有什么区别，仅仅是把观测站从地面搬到天上而已。事实上这两者有很大区别，地基测站仅随地球自转，而观测卫星则在低轨轨道上高速运行，其运动的力学规律和低轨目标卫星极其类似。当依据观测星对低轨目标卫星的测角资料计算初轨时，由于没有距离观测量，将出现两重解，真解是目标卫星轨道，假解是观测卫星自身轨道。

[0003] 如果照搬地基观测的初轨计算方法，对半长轴20000km以上的高轨目标，可以收敛到真解附近，这是由于高轨卫星运动速度较慢，和观测星的运动规律有较大差异，但是对于低轨目标卫星，由于运动力学规律的极度相似性，以及目标星相对观测星的高速动态性和距离观测量缺失，如果不对传统的测角资料定初轨运动方程加以改进，会使得计算程序自动收敛到假解附近，这就会严重影响定轨精度。

发明内容

[0004] 针对现有天基确定低轨卫星目标轨道技术会自动收敛到假解附近，严重影响定轨精度的问题，本发明的目的在于提供一种依据低轨目标卫星相对观测星测角资料定初轨的方法，该方法屏蔽传统的依据天基测角资料定初轨运动方程的假解，使程序向着真解的方向收敛，并达到较高的定轨精度。

[0005] 为了达到上述目的，本发明采用如下的技术方案：

[0006] 一种依据天基卫星测角资料确定低轨目标卫星初始轨道的方法，该方法采用加权测角资料确定低轨卫星初始轨道的运动方程，并利用该运动方程使得计算程序对所有目标收敛到真解。

[0007] 在本方案的优选方案中，所述方法首先利用高速运动低轨目标卫星相对于低轨观测卫星的测角资料，得到目标星相对观测星单位矢量 的条件方程；

[0008] 接着，在给条件方程加约束的前提下得到测距ρ0加约束的条件方程和

[0009] 最后，利用选择一定步长、给约束方程中的ρ0以及目标星的半长轴a0加一定权重的方法，得到迭代解作为下一步迭代的初值：位置 速度 进而迭代计算直到收敛得到低轨目标卫星初始轨道，并由此屏蔽假解，得到真解，确定目标星的位置和速度。

[0010] 进一步的，具体的实现步骤如下：

[0011] 步骤1，

[0012] 依据地基定初轨的单位矢量法，得到目标相对观测星的单位矢量

[0013]

[0014] 式中 为目标星相对于观测星的单位矢量，ρk为目标星相对于观测星的距离，观测量(αk,δk)为 在惯性坐标系下的赤经赤纬；

[0015] 步骤2，

[0016] 引进两个与 正交的单位矢量

[0017]

[0018] 单位矢量 由观测量赤经赤纬(αk,δk)计算，由此对每个观测量(αk,δk)产生下述两个加约束的条件方程：

[0019] 其中加权

[0020] 式中σ为迭代过程中得到的中误差；

[0021] 步骤3，

[0022] 增加对ρ0的约束条件；

[0023] 初始迭代只对测距ρ0加约束求解 对ρ0加约束条件方程为：

[0024] 加权1/wr,(wr＝100米)

[0025] 其解作为以下迭代计算的初值

[0026] 第二次迭代时，改为求解 的改正量

[0027]

[0028] 此时对每个观测量(αk,δk)形成下述两个条件方程：

[0029] (加权1/ρk*σ)

[0030] 式中 由 计算，而 由 计算；

[0031] 步骤4,

[0032] 依据 再对测距ρ0加约束条件，得到测距ρ0加约束的条件方程：

[0033]

[0034] ρ0的初值由圆轨道计算得到，将条件方程加上述约束条件的解作为下一步迭代的初值

[0035] 步骤5，

[0036] 在迭代计算后期对半长轴a0再加入一个约束条件，需要适当加权，以匹配ρ0，以提高定轨精度。

[0037] 利用本发明解决了对于低轨目标卫星，由于运动力学规律的相似，目标相对观测星的高速动态性和距离观测量缺失，使得计算程序自动收敛到假解附近的数学奇异问题，取得了判断解的真伪的有益效果，并有效的提高了定初轨精度。

附图说明

[0038] 以下结合附图和具体实施方式来进一步说明本发明。

[0039] 图1为目标卫星与观测星的几何关系图；

[0040] 图2为本发明的计算流程图；

[0041] 图3为距离观测星1000km左右的低轨目标的定初轨结果；

[0042] 图4为距离观测星1600km左右的低轨目标的定初轨结果；

[0043] 图5为距离观测星2000km左右的低轨目标的定初轨结果。

具体实施方式

[0044] 为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体图示，进一步阐述本发明。

[0045] 本发明利用定初轨的通用技术之一的单位矢量法，得到目标星相对观测星的单位矢量 (下标k为采样点的序号)的方程，其中单位矢量 方程的解有两个加约束的条件方程：

[0046]

[0047] 由此，屏蔽假解得到真解，确定目标星的位置和速度。其具体的实现过程如下：

[0048] 参见图1，其所示为目标卫星与观测星的几何关系图。由图可知，

[0049] 式1

[0050] 式中ρk为目标与测站的距离。

[0051] 由计算程序依据地基定初轨的通用技术之一单位矢量法，得到目标相对观测星的单位矢量 (如式2)，下式中(αk,δk)为观测量，表征 方向的惯性坐标系下的赤经赤纬。

[0052] 式2

[0053] 按照二体问题的动力学模型，tk时刻的 与t0时刻的 有如下关系

[0054] 式3

[0055] 式中，fk,gk是t0,tk, 的函数。

[0056] 由上述推导可得：

[0057] 式4

[0058] 为了消去未知量ρk，形成条件方程，引进两个与 正交的单位矢量 作为参考矢量：

[0059] 式5

[0060] 这两个矢量由观测量赤经赤纬(αk,δk)计算，分别用 点乘式4两边，消去了项。这样，对每个观测量(αk,δk)，将产生下述两个条件方程：

[0061] (加权1/ρkσ)  式6

[0062] 上式中σ为迭代过程中得到的中误差。

[0063] 本发明中的定初轨指的是依据观测星测得的测角资料确定初轨，目标是半长轴为9000km以下的低轨目标卫星。当低轨目标卫星的轨道高度与观测星相差不多时，出现了即使增加了对星间距离ρ0(下标0代表初始值)的约束条件仍不收敛的情形。由此，在单位矢量法中，同时增加对ρ0与a0的两个约束条件，对低轨目标卫星将有很好的收敛性。其过程如下：

[0064] 初始迭代仅对ρ0加约束，求解 对ρ0加约束的条件方程为：

[0065] 加权1/wρ,(wρ＝100米)  式7

[0066] 它的解作为以下迭代计算的初值

[0067] 第二次迭代时，改为求解 的改正量 即参考矢量法。

[0068] 式8

[0069] 分别用 点乘 (该公式为式4的改正量，即改正后与 与 相减得到的量值，式中带*的符号为每次的计算迭代值，为已知量)两边，此时对每个观测量(αk,δk)形成下述两个条件方程：

[0070] (加权1/ρk*σ)  式9

[0071] 式中 由 计算，而 由 计算。

[0072] 依据 以及图1所示的几何关系，再对测距ρ0加约束条件，将形成下述条件方程：

[0073] 加权1/wρ,(wρ＝100米)  式10

[0074] ρ0的初值由圆轨道计算得到，将条件方程加上述约束条件的解作为下一步迭代的初值

[0075] 这里涉及的圆轨道初轨算法，是指假定目标轨道为圆轨道，即rk＝r0＝a，这是一个对a进行搜索的算法。仅利用观测弧段的第一点(t1,α1,δ1)与最后一点(tN,αN,δN)或等价的由

[0076]

[0077] 两边平方得：

[0078]

[0079] 现在令r1＝a，可解出：

[0080]

[0081] 同理

[0082]

[0083]

[0084] 即如已知a可解出ρ1,ρN，于是得 可以计算出平均角速度：

[0085]

[0086] 另一方面，由n2a3＝μ，又可计算出一个n'，如果半长轴a的值正确，应该有[0087]

[0088] 对半长轴a进行搜索，即找到一个a，使|n*-n'|达最小，求得半长轴a后，就得到了与 由

[0089]

[0090]

[0091] 可解出 从而获得了初轨。上式中f1,g1,fN,gN是卫星状态t0,t1,tN, 的函数。

[0092] 从第四次迭代开始(即迭代后期)，对半长轴a0再加入一个约束条件(需要适当加权，以匹配ρ0)，它形成下述条件方程：

[0093] 加权1/wa,(wa＝100米)  式11

[0094] 式中：

[0095]

[0096] 式12

[0097] 该两式右端都代入迭代值 a0*计算，a0*由下式计算：

[0098]   式13

[0099] 上述方案中通过对ρ0与a0增加两个约束条件，使得低轨目标卫星将有很好的收敛性。以下具体介绍三种加约束的方法：

[0100] (1)第1种方法

[0101] 对历元时刻(需要确定初轨的时刻)的ρ0与a0同时加约束，其中对ρ0与a0的权重都取100米之倒数。

[0102] 对历元时刻半长轴a0进行搜索，一个a0对应一个ρ0，对a0的搜索空域确定为9000km到6600km。将分三层搜索，第一层搜索以7800km为中心，上下1200km，搜索步长为5km，第二层用上一层搜索的最佳点为中心，上下5km空域，步长为0.5km，最后找到最佳点a0与对应的ρ0作为约束条件。

[0103] (2)第2种方法

[0104] 对历元时刻的ρ0与a0同时加约束，其中对ρ0的加权为100米之倒数，而对a0的权重取为15km之倒数。

[0105] (3)第3种方法

[0106] 仅对历元时刻ρ0加约束，对它的权重为100米之倒数。

[0107] 后两种方法对a0的搜索起点都取为第一种方法对a0的搜索解，搜索范围为上下75km。分两层搜索，第一层搜索步长为3km，第二层搜索步长为0.3km。

[0108] 参见图2，其所示为基于上述三种加约束方法进行计算的流程图。由图，整个流程如下：

[0109] 1)首先计算第1种方法，对历元时刻(需要确定初轨的时刻)的ρ0与a0同时加约束，其中对ρ0与a0的权重都取100米之倒数，具体过程如上所述，此处不赘述。

[0110] 2)对少数目标不能收敛到真解的情况下，计算第3种方法，采取仅对ρ0加约束，对它的权重取为100米之倒数。对a0的搜索起点取第1种方法对a0的搜索解，搜索范围为上下75km，再分两层搜索，第一层搜索步长为3km，第二层搜索步长为0.3km，判断是否收敛，如果收敛，则计算结束，用它的解作为初轨；如果不收敛，则转入步骤3)。

[0111] 3)当第3种方法不收敛时，计算第2种方法，采用ρ0与a0同时加约束，其中对ρ0的加权为100米之倒数，而对a0的权重取为15km之倒数，判断是否收敛。如果收敛，则计算结束，用它的解作为初轨；如果不收敛，则转入步骤4。

[0112] 4)如果后两种方法都不收敛，而第1种方法观测值与计算值的差(o-c)的中误差小于50σ，则计算结束，用第1种方法的解作为初轨；

[0113] 5)如果后两种方法都不收敛，且第1种方法观测值与计算值的差(o-c)的中误差大于50σ，则用圆轨道算法计算初轨。圆轨道法假设卫星运动角速度相同，收敛条件为|n*-n'|≤ε，具体的方法如上所述。

[0114] 6)最后判断真假得到真解，计算结束。

[0115] 对于初轨计算问题，若观测量赤经、赤纬αk,δk为均值为0，方差为σ'的正态分布随机量，C.R下界为：

[0116] 式14

[0117] 初轨计算的参数θ为六维向量，可取为 也可取为轨道根数。初始状态θ、时间t已知的前提下，函数α(θ,t)和δ(θ,t)可求，利用正确无误差的θ和t，即可计算出初轨的C.R下界。由于初轨计算是非线性模型，无法达到C.R下界这一最佳精度，只能接近最佳精度。

[0118] 由C.R下界的表达式可知，它由两部分乘积组成，一部分为观测误差，一部分为定初轨法方程系数矩阵的逆。因此，判断某一方法是否为最佳精度算法，应有两个判据：

[0119] (1)初轨精度与观测误差成准线性关系；

[0120] (2)初轨精度已接近C.R下界。

[0121] 由C.R下界公式可以看出，对于不加约束条件的初轨计算，其精度影响因素如下：

[0122] (1)观测误差σ'：σ'越小，即测角精度越高，初轨精度越高；

[0123] (2)采样频率f：采样频率f越高，即相同时间内的采样数据点N越多，初轨精度越高；

[0124] (3)观测弧段时间跨度T：观测时间越长，即相同采样频率下的采样数据点N越多，初轨精度越高。

[0125] 如图3、图4、图5分别为距离观测星1000km左右、1600km左右、2000km左右低轨目标卫星的定初轨结果图。仿真基于的观测星轨道高度为660km，降交点地方时6:30的太阳同步轨道，光轴指向为轨道面法向负方向，仿真观测数据为目标相对于观测星的矢量在惯性空间中的测角资料。以下将提供实例给予本发明方案的进一步了解。

[0126] 实施例1

[0127] 卫星的轨道根数如下：

[0128] a＝7203.709km

[0129] i＝80°

[0130] Ω＝130.14°

[0131] e＝0.00187

[0132] λ＝M+ω＝272.59°

[0133] ρ0＝1526.9km

[0134] ρmin＝635.4km

[0135] 测角精度19角秒，数据479秒，240点，每2秒一个点，初轨精度仿真结果以及和C.R下界的关系见图3。

[0136] 实施例2

[0137] 目标星轨道根数如下：

[0138] a＝7311.066km

[0139] i＝70°

[0140] Ω＝143.93°

[0141] e＝0.0016

[0142] λ＝M+ω＝272.95°

[0143] ρ0＝1840.2km

[0144] ρmin＝1187.3km

[0145] 测角精度19角秒，数据258秒，129点，每2秒一个点，初轨精度仿真结果以及和C.R下界的关系见图4。

[0146] 实施例3

[0147] 卫星的轨道根数如下：

[0148] a＝7405.907km

[0149] i＝70°

[0150] Ω＝146.75°

[0151] e＝0.00133

[0152] λ＝M+ω＝275.11°

[0153] ρ0＝2241.6km

[0154] ρmin＝1523.4km

[0155] 测角精度19角秒，数据281秒，141点，每2秒一个点，初轨精度仿真结果以及和C.R下界的关系见图5。

[0156] 由仿真结果可见，通过本发明提供的方案能够将计算程序结果收敛在真解附近，而且和理论最佳解C.R下界差别不大，充分证明此方法的正确可行性。

[0157] 以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。