Enhanced wide-area multi-GNSS RTK and rapid static positioning in the presence of ionospheric disturbances

Ionospheric disturbances are the phenomena which adversely affect the performance of precise positioning. This holds true even for multi-constellation relative positioning supported with network-derived ionospheric corrections. In such scenario the unfavorable effect is caused by a poor accuracy of corrections, which, in turn, is driven by the deterioration of the spatial interpolation process. The positioning becomes even more challenging in a wide-area scenario with baselines over 100 km. In this paper, we assess the methodology which aims at reliable and accurate wide-area RTK and rapid static positioning in the presence of severe ionospheric conditions. The approach takes advantage of multi-constellation network ionospheric corrections and an algorithm which allows the elimination of the temporal variations of the ionospheric delay. The experimental evaluation was performed on the basis of multi-station RTK and static positioning using GPS, BDS and Galileo data collected at high latitudes during the ionospheric storm on August 25–26, 2018. The results confirmed the deterioration of the accuracy of the network ionospheric corrections and consequently a decline in the positioning performance with routine models such as ionosphere-float and ionosphere-weighted. On the other hand, the results obtained with the application of the developed methodology demonstrated a very distinctive improvement in the ambiguity resolution domain and thus proved the advantage over benchmark models. In this case, the developed methodology allowed up to 20% enhancement of the ambiguity success rate with respect to benchmark strategies.


Introduction
Over the last decades a number of GNSS positioning techniques have been developed. Depending on the classification, the methods may be referred to as absolute and relative, code-and phase observable-based or real time and postprocessed (Katsigianni et al. 2019). Nonetheless, all of these methods suffer from the impact of severe ionospheric activity (Kim and Tinin 2007;Tiwari et al. 2009;Bergeot et al. 2010). Such conditions lead to the significant growth of plasma content, which in the extreme cases may amplify the ionospheric delay of GNSS signals up to 100 m. More importantly, the storm-time ionosphere and consequently time series of satellite data are characterized by high temporal variability as well (Jin et al. 2017;Su et al. 2019).
Since the ionospheric irregularities and their impact on GNSS positioning are the subject to intensive interdisciplinary studies, a number of them have been recognized so far. Rishbeth (2000) characterized ionospheric equatorial anomaly, which is responsible for creating an electron concentration trough at geomagnetic equator as well as north and south crests at ~ 15° of geomagnetic latitude. Such phenomena are said to be responsible for adverse GNSS signal scintillations, which in turn cause the occurrence of cycle slips in phase data and distort the positioning performance. Other frequently investigated ionospheric phenomena are Traveling Ionospheric Disturbances (TID), which are defined as wave-like fluctuations in Total Electron Content (TEC) that propagate through the ionosphere (Hines 1960;Tsugawa et al. 2007;Deng et al. 2013). It is common to classify the TIDs with the size criterion . While the large-scale TIDs are believed to be triggered by geomagnetic storms, the medium-scale structures (MSTID) are induced by the atmospheric gravity waves in the neutral atmosphere as well as the auroral substorms (Hunsucker 1982;Hocke and Schlegel 1996). The modeling of MSTID with GNSS signals as well as their impact on precise positioning performance has been reported in several papers. The studies given by Saito et al. (1998), Otsuka et al. (2013), Chen et al. (2016), Savastano et al. (2017) who proved a high applicability of GNSS observables for MSTID characterization. On the other hand, although these structures are relatively weak (up to a few TEC units), they can significantly deteriorate the estimated coordinates, and thus a reliable GNSS positioning in such conditions is still considered as a challenging task (Sieradzki andPaziewski 2015, Hernández-Pajares et al. 2017). A noticeable GNSS signal degradation is especially produced by high latitude ionospheric irregularities over polar and auroral regions (Spogli et al. 2009;Sieradzki 2015;Cherniak and Zakharenkova 2016). The occurrence of the disturbances for this area is considered to be driven by the coupling between ionosphere, solar wind and the magnetosphere. This system, which is highly dependent on Sun's activity, is responsible for particle precipitation and a convection pattern at high latitudes. The former effect implicates the generation of small-and mediumscale irregularities, which are particularly intense at the auroral oval. The latter, in turn, makes it possible to feed polar ionosphere with plasma from middle latitude reservoir and a formation of polar patches (Carlson 2012;Sieradzki and Paziewski 2019). The occurrence of the above structures causes the extreme complexity of ionosphere and frequent scintillations of GNSS signals (Jiao et al. 2013;Prikryl et al. 2014). Due to this reason the positioning at high latitudes has to be considered as the worst-case and most challenging scenario. On the other hand, there is no doubt that testing of algorithms in such conditions should provide a reliable evaluation of their efficiency.
To overcome the negative impact of the ionospheric delay on GNSS positioning several approaches have been developed over the years. The first group of such algorithms implements and employs empirical ionospheric models, such as International Reference Ionosphere (IRI), Klobuchar, NeQuick and Neustrelitz, however, due to the high level of spatiotemporal generalization, they do not meet the requirements of rapid precise positioning even for quiet conditions. Another approach, mainly applied in post-processing applications, utilizes GNSS-based global ionospheric maps (Zhang et al. 2016). These maps, however, feature insufficient spatio-temporal resolution to capture the disturbances during high ionospheric activity (Park et al. 2016). Analyzing the other possibilities it should be noted that the benefit from the application of ionosphere-free linear combination does not offset the drawback of providing a float-only solution. Hence, a method which is currently the most predestined to provide reliable and precise solution is relative positioning supported with real-time network-derived ionospheric corrections (Rizos 2002). This algorithm allows the achievement of cm-level accuracy position even for the long-range baselines ). On the other hand, it is extremely difficult to provide accurate ionospheric corrections in the presence of the severe ionospheric disturbances Paziewski 2016;Jacobsen and Andalsvik 2016). To handle the impact of the ionospheric disturbances in RTK positioning, the researchers proposed to support RTK users with indices providing statistical information on expected residual ionospheric errors (Wanninger 2004). For such a task Lejeune and Warnant (2008) proposed to employ a network-based monitoring service. Other studies suggested the mitigation of the ionospheric scintillations by advanced stochastic modeling of observables (Park et al. 2017). Nevertheless, the reliable RTK positioning is still a demanding task under the disturbed ionosphere.
In this paper, we demonstrate the advantages and assess the methodology which aims at reliable wide-area multi-station RTK and rapid static positioning in the presence of severe ionospheric conditions. The presented approach utilizes multi-constellation network ionospheric corrections and an algorithm suited for the elimination of the temporal ionospheric delay variations-rate of TEC correction (RTC).
The study is organized as follows. In the following section we briefly describe the developed methodology for wide-area relative positioning under ionospheric disturbances. The subsequent section characterizes the ionospheric conditions during the period employed for the experimental evaluation of the methodology. Next, we present a performance assessment of the enhanced widearea positioning for different ionospheric conditions. Finally, conclusions are drawn in the last section.

