CN109541663A - The correcting method of appearance Multipath Errors is surveyed in a kind of GNSS positioning - Google Patents

Abstract

The invention discloses a method for correcting multi-path errors of GNSS positioning and attitude measurement. The trend surface analysis method is used to fit the multi-path spatial distribution characteristics in the sky grid, which improves the modeling ability of the multi-path modeling method based on spatial repeatability for high-frequency multi-path, and improves the positioning and attitude measurement accuracy; The requirements for the scale of the sky grid are improved, and the calculation efficiency is improved. This method is suitable for both static and dynamic scenes with unchanged multi-path environments.

Description

A Correction Method for Multipath Error of GNSS Positioning and Attitude Measurement

technical field

The invention relates to the field of satellite positioning and navigation, in particular to a method for correcting multi-path errors in GNSS positioning and attitude measurement.

Background technique

The modeling method based on multipath spatiotemporal repeatability can realize the real-time processing of multipath effects in GNSS positioning. Among them, the multi-path modeling methods based on spatial repeatability, such as spherical harmonic function model and multi-path hemisphere map (MHM) method, are not only suitable for static positioning mode, but also for dynamic carriers with invariant multi-path environments such as ships and aircraft. positioning, and therefore has significant advantages. But these methods have obvious limitations. The accuracy of the spherical harmonic model depends on the fitting order, and it is difficult to take into account the model accuracy and computational efficiency. Although the computational efficiency of the MHM method is greatly improved compared with the spherical harmonic model, it uses the average residual error in the grid as the correction value, ignoring the difference of multi-path values in the grid, resulting in poor high-frequency multi-path modeling. Fine mesh segmentation can ensure accuracy.

SUMMARY OF THE INVENTION

The purpose of the present invention is to provide a method for correcting multi-path errors in GNSS positioning and attitude measurement. This method uses the trend surface method to fit the spatial distribution characteristics of multi-paths in the sky map grid, which improves the modeling ability of high-frequency multi-paths, and at the same time reduces the requirements for grid scale to a certain extent and improves the calculation. efficiency.

The concrete technical scheme that realizes the object of the present invention is:

A method for correcting multi-path errors in GNSS positioning and attitude measurement, the method includes the following specific steps:

Step 1: Short Baseline Solution and Multipath Computation

(1) Establish a common clock single-difference model

Based on the common clock multi-antenna receiver to solve the short baseline, the inter-station single difference can eliminate the atmospheric delay, satellite clock difference and receiver clock difference, and the common clock single difference model is obtained:

where Δ represents the single-difference operator; λ is the carrier wavelength; φ ⁱ , ^{ri , N i and φ i mth are} _the observed value of the carrier phase, the geometric distance between the satellite and the antenna, the ambiguity parameter and the multipath effect, respectively; φ _UPD is UPD (uncalibrated phasedelay, phase delay) parameter;

(2) Solve the short baseline

Kalman filtering is used to solve the single-difference model; at k epochs, the matrix form of the observation equation is:

Z _k =H _k X _k +ε _k (2)

X=[dB _x dB _y dB _z VB _x VB _y VB _z -ΔN ¹ … -ΔN ^m Δφ _UPD ] ^T (3)

Among them, Z _k is the residual vector of the observation value, that is, the single-difference observation minus the initial geometric distance, the initial ambiguity and the initial UPD value; X _k is the state vector; [dB _x dB _y dB _z ] ^T represents the baseline vector increment; [VB _x VB _y VB _z ] ^T is the baseline velocity vector increment; ΔN ⁱ is the ambiguity parameter increment of the i-th satellite; H _k is the design matrix, l _x , _ly , and l _z are the satellite-receiver connection Cosine of the direction of the line; (X ⁱ , Y ⁱ , Z ⁱ ) are satellite coordinates; (X ⁰ , Y ⁰ , Z ⁰ ) are the approximate coordinates of the receiver; ρ ⁰ is the distance between the satellite and the approximate coordinates of the receiver ε _k is the observation noise; the superscript T represents the transpose of the matrix;

