Adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution

Abstract

The present invention discloses an adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution. The method does not make an assumption about the number of satellites that have faults or deviations, but modifies the Gaussian distribution assumption of errors in a traditional positioning model and converts it into a Gaussian mixture distribution, i.e., a non-Gaussian distribution. Parameters are solved by a maximum likelihood method based on real-time observation data, and the probability distribution of the deviation range of each satellite fault is adaptively calculated. The method can be used for identification and robust positioning under multiple faults, and can also be used for robust positioning and multipath elimination in a multipath environment. The present invention solves the shortcomings of existing RAIM and robust positioning methods. Currently, existing fault detection and identification methods have a large amount of computation, rely on a priori assumptions about the number of faults, and have a high misjudgment rate. The method provided by the present invention has a high fault detection success rate, a small amount of computation, high speed, high accuracy, and robustness.

Description

Technical Field

[0001] The invention belongs to the technical field of satellite navigation and positioning, and relates to an adaptive real-time multipath elimination and anti-error positioning method based on non-Gaussian distribution. Background Art

[0002] When using satellite navigation systems to locate a receiver, the pseudorange or carrier measurements observed by the receiver may contain propagation errors, such as electronic failures, errors in the satellite-broadcast ephemeris and clock, abnormal atmospheric delays, multipath effects, and receiver failures. These failures can further lead to large errors in the positioning results output by the user's receiver, which may adversely affect the user's own safety. With the rapid development and widespread use of global satellite navigation systems, real-time, rapid integrity fault monitoring and robust positioning methods are playing an increasingly important role.

[0003] RAIM (Receiver Autonomous Integrity Monitoring) involves real-time monitoring of individual satellite observations by the receiver itself. It leverages the redundancy created by multiple satellite observations received by the receiver to detect and identify faults. Its advantages are that it requires no external equipment, is low-cost, and easy to implement. However, its disadvantages are that once a faulty satellite is identified, it must be removed, increasing computation time. Furthermore, the reduced number of satellites can lead to reduced positioning accuracy.

[0004] Currently, widely used RAIM methods fall into two main categories: residual analysis and maximum solution separation. Residual analysis is based on least squares residuals. After least squares positioning, the sum of squared residuals (or studentized residuals) is calculated and compared with a threshold. However, the least squares method is very sensitive to faulty satellites. When a faulty satellite has a large error, the positioning result will be significantly offset, causing the estimated absolute value of the residual for the faulty satellite to decrease while the absolute value of the residual for non-faulty satellites to increase, ultimately leading to misjudgment and inability to accurately identify the true faulty satellite. The maximum solution separation method divides all observations into several subsets based on a certain method. Positioning is performed on each observation subset separately, and comparisons and judgments are made based on statistics such as the sum of squared residuals after positioning. For example, the optimal combination is the one with no faulty satellites. However, this method requires artificial assumptions about the number of faulty satellites and requires a separate solution for each combination, which is computationally intensive. Incorrect assumptions about the number of faulty satellites can also lead to misjudgments.

[0005] Compared to the RAIM method, the robust positioning method does not exclude faulty satellites. It is a robust positioning method that overcomes the disadvantage of the least squares method, where the residual sum of squares is very sensitive to large errors. By constructing a robust loss function, the receiver can use all satellite observations for positioning while reducing the impact of faulty observations on the positioning results. However, its disadvantage is that the error distribution assumption of the robust positioning method is thicker than the normal distribution. Therefore, when there are no faulty satellites, the positioning results obtained by the method will be larger than those of the least squares method and are often not unbiased estimates.

[0006] When satellite signals received by a receiver are obscured or reflected, the observed values ​​contain significant errors. This type of multipath error can be considered a special type of fault error: fault errors occur simultaneously in multiple satellite observations at the same time, resulting in a greater number of faults than other fault modes. Current RAIM fault detection methods are mostly designed for faults of no more than three or four. However, in multi-mode, multi-frequency positioning, the number of available satellites is large, and the number of multipath faults far exceeds three or four. Identifying and eliminating multipath errors is crucial for the magnitude of receiver positioning errors, and fault elimination methods are needed to address a greater number of faults. Summary of the Invention

[0007] The present invention aims to overcome the above-mentioned drawbacks and provide an adaptive real-time multipath elimination and robust positioning method based on a non-Gaussian distribution. The method does not make any assumptions about the number of satellites experiencing failures or deviations, but instead modifies the Gaussian distribution assumption of errors in the traditional positioning model and converts it into a Gaussian mixture distribution, i.e., a non-Gaussian distribution. Parameters are solved by a maximum likelihood method based on real-time observation data, and the probability distribution of the deviation range of each satellite failure is adaptively calculated. The method can be used for identification and robust positioning under multiple faults, as well as for robust positioning and multipath elimination in a multipath environment. The present invention overcomes the shortcomings of existing RAIM and robust positioning methods. Currently, existing fault detection and identification methods have a large amount of computation, rely on a priori assumptions about the number of faults, and have a high misjudgment rate. The method provided by the present invention has a high fault detection success rate, is fast and has low computational complexity, and is highly accurate and robust.

