Automated Discontinuity Detection and Reconstruction in Subsurface Environment of Mars Using Deep Learning: A Case Study of SHARAD Observation

: Machine learning (ML) algorithmic developments and improvements in Earth and planetary science are expected to bring enormous beneﬁts for areas such as geospatial database construction, automated geological feature reconstruction, and surface dating. In this study, we aim to develop a deep learning (DL) approach to reconstruct the subsurface discontinuities in the subsurface environment of Mars employing the echoes of the Shallow Subsurface Radar (SHARAD), a sounding radar equipped on the Mars Reconnaissance Orbiter (MRO). Although SHARAD has produced highly valuable information about the Martian subsurface, the interpretation of the radar echo of SHARAD is a challenging task considering the vast stocks of datasets and the noisy signal. Therefore, we introduced a 3D subsurface mapping strategy consisting of radar echo pre-processors and a DL algorithm to automatically detect subsurface discontinuities. The developed components the of DL algorithm were synthesized into a subsurface mapping scheme and applied over a few target areas such as mid-latitude lobate debris aprons (LDAs), polar deposits and shallow icy bodies around the Phoenix landing site. The outcomes of the subsurface discontinuity detection scheme were rigorously validated by computing several quality metrics such as accuracy, recall, Jaccard index, etc. In the context of undergoing development and its output, we expect to automatically trace the shapes of Martian subsurface icy structures with further improvements in the DL algorithm.


Introduction
Conventional computer algorithms in the geoscience field have been applied mainly for spatial data processing such as stereo analysis of topography, multispectral/hyperspectral analysis, and other remote sensing data compiling. These days, novel approaches combined with contemporary machine learning algorithms are changing the trend of computational geoscience. This is because several ground and remote sensors are producing gigantic geoscientific datasets, which cannot be managed effectively by manual interpretation in terms of precision and speed. The active dynamics of the terrestrial surface provide an optimal ground, where a massive data collection scheme and machine leaning play a significant role. In [1,2], a digital terrain model (DTM) of a wide area with very high resolution was constructed using an advanced stereo machine vision algorithm on the Advanced Land Observing Satellite (ALOS) Panchromatic Remote-Sensing Instrument for stereo mapping (PRISM) and WorldView satellite data respectively. research, employing MARSIS and SHARAD, is the methodology to effectively trace the discontinuities from subsurface signals.
Thus, in this study, we focused on the development of an automated algorithm for subsurface discontinuity detection using DL. The detected discontinuities were later utilized for subsurface reconstruction. Considering the high signal-to-noise ratio of SHARAD, we actively employed a deep-learning approach: specifically, a convolutional neural network along with add-on processors for noise reduction, registration between clutter simulation and radargram, and feature extraction.
The constructed algorithms were tested on Martian geomorphic features such as Lobate Debris Apron (LDA), Polar Layered Deposits, and mid-latitude potential ice layers beneath the Phoenix landing site, as described in Section 2. The technical details of the algorithms are stated in Section 3. The result and validation outcomes employing error metrics are presented in Sections 4 and 5, along with relevant discussion.

