A method for estimating the orbit determination accuracy of earth-moon space optical observation

Abstract

The application discloses a kind of lunar space optical observation orbit determination precision estimation method, solve the problem that orbit determination precision evaluation needs to be calculated by a large number of Monte Carlo simulation, causes the problem of calculation resource and time consumption, belong to orbit determination precision estimation field;Including: using circular restricted three-body problem to describe the gravitational model of lunar space, construct the dynamics equation of target in lunar space under conjunction coordinate system;Linearized dynamics equation and the differential equation that state transition matrix satisfies are obtained;Construct target imaging signal-to-noise ratio model and apparent magnitude model to obtain target observation standard, construct the Jacobian matrix of observation model according to the relative position relationship of target and sensor;Based on posterior Cramer-Lo bound theory, combined with the differential equation that state transition matrix satisfies and the Jacobian matrix of observation model, construct orbit determination precision estimation index model, the lower bound of position error and the lower bound of velocity error are used as precision estimation index in the application, reduce the amount of calculation, effectively reflect the observation orbit determination effect.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of orbit determination accuracy estimation, and relates to a method for estimating orbit determination accuracy of optical observation in the earth-moon space. BACKGROUND

[0002] The earth-moon space has become the latest destination of human spaceflight, and the American Artemis program will build a "Gateway" space station on the moon and has completed the Artemis I related task. The Chinese Chang'e series of tasks have completed the lunar landing and return, and are promoting the first lunar landing plan in the 21st century. With the commercialization of the space industry, space activities related to the moon have transitioned from government-led in the past to commercial companies. Recently, the moon lander "Blue Ghost" of the American private company "Firefly" Aerospace successfully landed in the "Mare Serenitatis" basin on the near side of the moon. This is the first completely successful soft landing on the moon by a commercial company in human history. It can be predicted that human activities in the earth-moon space will increase dramatically, and with them come the potential threat of uncleanable space debris and near-earth asteroids to human activities, so it is necessary to enhance the observation capability of the earth-moon space.

[0003] However, due to the radius of the earth-moon space being more than ten times that of the geostationary orbit, ground radar equipment is difficult to meet the observation of objects in this space due to power limitations and distance attenuation. In addition, due to the limitation of the moon's exclusion angle, ground-based optical equipment will have a period of time in a synodic month cycle during which it cannot observe the space around the moon. For this reason, many scholars have carried out research on space-based optical observation of the earth-moon space.

[0004] Currently, when performing simulation work, the state estimation of the target is affected by the measurement error itself, has certain randomness, and needs to be simulated by Monte Carlo, resulting in consumption of more computing resources and computing time; for example, the application number is 202211481656.0, the publication number is CN115905800A, and the invention name is a method and system for analyzing the accuracy of space target orbit determination based on optical observation. In the method and system, a plurality of sets of simulation data under different random noises are generated by using random numbers to analyze the accuracy of space target orbit determination, the target position estimation value is obtained by parameter solving, and the average position deviation is generated by comparing the actual position, and the process is repeated for statistical analysis.

[0005] The prior art needs a large amount of computing hardware resources and a long computing time, and the computing result needs to be processed manually. If the result does not meet the requirements, the input conditions need to be adjusted and the long-time repeated calculation needs to be performed again, which delays the research progress. SUMMARY

[0006] In order to solve the technical problem that the evaluation of orbit determination accuracy requires a large number of Monte Carlo simulation calculations that consume a lot of computing resources and computing time, the present invention proposes a method for estimating orbit determination accuracy through optical observation of Earth-Moon space. This method uses the industrial Internet of Things information perception technology and adopts the posterior Cramer-Rao lower bound as the basis for determining the accuracy of orbit determination. By converting statistical data and processing industrial information and data, the position accuracy estimation index and velocity accuracy estimation index for orbit determination calculation using the extended Kalman filter algorithm are obtained. This index is only related to the relative state between the target and the sensor and the error level of the sensor itself, and is not affected by the error in the simulation; long-term simulation can be performed. The method works to determine the position of a target in orbit for a long time; the a posteriori Cramer-Rao lower bound represents an unbiased estimate that can be achieved mathematically for orbit determination, and effectively reflects the observational orbit determination effect of the target; by establishing a target imaging signal-to-noise ratio model and an apparent magnitude model, the target visibility is jointly evaluated from the imaging signal-to-noise ratio requirement and the apparent magnitude model, assuming that the solar phase angle is less than or equal to 90° and the target is outside the celestial body exclusion angle, and then the a posteriori Cramer-Rao lower bound is used to evaluate the theoretical accuracy of the target orbit determination, and an estimated value of the calculation result is derived through a mathematical formula; the present invention has low requirements on computing hardware and consumes little computing time, and can flexibly change the input conditions to quickly obtain an estimate of the result.

