一种地月空间航天器的天基测定轨方法

摘要

本发明属于航天器轨道判定领域，尤其涉及一种地月空间航天器的天基测定轨方法。该方法包括：步骤1、获取待定轨航天器在初始历元时的初始状态以及初始误差；步骤2、在初始历元的下一历元时，获取低轨卫星与待定轨航天器之间的距离，根据该距离对所述初始状态以及所述初始误差进行扩展卡尔曼滤波处理，得到待定轨航天器在下一历元的状态以及误差；步骤3、以所述下一历元作为初始历元，以所述下一历元的状态作为初始状态，以所述下一历元的误差作为初始误差，重复执行步骤2，得到多组状态以及误差，将所有历元的状态以及误差进行整合得到待定轨航天器的轨迹。通过该方法能够达到地月空间航天器的快速和高精度定轨的效果。

权利要求

1.一种地月空间航天器的天基测定轨方法，其特征在于，包括：

步骤1、获取待定轨航天器在初始历元时的初始状态以及初始误差；

步骤2、在初始历元的下一历元时，获取通过GNSS确定的低轨卫星的位置和速度，获取通过星上测量设备确定的低轨卫星与待定轨航天器之间的距离，根据低轨卫星的位置和速度和该距离对所述初始状态以及所述初始误差进行扩展卡尔曼滤波处理，得到待定轨航天器在下一历元的状态以及误差；

步骤3、以所述下一历元作为初始历元，以所述下一历元的状态作为初始状态，以所述下一历元的误差作为初始误差，重复执行步骤2，得到多组状态以及误差，将所有历元的状态以及误差进行整合得到待定轨航天器的轨迹；

所述扩展卡尔曼滤波处理具体包括：

时间更新以及测量更新；

所述时间更新的具体过程为：

基于动力学模型，根据第i个历元ti的状态Xi对第i+1个历元ti+1的状态Xi+1进行预测；

所述测量更新的具体过程为：

利用第i+1个历元ti+1下的低轨卫星与所述待定轨航天器之间的距离以及该历元对应的通过GNSS确定的低轨卫星的位置和速度，对第i+1个历元ti+1的状态以及第i+1个历元ti+1的误差进行校正；

所述根据第i个历元ti的状态Xi对第i+1个历元ti+1的状态Xi+1进行预测具体为：其中，Φi+1|i＝ΦY(ti+1,ti)，Γi+1|i＝ΓY(ti+1,ti)，ΦY(ti+1,ti)为历元ti+1相对于历元ti的状态转移矩阵,ΓY(ti+1,ti)为历元ti的加速度过程噪声对历元ti+1的位置速度状态的转移矩阵，ui＝u(ti)＝(uxi,uyi,T T T Tuzi) 为加速度过程噪声，Y＝(r ,v ) ，r为航天器的位置矢量，v为航天器的速度矢量，为加速度过程噪声ui的协方差矩阵，其中的 和ρuz分别为在地心惯性系x、y和z坐标轴方向的加速度过程噪声标准差,X(ti+1，Xi)为数值积分后历元ti+1的未增加过程噪声影响的状态预测值，Pi为历元ti的状态协方差矩阵， ri为历元ti的待定轨航天器的位置矢量，vi为历元ti的待定轨航天器的速度矢量，Xi+1|i为根据ti历元对应的定轨结果Xi计算出的ti+1历元对应的状态预测值，Pi+1|i为ti+1历元的状态预测值Xi+1|i所对应的协方差矩阵，Y(ti)为从ti历元到ti+1历元的轨道预报过程中ti对应的位置速度状态矢量，Y(ti+1)为从ti历元到ti+1历元的轨道预报过程中ti+1历元对应的位置速度状态矢量；

所述利用第i+1个历元ti+1下的低轨卫星与所述待定轨航天器之间的距离，对第i+1个历元ti+1的状态以及第i+1个历元ti+1的误差进行校正的具体过程为：Xi+1＝Xi+1|i+Ki+1(zi+1‑h(Xi+1|i))

其中，Pi+1为ti+1历元定轨结果状态Xi+1的协方差矩阵；Ki+1是卡尔曼增益；Pi+1|i为ti+1历元的状态预测值Xi+1|i所对应的协方差矩阵； 是测量噪声标准差；Hi+1是星间距离测量值对状态的设计矩阵；I是单位矩阵，ti+1为时间，zi+1为星上测量设备在

