全方向长期伴飞卫星集群轨道设计与仿真*

韩耀昆，陈琪锋

(中南大学 航空航天学院·长沙·410083)

0 引 言

集群飞行是近年来空间技术的一个新兴概念，即用多个小微卫星代替单一卫星实现系统的功能[1]，每个卫星在空间上相互独立，通过无线组网进行交互。与传统单一卫星相比，集群的优势在于一旦某个卫星发生故障，对系统功能的影响较小，可以增强空间生存能力和快速反应能力，而单一卫星一旦故障，会丧失部分功能甚至完全无法工作[2-5]。除此之外，随着微小卫星和皮纳卫星等技术的不断发展成熟，集群任务已经成为可能，相较于普通卫星，可以大大降低成本[6-7]。但是集群也有一定的劣势，卫星集群实质上是将单一卫星实现的功能变为多颗小微卫星协同完成[6-10]，所以会增大协同管理和协同控制的难度，并且随着集群数量的增多，复杂性也会大大增加[10-12]。

集群任务不同于精确编队，不要求长期保持严格的构型和精确的相对位置保持，仅需要保持航天器松散地长期伴飞，不需要巨大的燃料消耗，主要依靠自身摄动保持构型[2-11]。而在航天动力学中，航天器要受到J2摄动、大气阻力、太阳光压等多种摄动因素的作用[12-15]，由于集群航天器中各航天器的初始轨道根数和面质比等参数都可能有一定的差异，因此每一个航天器之间的漂移(面内和面外)是不可避免的[2,12]，也就导致航天器的整体构型会发生变化。为使集群构型保持长期稳定，需要建立相对运动动力学模型，计算长期稳定的相对运动轨迹，Sedwick等[16]将J2摄动以力函数的形式引入到C-W方程中，Schweighart等[17]则给出了包含J2摄动的线性微分方程。参考轨道为椭圆轨道的情况则是被Lawden等[18]和 Tshauner等[19]解决，他们两人分别独自于1963年和1965年写出了线性化的相对运动方程并求出其解析解，因此该方程叫T-H方程。Carter等[20]则解决了T-H方程初始值奇异的问题。此外，也可以用数值积分法进行相对运动的研究，此方法最为直接，且由于计算过程不需要进行简化或线性化，故而精度较高，只是无法给出解析解，常用来进行验证。Schaub[6]等人通过对航天器平均轨道根数长期漂移速率的分析研究，提出以平均轨道要素的差异来抵消J2摄动的影响，并得到了J2摄动条件下相对运动轨道的初始条件。Gurfil[7]等人基于能量匹配建立了开普勒轨道下椭圆运动动力学模型，进而得出了相对有界周期运动的条件，并利用单冲量控制方法来维持相对有界。党朝辉[2]求解了LVLH坐标系下各坐标分量的距离上下界及星间距离的上下边界，并据此提出包络球和包络盒的模型。

本文提出了一种在J2摄动的影响下，在中心卫星周围任意给定初始方向上部署长期稳定伴飞卫星的轨道设计方法。通过在中心卫星周围不同方向上部署，实现卫星集群对中心卫星的全方向伴飞，并进行大规模集群的仿真验证。本文星群构型设计的目标是使星群分布于中心卫星的不同方向上形成全向包围分布并能够长期稳定运行的构型。主要技术难点在于求解在任意给定时刻处于中心卫星任意给定方向上，并且能够稳定伴随飞行的伴飞卫星轨道根数。

1 全方向长期伴飞卫星轨道设计方法

1.1 坐标系

如图1所示，分别定义两个坐标系。

地球为中心的惯性坐标系，简称地心惯性坐标系(Earth Centered Inertial, ECI)，用来描述卫星轨道运动。ECI坐标系原点在地心Oe，基本平面为J2000地球平赤道面，X轴由地心指向J2000春分点，Z轴沿地球自转轴从地心指向北极，Y轴与X，Z轴垂直且由右手规则确定。向量r=[X,Y,Z]T表示在地心惯性坐标系中的空间位置。

当地水平当地垂直坐标系(Local-Vertical-Local-Horizontal，LVLH)用来描述轨道上两颗近距离飞行的航天器之间的相对运动，为方便，一般将其中一个航天器称为中心卫星(Center Satellite)，另一个航天器称为伴飞卫星(Deputy Satellite)[6]。在下文中，下标(·)c和(·)d分别代表中心卫星C和伴飞卫星D。LVLH坐标系原点为中心卫星质心，x轴由地心指向中心卫星质心，z轴垂直于轨道平面且指向轨道角动量方向，y轴与x、z轴垂直且由右手规则确定。向量ρ=[x,y,z]T表示在LVLH坐标系中的空间位置。

轨道根数以向量形式描述为：

L=[a,e,i,ω,Ω,M]

(1)

其中a表示半长轴，e表示偏心率，i表示轨道倾角，ω表示近地点角距，Ω表示升交点赤经，M表示平近点角。

1.2 卫星长期伴飞的J2不变条件

将J2摄动加到经典开普勒轨道运动中会导致密切轨道根数的三种变化，即短期和长期振荡以及长期增长。短期项表现为轨道根数的振荡，但不会导致轨道漂移[6]。本文研究内容主要避免的是长期增长，这种增长用平均轨道根数来描述比用密切轨道根数更加简单[8]，本文以下公式如无说明均为平均轨道根数。

对于开普勒运动(即不存在摄动)，只有平近点角M是与时间相关的轨道元素。M增长的速率由平均轨道速率n给出，平均轨道速率也叫平均角速度，定义为n=2π/T，其中T为轨道周期。对于有界的两颗卫星之间的相对运动，平近点角变化速率必须相等，即差值为零：

(2)

式中因为平近点角的变化率取决于a，所以公式(2)可近似为一阶：

(3)

(4)

根据平均摄动法[8]，J2摄动下的升交点赤经、近地点角距和平近点角变化率为

(5)

(6)

(7)

为了设计长期稳定伴飞构型，采用J2不变的条件[6]：

(8)

其中θM为平均纬度幅角，公式(8)即表示升交点赤经和纬度幅角的漂移率相同，从而两颗卫星相对位置不会漂移。由公式(7)可知二者都是L，η，i的函数，中心卫星与伴飞卫星升交点赤经和平均纬度幅角变化率的差值为：

(9)

将公式(7)代入公式(9)并线性化和略去高阶项之后可以将上式简化为[6]：

(10)

可以看出，公式(10)是a、e、i三个轨道根数之间的关系，因而平均轨道根数Ω、ω和M不会直接导致J2引起的长期增长，所以这三个参数可以任意选取。

满足公式(10)条件的伴飞卫星在J2摄动的影响下可以实现对中心卫星长期伴飞，对于中心卫星为圆轨道和椭圆轨道均可用此方法实现。

1.3 全方向部署的轨道设计方法

要实现伴飞卫星对中心卫星的全方向部署，需要在任意方向都可以部署伴飞卫星，本节给出在任意方向部署的轨道设计方法，基本设计思路为利用轨道根数的差异在不同方向产生位移。公式(10)中的两个公式是半长轴、偏心率和轨道倾角三个轨道根数之间的关系，可以通过改变偏心率在LVLH坐标系中的x方向产生的位移，y方向的位移利用改变近地点角距和平近点角实现，z方向的位移则通过改变升交点赤经实现[16-17]。

在地心惯性坐标系中的中心卫星平均轨道根数为Lc=[ac,ec,ic,ωc,Ωc,Mc]，伴飞卫星的平均轨道根数为Ld=[ad,ed,id,ωd,Ωd,Md]。伴飞卫星与中心卫星轨道根数之差为δ=[δa,δe,δi,δΩ,δω,δM]。在LVLH坐标系中，定义伴飞卫星相对于中心卫星的位置为ρ0=[x0,y0,z0]T，

(11)

其中ε为ρ在xy平面的投影与x轴夹角，ψ为ρ与z轴的夹角，下标0表示初始时刻，如图2所示。

图2 LVLH坐标系中ε和ψ的定义Fig.2 Define ε and ψ in LVLH coordinate system

根据开普勒方程，得到中心卫星的Mc和ec以及偏近点角Ec之间的关系[8]：

Ec=Mc+ecsinEc

(12)

利用牛顿迭代法，即可求出Ec。

定义旋转矩阵Pc,Qc分别为[8]：

(13)

(14)

式中uc为中心卫星的纬度幅角，引入真近点角fc，则uc=fc+ωc，真近点角fc可由下式求得[8]，符号与平近点角一致：

(15)

从而得到中心卫星在地心惯性坐标系中的空间位置rc：

(16)

同样，伴飞卫星在地心惯性坐标系中的空间位置rd：

(17)

则在地心惯性坐标系中伴飞卫星相对于中心卫星的空间位置为r0=rd-rc。

r0与ρ0之间的关系为

r0=Tρ0

(18)

其中T为从LVLH坐标系到地心坐标系的旋转变换矩阵，为

(19)

根据rc和r0即可求得rd=r0+rc，则有

(20)

根据公式(8)求得δa,δe,δi,之后代入公式(19)，若伴飞卫星相对中心卫星的距离已知，即可求得指定方向伴飞卫星轨道根数，若未知，可得关于距离的表达式。

需要特别指出的是，中心卫星为圆轨道，或者中心卫星的偏心率ec<-δe时，直接计算得到ed<0，在航天轨道力学中无意义，因此δe要分情况讨论。

(1)δe为正

即伴飞卫星初始位置在中心卫星上方(相对于地球)，此时可以直接计算伴飞卫星的轨道根数。

(2)δe为负

即伴飞卫星初始位置在中心卫星下方，此时直接计算会导致伴飞卫星的偏心率为负值，不满足偏心率的定义，计算时也会出现错误，为了使伴飞卫星可以在参考卫星下方，需对伴飞卫星轨道根数进行处理。

