空间激光干涉引力波探测器由3颗航天器组成，3颗航天器以120°的相位差分布在同一圆轨道上，形成尺度十万公里到百万公里的等边三角形构型。为了实现高精度的引力波探测，探测器的等边三角形构型必须具有高度的稳定性。文献[11-12]分别采用粒子群算法和组合优化方法设计了满足等边三角形构型稳定性的天琴计划任务轨道。文献[13]采用两步优化策略对LISA编队的滞后角、呼吸角变化量和臂长变化率进行了数值优化，并分析了轨道周期偏差对构型稳定性的影响。文献[14]采用遗传算法对空间太极计划航天器任务轨道进行了优化设计。文献[15]以摄动轨道平均周期和平均偏心率为优化指标对ASTROD-GW(ASTROD optimized for Gravitation Wave detection)计划任务轨道进行了优化设计。入轨误差和摄动会造成航天器的实际轨道和名义轨道存在偏差，这种偏差会导致等边三角形构型破坏，无法满足任务要求。因此，必须通过轨道修正将航天器的轨道偏差控制在一定范围内。现有的文献尚未对空间激光干涉引力波探测器的轨道偏差修正问题进行研究。
针对空间激光干涉引力波探测器轨道修正问题，提出一种基于虚拟编队构型设计的航天器轨道修正方法。本文首先对航天器受到的主要摄动进行量级分析，建立平均轨道要素摄动方程。然后，根据探测器等边三角形构型稳定性约束对虚拟编队的构型进行设计，并考虑摄动的影响, 利用平均轨道要素摄动方程对平均半长轴进行小量偏置，求解平均轨道要素修正量和相应的轨道修正速度脉冲。最后，通过数值仿真验证该轨道修正方法的有效性。
1 平均轨道要素摄动方程 1.1 问题描述 天琴计划如图 1所示, 图中：S为太阳；E为地球；S1、S2、S3为3颗航天器。
 |
| 图 1 天琴计划轨道构型示意图 Fig. 1 Schematic of orbital configuration for TianQin |
| 图选项 |
该计划在距离地球约10万公里高度的同一圆轨道面内放置3颗航天器，各航天器相位角相差120°，形成等边三角形构型，轨道面法向指向白矮星双星系统RX J0806.3+1527，任务寿命为3~5年[9-12]。在地心赤道惯性坐标系下，天琴计划航天器名义轨道要素如表 1所示[11]。表中：a、e、i、Ω、ω和M分别为航天器的半长轴、偏心率、轨道倾角、升交点赤经、近地点幅角和平近点角。
表 1 天琴计划航天器名义轨道要素[11] Table 1 Nominal orbital elements of spacecraft for TianQin[11]
| 航天器 | a/km | e | i/(°) | Ω/(°) | ω/(°) | M/(°) |
| S1 | 100 000 | 0 | 74.39 | 211.58 | 0 | 0 |
| S2 | 100 000 | 0 | 74.39 | 211.58 | 0 | 120 |
| S3 | 100 000 | 0 | 74.39 | 211.58 | 0 | 240 |
表选项
定义两颗航天器之间的距离为臂长、视线相对速度为臂长变化率、2个臂长之间的夹角为呼吸角，天琴计划等边三角形构型稳定性指标如表 2所示[11]。表中：ΔL为臂长变化量；Δθ为呼吸角变化量；
表 2 天琴计划等边三角形构型稳定性指标[11] Table 2 Stability index of equilateral triangular configuration for TianQin[11]
| 参数 | ΔL/km | Δθ/(°) |  /(m·s-1) |
| 大小要求 | < 1 500 | < 1.5 | < 9 |
表选项
设ri(t)、vi(t)分别为t时刻航天器Si在地心赤道惯性坐标系下的位置、速度矢量，根据前述定义ΔL、Δθ、
 | (1) |
 | (2) |
式中：rij(t)=ri(t)-rj(t)；rik(t)=ri(t)-rk(t)；vij(t)=vi(t)-vj(t)；名义臂长
1.2 摄动方程 由表 1可知，天琴计划航天器的轨道高度约为10万公里，受到的主要摄动为地球J2摄动、日月引力摄动和太阳光压摄动。由于航天器采用无拖曳控制，可消除非保守力，因此太阳光压不予考虑[1]。
地球J2摄动加速度和地球中心引力加速度的大小之比满足：
 | (3) |
式中：Fo为地球中心引力加速度；FJ2为地球J2摄动加速度；地球二阶带谐系数J2=1.082 63×10-3；地球赤道半径rE=6 378.140 km；r为航天器到地心的距离。
日月引力摄动为第三体摄动，第三体摄动加速度和地球中心引力加速度的大小之比满足：
 | (4) |
