GNSS diffraction error elimination method based on regular hexagon grid

Abstract

The application provides a GNSS diffraction error elimination method based on a regular hexagon grid, comprising: acquiring observation data and broadcast ephemeris of GNSS reference stations and monitoring stations, calculating observation residual sequences of each satellite, and acquiring corresponding satellite elevation angles and azimuth angles; constructing a regular hexagon grid model covering an observation hemisphere, the regular hexagon grid model comprising a plurality of regular hexagon units with equal areas and isotropy, mapping the observation residual sequences into corresponding regular hexagon units according to the satellite elevation angles and azimuth angles, and forming a grid residual set; statistically analyzing the grid residual set in each regular hexagon unit, and calculating statistical characteristic parameters of the residuals; when the statistical characteristic parameters exceed a preset threshold, determining that the corresponding regular hexagon unit is contaminated by diffraction errors, generating a data quality control grid mask, and eliminating the observation values in the regular hexagon unit determined to be contaminated by diffraction errors in subsequent positioning calculation.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of satellite navigation and positioning, and particularly relates to a GNSS diffraction error elimination method based on a regular hexagonal grid. BACKGROUND

[0002] The Beidou / Global Navigation Satellite System (GNSS) high-precision positioning technology has become an important means for structural vibration and deformation monitoring due to its real-time, all-weather, high automation, and millimeter-level precision. However, in actual monitoring applications, Beidou / GNSS devices need to be placed on specific monitoring points, and there are inevitably environmental obstructions in the observation field of view. These obstructions will interfere with Beidou / GNSS signals, causing multipath effect and diffraction effect errors, which seriously affect the positioning accuracy and reliability, thereby reducing the usability of Beidou / GNSS in structural health monitoring.

[0003] To solve the above problems, existing research has proposed error mapping methods based on spatial grids, such as the Multipath Hemisphere Map (MHM) method. This method maps the observation residuals according to the elevation and azimuth angles of the satellites to a fixed rectangular grid, and averages the residuals within the grid to establish a site-related error correction model, thereby modeling and correcting the Beidou / GNSS positioning errors.

[0004] However, this method based on fixed rectangular grids has many defects. The rectangular grid structure itself has inherent deficiencies, the grid cell area will decrease sharply with the increase of the elevation angle, leading to sparse spatial sampling in the high-elevation area, making the modeling results unreliable; the anisotropic characteristics of the rectangular grid will also introduce directional bias, affecting the stability of parameter estimation; and the complex node adjacency relationship restricts the interpolation accuracy and smoothness. In addition, the traditional method is generally based on the assumption of zero mean residual, which makes it difficult to effectively identify and separate the diffraction errors that exhibit non-Gaussian distribution in statistical characteristics, resulting in distorted model. SUMMARY

[0005] The present application proposes a GNSS diffraction error elimination method based on a regular hexagonal grid, which solves the problem that the existing error mapping method based on a fixed rectangular grid cannot effectively identify and separate complex error patterns with non-Gaussian distribution, resulting in insufficient modeling accuracy and reliability.

[0006] To solve the above technical problems, the present application provides a GNSS diffraction error elimination method based on a regular hexagonal grid, comprising the following steps:

[0007] Step S1: Obtain the observation data and broadcast ephemeris of the Global Navigation Satellite System (GNSS) reference station and monitoring station, and calculate the observation residual sequence of each satellite to obtain the corresponding satellite elevation and azimuth angles;

[0008] Step S2: Constructing a regular hexagon grid model covering the observation hemispherical surface, the regular hexagon grid model comprising a plurality of regular hexagon cells with equal area and isotropy, mapping the observation residual sequence into the corresponding regular hexagon cell according to the satellite elevation angle and azimuth angle, forming a grid residual set;

[0009] Step S3: Performing statistical distribution analysis on the grid residual set in each regular hexagon cell, calculating the statistical characteristic parameters of the residual;

[0010] Step S4: When the statistical characteristic parameters exceed the preset threshold, determining that the corresponding regular hexagon cell is contaminated by diffraction error, generating a data quality control grid mask, and excluding the observation values in the regular hexagon cell determined to be contaminated by diffraction error in subsequent positioning calculation.

[0011] Preferably, the observation residual sequence of each satellite is calculated by using a double-difference relative positioning algorithm in step S1, and the corresponding satellite elevation angle and azimuth angle are obtained, including the following steps:

[0012] Step S11: Establishing a double-difference observation equation including multipath effect and diffraction error parameters:

[0013] ;

[0014] In the formula, is a double-difference operator; is the double-difference pseudo-range observation value of satellite , to the reference station and the monitoring station ; is the double-difference carrier phase observation value of satellite , to the reference station and the monitoring station ; is the wavelength of the carrier phase observation value of the mth frequency; is the double-difference geometric distance from the satellite to the receiver; is the speed of light in vacuum; is the double-difference carrier phase integer ambiguity parameter; , are the double-difference pseudo-range observation noise and double-difference carrier phase observation noise containing multipath effect and diffraction error information, respectively;

[0015] Step S12: Estimating the position parameters and ambiguity parameters of each satellite using Kalman filtering every epoch:

[0016] ;