ti+1历元获得的距离测量值,r代表待定轨航天器，s代表用于测量的低轨卫星， 为低轨卫星和待定轨航天器低轨卫星和待定轨航天器之间的几何距离，ε(ti+1)为测量噪声,h(Xi+1|i)表示通过状态预测值Xi+1|i和GNSS确定的低轨卫星的位置和速度计算出的ti+1历元的低轨卫星和待定轨航天器之间的距离预测值，h(X(ti+1),ti+1)表示通过ti+1历元的状态X(ti+1)计算出低轨卫星和待估航天器间距离的函数表达式；

说明书

一种地月空间航天器的天基测定轨方法

技术领域

[0001] 本发明属于航天器轨道判定领域，尤其涉及一种地月空间航天器的天基测定轨方法。

背景技术

[0002] 一般航天器的定轨方式主要有两种，一种是依靠全球导航卫星系统(GNSS)，另一种是依靠地面站测量。在距离地球几万公里以内的空间范围内，基于GNSS测量的定轨方式已经比较成熟，利用GNSS实现航天器的定轨是最常用的方法。绕地运行的航天器通过GNSS信号接收机接收GNSS星座发出的测量信号获得测量数据，然后对测量数据进行处理，计算出位置速度状态信息进而实现定轨。由于GNSS测量信号随着距离的增加而衰减，当航天器距离地球太远时无法利用GNSS实现定轨，所以远距离航天器的定轨一般采用其它方法。

[0003] 处于地月转移轨道、环月轨道及太阳系其他行星转移轨道的航天器距地球远，主要依靠地面站测量实现定轨，如果对定轨的精度要求不高，还可以采用天文导航的方式，通过测量距航天器较近的天体与遥远恒星的星光角距，利用星光角距变化与航天器位置变化的关系，实现自主定轨。

[0004] 在以往的月球探测任务和行星探测任务中，通过地面站测量实施定轨的方式已经被广泛使用，其可行性和可靠性已经被验证，但这种方式具有一定的局限性。航天器远离地球，速度较慢，地面站随地球自转同步旋转，当地面站被地球遮挡，位于相对于航天器的地球背面时，无法进行定轨，且持续时间可达十小时以上。另一方面，航天器的低速和地球的低转速导致航天器相对于地面站的距离和速度变化较为缓慢，因此测量敏感性较低，对于实时定轨方法，需要很长的时间才能收敛，定轨精度也不高。通过增加地面站的数量和覆盖范围虽然可以解决这些问题，但是为了满足未来地月空间的大量探索任务，需要建设覆盖全球的地面站，这显然具有较大的阻力，而且花费巨大，所以需要更好的解决方法。

发明内容

[0005] 本发明所要解决的技术问题是提供一种地月空间航天器的天基测定轨方法。

[0006] 本发明解决上述技术问题的技术方案如下：一种地月空间航天器的天基测定轨方法，包括：

[0007] 步骤1、获取待定轨航天器在初始历元时的初始状态以及初始误差；

[0008] 步骤2、在初始历元的下一历元时，获取低轨卫星与待定轨航天器之间的距离，根据该距离对所述初始状态以及所述初始误差进行扩展卡尔曼滤波处理，得到待定轨航天器在下一历元的状态以及误差；

[0009] 步骤3、以所述下一历元作为初始历元，以所述下一历元的状态作为初始状态，以所述下一历元的误差作为初始误差，重复执行步骤2，得到多组状态以及误差，将所有历元的状态以及误差进行整合得到待定轨航天器的轨迹。

[0010] 本发明的有益效果是：采用基于低轨卫星的天基测量方式，实现地月空间航天器的快速和高精度定轨，解决传统地基测量信号容易遮挡、定轨收敛慢的问题，定轨精度高于传统定轨方案。

[0011] 在上述技术方案的基础上，本发明还可以做如下改进。

[0012] 进一步，所述扩展卡尔曼滤波处理具体包括：

[0013] 时间更新以及测量更新。

[0014] 进一步，所述时间更新的具体过程为：

[0015] 基于动力学模型，根据第i个历元的状态对第i+1个历元的状态进行预测。

[0016] 进一步，所述测量更新的具体过程为：

[0017] 利用第i+1个历元下的低轨卫星与所述待定轨航天器之间的距离，对第i+1个历元的状态以及第i+1个历元的误差进行校正。

