Method for calculating code deviation based on ground-based GNSS (Global Navigation Satellite System) data

Abstract

The invention provides a method for calculating code deviation based on ground-based GNSS (Global Navigation Satellite System) data. The method comprises the following steps of: traversing files of a preprocessing module; performing satellite-by-satellite traversal by a preprocessing module; performing arc-section-by-arc traversal of the preprocessing module; traversing files of the post-processing module; performing satellite-by-satellite traversal of the post-processing module; performing arc-section-by-arc traversal on the post-processing module; and the resolving module is used for resolving the DCB by using a least square method. The beneficial effects of the invention are that through simultaneous processing of multiple satellites and multiple observation station data in a specified range, joint DCB resolving is carried out, more observation information is utilized, resolving precision and reliability are improved, a preprocessing module is added, cycle slip detection and pseudo-range smoothing modules are introduced, result robustness is enhanced, and a basis is provided for regional ionosphere modeling.

Description

Technical Field

[0001] The invention belongs to the field of GNSS applications, and in particular relates to a method for calculating code deviation based on ground-based GNSS data. Background Art

[0002] Code bias is a systematic error caused by hardware delay differences during the transmission (satellite) and reception (receiver) of GNSS signals, which manifests itself as a fixed deviation of pseudorange observations at different frequencies.

[0003] The inversion of Difference of Code Bias (DCB) is crucial to improving the positioning accuracy of the Global Navigation Satellite System (GNSS). By accurately estimating and correcting DCB, the accuracy of the ionospheric total electron content (TEC) model can be improved, thereby improving the effect of ionospheric delay correction, which is crucial to improving the accuracy of GNSS positioning, navigation and timing services. In addition, accurate DCB inversion helps to obtain high-precision ionospheric TEC data, which is of great significance for studying the temporal and spatial variation characteristics of the ionosphere and space weather events (such as ionospheric storms, magnetic storms, etc.). In high-precision GNSS applications (such as precise positioning, precise orbit determination, etc.), the presence of DCB will affect the measurement accuracy. Through DCB inversion and correction, its impact on the measurement results can be effectively reduced and the reliability of GNSS applications can be enhanced. Summary of the invention

[0004] In view of this, the present invention aims to propose a method for calculating code bias based on ground-based GNSS data to solve the problem that the intra-day variation of DCB is generally ignored in the ionospheric observation inversion program, thereby reducing the accuracy of TEC estimation.

[0005] To achieve the above object, the technical solution of the present invention is achieved as follows: A method for calculating code deviation based on ground-based GNSS data comprises the following steps: S1, file traversal of preprocessing module; S2, satellite-by-satellite traversal of the preprocessing module; S3, arc segment by arc segment traversal of the preprocessing module; S4, file traversal of post-processing module; S5, satellite-by-satellite traversal of the post-processing module; S6, arc segment by arc segment traversal of the post-processing module; S7, the solving module solves the DCB using the least square method.

[0006] Furthermore, in step S1, the file traversal of the preprocessing module includes: S11, file reading: read the navigation file, obtain the orbit radius sine term, satellite orbit eccentricity and perigee angular distance parameters, read the observation file, and obtain the specified observation type data; S12, SPP solution: Based on the read navigation files and observation files, perform pseudo-range single-point positioning calculation to obtain the coordinates of the navigation satellite and the altitude angle between the satellite and the receiver; S13, altitude angle & azimuth angle storage: the altitude angle information of different satellites is stored according to the information format, and the stored information includes the altitude angle and its corresponding satellite number and time.

[0007] Furthermore, in step S2, the satellite-by-satellite traversal of the preprocessing module includes: S21, satellite health status judgment: judge whether the SVH value of the current satellite is 0. If not, it is determined that the current satellite has not passed the preprocessing and jumps to the next satellite. Otherwise, perform subsequent processing; S22, observation value storage: based on step S11, store the pseudorange and carrier of the specified time period according to the time and observation value type; S23, cycle slip detection: Geometry-Free is used to detect cycle slips. The formula is as follows: (1) (2) (3) In the formula, is the sequence subscript, , are the wavelengths of bands 1 and 2 respectively, in meters. , They are the phase observations of bands 1 and 2, in cycles. When , it is determined that a cycle slip has occurred, and the observation value at the corresponding time is cleared to zero; S24, non-empty arc statistics: calculate the length of observation values ​​of different observation types that are not 0, and take their intersection. If the current arc length is less than 45, count the next arc, otherwise, perform subsequent processing.