Target Area
The first target area to test the discontinuity detection is LDA, which is regarded as an evidence of potential glacial processes [28] or ice debris transport in the ancient Martian environment [29]. It can be characterized as an alluvial-fan-like mass of debris, whose thickness varies between tens to hundreds of meters, that flowed from the mountain slope covered by a glacier or periglacial landforms roughly some million years ago [30]. Despite ongoing arguments regarding the formation mechanism of LDA [28][29][30][31], its geomorphic origin is still unclear. Therefore, the subsurface investigation of LDA employing SHARAD data is expected to provide useful clues for LDA structures. It is worth noting that LDAs are mostly distributed along with the dichotomy of the southern hemisphere and northern lowlands [32]. Therefore, the hypothesis that such highland features were presumably involved with glacial and periglacial processes has gained attention. Thus, the scenario established on this clear evidence is that climate change induced the melting of the glaciers, which allowed sediment flow down through the steep sides of LDAs and produced their fan-shaped features (see Figure 1c,d).
In this study, two sites, i.e., Euripus Mons LDA & an anonymous LDA (44.8 • S 105.2 • E) close to Euripus Mons, were chosen as the target for SHARAD subsurface feature detection. Based on previous studies, high ice content has already been identified [28] in those target areas. It has also been proposed that at least some of them might be composed of pure ice and, protected from sublimation by a thin debris cover. Recent modeling practices, assuming glacial deformation, have helped to reconstruct some of the inherent structural properties [33]. Moreover, similar modeling studies employing High Resolution Stereo Camera (HRSC) DEM and non-linear flow calculations have shown a potential involvement of dirty ice [34]. Considering the undergoing arguments over these target areas, the identification of subsurface structures by a deep-learning approach may imply essential clues about the origin of LDA and its deposits. Therefore, we choose two LDAs as our first target area to deploy the developed algorithms. In this study, two sites, i.e., Euripus Mons LDA & an anonymous LDA (44.8° 105.2° ) close to Euripus Mons, were chosen as the target for SHARAD subsurface feature detection. Based on previous studies, high ice content has already been identified [28] in those target areas. It has also been proposed that at least some of them might be composed of pure ice and, protected from sublimation by a thin debris cover. Recent modeling practices, assuming glacial deformation, have helped to reconstruct some of the inherent structural properties [33]. Moreover, similar modeling studies employing High Resolution Stereo Camera (HRSC) DEM and non-linear flow calculations have shown a potential involvement of dirty ice [34]. Considering the undergoing arguments over these target areas, the identification of subsurface structures by a deep-learning approach may imply essential clues about the origin of LDA and its deposits. Therefore, we choose two LDAs as our first target area to deploy the developed algorithms.
Possibly, the Martian deposits in the southern polar zone represent the most interesting history of the evolution of Mars as a solid planet. Several researchers have proposed that the multiple layers of the southern Martian pole represent a sequence of episodic climate changes [35,36]. Therefore, the Martian south polar layered deposits (SPLD) presumably have involvements with the other footprints of episodic environmental changes, as shown in the potential coastal line of the northern hemisphere [37,38], layered fluvial deposition [39,40] and a recently claimed underground lake [41]. The investigation of the subsurface structure of SPLD was proposed as a prime target of SHARAD and MARSIS campaigns, which revealed expected multiple subsurface layering as described in [42]. However, a high density of discontinuity points in SPLD induced severe difficulties in correctly inferring the detailed structures and locations of subsurface layers. Therefore, we choose a part of the SPLD structure as a test area for the constructed algorithm ( Figure 1c). The most challenging target area for the application of the developed algorithms is the Phoenix landing site (Figure 1d), where Possibly, the Martian deposits in the southern polar zone represent the most interesting history of the evolution of Mars as a solid planet. Several researchers have proposed that the multiple layers of the southern Martian pole represent a sequence of episodic climate changes [35,36]. Therefore, the Martian south polar layered deposits (SPLD) presumably have involvements with the other footprints of episodic environmental changes, as shown in the potential coastal line of the northern hemisphere [37,38], layered fluvial deposition [39,40] and a recently claimed underground lake [41]. The investigation of the subsurface structure of SPLD was proposed as a prime target of SHARAD and MARSIS campaigns, which revealed expected multiple subsurface layering as described in [42]. However, a high density of discontinuity points in SPLD induced severe difficulties in correctly inferring the detailed structures and locations of subsurface layers. Therefore, we choose a part of the SPLD structure as a test area for the constructed algorithm ( Figure 1c). The most challenging target area for the application of the developed algorithms is the Phoenix landing site (Figure 1d), where weak discontinuity caused by potential icy bodies are expected. There have been assumptions that the layer of permafrost is extended to the mid-latitude region of Mars [43], which were later confirmed from the excavation during the Phoenix mission [44]. However, such ice layers close to the surface were not the target of SHARAD and MARSIS. Although SHARAD data analysis on the Phoenix landing site identified the discontinuity (which may be an extension of the ice layer very close to the surface), the discontinuity signals are very weak. Due to this difficulty, the application to SHARAD data over the Phoenix landing site should be tested to prove the robustness of developed algorithms. [19], is the first spaceborne ground penetration RADAR (GPR) to map the upper portions of the crust of a solid planet or satellite with an extremely low-frequency electromagnetic wave. Although the demand to investigate the subsurface liquid and solid water in the footprint of Martian geological evaluation has been actively proposed [45], the observation priority of MARSIS was given on the Martian polar deposits, as they may play a crucial role in atmospheric circulation and evolution [46,47]. In contrast, the SHARAD instrument of the Mars Reconnaissance Orbiter, which was developed through collaboration with the Italian Space Agency and the Jet Propulsion Laboratory, was launched in 2005 and characterized by better vertical resolution but lower penetration capability. Due to such system specifications, the complex substructures of the Martian crust, such as polar deposits, pedestal craters [48], and LDA [28], have been successfully observed employing SHARAD. To understand the system characteristics and operational concepts of SHARAD, see Table 1 and Figure 2, respectively. Since this study aims to extract detailed information on subsurface layering utilizing ML algorithms, SHARAD became the prime data source rather than MARSIS. The potential application to MARSIS data with developed ML algorithms will be explored in a future study.

