CN107728180B - GNSS precision positioning method based on multi-dimensional particle filter deviation estimation - Google Patents

GNSS precision positioning method based on multi-dimensional particle filter deviation estimation Download PDF

Info

Publication number
CN107728180B
CN107728180B CN201710790196.2A CN201710790196A CN107728180B CN 107728180 B CN107728180 B CN 107728180B CN 201710790196 A CN201710790196 A CN 201710790196A CN 107728180 B CN107728180 B CN 107728180B
Authority
CN
China
Prior art keywords
particle
particles
value
gnss
equation
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.)
Active
Application number
CN201710790196.2A
Other languages
Chinese (zh)
Other versions
CN107728180A (en
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.)
Southwest Jiaotong University
Original Assignee
Southwest Jiaotong University
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 Southwest Jiaotong University filed Critical Southwest Jiaotong University
Priority to CN201710790196.2A priority Critical patent/CN107728180B/en
Publication of CN107728180A publication Critical patent/CN107728180A/en
Application granted granted Critical
Publication of CN107728180B publication Critical patent/CN107728180B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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/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
    • G01S19/43—Determining position using carrier phase measurements, e.g. kinematic positioning; using long or short baseline interferometry
    • G01S19/44—Carrier phase ambiguity resolution; Floating ambiguity; LAMBDA [Least-squares AMBiguity Decorrelation Adjustment] method

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)

Abstract

The invention discloses a GNSS precision positioning method based on multi-dimensional particle filter deviation estimation, which comprises the following steps: step 1: data preprocessing, namely importing a satellite ephemeris, a current epoch pseudo-range observation value and a current epoch phase observation value; step 2: establishing a GNSS double-difference observation equation, and obtaining a single epoch method equation or a method equation accumulated with a previous epoch after linearization; and step 3: correcting corresponding deviation in the method equation by using particles with a plurality of particle values, and solving the method equation; updating the particle filtering weight, and calculating the numerical value of the deviation decimal part between the phase systems and the root mean square of the particles according to the weighted particles; and 4, step 4: repeating the step 1-3, and outputting a corresponding deviation estimated value and a corresponding variance after filtering convergence; and 5: and correcting the deviation value in the observed value or the normal equation according to the estimated value of the error parameter, and fixing the integer ambiguity to realize precise positioning. The method can realize simultaneous estimation of a plurality of GNSS deviation parameters by particle filtering, and can realize precise relative positioning by using less observation satellites.

Description

