CN107065025A - 一种基于重力梯度不变量的轨道要素估计方法 - Google Patents

一种基于重力梯度不变量的轨道要素估计方法 Download PDF

Info

Publication number
CN107065025A
CN107065025A CN201710024671.5A CN201710024671A CN107065025A CN 107065025 A CN107065025 A CN 107065025A CN 201710024671 A CN201710024671 A CN 201710024671A CN 107065025 A CN107065025 A CN 107065025A
Authority
CN
China
Prior art keywords
mrow
msub
mtd
mfrac
mtr
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201710024671.5A
Other languages
English (en)
Other versions
CN107065025B (zh
Inventor
陈培
孙秀聪
张键
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN201710024671.5A priority Critical patent/CN107065025B/zh
Publication of CN107065025A publication Critical patent/CN107065025A/zh
Application granted granted Critical
Publication of CN107065025B publication Critical patent/CN107065025B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V7/00—Measuring gravitational fields or waves; Gravimetric prospecting or detecting

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
  • Navigation (AREA)

Abstract

一种基于重力场梯度不变量的轨道要素估计方法，其步骤如下：一：准备工作；二：理想重力梯度张量在东北天(East‑North‑Up，ENU)坐标系下的分解和特征值；三：求解测量历元J2模型重力场张量在ENU坐标系下的分解和特征值；四：利用J2重力梯度矩阵特征值求解各测量历元r和；五：利用和计算初始轨道要素；六：轨道要素平滑。通过以上步骤本发明创新的使用了卫星运行过程中所测量的重力梯度矩阵不变量信息，通过迭代求解出卫星到地心的几何距离与地心纬度，并以此作为观测量给出了以轨道根数中的半长轴、偏心率、轨道倾角、近地点俯角和真近点角为状态参数的测量方程；并创新地引入了轨道根数的摄动动力学方程，采用批处理最小二乘估计出初始历元的五个轨道根数；该方法具有自主程度高、抗干扰能力强、建设成本低等优点，而且在深空探测领域具有传统方法所不具有的优势。

Description

一种基于重力梯度不变量的轨道要素估计方法
技术领域
本发明提供一种基于重力场梯度不变量的轨道要素估计方法，它涉及一种对在轨卫星利用重力场梯度测量数据进行轨道要素估计的方法，属于导航技术领域。
背景技术
地球重力场是受地球重力作用的空间范围，而由于地球内部质量分布的不规则性，使得地球重力场不是按一个简单规律变化的力场。为了建立地球重力场模型，传统的物理大地测量学依据地球表面收集的观测数据形成和发展了一系列重力场理论和方法，如重力势、重力位函数、大地水准面等。由于观测手段的限制，这个阶段形成的地球重力场模型比较粗糙，并不能给出相对精确的地球重力场分布，使得利用地球重力场特性的应用受限。进入20世纪后半叶，利用新的测量技术(如电磁波测距、人造地球卫星定位***和甚长基线干涉测量(VLBI)等)和先进的重力测量仪器为建立更加精准的地球重力场模型提供了支持。如今，由于高精度的地球重力场模型的建立，利用重力梯度进行导航、定姿成为了众多国家、学者研究的热点。
自上世纪60年代开始，有学者便提出利用重力梯度仪辅助惯性导航。重力梯度辅助惯性导航，是指利用重力梯度值来修正惯性导航***(Inertial Navigation System，INS)随时间积累的定位误差，具有自主性强、隐蔽性好、不受地域和时域限制等优点。但由于惯性器件的固有缺陷(例如陀螺漂移和加速度计零偏置)使得重力场辅助导航***并不能从根本上解决INS误差积累的问题。而且，这些重力梯度导航研究中，主要应用场景局限于水下潜艇和航空飞行器，重力梯度导航没有成为一种独立的导航方法。
近年来，人们开始着眼于在轨卫星的重力梯度导航，为卫星的自主运行提供一种新方案。传统的卫星导航定轨方法主要是基于地面站遥测的定轨方法。然而，使用地面站对卫星进行定轨存在诸多限制，如建设成本、地面站故障、数据处理能力、数据传送能力、地域限制、政治因素等。对于目前拥有大量发射任务的我国来说，地面站的数量以及处理能力并不能完全满足当前需要。因此有必要提出一种新型的定轨方法来作为传统定轨方法的补充，重力梯度导航成为了避免以上弊端的有效方案之一，其具有建设成本低、完全不依赖外界提供的信息，星上自主性强，不易受外界因素干扰的优势。
当前，美国、俄罗斯以及中国等航天大国都提出了对于深层空间的探测计划。当探测器远离地球后，继续采用传统地面站遥测定轨会出现信号强度微弱，延迟严重，受到各种干扰等问题，将会给深空探测器实时定轨带来巨大的困难。而如果事先已建立了要探测行星的重力场模型，那么重力场梯度导航将为卫星提供自主的、实时的、可靠的导航、定轨结果。
综上，由于当前高精度地球重力场模型的建立，以及先进的重力测量仪器的出现为重力梯度导航提供了支持。随着研究的深入，使用重力梯度为航天器进行定轨成为了当前研究的热点。并且由于其相较于传统定轨方法所具有的优点与意义让其研究具有了必要性。本发明便提出了一种基于重力场梯度不变量的轨道要素估计方法，使用重力场梯度矩阵中的不变量信息估计出航天器轨道六要素中的五个要素。
发明内容
(一)发明目的：
本发明旨在提供一种新型的自主定轨、导航方法，该方法创新的使用了卫星运行过程中所测量的重力梯度矩阵不变量信息，通过迭代求解出卫星到地心的几何距离与地心纬度，并以此作为观测量给出了以轨道根数中的半长轴、偏心率、轨道倾角、近地点俯角和真近点角为状态参数的测量方程。为了提高轨道根数的估计精度，创新的引入了轨道根数的摄动动力学方程，采用批处理最小二乘估计出初始历元的五个轨道根数。该方法具有自主程度高、抗干扰能力强、建设成本低等优点，而且在深空探测领域具有传统方法所不具有的优势。本发明提出了一个完整的轨道根数估计方法流程，并为了减少过程中的计算难度，给出了较为简单的计算方法。
(二)技术方案
本发明一种基于重力场梯度不变量的轨道要素估计方法，其步骤如下。
步骤一：准备工作
地球引力势函数
在初步研究航天器运动规律时，把地球看作理想的球体，它对航天器产生的引力指向地心，其大小F与航天器的质量m成正比，而与地心至航天器的距离r的平方成反比：
上式中，F为航天器所受地球引力，G为地球引力常量，m为航天器质量，r为地心至航天器的几何距离。F除以m就是引力加速度g。把两个常数合并μ＝GM，则得到下式：
上式中，μ成为地球引力常数。
此情况下地球引力场成为中心引力场。把引力加速度表示成引力势的梯度，则中心引力场的引力势函数是：
实际的地球并不是球对称的，具有变平度、梨形和赤道变形等，因而引力加速度是距离r、纬度(地心纬度)和经度λ的函数，而r，λ是在地球固连坐标系Se中位置的球坐标。为了严格地而且方便地描述地球引力加速度的分布，把它表示成地球引力势函数U的梯度：
上式中grad(*)为求梯度算子。
以下给出，在只考虑地球非球形J2项影响下的引力场表达式：
上式中R＝6378.1363km，为地球赤道半径。
开普勒轨道根数
在惯性空间观察航天器轨道，需要定义轨道要素。轨道要素，又称轨道根数，是决定开普勒轨道的运动特征的一组常数。在考虑轨道摄动后，轨道根数不再是常数，可以作为状态变量。
轨道根数一共有6个，表示及含义如下：
a：椭圆轨道半场轴；
e：椭圆轨道偏心率；
Ω：升交点赤经；航天器轨道由南向北穿越赤道的点称为升交点B。赤道平面内，由春分点向东转到升交点B的角度称为升交点赤经，有效范围是0至360度；
i：轨道倾角；轨道平面与赤道平面的夹角。有效范围是0至180度；
tp：近地点时刻；航天器通过近地点的时刻，也可以用纪元时刻t0的真近点角θ(t0)，或偏近点角E(t0)，或平近点角M(t0)；
ω：近地点幅角；轨道平面内，有升交点转到近地点的角度。有效范围是0至360度。
他们的作用分别是：轨道半长轴a和偏心率e确定轨道的大小和形状；升交点赤经Ω和轨道倾角i确定轨道平面在惯性空间的取向；近地点幅角ω确定拱线在轨道平面内的位置；纪元时刻的平近点角M(t0)确定时间的起点。
轨道摄动方程
在摄动力作用下航天器的轨道不是开普勒轨道，其轨道要素是随时间变化的：a(t)，e(t)，Ω(t)，i(t)，…。
以下给出轨道摄动方程：
上式中，与前文已出现的符号具有相同含义；fr表示摄动加速度在地心轨道坐标系的径向分量；fu表示摄动加速度在地心轨道坐标系的横向分量；fh表示摄动加速度在地心轨道坐标系的副法向分量；p为半通径。
步骤二：理想重力梯度张量在东北天(East-North-Up，ENU)坐标系下的分解和特征值
理想球体重力场模型如公式(3)所示，式中r有如下表示：
上式中，x，y，z分别是卫星位置向量r在地球固连坐标系上的三个分量投影。
重力梯度张量可以看成是引力加速度的梯度，引力加速度又是引力场的梯度。因而直接在ENU坐标系下分解：
上式中，下标E，N，U表示Γ在ENU坐标系下分解的不同元素。
通过式(8)可以方便的得到重力梯度矩阵的三个特征值：
从式(9)可以知，三个特征值两两线性相关，只有一个特征量。
步骤三：求解测量历元J2模型重力场张量在ENU坐标系下的分解和特征值
J2模型重力场张量由式(5)给出。采用求梯度算子求解J2模型重力场张量在ENU坐标系下的分解：
其中，各个分量形式如下：
由式(10)求解矩阵的特征值λ1，λ2，λ3，特征值与的关系如下：
可以验证λ1，λ2，λ3有如下关系：
λ1+λ2+λ3＝0 (13)
因此，三个特征值中只有两个独立量。
步骤四：利用J2重力梯度矩阵特征值求解各测量历元r和
1)计算r初值
使用J2模型求解的重力梯度张量矩阵特征值(12)，和理想球体重力场模型特征值表达式(9)计算r初值。
式(12)中的特征值λ1和λ2任选其一皆可，带入式(9)，反解r则有如下初值表达式：
或
2)计算初值
由式(12)做如下计算：
从式(16)可以得到的表达式：
从式(17)可以得到的另一表达式：
以上的两种表达式，任选其一便可。对于式(18)与式(19)中余弦值的正负号在此不做讨论。
3)利用更新r
利用式(12)中的λ2表达式更新r，将r0和带入反解r3项，表达式如下：
4)利用r1更新
方法一：
用式(12)中的三个特征值做如下计算：
将r1和带入上式(23)，反解上式中的项，如下：
方法二：
使用式(21)和式(22)做如下计算：
将r1和带入上式(25)，反解上式中的项，如下：
步骤五：利用r和计算初始轨道要素
1)计算半长轴和偏心率
在所有计算历元的结果r中搜索最大值rmax和最小值rmin，则有以下关系：
ra＝rmax，rp＝rmin (27)
上式中，下标a表示远地点，下标p表示近地点。
则半长轴：
偏心率：
2)计算轨道倾角
在此假设，通过步骤四所求得的值的正负号已通过其他方法得知，在所有计算历元的结果中搜索最大值，
轨道倾角：
3)计算近地点幅角和初始真近点角
在2)中搜索到的rmin其对应的纬度记为则近地点幅角ω由下式计算：
初始历元真近点角θ0，由下式计算：
上式中，为第一个历元所计算出的纬度。
步骤六：轨道要素平滑
将初始历元五个轨道要素作为待估参数，如下：
x＝[a0 ex0 ey0 i0 u0] (33)
上式中，下标“0”表示初始历元参数。ex，ey，u表达式如下：
待估参数x的初值，由步骤五给出。
轨道动力学模型由式(6)描述，半通径由于J2项引起的摄动力加速度在地心轨道坐标系下的分解如下：
结合式(6)与式(35)，以及由步骤五提供的初值，便构成了已知初值的微分方程组求解问题，则通过数值积分方法可求得任意k历元时刻的五个轨道要素，xk＝[ak exk eyk ikuk]，则k历元时刻的观测方程如下：
上式中，为步骤四中所求得的各历元卫星到地心的几何距离与纬度，表示计算误差，下标“k”表示第k个历元。
简单地假设是一个高斯噪声，则k历元测量协方差矩阵如下：
上式中，σλ为重力梯度不变量误差，其误差大小要根据测量仪器的测量噪声与J2模型误差大小选取。
将各个历元的测量值z放置在一起写成如下形式：
上式中，N为测量历元的个数。各个历元测量协方差矩阵写成以下形式：
采用准西格玛最小二乘批处理滤波估计状态参数，其估计式如下：
上式中，下标j表示迭代次数，表示利用计算出的测量量z的计算值与敏度矩阵。
式(40)中，求解H阵难度较大，将Z写成关于x＝[a0 ex0 ry0 i0 u0]的隐函数表达，如下：
Z＝f(a0 ex0 ey0 i0 u0) (41)
通过普通解析法求偏导计算敏度矩阵H，计算量非常巨大。因此可以通过以下方法求解。
在第j次迭代完毕，得到状态参数估计值则进行如下采样：
上式中，L为x的维度，P0为先验协方差矩阵，表示P0的第i列平方根。先验协方差矩阵由下式给出：
将式(42)通过式(6)积分，和式(36)传递，计算测量值如下：
χi，j→γi，j，i＝0，…，L (44)
则密度矩阵H可以近似通过下式得出：
上式中，表示P0对角线元素的平方根。
通过上述步骤，提出了一种基于重力场梯度不变量的轨道要素估计方法。该方法采用重力梯度仪所测量的重力梯度信息，并利用理想球体重力场模型与带J2项的重力梯度模型求解重力梯度张量的特征值，给出了利用该特征值求解五个轨道要素的方法，为了提高求解精度引入轨道摄动动力学方程来平滑五个轨道要素。针对计算量大的部分，给出了相应的降低计算难度的方法。综上，该方法完全不依赖于外界信息的输入，具有很强的自主性，建设成本低在自主导航以及深孔探测领域具有很强的应用前景。若引入的更高精度的重力场将会进一步提高定轨精度，具有很强的潜力价值。
(三)优点
本发明提供的一种基于重力场梯度不变量的轨道要素估计方法方法的优点在于：
①本发明中，创新的使用了重力梯度测量矩阵中不变量信息，结合重力场模型给出了利用不变量信息计算轨道要素的方法。确定轨道要素过程中不需要额外引入其他外界信息，具有很强的自主性、抗干扰能力；
②本发明提出的方法，具有实现成本低的优势。不同于地面站定轨需要大量的资金投入进行***建设，本方法只需在卫星上搭载一台重力梯度测量仪便可实现自主定轨，实现成本低；
③本发明提出的方法，在深空探测领域具有巨大的应用前景。在掌握待勘测行星精确的重力场模型后，便可实现实时自主导航、定轨。
附图说明
图1是本发明所述方法流程图。
图2是本发明中利用J2重力梯度矩阵特征值求解各测量历元r和的流程图
具体实施方式
下面将结合附图1、图2和技术方案对本发明的具体实施过程做进一步的详细说明。
步骤一：理想重力梯度张量在东北天(ENU)坐标系下的分解和特征值
理想球体重力场模型如公式(3)所示，式中r有如下表示：
上式中，x，y，z分别是卫星位置向量r在地球固连坐标系上的三个分量投影。
重力梯度张量可以看成是引力加速度的梯度，引力加速度又是引力场的梯度。因而直接在ENU坐标系下分解：
上式中，下标E，N，U表示Γ在ENU坐标系下分解的不同元素。
通过式(47)可以方便的得到重力梯度矩阵的三个特征值：
从式(48)可以知，三个特征值两两线性相关，只有一个特征量。
该步骤对应附图1中的第二个方框。
步骤二：求解测量历元J2模型重力场张量在ENU坐标系下的分解和特征值
J2模型重力场张量由式(5)给出。通过求梯度计算求解J2模型重力场张量在ENU坐标系下的分解，得下式：
其中，各个分量形式如下：
由式(49)求解矩阵特征值λ1，λ2，λ3，特征值与的关系如下：
可以验证λ1，λ2，λ3有如下关系：
λ1+λ2+λ3＝0 (52)
因此，三个特征值中只有两个独立量。
该步骤对应附图1中第二个方框。
步骤三：利用J2重力梯度矩阵特征值求解各测量历元r和
1)计算r初值
使用J2模型求解的重力梯度张量矩阵特征值(51)，和理想球体重力场模型特征值表达式(48)计算r初值。
式(51)中的特征值λ1和λ2任选其一皆可，带入式(48)，反解r则有如下初值表达式：
或
该步骤对应附图2中第二个方框。
2)计算初值
由式(51)做如下计算：
从式(55)可以得到的表达式：
或，从式(56)可以得到的另一表达式：
以上的两种表达式，任选其一便可。对于式(57)与式(58)中余弦值的正负号在此不做讨论。
该步骤对应附图2中第三个方框。
3)利用更新r
利用式(51)中的λ2表达式更新r，将λ0和带入反解r3项，表达式如下：
该步骤对应附图2中第四个方框。
4)利用r1更新
方法一：
用式(51)中的三个特征值做如下计算：
将r1和带入上式(62)，反解上式中的项，如下：
方法二：
使用式(60)和式(61)做如下计算：
将r1和带入上式(64)，反解上式中的项，如下：
该步骤对应附图2中第五个方框。
5)判断几何距离r和纬度是否达到要求精度
如果达到符合精度则跳出循环，执行步骤四，否则执行3)。
该步骤对应附图2中的第六个方框。
步骤三对应附图1中的第四个方框。
步骤四：判断测量历元是否计算完毕
若所有历元计算完毕则执行步骤五，否则执行步骤二，计算下一历元重力梯度张量矩阵特征值。
该步骤对应附图1中的第五个方框。
步骤五：利用r和计算初始轨道要素
1)计算半长轴和偏心率
在所有计算历元的结果r中搜索最大值rmax和最小值rmin，则有以下关系：
ra＝rmax，rp＝rmin (66)
上式中，下标a表示远地点，下标p表示近地点。
则半长轴：
偏心率：
2)计算轨道倾角
在此假设，通过步骤四所求得的值的正负号已通过其他方法得知，在所有计算历元的结果中搜索最大值，
轨道倾角：
3)计算近地点幅角和初始真近点角
在2)中搜索到的rmin其对应的纬度记为则近地点幅角ω由下式计算：
初始历元真近点角θ0，由下式计算：
上式中，为第一个历元所计算出的纬度。
步骤六：轨道要素平滑
将初始历元五个轨道要素作为待估参数，如下：
x＝[a0 ex0 ey0 i0u0] (72)
上式中，下标“0”表示初始历元参数。ex，ey，u表达式如下：
待估参数x的初值，由步骤五给出。
轨道动力学模型由式(6)描述，半通径由于J2项引起的摄动力加速度在地心轨道坐标系下的分解如下：
结合式(6)与式(74)，以及由步骤五提供的初值，便构成了已知初值的微分方程组求解问题，则通过数值积分方法可求得任意k历元时刻的五个轨道要素，xk＝[ak exk eyk ikuk]，则k历元时刻的观测方程如下：
上式中，为步骤三中所求得的各历元卫星到地心的几何距离与纬度，表示计算误差，下标“k”表示第k个历元。
简单地假设是一个高斯噪声，则k历元测量协方差矩阵如下：
上式中，σλ为重力梯度不变量误差，其误差大小要根据测量仪器的测量噪声与J2模型误差大小选取。
将各个历元的测量值z放置在一起写成如下形式：
上式中，N为测量历元的个数。各个历元测量协方差矩阵写成以下形式：
采用准西格玛最小二乘批处理滤波估计状态参数，其估计式如下：
上式中，下标j表示迭代次数，表示利用计算出的测量量z的计算值与敏度矩阵。
式(79)中，求解H阵难度较大，将Z写成关于x＝[a0 ex0 ey0 i0 u0]的隐函数表达，如下：
Z＝f(a0 ex0 ey0 i0 u0) (80)
通过普通解析法求偏导计算敏度矩阵H，计算量非常巨大。因此可以通过以下方法求解。
在第j次迭代完毕，得到状态参数估计值则进行如下采样：
上式中，L为x的维度，P0为先验协方差矩阵，表示P0的第i列平方根。先验协方差矩阵由下式给出：
将式(81)通过式(6)积分，和式(75)传递，计算测量值如下：
χi，j→γi，j，i＝0，…，L (83)
则密度矩阵H可以近似通过下式得出：
上式中，表示P0对角线元素的平方根。
当最小二乘迭代达到精度要求时结束迭代输出最终平滑的轨道要素结果。
通过上述步骤，给出了一种基于重力场梯度不变量的轨道要素估计方法的完整流程。该流程采用重力梯度矩阵不变量信息，引入重力场模型与轨道摄动方程，得到五个轨道要素，并且简化了求解计算难度，降低了对计算机的性能要求。综上，该方法流程不依赖于外界信息的输入，在满足当前以及以后卫星等航天器自主运行的要求，建设成本远低于地面站建设，并且在未来的深空探测、自主导航领域具有光明的应用前景。