Methodology for enhanced multi-constellation wide-area relative positioning under high ionospheric activity
The enhanced methodology of positioning consists of several consecutive steps performed at the reference network and rover sides: In the first step the reference network is processed in a relative mode with station coordinates held fixed (step 1.a). Taking advantage of the resolved double differenced (DD) ambiguities (N), corresponding precise DD ionospheric delays (I) may be obtained for each of the processed baselines (step 1.b), employing a geometry-free linear combination of dual-frequency phase observables (Schaer 1999): where the superscripts i and j, and the subscripts k and l identify the satellites and stations, respectively, φ is the phase observables in cycles, λ denotes the wavelength, f corresponds to the signal frequency.
Such obtained precise DD ionospheric delays are decomposed to undifferenced values and interpolated using ordinary kriging method to an approximate position of the rover receiver (Wielgosz et al. 2003). This procedure takes advantage of the least squares estimation with parameter constraining to handle rank deficiency of the decomposition model. Finally, the absolute ionospheric corrections are converted to DD values and applied to DD observables for the rover site as described in detail in Paziewski and Wielgosz (2014). (1) Due to the existence of non-overlapping signal frequencies in the used constellations (GPS, Galileo and BDS), we take advantage of a loosely combined model for the integration of multi-constellation signals. As former studies show, both tightly and loosely combined models are characterized to exhibit similar performance in the case of unobstructed visibility of satellites . Due to a noticeable impact of the ionospheric delay, the methodology assumes not only support from networkderived ionospheric corrections but also parametrization of the residual slant ionospheric delays. In this context, the observational model is often termed as ionosphereweighted (Bock et al. 1986;Teunissen 1998;Kashani et al. 2007;Paziewski 2016). Hence the system of observation equation with intra-system double-differenced dual-frequency observables reads: in which T denotes zenith tropospheric delay, α is the coefficient of the tropospheric mapping function, P is the measured pseudorange, geometric distance is denoted as ρ and finally μ stands for the constant coefficient employed for conversion of the ionospheric delay on the first applied frequency into that on the second one, according to the equation: Corresponding vector of estimates reads: in which corrections to a priori coordinates were denoted as {δX, δY , δZ} , the zenith tropospheric delay estimates are given as δT , δI denotes a set of epoch-wise DD ionospheric delays and N stands for a set of phase ambiguities.
The ionosphere-weighted model takes advantage of constraining of the DD slant ionospheric delays (DD-SID) to the a priori values, which is done here by the introduction of the pseudo-observations with assigned variance factors. A similar approach is applied for zenith tropospheric delays as well as for station coordinates. As a result the combined vector of pseudo-observations aggregates on the left side of the estimated parameters and a priori values of the constrained parameters on the opposite. (2) (4) x = [δX, δY , δZ, δT , δI, N ] T , Page 4 of 16 Paziewski and Sieradzki Earth, Planets and Space (2020) 72:110 This well-recognized observational model is fed with corrections which eliminate temporal variability of the ionospheric delays, as validated by the authors in Sieradzki and Paziewski (2016). Taking advantage of rate of TEC equation, which is basically the time difference of the geometry-free linear combination ( L 4 ), we may obtain temporal change of the ionospheric delay between epoch i and j (Eq. 6): where geometry-free linear combination at selected epoch i, is formed of dual-frequency phase observables according to Eq. (7): A corresponding equation may be formed at consecutive epoch j, as well. Hence the t ij epoch difference of L 4 combination yields: Providing that there were no cycle slips between epochs t i and t j , the phase ambiguities are constant, which allows for the derivation of a temporal change of the ionospheric delay directly from equation: Such a derived increment of the ionospheric delay corresponds to a geometry-free linear combination but it may be easily converted to any signal frequency and applied as a correction to observables. As a result, the observables are free from the temporal variability of the ionospheric delay since updated observables of following epoch t j are subject to the ionospheric delay equal to epoch t i (Sieradzki and Paziewski 2016).
In a regular ionosphere-weighted model the ionospheric delays are considered as epoch-wise parameters. The main benefit of RTC is that the parametrized ionospheric delays are considered in the model (Eq. 2) as time-constant parameters, hence are accumulated estimates. RTC corrections cannot be considered as temporal interpolation in strict sense but fulfills an equivalent function, i.e., removes ionospheric fluctuation in time domain. As a reference epoch we adopted the epoch with the lowest STEC derived directly from L4 time series. The reason of that is positive enhancement of plasma for auroral irregularities. Consequently, the delay at the reference epoch is close to the background state of ionosphere, which is expected to be similar for all stations in network and will be eliminated for double-differenced observations. The algorithm is employed after the application of the network ionospheric corrections, therefore the residual ionospheric delays are not only constant but also reduced with respect to the original values.
Stochastic model incorporates weights corresponding to the linearized phase and code observables as well as the pseudo-observables of constrained parameters. Since the functional model is based on a relative approach, the former weight matrix must take into account the inter-satellite and inter-station correlations resulting from double differencing process as well as variances assigned to elevations and a priori zenith variances. The weight matrix of DD-SID pseudoobservables also takes into account the mathematical correlation and a priori variances of ionospheric delay assigned to the expected residual ionospheric error. Hence this approach is an advance to basic stochastic modeling which provides equal variances of all DD-SID pseudo-observables (Paziewski and Wielgosz 2014).

