A method for optimizing the feature selection of autonomous navigation sequence images during planetary soft landing

Abstract

The invention discloses an optimization selection method for autonomous navigation sequence image features in a planetary soft landing process, belonging to the field of spacecraft navigation guidance control technology, and comprising the following steps: S1, establishing discrete time state and observation error equations of visual-aided inertial navigation; S2, analyzing the minimum number of observations required for convergence of observable states according to the rank of the observability matrix of the discrete time system; S3, aiming at unknown environments, the invention designs an observability index, and determines the feature points with the highest observability under the condition of limiting the number of feature points; S4, constructing an observation strategy, and determining the switching timing of feature points by using the conclusion of the minimum number of observations described in S2 and the observability index described in S3.

Description

Technical Field

[0001] The invention relates to an autonomous navigation sequence image feature optimization selection method in a planet soft landing process, belonging to the technical field of spacecraft navigation guidance control. Background Art

[0002] In 2020, China and the United States successively launched Mars probes and implemented landing missions. Among them, high-precision autonomous navigation is the key to ensuring the smooth and accurate landing of the lander. The inertial measurement unit (IMU) is the basic autonomous navigation device carried by the probe, but the accumulation of navigation errors caused by integration is a defect of inertial navigation. Therefore, it is often combined with external measurement equipment to improve navigation accuracy.

[0003] In order to achieve precise planetary landings at the level of 10 to 100 meters, NASA and ESA began to study vision-based relative navigation technology at the end of the last century. It is feasible to use the navigation method of the sequence image-assisted inertial measurement unit in the landing process, and high-precision navigation results can be obtained. This autonomous navigation method based on the combination of sequence images and inertial measurement units is usually called visually assisted inertial navigation (VAIN).

[0004] At present, the observability research on vision-assisted navigation has made some progress, but there are still the following shortcomings:

[0005] (1) Existing papers are all aimed at specific observation conditions, and lack discussion on the number of observations of the same feature point. If the observability of different observation schemes can be analyzed, it can guide the selection of better observation schemes in engineering practice and obtain boundary conditions that meet observability requirements. The minimum number of observations required is a necessary condition to ensure the convergence of the navigation system and can be used to guide the design of observation strategies.

[0006] (2) Research on the observability of the VAIN problem mainly focuses on observable properties, which are used to explore the unobservable aspects of the system, but there is less research on the observability used to describe the accuracy of the navigation system. However, existing research methods are only applicable to the pre-flight design database or the post-landing analysis and filtering process, and cannot achieve online evaluation and guide feature point selection during the landing process. During the landing phase, the lander is limited by computing power. If the feature points that contribute more to navigation accuracy can be selected by calculating the observability online, higher navigation accuracy can be obtained while observing as few feature points as possible. Summary of the invention

[0007] The technical problem solved by the present invention is: to overcome the deficiencies of the prior art, systematically analyze the observability of the visually assisted inertial navigation method, and design an observability index and an observation strategy to guide the autonomous planning of feature points in the process.

[0008] The object of the present invention is achieved through the following technical solution: a method for optimizing the selection of features of autonomous navigation sequence images during a planetary soft landing process, comprising the following steps:

[0009] S1. Establish the discrete-time state equation and observation error equation of vision-assisted inertial navigation;

[0010] S2. Determine the minimum number of observations required for the observable state to converge based on the rank of the observability matrix of the discrete-time system;

[0011] S3, set the observability index, and determine the feature point with the highest observability under the condition of limiting the number of feature points;

[0012] S4. Construct an observation strategy, and use the minimum number of observations described in step S2 and the observability index described in step S3 to determine the switching timing of the feature points.

[0013] The specific process of step S1 is as follows:

[0014] 11) Calculate the continuous-time state equation of vision-assisted inertial navigation;

[0015] 12) Calculate the continuous time error state equation of the lander error state;

[0016] 13) Calculate the discrete time state error equation;

[0017] 14) Calculate the observation error equation for discrete time.

