空间引力波探测器轨道确定及仿真分析

李康康，谷德峰，安子聪，宋佳凝，朱炬波

（1.中山大学 天琴中心，珠海 519082；2.中山大学 物理与天文学院，珠海 519082）

引 言

2015年，地面引力波探测项目激光干涉引力波天文台（Laser Interferometer Gravitational-wave Observatory，LIGO）首次探测到双黑洞合并产生的引力波[1]，开启了观察黑洞合并的新窗口。受地面震动、引力梯度噪声以及干涉臂长受限的影响，地面引力波探测无法覆盖中低频段的引力波。空间引力波探测任务由于采用大尺度空间探测器编队，不受臂长尺度约束，探测频段更广，波源更加丰富[2]，为认识宇宙演化提供一种新的方式，是当前国际热门研究领域。其中LISA（Laser Interferometer Space Antenna）由美国国家航空航天局（National Aeronautics and Space Administration，NASA）和欧洲航天局（European Space Agency，ESA）共同提出[3]，是目前发展最为成熟、最具有代表性的空间引力波探测计划，预计于2034年发射。LISA由3颗探测器构成的大尺度分布式日心轨道编队，位于黄道内，编队中心约在地球与太阳构成的后掠角20°位置上，距离地球5～6.5千万km，由约为250万 km的干涉臂长构成等边三角形[4]。图1给出了LISA轨道在太阳系下的示意图。

图1 类LISA轨道结构示意图[5]Fig.1 Schematic of LISA-type orbit configuration[5]

LISA在发射入轨阶段和科学实验阶段对定轨精度以及臂长的测量精度有一定的要求[6]。一是来自星间激光干涉仪对波束指向的需求；二是来自空间引力波探测系统时间延迟干涉（Time Delay Interferometry，TDI）[7]压制激光频率噪声的需求[8]；三是来自编队构形高精度初始化和稳定性保持的需求[9]。

测控距离远、跟踪弧段少、精度要求高给LISA精密定轨（Precise Orbit Determination，POD）带来一定的挑战。本文以LISA为例，分析研究其探测器的精密定轨问题，为空间引力波探测计划的工程实践提供理论支撑和技术参考。LISA探测器距离地球较为遥远，难以采用低轨卫星的定轨、定位手段。而基于深空网的无线电测距测速定轨技术已成功应用于深空探测任务，如美国火星探测任务“好奇号”（Curiosity）[10]、“嫦娥三号”（Chang'E-3，CE-3）[11]、“嫦娥五号飞行试验器”（Chang'E-5T1，CE-5T1）[12]等。因此，提出采用DSN（Deep Space Network）和CSDN（Chinese Deep Space Network）对LISA探测器进行精密定轨。

尽管LISA早在20世纪90年代提出，但对LISA轨道确定的相关研究仍较少，其中文献[13]对早期的LISA轨道进行了初步定轨仿真，分析了基于DSN、甚长基线干涉测量（Very Long Baseline Interferometry，VLBI）、星间激光测量数据3种不同的组合定轨策略的轨道精度将地球当作质点，且假设测站位于地球中心，未考虑地球的自转，也没有模拟分析不同深空网、跟踪弧长、测量数据类型等因素的影响。为此，本文基于多深空网测量模式下的LISA精密定轨分析，考虑地球自转以及测站的空间位置，模拟生成不同测站的测距、测速数据，采用非线性加权最小二乘（Weighted Least Squares，WLS）和蒙特卡洛方法（Monte Carlo，MC）分析深空网的跟踪弧长、测量数据类型以及不同深空网测站数量及其分布对LISA精密轨道确定的影响。

1 场景环境建立

LISA精密定轨仿真分析中采用DSN和CDSN作为跟踪测站。其中，DSN的3个测站分别位于加利福尼亚戈尔德斯顿、澳大利亚堪培拉和西班牙马德里[14]。CDSN主要包括佳木斯深空测站、喀什深空测站、阿根廷萨帕拉市深空测站[15]。表1给出了上述6个测站的经纬度。

表1 测站的经纬度数据Table 1 Longitude and latitude parameters of stations