Characteristics of the ionospheric conditions during the validation experiment
The experimental validation of the methodology was performed with GNSS data collected on August 25 and 26, 2018 by permanent stations located at the northern high latitudes and marked with red dots in Fig. 1. We preceded this validation with the characterization of ionospheric conditions over the experimental area. The carefully chosen test period corresponds to the one of the strongest geomagnetic events in the recent years. The cause and the unexpected strength of this storm are still not fully explained, but it is thought that it was originated by a coronal mass ejection (CME) erupted from the Sun on August 20, 2018. The arrival of this slow-moving CME on August 25th caused the southward reorientation of interplanetary magnetic field (IMF), which in turn initiated the minor geomagnetic storm (G1) in the last hours Page 5 of 16 Paziewski and Sieradzki Earth, Planets and Space (2020) 72:110 of that day. The continuation of steady negative values of IMF Bz on August 26, permanently exceeding -10 nT, and a magnetic reconnection related to such conditions led to the intensification of the storm to G3 level. For the characterization of ionospheric conditions during the storm, we used relative STEC (slant total electron content). This GNSS-based parameter has been recently proved to be effective for the detection of the differentscale ionospheric disturbances (Sieradzki and Paziewski 2019). It is computed using arcs of geometry-free combination data (L4) in two steps. The initial step aims at cycle-slip repairing and is realized with the leveling of phase to data to code one according to the equation: where L 4(k) , L 4(k) and P 4(k) denote the leveled phase, phase, and code data at epoch k.
Equation (11) is applied to clean time series longer than 15 min. In the opposite case, a poor-quality phase data are substituted by the five-epoch moving average of −P 4(k) . The second part of the algorithm is the computation of relative STEC values defined as the difference between epoch-wise ionospheric information L 4(k) and a background STEC level. The latter factor is approximated with three-step iterative fitting of a fourth-order polynomial applied to clean arc of L 4(k) values. After each solution, we removed the epochs with differences higher than 1.5 TECu from further processing to stabilize the fitting at the background level of STEC variations.
The starting point of the analysis is an overall characterization of ionospheric conditions prevailing at northern high latitudes. For this purpose, we used the relative STEC values derived from approximately 180 stations located in the northern hemisphere. The summary of (11) L 4(k) = L 4(k) − median(L 4(k) + P 4(k) ) arc , these results is given in the bottom panel of Fig. 3. It shows the temporal variations of the relative STEC values as a function of the altitude-adjusted corrected geomagnetic (AACGM) latitude (Shepherd 2014) for day-and nightside ionosphere (6-18 and 18-6 of Magnetic Local Time, respectively). The positions of the ionospheric pierce points were computed using single layer model at 350 km attitude; however, they were prepared in a different way for auroral (< 77° of AACGM latitude) and polar regions. In the former case the original ACCGM latitude was used, whereas in the polar region we adopted the points resulting from the preliminarily projection on noon-midnight axis. Such distinction allows the depiction of short-term variations of the auroral oval as well as the convection of plasma across the cap. In order to clarify the occurrence of irregularities during the test period, the two top panels in Fig. 2 depict the variations of the selected geomagnetic indices, i.e., IMF Bz and SYM/H adopted from OMNI database service (http://omniw eb.gsfc.nasa.gov) as well as AE provided by World Data Center for Geomagnetism (http://wdc.kugi.kyoto -u.ac. jp). The latter parameter is a real-time product based on 10 observatories.
The results clearly confirm the dependence of ionospheric irregularity occurrence on geomagnetic conditions. The initial hours of the test period (0:00-16:00 UTC on August 25, 2018) were characterized by northwardly oriented IMF and consequently all of the high latitude ionosphere was quiet. Considering the experimental network area, marked with black line in Fig. 2, the relative STEC were continuously below 1 TECu. Thus, the positioning for this interval is free from the disturbing impact of ionosphere and consequently can be treated as a benchmark for different scenarios. The following hours (from 16:00 UTC on August 25 until 9:00 UTC on August 26), involving the main phase of the storm, were the extensively long period of negative IMF Bz. As expected it was characterized by the occurrence of different ionospheric structures. In terms of positioning performance the most relevant are small and medium irregularities in the auroral oval area. According to Fig. 2 their occurrence during this period was practically continuous with higher intensity observed on the nightside. The detected disturbances are coincident with rapid AE variations, which suggest the particle precipitation as their primary source. The STEC enhancement for these structures often exceeds 10 TECu, making GNSS positioning an extremely challenging task. Apart from the development of auroral oval, the results of relative STEC also depict other structures, the first of which were polar patches detected for cap area within the interval 20:00-23:00 UTC on August 25. Their occurrence seems to be related to the positive effect of the storm and increased amount Page 6 of 16 Paziewski and Sieradzki Earth, Planets and Space (2020) 72:110 of plasma on the dayside. The increase of high-latitude TEC for this period was also reported by Blagoveshchensky and Sergeeva (2020). The another consequence of geomagnetic conditions was the formation of a tongue of ionization extending from post-noon to pre-midnight sectors of MLT known as storm enhanced density (SED). Its signatures on both sides of ionosphere are clearly visible at about 55° of AACGM latitude (1:00-3:00 UTC on August 26). Since these phenomena are extensively large, their impact on positioning is expected to be different than auroral oval. Unfortunately, they occurred far from the test network and could not be analyzed. The ionospheric activity during the recovery phase of the storm confirms its strong dependence on geomagnetic conditions. As one can observe the more intense variability at auroral latitudes corresponds to two periods: 15:00-17:00 UTC and 19:00-21:00 UTC, which were characterized by southwardly oriented IMF and increased values of AE index. The following step is the more detailed analysis of ionospheric conditions at the experimental network. Thus, Fig. 3 provides the relative STEC from a rover site HETT and an outermost reference point-KEV2. For the clarification purpose it presents the results from the selected satellites (GPS and GLONASS) covering the entire test period. In this case we applied the elevation mask set to 10°, which is commonly adopted for GNSS positioning. The time series given in Fig. 3 are consistent with the global characterization of ionospheric activity presented above, both for the quiet prestorm period as well as for the following phases of the event. The one exception, observed between 6:00 and 10:00 UTC on August 26, is related to the location of the test network below the equatorward boundary of auroral oval. As demonstrated the disturbances over the experimental network typically reached 4-8 TECu.
Only for the period of 2:00-4:00 UTC on August 26 they exceeded 10 TECu. These strongest structures were observed only locally at dawn sector of Magnetic Local Time. Analyzing the particular periods there is no difference between the scale of ionospheric variability at specific satellites. This may be explained by the latitudinal range of the oval and indicates that for the specific intervals all of the GNSS data are affected by Page 7 of 16 Paziewski and Sieradzki Earth, Planets and Space (2020) 72:110 the occurrence of ionospheric disturbances. The same agreement applies to the time series derived for both stations. On the other hand, the results for both sites reveal strong differences in the relative STEC at particular epochs. In such conditions the degradation of GNSS positioning performed with common algorithms is rather expected. In order to depict the discrepancies in temporal variations of STEC, Fig. 4 presents data of two GPS satellites registered at all of the stations used in the positioning experiment. In particular, the time series of GPS PRN 7 correspond to the typical disturbed conditions at high latitudes, whereas these of GPS PRN 19 demonstrate the strongest structures detected during this event. The initial batch of results given in the top panel (20:00-20:45 UTC) depict the common pattern of small-scale irregularities occurring in the auroral zone for southwardly oriented IMF. Due to their size, the ionospheric pattern varies significantly between stations and results in differences at the level of 3-4 TECu. The following larger structure is basically observed for all the sites, but its signature varies between them. This is related to the lengths of the baselines in the network and the co-occurrence of the smaller ionospheric structures as well. The combined impact of these factors amplifies the observed discrepancies, which reach as high as 6 TECu. The case presented in the bottom panel has to be considered as the most challenging scenario for positioning, wherein the enhancement in the amount of plasma lead to the differences exceeding 10 TECu. Their maximal values are mainly a consequence of a high gradient of STEC at the edges of structures. In this case one can also notice a clear dependence of results on the distance between stations. The best examples of this effect are the values for the outermost reference station KEV2, which do not reveal the structures at 3:15 and 3:45 UTC, respectively. Finally, the selected time series demonstrate the high elevated data, i.e., about 65° and 55° at the middle epochs for GPS PRN 7 and PRN 19, respectively. It means that the discrepancies depicted here will propagate to all double-differenced observables assuming these satellites as pivot ones.

Experiment design and assessment of the network ionospheric corrections
In the experiment, we have used GNSS observations collected by high latitude stations of EUREF Permanent Network (EPN). Such dataset ensures harsh ionospheric conditions, which are advantageous for a comprehensive examination of the proposed method's performance. The stations offer the observations with 30 s sampling interval, which is not favorable for practical geodetic and surveying applications. However, as demonstrated in the previous section, the signatures of ionospheric structures have periods of at least several minutes, and thus we should not expect any significant benefit in a positioning performance when using 1-s data instead of 30 s. There is no doubt that the temporal variations of ionosphere are lower in such case, but the consecutive epochs with shorter interval would be affected by a similar ionospheric gradient between stations. As a result, the 1-s data bring little new information and enhance the solution only to a small extent. Figure 5 presents the experimental network used for the evaluation of the positioning. We have employed HETT station as a simulated rover receiver and four neighboring stations (KILP, KEV2, SOD3, KIR8) as the reference network. As a result the length of the baselines in the rover solution fits the range of 123-201 km and hence the positioning may be considered as wide-area. Moreover, the stations fulfill the requirements of multiconstellation signal collection thus GPS, Galileo and BDS signals were employed. As shown in Figs. 2 and 3, the date of the experiment ensures the periods with high and low ionospheric activity, making it feasible to validate the methodology in different ionospheric conditions. Figure 6 depicts the number of satellites tracked over the analyzed period as well as corresponding PDOP value. As can be seen in the figure, the total number of satellites ranges between 11 and 21 with PDOP always smaller than 1.5. This implies that the satellite geometry conditions may be considered good provided that multiconstellation signals are used.

Assessment of the network ionospheric corrections
The methodology of wide-area relative positioning takes advantage of network ionospheric corrections, which are spatially interpolated to the position of the rover receiver. In this case, there was no need for a temporal interpolation of the corrections since their sampling interval equaled to the sampling interval of the rover solution. It  is theoretically the best case scenario, since only spatial interpolation is applied without any processing in time domain. Unfortunately, taking into account the dynamics of the high latitude ionosphere and the size of the auroral irregularities, the application of the network ionospheric corrections, even with adequate time fineness, still does not guarantee a high efficiency of error mitigation process. We attribute it to a poor quality of the network ionospheric corrections, which, in turn, is driven by the deterioration of the spatial interpolation . This phenomena may be explained by a strong impact from small, locally produced, ionospheric irregularities. It is clear that the degradation of the corrections' accuracy will be even higher if they are characterized with time fineness inadequate to the ionospheric conditions,. This, however, strongly depends on the dynamics of the processes in the ionosphere (Kashani et al. 2007).
To verify these hypotheses, we have assessed the accuracy of network-derived ionospheric corrections with our experimental dataset during high and low ionospheric activity. Figure 7 presents the ionospheric correction residuals and true delays given for the baselines denoted in Fig. 5. Top panel corresponds to the true ionospheric delays which were obtained from a geometry-free linear combination (Eq. 1) after an application of integer ambiguities from the solution with station coordinates held fixed. Bottom panel corresponds to the ionospheric correction residuals defined as the differences between the interpolated network-derived corrections and the true ionospheric delays. Figure 7 clearly depicts the periods with high and low ionospheric activity (from 18:00 UTC on August 25, until 24:00 UTC on August 26 and 00:00-18:00 UTC on August 25, respectively). During the former period the DD delays reached up to 5 m, whereas the low activity period is characterized by the delays fitting the range of ± 0.5 m. What is more important, Fig. 7 illustrates a clear deterioration of the ionospheric correction accuracy during the ionosphere active period. Comparing both panels we see that residuals are only slightly lower than true ionospheric delays. This confirms that spatial interpolation of auroral oval irregularities was of poor performance. Table 1 summarizes the accuracy of the ionospheric corrections. As we can read from the table, during low ionosphere activity it was feasible to obtain high accuracy corrections. There were over 81% of the corrections with residuals lower than 5 cm and almost 96% lower than 10 cm. The standard deviation (STD) of the correction residuals equaled to 4.6 cm. Such accuracy was, however, not achievable during the disturbed ionosphere. As we can learn from Table 1, only 39% of the corrections were characterized with residuals fitting the range of ± 5 cm and 60% fitting the range of ± 10 cm. As shown in (Paziewski 2016) such accuracy of the corrections is a prerequisite for reliable and fast RTK positioning. Similarly, also STD of ionospheric corrections was at a significantly higher level with respect to the inactive ionosphere period and reached 31 cm. Such a number indicates a low accuracy of the corrections during the active ionosphere, which seems to prevent high-performance RTK positioning.

Performance assessment of enhanced wide-area positioning under disturbed and quiet ionosphere
This section evaluates the performance of the proposed enhanced positioning model. As the benchmark solutions we used two well-known models (#1 and #2), therefore eventually three strategies were tested:  1. Ionosphere-float model (IF) which parametrizes DD-SID, This model is not supported by any ionospheric corrections. 2. Ionosphere-weighted model (IW) which parametrizes the residual DD-SID after an application of the interpolated network ionospheric corrections. We consider this approach as superior to #1. 3. Enhanced wide-area positioning model, which takes advantage of modified IW model with RTC (IW-RTC).
The computations were performed for both rapid static and RTK modes using imitated short time sessions of 15-min-long and 30-s interval. This means that each session lasted only 30 epochs and we processed 192 sessions during the analyzed 2 days of the experiment. Table 2 summarizes other details of the data processing.
The above methodology is one of the common standards of relative positioning. The only uncertainty is related to the lower accuracy of broadcast ephemeris of BDS system, what is particularly true for geostationary (GEO) satellites (Montenbruck et al. 2015;Ouyang et al. 2019). However, considering the localization of the experiment in the northern Europe and applied elevation mask, the acquisition of GEO BDS signals was strongly limited. Since we have collected the observations only from a single low-elevated GEO BDS satellite, its impact on the positioning can be considered as negligible. Taking into account the error of broadcast ephemeris of 2 m and the range of employed baselines, we expect after Shi et al. (2017) that this may produce a deterioration in the coordinate domain of about 1 cm. This value cannot be neglected in the most precise applications. However such error should not importantly impact on the ambiguity resolution, which is a key performance domain in our analysis. Moreover, such impact of the ephemeris error is much lower than the temporal dynamics of the ionosphere.

Ambiguity resolution domain
As the indicators of the ambiguity resolution performance we used mean time to fix (TTF) and ambiguity success rate (ASR). ASR is here defined as a ratio of sessions with correctly fixed ambiguities with respect to the number of all sessions. TTF is a number of epochs required to obtain a solution with correct integer ambiguities.

Low ionospheric activity period
We begin the evaluation with a brief demonstration of models' performance during low ionospheric activity. Under such conditions all models, namely: IF, IW and IW-RTC provided fast converging solutions in both RTK and static modes. The detailed results proved that the models required up to about 2 epochs to obtain a correct integer ambiguity resolution. Specifically, in the RTK mode the mean TTFs were of 1.9 and 1.5 epochs for the benchmark IF and IW models, respectively, while for IW-RTC this statistics dropped to 1.4 epochs. As expected, the rapid static mode provided comparable high-performance results. In this case the ionosphere-float model required the longest time to reach a correct integer ambiguity resolution and thus a cm-level positioning, which was expressed by the highest mean TTF of 2.2 epochs. Although the TTF of IW and IW-RTC were close, it is still discernable that for the proposed model this statistics was slightly lower than for the ionosphere-weighted one being 1.4 and 1.6 epochs, respectively.
The results of ASR as a function of the session length confirmed a high performance of all analyzed models during the quiet ionospheric period. In specific, the IW and IW-RTC models applied to the RTK mode allow obtaining over 90% of sessions with correctly resolved  Paziewski and Sieradzki Earth, Planets and Space (2020) 72:110 ambiguities with just a single epoch of observations. In the case of the model deprived of any ionospheric corrections (IF), this parameter reached noticeably lower value of 82%. Correspondingly, after 30 epochs, IW and IW-RTC were characterized with ASR of 100%. Such uppermost ratio was not attainable for the IF RTK, nevertheless still a high value of 98.6% was reached. As expected, also in the rapid static mode all analyzed strategies allowed the achievement of a high ratio of sessions with correctly fixed ambiguities. After 10 epochs the mean ASR reached the level of 98.6% for IW-RTC, whereas for the benchmark strategies (IF and IW) these statistics were slightly lower of 95.8% and 97.2%, respectively. Finally, the mean ASR corresponding to 30 epoch long solution fitted the range of 98.6-100%, again with the highest value for the IW-RTC model. These results demonstrate that during quiet ionospheric conditions the advantage of the enhanced model over the ionosphere-weighted model is noticeable, but not significant. We attribute a high performance of the latter approach to the quality of the network corrections which allowed effective mitigation of the ionospheric delay. We may conclude with a high degree of certainty that such ionospheric conditions were not challenging and therefore it was feasible to handle them with all ionosphere modeling strategies. Hence a more comprehensive validation of the methodology is performed in the following subsection taking advantage of the harsh ionospheric conditions.

High ionospheric activity period
In Fig. 8 we depict the time to fix given for both RTK and rapid static modes during the presence of ionospheric disturbances. We can clearly notice that a much longer time was required to obtain a correct ambiguity resolution with respect to the period of quiet ionosphere. In many cases, especially during the most disturbed periods (from 18:00 UTC on August 25 until 6:00 UTC on August 26 and 15:00-24:00 UTC on August 26), it was not possible to converge to a correct fix within a 30 epoch time limit. These sessions are depicted with grey bars in the figure.
The values of mean TTF for the RTK mode reached up to 7.6 and 5.7 epochs for the IF and IW benchmark models, respectively. The bottom panels of Fig. 8 show the effect of a better ionospheric delay handling of the IW-RTC model. In the case of RTK mode the corresponding mean TTF equaled to 5.4 epochs. This value is slightly superior to that of IW, however the figure demonstrates that much more session converged within the 30 epoch long period with respect to benchmark strategies.
In analogy to RTK mode, high ionospheric activity constitutes a challenging scenario for a rapid static positioning. The right column panels of Fig. 8 illustrate that ionospheric disturbances are the main factors contributing to the deterioration of the ambiguity resolution. What follows from the figure is that even a lot more sessions did not converge to a fix within a 30 epoch time limit in the case of the benchmark strategies with respect to IW-RTC. We have discovered that the mean TTF of IW-RTC was superior to that of the benchmark models, as we obtained a clear drop from the level of 6.9 epochs down to 5.3 epochs.
A clear advantage of the enhanced model over benchmark strategies is justified by Fig. 9 which depicts ASR as a function of session's length during high ionospheric activity. In the RTK mode after 10 epochs ASRs are 53.4%, 55.2%, 72.4% and after 30 epochs are 75.0%, 71.6%, 88.8% for the IF, IW and IW-RTC models, respectively. Even more noticeable discrepancies between the strategies are reflected in these statistics given for the static solution. The bottom panel of Fig. 9 shows that after 10 epochs it was feasible to reach the mean IW-RTC ASR of 75.2% and only of 52.1% and 54.7% with the benchmark models IF and IW, respectively. After 30 epochs this parameter reached the level of 74.4% for the benchmark strategies and 94.0% for the IW-RTC model. The close ASR results obtained for the strategies IF and IW indicate that there was no benefit from the network corrections, which supported the latter model. We attribute such outcomes to the poor quality of such corrections. This may be related to the occurrence of small ionospheric irregularities as well as a high TEC gradient which may occur at the equatorward boundary of the auroral oval. The former factor cannot be mitigated with any spatial interpolation, whereas the latter may be only to a some extent. This however depends on the ionospheric conditions at particular reference stations. These results are basically consistent with the statistics of the ionospheric correction residuals presented in Sect. 4. On the contrary, the application of the IW-RTC model yields much better performance. Particularly, this strategy is characterized with up to 20% higher values of the ASR over the benchmark ones. A clear improvement observed for the IW-RTC model confirms that the impact of small ionospheric irregularities on GNSS positioning can be effectively reduced with the Rate of TEC Corrections.
The results of ASR as a function of session's length gives also the presumption on the models' performance in the case of application of 1 s sampling interval. With regard to the ionosphere-weighted model we should not expect a significant shortening of the time required for an ambiguity fixing since the temporal change of Paziewski and Sieradzki Earth, Planets and Space (2020) 72:110 the ionospheric delay for 1 s data is much lower than its spatial disturbances. In the case of enhanced model, the efficiency will be strongly dependent on the spatial size of the structures (their signatures in GNSS data). For smaller disturbances we can expect that the background level of the ionosphere will be faster identified with RTC and thus, the improvement in a convergence time is possible. In the opposite case the efficiency of RTC should be at the similar level as in our experiment.

Coordinate domain
In the previous section we have demonstrated a noticeable benefit from the IW-RTC model in terms of the ambiguity resolution performance during high ionospheric activity. This section investigates a potential impact on the coordinate domain. To this end, we analyzed the statistics such as the standard deviation of the topocentric coordinates, being a measure of position repeatability, as well as the mean coordinate residuals with respect to the ground truth position. The summary Fig. 8 Time to first fix for RTK (left column) and rapid static (right column) modes obtained in the consecutive sessions during high ionospheric activity given for IF (top panels), IW (middle panels) and IW-RTC models (bottom panels). Black bars depict sessions without a solution, grey ones correspond to the sessions which were not correctly fixed within the 30 epoch time limit Fig. 9 Ambiguity success rates as a function of sessions' length obtained in a cumulative RTK and rapid static mode solutions during high ionospheric activity Page 13 of 16 Paziewski and Sieradzki Earth, Planets and Space (2020) 72:110 of these results, given separately for float and fixed solutions, is provided in Table 3. It is expected that after the correct ambiguity resolution the differences in coordinate precision between the strategies may not be significant. This is due to the fact that the accuracy of the fixed solution is driven mainly by the precision of phase observations. This point is confirmed by the detailed results. In general, STDs of the fixed RTK and rapid static solutions obtained with the IW-RTC model are in the range of 12-16 mm, 9-10 mm for North and East components, respectively. The benchmark strategies IF and IW provided slightly lower coordinate repeatability which is reflected by STD fitting the range of 10-19 mm, 10-13 mm, respectively. The height component was obtained with a noticeably lower accuracy than that of the horizontal coordinates, since the STDs were in the range of 54-64 mm for the IW-RTC model and 63-100 mm for the benchmark strategies. Such level of accuracy was of course expected considering the wide-area scenario and relatively short sessions. Some improvements in this field may be expected after an application of more advanced tropospheric delay mitigation methods (Wielgosz et al. 2012). This issue was however not comprehensively investigated in this study, since we focused on the impact of ionosphere and its influence on the ambiguity resolution domain. Even so, detailed results show that the IW-RTC model outperforms the benchmark strategies, which is reflected, e.g., by lower STD of height component. For example, during high ionospheric activity this statistics dropped from 92 mm for IF RTK down to 54 mm for the IW-RTC.
All of the strategies allowed obtaining comparable mean coordinate residuals related to the ground true plane coordinates. These statistics fitted the range of 0-3 mm and 0-5 mm for North and East components, respectively. On the other hand, the results showed a clear benefit from the application of the IW-RTC model to the reduction of the height bias during the high ionospheric activity period. The application of IW-RTC model resulted in the drop of the mean height residual with regard to IF from 30 to 22 mm and from 32 to 29 mm for RTK and rapid static modes, respectively. If we inspect Table 3, we discover a foreseeable deterioration of the float solution precision during high ionospheric activity. For example, in the case of IW-RTC RTK, coordinate STD of the float solution during high ionospheric activity dropped approximately 3-4 times as compared to the undisturbed period. Fortunately, no deterioration was observed after correct ambiguity fixing and this solution provided a similar position accuracy for both disturbed and quiet periods.

Conclusions
This paper introduced and assessed the methodology aiming at a reliable wide-area multi-GNSS relative positioning under the presence of ionospheric disturbances. The enhanced approach takes advantage of both multiconstellation network ionospheric corrections and the RTC algorithm which eliminate the temporal variations of the ionospheric delay. As a consequence the residual DD ionospheric delays are reduced thanks to the application of the network corrections and, what is more important, may be estimated as constant parameters for the entire DD observational arc taking benefit from RTC.
The experimental evaluation of the methodology was performed on the basis of RTK and rapid static positioning with GPS, BDS and Galileo data collected at high latitudes during the ionospheric storm on August 25-26, 2018. The performance of the enhanced model was compared to the benchmark solutions, which in this case were the ionosphere-float and the ionosphere-weighted models.
The results confirmed the deterioration of the accuracy of the network ionospheric corrections and consequently a decline of the positioning performance with regular ionosphere-weighted model during active ionosphere. We may conclude that the support from network ionospheric corrections under the presence of the ionospheric disturbances is not enough for reliable precise positioning due to poor quality of these corrections. On the contrary, the results obtained with the application of the proposed methodology demonstrated a distinctive improvement in the ambiguity resolution domain and thus proved its superiority over the ionosphere-float and ionosphere-weighted models. More specifically, a clear advantage of the enhanced methodology over benchmark solutions was confirmed by the values of ambiguity success rate during high ionospheric activity. The IW-RTC ASR was about 20% higher with respect to that of the benchmark models. A similar level of improvement was observed for the rapid static mode, since after 30 epochs the ASR reached the level of 74.4% for the IF and IW models and 94.0% for IW-RTC.
Eventually, the IW-RTC model proved also better performance than IF and IW models in the coordinate domain. The most noticeable improvement was reflected in the precision of the height component. For example, in RTK mode this statistics dropped from 92 mm for IF down to 54 mm for the IW-RTC model during high ionospheric activity.