Technical field

[0001] The present invention is related to methods and systems for detecting disturbances in signals received from Global Navigation Satellite System (GNSS) satellites. In particular, the present invention is related to methods and systems for detecting ionospheric scintillations and other line-of-sight disturbances in such signals.

Background art

[0002] Global Navigation Satellite System (GNSS) is a general name for a satellite navigation system. A GNSS system comprises, among other things, a constellation of multiple navigation satellites, which are orbiting the earth and which transmit dedicated navigation signals also called ‘ranging signals’. Most known is the Global Positioning System (GPS) developed by the United States Government. Other systems include the Galileo system which is currently in development by the European Space Agency (ESA), the Russian GLONASS and the Chinese BeiDou system. The primary purpose of these systems is to provide position, velocity and time to any user on or near the Earth’s surface. The user determines his position in space and time by measuring his pseudorange to at least 4 satellites whose position and time are accurately determined, the range being defined as the distance between the user and a satellite.

[0003] Proper operation of high-rate GNSS receivers depends on high- quality phase and code observations. Especially the quality of the phase observations is important for the high-accuracy positioning modes RTK (Real-Time Kinematic) and PPP (Precise Point Positioning). Early detection of reduced quality measurements due to external circumstances is important to warn the user of the GNSS receiver that normal operation is blocked or limited. Two such circumstances are ionospheric scintillation and operation under foliage. In both cases the signal propagations along the lines-of-sight to the satellites are compromised. Because this may hold for a subset of all lines-of-sights, limited operation may still be possible if the compromised signals are discarded in the RTK or PPP positioning algorithms and only the signals of the remaining satellites are being used.

[0004] Each GNSS satellite sends signals at different carrier frequencies.

These carrier frequencies differ per constellation and in fact even per satellite block within a constellation. GPS for example uses two frequencies (L1 = 1575.42 MHz and L2 = 1227.60 MHz) in the older blocks II and three frequencies in the newer blocks I l-F and III (L1 , L2, and L5 = 1176.45 MHz). On these frequencies different signals may be encoded, e.g. L1CA, and L1 P(Y) on L1 , and L2P(Y) on L2 on older satellites and L1C and L2C on newer satellites. These signals can be tracked by a GNSS receiver resulting in code measurements, but also the carrier phases of the signals can be tracked resulting in phase measurements which are much more accurate than the code measurements. Despite an additional phase cycle bias, the phase measurements are therefore much more suitable for the precise positioning techniques. Because on older GPS satellites there is no L2C signal and the L2P(Y) code is encrypted, the L2 frequency is only accessible through a semi-codeless tracking technique. With this technique the resulting L2P phase measurement is correlated with the L1 phase measurement.

[0005] Measurement of scintillation is possible with dedicated GNSS receivers (ionospheric scintillation monitoring receivers, or ISM Rs) that output standardized scintillation indices based on the sampling of the signal amplitude and phase at a high rate (typically 50 Hz or more). Besides the need for operating at high rate, these receivers require a very stable clock reference to distinguish phase variations due to the ionosphere from phase variations due to the clock jitter. This comes with a significant cost, size and power consumption penalty. Standard GNSS receivers typically operate at lower rates, such as about 1 Hz, and use cheaper and smaller clocks.

[0006] US 2015/0253431 describes a method for detecting ionospheric scintillation by calculating a geometry free combination (GFC) parameter based on code phase pseudo-range measurements of two signals received from a same GNSS satellite having different frequencies. The GFC parameter is calculated by:

GFC = ( Φ ₁ — Φ ₂ ) / (μ*-1), wherein Φ ₁ and Φ ₂ are the phase measurements of the two signals, and μ* = f ₁ ² / f ₂ ² wherein fi and f2 represent the carrier frequencies of the two signals. The GFC parameter provides an estimation of the ionosphere signal delay for the signals. A dispersion of the GFC parameter, which can be represented by the variance of the GFC, is used for detecting the presence of scintillation. If a value of the dispersion exceeds a predetermined threshold value, ionosphere scintillation is detected.

[0007] One drawback of the above method is that it is slow, as it needs to gather a batch of data over a certain period of time before being able to determine the variance. This is because detrending needs to be applied to eliminate unknowns in the equations and compute a variance. This may pose problems in high-accuracy positioning modes of GNSS receivers, where early detection of compromised signals is critical for reliable operation. Another drawback is that computing a variance over a sample time can result in inaccurate estimates of the variance, in particular when the sample consists of both a low noise part and a high noise part, which would occur exactly when scintillation starts kicking in.

