A Method for Relative Attitude Estimation of a Failed Spacecraft under Computational Resource Constraints

Abstract

The present invention relates to the technical field of spacecraft navigation, and discloses a method for estimating the relative attitude of a failed spacecraft under computational resource constraints, including performing measurement data preprocessing, obtaining camera relative attitude observation data based on an ellipse recognition algorithm, and obtaining lidar relative attitude measurement data based on an iterative closest point algorithm; designing pseudo-measurement factors for the camera and lidar, and designing a variational integral prediction factor for the attitude increment of the failed spacecraft at adjacent moments; performing factor graph optimization based on a moving window to estimate the relative attitude of the failed spacecraft; which is beneficial to exploiting the momentum conservation characteristics of the failed spacecraft through the designed variational integral prediction factor, constructing the constraint relationship of the attitude increment at adjacent moments, and solving the problem that the estimation accuracy of the existing scheme depends on high-frequency measurement input; through the established full-state factor graph model of the failed spacecraft, fusing the low-frequency measurement data of the camera and lidar, and realizing the relative pose estimation of non-cooperative targets under low computing power load.

Description

Technical Field

[0001] The present invention relates to the technical field of spacecraft navigation, and more particularly to a method for estimating the relative attitude of a failed spacecraft under computational resource constraints. Background Art

[0002] The proximity operation of a failed spacecraft is an urgent need to ensure the healthy development of the space industry. Due to the significant non-cooperative characteristics of the failed target, the use of multiple vision sensors for fusion observation can significantly improve the detection and evaluation capabilities, but there are still many difficult challenges in observing and estimating it. On the one hand, the computational performance of the on-board platform is insufficient to meet the high-frequency sampling of the vision sensor, resulting in serious limitations in attitude estimation. On the other hand, the current multi-source fusion method based on Bayesian filtering cannot handle the missing measurement sources, and often ignores the spacecraft dynamics identification, resulting in the need for additional time registration when processing the asynchronous update of multi-measurement information, and cannot solve the problem of estimation divergence caused by missing measurement information. Therefore, under the condition of limited computational performance, designing a multi-vision sensor fusion scheme for the relative attitude estimation of a failed spacecraft has great theoretical and engineering value.

[0003] In the published document with the publication number CN111678512B, a method for determining the attitude of a satellite by combining a star sensor and a gyroscope based on a factor graph is disclosed. An incremental smoothing optimization method based on a factor graph is designed for the attitude determination by fusing the star sensor and the gyroscope, realizing the high-precision attitude determination of the body spacecraft itself. However, the method in this published document only considers its own attitude determination and cannot handle the relative attitude determination of a failed spacecraft lacking information interaction. In the published document with the publication number CN114114311A, a method for measuring the relative pose of a non-cooperative spacecraft based on multi-source information fusion is disclosed. A relative pose measurement method is designed for the fusion of a binocular camera and a lidar, realizing the determination of the relative pose information of a non-cooperative spacecraft. However, the estimation accuracy of the method in this published document is seriously affected by the observation noise, and is extremely dependent on the high-frequency input of the visual measurement information. Facing the complex and changeable space lighting environment and the severely limited on-board computing power resources, the actual engineering implementation accuracy is not high.

[0004] Most of the existing research methods for the problem of multi-sensor fusion estimation of spacecraft are based on the cooperative assumption of complete information interaction, and cannot handle the fusion estimation of multiple vision sensors under low computing power loads, affecting the execution of the proximity operation task of a failed spacecraft. Summary of the Invention

[0005] In order to overcome the above-mentioned defects of the prior art, the present invention provides a method for estimating the relative attitude of a failed spacecraft under computational resource constraints to solve the problems existing in the above-mentioned background art.

[0006] The present invention provides the following technical solution: A method for estimating the relative attitude of a failed spacecraft under computational resource constraints, comprising the following steps:

[0007] Step S01: Perform preprocessing of measurement data, obtain camera relative attitude observation data based on an ellipse recognition algorithm, and obtain lidar relative attitude measurement data based on the iterative closest point algorithm;