[0018] 进一步，所述根据第i个历元ti的状态Xi对第i+1个历元ti+1的状态Xi+1进行预测具体为：

[0019]

[0020] 其中，Φi+1|i＝ΦY(ti+1,ti)，Γi+1|i＝ΓY(ti+1,ti)，ΦY(ti+1,ti)为历元ti+1相对于历元ti的状态转移矩阵,ΓY(ti+1,ti)
为历元ti的加速度过程噪声对历元ti+1的位置速度状态的转移矩阵，ui＝u(ti)＝(uxi,uyi,T T T T
uzi) 为加速度过程噪声，Y＝(r ,v ) ，r为航天器的位置矢量，v为航天器的速度矢量，为加速度过程噪声ui的协方差矩阵，其中的ρ为加速度过程噪声标
准差,X(ti+1，Xi)为数值积分后历元ti+1的状态预测值，，Pi为历元ti的状态协方差矩阵，ri为历元ti的待定轨航天器的位置矢量，vi为历元ti的待定轨航天器的速度矢
量。

[0021] 进一步，所述利用第i+1个历元ti+1下的低轨卫星与所述待定轨航天器之间的距离，对第i+1个历元ti+1的状态以及第i+1个历元ti+1的误差进行校正的具体过程为：

[0022]

[0023] Xi+1＝Xi+1|i+Ki+1(zi+1‑h(Xi+1|i))

[0024]

[0025] 其中，Ki+1是卡尔曼增益； 是测量噪声标准差；Hi+1是星间距离测量值对状态的设计矩阵；I是单位矩阵， ti为时间，z(ti)为星间链路距离测量值，r代表待定轨航天器，s代表用于测量的低轨卫星， 为低轨卫星和待定轨航天器之间的几何距离，ε(ti)为测量噪声,h(Xi+1|i)表示通过状态预测值Xi+1|i计算出的ti+1历元的LEO导航星和待定轨航天器之间的距离预测值，

[0026]

[0027] 本专利受到中国科学院战略性先导科技专项资助，专项编号：XDA30010000。

附图说明

[0028] 图1为本发明一种地月空间航天器的天基测定轨方法的实施例提供的流程示意图；

[0029] 图2为本发明一种地月空间航天器的天基测定轨方法的实施例提供的LEO卫星天基测量的地月航天器定轨方案示意图；

[0030] 图3为本发明一种地月空间航天器的天基测定轨方法的实施例提供的DRO轨道地月旋转系示意图；

[0031] 图4为本发明一种地月空间航天器的天基测定轨方法的实施例提供的DRO轨道地心惯性系示意图；

[0032] 图5为本发明一种地月空间航天器的天基测定轨方法的实施例提供的RO轨道示意图；

[0033] 图6为本发明一种地月空间航天器的天基测定轨方法的实施例提供的DRO入轨转移轨道示意图；

[0034] 图7为本发明一种地月空间航天器的天基测定轨方法的实施例提供的站星测量精度与站星距离的关系示意图；

[0035] 图8为本发明一种地月空间航天器的天基测定轨方法的实施例提供的仅使用地基测量的低能耗DRO转移轨道定轨结果示意图；

[0036] 图9为本发明一种地月空间航天器的天基测定轨方法的实施例提供的仅使用LEO测量的低能耗DRO转移轨道定轨结果示意图；

[0037] 图10为本发明一种地月空间航天器的天基测定轨方法的实施例提供的仅使用地基测量的DRO停泊轨道定轨结果示意图；

[0038] 图11为本发明一种地月空间航天器的天基测定轨方法的实施例提供的仅使用LEO测量的DRO停泊轨道定轨结果示意图；

[0039] 图12为本发明一种地月空间航天器的天基测定轨方法的实施例提供的仿真使用的RO轨道及航天器的初始相位示意图；

[0040] 图13为本发明一种地月空间航天器的天基测定轨方法的实施例提供的仅使用地基测量的RO轨道定轨结果示意图；

[0041] 图14为本发明一种地月空间航天器的天基测定轨方法的实施例提供的仅使用LEO测量的RO轨道定轨结果示意图；

[0042] 图15为本发明一种地月空间航天器的天基测定轨方法的实施例提供的扩展卡尔曼滤波实现定轨示意图。

具体实施方式

[0043] 以下对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

[0044] 如图1所示，一种地月空间航天器的天基测定轨方法，包括：

