WO2016185500A1 - Method for forecasting ionosphere total electron content and/or scintillation parameters - Google Patents

Abstract

The present invention relates to a method of TEC and scintillation empirical forecasting, in particular short-term forecasting (seconds to minutes). The output of the method is necessary to feed mitigation algorithms aiming at improving accuracy on GNSS precise positioning techniques (RTK, NRTK, and PPP) under ionospheric harsh conditions.

Description

METHOD FOR FORECASTING IONOSPHERE TOTAL ELECTRON

CONTENT AND/OR SCINTILLATION PARAMETERS

The present invention concerns a method for forecasting ionosphere total electron content and/or scintillation parameters.

More in detail, the present invention relates to a method of TEC (Total Electron Content) and scintillation empirical forecasting, in particular short-term forecasting (seconds to minutes) . The output of the method is necessary to feed mitigation algorithms aiming at improving accuracy on GNSS precise positioning techniques (RTK, NRTK, and PPP) under ionospheric harsh conditions. See at the end of description for acronyms .

State of the art

The ionosphere is the single largest contributor to the GNSS error budget. Although the bulk of its effect on the propagation of GNSS signals can be generally modelled to a first order using multiple frequency range measurements and special algorithms, its state can be very erratic, depending on location, season, local time and prevailing solar activity (see e.g. Kim and Tinin, 2011, and references therein) . More often but not only around the peak of the 11-year solar cycle the ionosphere may become so disturbed as to lead to severe satellite signal degradation, affecting in particular real time high accuracy carrier phase based techniques such as RTK, NRTK and PPP.

This is due to the rapid fluctuations in the amplitude and phase of radio signals from GNSS satellites as they pass through small-scale plasma density irregularities in the disturbed ionosphere as sketched in Figure 1 (Wernik & Liu 1974; Kintner et al. 2001, 2009 and the references therein) . These rapid fluctuations are called "ionospheric scintillations" and are worldwide represented by two indices S4 (amplitude scintillation) and σ,, (phase scintillation) . To characterize scintillation other quantities are also used, i.e. the opposite of the slope (often noted "p") of the phase power spectral density in log-log axes, and the spectral strength (often noted "T") that is the detrended phase power spectral density at 1.0Hz (see for details Wernik et al., 2007).

The effects of a disturbed ionosphere are more prominent in regions at high and equatorial latitudes and this is well depicted by Basu et al (2002) in Fig. 2 where the occurrence of scintillation phenomenon is shown on the global scale for high solar conditions (a) and low solar conditions (b) in geomagnetic coordinates. The mechanisms responsible of the polar and equatorial scintillation are quite different each other from a physical point of view but their effects on GNSS based applications are equally disruptive. In Fig. 3, an example of the effects of strong scintillation is shown on the positioning at high latitude: phase scintillation has severely affected the 3D positioning accuracy on 30 October, 2003 at Ny- Alesund (78°55'N, 11°56'E), as shown in Alfonsi et al. (2006) .

Together with the scintillation indices, the ionospheric TEC and its sudden variation are also important indicators of the ionospheric scenario that may occur: polar patches at high latitudes and plasma bubbles at low-equatorial latitudes are TEC sudden fluctuations (in time and space) often accompanied by scintillation (see e.g., De Franceschi et al 2008; Alfonsi et al, 2011) . Moreover, these sudden variations are quite relevant as they exacerbate the so-called "ionospheric delay" that can adversely affect the GNSS signal propagation through the ionosphere. In this case, ionospheric delay cannot be fully eliminated by a linear combination of observables on two frequencies such as GPS LI and L2. The scintillation scenarios also depend on the longitudinal sector. For example, the Brazilian ionosphere can be considered as one of the most affected area where scintillation events can be severe, so Brazil represents the worst-case scenario for disruption to real time high accuracy GNSS based applications. In Brazil, ionospheric disturbances show a strong seasonal dependence (with equinoctial and summer solstice presenting stronger effects) during periods within and beyond the peaks of the solar cycles (Akala et al. 2011 and Muella et al. 2013) . Moreover, in Brazil the scintillation is a daily event mostly confined during post sunset (22 UT - 04 UT) .

In the art, these problems do not find a satisfactory solution, and problems remains as follows.

Observed scintillation activity is characterized by a considerable spatial and temporal variability, which depends on factors such as the frequency, zenith angle or angle between the ray path and the Earth's magnetic field. The effect of these factors can be accurately defined based on the scintillation theory. However, scintillation dependencies on local time, season, solar and magnetic activity have a stochastic character, meaning that there is no unique relationship between the strength and/or occurrence of scintillation and the particular agent. That is why it is so difficult to forecast the occurrence of scintillations, and therefore predict the impact of ionospheric disturbances on radio communications, navigation or positioning systems.

Actual models able to predict ionospheric scintillation can be divided into three families: analytical models, global cliraatological models and climatological models based on in situ data.

Most known analytical models are cited in Priyadarshsi et al., 2015. Most of these models suffers from the limitation of being limited to the geographical sectors and/or frequencies (historical modelling was for VHF not for GNSS frequencies) for which they have been developed (Priyadarshsi et al., 2015) . The most known and used climatological models are the WBMOD (WideBand MODel, Secan et al., 1995) and the GISM (Global Ionospheric Scintillation Model, Beniguel and Buonomo, 1999 and Beniguel, 2002) . Because of the climatological nature of such models, they cannot properly account for the actual conditions but rather give an average dependence of the propagation characteristics on geophysical conditions. As an example, the patchy character of the irregular structure of the low-latitude ionosphere is completely absent (Priyadarshsi et al., 2015), so they often fail in catching the detailed ionospheric morphology needed to feed mitigation algorithms aiming at improving accuracy on GNSS precise positioning techniques. Most known climatological models based on in-situ data are the Basu et al. model (1976) and the WAM model (Wernik et al., 2007). Both models have the limitation to be climatological ones and are limited in space and time by the in-situ data used in its construction, so their outputs is not useful to improve GNSS precise positioning accuracy.

