A Beidou satellite signal real-time dynamic positioning method, system, device and medium

Abstract

The application discloses a Beidou satellite signal real-time dynamic positioning method, system, device and medium based on carrier phase observations, and belongs to the technical field of Beidou satellite signal positioning, and comprises the following steps: identifying the observation values polluted by multipath from original carrier phase observations, extracting a first multipath characteristic component and a second multipath characteristic component from the polluted observation values of the first type of satellites and the second type of satellites, and performing error compensation on the original carrier phase observations based on the two types of characteristic components to obtain a first position coordinate; when the residual error calculated based on the first position coordinate presents a preset feature, the two types of characteristic components are corrected by using the observation data of the second type of satellites in a preset space region, and error compensation is performed based on the corrected characteristic components to obtain a second position coordinate. The application realizes the cooperative calibration of multipath characteristic components by utilizing the spatial complementarity of the Beidou mixed constellation.

Description

Technical Field

[0001] This invention relates to the field of BeiDou satellite signal positioning technology, and in particular to a method, system, device and medium for real-time dynamic positioning of BeiDou satellite signals. Background Technology

[0002] The BeiDou satellite system is a hybrid constellation consisting of geostationary orbit satellites and medium Earth orbit satellites, providing carrier phase observations for real-time dynamic positioning. These carrier phase observations are susceptible to multipath error, caused by environmental reflections around the receiver, which degrades positioning accuracy.

[0003] Geostationary orbit satellites maintain a quasi-stationary position during the observation period, exhibiting a constant bias in their multipath error. Medium Earth orbit satellites, however, remain in motion during the observation period, with their multipath error varying with their elevation angle. Current technologies process the multipath characteristics of satellites in different orbit types independently, lacking a collaborative calibration mechanism between them. When systematic biases exist in the extracted multipath feature components, the complementarity of hybrid constellations cannot be used to verify and correct these biases, limiting the effectiveness of error compensation.

[0004] Therefore, a method, system, equipment, and medium for real-time dynamic positioning using BeiDou satellite signals are needed. Summary of the Invention

[0005] The technical problem solved by this invention is: how to achieve coordinated calibration of multi-path characteristic components of satellites of different orbit types by utilizing the complementarity of geostationary orbit satellites and medium Earth orbit satellites in the BeiDou hybrid constellation in terms of space coverage, and at the same time, improve the real-time dynamic positioning accuracy by directionally correcting systematic deviations based on the residual tracing mechanism.

[0006] In view of the above-mentioned problems, the present invention provides a method, system, device and medium for real-time dynamic positioning of BeiDou satellite signals based on carrier phase observations, the specific technical solution of which is as follows: A real-time dynamic positioning method for BeiDou satellite signals includes the following steps: The observations contaminated by multipath are identified from the original carrier phase observations. The first multipath feature component is extracted from the contaminated observations of the first type of satellites, and the second multipath feature component is extracted from the contaminated observations of the second type of satellites. The original carrier phase observations are then error-compensated based on the first and second multipath feature components, and the first position coordinates are calculated. The residual is calculated based on the first position coordinates. When the residual exhibits a preset characteristic, the observation data of the second type of satellite in the preset space region is extracted from the contaminated observations of the second type of satellite. The first multipath feature component is corrected based on the observation data. The second multipath feature component is corrected based on the corrected first multipath feature component. The original carrier phase observation is error-compensated based on the corrected feature component and the second position coordinates are obtained and output. When the residual does not exhibit the preset feature, output the first position coordinates; Store the output position coordinates. If the distance between the current position and the stored position is less than a preset threshold, output the stored position coordinates.

[0007] Preferably, identifying multipath-contaminated observations from the original carrier phase observations includes: For each visible satellite, carrier phase observations at the first and second frequency points are selected and linearly combined to form inter-frequency combined observations; Collect the values ​​of the frequency combination observations over a continuous observation period, calculate the difference between the frequency combination observations of adjacent epochs, and obtain the time-domain rate of change sequence; Based on the relationship between the time-domain rate of change sequence and frequency, a dispersion relation test sequence is obtained; Determine whether there is a jump in the values ​​of the dispersion relation test sequence in adjacent epochs. If a jump occurs, it is determined that the satellite observation value is contaminated by multipath.

[0008] Preferably, the extraction of the first multipath feature component from the contaminated observations of the first type of satellite includes: The carrier phase observations of the first type of satellites are sequentially processed by first-order time difference processing, second-order time difference processing, and time domain accumulation processing. The center position value is extracted from the accumulation result as the first multipath feature component. The extraction of the second multipath feature component from contaminated observations of the second type of satellites includes: The carrier phase observations of the second type of satellite are processed by first-order time difference to obtain a second time difference sequence. Spline interpolation is performed on the representative values ​​of the elevation angle interval to establish a first spatial reference function. The first spatial reference function is subtracted from the second time difference sequence to obtain a first residual sequence. The coefficient parameters of the elevation angle polynomial function are obtained by least-squares fitting of the first residual sequence as the second multipath feature component.

[0009] Preferably, the step of performing error compensation on the original carrier phase observation value based on the first multipath feature component and the second multipath feature component and then calculating the first position coordinate includes: The observation space is divided into grids, and the number of first-type satellites and second-type satellites covered in each grid are counted. The coverage ratio is calculated as the first-type feature attribution degree and the second-type feature attribution degree. Based on the first type of feature attribution degree and the second type of feature attribution degree, the first multipath feature component and the second multipath feature component of each grid-covered satellite are weighted and combined to obtain the multipath error prediction value of the grid. For each visible satellite, the first compensated observation is formed by subtracting the multipath error prediction value of the corresponding grid from the original carrier phase observation value; Based on the first compensated observation, carrier phase differential positioning is calculated to obtain the first integer ambiguity and the first position coordinates.

[0010] Preferably, the extraction of observation data of the second type of satellite in a preset space region from the contaminated observations of the second type of satellite includes: The residual is mapped to a grid to form a grid residual distribution. When the difference between the residuals of adjacent grids is within a preset spatial correlation range, it is determined that the preset feature is presented. A first correlation line is obtained by linearly fitting the grid residuals and the first type of feature attribution degree. When the slope of the first correlation line exceeds a preset first slope threshold, grids with the first type of feature attribution degree exceeding the preset first attribution threshold are selected to form a first spatial direction set. A second correlation line is obtained by linearly fitting the grid residuals and the second type of feature attribution degree. When the slope of the second correlation line exceeds the preset second slope threshold, grids with the second type of feature attribution degree exceeding the preset second attribution threshold are selected to form a second spatial direction set. The intersection region of the first spatial direction set and the second spatial direction set is taken as the preset spatial region, and the carrier phase observation values ​​of the second type of satellite in the preset spatial region are extracted as the observation data.

[0011] Preferably, modifying the second multipath feature component based on the modified first multipath feature component includes: The observation data is processed by first-order time difference and the regional multipath feature value is statistically extracted. The difference between the regional multipath feature value and the statistical value of the first multipath feature component of the first type of satellite in the preset space region is calculated as the feature deviation value. The feature deviation value is superimposed on the first multipath feature component to obtain the corrected first multipath feature component. Substitute the center elevation angle of the preset spatial region into the elevation angle polynomial function to obtain the theoretical multipath feature value. Calculate the difference between the theoretical multipath feature value and the statistical value of the corrected first multipath feature component in the preset spatial region as the second type of feature adjustment amount. Add the second type of feature adjustment amount to the constant term coefficient of the elevation angle polynomial function to obtain the corrected second multipath feature component.

[0012] Preferably, the step of performing error compensation on the original carrier phase observation based on the corrected feature components, calculating and outputting the second position coordinates includes: Based on the corrected first multipath feature component and the corrected second multipath feature component, and combined with the first type of feature attribution degree and the second type of feature attribution degree, the corrected multipath error prediction value of each grid is calculated. For each visible satellite, the second compensated observation is formed by subtracting the corrected multipath error value of the corresponding grid from the original carrier phase observation value. Using the first integer ambiguity and the first position coordinates as initial estimates, the extended Kalman filter algorithm is used to iteratively solve the carrier phase double-difference observation equation established based on the second compensated observation value. The integer ambiguity is fixed and the second position coordinates are calculated and output.