[0008] In order to achieve the above-mentioned object of the invention, the present invention provides the following technical solutions:

[0009] An adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution comprises the following steps:

[0010] (1) Establish the positioning solution model of the receiver within the positioning period and convert it linearly to obtain the standardized positioning solution model Y = Hβ + ε, ε ~ N (0, σ 2 I), where Y=(y1,...,yn ) T is the standardized transformation of the known satellite observation matrix, H=(h1,...,h n ) T is the normalized transformation of the n×p-dimensional known observation positioning geometry matrix, n is the total number of satellites, β is the p×1-dimensional unknown parameter vector, and ε is the normalized transformation of the n×1-dimensional observation noise matrix, which obeys the normal distribution N(0,σ 2 I),σ 2 is the variance coefficient, I is the unit diagonal matrix;

[0011] (2) Use Gaussian mixture distribution to model the satellite observations and define the n-dimensional unknown vector Z as a latent variable, where each element Z i ∈{1,2,3},i=1,2,...,n, and obtain the Gaussian mixture distribution model of each satellite observation where γ k Z i The prior probability of

[0012] make represents the parameter vector to be solved consisting of all unknown parameters to be solved in the Gaussian mixture distribution model of each satellite observation, where

[0013] μ=(μ1,...,μ K ) T ,μ1=0,σ 2 =(σ 2 1,...,σ 2 K ) T ,γ=(γ1,...,γ K ) T ,K=3,γ3=1-γ1-γ2.

[0014] (3) Iteratively solve the parameter vector θ based on the EM algorithm to obtain the optimal value of the parameter vector;

[0015]

[0016] (4) According to the optimal value of the parameter vector to be solved Calculate the real-time protection level of positioning and compare it with the given alarm threshold to determine whether the positioning result is reliable;

[0017] (5) The optimal value of the parameter vector to be solved in Substitute the obtained standardized positioning solution model in step (1) to obtain the test statistic. Based on the comparison between the test statistic and the rejection region of the hypothesis test, it is determined whether each satellite has a fault that affects positioning.

[0018] Furthermore, in step (2), the Gaussian mixture distribution model of each satellite observation is established by the following method:

[0019] Define n-dimensional unknown vector Z as hidden variable, where each element Z i ∈{1,2,3},i=1,2,...,n,Z i =1 indicates observation noise ε i Does not contain deviations that cause unacceptable errors in receiver positioning, Z i =2 or Z i =3 respectively represent the observation noise ε i Contains varying degrees of deviation that can cause unacceptable errors in receiver positioning;

[0020] Assuming that different Z i Under the value, the satellite observation quantity y i Obey different Gaussian distributions, namely

[0021]

[0022]

[0023] where y i is the i-th element of the known satellite observation matrix Y, h i The transpose of the i-th row of the geometry matrix H for the known observation location, h i T Indicates h i The transpose of f k (y i ) is a given Z i Time i The conditional Gaussian distribution of Z i Unknown time i The mixed Gaussian density function f m (y i ) in the k-th Gaussian component;

[0024] When Z i When unknown, the satellite observation quantity y i For Gaussian mixture distribution, the Gaussian mixture distribution model is as follows:

[0025]

[0026] Among them, γ k Z i The prior probability, P(Z i =k) ​​=γ k .

[0027] Furthermore, in step (3), the parameter vector θ to be solved is iteratively solved based on the EM algorithm to obtain the optimal value of the parameter vector to be solved The method comprises the following steps:

[0028] (31) Calculate the expected log-likelihood function of the known satellite observation matrix Y:

[0029]

[0030]

[0031] in t is the number of iterations of the parameter vector θ to be solved;

[0032] (32) Set the initial value θ of the parameter vector θ to be solved (0) , at this time, take t = 0;

[0033] (33) According to θ (t) Substituting in, using the Bayesian method, we can calculate

[0034] (34) Update t to t+1 and change the value calculated in step (33) to Substitute the expected log-likelihood function Q(θ|θ (t-1) ,Y), solve Let Q(θ|θ (t-1) ,Y) is maximized, and we get

[0035] (35) As γ k Substitute the value of into the expected log-likelihood function Q(θ|θ (t-1) ,Y), solve Let Q(θ|θ (t-1) ,Y)maximization;

[0036] (36) Repeat steps (33), (34) and (35), each time updating t to t+1 before starting step (34), until the θ calculated in step (35) is (t) and θ (t-1) The difference is less than the given first threshold value, and the final iteration number t is obtained. (t) , the optimal value of the parameter vector to be solved

[0037] Furthermore, in step (32), the iterative initial value of the parameter vector θ to be solved is where β (0) Estimated from the state equation or state prediction in Kalman filtering, or determined using robust positioning methods, μ (0) ,(σ2 ) (0) and γ (0) , and determined based on historical experience data.

[0038] Furthermore, in step (35), solve Let Q(θ|θ (t-1) ,Y) maximization steps are as follows:

[0039] (351)Fixed(σ 2 ) (t) , find the value that makes Q(θ|θ (t-1) ,Y) reaches the maximum

[0040]

[0041]