上述变换处理如图3所示，点代表卫星，实线为中心卫星轨道，虚线为目标伴飞卫星轨道，未经处理之前的轨道为黑虚线轨道，要将黑虚线轨道变为红虚线轨道，则需要将轨道旋转180°，即近地点角距在原来的基础上加180°，而卫星从原来的远地点变为近地点，即平近点角在原来的基础上减去180°，这样即可实现伴飞卫星在中心卫星下方。为避免运算错误，此时计算时δe需要取绝对值。

图3 初始位置在中心卫星下方轨道变换示意图Fig.3 Orbital transformation of the initial position below the center satellite

2 仿真试验及结果分析

2.1 长期伴飞仿真验证

给定中心卫星C的初始平均轨道根数Lc(t0)=[7500 km,0,35°,0°,0°,0°]，设在初始时刻，伴飞卫星的初始方位角为ε0=50°，高度角为ψ0=60°，且|ρ0|=65km，设计伴飞卫星的轨道根数。采用公式(8)计算得到δe=0.004,δa=-0.713 m,δi=0.0025°，再用公式(18)计算得到δΩ=-0.0058,δω=0.0091，δM=0.000013，最后结合轨道根数的定义，得到地心惯性坐标系伴飞卫星初始轨道根数为Ld(t0)=[7499999.287 m,0.004,35.0025°,359.6677°,0.5214°,0.000745°]。

之后将平根数转换为密切轨道根数，进行数值仿真。其中地球非球形引力场摄动仅考虑J2摄动，轨道预报模型选取高精度的HPOP模块，重力模型采用JGM3，其中n=21，大气阻力模型采用1976-Standard模型。为验证长期伴飞的稳定性，仿真时间为60天。两航天器的参数为：质量mc=md=15kg，阻力系数CDc=CDd=2.2，最大截面积Sc=Sd=0.09 m2，其他参数皆采用默认参数，且所有卫星一致。60天内得相对距离变化情况为图4。可以看出距离始终在小范围内波动，可以实现长期稳定伴飞。

图4 2颗卫星60天内距离变化Fig.4 Distance change of 2 satellite within 60 days

2.2 卫星集群仿真验证

本节对大规模(100颗)卫星集群全方向伴飞构型的设计进行仿真验证。中心卫星轨道根数及其他参数同上例，随机生成100颗伴飞卫星相对中心卫星的位置向量，方位角随机选取保证初始全方位覆盖，初始距离在100km以内随机选取。

60天内相对距离变化如图5。在LVLH坐标系下表示的60天内不同阶段构型如图6，中间标注点代表中心卫星，周围其他点表示伴飞卫星。

图5 100颗卫星60天内距离变化Fig.5 Distance change of 100 satellites within 60 days

(a) 初始构型

(b) 20天构型

(c) 40天构型

(d) 60天构型图6 100颗卫星60天内构型变化Fig.6 Formation configuration of 100 satellites within 60 days

由图5和图6可以看出，60天内集群一直都可以保持环绕中心卫星的松散伴飞构型，距离尺度没有显著变化。

考虑到伴飞卫星之间可能发生碰撞，统计了所有伴飞卫星两两之间在运行过程中的最小距离，考虑到本文所研究的大规模卫星集群主要为微纳卫星或皮卫星，体积尺度在厘米级或分米级，设置卫星之间安全距离为100 m，统计结果如图7，为方便观看及统计，作二维图，图中横坐标表示伴飞卫星的个数，纵坐标表示伴飞卫星两两之间在运行过程中的最小距离，对于碰撞问题需要特别关注的是距离下界，故对纵坐标做对数化处理，所有横坐标为i(i=1,2，…，100)的点表示第i颗卫星与其它伴飞卫星在运行过程中的最小距离，为避免重复统计，对于第i颗卫星，只统计其与序号大于i的伴飞卫星之间运行过程中的最小距离。由图7可以看出伴飞卫星运行过程中共有10次在安全距离以下，从总运行时长看，运行总圈数约为800圈，碰撞概率为0.0125次/圈，碰撞概率很低。

图7 伴飞卫星60天内两两之间最小距离Fig.7 The minimum distance between deputy satellites within 60 days

本文的初步设计没有考虑星群卫星构型设计中碰撞避免的问题，因为所考虑的星群构型中卫星间的碰撞概率很低，这一方面是卫星间距离尺度较大，在100 km量级；另一方面是因为星群各卫星较均匀地部署在参考星的周期各方向上，并且利用J2不变轨道的特点，星群运行中各环绕星的相对运动关系比较稳定。在低碰撞概率情况下，在运行中根据状态进行碰撞预测，并提前采取规避可以作为一种处理方式。关于星群无碰撞构型设计问题，还需要依据轨道动力学规律进一步深入研究。

3 结 论

本文提出了在指定初始方向部署卫星并且可以实现长期稳定伴飞的一种伴飞卫星轨道设计方法，并将其应用于全向伴飞卫星集群的构型设计。通过对平均轨道根数的分析，利用J2不变的概念，保证了卫星在J2摄动的条件下相对漂移为零或很小；从几何关系上对轨道进行分析，得到相对位置与轨道根数之间的关系，进而找出在指定方向部署卫星的方法；用上述方法开展全向伴飞卫星集群构型设计，通过仿真验证了其有效性。