[0045] 步骤1，获取初始历元t0下，待定轨航天器的初始状态X0以及初始误差P0；

[0046] 步骤2，通过下一历元t1下的低轨卫星与所述待定轨航天器之间的距离，对所述初始状态X0以及所述初始误差P0进行扩展卡尔曼滤波处理，得到历元t1的定轨结果，所述定轨结果包括：状态X1以及误差P1；

[0047] 步骤3，重复步骤2，得到n个定轨结果，将n个所述定轨结果进行整合得到待定轨航天器的轨迹。

[0048] 在一些可能的实施方式中，采用基于低轨卫星的天基测量方式，实现地月空间航天器的快速和高精度定轨，解决传统地基测量信号容易遮挡、定轨收敛慢的问题，定轨精度高于传统定轨方案。

[0049] 需要说明的是，“测定轨”指的就是轨迹确定，而且是地月空间航天器的轨道。“天基”指的是相对于“地基”而言，地基定轨方法中使用的是位置已知的地面站。本专利中，定轨方法中使用的是位置速度状态已知的低轨(LEO)卫星，不是地面站；因为低轨卫星是运行在天上的，绕地球运动，所以叫做“天基”。本专利中提及的历元为一种航天领域的时间计量单位，且第一历元即对应的是历元t1，第二历元对应的是历元t2，以此例推，同理，第一历元的状态即为X1，第i历元的状态为Xi，以此类推，同样的也适用于误差。

[0050] 本专利的技术方案如图2所示，LEO卫星配置GNSS接收机接收GNSS导航卫星信号，在轨可实现m级定位精度和30ns级授时精度。LEO卫星和地月空间航天器均配备无线电收发设备，使得LEO卫星能与地月航天器建立链路，完成测距测速。地月航天器接收到测距信号和LEO卫星的轨道状态后，可在轨自行计算轨道信息。

[0051] 本专利提出一种利用低轨卫星测量数据和扩展卡尔曼滤波算法对地月空间航天器实施定轨的方法。该方法有两个关键点：一个是低轨卫星的状态和轨道已知时，仅使用距离测量数据对地月空间航天器进行定轨的原理；另一个是基于测量数据的扩展卡尔曼滤波定轨算法。下面着重介绍这两方面：

[0052] 基于低轨卫星距离测量数据实现定轨的基本原理

[0053] 利用GNSS测量实现低轨卫星定轨的方法已经比较成熟，精度较高，所以低轨卫星的轨道已确定。一颗低轨卫星与地月空间航天器之间建立测量链路(LEO星载设备发送测距信号，该测距信号含有发送时刻的时间信息，地月空间航天器星载接收机接收到测距信号，提取出发射时间，然后计算出发射与接收的时间差，将其乘以光速获得距离，该技术为成熟技术)，获得两者之间的距离随时间的变化序列(这里提到的序列是一个理论上的描述，每个历元一个距离数据，不同历元的数据合在一起称作序列。测量数据用于定轨过程，可以不记录，数据获取后直接发送给定轨程序；即便需要记录，简单地写入文件或内存即可，不需要特别的算法。)。与地面站测量不同，因为低轨卫星的运行速度非常快，绕地周期短，所以星间测量被地球遮挡的时间也短，星间距离变化相对较快，具有更好的测量敏感性。地月空间航天器初始状态不同时，轨道也不一样，测量链路获得的测量序列也不一样，因此当测量序列(不同时间的距离测量数据组成的序列)和初始状态已知时，可以唯一地确定地月空间航天器的轨道。(这一段文字的第一句中说：“低轨卫星的轨道由GNSS确定”，因此低轨卫星的位置速度状态是已知的，是前提条件，相当于导航星；需要定轨的是地月空间的其他航天器。这一句之前的文字是对定轨原理的描述，即为什么可以实现定轨，后续的才是具体的实现方法。)距离测量数据的时间序列与动力学模型(高精度轨道动力学模型，计算使用的模型是实测地球月球引力模型，在后面的实施例的“力学模型及积分器设置”表给出了具体的名称。)进行匹配(匹配的过程在扩展卡尔曼滤波中实现，输出结果是不同历元的航天器位置速度状态组成的轨道。扩展卡尔曼滤波的具体描述如图15所示，图中ti‑1为第i‑1个历元，ti为第i个历元，ti+1为第i+1个历元，虚线代表真实轨道，空心方块代表定轨结果，空心圆圈代表测量值，实心圆圈代表预测值，图中，从第i‑1历元的定轨结果到第i个历元的预测值之间进行时间更新，第i个历元的预测值到第i个历元的定轨结果之间进行测量更新，第i个历元的定轨结果到第i+1个历元的预测值之间进行时间更新，第i+1个历元的预测值到第i+1个历元的定轨结果之间进行测量更新，第i+1个历元的定轨结果到下一个历元的预测值之间进行时间更新，以此往复。)，使用扩展卡尔曼滤波的方法计算出最优匹配轨道，从而实现地月空间航天器的定轨。

