The Potsdam Open Source Radio Interferometry Tool (PORT)

The fundamental input data sources feeding the PORT processing chain (level 2 data analysis) are time series of group delays stored in the vgosDB ¹⁴ or legacy NGS ¹⁵ data formats. In addition to the delay information, meteorological data, such as temperature, pressure, and relative humidity, provided by sensors at the antenna sites, as well as hardware-related calibration data are included in the vgosDB/NGS file(s). Furthermore, PORT collects a priori EOP time series from the International Earth Rotation and Reference Systems Service (IERS) data archives and calculates reference point displacements at a station site. Atmospheric refraction causes delay and bending, which are considered in PORT as well. Mapping functions and atmospheric parameter models consider both effects at the same time. The a priori zenith delays are calculated using the Saastamoinen refraction model (Petit & Luzum 2010) inserting meteorological data from either in situ observations or from meteorological analysis fields of ECMWF. The a priori horizontal atmospheric gradients are calculated from the APG model (Böhm et al. 2013). Finally, general and session-specific processing options and runtime control variables are initialized from ASCII text files using a PORT-specific data format. If information on outlier data points, stored in a separate ASCII text file, is available, the corresponding observations are flagged and excluded from the subsequent analysis tasks.

Besides the functional code for VLBI data analysis, PORT includes code for operational, administrative and format conversion tasks. Monitoring of the IVS archives for the publication of novel or revised vgosDB/NGS data files is handled by a set of Python scripts. The scripts are intended to run regularly at predefined intervals. If new observations are found, the operator on duty is alerted by e-mail and optionally in terms of SMS messages to analyze the session. The operational processing toolkit of PORT includes a database which holds information on the analysis status of each individual VLBI session.

Table 1 lists the most relevant components of the functional model. The calculation takes into account ITRF2014 (Altamimi et al. 2016) station and ICRF3 (Charlot et al. 2020) source coordinates (other reference frames may be chosen as well), loading models, atmospheric mapping functions (Dill & Dobslaw 2013; Balidakis et al. 2018), and relativistic corrections according to the consensus model (Eubanks 1991; Titov et al. 2011) following the Conventions by the IERS (Petit & Luzum 2010). The JPL solar system ephemerides listed in Table 1 are supported out of the box; ephemerides format-compatible with the JPL ephemerides may be used as well. The observables are assumed to be uncorrelated and semi-empirically determined weights are applied.

In addition, the design matrix entries, consisting of the partial derivatives of the observables with respect to each parameter, are calculated. The derivatives are taken at the parameters' a priori values, since in geodetic VLBI the adjustments are assumed to be less than about 4 mas for radio source coordinates and less than about 2 m for station coordinates and thus the linear approximation of the functional model suffices for the desired accuracy (Eubanks 1991).

Before starting with the parameter estimation the user is advised to check the correlator report and auxiliary station-specific information. PORT supports these preprocessing tasks, first, by downloading the correlator report (if it exists) and by storing it in the PORT data directory and, second, by providing plots of meteorological data from sensors at the individual stations as well as phase calibration information recorded during the 24 hr measurement period.

Figures 2 and 3 show examples of these auxiliary data recorded at station ONSALA60 during session 14MAY10XA. Together with the correlator report, these plots are supposed to guide the analyst in his/her decisions during the subsequent analysis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The parameter estimation process starts by determining the difference between the observed and the modeled delay values—the reduced observation vector. All parameters are realized as linear splines (piecewise linear functions). The reference epoch is the 0:00 UTC within the session (24 hr sessions) or closest to the observation epochs (CONT sessions, intensive sessions) in order to be compatible with other space geodetic techniques parameter epochs, such as the ones from GNSS. The length of the spline segments is constant inside a session for each parameter and can be chosen by the analyst. Depending on the chosen segment length several parameters will be outputted by PORT until the first and the last epoch of the observations are covered within a parameter segment. Occasionally, stations do not observe for some time during a session owing to other commitments, for example, the participation in an Intensive session, unforgiving weather, or technical issues. To avoid the noise increase associated with setting up unknowns that are not supported by observations, PORT is capable of refraining from doing so.

For each EOP two parameters are estimated as part of the normal equation system, one at the epoch of the first scan, the other at the epoch of the last scan of the session. Thus, PORT allows to output EOP offsets as well as their rates. Optionally, EOP estimation at midnight epochs may be set up to facilitate the comparison with other space geodetic techniques. For clock parameters, 1 hourly segments are chosen. Accordingly, PORT outputs 25 parameters for (n_clk − 1) clocks in the network. One value at 0:00 between day 1 and day 2 and then forward and backward in time every integer hour until no more observation epochs are outside. The associated constraints are added in terms of pseudo-observations either by constraining a parameter to its a priori value with a certain weight (absolute constraint) or by constraining the rate of change of a parameter with a certain weight (relative constraint). If just one adjusted parameter for the entire session is targeted, tight absolute constraints can be applied. Depending on the type of the parameter and its anticipated fluctuation per time, the segment lengths listed in Table 2 are typically applied. In addition, the table lists the corresponding absolute and relative constraints, which have been obtained empirically. For ITRF2020, a set of specific piecewise linear functions intervals for clock and atmospheric delay estimation was used for each session, depending on the observation geometry. For stations with poor observation geometry, fewer knots were set up.

