A method for recovering the harmonic gravity field of Phobos' ellipsoid based on Doppler orbit data

Abstract

The present invention discloses a method for recovering the ellipsoidal harmonic gravity field of Phobos based on Doppler orbit data. For the ellipsoidal small celestial body Phobos, this method selects the ellipsoidal harmonic form with stronger applicability instead of the traditional spherical harmonic form. By processing the two-way velocity measurement data of the detector, the partial derivatives of the two-way Doppler velocity measurement value to the detector position and velocity vector at any time, and the partial derivatives of the perturbation acceleration to the detector position, velocity vector and the parameter to be estimated at any time are solved; the solved partial derivatives are formed into a design matrix; the single arc segment normal equation is determined; the multi-arc segment normal equation is integrated to solve the correction value of the parameter to be estimated, and the iteration is repeated until convergence, and finally the accurate solution value of the parameter to be estimated is obtained. Compared with the traditional spherical harmonic coefficients, the ellipsoidal harmonic coefficients can extract the gravity field information of the ellipsoidal small celestial body Phobos from the Doppler orbit data with faster convergence speed, fewer coefficients and higher accuracy.

Description

Technical Field

[0001] The invention belongs to the field of deep space exploration and computing technology, and in particular relates to a method for recovering the harmonic gravity field of the Phobos ellipsoid based on Doppler orbit data. Background Art

[0002] Phobos is the larger and closer of the two natural satellites of Mars, orbiting Mars at an altitude of 6000 km. However, its future is uncertain, and it is expected to collide with the surface of Mars or disintegrate into dust in the next 30 to 50 million years. This ellipsoidal object without an atmosphere is one of the darkest objects in the Solar System, with an albedo of only 0.071. These unique physical properties make it a focus of scientific research. Mars exploration missions have been carried out since the 1960s and have achieved many milestones. However, exploration of Phobos is relatively rare, as it is impossible to orbit Phobos in a stable manner due to the Hill sphere located inside Phobos. Therefore, exploration of Phobos is mainly carried out by flyby or co-orbiting with it in a quasi-satellite orbit or QSO (Quasi-Satellite Orbit). To date, all Mars exploration missions have conducted only one flyby of Phobos (in a hyperbolic orbit around Phobos). Recent work on the determination of the mass of Phobos and its gravity field (spherical harmonic coefficients up to second order) was based on data from Mars Express. Considering the safety of the spacecraft, the probe's flyby was carried out at a relatively large distance from Phobos, resulting in only a small deflection of the spacecraft's orbit caused by Phobos' gravity. Therefore, the characteristics of the trajectory deflection in the Doppler signal are close to the measurement noise.

[0003] Therefore, in this case, it is crucial to deploy data processing tools that can extract the full content of the signal. Phobos is not a round body with a diameter of 20 km, but a triaxial ellipsoid of approximately 27 km×22 km×18 km. Since spherical harmonics are more suitable for expressing the gravity field of nearly spherical objects, the gravity field of the ellipsoidal Phobos is not suitable for modeling using spherical harmonics.

[0004] Currently, there is still a lack of a method and application for recovering the ellipsoidal harmonic gravity field of Phobos based on the ellipsoidal harmonic model. Summary of the invention

[0005] In view of the above technical problems, the present invention aims to overcome the traditional spherical harmonic function expression of gravity field, and proposes a method for recovering the spherical harmonic gravity field of Phobos based on Doppler orbit data, which is closer to an ellipsoid than a sphere. This method can represent the gravity field of an ellipsoidal asteroid with fewer coefficients and smaller errors.

[0006] To achieve the above objectives, the technical solution of the present invention provides the following method:

[0007] A method for recovering the harmonic gravity field of the Phobos ellipsoid based on Doppler orbit data, the steps are as follows:

[0008] Step 1: Establish the Doppler two-way velocity observation equation, use the detector Doppler orbit data to measure the observation value, and obtain the two-way Doppler velocity value; solve the partial derivatives of the two-way Doppler velocity value with respect to the detector position and velocity vector at any time, and the partial derivatives of the perturbation acceleration with respect to the detector position, velocity vector and the parameters to be estimated at any time; form a design matrix with the solved partial derivatives; wherein the parameters to be estimated include the GM term and the ellipsoidal harmonic gravity potential coefficient ;

[0009] Step 2: Determine the single arc segment normal equation based on the partial derivative in step 1;

[0010] Step 3: Fusion of multi-segment normal equations to solve the estimated parameter correction values The initial value is modified and used as the new initial value to iterate and repeat steps 1 and 2 until convergence to obtain the accurate solution value of the parameter to be estimated.

[0011] The present invention also provides an electronic device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the program, the method for recovering the harmonic gravity field of the Phobos ellipsoid based on Doppler orbit data as described above is implemented.

[0012] Beneficial effects of the present invention:

[0013] The method of the present invention can more accurately restore the ellipsoidal harmonic gravity field of Phobos, taking into full account the characteristics of its ellipsoidal shape, and can provide a more accurate gravity field expression compared to the traditional spherical harmonic function form. This is of great significance for the gravity field research of Phobos and related scientific exploration.

[0014] The present invention is based on the two-way velocity measurement mode. In the Phobos exploration mission, the two-way velocity measurement observation value is used to achieve the precise orbit determination of the detector and the inversion solution of the ellipsoidal harmonic gravity field. Compared with the spherical harmonic gravity field, the ellipsoidal harmonic gravity field obtained by the present invention can be represented by fewer coefficients, has a smaller convergence area and a faster convergence speed. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 It is a schematic diagram of the dual-pass Doppler velocity measurement mode;

[0016] Figure 2 Schematic diagram of the convergence domain of harmonic functions. DETAILED DESCRIPTION

[0017] In order to better understand the technical solution of the present invention, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments, but the content of the present invention is not limited thereto.

[0018] The present invention proposes to express the gravity field of the ellipsoidal small celestial body Phobos in ellipsoidal harmonic form. Specifically, by processing the detector's two-way velocity measurement Doppler orbit tracking data, the detector position is precisely determined and the ellipsoidal harmonic form gravity field is solved, thereby achieving an improvement in the accuracy of gravity field solution.

[0019] The present invention aims at small celestial bodies of elongated or ellipsoidal shape, and uses ellipsoidal harmonics with fewer coefficients, smaller convergence area and faster convergence speed to represent their gravity fields. The convergence domain diagrams of the two harmonic functions are shown in Figure 2 .

[0020] The method provided by the present invention is as follows:

[0021] Step 1, establish the Doppler two-way velocity observation equation, use the detector Doppler orbit data to measure the velocity observation value, and obtain the two-way Doppler velocity value; solve the partial derivatives of the two-way Doppler velocity value to the detector position and velocity vector at any time, and the partial derivatives of the perturbation acceleration to the detector position, velocity vector and the parameters to be estimated at any time; form a design matrix with the solved partial derivatives; wherein the parameters to be estimated include the GM term and the ellipsoidal harmonic gravity potential coefficient .

[0022] Specific methods include:

[0023] The Doppler two-way velocity measurement observation equation is as follows:

[0024] (1)

[0025] in, , , , They are all one-way distance measurement values ​​of different links. and are the two-way ranging values ​​at the start and end of the integration period, respectively. is the two-way speed measurement value, is the Doppler frequency shift integration period, and the schematic diagram shows Figure 1 .

[0026] The two-way speed measurement observation value Right now ,right Satellite status at all times Right now The partial derivative of for:

[0027] (2)

[0028] (3)

[0029] , , They are The Cartesian coordinates of the satellite's position at the moment, , , They are Satellites at x , y , z Velocity components in three directions.

[0030] The two-way distance measurement value The partial derivative of the satellite state at time t is:

[0031] (4)

[0032] (5)

[0033] (6)

[0034] (7)