[0018] The continuous-time state equation of vision-assisted inertial navigation calculated in step 11) is:

[0019]

[0020] Where: r, v represent the position and velocity vector of the lander in the landing coordinate system; represents the derivatives of r and v; q represents the attitude quaternion from the landing coordinate system to the body coordinate system; C(q) represents the attitude transfer matrix corresponding to q; g represents the gravitational acceleration in the landing system; and Represents the acceleration and angular velocity output by the IMU;

[0021] The continuous time error state equation of the lander error state calculated in step 12) is specifically:

[0022] definition and are the errors of r and v respectively, δθ represents the attitude angle error; the error state vector is expressed as

[0023] Then the continuous time error state equation of the lander error state is expressed as:

[0024]

[0025] Where: express The derivative of ; the system noise is defined as n a 、n ω Respectively represent the acceleration drift noise and angular velocity drift noise of IMU; 0 3×3 ∈R 3 ×3 Represents a 3×3 zero matrix; I3∈R 3×3 represents the 3×3 identity matrix; assuming that n s feature points with known positions, referred to as known points, and n p feature points whose positions are unknown, referred to as unknown points, r pi Represents the i-th feature point p i The coordinates in the landing coordinate system {L}, is a known point, is an unknown point; the position state vector of the unknown point is expressed as The continuous-time state equation for unknown points is

[0026]

[0027] Where: Represents x p The derivative of

[0028] The specific process of calculating the discrete-time state error equation in step 13) is as follows:

[0029] The state to be estimated is a combination of the lander state and the feature point position The discrete-time state error equation from time i to time j is expressed as

[0030]

[0031] Where: the upper right corner subscript i represents the i-th moment, (j, i) represents from the i-th to the j-th moment; the lander state The state transition matrix Where: Δt represents the discrete time interval;

[0032] The specific process of calculating the discrete time observation error equation in step 14) is as follows:

[0033] Let ri represent the feature point p pointed by the lander in the landing system {L} i The vector, measurement quantity b e i Yes i The unit vector after rotating to this system {B}, the observation equation is:

[0034]

[0035] The measurement quantity at time k is the direction vector of each feature point z (k) The estimated value is Then (k) The measurement error is:

[0036]

[0037] The observation error equation in discrete time is:

[0038]

[0039] Where: η (k) represents the measurement noise;

[0040] The minimum number of observations of the observable state described in step S2 is 3 times.

[0041] The observability index described in step S3 is:

[0042]

[0043] Assume that p1 is a known point, p2 is an unknown point, and C is the position of the lander. In the observation profile Cp1p2, α∈[arcsin(z v / ||r||), min(π-β, π / 2)) represents r and r p2 The angle between them, r = -r1, z v represents the height of the lander, and the angle between the vectors r1 and r2 of the observed feature points is β∈(0,π-arcsin(z v / ||r||)), and the formation of a triangle requires α+β<π; δr represents the observability index; represents the position error bound when observing the last feature point; the attitude estimation error bound is P is the covariance matrix of the extended Kalman filter, P(7:9, 7:9) represents the sub-blocks of rows 7 to 9 and columns 7 to 9 in the P matrix; μ represents the expansion parameter.

[0044] The observation strategy described in step S4 is:

[0045] (1) Let α = arcsin(z v / ||r||);

[0046] (2) Solving β *(k) =arg min(δr(β,μ>0)), along β=β *(k) Direction observation feature points, filtering for 3 cycles;

[0047] (3) If the lander lands, exit this method; if the lander does not land, continue to solve β *(k) =argmin(δr(β,μ>0)), compare δr(β *(k) , μ>0) and δr(β (k，k-1) , μ=0); if δr(β *(k) ,μ>0)<δr(β (k，k-1) , μ=0), then execute step (4), otherwise execute step (5);

[0048] (4) Along β = β *(k) Directionally observe feature points, filter for 3 cycles, and return to step (3);

[0049] (5) Observe the old feature points, filter for 1 cycle, and return to step (3).

[0050] Compared with the prior art, the present invention has the following beneficial effects:

[0051] (1) Regarding the observability problem of visual-assisted inertial navigation, existing research results mostly focus on the analysis of unobservable directions, and lack research on the minimum number of observations. The present invention systematically discusses the minimum number of observations required to ensure convergence of the observable state when observing different numbers of feature points. The research results on the minimum number of observations can be used to guide the design of observation strategies, constraining the minimum number of observations of the same feature point when the lander extracts feature points from sequence images as navigation landmarks, thereby ensuring the convergence of the lander state estimate.

[0052] (2) Regarding the observability problem of visual-assisted inertial navigation, most existing research results are only applicable to the pre-flight design database or post-landing analysis and filtering process, which makes it difficult to evaluate online and guide the selection of feature points. The observability index designed by the present invention is only in the form of a single variable elementary function, which is suitable for autonomous planning and can be used to guide the selection of feature points under the condition of limiting the number of feature point observations, thereby improving navigation accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] Figure 1 Schematic diagram of observation geometry.

[0054] Figure 2 Schematic diagram of autonomous navigation simulation during landing.

[0055] Figure 3 Monte Carlo simulation for different numbers of feature points. 3(a) is the position error, 3(b) is the velocity error, and 3(c) is the attitude error.

[0056] Figure 4 are different observability index pairs. 4(a) is the position error, and 4(b) is the time consumed by the optimization algorithm. DETAILED DESCRIPTION

[0057] To make the objectives, technical solutions and advantages of the present invention more clear, the embodiments of the present invention will be further described in detail below with reference to the accompanying drawings.

[0058] Embodiment 1:

[0059] A method for optimizing and selecting features of autonomous navigation sequence images during a planetary soft landing process comprises the following steps:

[0060] (1) Establish the state equation and observation equation of vision-assisted inertial navigation; First, the continuous-time state equation of vision-assisted inertial navigation is:

[0061]

[0062] Where: r, v represent the position and velocity vector of the lander in the landing coordinate system; represents the derivatives of r and v; q represents the attitude quaternion from the landing coordinate system to the body coordinate system; C(q) represents the attitude transfer matrix corresponding to q; g represents the gravitational acceleration in the landing system; and Represents the acceleration and angular velocity output by the IMU;

[0063] definition and are the errors of r and v respectively, and δθ represents the attitude angle error. The error state vector is expressed as

[0064] Therefore, the continuous-time error state equation of the lander error state is expressed as:

[0065]

[0066] Where: express The derivative of ; the system noise is defined as n a 、n ω Respectively represent the acceleration drift noise and angular velocity drift noise of IMU; 0 3×3 ∈R 3×3 Represents a 3×3 zero matrix; I3∈R 3×3 represents the 3×3 identity matrix. Assume that n s feature points with known positions, referred to as known points, and n p feature points whose positions are unknown, referred to as unknown points, r pi Represents the i-th feature point p i The coordinates in the landing coordinate system {L}, is a known point, is an unknown point. The position state vector of the unknown point is expressed as The continuous-time state equation for unknown points is

[0067]

[0068] Where: Represents x p The derivative of .

[0069] The state to be estimated is a combination of the lander state and the feature point position The discrete-time state error equation from time i to time j is expressed as

[0070]

[0071] Where: the upper right corner subscript (i) represents the i-th moment, (j, i) represents from the i-th to the j-th moment; the lander state The state transition matrix

[0072]

[0073] Where: Δt represents the discrete time interval;

[0074]

[0075] Assume r i Indicates the landing system {L} from the lander to the feature point p i The vector, measurement quantity b e i Yes i After rotating to the unit vector of this system {B}, the observation equation is

[0076]

[0077] The measurement quantity at time k is the direction vector of each feature point z (k) The estimated value is Then (k) The measurement error is

[0078]

[0079] The discrete-time observation error equation is

[0080]

[0081] Where: η (k) represents the measurement noise;

[0082]

[0083]