模拟生成的观测数据为测距和测速信息，其中观测数据包含了系统误差和随机误差。考虑测距系统偏差主要由测距设备延迟、对流层带来的天顶和视向偏差、电离层带来的视向偏差以及时钟钟差和钟漂等因素组成；产生随机误差的主要是时标和频率源不稳定以及热噪声等因素，并假设其服从高斯分布。仿真分析中，根据文献[14]美国深空网的测站测量精度，本次仿真设置测距系统误差为2.055 m，其中测站设备时延2.0 m，对流层和电离层以及钟差钟漂等带来的误差0.055 m；测距随机噪声的标准差为0.6 m。测速误差仿真时只涉及到随机误差，观测噪声的标准差为0.03 mm/s。

采用文献[4]中给出的最新的、250万 km臂长的LISA设计轨道，轨道初始状态采用文献[5]中初始时间为2028年3月22日12:00:00 UTC，在J2000地心惯性系下（J2000 Earth-Centered Inertial，J2000 ECI）的轨道初值，如表2所示。仿真分析过程：①根据表2给出的LISA探测器初始状态，利用Adams-Cowell轨道积分器分别计算LISA 3个探测器的轨道信息，获得相应的轨道文件；②根据表1所示的测站信息，结合步骤①中得到的轨道文件，利用测距、测速的数学模型迭代计算探测器相对于测站的距离和速度的理论值，并在理论值的基础上添加设定的系统偏差和随机噪声，从而得到观测数据；③利用观测数据对LISA探测器轨道求解，分析不同仿真场景下LISA探测器的定轨精度，其中轨道精度由设计轨道与求解结果之间的偏差衡量。本文采用NUDTTK软件[16]进行仿真分析。

表2 J2000地心惯性系下的LISA探测器初始状态[5]Table 2 Initial state of LISA detector in J2000 ECI[5]

2 测量模型与仿真策略

本文采用深空网的双向测量模式，测距测量模型可表示为

测速信息可由测距模型的差分得到，具体的测速测量模型为

考虑系统误差和随机噪声等，定轨过程的观测方程最终可表示为

LISA可忽略地球和月球重力场[17]，在仿真观测数据以及定轨过程都不考虑其影响，因此既不会对仿真结果造成重大影响又可简化仿真过程以提升计算效率。另外，考虑卫星入轨会有位置和速度的误差，参考文献[4]，对表2中的设计轨道初值添加对应的初始误差后作为定轨的初始位置速度，即分别在初始位置的x，y，z方向上各添加10 km误差，初始速度的3轴方向各添加1 cm/s误差作为定轨的初始状态。采用非线性加权最小二乘算法估计轨道初始状态，根据测距和测速的观测噪声对测量数据进行加权定轨，并采用MC方法进行35次仿真，计算定轨位置和速度的均方根（Root Mean Square，RMS）误差的平均值评估定轨精度。定轨模型和参数设置总结于表3。

表3 仿真参数与定轨模型Table 3 Simulation parameters and orbit determination model

3 仿真与结果分析

本节从定轨弧长、测轨数据类型、测站分布及其数量等方面仿真分析LISA定轨的精度。

3.1 跟踪弧长对定轨精度的影响

由于观测数据中存在随机误差，一般随着定轨弧长增加，定轨精度会因逐渐平滑随机误差以及测站对探测器观测几何约束逐渐变强而提高。下面通过分析不同的定轨弧长对LISA轨道确定的影响，讨论上述仿真条件下LISA轨道的位置、速度估计精度。以DSN作为跟踪测站，定轨弧长分别取[2，3，5，7，10，15，20，25，30]d，采用测距/测速组合测量数据进行初始轨道的解算。

图2为多次MC仿真统计的LISA探测器的位置、速度RMS误差的平均值曲线，其中误差棒对应于统计的标准差。从图2中可以看出，定轨精度随着定轨弧长增加而提高，标准差随着定轨弧长增加而降低，可信度提高。LISA系统的3个探测器在同一弧长定轨精度不完全相同，这主要是因为即使采用同一定轨弧长，但由于探测器和测站的几何关系以及测量数据不完全相同、观测的数据量和随机误差都会有所不同等因素导致的，但定轨精度的量级是基本一致的。

图2 LISA探测器定轨误差（1σ）Fig.2 Orbit determination errors of LISA detector（1σ）

当定轨弧长达到20 d时，LISA 3星的位置误差均可达92 m，速度误差均可达3.3 mm/s；此后随着定轨弧长的增加，定轨精度变化不大。

