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

基于对偶四元数的空间引力波探测器稳定构型优化设计张锦绣1，张 谕2，王继河1，杨继坤2，逯振坤2，宋昱岐2（1.中山大学 航空航天学院，深圳 518107；2.中山大学 物理与天文学院，珠海 519082）摘 要：针对空间引力波探测器的轨道构型稳定性，研究了空间大尺度航天器编队稳定构型设计和优化方法。考虑日地月引力场对编队构型的影响，在日心坐标系下，将地球、月球和3颗航天器视为一个编队进行研究；基于对偶四元数建立自然/人造天体动力学模型，对有/无中心天体的空间引力波探测系统模型进行统一的描述；给出了基于遗传算法的待优化参数，并设计了优化目标函数。仿真结果表明：在不进行轨道控制情况下，1年内探测器满足轨道构型的稳定性要求。关键词：航天器编队；轨道构型；对偶四元数；遗传算法中图分类号：TN95 文献标识码：A 文章编号：2096-9287（2023）03-0268-09DOI:10.15982/j.issn.2096-9287.2023.20220085引用格式：张锦绣，张谕，王继河，等.基于对偶四元数的空间引力波探测器稳定构型优化设计J.深空探测学报（中英文），2023，10（3）：268-276.Reference format:ZHANG J X，ZHANG Y，WANG J H，et al.Optimal design of stable configuration of spacegravitational wave detector based on dual quaternionJ.Journal of Deep Space Exploration，2023，10（3）：268-276.引言为印证广义相对论，引力波探测成为当前国际物理和空间科学领域的研究热点。国内外各研究机构围绕空间引力波探测技术展开了深入研究，如：欧洲航天局（European Space Agency，ESA）的“激光干涉测量空间天线”（Laser Interferometer Space Antenna，LISA）/eLISA（evolved LISA）计划1-2，中国的“天琴计划”3、“太极计划”4。上述引力波探测计划大同小异，如图1所示，它们的空间引力波探测器均由3个航天器组成等边三角形构型，C位于构型平面的中心，表示地球或者虚拟的点，航天器围绕太阳周期性运动的同时也围绕C周期性运动。另外，它们都假设空间引力波会造成航天器搭载的检验质量相对位置发生改变，并通过测量该变化量来印证引力波的存在。检验质量之间的相对位置，是通过安装在航天器上的激光干涉测距系统测量得到的。为了保证高精度激光干涉测距系统可以长时间处于稳定对准测量状态，航天器编队需要运行在可以保持长期稳定构型的轨道上。不同的是，LISA/eLISA计划和“太极计划”，构型的中心C为虚拟点；“天琴计划”是以地球为中心天体，构型中心C为地球。空间引力波探测编队构型的稳定性主要由3个指标表征：臂长（任意两个探测器之间的距离）、呼吸角（每两个臂之间的角度）和星间视线相对速度。臂长决定了干涉仪的探测频带，不等臂长将会在激光干涉测量中引入激光频率噪声，影响核心科学载荷的设计；呼吸角的变化影响探测器的指向精度；星间视线相对速度影响超稳时钟和相位计等重要载荷的指标设计。探测器为了成功探测到引力波，需要将3个航天器之间的臂长、呼吸角以及相对速度等三角形轨道构型参数长时间保持在一定范围内。准确的航天器动力学模型和轨道优化算法是设计满足上述轨道构型的两大关键要素。对于如LISA/eLISA计划和“太极计划”这种轨道构型中心为虚拟点的航天器编队方案，多数学者以笛卡 收稿日期：2022-09-17 修回日期：2022-12-28基金项目：广东省基础与应用基础研究重大项目（2019B030302001）SC1SC3SC2CZsYsXsOs太阳 图 1 空间引力波探测器构型示意图Fig.1 Schematic diagram of space gravitational wave detector configuration第 10 卷 第 3 期深 空 探 测 学 报（中英文）Vol.10 No.32023 年 6 月Journal of Deep Space ExplorationJune 2023尔坐标或轨道要素为状态量，在日心坐标系下推导航天器间相对动力学方程，然后将其简化为线性的C-W方程进行轨道设计。但是，当探测器轨道构型中心有中心天体时（例如“天琴计划”），上述方式不再适用。这是因为航天器之间的臂长和航天器到中心天体之间的距离在同一量级，不满足C-W方程的简化条件。针对有中心天体的航天器编队情况，本文以对偶四元数为工具，在J2000日心坐标系下建立航天器之间的相对运动动力学模型。值得注意的是，本文所提方法同样适用于轨道构型中心为虚拟点的航天器编队方案。为了得到满足空间引力波探测任务要求的高精度轨道构型，通常需要采用优化算法对航天器运行轨道进行优化。文献5使用全局优化的禁忌搜索算法对LISA引力波探测轨道进行了优化设计，价值函数有18个自变量，找到了满足引力波探测要求的轨道，并给出了从地球停泊轨道进入引力波探测实验轨道的发射段和分离段的轨道设计。文献6分析了影响LISA编队构型的性能指标（呼吸角、滞后角、相对速度）相对于初始轨道要素的敏感性，将LISA轨道优化问题转化为两步级联单目标优化问题，先优化滞后角，后优化呼吸角，并分析了轨道周期偏差对构型稳定性的影响。针对“天琴计划”轨道构型优化设计，文献7对构型发生长周期漂移的轨道要素进行了分析，采用粒子群算法对3颗航天器的轨道半长轴进行优化，初步满足轨道构型的稳定性指标。文献8采用组合优化方法，讨论了由于轨道进动而导致探测器指向的缓慢长期漂移，设计了满足构型稳定性要求的任务轨道。虽然上述研究满足了引力波探测构型稳定性的需求，但是它们均是在已知初始轨道的基础上，对航天器的轨道根数进行优化，并且轨道优化的起始时刻是确定的。没有考虑初始轨道以外是否存在更加稳定的轨道，使得航天器满足引力波探测构型稳定性的需求。运行在甚高轨道的航天器会受到地球、月球和太阳等主要天体的引力影响，若在轨道优化过程中通过星历得到每一时刻地球、月球和太阳的位置，这严重影响了优化算法的计算效率。因此，有必要以太阳为中心建立惯性坐标系，为地球和月球建立简化但足够精确的动力学模型，从而加快优化速度。把地球、月球和航天器看作刚体来描述它们的运动规律，这实际上是描述六自由的运动问题。而对偶四元数9-13是目前描述刚体运动的最合适的数学工具，这是因为它将平移和旋转统一在同一框架下，同时考虑了二者之间的耦合影响，并且可以提供很高的计算效率。综上所述，本文在J2000日心坐标系下，研究日地月引力场中，空间引力波探测系统稳定构型的全局搜索优化方案。基于对偶四元数，在J2000日心坐标系下建立地球、月球和3个航天器的姿轨一体化动力学模型；分析构型稳定性的表征，确定优化的目标函数和待优化的变量；使用遗传算法对航天器轨道构型进行优化设计，使得地球、月球和3颗航天器在不进行轨道控制的前提下尽可能长时间地保持稳定的构型；通过数学仿真验证优化方法的有效性和普适性。1 日心系编队动力学表征模型 1.1 自然/人造天体动力学方程考虑日地月引力场对引力波探测器的影响，本文在日心系下，将自然天体（地球、月球）和人造天体（3个航天器）视为一个完整的编队。这样可以在建立编队动力学模型的时候，更加方便地对日地月引力进行表述。OsXsYsZs地球在J2000日心坐标系14下的运动学方程为 qe=12 qe wee(1)qe=qe+12qereeqeree wee=wee+(ree+weeree)wee其中：；是地固坐标系相对于日心坐标系的四元数；是地球在日心坐标系下的位置向量在地固坐标系下的表示；运算符 表示对偶四元数和四元数的乘法；速度旋量为，是地球的角速度在地固坐标系下的表示。OsXsYsZs地球在J2000日心坐标系下的动力学方程为 wee=M1eFeM1e(weeMe wee)(2)Fe=fe+eMe=medI/d+JemeJe其中：为作用于地球质心的对偶力旋量；表示地球的对偶惯量矩阵；表示地球的质量；表示地球的转动惯量。OsXsYsZs月球在J2000日心坐标系下的运动学方程为 qm=12 qm wmm(3)qm=qm+12qmrmmqmrmm wmm=wmm+(rmm+wmmrmm)wmm其中：；是月固坐标系相对于日心坐标系的四元数；是月球在日心坐标系下的位置向量在月固坐标系下的表示；速度旋量为，是月球的角速度在月固坐标系下的表示。OsXsYsZs月球在J2000日心坐标系下的动力学方程为第 3 期张锦绣，等：基于对偶四元数的空间引力波探测器稳定构型优化设计269 wmm=M1mFmM1m(wmmMm wmm)(4)Fm=fm+mMm=mmdI/d+JmmdJd其中：为作用于月球质心的对偶力旋量；表示月球的对偶惯量矩阵；表示月球的质量；表示月球的转动惯量。Oixiyizi以航天器的质心为坐标原点，建立本体坐标系，3个坐标轴方向分别与航天器的惯量主轴重合。航天器在日心坐标系下运动学方程表示为15 qi=12 qi wii(5)i=1,2,3 qi=qi+12qiriiqiirii wii=wii+viiwiivii=rii+wiirii其中：分别表示航天器SC1、SC2和SC3；对偶四元数，表示第 个航天器的本体坐标系相对于日心坐标系的四元数，表示航天器的位置矢量在本体坐标系下的表示；对偶速度旋量，表示航天器相对于太阳的角速度在本体坐标系下的表示，表示航天器相对于太阳的速度矢量在本体坐标系下的表示。动力学方程为Mi wii=Fii wiiMi wii(6)Mi=midI/d+JimiJiFiiFii=fis+0fisFii=fie+fim+fis+fids+igfiefimfisfidsig其中：表示对偶惯量矩阵，表示航天器的质量，表示航天器的转动惯量；表示航天器受到的对偶力。当图1中C点为虚拟点时，航天器以太阳为中心天体，若只考虑太阳引力，那么航天器所受到的对偶力为，表示太阳的引力。当图1中C点以地球为中心天体时，航天器受到地球、月球和太阳的引力和太阳光压，那么航天器所受到的对偶力为：，、和分别表示地球引力、月球引力、太阳引力、太阳光压和重力梯度力矩。fie=emiriirie?riirie?3+rie?rie?3fim=mmiriirim?riirim?3+rim?rim?3fis=smirii?rii?3ig=3emi(riirie)Ji(riirie)?riirie?5(7)e=398 600.441 90m=4 902.800 076s=132 712 440 040.944 00fids=PA(r/r3)AU2(1+)Pr其中：地心引力常数为 km3/s2，月心引力常数为 km3/s2；日心引力常数为 km3/s2；太阳光压为，表示太阳光压强，它与航天器到太阳的距离有关，A表示航天器表面积，AU表示太阳相对于地心的位置矢量，表示地球与太阳的距离，表示卫星表面材料的反射系数。1.2 编队成员间相对动力学方程OixiyiziOjxjyjzji,j=1,2,3,i,j本文主要研究航天器之间的相对运动，即两两航天器本体坐标系和（）之间的相对运动。利用乘性误差对偶四元数，该相对运动可表示为 qij=qj qi=qij+12qijriij(8)qi qjijqijijriij=rijriiijOixiyizi其中：和分别表示第 个和第 个航天器相对于日心坐标系的对偶四元数表示；为第 个航天器相对于第个航天器的姿态四元数；为第 个航天器相对于第 个航天器的相对位置矢量在坐标系下的分量。OixiyiziOjxjyjzj坐标系相对于坐标系的运动学方程为2 qij=qij wiij(9)wiij=wiij+(riij+wiijriij)OixiyiziwiijOixiyizi其中：为相对速度旋量在坐标系中的表示，为相对角速度在中的表示。航天器之间相对动力学方程为 wiij=wii qij wjj qij+wiij(qij wjj qij)(10)本研究关注的核心指标是探测器之间的相对运动参数，但由于编队中包含了自然天体，所以需要建立探测器和临近自然天体之间的相对运动关系。航天器相对于临近天体的相对运动方程为2 qiC=qiC wiiC(11)qiC=qC qi=qiC+12qiCriiC qi qCiqiCiriiC=riCriiiOixiyizi wiiC=wiiC+(riiC+wiiCriiC)其中：，和分别表示第 个航天器和天体相对于日心坐标系的对偶四元数表示；为第 个航天器相对于天体的姿态四元数；为第 个航天器相对于天体的相对位置矢量在坐标系下的分量。为相对速度旋量。航天器相对于临近天体的相对动力学方程为 wiiC=wii qiC wCC qiC+wiiC(qiC wCC qiC)(12)2 稳定构型优化设计 2.1 待优化参数d本文将进入预定轨道的起始时刻 作为优化参数之一，寻找使得航天器构型最稳定的时刻。270深空探测学报（中英文）2023年在日心坐标系下，确定了轨道优化的起始时刻后，此时地球、月球的初始位置和速度均可通过星历DE421获得，通过式（1）（4）可得到初始时刻之后地球和月球在日心坐标系下的位置和速度。航天器本体坐标系相对于日心坐标系的四元数表示为q=cos(2),sin(2)n=f(,)(13)nn=nx,ny,nz其中：为单位轴；为绕单位轴转过的角度；为方位角；为极角。单位轴由极坐标形式表示为nx=cos sinny=sin sinnz=cos(14),0,2其中：。rii=r,0,0Tr=p/(1+ecos f)p=a(1e2)aefvii=vr,vt,0Tvr=/a(1e2)esin fvt=/a(1e2)(1+ecos f)p、e、f在航天器轨道坐标系中，位置矢量为，其中，为轨道半长轴，为偏心率，为真近点角。速度矢量为，为径向速度,为横向速度。所以当航天器本体坐标系和质心轨道坐标系重合，每个航天器需要6个参数作为待优化参数，优化得到航天器的轨道。2.2 目标函数的确定CFCF1CF2根据表征等边三角形稳定性要素指标要求，目标函数由卫星之间的臂长变化量和呼吸角变化量两部分线性组成，并求最小值。CF1(t)=1c1wtt0L12(t)L2+L23(t)L2+L13(t)L2)dtCF2(t)=1c2wtt0(1(t)02+2(t)02+3(t)0)2dtCF(t)=w1CF1(t)+w2CF2(t)(15)w1、w2c1、c20=60L其中：为权重；为归一化常数；为卫星之间的臂长。干涉臂长和呼吸角的计算公式为Lij(t)=?rij(t)?i(t)=arccosrij(t)rik(t)?rij(t)?rik(t)(16)i,j,k=1,2,3i,j,krij(t)rik(t)rjk(t)其中：，、为两航天器之间相对位置向量。2.3 遗传算法优化方法根据前面的推导结果可以看出，对于空间引力波CF(t)探测系统，只要确定了构型的优化参数，就可以计算出地球、月球和航天器的位置和速度。取不同的优化参数，就可以得到不同的构型，构型稳定性问题可以抽象为求取目标函数全局最小值的问题。对于引力波探测器轨道的优化，目标函数有19个自变量，存在数量众多的局部极小值。使用简单的局部优化算法易陷入局部极小值点而无法找到满足要求的轨道，所以这里选择遗传算法作为全局优化算法。遗传算法只利用目标函数的信息，无需梯度等高价信息，适用于大规模、高度非线性的不连续的目标函数的优化，具有很强的通用性，是构型优化设计算法的一个理想选择。2.3.1 编码和适应度函数遗传算法的编码方式包括二进制编码、浮点数编码等，本文采用二进制编码方式，通过控制二进制的位数决定参数的精度。选择合适的适应度函数是设计遗传算法过程中的一个关键问题。本文以式（17）作为适应度函数f=11+CF(17)2.3.2 遗传算子1）选择算子通过选择算子选出种群中较优的个体形成一个新的种群，可使得种群中的个体不断逼近最优解。选择算子可采用常用的轮盘算法，设定个体被选中的概率如下Sj=fjnj=1fj(18)fjn其中：为个体 的适应度值；为种群数目。在上述选择方式下，个体的适应度值越大，被选择的概率也越大，则其基因构造被遗传到下一代的可能性越大。这种选择方式使得种群中的个体都存在被选中的概率，保证了较小个体的遗传性。2）交叉算子交叉是指按交叉概率从种群中选择部分个体，通过交换个体部分基因形成新的个体，从而形成新的种群。交叉算子可选择单点交叉算子。3）变异算子变异算子模拟的是生物进化过程中的基因突变现象，具体操作为：令变异概率为Pc，则变异的节点个数为nSPc，S为单个个体的节点数；根据生成的随机数r决定变异方向，选中某些个体上的某些节点，然后用新值替代这些节点上的原值。变异算子维持了群体第 3 期张锦绣，等：基于对偶四元数的空间引力波探测器稳定构型优化设计271的多样性，有利于防止早熟现象的出现。2.3.3 算法流程本文算法的计算步骤如下：1）生成初始种群；2）根据设计的适应度函数，计算该代种群所有个体的适应度；3）判断是否满足终止条件，若不满足则到步骤4，若满足则到步骤5；4）根据交叉概率和变异概率，经过选择算子、交叉算子和变异算子产生下一代种群，然后返回步骤3；5）选出最后一代种群中的最优个体，输出相应的最优解，算法终止。3 仿真算例与结果分析为了验证所提方法的普适性，以地心轨道和日心轨道空间引力波探测任务为应用背景，对构型优化问题进行仿真验证。3.1 地心轨道空间引力波探测系统构型优化仿真以地心轨道空间引力波探测系统为应用背景，对p、c1=3/104c2=3/400w1=0.5w2=0.5P=4.56106Nm2=0.21A=12 m2AU=149 597 870 691所提方法进行仿真验证。3个航天器各自围绕地球飞行，轨道高度为10万km，假设轨道偏心率和真近点角均为0，航天器轨道由4个参数唯一确定，参数范围如表1所示。轨道起始日期为2034年1月1日，在1年的时间内优化得到最优的构型。遗传算法的种群规模为100，交叉概率为0.7，变异概率为0.08，迭代次数为100代。其余参数分别为，。，km。将优化结果表示为位置速度后如表2所示。表 1 地心轨道编队待优化参数及取值范围Table 1 Geocentric orbit formation parameters to be optimizedand their value ranges待优化参数取值范围半通径pi/km99 900100 100旋转角度i/rad02方位角i/rad02极角i/rad02天数d/d0365.25 表 2 地心轨道编队优化结果Table 2 Geocentric orbital formations optimization results航天器x/kmy/kmz/kmvx/（kms1）vy/（kms1）vz/（kms1）SC192 047.119 927 807.151 427 585.698 40.776 7161.138 8331.443 746SC212 013.557 466 493.722 973 783.279 91.976 2680.051 6110.275 268SC379 297.944 730 902.071 652 579.608 91.200 7261.070 8871.181 493 L12L12L13L13L23L23图2是地心轨道编队航天器SC1、SC2和SC3之间的干涉臂长变化率曲线，图2中绿色曲线表示干涉臂长的变化率，最大值为0.14%；红色曲线表示干涉臂长的变化率，最大值为0.14%；蓝色曲线表示干涉臂长的变化率，最大值为0.18%。0501001502002503003504000.200.150.100.0500.050.100.150.20L12L13L23时间/d臂长变换率/%图 2 地心轨道编队航天器之间的干涉臂长变化率Fig.2 Geocentric orbit formation arm length variation between spacecraft112233图3是地心轨道编队航天器SC1、SC2和SC3轨道构型呼吸角的变化曲线，图中绿色曲线表示呼吸角，变化量最大值为0.14；红色曲线表示呼吸角，变化量最大值为0.084；蓝色曲线表示呼吸角，变化量最大值为0.13。59.8559.9059.9560.0060.0560.1060.15呼吸角/（）050100150200250300350400123时间/d 图 3 地心轨道编队呼吸角Fig.3 Geocentric orbit formation breathing angle272深空探测学报（中英文）2023年v1v2v3图4是地心轨道编队航天器SC1、SC2和SC3相对速度变化曲线，图中绿色曲线表示SC1和SC2之间的相对速度，变化量最大值为6.17 m/s；红色曲线表示SC1和SC3之间的相对速度，变化量最大值为4.56 m/s；蓝色曲线表示SC2和SC3之间的相对速度，变化量最大值为6.70 m/s。86420246相对速度/（ms1）050100150200250300350400v1v2v3时间/d 图 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方向参数L12/kmL13/kmL23/km1/（）2/（）3/（）v12/（ms1）v13/（ms1）v23/（ms1）径向均值173 671173 313173 51860.056 260.262 860.174 00.055 70.011 40.037 1标准差4804164300.234 80.278 30.243 11.653 92.287 12.586 0切向均值173 479173 261173 35260.025 960.151 260.078 30.029 90.006 90.022 1标准差2373062800.172 00.137 90.154 61.658 22.286 32.571 3法向均值173 515173 283173 35860.038 460.171 860.081 80.033 20.012 70.024 7标准差2913413070.192 10.168 10.170 91.648 52.286 32.558 3 表 4 速度误差对地心轨道编队构型的影响Table 4 Effect of velocity error on configuration of geocentric orbital formation方向参数L12/kmL13/kmL23/km1/（）2/（）3/（）v12/（ms1）v13/（ms1）v23/（ms1）径向均值173 274173 236173 21260.01360.03559.9990.003 370.008 230.015 4标准差94.69088.485151.1850.039 30.052 60.079 01.5322.2602.327切向均值173 274173 209173 21859.99760.03 460.0020.000 1400.000 280.016 2标准差109.917116.087156.5260.057 50.061 40.082 01.5302.2622.334法向均值173 284173 221173 15860.00460.04 059.9670.004 510.003 820.005 38标准差68.16577.058107.0580.031 40.037 90.047 91.5272.2612.329 3.2 日心轨道空间引力波探测系统构型优化仿真L=2.5106km以日心轨道空间引力波探测系统为应用背景，对所提方法进行验证。该方案中3颗航天器各自围绕太阳飞行，所构成等边三角形的臂长，编队三角形中心与地球相差约20。本算例中待优化参数范围如表5所示，遗传算法参数配置与地心轨道编队算例相同。将优化结果表示为位置速度后如表6所示。图5是日心轨道编队航天器SC1、SC2和SC3之间表 5 日心轨道编队待优化参数及取值范围Table 5 Heliocentric orbit formation parameters to beoptimized and their value ranges待优化参数取值范围单位pi0.9991.001AUi02radi02radi02rade00.01f02radd0365.25d第 3 期张锦绣，等：基于对偶四元数的空间引力波探测器稳定构型优化设计273L12L12L13L13L23L23的干涉臂长变化率曲线，图中绿色曲线表示干涉臂长的变化率，最大值为0.48%；红色曲线表示干涉臂长的变化率，最大值为0.25%；蓝色曲线表示干涉臂长的变化率，最大值为0.32%。0.40.30.20.100.10.20.30.40.5050100150200250300350400L12L13L23时间/d臂长变换率/%图 5 日心轨道编队航天器之间的臂长变化率Fig.5 Heliocentric orbit formation arm length variation between spacecraft112233图6是日心轨道编队航天器SC1、SC2和SC3轨道构型呼吸角的变化曲线，图中绿色曲线表示呼吸角，变化量最大值为0.36；红色曲线表示呼吸角，变化量最大值为0.38；蓝色曲线表示呼吸角，变化量最大值为0.26。59.659.759.859.96060.160.260.360.4050100150200250300350400呼吸角/（）123时间/d 图 6 日心轨道编队航天器呼吸角Fig.6 Heliocentric orbit formation spacecraft breathing angle v1v2图7是日心轨道编队航天器SC1、SC2和SC3相对速度变化曲线，图中绿色曲线表示SC1和SC2之间的相对速度，变化量最大值为1.82 m/s；红色曲线表示SC1和SC3之间的相对速度，变化量最大值为v31.05 m/s；蓝色曲线表示SC2和SC3之间的相对速度，变化量最大值为0.92 m/s。2.01.51.00.500.51.01.52.0相对速度/（ms1）050100150200250300350400v1v2v3时间/d 图 7 日心轨道编队航天器之间的相对速度Fig.7 Heliocentric orbit formation relative velocity between spacecraft 上述结果表明，在考虑了太阳引力作用下，日心轨道空间引力波探测系统在一年内的变化能满足引力波探测的需求。日心编队构型轨道频率变化缓慢，可以在较长时间保持稳定的构型。以优化得到的结果作为初始条件，用Monte-Carlo法分析入轨误差对日心轨道编队构型的影响。分别向航天器的径向、切向和法向位置增加标准差为200 km的位置误差，对构型的影响如表7所示。分别向航天器的径向、切向和法向速度增加标准差为1 cm/s的速度误差，对构型的影响如表8所示。位置误差和速度误差均符合正态分布。可以看出，径向位置误差和切向速度误差导致构型标准差的变化最大。若位置误差和速度误差同时存在，当速度误差为3 cm/s时，位置误差最大不超过90 km，否则构型将会发散。地心轨道编队和日心轨道编队的优化结果对比如表9所示，地心轨道编队的臂长变化率最大值为0.18%，呼吸角变化量最大值为0.14，相对速度的变化量最大值为6.70 m/s；日心轨道编队的臂长变化率最大值为0.48%，呼吸角变化量最大值为0.38，相对速度的变化量最大值为1.82 m/s。地心轨道编队和日心轨道编队的优化结果均满足各自的指标要求，入轨误差中对构型稳定性影响较大的均为径向位置误差和切向速度误差。表 6 日心轨道编队优化结果Table 6 Heliocentric orbital formations optimization results航天器x/kmy/kmz/kmvx/（kms1）vy/（kms1）vz/（kms1）SC1135 509 841.314 859 413 755.739 927 127 784.819 012.868 44424.492 13210.608 043SC2133 979 739.828 260 820 480.983 125 744 263.425 113.077 96324.520 13210.871 339SC3135 019 533.789 058 772 873.331 424 762 996.642 712.858 99624.798 79910.521 874274深空探测学报（中英文）2023年表 7 位置误差对日心轨道编队构型的影响Table 7 Effect of position error on configuration of heliocentric orbital formation方向参数L12/kmL13/kmL23/km1/（）2/（）3/（）v12/（ms1）v13/（ms1）v23/（ms1）径向均值2 506 1872 505 7222 494 52160.17760.21859.9740.4700.4380.756标准差13 57019 29213 0400.4630.4650.3153.1833.5012.465切向均值2 500 8272 500 1422 494 24760.06960.09659.8350.0430.0010.051标准差6 9413 4653 4640.1560.2160.1541.3660.7270.725法向均值2 501 0342 500 3762 494 24760.07360.10159.8410.0420.0060.051标准差6 3893 0572 8400.1480.2050.1401.2680.6690.611 表 8 速度误差对日心轨道编队构型的影响Table 8 Effect of velocity error on configuration of heliocentric orbital formation方向参数L12/kmL13/kmL23/km1/（）2/（）3/（）v12/（ms1）v13/（ms1）v23/（ms1）径向均值2 502 0112 501 0092 499 65160.11060.13659.8790.0430.0060.053标准差8 8865 1284 9540.1870.2550.1801.7831.0531.080切向均值2 509 3112 507 5012 494 46160.21760.31060.0340.4890.5150.940标准差16 76115 07814 4400.4310.4450.3664.1012.8052.821法向均值2 501 2052 500 4462 494 87760.08460.10959.8470.0450.0040.049标准差6 3953 1032 9100.1490.2060.1411.2770.7050.663 表 9 仿真算例数据对比Table 9 Comparision of simuilation results变化量地心轨道编队日心轨道编队优化结果参考指标8优化结果参考指标16臂长变化率/%L120.1410.4810L130.1410.2510L230.1810.3210呼吸角/（）10.140.20.36120.080.20.38130.130.20.261相对速度/（ms1）v16.17101.8210v24.56101.0510v36.70100.9210 4 结论本文以空间引力波探测器编队构型设计为研究背景，重点研究了探测器动力学建模和轨道优化算法。仿真结果表明，所提优化算法可以使引力波探测器在持续一年的时间里，保持构型的稳定性，其中，无中心天体的日心轨道探测器编队的臂长变化率保持在0.48%以内，呼吸角最大变化量保持在0.38%以内；地心轨道探测器编队的相对速度变化量保持在6.70 m/s以内。另外，径向位置误差和切向速度误差对编队构型稳定性的影响较大。上述结果说明，所提基于对偶四元数编队动力学模型，可以对探测器有无中心天体的情况进行统一描述，在保证动力学模型精度的同时提升模型的适用性。本文的研究成果能够为未来开展空间引力波探测提供建议与参考。参考文献STEBBINS R T.Rightsizing LISAJ.Classical and Quantum Gravity，2009，26（9）：094014.1SWEETSER T H.An end-to-end trajectory description of the LISAmissionJ.Classical and Quantum Gravity 2005，22（10）：S429.2LUO J，CHEN L S，DUAN H Z，et al.TianQin：a space-bornegravitational wave detectorJ.Class Quantum Gravity 2016;33（3）：19.3HU W R，WU Y L.The Taiji Program in Space for gravitational wavephysics and the nature of gravityJ.National Science Review,2017,4（5）：685-686.4XIA Y，LI G，HEINZEL G，RUEDIGER A，et al.Orbit design for theLaser Interferometer Space Antenna（LISA）J.Science ChinaPhysics，Mechanics and Astronomy，2010，53（1）：179-186.5CHIHANG Y，HAO Z.Formation flight design for a LISA-likegravitational wave observatory via Cascade optimizationJ.Astrodynamics 2019，3（2）：155-171.6万小波，张晓敏，黎明.天琴计划轨道构型长期漂移特性分析J.中国空间科学技术 2017;37（3）：110-116.WAN X B，ZHANG X M，LI M.Analysis of long-