# Planetary Radio Interferometry and Doppler Experiment (PRIDE) technique: A test case of the Mars Express Phobos fly-by

D. Duev, S. Pogrebenko, G. Cimò, G. Molera Calvés, T. Bocanegra Bahamón, L. Gurvits, M. Kettenis, J. Kania, V. Tudose, P. Rosenblatt, et al.

► To cite this version:

D. Duev, S. Pogrebenko, G. Cimò, G. Molera Calvés, T. Bocanegra Bahamón, et al. Planetary Radio Interferometry and Doppler Experiment (PRIDE) technique: A test case of the Mars Express Phobos fly-by. Astronomy & Astrophysics - A&A, 2016, 593, pp.A34. 10.1051/0004-6361/201628869. hal-02320320

HAL Id: hal-02320320

https://hal.science/hal-02320320v1

Submitted on 13 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Planetary Radio Interferometry and Doppler Experiment (PRIDE) technique: A test case of the Mars Express Phobos fly-by

D. A. Duev<fnref n="1" />,²,³, S. V. Pogrebenko², G. Cimò²,⁴, G. Molera Calvés⁵,², T. M. Bocanegra Bahamón⁶,²,⁷, L. I. Gurvits²,⁶, M. M. Kettenis², J. Kania⁸,², V. Tudose⁹, P. Rosenblatt¹⁰, J.-C. Marty¹¹, V. Lainey¹², P. de Vicente¹³, J. Quick¹⁴, M. Nickola¹⁴, A. Neidhardt¹⁵, G. Kronschnabl¹⁵, C. Ploetz¹⁵, R. Haas¹⁶, M. Lindqvist¹⁶, A. Orlati¹⁷, A. V. Ipatov¹⁸, M. A. Kharinov¹⁸, A. G. Mikhailov¹⁸, J. E. J. Lovell¹⁹, J. N. McCallum¹⁹, J. Stevens²⁰, S. A. Gulyaev²¹, T. Natush²¹, S. Weston²¹, W. H. Wang⁷, B. Xia⁷, W. J. Yang²², L.-F. Hao²³, J. Kallunki²⁴, and O. Witasse²⁵

*(Affiliations can be found after the references)*

Received 7 May 2016 / Accepted 31 May 2016

## ABSTRACT

*Context.* The closest ever fly-by of the Martian moon Phobos, performed by the European Space Agency's Mars Express spacecraft, gives a unique opportunity to sharpen and test the Planetary Radio Interferometry and Doppler Experiments (PRIDE) technique in the interest of studying planet-satellite systems.

*Aims.* The aim of this work is to demonstrate a technique of providing high precision positional and Doppler measurements of planetary spacecraft using the Mars Express spacecraft. The technique will be used in the framework of Planetary Radio Interferometry and Doppler Experiments in various planetary missions, in particular in fly-by mode.

*Methods.* We advanced a novel approach to spacecraft data processing using the techniques of Doppler and phase-referenced very long baseline interferometry spacecraft tracking.

*Results.* We achieved, on average, mHz precision (30 µm/s at a 10 s integration time) for radial three-way Doppler estimates and sub-nanoradian precision for lateral position measurements, which in a linear measure (at a distance of 1.4 AU) corresponds to ~50 m.

**Key words.** techniques: interferometric – techniques: miscellaneous – methods: data analysis – astrometry

## 1. Introduction

The Planetary Radio Interferometry and Doppler Experiments (PRIDE) project (Duev et al. 2012) initiated by the Joint Institute for VLBI ERIC (JIVE, Dwingeloo, The Netherlands) utilises Doppler and phase-referencing very long baseline interferometric (VLBI<fnref n="1" />) observations to provide high precision spacecraft state vector estimation. In this paper, we further advance a novel approach to spacecraft data processing developed within PRIDE.

On December 29, 2013, the European Space Agency's Mars Express (MEX) spacecraft made the closest ever fly-by of Phobos, one of the two Martian moons, just some 45 km from its surface (Witasse et al. 2014). PRIDE observations of MEX during this event, involving more than 30 radio telescopes spread around the globe, were carried out by our team on December 28–29, 2013 (PI Pascal Rosenblatt, Royal Observatory of Belgium; European VLBI Network (EVN)/Global VLBI experiment code GR035). These observations allow reconstructing MEX's trajectory in the vicinity of Phobos with a high accuracy, which will in turn help to put a better constraint on the geophysical parameters of Phobos, possibly shedding light on its origin. The PRIDE data processing technique has been specifically refined for the observations of MEX during this event to provide high precision positional and Doppler measurements. In particular, we describe here the improvements made to the correlator software

at JIVE (Keimpema et al. 2015) that allow efficient handling of such data, and we demonstrate the positive impact of these enhancements on the spacecraft position estimates obtained at the post-processing stage. The use of the Doppler and VLBI data in dynamical modelling of MEX motion to estimate the geophysical parameters related to the interior composition of this Mars moon is discussed in Rosenblatt et al. (in prep.).

The GR035 experiment has been conducted as a live end-to-end verification of the PRIDE technique which will be used in future planetary missions, in particular ESA's Jupiter Icy Satellites Explorer, JUICE (Grasset et al. 2013). In this paper, we present the technique, including the data processing algorithms. In Sect. 2, we describe the set-up of the GR035 experiment. Section 3 describes Doppler and VLBI data processing pipeline and presents the results of the experiment. Section 4 provides the reader with conclusions and outlook. The scientific evaluation of the experiment will be given in a separate paper (Rosenblatt et al., in prep.). The experiment GR035 was complementary to the nominal MEX Radio Science experiment MaRS (Pätzold et al. 2004, 2016).