表4中给出了LISA-1探测器在轨道径向（R）、切向（T）、法向（N）以及三维方向（3D）下的定轨位置精度。从表4中可以看出，径向定轨精度受定轨弧长的变化影响不大，在2.05 m附近波动，表现为测距数据的系统偏差；而切向和法向的位置精度会随着定轨弧长的变化逐渐提高，且法向的位置精度最差。这是由于测距信息提供的测站与探测器的几何距离主要体现在径向上，对切向和法向的约束较弱，需增加跟踪弧段、增加观测数据量来提高定位精度。

表4 不同跟踪弧长下LISA-1定轨位置误差（1σ）Table 4 LISA-1 position errors of orbit determination with different tracking arc length（1σ）

3.2 观测数据类型对定轨精度的影响

图3和图4为LISA-1探测器在不同观测数据类型（测距/测速、单独测距、单独测速）下的定轨位置误差。从图3和图4中可知：用测距、测速数据定轨的收敛情况与使用测距/测速数据定轨的收敛情况一致，根据图4中的定轨位置曲线可以归纳出仅使用测速数据定轨的精度较低，定轨位置误差大。由上述分析可知，定轨弧长达到20 d后，定轨精度随定轨弧长增加提升不再明显，故采用DSN 20 d的跟踪弧长分析测距、测速以及测速/测距融合3种测量方案的LISA定轨精度，以此讨论测距、测速数据类型对定轨精度的影响程度。

图3 LISA-1探测器在测距/测速、测距数据类型下的定轨误差（1σ）Fig.3 Orbit determination error of LISA-1 detector under ranging/velocity measurement and ranging data types（1σ）

图4 LISA-1探测器在测速条件下的定轨误差（1σ）Fig.4 Orbit determination error of LISA-1 detector under the condition of velocity measurement（1σ）

图5是利用定轨弧长为20 d并分别采用测距/测速、测距、测速数据定轨后的位置速度精度。从图5中易得知：同种定轨测量数据类型下，LISA 3星的定轨精度量级一致；采用测距/测速组合定轨比仅测距数据定轨精度有所提高；仅采用测速数据的定轨精度最差，位置误差大于10 km，速度误差大于1 m/s，无法满足高精度定轨的需求。

图5 不同测量数据类型下的定轨误差（1σ）Fig.5 Orbit determination errors with different types of measurement data（1σ）

以LISA-1为例，表5给出了20 d跟踪弧长下，不同定轨数据类型在径向、切向、法向和三维方向上的位置精度。由表5中数据可知，测距和测距/测速方案定轨中，径向误差约2.0 m，这与表4中的径向结果相吻合；而仅测速数据定轨方案的径向定位误差较大。对数据进一步分析可得：在20 d的定轨弧长下，采用测距/测速定轨精度比仅采用测距数据定轨位置精度提高32.23%，比仅采用测速定轨位置精度提高99.52%。可见，测距数据对定轨径向精度的贡献比较大，仅利用测速数据，径向定轨精度明显较差；联合测速数据，可适当提高飞行方向和法向的定轨精度。

表5 不同测量数据类型LISA-1的定轨结果（1σ）Table 5 Orbit determination results of LISA-1 with different types of measurement data（1σ）

3.3 多测控网组合方案对定轨精度的影响

除跟踪弧长和定轨数据类型外，不同的深空网因测站数量和位置分布不同也会影响定轨效果，下面分析DSN、CDSN以及DSN/CDSN 3种不同深空网定轨方案对LISA定轨精度的影响。同时，考虑到LISA系列探测器均分布在同一太阳轨道平面上，多数仅真近点角（true anomaly）不同，深空网对LISA 3星的覆盖情况基本一致，故本节仅以LISA-1为例进行分析。

3.3.1 深空网对探测器的跟踪观测情况

图6为DSN和CDSN跟踪7 d高度角的变化示意图，从图6中可看出，高度角最大不超过80°，且高度角的变化周期约为3 d，主要由LISA轨道的特点决定，因其位于地球后的太阳轨道，与地球的相对位置几乎保持不变，测站随着地球自转与探测器发生周期性变化。若设最小观测高度角为10°，则DSN/CDSN和DSN在跟踪20 d弧长下对LISA覆盖率分别达到70.24%和59.62%，而由于CDSN的测站分布不如DSN均匀，对LISA探测器的覆盖率为45.07%，略低于DSN。

