Method and apparatus for single epoch position bounding

Abstract

The invention relates to a method and apparatus for single epoch position bounding. A method for determining a protection level for a position estimate using a single epoch of GNSS measurements, the method comprising: specifying a prior probability density of states x P(x); specifying a system model h(x) relating states x to measured observations z; quantifying a quality metric q associated with the measurements; specifying a non-Gaussian residual error probability density model f(r|θ,q) and fitting the model parameters θ using a set of experimental data; defining a posterior probability density P(x|z,q,θ); estimating states x; and calculating the protection level by integrating the posterior probability density P(x|z,q,θ) over states x.

Description

TECHNICAL FIELD

[0001] The present disclosure relates generally to satellite-based navigation, and more particularly, to a method of determining a protection level of a position estimate using a single epoch of global navigation satellite system (GNSS) observation data and / or measurements from other related sensors and an operating device thereof. BACKGROUND

[0002] GNSS receivers receive satellite signals transmitted from one or more GNSS satellite constellations and use information contained in the satellite signals to estimate its position, etc. Generally, this estimation alone cannot provide sufficient accuracy and reliability due to, for example, multipath interference suffered by the satellite signals. The consequences of providing a position estimate with errors exceeding a tolerable range can be significant. Therefore, there is a need to determine a probability bound on the error associated with a position estimate.

[0003] The probability bound on the error associated with a position estimate can be determined by modeling errors of a set of GNSS observations. However, when the observations are not independent of each other, modeling of the observation errors is difficult, resulting in difficulty in determining the probability bound on the error associated with a position estimate. SUMMARY

[0004] According to some embodiments of the present disclosure, a method of determining a protection level of a position estimate using a single epoch of global navigation satellite system (GNSS) measurements is provided. The method includes pre-specifying a prior probability density of a state x P(x); pre-specifying a system model h(x) relating the state x to the measured observations z; quantifying a quality metric q associated with the measurements during operation; pre-specifying a non-Gaussian residual error probability density model f(r|0,q) and fitting the model parameters 0 offline from the prior; defining a posterior state probability density P(x|z,q,0) during operation; estimating the state x during operation; and computing the protection level by integrating the posterior probability density P(x|z,q,0) over the state x during operation.

[0005] Optionally, the system model can be specified as z = h(x) + r, where r is a residual error.

[0006] In some arrangements, quantifying the quality metric q associated with the measurements can include using a random sample consensus (RANSAC) technique to exclude one or more outliers of the measurements.

[0007] The non-Gaussian residual error probability density model can be a student-t distribution function. Specifying the residual error probability density model f(r|0,q) can further include using the student-t distribution function to identify a model of pseudoranges as:

[0008]

[0009] where r is a residual error, v is a degree of freedom parameter, σ is a scaling parameter, and Γ is a gamma function.

[0010] Optionally, the non-Gaussian residual error probability density model can be a student-t distribution function. Specifying the residual error probability density model f(r|θ, q) can further include identifying the model of the carrier phase using the student-t distribution function as:

[0011]

[0012] where r is a residual error, v is a degree of freedom parameter, σ is a scaling parameter, w is a weight, and Γ is a gamma function.

[0013] Fitting the model parameters θ using a set of experimental data can include defining the degree of freedom parameter and the scaling parameter using unconstrained student-t parameters; defining the unconstrained student-t parameters as a function of a mass metric; and finding a maximum likelihood value of the unconstrained student-t parameters.

[0014] In an example method, fitting the model parameters θ using a set of experimental data can include defining the degree of freedom parameter and the scaling parameter and a uniform weight parameter using unconstrained parameters; defining the unconstrained parameters as a function of a mass metric; and finding a maximum likelihood value of the unconstrained parameters.

[0015] Defining the posterior probability density can be based on:

[0016]

[0017] where i is a natural number.

[0018] Optionally, integrating the posterior probability density over states can include applying Markov Chain Monte Carlo (MCMC) numerical integration using a plurality of interacting chains with modified density functions.

[0019] In some arrangements, integrating the posterior probability density over the state can include: performing a first round of sampling by extracting a first set of samples from an unbiased distribution and estimating a first position of a quantile on the first set of samples; performing a second round of sampling by extracting a second set of samples from a first constrained distribution determined based on the first position of the quantile; estimating a second position of the quantile on the second set of samples, wherein the second position of the quantile is used to set a second constrained distribution for a third round of sampling; repeating the sampling using the constrained distribution set in the previous sampling until an nth round of sampling is completed, where n is a natural number greater than 2; determining whether the samples of the n rounds of sampling are sufficient to compute the protection level; and combining the samples of the n rounds of sampling with appropriate weights when the n rounds of samples are sufficient to compute the protection level.

[0020] The method can further include propagating the nearest single epoch bound forward to the current time using the delta phase error distribution or data obtained by the one or more sensors.

[0021] The propagating can further include modeling the error distribution on the delta phase using a non-Gaussian distribution function and numerically evaluating the bound on the delta position using Markov Chain Monte Carlo techniques or modeling the error distribution on the delta position using a Gaussian overshoot method.

[0022] Optionally, the state x can be defined as:

[0023]

[0024] where p is the position of the rover, c is the speed of light in vacuum, Δdt r is the difference in receiver clock bias between the rover and the base station, Δk p,r is the difference in receiver code instrument delay on the pseudorange, Δk L,r is the difference in receiver code instrument delay on the carrier phase, and ΔT zenith is the difference in zenith tropospheric delay.

[0025] According to some embodiments of the disclosure, there is provided a method of quantifying the quality of GNSS measurements. The method comprises the steps of: performing GNSS measurements in a window around an epoch of interest; determining changes in the GNSS measurements; identifying a consensus solution for the changes in position and clock bias; and characterizing the quality of the GNSS measurements by deviation from the consensus. The GNSS measurements can comprise at least one of phase measurements, pseudorange measurements, and Doppler measurements.

[0026] Optionally, identifying a consensus solution for the changes in position and clock bias can further comprise: modeling the GNSS measurement data using a system representing position and clock bias offsets relative to a first epoch of the window; and solving the system as a weighted linear least squares problem.

[0027] Solving the system as a weighted linear least squares problem can further comprise introducing a loss function on the residual errors as:

[0028]

[0029] where p(s) is a loss function, w i is a weight for each measurement, r i is a residual error for each measurement, and i is a natural number.

[0030] Optionally, the loss function can be a Huber loss function:

[0031]

[0032] where d represents a threshold parameter.

[0033] According to some embodiments of the present disclosure, there is also provided a device configured to determine a protection level of a position estimate using a single epoch of GNSS measurements. The device comprises a receiver configured to receive GNSS signals and to process the received signals to generate measurements, and a processor configured to: specify a prior probability density P(x) of a state x; specify a system model h(x) relating the state x to the measurements z; quantify a quality metric q associated with the measurements; specify a non-Gaussian residual error probability density model f(r|0,q) and fit the model parameters 0 offline from the prior; define a posterior probability density P(x|z,q,0); estimate the state x using, for example, the posterior probability density P(x|z,q,0); and compute the protection level by integrating the posterior probability density P(x|z,q,0) over the state x.

[0034] Optionally, the processor can be further configured to communicate with an external sensor to obtain position change data tracked by the sensor to determine the bound propagation.

