基于对偶四元数的空间引力波探测器稳定构型优化设计

张锦绣，张 谕，王继河，杨继坤，逯振坤，宋昱岐

（1.中山大学 航空航天学院，深圳 518107；2.中山大学 物理与天文学院，珠海 519082）

引 言

为印证广义相对论，引力波探测成为当前国际物理和空间科学领域的研究热点。国内外各研究机构围绕空间引力波探测技术展开了深入研究，如：欧洲航天局（European Space Agency，ESA）的“激光干涉测量空间天线”（Laser Interferometer Space Antenna，LISA）/eLISA（evolved LISA）计划[1-2]，中国的“天琴计划”[3]、“太极计划”[4]。

上述引力波探测计划大同小异，如图1所示，它们的空间引力波探测器均由3个航天器组成等边三角形构型，C位于构型平面的中心，表示地球或者虚拟的点，航天器围绕太阳周期性运动的同时也围绕C周期性运动。另外，它们都假设空间引力波会造成航天器搭载的检验质量相对位置发生改变，并通过测量该变化量来印证引力波的存在。检验质量之间的相对位置，是通过安装在航天器上的激光干涉测距系统测量得到的。为了保证高精度激光干涉测距系统可以长时间处于稳定对准测量状态，航天器编队需要运行在可以保持长期稳定构型的轨道上。不同的是，LISA/eLISA计划和“太极计划”，构型的中心C为虚拟点；“天琴计划”是以地球为中心天体，构型中心C为地球。

图1 空间引力波探测器构型示意图Fig.1 Schematic diagram of space gravitational wave detector configuration

空间引力波探测编队构型的稳定性主要由3个指标表征：臂长（任意两个探测器之间的距离）、呼吸角（每两个臂之间的角度）和星间视线相对速度。臂长决定了干涉仪的探测频带，不等臂长将会在激光干涉测量中引入激光频率噪声，影响核心科学载荷的设计；呼吸角的变化影响探测器的指向精度；星间视线相对速度影响超稳时钟和相位计等重要载荷的指标设计。探测器为了成功探测到引力波，需要将3个航天器之间的臂长、呼吸角以及相对速度等三角形轨道构型参数长时间保持在一定范围内。准确的航天器动力学模型和轨道优化算法是设计满足上述轨道构型的两大关键要素。

对于如LISA/eLISA计划和“太极计划”这种轨道构型中心为虚拟点的航天器编队方案，多数学者以笛卡尔坐标或轨道要素为状态量，在日心坐标系下推导航天器间相对动力学方程，然后将其简化为线性的C-W方程进行轨道设计。但是，当探测器轨道构型中心有中心天体时（例如“天琴计划”），上述方式不再适用。这是因为航天器之间的臂长和航天器到中心天体之间的距离在同一量级，不满足C-W方程的简化条件。针对有中心天体的航天器编队情况，本文以对偶四元数为工具，在J2000日心坐标系下建立航天器之间的相对运动动力学模型。值得注意的是，本文所提方法同样适用于轨道构型中心为虚拟点的航天器编队方案。

为了得到满足空间引力波探测任务要求的高精度轨道构型，通常需要采用优化算法对航天器运行轨道进行优化。文献[5]使用全局优化的禁忌搜索算法对LISA引力波探测轨道进行了优化设计，价值函数有18个自变量，找到了满足引力波探测要求的轨道，并给出了从地球停泊轨道进入引力波探测实验轨道的发射段和分离段的轨道设计。文献[6]分析了影响LISA编队构型的性能指标（呼吸角、滞后角、相对速度）相对于初始轨道要素的敏感性，将LISA轨道优化问题转化为两步级联单目标优化问题，先优化滞后角，后优化呼吸角，并分析了轨道周期偏差对构型稳定性的影响。针对“天琴计划”轨道构型优化设计，文献[7]对构型发生长周期漂移的轨道要素进行了分析，采用粒子群算法对3颗航天器的轨道半长轴进行优化，初步满足轨道构型的稳定性指标。文献[8]采用组合优化方法，讨论了由于轨道进动而导致探测器指向的缓慢长期漂移，设计了满足构型稳定性要求的任务轨道。