[0008] Step S02: Design pseudo-measurement factors for the camera and lidar, and design a variational integral prediction factor for the attitude increment of the failed spacecraft at adjacent moments;

[0009] Step S03: Perform factor graph optimization based on a moving window to estimate the relative attitude of the failed spacecraft.

[0010] Preferably, the ellipse recognition algorithm includes arc segment extraction, ellipse detection, and clustering fusion; first, perform Gaussian filtering on the original image once, then use the Canny operator to extract the image edges, divide the boundary points into positive and negative categories according to the phase of the boundary point gradient, detect the 8-neighborhood connectivity of the boundary points with the same gradient positivity and negativity and connect them into arc segments for extraction;

[0011] Design two indicators to screen the extracted arc segments: Indicator 1: The total length of the arc segment ; Indicator 2: The length of the short side of the arc segment boundary box , combined with the concavity and convexity of the alternative arc segments, the quadrant (Ⅰ, Ⅱ, Ⅲ, Ⅳ) to which they belong can be divided;

[0012] Based on the extracted arc segments, first combine the arc segments in a counterclockwise order of quadrants in pairs, expressed as: ((Ⅰ, Ⅱ), (Ⅱ, Ⅲ), (Ⅲ, Ⅳ), (Ⅳ, Ⅰ)), and eliminate the combinations that do not meet the position constraints; list the possible combinations (Ⅰ, Ⅱ, Ⅲ), (Ⅱ, Ⅲ, Ⅰ), (Ⅲ, Ⅳ, Ⅱ), (Ⅳ, Ⅰ, Ⅲ), and eliminate the combinations that do not meet the concentricity constraints.

[0013] Preferably, the acquisition method of the camera relative attitude measurement data is:

[0014] Perform clustering fusion through two-stage screening;

[0015] The first stage: fuse the ellipses with too close center positions;

[0016] The second stage, check the contour integrity and exclude the ellipses below the set threshold;

[0017] The contour integrity is the ratio of the arc segments detected in the image to the estimated elliptical edge. The distance between the centers of the ellipses to be fused is used as the clustering index, and the ellipse with the highest contour integrity is retained. By finally identifying the five-element parameters of the ellipse, the rotation matrix for affine transformation into a perfect circle is calculated, which is the relative attitude measurement data obtained from the camera image.

[0018] Preferably, the method for obtaining the relative attitude measurement data of the lidar is as follows:

[0019] The attitude measurement data is obtained based on the iterative closest point algorithm. Based on two point sets with a rotational correspondence relationship, the following optimization loss function is constructed: , where R represents the relative attitude, and E(R) represents the loss function of the relative attitude. represents the relative attitude input parameter for which the loss function obtains the minimum value; U is the known point set of the three-dimensional model of the failed spacecraft, V is the lidar measurement point set, F represents the Frobenius norm of the matrix, "s.t." indicates that the variables in the objective function on its left satisfy the constraints on the right, and I3 is the third-order identity matrix; , where n represents the number of feature points, and m = 1, 2, 3,..., n; u m is the m-th three-dimensional feature point in the known three-dimensional model of the failed spacecraft, and v m is the m-th three-dimensional feature point measured by the lidar;

[0020] The rotation matrix that minimizes the loss function after iterative optimization is the relative attitude measurement data obtained by the lidar.

[0021] Preferably, in step S02, the pseudo-measurements of the camera and the lidar are expressed as:

[0022] ; ;

[0023] Among them, and respectively represent the observations of the camera and the lidar at time k; and respectively represent the noises of the two vision sensors; and respectively represent the observation functions of the two sensors; x k represents the motion state of the failed spacecraft at time k;

[0024] The corresponding pseudo-measurement odometer factor is expressed as:

[0025] ; ;

[0026] Among them, denotes the cost function; f cam denotes the pseudo-measurement odometry factor corresponding to the camera, f LiDAR denotes the pseudo-measurement odometry factor corresponding to the lidar.