[0008] Furthermore, in step S3, the arc segment-by-arc traversal of the preprocessing module includes: S31, pseudorange smoothing: hatch is used to smooth the pseudoranges of the two specified observation types. The formula is as follows: (4) In the formula, For band The wavelength, in meters, For band Phase observation, in weeks, Indicates that the subscript is The original pseudorange observation value of Indicates that the subscript is The smoothed pseudorange observation value is is the weight value, which depends on the length of the non-empty arc segment; S32, puncture point calculation: Calculate the puncture point based on the altitude angle and azimuth angle. The formula is as follows: (6) (7) (8) (9) like >70° and ,or <-70° and , (10) otherwise, The value of is shown in (11): (11) In the formula, is the altitude angle, is the latitude of the puncture point, is the longitude of the puncture point, is the receiver latitude, is the receiver longitude, is the azimuth, is the angle formed by the receiver, the center of the earth, and the puncture point. is the angle between the center of the earth - the puncture point and the receiver - the satellite ; S33, elevation angle screening: Eliminate data with elevation angles less than 30°. If the current satellite and station still have data after the above step S31, it is determined that the current satellite and station have passed the preprocessing, and then the next station is preprocessed.

[0009] Furthermore, in step S4, the files of the post-processing module are traversed, and the steps are the same as those of step S1.

[0010] Furthermore, in step S5, the satellite-by-satellite traversal of the post-processing module includes: S51, observation value storage, the steps are the same as step S22; S52, cycle slip detection, the steps are the same as step S23; S53, non-empty arc segment statistics, the steps are the same as step S24.

[0011] Furthermore, in step S6, the arc segment-by-arc traversal of the post-processing module includes: S61, pseudorange smoothing, the steps are the same as step S31; S62, construct P4 matrix: each element of P4 is the VTEC value before DCB is eliminated, and the calculation formula is as follows: (5) In the formula, is the altitude angle, is the difference in pseudorange of the specified observation type after smoothing, , Bands , The frequency in Hz, is the radius of the Earth, is the puncture point height; S63, puncture point calculation, the steps are the same as step S32; S64. Construct coefficient matrix B: The coefficient matrix B is composed of spherical harmonic coefficients and DCB coefficient modules. The number of columns of spherical harmonic coefficients is the square of the order. The formula is as follows: (12) (13) In the formula, is the latitude of the puncture point, is the longitude of the puncture point, for Step The normalized Legendre coefficients are concluded in matrix B. If If it is not 0, then formula (13) follows formula (12); The number of columns of the DCB coefficient module is the sum of the number of satellites and the number of stations that have passed the preprocessing. It is composed of the satellite DCB and the receiver DCB. At the corresponding satellite and the corresponding station, the value of the element is as shown in the following formula, and the non-corresponding position is 0; (14); In the formula, is the altitude angle, , Bands , The frequency in Hz, is the speed of light, is the radius of the Earth, is the puncture point height.

[0012] Furthermore, in step S7, the solving module solves DCB using the least square method, including: After the matrix B and matrix P4 are constructed, the least squares method is used to solve the unknown DCB and spherical harmonic coefficients. The formula is as follows: (15); In the formula, To solve the result, we need to calculate the spherical harmonic coefficients, the satellite DCB and the receiver DCB. is the coefficient matrix, is the observation matrix.

[0013] Compared with the prior art, the method for calculating code deviation based on ground-based GNSS data described in the present invention has the following advantages: The present invention processes data from multiple satellites and multiple stations within a specified range simultaneously to perform joint DCB solution, utilizes more observation information, improves solution accuracy and reliability, adds a preprocessing module, introduces modules such as cycle slip detection and pseudorange smoothing, enhances the robustness of the results, and provides a basis for regional ionospheric modeling. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] The accompanying drawings constituting a part of the present invention are used to provide a further understanding of the present invention. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute an improper limitation of the present invention. In the accompanying drawings: Figure 1 This is a technical roadmap described in an embodiment of the present invention. DETAILED DESCRIPTION

