Patent Grant
11709281
This application is a 371 of international application of PCT application serial no. PCT/CN2020/072411, filed on Jan. 16, 2020, which claims the priority benefit of China application no. 201910428335.6, filed on May 22, 2019. The entirety of each of the above mentioned patent applications is hereby incorporated by reference herein and made a part of this specification.
The present invention relates to a satellite positioning method for Global Navigation Satellite Systems (GNSSs), and in particular to a high-precision point positioning method and device based on original GNSS observation values of a smartphone.
Smartphones are indispensable tools in people's lives nowadays, and GNSS (Global Navigation Satellite System) modules in smartphones have greatly improved the lives of modern people. In the development of GNSS navigation and positioning technology, navigation or positioning precision has always been a key problem that restricts GNSS navigation and positioning technology from being further applied in human production and life and playing a significant role, and the same is true for GNSS navigation and positioning of smartphones. Because GNSS navigation and positioning modules in smartphones have always been encapsulated in operating systems of smartphones, researchers can only acquire final positioning results to carry out application-level development. Therefore, relative to conventional geodesic receivers, the analysis of original observation values of smartphones and the research of positioning algorithms are very scarce.
Since 2016, Google has provided an interface for accessing original GNSS observation data in the Android N operating system, and scholars have begun to evaluate the quality of GNSS data received by smartphones and analyze the positioning performance. At present, domestic and foreign scholars' researches on using original GNSS observation values to achieve high-precision positioning of smartphones are mainly divided into three directions: (1) analyzing the quality of original GNSS observation value data of a mobile end; (2) comprehensively using various types of GNSS data to increase positioning precision by a filtering method; (3) using a pseudorange or carrier difference method to increase the positioning precision of a mobile end. All these researches apply existing navigation and positioning enhancement means into smartphones, and most of research objects are Android tablets. For smartphones used by ordinary people, the current positioning performance is poor. More importantly, there is no precedent to improve conventional positioning solutions according to the nature of original observation values of smartphones at present.
Objective of the invention: Based on the above background, the PPP (Precise Point Positioning) algorithm is improved to achieve the high-precision real-time point positioning of smartphones, aiming at the inconsistency between pseudorange observation values and carrier observation values of smartphones found in the test.
Technical solution: In order to achieve the aforementioned objective of the present invention, the present invention adopts the following technical solution.
A high-precision point positioning method based on a smartphone comprises the following steps:
(1) acquiring original observation values of a navigation and positioning module of the smartphone, such as GNSS pseudoranges and carrier phases;
(2) after preprocessing the data to decrease part of error influences, generating an uncombined observation value model from the original observation values according to an improved precise point positioning method based on an estimation of double clock biases;
(3) determining each satellite observation value weight according to a satellite elevation angle; and
(4) carrying out calculation by employing an improved static Kalman filtering to give a high-precision point positioning result.
In a preferred embodiment, in step (1), the API based on location service provided by the operating system of the smartphone is used to directly acquire part of original GNSS data, including time data, GNSS system type and carrier phase observation values, and a pseudorange is then calculated by a signal propagation time difference according to the time data among the original GNSS data.
In a preferred embodiment, in step (2), the decrease of error influence includes using a precise ephemeris and a precise clock bias file to eliminate an orbit error and a satellite clock bias for received observation values, using an ionospheric grid file to reduce ionospheric delay and regarding a multi-path effect as observation noise, and the generated uncombined observation value model is:
Pig=ρig+cd{tilde over (t)}P+dtropg+ϵPg
Φig=ρig+cd{tilde over (t)}Φ+dtropg+Ñig+ϵΦg
Pje=ρje+cd{tilde over (t)}P+cd{tilde over (t)}syse+dtrope+ϵPe
Φje=ρje+cd{tilde over (t)}Φ+cd{tilde over (t)}syse+dtrope+Ñje+ϵΦe
Pkc=ρkc+cd{tilde over (t)}P+cd{tilde over (t)}sysc+dtropc+ϵPc
Φkc=ρkc+cd{tilde over (t)}Φ+cd{tilde over (t)}sysc+dtropc+Ñkc+ϵΦc
in the formula, the superscripts g, e and c respectively represent the GPS system, the Galileo system and the BDS system, the subscripts i, j and k represent ith, jth and kth satellites, P and Φ on the left are respectively pseudorange and carrier observation values, ρ is the distance from a smartphone terminal to a satellite, c is light speed, d{tilde over (t)}P and d{tilde over (t)}Φ are respectively mobile end clock biases of pseudorange and carrier observation values, dtrop is a tropospheric delay, N is carrier integer ambiguity, and ϵP and ϵΦ are respectively residuals of pseudorange and carrier phase observation values; and d{tilde over (t)}syse and d{tilde over (t)}sysc respectively represent a time bias between the Galileo system and the GPS system and a time bias between the BDS system and the GPS system.
In a preferred embodiment, in step (3), the weight determination solution for determining the observation value weight of each satellite according to a satellite elevation angle is as follows:
when the satellite elevation angle is less than 10°, the observation value weight of the satellite is 0;
when the satellite elevation angle is greater than 10°, the observation value weight of the satellite is sin E, and E is a satellite elevation angle.
In a preferred embodiment, in step (4), improved static Kalman filtering is as follows:
for a single epoch, according to the uncombined observation value model described in step (2), the number of observation equations for a satellite is two, and if n1 GPS satellites, n2 Galileo satellites and n3 BDS satellites are observed in a certain epoch, n1+n2+n3≥8 is met; at this point, a parameter vector to be estimated in Kalman filtering is:
XT=[xyzcd{tilde over (t)}Pcd{tilde over (t)}Φcd{tilde over (t)}sysecd{tilde over (t)}syscdtropÑ1g . . . Ñn
in the formula, x, y and z are positional parameters of the smartphone; then the total number of observation equations is 2×(n1+n2+n3), the number of parameters to be estimated is n1+n2+n3+8, and the number of redundant observations is n1+n2+n3−8;
the observation equation coefficient matrix H of Kalman filtering is:
In the matrix H, the row number is 2×(n1+n2+n3), the column number is n1+n2+n3−8, the subscripts 1 to n represent satellite numbers, n=2×(n1+n2+n3), and the superscripts represent GNSS system types, that is, g, e and c respectively represent the GPS system, the Galileo system and the BDS system; the odd rows correspond to pseudorange observation values, and the even rows correspond to carrier phase observation values; the first three columns α, β and γ are satellite-mobile end direction cosines, the fourth and fifth columns are clock bias coefficients of pseudorange observation values and carrier observation values, the odd rows of the fourth column are 1, and the even rows of the fifth column are 1; the sixth and seventh columns are inter-system bias coefficients; when the observation values come from the Galileo satellites, the sixth column is 1; when the observation values come from the BDS satellites, the seventh column is 1; MF is tropospheric wet delay projection coefficients; the columns after the ninth column are carrier phase ambiguity coefficients, the odd rows are 0, the even rows of the i+8th column are 1, and i is satellite number; and in the process of filtering, the weight of each satellite observation value is determined in step (3).
Based on the same inventive concept, the present invention provides a high-precision point positioning device based on a smartphone, which comprises a memory, a processor and a computer program stored in the memory and capable of running in the processor. When loaded into the processor, the computer program implements the high-precision point positioning method based on a smartphone.
Beneficial effects: Based on the conventional positioning algorithm, the high-precision point positioning method based on a smartphone proposed by the present invention optimizes the uncombined PPP observation value model according to the unique nature of original GNSS observation values of the smartphone terminal, and can achieve sub-meter level positioning precision of ordinary smartphones, and the convergence speed is fast. Based on the high-precision positioning technology of the smartphone terminal, the present invention can provide a better user experience in urban positioning, on-board navigation, tourism transportation and so on for users. In the coming era of intelligent interconnection, location service with higher precision means an unpredictable application prospect. Low-cost location service is included in the Internet of Things, automated driving, smart city and other development directions of modern people's lives, so the present invention is of great significance.
The technical solution of the present invention will be further described in detail below with reference to the drawings.
In the first step, original observation values, such as GNSS pseudoranges and carrier phases, are acquired through the API based on location service provided by an Android system. Part of original GNSS data, including time data, a GNSS system type, carrier phases (the first row in Table 1), a pseudorange rate and a carrier-to-noise ratio, are first directly acquired by using the GnssMeasurement class and the GnssClock class in an android.location package in Android system development. The original data contents included in the GnssMeasurement class and the GnssClock class are shown as Table 1 and Table 2.
Pseudoranges are then calculated according to the time data among the original GNSS data, the essence of which is to calculate the pseudoranges by a signal propagation time difference:
P=(tRx−tTx)·c (1)
Wherein P represents a pseudorange, tRx represents the moment when a signal is received by the receiver, tTx represents the moment when a signal transmitted by a satellite is received, and c is light speed. tTx can be directly acquired by the getReceivedSvTimeNanos( ) method, and value is the number of seconds in the GPS week.
tRx=tRx
tRx
wherein tRx
In the second step, after the data are preprocessed to decrease part of error influences, an uncombined observation value model is generated from the original observation values according to an improved PPP algorithm based on an estimation of double clock biases.
In a preliminary experiment, the smartphone and the geodesic receiver are placed at the place for synchronous observation, original observation values are collected, and the original observation value file of the smartphone is converted into a Rinex file.
The different original satellite observation values of the same device are analyzed, finding that the differences between pseudorange change rates and carrier observation value change rates of different satellites at the same moment are equal, so it is considered that the pseudorange observation values and carrier observation values of the smartphone terminal have different clock biases during positioning.
According to the aforementioned nature, the positioning solution is improved to estimate two clock bias parameters during parameter calculation. Taking the GPS system as an example, the single-frequency equation is:
Pig=ρig+cd{tilde over (t)}P−cdTg+dorbg+dtropg+diong+dmult/Pg+ϵPg (4)
Φig=ρig+cd{tilde over (t)}Φ−cdTg+dorbg+dtropg+diong+Ñig+dmult/Φg+ϵPg (5)
In the formula, the superscript g represents the GPS system, the subscript i represents the ith satellites, Pig and Φig are pseudorange and carrier observation values, ρig is the distance from the smartphone terminal to a satellite, c is light speed, d{tilde over (t)}P and d{tilde over (t)}Φ are respectively clock biases of pseudorange and carrier observation values, dTg is a satellite clock bias, dtropg is a tropospheric delay, dorbg is a satellite orbit error, diong is an ionospheric delay, dmult/Pg and dmult/Φg are respectively multipath effect delays of pseudoranges and carriers, Ñig is carrier integer ambiguity, and ϵPg and ϵΦg are respectively residuals of pseudorange and carrier phase observation values.
A predicted precise ephemeris, a clock bias file and a predicted ionospheric grid file provided by the IGS (International GNSS Service) are downloaded in real time to eliminate satellite orbit errors and clock biases and decrease ionospheric delay errors, and the influence of errors such as the relativistic effect and the rotation of the earth is corrected and eliminated by a model. After the multipath effect is regarded as observation noise, the observation value equation is simplified to give an uncombined observation value model:
Pig=ρig+cd{tilde over (t)}P+dtropg+ϵPg (6)
Φig=ρig+cd{tilde over (t)}Φ+dtropg+Ñig+ϵΦg (7)
Pje=ρje+cd{tilde over (t)}P+cd{tilde over (t)}syse+dtrope+ϵPe (8)
Φje=ρje+cd{tilde over (t)}Φ+cd{tilde over (t)}syse+dtrope+Ñje+ϵΦe (9)
Pkc=ρkc+cd{tilde over (t)}P+cd{tilde over (t)}sysc+dtropc+ϵPc (10)
Φkc=ρkc+cd{tilde over (t)}Φ+cd{tilde over (t)}sysc+dtropc+Ñkc+ϵΦc (11)
In the formula, the pseudorange and carrier observation value of each GNSS system are on the left, the superscripts e and c respectively represent the Galileo system and the BDS system, the subscripts j and k represent the jth satellite and the kth satellite, and d{tilde over (t)}syse and d{tilde over (t)}sysc respectively represent a time bias between the Galileo system and the GPS system and a time bias between the BDS system and the GPS system.
In the third step, each satellite observation value weight is determined. In the method, the weight determination solution for determining the observation value weight of each satellite according to a satellite elevation angle is as follows:
In the formula, W is a weight, and E is an elevation angle between the smartphone terminal and a satellite.
In the fourth step, Kalman filtering is employed to carry out parameter estimation to give a high-precision positioning result.
The uncombined model shown in the second step is generated for the observation values. At this point, the number of observation equations for a satellite in a single epoch is two, and if n1 GPS satellites, n2 Galileo satellites and n3 BDS satellites are observed in a certain epoch, at this point, a parameter vector to be estimated in Kalman filtering is:
XT=[xyzcd{tilde over (t)}Pcd{tilde over (t)}Φcd{tilde over (t)}syscd{tilde over (t)}sysdtropÑ1g . . . Ñn
In the formula, x, y and z are positional parameters of the smartphone.
Then the total number of observation equations is 2×(n1+n2+n3), the number of parameters to be estimated is n1+n2+n3+8, and the number of redundant observations is n1+n2+n2−8.
The specific steps of Kalman filtering employed by the present invention are as follows.
(1) an initial system state value X0 and a variance P0 thereof are inputted. X is the parameter vector to be estimated in formula (13), and the subscript 0 represents initial moment, that is, initial values are set for the parameters to be estimated; at this point, the first three parameters are approximate coordinates (acquired by standard point positioning using pseudorange) of the smartphone, and the rest of the parameters are respectively pseudorange clock bias, carrier clock bias, inter-system bias, tropospheric wet delay and ambiguity of n1+n2+n3 satellites, which are all set as 0.
(2) further predicted values Xk,k−1 and a variance-covariance matrix Pk,k−1 (the subscript k is the number of epochs) are calculated:
Xk,k−1=Φk,k−1Xk−1 (14)
Pk,k−1=Φk,k−1Pk−1ΦTk,k−1+Γk,k−1Qk−1ΓTk,k−1 (15)
In the formula, Φk,k−1 is a state transition matrix (because the number of satellites may change, the nth-order Xk−1 of the previous epoch needs to be changed into the mth-order Xk of this epoch). Γk,k−1 is a system noise driving matrix, Qk−1 is a positive definite matrix of system errors (errors of the model), input is required for each epoch, and k,k−1Qk−1 k,k−1 can be directly regarded as a whole, defined as a mth-order symmetric matrix. It is the same as P0 during initial positioning. During positioning in subsequent epochs, the true values of three positional parameters and n ambiguity parameters should not be changed, so process noise is 0 (variance corresponding to system noise). Random models for smartphone clock bias, inter-system bias and tropospheric delay parameters are simulated by using random walk process models.
(3) a gain matrix Kk is calculated:
Kk=Pk,k−1HkT(HkPk,k−1HkT+Rk)−1 (16)
In the formula, Hk is a coefficient matrix for observation equations. For conventional static Kalman filtering, Hk in the method needs to be modified accordingly:
The row number of Hk is 2×(n1+n2+n3), the column number is n1+n2+n2+8, and in the formula, n=2×(n1+n2+n3); the subscripts represent satellite numbers, and the superscripts represent GNSS system types; the odd rows correspond to pseudorange observation values, and the even rows correspond to carrier phase observation values; the first three columns α, β and γ are satellite-mobile end direction cosines, and the calculation formula is shown as formula (18); the fourth and fifth columns are clock bias coefficients of pseudorange observation values and carrier observation values, the odd rows of the fourth column are 1, and the even rows of the fifth column are 1; the sixth and seventh columns are inter-system bias coefficients; when the observation values come from the Galileo satellites, the sixth column is 1; when the observation values come from the BDS satellites, the seventh column is 1; and in the eighth column, MF is tropospheric wet delay projection coefficients (the tropospheric dry delay is corrected by a Hopfield model, and the projection function is a Niell model). The ninth column to the n1+n2+n2+8th column are carrier phase ambiguity coefficients, the 2×ith row of the i+8th column is 1, and i is satellite number.
In the formula, xs, ys and zs are satellite coordinates (acquired through an ephemeris file), x0, y0 and z0 are approximate coordinates (acquired by standard point positioning using pseudorange) of the smartphone, and ρ0 is an approximate distance (calculated through satellite coordinates and approximate smartphone coordinates) from the smartphone terminal to a satellite.
In formula (16), Rk is an observation noise variance matrix:
In the formula, σP and σΦ are respectively zenithal standard deviations of pseudorange observation values and carrier observation values, and W is a weight of each satellite, i.e., each satellite observation value weight calculated by formula (12). In the process of Kalman filtering, the setting of variances of observation values is very important, and improper setting can easily cause filtering divergence, severely affecting the positioning result. In the present invention, the zenithal standard deviations of observation values are determined according to the quality of the original GNSS observation value data of the smartphone detected in the previous experiment. The zenithal standard deviations of pseudorange observation values and carrier phase observation values of the Xiaomi 8 smartphone used in the example are respectively set as 3.0 m and 0.2 m.
(4) the estimated filtering parameter value vector Xk of the current epoch and a corresponding variance-covariance matrix Pk are calculated:
Xk=Xk,k−1+Kk(Lk−HkXk,k−1) (20)
Pk=(I−KkHk)Pk,k−1(I−KkHk)T+KkRkKkT (21)
In the formula, Lk is an observation value vector, i.e., the pseudorange and carrier observation values of each GNSS system in formulas (6) to (11), and I is a unit matrix.
After calculation is completed, the first three items in the Xk vector are accurate coordinates of the smartphone calculated in the current epoch. If the current epoch is not the last epoch, then return to step (2) to carry out circular filtering calculation until a multi-epoch filtering solution is obtained in the end.
The time of the experiment was 2018.10.19, a Xiaomi 8 smartphone was used, and the test site was the GE01 control point in the campus of the Southeast University. Observations were carried out five times, about 6 minutes each time, and the sampling interval was 1 s. The accurate coordinates of the observation station had been obtained in advance by means of network RTK.
The method of the present invention was employed to carry out a positioning test on the Xiaomi 8 smartphone. The result is shown as
Errors of positioning in the five time periods were counted, as shown in Table 3. The average error in plane was 0.81 m, while the average error in elevation was 1.65 m. This result is better than the highest positioning precision that can be achieved currently by using ordinary smartphones.
The convergence time is defined as the time from the moment when positioning is started to the moment when the errors of positioning in both N and E directions are less than 1 m and the errors in subsequent epochs no longer exceed 1 m. The convergence times in the five time periods are shown as table 4. Convergence can be achieved within 30 s in each time period, indicating that the method can be applied to real-time positioning and provide a low-delay high-precision smartphone positioning result.
Based on the same inventive concept, an example of the present invention discloses a high-precision point positioning device based on a smartphone, which comprises a memory, a processor and a computer program stored in the memory and capable of running in the processor. When loaded into the processor, the computer program implements the high-precision point positioning method based on a smartphone. A test result indicates that by using the positioning method proposed by the present invention, a positioning precision higher than 0.9 m in plane and 1.7 m in elevation can be achieved for the smartphone.
| Number | Date | Country | Kind |
|---|---|---|---|
| 201910428335.6 | May 2019 | CN | national |
| Filing Document | Filing Date | Country | Kind |
|---|---|---|---|
| PCT/CN2020/072411 | 1/16/2020 | WO |
| Publishing Document | Publishing Date | Country | Kind |
|---|---|---|---|
| WO2020/233158 | 11/26/2020 | WO | A |
| Number | Name | Date | Kind |
|---|---|---|---|
| 20140002299 | Leandro et al. | Jan 2014 | A1 |
| Number | Date | Country |
|---|---|---|
| 104597465 | May 2015 | CN |
| 104714244 | Jun 2015 | CN |
| 107356947 | Nov 2017 | CN |
| 107356947 | Nov 2017 | CN |
| 107807373 | Mar 2018 | CN |
| 108363079 | Aug 2018 | CN |
| 109343090 | Feb 2019 | CN |
| 110275192 | Sep 2019 | CN |
| WO-2011034624 | Mar 2011 | WO |
| Entry |
|---|
| “International Search Report (Form PCT/ISA/210) of PCT/CN2020/072411,” dated Apr. 20, 2020, pp. 1-5. |
| Number | Date | Country | |
|---|---|---|---|
| 20220155465 A1 | May 2022 | US |
