基于二阶CW方程的引力波探测编队构形设计

焦博涵，党朝辉

（1.西北工业大学 航天学院，西安 710072；2.西北工业大学 航天飞行动力学技术重点实验室，西安 710072）

引 言

2015年9月14日激光干涉引力波天文台（Laser Interferometer Gravitational-Wave Observatory，LIGO）合作小组第一次探测到了引力波信号，这一发现填补了广义相对论实验验证的最后一块缺失的拼图，同时也揭开了人们利用引力波进行宇宙研究的序幕。而随着研究的深入，传统地面引力波探测装置由于地面噪声、臂长等限制的存在已越发不能满足研究所需的探测需求[1]。为此一系列采用卫星编队进行引力波探测的空间任务被陆续提出。欧洲航天局（European Space Agency，ESA）与美国国家航空航天局（National Aeronautics and Space Administration，NASA）于1993年最早提出了空间引力波探测计划LISA[2]项目。随后，日本、欧洲等也先后提出了类似的计划如DECIGO[3]、BBO[4]、ALIA[5]等。中国从2008年开始探讨空间引力波探测的可行性，先后开展了概念与方案设计、关键科学载荷研制等，并于2014年和2016年提出“天琴计划”[6]和“太极计划”[7]。

空间引力波探测任务是利用3个航天器构成空间正三角形编队，根据迈克尔逊激光干涉原理，通过测量相邻两个航天器间的臂长变化来表征引力波信号[8-9]。因此如何设计一个可以长期稳定保持的正三角形编队便成为了空间引力波探测任务的重要问题，对此相关学者进行了大量的研究并给出了3类设计方法[9]：①共轨星座方式，即将3个航天器均匀地布置在同一条围绕中心天体的圆形轨道上；②相对绕飞方式，即将3个航天器均匀地布置在圆参考轨道附近的相对圆形绕飞轨道上，需要额外说明的是：这类方法在具体进行设计时又可分为基于绝对运动的设计方法（通过确定一个航天器的绝对运动轨道并根据正三角形编队的几何关系，将其绕垂直轨道平面的方向旋转，得到另外两个航天器的轨道[10-12]）和基于相对运动的设计方法〔利用CW（Clohessy-Wiltshire）方程构造参考轨道附近的空间圆形相对轨道，并将3个航天器均匀布置在该相对轨道上从而构成三角形编队[12-13]〕；③三角平动点方式，即将3个航天器分别布置在日地圆形限制性三体模型中的L3、L4、L5共3个平动点上。针对上述3类设计方法，由于基于三角平动点方式的设计方法的工程实现难度较大，因此对于绝大多数的空间引力波探测任务通常采用前两类设计方法。但不论是共轨星座方式还是相对绕飞方式，其均为理想化二体引力场下（无摄动或线性引力场）的设计结果。因此会导致编队在高精度引力场运行过程中出现较为严重发散情况（具体表现在正三角形编队臂长、臂长变化率、呼吸角等指标的发散上）[14-16]，从而无法满足引力波探测任务所需的编队稳定性。