[0008] WO 2009/125011 describes a method for vector phase tracking a plurality of global positioning satellite carrier signals. In the method, the ionospheric drift is estimated by a linear minimum mean square error estimator using the square difference of carrier phase measurement and ionospheric drift for multiple epochs. Such an estimator suffers from the same drawbacks as indicated above, i.e. it is slow in that it needs to gather a batch of data over multiple epochs in order to provide an estimate.

Summary

[0009] There is therefore a need in the art for methods allowing to estimate satellite signal disturbances more reliably, more accurately and/or faster than prior art methods and without requiring the expensive infrastructure of current ISM Rs. There is also a need in the art for GNSS receivers, in particular non-ISMR GNSS receivers implementing such improved methods.

[0010] According to a first aspect of the invention, there is therefore provided a method of detecting a disturbance in a (carrier) signal of a GNSS satellite as set out in the appended claims.

[0011] A method as described herein comprises receiving from the GNSS satellite a first signal and a second signal, the second signal having a carrier frequency different from a carrier frequency of the first signal and measuring a carrier phase difference between the first signal and the second signal. According to an aspect of the invention, the method comprises determining an ionosphere delay rate. The ionosphere delay rate can be determined for a reference frequency, which can be the carrier frequency of the first signal or the second signal, or any other suitable frequency of a signal emitted by the GNSS satellite. The method further comprises determining a first indicator based on a quadratic expression comprising or consisting of a squared sum of a first term and a second term. The first term comprises or consists of the measured carrier phase difference and the second term comprises the determined ionosphere delay rate.

[0012] According to the present invention, a quadratic expression (i.e. an expression comprising a squared term) is used as basis to construct an indicator for detecting the disturbance. The indicator is advantageously a variance tracer, in particular a tracer for the thermal noise variance. The quadratic expression is advantageously an observable for the indicator. By using a quadratic expression comprising the squared sum of both the measured carrier phase difference between two or more signals and the ionosphere delay rate, a more accurate indicator/tracer is obtained, which advantageously provides a direct and unbiased indication of the thermal noise variance, hence allowing for more accurate detection. The quadratic expression is advantageously defined by ½ (ΔΦ ₁₂ + μİ _r Δt) ² , in which ΔΦ ₁₂ represents the measured carrier phase difference and i _r represents the determined ionosphere delay rate. Furthermore, the indicator so obtained provides faster convergence and eliminates the need to take batches of data taken over multiple epochs.

[0013] Advantageously, the first indicator is determined by a first Kalman filter, i.e. the first indicator is an estimated state obtained by the first Kalman filter. Advantageously, the quadratic expression is used in the first Kalman filter as an observable. By using a Kalman filter to determine the indicator or tracer, an even faster initialisation is obtained. Furthermore, the Kalman filter only uses input measurements of the current and previous epoch, and therefore eliminates any waiting times due to the use of batches of data over multiple epochs. Yet additionally, using a quadratic expression as observable in the Kalman filter allows to obtain a direct and unbiased estimation of the thermal noise variance, providing more accurate results compared to the indicators of the prior art, such as in US 2015/0253431 which use an estimate of the ionosphere delay rate variance instead, which in turn may be an indicator of the thermal noise variance.

[0014] Advantageously, the measured carrier phase difference is used as an input for determining the ionosphere delay rate. Advantageously, the ionosphere delay rate is determined by a second Kalman filter. By using a Kalman filter to determine the ionosphere delay rate, a more accurate result for the ionosphere delay rate can be obtained, which can furthermore be obtained faster. Advantageously, the second Kalman filter uses as input the measured carrier phase difference. Even more advantageously the output of the first Kalman filter is used to determine a variance in the second Kalman filter, for even more accurate results.

[0015] According to a second aspect of the invention, there is provided a method of determining navigation data based on GNSS navigation signals as set out in the appended claims. Methods according to the second aspect comprise checking a signal disturbance by performing the method according to the first aspect. For the signal in respect of which a disturbance is detected, the method comprises excluding a measurement of the signal, or applying a reduced weight to the measurement when determining the navigation data, wherein the measurement can refer to a carrier phase measurement and/or a code measurement.

[0016] According to a third aspect of the invention, there is provided a

GNSS receiver system, or an ionospheric scintillation monitoring device, as set out in the appended claims. The GNSS receiver system comprises a computing unit configured or programmed to carry out the method of the first aspect and possibly the method of the second aspect. The ionospheric scintillation monitoring device comprises a computing unit configured or programmed to carry out the method of the first aspect above.