[0042] e n =(1,...,1) T

[0043]

[0044]

[0045] Among them use or The latest value of the two, steps (351) and (352) need to be repeated multiple times, each time updating and Until the difference between the two updated values ​​is less than a given first threshold;

[0046] (352) Under the condition of (t-1) ,Y) reaches the maximum (σ 2 ) (t) ;

[0047]

[0048] Furthermore, in step (36), the first threshold value is determined as β and σ calculated by multiple experiments. 2 The average error is 5%.

[0049] Furthermore, in the step (4), the real-time protective level of the positioning includes a horizontal protective level HPL and a vertical protective level VPL;

[0050]

[0051] Among them, κ α is the variance expansion coefficient, which is set according to the empirical value or the actual false alarm and missed alarm probability. 11 , C 22 and C 33 are the optimal values ​​of the parameter vectors to be solved middle The variance of the 1st, 2nd, and 3rd elements of ;

[0052] The alarm threshold is a vector with a length of 2, including threshold values ​​in the horizontal direction and the vertical direction.

[0053] Furthermore, the The variance of each element of is the matrix The first p diagonal elements in ;

[0054] Among them A E =(H E T W E H E ) -1 H E T W E , Σ Y is the posterior covariance matrix of the known satellite observation matrix Y, which is a diagonal matrix. The i-th element in the diagonal element is:

[0055] where μ k , are the optimal values ​​of the parameter vectors to be solved in step (3) The value of the corresponding element in ,γ ik According to γ (t) The value of the corresponding element is calculated by step (33).

[0056] Furthermore, in step (5), the optimal value of the parameter vector to be solved is in Substitute into the standardized positioning solution model obtained in step (1) and calculate The resulting test statistic is: in It represents the median of the absolute values ​​of all residuals (i.e., the estimated values ​​of the errors).

[0057] Furthermore, in step (5), the rejection region of the hypothesis test is {T i ||T i |>T}, |T i |>T indicates that the i-th satellite has a fault that affects positioning; otherwise, it indicates that the i-th satellite does not have a fault that affects positioning;

[0058] The T is the third threshold value, T = Φ -1 (1-α / 2), T is the α quantile of the standard normal distribution, α is the significance level determined after comprehensively considering the false alarm and missed alarm rates, for example, it can be 5%, Φ -1 is the inverse of the standard normal cumulative distribution function.

[0059] Furthermore, in step (1), the positioning solution model of the receiver during the positioning period is in is the n×1-dimensional known observation quantity, n is the total number of satellites, is the n×p dimensional known observation positioning geometry matrix, β is the p×1 dimensional unknown parameter vector, is the observation error vector, which is an n×1 dimensional vector, Obey the normal distribution N(0,σ 2 ∑),σ 2 is the variance coefficient, ∑ is the covariance prior matrix of each satellite observation error,

[0060] The linearization conversion process is as follows:

[0061]

[0062] After linearization conversion I is the unit diagonal matrix;

[0063] We get Y=Hβ+ε,ε~N(0,σ 2 I);

[0064] in

[0065] Compared with the prior art, the present invention has the following beneficial effects:

[0066] (1) The present invention provides an adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution, which combines RAIM, multipath elimination, and positioning. It reduces the impact of faulty or biased satellite observations on positioning through adaptive weights based on fault probability and fault error size. Fault detection, identification, and elimination can be performed simultaneously with positioning calculations, with fewer computational steps and strong real-time performance.

[0067] (2) The adaptive real-time multipath elimination and anti-error positioning method based on non-Gaussian distribution of the present invention has an adaptive weight in positioning, which is more robust than other weighted least squares solution methods. Other weighted least squares methods are very sensitive to weight settings. When the weight settings deviate slightly from the actual ones, the positioning results will be greatly offset, and the faulty satellite cannot be successfully identified. The weight calculation of this method is based on real-time observation data and adaptive iterative calculation. It can adaptively adjust the weight of the faulty satellite according to the size of the fault deviation, thereby ensuring good positioning accuracy, fault detection success rate and robustness.

[0068] (3) The present invention provides an adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution. In terms of identifying faulty satellites, it does not only have two states, yes and no, but can also provide the probability of failure and the size of the failure range for each satellite based on observation data.

[0069] (4) The present invention provides an adaptive real-time multipath elimination and robust positioning method based on a non-Gaussian distribution, which does not require any prior assumptions about the number of faulty satellites. Other methods require separate observation subset calculations and comparisons under different assumptions about the number of faulty satellites, which results in many calculations and is particularly time-consuming when there are a large number of faulty or deviated satellites. This method modifies the Gaussian distribution assumption of errors in traditional positioning models, converting it to a Gaussian mixture distribution. This method directly estimates the probability of each satellite failing, eliminating the need to exclude faulty satellites and facilitating faster calculations. DETAILED DESCRIPTION

[0070] The following detailed description of the present invention will make the features and advantages of the present invention more clear and explicit.

[0071] The word “exemplary” is used exclusively herein to mean “serving as an example, example, or illustration.” Any embodiment described herein as “exemplary” is not necessarily to be construed as preferred or advantageous over other embodiments.