[0013] A real-time dynamic positioning system for BeiDou satellite signals, employing the real-time dynamic positioning method for BeiDou satellite signals as described above, includes: A multipath identification module is used to identify observations contaminated by multipath from the original carrier phase observations; The feature extraction module is used to extract the first multipath feature component and the second multipath feature component from the contaminated observations of the first type of satellite and the second type of satellite, respectively. The first positioning module is used to perform error compensation on the original carrier phase observation value based on the first multipath feature component and the second multipath feature component and calculate the first position coordinates. The feature correction module is used to correct the first multipath feature component and the second multipath feature component from the observation data of the second type of satellite in the preset space region when the residual calculated based on the first position coordinates presents a preset feature. The second positioning module is used to perform error compensation on the original carrier phase observation based on the corrected feature components, calculate the second position coordinates, and output them. The location management module is used to output the first location coordinates when the residual does not exhibit the preset feature, store the output location coordinates, and output the stored location coordinates when the distance between the current location and the stored location is less than a preset threshold.

[0014] A computer device includes a memory and a processor, the memory storing a computer program, and the processor executing the computer program to implement the steps of a real-time dynamic positioning method for BeiDou satellite signals as described above.

[0015] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of a real-time dynamic positioning method for BeiDou satellite signals as described above.

[0016] Compared with the prior art, the beneficial effects of the present invention are as follows: The observation space is divided into grids and feature attribution is established. A spatial direction tracing mechanism for the source of multi-path features is established, which can trace the feature type and contribution ratio of the residual source through the grid.

[0017] Linear fitting analysis is performed on the residuals and feature attribution degrees to analyze their correlation, trace the source of systematic bias, and determine the spatial direction range that needs to be corrected, thus avoiding the computational overhead of global correction.

[0018] By using observation data from medium Earth orbit satellites in a predetermined space region to correct the multipath characteristic components of geostationary orbit satellites, and then using the corrected multipath characteristic components of geostationary orbit satellites to correct the constant term coefficients of the elevation angle polynomial function of medium Earth orbit satellites, bidirectional collaborative calibration is achieved based on the characteristics of the BeiDou hybrid constellation, eliminating systematic deviations in multipath characteristic components and improving the error compensation effect of carrier phase observations. Attached Figure Description

[0019] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the accompanying drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. In all the drawings, similar elements or parts are generally identified by similar reference numerals. In the drawings, the elements or parts are not necessarily drawn to scale.

[0020] Figure 1 The following is a flowchart of a real-time dynamic positioning method for BeiDou satellite signals based on carrier phase observations, provided as an embodiment of the present invention.

[0021] Figure 2 The flowchart below shows the first position coordinate acquisition process of a real-time dynamic positioning method for BeiDou satellite signals based on carrier phase observations, provided as an embodiment of the present invention.

[0022] Figure 3 The flowchart below shows a second position coordinate acquisition process for a real-time dynamic positioning method for BeiDou satellite signals based on carrier phase observations, provided as an embodiment of the present invention. Detailed Implementation

[0023] 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, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0024] It should be understood that, when used in this specification and the appended claims, the terms "comprising" and "including" indicate the presence of the described features, integrals, steps, operations, elements and / or components, but do not exclude the presence or addition of one or more other features, integrals, steps, operations, elements, components and / or collections thereof.

[0025] It should also be understood that the terminology used in this specification is for the purpose of describing particular embodiments only and is not intended to limit the invention. As used in this specification and the appended claims, the singular forms “a,” “an,” and “the” are intended to include the plural forms unless the context clearly indicates otherwise.

[0026] It should also be further understood that the term "and / or" as used in this specification and the appended claims refers to any combination of one or more of the associated listed items and all possible combinations, and includes such combinations.

[0027] Example 1, referring to Figures 1-3 This is one embodiment of the present invention, which provides a real-time dynamic positioning method for BeiDou satellite signals based on carrier phase observations, including: S1: Identify multipath contaminated observations from the original carrier phase observations, extract the first multipath feature component from the contaminated observations of the first type of satellites, extract the second multipath feature component from the contaminated observations of the second type of satellites, perform error compensation on the original carrier phase observations based on the first and second multipath feature components, and calculate the first position coordinates. Step S1 includes the following steps, S11-S17: S11: Obtain raw carrier phase observations from the BeiDou satellite system receiver, identify satellites affected by multipath interference, and classify them into a first-class satellite set and a second-class satellite set according to their orbital type. The first-class satellites are geostationary orbit satellites, and the second-class satellites are medium Earth orbit satellites.

[0028] Specifically, the multi-frequency carrier phase observations of the current observation epoch are obtained from the receiver as the raw carrier phase observations. The raw carrier phase observations include the carrier phase measurement values ​​of each visible satellite at multiple frequency points and the corresponding satellite identifier. The satellite identifier is a number used to uniquely identify a satellite in the BeiDou satellite system.

[0029] For each visible satellite, carrier phase observations at the first and second frequency points are selected and linearly combined to form an inter-frequency combined observation. The inter-frequency combined observation eliminates error terms shared by the two frequency points, retaining both frequency-dependent and frequency-independent error terms.

[0030] Collect the values ​​of inter-frequency combined observations over a continuous observation period. Calculate the difference between inter-frequency combined observations at adjacent epochs to obtain the time-domain rate of change sequence.

[0031] It should be noted that inter-frequency combined observations include frequency-dependent error terms and frequency-independent error terms. The frequency-dependent error terms conform to the dispersion relation of electromagnetic waves propagating in the ionosphere, and their time-domain rate of change is inversely proportional to the square of the frequency. The frequency-independent error terms are determined by the geometry of the reflecting environment around the receiver, and their time-domain rate of change is unrelated to the frequency.

[0032] The dispersion relation test sequence is obtained by calculating the ratio of the time-domain rate of change sequence to the square of the frequency. When the dispersion relation test sequence remains constant over consecutive epochs, the inter-frequency combined observations are determined to conform to the dispersion relation. When the dispersion relation test sequence shows a jump over consecutive epochs, the inter-frequency combined observations are determined to not conform to the dispersion relation, and the satellite observations are contaminated by error terms independent of frequency.

[0033] Satellite identifiers of satellites contaminated by frequency-independent error terms are included in the set of satellite identifiers contaminated by multipath.

[0034] Satellite orbit types are identified based on satellite identifiers. Satellite identifiers with geostationary orbit types are classified into the first category of satellite identifiers, and satellite identifiers with medium Earth orbit types are classified into the second category of satellite identifiers.

[0035] S12: For each satellite in the first type of satellite set, extract the carrier phase observation value of this satellite from the original carrier phase observation value, and extract the offset part of this carrier phase observation value that remains constant during the observation period as the first multipath feature component of this satellite.

[0036] Specifically, the spatial position of the first type of satellite relative to the receiver remains quasi-stationary during the observation period. The quasi-stationary state is defined as a state in which the change in spatial position during the observation period is less than a preset first spatial threshold. In the quasi-stationary state, the geometric path difference between the direct path and the reflected path of signal propagation remains unchanged during the observation period, resulting in a constant offset characteristic of multipath error in the carrier phase observation values.

[0037] For each satellite in the first type of satellite set, the carrier phase observation values ​​of that satellite during the continuous observation period are extracted from the original carrier phase observation values ​​according to its satellite identifier, forming the first original sequence.

[0038] The first original sequence is subjected to first-order time-difference processing to obtain the first time-difference sequence. First-order time-difference processing extracts the difference between observations at adjacent epochs in the first original sequence. First-order time-difference processing eliminates components in the first original sequence that remain constant along the time axis, while retaining components that change along the time axis.