From epoch k to epoch k+1, the prediction equation of the state vector and its covariance matrix are:

X _k+1,k =A _k X _k (6)

The update formula of the state vector and its covariance matrix are:

X _k+1 =X _k+1,k +K _k+1 (Z _k+1 -A _k+1 X _k+1,k ) (8)

C _k+1 =C _k+1,k -K _k+1 A _k+1 C _k+1,k (9)

Among them, A _k is the state transition matrix, Q _k is the error covariance matrix of the state transition equation, P _k+1 is the variance matrix of the observation noise, K _k+1 is the gain matrix, R _k is the variance matrix of the observation noise, C _{k+1, k} and C _k+1 are the prediction variance matrix and the update variance matrix, respectively;

The floating-point solution Ni of the ambiguity parameter can be obtained by Kalman filtering, and the method of fixing the floating-point ambiguity is as follows: select a reference satellite ⁱ , and combine the floating-point ambiguity ^Ni and UPD parameters of the satellite into a new UPD parameter ; At this time, the new UPD parameter is actually the sum of UPD and the floating-point ambiguity of the reference satellite, while the floating-point ambiguity of other non-reference satellites becomes the difference between the original floating-point ambiguity and the floating-point ambiguity of the reference satellite; as follows The formula shows:

in, is the new UPD parameter obtained by merging; Ni is the floating-point ambiguity of the reference satellite ⁱ , k=1,2...i-1,i+i,...,n is the new floating-point ambiguity of satellite k after the UPD parameters are combined with the reference satellite floating-point ambiguity; the ambiguity is fixed in three steps: Firstly, ignoring the integer characteristics of ambiguity, the parameter estimation method is used to solve the baseline vector and the ambiguity floating-point solution; then the ambiguity floating-point solution and its covariance matrix are used to fix the ambiguity satellite by satellite according to its integer characteristics; Resolve the original observation equation by replacing the original observation equation to obtain the baseline vector solution after the ambiguity is fixed;

(3) Calculate the multipath based on the obtained solution

The posterior residual of the Kalmam filter solution is the multipath of the satellite:

Step 2: Establishment of a multi-path hemisphere model based on trend surface analysis

(1) Assign multipath data to the sky grid

Calculate the altitude and azimuth for the satellites that have obtained multipath data; first, calculate the observation vector [Δe Δn Δu] ^T from the user to the satellite in the station center coordinate system:

Among them, [Δx Δy Δz] ^T is the observation vector from the user to the satellite in the geocentric geo-fixed rectangular coordinate system; λ and φ are the longitude and latitude, respectively; Next, calculate the observation vector from the user to the satellite in the carrier coordinate system [Δx′ Δy′ Δz ′] ^T :

Among them, ψ, θ and γ are the carrier heading angle, pitch angle and roll angle respectively; the satellite altitude angle ele and azimuth angle azi are:

Divide the sky into different grids according to the altitude and azimuth, and assign the satellite multipath data corresponding to the altitude and azimuth to the grid;

(2) Build a multi-path trend surface model for each grid of the sky

First, model the primary and secondary trend surfaces for the multipath of the sky grid:

A trend surface:

Secondary trend surface:

where x _i , y _i represent the azimuth and elevation angles corresponding to the i- _th multipath value in the sky grid, a ₀ , a ₁ ,...an represent the trend surface fitting coefficients, is the fitting value of the trend surface model;

Secondly, carry out the R ² test and the significant F test for the appropriate fit; in which the coefficient of fit R ² is represented by the proportion of the regression square sum SSR to the total deviation square sum _{S T} ; the value of _R ² ranges from 0 to 1. The larger the value, the higher the fitting degree of the model; its calculation formula is:

where _zi is is the residual sum of squares

If R ² is greater than a predetermined threshold, the model is judged to pass the test;

If both models pass the test and then use the F test to select the best, the F test is used to test whether the trend surface equation is significant, and its calculation formula is:

where p represents the degrees of freedom of the regression sum of squares, given by Calculated, (np-1) is the degree of freedom of the remaining sum of squares, r is the number of trend surface analysis, n is the total number of residual values in the grid point; F obeys the F distribution with degrees of freedom (p, np-1), In the case of a given confidence level, by checking the F distribution table, the critical value F is obtained, so as to test which trend surface equation fitted is more significant under the confidence level α;

If none of them pass, the mean value of residuals in the grid is used to replace the trend surface model directly;

Step 3: Multipath correction in real time

Collect GNSS observation data in real time, calculate the altitude and azimuth angles of satellites in the carrier coordinate system, select the multi-path trend surface model corresponding to the sky map grid for multi-path estimation, and correct the satellite observation data.

The present invention combines the trend surface analysis method with the traditional MHM method. One of its advantages is that it follows the multi-path modeling method based on spatial repeatability and is suitable for dynamic positioning and navigation applications in which static and multi-path environments remain unchanged. The trend surface method fits the spatial distribution characteristics of multi-paths in the sky map grid, which improves the modeling ability of high-frequency multi-paths, and reduces the requirements for the grid scale and improves the computational efficiency.

Description of drawings

Fig. 1 is the flow chart of the present invention;

Fig. 2 is the power spectrum analysis comparison diagram of the method of the present invention and the model value of the multipath hemisphere;

Figure 3 is a comparison diagram of the baseline solutions of the two models under different grid resolutions.

Detailed ways

The present invention will be described in detail below with reference to the accompanying drawings and embodiments.

Example

The invention is suitable for multi-path error correction of static and multi-path environment-invariant dynamic carriers in GNSS positioning and navigation applications. Compared with the existing multi-path modeling methods based on spatial repeatability, the present invention improves the modeling capability of high-frequency multi-paths, reduces the requirement for grid scale, and improves computing efficiency.

In order to verify the present invention, a static short-baseline multi-path correction experiment is implemented, and the specific implementation steps include:

Step 1: Short Baseline Solution and Multipath Computation

(1) Establish a common clock single-difference model

Based on the common clock multi-antenna receiver to solve the short baseline, the inter-station single difference can eliminate the atmospheric delay, satellite clock difference and receiver clock difference, and the common clock single difference model is obtained:

where Δ represents the single-difference operator; λ is the carrier wavelength; φ ⁱ , ^{ri , N i and φ i mth are} _the observed value of the carrier phase, the geometric distance between the satellite and the antenna, the ambiguity parameter and the multipath effect, respectively; φ _UPD is a UPD (uncalibrated phasedelay, phase delay) parameter.

(2) Solve the short baseline

The single-difference model is solved using Kalman filtering. At epoch k, the matrix form of the observation equation is:

Z _k =H _k X _k +ε _k (2)

X=[dB _x dB _y dB _z VB _x VB _y VB _z -ΔN ¹ … -ΔN ^m Δφ _UPD ] ^T (3)

Among them, Z _k is the residual vector of the observation value, that is, the single difference observation minus the initial geometric distance, the initial ambiguity and the initial UPD value; X _k is the state vector; [dB _x dB _y dB _z ] ^T represents the baseline vector increment; [VB _x VB _y VB _z ] ^T is the baseline velocity vector increment; ΔN ⁱ is the ambiguity parameter increment of the i-th satellite; H _k is the design matrix, l _x , _ly , and l _z are the satellite-receiver connection Cosine of the direction of the line; (X ⁱ , Y ⁱ , Z ⁱ ) are satellite coordinates; (X ⁰ , Y ⁰ , Z ⁰ ) are the approximate coordinates of the receiver; ρ ⁰ is the distance between the satellite and the approximate coordinates of the receiver _εk is the observation noise. The superscript T represents the transpose of the matrix.

From epoch k to epoch k+1, the prediction equation of the state vector and its covariance matrix are:

X _k+1,k =A _k X _k (6)

The update formula of the state vector and its covariance matrix are:

X _k+1 =X _k+1,k +K _k+1 (Z _k+1 -A _k+1 X _k+1,k ) (8)

C _k+1 =C _k+1,k -K _k+1 A _k+1 C _k+1,k (9)

