13 m天线观测深空探测器数据的处理与分析

贺庆宝，陈祎豪，杨轩，鄢建国，李斐

(1.南方科技大学 地球与空间科学系，深圳518055；2.武汉大学 测绘遥感信息工程国家重点实验室，武汉430079)

1 引言

武汉大学13 m天线于2013年建成，主要应用于地球遥感卫星数据的接收，同时也肩负深空探测学科的建设、教学和科研工作[1]。为应用于深空探测，天线陆续配备了X波段制冷接收机、氢原子钟、数据采集系统，同时也对深空探测器进行了观测，如我国的月球探测器嫦娥三号(Chang’E-3,CE3)、欧空局及美国的火星探测器火星快车(Mars Express,MEX)和火星轨道勘察器(Mars Reconnaissance Orbiter,MRO)、美国的木星探测器Juno等。在观测时间安排上，我们在安装制冷接收机前后以及安装氢原子钟前后都观测了嫦娥三号，用于对比分析接收信号信噪比(signal to noise ratio,SNR)的强弱和接收频率精度的高低。同时我们在火星与太阳角距离较小和较大时都对火星探测器(MROMEX)进行了观测，用以分析等离子体对信号相位的影响。

接收信号信噪比的强弱直接决定信号能否被成功提取，同时也会影响后续提取频率与相位的精度。提取得到的多普勒频率可用来对探测器进行测定轨，也可用来做行星科学的研究，如重力场提取与行星大气的研究[2,3]；提取得到的相位抖动可以用来研究空间介质对电磁波影响的特性与规律[4,5]。随着我国深空探测项目的开展，如后续的探月项目嫦娥五号(Chang’E-5,CE5)、火星探测、小行星探测等，13 m天线不仅可用来做深空探测学科建设，将来也有可能应用在深空探测任务中。本文通过分析13 m天线近两年来观测深空探测器的数据，及介绍天线接收功率以及多普勒频率精度、相位抖动情况，为相关学科建设以及加入深空探测任务、空间科学研究等做参考与准备。

本文第2章主要介绍用13 m天线观测不同深空探测器的功率谱及制冷接收机安装前后的信噪比变化。第3章主要介绍多普勒频率的精度及稳定性。第4章主要介绍一种事后计算总相位的方法及相位抖动分析。第5章对全文做了总结。全文时间均为协调世界时(coordinated universal time,UTC)。

2 接收功率分析

目前深空探测任务的下行频率大部分使用的都是X波段频率，13 m天线配备X波段接收机能对大部分深空探测器进行观测。在深空探测中，探测器距离地球小于2×106km时，X波段下行频率限制在8 450～8 500 MHz；当探测器距离地球大于2×106km时，X波段下行频率限制在8 400～8 450 MHz[6]。我们以接收频率8 470 MHz为例，计算了13 m天线归一化的增益曲线，计算方法见参考文献[7]。图1a)展示了13 m天线的实物图，图1b)给出了天线在8 470 MHz的波束图。由图1b)中数据可知，13 m天线X波段的主波束宽度(两边各下降3 dB)为0.16°。

13 m天线X波段接收机的制冷装置在2017年12月2日进行安装，我们在安装前后分别对CE3进行了观测。图2a)与图2b)分别是2017年11月10日和2018年1月9日观测CE3得到载波通道信噪比的结果。在计算信噪比时，采用1 024点做快速傅里叶变换(fast Fourier transform,FFT)，然后1 s积分，具体方法见参考文献[7]。由图2结果可知，在接收机制冷装置安装后，信噪比提高了约3 dB。考虑CE3发射功率不变，信噪比提高3 dB意味着天线系统噪声温度降低了一半，说明制冷效果良好。

在制冷接收机安装之后，我们利用13 m天线对火星探测器MRO,MEX以及木星探测器Juno也进行了观测。火星离地球最近时约5.5×107km，13 m天线X波段主波束在此距离上对应的平面宽度约为1.5×105km，远大于目前火星探测器的轨道高度[8]。因此在观测火星探测器时，只需将天线对准火星即可。另外，我们在2018年1月10日也对木星探测器Juno进行了观测，此时木星距离地球约8.7×108km。由于Juno飞行轨道距离木星很远[9]，观测时我们根据轨道预报文件计算Juno的位置，然后将天线对准Juno进行观测。

图1 13 m天线的实物图以及X波段波束图

图2 制冷接收机安装前后CE3信号的信噪比