MARSIS, equipped in Mars Express
confirmed from the excavation during the Phoenix mission [44]. However, such ice layers close to the surface were not the target of SHARAD and MARSIS. Although SHARAD data analysis on the Phoenix landing site identified the discontinuity (which may be an extension of the ice layer very close to the surface), the discontinuity signals are very weak. Due to this difficulty, the application to SHARAD data over the Phoenix landing site should be tested to prove the robustness of developed algorithms. [19], is the first spaceborne ground penetration RADAR (GPR) to map the upper portions of the crust of a solid planet or satellite with an extremely lowfrequency electromagnetic wave. Although the demand to investigate the subsurface liquid and solid water in the footprint of Martian geological evaluation has been actively proposed [45], the observation priority of MARSIS was given on the Martian polar deposits, as they may play a crucial role in atmospheric circulation and evolution [46,47]. In contrast, the SHARAD instrument of the Mars Reconnaissance Orbiter, which was developed through collaboration with the Italian Space Agency and the Jet Propulsion Laboratory, was launched in 2005 and characterized by better vertical resolution but lower penetration capability. Due to such system specifications, the complex substructures of the Martian crust, such as polar deposits, pedestal craters [48], and LDA [28], have been successfully observed employing SHARAD. To understand the system characteristics and operational concepts of SHARAD, see Table 1 and Figure 2, respectively. Since this study aims to extract detailed information on subsurface layering utilizing ML algorithms, SHARAD became the prime data source rather than MARSIS. The potential application to MARSIS data with developed ML algorithms will be explored in a future study.

System Characteristics SHARAD Specification
Carrier frequency 20 MHz

System Characteristics SHARAD Specification
Carrier frequency 20 MHz Vertical resolution (if surface permittivity is 4.0) 7.5 m * Transmitter power 10 W Pulse length 85 µs Pulse repetition frequency 700/350 Hz Horizontal resolution (along track x cross track) 0.3-1 km × 3-6 km Estimated maximum penetration depth in the scenario of porous (30%) materials and underneath ice filled pore spaces below [20] 1500 m **