GNSS precision positioning method based on multi-dimensional particle filter deviation estimation
Technical Field
The invention relates to the technical field of satellite positioning systems and positioning measurement, in particular to a GNSS precision positioning method based on multi-dimensional particle filter deviation estimation.
Background
The application of a Global Navigation Satellite System (GNSS) in the field of precise positioning is more and more extensive; the GNSS observation values comprise tropospheric and ionospheric delay errors, FDMA signal inter-frequency bias of GLONASS, multi-system GNSS positioning inter-system bias and the like. In data processing for GNSS precision positioning, these errors in the phase observations need to be eliminated or estimated. More than one unknown parameter of error may be present in the GNSS phase observations at the same time, such as GLONASS inter-frequency bias and inter-system bias with other constellations may be present simultaneously in GNS multi-mode positioning where the GLONSS system participates. When the system combination frequency exceeds two groups, two or more than two intersystem deviation parameters exist; the traditional solution method uses a method based on least square method, and these parameters have certain correlation with the ambiguity parameters of GNSS, such as the phase inter-system bias and ambiguity parameters are completely correlated. At this time, equation rank deficiency or approximate rank deficiency occurs in the solution, and more observed values are needed to enhance the solution equation, or a priori values are adopted to eliminate/weaken the rank deficiency. Increasing the observation data will prolong the positioning time; the method using the prior value requires the prior value to have high reliability, needs to be determined in advance, and is difficult to handle the situation of error value change.
Disclosure of Invention
The invention provides a GNSS precision positioning method based on multi-dimensional particle filter deviation estimation, which can simultaneously estimate a plurality of error terms of the GNSS by using particle filter.
A GNSS precision positioning method based on multi-dimensional particle filter deviation estimation comprises the following steps:
step 1: preprocessing satellite navigation data, and importing a satellite ephemeris, a current epoch pseudo-range observation value and a current epoch phase observation value;
step 2: establishing an observed value equation containing a plurality of error parameters and linearizing to obtain a single epoch method equation or a previous epoch accumulation method equation;
and step 3: correcting corresponding deviation in the method equation by using particles with a plurality of particle values, and solving the method equation; fixing the ambiguity by an LAMBDA method, outputting a RATIO value, establishing a function related to the RATIO value, and updating the particle weight by using the function value; calculating the numerical value of the deviation decimal part between the phase systems and the root mean square of the particles according to the weighted particles;
and 4, step 4: repeating the steps 1-3, and outputting estimated values of unknown parameter vectors after filtering convergence, wherein the estimated values comprise estimated values of two or more error parameters;
and 5: and correcting the deviation value in the observed value or the normal equation according to the estimated value vector of the error parameter, and fixing the integer ambiguity to realize precise positioning.
Further, the specific process of step 3 is as follows:
(1) sampling in M dimensions produces an initial set of particles
Figure BDA0001398907340000021
For the kth moment, generating a particle set by a filtering result of the last moment; wherein x is a particle numerical value, w is a corresponding weight, N is the number of particles, i is 1, and 2 … N is a particle serial number; will be described below
Figure BDA0001398907340000022
Denoted as vector xi 0Corresponding to xi kFinger substitute
Figure BDA0001398907340000023
(2) If the error parameters contain the deviation parameters among the GNSS systems, calculating the root mean square of the corresponding particle dimensions; for each dimension, it is determined whether the root mean square is greater than a given threshold stdgroupIf the deviation is larger than the threshold value, cluster analysis is carried out on the particles, the clustered particles are integrated through the cluster analysis, and the step is skipped if the deviation of other types is larger than the threshold value;
(3) for each particle, correcting a corresponding deviation in a GNSS observation equation using the multi-dimensional value of the particle; resolving a law equation to obtain a floating point solution of an unknown quantity and a corresponding covariance matrix; fixing the ambiguity by an LAMBDA method, and outputting a RATIO value of the corresponding particle;
(4) establishing a likelihood function related to the RATIO value, updating the particle filtering weight by using the function value, and standardizing the particle weight as a new particle weight;
(5) calculating expected values of particles as estimates of unknown deviation vectors
Figure BDA0001398907340000024
Figure BDA0001398907340000025
Calculating the variance of the particles
Figure BDA0001398907340000026
Figure BDA0001398907340000027
(6) Judging whether the particle filter is convergent or not, and judging whether the root mean square is smaller than a set threshold std or notthdIf yes, outputting an estimated value of the phase deviation and the variance of the particles as an estimated value result;
(7) if the resampling condition is satisfied, resampling according to the updated weight;
(8) predicting particles at the next moment, and discretizing the resampled particles in real time:
Figure BDA0001398907340000028
wherein:
Figure BDA0001398907340000029
random noise added during discretization; and calculating the particle value of the next epoch, and transferring to the step 1.
Further, the process of integrating the clustered particles by clustering analysis in the step (2) is as follows:
s1: selecting two particles with the maximum and minimum particle values as initial particles, calculating the distance from other particles to the initial particles, and dividing the particles into two groups according to the distance;
s2: calculating the center of gravity g of two particle groups1，g2It is defined as follows:
Figure BDA0001398907340000031
wherein h is 1,2 is a group number, and Nh is the number of particles in each group; the gravity center distance of the particle group is as follows:
d＝|g1-g2|；
s3: judging that the difference between the distance d and the carrier wavelength lambda is smaller than eps, if | d-lambda | is smaller than eps, transferring one group of particles to the other group by adding or subtracting one wavelength on the dimension M to realize the combination of the two particle groups on the dimension M, wherein M is 1,2 …, M; wherein eps is generally a minimum;
s4: if the parameter vector is unknown
Figure BDA0001398907340000032
Is the phase intersystem deviation, the steps S1-S3 are repeated for the corresponding particle dimension value of the element.
Further, the particle filtering weight updating process in step (4) is as follows:
s11: establishing a functional relation between the likelihood function and RATIO:
Figure BDA0001398907340000033
in the formula: (RATIO) is a function of the RATIO value;
s12: according to the functional relationship established in step S11 and the RATIO value RATIO corresponding to the ith particleiCalculating likelihood function values of corresponding particles
Figure BDA0001398907340000034
Figure BDA0001398907340000035
S13: multiplying the likelihood function value with the weight of the corresponding particle to obtain the updated particle weight
Figure BDA0001398907340000036
Figure BDA0001398907340000037
S14: normalizing the weight of the particles, i.e. taking the ratio of the weight of each particle to the sum of all the weights of the particles as the new weight of the particles
Figure BDA0001398907340000038
Figure BDA0001398907340000039
Further, the resampling process in the step (7) is as follows:
s21: and accumulating the weight of the particles according to the serial numbers to obtain a cumulative distribution function value set of each particle:
Figure BDA0001398907340000041
s22: calculating the number of particles N requiredk+1：
Figure BDA0001398907340000042
In the formula:
Figure BDA0001398907340000043
n is the number of particles corresponding to the unit variance,
Figure BDA0001398907340000044
the number of particles is the minimum;
s23: generating a uniform or random cumulative distribution function value:
Figure BDA0001398907340000045
s24: comparing the cumulative distribution function value corresponding to the particle serial number with the uniform or random cumulative distribution function value in sequence; for m 1, i 1, e.g
Figure BDA0001398907340000046
Deleting the ith particle, i equals to i +1, otherwise copying the ith particle to a new particle set, and m equals to m + 1; until m is equal to Nk+1To obtain a new particle set of
Figure BDA0001398907340000047
S25: set the new set of particles to equal weight:
Figure BDA0001398907340000048
and obtaining a new particle set and a weight value.
Further, the establishment process of the single epoch method equation or the accumulation method equation with the previous epoch in the step 2 is as follows:
the pseudo-range non-difference observation equation of the GNSS system is as follows:
Figure BDA0001398907340000049
the GNSS system phase non-difference observation equation is as follows:
Figure BDA00013989073400000410
in the formula: i is the satellite serial number, a is the observation station serial number, P is the non-differential pseudo-range observed value of the GNSS satellite, phi is the non-differential phase observed value of the GNSS satellite, c is the light speed, delta taFor GNSS observation station receiver clock error, ρ is the distance between the observation station and the GNSS satellite, δ tiFor GNSS satellite clock error, di aFor pseudo-range hardware delay at the receiver, diPseudo-range hardware delay of GNSS satellite terminal, I is ionosphere delay error, T is troposphere delay error, epsilon is observation noise of pseudo-range observation value, mui aFor receiver-side phase hardware delay, muiFor GNSS satellite-side phase hardware delay, λiIs the carrier wavelength, N, of the ith satellitei aZeta is the observation noise of the phase observation value;
performing double-difference combination on the pseudo-range non-difference observation value and the phase non-difference observation value of the GNSS system, eliminating the clock difference of a satellite and the clock difference of a receiver, and correcting the delay error of an ionosphere and the delay error of a troposphere; obtaining a double-difference pseudo-range observation equation in the GNSS system:
Figure BDA0001398907340000051
the double-difference carrier phase observation equation in the GNSS system is as follows:
Figure BDA0001398907340000052
in the formula: s1Is any satellite system, swIs s is1Another satellite system other than the two, W is 2, 3, … W, where W is the number of combined satellite bands, b is the station number of the other station of the double-difference observation, j is the satellite number of the other GNSS satellite constituting the double-difference observation, d is the pseudorange intersystem bias, and μ is the phase intersystem bias, (k) is the phase intersystem biasj-ki)ΔγabThe frequency deviation is the frequency deviation when the GLONASS system FDMA observed value exists, k is the satellite number, and delta gamma is the frequency deviation rate;
the double-difference pseudo-range observation equation in the GNSS system and the double-difference carrier phase observation equation in the GNSS system can be converted into the following equations after linearization:
v＝Ax+Db+Cz+l
in the formula: x is a vector formed by coordinate components of measuring stations except ambiguity and inter-frequency deviation, b is a vector of single-difference ambiguity unknowns among receivers, z is an unknown vector containing a plurality of deviations, A, D and C are coefficient matrixes corresponding to the unknowns respectively, l is a constant term vector, P is a weight matrix, and v is an observed value residual error vector;
according to the linearized equation, a single epoch method equation or an addition equation with a previous epoch can be used:
Figure BDA0001398907340000053
the invention has the beneficial effects that:
(1) the method can realize simultaneous estimation of a plurality of deviation parameters of the GNSS by using particle filtering, including GLONSS phase inter-frequency deviation, multi-mode GNSS phase inter-system deviation, troposphere delay and the like, and realize GNSS precision positioning;
(2) in the GNSS multimode combined positioning, the method can utilize the ambiguity whole-cycle characteristic among GNSS systems, realize the solution of the deviation among phase systems when the satellite number is less, and achieve the aim of precise positioning.
Drawings
FIG. 1 is a schematic flow chart of the present invention.
FIG. 2 is a RAITO distribution diagram of GNSS two phase system bias estimation in the embodiment of the present invention.
Fig. 3 illustrates the convergence process of the full constellation in the embodiment of the present invention, in which fig. a to d sequentially show the particle distribution maps and the error estimation positions of the first epoch to the fourth epoch, the background is the slope distribution map of each epoch, the granular point is the particle position, and the pentagram is the deviation estimation position.
Fig. 4 shows PRN numbers for six satellites used in an embodiment of the invention.
FIG. 5 is a sequence number of an estimation result when the systematic offsets of two phases of the GNSS six satellites are estimated simultaneously according to the embodiment of the present invention.
FIG. 6 is a comparison of baseline solution results after correction of the offset between two phase systems in an embodiment of the present invention.
Detailed Description
The invention is further described with reference to the following figures and specific embodiments.
As shown in fig. 1, a GNSS precision positioning method based on multidimensional particle filter bias estimation includes the following steps:
step 1: preprocessing satellite navigation data, and importing a satellite ephemeris, a current epoch pseudo-range observation value and a current epoch phase observation value;
the preprocessing generally includes data format conversion, gross error detection and elimination, cycle slip detection and repair, and the like.
Step 2: and establishing an observed value equation containing a plurality of error parameters and linearizing to obtain a single epoch method equation or an accumulation method equation with the previous epoch.
The pseudo-range non-difference observation equation of the GNSS system is as follows:
Figure BDA0001398907340000061
the GNSS system phase non-difference observation equation is as follows:
Figure BDA0001398907340000062
wherein: i is the satellite serial number, a is the observation station serial number, P is the non-differential pseudo-range observed value of the GNSS satellite, phi is the non-differential phase observed value of the GNSS satellite, c is the light speed, delta taFor GNSS observation station receiver clock error, ρ is the distance between the observation station and the GNSS satellite, δ tiFor GNSS satellite clock error, di aFor pseudo-range hardware delay at the receiver, diPseudo-range hardware delay of GNSS satellite terminal, I is ionosphere delay error, T is troposphere delay error, epsilon is observation noise of pseudo-range observation value, mui aFor receiver-side phase hardware delay, muiFor GNSS satellite-side phase hardware delay, λiIs the carrier wavelength, N, of the ith satellitei aZeta is the observed noise of the phase observation for the whole-cycle ambiguity.
Performing double-difference combination on the pseudo-range of the GNSS system and the non-differential observation value of the phase of the GNSS system to obtain a double-difference pseudo-range observation equation in the GNSS system and a double-carrier-difference phase observation equation in the GNSS system:
Figure BDA0001398907340000071
Figure BDA0001398907340000072
in the formula: s1Is any satellite system, swFor another satellite system, W is 1,2, 3, … W, where W is the number of combined satellite bands, b is the station number of another station of the double-difference observation, j is the satellite number of another GNSS satellite constituting the double-difference observation, d is the pseudorange intersystem bias, and μ is the phase intersystem bias, which requires accurate estimation, (k) andj-ki)Δγabthe GLONASS systemInter-frequency bias at FDMA observations, k is the satellite number, Δ γ is the inter-frequency bias ratio, which is zero without GLONASS system FDMA observations, when the combined frequency is large, multiple inter-phase systematic biases, possibly GLONASS inter-frequency biases, and when the baseline is long, tropospheric and ionospheric delay biases, etc. need to be estimated simultaneously.
After linearization, a double-difference pseudo-range observation equation in the GNSS system and a double-carrier-difference phase observation equation in the GNSS system can be converted into:
v＝Ax+Db+Cz+l
in the formula: x is a vector formed by coordinate components of measuring stations except ambiguity and inter-frequency deviation, b is a vector of single-difference ambiguity unknowns between receivers, z is a vector of deviation unknowns, A, D and C are coefficient matrixes corresponding to the unknowns respectively, l is a constant term vector, P is a weight matrix, and v is an observation value residual vector;
according to the linearized equation, a single epoch method equation or an addition equation accumulated with a previous epoch can be obtained:
Figure BDA0001398907340000073
and step 3: correcting corresponding deviation in the method equation by using particles with a plurality of particle values, and solving the method equation; fixing the ambiguity by an LAMBDA method, outputting a RATIO value, establishing a function related to the RATIO value, and updating the particle weight by using the function value; calculating the numerical value of the deviation decimal part between the phase systems and the root mean square of the particles according to the weighted particles; if the root mean square of the particles is smaller than a certain threshold, the filtering is stopped, and the specific process is as follows:
(1) sampling in M dimensions produces an initial set of particles
Figure BDA0001398907340000074
For the kth moment, generating a particle set by a filtering result of the last moment; wherein x is a particle numerical value, w is a corresponding weight, N is the number of particles, i is 1, and 2 … N is a particle serial number; will be described below
Figure BDA0001398907340000081
Denoted as vector xi 0Corresponding to xi kFinger substitute
Figure BDA0001398907340000082
(2) If the error parameters contain the deviation parameters among the GNSS systems, calculating the root mean square of the corresponding particle dimensions; for each dimension, the following operations are performed: judging whether the corresponding root mean square is larger than a given threshold stdgroupIf the particle size is larger than the threshold value, performing cluster analysis on the particles, and dividing the particles into two groups; if the distance between the particle centers of gravity is approximately one carrier wavelength, then the two particle populations are merged, skipping this step for other deviations.
The process of integrating clustered particles by clustering analysis is as follows:
s1: selecting two particles with the maximum and minimum particle values as initial particles, calculating the distance from other particles to the initial particles, and dividing the particles into two groups according to the distance;
s2: calculating the center of gravity g of two particle groups1，g2It is defined as follows:
Figure BDA0001398907340000083
wherein h is 1,2 is a group number, and Nh is the number of particles in each group; the gravity center distance of the particle group is as follows:
d＝|g1-g2|；
s3: judging that the difference between the distance d and the carrier wavelength lambda is smaller than a smaller value eps, if | d-lambda | is smaller than eps, transferring one group of particles to the other group by adding or subtracting one wavelength on the dimension M to realize the combination of the two particle groups on the dimension M, wherein M is 1,2 …, M; where eps is a small value relative to the carrier wavelength;
s4: if the parameter loss is unknown
Figure BDA0001398907340000084
Is a phase intersystem deviation, thenThe steps S1-S3 are repeated for the particle dimension values corresponding to the element.
(3) For each particle, correcting a corresponding deviation in a GNSS observation equation using the multi-dimensional value of the particle; resolving a law equation to obtain a floating point solution of an unknown quantity and a corresponding covariance matrix; fixing the ambiguity by an LAMBDA method, and outputting a RATIO value of the corresponding particle;
(4) establishing a likelihood function related to the RATIO value, updating the particle filtering weight by using the function value, and standardizing the particle weight as a new particle weight;
the particle filter weight update process is as follows:
s11: establishing a functional relation between the likelihood function and RATIO:
Figure BDA0001398907340000085
in the formula: (RATIO) is a function of the RATIO value;
s12: according to the functional relationship established in step S11 and the RATIO value RATIO corresponding to the ith particleiCalculating likelihood function values of corresponding particles
Figure BDA0001398907340000091
Figure BDA0001398907340000092
S13: multiplying the likelihood function value with the weight of the corresponding particle to obtain the updated particle weight
Figure BDA0001398907340000093
Figure BDA0001398907340000094
S14: normalizing the weight of the particles, i.e. taking the ratio of the weight of each particle to the sum of all the weights of the particles as the new weight of the particles
Figure BDA0001398907340000095
Figure BDA0001398907340000096
(5) Calculating expected values of particles as estimates of unknown deviation vectors
Figure BDA0001398907340000097
Figure BDA0001398907340000098
Calculating the variance of the particles
Figure BDA0001398907340000099
Figure BDA00013989073400000910
(6) Judging whether the particle filter is convergent or not, and judging whether the root mean square is smaller than a set threshold std or notthdIf yes, outputting an estimated value of the phase deviation and the variance of the particles as an estimated value result;
(7) if the resampling condition is satisfied, resampling according to the updated weight;
the resampling process is as follows:
s21: and accumulating the weight of the particles according to the serial numbers to obtain a cumulative distribution function value set of each particle:
Figure BDA00013989073400000911
s22: calculating the number of particles N requiredk+1：
Figure BDA00013989073400000912
In the formula:
Figure BDA00013989073400000913
n is the number of particles corresponding to the unit variance,
Figure BDA00013989073400000914
the number of particles is the minimum;
s23: generating a uniform or random cumulative distribution function value:
Figure BDA0001398907340000101
s24: comparing the cumulative distribution function value corresponding to the particle serial number with the uniform or random cumulative distribution function value in sequence; for m 1, i 1, e.g
Figure BDA0001398907340000102
Deleting the ith particle, i equals to i +1, otherwise copying the ith particle to a new particle set, and m equals to m + 1; until m is equal to Nk+1To obtain a new particle set of
Figure BDA0001398907340000103
S25: set the new set of particles to equal weight:
Figure BDA0001398907340000104
and obtaining a new particle set and a weight value, and outputting the generated new particle set and the weight value to the next step.
(8) Predicting particles at the next moment, and discretizing the resampled particles in real time:
Figure BDA0001398907340000105
wherein:
Figure BDA0001398907340000106
random noise added during discretization; the particle value of the next epoch is calculated,and (5) transferring to the step 1.
And 4, step 4: repeating the steps 1-3, and outputting estimated values of unknown parameter vectors after filtering convergence, wherein the estimated values comprise estimated values of two or more error parameters;
and 5: and correcting the deviation value in the observed value or the normal equation according to the estimated value of the error parameter, and fixing the integer ambiguity to realize precise positioning.
The method of the invention is used for simultaneously estimating the deviation between two GNSS phase systems and processing a short baseline, wherein the observed value of the baseline comprises a GPS system L5, a Galieo system E5a and a QZSS system L5; the deviation parameters between two phase systems of the GPS L5-Galileo E5a and the GPS L5-QZSS L5 need to be estimated simultaneously; the RAITO distribution diagram of GNSS two phase system bias simultaneous estimation is shown in fig. 2, and it can be seen that the RATIO value can map the values of two parameters at the same time; the convergence process of estimating the offset between two phase systems using the method with all the observation satellites is shown in fig. 3, the particles can be successfully converged, and two parameter values are estimated, wherein the problem of the half cycle of the offset between the phase systems in one dimension (shown in fig. 3 c) is successfully solved (shown in fig. 3 d); simulating the situation that one six observation satellites come from three constellations, wherein the PRN numbers of the six used satellites are shown in FIG. 4; used for simulating the received satellite and giving the used data information; the estimated two intersystem bias results are shown in fig. 5, and it can be seen that a small amount of bias results from three GNSS constellations are used, and the method can still achieve correct estimation of the two intersystem bias; the baseline calculation result obtained finally is shown in fig. 6, the positioning result is obviously improved after the method of the invention is used, and the ambiguity fixed proportion is also improved from 19.3% to 99.7%.
The method can realize simultaneous estimation of a plurality of deviation parameters of the GNSS by using particle filtering, including GLONASS phase inter-frequency deviation, multi-mode GNSS phase inter-system deviation, troposphere delay and the like, and realize GNSS precision positioning; in the GNSS multimode combined positioning, the method can utilize the ambiguity whole-cycle characteristic among GNSS systems, realize the solution of the phase intersystem deviation when the satellite number is less, and achieve the purpose of precise positioning.
In this context: lambda is an abbreviation for Least-square AMBiguity decay Adjustment, referring to the Least squares AMBiguity Decorrelation Adjustment; the reliability inspection index of RATIO in the GNSS field when the ambiguity of the whole cycle is fixed; PRN is an abbreviation for Pseudo Random Noise, referring to Pseudo Random Noise.