## 2. Experiment GR035 set-up

The error in the a priori position of a phase calibrator directly affects the estimates of the MEX orbit parameters proportionally to the separation angle between the calibrator and the target. Therefore, in order to reach the best positional accuracy of MEX achievable with the modern ground-based VLBI – about a nanoradian (0.2 mas), which is translated to ~100 m at the orbit, phase calibrators are necessary within ~2°

<footnote n="1">¹ Very long baseline interferometry (VLBI) is a technique in which the signals from a network of radio telescopes, spread around the world, are combined to simulate one single telescope with a resolution far surpassing that of the individual telescopes.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Mars Express Phobos fly-by experiment field on the sky. A large black star denotes the field centre; three crosses denote the pointings used in the experiment ET027 to observe possible secondary calibrators; five smaller black stars denote the detected sources; the circle represents the approximate size of the primary beam of a 30 m antenna. The colour insets show calibrated images of three of the sources in the field obtained in the ET027 experiment.</figcaption></figure>

(Beasley & Conway 1995) of the MEX position at the observational epoch with absolute errors in their position not exceeding 0.2 mas.

In the case of the GR035 experiment, the timing of the observations and therefore the sky position of the events were, obviously, defined by MEX ballistics. The nearest bright source suitable as a primary phase calibrator, J1232–0224, with 0.16 mas error in RA, 0.27 mas error in Dec and 0.717 Jy total X-band flux density (Petrov 2015), happened to be at a fairly large angular distance of ~2.5° from the target. For this reason, on March 19, 2013, we performed VLBI observations of the fly-by event field with the EVN (experiment code ET027). We were able to identify three new sources, one of which – “CAL5” or J1243–0218 (0.133 Jy total flux density; see Fig. 1) – appeared to be suitable both as a secondary calibrator and as an arc stability test source to ensure the quality of the final astrometrical solution for MEX.

The observations of the fly-by event were carried out within the global VLBI experiment GR035. During the run on December 28–29, 2013, we observed three consecutive revolutions of MEX around Mars, each 7 h long, in order to provide a better coverage for further orbit reconstruction. The total duration of the experiment was 25 h. Telescopes that took part in the observations<fnref n="2" /> were split into two sub-arrays. The antennas in the first sub-array (see Fig. 2, stations depicted with triangles) observed in Doppler mode with a dual S/X-band (2/8.4 GHz) frequency set-up with long scans (20-min) to minimise the Doppler frequency detection noise. The second sub-array antennas (see Fig. 2, stations depicted with circles) used an X-band (8.4 GHz) set-up and operated in a phase-referencing VLBI mode with short scans (Ros et al. 1999), alternating between the target and the phase calibrators. The frequency set-up used by each of the stations is shown in Fig. 2. Standard VLBI recording equipment was used at all stations. Strong sources M87 and J1222+0413

with total X-band flux densities of 1.968 and 0.701 Jy, respectively (Petrov 2015), were used as fringe<fnref n="3" /> finders.

The Australian, New Zealand and eastern EVN stations began the tracking, then subsequently the western EVN and VLBA stations stepped in. This is illustrated in Fig. 3.

### 3. Doppler and VLBI data processing pipeline and results

We have advanced the generic spacecraft data processing pipeline developed within the scope of PRIDE (outlined in Fig. 4) to achieve a very high precision of MEX positional estimation, which is shown schematically in Fig. 5.

First, we performed narrowband processing of the single-dish open-loop data collected by all of the stations. We used the SWSpec/SCtracker/dPLL<fnref n="4" /> (Molera Calvés 2012; Molera Calvés et al. 2014) software package to obtain the topocentric Doppler detections. SWSpec extracts the raw data from the channel where the spacecraft carrier signal is expected to be recorded. Then it performs a window-overlapped add (WOLA) discrete Fourier transform (DTF) and time integration over the obtained spectra. The result is an initial estimate of the spacecraft carrier tone along the scan. The moving phase of the spacecraft carrier tone throughout the scan is modelled by performing an *n*-order frequency polynomial fit. SCtracker uses this initial fit to stop the phase of the carrier tone, allowing subsequent tracking, filtering and extraction of the carrier tone in narrower bands (from the initial 16 MHz channel bandwidth down to a 2 kHz bandwidth) using a second-order WOLA DFT-based algorithm of the Hilbert transform approximation. The Digital Phase-Locked-Loop (dPLL) performs high precision reiterations of the previous steps – time-integration of the overlapped spectra, phase polynomial fitting, and phase-stopping correction – on the 2 kHz bandwidth signal, using 20 000 FFT points and 10-s integration time. The output of the dPLL is the filtered down-converted signal and the final residual phase in the stopped band with respect to the initial phase polynomial fit. The bandwidth of the output detections is 20 Hz with a frequency spectral resolution of 2 mHz. The Doppler observable is obtained by adding the base frequency of the selected channel to the 10-s averaged carrier tone frequencies retrieved by the dPLL.

