GNSS measurement processing to identify modalities

Abstract

The GNSS measurements are processed to identify modalities. A method and apparatus for processing GNSS measurements is disclosed. The GNSS measurements include carrier phase measurements. A state vector is defined that includes a state variable. A posterior probability density of the state vector is obtained, which is based on a non-Gaussian residual error model of the GNSS measurements. A search of posterior probability densities is performed to identify a modal set of probability densities. Using the identified set of modalities, state information is inferred based on a posterior probability density. The searching includes converting the posterior probability density into a hybrid model containing a plurality of hybrid components, where each hybrid component is a multivariable distribution.

Description

Technical Field

[0001] The present invention relates to Global Navigation Satellite Systems (GNSS). In particular, the present invention relates to methods and apparatus for processing GNSS measurements. 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 in the accuracy 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 a posteriori probability density as accurately and efficiently as possible. In this regard, a challenge arises when the GNSS measurements include carrier phase measurements. Due to the integer ambiguities inherent in each carrier phase measurement, the a posteriori probability density (e.g., over position) is multimodal. Each distribution mode corresponds to a possible combination of integer values ​​for resolving (or "fixing") the ambiguities.

[0012] 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.

[0013] A method and apparatus for processing GNSS measurements are disclosed. The GNSS measurements include carrier phase measurements. A state vector is defined, which includes state variables. A posterior probability density of the state vector is obtained, which is based on a non-Gaussian residual error model of the GNSS measurements. A search of the posterior probability density is performed to identify a set of modes of the probability density. Using the identified set of modes, state information is inferred based on the posterior probability density. The search includes transforming the posterior probability density into a mixture model including a plurality of mixture components, wherein each mixture component is a multivariate distribution.

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

[0015] Obtaining the plurality of GNSS measurement results, wherein the plurality of GNSS measurement results include a plurality of carrier phase measurement results;

[0016] defining a state vector, the state vector comprising state variables;

[0017] 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 of respective ones of the GNSS measurement results, the one or more residual error models comprising at least one non-Gaussian model; and

[0018] Performing a search to identify a set of modes of the posterior probability density, wherein the step of searching comprises:

[0019] transforming the posterior probability density into a mixture model comprising a plurality of mixture components, wherein each mixture component is a multivariate distribution; and

[0020] identifying the set of modalities using the mixture model,

[0021] wherein each mode is associated with a corresponding one of the multivariate distributions; and

[0022] Using the identified set of modalities, state information is inferred based on the posterior probability density.

[0023] 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).

[0024] 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.

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

[0026] 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).

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

[0028] 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.

[0029] 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.

[0030] 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.

[0031] 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).

[0032] The step of identifying the set of modes optionally includes: 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.

[0033] 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 peak of each Gaussian approximation. The peak of each Gaussian approximation may be found analytically.

[0034] 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.

[0035] 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*.

[0036] 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.

[0037] 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.

[0038] 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.

[0039] 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.

[0040] 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.

[0041] 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.

[0042] 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.

[0043] 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.

[0044] 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.

[0045] 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.

[0046] 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.

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

[0048] The step of inferring may include numerically integrating the posterior probability density including the at least one non-Gaussian model, wherein the integration is based on the identified set of modalities.

[0049] The step of numerical integration optionally includes importance sampling based on the identified set of modes.

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

[0051] 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.

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

[0053] a signal processing unit configured to generate a plurality of GNSS measurements, wherein the plurality of GNSS measurements include a plurality of carrier phase measurements; and

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

[0055] obtaining the plurality of GNSS measurement results;

[0056] defining a state vector, the state vector comprising state variables;

[0057] 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 of respective GNSS measurement results in the GNSS measurement results, the one or more residual error models comprising at least one non-Gaussian model;

[0058] performing a search to identify a set of modes of the posterior probability density; and

[0059] Using the identified set of modalities, inferring state information based on the posterior probability density,

[0060] The searching steps include:

[0061] transforming the posterior probability density into a mixture model comprising a plurality of mixture components, wherein each mixture component is a multivariate distribution; and

[0062] identifying the set of modalities using the mixture model,

[0063] Therein, each mode is associated with a corresponding multivariate distribution in the multivariate distributions.

[0064] 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

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

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

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

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

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

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

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

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

[0073] 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.

[0074] 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).

[0075] Figure 11 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 1 As 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.

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

[0077] 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 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.

[0078] Figure 2 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.

[0079] 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.

[0080] 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.

[0081] 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:

[0082]

[0083] 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.

[0084] 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:

[0085]

[0086] 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.

[0087] 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.

[0088] In step 240, processor 120 searches for significant modes of the posterior probability density. Then, in step 250, processor 120 uses the information about the location of the significant modes determined in step 240 to infer state information (particularly location and protection level) based on the posterior probability density.

[0089] According to the present example, the inference in step 250 involves numerically integrating the posterior probability density (see step 252).

[0090] 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.

[0091] 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.

[0092] 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).

[0093] Figure 3The 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.

[0094] 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.

[0095] 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.

[0096] 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.

[0097] 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.

[0098] 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.

[0099] 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.

[0100] 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.

[0101] 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.

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

[0103] 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.

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

[0105] Notation and Problem Definition

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

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

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

[0109] 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.

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

[0111]

[0112] 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:

[0113]

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

[0115]

