Planetary gravity field inversion method

Abstract

The invention belongs to the cross technical field of satellite gravity mechanics, geodetics, orbit perturbation, space science and the like. The invention particularly relates to a planet gravity field inversion method. The planet gravity field inversion method comprises the following steps: S1, acquiring shadow region duration of a planet in an orbital period through orbital device load data; s2, acquiring an included angle beta (t) between a sun-planet vector and an orbiter plane through a conical earth shadow model or a cylindrical earth shadow model according to the shadow region duration time tau (t), and further acquiring an ascending node right ascension advancing amount omega (t); s3, solving the ascending node right ascension advancement amount omega (t) to obtain a long-term item and a long-period item; and S4, obtaining a low-order even-number term harmonic coefficient and a low-order odd-number term harmonic coefficient of the gravity field by using a linear perturbation method, an average perturbation method or an analytical method. According to the invention, the low-order gravity field can be measured only by using the sun sensitive load carried by the orbiter without radio measurement, and long-distance measurement can be realized.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of satellite gravimetry, geodesy, orbital perturbation and space science, and particularly relates to a planetary gravity field inversion method. BACKGROUND

[0002] Planetary gravity field measurement is one of the hotspots and frontiers in the field of deep space exploration, and is an important means to understand the material distribution and evolution of a planet. At present, the planetary gravity field is obtained by deep space exploration technologies such as radio ranging, radio velocity measurement and very long baseline interferometry (VLBI) to obtain the orbit position information of a spacecraft, for example, patent document CN114740541A, and the main advantage is high precision and suitable for long-distance tasks within hundreds of millions of kilometers. However, the radio signal is affected by the interference in the propagation path (such as the atmosphere and solar plasma). Taking Mars as an example, the distance between Mars and the Earth changes from 50 million to 500 million kilometers, and the signal strength decreases sharply with the distance according to the inverse square law, which worsens the ranging and velocity measurement accuracy. In addition, the global coverage of radio observation is difficult, and the accuracy of the non-conservative force model of the planet is poor, which makes it difficult to determine the global gravity field of the planet. Therefore, it is of great significance to directly measure the planetary gravity field information by using the orbiters' own load. SUMMARY

[0003] In view of the defects of the prior art, the application provides a new planetary gravity field inversion method, which does not rely on radio measurement results, but only relies on the load data of the orbiters, and aims to solve the problem of radio measurement error caused by long distance.

[0004] To achieve the above purpose, the application provides a planetary gravity field inversion method for obtaining low-order zonal harmonic coefficients of the gravity field, comprising the following steps: S1. In a circular orbit, the duration of the shadow zone of the planet in an orbit period is obtained by the load data of the orbiters. ; wherein, μ G is the gravitational constant of the planet, R P R is the average radius of the planet, h H is the orbital height, θ is the angular distance of the orbiters, and t is the running time; S2. According to the duration of the shadow zone τ (t), the angle between the sun-planet vector and the plane of the orbiters is obtained by a conical shadow model or a cylindrical shadow model β (t), and the ascending node precession Ω (t) is obtained. S3. The ascending node precession obtaining long-term terms and long-period terms ; S4. obtaining low-order even-degree zonal coefficients of the gravity field by linear perturbation method or average perturbation method according to long-term terms and long-period terms obtaining low-order odd-degree zonal coefficients of the gravity field by linear perturbation method or analytical method.

[0005] Preferably, in the step S1, the duration of the shadow zone of the planet in an orbital period is obtained by the method of STK simulation.

[0006] Preferably, the orbiter is a low-inclination orbiter, and in the step S2, the duration of the shadow zone is obtained by a conical shadow model, the formula of the conical shadow model being: β (t). .

[0007] Preferably, the orbiter is a high-inclination orbiter, and in the step S2, the duration of the shadow zone is obtained by a cylindrical shadow model, the formula of the cylindrical shadow model being: β (t).