[0054] 扩展卡尔曼滤波定轨方法

[0055] X＝(rT,vT)T为待估量，表示航天器的位置速度状态，其误差由状态协方差矩阵P表示。位置矢量r和速度矢量v如式所示，x、y和z为航天器位置在J2000地心惯性系下的坐标分量，vx、vy和vz为速度的坐标分量。协方差矩阵P如式所示，对角线元素为位置和速度的方差，非对角线元素为协方差。

[0056]

[0057]

[0058] 初始化：t0为初始历元，对应的初始状态及协方差矩阵为X0和P0，如式和式所示。

[0059] X0＝(x0 y0 z0 vx0 vy0 vz0)T

[0060]

[0061] 定轨流程：利用t1历元测量数据对X0和P0进行一次扩展卡尔曼滤波处理，得到t1历元的定轨结果X1和P1，然后利用t2历元测量数据对X1和P1再次进行一次扩展卡尔曼滤波处理，得到t2历元的定轨结果X2和P2，重复上述过程，计算出X3与P3、X4与P4、X5与P5等，X1、X2、X3等状态组成最终的定轨结果。

[0062] 基于扩展卡尔曼滤波的定轨方法是实时定轨方法，所以理论上只要有测量数据，此过程可以一直持续下去。工程任务中，根据实际情况选择结束条件。

[0063] 每一次扩展卡尔曼滤波处理有两个步骤：第一个步骤叫做“时间更新”，即利用前一历元状态(状态和状态协方差矩阵一起使用)使用高精度动力学模型对下一历元状态进行预测；第二个步骤叫“测量更新”，利用下一历元的测量数据对获得的预测状态进行校正，计算出下一历元的状态定轨结果。在“时间更新”阶段需用到轨道动力学模型，在“测量更新”阶段需要用到测量数据对待估状态的偏导数，因此先详细先介绍轨道动力学模型、测量模型和测量模型的线性化。

[0064] 1.轨道动力学模型

[0065] 地月系统内无机动自由航行的航天器的运动规律可以用一阶常微分方程表示，如式所示：

[0066]

[0067] 其中，r为航天器的位置矢量，v为航天器的速度矢量，rm为月球的位置矢量，rs为太阳的位置矢量，月球和太阳位置矢量由星历求得。μe为地球引力常数，μm为月球引力常数；μs为太阳引力常数；aa为大气阻力加速度；as为太阳辐射压力加速度；ae和am分别为地球和月球非球形摄动加速度。式中的向量均表示在J2000地心惯性坐标系中。令状态向量为Y＝T T T(r ,v) ，则上式可以简写为：

[0068]

[0069] 结合力学模型，使用数值积分的方式可以计算出轨道状态。任意历元t的状态相对于初始历元t0的状态转移矩阵定义为：

[0070]

[0071] 使用高精度轨道动力学模型的状态预测过程的线性化数学表达，表示由上一历元状态转移到下一历元状态的过程，状态转移矩阵乘以上一历元的状态等于下一历元的状态，即Y(t)＝ΦY(t,t0)Y(t0)。

[0072] 由于动力学模型有误差，因此引入加速度过程噪声u(t)＝(ux uy uz)T模拟未建模的动力学模型误差。每个历元的u(t)都是不同的，但是服从均值为零的正态分布。t0历元的加速度过程噪声对t历元状态的影响表示为转移矩阵ΓY：

[0073]

[0074] 将其进一步展开可以得到：

[0075]

[0076] 式中的Δt为数值积分步长，I3×3为单位矩阵。

[0077] 2.测量模型及其线性化

[0078] 对于单个星间链路(两个卫星)，通过测量获得低轨卫星和航天器之间的距离，其模型表示如下：

[0079]

[0080] 其中，t为时间，z(t)为星间链路距离测量值，r代表航天器，s代表用于测量的低轨卫星， 为低轨卫星和航天器之间的几何距离，ε(t)为测量噪声。