式中：Fε为第三体摄动加速度；με为第三天体的引力常数；地球引力常数μe=3.986 005×1014 m3/s2；rε为第三天体到地心的距离。
月球引力常数μm=4.902 8×1012 m3/s2；太阳引力常数μs=1.327 124 4×1020 m3/s2；取月球到地球的距离为地月平均距离，rm=384 403.9 km；取太阳到地球的距离为日地平均距离，rs=1 AU(AU为天文单位，指地球和太阳之间的平均距离)。根据式(3)和式(4)计算，得到天琴计划航天器受到的各摄动量级如下：
 | (5) |
式中：Fm为月球引力摄动加速度；Fs为太阳引力摄动加速度。
由式(5)可知，日月引力摄动的量级相近，远大于地球J2摄动。因此，忽略地球J2摄动，主要考虑日月引力摄动。对于第三体摄动，存在相应的摄动势函数，并且主要部分为二阶项[16]。略去高阶项，在地心赤道惯性坐标系下，航天器的平均摄动势函数可描述为[17]
 | (6) |
式中：
式(6)是采用平均化思想，由日月引力摄动下的二阶摄动势函数对航天器轨道平近点角和日月轨道平近点角进行积分平均化得到的。因此，式(6)不包含短周期项和与日月第三体运动相关的长周期项。
将式(6)代入拉格朗日摄动方程得到航天器平均轨道要素摄动方程为[17]
 | (7) |
式中：n=(μe/a3)1/2；η2=1-e2。
由式(7)可知，在日月引力摄动作用下，航天器平均轨道要素Ω、ω、M发生长期和长周期变化，e、i发生长周期变化，a不受影响。
2 轨道修正 2.1 虚拟编队构型设计 “虚拟编队”的概念由文献[18]提出。定义始终运行在名义轨道上的航天器为理想航天器，运行在实际轨道上的航天器为真实航天器。真实航天器与理想航天器的平均轨道要素之差可视为小量，以理想航天器为参考航天器，真实航天器为绕飞航天器，组成虚拟编队。天琴计划的3颗航天器分别与对应的理想航天器组成3个虚拟编队，如图 2所示。图中：S′1、S′2、S′3为理想航天器，构成理想的等边三角形构型；S1、S2、S3为真实航天器，构成真实的三角形构型；Lij为真实航天器Si和Sj构成的臂长；θi为真实航天器Si对应的呼吸角；真实航天器S1、S2、S3分别与理想航天器S′1、S′2、S′3组成3个虚拟编队，虚拟编队构型的包络半径为Δr。
 |
| 图 2 天琴计划虚拟编队示意图 Fig. 2 Schematic of virtual formation for TianQin |
| 图选项 |
虚拟编队构型包络半径Δr可由相对运动关系得到。设轨道修正后的真实航天器平均轨道要素为(a+δa, e+δe, i+δi, ω+δω, Ω+δΩ, M+δM)，在理想航天器轨道坐标系下修正后的真实航天器相对位置分量表达式为[19]
 | (8) |
式中：A、C与虚拟编队构型的大小和形状相关；B与虚拟编队构型的中心相关；t为时间; φ、ψ为相位角参数。
 | (9) |
值得注意的是，式(8)是基于二体模型和圆轨道得到的，针对密切轨道要素。但在一阶近似的情况下，式(8)在有摄动的情况下也可适用于平均轨道要素。
为保证绕飞轨道闭合，有
 | (10) |
为使理想航天器位于虚拟编队构型中心，有
 | (11) |
虚拟编队构型的包络半径Δr可表示为
 | (12) |
理想航天器运行于圆轨道时，将式(9)代入式(12)得
 | (13) |
由表 2可知，天琴计划对3颗航天器组成的等边三角形构型稳定性有严格的要求，其中臂长变化量、呼吸角变化量与虚拟编队构型大小直接相关，因此必须对虚拟编队构型的包络半径Δr进行设计。当真实航天器与理想航天器的平均轨道要素之差为小量时，真实航天器和理想航天器组成的虚拟编队构型包络半径远小于2个理想航天器之间的距离。此时，可将真实航天器投影到理想航天器所在平面，用投影后的平面几何关系近似真实航天器之间的空间几何关系，如图 3所示。图中：真实航天器S1、S2、S3投影到理想航天器所在平面后位于虚线圆内的任意位置；θmax为最大呼吸角；θmin为最小呼吸角。
 |
