Integration for GNSS measurement processing

Abstract

The invention relates to integration for GNSS measurement processing. A GNSS receiver and a method of processing GNSS measurements are disclosed. The method comprises defining (220) a state vector containing a state variable. The method further comprises obtaining (230) a posterior probability density of the state vector. State information is inferred based on the posterior probability density. The inference involves integrating a posterior probability density. The integration includes dividing (234) a domain of the posterior probability density into a plurality of slices along a dimension associated with one of the state variables; individually integrating (236) the respective slices to produce respective integrated results for the respective slices; and combining (238) the integrated results of the slices.

Description

Technical Field

[0001] The present invention relates to a global navigation satellite system (GNSS), and more particularly to a method and apparatus for processing GNSS measurements to infer state information by integrating a posterior probability density function. Background Art

[0002] GNSS positioning technology is well known in the art. Existing GNSS include the Global Positioning System (GPS), Galileo, GLONASS and the BeiDou Satellite Navigation System (BDS), also referred to herein as "BeiDou". Each GNSS includes a satellite constellation, also referred to in the art as a "space vehicle" (SV), which orbits the earth. Typically, each SV sends multiple satellite signals. These signals are received by a GNSS receiver, for which state information, such as a navigation solution, needs to be calculated. A GNSS receiver can use these signals to generate multiple ranging and / or Doppler measurements to derive information about the distance (and / or the rate of change of these distances) between the receiver and the corresponding satellite. When a sufficient number of measurements can be made, the position and / or velocity of the receiver can be calculated.

[0003] Various influences can introduce errors into the measurements. Errors in the measurements, in turn, can lead to errors in the resulting state estimates. Sources of error include, for example, multipath interference. This can be a particular problem for road vehicle navigation systems. In signal environments often encountered by road vehicles, multipath interference can be caused by other vehicles, buildings or vegetation reflecting or blocking satellite signals.

[0004] Other sources of error include changes in atmospheric delays as satellite signals pass through the troposphere or ionosphere. In addition to this "innocent" error source, there may also be malicious interference, such as jamming or spoofing signals. Jamming signals are sent with the intention of deliberately preventing a GNSS receiver from successfully receiving satellite signals. Spoofing involves sending signals that mimic real satellite signals with the intention of causing a GNSS receiver to mistakenly acquire and track the spoofed signals. An attacker can use this to intentionally introduce false information into the multilateration positioning calculation, resulting in a false position or velocity estimate.

[0005] Obviously, despite all of these sources of error, it is desirable to compute the state estimate as accurately as possible. However, it is inevitable that they will sometimes result in inaccurate state estimates. It is desirable to know whether the state estimate is accurate, and to quantify the expected accuracy. For example, in autonomous driving applications, it is important to be able to determine when the uncertainty in the vehicle's position is so large that the position estimate cannot be safely used to control the navigation of the vehicle.

[0006] In this context, the concept of a "warning limit" or "alarm limit" is sometimes used. This may be defined as the maximum permissible error (eg position error) that may be permitted without triggering a dangerous situation.

[0007] The concept of “integrity risk” is also used. This can be defined as the probability that the actual position error exceeds the warning limit in the absence of any warning signs at a given moment. Therefore, integrity risk describes the trust or confidence level of the GNSS position estimate.

[0008] The uncertainty associated with a position estimate can be quantified in terms of a “protection level”. This is defined as a statistical error limit, calculated such that the probability of a position error exceeding the limit is less than or equal to a target integrity risk. The protection level can be provided in terms of a radius, in which case the protection level defines a circle within which the probability of the vehicle being located is (quantifiably) high. The protection level can be used in conjunction with a warning limit to warn the user of a navigation application, or to trigger other mitigation measures in an automated system. Whenever the protection level exceeds a warning limit, the system issues a warning to the user or takes the necessary mitigation measures (as appropriate). However, there is always a small risk that the position estimate deviates from the true position by more than the reported protection level. One of the goals of an integrity monitoring system is to minimize this risk by estimating the position or error limit (e.g., protection level) as reliably as possible.

[0009] EP 3 859 397 A1 discloses a method for determining a protection level for a position estimate using a single epoch of GNSS measurement results. The method comprises: specifying a prior probability density P(x) of a state x; specifying a system model h(x) relating the state x to an observable z of the measurement result; quantifying a quality indicator q associated with the measurement result; specifying a non-Gaussian residual error probability density model f(r|θ,q) and fitting the model parameters θ using experimental data; and defining a posterior probability density P(x|z,q,θ); estimating the state x; and calculating the protection level by integrating the posterior probability density P(x|z,q,θ) over the state x. EP 3 859 397 A1 proposes a specific example of a Student's T distribution as a non-Gaussian residual error probability density model. Summary of the invention

[0010] As described above, it is desirable to ensure the integrity of the navigation solution by calculating reliable error limits that are tight enough to provide high availability of the system. More generally, it is desirable to infer state information as accurately as possible from GNSS measurements. A key aspect of achieving these goals is to model the probability density function of the error as accurately as possible in the set of GNSS measurements used to infer state information. Errors in modeling the probability density function may cause errors in the resulting inferred state information. Errors in the state limits themselves may manifest as too high or too low confidence in the state information. Both types of errors may have their own problems. For example, if the estimate provided for the protection level is too low, the risk of the true error exceeding the alarm limit will increase without the system realizing it. On the contrary, if the estimate provided for the protection level is too high, the system may issue too many false alarms.