为解决上述编队构形长时间运行的稳定性问题，有不同研究提出了若干解决方法。针对这些方法，从其所依据的模型来看，可将其分为基于Relative Orbital Elements（ROE）方程[17]的编队设计、基于绝对轨道动力学方程的编队设计以及基于二阶CW方程[10]的编队设计共3类。①基于ROE方程的编队设计：ROE方程是指基于相对轨道根数差的相对运动方程。清华大学蒋方华副教授团队[18]基于ROE方程，通过3个约束条件，将决策变量由18个减少到14个，并提出了一种变量区域自适应调整算法，实现了高效优化，而Joffre等[19]则进一步考虑了太阳系内主要天体的摄动影响。中国科学院空间应用工程与技术中心张皓研究员[20]在ROE方程的基础上提出了基于非奇异轨道根数的相对运动方程，并通过分析优化指标间的相关性将多目标优化问题转化为了多个单目标优化问题。②基于绝对轨道动力学方程：李卓[21]在二体模型的基础上，考虑多种摄动因素，采用辛积分方法实现了构形优化。易照华[22]着重关注LISA任务编队中心与地球同轨的特点，建立了共轨限制性三体问题，并推导了共轨限制性三体模型下的编队臂长的表达式，但却并未研究如何利用该表达式进行编队优化。③基于二阶CW方程的编队设计：Nayak[10]基于二阶CW方程推导了编队臂长与编队平面倾角的关系，发现可通过改变编队平面倾角来降低编队臂长的波动峰值。Pucacco[23]和Dhurandhar[24]在二阶CW方程的基础上考虑不同的摄动因素进行了进一步的优化。虽然Nayak等基于二阶CW方程对编队构形进行了优化设计，但其仍存在以下两点不足：①Nayak所得到的是一种具有二阶精度的CW方程周期解而并非二阶CW方程的解析解，因此并不能直接地解释构形发散地深层机理；②Nayak的方法计算较为复杂且不直观，无法解释更优稳定解的物理机制，且构型稳定性的优化结果仍然有限。

针对Nayak方法的局限性，本文提出了一种基于二阶CW方程的编队构形设计方法：基于摄动法推导了二阶CW方程的近似解析解，得益于这种解析解，成功解释了编队构形发散的主要原因，并依据分析结果给出了一种以航天器相位角为优化变量的构形优化方法。该方法相较于Nayak的方法具有更好的优化结果，且在优化臂长指标的基础上能够同时实现对于呼吸角指标的优化。

1 空间引力波探测编队标称构形设计

本小节首先回顾太极计划对于编队构形设计的主要要求。其次，给出描述相对运动的CW方程，并在此基础上通过构建空间圆形绕飞轨道实现标称构形的设计。

1.1 太极引力波探测计划

作为中国的空间引力波探测计划，太极计划于2016年由中国科学院正式提出，旨在探测频段位于0.1～1 mHz范围内的引力波源。太极任务由3颗完全相同的航天器组成，它们以 20◦的角度沿超前或落后地球的日心轨道运行。如图1所示，3颗航天器形成一个臂长约300万km的等边三角形。

图1 太极引力波探测计划Fig.1 Taiji gravitational wave detection program

与LISA计划类似，太极计划同样利用每颗航天器上所携带的一对激光测距干涉仪和两个重力参考传感器，基于正三角形编队构形形成三组迈克尔逊干涉仪实现对于引力波信号的探测。在实际工程中，若要实现上述原理的引力波探测任务，则需要对星间所形成的三组迈克尔逊干涉仪提出极高的精度要求，即需要对正三角形编队的整体稳定性提出一定的要求。

文献[21]指出了任务过程中对于编队的稳定性要求（示意图见图2），其具体指标要求见表1。在本文中，从简化问题的角度出发，将着重关注臂长变化这一指标。

表1 “太极”任务编队稳定性指标Table 1 Formation stability index of Taiji mission

图2 编队稳定性指标示意图Fig.2 Diagram of formation stability index

1.2 参考坐标系

本文采用如下参考坐标系用于描述空间引力波探测编队的运动：

1）日心惯性坐标系

本文日心惯性坐标系的建立参考国际天体参考系[25]（International Celestial Reference System），其坐标原点位于太阳的质心，X轴位于黄道面内并指向J2000春分点，Z轴与黄道面垂直，与地球公转速度方向一致，Y轴位于黄道面内，与X、Z轴构成右手坐标系。

2）Hill坐标系

Hill坐标系常用于描述相对运动，其坐标原点位于正三角形编队的虚拟中心，x轴位于虚拟中心的轨道平面内并由日心指向虚拟中心，z轴垂直于虚拟中心的轨道平面并指向角动量方向，y轴分别与另外两轴垂直从而构成右手坐标系。该坐标系也常称为当地垂直当地水平（LVLH）坐标系。

上述两坐标系如图3所示。

图3 坐标系示意图Fig.3 Diagram of coordinate system

1.3 编队运动方程

