A method for inverting planetary gravity fields
By measuring shadow time using orbiter payload data and combining it with Earth shadow models and spectral analysis, the problem of radio measurement errors in long-distance planetary exploration was solved, and high-precision planetary gravity field inversion was achieved.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- HUAZHONG UNIV OF SCI & TECH
- Filing Date
- 2025-08-07
- Publication Date
- 2026-07-31
AI Technical Summary
In long-distance planetary exploration, existing technologies using radio measurements to obtain planetary gravitational field information are subject to signal interference and insufficient model accuracy, making it difficult to determine the global gravitational field of a planet.
By using data from the orbiter's own payload, the duration of the planet's shadowed region during its orbital period is measured. Combined with a conical or cylindrical Earth shadow model, the angle between the Sun-planet vector and the orbiter plane is obtained. Using spectral analysis and perturbation methods, the low-order band harmonic coefficients of the gravitational field are retrieved.
It enables high-precision acquisition of planetary gravitational field information at long distances, reduces radio measurement errors, and improves the accuracy of gravitational field determination.
Smart Images
Figure CN120993509B_ABST
Abstract
Description
Technical Field
[0001] This application belongs to the interdisciplinary technical fields of satellite gravimetry, geodesy, orbital perturbation, and space science. Specifically, it relates to a method for inverting planetary gravitational fields. Background Technology
[0002] Planetary gravity field measurement is currently a hot topic and a frontier in deep space exploration, serving as a crucial means of understanding the composition and evolution of planets. Currently, planetary gravity fields are obtained through inversion of spacecraft orbital position information using deep space exploration technologies such as radio ranging, velocity measurement, and Very Long Baseline Interference (VLBI), as illustrated in patent document CN114740541A. The main advantages of this approach are high precision and suitability for long-distance missions within hundreds of millions of kilometers. However, radio signals are affected by interference along their propagation path (such as the atmosphere and solar plasma). Taking Mars as an example, its distance from Earth varies from 0.5 to 500 million kilometers. The inverse square law of distance causes a sharp decrease in signal strength with distance, deteriorating the accuracy of ranging and velocity measurements. Furthermore, the difficulty in achieving global coverage of radio observations and the poor accuracy of non-conservative planetary force models make determining the global gravity field of planets quite challenging. Therefore, directly measuring planetary gravity field information using the orbiter's own payload is of great significance. Summary of the Invention
[0003] To address the shortcomings of existing technologies, this application provides a novel method for planetary gravity field inversion, which relies solely on orbiter payload data without the aid of radio measurement results, aiming to solve the problem of radio measurement errors caused by long distances.
[0004] 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, m is the gravitational constant of the planet. R P The average radius of the planet. h For the orbital height, i The angle is the angular distance of the orbiter's operation, and t is the time of operation; S2. Based on the duration of the shaded area t (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. Oh (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.
[0005] Preferably, in step S1, the duration of the shadowed region of the planet within one orbital period is obtained. The method is STK simulation.
[0006] Preferably, the orbiter is a low-inclination orbiter, which is obtained in step S2 using a conical ground shadow model. β (t), the formula for the cone-shaped ground shadow model is: .
[0007] Preferably, 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:
[0008] in, R S The average radius of the sun. A This represents the average distance between the Sun and the planets.
[0009] Preferably, in step S2, the precession of the right ascension of the ascending node is obtained. Oh The formula for (t) is: ; in, e It is the angle between the obliquity of the ecliptic and the obliquity of the sun. G This is the true solar longitude of the zodiac. i This represents the track inclination angle.
[0010] Preferably, 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.
[0011] As a further preferred option, , ; in, The number of roots without photoperiod is the average root number. It is a short-period term. This is the precession of the right ascension of the ascending node caused by solar radiation pressure. This is the precession of the right ascension of the ascending node caused by atmospheric drag. r For planet-orbiter vectors, f For true near point angle, g The average angular velocity of the orbiter. oh The perigee argument, a For the semi-major axis, Y This refers to the acceleration component of the nonconservative force acceleration measured by the accelerometer along the normal direction of the orbital plane.
[0012] Preferably, 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.
[0013] As a further preferred option, n is 2 or 4, to obtain the low-order even-numbered terms of the gravitational field with harmonic coefficients. J 2 and J Method 4 is as follows:
[0014] ; in, G The gravitational constant is... i For the track inclination angle, M P Planetary mass, a It is the semi-major axis.
[0015] 1. Since the accuracy of the inverted gravity field is mainly affected by the noise level of the payload itself, as the measurement relies solely on the orbiter's own payload data, this application enables long-distance measurement. 2. Due to the use of shadow duration t (t) is a measurement technique, so it does not rely on radio observation; 3. Since the orbiter can carry accelerometers, it is preferable to use non-conservative force acceleration measured by accelerometers, after deducting the precession of solar radiation pressure and atmospheric drag on the right ascension of the ascending node. The impact; when using an accuracy of 10 -12 m / s 2 When an accelerometer of this magnitude is used as a solar-sensitive payload, it can measure non-conservative forces acting on the orbiter itself. This allows it to measure and subtract non-conservative force perturbations while simultaneously sensing the duration of the orbiter's shadow period. Therefore, it can measure and subtract the model at the same time, resulting in high accuracy. Attached Figure Description
[0016] Figure 1 This indicates the definition from the relevant perspective of this application, where fThe angle is the angle subtended by the Sun-Planet-Orbiter; β The angle is the angle between the Sun-Planet Vector and the orbital plane of the orbiter; i The angle is the angular distance of the orbiter's movement. S Planet-star vector, r This is the planet-orbiter vector.
[0017] Figure 2a Represents a cone-shaped shaded area model; Figure 2b Represents a cylindrical shaded area model; Figure 3 A schematic diagram illustrating the STK simulation of Example 1; Figure 4 This indicates the duration of the shaded area in Examples 1 and 2. t Simulation results; Figure 5 This refers to Example 1. t The power spectral density; Figure 6 Orbiters representing Embodiments 1 and 2 β Angle simulation results; Figure 7 Orbiters representing Embodiments 1 and 2 β Angular power spectral density; Figure 8 The sine value of the precession of the right ascension of the ascending node in Examples 1 and 2 is sin( Oh ); Figure 9 The sin(representation) of Examples 1 and 2 Oh Power spectral density; Figure 10 This describes the process of Embodiment 1 of this application. Detailed Implementation
[0018] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.
[0019] In the description of this application, it should be understood that the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined as "first" or "second" may explicitly or implicitly include one or more of that feature. In the description of this application, "multiple" means two or more, unless otherwise explicitly specified.
[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, m is the gravitational constant of the planet. R P The average radius of the planet. h For the orbital height, i 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 oh Right ascension of ascending node Oh 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: Figure 1 As shown; where, f The angle is the angle subtended by the Sun-Planet-Orbiter; βThe angle is the angle between the Sun-Planet Vector and the orbital plane of the orbiter; i The angle is the angular distance of the orbiter's movement. S Planet-star vector, r This is the planet-orbiter vector.
[0024] S2. Based on the duration of the shaded area t (t), obtain the angle between the Sun-Planet Vector and the orbiter plane. β (t); Figure 2 shows two shaded area models, respectively. Figure 2a Conical ground shadow model and Figure 2b Columnar ground shadow model.
[0025] When the orbiter is a high-inclination orbiter, it is obtained through a cylindrical ground shadow model. β (t), the formula for the cylindrical ground shadow model is:
[0026] in, R S The average radius of the sun. A This represents the average distance between the Sun and the planets. When the orbiter is a low-inclination orbiter It is always true, therefore it can be obtained through the cone-shaped ground shadow model. β (t), the formula for the cone-shaped ground shadow model is: .
[0027] Furthermore, through β (t) Obtain the precession of the right ascension of the ascending node Oh (t), specifically:
[0028] in, e It is the angle between the obliquity of the ecliptic and the obliquity of the sun. G This is the true solar longitude of the zodiac. i The inclination angle is the orbital inclination angle, both of which are known quantities.
[0029] S3. Precession momentum from right ascension at the ascending node, determined through spectral analysis. Obtain long-term items and long-period terms ; of which less than 10 -5 The Hz term represents the long-term term and its harmonics, which are... ; higher than 10 -5 A Hz signal is a long-period signal. ; , in, The number of roots without photoperiod is the average root number. It is a short-period term. This is the precession of the right ascension of the ascending node caused by solar radiation pressure. This refers to the precession of the right ascension of the ascending node caused by atmospheric drag. ; in, r For planet-orbiter vectors, f For true near point angle, g The average angular velocity of the orbiter. oh The perigee argument, a For the semi-major axis, e Eccentricity Y The acceleration component of the nonconservative force acceleration measured by the accelerometer is the acceleration component along the normal direction of the orbital plane; therefore, it can be measured and subtracted by the accelerometer. and The part.
[0030] S4. Obtain the low-order band harmonic coefficients of the gravitational field J n Where n is the order, which is a natural number; the higher the required precision, the higher the upper limit of n. Among them, the harmonic coefficients of the lower-order even-numbered terms of the gravitational field (i.e., n=2, 4, 6...) are determined according to the long-term term. The coefficients of the low-order odd-numbered terms of the gravitational field (i.e., n=3, 5, ...) are obtained through linear perturbation or average perturbation; the coefficients of the long-period terms of the gravitational field are obtained through long-period terms. It can be obtained using linear perturbation or analytical methods.
[0031] To obtain through the average perturbation method J 2 and J Taking 4 as an example, the harmonic coefficients of the low-order even-numbered terms of the gravitational field are obtained. J 2 and J Method 4 is as follows:
[0032] ; in, G The gravitational constant is... i For the track inclination angle, M P For the orbiter's level anterior angle, a For the semi-major axis, e Eccentricity; because e =0, the above expression can be further simplified to:
[0033] ; Other methods can be used to obtain the low-order band harmonic coefficients of other gravitational fields.
[0034] The following is an example.
[0035] Example 1 The planetary gravitational field inversion method based on orbiter shadowing time variations includes the following steps: Step 1: Measure the time the satellite spends in the shadow zone for each orbit using orbiter payload data.
[0036] like Figure 2a As shown, when using the cone shading model, the orbiter enters three phases within one orbital period: the illuminated region, the penumbra, and the umbra. We will utilize the orbiter's solar-sensitive payload data to record the time when the orbiter enters each phase. t ,in t 1 This indicates the time from the center of the illuminated area to entering the first penumbra. t 2 Indicates the time from the center of the illuminated area to entering the total shadow area. t 3 This indicates the time from the center of the illuminated area to entering the second penumbra. t 4 This indicates the time from the center of the illuminated area to the next illuminated area. R e is the radius of the planet. R 's' represents the radius of the star opposite the planet. For solar-sensitive payloads, there will be a noticeable jump signal at the moment of entering and exiting the shadow region. For example, an electrostatic accelerometer measures the non-conservative forces experienced by an orbiter in orbit. As the orbiter moves from the illuminated region into the penumbra, the solar radiation pressure it experiences gradually decreases, reaching zero upon entering the total shadow. Therefore, we need to record the moments when the payload measurement begins to decrease (increase) and stops decreasing (increase). This allows us to calculate the duration of the planet's entire shadow period during one orbital period. t A And time in the entire film zone t U :
[0037]
[0038] This application selects a circular track model for calculation. Under the circular track, the tracker moves at a constant speed. .in m is the gravitational constant of the planet.R P The average radius of the planet. h For the orbital height, i The angle represents the angular distance of the orbiter's movement, and t represents the time of travel. Since a circular orbit is chosen, it is assumed that... h It remains unchanged.
[0039] like Figure 3 As shown, we used STK aerospace simulation to simulate Mars as the celestial body. The pink dot in the figure represents the spatial position of the orbiter, which orbits in a circular orbit (i.e., with an eccentricity of 1 / 200°). e =0) revolves around Mars, ensuring that the equations (3) and (4) are correct. h Relatively stable; according to formulas (11) and (12), the lower the orbital altitude, the greater the orbital perturbation it experiences. Since the average radius of Mars is 3389.5 km, and the lowest orbital altitude is approximately 300 km, we selected a semi-major axis... a The orbital inclination is 3700 km. According to formula (11), the smaller the orbital inclination, the greater the corresponding change in precession, which means that the Earth's gravity will cause a greater orbital perturbation, which is beneficial to the calculation in this application. However, if the inclination is too small, the change in perturbation will be smaller, which will affect the accuracy of the estimation. Therefore, the orbital inclination is important. i Set to a low tilt angle (15 degrees); similarly, it can be seen from formulas (11) and (12) that the right ascension of the ascending node is... Oh Precession and perigee angle oh , near point angle M and the right ascension of the ascending node Oh The initial value is irrelevant. For ease of calculation, in this embodiment, we selected 0deg as the initial value. The simulation parameter settings are shown in Table 1, and the simulation results are as follows. Figure 3 As shown.
[0040] Table 1 Simulation Parameters
[0041] Simulations show that the orbiter spends a certain amount of time in the shadow zone throughout the year. Changes such as Figure 4 As shown.
[0042] Step 2: Estimate the right ascension precession of the orbiter's ascending node based on the time in the shaded area.
[0043] 2.1 Selection of Shadow Model For the cone-shaped ground shadow model, based on geometric relationships... and With the orbiter track β 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... Figure 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 Figure 4 The time variation of the orbiter in the shaded area, measured under the simulation conditions of Example 1, is given. t (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. t Variation with simulation time t (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 ( t =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 Figure 5 The simulation is given in the middle. t The power spectral density curve of (t) is plotted with frequency on the x-axis and duration on the y-axis. t 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, and its expression as a function of time is as follows:
[0049] in, e It is the angle between the obliquity of the ecliptic and the obliquity of the sun. G This is the true solar longitude of the zodiac. i This represents the track inclination angle.
[0050] The influence can be derived from the above formula. β There are two main factors that cause the angle to change, namely 𝛤 Changes (caused by planets orbiting the sun) and Oh Changes (caused by orbital perturbations). Therefore, we can... β Extracting the right ascension precession of the ascending node from angle (t) . Figure 6 and Figure 7 Simulations are given respectively. β The time-domain variation of the angle and the power spectral density curve, with frequency on the x-axis and frequency on the y-axis. β The power spectral density of the angle, its frequency domain exhibits the same characteristics as t Similar frequency components, but due to limitations in measurement methods, through t Received β Since it can only be sampled once per orbital period, the orbital frequency signal component is missing. Therefore, this application is suitable for observing long-term and long-period signals.
[0051] Step 3: Determine the gravitational field coefficients based on the precession of the right ascension of the ascending node.
[0052] For the perturbed two-body problem, the equations of motion of the orbiter can be expressed as:
[0053] in The distance from the orbiter to the center of the planet. Accelerate the orbiter. The orbiter is accelerated by the gravitational pull of the planet's center. This is the perturbation acceleration acting on the orbit, which can be further expressed as follows, based on the main perturbation sources of terrestrial planets:
[0054] In the above formula degree Represents the gradient. Represents the gravitational perturbation function of a planet that is not spherical. Represents the gravitational perturbation function of the third body. Represents the acceleration due to solar radiation perturbation. Represents the atmospheric drag perturbation acceleration, where It can be represented as:
[0055] in and This is the gravitational field coefficient. For harmonic terms, For the field, m For the number of times (order). n For degree; P nm To associate Legendre functions, P n For Legendre functions, f The latitude of the Earth's core. l Longitude l nm This is the phase angle. For terrestrial planets, the largest perturbation factor is... Harmonic terms (first-order minors), other harmonic terms and Tian harmonic terms are second-order minors, third-body gravitational perturbation. Solar radiation perturbation and atmospheric drag perturbation This is usually also considered a second-order minor quantity. For near-Earth satellites with a large surface-to-mass ratio, atmospheric drag perturbations can even reach the level of first-order minor quantities. These perturbations partially depend on model subtractions, but if the orbiter is equipped with an accelerometer, it can measure the solar radiation pressure experienced by the orbiter in orbit. and atmospheric drag Perturbation acceleration can further improve estimation accuracy. The non-conservative force acceleration measured by the accelerometer is expressed as... F ( X , Y , Z ) ,in X The acceleration component is along the direction of the orbiter's motion. Y The acceleration component is along the normal direction of the orbital plane. Z For along r The acceleration component in the direction of acceleration. This can be measured and subtracted using an accelerometer. and The contribution is given by the following formula:
[0056] in and These are the precession of the right ascension of the ascending node caused by solar radiation pressure and atmospheric drag, respectively. f The true anomaly angle. Since the method in this application does not directly measure the orbit, the main methods currently available for determining the gravitational field through orbital plane precession include the Kaula linear perturbation method and the average perturbation method. Their basic principle is to convert the perturbation of the satellite position into the perturbation of six orbital elements, and the movement of the orbital elements can be represented by a function of the perturbation potential, thereby establishing the relationship between the orbital perturbation and the perturbation potential.
[0057] Since the physical meaning of the average perturbation method is relatively clear, the following formulas will be derived from the average perturbation method. J 2 and J Formula 4.
[0058] By constraining the orbiter with a circular orbit, the perturbation caused by the gravitational field coefficient mainly causes the precession of the right ascension of the ascending node, according to the Lagrange perturbation equations:
[0059] The left side of the equation represents the six orbital elements. Taking the first derivative with respect to time. By Kaula's linear perturbation method, we can express the solution for the right ascension of the ascending node in equation (9) as:
[0060] in The number of roots without photoperiod; It is a long-term term; It is a long-period term; For short-period terms. Formula (12) divides the perturbation terms in the equation of motion into the following three categories: 1) without angular coordinates, only related to the semi-major axis. a eccentricity e and tilt angle i The relevant long-term terms; 2) with the perigee argument or ascending node right ascension Related long-period terms; 3) Angle of near the satellite The relevant short-period terms.
[0061] Gravitational perturbations are categorized into long-term, long-period, and short-period terms. These gravitational potentials cause different frequency components in the orbital element time series, such as... Figure 8 and Figure 9 As shown, the horizontal axis represents frequency, and the vertical axis represents the power spectral density at the right ascension of the ascending node. Spectral density values below 10 are... -5 The Hz term represents the long-term term and its harmonics, which are... ; higher than 10 -5 A Hz signal is a long-period signal. .
[0062] The perturbations caused by Earth's non-spherical gravity, with the long-term term accounting for the majority, are given below at the right ascension of the ascending node. The main long-term components of the first-order long-term solution and second-order long-term terms :
[0063]
[0064] in, G Let be the gravitational constant; since it is a circular orbit, the above equation can be simplified to:
[0065] ; Therefore, we obtain the precession of the right ascension of the ascending node using Equation 6. Then, through spectrum analysis (such as...) Figure 9 As shown), filter out the long-term term signals contained therein. The harmonic coefficients of the lower-order even-numbered terms of the gravitational field can then be calculated using the formula above. J 2. J 4… Similarly, by solving the first-order long-period terms… The band harmonic coefficients of the low-order odd-numbered terms of the gravitational field can be calculated. J 3. J 5…; This step can also utilize the Kaula linear perturbation method to analytically process the perturbation and solve for the gravitational field, avoiding complex numerical integration and achieving higher computational efficiency.
[0066] The overall process of Example 1 is summarized as follows: Figure 10 As shown.
[0067] Example 2 Repeat Example 1 using the same steps, except that the track inclination angle is changed. i Set to a high tilt angle (85 degrees); in step 2.2, use a conical ground shadow model: when hour, ; when hour
[0068] Among them, duration t Power spectral density, orbiter β Angle, orbiter β Angular power spectral density, sinusoidal value of the precession of the right ascension of the ascending node (sin( Oh ), and sin( Oh The simulation results of the power spectral density are also indicated by different colors. Figure 3-Figure 9 I won.
[0069] Those skilled in the art will readily understand that the above description is merely a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this application should be included within the protection scope of this application.
Claims
1. A planetary gravity field inversion method for retrieving low-degree band harmonic coefficients of a gravity field, characterized in that, Includes the following steps: S1. Under the circular orbit, the duration of the shadow zone of the planet in an orbital period is obtained by the orbital data of the orbiter payload ; wherein, μ Gp is the gravitational constant of the planet, R P Rp is the average radius of the planet, h H is the orbital height, θ Θ is the angular distance of the orbiter's run, t t is the time of the run; 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.