| 图 3 天琴计划虚拟编队构型几何关系示意图 Fig. 3 Schematic of geometric relationship of virtual formation configuration for TianQin |
| 图选项 |
由几何关系得
 | (14) |
式中：
 |
 |
 | (15) |
式(15)根据等边三角形构型臂长变化量和呼吸角变化量给出了虚拟编队构型包络半径Δr的约束条件。在不考虑摄动的情况下，虚拟编队的构型应同时满足式(10)、式(11)、式(13)、式(15)，求解可得修正后的真实航天器轨道要素。由于约束条件少于变量个数，存在一定的优化空间减小燃料消耗，可通过松散控制实现轨道修正[20]。
2.2 日月引力摄动对虚拟编队的影响及修正 当真实航天器和理想航天器的平均轨道要素之差为小量时，真实航天器和理想航天器的平均轨道要素变化率之差可用一阶变分表示[21]。对式(7)取一阶变分得
 | (16) |
由式(16)可知，在日月引力摄动作用下，δe、δi发生长周期变化，δΩ、δω、δM发生长期变化和长周期变化。
 | (17) |
式中：
 |
为使虚拟编队的虚拟中心相对于理想航天器不发生沿迹向漂移[19]，由式(11)得
 | (18) |
在偏心率为0的情况下，将式(17)代入式(18)化简得
 | (19) |
式中：
 |
根据式(19)对真实航天器平均半长轴进行小量偏置，可以抑制虚拟编队的虚拟中心的沿迹向漂移。式(19)与式(10)不一致的原因是式(10)是在二体模型下给出的，未考虑摄动的影响。考虑日月引力摄动后，修正后的平均轨道要素可由式(11)、式(13)、式(15)、式(19)求出。
2.3 轨道修正速度脉冲 设真实航天器平均轨道要素为(a+σa, e+σe, i+σi, Ω+σΩ, ω+σω, M+σM)，修正后为(a+δa, e+δe, i+δi, Ω+δΩ, ω+δω, M+δM)，修正量为(Δa, Δe, Δi, ΔΩ, Δω, ΔM)，则：Δa=δa-σa，Δe=δe-σe，Δi=δi-σi，ΔΩ=δΩ-σΩ，Δω=δω-σω，ΔM=δM-σM。
采用四脉冲控制实现轨道修正[22]。
1) 在纬度幅角u=uh处施加法向脉冲ΔVh，有
 | (20) |
 | (21) |
2) 在纬度幅角u=ur处施加径向脉冲ΔVr1，有
 | (22) |
 | (23) |
3) 在纬度幅角u=ur+π处同时施加迹向脉冲ΔVt和径向脉冲ΔVr2，有
 | (24) |
 | (25) |
式中：
3 数值仿真 假设由于入轨误差和摄动等因素，真实航天器和理想航天器的平均轨道要素偏差如表 3所示。
表 3 平均轨道要素偏差 Table 3 Deviation of mean orbital elements
| 航天器 | σa/km | σe | σi/(°) | σΩ/(°) | σω/(°) | σM/(°) |
| S1 | 7 | 0.001 | -0.04 | 0.03 | 0.03 | -0.05 |
| S2 | 5 | 0.001 | 0.04 | -0.03 | -0.02 | 0.04 |
| S3 | -8 | 0.001 | 0.05 | 0.02 | -0.04 | 0.03 |
表选项
利用平均轨道要素摄动方程式(7)对真实航天器的轨道进行数值仿真。在理想航天器轨道坐标系下，虚拟编队构型如图 4所示。等边三角形构型稳定性参数变化曲线如图 5所示。
 |
| 图 4 天琴计划虚拟编队构型(30 d) Fig. 4 Virtual formation configuration for TianQin (30 d) |
| 图选项 |
 |
