Vision-based autonomous orbit determination method for lunar aircraft

Abstract

The invention provides a vision-based autonomous orbit determination method for a lunar aircraft, and belongs to the technical field of spacecraft orbit determination. The method comprises the following steps: firstly, extracting a centroid point of a lunar surface annular mountain according to a lunar surface image and data, and constructing a geometric information database for registration; secondly, setting a preset orbit of the aircraft, adjusting a visual axis direction of the aircraft and a camera focal length of the aircraft, photographing the moon, registering a photographed image with a geometric information database to find corresponding real annular mountain information, and calculating a current attitude of the aircraft by using a camera attitude estimation method; and thirdly, filtering is carried out on a resolving result, so that orbit determination is autonomously completed in real time. According to the invention, the dependence on ground measurement and control systems such as a ground station can be reduced, and the autonomy and real-time performance of orbit determination when the lunar aircraft executes a daily navigation task are improved; meanwhile, orbit determination can be carried out by shooting only through light rays without emitting energy, so that the observation range is enlarged; and the reliability is high.

Description

Technical Field

[0001] The present invention belongs to the technical field of spacecraft orbit determination, relates to vision principles and camera pose estimation related technologies, and specifically relates to a vision-based autonomous orbit determination method for a lunar orbiting vehicle. Background Art

[0002] For lunar orbiters (such as lunar relay satellites and lunar orbiters), spacecraft orbit determination technology is crucial. Through orbit determination technology, the orbit of the lunar orbiter can be accurately designed to ensure that it covers the key areas of the moon and is in the best position through real-time adjustments, ensuring the continuity of the communication link and the efficiency of data transmission when performing lunar exploration missions. Spacecraft orbit determination technology has high requirements for accuracy and real-time performance, and should also consider multiple factors, such as the autonomy and reliability of the overall system and technical costs.

[0003] Traditional orbit determination technology for lunar spacecraft mainly relies on radio technology and ground measurement and control systems. The ground station on the earth processes radio tracking data to obtain the relative speed and position of the spacecraft relative to the ground station on the earth, and then guides the spacecraft to achieve orbit determination. After searching the existing literature, it was found that Jin-chao Xia, Ge-shi Tang, Yan-long Bu and Chun-lai Li et al. published "Lunar orbiter automotive navigation based on camera and laser altimeter" in the Proceedings of 2014 IEEE Chinese Guidance, Navigation and Control Conference in 2014. They analyzed the current orbit determination technology for lunar spacecraft and believed that although the traditional orbit determination technology relying on the ground measurement and control system still plays an important role, the traditional orbit determination technology has extremely high orbit determination cost, limited observation arc and unmeasurable back of the moon, and the autonomy and real-time performance need to be improved. Summary of the invention

[0004] The technical problem to be solved by the present invention is that when a lunar orbiter uses traditional orbit determination technology to determine its orbit, it mainly relies on the ground measurement and control system, resulting in high overall orbit determination costs, limited observation range, and the need for further improvement in autonomy and real-time performance, which has a huge impact on the lunar orbiter's lunar exploration mission. The main application scenario of the present invention is the flight process of a lunar orbiter performing daily navigation tasks.

[0005] In order to achieve the above object, the technical solution adopted by the present invention is:

[0006] A vision-based autonomous orbit determination method for a lunar circumnavigation vehicle. The method first uses the rich lunar surface images and data collected from multiple lunar exploration missions in the past to extract the centroid points of lunar craters and construct a geometric information database for registration. Secondly, the preset orbit of the aircraft is set, the direction of the aircraft's visual axis and the focal length of the aircraft's camera are adjusted, and simulated photos of the moon are taken. The corresponding real crater information is found by registering the captured images with the geometric information database, and the current aircraft's posture is calculated using camera posture estimation related technologies. Thirdly, the calculated results are filtered to complete the orbit determination autonomously and in real time.

[0007] Assume that a lunar world coordinate system with the center of the moon as the origin is established, the real geographical location coordinates of the main craters are recorded through the existing lunar surface images, different craters are numbered and the centroids of these craters are extracted. At the same time, a preset orbit of the aircraft is set, wherein the preset orbit is an elliptical orbit and is located in the XOY plane of the lunar world coordinate system, the center of the moon (i.e., the origin of the lunar world coordinate system) is located at the left focus of the preset orbit, and the initial position of the aircraft is the perimoon point of the preset orbit (i.e., the left vertex of the ellipse), then the vision-based autonomous orbit determination method of the lunar circumlunar aircraft is divided into three parts: a geometric information database establishment part, a lunar surface photography and result registration part, and a calculation and filtering part of the aircraft posture. The innovation of the present invention is to propose a new autonomous orbit determination method for the lunar circumlunar aircraft that combines visual principles with a camera posture estimation algorithm. The focus is on a method proposed in step 4 to match the lunar surface image taken by the lunar circumlunar aircraft with the actual lunar craters, and in step 5, the camera posture estimation algorithm is used to solve the posture of the aircraft through the two-dimensional geometric information of the image and the three-dimensional real geographic information. The method for autonomous orbit determination of a lunar orbiting vehicle comprises the following steps:

[0008] Part 1: Geometric information database establishment part;

[0009] Step 1: Store the numbering information of adjacent crater groups;

[0010] According to the real geographical location coordinates, calculate the N nearest craters to each crater except itself, record and store the numbers of these adjacent crater groups, then each crater corresponds to a numbering matrix containing a total of N+1 elements, and the first element of this numbering matrix is ​​its own number.

[0011] Step 2: Calculate and store the geometric information of each crater group;

[0012] Step 2.1, first, obtain the coordinates of the centroid of the two-dimensional image corresponding to each adjacent crater group. If the adjacent crater number matrix corresponding to the crater numbered k0 is [k0, k1, ..., k N ], assuming that the corresponding coordinates of the centroid of the two-dimensional image are in Indicates number k n The coordinates of the centroid of the two-dimensional image corresponding to the crater, represents the horizontal coordinate of the centroid of this two-dimensional image, Indicates the vertical coordinate of the centroid of this two-dimensional image. Since the two-dimensional plane coordinate system takes the upper left corner as the origin, for the convenience of calculation, the origin of the two-dimensional plane coordinate system is changed to At, the coordinates of the remaining two-dimensional image centroid points are adjusted to:

[0013]

