Asteroid surface pixel-level high-precision terrain construction method based on improved stereophotometry

Abstract

The invention relates to an asteroid surface pixel-level high-precision terrain construction method based on an improved stereophotometry, and the method comprises the following steps: obtaining asteroid multi-view image data which are obtained through shooting of a detector under different illumination conditions and visual angles; establishing a multi-view stereo photometric method variation model, reconstructing a terrain height field by the model through minimizing a photometric error, and introducing a multi-view stereo confidence constraint, a non-convex estimator and a terrain edge protection mechanism to construct a target function; the method comprises the following steps: converting a variable updating problem of a multi-view stereophotometry variation model into an independent nonlinear optimization problem, decomposing an objective function into a plurality of sub-problems through a Lagrange operator, and performing efficient solution by adopting a quasi-Newton method; in the iteration process, the terrain resolution is gradually optimized, and finally pixel-level high-precision terrain data are output. Compared with the prior art, the method has the advantages of being capable of adapting to different terrains to achieve high-precision terrain construction and the like.

Description

Technical Field

[0001] The present invention relates to the field of planetary surface terrain construction, and particularly to a high-precision terrain construction method for asteroid surface pixels based on an improved stereophotometry method. Background Art

[0002] High-resolution terrain is crucial for scientific exploration and sample collection activities in asteroid exploration missions. For the global digital elevation model, it plays an important role in revealing the geological origin and evolution process of asteroids. For example, by accurately measuring the volume and estimating the mass to infer the density of asteroids, and the global slope, elevation, and surface roughness maps provide quantitative data for in-depth understanding of the surface processes affecting the evolution of the regolith. The high-resolution elevation, slope, and the geometric height of boulders within the sampling ellipse in the sampling area are crucial for scientists to assess potential risks during the sampling process and ensure the safety and efficiency of the sampling activities.

[0003] Currently, the high-resolution three-dimensional terrain construction methods mainly include the binocular disparity method (BD), the perspective projection transform method (PPT), structure from motion (SfM), shape from shading (SfS), and stereophotoclinometry (SPC), etc. Among them, for the high-resolution three-dimensional terrain construction in the deep space exploration environment, currently, SfS and SPC are mainly used. The theoretical basis of SfS and SPC is the photometry method.

[0004] The typical photometric method for terrain construction is divided into two steps. First, a cost function is minimized to find the surface normal, and then, a second cost function is minimized to recover the surface height from the normal. However, such a process is vulnerable to the situation where outliers are non-integrable, making it impossible to solve the terrain. To address this problem, some existing techniques have attempted methods for directly solving the surface height, only calculating the surface normal as needed during the optimization process. The results show that the method of directly solving the surface height can, to a certain extent, solve the problem of non-integrability of the surface based on the gradient integration method. To improve the convergence speed, avoid local minima, and handle camera position and pointing errors, many scholars have used a multi-resolution method based on multiple images to solve the unknown terrain and albedo, etc. The literature "Improved methods of estimating shape from shading using the lightsource coordinate system" proposed a method for constructing the lunar surface terrain based on multi-image photometry. This method takes into account the lunar-Lambert reflectance model and image exposure, and preprocesses the image set through bundle adjustment to reduce the errors caused by camera pose information. However, this method does not constrain the terrain smoothness, so the obtained terrain has the problem of discontinuous changes. In the literature "Shape from shading using linear approximation", a very rough digital terrain model (DTM) jointly produced by the Lunar Orbiter Laser Altimeter (LOLA) and LROC (Lunar Reconnaissance Orbiter Camera) images was used as the input of the initial terrain, with a ground sampling distance of about 100 meters per pixel. Their final result has an elevation error of about 30 meters. This method directly discretizes the surface normal parameters through finite differences, and this local linearization strategy will amplify the noise interference when the input DTM resolution is low (100 meters / pixel).

[0005] In summary, the traditional photometric method estimates the gradient field of the object surface rather than the surface shape itself. Its essence is to obtain the surface shape by integrating the gradient field. However, due to the presence of noise such as shadows or highlights in the image, on the one hand, it will lead to a large error in the calculated gradient field, and on the other hand, the process from gradient integration to surface shape will inevitably encounter the situation where the selected integration path is non-integrable. The above problems may all cause the surface shape solving process to fall into a local optimal solution, and it is impossible to obtain accurate geometric shape information of the surface. In addition, for the existing high-resolution terrain construction methods, due to the lack of consideration of the possible "noise" in the image, such as the scattered shadow areas in the image, there are also many non-smooth details in the constructed terrain. Such characteristics are particularly obvious in the images of extraterrestrial celestial bodies in the form of rubble piles. Summary of the Invention

