MAPPING LANDSLIDES IN LUNAR IMPACT CRATERS USING CHEBYSHEV POLYNOMIALS AND DEM’S

: Geological slope failure processes have been observed on the Moon surface for decades, nevertheless a detailed and exhaustive lunar landslide inventory has not been produced yet. For a preliminary survey, WAC images and DEM maps from LROC at 100 m/pixels have been exploited in combination with the criteria applied by Brunetti et al. (2015) to detect the landslides. These criteria are based on the visual analysis of optical images to recognize mass wasting features. In the literature, Chebyshev polynomials have been applied to interpolate crater cross-sections in order to obtain a parametric characterization useful for classification into different morphological shapes. Here a new implementation of Chebyshev polynomial approximation is proposed, taking into account some statistical testing of the results obtained during Least-squares estimation. The presence of landslides in lunar craters is then investigated by analyzing the absolute values off odd coefficients of estimated Chebyshev polynomials. A case study on the Cassini


INTRODUCTION
Geological slope failure processes have been observed on the Moon surface for decades (Pike, 1971;Lindsay, 1976;Xiao et al., 2013), nevertheless a detailed and exhaustive lunar landslide inventory has not been produced yet.As a part of the collaboration in the 'Moon Mapping Project' between Italy and China (see Giommi et al., 2016), a working group has been established to focus on the detection, characterization and mapping of landslides in impact craters.In Brunetti et al. (2015), some criteria based on the visual analysis of optical images to recognize mass wasting features have been proposed.On the other hand, the need for a more objective automatic method is required in order to make the recognition process independent from the subjective interpretation and to carry out an exhaustive search.In Mahanti et al. (2014) the Chebyshev polynomials have been applied to interpolate crater cross-sectional profiles derived from a DEM in order to obtain a parametric characterization useful for classification into different morphological shapes.The presence of landslides is then recognized by analyzing the asymmetry of crater's profiles, which can be detected from the analysis of odd coefficients of the Chebyshev polynomials (Mahanti et al., 2015).

(*) corresponding author
In this paper the use of Chebyshev polynomials for detection and classification of landslides is continued.The aim here is to present the implemented methodology and to add new statistical tools for the analysis of Chebyshev interpolation.A case study consisting in the Cassini A crater has been considered.In this crater, a large slope failure has been detected, so that it may be used for validation the algorithms for the analysis of crater cross-sectional profiles.The LROC WACGLD100 DEM has been used for the initial set up of the methodology (resolution 100 m x 100 m), NAC (Narrow Angle Camera) images have been analyzed to investigate the lunar ground surface texture.