图3给出了13 m天线观测各个探测器信号得到的自功率谱。自功率谱是通过利用1 s数据做FFT然后平方所得，频率分辨率为1 Hz。图中的Pcal信号为相位校准信号。图3a)是2018年1月9日观测CE3得到的自功率谱；图3b)是2018年1月10日观测Juno得到的自功率谱；图3c)是2018年1月11日观测MRO得到的自功率谱，此时火星距离地球约2.8×108km；图3d)是2018年3月5日观测MEX得到的自功率谱，此时火星距离地球约2×108km。从图3可见，CE3信号最强，其次是MRO和MEX，最弱的是Juno，这主要是因为接收信号的功率与距离的平方成反比。尽管如此，我们还是成功接收并提取到了远在8.7×108km外的Juno探测器信号。

图3 13 m天线观测月球、火星、木星探测器信号的自功率谱

3 频率分析

3.1 接收频率精度分析

接收频率扣除探测器的发射频率即为多普勒频移，它可用来对探测器进行测定轨及开展行星科学研究；因此，接收频率精度的高低会影响测定轨及科学研究的结果。接收频率的精度主要受频率发射端与接收端基准频率的稳定性、随机噪声、空间介质等影响。为提高13 m天线基准频率的稳定性，13 m天线于2017年12月19日安装了氢钟，替代了之前所用的铷钟。理论上氢钟能提供更稳定的基准频率。图4给出了氢钟安装前后接收频率的精度对比，积分时间均为1 s。图4a)与图4b)分别是2017年12月10日和2018年1月9日观测CE3计算得到的接收频率。图中显示的是转换成基带之后的频率，真实的接收频率等于基带频率加上本振频率8 469 MHz。频率计算方法见参考文献[10]。图4c)是2017年12月10日CE3接收频率的9次多项式拟合残差，其均方根(root mean square,RMS)为42 mHz。不仅如此，2017年11月10日观测CE3的接收频率，其精度也约为42 mHz。可见在氢钟安装之前，13 m天线观测CE3的接收频率所能达到的精度即为42 mHz。图4d)是2018年1月9日即氢钟安装之后观测CE3的接收频率9次多项式拟合残差，其RMS为20 mHz。考虑CE3的发射信号都是由地面上行站提供频率基准，图4d)中接收频率精度提高一倍的主要原因则来自于13 m天线氢钟的安装。

图4 氢钟安装前后13 m天线观测CE3得到的接收频率

如第2章所提到，13 m天线除了观测CE3，还多次观测了火星探测器MRO和MEX，以及木星探测器Juno。我们对这些观测数据都进行了处理，计算了其接收频率，并分析其精度。图5给出了氢钟安装后不同探测器接收频率的拟合残差，积分时间均为1 s。图5a)给出了2018年1月10日观测Juno的频率拟合残差，RMS为0.067 Hz。当天观测时段Juno距离木星约7.93×106km。图中结果有一些波动，这可能是由空间介质或者信号穿过小行星带引起。图5b)给出了2018年1月11日观测MRO的频率拟合残差，RMS为0.389 Hz。当天观测时段MRO距离火星从约260 km变化至290 km。图中结果有明显的大幅度变化，周期约数百秒，这应主要由MRO当时离火星较近，受火星重力场影响较大而引起。图5c)给出了2018年3月5日观测MEX的频率拟合残差，RMS为0.021 Hz。当天观测时段MEX距离火星从约8 400 km变化到6 000 km。图中结果整体非常平稳，这应主要由此时飞行器轨道较高，对火星重力场高阶变化不敏感所致。图5d)给出了2019年9月28日观测MRO的频率拟合残差，RMS为0.024 Hz。当天观测时段MRO距离火星约260 km。图中结果没有类似图5b)中大幅度抖动，原因是此数据的时间跨度较短，重力场引起的长周期抖动已被高阶拟合去除。图中有很小幅度的波动，这可能是受空间等离子体的影响而引起，当天火星与太阳夹角为8.6°。

图5 13 m天线观测不同深空探测器的频率拟合残差

综上所述，13 m天线观测CE3以及火星探测器的接收频率精度约为20 mHz，且能以较高精度提取木星探测器的接收频率。

3.2 频率稳定性分析

前面给出了1 s积分的接收频率结果。除此之外，我们还计算了5 ms积分的接收频率，然后分别利用5 ms积分和1 s积分的结果做阿伦方差以分析其频率稳定性。

图6a)给出了2018年1月9日13 m天线观测CE3得到5 ms积分的接收频率拟合残差，其变化范围在±10 Hz左右，这主要是由短时间积分随机误差较大引起。图6b)是图6a)的部分放大图，为了看清其变化的周期，我们把横轴转换成了秒。由图6b)可知，其频率拟合残差存在明显的毫秒量级周期抖动。我们将图6a)结果做阿伦方差，得到结果如图6c)黑色圈线所示，红色点线为直线拟合结果(斜率为-1，代表随机噪声的理论变化趋势)。图6d)是阿伦方差的直线拟合残差，其结果显示时间间隔在毫秒量级时抖动较大，但随着时间间隔越长，抖动幅度越小，最后在时间间隔为1 s左右稳定下来。这种现象不仅存在于CE3观测数据中，也存在于火星和木星探测器的观测数据中。

