CN111125874B - 一种动平台高精度测轨预报方法 - Google Patents

Abstract

本发明涉及一种动平台高精度测轨预报方法：该方法首先建立轨道模型：空间目标高精度轨道运动摄动力包括地球质心引力、地球非球形引力摄动、大气阻尼摄动、光压摄动、日月引力摄动和潮汐摄动；弹道目标高精度轨道运动模型的摄动力包括地球质心引力、地球非球形引力摄动、大气阻尼摄动、火箭推力；然后对测量参数进行坐标变换，完成惯性坐标系与地固坐标系的转换；再进行力模型计算，将大气阻尼摄动、光压摄动、日月引力摄动和潮汐摄动等进行计算；最后使用改进的加权最小二乘微分算法，得到目标的速度和加速度参数，在此基础上分别利用弹道目标的短弧数据及空间目标的多天多圈数据采用高精度轨道确定算法进行定轨。

Description

一种动平台高精度测轨预报方法

技术领域

本发明涉及一种动平台高精度测轨预报方法。属于轨道测量技术领域。

背景技术

在实际的卫星飞行运动中，卫星的初始状态是未知的，并且其运动微分方程中的物理常量和模型本身也是带有误差的，从而导致积分的轨道与卫星真实轨道存在偏差。所以，为了获得尽可能接近真实轨道的积分轨道，必须对卫星进行跟踪观测，最后综合卫星跟踪的几何信息和卫星运动的动力信息来估计卫星初始状态及相关物理参数。

精密定轨是在低精度的参考轨道(简称初轨)的基础上，利用区域或全球跟踪站的观测数据对参考轨道予以改进。轨道改进的同时还可以根据需要解算整周模糊度、测站坐标、对流层延迟、地球自转参数、天线相位中心偏差、地球质心偏差等参数。

由于卫星的摄动力比较复杂，导致难以得到任意时刻卫星位置和速度。因此，在建立摄动力模型之后，必须求出摄动加速度对卫星坐标、速度以及力模型参数的偏导数。

在建立线性化的观测方程后，可根据实际情况选择合适的参数估计方法对卫星的初始状态及相关参数进行估值，再用估计的参数对卫星的运动方程积分即得到精密轨道。至此，卫星精密定轨问题就成了一个参数估计问题。

由于动平台在观测期间可能存在运动，为了保证精度需要在惯性坐标系下进行空间目标和弹道目标的测轨和预报，因此需要建立惯性坐标与地固坐标的时空转换关系。空间目标和弹道目标在运动过程中受多种摄动力的影响，为了得到高精度轨道参数都需要采用高精度轨道模型，在此基础上分别利用弹道目标的短弧数据及空间目标的多天多圈数据采用高精度轨道确定算法进行定轨。

发明内容

本发明所要解决的技术问题是针对上述现有技术提供一种算法简单，效率高，可靠性好，预报精度高的动平台精确测轨预报方法。

本发明解决上述问题所采用的技术方案为：一种动平台精确测轨预报方法，该方法首先建立轨道模型：空间目标高精度轨道运动摄动力包括地球质心引力、地球非球形引力摄动、大气阻尼摄动、光压摄动、日月引力摄动和潮汐摄动；弹道目标高精度轨道运动模型的摄动力包括地球质心引力、地球非球形引力摄动、大气阻尼摄动、火箭推力；然后进行坐标变换，完成惯性坐标系与地固坐标系的转换；再进行力模型计算，将大气阻尼摄动、光压摄动、日月引力摄动和潮汐摄动进行计算；最后使用改进的加权最小二乘微分算法，得到目标的速度和加速度参数，在此基础上分别利用弹道目标的短弧数据及空间目标的多天多圈数据采用高精度轨道确定算法进行定轨。

优选地，所述方法包含以下步骤：

步骤一、建立运动模型

包括空间目标高精度轨道运动模型和弹道目标高精度轨道运动模型，两种模型分列摄动力项；

步骤二、坐标转换和摄动力计算

2.1)坐标转换

惯性坐标与地固坐标的坐标转换计算公式如下：

其中：为目标在j2000历元平赤坐标系中的坐标，

为目标地固坐标系中的坐标，

P、N和GAST分别为岁差矩阵、章动矩阵和格林尼治恒星时；

2.2)力模型计算

2.2.1)地球引力

地球引力等于中心天体各部分质量引力作用之和，利用球谐函数展开将其表达成如下形式：

其中