[0007] The purpose of the present invention is specifically achieved through the following technical solutions:

[0008] The present invention discloses a method for estimating orbit determination accuracy by optical observation in Earth-Moon space, the method comprising:

[0009] Step 1: Use the circularly restricted three-body problem to describe the gravitational model of Earth-Moon space. Based on the gravitational model, construct the rendezvous coordinate system and the dynamic equations of the target in the Earth-Moon space under the rendezvous coordinate system.

[0010] Step 2: Perform Taylor series expansion on the dynamic equation at the known trajectory and linearize it by ignoring the high-order terms greater than or equal to the second order to obtain the linearized dynamic equation;

[0011] Step 3: Combine the dynamic equation and the linearized dynamic equation to obtain the differential equation satisfied by the state transfer matrix;

[0012] Step 4: Build a target imaging signal-to-noise ratio model using the CCD equation; use a Lambertian sphere with both diffuse and specular reflection as the target and build a magnitude model; use the target imaging signal-to-noise ratio calculated by the target imaging signal-to-noise ratio model as a target observation criterion for judging whether the target is visible when the signal-to-noise ratio is greater than a threshold, the solar phase angle in the magnitude model is less than or equal to 90°, and the target is outside the celestial exclusion angle;

[0013] Step 5: When the target is visible, the observation model and the Jacobian matrix of the observation model in the rendezvous coordinate system are constructed according to the relative position relationship between the target and the sensor;

[0014] Step 6: Based on the a posteriori Cramer-Rao lower bound theory, combined with the differential equation satisfied by the state transfer matrix and the Jacobian matrix of the observation model, an orbit determination accuracy estimation index model is constructed. The position error lower bound in the orbit determination accuracy estimation index model is used as the position accuracy estimation index for orbit determination calculation using the extended Kalman filter algorithm; the velocity error lower bound in the orbit determination accuracy estimation index model is used as the velocity accuracy estimation index for orbit determination calculation using the extended Kalman filter algorithm.

[0015] In step 1, the circular restricted three-body problem is used to describe the gravitational model of the Earth-Moon space. The method for constructing the rendezvous coordinate system based on the gravitational model is as follows:

[0016] In the Earth-Moon space, the Earth and the Moon are regarded as the two main celestial bodies, the target is regarded as the third body, the unit mass is the sum of the masses of the two main celestial bodies, the unit distance is the distance between the two main celestial bodies, and the period of the two main celestial bodies rotating around their common center of mass is 2π. The gravitational model of the Earth-Moon space described by the circular restricted three-body problem is obtained; the rendezvous coordinate system constructed based on the gravitational model is:

[0017] The origin of the conjunction coordinate system is the common center of mass of the Earth and the Moon. The x-axis points from the center of mass to the Moon and rotates with the Moon. The y-axis is determined by the Cartesian coordinate system. The z-axis is perpendicular to the plane of the celestial body's rotation. The mass of the Earth is , the mass of the moon is , the coordinates of the Earth are , the coordinates of the moon are , the target position coordinates are , where is the quality coefficient, where .

[0018] In step 1, the dynamic equation of the target in the Earth-Moon space in the rendezvous coordinate system is:

[0019] ;

[0020] Where, is the acceleration coordinate of the target, is the velocity coordinate of the target, is the distance from the target to the Earth, is the distance from the target to the moon; .

[0021] In step 2, the linearized dynamic equation is:

[0022] ;

[0023] Where, is a nonlinear function describing the system dynamics, is the state vector in the dynamics equation, ; is the initial state of the known trajectory, is the initial state vector of the known trajectory, , is the initial position coordinate of the target, is the initial velocity coordinate of the target, T is the transpose; is the Jacobian matrix, is the time variable, for and The state deviation between them.

[0024] In step 3, the differential equation satisfied by the state transfer matrix is:

[0025] ;

[0026] Where, is the derivative of the state transfer matrix, is the state transition matrix, is the initial time.

[0027] In step 4, the target imaging signal-to-noise ratio model constructed by the CCD equation is:

[0028] ;

[0029] Where, is the signal term, is the noise term; The number of signal electrons received by the CCD per second. is the integration time of CCD receiving photons; is the pixel occupied by the target image, is the number of signal electrons received by the CCD per second from the sky background. is the number of dark current electrons generated by the CCD per second, is the CCD signal read noise.

