Estimation method for entering disturbance compensation state of planet atmosphere

Abstract

The invention discloses a planetary atmosphere entering disturbance compensation state estimation method, and belongs to the technical field of deep space exploration. The implementation method of the invention comprises the following steps: establishing a planet atmosphere entry kinetic model, an uncertain parameter disturbance model and a navigation measurement model, wherein uncertain parameters comprise atmospheric density, ballistic coefficient and lift-drag ratio; the accelerometer and radio distance measurement and speed measurement combined navigation is adopted, so that the observability of the navigation system is improved. Aiming at the kinetic model considering the influence of the parameter uncertainty, designing a non-linear disturbance observer, and estimating the disturbance caused by the parameter uncertainty through the non-linear disturbance observer. Introducing an interference estimation value to correct the filtering dynamic model and the navigation measurement model, and introducing the interference estimation value into a filtering process to compensate the disturbance; and performing state estimation by adopting a rank filtering method. According to the method, parameter disturbance is observed, disturbance compensation is carried out in the filtering process, and the navigation state estimation precision under the uncertain condition is improved.

Description

Technical Field

[0001] The invention relates to a state estimation method, in particular to a planetary atmosphere entry disturbance compensation state estimation method, and belongs to the technical field of deep space exploration. Background Art

[0002] Planetary landing exploration and sample return are the most challenging planetary exploration activities. Carrying out planetary landing and sample return missions is of great significance for exploring the origin of life, studying the evolution of the earth, and improving aerospace technology. Safe and accurate landing on the surface of a planet is a prerequisite for successfully carrying out landing exploration missions and sample return. High-precision autonomous navigation and guidance control during the planetary entry, descent, and landing phases is the basis for achieving safe and accurate planetary landing.

[0003] High-precision autonomous navigation during planetary atmospheric entry is a prerequisite for active compensation and precise control, and is an important basis for achieving safe and precise landing on the surface of a planet. Due to the complex and changeable environment of planetary atmospheric entry and the limitations of existing technologies, there are generally uncertainties in parameters such as atmospheric density, ballistic coefficient, and lift-to-drag ratio during the planetary atmospheric entry process. These uncertainties will cause disturbances to the motion characteristics of the probe. Parameter uncertainty causes large model errors in the navigation model. When performing navigation filter estimation, it is necessary to consider the impact of parameter uncertainty and compensate or suppress the disturbance caused by parameter uncertainty. Otherwise, the accuracy of state estimation may decrease, or even filter divergence may occur. This will not only seriously affect the effect of guidance and control and reduce landing accuracy, but may also cause the probe to deviate from the predetermined trajectory, causing yaw, and ultimately leading to the failure of the entire landing mission. Summary of the invention

[0004] In order to solve the problem that parameter uncertainty during the planetary atmosphere entry process causes disturbances to the detector, thereby affecting the state estimation accuracy of the navigation system, the purpose of the present invention is to provide a planetary atmosphere entry disturbance compensation state estimation method, which first establishes a planetary atmosphere entry dynamics model, an uncertain parameter disturbance model and a navigation measurement model; uses accelerometer and radio ranging and speed measurement combined navigation to improve the observability of the navigation system. For the dynamic model considering the influence of parameter uncertainty, a nonlinear disturbance observer is designed, and the disturbance caused by parameter uncertainty is estimated by the nonlinear disturbance observer. The interference estimate value is introduced to correct the filtering dynamics model and the navigation measurement model to achieve compensation for the disturbance; and the rank filtering method is used for state estimation. The present invention improves the navigation state estimation accuracy under uncertainty conditions by observing the parameter disturbance and performing disturbance compensation during the filtering process.

[0005] The purpose of the present invention is achieved through the following technical solutions.

[0006] The planetary atmosphere entry disturbance compensation state estimation method disclosed in the present invention comprises the following steps:

[0007] Step 1: Establish a three-degree-of-freedom dynamic model for planetary atmosphere entry; consider the impact of parameter uncertainty on the motion characteristics of the probe and establish an uncertain parameter perturbation model; the uncertain parameters include atmospheric density, ballistic coefficient and lift-to-drag ratio; use an accelerometer and radio combined navigation scheme to establish a corresponding navigation measurement model.

[0008] X = [r T ,v T ] T is the state quantity, where the detector position vector r = [r x ,r y ,r z ] T , the detector velocity vector v = [v x ,v y ,v z ] T , the planetary atmosphere enters the three-degree-of-freedom dynamic model:

[0009]

[0010] Where σ is the roll angle, g(r) = -μ / ||r|| 2 is the local gravitational acceleration, μ is the gravitational constant of the planet; D and L are the drag and lift accelerations respectively, and the specific expression is

[0011]

[0012] Where v = ||v|| is the velocity, ρ is the real atmospheric density, S is the reference area, m s is the detector mass, C D and C L are the drag coefficient and lift coefficient, respectively, q d =ρv 2 / 2 is the dynamic pressure, β=m s / SC D is the true ballistic coefficient, η=L / D is the true lift-to-drag ratio;

[0013] Considering the influence of atmospheric density, ballistic coefficient, and lift-to-drag ratio uncertainty on atmospheric entry dynamics; at any time, the true atmospheric density value of the planet is expressed by the following formula

[0014] ρ=ρ * (1+Δ ρ ) (3)

[0015] Among them, ρ represents the real atmospheric density at the current moment, ρ *represents the nominal value of the atmospheric density model calculated by equation (4), Δ ρ It indicates the percentage deviation between the model nominal value and the true value;

[0016]

[0017] Among them, ρ0 is the reference atmospheric density, h is the current altitude of the detector, h0 is the reference altitude, and h s is the proportional height;

[0018] Similarly, the true ballistic coefficient β and the true lift-to-drag ratio η are expressed as follows:

[0019] β=β * (1+Δ β ) (5)

[0020] η=η * (1+Δ η ) (6)

[0021] Among them, β * and η * are the model nominal values ​​of ballistic coefficient and lift-to-drag ratio, Δ β and Δ η They represent the percentage deviations between the true and model nominal values ​​of the ballistic coefficient and lift-to-drag ratio, respectively;

[0022] The navigation equipment carried by the probe includes an accelerometer and a radio transceiver; the accelerometer measurement model is simplified as follows:

[0023]

[0024] in, is the acceleration information measured by the accelerometer, d is the interference term caused by parameter uncertainty, ε a is the accelerometer measurement noise, defined as zero-mean white noise; a B is the non-gravitational acceleration, calculated by the following formula:

[0025]

[0026] Assume that the position vector and velocity vector of the i-th radio beacon in the planetary inertial coordinate system are r Bi =[r xBi ,r yBi ,r zBi ] T , v Bi =[v xBi ,v yBi ,v zBi ] T , then the relative distance and relative speed between the detector and the i-th radio beacon are

[0027]

[0028] Therefore, the radio measurement model is expressed as:

[0029]

[0030] Where N is the number of radio beacons, ε R and ε V They are radio ranging noise and speed measurement noise, respectively, both defined as zero-mean white noise;

[0031] In summary, the navigation measurement model is expressed as:

[0032]

[0033] Where H(·) is a nonlinear vector function and ε is the measurement noise.

[0034] Step 2: Consider the influence of parameter uncertainty and rewrite the dynamic model established in step 1; on this basis, design a nonlinear disturbance observer; use the nonlinear disturbance observer to estimate the disturbance caused by parameter uncertainty and obtain the disturbance estimate

[0035] Considering the influence of parameter uncertainty, the dynamic equation (1) is rewritten as:

[0036]

[0037] Where X = [r T ,v T ] T , d is the disturbance term caused by parameter uncertainty, and the system control input u=σ;

[0038] The nonlinear disturbance observer is expressed as:

[0039]

[0040] in, is the disturbance estimate output by the observer, z is the internal variable of the nonlinear disturbance observer, Q=diag(λ1,λ2,λ3) is the observer output gain, λ1,λ2,λ3 are all constants greater than 0; the input state of the nonlinear disturbance observer, namely X and v, is obtained by dead reckoning based on the actual measurement value of the accelerometer;

[0041] The observation error of the nonlinear disturbance observer is defined as

[0042]

[0043] The candidate Lyapunov function is selected as follows to perform stability analysis on the nonlinear disturbance observer

[0044]

[0045] but

[0046]

[0047] According to equations (12), (13) and (14), we can get

[0048]

[0049] Pick Therefore, the observation error equation is

[0050]

[0051] Substituting equation (18) into equation (16), we can obtain

[0052]

[0053] Let λ=min(2λ1,2λ2,2λ3), we can get inequality (20)

[0054] Q T +Q≥λI (20)

[0055] Where I is the identity matrix, then

[0056]

[0057] Combining formula (15), we can get

[0058]

[0059] Lyapunov function V L is a function of time t, and the value of the Lyapunov function at time t is V L (t), then the solution of equation (22) is

[0060] V L (t)≤V L (0)exp(-λt) (23)

[0061] When t→∞, V L (t)→0, so Exponential convergence;

[0062] In summary, according to the Lyapunov theory, the designed disturbance observer is stable.

[0063] Step 3: Use the disturbance estimate output by the nonlinear disturbance observer in step 2 The rewritten dynamic model in step 2 and the navigation measurement model in step 1 are corrected, and the rank filtering method is used for state estimation, that is, the planetary atmosphere entry disturbance compensation state estimation is realized.

[0064] Step 3.1: Based on the state estimate at time k-1 And the error variance matrix P k-1 , get the state prediction at time k

[0065] Discretize equation (12) as follows:

[0066] X k =F(X k-1 ,d k-1 )+W (24)

[0067] Among them, X k is the state quantity at time k, X k-1 is the state quantity at time k-1; d k-1 is the interference term caused by parameter uncertainty at time k-1; F(·) is a nonlinear vector function; W is the system noise;

[0068] The disturbance estimate obtained by introducing the nonlinear disturbance observer at time k-1 is By modifying equation (12), the state prediction at time k is obtained as follows:

[0069]

[0070] Among them, m is the number of sampling layers, n is the dimension of the state quantity, and there are 2mn sample points in total; the intermediate variable X k / (k-1),i It is expressed as:

[0071]

[0072] Among them, χ k-1,i for The i-th sampling point of:

[0073]

[0074] in, represents the state estimate at time k-1, Represents the matrix P k-1 The i-th column vector of square roots; standard normal skewness Sampling correction coefficient r'=1;

[0075] Step 3.2: State prediction based on step 3.1 Get the measurement estimate at time k

[0076] Re-rank sampling, and get the i-th new sampling point as

[0077]

[0078] Among them, χ k / (k-1),i for The i-th sampling point; Represents the matrix P k / (k-1) The i-th column vector of the square root; the prediction error variance matrix P k / (k-1) for:

[0079]

[0080] in, is the covariance correction coefficient corresponding to the i-th sampling point, which can generally be taken as i=1,2,...,2mn;covariance weight coefficient Q W is the system noise variance matrix;

[0081] The disturbance estimate obtained by introducing the nonlinear disturbance observer at time k-1 is Modify equation (7); the measurement estimate at time k is:

[0082]

[0083] Among them, the intermediate variable Z k / (k-1),i It is expressed as:

[0084]

[0085] Step 3.3: State prediction based on steps 3.1 and 3.2 and measurement estimation Combined with the real measurement value Z of the navigation system k , get the state estimate at time k And the error variance matrix P k ;

[0086] The filter gain matrix at time k is:

[0087]

[0088] in

[0089]

[0090]

[0091] Among them, R εis the measurement noise variance matrix;

[0092] In summary, the state estimation at time k is Its error variance matrix P k They are

[0093]

[0094]

[0095] Among them, Z k is the true measurement value of the navigation system at time k, that is, the true measurement value of the accelerometer and radio.

[0096] Beneficial effects:

[0097] 1. The planetary atmosphere entry disturbance compensation state estimation method disclosed in the present invention introduces a nonlinear disturbance observer to observe the disturbance term caused by parameter uncertainty, and introduces the disturbance estimation value into the filtering process to compensate for the disturbance, thereby improving the state estimation accuracy under uncertainty conditions.

[0098] 2. The planetary atmosphere entry disturbance compensation state estimation method disclosed in the present invention adopts a rank filtering method to fuse multi-source information, and effectively simulates the probability distribution of the system state through a rank sampling method, thereby further improving the accuracy of navigation state estimation. BRIEF DESCRIPTION OF THE DRAWINGS

[0099] Figure 1 Flow chart of the state estimation method for planetary atmosphere entry disturbance compensation;

[0100] Figure 2 is the relationship between the estimated value of the three-axis interference and time, where: Figure 2 (a) is the relationship between the estimated value of X-axis interference and time. Figure 2 (b) is the relationship between the Y-axis interference estimation value and time. Figure 2 (c) is the relationship between the estimated value of Z-axis interference and time;

[0101] Figure 3 is the relationship between the three-axis interference estimation error and time, where Figure 3 (a) is the relationship between the X-axis interference estimation error and time. Figure 3 (b) is the relationship between the Y-axis interference estimation error and time. Figure 3 (c) is the relationship between the Z-axis interference estimation error and time;

[0102] Figure 4 is the relationship between the three-axis position estimation error and time, where Figure 4 (a) is the relationship between the X-axis position estimation error and time. Figure 4(b) is the relationship between the Y-axis position estimation error and time. Figure 4 (c) is the relationship between the Z-axis position estimation error and time;

[0103] Figure 5 is the relationship between the three-axis velocity estimation error and time, where Figure 5 (a) is the relationship between the X-axis velocity estimation error and time. Figure 5 (b) is the relationship between the Y-axis velocity estimation error and time. Figure 5 (c) is a graph showing the relationship between the Z-axis velocity estimation error and time. DETAILED DESCRIPTION

[0104] In order to better illustrate the purpose and advantages of the present invention, the invention is further described below in conjunction with an embodiment and corresponding drawings.

[0105] Taking the Mars Science Laboratory rover mission as the background, the Mars atmosphere entry navigation process with parameter uncertainty is simulated, the state of the rover is estimated by the method proposed in this invention, and the estimated state value is compared with the actual value to verify the effectiveness of the proposed method. The simulation parameters are shown in Table 1. The roll angle σ=0 during the entire atmospheric entry process. The atmospheric density uncertainty Δ ρ = 20%, uncertainty of ballistic coefficient Δ β =10%, lift-to-drag ratio uncertainty Δ η =10%. The nonlinear disturbance observer gain is set to Q=diag(150,100,100). The detector initial entry state is shown in Table 2.

[0106] Table 1 Simulation parameters

[0107]

[0108] Table 2 Detector enters initial state

[0109]

[0110] To simplify the simulation, the response frequency of the accelerometer and the radio is set to 100 Hz, the total simulation time is set to 250 s, and the simulation integration step is 0.01 s. The accelerometer measurement noise variance is [5×10 -8 ,5×10 -8 ,5×10 -8 ](m / s 2 ) 2 , the radio ranging noise variance is 400m 2 , the speed measurement noise variance is 0.4 (m / s) 2 The system noise variance matrix is ​​set to Q W=diag(10,10,10,0.1,0.1,0.1). The present invention uses three orbiters as radio beacons, and the initial states of the orbiters are shown in Table 3.

[0111] Table 3 Orbiter initial state

[0112]

[0113] like Figure 1 As shown, the planetary atmosphere entry disturbance compensation state estimation method disclosed in this embodiment has the following specific implementation steps:

[0114] Step 1: Establish a three-degree-of-freedom dynamic model for planetary atmosphere entry; consider the impact of parameter uncertainty on the motion characteristics of the probe and establish an uncertain parameter perturbation model; the uncertain parameters include atmospheric density, ballistic coefficient and lift-to-drag ratio; use an accelerometer and radio combined navigation scheme to establish a corresponding navigation measurement model.

[0115] X = [r T ,v T ] T is the state quantity, where the detector position vector r = [r x ,r y ,r z ] T , the detector velocity vector v = [v x ,v y ,v z ] T , the planetary atmosphere enters the three-degree-of-freedom dynamic model:

[0116]

[0117] Where σ is the roll angle, g(r) = -μ / ||r|| 2 is the local gravitational acceleration, μ is the gravitational constant of the planet; D and L are the drag and lift accelerations respectively, and the specific expression is

[0118]

[0119] Where v = ||v|| is the velocity, ρ is the real atmospheric density, S is the reference area, m s is the detector mass, C D and C L are the drag coefficient and lift coefficient, respectively, q d =ρv 2 / 2 is the dynamic pressure, β=m s / SC D is the true ballistic coefficient, η=L / D is the true lift-to-drag ratio;

[0120] Considering the influence of atmospheric density, ballistic coefficient, and lift-to-drag ratio uncertainty on atmospheric entry dynamics; at any time, the true atmospheric density value of the planet is expressed by the following formula

[0121] ρ=ρ * (1+Δ ρ ) (39)

[0122] Among them, ρ represents the real atmospheric density at the current moment, ρ * represents the nominal value of the atmospheric density model calculated by equation (4), Δ ρ It indicates the percentage deviation between the model nominal value and the true value;

[0123]

[0124] Among them, ρ0 is the reference atmospheric density, h is the current altitude of the detector, h0 is the reference altitude, and h s is the proportional height;

[0125] Similarly, the true ballistic coefficient β and the true lift-to-drag ratio η are expressed as follows:

[0126] β=β * (1+Δ β ) (41)

[0127] η=η * (1+Δ η ) (42)

[0128] Among them, β * and η * are the model nominal values ​​of ballistic coefficient and lift-to-drag ratio, Δ β and Δ η They represent the percentage deviations between the true and model nominal values ​​of the ballistic coefficient and lift-to-drag ratio, respectively;

[0129] The navigation equipment carried by the probe includes an accelerometer and a radio transceiver; the accelerometer measurement model is simplified as follows:

[0130]

[0131] in, is the acceleration information measured by the accelerometer, d is the interference term caused by parameter uncertainty, ε a is the accelerometer measurement noise, defined as zero-mean white noise; a B is the non-gravitational acceleration, calculated by the following formula:

[0132]

[0133] Assume that the position vector and velocity vector of the i-th radio beacon in the planetary inertial coordinate system are rBi =[r xBi ,r yBi ,r zBi ] T , v Bi =[v xBi ,v yBi ,v zBi ] T , then the relative distance and relative speed between the detector and the i-th radio beacon are

[0134]

[0135] Therefore, the radio measurement model is expressed as:

[0136]

[0137] Where N is the number of radio beacons, ε R and ε V They are radio ranging noise and speed measurement noise, respectively, both defined as zero-mean white noise;

[0138] In summary, the navigation measurement model is expressed as:

[0139]

[0140] Where H(·) is a nonlinear vector function and ε is the measurement noise.

[0141] Step 2: Consider the influence of parameter uncertainty and rewrite the dynamic model established in step 1; on this basis, design a nonlinear disturbance observer; use the nonlinear disturbance observer to estimate the disturbance caused by parameter uncertainty.

[0142] Considering the influence of parameter uncertainty, the dynamic equation (1) is rewritten as:

[0143]

[0144] Where X = [r T ,v T ] T , d is the disturbance term caused by parameter uncertainty, and the system control input u=σ;

[0145] The nonlinear disturbance observer is expressed as:

[0146]

[0147] in, is the disturbance estimate output by the observer, z is the internal variable of the nonlinear disturbance observer, Q=diag(λ1,λ2,λ3) is the observer output gain, λ1,λ2,λ3 are all constants greater than 0; the input state of the nonlinear disturbance observer, namely X and v, is obtained by dead reckoning based on the actual measurement value of the accelerometer.

[0148] The observation error of the nonlinear disturbance observer is defined as

[0149]

[0150] The candidate Lyapunov function is selected as follows to perform stability analysis on the nonlinear disturbance observer

[0151]

[0152] but

[0153]

[0154] According to equations (12), (13) and (14), we can get

[0155]

[0156] Pick Therefore, the observation error equation is

[0157]

[0158] Substituting equation (18) into equation (16), we can obtain

[0159]

[0160] Let λ=min(2λ1,2λ2,2λ3), we can get inequality (20)

[0161] Q T +Q≥λI (56)

[0162] Where I is the identity matrix, then

[0163]

[0164] Combining formula (15), we can get

[0165]

[0166] Lyapunov function V L is a function of time t, and the value of the Lyapunov function at time t is V L (t), then the solution of equation (22) is

[0167] VL (t)≤V L (0)exp(-λt) (59)

[0168] When t→∞, V L (t)→0, so Exponential convergence;

[0169] In summary, according to the Lyapunov theory, the designed disturbance observer is stable.

[0170] Step 3: Use the disturbance estimate output by the nonlinear disturbance observer in step 2 The equation (12) in step 2 and the navigation measurement model in step 1 are modified, and the rank filtering method is used for state estimation, that is, the planetary atmosphere entry disturbance compensation state estimation is realized.

[0171] Step 3.1: Based on the state estimate at time k-1 And the error variance matrix P k-1 , get the state prediction at time k

[0172] Discretize equation (12) as follows:

[0173] X k =F(X k-1 ,d k-1 )+W (60)

[0174] Among them, X k is the state quantity at time k, X k-1 is the state quantity at time k-1; d k-1 is the interference term caused by parameter uncertainty at time k-1; F(·) is a nonlinear vector function; W is the system noise;

[0175] The disturbance estimate obtained by introducing the nonlinear disturbance observer at time k-1 is By modifying equation (12), the state prediction at time k is obtained as follows:

[0176]

[0177] Among them, m is the number of sampling layers, n is the dimension of the state quantity, and there are 2mn sample points in total; the intermediate variable X k / (k-1),i It is expressed as:

[0178]

[0179] Among them, χ k-1,i for The i-th sampling point of:

[0180]

[0181] in, represents the state estimate at time k-1, Represents the matrix P k-1 The i-th column vector of square roots; standard normal skewness Sampling correction coefficient r'=1;

[0182] Step 3.2: State prediction based on step 3.1 Get the measurement estimate at time k

[0183] Re-rank sampling, and get the i-th new sampling point as

[0184]

[0185] Among them, χ k / (k-1),i for The i-th sampling point; Represents the matrix P k / (k-1) The i-th column vector of the square root; the prediction error variance matrix P k / (k-1) for:

[0186]

[0187] in, is the covariance correction coefficient corresponding to the i-th sampling point, which can generally be taken as i=1,2,...,2mn;covariance weight coefficient Q W is the system noise variance matrix;

[0188] The disturbance estimate obtained by introducing the nonlinear disturbance observer at time k-1 is Modify equation (7); the measurement estimate at time k is:

[0189]

[0190] Among them, the intermediate variable Z k / (k-1),i It is expressed as:

[0191]

[0192] Step 3.3: State prediction based on steps 3.1 and 3.2 and measurement estimation Combined with the real measurement value Z of the navigation system k , get the state estimate at time k And the error variance matrix P k ;

[0193] The filter gain matrix at time k is:

[0194]

[0195] in

[0196]

[0197]

[0198] Among them, R ε is the measurement noise variance matrix;

[0199] In summary, the state estimation at time k is Its error variance matrix P k They are

[0200]

[0201]

[0202] Among them, Z k is the true measurement value of the navigation system at time k, that is, the true measurement value of the accelerometer and radio.

[0203] The relationship between the estimated value of the three-axis disturbance obtained by the nonlinear disturbance observer and time is as follows: Figure 2 The simulation results show that in the early stage of the entry process (the first 50 seconds), the disturbance caused by parameter uncertainty is relatively small due to the relatively thin atmosphere. The disturbance increases significantly from 50 seconds to 200 seconds of the entry process. The relationship between the three-axis disturbance estimation error of the nonlinear disturbance observer and time is shown in Figure 3 As shown, the interference estimation error is on the order of 10 -3 m / s 2 ,During the period from 50 seconds to 200 seconds of the entry process, the relative error of interference estimation is maintained below 1%, with high estimation accuracy.

[0204] When there is parameter uncertainty, directly applying rank filtering will lead to poor state estimation accuracy, or even filter divergence, and the detector state cannot be correctly estimated. The interference estimate is used to correct the dynamic model and navigation observation model. The detector state estimation result is as follows: Figure 4 and Figure 5 As shown in the figure, the X-axis position estimation error gradually converges to within 10m, the Y-axis and Z-axis position estimation errors gradually converge to within 5m; the X-axis velocity estimation error converges to within 2m / s, and the Y-axis and Z-axis velocity estimation errors all converge to within 0.5m / s. It can be seen from the simulation results that the method proposed in the present invention can achieve high-precision state estimation under parameter uncertainty conditions.

[0205] The specific description above further illustrates the purpose, technical solutions and beneficial effects of the invention in detail. It should be understood that the above is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A planetary atmosphere entry disturbance compensation state estimation method, characterized in that: The following steps are included: Step 1: Establish a three-degree-of-freedom dynamic model for planetary atmosphere entry; consider the influence of parameter uncertainty on the motion characteristics of the probe and establish an uncertain parameter perturbation model; the uncertain parameters include atmospheric density, ballistic coefficient and lift-to-drag ratio; use accelerometer and radio combined navigation to establish a corresponding navigation measurement model; Step 2: Consider the influence of parameter uncertainty and rewrite the dynamic model established in step 1; on this basis, design a nonlinear disturbance observer; use the nonlinear disturbance observer to estimate the disturbance caused by parameter uncertainty and obtain the disturbance estimate Step 3: Use the disturbance estimate output by the nonlinear disturbance observer in step 2 The rewritten dynamic model in step 2 and the navigation measurement model in step 1 are corrected, and the rank filtering method is used for state estimation, that is, the planetary atmosphere entry disturbance compensation state estimation is realized.

2. The planetary atmosphere entry disturbance compensation state estimation method according to claim 1, characterized in that: The specific implementation method of step one is: X = [r T ,v T ] T is the state quantity, where the detector position vector r = [r x ,r y ,r z ] T , the detector velocity vector v = [v x ,v y ,v z ] T , the planetary atmosphere enters the three-degree-of-freedom dynamic model: Where σ is the roll angle, g(r) = -μ / ||r|| 2 is the local gravitational acceleration, μ is the gravitational constant of the planet; D and L are the drag and lift accelerations respectively, and the specific expression is Where v = ||v|| is the velocity, ρ is the real atmospheric density, S is the reference area, m s is the detector mass, C D and C L are the drag coefficient and lift coefficient, respectively, q d =ρv 2 / 2 is the dynamic pressure, β=m s / SC D is the true ballistic coefficient, η=L / D is the true lift-to-drag ratio; Considering the influence of atmospheric density, ballistic coefficient, and lift-to-drag ratio uncertainty on atmospheric entry dynamics; at any time, the true atmospheric density value of the planet is expressed by the following formula p=p * (1+D ρ ) (3) Among them, ρ represents the real atmospheric density at the current moment, ρ * represents the nominal value of the atmospheric density model calculated by equation (4), Δ ρ It indicates the percentage deviation between the model nominal value and the true value; Among them, ρ0 is the reference atmospheric density, h is the current altitude of the detector, h0 is the reference altitude, and h s is the proportional height; The true ballistic coefficient β and the true lift-to-drag ratio η are expressed as follows: β=β * (1+D β ) (5) the=the * (1+D η ) (6) Among them, β * and η * are the model nominal values ​​of ballistic coefficient and lift-to-drag ratio, Δ β and Δ η They represent the percentage deviations between the true and model nominal values ​​of the ballistic coefficient and lift-to-drag ratio, respectively; The navigation equipment carried by the probe includes an accelerometer and a radio transceiver; the accelerometer measurement model is simplified as follows: in, is the acceleration information measured by the accelerometer, d is the interference term caused by parameter uncertainty, ε a is the accelerometer measurement noise, defined as zero-mean white noise; a B is the non-gravitational acceleration, calculated by the following formula: Assume that the position vector and velocity vector of the i-th radio beacon in the planetary inertial coordinate system are r Bi =[r xBi ,r yBi ,r zBi ] T , v Bi =[v xBi ,v yBi ,v zBi ] T , then the relative distance and relative speed between the detector and the i-th radio beacon are The radio measurement model is expressed as: Where N is the number of radio beacons, ε R and ε V They are radio ranging noise and speed measurement noise, respectively, both defined as zero-mean white noise; The navigation measurement model is expressed as: Where H(·) is a nonlinear vector function and ε is the measurement noise.

3. The planetary atmosphere entry disturbance compensation state estimation method according to claim 2, characterized in that: The specific implementation method of step 2 is: Considering the influence of parameter uncertainty, the dynamic equation (1) is rewritten as: Where X = [r T ,v T ] T , d is the disturbance term caused by parameter uncertainty, and the system control input u=σ; The nonlinear disturbance observer is expressed as: in, is the disturbance estimate output by the observer, z is the internal variable of the nonlinear disturbance observer, Q=diag(λ1,λ2,λ3) is the observer output gain, λ1,λ2,λ3 are all constants greater than 0; the input state of the nonlinear disturbance observer, namely X and v, is obtained by dead reckoning based on the actual measurement value of the accelerometer.

4. The planetary atmosphere entry disturbance compensation state estimation method according to claim 3, characterized in that: The specific implementation method of step three is: Step 3.1: Based on the state estimate at time k-1 And the error variance matrix P k-1 , get the state prediction at time k Discretize equation (12) as follows: X k =F(X k-1 ,d k-1 )+W (14) Among them, X k is the state quantity at time k, X k-1 is the state quantity at time k-1; d k-1 is the interference term caused by parameter uncertainty at time k-1; F(·) is a nonlinear vector function; W is the system noise; The disturbance estimate obtained by introducing the nonlinear disturbance observer at time k-1 is By modifying equation (12), the state prediction at time k is obtained as follows: Among them, m is the number of sampling layers, n is the dimension of the state quantity, and there are 2mn sample points in total; the intermediate variable X k / (k-1),i It is expressed as: Among them, χ k-1,i for The i-th sampling point of: in, represents the state estimate at time k-1, Represents the matrix P k-1 The i-th column vector of square roots; standard normal skewness Sampling correction coefficient r'=1; Step 3.2: State prediction based on step 3.1 Get the measurement estimate at time k Re-rank sampling, and get the i-th new sampling point as Among them, χ k / (k-1),i for The i-th sampling point; Represents the matrix P k / (k-1) The i-th column vector of the square root; the prediction error variance matrix P k / (k-1) for: Among them, r i * is the covariance correction coefficient corresponding to the i-th sampling point, taking r i * =1, i = 1, 2, ..., 2mn; covariance weight coefficient Q W is the system noise variance matrix; The disturbance estimate obtained by introducing the nonlinear disturbance observer at time k-1 is Modify equation (7); the measurement estimate at time k is: Among them, the intermediate variable Z k / (k-1),i It is expressed as: Step 3.3: State prediction based on steps 3.1 and 3.2 and measurement estimation Combined with the real measurement value Z of the navigation system k , get the state estimate at time k And the error variance matrix P k ; The filter gain matrix at time k is: in Among them, R ε is the measurement noise variance matrix; In summary, the state estimation at time k is Its error variance matrix P k They are Among them, Z k is the true measurement value of the navigation system at time k, that is, the true measurement value of the accelerometer and radio.

5. The planetary atmosphere entry disturbance compensation state estimation method according to claim 4, characterized in that: The stability of the nonlinear disturbance observer in step 2 is verified by the following method: The observation error of the nonlinear disturbance observer is defined as The candidate Lyapunov function is selected as follows to perform stability analysis on the nonlinear disturbance observer but According to equations (12), (13) and (27), we can get Pick Therefore, the observation error equation is Substituting equation (31) into equation (29), we get Let λ=min(2λ1,2λ2,2λ3), and we get inequality (33) Q T +Q≥λI (33) Where I is the identity matrix, then Combining formula (28), we can get Lyapunov function V L is a function of time t, and the value of the Lyapunov function at time t is V L (t), then the solution of equation (35) is V L (t)≤V L (0)exp(-λt) (36) When t→∞, V L (t)→0, so Exponential convergence; According to Lyapunov theory, it is verified that the designed disturbance observer is stable.

Images