[0081] 在扩展卡尔曼滤波的“测量更新”定轨流程中，需要用到距离对状态的偏导数，因此需要对测量进行线性化，对测量模型求偏导得到：

[0082]

[0083] 其中，r(t)为航天器质心位置； 为航天器指向低轨卫星的单位向量(Light of Sight,LoS)。

[0084] 3.扩展卡尔曼滤波的“时间更新”

[0085] “时间更新”步骤即由前一个历元ti的状态Xi和协方差矩阵Pi，通过高精度轨道动力学模型，利用数值积分的方法计算推出下一个历元ti+1的状态预测值Xi+1|i及其对应的协方差矩阵Pi+1|i，计算方式如下：

[0086] Xi+1|i＝X(ti+1,X(ti)＝Xi)+Γi+1|iui

[0087]

[0088] 轨道动力学模型如式所示，未建模的轨道动力学误差由加速度过程噪声ui表示。状态转移矩阵计算方式如下：

[0089] Φi+1|i＝ΦY(ti+1,ti)

[0090] 过程噪声转移矩阵表示为：

[0091] Γi+1|i＝ΓY(ti+1,ti)

[0092] 在公式中，Qc是加速度过程噪声ui的协方差矩阵如式所示，它是一个三乘三的对角矩阵。对角线上的值表示加速度噪声的标准差，是固定值，根据未建模的动力学模型误差设置。增加加速度过程噪声能够防止状态协方差矩阵收敛过快进而导致定轨结果恶化，同时允许后续的测量数据持续影响滤波估计。在定轨仿真测试中，采用试错法选择以获得最优定轨结果。

[0093]

[0094] 4.扩展卡尔曼滤波的“测量更新”

[0095] 扩展卡尔曼滤波的第二个步骤是“测量更新”，利用ti+1历元的距离测量数据zi+1对来自时间更新的状态预测值Xi+1|i及其对应的协方差矩阵Pi+1|i进行状态修正处理，计算出ti+1历元的定轨结果Xi+1和协防差矩阵Pi+1：

[0096]

[0097] Xi+1＝Xi+1|i+Ki+1(zi+1‑h(Xi+1|i))

[0098]

[0099] 式中Ki+1是卡尔曼增益； 是测量噪声标准差；Hi+1是星间距离测量值对状态的设计矩阵(偏导数)；I是单位矩阵。设计矩阵Hi+1的计算方法如下：

[0100]

[0101] 状态协方差矩阵Pi+1表示估计出的状态的精度，也可以用来评估定轨过程的收敛性。

[0102] 实施例1，下面以三种地月空间典型轨道为例，利用本专利的方法对位于这些轨道上的航天器进行定轨计算，并对结果进行评估。

[0103] 一、三种轨道基本介绍

[0104] 第一种轨道是DRO轨道。航天器离地球数万公里之外，位于月球附近，在地月旋转系下绕月球运动，绕月球旋转的方向与月球公转的方向相反，在地心惯性系下绕地球旋转的方向与月球公转方向相同，而且月球的公转周期与航天器绕月球的周期成整数比，该比例为共振比。此时航天器所在的轨道为远距离逆行轨道(DRO)，典型的共振比为2:1、3:1和4:1。本专利仿真用到的是周期比2:1的DRO，其在地月旋转系和地心惯性系下的表示如图1所示，图中的DU表示地月平均距离，图3是DRO(2:1)在地月旋转系中的表示，从图中可以看出轨道距离月球相对较近，距离地球相对较远，形状近似椭圆。图4是DRO(2:1)在地心惯性系中的表示，轨道形状近似椭圆，但是地球在轨道的中心，不在焦点上。

[0105] 第二种是RO轨道。大椭圆轨道指偏心率很大，近地点较低，远地点很高的椭圆轨道。这种轨道数量众多，其中一类具有特殊性质，即月球公转周期与其轨道周期成整数比，这个性质使该类轨道受到关注。图5所示的是一种稳定的大椭圆轨道，左图为地月旋转系，右图为地心惯性系，共振比为3:1，本专利仿真使用此RO轨道。图中的DU表示地月平均距离。