[0072] To address the shortcomings of existing RAIM and robust positioning methods, the present invention proposes an adaptive, real-time multipath mitigation and robust positioning method based on a non-Gaussian distribution. This method performs multipath mitigation and multiple fault detection on observed quantities simultaneously with positioning, avoiding the need for other methods to perform fault detection followed by re-positioning, significantly reducing computational time. Furthermore, it adds state estimation of each satellite's fault presence, fault error magnitude, and fault probability. Iteratively estimates the posterior probability based on observed data, and solves the problem using an adaptive weighted least squares method based on the fault probability and fault error. Compared to other least squares positioning methods based on fixed weights, this method automatically adjusts satellite weights based on the satellite fault probability and error magnitude, exhibiting excellent adaptability to both the number of faulty satellites and error magnitude. While ensuring accuracy, it significantly reduces the impact of faulty satellites on positioning results, ensuring robustness. Furthermore, this method avoids the need to make assumptions about the number of faulty satellites and then perform positioning calculations and comparisons under different assumptions, instead directly solving the optimal expected logarithmic maximum likelihood based on the satellite fault probability during the calculation, significantly reducing computational effort. Finally, this method can also estimate the probability of failure of each satellite, overcoming the shortcomings of other methods that can only determine whether a satellite has failed but cannot estimate the failure probability.

[0073] like Figure 1 As shown, an adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution includes the following steps:

[0074] Step 1. In each positioning cycle, establish the receiver positioning solution model and standardize it.

[0075] The positioning equation is in is a known observation, including pseudorange or carrier observation, which is an n×1 dimensional vector, where n is the total number of satellites. is the known observation positioning geometry matrix, which is n×p dimensional, and β is the p×1 dimensional unknown parameter vector, which is the receiver parameter to be solved, including the three-dimensional position coordinates, the receiver clock error or velocity parameter to be solved, is the observation error vector, which is an n×1 dimensional vector, Obey the normal distribution N(0,σ 2 ∑),σ 2 is the unknown variance coefficient, ∑ represents the covariance prior matrix of the errors of each satellite observation given in advance. It is generally assumed that the observations of different satellites are independent of each other. ∑ is a positive definite diagonal matrix.

[0076] The value of ∑ can be determined based on multiple experiments in advance, such as a function related to the satellite pitch angle, carrier-to-noise ratio, etc., or can be simply set to a unit matrix.

[0077] In order to facilitate subsequent calculations, a linear transformation is applied to the receiver positioning solution model to standardize it, and a standardized receiver positioning solution model is obtained:

[0078]

[0079]

[0080] After transformation I is the unit diagonal matrix.

[0081] We get Y=Hβ+ε,ε~N(0,σ 2 I);

[0082] in

[0083] Step 2. Use Gaussian mixture distribution to model the observations of each satellite, which is equivalent to modeling the errors.

[0084] Define ε i is any element of ε, 1≤i≤n, assuming that the observation noise of each satellite is ε i Since the observation noise ε may be affected by multipath, signal shielding, satellite or receiver failure, etc., we define an n-dimensional unknown vector Z, where Z is a hidden variable and each element Z i ∈{1,2,3},i=1,2,...,n,Z i Different values ​​of represent different failure modes, such as no failure or deviation, mean shift, variance increase, etc. Specifically, Z i =1 indicates observation noise ε i Does not contain faults or deviations that cause large errors in receiver positioning, Z i =2 or Z i =3 respectively represent the observation noise ε i There are different degrees of deviations that lead to large errors in receiver positioning. In this case, there may be a drift in the mean or an increase in the variance.

[0085] Assuming that different Z i Under the value, the observation y i Obey different Gaussian distributions:

[0086]

[0087]

[0088] where y i is the i-th element of the known observation Y after normalization in step (1), h iis the transpose of the i-th row of the normalized known observation positioning geometry matrix H, h i T Indicates h i The transpose of .

[0089] When Z i When unknown, y i It is a Gaussian mixture distribution, and the mixed Gaussian distribution model of the observation quantity, that is, its density function is as follows:

[0090]

[0091] Among them, γ k Z i The prior probability, P(Z i =k) ​​=γ k .

[0092] The parameter vector to be solved is The last 7 parameters are the parameters of the Gaussian mixture distribution. It should be noted that we assume μ1 = 0, that is, the first component of the Gaussian mixture distribution has zero mean, and the Gaussian distribution f1(y i ), is the error distribution when there is no fault or deviation. γ3=1-γ1-γ2 is the redundant parameter.

[0093] make Represents a vector consisting of all unknown parameters to be solved, where the superscript T in the upper right corner represents the transpose of the vector.

[0094] μ=(μ1,...,μ K ) T ,μ1=0,σ 2 =(σ 2 1,...,σ 2 K ) T ,γ=(γ1,...,γ K ) T ,K=3,γ3=1-γ1-γ2.

[0095] When multiple multipath errors coexist, this method uses a non-Gaussian assumption for the error distribution. Specifically, it uses a Gaussian mixture distribution to model the errors, making it more suitable for the error distribution characteristics of multiple faults under multipath effects. Unlike other positioning methods, this method does not require a priori assumptions about the number of satellites experiencing faults or deviations. Other integrity troubleshooting methods perform separate positioning solutions for different numbers of faulty satellites, which results in a high computational load.

