IONOSPHERIC SCINTILLATION PREDICTION

Abstract

A method of generating scintillation prediction data comprises: •a) for a plurality of satellites (12A), and reference stations (10A, 10B, 10C, 10D) measuring phase scintillation data, satellite by satellite for each reference station during multiple epochs (t−2, t−1, t, t+k); •b) forecasting an expected phase scintillation value for each satellite and reference station for a period of at least 24 hours based on a cyclical prediction model; •c) for a given user location (20) and a given satellite, spatially interpolating the expected phase scintillation values of the plurality of reference stations to determine a predicted phase scintillation index; •d) repeating step c) for further satellites visible from the user location.

Claims

1. A method for generating ionospheric scintillation prediction data comprising: a) for a plurality of satellites, and reference stations, measuring phase scintillation data, satellite by satellite for each reference station during multiple epochs; b) forecasting an expected phase scintillation value for each satellite and reference station for a period of at least 24 hours based on a cyclical prediction model; c) for a given user location and a given satellite, spatially interpolating the expected phase scintillation values of the plurality of reference stations to determine a predicted phase scintillation index; and d) repeating step c) for satellites visible from the user location.

2. The method according to claim 1, wherein the method further comprises a step e) of converting the phase scintillation index for each satellite into a range error with respect to the user location.

3. The method according to claim 2, wherein the ionospheric scintillation prediction data comprises a Scintillation Dilution Of Precision (SDOP) index and step e) comprises computing the SDOP index by mapping the predicted phase scintillation indices of the satellites visible from the user location onto a position domain.

4. The method according to claim 1, wherein step b) takes place using a Holt-Winters prediction model.

5. The method according to claim 4, wherein step b) takes place using a multiplicative Holt-Winters prediction model, whereby for an expected phase scintillation value k epochs ahead, the sum of the local mean level and the trend over k epochs is multiplied by a factor indicative of the variation over the period.

6. The method according to claim 1, wherein step c) takes place by determining ionosphere pierce points (IPPs) for signals received from the given satellite for each of the reference stations

7. The method according to claim 6, further comprising performing linear spatial interpolation between the IPPs.

8. The method according to claim 6, further comprising determining the height of the IPPs by determining the time difference between scintillation start-up and sunset at ground level.

9. The method according to claim 1, further comprising providing the ionospheric scintillation prediction data to a user at the user location.

10. The method according to claim 1, wherein the ionospheric scintillation prediction data comprises loss-of-lock indicators and the method further comprises: e) for the plurality of satellites and reference stations, determining loss-of-lock occurrences for a carrier tracking loop of a Global Navigation Satellite System (GNSS) receiver at the reference station; f) forecasting expected loss-of-lock occurrences for each satellite and reference station for a period of at least 24 hours based on a cyclical prediction model; g) for the user location and given satellite, spatially interpolating the expected loss-of-lock occurrences of the plurality of reference stations to determine a loss-of-lock indicator for the given satellite at the user location.

11. The method according to claim 10, further comprising displaying the predicted loss-of-lock indicator to a user at the user location.

12. The method according to claim 1, further comprising determining a user location for a given user and generating the ionospheric scintillation prediction data on the basis of the determination.

13. The method according to claim 1, wherein the reference stations comprise multi-frequency receivers and step a) takes place by a linear combination of the dual frequency phase measurements.

14. The method according to claim 1, wherein the plurality of satellites are earth-orbiting satellites and the method further comprises estimating the phase scintillation information for ionosphere pierce points (IPPs) between the user location and geostationary satellites visible from the user location and transmitting the phase scintillation information to a user at the user location.

15. A system for generating ionospheric scintillation prediction data from Global Navigation Satellite System (GNSS) satellite signals and for providing the ionospheric scintillation prediction data to a user at a user location, the system comprising: a plurality of reference stations having multi-frequency receivers adapted to receive the GNSS satellite signals; a processing unit; and a memory storing instructions, which when executed by the processing unit causes the processing unit to: a) for a plurality of satellites, and the reference stations, measure phase scintillation data, satellite by satellite for each reference station during multiple epochs; b) forecast an expected phase scintillation value for each satellite and reference station for a period of at least 24 hours based on a cyclical prediction model; c) for a given user location and a given satellite, spatially interpolate the expected phase scintillation values of the plurality of reference stations to determine a predicted phase scintillation index; and d) repeating step c) for satellites visible from the user location.