[0017] Described herein are methods that detect carrier phase fluctuations on a satell ite-by-satellite basis. The methods are light on CPU and can run at low observation intervals, e.g. 1 second. As a result, methods as described herein can advantageously be used in standard GNSS (positioning I navigation) receivers, providing accurate results not only at higher sampling rates of 50 Hz, but even at lower sampling rates, e.g. 1 Hz or 10 Hz. The methods described herein generally do not make a distinction between scintillation and foliage as cause for the phase fluctuations. Methods as described herein are targeted to static or slow-moving rover GNSS receivers, but can also run on measurements of an RTK base receiver. In the latter case, the measurements from the RTK base can be used to steer operation of the rover positioning method.

Brief description of the figures

[0018] Aspects of the invention will now be described in more detail with reference to the appended drawings, wherein same reference numerals illustrate the same features and wherein:

[0019] Figure 1A represents a block diagram of a GNSS receiver system as described herein. Figure 1 B represents an operational diagram of the GNSS receiver system of Fig. 1A operating according to the present invention;

[0020] Figure 2 represents a flow diagram of methods according to the present invention;

[0021] Figure 3 represents a diagram of two Kalman filters implementing the Thermal noise filter (TNF) and the Ionosphere delay rate filter (IDRF) as described herein;

[0022] Figure 4 represents TNF results of PRN16 for a Singapore 2013 dataset after IDRF feedback to the TNF was applied according to the scheme of Fig. 3; [0023] Figure 5 represents TNF results of PRN16 for the same Singapore

2013 dataset of Fig. 4 before IDRF feedback to the TNF was applied; [0024] Figure 6 represents IDRF results of PRN16 for a Singapore 2013 dataset after IDRF feedback to the TNF was applied;

[0025] Figure 7A represents TNF estimates and model values under occurrence of a scintillation event. The dots show the root of the observation values;

[0026] Figure 7B represents the scintillation factor, i.e. the ratio between estimated and modelled values. A clear scintillation peak can be seen between 17.30 h and 18.00 h;

[0027] Figure 7C represents corresponding estimated ionosphere delay rate and the observations used;

[0028] Figure 8A represents IDRF results of a L1-L2C combination recording in Singapore 2013;

[0029] Figure 8B represents IDRF results of a L1-L2P combination recording in Singapore 2013 at a same time corresponding to the recording of Fig. 8A.

Detailed Description

[0030] Referring to Fig. 1A, a GNSS receiver system 100 comprises a receiver unit 10, referred to as a ‘receiver’, and an antenna 14 for receiving GNSS navigation signals 21 from a plurality of navigation satellites 20, 20’, 20”. The navigation satellites can all belong to a specific GNSS constellation such as though not limited to GPS, Galileo, GLONASS, etc. but alternatively they can belong to more than one specific GNSS constellation. The GNSS receiver unit 10 comprises a positioning module 11 and a tracking module 12 coupled to antenna 14. The tracking module 12 is configured to perform code and phase measurements on the navigation signals 21 received from satellites 20, 20’, 20”. The positioning module 11 is configured to process the code and phase measurements from tracking module 12 to compute navigational data, such as position, velocity, elevation, etc. The positioning module 11 can be configured to operate according to one or a plurality of operational modes, such as PPP or RTK.

[0031] The antenna 14 can be mounted in a fixed relationship to the receiver unit 10. Alternatively, the antenna 14 is movable with respect to the receiver unit 10, and can be coupled to the receiver unit 10 through a wired or wireless connection.

[0032] Positioning module 11 is operable to process the code and phase measurements determined by tracking module 12 and determine a navigation output, such as position, velocity, time, etc. Positioning module 11 feeds the navigation output to an output module 13 of the GNSS receiver system, such as a user interface, or any other suitable interface communicating with an external device. [0033] The navigation signals 21 comprise carrier signals modulated by pseudo-random binary codes that are unique to each satellite 20, 20’, 20”. The tracking module 12 measures the time difference between the signal emission by the satellite and its reception by the antenna 14, resulting in code measurements. This time difference, multiplied by the speed of light, is the pseudorange. The term “pseudorange” is used instead of “range” because the time difference is measured with non-ideal clocks and contains a clock offset term. To measure the time difference, the receiver synchronizes a local replica of the spreading code with the spreading code incoming from the satellite. The time difference is the delay that the receiver had to apply to its local replica to synchronize it with the incoming code.

[0034] To increase positioning accuracy, the tracking module 12 is configured to additionally measure the carrier phase of the signals received from satellites 20, 20’, 20”. This allows to obtain a higher positioning accuracy. The precision of the phase observables is typically about 1 percent of the carrier wavelength, i.e. ~0.01 x 20 cm = 2 mm. The precision of the code measurements is typically about 0.1 percent of the chip length, i.e. ~0.001 x 300 m = 30 cm.

[0035] In the context of the present description, we assume that satellites