GM为引力常数，r为地心向径，P_nm为缔合勒让德多项式，R_c为正常椭球的长半径，M_c为参考椭球子午圈曲率半径，δ_0m为地心纬度，φ和λ分别为探测器在星固坐标系下的纬度和经度，φ'和λ'分别为积分单元的经度和纬度，N_nm为归一化章动矩阵，C_nm和S_nm为归一化地球重力场模型系数，和/>即中心天体的引力场谐系数，采用EGM96地球非球形引力场模型；

2.2.2)日月引力摄动

第三体摄动对探测器相对于中心天体质心坐标下处在处质量为M_p的摄动天体运动产生的摄动加速度为：

其中

摄动天体其对探测器的直接引力：

摄动体对中心天体的质心引力加速度：

为探测器在中心天体质心坐标系中的坐标，

日月引力摄动的计算涉及日月星历的计算，采用DE405模型进行计算；

2.2.4)大气摄动

大气摄动的阻尼加速度由下式表达

其中C_d是大气阻尼系数，A为卫星有效横截面积，ρ为卫星所在位置的大气密度，和/>分别为卫星和大气的运动速度；

2.2.4)光压摄动

卫星受到的光压力可以近似为太阳光垂直照射一表面对表面产生的压力，即

其中C_R＝1+η，A为球形卫星截面积或卫星帆板的面积，为太阳相对卫星的位置矢量，太阳辐射压强度P_m与卫星到太阳的距离有关，在1个天文单位处太阳辐射压强为

P₀≈4.5605×10^-6N/m²；

2.2.5)潮汐摄动

只考虑固体潮，对于位置在处质量为M_p的摄动天体引起的固体潮，对应的谐系数可表示为

其中k_nm为积分式所取的阶数，M_p为摄动天体质量，M_c为地球质量，φ'和λ'分别为积分单元的经度和纬度；

步骤三、加权最小二乘微分改进

待估计的和大气阻尼参数β统称为状态量记作X，状态量满足的微分方程可以统一写成如下形式

X(t)＝f(X_k,t)

其中：X_k为需要估计的t_k时刻的状态量；

对弹道目标进行观测的随机采样数据称为观测量，记作Y，通过线性化转化，观测量与状态量的关系可表示为

y_i＝B_ix_k+V_i

其中V为测量误差，i＝1,...,l，l为观测资料的点数，B_i＝H_iΦ(t_k,t_i),H_i为测量矩阵；Φ(t_k,t_i)为状态转移矩阵；

则轨道根数的预测结果为：

其中W_i为对应y_i的权矩阵，可取为

与现有技术相比，本发明的优点在于：

1)本发明建立了动平台的预报模型，给出了坐标转换方程；对转换方程中的摄动力求解给出了计算公式，满足动平台轨道模型的高精度要求。

2)本发明设计采用加权最小二乘微分改进算法对方程进行求解，算法简单，实现便利。

附图说明

图1为本发明实例的原理框图

图2为本发明实例的高精度轨道模型。

图3为图1中实例的动平台高精度测轨道确定处理流程。

具体实施方式

以下结合附图实施例对本发明作进一步详细描述。

本实施例中的一种动平台高精度测轨预报方法，含以下步骤：

步骤1、建立运动模型

为了满足空间目标和弹道目标的高精度测轨预报需求，首先根据空间目标和弹道目标的运动特点建立运动模型，其中空间目标高精度轨道运动摄动力包括地球质心引力、地球非球形引力摄动、大气阻尼摄动、光压摄动、日月引力摄动和潮汐摄动；弹道目标高精度轨道运动模型的摄动力包括地球质心引力、地球非球形引力摄动、大气阻尼摄动、火箭推力，如图1所示。

步骤2、惯性坐标与地固坐标的转换和摄动力计算

步骤2.1、坐标转换

因为动平台测量所得参数与需要得到的轨道根数分属不同坐标系，故需要进行坐标变化，尤其涉及惯性坐标与地固坐标的坐标转换。

惯性坐标与地固坐标的坐标转换计算公式如下：

其中：为目标在j2000.0历元平赤坐标系中的坐标，

为目标地固坐标系中的坐标，

R_z为位置矢量的旋转角(旋转参数),

P、N和GAST分别为岁差矩阵、章动矩阵和格林尼治恒星时。

岁差矩阵采用经典三次旋转：

其中

θ_A＝2004.1917476"t-0.4269353"t²-0.0418251"t³-0.0000601"t⁴-0.0000001"t⁵

