CN113848570A - Self-adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution - Google Patents

Self-adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution Download PDF

Info

Publication number
CN113848570A
CN113848570A CN202110988347.1A CN202110988347A CN113848570A CN 113848570 A CN113848570 A CN 113848570A CN 202110988347 A CN202110988347 A CN 202110988347A CN 113848570 A CN113848570 A CN 113848570A
Authority
CN
China
Prior art keywords
positioning
satellite
theta
gaussian
value
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202110988347.1A
Other languages
Chinese (zh)
Inventor
由立华
毕波
陈启亮
高参
吴雨航
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Microelectronic Technology Institute
Mxtronics Corp
Original Assignee
Beijing Microelectronic Technology Institute
Mxtronics Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Microelectronic Technology Institute, Mxtronics Corp filed Critical Beijing Microelectronic Technology Institute
Priority to CN202110988347.1A priority Critical patent/CN113848570A/en
Publication of CN113848570A publication Critical patent/CN113848570A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/20—Integrity monitoring, fault detection or fault isolation of space segment
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/22—Multipath-related issues
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/23—Testing, monitoring, correcting or calibrating of receiver elements
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Security & Cryptography (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

The invention discloses a self-adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution, which does not make assumptions on the number of satellites with faults or deviations, but modifies the Gaussian distribution assumption of errors in the traditional positioning model to change the Gaussian distribution assumption into Gaussian mixed distribution, namely non-Gaussian distribution, solves parameters by a maximum likelihood method based on real-time observation data, calculates the probability distribution condition of the deviation range of each satellite fault in a self-adaptive manner, can be used for identification and robust positioning under multiple faults and robust positioning and multipath elimination under a multipath environment, overcomes the defects of the existing RAIM and robust positioning methods, has larger calculation amount, depends on the prior assumption of the number of the faults, has higher misjudgment rate, and has higher fault detection success rate, meanwhile, the method has small operand, high speed, high precision and robustness.

Description

Self-adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution
Technical Field
The invention belongs to the technical field of satellite navigation positioning, and relates to a non-Gaussian distribution based self-adaptive real-time multipath elimination and robust positioning method.
Background
When a satellite navigation system is used for positioning a receiver, measured values such as pseudoranges or carriers observed by the receiver may contain propagation errors, such as electronic faults, ephemeris and clock errors broadcast by satellites, abnormal atmospheric delay, multipath effects, and faults of the receiver, which may further cause a large error in a positioning result output by a user receiver, and further may adversely affect the safety of the user itself. With the rapid development and wide use of the global satellite navigation system, real-time and rapid integrity fault monitoring and robust positioning methods play an increasingly important role.
RAIM (Receiver Autonomous Integrity Monitoring) is that a Receiver detects observation values of each satellite in real time, and detects and identifies a fault by using redundancy formed by a plurality of satellite observation quantities received by the Receiver. The method has the advantages of no need of assistance of external equipment, low cost and easy realization, but has the defects that only the fault satellite can be removed after being identified, the operation time is increased, and the positioning precision is reduced due to the reduction of the number of the satellites.
At present, the widely used RAIM methods mainly include two major types, namely a residual error analysis method and a maximum solution separation method. The residual error analysis method is based on the least square residual error, after least square positioning, the residual error square sum or the student residual error is calculated and compared with a threshold value, however, the least square method is very sensitive to a fault star, when the fault star has a large error, a large deviation occurs in a positioning result, the estimated absolute value of the residual error of the fault star is reduced, meanwhile, the absolute value of the residual error of a non-fault star is increased, finally, misjudgment is caused, and the real fault star cannot be accurately identified. The maximum solution separation rule divides all the observed values into a plurality of groups of subset combinations according to a certain method, positions each group of observed subsets respectively, and compares and judges the observed subsets based on some statistics such as the square sum of residual errors after positioning, and the like, for example, the optimal combination is a combination without fault satellites, and the like.
Compared with the RAIM method, the robust positioning method does not eliminate the fault satellite, is a stable positioning method, overcomes the defect that the sum of squared residuals is very sensitive to large errors in the least square method, and reduces the influence of fault observed quantity on the positioning result while the receiver positions by utilizing all satellite observed quantities by constructing a stable loss function. However, the error distribution assumed by the robust positioning method is thicker than the normal distribution, so that when there is no satellite fault, the positioning result obtained by the robust positioning method has larger error than that obtained by the least square method, and is not always an unbiased estimation.
When satellite signals received by a receiver are shielded or reflected, the observed values contain large errors, and the multipath errors can be regarded as special fault errors, namely, the fault errors occur in a plurality of satellite observed values at the same time in the same observation moment, and compared with other fault modes, the number of faults is more. Most of the existing RAIM fault detection methods aim at the fault number not exceeding three or four, but in the multi-mode multi-frequency point positioning process, the available star number is large, and the multi-path fault number far exceeds three or four. The identification and elimination of the multipath error are key factors affecting the positioning error of the receiver, and a fault elimination positioning method aiming at more faults is needed.
Disclosure of Invention
The invention aims to overcome the defects, provides a non-Gaussian distribution based adaptive real-time multipath elimination and robust positioning method, does not make assumptions on the number of failed or deviated satellites, but modifies the Gaussian distribution assumption of errors in the traditional positioning model to change the Gaussian distribution into Gaussian mixed distribution, namely non-Gaussian distribution, solves parameters by a maximum likelihood method based on real-time observation data, adaptively calculates the probability distribution situation of the deviation range of each satellite fault, can be used for identification and robust positioning under multiple faults and robust positioning and multipath elimination under a multipath environment, overcomes the defects of the existing RAIM and robust positioning method, has larger calculation amount depending on the prior assumption of the number of the faults and higher misjudgment rate, and has higher fault detection success rate, meanwhile, the method has small operand, high speed, high precision and robustness.
In order to achieve the above purpose, the invention provides the following technical scheme:
a non-Gaussian distribution based self-adaptive real-time multipath elimination and robust positioning method comprises the following steps:
(1) establishing a positioning calculation model of the receiver in the positioning period, and linearly converting the positioning calculation model to obtain a standardized positioning calculation model Y which is H beta + epsilon, epsilon-N (0, sigma)2I) Wherein Y is (Y)1,...,yn)TFor a known normalized transformation of the observation matrix of each satellite, H ═ H1,...,hn)TStandardized transformation of a known observation positioning geometric matrix for N x p dimensions, N being the total number of satellites, beta being an unknown parameter vector for p x 1 dimensions, and epsilon being the standardized transformation of an observation noise matrix for N x 1 dimensions, obeying a normal distribution N (0, sigma)2I),σ2Is a variance coefficient, and I is a unit diagonal matrix;
(2) modeling each satellite observation using a mixture of Gaussian distributions, defining an n-dimensional unknown vector Z as a hidden variable, wherein each element Z is a variableiE {1,2,3}, i ═ 1, 2.. times, n, a gaussian mixture distribution model of each satellite observation is obtained
Figure BDA0003231578270000031
Wherein gamma iskIs ZiThe prior probability of (a) being,
Figure BDA0003231578270000032
order to
Figure BDA0003231578270000033
A parameter vector to be solved consisting of all unknown parameters to be solved in a Gaussian mixture distribution model representing the observed quantity of each satellite, wherein
μ＝(μ1,...,μK)T,μ1＝0,σ2＝(σ2 1,...,σ2 K)T,γ＝(γ1,...,γK)T,K＝3,γ3＝1-γ1-γ2。
(3) Iteratively solving a parameter vector theta to be solved based on an EM algorithm to obtain an optimal value of the parameter vector to be solved;
Figure BDA0003231578270000034
(4) according to the optimal value of the parameter vector to be solved
Figure BDA0003231578270000035
Calculating the real-time protective level of positioning, comparing the real-time protective level with a given alarm threshold, and judging whether the positioning result is reliable or not;
(5) the optimal value of the parameter vector to be solved
Figure BDA0003231578270000036
In (1)
Figure BDA0003231578270000037
And (3) substituting the standard positioning calculation model obtained in the step (1) to obtain test statistic, and judging whether each satellite has a fault affecting positioning or not according to comparison between the test statistic and the rejection region of the hypothesis test.
Further, in the step (2), a gaussian mixture distribution model of each satellite observation is established by the following method:
defining an n-dimensional unknown vector Z as a hidden variable, wherein each element Zi∈{1,2,3},i＝1,2,...,n，Zi1 denotes the observation noise εiWithout deviations, Z, causing unacceptable errors in the positioning of the receiveri2 or Zi3 denotes the observation noise εiContain varying degrees of bias that cause unacceptable errors in receiver positioning;
assuming at different ZiUnder the value, each satellite observed quantity yiSubject to different Gaussian distributions, i.e.
Figure BDA0003231578270000041
Figure BDA0003231578270000042
Wherein y isiFor the i-th element, h, of the known satellite observation matrix YiTransposing the ith row of the geometric matrix H for known observationsi TRepresents hiTransposing; f. ofk(yi) For a given ZiTime yiConditional Gaussian distribution of (a) also being ZiUnknown time yiMixed gaussian density function fm(yi) The kth gaussian component of (1);
when Z isiAt unknown time, each satellite observed quantity yiFor gaussian mixture distribution, the gaussian mixture distribution model is obtained as follows:
Figure BDA0003231578270000043
wherein, γkIs ZiA priori of, P (Z)i＝k)＝γk。
Further, in the step (3), the parameter vector θ to be solved is iteratively solved based on the EM algorithm to obtain the parameter to be solvedOptimal value of number vector
Figure BDA0003231578270000044
The method comprises the following steps:
(31) and calculating to obtain an expected log-likelihood function of each known satellite observation matrix Y:
Figure BDA0003231578270000045
Figure BDA0003231578270000051
wherein
Figure BDA0003231578270000052
t is the iteration number of the parameter vector theta to be solved;
(32) setting an iteration initial value theta of a parameter vector theta to be solved(0)Taking t as 0;
(33) according to theta(t)Substituting, and calculating by Bayesian method
Figure BDA0003231578270000053
(34) Updating t to t +1, calculating in step (33)
Figure BDA0003231578270000054
Substituting into the desired log-likelihood function Q (θ | θ)(t-1)Y), solving for
Figure BDA0003231578270000055
Make Q (theta | theta)(t-1)Y) maximization to obtain
Figure BDA0003231578270000056
(35) Will be provided with
Figure BDA0003231578270000057
As gammakValue substitution period ofThe log-likelihood function Q (theta | theta)(t-1)Y), solving for
Figure BDA0003231578270000058
Make Q (theta | theta)(t-1)Y) maximization;
(36) repeating steps (33), (34) and (35), each time updating t to t +1 before starting step (34), until θ calculated in step (35)(t)And theta(t-1)Is smaller than a given first threshold value, theta is obtained when the final iteration number t is obtained(t)Optimum value of parameter vector to be solved
Figure BDA0003231578270000059
Further, in the step (32), the iterative initial value of the parameter vector θ to be solved
Figure BDA00032315782700000510
Wherein beta is(0)Estimated from equations of motion or state prediction in Kalman filtering, or determined using robust positioning methods, mu(0)，(σ2)(0)And gamma(0)And determined from historical empirical data.
Further, in the step (35), solving
Figure BDA00032315782700000511
Make Q (theta | theta)(t-1)Y) the step of maximization is as follows:
(351) fixation (sigma)2)(t)Then, Q (theta | theta) is obtained(t-1)Y) maximization
Figure BDA0003231578270000061
Figure BDA0003231578270000062
Figure BDA0003231578270000063
en＝(1,...,1)T
Figure BDA0003231578270000064
Figure BDA0003231578270000065
Figure BDA0003231578270000066
Therein
Figure BDA0003231578270000067
Use of
Figure BDA0003231578270000068
Or
Figure BDA0003231578270000069
The latest value of the two, step 351 and step 352, needs to be repeated for a plurality of times, and each time the value is updated
Figure BDA00032315782700000610
And
Figure BDA00032315782700000611
until the difference between the two update values is less than a given first threshold value;
(352) in that
Figure BDA00032315782700000612
Under the condition (2), Q (theta | theta) is obtained(t-1)Y) maximization2)(t)；
Figure BDA00032315782700000613
Further, in the step (36), the first thresholdThe values were determined as beta and sigma calculated for multiple experiments25% of the average error.
Further, in the step (4), locating the real-time protective level comprises a horizontal protective level HPL and a vertical protective level VPL;
Figure BDA00032315782700000614
wherein, καSetting the variance expansion coefficient according to empirical value or actual false alarm and false alarm probability, C11，C22And C33Respectively are the optimal values of the parameter vectors to be solved
Figure BDA0003231578270000071
In
Figure BDA0003231578270000072
Variance of 1,2,3 elements of (a);
the alarm threshold is a vector with the length of 2 and comprises threshold values in the horizontal direction and the vertical direction.
Further, the
Figure BDA0003231578270000073
Is a matrix of variances of each element
Figure BDA0003231578270000074
The first p diagonal elements in (1);
Figure BDA0003231578270000075
wherein A isE＝(HE TWEHE)-1HE TWE，
Figure BDA0003231578270000076
ΣYThe posterior covariance matrix of each known satellite observation matrix Y is a diagonal matrix, and the ith element in the diagonal element is:
Figure BDA0003231578270000077
wherein muk,
Figure BDA0003231578270000078
Are all the optimal values of the parameter vectors to be solved in the step (3)
Figure BDA0003231578270000079
Value of the corresponding element in (1), γikAccording to gamma(t)The value of the corresponding element is calculated by step (33).
Further, in the step (5), the optimal value of the parameter vector to be solved is obtained
Figure BDA00032315782700000710
In (1)
Figure BDA00032315782700000711
Substituting into the standardized positioning calculation model obtained in the step (1), and calculating
Figure BDA00032315782700000712
The test statistics obtained were:
Figure BDA00032315782700000713
wherein
Figure BDA00032315782700000714
Representing the median of the absolute values of all residuals (i.e. estimates of the errors).
Further, in the step (5), the rejection field of the hypothesis test is { T }i||Ti|>T}，|Ti|>T shows that the ith satellite has a fault which influences the positioning, otherwise shows that the ith satellite does not have the fault which influences the positioning;
t is a third threshold value, and T is phi-1(1-alpha/2), T is an alpha quantile point of standard normal distribution, and alpha is a significance level determined after comprehensively considering false alarm rate and false alarm rate, for example, the significance level can be 5 percent and phi-1Is the inverse function of a standard normal cumulative distribution function.
Further, in the step (1), the positioning calculation model of the receiver in the positioning period is
Figure BDA00032315782700000715
Wherein
Figure BDA00032315782700000716
Known observations are of dimension n x 1, n is the total number of satellites,
Figure BDA00032315782700000717
a geometric matrix is positioned for known observation of dimension n x p, beta is an unknown parameter vector of dimension p x 1,
Figure BDA0003231578270000081
an observation error vector is an n x 1-dimensional vector,
Figure BDA0003231578270000086
obeying a normal distribution N (0, σ)2∑),σ2Is a variance coefficient, sigma is a covariance prior matrix of the observed quantity error of each satellite,
Figure BDA0003231578270000082
the linearization conversion process is as follows:
Figure BDA0003231578270000083
after linear conversion
Figure BDA0003231578270000084
I is a unit diagonal matrix;
obtaining Y ═ H beta + epsilon, epsilon-N (0, sigma)2I)；
Wherein
Figure BDA0003231578270000085
Compared with the prior art, the invention has the following beneficial effects:
(1) the invention relates to a non-Gaussian distribution based self-adaptive real-time multipath elimination and robust positioning method, which combines RAIM, multipath elimination and positioning together, reduces the influence of faults or satellite observation values containing deviation on positioning through self-adaptive weight based on fault probability and fault error magnitude, can detect, identify and eliminate faults while resolving the positioning, and has less operation steps and strong real-time performance;
(2) the self-adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution is more stable in self-adaptive weight during positioning than other weighted least square solving methods. Other weighted least square methods are very sensitive to weight setting, when the weight setting is slightly deviated from the actual value, the positioning result is greatly deviated, so that the fault satellite cannot be successfully identified, the weight calculation of the method is self-adaptive iterative calculation based on real-time observation data, and the weight of the fault satellite can be self-adaptively adjusted according to the fault deviation, so that good positioning precision, success rate and robustness of fault detection can be ensured;
(3) the invention relates to a non-Gaussian distribution based self-adaptive real-time multipath elimination and robust positioning method, which can give the probability of each satellite failing and the size of a failure range based on observation data instead of only two states of yes and no in the identification of a failed satellite.
(4) The invention relates to a non-Gaussian distribution based self-adaptive real-time multipath elimination and robust positioning method, which does not need to presume the number of failed satellites in advance. The other methods need to respectively carry out observation subset calculation and comparison under the assumption of different numbers of fault satellites, the calculation times are multiple, especially when the number of fault or deviation satellites is large, the time consumption is long, the method modifies the Gaussian distribution assumption of errors in the traditional positioning model, changes the Gaussian distribution assumption into Gaussian mixed distribution, directly estimates the probability of each satellite having faults, does not need to remove the fault satellites, and is high in calculation speed.
Detailed Description
The features and advantages of the present invention will become more apparent and appreciated from the following detailed description of the invention.
The word "exemplary" is used exclusively herein to mean "serving as an example, embodiment, or illustration. Any embodiment described herein as "exemplary" is not necessarily to be construed as preferred or advantageous over other embodiments.
Aiming at the defects of the existing RAIM and robust positioning method and the like, the invention provides a non-Gaussian distribution based self-adaptive real-time multipath elimination and robust positioning method, which can perform multipath elimination and observed quantity multi-fault detection while positioning, avoid the problem of re-positioning after fault detection after positioning by other methods, and greatly reduce the operation time. Secondly, the state estimation of whether each satellite has faults or not, the fault error magnitude and the fault probability is added, the iterative estimation of the posterior probability is carried out based on observation data, the solution is carried out by using an adaptive weighted least square method based on the fault probability and the fault error, and compared with other least square positioning methods based on fixed weight, the weight of the satellite can be automatically adjusted by the method according to the fault probability and the error magnitude of the satellite, the method has good adaptivity to the number and the error magnitude of fault satellites, the influence of the fault satellite on the positioning result is greatly reduced while the precision is ensured, and the robustness of the positioning result is considered. In addition, the method avoids the situation that other methods need to make assumptions on the number of the failed satellites and perform positioning calculation and comparison under different assumptions, solves the optimal expected log maximum likelihood directly on the basis of the satellite failure probability during calculation, and greatly reduces the calculation amount. Finally, the method can also estimate the probability of each satellite failure, and overcomes the defect that other methods can only judge whether the satellite fails and cannot estimate the failure probability.
As shown in fig. 1, a non-gaussian distribution based adaptive real-time multipath cancellation and robust positioning method includes the following steps:
step 1, in each positioning period, a receiver positioning calculation model is established and standardized.
Equation of locationIs composed of
Figure BDA0003231578270000091
Wherein
Figure BDA0003231578270000092
For known observations, including pseudorange or carrier observations, an n x 1 dimensional vector, n being the total number of satellites,
Figure BDA0003231578270000101
for the known observation positioning geometrical matrix, the measurement is n multiplied by p dimension, beta is p multiplied by 1 dimension unknown parameter vector, namely the parameters to be solved of the receiver, including three-dimensional position coordinates, the parameters to be solved of the receiver clock error or speed,
Figure BDA0003231578270000102
an observation error vector is an n x 1-dimensional vector,
Figure BDA0003231578270000103
obeying a normal distribution N (0, σ)2∑),σ2Sigma is an unknown variance coefficient, represents a covariance prior matrix of observation errors of each satellite given in advance, generally, the observations of different satellites are assumed to be independent of each other, sigma is a positive definite diagonal matrix,
Figure BDA0003231578270000104
the value of Σ may be determined from a number of experiments in advance, for example, a function related to the satellite pitch angle, the carrier-to-noise ratio, and the like, or may be simply set as an identity matrix.
In order to facilitate subsequent calculation, applying linear transformation to the receiver positioning calculation model to standardize the receiver positioning calculation model to obtain a standardized receiver positioning calculation model:
Figure BDA0003231578270000105
Figure BDA0003231578270000106
after transformation
Figure BDA0003231578270000107
And I is a unit diagonal matrix.
Obtaining Y ═ H beta + epsilon, epsilon-N (0, sigma)2I)；
Wherein
Figure BDA0003231578270000108
Step 2, modeling each satellite observation quantity by using Gaussian mixture distribution, which is also equivalent to error modeling
Definition of εiIs any element of epsilon, i is more than or equal to 1 and less than or equal to n, and the observation noise epsilon of each satellite is assumediIndependent of each other. Since the observation noise epsilon can be influenced by multipath, signal shielding, faults of a satellite or a receiver, and the like, an n-dimensional unknown vector Z is defined, wherein each element Z is an implicit variablei∈{1,2,3},i＝1,2,...,n，ZiRepresents different failure modes, such as no failure or deviation, mean shift, variance increase, and more specifically, Zi1 denotes the observation noise εiWithout faults or deviations, Z, causing large errors in the positioning of the receiveri2 or Zi3 denotes the observation noise εiContain different degrees of bias that cause large errors in receiver positioning, where there may be a shift in the mean or an increase in the variance;
assuming at different ZiUnder the value, the observed quantity yiObey different gaussian distributions:
Figure BDA0003231578270000111
Figure BDA0003231578270000112
wherein y isiIs the step (1)) The ith element, h, of normalized known observed quantity YiTransposing the ith row of the geometric matrix H for known observations after normalizationi TRepresents hiThe transposing of (1).
When Z isiWhen unknown, yiFor gaussian mixture distribution, the observed quantity mixture gaussian distribution model, i.e. its density function, is as follows:
Figure BDA0003231578270000113
wherein, γkIs ZiA priori of, P (Z)i＝k)＝γk。
The parameter vector to be solved is
Figure BDA0003231578270000114
The latter 7 parameters are those of the gaussian mixture distribution. It should be noted that, let us assume that10, i.e. a gaussian distribution f with zero mean of the first component of the gaussian mixture and unknown error1(yi) The error distribution is error distribution without failure or deviation. Gamma ray3＝1-γ1-γ2Are redundant parameters.
Order to
Figure BDA0003231578270000115
And representing all vectors formed by unknown parameters to be solved, wherein the upper right corner mark T represents the transposition of the vectors.
μ＝(μ1,...,μK)T,μ1＝0,σ2＝(σ2 1,...,σ2 K)T,γ＝(γ1,...,γK)T,K＝3,γ3＝1-γ1-γ2。
When a plurality of multipath errors exist simultaneously, the invention uses non-Gaussian assumption for error distribution, namely, Gaussian mixed distribution is adopted to model the errors, so that the error distribution is better suitable for error distribution characteristics of a plurality of faults under the multipath effect. The difference from other positioning methods is that no prior assumptions need be made about the number of satellites that have failed or are biased. Other integrity fault elimination methods respectively perform positioning calculation under different numbers of fault satellites, and the calculation amount is large.
The invention does not make assumptions on the number of the failed or deviated satellites, but modifies the Gaussian distribution assumption of errors in the traditional positioning model to change the Gaussian distribution assumption into Gaussian mixed distribution, namely non-Gaussian distribution, solves parameters by a maximum likelihood method based on real-time observation data, and adaptively calculates the probability distribution condition of the deviation range of each satellite failure, thereby avoiding the reduction of resolving efficiency caused by improper artificial prior assumption, ensuring that the positioning result has good precision and robustness, and simultaneously eliminating the deviation caused by the influence of failure or multipath and the like on the satellite observation value.
Meanwhile, unknown parameters such as standard errors, deviation sizes and the like of the fault satellites are increased, the standard errors and the deviation sizes of the fault satellites and the fault satellites are respectively estimated while positioning calculation is carried out by using all the satellites, the influence of the fault satellites is reduced, and the positioning accuracy and the stability are improved.
And 3, calculating an expected log-likelihood function.
By expecting the posterior distribution about Z to obtain an expected log-likelihood function of the observed value, i.e.
Figure BDA0003231578270000121
Wherein
Figure BDA0003231578270000122
Will be calculated in step 5.
t is the number of iterations and the initial value of the iteration is theta(0)See step 4.
Figure BDA0003231578270000123
Indicating a failure mode ziK, a posteriori probability calculated at the t-th iteration based on real-time observations.
The difference between the method and other positioning methods is that other anti-difference positioning methods mainly calculate the log-likelihood of satellites without faults under the assumption of different numbers of fault satellites, obviously reduce the weight of the satellites with faults and cause less satellites to be positioned due to the fact that all information cannot be utilized. There are also multipath cancellation methods that estimate the multipath error, but do not utilize the real-time positioning matrix for estimation, and introduce additional prediction error when the receiver is in motion.
When the likelihood function is calculated, the invention adds the hidden variable Z to indicate whether the satellite has faults and the probability, and combines the combined log likelihood function of (Y, Z)
Figure BDA0003231578270000124
Conditional expectation is taken for Z posterior distribution, and the log likelihood function logP (y) of each satellite is dynamically adjusted in each iteration according to the Z posterior distributioni|ZiK, theta), different weights are given to the log-likelihood functions of the satellites according to the probability of failure of each satellite, and the deviation can be calculated in a self-adaptive manner without artificial assumption, so that the real-time performance is strong.
And 4, setting an iteration initial value.
Figure BDA0003231578270000131
Wherein beta is(0)It can be estimated from equations of motion or state prediction in kalman filtering, or it can be determined using robust positioning methods.
Figure BDA0003231578270000132
The range size for the satellite observation value to be deviated can be generally selected according to experience values through multiple experiments, such as pseudo-range positioning setting
Figure BDA0003231578270000133
The unit is meter.
Figure BDA0003231578270000134
The standard deviation of different normal distributions of the observed value of each satellite under each fault mode can be selected according to experience values through multiple experiments, and the standard deviation can be selected as much as possible
Figure BDA0003231578270000135
The larger the distribution, the more even the distribution. Such as pseudorange location may be set
Figure BDA0003231578270000136
The unit is square meter.
Figure BDA0003231578270000137
And is
Figure BDA00032315782700001310
γ(0)The prior probability of no information for possible failure of the observed value of each satellite can be selected empirically according to the probability of failure in multiple experiments, for example, γ can be set(0)＝(0.6,0.2,0.2)T。
Step 5. based on theta(t-1)Computing
Figure BDA0003231578270000138
t 1, 2. t denotes an iterative update θ ═ βT,μT,(σ2)T,γT)TThe number of times. Each time θ is updated, the log-likelihood function Q (θ | θ) is made(t-1)And Y) increases.
Figure BDA0003231578270000139
A posteriori probability update based on the t iteration of the observation data indicating that the ith satellite is in the kth failure mode.
Can be calculated to obtain
Figure BDA0003231578270000141
The method is different from other robust positioning methods in that other methods cannot consider the fault probability and fault deviation of each satellite, directly eliminate the faulted satellites, cannot utilize all information, and only solve parameters such as positions and the like based on the satellites which do not have faults, and because the number of the satellites participating in positioning is reduced, errors are increased, or the deviation of the satellites is directly estimated, and cannot be based on actual observation data, additional prediction errors are possibly introduced.
The method of the invention is based on actual observation data, adaptively calculates the posterior probability of each satellite fault and the fault or multipath deviation range, uses the posterior probability and the deviation for equivalent weighting, and uses all observation data, thereby obtaining good positioning precision and taking into account robustness.
Step 6. calculating gamma(t). Will be provided with
Figure BDA0003231578270000142
Substituting Q (θ | θ!)(t-1)Y), calculating
Figure BDA0003231578270000143
Make Q (theta | theta)(t-1)Y) maximization, equivalent to maximization
Figure BDA0003231578270000144
Figure BDA0003231578270000145
Using Lagrange multiplier method, adding constraint condition
Figure BDA0003231578270000146
Can obtain the product
Figure BDA0003231578270000147
Step 7. calculating (beta, mu, sigma)2)(t). Will be provided with
Figure BDA0003231578270000148
Substituting Q (θ | θ!)(t-1)Y), calculating
Figure BDA0003231578270000149
Make Q (theta | theta)(t-1)Y) maximization, in two steps.
The first step is as follows: fixation (sigma)2)(t)Then, Q (theta | theta) is obtained(t-1)Y) maximization
Figure BDA00032315782700001410
Figure BDA00032315782700001411
Figure BDA00032315782700001412
en＝(1,...,1)T
Figure BDA0003231578270000151
Figure BDA0003231578270000152
Figure BDA0003231578270000153
Therein
Figure BDA0003231578270000154
Use of
Figure BDA0003231578270000155
Or
Figure BDA0003231578270000156
The most recent of the two.
The second step is that: in that
Figure BDA0003231578270000157
Under the condition (2), Q (theta | theta) is obtained(t-1)Y) maximization2)(t)。
Figure BDA0003231578270000158
Then, updating beta alternately(t)And (σ)2)(t)Up to beta(t)And (σ)2)(t)Is less than a given first threshold value from the update value before updating, completing beta(t)And (σ)2)(t)And (4) solving. It can be shown that this iterative calculation method enables β(t)And σ(t)Converge to a stable extreme point.
The first threshold value can be chosen here by a person based on a number of experiments, for example, based on beta and sigma calculated from a number of experiments2Is determined as 5% of the average error.
The weight setting method is different from other robust positioning methods in that other methods either directly eliminate the fault satellite, or set the weight of the fault satellite according to information such as carrier-to-noise ratio and the like irrelevant to the observation, or judge the fault satellite based on a least square positioning result which contains fault information and is not reliable any more.
From this step
Figure BDA0003231578270000159
It can be seen that the invention is actually an iterative dynamic adaptive weighted least squares fix, the weight setting takes into account the probability of each satellite failing and the standard error of the failure, the weight setting is updated during each iteration and is calculated based on the observation data, thus avoiding the adverse effect of the artificial assumption inaccuracy on the solution.
At the same time, from WkThe calculation shows that the influence of the fault satellite on the positioning is considered while the positioning is carried out, so that the adverse influence of the fault satellite can be eliminated while the positioning is carried out, the fault satellite does not need to be removed for positioning again, and the time consumption can be reduced.
Step 8, repeating steps 5-7, each time updating t to t +1 before starting step 5, until theta calculated in step 7(t)And theta(t-1)Until the difference is smaller than a given first threshold value. T at this time is the final iterative update calculation number, θ(t)Beta in (A) to (B)(t)Namely the solved value of the unknown parameter vector beta of the positioning period,
Figure BDA0003231578270000166
from the theoretical basis of the EM algorithm, it can be seen that using the above-mentioned EM algorithm to iterate the calculation, θ will be made(t)The sequence converged to a stable extremum. Experiments show that the convergence rate of the method is high, and iteration convergence to a stable value can be completed for 2-3 times generally. And performing parameter solution by using an EM (effective velocity) algorithm, and converting positioning solution into self-adaptive weighted least square solution based on satellite fault probability and fault error by calculating an expected log-likelihood function and deducing a formula without manually setting weight.
And 9, calculating real-time protective horizontal HPL and VPL of the current positioning, wherein HPL is horizontal protective horizontal and VPL is vertical protective horizontal, comparing the HPL with a given alarm threshold (a second threshold), judging whether the current positioning is successful, and outputting a final state.
Calculating the real-time protection level HPL and VPL of the current positioning, namely the position and other parameters under the given false alarm probability
Figure BDA0003231578270000161
Is the maximum possible error.
Figure BDA0003231578270000162
AE＝(HE TWEHE)-1HE TWE
Figure BDA0003231578270000163
Posterior covariance matrix sigma of YYIs a diagonal matrix, the ith element in the diagonal element is
Figure BDA0003231578270000164
Wherein muk,
Figure BDA0003231578270000165
Are all the optimal values of the parameter vectors to be solved in the step 8
Figure BDA0003231578270000171
I.e. the final value at which the iteration stops, gammaikAccording to gamma(t)Calculating the value of the corresponding element by step (33)
Figure BDA0003231578270000172
The first p diagonal elements are
Figure BDA0003231578270000173
Of each element (respectively with C)11,...,CppRepresentative), which can be used to calculate the horizontal protective horizontal HPL and the vertical protective horizontal VPL, etc.
Figure BDA0003231578270000174
Figure BDA0003231578270000175
καCan be based on empirical values orAnd the actual needed false alarm, false alarm probability and the like, for example, the false alarm and false alarm probability can be generally set to 5.
And comparing the VPL and the HPL with a given alarm threshold (a second threshold), if the VPL and the HPL are lower than the given alarm threshold (the second threshold), the positioning is successful, and if the VPL and the HPL are higher than the second threshold, the positioning result is unreliable.
Here, the second threshold may be selected according to experience or accuracy required by a specific application scenario, for example, pseudorange location may be set to 15 meters, or two different thresholds in the horizontal direction and the vertical direction may be set for VPL and HPL, respectively, where the second threshold is a vector with a length of 2.
The difference between the invention and other robust positioning methods is that other methods only consider the standard error of the fault-free satellite
Figure BDA0003231578270000176
Failure to take into account the presence of a failed satellite or the observation containing multipath errors,
Figure BDA0003231578270000177
the value of (A) is influenced to increase, and the invention adds two parameters when positioning
Figure BDA0003231578270000178
For estimating the standard error of a faulty satellite or multipath effect, which results in
Figure BDA0003231578270000179
The estimated value is more accurate, and the positioning result is higher in precision. At the same time, due to the addition of
Figure BDA00032315782700001710
The fault or multi-path error of the satellite is absorbed, so that the positioning result has good robustness.
And step 10, judging whether each satellite has faults affecting positioning.
Computing
Figure BDA00032315782700001711
Obtaining test statistic:
Figure BDA00032315782700001712
representing the median of the absolute values of all residuals, the rejection field of this hypothesis test is { T }i||Ti|>T}。|Ti|>And T indicates that the ith satellite has a fault or deviation, otherwise, indicates that the ith satellite observation value has no fault or deviation.
T is a third threshold value, and T ═ Φ may be set-1(1-. alpha./2). T is an alpha quantile point of standard normal distribution, and alpha is a significance level determined after comprehensively considering false alarm and false alarm rate. Phi-1Is the inverse function of a standard normal cumulative distribution function. A may generally be chosen empirically or artificially based on the tolerable level of the fault, e.g. 5% may be chosen.
The method is different from other methods in that the method combines fault and multipath deviation identification and positioning, reduces the influence of a fault satellite on positioning by self-adaptive weight based on deviation probability and deviation value, can detect and eliminate the fault at the same time of positioning calculation, does not need to perform positioning calculation again, has few operation steps, avoids the problem that other methods judge and eliminate the fault satellite after positioning and perform positioning calculation again, reduces the operation amount, accelerates the operation speed, improves the positioning accuracy and simultaneously considers the robustness.
Meanwhile, the invention judges the fault satellite differently from other methods, which has or does not have two states, and is embodied by the fault occurrence probability and the deviation, so that whether the fault exists can be flexibly judged under different thresholds, the provided information is richer, and more effective information can be fed back to other receiver modules.
Based on the existing theory, experiments show that the receiver fault detection and anti-multipath steady positioning method based on Gaussian mixture distribution can quickly identify the fault satellite and determine the influence of the fault satellite on positioning in a self-adaptive manner, is a simple, efficient and quick fault satellite identification method for a user-side receiver, and provides reliability guarantee for positioning, time service and navigation by utilizing GPS, Beidou and dual modes.
The invention has been described in detail with reference to specific embodiments and illustrative examples, but the description is not intended to be construed in a limiting sense. Those skilled in the art will appreciate that various equivalent substitutions, modifications or improvements may be made to the technical solution of the present invention and its embodiments without departing from the spirit and scope of the present invention, which fall within the scope of the present invention. The scope of the invention is defined by the appended claims.
Those skilled in the art will appreciate that those matters not described in detail in the present specification are well known in the art.