| 图 5 天琴计划等边三角形构型稳定性参数变化曲线 Fig. 5 Curves of stability parameters of equilateral triangular configuration for TianQin |
| 图选项 |
由图 4可知，在30 d内，虚拟编队构型存在明显的发散，构型半径量级由105 m变为106 m。由图 5可知，天琴计划航天器等边三角形构型臂长和呼吸角存在明显的长期漂移，臂长变化量和呼吸角变化量不满足任务要求。
根据第2节给出的航天器轨道修正方法，考虑日月引力摄动，求得一组平均轨道要素修正量如表 4所示，所需速度脉冲如表 5所示。由表 3、表 4可知，平均轨道要素修正量小于偏差量，是对探测器轨道偏差的部分修正；并且，半长轴存在小量偏置。
表 4 平均轨道要素修正量 Table 4 Correction value of mean orbital elements
| 航天器 | Δa/km | Δe | Δi/(°) | ΔΩ/(°) | Δω/(°) | ΔM/(°) |
| S1 | -7.001 4 | 0 | 0 | 0 | 0 | 0.011 9 |
| S2 | -4.998 6 | 0 | 0 | 0 | 0 | -0.011 9 |
| S3 | 8.001 8 | 0 | 0 | 0 | 0 | 0.004 6 |
表选项
表 5 轨道修正速度脉冲 Table 5 Velocity pulse for orbit correction
| 航天器 | (ΔVh, uh) | (ΔVr1, ur) | (ΔVt, ΔVr2, ur+π) |
| S1 | (0 m/s, -) | (-0.103 9 m/s, -90°) | (-0.069 9 m/s, -0.103 9 m/s, 90°) |
| S2 | (0 m/s, -) | (0.103 9 m/s, -90°) | (-0.049 9 m/s, 0.103 9 m/s, 90°) |
| S3 | (0 m/s, -) | (-0.040 2 m/s, -90°) | (0.079 9 m/s, -0.040 2 m/s, 90°) |
表选项
采用相同的仿真条件进行数值仿真。在理想航天器轨道坐标系下，轨道修正后的虚拟编队构型如图 6所示。相应的等边三角形构型稳定性参数变化曲线如图 7所示。
 |
| 图 6 轨道修正后的天琴计划虚拟编队构型(30 d) Fig. 6 Virtual formation configuration for TianQin after orbit correction (30 d) |
| 图选项 |
 |