[0096] The present invention does not make any assumptions about the number of satellites that have failed or deviated. Instead, it modifies the Gaussian distribution assumption of errors in the traditional positioning model and converts it into a Gaussian mixture distribution, i.e., a non-Gaussian distribution. Parameters are solved using the maximum likelihood method based on real-time observation data, and the probability distribution of the deviation range of each satellite failure is adaptively calculated. This avoids the reduction in solution efficiency caused by improper artificial prior assumptions, ensures that the positioning results are both accurate and robust, and eliminates deviations in satellite observation values ​​caused by failures or multipath.

[0097] At the same time, unknown parameters such as the standard error and deviation size of the faulty satellite are added. While using all satellites for positioning solution, the standard error and deviation size of the fault-free satellite and the faulty satellite are estimated respectively, which reduces the impact of the faulty satellite and improves positioning accuracy and robustness.

[0098] Step 3. Calculate the expected log-likelihood function.

[0099] By finding the expectation of the posterior distribution about Z, we can obtain the expected log-likelihood function of the observed value, that is,

[0100]

[0101] in Will be calculated in step 5.

[0102] t is the number of iterations, the initial value of iteration θ (0) For selection, see step 4. Indicates failure mode z i =k is the posterior probability calculated at the tth iteration based on the real-time observation data.

[0103] This invention differs from other positioning methods in that, under different assumptions about the number of faulty satellites, these robust positioning methods primarily calculate the log-likelihood for healthy satellites, significantly reducing the weight of faulty satellites. This failure to utilize all information leads to an underestimation of the number of satellites located. Some multipath mitigation methods also estimate multipath errors, but fail to utilize the real-time positioning matrix, introducing additional prediction errors when the receiver is in motion.

[0104] When calculating the likelihood function, the present invention adds a hidden variable Z to indicate whether the satellite has a fault and its probability. The joint log-likelihood function of (Y, Z) is converted to

[0105]

[0106] Take the conditional expectation of the posterior distribution of Z and dynamically adjust the log-likelihood function logp(y) of each satellite at each iteration according to the posterior distribution of Z. i |Z i=k,θ), which not only utilizes the information of all satellites but also assigns different weights to the log-likelihood function of each satellite according to the probability of failure of each satellite. It can adaptively calculate the size of its deviation without artificial assumptions and has strong real-time performance.

[0107] Step 4. Set the initial value of the iteration.

[0108] where β (0) It can be estimated based on the state equation of motion or state prediction in Kalman filtering, or it can be determined using robust positioning methods.

[0109] The range of satellite observation value deviation can be selected based on multiple experiments and experience. For example, it can be set in pseudo-range positioning. The unit is meter.

[0110] The standard deviation of the normal distribution of each satellite's observation value under each fault mode can generally be selected based on empirical values ​​through multiple experiments to make the The bias makes the distribution flatter. For example, when pseudo-range positioning, you can set The unit is square meters.

[0111] and γ (0) is the uninformative prior probability that the observation value of each satellite may fail. Generally, it can be selected based on the probability of failure in multiple experiments according to the empirical value. For example, γ (0) =(0.6,0.2,0.2) T .

[0112] Step 5. Based on θ (t-1) calculate t=1,2,.... t represents iterative update θ=(β T ,μ T ,(σ 2 ) T ,γ T ) T Each time θ is updated, the log-likelihood function Q(θ|θ (t-1) ,Y) value increases. represents the updated value of the posterior probability of the tth iteration based on the observed data that the ith satellite is in the kth failure mode.

[0113] It can be calculated

[0114] The difference between the present invention and other robust positioning methods is that other methods fail to consider the probability of failure of each satellite and the size of the failure deviation, directly eliminate the failed satellites, fail to utilize all information, and only solve parameters such as position based on satellites that have not failed. Due to the increase in the number of satellites involved in positioning, or directly estimate the size of the deviation, it is not based on actual observation data and may introduce additional prediction errors.

[0115] The method of the present invention is based on actual observation data, adaptively calculates the posterior probability of each satellite failure and the range of its failure or multipath deviation, uses the posterior probability and the deviation size for equivalent weighting, and utilizes all observation data, thereby achieving good positioning accuracy while taking into account robustness.

[0116] Step 6. Calculate γ (t) .Will Substitute Q(θ|θ (t-1) ,Y), find Let Q(θ|θ (t-1) ,Y) maximization, which is equivalent to maximizing

[0117]

[0118] Use the Lagrange multiplier method to add constraints Available

[0119] Step 7. Calculate (β, μ, σ 2 ) (t) .Will Substitute Q(θ|θ (t-1) ,Y), find Let Q(θ|θ (t-1) ,Y) maximization is carried out in two steps.

[0120] Step 1: Fix (σ 2 ) (t) , find the value that makes Q(θ|θ (t-1) ,Y) reaches the maximum

[0121]

[0122]

[0123] e n =(1,...,1) T