[0116] The proportionality constant S can be inferred from the fact that the integral of prob(θ|y,I) with respect to θ must equal 1. The goal is then to find a credible interval (L,U) such that the risk of the true value being outside the interval is α, and the left and right tails share this risk equally:

[0117]

[0118] Solution

[0119] 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:

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

[0121] 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.

[0122] 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 .

[0123] • 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.

[0124] The method locates important modes based on a systematic search and then performs importance sampling to find the integral. The knowledge of the modes can be used to help the importance sampling algorithm move between modes efficiently and effectively.

[0125] Transform the posterior into a mixture

[0126] The posterior probability distribution given by equation (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. 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.

[0127] 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).

[0128] Approximation using parcel distribution

[0129] The first step 310 in processing the integral is to define the phase measurement error model Approximation of , which is a parcel distribution:

[0130]

[0131] 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 modes; the actual evaluation of the integral via importance sampling will use the original version.

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

[0133]

[0134] Then formula (2) becomes

[0135]

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

[0137]

[0138] 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:

[0139]

[0140] Therefore, the marginalized posterior probability density becomes

[0141]

[0142] 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:

[0143]

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

[0145]

[0146] Changing the integration domain

[0147] 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:

[0148]

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

[0150]

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

[0152]

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

[0154] Forming a Mix

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

[0156]

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

[0158]

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

[0160]

[0161] 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.

[0162] Identify significant components

[0163] 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.

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

[0165]

[0166] 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:

[0167]

[0168] 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 Ω′.

[0169] 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:

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

[0171] 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:

[0172]

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

[0174]

[0175] 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:

[0176]

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

[0178]

[0179] 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 the proposed LAMBDA algorithm).

[0180] 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 .

[0181] 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.

[0182] As mentioned above, this information is then used to support inference—either numerical integration or integration by other means.

[0183] 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.

[0184] It should be understood that 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 distributions or non-parametric distributions) having the above general characteristics.

[0185] 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 steepest descent method) and / or for numerical integration. For example, suitable alternative nonlinear optimization algorithms include, but are not limited to, the Levenburg-Marquardt algorithm.

[0186] As mentioned above, although the examples focus primarily on the computation of state bounds (particularly position bounds in the form of protection levels), modal lookups are useful for computing state information other than (or in addition to) state bounds. In particular, modes may be useful for computing maximum likelihood estimates of the state itself.

[0187] 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 bound, this is not limiting. 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 the GNSS measurements may be to infer information about other state variables, such as 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, it may be desirable to explicitly identify such modes for the same reasons as required when attempting to infer state information related to position.

[0188] Generally speaking, Figures 2 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.

[0189] 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.

[0190] 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.

[0191] It should be understood that Figure 1The 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.

[0192] 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.

[0193] 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.

[0194] 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.

[0195] 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).

[0196] 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.

[0197] 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 measurements, wherein the plurality of GNSS measurements comprises a plurality of carrier phase measurements; defining (220) a state vector, the state vector comprising state variables; obtaining (230) a posterior probability density of 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; and A search (240) is performed to identify a set of modes of the posterior probability density, wherein the step of searching comprises: 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, wherein each mode is associated with a corresponding one of the multivariate distributions; and Using the identified set of modalities, state information is inferred (250) based on the posterior probability density.

2. The method according to claim 1, 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.

3. The method according to claim 1 or 2, 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.

4. The method according to claim 3, 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.

5. The method according to claim 3, 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.

6. The method according to claim 5, 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.

7. The method according to claim 6, 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).

8. The method according to claim 6 or 7, further comprising, for each identified candidate integer set: Another cost function is defined (448) 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 set of candidate integers; and A second local search is performed (449) on the other cost function to find a state vector associated with a local minimum of the cost function associated with the set of candidate integers.

9. The method according to claim 1 or 2, wherein: Before transforming the obtained posterior probability density into the mixture model, the obtained posterior probability density is defined as the product of: a plurality of first residual error models, wherein the plurality of first residual error models describe probability distributions of errors in corresponding pseudorange measurement results; a plurality of second residual error models describing probability distributions of errors in corresponding first carrier phase measurements; as well as a plurality of third residual error models, the plurality of third residual error models describing a probability distribution of errors in corresponding second carrier phase measurements, The first carrier phase measurement result includes a carrier phase measurement result for each visible GNSS constellation and for each GNSS frequency band, The second carrier phase measurement result includes 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.

10. The method according to claim 1 or 2, wherein: The state vector includes a phase deviation for each GNSS constellation and each GNSS frequency band.

11. The method according to claim 1 or 2, wherein: The step of inferring includes numerically integrating (252) the posterior probability density including the at least one non-Gaussian model, wherein the integration is based on the identified set of modalities.

12. The method according to claim 11, wherein: The step of performing numerical integration (252) includes performing importance sampling based on the identified set of modes.

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

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

15. A GNSS receiver (100), comprising: a signal processing unit (110) configured to generate a plurality of GNSS measurements, wherein the plurality of GNSS measurements include a plurality of carrier phase measurements; and 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, the state vector comprising state variables; obtaining (230) a posterior probability density of 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; performing (240) a search to identify a set of modes of the posterior probability density; and Using the identified set of modalities, state information is inferred (250) based on the posterior probability density, wherein the step of searching comprises: 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.

Images