[0014] in, Indicates number k n The coordinates of the centroid point of the two-dimensional image corresponding to the crater in the adjusted two-dimensional plane coordinate system (hereinafter referred to as the two-dimensional point).

[0015] Step 2.2, calculate the angle difference matrix of this adjacent crater group Modulus ratio matrix Ratio matrix to radius Specific:

[0016] First calculate the angle difference matrix, assuming but:

[0017]

[0018] in, Represented as a two-dimensional point The horizontal axis of Represented as a two-dimensional point The vertical coordinate of Represents a two-dimensional point The angle between the vector formed with the origin and the horizontal axis of the coordinate system; Represents a two-dimensional point The angle between the vector formed by the origin and the horizontal axis of the coordinate system; Δθ mn express and The angle difference between them; m means number k m The sequence of the crater numbered in the adjacent crater number matrix corresponding to the crater numbered k0; n represents the crater numbered k n The sequence of the crater in the adjacent crater number matrix corresponding to the crater numbered k0.

[0019] The angle difference matrix of point k0 is calculated as follows: Assume the number is k n The radius of the crater is Then the module length ratio matrix Ratio matrix to radius The calculation method is as follows:

[0020]

[0021] in, Represents a two-dimensional point The distance from the origin; Represents a two-dimensional point The distance from the origin; Represents a two-dimensional point and 2D points The ratio of the distance to the origin; Indicates number k m The radius of the crater; Indicates number k n The radius of the crater; Indicates that the number is k n The crater numbered k m ratio of the radius of the crater.

[0022] Then get

[0023] The three geometric information matrices (angle difference matrix, modulus ratio matrix, radius ratio matrix) of each centroid point are calculated and stored respectively, so as to establish a geometric information database.

[0024] Part II: lunar surface photography and result registration;

[0025] Step 3: Obtain real-time lunar images and solve related geometric information;

[0026] Assuming that the spacecraft is flying on a preset orbit, adjust the visual axis direction of the spacecraft camera to face the center of the moon, make the optical center of the camera on the spacecraft consistent with the projection point of the center of the moon on the two-dimensional pixel plane, and adjust the focal length so that the moon projection can be tangent to the short side of the image taken by the camera.

[0027] Take a picture with the camera and extract the centroid of the crater, and obtain the coordinates of the centroid of the crater closest to the center of the image (i.e., the optical center), marked as in represents the coordinates of the 2D point with sequence q in the center point group of the crater closest to the center of the image; Represents a two-dimensional point The horizontal axis of Represents a two-dimensional point The vertical coordinate of Sequence in the cluster of crater centroids closest to the image center.

[0028] The two-dimensional point coordinates after the origin is changed are generated by formula (1) The angle difference matrix is ​​obtained by formula (2) and formula (3): The modulus ratio matrix is ​​obtained from formula (4) and formula (5): The radius ratio matrix is ​​obtained from formula (6): After obtaining the three geometric information matrices, they are aligned with the geometric information database.

[0029] Step 4: align the geometric information and obtain the corresponding real geographic information;

[0030] Step 4.1: Use the geometric information matrix obtained in step 3 to align with the geometric information database in the first part. The alignment criteria are:

[0031] (1) Set the registration thresholds of three geometric information: angle difference matching threshold σ a =5, module length ratio matching threshold σ m =0.05, radius ratio matching threshold σ r =0.1;

[0032] (2) Take out the first element ξ1 of a certain geometric information matrix in the order of angle difference, modulus ratio and radius ratio, and traverse all geometric information matrices of the same type in the database. When there is an element ξ r satisfy:

[0033] ξ r -ξ1|≤σ w (7)

[0034] Then ξ1 is considered to match this number matrix and this number is recorded, where σ w The registration threshold representing this geometric information;

[0035] (3) Take the next element ξ2 and perform the same operation until the last element ξ is taken. n And complete this calculation process;

[0036] (4) Count the number of times this geometric information matrix is ​​successfully matched with each number matrix in the geometric information database, which should be a constant between 0 and n;

[0037] (5) The above process is used to perform the above operations on the other two geometric information matrices until all three geometric information matrices are calculated.

[0038] Step 4.2, using the results of step 4.1 to calculate the matching degree of each number matrix, the calculation criteria are:

[0039] (1) Calculate the matching weight of the angle difference, the modulus ratio and the radius ratio, denoted as δ A , δ M and δ R ;

[0040] (2) Count the number of successful matches of three different geometric information matrices with the same number in the geometric information database, recorded as N1, N2 and N3;

[0041] (3) Calculate the overall matching degree corresponding to this number

[0042]

[0043] Among them, N ad , N mr and N rr Respectively and The number of elements in ;

[0044] (4) Sort the overall matching degree of all numbers and obtain the number I with the highest matching degree h , then it is believed that The corresponding real crater number is I h Adjacent craters.

[0045] Step 4.3, finally, use the number I obtained in step 4.2 h The three-dimensional point coordinates of the real craters corresponding to the other two-dimensional centroid points in the captured image are obtained by using the adjacent crater groups and the existing digital maps.

[0046] Part 3: Calculation and filtering of aircraft posture;

[0047] After successfully obtaining the coordinates of the 2D centroid point in the captured image and the 3D point coordinates of the corresponding real crater, the coordinates of the aircraft itself can be solved using the relevant technology of camera pose estimation. The specific steps are as follows:

[0048] Step 5: Calculate the current aircraft posture;

[0049] The coordinate pair formed by the 2D centroid point calculated in the second part and the 3D point of the real crater is used as input data, and the current aircraft posture is calculated using the Efficient Perspective-n-Points (EPNP) algorithm. Assume u m =(μ m ,ν m )and For a set of 2D centroid points and 3D points, there are four control points And the corresponding weight value α m =[α m1 ,α m2 ,α m3 ,α m4 ] Τ satisfy:

[0050]

[0051] in, represents the three-dimensional point coordinates of the real crater in sequence m in the input data; j represents A sequence of a control point among the four control points; express The corresponding j-th control point, α mj express The corresponding weight value.

[0052] If point The corresponding point in the aircraft camera coordinate system is point The corresponding control points are Then it satisfies:

[0053]

[0054] in, express The corresponding three-dimensional coordinate point in the aircraft camera coordinate system; o represents A sequence of a control point among the four control points; express The corresponding o-th control point, α mo express The corresponding weight value.

[0055] Assume that the intrinsic parameter matrix of the aircraft camera is:

[0056]