20, 20’, 20” each emit at least two signals having different carrier frequencies f ₁ and f ₂ , i.e. the navigation signals 21 of each satellite 20, 20’, 20” comprise at least two signals with carrier frequencies f ₁ and f ₂ respectively. The carrier frequencies fi, f2 refer to carriers in different frequency bands at which the satellite emits navigation signals, e.g. the L1, L2 and L5 frequency bands for GPS satellites.

[0036] Referring to Fig. 1 B, the tracking module 12 is configured to perform a phase measurement on the signals, resulting in measured phases Φ ₁ and Φ ₂ for the two different carrier frequencies f ₁ and f ₂ respectively. Tracking module 12 determines for each satellite one or more code measurements P ₁ and P ₂ . Additionally, tracking module 12 determines for each satellite carrier phase measurements Φ ₁ and Φ ₂ for signals with carrier frequency f ₁ and f ₂ respectively. This is shown diagrammatically in Fig. 1B by the different boxes 200, 201, 202 representing the measurements made by tracking module 12 for satellites 20, 20’ and 20” respectively.

[0037] According to an aspect of the present invention, the carrier phase measurements are utilized to detect a disturbance in the signal of a respective satellite. To this end, receiver unit 10 comprises a disturbance detection module 15, which analyzes the phase measurements Φ ₁ and Φ ₂ on a satellite-by-satellite basis to detect whether the signal is disturbed. If a disturbance is detected, disturbance detection module 15 is configured to provide an indication to the positioning module 11 that the navigation signal of a particular satellite is disturbed. Based on the indication, the positioning module 11 is configured to either discard the phase measurements and/or code measurements of the particular satellite for which a disturbance is detected, or to apply a lower weighting value to these measurements when determining the navigation output.

[0038] The disturbance detection module 15 is implemented with a disturbance detection method according to aspects of the invention as described below. In the following, the disturbance detection method will be described considering the L1-L2 carrier phase differences as inputs, where L1 and L2 are the GPS L1 and L2 carrier frequencies. It will however be appreciated that methods according to the present invention can operate on any two or more carrier frequencies for any GNSS constellation, such as GLONASS, Galileo, BeiDou, QZSS, etc.

[0039] Firstly, some theoretical considerations on a tracer for ionospheric scintillation are given. The L1-L2 carrier phase difference Φ ₁₂ (k) = Φ ₁ (k) - Φ ₂ (k) at epoch t(k) depends on the following parameters: Φ ₁₂ (k) = — I ₁₂ (k) + A ₁₂ + mp ₁₂ (k) + b ₁₂ (k) + φ ₁₂ (k) + ε ₁₂ (k) , (1) where I ₁₂ (k) = I ₁ (k) - / ₂ (k) is the ionosphere delay difference, A ₁₂ the constant ambiguity term, mp ₁₂ (k) the phase multipath difference, b ₁₂ (k) the hardware delay difference, φ ₁₂ (k) the phase wind-up difference, and ε ₁₂ (k) the thermal noise.

Because with λ ₁ , λ ₂ , λ _r , referring to the wavelength of respectively the carrier frequency f ₁ , the carrier frequency f ₂ and a reference frequency r, the ionosphere delay difference can be expressed as

The reference frequency r can be the frequency of the first carrier, the second carrier or any other frequency. The selection of the reference frequency depends on the GNSS system. Typically, the first frequency is used, i.e. . However, for GLONASS, which uses different frequencies for different satellites, the nominal first frequency, i.e. the first frequency of a satellite with frequency number zero is used as the reference frequency. If the hardware delay rates, phase wind-up rates and multipath rates are neglected, and it is assumed that there are no cycle slips, the L1-L2 carrier phase time difference yields: Since the ionosphere delay changes slowly, it can be assumed that between two epochs (typically one second) the ionosphere delay rate is constant, yielding: where Δt(k) = t(k) - t(k - 1) is the time difference between the two epochs, and İ _r (k) is the ionosphere delay rate. Dividing by Δt(k) gives: which Pi, X., A. J. Mannucci, U.J. Lindqwister, and C.M. Ho (1997), Monitoring of global ionospheric irregularities using the worldwide GPS network, Geophysical Research Letters, Vol. 24, No. 18, pp. 2283-2286 call the Rate Of TEC (total electron content), abbreviated as ROT. They define the standard deviation of the ROT over a certain time, the Rate Of TEC Index (ROTI), as an indicator of scintillation. In both Pi et al. and US 2015/0253431 described above, the variances are driven by the L1-L2 thermal noise ε ₁₂ .

[0040] Conversely, in methods of the present invention, a different approach is taken. The variance of the thermal noise is taken as a tracer for scintillation and other line-of-sight disturbances.

[0041] From Eq. (5), it has been observed that squaring the term