LANDSLIDES IN IMPACT CRATERS
The visual recognition and mapping of lunar landslides takes advantage of the same criteria commonly used by geomorphologists to detect and map landslides on Earth (e.g., Guzzetti et al., 2012).In particular, remotely-sensed data have been demonstrated to offer unprecedented opportunities for mapping geo-hazards over large areas (Scaioni et al., 2014).Remote-sensing data are even more important when dealing with extraterrestrial bodies, because in such a case they are the major -if not unique -data source.Recently, some typical criteria that have been normally applied to map landslides on Earth, were recently used to identify and draw landslide maps in Valles Marineris, Mars (Brunetti et al., 2014).Figure 1 shows an example of a landslide mapped in a lunar crater (Brunetti et al., 2015).For the detection of the landslides inside the impact craters, the QuickMap™ web interface and the open source Java Mission-planning and Analysis for Remote Sensing (JMARS) software; a Geographic Information System (GIS) developed by the Arizona State University (http://jmars.asu.edu/) are used.In order to focus on likely post-impact failures (i.e., those non-related to impact cratering), complex craters are not included in the analysis.The reason is that complex craters usually have terraced crater walls that form following the collapse (e.g., Melosh, 1989) and prevent the recognition of post-impact landslides.

The concept
Landslides in impact craters occurred after the formation of the main cavity as a consequence of soil instability, Moonquakes, direct new impacts from meteorites or seismic waves generated by other impacts in the nearby.The analysis of the impact crater shapes, either in 2-D and 3-D, is then fundamental for understanding both impact and post-impact processes.In Mahanti et al. (2014) an analysis of methods introduced for quantitative morphological study of crater shape is reported together with a review of the related literature.On one side, methods based on the extraction of some crater measurable features (diameter and depth in different locations, circularity, slope, etc.) and qualitative characteristics (central peaks, surface texture, asymmetries, etc.) have been proposed (see, e.g., Pommerol et al., 2012).These methods rely on the assumption of a theoretical model that should be applied to compare groups of similar craters.Unfortunately, such methods do not offer a detailed description of the crater shape to be analyzed for the detection of the surface degradation processes (Mahanti et al., 2015).
On the other side, polynomial approximations have been used to describe crater elevation cross-sectional profiles.These methods can be classified as data-driven, since they do not need any a priori model to be assumed.Since the approximation of more complex shapes of the profiles can be done by simply increasing the order of the approximating polynomial, this solution is potentially efficient in the case of craters affected by soil degradation (see an example in Fig. 1).The approximation level depends on the degree of the adopted polynomials: the terms that are omitted give rise to the so called truncation error, whose magnitude is related to the specific implemented polynomials.
In Mahanti et al. (2014) the Chebyshev polynomials have been used (see Mason and Handscomb, 2010) for approximating craters' cross-sectional profiles.These are a series of orthogonal polynomials, each of them featuring a unique and uncorrelated shape with respect to any other members of the series.More details will be given in next subsection 3.2.The approximation of each profile is accomplished by considering a cross-sectional length as twice the crater diameter from rim-to-rim.This distance is the normalized between -1 and +1, because this is the domain of Chebyshev polynomials.In this interval, any arbitrary continuous function can by approximated (Gautschi, 2004).
In the case under consideration, the function to approximate is the discrete crater profile f(x), being x the sample direction.

Mathematical background
Among four different types of Chebyshev polynomials, in the abovementioned paper it is suggested to implement the so called Type I for approximating crater profiles.This is motivated by the greater simplicity of coefficients related to this representation.The formulation of polynomials' basis functions is based on a recursive series: where Tn(x) is the polynomial piece of order n.The pieces of order 0 and 1 are T0(x) = 1 and T1(x) = x, respectively.In order to approximate a function f(x), a linear combination p(x) of the basis functions is adopted: where M is the degree of the Chebyshev polynomial and Cn are the coefficients that modulates the amplitude of each basis component.Coefficients Cn are estimated on a Leastsquares basis in order to fit with real profile data, as discussed in subsection 3.3.The residual approximation error o(x M ) is equal to the sum of the missing terms after degree M that are not considered in the approximation.
As it results from Eq.'s (1) and ( 2), in the Chebyshev polynomial series even (symmetric w.r.t.vertical axis) and odd (anti-symmetric) basis functions alternatively appear.Consequently, the size of odd coefficients may express the degree of asymmetry of the approximated crater profile.

Relevant properties of Chebyshev polynomials
Several properties make the Chebyshev polynomials particularly efficient for approximating crater cross-sectional profiles.These could be summarized in four main points: 1. Orthogonality of basis functions and uncorrelation among the estimated coefficients: even though the total number of adopted coefficients in Eq. ( 2) is different, the estimated values of the lower order coefficients is the same.These properties are important because make independent the estimated coefficients from the specific estimation process, so that they can be compared in a meaningful way.Indeed, lower numbered coefficients Tn(x) have a larger impact in the approximation of the crater profile; 2. Chebyshev polynomials results in the smallest maximum deviation with respect to the interpolated function, i.e., the crater cross-sectional profile in such a case (see Mason and Handscomb, 2010); 3. Extrema values of pM(x) always occur at some specific positions on the reference axis (x = -1, 0, +1).This property makes easier to link the estimated polynomials' coefficients with the geometry of the crater; 4. Correlation among the lower order coefficients as well as some combinations of coefficients with important morphological properties of the crater and its surrounding (average crater profile elevation, local topographic gradient, crater depth, etc. -see Mahanti et al., 2014 for a complete list).This does not mean that morphological features can be directly measured, but that a set of objective numerical indicators can be obtained through a repeatable almost automatic process; and 5. Detection of asymmetry in the crater cross-sectional profile on the basis of the analysis of odd polynomials.
Taking advantage of these properties, Mahanti et al. (2014) demonstrated that lunar crater cross-sectional profiles could be approximated by using Chebyshev polynomials.Using the first 17 coefficients (M=16), it is possible to describe in a compact standard format the profile, and also to operate a classification of differ crater shapes.
In Mahanti et al. (2015) some initial results of the analysis of asymmetric components of Chebyshev polynomials demonstrated the possibility of detecting mass wasting processes on crater walls.

Estimation procedure
The estimation of the Chebyshev polynomial coefficients is usually carried out on the basis of Least-squares (LS) principle (see Teunissen, 2009).Sampled elevation data yi are considered as the observation to approximate and cast in vector y0.Eq. ( 2), which describes the relation between the interpolated function and the sampled coordinates x, are implemented as parametric linear function model: where c is the vector containing the unknown M+1 coefficients to estimate: and A the so called design matrix (size m×M, being m the number of points sampled on the crater cross-sectional profile): Even though the observations are usually given the same weight in the LS estimate, a weight matrix W (size m×m) may be also introduced for weighting specific points in the profile.For example, in the case the precision of elevation points might be uneven, some weights could be used to differentiate points on steep crater slopes, on rim edge, on the extremes of the cross-sectional profile, for instance.In the case all observations are expected to feature the same precision, the weight matrix becomes W=Im×m.
The standard formula for LS solution allows estimating the unknown coefficient vector c: After the estimate of the solution, the sigma naught can be worked out: where residuals are given as: A 2 statistical test on the estimated (see, e.g., Teunissen, 2006) is normally applied to check out the conformity of the adopted functional model to real data.Given the a priori precision of elevation points as , the following relationship is used for testing: If the 2 test on the estimated fails, three reasons have to be considered and discussed: 1.A non-appropriate number of coefficients has been considered, i.e., the adopted functional model is not adequate; 2. The a priori precision ( ) of elevation data in the crater profile is not correctly guessed; and 3.The presence of anomalies in observations has to be investigated, since these may be due to gross measurement errors or unmodelled systematic errors).To this purpose, either a visual analysis of residuals plotted over the estimated approximating polynomial, and a statistical test on the standardized residuals are operated (see, e.g., Teunissen, 2006).
A second analysis that is applied concerns the statistical significance of the estimated Chebyshev's coefficients.This property can be investigated after computing the covariance matrix of the estimated solution: which reports the estimated variances ( ) of polynomial coefficients on the main diagonal.Notice that the offdiagonal elements of the covariance matrix should be naught in such a case, since the coefficients are supposed to be uncorrelated among them.
A t Student test on each estimated coefficient is applied to check out whether it is significantly different from Ci=0: In the case a coefficient i has not passed this test, this is constrained to be Ci=0.Because of the orthogonality of the Chebyshev basis functions, to constrain one or more coefficients is equivalent to put them as zero in the final polynomial pM(x) that will be implemented for the approximation of the cross-sectional profile.