Introduction to SHARAD Data Products
SHARAD tracks on our target areas are identified (see Table 2) using the interfaces established on Planetary Data System (PDS, https://ode.rsl.wustl.edu/mars/indexproductsearch.aspx). The SHARAD tracks are footprints of the signals collected by the sensor. They consist of varying along-track lengths. The data obtained from the PDS comprises of a radargram image file, which is a real-valued image of the radar backscatter power, with the vertical and horizontal axis representing the time delay and along-track samples, respectively. The along-track samples have a spacing of 128 pixels per degree (nearly 460 m) and are associated with latitude-longitude pairs, spacecraft position, planetary radius, and other information contained in a geometry file. The round-trip delay samples in the vertical axis are uniform at 0.0375 microseconds per pixel. The coordinate system is planetocentric with longitude positive toward the east [49]. On-board pre-summing is usually performed with typical pre-summing factors like 4 or 8, to control data volume transmission to Earth [50]. The degree of pre-summing affects only the Doppler frequency bandwidth (in Hz) of the Doppler spectrum. Doppler processing narrows the along-track resolution by filtering out signals that return to the sensor with more than a specified Doppler shift. Echoes from areas perpendicular to the flight path at this band of Doppler components are characterized only by round-trip delay, so reflections from features far from the ground track can appear as apparent subsurface returns and are depicted as clutters [51].
Since the data processing algorithm is not feasible to implement into in-house S/W, the CO-SHARPS [52] processing boutique is employed for radar signal processing. We describe the CO-SHARPS processing parameters in Section 3.1.

Processing Method
The complete process, covering radargram preprocessing, discontinuity detection and validation in addition to subsequent analysis is shown in Figure 3. The flowchart can be summarized as follows: 1.
The SHARAD tracks for the target regions are identified by inspecting track locations/coverage.

2.
The processed and raw radargrams are obtained using the CO-SHARPS processing boutique.

3.
A bilateral filtering method is applied to the radargram to reduce the backscatter noise. 4.
The clutter-gram is co-registered with raw radargram using a Fourier-based phase correlation process.

5.
The labelled data is prepared by manually detecting discontinuities in the radargram. 6.
New features such as Euclidean norm, entropy, etc. are prepared for training the CNN model. 7.
The model is validated using various quality metrics. The discontinuity points are accumulated and processed into a subsurface discontinuity plane reconstruction algorithm.

5.
The labelled data is prepared by manually detecting discontinuities in the radargram. 6. New features such as Euclidean norm, entropy, etc. are prepared for training the CNN model. 7. The model is validated using various quality metrics. The discontinuity points are accumulated and processed into a subsurface discontinuity plane reconstruction algorithm. A detailed description of the methodology is provided in following sub-sections.

Pre-processing
After identifying the product IDs of tracks over proposed target areas, the CO-SHARPS processing boutique was employed for the purpose of clutter simulation. The clutter simulator available in the CO-SHARPS boutique is based on the JPL (Jet Propulsion Laboratory, Caltech) focussed processor and comprises products called QDA clutter-grams [53].
The processing parameters of the CO-SHARPS boutique are decided as follows. The QDA products have a fixed aperture (frames) of 4096. The pre-sum value, which varies the bandwidth or frame interval, is set as 16. A Hann window is selected to reduce sidelobes in the azimuth compression. The QDA products in this study are simulated using MOLA DEM with the ellipsoid as the equipotential reference surface, i.e., aeroid. The multi-look bandwidth is fixed to 15 Hz for the QDA products. The simulator also offers the depth-corrected radargram if a dielectric constant is specified during QDA processing. It is worth noting that the dielectric constant in the clutter simulator is assumed to be uniform over the entire subsurface for simplified processing ( = 3.15, in this study). The dielectric constant varies in both latitudinal and longitudinal direction according to the subsurface material properties and is an area of ongoing research [26,27]. However, a dielectric constant value of 3.15 has been used as a suitable approximation in past studies of the Euripus Mons region [54,55].
It is optional to adopt either the depth-converted radargram or the time-domain radargram. The process of depth conversion results in under-sampling of the signals and, hence, would lead to loss A detailed description of the methodology is provided in following sub-sections.

Pre-Processing
After identifying the product IDs of tracks over proposed target areas, the CO-SHARPS processing boutique was employed for the purpose of clutter simulation. The clutter simulator available in the CO-SHARPS boutique is based on the JPL (Jet Propulsion Laboratory, Caltech) focussed processor and comprises products called QDA clutter-grams [53].
The processing parameters of the CO-SHARPS boutique are decided as follows. The QDA products have a fixed aperture (frames) of 4096. The pre-sum value, which varies the bandwidth or frame interval, is set as 16. A Hann window is selected to reduce sidelobes in the azimuth compression. The QDA products in this study are simulated using MOLA DEM with the ellipsoid as the equipotential reference surface, i.e., aeroid. The multi-look bandwidth is fixed to 15 Hz for the QDA products. The simulator also offers the depth-corrected radargram if a dielectric constant is specified during QDA processing. It is worth noting that the dielectric constant in the clutter simulator is assumed to be uniform over the entire subsurface for simplified processing (ε = 3.15, in this study). The dielectric constant varies in both latitudinal and longitudinal direction according to the subsurface material properties and is an area of ongoing research [26,27]. However, a dielectric constant value of 3.15 has been used as a suitable approximation in past studies of the Euripus Mons region [54,55].
It is optional to adopt either the depth-converted radargram or the time-domain radargram. The process of depth conversion results in under-sampling of the signals and, hence, would lead to loss of useful information, potentially indicative of desired clutter. In this study, the time-domain radargram is utilized as it does not suppress the original signal. The selection of dielectric constant does not affect the DL algorithm in this study, as the model has been trained on time-domain radargram, and the depth conversion has been employed only for 3D-subsurface reconstruction.
Noise-reduction techniques are introduced to reduce the backscatter or background speckles in radargram images. It supports the identification of discontinuities in the data as the signal-to-noise ratio in the original radargram is not adequate. Several filtering methods such as log-Gabor filtering, block-matching, and wavelet shrinkage denoising, etc. [56] were applied for suppressing the backscatter noise. Finally, bilateral filtering was employed (refer Figure 4) as it provides a balanced outcome between the preservation of the original signals and the noise levels.
Noise-reduction techniques are introduced to reduce the backscatter or background speckles in radargram images. It supports the identification of discontinuities in the data as the signal-to-noise ratio in the original radargram is not adequate. Several filtering methods such as log-Gabor filtering, block-matching, and wavelet shrinkage denoising, etc. [56] were applied for suppressing the backscatter noise. Finally, bilateral filtering was employed (refer Figure 4) as it provides a balanced outcome between the preservation of the original signals and the noise levels.

Feature Space Construction
The clutter-simulated radargram is co-registered with the original radargram to remove vertical and horizontal offsets. Due to the nature of radargram signals, translation, rotation, and scaleinvariant image registration are essential. Therefore, we employed a fast Fourier transform (FFT)based registration method, which utilizes frequency domain analysis to find matching signals [57]. The performance of the FFT-based registration method is demonstrated in Figure 5.

Feature Space Construction
The clutter-simulated radargram is co-registered with the original radargram to remove vertical and horizontal offsets. Due to the nature of radargram signals, translation, rotation, and scale-invariant image registration are essential. Therefore, we employed a fast Fourier transform (FFT)-based registration method, which utilizes frequency domain analysis to find matching signals [57]. The performance of the FFT-based registration method is demonstrated in Figure 5. The labelled discontinuity is prepared manually for approximately twenty tracks from different regions for the training components of the machine learning algorithm (see Figure 6) The characteristics of "subsurface discontinuity" do not vary significantly for different regions, as per our experiences. Therefore, our deep learning model does not employ regionally selected training vectors. The construction of relevant feature spaces is highly essential for the performance of the ML algorithm. Several features are generated and tested with approaches such as normalized difference, mean, Gabor filter and image entropy, etc. Finally, a three-channel feature space is constructed using the radargram, registered clutter-gram, and the normalized difference of the radargram and cluttergram (refer to Figure 6). In order to enhance the robustness of the model and increase the number of The labelled discontinuity is prepared manually for approximately twenty tracks from different regions for the training components of the machine learning algorithm (see Figure 6) The characteristics of "subsurface discontinuity" do not vary significantly for different regions, as per our experiences. Therefore, our deep learning model does not employ regionally selected training Appl. Sci. 2020, 10, 2279 9 of 18 vectors. The construction of relevant feature spaces is highly essential for the performance of the ML algorithm. Several features are generated and tested with approaches such as normalized difference, mean, Gabor filter and image entropy, etc. Finally, a three-channel feature space is constructed using the radargram, registered clutter-gram, and the normalized difference of the radargram and clutter-gram (refer to Figure 6). In order to enhance the robustness of the model and increase the number of training samples, an additional training data set was created using varied orientation and affine transforms on the original dataset. The augmentation was performed by rotating the images through a range of −5 to +5 degrees, mirror/left-right flip, resampling and image zooming. The overall training dataset consisted of 2160 images of size 256 × 256. The labelled discontinuity is prepared manually for approximately twenty tracks from different regions for the training components of the machine learning algorithm (see Figure 6) The characteristics of "subsurface discontinuity" do not vary significantly for different regions, as per our experiences. Therefore, our deep learning model does not employ regionally selected training vectors. The construction of relevant feature spaces is highly essential for the performance of the ML algorithm. Several features are generated and tested with approaches such as normalized difference, mean, Gabor filter and image entropy, etc. Finally, a three-channel feature space is constructed using the radargram, registered clutter-gram, and the normalized difference of the radargram and cluttergram (refer to Figure 6). In order to enhance the robustness of the model and increase the number of training samples, an additional training data set was created using varied orientation and affine transforms on the original dataset. The augmentation was performed by rotating the images through a range of −5 to +5 degrees, mirror/left-right flip, resampling and image zooming. The overall training dataset consisted of 2160 images of size 256 × 256.

Machine Vision Algorithm
The convolutional neural network has been used for training of the model. The trained network can be divided into two parts (i.e., the down-sampling or contracting path and the up-sampling or expanding path) [58]. The architecture of the CNN model used in this research can be represented by Figure 7.
The down-sampling or contracting path is composed of blocks, which are comprised of 3 × 3 Convolution layer and a 2 × 2 Max Pooling layer with an activation function at each layer [50]. The input layer is normalized by adjusting and scaling the activations or batch normalization. The batch normalization improves the speed, performance, and stability of neural networks [59]. The trained deep learning model in this research draws its strength from making normalization a part of the model architecture and performing the normalization for each training mini batch. The network uses max-pooling layers, which reduce the height and width information while keeping the number of channels of the input matrix constant. These layers are intended to increase the resolution of the output. Pooling layers can work with different methods, including maximum, average, or median layers. The contracting path captures the context of the input image in order to perform segmentation. This coarse contextual information is transferred to the expanding path using skip connections.
The expanding path is composed of the block with 3 × 3 convolution layers, deconvolution layers and a concatenation operation on the cropped feature map from the contracting path. The purpose of expanding the path is to enable precise localization. The workflow for training and testing the image dataset is shown in Figure 8.

Machine Vision Algorithm
The convolutional neural network has been used for training of the model. The trained network can be divided into two parts (i.e., the down-sampling or contracting path and the up-sampling or expanding path) [58]. The architecture of the CNN model used in this research can be represented by Figure 7. The down-sampling or contracting path is composed of blocks, which are comprised of 3 × 3 Convolution layer and a 2 × 2 Max Pooling layer with an activation function at each layer [50]. The input layer is normalized by adjusting and scaling the activations or batch normalization. The batch normalization improves the speed, performance, and stability of neural networks [59]. The trained deep learning model in this research draws its strength from making normalization a part of the model architecture and performing the normalization for each training mini batch. The network uses max-pooling layers, which reduce the height and width information while keeping the number of channels of the input matrix constant. These layers are intended to increase the resolution of the output. Pooling layers can work with different methods, including maximum, average, or median layers. The contracting path captures the context of the input image in order to perform segmentation. This coarse contextual information is transferred to the expanding path using skip connections. The expanding path is composed of the block with 3x3 convolution layers, deconvolution layers and a concatenation operation on the cropped feature map from the contracting path. The purpose of expanding the path is to enable precise localization. The workflow for training and testing the image dataset is shown in Figure 8.

Evaluation Metric
Several evaluation metrics such as accuracy, precision, F-score, etc. have been used to evaluate the performance of the model. The accuracy depicts the ratio of the correctly classified samples to the total number of samples. Precision helps to evaluate how well the model classifies the relevant class with respect to the total number of predicted instances of that class. A similar metric called

Evaluation Metric
Several evaluation metrics such as accuracy, precision, F-score, etc. have been used to evaluate the performance of the model. The accuracy depicts the ratio of the correctly classified samples to the total number of samples. Precision helps to evaluate how well the model classifies the relevant class with respect to the total number of predicted instances of that class. A similar metric called "specificity" is also used.
Precision (P) = True Positive/(True Positive + False Positives) Speci f icity = True Negative/(True Negative + False Positives) Recall (also known as sensitivity) denotes the fraction of correctly predicted relevant instances (class y = 1) in the actual number of relevant instances. This gives us a measure of how accurately the model predicts with respect to the actual instances of that class.
These two terms, i.e., precision and recall, give a good estimate of the model. The ideal value for both is close to 1. To combine and evaluate these together, we use the metric F-score or F-measure. A high F-score is an indication of a good model.
The area under curve (AUC) of the receiver operating characteristics (ROC) curve is a statistical value that denotes the ability of a classifier to perform as its discrimination threshold (probability) is varied. An ideal AUC value is close to 1.
The Jaccard Index or IoU (Intersection over Union) is a metric is used to measure the similarities or overlap between predicted and actual sample sets A and B, respectively. It is defined as the cardinality or size of the intersection divided by the size of the union of the two sets. A value close to 1 is considered to be ideal, i.e., when both the sets are equal.

Results and Discussion
The training was performed on the tracks over Euripus Mons, which also contained the portion beyond the geographic extent of Euripus mons. Several new tracks, with almost no similarity to the other tracks in the region, were processed for testing. The results obtained are shown in Figure 9. Tables 3 and 4 show the confusion matrix and performance metrics over the respective target areas. The testing performed on the Euripus Mons tracks illustrate the best overall results as expected. It is important to note that a substantial proportion of the total labelled discontinuities belonged to the Euripus Mons region when compared to the Anonymous LDA and Polar region (refer to columns 1 and 3 of Table 3). The testing on the Anonymous LDA region also yielded fair results considering the true class size ratio. The model was able to predict 22 of the total 81 discontinuity pixels in the area. A few tracks in the polar region were also processed for testing. It can be inferred from Table 3 that there was no true class in the actual data for this region. However, the model predicted a small number of false positives. Despite the vast size of the Polar region, not many tracks which cover the Euripus Mons region also pass through it (refer Table 2).
The model, when tested on tracks from the Phoenix landing site, did not predict any pixel as discontinuity even though there were 67 manually detected true class pixels. While preparing the labelled dataset for this region, some tracks were omitted, as there was no/very little discontinuity in the region. Since the Phoenix landing site does not contain significant topographical variation, the discontinuities are not considerable. The test results were not very accurate, as expected in this region, indicating that geographically more distributed labelled data is required for training, to enhance the robustness of the model. However, it is important to note that the manually detected labelled data had only 1 of 3000 pixels as discontinuity pixels.
Moreover, the Phoenix landing site lies in the northern lowland of Mars, and the training dataset belongs to the southern hemisphere, which has a different geological origin. Due to significant variation in the subsurface strata of Euripus Mons, a southern hemispheric highland, and the Phoenix landing site, a primary northern lowland of Mars, the training may not be spatially robust. To overcome this problem, and improve the model accuracy, more training data must be prepared covering a wide array of the geographical locations of Mars. The model was trained for 20 epochs with a batch size of 16, making it approximately 135 iterations. The epoch refers to a complete cycle of training, i.e., one forward pass and one backward pass of all the training samples. The training and validation loss as well as accuracy plotted at each epoch are shown in Figure 10. The training loss is initially high and decreases monotonically in subsequent epochs, demonstrating that the model did not present a bias (or no-skill) and was trained adequately. The accuracy increases with the number of epochs. Although some undulations occur in training accuracy, however, they may be attributed to a random sampling of training data in each epoch.    It is essential to assess the outcomes of the model in the geological context, as well as in terms of quantitative metrics for further scientific applications. The efforts towards geological applications were accomplished using discontinuity plane reconstruction, together with contextual interpretation.
After collecting all discontinuity points from the sublayers, the points were clustered using the Balanced iterative reducing and clustering using hierarchies (BIRCH) algorithm [60] into a few initial groups according to their spatial distributions. Kriging was applied to extract interpolated sublayer planes and the Euclidian distances between discontinuity points and interpolated planes. The best fit underground sublayers were extracted using an iterative procedure to attain an optimal number of clusters and distances. In this research, manual intervention was required to accomplish this task; however, the development of an automated/semi-automated routine is under development. The above procedure to reconstruct the underlying discontinuity planes revealed the gridded subsurface structures of Euripus Mons and anonymous LDAs, as shown in Figure 11. We discovered three clear discontinuity sublayers underneath the Euripus Mons LDA, represented by the yellow, green, and blue face. It is essential to assess the outcomes of the model in the geological context, as well as in terms of quantitative metrics for further scientific applications. The efforts towards geological applications were accomplished using discontinuity plane reconstruction, together with contextual interpretation.
After collecting all discontinuity points from the sublayers, the points were clustered using the Balanced iterative reducing and clustering using hierarchies (BIRCH) algorithm [60] into a few initial groups according to their spatial distributions. Kriging was applied to extract interpolated sublayer planes and the Euclidian distances between discontinuity points and interpolated planes. The best fit underground sublayers were extracted using an iterative procedure to attain an optimal number of clusters and distances. In this research, manual intervention was required to accomplish this task; however, the development of an automated/semi-automated routine is under development. The above procedure to reconstruct the underlying discontinuity planes revealed the gridded subsurface structures of Euripus Mons and anonymous LDAs, as shown in Figure 11. We discovered three clear discontinuity sublayers underneath the Euripus Mons LDA, represented by the yellow, green, and blue face.   We also identified two discontinuity layers under anonymous LDA, as shown in Figure 11, and ascertained different morphologies in northern and southern aprons. We propose that such different morphologies of the two subsurface layers imply a certain temporal/contextual disconnection between paleoclimate events and progressive developments during intermediate periods, perhaps involving sublimation. This is a significant clue towards identifying the morphological evolution of LDAs. We did not observe three discontinuity layers, which follow the shape of surface topography, as claimed in a previous study [61]. This may be due to sparsity in the reconstructed discontinuity points. Moreover, the integration of several studies [25,26] for the estimation of dielectric constants will be considered for future discontinuity layer reconstruction algorithms. The depth of discontinuity layers will be adjusted to precisely depict subsurface structures using varying dielectric constants.
Finally, the assessment of our automated discontinuity detection approach guarantees sufficient accuracy to replace manual detection, given that proper training and pre-processing are preceded. Moreover, the reconstructed discontinuity planes over LDAs and the automated MV algorithms on SHARAD signals can provide meaningful geological information. The latter implies the potential of our method to mine significant clues regarding Martian evolution. We also identified two discontinuity layers under anonymous LDA, as shown in Figure 11, and ascertained different morphologies in northern and southern aprons. We propose that such different morphologies of the two subsurface layers imply a certain temporal/contextual disconnection between paleoclimate events and progressive developments during intermediate periods, perhaps involving sublimation. This is a significant clue towards identifying the morphological evolution of LDAs. We did not observe three discontinuity layers, which follow the shape of surface topography, as claimed in a previous study [61]. This may be due to sparsity in the reconstructed discontinuity points. Moreover, the integration of several studies [25,26] for the estimation of dielectric constants will be considered for future discontinuity layer reconstruction algorithms. The depth of discontinuity layers will be adjusted to precisely depict subsurface structures using varying dielectric constants.

Conclusion and Future Work
Finally, the assessment of our automated discontinuity detection approach guarantees sufficient accuracy to replace manual detection, given that proper training and pre-processing are preceded. Moreover, the reconstructed discontinuity planes over LDAs and the automated MV algorithms on SHARAD signals can provide meaningful geological information. The latter implies the potential of our method to mine significant clues regarding Martian evolution.

Conclusions and Future Work
Although the data analysis of orbital radars such as MARSIS and SHARAD are producing many compelling clues about Martian surface and climate evolution, their interpretation still depends on 1D profiling of signals, which are created manually by experienced interpreters. This manual profiling has reduced the potential for synthesized use of SHARAD and MARSIS data sets with other spatial information to extract scientific clues about the Martian surface and subsurface. Therefore, we have developed an automated 3D subsurface mapping procedure that enables us to automatically define surface discontinuity and reconstruct subsurface layers from SHARAD signals using ML/DL. We have applied the established scheme over several target areas such as LDAs in the Euripus Mons region, high altitude potential permafrost, and polar deposit. The developed method automatically identified several discontinuity layers with adequate accuracy. The detected discontinuity points were utilized to reconstruct subsurface discontinuity planes.
In comparison to the precedent studies aimed at bringing partially automated improvement for the analysis of SHARAD and MARSIS data sets such as point detection, auto registration, filtering, and reconstruction, our method depicts further potential to be developed as a comprehensive SHARAD and MARSIS procedure with improved accuracy. We are planning to improve the established algorithms as follows.
(a) Employment of HRSC & Context Camera (CTX) stereo DTMs for better clutter simulation. (b) Automated subsurface layer clustering based on more robust approaches, such as support vector machine (SVM) algorithms. (c) ML algorithm improvement, by preparing more labelled data from spatially varied regions to increase the robustness of the model.
We will publish these improvements in the public domain after the fully stabilized components are integrated with the MARSIS/SHARAD signal processor. However, for robust training, SHARAD tracks and discontinuity labels from several other areas will be required based on user feedback for the published software. The case will not be limited only to Mars, and may also apply to future planetary subsurface radars such as Radar for Europa Assessment and Sounding: Ocean to Near-surface (REASON) [62].
Author Contributions: V.G. was charged of major data processing and writing. S.K.G. developed ML algorithms and wrote corresponding text. J.K. designed the overall strategy and structured the draft. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.