During this experiment, three transmitting stations provided 24-h coverage: the 35-m ESTRACK station New Norcia (NNO) in Australia and the 70-m Deep Space Stations 63 (DSS-63) in Robledo (Spain) and 14 (DSS-14) in Goldstone (CA, USA). In order to estimate the Doppler noise $f_n$, the topocentric detections with a 10-s integration time for each station were differed with predicted three-way Doppler values, and the standard deviation of the result was calculated for each 2-min scan (see the histogram in Fig. 6). The resulting mean value of $f_n$ is 2.5 mHz, median –2.2 mHz, and mode (maximum of a fitted log-normal distribution) –1.7 mHz. The mode value translates to $0.5 \cdot c \cdot (f_n/f_0) = 30 \ \mu\text{m/s}$ in linear measure for the three-way Doppler, where $f_0 = 8.4 \text{ GHz}$, and *c* is the speed of light in a vacuum. This is comparable to the precision of the Doppler detections provided by the DSN and ESOC (see e.g. Tyler et al. 1992; Budnik et al. 2004).

<footnote n="2">² Additional information about particular antennas can be found by the two-letter codes shown in Fig. 2 in the databases of, e.g. the International VLBI Service at ftp://ivscc.gsfc.nasa.gov/pub/control/ns-codes.txt and the European VLBI Network at http://www.evlbi.org/user_guide/EVNstatus.txt</footnote>

<footnote n="3">³ The response of a VLBI system.</footnote>

<footnote n="4">⁴ Wagner, J., Molera Calvés, G., and Pogrebenko, S.V. 2009–2014, Metsähovi Software Spectrometer and Spacecraft Tracking tools, Software Release, GNU GPL, http://www.metsahovi.fi/en/vlbi/spec/index</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Experiment GR035: participating telescopes are denoted by their two-letter codes. Stations depicted with a triangle were engaged in the Doppler part of the experiment, with a circle in the phase-referencing VLBI part. Station colour denotes frequency set-up used: red – 8332–8444 MHz, 8 channels of 16 MHz USB (Upper Side Band); blue – 8396–8444 MHz, 4 channels of 16 MHz USB; orange – 8380–8428 MHz, 4 channels of 16 MHz USB; green – 2289–2305/8396–8412 MHz, 4 channels of 16 MHz USB.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Observational time ranges (UTC) of the telescopes that participated in the phase-referencing VLBI part. Green indicates the Australian, New Zealand and eastern EVN stations, blue the western EVN stations, and red the VLBA stations. Station codes as in Fig. 2. The dashed black vertical line denotes the time of the fly-by event as seen from the Earth.</figcaption></figure>

In order to process the VLBI data, streams from each station must be synchronised with a common base, usually the International Atomic Time TAI. The behaviour of station clocks is regularly checked against the GPS<fnref n="5" /> time scale tied to the TAI time. We examined these time series and chose Medicina (station code Mc) as the absolute reference station for the current experiment, which appeared to have the best long-term stability and the smallest absolute clock rate value around the date of the experiment. Clock parameters of the rest of the stations were referenced to Medicina using the fringe finder data<fnref n="6" />.

The spacecraft cross-correlation spectrum is smeared owing to the intrinsic change of the frequency (emitted by a spacecraft as it retransmits the signal in the two-/three-way Doppler

<figure><img src="image_5.png" /><figcaption>Fig. 4. Generic PRIDE data flow and processing pipeline. Spacecraft Doppler and delay observables obtained, respectively, at the narrow-band and post-processing steps, alongside the final astrometric solution, that can be used for a variety of scientific applications (Duev et al. 2012).</figcaption></figure>

regime) caused by a change in the relative velocity. To mitigate this frequency smearing, a Doppler phase correction must be applied to the spacecraft data. In VLBI, data from individual telescopes are reduced to a common phase centre, usually the geocentre. Therefore, to compute an empirical Doppler phase correction, we first reduced all the topocentric frequency detections $f_{tc}(t)$ to the geocentre (see Fig. 7) using the equation

$$f_{gc}(t) = f_{tc}(t - \tau_{gc}) \cdot \left(1 - \frac{d\tau_{gc}}{dt}\right), \qquad (1)$$

where $t$ is UTC time, and $\tau_{gc}$ and $d\tau_{gc}/dt$ are the total near-field VLBI signal delays and delay rates with respect to the geocentre.

<footnote n="5">⁵ Global Positioning System.</footnote>

<footnote n="6">⁶ This procedure is usually referred to as the “clock search”.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** GR035 data processing pipeline. The positional measurements and the Doppler detections are fed into a dynamical model of MEX motion (Rosenblatt et al., in prep.).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Doppler detection noise in mHz at X-band. Topocentric detections with a 10-s integration time for each station were differenced with predicted values, and the standard deviation of the result was calculated for each 2-min scan providing what is referred to as the measurements in the histogram.</figcaption></figure>

The resulting geocentric frequencies (consistent with each other and with the geocentric frequency prediction at a sub-mHz level) are subsequently averaged and integrated into phases on a per scan basis. This way, phases are reset to zero at the beginning of each scan. This is done to avoid numerical errors associated with a fairly wide dynamical range of changes in the phase.

The resulting phase correction is applied at the next processing step, the broadband correlation with the EVN software correlator at JIVE, SFXC (Keimpema et al. 2015). For the correlation of both the calibrator and the spacecraft data, we used the signal delay models described in Duev et al. (2012). These models have been implemented in a software package *pypride* (Python tools for PRIDE), whose output is compatible with the SFXC. The package is mostly written in the Python programming language with an extensive use of modules providing JIT-compilation to boost performance. The most computationally expensive sub-routines are written in Fortran. Most of the tasks are automated and parallelised. The package *pypride* calculates VLBI delays and the uvw-projections of baselines/Jacobians for far- and near-field sources and for space VLBI (for details, see Duev et al. 2012, 2015). The software can also be used to calculate Doppler frequency shift predictions for spacecraft observations.