[0035] in, The value is or . Indicates the time when the ground observation station transmits the uplink signal. At this time, the ground station position vector is ; Indicates the moment when the detector receives the uplink signal. At this time, the detector position vector is ; Indicates the moment when the ground station receives the downlink signal. At this time, the ground station position vector is ; Indicates the uplink distance; Indicates the downlink distance.

[0036] By establishing various precise perturbation force models and solving various perturbation accelerations, the numerical integration method can be used to calculate the initial acceleration of the detector. Satellite status at all times The integral solution is obtained at any time Satellite status at . Doppler data observation value at any time With theoretical value The residual between . Then we have:

[0037] (8)

[0038] Theoretical value The perturbation model can be expressed as the satellite state of the detector at any time and time function, so in The theoretical value in formula (2) Taylor first-order linear expansion gives:

[0039] (9)

[0040] in, is the two-way speed measurement observation value Right now For any time The partial derivative of the orbital vector at any time is obtained by observing the value at any time The partial derivative of the orbital vector is denoted by , which can be expressed as:

[0041] (10)

[0042] in, For any time The partial derivative matrix of the detector state vector at time t with respect to the state vector at the initial time and the parameters to be estimated.

[0043] For a time series observation, formula (9) can be expressed as:

[0044] (11)

[0045] in, y t for t The residual of the observation at time t, for t The residual of the moment;

[0046] The time series equation for all moments Combined together, it can be expressed as the Jacobian matrix The design matrix That is:

[0047] (12)

[0048] in, for Satellite position vector at time, for Satellite position vector at time, for Satellite velocity vector at time, are the parameters to be estimated (specifically the GM term and the ellipsoidal harmonic gravity potential coefficient ); for Satellite velocity vector at this moment; for Satellite velocity vector at time.

[0049] The first-order derivative of the design matrix is:

[0050] (13)

[0051] in, n p is the number of parameters to be estimated; d Expressed as a differential, dadr represents the partial derivative of the perturbation acceleration with respect to the satellite position vector, dadv represents the partial derivative of the perturbation acceleration with respect to the velocity vector, dadp represents the partial derivative of the perturbation acceleration with respect to the estimated parameter (including the perturbation acceleration with respect to the GM term and the ellipsoidal harmonic gravity potential coefficient partial derivatives of ).

[0052] The second-order derivative of the design matrix is:

[0053] (14)

[0054] for Satellite position vector at time For initial Satellite position vector at time The partial derivative of for Satellite position vector at time Parameters to be estimated The partial derivative of

[0055] By integrating the variational equations (13) and (14), we can obtain the detector position vector at any time.

[0056] Ellipsoid Harmonic Derivatives , and The calculation method is as follows. Usually the ellipsoidal harmonic gravity field can be expanded as:

[0057] (15)

[0058] in, is the gravity potential at the calculation point, is the gravitational constant, is the mass of the central celestial body, n is the order, and the superscript p is the degree; is the ellipsoidal harmonic gravity potential coefficient when the order is n and the degree is p; is the second kind of Lamé function with order n and degree p; is the first type of Lame function with order n and degree p. The Lame function is divided into four categories: K, L, M and N. , , is the ellipsoid coordinate of the calculation point, is the maximum semi-axis length of the reference ellipsoid.

[0059] It should be noted that when p is a superscript, it indicates the number of times.

[0060] In order to obtain partial derivatives conveniently, formula (15) is written as follows:

[0061] (16)

[0062] in, The sub-items can be obtained by integration:

[0063] (17)

[0064] Among them, the focal length ,focal length , a , b , c are the three semi-axis lengths of the reference ellipsoid, t express t node, n Is the order.

[0065] express t At the node, the order is n , the number of times is p， Ellipsoid coordinates Take The normalized first-kind Lamé function is calculated as: ; express t At the node, the order is n , the number of times is p The term function of express t At the node, the order is n , the number of times is p The partial function of , within this range t Samples are taken at even intervals. , and According to the sub-item , and By λ Replace with t Calculated.

[0066] Itemization ( i =1, 2, 3) is calculated according to the following formula:

[0067] (18)