| 图 7 轨道修正后的天琴计划等边三角形构型稳定性参数变化曲线 Fig. 7 Curves of stability parameters of equilateral triangular configuration for TianQin after orbit correction |
| 图选项 |
由图 6可知，经过轨道修正，在30 d内虚拟编队构型呈现为稳定的八字形空间封闭曲线，发散程度大大减小，构型半径量级保持为105 m。由图 7可知，天琴计划航天器等边三角形构型臂长和呼吸角的长期漂移得到有效抑制，臂长变化量和呼吸角变化量满足任务要求。
4 结论 对空间激光干涉引力波探测器轨道偏差修正问题进行了研究，提出了一种基于虚拟编队构型设计的航天器轨道修正方法。
1) 对探测器受到的主要摄动进行量级分析，建立了日月引力摄动下的航天器平均轨道要素摄动方程。
2) 根据探测器等边三角形构型稳定性指标对虚拟编队构型大小进行了设计。
3) 通过平均半长轴的小量偏置抑制了摄动带来的虚拟编队虚拟中心沿迹向漂移。求解出的航天器平均轨道要素修正量小于偏差量。轨道修正采用四脉冲控制实现。
4) 数值仿真结果表明，采用该方法对探测器轨道进行部分修正，有效抑制了探测器等边三角形构型臂长和呼吸角的长期漂移，满足任务要求，具有减少燃料消耗、延长任务寿命的潜力。
本文未考虑短周期项、与日月第三体运动相关的长周期项的影响，以及求解航天器平均轨道要素修正量时的燃料消耗问题。即，如何在完整力模型下实现燃料最优轨道修正有待进一步研究。
参考文献
| [1] | 罗子人, 白姗, 边星, 等. 空间激光干涉引力波探测[J]. 力学进展, 2013, 43(4): 415-447. LUO Z R, BAI S, BIAN X, et al. Space laser interferometry gravitational wave detection[J]. Advances in Mechanics, 2013, 43(4): 415-447. (in Chinese) |
| [2] | WEBER J. Observation of the thermal fluctuations of a gravitational-wave detector[J]. Physical Review Letters, 1966, 17(24): 1228-1230. DOI:10.1103/PhysRevLett.17.1228 |
| [3] | BLAIR D, JU L, ZHAO C N, et al. Gravitational wave astronomy:The current status[J]. Science China Physics Mechanics & Astronomy, 2015, 58(12): 120402. |
| [4] | ABBOTT B P, ABBOTT R, ABBOTT T D, et al. Observation of gravitational waves from a binary black hole merger[J]. Physical Review Letters, 2016, 116(6): 061102. DOI:10.1103/PhysRevLett.116.061102 |
| [5] | CORNELISSE J W. LISA mission and system design[J]. Classical and Quantum Gravity, 1996, 13(11A): A251-A258. DOI:10.1088/0264-9381/13/11A/034 |
| [6] | ARMANO M, AUDLEY H, BAIRD J, et al. Beyond the required LISA free-fall performance:New LISA pathfinder results down to 20 μHz[J]. Physical Review Letters, 2018, 120(6): 061101. DOI:10.1103/PhysRevLett.120.061101 |
| [7] | WU A M, NI W T. Deployment and simulation of the ASTROD-GW formation[J]. International Journal of Modern Physics D, 2013, 22(1): 1341005. DOI:10.1142/S0218271813410058 |
| [8] | JIN G.Program in space detection of gravitational wave in Chinese Academy of Sciences[C]//11th International LISA Symposium.Bristol: IOP Publishing, 2017, 840(1): 012009. |
| [9] | LUO J, CHEN L S, DUAN H Z, et al. TianQin:A space-borne gravitational wave detector[J]. Classical and Quantum Gravity, 2016, 33(3): 035010. DOI:10.1088/0264-9381/33/3/035010 |
| [10] | HU X C, LI X H, WANG Y, et al. Fundamentals of the orbit and response for TianQin[J]. Classical and Quantum Gravity, 2018, 35(9): 095008. DOI:10.1088/1361-6382/aab52f |
| [11] | 万小波, 张晓敏, 黎明. 天琴计划轨道构型长期漂移特性分析[J]. 中国空间科学技术, 2017, 37(3): 110-116. WAN X B, ZHANG X M, LI M. Analysis of long-period drift characteristics for orbit configuration of TianQin Mission[J]. Chinese Space Science and Technology, 2017, 37(3): 110-116. (in Chinese) |
| [12] | YE B B, ZHANG X, ZHOU M Y, et al. Optimizing orbits for TianQin[J]. International Journal of Modern Physics D, 2019, 28(9): 1950121. DOI:10.1142/S0218271819501219 |
| [13] | YANG C, ZHNAG H. Formation flight design for a LISA-like gravitational wave observatory via Cascade optimization[J]. Astrodynamics, 2019, 3(2): 155-171. DOI:10.1007/s42064-018-0042-9 |
| [14] | 唐文林.中国空间引力波探测计划卫星轨道设计的初步研究[D].北京: 中国科学院大学, 2014: 17-31. TANG W L, Preliminary study of the orbit design of the Chinese mission to detect the gravitational wave in space[D].Beijing: University of Chinese Academy of Sciences, 2014: 17-31(in Chinese). |
| [15] | WANG G, NI W T. Orbit optimization for ASTROD-GW and its time delay interferometry with two arms using CGC ephemeris[J]. Chinese Physics B, 2013, 22(4): 571-579. |
| [16] | 刘林, 汤靖师. 卫星轨道理论与应用[M]. 北京: 电子工业出版社, 2015: 227-243. LIU L, TANG J S. Theory and application of satellite orbit[M]. Beijing: Publishing House of Electronics Industry, 2015: 227-243. (in Chinese) |
| [17] | PRADO B D A, FERNANDO A. Third-body perturbation in orbits around natural satellites[J]. Journal of Guidance, Control, and Dynamics, 2003, 26(1): 33-40. |
| [18] | 杜耀珂, 杨盛庆, 完备, 等. 近地卫星严格回归轨道保持控制[J]. 航空学报, 2018, 39(12): 322449. DU Y K, YANG S Q, WAN B, et al. Strictly-regressive orbit maintenance control of near earth satellites[J]. Acta Aeronautica et Astronautica Sinica, 2018, 39(12): 322449. (in Chinese) |
| [19] | 孟云鹤. 航天器编队飞行导论[M]. 北京: 国防工业出版社, 2014: 28-32. MENG Y H. Introduction to spacecraft formation flying[M]. Beijing: National Defense Industry Press, 2014: 28-32. (in Chinese) |
| [20] | HOU X Y, ZHAO Y H, LIU L. Formation flying in elliptic orbits with the J2 perturbation[J]. Research in Astronomy and Astrophysics, 2012, 12(11): 1563. DOI:10.1088/1674-4527/12/11/010 |
| [21] | ROSCOE C W T, VADALI S R, ALFRIEND K T. Third-body perturbation effects on satellite formations[J]. Journal of the Astronautical Sciences, 2013, 60(3-4): 408-433. DOI:10.1007/s40295-015-0057-x |
| [22] | 曹喜滨, 张锦绣, 王峰. 航天器编队动力学与控制[M]. 北京: 国防工业出版社, 2013: 103-121. CAO X B, ZHANF J X, WANG F. The dynamics and control of spacecraft formation flying[M]. Beijing: National Defense Industry Press, 2013: 103-121. (in Chinese) |