图6 不同深空网对LISA-1的跟踪高度角变化Fig.6 Elevation angles of LISA-1 tracked by DSN and CDSN

表6给出了采用不同深空网进行跟踪20 d时，不同测站数量对LISA-1的观测情况。由于测站分布在不同半球，DSN和CDSN单独跟踪时，都不存在3个测站同时观测到探测器的情况，且单站观测率最高。而DSN/CDSN联合跟踪测量时，双站和3个测站共视率共占比44.87%，其中3站同时观测率为14.57%。

表6 测站共视率统计（LISA-1）Table 6 Observation rate statistics of simultaneous observation at the ground stations（LISA-1）

图7分别展示了DSN和CDSN对LISA-1在2028年3月22日12:00:00UTC —2028年3月29日12:00:00UTC期间跟踪7 d的测站跟踪时段情况。DSN和CDSN联合观测利于增加有效观测弧段，利于提高观测率。

图7 不同深空站跟踪弧段统计（LISA-1）Fig.7 Tracking arcs analysis of LISA-1 by different deep space stations

3.3.2 多网深空网定轨精度分析

根据上述对深空网的观测几何分析可知，不同的深空网观测的弧段、连续观测弧长以及测站几何约束能力各不相同，这都会对LISA的定轨精度产生影响。

图8中分别给出了基于DSN、CDSN以及DSN/CDSN 3种不同深空网下的LISA-1轨道确定的结果，其中测量信息采用测距/测速组合模式。图8中可以看出：在相同跟踪弧长的情况下，采用DSN/CDSN组合网络比仅采用DSN和CDSN的单深空网络的定轨精度明显提高。分析2 d的定轨弧长发现：CDSN定轨精度明显低于DSN，这是由于CDSN具有较大数据量，同时注意到，尽管CDSN的总观测覆盖率比DSN低，但CDSN双站共视率略高，此后基于CDSN的定轨精度会在一定程度上略高于基于DSN的定轨精度，但总体上两者的定轨精度均在相同的量级上变化；当定轨弧长达到7 d后，两者的定轨结果差异明显减小。

图8 在不同深空网下的初始定轨误差（1σ）Fig.8 Orbit determination errors of LISA-1 using different deep space nets（1σ）

联合DSN/CDSN，尤其在跟踪弧段较短的时候，可提高收敛速度，从图8中可看出，当定轨弧长在7 d内，多深空网对定轨精度提升效果较明显。另外，统计分析不同定轨弧长[2，3，5，7，10，15，20，25，30]d时的定轨结果，可得基于DSN/CDSN组合网络的LISA-1轨道的位置平均误差为157.312 2 m，速度平均误差为6.0 mm/s；与相同情况下基于DSN的定轨结果相比较，平均位置估计精度提高43.73%，平均速度估计精度提高43.40%。

4 结 论

本文对基于深空网的空间引力波探测器——LISA探测器定轨进行了仿真分析，采用无线电测距、测速方法对LISA进行跟踪观测，利用非线性加权最小二乘批处理、蒙特卡洛方法分析了跟踪弧长、测轨数据类型、测站数量和分布对LISA探测器精密定轨的影响。当定轨弧长达到20 d后，定轨精度的提高效果随跟踪弧长的增加变得缓慢。随后采用不同观测类型数据进行分析组合定轨策略：测距信息对定轨影响较大，且对轨道径向有着较强的约束；仅测速数据无法实现精密定轨；在DSN测控体系中，对深空探测器精密定轨以测速为主，只有非常少量的测距数据，测速辅以少量的测距数据，可以实现测距和测速的优势互补。最后，分析了DSN和CDSN对LISA-1的观测几何约束情况，DSN/CDSN联合观测可有效增加多站共视比率；统计[2，3，5，7，10，15，20，25，30]d定轨弧长精密定轨后的平均精度，DSN/CDSN位置平均误差为157.312 2 m，较DSN定轨的位置估计精度提高了43.73%。

本文的仿真方法和结果分析对天琴计划[18]、太极计划[19]等其它空间引力波探测计划探测器的精密轨道确定具有一定的参考意义。