Concerning the TEC parameter, to-date there is no method for forecasting it, to the knowledge of the Inventors .

It is the object of the present invention to provide a method and a system for forecasting ionosphere parameters which solves the problems and overcomes the drawbacks of the present invention.

It is subject-matter of the present invention a method and a system according to the enclosed claims.

Once forecasted values are obtained, they are used to feed prior art algorithms able to mitigate the ionospheric effects on precise positioning (e.g. Aquino et al, 2009 and reference therein) .

The invention will be now described by way of illustration but not by way of limitation, with particular reference to the drawings of the enclosed figures, wherein: - Figure 1 shows a scheme of the varying effects of scintillation on GNSS: electromagnetic wave travelling along a given raypath from satellite to receiver is influenced when passing through ionospheric irregularities, (from Kintner et al.,

2009) ;

- Figure 2 shows the occurrence of scintillation at solar maximum (a) and solar minimum (b) ; scintillation is most intense and most frequent in two bands, surrounding the magnetic equator, and at poleward latitudes (from Basu et al., 2002);

- Figure 3 shows (a) Phase scintillation index σ_Φ and (b) 3D positioning error measured on 30 October 2003 at Ny-Alesund, Norway (adapted from Alfonsi et al. 2006).

- Figure 4 shows a box diagram of the invention;

- Figure 5 shows temporal flow chart representing the idea of estimation of parameters and making the forecast. Present is denoted as tO (grey horizontal line), past N measurements are represented by thin horizontal black lines separated by sampling interval At and forecasting (thick horizontal black line at the bottom) is separated from the present by the forecasting horizon h, according to the invention;

- Figure 6 shows a sketch of adjacent triangles geometry and quantities, as calculated in a step of the method according to the invention;

- Figure 7 shows vTEC data used in the test run of the model according to the invention. Each line represents a satellite-receiver link;

- Figure 8 shows Delaunay triangulation of the computational domain with internal points numbering, in a test example according to the invention;

- Figure 9 shows a spectrum of the SVD values (from 1 to 112) plotted in decreasing magnitude;

- Figure 10 shows reconstructed velocity field in an application example; the light grey arrow shows velocity vector with magnitude of 100 m/s;

Figure 11 shows comparison of the forecasting (black circles) with measured vTEC data (grey circles) as a function of time for pierce point number 24 (for its position see Figure 8);

Figure 12 shows S4 indices corresponding to the data shown in Figure 7. Each line represents a satellite-receiver link;

Figure 13 shows a spectrum of the SVD values for the S4 forecasting, according to the invention; Figure 14 shows an example of S4 1-minute forecasting for each piercing point (x-axis) : input S4 (triangle), i.e. the initial condition at TO, actual S4 (cross) at TO+1 min, and forecasted S4 (plus) for T0+lmin, according to the invention; Figure 15 shows time profile (left plot) and corresponding distribution (right plot) of the difference between actual and forecasted value of TEC for the day 26/09/2013 according to the invention, box indicates where the ionospheric effects are expected to exacerbate;

Figure 16 shows RTK positioning errors (in meters) in North-South direction (dN, upper panels), East- West direction (dE, middle panels) and Up-Down direction (dU, bottom panels) obtained using IGS TEC map (left panels) and forecasted TEC map (right panels) for the day 26/09/2013 according to the invention;

Figure 17 shows a schematic block diagram of a positioning system according to the invention.

Detailed description of the embodiments of the invention

The past value of ionosphere parameters are known at points called ionospheric pierce points (IPP), which are intersections of the ray path of satellite signal between a certain receiver and a satellite with ionosphere at a certain height. The IPPs have its coordinates at the geographic coordinate system that we denote by r_k for the k-th. IPP.

Making reference to Fig. 4, the data needed for next steps (i.e. values of TEC and scintillation parameters as well as positions of corresponding IPPs) of the algorithm are acquired at initial stage 100. The flow chart of Fig. 4 shows two calculation channels. Both channels have separate input data flows (102 for TEC and 104 for scintillation parameters respectively) . TEC input data usually requires known calibration procedure 103 to eliminate satellite and receiver biases from GNSS observables.

The forecasting method distinguishes two channels too: one for TEC and the other for scintillation parameters. The channels differ in description of the quantities to be forecasted. In the TEC channel, the continuity equation in the conservative form is used while description of scintillation parameters uses the continuity equation with source term added, as described in the following. Both channels starts with the triangulation procedure 200 that gives a structure to the computational domain used in the subsequent steps of the algorithm.

Either scintillation parameters and TEC parameter can be provided alone to the mitigation algorithm, obtaining a mitigation on the GNSS prediction error. Of course, providing both sets of forecasted parameters will result in a much better mitigation.

TEC channel

The method exploits the transport theory for a scalar field. The basic equation of continuity for a scalar, provided the velocity field v is known, and is:

in a volume V with boundary dV .