[0039] The first time difference sequence is subjected to second-order time differencing to obtain a first second-order difference sequence. Second-order time differencing extracts the differences between adjacent epochs in the first time difference sequence. It eliminates components that exhibit linear changes along the time axis and retains components that exhibit non-linear changes along the time axis.

[0040] It should be noted that first-order time difference processing eliminates constant components, second-order time difference processing eliminates linearly varying components, and time-domain accumulation processing weakens random components and extracts systematic components.

[0041] The first second-order difference sequence is accumulated in the time domain to obtain the first accumulated sequence. The time-domain accumulation process involves epoch-by-epoch summation of the first second-order difference sequence along the time axis. This process reduces the random components of observation noise and extracts the systematic components of multipath error.

[0042] Values ​​representing the systematic components are extracted from the first cumulative sequence. These values ​​are obtained through statistical extraction, which involves calculating the center position of the first cumulative sequence within the observation period. These values ​​representing the systematic components are then used as the first multipath characteristic component of this satellite.

[0043] S13: For each satellite in the second type of satellite set, extract the carrier phase observation value of this satellite from the original carrier phase observation value, obtain the satellite elevation angle of this satellite, establish the functional relationship between the carrier phase observation value and the satellite elevation angle, and extract the parameters of this functional relationship as the second multipath feature component of this satellite.

[0044] Specifically, the second type of satellite remains in motion relative to the receiver during the observation period. This motion is defined as a change in spatial position exceeding a preset second spatial threshold during the observation period. In this motion state, the geometric path difference between the direct and reflected paths of signal propagation changes with the satellite's spatial position, causing multipath errors to exhibit spatial position-dependent characteristics in the carrier phase observations.

[0045] For each satellite in the second type of satellite set, the carrier phase observation values ​​for that satellite during the continuous observation period are extracted from the original carrier phase observation values ​​based on its satellite identifier, forming a second original sequence. Based on the satellite's broadcast ephemeris and receiver position, the satellite's elevation angle for the continuous observation period is calculated, forming an elevation angle sequence.

[0046] The second original sequence is subjected to first-order time difference processing to obtain the second time difference sequence.

[0047] The elevation angle sequence is divided into several elevation angle intervals according to their numerical range. For each elevation angle interval, the corresponding values ​​in the second time difference sequence are extracted. These values ​​are then statistically extracted to obtain the representative value of the interval, which is used to calculate the center position value of the numerical set.

[0048] It should be noted that changes in satellite geometric position and atmospheric delay exhibit spatial continuity within adjacent elevation angle intervals, defined as the difference between representative values ​​in adjacent intervals being less than a preset continuity threshold. Changes in multipath error exhibit spatial discontinuity, defined as the difference between representative values ​​in adjacent intervals being greater than a preset continuity threshold.

[0049] Spline interpolation is performed on the representative values ​​of each elevation angle interval to establish a first spatial reference function between the interval representative values ​​and the elevation angle. Spline interpolation establishes a continuous transition relationship between adjacent elevation angle intervals. The first spatial reference function characterizes the spatially continuous components of satellite geometric position variation and atmospheric delay variation.

[0050] The function values ​​of the first spatial reference function at each epoch are subtracted from the corresponding values ​​of the second time difference sequence to obtain the first residual sequence. The first residual sequence represents the multipath error variation components with spatial discontinuity.

[0051] A second spatial reference function is established between the first residual sequence and the elevation angle sequence. The second spatial reference function adopts the form of an elevation angle polynomial function.

[0052] The elevation angle polynomial function is expressed as: in, This represents the predicted multipath error value. For the satellite elevation angle, The coefficient of the constant term, The coefficient of the linear term, The coefficient of the quadratic term.

[0053] It should be noted that the elevation angle polynomial function is obtained by performing a least-squares fit between the first residual sequence and the elevation angle sequence. The least-squares fit involves adjusting three coefficient parameters to minimize the sum of squared deviations between the elevation angle polynomial function and the first residual sequence. (Constant term coefficients) This represents the numerical value of the multipath error at zero elevation angle, which provides a constant contribution to the multipath error across the entire elevation angle range. Linear term coefficient. Characterizing the linear rate of change of multipath error with elevation angle, when When the value is positive, the multipath error increases linearly with the elevation angle. When the value is negative, the multipath error decreases linearly with increasing elevation angle. (Quadratic coefficient) Characterizing the nonlinear rate of change of multipath error with elevation angle, when When the value is positive, the multipath error exhibits an upward convex curve characteristic as it changes with the elevation angle. When the value is negative, the multipath error exhibits a concave curve characteristic as it changes with the elevation angle. The technical reason for using a quadratic polynomial form is that: During the observation period, Type II satellites move from low to high elevation angles. The multipath error is affected by the geometry of the reflection environment, resulting in differences in the reflection path length at different elevation angle positions. This difference presents a nonlinear distribution pattern with a single peak or valley on the elevation angle axis, and the quadratic term can fully characterize this nonlinearity. The technical advantages of the elevation angle polynomial function compared to the piecewise averaging method are: the function is continuously differentiable across the entire elevation angle range, allowing for the calculation of predicted multipath error values ​​at any elevation angle position; the function's extrapolation capability enables prediction of multipath error even outside the elevation angle range covered during the observation period; and the compact representation of the three coefficient parameters reduces storage and computational overhead.

[0054] Extracting the three coefficient parameters of the elevation angle polynomial function and This serves as the second multipath feature component of the satellite.

[0055] S14: Based on the spatial positions of each satellite in the first and second type of satellite sets, the observation space is divided into grids according to the satellite azimuth and satellite elevation angles. For each grid, the affiliation relationship between this grid and the first and second type of satellites is established.

[0056] Specifically, based on the satellite broadcast ephemeris and receiver position, the azimuth and elevation angles of each satellite in the first type of satellite set are calculated, and the azimuth and elevation angles of each satellite in the second type of satellite set at each epoch are calculated.

[0057] The numerical range of satellite azimuth angle is divided into several azimuth angle intervals, and the numerical range of satellite elevation angle is divided into several elevation angle intervals. Each combination of azimuth angle interval and elevation angle interval forms a grid.

[0058] It should be noted that the purpose of grid partitioning is to establish a spatial direction tracing mechanism for multipath errors. The reflection environment around the receiver has different reflection characteristics in different spatial directions, leading to differences in multipath error characteristics across these directions. Grid partitioning discretizes the observation space into a finite number of directional units, establishing a correspondence between spatial directions and the sources of multipath features. This grid mechanism provides spatial tracing capability for subsequent positioning residual analysis. When a positioning residual occurs in a certain spatial direction, the grid can be used to trace back to the type and contribution ratio of the feature source in that direction, enabling directional correction.

[0059] For each grid cell, the coverage of this grid cell by the first type of satellite set and the second type of satellite set is statistically analyzed.

[0060] For each satellite in the first type of satellite set, determine whether its azimuth angle falls within the azimuth angle range of the grid and whether its elevation angle falls within the elevation angle range of the grid. When the conditions are met, record the satellite identifier of that satellite in the first coverage identifier set of this grid. Count the number of satellite identifiers in the first coverage identifier set, and use this count as the first coverage count of this grid.

[0061] For each satellite in the second type of satellite set, at each observation epoch, determine whether its azimuth angle falls within the azimuth angle interval of the grid and whether its elevation angle falls within the elevation angle interval of the grid. When the conditions are met, record the satellite identifier and corresponding epoch to the second coverage identifier set of this grid. Count the number of epochs in the second coverage identifier set as the second coverage count of this grid.

[0062] Calculate the sum of the first coverage quantity and the second coverage quantity to obtain the total coverage quantity of this raster. Calculate the ratio of the first coverage quantity to the total coverage quantity to obtain the first type feature classification degree of this raster. Calculate the ratio of the second coverage quantity to the total coverage quantity to obtain the second type feature classification degree of this raster.

[0063] It should be noted that feature attribution quantifies the composition of multipath feature sources for the spatial direction corresponding to a grid. The first type of feature attribution reflects the proportion of multipath features contributed by quasi-stationary satellites in that spatial direction, while the second type reflects the proportion contributed by moving satellites. The technical role of feature attribution is reflected in subsequent positioning residual analysis: when a residual appears in the spatial direction corresponding to a grid, by comparing the numerical relationship between the first and second type feature attributions for that grid, the correlation between the residual in that grid direction and the multipath features of the first and second type satellites can be determined, providing directional guidance for the subsequent orientation correction of feature components.

[0064] The first type of feature attribution degree and the second type of feature attribution degree are used as the feature attribution weight vector for this grid.

[0065] S15: For each grid cell, extract the corresponding multipath feature components based on the satellites covering the grid cell, and calculate the multipath error prediction value of the grid cell based on the feature attribution weight vector of the grid cell.

[0066] Specifically, based on the first set of coverage identifiers of the raster, the first multipath feature components of the corresponding satellites are extracted to form a first feature component set. Statistical extraction is performed on the first feature component set to obtain a first type of feature reference value. The statistical extraction is used to calculate the center position value of the first feature component set.

[0067] Based on the second coverage identifier set of the raster, the second multipath feature components of the corresponding satellites are extracted to form a second feature component set. The center elevation angle of the raster is calculated, which is the midpoint value of the elevation angle interval of the raster. For each second multipath feature component in the second feature component set, the center elevation angle is substituted into its corresponding elevation angle polynomial function to calculate the function value of the elevation angle polynomial function at the center elevation angle, forming a function value set. Statistical extraction is performed on the function value set to obtain the second type of feature reference value.

[0068] Using the feature assignment weight vector of the raster, a weighted combination calculation is performed on the first type of feature reference value and the second type of feature reference value. The weighted combination calculation is the product of the first type of feature reference value and the first type of feature assignment degree plus the product of the second type of feature reference value and the second type of feature assignment degree. The result of the weighted combination calculation is used as the multipath error prediction value for this raster.

[0069] S16: For each visible BeiDou satellite in the current observation epoch, find the corresponding grid according to the spatial position of the satellite, extract the multipath error prediction value of the grid to compensate the carrier phase observation value of the satellite, and combine the compensated carrier phase observation values ​​of all satellites to form the first compensated observation value.

[0070] Specifically, calculate the satellite's azimuth and elevation angles. Determine the azimuth interval to which the satellite belongs based on the azimuth angle, and the elevation interval to which it belongs based on the elevation angle. Determine the grid cell to which the satellite belongs based on the azimuth and elevation intervals. Extract the multipath error prediction value for this grid cell.

[0071] The compensated carrier phase observations of this satellite are obtained by subtracting the multipath error prediction from the original carrier phase observations. These compensated carrier phase observations from all visible BeiDou satellites at the current observation epoch are then combined to form the first compensated observation.

[0072] S17: Based on the first compensated observation, perform carrier phase differential positioning calculation, fix the integer ambiguity and calculate the receiver position to obtain the first integer ambiguity and the first position coordinates.

[0073] Specifically, among all visible BeiDou satellites in the current observation epoch, the satellite ranked first according to its elevation angle is selected as the reference satellite. The difference between the compensated carrier phase observations of other satellites and the reference satellite is calculated to form a single-difference observation.

[0074] A reference station at a known location is selected as the reference station, and its carrier phase observations are obtained. The difference between the single-difference observations of the current receiver and those of the reference station is calculated to form double-difference observations. A carrier phase differential observation equation based on the double-difference observations is then established.

[0075] A Kalman filter algorithm is used to iteratively solve the carrier phase differential observation equation, obtaining a floating-point solution for integer ambiguity. Integer search processing is then performed on the floating-point solution to obtain a set of integer candidate solutions. The reliability of the integer candidate solution set is verified by calculating the ratio of the optimal integer solution to the second-best integer solution. If the ratio is greater than a preset reliability threshold, the optimal integer solution is considered to meet the reliability condition. The optimal integer solution that meets the reliability condition is fixed as the first integer ambiguity.

[0076] The carrier phase differential observation equation is solved based on the first integer ambiguity to obtain the receiver position, which is used as the first position coordinate.

[0077] S2: Calculate the residual based on the first position coordinates. When the residual exhibits a preset characteristic, extract the observation data of the second type of satellite in the preset space region from the contaminated observation values ​​of the second type of satellite. Correct the first multipath feature component based on the observation data. Correct the second multipath feature component based on the corrected first multipath feature component. Perform error compensation on the original carrier phase observation value based on the corrected feature component and solve to obtain the second position coordinates and output them. Step S2 includes the following steps, S21-S25: S21: Calculate the residual based on the first position coordinates.

[0078] Specifically, based on the first position coordinates and the receiver's three-dimensional coordinates, combined with the satellite position information in the satellite broadcast ephemeris, the geometric distances between each visible satellite and the receiver at the current observation epoch are calculated. Based on the geometric distances, satellite clock corrections, receiver clock estimations, tropospheric delay models, and ionospheric delay models, the theoretical carrier phase values ​​for each satellite are calculated.

[0079] The carrier phase residuals of each satellite are obtained by subtracting the corresponding theoretical carrier phase values ​​from the original carrier phase observations. These residuals are then used as the final residuals.

[0080] S22: When the residuals exhibit preset characteristics, extract the observation data of the second type of satellite in the preset space region from the contaminated observations of the second type of satellite.

[0081] Specifically, the residuals of each satellite are mapped to the corresponding grid. The grid to which each satellite belongs is determined based on its azimuth and elevation angles, and the residuals are assigned to the corresponding grids to form a grid residual distribution.

[0082] Analyze the spatial correlation characteristics of the raster residual distribution. Statistically count the residual values ​​of adjacent raster cells and calculate the difference between their residuals. When the difference between the residuals of adjacent raster cells is within a preset spatial correlation range, the raster residual distribution is determined to exhibit spatial correlation, and spatial correlation is a preset feature.

[0083] It should be noted that the spatial correlation characterization residuals exhibit a systematic distribution pattern in space, indicating that there is a systematic bias in the extraction of multi-path feature components. By analyzing the correlation between the residuals and the feature attribution degree, the source of the systematic bias can be traced and targeted correction can be achieved.

[0084] When the grid residual distribution exhibits spatial correlation, it is determined that the residual exhibits a preset characteristic.

[0085] Based on the feature assignment weight vector of each raster, the correlation between raster residuals and the first type of feature assignment degree is analyzed. The first type of feature assignment degree values ​​and raster residual values ​​of each raster are extracted to establish an assignment degree-residual data pair set. A linear fit is performed on the assignment degree-residual data pair set to obtain the first correlation line. The linear fit involves determining the line parameters using the least squares method to minimize the sum of squared deviations between the line and the data pair set. The slope of the first correlation line is calculated. When the slope of the first correlation line is greater than a preset first slope threshold, a positive correlation is determined between the raster residuals and the first type of feature assignment degree.

[0086] Extract the second-type feature attribution values ​​and raster residual values ​​for each raster, and establish an attribution-residual data pair set. Perform linear fitting on the attribution-residual data pair set to obtain a second correlation line. Calculate the slope of the second correlation line. When the slope of the second correlation line is greater than a preset second slope threshold, it is determined that there is a positive correlation between the raster residual and the second-type feature attribution.

[0087] It should be noted that the positive correlation between raster residuals and feature attribution degrees indicates that the systematic components of the residuals are caused by the extraction bias of the corresponding multipath feature components. By tracing the source of the residuals through the attribution degree, the type and spatial direction range of the feature components that need to be corrected can be determined.

[0088] When there is a positive correlation between the raster residual and the first type of feature attribution degree, rasters with first type of feature attribution degrees greater than a preset first attribution threshold are selected from the raster set to form a first target raster set. The spatial direction range corresponding to each raster in the first target raster set is extracted to form a first spatial direction set.

[0089] When there is a positive correlation between the raster residual and the second-type feature attribution degree, rasters with second-type feature attribution degrees greater than a preset second attribution threshold are selected from the raster set to form a second target raster set. The spatial direction range corresponding to each raster in the second target raster set is extracted to form a second spatial direction set.

[0090] A preset spatial region is determined based on the first set of spatial directions and the second set of spatial directions. The preset spatial region is the intersection of the first set of spatial directions and the second set of spatial directions, and the intersection region is the range of spatial directions that simultaneously belong to the first set of spatial directions and the second set of spatial directions.

[0091] It should be noted that the preset space region is covered by both Type I and Type II satellites. Type I satellites maintain a quasi-stationary state in this region, providing a space reference, while the trajectories of Type II satellites pass through this region, providing dynamic sampling. By extracting measured multipath features from the observation data of Type II satellites in this region and comparing them with the multipath features of Type I satellites in the region, the accuracy of the first multipath feature component can be verified and corrected.

[0092] For each satellite in the second type of satellite set, based on its satellite identifier, extract the carrier phase observation values ​​and corresponding azimuth and elevation angles for that satellite during the continuous observation period from the original carrier phase observation values. Determine whether the satellite azimuth and elevation angles at each epoch fall within a preset spatial region. When the condition is met, extract the carrier phase observation values ​​for that epoch to form a regional observation set. Use this regional observation set as the observation data for the second type of satellites within the preset spatial region.

[0093] S23: Correct the first multipath feature component based on the observation data.

[0094] Specifically, for each satellite in the second type of satellite set, observation values ​​are extracted from its regional observation set. First-order time difference processing is then performed on the regional observation set to obtain a regional time difference sequence.

[0095] Statistical extraction is performed on the regional time difference sequence to obtain regional multipath feature values. These extracted values ​​are used to calculate the center position of the regional time difference sequence. The regional multipath feature values ​​characterize the measured multipath error characteristics of the second type of satellite in a predetermined space region.

[0096] For satellites covering a predetermined space region within the first type of satellite set, their first multipath feature components are extracted to form a first type of regional feature set. Statistical extraction is then performed on the first type of regional feature set to obtain reference values ​​for the first type of region.

[0097] The difference between the regional multipath characteristic value of the second type of satellite and the regional reference value of the first type of satellite is calculated to form the characteristic deviation value.

[0098] It should be noted that both Type I and Type II satellites observe the same reflection environment within a preset space region, and theoretically, their multipath errors should be consistent. The existence of the feature deviation value indicates a systematic bias in the extraction of the first multipath feature component. This bias arises because the spatial position of the Type I satellite is fixed, and its multipath feature extraction relies on long-term statistical analysis. During the statistical process, temporal changes in the reflection environment or spatial drift of the receiver position accumulate as systematic bias. Temporal changes include changes in the physical state of the reflector surface and changes in the geometric structure of the surrounding environment, while spatial drift includes loosening of the receiver hardware and offset of the measurement reference point. By superimposing the feature deviation value onto the first multipath feature component, the multipath features of the Type I satellite in the preset space region can be calibrated to the measured values ​​of the Type II satellite, eliminating the systematic bias. Since the systematic bias remains constant during the observation period, after calibration for the Type I satellites covering the preset space region, the systematic bias in the multipath features of all spatial directions covered by these satellites is simultaneously eliminated.

[0099] For each satellite in the first type of satellite set covering a predetermined space region, the feature deviation value is superimposed on its first multipath feature component to obtain the corrected first multipath feature component. The superposition operation involves adding the value of the first multipath feature component to the value of the feature deviation value.

[0100] S24: Correct the second multipath feature component based on the corrected first multipath feature component.

[0101] Specifically, for satellites covering a preset space region in the first type of satellite set, the corrected first multipath feature component is extracted to form a corrected first type of region feature set. Statistical extraction is performed on the corrected first type of region feature set to obtain the corrected first type of region reference value.

[0102] For each satellite in the second type of satellite set, its second multipath feature component is extracted. The second multipath feature component is the coefficient parameter of the elevation angle polynomial function, the form of which has been defined in step S13.

[0103] Calculate the center elevation angle of the preset space region. Substitute the center elevation angle into the elevation angle polynomial function, calculate the function value, and use it as the theoretical multipath characteristic value of the second type of satellite in the preset space region.

[0104] The difference between the theoretical multipath characteristic value of the second type of satellite in the preset space region and the corrected reference value of the first type of region is calculated as the adjustment amount of the second type of characteristic.

[0105] It should be noted that the spatial position of the first type of satellite in the preset space region remains unchanged during the observation period, and its multipath error is a constant value, characterized by the corrected first multipath feature component. The spatial position of the second type of satellite in the preset space region moves from a low elevation angle to a high elevation angle over time, and its multipath error varies with the elevation angle, characterized by an elevation angle polynomial function. When the second type of satellite moves to the center elevation angle position of the preset space region, the first and second type of satellites are in the same spatial direction, observing the same reflection environment, and their multipath errors should theoretically be equal. The adjustment amount of the second feature reflects the difference in their multipath errors. The technical reason for this difference is that the extraction of the second multipath feature component is based on fitting observation data across the entire elevation angle range. The fitting process pursues global optimization rather than local optimization, leading to the accumulation of fitting errors at local positions within the preset space region. This difference is addressed by adjusting the constant term coefficients of the elevation angle polynomial function. This allows the function value at the center elevation angle of a preset spatial region to be equal to the corrected first multipath characteristic component, achieving local calibration. The coefficients of the first-order terms are preserved. and quadratic coefficient The unchanging technical reason is: and The linear and nonlinear rates of change of multipath error with elevation angle are characterized. These rates are determined by the trajectory of the second-type satellite and the spatial geometry of the reflection environment, and are not affected by single-point errors in the preset space region. The technical innovation of adjusting only the constant term coefficients lies in: utilizing the spatial reference characteristics of the first-type satellite in the preset space region to calibrate the multipath error benchmark value of the second-type satellite in that region, while maintaining the rate of change characteristics of the second-type satellite across the entire elevation angle range, thus achieving collaborative calibration based on the characteristics of the BeiDou hybrid constellation.

[0106] For each satellite in the second type of satellite set, extract the constant term coefficients of its second multipath feature component. coefficient of the first term and quadratic coefficient .

[0107] The coefficients of the corrected constant term are calculated according to the following rules: in, The constant term coefficients of the second multipath characteristic component are... This is the adjustment amount for the second type of feature. These are the coefficients of the corrected constant term.

[0108] It should be noted that the calculation rule for the corrected constant term coefficients is achieved by directly superimposing the second type of feature adjustment amount onto the original constant term coefficients. The technical implementation details of this rule are as follows: when the value of the second type of feature adjustment amount is positive, it indicates that the theoretical multipath characteristic value of the second type of satellite in the preset space region is less than the corrected reference value for the first type of region, and the constant term coefficients are adjusted upwards; when the value of the second type of feature adjustment amount is negative, it indicates that the theoretical multipath characteristic value of the second type of satellite in the preset space region is greater than the corrected reference value for the first type of region, and the constant term coefficients are adjusted downwards. This superposition operation causes the function values ​​of the elevation angle polynomial function to increase or decrease synchronously at all elevation angle positions. This rule achieves a translation of the function curve along the vertical axis. The technical effects are as follows: the multipath error prediction values ​​of the second-type satellites in the preset space region are consistent with the corrected measured values ​​of the first-type satellites, eliminating prediction bias in the preset space region; the multipath error prediction values ​​of the second-type satellites in non-preset space regions maintain their original linear and nonlinear rates of change, ensuring that the prediction accuracy in other space directions is not affected by local calibration; the multipath characteristics of the first and second-type satellites are consistent in the preset space region, and systematic biases in all space directions are eliminated synchronously. This rule is based on the unique characteristics of the BeiDou hybrid constellation: the first-type satellites provide space reference anchor points, the second-type satellites provide space rates of change, and the space reference is unified through constant term adjustment, while motion characteristics are preserved through rate of change maintenance.

[0109] Keep the coefficient of the first term and quadratic coefficient The value remains unchanged. The coefficient of the corrected constant term will be... coefficient of the first term and quadratic coefficient As the corrected second multipath feature component.

[0110] S25: Based on the corrected characteristic components, perform error compensation on the original carrier phase observation values, calculate the second position coordinates, and output them.

[0111] Specifically, for each grid cell, the corrected first multipath feature component of the corresponding satellite is extracted based on the grid cell's first coverage identifier set. The corrected first multipath feature component is then statistically extracted to obtain the corrected first-class feature reference value.

[0112] Based on the second coverage identifier set of the raster, the corrected second multipath feature components of the corresponding satellites are extracted. For each corrected second multipath feature component, the center elevation angle of the raster is substituted into the corrected elevation angle polynomial function to calculate the function value. Statistical extraction is performed on the function values ​​to obtain the corrected second-type feature reference value.

[0113] Using the feature assignment weight vector of the raster, a weighted combination is performed on the corrected first-class feature reference value and the corrected second-class feature reference value to obtain the corrected multipath error prediction value of the raster. The weighted combination calculation is the product of the corrected first-class feature reference value and the first-class feature assignment degree plus the product of the corrected second-class feature reference value and the second-class feature assignment degree.

[0114] For each visible BeiDou satellite at the current observation epoch, the satellite's grid is determined based on its azimuth and elevation angles, and the corrected multipath error prediction value for that grid is extracted. The corrected multipath error prediction value is subtracted from the original carrier phase observation values ​​of the satellite to obtain the second compensated carrier phase observation value. The second compensated carrier phase observation values ​​of all satellites are combined to form the second compensated observation value.

[0115] Among all visible BeiDou satellites, the satellite ranked first according to its elevation angle is selected as the reference satellite. The difference between the second compensated carrier phase observations of other satellites and the reference satellite is calculated to form the second single-difference observation. A reference station with a known location is selected as the reference station, and the difference between the second single-difference observation of the current receiver and the reference station is calculated to form the second double-difference observation.

[0116] Establish the carrier phase double-difference observation equation.

[0117] The carrier phase double-difference observation equation is expressed as: in, For satellite and satellite The epoch The second double-difference observation, For satellite and satellite The epoch The double difference geometric distance, For carrier wavelength, For satellite and satellite Double-difference integer ambiguity, For satellite and satellite The epoch Double-difference observation noise.

[0118] It should be noted that the carrier phase double-difference observation equation establishes the mathematical relationship between the second double-difference observation value and the unknown parameters. This equation is the fundamental observation model for carrier phase differential positioning. Double-difference geometric distance. The calculation depends on the current receiver's three-dimensional coordinates, the base station's three-dimensional coordinates, and the satellite's... 3D coordinates and satellites The three-dimensional coordinates of the receiver are given, where the three-dimensional coordinates of the current receiver are unknown parameters. Double-difference geometric distance is calculated by determining the distance between the current receiver and the satellite. Euclidean distance, base station and satellite Euclidean distance, current receiver and satellite Euclidean distance, base station and satellite The Euclidean distance is then used to obtain the difference combination. Double-difference integer ambiguity. This is an integer unknown parameter characterizing the integer uncertainty in the carrier phase observations; this parameter remains constant throughout the observation period. Double-difference observation noise. This includes residual multipath error, receiver thermal noise, residual atmospheric delay, and model error. The second double-difference observation is calculated based on the second compensated carrier phase observation, which subtracts the corrected multipath error prediction. Compared to the uncorrected multipath error prediction, the corrected multipath error prediction eliminates the systematic biases of the first and second multipath characteristic components. The elimination of systematic biases reduces the observation noise of the second double-difference observation. The statistical characteristics of the first double-difference observation are superior to those of the first double-difference observation, as the mean of the observation noise is closer to zero and the variance of the observation noise is reduced. This improvement in the statistical characteristics of the observation noise directly affects the solution accuracy of the carrier phase double-difference observation equation, thus reducing the double-difference integer ambiguity. The floating-point solution estimation accuracy is improved, the success rate of integer fixation is improved, and the estimation accuracy of the current receiver's three-dimensional coordinates is improved. The technical advantages of using the double-difference form of the carrier phase double-difference observation equation are as follows: Double-difference processing eliminates the common-mode error of satellite clock error and receiver clock error, eliminating the need to estimate clock error parameters; Double-difference processing reduces the common-mode error of tropospheric delay and ionospheric delay. The degree of common-mode error reduction depends on the baseline length between the reference station and the current receiver. When the baseline length is less than a preset baseline threshold, the common-mode error reduction effect meets the positioning accuracy requirements; Double-difference processing retains the integer characteristics of integer ambiguity, and the estimation accuracy of the current receiver's three-dimensional coordinates is improved after the integer is fixed.

[0119] The first integer ambiguity is taken as the double-difference integer ambiguity. The initial estimate is obtained by using the first position coordinate as the initial estimate of the current receiver's three-dimensional coordinates. An extended Kalman filter algorithm is used to iteratively solve the carrier phase double-difference observation equation.

[0120] It should be noted that the Extended Kalman Filter (EKF) algorithm establishes state equations and observation equations, and achieves epoch-wise optimization of parameter estimates through recursive calculations of state prediction and observation updates. The state equations characterize the temporal evolution of the state vector in adjacent epochs, and the state vector contains the current receiver's 3D coordinates and double-difference integer ambiguities. For carrier-phase differential positioning applications, the change of the current receiver's 3D coordinates in adjacent epochs is determined by the receiver's motion velocity, while the double-difference integer ambiguities remain unchanged in adjacent epochs. The state equations describe the temporal propagation of the state vector through a state transition matrix, the form of which is determined based on the receiver's motion model. The observation equations are linearized forms of the carrier-phase double-difference observation equations with respect to the state vectors. The observation matrix is ​​obtained by taking the partial derivatives of the carrier-phase double-difference observation equations with respect to the current receiver's 3D coordinates and double-difference integer ambiguities. The state prediction step of the EKF algorithm calculates the predicted state value for the current epoch based on the state equations and the state estimate from the previous epoch. The observation update step calculates the Kalman gain based on the observation equations, the second double-difference observation value for the current epoch, and the predicted state value. The state estimate is obtained by correcting the predicted state value using the Kalman gain. The calculation of Kalman gain depends on the covariance matrix of state prediction and the covariance matrix of observation noise. The observation noise covariance matrix reflects the accuracy of the second double-difference observation. The observation noise variance of the second double-difference observation is smaller than that of the first double-difference observation, and the diagonal elements of the observation noise covariance matrix are smaller. This results in a larger weight for observation updates in the Kalman gain and a greater correction magnitude from the state estimate to the observation. The technical advantages of using the first integer ambiguity and the first position coordinates as the initial estimate are: the initial estimate is close to the true value, the initial error of state prediction is small, and the convergence speed of the extended Kalman filter algorithm is faster; the initial estimate is obtained based on the first compensated observation, which already includes some multipath error compensation effects, avoiding convergence instability caused by solving from scratch. The technical advantages of the extended Kalman filter algorithm compared to the batch least squares algorithm are: recursive epoch-by-epoch update, low computational complexity, suitable for real-time dynamic positioning applications; fusion of multi-epoch observation information, resulting in high stability of the state estimate; and real-time output of the current epoch's position estimate, meeting real-time positioning requirements.

[0121] After multiple epochs of state prediction and observation update iterations, the double-difference integer ambiguity is obtained. The floating-point solution and the estimated three-dimensional coordinates of the current receiver.

