Moon satellite formation orbit determination method and system

Abstract

The invention provides a lunar satellite formation orbit determination method and system. The method comprises the following steps: S1, establishing a lunar satellite high-precision dynamic model; s2, establishing a lunar satellite laser ranging and inter-satellite ranging mathematical model; s3, measuring the distance, namely the absolute position, of the A-type satellite relative to the ground station through laser ranging, and measuring the relative position of the B-type satellite relative to the A-type satellite through inter-satellite ranging; s4, preprocessing the absolute position and relative position data of the satellite measured by laser ranging and inter-satellite ranging, and generating a format data file required by orbit calculation; and S5, resolving the position coordinates of each satellite in the monthly formation through the format data file. According to the method, laser ranging and inter-satellite ranging are combined for orbit determination, precise orbit determination of the lunar satellite formation is achieved, compared with a traditional orbit determination mode, the precision can be improved from the hundred-meter level to the meter level, the requirements for the arrangement number of ground stations and set distribution are reduced, and the measurement and control pressure of existing ground stations is relieved.

Description

Technical Field

[0001] The present invention relates to the fields of on-orbit spacecraft measurement and control technology and satellite orbit calculation technology, and in particular to a method and system for determining the orbit of a lunar satellite formation. Background Art

[0002] With the development of aerospace technology, human beings are exploring the moon faster and faster, and the demand for various space missions is increasing. However, ground stations cannot cover all areas of the lunar surface without interruption and blind spots, so my country has proposed to establish a lunar satellite constellation. At present, the traditional lunar satellite orbit determination technology is ground-based radio observation technology, including S / X ranging and velocity measurement and Very Long Baseline Interferometry (VLBI) measurement technology, which has been fully applied in lunar exploration projects at home and abroad. However, the accuracy of radio orbit determination technology is low, only from kilometers to hundreds of meters; and the number of satellites in orbit is increasing. Relying entirely on ground measurement and control stations will cause the ground station to be overburdened and make the satellite system very fragile. With the development of lunar exploration technology in the new era, it is very necessary to expand the means of precise orbit determination of Earth-Moon space probes. In summary, in view of the needs of establishing a lunar constellation and the limitations of ground measurement and control resources and measurement and control accuracy, new high-precision orbit determination means need to be introduced.

[0003] As an emerging method, inter-satellite ranging technology determines orbits by measuring the distance and relative position between satellites, reducing dependence on ground stations. However, single inter-satellite ranging can only determine the relative position between satellites, and is limited by error accumulation and overall rotation of the constellation, requiring other high-precision measurement methods for compensation and correction.

[0004] Lunar laser ranging is currently the most accurate method of measuring the distance between the Earth and the Moon. With its high precision and long distance, laser ranging has become an ideal supplement to the intersatellite ranging system.

[0005] The invention patent with the existing publication number CN114935947B discloses a method and electronic equipment for satellite formation holding control, including: calculating the relative orbital elements according to the absolute orbital elements of the current master and slave satellites; calculating the nominal relative orbital elements according to the decay rate of the relative orbital elements; calculating the current formation orbital element holding error according to the current relative orbital elements and the nominal relative orbital elements; calculating the ignition latitude angle and pulse size of the pulse with the goal of minimizing fuel consumption according to the current formation orbital element holding error; and calculating the continuous thrust ignition start and end time at each pulse using the propulsion system parameters based on the ignition latitude angle and pulse size of the pulse obtained by the calculation. Summary of the invention

[0006] In view of the defects in the prior art, the present invention provides a method and system for determining the orbit of a lunar satellite formation.

[0007] According to a method and system for determining the orbit of a lunar satellite formation provided by the present invention, the scheme is as follows: In a first aspect, a method for determining the orbit of a lunar satellite formation is provided, the method comprising: dividing the lunar orbit satellites in the lunar satellite formation into Class A satellites and Class B satellites; the Class A satellites are equipped with laser corner reflectors and cooperate with ground laser stations to carry out laser ranging tests; the Class A satellites are equipped with integrated communication devices and cooperate with Class B satellites to carry out inter-satellite measurement technology tests; the Class B satellites are equipped with communication and navigation payloads and cooperate with Class A satellites to carry out inter-satellite measurement technology tests; the specific method for determining the orbit comprises: Step S1: Establishing a high-precision dynamic model of the lunar satellite; Step S2: Establishing mathematical models of lunar satellite laser ranging and inter-satellite ranging; Step S3: measuring the distance, i.e., the absolute position, of the Class A satellite relative to the ground station by laser ranging, and measuring the relative position of the Class B satellite relative to the Class A satellite by inter-satellite ranging; Step S4: pre-processing the absolute position and relative position data of the satellite measured by laser ranging and inter-satellite ranging to generate a data file in a format required for orbit calculation; Step S5: Calculate the position coordinates of each satellite in the lunar formation using the formatted data file.

[0008] Preferably, the step S1 includes: for a lunar satellite, the coordinate system selects the epoch moon-center inertial coordinate system, and the high-precision dynamic model of the lunar satellite is as follows:

[0009] in, for The differential of is the state vector of the state model; They are the position and velocity of the lunar satellite in the X, Y, and Z directions respectively; is the system nonlinear continuous state transfer function of the state model; Represents system noise.

[0010] Preferably, the step S2 comprises: The lunar satellite laser ranging observation is the distance from the laser station to the satellite. Suppose the original observation distance is , the observation equation is:

[0011] In the formula, The ranging error caused by the change of station position due to the earth's tides; The distance measurement error caused by the refraction effect of light in the atmosphere; The distance measurement error caused by the relativistic effect of light in the gravitational field; is the deviation of the laser reflection point on the satellite surface from the center of mass; is the system delay error of the station; Intersatellite ranging is a two-way one-way distance measurement in a time-division system. The ranging value contains the distance and clock difference information of the two satellites. By summing and subtracting the two ranging values, the clock difference and distance information of the two satellites are decoupled. The observation equation for intersatellite two-way ranging is:

[0012] In the formula, for The relative distance between satellites A and B at time, i.e. the observed value of inter-satellite ranging; They are the transmission delay and reception delay of Class A satellite respectively; They are the transmission delay and reception delay of Class B satellite respectively; It is the error in one-way ranging, including the error caused by the phase center deviation of the satellite antenna and the relativistic effect.

[0013] Preferably, the step S4 pre-processes the absolute position and relative position data of the satellite by using observation equations, including: 1) Space-time conversion: To ensure the accuracy of error correction, all calculations must be performed in the solar system barycenter coordinate system; 2) System error correction: In laser ranging, the system error is the delay error of the laser ranging system of the ground station; in inter-satellite ranging, the system error is the transmission and reception delay of Class A satellites and Class B satellites; 3) Other error corrections: Tropospheric refraction correction: Based on the temperature, humidity and pressure data measured at the station, the tropospheric refraction correction is completed for satellite ranging and velocity measurement through the atmospheric refraction correction model; Relativistic delay correction: Based on the distance from the satellite to each celestial body, the bending of light caused by the gravitational field of the celestial body is calculated to complete the relativistic delay correction; Earth tide correction: Based on the solid tide, ocean tide, atmospheric load tide generated by the external gravitational force on the earth, and the solid polar tide and ocean polar tide parameters caused by the centrifugal disturbance caused by the earth's rotation, the coordinate deformation variables of the ground station are calculated to complete the earth tide correction; Center of mass offset correction: The center of mass offset correction is completed according to the distance from the on-board corner reflector to the satellite center of mass; Satellite antenna phase center correction: correct the satellite antenna phase offset; 4) Data format conversion: After the orbit measurement data is corrected and processed, the format conversion is performed to generate data content in the format required for orbit determination calculation.