ΔΦ ₁₂ (k) + μ İ _r (k)Δt(k) gives a variance tracer. To this end, an observable is formed by

In Eq. (7) the factor ½ is optional. It is introduced here to obtain a direct observable of the thermal noise variance because measurements of two (consecutive) epochs will be used to determine the input noise, as will be described below. To determine the observable of Eq. (7), two factors need to be determined, viz. the carrier phase difference at epoch k: ΔΦ ₁₂ (k) and the ionosphere delay rate at epoch k: İ _r (k). From Eq. (3) it can be seen that μ represents the relation between the delay rates at different frequencies, making Eq. (7) independent of the choice for a reference frequency. In aspects of the present invention, a disturbance is detected on the basis of the observable of Eq. (7).

[0042] In the context of the present disclosure the term carrier phase difference ΔΦ ₁₂ refers to a difference of the carrier phase between two signals having different carrier frequencies f ₁ and f ₂ and between two (consecutive) epochs, e.g. k and k-1. Such carrier phase difference may be referred to as carrier phase double difference.

[0043] In the context of the present disclosure, the ionosphere delay rate can refer to a change in ionosphere delay along the line of sight. The ionosphere delay rate can have both a temporal and spatial component, and can include either one or both the ionosphere-over-time changes and the line-of-sight changes. No Taylor assumption is generally required.

[0044] Referring to Fig. 2, methods according to the present invention comprise a first step 101 of receiving from at least one, preferably each one of the GNSS satellites 20, 20’, 20” a first signal and a second signal having a carrier frequency different from a carrier frequency of the first signal. The first and the second signals are received by antenna 14 of GNSS receiver system 100. In a second step 102, a carrier phase difference between the first signal and the second signal is measured. Advantageously, step 102 comprises measuring a first carrier phase Φ ₁ of the first signal and a second carrier phase Φ ₂ of the second signal and taking the difference ΔΦ ₁₂ . The carrier phases are measured by tracking module 12.

[0045] In a third step 103, an ionosphere delay rate i _r for a reference frequency is determined. As will be described further below, advantageously, the ionosphere delay rate is determined on the basis of the measured carrier phase difference. Step 103 can comprise determining a thermal noise variance. In this case, the ionosphere delay rate is determined on the basis of the measured carrier phase difference and the determined thermal noise variance.

[0046] In a fourth step 104, a first indicator based on a quadratic expression is determined. The quadratic expression comprises a squared sum of a first term and a second term, wherein the first term comprises the measured carrier phase difference of step 102 and the second term comprises the ionosphere delay rate determined in step 103. Advantageously, the first indicator refers to the variance tracer, determined according to Eq. (7). This variance tracer is representative of the thermal noise variance.

[0047] In a fifth step 105, a test is performed on the basis of the indicator determined in step 104. If the test meets a predetermined requirement, e.g. its outcome is true, the method can determine that a disturbance is detected. If the test does not meet the requirement, the method can determine that a disturbance is not detected. The test can comprise comparing the indicator value to a threshold. If the indicator value exceeds the threshold, the test outcome can be true. Advantageously, a plurality of threshold levels are provided, each one corresponding to a different severity level of the disturbance. By comparing the indicator value with the plurality of threshold levels, the disturbance can be categorized into different severity levels.

[0048] It will be appreciated that steps 101-105 are repeated at each epoch. Additionally, steps 103 and 104 can be carried out simultaneously or in a different order than depicted in Fig. 2.

[0049] Referring again to Fig. 1 B, the positioning module 11 , which can comprise an RTK or PPP positioning engine, takes input from carrier phase measurements Φ ₁ and Φ ₂ , and code measurements Pi and P2 from n satellites (only measurements 200, 201, 202 for three satellites are shown) determined in tracking module 12. For each satellite the carrier phase measurements are input to the disturbance detection module 15, which is configured to determine a variance tracer, in particular a thermal noise variance as described herein. The dashed lines illustrate two possible outcomes of the test. If the variance tracer exceeds a threshold, the carrier phase measurement input to positioning module 11 may be blocked, or may be given reduced importance level (e.g. by appropriate weighting) or a warning may be given. If the thermal noise variance does not exceed a threshold, the phase measurement input is not blocked.

[0050] It will be appreciated that in step 105 multiple levels of thresholds may be defined, and based thereon, the carrier phase measurements may be (i) blocked, or (ii) may be reduced in importance level (e.g. by appropriate weighting) or (iii) a warning may be given but the carrier phase measurement may not be blocked. In addition, or alternatively to blocking the carrier phase measurement, the code measurements may be blocked based on the outcome of the test of step 105.

