航天器相对位置增量定轨方法
【专利摘要】航天器相对位置增量定轨方法，属于航天器轨道确定技术领域。本发明是为了解决现有基于测量历元差分的定轨方法中，需要采用其他定轨方法获得初始位置信息的问题。本发明所述的航天器相对位置增量定轨方法，首先，确定航天器在三维空间中运动的轨道平面；然后，将三维相对位置矢量转换为二维平面内航天器相对位置坐标，在该二维平面内采用最小二乘原理进行航天器运动轨道形状的拟合，采用拉格朗日乘子法求解拟合矩阵的特征值；最后，根据拟合矩阵的特征值判断航天器的轨道类型，并根据中心天体的位置，获得航天器在惯性参考系下的绝对位置信息，实现航天器的初轨确定。本发明适用于对航天器实现定轨。
【专利说明】
航天器相对位置増量定轨方法
技术领域
[0001 ]本发明属于航天器轨道确定技术领域。
【背景技术】
[0002] 航天器轨道确定的任务就是应用少量的观测数据，快速确定航天器两体问题下的 轨道要素。观测数据由于测站数目的不同，测量量(如:位置、速度、角度)属性的差别，相应 地定轨方法也有所不同。其中位置观测量可通过全球卫星导航系统或地面雷达获得，测量 得到的位置信息均可通过简单的转换获得航天器在地球惯性坐标系下的位置矢量或方向 矢量，采用三个地心位置矢量可由吉布斯方法定轨，采用两个位置矢量及他们之间的时间 间隔定轨的兰伯特问题。在脉冲星导航系统及全球卫星导航系统中，利用观测数据的历元 差分方法可有效消除测量部分的公共误差提供更精确的测量信息，然而在基于测量历元差 分的导航方法中，初始位置信息的估计往往先采用地心位置矢量通过其他定轨方法给出。 为此使得基于历元差分相对位置信息导航算法无法独立应用，需依靠其他定轨算法提供可 靠的初始位置估计后方可进行。
【发明内容】
[0003] 本发明是为了解决现有基于测量历元差分的定轨方法中，需要采用其他定轨方法 获得初始位置信息的问题，现提供航天器相对位置增量定轨方法。
[0004] 航天器相对位置增量定轨方法，该方法为：
[0005] 首先，确定航天器在三维空间中运动的轨道平面；
[0006] 然后，将三维相对位置矢量转换为二维平面内航天器相对位置坐标，在该二维平 面内采用最小二乘原理进行航天器运动轨道形状的拟合，并采用拉格朗日乘子法求解拟合 矩阵的特征值；
[0007] 最后，根据拟合矩阵的特征值判断航天器的轨道类型，若航天器轨道为椭圆轨道， 则中心天体位于椭圆的一个焦点处，结合开普勒时间方程与测量相对应的时间序列确定中 心天体的位置;若航天器轨道为圆轨道，则中心天体在轨道的中心;根据中心天体的位置， 获得航天器在惯性参考系下的绝对位置信息，实现航天器的初轨确定。
[0008] 本发明所述的航天器相对位置增量定轨方法，是一种直接使用历元差分测量的相 对距离增量信息进行定轨的方法，仅依靠航天器相对于空间任意未知位置的相对位置增量 及相应的观测时间。该方法使得历元差分测量信息在导航应用时不依赖与其他测量信息， 提高了基于历元差分测量信息的导航系统的自主性，本发明定轨方法简单，易于实现。
【附图说明】
[0009] 图1为基于历元差分的航天器相对位置示意图；
[0010] 图2为轨道平面的各参考系示意图；
[0011] 图3为本发明所述的航天器相对位置增量定轨方法的流程图。
【具体实施方式】
【具体实施方式】 [0012] 一:参照图3具体说明本实施方式，本实施方式所述的航天器相对位 置增量定轨方法，该方法为：
[0013] 首先，确定航天器在三维空间中运动的轨道平面；
[0014] 然后，将三维相对位置矢量转换为二维平面内航天器相对位置坐标，在该二维平 面内采用最小二乘原理进行航天器运动轨道形状的拟合，并采用拉格朗日乘子法求解拟合 矩阵的特征值；
[0015] 最后，根据拟合矩阵的特征值判断航天器的轨道类型，若航天器轨道为椭圆轨道， 则中心天体位于椭圆的一个焦点处，结合开普勒时间方程与测量相对应的时间序列确定中 心天体的位置;若航天器轨道为圆轨道，则中心天体在轨道的中心;根据中心天体的位置， 获得航天器在惯性参考系下的绝对位置信息，实现航天器的初轨确定。
【具体实施方式】 [0016] 二:参照图1具体说明本实施方式，本实施方式是对一 所述的航天器相对位置增量定轨方法作进一步说明，本实施方式中，确定航天器在三维空 间中运动的轨道平面的具体方法为：
[0017]测量航天器在地心惯性坐标系下历元时刻与之间航天器移动的相对位置矢 量Δ r (t j-i，t j)，并根据该相对位置矢量Δ r (t j-i，t j)获得航天器的初始位置r (to)，其中j = 1,2,3,...；
[0018]根据相对位置矢量ArUj-i,!:」），采用最小二乘原理确定单位角动量矢量h;
[0020]其中，hx、hy、hz分别为单位角动量矢量h在空间三维参考系中的X轴、y轴、z轴的分 量，Arx(tj-1，tj)、Ary(tj-1，tj)、Arz(tj-1，tj)分别为相对位置矢量 Δ r(tj-1，tj)在三维参 考系中x轴、y轴、z轴的分量，
[0021 ]根据下式利用单位角动量矢量h获得轨道倾角i和升交点赤经Ω，
[0024]利用轨道倾角i和升交点赤经Ω确定三维轨道平面坐标。
[0025]图1描述了相对位置测量信息，r为航天器在地心惯性坐标系下相对于地心的位置 矢量，Δ r(tj-i，tj)和△ r(tj，tj+i)分别为历元时刻（tj-i，tj)和（tj，tj+i)间航天器移动的相 对位置矢量，定轨就是根据一系列的测量信息ArUntJ及t」(j = l，2,3，...）获得航天器 的初始位置r (to)及速度或轨道六根数。
[0026]设相对位置Ar为地心惯性J2000参考系下的矢量。根据向量间关系：任意两个矢 量差乘即可获得轨道平面的法向量，即单位角动量矢量，这里考虑测量噪声的存在，利用所 有观测量，采用最小二乘原理确定单位角动量矢量h。上述式(1)中可通过求解目标函数：
[0028]来获得约束条件G(hx，hy，hz)下的最小值。并根据非线性优化理论，选取适当的优 化算法(如牛顿法、最速下降法等)获得单位角动量矢量h。
【具体实施方式】 [0029] 三:本实施方式是对二所述的航天器相对位置增量定 轨方法作进一步说明，本实施方式中，将三维相对位置矢量转换为二维平面内航天器相对 位置坐标的具体方法为：
[0030] 设参考系Pxyz以航天器的初始位置为原点，z轴垂直于轨道平面，x，y轴由地球 J2000参考系根据坐标转换矩阵Q Xx获得；
[0031 ]根据轨道倾角i和升交点赤经Ω获得航天器相对位置在参考系Pxyz中的矢量&:
[0032] xk = Qxxmk (4)
[0033]
为k时刻航天器位置相对于初始位置的矢量，QXx为坐标转 换矩阵，具体公式如下：
[0035]利用矢量a在x，y轴上的投影获得航天器各点在二维平面内的位置坐标(xk，yk)。
[0036] 图2给出了在轨道平面内各参考系之间的关系，其中参考系Pxyz以航天器未知的初 始位置为原点，z轴垂直于轨道平面，x，y轴由地球J2000参考系根据坐标转换矩阵Q Xx获得。 参考系0XW以椭圆轨道中心为原点，z轴垂直于轨道平面，x，y轴分别沿椭圆长轴和短轴方 向，满足右手定则。参考系〇 ' xyz以地球中心为原点，为轨道平面内的近焦点坐标系。
【具体实施方式】 [0037] 四：本实施方式是对一所述的航天器相对位置增量定 轨方法作进一步说明，本实施方式中，采用最小二乘原理进行航天器运动轨道形状拟合的 具体方法为：
[0038]设平面内任意轨道的表达式为：
[0040] 其中，U=[A，B，C，D，E，M]，X=[x2，Xy，y2， X，y，l]，A、B、C、D、E、M 均为多项式各项的 系数，X和y分别为二维平面参考系Pxyz中轨道上各点坐标参数，
[0041] 根据最小二乘原理将公式(6)转换为以下形式：
[0043] 式中；，可由矢量xk在x，y轴上的投影获得航天器各点在 二维平面内的位置坐标(Xk，yk)计算获得。
[0044] 本实施方式公式(6)中含有的：
[0045] 4AC-B2=1 (7)
[0046] 公式(7)是为了避免全零解及抛物线、双曲线情况，所增加的约束条件。
[0047] 平面椭圆拟合问题可通过Hough变换、Kalman滤波及最小二乘法等求解，考虑到计 算的复杂度与鲁棒性，本实施方式中选择最小二乘法求解椭圆拟合问题。
[0048]【具体实施方式】五:本实施方式是对【具体实施方式】四所述的航天器相对位置增量定 轨方法作进一步说明，本实施方式中，采用拉格朗日乘子法求解拟合矩阵的特征值的具体 方法为：
[0049]将公式(8)改写成矩阵形式：
[0053]采用拉格朗日乘子求解式(9)，则有：
[0055] 上式中W = VTV记为散布矩阵，采用拉格朗日乘子求解矩阵W特征值(λι，υ〇,其中i =1，2, . . .，6，M为特征值，U为特征向量，最小的正实数解即为轨道拟合问题的解，即获得 参数U的估计值。
[0056] Fitzgibbon(菲茨吉本)给出了基于最小二乘原理的轨道拟合问题的解，在此基础 上，Hal ir (海立)进一步给出了快速非迭代的椭圆拟合算法，Nikolai (尼古拉)在Matlab的 CodeExchange (程序交流论坛）中公布了稳定可靠的拟合程序。上述对式(9)的求解就基于 Nikolai的程序。将式(6)与椭圆几何参数模型相比较，从而可根据计算得到的参数A获得图 2中参数xo，yo，a，lx5：a。
【具体实施方式】 [0057] 六:本实施方式是对五所述的航天器相对位置增量定 轨方法作进一步说明，本实施方式中，根据拟合矩阵的特征值判断航天器的轨道类型的具 体方法为：
[0058]在无噪声的情况下参数A满足如下条件：
[0060]考虑噪声的影响，设置门限eps用以区分圆轨道情况，
[0061] 将参数A，B，C代入下式，
[0062] (A/C-1)2+B2<eps (22)
[0063] 若结果满足式(22)，则轨道为圆轨道，
[0064] 若结果不满足式(22)，则轨道为椭圆轨道。
[0065] 本实施方式中，由于圆轨道可视为偏心率为零时的一种特殊的椭圆轨道，因此才 有式(21)的条件。
【具体实施方式】 [0066] 七:本实施方式是对六所述的航天器相对位置增量定 轨方法作进一步说明，本实施方式中，航天器轨道为椭圆轨道时，航天器初轨确定的具体方 法为：
[0067] 根据以下椭圆的几何参数方程确定椭圆在参考系Pxyz下的几何参数α、χο、y〇、a和 b：
[0070] 利用下式获得偏心率e:
[0072] 设参考系0x，y，z，以椭圆轨道中心为原点，z'轴垂直于轨道平面，x'，y'轴分别与椭 圆长轴方向和短轴方向同向，将椭圆轨道上各点坐标转换到参考系OxW下，满足右手定 则，有：
[0073] pk = Qxx' (xk+Sxx1) (14)
[0074] 其中，讲为航天器相对于参考系0x，y，z，原点的位置矢量，Q xx，为坐标转换矩阵，Sxx， 为坐标平移向量，解得：
[0076]设地球位于椭圆轨道的右焦点处，根据位置矢量Pk在参考系Ox，y，Z，下各方向分量 获得轨道各点的真近点角Qk :
[0078]将真近点角0k代入到开普勒时间方程中，获得平近点角Me(k):
[0080]根据平近点角与时间的关系
获得相应位置相对于近地点时刻的 时间t ' k，其中以近地点时刻为时间起点，若t ' k与测量量记录的时间tk具有相同的变化趋 势，则地球位于椭圆轨道的右焦点处，且仅有常值偏差;若t 'k与tk变化趋势不符，则地球位 于椭圆轨道的左焦点处，
[0081 ]根据地球位于椭圆轨道的位置获得近地点幅角ω和真近点角θ0:
[0084]
,.q〇 = XQsinc〇-y()cosc〇，（pQ，q())为航天器 初始位置在近焦点坐标系0 ' xyZ下的坐标；
[0085] 即能够获得航天器在地心J2000参考系下相对地心的绝对位置矢量rk，
[0086] + (2〇)
[0087]其中，化e为坐标转换矩阵的转置矩阵，为坐标平移向量，
【具体实施方式】 [0089] 八:本实施方式是对六所述的航天器相对位置增量定 轨方法作进一步说明，本实施方式中，航天器轨道为圆轨道时，航天器初轨确定的具体方法 为：
[0090] 利用下式确定圆轨道的半径Rc及中心位置坐标0c:
[0092] 利用圆轨道在平面参考系下的中心位置坐标0。，确定航天器各点在地心J2000参 考系下航天器的绝对位置矢量rk:
[0093] 24
[0094] 根据角动量矢量h获得位置矢量V k的速度信息Vk:
[0096]其中，μ为中心天体的引力常数。
[0097] 本实施方式中所述的引力常数μ，对于地球一般取y=398600km3/s2。本实施方式对 圆轨道确定了轨道各点的绝对位置及速度信息，即完成轨道的确定。
[0098]本发明讨论了根据相对位置矢量确定航天器椭圆、圆轨道的方法。下面对上述定 轨算法分别采用椭圆轨道和圆轨道进行检验。设观测的相对位置矢量的信噪比为20dB，每 15s获得一次观测信息，总观测数据为200时，分别进行1000次蒙特卡洛仿真，结果如表1和 表2所示。
[0099]表1 ST5-A卫星椭圆轨道定轨结果
[0101] 表2 PICARD卫星圆轨道定轨结果
[0103]仿真结果验证了本发明所述方法的可行性。
【主权项】
1. 航天器相对位置增量定轨方法，其特征在于，该方法为： 首先，确定航天器在Ξ维空间中运动的轨道平面； 然后，将Ξ维相对位置矢量转换为二维平面内航天器相对位置坐标，在该二维平面内 采用最小二乘原理进行航天器运动轨道形状的拟合，并采用拉格朗日乘子法求解拟合矩阵 的特征值； 最后，根据拟合矩阵的特征值判断航天器的轨道类型，若航天器轨道为楠圆轨道，则中 屯、天体位于楠圆的一个焦点处，结合开普勒时间方程与测量相对应的时间序列确定中屯、天 体的位置;若航天器轨道为圆轨道，则中屯、天体在轨道的中屯、；根据中屯、天体的位置，获得 航天器在惯性参考系下的绝对位置信息，实现航天器的初轨确定。2. 根据权利要求1所述的航天器相对位置增量定轨方法，其特征在于，确定航天器在Ξ 维空间中运动的轨道平面的具体方法为： 测量航天器在地屯、惯性坐标系下历元时刻与tw之间航天器移动的相对位置矢量A r (tJ-l，tJ)，并根据该相对位置矢量Δr(tw，tJ)获得航天器的初始位置r(to)，其中j = l，2， 3，. . . s 根据相对位置矢量Δ ;r(tj-i，tj)，采用最小二乘原理确定单位角动量矢量h;Cl) 其中，hx、hy、hz分别为单位角动量矢量h在空间Ξ维参考系中的X轴、y轴、z轴的分量，Δ Γχ(t j-i，t j)、Δ。(t j-i，t j)、Δ。(t j-i，t j)分别为相对位置矢量 Δ r (t j-i，t j)在；维参考系中 X轴、y轴、z轴的分量， 根据下式利用单位角动量矢量h获得轨道倾角i和升交点赤经Ω，、'3'/ 利用轨道倾角i和升交点赤经Ω确定Ξ维轨道平面坐标。3. 根据权利要求2所述的航天器相对位置增量定轨方法，其特征在于，将Ξ维相对位置 矢量转换为二维平面内航天器相对位置坐标的具体方法为： 设参考系PxyzW航天器的初始位置为原点，Z轴垂直于轨道平面，x，y轴由地球J2000参 考系根据坐标转换矩阵Qxx获得； 根据轨道倾角i和升交点赤经Ω获得航天器相对位置在参考系Pxyz中的矢量Xk: xk = Qxxmk (4) 式中，为k时刻航天器位置相对于初始位置的矢量，Qxx为坐标转换矩 阵，具体公式如下：(5) 利用矢量a在x，y轴上的投影获得航天器各点在二维平面内的位置坐标(xk，yk)。4. 根据权利要求1所述的航天器相对位置增量定轨方法，其特征在于，采用最小二乘原 理进行航天器运动轨道形状拟合的具体方法为： 设平面内任意轨道的表达式为：(飯 其中，U=[A，B，C，D，E，M]，X=[χ2，巧，y2，χ，y，l]，A、B、C、D、E、M均为多项式各项的系数， X和y分别为二维平面参考系Pxyz中轨道上各点坐标参数， 根据最小二乘原理将公式(6)转换为W下形式：5. 根据权利要求4所述的航天器相对位置增量定轨方法，其特征在于，采用拉格朗日乘 子法求解拟合矩阵的特征值的具体方法为： 将公式(8)改写成矩阵形式：上式中W=vTy记为散布矩阵，采用拉格朗日乘子求解矩阵W特征值(λι，υ〇,其中i = l， 2,...，6，λι为特征值，Ui为特征向量，最小的正实数解即为轨道拟合问题的解，即获得参数U 的估计值。6. 根据权利要求5所述的航天器相对位置增量定轨方法，其特征在于，根据拟合矩阵的 特征值判断航天器的轨道类型的具体方法为： 在无噪声的情况下参数A满足如下条件：^21 ) 考虑噪声的影响，设置口限eps用W区分圆轨道情况， 将参数A，B，C代入下式， (A/C-1)2+B2<eps (22) 若结果满足式(22)，则轨道为圆轨道， 若结果不满足式(22)，则轨道为楠圆轨道。7.根据权利要求6所述的航天器相对位置增量定轨方法，其特征在于，航天器轨道为楠 圆轨道时，航天器初轨确定的具体方法为： 根据W下楠圆的几何参数方程确定楠圆在参考系Pxyz下的几何参数〇、^〇、7〇、3和6:设参考系OxVz'W楠圆轨道中屯、为原点，Z'轴垂直于轨道平面，x'，y'轴分别与楠圆长 轴方向和短轴方向同向，将楠圆轨道上各点坐标转换到参考系Ox'yV下，满足右手定则，有： pik = Qxx'（Xk+Sxx'） （14) 其中，Pk为航天器相对于参考系Ox'yV原点的位置矢量，Qxx'为坐标转换矩阵，Sxx'为坐标 平移向量，解得：设地球位于楠圆轨道的右焦点处，根据位置矢量Pk在参考系Ox'yV下各方向分量获得轨 道各点的真近点角9k:根据平近点角与时间的关系= 获得相应位置相对于近地点时刻的时间 t'k，其中W近地点时刻为时间起点，若t'k与测量量记录的时间tk具有相同的变化趋势，贝U 地球位于楠圆轨道的右焦点处，且仅有常值偏差;若t ' k与tk变化趋势不符，则地球位于楠圆 轨道的左焦点处， 根据地球位于楠圆轨道的位置获得近地点幅角ω和真近点角θ〇:其中，扔而.跑S货一於础η (0-如2护，qo = XOSin ω -yocos ω，（ΡΟ, qo)为航天器初始位 置在近焦点坐标系0胃下的坐标； 即能够获得航天器在地屯、J2000参考系下相对地屯、的绝对位置矢量η，(如) 其中，溪;:为坐标转换矩阵谷的转置矩阵，式,Ε为坐标平移向量，8.根据权利要求6所述的航天器相对位置增量定轨方法，其特征在于，航天器轨道为圆 轨道时，航天器初轨确定的具体方法为： 利用下式确定圆轨道的半径Rc及中屯、位置坐标Oc:(23》 利用圆轨道在平面参考系下的中屯、位置坐标Oc，确定圆轨道时航天器各点在地屯、J2000 参考系下航天器的绝对位置矢量k:(24) 根据角动量矢量h获得位置矢量r/ k的速度信息vk:as) 其中，μ为中屯、天体的引力常数。
【文档编号】G01C21/24GK106092099SQ201610388810
【公开日】2016年11月9日
【申请日】2016年6月2日 公开号201610388810.8, CN 106092099 A, CN 106092099A, CN 201610388810, CN-A-106092099, CN106092099 A, CN106092099A, CN201610388810, CN201610388810.8
【发明人】徐国栋, 宋佳凝, 姜昆, 董立珉, 徐振东, 张兆祥
【申请人】哈尔滨工业大学