16.-17. (canceled)

18. A non-transitory computer readable medium comprising instructions which when executed by a processing unit, causes the processing unit to: a. for a plurality of satellites, and the reference stations, measure phase scintillation data, satellite by satellite for each reference station during multiple epochs; b. forecast an expected phase scintillation value for each satellite and reference station for a period of at least 24 hours based on a cyclical prediction model; c. for a given user location and a given satellite, spatially interpolate the expected phase scintillation values of the plurality of reference stations to determine a predicted phase scintillation index; and d. repeating step c) for satellites visible from the user location.

19. The non-transitory computer readable medium of claim 18, further comprising instructions, which when executed by the processing unit, causes the processing unit to map ionospheric scintillation prediction data into a position domain on a user display.

20. The non-transitory computer readable medium of claim 18, further comprising instructions, which when executed by the processing unit, causes the processing unit to: enabling a user at a user location to receive ionospheric scintillation prediction data; and view the ionospheric scintillation prediction data on a user display at the user location.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0031] The features and advantages of the invention will be appreciated upon reference to the following drawings of a number of exemplary embodiments, in which:

[0032] FIG. 1 shows the geometry of the thin-layer ionosphere;

[0033] FIG. 2 shows a plot of daily repetition of scintillation along GPS PRN 12 trajectory from 24 Nov. to 4 Dec. 2013 in Recife;

[0034] FIG. 3 shows measurement of the time series of σ.sub.φ60 over a plurality of epochs for a given satellite and a given reference station;

[0035] FIG. 4 shows spatial interpolation for each ionospheric patch between a satellite and the reference stations;

[0036] FIG. 5 shows a given user location with the predicted phase scintillation indices of the visible satellites at the user location mapped onto the position domain;

[0037] FIG. 6 shows a plot of observed and predicted Loss-Of-Lock for 27th November 2013;

[0038] FIG. 7 shows an illustration indicating how the height of irregularities can be computed based on the geometry of the earth and the start time of scintillation;

[0039] FIG. 8 provides a plot of L-band scintillation measurement at IPPs for reference stations in South America and predictions for a user location;

[0040] FIG. 9 shows a flow chart of processing steps for the scintillation forecast of the invention; and

[0041] FIG. 10 shows time series of measured and predicted SD OP for a user location from 30 Nov. to 3 Dec. 2013.

DESCRIPTION OF ILLUSTRATIVE EMBODIMENTS

[0042] In order to better understand the context of the invention, a brief introduction to the calculation of ionospheric effects is provided, on which the invention is based. Due to the Earth's rotation, the ionosphere is a highly dynamic medium and the electron density can vary significantly from minute to minute at a given location resulting in temporal and spatial variations of TEC. FIG. 1 shows the geometry of the thin-layer ionosphere I at height H above the Earth E, including a satellite 12 visible from a user location 20. The line of sight between the satellite 12 and user location 20 passes the ionosphere I at the ionospheric pierce point (IPP). The satellite zenith angle ξ at the IPP is thus distinct from the satellite zenith angle ξ′ measured at the user location 20. For simplicity and to exclude height dependency of TEC, a thin-layer approximation is used for the ionosphere and the slant TEC is mapped onto the vertical TEC using a simple cosine function:

VTEC(t,φ,λ,H)=m(ξ).Math.TEC(t,{right arrow over (ρ)}(t))   (1))

{right arrow over (ρ)}(t)={right arrow over (r)}.sup.s(t)−{right arrow over (r)}.sub.R   (2)

Where φ and λ are geographic latitude and longitude of the ionospheric pierce point (IPP), H is the height of the thin ionosphere layer, ξ is the satellite zenith angle at the IPP, {right arrow over (r)}.sup.s and {right arrow over (r)}.sub.R stand respectively for geocentric position vectors of satellite and receiver.