FIRST EXPERIMENTAL RESULTS
The procedure for estimating the Chebyshev approximation of a lunar crater cross-sectional profile has been tested on the Cassini A crater.In next subsection 4.1 a description of this crater will be given.The choice has been made because of the particular shape of Cassini A, which is visible in Figure 2. A large slope failure invaded the Eastern lobe of the crater, resulting in an elongation of the original circular shape after the impact.Consequently, this crater is a perfect case study to prove algorithms for landslide recognition and analysis.

Cassini A crater
Cassini A (50.5° N, 4.8° E) is an impact crater located in Palus Nebularum, in the Eastern part of Mare Imbrium.It belongs to a family of 15 craters (Fig. 3), all entitled to Cassini but distinguished with different letters.Cassini A, whose diameter ranges approx.15 km, is contained in the bigger Cassini crater (average diameter 57 km).The internal area of Cassini crater was filled by lava that created a flat surface.
As already mentioned at the beginning of this Section, this crater has been chosen due to the presence of a large landslide inside.This slope failure can be easily detected with a visual interpretation of optical images, as reported in Brunetti et al. (2015).On the other hand, if the crosssectional profile W-E passes through the landslide body, the cross-section S-N is much less altered.Consequently, Cassini A also allows the comparison between crosssectional profiles with and without landslides.

Data sets
The digital elevation model (DEM) adopted for the extraction of cross-sectional profile is WACGLD100 (Scholten et al., 2012).This DEM was obtained from the photogrammetric process of images recorded by LROC Wide Angle Camera -WAC (Robinson et al., 2010).The resolution of this DEM is 100 m x100 m, while the average accuracy of elevation is better than ±20 m.In lunar nearside maria, this accuracy is even better than ±10 m.In order to analyze the surface texture, images from LROC Narrow Angle Camera -NAC have been used.The GSD of NAC images may reach 0.5 m.In the future development of this research, the use of Chang'E-2 images recorded by the CCD camera onboard is devised, considering a GSD between 1.5 and 7 m.

Cross-sectional profile approximation process
Four cross-sectional profiles have been interpolated from WACGLD100 DEM using bilinear interpolation.As it can be seen in Figure 2, all these profiles pass through the same points and are aligned along directions: W-E, S-N, SW-NE, NW-SE.The length of each section has been extended beyond both rims of a quantity that is approximately 30% of the rim-to-rim distance.This length is shorter than the one adopted in Mahanti et al. (2014), who used profiles twice the rim-to-rim distance to cope with a large group of craters.The solution implemented here is motivated by the need of tailoring the approximation to model the profiles of Cassini A, rather than defining a standard profile length to be applied in general.Resulting cross-sectional profiles feature a spatial resolution steps of 200 m, a total length (the same for all directions) of 25 km, and a total number of 127 elevation points.Considering that the interpolation from the grid DEM to the cross-sectional profiles may degrade the accuracy of elevation points, we assumed an accuracy of ±10 m for them.
The approximation with Chebyshev polynomials was operated by considering M=26 terms, corresponding to a maximum degree of 25.Also in this case, the number of coefficients suggested by Mahanti et al. (2014), i.e., M=16, has been changed in order to better model the shape of Cassini A crater.Furthermore, testing the statistical significance of the higher order coefficients was a point be focused in this research.According to this concept, the crosssectional profiles shown in Figure 5 have been obtained.A first analysis concerned the statistical testing of the results.In all four section, the test (9) on the estimated variance of unary weight ( ) failed.This occurred also when the a priori precision has been assumed to be ±20 m.Then a second test has concerned the statistical significance of the estimated Chebishev coefficients (11).In such a case, In the W-E profile, three coefficients have been found to be non significant.Among these, only one is included in the subset with M 16 (i.e., the one proposed in Mahanti et al., 2014).In other profiles, a number between 8-9 coefficients has not passed test (11).Among these rejected coefficients, 3 of them feature M 16 (results are shown in Table 4).By removing the non-significant coefficients Ci from the final polynomial adopted for the approximation of the crosssectional profiles, an increase of the estimated has been found, as expected due to the higher truncation error.On the other hand, the RMS (Root Mean Square) of the estimated for the four profiles before and after removing nonsignificant Ci increased from 27.9 m to 30.1 m (+2.2 m, corresponding to 7.8% of the initial RMS).This difference is relatively small when compared with the assumed a priori precision of the elevation data.Thus, removing these coefficients does not influence so much the quality of the approximation process.On the other hand, in only one case a low-order Ci (M 4) has not passed test (11).This has occurred for coefficient C3 of S-N profile.This result is particularly meaningful since Mahanti et al. (2014) relates this coefficient to the lack of asymmetry in the profile.This property can be confirmed by looking at the S-N profile reported in Figure 5. Next step of the analysis has been focused on the standard residuals after LS estimate, which are expected to be normally distributed with zero mean.As reported in Table 4, the fraction of points featuring large standard residual (larger than 2 and 3, respectively), resulted to be quite low to address the presence of a significant fraction of gross errors.On the other hand, the analysis of conformity of such standard residuals w.r.to the standard normal distribution has highlighted the presence of some unmodelled systematic errors.By looking at the profiles in Figure 5, this effect may be related to the presence of local short-length dunes (less than 1 km of horizontal extension) in the residuals.Red circles point out where these dunes are located.On the other hand, the frequency of these dunes is shorter that the critical Nyquist frequency, so that they are not due to ambiguous sampling.Further investigation about this outcome are reported in subsection 4.4.Table 4. -Results of some statistical analysis on the results of approximation using Chebishev polynomials.

Analysis of odd Chebishev coefficients
As already addressed in subsection 3.3, the presence of large odd coefficients is correlated to the general asymmetry of the cross-sectional profile.Asymmetry can be then linked to mass movement from the slopes to the bottom of the crater, as in the case of Cassini A crater.To investigate this property, the absolute value of standardized odd coefficients has been compared for all four profiles: The use of this quantity allows one to keep into consideration also the uncertainty of the estimated coefficient.On the other hand, neglecting the sign does not result in the loss of information for the purpose of seeking for the asymmetry.The standardized odd coefficients for the four cross-sectional profiles of Cassini A are plotted in Figure 6.Here having i>13 have been omitted because of their small size.As it can be seen by comparing Figures 5 and 6, two extreme cases can be pinpointed.The largest values of can be found in correspondence of the W-E profile that is clearly asymmetric.Smallest values are related to the S-N profile, featuring a high degree of symmetry.Indeed, in Figure 2 this profile does not intersect the landslide body.On the contrary, both diagonal profiles do that, and this fact can be recognized in the plot of absolute standardized coefficients.In the case of SW-NE profile, a very large value for is reported.In the case of NW-SE profile, the coefficients with odd order between 3 and 9 show values that are higher than the average value.

