Ionospheric Modelling using GPS to Calibrate the MWA. II: Regional Ionospheric Modelling using GPS and GLONASS to Estimate Ionospheric Gradients

Abstract

We estimate spatial gradients in the ionosphere using the Global Positioning System and GLONASS (Russian global navigation system) observations, utilising data from multiple Global Positioning System stations in the vicinity of Murchison Radio-astronomy Observatory. In previous work, the ionosphere was characterised using a single-station to model the ionosphere as a single layer of fixed height and this was compared with ionospheric data derived from radio astronomy observations obtained from the Murchison Widefield Array. Having made improvements to our data quality (via cycle slip detection and repair) and incorporating data from the GLONASS system, we now present a multi-station approach. These two developments significantly improve our modelling of the ionosphere. We also explore the effects of a variable-height model. We conclude that modelling the small-scale features in the ionosphere that have been observed with the MWA will require a much denser network of Global Navigation Satellite System stations than is currently available at the Murchison Radio-astronomy Observatory.

1 INTRODUCTION

The Earth’s ionosphere has a significant effect upon radio astronomy observations, in particular at low radio frequencies. Below a critical frequency, the ionosphere is opaque (Rawer 1993; Wilson, Rohlfs, & Huettemeister, 2014) and can radiate (Davies, 1990). Refractive effects manifest themselves as apparent position shifts of celestial radio sources (Wilson et al., 2014). The ionosphere is dispersive (Davies, 1990) and causes Faraday Rotation of radio waves (Wilson et al., 2014). The ionosphere can also cause diffractive effects (Davies 1990).

The ionosphere is of great interest as a target for research for many reasons and has been the subject of detailed study for decades. See Davies (1990) and Rawer (1993) for general reviews of the ionosphere and Thompson, Moran, & Swenson Jr (2008) and Wilson et al. (2014), for example, for the connection between ionospheric research and radio astronomy.

With a new generation of wide-field low radio frequency telescopes now in operation, including the Murchison Widefield Array (MWA) (Tingay et al. 2013), Low-Frequency Array (LOFAR) (van Haarlem et al. 2013), Precision Array for Probing the Epoch of Reionisation (PAPER) (Parsons et al. 2010), and the Long Wavelength Array (LWA) (Ellingson et al. 2009), interest in the effect of the ionosphere in radio astronomy is greatly renewed. This is for two reasons: first, the new generation of radio telescopes have the ability to probe the ionosphere in unprecedented detail (Loi et al. 2015b). Second, because as the new generation of radio telescopes are designed and built, calibration of the effects of the ionosphere become more challenging and characterising its effects radio astronomy observations is critical.

The MWA, the low frequency precursor for the Square Kilometre Array (SKA) located in Western Australia, has 128 aperture array elements (called tiles), has the maximum baseline length of ~ 3 km, and an extreme wide field-of-view (FoV) capability (25° full width at half maximum at 150 MHz) (Tingay et al. 2013). These characteristics place the MWA in a regime where different paths through the ionosphere are observed for different sources across the FoV, but the effect of the ionosphere on each array element can be assumed to be the same (Lonsdale et al. 2009). However, future instruments (e.g. the low frequency component of the SKA Hall 2005) will have much longer baselines, necessitating not only a direction-dependent calibration, but also a different solution for each interferometer element, a far more difficult and computational intensive problem to solve.

This motivates us to look at Global Satellite Navigation Systems (GNSS) as a possible source of information on the ionosphere, not only as a direct source of information for calibration (perhaps a low-resolution model that can reduce the parameter space to be searched), but also both for climatology (understanding the range of ionospheric conditions in a statistical sense), and for identifying whether conditions prevailing during a particular observations were favourable for radio astronomy (without taking the much more route of determining this from the radio telescope data itself). In previous work (Arora et al. 2015), we undertook an initial study of refractive effects due to the ionosphere, as observed by the MWA, and compared them with independent measurements using the Global Positioning System (GPS). Bulk ionospheric gradients causing the refractive effects observed with the MWA were found to agree well with those estimated from GPS observables. The results presented in Arora et al. (2015) establish a methodology and show that ionospheric information can plausibly be obtained from GPS, to help calibrate the MWA.

The research presented here aims to build on our earlier work. Before, ionospheric modelling was performed by using data from a single GPS station for any given ionospheric solution. To capture the ionospheric behaviour on finer spatial scales, additional data is required. To this end, we incorporate the GLONASS satellite system into our analysis. GLONASS currently has 24 active satellites in orbit. Further, we now upgrade our ‘single-station’ analysis to a ‘multi-station’ analysis, whereby each ionospheric solution is calculated using data from multiple receiving stations.

Finally, we explore the effectiveness of various relaxations of the single-layer model for ionospheric modelling. Methods to include spatial and temporal variations into the height of the single-layer model are discussed.

This paper is organised as follows: Section 2 presents the GNSS data pre-processing methodology. In Section 3, the GLONASS system overview and a combined GPS and GLONASS observation model are presented. Further, the effect of the single-layer model height on the estimated ionosphere coefficients and methods to incorporate the variation in single-layer model height are presented. The multi-station approach to estimate ionosphere coefficients using GPS and GLONASS is presented in Section 3. Section 4 presents the summary of obtaining ionosphere gradients from MWA observations as a function of position shifts. The results from the multi-station approach are presented and discussed in Section 5. The paper is concluded in Section 6, where we discuss future directions for this work.

2 GNSS DATA PRE-PROCESSING

In GNSS data pre-processing, discontinuities in the phase and code observables are identified and repaired (Lichtenegger & Hofmann-Wellenhof 1990; Remondi 1985, among others). The uncertainties in the GNSS observables are phase cycle-slips and jumps, and multi-path effects. In our previous work (Arora et al. 2015), pre-processing of the GNSS observables was applied; however, it is discussed in detail here for the first time.

2.1. Cycle-slip detection and repair

When a receiver tracks a satellite, the integer and fractional number of cycles (total number of wavelengths at the GPS frequency) between the receiver and the satellite are recorded as phase observables. However, the initial number of phase cycles, upon first acquisition of the satellite signal, remain ambiguous. The ambiguities present in the phase observables remain constant for a complete satellite pass, unless a cycle-slip occurs. Cycle-slip can occur for a number of reasons including, temporary blocking of the GNSS signal by a physical obstruction, multi-path effects, high ionospheric activity, and low signal-to-noise ratio (Hofmann, Lichtenegger, & Wasle 2008). It is important to account for cycle-slips in order to ensure the continuity of carrier phase data on which high precision GNSS applications are dependent.

There are three stages for pre-processing of cycle-slips. First, the cycle-slip is detected, second, its magnitude is quantified, and third, it is flagged or accounted for in the observables. The generic approach to cycle-slip detection is by forming linear combinations of observables. During pre-processing using the BERNESE software (Dach et al. 2007), a combination of the phase and code observables is formed, also known as Melbourne–Wübbena combination (Melbourne 1985; Wübbena 1985), which allows detection of cycle-slips. However, the noise of the observable in Melbourne–Wübbena combination is driven by the noise of the code observable, the code observables are found to have a precision of 25 cm or worse. Other combinations, namely, the ‘wide-lane’ (Hofmann et al. 2008) and ‘ionosphere-free’ (Hofmann et al. 2008) combinations, although driven by the very precise phase observables, have noise of 5.7 and 3.0 times the original observables, respectively (Dach et al. 2007).

The Geometry-Free combination can also be used to detect cycle-slips (Vaclavovic & Dousa 2015), the Geometry-Free phase observable, used to detect cycle-slips, is denoted as L4. The noise of the L4 observable is $\sqrt{2}$ times the precision of the phase observables, lower than all of the earlier mentioned linear combinations. In our approach, cycle-slips are detected by using Geometry-Free combination of observables.

For a Geometry-Free combination, the new phase observables (L4) has the constant instrumental term, which has ambiguities and other biases, and the ionospheric error. The time difference of the L4 observables, L4(t) − L4(t − 1), can be used to eliminate all other terms except the variable part of the ionospheric error. Any unusual variation in the ionosphere can be easily flagged for a cycle-slip. Following Vaclavovic & Dousa (2015), the expression for cycle-slip detection, using Geometry-Free phase observables, is given as follows:

(1) $$\begin{equation} | {\rm L4}(t) - {\rm L4}(t-1)| > k \cdot \sigma _{{\rm L4}} + \Delta I_{\text{max}}, \end{equation}$$

where L4 is the Geometry-Free observable formed from GNSS phase observables L1 and L2, t is the observation time, k is a scaling factor, $\sigma _{{\rm L4}}$ is the precision of the ${\rm L4}$ observable, and $I_{\text{max}}$ is the maximal ionospheric delay. Following Vaclavovic & Dousa (2015), $\Delta I_{\text{max}}$ is chosen to be 0.4 $\text{m h}^{-1}$ and the factor k as 4.

Once the cycle-slip is detected, the hypothesis can further be tested by using the Geometry-Free code observable (P4) for a sufficient number of epochs (Teunissen & Kleusberg 1998).

The cycle-slip can be repaired by estimating the time propagation of the ionosphere from L4 observables, using the information before and after the slip. The L4 observables being precise, are capable of sensing ionospheric variations as small as ~ 0.04 Total Electron Content Unit (TECU) (at GPS frequencies, 1 TECU = 10¹⁶ electrons m⁻²). Considering the location of GNSS receivers we are using (far below the equatorial anomaly), extreme ionospheric variations are not expected. However, ionospheric variations of the order of 1 TECU every 30 s can be easily accounted for with this algorithm.

The above algorithm, however, is not capable of differentiating whether the slip occurred on L1 or L2 phase observable. Since our software makes use of L4 observables for estimating the ionosphere and other unknowns, estimating cycle-slips at individual frequencies is not a necessary.

3 IONOSPHERIC MODELLING USING GPS AND GLONASS

3.1. GLONASS system overview

The GLONASS system, currently has 24 operational satellites in its constellation, which continuously transmit dual frequency data centred at frequencies L1₀ = 1602.0 MHz and L2₀ = 1246.0 MHz. Each GLONASS satellite transmits on a different frequency using a 15-channel Frequency Division Multiple Access (FDMA) technique. The frequency for each channel is given by L1 _n = L1₀ + n × (9/16) MHz and L2 _n = L2₀ + n × (7/16) MHz, where n is the GLONASS channel number, n = −7, . . ., 0, . . ., 6. The GLONASS channel number for each satellite can be obtained from the GLONASS broadcast ephemerides file.

The GPS and GLONASS systems transmit in different reference times, which needs to be compensated in the code observables, while realising a common time system for processing the data. However, if the Geometry-Free combination is used, as in our work, the code observable correction is compensated and not required.

3.2. GPS and GLONASS observation model

The Geometry-Free GPS observation model is discussed in detail in Arora et al. (2015). We recall it below and append the GLONASS observation model as follows:

(2) $$\begin{equation} \displaystyle E(\Phi _{r,21}^{Gs}) = \displaystyle \Phi _{r,1}^{Gs}-\Phi _{r,2}^{Gs} = - \iota _{r,21}^{Gs} + \mathrm{C}_{r}^{Gs}, \end{equation}$$

(3) $$\begin{equation} \displaystyle E(P_{r,21}^{Gs}) = \displaystyle P_{r,1}^{Gs}-P_{r,2}^{Gs} = \iota _{r,21}^{Gs} + c \cdot (d_{r,21} - d^{Gs}_{,21}), \end{equation}$$

(4) $$\begin{equation} \displaystyle E(\Phi _{r,21}^{Rs}) = \displaystyle \Phi _{r,1}^{Rs}-\Phi _{r,2}^{Rs} = - \frac{\mu _{21}^{R}}{\mu _{21}} \iota _{r}^{s} + \mathrm{C}_{r}^{Rs}, \end{equation}$$

(5) $$\begin{equation} \displaystyle E(P_{r,21}^{Rs}) = \displaystyle P_{r,1}^{Rs}-P_{r,2}^{Rs} = \frac{\mu _{21}^{R}}{\mu _{21}} \iota _{r}^{s} + c \cdot d^{Rs}_{r,21}, \end{equation}$$

where E(·) is the expectation operator, Φ is the phase observable, P is the code observable, subscript _r indicates receiver, ₁, ₂, and ₂₁ indicates GNSS frequency/frequency combinations corresponding to phase (or code) observables, L1 (or C1), L2 (or P2), and L4 (or P4), respectively. Superscripts ^Gs , ^Rs indicate GPS and GLONASS satellites, respectively, c · (d _{r, 21}) and c(d^Gs _{, 21}) are the GPS receiver and satellite differential code biases (DCBs), respectively, μ₂₁ is the GPS frequency coefficient given as, μ₂₁ = μ₁ − μ₂ and $\mu _{1} = \frac{1}{f_{1}^{2}}$ , $\mu _{2} = \frac{1}{f_{2}^{2}}$ , f ₁ and f ₂ are GPS frequencies at L1 and L2. Similarly, the GLONASS frequency coefficient is given by μ ^R ₂₁.

The instrumental biases and other unknowns are estimated for each GNSS receiver using the method of least squares with a Kalman filter, as described in Arora et al. (2015). The precision of the time-constant parameters propagate as the inverse of the square-root of the number of epochs (n), $\sigma = 1/\sqrt{n}$ . For a continuous satellite arc, n is very large, ~ 100 or more. This results in a very precise estimation of time-constant parameters.

The line-of-sight Total Electron Content (TEC) between receiver and the satellite, is also known as Slant TEC (STEC). STEC can be retrieved from L4 phase observables (Φ _{r, 21}). In our work, we retrieve the STEC for both GPS and GLONASS satellites, by substituting for the time-constant parameters for L4 observables, estimated using the single-station approach. Figure 1 presents the retrieved STEC for MRO1 and MEDO Geoscience Australia (GA) GNSS stations on DOY 062, refer Figure 6 and Table 1 for details of all the stations.