[0017] wherein, is the pseudo-range and carrier phase observation at epoch ; is the design matrix at epoch ; is the state vector containing position parameters and ambiguity parameters at epoch ; is the residual vector at epoch ; denotes the coefficient matrix of the state transition equation between epoch and epoch ; is the normal white noise with mean zero and covariance matrix ;

[0018] Step S13: solving the formula in step S12 to obtain the float solution of each ambiguity parameter and the standard deviation , calculating the ambiguity fixing success rate of the ambiguity parameter, sorting the ambiguity parameters in descending order of the ambiguity fixing success rate, and sequentially substituting the ambiguity parameters into the ambiguity fixing criterion in order, and fixing the current ambiguity parameter if the ambiguity parameter meets the ambiguity fixing criterion; the expressions of the ambiguity fixing success rate and the ambiguity fixing criterion are:

[0019] ;

[0020] ;

[0021] ;

[0022] ;

[0023] wherein, is the ambiguity fixing success rate; is the float ambiguity closest to the integer; is the ambiguity fixing criterion; is the complementary error function; is the summation index variable; is the threshold value of the ambiguity fixing criterion;

[0024] Step S14: substituting the fixed ambiguity parameter into the formula in step S12 to obtain the carrier phase double difference residual of each satellite at each epoch, forming the observation residual sequence, and outputting the satellite elevation angle and azimuth angle.

[0025] Preferably, the step S2 adopts a layer-vertex-edge coding mode to construct a regular hexagonal grid model covering the observed hemispherical surface, comprising the following steps:

[0026] Step S201: defining the regular hexagonal unit at the center of the hemispherical surface as the first layer, and the layer coding is 1, and the number of layers extending outward is sequentially increased;

[0027] Step S202: defining the vertex coding of each layer, starting from the regular hexagonal unit at the top left corner of the current layer and coding clockwise;

[0028] Step S203: defining the edge coding of each layer, for the regular hexagonal unit that is not a vertex, the regular hexagonal unit is attributed to the nearest vertex in the counterclockwise direction, and the edge coding is performed clockwise within the current vertex group.

[0029] Preferably, in the step S2, the observed residual sequence is mapped into the corresponding regular hexagonal unit, comprising the following steps: calculating the geometric distance from the sky projection point corresponding to the observed residual to the center point of each adjacent regular hexagonal unit; and when and only when the geometric distance is less than or equal to the circumscribed circle radius of the regular hexagonal unit, and the geometric distance is the minimum value among all adjacent regular hexagonal units, the current observed residual is assigned to the corresponding regular hexagonal unit.

[0030] Preferably, the polar coordinates of the center point of the regular hexagonal unit are calculated, comprising the following steps:

[0031] Step S211: establishing a Cartesian coordinate system with the origin of the center point of the regular hexagonal unit, and defining a reference point :

[0032] ;

[0033] In the formula, , are the horizontal and vertical coordinates of ; is the layer coding number of the current regular hexagonal unit; is the side length of the regular hexagon;

[0034] Step S212: using the six-fold rotational symmetry of the regular hexagonal grid, establishing a reference point coordinate system, for the regular hexagonal unit with layer coding , vertex coding , and edge coding , the polar coordinate radius of the center point is:

[0035] ;

[0036] ;

[0037] ;

[0038] In the above formula, is an interpolation coefficient; is a normalization function, which describes the variation of the radius of the center point relative to the maximum radius on the edge of the hexagonal grid; is the length of the edge of the regular hexagonal unit;

[0039] Step S213: In combination with the polar coordinate radius and angle formula, the polar coordinates of the center point of the regular hexagonal unit are obtained as follows:

[0040] ;

[0041] ;

[0042] ;

[0043] ;

[0044] In the formula, is the angle of the center point of the regular hexagonal unit; is the base angle of the vertex ; is the offset angle of the base angle.

[0045] Preferably, in step S3, the probability density curve of the residual is fitted by using a kernel density estimation method, and the statistical characteristic parameters are calculated, and the expression of the probability density curve of the residual fitted by using the kernel density estimation method is as follows:

[0046] ;

[0047] ;

[0048] ;

[0049] In the formula, is a kernel density estimation formula; is the sample size; is the width of the kernel function; is the kernel function; denotes the independent variable of the probability density function; denotes the i-th residual sample value; is the sample standard deviation.

[0050] Preferably, the statistical characteristic parameters in step S3 include a relative intensity of a secondary peak and a bimodal separation degree for characterizing the bimodal distribution characteristics, and the calculation of the relative intensity of the secondary peak and the bimodal separation degree includes the following steps: identifying local maximum points in the probability density curve, determining a main peak height if the number of detected peak values is not less than 2 and the corresponding positions , and a secondary peak height and the corresponding positions ;

[0051] The relative intensity of the secondary peak is calculated as follows:

[0052] ;

[0053] The bimodal separation degree is calculated as follows:

[0054] .

[0055] Preferably, the statistical characteristic parameters include a heavy tail ratio for characterizing the degree of tailing , and the expression of the heavy tail ratio is as follows:

[0056] ;

[0057] ;

[0058] ;

[0059] ;