[0068] in, is the coefficient of the Lamé function polynomial, i is 1, 2, 3; the values ​​of m are shown in Table 1.

[0069] Itemization It can be calculated according to Table 1:

[0070] Table 1 Different categories Expression, item

[0071]

[0072] Polynomial coefficients The value varies depending on the classification of the Lamé function, and is usually solved by a recursive relationship:

[0073] (19)

[0074] Among them, the parameters , and The values ​​can be obtained according to Tables 2 to 5, where p is the eigenvalue of the corresponding matrix after formula (19) is written in matrix form.

[0075] To simplify the formula, let , , .

[0076] Table 2 K-type Lamé function polynomial coefficients solution recursive relationship ( )

[0077] Table 3. Recursive relationship for solving polynomial coefficients of L-type Lamé functions ( )

[0078]

[0079] Table 4. Recursive relationship for solving polynomial coefficients of M-type Lamé functions ( )

[0080]

[0081] Table 5 Recursive relations for solving polynomial coefficients of N types of Lamé functions ( )

[0082]

[0083] According to the sub-items in Table 2-5 , and The coefficients obtained from the recurrence relation are written in matrix form as follows:

[0084] (20)

[0085] The above formula satisfies: , so the problem of solving the coefficients of the Lamé function polynomial can be simplified to the problem of solving the eigenvalues ​​and eigenvectors of the matrix. To facilitate the solution of the eigenvalues and the eigenvector , is the constituent element of the eigenvector (same as in formula 19 b i ), the present invention transforms the matrix in formula 20 to obtain , written in the form of a tridiagonal matrix, that is:

[0086] (twenty one)

[0087] Among them, S is the transformed tridiagonal matrix, D is the diagonal matrix related to the matrix transformation, represents the inverse matrix of matrix D, , , , is the parameter item after matrix transformation (the subscript can be 0, 1, ..., r), and the eigenvector corresponding to the matrix is V t Therefore, the eigenvalues ​​of the two matrices before and after the transformation are the same, and the eigenvectors of the original matrix before the required transformation can be obtained according to Calculated.

[0088] Itemization It can be calculated according to Table 6.

[0089] Table 6: Items in different categories Expression

[0090]

[0091] in, , , The Cartesian coordinates of the calculation point.

[0092] The first-order derivative of the gravitational potential with respect to time is the perturbation acceleration caused by the gravitational field of the celestial body. for:

[0093] (twenty two)

[0094] in, Perturbation acceleration , P is the rectangular coordinate With ellipsoidal coordinates The partial derivative matrix P and the gravity potential on the ellipsoidal coordinates The desired partial derivative The specific solution formula of is shown in the supplementary content.

[0095] The second-order derivative of the gravitational potential with respect to time, that is, the partial derivative of the perturbation acceleration with respect to the position of the detector at any time (Right now ) is:

[0096]

[0097] in, is the item function, ; i Represents the ellipsoid coordinate subscript, with values ​​of 1, 2, 3; , , and All are item functions; The order is n , the number of times is p The second kind of Lamé function when Indicates the maximum semi-axis length of the reference ellipsoid;

[0098] It is divided into four types of Lamé functions; , They are all partial functions, which can be calculated by formula 25, and the only difference is the subscript; and Represents the rectangular coordinates of the calculation point, i takes values ​​of 1, 2, 3; j takes values ​​of 1, 2, 3.

[0099] Itemized Function It can be calculated according to the following formula:

[0100] (twenty four)

[0101] in, To calculate the partial derivatives of the point ellipsoidal coordinates with respect to the rectangular coordinates, i takes values ​​of 1, 2, 3; m takes values ​​of 1, 2, 3.

[0102] The specific steps for solving this second-order partial derivative are detailed in the supplementary instructions.

[0103] Perturbation acceleration is one of the parameters to be estimated, the ellipsoidal harmonic gravity potential coefficient The partial derivative of for:

[0104] (25)

