Earth rotation parameter determination method and system based on GNSS data stream

Abstract

The invention provides an earth rotation parameter determination method and system based on a GNSS data stream, and the method comprises the steps: achieving the purpose of increasing the estimation of earth rotation parameters in a process of achieving precise orbit determination through employing square root information filtering, regarding the acceleration values of polar shift and UT1 change as white noise in time updating, and forming a state equation; in the measurement updating process, a UT1 value in an ultra-fast ERP forecasting product provided by the IGS is added to serve as virtual observation to form a virtual observation equation, so that UT1 parameters and orbit parameters are decoupled and can be estimated, stable estimation of real-time earth rotation parameters is achieved, the precision is superior to that of earth rotation parameter forecasting products provided by the IGS and IERS, and the accuracy of the earth rotation parameter forecasting products provided by the IGS and the IERS is improved. And meanwhile, a precise track product with higher precision is obtained, and the real-time high-precision ERP application requirement can be better met.

Description

Technical Field

[0001] The present 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 stream. Background Art

[0002] The Earth Rotation Parameters (ERP) 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). They are necessary physical parameters for the conversion between the Terrestrial Reference Frame (TRF) and the Celestial Reference Frame (CRF). They are of great significance to the establishment and maintenance of the Earth Reference Frame, deep space exploration, precise orbit determination of navigation satellites, climate change and other applications.

[0003] VLBI, GNSS, SLR, DORIS and other space geodetic technologies are currently the most accurate means of measuring ERP. Due to the different sampling intervals of different observation technologies, as well as the differences in data processing accuracy and complexity, there is a delay of several hours to several days in the accurate measurement of ERP. Real-time ERP is usually obtained by forecasting, and the ultra-rapid (IGU) ERP products provided by the International GNSS Service (IGS) have a forecast part that is 4 times less accurate than the calculation part. Currently, IGS provides GNSS real-time observation data streams from hundreds of globally distributed stations, with an update frequency of seconds, and broadcasts using the RTCM Internet Transport Protocol (NTRIP), which can realize real-time estimation and update of GNSS satellite orbits, clock errors and ERP products.

[0004] Therefore, it is necessary to propose a new scheme for real-time estimation of the earth's rotation parameters to obtain real-time ERP products with higher accuracy. Summary of the invention

[0005] The present invention provides a method and system for determining the earth rotation parameters based on GNSS data stream. The ERP is determined by determining the transformation matrix between the earth reference system and the celestial reference system using GNSS technology. The reference point coordinates in the earth reference system are determined by the ground station, while the satellite coordinates in the celestial reference system are determined by the dynamic orbit of the satellite. Therefore, the solution of station-satellite pair networking can realize the synchronous solution of parameters such as satellite orbit, clock error, ambiguity, troposphere and ERP. In a first aspect, the present invention provides a method for determining earth rotation parameters based on GNSS data stream, comprising: Step 1: Obtain the initial state of each station satellite and the parameters to be estimated based on the collected GNSS real-time observation data stream; Step 2: Get the initial value of the orbital position and perform numerical integration on the orbit to obtain the state transfer parameters of the previous and next epochs; Step 3: By updating the state of each parameter in various state transition models, an updated information matrix is ​​obtained; Step 4: Construct the observation equation based on the initial states of the satellites at each station and the parameters to be estimated; Step 5: Use the preset model to correct the observation value error in the observation equation and obtain the error equation after standardized residual; Step 6: Fuse the state quantity prior information and observation information, measure and update the information matrix at the next moment; Step 7: After the filter converges, the non-difference ambiguity fixing method is used to quickly fix the ambiguity parameters, and the ambiguity fixed value is added to the observation equation solved by the filter for constraint. Otherwise, skip this step and execute step 8; Step 8: Solve the information equation to obtain the corrected values ​​of all parameters at the next moment, and superimpose them with the initial values ​​after time update to obtain the estimated values ​​of all parameters at the current moment; 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.

[0006] According to a method for determining earth rotation parameters based on GNSS data stream provided by the present invention, step 1 comprises: The coordinates of each station in the earth reference system are regarded as constants. The implementation method is to constrain them to the specified GNSS data format solution that is updated regularly, and obtain the parameters to be estimated, including the earth's rotation parameters, satellite position, velocity, light pressure model parameters, satellite clock error, station clock error, ambiguity, station zenith tropospheric delay and inter-system bias; The Earth 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 change in day length, denoted as .

[0007] According to a method for determining earth rotation parameters based on GNSS data stream provided by the present invention, step 2 comprises: Determine the satellite's equation of motion in the inertial coordinate system:

[0008] in, is the position vector of the satellite's center of mass, is the satellite speed, is the dynamic parameter to be estimated in the satellite dynamics equation, , and They are the conservative force, non-conservative force and unmodeled empirical perturbation force acting on the satellite; The equations of motion are solved using numerical integration methods.

[0009] According to a method for determining earth rotation parameters based on GNSS data stream provided by the present invention, step 3 comprises: For the Earth's rotation parameters, the acceleration values ​​of polar motion and UT1 change are regarded as white noise, and the state transfer matrix of the previous and next epochs is derived based on the white noise:

[0010] in and represents two epochs, is the state transition matrix, Can characterize process noise The covariance matrix of , then:

[0011]

[0012] in represents the third-order identity matrix, is a 3rd order zero matrix, is the time interval between the previous and next epochs, and we can derive:

[0013] in is a 3rd order diagonal matrix, The diagonal elements are the variance values ​​of the noise per unit time process of the polar motion rate and day length change in the x and y directions, respectively.

[0014] According to a method for determining earth rotation parameters based on GNSS data stream provided by the present invention, step 4 comprises: Get the raw GNSS observations from the real-time streaming data and compose the ionosphere-free combined observations after gross error detection:

[0015] in, and are pseudorange and phase observations, respectively, is the geometric distance between the station and the satellite, and are the station and satellite clock errors, is the tropospheric delay, is the fuzziness, is the wavelength of the ionospheric-free phase combination, and is the unmodeled residual.

[0016] According to a method for determining earth rotation parameters based on GNSS data stream provided by the present invention, step 5 comprises: For the Earth's rotation parameters, we have:

[0017] in is the rotation matrix, , , and They represent precession, nutation, Earth rotation and polar motion, respectively. is the site coordinate in the earth reference system, which is considered as a constant here. is the satellite coordinate in the celestial coordinate system. Linearizing the above formula yields:

[0018] in , ,and They are , and The initial value of ERP is discussed only. , then:

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

[0020] in and They are the ERP parameters and The initial value of is the Earth's rotation angle, which can be calculated as follows:

[0021] in is the Julian day corresponding to UT1 time; For the rate of pole motion and the change in day length, add the derivative with respect to time ; Based on the UT1 forecast value provided by the IGS service organization or the IERS service organization, the virtual observation equation added for the UT1 value is:

[0022] in is the UT1 estimate obtained by time update, The UT1 value representing the external constraint adopts the ERP forecast product provided by IGS or IERS.

[0023] In a second aspect, the present invention further provides a system for determining earth rotation parameters based on GNSS data stream, comprising: The data acquisition module is used to obtain the initial state of each station satellite and the parameters to be estimated 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 numerically integrate the orbit to obtain the state transfer parameters of the previous and next epochs; The time update module is used to update the state of each parameter in various state transfer models to obtain the updated information matrix; A construction module is used to construct the observation equation based on the initial states of the satellites at each station and the parameters to be estimated; A correction module is used to correct the observation value error in the observation equation using a preset model to obtain an error equation after standardized residual; The fusion module is used to fuse the state quantity prior information and the observation information, and measure and update the information matrix at the next moment; The constraint module is used to use the non-difference ambiguity fixing method to quickly fix the ambiguity parameters after the filter converges, and add the ambiguity fixed value to the observation equation solved by the filter 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 correction values ​​of all parameters at the next moment, and superimpose them with the initial values ​​after time update to obtain the estimated values ​​of all parameters at the current moment; 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 repeatedly executed to the corresponding execution steps of the judgment module. If yes, the process ends.

[0024] In a third aspect, the present invention also provides an electronic device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the program, the method for determining the earth's rotation parameters based on the GNSS data stream as described in any one of the above is implemented.

[0025] In a fourth aspect, the present invention further provides a non-transitory computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements any of the above-described methods for determining the earth's rotation parameters based on the GNSS data stream.

[0026] In a fifth aspect, the present invention further provides a computer program product, comprising a computer program, which, when executed by a processor, implements any of the above-described methods for determining the earth's rotation parameters based on the GNSS data stream.

[0027] The method and system for determining the earth rotation parameters based on GNSS data stream provided by the present invention realizes real-time estimation and update of the earth rotation parameters through the GNSS real-time observation data stream, and can obtain high-precision ERP products 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. The present invention effectively solves the pain point that high precision and low latency cannot be achieved at the same time in the traditional earth rotation parameter determination method, and can better serve the application needs of high-precision real-time ERP. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] In order to more clearly illustrate the technical solutions in the present invention or the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0029] Figure 1 It is one of the flow charts of the method for determining the earth rotation parameters based on the GNSS data stream provided by the present invention; Figure 2 This is the second flow chart of the method for determining the earth rotation parameters based on the GNSS data stream provided by the present invention; Figure 3 It is a structural schematic diagram of the earth rotation parameter determination system based on GNSS data stream provided by the present invention; Figure 4 It is a structural schematic diagram of the electronic device provided by the present invention. DETAILED DESCRIPTION

[0030] In order to make the purpose, technical solution and advantages of the present invention clearer, the technical solution of the present invention will be clearly and completely described below in conjunction with the drawings of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0031] In view of the limitations in the prior art, the present invention uses GNSS technology to determine ERP by determining the transformation matrix between the earth reference system and the celestial reference system. In the process of implementing precise orbit determination using station-satellite networking, the synchronous solution of ERP can be achieved. The technical solution of the present invention is further explained below in conjunction with the accompanying drawings and specific embodiments.

[0032] Figure 1 is one of the flow charts of the method for determining the earth rotation parameters based on the GNSS data stream provided by an embodiment of the present invention, such as Figure 1 As shown, including: Step 1: Obtain the initial state of each station satellite and the parameters to be estimated based on the collected GNSS real-time observation data stream; Step 2: Get the initial value of the orbital position and numerically integrate the orbit to obtain the state transfer parameters of the previous and next epochs; Step 3: By updating the state of each parameter in various state transition models, an updated information matrix is ​​obtained; Step 4: Construct the observation equation based on the initial states of the satellites at each station and the parameters to be estimated; Step 5: Use the preset model to correct the observation value error in the observation equation and obtain the error equation after standardized residual; Step 6: Fuse the state quantity prior information and observation information, measure and update the information matrix at the next moment; Step 7: After the filter converges, the non-difference ambiguity fixing method is used to quickly fix the ambiguity parameters, and the ambiguity fixed value is added to the observation equation solved by the filter for constraint. Otherwise, skip this step and execute step 8; Step 8: Solve the information equation to obtain the corrected values ​​of all parameters at the next moment, and superimpose them with the initial values ​​after time update to obtain the estimated values ​​of all parameters at the current moment; 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.

[0033] Specifically, combined Figure 2 The logical route shown includes: Step 1, initialization, to prepare and preprocess GNSS data, including decoding and transferring of real-time streams, preparation of table files that various error models depend on, and data reading. In this way, the initial state of each station satellite and the parameters to be estimated are obtained to provide initial values ​​for filter estimation. The coordinates of the station in the earth reference system are regarded as constants, and are strongly constrained to the regularly updated SINEX solution or the post-precision single-point positioning static solution. The parameters to be estimated include the earth's rotation parameters, satellite position, speed, light pressure model parameters, satellite clock error, station clock error, ambiguity, station zenith tropospheric delay, and inter-system bias when GNSS multi-system joint processing is performed. The earth's rotation parameters include polar shift (x and y directions), polar shift rate (x and y directions), the difference between UT1 and UTC, and the change in day length, a total of 6 parameters, denoted as .

[0034] Step 2, orbital integration, obtain the initial value of the orbital position and the state transfer parameters of the previous and next epochs. The satellite's motion equation in the inertial coordinate system can be described by a set of first-order differential equations:

[0035] in, is the position vector of the satellite's center of mass, is the satellite speed, It is the dynamic parameter to be estimated in the satellite dynamics equation. , and They are the conservative force, non-conservative force and unmodeled empirical perturbation force acting on the satellite. Usually, the numerical integration method can be used to solve the differential equation, usually using a combination of the RKF single-step method and the Adams multi-step method.

[0036] Step 3, time update, for the earth's rotation parameters, the acceleration values ​​of polar shift and UT1 change are regarded as white noise, and the state transfer matrix of the previous and next epochs is derived. That is, the earth's rotation parameters of the previous and next epochs satisfy the following formula:

[0037] in and represents two epochs, is the state transition matrix, Can characterize process noise The covariance matrix of , then:

[0038]

[0039] in represents the third-order identity matrix, is a 3rd order zero matrix, is the time interval between the previous and next epochs. Therefore, it can be deduced that:

[0040] in It is a 3rd order diagonal matrix, and the diagonal elements are the variance values ​​of the polar motion rate (x and y directions) and the noise of the unit time process of the day length change. By performing time series analysis on the ERP sequence, when the correlation time is 300 seconds, Set to 1e-6 , 1e-6 and 1e-6 .

[0041] Step 4: Initialize the observation equation, read the GNSS observations, and construct the observation equation. Get the original GNSS observations from the real-time streaming data, and form the ionosphere-free combined observations after gross error detection:

[0042] in, and are pseudorange and phase observations, respectively, is the geometric distance between the station and the satellite, and are the station and satellite clock errors, is the tropospheric delay, is the fuzziness, is the wavelength of the ionospheric-free phase combination, and is the unmodeled residual.

[0043] Step 5: Form the error equation. Use common models to correct the observation error, and then complete the linearization of the observation equation for the parameters in step 1. For the Earth's rotation parameters, we have:

[0044] in is the rotation matrix, , , and They represent precession, nutation, Earth rotation and polar motion, respectively. is the site coordinate in the earth reference system, which is considered as a constant here. is the satellite coordinate in the celestial coordinate system. Linearizing formula (7) yields:

[0045] in , ,and They are , and Here we only discuss the ERP parameters. , then:

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

[0047] in and They are the ERP parameters and The initial value of is the Earth's rotation angle, which can be calculated as follows:

[0048] in is the Julian day corresponding to UT1 time. For the polar motion rate and day length change, we only need to add the differential of time That's it.

[0049] In addition, for the UT1 value in the earth rotation parameter, since it is coupled with the orbital parameters in GNSS orbit determination and cannot be estimated, the UT1 forecast 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:

[0050] in is the UT1 estimate obtained by time update, The UT1 value representing the external constraint can be the ERP forecast product provided by IGS or IERS.

[0051] Step 6: Non-differenced ambiguity fixation. After the filter converges, the non-differenced ambiguity fixation is achieved using the solved wide and narrow lane uncalibrated phase delay (UPD) information. The fixed ambiguity is added as constraint information to the observation equations to finally obtain the ambiguity fixed solution. Otherwise, skip this step.

[0052] Step 7: Measurement update. The state quantity prior information is integrated with the observation information, and the measurement update is obtained. The information matrix of the moment.

[0053] Step 8: Parameter recovery. Solve the information equation in step 7 to obtain The corrected values ​​of all parameters at the moment are added to the initial values ​​after the time update to obtain the estimated values ​​of all parameters at that moment.

[0054] 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, then end normally.

[0055] In one embodiment, a total of 32 days of GPS data from 120 globally distributed stations with real-time data streams from 309 to 340 in 2023 are used to conduct an experiment on real-time ERP estimation. The ERP parameters are updated in the real-time filtering orbit determination network solution, and the data processing interval is 300s. In order to ensure the accuracy of ambiguity fixation, non-difference ambiguity fixation is performed 2 days after the filter is started. After the ambiguity is fixed, a more accurate ERP valuation can be obtained, and the time interval of the ERP valuation sequence is 300s. Due to the time consumption of data decoding and data solution, the time delay of ERP valuation is about 60s.

[0056] The IERS 20C04 post-precision product is used as a reference value to evaluate the accuracy of the above ERP solution sequence. The IERS20C04 product is currently recognized as the highest accuracy ERP product, with a time delay of about 1 month, providing ERP values ​​at 0:00 and 12:00 within a day, while the sampling interval of the ERP real-time estimation product sequence is 300s. In order to avoid the accuracy loss of the reference product due to interpolation, the real-time estimated ERP solution sequence is thinned to the IERS 20C04 product time for accuracy evaluation. For comparative analysis, the accuracy of the forecast products (IERS eopc04_extended) available in real time from IGS and IERS in the same period is also evaluated in the same way. The results are shown in Table 1.

[0057] Table 1 Accuracy (RMS) statistics of IERS and IGU forecast products and ERP product series calculated by real-time filtering

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

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

[0060] Table 2 Statistical accuracy (RMS) of 30-day GPS real-time filtered precise orbit determination products under 3 ERP processing strategies (compared with the final precise products of IGS)

[0061] It can be seen that the results of using the IGU predicted ERP product directly for orbit determination in Scheme 1 are significantly worse, with an RMS value of 7.9 cm in the three-dimensional direction. Comparing Scheme 2 with Scheme 3, the method of this patent can obtain orbit results with higher accuracy, and the three-dimensional orbit accuracy is improved by more than 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 filtering precision orbit determination.

[0062] The earth rotation parameter determination system based on GNSS data stream provided by the present invention is described below. The earth rotation parameter determination system based on GNSS data stream described below and the earth rotation parameter determination method based on GNSS data stream described above can refer to each other.

[0063] Figure 3 is a schematic diagram of the structure of the earth rotation parameter determination system based on GNSS data stream provided by an embodiment of the present invention, such as Figure 3 As shown, it includes: a data acquisition module 31, a track 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: The data acquisition module 31 is used to obtain the initial state of each station satellite and the parameters to be estimated according to the collected GNSS real-time observation data stream; The orbit integration module 32 is used to obtain the initial value of the orbit position and numerically integrate the orbit to obtain the state transfer parameters of the previous and next epochs; The time updating module 33 is used to update the state of each parameter in various state transition models to obtain an updated information matrix; The construction module 34 is used to construct the observation equation based on the initial states of the satellites at each measuring station and the parameters to be estimated; The correction module 35 is used to correct the observation value error in the observation equation using a preset model to obtain an error equation after standardized residual; The fusion module 36 is used to fuse the state quantity prior information and the observation information, and measure and update the information matrix at the next moment; The constraint module 37 is used to use the non-difference ambiguity fixing method to quickly fix the ambiguity parameters after the filter converges, and add the ambiguity fixed value to the observation equation solved by the filter for constraint, otherwise skip this step and execute the corresponding execution steps of the solution module; The solving module 38 is used to solve the information equation to obtain the corrected values ​​of all parameters at the next moment, and to superimpose them with the initial values ​​after time update to obtain the estimated values ​​of all parameters at the current moment; The judgment module 39 is used to judge whether all epochs have been processed. If not, the initial state of the parameters is updated, and the orbit integration module is repeatedly executed to the corresponding execution steps of the judgment module. If yes, the process ends.

[0064] Figure 4 An example of a physical structure diagram of an electronic device is shown in FIG. Figure 4As shown, the electronic device may include: a processor (processor) 410, a communication interface (Communications Interface) 420, a memory (memory) 430 and a communication bus 440, wherein the processor 410, the communication interface 420, and the memory 430 communicate with each other through the communication bus 440. The processor 410 can call the logic instructions in the memory 430 to execute the method for determining the earth rotation parameters based on the GNSS data stream, and the method includes: Step 1: Obtain the initial state of each station satellite and the parameter to be estimated according to the collected GNSS real-time observation data stream; Step 2: Obtain the initial value of the orbital position and the numerical integration of the orbit to obtain the state transfer parameters of the previous and next epochs; Step 3: Obtain the updated information matrix by updating the state of each parameter in various state transfer models; Step 4: Construct an observation equation based on the initial state of each station satellite and the parameter to be estimated; Step 5: Use the preset model to correct the observation value error in the observation equation to obtain the standard The error equation after the residual is normalized; Step 6: Fuse the prior information of the state quantity with the observation information, and measure and update the information matrix at the next moment; Step 7: After the filter converges, the non-difference ambiguity fixing method is used to realize the rapid fixation of the ambiguity parameters, and the ambiguity fixed value is added to the observation equation solved by the filter for constraint, otherwise skip this step and execute step 8; Step 8: Solve the information equation to obtain the correction value of all parameters at the next moment, and superimpose it with the initial value after the time update to obtain the estimation value of all parameters at the current moment; 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.

[0065] In addition, the logic instructions in the above-mentioned memory 430 can be implemented in the form of a software functional unit and can be stored in a computer-readable storage medium when it is sold or used as an independent product. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art or the part of the technical solution, can be embodied in the form of a software product, and the computer software product is stored in a storage medium, including a number of instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to perform all or part of the steps of the method described in each embodiment of the present invention. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), disk or optical disk, etc. Various media that can store program codes.

[0066] The device embodiments described above are merely illustrative, wherein the units described as separate components may or may not be physically separated, and the components displayed as units may or may not be physical units, that is, they may be located in one place, or they may be distributed on multiple network units. Some or all of the modules may be selected according to actual needs to achieve the purpose of the scheme of this embodiment. Ordinary technicians in this field can understand and implement it without paying creative labor.

[0067] Through the description of the above implementation methods, those skilled in the art can clearly understand that each implementation method can be implemented by means of software plus a necessary general hardware platform, and of course, can also be implemented by hardware. Based on this understanding, the above technical solution is essentially or the part that contributes to the prior art can be embodied in the form of a software product, and the computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, a disk, an optical disk, etc., including a number of instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to execute the methods described in each embodiment or some parts of the embodiments.

[0068] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for determining earth rotation parameters based on GNSS data stream, characterized in that: include: Step 1: Obtain the initial state of each station satellite and the parameters to be estimated based on the collected GNSS real-time observation data stream; Step 2: Get the initial value of the orbital position and numerically integrate the orbit to obtain the state transfer parameters of the previous and next epochs; Step 3: By updating the state of each parameter in various state transition models, an updated information matrix is ​​obtained; Step 4: Construct the observation equation based on the initial states of the satellites at each station and the parameters to be estimated; Step 5: Use the preset model to correct the observation value error in the observation equation and obtain the error equation after standardized residual; Step 6: Fuse the state quantity prior information and observation information, measure and update the information matrix at the next moment; Step 7: After the filter converges, the non-difference ambiguity fixing method is used to quickly fix the ambiguity parameters, and the ambiguity fixed value is added to the observation equation solved by the filter for constraint. Otherwise, skip this step and execute step 8; Step 8: Solve the information equation to obtain the corrected values ​​of all parameters at the next moment, and superimpose them with the initial values ​​after time update to obtain the estimated values ​​of all parameters at the current moment; 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 the 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 implementation method is to constrain them to the specified GNSS data format solution that is updated regularly, and obtain the parameters to be estimated, including the earth's rotation parameters, satellite position, velocity, light pressure model parameters, satellite clock error, station clock error, ambiguity, station zenith tropospheric delay and inter-system bias; The Earth 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 change in day length, denoted as .

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

4. The method for determining the earth rotation parameters based on GNSS data stream according to claim 1, characterized in that: Step 3 includes: For the Earth's rotation parameters, the acceleration values ​​of polar motion and UT1 change are regarded as white noise, and the state transfer matrix of the previous and next epochs is derived based on the white noise: in and represents two epochs, is the state transition matrix, Can characterize process noise The covariance matrix of , then: in represents the third-order identity matrix, is a 3rd order zero matrix, is the time interval between the previous and next epochs, and we can derive: in is a 3rd order diagonal matrix, The diagonal elements are the variance values ​​of the noise per unit time process of the polar motion rate and day length change in the x and y directions, respectively.

5. The method for determining the earth rotation parameters based on GNSS data stream according to claim 1, characterized in that: Step 4 includes: Get the raw GNSS observations from the real-time streaming data and compose the ionosphere-free combined observations after gross error detection: in, and are pseudorange and phase observations, respectively, is the geometric distance between the station and the satellite, and are the station and satellite clock errors, is the tropospheric delay, is the fuzziness, is the wavelength of the ionospheric-free phase combination, and is the unmodeled residual.

6. The method for determining the earth rotation parameters based on GNSS data stream according to claim 1, characterized in that: Step 5 includes: For the Earth's rotation parameters, we have: in is the rotation matrix, , , and They represent precession, nutation, Earth rotation and polar motion, respectively. is the site coordinate in the earth reference system, which is considered as a constant here. is the satellite coordinate in the celestial coordinate system. Linearizing the above formula yields: in , ,and They are , and The initial value of ERP is discussed only. , then: in Represents the rotation matrix The initial value of the partial derivative on the right side of the equation can be expressed as: in and They are the ERP parameters and The initial value of is the Earth's rotation angle, It can be calculated by the following formula: in is the Julian day corresponding to UT1 time; For the rate of pole motion and the change in day length, add the derivative with respect to time ; Based on the UT1 forecast value provided by the IGS service organization or the IERS service organization, the virtual observation equation added for the UT1 value is: in is the UT1 estimate obtained by time update, The UT1 value representing the external constraint adopts the ERP forecast product provided by IGS or IERS.

7. A system for determining earth rotation parameters based on GNSS data stream, characterized in that: include: The data acquisition module is used to obtain the initial state of each station satellite and the parameters to be estimated 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 numerically integrate the orbit to obtain the state transfer parameters of the previous and next epochs; The time update module is used to update the state of each parameter in various state transfer models to obtain the updated information matrix; A construction module is used to construct the observation equation based on the initial states of the satellites at each station and the parameters to be estimated; A correction module is used to correct the observation value error in the observation equation using a preset model to obtain an error equation after standardized residual; The fusion module is used to fuse the state quantity prior information and the observation information, and measure and update the information matrix at the next moment; The constraint module is used to use the non-difference ambiguity fixing method to quickly fix the ambiguity parameters after the filter converges, and add the ambiguity fixed value to the observation equation solved by the filter 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 correction values ​​of all parameters at the next moment, and superimpose them with the initial values ​​after time update to obtain the estimated values ​​of all parameters at the current moment; 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 repeatedly executed to the corresponding execution steps of the judgment module. If yes, the process ends.

8. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: When the processor executes the program, the method for determining the earth rotation parameters based on the GNSS data stream as described in any one of claims 1 to 6 is implemented.

9. A non-transitory computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the method for determining the earth rotation parameters based on the GNSS data stream as claimed in any one of claims 1 to 6 is implemented.

10. A computer program product, comprising a computer program, characterized in that When the computer program is executed by a processor, the method for determining the earth rotation parameters based on the GNSS data stream as claimed in any one of claims 1 to 6 is implemented.

Images