鉴于此抖动普遍存在于观测数据中，说明其原因来自于13 m观测台站，而最有可能产生此抖动的环节则是由氢钟提供的混频信号。经与氢钟设计方交流后得知，氢钟毫秒量级频率稳定度取决于氢钟锁相接收机所采用恒温晶振的等级，武汉大学13 m天线的氢钟所采用的是常规瑞士oscilloquartz的8789晶振，其毫秒量级稳定度较差，从而造成了抖动现象。

另外，我们利用相同观测数据得到1 s积分的频率结果做阿伦方差，结果如图7所示。由于观测数据长度有限，我们做阿伦方差所使用的时间间隔只设置从1 s到100 s。由图7可知，13 m天线观测CE3的接收频率1 s稳定度约为5.6×10-12。这与利用上海佘山25 m天线进行CE3天地对接得到的频率稳定度(7.0×10-12)基本相当[11]。

图6 2018年1月9日13 m天线观测CE3得到5 ms积分接收频率的拟合残差以及阿伦方差

图7 13 m天线观测CE3得到1 s积分接收频率的阿伦方差

综上所述，13 m天线的氢钟在毫秒量级的频率稳定度较差，但1 s稳定度较高(10-12量级)。目前深空探测的数据处理基本使用1 s甚至更长时间积分以提高SNR，因此13 m天线氢钟的频率稳定度能够满足目前深空探测的需求。

4 总相位计算与分析

深空探测器发射的电磁波信号先穿过空间介质(包括空间等离子体、地球大气电离层)，然后到达地面。与频率相比，电磁波总相位对空间介质的影响更敏感。接收信号总相位的计算一般与频率计算同时进行，如利用锁相环计算频率和总相位。但此方法有其局限性，如锁相环在信噪比较低时容易失锁。由于13 m天线接收深空探测器信号的信噪比一般较低，为此我们提出一种事后计算总相位的方法，能适用于极低信噪比情况。我们用此方法分别计算了CE3以及火星探测器观测数据的总相位并做了分析。

4.1 总相位计算方法

此方法的总体思路为：先计算得出1 s积分的接收频率并做多项式拟合；然后拿拟合频率构建混频函数并与接收信号混频，得到残余相位；最后把残余相位与混频相位相加，得到总相位。下面我们将做具体说明。

(1)计算频率并做拟合

我们先利用参考文献[10]提出的方法计算得到当天接收信号1 s积分的频率，然后再回到原始数据，并把原始数据分为若干段，每段数据的时间跨度为5 s，后面每次计算5 s数据的总相位。每一段5 s数据的信号表达式可表示为：

φtotal表示为：

其中，f0表示这段数据起始点的频率，t为时间，变化范围为0～5 s，θ为初始相位。由于有多普勒效应且时间跨度为5 s，我们设定其频率满足二次多项式变化，a和b分别是二次和一次多项式的系数。如此，将式(2)对时间t求导再除以2π，则可得到理论频率的表达式：

此时我们在之前已计算得到的频率中找到相同时间段的结果，并做二次多项式拟合。拟合后的频率可表示为：

其中，p1,p2,p3为二次拟合后得到的多项式系数，则有p1≈3a,p2≈2b,p3≈f0。

(2)混频并平均

此时我们利用前面频率拟合结果来构建两路正交的混频函数，并与此5 s原始数据做混频处理。混频运算的目的是为了将频率搬移至0 Hz附近，以便最后提取残余相位。混频之后，我们得到两路信号，即I路与Q路。混频运算公式如下：

其中，两路混频信号的相位表示为：

混频信号相位的参数设计主要由式(3),(4)得来，目的是尽可能地逼近真实相位。

混频运算不仅把原始数据的频率搬移至低频，同时也把频率搬移到了高频。为了消除其高频频率，我们将混频后得到的I路与Q路信号分别进行平均处理即积分。采用多长时间积分可根据信号的强弱来决定，对于极弱信号，可使用1 s甚至更长时间积分，以保证相位能够被正确提取，这也是此方法的优势所在。平均之后的结果可表示为：

其中，

此时，我们得到了两路超低频信号，两路信号的相位约等于原始信号的初始相位，即Δφ≈θ。

(3)得到残余相位和总相位

在混频并平均处理后，我们将两路信号进行反正切运算，得到残余相位：

之后把残余相位与混频相位叠加，即得到每一段数据(5 s)的总相位。最后再把每段数据的总相位通过调整2π模糊度以连接起来，即可得到整段观测数据的总相位。

4.2 总相位计算结果与分析