[0043] The total VTEC variation along a satellite trajectory can be shown to be as follows:

[00001]$\begin{array}{cc}\mathrm{dVTEC}=\left(\frac{\partial \mathrm{VTEC}}{\partial t}+\frac{\partial \mathrm{VTEC}}{\partial \varphi }.Math.\frac{\partial \varphi }{\partial t}+\frac{\partial \mathrm{VTEC}}{\partial \lambda }.Math.\frac{\partial \lambda }{\partial t}\right).Math.\mathrm{dt}& \left(3\right)\end{array}$

Where

[0044] [00002]$\frac{\partial \varphi }{\partial t}.Math..Math.\mathrm{and}.Math..Math.\frac{\partial \lambda }{\partial t}$

stand for projected satellite velocity on the thin-ionosphere layer respectively along latitude and longitude.

[0045] For time series analysis, the pattern of VTEC observations can be decomposed into regular (i.e. trend: long-term and seasonal variations due to Earth daily rotation and annual revolution) and irregular (short term and random-like) variations. The regular variation may be a few TECU and the irregular variation may be as little as a tenth of a TECU.

[0046] GNSS signal scintillation is due to the irregular and rapid fluctuations of TEC in space. In order to measure scintillation, it is required that TEC is detrended. A simple way of detrending is between-epoch differencing of TEC observations. Assuming a 1 Hz or higher sampling rate (dt≦1 s) and a zero time derivative of VTEC yields:

[00003]$\begin{array}{cc}\mathrm{dVTEC}\approx \frac{\partial \mathrm{VTEC}}{\partial \varphi }.Math.d.Math..Math.\varphi +\frac{\partial \mathrm{VTEC}}{\partial \lambda }.Math.d.Math..Math.\lambda =\stackrel{.fwdarw.}{\nabla }.Math.\mathrm{VTEC}.Math..Math.\mathrm{with}.Math..Math.\stackrel{.fwdarw.}{\nabla }=\left[d.Math..Math.\varphi ,d.Math..Math.\lambda \right]& \left(4\right)\end{array}$

where {right arrow over (∇)} TEC denotes the gradient of VTEC. Between-epoch differencing of VTEC observations (with sampling rate≦1 Hz) leads to a detrended VTEC time series (dVTEC), which is more likely to consist of only the irregular and random-like fluctuations of VTEC along a satellite trajectory. dVTEC is the rate of TEC which is also referred to as ROTI.

Phase Scintillation and Loss-of-Lock

[0047] Any changes in TEC along a signal path result in a phase change of the signal. The following equation is used to convert dVTEC to the GNSS phase variation:

[00004]$\begin{array}{cc}d.Math..Math.\phi \left(t\right)=\frac{40.3}{c.Math.f}.Math.\mathrm{dVTEC}\left(t\right)& \left(5\right)\end{array}$

where f is the relevant GNSS signal frequency and c is the speed of light. Rapid fluctuation of TEC leads to phase scintillation of the signal. In case of severe phase scintillation, the carrier tracking loop in a GNSS receiver may have problems to track. For a receiver tracking an L-band signal, a rapid change of only 1 radian of phase in a time interval equal to the inverse of the receiver bandwidth can be enough to cause problems for the receiver's tracking loop. In this case, the signal amplitude also fades enough to result in Loss-Of-Lock and consequently a phase discontinuity or cycle slip. The fading of a signal due to ionosphere disturbances is called amplitude scintillation. In the following, one radian is used as a threshold to detect cycle slips in the time series of dφ(t). Note that one radian corresponds to ˜0.2 TECU (and ˜3 cm range error for the GPS L1 frequency). The number of cycle slips per minute is used as an index (denoted by NLL.sub.60) to measure the severity of the amplitude scintillation along a satellite trajectory. In case of no amplitude scintillation, NLL.sub.60 will be equal to zero; in the most extreme case NLL.sub.60 can never exceed 60.

[0048] To measure phase scintillation, the standard deviation of dφ over every 60 seconds is calculated (denoted by σ.sub.φ60) excluding the cycle slips. Note that the phase scintillation index will be 0<σ.sub.φ60<1 radians.