[0008] wherein, R S is the average radius of the sun, A is the average distance between the sun and the planet.

[0009] Preferably, in the step S2, the amount of precession of the right ascension of the ascending node is obtained, and the formula of the amount of precession of the right ascension of the ascending node is: Ω (t). ; wherein, ε is the obliquity of the ecliptic; Г is the true solar longitude of the ecliptic, i is the orbital inclination.

[0010] Preferably, in the step S3, the effect of the solar radiation pressure and the atmospheric drag on the amount of precession of the right ascension of the ascending node is deducted by the non-conservative force acceleration measured by the accelerometer.

[0011] As a further preferred, , ; wherein, is the average root number without perturbation, is the short-period term, ​​Ω is the mean motion of the ascending node due to solar radiation pressure, Ω is the mean motion of the ascending node due to atmospheric drag, r Ω is the planet-orbit vector, f Ω is the true anomaly, ζ Ω is the mean angular velocity of the orbit, ω Ω is the argument of perigee, a Ω is the semi-major axis, Y Ω is the acceleration component of the non-conservative force acceleration measured by the accelerometer along the normal of the orbit plane.

[0012] Preferably, in the step S4, the low-order even zonal harmonic coefficients of the gravity field are obtained by the average perturbation method J n wherein n is an even number.

[0013] As further preferred, n is 2 and 4, the low-order even zonal harmonic coefficients of the gravity field are obtained by the method J 2 and J 4, specifically as follows:

[0014] ; wherein, G G is the gravitational constant, i i is the orbit inclination, M P M is the planet mass, a a is the semi-major axis.

[0015] 1. The precision of the inversion of the gravity field is mainly affected by the noise level of the payload itself, only by means of the measurement of the payload data of the orbit, so that the present application can realize long-distance measurement; 2. Since the technical means of the duration of the shadow zone τ (t) measurement is adopted, it is not dependent on radio observation; 3. Since the orbit can be equipped with an accelerometer, the influence of the solar radiation pressure and the atmospheric drag on the ascending node right ascension is deducted by the non-conservative force acceleration measured by the accelerometer; when an accelerometer with a precision of 10 -12 m / s 2 order is used as a solar sensitive payload, since the accelerometer can measure the non-conservative force suffered by the orbit, it can measure and deduct the non-conservative force perturbation while sensing the duration of the shadow zone, so that the model can be measured and deducted at the same time, with high precision. BRIEF DESCRIPTION OF DRAWINGS

[0016] Fig. 1 The relevant angle definitions of the present application are shown in the figure, wherein φis the angle between the sun-planet-orbiters vector and the orbiters orbital plane; β is the angle between the sun-planet-orbiters vector and the orbiters orbital plane; θ is the angular distance the orbiters are running, S is the planet-star vector, r is the planet-orbiters vector.

[0017] Fig. 2a represents the conical shadow zone model; Fig. 2b represents the cylindrical shadow zone model; Fig. 3 represents the schematic diagram of the STK simulation of embodiment 1; Fig. 4 represents the shadow zone duration of embodiment 1 and embodiment 2 τ simulation results; Fig. 5 represents the power spectral density of τ of embodiment 1; Fig. 6 represents the orbiters β angle simulation results of embodiment 1 and embodiment 2; Fig. 7 represents the power spectral density of orbiters β angle of embodiment 1 and embodiment 2; Fig. 8 represents the sine value of the ascending node precession of embodiment 1 and embodiment 2 Ω ); Fig. 9 represents the power spectral density of Ω ) of embodiment 1 and embodiment 2; Fig. 10 represents the flow of embodiment 1 of the present application. DETAILED DESCRIPTION

[0018] In order to make the purposes, technical solutions and advantages of the present application clearer, further detailed description will be made to the present application in combination with the drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application, and are not used to limit the present application.