[0035] According to some embodiments of the present disclosure, there is also provided a non-transitory computer readable medium having stored therein instructions which, when executed by a processor, perform a method of determining a protection level of a position estimate using a single epoch of GNSS measurements, the method comprising the steps of: specifying a prior probability density P(x) of a state x; specifying a system model h(x) relating the state x to the measurements z; quantifying a quality metric q associated with the measurements; specifying a non-Gaussian residual error probability density model f(r|0,q) and fitting the model parameters 0 offline from the prior; defining a posterior probability density P(x|z,q,0); estimating the state x using, for example, the posterior probability density P(x|z,q,0); and computing the protection level by integrating the posterior probability density P(x|z,q,0) over the state x. Attached Figure Description

[0036] Figure 1 This is a block diagram of an apparatus consistent with some embodiments of this disclosure.

[0037] Figure 2 This is a flowchart illustrating a method for determining the protection level of a location estimate that is consistent with some embodiments of this disclosure.

[0038] Figure 3 This is a flowchart illustrating a method for identifying outliers in GNSS measurements that is consistent with some embodiments of this disclosure.

[0039] Figure 4 This is a flowchart illustrating a method for quantifying quality metrics associated with GNSS measurements that is consistent with some embodiments of this disclosure.

[0040] Figure 5 A comparison is shown of the residual error distribution of the observations, the Gaussian distribution fitted to the residual error distribution of the observations, and the student-t distribution fitted to the residual error distribution of the observations, consistent with some embodiments of this disclosure.

[0041] Figure 6 This is a schematic diagram illustrating a multi-peak sampling method using parallel interaction chains, consistent with some embodiments of this disclosure.

[0042] Figure 7 This is a flowchart illustrating a method for performing umbrella sampling consistent with some embodiments of this disclosure.

[0043] Figure 8A and Figure 8B The Stanford plot shown illustrates the bounding half-width and true error obtained from experimental data in the along-track direction, consistent with some embodiments of this disclosure. Figure 8C and Figure 8D The Stanford plot shows the bounding half-width and true error obtained from the experimental data in the cross-track dimension. Figure 8E and Figure 8F The Stanford plot shows the bounding half-width and true error obtained from the experimental data in the vertical dimension. Detailed Implementation

[0044] Reference will now be made in detail to exemplary embodiments, examples of which are illustrated in the accompanying drawings. The following description refers to the accompanying drawings, wherein, unless otherwise indicated, the same reference numerals in the different drawings denote the same or similar elements. The implementations set forth in the following description of the exemplary embodiments do not represent all implementations consistent with this disclosure. Rather, they are merely examples of systems, apparatuses, and methods consistent with aspects of this disclosure set forth in the appended claims.

[0045] A GNSS receiver receives satellite signals transmitted from one or more GNSS satellite constellations through an antenna and uses information contained in the satellite signals to estimate its position. Many GNSS applications require a certain level of accuracy and reliability of the estimated position. For example, for autonomous vehicles that rely on trustworthy global position information, the consequences of providing a position estimate (or any other parameter related to the estimate) with an error exceeding a tolerable range can be significant.

[0046] The error associated with a position estimate along a particular direction can be quantified using the concept of a "protection level" in that direction. A protection level can be defined as a probabilistic bound on the error of a position estimate, specifying a bounding region and suggesting that the true position is outside that region with a probability lower than a given value. The region can be provided in the form of an interval, a circle, or a rectangle, or any other shape. For example, a protection level can be provided as a circle around a horizontal position estimated by a GNSS receiver. The protection level can be used to determine when to alert a user that the receiver is unable to provide a position estimate with sufficient accuracy. For example, if the protection level exceeds a predetermined alert limit, an alarm can be triggered so that the client system can make appropriate adjustments. The ability of a system to provide a warning to a user when the system position estimate is not reliable can be described using the concept of "integrity." If the protection level does not exceed the alert limit, the user can have confidence in the accuracy of the position estimate.

[0047] Determining a protection level from a set of GNSS observables (e.g., pseudoranges or carrier phases or Doppler) is a statistical problem that depends primarily on the error probability distribution of the observables. Given a set of known error distributions of the observables, Bayes' method can be used to compute the posterior probability density (i.e., the probability density over the position). The protection level can then be determined by integrating the posterior probability density.

[0048] However, when the observables are statistically dependent on each other, the observable errors need to be modeled in a joint distribution. This makes the error distribution more difficult to model and the posterior probability density more difficult to determine, resulting in difficulty in computing the protection level.

[0049] Embodiments of the present disclosure provide a method for determining a protection level of a position estimate using a single epoch of GNSS measurements. An epoch can be represented as a set of measurements from a GNSS receiver at a single point in time. In embodiments, a prior probability density of states is specified; a system model linking states to measurements of observables is specified; and a non-Gaussian measurement error probability density model is specified based on offline-determined fitting model parameters. In embodiments, a quality metric is derived and outliers of measurements are excluded or downweighted using a window-based approach or RANSAC techniques. A posterior probability density is defined and the protection level associated with the position estimate is determined by integrating the posterior probability density using Markov Chain Monte Carlo (MCMC) or importance sampling methods. In embodiments, limit propagation is implemented using a delta phase error distribution or data obtained by a sensor.

[0050] The embodiments disclosed herein have one or more technical effects. By considering only the observables from a single epoch of GNSS measurements, a major source of non-independence (e.g., time correlation) is removed and error distributions on individual signals are modeled as one-dimensional probability density functions separately, resulting in a fast and rigorous determination of the protection level. The use of a non-Gaussian error probability density model allows for proper consideration of tail probability densities, ensuring accuracy of the processing. Excluding or downweighting outliers of measurements allows for enhanced accuracy and consistency of the determination. Integrating the posterior probability density using multiple MCMC chains with modified density functions allows for consideration of tail probability densities on multi-modal distributions, further enhancing accuracy of the processing. Performing limit propagation allows the system to maintain high availability and low latency.

[0051] Figure 1 is a block diagram of an apparatus 100 consistent with some embodiments of the present disclosure. Referring to Figure 1 , the apparatus 100 can be installed in a mobile vehicle or in a fixed location. The apparatus 100 can take any form including, but not limited to, a laptop computer, a GNSS receiver, a wireless terminal including a mobile phone, a wireless handheld device, or a wireless personal device, or any other form. The apparatus 100 includes an antenna 102, a receiver 104 coupled to the antenna 102, a processor 106, a memory 108, a local clock 110, and an input / output device 112.

[0052] The antenna 102 is configured to receive GNSS signals. The GNSS signals can be received from a single satellite or multiple satellite signals transmitted in multiple frequency bands from multiple satellites, respectively. The GNSS signals can also include signals originating from one or more virtual sources that are derived from reflected and / or scattered satellite signals. However, the signals received by the antenna 102 are not limited to satellite signals and can be any electromagnetic waves transmitted from any source (e.g., wireless cellular signals). The antenna 102 can be any type of antenna such as a patch antenna, a helical antenna, a crossbow antenna, orthogonally placed monopole antennas, etc. The antenna 102 can be an antenna array.

[0053] The receiver 104 is coupled to the antenna 102 and is configured to receive GNSS signals via the antenna 102. The GNSS signals received by the antenna 102 can be transmitted to the receiver 104 via a coaxial RF cable or any other cable suitable for transmitting RF signals. In embodiments, the receiver 104 can be part of a transceiver modem that includes a transmitter configured to transmit data to an external device.