[0027] Preferably, if the measurement data satisfies the Gaussian distribution, the cost function of the odometry factor is expressed as:

[0028] ;

[0029] where is the error function, expressed as the error between the sensor measurement value and the state prediction value, i.e., ; is the Mahalanobis distance of the error function, is the measurement noise variance matrix of the corresponding sensor, denotes the mileage.

[0030] Preferably, in the step S02, the variational integration prediction factor is obtained as follows:

[0031] Apply variational integration to the rigid body rotation Euler equation, where is the rotation increment at the current moment and

[0032] ;

[0033] where denotes the vector part of the attitude quaternion at time k; the rotation increment is defined as ; denotes quaternion multiplication, is 's conjugate quaternion; denotes the moment of inertia of the failed spacecraft; denotes the cross product matrix of three-dimensional vectors;

[0034] Let , then:

[0035] ;

[0036] Design the variational integration prediction factor as:

[0037] .

[0038] Preferably, in the step S03, the factor acquisition process in the factor graph is as follows:

[0039] Solve the maximum a posteriori probability distribution of the state quantity according to the observed quantity:

[0040] ;

[0041] Among them, represents the posterior probability of the state quantity determined according to the observed quantity; represents the joint probability of the observed quantity and the state quantity; represents the prior probability of the state quantity before there is no measurement quantity information; represents the probability of the occurrence of the measurement quantity, which is also called the marginal probability; represents the probability of the measurement quantity in the case where the state quantity has occurred, which is also called the likelihood function;

[0042] ; ;

[0043] Among them, represents the state quantity in the case where it has occurred, the measurement quantity the likelihood probability of occurrence; represents the state from the previous moment changing to the state at the next moment the transition probability; represents the prior probability of the initial state ;

[0044] ; Among them represents the state quantity in the case where it has occurred, the measurement quantity the likelihood probability of occurrence, represents the measurement quantity other than the IMU in the whole system; z j represents the measurement information obtained through the camera or lidar and processed into an odometer; Therefore, the posterior probability of the state inferred from the observed quantity is:

[0045] ;

[0046] (1);

[0047] Among them, represents the IMU measurement information at time; Each term in the above formula corresponds to a factor in the factor graph.

[0048] Preferably, in the step S03, the specific method for optimizing the factor graph is:

[0049] After the factor graph construction is completed, the optimization of the entire graph is carried out; the states and constraints at all times are transformed into a least squares problem, and the final optimized estimation result is obtained through iterative solution by a numerical optimization method; to avoid the huge computational burden caused by the synchronous optimization of all states, the sliding window method is used to limit the optimization to time to At a moment, the state from the initial moment to this moment is regarded as a priori factor and added to the figure;

[0050] Taking the negative logarithm of the posterior probability in Equation (1), the optimal state solution is:

[0051] ;

[0052] Wherein, represents the optimal state that minimizes the cost function.

[0053] Preferably, in step S03, the relative attitude of the failed spacecraft is estimated by obtaining the minimum value of the following formula, and the formula is expressed as:

[0054] .

[0055] The technical effects and advantages of the present invention:

[0056] (1) By providing step S01, the present invention is conducive to effectively improving the robustness of visual attitude recognition through the established method based on arc segment detection, through three stages of arc segment extraction, ellipse detection and clustering fusion.

[0057] (2) By providing step S02, the present invention is conducive to mining the momentum conservation characteristics of the failed spacecraft by designing the variational integral prediction factor, constructing the constraint relationship of the attitude increment between adjacent moments, and solving the problem that the estimation accuracy of the existing scheme depends on high-frequency measurement input.

[0058] (3) By providing step S03, the present invention is conducive to realizing the relative pose estimation of non-cooperative targets under low computing power load by integrating the low-frequency measurement data of cameras and lidar through the established full-state factor graph model of the failed spacecraft. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 is a flowchart of the method for estimating the relative attitude of a failed spacecraft under the computing resource constraint of the present invention.

[0060] Figure 2 is a flowchart of the ellipse recognition algorithm of the present invention.