[0030] In step 4, the magnitude model is:

[0031] ;

[0032] Where, is the apparent magnitude, is the apparent magnitude of the Sun, is the target reflection cross-sectional area, is the average reflectivity, are the coefficients of diffuse reflection and specular reflection, is the phase angle function of diffuse reflection, is the phase angle function of specular reflection, is the solar phase angle, i.e. the angle between the sun, target and sensor, is the distance between the target and the sensor.

[0033] In step 4, in the target observation standard, the solar phase angle is: ; when When the angle is greater than 90°, imaging is impossible;

[0034] ; is the angle between celestial body, sensor and target, where celestial body is the sun, earth or moon; if If the exclusion angle is smaller than that of the celestial body, imaging cannot be performed;

[0035] Where, is the vector of the target pointing to the sensor, is the target vector pointing to the sun, is the vector pointing from the sensor to the target, is the vector pointing the sensor toward the sun, earth, or moon, To obtain the vector modulus length operation.

[0036] In step 5, the observation model constructed in the conjunction coordinate system is:

[0037] ;

[0038] Where, for The target azimuth change rate over time, for The target pitch angle change rate over time; is the target azimuth measured by the sensor in the conjunction coordinate system, The angle between the projection on the xy plane and the positive direction of the x-axis; is the target pitch angle measured by the sensor in the conjunction coordinate system, and The angle of the projection on the xy plane; , The vector pointing from the sensor to the target The coordinates of for About time variables The derivative of , The velocity coordinate of the sensor pointing to the target; for The module length, ;

[0039] Jacobian matrix of the constructed observation model for:

[0040] .

[0041] In step 6, the orbit determination accuracy estimation index model is:

[0042] ;

[0043] in, is the lower bound of the position error, ; is the lower bound of the velocity error ;

[0044] Where, is the absolute position error, is the absolute velocity error, To find the expected operation, is an unbiased estimate of the target position state, is an unbiased estimate of the target velocity state, is the posterior Clamer-Rao lower bound matrix, which is a 6×6 matrix. for The first row and first column of The estimated lower bound of for The second row and second column of The estimated lower bound of for The third row and third column of The estimated lower bound of for The 4th row and 4th column of The estimated lower bound of for The 5th row and 5th column of The estimated lower bound of for The 6th row and 6th column of The estimated lower bound of .

[0045] The beneficial effects of the present invention are:

[0046] (1) The present invention processes industrial information and data through mathematical formulas and performs conversion statistics. This method is only related to the relative state between the target and the sensor and the error level of the sensor itself, and is not affected by the errors in the simulation. It solves the technical problem that the evaluation of orbit determination accuracy requires a large number of Monte Carlo simulations, which consumes a lot of computing resources and time.

[0047] (2) It can conduct long-term simulation work to determine the long-term target position in orbit;

[0048] (3) Through the information perception technology of the Industrial Internet of Things, the posterior Cramer-Rao lower bound is used to represent the unbiased estimate value that can be achieved mathematically, which effectively reflects the observation and orbit determination effect of the target.

[0049] (4) By establishing a target imaging signal-to-noise ratio model and a visual magnitude model, the target visibility is jointly evaluated from the imaging signal-to-noise ratio requirement and the solar phase angle less than or equal to 90° and the target is outside the celestial body exclusion angle in the visual magnitude model. Then, the theoretical accuracy of the target orbit determination is evaluated using the posterior Cramer-Rao lower bound. The technical solution of giving an estimated value of the calculation result by mathematical formula deduction has low requirements for computing hardware, does not rely on a large number of statistical calculations, consumes less computing time, can flexibly change the input conditions and quickly obtain an estimate of the result, significantly reducing computing consumption. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] The present invention will be described in further detail below with reference to the accompanying drawings and examples.

[0051] Figure 1 Schematic diagram of the relative positions of the Earth, the Moon, and a target in a conjunction coordinate system provided by an embodiment of the present invention.

[0052] Figure 2 Schematic diagram of the solar phase angle provided by an embodiment of the present invention.

[0053] Figure 3 This is a schematic diagram of the celestial body exclusion angle provided by an embodiment of the present invention.

[0054] Figure 4 It is a schematic diagram of azimuth and elevation angles in a conjunction coordinate system provided by an embodiment of the present invention.

[0055] Figure 5 Schematic diagram of a target and sensor orbit determination scenario provided by an embodiment of the present invention.

[0056] Figure 6 3 is a schematic diagram comparing the orbit determination result and the orbit determination accuracy estimation result provided by an embodiment of the present invention.