[0122] An integer search is performed on the floating-point solutions of the double-difference integer ambiguity to obtain a set of integer candidate solutions. The integer search searches the integer space near the floating-point solutions to find the integer combination that minimizes the sum of squared residuals of the carrier phase double-difference observation equation. A ratio test is performed on the set of integer candidate solutions to calculate the ratio between the optimal and second-best integer solutions. The criterion for the ratio test is: if the ratio is greater than a preset reliability threshold, the optimal integer solution is deemed to meet the reliability condition, and the optimal integer solution is fixed as the second integer ambiguity.

[0123] Substituting the fixed value of the second integer cycle ambiguity into the carrier phase double-difference observation equation, the final solution of the current receiver's three-dimensional coordinates is obtained, yielding the second position coordinates. The second position coordinates are then output.

[0124] S3: When the residual does not exhibit the preset characteristics, output the first position coordinates; S4: Store the output position coordinates. In response to the current position being less than a preset threshold, output the stored position coordinates.

[0125] Specifically, when step S2 outputs the second position coordinates, the second position coordinates are stored in the position coordinate memory. When step S3 outputs the first position coordinates, the first position coordinates are stored in the position coordinate memory.

[0126] In subsequent observation epochs, the spatial distance between the currently output position coordinates and the stored position coordinates is calculated. When the spatial distance is less than a preset distance threshold, the stored position coordinates are directly output without needing to repeat the multipath feature extraction and error compensation in steps S1 and S2.

[0127] Example 2, an embodiment of the present invention, provides a real-time dynamic positioning method for BeiDou satellite signals based on carrier phase observations, based on the previous embodiment, including: In step S11, the determination of constancy and jump in the dispersion relation test sequence is based on the statistical characteristics of the sequence itself. The dispersion relation test sequence is statistically analyzed over a continuous observation period, and the mean and standard deviation of the sequence are calculated. When the change in sequence values ​​between adjacent epochs is less than the standard deviation, it is determined to be constant; when the change in sequence values ​​between adjacent epochs is greater than the sum of the mean and the standard deviation, it is determined to be a jump.

[0128] In steps S12 and S13, the preset first spatial threshold and the preset second spatial threshold are determined based on the statistical distribution of satellite spatial position changes. The spatial position changes of each satellite in the first type of satellite set within adjacent epochs of the observation period are statistically analyzed, and the upper quantile of the change distribution is extracted as the criterion for determining the first spatial threshold. Similarly, the spatial position changes of each satellite in the second type of satellite set within adjacent epochs of the observation period are statistically analyzed, and the lower quantile of the change distribution is extracted as the criterion for determining the second spatial threshold. Satellites in the first type of satellite set whose spatial position changes are less than the first spatial threshold are classified as quasi-stationary; satellites in the second type of satellite set whose spatial position changes are greater than the second spatial threshold are classified as moving.

[0129] In step S13, the preset continuity threshold is determined based on the statistical distribution of the differences between representative values ​​of intervals. The difference between the representative values ​​of all adjacent elevation angle intervals is calculated, and statistical analysis is performed on the differences to calculate the median and interquartile range of the difference distribution. Spatial continuity is determined when the difference between the representative values ​​of adjacent intervals is less than the sum of the median and the interquartile range; spatial discontinuity is determined when the difference between the representative values ​​of adjacent intervals is greater than the sum of the median and the interquartile range.

[0130] In step S14, the number of azimuth and elevation angle intervals is determined based on the spatial uniformity of satellite coverage. The distribution of the first and second type of satellite sets along the azimuth and elevation angle axes is statistically analyzed, and the spatial density of satellite coverage is calculated. The number of azimuth and elevation angle intervals is adjusted to minimize the variance of coverage density for each grid, achieving a spatially uniform distribution of grid coverage. Coverage density is calculated by dividing the total number of grids by the spatial area of ​​the grid.

[0131] In step S22, the preset spatial correlation range is determined based on the spatial variability characteristics of the raster residuals. The residual differences between all adjacent raster pairs are calculated, and a variogram analysis is performed on these residual differences. The spatial distance corresponding to when the variogram reaches the sill value is extracted as the spatial correlation range. The spatial correlation range characterizes the critical distance at which the residuals transition from correlation to independence in space.

[0132] In step S22, the preset first and second attribution thresholds are determined based on the quantiles of the feature attribution degree distribution. Statistical distribution analysis is performed on the first type of feature attribution degree for all rasters, and the upper tertiary of the distribution is extracted as the first attribution threshold. Statistical distribution analysis is also performed on the second type of feature attribution degree for all rasters, and the upper tertiary of the distribution is extracted as the second attribution threshold. The upper tertiary represents the critical value where the value of the attribution degree of this type of feature is within the top quartile of all rasters.

[0133] In step S22, the preset first slope threshold and the preset second slope threshold are determined based on the significance test of linear fitting. Linear regression analysis is performed on the set of membership-residual data, and the standard error of the regression coefficients and the t-statistic are calculated. When the absolute value of the t-statistic is greater than the critical value corresponding to the significance level, the slope is considered significantly different from zero. When the regression coefficient is positive and the significance test passes, a positive correlation is considered between the membership degree and the residuals.

[0134] In steps S17 and S25, the preset reliability threshold is determined based on the statistical test theory of fixed integer ambiguity. The theoretical basis for the ratio test is that the likelihood ratio of the optimal solution and the suboptimal solution follows an F distribution, and the ratio threshold is determined based on the critical value of the F distribution at a given confidence level. The preset baseline threshold is determined based on the spatial decorrelation distance of atmospheric delay. By analyzing the decay law of the correlation coefficient between the reference station and the current receiver with the baseline length, the baseline length corresponding to the correlation coefficient dropping to the preset correlation level is extracted as the baseline threshold.

[0135] In step S4, the preset distance threshold is determined based on the spatial stability range of the multipath characteristics. The changes in the multipath characteristics of the first type of satellites during the observation period are statistically analyzed. When a change in receiver position causes a change in the observation geometry of the first type of satellites, the multipath characteristics change accordingly. By analyzing the relationship between the sky coverage angle interval of the first type of satellites and the amount of change in receiver position, the upper limit of the position change that makes the change in multipath characteristics less than the feature statistical error is extracted as the distance threshold.

[0136] Example 3 is an embodiment of the present invention, which provides a real-time dynamic positioning system for BeiDou satellite signals based on carrier phase observations, comprising: A multipath identification module is used to identify observations contaminated by multipath from the original carrier phase observations; The feature extraction module is used to extract the first multipath feature component and the second multipath feature component from the contaminated observations of the first type of satellite and the second type of satellite, respectively. The first positioning module is used to perform error compensation on the original carrier phase observation value based on the first multipath feature component and the second multipath feature component, and calculate the first position coordinates. The feature correction module is used to correct the first multipath feature component and the second multipath feature component from the observation data of the second type of satellite in the preset space region when the residual calculated based on the first position coordinates presents a preset feature. The second positioning module is used to perform error compensation on the original carrier phase observation based on the corrected feature components, calculate the second position coordinates, and output them. The position management module is used to output the first position coordinates when the residual does not exhibit the preset characteristics, store the output position coordinates, and output the stored position coordinates when the distance between the current position and the stored position is less than a preset threshold.

[0137] This embodiment also provides an electronic device applicable to a real-time dynamic positioning method for BeiDou satellite signals based on carrier phase observations, comprising: a memory and a processor; the memory is used to store computer-executable instructions, and the processor is used to execute the computer-executable instructions to realize the real-time dynamic positioning method for BeiDou satellite signals based on carrier phase observations as proposed in the above embodiment.

[0138] This embodiment also provides a storage medium storing a computer program that, when executed by a processor, implements a real-time dynamic positioning method for BeiDou satellite signals based on carrier phase observations as proposed in the above embodiment.

[0139] The storage medium proposed in this embodiment belongs to the same inventive concept as the method for real-time dynamic positioning of Beidou satellite signals based on carrier phase observations proposed in the above embodiments. Technical details not described in detail in this embodiment can be found in the above embodiments, and this embodiment has the same beneficial effects as the above embodiments.