虽然上述研究满足了引力波探测构型稳定性的需求，但是它们均是在已知初始轨道的基础上，对航天器的轨道根数进行优化，并且轨道优化的起始时刻是确定的。没有考虑初始轨道以外是否存在更加稳定的轨道，使得航天器满足引力波探测构型稳定性的需求。运行在甚高轨道的航天器会受到地球、月球和太阳等主要天体的引力影响，若在轨道优化过程中通过星历得到每一时刻地球、月球和太阳的位置，这严重影响了优化算法的计算效率。因此，有必要以太阳为中心建立惯性坐标系，为地球和月球建立简化但足够精确的动力学模型，从而加快优化速度。把地球、月球和航天器看作刚体来描述它们的运动规律，这实际上是描述六自由的运动问题。而对偶四元数[9-13]是目前描述刚体运动的最合适的数学工具，这是因为它将平移和旋转统一在同一框架下，同时考虑了二者之间的耦合影响，并且可以提供很高的计算效率。

综上所述，本文在J2000日心坐标系下，研究日–地–月引力场中，空间引力波探测系统稳定构型的全局搜索优化方案。①基于对偶四元数，在J2000日心坐标系下建立地球、月球和3个航天器的姿轨一体化动力学模型；②分析构型稳定性的表征，确定优化的目标函数和待优化的变量；③使用遗传算法对航天器轨道构型进行优化设计，使得地球、月球和3颗航天器在不进行轨道控制的前提下尽可能长时间地保持稳定的构型；④通过数学仿真验证优化方法的有效性和普适性。

1 日心系编队动力学表征模型

1.1 自然/人造天体动力学方程

考虑日地月引力场对引力波探测器的影响，本文在日心系下，将自然天体（地球、月球）和人造天体（3个航天器）视为一个完整的编队。这样可以在建立编队动力学模型的时候，更加方便地对日地月引力进行表述。

地球在J2000日心坐标系[14]OsXsYsZs下的运动学方程为

地球在J2000日心坐标系OsXsYsZs下的动力学方程为

其中：=fe+ετe为作用于地球质心的对偶力旋量；=medI/dε+εJe表示地球的对偶惯量矩阵；me表示地球的质量；Je表示地球的转动惯量。

月球在J2000日心坐标系OsXsYsZs下的运动学方程为

月球在J2000日心坐标系OsXsYsZs下的动力学方程为

其中：=fm+ετm为作用于月球质心的对偶力旋量；=mmdI/dε+εJm表示月球的对偶惯量矩阵；md表示月球的质量；Jd表示月球的转动惯量。

以航天器的质心为坐标原点，建立本体坐标系Oi−xiyizi，3个坐标轴方向分别与航天器的惯量主轴重合。航天器在日心坐标系下运动学方程表示为[15]

其中：i=1,2,3分别表示航天器SC1、SC2和SC3；对偶四元数qi表示第i个航天器的本体坐标系相对于日心坐标系的四元数，rii表示航天器的位置矢量在本体坐标系下的表示；对偶速度旋量=wii+εvii，wii表示航天器相对于太阳的角速度在本体坐标系下的表示，表示航天器相对于太阳的速度矢量在本体坐标系下的表示。

动力学方程为

其中：=midI/dε+εJi表示对偶惯量矩阵，mi表示航天器的质量，Ji表示航天器的转动惯量；表示航天器受到的对偶力。当图1中C点为虚拟点时，航天器以太阳为中心天体，若只考虑太阳引力，那么航天器所受到的对偶力为fsi表示太阳的引力。当图1中C点以地球为中心天体时，航天器受到地球、月球和太阳的引力和太阳光压，那么航天器所受到的对偶力为：分别表示地球引力、月球引力、太阳引力、太阳光压和重力梯度力矩。

其中：地心引力常数为µe=398 600.441 90 km3/s2，月心引力常数为µm=4 902.800 076 km3/s2；日心引力常数为µs=132 712 440 040.944 00 km3/s2；太阳光压为P⊙表示太阳光压强，它与航天器到太阳的距离有关，A表示航天器表面积，r⊙表示太阳相对于地心的位置矢量，AU表示地球与太阳的距离，δ 表示卫星表面材料的反射系数。

1.2 编队成员间相对动力学方程

本文主要研究航天器之间的相对运动，即两两航天器本体坐标系Oi−xiyizi和Oj−xjyjzj（i,j=1,2,3,i≠j）之间的相对运动。利用乘性误差对偶四元数，该相对运动可表示为

坐标系Oi−xiyizi相对于坐标系Oj−xjyjzj的运动学方程为

