利用大气阻力的低能耗高连续低轨星座构型保持方法

摘要

针对传统低轨星座构型维持中频繁机动导致的寿命缩减以及性能下降问题，提出了一种利用大气阻力的低轨星座相对构型自主保持方法。首先通过理论分析明确了大气阻力和卫星轨道相位之间的关系，而后在卫星能源、姿态指向等约束下，设计了通过旋转帆板改变卫星迎风面积的星座构型自主保持方法，并实现了算法的仿真验证。仿真结果表明，通过改变帆板迎风面积，可使卫星相对相位偏差量控制在0.1 °以内，相对半长轴差控制在40 m以内，卫星一年内无需进行主动相位控制，整星寿命期间燃耗节省至少13.284 kg。该方法能有效降低卫星燃料消耗，增加卫星整体服务连续性。

关键词

低轨星座; 相位保持; 大气阻力; 太阳帆板

当前，伴随着众多低轨宽带通信星座计划的发布，国际太空市场异常火爆。协同控制技术是实现低轨卫星星座应用的关键技术，而星座构形控制是协同控制研究的重要内容。由于地球扁率、日月第三体摄动、光压摄动以及大气摄动等对低轨星座几何结构的影响，导致星座工作性能的下降。卫星轨道控制系统必须进行构型维持，以达成星座构型长期稳定的目的^［1］。然而在推力器点火进行轨道控制过程中，卫星需要暂停服务，进而短暂地影响星座整体工作性能。此外，频繁地进行主动构型控制需要卫星与地面站频繁的传达指令， 这不仅加重工作人员负担以及运营成本，也会影响星座的整体可靠性与寿命^［2］。

目前，以星链（Starlink）、一网（OneWeb）、铱星（Iridium）等为代表的大型低轨星座，均通过推力器变轨调节卫星半长轴，进而实现卫星的相位控制。Iridium先将卫星送入高度较低的入轨轨道，再使卫星同时抬高轨道半长轴并改变轨道倾角，在进入工作轨道之后，保持升交点赤经漂移速率。并在入轨后保持各个卫星的平半长轴与平倾角，使得卫星升交点赤经漂移速率相同，以保证长期的星座构型稳定，其单星控制频率约为27 d。OneWeb则是在进入工作轨道后，对各星交替采用升轨控制与降轨控制的方法进行相位控制，其单星维持频率约为15 d。Starlink采用升高或降低半长轴来维持星座控制，再通过卫星升轨以维持轨道高度。但是由于Starlink星座的轨道高度分布在300 km至1 000 km，受大气阻力影响较大，卫星的平均半长轴每日衰减量大约为10 m。此外，同样由于大气阻力影响，星座进行维持控制时振荡幅度较大，卫星降轨后需要马上进行升轨来维持高度，这导致星座的半长轴控制较为频繁^［3］。因此，需要研究一种能够达成星座构型自主稳定的方法。

卫星星座结构自稳定是指在一定的容许范围内，卫星不需要主动使用推力器进行轨道控制，并且能使星座几何结构长期稳定的一种状态^［4］。对于常见的步行者（Walker）星座，由于其中卫星的轨道高度、偏心率和倾角都相同，因此在地球扁率摄动影响下的长期轨道变化是相同的。蒋虎^［5］对低轨卫星主要摄动源进行了评估，并获得不同摄动力的量级。对于低轨星座而言，大气摄动是造成星座构型恶化的主要因素之一。黄卫东等^［6］首先分析了大气摄动对星座构形发散的影响，并建立了一种降低大气摄动构形影响的补偿方法。李玖阳等^［7］在低轨Walker星座中使用两次构型偏置的方法，减少了漂移率。刘奇等^［8］提出了一种基于极限环的高精度站位保持方法，能够在卫星定轨数据精度不足，数据采样间隔较大的情况下实现星座高精度站位保持。孙俞等^［9］则是通过分析卫星面质比与太阳活动以及大气阻力之间的关系，提出了一种针对星座长期构型维持的自适应绝对站位保持方法，有效减少了控制频率。

在大多数情况下，大气阻力通常被视作必须被克服的扰动，尤其是对于绝对构型维持策略来说，大气阻力必须通过推力器点火加以消除。但是在特定条件下，大气阻力也可以用来实现卫星轨道的控制。Omar^［10］通过基于仿真软件开发的数值积分器验证了气动力可以有效用于轨道维持，但是对于轨道机动的能力较弱。Dutoit^［11］验证了利用大气阻力进行低轨星座构型控制的理论可行性。近年来，大气阻力被广泛应用于低轨卫星飞行编队控制，Ivanov^［12］ 和Larbi^［13］对利用大气阻力进行编队构型维持的问题，研究了基于线性二次型调节器的相对运动控制算法以及前馈和反馈PID（比例、微分、积分）的组合控制律。Lambert等^［14］验证了利用大气阻力进行低轨卫星飞行编队构型维持的可行性与有效性。Traub等^［15］研究了新型空气动力学材料结合大气呼吸电推进系统进行大气摄动补偿方法。赵晓峰等^［16］针对伴飞卫星的长期稳定，提出了一种利用改变姿态，调整面质比来进行姿轨控一体化设计。针对低轨卫星以及同步轨道卫星，Kumar等^［17］提出了一种利用环境摄动，即大气阻力与太阳光压进行编队构型维持方法。而Sun等^［18］则提出了一种由六块空气动力板为执行器的姿态轨道耦合控制模型。Shao等^［19］基于五块空气动力板，提出了基于非线性李雅普诺夫的反馈控制律。不过由于编队飞行与星座维持的任务目标不同，对构型的要求也不同。针对星座构型保持问题，林西强等^［20］验证了气动星座控制的可行性，并研究了其最优控制问题。杨维廉^［21］研究了星座同轨道面内利用空气阻力进行半长轴控制的方法。向开恒等^［22］在利用大气阻力的基础上，利用卫星最大、最小弹道系数进行星座构型维持。Foster等^［23］提出了一种对部署在同一轨道上的无推进大型卫星星座进行差分阻力控制的方法，入轨时将卫星送进已提前偏置的相位，通过最优化最大与最小弹道系数的控制率，成功实现相位调整。尽管上述办法通过调整卫星姿态以及使用卫星面质比的极值作为控制量在理论上可行，但是对在低轨星座应用中广泛存在的帆板能源损耗、帆板转动能力以及三轴稳定卫星姿态指向的约束考虑不足。