Among them, A _k is the state transition matrix, Q _k is the error covariance matrix of the state transition equation, P _k+1 is the variance matrix of the observed noise, K _k+1 is the gain matrix, R _k is the variance matrix of the observed noise, C _{k+1, k} and C _k+1 are the prediction variance matrix and the update variance matrix, respectively.

The floating-point solution Ni of the ambiguity parameter can be obtained by Kalman filtering, and the method of fixing the floating-point ambiguity is as follows: select a reference satellite ⁱ , and combine the floating-point ambiguity ^Ni and UPD parameters of the satellite into a new UPD parameter ; At this time, the new UPD parameter is actually the sum of UPD and the floating-point ambiguity of the reference satellite, while the floating-point ambiguity of other non-reference satellites becomes the difference between the original floating-point ambiguity and the floating-point ambiguity of the reference satellite; as follows The formula shows:

in, is the new UPD parameter obtained by merging; Ni is the floating-point ambiguity of the reference satellite ⁱ , k=1,2...i-1,i+i,n is the new floating-point ambiguity of satellite k after the UPD parameter is combined with the reference satellite floating-point ambiguity; the ambiguity is fixed in three steps: first ignore Integer characteristics of ambiguity, the parameter estimation method is used to solve the baseline vector and ambiguity floating-point solution; then the ambiguity floating-point solution and its covariance matrix are used to fix the ambiguity satellite by satellite according to its integer characteristics; The original observation equation is re-solved to obtain the baseline vector solution after the ambiguity is fixed.

(3) Calculate the multipath based on the obtained solution

The posterior residual of the Kalmam filter solution is the multipath of the satellite:

Step 2: Establishment of a multi-path hemisphere model based on trend surface analysis

(1) Assign multipath data to the sky grid

Elevation and azimuth are calculated for satellites for which multipath data has been acquired. First, calculate the observation vector [Δe Δn Δu] ^T from the user to the satellite in the station center coordinate system:

Among them, [Δx Δy Δz] ^T is the observation vector from the user to the satellite in the geocentric geo-fixed rectangular coordinate system; λ and φ are the longitude and latitude, respectively; Next, calculate the observation vector from the user to the satellite in the carrier coordinate system [Δx′ Δy′ Δz ′] ^T :

Among them, ψ, θ and γ are the carrier heading angle, pitch angle and roll angle respectively; the satellite altitude angle ele and azimuth angle azi are:

Divide the sky into different grids according to the altitude and azimuth angles, and the division unit can be in degrees (such as 1°*1°, azimuth angle×altitude angle); you can also consider the equal area division method, that is, take a certain area as the unit Perform grid division (by keeping the interval of the height angle fixed, the size of the azimuth angle interval is different to keep the area of each grid point consistent: for example, when the grid point area is M°*N°, the height angle is between 0° and N°. M grid points are divided between, and M*COS (N/360*2pi) grid points are divided between N° and 2N° for the height angle. Then, the satellite multipath data corresponding to the altitude and azimuth are allocated to the grid.

(2) Build a multi-path trend surface model for each grid of the sky

First, model the primary and secondary trend surfaces for the multipath of the sky grid:

A trend surface:

Secondary trend surface:

where x _i , y _i represent the azimuth and elevation angles corresponding to the i- _th multipath value in the sky grid, a ₀ , a ₁ ,...an represent the trend surface fitting coefficients, is the fitting value of the trend surface model;

Secondly, carry out the R ² test and the significant F test for the appropriate fit; in which the coefficient of fit R ² is represented by the proportion of the regression square sum SSR to the total deviation square sum _{S T} ; the value of _R ² ranges from 0 to 1, the larger the value, the higher the fit of the model; its calculation formula is

where _zi is is the residual sum of squares

If R ² is larger than a certain threshold, the model is judged to pass the test.

If both models pass the test and then use the F test to select the best, the F test is used to test whether the trend surface equation is significant, and its calculation formula is:

In the formula, p represents the degree of freedom of the regression sum of squares, which is given by Calculated, (np-1) is the degree of freedom of the remaining sum of squares, r is the number of trend surface analysis, n is the total number of residual values in the grid point; F obeys the F distribution with degrees of freedom (p, np-1), In the case of a given confidence level, the critical value F can be obtained by checking the F distribution table, so as to test which trend surface equation fitted is more significant at the confidence level α.

If none of them pass, the mean value of residuals in the grid is used to replace the trend surface model directly;

Step 3: Multipath correction in real time

Collect GNSS observation data in real time, calculate the altitude and azimuth of satellites in the carrier coordinate system, select the multi-path trend surface model corresponding to the sky map grid for multi-path estimation, and correct the satellite observation data.

By comparing the method of the present invention with the existing multi-path modeling methods based on spatial repeatability (such as multi-path hemisphere map, MHM), the method of the present invention is better than MHM in the correction effect of high-frequency multi-path. Figure 2 is the result of extracting the model value of one of the satellites (taking G06 as an example) for power spectrum analysis after using the multipath models of two methods. As shown in the figure, the power of the method of the present invention is higher than that of the MHM model in the entire frequency band , indicating that the method of the present invention contains more abundant multi-path changes than the MHM model, which can not only perform better correction of low-frequency multi-path, but also improve the correction performance of high-frequency multi-path.

The sky is divided into grids of 1°*1°, 2°*2°, 5°*5°, and 9°*9° respectively. The comparison results of the baseline solutions of the two methods are shown in Figure 3. The results show that with the sky As the grid increases, the baseline solution of the method of the present invention is more stable, indicating that the method of the present invention has lower requirements on the fineness of the sky grid, so the algorithm is more efficient.

Claims (1)

1. a correction method for GNSS positioning attitude measurement multi-path error, is characterized in that, the method comprises the following concrete steps: Step 1: short baseline solution and multi-path calculation (1) Establish a common clock single-difference model Based on the common clock multi-antenna receiver to solve the short baseline, the inter-station single difference can eliminate the atmospheric delay, satellite clock difference and receiver clock difference, and the common clock single difference model is obtained: where Δ represents the single-difference operator; λ is the carrier wavelength; φ ⁱ , ^{ri , N i and φ i mth are} _the observed value of the carrier phase, the geometric distance between the satellite and the antenna, the ambiguity parameter and the multipath effect, respectively; φ _UPD is the phase delay UPD parameter; (2) Solve the short baseline Kalman filtering is used to solve the single-difference model; at k epochs, the matrix form of the observation equation is: Z _k =H _k X _k +ε _k (2) X=[dB _x dB _y dB _z VB _x VB _y VB _z -ΔN ¹ … -ΔN ^m Δφ _UPD ] ^T (3) Among them, Z _k is the residual vector of the observation value, that is, the single-difference observation minus the initial geometric distance, the initial ambiguity and the initial UPD value; X _k is the state vector; [dB _x dB _y dB _z ] ^T represents the baseline vector increment; [VB _x VB _y VB _z ] ^T is the baseline velocity vector increment; ΔN ⁱ is the ambiguity parameter increment of the i-th satellite; H _k is the design matrix, l _x , _ly , and l _z are the satellite-receiver connection Cosine of the direction of the line; (X ⁱ , Y ⁱ , Z ⁱ ) are satellite coordinates; (X ⁰ , Y ⁰ , Z ⁰ ) are the approximate coordinates of the receiver; ρ ⁰ is the distance between the satellite and the approximate coordinates of the receiver ε _k is the observation noise; the superscript T represents the transpose of the matrix; From epoch k to epoch k+1, the prediction equation of the state vector and its covariance matrix are: X _k+1,k =A _k X _k (6) The update formula of the state vector and its covariance matrix are: X _k+1 =X _k+1,k +K _k+1 (Z _k+1 -A _k+1 X _k+1,k ) (8) C _k+1 =C _k+1,k -K _k+1 A _k+1 C _k+1,k (9) Among them, A _k is the state transition matrix, Q _k is the error covariance matrix of the state transition equation, P _k+1 is the variance matrix of the observation noise, K _k+1 is the gain matrix, R _k is the variance matrix of the observation noise, C _{k+1, k} and C _k+1 are the prediction variance matrix and the update variance matrix, respectively; The floating-point solution Ni of the ambiguity parameter can be obtained by Kalman filtering, and the method of fixing the floating-point ambiguity is as follows: select a reference satellite ⁱ , and combine the floating-point ambiguity ^Ni and UPD parameters of the satellite into a new UPD parameter ; At this time, the new UPD parameter is actually the sum of UPD and the floating-point ambiguity of the reference satellite, while the floating-point ambiguity of other non-reference satellites becomes the difference between the original floating-point ambiguity and the floating-point ambiguity of the reference satellite; as follows The formula shows: in, is the new UPD parameter obtained by merging; Ni is the floating-point ambiguity of the reference satellite ⁱ , is the new floating-point ambiguity of satellite k after the UPD parameters and the reference satellite floating-point ambiguity are combined; the ambiguity is fixed through three steps: first, ignoring the integer characteristics of the ambiguity, using the parameter estimation method to solve the baseline vector and the ambiguity Floating-point solution; then use the ambiguity floating-point solution and its covariance matrix to fix the ambiguity satellite by satellite according to its integer characteristics; re-solve the integer ambiguity back to the original observation equation to obtain the baseline vector solution after the ambiguity is fixed ; (3) Calculate the multipath based on the obtained solution The posterior residual of the Kalmam filter solution is the multipath of the satellite: Step 2: Establishment of a multi-path hemisphere model based on trend surface analysis (1) Assign multipath data to the sky grid Calculate the altitude and azimuth for the satellites that have obtained multipath data; first, calculate the observation vector [Δe Δn Δu] ^T from the user to the satellite in the station center coordinate system: Among them, [Δx Δy Δz] ^T is the observation vector from the user to the satellite in the geocentric geo-fixed rectangular coordinate system; λ and φ are the longitude and latitude, respectively; Next, calculate the observation vector from the user to the satellite in the carrier coordinate system [Δx′ Δy′ Δz ′] ^T : Among them, ψ, θ and γ are the carrier heading angle, pitch angle and roll angle respectively; the satellite altitude angle ele and azimuth angle azi are: Divide the sky into different grids according to the altitude and azimuth, and assign the satellite multipath data corresponding to the altitude and azimuth to the grid; (2) Build a multi-path trend surface model for each grid of the sky First, model the primary and secondary trend surfaces for the multipath of the sky grid: A trend surface: Secondary trend surface: where x _i , y _i represent the azimuth and elevation angles corresponding to the i- _th multipath value in the sky grid, a ₀ , a ₁ ,...an represent the trend surface fitting coefficients, is the fitting value of the trend surface model; Secondly, carry out the R ² test and the significant F test for the appropriate fit; in which the coefficient of fit R ² is represented by the proportion of the regression square sum SSR to the total deviation square sum _{S T} ; the value of _R ² ranges from 0 to 1. The larger the value, the higher the fitting degree of the model; its calculation formula is: where _zi is is the residual sum of squares If R ² is greater than a predetermined threshold, the model is judged to pass the test; If both models pass the test and then use the F test to select the best, the F test is used to test whether the trend surface equation is significant, and its calculation formula is: where p represents the degrees of freedom of the regression sum of squares, given by Calculated, (np-1) is the degree of freedom of the remaining sum of squares, r is the number of trend surface analysis, n is the total number of residual values in the grid point; F obeys the F distribution with degrees of freedom (p, np-1), In the case of a given confidence level, by checking the F distribution table, the critical value F is obtained, so as to test which trend surface equation fitted is more significant under the confidence level α; If none of them pass, the mean value of residuals in the grid is used to replace the trend surface model directly; Step 3: Multipath correction in real time Collect GNSS observation data in real time, calculate the altitude and azimuth angles of satellites in the carrier coordinate system, select the multi-path trend surface model corresponding to the sky map grid for multi-path estimation, and correct the satellite observation data.