At present, efforts are being made to incorporate VieVS@GFZ's Kalman Filter code (Nilsson et al. 2015; Soja et al. 2015, 2016; Soja 2016) into PORT as well. In the current PORT version the parameters are estimated with least squares adjustment (Koch 1999; Niemeier 2001; Teke et al. 2012). The adjustment is performed in two iterations. The parameters estimated in the first solution are clock polynomials of arbitrary degree (usually not higher than second) for each station—except for one station. The coefficients of this reference station are defined to be zero and, hence, this clock serves as the reference for the other clocks in the network. The clock adjustments may exceed hundreds of microseconds, in extreme cases even one second. In addition, harmonic clock variations, baseline clock offsets and one zenith delay parameter per station are estimated, their contributions are removed and the final adjustment estimates the remaining parameters including again clock offsets, drifts and drift rates.

In Table 2 the parameters estimated by PORT during single session analysis are listed. Apart from EOP, station positions, radio source coordinates, and tropospheric delay parameters (Balidakis et al. 2018), station clock offsets and drifts may be determined. In addition, PORT estimates baseline clock offsets, constant delays on specific baselines, which absorb instrumental effects (e.g., due to failure of individual frequency channels because of radio-frequency interference) on the corresponding baseline. The complete mathematical model of VLBI data analysis is for example described in Heinkelmann (2008).

The estimation results are presented to the analyst in several ways. First, following the approach familiar from the VieVS software package (Böhm et al. 2012, 2018), the residuals v

$$\begin{eqnarray}&&{\boldsymbol{v}}\equiv {\boldsymbol{A}}\cdot {\boldsymbol{x}}-{\boldsymbol{l}}\end{eqnarray} \tag{ 1 }$$

as well as their uncertainties are displayed as time series on screen for visual inspection. Here, A denotes the design matrix, x the adjustments for the unknown parameters, and l the reduced observation vector. The user may choose to restrict the plotted data by selecting individual baselines, stations, or radio sources.

Figures 4 and 5 illustrate the interactive procedure. Initially, the analyst can check the residuals from the first estimation stage (see Section 3.4) for large clock breaks. Representing residuals station-wise allows the user to identify systematic errors related to hardware issues at the selected station. For example, the time series of residuals obtained from the analysis of session 18OCT01XA at station WETTZELL, shown in Figure 4, exhibits a discontinuity of about 1.7 × 10⁹ cm, corresponding to a clock jump of about 55 ms at 58393.6294 MJD (2018 October 2 15:06:20 UTC). Using the mouse pointer a marker (red line) is placed at that time to insert a clock break. Every time a clock break is added at a given station, the coefficients of the affected clock polynomial are estimated anew. After repeating the estimation cycle, the residuals from the final estimation stage are inspected and outlier events are flagged, see Figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Second, to support the analyst in optimizing station-, baseline-, and/or radio source-based weights or even to removing a particular station altogether, a table detailing the goodness of the fit (χ²), the weighted root mean square (rms) residual deviations, the weighted mean values, as well as the (effective) number of observations per station, baseline, and radio source is produced and printed on screen.

Figure 6 displays an example for this type of table obtained as a result of processing session 06MAR23XE using PORT's current runtime option and processing control settings, denoted as data analysis version "2020b." These settings include information on, among others, selected models (see Table 1), outlier removal methods, elevation cut-off angles, vgosDB quality control flags, threshold values of the minimum number of stations, sources and/or baselines required, parameterizations of the station clock errors, application of constraints (pseudo-observations), and additional uncertainties to be added to the a priori measurement uncertainties. In the first part the station- and baseline-specific statistics are listed, including the goodness of the fit, in the following denoted "chi-squared" (and also "chi∧2" in Figure 6),

$$\begin{eqnarray}&&{\chi }^{2}\equiv \displaystyle \frac{{{\boldsymbol{v}}}^{T}\cdot {\boldsymbol{P}}\cdot {\boldsymbol{v}}}{n-{n}_{\mathrm{par}}},\end{eqnarray} \tag{ 2 }$$

the weighted mean value of the residuals

$$\begin{eqnarray}&&\mathrm{WM}\equiv \langle {\boldsymbol{v}}\rangle \equiv \displaystyle \frac{{{\boldsymbol{e}}}^{T}\cdot {\boldsymbol{P}}\cdot {\boldsymbol{v}}}{{{\boldsymbol{e}}}^{T}\cdot {\boldsymbol{P}}\cdot {\boldsymbol{e}}},\end{eqnarray} \tag{ 3 }$$