[0061] Figure 3 is the full-state factor graph model of the relative navigation of the failed spacecraft of the present invention.

[0062] Figure 4 is a curve graph of the attitude estimation error. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0063] The technical solutions of the present invention will be clearly and completely described below in conjunction with the accompanying drawings of the present invention. In addition, the forms of the various structures described in the following embodiments are merely examples, and a method for estimating the relative attitude of a failed spacecraft under computational resource constraints according to the present invention is not limited to the various structures described in the following embodiments. All other embodiments obtained by those of ordinary skill in the art without creative efforts fall within the scope of protection of the present invention.

[0064] As Figure 1 shown, the present invention provides a method for estimating the relative attitude of a failed spacecraft under computational resource constraints, including the following steps:

[0065] Step S01: Perform preprocessing of measurement data, obtain camera relative attitude observation data based on an ellipse recognition algorithm, and obtain lidar relative attitude measurement data based on an iterative closest point algorithm; the measurement data is the relative attitude measurement data;

[0066] The visual camera carried on the servicing spacecraft is used to photograph the circular features on the surface of the failed spacecraft, detect the ellipse parameters from a single-frame image, and calculate the relative attitude measurement data during the process of affine-transforming the ellipse into a perfect circle; the circular features include, but are not limited to, thruster nozzles and docking rings, etc., and the ellipse parameters include, but are not limited to, the centroid pixel coordinates, semi-major axis, semi-minor axis, and ellipse inclination angle; the relative attitude measurement data includes camera relative attitude measurement data and lidar relative attitude measurement data;

[0067] Step S02: Design pseudo-measurement factors for the camera and lidar, and design a variational integral prediction factor for the attitude increment of the failed spacecraft at adjacent moments; design odometer factors for the camera image information and lidar point cloud information respectively, and then add the pseudo-measurement odometer factors to the full-state factor graph;

[0068] Step S03: Define all states from the initial moment to the k-th moment, and define the camera, lidar, and variational integral predictor as observations of the system state; perform factor graph optimization based on a moving window to effectively reduce the computational resource requirements and achieve accurate estimation of the relative attitude of the failed spacecraft during proximity operations.

[0069] In this embodiment, it should be specifically noted that the method for obtaining the ellipse parameters is specifically as follows:

[0070] As Figure 2As shown in the figure, an ellipse recognition algorithm is used, including arc segment extraction, ellipse detection and cluster fusion; first, a Gaussian filter is performed on the original image, and then the Canny operator is used to extract the image edge. According to the phase of the boundary point gradient, it can be divided into positive and negative categories. Then, the connectivity of the 8-neighborhood boundary points with the same gradient positivity is detected and connected into arc segments for extraction; two indicators are designed to screen the extracted arc segments: Indicator 1: Total length of the arc segment , this indicator is to filter out false detections caused by noise; Indicator 2: The short side length of the arc segment bounding box , this indicator ensures that the arc segment has enough curvature to form an ellipse, rather than a straight line segment; combined with the concavity and convexity of the candidate arc segment, the quadrant to which it belongs can be divided (Ⅰ, Ⅱ, Ⅲ, Ⅳ);

[0071] Based on the extracted arc segments, first combine the arc segments in pairs in the counterclockwise order of the quadrants, expressed as: ((Ⅰ, Ⅱ), (Ⅱ, Ⅲ), (Ⅲ, Ⅳ), (Ⅳ, Ⅰ)), and eliminate the combinations that do not meet the position constraints; the position constraints are expressed as the constraints on the arc segment positions; taking the combination (Ⅰ, Ⅱ) as an example, the right end of arc segment Ⅱ must be located to the left of the left end of arc segment Ⅰ; taking the combination (Ⅱ, Ⅲ) as an example, the lower end of arc segment Ⅱ must be located above the upper end of arc segment Ⅲ; on this basis, the possible Combinations (Ⅰ, Ⅱ, Ⅲ), (Ⅱ, Ⅲ, Ⅰ), (Ⅲ, Ⅳ, Ⅱ), (Ⅳ, Ⅰ, Ⅲ), and eliminate combinations that do not satisfy the concentricity constraint; the concentricity constraint utilizes the property that the trajectory of the midpoints of the parallel chords of an ellipse passes through the center of a circle. Two arc segments can determine a center of a circle, so the position error of the two centers determined by three arc segments must be within a certain range, and then the five parameters of the ellipse are estimated.

[0072] In this embodiment, it should be specifically explained that the method for obtaining the camera relative posture measurement data is:

[0073] Cluster fusion is performed through two-stage screening;

[0074] Stage 1: Merge ellipses whose centers are too close to each other;

[0075] In the second stage, the contour completeness is checked and ellipses below the set threshold are excluded;

[0076] The contour completeness is the ratio of the arc segments detected in the image to the edge of the estimated ellipse. The lower the ratio, the higher the probability of mismatch. The distance between the centers of the ellipses to be fused is used as the clustering index, and the ellipse with the highest contour completeness is retained. By finally identifying the five element parameters of the ellipse, the rotation matrix of its affine transformation into a perfect circle is calculated, which is the relative posture measurement data obtained by the image.

[0077] In this embodiment, it should be specifically explained that the laser radar relative attitude measurement data is obtained in the following manner:

[0078] The point cloud data of the failed spacecraft is obtained by the lidar carried on the service spacecraft, and the attitude measurement data is obtained based on the iterative closest point algorithm; based on two point sets with a rotational correspondence relationship, the following optimization loss function is constructed: , where R represents the relative attitude, and E(R) represents the loss function of the relative attitude, represents the relative attitude input parameter for obtaining the minimum value of the loss function; U is the known point set of the three-dimensional model of the failed spacecraft, V is the lidar measurement point set, F represents the Frobenius norm of the matrix, "s.t." means that the variables in the objective function on its left satisfy the constraints on the right, and I3 is the third-order identity matrix; , where n represents the number of feature points, m = 1, 2, 3,..., n; u m is the m-th three-dimensional feature point in the known three-dimensional model of the failed spacecraft, and v m is the m-th three-dimensional feature point measured by the lidar;

[0079] The rotation matrix that minimizes the loss function after iterative optimization is the relative attitude measurement data obtained by the lidar.

[0080] In this embodiment, it should be specifically noted that the pseudo-measurements of the camera and the lidar in step S02 are expressed as:

[0081] ; ;

[0082] Among them, , respectively represent the observations of the camera and the lidar at time k; , respectively represent the noises of the two vision sensors; , respectively represent the observation functions of the two sensors; x k represents the motion state of the failed spacecraft at time k;

[0083] The corresponding pseudo-measurement odometry factor is expressed as:

[0084] ; ;

[0085] Among them, represents the cost function; f cam represents the pseudo-measurement odometry factor corresponding to the camera, and f LiDAR represents the pseudo-measurement odometry factor corresponding to the lidar; if the measurement data all satisfy the Gaussian distribution, the cost function of the odometry factor is expressed as:

[0086] ;

[0087] Among them, is the error function, which is expressed as the error between the sensor measurement value and the state prediction value, that is ; is the Mahalanobis distance of the error function, is the measurement noise variance matrix of the corresponding sensor, represents the mileage.

[0088] In this embodiment, it should be specifically noted that in step S03, considering the relative navigation and positioning problem of the failed spacecraft, all states from the initial time to time k are defined as , and the camera, lidar, and variational integration predictor are all defined as observations of the system state ;

[0089] The specific method of variational integration prediction is as follows:

[0090] Considering that the failed spacecraft is freely tumbling in space, without considering the external torque and disturbance torque, and conserving angular momentum and momentum during the proximity operation of the servicing spacecraft; applying variational integration to the rigid body rotation Euler equation, the rotation increment at the current time and the rotation increment at the next time are used to construct an equality constraint:

[0091] ;

