Planetary atmospheric entry perturbation compensation state estimation method

Abstract

The planet atmosphere entry disturbance compensation state estimation method disclosed by the application belongs to the technical field of deep space exploration. The method of the application is as follows: a planet atmosphere entry dynamics model, an uncertain parameter disturbance model and a navigation measurement model are established, the uncertain parameters including atmosphere density, ballistic coefficient and lift-drag ratio; an accelerometer and a radio ranging and speed measurement combined navigation are adopted to improve the observability of the navigation system. For the dynamics model considering the influence of parameter uncertainty, a nonlinear disturbance observer is designed to estimate the disturbance caused by parameter uncertainty through the nonlinear disturbance observer. The disturbance estimation value is introduced to correct the filtering dynamics model and the navigation measurement model, and the disturbance estimation value is introduced into the filtering process to compensate the disturbance; a rank filtering method is adopted for state estimation. The application improves the navigation state estimation precision under the condition of uncertainty by observing the parameter disturbance and compensating the disturbance in the filtering process.

Description

Technical Field

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

[0002] Planetary landing and sample return missions are among the most challenging of planetary exploration activities. These missions are crucial for exploring the origins of life, studying Earth's evolution, and advancing space technology. Safe and precise landing on a planetary surface is a prerequisite for successful landing and sample return missions. High-precision autonomous navigation and guidance control during the planetary entry, descent, and landing phases are fundamental to achieving safe and precise planetary landings.

[0003] High-precision autonomous navigation during planetary atmospheric entry is a prerequisite for active compensation and precise control, and is crucial for achieving a safe and precise landing on a planetary surface. Due to the complex and variable environment of planetary atmospheric entry, coupled with the limitations of existing technologies, the process is subject to widespread uncertainty in parameters such as atmospheric density, ballistic coefficient, and lift-to-drag ratio. These uncertainties can perturb the probe's kinematic characteristics. Parameter uncertainty leads to significant errors in the navigation model. The impact of parameter uncertainty must be considered during navigation filter estimation, compensating for or suppressing the perturbations caused by parameter uncertainty. Failure to do so can lead to a decrease in state estimation accuracy and even filter divergence. This can severely impact guidance and control effectiveness, reducing landing accuracy, and potentially causing the probe to deviate from its intended trajectory, resulting in yaw and ultimately failure of the landing mission. Summary of the Invention

[0004] To address the problem of parameter uncertainty causing disturbances to the detector during planetary atmospheric entry, which in turn affects the state estimation accuracy of the navigation system, the present invention aims to provide a planetary atmospheric entry disturbance compensation state estimation method. This method first establishes a planetary atmospheric entry dynamics model, an uncertain parameter disturbance model, and a navigation measurement model; then uses accelerometers and radio ranging and velocity combined navigation to improve the observability of the navigation system. A nonlinear disturbance observer is designed for the dynamics model that considers the influence of parameter uncertainty, and the disturbance caused by parameter uncertainty is estimated using the nonlinear disturbance observer. Disturbance estimates are introduced to correct the filtering dynamics model and navigation measurement model to achieve disturbance compensation; and a rank filtering method is used for state estimation. The present invention improves the navigation state estimation accuracy under uncertainty conditions by observing parameter disturbances 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 probe's motion characteristics 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 three-degree-of-freedom dynamic model of the planetary atmosphere entering is:

[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 true 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, 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 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), Δ ρ 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 height of the detector, h0 is the reference height, and h s is the proportional height;

[0018] Similarly, the true ballistic coefficient β and 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 Δ η represent the percentage deviations between the true and model nominal values ​​of ballistic coefficient and lift-to-drag ratio, respectively;

[0022] The navigation equipment carried by the probe includes an accelerometer and a radio transceiver. The simplified accelerometer measurement model is 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 as follows:

[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 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 value.

[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, and λ1,λ2,λ3 are all constants greater than 0. The input states of the nonlinear disturbance observer, namely X and v, are obtained by dead reckoning based on the actual measurement values ​​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 formula (12), (13) and (14), we can get

[0048]

[0049] Pick Therefore, the observation error equation is

[0050]

[0051] Substituting formula (18) into formula (16) we can get

[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 As a function of time t, the value of the Lyapunov function at time t is V L (t), then the solution of formula (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 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 of step 2 and the navigation measurement model of 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: Estimation of the state at time k-1 and the error variance matrix P k-1 , get the state prediction at time k

[0065] Discretize Equation (12) as:

[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 Expressed as:

[0071]

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

[0073]

[0074] in, represents the state estimate at time k-1, Represents the matrix P k-1 The i-th column vector of the square root; the 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, the i-th new sampling point is

[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 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 and 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 estimate 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 Flowchart of the state estimation method for planetary atmospheric 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 X-axis interference estimation value 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 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] Using the Mars Science Laboratory rover mission as a backdrop, a simulation of the Mars atmospheric entry navigation process was conducted under conditions of parameter uncertainty. The proposed method was used to estimate the rover's state and compare the estimated state with the actual value to verify the effectiveness of the proposed method. The simulation parameters are shown in Table 1. Throughout the atmospheric entry process, the roll angle σ = 0. 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 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 is specifically implemented in the following steps:

[0114] Step 1: Establish a three-degree-of-freedom dynamic model for planetary atmosphere entry; consider the impact of parameter uncertainty on the probe's motion characteristics 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 three-degree-of-freedom dynamic model of the planetary atmosphere entering is:

[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 true 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, 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 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), Δ ρ 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 height of the detector, h0 is the reference height, and h s is the proportional height;

[0125] Similarly, the true ballistic coefficient β and 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 Δ η represent the percentage deviations between the true and model nominal values ​​of ballistic coefficient and lift-to-drag ratio, respectively;

[0129] The navigation equipment carried by the probe includes an accelerometer and a radio transceiver. The simplified accelerometer measurement model is 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 as follows:

[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 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, and λ1,λ2,λ3 are all constants greater than 0. The input states of the nonlinear disturbance observer, namely X and v, are obtained by dead reckoning based on the actual measurement values ​​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 formula (12), (13) and (14), we can get

[0155]

[0156] Pick Therefore, the observation error equation is

[0157]

[0158] Substituting formula (18) into formula (16) we can get

[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 As a function of time t, the value of the Lyapunov function at time t is V L (t), then the solution of formula (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 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: Estimation of the state at time k-1 and the error variance matrix P k-1 , get the state prediction at time k

[0172] Discretize Equation (12) as:

[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 Expressed as:

[0178]

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

[0180]

[0181] in, represents the state estimate at time k-1, Represents the matrix P k-1 The i-th column vector of the square root; the 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, the i-th new sampling point is

[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 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, we can get the state estimation at time k and 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 three-axis disturbance estimation value 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 into the entry process. The relationship between the three-axis disturbance estimation error of the nonlinear disturbance observer and time is shown as follows: 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%, which has a 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, and 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. The simulation results show that the proposed method can achieve high-precision state estimation under conditions of parameter uncertainty.

[0205] The above specific description further illustrates the purpose, technical solutions and beneficial effects of the invention in detail. It should be understood that the above description 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 by: The following steps are included: Step 1: Establish a three-degree-of-freedom dynamic model for planetary atmosphere entry; consider the impact of parameter uncertainty on the probe's motion characteristics and establish an uncertain parameter perturbation model; the uncertain parameters include atmospheric density, ballistic coefficient, and lift-to-drag ratio; use accelerometer and radio integrated 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 value. Step 3: Use the disturbance estimate output by the nonlinear disturbance observer in step 2 The rewritten dynamic model of step 2 and the navigation measurement model of 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 three-degree-of-freedom dynamic model of the planetary atmosphere entering is: 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 true 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, 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 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), Δ ρ 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 height of the detector, h0 is the reference height, and h s is the proportional height; The true ballistic coefficient β and 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 Δ η represent the percentage deviations between the true and model nominal values ​​of ballistic coefficient and lift-to-drag ratio, respectively; The navigation equipment carried by the probe includes an accelerometer and a radio transceiver. The simplified accelerometer measurement model is 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 as follows: 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 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, and λ1,λ2,λ3 are all constants greater than 0. The input states of the nonlinear disturbance observer, namely X and v, are obtained by dead reckoning based on the actual measurement values ​​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: Estimation of the state at time k-1 and the error variance matrix P k-1 , get the state prediction at time k Discretize Equation (12) as: 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 Expressed as: Among them, χ k-1,i for The i-th sampling point: in, represents the state estimate at time k-1, Represents the matrix P k-1 The i-th column vector of the square root; the 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, the i-th new sampling point is 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, and 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 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 and 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 formula (12), (13) and formula (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 As a function of time t, the value of the Lyapunov function at time t is V L (t), then the solution of formula (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