[0060] In the formula, Q is a quartile range; , P25 and P75 are the 25th percentile and the 75th percentile, respectively; ,is a median absolute deviation; indicates taking the median; represents the center position of the residual distribution; is a residual vector, the residuals in the residual vector are arranged in ascending order; is the total number of observation values; is the kth residual value after sorting. Preferably, the statistical characteristic parameters in step S4 at least include a skewness for characterizing the asymmetry of the distribution, a heavy tail ratio for characterizing the degree of tailing, a relative intensity of a secondary peak and a bimodal separation degree for characterizing the bimodal distribution characteristics; and a diffraction error decision rule is established according to the comparison result of the statistical characteristic parameters and a preset threshold, and the diffraction error decision rule includes:

[0061] Preferably, the statistical characteristic parameters in step S3 include a relative intensity of a secondary peak and a bimodal separation degree for characterizing the bimodal distribution characteristics, and the calculation of the relative intensity of the secondary peak and the bimodal separation degree includes the following steps: identifying local maximum points in the probability density curve, determining a main peak height if the number of detected peak values is not less than 2

[0062] (1) the absolute value of the skewness is greater than 0.34;

[0063] (2) the heavy tail ratio is greater than 1.6;

[0064] (3) the relative intensity of the secondary peak is greater than 0.24, and the bimodal separation degree is greater than 2.7 times the median absolute deviation;

[0065] (4) the shift index is greater than 1.9, and the shift index is the ratio of the median of the residual error to the median absolute deviation.

[0066] Preferably, the generation of the data quality control grid mask in step S4 rejects the observation values in the regular hexagonal unit determined to be contaminated by diffraction error in subsequent positioning solution, comprising the following steps:

[0067] Step S41: generate the data quality control grid sequence of each monitoring station, mark the decision value of the regular hexagonal unit contaminated by diffraction error as 0, and mark the decision value of the regular hexagonal unit not contaminated as 1;

[0068] Step S42: in the real-time positioning process, the corresponding regular hexagonal unit decision value is retrieved according to the elevation angle and azimuth angle of the satellite at the current epoch, and if the decision value is 0, the observation value of the satellite is rejected for positioning solution.

[0069] The beneficial effects of the present application at least include:

[0070] 1. The regular hexagonal grid has isotropic characteristics, avoiding the directional deviation caused by anisotropy of the traditional rectangular grid, and the area of the regular hexagonal grid is equal, solving the problem of sparse sampling in the high elevation angle area of the rectangular grid, and improving the reliability and accuracy of modeling;

[0071] 2. The statistical characteristic parameters of the residual error are calculated by statistical distribution analysis and compared with the preset threshold value, which can effectively identify and separate complex diffraction errors of non-Gaussian distribution, overcome the defects of traditional methods based on the assumption of zero mean value of residual error, and improve the accuracy of error rejection and the adaptability of the model;

[0072] 3. The data quality control grid mask is generated, which can reject the observation values contaminated by diffraction error in subsequent positioning solution, further optimize the positioning result, improve the positioning accuracy and reliability, and is especially suitable for structural health monitoring in complex shielding environment. BRIEF DESCRIPTION OF DRAWINGS

[0073] Figure 1 The method flowchart of the embodiment of the present application is shown in the figure;

[0074] Figure 2 The coding diagram of the regular hexagonal unit grid in the embodiment of the present application is shown in the figure;

[0075] Figure 3 This is a schematic diagram of the station distribution and observation environment according to an embodiment of the present invention;

[0076] Figure 4 This is a residual sky map grid model of the monitoring station in an embodiment of the present invention;

[0077] Figure 5 The diffraction error identification residual sky grid model of the monitoring station in this embodiment of the invention;

[0078] Figure 6 This is a sky map showing the data quality control of the monitoring station in an embodiment of the present invention. Detailed Implementation

[0079] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the protection scope of the present invention.

[0080] like Figure 1 As shown, this embodiment of the invention provides a GNSS diffraction error elimination method based on a regular hexagonal grid, including the following steps:

[0081] Step S1: Obtain observation data and broadcast ephemeris from GNSS reference stations and monitoring stations of the Global Navigation Satellite System. Use the double-difference relative positioning algorithm to calculate the observation residual sequence of each satellite and obtain the corresponding satellite elevation angle and azimuth angle.

[0082] Specifically, for BeiDou / GNSS data from base stations and monitoring stations in practical positioning applications, 10 consecutive days of observation data and broadcast ephemeris files are selected. Then, a post-hoc double-difference relative positioning algorithm is used to solve the data, outputting the double-difference residual sequence, elevation angle, and azimuth angle for each satellite. This includes the following steps:

[0083] Step S11: For the base station and monitoring station Based on the carrier phase and pseudorange observations, construct the double-difference observation equation:

[0084] ;

[0085] In the formula, It is a double difference operator; For satellites at frequency m , For the base station and monitoring station The double-difference pseudorange observations; for the satellite , of the reference station and the monitoring station in units of cycles; for the carrier phase observation value of the mth frequency for the double-difference geometric distance from the satellite to the receiver; for the speed of light in vacuum; for the double-difference carrier phase integer ambiguity parameter in units of cycles; , respectively for the double-difference pseudo-range observation noise and the double-difference carrier phase observation noise, which contain the multipath effect and diffraction error information.