[0015] It should be noted that, in the absence of conflict, the embodiments of the present invention and the features in the embodiments may be combined with each other.

[0016] In the description of the present invention, it should be understood that the terms "center", "longitudinal", "lateral", "up", "down", "front", "back", "left", "right", "vertical", "horizontal", "top", "bottom", "inside", "outside" and the like indicate positions or positional relationships based on the positions or positional relationships shown in the accompanying drawings, and are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore cannot be understood as limiting the present invention. In addition, the terms "first", "second", and the like are only used for descriptive purposes, and cannot be understood as indicating or implying relative importance or implicitly indicating the number of technical features indicated. Thus, features defined as "first", "second", and the like may explicitly or implicitly include one or more of the features. In the description of the present invention, unless otherwise specified, "multiple" means two or more.

[0017] In the description of the present invention, it should be noted that, unless otherwise clearly specified and limited, the terms "installed", "connected", and "connected" should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection, or it can be indirectly connected through an intermediate medium, or it can be the internal communication of two components. For ordinary technicians in this field, the specific meanings of the above terms in the present invention can be understood by specific circumstances.

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

[0019] like Figure 1 The technical roadmap shown in the figure is a method for calculating code deviation based on ground-based GNSS data, comprising the following steps: S1, file traversal of preprocessing module; S2, satellite-by-satellite traversal of the preprocessing module; S3, arc segment by arc segment traversal of the preprocessing module; S4, file traversal of post-processing module; S5, satellite-by-satellite traversal of the post-processing module; S6, arc segment by arc segment traversal of the post-processing module; S7, the solving module solves the DCB using the least square method.

[0020] The present invention reads the navigation file data of ground-based GNSS (Beidou, GPS, GLONASS and GALILEO) and calculates the elevation and azimuth of each navigation satellite based on the approximate coordinates of the ground station. The observation file data is read one by one, the observation values ​​of the corresponding observation type are extracted, the cycle slip is eliminated, and the phase smoothed pseudorange is used for each continuous arc segment of each satellite. After the elevation angle is screened for each epoch, the VTEC containing the DCB error is stored in Matrix, at the same time, calculate the coordinates of the puncture points corresponding to each epoch, and on this basis calculate the various spherical harmonic coefficients and DCB coefficients, and finally solve them using the least squares method , and obtain the spherical harmonic coefficients and DCB. The specific process of the present invention is as follows: a) Preprocessing module 1. File reading Read the navigation file to obtain parameters such as orbit radius sine term, satellite orbit eccentricity and perigee angular distance. Read the observation file to obtain data of the specified observation type.

[0021] 2.SPP solution Based on the read navigation files and observation files, pseudo-range single-point positioning calculation is performed to obtain the coordinates of the navigation satellite and the altitude angle between the satellite and the receiver.

[0022] 3. Altitude & azimuth storage The altitude angle information of different satellites is stored in a certain information format, and the stored information includes the altitude angle and its corresponding satellite number and time.

[0023] 4. Satellite health status assessment Determine whether the SVH value of the current satellite is 0. If not, determine that the current satellite has not passed the preprocessing and jump to the next satellite; otherwise, perform subsequent processing.

[0024] 5. Observation storage Based on the file reading, the pseudorange and carrier of the specified time period are stored according to the time and observation type.

[0025] 6. Cycle slip detection Cycle slips will affect the accuracy of DCB inversion, so the data with cycle slips should be eliminated. The present invention uses Geometry-Free (GF) to detect cycle slips. The formula is as follows: (1) (2) (3) In the formula, is the sequence subscript, , is the wavelength of bands 1 and 2, in meters, , are the phase observations of band 1 and band 2, respectively, in cycles. When , it can be considered that a cycle slip has occurred, and the observation value at the corresponding moment is reset to zero.