航天器之间相对动力学方程为

本研究关注的核心指标是探测器之间的相对运动参数，但由于编队中包含了自然天体，所以需要建立探测器和临近自然天体之间的相对运动关系。

航天器相对于临近天体的相对运动方程为

航天器相对于临近天体的相对动力学方程为

2 稳定构型优化设计

2.1 待优化参数

本文将进入预定轨道的起始时刻d作为优化参数之一，寻找使得航天器构型最稳定的时刻。

在日心坐标系下，确定了轨道优化的起始时刻后，此时地球、月球的初始位置和速度均可通过星历DE421获得，通过式（1）～（4）可得到初始时刻之后地球和月球在日心坐标系下的位置和速度。

航天器本体坐标系相对于日心坐标系的四元数表示为

其中：n为单位轴；σ 为绕单位轴转过的角度；β为方位角；φ为极角。单位轴n=[nx,ny,nz]由极坐标形式表示为

其中：β,φ ∈[0,2π]。

在航天器轨道坐标系中，位置矢量为rii=[r,0,0]T，其中r=p/(1+ecosf)，p=a(1−e2)，a为轨道半长轴，e为偏心率，f为真近点角。速度矢量为vii=[vr,vt,0]T，为径向速度,(1+ecosf)为横向速度。所以当航天器本体坐标系和质心轨道坐标系重合，每个航天器需要6个参数p、σ、β、φ、e、f作为待优化参数，优化得到航天器的轨道。

2.2 目标函数的确定

根据表征等边三角形稳定性要素指标要求，目标函数CF由卫星之间的臂长变化量CF1和呼吸角变化量CF2两部分线性组成，并求最小值。

其中：w1、w2为权重；c1、c2为归一化常数；θ0=60◦；L为卫星之间的臂长。

干涉臂长和呼吸角的计算公式为

其中：i,j,k=1,2,3，i≠j≠k，rij(t)、rik(t)、rjk(t) 为两航天器之间相对位置向量。

2.3 遗传算法优化方法

根据前面的推导结果可以看出，对于空间引力波探测系统，只要确定了构型的优化参数，就可以计算出地球、月球和航天器的位置和速度。取不同的优化参数，就可以得到不同的构型，构型稳定性问题可以抽象为求取目标函数CF(t)全局最小值的问题。

对于引力波探测器轨道的优化，目标函数有19个自变量，存在数量众多的局部极小值。使用简单的局部优化算法易陷入局部极小值点而无法找到满足要求的轨道，所以这里选择遗传算法作为全局优化算法。遗传算法只利用目标函数的信息，无需梯度等高价信息，适用于大规模、高度非线性的不连续的目标函数的优化，具有很强的通用性，是构型优化设计算法的一个理想选择。

2.3.1 编码和适应度函数

遗传算法的编码方式包括二进制编码、浮点数编码等，本文采用二进制编码方式，通过控制二进制的位数决定参数的精度。

选择合适的适应度函数是设计遗传算法过程中的一个关键问题。本文以式（17）作为适应度函数

2.3.2 遗传算子

1）选择算子

通过选择算子选出种群中较优的个体形成一个新的种群，可使得种群中的个体不断逼近最优解。选择算子可采用常用的轮盘算法，设定个体被选中的概率如下

其中：f′为个体j的适应度值；n为种群数目。在上述选择方式下，个体的适应度值越大，被选择的概率也越大，则其基因构造被遗传到下一代的可能性越大。这种选择方式使得种群中的个体都存在被选中的概率，保证了较小个体的遗传性。

2）交叉算子

交叉是指按交叉概率从种群中选择部分个体，通过交换个体部分基因形成新的个体，从而形成新的种群。交叉算子可选择单点交叉算子。

3）变异算子

变异算子模拟的是生物进化过程中的基因突变现象，具体操作为：令变异概率为Pc，则变异的节点个数为[nSPc]，S为单个个体的节点数；根据生成的随机数r决定变异方向，选中某些个体上的某些节点，然后用新值替代这些节点上的原值。变异算子维持了群体的多样性，有利于防止早熟现象的出现。

2.3.3 算法流程

本文算法的计算步骤如下：

1）生成初始种群；

2）根据设计的适应度函数，计算该代种群所有个体的适应度；

3）判断是否满足终止条件，若不满足则到步骤4，若满足则到步骤5；

