一种高面质比航天器轨道动力学分析方法

技术领域

本发明涉及一种航天器轨道动力学分析方法，属于航天器轨道动力学领域。

背景技术

1957年10月4日，苏联成功地发射了第一颗人造卫星，开创了空间时代，随着人类 空间领域的不断发展，大量面向不同任务要求的航天器被送入太空。然而在人类进入空间 的初期，只是利用了近地的环地球轨道，完成较单一的使命。近几十年来，随着航天技术的 迅速发展和人类社会的不断进步，产生了各种各样的航天应用和研究领域，同时也产生了 各种各样的航天器，高面质比航天器就是其中的一种，如太阳帆、充气卫星以及目前正在研 发的芯片卫星等都属于高面质比航天器。

近些年来将太阳帆用于深空探测等任务的研究被广泛进行，相对于传统的航天器 利用化学燃料产生推力，太阳帆是利用太阳光子撞击帆面后发生动量交换的原理产生太阳 光压力推动其做各种轨道机动运行。因此，太阳帆可以源源不断的获得推力而不受化学燃 料的限制，同时由于太阳帆不需要携带大量的燃料，这样就可以减小发射质量，极大的降低 发射成本。除了被广泛的用于深空探测的研究，太阳帆也在不断的被应用于地球轨道任务 的研究，例如将太阳帆用于地球磁场尾迹探测的任务。由于太阳帆显著的高面质比特性，因 此可以利用太阳光压力实现轨道拱线的太阳同步进动，这样就可以被动的长时间处于地球 磁场尾迹中，进行长时间的观察，这些都说明研究高面质比飞行器有着极其重要的意义。

为了更加准确的理解这些航天器的在轨运行状态，需要对其轨道动力学重新进行 分析研究。与传统意义上的航天器不同，由于其显著高的面质比特性会使得太阳光压力和 大气阻力等摄动力对其轨道产生极其重要的影响，因此就需要重新认识该类航天器的轨 道。由于太阳光压力和大气阻力都是与面质比成比例的，因此在该类航天器的动力学建模 过程中必须考虑太阳光压力和大气阻力的作用效果。传统的航天器都是以开普勒轨道为基 础，并利用控制系统装置消除这些摄动力对轨道的影响。然而在高面质比航天器轨道设计 中并不是主动的抵消太阳光压力和大气阻力等摄动力对轨道的影响，而是主动利用这些摄 动力获得新的任务轨道。另一方面，传统的轨道动力学分析都是以高斯摄动方程为基础，然 而对于高面质比航天器，由于太阳光压力和大气阻力的复杂性，会使得高斯摄动方程只能 用于特殊情况下的轨道分析。当前对于高面质比航天器的轨道动力学分析方法主要是基于 高斯摄动方程的分析方法，该方法形式简单、意义明确，但是通常需要借助轨道平均技术对 其进行求解，这会使得计算精度降低同时随着摄动力的增多对于该方程的求解会变得十分 的困难。

发明内容

为了克服现有技术的不足，本发明提供一种基于哈密尔顿原理的高面质比航天器 动力学分析方法，能够提高计算精度和计算效率。

本发明解决其技术问题所采用的技术方案包括以下步骤：

(1)计算地球惯性坐标系下地球的引力势能其中， φ代表地心纬度，P_n代表第n阶勒让德多项式，R代表地球的半径，μ代表引力系数，J_n代表带 谐系数，r代表航天器所在位置处的地心距；

(2)计算太阳光压力其中代表太阳光压力，A代表航天器的横截面积，P代表太阳光压强度，代表太阳光线的单位矢 量，代表沿着表面正法线方向的单位矢量。ρ_a表示被吸收系数，ρ_d表示漫反射系数，ρ_s表示 完全反射系数；

(3)建立地球阴影坐标系原点o位于地球质心，轴沿太阳光线方向，轴 垂直于航天器所在轨道平面，轴由右手定则确定；若航天器在坐标系下的位置坐标 (x′,y′,z′)满足x′＞0且R代表地球半径，则判断航天器处于地球阴影中；

(4)计算大气阻力其中，C_D代表阻力系数，ρ代表大气密 度，代表航天器速度，A_drag代表航天器沿速度方向的横截面积；

(5)根据高面质比航天器所受到的动能T和势能V得到拉格朗日函数L＝T-V；选取 广义坐标q，求得相应的广义动量得到哈密尔顿函数其中p_i、q_i分别 表示广义坐标、动量的第i个分量；