在不考虑摄动力与控制力的情况下，若假设：①参考航天器的运动轨道为圆形轨道；②从航天器与参考航天器之间的距离远小于其各自的轨道半径，则二体引力场中的线性相对运动动力学方程可在Hill坐标系中表示为：

在不考虑摄动力与控制力的情况下，由1.1节可知：①太极计划的编队虚拟中心位于偏心率接近为0的地球公转轨道；②编队虚拟中心距离太阳的距离为1 AU远大于3×106km的编队尺度。因此针对太极计划编队的上述两个特点，可采用如下建立在Hill坐标系中的CW方程来描述整个编队的运动[26]

其中：x0、y0、z0、、、分别表示t=0时刻相对位置和相对速度的值。

1.4 基于CW方程的标称构形设计

在长时间的任务周期中，为了保持编队整体的有界性，必须消除式（2）中的长期项，由此可得到CW方程的周期性运动条件为

将该条件代入到式（2）中便可得到CW方程的周期解为

正三角形编队构形设计的关键在于空间圆形相对轨道的构建。而所谓空间圆形是指从航天器相对于参考航天器间的距离为一定值，因此有

进而将式（4）代入便可得空间圆形编队的初始条件为

由此若要获得正三角形编队，则只需要在上式所确定的空间圆形相对轨道上布置三个相位角相差120◦的航天器即可，即

而式（6）中z方向与x方向在相对位置和相对速度上的正负比例关系则代表了顺时针旋转和逆时针旋转两类构形模式[18-19]。

综上，式（6）和（7）即为基于CW方程所设计的正三角形标称构形。

2 二阶CW方程及其近似解析解

在上一节中介绍了描述相对运动的CW方程，其是在非线性相对运动动力学的基础上，将非线性的引力梯度项进行Taylor级数展开并保留一阶项所得到的线性方程。若对非线性的引力梯度项的Taylor级数保留至二阶项，则可获得式（8）所示的二阶CW方程

由于二阶CW方程相较于CW方程保留了一部分二阶非线性项，因此其具有更高的精度，能够体现出更多的非线性规律。为了充分利用二阶CW方程中所蕴含的非线性信息，望求取式（8）的解析解。因式（8）为一组非线性微分方程组，无法直接写出其通解的形式，因此考虑使用摄动法求解式（8）的近似解析解。

摄动法（或称小参数法）是一种将非线性因素作为对线性系统的摄动，从而在线性解基础上寻求非线性系统近似解的方法[27]。据摄动法核心思想式（8）的解可以写为

其中：ε即为摄动法中的小参数，通过选择保留ε的不同阶次项即可获得不同精度的近似解析解。对于本文所要研究的问题，为避免使问题过度复杂，本文仅选择保留到ε的一次项。

对于式（8）可选取ε=n2/a作为小参数，此时根据式（9），在保留ε的一次项基础上将近似解析解的形式代入到式（8）中，则可将式（8）的求解问题转化为如下两个微分方程组的求解问题

对于式（10）的方程组而言，它的解即为CW方程的解析解。进一步，将式（10）的解代入到式（11）微分方程组的右端，并根据线性常微分方程理论，可求得式（11）的解为

其中：α0∼α8、β0∼β7、γ0∼γ6分别为与x0、y0、z0、、、有关的代数表达式，其具体形式为

进而，根据式（2）、式（12）可得到二阶CW方程的近似解析解为

为了直观地说明二阶CW方程在精度上相较于CW方程的优越性，分别采用CW方程的解析解〔式（2）〕、二阶CW方程的近似解析解〔式（36）〕以及二体非线性方程，以标称构形设计结果为初始条件进行演化计算，并分别绘制CW方程解析解、二阶CW方程近似解析解的演化结果与二体方程的演化结果的差的变化情况，其结果如图4所示：

图4 CW方程与二阶CW方程的精度比较Fig.4 Accuracy comparison between CW equation and second-order CW equation

可以看出，二阶CW方程相较于CW方程具有更好的精度，也因此可以充分利用二阶CW方程的这一优势进行编队设计。