[0019] In the description of the present application, it should be understood that the terms "first", "second" are only used for the purpose of description, and cannot be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Therefore, the features limited by "first", "second" can explicitly or implicitly include one or more of the features. In the description of the present application, the meaning of "multiple" is two or more, unless otherwise specifically limited.

[0020] Furthermore, throughout this specification, references to "an embodiment"; "an embodiment," "an example," or similar language indicate that a particular feature, structure, or characteristic described in connection with that embodiment is included in at least one embodiment of this application. Therefore, the appearance of the phrase "in one embodiment;" throughout this specification, and similar language, may, but not necessarily, refer to the same embodiment.

[0021] In light of the above background, this application proposes a method for determining a planet's gravitational field by measuring the duration of shadowing in an orbiter's shadow region, with a focus on obtaining low-order gravitational fields. J 2. J Information on fourth-order coefficients. Lower-order gravitational field coefficients cause orbital perturbations in satellites. For circular orbit satellites, this primarily manifests as the precession of the right ascension of the ascending node, causing variations in the duration of shadows per orbit. Therefore, by measuring this variation, lower-order gravitational field coefficients can be estimated. Since the spacecraft's entry into and exit from shadow areas alters the non-conservative forces acting on it, the resulting changes in shadow duration can be measured using accelerometers. Additionally, based on the changes in sunlight and temperature during spacecraft entry into and exit from shadow areas, the resulting changes in shadow duration can also be obtained using star sensors, temperature sensors, or accelerometers.

[0022] To achieve the above objectives, this application provides a method for inverting planetary gravitational fields to obtain low-order band harmonic coefficients of the gravitational field, comprising the following steps: S1. Under circular orbit, obtain the duration of the planet's shadowed region within one orbital period using orbiter payload data. ;in, μ The gravitational constant of a planet, R P The average radius of the planet. h For the orbital height, θ The angle is the angular distance of the orbiter's movement, and t is the time of movement. Since a circular orbit is chosen, it is assumed that... h Keep it unchanged; then by setting the semi-major axis a eccentricity e Track inclination i Perigeal argument ω Right ascension of ascending node Ω Initial value, mean anterior angle M Using initial values ​​and methods such as STK simulation, the duration of the shadowed region of a planet within one orbital period is obtained. .

[0023] Relevant angles are defined as follows: Fig. 1 As shown; where, φ The angle is the angle subtended by the Sun-Planet-Orbiter; βis the angle between the sun-planet vector and the orbiter plane; θ is the angular distance of the orbiter, S is the planet-star vector, r is the planet-orbiter vector.

[0024] S2. According to the shadow zone duration τ (t), the angle between the sun-planet vector and the orbiter plane β (t) is obtained; Two shadow zone models are given in Fig. 2, which are Fig. 2a the conical shadow model and Fig. 2b the cylindrical shadow model.

[0025] When the orbiter is a high-inclination orbiter, the angle β (t) is obtained by the cylindrical shadow model, whose formula is:

[0026] wherein, R S is the average radius of the sun, A is the average distance between the sun and the planet; When the orbiter is a low-inclination orbiter, always holds, so the angle β (t) can be obtained by the conical shadow model, whose formula is: .

[0027] Further, the precession of the ascending node right ascension β (t) is obtained by Ω (t), specifically:

[0028] wherein, ε is the obliquity of the ecliptic; Г is the true sun longitude of the ecliptic, i is the orbital inclination, all of which are known quantities.

[0029] S3. The long-term term and the long-period term are obtained from the precession of the ascending node right ascension by spectral analysis; wherein the signals lower than 10 -5 Hz are the long-term term and its multiple frequency signals, that is ; the signals higher than 10 -5 Hz are the long-period term signals ; , wherein, is the mean root number of the camera, is the short-period term, is the precession of the ascending node due to solar radiation pressure, is the precession of the ascending node due to atmospheric drag; ; wherein, r is the planet-orbit vector, f is the true anomaly, ζ is the mean angular velocity of the orbit, ω is the argument of perigee, a is the semi-major axis, e is the eccentricity, Y is the acceleration component of the non-conservative force measured by the accelerometer along the normal direction of the orbital plane; thus, the acceleration component of the non-conservative force along the normal direction of the orbital plane can be measured by the accelerometer and subtracted from and .