z_A＝-2.5976176"+2306.0803226"t+1.0947790"t²+0.0182273"t³+0.0000470"t⁴-0.0000003"t⁵

t表示起算历元时刻到目标历元时刻的懦略世纪数；

R_z、R_y为位置矢量的旋转角(旋转参数)

章动矩阵的计算方法为：

N＝R₁(-ε_A-Δε)R_z(-Δψ)R_x(ε_A)

其中

ε_A＝84381.406"-46.836769t-0.0001831t²+0.00200340t³-5.76×10^-7t⁴-4.34×10^-8t⁵

ε_A为黄赤交角，

Δε为交角章动，

ΔΨ为黄经章动，

R₁、R_z、R_x为位置矢量的旋转角(旋转参数)

IAU2000A岁差章动模型包含了678项日月章动和687项行星章动，计算表达式如下：

对于日月章动部分有：

其中A_i、A_i'、A_i”、A_i"'为黄经章动系数，B_i、B_i'、B_i"、B_i"'为交角章动系数，α_i只与5个日月幅角相关。对于行星章动部分有：

其中A_i、A_i”为黄经章动系数，B_i、B_i"为交角章动系数，α_i与日月及行星相关的14个幅角相关。格林尼治恒星时GAST可表示为：GAST＝GMST+EE

GMST，EE分别为格林尼治平恒星时和二分差，它们的计算方法如下：

GMST＝θ(UT1)+0".014506"+4612.15739966"t+1.39667721t²

-0.00009344"t³+0.00001882"t⁴

其中θ(UT1)为世界时，Δψ为黄经章动，ε_A为真黄赤交角，C_i为地球重力场模型系数，Ω为速度矢量。

对空间目标精密定轨来说，严格按照上述转换关系计算，对于弹道导弹可以不考虑岁差和章动，只考虑与地球自转相关的恒星时GAST即可。

步骤2.2、力模型计算

根据动平台测量得到的测量参数，经坐标变换等计算得到目标的初始轨道根数参数，再利用力学修正模型进行修正，从而提高轨道预报精度。

步骤2.2.1地球引力

地球引力包含地球质心引力及地球非球形引力摄动。对于自然界的大行星来说，其形状都不是严格的球形，内部的密度分布也是不均匀的，这使得中心天体对环绕型探测器的引力作用不严格等于其质心引力，而是等于中心天体各部分质量引力作用之和即引力位函数为：

其中

G为万有引力，

ρ(s)为质点所在平行圈半径，

其中为探测器的位置，/>为地球内部的某一点的位置，上述积分应当对整个中心天体体积进行积分。中心天体对探测器的引力作用可表示成：

其中gradV为引力位梯度，表示在径向、横向和法向的引力。

由于中心天体的质量分布是无法准确知道的，上式虽然严格却并不能直接用来计算中心天体的引力势，一般是通过无穷级数展开来表示，一种常用的方法是利用球谐函数展开将其表达成如下形式：

其中

GM为引力常数，r为地心向径，P_nm为缔合勒让德多项式，R_c为正常椭球的长半径，M_c为参考椭球子午圈曲率半径，δ_0m为地心纬度，φ和λ分别为探测器在星固坐标系下的纬度和经度，φ'和λ'分别为积分单元的经度和纬度，N_nm为归一化章动矩阵，C_nm和S_nm为归一化地球重力场模型系数，和/>即中心天体的引力场谐系数，球谐参数的计算与特定的重力场模型相关，这里采用EGM96地球非球形引力场模型。

步骤2.2.2日月引力摄动

除了中心天体的引力外，环绕型探测器还受到其它天体的引力作用，对于地球卫星来说主要的摄动体是太阳和月球。由于距离较远，摄动天体对环绕型探测器的引力作用可以考虑成质点引力，在中心天体质心坐标下处在处质量为M_p的摄动天体其对探测器的直接引力作用可以表示成：

其中为探测器在中心天体质心坐标系中的坐标。除了对探测器的引力作用外，摄动体对中心天体的质心还有一引力加速度：

故第三体摄动对探测器相对于中心天体质心运动产生的摄动加速度为：

日月引力摄动的计算涉及日月星历的计算，采用DE405模型进行计算。

步骤2.2.3大气摄动

严格的计算大气摄动需要考虑卫星表面各部分与大气的作用力，对于一般形状不规则的卫星来说这是很困难的，除此之外大气摄动还与大气状态相关，而地球大气目前还存在较多不确定因素，所以对于一般的卫星来说严格的计算大气摄动力比较困难。大气摄动占主导作用的是大气的阻尼作用，作为近似可以只考虑卫星受到的阻尼作用而忽略大气的升力作用，相应的阻尼加速度由下式表达