[0105] in, is the perturbation acceleration, is the N-body perturbation force, The order is n , the number of times is p The second kind of Lamé function when represents the maximum semi-axis length of the reference ellipsoid, is the item function Calculate the rectangular coordinates of the point The specific solution formula is shown in the supplementary instructions.

[0106] The partial derivative of the perturbation acceleration with respect to the GM term of one of the parameters to be estimated is for:

[0107] (26)

[0108] in, represents the partial derivative of the gravity potential V with respect to the rectangular coordinates of the calculation point, i represents the subscript of the rectangular coordinates of the calculation point, and its value is 1, 2, 3; is the gravitational constant, is the mass of the central object.

[0109] It should be noted that when combining the method equations, it is necessary to pay attention to the order in which the estimated parameters are stored and the order in which the estimated parameters are stored. Keep the storage order consistent. For example: dadp First store the partial derivative of the perturbation acceleration to the GM term and then store the partial derivative of the perturbation acceleration to the ellipsoidal harmonic gravity potential coefficient. When storing, you need to store the GM term prior value first, and then store the ellipsoidal harmonic gravity potential coefficient.

[0110] Step 2, determine the single arc segment normal equation according to the partial derivative in step 1;

[0111] According to the linear unbiased minimum variance, the normal equation can be determined as:

[0112] (27)

[0113] in, Represents the observed value at any time The partial derivative of the orbital vector when , for The transposed matrix of is the weight matrix, representing the accuracy difference between observations. , represents the noise of two-way speed measurement observations; is the correction value of the parameter to be estimated, is the prior value of the parameter to be estimated, is the residual of the observation at time t, is the prior error matrix of the parameters to be estimated, generally written in the form of a diagonal matrix:

[0114] (28)

[0115] Substitute the partial derivatives obtained in step 1 into formula (27) and replace the parameters to be estimated Local Parameters and global parameters , formula (27) can be written as a single arc segment normal equation:

[0116] (29)

[0117] in, is a local parameter The normal equation matrix, For global parameters The normal equation matrix, and are transposed to each other, representing the connection matrix of local parameters and global parameters respectively, z 1. z 2 corresponds to the part on the right side of the equal sign in formula 27, and represents the constituent elements on the right side of the equal sign of the single arc segment normal equation in formula 29.

[0118] Step 3: Fusion of multi-segment normal equations to solve the correction values ​​of the parameters to be estimated The initial value is corrected and used as the new initial value to iterate and repeat step 1 and step 2 until the relative error of the root mean square of the residuals of two adjacent iterations is less than 2%, and the accurate solution value of the estimated parameter is obtained.

[0119] The normal equation after multi-arc segment fusion is:

[0120] (30)

[0121] The superscripts represent the numbers of the merged arcs. z 1 is the right side of the equal sign of the local parameter solution formula corresponding to the linear unbiased least squares method , represents the local parameter prior value, z 2 is the right side of the equal sign of the global parameter solution formula corresponding to the linear unbiased least squares method , The superscript n indicates different arc segments. Represents the global parameter prior value.

[0122] Additional notes:

[0123] 1. The components of the partial derivative matrix P of the calculated point ellipsoid coordinates to the Cartesian rectangular coordinates can be solved according to the following formula:

[0124] (A.1)

[0125] (A.2)

[0126] (A.3)

[0127] (A.4)

[0128] in, To calculate the partial derivatives of the ellipsoidal coordinates of a point with respect to the rectangular coordinates, i=1, 2, 3; j=1, 2, 3.

[0129] 2. To facilitate the calculation of partial derivatives, the ellipsoidal harmonic expansion can be written as:

[0130] (A.5)

[0131] in,

[0132] ,

[0133] So the gravitational potential V has an influence on the rectangular coordinates of the calculation point. The partial derivatives of can be obtained through the ellipsoidal coordinates Partial derivatives with respect to rectangular coordinates To express:

[0134] (A.6)

[0135] i The possible values ​​are 1, 2, or 3.

[0136] Itemized Function Calculate the rectangular coordinates of the point The partial derivative of is:

[0137] (A.7)

[0138] (A.8)

