SIMULATION EXPERIMENT ON LANDING SITE SELECTION USING A SIMPLE GEOMETRIC APPROACH

Safe landing is an important part of the planetary exploration mission. Even fine scale terrain hazards (such as rocks, small craters, steep slopes, which would not be accurately detected from orbital reconnaissance) could also pose a serious risk on planetary lander or rover and scientific instruments on-board it. In this paper, a simple geometric approach on planetary landing hazard detection and safe landing site selection is proposed. In order to achieve full implementation of this algorithm, two easy-to-compute metrics are presented for extracting the terrain slope and roughness information. Unlike conventional methods which must do the robust plane fitting and elevation interpolation for DEM generation, in this work, hazards is identified through the processing directly on LiDAR point cloud. For safe landing site selection, a Generalized Voronoi Diagram is constructed. Based on the idea of maximum empty circle, the safest landing site can be determined. In this algorithm, hazards are treated as general polygons, without special simplification (e.g. regarding hazards as discrete circles or ellipses). So using the aforementioned method to process hazards is more conforming to the real planetary exploration scenario. For validating the approach mentioned above, a simulated planetary terrain model was constructed using volcanic ash with rocks in indoor environment. A commercial laser scanner mounted on a rail was used to scan the terrain surface at different hanging positions. The results demonstrate that fairly hazard detection capability and reasonable site selection was obtained compared with conventional method, yet less computational time and less memory usage was consumed. Hence, it is a feasible candidate approach for future precision landing selection on planetary surface.  Corresponding author


INTRODUCTION
Safe planetary landing is crucial to the space exploration mission.It largely depends on precisely and efficiently detecting and avoiding potential hazardous obstacles (like boulders, craters and steep slopes etc.), analysing the landing area in real time during the terminal phase of powered descent and selecting the safe landing site for the spacecraft.
Before the 1990s, the planetary surface landing missions were implemented basically in the absence of much prior knowledge of the targeted body (Bennett, 1972;Epp and Smith, 2007).There were several Mars landers been launched during recent 20 years (Ball et al., 2007).Up to now, among all these unmanned probes, for the first time, the automatic hazard detection and avoidance was realized during China's Chang'e III mission.
With the rapid development of deep space exploration activities, the future science-driven mission (like asteroid or comet exploration) will require soft landing within the potential high scientific value region (Furfaro et al., 2012;Kawaguchi, 2013;Rodgers et al., 2016;Witte et al., 2016).It happens frequently that the terrain conditions of these area are more complicated, this not only poses a great challenge to the spacecraft Guidance, Navigation, & Control Systems, but also requests more higher demand to the performance of landing obstacle detection and the efficiency of the landing site selection.Up to now, using multiple sensors data for hazard detection have been proposed (Brady et al., 2009;Neveu et al., 2015), including passive optical image data (Bajracharya, 2002;Cheng et al., 2001;Cohanim et al., 2013;Huertas et al., 2006;Mahmood and Saaj, 2015;Matthies et al., 2008;Woicke and Mooij, 2016;Yan et al., 2013), LiDAR data (Amzajerdian et al., 2013;Chakroborty et al., 2009;de Lafontaine et al., 2006;Johnson et al., 2002) and Radar data (Pollard et al., 2003), and so forth.
Nowadays, Several trials have been carried out during the course of NASA's ALHAT project to validate the hazard detection method using LiDAR sensor on-board a helicopter (Carson et al., 2013a;Carson et al., 2015;Carson et al., 2013b;Epp et al., 2014;Epp and Smith, 2007;Trawny et al., 2015;Villalpando et al., 2013).Likewise, NASA proposed to develop the 3D Imaging Flash LiDAR system to hazard detection for future planetary landing missions (Amzajerdian et al., 2011).Based on active LiDAR methods, JPL utilized the method based on active LiDAR for the avoidance of obstacles (Johnson et al., 2008).ESA also developed its own 3D Imaging Flash LiDAR system for the 2018 Lunar mission (Kerr et al., 2013;Parreira et al., 2013).
For safe landing site selection, the safe index of the landing site was proposed in literatures (Johnson et al., 2008;Simões et al., 2012;Simoes et al., 2009), but it was not reported on how to locate the safe landing spot.Shao et al. (2008) presented a safe landing site selection approach based on computational geometry and genetic algorithm, but obstacles were only treated as discrete circles or ellipses, therefore the suboptimal largest circles appeared sometimes.Considering the factors of landing safety, fuel consumption, and scientific returns, Serrano established the framework of multi-source data fusion for landing site selection based on bayesian theory (Serrano, 2006).Also taking account of landing safety, fuel consumption and touchdown performance, Cui et al. (2016) constructed safety index to select optimal safe landing spot.Cohanim et al. (2012) proposed the computing architecture using multiple weighted cost distance.Using this method, the suitable position of landing site can be determined, so it has high scalability and referable value.Ivanov et al. (2013) presented a framework for safe landing area selection based on probability reasoning, better dealing with the influence on landing site selection from the sensor noise and navigation positioning error.Method was proposed respectively in literature (Howard and Seraji, 2004;Seraji and Serrano, 2009;Serrano et al., 2005) from the angle of terrain classification based on multi-sensor data fusion and safe landing area selection.Câ mara et al. (2015) summarized the multi-sensor data fusion strategies for hazard detection and safe site selection for planetary and small body landings.Witte (2013) put forward a framework on hazard detection and avoidance maneuver on the basis of stochastic modelling.
This study focuses on the simple geometric approach for planetary landing hazard detection and safe landing site selection.The method has two potential benefits.On one hand, hazards is identified directly through the processing on LiDAR point cloud, there is no need to do the robust plane fitting and interpolation for Hazard DEM creation and analysis.On the other hand, hazards are viewed as discrete polygons, which is more in line with the actual situation, and thus more suitable landing site can been obtained.
The rest of this paper is organized as follows.The proposed methodology will be elaborated in Section 2. The preliminary results are presented and discussed in Section 3. Finally, this study is summarized in Section 4.

METHODOLOGY
The methodology in this paper consists of three steps, which are pre-processing step, landing hazard detection, and safe landing site selection.In the following sub-sections, details on steps of the proposed methodology are described.

Pre-processing
Due to the existing of various noises during the processing of the scanning, we applied the k-nearest neighbour algorithm to denoise the point cloud data.
The LiDAR point cloud has almost the same resolution and position accuracy over the entire experimental scene.So in this experiment, after denoising of the LiDAR data, the 3D points lied in the experimental field should be project to the 2D reference grid with side length s.Thus slope index and roughness index could be computed for the subsequent hazard detection.Suppose that P=(X, Y, Z) is a 3D point from the point cloud, where X-Y plane is parallel to the ground (grid), Z is elevation value of P above the ground.Then we can use equation ( 1) to calculate the corresponding grid coordinate of point P.
To do the calculation in equation ( 1), first we should discretize the whole area into equal size grids.Towards discretion of 2D space for LiDAR point cloud, average point space is a key index, namely nominal pulse space (NPS), with its inverse nominal pulse density (NPD).NPD is typically expressed as pulses per square meter (pls/m 2 ).Hence via the formula NPS = 1/ NPD from literature (Heidemann, 2012), NPS can be obtained.In the case of 6,000,000 points within the area of 15×4.2m 2 , the average point space (that is to say NPS) is about 3 mm.With regard to the computing the slope index and roughness index, it is vital to consider the number of sampled point lied in each reference grid.Normally we assume that the terrain point elevation obeys a normal distribution in each grid, actually we often use t-distribution to describe it.When the number of points exceeds 120, t-distribution is close to the normal distribution, thus it can be used to approximately substitute normal distribution.Hence, we use more than 10 times of NPS as the side length of the discretized grid (in this paper, side length is 4 cm) to ensure that there has sufficient points to computing the terrain characteristic.

Hazard detection based on simple geometric metric
In this paper, two metrics, namely slope index and roughness index, will be applied to evaluate the terrain surface for simulating hazard detection during the spacecraft powered descent process.
For each grid, based on equation (2), we employ 2D sliding window to calculate slope index, the result is expressed in degrees.
slope(i, j) atan((maxh-minh)/LanderSize) 180/π  (2) Where slope(i,j): the slope of grid(i,j) maxh: maximum elevation value of points within the lander size region centered at grid(i,j) minh: minimum elevation value of points within the lander size region centered at grid(i,j) LanderSize: the assumed lander size in our research atan: denotes the arctangent function π: denotes the circular constant Similar to slope index, we employ the variance of the point elevations within the same region to compute roughness index, the formula is defined as follows.
Where roughness(i,j): the roughness of grid(i,j) hk: the k th point elevation within the same region as computing slope index h : the mean elevation of all points within the same region as computing slope index n: the number of the points within the same region as computing slope index tolh: the maximum height above the ground of the hazard the lander can endured When there are insufficient points (such as, below 30), then NoData is assigned to the two indices for the corresponding grid.
After the acquisition of the slope map and roughness map, binary threshold segmentation on the aforementioned two types of map should be done according to landing demand on terrain condition (slope and roughness index).Do the OR operation on the two segmented map followed by morphological operations (such as open and close operation) and vectorization, the potential obstacles could be extracted.

Safe landing site selection based on Generalized Voronoi Diagram
Inspired by the idea of Maximum Empty Circle usually used in robot motion planning, we employed generalized voronoi diagram to carry out the safe landing site selection under the condition on the hazards distribution.Here hazards are treated as the dispersed polygons.Our work is to find the circle with the largest radius between the polygons (not intersect with any polygons).
To construct the generalized voronoi diagram, here the approximated method was used.In more detail, the method is as follows.First, each edge from each polygon will be discretized with a specified step  , resulting in so many points on the polygons boundaries.Then these scattered points will be used to computing voronoi diagram.Next those voronoi edges which intersect with the polygons and the envelope of the LiDAR point cloud should be removed.Only those skeleton lines of voronoi diagram could be left.Finally via calculating the intersect points of these skeleton line, the candidate site with the maximum empty circle can be gained.
For the ultimately quantitatively evaluating the most suitable landing site, a performance index (Takahashi et al., 2013) incorporating slope index, roughness index and the radius of the inscribed circle on each candidate site will be computed.Finally, the site with the lowest performance index is the recommend landing point.The formula for defining the performance index is expressed in equation ( 4).With respect to defining the weight of each single index, we adopt the following rules.For the weight on slope, we use the reciprocal of the permissible slope value (here is 15°).Owing to the higher sensitivity characteristics of roughness index to the landing safety under lower roughness, after multiple tests, we introduce the formula Wr=e -1/(10*r) to determine the weight on roughness index where symbol "e" and "r" denote Euler's number and roughness value respectively.Similar to calculate the weight of roughness index, we use formula Wd = 0.8 d/LanderSize to determine the weight on radius of inscribed circle considering the rapidly dropping characteristic on landing safety along with the increasing of the distance between the candidate site and the nearest obstacle.In a word, our strategy is that the slope and roughness of selecting the landing site is as small as possible, meanwhile the distance departing the hazard is as far as possible.In addition, given equivalent slope and roughness, the algorithm prefers the site with more inscribed circle radius.Furthermore, appearing the n site with all three equivalent indexes, we should take into account other factors, such as the fuel cost, the site's scientific value, etc.These has gone beyond the research of this paper.

EXPERIMENTAL RESULTS
In this research, we used a RIEGL VZ-400 commercial scanner mounted on the rail above the simulated lunar surface to acquire the 3D terrain point cloud at 3 hanging positions with the same height above the ground.These point cloud from different stations then can be merged together for terrain analysis and hazard detection.
After grid processing, we use the aforementioned method to generate the slope map and roughness map, presented in Figure 1 and Figure 2 respectively.From the result, we can see that the vast majority of hazards can be detected via the proposed approach regardless of craters, boulders, hilly regions, but the boundary of extracted hazards could be slightly larger than the actual size when use slope index.It may be caused by less using the statistics from the points around the corresponding grid.It only employs the maximum and minimum elevation value to compute, not consider the spatial distribution of these points, so sometimes the slope degrees could be exaggerated.On the other hand, the larger grid size (4 cm) will also caused the uncertainty of obstacle boundary.After all, the size of some rocks in the experimental field only has several centimetres.Hence, in our future work, this drawback will be averted.Similar situation has less appeared on roughness index, especially for boulder detection.The detected hazards are shown as the red polygons in Figure 3.In the stage of safe landing site selection, via the method described in section 2.3, we subdivide the whole area as the blue lines in Figure 3 shows.Then, by searching all points intersected of 3 blue lines, we select 10 candidate sites according to the distance to its nearest obstacle.Each of these 10 selected candidate sites is surrounded by a circle with green color, respectively.Next we calculate the performance index for each site using the method mentioned in section 2.3.By accumulating the sub-index scores listed in Table 1, we know that the final recommend landing spot is located at site 5, its total score is 0.983, the lowest value.1.Each site's sub-index value and its corresponding score

Site
The calculation for safe landing site selection occupied most of the time of the whole process.The reason is that removing redundant voronoi edge and searching candidate sites could consume too much computing resources.Actually, removing additional edges only promoted visual effect, but had no help on site selection.So in practical applications, we will not do that.Towards searching candidate sites, by increasing the step ε , we can induce the points during the polygons discretization resulting in the inducing the edges of voronoi diagram.Indeed, we can even only use the vertices of polygon to computing the voronoi diagram.Supplemented by some parallel algorithms, our methodology could satisfy the real time processing demand on planetary landing hazard detection and safe landing site selection.

CONCLUSION
In this paper, a simple geometric methodology for planetary landing hazard detection and safe landing site selection is proposed.And the validation on this method using the simulated planetary surface terrain data is also implemented.
During the whole experimental process, we find that the error introduced in each step could lead to the uncertainty to a certain extent.So in future, we would carry out this work considering the uncertainty contained in this sensor data.This will certainly promote the quality and reliability of the landing mission.
index on corresponding candidate site r: roughness index on corresponding candidate site d: radius of inscribed circle on corresponding candidate site Ws: the weight of slope on PI Wr: the weight of roughness on PI Wd: the weight of distance on PI

Figure 1 .
Figure 1.Hazard presented in slope map

Figure 3 .
Figure 3. Final result on hazard and safe landing site selection