(6)基于建立的动力学模型，分别仿真太阳光压和大气阻力对处于不同初始轨道 条件下的高面质比航天器轨道的影响，得到高面质比航天器在一个轨道周期内轨道要素的 变化，得出轨道要素变化与初始轨道条件之间的关系，进而分析太阳光压和大气阻力对不 同初始轨道条件下的高面质比航天器的影响；

(7)基于建立的动力学模型，仿真太阳光压和大气阻力对于不同初始轨道条件下 的高面质比航天器的偏心率、半长轴以及近地点俯角的影响，从仿真图中搜索满足偏心率 和半长轴在一个周期内总的变化为零而近地点俯角的变化等于太阳光线进动角度的线，从 而得到平衡轨道所需满足的初始轨道条件。

本发明的有益效果是：基于哈密尔顿原理对高面质比航天器进行了轨道动力学建 模，然后对其轨道特性进行了分析。在该过程中同时考虑了太阳光压、大气阻力以及地球扁 率三种摄动因素对航天器的影响。同时为了更加真实的反映太阳光压的作用效果，在模型 中还考虑了地球阴影，当航天器处于地球阴影中时其所受到的太阳光压力为零。与基于高 斯摄动方程的方法相比，本发明不需要进行轨道平均处理，适应范围更加广泛。最后本发明 所用到的方法也可以用来研究绕其他星体运行的高面质比航天器的轨道特性分析。

附图说明

图1是绕地运行的航天器在地心惯性坐标系下的示意图；

图2是太阳在黄道面内的位置示意图；

图3是地球阴影示意图；

图4是只考虑太阳光压时轨道半长轴在一个周期内变化的示意图；

图5是只考虑太阳光压时轨道偏心率在一个周期内变化的示意图；

图6是只考虑太阳光压时轨道近地点俯角在一个周期内变化的示意图；

图7是只考虑大气阻力时轨道半长轴在一个周期内变化的示意图；

图8是只考虑大气阻力时轨道偏心率在一个周期内变化的示意图；

图9是在太阳光压和大气阻力作用下的平衡轨道示意图；

图10是本发明的分析流程图。

具体实施方式

本发明解决其技术问题所采用的技术方案包括以下步骤：

(1)建立动力学模型，具体步骤如下：

步骤一：参考《远程火箭与卫星轨道力学基础》，可知在地球惯性坐标系下地球的 引力势能函数模型可以表示为其中U代表地球引力势场，φ 代表地心纬度，P_n代表第n阶勒让德多项式，R代表地球的半径，μ代表引力系数，J_n代表带谐 系数，r代表航天器所在位置处的地心距。

步骤二：参考BongWie.SolarSailAttitudeControlandDynamics,Part1 (JournalofGuidance,Control,andDynamics,Vol.27,No.4,July-August2004)，可以 建立太阳光压力模型。假设系数ρ_a表示被吸收系数，ρ_d表示漫反射系数，ρ_s表示完全反射系 数，这些参数的取值与航天器表面材料等性质相关。此时太阳光压力就可以表示为

其中代表太阳光压力，A代表 航天器的横截面积，P代表太阳光压强度，对于近地轨道航天器可以近似认为P＝4.65×10^-6N/m²，代表太阳光线的单位矢量，代表沿着表面正法线方向的单位矢量。

步骤三：建立地球阴影模型。首先建立地球阴影坐标系原点o位于地球质 心，轴沿太阳光线方向，轴垂直于航天器所在轨道平面，轴由右手定则确定。根据几 何关系，可以确定在该坐标系下地球阴影满足如下条件：(1)x′＞0；(2)其 中(x′,y′,z′)为航天器在坐标系下的位置坐标，R代表地球半径。根据此条件即可判 断航天器此时是否处于地球阴影中。

步骤四：建立大气阻力模型。以大气分子撞击卫星表面建立阻力模型，可以近似的 认为入射能量被完全吸收，产生的阻力为其中Drag代表大气阻 力，C_D代表阻力系数，ρ代表大气密度，代表航天器速度，A_drag代表航天器沿速度方向的横 截面积。