本文首先对低轨道摄动环境进行了研究，分析了其轨道变化特点，并论述了选择相对构型作为低轨星座构型维持的合理性；此外，出于卫星服务连续性以及卫星能源供给的考量，研究了能源损耗，帆板转动能力以及卫星姿态约束的问题；之后推导了卫星帆板迎风面积，半长轴以及角速度之间的关系；在此基础上，建立了控制策略流程并实现算法仿真验证；最后结果表明，通过高频次，小幅度调整帆板迎风面积，改变卫星面质比，从而改变轨道漂移速率，能够实现较高精度的相位保持，并且有效降低了主动控制频率，增加了星座服务连续性。对低轨星座整体性能以及寿命提升具有十分重要的意义。

1 低轨道摄动分析

卫星运行在低轨道上存在4种摄动力的作用：地球形状非球形引力、日月引力、太阳光压引起的摄动以及大气阻力。在摄动力的持续作用下，卫星无法根据二体运动规律运行，其轨道的半长轴、偏心率、轨道倾角、升交点赤经也在不断变化，从而使卫星轨道产生漂移。因此，需要对低轨卫星进行摄动分析，并根据其轨道变化特点，对卫星的轨道位置进行保持控制。表1为经过96 h轨道预报后，轨道高度在1 000 km的航天器所受的各种轨道摄动数量级。表中LEO （low earth orbit） 表示低地球轨道，由表可见，对于轨道高度1 000 km以及以下的卫星，除地球引力以外，地球非球形J₂摄动与大气阻力为主要摄动源，摄动量级分别为10⁷与10⁵，其他摄动为小量，并且随着越来越多低轨卫星卫星轨道高度的进一步降低，大气阻力摄动会随着高度降低而产生显著效果，而日月第三体引力摄动则随着轨道高度降低而减弱^［24］。因此本文将在上述4种摄动力之中更进一步分析地球非球形摄动与大气阻力摄动对卫星轨道的影响。

表1  低轨卫星受摄动力加速度量级

Tab. 1  Magnitude of perturbed dynamic acceleration of LEO satellite ( m )

摄动力	加速度量级
二阶非球形引力	107
大气阻力摄动力	105
四阶非球形引力	103
六阶非球形引力	103
月球引力	102
太阳引力	102
固体潮汐摄动力	102
太阳光压摄动力	102

1.1　地球非球形摄动

由地球非球形导致的摄动通常是影响卫星轨道的主要摄动力，并且随着轨道高度降低，地球非球形摄动的影响增加。对于轨道高度1 000 km以下的低轨卫星，地球中心引力与地球非球形J₂项摄动对卫星轨道产生的影响最大。它们对卫星各轨道根数的长期变化率影响表示为