Figure 1. Retrieved STEC for the MRO1 and MEDO Geoscience Australia (GA) GNSS stations on DOY 062, year 2014. (a) STEC for MRO1—GPS only. (b) STEC for MEDO—GPS+GLONASS.

Table 1. Description of the selected GA GPS/GNSS stations and the MWA. Data were available for all the four observing sessions, DOY 062, 063, 065, and 075, year 2014. The acronyms given under GNSS, G, and GR stand for GPS only and GPS+GLONASS, respectively.

^a Partial data available, from 00:00:00 UTC to 18:07:00 UTC on DOY 062, year 2014.

3.3. Single-station versus multi-station approach

In a single-station approach, the ionospheric coefficients and time constant parameters, namely the frequency dependent receiver and satellite biases on code observable, also known as DCBs, and the constant phase term (containing the ambiguities and other biases) are estimated using the observables from a single-station. Due to the limited number of observations, the ionospheric coefficients are estimated at an interval of 10 min to allow sufficient robustness. The ionosphere gradients show artificially high levels of variation as a satellite sets and rises, again due to the limited number of observations used.

For multi-station ionosphere modelling, the GNSS model can be designed such that all the parameters (ionosphere and other time-constant biases) are estimated in a multi-station mode. This approach adds constraints to the time-constant satellite-specific bias parameters, namely, the GPS satellite DCBs. However, the multi-station approach must account for different number of satellites being visible, for different receivers at any given time. Another feasible approach is to consider only satellites that are visible to all receivers at a given time; however, this can result in loss of information. In this work, the multi-station modelling is performed using the retrieved STEC for each GNSS receiver, which eliminates the need to estimate any time constant parameter.

In our work, the retrieved STEC is derived from the precise phase observables only, in contrast to the generic approach of using phase-smoothed-code observables (Gao et al. 2002; Chevalier et al. 2013). In the phase-smoothed-code approach, the phase as well as the code observables are used to retrieve the STEC, by substituting for the receiver and satellite DCBs. Also that, the noise of retrieved STEC from the phase-smoothed-code approach is driven by the noise of the code observable. By using the phase observables to retrieve STEC, the noise of the STEC is driven by the noise of the phase observable. In our approach, the time constant biases in the phase observable, specific to each receiver and satellite are estimated during the single-station approach, with sufficient precision that complements the noise of the phase observable. These constant phase biases are then subtracted to retrieve the STEC for each receiver-satellite pair. By using retrieved STEC, the ionospheric coefficients can be estimated at a higher temporal resolution, in our work, the ionospheric coefficients are estimated every 2 min.

3.4. Effective height of the ionospheric layer

For ionospheric modelling using a single-layer model, the ionospheric electron density is assumed to be concentrated at a fixed height, H _ion. H _ion, is used to compute the obliquity factor (mapping function) and the coordinates of the Ionospheric Pierce Point (IPP).

To understand the effect of H _ion on estimated ionospheric coefficients and receiver DCBs, H _ion was varied between 350 and 550 kms in steps of 50 km. GPS observables for GA station MRO1 for DOY 062, year 2014 were used for this analysis. GPS data were processed using a single-station, single-layer ionospheric modelling approach. Figure 2(a) presents the estimated values of VTEC (Vertical TEC) for different values of H _ion. The differences in VTEC lie between ~ 0.5 and ~ 1 TECU at different times during the day. This effect is absorbed by the receiver DCBs, the receiver DCBs are affected by an amount corresponding to the maximum constant difference in VTEC over 24 h, that is ~ 0.5 TECU [Figure 2(b)]. Hence, selection of the value of H _ion plays a significant role in ionosphere modelling and has an effect on the estimated receiver DCBs. The ionosphere gradients, however, do not seem to be significantly affected when the height is varied by a constant value, refer Figure 4. It is important to incorporate the spatial variations, if they are significant, to the ionospheric single-layer height in order to observe any significant change in the gradients.

Figure 2. Effect on estimated VTEC and receiver DCBs by the choice of H _ion, H _ion is varied between 350 to 550 km in steps of 50 km. Note the average precision of VTEC is ~ 0.03 TECU. (a) VTEC as a function of H _ion. (b) Receiver DCB as a function of H _ion.

A more realistic representation of the single-layer height would be to account for the effective height of the ionospheric layer, H _eff. H _eff is the height at which the electron density reaches its median value, and is a function of location and time. H _eff can be deduced from the ionospheric profiles. The ionospheric profiles can be obtained from empirical models like the International Reference Ionosphere (IRI) (Bilitza et al. 2014). However, for IRI model, the maximum height for the ionospheric profiles are limited to 2000 km. Ionospheric–Plasmaspheric models like the Parameterised ionospheric model (PIM) (Daniell et al. 1995) and extension of IRI to plasmasphere (IRI-Plas) (Gulyaeva & Bilitza 2012; Gulyaeva et al. 2013) model the ionosphere up to plasmaspheric altitudes (above 20 000 km). In our work, we make use of IRI-Plas model, available as an external software, to generate ionospheric profiles.