步骤五：建立高面质比航天器的轨道动力学模型。首先根据高面质比航天器所受 到的动能T和势能V得到拉格朗日函数L，其表达式为L＝T-V。然后选取广义坐标q，并求得相 应的广义动量最后就可以得到哈密尔顿函数其中p_i,q_i分别表示广 义坐标/动量的第i个分量。

需要注意的是，在仿真计算时首先需要将初始的轨道六要素转化为广义坐标，当 仿真完成后又需要将广义坐标转化为轨道六要素，从而可以对轨道特性进行分析。

(2)分析摄动力对高面质比航天器轨道特性的影响，具体步骤如下：

步骤一：分析太阳光压、大气阻力对于高面质比航天器轨道的影响。基于建立的动 力学模型，分别仿真太阳光压和大气阻力对处于不同初始轨道条件下的高面质比航天器轨 道的影响。本发明选择一个轨道周期作为积分间隔，从而得到高面质比航天器在一个轨道 周期内轨道要素的变化。对仿真图进行分析，可以得出轨道要素变化与初始轨道条件之间 的关系，进而可以分析太阳光压和大气阻力对不同初始轨道条件下的高面质比航天器的影 响。

步骤二：分析高面质比航天器的平衡轨道。本发明假设航天器的运行轨道位于黄 道面内，基于建立的动力学模型，仿真了太阳光压和大气阻力对于不同初始轨道条件下的 高面质比航天器的偏心率、半长轴以及近地点俯角的影响。对仿真图进行分析，从图中可以 搜索到满足偏心率和半长轴在一个周期内总的变化为零而近地点俯角的变化等于太阳光 线进动角度的线，从而得到平衡轨道所需满足的初始轨道条件。

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施 例。

首先建立高面质比航天器在各种摄动力作用下的轨道动力学模型，具体过程如 下：

如图1所示，参考《远程火箭与卫星轨道力学基础》可得在地球惯性坐标系下的地 球引力势能函数模型，表示形式如下：

其中U代表地球引力势能，φ代表地心纬度，P_n代表第n阶勒让德多项式，R代表地 球半径，μ表示引力系数，J_n表示带谐系数，r代表航天器所在位置处的地心距。

在本发明中关于地球扁率的描述只取到J₂项，那么地球的引力势函数就可以简化 为：

其中(x,y,z)表示航天器在地心惯性坐标系下的位置坐标。

如图2所示，建立太阳光压力模型。假设系数ρ_a表示被吸收系数，ρ_d表示漫反射系 数，ρ_s表示完全反射系数，这些参数的取值与航天器表面材料等性质相关。此时太阳光压 力就可以表示为：

其中代表太阳光压力，A代表航天器的横截面积，P代表太阳光压强度，对于近 地轨道航天器可以近似认为P＝4.65×10^-6N/m²，表示太阳光线的单位矢量，表示沿着 表面正法线方向的单位矢量。

航天器表面对太阳光的反射比较复杂，因此在讨论太阳光压力对航天器轨道的影 响时，可以近似地认为太阳光压力的方向与太阳光的入射方向一致。同时本发明中假设航 天器垂直于太阳光线的面积保持不变，对上述表达式进行化简。于是就可以得到太阳光压 力的近似表达式：

F_SRP＝p_SRc_RA

其中p_SR为太阳光压力，在地球轨道附近约等于4.56×10^-6N/m²，c_R为反射系数，A为 航天器的横截面积。

因此就可以得到如下表达式：

其中代表太阳光压力，分别代表沿地心惯性系三个坐标轴 x,y,z正方向的单位矢量，λ代表太阳光线与地心惯性坐标系x轴之间的夹角，ε代表黄赤交 角。

则太阳光压在地心惯性坐标系下的三个分量分别为：

其中F_SRPx,F_SRPy,F_SRPz分别代表太阳光压力沿地心惯性坐标系x,y,z轴方向的分量。

如图3所示，建立地球阴影模型。假设太阳位于无穷远处，那么就可以忽略太阳的 视差，因此阴影区域就为一个圆柱形，它的半径等于地球的半径。首先建立地球阴影坐标系 原点o位于地球质心，轴沿太阳光线方向，轴垂直于轨道平面，轴由右手定则 确定。根据几何关系，可知在该坐标系下地球阴影满足如下条件：

(1)x′＞0；

(2)

其中(x′,y′,z′)为航天器在坐标系下的位置坐标，R为地球半径。

因此，就可以得到下面的表达式：