[0124]

[0125]

[0126] Among them use or The more recent of the two values.

[0127] Step 2: In Under the condition of (t-1) ,Y) reaches the maximum (σ 2 ) (t) .

[0128]

[0129] Then, update β alternately (t) and (σ 2 ) (t) , until β (t) and (σ 2 ) (t) The difference between the updated value and the value before the update is less than the given first threshold value, completing β (t) and (σ 2 ) (t) It can be proved that this iterative calculation method can make β (t) and σ (t) Converges to a stable extreme point.

[0130] Here, the first threshold value can be set manually based on multiple experiments, for example, β and σ calculated based on multiple experiments 2 The average error is determined by 5%.

[0131] The weight setting of the present invention is different from other robust positioning methods in that other methods either directly eliminate the faulty satellite, or set the weight of the faulty satellite based on information irrelevant to the current observation, such as the carrier-to-noise ratio, or judge the faulty satellite based on the least squares positioning result that contains fault information and is no longer reliable.

[0132] From this step It can be seen that the present invention is actually an iterative dynamic adaptive weighted least squares positioning. The setting of weights takes into account the probability of failure of each satellite and the standard error of the failure. The weight setting is updated at each iteration and is calculated based on the observed data, which avoids the adverse effects of inaccurate human assumptions on the solution.

[0133] At the same time, from W k It can be seen from the calculation that the present invention takes into account the influence of the faulty satellite on positioning while positioning, so it can eliminate the adverse influence of the faulty satellite while positioning, and there is no need to remove the faulty satellite and position again, which can reduce time consumption.

[0134] Step 8. Repeat steps 5-7, each time updating t to t+1 before starting step 5, until the θ calculated in step 7 is (t) and θ (t-1) Until the difference is less than the given first threshold value. At this time, t is the final number of iterative update calculations, θ (t) β in (t) That is the value of the unknown parameter vector β of this positioning cycle.

[0135] According to the theoretical basis of the EM algorithm, the above EM algorithm can be used to iteratively calculate θ (t) The sequence converges to a stable extreme point. Experiments show that the method of the present invention converges quickly, typically completing 2-3 iterations to reach a stable value. The EM algorithm is used for parameter solution. By calculating the expected log-likelihood function and deriving formulas, the positioning solution is transformed into an adaptive weighted least squares solution based on satellite failure probability and failure error, eliminating the need for manual weight setting.

[0136] Step 9. Calculate the real-time protection level HPL and VPL of this positioning, where HPL is the horizontal protection level and VPL is the vertical protection level. Compare them with the given alarm threshold (second threshold value) to determine whether this positioning is successful and output the final status.

[0137] Calculate the real-time protection level HPL and VPL of this positioning, that is, under a given false alarm probability, the position and other parameters The maximum possible error.

[0138]

[0139] A E =(H E T W E H E ) -1 H E T W E

[0140]

[0141] The posterior covariance matrix Σ of Y Y , is a diagonal matrix, and the i-th element in the diagonal element is

[0142] where μ k , These are the optimal values ​​of the parameter vectors to be solved in step 8 The value of the corresponding element in , that is, the final value when the iteration stops, γ ikAccording to γ (t) Calculate the value of the corresponding element by step (33)