[0057] Figure 7 3 is a schematic diagram comparing the orbit determination velocity result and the orbit determination velocity accuracy estimation result provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0058] In order to enable those skilled in the art to better understand the solutions of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the embodiments described are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts should fall within the scope of protection of the present invention.

[0059] An embodiment of the present invention provides a method for estimating orbit determination accuracy by optical observation in Earth-Moon space, the method comprising:

[0060] Step 1: Use the circularly restricted three-body problem to describe the gravitational model of Earth-Moon space. Based on the gravitational model, construct the rendezvous coordinate system and the dynamic equations of the target in the Earth-Moon space under the rendezvous coordinate system.

[0061] This step uses the circularly restricted three-body problem to describe the gravitational model of Earth-Moon space. The equations of motion describe the motion of the target and sensor, enabling a description of how their positions and velocities change over time. In the circularly restricted three-body problem, the two primary celestial bodies, with their greater mass, move in a circular motion under the mutual gravitational pull of each other. The third celestial body, with its much smaller mass, is affected by the gravitational pull of the two primary bodies but has no effect on their motion. In Earth-Moon space, the two primary celestial bodies are the Earth and the Moon, and the target is the third body.

[0062] Step 2: Perform Taylor series expansion on the dynamic equation at the known trajectory and linearize it by ignoring the high-order terms greater than or equal to the second order to obtain the linearized dynamic equation;

[0063] Step 3: Combine the dynamic equation and the linearized dynamic equation to obtain the differential equation satisfied by the state transfer matrix;

[0064] Step 4: Build a target imaging signal-to-noise ratio model using the CCD equation; use a Lambertian sphere with both diffuse and specular reflection as the target and build a magnitude model; use the target imaging signal-to-noise ratio calculated by the target imaging signal-to-noise ratio model as a target observation criterion for judging whether the target is visible when the signal-to-noise ratio is greater than a threshold, the solar phase angle in the magnitude model is less than or equal to 90°, and the target is outside the celestial exclusion angle;

[0065] Step 5: When the target is visible, the observation model and the Jacobian matrix of the observation model in the rendezvous coordinate system are constructed according to the relative position relationship between the target and the sensor;

[0066] Step 6: Based on the a posteriori Cramer-Rao lower bound theory, combined with the differential equation satisfied by the state transfer matrix and the Jacobian matrix of the observation model, an orbit determination accuracy estimation index model is constructed. The position error lower bound in the orbit determination accuracy estimation index model is used as the position accuracy estimation index for orbit determination calculation using the extended Kalman filter algorithm; the velocity error lower bound in the orbit determination accuracy estimation index model is used as the velocity accuracy estimation index for orbit determination calculation using the extended Kalman filter algorithm.

[0067] In step 1, the circular restricted three-body problem is used to describe the gravitational model of the Earth-Moon space. The method for constructing the rendezvous coordinate system based on the gravitational model is as follows:

[0068] In the Earth-Moon space, the Earth and the Moon are regarded as the two main celestial bodies, the spacecraft is regarded as the third body, the unit mass is the sum of the masses of the two main celestial bodies, the unit distance is the distance between the two main celestial bodies, and the period of the two main celestial bodies rotating around their common center of mass is 2π. The gravitational model of the Earth-Moon space described by the circular restricted three-body problem is obtained; Figure 1 As shown, the conjunction coordinate system constructed based on the gravitational model is:

[0069] The origin of the conjunction coordinate system is the common center of mass of the Earth and the Moon. The x-axis points from the center of mass to the Moon and rotates with the Moon. The y-axis is determined by the Cartesian coordinate system. The z-axis is perpendicular to the plane of the celestial body's rotation. The mass of the Earth is , the mass of the moon is , the coordinates of the Earth are , the coordinates of the moon are , the target position coordinates are , where is the quality coefficient, where .

[0070] In step 1, the dynamic equation of the target in the Earth-Moon space in the rendezvous coordinate system is:

[0071] ;

[0072] Where, is the acceleration coordinate of the target, is the velocity coordinate of the target, is the distance from the target to the Earth, is the distance from the target to the moon; .

[0073] In step 2, the linearized dynamic equation is:

[0074] ;

[0075] Where, is a nonlinear function describing the system dynamics, is the state vector in the dynamics equation, ; is the initial state of the known trajectory, is the initial state vector of the known trajectory, , is the initial position coordinate of the target, is the initial velocity coordinate of the target, T is the transpose; is the Jacobian matrix, is the time variable, for and The linearized dynamic equation is derived as follows:

[0076] The original dynamic equation of the target in the Earth-Moon space in the conjunction coordinate system can be simplified as:

[0077] ;

[0078] Assume there is a known trajectory , then the state vector and The deviation between them is:

[0079] ;

[0080] Will exist Performing Taylor series expansion around and ignoring higher-order terms greater than or equal to second order, we can obtain:

[0081] ;

[0082] in, is the Jacobian matrix denoted as ;The linearized dynamic equation is: .

[0083] In step 3, the differential equation satisfied by the state transfer matrix is:

[0084] ;

[0085] Where, is the derivative of the state transfer matrix, is the state transition matrix, is the initial time. The derivation process of the differential equation satisfied by the state transfer matrix is ​​as follows:

[0086] Substituting the linearized dynamic equation into the original dynamic equation yields:

[0087] ;

[0088] Because the trajectory is known satisfy ;

[0089] therefore, ;

[0090] Then the linearized deviation dynamic equation is: ;

[0091] State transition matrix Describes state deviation From the initial time To time variable The linear transformation of is defined as:

[0092] ;right Taking the derivative, we get ;

[0093] Then, substitute the linearized deviation dynamic equation into the above equation and replace the left side of the equal sign ,get:

[0094] ;

[0095] Substitute the definition of the state transfer matrix into the above formula and replace the left side of the equal sign ,get:

[0096] ;

[0097] because Can be the initial deviation at any time, so the two sides of the equal sign It can be eliminated, and then the quantities on both sides of the equal sign are swapped to obtain the differential equation satisfied by the state transfer matrix as follows:

[0098] ;

[0099] According to the definition of the state transfer matrix, when At the initial moment, the state transfer matrix satisfies:

[0100] ;

[0101] That is, at the initial moment, the state deviation does not change, so the initial moment state transfer matrix is the identity matrix . And for any moment The state transition matrix , can be solved by numerically integrating the differential equation of the state transfer matrix. In the present invention, the state transfer matrix Used in the subsequent accuracy estimation process.

[0102] In step four, the target imaging signal-to-noise ratio model constructed by the CCD equation is as follows:

[0103]

[0104] In the formula, is a signal term, is a noise term; is the number of signal electrons received by the CCD per second, is the integration time of the CCD receiving photons; is the number of pixels occupied by the target image, is the number of signal electrons received by the CCD per second, is the number of dark current electrons generated by the CCD per second, is the CCD signal reading noise.The derivation process of the target imaging signal-to-noise ratio model constructed by the CCD equation is as follows:

[0105] By constructing a high-fidelity space-based optical observation model, the visibility of the sensor to the target during space-based optical observation is reflected, and the sensor can only implement measurement when the target is visible. The space-based visible light detection payload generally uses a CCD camera, and the signal-to-noise ratio (SNR) is an important indicator affecting the capability of the payload and is directly related to the photon energy received by the detection system. Since space-based observation is less disturbed by various factors, the CCD equation used in the present application is as follows:

[0106]

[0107] In the formula, is a signal term, is a noise term, is the number of signal electrons generated by the target, is the number of pixels occupied by the target image, is the number of signal electrons generated by the sky background, is the dark current generated by the CCD itself, is the CCD signal reading noise. Considering that the integration time of the CCD receiving photons is , in the above formula, , , are respectively as follows:

[0108]

[0109] In the formula, and are respectively the number of signal electrons received by the CCD per second, is the number of dark current electrons generated by the CCD itself per second, and the target imaging signal-to-noise ratio model constructed by the CCD equation is as follows: ​​

[0110] ;

[0111] Where, and The CCD manufacturer provides the relevant hardware parameters, and 、 Determined by the observed target and the target sky background, where

[0112] ;

[0113] ;

[0114] Where, is the luminous flux density of a target with a visible magnitude of 0, is the atmospheric penetration rate, is the optical transmittance of the system, is the transmittance of the secondary mirror, is the telescope aperture, is the effective focal length of the telescope, is the average quantum efficiency of the CCD, is the pixel size, is the sky background brightness.

[0115] In step 4, the magnitude model is:

[0116] ;

[0117] Where, is the apparent magnitude, is the apparent magnitude of the Sun, is the target reflection cross-sectional area, is the average reflectivity, are the coefficients of diffuse reflection and specular reflection, is the phase angle function of diffuse reflection, is the phase angle function of specular reflection, is the solar phase angle, i.e. the angle between the sun, target and sensor, is the distance between the target and the sensor.

[0118] ;

[0119] .

