A method and system for determining Earth rotation parameters based on GNSS data stream

Abstract

This invention provides a method and system for determining Earth rotation parameters based on GNSS data streams. The method involves adding estimation of Earth rotation parameters during the precise orbit determination process using square root information filtering. In time updates, the acceleration values ​​of polar motion and UT1 variation are treated as white noise to form a state equation. In measurement updates, UT1 values ​​from the ultra-fast ERP forecast product provided by IGS are added as virtual observations to form a virtual observation equation. This decouples the UT1 parameters from the orbital parameters, allowing them to be estimated. This achieves stable estimation of real-time Earth rotation parameters with accuracy superior to Earth rotation parameter forecast products provided by IGS and IERS, while also obtaining more accurate precise orbital products, better serving the real-time, high-precision ERP application requirements.

Description

Technical Field

[0001] This invention relates to the field of satellite measurement technology, and in particular to a method and system for determining Earth rotation parameters based on GNSS data streams. Background Technology

[0002] Earth Rotation Parameters (ERPs) include polar motion (PM), polar motion rate, Universal Time (UT1), and Length-of-Day (LOD). Together with precession and nutation, they constitute the Earth Orientation Parameter (EOP). These are essential physical parameters for converting between the Terrestrial Reference Frame (TRF) and the Celestial Reference Frame (CRF), and are of great significance for applications such as the establishment and maintenance of the Earth Reference Frame, deep space exploration, precise orbit determination of navigation satellites, and climate change analysis.

[0003] VLBI, GNSS, SLR, and DORIS are currently the most accurate techniques for determining ERP (Earthquake Performance). However, due to differences in sampling intervals and data processing precision and complexity among these different observation techniques, the accurate determination of ERP can be delayed by hours to days. Real-time ERP acquisition typically relies on forecasting. In the International GNSS Service (IGS) Ultra-rapid Products (IGU) ERP products, the forecasting component has a four-fold reduction in accuracy compared to the calculated component. Currently, IGS provides real-time GNSS observation data streams from hundreds of globally distributed stations, with second-level update frequencies, and broadcasts using the RTCM (Network Transmission Protocol Transmitted via IP) protocol. This allows for real-time estimation and updates of GNSS satellite orbits, clock errors, and ERP products.

[0004] Therefore, a new scheme is needed to estimate the Earth's rotation parameters in real time in order to obtain a more accurate real-time ERP product. Summary of the Invention

[0005] This invention provides a method and system for determining Earth rotation parameters (ERPs) based on GNSS data streams. The ERP determination using GNSS technology is achieved by determining the transformation matrix between the Earth reference frame and the celestial reference frame. The coordinates of the reference point in the Earth reference frame are determined by the ground station, while the satellite coordinates in the celestial reference frame are determined by the satellite's dynamic orbit. Therefore, using a station-satellite pair network for calculation allows for the simultaneous calculation of parameters such as satellite orbit, clock error, ambiguity, troposphere, and ERP.

[0006] In a first aspect, the present invention provides a method for determining Earth rotation parameters based on GNSS data streams, comprising:

[0007] Step 1: Obtain the initial state of satellites and parameters to be estimated at each station based on the collected GNSS real-time observation data stream;

[0008] Step 2: Obtain the initial value of the orbit position and perform numerical integration on the orbit to obtain the state transition parameters of the preceding and following epochs;

[0009] Step 3: Update the parameters in various state transition models to obtain the updated information matrix;

[0010] Step 4: Based on the initial state of the satellites and the parameters to be estimated at each station, construct the observation equations;

[0011] Step 5: Correct the observation error in the observation equation using the preset model to obtain the error equation after standardization of residuals;

[0012] Step 6: Fuse prior information of state variables with observation information, and measure and update the information matrix for the next time step;

[0013] Step 7: After the filtering converges, the ambiguity parameter is quickly fixed using the non-differential ambiguity fixing method. The fixed ambiguity value is added to the observation equation of the filtering solution for constraint. Otherwise, skip this step and proceed to step 8.