[0051] In one embodiment, the observable of Eq. (7) is used as a direct indication of the thermal noise variance. A disturbance can then be detected if the value of y _k meets a predetermined requirement, e.g. if it exceeds a threshold value. However, the observable of Eq. (7) can be noisy, and appropriate noise reduction techniques can be used. One advantageous noise reduction technique is Kalman filtering.

[0052] Advantageously, the observable of Eq. (7) is used in a Kalman filter to determine the thermal noise variance. In the example below, this Kalman filter is named the Thermal noise filter (TNF) and can be set up as follows.

[0053] A Kalman filter takes observations y _k at epoch k and its variancecovariance matrix Q _yk as input. These observations are related to the unknown states x _k by a linear relation H _k (the design matrix) as: E[y _k } = H _k x _k . (8)

Eq. (8) is called the observation equation, which will be described explicitly for the TNF below. The Kalman filter also takes the states of a previous epoch as input. The goal of a Kalman filter is to estimate the states x _k ; the estimated (filtered) states are denoted by

[0054] The equations for a Kalman filter can be found in many text books.

In the present description, a notation is used that is close to Teunissen, P.J.G., Dynamic data processing: Recursive least-squares, Series on mathematical geodesy and positioning, VSSD, ISBN 9065622179, 9789065622174, Delft. An overview of the equations relating to the TNF is shown on the left hand side of Fig. 3.

[0055] A Kalman filter consist of two steps: a prediction (time update), and a measurement update step. The following matrices are shown in Fig. 3: Φ _k is the transition matrix, is the system noise matrix, is the variance-covariance matrix of the predicted states is the variance-covariance matrix of the filtered states of the previous epoch, K _k is the Kalman gain matrix, and is the variance-covariance matrix of the predicted residuals v _k .

The thermal noise filter (TNF)

[0056] Eq. (5) shows the following thermal noise term:

The observable from Eq. (7) has expected values:

In Eq. (10), E[. } represents the mathematical expectation and D[. } the mathematical dispersion.

[0057] A suitable noise model for the mathematical dispersion can be set up. By way of example, two possible noise models are given below: and where is the thermal noise variance of Φ ₁₂ . The model of Eq. (12) is suitable for the GPS L2P signal which has correlations between the L1 and L2 carrier phase observations. The model of Eq. (11) is suitable for most other GNSS signals as long as Δt(k) is sufficiently large.

[0058] The covariance matrices have model parameters α(k) and correlation coefficient ρ(k). These can typically be modelled as function of elevation angle є(k) or carrier-to-noise ratio C/N ₀ . Simple piece-wise linear functions are sufficient and can be found when tracker covariance information is available from the receiver by a calibration with a data set under non-scintillation open-sky conditions. The α(k) parameter however will vanish from the observation equations below, and no calibration is required.

[0059] For the correlation coefficient ρ(k) the following model can be used for GNSS receivers of Septentrio, Belgium: and other models can be easily set up for other receivers as indicated above.

[0060] Using the model of Eq. (11), the following observation equation can be derived from Eq. (10) for the Kalman filter (TNF):

[0061] Using instead the model of Eq. (12), the following observation equation can be derived from Eq. (10) for the Kalman filter (TNF):

[0062] In both cases, the observable y _k of Eq. (7) is used as input for the

TNF, which estimates the scalar state .

[0063] In order to operate the TNF, information for the ionosphere delay rate İ _r (k) is required. The ionosphere delay rate can be determined according to known methods. By way of example, the ionosphere delay rate can be estimated by the positioning module 11 , e.g. if the ionosphere delay rates are states in the positioning algorithm implemented in positioning module 11 , or more likely if the ionosphere delays themselves are states in the positioning algorithm and the ionosphere delay rates are derived from them. Alternatively, according to an aspect of the present invention, the following relation, derived from Eq. (5), is used for determining the ionosphere delay rate:

It will be appreciated that only the carrier phase difference ΔΦ ₁₂ (k) needs to be measured I determined for determining the ionosphere delay rate. From Eq. (7) it can be seen that the thermal noise variance is advantageously determined based on only one type of input or measurement, i.e. based on the carrier phase difference ΔΦ ₁₂ (k) only.

[0064] Advantageously, a second Kalman filter, referred to herein as the

Ionosphere Delay Rate Filter (IDRF), is used for determining İ _r (k). As shown in Fig. 3, the filtered observations of the IDRF are fed to the TNF. Vice versa, the output of the TNF is used to compute the (scalar) variance-covariance matrix of the IDRF, i.e. two Kalman filters are used to determining the thermal noise variance. It will be appreciated that although the present description refers to the variance-covariance matrix, the matrix is in fact a scalar in the instant case.

The ionosphere delay rate filter (IDRF)

[0065] This second Kalman filter takes the right hand side of Eq. (16) as the observable y _k :