Claims (5)

1. A GNSS precision positioning method based on multi-dimensional particle filter deviation estimation is characterized by comprising the following steps:
step 1: preprocessing satellite navigation data, and importing a satellite ephemeris, a current epoch pseudo-range observation value and a current epoch phase observation value;
step 2: establishing an observed value equation containing a plurality of error parameters and linearizing to obtain a single epoch method equation or a previous epoch accumulation method equation;
and step 3: correcting corresponding deviation in the method equation by using particles with a plurality of particle values, and solving the method equation; fixing the ambiguity by an LAMBDA method, outputting a RATIO value, establishing a function related to the RATIO value, and updating the particle weight by using the function value; according to the weighted particles, calculating the numerical value and the root mean square of the particles of the deviation decimal part between the phase systems, specifically:
step (1): sampling in M dimensions produces an initial set of particles
Figure FDA0002700641660000011
For the kth moment, generating a particle set by a filtering result of the last moment; wherein x is a particle numerical value, w is a corresponding weight, N is the number of particles, i is 1, and 2 … N is a particle serial number; will be described below
Figure FDA0002700641660000012
Denoted as vector xi 0Corresponding to xi kFinger substitute
Figure FDA0002700641660000013
Step (2): if the error parameters contain the deviation parameters among the GNSS systems, calculating the root mean square of the corresponding particle dimensions; for eachDimension, judging whether the root mean square is larger than a given threshold stdgroupIf the deviation is larger than the threshold value, cluster analysis is carried out on the particles, the clustered particles are integrated through the cluster analysis, and the step is skipped if the deviation of other types is larger than the threshold value;
and (3): for each particle, correcting a corresponding deviation in a GNSS observation equation using the multi-dimensional value of the particle; resolving a law equation to obtain a floating point solution of an unknown quantity and a corresponding covariance matrix; fixing the ambiguity by an LAMBDA method, and outputting a RATIO value of the corresponding particle;
and (4): establishing a likelihood function related to the RATIO value, updating the particle filtering weight by using the function value, and standardizing the particle weight as a new particle weight;
and (5): calculating expected values of particles as estimates of unknown deviation vectors
Figure FDA0002700641660000014
Figure FDA0002700641660000015
Calculating the variance of the particles
Figure FDA0002700641660000016
Figure FDA0002700641660000017
And (6): judging whether the particle filter is convergent or not, and judging whether the root mean square is smaller than a set threshold std or notthdIf yes, outputting an estimated value of the phase deviation and the variance of the particles as an estimated value result;
and (7): if the resampling condition is satisfied, resampling according to the updated weight;
and (8): predicting particles at the next moment, and discretizing the resampled particles in real time:
Figure FDA0002700641660000021
wherein:
Figure FDA0002700641660000022
random noise added during discretization; calculating the particle value of the next epoch, and transferring to the step (1);
and 4, step 4: repeating the steps 1-3, and outputting estimated values of unknown parameter vectors after filtering convergence, wherein the estimated values comprise estimated values of two or more error parameters;
and 5: and correcting the deviation value in the observed value or the normal equation according to the estimated value of the error parameter, and fixing the integer ambiguity to realize precise positioning.
2. The GNSS fine positioning method based on multi-dimensional particle filter bias estimation as claimed in claim 1, wherein the process of integrating clustered particles by cluster analysis in step (2) is as follows:
s1: selecting two particles with the maximum and minimum particle values as initial particles, calculating the distance from other particles to the initial particles, and dividing the particles into two groups according to the distance;
s2: calculating the center of gravity g of two particle groups1，g2It is defined as follows:
Figure FDA0002700641660000023
wherein h is 1,2 is a group number, and Nh is the number of particles in each group; the gravity center distance of the particle group is as follows:
d＝|g1-g2|；
s3: judging a difference value eps between the distance d and the carrier wave wavelength lambda, if | d-lambda | is less than eps, transferring one group of particles to the other group by adding or subtracting one wavelength on a dimension M to realize the combination of the two groups of particles on the dimension M, wherein M is 1,2 …, M;
s4: if the parameter loss is unknown
Figure FDA0002700641660000024
Is the phase intersystem deviation, the steps S1-S3 are repeated for the corresponding particle dimension value of the element.
3. The GNSS precision positioning method based on multidimensional particle filter bias estimation as recited in claim 1, wherein the particle filter weight updating process in step (4) is as follows:
s11: establishing a functional relation between the likelihood function and RATIO:
Figure FDA0002700641660000025
in the formula: (RATIO) is a function of the RATIO value;
s12: according to the functional relationship established in step S11 and the RATIO value RATIO corresponding to the ith particleiCalculating likelihood function values of corresponding particles
Figure FDA0002700641660000031
Figure FDA0002700641660000032
S13: multiplying the likelihood function value with the weight of the corresponding particle to obtain the updated particle weight
Figure FDA0002700641660000033
Figure FDA0002700641660000034
S14: normalizing the weight of the particles, i.e. taking the ratio of the weight of each particle to the sum of all the weights of the particles as the new weight of the particles
Figure FDA0002700641660000035
Figure FDA0002700641660000036
4. The GNSS fine positioning method based on multi-dimensional particle filter bias estimation according to claim 1, wherein the resampling process in step (7) is as follows:
s21: and accumulating the weight of the particles according to the serial numbers to obtain a cumulative distribution function value set of each particle:
Figure FDA0002700641660000037
s22: calculating the number of particles N requiredk+1：
Figure FDA0002700641660000038
In the formula:
Figure FDA0002700641660000039
n is the number of particles corresponding to the unit variance,
Figure FDA00027006416600000310
the number of particles is the minimum;
s23: generating a uniform or random cumulative distribution function value:
Figure FDA00027006416600000311
s24: comparing the cumulative distribution function value corresponding to the particle serial number with the uniform or random cumulative distribution function value in sequence; for m 1, i 1, e.g
Figure FDA00027006416600000312
Deleting the ith particle, i equals to i +1, otherwise copying the ith particle to a new particle set, and m equals to m + 1; until m is equal to Nk+1To obtain a new particle set of
Figure FDA00027006416600000313
S25: set the new set of particles to equal weight:
Figure FDA0002700641660000041
and obtaining a new particle set and a weight value.
5. The GNSS fine positioning method based on multi-dimensional particle filter bias estimation according to claim 1, wherein the establishment of the single epoch method equation or the cumulative method equation with the previous epoch in step 2 is as follows:
the pseudo-range non-difference observation equation of the GNSS system is as follows:
Figure FDA0002700641660000042
the GNSS system phase non-difference observation equation is as follows:
Figure FDA0002700641660000043
in the formula: i is the satellite serial number, a is the observation station serial number, P is the non-differential pseudo-range observed value of the GNSS satellite, phi is the non-differential phase observed value of the GNSS satellite, c is the light speed, delta taFor GNSS observation station receiver clock error, ρ is the distance between the observation station and the GNSS satellite, δ tiFor GNSS satellite clock error, di aFor pseudo-range hardware delay at the receiver, diPseudo-range hardware delay for GNSS satellite, I is electricityAn error in the delayer delay, T the tropospheric delay error, ε the observation noise of the pseudorange observations, μi aFor receiver-side phase hardware delay, muiFor GNSS satellite-side phase hardware delay, λiIs the carrier wavelength, N, of the ith satellitei aZeta is the observation noise of the phase observation value;
performing double-difference combination on the pseudo-range non-difference observed value of the GNSS system and the phase non-difference observed value of the GNSS system to obtain double-difference pseudo-range and phase observation equations in the GNSS system, wherein the double-difference pseudo-range and phase observation equations are respectively as follows:
Figure FDA0002700641660000044
Figure FDA0002700641660000045
in the formula: s1Is any satellite system, swW is 1,2, 3, … W, where W is the number of satellite systems, b is the station number of the other station of the double-difference observation, j is the satellite number of the other GNSS satellite constituting the double-difference observation, d is the pseudorange intersystem bias, μ is the phase intersystem bias, (k) is the phase intersystem biasj-ki)ΔγabThe frequency deviation is the frequency deviation when the GLONASS system FDMA observed value exists, k is the satellite number, and delta gamma is the frequency deviation rate;
the double-difference pseudo-range observation equation in the GNSS system and the double-difference carrier phase observation equation in the GNSS system can be converted into the following equations after linearization:
v＝Ax+Db+Cz+l
in the formula: x is a vector formed by coordinate components of measuring stations except ambiguity and inter-frequency deviation, b is a vector of single difference ambiguity unknowns among receivers, z is an unknown vector containing a plurality of deviations to be estimated, A, D and C are coefficient matrixes corresponding to the unknowns respectively, l is a constant term vector, P is a weight matrix, and v is an observed value residual error vector;
according to the linearized equation, a single epoch method equation or an addition equation with the previous epoch can be obtained:
Figure FDA0002700641660000051
CN201710790196.2A 2017-09-05 2017-09-05 GNSS precision positioning method based on multi-dimensional particle filter deviation estimation Active CN107728180B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710790196.2A CN107728180B (en) 2017-09-05 2017-09-05 GNSS precision positioning method based on multi-dimensional particle filter deviation estimation

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710790196.2A CN107728180B (en) 2017-09-05 2017-09-05 GNSS precision positioning method based on multi-dimensional particle filter deviation estimation