[0106] 第三种是DRO入轨转移轨道。航天器从地球附近出发转移至DRO不采用直接转移方式，而是使用低能耗转移轨道，这种方式不需要消耗太多燃料，但代价是很长的转移时间。低能耗转移轨道如图6所示，位于该轨道上的航天器在地球、太阳和月球三体作用下，逐渐运行到DRO轨道插入点，经过机动进入DRO轨道，在此过程中，航天器距离地球最远可达一百万公里以上。

[0107] 二、链路和轨道设置。

[0108] 星间测量的噪声设为0.5m，系统差为5m。地面站位于北京密云，测距误差随着站星距离的变化关系，如图7所示。拟合为星地距离的函数，y为测量噪声(单位为m)，x为站星距离(单位为km)，可见星地距离50万km，测距精度优于1m，星地距离100万km，测距精度优于3m，星地距离150万km，测距精度为5m。星地测量系统差设为5m。

[0109]

[0110] DRO和RO轨道初始时间设置为2023年1月1日零时(UTC时间)，初始状态如表1所示(表示在J2000坐标系下)，力学模型如表2所示。LEO轨道是高度为500km的太阳同步轨道，初始轨道状态如表3所示(表示在J2000坐标系下)，力学模型如表2所示。

[0111] 表1 DRO和RO轨道的初始时刻状态(J2000坐标系)

[0112]

[0113] 表2力学模型及积分器设置

[0114]

[0115]

[0116] 表3 LEO轨道的初始时刻状态(J2000坐标系)

[0117] 参数 LEO半长轴/km 6878.137
偏心率 0
轨道倾角/° 97.4065
近地点幅角/° 0
升交点赤经/° 10.3886
真近点角/° 0
绕地/月周期 1.5小时

[0118] 假设低能耗DRO转移入轨轨道(如图3所示)的发射(星箭分离)时刻10Feb 202315:33:36.094，卫星进入大椭圆轨道，飞行一圈后，近地点机动时刻16Feb 202323:48:03.018(UTC)，机动脉冲为35.05m/s，深空机动时刻24Mar 202317:21:30.665(UTC)，深空机动脉冲
6.32m/s，到达时刻4Jun 202323:14:11.872(UTC)，入轨脉冲63.46m/s，飞行时间107.98天，卫星脉冲消耗104.83m/s，如表4所示。低能耗DRO转移入轨轨道的力学模型设置如表5所示，太阳与行星的实际位置根据JPL的星历文件DE430插值得到。

[0119] 表4低能耗DRO转移轨道机动

[0120]

[0121]

[0122] 表5低能耗DRO转移入轨轨道的力学模型

[0123]

[0124] 低轨卫星通过接收GNSS卫星实现定轨，位置精度优于1m，速度精度达到亚mm/s，仿真中低轨卫星的精度设置如表6所示。

[0125] 表6低轨卫星精度

[0126]

[0127] 基于低轨卫星的天基定轨参数设置如表7所示，力学模型如表2所示。

[0128] 表7基于低轨卫星的天基定轨参数设置

[0129]

[0130] 假设地面站位于北京密云，位置精确已知。地基定轨参数设置如表8所示，力学模型仍采用表2所列模型。

[0131] 表8基于地基的定轨参数设置

[0132]

[0133]

[0134] 定轨仿真结果。

[0135] 下面仿真对比基于低轨卫星的天基定轨结果和地基定轨结果。每种测量方式的结果分两种情形，假设测量是连续的(地球遮挡是需要考虑的)和假设每8小时弧段包含2小时测量。

[0136] 低能耗DRO转移轨道的定轨结果

[0137] 低能耗DRO转移轨道上的航天器在转移途中实施了3次机动，定轨过程过程模拟了含有误差的脉冲，误差为实际脉冲大小的20％。

[0138] 仅使用地基测量时，收敛的判定标准为位置残差小于10km，定轨结果如图8所示。

[0139] 仅使用LEO测量时，收敛的判定标准为位置残差小于1km，定轨结果如图9所示。

[0140] 定轨统计结果如表9所示。(MC)

[0141] 表9低能耗DRO转移轨道定轨结果

[0142]

[0143]

[0144] DRO轨道的定轨结果

[0145] DRO轨道定轨收敛的判定标准为位置残差小于100m。

[0146] 仅使用地基测量的定轨结果如图10所示。

[0147] 仅使用LEO测量的定轨结果如图11所示。

[0148] 定轨统计结果如表10所示。(MC)

[0149] 表10 DRO停泊轨道定轨结果

[0150]

[0151] RO轨道的定轨结果