[0011] It is also desirable to perform the integration of the posterior probability density as accurately and efficiently as possible. In particular, for safety-critical applications, it is desirable to accurately quantify tail risk. "Tail risk" refers to extreme events that occur with low probability - that is, events that lie in the tail of the posterior probability density. It may be desirable to integrate the tail of the density function in order to quantify the total probability of "extreme" events occurring. In some examples, it may be desirable to integrate the tail (defined by one or more specific thresholds) in order to determine the total probability of events occurring outside the thresholds. In other examples, it may be desirable to integrate the tail in order to find one or more thresholds that are consistent with a predefined total probability of extreme events. For example, it may be desirable to answer questions such as "If we accept 10 extreme events to occur per hour, what is the total probability of extreme events?" -7 errors (i.e., 10 -7 What level of protection should we choose to achieve this goal?

[0012] The difficulty arises because typical methods of numerically integrating the probability density function describing GNSS measurements tend to be more accurate in the region of the density function where the probability is highest. The inventors have recognised that this approach may not give accurate results when integrating over the tails of the probability density.

[0013] In this regard, particular challenges arise when GNSS measurements include carrier phase measurements. Carrier phase measurements have the potential to significantly improve the accuracy of position estimates. However, due to the integer ambiguities inherent in individual carrier phase measurements, the posterior probability density (e.g., over position) is multimodal. Each mode of the distribution corresponds to a possible combination of integer values ​​for resolving (or "fixing") the ambiguity. Integrating over the tails of these multimodal probability distributions presents special challenges to the integration methods that have been used for this purpose in the past.

[0014] The inventors have recognized that it may be advantageous to explicitly identify the modes. Information about the modes may be exploited to improve the efficiency of the integration (which may involve numerical integration), for example by enabling more targeted random sampling of the probability density. This approach may be of particular benefit in the case of non-Gaussian distributions (which allows for more accurate modeling of the measurement errors). Existing algorithms for handling integer ambiguities in carrier phase measurements tend to assume a Gaussian distribution of the errors. Such algorithms allow for computational efficiency at the expense of less accurate modeling of the error distribution.

[0015] The present invention discloses a GNSS receiver and a method for processing GNSS measurement results. The method includes defining a state vector, the state vector including state variables. The method also includes obtaining a posterior probability density of the state vector. State information is inferred based on the posterior probability density. The inference involves integrating the posterior probability density. The integration includes dividing the domain of the posterior probability density into multiple slices along a dimension associated with one of the state variables; integrating each slice separately to produce a corresponding integration result for each slice; and combining the integration results of the slices.

[0016] According to one aspect of the present disclosure, a method for processing a plurality of GNSS measurement results is provided, the method comprising the following steps:

[0017] obtaining the plurality of GNSS measurement results;

[0018] defining a state vector, wherein the state vector includes state variables;

[0019] obtaining a posterior probability density of the state vector, wherein the posterior probability density is based on one or more residual error models, the one or more residual error models describing a probability distribution of errors for respective ones of the GNSS measurements, the one or more residual error models comprising at least one non-Gaussian model, wherein the posterior probability density has a plurality of dimensions, each dimension corresponding to one of the state variables; and

[0020] Inferring state information based on the posterior probability density, wherein the step of inferring comprises integrating the posterior probability density,

[0021] The step of integrating includes, for one of the state variables:

[0022] partitioning the domain of the posterior probability density into a plurality of slices along a dimension associated with the one state variable;

[0023] Integrate each slice of the plurality of slices individually to generate a corresponding integration result for each slice; and

[0024] The integration results of the slices are combined.

[0025] The one state variable on which the domain is divided into slices may be referred to as a “state of interest.” The state variable may include one or more position variables. The state of interest may be one of these one or more position variables.

[0026] The state variables may include two or more position variables, and the method may include integrating the posterior probability density of each position variable, wherein the integration includes, for each of the position variables (individually): dividing the domain of the posterior probability density into multiple slices along a dimension associated with the position variable; and integrating each slice individually.

[0027] In this way, the integral can be divided into slices along different dimensions in order to infer state information about the variables at each location.

[0028] In some examples, the position variables may include any one, any two, or all three of the following: along-track position, cross-track position, and vertical position. The three position variables form an orthogonal system, where the along-track position is aligned with the current direction of travel. When the GNSS receiver is stationary, the along-track position may be aligned with the last estimated direction of travel.

[0029] The number of slices may be at least 20, 50, 100, or 150. The number of slices may be less than 1000, 750, 500, 300, or 250. All combinations of these upper and lower limits are disclosed herein. In some examples, the number of slices is approximately 200 (optionally exactly 200).

[0030] Each slice may have a finite width in the dimension of interest (the dimension associated with the one state variable). That is, each slice may be associated with a bounded integration domain defined in the dimension of interest. Each slice may be relatively narrow compared to the overall integration domain spanned by the set of slices.

[0031] Each slice may have a width in the dimension of interest that is selected based on the width of a peak in the posterior probability density. The peak of the width used may be the peak associated with the maximum likelihood state. The slice width may be set proportionally to the peak width. In particular, the slice width may be set equal to a predetermined fraction or proportion (less than or equal to 100%) of the peak width, such as half the peak width. The peak width may be defined as the full width at half maximum (FWHM).

[0032] Individual slices may be unbounded or loosely bounded in dimensions other than the dimension of interest (i.e., dimensions associated with state variables other than the one state variable). Here, "unbounded" means that for the "other" dimensions, the integration is over the entire space. "Loosely" bounded means that the bounds on the "other" dimensions are much wider than the bounds on the dimension of interest. Individual slices may extend across the entire integration domain over which state information is to be inferred in the "other" dimensions. All slices may have the same (wide) bounds on dimensions other than the dimension of interest, or may be similarly unbounded.

[0033] The set of slices may span a continuous domain in the dimension of interest. The width of the domain may be selected based on the width of the peak in the second a posteriori probability density. The second a posteriori probability density may be calculated based only on pseudorange measurements and not on any carrier phase measurements. The peak in the second a posteriori probability density is typically wider than the peak in the a posteriori probability density based on all measurements (i.e., also including carrier phase measurements).

[0034] The integration may include (or consist of) integrating a plurality of disturbance states, wherein the disturbance states include all states in the state vector except a state of interest (eg, a position variable of interest).

[0035] The plurality of GNSS measurements may include a plurality of carrier phase measurements.

[0036] Integration may include numerical integration, but is not limited thereto.

[0037] The plurality of GNSS measurements may include pseudorange measurements. The a posteriori probability density may be constrained by the pseudorange measurements.

[0038] The at least one non-Gaussian model may include any one or any combination of two or more of the following:

[0039] Student’s T distribution, to model pseudorange measurements;

[0040] Student's T distribution, which models carrier range measurements (specifically, carrier range measurements whose integer ambiguities have been fixed); and

[0041] Periodic probability density,models range measurements (especially floating point range measurements,,for which integer ambiguities remain unfixed).

[0042] State variables may also include speed variables and / or time variables.

[0043] A single residual error model may describe the probability distribution of one or more GNSS measurements. That is, the error in one GNSS measurement may be modeled as a univariate distribution, or the errors in two or more GNSS measurements may be jointly modeled as a multivariate distribution. Each GNSS measurement is modeled by one of the (one or more) residual error models. Thus, the number of residual error models is less than or equal to the number of GNSS measurements. In one extreme case, each GNSS measurement may be modeled by a separate (univariate) residual error model. In another extreme case, all GNSS measurements may be modeled together by a single (multivariate) residual error model.

[0044] Of particular interest for this disclosure are one or more residual error models that describe the probability distribution of the errors in those carrier phase measurements.

[0045] The residual error model may be pre-generated for use in the inference method. The residual error model may be generated from empirical data including GNSS measurements, quality indicators associated with those GNSS measurements, and residual errors associated with those GNSS measurements. In general, such data may be referred to as empirical data or training data.

[0046] Obtaining the posterior probability density may include calculating the posterior probability density. Calculating the posterior probability density may include multiplying the probability distribution of the error by the prior probability density of the state vector.

[0047] The posterior probability density may be calculated based on the GNSS measurements from a single measurement epoch.

[0048] The state variable may include a carrier phase ambiguity. The carrier phase ambiguity may include a floating point carrier phase ambiguity value.

[0049] The inferred state information may include at least one of: a state estimate; and an error bound on the state estimate.

[0050] The state estimate may be a position estimate. The error bound may be an error bound for the position estimate.

[0051] The error bounds on the state vector define probability bounds on the range of the individual state variables. In some cases, bounds can be calculated for state variables (e.g., position variables) in addition to estimating the state variables themselves. In other cases, bounds can be calculated without estimating the state variables themselves.

[0052] Calculating estimates of one or more state variables of the state vector may include calculating maximum likelihood estimates of the state variables. Similarly, calculating error bounds may include finding the posterior probability density of the state variables to calculate bounds on these variables.

[0053] The error bound of the position estimate may include, or consist of, a protection level.

[0054] The step of integrating each slice individually may include, for each slice, performing at least one of the following: importance sampling; MCMC method; and approximating the posterior probability density using a mathematical model that can perform analytical integration.

[0055] The first two alternatives are examples of numerical integration techniques. Various numerical techniques can be used to improve the efficiency of numerical integration, which are based on the knowledge of the identified set of modes. The latter alternative does not use numerical integration. The mathematical model used to approximate the posterior probability density can be a parametric model, where the parameters of the model are based on the identified set of modes. Suitable parametric models include (but are not limited to) mixed models, such as Gaussian mixture models (GMM).

[0056] Optionally, the step of integrating comprises, for each slice: performing a search in the slice to identify a set of modes of the posterior probability density, and wherein the integrating is based on the identified set of modes.

[0057] The method may further include outputting a list of modalities. Optionally, the method may include outputting each modality and an associated probability.

[0058] The step of integrating each slice individually may include, for each slice, performing at least one of the following: importance sampling based on the identified modal set; MCMC method based on the identified modal set; and approximating the posterior probability density using a mathematical model that can be analytically integrated, wherein the mathematical model is based on the identified modal set.

[0059] The search may be a systematic search. Here, "systematic" means that the search is deterministic, rather than random / random. It also means that the search is complete (in other words, exhaustive) within certain specified limits. Thus, the search is guaranteed to find all modalities that meet the specified criteria.

[0060] In some examples, the search includes a recursive depth-first search of the set of carrier phase ambiguities, meaning that a search for one ambiguity triggers a search for the next ambiguity until the complete set of candidate ambiguities is solved.

[0061] In some examples, the search includes the steps of transforming the posterior probability density into a mixture model comprising multiple mixture components, wherein each mixture component is a multivariate distribution; and using the mixture model to identify the set of modalities, wherein each modality is associated with a corresponding multivariate distribution in the multivariate distributions.

[0062] The inventors have recognized that, by transforming the posterior probability density into a hybrid model and using the hybrid model to identify the modal set, rather than directly searching the original form of the posterior probability density, the modality can be efficiently identified. In some examples, according to the original form of the posterior probability density, the transformation may include approximating the posterior probability density by one or more different functions. However, approximating the posterior probability density in this way has been found to be advantageous because it can allow the modality to be found more easily without affecting the overall quality of the inference. This is because state information can still be inferred based on the posterior probability density in its original form (i.e., rather than in its approximated form). In other words, the approximated form is only used to help find the modality. The modality is then used together with the original form of the posterior probability density for inference (including, for example, any numerical integration).

[0063] The mixture model may include a sum of the mixture components. The step of identifying the set of modes using the mixture model may include, for each of the plurality of mixture components, identifying a mode for each mixture component.

[0064] The multivariate distributions that form each mixture component may include the product of univariate distributions.

[0065] Each mixture component may comprise (or consist of) a linear distribution. In particular, it may comprise (or consist of) a product of linear distributions. Here, the term "linear" distribution is used to denote a distribution that is not a parcel distribution. In one non-limiting example, each mixture component comprises a product of Student's T distributions (e.g., one Student's T distribution per carrier phase measurement).