[0014] Step 8: Solve the information equation to obtain the corrected values ​​of all parameters at the next time step, and superimpose them with the initial values ​​after the time update to obtain the estimated values ​​of all parameters at the current time step;

[0015] Step 9: Determine whether all epochs have been processed. If not, update the initial state of the parameters and repeat steps 2 to 9. If yes, end the process.

[0016] According to the method for determining Earth rotation parameters based on GNSS data stream provided by the present invention, step 1 includes:

[0017] The coordinates of each station in the Earth reference system are regarded as constants. The method is to constrain them to a specified GNSS data format solution that is updated periodically, and obtain the parameters to be estimated, including Earth rotation parameters, satellite position, velocity, light pressure model parameters, satellite clock error, station clock error, ambiguity, station zenith tropospheric delay and inter-system deviation.

[0018] The Earth's rotation parameters include polar motion in the x and y directions, polar motion rates in the x and y directions, the difference between UT1 and UTC, and the variation in day length, denoted as... .

[0019] According to the method for determining Earth rotation parameters based on GNSS data stream provided by the present invention, step 2 includes:

[0020] Determine the equations of motion of the satellite in the inertial coordinate system:

[0021]

[0022] in, It is the position vector of the satellite's center of mass. It's the satellite's speed. These are the dynamic parameters to be estimated in the satellite dynamics equations. , and These are the conservative forces, non-conservative forces, and unmodeled empirical perturbations acting on the satellite;

[0023] The equations of motion are solved using numerical integration.

[0024] According to the method for determining Earth rotation parameters based on GNSS data stream provided by the present invention, step 3 includes:

[0025] For the Earth's rotation parameters, the acceleration values ​​of polar motion and UT1 variation are treated as white noise, and the state transition matrix of the previous and next epochs is derived based on the white noise:

[0026]

[0027] in and Representing two epochs, It is the state transition matrix. It can characterize process noise covariance matrix Then we have:

[0028]

[0029]

[0030] in Represents a 3rd order identity matrix. It is a 3rd order zero matrix. Given the time interval between consecutive epochs, we can derive:

[0031]

[0032] in It is a 3x3 diagonal matrix. The diagonal elements are the polar migration velocities in the x and y directions and the variance of the diurnal variation noise per unit time, respectively.

[0033] According to the method for determining Earth rotation parameters based on GNSS data stream provided by the present invention, step 4 includes:

[0034] Raw GNSS observations are obtained from real-time streaming data and combined into ionospherically-free composite observations after gross error detection:

[0035]

[0036] in, and These are pseudorange and phase observations, respectively. It is the geometric distance between the station and the satellite. and These are the station and satellite clock biases, It is a tropospheric delay. For ambiguity, The wavelength of the phase combination without ionosphere and These are unmodeled residuals.

[0037] According to the method for determining Earth rotation parameters based on GNSS data stream provided by the present invention, step 5 includes:

[0038] Regarding the Earth's rotation parameters, we have:

[0039]

[0040] in It is a rotation matrix. , , and These represent precession, nutation, Earth's rotation, and polar motion, respectively. These are the station coordinates in the Earth's reference frame, and are considered constants here. These are satellite coordinates in the celestial coordinate system. Linearizing the above formula yields:

[0041]

[0042] in , ,and They are , and The initial values ​​are discussed only for ERP parameters, denoted as... Then we have:

[0043]

[0044] in Representing the rotation matrix Given the initial value of , the partial derivative on the right side of the equation can be expressed as:

[0045]

[0046] in and These are the ERP parameters and initial value, It is the Earth's rotation angle, which can be calculated by the following formula:

[0047]

[0048] in It is the Julian Day corresponding to UT1 time;

[0049] For polar migration rate and diurnal variation, increase the derivative with respect to time. ;

[0050] Based on the UT1 forecast values ​​provided by the IGS or IERS service organizations, the virtual observation equation added for the UT1 values ​​is as follows:

[0051]

[0052] in The UT1 value obtained from the time update. The UT1 value represents the external constraint and is obtained using the ERP forecast product provided by IGS or IERS.

[0053] Secondly, the present invention also provides a system for determining Earth rotation parameters based on GNSS data streams, comprising:

[0054] The data acquisition module is used to acquire the initial state of satellites and parameters to be estimated at each station based on the collected GNSS real-time observation data stream.

[0055] The orbit integration module is used to obtain the initial value of the orbit position and integrate the orbit values ​​to obtain the state transition parameters of the preceding and following epochs;

[0056] The time update module is used to update the parameters in various state transition models to obtain the updated information matrix.

[0057] The module is used to construct observation equations based on the initial states of satellites and parameters to be estimated at each station;

[0058] The correction module is used to correct the observation error in the observation equation using a preset model, and obtain the error equation after standardization of residuals.

[0059] The fusion module is used to fuse prior information of state variables with observation information, and measure and update the information matrix for the next time step.

[0060] The constraint module is used to quickly fix the ambiguity parameters after the filtering converges using the non-differential ambiguity fixing method. The fixed ambiguity value is added to the observation equation of the filtering solution for constraint. Otherwise, this step is skipped and the execution steps corresponding to the solution module are executed.

[0061] The solution module is used to solve the information equation to obtain the corrected values ​​of all parameters at the next time step, and then superimpose them with the initial values ​​after the time update to obtain the estimated values ​​of all parameters at the current time step.

[0062] The judgment module is used to determine whether all epochs have been processed. If not, the initial state of the parameters is updated, and the orbit integration module is executed repeatedly until the corresponding execution step of the judgment module is reached. If yes, the process ends.

[0063] Thirdly, the present invention also provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the method for determining Earth rotation parameters based on GNSS data stream as described above.

[0064] Fourthly, the present invention also provides a non-transitory computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the method for determining Earth rotation parameters based on GNSS data streams as described above.

[0065] Fifthly, the present invention also provides a computer program product, including a computer program that, when executed by a processor, implements the method for determining Earth rotation parameters based on GNSS data stream as described above.

[0066] The method and system for determining Earth rotation parameters based on GNSS data stream provided by this invention achieve real-time estimation and updating of Earth rotation parameters through real-time GNSS observation data streams, resulting in a high-precision ERP product with high update frequency and low latency. The update interval and time delay can reach the minute level, and the accuracy is nearly doubled compared to the ERP forecast product provided by IGU. This invention effectively solves the pain point of traditional Earth rotation parameter determination methods where high accuracy and low latency are mutually exclusive, and can better serve the application needs of high-precision real-time ERP. Attached Figure Description

[0067] To more clearly illustrate the technical solutions in this invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0068] Figure 1 This is one of the flowcharts illustrating the method for determining Earth rotation parameters based on GNSS data stream provided by the present invention;

[0069] Figure 2 This is the second flowchart of the method for determining Earth rotation parameters based on GNSS data stream provided by the present invention;

[0070] Figure 3 This is a schematic diagram of the Earth rotation parameter determination system based on GNSS data stream provided by the present invention;

[0071] Figure 4 This is a schematic diagram of the structure of the electronic device provided by the present invention. Detailed Implementation

[0072] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this 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 this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

[0073] To address the limitations of existing technologies, this invention utilizes GNSS technology to determine the ERP by determining the transformation matrix between the Earth reference system and the celestial reference system. During the precise orbit determination process using a satellite-station network, the ERP can be calculated synchronously. The technical solution of this invention will be further explained below with reference to the accompanying drawings and specific embodiments.

[0074] Figure 1This is one of the flowcharts illustrating the method for determining Earth rotation parameters based on GNSS data streams provided in this embodiment of the invention. Figure 1 As shown, it includes:

[0075] Step 1: Obtain the initial state of satellites and parameters to be estimated at each station based on the collected GNSS real-time observation data stream;

[0076] Step 2: Obtain the initial value of the orbit position and integrate the orbital values ​​to obtain the state transition parameters for the preceding and following epochs;

[0077] Step 3: Update the parameters in various state transition models to obtain the updated information matrix;

[0078] Step 4: Based on the initial state of the satellites and the parameters to be estimated at each station, construct the observation equations;

[0079] Step 5: Correct the observation error in the observation equation using the preset model to obtain the error equation after standardization of residuals;

[0080] Step 6: Fuse prior information of state variables with observation information, and measure and update the information matrix for the next time step;

[0081] Step 7: After the filtering converges, the ambiguity parameter is quickly fixed using the non-differential ambiguity fixing method. The fixed ambiguity value is added to the observation equation of the filtering solution for constraint. Otherwise, skip this step and proceed to step 8.

[0082] Step 8: Solve the information equation to obtain the corrected values ​​of all parameters at the next time step, and superimpose them with the initial values ​​after the time update to obtain the estimated values ​​of all parameters at the current time step;

[0083] Step 9: Determine whether all epochs have been processed. If not, update the initial state of the parameters and repeat steps 2 to 9. If yes, end the process.

[0084] Specifically, combined Figure 2 The logical route shown includes:

[0085] Step 1, Initialization, involves GNSS data preparation and preprocessing, including real-time stream decoding and transfer, preparation and data reading of tables required by various error models. This obtains the initial state of satellites and parameters to be estimated at each station, providing initial values ​​for filtering estimation. The coordinates of the stations in the Earth reference system are considered constants, strongly constrained to periodically updated SINEX solutions or post-hoc precise single-point positioning static solutions. Parameters to be estimated include Earth rotation parameters, satellite positions and velocities, light pressure model parameters, satellite clock errors, station clock errors, ambiguities, and station zenith tropospheric delay. In GNSS multi-system joint processing, inter-system biases are also included. Earth rotation parameters include polar motion (x and y directions), polar motion rate (x and y directions), the difference between UT1 and UTC, and diurnal variation, totaling six parameters, denoted as... .

[0086] Step 2: Orbit integration to obtain the initial value of the orbital position and the state transition parameters for the preceding and following epochs. The equations of motion of the satellite in the inertial coordinate system can be described by a set of first-order differential equations:

[0087]

[0088] in, It is the position vector of the satellite's center of mass. It's the satellite's speed. These are the dynamic parameters to be estimated in the satellite dynamics equations. , and These are the conservative forces, non-conservative forces, and unmodeled empirical perturbations acting on the satellite. Typically, numerical integration methods can be used to solve this differential equation, often employing a combination of the RKF single-step method and the Adams multi-step method.

[0089] Step 3, Time Update: For the Earth's rotation parameters, the acceleration values ​​of polar motion and UT1 variation are treated as white noise to derive the state transition matrix for the preceding and following epochs. That is, the Earth's rotation parameters for the preceding and following epochs satisfy the following formula:

[0090]

[0091] in and Representing two epochs, It is the state transition matrix. It can characterize process noise covariance matrix Then we have:

[0092]

[0093]

[0094] in Represents a 3rd order identity matrix. It is a 3rd order zero matrix. Let be the time interval between consecutive epochs. Therefore, we can derive:

[0095]

[0096] in It is a 3x3 diagonal matrix, with diagonal elements representing the polar migration velocities (x and y directions) and the variance of the diurnal variation noise per unit time, respectively. Time series analysis of the ERP series shows that at a correlation time of 300 seconds... Set to 1e-6 1e-6 and 1e-6 .

[0097] Step 4: Initialize the observation equations. Read GNSS observations and construct the observation equations. Obtain raw GNSS observations from real-time streaming data, and after gross error detection, assemble ionospherically-free composite observations:

[0098]

[0099] in, and These are pseudorange and phase observations, respectively. It is the geometric distance between the station and the satellite. and These are the station and satellite clock biases, It is a tropospheric delay. For ambiguity, The wavelength of the phase combination without ionosphere and These are unmodeled residuals.

[0100] Step 5: Formulate the error equations. Correct the observation errors using common models, and then linearize the observation equations for the parameters from Step 1. For the Earth's rotation parameters, we have:

[0101]

[0102] in It is a rotation matrix. , , and These represent precession, nutation, Earth's rotation, and polar motion, respectively. These are the station coordinates in the Earth's reference frame, and are considered constants here. These are the satellite coordinates in the celestial coordinate system. Linearizing formula (7), we get:

[0103]

[0104] in , ,and They are , and The initial values ​​are then set. Here, we will only discuss the ERP parameters, and denote them as follows: Then we have:

[0105]

[0106] in Representing the rotation matrix The initial value of . The partial derivative on the right side of the equation can be expressed as:

[0107]

[0108] in and These are the ERP parameters and initial value, It is the Earth's rotation angle, which can be calculated by the following formula:

[0109]

[0110] in This is the Julian day corresponding to UT1 time. For the variations in polar migration rate and day length, only the derivative with respect to time needs to be added. That's it.

[0111] Furthermore, regarding the UT1 value in the Earth's rotation parameters, it is impossible to estimate due to its coupling with orbital parameters during GNSS orbit determination. Therefore, the UT1 prediction value provided by the International GNSS Service (IGS) or the International Earth Rotation Service (IERS) is added as a virtual observation to form a virtual observation equation, making the UT1 parameter estimable. The added virtual observation equation is as follows:

[0112]

[0113] in The UT1 value obtained from the time update. The UT1 value, representing external constraints, can be obtained from ERP forecast products provided by IGS or IERS.

[0114] Step 6: Fixing Undifferentiated Ambiguities. After the filtering converges, the undifferentiated ambiguities are fixed using the calculated wide-lane and narrow-lane uncalibrated phase delay (UPD) information. The fixed ambiguities are then added as constraint information to the observation equations to obtain the final ambiguity-fixed solution. Otherwise, this step is skipped.

[0115] Step 7, Measurement Update. The prior information of the state variables is fused with the observation information to obtain the measurement update. Information matrix at any given time.

[0116] Step 8, parameter recovery. Solve the information equation in step 7 to obtain... The corrected values ​​for all parameters at time point are then superimposed onto the initial values ​​after the time update to obtain the estimated values ​​for all parameters at that time point.

[0117] Step 9: Determine if all epochs have been processed. If not, update the initial parameter state and repeat steps 2 to 9; if yes, terminate normally.

[0118] In one embodiment, an experiment was conducted to estimate the ERP in real time using 32 days of GPS data from 120 globally distributed stations with real-time data streams, spanning days 309 to 340 of 2023. The ERP parameters were updated during real-time filtering and orbit determination, with a data processing interval of 300 seconds. To ensure the accuracy of ambiguity fixing, unequal ambiguity fixing was performed two days after filtering began. Ambiguity fixing yielded a more accurate ERP estimate, with the ERP estimate sequence time interval being 300 seconds. Due to the time required for data decoding and computation, the time delay for ERP estimation was approximately 60 seconds.

[0119] The accuracy of the aforementioned ERP solution sequence was evaluated using the IERS 20C04 post-hoc refined product as a reference. The IERS 20C04 product is currently recognized as the most accurate ERP product, with a time delay of approximately one month, providing ERP values ​​at 0:00 and 12:00 within a day. The sampling interval of the real-time ERP estimation product sequence is 300 seconds. To avoid accuracy loss of the reference product due to interpolation, the real-time estimated ERP solution sequence was thinned to the IERS 20C04 product time and the difference was used for accuracy evaluation. For comparative analysis, the accuracy of the real-time available forecast products (IERS eopc04_extended) of IGS and IERS during the same period was also evaluated using the same method. The results are shown in Table 1.

[0120] Table 1. Statistics on the accuracy (RMS) of IERS and IGU forecast products and ERP product sequences calculated by real-time filtering.

[0121]

[0122] As can be seen, the ERP sequence obtained by real-time filtering is... , UT1 , and The RMS values ​​are 80. 73 , twenty three , 238 , 219 and 44 Compared to real-time IGS and IERS forecast products, it has a significant advantage in accuracy.

[0123] In one embodiment, the accuracy of the ERP parameters affects the calculation accuracy of the real-time orbit product in real-time filtered precise orbit determination. This example compares the results of real-time filtered precise orbit determination under different ERP processing strategies and designs three experimental processing schemes. Scheme 1: Do not estimate the ERP parameters and directly use the ultra-fast forecast product (IGU) provided by IGS; Scheme 2: Do not estimate the ERP parameters and directly use the post-precision product of IERS 20C04; Scheme 3: Estimate the ERP parameters synchronously according to the method of this patent. The GNSS real-time stream data used in this example is the same as in Example 1, and the post-precision orbit product released by IGS is used as a reference to evaluate the calculation accuracy of the real-time orbit. Table 2 shows the RMS statistics of the real-time filtered precise orbit determination solution products of the three schemes in the tangential (A), normal (C), radial (R), and three-dimensional (3D) directions.

[0124] Table 23 shows the statistical accuracy (RMS) of 30-day GPS real-time filtering precision orbit determination products under 23 ERP processing strategies (compared with IGS final precision products).

[0125]

[0126] As can be seen, the results of Scheme 1, which directly uses the IGU-predicted ERP product for orbit determination, are significantly poor, with an RMS value of 7.9 cm in the three-dimensional direction. Comparing Scheme 2 and Scheme 3, the method of this patent can obtain more accurate orbit results, with a three-dimensional orbit accuracy improvement of over 1 cm. In summary, it can be seen that the real-time ERP estimation method of this patent can effectively improve the solution accuracy of real-time filtered precision orbit determination.

[0127] The Earth rotation parameter determination system based on GNSS data stream provided by this invention is described below. The Earth rotation parameter determination system based on GNSS data stream described below can be referred to in correspondence with the Earth rotation parameter determination method based on GNSS data stream described above.

[0128] Figure 3This is a schematic diagram of the structure of the Earth rotation parameter determination system based on GNSS data stream provided in an embodiment of the present invention, as shown below. Figure 3 As shown, it includes: a data acquisition module 31, an orbit integration module 32, a time update module 33, a construction module 34, a correction module 35, a fusion module 36, a constraint module 37, a solution module 38, and a judgment module 39, wherein:

[0129] The data acquisition module 31 is used to acquire the initial state of satellites and parameters to be estimated at each station based on the collected GNSS real-time observation data stream;

[0130] The orbit integration module 32 is used to obtain the initial value of the orbit position and integrate the orbit values ​​to obtain the state transition parameters of the preceding and following epochs;

[0131] The time update module 33 is used to update the parameters in various state transition models to obtain the updated information matrix;

[0132] Module 34 is used to construct observation equations based on the initial states of satellites and parameters to be estimated at each station;

[0133] The correction module 35 is used to correct the observation error in the observation equation using a preset model, and obtain the error equation after standardization of residuals;

[0134] The fusion module 36 is used to fuse prior information of state variables with observation information, and to measure and update the information matrix for the next time step.

[0135] The constraint module 37 is used to quickly fix the ambiguity parameters using the non-differential ambiguity fixing method after the filtering converges. The fixed ambiguity value is added to the observation equation of the filtering solution for constraint. Otherwise, this step is skipped and the execution step corresponding to the solution module is executed.

[0136] Solver module 38 is used to solve the information equation to obtain the corrected values ​​of all parameters at the next time step, and then superimpose them with the initial values ​​after the time update to obtain the estimated values ​​of all parameters at the current time step.

[0137] The judgment module 39 is used to determine whether all epochs have been processed. If not, the initial state of the parameters is updated and the orbit integration module is executed repeatedly until the corresponding execution step of the judgment module is reached. If yes, the process ends.

[0138] Figure 4 An example is a schematic diagram of the physical structure of an electronic device, such as... Figure 4As shown, the electronic device may include: a processor 410, a communication interface 420, a memory 430, and a communication bus 440. The processor 410, communication interface 420, and memory 430 communicate with each other via the communication bus 440. The processor 410 can call logical instructions in the memory 430 to execute a method for determining Earth rotation parameters based on GNSS data streams. This method includes: Step 1: Obtaining the initial states of satellites and parameters to be estimated at each station based on the collected real-time GNSS observation data stream; Step 2: Obtaining initial orbital position values ​​and integrating the orbital values ​​to obtain state transition parameters for previous and subsequent epochs; Step 3: Updating the state of each parameter in various state transition models to obtain an updated information matrix; Step 4: Constructing observation equations based on the initial states of satellites and parameters to be estimated at each station; Step 5: Correcting observation errors in the observation equations using a preset model to obtain a standard... Step 6: Fuse the prior information of the state variables with the observation information, and measure and update the information matrix for the next time step; Step 7: After the filtering converges, use the non-differential ambiguity fixing method to quickly fix the ambiguity parameters, and add the fixed ambiguity values ​​to the observation equation of the filtering solution for constraint. Otherwise, skip this step and proceed to Step 8; Step 8: Solve the information equation to obtain the correction values ​​of all parameters for the next time step, and superimpose them with the initial values ​​after the time update to obtain the estimated values ​​of all parameters for the current time step; Step 9: Determine whether all epochs have been processed. If not, update the initial state of the parameters and repeat steps 2 to 9. If yes, end the process.

[0139] Furthermore, the logical instructions in the aforementioned memory 430 can be implemented as software functional units and, when sold or used as independent products, can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, or the part that contributes to the prior art, or a part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0140] The device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Those skilled in the art can understand and implement this without any creative effort.

[0141] Through the above description of the embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus necessary general-purpose hardware platforms, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solutions, in essence 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 ROM / RAM, magnetic disk, optical disk, etc., and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute the methods described in the various embodiments or some parts of the embodiments.

[0142] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for determining Earth's rotation parameters based on GNSS data stream, characterized in that, include: Step 1: Obtain the initial state of satellites and parameters to be estimated at each station based on the collected GNSS real-time observation data stream; Step 2: Obtain the initial value of the orbit position and integrate the orbital values ​​to obtain the state transition parameters for the preceding and following epochs; Step 3: By updating the parameters in various state transition models, the updated information matrix is ​​obtained, including: Regarding the Earth's rotation parameters, the acceleration values ​​of polar motion and UT1 variation are treated as white noise. Based on the white noise, the state transition matrix of the preceding and following epochs is derived. The Earth's rotation parameters of the preceding and following epochs satisfy the following formula: in and Representing two epochs, It is the state transition matrix. Used to characterize process noise covariance matrix ; Step 4: Based on the initial states of satellites and parameters to be estimated at each station, construct the observation equations, including: Raw GNSS observations are obtained from real-time streaming data and combined into ionospherically-free composite observations after gross error detection: in, and These are pseudorange and phase observations, respectively. It is the geometric distance between the station and the satellite. and These are the station and satellite clock biases, It is a tropospheric delay. For ambiguity, The wavelength of the phase combination without ionosphere and These are unmodeled residuals; Step 5: Correct the observation error in the observation equation using the preset model to obtain the error equation after standardization of residuals; Based on the UT1 forecast values ​​provided by the IGS or IERS service organizations, the virtual observation equation added for the UT1 values ​​is as follows: in The UT1 value obtained from the time update. The UT1 value, representing external constraints, is obtained using ERP forecast products provided by IGS or IERS. Step 6: Fuse prior information of state variables with observation information, and measure and update the information matrix for the next time step; Step 7: After the filtering converges, the ambiguity parameter is quickly fixed using the non-differential ambiguity fixing method. The fixed ambiguity value is added to the observation equation of the filtering solution for constraint. Otherwise, skip this step and proceed to step 8. Step 8: Solve the information equation to obtain the corrected values ​​of all parameters at the next time step, and superimpose them with the initial values ​​after the time update to obtain the estimated values ​​of all parameters at the current time step; Step 9: Determine whether all epochs have been processed. If not, update the initial state of the parameters and repeat steps 2 to 9. If yes, end the process.

2. The method for determining Earth rotation parameters based on GNSS data stream according to claim 1, characterized in that, Step 1 includes: The coordinates of each station in the Earth reference system are regarded as constants. The method is to constrain them to a specified GNSS data format solution that is updated periodically, and obtain the parameters to be estimated, including Earth rotation parameters, satellite position, velocity, light pressure model parameters, satellite clock error, station clock error, ambiguity, station zenith tropospheric delay and inter-system deviation. The Earth's rotation parameters include polar motion in the x and y directions, polar motion rates in the x and y directions, the difference between UT1 and UTC, and the variation in day length, denoted as... .

3. The method for determining Earth rotation parameters based on GNSS data stream according to claim 1, characterized in that, Step 2 includes: Determine the equations of motion of the satellite in the inertial coordinate system: in, It is the position vector of the satellite's center of mass. It's the satellite's speed. These are the dynamic parameters to be estimated in the satellite dynamics equations. , and These are the conservative forces, non-conservative forces, and unmodeled empirical perturbations acting on the satellite; The equations of motion are solved using numerical integration.

4. The method for determining Earth rotation parameters based on GNSS data stream according to claim 1, characterized in that, Step 3 includes: in Represents a 3rd order identity matrix. It is a 3rd order zero matrix. Given the time interval between consecutive epochs, we can derive: in It is a 3x3 diagonal matrix. The diagonal elements are the polar migration velocities in the x and y directions and the variance of the diurnal variation noise per unit time, respectively.

5. The method for determining Earth rotation parameters based on GNSS data stream according to claim 1, characterized in that, Step 5 includes: Regarding the Earth's rotation parameters, we have: in It is a rotation matrix. , , and These represent precession, nutation, Earth's rotation, and polar motion, respectively. These are the station coordinates in the Earth's reference frame, and are considered constants here. These are satellite coordinates in the celestial coordinate system. Linearizing the above formula yields: in , ,and They are , and The initial values ​​are discussed only for ERP parameters, denoted as... Then we have: in Representing the rotation matrix Given the initial value of , the partial derivative on the right side of the equation can be expressed as: in and These are the ERP parameters and initial value, It is the Earth's rotation angle. It can be calculated using the following formula: in It is the Julian Day corresponding to UT1 time; For polar migration rate and diurnal variation, increase the derivative with respect to time. .

6. A system for determining Earth rotation parameters based on GNSS data stream, based on the method for determining Earth rotation parameters based on GNSS data stream according to any one of claims 1 to 5, characterized in that, include: The data acquisition module is used to acquire the initial state of satellites and parameters to be estimated at each station based on the collected GNSS real-time observation data stream. The orbit integration module is used to obtain the initial value of the orbit position and integrate the orbit values ​​to obtain the state transition parameters of the preceding and following epochs; The time update module is used to update the parameters in various state transition models to obtain the updated information matrix. The module is used to construct observation equations based on the initial states of satellites and parameters to be estimated at each station; The correction module is used to correct the observation error in the observation equation using a preset model, and obtain the error equation after standardization of residuals. The fusion module is used to fuse prior information of state variables with observation information, and measure and update the information matrix for the next time step. The constraint module is used to quickly fix the ambiguity parameters after the filtering converges using the non-differential ambiguity fixing method. The fixed ambiguity value is added to the observation equation of the filtering solution for constraint. Otherwise, this step is skipped and the execution steps corresponding to the solution module are executed. The solution module is used to solve the information equation to obtain the corrected values ​​of all parameters at the next time step, and then superimpose them with the initial values ​​after the time update to obtain the estimated values ​​of all parameters at the current time step. The judgment module is used to determine whether all epochs have been processed. If not, the initial state of the parameters is updated, and the orbit integration module is executed repeatedly until the corresponding execution step of the judgment module is reached. If yes, the process ends.

7. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the method for determining Earth rotation parameters based on GNSS data stream as described in any one of claims 1 to 5.

8. A non-transitory computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the method for determining Earth rotation parameters based on GNSS data stream as described in any one of claims 1 to 5.

9. A computer program product, comprising a computer program, characterized in that, When the computer program is executed by the processor, it implements the method for determining Earth rotation parameters based on GNSS data stream as described in any one of claims 1 to 5.

Images