3 编队构形的分析与优化

在第2节中推导并给出了二阶CW方程的近似解析解表达式，本节将基于该近似解析解针对空间引力波探测编队构形的设计与优化问题进行分析。

3.1 完美圆形绕飞轨道不存在性证明

值得说明的是，对于以CW方程为基础的正三角形编队构形优化设计问题，另一种描述为：在不同精度的引力场中是否存在完美的圆形相对绕飞轨道[9]。对于这个问题，尽管是最为简单的二体引力场也是难以回答的，而在获得能够表征一定二体引力场非线性特征的二阶CW方程的基础上，本节尝试在二阶CW方程下回答这一问题。

同基于CW方程的标称构形思路一致，为了构成稳定的空间圆形构形则需要消除式（36）中的长期项和漂移项，这要求如下关系式成立

至此在满足式（37）的约束下得到了二阶CW方程的周期解，进一步需要结合式（5）给出二阶CW方程下形成空间圆形编队所需要满足的条件。先将式（37）代入式（36），随后再将式（36）代入到式（5）中，忽略关于ε的二次项，且为保证最后的结果为一常值需令sin2nt和cos2nt外其余项的系数均为0，由此有下列各式成立

式（37）、（38）的意义在于当两式中的所有表达式同时成立时才可以形成空间圆形编队。此外，由式（37）、（38）可以看出，二阶CW方程的空间圆形编队条件是在CW方程空间圆形编队条件的基础上又增添了若干约束，当忽略关于ε的一次项时，两式便退化为CW方程的空间圆形编队条件，进一步验证了二阶CW方程的正确性。

下面利用式（38）中的表达式进行线性组合可得到如下表达式

进一步化简可以推出

若要式（40）成立则需要满足z0==0，此时考虑到ε为一近似为0的小量，则若式（38）中的η0成立则需要满足x0≈≈0，而这种情况显然不会出现，这说明在考虑 η0、η1、η2、η3、η44个表达式均成立的假设下推出了矛盾的结果，因此可以说明 η0、η1、η2、η3、η4不可能同时成立，而这则说明在二阶CW方程下不存在严格精确的空间圆形绕飞轨道。

3.2 基于二阶CW方程的构形发散性分析

由3.1节可知，在二阶CW方程下不存在严格精确的空间圆形绕飞轨道，因此进一步可以说明在二阶CW方程的精度下无法通过完全解析的方法实现对于正三角形编队的设计。为此，较为自然且直接的思路是如何利用已有的有限解析结果来辅助编队的优化设计。以此为目的，本小节旨在通过二阶CW方程分析基于CW方程的标称构形设计结果在二体非线性引力场下的发散机理，从而为后续的编队优化设计提供依据。

根据式（6），首先考虑周期性条件以及无漂移条件，即将y˙0=−2nx0和y0=/n代入式（36），可得

通过比较式（41）和式（36），并观察包含t的项的变化，可以发现：标称构形的周期性条件消除了二阶CW方程近似解析解中的x方向的发散项，但并未消除y方向的发散项；标称构形的无漂移条件消除了二阶CW方程近似解析解中x和y方向的主要漂移项。

进一步，分析式（6）中x方向与z方向在相对位置、相对速度上的关系，不丧失一般性可将z0=和代入到式（41）中。结果显而易见，同样从包含t的项的变化来看，上述两条件并未对式（41）的形式造成影响。但通过观察y方向发散项系数 β1的具体表达式〔式（42）〕可以发现：除x0==0的情况外，不论何种相对初始状态，其构形在二阶CW方程下均会发散，然而根据式（6）中x0和x的表达式，不会出现x0==0的情况。综合上述分析可以发现基于CW方程的标称构形设计结果在二阶CW方程下必然会出现发散的情况，即二体引力场中的非线性部分会破坏标称构形的周期条件，因此在进行编队优化设计的过程中应着重考虑周期条件的修正

3.3 基于二阶CW方程的编队优化设计

在分析基于CW方程的标称构形设计结果在二阶CW方程下的发散机理的基础上，本节主要研究如何通过对已有的标称设计结果进行修正，从而实现编队的优化设计。

1）周期匹配条件与能量匹配条件

根据3.2节的分析结果，二体引力场下基于CW方程的标称构形设计结果的发散原因是因为二体引力场中的非线性部分破坏了标称构形的周期条件。因此在研究编队的优化设计问题之前，需要针对二体引力场下的相对运动周期性条件进行讨论。

对于二体引力场下的相对运动，其周期性相对运动条件存在两种形式上的表达[28]：一种是由轨道根数所描述的周期匹配条件，要求编队中的从航天器和参考航天器具有相同的半长轴；另一种是由笛卡尔坐标所描述的能量匹配条件，要求从航天器的相对运动初始状态满足式（43），其中：R表示参考航天器的轨道半径。

为了说明上述两个周期性相对运动条件的有效性，图5和图6分别给出了基于周期匹配条件和能量匹配条件修正后太极计划编队臂长的演化情况（注：SC1、SC2、SC3的相位角依次为0◦、120◦、240◦）。

图6 二体引力场下基于能量匹配条件修正后的臂长演化情况Fig.6 Evolution of arm length based on modified energy matching condition in two-body gravitational field

从整体来看，经周期匹配条件和能量匹配条件修正后，与图7相比编队臂长在10年任务周期内不再呈现发散趋势，其演化表现出明显的周期性特征。但具体来看，经周期匹配条件修正后，航天器1与航天器2以及航天器1与航天器3之间的臂长变化最大值已经不满足太极计划的臂长要求；而经能量匹配修正后，3颗航天器两两之间的臂长变化最大值均能够满足任务要求，且臂长波动范围能够保持在较细小的范围之内。

图7 二体引力场下标称构形臂长演化情况Fig.7 Evolution of arm length of nominal configuration in two-body gravitational field

尽管周期匹配条件和能量匹配条件仅在表达方式上有所不同，但其在应用于编队设计的过程中却表现出了不同的效果。通过简单分析可以发现，造成这种现象的原因与标称构形的设计方法有着紧密的联系：由于本文所采用的是基于CW方程的相对运动设计方法，所得到的设计结果均基于笛卡尔坐标所表示。对于周期匹配修正而言，虽然其修正了由于非线性所破环的周期条件，但同时在由笛卡尔坐标向轨道根数进行相互转换的过程中同样破环了式（6）中的其它条件；但对于能量匹配修正，由式（43）可以看出，其仅对式（6）中的周期条件进行了修正而并未对其他条件造成影响。因此，通过上述分析可以得到如下结论：在基于CW方程所设计的标称构形的基础上进行优化设计时，能量匹配条件是一项行之有效的约束条件。

2）编队优化问题的建立

结合上述针对相对运动周期性条件的分析结果，可以进一步构建基于二阶CW方程的编队优化模型。具体地，在优化指标方面，可直接基于式（36）的近似解析解构造如下的臂长变化量的解析表达式

其中：l0表示编队的标称臂长；下标i和j表示航天器的编号，满足i=1,2,3，j=1,2,3且i≠j。与现有文献中常见的直接利用数值方法求解非线性方程进而构造臂长的方法相比，解析的臂长表达式虽然在精度上略有不足，但其在计算效率上却具有显著优势。在约束方面，根据上文对于相对运动周期性条件的分析结果，可直接采用式（43）作为约束条件。

在优化变量方面，通过分析图5和图6的结果可以发现，航天器1/航天器2与航天器1/航天器3的臂长变化始终处于同一量级，且大于航天器2/航天器3的臂长。从几何的角度来看，说明演化过程中的编队构形更趋向于等腰三角形的形状而非正三角形。受这一现象启发，在优化过程中可针对式（6）中的航天器相位角α进行修正（如图8所示），以期望使构形更趋向于正三角形。

图8 优化变量的几何意义Fig.8 Geometric meaning of optimization variables

综上，最终得到基于二阶CW方程的编队优化模型为