[0086] Step S12: write the above formula into a matrix form

[0087] ;

[0088] wherein, is the current epoch; is the observation value vector of the epoch , i.e. the pseudo-range and carrier phase observation values; is the design matrix of the epoch ; is the state vector of the epoch , which contains the position parameters and ambiguity parameters; is the residual vector of the epoch .

[0089] The following formula is constructed to build the standard Kalman filtering equation, and each position parameter (X, Y, Z), ambiguity parameter and its standard deviation is estimated epoch by epoch:

[0090] ;

[0091] wherein, represents the coefficient matrix of the state transition equation between the epoch and the epoch ; is a normal white noise with a mean of zero and a covariance matrix .

[0092] In the estimation process, the ambiguity of one satellite arc segment is set as one parameter for estimation, and after the parameter estimation, the floating point solution of each ambiguity parameter and its standard deviation are obtained.

[0093] Step S13: the obtained floating point ambiguity solution and its standard deviation The ambiguity fixing success rate is brought into the following formula, the ambiguity parameters are sorted in descending order of ambiguity fixing success rate, the ambiguity parameters are sequentially substituted into the ambiguity fixing criterion in order, and if the ambiguity parameters meet the ambiguity fixing criterion, the current ambiguity parameters are fixed. The expression of the ambiguity fixing success rate and the ambiguity fixing criterion is:

[0094] ;

[0095] ;

[0096] ;

[0097] ;

[0098] In the formula, is the ambiguity fixing success rate; is the floating-point ambiguity is the closest integer; is the ambiguity fixing criterion; is the complementary error function; is the summation index variable, which is a positive integer starting from 1; is the threshold value of the ambiguity fixing criterion, which usually represents the confidence level.

[0099] The formula calculates the reliability index of ambiguity fixing by accumulating the complementary error function values under different integer offsets When the calculated value is greater than or equal to , it is considered that the current ambiguity parameters meet the fixing criterion and can be fixed.

[0100] When the criterion exceeds 99.9%, it is considered that the ambiguity can be reliably fixed, where 99.9% is an empirical value commonly used for ambiguity fixing.

[0101] Step S14: After ambiguity fixing, the fixed ambiguity is substituted back into the Kalman filter equation for update calculation to obtain the monitoring station position parameters (X, Y, Z) after ambiguity fixing, and output the double difference residual sequence, signal-to-noise ratio sequence, satellite elevation angle and azimuth angle sequence of each satellite and each epoch.

[0102] Step S2: Construct a regular hexagonal grid model covering the observation hemisphere, the regular hexagonal grid model includes a plurality of regular hexagonal cells with equal area and isotropy, map the observation residual sequence into the corresponding regular hexagonal cell according to the satellite elevation angle and azimuth angle, and form a grid residual set.

[0103] Specifically, a regular hexagonal grid model covering the observation hemisphere is constructed using a layer-vertex-edge encoding method, as shown inFigure 2 The method comprises the following steps:

[0104] Step S201: defining the regular hexagonal element at the center of the hemisphere as the first layer, and the layer code is 1, and then the layer code of all grids of each layer is increased by 1;

[0105] Step S202: defining the vertex code of the first layer as 1, and then sequentially coding the vertices of the quasi-hexagon in a clockwise direction from the regular hexagon at the top left corner of each layer;

[0106] From the third layer, the quasi-hexagon appears a non-vertex edge, and the vertex code of the hexagon on the non-vertex edge is consistent with the vertex code of the vertex closest to it in the counterclockwise direction;

[0107] Step S203: defining the edge code at each layer vertex as 1;

[0108] From the third layer, a non-vertex edge appears, and all hexagons with the same vertex code are divided into a group, and the layer code is sequentially coded in a clockwise direction from the vertex.

[0109] Mapping the observation residual sequence into the corresponding regular hexagonal element, specifically comprising the following steps:

[0110] Step S211: calculating the center point of each regular hexagonal grid in polar coordinates;

[0111] Step S212: establishing a Cartesian coordinate system with the center point of the regular hexagonal grid as the origin O(0, 0), and for the nth layer (n = 1, 2, 3, …), defining the reference point , which corresponds to the first vertex (v = 1) of the current layer and is located on the first edge (e = 1):

[0112] ;

[0113] In the formula, , are the horizontal and vertical coordinates of , respectively; is the side length of the regular hexagon.

[0114] Step S213: using the six-fold rotational symmetry of the hexagonal grid, i.e., the grid structure coincides with itself after rotating 60° by an integer multiple around the origin, and the positions of other vertices are derived from the reference point by rotation transformation.

[0115] For any , the vertex code , and the edge code , the polar coordinate radius r of the hexagonal center point is: ​

[0116] ;

[0117] ;

[0118] ;

[0119] In the above formula, is an interpolation coefficient; is a normalization function, which describes the variation of the radius of the center point relative to the maximum radius on the edge of the hexagonal grid; is the edge length of the regular hexagonal unit;

[0120] For , the center point angle is:

[0121] ;

[0122] ;

[0123] ;