其中C_d是大气阻尼系数，A为卫星有效横截面积，ρ为卫星所在位置的大气密度，和/>分别为卫星和大气的运动速度。

大气密度随高度、经度、纬度、地方时、季节、年份、太阳辐射流量及地磁流量等因素的影响，如果可以获取太阳辐射流量和地磁指数等空间天气数据，可以采用考虑这些参数的大气密度模型，提高大气密度的计算精度。

步骤2.2.4光压摄动

卫星受到的光压力可以近似为太阳光垂直照射一表面对表面产生的压力，即

其中C_R＝1+η，A为球形卫星截面积或卫星帆板的面积，为太阳相对卫星的位置矢量，太阳辐射压强度P_m与卫星到太阳的距离有关，在1个天文单位处太阳辐射压强为

P₀≈4.5605×10^-6N/m²

对于一般的精度需求，上面的计算方法可以较好的反应光压摄动的主要影响，由于卫星有效截面及卫星表面的反射系数的计算都存在较大误差，在定轨计算中一般会把C_R作为待估参数。

步骤2.2.5潮汐摄动

由于中心天体不是刚体，第三体的引力作用除了前面介绍的直接摄动作用外还会引起中心天体形变从而使中心天体的引力位发生变化，这种摄动影响即潮汐摄动。潮汐摄动按中心天体不同部分对第三体引力的响应又分固体潮、海潮和大气潮，其中固体潮影响最大，对于我们的精度需求来说只需要考虑固体潮就够了，固体潮的影响可表示为一附加引力位，对于位置在处质量为M_p的摄动天体引起的固体潮，对应的谐系数可表示为

其中k_nm为积分式所取的阶数，M_p为摄动天体质量，M_c为地球质量，φ'和λ'分别为积分单元的经度和纬度。

空间目标精密定轨一般需要考虑地球质心引力、地球非球形引力摄动、大气摄动、光压摄动和潮汐摄动，弹道导弹被动段一般只需要考虑地球质心引力、地球非球形引力摄动及大气阻尼，而且还可以进行简化考虑。将所有力模型枚举项相加作为轨道根数的修正项，提供精确的轨道根数。

步骤3、加权最小二乘微分改进

轨道预报就是按照给定的初始状态矢量对卫星轨道方程进行积分，从而得到任意时刻卫星状态矢量的过程，即在建立了运动模型基础上，弹道导弹目标和空间目标的精密定轨问题即如何由各种观测资料估计弹道目标轨道量及一些参数的问题。精密定轨的实质就是利用观测数据，结合动力学模型获得卫星状态量的最佳估计值，待估计的和大气阻尼参数β统称为状态量记作X，它是7维向量，状态量满足的微分方程可以统一写成如下形式

其中F为n维向量函数，X_k为需要估计的t_k时刻的状态量。该方程的解可通过数值积分或构造分析解获得，一般可以表示成如下形式

X(t)＝f(X_k,t)

对弹道目标进行观测的随机采样数据(距离、方位和仰角)称为观测量，记作Y，它是m(m≥1)维向量，与状态量的关系可表示为

Y(t)＝G(X(t),t)+V

其中V为测量误差。一般我们有一系列的随机测量数据，于是有

Y_i＝G(X(t_i),t_i)+V_i＝G(f(X_k,t_i),t_i)+V_i,i＝1,...,l

其中l为观测资料的点数，通常有m×l＞n。由前面的分析可知，弹道目标定轨问题实际上归结为常微分方程边值问题，其残差为：

其中：n为迭代次数。

通过线性化，定轨对应的非线性问题转化为线性方程组的求解问题，这里将其改写为如下形式

y_i＝B_ix_k+V_i

其中

B_i＝H_iΦ(t_k,t_i)，H_i为测量矩阵；Φ(t_k,t_i)为状态转移矩阵；

一般假定V_i为测量噪声，均值为零，且有方差矩阵

方程的求解可以引用线性理论中的最优估计理论，这里采用加权最小二乘估计，即

其中W_i为y_i对应的权矩阵，可取为由/>给出

则轨道根数的预测结果为：

定轨计算流程如图2所示。

实例：

于2015年11月11日10时至18日10时共7天时间内对选定的7颗激光标校星(其中有两颗为双星)进行了连续观测。通过任务计划引导截获模式对7颗星所有过境弧段进行跟踪。雷达获取目标数据后，进行轨道改进，利用改进获得的精密轨道根数进行轨道预报以考核定轨及预报精度，卫星定轨精度如下：