Though H _eff can be deduced from the ionospheric profiles, there is no direct source of information of this parameter, regarding its spatial variation. H _eff can also be related to hmF2, hmF2 is the height at which maximum ionisation is reached, which lies in the F2 region. The spatial variation of hmF2 can be obtained from the global hmF2 maps generated by IRI-Plas (refer Figure 3), available for download ¹ . The temporal and spatial resolution of hmF2 global maps is 1 h and 5°/2.5° in longitude and latitude, respectively.

Figure 3. HmF2 global map from IRI-Plas model for DOY 062, 2014.

Figure 4. Effect on estimated ionosphere gradients by the choice of H _ion, H _ion is varied between 350 to 550 km in steps of 50 km. (a) EW gradient v/s HION. (b) NS gradient v/s HION.

The effective height, H _eff, though is different from hmF2, however can be related to hmF2. Ionospheric profiles generated using IRI-Plas model show that H _eff is greater than hmF2, refer Figure 5. It can be noted from Figure 5, the difference of H _eff and hmF2 increases greatly during the night time, since the electron density in the F2 region decreases and hence H _eff increase significantly due to the plasmaspheric electron density. A discussion on H _eff and hmF2 can be found in Komjathy & Langley (1996). For generating global TEC maps, a variable height of the single layer model has been used for each of the ground stations (Komjathy et al., 1998). In our work, since we make use of regional GNSS network to determine ionospheric gradients, variable height for each receiver-satellite pair is computed.

Figure 5. Temporal variation of hmF2 and H _eff obtained from IRI-Plas ionosphere profiles at the Taylor series expansion point.

To make use of the effective height, H _eff, as an input to a single-layer model, a method for determining the H _eff for all the IPPs is required. For sufficiently small regional scales, as in our work, the offset between H _eff and hmF2 can be computed at the central location of the network. The only significant variation between H _eff and hmF2 is due to the solar time, which can be compensated for different IPPs within the network, as an argument of longitude of the IPP. In our work, a constant offset between H _eff and hmF2 is computed for the Taylor series expansion point. The Taylor series expansion point is the MWA look direction, and remains constant throughout the day (entire observation period). This offset is then applied to the hmF2 values corresponding to satellite IPPs, obtained from the hmF2 global maps, refer Figure 3.

3.5. GNSS multi-station modelling using retrieved STEC

The GA GNSS stations in the near vicinity of MRO are selected to perform regional modelling. The data from selected GA stations, namely, MRO1, WILU, MTMA, YAR3, MEDO, GASC, and TOMP were used for modelling. Figure 6 presents the location of the selected GNSS stations. The details of each GNSS station are given in Table 1.

Figure 6. Selected GNSS station locations from Geoscience Australia’s network (red), MWA location (blue), and MWA IPP (green) for the four MWA observation nights (DOY 062, 063, 065, and 075 marked by 1 to 4, respectively).

The retrieved STEC for all the GPS and GLONASS satellites are used for regional ionospheric modelling. For a single-layer model assumption, STEC can be related to the VTEC using a simple mapping function at a fixed height, refer Figure 7. The obliquity factor or the mapping function, F^s , is discussed in detail in our earlier paper (Arora et al. 2015), the geometry of the obliquity factor is illustrated in Figure 7. We briefly summarise the mapping function equation as follows:

(6) $$\begin{eqnarray} \left. \begin{array}{rcl}\displaystyle STEC &=& \displaystyle VTEC \cdot F^{s} \\ \\ \displaystyle F^{s} &=& \displaystyle \frac{1}{\cos (z^{\prime })} =\frac{1}{\sqrt{1-\sin ^{2} z^{\prime }}} \\ \\ \displaystyle \sin z^{\prime } &=& \displaystyle \frac{R_{e}}{R_{\text{e}}+H_{{\rm ion}}} \sin (z) \end{array}\right\rbrace, \end{eqnarray}$$

where z′ is the zenith angle at the IPP, $R_{\text{e}}$ is the mean radius of the Earth, considered to be 6 371 km and assuming a spherical Earth, H _ion is the height at the sub-ionospheric point, and z is the zenith angle of the satellite as seen by the receiver.

Figure 7. Ionosphere single-layer model representation.

Two models are considered for the single-layer height; first, assuming a fixed height of 450 km (H _ion); second, the spatially and temporally varying height inferred from the IRI-Plas model is used (H _eff).

