GNSS positioning system and method using multiple processing threads

A plurality of GNSS satellite signals feeds multiple signal processing engines, each operating in certain processing mode including carrier smoothed pseudorange positioning, precise point positioning (PPP), pseudorange differential (DGNSS), carrier phase differential (RTK). Each processing engine (or processing thread of the same engine) runs the same unified numerical algorithm and uses the same or different sets of parameters. All engines can use the same set of signals, or the set of signals can be split into non-intersecting subsets, or the set of signals can be split into the overlapping subsets. Each engine produces estimates of certain parameters, namely carrier phase ambiguities and ionospheric delays for each satellite. These estimates are then combined into a resulting estimate which in turn is used for calculation of the final position reported by the receiver.

BACKGROUND OF THE INVENTION

Field of the Invention

The present invention relates to a satellite navigation systems, and, more particularly, to more accurate determination of phase ambiguity in satellite navigation signals using multiple processing threads.

Description of the Related Art

The number of such signals can exceed several dozen. The following fundamental set of observables is used:

where the following notations are used (see [1, Chapter 6, 7]):k is the index of the station;p is the number of the satellite (out of K total satellites from which signals are received). It is assumed that all satellites available for tracking in the receiver station k are ordered and this ordering number includes the satellite system. For example, satellites with numbers from 1 to 5 are GPS, satellites with numbers from 6 to 11 are GLONASS, satellites with numbers from 12 to 15 are BEIDOU, etc.;b is the frequency band (for example L1, LP2, L2C, L5, E6, B1, and so on);(p,b) the pair of indices indicating the signal of the satellite p at the frequency band b;t is the current time; the term ‘epoch’ is also used to denote the current discrete time instant;c is the light speed;fbpis the frequency of the signal corresponding to the satellite p and the frequency band bdtk(t) is the current clock bias of the station k;dtp(t) is the current clock bias of the satellite p;Ik,bp(t) is the ionospheric delay affecting the signal (p, b) received by the station k. Thus Ik,L1p(t) is related to the L1 band. Basically, the ionospheric delay depends on the position of the station, position of the satellite, frequency of the signal, and the total electronic content (TEC) corresponding to the time instant t;Tkp(t) is the tropospheric delay. In contrast to the ionospheric delay it doesn't depend on the signal frequency and is called ‘non-dispersive’ delay.Pk,bp(t) and φk,bp(t) are pseudorange and phase measurements respectively;Mk,b,Ppand Mk,b,φpare code and phase multipath, affecting pseudorange and carrier phase measurements respectively;

Nk,bp(tCS,k,bp) carrier phase ambiguity corresponding to the signal (p,b) received by the station k. Note that it corresponds to the last detected cycle slip and therefore it explicitly dependents on the time instant tCS,k,bp, when the cycle slip was detected. The carrier phase ambiguity remains unchanged until the cycle slip occurs.ρkp(t) is the true topocentric range between the satellite and the receiver station. The measurement equations (1) and (2) relate this quantity with pseudorange and carrier phase observables, respectively;Quantities Db,φp, Db,Pp, εk,b,Pp(t) and εk,b,φpdenote hardware biases and noise. D stands for satellite-related bias, while d denotes the receiver-related bias. Biases reflect a systematic hardware component of the measurement error. They are constant or slow varying, in contrast to the noise component of the error having the stochastic nature. Satellite position error is part of the slow varying bias. It is also called ‘the ephemerides’ error.

Thus, the receiver position is measured by the pseudorange and carrier phase observables for the plurality of satellites. Error components, including biases and noise, affecting the observable equations (1) and (2), prevent a direct solution for the receiver antenna position.

Carrier phase measurements are much more precise, compared to the pseudorange measurements, since the carrier phase noise has standard deviation in the centimeter or even millimeter range, while the standard deviation of the pseudorange measurements is usually of the meter or decimeter range. On the other hand, the carrier phase measurement is affected by the carrier phase ambiguity, which is an unknown integer valued quantity.

Thus, elimination of measurement errors is necessary for precise positioning. To achieve a high precision in position determination, different methods of errors mitigation are applied. For example, tropospheric errors can be precisely modeled and compensated in observables of equations (1) and (2). Ionospheric errors can be estimated along with other unknowns. Noise is easily filtered.

Errors common to two receivers, like clock and hardware biases of the satellite, can be compensated in a difference between two receivers. Usually one of receivers occupies a known position, while an antenna of another receiver is attached to the object to be located. The first receiver is called ‘the base’ while another receiver is called ‘the rover’. The processing mode involving calculation of the across-receiver difference (also called the ‘first difference’) is referred to as differential GNSS processing or DGNSS. The DGNSS processing is performed in real time and includes not only pseudoranges but also carrier phase observables, and is referred to as real time kinematic (RTK) processing.