the weighted root mean square of the residuals

$$\begin{eqnarray}&&\mathrm{WRMS}\equiv \sqrt{\displaystyle \frac{{{\boldsymbol{v}}}^{T}\cdot {\boldsymbol{P}}\cdot {\boldsymbol{v}}}{{{\boldsymbol{e}}}^{T}\cdot {\boldsymbol{P}}\cdot {\boldsymbol{e}}}},\end{eqnarray} \tag{ 4 }$$

the weighted standard deviation of the residuals

$$\begin{eqnarray}&&\mathrm{WSTD}\equiv \sqrt{\displaystyle \frac{{({\boldsymbol{v}}-\langle {\boldsymbol{v}}\rangle )}^{T}\cdot {\boldsymbol{P}}\cdot ({\boldsymbol{v}}-\langle {\boldsymbol{v}}\rangle )}{{{\boldsymbol{e}}}^{T}\cdot {\boldsymbol{P}}\cdot {\boldsymbol{e}}}},\end{eqnarray} \tag{ 5 }$$

the effective numbers of observations

$$\begin{eqnarray}&&{n}_{\mathrm{eff}}\equiv \displaystyle \frac{{{\boldsymbol{e}}}^{T}\cdot {\boldsymbol{P}}\cdot {\boldsymbol{e}}}{\max \left({\boldsymbol{P}}\right)}\end{eqnarray} \tag{ 6 }$$

denoted "n_eff" in Figure 6, and the Student t-distribution (Student 1908), in Figure 6 listed as "t_cdf." Here, n denotes the actual number of observations (listed as "N" in Figure 6), n_par the number of parameters, P the weight matrix for the observations, e ≡ (1, 1,...,1)^T , and $\max \left({\boldsymbol{P}}\right)$ is the largest element of P . These statistics guide the analyst in the decision process which station to down- or up-weight and/or whether baseline clock offsets need to be admitted. In the second part of the table the resulting adjustments for station positions, radio source positions, and Earth orientation parameters are listed.

Zoom In Zoom Out Reset image size

On the basis of the information described above, the analyst can enter or modify corresponding weights, activate baseline offsets, and remove/deselect stations and/or radio sources. A weighted Student t-test triggers a warning message if mean observable residuals at individual stations, baselines, and/or radio sources significantly differ from the value expected on the basis of the null hypothesis and exhibit statistical anomalies (see, e.g., Bevington 1992; Press et al. 1992). Based on this information the estimation/editing/review procedure is iterated until acceptable overall goodness of fit (typically, 0.1 < χ² < 10) and weighted residual rms values are achieved.

Besides the interactive mode discussed in Sections 3.3 and 3.5, PORT supports a non-interactive batch mode. Sessions without (or with already known) clock breaks may be processed in this unsupervised automatic mode where outliers are identified using Baarda's data snooping method, an outlier detection based on statistical tests (Baarda 1967, 1968; Rofatto et al. 2020). For clarity it is emphasized that in automatic mode—just as in interactive mode—PORT estimates parameters and assembles normal equations as well. The speedy analysis of large number of sessions enabled by automatic mode allows for efficient tests and intercomparisons of diverse models (Table 1).

Once the analyst has reviewed the residuals and verified their statistics, the estimated parameters, the a priori information employed during the estimation procedure, as well as the normal equation systems prepared for the global least squares adjustment are exported, stored in the Solution INdependent EXchange (SINEX) format. ¹⁶ Finally, the user has the option to upload the SINEX file to the IVS data servers using a provided Python script. In addition, these data are written to disk in PORT-specific binary file format for subsequent processing. Another set of output files is generated from the selected options and outlier. ¹⁷

Verification and validation of the estimated parameters are an important aspect of geodetic data analysis and one of the key PORT development objectives is to provide the analyst with suitable tools to visualize the results. Furthermore, PORT supports the intercomparison with the analysis results provided by other IVS analysis centers; its quality control tool set allows for the creation of graphical representations of deviations between two or more SINEX-formatted solutions. The comparisons include overall statistics such as goodness of fits and weighted residual rms, and subset statistics of individual EOP variables, stations, or radio sources. For ease of access the images are displayed through a set of static HTML pages.

Figures 7 and 8 exemplarily illustrate statistical results for GFZ's VLBI contribution to ITRF2020 ¹⁸ (K. Balidakis et al. 2021, in preparation). Figure 7 displays the temporal evolution of weighted rms corresponding to the PORT analysis of more than 6500 sessions recorded between 1979 August and 2020 December. The histogram of the deviations between the absolute value of UT1−UTC extracted from PORT's "gfz2020b" and NASA Goddard's "gsf2019a" solutions are plotted in Figure 8. While perfect agreement has not yet been achieved, the two data sets appear to be consistent on the 30–100 μs level.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size