[0084] (2) According to the rank of the observability matrix of the discrete-time system, the minimum number of observations required for the observable state to converge is analyzed.

[0085] Without considering the noise, according to the error equation of the discrete time system:

[0086]

[0087] After N observations starting from any k time, the observability matrix of the system can be expressed as

[0088]

[0089] When the real state x is brought into The system is observable when it is full rank. is full rank, Full rank is equivalent to Full rank.

[0090]

[0091]

[0092] Case 1: n s =0,n p ≥1

[0093] There are four unobservable directions in the system, and it is necessary to determine when the rank of the observability matrix reaches (9+3n p )-4=5+3n p When the number of observations N = 2, After row-column transformation, it becomes

[0094]

[0095] It can be seen The rank condition is not satisfied. When the number of observations increases to N = 3, After row-column transformation, it becomes

[0096]

[0097] It can be seen at this time The dimension of the right null space is 4, corresponding to the 4 unobservable directions of the system, namely the three-axis position and the yaw angle. At the same time, except for these 4 unobservable directions, there are no other unobservable directions. This means that when observing several unknown points, at least 3 consecutive observations can make the system states other than the three-axis position and the yaw angle all observable.

[0098] Case 2: n s =1,n p ≥0

[0099] There is one unobservable direction, and we need to determine when the rank of the observability matrix reaches (9+3n p )-1=8+3n p When the number of observations N = 2, After row-column transformation, it becomes

[0100]

[0101] It can be seen The rank condition is not satisfied. When the number of observations increases to N = 3, After row-column transformation, it becomes

[0102]

[0103] It can be seen at this time The dimension of the right null space is 1, which corresponds to the unobservable direction of the system, namely the yaw angle. At the same time, except for this unobservable direction, there are no other unobservable directions. This means that when observing one known point and any unknown points, at least three consecutive observations can make the state of the system except the yaw angle observable.

[0104] Case 3: n s >1,n p ≥0

[0105] When the lines between all known points are parallel to g, there is one unobservable direction: the yaw angle. It is necessary to determine when the rank of the observability matrix reaches (9+3n p )-1=8+3n p ; When there is a line between known points that is not parallel to g, the whole state is observable. It is necessary to determine when the rank of the observability matrix reaches 9+3n p When the number of observations N = 2, After row-column transformation, it becomes

[0106]

[0107] It can be seen The rank condition is not satisfied. When the number of observations increases to N = 3, After row-column transformation, it becomes

[0108]

[0109] When the lines between all known points are parallel to g, the unobservable direction of the system is the yaw angle. The dimension of the right null space is 1, which corresponds to the unobservable direction of the system. This means that when observing at least 2 known points and any number of unknown points, and the lines between all known points are parallel to g, after at least 3 consecutive observations, the state of the system except the yaw angle can be observed. As mentioned in case 1.4, the distribution of such characteristic points can be ignored. When there are lines between 2 known points that are not parallel to g, the system has no unobservable direction. At this time, It is full rank. This means that when observing at least 2 known points and any number of unknown points, and there is a line between the 2 known points that is not parallel to g, the state of the system can be made fully observable after at least 3 consecutive observations.

[0110] Combining the observability analysis of discrete systems in the above three cases and the observability analysis results of continuous systems in the above four cases, we can draw the following conclusion: the minimum number of observations of the observable state is 3 times.

[0111] (3) For unknown environments, the present invention designs an observability index to guide the lander to autonomously track and photograph the feature points with the highest observability under the condition of limiting the number of feature points.

[0112] The observability of the system can only indicate that the system state converges to the true value, but it cannot measure the size of the estimation error. When computing resources are limited on the lander, the desired goal is to achieve higher navigation accuracy while tracking fewer feature points. Therefore, when the number of feature points is limited, it is necessary to rationally plan and select the observed feature points to improve navigation accuracy. Observability is an important indicator that describes the convergence speed and accuracy of the system. The greater the observability, the higher the navigation accuracy. The present invention will construct the observability of the lander VAIN navigation system.

[0113] First, the following assumptions are made:

[0114] Assumption 1: During feature point planning, all unknown feature points are in the XY plane of the landing system.

[0115] When the landing environment is unknown, the lander needs to autonomously plan the feature points to be observed. When the state of the feature points has three-dimensional coordinates, it is difficult to optimize the selection of feature points in the three-dimensional unknown space; and when the flight altitude of the lander is high, if the height of the feature points is considered to be much smaller than the flight altitude of the lander, the surface of the planet is regarded as a horizontal plane, and the feature points are distributed on the horizontal plane. Therefore, based on assumption 1, the present invention designs a method for planning feature points in a horizontal plane to guide the observation and filtering of the lander in three-dimensional space.

[0116] Based on observability analysis, under landing conditions with unknown environments, if only unknown points are observed, both position and yaw are unobservable. However, if the origin of the landing system {L} is established at the landing point and the landing point is continuously observed, the landing point can be used as a known point in the landing system, and the position can be observable. However, it is impossible to observe two known points in an unknown environment. The following discusses the observation mode of observing only two feature points, that is, on the basis of observing the landing point as a known point p1, an additional unknown point p2 is observed. Under this observation mode, an observability index is designed to plan the feature points to be observed.

[0117] If we only care about the error in the estimated position of the lander, the observation matrix takes the items related to the position state of the lander:

[0118]

[0119] Traditional observability indicators are often functions of the observation matrix. When the estimated value of the lander position and attitude is known, H r1 is fixed, H r2 It is a function of the unknown point position r2. It can be seen that traditional observability indicators are all functions of r2, which means that if the relative position relationship between the unknown point and the lander can be reasonably planned, the observability can be increased, thereby improving the navigation accuracy.

[0120] The position of the lander is represented by C, such as Figure 1 As shown, in the observation profile Cp1p2, α∈[arcsin(z v / ||r||), min(π-β, π / 2)) represents r and r p2 The angle between (r = -r1), z v represents the height of the lander, and the angle between the vectors r1 and r2 of the observed feature points is β∈(0,π-arcsin(z v / ||r||)), and forming a triangle requires α+β<π.

[0121] Depend on Figure 1 It can be seen that the error in the lander position estimation is caused by the attitude estimation error and the position error of p2. The attitude estimation error is bounded by P is the covariance matrix of the extended Kalman filter, and P(7:9, 7:9) represents the sub-blocks of rows 7 to 9 and columns 7 to 9 in the P matrix. The position error bound of p2 first inherits the position estimation error when observing the previous feature point. Taking into account the lander state recursion error and b The measurement noise of e2 brings new errors, and the position error of p2 at time k is bounded by μ represents the dilation parameter.

[0122] The position error δr along the r1 direction is used as the observability index. δr is divided into two parts. The first part is The projection in the r1 direction is expressed as

[0123]

[0124] The second part is the influence of the posture error when observing p2, which is expressed as

[0125]

[0126] In summary, the observability index describing the position estimation error of the lander in the r direction is obtained:

[0127]

[0128] The independent variables of δr are α and β, and δr can be reduced by optimizing α and β. That is, under assumption 1, the observability can be improved by planning the position of p2 on the ground plane.

[0129] (4) An observation strategy was constructed, using the conclusion of the minimum number of observations described in S2 and the observability index described in S3 to guide the autonomous switching of feature points during the landing process.

[0130] First, we discuss the effect of α on observability when β is constant. For any given β, the distance from p2 to the lander is is about α∈[arcsin(z v / ||r||), min(π-β,π / 2)), when α=arcsin(z v / ||r||), the observation section is perpendicular to the ground, and when α=min(π-β,π / 2), the observation section is parallel to the ground. Estimate the feature point coordinates r p2 When, under Assumption 1, it is transformed into estimating only The corresponding observation matrix can be expressed as p2 Rewritten as:

[0131]

[0132] Describing with GDOP Estimated accuracy of:

[0133]

[0134] Where: r 2z =[0 0 1]r2. The more hours The higher the estimation accuracy is, the higher the ground surface is. 2z | is equal to the height of the lander, so the closer the distance ||r2|| from the feature point to the lander is, The higher the estimation accuracy, the better. The estimation accuracy is the highest. According to the analysis in the previous article, we only need to make the observation section perpendicular to the ground, so in Algorithm 1, we take α = arcsin(z v / ||r||).

[0135] When α = arcsin(z v / ||r||), δr is a single variable elementary function about β∈(0,π-α), and the minimum point can be easily solved using optimization tools. From this, we can get the optimization algorithm for the observation angle β that minimizes δr when observing a known point and a new feature point in the observation profile perpendicular to the ground: exist Find the minimum point β of δr * =arg min(δr(β, μ>0)). δr is an elementary function, and the algorithm complexity is O(1). The traditional observability index mentioned above involves complex matrix inversion and singular value calculations, and the algorithm complexity is O(n 3 ). Therefore, the calculation speed may be faster when using the observability proposed in this paper to plan feature points.

[0136] Under the observability index of δr, when μ>0, the β obtained by optimizing δr(β, μ>0) *(k) is the observation angle of the optimal configuration when observing a new feature point; if we continue to track and observe the feature point at the previous moment, then μ = 0, and the corresponding observation angle β (k，k-1) It can be calculated through geometric relationships. At this time, the observability index δr(β (k，k-1) , μ=0). Because the value of parameter μ changes, δr(β (k，k-1) ,μ>0)>δr(β (k，k-1) , μ=0), even if the β *(k) It can make δr(β *(k) ,μ>0)<δr(β (k，k-1) , μ>0), but it may not necessarily make δr(β*(k) ,μ>0)<δr(β (k，k-1) , μ=0). Therefore, in Algorithm 1, by comparing δr(β *(k) , μ>0) and δr(β (k，k-1) , μ=0), we can compare the observability of the new feature point with that of the old feature point. If the observability of the old feature point is large, continue to track and observe; if the new feature point can bring greater observability, switch to β *(k) In order to ensure the convergence of the observable state in Algorithm 1, according to the analysis results of the minimum number of observations, after observing a new feature point, it is necessary to observe continuously for at least 3 filtering cycles.

[0137] Therefore, the observation strategy designed is as follows:

[0138] (1) Let α = arcsin(z v / ||r||);

[0139] (2) Solving β *(k) =arg min(δr(β,μ>0)), along β=β *(k) Direction observation feature points, filtering for 3 cycles;

[0140] (3) If the lander lands, exit this method; if the lander does not land, continue to solve β *(k) =argmin(δr(β,μ>0)), compare δr(β *(k) , μ>0) and δr(β (k，k-1) , μ=0); if δr(β *(k) ,μ>0)<δr(β (k，k-1) , μ=0), then execute step (4), otherwise execute step (5);

[0141] (4) Along β = β *(k) Directionally observe feature points, filter for 3 cycles, and return (3).

[0142] (5) Observe the old feature points, filter for 1 cycle, and return to (3).

[0143] The method of Example 1 is applied, and the landing coordinates are selected as the landing site of Tianwen-1: 109.9°E, 25.1°N. The landing guidance rate of Apollo landing on the moon is used as the landing guidance algorithm. The landing phase lasts 198 seconds. The initial position of the lander in the landing coordinate system is [3300 5870 6570] T m, initial speed is [-56 -66 -90] T m / s, initial attitude is [45 1800] T°. First, verify the observation results under specific initial conditions. When observing the origin of the landing system and an unknown point, assume that the initial position estimation error is 500m, the initial velocity estimation error is 3m / s, the initial attitude angle estimation error is 1°, and the line of sight angle measurement error is 0.01rad (1σ). The landing process is as follows Figure 2 As shown in the figure, the estimated trajectory of the lander gradually converges to the true value. The red asterisks indicate the feature points observed at each filtering. The blue dotted line connects the lander and the feature points, indicating the observation line of sight at each moment. It can be seen that the feature points observed during the landing process are frequently switched under the premise of meeting the observability conditions, which verifies the effectiveness of the observation strategy.