[0054] The receiver 104 can perform operations on the received GNSS signals such as amplification, filtering, mixing, and digitization. The receiver 104 can also process the digitized signals and generate measurements such as pseudorange measurements and carrier phase measurements. The receiver 104 can also determine measurement quality indicators (also referred to as quality metrics). The quality metrics can provide information about the quality of the environment in which the measurements were performed. Examples of quality metrics include (but are not limited to) carrier-to-noise density ratio, satellite elevation angle, time-to- lock, time-to-code, and cycle slip. The accuracy of the measurements made by the receiver 104 can be related to the quality metrics associated with the measurements. The receiver 104 is also configured to communicate with the clock 110.

[0055] The processor 106 can include one or more special purpose processing units, application specific integrated circuits (ASICs), field programmable gate arrays (FPGAs), or various other types of processors or processing units. The processor 106 can receive the results of the pseudorange and carrier phase measurements from the receiver 104 and the quality metrics associated with the measurements and also process the information to estimate the current position of the device 100 and a protection level of the position estimate. The processor 106 can also provide its own measurement quality metrics. In embodiments, the processor 106 can include a navigation filter (e.g., a Kalman filter or a recursive LS filter) for determining the protection level. The processor 106 can be coupled to external motion sensors such as accelerometers, gyroscopes, or wheel speed sensors to obtain position tracking information from the motion sensors and aid in determining the protection level. The processor 106 is also configured to communicate with the input / output device 112 and the memory 108.

[0056] In embodiments, the receiver 104 can include a built-in processor (not shown) that performs all or part of the functions of the processor 106. In embodiments, the built-in processor of the receiver 104 can be a front-end processor that controls signal processing in the receiver 104, and the processor 106 can be a back-end processor that performs further computations based on the signal processing in the receiver 104. In embodiments, the processor 106 can assign computational tasks to a remote computer (not shown) such that the remote computer performs a portion of the computations and sends the results of the computations to the processor 106.

[0057] The memory 108 can be any type of computer readable memory, including volatile or non-volatile memory devices or a combination thereof. The memory 108 can store information related to the identity of the device 100 and GNSS signals received by the antenna 102. The memory 108 can also store post-processed signals. The memory 108 can also store pseudorange and carrier phase measurements and quality metrics associated with these measurements and measurements from other sensors. The memory 108 can also store computer readable program instructions, mathematical models, and algorithms for signal processing in the receiver 104 and computations in the processor 106. The memory 108 can also store computer readable program instructions for execution by the processor 106 to operate the device 100.

[0058] The local clock 110 provides the time of the local in which the device 100 is located. The local clock 110 can be used to determine the time of arrival of signals for pseudorange measurements and position estimation.

[0059] The input / output device 112 can be used to communicate the results of signal processing and computations to a user or another device. The input / output device 112 can include a user interface that includes a display and input devices to send user commands to the processor 106. The display can be configured to display signal reception status at the device 100, data stored at the memory 108, signal processing status, and computation results, among others. For example, the display can show the results of determining a protection level to the user so that the user can better understand the position of the device 100. The display can include, but is not limited to, a cathode ray tube (CRT), a liquid crystal display (LCD), a light-emitting diode (LED), a gas plasma display, a touch screen, or other image projection device used to display information to a user. The input devices can be any type of computer hardware device for receiving data and control signals from a user. The input devices can include, but are not limited to, a keyboard, a mouse, a scanner, a digital camera, a joystick, a trackball, cursor direction keys, a touch screen monitor, or an audio / video commander, among others. The input / output device 112 can also include a machine interface such as an electrical bus connection or a wireless communication link.

[0060] In embodiments, the mathematical models and algorithms used to determine the protection level of a position estimate can be based on Bayes' theorem.

[0061] Bayes' theorem provides a way to determine the probability density of an event based on known prior information about the event. Bayes' theorem can be expressed using the following equation:

[0062]

[0063] where P(A|B) is the conditional probability of event A given that B is true, P(B|A) is the conditional probability of event B given that A is true, and P(A) and P(B) are the probabilities of A and B, respectively, considered independently.

[0064] When a set of continuous domain measurements are applied to an inference problem, Bayes' formula can be converted to the following equation:

[0065]

[0066] where P(state|data) is the probability density of a particular state given a set of observations, P(data|state) is the probability density of a particular set of observations given a state, P(state) is the prior probability density of a state, and P(data) is the probability density of a set of observations. Information about a state can be inferred from the observables (e.g., pseudoranges or carrier phases), and P(state|data) corresponds to the posterior probability density that needs to be determined in order to compute the protection level of a position estimate. P(data) can be treated as an unknown normalization factor, and can be inferred using the fact that the integral of the posterior probability density P(state|data) over all states is equal to 1. P(data|state) is closely related to the measurement error probability distribution and can be specified by a mathematical model as described below.

[0067] In embodiments, the above equation (2) is applied to GNSS measurements. The state (x) can include, but is not limited to, position, clock bias, instrument bias, or atmospheric parameters. The clock bias can include receiver clock bias and satellite clock bias. The GNSS measurements can include a set of observations such as pseudorange measurements and carrier phase measurements. The measurements can also include a set of quality metrics q (metadata) that indicate information about the quality of the observations, such as the carrier noise density ratio.

[0068] A set of observations z can be split into two components: a system model h(x) that relates the state x to the observables in the set of observations z, and a random measurement error component (residual) r:

[0069] z = h(x) + r; (3)

[0070] The probability density of the residual r can be specified by a function f joint(r|θ, q) defines the probability density of a set of specified residuals r given error model parameters θ and a set of quality metrics q. Returning to equation (2), the posterior probability density P(x|z, q, θ) can be expressed as:

[0071]

[0072] With the denominator P(z|q, θ) as an unknown normalization constant, and r = z - h(x), the posterior probability density P(x|z, q, θ) can be expressed as follows:

[0073] P(x|z, q, θ) ∝ P(z|x, q, θ)P(x) (5)

[0074] P(x|z, q, θ) ∝ f joint (z - h(x)|θ, q)P(x) (6)

[0075] In an implementation, in determining the posterior probability density P(x|z, q, θ), the observations from a single measurement epoch are considered. For a single measurement epoch, the measurements can be considered independent since there is no temporal correlation between the measurements. Thus, the joint probability density f joint (r|θ, q) can be computed as a simple product of the individual observation probability densities f(r|θ, q):

[0076] f joint (z - h(x)|θ, q) = Π i f(z i -h i (x)|θ, q i ) (7)

[0077] Thus, equation (6) can be simplified as follows:

[0078] P(x|z, q, θ) ∝ P(x)Π i f(z i -h i (x)|θ, q i ) (8)

[0079] In this way, the measurement error distribution is simplified to a one-dimensional function f(r|θ, q) of the residuals r, which can be determined using experimental data. The protection level of the position estimate can also be further determined from the posterior probability density P(x|z, q, θ), for example, by integrating the posterior probability density P(x|z, q, θ) over the state x.

[0080] Figure 2 is a flowchart illustrating a method 200 of determining a protection level of a position estimate consistent with some embodiments of the present disclosure. The method can be performed by a processor, such as a processor 102 of a positioning device 100, or a processor 202 of a server 200, or a processor 302 of a client device 300. Figure 1device 100. Referring to FIG. 2, the method 200 includes a step S210 of specifying prior information about the state x. For example, the prior probability density P(x) of the state x can be specified using knowledge about the state x before taking measurements. In embodiments, specifying the prior information can be performed before taking measurements. Figure 2

[0081] The method 200 includes a step S220 of specifying a system model h(x) relating the state x to the measured observable. The observable can include the pseudorange R and the carrier phase Φ.

[0082] In embodiments, a mathematical model of the pseudorange R can be expressed as follows:

[0083] R = p + c(dt r - dt s ) + T + a f STEC + k P,r + k P,s + m P + e P (9)