[0057] Among them, f is the focal length of the camera, w is the width of the camera, and h is the height of the camera. According to the optical imaging law of the camera, we can get:

[0058]

[0059] Among them, d m Indicate point The depth of The Z value in u mRepresents the coordinates of the two-dimensional centroid point with a sequence of m in the input data; R represents the rotation matrix of the aircraft camera; t represents the translation vector of the aircraft camera.

[0060] Combining formula (9), formula (10), formula (11), and formula (12) and simplifying them, we get:

[0061]

[0062] in, They are represented as control points C o c The x, y, and z values ​​of

[0063] If n points are combined, we can get the system of equations:

[0064] M·x=0 (14)

[0065] Among them, M represents the parameter matrix; x represents The column vector formed by the coordinates of the four control points in .

[0066] By solving the equations in formula (14), the coordinates of the control points in the camera coordinate system can be obtained: Then solve the point After completing the above calculations, the rotation matrix and translation vector of the aircraft camera are solved based on the relationship between n pairs of three-dimensional points in the lunar world coordinate system and the three-dimensional points in the aircraft camera coordinate system. The results are corrected using the existing Gauss-Newton correction method in the field to finally obtain the positioning coordinates.

[0067] Step 6: Filtering of positioning results;

[0068] Step 6.1, since the application scenario of the present invention is that the lunar spacecraft performs daily navigation missions, it can be regarded as having collected the orbital operation data of previous navigations. Although the orbits of different navigations of the same daily navigation mission will change, the orbits of each navigation will not differ too much. Although the navigation orbits of the same daily navigation mission are different in different navigations, the relevant orbital operation data are still of reference significance to the orbit determination results of this navigation.

[0069] Define the relative angle moment T αis: a non-fixed navigation time value, when the spacecraft flies on different navigation orbits (the navigation starting point of the spacecraft is at the perigee of the navigation orbit, and only the case of sailing around the moon is considered), the navigation time required when the angle between the line connecting the spacecraft and the center point of the preset orbit of this navigation is α. Even in multiple different navigations, as long as the relative angle is the same, the angle between the line connecting the spacecraft and the center point of the preset orbit of this navigation is the same in each navigation.

[0070] Set a relative angle moment value We can get:

[0071]

[0072] Among them, A n Indicates the number of times this daily navigation mission has been executed; A k Indicates the number of sailings under this daily sailing mission; Indicates that the relative angle moment is When the aircraft and the A k The actual angle of the line connecting the center points of the preset orbit of the voyage; α r Indicates that the relative angle moment is When the aircraft and the A k The ideal angle of the line connecting the center points of the preset orbit of the voyage; Indicates that the relative angle moment is When the aircraft and the A k The actual distance between the center points of the preset orbit of the voyage; Indicates that the relative angle moment is When the aircraft and the A k The ideal distance of the center point of the preset orbit of the voyage; Indicates that the relative angle moment is When the aircraft is in A k The relative error of the angle between the line connecting the center point of the previous voyage and the preset orbit of this voyage; Indicates that the relative angle moment is When the aircraft is in A k The relative error between the distance from the center point of the preset orbit of this voyage and the previous voyage.

[0073] By using formula (15) and formula (16), the relative angle moment is calculated as The calculation method of the relative error of the angle between the aircraft and the center point of the preset navigation track during each past navigation and the relative error of the distance between the aircraft and the center point of the preset navigation track is:

[0074]

[0075] in, Indicates that the relative angle moment is The relative error mean of the angle between the aircraft and the center point of the preset orbit of this voyage during each previous voyage. Indicates that the relative angle moment is The relative error average of the distance between the aircraft and the center point of the preset orbit of this voyage during each previous voyage is .

[0076] Assume that the relative angle moment corresponding to the aircraft positioning coordinates solved in step 5 is When satisfied:

[0077]

[0078] in, Indicates the current relative angle moment (i.e. ) under the condition of the aircraft, the actual angle between the position of the aircraft and the center point of the preset orbit of this voyage; Indicates the current relative angle moment (i.e. ) under the following conditions: the actual distance between the position of the aircraft and the center point of the preset orbit of this voyage; Indicates the current relative angle moment (i.e. ) is the actual distance between the position of the aircraft and the center point of the preset orbit of this voyage. If the positioning result is within the acceptable error range, proceed to step 6.2; if it does not satisfy formula (19) and formula (20), the positioning result is considered unacceptable and steps 3 to 5 are repeated.

[0079] Step 6.2, use the current positioning moment data to predict the correct position that should be reached at the next moment. Since the trajectory will not change suddenly during the same flight, there should be no abrupt change between the positioning result of this time and the positioning result of the previous moment. Assume that in this voyage, the relative angle moment is The actual flight time corresponding to The predicted position of the aircraft at the next flight time is:

[0080]

[0081] Among them, x pre The x value represents the predicted position of the aircraft at the next flight time; y pre The y value represents the predicted position of the aircraft at the next flight time; z preIndicates the z value of the predicted position of the aircraft at the next flight time; a is the major semi-axis of the preset orbit of this flight; b is the minor semi-axis of the preset orbit of this flight; Indicates the actual sailing time corresponding to the next sailing time; T n Indicates the cycle of completing the navigation on the preset track of this navigation.

[0082] Step 6.3, all the positions passed by this voyage are fitted by using the unconstrained nonlinear least squares method related algorithm (such as fminsearch function, etc.) to obtain the fitted orbit, and compared with the preset orbit of the spacecraft this voyage set at the beginning, and three relative errors (relative error of the orbit semi-major axis, relative error of the orbit semi-minor axis, and relative error of the orbit focus) are obtained by calculating the orbit parameters. If the three relative errors are all lower than 0.5%, it can be considered that the fitted orbit has a high degree of fit with the preset orbit, and the orbit determination result is considered to be in a reasonable range and accepted, and the orbit determination is completed. If any of the three relative errors is higher than 0.5%, it is considered that the orbit determination has failed and the orbit determination result is not accepted.

[0083] The specific process of calculating the relative error is:

[0084] (1) Calculate the major semi-axis, minor semi-axis and focal length of the fitted orbit of this voyage, as well as the coordinates of the center point of the fitted orbit;

[0085] (2) Calculate the relative error of the semi-major axis (true value is a), the relative error of the semi-minor axis (true value is b), and the relative error of the focal length (true value is );