[0092] Among them, represents the vector part of the attitude quaternion at time k; the rotation increment is defined as , represents quaternion multiplication, is 's conjugate quaternion; represents the moment of inertia of the failed spacecraft; represents the cross product matrix of three-dimensional vectors;

[0093] Let , then:

[0094] ;

[0095] The variational integration prediction factor is designed as:

[0096] ;

[0097] When using the numerical optimization method to solve the variable that minimizes the error function, the Jacobian matrix still needs to be calculated. To further reduce the computational complexity, the Jacobian matrix is given in an analytical form:

[0098] .

[0099] In this embodiment, it should be specifically noted that the state estimation problem is actually to solve the maximum a posteriori probability distribution of the state quantity according to the observed quantity:

[0100] ;

[0101] Among them, represents the posterior probability of the state quantity determined according to the observed quantity; represents the joint probability of the observed quantity and the state quantity; represents the prior probability of the state quantity before there is no measurement information; represents the probability of the occurrence of the measurement quantity, which is also called the marginal probability; represents the probability of the measurement quantity when the state quantity has occurred, which is also called the likelihood function;

[0102] ; ;

[0103] Among them, represents the state quantity when the measurement quantity has occurred; represents the transition probability from the state at the previous moment to the state at the next moment; represents the prior probability of the initial state ;

[0104] ; Among them represents the state quantity when the measurement quantity has occurred, represents the measurement quantity other than the IMU in the whole system; z j represents the measurement information obtained by the camera or lidar and processed into an odometer; Therefore, the posterior probability of the state inferred from the observed quantity is:

[0105] ;

[0106] (1);

[0107] Among them, represents the IMU measurement information at the Figure 3 moment; Each term in the above formula corresponds to Figure 3As shown; after the factor graph construction is completed, further optimization of the entire graph is carried out, which is also called smoothing; the states and constraints at all times are transformed into a least squares problem, and iterative solutions are obtained through numerical optimization methods such as the Gauss-Newton method and the LM method to finally obtain the optimized estimation result; to avoid the huge computational burden caused by the synchronous optimization of all states, a sliding window method is adopted to limit the optimization to time to time, and the states from the initial time to time are regarded as prior factors and added to the graph;

[0108] Taking the negative logarithm of the posterior probability in Equation (1), the optimal state solution is:

[0109] ;

[0110] where represents the optimal state that minimizes the cost function; the relative attitude estimation problem of the final failed spacecraft is described as the following minimum value solving problem: ;

[0111] Finally, a relative attitude estimation experiment is carried out on the ground dual-manipulator semi-physical platform. The camera and lidar are placed on the service spacecraft, and relative measurement data of the target spacecraft are collected during the proximity operation. In addition, eight Vicon infrared cameras are installed on the periphery of the experimental scene, and together with the retroreflective target balls placed on the surface of the spacecraft in advance, accurate true values of the relative attitude are obtained. The initial values of the experimental scene are set as follows: the target attitude quaternion [-0.0041, -0.7067, -0.0053, 0.7075], the target relative position [7.03316, 0.03439, 1.90147] m, the relative angular velocity [0.5, 1.5, 1] * 10-3 rad / s, and the relative linear velocity [-1, 1, -1] * 10-2 m / s. The moment of inertia of the target spacecraft is set as follows:

[0112] ;

[0113] The ellipse of the recognized circular nozzle is converted into a perfect circle through affine transformation, so as to obtain the attitude measurement value of the camera. The image processing time for each frame is about 250 milliseconds, so the camera measurement frequency is 4 Hz. Although the acquisition frequency of repeated scans can reach tens of thousands of Hz, considering the computational efficiency of ICP, the actual LiDAR measurement frequency is 0.5 Hz.