[0014] Preferably, calculating the position coordinates of each satellite in the lunar formation in step S5 includes: Step S5.1: Calculate the initial orbit using orbit measurement data combined with satellite dynamics equations; Step S5.2: perform orbital integration using the initial orbital information to calculate the reference orbit and state transfer matrix; Step S5.3: Linearize the observation equation according to the reference orbit and the state transfer matrix; Step S5.4: Perform linear optimal estimation on the linearized observation equation to obtain solution parameters; Step S5.5: Calculate theoretical observation values ​​according to the solution parameters; Step S5.6: Calculate the residuals, remove the observations whose residuals exceed the threshold, and obtain the orbit determination results of each satellite.

[0015] Preferably, the step S5.4 comprises: The state differential equation of satellite motion is expressed as:

[0016] in,

[0017] In the formula, and Indicates the state quantity to be estimated; represents the kinetic parameters; Indicates the state quantity to be estimated at the initial moment of the orbit determination arc segment; Indicates the initial time of the orbit determination arc; Indicates the state quantity of the orbit determination arc at the initial moment; Indicates satellite position; Indicates the satellite speed; represents the satellite acceleration; represents other parameters to be estimated in the kinetic model; In the above formula In the reference state Expand and remember , after ignoring the higher-order terms, it can be expressed as a linear equation:

[0018] In the formula, let , ; Its solution is expressed as: ; Solved ; in, represents the state transfer matrix, represents the derivative of the state transfer matrix; represents the identity matrix; Indicates the reference state; Indicates the difference between the state to be estimated and the reference state; Indicates time; represents the difference derivative between the estimated state and the reference state; express The difference between the state to be estimated and the reference state at any moment; Satellite in Observable quantity at time , represents the measurement noise; Indicates that the satellite is The state vector at the moment; Indicates the observation time; Represents observation data The corresponding truth value; i Represents the i-th data; Expand the above formula at the reference state and consider only the first-order terms:

[0019] In the formula, represents the observed partial derivative of the observed quantity with respect to the state quantity at the observed epoch, represents the observed partial derivative of the improved epoch state quantity; represents the actual observed value; Get the linear equation ; in, Indicates the state quantity to be estimated; Represents random error.

[0020] Preferably, the step S5.5 comprises: Solution The best estimate of is the least square method for parameter estimation; the weight matrix of the observed quantity is , and the estimate is obtained based on the linear unbiased minimum variance estimate:

[0021] in, express The transpose of k Represents the amount of observation data.

[0022] In a second aspect, a system for determining the orbit of a lunar satellite formation is provided, the system comprising: dividing the lunar orbit satellites in the lunar satellite formation into Class A satellites and Class B satellites; the Class A satellites are equipped with laser corner reflectors to cooperate with ground laser stations to carry out laser ranging tests; the Class A satellites are equipped with integrated communication machines to cooperate with Class B satellites to carry out inter-satellite measurement technology tests; the Class B satellites are equipped with communication and navigation payloads to cooperate with Class A satellites to carry out inter-satellite measurement technology tests; the orbit determination specific system comprises: Module M1: Establish a high-precision dynamic model of the lunar satellite; Module M2: Establish mathematical models for lunar satellite laser ranging and inter-satellite ranging; Module M3: Measure the distance, i.e., absolute position, of a Class A satellite relative to a ground station through laser ranging, and measure the relative position of a Class B satellite relative to a Class A satellite through intersatellite ranging; Module M4: Preprocess the absolute and relative position data of the satellite measured by laser ranging and inter-satellite ranging to generate data files in the format required for orbit calculation; Module M5: Calculate the position coordinates of each satellite in the lunar formation through formatted data files.

[0023] Preferably, the module M1 includes: for a lunar satellite, the coordinate system selects the epoch moon-center inertial coordinate system, and the high-precision dynamic model of the lunar satellite is as follows:

[0024] in, for The differential of is the state vector of the state model; They are the position and velocity of the lunar satellite in the X, Y, and Z directions respectively; is the system nonlinear continuous state transfer function of the state model; represents the system noise; The module M2 comprises: The lunar satellite laser ranging observation is the distance from the laser station to the satellite. Suppose the original observation distance is , the observation equation is:

[0025] In the formula, The ranging error caused by the change of station position due to the earth's tides; The distance measurement error caused by the refraction effect of light in the atmosphere; The distance measurement error caused by the relativistic effect of light in the gravitational field; is the deviation of the laser reflection point on the satellite surface from the center of mass; is the system delay error of the station; Intersatellite ranging is a two-way one-way distance measurement in a time-division system. The ranging value contains the distance and clock difference information of the two satellites. By summing and subtracting the two ranging values, the clock difference and distance information of the two satellites are decoupled. The observation equation for intersatellite two-way ranging is:

[0026] In the formula, for The relative distance between satellites A and B at time, i.e. the observed value of inter-satellite ranging; They are the transmission delay and reception delay of Class A satellite respectively; They are the transmission delay and reception delay of Class B satellite respectively; It is the error in one-way ranging, including the error caused by the phase center deviation of the satellite antenna and the relativistic effect; The module M4 pre-processes the absolute position and relative position data of the satellite through observation equations, including: 1) Space-time conversion: To ensure the accuracy of error correction, all calculations must be performed in the solar system barycenter coordinate system; 2) System error correction: In laser ranging, the system error is the delay error of the laser ranging system of the ground station; in inter-satellite ranging, the system error is the transmission and reception delay of Class A satellites and Class B satellites; 3) Other error corrections: Tropospheric refraction correction: Based on the temperature, humidity and pressure data measured at the station, the tropospheric refraction correction is completed for satellite ranging and velocity measurement through the atmospheric refraction correction model; Relativistic delay correction: Based on the distance from the satellite to each celestial body, the bending of light caused by the gravitational field of the celestial body is calculated to complete the relativistic delay correction; Earth tide correction: Based on the solid tide, ocean tide, atmospheric load tide generated by the external gravitational force on the earth, and the solid polar tide and ocean polar tide parameters caused by the centrifugal disturbance caused by the earth's rotation, the coordinate deformation variables of the ground station are calculated to complete the earth tide correction; Center of mass offset correction: The center of mass offset correction is completed according to the distance from the on-board corner reflector to the satellite center of mass; Satellite antenna phase center correction: correct the satellite antenna phase offset; 4) Data format conversion: After the orbit measurement data is corrected and processed, the format conversion is performed to generate data content in the format required for orbit determination calculation.

[0027] Preferably, the module M5 calculates the position coordinates of each satellite in the lunar formation, including: Module M5.1: Calculate the initial orbit using orbit measurement data combined with satellite dynamics equations; Module M5.2: Use the initial orbital information to perform orbital integration and calculate the reference orbit and state transfer matrix; Module M5.3: Linearize the observation equations based on the reference orbit and state transfer matrix; Module M5.4: Perform linear optimal estimation on the linearized observation equation to obtain the solution parameters; Module M5.5: Calculate theoretical observations based on solution parameters; Module M5.6: Calculate the residuals, remove the observations whose residuals exceed the threshold, and obtain the orbit determination results for each satellite; The module M5.4 includes: The state differential equation of satellite motion is expressed as:

[0028] in,

[0029] In the formula, and Indicates the state quantity to be estimated; represents the kinetic parameters; Indicates the state quantity to be estimated at the initial moment of the orbit determination arc segment; Indicates the initial time of the orbit determination arc; Indicates the state quantity of the orbit determination arc at the initial moment; Indicates satellite position; Indicates the satellite speed; represents the satellite acceleration; represents other parameters to be estimated in the kinetic model; In the above formula In the reference state Expand and remember , after ignoring the higher-order terms, it can be expressed as a linear equation:

[0030] In the formula, ; Its solution is expressed as: ; Solved ; in, represents the state transfer matrix, represents the derivative of the state transfer matrix; represents the identity matrix; Indicates the reference state; Indicates the difference between the state to be estimated and the reference state; Indicates time; represents the difference derivative between the estimated state and the reference state; express The difference between the state to be estimated and the reference state at any moment; Satellite in Observable quantity at time , represents the measurement noise; Indicates that the satellite is The state vector at the moment; Indicates the observation time; Represents observation data The corresponding truth value; i Represents the i-th data; Expand the above formula at the reference state and consider only the first-order terms:

[0031] In the formula, represents the observed partial derivative of the observed quantity with respect to the state quantity at the observed epoch, represents the observed partial derivative of the improved epoch state quantity; represents the actual observed value; Get the linear equation ; in, Indicates the state quantity to be estimated; represents random error; The module M5.5 includes: Solution The best estimate of is the least square method for parameter estimation; the weight matrix of the observed quantity is , and the estimate is obtained based on the linear unbiased minimum variance estimate:

[0032] in, express The transpose of k Represents the amount of observation data.

[0033] Compared with the prior art, the present invention has the following beneficial effects: 1. The present invention uses laser ranging and intersatellite ranging to jointly determine orbits, achieving precise orbit determination of the lunar satellite formation. Compared with traditional orbit determination methods, the accuracy can be improved from hundreds of meters to meters; 2. The present invention utilizes inter-satellite ranging information to reduce the number and distribution requirements of ground stations, thereby alleviating the measurement and control pressure of existing ground stations; 3. The present invention introduces laser ranging information, which effectively solves the problem of overall rotation of the constellation.

[0034] Other beneficial effects of the present invention will be explained in the specific implementation manner through the introduction of specific technical features and technical solutions. Through the introduction of these technical features and technical solutions, those skilled in the art should be able to understand the beneficial technical effects brought about by the technical features and technical solutions. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Other features, objects and advantages of the present invention will become more apparent from the detailed description of non-limiting embodiments made with reference to the following drawings: Figure 1 It is the overall flow chart of the present invention; Figure 2 This is a measurement model diagram of the system of the present invention; Figure 3 This is a data preprocessing flow chart of the present invention; Figure 4 This is a flow chart of the track improvement process of the present invention. DETAILED DESCRIPTION

[0036] The present invention is described in detail below in conjunction with specific embodiments. The following embodiments will help those skilled in the art to further understand the present invention, but are not intended to limit the present invention in any form. It should be noted that, for those of ordinary skill in the art, several changes and improvements can also be made without departing from the concept of the present invention. These all belong to the protection scope of the present invention.

[0037] The embodiment of the present invention provides a method for determining the orbit of a lunar satellite formation, which conducts joint observations of laser ranging and inter-satellite ranging to achieve centimeter-level precise orbit determination of a lunar satellite formation, and can greatly reduce the number of ground tracking stations and lower the requirements for the geometric distribution of ground stations, thus becoming a means that can be considered to solve the current problems of low orbit determination accuracy of lunar orbit satellites and limited resources of ground measurement and control stations.

[0038] Specifically, the lunar orbiting satellites in the lunar satellite formation are divided into Class A satellites and Class B satellites; Class A satellites are equipped with laser corner reflectors to cooperate with ground laser stations to carry out laser ranging tests; Class A satellites are equipped with integrated communication equipment to cooperate with Class B satellites to carry out inter-satellite measurement technology tests; Class B satellites are equipped with communication and navigation payloads to cooperate with Class A satellites to carry out inter-satellite measurement technology tests.

[0039] The ground laser station adjusts the telescope pointing according to the satellite orbit prediction information, transmits laser to the corner reflector on the A-type satellite, receives the laser reflected by the corner reflector, and records the time of laser transmission and reception. The position and velocity of the main celestial bodies in the solar system are obtained according to the DE430 ephemeris.

[0040] Reference Figure 1 As shown, the specific method for determining the orbit includes: Step S1: Establishing a high-precision dynamic model of the lunar satellite; High-precision dynamic models include: lunar non-spherical perturbation model, third-body gravity model caused by the sun and the earth, lunar solid tide perturbation model, lunar physical libration model, solar radiation pressure model, jet unloading, etc.

[0041] The non-spherical perturbation of the moon is the effect of the change of the gravitational field caused by the uneven distribution of the moon's mass on the orbiter's motion; The third-body gravitational perturbation caused by the sun and the earth is the additional disturbance effect of the sun's and the earth's gravity on the lunar orbiter. The lunar solid tidal perturbation is the influence of the periodic change of the gravitational field caused by the tidal deformation of the moon due to the gravitational effect of the earth on the orbiter. The lunar physical libration perturbation is the disturbance of the orbiter caused by the change of the gravitational field caused by the libration during the rotation of the moon. Solar light pressure perturbation is the motion disturbance caused by the tiny pressure exerted by solar radiation photons on the orbiter surface. Jet unloading is the instantaneous effect of the reaction force on the trajectory of the orbiter when it uses the jet device to unload the flywheel angular momentum; For the lunar satellite, the coordinate system is the lunar center inertial coordinate system. The high-precision dynamic model of the lunar satellite is as follows:

[0042] in, for The differential of is the state vector of the state model; They are the position and velocity of the lunar satellite in the X, Y, and Z directions respectively; is the system nonlinear continuous state transfer function of the state model; Represents system noise.

[0043] Step S2: Establishing mathematical models of lunar satellite laser ranging and inter-satellite ranging; Perform data preprocessing on the laser ranging distance information to generate data files in the format required for orbit calculation; Data preprocessing includes: system delay correction, tropospheric refraction correction, relativistic delay correction, earth tide correction, corner reflector center of mass offset correction, time and space correction, and data format conversion; the data file format used for orbit calculation is CRD (Consolidated Prediction Format).

[0044] The lunar satellite laser ranging observation is the distance from the laser station to the satellite. Suppose the original observation distance is , the observation equation is:

[0045] In the formula, The ranging error caused by the change of station position due to the earth's tides; The distance measurement error caused by the refraction effect of light in the atmosphere; The distance measurement error caused by the relativistic effect of light in the gravitational field; is the deviation of the laser reflection point on the satellite surface from the center of mass; is the system delay error of the station; Intersatellite ranging is a two-way one-way distance measurement in a time-division system. The ranging value contains the distance and clock difference information of the two satellites. By summing and subtracting the two ranging values, the clock difference and distance information of the two satellites are decoupled. The observation equation for intersatellite two-way ranging is:

[0046] In the formula, for The relative distance between satellites A and B at time, i.e. the observed value of inter-satellite ranging; They are the transmission delay and reception delay of Class A satellite respectively; They are the transmission delay and reception delay of Class B satellite respectively; It is the error in one-way ranging, including the error caused by the phase center deviation of the satellite antenna and the relativistic effect.

[0047] Step S3: measuring the distance, i.e., the absolute position, of the Class A satellite relative to the ground station by laser ranging, and measuring the relative position of the Class B satellite relative to the Class A satellite by inter-satellite ranging; The N satellites in the satellite formation are respectively combined with satellites in the satellite formation except the current satellite to obtain N×(N-1) satellite combinations; each satellite combination includes a current satellite and a link-building satellite, and the N satellites in the satellite formation are all used as current satellites, and one current satellite corresponds to N-1 satellite combinations; N>1.

[0048] Step S4: pre-processing the absolute position and relative position data of the satellite measured by laser ranging and inter-satellite ranging to generate a data file in a format required for orbit calculation; The satellite's absolute and relative position data are preprocessed through observation equations, including: 1) Space-time conversion: To ensure the accuracy of error correction, all calculations must be performed in the solar system barycenter coordinate system; 2) System error correction: In laser ranging, the system error is the delay error of the laser ranging system of the ground station; in inter-satellite ranging, the system error is the transmission and reception delay of Class A satellites and Class B satellites; 3) Other error corrections: Tropospheric refraction correction: Based on the temperature, humidity and pressure data measured at the station, the tropospheric refraction correction is completed for satellite ranging and velocity measurement through the atmospheric refraction correction model; Relativistic delay correction: Based on the distance from the satellite to each celestial body, the bending of light caused by the gravitational field of the celestial body is calculated to complete the relativistic delay correction; Earth tide correction: Based on the solid tide, ocean tide, atmospheric load tide generated by the external gravitational force on the earth, and the solid polar tide and ocean polar tide parameters caused by the centrifugal disturbance caused by the earth's rotation, the coordinate deformation variables of the ground station are calculated to complete the earth tide correction; Center of mass offset correction: The center of mass offset correction is completed according to the distance from the on-board corner reflector to the satellite center of mass; Satellite antenna phase center correction: correct the satellite antenna phase offset; 4) Data format conversion: After the orbit measurement data is corrected and processed, the format conversion is performed to generate a data file in the format required for orbit determination calculation.

[0049] Step S5: Calculate the position coordinates of each satellite in the lunar formation using the formatted data file.

[0050] The step S5 calculates the position coordinates of each satellite in the lunar formation, including: Step S5.1: Calculate the initial orbit using orbit measurement data combined with satellite dynamics equations; Step S5.2: perform orbital integration using the initial orbital information to calculate the reference orbit and state transfer matrix; Step S5.3: Linearize the observation equation according to the reference orbit and the state transfer matrix; Step S5.4: Perform linear optimal estimation on the linearized observation equation to obtain solution parameters; The state differential equation of satellite motion is expressed as:

[0051] in,

[0052] In the formula, and Indicates the state quantity to be estimated; represents the kinetic parameters; Indicates the state quantity to be estimated at the initial moment of the orbit determination arc segment; Indicates the initial time of the orbit determination arc; Indicates the state quantity of the orbit determination arc at the initial moment; Indicates the satellite position; Indicates the satellite speed; represents the satellite acceleration; represents other parameters to be estimated in the kinetic model; In the above formula In the reference state Expand and remember , after ignoring the higher-order terms, it can be expressed as a linear equation:

[0053] In the formula, ; Its solution is expressed as: ; Solved ; in, represents the state transfer matrix, represents the derivative of the state transfer matrix; represents the identity matrix; Indicates the reference state; Indicates the difference between the state to be estimated and the reference state; Indicates time; represents the difference derivative between the estimated state and the reference state; express The difference between the state to be estimated and the reference state at any moment; Satellite in Observable quantity at time , represents the measurement noise; Indicates that the satellite is The state vector at the moment; Indicates the observation time; Represents observation data The corresponding truth value; i Represents the i-th data; Expand the above formula at the reference state and consider only the first-order terms:

[0054] In the formula, represents the observed partial derivative of the observed quantity with respect to the state quantity at the observed epoch, represents the observed partial derivative of the improved epoch state quantity; represents the actual observed value; Get the linear equation ; in, Indicates the state quantity to be estimated; Represents random error.

[0055] Step S5.5: Calculate theoretical observation values ​​according to the solution parameters; This step S5.5 includes: Solution The best estimate of is the least square method for parameter estimation; the weight matrix of the observed quantity is , and the estimate is obtained based on the linear unbiased minimum variance estimate:

[0056] in, express The transpose of k Represents the amount of observation data.

[0057] Step S5.6: Calculate the residuals, remove the observations whose residuals exceed the threshold, and obtain the orbit determination results of each satellite.

[0058] In the present invention, for Class A satellites: Using orbit measurement data and satellite dynamics equations, the initial satellite orbit is calculated based on the Laplace method. In the orbit determination arc, the observation values ​​are used to determine the orbit and the orbit parameters of each satellite are obtained by calculation; Calculate theoretical observation values ​​of the observation values ​​according to the solution parameters; The residuals are calculated based on the observed values ​​and the theoretical observed values, and the observed values ​​whose residuals exceed a threshold are eliminated.

[0059] In the orbit determination arc, the observation values ​​are used to determine the orbit and obtain the solution parameters, including: Using the initial orbit information to perform orbital integration, a reference orbit and state transfer matrix sampled at certain time intervals are obtained; Linearizing the observation equation according to the reference orbit and the state transfer matrix; The linear optimal estimation method is used to solve the observation equation after linearization to obtain the solution parameters.

[0060] The observation values ​​include laser station data; and the theoretical observation values ​​include theoretical observation values ​​of the laser station data.

[0061] For Class B satellites: According to the initial orientation elements and intersatellite measurement data provided by the A-type satellite, combined with the satellite dynamics equation, the satellite initial orbit is calculated based on the Laplace method; In the orbit determination arc, the observation values ​​are used to determine the orbit and the orbit parameters of each satellite are obtained by calculation; Calculate theoretical observation values ​​of the observation values ​​according to the solution parameters; The residuals are calculated based on the observed values ​​and the theoretical observed values, and the observed values ​​whose residuals exceed a threshold are eliminated.

[0062] In the orbit determination arc, the observation values ​​are used to determine the orbit and obtain the solution parameters, including: Using the initial orbit information to perform orbital integration, a reference orbit and state transfer matrix sampled at certain time intervals are obtained; The observation equation is linearized according to the reference orbit and the state transfer matrix; The linear optimal estimation method is used to solve the observation equation after linearization to obtain the solution parameters.

[0063] The observed values ​​include inter-satellite ranging data; and the theoretical observed values ​​include theoretical observed values ​​of the inter-satellite ranging data.

[0064] Each satellite in the satellite formation is used as the local satellite for measurement update and time update to complete the orbit determination of each satellite.

[0065] The present invention also provides a system for determining the orbit of a lunar satellite formation. The system can be implemented by executing the process steps of the method for determining the orbit of a lunar satellite formation, that is, a person skilled in the art can understand the method for determining the orbit of a lunar satellite formation as a preferred implementation of the system for determining the orbit of a lunar satellite formation. The system includes: dividing the lunar orbit satellites in the lunar satellite formation into Class A satellites and Class B satellites; Class A satellites are equipped with laser corner reflectors to carry out laser ranging tests in cooperation with ground laser stations; Class A satellites are equipped with integrated communication machines to carry out inter-satellite measurement technology tests in cooperation with Class B satellites; Class B satellites are equipped with communication and navigation payloads to carry out inter-satellite measurement technology tests in cooperation with Class A satellites; the specific system for determining the orbit includes: Module M1: Establish a high-precision dynamic model of the lunar satellite; Specifically, module M1 includes: for the lunar satellite, the coordinate system selects the lunar center inertial coordinate system of the epoch, and the high-precision dynamic model of the lunar satellite is as follows:

[0066] in, for The differential of is the state vector of the state model; They are the position and velocity of the lunar satellite in the X, Y, and Z directions respectively; is the system nonlinear continuous state transfer function of the state model; represents the system noise; Module M2: Establish mathematical models for lunar satellite laser ranging and inter-satellite ranging; Module M2 includes: The lunar satellite laser ranging observation is the distance from the laser station to the satellite. Suppose the original observation distance is , the observation equation is:

[0067] In the formula, The ranging error caused by the change of station position due to the earth's tides; The distance measurement error caused by the refraction effect of light in the atmosphere; The distance measurement error caused by the relativistic effect of light in the gravitational field; is the deviation of the laser reflection point on the satellite surface from the center of mass; is the system delay error of the station; Intersatellite ranging is a two-way one-way distance measurement in a time-division system. The ranging value contains the distance and clock difference information of the two satellites. By summing and subtracting the two ranging values, the clock difference and distance information of the two satellites are decoupled. The observation equation for intersatellite two-way ranging is:

[0068] In the formula, for The relative distance between satellites A and B at time, i.e. the observed value of inter-satellite ranging; They are the transmission delay and reception delay of Class A satellite respectively; They are the transmission delay and reception delay of Class B satellite respectively; It is the error in one-way ranging, including the error caused by the phase center deviation of the satellite antenna and the relativistic effect.

[0069] Module M3: Measure the distance, i.e., absolute position, of a Class A satellite relative to a ground station through laser ranging, and measure the relative position of a Class B satellite relative to a Class A satellite through intersatellite ranging; Module M4: Preprocess the absolute and relative position data of the satellite measured by laser ranging and inter-satellite ranging to generate data files in the format required for orbit calculation; Module M4 pre-processes the absolute and relative position data of the satellite through observation equations, including: 1) Space-time conversion: To ensure the accuracy of error correction, all calculations must be performed in the solar system barycenter coordinate system; 2) System error correction: In laser ranging, the system error is the delay error of the laser ranging system of the ground station; in inter-satellite ranging, the system error is the transmission and reception delay of Class A satellites and Class B satellites; 3) Other error corrections: Tropospheric refraction correction: Based on the temperature, humidity and pressure data measured at the station, the tropospheric refraction correction is completed for satellite ranging and velocity measurement through the atmospheric refraction correction model; Relativistic delay correction: Based on the distance from the satellite to each celestial body, the bending of light caused by the gravitational field of the celestial body is calculated to complete the relativistic delay correction; Earth tide correction: Based on the solid tide, ocean tide, atmospheric load tide generated by the external gravitational force on the earth, and the solid polar tide and ocean polar tide parameters caused by the centrifugal disturbance caused by the earth's rotation, the coordinate deformation variables of the ground station are calculated to complete the earth tide correction; Center of mass offset correction: The center of mass offset correction is completed according to the distance from the on-board corner reflector to the satellite center of mass; Satellite antenna phase center correction: correct the satellite antenna phase offset; 4) Data format conversion: After the orbit measurement data is corrected and processed, the format conversion is performed to generate data content in the format required for orbit determination calculation.

[0070] Module M5: Calculate the position coordinates of each satellite in the lunar formation through formatted data files.

[0071] Module M5 calculates the position coordinates of each satellite in the lunar formation, including: Module M5.1: Calculate the initial orbit using orbit measurement data combined with satellite dynamics equations; Module M5.2: Use the initial orbital information to perform orbital integration and calculate the reference orbit and state transfer matrix; Module M5.3: Linearize the observation equations based on the reference orbit and state transfer matrix; Module M5.4: Perform linear optimal estimation on the linearized observation equation to obtain the solution parameters; The module M5.4 includes: The state differential equation of satellite motion is expressed as:

[0072] in,

[0073] In the formula, and Indicates the state quantity to be estimated; represents the kinetic parameters; Indicates the state quantity to be estimated at the initial moment of the orbit determination arc segment; Indicates the initial time of the orbit determination arc; Indicates the state quantity of the orbit determination arc at the initial moment; Indicates the satellite position; Indicates the satellite speed; represents the satellite acceleration; represents other parameters to be estimated in the kinetic model; In the above formula In the reference state Expand and remember , after ignoring the higher-order terms, it can be expressed as a linear equation:

[0074] In the formula, ; Its solution is expressed as: ; Solved ; in, represents the state transfer matrix, represents the derivative of the state transfer matrix; represents the identity matrix; Indicates the reference state; Indicates the difference between the state to be estimated and the reference state; Indicates time; represents the difference derivative between the state to be estimated and the reference state; express The difference between the state to be estimated and the reference state at any moment; Satellite in Observable quantity at time , represents the measurement noise; Indicates that the satellite is The state vector at the moment; Indicates the observation time; Represents observation data The corresponding truth value; i Represents the i-th data; Expand the above formula at the reference state and consider only the first-order terms:

[0075] In the formula, represents the observed partial derivative of the observed quantity with respect to the state quantity at the observed epoch, represents the observed partial derivative of the improved epoch state quantity; represents the actual observed value; Get the linear equation ; in, Indicates the state quantity to be estimated; represents random error; Module M5.5: Calculate theoretical observations based on solution parameters; The M5.5 module includes: Solution The best estimate of is the least square method for parameter estimation; the weight matrix of the observed quantity is , and the estimate is obtained based on the linear unbiased minimum variance estimate:

[0076] in, express The transpose of k Represents the amount of observation data.

[0077] Module M5.6: Calculate the residuals, remove the observations whose residuals exceed the threshold, and obtain the orbit determination results for each satellite.

[0078] Next, the present invention will be described in more detail.

[0079] The present invention provides a method for determining the orbit of a lunar satellite formation, which specifically includes: Precise orbit determination is the process of fitting satellite orbits using precise dynamics using various types of observation data containing observation errors. Obviously, the accuracy of orbit determination calculations is largely limited by the accuracy of measurement data, the accuracy of observation models, and the accuracy of dynamic models, in addition to observation geometry. Precise orbit determination involves the processing of various types of observation data, including precise observation modeling and correction of various measurement errors. Satellites of different orbit types are subject to different dynamic influences, so the dynamic models considered in orbit determination calculations are also different. In addition, both the dynamic model and the observation model have certain errors, especially the dynamic model part. In summary, in the process of orbit calculation, it is necessary to select a suitable mathematical model and optimize the measurement data. This will improve the accuracy of the calculation results.

[0080] Figure 1 This is a flow chart of an embodiment of a method for combined orbit determination of lunar satellite formation laser ranging and inter-satellite ranging according to the present invention. Figure 2 This is a system measurement model diagram related to the joint orbit determination method of the lunar satellite formation laser ranging and inter-satellite ranging of the present invention. Figure 1 and Figure 2 The laser ranging and intersatellite ranging combined orbit determination method comprises: Step 1: Establish a high-precision dynamic model of the lunar satellite.

[0081] For the lunar satellite, the coordinate system is the lunar center inertial coordinate system of epoch (J2000.0), considering the gravitational forces of the sun, moon, earth and other celestial bodies on the lunar satellite, the lunar solid tidal perturbation, the lunar physical libration, solar radiation pressure, jet unloading and high-precision ephemeris, the high-precision dynamic model of the lunar satellite is as follows:

[0082] The above formula can be simplified as:

[0083] in, for The differential of is the state vector of the state model; are the position and velocity of the lunar satellite in the X, Y, and Z directions respectively; is the system nonlinear continuous state transfer function of the state model; represents the system noise; represents the vector from the center of the moon to the satellite; Represents the vector from the center of the Earth to the satellite; heliocentric-to-satellite vector; are the gravitational constants of the Sun, Moon, and Earth respectively; is the vector from the heliocenter to the satellite; is the vector from the center of the moon to the center of the earth in the geocentric coordinate system; is the vector from the center of the moon to the center of the sun; is the coordinate of the moon's position in the solar mass center coordinate system; is the coordinate of the earth's position in the geocentric coordinate system; is the coordinate of the sun's position in the sun's mass center coordinate system; are the system noise respectively.

[0084] Step 2: Establish mathematical models for lunar satellite laser ranging and intersatellite ranging.

[0085] The lunar satellite laser ranging observation is the distance from the laser station to the satellite. Suppose the original observation distance is , the observation equation is:

[0086] In the formula, The ranging error caused by the change of station position due to the earth's tides; The distance measurement error caused by the refraction effect of light in the atmosphere; The distance measurement error caused by the relativistic effect of light in the gravitational field; is the deviation of the laser reflection point on the satellite surface from the center of mass; is the system delay error of the station; Intersatellite ranging is a two-way one-way distance measurement in a time-division system. The ranging value contains the distance and clock difference information of the two satellites. By summing and subtracting the two ranging values, the clock difference and distance information of the two satellites are decoupled. The observation equation for intersatellite two-way ranging is:

[0087] In the formula, for The relative distance between satellites A and B at time, i.e. the observed value of inter-satellite ranging; They are the transmission delay and reception delay of Class A satellite respectively; They are the transmission delay and reception delay of Class B satellite respectively; It is the error in one-way ranging, including the error caused by the phase center deviation of the satellite antenna and the relativistic effect.

[0088] Step 3: The ground laser station conducts a laser ranging test to obtain the distance information of Class A satellite relative to the ground station, which is the absolute position information; and uses inter-satellite two-way ranging to provide the relative position information of Class B satellite relative to Class A satellite.

[0089] Step 4: Use the observation equation in step 2 to preprocess the original position information obtained in step 3: The specific steps of data preprocessing are as follows: Figure 3 As shown, including: (1) Space-time conversion: In order to ensure the accuracy of error correction, all calculations must be performed in the solar system barycenter coordinate system; (2) System error correction: In laser ranging, the system error is the delay error of the laser ranging system of the ground station; in intersatellite ranging, the system error is the transmission and reception delay of the A-type satellite and the B-type satellite; (3) Error correction: 1. Tropospheric refraction correction; Based on the temperature, humidity and pressure data measured at the measuring station, the tropospheric refraction correction is completed for the satellite ranging and velocity measurement through the atmospheric refraction correction model.

[0090] 2. Relativistic delay correction; Based on the distance from the satellite to each celestial body, the bending of light caused by the gravitational field of the celestial body is calculated to complete the relativistic time delay correction.

[0091] 3. Earth tide correction; Based on the solid tide, ocean tide, atmospheric load tide generated by the external gravitational force on the earth, and the parameters of solid polar tide and ocean polar tide caused by the centrifugal disturbance caused by the rotation of the earth, the coordinate deformation variables of the ground station are calculated to complete the earth tide correction.

[0092] 4. Correction of center of mass offset; The center of mass offset correction is completed based on the distance from the on-board corner reflector to the satellite's center of mass.

[0093] 5. Satellite antenna phase center correction; Correction for satellite antenna phase offset.

[0094] (4) Data format conversion: After the orbit measurement data is corrected and processed, the format conversion is performed to generate data content in the format required for orbit determination calculation.

[0095] Step 5: Calculate the position coordinates of each satellite in the formation, such as Figure 4As shown, the following steps are included: (1) Initial orbit calculation: using orbit measurement data to calculate and generate the satellite's initial orbit information for orbit improvement: Based on the Laplace method, the ranging data is used to complete the satellite's initial orbit calculation. The dynamic model used in the initial orbit calculation is the model in step 1, which is suitable for satellites in lunar orbit.

[0096] (2) Use the initial orbital information to perform orbital integration and calculate the reference orbit and state transfer matrix; Calculating Satellite Orbits Using Numerical Integration Methods

[0097] in, Respectively expressed in and Satellite position vector at this moment; Respectively expressed in and Satellite velocity vector at this moment; represents the time step; represents the higher-order error term; Using the orbital integration method, the initial orbital parameters are taken as initial values, and the trajectory obtained by integration calculation is called the reference orbit; The state differential equation of satellite motion is expressed as:

[0098] in, ; In the formula, and Indicates the state quantity to be estimated; represents the kinetic parameters; Indicates the state quantity to be estimated at the initial moment of the orbit determination arc segment; Indicates the initial time of the orbit determination arc; Indicates the state quantity of the orbit determination arc at the initial moment; Indicates the satellite position; Indicates the satellite speed; represents the satellite acceleration; represents other parameters to be estimated in the kinetic model; (3) Linearize the observation equation according to the reference orbit and state transfer matrix; In the above formula In the reference state Expand and remember , after ignoring the higher-order terms, it can be expressed as a linear equation:

[0099] In the formula, ; Its solution is expressed as: ; The solution of the state equation is ; in, represents the state transfer matrix, represents the derivative of the state transfer matrix; represents the identity matrix; Indicates the reference state; Indicates the difference between the state to be estimated and the reference state; Indicates time; represents the difference derivative between the state to be estimated and the reference state; express The difference between the state to be estimated and the reference state at any moment; (4) Perform linear optimal estimation on the linearized observation equation to obtain the solution parameters; Satellite in Observable quantity at time , represents the measurement noise; Indicates that the satellite is The state vector at the moment; Indicates the observation time; Represents observation data The corresponding truth value; i Represents the i-th data; Expand the above formula at the reference state and consider only the first-order terms:

[0100] In the formula, represents the observed partial derivative of the observed quantity with respect to the state quantity at the observed epoch, represents the observed partial derivative of the improved epoch state quantity; represents the actual observed value; Get the linear equation ; in, Indicates the state quantity to be estimated; Represents random error.

[0101] Step 5: Calculate theoretical observation values ​​based on the solution parameters.

[0102] Solution The best estimate of is the least square method for parameter estimation; the weight matrix of the observed quantity is , according to the linear unbiased minimum variance estimation, we can get the estimate of the theoretical observation value:

[0103] in, express The transpose of k Represents the amount of observation data.

[0104] Step 6: Calculate the residual based on the difference between the observed value and the theoretical value, and eliminate the observations whose residual exceeds the threshold.

[0105] Both the state equation and the observation equation are the results of linear approximation. The errors caused by the nonlinear part are solved through continuous iteration, and the observations whose residuals exceed the threshold are eliminated.

[0106] Each satellite in the satellite formation is used as the local satellite for measurement update and time update to obtain the orbit determination result of each satellite.

[0107] The embodiments of the present invention provide a method and system for determining the orbit of a lunar satellite formation, which uses laser ranging and inter-satellite ranging to jointly determine the orbit, thereby realizing precise orbit determination of a lunar satellite formation. Compared with traditional orbit determination methods, the accuracy can be improved from hundreds of meters to meters, and the requirements for the number and distribution of ground stations are reduced, thereby alleviating the measurement and control pressure of existing ground stations.

[0108] Those skilled in the art know that, in addition to realizing the system and its various devices, modules, and units provided by the present invention in a purely computer-readable program code, it is entirely possible to realize the same functions in the form of logic gates, switches, application-specific integrated circuits, programmable logic controllers, and embedded microcontrollers by logically programming the method steps. Therefore, the system and its various devices, modules, and units provided by the present invention can be considered as a hardware component, and the devices, modules, and units included therein for realizing various functions can also be regarded as structures within the hardware component; the devices, modules, and units for realizing various functions can also be regarded as both software modules for realizing the method and structures within the hardware component.

[0109] The above describes the specific embodiments of the present invention. It should be understood that the present invention is not limited to the above specific embodiments, and those skilled in the art can make various changes or modifications within the scope of the claims, which does not affect the essence of the present invention. In the absence of conflict, the embodiments of the present application and the features in the embodiments can be combined with each other arbitrarily.

Claims

1. A method for determining the orbit of a lunar satellite formation, characterized in that: The lunar orbit satellites in the lunar satellite formation are divided into Class A satellites and Class B satellites; the Class A satellites are equipped with laser corner reflectors to cooperate with ground laser stations to carry out laser ranging tests; the Class A satellites are equipped with integrated communication devices to cooperate with Class B satellites to carry out inter-satellite measurement technology tests; the Class B satellites are equipped with communication and navigation payloads to cooperate with Class A satellites to carry out inter-satellite measurement technology tests; the specific method for determining the orbit includes: Step S1: Establishing a high-precision dynamic model of the lunar satellite; Step S2: Establishing mathematical models of lunar satellite laser ranging and inter-satellite ranging; Step S3: measuring the distance, i.e., the absolute position, of the Class A satellite relative to the ground station by laser ranging, and measuring the relative position of the Class B satellite relative to the Class A satellite by inter-satellite ranging; Step S4: pre-processing the absolute position and relative position data of the satellite measured by laser ranging and inter-satellite ranging to generate a data file in a format required for orbit calculation; Step S5: Calculate the position coordinates of each satellite in the lunar formation using the formatted data file.

2. The method for determining the orbit of a lunar satellite formation according to claim 1, characterized in that: The step S1 includes: for a lunar satellite, the coordinate system selects the lunar center inertial coordinate system of the epoch, and the high-precision dynamic model of the lunar satellite is as follows: in, for The differential of is the state vector of the state model; They are the position and velocity of the lunar satellite in the X, Y, and Z directions respectively; is the system nonlinear continuous state transfer function of the state model; Represents system noise.

3. The method for determining the orbit of a lunar satellite formation according to claim 1, characterized in that: The step S2 comprises: The lunar satellite laser ranging observation is the distance from the laser station to the satellite. Suppose the original observation distance is , the observation equation is: In the formula, The ranging error caused by the change of station position due to the earth's tides; The distance measurement error caused by the refraction effect of light in the atmosphere; The distance measurement error caused by the relativistic effect of light in the gravitational field; is the deviation of the laser reflection point on the satellite surface from the center of mass; is the system delay error of the station; Intersatellite ranging is a two-way one-way distance measurement in a time-division system. The ranging value contains the distance and clock difference information of the two satellites. By summing and subtracting the two ranging values, the clock difference and distance information of the two satellites are decoupled. The observation equation for intersatellite two-way ranging is: In the formula, for The relative distance between satellites A and B at time, i.e. the observed value of inter-satellite ranging; They are the transmission delay and reception delay of Class A satellite respectively; They are the transmission delay and reception delay of Class B satellite respectively; It is the error in one-way ranging, including the error caused by the phase center deviation of the satellite antenna and the relativistic effect.

4. The method for determining the orbit of a lunar satellite formation according to claim 3, characterized in that: The step S4 pre-processes the absolute position and relative position data of the satellite through the observation equation, including: 1) Space-time conversion: To ensure the accuracy of error correction, all calculations must be performed in the solar system barycenter coordinate system; 2) System error correction: In laser ranging, the system error is the delay error of the laser ranging system of the ground station; in inter-satellite ranging, the system error is the transmission and reception delay of Class A satellites and Class B satellites; 3) Other error corrections: Tropospheric refraction correction: Based on the temperature, humidity and pressure data measured at the station, the tropospheric refraction correction is completed for satellite ranging and velocity measurement through the atmospheric refraction correction model; Relativistic delay correction: Based on the distance from the satellite to each celestial body, the bending of light caused by the gravitational field of the celestial body is calculated to complete the relativistic delay correction; Earth tide correction: Based on the solid tide, ocean tide, atmospheric load tide generated by the external gravitational force on the earth, and the solid polar tide and ocean polar tide parameters caused by the centrifugal disturbance caused by the earth's rotation, the coordinate deformation variables of the ground station are calculated to complete the earth tide correction; Center of mass offset correction: The center of mass offset correction is completed according to the distance from the on-board corner reflector to the satellite center of mass; Satellite antenna phase center correction: correct the satellite antenna phase offset; 4) Data format conversion: After the orbit measurement data is corrected and processed, the format conversion is performed to generate a data file in the format required for orbit determination calculation.

5. The method for determining the orbit of a lunar satellite formation according to claim 1, characterized in that: The step S5 of calculating the position coordinates of each satellite in the lunar formation includes: Step S5.1: Calculate the initial orbit using orbit measurement data combined with satellite dynamics equations; Step S5.2: perform orbital integration using the initial orbital information to calculate the reference orbit and state transfer matrix; Step S5.3: Linearize the observation equation according to the reference orbit and the state transfer matrix; Step S5.4: Perform linear optimal estimation on the linearized observation equation to obtain solution parameters; Step S5.5: Calculate theoretical observation values ​​according to the solution parameters; Step S5.6: Calculate the residuals, remove the observations whose residuals exceed the threshold, and obtain the orbit determination results of each satellite.

6. The method for determining the orbit of a lunar satellite formation according to claim 5, characterized in that: The step S5.4 comprises: The state differential equation of satellite motion is expressed as: in, ; In the formula, and Indicates the state quantity to be estimated; represents the kinetic parameters; Indicates the state quantity to be estimated at the initial moment of the orbit determination arc segment; Indicates the initial time of the orbit determination arc; Indicates the state quantity of the orbit determination arc at the initial moment; Indicates satellite position; Indicates the satellite speed; represents the satellite acceleration; represents other parameters to be estimated in the kinetic model; In the above formula In the reference state Expand and remember , after ignoring the higher-order terms, it can be expressed as a linear equation: In the formula, ; Its solution is expressed as: ; Solved ; in, represents the state transfer matrix, represents the derivative of the state transfer matrix; represents the identity matrix; Indicates the reference state; Indicates the difference between the state to be estimated and the reference state; Indicates time; represents the difference derivative between the estimated state and the reference state; express The difference between the state to be estimated and the reference state at any moment; Satellite in Observable quantity at time , represents the measurement noise; Indicates that the satellite is The state vector at the moment; Indicates the observation time; Represents observation data The corresponding truth value; i Represents the i-th data; Expand the above formula at the reference state and consider only the first-order terms: In the formula, represents the observed partial derivative of the observed quantity with respect to the state quantity at the observed epoch, represents the observed partial derivative of the improved epoch state quantity; represents the actual observed value; Get the linear equation ; in, Indicates the state quantity to be estimated; Represents random error.

7. The method for determining the orbit of a lunar satellite formation according to claim 5, characterized in that: The step S5.5 comprises: Solution The best estimate of is the least square method for parameter estimation; the weight matrix of the observed quantity is , and the estimate is obtained based on the linear unbiased minimum variance estimate: in, express The transpose of k Represents the amount of observation data.

8. A lunar satellite formation orbit determination system, characterized in that: include: The lunar orbit satellites in the lunar satellite formation are divided into Class A satellites and Class B satellites; the Class A satellites are equipped with laser corner reflectors to cooperate with ground laser stations to carry out laser ranging tests; the Class A satellites are equipped with integrated communication devices to cooperate with Class B satellites to carry out inter-satellite measurement technology tests; the Class B satellites are equipped with communication and navigation payloads to cooperate with Class A satellites to carry out inter-satellite measurement technology tests; the orbit determination specific system includes: Module M1: Establish a high-precision dynamic model of the lunar satellite; Module M2: Establish mathematical models for lunar satellite laser ranging and inter-satellite ranging; Module M3: Measure the distance, i.e., absolute position, of a Class A satellite relative to a ground station through laser ranging, and measure the relative position of a Class B satellite relative to a Class A satellite through intersatellite ranging; Module M4: Preprocess the absolute and relative position data of the satellite measured by laser ranging and inter-satellite ranging to generate data files in the format required for orbit calculation; Module M5: Calculate the position coordinates of each satellite in the lunar formation through formatted data files.

9. The orbit determination system for a lunar satellite formation according to claim 8, characterized in that: The module M1 includes: for the lunar satellite, the coordinate system selects the lunar center inertial coordinate system of the epoch, and the high-precision dynamic model of the lunar satellite is as follows: in, for The differential of is the state vector of the state model; They are the position and velocity of the lunar satellite in the X, Y, and Z directions respectively; is the system nonlinear continuous state transfer function of the state model; represents the system noise; The module M2 comprises: The lunar satellite laser ranging observation is the distance from the laser station to the satellite. Suppose the original observation distance is , the observation equation is: In the formula, The ranging error caused by the change of station position due to the earth's tides; The distance measurement error caused by the refraction effect of light in the atmosphere; The distance measurement error caused by the relativistic effect of light in the gravitational field; is the deviation of the laser reflection point on the satellite surface from the center of mass; is the system delay error of the station; Intersatellite ranging is a two-way one-way distance measurement in a time-division system. The ranging value contains the distance and clock difference information of the two satellites. By summing and subtracting the two ranging values, the clock difference and distance information of the two satellites are decoupled. The observation equation for intersatellite two-way ranging is: In the formula, for The relative distance between satellites A and B at time, i.e. the observed value of inter-satellite ranging; They are the transmission delay and reception delay of Class A satellite respectively; They are the transmission delay and reception delay of Class B satellite respectively; It is the error in one-way ranging, including the error caused by the phase center deviation of the satellite antenna and the relativistic effect; The module M4 pre-processes the absolute position and relative position data of the satellite through observation equations, including: 1) Space-time conversion: To ensure the accuracy of error correction, all calculations must be performed in the solar system barycenter coordinate system; 2) System error correction: In laser ranging, the system error is the delay error of the laser ranging system of the ground station; in inter-satellite ranging, the system error is the transmission and reception delay of Class A satellites and Class B satellites; 3) Other error corrections: Tropospheric refraction correction: Based on the temperature, humidity and pressure data measured at the station, the tropospheric refraction correction is completed for satellite ranging and velocity measurement through the atmospheric refraction correction model; Relativistic delay correction: Based on the distance from the satellite to each celestial body, the bending of light caused by the gravitational field of the celestial body is calculated to complete the relativistic delay correction; Earth tide correction: Based on the solid tide, ocean tide, atmospheric load tide generated by the external gravitational force on the earth, and the solid polar tide and ocean polar tide parameters caused by the centrifugal disturbance caused by the earth's rotation, the coordinate deformation variables of the ground station are calculated to complete the earth tide correction; Center of mass offset correction: The center of mass offset correction is completed according to the distance from the on-board corner reflector to the satellite center of mass; Satellite antenna phase center correction: correct the satellite antenna phase offset; 4) Data format conversion: After the orbit measurement data is corrected and processed, the format conversion is performed to generate a data file in the format required for orbit determination calculation.

10. The orbit determination system for a lunar satellite formation according to claim 8, characterized in that: The module M5 calculates the position coordinates of each satellite in the lunar formation, including: Module M5.1: Calculate the initial orbit using orbit measurement data combined with satellite dynamics equations; Module M5.2: Use the initial orbital information to perform orbital integration and calculate the reference orbit and state transfer matrix; Module M5.3: Linearize the observation equations based on the reference orbit and state transfer matrix; Module M5.4: Perform linear optimal estimation on the linearized observation equation to obtain the solution parameters; Module M5.5: Calculate theoretical observations based on solution parameters; Module M5.6: Calculate the residuals, remove the observations whose residuals exceed the threshold, and obtain the orbit determination results for each satellite; The module M5.4 includes: The state differential equation of satellite motion is expressed as: in, In the formula, and Indicates the state quantity to be estimated; represents the kinetic parameters; Indicates the state quantity to be estimated at the initial moment of the orbit determination arc segment; Indicates the initial time of the orbit determination arc; Indicates the state quantity of the orbit determination arc at the initial moment; Indicates satellite position; Indicates the satellite speed; represents the satellite acceleration; represents other parameters to be estimated in the kinetic model; In the above formula In the reference state Expand and remember , after ignoring the higher-order terms, it can be expressed as a linear equation: In the formula, ; Its solution is expressed as: ; Solved ; in, represents the state transfer matrix, represents the derivative of the state transfer matrix; represents the identity matrix; Indicates the reference state; Indicates the difference between the state to be estimated and the reference state; Indicates time; represents the difference derivative between the estimated state and the reference state; express The difference between the state to be estimated and the reference state at any moment; Satellite in Observable quantity at time , represents the measurement noise; Indicates that the satellite is The state vector at the moment; Indicates the observation time; Represents observation data The corresponding truth value; i Represents the i-th data; Expand the above formula at the reference state and consider only the first-order terms: In the formula, represents the observed partial derivative of the observed quantity with respect to the state quantity at the observed epoch, represents the observed partial derivative of the improved epoch state quantity; represents the actual observed value; Get the linear equation ; in, Indicates the state quantity to be estimated; represents random error; The module M5.5 includes: Solution The best estimate of is the least square method for parameter estimation; the weight matrix of the observed quantity is , and the estimate is obtained based on the linear unbiased minimum variance estimate: in, express The transpose of k Represents the amount of observation data.

Images