[0066] The modality set comprises at least two modalities, and optionally comprises at least 10, 100 or 1000 modalities.

[0067] The a posteriori probability density can be obtained in raw form, which includes the periodic distribution of the residual errors in the individual carrier phase measurements.

[0068] The step of transforming the posterior probability density optionally includes: defining a package distribution for each carrier phase measurement result; and transforming the package distribution into the mixed model, wherein each mode is associated with an integer set a, and each integer is associated with a corresponding carrier phase measurement result in the carrier phase measurement results, wherein each integer indicates the number of cycles in the corresponding package distribution.

[0069] In some examples, the posterior probability density can be obtained in a form that already contains these parcel distributions. In other examples, a parcel distribution is defined that approximates the obtained posterior probability density.

[0070] Inferring the state information may include performing an integral. The transformation may also include changing the integration domain of the integral. In particular, the integration domain may be changed such that the integration is performed over a set of real numbers (rather than over a period of the carrier phase state).

[0071] The step of identifying the set of modes may include: defining a floating-point cost function based on the mixture model by relaxing the constraint that each value a indicating the number of cycles in the corresponding carrier phase measurement result is an integer; performing a first local search of the floating-point cost function to find a floating-point valued state vector associated with a local minimum of the cost function; approximating the floating-point cost function in the region of the local minimum by a multivariate Gaussian distribution; and identifying the set of modes based on the approximation using a multivariate Gaussian distribution.