[0006] The object of the present invention is to provide a method for constructing a high-precision terrain at the pixel level on the surface of an asteroid based on an improved stereo-photometric method, which directly solves the terrain under a variational framework by combining low-resolution terrain and high-resolution image sets, introduces a non-convex estimator to enhance the photometric robustness constraint and a multi-view stereo confidence constraint to enhance the robustness of the method to noise, and introduces a terrain edge protection operation to solve the problem of over-smoothing of the terrain edge in the current terrain optimization method, thereby improving the accuracy of the constructed terrain.

[0007] The object of the present invention can be achieved by the following technical solutions:

[0008] A method for constructing a high-precision terrain at the pixel level on the surface of an asteroid based on an improved stereo-photometric method, comprising the following steps:

[0009] Data acquisition: Obtain multi-view image data of the asteroid, where the image data is obtained by the detector under different illumination conditions and viewing angles;

[0010] Model construction: Establish a variational model of multi-view stereo-photometry, which reconstructs the terrain height field by minimizing the photometric error, and introduce a multi-view stereo confidence constraint, a non-convex estimator and a terrain edge protection mechanism to construct the objective function;

[0011] Model solution: Convert the variable update problem of the multi-view stereo-photometry variational model into an independent non-linear optimization problem, decompose the objective function into multiple sub-problems through the Lagrangian operator, and use the quasi-Newton method for efficient solution; during the iteration process, gradually optimize the terrain resolution, and finally output high-precision terrain data at the pixel level.

[0012] The objective function is expressed as:

[0013] F(z,a) = D image (z) + υ·D sfs (z,a) + α1·L edge (z) + α2·L dem (z)

[0014] This objective function consists of two data terms and two regularization terms. Among them, D image is the multi-view stereo data term, representing the multi-view stereo confidence constraint, which is constructed by using a depth parameterized model and the photometric consistency of the projected surface points; D sfs is the stereo-photometric data term, which is constructed by using the lunar-Lambert photometric model and the difference between the image rendered based on the photometric model and the real image; L egde represents the terrain edge regularization term, which is used to constrain the continuity of the terrain edge; L demRepresents the terrain regularization term, which is used to constrain the correctness of the terrain geometric structure; ν, α1, and α2 are weights; z is the depth corresponding to the terrain in the camera space coordinate system, and a is the albedo.

[0015] Combining multi-view stereo and photometric constraints to optimize the low-resolution terrain into a terrain with pixel-level accuracy. The multi-view stereo data term D image Is expressed as:

[0016]

[0017] Where, Defines the brightness difference between the corresponding points in the reference image I0(x,y) and other views I k (x,y), Defines the gradient difference between the corresponding points in the reference image I0(x,y) and other views I k (x,y); (x,y) represents the pixel position in the image, k represents the view index, n represents the total number of views, Ω0 represents the reference image domain, J(·) represents the Jacobian matrix, and γ represents the weighting factor, which is used to balance the contributions of the two terms to D image Of, Represents the L2 norm, Represents the Frobenius norm, and Γ D (·) represents the Geman-McClure non-convex estimator.

[0018] The stereo photometric data term D sfs Is expressed as:

[0019]

[0020] Where, k represents the view index, I represents the image, Ω0 represents the reference image domain, a(x,y) represents the albedo corresponding to the pixel at (x,y) in the image, z(x,y) represents the optimized high-resolution terrain, Represents the L2 norm, R k Represents the lunar-Lambert photometric model, T k Represents the parameter related to the image exposure, and Ψ(·) is the Cauchy non-convex estimator.

[0021] The terrain edge regularization term L egde Is expressed as:

[0022]

[0023] Where, z(x,y) represents the optimized high-resolution terrain, Ω0 represents the reference image domain, Represents the Frobenius norm, Ψ(·) is the Cauchy non-convex estimator, and Hess(·) represents the smoothness constraint on the terrain edge, which is defined as:

[0024]

[0025] where z xx represents the second derivative of the terrain z in the x direction, and z yy represents the second derivative of the terrain z in the y direction, and z xy represents the gradient of the terrain z in the x direction, and z x is the derivative in the y direction.

[0026] The terrain regularization term L dem is defined as:

[0027]

[0028] where z0(x, y) represents the initial low-resolution terrain, z(x, y) represents the optimized high-resolution terrain, Ω0 represents the reference image domain, represents the L2 norm.

[0029] The Geman-McClure non-convex estimator is defined as:

[0030]

[0031] where s represents the independent variable of the function input, and ε is defined as a positive constant close to 0 to ensure the differentiability of the function at s = 0.

[0032] The Cauchy non-convex estimator is defined as:

[0033] ψ(s) = λ 2 (1 + log(1 + s 2 / λ 2 ))

[0034] where s represents the independent variable of the function input, and λ is an empirical parameter.

[0035] In the model solution, by introducing an auxiliary vector field to directly solve for the terrain z, the solution processes of the terrain and albedo are expressed as an equivalent constrained optimization problem:

[0036]