其中ε代表黄赤交角，(x,y,z)代表航天器在地心惯性坐标系中的坐标，λ代表太阳 光线与地心惯性坐标系x轴之间的夹角。

根据上面的坐标之间的转换，地心惯性坐标系下的地球阴影区域的表达式为：

接着，我们建立大气阻力模型。以大气分子撞击卫星表面建立阻力模型，可以近似 的认为入射能量被完全吸收，产生的阻力为：

其中Drag代表大气阻力，C_D代表阻力系数，ρ代表大气密度，代表航天器的速度， A_drag代表航天器沿速度方向的横截面积，本发明中假设沿速度方向的横截面积近似于航天 器的横截面积即A_drag＝A。在地心惯性坐标系中，大气阻力沿三个轴方向的分量可以表示 为：

其中D_x,D_y,D_z分别代表大气阻力沿地心惯性坐标系x,y,z轴方向的分量，v_x,v_y,v_z代表航天器速度沿地心惯性坐标系x,y,z轴方向的分量。在本发明中采用指数形式的大 气密度模型即其中h₀表示参考高度，ρ₀表示参考高度处的大气密度，H表示 标称高度。这里取h₀＝600km,ρ₀＝1.454×10^-13kg/m³，H＝71.835km。

得到上述几个模型后，我们就可以建立高面质比航天器的轨道动力学模型。首先 根据高面质比航天器的动能T和势能V得到拉格朗日函数L，其表达式为：

L＝T-V

然后选取航天器在地心惯性系中的位置坐标(x,y,z)为广义坐标q，即q＝(x,y,z )^T，并求得相应的广义动量：

其中q,p分别代表系统的广义坐标和广义动量。

最后就可以得到整个系统的哈密尔顿函数：

其中H表示哈密尔顿函数，q_i,p_i分别表示广义坐标和广义动量的第i个分量。

经过推导得到T,V,L,H的具体表达式如下：

其中(x,y,z)代表航天器此时在地心惯性坐标系下的三个位置坐标分量，m代表航 天器的质量，R为地球半径。

将哈密尔顿函数分别对广义坐标和广义动量求偏导，可得正则方程为

其中F_i表示系统中的非保守力。和具体表达式如下：

其中F_SRPx，F_SRPy，F_SRPz分别代表太阳光压力沿地心惯性坐标系x,y,z轴方向的分量， D_x，D_y，D_z分别代表大气阻力沿地心惯性坐标系x,y,z轴方向的分量。

在仿真开始时需要将轨道六要素转化为广义坐标。

在求解正则方程过程中所需的输入量为广义坐标动量(x,y,z,p_x,p_y,p_z)，但是通 常是使用经典的轨道六要素描述航天器在轨运行的状态，即半长轴(a)或近地点高度(h_p)、 轨道倾角(i)、升交点赤经(Ω)、近地点俯角(ω)、偏心率(e)、真近点角(f)或偏近点角(E) 或平近点角(M)。

当给定航天器的轨道六要素，则由轨道力学知识可知：

其中r代表航天器此时与地心的距离，a代表半长轴，h_p代表轨道近地点高度，e代 表偏心率，f代表真近点角。

已知在地心惯性坐标系下航天器的位置矢径因此通过坐标转换 就可以得到航天器在地心惯性坐标下的三个位置坐标分量x,y,z：

其中代表绕x轴的坐标旋转矩阵，代表绕y轴 的坐标旋转矩阵，代表绕z轴的坐标旋转矩阵，i代表轨道倾角， Ω代表真近点角，ω代表近地点俯角。将上式进行化简，就可以得到如下分量形式：

x＝r(cosΩcos(ω+f)-sinΩsin(ω+f)cosi)

y＝r(sinΩcos(ω+f)+cosΩsin(ω+f)cosi)

z＝r(sin(ω+f)sini)

由轨道力学知识可知航天器的速度可以表示为如下形式：

其中vr表示航天器的径向速度分量，vf表示切向速度分量。

通过与前边相同的坐标转换就可以得到航天器速度矢量在地心惯性坐标系下的 三个分量：

进一步化简可以得到：

于是，三个广义动量(px,py,pz)就可以表示为：

当仿真完成后又需要将广义坐标转化为轨道六要素，从而可以对轨道特性进行分 析。

由广义坐标/动量可以表示得到航天器此时的位置矢量和速度矢量：

其中表示航天器的位置矢量，表示航天器的速度矢量。