The technical concept of the method is to reconstruct the velocity field v from TEC measurements and subsequently evolve the scalar field / according to equation (1) with the desired time resolution using an appropriate numerical integration scheme. In our case, the velocity field v is the velocity of the integrated electron density (TEC) along the line of sight connecting the receiver and the satellite, while / is the scalar field of TEC. The reconstruction of the velocity field v is performed by fitting it to the time changes of TEC field using recent data (i.e. using TEC data to derive velocities by equation (1))). Special hypothesis can be made for the values of scalar and velocity fields at the boundary of the volume (e.g. one can suppose that the flux at the boundary of the total volume is zero) . It is here to be understood that the volume is a generic name for a element of space, therefore if the equation is integrated over a surface, V will be a surface area.

According to an aspect of the invention, performing the forecasting for TEC channel may consist of the following steps:

1. discretizing the space into Delaunay triangles, hereafter "triangulation" (block 200 on the method flow chart);

2. forming the approximation to the equation (1) as a budget (i.e. computation of the net flux) of conserved quantity (TEC) between triangles forming triangulation, approximating TEC piecewise linearly on the triangular grid and considering velocity constant over a triangle (301 on the method flow chart);

3. performing SVD of the resulting matrices (block 401 on the method flow chart) ;

4. solution for the velocity field is sought by applying suitable pseudo-inverse operator to the temporal changes of the field f (block 501 on the method flow chart) ;

5. making the forecasting step by applying suitable time integration method to the equation of continuity keeping the velocity field constant over the forecasting horizon (TO to T0+horizon, block 601 on the method flow chart) .

The above steps are further detailed in the following.

1. Triangulation domain is performed by suitable triangulation algorithm (for example Delaunay triangulation) . The output of this step of the algorithm is a set of triangles

Δ_k

where K is the total number of triangles. The underlying technical concept of this step is assimilating the region of the ionosphere, where in the scalar field is to be calculated, to a plane (at an altitude value chosen in the range of 350-400 km typically) , and triangulation keeps locality of the solution of the general equation.

2. The triangulation Δ allows to discretize equation (1) into the form:

wherein we assume the velocity v_k as constant in the triangle Δ_k . In fact, this formula holds also for any type of discretization, be it by squares or other geometrical forms, and even with discretization by areas which are not equal to each other.

The vertices of the triangles are the pierce points of the above plane of the ionosphere with the path between transmitter and receiver. Let us denote the Δ_k ' s triangle vertices by: r_kl, r_k2, r_k3 , orientated counter-clockwise and the corresponding values of a scalar / at these points: f_kv f_kZ, f_k3 . To the edge- vectors (vectors of the sides of triangles) we assign the index of the opposite vertex, by u_kl = r_fe3— r_k2,.... The function /, representing the vTEC scalar field, is then approximated piecewise linearly over the triangulation. The velocities on triangles A

A_l2, adjacent to A

on its respective sides 1,2,3, form the set v_llt v_l2, v_a. The geometry of adjacent triangles and quantities used in the following formulation are shown in Fig. 6.

At a given time t, the total flux

£ for each triangle is given by the sum of the fluxes Φ

across each i-th triangle boundary that can be approximated as the sum of contributions from adjacent triangles :

where u_ki is the vector perpendicular to the vector

(the flux is assumed as positive when going out from the triangle) , and is the sum of the scalar at

points forming particular vector at time instant t.

The contribution of each in equation (3) could be

rewritten as:

where superscripts t and t + At denote consecutive time instants separated by the time interval At (the values of

the time series before the starting point of the forecast are used in this velocities calculation step) ,

is the standard measure (area) of triangle Δ_k and

we can write (4) as:

where in . Such equations are composed for