[0037] where θ is the Lagrangian operator, represents the gradient of the terrain z, a represents the albedo, represents the terrain, albedo, and Lagrangian operator obtained by the optimized solution.

[0038] In the model solution, 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) , assuming the initial value of the reflectivity is 1, i.e., uniform white, at the iteration step (t), each variable is updated according to the following steps:

[0039]

[0040] where k represents the view index, I0(x, y) is the reference image, and I k (x, y) is other views, represents the L2 norm, represents the Frobenius norm, Γ D (·) represents the Geman-McClure non-convex estimator, R k represents the moon-Lambert photometric model, T k represents the parameter related to image exposure, Ψ(·) is the Cauchy non-convex estimator, γ, β, and κ all represent weight parameters, defined as empirical constants, Hess(·) represents the smoothness constraint on the terrain edge, z0 represents the initial low-resolution terrain; the update of the Lagrangian operator θ is a series of independent non-linear optimization problems, solved using the implementation of the L-BFGS method; the update of the dual operator u is a non-linear optimization sub-problem, independently solved using the L-BFGS method for each pixel; the update of the terrain z uses the conjugate gradient of the normal equation to solve this problem; the solution of the albedo a is a linear least squares problem, solved using the gradient descent method.

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

[0042] (1) In the traditional two-step method, the normal is optimized first and then the terrain is integrated. It is easy to fail in solving due to the non-integrability of the gradient field. The present invention constructs a variational model of multi-view stereo photometry, directly solves the terrain height field as an optimization variable, bypassing the gradient integration step. At the same time, a multi-view stereo data term is introduced, and the photometric consistency of multi-view images is used to directly constrain the terrain height, rather than indirectly relying on normal integration. The problem of non-integrability of the gradient field is avoided, and the terrain geometry is directly recovered from the image data, improving the stability and accuracy of terrain solving.

[0043] (2) The existing methods do not constrain the terrain smoothness, resulting in discontinuities in the reconstructed terrain. The present invention introduces a terrain regularization term and an edge protection mechanism, which retain edge details while constraining the global smoothness of the terrain, solving the problem of over-smoothing or discontinuity in the terrain mutation area in the traditional method.

[0044] (3) Traditional methods are sensitive to image noise, resulting in distorted terrain details. The present invention introduces multi-view stereo confidence constraints, combines brightness differences and gradient differences, and uses multi-view cross-validation to exclude single-view noise interference. A non-convex estimator is introduced to enhance the adaptability to areas with abnormal illumination, and the accuracy of terrain details can still be maintained in the presence of noise such as shadows and highlights, improving the reconstruction robustness under complex illumination conditions. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0046] Figure 2 is a schematic diagram of 1D data for integrating the noise gradient field by least squares and sparse representation on a regular grid;

[0047] Figure 3 is a schematic diagram of the experimental area in one embodiment;

[0048] Figure 4 are the terrains of the experimental area with OLA80cm and SPC40cm resolutions and the terrains with 40cm and 25cm resolutions generated by the method of the present invention in one embodiment. Among them, (4a) and (4c) are the terrains with OLA80cm and SPC40cm resolutions respectively, (4b) and (4d) are the terrains with 40cm and 25cm resolutions constructed by the present invention respectively, and (4e)-(4h) are the enlarged detail views of (4a)-(4d) respectively;

[0049] Figure 5 is the terrain profile line analysis of the experimental area with OLA80cm and SPC40cm resolutions and the terrains with 40cm and 25cm resolutions generated by the method of the present invention in one embodiment. Among them, (5a) is a broken line graph of the terrain profile line 1 corresponding to the four terrains, (5b) is a schematic diagram of the image and the terrain profile line 1 (solid line) and the terrain profile line 2 (dashed line), (5c) is a broken line graph of the terrain profile line 2 corresponding to the four terrains, (5d) and (5e) are the partial enlarged views of (5a), and (5f) and (5g) are the partial enlarged views of (5c). DETAILED DESCRIPTION OF THE EMBODIMENTS

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

[0051] This embodiment provides a method for constructing a high-precision terrain at the pixel level of the asteroid surface based on an improved stereo photometric method, as Figure 1 shown, including the following steps:

[0052] S1, Data acquisition: Obtain multi-view image data of asteroids, where the image data is captured by a detector under different lighting conditions and viewpoints.

[0053] S2, Model construction: Establish a variational model for multi-view stereo photometry. The model reconstructs the terrain height field by minimizing the photometric error, and introduces multi-view stereo confidence constraints, non-convex estimators, and terrain edge protection mechanisms to construct the objective function.

[0054] Multi-view photometry takes an initial low-resolution terrain and multiple high-resolution images of corresponding regions as inputs, and optimizes the terrain based on the principle of photometry. The objective function for constructing the multi-view photometry model is as follows:

[0055]