[0139] Itemized Function Ellipsoid coordinates of the calculated point The first-order partial derivative of is:

[0140] (A.9)

[0141] in, , , i The value is 1, 2, 3. .

[0142] Itemized Function Ellipsoid coordinates of the calculated point The first-order partial derivative of is:

[0143] (A.10)

[0144] Substituting the above partial derivative into formula A.6, we can get The first-order partial derivative of the axial gravity potential (i.e., acceleration ) is:

[0145] (A.11)

[0146] Itemized Function Calculate the rectangular coordinates of the point The partial derivative of As shown in Table A.1:

[0147] Table A.1 Items of four types of Lamé functions Partial derivatives with respect to Cartesian coordinates

[0148]

[0149] 3. Gravity edge Axis pair The second-order partial derivative of )for:

[0150] (A.12)

[0151] The sub-item function Calculate the rectangular coordinates of the point The partial derivative of It can be expressed as:

[0152] (A.13)

[0153] Itemized Function Second-order partial derivatives of the rectangular coordinates of the calculation point It can be solved from Table A.2-A.5.

[0154] Table A.2 K-type Lamé function sub-items Second-order partial derivatives of the Cartesian coordinates of the calculation point

[0155]

[0156] Table A.3L-type Lamé function sub-items Second-order partial derivatives of the Cartesian coordinates of the calculation point

[0157]

[0158] Table A.4 Class M Lamé functions Second-order partial derivatives of the Cartesian coordinates of the calculation point

[0159]

[0160] Table A.5 N-type Lamé function sub-items Second-order partial derivatives of the Cartesian coordinates of the calculation point

[0161]

[0162] The method provided by the present invention realizes that for small celestial bodies of slender or ellipsoidal shape, the gravity field thereof is represented by ellipsoidal harmonics with fewer coefficients, smaller convergence area and faster convergence speed, which is different from the traditional spherical harmonic expansion form to express the gravity field. The present invention solves the first-order Phobos ellipsoidal harmonic gravity field model by processing the radio data of the ESA Mars Express orbit. Compared with the spherical harmonic gravity field model, the solution accuracy of the Phobos GM term is improved by 3-4 times.

[0163] The present invention also provides an electronic device, including a readable storage medium, on which a computer program is stored. When the computer program is executed, a method for recovering the harmonic gravity field of the Phobos ellipsoid based on Doppler orbit data as described above is implemented.

[0164] In specific implementation, the method proposed in the technical solution of the present invention can be implemented by technical personnel in this field using computer software technology to realize automatic operation process. System devices for implementing the method, such as computer-readable storage media storing the corresponding computer program of the technical solution of the present invention and computer equipment running the corresponding computer program, should also be within the scope of protection of the present invention.

[0165] The above description is only a preferred specific implementation manner of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications, equivalent substitutions and improvements made by any technician familiar with the technical field within the technical scope disclosed by the present invention should be included in the scope of protection of the invention.

Claims

1. A method for recovering the harmonic gravity field of the Phobos ellipsoid based on Doppler orbit data, characterized in that: Here are the steps: Step 1: Establish the Doppler two-way velocity observation equation, use the detector Doppler orbit data to measure the velocity observation value, and obtain the two-way Doppler velocity value; solve the partial derivatives of the two-way Doppler velocity value to the detector position and velocity vector at any time, and the partial derivatives of the perturbation acceleration to the detector position, velocity vector and the estimated parameters at any time; form a design matrix with the solved partial derivatives; wherein the estimated parameters include the GM term and the ellipsoidal harmonic gravity potential coefficient ; Step 2: Determine the single arc segment normal equation based on the partial derivative in step 1; Step 3: Fusion of multi-segment normal equations to solve the estimated parameter correction values The initial value is corrected and used as the new initial value to iterate and repeat steps 1 and 2 for iterative calculation until convergence to obtain the accurate solution value of the parameter to be estimated.

2. The method for recovering the harmonic gravity field of the Phobos ellipsoid based on Doppler orbit data according to claim 1, characterized in that: In step 1, the two-way Doppler velocity value is calculated by counting the time interval t c The integrated Doppler counts within the range are obtained, and the measured rate of change of distance is modeled as the difference between the round-trip distances at the end and starting points of the counting time.

3. The method for recovering the harmonic gravity field of the Phobos ellipsoid based on Doppler orbit data according to claim 2, characterized in that: The calculation method of the Doppler two-way velocity observation equation is as follows: Ground station The detector transmits uplink signals at all times. Receive and send downlink signals to the ground station at all times. Successfully receive downlink signals at all times; to The signal transmission process of the moment is the same; a single speed measurement model iterates the light travel time; is the time it takes for the first uplink signal to be transmitted from the ground station to the detector, It is the time for the first downlink signal to be transmitted from the detector to the ground station; is the time it takes for the first uplink signal to be transmitted from the ground to the detector, It is the time for the first downlink signal to be transmitted from the detector to the ground station; The Doppler two-way velocity measurement observation equation is as follows: in, , , , They are all one-way distance measurement values ​​of different links. and are the two-way ranging values ​​at the start and end of the integration period, respectively. is the two-way speed measurement observation value, is the Doppler frequency shift integration period; The two-way speed measurement observation value Right now ,right Satellite status at all times Right now The partial derivative of for: , , They are The Cartesian coordinates of the satellite's position at the moment, , , They are Satellites at x , y , z Velocity components in three directions; The two-way distance measurement value right The partial derivative of the satellite state at time is: in, The value is or ; Indicates the time when the ground observation station transmits the uplink signal. At this time, the ground station position vector is ; Indicates the moment when the detector receives the uplink signal. At this time, the detector position vector is ; Indicates the moment when the ground station receives the downlink signal. At this time, the ground station position vector is ; Indicates the uplink distance; Indicates the downlink distance.

4. The method for recovering the harmonic gravity field of the Phobos ellipsoid based on Doppler orbit data according to claim 1, characterized in that: In step 1, the method for calculating the partial derivative of the perturbation acceleration with respect to the detector position at any time is as follows: Taking the first-order derivative of the ellipsoidal harmonic gravity field formula with respect to time, we can obtain the perturbation acceleration caused by the celestial gravity field. ; P is the rectangular coordinate With ellipsoidal coordinates The partial derivative matrix of ; The second-order derivative of the gravitational potential with respect to time, that is, the partial derivative of the perturbation acceleration with respect to the position of the detector at any time, is expressed as or , the formula is: in, For function Ellipsoid coordinates of the calculated point The partial derivative of , The value is a rectangular coordinate Time-division term function The function value of i =1,2,3; ,in, , , and All are item functions; Itemized Function Calculated by the following formula: Among them, m=1,2,3.

5. The method for recovering the harmonic gravity field of the Phobos ellipsoid based on Doppler orbit data according to claim 4, characterized in that: In step 1, the ellipsoidal harmonic gravity field formula is: in, V is the gravity position, is the gravitational constant, is the mass of the central celestial body, n is the order, and the superscript p is the degree; is the ellipsoidal harmonic gravity potential coefficient when the order is n and the degree is p; Partial Function By integrating the form: Among them, the focal length ;focal length ; a , b , c are the three semi-axis lengths of the reference ellipsoid; here t represents the integration node, The value is taken according to the specified interval within the range. n is the order; The calculation of is similar, is the maximum semi-axis length of the reference ellipsoid; express t At the node, the order is n , the number of times is p， Ellipsoid coordinates Take The normalized Lamé function of the first kind when ; Partial Function Calculated using the following formula: in, are the coefficients of the Lamé function polynomial, i The values ​​are 1, 2, and 3; the value of m is different depending on the type of Lamé function; Itemized Function The values ​​of are as follows: When the classification is K, the order is 1 if it is an even number and ; When classified as L, the order is even. , if the order is odd, then ; When classified as M, the order is even. , if the order is odd, then ; When the classification is N categories, the order is even. , if the order is odd, then ; in, , , The Cartesian coordinates of the calculation point.