[0030] S4. Obtaining low-order even-harmonic coefficients of the gravity field J n ; wherein n is the order, which is a natural number; the higher the accuracy to be solved, the higher the upper limit of n; wherein the low-order even-harmonic coefficients of the gravity field (i.e. n = 2, 4, 6,...) are obtained according to the long-term term by using the linear perturbation method or the average perturbation method; the low-order odd-harmonic coefficients of the gravity field (i.e. n = 3, 5,...) are obtained according to the long-period term by using the linear perturbation method or the analytical method.

[0031] Taking the low-order even-harmonic coefficients J 2 and J 4 obtained by the average perturbation method as examples, the method for obtaining the low-order even-harmonic coefficients J 2 and J 4 is as follows:

[0032] ; wherein, G is the gravitational constant, i is the orbital inclination, M P is the orbit plane anomaly, a is the semi-major axis, e is the eccentricity; since e = 0, the above formula can be further simplified as:

[0033] ; Other low order zonal coefficients of the gravitational field can be obtained by other methods.

[0034] The following is an example.

[0035] Example 1 The method of retrieving the gravitational field of a planet based on the variation of the shadow time of an orbiter contains the following steps: Step 1: Measure the time of the satellite in the shadow region at each orbit using the data of the payload of the orbiter.

[0036] As shown in Fig. 1, when using the conical shadow model, the orbiter is in three stages in one orbit period, i.e. the illumination region, the penumbra and the umbra. We want to record the time of the orbiter entering each stage using the payload data sensitive to the sun in the orbiter Fig. 2a τ where τ 1 represents the time from the center of the illumination region to the entry into the first penumbra, τ 2 represents the time from the center of the illumination region to the entry into the umbra, τ 3 represents the time from the center of the illumination region to the entry into the second penumbra, τ 4 represents the time from the center of the illumination region to the entry into the next illumination region, R e is the radius of the planet, and R s is the radius of the star relative to the planet. The payload sensitive to the sun will have a clear jump signal at the moment of entering and leaving the shadow region, such as the electrostatic accelerometer, which can measure the non-conservative force received by the orbiter in orbit. When the orbiter enters the penumbra from the illumination region, the solar radiation pressure begins to gradually decrease, and becomes 0 after entering the umbra. Therefore, we need to record the time when the payload measurement begins to decrease (increase) and stops decreasing (increasing). The duration of the orbiter in the entire shadow region in one period of the planet τ A and the time of the orbiter in the umbra τ U :

[0037]

[0038] This application selects the circular orbit model for calculation. Under the circular orbit, the orbiter moves at a constant speed, i.e. where μ is the gravitational constant of the planet,​R P is the mean radius of the planet, h is the orbit altitude, θ is the angular distance of the orbiter, and t represents the time of the operation. Since the selected orbit is a circular orbit, it is considered that h is constant.

[0039] As shown in Fig. 3 , we use STK aviation simulation to simulate the orbiter around Mars as the central celestial body. The pink dot in the figure represents the spatial position of the orbiter, which rotates around Mars in a circular orbit (i.e., the eccentricity e = 0), ensuring the relative stability of h in formulas (3) and (4); according to formulas (11) and (12), the lower the orbit altitude of the orbiter, the greater the orbit perturbation it receives, and the average radius of Mars is 3389.5 km, and the lowest altitude of the orbiter is about 300 km, so the semi-major axis a selected by us is 3700 km; according to formula (11), the smaller the orbit inclination, the greater the change in the precession, which means that the Earth's gravity will cause greater orbit perturbation, which is beneficial to the calculation of the present application, but at the same time, if the inclination is too small, the change amplitude of the perturbation will be small, which will affect the estimation accuracy, so the orbit inclination i is set to a low inclination (15 deg); it can also be seen from formulas (11) and (12) that the precession of the ascending node right ascension Ω is independent of the initial value of the argument of perigee ω , the mean anomaly M , and the ascending node right ascension Ω In order to facilitate calculation, we selected 0 deg as the initial value in this embodiment. The simulation parameter settings are shown in Table 1, and the simulation results are shown in Fig. 3 .