按照轨道力学知识，在得到航天器的位置矢量和速度矢量后，可以按下面的步 骤可以分别计算出所有的轨道要素：

其中μ代表引力系数。

上面关于真近点角的求解存在奇异性，因此可以按下面的方法解决:如果则0＜f＜π，否则π＜f＜2π。

由位置矢量和速度矢量可以计算得到航天器此时的角动量：

因此根据角动量h就可以得到轨道倾角：

其中代表沿地心惯性坐标系z轴方向的单位矢量，h_z表示轨道角动量在地心惯 性坐标系z轴方向的分量。

需要注意的是，本发明中不考虑逆行轨道即0≤i≤π，因此由上面的公式计算得到 的轨道倾角就不会发生奇异。

由角动量可以计算得到轨道节线矢量:

其中代表轨道节线矢量，k代表沿地心惯性坐标系z轴方向的单位矢量，代表 轨道角动量矢量。

根据轨道节线矢量就可以计算得到升交点赤经：

其中代表沿地心惯性坐标系x轴方向的单位矢量，n_x代表轨道节线矢量在地心惯 性坐标系x轴方向的分量。

类似的，求解反三角函数出现的角度奇异性可以按下面的方法解决：如果 否则π＜f＜2π。接着可以得到近地点俯角：

如果则0＜ω＜π，否则π＜ω＜2π。

分析摄动力对高面质比航天器轨道特性的影响，具体步骤如下：

分析实例一，为了研究太阳光压力和大气阻力对高面质比航天器的影响效果，基 于建立的动力学模型，分别仿真太阳光压和大气阻力对处于不同初始轨道条件下的高面质 比航天器轨道的影响。在仿真过程中选择的航天器的初始轨道要素为：近地点高度h_p＝ 669.42km，轨道倾角i＝23.4°，偏心率e的变化范围为0.01～0.8，近地点俯角的变化范围为 0～2π同时航天器的面质比为A/m＝32.24m²/kg，仿真结果如图4-8所示，图4-6表示的是在 单独太阳光压力作用下，轨道半长轴、偏心率以及近地点俯角在一个轨道周期的变化。对图 4-6分析可知当高面质比航天器的初始近地点俯角在0°-180°之间时(初始太阳光线角为 0°)半长轴增大，偏心率减小；反之，当初始近地点俯角在180°-360°之间时半长轴减小，偏 心率增大。这是因为当初始近地点俯角在0°-180°之间时航天器在远离太阳的半个轨道获 得的能量大于靠近太阳时消耗的能量，因此航天器的总能量增加，半长轴增加；反之，当初 始近地点俯角在180°-360°之间时航天器的总能量是减小的，所以半长轴减小。图7-8表示 的是在单独大气阻力作用下，轨道半长轴和偏心率在一个周期内的变化。对图7-8分析可知 无论初始条件如何变化，轨道半长轴和偏心率都是减小的，这是因为大气阻力为耗散力，因 此高面质比航天器的能量总是在变小，所以轨道半长轴和偏心率一直在减小。

分析实例二，由前边的分析可知，当航天器受到太阳光压和大气阻力摄动时轨道 会发生很大的变化。为了保证航天器能够在轨长时间运行，就必须寻找一种平衡轨道，即在 这两种摄动力作用下的航天器轨道相对于太阳光线保持不变。本发明假设航天器的运行轨 道位于黄道面内，此时根据前边的分析可知这两种摄动力只会改变轨道面内要素：偏心率、 半长轴以及近地点俯角。基于建立的动力学模型，仿真了太阳光压和大气阻力对于不同初 始轨道条件下的高面质比航天器的偏心率、半长轴以及近地点俯角的影响，仿真结果如图9 所示。对图9进行分析可知，在图中存在三个曲面分别代表：Δe＝0即偏心率在一个轨道周 期内的总变化为0，Δa＝0即半长轴在一个轨道周期内的总的变化为0，Δω＝Δλ即近地点 俯角在一个轨道周期的进动角度等于太阳光线在一个周期的进动角度。于是可得处于这三 个曲面交线上的点即表示满足平衡轨道要求的轨道要素。当高面质比航天器的初始轨道要 素位于这条线上时，其轨道半长轴、偏心率在一个周期内的总的变化为0，仅有近地点俯角 随着太阳的进动而同步变化，因此轨道相对于太阳光线保持不变。