$$\left\{\begin{array}{l}\dot{a}=\mathrm{0}\\ \dot{e}=\mathrm{0}\\ \dot{i}=\mathrm{0}\\ \dot{\mathrm{\Omega }}=-\frac{\mathrm{3}{J}_{\mathrm{2}}{R}_{\mathrm{e}}^{\mathrm{2}}}{\mathrm{2}{p}^{\mathrm{2}}}n\mathrm{c}\mathrm{o}\mathrm{s}i\\ \dot{\lambda }=\frac{\mathrm{3}{J}_{\mathrm{2}}{R}_{\mathrm{e}}^{\mathrm{2}}}{\mathrm{2}{p}^{\mathrm{2}}}n\left[\left(\mathrm{2}-\frac{\mathrm{5}}{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}i\right)-\sqrt[]{\mathrm{1}-e^{\mathrm{2}}}·\left(\mathrm{1}-\frac{\mathrm{3}}{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}i\right)\right]\end{array}\right.$$

（1）

式中：a为半长轴；e为偏心率；i为轨道倾角；Ω为升交点赤经；$\lambda$为沿迹角；$p=a\left(\mathrm{1}-e^{\mathrm{2}}\right)$；$n=\sqrt[]{\frac{\mu }{a^{\mathrm{3}}}}$；R_e为地球半径；μ为地球引力常数。

因此地球引力与非球形J₂项摄动对卫星半长轴、偏心率以及倾角无长期变化影响，但是对升交点赤经与沿迹角的变化影响明显。

1.2　太阳光压摄动

太阳光压主要对卫星的偏心率产生影响，其影响大小主要由面质比决定。作用在航天器单位表面积的辐射压力为

$$\mathrm{d}F=p\mathrm{d}s\left|\mathrm{c}\mathrm{o}\mathrm{s}\alpha \right|\left[\left(\mathrm{1}-c\right){\mathit{u}}_{\mathrm{I}}-c\mathrm{\text{'}}{\mathit{u}}_{\mathrm{F}}\right]$$

（2）

式中：$p$为太阳光压强度；$\alpha$为太阳光入射角；${\mathit{u}}_{\mathrm{I}}$、${\mathit{u}}_{\mathrm{F}}$分别为入射光、反射光方向的单位矢量；$c$和$c\mathrm{\text{'}}$分别为表面吸收率和反射率。太阳光压导致的轨道变化对中高轨道卫星的轨道预报与轨道控制较为重要，但是太阳光压摄动虽然对低地球轨道（LEO）卫星影响可以忽略不计，光压摄动产生的数量级只相当于地球非球形引力J₈项产生的加速度的二阶小量。因此在进行低轨航天器的轨道设计时，一般可以忽略光压摄动。

1.3　日月第三体引力摄动

第三体引力摄动加速度为

$${\mathit{a}}_{\mathrm{T}}=-{\mu }_{\mathrm{s}\mathrm{l}}\left(\frac{\mathit{\Delta }}{{\mathit{\Delta }}^{\mathrm{3}}}+\frac{{\mathit{r}}_{\mathrm{s}\mathrm{l}}}{{\mathit{r}}_{\mathrm{s}\mathrm{l}}^{\mathrm{3}}}\right),\Delta =\mathit{r}-{\mathit{r}}_{\mathrm{s}\mathrm{l}}$$

（3）

式中：${\mu }_{\mathrm{s}\mathrm{l}}$表示第三体的引力系数；$\mathit{r}$和${\mathit{r}}_{\mathrm{s}\mathrm{l}}$分别表示航天器和第三体在J2000系中的位置矢量；Δ为航天器相对太阳/月球中心的位置矢量。

由于日月第三体引力的存在，地球表面会因弹性形变而产生固体潮与大气潮，也会对航天器轨道产生影响。其中最大的影响来自固体潮，但是随着航天器轨道降低，日、月引力摄动就越小，对于低轨卫星，固体潮的影响与月球引力基本上处于一个数量级，第三体引力会引发轨道短周期项变化，这种高频率小幅波动不会对轨道有长期影响，因此可以不考虑日、月引力的的短周期项影响。但是太阳引力可能会引起轨道倾角的长周期变化，对低轨卫星的轨道稳定性会产生一定影响。进行轨道控制时需要考虑轨道倾角的维持。

1.4　大气阻力摄动

在近地轨道上的卫星在长时间高速运行于高层大气中时，大气阻力的积累会导致轨道衰减。大气阻力加速度公式为

$$a_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}=-\frac{\mathrm{1}}{\mathrm{2}}{C}_{\mathrm{d}}\frac{S}{m}\rho v$$

（4）

式中：C_d为阻力系数；$\rho$为航天器所在高度大气密度；S为卫星的迎风面积；m为卫星的质量；v为卫星相对大气的速度。通过调整帆板指向，可改变卫星迎风面积S，进而改变卫星所受的大气阻力摄动。

2 低轨星座构型维持方法

近地轨道卫星主要受到非球形引力摄动与大气阻力摄动。轨道倾角主要受到太阳光压摄动影响，产生长期摄动项，幅值0.01 （°）∙年^－1，因为轨道倾角对观测无影响，故不进行轨道倾角控制。轨道偏心率主要受到太阳光压的短周期项影响，无长期项影响。由图1可见，对于LEO星座，大气阻力摄动对卫星的半长轴产生长期摄动项影响，造成轨道衰减，特别是当星座内卫星面质比不同时，半长轴衰减不一致，进而改变相位差，导致星座构型恶化。因此需要进行轨道半长轴维持，以确保寿命期内轨道高度与相对相位满足任务要求。

图1  一年中卫星半长轴演化情况

Fig. 1  One-year evolution of satellite semi-axis

2.1　构型维持策略选择

当前，对于星座构型维持一般有绝对构型维持和相对构型维持两种思路。在选择相应策略之前，需要对低轨环境下主要摄动项对整体星座构型的影响进行研究与分析。卫星间相对位置变化则反映了星座整体结构的空间几何的相对变化。下面主要给出星座的相对变化。

卫星的相对位置关系可以用卫星轨道半长轴、偏心率、倾角以及卫星相位和升交点位置关系等来描述。从第1.1节对地球非球形摄动的分析得知，地球扁率摄动不会引起卫星轨道半长轴、偏心率和倾角的长期变化，但会导致卫星相位和升交点赤经的长期变化，因此可以用卫星之间的相位和升交点赤经的变化来描述星座结构的空间变化。

由式（1）可知地球非球形摄动对卫星轨道的长期影响，对于星座中任意的两颗星i和j，由轨道摄动模型得地球扁率摄动引起的卫星轨道面之差长期变化率$\mathrm{\Delta }{\Omega }_{\mathrm{1}}$和相位差长期变化率$\mathrm{\Delta }{\lambda }_{\mathrm{1}}$为

$$\mathrm{\Delta }{\Omega }_{\mathrm{1}}=-\frac{\mathrm{3}{J}_{\mathrm{2}}}{\mathrm{2}}(\frac{n_{\mathrm{i}}\mathrm{c}\mathrm{o}\mathrm{s}{i}_{\mathrm{i}}}{p_{\mathrm{i}}^{\mathrm{2}}}-\frac{n_{\mathrm{j}}\mathrm{c}\mathrm{o}\mathrm{s}{i}_{\mathrm{j}}}{p_{\mathrm{j}}^{\mathrm{2}}})$$

（5）

$$\mathrm{\Delta }{\lambda }_{\mathrm{1}}=\frac{\mathrm{3}{J}_{\mathrm{2}}}{\mathrm{2}}\left\{\begin{array}{l}\frac{n_{\mathrm{i}}}{p_{\mathrm{j}}^{\mathrm{2}}}\left[(\mathrm{2}-\frac{\mathrm{5}}{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}{i}_{\mathrm{i}})-(\mathrm{1}-\frac{\mathrm{3}}{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}{i}_{\mathrm{i}})\sqrt[]{\mathrm{1}-e_{\mathrm{i}}^{\mathrm{2}}}\right]\\ -\frac{n_{\mathrm{i}}}{p_{\mathrm{j}}^{\mathrm{2}}}\left[(\mathrm{2}-\frac{\mathrm{5}}{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}{i}_{\mathrm{j}})-(\mathrm{1}-\frac{\mathrm{3}}{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}{i}_{\mathrm{j}})\sqrt[]{\mathrm{1}-e_{\mathrm{j}}^{\mathrm{2}}}\right]\end{array}\right\}$$

（6）

式（5）—（6）中：$\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{1}}$为两卫星的轨道面赤经差的长期变化率；$\mathrm{\Delta }{\lambda }_{\mathrm{1}}$为两卫星相位差的长期变化率，通过上述两个关键变量可以描述星座整体构型的变化情况。

由式（5）—（6）可见，对于卫星半长轴、偏心率和倾角均相同的星座，地球扁率摄动不会引起卫星之间的相位差和升交点赤经差的长期变化，这说明地球扁率的长期不会引起这类星座结构的空间几何变化。

因此，对轨道高度、偏心率和倾角都相同的卫星，由地球扁率摄动导致的升交点赤经变化量均相同，在对由同种轨道类型的卫星构成的星座，升交点赤经变化量带来的升交点赤经长期漂移会引起星座整体性漂移，对于区域星座，这会导致其服务区域的系统性漂移。但是对于全球星座的覆盖以及服务不会产生影响。此外，根据云朝明等^［25］对不同高度情况下的大型星座控制频率的研究，高度较低且数量较少的星座比较适合采用相对构型维持方法。相比之下，绝对构型保持虽然简单，可靠，但是任务指标与控制基准和最大允许偏差转化时并不完全等价，在实际构型保持中若按照某一固定允许误差作为判别依据，要么偏于保守，浪费燃料，要么构型已经破坏，却没做出相应的调整。加之绝对构型控制策略的缺点较为明显，对发射入轨精度要求较高；卫星燃耗无法准确预估，进而对卫星寿命设计有不利影响；控制效率低，且容易造成轨道面相位混乱。相较于绝对构型维持方法，相对构型维持方法控制效率更高；对卫星发射入轨精度要求相对较低；星座全球覆盖性基本固定，利于星座性能分析。

因此本文依据上述讨论，决定采用相对相位维持的策略，通过调整轨道半长轴来维持所有卫星之间的相对位置关系。尽管卫星的绝对经度产生了漂移，但是星座整体构型能够维持稳定。

2.2　姿态与能源约束

2.2.1　卫星姿态约束

卫星可动部件的运动引起的干扰力矩，如太阳帆板驱动机构转动、定向天线运动等，都会影响卫星姿态控制精度。而这些运动部件的挠性振动与卫星姿态运动之间的相互耦合，会直接影响卫星对姿态稳定度的要求。

卫星的太阳帆板面积远大于卫星本体，且由电机驱动以满足对日定向运动。电机产生的力矩不仅驱动帆板旋转，也会反作用于卫星本体，进而对卫星姿态产生扰动。不仅如此，帆板的大幅度转动还会改变卫星质量特性，导致系统惯量出现不确定性，进而影响卫星的姿态控制性能。因此为尽量减小帆板转动对姿态的影响，设计卫星迎风面积控制方法时考虑帆板旋转角速度约束，并且尽可能减小旋转幅度。

图2中，A代表卫星本体，B为太阳帆板。C₀为卫星质心，C₂为帆板质心，C₀是轨道坐标系e1的坐标原点，e1_z轴指向地心，e1_x轴指向卫星速度方向，并在轨道平面内与e1_z轴垂直。

图2  卫星坐标系示意图

Fig. 2  Schematic diagram of satellite coordinate system

坐标系e2为星体固连坐标系，在无姿态误差时，e2与e1完全重合。滚动角、俯仰角和偏航角描述e2相对于e1的转动。

坐标系e3是帆板的固连坐标系，帆板转动时，设α为转角，可得坐标系e2与e3之间转换关系：

$${\mathit{e}}_{\mathrm{3}}={\mathit{A}}_{\mathrm{32}}\left(\alpha \right){\mathit{e}}_{\mathrm{2}}$$

（7）

其中：${\mathit{A}}_{\mathrm{32}}\left(\alpha \right)$为坐标转换矩阵。

在此简化，假定帆板完全刚性，可将卫星动力学方程表示为

$$\left[{\mathit{J}}_{\mathrm{a}}+{\mathit{J}}_{\mathrm{b}\mathrm{0}}+{\mathit{A}}_{\mathrm{23}}\left(\alpha \right){\mathit{J}}_{\mathrm{b}}{\mathit{A}}_{\mathrm{32}}\left(\alpha \right)\right]\ddot{\Phi }\left(t\right)+{\mathit{A}}_{\mathrm{23}}\left(\alpha \right){\mathit{J}}_{\mathrm{b}}\ddot{\alpha }\left(t\right)=\mathit{U}\left(t\right)+{\mathit{{\rm T}}}_{\mathrm{d}\mathrm{a}}\left(\alpha \right)-{\mathit{G}}_{\mathrm{a}}$$

（8）

$${\mathit{D}}^{\mathrm{T}}{\mathit{J}}_{\mathrm{b}}{\mathit{A}}_{\mathrm{32}}\left(\alpha \right)\ddot{\Phi }\left(t\right)+{\mathit{D}}^{\mathrm{T}}{\mathit{J}}_{\mathrm{b}}\mathrm{\Delta }\ddot{\alpha }\left(t\right)={\mathit{D}}^{\mathrm{T}}{\mathit{T}}_{\mathrm{c}}\left(t\right)+{\mathit{D}}^{\mathrm{T}}{\mathit{T}}_{\mathrm{d}\mathrm{b}}\left(\alpha \right)-{\mathit{D}}^{\mathrm{T}}{\mathit{G}}_{\mathrm{b}}$$

（9）

式中：${\mathit{J}}_{\mathrm{a}}$为卫星在星体固连坐标系下惯量矩阵；${\mathit{J}}_{\mathrm{b}}$为帆板在帆板的固连坐标系下的惯量矩阵；${\mathit{J}}_{\mathrm{b}\mathrm{0}}$为$U\left(t\right)$帆板质心C₂在e₂下的惯量矩阵；${\mathit{G}}_{\mathrm{a}}$，${\mathit{G}}_{\mathrm{b}}$为轨道耦合力矩；${\mathit{T}}_{\mathrm{c}}\left(t\right)$为帆板输出力矩；${\mathit{T}}_{\mathrm{d}\mathrm{a}}\left(\alpha \right)$，${\mathit{T}}_{\mathrm{d}\mathrm{b}}\left(\alpha \right)$为外干扰力矩。

由式（9）解出代入（8）可以得到卫星姿态动力学方程：

$$\left[\mathit{J}+\mathrm{\Delta }\mathit{J}\left(\alpha \right)\right]\ddot{\Phi }\left(t\right)=\mathit{U}\left(t\right)-\mathit{H}\left(\alpha \right){\mathit{T}}_{\mathrm{c}}\left(\alpha \right)-\mathit{G}\left(\alpha \right)$$

（10）

式中：$\mathit{H}\left(\alpha \right)={\left({\mathit{D}}^{\mathrm{T}}{\mathit{J}}_{\mathrm{b}}\mathit{D}\right)}^{-\mathrm{1}}{\mathit{A}}_{\mathrm{23}}\left(\alpha \right){\mathit{J}}_b\mathit{D}{\mathit{D}}^{\mathrm{T}}$；$\mathit{J}={\mathit{J}}_{\mathrm{a}}+{\mathit{J}}_{\mathrm{b}\mathrm{0}}$；$\mathit{J}={\mathit{J}}_{\mathrm{a}}+{\mathit{J}}_{\mathrm{b}\mathrm{0}}$；$\mathrm{\Delta }\mathit{J}\left(\alpha \right)={\mathit{A}}_{\mathrm{23}}\left(\alpha \right){\mathit{J}}_{\mathrm{b}}{\mathit{A}}_{\mathrm{32}}\left(\alpha \right)-\mathit{H}\left(\alpha \right){\mathit{J}}_{\mathrm{b}}{\mathit{A}}_{\mathrm{32}}\left(\alpha \right)$；$\mathit{G}={\mathit{G}}_{\mathrm{a}}+\mathit{H}\left(\alpha \right){\mathit{G}}_{\mathrm{b}}$

由姿态动力学方程可知，帆板转角会引起卫星姿态扰动。

此外帆板转动引起的质量特性变化也需要考虑。因此需要对帆板转动运动及其帆板根铰链进行受力分析。在三维空间下，根据经典力学和矢量运算的相关知识可以进行如下运算。

卫星角速度乘以卫星质心到帆板质心矢量，可得出帆板线速度矢量，即为

$$\mathit{V}=\mathit{\omega }\times {\mathit{L}}_{\mathrm{1}}$$

（11）

式中：$\mathit{\omega }$为某时刻卫星的角速度矢量；$\mathit{V}$为太阳翼在质心处速度矢量；L₁为卫星质心到帆板质心矢量。

由线速度与角速度可得帆板质心处沿切线方向加速度，可表示为

$$\mathit{a}=\mathit{\omega }\times \mathit{V}=\mathit{\omega }\times \mathit{\omega }\times {\mathit{L}}_{\mathrm{1}}$$

（12）

之后根据力矩公式：

$$\mathit{T}=\mathit{F}\times {\mathit{L}}_{\mathrm{2}}=m\cdot \mathit{a}\cdot {\mathit{L}}_{\mathrm{2}}$$

（13）

将式（12）代入力矩公式，可得：

$$\mathit{T}=\mathit{F}\times \mathit{L}=m\times \mathit{\omega }\times \mathit{\omega }\times {\mathit{L}}_{\mathrm{1}}\times {\mathit{L}}_{\mathrm{2}}$$

（14）

式中：$\mathit{F}$为太阳翼在质心处的向心力矢量，$\mathit{T}$为帆板根铰链在帆板向心力矢量作用下的力矩矢量，L₂为帆板质心到帆板铰链连接处矢量。

通过上述对卫星姿态性能以及帆板转动时质量特性两方面的结合分析，根据实际情况，一般帆板转速不超过0.5 °∙s^－1 。

2.2.2　卫星能源约束

在卫星在轨运行期间，卫星所需的能量均需依靠太阳帆板提供。由于太阳帆板面积恒定，因此其所能提供的能量大小直接取决于帆板法线和太阳矢量的夹角^［26］。

三轴稳定卫星为保持天线面的对地指向和太阳对帆板的直射，采取了动态偏航控制，图3描述了三轴稳定卫星如何偏航控制，以质心轨道坐标系为基准，利用欧拉角进行定位，即用X轴、Y轴、Z轴来描述定位。其中，X轴为滚动轴，Y轴为俯仰轴，Z轴为偏航轴。当太阳高度角β较大时，太阳帆板虽然可以沿着特定方向旋转以跟踪太阳，且能满足精度要求，但无法满足整星能源供电要求。

图3  三轴稳定卫星偏航机动示意图

Fig. 3  Sketch of three axis stabilized satellite yaw maneuver

因此面向工程实际，为保证卫星有足够的能源供给，在卫星设计时均留有一定的能源余量，因而，在整体的能源需求保证下，帆板在一定转角范围内的转动不会影响卫星平稳运行。实际能源吸收c与帆板法向与太阳矢量之间夹角α之间的关系为

$$c=E·\mathrm{c}\mathrm{o}\mathrm{s}\alpha$$

（15）

式中：E为帆板完全朝向太阳矢量时能源吸收量，在一般工程约束下；α取±5°之间。但是由于帆板需要同时兼顾利用大气阻力进行相位维持的功能，故此设计了一个能源损耗与α角之间的关系（表2）。

表2  能源损耗与α角关系表

Tab 2  Relation between energy loss and α angle

$\alpha$/（°）	能源损耗/％
18.19	5
25.84	10
36.87	20
45.57	30

因此在设计帆板旋转控制率时，需要考虑能源损耗，根据工程经验，将能源损耗30 ％时，$\alpha$角为45.57°作为帆板转角的最大约束。

2.3　轨道保持策略

首先建立半长轴与迎风面积之间的数学关系。用轨道半长轴表示卫星角速度：

$$n=\sqrt[]{\frac{\mu }{a^{\mathrm{3}}}}$$

（16）

卫星角加速度：

$$\dot{n}=-\frac{\mathrm{3}}{\mathrm{2}}n\frac{\dot{a}}{a}$$

（17）

半长轴变化量与迎风面积关系：

$$\dot{a}=-C_d\frac{S}{m}n{a}^{\mathrm{2}}\rho$$

（18）

式中：$C_{\mathrm{d}}$为阻力系数；$\rho$为航天器所在高度大气密度；$S$为卫星的迎风面积；$m$为卫星的质量。

由式（16）与（17）可知，卫星角加速度可表示：

$$\dot{n}=\frac{\mathrm{3}}{\mathrm{2}}{C}_{\mathrm{d}}\frac{S}{m}\frac{\mu \rho }{a^{\mathrm{2}}}$$

（19）

两星在t时间后转速差可表示为

$$\mathrm{\Delta }n={\dot{n}}_{\mathrm{1}}t-{\dot{n}}_{\mathrm{2}}t$$

（20）

因此在经过t时间后，两星相位差可以表示为

$$\mathrm{\Delta }u=\mathrm{\Delta }{u}_{\mathrm{1}}+\frac{\mathrm{1}}{\mathrm{2}}({\dot{n}}_{\mathrm{1}}-{\dot{n}}_{\mathrm{2}})t^{\mathrm{2}}$$

（21）

其中：$\mathrm{\Delta }{u}_{\mathrm{1}}$为两星初始相位差。

将式（20）代入式（21），可以发现决定相位差的关键变量为角速度差，二者关系为

$$\mathrm{\Delta }u=\mathrm{\Delta }{u}_{\mathrm{1}}+\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }nt$$

（22）

因此可以将$\mathrm{\Delta }n$作为控制相位差的主要变量，之后建立迎风面积与$\mathrm{\Delta }n$之间的关系，将式（19）代入式（20）可得到S2星迎风面积与两星转速差关系：

$$S_{\mathrm{2}}=\frac{a_{\mathrm{2}}^{\mathrm{2}}(\mathrm{3}t{C}_{\mathrm{d}}\mu \rho S_{\mathrm{1}}-\mathrm{2}\mathrm{\Delta }nm{a}_{\mathrm{1}}^{\mathrm{2}})}{\mathrm{3}T{C}_{\mathrm{d}}\mu \rho a_{\mathrm{1}}^{\mathrm{2}}}$$

（23）

式中：$S_{\mathrm{z}\mathrm{h}}\le S_{\mathrm{2}}\le S_{\mathrm{z}\mathrm{h}}+S_{\mathrm{F}}$ ；${\mathrm{S}}_{\mathrm{z}\mathrm{h}}$为卫星主体迎风面积；$S_{\mathrm{F}}$为卫星帆板迎风面积。

此外需要定义卫星有效迎风面积S与帆板迎角θ的关系，迎角θ如图4所示。

图4  帆板迎角定义图

Fig. 4  Definition of solar panel incidence

有效迎风面积由卫星主体迎风面积加帆板有效迎风面积组成，帆板厚度忽略不计，表达如下：

$$S=S_{\mathrm{z}\mathrm{h}}+S_{\mathrm{F}}\mathrm{s}\mathrm{i}\mathrm{n}\theta$$

（24）

其中：$\mathrm{0}°\le \theta \le \mathrm{180}°$ 。

2.4　控制策略步骤

在以上数学关系的基础上，帆板控制策略流程如图5所示。以初始时刻处于不同相位、不同迎风横截面积但位于同一轨道面的双星为例。

图5  相位保持策略流程

Fig. 5  Process of phase keeping strategy

（1）获取两星实时轨道数据后，进行高精度轨道外推，外推时长为一个轨道周期。由于高精度轨道外推不是本文研究重点，特此对高精度轨道外推模块进行简要说明。

轨道外推模块用于生成卫星高精度位置速度预报及其轨道状态转移矩阵。其首先分别构建轨道信息预报模型及轨道状态转移矩阵模型。而后，利用高精度轨道摄动模块得到的卫星加速度信息及轨道确定模块输出的初始轨道信息，采用数值积分方法得到所需的高精度卫星轨道外推信息及状态转移矩阵信息。

（2）计算两星下一周期角速度差

两星在T时间后的角速度差为

$$\mathrm{\Delta }n=\left({\dot{n}}_{\mathrm{1}}-{\dot{n}}_{\mathrm{2}}\right)T$$

（25）

（3）判断是否进行相位控制

当两星角速度差为0时，S2星不进行控制。当两星角速度差不为0时，则进行角速度控制。

（4）判断角速度变化

首先计算两星偏移相位差，作为输入量。由式（17）可知迎风面积与角加速度之间为正相关，如果S1星角速度大于S2星角速度，则需增加S2迎风面积。如果S1星角速度小于S2星角速度，则减小S2迎风面积。

（5）确定控制量

本方法希望能尽可能减少帆板的转动幅度。由于有效迎风面积的约束，需要确定进行控制的时间范围。相位控制采用高频次小量控制，根据轨道外推的时长进行控制，即卫星一轨周期。将$T$代入式（23）得到$S_{\mathrm{2}}$，即为下一轨内修正相位差所需的新的迎风面积。

（6）判断$S^{\text{'}}$是否落在约束区间内，如果$S^{\text{'}}\le \left[S_{\mathrm{z}\mathrm{h}}\right]$则卫星在这一轨运行时有效迎风面积等于卫星本体面积，如果$S^{\text{'}}\ge \left[S_{\mathrm{z}\mathrm{h}}+S_F\right]$时，卫星在这一轨运行时有效迎风面积等于卫星本体加最大帆板面积。当$S^{\text{'}}\in \left[S_{\mathrm{z}\mathrm{h}},S_{\mathrm{z}\mathrm{h}}+S_F\right]$时，进行下一步迎角计算。否则$S^{\text{'}}$保持当前大小，直到下一个轨道预报周期重新进行计算。

（7）确定帆板迎角

根据式（24），求出新的帆板迎角。继续获取实时数据，重复上述步骤。

2.5　燃耗分析

化学推力通常是给被控卫星一个瞬时速度增量，以达到维持其轨道的效果。如果采用传统的主动构型维持方法，不仅会因为推力器点火导致LEO星座服务中断，整体性能下降，也会进一步增加燃耗，此部分将分析采用主动维持方法所需的燃耗。

卫星相位变化与半长轴变化对应关系：

$$\mathrm{\Delta }u=-\frac{\mathrm{3}n}{\mathrm{2}a}\mathrm{\Delta }a\mathrm{\Delta }t$$

（26）

卫星在运行一轨后的半长轴变化量为

$$\mathrm{\Delta }a=-\mathrm{2}\mathrm{\pi }{C}_{\mathrm{d}}\frac{S}{m}{a}^{\mathrm{2}}\rho$$

（27）

由式（28）可以推算出卫星一轨的相位变化量：

$$\mathrm{\Delta }u=-\frac{\mathrm{6}{\mathrm{\pi }}^{\mathrm{2}}{C}_{\mathrm{d}}S\rho a}{m}$$

（28）

一个轨道周期后两星相对相位差可表示为

$$\mathrm{\Delta }u=\mathrm{\Delta }{u}_{\mathrm{1}}-\mathrm{\Delta }{u}_{\mathrm{2}}$$

（29）

进而计算出最大允许漂移量$\psi$与$\mathrm{\Delta }u$之间的倍数关系，即为卫星绕地圈数。将圈数与轨道周期相乘即可获得相位漂移临界时长$t$。

由式（26）可知：

$$\mathrm{\Delta }a=-\frac{\mathrm{2}a\mathrm{\Delta }u}{\mathrm{3}n\mathrm{\Delta }t}$$

（30）

式中：$\mathrm{\Delta }u=\psi$。由此可得卫星经过$\mathrm{\Delta }t$时间后的半长轴变化量$\mathrm{\Delta }a$。由于需要采用主动构型维持来修正半长轴变化量$\mathrm{\Delta }a$，推力器点火产生的速度增量$\mathrm{\Delta }V$可表示为

$$\mathrm{\Delta }V=\frac{n}{\mathrm{2}}\mathrm{\Delta }a$$

（31）

燃耗$\mathrm{\Delta }m$与速度增量$\mathrm{\Delta }V$关系为

$$\mathrm{\Delta }m=m_{\mathrm{0}}\cdot {\left(\mathrm{1}-e\right)}^{-\frac{\mathrm{\Delta }v}{I_{\mathrm{s}\mathrm{p}}\cdot \mathrm{9.8}}}$$

（32）

式中：m₀为卫星初始质量；I_sp为卫星推力器比冲。

由此可得卫星采用主动维持时所需的燃耗，而这部分燃耗可以通过旋转帆板，调整迎风面积来节省。

3 仿真分析

对同轨道面内两颗卫星S1星，S2星的相位保持情况进行仿真，周期为一年，步长60 s，使用化学推力器。两星初始轨道参数与仿真环境见表3。

表3  LEO仿真环境

Tab. 3  LEO simulation environment

仿真参数	数值
航天器质量/kg	20
轨道高度/km	1 000
轨道倾角/（°）	45
空气阻力系数/常数	2.2
卫星主体迎风面积/㎡	0.5
帆板最大迎风面积/㎡	10
最大相位差容许偏差/（°）	±2.0
S1星初始真近点角/（°）	0.0
S2星初始真近点角/（°）	30.0
推力器推力/N	5

该仿真过程使用高精度轨道外推模型，主要考虑了非球形引力摄动、日月第三体引力摄动以及大气阻力摄动，仿真结果如图6—8所示。为期一年的仿真表明，通过调整帆板迎风面积进行相位保持能够实现稳定的相位保持，一年中两星相位差由初始额定值30.0°经过260 d到达最大值，但小于29.9°，在仿真周期一年时间内，两星相位差始终不超过0.1°，远小于本文设定的±2.0°的最大容许偏差，并在后期逐渐收敛并靠近初始相位差。此外两星的半长轴差最大值约为31 m，且始终不超过40 m，倾角演化趋势基本相同，始终未超过$\pm \mathrm{6}\times {\mathrm{10}}^{-\mathrm{5}}$ °，因此该方法能够实现LEO星座的相对构型控制。

图6  两星相位偏差变化曲线

Fig. 6  Phase deviation variation of two satellites

图7  两星半长轴偏差变化曲线

Fig. 7  Variation of semi-major axis deviation of two satellites

图8  两星倾角差变化曲线

Fig. 8  Variation of inclination difference between two satellites

由图9可见，该方法使帆板小幅度旋转，根据实际控制需求在每一个周期旋转一定角度，从而产生新的迎风面积，而不是简单的调整帆板至最大或者最小迎风面积。

图9  S2星迎风面积变化曲线

Fig.9  Variation of solar panel incidence of S2 satellites

在卫星相对构型维持情况下，通过初步燃耗分析可知，当最大相位容许偏差为$\pm$2.0°时，在无控情况下两星相对相位差在整个寿命期间会有5次漂出允许范围，单次主动相位控制所需的燃耗为2.66 kg，卫星寿命期间总共消耗燃料13.28 kg。考虑到卫星需修正入轨偏差以及更复杂情况下进行主动相位保持，真实燃耗会大于初步分析结果。

4 结 论

本文提出了一种利用大气阻力进行LEO星座构型自主维持的方法，在充分考虑帆板转动对卫星质量分布的影响，能源影响以及姿态指向的基础上，在同一轨道面上有效实施了一年内卫星相位的稳定控制，使得两颗卫星的相位偏差量稳定在0.1°以内，完全符合预期±2.0 °的目标。两颗卫星的半长轴差可以控制在40 m以内。两星倾角差也不超过$\pm \mathrm{6}\times {\mathrm{10}}^{-\mathrm{5}}°$，可见该方法能够稳定维持星座构型。相较于主动进行相位保持的相对构型维持方法，利用帆板迎风面进行相位保持可以在卫星寿命周期内减少至少13.284 kg燃耗，且一年之内无需进行主动保持。可见该方法能够有效减少主动控制频率，减少燃耗，增强星座服务连续性，这对于低轨星座在通信宽带以及增强导航方面的应用都有较大意义。

作者贡献声明

郑 奕：研究构思，参与实验，数据处理，撰写论文。

林宝军：研究构思，学术指导，论文修改，项目支持。

刘迎春：研究构思，学术指导，论文修改。

林 夏：学术指导，论文修改，项目支持。