Method for constructing photometric high-resolution terrain navigation landmark only constrained by single image

Abstract

The invention relates to a photometric high-resolution terrain navigation landmark construction method only constrained by a single image. The method comprises the following steps: acquiring a single high-resolution image and low-resolution initial terrain data of a target asteroid surface; under a Bayesian framework, a terrain and albedo inversion problem is modeled into a posterior probability maximization problem, a Bayesian model is constructed, a likelihood item in the model is used for describing a relationship among luminosity information of an image, the terrain and the albedo, an inversion result is ensured to be consistent with an observation image, and the accuracy of the inversion result is improved. A local terrain smoothness constraint and a local albedo consistency hypothesis are introduced into a priori item so as to solve the ill-conditioned problem under the constraint of a single image; converting the Bayesian model into an optimization problem under a variational framework, and iteratively optimizing the terrain and albedo by adopting an alternating direction multiplier method; and constructing a terrain navigation landmark based on the high-resolution terrain and albedo data obtained by inversion. Compared with the prior art, the method has the advantages of being capable of achieving fine three-dimensional terrain and albedo inversion and the like.

Description

Technical Field

[0001] The present invention relates to the technical field of navigation landmark construction, and more particularly to a photometric high-resolution terrain navigation landmark construction method with only single-image constraint. Background Art

[0002] The deep space environment is complex, the long-distance observation conditions are limited, and the asteroid exploration tasks are complex and refined, all of which urgently rely on higher-precision visual positioning technologies. The refined asteroid terrain navigation landmarks are the primary prerequisite and basic information guarantee for supporting the asteroid to achieve attachment and sampling return detection visual navigation. In view of the major demand for the construction of fine terrain navigation landmarks on the asteroid surface, a very feasible and reliable method is visual navigation based on terrain navigation landmarks.

[0003] The visual navigation method based on terrain navigation landmarks was first applied in the Hayabusa-1 mission. This visual navigation method is called OBIRON (On-Board Image Registration for Optical Navigation). Different from the matching of navigation landmarks and navigation maps based on graphic features, the OBIRON method does not rely on a complete navigation map of the landing area. It comes with a set of terrain-based navigation landmarks. Each landmark includes a small piece of surface elevation grid and an albedo map. OBIRON does not directly match the landing image with the navigation map, but matches the image block features predicted based on the navigation landmarks with the navigation camera image. Among them, the predicted image block is rendered according to the prior known detector state and navigation camera model. Subsequently, the correlation matching method is used to predict the landmark positions of the image and the observed image, and the difference is calculated between the predicted position and the actual position of each landmark image, so as to obtain the least squares solution of the spacecraft pose. By iteratively solving and approaching the asteroid, the navigation accuracy and landmark position accuracy are continuously improved. The SPC technology is a stereophotometric modeling method specifically proposed by Gaskell et al. for the construction of navigation landmarks. SPC plays a crucial role in asteroid exploration missions and has been successfully applied to the three-dimensional reconstruction of many asteroid exploration missions such as Dawn, Rosetta, Hayabusa2, and Osiris. It can process multi-resolution images and obtain an accuracy comparable to that of the highest-resolution images. The images used in this method require multiple images to simultaneously meet specific observation combination conditions and lighting conditions.

[0004] Considering the resolution constraints between the initial terrain and the corresponding regional images, as well as the constraints of lighting and observation conditions, in many scenarios, there is no image set but only a single image that can be used to optimize the terrain. For example, taking the landing mission of the OSIRIS-Rex mission as an example, before the official landing, it arranged two rehearsal previews, that is, the detector flew along the preset landing trajectory for rehearsal. This process will obtain images with higher resolution and observation and lighting conditions very similar to those during the actual landing. Obviously, the resolution, lighting, and observation conditions of these images are unique. However, there is currently a lack of a method for constructing photometric navigation landmarks based on a single image constraint. Because compared with the terrain construction method based on multiple image constraints, the photometric construction of terrain and albedo based on a single image constraint is essentially an ill-posed problem. Summary of the Invention

[0005] The purpose of the present invention is to provide a method for constructing a photometric high-resolution terrain navigation landmark with only a single image constraint, using the assumption of "constant albedo in the local area" to achieve the construction of a high-resolution terrain navigation landmark under the mutual constraint of a single high-resolution image and a low-resolution terrain. In the constructed Bayesian framework, fine three-dimensional terrain and albedo inversion are realized.

[0006] The purpose of the present invention can be achieved through the following technical solutions:

[0007] A method for constructing a photometric high-resolution terrain navigation landmark with only a single image constraint, comprising the following steps:

[0008] Data acquisition and preprocessing: Obtain a single high-resolution image of the surface of the target asteroid and low-resolution initial terrain data;

[0009] Establish a terrain and albedo inversion model under the Bayesian model: Under the Bayesian framework, model the inversion problem of terrain and albedo as a problem of maximizing the posterior probability, and construct a Bayesian model considering the prior of local albedo consistency. The Bayesian model includes a likelihood term and a prior term. Among them, the likelihood term is used to describe the relationship between the photometric information of the image and the terrain and albedo, ensuring that the inversion result is consistent with the observed image. The prior term introduces local terrain smoothness constraints and the assumption of local albedo consistency to solve the ill-posed problem under a single image constraint;

[0010] Terrain and albedo inversion: Convert the Bayesian model into an optimization problem under the variational framework, and use the alternating direction method of multipliers to iteratively optimize the terrain and albedo, gradually approaching the optimal solution. In each iteration, update the terrain and albedo respectively, and introduce Lagrange multipliers to handle the constraint conditions. The constraint conditions include local albedo consistency constraints and terrain smoothness constraints;

[0011] Navigation landmark construction: Based on the high-resolution terrain and albedo data obtained by inversion, navigation landmarks are constructed. The navigation landmarks are defined by landmark positions, reference coordinate systems, grid sizes, and relative elevation and albedo maps. Each navigation landmark represents a local area on the surface of the target celestial body.

[0012] The specific construction process of the Bayesian model considering the prior of local albedo consistency is as follows:

[0013] Use photometric method to recover high-resolution terrain from low-resolution initial terrain. Based on the photometric constraints of a single high-resolution image, invert the high-resolution terrain from the low-resolution initial terrain. Among them, the high-frequency information required in the high-resolution terrain inversion corresponds to the terrain details constrained by the high-resolution image, and the low-frequency information corresponds to the geometric structure constrained by the initial terrain;

[0014] According to Bayes' rule, recover the high-resolution terrain by maximizing the posterior distribution of a single high-resolution image and the low-resolution terrain, and construct a Bayesian model based on the posterior distribution.

[0015] Given a low-resolution terrain, denote its elevation as z0, and assume a high-resolution terrain z, and The relationship between the low-resolution terrain and the high-resolution terrain is expressed as:

[0016] z0 = Dz + ε z

[0017] where ε z represents a random process; D: Ω HR →Ω LR represents a linear downsampling operation operator;

[0018] Use high-resolution images to constrain the process of generating high-resolution from low-resolution terrain. Based on photometric method, according to the SFS principle, establish the relationship between the image and the albedo model and minimize it to obtain the normal vector of any point in the image:

[0019]

[0020] where the normal vector is a function of the slope of the current point, that is, n(x, y) = g(z x (x, y), z y (x, y)), g(·) is a function expression, z x (x, y) and z y(x, y) represents the gradient of z(x, y), where z(x, y) is the elevation of pixel (x, y), I(x, y) represents the image gray value of pixel (x, y), and R(x, y) represents the reflectivity corresponding to pixel (x, y). The reflectivity is calculated based on the lunar-Lambert reflectivity model using the observation and illumination conditions at the time of the current image capture; the slope characterizes the magnitude of the elevation gradient.

[0021] The posterior distribution is expressed as:

[0022]

[0023] where, I represents the image, a represents the albedo, and P represents the probability density function.

[0024] The construction process of the likelihood term includes the following steps:

[0025] Assume that the image I and the low-resolution terrain z0 follow independent distributions, then:

[0026] P(z0, I|z, a) = P(z0|z, a)P(I|z, a)

[0027] Assume that the low-resolution terrain z0 is independent of the surface albedo a, then:

[0028] P(z0, I|z, a) = P(z0|z)P(I|z, a)

[0029] Assume the noise that is, it follows a Gaussian distribution with a mean of 0 and a variance of σ z then:

[0030]

[0031] where, represents the two-norm operator, and D(·) represents the linear downsampling operation;

[0032] Construct the relationship between the image I and the reflectivity R: Considering the self-shadow characteristic of the terrain, express the reflectivity R as a function of the normal vector corresponding to the current terrain, then:

[0033] I(x, y) = Λ × a(x, y) × max{R(n(x, y), s(x, y), v(x, y)), 0} + ε I

[0034] Among them, (x, y) represents the pixel index in the image. The max{·} operation constrains the self-shadow characteristics, that is, when the object is occluded or should not receive light, the self-shadow factor is zero; s(x, y) and v(x, y) are the unit vectors of the incident light and the outgoing light at (x, y) respectively, n(x, y) is the normal vector at (x, y), Λ represents the parameter related to the image exposure. There are many factors affecting Λ, and it is usually regarded as a constant during the terrain optimization process. a(x, y) represents the albedo at (x, y), and ε I represents the image noise, which is obtained through pre-calibration;

[0035] Before estimating the terrain and albedo, Λ is initialized as a constant, that is:

[0036]

[0037] Among them, x and y respectively define the integration intervals in the row and column directions of the image;

[0038] Assume That is, it follows a Gaussian distribution with a mean of 0 and a variance of σ I , then there is:

[0039]

[0040] When the incident light vector, the outgoing light vector, and the camera internal parameters are known, whether it is orthographic projection or perspective projection, n(x, y) is essentially a function of, then:

[0041]

[0042] The construction process of the prior term includes the following steps:

[0043] Assume that the terrain and reflectivity are independent of each other, then there is:

[0044] P(z, a) = P(z)P(a)

[0045] For the constraint of the terrain z, a second-order regularization is used to constrain the smoothness of the local terrain, and the terrain optimization is realized by combining the edge protection characteristics of the second-order coupling model and the directional image information:

[0046]

[0047] Among them, G is the Laplacian kernel function, and w z represents the bandwidth of the kernel function, which is used to control the size of the smoothness. Δz represents the second derivative of the terrain z, and μ is used to control the dispersion degree of the distribution, which is set as a constant;

[0048] For the constraint on the albedo \(a\), the local region albedo consistency constraint is adopted, that is, it is assumed that each superpixel has a consistent albedo characteristic, and it is realized by using the superpixel segmentation method. For the superpixel segmentation result, it is represented by an indicator function, that is:

[0049]

[0050] If the above process is considered to follow the Pott model, then there is:

[0051]

[0052] where the parameter \(\lambda>0\), \(\Omega\) i \(\in\Omega\) HR , \(i\) is the superpixel index, \(n\) is the total number of superpixels, represents the indicator function, \(|\cdot|\) represents the distance defined by the superpixel, and \(\lambda\) is used to control the discreteness of this distribution and is set as a constant.

[0053] The optimal estimation of the terrain and albedo is achieved by maximizing the posterior distribution, that is, by minimizing its negative logarithm. Then, the following energy function is constructed:

[0054]

[0055] where \(D\) I (z,a) expresses the assumption based on photometry, that is, the relationship between pixel intensity, albedo and reflectance model; represents the constraint of the initial terrain; \(L\) z (z) and \(L\) a (a) are the regularization terms for the terrain and albedo respectively; \(\beta\), \(\gamma\) z and \(\gamma\) a correspond to the weights of each data term and regularization term respectively;

[0056]

[0057]

[0058] Express the energy function in the form of an optimization problem under the variational framework, and obtain:

[0059]

[0060] By introducing an auxiliary vector field to directly solve for \(z\), the above optimization problem under the variational framework is re-expressed as an equivalent constrained optimization problem:

[0061]

[0062] The iterative optimization of terrain and albedo using the alternating direction multiplier method includes the augmented Lagrangian minimization process of the variables to be solved, z and a, and the gradient descent process for the dual variable u:

[0063]

[0064] where u is the dual operator defined by the gradient descent of z, and κ represents the step size of the gradient descent defined in the optimization method;

[0065] Assume that the estimated values of the variables in the current step are a (k) , θ (k) , z (k) and the dual operator u defined based on its gradient change (k) . At the iteration step (k), the variables are updated according to the following steps:

[0066]

[0067] The navigation landmark includes:

[0068] Landmark position: The center of the navigation landmark is defined as the landmark, and its position vector relative to the celestial-fixed coordinate system of the target celestial body is defined as r lm ;

[0069] Reference coordinate system The navigation landmark coordinate system takes the landmark as the origin, the Z-axis is aligned with the landmark position vector r lm , the X-axis points east, and the Y-axis points north;

[0070] Grid size and resolution: The grid size is defined as the number of grids of the navigation landmark, and the grid resolution is defined as the resolution of each grid cell;

[0071] Relative elevation z': The relative elevation is defined as the elevation of z in the plane grid coordinates , obtained by z' = F(z). The relative elevation value is represented as an N×N matrix, where z represents the elevation in the celestial-fixed coordinate system, and F(·) represents the rotation matrix from the celestial-fixed coordinate system to the plane grid coordinates ;

[0072] Relative albedo a': The relative albedo is defined as the albedo in the plane grid coordinate system in the reference coordinate system , which is jointly determined by the relative elevation, the image, and the illumination conditions corresponding to the image used.

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

[0074] (1) The present invention introduces a local albedo consistency constraint to achieve the construction of high-resolution terrain navigation landmarks under the constraint of only a single high-resolution image. This method is implemented in a Bayesian framework. The photometric estimation is described by the likelihood term, and the prior term describes the constraints on local terrain smoothness and the assumption of local albedo consistency, avoiding ill-posed problems. For the solution of this Bayesian model, it is transformed into an optimization problem under the variational framework, and the alternating direction multiplier method is used to solve this non-convex constrained optimization problem, avoiding local non-integrable problems.

[0075] (2) The terrains of the same area constructed by the present invention using a single image under different illumination and observation conditions have good consistency, especially in areas without boulders. The evaluation based on the correlation score and position error also shows the reliability of the navigation landmarks constructed by the present invention based on a single image. BRIEF DESCRIPTION OF THE DRAWINGS

[0076] Figure 1 is a flowchart of the method of the present invention;

[0077] Figure 2 is a schematic diagram of the definition of terrain landmarks and related attributes. Among them, (2a) is a schematic diagram of the relationship between the terrain landmark coordinate system, the star-fixed coordinate, and the camera coordinate system, (2b) is the relative elevation, and (2c) is the relative albedo;

[0078] Figure 3 is an image block with different illumination and observation conditions in an embodiment. Among them, (3a)-(3d) correspond to images a-d respectively;

[0079] Figure 4 is a 25-cm resolution terrain constructed based on the SPC 40-cm resolution terrain and images a-d in an embodiment. Among them, (4a)-(4d) correspond to the terrains constructed from images a-d respectively, and (4e)-(4h) are the enlarged detail views corresponding to (4a)-(4d);

[0080] Figure 5 is an analysis result diagram of the terrain profile line 1 of the 25-cm resolution terrain constructed based on the SPC 40-cm resolution terrain and images a-d in an embodiment. Among them, (5a) is a terrain polyline diagram of the terrain profile line 1, (5b) is a schematic diagram of the image and the profile line, (5c) is a box plot of the terrain constructed from images a-d for the terrain profile line 1, (5d) and (5e) are local enlarged views of (5a), and (5f) and (5g) are local enlarged views of (5b);

[0081] Figure 6Rendering result diagrams of navigation landmarks constructed based on image d under the lighting and observation conditions of images a, b, and c in one embodiment, where (6a) is image d, (6b) is the rendering and projection result of the terrain constructed from image d under the lighting and observation conditions corresponding to image a; (6c) is the rendering and projection result of the terrain constructed from image d under the lighting and observation conditions corresponding to image b; (6d) is the rendering and projection result of the terrain constructed from image d under the lighting and observation conditions corresponding to image c;

[0082] Figure 7 For the rendering results of navigation landmarks constructed based on a single image under different lighting and observation conditions in one embodiment, as well as the correlation score and correlation position error diagrams with the corresponding navigation camera images. Detailed implementation manners

[0083] The present invention will be described in detail below with reference to the accompanying drawings and specific embodiments. This embodiment is implemented on the premise of the technical solution of the present invention, and gives detailed implementation manners and specific operation processes, but the protection scope of the present invention is not limited to the following embodiments.

[0084] Embodiment 1

[0085] This embodiment provides a photometric high-resolution terrain navigation landmark construction method with only a single image constraint, as Figure 1 shown, including the following steps:

[0086] S1, Data acquisition and preprocessing: Acquire a single high-resolution image of the target asteroid surface and low-resolution initial terrain data.

[0087] S2, Establish a terrain and albedo inversion model under the Bayesian model: Under the Bayesian framework, model the inversion problem of terrain and albedo as a problem of maximizing the posterior probability, and construct a Bayesian model considering the prior of local albedo consistency. The Bayesian model includes a likelihood term and a prior term. Among them, the likelihood term is used to describe the relationship between the photometric information of the image and the terrain and albedo, ensuring that the inversion result is consistent with the observed image, and the prior term introduces local terrain smoothness constraints and local albedo consistency assumptions to solve the ill-posed problem under a single image constraint.

[0088] In this embodiment, the specific construction process of the Bayesian model considering the prior of local albedo consistency is as follows:

[0089] S21, Recover the high-resolution terrain from the low-resolution initial terrain using photometry. Based on the photometric constraint of a single high-resolution image, invert the high-resolution terrain from the low-resolution initial terrain. Among them, the high-frequency information required for high-resolution terrain inversion, corresponding to the terrain details, is constrained by the high-resolution image, and the low-frequency information, corresponding to the geometric structure, is constrained by the initial terrain;

[0090] S22. According to Bayes' rule, the high-resolution terrain is restored by maximizing the posterior distributions of a single high-resolution image and the low-resolution terrain, and a Bayesian model is constructed based on the posterior distributions.