[0124] wherein, is the actual angle of the vertex ; is the offset angle of the base angle.

[0125] Combined with the radius and angle formulas, the complete polar coordinate expression formula of the center point of the hexagon is obtained:

[0126] .

[0127] Step S213: After obtaining the center point coordinates of each grid, the geometric distance from the sky projection point corresponding to the observation residual to the center points of the adjacent regular hexagonal units is calculated; only when the geometric distance is less than or equal to the circumscribed circle radius of the regular hexagonal unit, and the geometric distance is the minimum value among all adjacent regular hexagonal units, the current observation residual is assigned to the corresponding regular hexagonal unit:

[0128] ;

[0129] In the formula, is the geometric distance between the sky projection point corresponding to the jth observation residual and the center point of the ith regular hexagonal unit; is the circumscribed circle radius of the regular hexagonal unit; is the number of layers of the hexagon.

[0130] Step S3: Perform statistical distribution analysis on the grid residual set in each regular hexagon unit, fit the probability density curve of the residual using the kernel density estimation method, and calculate the statistical characteristic parameters; the statistical characteristic parameters at least include: skewness for representing the asymmetry of the distribution, heavy tail ratio for representing the degree of tailing, and secondary peak relative intensity and bimodal separation degree for representing the bimodal distribution characteristics. Specifically, the following steps are included:

[0131] Step S31: Perform kernel density estimation on the residual in each grid, and the kernel density estimation formula is:

[0132] ;

[0133] ;

[0134] ;

[0135] In the formula, is the kernel density estimation formula; is the sample size; is the bandwidth, which determines the width of the kernel function; is the kernel function; represents the independent variable of the probability density function, i.e. the point whose probability density value is to be estimated, which is a continuous variable representing the value that the residual can take; represents the i-th residual sample value, and the residual refers to the difference between the observed value and the model predicted value, and here the residual sample is the original data point used to estimate the probability density function; is the sample standard deviation.

[0136] Step S32: Based on the kernel density estimation result, identify the local maximum points in the density curve. To exclude the pseudo-peak caused by noise, set the minimum peak prominence to 5% of the maximum density value. If the number of detected peaks is not less than 2, further calculate the secondary peak relative intensity B and the bimodal separation degree to quantify the bimodal characteristics:

[0137] ;

[0138] In the formula, is the main peak height; is the secondary peak height; , are the positions of the main peak and the secondary peak respectively. The secondary peak relative intensity B reflects the height ratio of the secondary peak to the main peak, and is used to measure the relative significance of the secondary peak in the bimodal structure. A higher B value indicates that the secondary peak is more significant, and the bimodal structure is more obvious. The bimodal separation degree D represents the horizontal distance between the two peaks, and is used to measure the separation degree of the bimodal in the numerical range. A larger D value means that the bimodal is more dispersed in the numerical value, and the bimodal structure is clearer.

[0139] Step S33: Quantify the asymmetry and systematic bias of the residual distribution using skewness, median, and absolute deviation of the median for quantitative evaluation.

[0140] ;

[0141] In the formula, Skewness; represents the sample mean, which is the average value of the residual vector R; s represents the sample standard deviation, which is used to measure the dispersion of the residuals.

[0142] Skewness The value of can reflect the direction and degree of skewness of the residual distribution. When When positive, it indicates a right-skewed distribution, meaning the right tail of the data is longer; when... A negative skewness indicates a left-skewed distribution, meaning the data has a longer left tail. The larger the absolute value of the skewness, the more significant the asymmetry of the distribution.

