[发明专利]一种基于重力场梯度不变量的轨道要素估计方法

摘要

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

主权项

1.一种基于重力场梯度不变量的轨道要素估计方法，其特征在于：其步骤如下：步骤一：准备工作1)地球引力势函数在初步研究航天器运动规律时，把地球看作理想的球体，它对航天器产生的引力指向地心，其大小F与航天器的质量m成正比，而与地心至航天器的距离r的平方成反比：上式中，F为航天器所受地球引力，G为地球引力常量，m为航天器质量，r为地心至航天器的几何距离，F除以m就是引力加速度g，把两个常数合并μ＝GM，则得到下式：上式中，μ成为地球引力常数；此情况下地球引力场成为中心引力场，把引力加速度表示成引力势的梯度，则中心引力场的引力势函数是：实际的地球并不是球对称的，具有扁平度、梨形和赤道变形，因而引力加速度是距离r、纬度和经度λ的函数，而r，λ是在地球固连坐标系S_e中位置的球坐标，为了严格地而且方便地描述地球引力加速度的分布，把它表示成地球引力势函数U的梯度：上式中grad(*)为求梯度算子；以下给出，在只考虑地球非球形J₂项影响下的引力场表达式：上式中Re＝6378.1363km，为地球赤道半径；2)开普勒轨道根数在惯性空间观察航天器轨道，需要定义轨道要素；轨道要素，又称轨道根数，是决定开普勒轨道的运动特征的一组常数；在考虑轨道摄动后，轨道根数不再是常数，可以作为状态变量；轨道根数一共有6个，表示及含义如下：a：椭圆轨道半场轴；e：椭圆轨道偏心率；Ω：升交点赤经；航天器轨道由南向北穿越赤道的点称为升交点B；赤道平面内，由春分点向东转到升交点B的角度称为升交点赤经，有效范围是0至360度；i：轨道倾角；轨道平面与赤道平面的夹角，有效范围是0至180度；t_p：近地点时刻；航天器通过近地点的时刻；ω：近地点幅角；轨道平面内，有升交点转到近地点的角度，有效范围是0至360度；他们的作用分别是：轨道半长轴a和偏心率e确定轨道的大小和形状；升交点赤经Ω和轨道倾角i确定轨道平面在惯性空间的取向；近地点幅角ω确定拱线在轨道平面内的位置；纪元时刻的平近点角M(t₀)确定时间的起点；3)轨道摄动方程在摄动力作用下航天器的轨道不是开普勒轨道，其轨道要素是随时间变化的：a(t)，e(t)，Ω(t)，i(t)，t_p，ω(t)；以下给出轨道摄动方程：上式中，与前文已出现的符号具有相同含义；f_r表示摄动加速度在地心轨道坐标系的径向分量；f_u表示摄动加速度在地心轨道坐标系的横向分量；f_h表示摄动加速度在地心轨道坐标系的副法向分量；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)求解矩阵的特征值λ₁,λ₂,λ₃，特征值与r,的关系如下：可以验证λ₁，λ₂，λ₃有如下关系：λ₁+λ₂+λ₃＝0  (13)因此，三个特征值中只有两个独立量；步骤四：利用J2重力梯度矩阵特征值求解各测量历元r和1)计算r初值使用J2模型求解的重力梯度张量矩阵特征值(12)，和理想球体重力场模型特征值表达式(9)计算r初值；式(12)中的特征值λ₁和λ₂任选其一皆可，带入式(9)，反解r则有如下初值表达式：或2)计算初值由式(12)做如下计算：从式(16)可以得到的表达式：从式(17)可以得到的另一表达式：以上的两种表达式，任选其一便可；3)利用更新r利用式(12)中的λ₂表达式更新r，将r₀和带入反解r³项，表达式如下：4)利用r₁更新方法一：用式(12)中的三个特征值做如下计算：将r₁和带入上式(23)，反解上式中的项，如下：方法二：使用式(21)和式(22)做如下计算：将r₁和带入上式(25)，反解上式中的项，如下：步骤五：利用r和计算初始轨道要素1)计算半长轴和偏心率在所有计算历元的结果r中搜索最大值r_max和最小值r_min，则有以下关系：r_a＝r_max,r_p＝r_min  (27)上式中，下标a表示远地点，下标p表示近地点；则半长轴：偏心率：2)计算轨道倾角在此假设，通过步骤四所求得的值的正负号已得知，在所有计算历元的结果中搜索最大值，轨道倾角：3)计算近地点幅角和初始真近点角在2)中搜索到的r_min其对应的纬度记为则近地点幅角ω由下式计算：初始历元真近点角θ₀，由下式计算：上式中，为第一个历元所计算出的纬度；步骤六：轨道要素平滑将初始历元五个轨道要素作为待估参数，如下：x＝[a₀ e_x0 e_y0 i₀ u₀]  (33)上式中，下标“0”表示初始历元参数，e_x,e_y,u表达式如下：e_x＝ecosωe_y＝esinω  (34)u＝ω+θ待估参数x的初值，由步骤五给出；轨道动力学模型由式(6)描述，半通径由于J2项引起的摄动力加速度在地心轨道坐标系下的分解如下：结合式(6)与式(35)，以及由步骤五提供的初值，便构成了已知初值的微分方程组求解问题，则通过数值积分方法可求得任意k历元时刻的五个轨道要素，x_k＝[a_k e_xk e_yk i_k u_k]，则k历元时刻的观测方程如下：上式中，为步骤四中所求得的各历元卫星到地心的几何距离与纬度，Δr,表示计算误差，下标“k”表示第k个历元；简单地假设Δr,是一个高斯噪声，则k历元测量协方差矩阵如下：上式中，σ_λ为重力梯度不变量误差，其误差大小要根据测量仪器的测量噪声与J2模型误差大小选取；将各个历元的测量值z放置在一起写成如下形式：上式中，N为测量历元的个数；各个历元测量协方差矩阵写成以下形式：采用准西格玛最小二乘批处理滤波估计状态参数，其估计式如下：上式中，下标j表示迭代次数，H_j表示利用计算出的测量量z的计算值与敏度矩阵；式(40)中，求解H阵难度较大，将Z写成关于x＝[a₀ e_x0 e_y0 i₀ u₀]的隐函数表达，如下：Z＝f(a₀ e_x0 e_y0 i₀ u₀)  (41)通过普通解析法求偏导计算敏度矩阵H,计算量非常巨大，因此可以通过以下方法求解；在第j次迭代完毕，得到状态参数估计值则进行如下采样：上式中，L为x的维度，P₀为先验协方差矩阵，表示P₀的第i列平方根；先验协方差矩阵由下式给出：将式(42)通过式(6)积分，和式(36)传递，计算测量值如下：χ_i,j→γ_i,j,i＝0,...,L  (44)则敏度矩阵H可以近似通过下式得出：上式中，表示P₀第i个对角线元素的平方根,i＝1,2,3,4,5。