[0143] The first p diagonal elements of The variance of each element (respectively C 11 ,...,C pp It can be used to calculate the horizontal protective level HPL and vertical protective level VPL, etc.

[0144]

[0145]

[0146] κ α It can be set based on experience or the false alarm and missed alarm probabilities required in practice. For example, it can generally be set to 5.

[0147] Compare VPL and HPL with the given alarm threshold (second threshold value). If VPL and HPL are lower than the given alarm threshold (second threshold), it indicates that the positioning is successful. If they are higher than the second threshold, it indicates that the positioning result is unreliable.

[0148] The second threshold value can be set based on experience or the accuracy required by the specific application scenario. For example, the pseudorange positioning can be set to 15 meters. Alternatively, two different threshold values ​​can be set for the horizontal and vertical directions of VPL and HPL, respectively. In this case, the second threshold value is a vector of length 2.

[0149] The difference between the present invention and other robust positioning methods is that other methods only consider the standard error of the fault-free satellites. When the presence of a faulty satellite or the observations contain multipath errors, The value of will be affected and increase, and the present invention adds two parameters during positioning Used to estimate the standard error of faulty satellites or multipath effects, which makes The estimated value of is more accurate, which also makes the positioning result more accurate. It absorbs satellite failures or multipath errors, making the positioning results more robust.

[0150] Step 10: Determine whether each satellite has a fault that affects positioning.

[0151] calculate The test statistic is:

[0152] Represents the median of all residual absolute values. The rejection region of this hypothesis test is {T i ||T i|>T}. |T i |>T indicates that there is a fault or deviation in the i-th satellite, otherwise it means that there is no fault or deviation in the observation value of the i-th satellite.

[0153] T is the third threshold value, which can be set to T = Φ -1 (1-α / 2). T is the α quantile of the standard normal distribution, and α is the significance level determined after comprehensively considering the false alarm and missed alarm rates. Φ -1 It is the inverse function of the standard normal cumulative distribution function. α can generally be selected based on experience or the tolerance of failure, such as 5%.

[0154] The difference between the present invention and other methods lies in that the present invention combines fault and multipath deviation identification with positioning, reduces the influence of faulty satellites on positioning through adaptive weights based on deviation probability and deviation value, and can perform fault detection and elimination at the same time as positioning and solving, without the need for re-positioning and solving. The calculation steps are few, avoiding the need for other methods to judge and eliminate faulty satellites after positioning and then re-position and solve, reducing the amount of calculation, speeding up the operation speed, and improving positioning accuracy while taking into account robustness.

[0155] At the same time, the judgment of the faulty satellite by the present invention is different from the two states of presence or absence of other methods. Instead, it is reflected by the probability of fault occurrence and the size of the deviation. It can flexibly judge whether there is a fault under different thresholds, provide richer information, and also feedback more effective information to other receiver modules.

[0156] Based on existing theories, the present invention demonstrates through experiments that the receiver fault detection and robust positioning method against multipath based on Gaussian mixture distribution can quickly identify faulty satellites and adaptively determine the impact of the faulty satellite on positioning. This method is a simple, efficient and fast faulty satellite identification method for user-side receivers, and provides reliability assurance for the use of GPS, Beidou and dual-mode joint positioning, timing and navigation.

[0157] The present invention has been described in detail above with reference to specific embodiments and exemplary examples. However, these descriptions should not be construed as limiting the present invention. Those skilled in the art will appreciate that various equivalent substitutions, modifications, or improvements may be made to the technical solutions and implementations of the present invention without departing from the spirit and scope of the present invention, all of which fall within the scope of the present invention. The scope of protection of the present invention shall be determined by the appended claims.

[0158] The contents not described in detail in the specification of the present invention belong to the common knowledge of those skilled in the art.

Claims

1. An adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution, characterized in that: The following steps are involved: (1) Establish the positioning solution model of the receiver within the positioning period and convert it linearly to obtain the standardized positioning solution model Y = Hβ + ε, ε ~ N (0, σ 2 I), where Y=(y1,...,y n ) T is the standardized transformation of the known satellite observation matrix, H=(h1,...,h n ) T is the normalized transformation of the n×p-dimensional known observation positioning geometry matrix, n is the total number of satellites, β is the p×1-dimensional unknown parameter vector, and ε is the normalized transformation of the n×1-dimensional observation noise matrix, which obeys the normal distribution N(0,σ 2 I),σ 2 is the variance coefficient, I is the unit diagonal matrix; (2) Use Gaussian mixture distribution to model the satellite observations and define the n-dimensional unknown vector Z as a latent variable, where each element Z i ∈{1,2,3},i=1,2,...,n, and obtain the Gaussian mixture distribution model of each satellite observation where γ k Z i The prior probability of make represents the parameter vector to be solved consisting of all unknown parameters to be solved in the Gaussian mixture distribution model of each satellite observation, where μ=(μ1,...,μ K ) T ,μ1=0,σ 2 =(s 2 1,...,s 2 K ) T ,γ=(γ1,...,γ K ) T ,K=3,γ3=1-γ1-γ2; (3) Iteratively solve the parameter vector θ based on the EM algorithm to obtain the optimal value of the parameter vector (4) According to the optimal value of the parameter vector to be solved Calculate the real-time protection level of positioning and compare it with the given alarm threshold to determine whether the positioning result is reliable; (5) The optimal value of the parameter vector to be solved in Substitute the obtained standardized positioning solution model in step (1) to obtain the test statistic. Based on the comparison between the test statistic and the rejection region of the hypothesis test, it is determined whether each satellite has a fault that affects positioning.

2. The method for adaptive real-time multipath elimination and robust positioning based on non-Gaussian distribution according to claim 1, characterized in that: In step (2), the Gaussian mixture distribution model of each satellite observation is established by the following method: Define n-dimensional unknown vector Z as hidden variable, where each element Z i ∈{1,2,3},i=1,2,...,n,Z i =1 indicates observation noise ε i Does not contain deviations that cause unacceptable errors in receiver positioning, Z i =2 or Z i =3 respectively represent the observation noise ε i Contains varying degrees of deviation that can cause unacceptable errors in receiver positioning; Assuming that different Z i Under the value, the satellite observation quantity y i Obey different Gaussian distributions, namely where y i is the i-th element of the known satellite observation matrix Y, h i The transpose of the i-th row of the geometry matrix H for the known observation location, h i T Indicates h i The transpose of f k (y i ) is a given Z i Time i The conditional Gaussian distribution of Z i Unknown time i The mixed Gaussian density function f m (y i ) in the k-th Gaussian component; When Z i When unknown, the satellite observation quantity y i For Gaussian mixture distribution, the Gaussian mixture distribution model is as follows: Among them, γ k Z i The prior probability, P(Z i =k) ​​=γ k .

3. The method for adaptive real-time multipath elimination and robust positioning based on non-Gaussian distribution according to claim 1, characterized in that: In step (3), the parameter vector θ to be solved is iteratively solved based on the EM algorithm to obtain the optimal value of the parameter vector to be solved The method comprises the following steps: (31) Calculate the expected log-likelihood function of the known satellite observation matrix Y: in t is the number of iterations of the parameter vector θ to be solved; (32) Set the initial value θ of the parameter vector θ to be solved (0) , at this time, take t = 0; (33) According to θ (t) Substitute in and use the Bayesian method to calculate (34) Update t to t+1 and change the value calculated in step (33) to Recorded as Substitute the expected log-likelihood function Q(θ|θ (t-1) ,Y), solve Let Q(θ|θ (t-1) ,Y) is maximized, and we get (35) As γ k , substitute the expected log-likelihood function Q(θ|θ (t-1) ,Y), solve Let Q(θ|θ (t-1) ,Y)About Maximization; (36) Repeat steps (33), (34) and (35), each time updating t to t+1 before starting step (34), until the θ calculated in step (35) is (t) and θ (t-1) The difference is less than the given first threshold value, and the final iteration number t is obtained. (t) , the optimal value of the parameter vector to be solved 4. The method for adaptive real-time multipath elimination and robust positioning based on non-Gaussian distribution according to claim 3, characterized in that: In the step (32), the initial value of the iteration of the parameter vector θ to be solved is where β (0) Estimated from the state equation or state prediction in Kalman filtering, or determined using robust positioning methods, μ (0) ,(σ 2 ) (0) and γ (0) Determined based on historical experience data.

5. The method for adaptive real-time multipath elimination and robust positioning based on non-Gaussian distribution according to claim 3, characterized in that: In step (35), solve Let Q(θ|θ (t-1) ,Y) maximization steps are as follows: (351)Fixed(σ 2 ) (t) , find the value that makes Q(θ|θ (t-1) ,Y) reaches the maximum e n ＝(1,...,1) T Among them use or The latest value of the two, steps (351) and (352) need to be repeated multiple times, each time updating and Until the difference between the two updated values ​​is less than a given first threshold; (352) Under the condition of (t-1) ,Y) reaches the maximum (σ 2 ) (t) ; 6. The method for adaptive real-time multipath elimination and robust positioning based on non-Gaussian distribution according to claim 3, characterized in that: In step (36), the first threshold value is determined as β and σ calculated by multiple experiments. 2 The average error is 5%.