6. The method for recovering the harmonic gravity field of the Phobos ellipsoid based on Doppler orbit data according to claim 5, characterized in that: In step 1, the partial derivative of the perturbation acceleration with respect to the estimated parameter is calculated. dadp The method is as follows: The partial derivative of the perturbation acceleration with respect to the estimated parameter GM is: in, represents the partial derivative of the gravity potential V with respect to the rectangular coordinates of the calculation point, i represents the subscript of the rectangular coordinates of the calculation point, and its value is 1, 2, 3; Calculate the harmonic gravity potential coefficient of the ellipsoidal acceleration The partial derivative of is: in, is the perturbation acceleration, is the N-body perturbation force, The order is n , the number of times is p The second kind of Lamé function when represents the maximum semi-axis length of the reference ellipsoid, is the item function Calculate the rectangular coordinates of the point The partial derivative of There are four types of Lamé functions.

7. The method for recovering the harmonic gravity field of the Phobos ellipsoid based on Doppler orbit data according to claim 1, characterized in that: The method for obtaining the design matrix is ​​as follows: The two-way speed observation value Right now For any time t Satellite status The partial derivative of , from any time When mapped to the initial The partial derivative of the time, that is, the observed value at any time The partial derivative of the orbital vector , expressed as: in, for Satellite status at all times; is the design matrix, representing any time The partial derivative matrix of the detector state vector at time t to the state vector at the initial time and the parameters to be estimated; in, for Satellite position vector at time, for Satellite position vector at time, for Satellite velocity vector at time, is the parameter to be estimated; for Satellite velocity vector at this moment; for Satellite velocity vector at this moment; Design Matrix The first-order derivative is: in, n p is the number of parameters to be estimated; dadr represents the partial derivative of the perturbation acceleration with respect to the detector position at any time; dadv represents the partial derivative of the perturbation acceleration with respect to the velocity vector; dadp Represents the partial derivative of the perturbation acceleration with respect to the estimated parameters, including the perturbation acceleration with respect to the GM term and the ellipsoidal harmonic gravity potential coefficient The partial derivative of The second-order derivative of the design matrix is: for Satellite position vector at time For initial Satellite position vector at time The partial derivative of for Satellite position vector at time Parameters to be estimated The partial derivative of Design Matrix The position vector of the detector at any time can be obtained by integrating the first-order derivative and the second-order derivative.

8. The method for recovering the harmonic gravity field of the Phobos ellipsoid based on Doppler orbit data according to claim 1, characterized in that: In step 2, the normal equation is determined according to the linear unbiased minimum variance: in, is the weight matrix, representing the accuracy difference between observations, , represents the observation noise; Represents the observed value at any time The partial derivative of the orbital vector when , for The transposed matrix of is the correction value of the parameter to be estimated, is the prior value of the parameter to be estimated, is the prior error matrix of the parameters to be estimated, for The residual of the observation at time; Substitute the partial derivatives obtained in step 1 into the normal equation and replace the estimated parameters Local Parameters and global parameters , the single arc segment normal equation can be expressed as: in, is a local parameter The normal equation matrix, For global parameters The normal equation matrix, and are transposed to each other, representing the connection matrix between local parameters and global parameters, z 1 is the right side of the equal sign of the local parameter solution formula corresponding to the linear unbiased least squares method , represents the local parameter prior value, z 2 is the right side of the equal sign of the global parameter solution formula corresponding to the linear unbiased least squares method , Represents the global parameter prior value.

9. The method for recovering the harmonic gravity field of the Phobos ellipsoid based on Doppler orbit data according to claim 8, characterized in that: In step 3, the normal equation after fusion of multiple arc segments is: in, The superscripts in the middle represent different arc segments, and the superscript n represents the number of the merged arc segment.

10. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: When the processor executes the program, the method for recovering the harmonic gravity field of the Phobos ellipsoid based on Doppler orbit data as claimed in any one of claims 1 to 9 is implemented.

Images