本发明提供一种基于大气阻力加速度测量的近地轨道确定方法，它涉及一种通过测量大气阻力加速度，对近地轨道飞行器进行轨道确定的方法，属于导航技术领域。
背景技术：
航天器的自动化、自主化、智能化是未来航天器发展的重要方向。航天器的自主运行和自我管理，一方面可以大大减少对地面系统的依赖，在降低任务运营成本的同时，提高作业效率；另一方面可以在很大程度上增强和改善航天器的功能和性能。航天器自主导航是实现航天器自主化的关键技术之一。对于近地轨道运行的航天器，基于全球导航卫星系统(Global Navigation Satellite System，GNSS)的导航是当前研究较多的自主导航方式。由于具备高精度、全天候和实时三维定位等优点，GNSS系统，尤其是全球定位系统(Global Positioning System，GPS)，已广泛应用于各类低轨卫星任务(对地成像、海洋测高、重力场测量、编队飞行等)。但是，GNSS系统的运行依赖于昂贵的地面监控和维护设施；且GNSS的信号较弱，容易受到干扰和欺骗，无法摆脱电子武器的威胁。因此，寻找一种可以在航天器上独立使用、不受外界电磁环境干扰、成本低廉易于实现、不依赖于其他航天器提供的定轨信息的自主导航方式，在经济、军事、科学研究等方面都具有重要的意义，例如用于军事侦察卫星的自主导航、类地球行星的探测。
由于在低地球轨道运行的航天器会受到大气的阻力的影响，因此利用静电加速度计可以对航天器所受到的大气阻力加速度进行精确的测量。静电加速度计是空间引力物理实验、卫星重力测量的关键载荷之一，其主要功能是测量航天器非保守力或者引力梯度信息。在各种类型的静电加速度计中，基于电容传感技术和静电反馈控制技术的静电加速度计具有分辨率高、灵敏度高、同时测量多个自由度加速度等优势。搭载其的航天器当受到大气阻力和太阳光压以及其它非保守力时，该静电加速度计与航天器会发生相对位移，根据位移变化结合相关技术来测量航天器所受到的非保守力。美国、欧洲目前已先后发射了CHAMP、GRACE、GOCE等地球重力场测量卫星，分别搭载了STAR、SuperSTAR、GRADIO静电加速度计，成功对全球重力场进行了测量。其中，GRADIO静电加速度计达到了的噪声密度水平，完全具备对大气阻力加速度进行测量的能力。
本专利提出的利用静电加速度计测量得到的大气阻力加速度进行轨道确定的方法是一种新型的航天器自主导航方式。这种新型的定轨方法只利用航天器自身获取到的大气阻力加速度信息和姿态信息，无需外界的辅助。由于在使用该方法进行定轨的过程中，只用到了大气阻力加速度矢量的方向而没有使用模值，因此这种方法不依赖于具体的大气模型，选用不同的地球大气密度模型对定轨精度没有任何影响。这有效地规避了地球大气模型的不确定度对定轨精度造成的一系列影响。
在利用大气阻力加速度测量进行轨道确定的过程中，需要使用星敏感器对当前航天器的进行姿态确定。相比于太阳敏感器、磁强计、地平仪和陀螺仪等其他常见姿态测量设备而言，星敏感器姿态测量精度比较高，抗干扰能力比较强，是目前航天器上最主要的姿态测量仪器。
综上，基于大气阻力加速度测量的近地轨道确定方法具有巨大的研究与应用价值，它不依赖于GNSS卫星提供的定轨信息，通过将航天器自身采集的大气阻力加速度和姿态信息转化为当前的位置速度信息，实现了航天器的完全自主的轨道确定。
技术实现要素：
(一)发明目的：本发明创新性地利用航天器上搭载的三轴静电加速度计测得的大气阻力加速度，并借助星敏感器获取的航天器当前姿态信息，对低地球轨道运行的航天器进行完全自主的轨道确定，不依赖于具体的地球大气模型。与传统的基于GNSS的轨道确定方法相比，这种新型的轨道确定方法广泛适用于低轨航天器，如CubeSat、微纳卫星、对地成像卫星、气象卫星等，其需要的所有信息都可以由航天器自主获取，不依赖于外部对象，极大地提高了自主运行的航天器的抗干扰性，实现成本低廉，具有广泛的应用前景。
(二)技术方案
本发明是一种基于大气阻力加速度测量的近地轨道确定方法，其步骤如下：
步骤一：准备工作
地心惯性坐标系的定义：地心惯性坐标系的原点在地球中心，x轴在地球赤道平面内，指向春分点，z轴垂直于赤道平面指向北极，y轴按右手法则确定。
地心赤道旋转坐标系的定义：地心赤道旋转坐标系的原点在地球中心，x轴沿赤道平面于格林尼治(Greenwich)子午面的相交线，z轴垂直于赤道平面，y轴按右手法则确定。
航天器第二轨道坐标系的定义：航天器第二轨道坐标系的原点在航天器质心，z轴指向地心，y轴垂直于轨道平面且与轨道动量矩相反，x轴在轨道平面内指向前。该坐标系是对地三轴稳定航天器姿态的参考基准。
航天器本体坐标系的定义：航天器本体坐标系的固连在航天器本身，坐标原点位于航天器质心，x，y，z三个轴在设计过程中进行定义，对于三轴姿态稳定的航天器，在没有姿态误差的情况下，其本体坐标系的三个轴与第二轨道坐标系重合。
航天器上搭载的三轴加速度计的x，y，z三个测量轴分别与航天器本体坐标系的x，y，z三个轴重合，这样三轴加速度计测量得到的加速度分量值就是在航天器本体坐标系下的分量值，省去了一步坐标转换过程。三轴加速度计测量得到的是航天器在空间中所受到的非保守力，包括大气阻力、太阳光压力、航天器发动机推力等。因此，在利用大气阻力加速度测量进行定轨时，航天器应在近地轨道飞行且所有发动机处于关机状态，以保证航天器所受到的所有非保守力中，除了大气阻力以外的其他非保守力都可以看作是无穷小量而忽略不计。
星敏感器通过感光元件对星图的识别，可以获取当前航天器相对惯性空间的姿态信息，并用坐标转换矩阵将其表示出来。
轨道预报器是指已知航天器在某一时刻的位置和速度矢量{r，v}，对航天器任意时间间隔Δt之后的位置和速度矢量进行计算的一种工具。常见的轨道预报器有二体轨道预报器、高精度轨道预报器等。
步骤二：大气阻力加速度测量
通过三轴加速度计对航天器所受到的非保守力进行测量，获得航天器所受非保守力在航天器本体坐标系下的分量列阵：
上式中，abx表示航天器所受到的非保守力在航天器本体坐标系的x轴方向的分量，同样地aby，abz分别为航天器所受到的非保守力在航天器本体坐标系的y，z两个轴方向的分量。在第一次测量之后，每隔一段相等时间Δt(间隔时间远小于轨道周期)测量一次，从而获得n个非保守力加速度测量值{ab1，ab2，...，abn}，n≥4。
通过星敏感器的测量信息，获得航天器本体坐标系到地球惯性坐标系的坐标转换矩阵Lib，则在地心赤道惯性坐标系下，航天器所受非保守力的分量列阵如下：
上式中，aix表示航天器所受到的非保守力在地心赤道惯性坐标系的x轴方向的分量，同样地aiy，aiz分别为航天器所受到的非保守力在地心赤道惯性坐标系的y，z两个轴方向的分量。将测得的n个非保守力加速度测量值按照以上步骤转换到地心赤道惯性坐标系下，可得{ai1，ai2，...，ain}，n≥4。
步骤三：非保守力加速度矢量单位化
由于大气模型具有一定的不确定性，为避免使用大气模型造成较大的计算误差，这里将加速度计测得的航天器所受到的非保守力加速度单位化，仅利用其中的加速度方向信息，单位化的过程如下：
将地心赤道惯性坐标系下的n个非保守力加速度测量值分别按照公式(3)进行单位化可得应注意的是，由于大气随地球的转动效应，后续步骤需要对测量得到的带有上标的非保守力加速度单位矢量进行校正，从而得到航天器速度矢量的单位矢量并将其代入定轨方程进行计算。
步骤四：给出轨道确定初始值
在圆轨道假设下，使用即可确定基于大气阻力加速度测量的轨道确定初始值。首先计算真近点角差：
其中θ表示航天器轨道的真近点角，Δθ表示两个真近点角之间的差值，表示未经过校正的非保守力加速度单位矢量。
计算轨道半长轴：
其中地球引力常数为μ＝GM＝3.98×1014，G为万有引力常量，其值为6.67×10-11N·m2·kg-2，M为地球质量，其值为5.965×1024kg。
计算航天器在轨速度的模值：
由此可以得到轨道确定过程中航天器速度矢量的初值：
下面计算轨道动量矩的一个方向向量：
由此可以得到轨道确定过程中航天器位置矢量的初值：
其中符号“×”表示矢量的叉乘。
步骤五：构建轨道确定方程
利用前述步骤获取的在地心赤道惯性坐标系下表示的n个非保守力加速度测量值的单位向量构建轨道确定方程如下：
其中uin为通过修正非保守力加速度获得的航天器速度矢量的单位矢量，vn＝P(r1，v1，(n-1)Δt)，P为航天器轨道预报器(二体轨道预报器或者其他更高精度的轨道预报器)。令
则公式(10)可以改写为：
u＝H(r1，v1，Δt，n) (12)
其中H是关于r1，v1，Δt，n的函数。
步骤六：求解轨道确定方程
下面使用数值方法求解方程(12)，由于步骤四已经给出相对准确的轨道确定初值，因此可以使用牛顿法等数值方法对该方程进行求解。
首先，利用步骤三中得到的{r1，v1}，通过二体轨道预报器计算Δt，2Δt，...，(n-1)Δt时间间隔后的航天器轨道位置和速度矢量{r2，v2}，{r3，v3}...{rn，vn}。利用这组轨道数据计算地球大气自转项，修正航天器受到的非保守力加速度，从而得到航天器速度矢量的单位矢量如下：
其中ωEi为当前时刻下在地心惯性坐标系下表示的地球自转角速度矢量，可以通过下面的公式得到：
ωEi＝LifωE (14)
其中ωE＝[0 0 ωE]T为在地心赤道旋转坐标系下表示的地球自转角速度矢量，ωE为地球自转角速度，其值为7.29×10-5rad/s。
其次，将方程组(12)代入牛顿法公式进行迭代计算。使用牛顿法对该方程进行求解的关键在于对H函数导数的求解。H函数的导数往往无法解析求解，因此这里以二体轨道预报器为例，提供一种求解H导数的方法：首先将{r1，v1}整合为一个矢量
然后在p矢量的六个分量上分别加一定的偏移量，得到：
之后分别以p1，p11，p12，p13，p14，p15，p16为初始轨道状态，利用二体轨道预报器预报Δt，2Δt，...，(n-1)Δt时间间隔后航天器的速度矢量，相应地可以得到{v1，v2，v3，...，vn}，{v11，v21，v31，...，vn1}，{v12，v22，v32，...，vn2}，{v13，v23，v33，...，vn3}，{v14，v24，v34，...，vn4}，{v15，v25，v35，...，vn5}，{v16，v26，v36，...，vn6}。那么由此可以计算中的各元素如下：
则构建迭代式如下：
其中表示的伪逆。迭代终止条件是pk+1矢量中的位置分量与pk矢量中的位置分量差的模值足够小于1m。
得到上述pk+1矢量后，将该矢量中的位置和速度分量重新设为{r1，v1}，代回到公式(13)重新修正非保守力加速度单位矢量，再重新利用牛顿法迭代求解方程(12)。如此重复计算，迭代终值条件是第(j+1)步得到的r1矢量和第j步得到的r1矢量差的模值小于1m。
通过以上步骤，提出了一种基于大气阻力加速度测量的近地轨道确定方法；该方法仅使用航天器自身搭载的三轴加速度计和星敏感器，获取几组航天器在轨运行时所受到的非保守力加速度值和航天器本体相对于地心惯性坐标系的坐标转换矩阵，利用该坐标转换矩阵将航天器所受非保守力加速度转化到地心惯性坐标系下表示，在圆轨道假设下获取定轨迭代计算初值，考虑到大气转动效应对测得的非保守力加速度单位矢量进行修正得到航天器速度矢量的单位矢量，最后利用轨道预报器和数值方法对定轨方程进行求解从而实现高精度的轨道确定。
该方法可以在航天器上独立使用、不受外界电磁环境干扰、成本低廉易于实现、不依赖于高精度地球大气模型、不依赖于其他航天器提供的定轨信息，在经济、军事、科学研究等方面都具有重要的意义，可用于军事侦察卫星的自主导航、类地球行星的探测。与传统的基于GNSS的轨道确定方法相比，这种新型的轨道确定方法广泛适用于低轨航天器，如CubeSat、微纳卫星、对地成像卫星、气象卫星等。
(三)优点
本发明提供的一种基于大气阻力加速度测量的近地轨道确定方法，相较于传统的轨道确定方法的优点在于：
1)本发明创新性地应用了三轴静电加速度计获取的大气阻力加速度信息，实现了航天器完全自主定轨，避免了由于外界电磁环境收到干扰而产生的定轨误差。
2)本发明提出的定轨方法广泛适用于任何在近地轨道运行的航天器，如CubeSat、微纳卫星、对地成像卫星、气象卫星等。它不依赖于高精度的地球大气模型，有效地避免了由地球大气模型的不确定度在定轨过程中产生的误差。
3)本发明提出的定轨方法实现成本低廉，只需航天器自身搭载三轴加速度计和星敏感器即可实现高精度轨道确定，无需使用其他GNSS卫星提供的定轨服务，对于低成本航天器而言非常易于实现，具有广泛的应用前景。
附图说明
图1是本发明所述方法流程图
具体实施方式
下面将结合图1和技术方案对本发明的具体实施过程做进一步的详细说明。
本发明一种基于大气阻力加速度测量的近地轨道确定方法，见图1所示，其步骤如下：
步骤一：大气阻力加速度测量
通过三轴加速度计对航天器所受到的非保守力进行测量，获得航天器所受非保守力在航天器本体坐标系下的分量列阵：
上式中，abx表示航天器所受到的非保守力在航天器本体坐标系的x轴方向的分量，同样地aby，abz分别为航天器所受到的非保守力在航天器本体坐标系的y，z两个轴方向的分量。在第一次测量之后，每隔一段相等时间Δt(间隔时间远小于轨道周期)测量一次，从而获得n个非保守力加速度测量值{ab1，ab2，...，abn}，n≥4。
通过星敏感器的测量信息，获得航天器本体坐标系到地球惯性坐标系的坐标转换矩阵Lib，则在地心赤道惯性坐标系下，航天器所受非保守力的分量列阵如下：
上式中，aix表示航天器所受到的非保守力在地心赤道惯性坐标系的x轴方向的分量，同样地aiy，aiz分别为航天器所受到的非保守力在地心赤道惯性坐标系的y，z两个轴方向的分量。将测得的n个非保守力加速度测量值按照以上步骤转换到地心赤道惯性坐标系下，可得{ai1，ai2，...，ain}，n≥4。
该步骤如图1中第一个方框所示。
步骤二：非保守力加速度矢量单位化
由于大气模型具有一定的不确定性，为避免使用大气模型造成较大的计算误差，这里将加速度计测得的航天器所受到的非保守力加速度单位化，仅利用其中的加速度方向信息，单位化的过程如下：
将地心赤道惯性坐标系下的n个非保守力加速度测量值分别按照公式(3)进行单位化可得应注意的是，由于大气随地球的转动效应，后续步骤需要对测量得到的带有上标的非保守力加速度单位矢量进行校正，从而得到航天器速度矢量的单位矢量并将其代入定轨方程进行计算。
该步骤如图1中第二个方框所示。
步骤三：给出轨道确定初始值
在圆轨道假设下，使用即可确定基于大气阻力加速度测量的轨道确定初始值。首先计算真近点角差：
其中θ表示航天器轨道的真近点角，Δθ表示两个真近点角之间的差值，表示未经过校正的非保守力加速度单位矢量。
计算轨道半长轴：
其中地球引力常数为μ＝GM＝3.98×1014，G为万有引力常量，其值为6.67×10-11N·m2·kg-2，M为地球质量，其值为5.965×1024kg。
计算航天器在轨速度的模值：
由此可以得到轨道确定过程中航天器速度矢量的初值：
下面计算轨道动量矩的一个方向向量：
由此可以得到轨道确定过程中航天器位置矢量的初值：
其中符号“×”表示矢量的叉乘。
该步骤如图1中第三个方框所示。
步骤四：构建轨道确定方程
利用前述步骤获取的在地心赤道惯性坐标系下表示的n个非保守力加速度测量值的单位向量构建轨道确定方程如下：
其中uin为通过修正非保守力加速度获得的航天器速度矢量的单位矢量，vn＝P(r1，v1，(n-1)Δt)，P为航天器轨道预报器(二体轨道预报器或者其他更高精度的轨道预报器)。令
则公式(33)可以改写为：
u＝H(r1，v1，Δt，n) (35)
其中H是关于r1，v1，Δt，n的函数。
该步骤如图1中第四个方框所示。
步骤五：航天器速度单位矢量误差修正
利用步骤三中得到的{r1，v1}，通过二体轨道预报器计算Δt，2Δt，...，(n-1)Δt时间间隔后的航天器轨道位置和速度矢量{r2，v2}，{r3，v3}...{rn，vn}。利用这组轨道数据计算地球大气自转项，修正航天器受到的非保守力加速度，从而得到航天器速度矢量的单位矢量如下：
其中ωEi为当前时刻下在地心惯性坐标系下表示的地球自转角速度矢量，可以通过下面的公式得到：
ωEi＝LifωE (37)
其中ωE＝[0 0 ωE]T为在地心赤道旋转坐标系下表示的地球自转角速度矢量，ωE为地球自转角速度，其值为7.29×10-5rad/s。
该步骤如图1中第五个方框所示。
步骤六：牛顿法求解定轨方程
下面使用数值方法求解方程(35)，由于步骤三已经给出相对准确的轨道确定初值，因此可以使用牛顿法等数值方法对该方程进行求解。
其次，将方程组(35)代入牛顿法公式进行迭代计算。使用牛顿法对该方程进行求解的关键在于对H函数导数的求解。H函数的导数往往无法解析求解，因此这里以二体轨道预报器为例，提供一种求解H导数的方法：首先将{r1，v1}整合为一个矢量
然后在p矢量的六个分量上分别加一定的偏移量，得到：
之后分别以p1，p11，p12，p13，p14，p15，p16为初始轨道状态，利用二体轨道预报器预报Δt，2Δt，...，(n-1)Δt时间间隔后航天器的速度矢量，相应地可以得到{v1，v2，v3，...，vn}，{v11，v21，v31，...，vn1}，{v12，v22，v32，...，vn2}，{v13，v23，v33，...，vn3}，{v14，v24，v34，...，vn4}，{v15，v25，v35，...，vn5}，{v16，v26，v36，...，vn6}。那么由此可以计算中的各元素如下：
则构建迭代式如下：
其中表示的伪逆。迭代终止条件是pk+1矢量中的位置分量与pk矢量中的位置分量差的模值小于1m。
得到上述pk+1矢量后，将该矢量中的位置和速度分量重新设为{r1，v1}，代回到公式(36)重新修正非保守力加速度单位矢量，再重新利用牛顿法迭代求解方程(35)。如此重复计算，迭代终值条件是第(j+1)步得到的r1矢量和第j步得到的r1矢量差的模值小于1m。
该步骤如图1中第六个方框所示。
通过以上步骤，提出了一种基于大气阻力加速度测量的近地轨道确定方法；该方法仅使用航天器自身搭载的三轴加速度计和星敏感器，获取几组航天器在轨运行时所受到的非保守力加速度值和航天器本体相对于地心惯性坐标系的坐标转换矩阵，利用坐标转换矩阵将航天器所受非保守力加速度转化到地心惯性坐标系下表示，在圆轨道假设下获取定轨迭代计算初值，考虑到大气转动效应对测得的非保守力加速度单位矢量进行修正得到航天器速度矢量的单位矢量，最后利用轨道预报器和数值方法对定轨方程进行求解从而实现高精度的轨道确定。
该方法可以在航天器上独立使用、不受外界电磁环境干扰、成本低廉易于实现、不依赖于高精度地球大气模型、不依赖于其他航天器提供的定轨信息，在经济、军事、科学研究等方面都具有重要的意义，可用于军事侦察卫星的自主导航、类地球行星的探测。与传统的基于GNSS的轨道确定方法相比，这种新型的轨道确定方法广泛适用于低轨航天器，如CubeSat、微纳卫星、对地成像卫星、气象卫星等。