7. The method for adaptive real-time multipath elimination and robust positioning based on non-Gaussian distribution according to claim 1, characterized in that: In the step (4), the real-time protective level of the positioning includes the horizontal protective level HPL and the vertical protective level VPL; Among them, κ α is the variance expansion coefficient, which is set according to the empirical value or the actual false alarm and missed alarm probability. 11 , C 22 and C 33 are the optimal values ​​of the parameter vectors to be solved middle The variance of the 1st, 2nd, and 3rd elements of ; The alarm threshold is a vector with a length of 2, including threshold values ​​in the horizontal direction and the vertical direction.

8. The method for adaptive real-time multipath elimination and robust positioning based on non-Gaussian distribution according to claim 7, characterized in that: described The variance of each element of is the matrix The first p diagonal elements in ; Among them A E =(H E T W E H E ) -1 H E T W E , Σ Y is the posterior covariance matrix of the known satellite observation matrix Y, which is a diagonal matrix. The i-th element in the diagonal element is: in are the optimal values ​​of the parameter vectors to be solved in step (3) The value of the corresponding element in ,γ ik According to γ (t) The value of the corresponding element is calculated by step (33).

9. The method for adaptive real-time multipath elimination and robust positioning based on non-Gaussian distribution according to claim 1, characterized in that: In step (5), the optimal value of the parameter vector to be solved is in Substituting into the standardized positioning solution model obtained in step (1), the test statistic obtained is: in represents the median of the absolute values ​​of all residuals, 10. The method for adaptive real-time multipath elimination and robust positioning based on non-Gaussian distribution according to claim 1, characterized in that: In step (5), the rejection region of the hypothesis test is {T i ||T i |>T},|T i |>T indicates that the i-th satellite has a fault that affects positioning; otherwise, it indicates that the i-th satellite does not have a fault that affects positioning; The T is the third threshold value, T=Φ -1 (1-α / 2), T is the α quantile of the standard normal distribution, α is the significance level determined after comprehensively considering the false alarm and missed alarm rates, Φ -1 is the inverse of the standard normal cumulative distribution function.

11. The method for adaptive real-time multipath elimination and robust positioning based on non-Gaussian distribution according to claim 1, characterized in that: In step (1), the positioning solution model of the receiver during the positioning period is in is the n×1-dimensional known observation quantity, n is the total number of satellites, is the n×p dimensional known observation positioning geometry matrix, β is the p×1 dimensional unknown parameter vector, is the observation error vector, which is an n×1 dimensional vector, Obey the normal distribution N(0,σ 2 ∑),σ 2 is the variance coefficient, ∑ is the covariance prior matrix of each satellite observation error, The linearization conversion process is as follows: After linearization conversion I is the unit diagonal matrix; Get Y=Hβ+ε,ε~N(0,σ 2 I); in

Images