Changes in spacecraft signal spectrum over time are due to different transmission modes used during a communication session (see Fig. 8). A straightforward approach to the correlation of such data at typical resolutions used in VLBI works only if the

so-called data-bands are present in the spacecraft spectrum (see Fig. 9, left – characteristic “bumps” around the carrier and the first sub-carriers). However, this approach effectively narrows the bandwidth by a factor of ~3, which results in higher noise in the group delay estimates at the next processing step. More importantly, however, these spectral features are not constant over time and may completely disappear for extended intervals. To overcome these difficulties, we realised an approach that makes use of the ranging tones (sub-carriers) present in the S/C spectrum most of the time; the phases of these tones are directly related to the carrier phase as they are all synthesised from the same reference signal. This allows most of the available bandwidth to be used. The Doppler phase correction stops the cross-correlation spectrum drift (see Fig. 10). Individual sub-carrier lines are clearly seen only at a sufficiently large spectral resolution owing to their intrinsic narrowness. Therefore, we first correlated data on several baselines with larger telescopes of the array using a very high spectral resolution and derived a spectral mask leaving one to four spectral points to a line depending on its width at that particular resolution. For an initial mask approximation, we used a peak identification approach incorporating continuous wavelet transform-based pattern matching (Du et al. 2006), after which the resulting mask was inspected and corrected by hand if necessary. To derive the optimal spectral resolution, we tried fringe-fitting the results of correlation at different resolutions on several baselines (see Fig. 11 for examples). The standard deviation of the 1-s integrated group delay estimates suggested an optimal spectral resolution value of $2^{19} = 524\,288$ points. Presumably, numerical effects come into play at higher resolutions, preventing a further increase in precision. Finally, the amplitudes of individual filtered lines are normalised to unity, while the phases are kept intact. This provides additional improvement in the precision of group delay estimation by making the fringes more pronounced in the lag domain.

We need to point out that for ~50% of the time during GR035, only the carrier line was present in the spacecraft spectrum (see Fig. 9, bottom; the mask here consists of a few points around the carrier). In this case, preserving this only feature in the spectrum does not allow group delay estimation owing to an extremely narrow effective bandwidth; however the phase may still be accurately extracted and used.

Correlation at such a high spectral resolution results in a massive amount of data. Therefore, in order to reduce the latter to a manageable level, the resulting spectra are compressed to a resolution of 256 spectral points. This is achieved by transforming the spectra into the lag domain, where 512 points around the central lag are cut out. These spectra are then transformed back into the frequency domain (see Fig. 12). The secondary peaks seen in Fig. 12 (right) around the central (true) fringe are the result of the compression operation at the previous stage, which is equivalent to convolving the comb-like spectra<fnref n="7" /> in frequency domain with a Fourier image of a set of rectangular windows, and can easily confuse the fringe fitting algorithm. However, after performing the clock search, we know that the fringe maximum in the lag domain lies within several lags around the zeroth lag, and so we can eliminate these “out-of-the-maximum” peaks by applying a squared cosine-window filter (shown by a black line in Fig. 12, right). Our tests have shown that fringe fitting the output of this compression procedure yields the same

<footnote n="7">⁷ A spectrum with filtered and normalised spectral lines resembles a comb.</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Topocentric frequency detections (**top**, Hz) and frequency detections reduced to a common phase centre – geocentre (**bottom**, Hz). Jumps in frequency at 03:30 and 11:30 UT are due to the uplink frequency changes at transmitting ground stations. Mean geocentric frequency was converted to phase and applied to the spacecraft signal at correlation to avoid a frequency smearing at high spectral resolution. Station two-letter codes as in Fig. 2.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** MEX signal spectrum types. “Comb” denotes the time intervals when ranging tones were present in the spectrum; “Carrier” when only the carrier was present in the spectrum.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Zoom into the carrier line without (**left panel**) and with the Doppler phase correction (**right panel**). 10 s integration time, $2^{15} = 32\,768$ points spectral resolution, baseline Hh-Ww, scan 135. 23:22–23:24 UTC, December 28, 2013.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Averaged amplitude spectrum after applying the Doppler phase correction in arbitrary units, $2^{14} = 16\,384$ points spectral resolution, baseline T6-Sv, scan 209 (**top**) and 191 (**bottom**). Only the carrier line was present in the spectrum in the second case, as was the case for ~50% of the time during GR035. The spectral mask is shown in orange dots.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Optimal FFT size to perform signal filtration as characterised by the standard deviation of the group delay estimates obtained after fringe-fitting as a function of the spectral resolution used at correlation. Spectral resolution ranges from $2^{10} = 1024$ to $2^{21} = 2\,097\,152$ points. Baselines T6-Sv, T6-Mc, Sv-Mc.</figcaption></figure>

precision of group delay estimates as in the case of the original non-compressed data.

The spectrum filtration and compression described above were first implemented in Python and thoroughly tested before being incorporated into the SFXC correlator. This implementation is integrated in the SFXC standard spectral averaging code
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Spectrum compression from 4096 (**left panel**) down to 256 (**right panel**) spectral points, and filtration in lag domain. The solid black line in the right panel shows the output filter profile. One-second integrated spectra are shown. Scan 209, baseline T6-Sv. The dashed black lines denote the central lag.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** *Left panel:* full time range uv-coverage for all sources. *Middle panel:* CLEAN’ed map of MEX integrated over the time range from 29/12 11:30–13:30 UTC. *Right panel:* self-calibrated CLEAN’ed map of the primary calibrator source J1232–0224 to the same scale and integrated over the same time range.</figcaption></figure>