Another sources of errors partially eliminated by across-receiver differences are ionospheric delay and ephemerides error. The closer the rover is to the base, the better is compensation of the ionospheric and ephemerides error.

For two stations k and l the across-receiver differences of pseudorange and carrier phase measurements can be written as

Another way for error mitigation includes using a precise satellites clock and precise ephemerides. They are available through a dedicated real time service. Precise point positioning (PPP) allows to achieve the centimeter level position with only one rover receiver, provided satellite clock and precise satellite position are available. The base station is not necessary in this case.

Finally, if neither base station, nor precise clock and ephemerides are available, the quality of the standalone position can be improved if carrier phase ambiguity and ionospheric delay are estimated, along with position, using broadcast ephemerides. The corresponding processing mode is equivalent smoothing of pseudoranges using carrier phase measurements or, in short, carrier phase smoothing of code pseudoranges, see [2]. Whatever processing mode is used, the linearization and filtering algorithms are used for recursive estimation of unknown position, carrier phase ambiguity, and ionospheric delay.

Two last quantities are undifferences for carrier phase smoothing and PPP processing modes. For DGNSS and RTK processing modes the carrier phase ambiguity and ionospheric delay are across-receiver differenced. The general form of the linearized navigation model has the following form (see [1, Chapt. 7])
bP(t)=Adx(t)+eξ(t)+Γi(t)+dP(5)
bφ(t)=Λ−1Adx(t)+Λ−1eξ(t)+n−Λ−1Γi(t)+dφ(6)

Let n be total number of satellite signals, including different satellite systems, different satellites, different frequency bands.

In the following description, all vectors are represented by columns, and the superscript symbolTdenotes the matrix transpose. RNis the N-dimensional Euclidean space. Given a linearization point x0(t)∈R3, notations used in equations (5) and (6) are as follows:bP(t)∈Rnis the n-dimensional vector of pseudorange residuals calculated at the linearization point;bφ(t)∈Rnis the n-dimensional vector of carrier phase residuals calculated at the linearization point;e=(1, 1, . . . , 1)T∈Rnis the vector consisting of all “ones”;dx(t)∈R3is the correction to the linearization point. Thus, the corrected position is calculated as x(t)=x0(t)+dx(t);ξ(t) is the arbitrary varying rover clock shift; it is undifferenced for the standalone and PPP processing modes, and it is across-receiver differenced in the DGNSS and RTK processing modes;Λ is the n-dimensional diagonal matrix with wavelengths λbp=c/fbpin the main diagonal. Each wavelength corresponds to the specific signal (p,b);A is the n×3 matrix of directional cosines;Γ is the n×n diagonal matrix with quantities (fL1p/fbp)2in the main diagonal;i(t)∈Rnis the vector of ionospheric delays related to the L1 frequency band. Ionospheric delays are undifferenced for the standalone and PPP processing modes, and they are across-receiver differenced in the DGNSS and RTK processing modes;n∈Rnis the vector of carrier phase ambiguities related to the L1 frequency band. Ambiguities are undifferenced for the standalone and PPP processing modes, and they are across-receiver differenced in the DGNSS and RTK processing modes;dp∈Rnand dφ∈Rnare vectors of pseudorange and carrier phase receiver hardware biases.

Consideration of pseudorange hardware biases leads to a necessity to consider the plurality of signals the receiver is able to track. In the case of a multi-frequency and multi-system receiver supporting the following bands:

the signals (L1 GPS, L1 Galileo, L1 SBAS, L1 QZSS), (L2 GPS, L2 QZSS), (L5 GPS, E5a Galileo, L5 SBAS, L5 QZSS), (E6 Galileo, E6 QZSS), respectively, can share the same hardware channel and therefore will be affected by the same hardware bias, as noted in [1, Chapt. 7]. Note that the biases vector dPand the clock shift variable ξ(t) appear as a sum in equation (5). This means that one of the biases, say dL1,G/E/S/Q,P, can be combined with ξ(t), while others can be replaced with their difference with dL1,G/E/S/Q,P. Thus, new bias variables appear:
η1=L2,G/Q,P−dL1,G/E/S/Q,P,η2=dL1,R,P−dL1,G/E/S/Q,P,η3=dL2,R,P−dL1,G/E/S/Q,Pη4=dL5,G/E/S/Q,P−dL1,G/E/S/Q,P,η5=dE5b,E,P−dL1,G/E/S/Q,P,η6=dE6,E/Q,P−dL1,G/E/S/Q,Pη7=dB1,B,P−dL1,G/E/S/Q,P,η8=dB2,B,P−dL1,G/E/S/Q,P,η9=dB3,B,P−dL1,G/E/S/Q,P(7)