[0072] The first local search may include solving a nonlinear optimization problem. The first local search may involve using an algorithm based on Gauss-Newton (e.g., Levenburg-Marquardt algorithm, or a mixture between Newton's method and the steepest descent method, as described in Section 2.2 (Formula 2.11) of K. Madsen, HB Nielsen, O. Tingleff, "METHODS FOR NON-LINEAR LEAST SQUARESPROBLEMS" in 2004). The first local search may be initialized by performing Gaussian approximation on each univariate distribution and finding the minimum of the resulting approximate cost function. The minimum of the approximate cost function may be found analytically.

[0073] The first local search can be designed to find the global minimum of the cost function. However, it may not be guaranteed to always find the global minimum. Sometimes it may find another local minimum close to the global minimum.

[0074] As described above, each mode can be associated with a set of integers a, each integer being associated with a corresponding carrier phase measurement. Each integer can indicate the number of cycles in the corresponding parcel distribution. In the approximation of the multivariate Gaussian distribution, the integer a can be replaced by a real number a*.

[0075] The approximation may be expressed in the form of a negative log-likelihood using a multivariate Gaussian distribution, wherein the multivariate Gaussian distribution is modeled as a quadratic function.

[0076] The step of identifying the set of modes based on an approximated multivariate Gaussian distribution may include characterizing the approximated multivariate Gaussian distribution in the region of the local minimum by a matrix, the matrix being one of: a covariance matrix and a Hessian matrix.

[0077] The step of identifying the modal set based on approximation using a multivariate Gaussian distribution optionally also includes: defining a weighted square norm based on the matrix based on the state vector associated with the local minimum and the matrix; and identifying a candidate integer set a whose distance from the local minimum is less than a threshold value based on the weighted square norm.

[0078] The step of identifying the candidate integer set a may use any suitable algorithm, including but not limited to the Least Squares Ambiguity Decorrelation Adjustment (LAMBDA) algorithm.

[0079] The step of identifying a set of candidate integers may include identifying a maximum likelihood set of candidate integers; finding a first weighted square norm associated with the maximum likelihood set of candidate integers; and identifying other sets of candidate integers having weighted square norms within a predetermined search threshold of the first weighted square norm.

[0080] Identifying the maximum likelihood candidate set may include running the LAMBDA algorithm for a first time.Identifying the other candidate sets may include running the LAMBDA algorithm for a second time.

[0081] The method may also include, for each identified candidate integer set: defining another cost function based on the floating point cost function, wherein, in the another cost function, each value a indicating the number of cycles in the corresponding carrier phase measurement result is constrained to be a corresponding integer in the candidate integer set; and performing a second local search on the another cost function to find a state vector associated with a local minimum of the cost function, the cost function being associated with the candidate integer set.

[0082] Each second local search may include solving a nonlinear optimization problem. The second local search may involve using an algorithm based on Gauss-Newton (e.g., the Levenburg-Marquardt algorithm as described above, or a mixture between Newton's method and the steepest descent method). The second local search may be designed to find a local minimum of the cost function. The local minimum may be used to characterize the corresponding mode. The local minimum of the cost function may correspond to a local maximum of the posterior probability density.

[0083] The second local search can also produce a matrix that characterizes the curvature of the posterior probability density in the region of the local minimum. The matrix can be one of the following: a covariance matrix and a Hessian matrix.

[0084] Before transforming the obtained a posteriori probability density to the hybrid model, the obtained a posteriori probability density can be defined as the product of: multiple first residual error models, which describe the probability distribution of errors in corresponding pseudorange measurement results; multiple second residual error models, which describe the probability distribution of errors in corresponding first carrier phase measurement results; and multiple third residual error models, which describe the probability distribution of errors in corresponding second carrier phase measurement results, wherein the first carrier phase measurement results include one carrier phase measurement result for each visible GNSS constellation and for each GNSS frequency band, the second carrier phase measurement results include the remaining carrier phase measurement results, the second residual error model does not include a phase deviation term, and the third residual error model includes a phase deviation term.

[0085] According to this form of the a posteriori probability density, one carrier phase measurement from each frequency band of each GNSS constellation is assigned to the first plurality and the other carrier phase measurements from that frequency band of that GNSS constellation are assigned to the second plurality.

[0086] The state vector may include a phase deviation for each GNSS constellation and each GNSS frequency band.

[0087] The system search may include performing a recursive search with certain recursion levels on each carrier phase measurement result, each recursion level optionally including: defining one or more candidate integer ambiguity values ​​associated with the carrier phase measurement result of that level; and testing each candidate integer ambiguity value.

[0088] Defining the one or more candidate integer ambiguity values ​​may include identifying a range of integer values ​​having a posterior probability density above a predetermined threshold.

[0089] The state vector may include a carrier phase ambiguity, and testing each candidate integer ambiguity value may include: fixing the candidate integer ambiguity value; and performing a local search on remaining carrier phase ambiguities and other state variables of the state vector.

[0090] Local searches may include solving nonlinear optimization problems.

[0091] Fixing the candidate integer ambiguity values ​​optionally includes replacing the periodic probability density model of the carrier phase measurements with a unimodal probability density model.

[0092] Testing each candidate integer ambiguity value optionally includes: finding a probability density associated with the candidate integer ambiguity value; and comparing the probability density to a threshold value.

[0093] The method may include, in response to the probability density exceeding a threshold, continuing recursively to a next level.

[0094] The method may include, in response to the probability density not exceeding a threshold, terminating the recursion at the current level and returning to a previous level.

[0095] The method optionally includes, at a final level of the recursion, identifying modes of the posterior probability density for respective candidate integer ambiguity values, wherein respective modes are associated with a combination of candidate integer ambiguity values ​​from all levels.

[0096] A computer program is also provided, the computer program comprising computer program code, the computer program code being configured to cause one or more processors to perform all the steps of the method as described above when the computer program is run on the one or more processors.

[0097] The one or more processors may comprise or consist of one or more processors of a GNSS receiver.The computer program may be stored on a computer-readable storage medium (optionally non-transitory).

[0098] A GNSS receiver is also provided, the GNSS receiver comprising:

[0099] a signal processing unit configured to generate a plurality of GNSS measurements; and

[0100] at least one processor, the at least one processor being configured to:

[0101] obtaining the plurality of GNSS measurement results;

[0102] defining a state vector, wherein the state vector includes state variables;

[0103] obtaining a posterior probability density of the state vector, wherein the posterior probability density is based on one or more residual error models, the one or more residual error models describing a probability distribution of errors for respective ones of the GNSS measurements, the one or more residual error models comprising at least one non-Gaussian model, wherein the posterior probability density has a plurality of dimensions, each dimension corresponding to one of the state variables; and

[0104] Inferring state information based on the posterior probability density, wherein the step of inferring comprises integrating the posterior probability density,

[0105] The step of integrating includes, for one of the state variables:

[0106] partitioning the domain of the posterior probability density into a plurality of slices along a dimension associated with the one state variable;

[0107] Integrate each slice of the plurality of slices individually to generate a corresponding integration result for each slice; and

[0108] The integration results of the slices are combined.

[0109] The GNSS receiver may further include an RF front end. The RF front end may be configured to receive GNSS signals via an antenna. The signal processing unit may be configured to perform GNSS measurements (range measurements and optional pseudorange measurements) on the GNSS signals received by the RF front end. BRIEF DESCRIPTION OF THE DRAWINGS

[0110] The present invention will now be described by way of example with reference to the accompanying drawings, in which:

[0111] Figure 1 is a schematic block diagram of a GNSS receiver according to an example;

[0112] Figure 2A is based on an example, exemplified by Figure 1 A flowchart of a method performed by a GNSS receiver;

[0113] Figure 2B This example is applicable to Figure 2A A flowchart of an example inference method;

[0114] Figure 2C is a flow chart illustrating a method of integrating a slice according to an example;

[0115] Figure 3 is a flowchart illustrating one way to transform the posterior probability density into a mixture model, which is applicable to Figure 2C Methods;

[0116] Figure 4 is a flowchart illustrating one way of identifying a modality set using a mixture model, which is applicable to Figure 2C Methods;

[0117] Figure 5 is exemplified in Figure 4 A flowchart of a method for identifying modes based on approximation using a Gaussian distribution; and

[0118] Figure 6 is a flowchart illustrating one way to identify a set of candidate integers, which is applicable to Figure 5method.

[0119] It should be noted that these drawings are schematic diagrams only and are not drawn to scale. DETAILED DESCRIPTION

[0120] Reference will now be made in detail to examples according to the present disclosure illustrated in the accompanying drawings. The described examples should not be construed as being limited to the description given in this section. Other examples may have different forms.

[0121] In order to understand the following description, it is assumed that the reader is familiar with the general principles described in EP 3 859 397 A1 (which is incorporated herein by reference).

[0122] Figure 1 1 is a schematic block diagram of an apparatus according to an example. The apparatus includes a GNSS antenna 101 and a GNSS receiver 100. The GNSS antenna 101 is configured to receive GNSS signals. It can be configured to receive GNSS signals from a single GNSS constellation (e.g., GPS), or it can be configured to receive GNSS signals from multiple constellations (e.g., GPS, Galileo, GLONASS and / or BeiDou). The GNSS receiver 100 includes an RF front end 105, a signal processing unit 110, a processor 120 and a memory 130. The RF front end 105 is configured to receive GNSS signals via the GNSS antenna 101 and output them to the signal processing unit 110. The RF front end 105 is configured to down-convert and digitize satellite signals received via the antenna 101. The RF front end essentially conditions the signal for subsequent signal processing. Other typical tasks performed by the front end include filtering and amplification. The satellite signals received at the RF front end 105 via the antenna 101 include at least one ranging signal for each of the plurality of satellites, such as a GPS L1 C / A signal. The signal processing unit 110 is configured to track the received GNSS signals (e.g., frequency, delay (code phase), and carrier phase) and generate GNSS measurements from the received GNSS signals. The processor 120 is configured to process the GNSS measurements obtained from the signal processing unit 110. Although it should be understood that there may be more than one processor within the GNSS receiver 100 to implement the method according to the present disclosure, for the purpose of this specification, it is assumed that there is only one processor 120, such as Figure 1As shown. In this example, the processor implements a navigation filter 122 as described in EP 3 859 397 A1 for the calculation of single epoch position limits. At each of a plurality of time increments (epochs), the navigation filter 122 estimates the current value of the state vector of the state variables, optionally together with their associated uncertainties. In the context of positioning, the state variables estimated by the navigation filter 122 typically include position and time variables, as well as optional speed and other variables. The memory 130 communicates with the processor 120. The memory 130 is configured to store software / firmware to be executed by the processor 120. The software / firmware is configured to control the processor 120 to perform a processing method according to the example. The memory can also be configured to store data used as input to the processor 120 and / or store data output by the processor 120.

[0123] Figure 2A A method performed by the processor 120 according to an example is illustrated.

[0124] In step 210, the processor 120 obtains a plurality of GNSS measurements from the signal processing unit 110. The GNSS measurements are all obtained from the same epoch, that is, the measurements are all made by the signal processing unit 110 at substantially the same time. The GNSS measurements include at least carrier phase measurements. In the present example, they also include pseudorange measurements. Each measurement (whether a pseudorange measurement or a carrier phase measurement) is associated with a particular GNSS signal transmitted by a particular satellite. In the present example, the plurality of GNSS measurements include pseudorange and carrier phase measurements of L1 and L2 signals of several GPS satellites. Typically, in order to calculate a position and time solution from only the GNSS measurements, measurements of at least four signals from four corresponding satellites are obtained. In general, the measurements may be associated with different satellites of the same constellation or satellites of different constellations.

[0125] Figure 2A The objective of the method is to infer state information from GNSS measurements. The state information may include values ​​of state variables, or bounds on these state variables, or both. According to the present example, the state information of primary interest is the position of the GNSS receiver and bounds on the position of the GNSS receiver. The position bounds are protection levels. In step 220, the processor 120 defines a state vector comprising state variables. The state variables include (among other things) a position variable. The state vector also includes state variables corresponding to integer ambiguities associated with corresponding carrier phase measurements.

[0126] In step 230, the processor 120 obtains the a posteriori probability density of the state vector. The a posteriori probability density is based on one or more residual error models, which describe the probability distribution of errors in each GNSS measurement result. The residual error model includes at least one non-Gaussian model.

[0127] According to this example, the errors in each GNSS measurement result are individually modeled as a univariate probability distribution. However, it should be understood that this is not necessary. In other examples, some or all GNSS measurements can be modeled together as a joint multivariate probability distribution. Traditionally, residual error models tend to assume that the errors have a Gaussian distribution. In practice, this assumption is often invalid or unreliable; however, it is chosen because it can allow the calculations required to infer state information to be greatly simplified. In an example according to the present disclosure, some (if not all) residual errors are modeled by non-Gaussian distributions. By relaxing the Gaussian assumption, a residual error model that better approximates the true distribution of the error can be designed. The simplified calculations associated with the Gaussian assumption are no longer applicable, which may increase the computational complexity; however, the example according to the present disclosure alleviates this situation by providing a technique to reduce the computational burden of the calculations associated with the non-Gaussian case. It should be understood that different non-Gaussian models may be applicable to different measurements and different application scenarios. The scope of the present disclosure is not limited to any particular non-Gaussian model (or model set). However, for the sake of clarity and completeness, some exemplary non-Gaussian error models will be described below.

[0128] The residual error model is predefined and stored in the memory 130. The residual error model for each pseudorange comprises a Student's T distribution, as described in EP 3 859 397 A1:

[0129]

[0130] Here, r is the residual, v is the degrees of freedom parameter, σ is a scaling parameter that defines the width of the core distribution, and Γ is the gamma function. This distribution may allow for better modeling of error distributions with heavy tails compared to a Gaussian (normal) distribution. The inventors have observed that the distribution of true pseudorange errors often has this type of heavy tails.

[0131] The residual error model of each carrier phase measurement result includes a periodic probability distribution. In particular, according to this example, the error model of each carrier phase is a parameter distribution of the following form:

[0132]

[0133] where r is the phase residual, σ is the scale parameter of the t-distribution (controlling the width of the peak), v is the degrees of freedom parameter of the t-distribution, and a is a normalization term chosen so that the total integrated probability over the period range (-0.5, 0.5) is 1.

[0134] It should be understood that the parameters in these error distributions can be determined empirically so that the model approximates the true error distribution as closely as possible.

[0135] In step 232, processor 120 infers state information based on the posterior probability density. This involves integrating the posterior probability density.

[0136] Figure 2B The inference in step 232 is illustrated in more detail. To simplify the explanation of the method (but without loss of generality), it will be assumed that the state information of interest is a 1-D position bound in the cross-track direction. For example, when a vehicle is traveling on a multi-lane highway, such a position bound is useful when trying to infer the lane in which the vehicle is traveling. It should be understood that if state information related to other position dimensions is required, then the above can be repeated for each dimension of interest. Figure 2B method.

[0137] In step 234, the processor 120 divides the integration domain in the dimension of interest into a plurality of slices. This means that each slice is associated with a relatively narrow integration domain defined in the dimension of interest (ie, the cross-track dimension in this example).

[0138] The slices can be selected adaptively. For example, according to the present implementation, two heuristic methods are used to select the slices. The first heuristic method determines the total width of the domain that the set of slices will cover. The second heuristic method determines the width of each individual slice.

[0139] First, a position solution is calculated (at least in the cross-track direction) using only pseudoranges. This is achieved by finding the maximum likelihood solution for the posterior probability density based only on the pseudorange measurements. The maximum likelihood solution is associated with a peak (local maximum) in the posterior based on pseudoranges. (The local maximum should also be the global maximum of the posterior; however, in some cases, the global maximum may theoretically be missing.) The processor 120 calculates the probability ("peak") associated with the maximum likelihood solution. It then finds the boundary in the cross-track direction where the posterior probability drops below a probability threshold. The probability threshold is defined at a predetermined increment (or "depth") below the peak. The boundary can be found by numerical methods. In this example, the Newton-Raphson method is used. The boundary found in this way determines the total width of the domain that the slice set will cover.

[0140] The second heuristic method is based on the "full" a posteriori probability density, i.e., the multimodal a posteriori probability based on the pseudorange measurements and the carrier phase measurements. The processor 120 finds the maximum likelihood solution based on the a posteriori probability density. This corresponds to the highest peak in the multimodal distribution. The processor 120 determines the width of the peak, in particular the full width at half maximum (FWHM). It then selects the slice width based on the determined FWHM value. In this example, the width of each slice is set to half the FWHM value.

[0141] The domain width determined by the first heuristic and the slice width determined by the second heuristic together implicitly determine the number of slices to be used. The domain width is divided into non-overlapping adjacent slices of the calculated slice widths.

[0142] This slice configuration strategy helps ensure that important parts of the posterior probability density are not missed, while avoiding the surge in slices (which may lead to an excessive increase in computational complexity). Of course, it should be understood that this slice arrangement is not restrictive. In other examples, other strategies can be adopted.

[0143] In step 236, each slice is integrated separately. Figure 2C A detailed example of this step is described. In step 238, the results of integrating the individual slices are combined by integrating these results to give an overall integrated result. The latter integration can be performed by any suitable numerical method. According to this example, a cumulative density function is formed from the slices via a numerical cumulative integration (particularly a cumulative trapezoidal numerical integration). Using the trapezoidal rule in this way provides a fairly accurate and computationally efficient way to combine the results of integrating the individual slices.

[0144] In fact, the integration of the slices is performed twice. The first integration is performed on all slices to determine the total integrated result for the entire space. The results from the individual slices can then be normalized by dividing each of the results by the total integrated result. Alternatively, as shown in this example, the first integration can be solved as a cumulative integral, as described above, to form a cumulative density function. The cumulative density function can then be normalized by dividing by the total integrated result. In either case, normalization ensures that the integral over the entire space is equal to one (as expected from a probability density function).

[0145] The second integration uses the normalized result (in this case, the normalized cumulative density function) as input. This second integration can be performed in various ways, depending on the state information required. In the first example, the goal might be to calculate the risk of having a position error greater than a specified limit (e.g., greater than 1 meter). To calculate this, the slice can be integrated (using the cumulative density function) between the limits of negative infinity and -1 meter, and between the limits of 1 meter and infinity. In the second example, a risk can be specified, and the goal might be to calculate the corresponding limit that achieves the specified risk. For example, the goal might be to calculate the risk of achieving 10 per hour. -7 The position bounds of the specified risk of error. In this case, the cumulative density function can be examined to determine the limits x such that the integral between negative infinity and -x and the integral between +x and infinity equals 10 -7 . Note, however, that the limits need not be symmetric in this way. We can do this by taking half the risk (0.5x10 -7 ) are assigned to each tail to calculate the independent upper limit U and lower limit L respectively.

[0146] Figure 2C One method for performing the integration of step 236 (for a single slice) is shown, as used in this example. This method is repeated for each slice.

[0147] First, in step 240, processor 120 searches for a significant mode of the posterior probability density (within the current slice). Then, in step 252, processor 120 numerically integrates the posterior probability density over the slice using information about the location of the significant mode (as determined in step 240).

[0148] Locating the modes of the posterior probability by preprocessing can help improve the accuracy and computational efficiency of the integration. For example, it can avoid the need to rely on numerical integration techniques to find and integrate over the relevant modes (which can be computationally expensive when using stochastic numerical integration methods and may miss modes depending on the random sample chosen). Knowing the location of the modes can benefit various integration methods. Suitable methods include: importance sampling based on the set of identified modes; or MCMC methods based on the set of identified modes. Other relevant integration methods may include the method described by Pompe et al. (Emilia Pompe, Chris Holmes, Krzysztof “A Framework for Adaptive MCMC Targeting Multimodal Distributions”, January 2019, https: / / arxiv.org / pdf / 1812.02609v2), and those cited in the introduction of that work.

[0149] Another possibility is to approximate the posterior probability density using a mathematical model that can be analytically integrated, where the approximation is based on the set of modes identified. For example, the posterior can be approximated with a Gaussian mixture model (GMM), where one or more mixture components are centered around each mode. This can avoid the need for numerical integration entirely, while still avoiding the need to assume a Gaussian distribution for the residual error.

[0150] The search in step 240 involves two stages. The first stage transforms the posterior probability density into a mixture model (see step 242). The mixture model consists of multiple mixture components, each of which is a multivariate distribution. The second stage uses the mixture model to identify the modal set (see step 244).

[0151] Figure 3 The transformation of step 242 is illustrated in more detail. First, in step 310, processor 120 defines parcel distributions for each carrier phase measurement result. In this example, these parcel distributions involve approximations of the periodic distributions of each carrier phase measurement result in the original a posteriori probability density. However, in other examples, the original posteriori may already be based on the parcel distribution of the carrier phase. In this case, no approximation is required - the parcel distribution in the original posteriori can be used to directly generate the mixture model.

[0152] In step 320, the domain of the numerical integration is transformed. In its original form, the domain of integration for each carrier phase deviation state is within a single period. However, by extracting the observations associated with the phase deviation terms from the definition of the posterior probability density, the domain of integration can be changed from this single period to a set of real numbers. With this transformation, all variables to be integrated can have the same domain - that is, the integration is performed over the set of real numbers for all variables.

[0153] In step 330, the numerical integration of the posterior probability density is reshaped as an integration over the mixture model. This is possible because the posterior is defined as the product of and. This definition is restated as a sum of products. Each such product is formed by a "linear" (i.e., non-wrapped) distribution. In this example, each of these "product" distributions is a product of univariate Student's T distributions.

[0154] After the transformation to the hybrid model in step 242, there follows a step 244 of using the hybrid model to identify the modality set. Figure 4 This is described in more detail.

[0155] In step 410, the processor 120 defines a cost function that will be used to search for a maximum likelihood solution for the state vector. The modal set can then be located near the maximum likelihood solution. The cost function is defined and the maximum likelihood solution is found by relaxing the constraint that the number of cycles in each carrier phase measurement result is an integer. Therefore, the cost function is recorded as a "floating point cost function". It is defined as a negative log-likelihood value based on the form of a single component in the mixture component.

[0156] In step 420, the processor performs a local search on the floating-point cost function to find a floating-point valued state vector associated with the minimum of the cost function. The purpose of the search is to find the global minimum, but mathematically, it is not guaranteed to be found in every case. This will depend on the shape of the cost function in each given case. In some cases, the search may (inadvertently) find another local minimum that is not the global minimum. The local search involves solving a nonlinear optimization problem. In this example, the local search is based on a hybrid between Newton's method and the steepest descent method, as described in K. Madsen, HB Nielsen, O. Tingleff, "METHODS FOR NON-LINEAR LEAST SQUARESPROBLEMS", Section 2.2 (Equation 2.11) in 2004.

[0157] The local search can be initialized in various ways. According to the present implementation, the mixture component on which the cost function is based is the product of univariate distributions. To initialize the search, each of these univariate distributions is approximated by a univariate Gaussian. The cost function then takes the form of a multivariate quadratic function, and the minimum cost solution can be found analytically, providing the starting coordinates for the search.

[0158] After finding the minimum value in step 420, the processor 120 proceeds to step 430. In this step, the processor 120 approximates the floating point cost function near the found minimum value by a multivariate Gaussian distribution. Since the cost function is expressed as a negative log-likelihood value, the approximation takes the form of a quadratic function. In the subsequent step 440, the processor 120 identifies the modal set based on the approximation using the multivariate Gaussian distribution. Now, reference will be made to Figure 5 Describe this step in more detail.

[0159] In step 441, the processor 120 characterizes the approximating Gaussian distribution by its Hessian matrix, which happens to appear as a byproduct of the local search algorithm. (The covariance matrix can be used equivalently because it is the inverse of the Hessian.) In step 442, the processor 120 defines a weighted square norm based on the Hessian / covariance matrix. In step 443, the processor 120 identifies a set of candidate integers whose distance from the local minimum found in step 420 is less than a threshold distance (according to the weighted square norm). These candidate integer sets a are possible solutions that "fix" the ambiguity in the carrier phase state to integer values. Each integer set is associated with a corresponding mode of the posterior probability density. However, the step of determining the integer has not yet defined the location of the mode in the state space. In order to locate the mode, the processor 120 defines an integer constraint cost function in step 448. This corresponds to the floating point cost function previously defined in step 410, but the integer constraint is reapplied. In step 449, processor 120 performs a local search on the cost function defined in step 448 to find the coordinates of the local minimum - in other words, the state vector that characterizes the location of the corresponding local peak in the posterior probability density. This defines the modal location in state space. Local search 449 is again performed using a hybrid between Newton's method and the steepest descent method. The optimization algorithm again produces a Hessian (or covariance) matrix as output, which characterizes the curvature of the cost function in the neighborhood of the local minimum (mode). This is useful because it characterizes the sharpness of the peak in the posterior probability density associated with that particular mode.

[0160] For each candidate integer set found in step 443, steps 448 and 449 are repeated. Figure 6 Describe the details of this step.

[0161] In the present implementation, the search for candidate integer sets is based on the LAMBDA algorithm (Teunissen P, (1995) The least-squares ambiguity decorrelation adjustment: a method for fast GPS integer ambiguity estimation, J Geod 70 (1-2): 65-82). The LAMBDA algorithm enforces integer constraints. It is run once (in step 610) to identify the maximum likelihood candidate integer set. In step 620, the weighted square norm is calculated for the maximum likelihood integer set (as previously defined in step 442). Then a second instance of the LAMBDA algorithm is run in step 630. This time, the algorithm is configured not to search for a single best integer set, but to find all integer sets whose weighted square norm differs from the weighted square norm of the best integer set by less than a search threshold. The search threshold is set empirically to ensure loose detection, i.e., all significant modes are found, and some additional, less significant modes are found slightly by mistake.

[0162] Already referred to Figures 2A to 6 The flowchart of describes the method, and we now provide the mathematical derivation for completeness.

[0163] Notation and Problem Definition

[0164] The data y contains linear observations (e.g. pseudorange measurements) and phase observations (expressed in cycles, y jk ∈[0,1]):

[0165] y i =h i (θ,z)+∈ i

[0166] y jk =h jk (θ,z)+φ j +∈ jk

[0167] There are n linear observations. The phase observations are subject to the unknown phase deviation φ j There are m different phase deviations, and n j The phase observations are affected by each bias. The model h relates states θ and z to the observations, where θ is the state of interest and z contains p linear (i.e., non-phase) interfering states. The linear observations are affected by the measurement noise ∈ i The influence of i (x), The phase observations are affected by the measurement noise ∈ jk The influence of , its probability density function is x∈[0,1]. jk (θ,z),φ j and ∈ jk They all represent phases in the range [0, 1], and the modulus of the sum of these terms is 1.

[0168] Given the observation result y and prior information I, we first use Bayes’ theorem to infer the states θ, z, and φ:

[0169]

[0170] Assuming that the prior term prob(θ,z,φ|I) is uniform, the posterior probability distribution of the state can be inferred based on the proportional constant S:

[0171]

[0172] The linear interference state z and the phase interference state φ can be marginalized:

[0173]

[0174] Note that when integrating within each slice, position variables other than the position in the dimension of interest are considered as interference states.

[0175] The proportionality constant S can be inferred from the fact that the integral of prob(θ|y,I) with respect to θ must equal 1. This is accomplished by integrating over the entire space (i.e., all slices) and using the result to normalize the density function (thus determining the value of S). This integration refers to the first integration of step 238 described above.

[0176] Then, we assume that the state information to be determined is the level of protection for a given risk. Therefore, the goal is to find a credible interval (L,U) such that the risk of the true value exceeding the interval is α, and the left and right tails share the risk equally:

[0177]

[0178] These integrals are solved in the second integral of the above step 238. That is, in the final step, when the interval is determined, they are reduced to numerical integration over the slices.

[0179] Solution

[0180] The main challenge in finding the location bounds (L,U) is that the posterior probability distribution defined in equation (1) is difficult to numerically integrate. This difficulty arises from several factors:

[0181] Typically there are about 40 states in total, including the state of interest θ, the interfering linear state z, and the phase deviation φ.

[0182] Due to the periodicity of the phase observations, the problem is highly multimodal. There may be thousands of significant modes, scattered in unpredictable locations in the high-dimensional state space.

[0183] When forming position bounds, we are concerned with quantifying the risk of rare events, with a typical acceptable risk α of about 10 per hour. -7 .

[0184] • For implementation in embedded real-time automotive systems (for example), it is desirable to be able to solve for position bounds in seconds on computing hardware suitable for mass-market applications.

[0185] These difficulties are partially addressed by splitting the integral into separate integrals for each of the multiple slices. The method is also based on a systematic search to locate the important modes (in each slice), and then importance sampling is performed in each slice to find the integral for that slice. Modal knowledge can be used to help the importance sampling algorithm move between modes efficiently and effectively.

[0186] Transform the posterior into a mixture

[0187] The posterior probability distribution given in (1) is highly multimodal, and the location of the modes in the high-dimensional state space is not obvious. Our strategy is to transform this single complex distribution into a mixture of simpler distributions in each slice. Each component in this mixture will generate a mode (or possibly more than one mode) in the posterior, and the task of locating the mode in the posterior becomes the simpler task of locating the maximum of each mixture component.

[0188] This transformation is divided into three main stages. The first stage (see Figure 3 Step 310 in (c) defines a parcel distribution that can be used to replace the existing phase error distribution; this approximation makes subsequent operations possible. Figure 3 Step 320 in Change to This is useful because this is the domain of our mixture components. The final stage (see Figure 3 Step 330 in ) reorders the terms into a mixed form (i.e., the sum of multiple simpler distributions).

[0189] Approximation using parcel distribution

[0190] The first step 310 in processing the integration is to define (for each slice) the phase measurement error model Approximation of , which is a parcel distribution:

[0191]

[0192] where f jk (x) is a linear distribution, In this implementation, the function f jk (x) is the Student's T distribution, as described above. This approximation does not need to be particularly accurate, since it is only used to locate and characterize the mode in each slice; the actual evaluation of the integral via importance sampling in each slice will use the original version.

[0193] It is helpful to consider separately the terms that depend on the phase deviation φ. Define A(θ,z) and B(θ,z,φ) so that:

[0194]

[0195] Then formula (2) becomes

[0196]

[0197] The wrapped phase measurement error distribution can be substituted into the phase term:

[0198]

[0199] This sets the stage for the rest of the transformations, as combinations of integration, product, and sum operations can be manipulated into the form we want (i.e., mixed). For convenience, we will define a notation for this term:

[0200]

[0201] Therefore, the marginalized posterior probability density for a given slice becomes

[0202]

[0203] Since each f in formula (3) jk (x) depends on only one phase deviation φ j , so the multivariable integral can be turned into a product of univariate integrals:

[0204]

[0205] We can then define the contributions from the individual phase deviations:

[0206]

[0207] Changing the integration domain

[0208] Each phase deviation φ j must have at least one associated observation (i.e., for all j, n j≥1). We will take the first of these observations (y j0 ) and extract it from the principal product term:

[0209]

[0210] We can then change the domain of integration (see Figure 3 ), which replaces the sum of the first observations of the individual phase deviations:

[0211]

[0212] In this form, a jk Only in k∈[1,n j -1], so the total number of these terms is:

[0213]

[0214] We can combine all of these into a single object

[0215] Forming a Mix

[0216] Using the distributive property, we can write the product of sums as a sum of products:

[0217]

[0218] Here, a j The collection of all a for a given j jk Item. j0 Put the terms into the sum and you get:

[0219]

[0220] We can substitute it into equation (4) and rearrange it to obtain:

[0221]

[0222] This is Figure 3 Step 330 in . The individual univariate distributions within each product term are all linear (i.e., non-wrapped) distributions. In this example, they are Student's T distributions; of course, other suitable error distribution models may be used instead of this type.

[0223] Identify significant components

[0224] Equation (5) is a mixture containing an infinite number of components, so it is not immediately obvious that the transformation will help. However, most of the mass of the distribution is usually contained in a small number of components, so the next task is to identify these significant components.

[0225] The unnormalized posterior of a single mixture component in equation (5) is:

[0226]

[0228] Identifying significant components means finding a feasible set of components Ω such that the total mass of the individual components is above a certain threshold t:

[0229]

[0230] While this is a well-specified problem, it is computationally expensive to evaluate. Instead, we approximate (6) with a multivariate normal distribution and use it to identify the approximately feasible set Ω′.

[0231] Bundling θ, z, φ, and a into a vector β, and temporarily relaxing the integer constraint on a, we obtain the nonlinear function g(β), The cost function c(β) is defined as the negative log-likelihood:

[0232] c(β)=-log(g(β))

[0233] This is the above mentioned Figure 4 The "floating point" cost function associated with step 410 in . As previously described, the minimum of this cost function is found in step 420. The minimum cost solution is defined as β * , the cost function can be approximated by the second-order Taylor expansion:

[0234]

[0235] Where H is c(β) at β = β * The Hessian matrix at .

[0236]

[0237] Please note that due to β * is a local minimum, so the gradient term is omitted in equation (7), so the gradient of c(β) is zero. Using this approximation, the cost function (negative log-likelihood) is a simple quadratic function, so the likelihood is a function with location parameter μ = β * and covariance ∑ = H -1 Multivariate normal distribution of In other words, the cost function is at its minimum value β * The vicinity is approximated by a multivariate Gaussian (step 430). The Gaussian distribution is characterized by the Hessian matrix H (see step 441). Then, by * and ∑ are truncated to the part corresponding to the a parameter (a * and∑ a), marginalization over θ, z, and φ becomes trivial:

[0238]

[0239] At this point, we can enforce the constraint that the values ​​in a should be integers, thus obtaining an approximate feasible set:

[0240]

[0241] The weighted square norm The search threshold is χ 2 The weighted square norm here is the weighted square norm described above in conjunction with step 442. Numerically, this is a problem that has been well explored in the context of integer ambiguity resolution, and efficient algorithms are available (e.g., Figure 6 610 and 630 of ). It should be appreciated that the LAMBDA algorithm is used here in a different context than its originally intended use. The original use of LAMBDA was to "fix" ambiguities, or more specifically to identify the maximum likelihood ambiguity set and then quantify the probability of an error (via a ratio test). In our case, we repurpose LAMBDA to provide an efficient search for important ambiguity sets (i.e., sets with non-negligible probabilities). Thus, in this example, the LAMBDA algorithm is used to find more than two (in fact, much more than two) ambiguity sets.

[0242] This can be done by identifying the maximum likelihood component (see step 610), finding its weighted square norm (see step 620), and then 2 Set at a fixed depth below this value to conveniently choose a reasonable search threshold χ 2 .

[0243] All that remains is to "peak" at the local maximum of the posterior probability density (or the local minimum of the integer constrained cost function) using a local search (steps 448 and 449). The final output of this local search is the state vector associated with the mode, and the Hessian matrix that characterizes the curvature of the cost function at that location.

[0244] As mentioned above, this information is then used to support inference—in this case, numerical integration of the slices.

[0245] It should be understood that the scope of the present disclosure is not limited to the above examples. Based on the above description, many variations will be apparent to those skilled in the art.

[0246] The specific parametric distribution (eg, Student's T distribution) used in the above examples is not limiting. Those skilled in the art may select other distributions (parametric or non-parametric) having the above general properties.

[0247] Similarly, different numerical methods may be used - for example, for local search (i.e., using a nonlinear optimizer other than a hybrid between Newton's method and the method of steepest descent) and / or for numerical integration. For example, suitable alternative nonlinear optimization algorithms include, but are not limited to, the Levenburg-Marquardt algorithm. In the example above, the numerical integration of the slices is performed using the trapezoidal rule. This is not required - other integration techniques may be applied.

[0248] As mentioned above, although the examples focus primarily on the computation of state bounds (particularly position bounds in the form of protection levels), slicing and modal lookups may be useful for computing state information other than (or in addition to) state bounds. For example, slicing / modal can be used to compute a maximum likelihood estimate of the state itself.

[0249] Although in the above examples, the state vector includes a position variable, and the state information to be inferred relates to a position estimate and / or a position limit, this is not restrictive. In other examples, the state vector may not include a position variable and / or the state information to be inferred may relate to other state variables. For example, the position of the GNSS receiver may be known (or actually fixed - as may be the case with a terrestrial reference station). The ultimate goal of processing GNSS measurements may be to infer information about other state variables, such as time estimates at the GNSS receiver, atmospheric parameters, errors in satellite orbits, or errors in satellite clocks. In these examples, integer ambiguities in the carrier phase measurements may still result in modes in the posterior probability density of the state variables of interest. Therefore, for the same reasons as required when trying to infer state information related to position, it may be desirable to explicitly identify such modes. Similarly, for the same reasons as before, slicing the integral may be desirable - particularly to improve the accuracy of the integral in the "tail" region of the posterior probability density.

[0250] Generally speaking, Figures 2A to 6 In the flowcharts of FIG. 1 , arrows between steps do not necessarily indicate a causal relationship between the steps. They merely indicate an exemplary order in which the steps may be performed. The method steps may be performed in an order different from the exemplary order shown in the figure.

[0251] In the above examples, each mode is characterized by the coordinates of the corresponding local maximum of the posterior probability density in the state space. Using a unique point in the state space to characterize each mode provides a compact and convenient way to specify the mode. However, this is not the only way to specify them. In other examples, one or more modes can be characterized by defining a finite region of the state space corresponding to the mode (alternatively or in addition to using the coordinates of the local maximum). For example, the region of the state space can be defined by a bounding volume, and the corresponding local maximum of the posterior probability density is located within the bounding volume. The bounding volume can be defined by specifying the upper and lower limits of the range of each state variable. This defines a box in the state space. In some cases, the bounding volume can be selected so that the posterior probability density decreases monotonically from the local maximum within the bounding volume.

[0252] In the above example, the LAMBDA algorithm is used to find a set of integer ambiguity candidates, i.e., to fix the integer ambiguities in the corresponding carrier phases. This is a convenient, efficient, effective and widely used algorithm; however, it is not mandatory to use it.

[0253] An example of a systematic search of modalities is described above. It should be understood that other systematic search methods may be used in its place. A suitable example of another algorithm for systematic search is described in European Patent Application No. 23168078.6 filed by the present applicant on April 14, 2023. EP 23168078.6 is hereby incorporated by reference into the present disclosure.

[0254] It should be understood that Figure 1 The various components shown in can be implemented in hardware, software, or a mixture of both. In addition, in a given implementation, some components can be grouped together or implemented separately. In the present implementation, block 105 is implemented entirely in hardware, block 110 is partially implemented in hardware, and the remaining components (downstream in the signal processing chain) are implemented in software. In some examples, navigation filter 122 can be implemented in a separate software module that runs on a processor different from processor 120. Other implementations are also possible that divide and distribute various functions in different ways between hardware and software or between different hardware components, software modules, and / or processors running software.

[0255] In the claims, any reference signs placed in brackets shall not be construed as limiting the claims. The word "comprising" does not exclude the presence of elements or steps other than those listed in the claim. However, when the word "comprising" is used, this also discloses as a special case the possibility that the listed elements or steps are exhaustive - that is, the device or method may consist of only these elements or steps. The word "one" before an element does not exclude the presence of multiple such elements. The embodiment may be implemented by means of hardware comprising several different elements. In a device claim that lists several means, several of these means may be implemented by the same hardware. The fact that certain measures are described in different dependent claims does not in itself indicate that the combination of these measures cannot be fully utilized. In addition, in the appended claims, a list including "at least one of A, B and C" should be interpreted as (A and / or B) and / or C.

[0256] In flowcharts, abstracts, claims, and descriptions related to methods, the order in which the steps are listed is generally not intended to limit the order in which the steps are performed. The steps may be performed in an order different from that shown (unless otherwise specified, or a subsequent step depends on the product of a previous step). However, in some cases, the order in which the steps are described may reflect a preferred order of operation.

[0257] Furthermore, in general, various embodiments may be implemented in hardware or dedicated circuits, software, logic, or any combination thereof. For example, some aspects may be implemented in hardware, while other aspects may be implemented in firmware or software, which may be executed by a controller, microprocessor, or other computing device, although these are not limiting examples. Although various aspects described herein may be illustrated and described in block diagrams, flow charts, or using other graphical representations, it is fully understood that the blocks, devices, systems, techniques, or methods described herein may be implemented in hardware, software, firmware, dedicated circuits or logic, general-purpose hardware or controllers or other computing devices, or combinations thereof, as non-limiting examples.

[0258] The embodiments described herein can be implemented by computer software, which can be executed by a data processor of a device, such as in a processor entity, or by hardware, or by a combination of software and hardware. In addition, it should be noted in this regard that any block of the logic flow in the figure can represent a program step, or interconnected logic circuits, blocks and functions, or a combination of program steps and logic circuits, blocks and functions. The software can be stored on a physical medium such as a memory chip or a memory block implemented in a processor, a magnetic medium (such as a hard disk or a floppy disk), and an optical medium (such as a DVD and its data variants or a CD).

[0259] The memory may be of any type suitable for the local technical environment and may be implemented using any suitable data storage technology, such as semiconductor-based memory devices, magnetic memory devices and systems, optical memory devices and systems, fixed memory, and removable memory. The data processor may be of any type suitable for the local technical environment and may include one or more of a general-purpose computer, a special-purpose computer, a microprocessor, a digital signal processor (DSP), an application-specific integrated circuit (ASIC), a gate-level circuit, and a processor based on a multi-core processor architecture, as non-limiting examples.

[0260] The embodiments discussed herein may be implemented in various components, such as integrated circuit modules. The design of integrated circuits is typically a highly automated process. Complex and powerful software tools are available to convert logic-level designs into semiconductor circuit designs that can be etched and formed on semiconductor substrates.

Claims

1. A method for processing a plurality of GNSS measurement results, the method comprising the following steps: Obtaining (210) the plurality of GNSS measurement results; defining (220) a state vector, wherein the state vector comprises state variables; obtaining (230) a posterior probability density for the state vector, wherein the posterior probability density is based on one or more residual error models describing a probability distribution of errors for respective ones of the GNSS measurements, the one or more residual error models comprising at least one non-Gaussian model, wherein the posterior probability density has a plurality of dimensions, each dimension corresponding to one of the state variables; and Inferring (232) state information based on the posterior probability density, wherein the step of inferring comprises integrating the posterior probability density, The step of integrating includes, for one of the state variables: partitioning (234) the domain of the posterior probability density into a plurality of slices along a dimension associated with the one state variable; integrating (236) each slice of the plurality of slices individually to produce a corresponding integration result for each slice; and The integration results of the slices are combined (238).

2. The method according to claim 1, wherein: The inferred state information includes at least one of the following: State estimation; and An error bound on the state estimate.

3. The method according to claim 1 or 2, wherein: The step of integrating each slice individually (236) includes, for each slice, performing at least one of the following: Importance sampling; MCMC methods; as well as The posterior probability density is approximated using a mathematical model that can be analytically integrated.

4. The method according to claim 1 or 2, wherein: The integration steps include, for each slice: performing a search (240) in the slice to identify the set of modes of the posterior probability density, And wherein the integration is based on the identified set of modes.

5. The method according to claim 4, wherein: The step of integrating each slice individually (236) includes, for each slice, performing at least one of the following: Importance sampling based on the identified set of modalities; MCMC methods based on the identified set of modes; and The a posteriori probability density is approximated using a mathematical model that can be analytically integrated, wherein the mathematical model is based on the identified set of modes.

6. The method according to claim 4, wherein: The search is a system search.

7. The method according to claim 4, wherein: The search includes the following steps: transforming (242) the posterior probability density into a mixture model comprising a plurality of mixture components, wherein each mixture component is a multivariate distribution; and identifying (244) the set of modalities using the mixture model, Therein, each mode is associated with a corresponding multivariate distribution in the multivariate distributions.

8. The method according to claim 7, wherein: The step of transforming (242) the posterior probability density comprises: defining (310) a packet distribution for each carrier phase measurement; and transforming (330) the parcel distribution into the hybrid model, wherein each mode is associated with a set of integers a, and each integer is associated with a corresponding carrier phase measurement result in the carrier phase measurement results, Wherein, each integer indicates the number of cycles in the corresponding package distribution.

9. The method according to claim 7 or 8, wherein: The step of identifying (244) the set of modalities comprises: defining (410) a floating point cost function based on the mixture model by relaxing the constraint that each value a indicating the number of cycles in the corresponding carrier phase measurement is an integer; performing (420) a first local search of the floating point cost function to find a floating point valued state vector associated with a local minimum of the cost function; Approximating (430) the floating point cost function in the region of the local minimum via a multivariate Gaussian distribution; and The set of modes is identified (440) based on approximation using a multivariate Gaussian distribution.

10. The method according to claim 9, wherein: The approximation is expressed in the form of a negative log-likelihood using a multivariate Gaussian distribution, where the multivariate Gaussian distribution is modeled as a quadratic function.

11. The method according to claim 9, wherein: The step of identifying (440) the set of modes based on an approximation with a multivariate Gaussian distribution comprises characterizing (441) the approximation with a multivariate Gaussian distribution in the region of the local minimum by a matrix, the matrix being one of: a covariance matrix and a Hessian matrix.

12. The method according to claim 11, wherein: The step of identifying (440) the set of modes using a multivariate Gaussian distribution based on approximation also includes: based on a state vector associated with the local minimum and the matrix, defining (442) a weighted square norm based on the matrix; and A set of candidate integers a having a distance from the local minimum less than a threshold value is identified (443) based on the weighted square norm.

13. The method according to claim 12, wherein: The step of identifying (443) a set of candidate integers comprises: identifying (610) a maximum likelihood set of candidate integers; finding (620) a first weighted square norm associated with the maximum likelihood set of candidate integers; and Other candidate sets of integers having weighted squared norms within a predetermined search threshold of the first weighted squared norm are identified (630).

14. A computer program product comprising a computer program, the computer program comprising computer program code, the computer program code being configured to cause one or more processors to perform all the steps of the method according to any one of claims 1 to 13 when the computer program is run on the one or more processors.

15. A GNSS receiver (100), comprising: a signal processing unit (110), the signal processing unit being configured to generate a plurality of GNSS measurement results; as well as At least one processor (120), the at least one processor being configured to: Obtaining (210) the plurality of GNSS measurement results; defining (220) a state vector, wherein the state vector comprises state variables; obtaining (230) a posterior probability density for the state vector, wherein the posterior probability density is based on one or more residual error models describing a probability distribution of errors for respective ones of the GNSS measurements, the one or more residual error models comprising at least one non-Gaussian model, wherein the posterior probability density has a plurality of dimensions, each dimension corresponding to one of the state variables; and Inferring (232) state information based on the posterior probability density, wherein the step of inferring comprises integrating the posterior probability density, The step of integrating includes, for one of the state variables: partitioning (234) the domain of the posterior probability density into a plurality of slices along a dimension associated with the one state variable; integrating (236) each slice of the plurality of slices individually to produce a corresponding integration result for each slice; and The integration results of the slices are combined (238).

Images