Predictability of the Scintillation

[0049] The electron density irregularities in the equatorial region take place between sunset and midnight with activity occasionally continuing until dawn. The irregularities are more likely to repeat daily after sunset at a given location in the equatorial region.

[0050] In order to show daily repetition of the irregularities in space, the phase scintillation index σ.sub.φ60, determined at Recife for PRN 12 from 24 Nov. to 4 Dec. 2013, is plotted in FIG. 2, whereby the solid circles show the IPPs with severe phase scintillation (0.7≦σ.sub.φ60 radians). It can clearly be seen in FIG. 2 that scintillation varies for different ionospheric patches.

[0051] Time series measurements of σ.sub.φ60 and NLL.sub.60 determined at Recife for PRN 12 from 24 Nov. to 4 Dec. 2013 have shown that repetition of scintillation in time takes place from one day to the next.

Prediction of Phase Scintillation

[0052] The repeatability of the scintillation in time and in space provides an opportunity to statistically predict the scintillation and provide a real-time scintillation forecast system.

[0053] For predicting scintillation effects on positioning at a given user location, the following steps are performed: [0054] 1) Temporal prediction of phase scintillation index satellite by satellite at the reference stations, [0055] 2) Spatial interpolation of the predicted phase scintillation index between reference stations for the user location, [0056] 3) Mapping phase scintillations of visible satellites at the user location into the position domain.

[0057] 1) Temporal Prediction of Phase Scintillation

[0058] As discussed above and on the basis of measurements it has been found that there are scintillation periods after sunset which almost repeat every day due to the earth's rotation. Assuming a day to be a full period, the time series of σ.sub.φ60 contains a repeated trend and a seasonal pattern. There are two seasons in the course of the day: scintillation and non-scintillation hours.

[0059] A Holt-Winters method is used for forecasting phase scintillation (σ.sub.φ60) for the next 24 hours. The Holt-Winters method is a statistical forecast method that captures appropriately the seasonality of data. There are two Holt-Winters models: the Additive model and the Multiplicative model. The main distinction between additive and multiplicative models is seasonality. The additive model is appropriate if the seasonal pattern is almost constant through the time series. If the seasonal pattern is proportional to the level (or local mean) of data, the multiplicative model is preferred.

[0060] Data analysis of the time series of σ.sub.φ60 for a given location has shown that the time series of σ.sub.φ60 contains a constant trend but the level of data is not constant. It can be considered that the seasonal pattern in the course of day is proportional to the level of data. Therefore, the multiplicative Holt-Winters model for forecasting the phase scintillation is used.

[0061] For a measured time series of σ.sub.φ60, denoted by X(t.sub.1), X(t.sub.2), . . . , X(t.sub.n), to forecast X(t.sub.n|k), the forecast made at time t.sub.n, k steps ahead, is denoted by {circumflex over (X)}.sub.n(t.sub.n+k). The mathematical equations of the multiplicative Holt-Winters model are written as:

L(t.sub.n)=α.Math.(X(t.sub.n)/I(t.sub.n−p))+(1−α).Math.(L(t.sub.n 1)+T(t.sub.n 1))   (6)

T(t.sub.n)=β.Math.(L(t.sub.n)+L(t.sub.n−1))+(1−β).Math.T(t.sub.n−1)   (7)

S(t.sub.n)=γ.Math.(X(t.sub.n)/L(t.sub.n))+(1−γ).Math.S(t.sub.n−p)   (8)

where L(t.sub.n), T(t.sub.n), and S(t.sub.n) stand respectively for the local mean level, trend, and seasonal factor at time t.sub.n and p is the number of observation epochs in the seasonal cycle. The smoothing parameters α, β, and γ are constants between 0 and 1 and used for updating the level, trend, and seasonal factor. The forecast made at time t.sub.n for k=1, 2, . . . , p is given by

{circumflex over (X)}.sub.n(t.sub.n+k)=(L(t.sub.n)+k.Math.T(t.sub.n)).Math.S(t.sub.n+k−p)   (9)