[0152] RO轨道定轨收敛的判定标准为位置残差小于100m。

[0153] RO轨道的不同位置与地球间的距离相差较大，初始位置的选择可能会影响定轨的收敛时间，仿真过程中使用的RO轨道如图12所示。

[0154] 仅使用地基测量的定轨结果如图13所示。

[0155] 仅使用LEO测量的定轨结果如图14所示。

[0156] 定轨统计结果如表11所示。

[0157] 表11 RO轨道定轨结果

[0158]

[0159] 结果分析。

[0160] 定轨结果显示，仅使用地面站测量时，定轨过程收敛所需要的时间为几十甚至上百个小时，而且对速度脉冲误差十分敏感，当速度脉冲误差较大时，仅靠地面站测量，定轨过程很难收敛。基于低轨卫星持续测量进行定轨时，三种不同轨道上的航天器定轨需要的收敛时间都在两三个小时左右，相对于地面站来说有着很大改善。定轨精度方面，对于三个不同轨道的航天器，使用两种测量弧长策略的定轨结果均显示基于低轨卫星测量的定轨精度好于基于地面站测量的定轨精度，这在非周期DRO转移入轨轨道上尤为明显，定轨精度差距在一个数量级左右。从定轨结果来看，基于低轨卫星测量的定轨方法，成功解决了基于地面站测量进行定轨时地面站被长时间遮挡及测量敏感性较低导致的收敛时间较长和精度较差的问题，而且低轨卫星的部署比全球范围建设大量地面站更方便实现，大量的低轨卫星在未来有可能组成类似GPS的天基测控网以满足更多的地月空间探索任务。

[0161] 优选地，在上述任意实施例中，所述扩展卡尔曼滤波处理具体包括：

[0162] 时间更新以及测量更新。

[0163] 优选地，在上述任意实施例中，所述时间更新的具体过程为：

[0164] 基于动力学模型，根据历元ti的状态Xi对历元ti+1的状态Xi+1进行预测，其中，i为不大于n的自然数。

[0165] 优选地，在上述任意实施例中，所述测量更新的具体过程为：

[0166] 利用历元ti+1下的低轨卫星与所述待定轨航天器之间的距离，对状态Xi+1以及误差Pi+1进行校正。

[0167] 优选地，在上述任意实施例中，所述根据历元ti的状态Xi对历元ti+1的状态Xi+1进行预测具体为：

[0168]

[0169] 其中，Φi+1|i＝ΦY(ti+1,ti)，Γi+1|i＝ΓY(ti+1,ti)，ΦY(ti+1,ti)为历元ti+1相对于历元ti的状态转移矩阵,ΓY(ti+1,ti)为
历元ti的加速度过程噪声对历元ti+1的位置速度状态的转移矩阵，ui＝u(ti)＝(uxi,uyi,uziT T T T
) 为加速度过程噪声，Y＝(r ,v ) ，r为航天器的位置矢量，v为航天器的速度矢量，为加速度过程噪声ui的协方差矩阵，其中的ρ为加速度过程噪声标
准差。

[0170] 优选地，在上述任意实施例中，所述利用历元ti+1下的低轨卫星与所述待定轨航天器之间的距离，对状态Xi+1以及误差Pi+1进行校正的具体过程为：

[0171]

[0172] Xi+1＝Xi+1|i+Ki+1(zi+1‑h(Xi+1|i))

[0173]

[0174] 其中，Ki+1是卡尔曼增益； 是测量噪声标准差；Hi+1是星间距离测量值对状态的设计矩阵；I是单位矩阵， ti为时间，z(ti)为星间链路距离测量值，r代表待定轨航天器，s代表用于测量的低轨卫星， 为低轨卫星和待定轨航天器之间的几何距离，ε(ti)为测量噪声,

[0175]

[0176] Hi+1为设计矩阵。

[0177] 读者应理解，在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。

[0178] 在本申请所提供的几个实施例中，应该理解到，所揭露的装置和方法，可以通过其它的方式实现。例如，以上所描述的方法实施例仅仅是示意性的，例如，步骤的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个步骤可以结合或者可以集成到另一个步骤，或一些特征可以忽略，或不执行。

[0179] 上述方法如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分，或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read‑OnlyMemory)、随机存取存储器(RAM，RandomAccessMemory)、磁碟或者光盘等各种可以存储程序代码的介质。

[0180] 以上，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。