[0026] 7. Non-empty arc segment statistics Calculate the lengths of non-zero observation values ​​of different observation types and take their intersection. If the length of the current arc segment is less than 45, count the next arc segment; otherwise, perform subsequent processing.

[0027] 8. Pseudorange smoothing Since the quality of pseudorange is affected by noise and multipath, and the time series is relatively discrete, hatch is used to smooth the pseudoranges of the specified two observation types. The calculation formula is as follows: (4) In the formula, For band The wavelength, in meters, For band Phase observation, in weeks, Indicates that the subscript is The original pseudorange observation value of Indicates that the subscript is The smoothed pseudorange observation value is is the weight value, which depends on the length of the non-empty arc segment.

[0028] 9. Calculate the puncture point coordinates The puncture point is calculated based on the altitude angle and azimuth angle, as shown in formulas (6)-(11): (6) (7) (8) (9) like >70° and ,or <-70° and

[0029] (10) otherwise, (11) 10. Altitude angle screening Data with an altitude angle less than 30° were eliminated.

[0030] If the current satellite and station still have data after the above nine steps, it means that the satellite and the current station have passed the preprocessing, and then the next station is preprocessed.

[0031] b) Post-processing module The satellites and stations that fail the preprocessing are eliminated and the next step of processing is carried out.

[0032] 1. File reading Read the navigation file to obtain the orbit radius sine term, satellite orbit eccentricity, perigee angle, etc. Read the observation file to obtain the specified observation type data.

[0033] 2.SPP solution Based on the read navigation files and observation files, pseudo-range single-point positioning calculation is performed to obtain the coordinates of the navigation satellite and the altitude angle between the satellite and the receiver.

[0034] 3. Altitude & azimuth storage The altitude angle information of different satellites is stored in a certain information format, and the stored information includes the altitude angle and its corresponding satellite number and time.

[0035] 4. Observation storage Based on the file reading, the pseudorange and carrier of the specified time period are stored according to the time and observation type.

[0036] 5. Cycle slip detection Cycle slips will affect the accuracy of DCB inversion, so the data with cycle slips should be eliminated. This paper adopts Geometry-Free (GF) to detect cycle slips, and the formulas are shown in (1) to (3).

[0037] when When , it can be considered that a cycle slip has occurred, and the observation value at the corresponding moment is reset to zero.

[0038] 6. Non-empty arc segment statistics Calculate the lengths of non-zero observation values ​​of different observation types and take their intersection. If the length of the current arc segment is less than 45, count the next arc segment; otherwise, perform subsequent processing.

[0039] 7. Pseudorange smoothing Since the pseudorange quality is affected by noise and multipath, the time series is relatively discrete, so hatch is used to smooth it. The calculation formula is shown in (4).

[0040] 8. Construct P4 matrix The various elements of P4 are the VTEC values ​​before DCB is removed, and the calculation formula is as follows: (5) In the formula, is the altitude angle, is the difference in pseudorange of the specified observation type after smoothing, , Band , The frequency in Hz, is the radius of the Earth, is the puncture point height.

[0041] 9. Calculate the puncture point coordinates The puncture point is calculated based on the altitude angle and azimuth angle, as shown in formulas (6)-(11).

[0042] 10. Construct coefficient matrix B The coefficient matrix B is composed of spherical harmonic coefficients and DCB coefficient modules. The number of columns of spherical harmonic coefficients is the square of the order, as shown in formulas (12)-(13): (12) (13) In the formula, is the latitude of the puncture point, is the longitude of the puncture point, for Step The normalized Legendre coefficients are concluded in matrix B. If If it is not 0, then (13) follows (12).

[0043] The number of columns of the DCB coefficient module is the sum of the number of satellites and the number of stations that have passed the preprocessing. It is composed of satellite DCB and receiver DCB. At the corresponding satellites and stations, the value of the element is as shown in formula (14), and it is 0 at the non-corresponding locations.

[0044] (14) c) Solving module After traversing all files and satellites, the matrix B and matrix P4 are constructed, and the unknown DCB and spherical harmonic coefficients are solved using the least squares method, as shown in formula (15): (15).

[0045] Advantages of the present invention: In order to solve the problem that the ionospheric observation inversion program generally ignores the intraday variation of DCB, thus reducing the accuracy of TEC estimation, a joint DCB solution is performed by simultaneously processing data from multiple satellites and multiple stations within a specified range. More observation information is used to improve the solution accuracy and reliability. At the same time, a preprocessing module is added, and modules such as cycle slip detection and pseudorange smoothing are introduced to enhance the robustness of the results, providing a basis for regional ionospheric modeling.

[0046] Example 1 Taking the observation file of the 364th day of 2024 as an example, the observation values ​​from 0 to 2 are selected for DCB solution. There were 90 ground observation files on that day, and all GPS navigation satellites were selected. The order of the spherical harmonic function was 15. After preprocessing, after S1, S2, and S3 preprocessing, there were 2 GPS satellites and 3 ground observation files that failed the test.

[0047] After S4, S5, and S6, the matrices B and P4 are generated, where the size of matrix B is 103568*342 and the size of P4 is 103568*1. On this basis, the matrix to be solved is initialized , the size is 342*1, of which the first 225 parameters are to be determined, the 226th to 255th parameters are the DCB of the satellite, and the 256th to 342nd parameters are the DCB of the receiver.

[0048] On this basis, the solution module uses the least square method to solve DCB, and the results are shown in Tables 1 to 12 below: Table 1 ; Table 2 ; Table 3 ; Table 4 ; Table 5 ; Table 6 ; Table 7 ; Table 8 ; Table 9 ; Table 10 ; Table 11 ; Table 12 .

[0049] use Calculate the residual, where The modulus is 718, and the sum of the satellite and receiver DCB is 6.4005e-09.

[0050] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for calculating code deviation based on ground-based GNSS data, characterized in that: The steps include: S1, file traversal of preprocessing module; S2, satellite-by-satellite traversal of the preprocessing module; In step S2, the satellite-by-satellite traversal of the preprocessing module includes: S21, satellite health status judgment: judge whether the SVH value of the current satellite is 0. If not, it is determined that the current satellite has not passed the preprocessing and jumps to the next satellite. Otherwise, perform subsequent processing; S22, observation value storage: store the pseudorange and carrier of the specified time period according to the time and observation value type; S23, cycle slip detection: Geometry-Free detection of cycle slips; S3, arc segment by arc segment traversal of the preprocessing module; In step S3, the arc segment by arc segment traversal of the preprocessing module includes: S31, pseudorange smoothing: hatch is used to smooth the pseudoranges of the two specified observation types; S33, elevation angle screening: remove data with elevation angle less than 30°. If the current satellite and station still have data after the above step S31, it is determined that the current satellite and station have passed the preprocessing, and then the next station is preprocessed. S4, file traversal of post-processing module; S5, satellite-by-satellite traversal of the post-processing module; S6, arc segment by arc segment traversal of the post-processing module; In step S6, the arc segment by arc segment traversal of the post-processing module includes: S62, constructing a P4 matrix: each element of P4 is a VTEC value without removing DCB; S64, construct coefficient matrix B: the coefficient matrix B is composed of spherical harmonic coefficients and DCB coefficient modules, and the number of columns of the spherical harmonic coefficients is the square of the order; The number of columns in the DCB coefficient module is the sum of the number of satellites and stations that have passed preprocessing; S7, the solution module uses the least square method to solve the DCB; In step S7, the solving module solves the DCB using the least square method, including: After the coefficient matrix B and the matrix P4 are constructed, the least squares method is used to solve the unknown DCB and the spherical harmonic coefficients.

2. The method for calculating code deviation based on ground-based GNSS data according to claim 1, characterized in that: In step S1, the file traversal of the preprocessing module includes: S11, file reading: read the navigation file, obtain the orbit radius sine term, satellite orbit eccentricity and perigee angular distance parameters, read the observation file, and obtain the specified observation type data; S12, SPP solution: Based on the read navigation files and observation files, perform pseudo-range single-point positioning calculation to obtain the coordinates of the navigation satellite and the altitude angle between the satellite and the receiver; S13, altitude angle & azimuth angle storage: the altitude angle information of different satellites is stored according to the information format, and the stored information includes the altitude angle and its corresponding satellite number and time.

3. The method for calculating code deviation based on ground-based GNSS data according to claim 1, characterized in that: In step S23, the formula for cycle slip detection is as follows: (1) (2) (3) In the formula, is the sequence subscript, , are the wavelengths of bands 1 and 2 respectively, in meters. , They are the phase observations of bands 1 and 2, in cycles. When , it is determined that a cycle slip has occurred and the observation value at the corresponding time is cleared to zero.

4. The method for calculating code deviation based on ground-based GNSS data according to claim 1, characterized in that: In step S2, the satellite-by-satellite traversal of the pre-processing module further includes: S24, non-empty arc statistics: calculate the length of observation values ​​of different observation types that are not 0, and take their intersection. If the current arc length is less than 45, count the next arc, otherwise, perform subsequent processing.

5. The method for calculating code deviation based on ground-based GNSS data according to claim 1, characterized in that: In step S31, the formula for pseudorange smoothing is as follows: (4) In the formula, For band The wavelength, in meters, For band Phase observation, in weeks, Indicates that the subscript is The original pseudorange observation value of Indicates that the subscript is The smoothed pseudorange observation value is is the weight value, which depends on the length of the non-empty arc segment.

6. The method for calculating code deviation based on ground-based GNSS data according to claim 1, characterized in that: In step S3, the arc segment-by-arc traversal of the preprocessing module further includes: S32, puncture point calculation: Calculate the puncture point based on the altitude angle and azimuth angle. The formula is as follows: (6) (7) (8) (9) like >70° and ,or <-70° and , (10) otherwise, The value of is shown in (11): (11) In the formula, is the altitude angle, is the latitude of the puncture point, is the longitude of the puncture point, is the receiver latitude, is the receiver longitude, is the azimuth, is the angle formed by the receiver, the center of the earth, and the puncture point. is the angle between the center of the earth - the puncture point and the receiver - the satellite .

7. The method for calculating code deviation based on ground-based GNSS data according to claim 2, characterized in that: In step S4, the files of the post-processing module are traversed, and the steps are the same as those of step S1.

8. A method for calculating code deviation based on ground-based GNSS data according to claim 3 or 4, characterized in that: In step S5, the satellite-by-satellite traversal of the post-processing module includes: S51, observation value storage, the steps are the same as step S22; S52, cycle slip detection, the steps are the same as step S23; S53, non-empty arc segment statistics, the steps are the same as step S24.

9. A method for calculating code deviation based on ground-based GNSS data according to claim 5 or 6, characterized in that: In step S6, the arc segment-by-arc segment traversal of the post-processing module further includes: S61, pseudorange smoothing, the steps are the same as step S31; S63, puncture point calculation, the steps are the same as step S32; In step S62, the calculation formula for constructing the P4 matrix is ​​as follows: (5) In the formula, is the altitude angle, is the difference in pseudorange of the specified observation type after smoothing, , Band , The frequency in Hz, is the radius of the Earth, is the puncture point height; In step S64, the formula for constructing the coefficient matrix B is as follows: (12) (13) In the formula, is the latitude of the puncture point, is the longitude of the puncture point, for Step The normalized Legendre coefficients are then concluded. In the coefficient matrix B, if If it is not 0, then formula (13) follows formula (12); The DCB coefficient module column number is composed of satellite DCB and receiver DCB. At the corresponding satellite and corresponding station, the value of the element is as shown in the following formula, and the non-corresponding position is 0; (14)； In the formula, is the altitude angle, , Band , The frequency in Hz, is the speed of light, is the radius of the Earth, is the puncture point height.

10. The method for calculating code deviation based on ground-based GNSS data according to claim 1, characterized in that: In step S7, the solution module uses the least square method to solve the DCB formula as follows: (15)； In the formula, To solve the result, we need to calculate the spherical harmonic coefficients, the satellite DCB and the receiver DCB. is the coefficient matrix, is the observation matrix.

Images