[0143] Step S34: Use the median to measure the central location of the residual distribution, arrange the residuals in ascending order (R(1)≤R(2)≤⋯≤R(n), select the median value according to the parity of the sample size, and finally output the median med:

[0144] ;

[0145] In the formula, R is the residual vector, which contains n observations; This is the k-th residual value after sorting; This represents the total number of observations.

[0146] Compared to the mean, the median is more robust to extreme values ​​and less susceptible to outliers.

[0147] Step S35: Measure the dispersion of the residual data using the absolute deviation of the median (MAD):

[0148] ;

[0149] In the formula, The mathematical operator or calculation function for the median is: given a set of values, median means to take the median, that is, the value in the middle position when the data is sorted in ascending order. If the number of data is even, the average of the two middle numbers is taken. is the residual between the i-th observation and the model prediction; 1.4826 is a constant factor used to convert MAD into a measure with the same dimensions as the standard deviation.

[0150] MAD can effectively reflect the dispersion degree of data by calculating the median of the absolute difference of each data point from the median, while avoiding the influence of extreme values on the dispersion degree estimate. By calculating the above statistical quantities, this study can comprehensively judge the skew direction and degree of residual distribution.

[0151] Step S36: Quantify the systematic bias using the ratio of median to MAD :

[0152] ;

[0153] This index reflects the deviation of the central position from the dispersion degree, and the larger the value, the more significant the systematic bias.

[0154] Step S37: Quantify the heavy-tailed characteristics of the residual distribution by calculating the ratio of interquartile range (IQR) to MAD, is an index that describes the dispersion degree of the middle 50% of observations in the data, defined as:

[0155] ;

[0156] where, , are the 25th and 75th percentiles, respectively.

[0157] Calculate the ratio of to MAD, heavy-tailed ratio (T):

[0158] ;

[0159] This ratio is used to measure the degree of heavy tails of the data. If T is significantly greater than 1.6, it indicates that the residual distribution has heavy-tailed characteristics, i.e., the tail is thicker than the normal distribution, and the probability of extreme values is higher.

[0160] Step S4: Establish a diffraction error decision rule based on the comparison results of the statistical characteristic parameters and the preset threshold values, and when the statistical characteristic parameters satisfy the diffraction error decision rule, determine that the positive hexagonal unit is contaminated by diffraction error, generate a data quality control grid mask, and exclude the observation values in the positive hexagonal unit determined to be contaminated by diffraction error in subsequent positioning calculation.

[0161] Specifically, based on the above statistical characteristics, the following decision rule is established to identify diffraction error:

[0162] Condition 1: The absolute value of skewness is greater than 0.34;

[0163] Condition 2: Heavy-tailed ratio is greater than 1.6;

[0164] Condition 3: the secondary peak relative intensity is greater than 0.24, and the bimodal separation degree is greater than 2.7 times of the median absolute deviation;

[0165] Condition 4: the offset index is greater than 1.9, and the offset index is the ratio of the residual median to the median absolute deviation.

[0166] The absolute value of skewness exceeding 0.34 indicates that the distribution has obvious asymmetry, the heavy tail ratio exceeding 1.6 indicates that the tail of the distribution is significantly heavier than the normal distribution, the secondary peak relative intensity exceeding 0.24 and the bimodal separation degree being greater than 2.7 times of the MAD indicate that there is a significant bimodal structure, and the offset index exceeding 1.9 indicates that there is a significant systematic offset. If any of the above conditions is met, it is determined that the residual distribution has diffraction error.

[0167] After all the grid judgments are completed, a data quality control grid sequence of each monitoring station is generated, and the decision value of the regular hexagonal unit contaminated by the diffraction error is marked as 0, and the decision value of the regular hexagonal unit not contaminated is marked as 1.

[0168] In the real-time positioning process, the corresponding regular hexagonal unit decision value is retrieved according to the elevation angle and the azimuth angle of the satellite at the current epoch, and if the decision value is 0, the observation value of the satellite is refused to be used for positioning calculation.

[0169] The method of the embodiment of the application adopts a GNSS post-difference relative positioning algorithm for positioning calculation, can use the observation value of the entire ambiguity arc segment to fix the ambiguity, is beneficial to retaining the observation value residual affected by the shielding, and outputs a double-difference residual sequence and an elevation angle and azimuth angle sequence.

[0170] The constructed regular hexagonal unit grid has more consistent adjacency, that is, the 6 adjacent units are adjacent to the center grid with edges, and the distance from the center grid is the same. This feature makes the hexagon have certain advantages in processing the nearest neighbor, moving path and other neighborhood processing problems. Under the condition of the same area, the hexagon is closer to the circle, and the grid structure is more compact. Compared with the quadrilateral structure, under the same data amount, the hexagon has higher data precision. Moreover, the hexagon is more isotropic, and has more advantages in spatial field modeling.

[0171] In the field of deformation monitoring, the baseline of Beidou / GNSS is generally short, and the spatial distribution characteristics of the residuals of each station are similar under the condition of no obstruction. In the area affected by multipath effect and diffraction effect, the residuals in the obstructed area will show non-normal distribution characteristics. According to this characteristic, the method proposes to construct a spatial grid to calculate the statistical information of the residuals and determine a set of residual spatial distribution models in an ideal state. Since the obstructed positions of different stations may be different, this method can better reflect the signal-to-noise ratio spatial distribution model in the unobstructed environment facing the area. In the data processing process, only the observation values in the area not disturbed by obstruction are used for positioning calculation, and high-precision and high-reliability positioning results will be obtained.

[0172] Embodiment one

[0173] The test data is collected from a real running slope monitoring system of a stone field. As shown in Figure 3 , a reference station (JZ01) is set in the system, located on the roof of the data center building, and 10 monitoring stations (WY01-WY10) are respectively located on the mining section of the stone mine.

[0174] Each station is equipped with a NET10 Plus receiver and a matching antenna of ComSp. The receiver is set to receive GPS / BDS-2 / BDS-3 / GALILEO satellite system observation data, and the sampling frequency is set to 5s. The data is stored and calculated every 4 hours as a time period. In the project, 4G wireless network transmission is used for data transmission, which is directly transmitted to the corresponding port of the server through the external network port mapping, and the GNSS data management software is responsible for receiving and storing. The present embodiment selects a representative WY08 station as the test data, and the station information is shown in Figure 4 .

[0175] The test data used in this embodiment is the Beidou-2 / Beidou-3 / GPS / Galileo observation data from May 17, 2022 to June 5, 2022, with a sampling frequency of 5s. Figure 5 The diffraction error identification residual sky grid model for the WY08 monitoring station can be seen from the figure, and the residuals show a smooth decreasing characteristic from high elevation angle to low elevation angle, and there is no obvious signal-to-noise ratio mutation phenomenon caused by obstruction. Figure 4 The residual sky grid model for the WY08 monitoring station is shown in Figure 6 The available data area generated after diffraction error identification and elimination of the WY08 station measured residual sky grid is shown in the figure. It can be seen from the figure that the edge part of the residuals is obviously lower than the high-quality residual grid model, and the processed Figure 5 Figure 6 ​The available data area of the color region is divided into a plurality of regular hexagonal grids, and a part of the area is excluded by the slope shelter. In the data processing, if the observation value is located in the color region, the observation value is involved in the calculation, and if the observation value is located in the blank region, the observation value is not involved in the calculation. It can be seen that the GNSS diffraction error detection and elimination method of the regular hexagonal grid division effectively solves the problem of data quality decline caused by the shelter.

[0176] The technical features of the above embodiments can be combined in any manner. In order to make the description simple, all possible combinations of the technical features in the above embodiments are not described, only the preferred embodiments of the present application are expressed, and the description is more specific and detailed, but it should not be understood as a limitation on the scope of the present application. As long as the combination of these technical features does not exist, it should be considered as the scope of the present application.

[0177] It should be noted that, for those skilled in the art, without departing from the concept of the present application, a number of modifications and improvements can be made, which are within the scope of the present application. Therefore, the protection scope of the present application should be subject to the appended claims.

Claims

1. A GNSS diffraction error elimination method based on a regular hexagonal lattice, characterized in that, Includes the following steps: Step S1: Obtain observation data and broadcast ephemeris from GNSS reference stations and monitoring stations of the Global Navigation Satellite System, calculate the observation residual sequence of each satellite, and obtain the corresponding satellite elevation angle and azimuth angle; Step S2: Construct a regular hexagonal grid model covering the observation hemisphere. The regular hexagonal grid model includes multiple regular hexagonal cells with equal area and isotropic properties. Based on the satellite elevation angle and azimuth angle, map the observation residual sequence to the corresponding regular hexagonal cells. Calculate the geometric distance from the sky projection point corresponding to the observation residual to the center point of each adjacent regular hexagonal cell. If and only if the geometric distance is less than or equal to the radius of the circumcircle of the regular hexagonal cell, and the geometric distance is the minimum value among all adjacent regular hexagonal cells, assign the current observation residual to the corresponding regular hexagonal cell to form a grid residual set. The calculation of the polar coordinates of the center point of the regular hexagonal unit includes the following steps: Calculating the polar coordinates of the center point of the regular hexagonal unit involves calculating the geometric distance from the sky projection point corresponding to the observation residual to the center points of adjacent regular hexagonal units. Step S211: Establish a Cartesian coordinate system with the origin of the center point of the regular hexagonal unit, and define a reference point in the Cartesian coordinate system. : ； In the formula, , They are respectively The x and y coordinates; This represents the layer code number of the current regular hexagonal cell; The side length is that of a regular hexagon; Step S212: Utilize the six-fold rotational symmetry of the regular hexagonal mesh to establish a reference point coordinate system. For the layer encoding as... Vertex encoding is The edge is encoded as The polar coordinate radius of the center point of the regular hexagonal element. for: ； ； ； In the above formula, These are the interpolation coefficients; As a normalization function, it describes how the radius of the center point varies relative to the maximum radius on the edges of a hexagonal grid. Let be the side length of the regular hexagonal unit; Step S213: Combining the polar coordinate radius and angle formulas, the polar coordinates of the center point of the regular hexagonal unit are obtained as follows: ； ； ； ； In the formula, The angle of the center point of the regular hexagonal unit; As vertices The basic perspective; The offset angle from the base angle; Step S3: Perform statistical distribution analysis on the set of grid residuals within each regular hexagonal cell, and calculate the statistical characteristic parameters of the residuals; Step S4: When the statistical feature parameter exceeds the preset threshold, it is determined that the corresponding regular hexagonal cell is contaminated by diffraction error. A data quality control grid mask is generated, and the observation values ​​in the regular hexagonal cells that are determined to be contaminated by diffraction error are removed in the subsequent positioning calculation.

2. The GNSS diffraction error elimination method based on a regular hexagonal grid according to claim 1, characterized in that: Step S1 uses a double-difference relative positioning algorithm to calculate the observation residual sequence of each satellite and obtain the corresponding satellite elevation angle and azimuth angle, including the following steps: Step S11: Establish the double-difference observation equation, including multipath effects and diffraction error parameters: ； In the formula, It is a double difference operator; For frequency m Up to satellite , For the base station and monitoring station The double-difference pseudorange observations; For frequency m Up to satellite , For the base station and monitoring station The double-difference carrier phase observations; For the first m Wavelength of carrier phase observation at each frequency; This represents the geometric distance from the satellite to the receiver after double difference. The speed of light in a vacuum; The parameter is the integer ambiguity parameter for the double-difference carrier phase. , These are, respectively, the double-difference pseudorange observation noise and the double-difference carrier phase observation noise, which include multipath effects and diffraction error information; Step S12: Estimate the position parameters and ambiguity parameters of each satellite epoch-by-epoch using Kalman filtering: ； In the formula, For the calendar The observed values ​​of pseudorange and carrier phase; For the calendar The design matrix; For the calendar A state vector containing position parameters and ambiguity parameters; For the calendar The residual vector; Represents the epoch He Liyuan The coefficient matrix of the state transition equations between them; The mean is zero and the covariance matrix is Normal white noise; Step S13: Solve the formula in step S12 to obtain the floating-point solutions for each ambiguity parameter. and standard deviation The ambiguity fixation success rate of the ambiguity parameters is calculated. The ambiguity parameters are then sorted in descending order of their fixation success rate. Each ambiguity parameter is sequentially substituted into the ambiguity fixation criterion. If a ambiguity parameter satisfies the ambiguity fixation criterion, the current ambiguity parameter is fixed. The expressions for the ambiguity fixation success rate and the ambiguity fixation criterion are: ； ； ； ； In the formula, Fixed success rate for ambiguity; For floating-point ambiguity The closest integer; Establish a fixed criterion for ambiguity; It is a complementary error function; For summation index variables; The threshold for a fixed ambiguity criterion; Step S14: Substitute the fixed ambiguity parameters back into the formula in step S12 to obtain the carrier phase double difference residuals of each satellite at each epoch, form the observation residual sequence, and output the satellite elevation angle and azimuth angle.

3. The GNSS diffraction error elimination method based on a regular hexagonal grid according to claim 1, characterized in that: Step S2 uses a layer-vertex-edge encoding method to construct a regular hexagonal grid model covering the observed hemisphere, including the following steps: Step S201: Define the regular hexagonal unit at the center of the hemispherical surface as the first layer, with the layer code 1, and the number of layers extending outwards increases sequentially; Step S202: Define the vertex encoding for each layer, starting from the hexagonal cell at the top left of the current layer and encoding clockwise; Step S203: Define the edge encoding for each layer. For non-vertex regular hexagonal cells, assign the non-vertex regular hexagonal cells to the nearest vertex in the counterclockwise direction, and perform edge encoding clockwise within the current vertex group.

4. The GNSS diffraction error elimination method based on a regular hexagonal grid according to claim 1, characterized in that: In step S3, the probability density curve of the residuals is fitted using the kernel density estimation method, and statistical characteristic parameters are calculated. The expression for the probability density curve of the residuals fitted using the kernel density estimation method is as follows: ； ； ； In the formula, Here is the formula for estimating kernel density; For sample size; The width of the kernel function; For kernel functions; The independent variable represents the probability density function; Indicates the first i One residual sample value; This represents the sample standard deviation.

5. The GNSS diffraction error elimination method based on a regular hexagonal grid according to claim 4, characterized in that: The statistical characteristic parameters mentioned in step S3 include the relative intensity of the secondary peak and the separation of the two peaks, which are used to characterize the bimodal distribution. Calculating the relative intensity of the secondary peak and the separation of the two peaks includes the following steps: identifying local maxima points in the probability density curve; if the number of detected peaks is not less than two, then determining the height of the primary peak. and corresponding positions and the height of the second peak and corresponding positions ; Calculate the relative intensity of the secondary peak: ； Calculate the bimodal resolution: 。 6. The GNSS diffraction error elimination method based on a regular hexagonal grid according to claim 4, characterized in that: The statistical characteristic parameters include the heavy-tail ratio, which characterizes the degree of tailing. The heavy-tail ratio The expression is: ； ； ； ； In the formula, Interquartile range; , These are the 25th and 75th percentiles, respectively; This represents the absolute deviation of the median. This indicates taking the median; This represents the central location of the residual distribution; For the residual vector, The residuals are arranged in ascending order; The total number of observations; For the sorted number k Each residual value.

7. The GNSS diffraction error elimination method based on a regular hexagonal grid according to claim 6, characterized in that: The statistical characteristic parameters mentioned in step S4 include at least: skewness for characterizing distribution asymmetry, heavy-tail ratio for characterizing the degree of tailing, relative intensity of the secondary peaks and separation of the bimodal distribution for characterizing the bimodal distribution characteristics; based on the comparison results of the statistical characteristic parameters with preset thresholds, a diffraction error decision rule is established, the diffraction error decision rule including: (1) The absolute value of the skewness is greater than 0.34; (2) The heavy-tail ratio is greater than 1.6; (3) The relative intensity of the secondary peak is greater than 0.24, and the separation of the two peaks is greater than 2.7 times the absolute deviation of the median; (4) The offset index is greater than 1.9, where the offset index is the ratio of the median of the residuals to the absolute deviation of the median.

8. The GNSS diffraction error elimination method based on a regular hexagonal grid according to claim 1, characterized in that: The step S4, which generates a data quality control grid mask to remove observations from hexagonal cells deemed to be contaminated by diffraction errors in subsequent localization calculations, includes the following steps: Step S41: Generate the data quality control grid sequence for each monitoring station, mark the decision value of the regular hexagonal cell contaminated by diffraction error as 0, and mark the decision value of the uncontaminated regular hexagonal cell as 1; Step S42: During real-time positioning, retrieve the corresponding hexagonal cell decision value based on the elevation angle and azimuth angle of the current epoch satellite. If the decision value is 0, refuse to use the satellite's observation value for positioning calculation.

Images