The VTEC is further modelled using a Taylor series polynomial expansion, the expansion point being the MWA IPP. A second order polynomial function is used to model the VTEC, given as follows:

(7) $$\begin{eqnarray} VTEC (\varphi _{\text{m}},s) &=& VTEC_{0} + (\varphi _{\text{m}} - \varphi _{\text{m}_{0}})f^{\prime } \varphi + (s - s_{0})f^{\prime } s \nonumber\\ & & + \,(\varphi _{\text{m}} - \varphi _{\text{m}_{0}})^{2}f^{\prime \prime } \varphi _{\text{m}} \varphi _{\text{m}} + (s - s_{0})^{2}f^{\prime \prime } s s \nonumber \\ && +\, (\varphi _{\text{m}} - \varphi _{\text{m}_{0}}) (s - s_{0}) f^{\prime \prime } \varphi _{\text{m}} s. \end{eqnarray}$$

The Sun fixed longitude, s, is related to the local solar time (LT) as $s = \lambda _{\text{m}} + LT - \pi$ , where $\lambda _{\text{m}}$ is the geomagnetic longitude at IPP, LT is in radians, VTEC ₀ is the VTEC at the central location in the network, and $ f^{\prime } s, \ f^{\prime } \varphi _{\text{m}}, \ f^{\prime \prime } s s, \ f^{\prime \prime } \varphi _{\text{m}} \varphi _{\text{m}}, \ f^{\prime \prime } \varphi _{\text{m}} s$ are the first- and second-order derivatives of VTEC along the Sun fixed longitude and latitude, respectively.

Figure 8 presents a snapshot of the satellite IPPs for all the GA GNSS stations and MWA IPP locations.

Figure 8. A snapshot of MWA IPP (blue), GA GNSS stations (red), and satellite IPPs for 5 min (10 epochs) during MWA observations (gray) in earth-fixed reference frame, on DOY 062, 2014. The MWA IPP is considered for the Taylor series expansion point.

4 MWA IONOSPHERE

