Moon antarctic terrain reconstruction method based on prior terrain information guiding stereo matching

Abstract

The invention relates to a lunar Antarctic terrain reconstruction method based on prior terrain information guiding stereo matching. The method comprises the following steps: adaptively selecting a seed point in a shadow area through image intensity analysis; identifying homogeneous pixels around each seed point according to distance and intensity differences, and obtaining a prior parallax of an introduction point based on prior terrain information; the corresponding cost is simulated according to the prior parallax to replace the original fuzzy cost, so that the cost volume is optimized, and the matching difficulty region is guided to obtain the more accurate cost; semi-global matching cost aggregation is carried out, and a disparity map is obtained by carrying out disparity selection by adopting a winner communication strategy; and performing triangulation processing on the disparity map to obtain a point cloud, and performing point cloud interpolation to generate a terrain. Compared with the prior art, by introducing low-resolution terrain prior information, accurate parallax estimation at the shadow of the image can be realized, so that the accuracy and the integrity of terrain reconstruction of the south moon are improved.

Description

Technical Field

[0001] The invention relates to the technical field of lunar photogrammetry, and in particular to a lunar south pole terrain reconstruction method based on prior terrain information guided stereo matching. Background Art

[0002] With the continuous advancement of deep space exploration technology, the lunar South Pole has become a hot spot for exploration due to its unique scientific value and potential resources. There may be a large amount of water ice in the lunar South Pole region, which is crucial for the construction of future lunar bases and deep space exploration missions. Therefore, the acquisition of high-precision lunar South Pole terrain data is of great significance for planning landing sites, evaluating resource distribution, and formulating scientific exploration plans.

[0003] High-precision terrain acquisition relies on photogrammetry, which can build a 3D model from images. Dense matching is a key step in photogrammetry, which involves finding points of the same name between two or more images, thereby obtaining a dense 3D point cloud through forward intersection. The dense matching process places strict requirements on the accuracy and robustness of the matching algorithm, especially when facing the complex and changeable terrain conditions of the lunar South Pole. The lunar South Pole has a rugged terrain with large elevation differences and is dotted with large and small impact craters. Due to the low solar altitude angle, the South Pole region is often covered by long-term shadows, resulting in almost no effective information in some areas in the image. This lack of information further exacerbates the difficulty of image registration and 3D reconstruction, especially when it comes to large-scale terrain mapping. Therefore, it is urgent to develop a dense matching algorithm that can be robust in complex lighting and multi-shadow environments. This is not only of great significance for terrain mapping in the lunar South Pole, but also provides important technical support for scientific exploration missions such as water ice distribution, resource exploration and environmental monitoring. Summary of the invention

[0004] The purpose of the present invention is to overcome the defects of the above-mentioned prior art and provide a method for reconstructing the lunar South Pole terrain based on stereo matching guided by prior terrain information, which can achieve accurate parallax estimation in image shadows, thereby improving the accuracy and completeness of the lunar South Pole terrain reconstruction.

[0005] The object of the present invention can be achieved by the following technical solution: A method for reconstructing lunar South Pole terrain based on stereo matching guided by prior terrain information, comprising the following steps:

[0006] S1, adaptively select seed points in the shadow area through image intensity analysis;

[0007] S2, identifying homogeneous pixels around each seed point according to the distance and intensity difference, and obtaining the prior disparity of the introduced point based on the prior terrain information;

[0008] S3, simulate the corresponding cost based on the prior disparity to replace the original blurred cost, so as to optimize the cost volume and guide the matching difficult area to obtain a more accurate cost;

[0009] S4, performing semi-global matching cost aggregation, adopting a winner-takes-all strategy to select disparity and obtain a disparity map;

[0010] S5. Triangulate the disparity map to obtain a point cloud, and perform point cloud interpolation to generate terrain.

[0011] Furthermore, the step S1 includes the following steps:

[0012] S11, obtain high-resolution LROC-NAC images and low-resolution terrain data, and perform image preprocessing to obtain updated external and epipolar images;

[0013] S12. Perform image intensity analysis on the epipolar image and select seed points.

[0014] Further, the step S11 includes the following process: obtaining an initial external square from a high-resolution LROC-NAC image, and performing feature matching on different high-resolution LROC-NAC images to obtain connection points;

[0015] Based on the connection points and low-resolution terrain data, the initial external square is adjusted to obtain an updated external square;

[0016] Generates an epipolar image consisting of the left and right images based on the tie points.

[0017] Furthermore, the step S12 includes the following process:

[0018] S121, performing intensity histogram calculation on the epipolar image to obtain a shadow determination threshold;

[0019] S122, binarizing the image according to the shadow determination threshold to obtain a preliminary shadow mask;

[0020] S123, performing Gaussian blur and morphological opening and closing operations on the preliminary shadow mask to increase the continuity of the shadow area;

[0021] S124, binarizing the mask and removing the shadow area that is too small to obtain a final shadow mask;

[0022] S125. Perform uniform sampling on the final shadow mask to obtain seed points.

[0023] Furthermore, the shadow determination threshold is specifically the intensity value of the first trough in the epipolar image intensity histogram.

[0024] Furthermore, the step S2 comprises the following steps:

[0025] S21, based on the seed point, obtaining the object point in the three-dimensional space and back-projecting it onto the epipolar image to obtain the prior disparity of a single pixel;

[0026] S22, calculating the distance and grayscale value between different pixels, and determining homogeneous pixels around the seed point;

[0027] S23. Obtain the prior disparity of homogeneous pixels according to the prior disparity of a single pixel.

[0028] Furthermore, the specific process of step S22 is as follows:

[0029] With the seed point as the center, define an initial window of size n*n, and calculate the weights of all pixels in the window. The formula is as follows:

[0030]

[0031] Among them, (i, j) and (xL, yL) represent the pixel coordinates of the seed point and the current evaluation pixel coordinates, I(i, j) and I(xL, yL) represent their corresponding intensity values, σ xy is the influence weight of spatial distance, σ l is the impact weight of intensity difference;

[0032] Combined with the given threshold γ, the homogeneous pixels around the seed point are determined by the following formula:

[0033] H(x L ,y L )={(i,j)∈N(x L ,y L )|W (i,j) ≥γ}

[0034] Among them, H(x L ,y L ) represents a set of similar pixels around the seed point, N(x L ,y L ) represents the pixel (x L ,y L ) is the neighborhood pixel with a window size of n*n as the center.

[0035] Furthermore, the specific process of step S3 is as follows:

[0036] Given the disparity of the introduced seed point, the simulation cost is obtained using the generalized Euclidean distance cost function according to the distance between each disparity and the seed point disparity within the disparity search range. The specific formula is as follows:

[0037]

[0038] Among them, dpx and dpy are the disparities of the current seed point, dx and dy are the candidate disparities within the disparity search range, parameter A is used to control the overall cost value to prevent the value from being too large, and parameter α is used to control the smoothness of the cost function;

[0039] The surface elevation error value at each seed point is directly used as the α value to adaptively obtain the cost curve of each seed point according to its reliability, and the cost of the introduced seed point replaces the original fuzzy cost value.

[0040] Furthermore, the step S4 comprises the following steps:

[0041] S41, obtaining aggregated costs by semi-global matching cost aggregation; adopting a dynamic programming strategy to aggregate the costs of each pixel from eight directions, namely, from the upper left direction to the lower right direction, to obtain an aggregated cost space;

[0042] S42, using a winner-takes-all strategy to obtain parallax;

[0043] S43, performing left-right consistency detection on the disparity map to eliminate erroneous disparity in the occluded area;

[0044] S44: Use median filtering and texture filtering to smooth the disparity map.

[0045] Furthermore, the step S5 specifically utilizes the updated exterior orientation element information of the image, calculates the same-name points on the original left and right images according to the disparity map, and performs forward intersection based on the midpoint method to obtain the object point cloud, and then uses inverse distance weighted interpolation to reconstruct the terrain for the object point cloud.

[0046] Compared with the prior art, the present invention has the following advantages:

[0047] The present invention firstly selects seed points in the shadow area adaptively by analyzing the image intensity; then identifies the homogeneous pixels around each seed point according to the distance and intensity difference; then obtains the prior disparity of the introduced point based on the prior terrain information, and simulates the corresponding cost to replace the original blurred cost according to the prior disparity, optimizes the cost volume, performs semi-global matching cost aggregation, adopts the winner-takes-all strategy to select the disparity to obtain the disparity map; finally, triangulates the obtained disparity map to obtain the point cloud, and performs point cloud interpolation to generate the terrain. Thus, by introducing the prior terrain information to guide the stereo matching process, the matching robustness can be enhanced, especially in the shadowy areas of the lunar South Pole; in addition, the proposed adaptive seed point selection strategy can adapt to images with different characteristics, accurately select seed points in the matching difficult areas, and thus more effectively improve the reconstruction accuracy of these areas.

[0048] The present invention guides stereo matching based on prior terrain information, obtains the prior parallax of seed points by introducing low-resolution terrain back projection, and simulates the corresponding cost based on the prior parallax to replace the original blurred cost, thereby combining the prior cost with the original semi-global matching cost, effectively improving the matching effect in the multi-shadow area of ​​the lunar South Pole image.

[0049] The present invention proposes a seed point selection algorithm based on image intensity information, and sets the shadow judgment threshold according to the image intensity histogram, which can effectively obtain the shadow mask corresponding to the image, thereby more accurately selecting the seed points located in the shadow area, and effectively introducing prior information in the difficult matching area.

[0050] The present invention obtains the seed point cost based on prior terrain information, and controls the fluctuation trend of the prior cost by combining the reliability of low-resolution terrain elevation information, thereby enhancing the cost robustness and reliably guiding the propagation of effective information. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] Figure 1 It is a schematic diagram of the method flow of the present invention;

[0052] Figure 2 It is a schematic diagram of the application process of the embodiment;

[0053] Figure 3 It is a local area display of the LROC NAC image in the embodiment;

[0054] Figure 4 is the seed point distribution selected in the embodiment;

[0055] Figure 5 The cost obtained by simulating the generalized Euclidean distance cost function in the embodiment;

[0056] Figure 6 A disparity map generated in an embodiment;

[0057] Figure 7a and 7b These are the terrains generated by the method of the present invention and the terrain maps generated by the semi-global matching method in the embodiments. DETAILED DESCRIPTION

[0058] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0059] Example

[0060] like Figure 1 As shown, a method for reconstructing lunar South Pole terrain based on stereo matching guided by prior terrain information includes the following steps:

[0061] S1, adaptively select seed points in the shadow area through image intensity analysis;

[0062] S2, identifying homogeneous pixels around each seed point according to the distance and intensity difference, and obtaining the prior disparity of the introduced point based on the prior terrain information;

[0063] S3, simulate the corresponding cost based on the prior disparity to replace the original blurred cost, so as to optimize the cost volume and guide the matching difficult area to obtain a more accurate cost;

[0064] S4, performing semi-global matching cost aggregation, adopting a winner-takes-all strategy to select disparity and obtain a disparity map;

[0065] S5. Triangulate the disparity map to obtain a point cloud, and perform point cloud interpolation to generate terrain.

[0066] This embodiment applies the above solution, such as Figure 2 As shown, the main contents are:

[0067] First, adaptively select seed points in shadow areas through image intensity analysis;

[0068] 101. Calculate the intensity histogram of the image and obtain the intensity value of the first trough as the shadow determination threshold;

[0069] 102. Binarize the image according to the shadow determination threshold to obtain a preliminary shadow mask;

[0070] 103. Perform Gaussian blur and morphological opening and closing operations on the preliminary shadow mask to increase the continuity of the shadow area;

[0071] 104. Binarize the mask and remove the shadow area that is too small to obtain a final shadow mask.

[0072] 105. Perform uniform sampling on the final shadow mask to obtain seed points;

[0073] Figure 3 This is a local area display of the Lunar Reconnaissance Orbiter Camera-Narrow Angle Camera (LROC-NAC) image in this embodiment. The image size is 1029*1036 pixels. The image intensity is analyzed by histogram to obtain a preliminary shadow mask, and the shadow mask is subjected to morphological operations to obtain the final shadow mask. Uniform seed points are selected on the shadow mask at intervals of 5 pixels. Figure 4 is the seed point distribution selected in the embodiment, from Figure 4 It can be seen that the selected seed points are basically located in the shadow area of ​​the image and are evenly distributed in the shadow area.

[0074] Second, identify homogeneous pixels around each seed point based on distance and intensity differences;

[0075] After obtaining the seed point, define an initial window of size n*n with the seed point as the center, and calculate the weights of all pixels in the window. The formula is as follows:

[0076]

[0077] Among them, (i, j) and (xL, yL) represent the pixel coordinates of the seed point and the current evaluation pixel coordinates, respectively, and I(i, j) and I(xL, yL) represent their corresponding intensity values. xy is the influence weight of spatial distance, σ L is the influence weight of the intensity difference. Given a threshold γ, the homogeneous pixels around the seed point can be determined by the following formula.

[0078] H(x L ,y L )={(i,j)∈N(x L ,y L )|W (i,j) ≥γ}

[0079] Among them, H(x L ,y L ) represents a set of similar pixels around the seed point, N(x L ,y L ) represents the pixel (x L ,y L ) is the neighborhood pixel with a window size of n*n as the center.

[0080] This step is mainly used to expand the influence range of the seed point. The selected homogeneous pixels have similar disparity with the seed point, so they can also be included as prior information for cost modification.

[0081] 3. Obtain the prior disparity of the introduction point based on the prior terrain information, and simulate the corresponding cost based on the prior disparity to replace the original blurred cost, so as to optimize the cost volume and guide the matching difficult area to obtain a more accurate cost.

[0082] According to the prior terrain information, the back-projection point on the left and right images can be obtained, thereby obtaining the disparity value of the point. Given the disparity of the introduced seed point, the simulation cost can be obtained by using the generalized Euclidean distance cost function based on the distance between each disparity and the disparity of the seed point within the disparity search range. The specific formula is as follows:

[0083]

[0084] Among them, dpx and dpy are the disparities of the current seed point, dx and dy are the candidate disparities within the disparity search range, parameter A is used to control the overall cost value to prevent the value from being too large, and parameter α is used to control the smoothness of the cost function. Considering that the disparity of the seed point is obtained by low-resolution terrain back projection, there is a certain uncertainty, and this uncertainty varies due to factors such as sampling position and elevation accuracy. In the process of prior cost calculation, the surface elevation error of each grid from the low-resolution prior terrain model is combined to determine the weight of each prior point. Specifically, the surface elevation error value at each seed point is directly used as the α value. In this way, the cost curve of each seed point is adaptively obtained according to its reliability. The cost of the introduced seed point is used to replace the original fuzzy cost value, thereby realizing the introduction of robust cost value to guide the optimization of the entire cost volume. Figure 5 is the cost obtained by simulating the generalized Euclidean distance cost function in the embodiment. In this embodiment, the low-resolution prior terrain uses the laser altimetry LOLADEM with a resolution of 20 m / pixel. The seed point parallax is dpx=4, dpy=4, A=1, α=0.05. It can be seen from the cost curve that the cost at the seed point is the smallest, while the cost values ​​of other parallax positions gradually increase with the increase of the parallax distance from the seed point.

[0085] Fourth, carry out semi-global matching cost aggregation, and adopt the winner-takes-all strategy to perform disparity selection to obtain the disparity map; specifically, first obtain the aggregated cost through semi-global matching cost aggregation; adopt the dynamic programming strategy to aggregate the cost of each pixel from eight directions (upper left to lower right) to obtain the aggregated cost space.

[0086] After that, the winner-takes-all strategy is adopted to obtain the disparity, and then the left-right consistency of the disparity map is checked, the erroneous disparity in the occluded area is eliminated, and finally the disparity map is smoothed using median filtering and texture filtering.

[0087] Figure 6 It is a disparity map generated in the embodiment. It can be seen from the disparity map result that the disparity result is relatively smooth, and a smooth and accurate disparity can also be obtained in the corresponding image shadow area.

[0088] Fifth, using the image's exterior orientation element information, the same-name points on the original left and right images are calculated based on the obtained disparity map, and forward intersection based on the midpoint method is performed to obtain the object point cloud. After obtaining the object point cloud, the terrain is obtained by using inverse distance weighted interpolation.

[0089] Figure 7a is the terrain of the area. It can be seen from the terrain results that the method proposed in the present invention can effectively obtain accurate disparity estimation for the shadow area, thereby generating high-precision terrain. Figure 7a) and the terrain generated by the original semi-global matching method ( Figure 7b ) were quantitatively analyzed and compared with the public LOLA DEM. The results showed that the average absolute error of the terrain obtained by the method of the present invention was 1.54 meters, and the root mean square error was 1.82 meters, while the average absolute error of the terrain obtained by the semi-global matching method was 1.75 meters, and the root mean square error was 2.27 meters. This proves that the method of the present invention effectively improves the accuracy of terrain reconstruction.

Claims

1. A method for reconstructing lunar South Pole terrain based on stereo matching guided by prior terrain information, characterized in that: The following steps are involved: S1, adaptively select seed points in the shadow area through image intensity analysis; S2, identifying homogeneous pixels around each seed point according to the distance and intensity difference, and obtaining the prior disparity of the introduced point based on the prior terrain information; S3, simulate the corresponding cost based on the prior disparity to replace the original blurred cost, so as to optimize the cost volume and guide the matching difficult area to obtain a more accurate cost; S4, performing semi-global matching cost aggregation, adopting a winner-takes-all strategy to select disparity and obtain a disparity map; S5. Triangulate the disparity map to obtain a point cloud, and perform point cloud interpolation to generate terrain.

2. The method for reconstructing lunar South Pole terrain based on stereo matching guided by prior terrain information according to claim 1, characterized in that: The step S1 comprises the following steps: S11, obtain high-resolution LROC-NAC images and low-resolution terrain data, and perform image preprocessing to obtain updated external and epipolar images; S12. Perform image intensity analysis on the epipolar image and select seed points.

3. The method for reconstructing lunar South Pole terrain based on stereo matching guided by prior terrain information according to claim 2, characterized in that: The step S11 includes the following processes: obtaining an initial external square from a high-resolution LROC-NAC image, and performing feature matching on different high-resolution LROC-NAC images to obtain connection points; Based on the connection points and low-resolution terrain data, the initial external square is adjusted to obtain an updated external square; Generates an epipolar image consisting of the left and right images based on the tie points.

4. The method for reconstructing lunar South Pole terrain based on stereo matching guided by prior terrain information according to claim 2, characterized in that: The step S12 includes the following process: S121, performing intensity histogram calculation on the epipolar image to obtain a shadow determination threshold; S122, binarizing the image according to the shadow determination threshold to obtain a preliminary shadow mask; S123, performing Gaussian blur and morphological opening and closing operations on the preliminary shadow mask to increase the continuity of the shadow area; S124, binarizing the mask and removing the shadow area that is too small to obtain a final shadow mask; S125. Perform uniform sampling on the final shadow mask to obtain seed points.

5. The method for reconstructing lunar South Pole terrain based on stereo matching guided by prior terrain information according to claim 4, characterized in that: The shadow determination threshold is specifically the intensity value of the first trough in the epipolar image intensity histogram.

6. The method for reconstructing lunar South Pole terrain based on stereo matching guided by prior terrain information according to claim 2, characterized in that: The step S2 comprises the following steps: S21, based on the seed point, obtaining the object point in the three-dimensional space and back-projecting it onto the epipolar image to obtain the prior disparity of a single pixel; S22, calculating the distance and grayscale value between different pixels, and determining homogeneous pixels around the seed point; S23. Obtain the prior disparity of homogeneous pixels according to the prior disparity of a single pixel.

7. The method for reconstructing lunar South Pole terrain based on stereo matching guided by prior terrain information according to claim 6, characterized in that: The specific process of step S22 is as follows: With the seed point as the center, define an initial window of size n*n, and calculate the weights of all pixels in the window. The formula is as follows: Among them, (i, j) and (xL, yL) represent the pixel coordinates of the seed point and the current evaluation pixel coordinates, I(i, j) and I(xL, yL) represent their corresponding intensity values, σ xy is the influence weight of spatial distance, σ L is the impact weight of intensity difference; Combined with the given threshold γ, the homogeneous pixels around the seed point are determined by the following formula: H(x L ,y L )={(i,j)∈N(x L ,y L )|W (i,j) ≥γ} Among them, H(x L ,y L ) represents a set of similar pixels around the seed point, N(x L ,y L ) represents the pixel (x L ,y L ) is the neighborhood pixel with a window size of n*n as the center.

8. The method for reconstructing lunar South Pole terrain based on stereo matching guided by prior terrain information according to claim 3, characterized in that: The specific process of step S3 is: Given the disparity of the introduced seed point, the simulation cost is obtained using the generalized Euclidean distance cost function according to the distance between each disparity and the seed point disparity within the disparity search range. The specific formula is as follows: Among them, dpx and dpy are the disparities of the current seed point, dx and dy are the candidate disparities within the disparity search range, parameter A is used to control the overall cost value to prevent the value from being too large, and parameter α is used to control the smoothness of the cost function; The surface elevation error value at each seed point is directly used as the α value to adaptively obtain the cost curve of each seed point according to its reliability, and the cost of the introduced seed point replaces the original fuzzy cost value.

9. The method for reconstructing lunar South Pole terrain based on stereo matching guided by prior terrain information according to claim 8, characterized in that: The step S4 comprises the following steps: S41, obtaining aggregated costs by semi-global matching cost aggregation; adopting a dynamic programming strategy to aggregate the costs of each pixel from eight directions, namely, from the upper left direction to the lower right direction, to obtain an aggregated cost space; S42, using a winner-takes-all strategy to obtain parallax; S43, performing left-right consistency detection on the disparity map to eliminate erroneous disparity in the occluded area; S44: Use median filtering and texture filtering to smooth the disparity map.

10. The method for reconstructing lunar South Pole terrain based on stereo matching guided by prior terrain information according to claim 3, characterized in that: The step S5 specifically utilizes the updated exterior orientation element information of the image, calculates the same-name points on the original left and right images according to the disparity map, and performs forward intersection based on the midpoint method to obtain the object point cloud, and then uses inverse distance weighted interpolation to reconstruct the terrain for the object point cloud.

Images