Publications (2)

Publication Number Publication Date
CN107728180A CN107728180A (en) 2018-02-23
CN107728180B true CN107728180B (en) 2021-01-29

Family

ID=61204906

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710790196.2A Active CN107728180B (en) 2017-09-05 2017-09-05 GNSS precision positioning method based on multi-dimensional particle filter deviation estimation

Country Status (1)

Country Link
CN (1) CN107728180B (en)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108919321B (en) * 2018-05-18 2019-05-10 长安大学 A kind of GNSS positioning Detection of Gross Errors method based on trial and error method
CN108955851B (en) * 2018-07-12 2020-07-17 北京交通大学 Method for determining GNSS error by using INS and DTM
CN110568464B (en) * 2019-06-19 2023-10-10 航天信息股份有限公司 BDS/GNSS multimode chip-based precise positioning method and BDS/GNSS multimode chip-based precise positioning device
CN110988935B (en) * 2019-11-25 2021-11-12 重庆市地理信息和遥感应用中心(重庆市测绘产品质量检验测试中心) Multi-system combination precision positioning method based on receiver-side deviation clustering optimization
CN111596321B (en) * 2020-05-29 2022-04-01 武汉大学 Multi-GNSS multi-path error star day filtering method and system using non-difference correction
CN113064195B (en) * 2021-03-16 2023-09-01 西南交通大学 High-precision low-computation carrier attitude measurement method utilizing geometric characteristics of multiple antennas
CN112799096B (en) * 2021-04-08 2021-07-13 西南交通大学 Map construction method based on low-cost vehicle-mounted two-dimensional laser radar
CN113504557B (en) * 2021-06-22 2023-05-23 北京建筑大学 Real-time application-oriented GPS inter-frequency clock difference new forecasting method
CN113536983B (en) * 2021-06-29 2023-12-12 辽宁工业大学 Petroleum pipeline stealing positioning method based on P-RLS adaptive filtering time delay estimation
CN113589349A (en) * 2021-09-30 2021-11-02 中国地质大学（武汉） GNSS real-time high-precision single-difference attitude measurement method based on particle filtering
CN114488227B (en) * 2022-01-26 2023-10-20 西南交通大学 Multipath error correction method based on spatial correlation

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105723185A (en) * 2013-11-25 2016-06-29 高通股份有限公司 Motion state based mobile device positioning

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105607106B (en) * 2015-12-18 2018-08-21 重庆邮电大学 A kind of low-cost and high-precision BD/MEMS fusions attitude measurement method
KR20170084896A (en) * 2016-01-13 2017-07-21 한국전자통신연구원 Apparatus for estimating indoor position information of user base on paticle filters method and using the same
CN106249256B (en) * 2016-07-08 2018-08-14 辽宁工程技术大学 Real-time GLONASS phase deviation estimation methods based on particle swarm optimization algorithm

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105723185A (en) * 2013-11-25 2016-06-29 高通股份有限公司 Motion state based mobile device positioning