[0040] Table 1 Simulation parameters

[0041] After simulation, the duration of the orbiter in the shadow area in a year changes as shown in Fig. 4 .

[0042] Step 2: Estimate the precession of the ascending node right ascension of the orbiter according to the shadow area time.

[0043] 2.1 Selection of shadow model For the conical shadow model, according to the geometric relationship, and and the orbit β The relationship between the angles is:

[0044]

[0045] in A This represents the average distance between the Sun and the planets. R S This represents the average radius of the Sun. Since our simulation uses a low-Earth orbit satellite, we can approximate sunlight as parallel light for low-Earth orbit satellites, i.e., the average distance between the Sun and planets. A Infinity ,at this time The duration of the penumbra is negligible, and the cone model degenerates into a cylindrical model, such as... Fig. 2b As shown. Formulas (3) and (4) under the cylindrical model can be rewritten as:

[0046] Using the formula above, we measure the change in duration. Calculator orbiter orbit Changes in angle .

[0047] refer to Fig. 4 The time variation of the orbiter in the shaded area, measured under the simulation conditions of Example 1, is given. τ (t), where the horizontal axis represents the simulation runtime t, and the vertical axis represents the time the satellite is in the shadow region as measured by the orbiter payload. τ Variation with simulation time τ (t). According to formula (5), when Time (i.e., the orbit in this simulation) β When the angle exceeds 66.36°, the orbiter will undergo a full illumination phase ( τ =0), at this time the relationship of formula (5) no longer holds, but since we selected a low-inclination tracker in the simulation, according to formula (6), as long as the track inclination angle is less than 34°, β The maximum angle can never reach 66.36°, therefore the orbiter will not experience a fully illuminated phase. Our measurements are continuous, therefore in Fig. 5 The simulation is given in the middle. τ The power spectral density curve of (t) is plotted with frequency on the x-axis and duration on the y-axis. τ The power spectral density of (t) contains various frequency information, which are given by various perturbation effects.

[0048] 2.2 Calculation of the precession of the right ascension of the ascending node β The angle is the angle between the Sun-planet vector and the orbiter plane, which varies with time and can be expressed as follows:

[0049] where, ε is the obliquity of the ecliptic; Г is the true Sun longitude in ecliptic, i is the orbital inclination.

[0050] From the above formula, it can be seen that the influence of β The angle changes mainly due to two factors, respectively the change (caused by the revolution of the planet around the sun) and Ω the change (caused by the orbital perturbation). Therefore, we can extract the precession of the ascending node right ascension β (t) from the angle . Fig. 6 and Fig. 7 give the time-domain variation and power spectral density curve of the simulated β angle, the horizontal coordinate is frequency, and the vertical coordinate is β The power spectral density of the angle shows similar frequency components as τ , but due to the limitation of the measurement method, the τ obtained by β can only be sampled once every orbital period, and the orbital frequency signal component is missing, so the present application is suitable for observing long-term signals and long-period signals.

[0051] Step three: determine the gravity field coefficient according to the precession of the ascending node right ascension.

[0052] For the two-body problem, the motion equation of the orbiter can be expressed as:

[0053] where is the distance from the orbiter to the center of the planet, is the acceleration of the orbiter, is the acceleration of the orbiter due to the central gravity of the planet, is the perturbation acceleration of the orbit, which can be further expressed as follows according to the main perturbation sources of the Earth-like planet:

[0054] In the above formula, grad represents the gradient, represents the planetary non-spherical gravitational perturbation function, represents the third body gravitational perturbation function, represents the solar radiation pressure perturbation acceleration, represents the atmospheric drag perturbation acceleration, where can be represented as:

[0055] where and are the gravitational field coefficients, is the zonal harmonic, is the sectorial harmonic, m is the order, n is the degree; P nm is the associated Legendre function, P n is the Legendre function, φ is the geocentric latitude, λ is the longitude, λ nm is the phase angle. For a planet like Earth, the largest perturbation factor is the zonal harmonic (first order small quantity), other zonal harmonics and sectorial harmonics are second order small quantities, the third body gravitational perturbation , the solar radiation pressure perturbation and the atmospheric drag perturbation are usually also second order small quantities. For a near-Earth satellite with a large mass ratio, the atmospheric drag perturbation can even reach the first order small quantity, these perturbation parts depend on the model subtraction, but if the orbiter is equipped with an accelerometer, the solar radiation pressure and atmospheric drag perturbation accelerations experienced by the orbiter in orbit can be measured to further improve the estimation accuracy. The non-conservative force acceleration measured by the accelerometer can be represented as F ( X , Y , Z ), where X is the acceleration component along the direction of motion of the orbiter, Y is the acceleration component along the normal direction of the orbital plane, Z is the acceleration component along the direction. The contributions of r and can be measured by the accelerometer and subtracted, which is given by:

[0056] where and ​respectively, are the precessions of the ascending node due to solar radiation pressure and atmospheric drag, f is the true anomaly. Since the method of the present application does not directly measure the orbit, the main methods currently available for determining the gravity field from the orbit plane precession are the Kaula linear perturbation method and the average perturbation method. The basic principle of these methods is to convert the perturbation of the satellite position into the perturbation of the six orbital elements, and the motion of the orbital elements can be expressed by the function of the disturbing potential, thereby establishing the relationship between the orbit perturbation and the disturbing potential.

[0057] Since the average perturbation method has a clearer physical meaning, the following formula will give the formula for obtaining J 2 and J 4 from the average perturbation method.

[0058] By constraining the orbit with a circular orbit, the perturbation caused by the gravity field coefficients mainly causes the precession of the ascending node. From the Lagrange perturbation motion equation:

[0059] where the left side of the equation is the first order derivative of the six orbital elements with respect to time. From the Kaula linear perturbation method, we can express the ascending node right ascension solution in formula (9) as:

[0060] where is the unperturbed average element; is the long-term term; is the long-period term; is the short-period term. Formula (12) divides the perturbation term in the motion equation into the following three categories: 1) the long-term term that does not contain the angular coordinate and is only related to the semi-major axis a , eccentricity e and inclination i ; 2) the long-period term related to the argument of perigee or the ascending node right ascension ; 3) the short-period term related to the satellite mean anomaly .

[0061] The perturbation caused by the gravity field is divided into long-term terms, long-period terms and short-period terms. These gravitational potentials cause different frequency components in the time series of the orbital elements, as shown in Fig. 8 and Fig. 9 , where the horizontal coordinate is the frequency and the vertical coordinate is the power spectral density of the ascending node right ascension. The signals below 10 -5 Hz are the long-term terms and their multiple frequency signals, that is ; the signals above 10 -5 Hz are the long-period term signals .

[0062] The perturbation caused by the Earth's non-spherical gravity, where the long-term term dominates, is given by The main long-term term component of the ascending node right ascension and the second order long-term term are given by

[0063]

[0064] where G is the gravitational constant; and since it is a circular orbit, the above equation can be simplified to G

[0065] ; Therefore, after we measure the ascending node right ascension precession by equation 6, we can filter out the long-term term signal Fig. 9 contained in it by spectral analysis, and then calculate the low order even zonal harmonic coefficients 2, J 4… of the gravity field by the above equation. Similarly, by solving the first order long period term J , we can calculate the low order odd zonal harmonic coefficients 3, J 5… of the gravity field; this step can also be used to solve the gravity field by the Kaula linear perturbation method analytically, avoiding the complex numerical integration, and thus more efficient. J The overall flow of Example 1 is summarized in