[0114] By fusing the measurement data of the camera and the lidar, the attitude estimation error is as Figure 4As shown in the figure. The present invention uses the classical multiplicative extended Kalman filter method (MEKF) in the field of spacecraft attitude estimation for comparison to highlight the accuracy of the estimation method proposed by the present invention. From Figure 4 It can be seen that the average attitude estimation error of the proposed method is 0.2264 degrees. Compared with the steady-state estimation error of MEKF which is 1.0829 degrees, the attitude estimation accuracy has been significantly improved. In addition, all filtering methods require a certain convergence time to ensure good estimation performance. The attitude estimation error of MEKF stabilizes at about 1 degree after 80 seconds. In contrast, the proposed factor-graph-based fusion estimation method adopts global optimization, does not require a convergence time, and maintains a stable and low estimation error throughout the proximity operation estimation process.

[0115] Finally: The above are only the preferred embodiments of the present invention and are not used to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

[0116] The above is only the specific implementation manner of the present application, but the protection scope of the present application is not limited thereto. Any person skilled in the art can easily think of changes or replacements within the technical scope disclosed in the present application, and all should be covered within the protection scope of the present application. Therefore, the protection scope of the present application shall be subject to the protection scope of the claims.

Claims

1. A method for estimating the relative attitude of a failed spacecraft under computational resource constraints, characterized in that: It includes the following steps: Step S01: Perform preprocessing of measurement data, obtain camera relative pose observation data based on the ellipse recognition algorithm, and obtain lidar relative pose measurement data based on the iterative closest point algorithm; Step S02: Design pseudo-measurement factors for the camera and lidar, and design a variational integral prediction factor for the attitude increment of the failed spacecraft at adjacent moments; Step S03: Perform factor graph optimization based on a moving window to estimate the relative pose of the failed spacecraft; In the said step S02, the pseudo-measurements of the camera and lidar are expressed as: ； ； Among them, and respectively represent the observations of the camera and the lidar at time k; and respectively represent the noises of the two vision sensors; and respectively represent the observation functions of the two sensors; x k represents the motion state of the failed spacecraft at time k; The corresponding pseudo-measurement odometer factor is expressed as: ； ； Among them, represents the cost function; f cam represents the pseudo-measurement odometry factor corresponding to the camera, f LiDAR represents the pseudo-measurement odometry factor corresponding to the lidar; If the measurement data satisfies the Gaussian distribution, the cost function of the odometer factor is expressed as: ； Among them, is the error function, expressed as the error between the sensor measurement value and the state prediction value, that is ; is the Mahalanobis distance of the error function, is the measurement noise variance matrix of the corresponding sensor, represents mileage; In the said step S02, the acquisition method of the variational integral prediction factor is: Applying variational integration to the Euler equations of rigid body rotation for the rotation increment at the current moment and the rotation increment at the next moment Construct an equality constraint: ； Among them, represents the vector part of the attitude quaternion at time k; the rotation increment is defined as , represents quaternion multiplication, is the conjugate quaternion of; represents the moment of inertia of the failed spacecraft; represents the cross product matrix of three-dimensional vectors; Let , then: ； Design the variational integral prediction factor as: 。 2. A method for estimating the relative attitude of a failed spacecraft under computational resource constraints according to claim 1, characterized in that: The said ellipse recognition algorithm includes arc segment extraction, ellipse detection, and clustering fusion; first perform a Gaussian filter on the original image, then use the Canny operator to extract the image edges, divide them into positive and negative categories according to the phase of the boundary point gradient, detect the 8-neighborhood connectivity of the boundary points with the same gradient positivity and negativity and connect them into arc segments for extraction; Design two metrics to screen the extracted arcs: Metric 1: Total arc length ; Metric 2: Short side length of the arc bounding box , combined with the concavity and convexity of the candidate arcs, the quadrant (Ⅰ, Ⅱ, Ⅲ, Ⅳ) to which they belong can be determined; Based on the extracted arcs, first combine the arcs pairwise in counterclockwise order of quadrants, expressed as: ((Ⅰ, Ⅱ), (Ⅱ, Ⅲ), (Ⅲ, Ⅳ), (Ⅳ, Ⅰ)), and eliminate the combinations that do not meet the position constraints; list the possible combinations (Ⅰ, Ⅱ, Ⅲ), (Ⅱ, Ⅲ, Ⅰ), (Ⅲ, Ⅳ, Ⅱ), (Ⅳ, Ⅰ, Ⅲ), and eliminate the combinations that do not meet the concentricity constraints.