[0084] where p is the geometric distance, c is the speed of light in vacuum, dt r and dt s are the receiver and satellite clock biases with respect to the GNSS time scale, T is the tropospheric delay, a f STEC is the ionospheric delay (STEC is the slant total electron count integrating the electron density along the signal path, a f is a frequency-dependent conversion factor relating STEC to delay), k P,r and k P,s are the receiver and satellite code instrument delays depending on code and frequency, m P is the influence of multipath interference on the pseudorange, e P is the influence of receiver noise on the pseudorange.

[0085] In embodiments, a mathematical model of the carrier phase Φ can be expressed as follows:

[0086] Φ = p + c(dt r - dt s ) + T - a f STEC + k L,r - k L,s + λn + λw + m L + e L (10)

[0087] In addition to the terms as defined above, k L,r and k L,s are the receiver and satellite code instrument delays, respectively, m L ​This refers to the effect of multipath interference on the carrier phase, where λ is the wavelength of the transmitted GNSS signal, n is the integer ambiguity, and w is the circular polarization winding (wind-up). L This refers to the effect of receiver noise on carrier phase. Carrier phase measurement can be a traditional single-frequency measurement or a combination of carrier phase measurements obtained from signals of different frequencies (e.g., wideband combination).

[0088] In this implementation, data obtained from one or more physical or virtual base stations, correction services, or enhancement systems can be used to calculate corrections for the pseudorange R and the carrier phase Φ. Specifically, corrections from base stations are considered; the base stations can be reference devices whose exact locations are known. This correction removes satellite bias terms and reduces the effects of ionospheric and tropospheric delays, as well as mitigating errors in clocking and orbit. Corrections based on base station data can be expressed as follows:

[0089] R correction =R base -ρ predicted (11)

[0090] Φ correction =Φ base -ρ predicted / λ (12) where R base and Φ base These are the pseudorange and carrier phase observed from the base station, R. correction and Φ correction These are corrections for pseudorange and carrier phase, respectively, ρ predicted It is the predicted geometric distance given the known location of a base station.

[0091] Based on the rover (e.g., Figure 1 The correction of the data obtained by the device 100 is then calculated as follows:

[0092] R corrected =R rover -R correction (13)

[0093] Φ corrected =Φ rover -Φ correction (14)

[0094] The corrected pseudorange observable can be expressed as follows:

[0095] R corrected =ρ+cΔdt r +ΔT+α f ΔSTEC+Δk P,r +m P,rover -m P,base +∈ P (15)

[0096] where, in addition to the terms defined above, Δdt r is the difference in receiver clock bias between the rover and base stations (Δdt r = dt rover - dt base ), ΔT is the difference in tropospheric delay, a f ΔSTEC is the difference in ionospheric delay, Δk P,r is the difference in receiver code instrument delay, m P,rover is the effect of multipath on the pseudorange measured at the rover, m P,base is the effect of multipath on the pseudorange measured at the base station. Both rover and base station receiver noise have been subsumed into ∈ P .

[0097] Similarly, the corrected phase observable can be represented as:

[0098] Φ corrected = p + cΔdt r + AT - a f ΔSTEC + Δk L,r + l(n rover - n base ) + lAw + m L,rover - m L,base + e L (16)

[0099] where, in addition to the terms defined above, Δk L,r is the difference in receiver code instrument delay, m L,rover is the effect of multipath on the carrier phase measured at the rover, m L,base is the effect of multipath on the carrier phase measured at the base station, n is the integer ambiguity, Δw is the difference in circular polarization windup. Both rover and base station receiver noise have been subsumed into ∈ L .

[0100] In an implementation, the portion of the circular polarization windup term due to satellite rotation can be canceled out, and for single epoch analysis, the windup due to rover rotation can be absorbed by the instrument bias term Δk L,r and thus ignored. If only short baselines (<20 km) are used, the ionospheric term a f ΔSTEC can also be ignored. The tropospheric delay is modeled using a dip factor M(E) that depends only on the satellite elevation angle E. Thus, the final observable model can be simplified as follows:

[0101] R corrected = p + cΔdt r + AT zenith M(E) + Δk P,r + m P,rover-m P,base +∈ P (17)

[0102] Φ corrected =ρ+cΔdt r +ΔT zenith M(E)+Δk L,r +λ(n rover -n base )+m L,rover -m L,base +∈ L (18) In addition to the terms defined above, ΔT zenith It is the difference in zenith tropospheric delay.

[0103] Thus, the observables (pseudorange R and carrier phase Φ) are represented as functions of state x. In this implementation, correction of the observables may be optional. In this implementation, specifying the system model may be performed before measurement.

[0104] In the implementation, state x can be defined as the following state vector:

[0105]

[0106] In addition to the terms defined above, p represents the position of the rover. Instrument deviation Δk P,r and Δk L,r Each GNSS signal type includes a bias; however, by definition, the pseudorange bias on a GNSS is taken as zero and is therefore omitted in the state vector. In the implementation, GNSS includes GPS, Galileo (a GNSS created by the European Union), and GLONASS (a GNSS created by Russia); therefore, a typical epoch has 10 states in x (three in p, one for Δdt). r Two pseudorange biases Δk used for Galileo and GLONASS P,r Three specific phase deviations Δk for GNSS L,r And a for ΔT zenith ).