[0066] . Fig. 10 Example 2

[0067] Repeat Example 1 with the same steps, except that the orbit inclination is set to a high inclination (85 deg); and in step 2.2, a conical eclipse model is used: i when , ; when ,

[0068] where the duration τ , power spectral density, orbiter β angle, orbiter β angle power spectral density, sine of the ascending node right ascension precession sin( Ω ), and sin(​​Ω The simulation results of the power spectral density are also marked in different colors in Fig. 3 - Fig. 9 .

[0069] Those skilled in the art easily understand that the above description is only the preferred embodiment of the present application, and is not used to limit the present application, any modification, equivalent replacement and improvement made within the spirit and principle of the present application should be included in the protection scope of the present application.

Claims

1. A method for inverting a planetary gravitational field, used to obtain low-order band harmonic coefficients of the gravitational field, characterized in that, Includes the following steps: S1. Under circular orbit, obtain the duration of the planet's shadowed region within one orbital period using orbiter payload data. ;in, μ The gravitational constant of a planet, R P The average radius of the planet. h For the orbital height, θ The angle is the angular distance of the orbiter's movement. t For runtime; S2. Based on the duration of the shaded area τ (t) The angle between the Sun-Planet Vector and the Orbiter Plane is obtained using a cone-shaped or cylindrical ground shadow model. β (t), and thus obtain the right ascension precession of the ascending node. Ω (t); S3. Precession momentum from right ascension at the ascending node, determined through spectral analysis. Obtain long-term items and long-period terms ; S4. Based on the long-term item The harmonic coefficients of the low-order even-numbered terms of the gravitational field are obtained through linear perturbation or average perturbation methods; based on the long-period terms... The harmonic coefficients of the low-order odd terms of the gravitational field can be obtained by using the linear perturbation method or the analytical method.

2. The planetary gravitational field inversion method as described in claim 1, characterized in that, In step S1, the duration of the shadowed region of the planet within one orbital period is obtained. The method is STK simulation.

3. The planetary gravitational field inversion method as described in claim 1, characterized in that, The orbiter is a low-inclination orbiter, which is obtained in step S2 using a cone-shaped ground shadow model. β (t), the formula for the cone-shaped ground shadow model is: 。 4. The planetary gravitational field inversion method as described in claim 1, characterized in that, The orbiter is a high-inclination orbiter, which is obtained in step S2 using a cylindrical ground shadow model. β (t), the formula for the cylindrical ground shadow model is: in, R S The average radius of the sun. A This represents the average distance between the Sun and the planets.

5. The planetary gravitational field inversion method as described in claim 1, characterized in that, In step S2, the precession of the right ascension of the ascending node is obtained. Ω The formula for (t) is: ； in, ε It is the angle between the obliquity of the ecliptic and the obliquity of the sun. Г This is the true solar longitude of the zodiac. i This represents the track inclination angle.

6. The planetary gravitational field inversion method as described in claim 1, characterized in that, In step S3, the non-conservative force acceleration measured by the accelerometer is subtracted from the precession of solar radiation pressure and atmospheric drag on the right ascension of the ascending node. The impact.

7. The planetary gravitational field inversion method as described in claim 1, characterized in that, In step S4, the low-order even-numbered band harmonic coefficients of the gravitational field are obtained by the averaging perturbation method. J n , where n is an even number.

8. The planetary gravitational field inversion method as described in claim 7, characterized in that, With n = 2 and 4, we obtain the band harmonic coefficients of the low-order even-numbered terms of the gravitational field. J 2 and J Method 4 is as follows: ； in, G The gravitational constant is... i For the track inclination angle, M P Planetary mass, a It is the semi-major axis.

Images