3. A method for estimating the relative attitude of a failed spacecraft under computing resource constraints according to claim 2, characterized in that: The acquisition method of the camera relative pose measurement data is: Perform clustering fusion through two-stage screening; The first stage: fuse the ellipses with too close center positions; The second stage: check the contour integrity and exclude the ellipses below the set threshold; The said contour integrity is the ratio of the arc segments detected in the image to the estimated ellipse edge. Use the distance between the centers of the ellipses to be fused as the clustering index, and retain the ellipse with the highest contour integrity; by finally identifying the five-element parameters of the ellipse, calculate the rotation matrix for its affine transformation into a perfect circle, which is the relative pose measurement data obtained from the camera image.

4. A method for estimating the relative attitude of a failed spacecraft under computing resource constraints according to claim 3, characterized in that: The acquisition method of the lidar relative pose measurement data is: Obtain attitude measurement data based on the Iterative Closest Point (ICP) algorithm; based on two point sets with a rotational correspondence relationship, construct the following optimization loss function: , where R represents the relative attitude, and E(R) represents the loss function of the relative attitude, represents the relative attitude input parameter for which the loss function obtains the minimum value; U is the known point set of the three-dimensional model of the failed spacecraft, V is the point set measured by the lidar, F represents the Frobenius norm of the matrix, "s.t." indicates that the variables in the objective function on the left side satisfy the constraints on the right side, and I3 is the third-order identity matrix; , where n represents the number of feature points, and m = 1, 2, 3, …, n; u m is the m-th three-dimensional feature point in the known three-dimensional model of the failed spacecraft, and v m is the m-th three-dimensional feature point measured by the lidar; The rotation matrix that minimizes the loss function after iterative optimization is the relative pose measurement data obtained from the lidar.

5. A method for estimating the relative attitude of a failed spacecraft under computational resource constraints according to claim 4, characterized in that: In the said step S03, the factor acquisition process in the factor graph is: Solve the maximum a posteriori probability distribution of the state quantity according to the observed quantity: ； Among them, represents the posterior probability of the state quantity determined according to the observed quantity; represents the joint probability of the observed quantity and the state quantity; represents the prior probability of the state quantity before there is no measurement quantity information; represents the probability of the occurrence of the measurement quantity, which is also called the marginal probability; represents the probability of the measurement quantity when the state quantity has occurred, which is also called the likelihood function; ； ； Among them, represents the state quantity the measured quantity when it has occurred the likelihood probability of occurrence; represents the state from the previous moment changing to the state at the next moment the transition probability; represents the initial state the prior probability; ; where represents a state quantity the measurement quantity when it has occurred the likelihood probability of occurrence, represents the measurement quantity in the entire system except for the IMU; z j represents the measurement information obtained by the camera or lidar and processed as an odometer; therefore, the a posteriori probability of the state inferred from the observed quantity is: ； （1）； Among them, represents the IMU measurement information at the moment; each term in the above formula corresponds to a factor in the factor graph.

6. A method for estimating the relative attitude of a failed spacecraft under computational resource constraints according to claim 5, characterized in that: In the said step S03, the specific method of factor graph optimization is: After the factor graph construction is completed, the entire graph is optimized; the states and constraints at all times are transformed into a least squares problem, and the final optimized estimation result is iteratively solved through numerical optimization methods; to avoid the huge computational burden caused by the synchronous optimization of all states, a sliding window method is used to limit the optimization to time to time, and the states from the initial time to time are regarded as prior factors and added to the graph; Take the negative logarithm of the posterior probability in formula (1), and the optimal state solution is: ； Among them, represents the optimal state that minimizes the cost function.

7. A method for estimating the relative attitude of a failed spacecraft under computational resource constraints according to claim 6, characterized in that: In the said step S03, the relative pose of the failed spacecraft is estimated by obtaining the minimum value of the following formula, and the formula is expressed as: 。

Images