Analysis of short-length dunes in crater profiles
As described in subsection 4.3, some short-length dunes have been found in the residuals between cross-sectional profile data and values in corresponding positions as predicted using estimated Chebyshev polynomials.First of all, two possible causative reasons of this result have been excluded.The first reason entails some artifacts possibly coming out during the interpolation from WACGLD100 DEM to the cross-sectional profiles.In fact, the lower spatial resolution of the latter (1 point/200m) with respect to the former (100 m x 100 m) prevents from such an effect.The second could be related to some systematic errors occurred during the photogrammetric processing of WAC imagery.On the other hand, there is not any spatial correlations between these observed dunes.It is then quite evident that these effects are probably due to the local topography, which could not be correctly modelled using Chebyshev polynomials.Two other factors have been considered.First, the excursion from the lowest to the highest points of each dune ranges in the order of a few tens of metres up to approximately 120 m.In a large fraction of cases, the local dunes feature a lack of mass in the upper part, and an almost equivalent accumulation area in the lower part.The interpretation of these phenomena could be then related to local small failures that mainly occurred in higher slopes, and outside the main accumulation body of the big landslide that can be recognized in the W-E direction.A last attempt to better understand this process has been done by looking at LROC-NAC (Narrow Angle Camera) images over corresponding areas.This analysis has not been accomplished yet in exhaustive manner.Three examples are shown in Figure 7, demonstrating that in the areas where short-length dunes have been found, images show different kinds of ground texture.These differences might be related to local mass movements.

CONCLUSIONS AND FUTURE WORK
This paper presented the application of Chebyshev polynomial for approximation lunar crater cross-sectional profiles and extracting some property (i.e., asymmetry) that could be related to mass wasting processes.The analysis of asymmetry in the cross-sectional profiles has confirmed to be a useful to detect large mass wasting processes.Here a method based on the analysis of absolute standardized odd coefficients has been successfully applied.On the other hand, some criteria to establish selective threshold have to be developed in order to make the discrimination of asymmetry more objective.Also the application of other statistical indicators of asymmetry (like skewness, see Crosilla et al, 2013) are worth to be applied.In the future, data from the Chinese mission Chang'E 2 will be also considered due to their similar spatial resolution (Zhao et al., 2011).Further studies will include: the measurement of the landslide volume; the analysis of relationships between landslides and characteristics of the hosting craters as well as the surrounding terrain; the lithological and mineralogical characterization of surfaces using multispectral data acquired by the sensor IIM (Imaging Interferometer Spectrometer) mounted on Chang'E 1 (Sun et al., 2005), comparing the data with available spectral libraries.

Figure 1 .
Figure 1.-Example of a landslide mapped in Cassini A crater (see Subsect. 4.1).The blue circle approximates the crater rim; purple and green shaded areas are the landslide scarp and deposit, respectively.(Credit: NASA/Goddard Space Flight Center/Arizona State University).

Figure 2 .
Figure 2. -Mosaic of tiles from LROC-NAC camera showing the Cassini A crater. Green dotted lines represent the position of cross-sectional profiles that have been analyzed using Chebyshev polynomials.Crosses point out the position of those short-length dunes that can be recognized in the residuals after approximation (see Subsect. 4.4).

Figure 3 .
Figure 3. -Mosaic of tiles from LROC-WAC camera showing the area with the most Cassini craters.

Figure 5
Figure 5. -Cross-sectional profiles of Cassini A crater; red circles highlight some short-length dunes found in residuals.

Figure 6 .
Figure 6.-Results of odd standardized coefficient ( ) analysis for cross-sectional profiles of Cassini A crater.

Figure 7 .
Figure 7. -Three portions of rectified LROC-NAC images showing some details in correspondence of short-length dunes detected in profile of Cassini A crater.On the left: W-E profile, area denoted as 'High left' in Fig. 5; on the centre: W-E profile, area denoted as 'High right'; on the right: N-S profile, area denoted as 'High right