The corresponding variance is: where is the design matrix of the thermal noise filter from Eq. (14) or (15). Substituting the estimated from the thermal noise filter yields: (19)

The slant ionosphere delay rate İ _r (k) can be modeled as the only state. Alternatively, a more elaborate model may better serve the purpose of estimating İ _r (k) If the ionosphere delay is written as a product of a zenith ionosphere delay / ^z (k) and a mapping function m(є(k)) as function of elevation angle є(k), as e.g. described in Klobuchar, J.A., Ionosphere time-delay algorithm for single-frequency GPS users, Transactions on Aerosphere and Electronic Systems, 23(3), pp. 325-331 , then the following observation equation is obtained: where ṁ(є(k)) is the time-derivative of m(є(k)) and is the zenith ionosphere delay rate.

[0066] The filtered state forms an estimate of İ _r (k) which is input to the thermal noise filter (TNF). Fig. 3 shows the input ΔΦ ₁₂ (k) is used as input for both filters in Eq. (14, 15) and Eq. (17). The from the IDRF is taken in the next epoch as input for observation of the TNF according to Eq. (7):

[0067] In the IDRF, as shown in Fig. 3, from Eq. (19) takes input from

Q _k and the estimated thermal noise variance It will be appreciated that in the IDRF, is not strictly required, but may be used for statistical tests, the outcomes of which may be used to steer this filter. Usually not output in a Kalman filter, but in this embodiment added in the measurement update of the IDRF are the filtered observations .

[0068] It will be convenient to note that two representation forms of the measurement update equations are possible. In Fig. 3 both representation forms are used. One representation form of the measurement update equations is used on the left-hand side of Fig 3. The other representation form is used on the right-hand side (the first three equations differ). It will however be appreciated that either one of the two representation forms could be used in either filter. The matrices differ for both filters, but they are written here in the same notation.

[0069] The combination of two Kalman filters, i.e. one for determining the thermal noise variance and one for determining/updating the ionosphere delay rate allows to further decrease convergence time and improves response time of the disturbance detection compared to when only the TNF would be used.

[0070] The thermal noise variance under non-scintillation open-sky circumstances can be modelled as function of elevation angle or carrier-to-noise ratio C/N ₀ using the same kind of calibration as indicated above. As detector for scintillation one can then use the ratio of the square root of the estimated thermal noise variance and the modelled value. Under the above conditions this ratio, which is referred to herein as scintillation factor (SF) will then be around 1. A threshold level for the scintillation factor (e.g. 3--5) can then be defined above which scintillation is considered significant.

[0071] The above method can be extended to using carrier phase measurements taken from three different carrier frequencies, e.g. L1, L2, and L5 for GPS. The carrier phase measurement of the third frequency will add an extra observable in the IDRF, obtaining a 2-by-1 observable vector comprised of Eq. (17) and ΔΦ ₁₅ (k)/-μ ₁₅ Δ t(k), where μ ₁₅ is defined as in Eq. (3) based on the carrier wavelengths L1 and L5 (while p being e.g. defined based on L1 and L2). The design matrix is extended to of which the top row is as in Eq. (20). The TNF can be extended from one to three observables, adding to Eq. (7) the observables and The design matrix will depend on the noise model for the three frequencies.

Experiments

[0072] The TNF and IDRF Kalman filters were tested on data from a

GNSS receiver (PolaRx4, Septentrio, Belgium). Bounding was applied in both filters, meaning that values larger than a threshold are set to that threshold value. This allows the filters to stay within a normal range of values. In the IDRF a w-test was applied - The measurement update was discarded if this value exceeded the threshold. Table 1 sets out the filter components that were used for the TNF and IDRF.

Table 1 : TNF and IDRF filter components used in the experiments

[0073] Figs. 4 and 6 show the performance of the two filters of Fig. 3 under normal non-scintillation circumstances. In the figures, the Model M represents the values determined by Eq. (24), while the Estimate E represents the results of the TNF. The observations O represent the observable y _k of the TNF. Fig. 4 shows the results of the TNF in the upper graph, and the scintillation factor in the lower graph. It can be seen that the estimate E remains close to the model M and the scintillation factor remains quite stable about the value 1 , and no scintillation events occur in this time frame. Fig. 6 shows the results of the IDRF. From the figures it can be seen that despite the noisy behaviour of the observations O, the Kalman filters are capable or estimating the thermal noise and ionosphere delay rate.