the whole area at N subsequent time instants (from to giving the matrix equation:

where M is a block matrix:

Matrix M has KN rows and 2K columns (as vectors u_ki have two components: u_kix and u_kiy ) . The corresponding column vector Af has KN entries:

and v is a column vector with 2K entries;

3. Matrix M could be factorized by using the so called SVD (Strang, 1998):

wherein V is a unitary 2K X 2K matrix (K is the number of triangles in the domain), U is the unitary NKxNK matrix (N is the number of time levels (At's) of the vTEC field to which the velocity field is fitted) and ∑ is a diagonal NKx2K matrix of singular values

indicates the complex conjugate)

4. Using SVD, we can apply the so called pseudo- inverse to obtain regularized solution to equation (7) :

Here

is a diagonal matrix with entries

The SVD regularization scheme replaces with zero the small singular values, i.e. small

as they introduce large errors in the solution compromising the stability of the solution to equation (7) . This is possible only when the SVD spectrum presents a discontinuity in the When the discontinuity is present, the behaviour of the field is fluid-like, and the continuity equation does not need extra terms (i.e., production and loss terms) . See below for an example (Fig. 9) . When this is not possible, a source term has to be added (see scintillation case) .

5. To make the forecast we use the equations (6) with estimated velocity field v (keeping it constant over a forecasting horizon h, see Fig. 5) and integrate in time the resulting equations:

Although steps 1-5 are only an aspect of the invention, they are particularly convenient in that the locality of solution of (1) is kept and the approximations are such that the calculation is fast and reliable. It is to be stressed that the invention method is realisable only when there is a suitable infrastructures of GNSS receivers, so that we have sufficient data to solve the equations. Prior art method do not exploit such infrastructure but model the ionosphere phenomena as such. Scintillation parameters channel

The forecasting of scintillation parameters channel consists essentially of the same steps as for the TEC channel except in the 2nd step we modify the continuity equation to the form:

wherein i.e. the "source term", is the density of

production (p) minus losses (/) rate. So, now the forecasting for the scintillation parameters channel can be casted into following steps:

1. discretizing the space into Delaunay triangles (block 200 on the method flow chart) ;

2. form the approximation to the equation (15) with the additional source term as a budget of a quantity between triangles forming triangulation approximating scintillation parameters piecewise linearly on the triangular grid and considering velocity and the source term constant over a triangle (block 302 on the method flow chart);

3. performing SVD of the resulting matrices (block 402 on the method flow chart) ;

4. solution for the velocity and the source term fields is sought by applying suitable pseudo- inverse operator to the temporal changes of the field / (block 502 on the method flow chart) ;

5. making the forecasting step by applying suitable time integration method to the equation of continuity keeping the velocity field and the source term constant over the forecasting horizon (block 602 on the method flow chart) .

1. Triangulation domain is performed by suitable triangulation algorithm (for example Delaunay triangulation) . The output of this step of the algorithm is a set of triangles 4

where K is the number of triangles.

2. When we add the source term for each triangle of the triangulation, Equation (3) changes to:

where ¾ is the source term for the k-th triangle. Now we can write the system of equations for estimating both the velocity field v and the source term in a compact form as:

where M is the block matrix (similarly to the case TEC) :

in which each block is now:

At the k-th position now the term μ(Δ_¾) appears which is the area of the k-th triangle, since the Equation (16) includes also the source term contributing to the total flux of the scalar field and related to the k-th triangle. The vector s takes the following form:

3. Matrix M could be factorized by SVD (Strang,

1998) :

wherein V is a unitary 3K X 3K matrix (K is the number of triangles in the domain), U is the unitary

matrix (N is the number of time levels, Δt's, of the scintillation parameters field to which the velocity field and source term are fitted) and ∑ is the diagonal NKx3K matrix of singular values

4. Using SVD, we can apply the so called pseudo- inverse to obtain regularized solution to equation

(17) :

Here ∑⁺ is a diagonal matrix with entries

The SVD regularization scheme replaces with zero the small singular values, i.e. small σ_i, as they introduce large errors in the solution compromising the stability of the solution.

5. To make the forecast, we use estimated velocity field v and source terms n_k (as a function of time, i.e. keeping them constant over forecasting horizon h) and integrate in time the resulting equations:

The above algorithm, also in its general form, is suitable for short term prediction, if the values of f are forecasted on the basis of an acquired time series of values of /. If instead the predicted values are added to time series and the new series is considered for the forecasting of a subsequent time value of /, then the forecasting can be extended for longer term.

Example of calculation for a small area

Example of computation

Test of the model run is here shown for both the TEC and scintillation parameters channels. Test has been conducted on the data taken on the 1 November 2011, a day characterized by strong scintillation regime. Data have been acquired by the CIGALA/CALIBRA network of GNSS receivers for scintillation (http: //is- cigala-calibra . fct . unesp.br/is/index . php) .

CIGALA/CALIBRA network is owned by the Brazilian Universidade Estadual Paulista "Julio de Mesquita Filho" and a summary of the receiver location, name and geographic coordinates of the scintillation receivers used in the model tests is in Table 1.

Table 1. List of the name, location and geographic coordinates of the scintillation receivers used to test the model.

TEC channel example

Fig. 7 shows the TEC data mapped to equivalent vertical value (vTEC) used for computations (model input) in which colours denote records for different pairs satellite-receiver. Computed ionospheric pierce points form our computational domain for which Delaunay triangulation (Delaunay, 1934) was created. The triangulation was formed from 33 points and has 56 triangles. The triangulation points are then divided into sets of boundary and internal points respectively (the triangulation with numbering of internal points is shown in Fig. 8) .

To reconstruct the velocity field, model was set to use 10 time levels (from To (starting point) to To+10min (final point)), with a 1 minute resolution. The SVD spectrum of the matrix M is depicted in Fig. 9, in which the logarithm of the singular values of the M matrix are plotted in decreasing magnitude. There is a notable jump after the singular value no. 56, which also corresponds exactly to the number of triangles in the computation domain. The TEC velocity field resulting from the regularization of the solution is shown in Fig. 10.

As an example of the forecasting results, comparison between forecasting (black circles) with the measured vTEC data (grey circles) as a function of time for a selected pierce point is shown in Fig. 11. The pierce point number is the 24 and their position is in Fig. 8. Forecasting horizon is 20 minutes.

Scintillation parameters channel example

As an example of forecasting of scintillation parameters, here we show the model test of S4 forecasting for the same day of the TEC forecasting, i.e. 1 November 2011, by using the same CIGALA/CALIBRA network data. Fig. 12 shows the S4 data used for computations (model input) in which colours denote different satellite-receiver records and that correspond to the vTEC value presented in Fig. 7.

The SVD spectrum of the matrix M, which now includes the source term, is depicted in Fig. 13, in which the logarithm of the singular values of the M matrix, σi , are plotted (from 1 to 141) in decreasing magnitude. There is a notable jump after the singular value no. 47, which also corresponds exactly to the number of triangles in our spatial domain. We keep only these highest 47 singular values in creating pseudo inverse that is the standard procedure in regularizing ill-posed matrices by SVD (see Strang, 1998) . The total number of singular values (/ = 141), is 3 times the number of triangles due to the contribution of the two components of the velocity field plus the source term for each triangle. The presence of the jump in the spectrum means that now the rank of the matrix M can be well defined.

As an example of S4 1-minute forecasting, Fig. 14 shows for each piercing point the input S4 (triangle), i.e. the initial condition at TO, measured S4 (cross) at TO+1 min, and forecasted S4 (plus) for T0+lmin.

Test and statistical validation under strong scintillation scenario

The model prediction capability is tested by means of the difference (Δ₀) between the actual (real) and the forecasted (modelled) values (Q) for each of the forecasted quantities, i.e. TEC and scintillation parameters (S4, σ_φ ,p and T, see e.g. Bougard et al., 2011) . Statistical analyses on the five Δ₀ variables have been carried out to assess the overall performance. The standard deviation (σ_ΔQ) of each Δ_Q distribution, obtained by considering GNSS data for 5 days under severe scintillation conditions, is the indicator of the average performance of the forecasting model. An example is in Fig. 15. The overall performance of the prediction is summarized in the Table 2.

Table 2. Error in the forecasted values of TEC and ionospheric parameters for different confidence levels.

In the Table 2, P indicates the confidence level or probability (here from 68% to 99% ) , associated to a given σ_ΔQ which is assumed to be the forecasting error. In the table, such error varies (depending on the desired confidence level) from about one order of magnitude in the case of S4, σ_φ, p, TEC to about 2 orders of magnitude for the T parameter. Considering typical values of the parameters during strong scintillation events and 99% of confident level, the relative error associated to the forecasted TEC is less than 5 % and the one associated to the scintillation parameters is about 10% - 15%.

GNSS precise positioning techniques and mitigation algorithms General description of GNSS precise positioning techniques

GNSS carrier phase based positioning techniques can provide a much higher precision and accuracy than their code based counterparts, and nowadays represent the high accuracy GNSS positioning techniques of choice. They have been widely used to support many high accuracy applications such as precision agriculture, constructions, land management and geodesy/land surveying.

The main carrier phase based positioning techniques, namely RTK/NRTK and PPP, will be briefly introduced and reviewed with particular emphasis on the impact of ionospheric TEC disturbances and scintillation. A detailed description of the main precise positioning techniques can be found in Yang et al. 2013 and references therein.

Real Time Kinematic (RTK)

RTK is one of the most widely adopted GNSS high accuracy positioning techniques . With the support of a nearby reference station (RS) with a known location, and a reliable wireless communication system, the common measurement errors/biases between the mobile "rover" (the "moving" station for which the precise position is needed, e.g. mounted on a car, tractor or plane) and the RS can be cancelled through differencing, and the rover positioning solution can reach the level of a few centimetres accuracy in realtime when using the carrier phase measurements. However, in traditional single RS based RTK positioning, the rover-RS distance is limited to within 10 to 20km. Beyond this boundary, the distance- dependent errors, which mainly consist of atmosphere related errors, may increase rapidly and lead to accuracy degradation. To solve this problem, Network based RTK (NRTK) was developed in the late 1990s as an evolutionary technique. It utilizes the measurements gathered from a network of Continuously Operating Reference Stations, interpolates and mitigates the distance-dependent errors at the rover location. NRTK can increase significantly the rover-RS distance whilst maintaining accuracy, depending on the local ionosphere and troposphere variation conditions. However, the ionospheric condition in equatorial regions, i.e. in Brazil, is very active so that even short distance of only 20km may collapse the inter-station correlation (Galera et al., 2013).

Precise Point Positioning (PPP)

Precise Point Positioning (PPP) has studied since late '90's as an alternative technique to RTK/NRTK and it is also based on carrier phase measurements. Combining the rover dual-frequency carrier phase measurements with external precise satellite orbits and clocks products, PPP is able to provide position solutions at centimetre to decimetre level, even less than 1 cm-level positioning in static mode.

Differently from RTK/NRTK, PPP uses a zero- difference technique, which does not require access to the observations of one or more reference stations with known coordinates. Zero-difference brings some advantages to PPP. First, the spatial operating range limitation of the differential techniques is overcome, as well as the need for simultaneous observations at both user and reference receivers. This leads to reduced infrastructure labour and equipment costs, as well as simplified operational logistics (Gao, 2006) . Second, PPP provides an absolute positioning solution in a global reference frame, instead of the relative solution to the reference station in a local reference frame, which reduces the positioning reliance on the reference station coordinates quality (Haasdyk and Janssen 2012) . Third, PPP only requires precise orbit and clock information, which can be computed from the measurements in a relatively sparse reference station network (typically hundreds of km inter-station distance) and theoretically, the same orbit and clock product can serve all PPP users distributed globally. The main criticism to it is the long convergence time required to achieve centimetre level performances (typically in the order of tens of minutes) . This largely jeopardizes productivity, and is the bottleneck for PPP to be applied in real time. Second, but important for the invention application, because PPP does not use a differencing technique, no ionospheric delay can be relatively cancelled out. Therefore, dual- frequency measurements are required to remove the first order effect of the ionosphere. However, if the ionosphere TEC disturbance and scintillation is strong, the ionosphere noise will dramatically increase (as largely described in section 1) extending the PPP convergence time.

Strategies to mitigate the ionospheric effect on GNSS precise positioning techniques

There are different approaches that can be used for mitigation using independent ionospheric information. These approaches can be divided into two main groups, one is mitigating scintillation at the observable level and the other at the positioning level. In the former group, the main idea is to provide 'clean observables' to the positioning solution, based on externally monitored or predicted scintillation information. That is informing users if their observations are contaminated by scintillation, which observations are involved, and how severe the situation is (see Fig. 1) . The user can then decide how to use these observations inside PVT (position, velocity and time) engine at GNSS receiver level. This group includes the screening and the weighting of contaminated observations. In the latter group, i.e. the mitigation at the positioning level, the main idea is to adapt the PVT engine itself, applying different strategies/parameters according to the monitored/forecasted scintillation/TEC level. The approaches in this group include cycle slip detection/correction and ambiguity resolution with an adaptive critical value.

Both these mitigation strategies take advantage from the forecasting model proposed in the invention as the forecasted ionospheric parameters can be used for the real-time applications of the mitigation algorithms. For example, scintillation forecasted parameters can be used to "clean" the observables from GNSS before using them to compute position (see above) . Moreover, forecasted TEC values can be used to estimate ionospheric delay to be ingested by the "rover" receiver in order to compute its precise position (see below example) . Screening

The simplest and most instinctive approach to mitigate scintillation effects at observable level is satellite screening. In the satellite-screening scheme, the scintillation indices at the rover location are monitored and checked against a pre-defined threshold. The observations with high indices values are isolated and screened out from the receiver PVT engine.

Weighting

Different from screening, in the weighting scheme the contaminated observations are utilized, rather than excluded, from the PVT solution. Each observation will be weighted in the positioning engine based on a scheme that depends on its scintillation level. If the weight assigned is really low, the observation will be effectively screened out from the positioning solution. The weighing scheme can be achieved by using the variance of the phase/code tracking jitter of specific observations (Aquino et al, 2009) . In case of scintillations, the measurement variance (the tracking jitter variance) can be estimated from the Conker model (Conker et al., 2003), which requires scintillation information as input. If the scintillation can be predicted for an arbitrary user location, the user will be able to apply the Conker model to compute the measurement variances for the PVT engine. Example of application of the invention on GNSS precise positioning

In RTK processing, atmospheric effects including ionospheric and tropospheric delays are cancelled out by the differencing techniques if a rover receiver is sufficiently close to a reference station (i.e., the atmospheric conditions affecting both receivers are similar) . But in long baseline RTK, which is the case of the network available in Brazil, and in PPP, external information about the atmosphere, the ionosphere in particular, is required, especially if the changeable conditions encountered in equatorial regions are considered. The ionospheric delay, one of the most significant GNSS error sources, can be quantified by the TEC due to their linear relationship. Thus, accurately determined TEC can be converted to the ionospheric delay on a GNSS ray path. In addition, the TEC gradient, either spatial or temporal, is also useful ionospheric information for positioning. Therefore, a regional TEC and TEC gradient map can represent valuable a priori information to assist GNSS positioning.

To show the improvement that is obtained by using the invention output to feed mitigation techniques, a test on a long baseline RTK for two stations in the Sao Paulo state (115 km) has been performed (Park et al., 2015) . Two different TEC maps have been used to feed the RTK algorithm: the one obtained by the model object of the invention and the one provided by IGS (GIM map, http://iono.jpl.nasa.gov/gim.html). Data for the selected baseline was processed on a day when strong scintillation was observed (29 September 2013) .

It is clear from Fig. 16 that the use of the maps obtained by the method subject-matter of the invention, instead of the one available from IGS, dramatically reduces the positioning errors in long-baseline RTK technique: 3D RMS calculated over the whole day of observations reduces from 1.4 m to 0.4 m. System of GNSS positioning according to an aspect of the invention

Making reference to Fig. 17, a GNSS positioning system 100 for the positioning of a rover 20 in an interest area 60 comprises:

- a GNSS network comprising scintillation and total electron content monitor receivers 10,

- a central unit 30, comprising a computer logic having a computer program installed on it and configured to execute the method of the invention;

- a communication network 40 between said central unit 30 and said receivers 10;

- a communication network 50 between said central unit 30 and said rover 20.

The central unit may be provided with means configured to execute the mitigation above described.

7. Fields of applications of GNSS precise positioning techniques

7.1 Precision agriculture (PA)

Although the concept of PA was introduced at the beginning of the last century, the technology available so far was not able to put it in practice. With the goal of increasing the field management yield, for almost two decades the concept of PA has evolved and is now treated as an agricultural practice which uses technology information based on the principle of variability of soil and climate. Based on specific geo- referenced data, it deploys the process of agricultural automation, dosing fertilizers and pesticides differently. The technological elements that highly contributed to the development of this concept were microprocessor and GNSS, which are coupled and integrated on the harvesters, seeders and other implements, allowing the data collection, cumulative tabulation and dosed application with localized inputs. Within this context, one of the fastest growing segments is the localized application of pesticides, and fertilizers as well as the autopilot that enables work for planting and harvesting every 24 hours. However, in the equatorial region, problems in GNSS signal quality, due to interference from the ionosphere, have impeded the operation as recommended by manufacturers. The result has been the complete stoppage of the machines during certain periods of the year and at certain times (post-sunset hours) of day causing considerable economical losses and discrediting to the technology.

7.2 Civil Aviation

In support to air navigation, reference stations must have their coordinates determined very accurately. It is also necessary an effective control of system operation with immediate warning information to users in case of any failure (integrity) .

To meet these requirements, complementary systems to GNSS have been designed and are currently in use. The International Civil Aviation Organization has classified these systems into two types, which can be applied in isolation or combined:

1. SBAS (space based augmentation systems, such as EGNOS and WAAS - cover large areas) this system is used to complement other satellite systems, e.g., GPS and/or GLONASS;

2. GBAS (ground based augmentation systems - local) - provide localized support such as in the vicinity of airports.

Due to the extreme ionospheric conditions observed in equatorial regions (i.e. in Brazil), the use of GBAS correction is encouraged because the SBAS is ineffective.

An example is the GaleSo International Airport in Rio de Janeiro where a GBAS system purchased from Honeywell Aerospace is in the process of certification, which was installed and is currently being tested. If the system meets the expectations of the authorities, it should be adopted in other airports in Brazil.

Therefore, new algorithms that enable GNSS high accuracy positioning techniques and systems to mitigate the effects of the ionosphere, like the forecasting invention method, may contribute to improving the scenario for the use of GNSS and SBAS (EGNOS) in Brazilian civil aviation.

In conclusion, no method besides the proposed invention is currently able to give short-term forecast (from seconds to minutes) or longer forecast of the basic quantities describing ionospheric and signal propagation conditions in particular under severe scintillation environment. These quantities are necessary to feed mitigation algorithms to improve the performance of the real-time GNSS precise positioning techniques. The present invention method provides a deterministic forecasting, contrary to the traditional models based on statistical approaches: the computation of the velocity of the scalar field plays a keyrole in the parameters forecasting.

List of Acronyms

References

Akala, A. 0., et al. "Ionospheric foF 2 variability at equatorial and low latitudes during high, moderate and low solar activity." Indian Journal of Radio & Space Physics 40.3 (2011): 124-129.

Alfonsi L, De Franceschi G, Romano V, AquinoM and Dodson A (2006) Positioning errors during the space weather event of October 2003. Location, ISSN 0973- 4627, 1, issue 5.

Alfonsi, L., Spogli, L., Tong, J. R., De Franceschi, G., Romano, V., Bourdillon, A., ... & Mitchell, C. N. (2011) . GPS scintillation and TEC gradients at equatorial latitudes in April 2006. Advances in Space Research, 47(10), 1750-1757.

Aquino, M. , et al. "Improving the GNSS positioning stochastic model in the presence of ionospheric scintillation." Journal of Geodesy 83.10 (2009): 953- 966.

Basu S, Basu S, Khan BK (1976) Model of equatorial scintillation from in situ measurements. Radio Sci 11:821-832

Basu, S., K.M. Groves, Su. Basu and P.J. Sultan, 2002, Specification and forecasting of scintillations in communication/navigation links: current status and future plans, J. Atmos . Solar-Terr. Phys. 64, 1745-1754. Beniguel Y (2002) Global ionospheric propagation model (GIM) : a propagation model for scintillations of transmitted radio signals. Radio Sci 37:RS1032. doi: 10.1029/2000RS002393

Beniguel Y, Buonomo S (1999) A multiple phase screen propagation model to estimate fluctuations of transmitted signals. Phys Chem Earth (C) 24:333-338

Bougard, B. , Sleewaegen, J. M. , Spogli, L., Sreeja, V.

V., & Galera Monico, J. F. (2011). CIGALA: challenging the solar maximum in Brazil with PolaRxS. Proceeding of the ION GNSS, 13231332.

Delaunay B. (1934), Sur la Sphere Vide. Bullettin of the Academy of Sciences USSR, VII, 793-800.

Conker, R. S., M. B. El-Arini, C. J. Hegarty, and T.

Hsiao, Modeling the effects of ionospheric scintillation on GPS/ Satellite-Based Augmentation

System availability, Radio Sci., 38(1), 1001, doi:10.1029/2000RS002604, 2003.

De Franceschi, G., L. Alfonsi, V. Romano, M. Aquino, A. Dodson, C.N. Mitchell, P. Spencer, A.W. Wernik (2008), Dynamics of high-latitude patches and associated small- scale irregularities during the October and November 2003 storms, J. Atmos. Solar-Terr. Phys, 70(6), 879- 888, doi: 10.1016/j . jastp.2007.05.018 Galera Monico

F. , Paulo de 0. Camargo, Vinicius A. Stuani Pereira, Bruno C. Vani, Daniele B. Marra Alves, Marcio Aquino (2013), CALIBRA Dl.l, User performance review. (http: //www. calibra- ionosphere.net/calibra/documents/calibradll .pdf)

Gao Y. (2006) Precise Point Positioning and its challenges, Inside GNSS, 1(8):16-18.

Haasdyk, J. and V. Janssen (2012), Choosing the best path: Global to national coordinate transformations, Coordinates, 8(2):10-16.

Kim, B. C, and M. V. Tinin. "Potentialities of multifrequency ionospheric correction in Global Navigation Satellite Systems." Journal of Geodesy 85.3 (2011) : 159-169

Kintner, P. M., H. Kil, T. L. Beach, and E. R. de Paula (2001), Fading Timescales Associated with GPS Signals and Potential Consequences, Radio Science, 36(4), 731- 743.

Kintner, P. M. , T. Humphreys, and J. Hinks. "GNSS and ionospheric scintillation." Inside GNSS 4.4 (2009): 22- 30.

Muella, Marcio TAH, Eurico R. de Paula, and Alan A. Monteiro. "Ionospheric scintillation and dynamics of Fresnel-scale irregularities in the inner region of the equatorial ionization anomaly." Surveys in Geophysics 34.2 (2013) : 233-251.

Park J., S. Veettil, M. Aquino, L. Yang, C. Cesaroni , Mitigation of Ionospheric Effects on GNSS RTK at Low Latitudes, accepted for presentation at ION GNSS, Tampa (Florida, USA) September 2015

(https : //www. ion. org/gnss/sessions . cfm?sessionID=398) Priyadarshi, S. (2015) . A Review of Ionospheric Scintillation Models. Surveys in Geophysics, 36(2), 295-324.

Secan JA, Bussey RM, Fremouw EJ, Basu S (1995) An improved model of equatorial scintillation. Radio Sci 30:607-617 Strang G., Introduction to Linear Algebra. Wellesley, MA: Wellesley Cambridge Press, 1998.

Wernik AW, Alfonsi L, Materassi M (2007) Scintillation modeling using in situ data. Radio Sci 42. doi:10.1029/2006RS003512

Wernik, A. W. and C. H. Liu (1974), Ionospheric irregularities causing scintillations of GHz frequency radio signalsJournal of Atmospheric and Solar- Terrestrial Physics, 36:871-879.

Yang L., Jihye Park, M Aquino, A Dodson, S.Veettil,

(2014). CALIBRA D3.1 Algorithms review.

(http: //www. calibra- ionosphere .net/calibra/documents/calibrad31.pdf) .

In the foregoing, preferred embodiments have been described and variations to the invention have been proposed, but it is to be understood that those skilled in the art will be able to modify and make changes to the invention without thereby falling outside the relevant scope of protection, as defined by the enclosed claims.

Claims

1. Method for the forecasting of at least a ionosphere scintillation parameter and/or a ionospheric total electron content using data from a GNSS network comprising scintillation and/or total electron content monitor receivers covering an interest area, the method comprising the execution of the following steps:

AO. acquiring real-time GNSS data from said monitor receivers along a predetermined time period, the data corresponding to a set of satellites in the fields of view of said monitor receivers;

Al . providing the ionospheric total electron content and/or said at least a ionosphere parameter values at Ionospheric Pierce Points based on said GNSS data, along said predetermined time period over a ionospheric domain Δ defined by the section of said fields of views with the ionosphere;

the method being characterized in that the following steps are executed:

B. subdividing said ionospheric domain Δ into domain portions Δκ;

C. solving for each domain portion Δ_k the transport equation:

Wherein dV is the measure of the domain portion, v_k is the velocity field of a scalar field /* in the domain portion Δ_k, which represents either total electron content parameter or a scintillation parameter at time t, said transport equation being solved for v_k and n_k which are assumed constant in Δ_κ, the value of the scalar field /' being taken as constant over the domain portion Δ* and determined on the basis of the ionospheric total electron content and/or the at least a ionosphere parameter values of step Al, wherein n_k - 0 if _ft is the total electron content parameter, the time derivative in the left-hand side of the transport equation being calculated on the basis of the values of f^l in a portion of said time period;

D; solving again for each portion Δκ the transport equation of step C to obtain the scalar field f^T at time T by using the velocity field t?_fc calculated in step C, wherein T is a time point after said time period.

2. Method according to claim 1, characterised in that said domain Δ is an area defined by the section of said fields of views with the ionosphere at an altitude value comprised between 350 and 400 km over the earth's surface.

3. Method according to claim 2, characterised in that said altitude value is around 350 km.

4. Method according to any claim 1 to 3, characterised in that each domain portion A_k is a triangle according to the Delaunay triangulation applied to said Ionospheric Pierce Points.

5. Method according to claim 4, characterised in that in step AO the time period is subdivided into time steps At, 2At, ... nAt, NAt, wherein n is the counter of the time step ranging from 1 to positive integer N, and in step C the transport equation for each domain portion gives rise to the following compact matrix equation system:

wherein M is a block matrix:

in which n-tb. block is:

Wherein and are the triangles adjacent to the

triangle A the term γ is the area of the

triangle Δ

_k in the case of said at least a scintillation parameter or 0 in the case of ionospheric total electron content; is the value of scalar field / at

time step nAt and ki is the index indicating the i-th side of the triangle; is the vector perpendicular

to i-th side of

triangle and assumed as positive when going out from such triangle; wherein vector s takes the following form:

And wherein:

With

Wherein V/,· is the velocity of the scalar field in the triangle adjacent to Ay on the i-th side of A*.

6. Method according to claim 5, characterised in that said block matrix is regularised according to Singular Value Decomposition.

7) Computer program, comprising code means configured in such a way that, when run on a computer, they perform the method of one of previous claims.

8. Method of GNSS precise positioning, using a GNSS network comprising scintillation and TEC monitor receivers covering an interest area, the method being characterised in that the following steps are executed: MP1. obtaining GNSS data from said receivers;

MP2. forecasting at least a scintillation parameter and/or total electron content values using the method of claim 1;

MP3. Providing GNSS precise positioning by either:

MP3_A. screening or weighting the GNSS data based on the output of step MP2, and calculating, by a receiver, the positioning based on the result of the screening or weighting;

MP3_B. calculating, by the receiver, the positioning based on the forecasted parameters of step MP2.

9) System (100) of GNSS precise positioning of a rover (20), comprising:

- a GNSS network comprising scintillation and total electron content monitor receivers (10) covering an interest area (60),

- a central unit (30) , comprising a computer logic having the program of claim 7 installed on it; - a communication network (40) between said central unit (30) and said receivers (10);

- a communication network (50) between said central unit (30) and said rover (20);

the central unit being provided with means configured to execute step MP3 of the method of claim 8.

Barzano & Zanardo Roma S.p.A.