This representation can be referred to as establishing a bias datum.

In one possible embodiment, the linearized equations (5) can now be expressed in the form
bP(t)=Adx(t)+eξ(t)+Γi(t)+Wηη  (8)

The bias vector η has the appropriate dimension mη. Note again that in this exemplary embodiment we follow notations introduced in [1], which is incorporated herein by reference in its entirety.

The bias vector η is three-dimensional (mη=3) for dual-band and dual-system GPS/GLONASS receivers, as only biases η1, η2, η3are presented among all possible biases listed in (7). In the case of the multi-band, multi-system receiver, the dimension of the vector η can be large. It is one-dimensional in the case of dual-band GPS-only receivers or single band GPS/GLONASS receivers.

The Wηis referred to as bias allocation matrix and has dimensions n×mη. It allocates a single bias, or none, to a certain signal. No bias is allocated to the signal corresponding to the GPS, Galileo, SBAS, or QZSS systems and b=L1 because we combined the bias dL1,G/E/S/Q,Pwith the clock bias ξ(t). In this case, the row of Wηconsists of zeroes.

Consider, for example, a dual-band GPS/GLONASS receiver. Suppose it tracks six GPS satellites and six GLONASS satellites. The total number of dual-band signals is 24. Let the signals be ordered in the following way: six GPS L1, six GPS L2, six GLONASS L1, and six GLONASS L2 signals. The biases allocation matrix presented in the linearized single difference pseudorange equation (8) takes the following form:

Further, the real-valued carrier phase ambiguities (also called float ambiguities) are combined with biases dφ, while pseudorange biases are considered as a real valued constant unknown parameter. Thus, after combination of carrier phase ambiguities with carrier phase biases, the equation (6) takes the form
bφ(t)=Λ−1Adx(t)+Λ−1eξ(t)+n−Λ−1Γi(t)  (10)

Note that the noise component is omitted in equations (8), (10) for the sake of brevity.

Accordingly, there is a need in the art for a more accurate determination of phase to, thereby, enable more accurate determination of position.

SUMMARY OF THE INVENTION

Accordingly, the present invention is directed to a plurality of GNSS satellite signal feeds to multiple signal processing engines, with each operating in a specific processing mode, such as carrier smoothed pseudorange positioning, precise point positioning (PPP), pseudorange differential (DGNSS), and carrier phase differential (RTK). Each processing engine (or processing thread of the same engine) runs the same numerical algorithm and uses the same or different sets of parameters. All engines can use the same set of signals, or the set of signals can be split into non-intersecting subsets, or the set of signals can be split into the overlapping subsets. Each engine produces estimates of certain parameters, namely carrier phase ambiguities and ionospheric delays for each satellite. These estimates are then combined into resulting estimate which in turn is used for calculation of the final position reported by the receiver.

DETAILED DESCRIPTION OF EMBODIMENTS OF THE INVENTION

Thus, equations (8) and (10) above are the basis for construction of the unified numerical scheme covering the following processing modes:

Standalone with carrier phase smoothing of pseudoranges,

As explained, different processing modes use the same set of equations presented in the generalized form. Therefore, the general purpose unified numerical scheme of the processing engine can be used. Each processing mode is thus performed by the certain embodiment of the same processing engine. The same source code, such as a C++ source code template in one embodiment, can be used to generate multiple processing threads working in parallel. We call them threads, since they can be implemented as parallel processing threads being governed by an appropriate operating system, such as LINUX, RTEMS, or eCos.

Each processing mode can be presented by one or more multiple threads, running the same algorithm with different parameters sets reflecting different assumptions about errors to be modeled or different statistic assumptions about stochastic noise component of the error. Implementation of this idea as it applies to only RTK is the subject of the patent [3]. In the present discussion, we explore a more general idea using multiple threads for support of all possible processing modes (and not only RTK) and describe more general voting procedure allowing to combine phase ambiguities and ionosphere delay estimates produced by different threads, not just final RTK positions produced by different threads. This allows for deeper and tighter integration of results produced by different parallel threads. It also allows for seamless transition between processing modes, i.e., standalone, PPP, DGNSS, RTK float, and RTK fixed.

Following the idea of using the generalize processing scheme, all unknowns are grouped into following groups:

Arbitrary varying variables dx(t) and ξ(t),

Slow varying variables subjected to certain dynamic model i(t),

Constant (time invariant) variables n and η.