[0074] Fig. 5 shows the TNF results for the same data as Fig. 4, where the thermal noise estimates have not been corrected by the ionosphere delay rates, i.e. a fixed value of the ionosphere delay rate of zero was applied (the IDRF was not used). It can be seen in the lower graph of Fig. 5 that the scintillation factor shows a higher variation. This plot illustrates the advantage of using the IDRF in methods described herein, since in Fig. 4 it removes the bump around 6.00-8.00 h which is visible in Fig. 5. [0075] Figs. 7A shows the TNF results and Fig. 7C the IDRF results of another measurement set including a scintillation event. These plots illustrate a successful detection of the scintillation event around 17.30-18.00 h.

[0076] Figs. 8A-B show the IDRF results for a satellite with both L2C and

L2P signals. The observation models of Eqs. (14) and (15) were used. Despite the lower noise of the L1-L2P ionosphere delay rate observations than those of L1-L2C, the Kalman filter results of the two signals are comparable, and the same ionosphere delay rate trend can be seen in both signals, evidencing the robustness of methods as described herein.

[0077] Aspects of the present disclosure are described in the following alphanumerically ordered clauses.

A1. Method of detecting a disturbance in a signal (21) of a GNSS satellite (20), comprising: receiving from the GNSS satellite at least a first signal and a second signal having a carrier frequency different from a carrier frequency of the first signal, measuring a carrier phase difference (ΔΦ ₁₂ ) between at least the first signal and the second signal, determining an ionosphere delay rate, determining a first indicator based on a quadratic expression comprising a squared sum of a first term and a second term, wherein the first term comprises the measured carrier phase difference and the second term comprises the ionosphere delay rate, detecting the disturbance if the first indicator meets a predetermined requirement.

A2. Method of clause A1 , wherein the first indicator corresponds to a thermal noise variance. A3. Method of clause A1 or A2, wherein the quadratic expression is given by: , or a scaled version thereof, wherein: ΔΦ ₁₂ represents the measured carrier phase difference, μ represents wherein λ ₁ and λ ₂ are the wavelengths of the carrier frequencies of the first signal and the second signal, respectively, and λ _r is the wavelength of a reference frequency, i _r represents the ionosphere delay rate for the reference frequency, and

At represents a time interval of discretization.

A4. Method of any one of the preceding clauses, wherein the disturbance is detected if the first indicator exceeds a predetermined threshold.

A5. Method of any one of the preceding clauses, wherein the first indicator is determined by a first Kalman filter.

A6. Method of clause A5, wherein the quadratic expression is used in the first Kalman filter as an observable (y _fe ).

A7. Method of any one of the preceding clauses, comprising using the measured carrier phase difference as an input for determining the ionosphere delay rate.

A8. Method of clause A7, comprising determining a second expression as: or a scaled version thereof, wherein: ΔΦ ₁₂ represents the measured carrier phase difference, p represents , wherein λ ₁ and λ ₂ are the wavelengths of the first signal and the second signal, respectively, and λ _r is the wavelength of a reference frequency,

Δt represents a time interval of discretization, and using the second expression as an input for determining the ionosphere delay rate.

A9. Method of clause A7 or A8, wherein the ionosphere delay rate is determined by a second Kalman filter, wherein the measured carrier phase difference is used as an input to the second Kalman filter. A10. Method of clause A9 in conjunction with clause A5 or A6, wherein an estimated state of the second Kalman filter is used for determining an observable of the first Kalman filter.

A11. Method of clause A10, wherein an output of the first Kalman filter is used to determine a covariance in the second Kalman filter.

A12. Method of determining navigation information based on GNSS navigation signals, comprising: receiving GNSS navigation signals from a plurality of GNSS satellites, the GNSS navigation signals comprising, for each of the plurality of GNSS satellites at least a first signal and a second signal having a carrier frequency different from a carrier frequency of the first signal, for each of the GNSS satellites, performing carrier phase measurements on the first signal and the second signal, performing the method of any one of the preceding clauses using the carrier phase measurements to detect a disturbance in signals of the respective GNSS satellite, if the disturbance is detected, determining the navigation information comprising reducing a weight of or discarding a measurement of the signals in respect of which the disturbance is detected.

A13. Method of clause A12, wherein reducing a weight of or discarding a measurement of the signals in respect of which the disturbance is detected comprises reducing a weight of or discarding one or more of a carrier phase measurement and a code measurement.

A14. GNSS receiver system (100) comprising an antenna (14) for receiving a first signal and a second signal of a GNSS satellite and computing means (10) for processing the first signal and the second signal, wherein the computing means is configured to carry out the method of any one of the clauses A1 to A11, preferably comprising a navigation unit (11) configured to carry out the method of clause A12 or A13.

A15. Ionospheric scintillation monitoring device, comprising an antenna for receiving a first signal and a second signal of a GNSS satellite and a signal tracking unit configured to carry out the method of any one of the clauses A1 to A11.