[0091] Given a low-resolution terrain, denote its elevation as z0, and assume a high-resolution terrain z, and that is, the high-resolution terrain can be downsampled to the low-resolution terrain. The relationship between the above low-resolution terrain and the high-resolution terrain is expressed as:

[0092] z0 = Dz + ε z (1)

[0093] where ε z represents a random process, such as measurement error, noise, and assume that is, it follows a Gaussian distribution with a mean of 0 and a variance of σ; D: Ω HR →Ω LR represents a linear downsampling operation operator.

[0094] Limited by the resolution of the initial terrain and the diversity of interpolation methods, it is ill-posed to sample from the low-resolution terrain to the high-resolution terrain only by interpolation. Therefore, a method must be found to constrain this problem and handle random noise. In fact, according to the principle of photometric measurement, for an image, the pixel intensity of each pixel in the image is related to the local terrain. Therefore, the high-resolution image can be used to constrain the process of generating the high-resolution terrain from the low-resolution terrain.

[0095] Photometry is a typical method for reconstructing the terrain shape of an object using an image. It aims to restore the shape of an object from the grayscale or color image of a single scene. At the same time, under some specific constraints, it can also invert the albedo information related to the object material from the image.

[0096] In this embodiment, the high-resolution image is used to constrain the process of generating the high-resolution terrain from the low-resolution terrain. Based on photometry, according to the SFS principle, by establishing the relationship between the image and the albedo model and minimizing it, the normal vector of any point in the image can be obtained:

[0097]

[0098] where the normal vector is a function of the slope of the current point, that is, n(x, y) = g(z x (x, y), z y (x, y)), g(·) is a function expression, z x (x, y) and z y(x, y) represents the gradient of z(x, y), where z(x, y) is the elevation of pixel (x, y), I(x, y) represents the image gray value of pixel (x, y), and R(x, y) represents the reflectivity corresponding to pixel (x, y). The reflectivity is calculated based on the lunar-Lambert reflectance model using the observation and illumination conditions at the time of the current image capture. However, it should be noted that the slope value obtained here represents the magnitude of the elevation gradient but cannot uniquely represent the direction of the slope. Therefore, the local shape is ambiguous at the singularities. To solve this ambiguity problem, additional constraints must be introduced.

[0099] Given the above brief discussion of terrain refinement based on interpolation methods and shape recovery from shadows using a single image, it can be concluded that in the context of corresponding images and low-resolution terrain, the high-frequency information required to achieve detail-preserving terrain refinement can be provided by photometric data. Similarly, the low-frequency information required to recover the shape from shadows can be conveyed by geometric data. The above process can be achieved by maximizing the posterior distribution of a single high-resolution image and low-resolution terrain. According to Bayes' rule, the posterior distribution can be expressed as:

[0100]

[0101] where, I represents the image, a represents the albedo, and P represents the probability density function. Assuming that each variable follows an independent distribution, in the proposed method, the idea of David Mumford will be adopted (see the literature: Strekalovskiy E, Chambolle A, Cremers D. A convex representation for the vectorial Mumford-Shah functional [C] / / 2012 IEEE Conference on Computer Vision and Pattern Recognition. IEEE, 2012: 1712-1719.), and an optimized variational model will be derived from the posterior distribution formula (3). Based on this optimized model, variational methods will be used to obtain the high-resolution terrain and albedo.

[0102] (1) Likelihood term

[0103] Assuming that the image I and the low-resolution terrain z0 follow independent distributions, then:

[0104] P(z0, I|z, a) = P(z0|z, a)P(I|z, a) (4)

[0105] Assuming that the low-resolution terrain z0 is independent of the surface albedo a, then:

[0106] P(z0, I|z, a) = P(z0|z)P(I|z, a) (5)

[0107] Also according to Equation (1), the noise i.e., it follows a Gaussian distribution with a mean of 0 and a variance of σ z then there is:[[]]

[0108]

[0109] where represents the two-norm operator, and D(·) represents the linear downsampling operation;

[0110] Furthermore, establish the relationship between the image I and the reflectance R: Considering the self-shadow characteristic of the terrain, express the reflectance R as a function of the normal vector corresponding to the current terrain, then there is:[[]]

[0111] I(x, y) = Λ × a(x, y) × max{R(n(x, y), s(x, y), v(x, y)), 0} + ε I (7)

[0112] where (x, y) represents the pixel index in the image, the max{·} operation restricts the self-shadow characteristic, that is, when the object is blocked or should not receive light, the self-shadow factor is zero; s(x, y) and v(x, y) are the unit vectors of the incident light and the outgoing light at (x, y) respectively, n(x, y) is the normal vector at (x, y), Λ represents the parameter related to the image exposure, and there are many factors affecting Λ. During the terrain optimization process, it is usually regarded as a constant, a(x, y) represents the albedo at (x, y), and ε I represents the image noise and is obtained through pre-calibration.

[0113] Before estimating the terrain and albedo, Λ is initialized to a constant, that is:[[]]

[0114]

[0115] where x and y respectively define the integration intervals in the row and column directions of the image;

[0116] Also assume i.e., it follows a Gaussian distribution with a mean of 0 and a variance of σ I then there is:[[]]

[0117]

[0118] When the incident light, the outgoing light vector, and the camera internal parameters are known, whether it is orthographic projection or projective projection, n(x, y) is essentially a function of, then Equation (9) above can also be expressed as:[[]]

[0119]

[0120] (2) Prior term

[0121] Assume that the terrain and albedo are independent of each other (the independence of albedo from the terrain comes from the Lambertian assumption), then we have:

[0122] P(z,a) = P(z)P(a) (11)

[0123] Regarding the constraint on the terrain z, the smoothness of the local terrain is mainly considered. The smoothness constraint on the terrain is necessary, which is conducive to obtaining a correct terrain. Therefore, the present invention also constructs a smoothness constraint on the terrain. The smoothness constraint based only on the first-order derivative of the terrain does not have edge protection and is prone to over-smoothing of the edges, resulting in inappropriate processing of the edges of stones. This embodiment uses second-order regularization, which combines the edge protection characteristics of the second-order coupling model and the directional image information to optimize the terrain.

[0124]

[0125] Where G is the Laplacian kernel function, w z represents the bandwidth of the kernel function, which is used to control the magnitude of smoothness, Δz represents the second-order derivative of the terrain z, and μ is used to control the dispersion degree of the distribution, which is set as a constant.

[0126] Regarding the constraint on the albedo a, the local area albedo consistency constraint is adopted. In the case of only one image, unless additional information is provided, such as the integrability boundary assumption or the statistical results of other studies, solving the albedo and reflectivity simultaneously in one observation is a mathematically ill-posed problem.

[0127] During the process of optimizing the terrain, the assumptions about the albedo include: assume that in the current image, the albedo is 1.0. This assumption is reasonable for optimizing the terrain of a small local area on the large-scale planetary surface, because it is considered that the substances on the surface of the local area do not differ much, so the albedo does not change much; another assumption adopted by various studies is "constant albedo in the local area", that is, it is assumed that the albedo of a point on the surface is highly similar to its neighborhood. Similarly, this assumption is logical and also holds for asteroids. This embodiment adopts the second assumption, that is, the assumption of constant albedo in the local area. The difference is that this embodiment not only considers the spatial distribution of pixels, but also comprehensively considers the spectral reflection characteristics of pixels. Therefore, the method of superpixel segmentation is used for albedo consistency constraint, that is, it is assumed that each superpixel has consistent albedo characteristics.

[0128] For the superpixel segmentation results, they can be represented by an indicator function, i.e.:

[0129]

[0130] If it is considered that the above process follows the Pott model, then there is:

[0131]

[0132] where the parameter λ > 0, Ω i ∈Ω HR , i is the superpixel index, n is the total number of superpixels, denotes the indicator function, |·| represents the distance defined by the superpixel, and λ is used to control the discreteness of this distribution, which is set as a constant.

[0133] S3. Terrain and albedo inversion: Convert the Bayesian model into an optimization problem under the variational framework, and use the alternating direction multiplier method to iteratively optimize the terrain and albedo, gradually approaching the optimal solution. In each iteration, update the terrain and albedo respectively, and introduce Lagrange multipliers to handle the constraint conditions, where the constraint conditions include local albedo consistency constraints and terrain smoothness constraints.

[0134] According to the above analysis, the estimation of the optimized terrain and albedo can ultimately be achieved by maximizing the posterior distribution formula (3), or rather, minimizing its negative logarithm, then the following energy function can be constructed:

[0135]

[0136] where D I (z, a) expresses the photometric-based assumption, that is, the relationship between pixel intensity, albedo, and the reflectance model; represents the constraint of the initial terrain; to solve the ambiguity between unknowns, regularization terms for the terrain and albedo are considered respectively, that is, L z (z) and L a (a); β, γ z and γ a correspond to the weights of each data term and regularization term respectively. The specific mathematical expressions of the above items are as follows:

[0137]

[0138] where, in equation (17), according to different projection methods, the expression of the normal vector will be different. In this section, the orthographic projection method is adopted, and at this time, the normal vector can be expressed as:

[0139]

[0140] where z xand z y respectively represent the gradients of the terrain z in two orthogonal directions.

[0141] In addition, another key element in Equation (16) is the reflectivity model. Here, the Lunar-Lambert reflectivity model is adopted.

[0142] Express the energy function of Equation (16) in the form of an optimization problem under the variational framework, and obtain:

[0143]

[0144] Obviously, the variational problem in the above equation is not only non-convex, but also the gradient of the terrain z is non-linearly expressed. In this embodiment, z will be directly solved instead of using the method based on gradient integration. By introducing an auxiliary vector field to achieve this, re-express the optimization problem under the above variational framework as an equivalent constrained optimization problem:

[0145]

[0146] To solve the above non-convex constrained optimization problem, use ADMM (Alternating Direction Method of Multipliers) to solve it. This can be reduced to an iterative process, including the augmented Lagrangian minimization process of the variables z and a to be solved, and the gradient descent process of the dual variable u:

[0147]

[0148] where u is the dual operator of z defined based on gradient descent, and κ represents the step size of gradient descent defined in the optimization method.

[0149] Assume that the estimated values of each variable in the current step are a (k) , θ (k) , z (k) and the dual operator u defined based on its gradient change (k) , at the iteration step (k), each variable is updated according to the following steps:

[0150]

[0151]

[0152] S4. Navigation landmark construction: Based on the retrieved high-resolution terrain and albedo data, construct terrain navigation landmarks.

[0153] The terrain navigation landmark is defined as a local area on the surface of the target celestial body with a specific size, represented as a grid with a specific size. The construction of the terrain landmark is completed during the ground work before the landing mission. Specifically, the navigation landmark is defined by the landmark position, the reference coordinate system, the grid size, and the relative elevation and relative albedo maps.

[0154] (i) Landmark position: The center (origin) of the navigation landmark is defined as the landmark, and the position vector of it relative to the celestial-fixed coordinate system of the target celestial body is defined as r lm ; as shown in (2a) of Figure 2 . It should be noted that the landmark itself is not necessarily a prominent feature, such as a boulder or a crater. The relative topography and albedo of the surrounding terrain are the basis of the navigation landmark required in the method.

[0155] (ii) Reference coordinate system The navigation landmark coordinate system takes the landmark as the origin, the Z-axis is aligned with the landmark position vector r lm , the X-axis points east, and the Y-axis points north.

[0156] (iii) Grid size and resolution: The grid size is defined as the number of grids of the navigation landmark (for example: N×N), and the grid resolution is defined as the resolution of each grid cell (for example: meters / grid). It should be noted that the grid size is not the same as the pixel size. The ground resolutions of different terrain landmarks in the landmark database are different.

[0157] (iv) Relative elevation z’: The relative elevation is defined as the elevation of z in the plane grid coordinates , obtained by z’ = F(z). The relative elevation value is represented as an N×N matrix, where z represents the elevation in the celestial-fixed coordinate system, and F(·) represents the rotation matrix from the celestial-fixed coordinate system to the plane grid coordinates , as shown in (2b) of Figure 2 .

[0158] (v) Relative albedo a’: The relative albedo is defined as the albedo in the plane grid coordinate system in the reference coordinate system , which is jointly determined by the relative elevation, the image, and the illumination conditions corresponding to the used image. It should be noted that it describes the surface material brightness relative to all other elements in this navigation landmark. The relative albedo cannot be used as a description of the inherent surface brightness, as shown in (2c) of Figure 2 .

[0159] Example 2

[0160] This example conducts experimental verification on the method of Example 1.

[0161] (1) Experimental data

[0162] The experimental data used in this embodiment are the images of the Bennu asteroid taken during the OSIRIS-Rex exploration mission, as shown in (3a)-(3d) of Figure 3 . These image blocks correspond to different lighting and observation conditions, and the specific information is shown in Table 1. These images are named Image a, Image b, Image c, and Image d respectively. The resolution of the images is approximately 25 cm. Figure 3 The direction vectors of the incident and outgoing light rays of the experimental images in Table 1

[0163] Table 1

[0164]

[0165] (2) Evaluation metrics

[0166] Topographies of the same area are constructed using images under different lighting and observation conditions. Generally speaking, if the proposed method is robust, then the topographies and albedos constructed using images under different lighting and observation conditions are almost the same. Therefore, the method proposed in the present invention is evaluated by comparing these topographies with each other. To evaluate the consistency or difference between different topographies, a quantitative evaluation method based on boxplot metrics is proposed.

[0167] A boxplot is an effective statistical tool for analyzing data dispersion. The important metrics it includes are the maximum value, minimum value, median, upper and lower quartiles, and outliers of a set of data. Based on these metrics, we propose to use the Mean Inter - quartile Range (MIQR) to evaluate the differences between topographies. MIQR is defined as

[0168]

[0169] where Q1 represents the first quartile, also known as the "lower quartile", which is equal to the 25% number after arranging all the values in the sample from small to large, and Q3 represents the third quartile, also known as the "upper quartile", which is equal to the 75% number after arranging all the values in the sample from small to large, and N represents the number of groups in the sample. In terms of the display of the boxplot, the Inter - quartile Range (IQR) reflects the size of the box body and represents the central tendency and distribution of the data, that is, the smaller the IQR, the more concentrated the data distribution within the current group.

[0170] (3) Experimental results

[0171] The input topography is the topography data with a resolution of 40 cm of SPC, and the input images are the above 4 images. The topography constructed by the proposed method is as shown in Figure 4 . Figure 4 as shown. Figure 4The coordinates of the terrain shown in (4a)-(4d) are in the fixed coordinate system of Bennu. To better display the details of terrain generation for comparison, the generated topographic maps are enlarged, as shown in Figure 4 (4e)-(4f) respectively. Judging from the visual results, the terrains generated based on these images are very close. The areas with differences are mainly located in the regions where the boundaries between the stones and their shadows are, as shown by the red circles in Figure 4 . Such results are reasonable because due to different lighting and observation conditions, the directions and sizes of the shadows of the stones in the images are different, which leads to different gray values of the corresponding regions in the images. However, since the initial terrains input are the same, according to the processing method for the shadow regions proposed, the initial terrain will play a dominant role at this time. Therefore, even though there are differences in pixel gray values, they still do not differ much from each other.

[0172] Furthermore, the differences between the four constructed terrains are further quantitatively evaluated through profile lines. The first profile line is selected in a relatively flat area of terrain change, as shown by the solid line area in Figure 5 (5b). The front part of this profile line mainly passes through flat areas. As shown in Figure 5 (5a), among the first 62 sampling points, the changes of the four terrains are very close. In the local enlarged view shown in Figure 5 (5d), it can be seen that for this part of the area, the maximum difference between the terrains does not exceed 0.3 m. According to the box plot shown in Figure 5 (5f), it can be seen that the maximum difference in IQR does not exceed 0.18 m, and for each data point, there are no outliers. The MIQR of the first 62 sampling points is shown in Table 2 as 0.1606. All of the above indicators illustrate the consistency of the four terrains. It also shows that for terrains constructed from images with different lighting and observation conditions in flat areas, the differences are not significant. The tail of this profile line passes through a boulder. Due to the influence of lighting, the texture distributions of the images around this boulder are somewhat different in the four images. As shown in Figure 5 (5a), starting from the 63rd sampling point, the differences between the four terrains begin to increase, and the maximum terrain deviation is 1.8 m, which occurs between image a and image d. According to the images shown in Figure 4 , it can be seen that the "stone-shadow" texture distributions of these two images are quite different. According to the principle of the method, in the shadow region, its initial terrain and its neighborhood will play a major role in constructing the terrain. Therefore, due to different shadow-covered areas and large terrain changes, when optimizing the new terrain, the input initial terrains and their neighborhoods are quite different. Therefore, the newly constructed terrains are quite different. According to Figure 5(5g) It can be seen that the boxes corresponding to the 71st and 72nd data are the largest, approximately 1.6 m. The MIQR of the 63rd - 77th sampling points is shown in Table 2 as 1.0693. However, in terms of outliers, there are still no outliers in these data, indicating the consistency of the terrain constructed by the four images within a certain range.

[0173] Analysis of MIQR Results of Profile Line 1 and Profile Line 2

[0174]

[0175]

[0176] This embodiment further discusses the effectiveness of the navigation landmarks generated based on a single image. As Figure 6 shown, the rendering results of the navigation landmarks constructed based on Image d under the lighting and observation conditions of Image a, Image b, and Image c are presented. It can be seen that under different lighting and observation conditions, image blocks with texture distributions consistent with the corresponding navigation camera images can be rendered based on this terrain. This demonstrates the accuracy of the relative elevation of the terrain constructed by the method proposed in the present invention. It should be noted that during the process of constructing the terrain based on a single image, the shadow areas in the image are ineffective information for terrain optimization. Therefore, the construction of the terrain in the shadow areas will be affected. Therefore, when constructing the terrain based on a single image, an image with good lighting and observation angles should be selected.

[0177] One of the important methods for evaluating the correctness of navigation landmarks is to compare the rendered navigation landmarks with the features in the corresponding images. Specifically, the navigation landmarks are rendered according to the lighting and observation geometric conditions of each navigation camera image, and then "image correlation" is performed with the image. The results include the correlation score and the error between the detected correlation peak position and the expected feature position in the image, that is, the correlation position deviation. Theoretically, for a high-quality navigation landmark, the main peak in its correlation results has a high correlation score, which is far from the correlation scores of any secondary peaks or false peaks. Generally speaking, if the correlation position deviation of a ground object is very small, it means that the position of the navigation landmark on the surface is accurate. As Figure 7 shown, the correlation results of the navigation landmarks generated by the proposed method and the navigation camera image set are presented. Here, the lighting and observation conditions of each image are from the ephemeris data matched with the corresponding image. According to Figure 7 the results in, it can be seen that in the correlation results, the correlation scores of the main peaks of most correlation results are greater than the threshold th cor , and the average main peak position deviation is less than th loc . Here, the threshold th cor = 4.5, th loc= 1.5 pixels, referring to the parameters in the initial task test of NFT. For instances of related errors, the score related to the main peak is less than th cor , and the position deviation is extremely large. We believe that this kind of error should be due to the fact that under a certain lighting condition, there will be large shadows projected by boulders in the navigation camera image, while there is none in the navigation landmark rendering result, which leads to related errors. Another possible reason is that the features are distributed at the edge of the navigation camera image. At this time, the feature texture projected by the navigation landmark into the image is relatively distorted, which is not conducive to accurate correlation. Such results indicate that during the actual process of constructing navigation landmark data, it is necessary to specifically evaluate the adaptability of navigation landmarks to various possible lighting conditions and the position distribution of navigation landmarks in the image.

[0178] In summary, the present invention proposes a photometric high-resolution navigation landmark construction method based on a single image constraint under a Bayesian framework. It is achieved by maximizing the posterior distributions of a single high-resolution image and a low-resolution terrain. The introduced constraints include the smoothness of the local terrain and the consistency of the local albedo, which solves the ill-posed problem of constructing the terrain with only a single image. The experimental results show that the proposed algorithm has good consistency for the terrain of the same area constructed using single images under different lighting and observation conditions, especially for flat areas. Based on the evaluation of the correlation score and position error, it shows that the navigation landmarks constructed based on a single image in the proposed method are also reliable, but their accuracy is slightly worse than that of the navigation landmarks constructed under the constraint of multiple images.

[0179] The preferred specific embodiments of the present invention have been described in detail above. It should be understood that those of ordinary skill in the art can make many modifications and variations based on the concept of the present invention without creative labor. Therefore, all technical solutions that can be obtained by those skilled in the art in the technical field based on the concept of the present invention through logical analysis, reasoning, or limited experiments on the basis of the prior art should be within the protection scope determined by the claims.

Claims

1. A photometric high-resolution terrain navigation landmark construction method with only single-image constraint, characterized in that It includes the following steps: Data acquisition and preprocessing: Obtain a single high-resolution image of the target asteroid's surface and low-resolution initial terrain data; Establish a terrain and albedo inversion model under the Bayesian model: Under the Bayesian framework, model the inversion problem of terrain and albedo as a problem of maximizing the posterior probability, and construct a Bayesian model considering the prior of local albedo consistency. The Bayesian model includes a likelihood term and a prior term. Among them, the likelihood term is used to describe the relationship between the photometric information of the image and the terrain and albedo, ensuring that the inversion result is consistent with the observed image. The prior term introduces local terrain smoothness constraints and local albedo consistency assumptions to solve the ill-posed problem under the constraint of a single image; Terrain and albedo inversion: Convert the Bayesian model into an optimization problem under the variational framework, and use the alternating direction multiplier method to iteratively optimize the terrain and albedo, gradually approaching the optimal solution. In each iteration, update the terrain and albedo respectively, and introduce Lagrange multipliers to handle the constraint conditions, where the constraint conditions include local albedo consistency constraints and terrain smoothness constraints; Navigation landmark construction: Based on the inverted high-resolution terrain and albedo data, construct terrain navigation landmarks. The navigation landmarks are defined by landmark positions, reference coordinate systems, grid sizes, and relative elevation and albedo maps. Each navigation landmark represents a local area on the surface of the target celestial body.

2. A method for constructing a photometric high-resolution terrain navigation landmark with only single-image constraint according to claim 1, characterized in that The specific construction process of the Bayesian model considering the prior of local albedo consistency is as follows: Use photometry to recover the high-resolution terrain from the low-resolution initial terrain. Based on the photometric constraints of a single high-resolution image, invert the high-resolution terrain from the low-resolution initial terrain. Among them, the high-frequency information required in the high-resolution terrain inversion corresponds to the terrain details constrained by the high-resolution image, and the low-frequency information corresponds to the geometric structure constrained by the initial terrain; According to the Bayesian rule, recover the high-resolution terrain by maximizing the posterior distribution of the single high-resolution image and the low-resolution terrain, and construct a Bayesian model based on the posterior distribution.

3. A method for constructing a photometric high-resolution terrain navigation landmark with only single-image constraint according to claim 2, characterized in that Given a low-resolution terrain with its elevation denoted as z0, and assume a high-resolution terrain z, and the relationship between the low-resolution terrain and the high-resolution terrain is expressed as: z0 = Dz + ε z where ε z represents a random process; D: Ω HR → Ω LR denotes a linear downsampling operation operator; Use the high-resolution image to constrain the process of generating the high-resolution from the low-resolution terrain. Based on photometry, according to the SFS principle, establish the relationship between the image and the albedo model and minimize it to obtain the normal vector of any point in the image: Among them, the normal vector is a function of the slope of the current point, that is, n(x, y) = g(z x (x, y), z y (x, y)), g(·) is a function expression, z x (x, y) and z y (x, y) represents the gradient of z(x, y), z(x, y) is the elevation of the pixel (x, y), I(x, y) represents the image gray value of the pixel (x, y), R(x, y) represents the reflectivity corresponding to the pixel (x, y), and the reflectivity is calculated based on the lunar-Lambert reflectivity model using the observation and illumination conditions at the time of shooting the current image; the slope characterizes the gradient magnitude of the elevation.

4. A method for constructing a photometric high-resolution terrain navigation landmark with only single-image constraint according to claim 3, characterized in that The posterior distribution is expressed as: Among them, I represents an image, a represents albedo, and P represents a probability density function.

5. A method for constructing a photometric high-resolution terrain navigation landmark with only single-image constraint according to claim 4, characterized in that, The construction process of the likelihood term includes the following steps: Assume that the image I and the low-resolution terrain z0 follow independent distributions, then: P(z0,I|z,a) = P(z0|z,a)P(I|z,a) Assume that the low-resolution terrain z0 is independent of the surface albedo a, then: P(z0,I|z,a) = P(z0|z)P(I|z,a) Assumed noise That is, it follows a Gaussian distribution with a mean of 0 and a variance of σ z Then there is: Among them, represents the two-norm operator, and D(·) represents the linear downsampling operation; Construct the relationship between the image I and the reflectance R: Considering the self-shadow characteristics of the terrain, express the reflectance R as a function of the normal vector corresponding to the current terrain, then: I(x,y) = Λ × a(x,y) × max{R(n(x,y), s(x,y), v(x,y)), 0} + ε I Among them, (x, y) represents the pixel index in the image. The max{·} operation constrains the self-shadow characteristic, that is, when the object is occluded or should not receive light, the self-shadow factor is zero; s(x, y) and v(x, y) are the unit vectors of the incident light and the outgoing light at (x, y) respectively, n(x, y) is the normal vector at (x, y), Λ represents the parameter related to the image exposure, a(x, y) represents the albedo at (x, y), and ε I represents the image noise, which is obtained through pre-calibration; Before estimating the terrain and albedo, Λ is initialized to a constant, that is: Among them, x and y respectively define the integration intervals in the row and column directions of the image; Hypothesis That is, it follows a Gaussian distribution with a mean of 0 and a variance of σ I Then we have: When the incident light ray, the outgoing light ray vector, and the camera internal parameters are known, whether it is orthographic projection or perspective projection, n(x, y) is essentially a function of, then:

6. A method for constructing a photometric high-resolution terrain navigation landmark with only single-image constraint according to claim 5, characterized in that, The construction process of the prior term includes the following steps: Assume that the terrain and the reflectance are independent of each other, then: P(z,a) = P(z)P(a) For the constraint of the terrain \(z\), a second-order regularization is adopted to constrain the smoothness of the local terrain, and terrain optimization is achieved by combining the edge protection characteristics of the second-order coupling model and the directional image information: where G is the Laplacian kernel function, w z represents the bandwidth of the kernel function, which is used to control the magnitude of smoothness, Δz represents the second derivative of the terrain z, and μ is used to control the dispersion degree of the distribution and is set as a constant; For the constraint of the albedo \(a\), a local area albedo consistency constraint is adopted, that is, it is assumed that each superpixel has a consistent albedo characteristic, which is realized by using the method of superpixel segmentation. For the superpixel segmentation result, it is represented by an indicator function, that is: Assuming that the above process follows the Pott model, then we have: where the parameter λ > 0, and Ω i ∈ Ω HR , i is the superpixel index, n is the total number of superpixels, denotes the indicator function, |·| represents the distance defined by the superpixel, and λ is used to control the discreteness of the distribution and is set as a constant.

7. A method for constructing a photometric high-resolution terrain navigation landmark with only single-image constraint according to claim 6, characterized in that, The optimal estimation of the terrain and albedo is achieved by maximizing the posterior distribution, that is, by minimizing its negative logarithm. Then, the following energy function is constructed: Among them, D I (z,a) expresses the photometric-based assumption of the relationship between pixel intensity, albedo, and reflectance model; represents the constraint of the initial terrain; L z (z) and L a (a) are regularization terms for terrain and albedo respectively; β, γ z and γ a correspond to the weights of each data term and regularization term respectively; 8. A method for constructing a photometric high-resolution terrain navigation landmark with only single-image constraint according to claim 7, characterized in that, Express the energy function in the form of an optimization problem under the variational framework, and we get: By introducing an auxiliary vector field to directly solve for z, the optimization problem under the above variational framework is reformulated as an equivalent constrained optimization problem:

9. A method for constructing a photometric high-resolution terrain navigation landmark with only single-image constraint according to claim 8, characterized in that The iterative optimization of the terrain and albedo using the alternating direction multiplier method includes the augmented Lagrangian minimization process of the variables \(z\) and \(a\) to be solved, and the gradient descent process of the dual variable \(u\): where \(u\) is the dual operator defined by \(z\) based on gradient descent, and \(\kappa\) represents the step size of gradient descent defined in the optimization method; Assume that the estimated values of the variables in the current step are a (k) , θ (k) , z (k) and the dual operator u defined based on its gradient change (k) . At the iteration step (k), the variables are updated according to the following steps:

10. A method for constructing a photometric high-resolution terrain navigation landmark with only single-image constraint according to claim 1, characterized in that, The navigation landmarks include: Landmark position: The center of the navigation landmark is defined as the landmark, and its position vector relative to the celestial-fixed coordinate system of the target celestial body is defined as r lm ; Reference coordinate system The navigation landmark coordinate system has the landmark as the origin, and the Z-axis is aligned with the landmark position vector r lm The X-axis points east and the Y-axis points north; Grid size and resolution: The grid size is defined as the number of grids of the navigation landmark, and the grid resolution is defined as the resolution of each grid cell; Relative elevation z': The relative elevation is defined as the elevation in the plane grid coordinates of z, obtained by z' = F(z). The relative elevation value is represented as an N×N matrix, where z represents the elevation in the star-fixed coordinate system, and F(·) represents the rotation matrix from the star-fixed coordinate system to the plane grid coordinates ; ​ Relative albedo a': The relative albedo is defined as the albedo in the grid coordinate system of the midplane in the reference coordinate system and is jointly determined by the relative elevation, the image, and the illumination conditions corresponding to the image used.

Images