[0086] (3) Evaluate whether the three relative errors (relative error of the orbit semi-major axis, relative error of the orbit semi-minor axis, and relative error of the orbit focus) are all lower than 0.5%. If this condition is met, the orbit determination is considered successful; if not, the orbit determination is considered a failure.

[0087] The beneficial effects of the present invention are:

[0088] (1) The present invention uses existing scientific data and images to build a database and uses the captured images to perform real-time orbit determination, thereby reducing the reliance on ground measurement and control systems such as ground stations, improving the autonomy and real-time performance of orbit determination of the lunar orbiter when performing routine navigation tasks, and reducing resource consumption;

[0089] (2) At the same time, the present invention utilizes the characteristics of the camera's passive sensor in areas such as the far side of the moon where traditional orbit determination technology has a poor effect. It does not need to emit energy and only needs light to take pictures for orbit determination, thereby increasing the observation range.

[0090] (3) In addition, the present invention can fully utilize the characteristics of the lunar surface, such as the lack of atmospheric obstruction and the slow changes of the lunar surface, so that this vision-based orbit determination method has a higher reliability. BRIEF DESCRIPTION OF THE DRAWINGS

[0091] Figure 1 The present invention is a specific flow chart of the lunar surface photography and result registration part of the vision-based autonomous orbit determination method for a lunar circumnavigation vehicle. In the figure: 1 is the camera of the lunar circumnavigation vehicle, 2 is the lunar surface image obtained by taking pictures after attitude determination and focusing, 3 is the geometric information obtained by extracting the adjacent crater mass center point group closest to the image center from the image, from top to bottom are the angle difference, the modulus ratio and the radius ratio, 4 is three geometric information databases, 5 is the number corresponding to the database information matrix matched to each geometric information element, 6 is the matching summary result of each information matrix number in the database, 7 is the matching degree of each number under the three geometric information, 8 is the matrix number with the highest matching degree obtained by calculating the total matching degree, and 9 is the position information of the real adjacent crater group corresponding to the obtained number.

[0092] Figure 2 This is a specific flow chart of the calculation and filtering part of the aircraft posture of the vision-based lunar orbiting aircraft autonomous orbit determination method of the present invention. In the figure: 1 is the real position information corresponding to the group of adjacent crater mass center points closest to the image center obtained by matching, 2 is the crater position information corresponding to the remaining mass center points of the image obtained by calculation, 3 is the two-dimensional position information of the mass center point in the image, 4 is the aircraft camera posture solver, 5 is the initial solution result of the aircraft camera posture, 6 is the filter of the aircraft camera posture solution result, 7 is the final acceptable solution result, and 8 is the prediction result of the aircraft position at the next moment. DETAILED DESCRIPTION

[0093] The present invention is further described below in conjunction with specific embodiments.

[0094] Assuming that relevant data reserves are available, the number of each crater is recorded and the geographical location information of its center of mass is known, the vision-based autonomous orbit determination method for lunar spacecraft is divided into three parts: the establishment of a geometric information database, the lunar surface photography and result registration part, and the calculation and filtering of the spacecraft's posture.

[0095] Part 1: Geometric information database establishment part;

[0096] Step 1: Store the numbering information of adjacent crater groups;

[0097] The specific calculation method is as follows:

[0098] (1) Obtain the three-dimensional coordinates of the centroid of the real crater numbered 1, use the K-Nearest Neighbor (KNN) algorithm to solve the nine other centroids with the closest Euclidean distance to the centroid of this crater, and obtain the numbers of these nine craters. Store the number 1 and them into a number matrix. This number matrix is ​​considered to be the group of centroids of the adjacent craters corresponding to the three-dimensional coordinates of the centroid of the real crater numbered 1. At this time, there are a total of 10 crater numbers in this matrix, and the first number is 1;

[0099] (2) Obtain the three-dimensional coordinates of the centroid of the real crater numbered 2, and repeat the above operation until the three-dimensional coordinates of the centroid of the real crater numbered 2 are completed;

[0100] (3) The adjacent crater centroid point groups corresponding to each number are stored in order.

[0101] Step 2: Calculate and store the geometric information of each crater group;

[0102] Step 2.1, after obtaining the adjacent crater centroid point group corresponding to each numbered real crater centroid point according to step 1, calculate the angle difference matrix between each adjacent crater centroid point group, as follows:

[0103] (1) Obtain the coordinates of the centroid points of the two-dimensional image corresponding to each adjacent crater centroid point group from the reserve data. If the adjacent crater number matrix corresponding to the centroid point of the crater numbered J0 is [J0, J1, ..., J9], assume that the coordinates of the centroid point of the two-dimensional image corresponding to it are u Ji =(μ Ji ,ν Ji ),i=0,1,...,9, where u Ji Indicates the number is J i The coordinates of the centroid of the two-dimensional image corresponding to the centroid of the real crater, μ Ji represents the horizontal coordinate of the centroid of this two-dimensional image, ν Ji Indicates the vertical coordinate of the centroid of this two-dimensional image, i indicates the number J i The sequence of the adjacent crater number matrix corresponding to the crater center point numbered J0;

[0104] (2) According to formula (1), u Ji =(μ Ji ,ν Ji ), the origin of the two-dimensional plane coordinate system where i=0,1,...,9 is changed to u J0 At, the coordinates of the centroid points of the remaining two-dimensional point images (hereinafter referred to as two-dimensional point coordinates) are obtained;

[0105] (3) Combine the calculated two-dimensional point coordinates with formula (2) and formula (3) to calculate the angle difference matrix of the adjacent crater group: The matrix contains 36 elements in total;

[0106] Step 2.2, calculate the modulus ratio matrix and radius ratio matrix between the centroids of each adjacent crater group, as follows:

[0107] (1) Combine the two-dimensional point coordinates calculated in step 2.1 with formula (4) and formula (5) to calculate the modulus ratio matrix The matrix contains a total of 72 elements;

[0108] (2) Obtain the radius of each crater corresponding to the group of adjacent crater centroids in step 2.1 Combined with formula (6), the radius ratio matrix is ​​calculated The matrix contains a total of 90 elements.

[0109] The three geometric information matrices (angle difference matrix, modulus ratio matrix, radius ratio matrix) of each centroid are calculated and stored respectively, so as to establish a geometric information database. At this time, the number of each crater (or the real crater centroid) not only corresponds to itself, but also corresponds to a crater centroid point group including 10 crater centroids and three geometric information matrices.

[0110] Part II: lunar surface photography and result registration;

[0111] Step 3: Obtain real-time lunar images and solve related geometric information;

[0112] Assuming that the aircraft is traveling on a preset orbit, the visual axis direction and focal length of the aircraft camera are adjusted before photographing the lunar surface, as follows:

[0113] (1) Adjust the visual axis direction of the aircraft camera so that it faces the center of the moon, ensuring that the optical center of the camera on the aircraft is consistent with the projection point of the center of the moon on the two-dimensional pixel plane. At this time, the unit vector of the visual axis direction is the unit vector of the line vector connecting the aircraft camera and the center of the moon (the direction points to the center of the moon);

[0114] (2) Adjust the focal length of the aircraft camera so that the two-dimensional projection of the moon is tangent to the short side of the image. If the three-dimensional coordinates of the center of mass of the highest two-dimensional image and the center of mass of the lowest two-dimensional image in the image (referred to as the highest three-dimensional point coordinates and the lowest three-dimensional point coordinates) are and in, represents the coordinates of the highest three-dimensional point, Represents the x value, y value, and z value of the highest three-dimensional point coordinates, P lw represents the coordinates of the lowest three-dimensional point, Represents the x, y, and z values ​​of the lowest three-dimensional point coordinates. If the external parameter matrix of the aircraft camera [R m |t m ]for:

[0115]

[0116] Among them, R m is the rotation matrix of the aircraft camera at this time; t m is the translation vector of the aircraft camera at this time. Then the camera focal length f0 should satisfy:

[0117]

[0118] Where h0 is the height of the aircraft camera; ξ h is the high offset; l is a low offset. h With ξ l The calculation method is:

[0119]

[0120] Among them, all variables come from formula (22) and are constants;

[0121] (3) After adjusting the visual axis direction and focal length of the aircraft camera, take a photo of the lunar surface and extract the centroid of the crater. Obtain the coordinates of the centroid of the crater closest to the image center (i.e., the optical center) and mark them as Among them, Indicates the number is G l The coordinates of the centroid of the two-dimensional image corresponding to the centroid of the real crater, represents the horizontal coordinate of the centroid of this two-dimensional image, Represents the vertical coordinate of the centroid of this two-dimensional image, l represents the number G l The sequence of the adjacent crater number matrix corresponding to the crater center point numbered G0. The angle difference matrix is ​​obtained by formula (1), formula (2), and formula (3): The modulus length ratio matrix is ​​obtained by formula (1), formula (4), and formula (5): From formula (1) and formula (6), we can get the radius ratio matrix: As attached Figure 1 As shown in 3, after obtaining the three geometric information matrices, they are aligned with the geometric information database.

[0122] Step 4: align the geometric information and obtain the corresponding real geographic information;

[0123] Step 4.1: Using the data obtained in step 3 and Compare with the geometric information database in the first part, and count and The matrix number of each element in the database that is successfully matched, as shown in the attached Figure 1 To reduce the amount of calculation, when an element in a geometric information matrix in the geometric information database already meets the matching condition, the other elements in the geometric information matrix will no longer be matched. After the statistics are completed, the number of successful matches for each number in the geometric information database is calculated, as shown in the attached figure. Figure 1 As shown in 6, the number of times will not exceed the number of elements of the geometric information matrix obtained in step 2. For example, the number of times a certain numbered angle difference matrix in the geometric information database is successfully matched will not exceed 36.

[0124] Step 4.2, using the result of step 4.1 to calculate the matching degree of each number matrix. The matching weight δ of this angle difference, modulus ratio and radius ratio A , δ M and δ R for:

[0125]

[0126] The matching weight is used to calculate the three geometric information matching degrees (angle difference matching degree, modulus length ratio matching degree, radius ratio matching degree) corresponding to each number in the geometric information database, as shown in the attached figure. Figure 1 As shown in 7.

[0127] Step 4.3, sort the overall matching degree of all numbers and obtain the number I with the highest matching degree h , as attached Figure 1 As shown in 8. The corresponding real crater is numbered I h The group of centroids of adjacent craters is numbered I h The real position information of the centroid points of the adjacent craters is shown in the attached Figure 2 Then, using this adjacent crater centroid point group and the existing digital map, the three-dimensional coordinates of the real crater centroid points corresponding to the remaining two-dimensional centroid points in the captured image are obtained, as shown in the attached figure. Figure 2 As shown in 2.

[0128] Part 3: Calculation and filtering of aircraft posture;

[0129] Step 5: Calculate the current aircraft posture;

[0130] After successfully obtaining the coordinates of the two-dimensional centroid point in the captured image and its corresponding three-dimensional centroid point coordinates, the coordinates of the aircraft itself are solved using the relevant technology of camera pose estimation. The specific steps are as follows:

[0131] (1) The acquired 2D centroid coordinates and the real 3D point coordinate pair are used as initial data and input into the aircraft camera pose solver, as shown in the attached figure. Figure 2 As shown in 4. The four control points and weight values ​​corresponding to each real three-dimensional coordinate point in the lunar world coordinate system are solved by the aircraft camera pose solver. After the solution is completed, the unknown parameter point depth d is eliminated using formula (9), formula (10), formula (11), and formula (12). m , and obtain the coordinate equation in the corresponding aircraft camera coordinate system. Multiple sets of coordinate equations are combined to obtain a linear equation system in the form of formula (14), where M is the parameter matrix. Solving the right null space x of M obtains the coordinates of the control points in the aircraft camera coordinate system, and then the three-dimensional coordinates in the corresponding aircraft camera coordinate system are obtained by linearly combining the coordinates of the four control points and the weight values.

[0132] (2) Based on multiple sets of three-dimensional points in the lunar world coordinate system 3D point with aircraft camera coordinate system Relationship:

[0133] P i w =R·P i c +t (26)

[0134] Among them, P i w represents the coordinates of a three-dimensional point in the lunar world coordinate system; R represents the rotation matrix of the aircraft camera; P i c Indicates P i w The corresponding three-dimensional point coordinates in the aircraft camera coordinate system; t is the translation vector of the aircraft camera.

[0135] The iterative closest point (ICP) method is used to solve the aircraft camera rotation matrix R and translation vector t, and the results are corrected by the existing Gauss-Newton method in the field to obtain the preliminary positioning coordinates Τ p =(x p ,y p ,z p ), where T p Represents the initial positioning coordinates, (x p ,y p ,zp ) indicates T p The x value, y value and z value of Figure 2 As shown in 5.