[0062] Note that the time interval of the time series of σ.sub.φ60 is 1 minute because the phase scintillation index σ.sub.φ60 is computed every minute (i.e. p=1440 for the seasonal cycle as one day contains 1440 minutes).

[0063] FIG. 3, shows a reference station 10A observing an orbiting satellite 12A during a number of epochs (t−2) to (t+k). The time series of σ.sub.φ60 is measured over a plurality of epochs for the satellite 12A, with the observations taking place through successive IPPs 14A(t−2) to 14A(t+k) defined within successive ionospheric patches 16(t−2) to 16(t+k). This is repeated for all further reference stations within the region.

[0064] 2) Spatial Interpolation of Phase Scintillation (Kriging Method)

[0065] In order to compute the phase scintillation index between a given user and a visible satellite, spatial interpolation is required. The spatial interpolation is carried out individually for each ionospheric patch (between a satellite and the reference stations), epoch by epoch, and is shown in FIG. 4 using the Kriging interpolation method. Kriging is a linear spatial interpolation method that takes the spatial correlation into account. It is also a best linear unbiased predictor if the covariance function is correctly selected. The covariance function is chosen as a simple linear function of the distance; i.e. c.sub.ij=L−l.sub.ij where c.sub.ij and l.sub.ij stand respectively for covariance and distance between IPPs of i and j, and L is an arbitrary constant that should be greater than the maximum distance between IPPs According to FIG. 4, a regional network of reference stations 10A to 10D receives signals from satellite 12A and provides the observations to a processing unit 22 at a Network Control Center (NCC) to process the data to compute phase scintillation indices for reference stations for the satellite at epoch (t). Then the processing unit 22 at time t can predict the phase scintillation indices for t+k epochs ahead for the reference stations for the given satellite. The processing unit 22 uses the Kriging interpolation method to calculate a prediction of phase scintillation index of IPP between a given user location and the satellite for t+k epochs ahead. The predicted scintillation index for the user can be either sent to the user by email or shown on the website.

[0066] 3) Scintillation Effect on GNSS Positioning

[0067] In order to measure the phase scintillation effect on GNSS positioning at a given user location the predicted phase scintillation indices of the visible satellites at the user location are mapped onto the position domain, as shown in FIG. 5. To do this, first the predicted phase scintillation index, between user u and satellite s, is converted to the range error dσ.sub.u.sup.s in meters, which for L1 is:

dρ(t).sub.u.sup.s=λ.Math.{circumflex over (σ)}.sub.φ60(t).sub.u.sup.s   (10)

where λ is the GNSS signal wavelength (0.19m for GPS L1). Then, the mathematical model of standard single point positioning is used to compute positioning error:

{right arrow over (r)}.sub.u(t)=(A.sup.TWA).sup.−1A.sup.TWY.sub.u and SDOP(t)=∥{right arrow over (r)}.sub.u(t)∥   (11)

Y.sub.u=[dρ(t).sub.u.sup.1, dρ(t).sub.u.sup.2, . . . , dρ(t).sub.u.sup.n].sup.T   (12)

where the design matrix A, weight matrix W and observation vector Y are time-dependent, n is number of visible satellites. In the weight matrix, the weight of a satellite is defined as a function of elevation angle of the satellite. We call the computed positioning error Scintillation Dilution Of Precision and denote it by SDOP. An ionosphere free solution e.g. for L1 and L2 would increase the SDOP index by a factor of 3.

[0068] The SDOP is an index that can be used as a measure of scintillation effect on precise positioning at the user location. It depends on the geometry of the visible satellites and phase scintillation status at the satellites' ionospheric pierce points. When there is no scintillation, dρ(t).sub.u.sup.s values are around 3 mm, resulting in SDOP to be at the level of a few mm. In case of scintillation, dρ(t).sub.u.sup.s values corresponding to the satellites affected by scintillation are increased, resulting in SDOP values that can be as large as several cm. An increase in SDOP leads to low accuracy of the ambiguity float solution and, consequently, to problems for ambiguity resolution in precise positioning models.

Prediction of Loss-of-Lock