4）根据交叉概率和变异概率，经过选择算子、交叉算子和变异算子产生下一代种群，然后返回步骤3；

5）选出最后一代种群中的最优个体，输出相应的最优解，算法终止。

3 仿真算例与结果分析

为了验证所提方法的普适性，以地心轨道和日心轨道空间引力波探测任务为应用背景，对构型优化问题进行仿真验证。

3.1 地心轨道空间引力波探测系统构型优化仿真

以地心轨道空间引力波探测系统为应用背景，对所提方法进行仿真验证。3个航天器各自围绕地球飞行，轨道高度为10万km，假设轨道偏心率和真近点角均为0，航天器轨道由4个参数p、σ、β、φ唯一确定，参数范围如表1所示。轨道起始日期为2034年1月1日，在1年的时间内优化得到最优的构型。遗传算法的种群规模为100，交叉概率为0.7，变异概率为0.08，迭代次数为100代。其余参数分别为c1=3/104，c2=3/400，w1=0.5，w2=0.5。P⊙=4.56×10−6Nm2，δ=0.21，A=12m2，AU=149 597 870 691km。将优化结果表示为位置速度后如表2所示。

表1 地心轨道编队待优化参数及取值范围Table 1 Geocentric orbit formation parameters to be optimized and their value ranges

表2 地心轨道编队优化结果Table 2 Geocentric orbital formations optimization results

图2是地心轨道编队航天器SC1、SC2和SC3之间的干涉臂长变化率曲线，图2中绿色曲线表示干涉臂长L12的变化率，最大值∆L12为0.14%；红色曲线表示干涉臂长L13的变化率，最大值∆L13为0.14%；蓝色曲线表示干涉臂长L23的变化率，最大值∆L23为0.18%。

图2 地心轨道编队航天器之间的干涉臂长变化率Fig.2 Geocentric orbit formation arm length variation between spacecraft

图3是地心轨道编队航天器SC1、SC2和SC3轨道构型呼吸角的变化曲线，图中绿色曲线表示呼吸角θ1，变化量最大值 ∆θ1为0.14°；红色曲线表示呼吸角θ2，变化量最大值∆ θ2为0.084°；蓝色曲线表示呼吸角θ3，变化量最大值∆ θ3为0.13°。

图3 地心轨道编队呼吸角Fig.3 Geocentric orbit formation breathing angle

图4是地心轨道编队航天器SC1、SC2和SC3相对速度变化曲线，图中绿色曲线表示SC1和SC2之间的相对速度，变化量最大值∆v1为6.17 m/s；红色曲线表示SC1和SC3之间的相对速度，变化量最大值∆v2为4.56 m/s；蓝色曲线表示SC2和SC3之间的相对速度，变化量最大值 ∆v3为6.70 m/s。

图4 地心轨道编队航天器之间的相对速度Fig.4 Geocentric orbit formation relative velocity between spacecraft

上述结果表明，在考虑了日、地、月引力和太阳光压摄动的情况下，地心空间引力波探测系统的轨道在一年内的变化能满足引力波探测的需求。随着时间的推移，构型逐渐发散，需要通过轨道机动将航天器编队配置恢复到初始状态。

以优化得到的结果作为初始条件，用Monte-Carlo法分析入轨误差对地心轨道编队构型的影响。分别向航天器的径向、切向和法向位置增加标准差为5 m的位置误差，对构型的影响如表3所示。分别向航天器的径向、切向和法向速度增加标准差为0.1 mm/s的速度误差，对构型的影响如表4所示。位置误差和速度误差均符合正态分布。可以看出，径向位置误差和切向速度误差导致构型标准差的变化最大。因此，航天器的入轨误差中径向位置误差和切向速度误差对构型稳定性的影响最大。若位置误差和速度误差同时存在，当速度误差是1 mm/s时，位置误差超过5 m 时构型将发散。

表3 位置误差对地心轨道编队构型的影响Table 3 Effect of position error on configuration of geocentric orbital formation

表4 速度误差对地心轨道编队构型的影响Table 4 Effect of velocity error on configuration of geocentric orbital formation

3.2 日心轨道空间引力波探测系统构型优化仿真

以日心轨道空间引力波探测系统为应用背景，对所提方法进行验证。该方案中3颗航天器各自围绕太阳飞行，所构成等边三角形的臂长L=2.5×106km，编队三角形中心与地球相差约20°。本算例中待优化参数范围如表5所示，遗传算法参数配置与地心轨道编队算例相同。将优化结果表示为位置速度后如表6所示。