[0136] Step 6: Filtering of positioning results;

[0137] Step 6.1, input the preliminary positioning coordinates into the filter, as shown in the attached Figure 2 The filter stores the relevant orbit data and the positioning status of the aircraft at each moment during this voyage, and uses the past data and motion status to filter the positioning result at this moment. Assume that the relative angle moment corresponding to the positioning result is The coordinates of the preset orbit center of this voyage (not the center of the moon, the center of the moon is the preset orbit focus) in the lunar world coordinate system are The coordinates of the preset orbit perigee (i.e. the preset orbit vertex) in the lunar world coordinate system are in, Indicates the coordinates of the preset orbit center of this voyage in the lunar world coordinate system. express The x-value, y-value, and z-value of Indicates the coordinates of the perigee of the preset orbit of this voyage in the lunar world coordinate system. express The x-value, y-value, and z-value of

[0138] At this time (the relative angle is ), the distance between the aircraft and the center of the preset track and the angle between the aircraft and the center of the preset track are:

[0139]

[0140] Among them, d1 represents the distance between the spacecraft and the center of the preset orbit at this time; v1 represents the vector formed by the center of the preset orbit and the perigee in the lunar world coordinate system (the direction points to the perigee of the preset orbit); v2 represents the vector formed by the center point of the preset orbit and the preliminary positioning coordinates in the lunar world coordinate system (the direction points to the preliminary positioning coordinates); θ1 represents the angle between the spacecraft and the center of the preset orbit.

[0141] When the relative angle is The ideal angle between the aircraft and the center point of the preset orbit corresponding to this voyage is α p , the corresponding ideal distance between the aircraft and the center point of the preset orbit is:

[0142]

[0143] in, It represents the ideal distance between the aircraft and the center point of the preset orbit corresponding to this voyage; a1 represents the major semi-axis of the preset orbit of this voyage; b1 represents the minor semi-axis of the preset orbit of this voyage.

[0144] If d1 and θ1 satisfy the conditions:

[0145]

[0146] in, Indicates that the relative angle moment is The relative error mean of the distance between the aircraft and the center point of the preset orbit of this voyage during each previous voyage; Indicates that the relative angle moment is When d1 and θ1 do not meet this condition, the positioning result is considered unacceptable and steps 3 to 5 are repeated.

[0147] Step 6.2, use the current positioning moment data to predict the correct position that should be reached at the next moment. Since the navigation track will not change suddenly during the same voyage, there should be no abrupt change between the positioning result of this time and the positioning result of the previous moment. Assume that in this voyage, the relative angle moment is The actual flight time corresponding to (i.e. the actual flight time corresponding to the initial positioning coordinates of this time), if the actual flight time at the next moment is Then the predicted position of the aircraft at the next flight time can be calculated by formula (21), as shown in the following figure: Figure 2 As shown in 8.

[0148] Step 6.3, use the fminsearch function to fit the trajectory of all the positions passed by the aircraft in this voyage to obtain the fitted orbit, and compare the fitting result with the preset orbit of the aircraft set at the beginning, calculate the orbital major semi-axis, orbital minor semi-axis and orbital focal length of the fitted orbit, and calculate the relative error of the major semi-axis, the relative error of the minor semi-axis and the relative error of the orbital focal length between the fitted orbit and the preset orbit (the true values ​​are the orbital major semi-axis, orbital minor semi-axis and orbital focal length of the preset orbit respectively), and make a final judgment on the orbit determination result of this time, as shown in the attached figure. Figure 2As shown in 7, when the three relative errors (relative error of the major axis, relative error of the minor axis and relative error of the orbit focal length) are all less than 0.5%, it is considered that the fitted orbit has a high degree of fit with the preset orbit, and the orbit determination result is considered to be in a reasonable range and accepted, completing the orbit determination. If any of the three relative errors is higher than 0.5%, the orbit determination is considered to have failed and the orbit determination result is not accepted.

[0149] The above-described embodiments merely express the implementation methods of the present invention, but they cannot be understood as limiting the scope of the patent of the present invention. It should be pointed out that for those skilled in the art, several modifications and improvements can be made without departing from the concept of the present invention, which all belong to the protection scope of the present invention.

Claims

1. A vision-based autonomous orbit determination method for a lunar orbiting vehicle, characterized in that: The autonomous orbit determination method for a lunar circumnavigation vehicle first extracts the centroid of lunar craters based on lunar surface images and data and constructs a geometric information database for registration; Secondly, set the preset orbit of the spacecraft, adjust the direction of the spacecraft's visual axis and the focal length of the spacecraft's camera, take pictures of the moon, and find the corresponding real crater information by aligning the captured images with the geometric information database, and use the camera pose estimation method to calculate the current position of the spacecraft; thirdly, filter the solution results to complete the orbit determination autonomously and in real time.

2. The method for autonomous orbit determination of a lunar orbiter based on vision according to claim 1, characterized in that: The autonomous orbit determination method for a lunar circumnavigation vehicle assumes that a lunar world coordinate system with the center of the moon as the origin is established, the real geographical location coordinates of the main craters are recorded through existing lunar surface images, different craters are numbered and the centroid points of the craters are extracted; at the same time, a preset orbit of the spacecraft is set, wherein the preset orbit is an elliptical orbit and is located in the XOY plane of the lunar world coordinate system, the center of the moon is located at the left focus of the preset orbit, and the initial position of the spacecraft is the perigee of the preset orbit, that is, the left vertex of the ellipse; the vision-based autonomous orbit determination method for a lunar circumnavigation vehicle is divided into three parts: a geometric information database establishment part, a lunar surface photography and result registration part, and a spacecraft posture solution and filtering part; The method for autonomous orbit determination of a lunar orbiting vehicle comprises the following steps: Part 1: Geometric information database establishment part; Step 1: Store the numbering information of adjacent crater groups; Calculate the N nearest craters to each crater except itself according to the real geographical coordinates, record and store the numbers of these adjacent crater groups, then each crater corresponds to a number matrix containing N+1 elements in total, and the first element of this number matrix is ​​its own number; Step 2: Calculate and store the geometric information of each crater group; Step 2.1, first, obtain the coordinates of the centroid of the two-dimensional image corresponding to each adjacent crater group; if the adjacent crater number matrix corresponding to the crater numbered k0 is [k0, k1, ..., k N ], assuming that the corresponding coordinates of the centroid of the two-dimensional image are in Indicates number k n The coordinates of the centroid of the two-dimensional image corresponding to the crater, represents the horizontal coordinate of the centroid of this two-dimensional image, Represents the vertical coordinate of the centroid of this two-dimensional image; Change the origin of the 2D plane coordinate system to At, the coordinates of the remaining two-dimensional image centroid points are adjusted to: in, Indicates number k n The coordinates of the centroid point of the two-dimensional image corresponding to the crater in the adjusted two-dimensional plane coordinate system (hereinafter referred to as the two-dimensional point); Step 2.2, calculate the angle difference matrix of this adjacent crater group Modulus ratio matrix Ratio matrix to radius Calculate and store three geometric information matrices of each centroid point respectively to establish a geometric information database, wherein the three geometric information matrices include an angle difference matrix, a modulus ratio matrix, and a radius ratio matrix; Part II: lunar surface photography and result registration; Step 3: Obtain real-time lunar images and solve related geometric information; Assuming that the spacecraft is traveling on a preset orbit, the optical center of the camera on the spacecraft is aligned with the projection point of the moon's spherical center on the two-dimensional pixel plane, and the focal length is adjusted so that the moon's projection is tangent to the short side of the image taken by the camera; Take a picture with the camera and extract the centroid of the crater. Get the coordinates of the centroid of the crater closest to the center of the image and mark them as in represents the coordinates of the 2D point with sequence q in the center point group of the crater closest to the center of the image; Represents a two-dimensional point The horizontal axis of Represents a two-dimensional point The vertical coordinate of Sequence in the cluster of crater centroids closest to the image center; Obtain three geometric information matrices according to step 2; Step 4: align the geometric information and obtain the corresponding real geographic information; Step 4.1, using the geometric information matrix obtained in step 3 to align with the geometric information database in the first part; Step 4.2, using the results of step 4.1 to calculate the matching degree of each number matrix, the calculation criteria are: (1) Calculate the matching weight of the angle difference, the modulus ratio and the radius ratio, denoted as δ A , δ M and δ R ; (2) Count the number of successful matches of three different geometric information matrices with the same number in the geometric information database, recorded as N1, N2 and N3; (3) Calculate the overall matching degree corresponding to this number Among them, N ad , N mr and N rr Respectively and The number of elements in ; (4) Sort the overall matching degree of all numbers and obtain the number I with the highest matching degree h , then it is believed that The corresponding real crater number is I h The adjacent crater group; Step 4.3, using the number obtained in step 4.2 as I h The adjacent crater groups and the digital map obtain the three-dimensional point coordinates of the real crater corresponding to the other two-dimensional centroid points in the captured image; Part 3: Calculation and filtering of aircraft posture; The coordinates of the aircraft itself are solved using the relevant technology of camera pose estimation. The specific steps are as follows: Step 5: Calculate the current aircraft posture; The coordinate pair formed by the two-dimensional centroid point calculated in the second part and the three-dimensional point of the real crater is used as input data, and the current aircraft posture is calculated using an efficient n-point perspective algorithm; Step 6: Filter the positioning results; Step 6.1, define the relative angle moment T α is: a non-fixed flight time value, which is the flight time required when the line connecting the aircraft and the center point of the preset orbit of this flight passes through an angle of α when the aircraft flies on different navigation orbits; even in multiple different navigations, as long as the relative angle moment is the same, the angle passed by the line connecting the aircraft and the center point of the preset orbit of this flight is the same in each navigation; Set a relative angle time value T αr , we can get: Among them, A n Indicates the number of times this daily navigation mission has been executed; A k Indicates the number of sailings under this daily sailing mission; Indicates that the relative angle moment is When the aircraft and the A k The actual angle of the line connecting the center points of the preset orbit of the voyage; α r Indicates that the relative angle moment is When the aircraft and the A k The ideal angle of the line connecting the center points of the preset orbit of the voyage; Indicates that the relative angle moment is When the aircraft and the A k The actual distance between the center points of the preset orbit of the voyage; Indicates that the relative angle moment is When the aircraft and the A k The ideal distance of the center point of the preset orbit of the voyage; Indicates that the relative angle moment is When the aircraft is in A k The relative error of the angle between the line connecting the center point of the previous voyage and the preset orbit of this voyage; Indicates that the relative angle moment is When the aircraft is in A k The relative error of the distance between the center point of the preset orbit of this voyage and the previous voyage; By using formula (15) and formula (16), the relative angle moment is calculated as The calculation method of the relative error of the angle between the aircraft and the center point of the preset navigation track during each past navigation and the relative error of the distance between the aircraft and the center point of the preset navigation track is: in, Indicates that the relative angle moment is The relative error mean of the angle between the aircraft and the center point of the preset orbit of this voyage during each previous voyage. Indicates that the relative angle moment is The relative error mean of the distance between the aircraft and the center point of the preset orbit of this voyage during each previous voyage; Assume that the relative angle moment corresponding to the aircraft positioning coordinates solved in step 5 is When satisfied: in, Indicates the current relative angle moment during this voyage The actual angle between the position of the aircraft and the center point of the preset orbit of this voyage; Indicates the current relative angle moment during this voyage The actual distance between the position of the aircraft and the center point of the preset orbit of this voyage; Indicates the current relative angle moment during this voyage If the actual distance between the position of the aircraft and the center point of the preset orbit of this voyage is 1, the positioning result is considered to be within the acceptable error range, and step 6.2 is performed; if formula (19) and formula (20) are not satisfied, the positioning result is considered to be unacceptable, and the operations of steps 3 to 5 are repeated; Step 6.2, use the current positioning time data to predict the correct position that should be reached at the next moment; assuming that in this voyage, the relative angle moment is The actual flight time corresponding to The predicted position of the aircraft at the next flight time is: Among them, x pre The x value represents the predicted position of the aircraft at the next flight time; y pre The y value represents the predicted position of the aircraft at the next flight time; z pre Indicates the z value of the predicted position of the aircraft at the next flight time; a is the major semi-axis of the preset orbit of this flight; b is the minor semi-axis of the preset orbit of this flight; Indicates the actual sailing time corresponding to the next sailing time; T n Indicates the cycle of completing the voyage on the preset track of this voyage; Step 6.3, perform trajectory fitting on all positions passed by this voyage to obtain a fitted orbit, and compare it with the preset orbit of the spacecraft for this voyage set at the beginning, and obtain three relative errors by calculating the orbit parameters; evaluate whether the three relative errors meet the conditions, if they meet the conditions, the orbit determination is considered to be successful, if not, the orbit determination is considered to have failed.

3. The vision-based autonomous orbit determination method for a lunar orbiting vehicle according to claim 2, characterized in that: The step 2.2 is as follows: First calculate the angle difference matrix, assuming but: in, Represented as a two-dimensional point The horizontal axis of Represented as a two-dimensional point The vertical coordinate of Represents a two-dimensional point The angle between the vector formed with the origin and the horizontal axis of the coordinate system; Represents a two-dimensional point The angle between the vector formed by the origin and the horizontal axis of the coordinate system; Δθ mn express and The angle difference between them; m means number k m The sequence of the crater numbered in the adjacent crater number matrix corresponding to the crater numbered k0; n represents the crater numbered k n The sequence of the crater numbered in the adjacent crater number matrix corresponding to the crater numbered k0; The angle difference matrix of point k0 is calculated as follows: Assume the number is k n The radius of the crater is Then the module length ratio matrix Ratio matrix to radius The calculation method is as follows: in, Represents a two-dimensional point The distance from the origin; Represents a two-dimensional point The distance from the origin; Represents a two-dimensional point and 2D points The ratio of the distance to the origin; Indicates number k m The radius of the crater; Indicates number k n The radius of the crater; Indicates that the number is k n The crater numbered k m The ratio of the radius of the crater to that of the crater; and then we get 4. The vision-based autonomous orbit determination method for a lunar orbiting vehicle according to claim 3, characterized in that: In step 3, the method for obtaining three geometric information matrices is: The two-dimensional point coordinates after the origin is changed are generated by formula (1) The angle difference matrix is ​​obtained by formula (2) and formula (3): The modulus ratio matrix is ​​obtained from formula (4) and formula (5): The radius ratio matrix is ​​obtained from formula (6):

5. The vision-based autonomous orbit determination method for a lunar orbiting vehicle according to claim 2, characterized in that: In step 4.1, the registration criteria are: (1) Set the registration thresholds of three geometric information: angle difference matching threshold σ a =5, module length ratio matching threshold σ m =0.05, radius ratio matching threshold σ r =0.1; (2) Take out the first element ξ1 of a certain geometric information matrix in the order of angle difference, modulus ratio and radius ratio, and traverse all geometric information matrices of the same type in the database. When there is an element ξ r satisfy: x r -ξ1|≤σ w (7) Then ξ1 is considered to match this number matrix and this number is recorded, where σ w The registration threshold representing this geometric information; (3) Take the next element ξ2 and perform the same operation until the last element ξ is taken. n And complete this calculation process; (4) Count the number of times this geometric information matrix is ​​successfully matched with each number matrix in the geometric information database, which should be a constant between 0 and n; (5) The above process is used to perform the above operations on the other two geometric information matrices until all three geometric information matrices are calculated.

6. The vision-based autonomous orbit determination method for a lunar orbiting vehicle according to claim 2, characterized in that: The method for calculating the current aircraft posture in step 5 is: Assumptions For a set of 2D centroid points and 3D points, there are four control points And the corresponding weight value α m =[α m1 ,α m2 ,α m3 ,α m4 ] Τ satisfy: in, represents the three-dimensional point coordinates of the real crater in sequence m in the input data; j represents A sequence of a control point among the four control points; express The corresponding j-th control point, α mj express The corresponding weight value; If point The corresponding point in the aircraft camera coordinate system is point The corresponding control points are Then it satisfies: in, express The corresponding three-dimensional coordinate point in the aircraft camera coordinate system; o represents A sequence of a control point among the four control points; express The corresponding o-th control point, α mo express The corresponding weight value; Assume that the intrinsic parameter matrix of the aircraft camera is: Among them, f is the focal length of the camera, w is the width of the camera, and h is the height of the camera. According to the optical imaging law of the camera, we can get: Among them, d m Indicate point The depth of The Z value in u m represents the coordinates of the two-dimensional centroid point with a sequence of m in the input data; R represents the rotation matrix of the aircraft camera; t represents the translation vector of the aircraft camera; Combining formula (9), formula (10), formula (11), and formula (12) and simplifying them, we get: in, Represented as control points The x, y and z values ​​of If n points are combined, we get the system of equations: M·x=0 (14) Among them, M represents the parameter matrix; x represents The column vector formed by the coordinates of the four control points in ; By solving the equations in formula (14), we can obtain the coordinates of the control points in the camera coordinate system: Then solve the point After completing the above calculations, the rotation matrix and translation vector of the aircraft camera are solved based on the relationship between n pairs of three-dimensional points in the lunar world coordinate system and the three-dimensional points in the aircraft camera coordinate system. The results are corrected using the Gauss-Newton correction method to obtain the positioning coordinates.

7. The vision-based autonomous orbit determination method for a lunar orbiting vehicle according to claim 2, characterized in that: In step 6.3, all positions passed by this voyage are fitted with trajectories using an unconstrained nonlinear least squares correlation algorithm to obtain a fitted orbit.

8. The vision-based autonomous orbit determination method for a lunar orbiting vehicle according to claim 2, characterized in that: In step 6.3, the three relative errors include the relative error of the orbit semi-major axis, the relative error of the orbit semi-minor axis, and the relative error of the orbit focus; if the three relative errors are all lower than 0.5%, it is considered that the fitted orbit has a high degree of fit with the preset orbit, and the orbit determination result is considered to be within a reasonable range and accepted to complete the orbit determination; if any of the three relative errors is higher than 0.5%, it is considered that the orbit determination has failed and the orbit determination result is not accepted.

9. The vision-based autonomous orbit determination method for a lunar orbiting vehicle according to claim 2, characterized in that: In step 6.3, the specific process of calculating the relative error is: (1) Calculate the major semi-axis, minor semi-axis and focal length of the fitted orbit of this voyage, as well as the coordinates of the center point of the fitted orbit; (2) Calculate the relative error of the semi-major axis, the relative error of the semi-minor axis, and the relative error of the orbit focal length between the fitted orbit of this voyage and the preset orbit respectively; define the relative error of the semi-major axis as the true value a, the semi-minor axis as the true value b, and the relative error of the orbit focal length as the true value

Images