and allows the arbitrary spectral and window filters to be specified as appropriately sized vectors.

The output cross-correlation spectra in the SFXC correlator-specific format were converted into the Measurement Set format (Kemball & Wieringa 2000) and into the FITS IDI files for further processing.

To derive the displacements of MEX from its a priori position in the post processing analysis of the data, we employed two different independent approaches: imaging and solving the astrometric measurement equation.

The first approach was realised employing a commonly used VLBI data reduction package AIPS (Greisen 2003). The FITS-files are loaded into the AIPS file system using a task<fnref n="8" /> FITLD. After a preliminary data inspection and editing, initial calibration is applied. This includes bandpass (AIPS task BPASS) calibration using the source J1222+0413 and antenna calibration (task ANTAB). For the stations that did not provide system temperature measurements during the experiment, nominal values were used to calibrate the antenna. To correct the delays and rates of the phase referencing calibrator J1232–0224, fringe fitting is performed with the task FRING. Next, a procedure called

“self-calibration” is applied to the calibrator using the CALIB task. During self-calibration, phase corrections for the antennas are calculated based on a model of the source. We started with a point-source model. Then a CLEAN’ed (Högbom 1974) map of the calibrator is produced using the AIPS task IMAGR. The resulting map is used as a new model. An example CLEAN’ed self-calibrated map of the source J1232–0224 is shown in Fig. 13, right panel. When we are satisfied with the map and with the calibration after a number of iterations of IMAGR and CALIB, the resulting phase corrections are applied to the spacecraft. At this stage it is possible to make an image of the spacecraft. The process of self-calibration fixes the centre of the map to the nominal a priori position<fnref n="9" />. However, the position of the spacecraft on

<footnote n="9">⁹ For 80% of the calibrators, the position is known to an accuracy better than 3 mas. According to Shu et al. (2016), 1167 calibrators were known within 7.5° of the ecliptic band by May 2016, and their number is growing. The median accuracy of their positions is 0.45 mas. A source is considered a calibrator if its median correlated flux density at baselines longer than 5000 km is above 30 mJy at 8 GHz. A dedicated observing program for improving positions of all known calibrators to a level of 0.3 mas is underway (Shu et al. 2016). Potentially, the accuracy can be further improved to reach a level of 0.1 mas if the necessary resources are allocated.</footnote>

<footnote n="8">⁸ In the AIPS environment, separate sub-programs are called “tasks”.</footnote>
<figure><img src="image_15.png" /><figcaption>Fig. 14. Displacements from the a priori lateral position of MEX as a function of time, measured using the imaging approach. Displacements in Right Ascension (mas) are shown in blue, in Declination (mas) in red. 2-min integration time; net time on target ~5.5h. December 28–29, 2013. The apparent jump in position near 01:00:00 UT is discussed in Sect. 3.</figcaption></figure>

the final image will have a shift from the centre of the map (defined by the spacecraft a priori position) due to the errors in the spacecraft ephemeris. The position measured on the map provides the actual coordinates of the spacecraft for each solution interval. The measurements are done with the task JMFIT by fitting a 2D Gaussian to the peak on the CLEAN’ed spacecraft map located within a box set using the image integrated over the full time range for the current sub-array (which efficiently suppresses the side lobes and any quickly changing variations in phase). An example image of MEX is shown in Fig. 13 (middle panel).

The described pipeline was automated with a ParselTongue script (Kettenis et al. 2006). When imaging the scans where only the carrier line was present in the MEX spectrum, we used only a few frequency channels around the carrier in a manner analogous to that used in spectral-line VLBI.

The alternative approach we used for estimating the MEX positional displacement (which we realised within the *pypride* software package) is similar to the approach used in geodetic VLBI and is based on solving the measurement equation for each epoch *t*.

$$ \Delta\phi|_{t} = (J \cdot \Delta\alpha)|_{t}, \qquad (2) $$

where $\Delta\phi$ is a vector of differential MEX carrier line phases on baselines, *J* is a matrix containing near-field analogues of uv-projections of baselines (see Duev et al. 2012), and $\Delta\alpha$ is the vector of corrections to the a priori lateral position of the spacecraft. The phases $\Delta\phi$ are subject to a 2π-ambiguity, which means that the corresponding phase delays $\tau_{ph} = \varphi/\omega_0$ ($\omega_0 = 2\pi f_0$) may have a bias of several cycles of ~120 ps for observations at X-band. In order to solve for this ambiguity, we employed the following approach.

First, the calibrator data are fringe-fitted and self-calibrated providing calibration group delays ($\tau_{gr} = d\phi/d\omega$, $\omega = 2\pi f$) and phases, which are applied to the MEX data. The rms error of the calibrator group delay estimates are used to set baseline weights. Then, the phase of the MEX carrier line is extracted and unwrapped. A naive approach to unwrapping does not work in most cases owing to large gaps and uneven spacing in time, and because the phase slope is not constant. Therefore, we made use of a wrapped Kalman smoother, which is based on a wrapped Kalman filter (WKF) algorithm described in Traa & Smaragdis (2013, see Appendix A). The calibrated MEX data are subsequently fringe-fitted. This yields residual group delays for the