表5 日心轨道编队待优化参数及取值范围Table 5 Heliocentric orbit formation parameters to be optimized and their value ranges

表6 日心轨道编队优化结果Table 6 Heliocentric orbital formations optimization results

图5是日心轨道编队航天器SC1、SC2和SC3之间的干涉臂长变化率曲线，图中绿色曲线表示干涉臂长L12的变化率，最大值∆L12为0.48%；红色曲线表示干涉臂长L13的变化率，最大值∆L13为0.25%；蓝色曲线表示干涉臂长L23的变化率，最大值∆L23为0.32%。

图5 日心轨道编队航天器之间的臂长变化率Fig.5 Heliocentric orbit formation arm length variation between spacecraft

图6是日心轨道编队航天器SC1、SC2和SC3轨道构型呼吸角的变化曲线，图中绿色曲线表示呼吸角θ1，变化量最大值∆θ1为0.36°；红色曲线表示呼吸角θ2，变化量最大值∆ θ2为0.38°；蓝色曲线表示呼吸角θ3，变化量最大值∆ θ3为0.26°。

图6 日心轨道编队航天器呼吸角Fig.6 Heliocentric orbit formation spacecraft breathing angle

图7是日心轨道编队航天器SC1、SC2和SC3相对速度变化曲线，图中绿色曲线表示SC1和SC2之间的相对速度，变化量最大值∆v1为1.82 m/s；红色曲线表示SC1和SC3之间的相对速度，变化量最大值∆v2为1.05 m/s；蓝色曲线表示SC2和SC3之间的相对速度，变化量最大值∆v3为0.92 m/s。

图7 日心轨道编队航天器之间的相对速度Fig.7 Heliocentric orbit formation relative velocity between spacecraft

上述结果表明，在考虑了太阳引力作用下，日心轨道空间引力波探测系统在一年内的变化能满足引力波探测的需求。日心编队构型轨道频率变化缓慢，可以在较长时间保持稳定的构型。

以优化得到的结果作为初始条件，用Monte-Carlo法分析入轨误差对日心轨道编队构型的影响。分别向航天器的径向、切向和法向位置增加标准差为200 km的位置误差，对构型的影响如表7所示。分别向航天器的径向、切向和法向速度增加标准差为1 cm/s的速度误差，对构型的影响如表8所示。位置误差和速度误差均符合正态分布。可以看出，径向位置误差和切向速度误差导致构型标准差的变化最大。若位置误差和速度误差同时存在，当速度误差为3 cm/s时，位置误差最大不超过90 km，否则构型将会发散。

表7 位置误差对日心轨道编队构型的影响Table 7 Effect of position error on configuration of heliocentric orbital formation

表8 速度误差对日心轨道编队构型的影响Table 8 Effect of velocity error on configuration of heliocentric orbital formation

地心轨道编队和日心轨道编队的优化结果对比如表9所示，地心轨道编队的臂长变化率最大值为0.18%，呼吸角变化量最大值为0.14°，相对速度的变化量最大值为6.70 m/s；日心轨道编队的臂长变化率最大值为0.48%，呼吸角变化量最大值为0.38°，相对速度的变化量最大值为1.82 m/s。地心轨道编队和日心轨道编队的优化结果均满足各自的指标要求，入轨误差中对构型稳定性影响较大的均为径向位置误差和切向速度误差。

表9 仿真算例数据对比Table 9 Comparision of simuilation results

4 结 论

本文以空间引力波探测器编队构型设计为研究背景，重点研究了探测器动力学建模和轨道优化算法。仿真结果表明，所提优化算法可以使引力波探测器在持续一年的时间里，保持构型的稳定性，其中，无中心天体的日心轨道探测器编队的臂长变化率保持在0.48 %以内，呼吸角最大变化量保持在0.38%以内；地心轨道探测器编队的相对速度变化量保持在6.70 m/s以内。另外，径向位置误差和切向速度误差对编队构型稳定性的影响较大。上述结果说明，所提基于对偶四元数编队动力学模型，可以对探测器有无中心天体的情况进行统一描述，在保证动力学模型精度的同时提升模型的适用性。本文的研究成果能够为未来开展空间引力波探测提供建议与参考。