[0140] Based on the above description of the implementation methods, those skilled in the art can clearly understand that the present invention can be implemented using software and necessary general-purpose hardware, and of course, it can also be implemented using hardware, but in many cases the former is a better implementation method. Based on this understanding, the technical solution of the present invention, or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as a computer floppy disk, read-only memory (ROM), random access memory (RAM), flash memory, hard disk, or optical disk, etc., including several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute the methods of the various embodiments of the present invention.

[0141] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A real-time dynamic positioning method for BeiDou satellite signals, characterized in that, Includes the following steps: The observations contaminated by multipath are identified from the original carrier phase observations. The first multipath feature component is extracted from the contaminated observations of the first type of satellites, and the second multipath feature component is extracted from the contaminated observations of the second type of satellites. The original carrier phase observations are then error-compensated based on the first and second multipath feature components, and the first position coordinates are calculated. The residual is calculated based on the first position coordinates. When the residual exhibits a preset characteristic, the observation data of the second type of satellite in the preset space region is extracted from the contaminated carrier phase observation value of the second type of satellite. The first multipath feature component is corrected based on the observation data. The second multipath feature component is corrected based on the corrected first multipath feature component. The original carrier phase observation value is error compensated based on the corrected feature component and the second position coordinates are obtained and output. When the residual does not exhibit the preset feature, output the first position coordinates; Store the output position coordinates. If the distance between the current position and the stored position is less than a preset threshold, output the stored position coordinates.

2. The real-time dynamic positioning method for BeiDou satellite signals according to claim 1, characterized in that, The process of identifying multipath-contaminated observations from the original carrier phase observations includes: For each visible satellite, carrier phase observations at the first and second frequency points are selected and linearly combined to form inter-frequency combined observations; Collect the values ​​of the frequency combination observations over a continuous observation period, calculate the difference between the frequency combination observations of adjacent epochs, and obtain the time-domain rate of change sequence; Based on the relationship between the time-domain rate of change sequence and frequency, a dispersion relation test sequence is obtained; Determine whether there is a jump in the values ​​of the dispersion relation test sequence in adjacent epochs. If a jump occurs, it is determined that the current satellite observation value is contaminated by multipath.

3. The real-time dynamic positioning method for BeiDou satellite signals according to claim 1, characterized in that, The extraction of the first multipath feature component from the contaminated observations of the first type of satellites includes: The carrier phase observations of the first type of satellites are sequentially processed by first-order time difference processing, second-order time difference processing, and time domain accumulation processing. The center position value is extracted from the accumulation result as the first multipath feature component. The extraction of the second multipath feature component from contaminated observations of the second type of satellites includes: The carrier phase observations of the second type of satellite are processed by first-order time difference to obtain a second time difference sequence. Spline interpolation is performed on the representative values ​​of the elevation angle interval to establish a first spatial reference function. The first spatial reference function is subtracted from the second time difference sequence to obtain a first residual sequence. The coefficient parameters of the elevation angle polynomial function are obtained by least-squares fitting of the first residual sequence as the second multipath feature component.

4. The real-time dynamic positioning method for BeiDou satellite signals according to claim 1, characterized in that, The step of performing error compensation on the original carrier phase observation value based on the first multipath feature component and the second multipath feature component and then calculating the first position coordinate includes: The observation space is divided into grids, and the number of first-type satellites and second-type satellites covered in each grid are counted. The coverage ratio is calculated as the first-type feature attribution degree and the second-type feature attribution degree. Based on the first type of feature attribution degree and the second type of feature attribution degree, the first multipath feature component and the second multipath feature component of each grid-covered satellite are weighted and combined to obtain the multipath error prediction value of the grid. For each visible satellite, the first compensated observation is formed by subtracting the multipath error prediction value of the corresponding grid from the original carrier phase observation value; Based on the first compensated observation, carrier phase differential positioning is calculated to obtain the first integer ambiguity and the first position coordinates.

5. The real-time dynamic positioning method for BeiDou satellite signals according to claim 4, characterized in that, The extraction of observation data of the second type of satellite in the preset space region from the contaminated observation data of the second type of satellite includes: The residual is mapped to a grid to form a grid residual. When the difference between the residuals of adjacent grids is within a preset spatial correlation range, it is determined that the preset feature is presented. A first correlation line is obtained by linearly fitting the grid residuals and the first type of feature attribution degree. When the slope of the first correlation line exceeds a preset first slope threshold, grids with the first type of feature attribution degree exceeding the preset first attribution threshold are selected to form a first spatial direction set. A second correlation line is obtained by linearly fitting the grid residuals and the second type of feature attribution degree. When the slope of the second correlation line exceeds the preset second slope threshold, grids with the second type of feature attribution degree exceeding the preset second attribution threshold are selected to form a second spatial direction set. The intersection region of the first spatial direction set and the second spatial direction set is taken as the preset spatial region, and the carrier phase observation values ​​of the second type of satellite in the preset spatial region are extracted as the observation data.

6. The real-time dynamic positioning method for BeiDou satellite signals according to claim 5, characterized in that, The second multipath feature component is modified based on the modified first multipath feature component, including: The observation data is processed by first-order time difference and the regional multipath feature value is statistically extracted. The difference between the regional multipath feature value and the statistical value of the first multipath feature component of the first type of satellite in the preset space region is calculated as the feature deviation value. The feature deviation value is superimposed on the first multipath feature component to obtain the corrected first multipath feature component. Substitute the center elevation angle of the preset spatial region into the elevation angle polynomial function to obtain the theoretical multipath feature value. Calculate the difference between the theoretical multipath feature value and the statistical value of the corrected first multipath feature component in the preset spatial region as the second type of feature adjustment amount. Add the second type of feature adjustment amount to the constant term coefficient of the elevation angle polynomial function to obtain the corrected second multipath feature component.

7. The real-time dynamic positioning method for BeiDou satellite signals according to claim 6, characterized in that, The step of performing error compensation on the original carrier phase observation based on the corrected feature components, calculating and obtaining the second position coordinates, and outputting the following includes: Based on the corrected first multipath feature component and the corrected second multipath feature component, and combined with the first type of feature attribution degree and the second type of feature attribution degree, the corrected multipath error prediction value of each grid is calculated. For each visible satellite, the second compensated observation is formed by subtracting the corrected multipath error value of the corresponding grid from the original carrier phase observation value. Using the first integer ambiguity and the first position coordinates as initial estimates, the extended Kalman filter algorithm is used to iteratively solve the carrier phase double-difference observation equation established based on the second compensated observation value. The integer ambiguity is fixed and the second position coordinates are calculated and output.

8. A real-time dynamic positioning system for BeiDou satellite signals, characterized in that, The method for real-time dynamic positioning of BeiDou satellite signals as described in any one of claims 1-7 includes: A multipath identification module is used to identify observations contaminated by multipath from the original carrier phase observations; The feature extraction module is used to extract the first multipath feature component and the second multipath feature component from the contaminated observations of the first type of satellite and the second type of satellite, respectively. The first positioning module is used to perform error compensation on the original carrier phase observation value based on the first multipath feature component and the second multipath feature component and calculate the first position coordinates. The feature correction module is used to correct the first multipath feature component and the second multipath feature component from the observation data of the second type of satellite in the preset space region when the residual calculated based on the first position coordinates presents a preset feature. The second positioning module is used to perform error compensation on the original carrier phase observation based on the corrected feature components, calculate and output the second position coordinates; The location management module is used to output the first location coordinates when the residual does not exhibit the preset feature, store the output location coordinates, and output the stored location coordinates when the distance between the current location and the stored location is less than a preset threshold.

9. A computer device, characterized in that, It includes a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps of the BeiDou satellite signal real-time dynamic positioning method according to any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, It stores a computer program, which, when executed by a processor, implements the steps of a real-time dynamic positioning method for BeiDou satellite signals according to any one of claims 1 to 7.

Images