[0069] As mentioned, deep signal fading caused by scintillation can lead to loss-of-lock of the carrier tracking loop in a GNSS receiver. It should be noticed that loss-of-lock depends also on quality and setup of receiver and antenna. For example, if a long cable is used for connection between antenna and receiver, it is more likely to experience loss-of-lock than the case of a short-cable. Because of this, it is difficult, if not impossible, to fully predict loss-of-lock at a given user location. Therefore, we only specify satellites which will have a chance for loss-of-lock in the next 24 hours.

[0070] In commercial GNSS reference networks, high quality receiver-antenna setups are generally used. If a reference station experiences loss-of-lock with respect to a satellite, then a user in the same area of the reference station will most likely also experience loss-of-lock.

[0071] NLL.sub.60 (the number of cycle slips per minute) is used as a measure of chance for loss-of-lock. In order to predict NLL.sub.60 for the next 24 hours, a loss-of-lock probability for the next 24 hours is assumed to be exactly the same as for the last 24 hours i.e. a simple unitary forecast model. Therefore, the measured time series of NLL.sub.60(t).sub.r.sup.s at reference station r and satellite s is copied from the last 24 h to the next 24 h, e.g.

N{circumflex over (L)}L.sub.60(t+24 h).sub.r.sup.s=NLL.sub.60(t).sub.r.sup.s.   (13)

[0072] Then, the Kriging method is used for interpolating the predicted N{circumflex over (L)}L.sub.60 values epoch by epoch for the given user location.

[0073] If the summation of interpolated N{circumflex over (L)}L.sub.60 over an hour is not zero then there is a probability of loss-of-lock during this hour. In FIG. 6, as an example, the chance of loss-of-lock at the user location is visualized for all satellites for 27th November 2013. In this figure, the past predictions are shown on the left of the figure with the future predictions shown on the right. The numbers indicate the actual number of observed cycle-slips during an hour as measured at the reference stations. This assists in determining the accuracy of past loss-of-lock prediction.

Adaptation Height of Scintillation

[0074] In the equatorial region, the irregularities of electron density (or plasma bubbles) mostly are generated after sunset at height H of the F-layer of ionosphere above the equator E. As shown in FIG. 7, for the sun at S, sunset is at point A on the ground and point B at the F-layer. The height of irregularities can be computed based on the start time of scintillation as follows:

[00005]$\begin{array}{cc}H=\frac{1-\cos \left(a\right)}{\cos \left(a\right)}.Math.R& \left(14\right)\end{array}$

where H is the height of the ionosphere and R the radius of earth. The difference between longitudes of A and B is denoted by a which corresponds to the time distance between the sunset on the ground and the start of scintillation.

[0075] In the evaluation of the IPPs the computed height H is used, being the height of the ionosphere thin-layer. It helps to adapt daily the ionospheric pierce points based on the height of the irregularities. This is then used in the subsequent Kriging and for the plotting of the scintillation map.

L-Band Scintillation

[0076] The geostationary satellites, which are at an almost fixed height and longitude above the equator, are used to provide GNSS corrections to users. These satellites use L-band signals for their communication. In case of severe scintillation, users may lose lock to these correction signals.

[0077] L-band scintillation values can be provided for all visible geostationary satellites at the user location. The time series of scintillation index is computed for the IPPs between the user and the L-band satellites. It will be understood that since the geostationary satellites do not generally transmit data sufficient to calculate ionospheric scintillation indices directly, interpolation within a patch to determine L-band scintillation data based on indices generated from the earth-orbiting satellites is a convenient way of generating valuable data. FIG. 8 provides a plot of L-band scintillation measurement at IPPs for reference stations in South America and predictions for a user location.

Flow Chart of Scintillation Forecast

[0078] To better illustrate the process according to an embodiment of the invention a flow chart of processing steps for forecasting scintillation is shown in FIG. 9. In the illustrated example shown in FIG. 9, the steps are carried out at the Network Control Center (NCC) for a dual frequency GNSS system. It will nevertheless be understood that certain steps may alternatively be carried out locally at the user location. In a first step 101, data from the reference GNSS stations is transmitted to the NCC, where it is received by the processing unit 22 for subsequent processing. The data includes the code and phase observables and Signal-to-noise ratio (SNR) at the two frequencies. In step 102, TEC values for the respective IPPs are calculated for each satellite and each reference station using the phase geometry-free linear combination between the two frequencies, with the slant TEC being converted to VTEC by equations (1) and (2) above.