In order to validate our modelled ionosphere against radio astronomy data, precisely the same data were used as in our previous work. Please see Wayth et al. (2015) and Hurley-Walker et al. (2016, submitted) for a full description of the GaLactic and Extragalactic All-sky Murchison Widefield Array (GLEAM) survey, and Arora et al. (2015, Section 4 for a detailed description of how these MWA data may be compared with GPS data. For clarity, the essentials are presented again below.

The observations we have used, operate in drift-scan mode, where the telescope remains pointed at a single point on the meridian as radio sources drift past. The telescope also cycles through five frequency bands spanning the range 73–230 MHz, with each band observed for 2 min in every 10 min. Each 2-min observation is then processed to produce one image for each of four neighbouring sub-bands which make up the full instantaneous bandwidth of 30.72 MHz. Each of these images will typically contain many hundreds of radio sources, most of which are unresolved, and almost all of which are already known from previous surveys.

A catalogue of sources has been compiled which are expected to be bright and unresolved with the MWA (Harvey et al. in preparation). For the brightest 100 sources in each image, an elliptical Gaussian is fitted to a small subset of pixels corresponding to the a priori location of that source. This gives us the location of the source on the sky at each of four different frequencies. By fitting the change in location of the source as a function of λ², we can determine the magnitude and direction of ionospheric refraction in a way that is blind to errors in the a prior location of the source, or instrumental or imaging effects. The ionospheric shift is then averaged over all sources within the FoV.

The final result is therefore a time-series of ionosphere gradients, in both north–south (NS) and east–west (EW) directions, averaged over all sources within the MWA FoV. The centre of this FoV is taken to be the MWA pierce point. Only the lowest frequency band was used, giving a time resolution of 10 min.

5 RESULTS AND DISCUSSION

5.1. Comparisons of ionosphere gradients—single-station v/s multi-station approach

The gradients obtained from MWA observations in right ascension (EW direction) and declination (NS direction), were computed from the ionospheric coefficients given in equation (7). The ionospheric coefficients were estimated by considering the centre of MWA FoV as the expansion point. The gradients can further be calculated by considering a latitude/longitude separation, on scales similar to the MWA FoV, around the expansion point.The MWA would see a FoV of width ~ 200 km at an altitude of 450 km, and the MWA derived gradients represent the bulk shift over the entire FoV. The ionospheric coefficients obtained using GPS and GLONASS observations were used to compute the EW and NS gradients, over 1° of latitude/longitude (1° ≅ 100 km) either side of the expansion point. The variation of latitude/longitude separation was carried for 0.5° to 10° in order to understand its effect on the computed gradients, refer Figure 9. The gradients computed from GNSS, using the latitude/longitude separation considerations, that closely represent the spatial scales of the MWA FoV, that is 0.5° to 2°, though exhibit some variation, however, the effect is modest, refer Figure 9.

Figure 9. Effect on estimated ionosphere gradients by the choice of the latitude/longitude separation, the variation for the latitude/longitude was done between 0.5° to 10°. (a) EW gradient v/s longitudinal seperation. (b) NS gradient v/s latitudinal seperation.

The MWA derived gradients are plotted (in green) along GNSS observed gradients for comparison. The ionospheric gradients, in EW and NS directions, were estimated for 2014 March 03 (DOY 062), 2014 March 04 (DOY 063), 2014 March 06 (DOY 065), and 2014 March 16 (DOY 075) using both single-station (Arora et al. 2015) and multi-station approaches (refer to Figures 10 and 11). In each of the Figures 10 and 11, two different ionosphere gradients are estimated, first, considering a fixed height of the ionospheric layer (H _ion) at 450 km, plotted as blue line. Second, the height is assumed to vary in space and time, derived from IRI-Plas model (H _eff), indicated by red curve.

Figure 10. EW ionosphere gradients observed from GNSS data [blue(H _ion) and red(H _eff)] and the MWA (green) using single-station approach, GPS only (left column) and multi-station approach, GPS+GLONASS (right column) on DOY 062, 063, 065, and 075, year 2014. Note the average precision of EW gradients is ~ 0.07 × 10⁻⁵ and ~ 0.03 × 10⁻⁵ for single-station and multi-station approach, respectively. (a) Single-station—GPS only, DOY 062. (b) Multi-station—GPS+GLONASS, DOY 062. (c) Single-station—GPS only, DOY 063. (d) Multi-station—GPS+GLONASS, DOY 063. (e) Single-station—GPS only, DOY 065. (f) Multi-station—GPS+GLONASS, DOY 065. (g) Single-station—GPS only, DOY 075. (h) Multi-station—GPS+GLONASS, DOY 075.

Figure 11. NS ionosphere gradients observed from GNSS data (blue) and the MWA (green) using single-station approach, GPS only (left column) and multi-station approach, GPS+GLONASS (right column) on DOY 062, 063, 065, and 075, year 2014. Note the average precision of NS gradients is ~ 0.05 × 10⁻⁵ and ~ 0.03 × 10⁻⁵ for single-station and multi-station approach, respectively. (a) Single-station—GPS only, DOY 062. (b) Multi-station—GPS+GLONASS, DOY 062. (c) Single-station—GPS only, DOY 063. (d) Multi-station—GPS+GLONASS, DOY 063. (e) Single-station—GPS only, DOY 065. (f) Multi-station—GPS+GLONASS, DOY 065. (g) Single-station—GPS only, DOY 075. (h) Multi-station—GPS+GLONASS, DOY 075.

Table 2 presents the EW ( $r_{\text{EW}}$ ) and NS ( $r_{\text{NS}}$ ) gradient correlation between the GNSS and MWA observed gradients for the single-station and multi-station approach for all the days of observations.

Table 2. Correlation between the GNSS and MWA observed gradients in EW ( $r_{\text{EW}}$ ) and NS ( $r_{\text{NS}}$ ) components, its standard error ( $\sigma _{\text{r}}$ ) using single-station approach and multi-station approach.

The correlation between the GPS and MWA EW gradients are identical within the errors between the single-station and multi-station approaches for two of the days, refer Table 2. For the remaining two days, the EW gradient correlation was found to be significantly better using the multi-station approach. The NS gradient correlation was found to be significantly better for three of the four days using the multi-station approach. The general trend showed that the EW and NS gradients show better correlations with the MWA observed gradients when estimated using a multi-station approach rather than single-station approach.

The single-station approach is limited by the number of observations to constrain the gradients, hence the gradients appear to be noisy. This is however not the case with the multi-station approach. By using the model values for ionospheric shell height (H _eff) varying in space and time, a curvature is defined for the ionospheric layer, hence the gradients appear to have a steeper slope (Figures 10 and 11), as compared to fixed ionospheric shell height (H _ion).

6 CONCLUSIONS

In this work, we have explored several incremental improvements on our previous work (Arora et al. 2015) including: (1) the addition of GLONASS data to augment GPS data; (2) the development of a multi-station ionospheric solution rather than a single-station solution; and (3) the sensitivity of our analysis to varying height for a single-layer ionospheric model. This work is designed to explore the most effective future directions for the development of ionospheric modelling to support calibration of the MWA and other future instruments.

The height of the single-layer model is seen to play a significant role in the estimated ionosphere coefficients. The estimated ionosphere coefficients and receiver DCBs were seen to have a common minimum offset while the height of the single layer (H _ion) was varied. While a variable height of the single layer was incorporated using IRI-Plas model, the gradients appear to have a steeper slope. The increased steepness in the slope of the gradients could be due to curvature incorporated in the height of the single-layer model (H _eff) by considering spatial and temporal variation.

For a single-layer model, the height of the ionospheric layer is therefore an important parameter which influences the estimated coefficients. Also, the modelled VTEC is limited by the obliquity factor used to map the STEC to VTEC. Since the effective height of the ionosphere is known to vary both temporally and spatially, it is important to model the ionosphere using a three-dimensional spatial model. We conclude that the future work should focus on the construction of a three-dimensional (or multi-layer) model for the ionosphere.

We have also found that the addition of GLONASS data to GPS data, and the use of a multi-station solution rather than a single-station solution, gives better results than our original work. The gradients from the multi-station approach were estimated at a higher time resolution (2 min) in comparison to single-station approach. Also, due to the large number of observations used to estimate gradients in a multi-station approach, the gradients seem to have a smoother temporal variation. For all the selected days of MWA observations, the correlation between MWA and GNSS estimated gradients was found to be identical within errors or higher with multi-station approach as compared to single-station approach.

The ionospheric modelling performed using GPS and GLONASS observations was also able to capture the spatial variations of the gradients, refer Figure 9. This encourages us towards deriving the higher order effects of the gradients estimated using GNSS. Future work will focus on deriving higher order effects in the gradients for various sources within MWA FoV and GNSS observations.

The NS gradients, estimated using the multi-station approach, agreed with the MWA observed gradients. The EW gradients had a better correlation than single-station approach; however, did not seem to correlate as well as the NS gradients. The current distribution of GNSS receiving stations, while demonstrated to be successful in characterising large-scale ionospheric features and validating our technical approach, is inadequate for advanced modelling.

The MWA with its wide FoV imaging capability, sees a position offset due to the ionosphere for each of the sources in its FoV. This capability of the MWA has been exploited to detect small-scale structures in the ionosphere (Loi et al. 2015a, 2015b, 2015c). The spatial scales of the structures are around 10–100 km. These are precisely the spatial scales that will need to be characterised if future longer baseline instruments are to be calibrated.

GPS satellites are capable of providing ionospheric information; however, the density of the pierce points is far too low to probe these small scales for the datasets currently available for the MRO. For a GPS receiver near the MRO, only 5–8 satellites are visible above horizon at any given time. The number of measurements can be increased to 10–15 satellites by including data from GLONASS satellites. Future Global Navigation Satellite Systems (GNSS), namely, BeiDou (China) and Galileo (European Union) are expected to be operational around the year 2020 (UN 2010). Regional satellite systems such as the Quasi-Zenith Satellite System (QZSS, Japan) (UN 2010), currently under development, will have four satellites, all of which pass over the MRO. In this scenario, 20–30 satellites would be above the horizon at the MRO at any given time, with an average separation of 200 km on the ionosphere.

We have established a methodology to obtain ionospheric information with the current GNSS infrastructure around the MRO. This methodology could be exploited to derive ionospheric corrections for the future LOFAR, like the SKA-low, where the direction-dependent effects become more dominant and deriving ionospheric information from the astronomical observations may not be feasible. The GNSS-derived ionospheric information can also be used for climatology. This could be useful in designing future instruments, devising calibration strategies, and for selecting data post-observation (to avoid wasting effort on data which is too badly corrupted by the ionosphere). The current GNSS infrastructure which limits the spatial resolution of the ionospheric corrections can be improved by deploying additional GNSS receivers.

To measure the ionosphere on scales of < 10 km, a dense cluster of GNSS receivers (5–10 receivers), with baselines as small as ~ 5 km would need to be installed. This would allow the ionospheric gradients, rather than just the STEC, to be determined towards each satellite. However, it is not guaranteed that a GNSS satellite would be in the MWA FoV at all times. Hence, another small cluster of GNSS stations would need to be deployed strategically. Deploying the further cluster at a distance of ~ 100 km would fill in the gaps between existing satellites. We call this approach the ‘cluster of clusters’ approach. For further densification, a cluster of GNSS stations at the median of existing clusters ( ~ 50 km) could be deployed.

The sparse population around the MRO and the lack of remote power and communication infrastructure constrains the possible locations for such a cluster. Locations indicated in Figure 12 are plausible cluster locations given that they are existing communities and homesteads likely to have the necessary infrastructure.

Figure 12. Current (green) and proposed (red) GNSS station locations in vicinity of MRO. The MWA location is marked in blue.

Future work will evaluate the expansion of GNSS stations around MRO, in view of generating a three-dimensional ionospheric model, to meet the ionosphere calibration requirements for future MRO instruments.