[0120] In step 4, in the target observation standard, such as Figure 2 As shown, the solar phase angle is: ; when When the angle is greater than 90°, imaging is impossible;

[0121] like Figure 3As shown, ; is the angle between celestial body, sensor and target, where celestial body is the sun, earth or moon; if If the exclusion angle is smaller than that of the celestial body, imaging cannot be achieved;

[0122] Where, is the vector of the target pointing to the sensor, is the target vector pointing to the sun, is the vector pointing from the sensor to the target, is the vector pointing the sensor toward the sun, earth, or moon, To obtain the vector modulus length operation;

[0123] The visibility of the sensor to the target is determined according to the target observation standard. Only when the target is visible can the target be observed, thereby generating observation data.

[0124] When the target is visible, the sensor measures the target and generates observation data, which is calculated by the measurement model. Therefore, it is necessary to establish a measurement model, such as Figure 4 As shown in Figure 2, in step 5, the observation model constructed is:

[0125] ;

[0126] Where, for The target azimuth change rate over time, for The target pitch angle change rate over time; is the target azimuth measured by the sensor in the conjunction coordinate system, The angle between the projection on the xy plane and the positive direction of the x-axis; is the target pitch angle measured by the sensor in the conjunction coordinate system, and The angle of the projection on the xy plane; , The vector pointing from the sensor to the target The coordinates of for About time variables The derivative of , The velocity coordinate of the sensor pointing to the target; for The module length, ;

[0127] Jacobian matrix of the constructed observation model for:

[0128] ;

[0129] ;

[0130] Where, .

[0131] In step 6, the orbit determination accuracy estimation index model is:

[0132] ;

[0133] Among them, As the lower bound of the position error; As the lower bound of the velocity error;

[0134] Where, is the absolute position error, is the absolute velocity error, To find the expected operation, is an unbiased estimate of the target position state, is an unbiased estimate of the target velocity state, is the posterior Clamer-Rao lower bound matrix, which is a 6×6 matrix. for The first row and first column of The estimated lower bound of for The second row and second column of The estimated lower bound of for The third row and third column of The estimated lower bound of for The 4th row and 4th column of The estimated lower bound of for The 5th row and 5th column of The estimated lower bound of for The 6th row and 6th column of The derivation process of the orbit determination accuracy estimation index model is as follows:

[0135] To evaluate orbit determination accuracy, observation data containing random noise is often used to perform orbit calculations, and the calculated results are compared with the true orbit values. The Monte Carlo method is then used to minimize the effects of randomness. This method is beneficial for evaluating results in specific orbit determination scenarios and can effectively simulate the effects of noise during actual observations. However, it suffers from high computational resource consumption and randomness in the evaluation results. Therefore, a deterministic metric is needed to evaluate orbit determination accuracy. The Posterior Cramer-Rao Lower Bound (PCRLB) provides an unbiased lower bound for the target error covariance. This method is used in this paper as the basis for evaluating orbit determination accuracy. The expression is:

[0136] ;

[0137] In the formula, the subscript k represents the ordinal number of the current calculation, is the error covariance matrix of the k-th step, is the state variable, is the estimated value of the variable, is the posterior Fisher information matrix, That is the posterior Cramer-Rao lower bound matrix of the kth step, and its recursive calculation process is:

[0138] ;

[0139] For the orbit determination process using the extended Kalman filter, the differential equation of the state transfer matrix is ​​substituted into , and the Jacobian matrix of the observation model , introducing the measurement error matrix , process noise matrix , the elements in the above formula can be expressed as:

[0140] ;

[0141] because 、 is a diagonal matrix. Substituting the above formula into the recursive calculation process formula, we can get:

[0142] ;

[0143] According to the matrix inversion lemma (Woodbury matrix identity), we can get:

[0144] ; simplifies to:

[0145] ;

[0146] Assumptions represents the observation matrix of the i-th sensor to the target at the k-th step. The case where the same target is observed by multiple different sensors can be expressed as:

[0147] ;

[0148] In orbit determination, As a 6×6 matrix, the absolute position error can be obtained , absolute speed error The relationship with the posterior Cramer-Rao lower bound matrix can be expressed as:

[0149] ;

[0150] definition 、 They represent the lower bounds of position error and velocity error, respectively, which are used as the orbit determination accuracy estimation indicators adopted by the present invention.

[0151] In order to explain the technical solution of the present invention, a specific example is provided for illustration:

[0152] Assume a scenario for the orbit determination of the Near Direct Halo Orbit (NRHO) in Earth-Moon space. Figure 5 As shown in the figure, the black track represents the target's motion trajectory, the tracks of other colors are the sensor's motion trajectory, the black dot is the sensor's initial position, and the gray dot is the moon.

[0153] Figure 5 The motion of the sensor and the target follows the dynamic equations, and the visibility of the target relative to the sensor complies with the target observation standard for judging whether the target is visible.

[0154] When the target is visible relative to the sensor, simulated observation data can be generated based on the measurement model, including the azimuth of the target relative to the sensor. , pitch angle , azimuth change rate and pitch angle change rate According to the observation data, the current common extended Kalman filter model can be used to obtain the orbit position and velocity results of the target, such as Figure 6 and Figure 7 Shown as solid line.

[0155] Using the orbit determination accuracy estimation index model, the position and velocity accuracy estimation indicators are calculated and , we can get the orbit position and velocity accuracy estimation effect as follows Figure 6 and Figure 7 As shown by the dotted line. Figure 6 and Figure 7It can be found that the accuracy estimation result has a good inclusion relationship with the error in the actual orbit determination process, that is, the actual orbit determination error is within the range of the accuracy estimation result, which shows the effectiveness of the accuracy estimation.

[0156] The beneficial effects of the embodiments of the present invention are:

[0157] (1) The present invention processes industrial information and data through mathematical formulas and performs conversion statistics. This method is only related to the relative state between the target and the sensor and the error level of the sensor itself, and is not affected by the errors in the simulation. It solves the technical problem that the evaluation of orbit determination accuracy requires a large number of Monte Carlo simulations, which consumes a lot of computing resources and time.

[0158] (2) It can conduct long-term simulation work to determine the long-term target position in orbit;

[0159] (3) Through the information perception technology of the Industrial Internet of Things, the posterior Cramer-Rao lower bound is used to represent the unbiased estimate value that can be achieved mathematically, which effectively reflects the observation and orbit determination effect of the target.

[0160] (4) By establishing a target imaging signal-to-noise ratio model and a visual magnitude model, the target visibility is jointly evaluated from the imaging signal-to-noise ratio requirement and the solar phase angle less than or equal to 90° and the target is outside the celestial body exclusion angle in the visual magnitude model. Then, the theoretical accuracy of the target orbit determination is evaluated using the posterior Cramer-Rao lower bound. The technical solution of giving an estimated value of the calculation result by mathematical formula deduction has low requirements for computing hardware, does not rely on a large number of statistical calculations, consumes less computing time, can flexibly change the input conditions and quickly obtain an estimate of the result, significantly reducing computing consumption.

[0161] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in the present invention should be included in the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be based on the scope of protection of the claims.

Claims

1. A method for estimating orbit determination accuracy by optical observation in Earth-Moon space, characterized in that: The method includes: Step 1: Use the circularly restricted three-body problem to describe the gravitational model of Earth-Moon space. Based on the gravitational model, construct the rendezvous coordinate system and the dynamic equations of the target in the Earth-Moon space under the rendezvous coordinate system. Step 2: Perform Taylor series expansion on the dynamic equation at the known trajectory and linearize it by ignoring the high-order terms greater than or equal to the second order to obtain the linearized dynamic equation; Step 3: Combine the dynamic equation and the linearized dynamic equation to obtain the differential equation satisfied by the state transfer matrix; Step 4: Build a target imaging signal-to-noise ratio model using the CCD equation; use a Lambertian sphere with both diffuse and specular reflection as the target and build a magnitude model; use the target imaging signal-to-noise ratio calculated by the target imaging signal-to-noise ratio model as a target observation criterion for judging whether the target is visible when the signal-to-noise ratio is greater than a threshold, the solar phase angle in the magnitude model is less than or equal to 90°, and the target is outside the celestial exclusion angle; Step 5: When the target is visible, the observation model and the Jacobian matrix of the observation model in the rendezvous coordinate system are constructed according to the relative position relationship between the target and the sensor; Step 6: Based on the a posteriori Cramer-Rao lower bound theory, combined with the differential equation satisfied by the state transfer matrix and the Jacobian matrix of the observation model, an orbit determination accuracy estimation index model is constructed. The position error lower bound in the orbit determination accuracy estimation index model is used as the position accuracy estimation index for orbit determination calculation using the extended Kalman filter algorithm; the velocity error lower bound in the orbit determination accuracy estimation index model is used as the velocity accuracy estimation index for orbit determination calculation using the extended Kalman filter algorithm.