其中：[t0,tf]表示任务周期；δα表示优化变量的边界值；非线性函数f代表能量匹配条件，即式（43）。通过求解式（45）的优化问题，可进一步计算出修正后正三角形编队构形条件

4 仿真校验与结果分析

本节以太极计划为任务背景，对3.3节中所构建的优化问题〔式（45）〕进行验证。根据式（45）可知，其本质上是一个双层优化问题：内层是以时间为优化变量，以臂长变化量最大为性能指标的优化问题；外层是以航天器相位角为优化变量，以任务周期内最大臂长变化量最小为性能指标的优化问题。对于内层优化问题，根据图5～图7可以看出：对于不同的初始条件，臂长的变化均存在一定频率的波动现象，对于优化问题而言其意味着存在若干的局部极值点，因此对于内层优化问题应采用全局优化算法进行求解，本文基于MATLAB中的GlobalSearch函数进行求解。对于外层优化问题，由于考虑∆α1=0，其实质上退化为了双变量优化问题，同时考虑到指标函数并非连续函数，因此采用模式搜索算法进行求解。

图9给出了经优化修正后二体引力场下编队臂长的演化情况。通过对比图6可以发现：从总体来看，经优化修正后相比于基于能量匹配的修正方法，3个臂长的平均最大变化量由0.42%降低至0.32%，而3个臂长的最大变化量由0.62%降低至0.44%，且相较于Nayak[10]将臂长最大变化量优化为1%的结果，可以看出本文所提出方法的优势。具体来看，经优化修正后航天器2和航天器3之间臂长的变化幅值得到了减小，而航天器1和航天器2以及航天器1和航天器3之间的臂长变化幅值的减小量并不明显，而这一现象则正好印证了式（45）中选取航天器相位角作为优化变量的有效性。

图9 二体引力场下通过优化修正后的臂长演化情况Fig.9 Evolution of arm length after optimization correction in two-body gravitational field

除考量编队的臂长变化情况外，图10和图11分别给出了基于能量匹配修正和优化修正后编队呼吸角在10年任务周期内的演化情况。可以看出，经优化修正后相较于基于能量匹配的修正方法，编队呼吸角的平均最大变化量由0.79%降低至0.6%，而呼吸角的最大变化量由0.92%降低至0.65%。虽然在式（45）的优化问题中仅将编队臂长最大变化量作为优化指标，但通过修正航天器的相位角将编队构形由等腰三角形修正为正三角形，进而使得呼吸角的变化幅值同样得到了抑制，而这也体现出选择航天器相位角作为优化变量的一项优势。

图10 二体引力场下基于能量匹配条件修正后的呼吸角演化情况Fig.10 Evolution of breathing angel length based on modified energy matching condition in two-body gravitational field

图11 二体引力场下通过优化修正后的呼吸角演化情况Fig.11 Evolution of breathing angel length after optimization correction in two-body gravitational field

5 结 论

本文研究了基于二阶CW方程的空间引力波探测正三角形编队的构形设计问题，利用摄动法给出了二阶CW方程的近似解析解，并基于该解析解分析了基于CW方程设计的标称构形的发散机理，构造了编队构形优化问题，最终得到如下结论：

1）在二阶CW方程的模型精度下，不存在完美的圆形相对绕飞轨道；

2）基于CW方程所设计的标称构形在二体非线性引力场下发散的主要原因是：二体引力场中的非线性部分破坏了标称构形条件中的周期条件；

3）针对基于CW方程所设计的标称构型，能量匹配修正相较于周期匹配修正具有更好的效果；

4）通过构造基于二阶CW方程的编队优化模型，在二体引力场下，将编队臂长10年内的平均误差降低至0.32%，最大误差降低至0.44%，同时通过考虑三个臂长变化的均衡性，以航天器的相位角为优化变量，在仅以臂长为优化指标的基础上实现了对呼吸角的优化。

本文的有关分析方法与仿真结果对太极引力波探测任务的编队设计问题具有一定的借鉴意义。