To estimate all three groups of variables in real time, the recursive estimator must be applied. It can be implemented as a sort of Kalman filter, particle filter, another recursive estimator which is convenient for the real time implementation. In one exemplary embodiment, the general purpose processing thread can run the recursive estimator described in [1, Chapt. 3.5].

To specify dynamic model that governs variation of i(t) let us suppose that the vector i(t) is restricted by equations

with Δt being the across-epoch time difference, and τiis the correlation time reflecting the rate of variation of the ionospheric delay in time. A typical value for τiis 600 seconds, provided no ionospheric scintillations happened. Otherwise, if detection criteria (not described here) signal a presence of fast variation of ionosphere, the instant value of τican be temporary scaled down until the fast variations end. The white noise εi(t) has a variance σi2, which provides the variance of the across-receiver ionosphere satisfying the condition
|ip(t)|≈s×10−6×∥baseline∥  (13)

with a certain scale factor s, which takes values 1 to 5, depending on the current solar activity. It follows from (11) that ∥ip(t)∥2=∥γiip(t−1)+εi(t)∥2. Then, assuming that the stochastic process ip(t) is stationary and ip(t) does not depend on εi(t), we take a mean value of both sides of the last equality. We obtain
σi2=(1−γi2)×(s×10−6×∥baseline∥)2(14)

Let us define the measurement and dynamic models. The measurement model combines equations (8) and (10).

The group of constant variables n and η form the vector

The group of arbitrary varying variables dx(t) and ξ(t) form the vector

Then the algorithm described in [1, Table 3.5.1] can be applied.

Definition of Multiple Processing Threads

The same general purpose numerical scheme or the engine template is used to generate multiple instances as shown inFIG. 1. In this exemplary implementation, several Standalone, several PPP and several RTK engines are implemented.

Each engine (indexed by the variable j ranging from 1 to m) calculates its own vector of carrier phase ambiguities njand its own vector of ionospheric delays estimate ij(t). Each estimate is accompanied by its covariance matrix Cnjand Cijrespectively. The resulting ambiguity and ionospheric delay vectors can be calculated in the following two ways:

A. Having covariance matrices Cnjand Cijwe select the thread corresponding to the minimum quantity tr(Cnj)+tr(Cij) where tr(·) denotes the trace of the matrix.

B. Weighted estimates are calculated according to the formulas

Estimates of the vectors n and i(t) produced in the way A or B according to (17) and (18) are then substituted into the equation (10), which is finally solved for time-varying variables

x_⁡(t)=(d⁢⁢x⁡(t)ξ⁡(t))
using minimization of least squares, according to the expression
x(t)=Ā+(bφ(t)−n+Λ−1Γi(t))  (19)
Ā=Λ−1[A|e]  (20)

where the symbol Ā+denotes the matrix pseudoinverse.

Using different sets of signals in the processing threads and using different sets of the processing parameters allows introduction of randomization into the set of a priori hypotheses placed into the measurement models, which in turn allows for faster ambiguity resolution and better averaging of the measurement errors.

Computational Complexity

The computational complexity of the processing thread can be roughly estimated as
T1=Cn3(21)

where C is a constant specific to the certain processing algorithm and n is number of signals. This estimate follows from the fact that most complex part of the numerical scheme is the matrix inverse of Cholesky factorization, or another factorization.

If all threads use the same set of signals, then the overall computational cost is estimated as
Tm=mT1(22)

Assume now that the set of n signals is divide into m threads. The resulting complexity can then be estimated as

which means that the computational complexity can be even reduced if the set of signals is appropriately distributed among processing threads.

In a general case, the set of signals is not necessarily split into m not intersecting subsets. Instead, subsets can overlap. Thus, resulting computational costTcan be estimated as

As a further option, each engine thread can use different settings of parameters, even if they run the same processing mode. Also, each engine thread can use different sets of signals. The processing threads can use different a priori assumptions and hypotheses about:

a) state of the ionosphere

b) level of the local multipath conditions

c) level of the measurement noise

The voting procedure judges which assumption is likely true. Alternatively, a weighting procedure aggregates all or some of the assumptions.

FIG. 2illustrates an exemplary navigation receiver on which the invention may be implemented. As will be appreciated by one of ordinary skill in the art, the various blocks shown inFIGS. 1-2can be implemented as discrete components, as an ASIC (or multiple ASICs) and/or as software running on a processor.

Having thus described the different embodiments of a system and method, it should be apparent to those skilled in the art that certain advantages of the described method and apparatus have been achieved. It should also be appreciated that various modifications, adaptations, and alternative embodiments thereof may be made within the scope and spirit of the present invention. The invention is further defined by the following claims.

REFERENCES (INCORPORATED HEREIN BY REFERENCE IN THEIR ENTIRETY)