我们利用上述方法计算了13 m天线于2018年1月9日观测CE3的总相位。图8a)和b)分别给出了计算第一段5 s原始数据时得到的混频相位以及残余相位结果，其中残余相位结果使用的是5 ms积分。残余相位值理论上约等于信号的初始相位，短时间内应稳定不变，但图8b)中结果存在明显毫秒量级的周期抖动，抖动幅度约为1 rad，这主要是由13 m天线氢钟提供的本振频率在毫秒量级的稳定性比较差造成。图8a)与b)结果之和即为总相位。

我们把每一段5 s数据的总相位连接起来，并做9次多项式拟合，得到拟合残差，其结果如图9a)所示。图中整体的变化范围在±1.5 rad内，主要原因如前所述，为氢钟引起。另外我们还根据同样的计算步骤得到了1 s积分的总相位结果，并做9次多项式拟合，拟合残差如图9b)所示。图9b)中结果并无类似图8b)中的大幅度抖动，其变化范围在±0.5 rad内，此量级恰好反映地球大气电离层的影响。由此也说明，尽管13 m天线氢钟在毫秒量级稳定性差，但总相位1 s积分的结果基本不受此影响。

图8 2018年1月9日13 m观测CE3数据计算总相位的中间值

图9 2018年1月9日13 m天线观测CE3的总相位拟合残差

除了计算CE3数据的总相位之外，我们还计算了13 m天线观测火星探测器(MROMEX)信号的总相位。火星距离地球遥远，其绕飞的探测器发送回来的电磁波信号需穿过长距离的空间等离子体区域，而该过程会引起电磁波相位的波动，特别是在火星与太阳角距离较小时。图10是从地球视角测量火星与太阳的角距离示意图，图中的β即为两者的角距离。

我们计算了13 m天线于2018年3月5日观测火星探测器MEX的总相位，积分时间为1 s，其拟合残差结果如图11a)中蓝色点线所示，RMS为0.18 rad。当天火星与太阳的夹角为81.8°，信号的传播路径距太阳较远，受等离子体的影响相对较小。另外我们还计算了2019年9月28日13 m天线观测火星探测器MRO的总相位，积分时间为1 s。当天火星距离地球约4×108km，且接收机的制冷设备因检修而暂停工作，造成接收信号的信噪比非常低。我们利用前面提出的方法依然计算得到了其总相位，它的拟合残差如图11a)中红色点线所示，RMS为1.73 rad，远大于图中蓝色点线的波动幅度。这主要是因为当天火星与太阳夹角较小，约为8.6°，信号传播路径靠近太阳，穿过较高浓度的等离子体区域，使得电磁波相位产生较大波动。

图10 火星与太阳的角距离关系图

图11 火星探测器的总相位拟合残差及功率谱

图11a)中没有选用两次MRO的观测数据做比较，而是选用了MEX和MRO的数据，是因为之前几次观测MRO结果均显示其受重力场变化的影响较大，如图5b)所示，不利于用来分析等离子体的影响。2018年3月5日观测MEX结果和2019年9月28日观测MRO结果均显示其受重力场变化的影响极小，如图5c)和d)所示，且考虑MEX及MRO下行频率的标准都由地面台站氢钟提供，因此它们总相位的抖动主要由空间介质引起，变化幅度的不同则主要与火星和太阳夹角有关。

我们把拟合残差做了功率谱分析，结果如图11b)所示。图中在x轴频率小于10-2.4Hz区域主要由测频精度引起，频率大于10-1Hz区域主要由装置随机噪声引起，频率在10-2.4～10-1之间区域主要由空间介质引起。我们把x轴频率在10-2.4～10-1之间的数据做直线拟合，得到红色点线和蓝色点线拟合直线的斜率分别为-2.2和-2.3。Woo等人[12]利用S和X双频的相位数据研究等离子体时获得直线拟合斜率为-2.45，这与图11b)的结果较接近，可见利用13 m天线的观测数据也可用来做等离子体的相关研究。

5 总结

文中分析了近两年来13 m天线观测深空探测器的数据，得到接收信号的信噪比在制冷接收机安装后提高了3 dB，其接收频率精度在氢钟安装之后提高了约1倍，达到0.02 Hz。分析结果显示13 m天线氢钟在毫秒量级的稳定度较差，但1 s的稳定度较高(10-12量级)，能够满足目前深空探测的需求。文中还提出一种事后计算总相位的方法，此方法能够适应于极低信噪比情况。在火星与太阳较近时，其火星探测器观测数据的总相位受等离子体影响引起较大幅度抖动。整体分析结果表明，13 m天线除了可提取高精度多普勒数据做测定轨之外，也能利用总相位抖动做空间介质等科学研究。

致谢

感谢曲春凯、金炜桐、邓青云、郭茜和冯鹏等人利用13 m天线收集深空探测器的数据。感谢王云博和杨永章的建议与帮助。