Also Published As

Publication number Publication date
CN107728180A (en) 2018-02-23

Similar Documents

Publication Publication Date Title
CN107728180B (en) GNSS precision positioning method based on multi-dimensional particle filter deviation estimation
CN107678050B (en) GLONASS phase inter-frequency deviation real-time tracking and precise estimation method based on particle filtering
CN107728171B (en) Particle filter based real-time tracking and precise estimation method for deviation between GNSS phase systems
CN111045034B (en) GNSS multi-system real-time precise time transfer method and system based on broadcast ephemeris
Bisnath Efficient, automated cycle-slip correction of dual-frequency kinematic GPS data
CN1864078B (en) Method for using three GPS frequencies to resolve carrier-phase integer ambiguities
EP1654559B1 (en) Method for generating clock corrections for a wide-area or global differential gps system
US6127968A (en) On-the-fly RTK positioning system with single frequency receiver
CN108254773B (en) Real-time clock error resolving method of multiple GNSS
CN108873029B (en) Method for realizing clock error modeling of navigation receiver
CN109143298B (en) Beidou and GPS observation value cycle slip detection and restoration method, equipment and storage equipment
CN109061694B (en) Low-orbit navigation enhanced positioning method and system based on GNSS clock correction fixation
CN110764127B (en) Relative orbit determination method for formation satellite easy for satellite-borne on-orbit real-time processing
CA2557984A1 (en) Method for back-up dual-frequency navigation during brief periods when measurement data is unavailable on one of two frequencies
CN107966722B (en) GNSS clock error resolving method
CN110346823B (en) Three-frequency ambiguity resolving method for Beidou precise single-point positioning
CN112146557A (en) GNSS-based real-time bridge deformation monitoring system and method
CN114935770B (en) Method and device for accelerating precision single-point positioning convergence speed by multiple calendars
CN113204042A (en) Multi-constellation combined train positioning method based on precise single-point positioning
CN111323796B (en) GNSS receiver high-sampling clock error resolving method
CN111123322B (en) Method, system, medium and device for preprocessing observed value real-time data of satellite-borne GNSS receiver
CN114994727A (en) Equipment for realizing high-precision time calibration and satellite positioning
CN110568464A (en) BDS/GNSS (broadband navigation satellite system/global navigation satellite system) multi-mode chip-based precision positioning method and device
Liu et al. Generating GPS decoupled clock products for precise point positioning with ambiguity resolution
CN110988935B (en) Multi-system combination precision positioning method based on receiver-side deviation clustering optimization

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
GR01 Patent grant
GR01 Patent grant