[0056] where μ and λ are the coefficients of two regularization terms respectively, μ > 0, λ > 0; I k (z(x,y)) represents the pixel intensity at the corresponding position on the k-th image obtained by projecting the terrain z(x,y) through the camera model corresponding to the k-th image. Generally speaking, the position projected onto the image is not exactly an integer multiple of the pixel grid. Therefore, the acquisition of pixel intensity usually involves an interpolation process; T k represents the exposure of the k-th image, and there are many factors affecting T k . During the terrain optimization process, it is usually regarded as a constant; a(x,y) represents the albedo that depends on the terrain; R k (z(x,y)) represents the reflectivity calculated based on the reflectivity model using the observations and lighting conditions during the shooting of the k-th image; ▽ represents the second derivative solver.

[0057] In the above objective function, the first term is the photometric constraint, which is the core part of photometry. The method with only photometric constraint is often unstable and will fall into local minima. Therefore, consider the second constraint term, weighted by μ, which is usually called the smoothness constraint. It constrains the obtained terrain to be smooth, and the weight μ is used to control the smoothness. As we all know, SfS can construct high-resolution terrain details, but its drawback is the drift in the spatial scale. Therefore, add the third term to make the result of terrain optimization not deviate too far from the initial terrain z0(x,y).

[0058] Thus, the improved objective function of the present invention is expressed as:

[0059] F(z,a) = D image (z) + υ·D sfs (z,a) + α1·L edge (z) + α2·L dem (z) (2)

[0060] The objective function consists of two data terms and two regularization terms. Among them, D image is the multi-view stereo data term, representing the multi-view stereo confidence constraint. It is constructed by using the depth parameterization model proposed in "Dense depth map reconstruction: A minimization and regularization approach which preserves discontinuities" and utilizing the photometric consistency (i.e., brightness constancy) of the projected surface points; D sfs is the stereo photometric data term, which is constructed by using the moon-Lambert photometric model and utilizing the difference between the image rendered based on the photometric model and the real image; L egde represents the terrain edge regularization term, which is used to constrain the continuity of the terrain edge; L dem represents the terrain regularization term, which is used to constrain the correctness of the terrain geometric structure; ν, α1, and α2 are weights; z is the depth corresponding to the terrain in the camera space coordinate system, and a is the albedo.

[0061] 1. Multi-view stereo data term

[0062] Combining multi-view stereo and photometric constraints, the low-resolution terrain is optimized to a terrain with pixel-level accuracy. The multi-view stereo data term D image is expressed as:

[0063]

[0064] Among them, defines the brightness difference between the corresponding points in the reference image I0(x,y) and other views I k (x,y), defines the gradient difference between the corresponding points in the reference image I0(x,y) and other views I k (x,y); (x,y) represents the pixel position in the image, k represents the view index, n represents the total number of views, Ω0 represents the reference image domain, J(·) represents the Jacobian matrix, γ represents the weighting factor, which is used to balance the contributions of the two terms to D image , represents the L2 norm, represents the Frobenius norm, Γ D (·) represents the Geman-McClure non-convex estimator.

[0065] This embodiment introduces a multi-view stereo confidence constraint, aiming to reduce the participation of images that are too bright or too dark in the local neighborhood centered at the current point in the terrain construction of the current point. When an individual image differs significantly from the average image intensity level of the image set, it indicates that the image intensity of the current point in this image is not only determined by itself and the corresponding terrain, but is more likely to be affected by other factors. According to the photometric principle, such areas are not conducive to the accurate construction of the terrain of the current point. Data item D image is defined as in Equation (3). In addition to directly constraining the image intensity of the current point, its neighborhood constraint should also be considered to ensure the constancy of the gradient of the current point. Therefore, the Jacobian coefficient J is used to constrain the gradient change. In Equation (3), the Geman-McClure non-convex estimator is introduced, and its definition is:

[0066]

[0067] where s represents the independent variable of the function input, and ε is defined as a positive constant close to 0 to ensure that the function is differentiable at s = 0.

[0068] 2. Stereo photometric data item

[0069] The stereo photometric data item D sfs is expressed as:

[0070]

[0071] where k represents the view index, I represents the image, Ω0 represents the reference image domain, a(x, y) represents the albedo corresponding to the pixel at (x, y) in the image, z(x, y) represents the optimized high-resolution terrain, represents the L2 norm, R k represents the lunar-Lambert photometric model, T k represents the parameter related to the image exposure, and Ψ(·) is the Cauchy non-convex estimator.

[0072] For Equation (1), it uses the least squares method to constrain the consistency between the reflectivity and the intensity of the corresponding pixels in the image. In addition to the least squares method, the sparsity assumption is also a method based on residual minimization. The residuals defined by the least squares method are all "small", while most of the residuals in the sparsity assumption are "zero". The former is usually robust to "noise", while the latter is robust to "outliers". Here, "noise" and "outliers" are collectively referred to as "outliers". Then, when constraining the photometric using multiple images, due to image quality, registration quality, etc., such "outliers" obviously exist. For example, due to the shadow effect, the pixel intensity in the corresponding area of the image is unreliable; the noise in the terrain causes incorrect estimation of the normal vector, which in turn leads to unreliable calculation of the reflectivity, etc. To solve the above problems, this embodiment proposes to use an estimator that is robust to outliers, which combines the robustness of the least squares method to noise and the robustness of sparsity to discontinuities. As Figure 2 shown, it is a schematic diagram of the integration principle of the least squares and sparse representation. The robustness of the least squares method to noise comes from the quadratic behavior near 0, which ensures that "small" residuals are considered "good" estimates. However, this quadratic method will have such problems at ±∞, that is, it will generate "large" residuals, which will be over-penalized. The estimator based on the sparsity assumption has the opposite characteristics, that is, it treats high residuals (discontinuities) exactly the same as low residuals, ensuring that discontinuities are not over-penalized, while low residuals (noise) will be over-penalized. Therefore, a good estimator is quadratic around 0, but sparse linear around ±∞. Obviously, only non-convex estimators have these two characteristics. Denote the non-convex estimator as Ψ.

[0073] Classical non-convex estimators Ψ include the Geman-McClure estimator, the Welsh estimator, and the Cauchy estimator, etc., which are respectively defined as:

[0074]

[0075] ψ w (x) = λ 2 (1 - exp(-x 2 / λ 2 )) (7)

[0076] ψ C (x) = λ 2 (1 + log(1 + x 2 / λ 2 )) (8)

[0077] Among them, λ is an empirical parameter. This embodiment uses the Cauchy estimator to make the photometric term robust to outliers, then there is,

[0078] ψ(s) = λ 2(1 + log(1 + s 2 / λ 2 )) (9)

[0079] where s represents the independent variable of the function input, and λ is an empirical parameter.

[0080] 3. Topographic edge regularization term

[0081] The topographic edge regularization term L egde is expressed as:

[0082]

[0083] where z(x, y) represents the optimized high-resolution terrain, Ω0 represents the reference image domain, represents the Frobenius norm, Ψ(·) is a Cauchy non-convex estimator, and Hess(·) represents the smoothness constraint on the topographic edge, defined as:

[0084]

[0085] where z xx represents the second-order derivative of the terrain z in the x direction, z yy represents the second-order derivative of the terrain z in the y direction, z xy represents the gradient of the terrain z in the x direction, and z x represents the derivative of the gradient of z in the y direction.

[0086] In Equation (1), isotropic second-order regularization is used for processing the topographic edge, but it usually has the situation of over-smoothing the edge. Especially for terrains with drastic changes such as boulders, such a regularization process often erases the topographic details at the edge. According to the discussion of the non-convex estimator in the above 3, a Cauchy non-convex estimator is introduced in this embodiment to achieve edge protection.

[0087] 4. Topographic regularization term

[0088] The topographic regularization term L dem is defined as:

[0089]

[0090] where z0(x, y) represents the initial low-resolution terrain, z(x, y) represents the optimized high-resolution terrain, Ω0 represents the reference image domain, represents the L2 norm.

[0091] S3. Model Solving: The variable update problem of the multi-view stereo photometric variational model is transformed into an independent non-linear optimization problem. The objective function is decomposed into multiple sub-problems through the Lagrangian operator, and the quasi-Newton method is used for efficient solution. During the iteration process, the terrain resolution is gradually optimized, and finally high-precision terrain data at the pixel level is output.

[0092] Obviously, the variational problem in Equation (2) is not only non-convex, but also the expression of the gradient ▽z of the terrain z is non-linear. In this embodiment, instead of first optimizing the normal and then fitting an integrable surface to solve the terrain in a two-step method, z will be directly solved. Therefore, the global and non-linear problem is transformed into a continuous global-linear and non-linear-local problem. By introducing an auxiliary vector field to directly solve the terrain z, the solution processes of the terrain and albedo are expressed as an equivalent constrained optimization problem:

[0093]

[0094] where θ is the Lagrangian operator, represents the gradient of the terrain z, a represents the albedo, represents the terrain, albedo, and Lagrangian operator obtained by optimization and solution.

[0095] 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) . Assume that the initial value of the reflectivity is 1, that is, it is uniformly white. At the iteration step (t), each variable is updated according to the following steps:

[0096]

[0097] where γ, β, and κ all represent weight parameters, defined as empirical constants. The update of the Lagrangian operator θ is a series of independent non-linear optimization problems, which are solved using the implementation of the L-BFGS method; the update of the dual operator u is a non-linear optimization sub-problem, which is independently solved using the L-BFGS method in each pixel; the update of the terrain z uses the conjugate gradient of the normal equation to solve this problem; the solution of the albedo a is a linear least squares problem, which is solved using the gradient descent method.

[0098] This embodiment conducts experimental verification on the above method.

[0099] The experimental data are the images of the asteroid Bennu taken during the OSIRIS-Rex exploration mission and the high-precision topographic data of the global or local areas released. Among them, the image data include PolyCam images and navigation camera images. These images have different image resolutions and correspond to different lighting and observation conditions. It should be noted that according to the principle and experience of constructing topography by photometry, the images used to construct high-precision topography need to meet specific conditions. In this embodiment, the lighting and observation conditions that the images should meet are: (i) the incident angle is preferably about 45° (acceptable range is 30 - 50°, limit is 0 - 60°); (ii) the exit angle is preferably 0° (acceptable range is 0 - 20°, limit is 0 - 60°). The topographic data are the local topographic data released by NASA. Both the image data and the topographic data can be downloaded from the PDS Earth Science Node. In addition, the ephemeris data used in the proposal can also be downloaded from the PDS Earth Science Node. These images and topographic data will be used as the input data of the proposed method and the data for accuracy comparison. The experimental area is the actual landing area "Nightingale" in the Bennu exploration mission, and the longitude and latitude coordinates of the center position are (42°, 56°). The characteristics of this area are that the center position is flat, the surrounding terrain is rugged, and it is covered with boulders as Figure 3 shown.

[0100] To evaluate the applicability of the present invention to different topographic characteristics and the characteristics of the method, for the experimental area, topographies and albedos with resolutions of 40 cm and 25 cm were respectively constructed based on the proposed method. Among them, the input data are the global topographic data with resolutions of 80 cm and SPC40 cm globally produced by NASA using OLA (OSIRIS-REx Laser Altimeter) data. The specific input and output data are shown in Table 1.

[0101] Table 1 General situation of input and output data of the method of the present invention

[0102]

[0103] The comparison of the topographic construction results in this embodiment is divided into two parts: (i) comparison with the topographic data released by NASA; (ii) comparison with other typical topographic construction methods. The quantitative topographic accuracy evaluation indexes used in the experiment are: Root Mean Square Error (RMSE), Max Absolute Error (Max-AE), and Mean Absolute Error (Mean-AE).

[0104] As Figure 3As shown in the image, the characteristics of this area are that there is a flat area in the middle with a diameter of about 8m, which is the actual landing point of OSIRIS-Rex. However, "giant boulders" are distributed around it. Compared with the flat area, the elevation changes of these "giant boulders" are very obvious. From these images, it can be found that some stones show obvious albedo contrast compared with their surrounding environment, while some stones have little difference in albedo from their surrounding environment. Accurately solving the terrain and albedo is very challenging.

[0105] As Figure 4 shown in (4a) and (4c), the terrain with OLA 80cm resolution and SPC 40cm resolution of the selected area respectively are the terrain data input for constructing the navigation landmarks with 40cm resolution and 25cm resolution in the present invention. As Figure 4 shown in (4b) and (4d), the terrains with 40cm resolution and 25cm resolution generated by the present invention. Figure 4 The coordinates of the terrains shown in (4a)-(4d) are in the fixed coordinate system of the asteroid Bennu. To more clearly show the details of each terrain, as Figure 4 shown in (4e)-(4h) are respectively Figure 4 the enlarged views of the terrains in (4a)-(4d). As shown by the red circles in the figure, whether it is the terrain with 40cm resolution or 25cm resolution, compared with the input terrain data, the terrain constructed by the present invention is smoother. Especially in the areas with large terrain undulations. As in the red circle, it is the "giant boulder" area, and the OLA data shows obvious "cliff-like" terrain changes. Although the terrain constructed by the present invention also has the "cliff-like" phenomenon, it is relatively smoother. As in the purple circle, it is also a stone area, but compared with the "giant boulder" area, the terrain change is more gentle. In this case, the OLA data still has the "cliff-like" terrain phenomenon, while the present invention processes the corresponding terrain more continuously. On the one hand, as the resolution increases, the terrain change should indeed be more continuous. Generally speaking, compared with the SPC 40cm resolution terrain and the OLA 80cm resolution terrain, the terrain with 40cm resolution is smoother than the terrain with 80cm resolution. On the other hand, the terrain change should also be continuous in reality, which shows that the present invention processes the terrain edge more reasonably.

[0106] Comparing the terrain with 40cm resolution constructed by the present invention with the SPC 40cm resolution terrain, visually, the restoration of the terrain details of both is relatively similar, indicating the correctness of the terrain construction method of the present invention. Compared with the OLA 80cm resolution terrain, many terrain details are constructed on the 40cm resolution terrain. The terrain with 25cm resolution constructed by the present invention based on the OLA 40cm resolution terrain is as Figure 4As shown in (4h), compared with the 40cm resolution terrain, the 25cm resolution terrain has more terrain details, such as Figure 4 As shown in the green circle in the middle. More terrain details mean that the rendered image has more texture details, which is obviously conducive to the precise matching of subsequent navigation landmarks with the corresponding resolution image, and thus conducive to accurate navigation positioning.

[0107] In order to further quantitatively evaluate the accuracy of the terrain constructed by the proposed method, two profile lines were drawn according to the texture of the image, and the change characteristics of the above four terrains were quantitatively analyzed through the profile lines. For the convenience of representation, the method proposed in this invention is named "Pro-SFS". Figure 5 (5a) and (5c) of Figure 5 The yellow solid and dashed topographic profiles in (5b) are shown in Figure 5. For areas where the profiles have typical changes, they are further enlarged, such as Figure 5 The first profile line is selected in the area where the terrain changes are relatively flat. It can be seen that the terrain changes of the four terrains in this area are very similar. Compared with the 40cm resolution terrain, the 25cm resolution terrain of the present invention restores more detailed terrain changes. Two typical sub-areas are selected, respectively. Figure 5 As shown in the yellow and green boxes in (5d) and (5e), ​​the OLA80cm terrain is quite different from the other three terrains, and the OLA80cm terrain is 0.2m higher than the other three terrains on average. The second profile line passes through a "boulder" area. It can be seen that the terrain changes of the four are generally very consistent, and the areas with the largest differences occur where the terrain changes dramatically, corresponding to the image, which is the edge of the stone. Figure 5 As shown in (5f) and (5g), according to the change of the profile line, the topography change trends of 40 cm and 25 cm of the present invention are consistent, while the topography change trends of OLA and SPC are consistent. Figure 5 In (5g), the OLA80cm terrain is about 0.2m higher than the SPC40cm terrain, and the SPC40cm terrain is about 0.5m higher than the terrain of the present invention. These phenomena illustrate that the processing of the edge by the present invention is to make it tend to be smooth. Taking SPC40cm as a reference standard and comparing the other three terrains with it, the quantitative evaluation results obtained are shown in Table 2. According to Table 2, it can be seen that the four terrains are very close, which illustrates the correctness of the present invention. The absolute root mean square error (RMSE) of the terrain profile is 0.3 meters, and the maximum error is about 1.2240 meters; the larger terrain difference originates from the edge area of ​​the "boulder".

[0108] Quantitative evaluation of the terrain with OLA 80 cm and SPC 40 cm resolutions in the experimental area and the topographic profile lines with 40 cm and 25 cm resolutions generated by the present invention

[0109]

[0110] The present invention proposes a method for constructing a high-precision terrain at the pixel level based on multi-view photometry under a variational framework. Compared with the traditional terrain construction method based on multi-view photometry, the present invention respectively introduces multi-view stereo confidence constraints to ensure the effectiveness of the current image intensity for constructing the terrain at the current point; introduces a non-convex estimator to enhance the robustness of photometric estimation; and introduces a convex estimator to enhance the terrain edge constraint to reduce the unevenness of the terrain edge. Experimental results show that the terrain with a 40 cm resolution constructed by the present invention is similar to the terrain with a 40 cm resolution constructed by SPC. In contrast, the edge of the terrain constructed by the present invention is smoother and has more terrain details.

[0111] 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 efforts. 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 shall fall within the protection scope determined by the claims.

Claims

1. An improved stereo photometric method for constructing high-precision terrain at the pixel level on the asteroid surface, characterized in that Including the following steps: Data acquisition: Acquire multi-view image data of asteroids, where the image data is obtained by a detector shooting under different lighting conditions and viewpoints; Model construction: Establish a variational model of multi-view stereo photometry. The model reconstructs the terrain height field by minimizing the photometric error, and introduces multi-view stereo confidence constraints, non-convex estimators, and terrain edge protection mechanisms to construct the objective function; Model solution: Convert the variable update problem of the multi-view stereo photometry variational model into an independent non-linear optimization problem. Decompose the objective function into multiple sub-problems through the Lagrangian operator, and use the quasi-Newton method for efficient solution; During the iteration process, gradually optimize the terrain resolution, and finally output pixel-level high-precision terrain data.

2. The high-precision terrain construction method for asteroid surface pixels based on improved stereophotometry according to claim 1, wherein The objective function is expressed as: F(z,a) = D image (z) + υ·D sfs (z,a) + α1·L edge (z) + α2·L dem (z) The objective function consists of two data terms and two regularization terms. Among them, D image is the multi-view stereo data term, representing the multi-view stereo confidence constraint. It is constructed by using a depth parameterized model and the photometric consistency of the projected surface points. D sfs is the stereo photometric data term, which is constructed by using the lunar-Lambert photometric model and the difference between the image rendered based on the photometric model and the real image. L egde represents the terrain edge regularization term, which is used to constrain the continuity of the terrain edge. L dem represents the terrain regularization term, which is used to constrain the correctness of the terrain geometric structure. ν, α1, and α2 are weights; z is the depth corresponding to the terrain in the camera space coordinate system, and a is the albedo.

3. A method for constructing a high-precision terrain at the pixel level on the surface of an asteroid based on an improved stereo photometry method according to claim 2, characterized in that, Combining multi-view stereo and photometric constraints to optimize low-resolution terrain into terrain with pixel-level accuracy, the multi-view stereo data item D image is expressed as: Among them, defines the luminance difference between the corresponding points of the reference image I0(x, y) and other views I k (x, y), defines the gradient difference between the corresponding points of the reference image I0(x, y) and other views I k (x, y); (x, y) represents the pixel position in the image, k represents the view index, n represents the total number of views, Ω0 represents the reference image domain, J(·) represents the Jacobian matrix, γ represents the weighting factor used to balance the contributions of the two terms to D image , represents the L2 norm, represents the Frobenius norm, Γ D (·) represents the Geman-McClure non-convex estimator.

4. A method for constructing a high-precision terrain at the pixel level on the asteroid surface based on an improved stereo photometry method according to claim 2, characterized in that The three-dimensional photometric data item D sfs is expressed as: where k represents the view index, I represents the image, Ω0 represents the reference image domain, a(x, y) represents the albedo corresponding to the pixel at (x, y) in the image, z(x, y) represents the optimized high-resolution terrain, represents the L2 norm, R k represents the lunar-Lambert photometric model, T k represents the parameter related to the image exposure, and Ψ(·) is the Cauchy non-convex estimator.

5. A method for constructing a high-precision terrain at the pixel level on the surface of an asteroid based on an improved stereo photometry method according to claim 2, characterized in that, The terrain edge regularization term L egde is expressed as: where \(z(x,y)\) represents the optimized high-resolution terrain, \(\Omega_0\) represents the reference image domain, represents the Frobenius norm, \(\Psi(\cdot)\) is the Cauchy non-convex estimator, and \(\text{Hess}(\cdot)\) represents the smoothness constraint on the terrain edges, defined as: where z xx represents the second derivative of the terrain z in the x - direction, z yy represents the second derivative of the terrain z in the y - direction, z xy represents the gradient of the terrain z in the x - direction, and z x represents the derivative of z in the y - direction.

6. The method for constructing a high-precision terrain at the pixel level on the asteroid surface based on the improved stereo photometry method according to claim 2, wherein The terrain regularization term L dem is defined as: Among them, z0(x, y) represents the initial low-resolution terrain, z(x, y) represents the optimized high-resolution terrain, and Ω0 represents the reference image domain. represents the L2 norm.

7. A method for constructing a high-precision terrain at the pixel level of an asteroid surface based on an improved stereo photometry method according to claim 3, characterized in that, The Geman-McClure non-convex estimator is defined as: Where s represents the independent variable of the function input, and ε is defined as a positive constant close to 0 to ensure that the function is differentiable at s = 0.

8. A method for constructing a high-precision terrain at the pixel level on the asteroid surface based on an improved stereo photometry method according to claim 4 or 5, characterized in that, The Cauchy non-convex estimator is defined as: ψ(s) = λ 2 (1 + log(1 + s 2 / λ 2 )) Where s represents the independent variable of the function input, and λ is an empirical parameter.

9. A method for constructing a high-precision terrain at the pixel level on the asteroid surface based on an improved stereo photometry method according to claim 2, characterized in that, In the solution of the model, by introducing an auxiliary vector field θ:= to directly solve for the terrain z, the solution processes of the terrain and albedo are expressed as an equivalent constrained optimization problem: where θ is the Lagrange operator, represents the gradient of the terrain z, and a represents the albedo, represents the terrain, albedo, and Lagrange operator obtained by optimized solution.

10. A high-precision terrain construction method for asteroid surface pixels based on improved stereo photometry according to claim 9, characterized in that In the model solution, it is assumed 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) . It is assumed that the initial value of the reflectivity is 1, that is, it is uniform white. At the iteration step (t), the variables are updated according to the following steps: where k represents the view index, I0(x,y) is the reference image, and I k (x,y) is other views, represents the L2 norm, represents the Frobenius norm, Γ D (·) represents the Geman-McClure non-convex estimator, R k represents the moon-Lambert photometric model, T k represents the parameters related to image exposure, Ψ(·) is the Cauchy non-convex estimator, γ, β, and κ all represent weight parameters, defined as empirical constants, Hess(·) represents the smoothness constraint on the terrain edge, z0 represents the initial low-resolution terrain; the update of the Lagrangian operator θ is a series of independent non-linear optimization problems, solved using the implementation of the L-BFGS method; the update of the dual operator u is a non-linear optimization sub-problem, independently solved in each pixel using the L-BFGS method; the update of the terrain z uses the conjugate gradient of the normal equation to solve this problem; the solution of the albedo a is a linear least squares problem, solved using the gradient descent method.

Images