本发明属于低轨卫星精密定轨技术领域，涉及一种低轨卫星大气阻力摄动建模和计算方法。
背景技术：
卫星技术的广泛应用，对卫星轨道确定精度提出了越来越高的要求。卫星定轨精度主要取决于两方面的因素，一是测量数据的质量，二是采用的定轨策略和方法。在高精度测量系统中,卫星动力学模型精度往往难以满足精度需求，需要在精密定轨过程中增加估计大气阻力因子、太阳光压系数，或引入周期性的经验力，甚至采用简化动力学方法来补偿力模型的不足；另一方面，当跟踪设备误差较大时,又需要借助高精度的卫星动力学模型，在定轨过程中增加估计测量数据的设备系统差，对测量设备进行标校。因此，实际定轨中必须综合考虑测量精度和动力学模型精度。
低轨(lowearthorbit,leo)卫星轨道运动受大气阻力摄动影响较大,因此在leo精密定轨过程中,通常需要求解的参数包括卫星的初始位置和速度、大气阻力因子cd，以及跟测站弧段相关的参数系统差。式(1)给出了传统大气阻力摄动加速度计算模型，主要参数包括大气阻力因子cd，大气密度ρ、迎风面积a、卫星质量m、以及相对速度v。传统的定轨方法中，全弧段一般仅解算一个大气阻力因子cd，用以补偿大气阻力摄动计算误差。这种方法仅适用于弧段适中且精度要求不高的情况，但对高精度leo定轨任务由于目前的大气密度模型都不十分精确，特别是在太阳活动和地磁指数活动比较剧烈(即磁暴)的时候，大气阻力在整个弧段存在较大变化，用这种方法计算的大气阻力摄动精度较差，定轨精度受到很大限制。
以我国“xx四号”02星为例，该星位于200km高的近圆轨道上，大气阻力摄动是主要的定轨误差源。如果定轨弧段选取太长，则在积分运动方程时，大气阻力误差的累积效应将很明显；如果定轨弧段太短，则参加定轨的测量数据太少，同样难以精确定轨。因此，对该星的定轨弧长取为0.8～3天。“xx四号”02星采用gps数据和usb数据定轨，并对测站usb外测设备进行标校。表1中给出了卫星在第90圈和91圈的两组实测数据算例，算例a采用1.5天数据，算例b采用0.8天数据，定轨期间发生了磁暴。
表1部分站usb测量系统差和随机差标校值(90和91圈算例)
表1结果显示，两组算例给出的活动站和渭南站usb测量系统差存在着明显差异，表明磁暴的发生给轨道摄动(主要是大气阻力摄动)的高精度建模和计算带来了困难，从而影响了测量系统的标校。因此，大气阻力摄动的高精度建模和计算(特别是在发生磁暴的情况下)，对于低轨卫星的精密定轨和一些科学应用，具有重要意义。
技术实现要素：
本发明的目的是提供一种低轨卫星大气阻力摄动建模和计算方法，解决了现有技术中存在的单一大气阻力因子难以精确描述大气变化引起的大气阻力摄动计算误差的问题。
本发明所采用的技术方案是，一种低轨卫星大气阻力摄动建模和计算方法，具体按照以下步骤实施：
步骤1、建立卫星的box-wing模型；
步骤2、建立分段线性大气阻力摄动模型；
步骤3、根据步骤1建立的box-wing模型和步骤2建立的分段线性大气阻力摄动模型，计算t时刻折线形大气阻力摄动加速度。
本发明的特点还在于：
步骤1具体为，将卫星主体部分和太阳帆板部分简化为简单的几何模型盒和翼，即box-wing，卫星本体简化为长方体，太阳帆板简化为长方形的面板，根据地面标定的卫星轮廓的几何尺寸，卫星各表面的面积表示为：ai(i＝1，…，nface)。
步骤2具体为，将整个定轨弧段[t0，tf]等分为n个弧段，则有n 1个节点，设各节点的大气阻力因子初值为cd(ti)，(i＝1，2，…，n 1)，利用一次牛顿插值多项式，计算任意两个节点之间的折线型大气阻力因子cd(t)：
步骤3中t时刻折线形大气阻力摄动加速度为：
其中，ai为卫星第i个平面的面积；γi为卫星相对速度矢量v与第i个平面法向ni的夹角,ni由定义的卫星姿态确定；相对速度矢量
步骤3还包括分段线性大气阻力因子估计。
分段线性大气阻力因子估计具体为：采用折线形大气阻力因子，实际的待估状态向量为：
卫星运动方程：
非线性观测方程：
z＝h(y，t)(5)
精密定轨的过程就是在给定初值y0的条件下，利用最小二乘方法联合求解方程组(4)、(5)，得到新的状态估值
本发明的有益效果是：通过本发明一种低轨卫星大气阻力摄动建模和计算方法计算卫星的大气阻力摄动，能够克服传统的单一大气阻力因子难以精确描述大气变化引起的大气阻力摄动计算误差，有效提高定轨精度，为低轨卫星的精密定轨，特别是在磁暴等特殊条件下的精密定轨提供了技术支持。
附图说明
图1是本发明一种低轨卫星大气阻力摄动建模和计算方法的leo卫星定轨计算流程图；
图2是本发明一种低轨卫星大气阻力摄动建模和计算方法的leo卫星box-wing模型示意图。
图3是本发明一种低轨卫星大气阻力摄动建模和计算方法的折线形大气阻力因子示意图。
具体实施方式
下面结合附图和具体实施方式对本发明进行详细说明。
本发明一种低轨卫星大气阻力摄动建模和计算方法，具体按照以下步骤实施：
步骤1、建立卫星的box-wing模型；
低轨卫星大气阻力的计算与卫星的迎风面积相关，传统的卫星常数中一般给定一个参考迎风面积，然而，由于卫星在轨姿态存在变化，其真实迎风面积并不固定，基于卫星的结构分析，将卫星的主体部分和太阳帆板部分简化成简单的几何模型——盒和翼(“box-wing”)，如图2所示，即本体简化成长方体、太阳帆板简化为长方形的面板，因此，可以根据地面标定的卫星轮廓的几何尺寸，简单计算卫星各表面的面积ai(i＝1，…，nface)
步骤2、建立分段线性大气阻力摄动模型；
将整个定轨弧段[t0，tf]等分为n个弧段，则有n 1个节点，给定各节点的大气阻力因子初值cd(ti)，(i＝1，2，…，n 1)，利用一次牛顿插值多项式，计算任意两个节点之间的折线型大气阻力因子cd(t)，如图3所示：
这样在整个弧段上计算的大气阻力因子呈分段线性，分段线性的大气阻力因子比直线型的大气阻力因子能够更加准确地反映大气密度快速波动引起的大气阻力变化；
步骤3、根据步骤1建立的box-wing模型和步骤2建立的分段线性大气阻力摄动模型，计算t时刻折线形大气阻力摄动加速度；
t时刻折线形大气阻力摄动加速度为：
其中，ai为卫星第i个平面的面积；γi为卫星相对速度矢量v与第i个平面法向ni的夹角,ni由定义的卫星姿态确定；相对速度矢量
卫星定轨过程是对状态参数迭代改进的过程，大气阻力摄动计算是其中的主要步骤之一，为了更加精确地计算大气阻力摄动加速度、估计卫星轨道，一般在定轨过程中将各节点的各节点的大气阻力因子cd(ti)也作为待估参数进行估计，因此，精密定轨过程中的大气阻力摄动加速度计算过程实际上还包括了大气阻力因子cd(ti)的估计过程，采用折线形大气阻力因子，实际的待估状态向量为：
卫星运动方程：
非线性观测方程：
z＝h(y，t)(5)
精密定轨的过程就是在给定初值y0的条件下，利用最小二乘方法联合求解方程组(4)、(5)，得到新的状态估值
实施例
将本发明一种低轨卫星大气阻力摄动建模和计算方法具体与卫星的精密定轨过程相结合进行说明，如图1所示。
1、测量数据预处理
收集地面跟踪设备对低轨卫星的测量数据z，对选取的定轨弧段上的测量数据进行预处理，剔除测量数据中的粗差和野值得到定轨观测值
2、设置初始概略轨道y*、各节点的大气阻力因子cd,i：
卫星的初始概略轨道由卫星的初定轨过程获得，一般精度在1km左右甚至更差，卫星的初始概略轨道用作精密定轨迭代过程的卫星轨道状态参数的初始值，利用各节点cd,i的初值一般可设为2.0；
3、建立卫星轨道动力学模型
建立低轨卫星的轨道动力学模型，包括二体运动、地球非球形摄动、潮汐摄动、三体引力摄动、大气阻力摄动、太阳光压摄动、相对论效应摄动以及经验加速度摄动；其中大气阻力摄动按照第三节中方法建立，包括box-wing模型、分段线性阻力因子以及式(2)表示的各表面的阻力之和；
4、建立观测方程
根据观测数据类型建立相应的观测方程
5、逐一处理观测数据
(1)读入一条观测数据，记录观测值
(2)基于星轨道动力学模型和轨道初值积分到观测时刻ti，得到参考状态y(ti)；
(3)由观测方程计算参考状态处的理论观测值h[y(ti)]，并计算观测数据改进量
(4)计算相应的法方程系数
(5)重复(1)～(4)步，直到所有观测数据处理完毕，从而构成法方程：
6、轨道改进
利用最小二乘法对(4)式求解可解算出卫星轨道状态的改进量δy，则新的轨道状态和各节点的大气阻力因子cd,i即为y＝y* δy；
轨道改进的过程需要进行迭代，将y作为新的参考轨道重复2～5步直到收敛。
1.一种低轨卫星大气阻力摄动建模和计算方法，其特征在于，具体按照以下步骤实施：
步骤1、建立卫星的box-wing模型；
步骤2、建立分段线性大气阻力摄动模型；
步骤3、根据步骤1建立的box-wing模型和步骤2建立的分段线性大气阻力摄动模型，计算t时刻折线形大气阻力摄动加速度。
2.根据权利要求1所述的一种低轨卫星大气阻力摄动建模和计算方法，其特征在于，所述步骤1具体为，将卫星主体部分和太阳帆板部分简化为简单的几何模型盒和翼，即box-wing，卫星本体简化为长方体，太阳帆板简化为长方形的面板，根据地面标定的卫星轮廓的几何尺寸，卫星各表面的面积表示为：ai(i＝1，…，nface)。
3.根据权利要求1所述的一种低轨卫星大气阻力摄动建模和计算方法，其特征在于，所述步骤2具体为，将整个定轨弧段[t0，tf]等分为n个弧段，则有n 1个节点，设各节点的大气阻力因子初值为cd(ti)，(i＝1，2，…，n 1)，利用一次牛顿插值多项式，计算任意两个节点之间的折线型大气阻力因子cd(t)：
4.根据权利要求1所述的一种低轨卫星大气阻力摄动建模和计算方法，其特征在于，所述步骤3中t时刻折线形大气阻力摄动加速度为：
其中，ai为卫星第i个平面的面积；γi为卫星相对速度矢量v与第i个平面法向ni的夹角，ni由定义的卫星姿态确定；相对速度矢量
5.根据权利要求1所述的一种低轨卫星大气阻力摄动建模和计算方法，其特征在于，所述步骤3还包括分段线性大气阻力因子估计。
6.根据权利要求5所述的一种低轨卫星大气阻力摄动建模和计算方法，其特征在于，所述分段线性大气阻力因子估计具体为：采用折线形大气阻力因子，实际的待估状态向量为：
卫星运动方程：
非线性观测方程：
z＝h(y，t)(5)
精密定轨的过程就是在给定初值y0的条件下，利用最小二乘方法联合求解方程组(4)、(5)，得到新的状态估值