scans when the sub-carriers are present in its spectrum. An SVM<fnref n="10" />-based unsupervised machine learning algorithm is used to identify and flag outliers in group delay estimates. An optimal fit of phase delays $\tau_{ph}$ to group delays is found by minimizing the squared error defined as $\sum_n \sum_i ((\tau_{ph}[i] + 2\pi n) - \tau_{gr}[i])^2$ over time periods when both phase and group delays are available with respect to $n \in \mathbb{N}$. This yields the number of phase cycles $n_m$ in question providing a solution to the 2π-ambiguity problem, also for the time intervals when no group delay data are available. If there is such a long gap in time in the data that the Kalman smoothing procedure fails to correctly unwrap the phase, the group delays are automatically split into an appropriate number of clusters using a DBSCAN algorithm (Ester et al. 1996), and then the phase delays are fitted to the corresponding clusters. This procedure yields an “unambiguous” phase to be used in the astrometric Eq. (2) solution.

A proper transition from group delays<fnref n="11" />, which are commonly used in VLBI astrometry, allows the Eq. (2) solution error to be reduced by an order of magnitude.

Equation (2) is in most cases overdetermined, therefore we used the singular value decomposition of the matrix *J* when solving it. The resulting angular corrections $\Delta\alpha$ are not displacements in Right Ascension and Declination per se (because they are defined for sources at infinity), but the angular displacement of the vector from geocentre to the target at a given epoch (see Duev et al. 2012).

In both approaches described above, the solution interval was set to the MEX scan length (2 min in most cases); for each target scan, two adjacent calibrator scans were used to perform calibration. In the processing, we used an elevation cut-off angle of 20°. Three stations were set as reference stations: Bd (from the start of the experiment until ~December 29, 2013 00:30 UTC), Ys (from ~01:00 until ~07:30 UTC), and Pt (from ~08:00 UTC to the end of the experiment). In the first part of the observations, the largest telescope of that sub-array (see Fig. 3, stations denoted in green), T6, had to be dropped because of phase stability problems at the station. In addition, most of the stations of the sub-array did not provide satisfactory system temperature $T_{sys}$.

<footnote n="10">¹⁰ Support vector machines.</footnote>

<footnote n="11">¹¹ The spacecraft signal is bandwidth limited to about 12 MHz in our case (see Fig. 9), which sets a certain limit on the group delay estimation accuracy.</footnote>
measurements, which resulted in poor data calibration and a consequent astrometric solution bifurcation. In Fig. 14 (before ~December 29, 2013 00:30 UTC) the solution nearest to the phase centre is shown. However, the solutions for the second and third sub-arrays with Ys and Pt as reference stations (Fig. 14 after ~December 29, 2013 01:00 UTC) show no bifurcation. Finally, the results of the imaging and “geodetic” approaches are consistent at a level of ~10 micro-arcseconds for the best-calibrated time range (from ~December 29, 2013 01:00 UTC onwards). The median $3\sigma$ formal error values for the full time range are 0.034 mas for RA and 0.058 mas for Dec, which translates into ~35 and 60 m at the orbit, respectively (MEX was at a distance of ~1.4 AU during the experiment).

We checked the stability of our array of telescopes by performing imaging of the source CAL5 (J1243–0218) using all available observations and the data of different sub-arrays. The derived astrometric position in all cases appeared to be the same as in the experiment ET027 with an uncertainty of less than 0.1 mas. At the same level of uncertainty, the coordinates are consistent with the Radio Fundamental Catalogue (Petrov 2015) values (RA $12^h43^m52.487864^s$; Dec $–02^d18’38.401056”$).

## 4. Conclusions and outlook

In this work, we were able to measure the lateral position and radial Doppler of the MEX spacecraft with a precision of about 50 m and 30 µm/s, respectively. This is comparable to what has been reported for other spacecraft (e.g. Jones et al. 2015). These measurements are used by our collaborators from the Royal observatory of Belgium (ROB) and the French National Centre for Space Studies (CNES) in the dynamical modelling of Mars Express motion aimed at estimating the geophysical parameters of Phobos (see Rosenblatt et al., in prep.).

The offsets in the estimated MEX RA and Dec of ~1 mas from the a priori orbital position seen in Fig. 14, although comparable to the formal orbit determination (OD) error budget, require further investigation and calibration. Similar results have been reported in the past for other spacecraft (e.g. Lanyi et al. 2005). To calibrate these offsets, observations of multiple spacecraft/calibrator pairs with telescope networks of different configurations will be required. In addition, models of signal delay caused by propagation effects must be improved. These are necessary since the systematic offsets most likely result from a combination of a number of error sources, but are dominated by the uncertainty in the a priori spacecraft position and propagation effects.

PRIDE was selected by ESA as one of the experiments of its L-class JUpiter ICy moons Explorer mission (JUICE; Grasset et al. 2013). The spacecraft data acquisition, processing, and analysis pipelines developed in this work will create a basis for implementation of PRIDE-JUICE.