2. The method according to claim 1, wherein In step 1, the circular restricted three-body problem is used to describe the gravitational model of the Earth-Moon space. The method for constructing the rendezvous coordinate system based on the gravitational model is as follows: In the Earth-Moon space, the Earth and the Moon are regarded as the two main celestial bodies, the target is regarded as the third body, the unit mass is the sum of the masses of the two main celestial bodies, the unit distance is the distance between the two main celestial bodies, and the period of the two main celestial bodies rotating around their common center of mass is 2π. The gravitational model of the Earth-Moon space described by the circular restricted three-body problem is obtained; the rendezvous coordinate system constructed based on the gravitational model is: The origin of the conjunction coordinate system is the common center of mass of the Earth and the Moon. The x-axis points from the center of mass to the Moon and rotates with the Moon. The y-axis is determined by the Cartesian coordinate system. The z-axis is perpendicular to the plane of the celestial body's rotation. The mass of the Earth is , the mass of the moon is , the coordinates of the Earth are , the coordinates of the moon are , the target position coordinates are , where is the quality coefficient, where .

3. The method according to claim 2, wherein In step 1, the dynamic equation of the target in the Earth-Moon space in the rendezvous coordinate system is: ； Where, is the acceleration coordinate of the target, is the velocity coordinate of the target, is the distance from the target to the Earth, is the distance from the target to the moon; .

4. The method according to claim 3, wherein In step 2, the linearized dynamic equation is: ； Where, is a nonlinear function describing the system dynamics, is the state vector in the dynamics equation, ; is the initial state of the known trajectory, is the initial state vector of the known trajectory, , is the initial position coordinate of the target, is the initial velocity coordinate of the target, T is the transpose; is the Jacobian matrix, is the time variable, for and The state deviation between them.

5. The method according to claim 4, wherein In step 3, the differential equation satisfied by the state transfer matrix is: ； Where, is the derivative of the state transfer matrix, is the state transition matrix, is the initial time.

6. The method according to claim 5, wherein In step 4, the target imaging signal-to-noise ratio model constructed by the CCD equation is: ； Where, is the signal term, is the noise term; The number of signal electrons received by the CCD per second. is the integration time of CCD receiving photons; is the pixel occupied by the target image, is the number of signal electrons received by the CCD per second from the sky background. is the number of dark current electrons generated by the CCD per second, is the CCD signal read noise.

7. The method according to claim 6, wherein In step 4, the magnitude model is: ； Where, is the apparent magnitude, is the apparent magnitude of the Sun, is the target reflection cross-sectional area, is the average reflectivity, are the coefficients of diffuse reflection and specular reflection, is the phase angle function of diffuse reflection, is the phase angle function of specular reflection, is the solar phase angle, i.e. the angle between the sun, target and sensor, is the distance between the target and the sensor.

8. The method according to claim 7, wherein In step 4, in the target observation standard, the solar phase angle is: ; when When the angle is greater than 90°, imaging is impossible; ; is the angle between celestial body, sensor and target, where celestial body is the sun, earth or moon; if If the exclusion angle is smaller than that of the celestial body, imaging cannot be performed; Where, is the vector of the target pointing to the sensor, is the target vector pointing to the sun, is the vector pointing from the sensor to the target, is the vector pointing the sensor toward the sun, earth, or moon, To obtain the vector modulus length operation.

9. The method according to claim 8, wherein In step 5, the observation model constructed in the conjunction coordinate system is: ； Where, for The target azimuth change rate over time, for The target pitch angle change rate over time; is the target azimuth measured by the sensor in the conjunction coordinate system, The angle between the projection on the xy plane and the positive direction of the x-axis; is the target pitch angle measured by the sensor in the conjunction coordinate system, and The angle of the projection on the xy plane; , The vector pointing from the sensor to the target The coordinates of for About time variables The derivative of , The velocity coordinate of the sensor pointing to the target; for The module length, ; Jacobian matrix of the constructed observation model for: 。 10. The method according to claim 9, wherein In step 6, the orbit determination accuracy estimation index model is: ； in, is the lower bound of the position error, ; is the lower bound of the velocity error ; Where, is the absolute position error, is the absolute velocity error, To find the expected operation, is an unbiased estimate of the target position state, is an unbiased estimate of the target velocity state, is the posterior Clamer-Rao lower bound matrix, which is a 6×6 matrix. for The first row and first column of The estimated lower bound of for The second row and second column of The estimated lower bound of for The third row and third column of The estimated lower bound of for The 4th row and 4th column of The estimated lower bound of for The 5th row and 5th column of The estimated lower bound of for The 6th row and 6th column of The estimated lower bound of .

Images