Claims (11)

1. A non-Gaussian distribution based adaptive real-time multipath elimination and robust positioning method is characterized by comprising the following steps:
(1) establishing a positioning calculation model of the receiver in the positioning period, and linearly converting the positioning calculation model to obtain a standardized positioning calculation model Y which is H beta + epsilon, epsilon-N (0, sigma)2I) Wherein Y is (Y)1,...,yn)TFor a known normalized transformation of the observation matrix of each satellite, H ═ H1,...,hn)TStandardized transformation of a known observation positioning geometric matrix for N x p dimensions, N being the total number of satellites, beta being an unknown parameter vector for p x 1 dimensions, and epsilon being the standardized transformation of an observation noise matrix for N x 1 dimensions, obeying a normal distribution N (0, sigma)2I),σ2Is a variance coefficient, and I is a unit diagonal matrix;
(2) modeling each satellite observation using a mixture of Gaussian distributions, defining an n-dimensional unknown vector Z as a hidden variable, wherein each element Z is a variableiE {1,2,3}, i ═ 1, 2.. times, n, a gaussian mixture distribution model of each satellite observation is obtained
Figure FDA0003231578260000011
Wherein gamma iskIs ZiThe prior probability of (a) being,
Figure FDA0003231578260000012
order to
Figure FDA0003231578260000013
And a parameter vector to be solved consisting of all unknown parameters to be solved in the Gaussian mixture distribution model representing each satellite observed quantity, wherein mu is (mu)1,...,μK)T,μ1＝0,σ2＝(σ2 1,...,σ2 K)T,γ＝(γ1,...,γK)T,K＝3,γ3＝1-γ1-γ2；
(3) Iteratively solving the parameter vector theta to be solved based on the EM algorithm to obtain the optimal value of the parameter vector to be solved
Figure FDA0003231578260000014
Figure FDA0003231578260000015
(4) According to the optimal value of the parameter vector to be solved
Figure FDA0003231578260000016
Calculating the real-time protective level of positioning, comparing the real-time protective level with a given alarm threshold, and judging whether the positioning result is reliable or not;
(5) the optimal value of the parameter vector to be solved
Figure FDA0003231578260000021
In (1)
Figure FDA0003231578260000022
And (3) substituting the standard positioning calculation model obtained in the step (1) to obtain test statistic, and judging whether each satellite has a fault affecting positioning or not according to comparison between the test statistic and the rejection region of the hypothesis test.
2. The non-gaussian based adaptive real-time multipath cancellation and robust positioning method according to claim 1, wherein in the step (2), the gaussian mixture distribution model of each satellite observation is established by the following method:
defining an n-dimensional unknown vector Z as a hidden variable, wherein each element Zi∈{1,2,3},i＝1,2,...,n，Zi1 denotes the observation noise εiWithout deviations, Z, causing unacceptable errors in the positioning of the receiveri2 or Zi3 denotes the observation noise εiContain varying degrees of bias that cause unacceptable errors in receiver positioning;
assuming at different ZiUnder the value, each satellite observed quantity yiSubject to different Gaussian distributions, i.e.
Figure FDA0003231578260000023
Figure FDA0003231578260000024
Wherein y isiFor the i-th element, h, of the known satellite observation matrix YiTransposing the ith row of the geometric matrix H for known observationsi TRepresents hiTransposing; f. ofk(yi) For a given ZiTime yiConditional Gaussian distribution of (a) also being ZiUnknown time yiMixed gaussian density function fm(yi) The kth gaussian component of (1);
when Z isiAt unknown time, each satellite observed quantity yiFor gaussian mixture distribution, the gaussian mixture distribution model is obtained as follows:
Figure FDA0003231578260000025
wherein, γkIs ZiA priori of, P (Z)i＝k)＝γk。
3. The method according to claim 1, wherein 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
Figure FDA0003231578260000026
The method comprises the following steps:
(31) and calculating to obtain an expected log-likelihood function of each known satellite observation matrix Y:
Figure FDA0003231578260000031
wherein
Figure FDA0003231578260000032
t is the iteration number of the parameter vector theta to be solved;
(32) setting an iteration initial value theta of a parameter vector theta to be solved(0)Taking t as 0;
(33) according to theta(t)Substituting, calculating by Bayesian method
Figure FDA0003231578260000033
(34) Updating t to t +1, calculating in step (33)
Figure FDA0003231578260000034
Substituting into the desired log-likelihood function Q (θ | θ)(t-1)Y), solving for
Figure FDA0003231578260000035
Make Q (theta | theta)(t-1)Y) maximization to obtain
Figure FDA0003231578260000036
(35) Will be provided with
Figure FDA0003231578260000037
As gammakSubstituting the value of (a) into the desired log-likelihood function Q (θ | θ [ ])(t-1)Y), solving for
Figure FDA0003231578260000038
Make Q (theta | theta)(t-1)Y) in connection with
Figure FDA0003231578260000039
Maximization;
(36) repeating steps (33), (34) and (35), each time updating t to t +1 before starting step (34), until θ calculated in step (35)(t)And theta(t-1)Is smaller than a given first threshold value, theta is obtained when the final iteration number t is obtained(t)Optimum value of parameter vector to be solved
Figure FDA00032315782600000310
4. The non-Gaussian distribution-based adaptive real-time multipath elimination and robust positioning method according to claim 3, wherein in the step (32), the iterative initial value of the parameter vector θ to be solved
Figure FDA00032315782600000311
Wherein beta is(0)Estimated from equations of motion or state prediction in Kalman filtering, or determined using robust positioning methods, mu(0)，(σ2)(0)And gamma(0)Determined from historical empirical data.
5. The non-Gaussian distribution based adaptive real-time multipath cancellation and robust positioning method according to claim 3, wherein in the step (35), the solution is
Figure FDA0003231578260000041
Make Q (theta | theta)(t-1)Y) the step of maximization is as follows:
(351) fixation (sigma)2)(t)Then, Q (theta | theta) is obtained(t-1)Y) maximization
Figure FDA0003231578260000042
Figure FDA0003231578260000043
Figure FDA0003231578260000044
en＝(1,...,1)T
Figure FDA0003231578260000045
Figure FDA0003231578260000046
Figure FDA0003231578260000047
Therein
Figure FDA0003231578260000048
Use of
Figure FDA0003231578260000049
Or
Figure FDA00032315782600000410
The latest value of the two, step 351 and step 352, needs to be repeated for a plurality of times, and each time the value is updated
Figure FDA00032315782600000411
And
Figure FDA00032315782600000412
until the difference between the two update values is less than a given first threshold value;
(352) in that
Figure FDA00032315782600000413
Under the condition (2), Q (theta | theta) is obtained(t-1)Y) maximization2)(t)；
Figure FDA00032315782600000414
6. The non-Gaussian distribution based adaptive real-time multipath cancellation and robust positioning method of claim 3, wherein in the step (36), the first threshold is determined as beta and sigma calculated by multiple experiments25% of the average error.
7. The non-gaussian based adaptive real-time multi-path cancellation and robust positioning method according to claim 1, wherein in the step (4), the positioning real-time protection level comprises a horizontal protection level HPL and a vertical protection level VPL;
Figure FDA0003231578260000051
wherein, καSetting the variance expansion coefficient according to empirical value or actual false alarm and false alarm probability, C11，C22And C33Respectively are the optimal values of the parameter vectors to be solved
Figure FDA0003231578260000052
In
Figure FDA0003231578260000053
Variance of 1,2,3 elements of (a);
the alarm threshold is a vector with the length of 2 and comprises threshold values in the horizontal direction and the vertical direction.
8. The method of claim 7, wherein the adaptive real-time multipath cancellation and robust positioning method based on non-Gaussian distribution
Figure FDA0003231578260000054
Is a matrix of variances of each element
Figure FDA0003231578260000055
The first p diagonal elements in (1);
Figure FDA0003231578260000056
wherein A isE＝(HE TWEHE)-1HE TWE，
Figure FDA0003231578260000057
ΣYThe posterior covariance matrix of each known satellite observation matrix Y is a diagonal matrix, and the ith element in the diagonal element is:
Figure FDA0003231578260000058
wherein muk,
Figure FDA0003231578260000059
Are all the optimal values of the parameter vectors to be solved in the step (3)
Figure FDA00032315782600000510
Value of the corresponding element in (1), γikAccording to gamma(t)The value of the corresponding element is calculated by step (33).
9. The non-Gaussian distribution-based adaptive real-time multipath elimination and robust positioning method as claimed in claim 1, wherein in the step (5), the optimal value of the parameter vector to be solved is determined
Figure FDA00032315782600000511
In (1)
Figure FDA00032315782600000512
Substituting the standard positioning calculation model obtained in the step (1) to obtain test statistic as follows:
Figure FDA00032315782600000513
wherein
Figure FDA00032315782600000514
Represents the median of the absolute values of all residuals,
Figure FDA00032315782600000515
10. the non-Gaussian distribution-based adaptive real-time multipath elimination and robust positioning method as claimed in claim 1, wherein in the step (5), the rejection domain of the hypothesis testing is { T }i||Ti|＞T}，|TiIf the value is greater than T, the fault influencing the positioning exists in the ith satellite, otherwise, the fault influencing the positioning does not exist in the ith satellite;
t is a third threshold value, and T is phi-1(1-alpha/2), T is an alpha quantile point of standard normal distribution, alpha is a significance level determined after comprehensively considering false alarm rate and false alarm rate, and phi-1Is the inverse function of a standard normal cumulative distribution function.
11. The non-Gaussian distribution-based adaptive real-time multipath elimination and robust positioning method according to claim 1, wherein in the step (1), the positioning solution model of the receiver in the positioning period is
Figure FDA0003231578260000061
Wherein
Figure FDA0003231578260000062
Known observations are of dimension n x 1, n is the total number of satellites,
Figure FDA0003231578260000063
a geometric matrix is positioned for known observation of dimension n x p, beta is an unknown parameter vector of dimension p x 1,
Figure FDA0003231578260000064
an observation error vector is an n x 1-dimensional vector,
Figure FDA0003231578260000065
obeying a normal distribution N (0, σ)2∑),σ2Is a variance coefficient, sigma is a covariance prior matrix of the observed quantity error of each satellite,
Figure FDA0003231578260000066
the linearization conversion process is as follows:
Figure FDA0003231578260000067
after linear conversion
Figure FDA0003231578260000068
I is a unit diagonal matrix;
obtaining Y ═ H beta + epsilon, epsilon-N (0, sigma)2I)；
Wherein
Figure FDA0003231578260000069
CN202110988347.1A 2021-08-26 2021-08-26 Self-adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution Pending CN113848570A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110988347.1A CN113848570A (en) 2021-08-26 2021-08-26 Self-adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110988347.1A CN113848570A (en) 2021-08-26 2021-08-26 Self-adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution

Publications (1)

Publication Number Publication Date
CN113848570A true CN113848570A (en) 2021-12-28

Family

ID=78976377

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110988347.1A Pending CN113848570A (en) 2021-08-26 2021-08-26 Self-adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution

Country Status (1)

Country Link
CN (1) CN113848570A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116068593A (en) * 2023-01-28 2023-05-05 中国铁建电气化局集团有限公司 Satellite positioning weight calculation method, device, equipment and medium based on Bayes

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116068593A (en) * 2023-01-28 2023-05-05 中国铁建电气化局集团有限公司 Satellite positioning weight calculation method, device, equipment and medium based on Bayes

Similar Documents

Publication Publication Date Title
RU2559842C2 (en) Method and apparatus for detecting and eliminating multiple failures of gnss system satellites
US8106823B2 (en) Method of operating a satellite navigation receiver
WO2021022251A1 (en) System and method for gaussian process enhanced gnss corrections generation
US7956802B1 (en) Integrity-optimized receiver autonomous integrity monitoring (RAIM) for vertical integrity monitoring
CN110007317B (en) Star-selection optimized advanced receiver autonomous integrity monitoring method
US11480690B2 (en) System and method for satellite positioning
US20220107427A1 (en) System and method for gaussian process enhanced gnss corrections generation
RU2510529C2 (en) Method of determining navigation parameters for carrier and hybridisation device associated with kalman filter bank
CN113466903B (en) Partial ambiguity fixing algorithm considering observed value system error
CN115047496B (en) Synchronous multi-fault detection method for GNSS/INS integrated navigation satellite
US9298532B2 (en) Device and method for determining a physical quantity
CN115267855B (en) Abnormal value detection method and differential positioning method in GNSS-INS tight combination
CN115420284B (en) Fault detection and identification method for combined navigation system
CN114417552A (en) Ambiguity confirming method, storage medium and electronic equipment
CN113848570A (en) Self-adaptive real-time multipath elimination and robust positioning method based on non-Gaussian distribution
Green et al. Position-domain integrity analysis for generalized integer aperture bootstrapping
CN116859415A (en) Quick, stable and high-precision multi-fault satellite identification and positioning method
CN112083463B (en) Method and device for detecting whether ambiguity is fixed correctly and positioning terminal
CN115561782B (en) Satellite fault detection method in integrated navigation based on odd-even vector projection
CN112526549B (en) Method and system for identifying integrity fault of foundation enhancement system
CN115728793A (en) Precise single-point positioning gross error detection and processing method based on DIA theory
Wang et al. solution separation-based integrity monitoring for integer ambiguity resolution-enabled GNSS positioning
Ober Ways to improve RAIM/AAIM availability using position domain performance computations
CN111505670A (en) Multipath detection and suppression method and system using dual antennas
CN115390099B (en) Maximum value joint chi-square fault elimination method based on odd-even vector projection

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination