For future lunar exploration and planetary missions, the digital elevation model (DEM) of the target object may not be well prepared before the mission, so developing a new robust crater detection algorithm (CDA) without prepared high-precision DEM is needed to meet the requirements of a high-reliability and high-precision detection and navigation system. In this paper, we presented a new robust lunar CDA method based on maximum entropy threshold segmentation. By calculating the entropy distribution of the ternary image, the threshold for retaining the maximum amount of image information is selected adaptively, a variety of evaluation indicators are proposed, and a multiple-indicator constraint matrix is constructed to realize the extraction and fitting of the craters. The proposed method has the following advantages: (1) it has strong robustness and is capable of extracting complete craters under multiple illumination conditions, which makes it suitable for the extraction of large-scale planetary and lunar images; (2) the extracted crater edges are clear and complete and do not merge with the surrounding environment edge; and (3) it avoids the problem of parameter sensitivity that is present in a traditional CDA algorithm. The proposed method was verified using an image taken by the Chang’e-2 lunar probe, and a comparison with the traditional method based on morphology and adaptive Canny edge detection shows that the number of craters detected increases by more than 35%, while the computational efficiency is improved by more than 40%.

1. Introduction

Planetary landing exploration is a key approach to study the formation and evolution of the solar system, the development and utilization of space resources, asteroid impact defense, and other major scientific issues. In the past decades, many countries have carried out a series of successful exploration missions for the Moon (Apollo [1, 2], Chang’e [37]), Mars (Tianwen-1 [8], Curiosity [9], and Opportunity [10]), asteroids (Rosetta [11]), and other celestial bodies. In general, terrain with high scientific value is more complex, and landing is more difficult; therefore, crater detection algorithms are important for lunar exploration and planetary missions. Chang’e-4 is the only exploration mission in a complex area that has successfully landed on the far side or back of the Moon [12, 13]. Compared with the traditional lunar sea area on the front of the Moon, the terrain on the back of the Moon is more complex and lacks high-precision terrain elevation data. Furthermore, traditional communication technology based on deep space networks presents greater risks for a mission due to communication delay problems [14], and therefore, algorithms need to be more robust.

To overcome the limitation of the traditional deep space network, it is necessary to develop a robust autonomous detection method for lunar craters under any illumination and construct autonomous hazard detection and avoidance (HDA) technology based on optical information for accurate navigation and autonomous landing-point selection for the lunar surface high-precision descent landing guidance, navigation, and control (GNC) system.

To detect craters, many crater detection algorithms (CDAs) have been studied by researchers all over the world. Some methods have focused on positive detection rates, while other methods pay more attention to speed or robustness. Studies by Stepinski et al., Delatte et al., Silburt et al., and Emami et al. [1518] used machine learning and convolutional neural networks (CNNs) separately to process digital elevation models (DEMs). Studies by Wang et al. [19] presented a high generalization performance transfer learning-based CDA. Some traditional methods are also used in this field [20], such as the Hough transform, Canny edge detection [21, 22], cascade decision forest [23], solar illumination direction, and texture analysis [24]. Jia et al. extracted more all craters larger than 200 m of Chang’e-5 landing area based on digital orthophoto map generated from more than 700 Lunar Reconnaissance Orbiter Camera, Narrow Angle Camera [25]. Solarna et al. proposed a marked point process CDA method using graph cuts and decimated wavelets [26]. Bandeira et al. detected subkilometer craters and achieved a high detection rate through texture analysis and shape information [27]. Galloway et al. [28] used a Hough transform and Canny edge detection for processing high-resolution images, which required a large number of calculations. Lee and Hogan [29] proposed an automated CDA with human-level performance by multiple neural networks, but the method used not only DEM but also thermal infrared imagery, so the method did not have universal applicability. Francis et al. [30] presented a dataset of craters for CDA based on deep learning. Cui et al. [31] directly used the crater curve matching for the planetary landing navigation system and obtained superior accuracy. Chen et al. [32] presented a CDA based on terrain analysis and mathematical morphology applied to DEM with a resolution of 100 m. Lee et al. [33] presented a CNN-based CDA, but the method needs to determine detection probabilities in advance for mission planning.

However, many existing CDAs have different limitations. Methods based on machine learning and CNNs require suitable train sets [34]. The Hough transform has poor adaptability because only circular craters can be detected. Furthermore, texture analysis cannot be adjusted to a variety of lighting conditions. Therefore, these methods have problems with robustness or speed; CDA based on image or optical information still have research value. In the view of CDA using unsupervised algorithms and CNNs, Emami and other researchers made a detailed research and comparative analysis [35]. From the study of Emami et al. [35], we can conclude that CDA based on convex grouping and interest points performed better in combining with CNNs for hypothesis verification. Some research results indicate that CDAs based on optical information processing still have great research values and have application prospects in the field of navigation. Maass and other researchers [36] used mainly image segmentation method to detect craters and applied research results to CNAV, which achieved good performance. Furthermore, Maass [37] firstly investigated image illumination direction in segmentation-based CDA for spacecraft navigation and then pointed that the navigation estimator can be considered robust and accurate enough in proposed CDA.

To meet the needs of reliable on-orbit real-time navigation and HDA systems, a robust and fast CDA is proposed in this paper. First, the gray image is filtered by an adaptive Gaussian filter to remove image noise while retaining edge information. Then, the image is divided into highlighted, shadowed, and average gray areas via maximum entropy threshold segmentation to reduce the sensitivity of the results to light sources and to improve the reliability of Canny edge detection. Finally, a multiple-indicator constraint method is proposed to match two edges that belong to one crater, and an ellipse is used to fit the two edges.

The structure of this paper is as follows. In Section 2, an image ternary segmentation method based on a maximum entropy threshold is proposed to deal with the disadvantages of traditional and adaptive Canny detection in CDAs; in Section 3, a multiple-indicator constraint matrix is constructed to accomplish crater edge matching and ellipse fitting; in Section 4, simulations are presented, and a lunar image is used as an example to demonstrate the detection quality and computational cost of the proposed model; finally, our conclusions are drawn in Section 5.

2. Crater Edge Detection Based on Maximum Entropy Threshold

In this study, craters were detected by calculating the gradient of the image, the maximum entropy threshold was used to segment the image before Canny edge detection, and then the craters could be detected after the clear edges were extracted. Compared to traditional methods, the crater edges extracted by our method are clear and obvious and do not merge with the surrounding environment; they do not need morphological processing such as dilation and erosion and can be directly used for crater edge matching and ellipse fitting.

2.1. Adaptive Canny Edge Detection

The craters on the lunar surface form highlighted and shadowed areas that are quite different from the surrounding environment and form obvious gray gradient edges under direct sunlight. Therefore, craters can be effectively detected by edge detection. There are many methods of edge detection, such as the Roberts operator, Prewitt operator, and Canny operator, but the Roberts and Prewitt operators have two disadvantages. The first disadvantage is that they are easily affected by noise and can detect multiple edges; the other disadvantage is that the scale of these operators is fixed, which is not conducive to detecting edges of different scales. Compared with the Roberts and Prewitt operators, the Canny edge detection can effectively suppress multiresponse edges and can improve the accuracy of edge detection via nonmaximum suppression. Therefore, in this study, we used the Canny edge detection.

First, for denoising the image and preserving the edge information of craters, a Gaussian filter is commonly used. The Sobel operator is then used to calculate the gradient of the image as follows: where are Sobel operators along the -axis and -axis, respectively. The norm and direction of the image gradient can be calculated as follows: where is the convolution window of the Canny edge detection. Finally, the strong and weak edges are extracted by nonmaximum suppression and dual-threshold detection. Dual-threshold, high-threshold , and low-threshold are the extraction standards of crater edges. The edges with gradient intensity bigger than high-threshold are sure to be crater edges, called strong crater edges, and the edges with gradient intensity smaller than low-threshold are sure to be noncrater edges. The edges whose gradient intensity lies between high-threshold and low-threshold are classified as weak crater edges or nonedges by their connectivity. As shown in Figure 1, the gradient intensity of edge A is bigger than , so edge A is sure to be a crater edge. The gradient intensity of crater B is between and , and edge B is connected to edge A, so it is a weak crater edge. Edge C is not connected to any crater edges, so it is not a crater edge.

It can be concluded that the selection of dual-thresholds, high-threshold , and low-threshold will seriously affect the effect of the Canny edge detection. Various values of and were selected to analyze a 7 m resolution digital orthophoto map (DOM) image of the lunar surface taken by Chang’e-2. The influence of different and values is shown in Figure 2. Figure 2(a) shows the original lunar image including regions A and B; Figure 2(b) shows the original region A image and different Canny edge detection results for region A with different ; Figure 2(c) shows the original region B image and different Canny edge detection results for region B with different ; Figure 2(d) shows different Canny edge detection results for region B with different .

We can infer from the results that the crater detection algorithm based on Canny has the following two disadvantages due to the parameter sensitivity: (i) the edges detected by Canny usually connect with surrounding edges, as shown in Figure 2(b); (ii) it is hard to choose an appropriate dual-threshold autonomously under different light conditions, which cannot meet the demand of a lunar exploration mission. To address these problems, He [38] proposed an adaptive dual-threshold selection method in the Canny edge detection algorithm, i.e., through the histogram statistics of the gray image after nonmaximum suppression, determine the gradient corresponding to the maximum number of pixels and the gradient amplitude mean square deviation of the maximum number of pixels as follows: where is the number of histogram bins and is the number of pixels that fall into bin whose gray value is .

The high-threshold is as follows:

Remove the pixels whose gray value is stronger than , and recalculate and using the steps above. The low-threshold is as follows:

This adaptive dual-threshold selection method was used to detect crater edges, as shown in Figure 2(a), and the results are given in Figure 3. For the same shooting area, by selecting different regions to analyze, the dual-thresholds selected by the traditional Canny edge detection method are quite different, which were easily affected by the camera field of view and the shooting area.

First, dual-thresholds were chosen to extract crater edges, when only region A or region B is analyzed, and the corresponding values are different from the of the whole image, as shown in Table 1, which leads to different results; this is not stable and robust, as shown in Figures 3(a) and 3(b).

Figure 3(a) shows the Canny edge detection results for the whole image; Figure 3(b) shows the Canny edge detection results for regions A and B, where (i) is the enlargement of region A in Figure 3(a), (ii) is the Canny edge detection result for region A alone, (iii) is the enlargement of region B in Figure 3(a), and (iv) is the Canny edge detection result for region B alone.

From Table 1 and Figure 3 results, we can infer that although the adaptive dual-threshold selection method in the Canny edge detection algorithm by Jiang [38] can select dual-threshold autonomously and extract the crater edges, the edges are usually broken and discontinuous and are sometimes connected with the surrounding area, which is not conducive to further processing. In those circumstances, the edges of craters cannot be extracted effectively for ellipse fitting; therefore, an image ternary segmentation method based on the maximum entropy threshold is proposed, which can effectively solve this problem.

2.2. Image Ternary Segmentation Method Based on Maximum Entropy Threshold

Under sunlight, the shadowed area, the highlighted area, and the surrounding environment can be seen on the lunar surface. The ternary segmentation method based on maximum entropy threshold is proposed to distinguish these three parts. After obtaining a filtered image, examining all gray colors, and calculating image entropy, the image can be ternary segmented by two thresholds, , which correspond to maximum entropy.

The entropy of the image is where is the probability density of the pixels whose gray value is .

We select to segment the image and divide the image into shadowed areas , background areas , and highlighted areas . The corresponding cumulative probability is as follows:

is satisfied as follows:

The entropies of those three areas are as follows:

The total entropy of the image is the sum of as follows:

We obtain a maximum by keeping the image information via selecting appropriate values of .

By choosing different values to segment the image, as shown in Figure 4(a), this leads to the distribution of the segmented image entropy, as shown in Figure 5. By finding the thresholds corresponding to the maximum entropy, the ternary segmentation image result based on the maximum entropy threshold is as shown in Figure 4(b).

We then detect the edges of craters after getting a ternary segmentation image. For the two different thresholds , the Canny edge detection produces results similar to those shown in Figure 6. From the results, we can infer that this method is robust in terms of dual-thresholds . No matter how the dual-thresholds are chosen, clear and continuous edges can be detected after image ternary segmentation, which helps to match semicrater edges and ellipse fitting in the next step.

3. Crater Edge Matching and Ellipse Fitting Based on Multiple-Indicator Constraints

After image ternary segmentation and crater edge detection, we need to match pairs of edges that belong to the same crater. To achieve this goal, we first use a series of feature points to represent the semicrater edge, as Figure 7 shows.

One semicrater edge can be represented by four feature points; they are the vertices and , geometric center , and midpoint of vertices . Then, the orientation vector of the semicrater edge can be defined as follows:

Taking into account the length, distance, orientation, brightness, and other constraints comprehensively, we match two semicrater edges by constructing multiple constraints . When , two edges match best, and when , two edges match worst; they do not match at all. The constraints are constructed as follows: (1)Normalized length similarity constraint

The lengths of two matched edges should not be widely different. The bigger the constraint, the better the two semicrater edges match. where are pixel lengths of the edges to be matched (2)Brightness constraint

Under sunlight, two semicrater edges belonging to the same crater should be in different highlighted and shadowed areas. The gray value of the geometric center is used to distinguish the two semicrater edges to be matched as follows: where are the two geometric centers’ gray values of the edges to be matched. indicates that two edges belong to the highlighted area and the shadowed area, respectively (3)Normalized orientation constraint

Two semicrater edges belonging to the same crater should be located in opposite directions, which means two orientation vectors are at an obtuse angle. The bigger the constraint is, the better the two semicrater edges match. (4)Distance constraint

Two semicrater edges belonging to the same crater should not be far apart. We use the geometric mean to measure the distance between two semicrater edges. The distance is divided into horizontal distance and vertical distance as follows: where is the sunlight vector and and are scale parameters. In general, the distance between two semicrater edges belonging to the same crater is shorter than half the perimeter of the crater. And the geometric mean is usually shorter than half the perimeter of the crater. The relationship of and can be written as follows:

In this paper, we chose (5)Normalized shape constraint

According to the above constraints , to match the edges of craters, there are two matching results, that is, X-shape results and O-shape results; O-shape results are real craters, and X-shape results usually come from stone shadow influence. To avoid this influence, the vector from the center of the shadowed area to the center of the highlighted area should be consistent with the direction of the sunlight vector as follows:

For any two semicrater edges, the multiple-indicator constraint can be calculated to represent the similarity between them as follows:

For semicrater edges, we construct a matching similarity multiple-indicator constraint matrix as follows:

is a symmetric matrix with diagonal elements that are zero; means the similarity between edges and .

First, we choose two largest elements in , which are located in coordinates and . Then, we remove the other elements in whose positions are row and column , and row and column . We name the new multiple-indicator constraint matrix . Then, we choose the new two largest elements in and repeat the above steps until the largest element is lower than the matching threshold . After the semicrater edge matching is complete, we ellipse fit pairs of edges that belong to one crater.

Craters can be described mathematically as follows:

We use the least-squares method (LSM) to calculate . Since the projections of the craters are ellipse, it is necessary to exclude the case for which has nonpositive eigenvalues. After ellipse fitting, the candidate for the crater center is the following:

4. Results and Discussion

Crater detection is basic to optical navigation, crater reconstruction, and HDA. To verify the effectiveness and robustness of the method proposed, two experiments are carried out in this article. The first experiment verifies that the proposed CDA works well and has robustness under different illumination directions for the same area. The second experiment verifies that compared with the traditional method, the proposed CDA can detect more craters and has lower computational efficiency.

4.1. Simulation Verification on Sand Table under Different Illumination Directions

Crater detection is basic to optical navigation, crater reconstruction, HAD, and landing site selection. A simulated lunar sand table is constructed to verify the robustness of the proposed CDA. The parallel light source is set at different positions of the sand table to simulate different sunlight illumination directions. The environment around the sand table is processed to avoid the influence on the detection of sand table craters. The experimental conditions are shown in Table 2.

Hardware-in-loop simulation experiment is constructed by a front-facing camera, sand table, parallel light source, and data switching exchanges. Figure 8 shows the camera and the lunar sand box used in the experiment. Figure 9 shows the hardware-in-loop simulation environment.

The parallel light sources are, respectively, set on the left and right sides of the lunar simulation sand table, with a height of about one meter, to simulate the different light directions of the Sun. After Gaussian smoothing of the sand table images under different light conditions, maximum entropy threshold segmentation is used to form the different highlight and shadow areas of the crater, and then the craters are detected. Detected results are as shown in Figure 10.

Figure 10 shows the following: (a) parallel light source is set on the left of sand table; (b) maximum entropy threshold segmentation of (a); (c) detected craters of (a); (d) parallel light source is set on the right of sand table; (e) maximum entropy threshold segmentation of (d); and (e) detected craters of (d).

We can make some discussions from Figure 10. (1)Seven craters are detected in each image. Crater No. 1~6 are all detected in two images because their shadows are easily segmented from the surrounding environment. Crater No. 7 and crater No. 8 are all covered by the shadow of crater No. 2, so the proposed CDA can detect only one of them in each direction of sunlight. This is caused by the shape and distribution of the craters on the sand table, not by our proposed method(2)Although crater No. 1~6 are all detected in two images, the ellipse shapes of the craters are different. That is because shadows of the same craters under different sunlight directions are usually quite different as shown in Figures 10(b) and 10(e)(3)CDA proposed can work well under different illumination directions, which means our CDA method has strong robustness.

4.2. Simulation Verification in Lunar Image

According to the lunar surface image data released by Chang’e series of probes, we can obtain optical images of the lunar surface with different resolutions. In this paper, a DOM (7 m) gray image, detected by the Chang’e-2 lunar probe, named CE2_GRAS_DOM_07m_C101_66N172W_A.tif was selected as the experimental data to extract large and small craters. We can download it from https://moon.bao.ac.cn/searchOrder_dataSearchData.search. The image information and PC parameters are shown in Table 3.

There were 1337 craters counted in Chang’e-2 DOM (7 m) image, and the number 1337 was used as the ground truth (Table 4). To detect craters fast and accurately, a sliding window sized is used to detect small craters.

At a matching threshold of , 1170 craters can be detected in this image without false positives. The detection results 1170 craters are shown in Figure 11. The image resolution is 7 m per pixel, and we can calculate the average diameter of each detected crater. For different sizes of crater, our method has different extraction rates; the bigger the crater, the higher the extraction rate. Through analysis, the main reasons for the decline in the extraction rate of smaller craters are as follows: the slope of lunar surface with small craters is too small, the changes of light and darkness are not obvious, and there is a large degree of overlap with the surrounding craters.

To make the detection results more obvious, Figure 12 shows the ground truth of all craters by manual counting. The red circles represent the center of the undetected craters, and the blue circles represent the crater detected by the method proposed.

Table 5 shows the extracted information for craters whose average diameters are larger than 700 m. The recorded information includes the average diameter, center coordinates, and direction angle ().

Many researchers [39, 40] proposed different parameters to evaluate the performance of the CDA. These methods work by counting craters detected in different states, i.e., true positives (TPs) corresponding to real detected craters and false positives (FPs) corresponding to detected craters that do not exist. A false negative (FN) represents the number of existing craters that are not detected. After obtaining the TP, FP, and FN, the quality factors can be calculated, which include the detection percentage (), the quality percentage (), the branching factor (), and the precision () as follows:

shows that our method proposed is better than random guess; , , and can help us to choose the best threshold.

Different thresholds were chosen, i.e., from 0.1 to 0.9, and then the detected crater numbers and quality factors were calculated using our proposed method. The results are shown in Table 6. The symbol represents the optimal value.

The detection percentage curve as shown in Figure 13 was drawn using the data in Table 6.

We can make the following conclusions from Table 6 and Figure 13. (1)When , although the number of detected craters is even more than manual counting ground truth, TP and detection percentage do not reach 100%. There are two reasons, the first reason is that some structures of old craters are broken, and the other reason is that some semicrater edges mismatch to other craters(2)When the matching threshold increases, the constraint for matching semicrater edges becomes stronger, and therefore, the number of extracted craters decreases; TP and FP decrease at the same time(3)When the matching threshold increases, FN increases along with it, because of irregularly shaped, incomplete crater edges, or craters that did not have an obvious change in gray values(4)The CDA achieves the best result only by choosing the appropriate threshold. When the matching threshold increases, quality factors and become better but the number of detected craters decreases. When , there is no FP which matters with navigation, and and are all the best value, which means is an appropriate threshold(5)Design an evaluation parameter to comprehensively consider various quality factors.

When , , which is larger than other parameters that correspond to other .

The area under the curve (AUC) in Figure 13 could be calculated to evaluate the performance of the proposed method; the higher the AUC, the better the proposed method. (i) means that the proposed CDA is perfect, and it can detect all craters without any FP and FN(ii) means the proposed CDA is better than random guess. If the matching threshold is set appropriately, the proposed method can work very well(iii) means that the proposed method works the same with random guess(iv) means that the proposed method works even worse than random guess.

The AUC of the proposed method is 0.9192, which is much higher than 0.5, which means that our method works well and has robustness for crater detection.

Table 7 shows a comparison of the quality factors of our proposed method with He’s CDA [38]. For the same lunar surface picture, our method obtained a better detection percentage, quality percentage, and branching factor, due to a larger number of TPs, smaller number of FPs, and smaller number of FNs.

Compared to He’s method and other traditional methods, the proposed method improves TP by 35.3%, and FP and FN are much smaller than He’s method. Quality factors are also better. It can be concluded that the proposed method is a more effective large-scale image CDA.

In Table 8, we compared the quality factors using a published lunar crater database [24, 41]; there are 4333 craters at the landing site.

By dividing the crater detection algorithm into substeps (image preprocessing, ternary segmentation, adaptive dual-threshold selection, edge detection, edge matching, and ellipse fitting), we can compare the difference between the proposed method and He’s method. Therefore, we selected a subgraph in Figure 11 of , calculated the time for each step, and analyzed the computational cost; the result is shown in Figure 14 and Table 9.

Compared with He’s method, the computational cost of the proposed CDA is significantly reduced. The time for image preprocessing is basically the same; He’s method takes a long time to adaptively select two thresholds, while ternary segmentation based on the maximum entropy threshold, as proposed in this paper, takes a short time. In addition, as compared with He’s method, using our method, the edges detected are continuous, smooth, and without broken edges, and the time for edge matching is also shorter. Because our method has less mismatching and strong robustness, the ellipse fitting time is also slightly less than for He’s method. Overall, compared to He’s method, the computational efficiency of He’s method can be improved by 40.1%.

5. Conclusions

In this paper, a robust crater detection algorithm based on a maximum entropy threshold was proposed. The proposed method can effectively improve the robustness of crater detection algorithms for Mars, the Moon, asteroids, and other celestial bodies and can enhance the quantity and reliability of crater extraction. It can also be used for the vision navigation tasks of Mars, the Moon, asteroids, and other celestial bodies, hazard avoidance, and path planning. In our method, the entropy distribution is first calculated using two thresholds of the detected image, and the maximum entropy threshold of the image is calculated to make the image ternary. Then, the edge of the ternary image is detected to complete the edge extraction of the crater; a multiple-indicator constraint matrix is constructed, and the extracted edge of the crater is matched and fitted by comprehensively considering the length similarity, brightness, and other normalized factors, so as to complete the extraction of the crater.

A simulation lunar sand table and a lunar surface image captured by the Chang’e probe were used to verify the proposed algorithm. The simulation results show that, compared with traditional crater detection algorithms, the extraction results and quality factors of the proposed method are significantly improved, the robustness is good, TPs increase by more than 30%, and the computational efficiency is improved by more than 40%. The method therefore has good development prospects for exploration missions to Mars, the Moon, asteroids, and other celestial bodies.

Data Availability

In this paper, a DOM (7 m) gray image, detected by the Chang’e-2 lunar probe, named CE2_GRAS_DOM_07m_C101_66N172W_A.tif was selected as the experimental data to extract large and small craters. We can download it from https://moon.bao.ac.cn/searchOrder_dataSearchData.search.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This work was supported by the Fourth Batch of Manned Space Preresearch Projects (Grant No. 18123060201).