[0144] The following Monte Carlo simulation verification is performed. Assume that the initial position estimation error range is ±1000m, the initial velocity estimation error range is ±10m / s, and the initial attitude estimation error range is ±2°. They obey uniform distribution. After 200 simulations, the 3σ envelope curve of the navigation error is as follows: Figure 3 shown. Figure 3 The navigation error is shown when observing different numbers of feature points.

[0145] from Figure 3 It can be seen that when n s =0 and n p = 1, the position and yaw of the lander do not converge; when n s =1 and n p = 1, the position of the lander can converge, but the yaw still does not converge; when n s =2 and n p =1, all states can converge.

[0146] In order to verify the superiority of algorithm 1 in improving navigation accuracy, algorithm 1 is compared with traditional observability indicators such as error ellipsoid, geometric precision factor and condition number. The 3σ envelope curve of navigation error after 200 simulations is shown in the figure below: Figure 4 (a) shown.

[0147] Figure 4 (a) shows the navigation position error results when using these three traditional observability indicators and the indicators in Algorithm 1 respectively. It can be seen that the position error after using the indicators in Algorithm 1 is smaller than that of other indicators, and the navigation accuracy is higher.

[0148] In addition, in order to verify the rapidity of online application of Algorithm 1 on the lander, the time consumption of optimizing and solving four indicators is compared. The running environment of the simulation program is Matlab R2020a, Intel(R)Core(TM)i7-10750H CPU, 16GBRAM. For the sake of fairness, the four indicators are all used to plan feature points in the observation profile perpendicular to the horizontal plane, so that they are all planning problems of single variable functions, the difference is that their respective objective functions are different. And the golden section method is used to optimize these four planning problems. The distribution of time consumption for each indicator is shown in the figure below. Figure 4 (b) is shown. From top to bottom, each horizontal line represents the maximum value, upper quartile, median, lower quartile, and minimum value. It can be seen that since the observability index in Algorithm 1 has a smaller algorithm complexity than the other three traditional indexes and a higher computational efficiency during optimization, it is significantly better than the other three methods in terms of rapid solution and is easier to apply online on the lander.

[0149] In summary, through the above embodiments, the feasibility and effectiveness of the method for optimizing the selection of features of autonomous navigation sequence images in the planetary soft landing process proposed by the present invention have been verified.

[0150] Although the present invention has been disclosed as above in the form of a preferred embodiment, it is not intended to limit the present invention. Any person skilled in the art may make possible changes and modifications to the technical solution of the present invention by using the methods and technical contents disclosed above without departing from the spirit and scope of the present invention. Therefore, any simple modifications, equivalent changes and modifications made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solution of the present invention shall fall within the protection scope of the technical solution of the present invention.

Claims

1. A method for optimizing the features of autonomous navigation sequence images during a planetary soft landing process, characterized in that: The steps include: S1. Establish the discrete-time state equation and observation error equation of vision-assisted inertial navigation; S2. Determine the minimum number of observations required for the observable state to converge based on the rank of the observability matrix of the discrete-time system; S3, set the observability index, and determine the feature point with the highest observability under the condition of limiting the number of feature points; The observability index described in step S3 is: Assume that p1 is a known point, p2 is an unknown point, and C is the position of the lander. In the observation profile Cp1p2, α∈[arcsin(z v / ||r||), min(π-β,π / 2)) represents r and r p2 The angle between them, r = -r1, z v represents the height of the lander, and the angle between the vectors r1 and r2 of the observed feature points is β∈(0,π-arcsin(z v / ||r||)), and the formation of a triangle requires α+β<π; δr represents the observability index; Represents the position error bound when observing the last feature point; ε is the attitude estimation error bound; μ represents the expansion parameter; S4, constructing an observation strategy, using the minimum number of observations in step S2 and the observability index in step S3 to determine the switching timing of the feature points; The observation strategy described in step S4 is: (1) Let α = arcsin(z v / ||r||); (2) Solving β *(k) =argmin(δr(β,μ>0)), along β=β *(k) Direction observation feature points, filtering for 3 cycles; (3) If the lander lands, exit this method; if the lander does not land, continue to solve β *(k) =argmin(δr(β,μ>0)), compare and δr(β (k,k-1) ,μ=0); if Then execute step (4), otherwise execute step (5); (4) Along β = β *(k) Directionally observe feature points, filter for 3 cycles, and return to step (3); (5) Observe the old feature points, filter for 1 cycle, and return to step (3).

2. The method for optimizing and selecting features of autonomous navigation sequence images during a planetary soft landing process according to claim 1, characterized in that: The specific process of step S1 is as follows: 11) Calculate the continuous-time state equation of vision-assisted inertial navigation; 12) Calculate the continuous time error state equation of the lander error state; 13) Calculate the discrete time state error equation; 14) Calculate the observation error equation for discrete time.

3. The method for optimizing and selecting features of autonomous navigation sequence images during a planetary soft landing process according to claim 2, characterized in that: The continuous-time state equation of vision-assisted inertial navigation calculated in step 11) is: Where: r, v represent the position and velocity vector of the lander in the landing coordinate system; represents the derivatives of r and v; q represents the attitude quaternion from the landing coordinate system to the body coordinate system; C(q) represents the attitude transfer matrix corresponding to q; g represents the gravitational acceleration in the landing system; and Represents the acceleration and angular velocity output by the IMU; 4. The method for optimizing and selecting features of autonomous navigation sequence images during a planetary soft landing process according to claim 3, characterized in that: The continuous time error state equation of the lander error state calculated in step 12) is specifically: definition and are the errors of r and v respectively, δθ represents the attitude angle error; the error state vector is expressed as Then the continuous time error state equation of the lander error state is expressed as: Where: express The derivative of ; the system noise is defined as n a 、n ω Respectively represent the acceleration drift noise and angular velocity drift noise of IMU; 0 3×3 ∈R 3×3 Represents a 3×3 zero matrix; I3∈R 3×3 represents the 3×3 identity matrix; assuming that n s feature points with known positions, referred to as known points, and n p feature points whose positions are unknown, referred to as unknown points, r pi Represents the i-th feature point p i The coordinates in the landing coordinate system {L}, is a known point, is an unknown point; the position state vector of the unknown point is expressed as The continuous-time state equation for unknown points is Where: Represents x p The derivative of .

5. The method for optimizing and selecting features of autonomous navigation sequence images during a planetary soft landing process according to claim 4, characterized in that: The specific process of calculating the discrete-time state error equation in step 13) is as follows: The state to be estimated is a combination of the lander state and the feature point position The discrete-time state error equation from time i to time j is expressed as Where: the upper right corner subscript i represents the i-th moment, (j,i) represents from the i-th to the j-th moment; the lander state The state transition matrix Where: Δt represents the discrete time interval; 6. The method for optimizing and selecting features of autonomous navigation sequence images during a planetary soft landing process according to claim 5, characterized in that: The specific process of calculating the discrete time observation error equation in step 14) is as follows: Assume r i Indicates the landing system {L} from the lander to the feature point p i The vector, measurement quantity b e i Yes i The unit vector after rotating to this system {B}, the observation equation is: The measurement quantity at time k is the direction vector of each feature point z (k) The estimated value is Then (k) The measurement error is: The observation error equation in discrete time is: Where: η (k) represents the measurement noise; 7. The method for optimizing and selecting features of autonomous navigation sequence images during a planetary soft landing process according to claim 1, characterized in that: The minimum number of observations of the observable state described in step S2 is 3 times.

8. The method for optimizing and selecting features of autonomous navigation sequence images during a planetary soft landing process according to claim 1, characterized in that: Said Where P is the covariance matrix of the extended Kalman filter, and P(7:9,7:9) represents the sub-blocks of rows 7 to 9 and columns 7 to 9 in the P matrix.

Images