从上面定轨结果来看，使用该算法不同高度卫星精度可达300m。

除上述实施例外，本发明还包括有其他实施方式，凡采用等同变换或者等效替换方式形成的技术方案，均应落入本发明权利要求的保护范围之内。

Claims (1)

1.一种动平台精确测轨预报方法，其特征在于：该方法首先建立轨道模型：空间目标高精度轨道运动摄动力包括地球质心引力、地球非球形引力摄动、大气阻尼摄动、光压摄动、日月引力摄动和潮汐摄动；弹道目标高精度轨道运动模型的摄动力包括地球质心引力、地球非球形引力摄动、大气阻尼摄动、火箭推力；然后进行坐标变换，完成惯性坐标系与地固坐标系的转换；再进行力模型计算，将大气阻尼摄动、光压摄动、日月引力摄动和潮汐摄动进行计算；最后使用改进的加权最小二乘微分算法，得到目标的速度和加速度参数，在此基础上分别利用弹道目标的短弧数据及空间目标的多天多圈数据采用高精度轨道确定算法进行定轨，所述方法包含以下步骤：

步骤一、建立运动模型

包括空间目标高精度轨道运动模型和弹道目标高精度轨道运动模型，两种模型分列摄动力项；

步骤二、坐标转换和摄动力计算

2.1)坐标转换

惯性坐标与地固坐标的坐标转换计算公式如下：

其中：为目标在j2000历元平赤坐标系中的坐标，

为目标地固坐标系中的坐标，

P、N和GAST分别为岁差矩阵、章动矩阵和格林尼治恒星时；

2.2)力模型计算

2.2.1)地球引力

地球引力等于中心天体各部分质量引力作用之和，利用球谐函数展开将其表达成如下形式：

其中

GM为引力常数，r为地心向径，P_nm为缔合勒让德多项式，R_c为正常椭球的长半径，M_c为参考椭球子午圈曲率半径，δ_0m为地心纬度，φ和λ分别为探测器在星固坐标系下的纬度和经度，φ'和λ'分别为积分单元的经度和纬度，N_nm为归一化章动矩阵，C_nm和S_nm为归一化地球重力场模型系数，和/>即中心天体的引力场谐系数，采用EGM96地球非球形引力场模型；

2.2.2)日月引力摄动

第三体摄动对探测器相对于中心天体质心坐标下处在处质量为M_p的摄动天体运动产生的摄动加速度为：

其中

摄动天体其对探测器的直接引力：

摄动体对中心天体的质心引力加速度：

为探测器在中心天体质心坐标系中的坐标，

日月引力摄动的计算涉及日月星历的计算，采用DE405模型进行计算；

2.2.4)大气摄动

大气摄动的阻尼加速度由下式表达

其中C_d是大气阻尼系数，A为卫星有效横截面积，ρ为卫星所在位置的大气密度，和/>分别为卫星和大气的运动速度；

2.2.4)光压摄动

卫星受到的光压力近似为太阳光垂直照射一表面对表面产生的压力，即

其中C_R＝1+η，A为球形卫星截面积或卫星帆板的面积，为太阳相对卫星的位置矢量，太阳辐射压强度P_m与卫星到太阳的距离有关，在1个天文单位处太阳辐射压强为

P₀≈4.5605×10^-6N/m²；

2.2.5)潮汐摄动

只考虑固体潮，对于位置在处质量为M_p的摄动天体引起的固体潮，对应的谐系数表示为

其中k_nm为积分式所取的阶数，M_p为摄动天体质量，M_c为地球质量，φ'和λ'分别为积分单元的经度和纬度；

步骤三、加权最小二乘微分改进

待估计的和大气阻尼参数β统称为状态量记作X，状态量满足的微分方程统一写成如下形式

X(t)＝f(X_k,t)

其中：X_k为需要估计的t_k时刻的状态量；

对弹道目标进行观测的随机采样数据称为观测量，记作Y，通过线性化转化，观测量与状态量的关系表示为

y_i＝B_ix_k+V_i

其中V为测量误差，i＝1,...,l_，l为观测资料的点数，B_i＝H_iΦ(t_k,t_i),H_i为测量矩阵；Φ(t_k,t_i)为状态转移矩阵；

则轨道根数的预测结果为：

其中W_i为对应y_i的权矩阵，取为