*Acknowledgements.* We would like to express our gratitude to the anonymous referee. The authors acknowledge the EC 7th Framework Programme (FP7/2008-2017) project ESPaCE (grant agreement #263466). Mars Express is a mission of the European Space Agency. Information about Mars Express telecommunication were provided by the Mars Express Project. The European VLBI Network is a joint facility of independent European, African, Asian, and North American radio astronomy institutes. Scientific results from data presented in this publication are derived from the following EVN project code: GR035. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Australia Telescope Compact Array is part of the Australia Telescope National Facility which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. This study made use of data collected through the AuScope initiative. AuScope Ltd is funded under the

National Collaborative Research Infrastructure Strategy (NCRIS), an Australian Commonwealth Government Programme. The authors would like to thank the personnel of the participating stations. R. M. Campbell, A. Keimpema, P. Boven (JIVE), M. Pätzold (University of Cologne), B. Häusler (University of BW Munich), and D. Titov (ESA/ESTEC) provided important support to various components of the project. G. Cimó acknowledges the Horizon 2020 project AS-TERICS. T. Bocanegra Bahamon acknowledges the NWO-ShAO agreement on collaboration in VLBI. J. Kania acknowledges the ASTRON/JIVE Summer Studentship programme. P. Rosenblatt is financially supported by the Belgium PRODEX program managed by the European Space Agency in collaboration with the Belgian Federal Science Policy Office.

## References

Beasley, A. J., & Conway, J. E. 1995, in Very Long Baseline Interferometry and the VLBA, eds. J. A. Zensus, P. J. Diamond, & P. J. Napier, ASP Conf. Ser., 82, 327

Budnik, F., Morley, T. A., & MacKenzie, R. A. 2004, in 18th International Symposium on Space Flight Dynamics, ESA SP, 548, 387

Du, P., Kibbe, W. A., & Lin, S. M. 2006, Bioinformatics, 22, 2059

Duev, D. A., Molera Calvés, G., Pogrebenko, S. V., et al. 2012, A&A, 541, A43

Duev, D. A., Zakhvatkin, M. V., Stepanyants, V. A., et al. 2015, A&A, 573, A99

Ester, M., Kriegel, H.-P., Sander, J., & Xu, X. 1996, in Second International Conference on Knowledge Discovery and Data Mining (KDD-96) (AAAI Press), 226

Grasset, O., Dougherty, M. K., Coustenis, A., et al. 2013, Planet. Space Sci., 78, 1

Greisen, E. W. 2003, Information Handling in Astronomy – Historical Vistas, 285, 109

Högbom, J. A. 1974, A&AS, 15, 417

Jones, D. L., Folkner, W. M., Jacobson, R. A., et al. 2015, AJ, 149, 28

Keimpema, A., Kettenis, M. M., Pogrebenko, S. V., et al. 2015, Exper. Astron., 39, 259

Kemball, A. J., & Wieringa, M. H., eds. 2000, MeasurementSet definition, ver. 2.0

Kettenis, M., van Langevelde, H. J., Reynolds, C., & Cotton, B. 2006, in Astronomical Data Analysis Software and Systems XV, eds. C. Gabriel, C. Arviset, D. Ponz, & S. Enrique, ASP Conf. Ser., 351, 497

Lanyi, G., Border, J., Benson, J., et al. 2005, Interplanetary Network Progress Report, 162, 1

Molera Calvés, G. 2012, Ph.D. Thesis, Aalto University, Pub. No. 42/2012

Molera Calvés, G., Pogrebenko, S. V., Cimò, G., et al. 2014, A&A, 564, A4

Pätzold, M., Neubauer, F. M., Carone, L., et al. 2004, in Mars Express: the Scientific Payload, eds. A. Wilson, & A. Chicarro, ESA SP, 1240, 141

Pätzold, M., Häusler, B., Tyler, G., et al. 2016, Planetary and Space Science, 127, 44

Petrov, L. 2015, Radio Fundamental Catalog, http://astrogeo.org/rfc, accessed: 2015-08-01

Ros, E., Marcaide, J. M., Guirado, J. C., et al. 1999, A&A, 348, 381

Shu, F., Petrov, L., Jiang, W., et al. 2016, ArXiv e-prints [arXiv:1605.07036]

Traa, J., & Smaragdis, P. 2013, IEEE Signal Processing Letters, 20, 1257

Tyler, G. L., Balmino, G., Hinson, D. P., et al. 1992, J. Geophys. Res., 97, 7759

Witasse, O., Duxbury, T., Chicarro, A., et al. 2014, Planet. Space Sci., 102, 18

<footnote n="1">1. California Institute of Technology, 1200 E California Blvd, Pasadena, CA 91125, USA
e-mail: duev@caltech.edu</footnote>

<footnote n="2">2. Joint Institute for VLBI ERIC, PO Box 2, 7990 AA Dwingeloo, The Netherlands</footnote>

<footnote n="3">3. Sternberg Astronomical Institute, Lomonosov Moscow State University, Universitetsky av. 13, 119991 Moscow, Russia</footnote>

<footnote n="4">4. ASTRON, The Netherlands Institute for Radio Astronomy, Postbus 2, 7990 AA, Dwingeloo, The Netherlands</footnote>

<footnote n="5">5. Aalto University, School of Electrical Engineering, Department of Radio Science and Engineering, 02150 Espoo, Finland</footnote>

<footnote n="6">6. Department of Astrodynamics and Space Missions, Delft University of Technology, 2629 HS Delft, The Netherlands</footnote>

<footnote n="7">7. Shanghai Astronomical Observatory, 80 Nandan Road, 200030 Shanghai, PR China</footnote>

<footnote n="8">8. Carnegie Mellon University, 5000 Forbes Ave, Pittsburgh, PA 15213, USA</footnote>
⁹ Institute for Space Sciences, Atomistilor 409, PO Box MG-23, 077125 Bucharest-Magurele, Romania

¹⁰ Royal Observatory of Belgium, Ringlaan 3, 1180 Brussels, Belgium

¹¹ CNES/GRGS, OMP 14 avenue Édouard Belin 31400 Toulouse, France

¹² IMCCE, Observatoire de Paris, PSL Research University, CNRS-UMR 8028 du CNRS, UPMC, Lille-1, 77 Av. Denfert-Rochereau, 75014 Paris, France

¹³ Observatorio de Yebes (IGN), Apartado 148, 19180 Guadalajara, Spain

¹⁴ Hartebeesthoek Radio Astronomy Observatory, 1740 Krugersdorp, South Africa

¹⁵ Federal Agency for Cartography and Geodesy, Geodetic Observatory of Wettzell, 60598 Frankfurt Am Main, Germany

¹⁶ Department of Earth and Space Sciences, Chalmers University of Technology, Onsala Space Observatory, 439 92 Onsala, Sweden

¹⁷ National Institute for Astrophysics, Radio Astronomy Institute, Radio Observatory Medicina, 75500 Medicina, Italy

¹⁸ Institute of Applied Astronomy, Russian Academy of Sciences, Kuzuzova Embankment 10, 191187 Saint-Petersburg, Russia

¹⁹ School of Physical Sciences, University of Tasmania, Private Bag 37, 7001 Hobart, Australia

²⁰ CSIRO Astronomy and Space Science, Australia Telescope National Facility, Narrabri NSW 2390, Australia

²¹ Institute for Radio Astronomy and Space Research, Auckland University of Technology, Private Bag 92006, 1142 Auckland, New Zealand

²² Xinjiang Astronomical Observatory, Chinese Academy of Sciences, 830011 Urumqi, PR China

²³ Yunnan Astronomical Observatory, Chinese Academy of Sciences, 650011 Kunming, PR China

²⁴ Metsähovi Radio Observatory, Aalto University, 02540 Kylmälä, Finland

²⁵ European Space Agency, ESA/ESTEC Scientific Support Office, 2200AG Noordwijk, The Netherlands
# Appendix A: Wrapped Kalman smoother

The wrapper Kalman filter (WKF) is a Kalman filter for which the filtered state distribution $P_{WN}$ is represented by a wrapped Gaussian (Traa & Smaragdis 2013):

$$P_{WN}(\varphi | \mu, \sigma^2) = \frac{1}{\sigma\sqrt{2\pi}} \sum_{l=-\infty}^{\infty} \exp\left[-\frac{(\varphi - (\mu + 2\pi l))^2}{2\sigma^2}\right], \quad \varphi \in \mathbb{S}^1 \qquad (A.1)$$

The wrapped Gaussian distribution results from mapping a normally distributed random variable $\gamma \sim N(\mu, \sigma^2)$, $\varphi \in \mathbb{R}^1$ onto a unit circle $\mathbb{S}^1$:

$$\varphi = \psi(\gamma) = \text{mod}(\gamma + \pi, 2\pi) - \pi \qquad (A.2)$$

The filtering algorithm is summarised below:

**Predict :**

$$\hat{\mathbf{z}}_{t}^{-} = A\hat{\mathbf{z}}_{t-1}^{-}$$

$$\hat{\mathbf{z}}_{t}^{-}[1] = \psi(\hat{\mathbf{z}}_{t-1}^{-}[1])$$

$$\hat{\Sigma}_{t}^{-} = A\hat{\Sigma}_{t-1}^{-}A^{T} + \Sigma_{v}$$

**Correct :**

$$K_t = \frac{\hat{\Sigma}_t^{-} B^T}{B \hat{\Sigma}_t^{-} B^T + \sigma_w^2}$$

$$g_t = \sum_{l=-1}^{1} ((\varphi_t + 2\pi l) - \hat{\mathbf{z}}_t^{-}[1])\eta_{t,l}$$

$$\hat{\mathbf{z}}_t = \hat{\mathbf{z}}_t^{-} + K_t g_t$$

$$\hat{\Sigma}_t = (I - K_t B) \hat{\Sigma}_t^{-} \qquad (A.3)$$

where $\mathbf{z}_t = \begin{bmatrix} \varphi_t \\ \dot{\varphi}_t \end{bmatrix}$ is the system's state, $\mathbf{z}_t[1] = \varphi_t$, $A = \begin{bmatrix} 1 & dt \\ 0 & 1 \end{bmatrix}$

is the linearised state transition matrix, $dt$ is the time interval between the $(t-1)^{th}$ and $(t)^{th}$ epochs, $\Sigma_v$ is the state covariance matrix, $\Sigma_t$ is the state covariance matrix estimate, $B = [1 \ 0]$ is the observation matrix, $\sigma_w$ is the observation variance, $I$ is an identity matrix, and $\eta_{t,l} = N(\varphi_t + 2\pi l | \hat{\mathbf{z}}_t^{-}[1], \sigma_w^2)/ \sum_{m=-\infty}^{\infty} N(\varphi_t + 2\pi m | \hat{\mathbf{z}}_t^{-}[1], \sigma_w^2)$ represents the probability of a replicate.

First, the filter is run on the phase data turned “backwards” in time running from $t_N$ to $t_1$, where $N$ is the number of data points. The system state $\mathbf{z}$ is initialised as

$$\mathbf{z}_0 = \begin{bmatrix} \varphi_N \\ 0 \end{bmatrix} \qquad (A.4)$$

The output of this filtering procedure at $t_0$ is used as the initial condition for running the Kalman filter “forward” in time. Usually, it is enough to run the filter backward and forward once to get a robust and reliable result, but if several iterations are needed, the output of the forward-run filter at $t_N$ is used to update the initial condition for the backward-run.