Claims (1)

1.一种基于重力场梯度不变量的轨道要素估计方法，其特征在于：其步骤如下：
步骤一：准备工作
地球引力势函数
在初步研究航天器运动规律时，把地球看作理想的球体，它对航天器产生的引力指向地心，其大小F与航天器的质量m成正比，而与地心至航天器的距离r的平方成反比：
<mrow> <mi>F</mi> <mo>=</mo> <mfrac> <mrow> <mi>G</mi> <mi>M</mi> <mi>m</mi> </mrow> <msup> <mi>r</mi> <mn>2</mn> </msup> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
上式中，F为航天器所受地球引力，G为地球引力常量，m为航天器质量，r为地心至航天器的几何距离。F除以m就是引力加速度g。把两个常数合并μ＝GM，则得到下式：
<mrow> <mi>g</mi> <mo>=</mo> <mfrac> <mi>&amp;mu;</mi> <msup> <mi>r</mi> <mn>2</mn> </msup> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
上式中，μ成为地球引力常数。
此情况下地球引力场成为中心引力场。把引力加速度表示成引力势的梯度，则中心引力场的引力势函数是：
<mrow> <mi>U</mi> <mo>=</mo> <mfrac> <mi>&amp;mu;</mi> <mi>r</mi> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>
实际的地球并不是球对称的，具有变平度、梨形和赤道变形等，因而引力加速度是距离r、纬度(地心纬度)和经度λ的函数，而r，λ是在地球固连坐标系Se中位置的球坐标。为了严格地而且方便地描述地球引力加速度的分布，把它表示成地球引力势函数U的梯度：
上式中grad(*)为求梯度算子。
以下给出，在只考虑地球非球形J2项影响下的引力场表达式：
上式中R＝6378.1363km，为地球赤道半径。
开普勒轨道根数
在惯性空间观察航天器轨道，需要定义轨道要素。轨道要素，又称轨道根数，是决定开普勒轨道的运动特征的一组常数。在考虑轨道摄动后，轨道根数不再是常数，可以作为状态变量。
轨道根数一共有6个，表示及含义如下：
a：椭圆轨道半场轴；
e：椭圆轨道偏心率；
Ω：升交点赤经；航天器轨道由南向北穿越赤道的点称为升交点B。赤道平面内，由春分点向东转到升交点B的角度称为升交点赤经，有效范围是0至360度；
i：轨道倾角；轨道平面与赤道平面的夹角。有效范围是0至180度；
tp：近地点时刻；航天器通过近地点的时刻，也可以用纪元时刻t0的真近点角θ(t0)，或偏近点角E(t0)，或平近点角M(t0)；
ω：近地点幅角；轨道平面内，有升交点转到近地点的角度。有效范围是0至360度。
他们的作用分别是：轨道半长轴a和偏心率e确定轨道的大小和形状；升交点赤经Ω和轨道倾角i确定轨道平面在惯性空间的取向；近地点幅角ω确定拱线在轨道平面内的位置；纪元时刻的平近点角M(t0)确定时间的起点。
轨道摄动方程
在摄动力作用下航天器的轨道不是开普勒轨道，其轨道要素是随时间变化的：a(t)，e(t)，Ω(t)，i(t)，…。
以下给出轨道摄动方程：
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mfrac> <mrow> <mi>d</mi> <mi>a</mi> </mrow> <mrow> <mi>d</mi> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mfrac> <mrow> <mn>2</mn> <msup> <mi>a</mi> <mn>2</mn> </msup> </mrow> <msqrt> <mrow> <mi>&amp;mu;</mi> <mi>p</mi> </mrow> </msqrt> </mfrac> <mrow> <mo>&amp;lsqb;</mo> <mrow> <mi>e</mi> <mi> </mi> <msub> <mi>sin&amp;theta;f</mi> <mi>r</mi> </msub> <mo>+</mo> <mrow> <mo>(</mo> <mrow> <mn>1</mn> <mo>+</mo> <mi>cos</mi> <mi>&amp;theta;</mi> </mrow> <mo>)</mo> </mrow> <msub> <mi>f</mi> <mi>u</mi> </msub> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mrow> <mi>d</mi> <mi>e</mi> </mrow> <mrow> <mi>d</mi> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <msqrt> <mfrac> <mi>p</mi> <mi>&amp;mu;</mi> </mfrac> </msqrt> <mrow> <mo>{</mo> <mrow> <msub> <mi>sin&amp;theta;f</mi> <mi>r</mi> </msub> <mo>+</mo> <mrow> <mo>&amp;lsqb;</mo> <mrow> <mrow> <mo>(</mo> <mrow> <mn>1</mn> <mo>+</mo> <mfrac> <mi>r</mi> <mi>p</mi> </mfrac> </mrow> <mo>)</mo> </mrow> <mi>cos</mi> <mi>&amp;theta;</mi> <mo>+</mo> <mfrac> <mrow> <mi>e</mi> <mi>r</mi> </mrow> <mi>p</mi> </mfrac> </mrow> <mo>&amp;rsqb;</mo> </mrow> <msub> <mi>f</mi> <mi>u</mi> </msub> </mrow> <mo>}</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mrow> <mi>d</mi> <mi>i</mi> </mrow> <mrow> <mi>d</mi> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mfrac> <mrow> <mi>r</mi> <mi> </mi> <mi>cos</mi> <mrow> <mo>(</mo> <mrow> <mi>&amp;omega;</mi> <mo>+</mo> <mi>&amp;theta;</mi> </mrow> <mo>)</mo> </mrow> </mrow> <msqrt> <mrow> <mi>&amp;mu;</mi> <mi>p</mi> </mrow> </msqrt> </mfrac> <msub> <mi>f</mi> <mi>h</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mrow> <mi>d</mi> <mi>&amp;omega;</mi> </mrow> <mrow> <mi>d</mi> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mfrac> <mn>1</mn> <mi>e</mi> </mfrac> <msqrt> <mfrac> <mi>p</mi> <mi>&amp;mu;</mi> </mfrac> </msqrt> <mrow> <mo>&amp;lsqb;</mo> <mrow> <mo>-</mo> <msub> <mi>cos&amp;theta;f</mi> <mi>r</mi> </msub> <mo>+</mo> <mrow> <mo>(</mo> <mrow> <mn>1</mn> <mo>+</mo> <mfrac> <mi>r</mi> <mi>p</mi> </mfrac> </mrow> <mo>)</mo> </mrow> <msub> <mi>sin&amp;theta;f</mi> <mi>u</mi> </msub> <mo>-</mo> <mfrac> <mrow> <mi>e</mi> <mi>r</mi> </mrow> <mi>p</mi> </mfrac> <mfrac> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <mrow> <mi>&amp;omega;</mi> <mo>+</mo> <mi>&amp;theta;</mi> </mrow> <mo>)</mo> </mrow> </mrow> <mrow> <mi>tan</mi> <mi> </mi> <mi>i</mi> </mrow> </mfrac> <msub> <mi>f</mi> <mi>h</mi> </msub> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mrow> <mi>d</mi> <mi>M</mi> </mrow> <mrow> <mi>d</mi> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <msqrt> <mfrac> <mi>&amp;mu;</mi> <msup> <mi>a</mi> <mn>3</mn> </msup> </mfrac> </msqrt> <mo>+</mo> <mfrac> <msqrt> <mrow> <mn>1</mn> <mo>-</mo> <msup> <mi>e</mi> <mn>2</mn> </msup> </mrow> </msqrt> <mi>e</mi> </mfrac> <msqrt> <mfrac> <mi>p</mi> <mi>&amp;mu;</mi> </mfrac> </msqrt> <mrow> <mo>&amp;lsqb;</mo> <mrow> <mrow> <mo>(</mo> <mrow> <mi>cos</mi> <mi>&amp;theta;</mi> <mo>-</mo> <mfrac> <mrow> <mn>2</mn> <mi>e</mi> <mi>r</mi> </mrow> <mi>p</mi> </mfrac> </mrow> <mo>)</mo> </mrow> <msub> <mi>f</mi> <mi>r</mi> </msub> <mo>-</mo> <mrow> <mo>(</mo> <mrow> <mn>1</mn> <mo>+</mo> <mfrac> <mi>r</mi> <mi>p</mi> </mfrac> </mrow> <mo>)</mo> </mrow> <msub> <mi>sin&amp;theta;f</mi> <mi>u</mi> </msub> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>
上式中，与前文已出现的符号具有相同含义；fr表示摄动加速度在地心轨道坐标系的径向分量；fu表示摄动加速度在地心轨道坐标系的横向分量；fh表示摄动加速度在地心轨道坐标系的副法向分量；p为半通径。
步骤二：理想重力梯度张量在东北天(East-North-Up，ENU)坐标系下的分解和特征值
理想球体重力场模型如公式(3)所示，式中r有如下表示：
<mrow> <mi>r</mi> <mo>=</mo> <msqrt> <mrow> <mo>(</mo> <mrow> <msup> <mi>x</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>y</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>z</mi> <mn>2</mn> </msup> </mrow> <mo>)</mo> </mrow> </msqrt> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>
上式中，x，y，z分别是卫星位置向量r在地球固连坐标系上的三个分量投影。重力梯度张量可以看成是引力加速度的梯度，引力加速度又是引力场的梯度。因而直接在ENU坐标系下分解：
<mrow> <mover> <mi>T</mi> <mo>&amp;OverBar;</mo> </mover> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mover> <mi>T</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>E</mi> <mi>E</mi> </mrow> </msub> </mtd> <mtd> <msub> <mover> <mi>T</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>E</mi> <mi>N</mi> </mrow> </msub> </mtd> <mtd> <msub> <mover> <mi>T</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>E</mi> <mi>U</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mover> <mi>T</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>N</mi> <mi>E</mi> </mrow> </msub> </mtd> <mtd> <msub> <mover> <mi>T</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>N</mi> <mi>N</mi> </mrow> </msub> </mtd> <mtd> <msub> <mover> <mi>T</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>N</mi> <mi>U</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mover> <mi>T</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>U</mi> <mi>E</mi> </mrow> </msub> </mtd> <mtd> <msub> <mover> <mi>T</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>U</mi> <mi>N</mi> </mrow> </msub> </mtd> <mtd> <msub> <mover> <mi>T</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>U</mi> <mi>U</mi> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mfrac> <mi>&amp;mu;</mi> <msup> <mi>r</mi> <mn>3</mn> </msup> </mfrac> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>2</mn> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow>
上式中，下标E，N，U表示Γ在ENU坐标系下分解的不同元素。
通过式(8)可以方便的得到重力梯度矩阵的三个特征值：
<mrow> <msub> <mi>&amp;lambda;</mi> <mn>1</mn> </msub> <mo>=</mo> <mo>-</mo> <mfrac> <mi>&amp;mu;</mi> <msup> <mi>r</mi> <mn>3</mn> </msup> </mfrac> </mrow>
<mrow> <msub> <mi>&amp;lambda;</mi> <mn>2</mn> </msub> <mo>=</mo> <mo>-</mo> <mfrac> <mi>&amp;mu;</mi> <msup> <mi>r</mi> <mn>3</mn> </msup> </mfrac> </mrow>
<mrow> <msub> <mi>&amp;lambda;</mi> <mn>3</mn> </msub> <mo>=</mo> <mfrac> <mrow> <mn>2</mn> <mi>&amp;mu;</mi> </mrow> <msup> <mi>r</mi> <mn>3</mn> </msup> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow> 2
从式(9)可以知，三个特征值两两线性相关，只有一个特征量。
步骤三：求解测量历元J2模型重力场张量在ENU坐标系下的分解和特征值
J2模型重力场张量由式(5)给出。采用求梯度算子求解J2模型重力场张量在ENU坐标系下的分解：
<mrow> <mover> <mi>T</mi> <mo>&amp;OverBar;</mo> </mover> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mover> <mi>T</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>E</mi> <mi>E</mi> </mrow> </msub> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <msub> <mover> <mi>T</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>N</mi> <mi>N</mi> </mrow> </msub> </mtd> <mtd> <msub> <mover> <mi>T</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>N</mi> <mi>U</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <msub> <mover> <mi>T</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>U</mi> <mi>N</mi> </mrow> </msub> </mtd> <mtd> <msub> <mover> <mi>T</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>U</mi> <mi>U</mi> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow>
其中，各个分量形式如下：
由式(10)求解矩阵的特征值λ1，λ2，λ3，特征值与r，的关系如下：
可以验证λ1，λ2，λ3有如下关系：
λ1+λ2+λ3＝0 (13)
因此，三个特征值中只有两个独立量。
步骤四：利用J2重力梯度矩阵特征值求解各测量历元r和
1)计算r初值
使用J2模型求解的重力梯度张量矩阵特征值(12)，和理想球体重力场模型特征值表达式(9)计算r初值。
式(12)中的特征值λ1和λ2任选其一皆可，带入式(9)，反解r则有如下初值表达式：
<mrow> <msub> <mi>r</mi> <mn>0</mn> </msub> <mo>=</mo> <msup> <mrow> <mo>(</mo> <mrow> <mo>-</mo> <mfrac> <mi>&amp;mu;</mi> <msub> <mi>&amp;lambda;</mi> <mn>2</mn> </msub> </mfrac> </mrow> <mo>)</mo> </mrow> <mrow> <mn>1</mn> <mo>/</mo> <mn>3</mn> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo> </mrow> </mrow> 3
或
<mrow> <msub> <mi>r</mi> <mn>0</mn> </msub> <mo>=</mo> <msup> <mrow> <mo>(</mo> <mrow> <mo>-</mo> <mfrac> <mi>&amp;mu;</mi> <msub> <mi>&amp;lambda;</mi> <mn>1</mn> </msub> </mfrac> </mrow> <mo>)</mo> </mrow> <mrow> <mn>1</mn> <mo>/</mo> <mn>3</mn> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>15</mn> <mo>)</mo> </mrow> </mrow>
2)计算初值
由式(12)做如下计算：
从式(16)可以得到的表达式：
从式(17)可以得到的另一表达式：
以上的两种表达式，任选其一便可。对于式(18)与式(19)中余弦值的正负号在此不做讨论。
3)利用更新r
利用式(12)中的λ2表达式更新r，将r0和带入反解r3项，表达式如下：
4)利用r1更新
方法一：
用式(12)中的三个特征值做如下计算：
将r1和带入上式(23)，反解上式中的项，如下：
方法二：
使用式(21)和式(22)做如下计算：
将r1和带入上式(25)，反解上式中的项，如下：
步骤五：利用r和计算初始轨道要素
1)计算半长轴和偏心率
在所有计算历元的结果r中搜索最大值rmax和最小值rmin，则有以下关系：
ra＝rmax，rp＝rmin (27)
上式中，下标a表示远地点，下标p表示近地点。
则半长轴：
<mrow> <mi>a</mi> <mo>=</mo> <mfrac> <mrow> <msub> <mi>r</mi> <mi>a</mi> </msub> <mo>+</mo> <msub> <mi>r</mi> <mi>p</mi> </msub> </mrow> <mn>2</mn> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>28</mn> <mo>)</mo> </mrow> </mrow>
偏心率：
<mrow> <mi>e</mi> <mo>=</mo> <mfrac> <mrow> <msub> <mi>r</mi> <mi>a</mi> </msub> <mo>-</mo> <msub> <mi>r</mi> <mi>p</mi> </msub> </mrow> <mrow> <msub> <mi>r</mi> <mi>a</mi> </msub> <mo>+</mo> <msub> <mi>r</mi> <mi>p</mi> </msub> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>29</mn> <mo>)</mo> </mrow> </mrow>
2)计算轨道倾角
在此假设，通过步骤四所求得的值的正负号已通过其他方法得知，在所有计算历元的结果中搜索最大值，
轨道倾角：
3)计算近地点幅角和初始真近点角
在2)中搜索到的rmin其对应的纬度记为则近地点幅角ω由下式计算：
初始历元真近点角θ0，由下式计算：
上式中，为第一个历元所计算出的纬度。
步骤六：轨道要素平滑
将初始历元五个轨道要素作为待估参数，如下：
<mrow> <mi>x</mi> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>a</mi> <mn>0</mn> </msub> </mtd> <mtd> <msub> <mi>e</mi> <mrow> <mi>x</mi> <mn>0</mn> </mrow> </msub> </mtd> <mtd> <msub> <mi>e</mi> <mrow> <mi>y</mi> <mn>0</mn> </mrow> </msub> </mtd> <mtd> <msub> <mi>i</mi> <mn>0</mn> </msub> </mtd> <mtd> <msub> <mi>u</mi> <mn>0</mn> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>33</mn> <mo>)</mo> </mrow> </mrow>
上式中，下标“0”表示初始历元参数。ex，ey，u表达式如下：
<mrow> <mtable> <mtr> <mtd> <mrow> <msub> <mi>e</mi> <mi>x</mi> </msub> <mo>=</mo> <mi>e</mi> <mi> </mi> <mi>cos</mi> <mi>&amp;omega;</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>e</mi> <mi>y</mi> </msub> <mo>=</mo> <mi>e</mi> <mi> </mi> <mi>sin</mi> <mi>&amp;omega;</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>u</mi> <mo>=</mo> <mi>&amp;omega;</mi> <mo>+</mo> <mi>&amp;theta;</mi> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>34</mn> <mo>)</mo> </mrow> </mrow>
待估参数x的初值，由步骤五给出。
轨道动力学模型由式(6)描述，半通径由于J2项引起的摄动力加速度在地心轨道坐标系下的分解如下：
<mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>f</mi> <mi>r</mi> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>f</mi> <mi>u</mi> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>f</mi> <mi>h</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mfrac> <mrow> <mn>3</mn> <msub> <mi>&amp;mu;J</mi> <mn>2</mn> </msub> <msup> <mi>R</mi> <mn>2</mn> </msup> </mrow> <mrow> <mn>2</mn> <msup> <mi>r</mi> <mn>4</mn> </msup> </mrow> </mfrac> <mrow> <mo>(</mo> <mrow> <mn>3</mn> <msup> <mi>sin</mi> <mn>2</mn> </msup> <mi>i</mi> <mi> </mi> <msup> <mi>sin</mi> <mn>2</mn> </msup> <mi>u</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>-</mo> <mfrac> <mrow> <mn>3</mn> <msub> <mi>&amp;mu;J</mi> <mn>2</mn> </msub> <msup> <mi>R</mi> <mn>2</mn> </msup> </mrow> <mrow> <mn>2</mn> <msup> <mi>r</mi> <mn>4</mn> </msup> </mrow> </mfrac> <msup> <mi>sin</mi> <mn>2</mn> </msup> <mi>i</mi> <mi> </mi> <mi>sin</mi> <mrow> <mo>(</mo> <mrow> <mn>2</mn> <mi>u</mi> </mrow> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>-</mo> <mfrac> <mrow> <mn>3</mn> <msub> <mi>&amp;mu;J</mi> <mn>2</mn> </msub> <msup> <mi>R</mi> <mn>2</mn> </msup> </mrow> <mrow> <mn>2</mn> <msup> <mi>r</mi> <mn>4</mn> </msup> </mrow> </mfrac> <mi>sin</mi> <mn>2</mn> <mi>i</mi> <mi> </mi> <mi>sin</mi> <mi>u</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>35</mn> <mo>)</mo> </mrow> </mrow>
结合式(6)与式(35)，以及由步骤五提供的初值，便构成了已知初值的微分方程组求解问题，则通过数值积分方法可求得任意k历元时刻的五个轨道要素，则k历元时刻的观测方程如下：
上式中，为步骤四中所求得的各历元卫星到地心的几何距离与纬度，Δr，表示计算误差，下标“k”表示第k个历元。
简单地假设Δr，是一个高斯噪声，则k历元测量协方差矩阵如下：
上式中，σλ为重力梯度不变量误差，其误差大小要根据测量仪器的测量噪声与J2模型误差大小选取。
将各个历元的测量值z放置在一起写成如下形式：
<mrow> <mi>Z</mi> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>z</mi> <mn>1</mn> </msub> </mtd> </mtr> <mtr> <mtd> <mtable> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> </mtable> </mtd> </mtr> <mtr> <mtd> <msub> <mi>z</mi> <mi>N</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>38</mn> <mo>)</mo> </mrow> </mrow>
上式中，N为测量历元的个数。各个历元测量协方差矩阵写成以下形式：
采用准西格玛最小二乘批处理滤波估计状态参数，其估计式如下：
<mrow> <msub> <mover> <mi>x</mi> <mo>^</mo> </mover> <mrow> <mi>j</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <msub> <mover> <mi>x</mi> <mo>^</mo> </mover> <mi>j</mi> </msub> <mo>+</mo> <msup> <mrow> <mo>(</mo> <mrow> <msubsup> <mi>H</mi> <mi>j</mi> <mi>T</mi> </msubsup> <msup> <mi>R</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <msub> <mi>H</mi> <mi>j</mi> </msub> </mrow> <mo>)</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <msubsup> <mi>H</mi> <mi>j</mi> <mi>T</mi> </msubsup> <msup> <mi>R</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>&amp;lsqb;</mo> <mrow> <mi>Z</mi> <mo>-</mo> <msub> <mover> <mi>Z</mi> <mo>&amp;OverBar;</mo> </mover> <mi>j</mi> </msub> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>40</mn> <mo>)</mo> </mrow> </mrow>
上式中，下标j表示迭代次数，Hj表示利用计算出的测量量z的计算值与敏度矩阵。
式(40)中，求解H阵难度较大，将Z写成关于的隐函数表达，如下：
<mrow> <mi>Z</mi> <mo>=</mo> <mi>f</mi> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <msub> <mi>a</mi> <mn>0</mn> </msub> </mtd> <mtd> <msub> <mi>e</mi> <mrow> <mi>x</mi> <mn>0</mn> </mrow> </msub> </mtd> <mtd> <msub> <mi>e</mi> <mrow> <mi>y</mi> <mn>0</mn> </mrow> </msub> </mtd> <mtd> <msub> <mi>i</mi> <mn>0</mn> </msub> </mtd> <mtd> <msub> <mi>u</mi> <mn>0</mn> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>41</mn> <mo>)</mo> </mrow> </mrow>
通过普通解析法求偏导计算敏度矩阵H，计算量非常巨大。因此可以通过以下方法求解。
在第j次迭代完毕，得到状态参数估计值则进行如下采样：
<mrow> <mtable> <mtr> <mtd> <mrow> <msub> <mi>&amp;chi;</mi> <mrow> <mn>0</mn> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>=</mo> <msub> <mover> <mi>x</mi> <mo>^</mo> </mover> <mi>j</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&amp;chi;</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>=</mo> <msub> <mover> <mi>x</mi> <mo>^</mo> </mover> <mi>j</mi> </msub> <mo>+</mo> <msub> <mrow> <mo>(</mo> <msqrt> <msub> <mi>P</mi> <mn>0</mn> </msub> </msqrt> <mo>)</mo> </mrow> <mi>i</mi> </msub> <mo>,</mo> <mi>i</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mi>L</mi> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>42</mn> <mo>)</mo> </mrow> </mrow>
上式中，L为x的维度，P0为先验协方差矩阵，表示P0的第i列平方根。先验协方差矩阵由下式给出：
<mrow> <msub> <mi>P</mi> <mn>0</mn> </msub> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msup> <mrow> <mo>(</mo> <msup> <mn>10</mn> <mrow> <mo>-</mo> <mn>6</mn> </mrow> </msup> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <msup> <mrow> <mo>(</mo> <msup> <mn>10</mn> <mrow> <mo>-</mo> <mn>6</mn> </mrow> </msup> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <msup> <mrow> <mo>(</mo> <msup> <mn>10</mn> <mrow> <mo>-</mo> <mn>6</mn> </mrow> </msup> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <msup> <mrow> <mo>(</mo> <msup> <mn>10</mn> <mrow> <mo>-</mo> <mn>6</mn> </mrow> </msup> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <msup> <mrow> <mo>(</mo> <msup> <mn>10</mn> <mrow> <mo>-</mo> <mn>6</mn> </mrow> </msup> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>43</mn> <mo>)</mo> </mrow> </mrow>
将式(42)通过式(6)积分，和式(36)传递，计算测量值如下：
χi，j→γi，j，i＝0，...，L (44)
则密度矩阵H可以近似通过下式得出：
<mrow> <msub> <mi>H</mi> <mi>j</mi> </msub> <mo>&amp;cong;</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mfrac> <mrow> <msub> <mi>&amp;gamma;</mi> <mrow> <mn>1</mn> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>&amp;gamma;</mi> <mrow> <mn>0</mn> <mo>,</mo> <mi>j</mi> </mrow> </msub> </mrow> <msub> <mrow> <mo>(</mo> <msqrt> <msub> <mi>P</mi> <mn>0</mn> </msub> </msqrt> <mo>)</mo> </mrow> <mn>11</mn> </msub> </mfrac> </mtd> <mtd> <mfrac> <mrow> <msub> <mi>&amp;gamma;</mi> <mrow> <mn>2</mn> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>&amp;gamma;</mi> <mrow> <mn>0</mn> <mo>,</mo> <mi>j</mi> </mrow> </msub> </mrow> <msub> <mrow> <mo>(</mo> <msqrt> <msub> <mi>P</mi> <mn>0</mn> </msub> </msqrt> <mo>)</mo> </mrow> <mn>22</mn> </msub> </mfrac> </mtd> <mtd> <mfrac> <mrow> <msub> <mi>&amp;gamma;</mi> <mrow> <mn>3</mn> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>&amp;gamma;</mi> <mrow> <mn>0</mn> <mo>,</mo> <mi>j</mi> </mrow> </msub> </mrow> <msub> <mrow> <mo>(</mo> <msqrt> <msub> <mi>P</mi> <mn>0</mn> </msub> </msqrt> <mo>)</mo> </mrow> <mn>33</mn> </msub> </mfrac> </mtd> <mtd> <mfrac> <mrow> <msub> <mi>&amp;gamma;</mi> <mrow> <mn>4</mn> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>&amp;gamma;</mi> <mrow> <mn>0</mn> <mo>,</mo> <mi>j</mi> </mrow> </msub> </mrow> <msub> <mrow> <mo>(</mo> <msqrt> <msub> <mi>P</mi> <mn>0</mn> </msub> </msqrt> <mo>)</mo> </mrow> <mn>44</mn> </msub> </mfrac> </mtd> <mtd> <mfrac> <mrow> <msub> <mi>&amp;gamma;</mi> <mrow> <mn>5</mn> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>&amp;gamma;</mi> <mrow> <mn>0</mn> <mo>,</mo> <mi>j</mi> </mrow> </msub> </mrow> <msub> <mrow> <mo>(</mo> <msqrt> <msub> <mi>P</mi> <mn>0</mn> </msub> </msqrt> <mo>)</mo> </mrow> <mn>55</mn> </msub> </mfrac> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>45</mn> <mo>)</mo> </mrow> </mrow>
上式中，表示P0对角线元素的平方根。
CN201710024671.5A 2017-01-13 2017-01-13 一种基于重力场梯度不变量的轨道要素估计方法 Active CN107065025B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710024671.5A CN107065025B (zh) 2017-01-13 2017-01-13 一种基于重力场梯度不变量的轨道要素估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710024671.5A CN107065025B (zh) 2017-01-13 2017-01-13 一种基于重力场梯度不变量的轨道要素估计方法

Publications (2)

Publication Number Publication Date
CN107065025A true CN107065025A (zh) 2017-08-18
CN107065025B CN107065025B (zh) 2019-04-23

Family

ID=59599280

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710024671.5A Active CN107065025B (zh) 2017-01-13 2017-01-13 一种基于重力场梯度不变量的轨道要素估计方法

Country Status (1)

Country Link
CN (1) CN107065025B (zh)

Cited By (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108415244A (zh) * 2017-12-28 2018-08-17 中国空间技术研究院 一种基于迭代调节的多自由度静电悬浮***几何对称性逼近方法
CN108548542A (zh) * 2018-07-13 2018-09-18 北京航空航天大学 一种基于大气阻力加速度测量的近地轨道确定方法
CN108919370A (zh) * 2018-07-25 2018-11-30 赛特雷德（重庆）科技有限公司 一种基于引力场测量的定位装置及方法
CN109738919A (zh) * 2019-02-28 2019-05-10 西安开阳微电子有限公司 一种用于gps接收机自主预测星历的方法
CN110053788A (zh) * 2019-03-15 2019-07-26 中国西安卫星测控中心 一种考虑复杂摄动的星座长期保持控制频次估计方法
CN110378012A (zh) * 2019-07-16 2019-10-25 上海交通大学 一种考虑高阶重力场的严格回归轨道设计方法
CN110442831A (zh) * 2019-07-31 2019-11-12 中国人民解放军国防科技大学 基于非线性偏差演化的空间非合作目标天基搜索方法
CN110705022A (zh) * 2019-08-30 2020-01-17 中国矿业大学 一种稀疏球面径向基函数局部重力场建模方法
CN110967041A (zh) * 2019-12-18 2020-04-07 自然资源部国土卫星遥感应用中心 基于张量不变理论的卫星引力梯度数据精度的验证方法
CN111505679A (zh) * 2020-04-20 2020-08-07 中国科学院国家空间科学中心 一种基于星载gnss的leo初轨确定方法
CN111721303A (zh) * 2020-07-01 2020-09-29 中国人民解放军陆军装甲兵学院 基于引力场的航天器惯性导航方法、***、介质及设备
CN112762924A (zh) * 2020-12-23 2021-05-07 北京机电工程研究所 基于重力梯度-地形异源数据匹配的导航定位方法
CN113433596A (zh) * 2021-06-25 2021-09-24 中国船舶重工集团公司第七0七研究所 一种基于空间域的重力梯度动态测量滤波方法
CN115327653A (zh) * 2022-08-15 2022-11-11 自然资源部国土卫星遥感应用中心 一种基于张量不变理论的卫星引力梯度粗差探测方法
CN115792982A (zh) * 2022-11-07 2023-03-14 北京航空航天大学合肥创新研究院(北京航空航天大学合肥研究生院) 北斗导航卫星广播星历参数拟合方法及存储介质
CN115808881A (zh) * 2023-01-21 2023-03-17 中国科学院数学与***科学研究院 一种无拖曳卫星在轨质量估计方法和自适应控制方法
CN116108319A (zh) * 2023-04-10 2023-05-12 中国人民解放军32035部队 一种恒定推力模式持续机动卫星的轨道预报方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2010132111A2 (en) * 2009-05-15 2010-11-18 Lockheed Martin Corporation Improved gravity sensing instrument
CN104061932A (zh) * 2014-06-10 2014-09-24 中国空间技术研究院 一种利用引力矢量和梯度张量进行导航定位的方法
CN105138005A (zh) * 2015-06-16 2015-12-09 北京航空航天大学 一种基于星间距离的相对轨道要素确定方法
CN106092099A (zh) * 2016-06-02 2016-11-09 哈尔滨工业大学 航天器相对位置增量定轨方法
CN106202640A (zh) * 2016-06-28 2016-12-07 西北工业大学 日‑地三体引力场中的晕轨道航天器偏置轨道设计方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2010132111A2 (en) * 2009-05-15 2010-11-18 Lockheed Martin Corporation Improved gravity sensing instrument
CN104061932A (zh) * 2014-06-10 2014-09-24 中国空间技术研究院 一种利用引力矢量和梯度张量进行导航定位的方法
CN105138005A (zh) * 2015-06-16 2015-12-09 北京航空航天大学 一种基于星间距离的相对轨道要素确定方法
CN106092099A (zh) * 2016-06-02 2016-11-09 哈尔滨工业大学 航天器相对位置增量定轨方法
CN106202640A (zh) * 2016-06-28 2016-12-07 西北工业大学 日‑地三体引力场中的晕轨道航天器偏置轨道设计方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
PEI CHEN 等: "Gravity Gradient Tensor Eigendecomposition for Spacecraft Positioning", 《JOURNAL OF GUIDANCE，CONTROL，AND DYNAMICS》 *
于锦海 等: "利用引力梯度不变量解算的GOCE引力场模型", 《中国科学：地球科学》 *
李建成 等: "由GOCE引力梯度张量不变量确定卫星重力模型的半解析法", 《武汉大学学报·信息科学版》 *

Cited By (27)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108415244B (zh) * 2017-12-28 2020-12-18 中国空间技术研究院 一种多自由度静电悬浮***几何对称性逼近方法
CN108415244A (zh) * 2017-12-28 2018-08-17 中国空间技术研究院 一种基于迭代调节的多自由度静电悬浮***几何对称性逼近方法
CN108548542A (zh) * 2018-07-13 2018-09-18 北京航空航天大学 一种基于大气阻力加速度测量的近地轨道确定方法
CN108548542B (zh) * 2018-07-13 2021-09-28 北京航空航天大学 一种基于大气阻力加速度测量的近地轨道确定方法
CN108919370A (zh) * 2018-07-25 2018-11-30 赛特雷德（重庆）科技有限公司 一种基于引力场测量的定位装置及方法
CN108919370B (zh) * 2018-07-25 2019-11-29 赛德雷特(珠海)航天科技有限公司 一种基于引力场测量的定位装置及方法
CN109738919A (zh) * 2019-02-28 2019-05-10 西安开阳微电子有限公司 一种用于gps接收机自主预测星历的方法
CN110053788A (zh) * 2019-03-15 2019-07-26 中国西安卫星测控中心 一种考虑复杂摄动的星座长期保持控制频次估计方法
CN110378012A (zh) * 2019-07-16 2019-10-25 上海交通大学 一种考虑高阶重力场的严格回归轨道设计方法
CN110378012B (zh) * 2019-07-16 2021-07-16 上海交通大学 一种考虑高阶重力场的严格回归轨道设计方法、***及介质
CN110442831A (zh) * 2019-07-31 2019-11-12 中国人民解放军国防科技大学 基于非线性偏差演化的空间非合作目标天基搜索方法
CN110705022A (zh) * 2019-08-30 2020-01-17 中国矿业大学 一种稀疏球面径向基函数局部重力场建模方法
CN110967041A (zh) * 2019-12-18 2020-04-07 自然资源部国土卫星遥感应用中心 基于张量不变理论的卫星引力梯度数据精度的验证方法
CN110967041B (zh) * 2019-12-18 2021-09-14 自然资源部国土卫星遥感应用中心 基于张量不变理论的卫星引力梯度数据精度的验证方法
CN111505679A (zh) * 2020-04-20 2020-08-07 中国科学院国家空间科学中心 一种基于星载gnss的leo初轨确定方法
CN111505679B (zh) * 2020-04-20 2022-05-03 中国科学院国家空间科学中心 一种基于星载gnss的leo初轨确定方法
CN111721303A (zh) * 2020-07-01 2020-09-29 中国人民解放军陆军装甲兵学院 基于引力场的航天器惯性导航方法、***、介质及设备
CN112762924A (zh) * 2020-12-23 2021-05-07 北京机电工程研究所 基于重力梯度-地形异源数据匹配的导航定位方法
CN112762924B (zh) * 2020-12-23 2023-07-14 北京机电工程研究所 基于重力梯度-地形异源数据匹配的导航定位方法
CN113433596B (zh) * 2021-06-25 2022-09-16 中国船舶重工集团公司第七0七研究所 一种基于空间域的重力梯度动态测量滤波方法
CN113433596A (zh) * 2021-06-25 2021-09-24 中国船舶重工集团公司第七0七研究所 一种基于空间域的重力梯度动态测量滤波方法
CN115327653A (zh) * 2022-08-15 2022-11-11 自然资源部国土卫星遥感应用中心 一种基于张量不变理论的卫星引力梯度粗差探测方法
CN115792982A (zh) * 2022-11-07 2023-03-14 北京航空航天大学合肥创新研究院(北京航空航天大学合肥研究生院) 北斗导航卫星广播星历参数拟合方法及存储介质
CN115792982B (zh) * 2022-11-07 2023-10-20 北京航空航天大学合肥创新研究院(北京航空航天大学合肥研究生院) 北斗导航卫星广播星历参数拟合方法及存储介质
CN115808881A (zh) * 2023-01-21 2023-03-17 中国科学院数学与***科学研究院 一种无拖曳卫星在轨质量估计方法和自适应控制方法
CN116108319A (zh) * 2023-04-10 2023-05-12 中国人民解放军32035部队 一种恒定推力模式持续机动卫星的轨道预报方法
CN116108319B (zh) * 2023-04-10 2023-08-11 中国人民解放军32035部队 一种恒定推力模式持续机动卫星的轨道预报方法

Also Published As

Publication number Publication date
CN107065025B (zh) 2019-04-23

Similar Documents

Publication Publication Date Title
CN107065025B (zh) 一种基于重力场梯度不变量的轨道要素估计方法
CN109556632B (zh) 一种基于卡尔曼滤波的ins/gnss/偏振/地磁组合导航对准方法
Chatfield Fundamentals of high accuracy inertial navigation
CN101344391B (zh) 基于全功能太阳罗盘的月球车位姿自主确定方法
CN109556631B (zh) 一种基于最小二乘的ins/gnss/偏振/地磁组合导航***对准方法
CN102353378B (zh) 一种矢量形式信息分配系数的组合导航***自适应联邦滤波方法
CN106871928A (zh) 基于李群滤波的捷联惯性导航初始对准方法
CN112161632B (zh) 一种基于相对位置矢量测量的卫星编队初始定位方法
CN105953795A (zh) 一种用于航天器表面巡视的导航装置及方法
Zhang et al. A novel INS/USBL integrated navigation scheme via factor graph optimization
Wang et al. Absolute navigation for Mars final approach using relative measurements of X-ray pulsars and Mars orbiter
CN107270937A (zh) 一种离线小波降噪快速初始对准方法
Xinlong et al. An autonomous navigation scheme based on geomagnetic and starlight for small satellites
CN105988129A (zh) 一种基于标量估计算法的ins/gnss组合导航方法
Fu et al. Unified all-Earth navigation mechanization and virtual polar region technology
Liu et al. Adaptive central difference Kalman filter with unknown measurement noise covariance and its application to airborne POS
Chen et al. Gravity gradient tensor eigendecomposition for spacecraft positioning
CN103017773B (zh) 一种基于天体表面特征和天然卫星路标的环绕段导航方法
Chen et al. Observability analysis for orbit determination using spaceborne gradiometer
Siouris Gravity modeling in aerospace applications
Kinatas et al. QUEST aided EKF for attitude and rate estimation using modified Rodrigues parameters
Pei et al. Autonomous orbit determination using epoch-differenced gravity gradients and starlight refraction
CN107702718A (zh) 一种基于瞬间可观测度模型的机载pos机动优化方法与装置
Liu et al. Absolute navigation and positioning of mars rover using gravity-aided odometry
Rosen Analysis of hybrid satellite-to-satellite tracking and quantum gravity gradiometry architecture for time-variable gravity sensing missions

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant