Method and system for determining a lunar satellite formation orbit

Abstract

The application provides a method and system for determining the formation orbit of a lunar satellite, comprising the following steps: S1, establishing a high-precision dynamics model of the lunar satellite; S2, establishing a laser ranging and inter-satellite ranging mathematical model of the lunar satellite; S3, measuring the distance of a satellite of type A relative to a ground station, i.e. the absolute position, by laser ranging, and measuring the relative position of a satellite of type B relative to the satellite of type A by inter-satellite ranging; S4, preprocessing the absolute position and relative position data of the satellites measured by laser ranging and inter-satellite ranging to generate a data file in a format required for orbit calculation; and S5, calculating the position coordinates of each satellite in the lunar satellite formation by using the data file. The application uses laser ranging and inter-satellite ranging to jointly determine the orbit, realizes precise orbit determination of the lunar satellite formation, and compared with the traditional orbit determination mode, the precision can be improved from the hundred-meter level to the meter level, and the requirements for the number and distribution of the ground stations are reduced, and the pressure of the existing ground stations is reduced.

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 advancement of space technology, humanity's lunar exploration efforts are accelerating, and the demand for various space missions is increasing. However, ground stations cannot provide uninterrupted, blind-spot coverage of every region of the lunar surface. Therefore, my country has proposed establishing a lunar satellite constellation. Currently, traditional lunar satellite orbit determination technologies rely on ground-based radio observations, including S / X ranging and velocity measurements and Very Long Baseline Interferometry (VLBI). These technologies have been extensively utilized in lunar exploration projects both domestically and internationally. However, radio orbit determination offers relatively low accuracy, ranging from kilometers to hundreds of meters. Furthermore, with the increasing number of satellites in orbit, relying solely on ground-based tracking and control (TT&C) stations would place an excessive burden on ground stations and render the satellite system extremely vulnerable. With the development of new-era lunar exploration technology, it is imperative to expand the means for precise orbit determination of cis-lunar space probes. In summary, to meet the demands of establishing a lunar constellation and the limitations of ground-based tracking and control (TT&C) resources and accuracy, new high-precision orbit determination methods are needed.

[0003] Intersatellite ranging technology, as an emerging method, determines orbits by measuring the distance and relative positions between satellites, reducing reliance on ground stations. However, intersatellite ranging alone can only determine the relative positions of satellites and is limited by error accumulation and overall constellation rotation, 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 existing invention patent with publication number CN114935947B discloses a method and electronic equipment for satellite formation holding control, including: calculating the relative orbital elements based on the absolute orbital elements of the current master and slave satellites; calculating the nominal relative orbital elements based on the decay rate of the relative orbital elements; calculating the current formation orbital element holding error based on the current relative orbital elements and the nominal relative orbital elements; calculating the ignition latitude argument and pulse size of the pulse based on the current formation orbital element holding error with the goal of minimizing fuel consumption; and calculating the continuous thrust ignition start and end time at each pulse using the propulsion system parameters based on the calculated ignition latitude argument and pulse size of the pulse. Summary of the Invention

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

[0007] According to the present invention, a method and system for determining the orbit of a lunar satellite formation is provided, and the scheme is as follows:

[0008] In a first aspect, a method for determining the orbit of a lunar satellite formation is provided. The method includes: 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 orbit determination method includes:

[0009] Step S1: Establish a high-precision dynamic model of the lunar satellite;

[0010] Step S2: Establishing mathematical models for lunar satellite laser ranging and inter-satellite ranging;

[0011] 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;

[0012] Step S4: Preprocessing the absolute and relative position data of the satellites measured by laser ranging and inter-satellite ranging to generate a data file in the format required for orbit calculation;

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

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

[0015]

[0016] 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 system noise.

[0017] Preferably, step S2 includes:

[0018] The lunar satellite laser ranging observation quantity is the distance from the laser station to the satellite. Let the original observation distance be , the observation equation is:

[0019]

[0020] Where, The ranging error caused by the change of station position due to the earth's tides; 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;

[0021] 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:

[0022]

[0023] Where, for The relative distance between satellites A and B at time instant, i.e. the observed value of inter-satellite ranging; are the transmission delay and reception delay of Class A satellite respectively; 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.

[0024] Preferably, the step S4 pre-processes the absolute position and relative position data of the satellite using observation equations, including:

[0025] 1) Space-time conversion: To ensure the accuracy of error correction, all calculations must be performed in the solar system barycenter coordinate system;

[0026] 2) System error correction: In laser ranging, the system error is the delay error of the ground station laser ranging system; in inter-satellite ranging, the system error is the transmission and reception delay of Class A and Class B satellites;

[0027] 3) Other error corrections:

[0028] Tropospheric refraction correction: Based on the temperature, humidity and pressure data measured at the observation station, the atmospheric refraction correction model is used to correct the tropospheric refraction of the satellite ranging and velocity measurements.

[0029] Relativistic delay correction: Based on the distance between the satellite and each celestial body, the bending of light caused by the gravitational field of the celestial body is calculated to complete the relativistic delay correction;

[0030] 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 extreme tide and ocean extreme 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;

[0031] Center of mass offset correction: The center of mass offset correction is completed based on the distance between the on-board corner reflector and the satellite center of mass;

[0032] Satellite antenna phase center correction: correct the satellite antenna phase offset;

[0033] 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.

[0034] Preferably, calculating the position coordinates of each satellite in the lunar formation in step S5 includes:

[0035] Step S5.1: Calculate the initial orbit using orbit measurement data combined with satellite dynamics equations;

[0036] Step S5.2: Using the initial orbital information to perform orbital integration, calculate the reference orbital and state transfer matrix;

[0037] Step S5.3: Linearize the observation equation according to the reference orbit and state transfer matrix;

[0038] Step S5.4: Perform linear optimal estimation on the linearized observation equation to obtain solution parameters;

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

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

[0041] Preferably, the step S5.4 includes:

[0042] The state differential equation of satellite motion is expressed as:

[0043]

[0044] in,

[0045] Where, and Represents 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 arc segment at the initial moment; Indicates satellite position; Indicates the satellite speed; represents the satellite acceleration; Represents other parameters to be estimated in the dynamic model;

[0046] In the above formula In the reference state Expand and remember , after omitting the higher-order terms, it can be expressed as a linear equation:

[0047]

[0048] In the formula, let , ;

[0049] Its solution is expressed as: ;

[0050] Solution ;

[0051] in, represents the state transition 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;

[0052] Satellite in Observable quantity at a moment , 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;

[0053] Expand the above formula at the reference state and consider only the first-order terms:

[0054]

[0055] Where, represents the observed partial derivative of the observed quantity with respect to the state quantity at the observation epoch, represents the observed partial derivative of the improved epoch state quantity; represents the actual observed value;

[0056] Get the linear equation ;

[0057] in, Represents the state quantity to be estimated; Represents random error.

[0058] Preferably, the step S5.5 includes:

[0059] Solution The best estimate of is obtained by using the least squares method for parameter estimation; the weight matrix of the observation is recorded as , and the estimated value is obtained based on the linear unbiased minimum variance estimate:

[0060]

[0061] in, express The transpose of k Indicates the amount of observation data.

[0062] In a second aspect, a system for determining the orbit of a lunar satellite formation is provided. The system includes: 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 orbit determination system specifically includes:

[0063] Module M1: Establish a high-precision dynamic model of the lunar satellite;

[0064] Module M2: Establish mathematical models for lunar satellite laser ranging and intersatellite ranging;

[0065] Module M3: measures the distance, i.e., the absolute position, of a Class A satellite relative to a ground station through laser ranging, and measures the relative position of a Class B satellite relative to a Class A satellite through inter-satellite ranging;

[0066] Module M4: Pre-processes the absolute and relative position data of satellites measured by laser ranging and inter-satellite ranging to generate data files in the format required for orbit calculation;

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

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

[0069]

[0070] 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;

[0071] The module M2 includes:

[0072] The lunar satellite laser ranging observation quantity is the distance from the laser station to the satellite. Let the original observation distance be , the observation equation is:

[0073]

[0074] Where, The ranging error caused by the change of station position due to the earth's tides; 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;

[0075] 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:

[0076]

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

[0078] The module M4 pre-processes the absolute and relative position data of the satellite through observation equations, including:

[0079] 1) Space-time conversion: To ensure the accuracy of error correction, all calculations must be performed in the solar system barycenter coordinate system;

[0080] 2) System error correction: In laser ranging, the system error is the delay error of the ground station laser ranging system; in inter-satellite ranging, the system error is the transmission and reception delay of Class A and Class B satellites;

[0081] 3) Other error corrections:

[0082] Tropospheric refraction correction: Based on the temperature, humidity and pressure data measured at the observation station, the atmospheric refraction correction model is used to correct the tropospheric refraction of the satellite ranging and velocity measurements.

[0083] Relativistic delay correction: Based on the distance between the satellite and each celestial body, the bending of light caused by the gravitational field of the celestial body is calculated to complete the relativistic delay correction;

[0084] 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 extreme tide and ocean extreme 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;

[0085] Center of mass offset correction: The center of mass offset correction is completed based on the distance between the on-board corner reflector and the satellite center of mass;

[0086] Satellite antenna phase center correction: correct the satellite antenna phase offset;

[0087] 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.

[0088] Preferably, the module M5 calculates the position coordinates of each satellite in the lunar formation, including:

[0089] Module M5.1: Calculate the initial orbit using orbital measurement data combined with satellite dynamics equations;

[0090] Module M5.2: Use the initial orbital information to perform orbital integration and calculate the reference orbit and state transfer matrix;

[0091] Module M5.3: Linearize the observation equations based on reference orbits and state transfer matrices;

[0092] Module M5.4: Perform linear optimal estimation on the linearized observation equation to obtain the solution parameters;

[0093] Module M5.5: Calculate theoretical observations based on solution parameters;

[0094] Module M5.6: Calculate the residuals, remove observations with residuals exceeding the threshold, and obtain the orbit determination results for each satellite;

[0095] The module M5.4 includes:

[0096] The state differential equation of satellite motion is expressed as:

[0097]

[0098] in,

[0099] Where, and Represents 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 arc segment at the initial moment; Indicates satellite position; Indicates the satellite speed; represents the satellite acceleration; Represents other parameters to be estimated in the dynamic model;

[0100] In the above formula In the reference state Expand and remember , after omitting the higher-order terms, it can be expressed as a linear equation:

[0101]

[0102] Where, ;

[0103] Its solution is expressed as: ;

[0104] Solution ;

[0105] in, represents the state transition 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;

[0106] Satellite in Observable quantity at a moment , 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;

[0107] Expand the above formula at the reference state and consider only the first-order terms:

[0108]

[0109] Where, represents the observed partial derivative of the observed quantity with respect to the state quantity at the observation epoch, represents the observed partial derivative of the improved epoch state quantity; represents the actual observed value;

[0110] Get the linear equation ;

[0111] in, Represents the state quantity to be estimated; represents random error;

[0112] The module M5.5 includes:

[0113] Solution The best estimate of is obtained by using the least squares method for parameter estimation; the weight matrix of the observation is recorded as , and the estimated value is obtained based on the linear unbiased minimum variance estimate:

[0114]

[0115] in, express The transpose of k Indicates the amount of observation data.

[0116] Compared with the prior art, the present invention has the following beneficial effects:

[0117] 1. This invention uses laser ranging and intersatellite ranging to jointly determine orbits, achieving precise orbit determination for a lunar formation of satellites. Compared with traditional orbit determination methods, the accuracy can be improved from hundreds of meters to meters.

[0118] 2. The present invention utilizes inter-satellite ranging information to reduce the requirements for the number and distribution of ground stations, thus alleviating the measurement and control pressure of existing ground stations;

[0119] 3. The present invention introduces laser ranging information, which effectively solves the problem of overall rotation of the constellation.

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

[0121] Other features, objects and advantages of the present invention will become more apparent upon reading the detailed description of non-limiting embodiments with reference to the following drawings:

[0122] Figure 1 It is the overall flow chart of the present invention;

[0123] Figure 2 This is a diagram of the measurement model of the system of the present invention;

[0124] Figure 3 This is a data preprocessing flow chart of the present invention;

[0125] Figure 4 This is a flow chart of the track improvement process of the present invention. DETAILED DESCRIPTION

[0126] The present invention will be described in detail below with reference to specific embodiments. The following examples 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 skilled in the art, several changes and improvements can be made without departing from the scope of the present invention. These all fall within the scope of protection of the present invention.

[0127] An 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. It can also significantly reduce the number of deployed ground tracking stations and lower the requirements for the geometric distribution of ground stations, becoming a viable means to address the current problems of low orbit determination accuracy of lunar-orbiting satellites and resource limitations of ground tracking stations.

[0128] 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.

[0129] The ground-based laser station adjusts the telescope's pointing direction based on satellite orbit predictions, transmits laser light toward and receives reflected laser light from the corner reflectors on Class A satellites, and records the times of laser emission and reception. The position and velocity of major celestial bodies in the solar system are obtained using the DE430 ephemeris.

[0130] Reference Figure 1 As shown, the specific method for orbit determination includes:

[0131] Step S1: Establish a high-precision dynamic model of the lunar satellite;

[0132] 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.

[0133] The non-spherical perturbation of the moon is the effect of the gravitational field change caused by the uneven distribution of the moon's mass on the orbiter's motion;

[0134] The third-body gravitational perturbation caused by the Sun and the Earth is the additional disturbance effect of the Sun's and Earth's gravitational pull on the lunar orbiter.

[0135] Lunar solid tidal perturbations are the effects of periodic changes in the gravitational field caused by the tidal deformation of the moon due to the gravitational pull of the Earth on the orbiter.

[0136] 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.

[0137] Solar pressure perturbation is the motion disturbance caused by the tiny pressure exerted by solar radiation photons on the orbiter surface.

[0138] Jet unloading is the instantaneous effect of the reaction force on the orbiter's trajectory when the orbiter uses the jet device to unload the flywheel angular momentum;

[0139] For lunar satellites, the coordinate system is the epoch moon-center inertial coordinate system. The high-precision dynamic model of lunar satellites is as follows:

[0140]

[0141] 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 system noise.

[0142] Step S2: Establishing mathematical models for lunar satellite laser ranging and inter-satellite ranging;

[0143] Perform data preprocessing on the laser ranging distance information to generate data files in the format required for orbit calculation;

[0144] 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).

[0145] The lunar satellite laser ranging observation quantity is the distance from the laser station to the satellite. Let the original observation distance be , the observation equation is:

[0146]

[0147] Where, The ranging error caused by the change of station position due to the earth's tides; 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;

[0148] 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:

[0149]

[0150] Where, for The relative distance between satellites A and B at time instant, i.e. the observed value of inter-satellite ranging; are the transmission delay and reception delay of Class A satellite respectively; 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.

[0151] 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;

[0152] The N satellites in the satellite formation are respectively combined with satellites in the satellite formation other than the current satellite to obtain N×(N-1) satellite combinations; each satellite combination includes the current satellite and the link establishment satellite. The N satellites in the satellite formation are all used as the current satellite, and one current satellite corresponds to N-1 satellite combinations; N>1.

[0153] Step S4: Preprocessing the absolute and relative position data of the satellites measured by laser ranging and inter-satellite ranging to generate a data file in the format required for orbit calculation;

[0154] The satellite's absolute and relative position data are preprocessed using observation equations, including:

[0155] 1) Space-time conversion: To ensure the accuracy of error correction, all calculations must be performed in the solar system barycenter coordinate system;

[0156] 2) System error correction: In laser ranging, the system error is the delay error of the ground station laser ranging system; in inter-satellite ranging, the system error is the transmission and reception delay of Class A and Class B satellites;

[0157] 3) Other error corrections:

[0158] Tropospheric refraction correction: Based on the temperature, humidity and pressure data measured at the observation station, the atmospheric refraction correction model is used to correct the tropospheric refraction of the satellite ranging and velocity measurements.

[0159] Relativistic delay correction: Based on the distance between the satellite and each celestial body, the bending of light caused by the gravitational field of the celestial body is calculated to complete the relativistic delay correction;

[0160] 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 extreme tide and ocean extreme 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;

[0161] Center of mass offset correction: The center of mass offset correction is completed based on the distance between the on-board corner reflector and the satellite center of mass;

[0162] Satellite antenna phase center correction: correct the satellite antenna phase offset;

[0163] 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.

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

[0165] Calculating the position coordinates of each satellite in the lunar formation in step S5 includes:

[0166] Step S5.1: Calculate the initial orbit using orbit measurement data combined with satellite dynamics equations;

[0167] Step S5.2: Using the initial orbital information to perform orbital integration, calculate the reference orbital and state transfer matrix;

[0168] Step S5.3: Linearize the observation equation according to the reference orbit and state transfer matrix;

[0169] Step S5.4: Perform linear optimal estimation on the linearized observation equation to obtain solution parameters;

[0170] The state differential equation of satellite motion is expressed as:

[0171]

[0172] in,

[0173] Where, and Represents 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 arc segment at the initial moment; Indicates satellite position; Indicates the satellite speed; represents the satellite acceleration; Represents other parameters to be estimated in the dynamic model;

[0174] In the above formula In the reference state Expand and remember , after omitting the higher-order terms, it can be expressed as a linear equation:

[0175]

[0176] Where, ;

[0177] Its solution is expressed as: ;

[0178] Solution ;

[0179] in, represents the state transition 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;

[0180] the satellite at the time the observation at the time , represents the measurement noise; represents the state vector of the satellite at the time the observation at the time represents the observation time; represents the observation data corresponding true value; i represents the i-th data;

[0181] Expanding the above formula at the reference state, only considering the first order term has:

[0182]

[0183] wherein, represents the observation partial derivative of the observation quantity to the state quantity at the observation epoch, represents the observation partial derivative to the improved epoch state quantity; represents the actual observation value;

[0184] The linear equation is obtained;

[0185] wherein, represents the state quantity to be estimated; represents the random error.

[0186] Step S5.5: calculating the theoretical observation value according to the solved parameters;

[0187] This step S5.5 includes:

[0188] solving the best estimate value of , and using the least square method to estimate the parameters; the weight matrix of the observation quantity is , and the estimate value is obtained according to the linear unbiased minimum variance estimation:

[0189]

[0190] wherein, represents the transpose of ; and k represents the observation data quantity.

[0191] Step S5.6: calculating the residual error, eliminating the observation value whose residual error exceeds the threshold value, and obtaining the orbit determination result of each satellite.

[0192] In the application, for the A-type satellite:

[0193] The initial orbit of the satellite is calculated based on the Laplace method by using the orbit measurement data and combining the satellite dynamics equation;

[0194] In the orbit determination arc, the observation values ​​are used to determine the orbit and calculate the orbit parameters of each satellite;

[0195] Calculate theoretical observation values ​​of the observation values ​​according to the solution parameters;

[0196] Residuals are calculated based on the observed values ​​and the theoretical observed values, and observations with residuals exceeding a threshold are eliminated.

[0197] In the orbit determination arc, the observation values ​​are used to determine the orbit and obtain the solution parameters, including:

[0198] Using the initial orbit information to perform orbital integration, we can obtain the reference orbit and state transfer matrix sampled at certain time intervals.

[0199] Linearizing the observation equation according to the reference orbit and the state transfer matrix;

[0200] The linear optimal estimation method is used to solve the observation equation after linearization to obtain the solution parameters.

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

[0202] For Class B satellites:

[0203] Based on the initial orientation elements and intersatellite measurement data provided by Class A satellites and the satellite dynamics equation, the initial satellite orbit is calculated based on the Laplace method;

[0204] In the orbit determination arc, the observation values ​​are used to determine the orbit and calculate the orbit parameters of each satellite;

[0205] Calculate theoretical observation values ​​of the observation values ​​according to the solution parameters;

[0206] Residuals are calculated based on the observed values ​​and the theoretical observed values, and observations with residuals exceeding a threshold are eliminated.

[0207] In the orbit determination arc, the observation values ​​are used to determine the orbit and obtain the solution parameters, including:

[0208] Using the initial orbit information to perform orbital integration, we can obtain the reference orbit and state transfer matrix sampled at certain time intervals.

[0209] The observation equation is linearized according to the reference orbit and state transfer matrix;

[0210] The linear optimal estimation method is used to solve the observation equation after linearization to obtain the solution parameters.

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

[0212] Each satellite in the satellite formation is measured and updated as the current satellite, and the orbit determination of each satellite is completed.

[0213] The application also provides a lunar-orbiting satellite formation orbit determination system, which can be realized by performing the flow steps of the lunar-orbiting satellite formation orbit determination method, i.e., the lunar-orbiting satellite formation orbit determination method can be understood by those skilled in the art as a preferred embodiment of the lunar-orbiting satellite formation orbit determination system. The system comprises: dividing the lunar-orbiting satellites in the lunar-orbiting satellite formation into A-type satellites and B-type satellites; configuring the A-type satellites with laser corner reflectors to carry out laser ranging tests in cooperation with ground laser stations; configuring the A-type satellites with integrated communication machines to carry out inter-satellite measurement technology tests in cooperation with the B-type satellites; configuring the B-type satellites with navigation and control loads to carry out inter-satellite measurement technology tests in cooperation with the A-type satellites; and the orbit determination system comprises:

[0214] Module M1: establishing a lunar-orbiting satellite high-precision dynamics model;

[0215] Specifically, module M1 comprises: for the lunar-orbiting satellite, selecting a coordinate system as a lunar center inertial coordinate system, and the lunar-orbiting satellite high-precision dynamics model is as follows:

[0216]

[0217] wherein, is the differential of ; is a state vector of the state model; are the positions and velocities of the lunar satellite in the X, Y, and Z directions, respectively; is a system nonlinear continuous state transition function of the state model; represents system noise;

[0218] Module M2: establishing a lunar-orbiting satellite laser ranging and inter-satellite ranging mathematical model;

[0219] Module M2 comprises:

[0220] The lunar satellite laser ranging observation quantity is the distance from the laser station to the satellite, and the original observation distance is set as , and the observation equation is:

[0221]

[0222] wherein, is the ranging error caused by the position change of the station due to the earth tide; is the ranging error caused by the refraction effect of light in the atmosphere; is the ranging 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;

[0223] 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:

[0224]

[0225] Where, for The relative distance between satellites A and B at time instant, i.e. the observed value of inter-satellite ranging; are the transmission delay and reception delay of Class A satellite respectively; 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.

[0226] Module M3: measures the distance, i.e., the absolute position, of a Class A satellite relative to a ground station through laser ranging, and measures the relative position of a Class B satellite relative to a Class A satellite through inter-satellite ranging;

[0227] Module M4: Pre-processes the absolute and relative position data of satellites measured by laser ranging and inter-satellite ranging to generate data files in the format required for orbit calculation;

[0228] Module M4 pre-processes the absolute and relative position data of the satellite through observation equations, including:

[0229] 1) Space-time conversion: To ensure the accuracy of error correction, all calculations must be performed in the solar system barycenter coordinate system;

[0230] 2) System error correction: In laser ranging, the system error is the delay error of the ground station laser ranging system; in inter-satellite ranging, the system error is the transmission and reception delay of Class A and Class B satellites;

[0231] 3) Other error corrections:

[0232] Tropospheric refraction correction: Based on the temperature, humidity and pressure data measured at the observation station, the atmospheric refraction correction model is used to correct the tropospheric refraction of the satellite ranging and velocity measurements.

[0233] Relativistic delay correction: Based on the distance between the satellite and each celestial body, the bending of light caused by the gravitational field of the celestial body is calculated to complete the relativistic delay correction;

[0234] 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 extreme tide and ocean extreme 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;

[0235] Center of mass offset correction: The center of mass offset correction is completed based on the distance between the on-board corner reflector and the satellite center of mass;

[0236] Satellite antenna phase center correction: correct the satellite antenna phase offset;

[0237] 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.

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

[0239] Module M5 calculates the position coordinates of each satellite in the lunar formation, including:

[0240] Module M5.1: Calculate the initial orbit using orbital measurement data combined with satellite dynamics equations;

[0241] Module M5.2: Use the initial orbital information to perform orbital integration and calculate the reference orbit and state transfer matrix;

[0242] Module M5.3: Linearize the observation equations based on reference orbits and state transfer matrices;

[0243] Module M5.4: Perform linear optimal estimation on the linearized observation equation to obtain the solution parameters;

[0244] The M5.4 module includes:

[0245] The state differential equation of satellite motion is expressed as:

[0246]

[0247] in,

[0248] Where, and Represents 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 arc segment at the initial moment; Indicates satellite position; Indicates the satellite speed; represents the satellite acceleration; Represents other parameters to be estimated in the dynamic model;

[0249] In the above formula In the reference state Expand and remember , after omitting the higher-order terms, it can be expressed as a linear equation:

[0250]

[0251] Where, ;

[0252] Its solution is expressed as: ;

[0253] Solution ;

[0254] in, represents the state transition 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;

[0255] Satellite in Observable quantity at a moment , 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;

[0256] Expand the above formula at the reference state and consider only the first-order terms:

[0257]

[0258] Where, represents the observed partial derivative of the observed quantity with respect to the state quantity at the observation epoch, represents the observed partial derivative of the improved epoch state quantity; represents the actual observed value;

[0259] Get the linear equation ;

[0260] in, Represents the state quantity to be estimated; represents random error;

[0261] Module M5.5: Calculate theoretical observations based on solution parameters;

[0262] The M5.5 module includes:

[0263] Solution The best estimate of is obtained by using the least squares method for parameter estimation; the weight matrix of the observation is recorded as , and the estimated value is obtained based on the linear unbiased minimum variance estimate:

[0264]

[0265] in, express The transpose of k Indicates the amount of observation data.

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

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

[0268] The present invention provides a method for determining the orbit of a lunar satellite formation, which specifically includes:

[0269] Precise orbit determination involves using precise dynamics to fit satellite orbits using various observational data containing observational errors. Obviously, the accuracy of orbit determination calculations is significantly limited by the accuracy of the measurement data, the accuracy of the observation model, and the accuracy of the dynamic model, in addition to the observational geometry. Precise orbit determination involves processing various types of observational data, including accurate observational modeling and correction for various measurement errors. Satellites in different orbital types are subject to different dynamic influences, necessitating different dynamic models to be considered in orbit determination calculations. Furthermore, both the dynamic model and the observation model contain certain errors, with the dynamic model being particularly significant. Therefore, the orbit calculation process requires the selection of an appropriate mathematical model and the optimization of the measurement data to improve the accuracy of the calculation results.

[0270] Figure 1 This is a flow chart of an embodiment of the 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 diagram of the system measurement model involved in the combined orbit determination method of lunar satellite formation laser ranging and inter-satellite ranging. Figure 1 and Figure 2 The laser ranging and intersatellite ranging combined orbit determination method includes:

[0271] Step 1: Establish a high-precision dynamic model of the lunar satellite.

[0272] For lunar satellites, the coordinate system is the lunar-centered inertial coordinate system of the epoch (J2000.0). Taking into account factors such as the gravitational forces of the Sun, Moon, Earth, and other celestial bodies on the lunar satellite, the lunar solid tide perturbation, the lunar physical libration, solar radiation pressure, jet unloading, and high-precision ephemeris, the high-precision dynamic model of the lunar satellite is shown below:

[0273]

[0274] The above formula can be simplified as:

[0275]

[0276] 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 barycenter 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 solar barycenter coordinate system; are the system noise respectively.

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

[0278] The lunar satellite laser ranging observation quantity is the distance from the laser station to the satellite. Let the original observation distance be , the observation equation is:

[0279]

[0280] Where, The ranging error caused by the change of station position due to the earth's tides; 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;

[0281] 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:

[0282]

[0283] Where, for The relative distance between satellites A and B at time instant, i.e. the observed value of inter-satellite ranging; are the transmission delay and reception delay of Class A satellite respectively; 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.

[0284] Step 3: The ground laser station conducts a laser ranging test to obtain the distance information of the 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 the Class B satellite relative to the Class A satellite.

[0285] Step 4: Use the observation equation in step 2 to preprocess the original position information obtained in step 3:

[0286] The specific steps of data preprocessing are as follows: Figure 3 Shown, including:

[0287] (1) Space-time conversion:

[0288] To ensure the accuracy of error correction, all calculations must be performed in the solar system barycenter coordinate system;

[0289] (2) System error correction:

[0290] 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;

[0291] (3) Error correction:

[0292] 1. Tropospheric refraction correction;

[0293] 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.

[0294] 2. Relativistic delay correction;

[0295] Based on the distance between the satellite and each celestial body, the bending of light caused by the gravitational field of the celestial body is calculated to complete the relativistic delay correction.

[0296] 3. Earth tide correction;

[0297] Based on the solid tide, ocean tide, atmospheric load tide generated by the external gravitational force on the earth, and the parameters of the solid extreme tide and ocean extreme 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.

[0298] 4. Correction of center of mass offset;

[0299] The center of mass offset correction is completed based on the distance between the on-board corner reflector and the satellite center of mass.

[0300] 5. Satellite antenna phase center correction;

[0301] Corrects the satellite antenna phase offset.

[0302] (4) Data format conversion:

[0303] After the orbit measurement data is corrected and processed, the format conversion is performed to generate the data content in the format required for orbit determination calculation.

[0304] Step 5: Calculate the position coordinates of each satellite in the formation, such as Figure 4 As shown, the following steps are included:

[0305] (1) Initial orbit calculation: using orbit measurement data to calculate and generate the satellite's initial orbit information for orbit improvement:

[0306] Based on the Laplace method, the ranging data is used to calculate the initial satellite orbit. The dynamic model used in the initial orbit calculation is the model in step 1, which is applicable to satellites in lunar orbit.

[0307] (2) Use the initial orbital information to perform orbital integration and calculate the reference orbit and state transfer matrix;

[0308] Calculating Satellite Orbits Using Numerical Integration Methods

[0309]

[0310] 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 high-order error term;

[0311] Using the orbital integration method, with the initial orbital parameters as the initial values, the trajectory obtained by integration calculation is called the reference orbit;

[0312] The state differential equation of satellite motion is expressed as:

[0313]

[0314] in, ;

[0315] Where, and Represents 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 arc segment at the initial moment; Indicates satellite position; Indicates the satellite speed; represents the satellite acceleration; Represents other parameters to be estimated in the dynamic model;

[0316] (3) Linearize the observation equation according to the reference orbit and state transfer matrix;

[0317] In the above formula In the reference state Expand and remember , after omitting the higher-order terms, it can be expressed as a linear equation:

[0318]

[0319] Where, ;

[0320] Its solution is expressed as: ;

[0321] The solution of the state equation is ;

[0322] in, represents the state transition 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;

[0323] (4) Perform linear optimal estimation on the linearized observation equation to obtain the solution parameters;

[0324] Satellite in Observable quantity at a moment , 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;

[0325] Expand the above formula at the reference state and consider only the first-order terms:

[0326]

[0327] Where, represents the observed partial derivative of the observed quantity with respect to the state quantity at the observation epoch, represents the observed partial derivative of the improved epoch state quantity; represents the actual observed value;

[0328] Get the linear equation ;

[0329] in, Represents the state quantity to be estimated; Represents random error.

[0330] Step 5: Calculate the theoretical observation value based on the solution parameters.

[0331] Solution The best estimate of is obtained by using the least squares method for parameter estimation; the weight matrix of the observation is recorded as , the estimated value of the theoretical observation is obtained according to the linear unbiased minimum variance estimate:

[0332]

[0333] in, express The transpose of k Indicates the amount of observation data.

[0334] 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.

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

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

[0337] An embodiment of the present invention provides a method and system for determining the orbit of a lunar satellite formation. It 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.

[0338] Those skilled in the art will appreciate that, in addition to implementing the system and its various devices, modules, and units provided by the present invention in purely computer-readable program code, it is entirely possible to implement the same functions of the system and its various devices, modules, and units provided by the present invention 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 a hardware component, and the devices, modules, and units included therein for implementing various functions can also be considered as structures within the hardware component; the devices, modules, and units for implementing various functions can also be considered as both software modules implementing the method and structures within the hardware component.

[0339] The above describes specific embodiments of the present invention. It should be understood that the present invention is not limited to the specific embodiments described above, and those skilled in the art may make various changes or modifications within the scope of the claims, which do not affect the essence of the present invention. The embodiments of this application and the features in the embodiments may be combined with each other in any manner unless there is a conflict.

Claims

1. A method for determining the orbit of a lunar satellite formation, characterized in that: The lunar orbit satellites in the lunar orbit satellite formation are divided 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 equipment 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 includes: Step S1: Establish a high-precision dynamic model of the lunar satellite; Step S2: Establishing mathematical models for 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: Preprocessing the absolute and relative position data of the satellites measured by laser ranging and inter-satellite ranging to generate a data file in the 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, wherein: The step S1 includes: for a lunar satellite, the coordinate system is selected as the epoch moon-center inertial coordinate system, and the high-precision dynamic model of the lunar satellite is as follows: in, Differentiation of; X = [xyzv x v y v z ] T is the state vector of the state model; x, y, z, v x , v y ,v z are the position and velocity of the lunar satellite in the X, Y, and Z directions respectively; f(X, t) is the system nonlinear continuous state transfer function of the state model; W(t) represents the system noise.

3. The method for determining the orbit of a lunar satellite formation according to claim 1, wherein: The step S2 includes: the lunar satellite laser ranging observation quantity is the distance from the laser station to the satellite, assuming the original observation distance is ρ′, the observation equation is: ρ0=ρ′-(Δρ tide +Dr atm +Dr rel +Dr mc +Dr sys ) Where Δρ tide The ranging error caused by the change of station position due to the earth tide; Δρ atm The distance measurement error caused by the refraction effect of light in the atmosphere; Δρ rel The distance measurement error caused by the relativistic effect of light in the gravitational field; Δρ mc is the deviation of the laser reflection point on the satellite surface from the center of mass; Δρ sys 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: Where ρ(t0) is the relative distance between satellites A and B at time t0, i.e., the observed value of inter-satellite ranging; are the transmission delay and reception delay of Class A satellite respectively; are the transmission delay and reception delay of Class B satellite respectively; Δρ corr 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, wherein: The step S4 pre-processes the absolute position and relative position data of the satellite using 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 ground station laser ranging system; in inter-satellite ranging, the system error is the transmission and reception delay of Class A and Class B satellites; 3) Other error corrections: Tropospheric refraction correction: Based on the temperature, humidity and pressure data measured at the observation station, the atmospheric refraction correction model is used to correct the tropospheric refraction of the satellite ranging and velocity measurements. Relativistic delay correction: Based on the distance between the satellite and 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 extreme tide and ocean extreme 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 based on the distance between the on-board corner reflector and 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, wherein: 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: Using the initial orbital information to perform orbital integration, calculate the reference orbital and state transfer matrix; Step S5.3: Linearize the observation equation according to the reference orbit and 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 ​​based on 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, wherein: The step S5.4 includes: the state differential equation of satellite motion is expressed as: in, Where X and represents the state quantity to be estimated; F represents the dynamic parameter; X(t0) represents the state quantity to be estimated at the initial time of the orbit determination arc; t0 represents the initial time of the orbit determination arc; X0 represents the state quantity at the initial time of the orbit determination arc; r represents the satellite position; represents the satellite velocity; a represents the satellite acceleration; p represents other parameters to be estimated in the dynamic model; Expand X in the above equation at the reference state X*, record x=XX*, and omit the high-order terms to express it as a linear equation: Where, The solution is: x(t) = Φ(t, t0)x(t0); Solution Among them, Φ(t, t0) represents the state transfer matrix, represents the derivative of the state transfer matrix; I represents the identity matrix; X* represents the reference state; x(t) represents the difference between the state to be estimated and the reference state; t represents time; represents the derivative of the difference between the state to be estimated and the reference state; x(t0) represents the difference between the state to be estimated and the reference state at time t0; Satellite in t i The observed quantity Y at time i =G(t i , X i )+ε i , ε i represents the measurement noise; X i Indicates that the satellite is at t i State vector at time t i Indicates the observation time; G() indicates the observation data Y i The corresponding true value; i represents the i-th data; Expand the above formula at the reference state and consider only the first-order terms: Where, represents the observed partial derivative of the observed quantity with respect to the state quantity at the observed epoch, H represents the observed partial derivative of the state quantity with respect to the improved epoch; Y represents the actual observed value; The linear equation y = Hx0 + ε is obtained; Among them, x0 represents the state quantity to be estimated; ε represents the random error.

7. The method for determining the orbit of a lunar satellite formation according to claim 6, wherein: The step S5.5 includes: solving the best estimate of x0, using the least squares method to estimate the parameters; denoting the weight matrix of the observed quantity as P, and obtaining the estimate according to the linear unbiased minimum variance estimation: Among them, H T represents the transpose of H; 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 and cooperate with ground laser stations to carry out laser ranging tests; the Class A satellites are equipped with integrated communication equipment 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 orbit determination system includes: Module M1: Establishing a high-precision dynamic model of the lunar satellite; Module M2: Establish mathematical models for lunar satellite laser ranging and intersatellite ranging; Module M3: measures the distance, i.e., the absolute position, of a Class A satellite relative to a ground station through laser ranging, and measures the relative position of a Class B satellite relative to a Class A satellite through inter-satellite ranging; Module M4: Pre-processes the absolute and relative position data of satellites 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 lunar satellite formation orbit determination system according to claim 8, characterized in that: The module M1 includes: for the lunar satellite, the coordinate system is the epoch moon-center inertial coordinate system, and the high-precision dynamic model of the lunar satellite is as follows: in, is the differential of X(t); X=[xyzv x v y v z ] T is the state vector of the state model; x, y, z, v x , v y , v z are the position and velocity of the lunar satellite in the X, Y, and Z directions respectively; f(X, t) is the system nonlinear continuous state transfer function of the state model; W(t) represents the system noise; The module M2 includes: The lunar satellite laser ranging observation quantity is the distance from the laser station to the satellite. Let the original observation distance be ρ′, and the observation equation is: ρ0=ρ′-(Δρ tide +Dr atm +Dr rel +Dr mc +Dr sys ) Where Δρ tide The ranging error caused by the change of station position due to the earth tide; Δρ atm The distance measurement error caused by the refraction effect of light in the atmosphere; Δρ rel The distance measurement error caused by the relativistic effect of light in the gravitational field; Δρ mc is the deviation of the laser reflection point on the satellite surface from the center of mass; Δρ sys 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: Where ρ(t0) is the relative distance between satellites A and B at time t0, i.e., the observed value of inter-satellite ranging; are the transmission delay and reception delay of Class A satellite respectively; are the transmission delay and reception delay of Class B satellite respectively; Δρ corr 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 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 ground station laser ranging system; in inter-satellite ranging, the system error is the transmission and reception delay of Class A and Class B satellites; 3) Other error corrections: Tropospheric refraction correction: Based on the temperature, humidity and pressure data measured at the observation station, the atmospheric refraction correction model is used to correct the tropospheric refraction of the satellite ranging and velocity measurements. Relativistic delay correction: Based on the distance between the satellite and 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 extreme tide and ocean extreme 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 based on the distance between the on-board corner reflector and 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 lunar satellite formation orbit determination system according to claim 8, characterized in that: Module M5 calculates the position coordinates of each satellite in the lunar formation, including: Module M5.1: Calculate the initial orbit using orbital 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 reference orbits and state transfer matrices; 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 observations with residuals exceeding 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, Where X and represents the state quantity to be estimated; F represents the dynamic parameter; X(t0) represents the state quantity to be estimated at the initial time of the orbit determination arc; t0 represents the initial time of the orbit determination arc; X0 represents the state quantity at the initial time of the orbit determination arc; r represents the satellite position; represents the satellite velocity; a represents the satellite acceleration; p represents other parameters to be estimated in the dynamic model; Expand X in the above equation at the reference state X*, record x=XX*, and omit the high-order terms to express it as a linear equation: Where, The solution is: x(t) = Φ(t, t0)x(t0); Solution Among them, Φ(t, t0) represents the state transfer matrix, represents the derivative of the state transfer matrix; I represents the identity matrix; X* represents the reference state; x(t) represents the difference between the state to be estimated and the reference state; t represents time; represents the derivative of the difference between the state to be estimated and the reference state; x(t0) represents the difference between the state to be estimated and the reference state at time t0; Satellite in t i The observed quantity Y at time i =G(t i , X i )+ε i , ε i represents the measurement noise; X i Indicates that the satellite is at t i State vector at time t i Indicates the observation time; G() indicates the observation data Y i The corresponding true value; i represents the i-th data; the above formula is expanded at the reference state, considering only the first-order terms: Where, represents the observed partial derivative of the observed quantity with respect to the state quantity at the observed epoch, H represents the observed partial derivative of the state quantity with respect to the improved epoch; Y represents the actual observed value; The linear equation y = Hx0 + ε is obtained; Among them, x0 represents the state quantity to be estimated; ε represents the random error; The module M5.5 includes: To find the best estimate of x0, the least squares method is used for parameter estimation; the weight matrix of the observation is denoted as P, and the estimate is obtained according to the linear unbiased minimum variance estimation: Among them, H T represents the transpose of H; k represents the amount of observation data.

Images