[0107] Method 200 includes step S230 of quantifying a quality metric q associated with the measurement. In implementations, some quality metrics (e.g., carrier-to-noise ratio, the effect of multipath interference, and satellite elevation angle) may be determined by a GNSS device (e.g., Figure 1Some other quality metrics (e.g., more information of GNSS signals) can be obtained by checking the consistency between different satellites. In embodiments, the quantified quality metric can be performed during the operation (e.g., during the measurement). When checking the consistency between different satellites, the unused quantities (e.g., the changes of pseudoranges and carrier phases close to the epoch of interest) can be the focus of the check.

[0108] In embodiments, some measurements that affect the accuracy and consistency of the method are excluded or down-weighted. For example, some signals that are severely affected by multipath interference are identified and excluded. Some non-line-of-sight signals that can cause large errors are also excluded.

[0109] In embodiments, the Random Sample Consensus (RANSAC) technique is used to exclude measurements that adversely affect the accuracy and consistency of the method. For example, a set of delta phase measurements (including signals in both frequency bands Ll and L2) within a one-second interval can be the input of the RANSAC technique. Only the signals that are phase-locked for at least one second can participate in this method, so all other signals are immediately excluded. Within a one-second short interval, the instrument bias Ak L,r , the tropospheric delay AT zenith M(E) and the base station multipath m L,base are extremely small and can be ignored. Then, the change of phase ΔΦ corrected is given by

[0110] ΔΦ corrected = Δρ + c ΔΔdt r + λ(Δn rover - Δn base ) + Δm L,rover + ε L (20)

[0111] where Δρ is the change of range within a one-second interval, ΔΔdt r is the change of residual clock bias, Δn rover and Δn base are the changes of ambiguities (i.e., cycle slips), Δm L,rover is the change of multipath interference contribution.

[0112] Some signals can be severely affected by non-line-of-sight (NLOS) signals and multipath interference, where the multipath interference term Δm L,rover can be large. Some signals can be severely affected by cycle slips, so the cycle slip term λ(Δn rover - Δn base ) can be large. The signals with a multipath interference term Δm L,rover larger than a threshold multipath interference or the cycle slip term λ(Δn rover - Δn base) signals greater than the threshold cycle slip can be identified as outliers and excluded.

[0113] Figure 3 is a flowchart illustrating a method 300 of identifying outliers of GNSS measurements, consistent with some embodiments of the present disclosure. Referring to Figure 3 , the method 300 includes a step S310 of selecting four random satellites and randomly selecting Ll band or L2 band delta phase measurements for each satellite.

[0114] The method 300 includes a step S320 of solving for changes in position and changes in residual clock bias using the four delta phase measurements. This is solvable because the system is a linear system with four measurements and four unknowns.

[0115] The method 300 includes a step S330 of predicting delta phases for all satellites and bands using the solution in step S320. Based on the solution obtained in step S320, delta phases for all satellites and bands can be predicted. The method 300 includes a step S340 of determining differences between predicted signals and actual signals and counting the number of signals with differences less than a threshold difference. Signals with differences between predicted values and actual values less than the threshold difference are considered normal values of measurements.

[0116] The method 300 includes a step S350 of repeating steps 310-340 for a certain number of iterations. The number of iterations can be determined based on a confidence level that ensures effective convergence to a solution. The method 300 includes a step S360 of identifying measurements with the most number of normal values. Step S360 is performed after all iterations in step S350 are completed.

[0117] The method 300 includes a step S370 of determining whether the highest number of normal values is greater than a threshold number of normal values (e.g., a predetermined minimum number of normal values). If the determination result is “yes”, the method 300 performs a step S380 of forming a final output of normal values of measurements. On the other hand, if the determination result is “no”, the method 300 performs a step S390 of excluding all signals as outliers of measurements. In an embodiment, the number of iterations in step S350 can be 100, the threshold difference in step S340 is 0.02 m, and the threshold number of normal values in step S370 is 15.

[0118] In this way, RANSAC technique is used to identify and exclude outliers of measurements, thereby enhancing accuracy and consistency of measurements.

[0119] Referring back to Figure 2In step S230 of method 200, instead of using the RANSAC technique, a window-based technique can be used to obtain the quality metric of the derivation. The window-based quality metric can be specified by looking for agreement in both the delta phase and the delta pseudorange over a small time window, e.g., up to four seconds long, as described below with reference to Figure 4 The detailed description.

[0120] Figure 4 is a flowchart illustrating a method 400 of quantifying a quality metric associated with GNSS measurements, consistent with some embodiments of the present disclosure. Referring to Figure 4 , method 400 includes a step S410 of defining a window around an epoch of interest. The width of the window is configurable. In embodiments, the window width can be a half-width of two seconds. The window can be symmetric and can have a full-width of four seconds, and can contain data from, e.g., five epochs (i.e., a number of samples n samples = 5). Using epochs after the epoch of interest means that when implemented in a real-time system, the algorithm has a latency equal to the window half-width.

[0121] Method 400 includes a step S420 of determining signals that are present at all of the specified epochs (e.g., five epochs) and that have phase-locked continuity for all of the epochs. The total number of signals n sig is the sum of the number of signals in the L1 band n L1 and the number of signals in the L2 band n L2 (n sig = n L1 + n L2 ). To be included, the signals need to be present and have phase-locked continuity for all of the epochs. Signals that do not satisfy these conditions can be excluded.

[0122] Method 400 includes a step S430 of obtaining phase and pseudorange measurement data. The total number of phase measurements is given by n samples x (n L1 + n L2 ). Similarly, the total number of pseudorange measurements is n samples x (n L1 + n L2 ), providing a total number of data points m = 2n samples x (n L1 + n L2 ). In embodiments, instead of phase and pseudorange measurements, Doppler measurements are performed in the window around the epoch of interest, and the obtaining in step S430 is obtaining Doppler measurement data.

[0123] Method 400 includes step S440 of modeling the measurement data by means of a system representing the position and clock offset relative to the first epoch of the window, together with fixed offset terms for phase and pseudorange in each frequency band. Therefore, the total number of free parameters is n = 4 (n samples -1)+4.

[0124] Method 400 includes step S450, for example, solving the system using a weighted linear least squares problem. For example, m data points are placed in a column vector y, a column vector x is defined to hold n free parameters, and an m×n matrix A is defined for the linear model y = Ax. Next, column vectors of weights w are defined, each measuring one weight, with pseudoranges assigned to the weights 1 / σ. PR The phase is given a weight of 1 / σ phase For a given state x, the residual r can be expressed as follows:

[0125] r = y - Ax (21)

[0126] untie The best-fit linear unbiased estimator is:

[0127]

[0128] This problem can be transformed into its standard form:

[0129] A′=diag(w)A (23-1)

[0130] y′=diag(w)y (23-2)

[0131] And for Solve

[0132]

[0133] We will solve equation (22) numerically.

[0134] In the implementation, the loss function on the residuals can be defined as follows:

[0135]

[0136] Where ρ(s) is the loss function.

[0137] For example, the Huber loss function can be introduced as follows:

[0138]

[0139] The threshold parameter δ is set to 2. This has the effect of subtracting the weight of outliers in the measurement. A robust Gauss-Newton algorithm can then be used to find the best-fit solution. Steps S440 and S450 correspond to identifying a consensus solution, e.g., a consensus on the changes in position and clock bias. Deviations from the consensus can be used to characterize the quality of GNSS measurements as described below.

[0140] Method 400 ends with step S460 of computing the root mean square (RMS) of the residuals for each signal to derive a quality metric. This corresponds to determining the deviation from the consensus. Once the solution of equation (22) is found, the residuals can be computed. For each signal, the RMS of the residuals can be computed, and this forms the derived quality metric. Phase and pseudorange metrics can be computed separately. In an implementation, these metrics can be used to exclude signals. As an example, the parameters of method 400 can be set as follows: the maximum allowed pseudorange RMS is 3 m; the maximum allowed phase RMS is 0.05 cycles; and at least 15 signals must pass both checks, otherwise the entire epoch is excluded.

[0141] In this way, the outliers of the measurements are de-weighted, and both the delta pseudorange and delta phase are used to obtain the derived quality metric, without the need for a large number of iterations to converge, resulting in enhanced efficiency and accuracy of the processing.

[0142] Referring back to Figure 2 , method 200 includes step S240 of specifying a non-Gaussian error probability density model f(r|0,q) and fitting the error probability density model parameters 0 using experimental data. In an implementation, specifying a non-Gaussian error probability density model and fitting the model parameters can be performed prior to the measurements (e.g., a priori offline). In an implementation, the non-Gaussian model is a student-t distribution model, and the pseudorange errors are modeled as a student-t distribution, where the probability density is represented as:

[0143]

[0144] where r is the residual, v is the degrees of freedom parameter, s is the scaling parameter defining the width of the core distribution, and G is the gamma function.

[0145] The parameters v and s are not constants, and can depend on the quality metric q. In an implementation, three quality metrics can be used: satellite elevation angle (el), carrier noise density ratio (cno), and time to lock (t). To simplify the fitting process, the basic student-t parameters s and v are first defined in terms of quantities that are valid for any real number (referred to as unconstrained intermediate parameters s' and v'):

[0146] s = log(l + e (σ′) ) (28)

[0147] v = log(l + e (v′) ) (29)

[0148] Using these unconstrained intermediate parameters σ' and v' allows the model fitting process to be treated as an unconstrained optimization problem, thereby avoiding the need to represent the constraint that both σ and v must be positive numbers.

[0149] The unconstrained student-t parameters are then defined as functions of the quality measures:

[0150]

[0151] where X0is a 2 x 1 matrix and X1is a 2 x 2 matrix. This model has a total of six free parameters: two in X0and four in X1. All of these unconstrained parameters are bundled into a vector θ and all of the quality measures are bundled into q, with the probability density of the residuals r i with quality measures q i being written as:

[0152] f pr (r i | θ, q i ) (31)

[0153] The model is identified by finding the maximum likelihood value of the nine parameters (six free parameters plus residuals r i , vector θ and quality measures q i ) given a set of measured observations. Using Bayes' theorem:

[0154] The maximum likelihood model is the model that maximizes P(model | data). Using θ to represent the model described in equation (32), using r to represent the observed residuals, and using q to represent the associated quality measures, equation (32) becomes:

[0155]

[0156] The maximum likelihood value of θ is denoted as and is given by

[0157]

[0158] This is a maximization problem with respect to θ. P(r | q) is a normalization factor and can be ignored. Equation (35) then becomes:

[0159]

[0160] The prior model parameters P(θ) are assumed to be uniform and can therefore be omitted. Equation (36) can therefore be simplified to:

[0161]

[0162] The entire data set fpr The combined probability density of (r|0, q) can be evaluated as the product of all individual probability densities of the residuals:

[0163]

[0164] In an implementation, instead of evaluating the entire data set directly for a numerical optimization, the log-likelihood is used and equation (39) becomes:

[0165]

[0166] Then, the maximum likelihood model can be found via standard numerical optimization In an implementation, separate models can be created for moving and stationary epochs, as the pseudorange error distribution has different characteristics in these two cases.

[0167] In an implementation, the carrier phase residuals are modeled as a student-t distribution, similar to the pseudorange model described above. In the phase residual model, the error distribution can be modeled as a student-t plus a uniform mixture (bounded by the range ±0.5). Since here the observable is a phase, errors greater than 0.5 cycles do not occur.

[0168] The basic probability density of the mixture is as follows:

[0169]

[0170] As with the pseudorange, the basic parameters w, s and v are not constants, but functions of the quality metric. First, the basic parameters are defined according to the unconstrained values:

[0171] s = log(l + e (σ′) ) (42)

[0172] v = log(l + e (v′) ) (43)

[0173]

[0174] The mapping between w' and w allows w' to take any value while constraining w between 0 and 1. Then, the unconstrained parameters are defined according to the quality metric. Here, only satellite elevation (el) information is used:

[0175]

[0176] The remaining unconstrained parameter w' is considered to be independent of the satellite elevation angle. In total, there are five free parameters, two in X0, two in X1, and one for w'. The model fitting process is then similar to the approach used for pseudoranges. Since the motion of the vehicle does not significantly affect the phase error distribution, a single model covers both stationary and moving epochs.

[0177] For performance optimization, one can approximate the uniform student mixture by taking the maximum of the two components of f(r phase (r) of the maximum of the two components of f(r

[0178]

[0179] Similarly to the pseudoranges described above, equation (46) can also be evaluated using the log-likelihood.

[0180] Figure 5 A comparison of the observed residual error distribution (410), a Gaussian distribution fit to the observed residual error distribution (420), and a student-t distribution fit to the observed residual error distribution (430) are shown, consistent with some embodiments of the present disclosure. The observed residual error distribution is a pseudorange residual error distribution obtained using a receiver set in a moving vehicle, where the satellite (Galileo) elevation angle is between 10-30 degrees. As shown, the Gaussian distribution fit (420) is a poor fit to the observed residual error distribution (410), as the tails of the observed residual error distribution (410) are much more severe than the Gaussian distribution fit. On the other hand, the student-t distribution fit (430) is an improved fit to the observed residual error distribution (410), indicating that the student-t distribution is a better choice for modeling the residual error distribution than the Gaussian distribution. Figure 5

[0181] Referring back to Figure 2 , the method 200 includes a step S250 of defining an unnormalized posterior probability density over the state x. In embodiments, defining the unnormalized posterior probability density can be performed during operation. In embodiments, the unnormalized posterior probability density is defined using equation (8) above. In the case where the mathematical form of the error distribution is defined and the error model parameters θ are fitted to a set of experimental data, the likelihood of a particular set of measurements can be computed given the state x:

[0182] P(z | θ, x) = Π i f(r i | θ, q i ) (47)

[0183] = Π i f(z i -h i (x) | θ, q i ​) (48) Here, h i (x) is the state of the given system, the predicted observation z i is the system model. With P(z) as the unknown normalisation constant, the final posterior probability density is given by:

[0184] P(x|θ,z,q)∝P(x)∏ i f(z i -h i (x)|θ,q i ) (49) where P(x) is the prior state probability. In embodiments, P(x) is considered uniform at all states except for ΔΤ zenith which is considered to have a Gaussian distribution centred on zero with a standard deviation of 5 cm.

[0185] The method 200 comprises a step S260 of estimating the state x and calculating the protection level by integrating the posterior probability density over the state x. The state x can be estimated, for example, using the posterior probability density P(x|z,q,θ). The integration can be done numerically, for example, using Markov Chain Monte Carlo (MCMC) or importance sampling methods. In embodiments, estimating the state x and calculating the protection level can be performed during operation.

[0186] In embodiments, the integration of the posterior probability density over the state x is performed based on MCMC but using multiple chains with modified density functions to give accurate estimates of the tail densities over multimodal distributions. Multimodality can be characterised by the presence of multiple distinct modes of distribution. For example, in a first phase, a set of parallel interacting chains are used to address the problem of sampling from a multimodal posterior. Instead of extracting samples directly from the posterior probability density, a modified probability density function is defined such that the σ parameter of the phase probability distribution (equation (42)) can be varied to:

[0187] σ = max(log(l + e (σ′) ), σ min ) (50)

[0188] This has the effect of broadening the modes of the posterior probability density. If σ min is set to a large value (e.g. 0.5 m), the modes merge and the distribution becomes unimodal. For example, a set of 10 MCMC chains are executed in parallel, each with a different σ min value. These chains can be viewed as forming a chain stack, with a single mode chain that is greatly broadened at the top, and un-broadened chains (σ minbottom. Adjacent chains can interact via replica exchange, i.e., the states of two adjacent chains are sometimes exchanged using the common Metropolis criterion. At the end of the sampling, only the states of the bottom (unwiden) chain are used. This bottom chain can accurately sample a multimodal posterior density.

[0189] Figure 6 is a schematic diagram illustrating a multimodal sampling method using parallel interacting chains consistent with some embodiments of the present disclosure. Ten MCMC chains are executed in parallel (for clarity, Figure 6 only three chains are shown), each MCMC chain having a different σ min value. The epoch has a bimodal posterior with modes at -1.5m and 0m. The unmodified chain (σ min = 0) can move between the two modes at -1.5m and 0m by interacting with the above gradually widened chains. At the highest level (σ min = 0.5), the distribution is unimodal. Thus, using multiple chains with modified density functions, the tail densities on a multimodal distribution can be accurately estimated.

[0190] In embodiments, instead of directly drawing a large number of samples from the posterior probability distribution, importance sampling is used. For example, umbrella sampling is used to draw samples biased towards regions of interest (e.g., the tails of the distribution). Figure 7 is a flowchart illustrating a method 700 of performing umbrella sampling consistent with some embodiments of the present disclosure. As Figure 7 shown, the method 700 includes a step S710 of performing a first round of sampling by drawing a first set of samples from an unbiased distribution and estimating a first position of a quantile on the first set of samples. For example, the quantile can be set to a median value (e.g., 0.2). This quantile will be used to constrain the sampling pool in a second round of sampling.

[0191] The method 700 includes a step S720 of performing a second round of sampling by drawing a second set of samples from a first constrained distribution determined based on the first position of the quantile. The first constrained distribution can be set to the tail of the unbiased distribution. The method 700 includes a step S730 of estimating a second position of the quantile on the second set of samples. The second position of the quantile is used to set a second constrained distribution for a third round of sampling.

[0192] The method 700 includes a step S740 of repeating sampling using the constrained distribution set in the previous sampling until an nth round of sampling is complete. For example, six rounds of sampling can be used. For the nth set of samples, the constrained distribution set in the (n-1)th sampling forces the sampling of only the tail of the q nFor example, for quantile value 0.2, only 0.2 fraction of the tail is sampled in the first round, only 0.04 fraction of the tail is sampled in the second round, only 0.008 fraction of the tail is sampled in the third round, and so on.

[0193] The method 700 comprises a step S750 of determining whether the sampling is sufficient for the required protection level calculation. If the determination in step S740 is "no", the method repeats from step S740 and extracts another round of samples. If, on the other hand, the determination in step S740 is "yes", in step S760 all samples are combined with the appropriate weights and the final protection level is determined.

[0194] In an embodiment, a system comprising a receiver (e.g., device 100 in Figure 1 The receiver is set in a vehicle and the set of experimental data comprises data obtained from the vehicle driving on highways for about 25 hours. Measurement residuals are obtained using a real-time kinematic (RTK) method using local (e.g., within 20 km) reference stations and fixing L1 and L2 integer ambiguities. The same data set is used to fit the error model parameters θ and test the protection level calculation. The data comprises measurements from GPS, Galileo and GLONASS. The phase observations are wide-lane, not available on GPS satellites not transmitting the civilian L2 signal. In line with some embodiments of the present disclosure, Figure 8A and Figure 8B Stanford plots showing the limit half-widths and the true errors obtained from experimental data in the along-track direction, Figure 8C and Figure 8D Stanford plots showing the limit half-widths and the true errors obtained from experimental data in the cross-track dimension, Figure 8E and Figure 8F Stanford plots showing the limit half-widths and the true errors obtained from experimental data in the vertical dimension. Figure 8A and Figure 8B The same data is shown, but with different axis limits. Similarly, Figure 8C and Figure 8D The same data is shown, but with different axis limits, Figure 8E and Figure 8F The same data is shown, but with different axis limits. Figures 8A to 8F Protection levels and actual errors of the test data are shown, with the integrity risk set to one event per 25 hours of data. The availability of the cross-track limits is shown in Table 1 below.

[0195] Table 1

[0196]

[0197] Here, the availability only considers epochs for which sufficient GNSS signals are available to compute the bounds.

[0198] The single-epoch position bound algorithm described above forms the core component of the integrity method. The algorithm can also be used to determine the full (multi-modal) posterior distribution, which can be used for other applications, not just for bounds or position estimates. In embodiments, the bound propagation (e.g., single-epoch bounds) is evaluated at a low rate (e.g., once every 10 seconds), and in parallel, the most recent single-epoch bounds and position are propagated forward to the current time using delta phase observations and sensor data. Mathematically, the bounds at the current time are the sum of the bounds at some previous epoch plus the bounds on the change in position since that time. To maintain integrity, there is a risk of sharing misleading information between the absolute bound computation and the propagation computation. Therefore, the propagation is done with a similar level of integrity as the main bound computation. In this way, the system can maintain high availability, low latency, and reduced computational load.

[0199] In embodiments, a prototype propagation scheme is used. The prototype propagation scheme is based on two different approaches: (1) modeling the error distribution on the delta phase using a non-Gaussian distribution (e.g., student-t distribution) and numerically evaluating the bounds on the delta position using MCMC; and (2) modeling the error distribution using a Gaussian over-bound method, which is used in the aviation domain by the Receiver Autonomous Integrity Monitoring (RAIM) algorithm.

[0200] In embodiments, sensor-based propagation can be used with GNSS signals or when GNSS signals are not available. For example, measurements made by motion sensors (e.g., accelerometers, gyroscopes, wheel speed sensors, or any other sensors placed in the vehicle) can be used for the measurements. In sensor-based propagation, an extended Kalman filter can be used to track the changes in position. Prior to a power outage, GNSS signals can be used to estimate the state of the system (e.g., orientation and sensor biases).

[0201] In embodiments, one or more model parameters can be added if multipath interference is severe. The model parameters can be estimated simultaneously with the position. In embodiments, instead of a local reference station, a State Space Representation (SSR) correction service can be used. The error associated with the SSR correction service can be treated as a coordinate shift in the solution, so the bound validity can be restored by enlarging the limits to allow for the possible size of the shift. Alternatively, the SSR error can be explicitly modeled.

[0202] The computer readable storage medium can be a tangible device that can store instructions for execution by a processor. The computer readable storage medium can be, for example but not limited to, an electronic storage device, a magnetic storage device, an optical storage device, an electromagnetic storage device, a semiconductor storage device, or any suitable combination of the foregoing. A non-exhaustive list of more specific examples of the computer readable storage medium includes the following: a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), a static random access memory (SRAM), a portable compact disc read-only memory (CD-ROM), a digital versatile disk (DVD), a memory stick, a floppy disk, a mechanically encoded device such as punch-cards or raised structures in a groove having instructions recorded thereon, and any suitable combination of the foregoing. A computer readable storage medium, as used herein, is not to be construed as being transitory signals per se, such as radio waves or other freely propagating electromagnetic waves.

[0203] Computer readable program instructions described herein can be downloaded to, and executed by, a computer system from a computer readable storage medium or to an external computer or external memory device or electronic storage device external to a computer. The computer readable storage medium can be, for example, a floppy disk, a CD-ROM, an optical disc storage, or any other medium that can be read by a computer system. The computer readable storage medium can also be, for example, a memory device or electronic storage device such as, for example, a semiconductor memory device, a flash memory, such as a solid state memory device, or a hard disk drive.

[0204] The flow diagrams and the block diagrams in the drawings are used to describe various implementations. It will be understood that each block of the flow diagrams and the block diagrams, and combinations of blocks in the flow diagrams and the block diagrams, can be implemented by computer readable program instructions. Such instructions can be executed by one or more processors of a general- purpose computer. The computer readable program instructions can be stored on a computer readable storage medium that can be implemented in any form of a tangible device.

[0205] It will be understood that the described implementations are not mutually exclusive, and elements, components, materials or steps described in connection with one example implementation can be combined with or removed from other implementations as appropriate to achieve the desired design goals.

[0206] Reference to "some implementations", or "some example implementations", means that a particular feature, structure, or characteristic described in connection with the implementation can be included in at least one implementation. The appearances of the phrase "one implementation", "some implementations" or "another implementation" in various places in the specification are not necessarily all referring to the same implementation, nor are they necessarily mutually exclusive or alternative implementations.

[0207] It should be understood that the steps of the example methods set forth herein are not necessarily required to be performed in the order described, and the order of the steps of these methods should be understood as merely an example. For example, depending on the functionality of the method set forth, two described steps can be executed substantially concurrently or the steps can sometimes be executed in the reverse order depending upon the functionality of each step. Additionally, although each of the methods is described as a series of steps, unless explicitly stated, the steps need not be performed in that exact order, nor do each and every step need necessarily be performed. It will also be appreciated that one or more steps could be at least partially implemented by one or more processors or computers.

[0208] As used in the present disclosure, the word "exemplary" is used herein to mean serving as an example, instance, or illustration. Any aspect or design described herein as "exemplary" is not necessarily to be construed as preferred or advantageous over other aspects or designs. Rather, use of the word exemplary is intended to present concepts in a concrete manner.

[0209] As used in the present disclosure, the term "or" encompasses any possible combination unless specifically stated otherwise or clear from context. For example, if a database is described as including A or B, then the database can include A or B or A and B unless specifically stated otherwise or clear from context. As a second example, if a database is described as including A, B, or C, then the database can include A or B or C or A and B or A and C or B and C or A and B and C unless specifically stated otherwise or clear from context.

[0210] Also, unless otherwise indicated herein, or in the appended claims, the singular forms "a", "an", and "the" are used herein not to denote count uniqueness but rather to denote "one or more."

[0211] Unless specifically stated otherwise, and as can be apparent from the context, quantities, values and ranges should be understood as being approximate, unless exact quantities or values are called for.

[0212] Although the elements in the accompanying method claims, if any, are recited in a particular order, unless otherwise dictated by the claim language, the order of elements in the claims is not necessarily intended to imply a specific order in which the elements are to be performed.

[0213] It will be understood that certain features of the disclosure, which are, for clarity, described in the context of separate embodiments, can also be provided in combination in a single embodiment. Conversely, various features of the disclosure, which are, for brevity, described in the context of a single embodiment, can also be provided separately or in any suitable subcombination or as suitable in any other described embodiment of the disclosure. The disclosure contemplates that in some embodiments, the features of the disclosure can be provided in combination while in other embodiments, the features may be provided separately or in any suitable subcombination.

[0214] It will also be understood that, without departing from the scope of the appended claims, those skilled in the art may make various modifications, substitutions, and variations to the details, materials, and arrangements of the components described and shown to illustrate the nature of the described embodiments. Therefore, the appended claims cover all such substitutions, modifications, and variations falling within the terms of the claims.

Claims

1. A method of determining a protection level for a position estimate using only a single epoch of Global Navigation Satellite System, GNSS, measurements, the method comprising the steps of: pre-specifying a prior probability density P(x) of a state x; pre-specifying a system model h(x) relating the state x to an observable z of the measurements; quantifying a quality metric q associated with the measurements during operation; For each observable z i Pre-specify the non-Gaussian residual error probability density model f(r) i |θ,q i ), and fit the model parameters θ determined offline prior, where r i It is with each of the observables z i The associated residual error, and r = zh(x); defining a posterior probability density P(x|z,q,0) during operation, wherein and i is a natural number that is an index for the observable z; estimating the state x during operation; and computing the protection level during operation by integrating the posterior probability density P(x|z,q,0) over the state x.

2. The method of claim 1, wherein, The step of quantifying the quality metric q associated with the measurements comprises using a Random Sample Consensus, RANSAC, technique to exclude one or more outliers of the measurements.

3. The method of claim 1 or 2, wherein, The non-Gaussian residual error probability density model is a student-t distribution function, and is given for each observable z i specifying the residual error probability density model f(r i |θ,q i ) further comprises: The model of pseudorange is identified using the student-t distribution function as: where r is the residual error, v is a degree of freedom parameter, s is a scaling parameter, and G is a gamma function.

4. The method of claim 1 or 2, wherein, The non-Gaussian residual error probability density model for each observable z is a student-t distribution function, and is specified as f(r i |θ,q i ) further comprises: The model of carrier phase is identified using the student-t distribution function as: where r is the residual error, v is a degree of freedom parameter, s is a scaling parameter, w is a weight parameter, and G is a gamma function.

5. The method of claim 3, wherein, The step of fitting the model parameters 0 using a set of experimental data comprises: defining the degree of freedom parameter and the scaling parameter using an unconstrained student-t parameter; defining the unconstrained student-t parameter as a function of the quality metric; and finding a maximum likelihood value of the unconstrained student-t parameter.

6. The method of claim 4, wherein, The step of fitting the model parameters 0 using a set of experimental data comprises: defining the degree of freedom parameter, the scaling parameter, and the weight parameter using an unconstrained student-t parameter; defining the unconstrained student-t parameter as a function of the quality metric; and finding a maximum likelihood value of the unconstrained student-t parameter.

7. The method of claim 1 or 2, wherein, The step of integrating the posterior probability density over the state comprises: applying Markov Chain Monte Carlo, MCMC, numerical integration using a plurality of interacting chains with modified density functions.

8. The method of claim 1 or 2, wherein, The step of integrating the posterior probability density over the state comprises: performing a first round of sampling by extracting a first set of samples from an unbiased distribution and estimating a first position of a quantile over the first set of samples; performing a second round of sampling by extracting a second set of samples from a first constrained distribution determined based on the first position of the quantile; estimating a second position of the quantile over the second set of samples, wherein the second position of the quantile is used to set a second constrained distribution for a third round of sampling; repeating the sampling using the constrained distribution set in a previous round of sampling until an nth round of sampling is completed, wherein n is a natural number greater than 2; determining whether the samples of n rounds of sampling are sufficient to compute the protection level; and combining the n rounds of sampled samples with appropriate weights when the n rounds of samples are sufficient to compute the protection level.

9. The method of claim 1 or 2, further comprising the step of: propagating the most recent single epoch bound forward to the current time using the delta phase error profile or data obtained by one or more sensors.

10. The method of claim 9, wherein, The step of propagating further comprises: modeling the error profile on delta phase using a non-Gaussian distribution function and numerically evaluating the bound on delta position using Markov Chain Monte Carlo techniques; or modeling the error profile on delta position using a Gaussian overshoot method.

11. The method of claim 1 or 2, wherein, The state x is defined as: where p is the position of the rover, c is the speed of light in vacuum, At r is the difference in receiver clock bias between the rover and base station, Ak P,r is the difference in receiver code instrument delay on the pseudorange, Ak L,r is the difference in receiver code instrument delay on the carrier phase, and AT zenith is the difference in zenith tropospheric delay.

12. A device configured to determine a protection level of a position estimate using only a single epoch of GNSS measurements, the device comprising: a receiver configured to receive GNSS signals and process the received signals to make measurements; and a processor configured to: specify a prior probability density P(x) of a state x; specify a system model h(x) relating the state x to the measured observable z; quantify a quality metric q associated with the measurements; for each observable z i specifying a non-Gaussian residual error probability density model f(r i | θ, q i ) and fitting the model parameters θ using a set of experimental data, where r i is the residual error associated with said each observable z i and r = z - h(x); define a posterior probability density P(x|z, q, 0) where and i is a natural number that is an index for the observable z; estimate the state x; and compute the protection level by integrating the posterior probability density P(x|z, q, 0) over the state x.

13. The apparatus of claim 12, wherein, The processor is further configured to: communicate with external sensors to obtain position change data tracked by the sensors to determine bound propagation.

14. A non-transitory computer readable medium having instructions stored thereon that, when executed by a processor, perform a method of determining a protection level of a position estimate using only a single epoch of GNSS measurements, the method comprising the steps of: specifying a prior probability density P(x) of a state x; specifying a system model h(x) relating the state x to the measured observable z; quantifying a quality metric q associated with the measurements; for each observable z i specifying a non-Gaussian residual error probability density model f(r i | θ, q i ), where r i is a residual error associated with said each observable z i , and r = z - h(x), model parameters θ being determined a priori offline; defining a posterior probability density P(x|z, q, 0) where and i is a natural number that is an index for the observable z; and estimating the state x and computing the protection level by integrating the posterior probability density P(x|z, q, 0) over the state x.

Images