[0079] In step 104, the TEC data is detrended by between-epoch differencing to compute the rate of TEC according to equation (4) and then the scintillation data (dφ) computed using equation (5). The scintillation data is then subject to various filtering steps at 106, primarily to remove station-dependent noise from the data. The station dependent noise may include local interference and other effects.

[0080] At step 108, the phase scintillation index is generated by calculating, the standard deviation of the change of phase dφ over every 60 seconds (denoted by σ.sub.φ60) according to equation (5). This data is then used for forecasting the expected phase scintillation index for the period 24 hours ahead using the Holt-Winters method according to equation (9) in step 110. This is performed satellite-by-satellite for each of the satellites observed by a given reference station and for each reference station.

[0081] In parallel to steps 108 and 110, loss of lock data, NLL.sub.60 (the number of cycle slips per minute) is generated in step 112 for each of the reference stations and satellites. Using equation (13) a predicted loss of lock for the next 24 hours is defined.

[0082] The data from both steps 110 and 112 is spatially interpolated in step 114 using Kriging to provide data specific for a given user location, based on the respective IPP for that user location and observations of a given satellite.

[0083] In step 116, the range error computed according to equation (10) is used to determine the SDOP for the user location according to equation (11). This can be plotted as a time series for a given period e.g. past and future 24 hour period on a user display (not shown). Additionally, at step 120, a regional map can be displayed for relevant data, including the phase scintillation index, the L-band scintillation values at the IPPs of the geostationary satellites and the loss-of lock probability.

Results of the Scintillation Forecast

[0084] In order to show performance of the Holt-Winters method, six reference stations in South America were used to predict scintillation of a user location: Lat. 11° S and Lon. 36° W. Data was processed from 27 to 30 Nov. 2013 when scintillation was severe.

[0085] The time series of σ.sub.φ60 between reference station Recife and five GPS satellites from 27 to 29 November was measured and the predicted σ.sub.φ60 calculated. The selected five satellites were affected by scintillation for a few hours. The day to day variability of the scintillation level and duration was observable in the data.

[0086] Using the predicted σ.sub.φ60 at the reference stations, the SDOP of the user was computed for 4 days from 30 Nov. to 3 Dec. 2013. The time series of the measured and predicted SDOP (at UTC 16:00) are shown in FIG. 10. The black curve over the last 24 hours is the measured SDOP which becomes the predicted SDOP for the next 24 hours. To control the quality of the prediction, the previously predicted SDOP for the last 24 hours is shown in light shading. As seen the timing of the predicted SDOP for the severe scintillation hours fit well to the measured SDOP. The level of predicted SDOP could be improved by optimization of the smoothing parameters in the Holt-Winters method based on further analysis of the observations over the previous 24 hours.

[0087] The repeatability of the small-scale irregularities in time and space allows scintillation to be predicted using the disclosed near real-time scintillation forecast system. Testing of the scintillation forecast for low-latitude regions (e.g. West-Africa, South-America and India) has shown promising results that are subject to further optimization. Statistical analysis of the prediction quality resulted in 70% correct prediction of severe scintillation.

[0088] Because of the seasonality pattern in the scintillation data, the Holt-Winters method is an effective method for predicting the phase scintillation in time. The Kriging interpolator is an appropriate method for spatial prediction of the phase scintillation. Combination of the Holt-Winters method and Kriging provides a tool for statistically predicting the scintillation for the next 24 h for any given location inside or outside the reference network. It will equally be understood that the same process may be used for predicting over periods of multiple days, should this be required.

[0089] Thus, the invention has been described by reference to the embodiments discussed above. It will be recognized that these embodiments are susceptible to various modifications and alternative forms well known to those of skill in the art without departing from the spirit and scope of the invention. Accordingly, although specific embodiments have been described, these are examples only and are not limiting upon the scope of the invention.