A New Method for Automatically Tracing Englacial Layers from MCoRDS Data in NW Greenland †

Englacial layering reflects ice dynamics within the ice bodies, which improves understanding of ice flow variation, past accumulation rates and vertical flows transferring between the surface and the underlying bedrock. The internal layers can be observed by using Radar Echo Sounding (RES), such as the Multi-channel Coherent Radar Depth Sounder (MCoRDS) used in NASA’s Operation IceBridge (OIB) mission. Since the 1960s, the accumulation of the RES data has prompted the development of automated methods to extract the englacial layers. In this study, we propose a new automated method that combines peak detection methods, namely the CWT-based peak detection or the Automatic Phase Picker (APP), with a Hough Transform (HT) to trace boundaries of englacial layers. For CWT-based peak detection, we test it using two different wavelets. The proposed method is tested with twelve MCoRDS radio echograms, which are acquired south of the Northern Greenland Eemian (NEEM) ice drilling site, where the folding of ice layers was observed. The method is evaluated in comparison to the isochrones that were extracted in an independent study. In comparison, the proposed new automated method can restore more than 70% of the englacial layers. This new automated layer-tracing method is publicly available on github.


Introduction
The radio echo sounding (RES) technique has been used in glaciology since the 1960s to map the bed topography and layering of the ice sheets in Greenland and Antarctica [1].RES operates at low frequencies that are usually tens of megahertz, to achieve penetration depths down to thousands of meters into the ice bodies.Through this technique, pulsed waves are transmitted into the ice body and reflections from interfaces in the ice body are detected, which mark changes in the dielectric properties between adjacent layers.Analysis of these internal reflections, which are thought to be isochronous except for the deepest layers, has shed lights on many aspects of ice dynamics, such as past accumulation rates [2,3], ice flow changes [4,5], bed variations and its relationship with surface conditions [6,7].
Over the last five decades, radar sounders have been collecting large volumes of data, which have revealed internal ice structures of the Greenland Ice Sheet (GrIS) and Antarctica Ice Sheet (AIS).RES data have been widely applied for mapping ice thickness, bed topography and the internal structure of these large ice bodies [8][9][10][11].The Operation IceBridge (OIB) mission is the largest airborne survey of the Earth's polar ice ever flown.It is partially carried out by NASA to fill the temporal gap (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) between the termination of the Ice, Cloud and Land Elevation Satellite (ICESat) mission, and the launch of ICESat-2.The MCoRDS developed by the Center for Remote Sensing of Ice Sheets (CreSIS) is one of the five radar instruments [12].Compared to the snow radar and the Ku-band radar, the MCoRDS works in the frequency range from 150 MHz to 195 MHz, leading to a penetration depth of more than 3 km in the GrIS.
RES data are very useful in studies of bedrock elevations, in which only the bottom layers imaged by RES are used.Apart from this, a large number of englacial layers, which are thought to be isochrones, are also observed in RES data.However, analysis of the englacial layers by using a large volume of RES data is limited by the low efficiency associated with manual tracing.The englacial layers usually exhibit less electromagnetic contrast than surface and bed layers; therefore, tracing englacial layers can be a laborious and time-consuming manual task.For instance, tracing 20 layers in a 20,000 km RES survey requires a skilled operator 10 years to complete [13].The increasing volume of RES data has motivated the development of automated or semi-automated methods to trace the internal reflections.Most of the layer-tracing methods involve manual selection of seed points and are based on the estimation of slope angles, which are used for propagating the seed points [14].The propagation is called peak-following.Sime et al. proposed a method for deriving layer slope using a binary thresholding system but it struggles with noisy and discontinuous data [13].Paton proposed two automated methods for tracing the internal reflections, one of which produced an initial estimation of layer positions and the other, known as the SNAKE, traced a layer from the initial estimation [14], but this method produced many false positive results for low-amplitude reflectors [15].MacGregor et al. presented two phase-based methods to infer the layer slope and to estimate the ice ages of traced englacial layers by using ice drilling core data [16]; however, this method required complex radar data.In the most commonly applied peak-following method, one critical problem is that the layer-tracing is liable to steer into neighboring layers in regions where the visible stratigraphy deteriorates [14].Instead of peak-following, Onana et al. and Holschuh et al. proposed to use the Radon Transform (RT) to trace layers with an assumption that a layer segment over a short distance is a straight line.Their methods were tested with Ku-band radar and MCoRDS data respectively [15,17].Another tracing scheme was tested with three frames of MCoRDS data near NEEM ice core where the folding of ice was observed [18].This tracing scheme combined the RT and Fahnestock et al.'s method, which was based on cross-correlation and peak-following [19].The layers traced by using this tracing scheme were discontinuous because (1) there was no way to discriminate small peaks that might be caused by instrumental or environment noises, and (2) the radar signal attenuation caused inaccurate estimation of slope angle.Besides, correlation calculation of too many peak-pairs limited the efficiency.
Challenges in the task of layer-tracing are to separate closely spaced layers, to bridge gaps such as disconnected layers due to no radar returns and to tackle poor echograms quality in some regions [14].In this study, we propose a new automated method for layer-tracing, which combines the CWT-based or the APP-based peak detection and the Hough Transform (HT) to improve the reliability and efficiency of the layer-tracing task.The proposed method avoids the manual selection of seed points by evaluating the peak prominence with wavelet coefficients or short-term to long-term ratios.The following layer-tracing is arranged in a hierarchical manner, which is consistent with the process of human visualization.We test the proposed method with radio echograms from the pre-IceBridge and the OIB missions.The test data and its geographical coverage are introduced in Section 2. A detailed description of the proposed method for tracing the englacial layers is presented in Section 3. Following in Section 4, the results are presented and compared to isochrones extracted in an independent study by MacGregor et al. [16].Finally, a summary of the benefits and disadvantages of the proposed method is discussed and concluded in Sections 5 and 6.The open source codes for the proposed method can be found at https://github.com/xiongsiting/layer-tracing.

NEEM Ice Drilling Core Site
The study site is located in northwestern Greenland and is south of the NEEM site.This region is an area of great interest within the OIB mission as it is located on an ice divide and is less than 100 km to the englacial tunnel revealed within the Petermann glaciers [20].In addition, disrupted layers are observed in the RES data in this area [21,22].Folding of the internal layers close to the bed was observed in the ice core and RES data [22].To test the proposed method for tracing englacial layers, we investigate 12 transects of MCoRDS data, the flight tracks of which are shown in Figure 1.The geographic information of these transects is listed in Table 1.

NEEM Ice Drilling Core Site
The study site is located in northwestern Greenland and is south of the NEEM site.This region is an area of great interest within the OIB mission as it is located on an ice divide and is less than 100 km to the englacial tunnel revealed within the Petermann glaciers [20].In addition, disrupted layers are observed in the RES data in this area [21,22].Folding of the internal layers close to the bed was observed in the ice core and RES data [22].To test the proposed method for tracing englacial layers, we investigate 12 transects of MCoRDS data, the flight tracks of which are shown in Figure 1.The geographic information of these transects is listed in Table 1.

MCoRDS L1B Data
The MCoRDS is a nadir-looking, five-channel, monostatic radar system [12], which has the capability to penetrate up to several kilometers into the polar ice sheets.Before 2010, the radar system operated at a central frequency of 150 MHz, and a bandwidth of 17 MHz with a peak transmit power of 200 W. From 2010 onwards, the central frequency was changed to 195 MHz with a bandwidth of 30 MHz and with increasing peak transmit power up to 1200 W [12,25].
The MCoRDS data are in the form of radio echograms, in which the x-axis is geographical locations (latitude and longitude) and the y-axis denotes the two-way travel time or depth (assuming a single dielectric constant of 3.15 for pure ice).Each column of one radio echogram is a two-way

MCoRDS L1B Data
The MCoRDS is a nadir-looking, five-channel, monostatic radar system [12], which has the capability to penetrate up to several kilometers into the polar ice sheets.Before 2010, the radar system operated at a central frequency of 150 MHz, and a bandwidth of 17 MHz with a peak transmit power of 200 W. From 2010 onwards, the central frequency was changed to 195 MHz with a bandwidth of 30 MHz and with increasing peak transmit power up to 1200 W [12,25].
The MCoRDS data are in the form of radio echograms, in which the x-axis is geographical locations (latitude and longitude) and the y-axis denotes the two-way travel time or depth (assuming a single dielectric constant of 3.15 for pure ice).Each column of one radio echogram is a two-way travel time signal recording a sequence of internal reflections that are characterized as the amplitude peaks.The data are released in two versions [26,27], the first of which has a temporal coverage from 16 October 2009 to 6 November 2012, and the second from 12 October 2012 to the present.Data acquired on each day are divided into segments.A segment is a contiguous dataset in which the radar settings do not change and each segment is broken into frames, each of which is approximately 50 km long.Meta information including the elevation of the platform, the latitude and longitude of the flight track, surface and bed layers are also released along with each frame of the radio echogram.The 12 transects of radio echograms were acquired for the period of 2011-2014, detailed information of which is listed in Table 1.

Ancillary Datasets
In 2015, MacGregor et al. completed a comprehensive study of tracing englacial layers of the entire Greenland by using the Pre-IceBridge and the early IceBridge mission data from 1993 to 2013.In addition, they related the traced englacial layers to ice age estimated from ice drilling core data.In their work, the layer-tracing task was realized by applying their phase-based methods to the complex radar data.For data acquired prior to 2006, they used the Automatic Radio Echo Sounding Processing (ARESP) algorithm proposed by Sime et al. in 2011 or manual tracing.The seed points were manually selected and then they are used to propagate layer-tracing.By propagating the englacial layers observed in ice cores, they traced the isochrones to the areas in the GrIS where there were MCoRDS data.These isochrones are publicly available in the National Snow and Ice Data Center (NSIDC) database as a product called Radiostratigraphy and Age Structure of the Greenland Ice Sheet (RRRAG) [28].
The phase-based methods were used in their study to calculate the layer-slope.The following layer-tracing was basically a peak following method.The phase-based methods require complex radar data in which the phase information is preserved; however, the MCoRDS data are released as amplitude in the NSIDC database.In this study, these isochrones are used for validating englacial layers traced by our proposed method.

MCoRDS Radio Echograms of GrIS and an Automated Tracing Strategy
Radio echograms acquired by the MCoRDS show deeper englacial layers down to a depth of more than 2 km, which is different from radio echograms acquired by higher frequency radars such as the Ku-band radar mapping fern and snow [17].The entire Greenland ice column consists of ice from the Holocene, glacial ice and Eemian interglacial ice [22].An example of an echogram is shown in Figure 2. Continuous ice layers are well preserved in the upper ∼2.2 km of the ice sheet, which corresponds to ice from the Holocene and includes the last glacial.Few continuous layers and the disrupted layers imply ice folding below ∼2.2 km, which is related to the Eemian interglacial ice [22].The radar signal below the bed interface is below the noise level; therefore, it is considered invalid for indicating any interfaces or features.The drop of the radar amplitude at the bottom of the radio echogram may result from the instrumental issues as this is not observed in data from the pre-IceBridge mission.In the pre-processing, we remove the bottom part (~300 rows) of the radio echograms, which contains a sharp drop in the radar amplitude.According to previous studies, tracing englacial layers from radio echograms can be generally separated into three steps, which are: (1) selecting reliable seed points; (2) calculating layer slope angles and (3) propagating layers from seed points with the calculated slope angles.Prominent englacial layers manifest in radio echograms with high amplitude contrast to surrounding areas.Thus, Peak detection can be used to extract the points that belong to layers.However, the radar return strength between the two dashed lines in Figure 2c diminishes with depth due to beam divergence and energy losses in the medium [19].In this case, high amplitudes do not always guarantee real peaks since noise is mingled with informative signals and shows up as high amplitude peaks as well [19].Conversely, low amplitude peaks can still be real, such as peaks located between the lower dashed lines and the bed layers in Figure 2c.The CWT-based method can address this problem without de-trending the signal prior to peak detection, which may lead to inconsistent peak detection [29].
In this study, we propose a new method of using the CWT-based peak detection [29], or APPbased peak detection [30], and the HT [31,32] to trace the englacial layers from the MCoRDS radio echograms.A processing flowchart of this new automated method is shown in Figure 3. First, an image with amplitude peaks is produced by using either a CWT-based or APP-based peak detection.In Section 3.2, the production of this peak image by these two algorithms is described in detail.In the second step, seed points are selected by applying a threshold, which is calculated by using a lognormal fitting to the generated image.Then, the Hough Transform (HT) is applied to the peak image to calculate the slope angles during the tracing process, which starts from the selected seed points.The tracing procedure is described in Section 3.4.Finally, a post-processing process is applied to the traced englacial layers to connect the disconnected layer segments and to geocode the traced englacial layers.The post-processing is described in Section 3.5.According to previous studies, tracing englacial layers from radio echograms can be generally separated into three steps, which are: (1) selecting reliable seed points; (2) calculating layer slope angles and (3) propagating layers from seed points with the calculated slope angles.Prominent englacial layers manifest in radio echograms with high amplitude contrast to surrounding areas.Thus, Peak detection can be used to extract the points that belong to layers.However, the radar return strength between the two dashed lines in Figure 2c diminishes with depth due to beam divergence and energy losses in the medium [19].In this case, high amplitudes do not always guarantee real peaks since noise is mingled with informative signals and shows up as high amplitude peaks as well [19].Conversely, low amplitude peaks can still be real, such as peaks located between the lower dashed lines and the bed layers in Figure 2c.The CWT-based method can address this problem without de-trending the signal prior to peak detection, which may lead to inconsistent peak detection [29].
In this study, we propose a new method of using the CWT-based peak detection [29], or APP-based peak detection [30], and the HT [31,32] to trace the englacial layers from the MCoRDS radio echograms.A processing flowchart of this new automated method is shown in Figure 3. First, an image with amplitude peaks is produced by using either a CWT-based or APP-based peak detection.In Section 3.2, the production of this peak image by these two algorithms is described in detail.In the second step, seed points are selected by applying a threshold, which is calculated by using a lognormal fitting to the generated image.Then, the Hough Transform (HT) is applied to the peak image to calculate the slope angles during the tracing process, which starts from the selected seed points.The tracing procedure is described in Section 3.4.Finally, a post-processing process is applied to the traced englacial layers to connect the disconnected layer segments and to geocode the traced englacial layers.The post-processing is described in Section 3.5.

CWT-Based Peak Detection
As mentioned above, profiles from radio echograms can be regarded as two-way travel time signals recording the internal reflections.The reflections occur at the boundaries between two layers of different dielectric properties.These boundaries are amplitude peaks on the radar signal.The goal of tracing layers on the radio echogram is to find the peak signals and to connect them to continuous layers.In this paper, we demonstrate a CWT-based peak detection to extract the peaks in every radar profile while preserving the overall layering feature in the radio echogram.
The CWT-based peak detection has been devised to identify peaks with different scales and amplitudes [29,33].It can preserve the location and frequency at the same time.Therefore, when analyzing the profile from radio echograms, it can produce a series of wavelet coefficients, which represent the correlation between the radar signals and the selected filter.The two-way travel time radar signals, ( ), can be convolved with a wavelet, ( ).The derived wavelet coefficients can be expressed as Equation (2), where a ∈ ℝ * is the scale of the selected wavelet and b ∈ ℝ is the translational value.
Applying the CWT to a radar signal with varying scales produces a scalogram, which shows the wavelet coefficients over time (or depth) and scales.Since no prior knowledge about the location of peaks is known, the CWT wavelet moves along the radar signal, step by step, which means the translational value, , is set to be 1.The two factors affecting the coefficient are the type of wavelet and the scale.
Two types of the wavelet, namely Morlet (MORL) and Mexican hat (MEXH) are analyzed with an example profile (1500th column) from radio echograms of 20110329_02_020.The derived

CWT-Based Peak Detection
As mentioned above, profiles from radio echograms can be regarded as two-way travel time signals recording the internal reflections.The reflections occur at the boundaries between two layers of different dielectric properties.These boundaries are amplitude peaks on the radar signal.The goal of tracing layers on the radio echogram is to find the peak signals and to connect them to continuous layers.In this paper, we demonstrate a CWT-based peak detection to extract the peaks in every radar profile while preserving the overall layering feature in the radio echogram.
The CWT-based peak detection has been devised to identify peaks with different scales and amplitudes [29,33].It can preserve the location and frequency at the same time.Therefore, when analyzing the profile from radio echograms, it can produce a series of wavelet coefficients, which represent the correlation between the radar signals and the selected filter.The two-way travel time radar signals, x(t), can be convolved with a wavelet, ψ(t).The derived wavelet coefficients can be expressed as Equation (2), where a ∈ R + * is the scale of the selected wavelet and b ∈ R is the translational value.
Applying the CWT to a radar signal with varying scales produces a scalogram, which shows the wavelet coefficients over time (or depth) and scales.Since no prior knowledge about the location of peaks is known, the CWT wavelet moves along the radar signal, step by step, which means the translational value, b, is set to be 1.The two factors affecting the coefficient are the type of wavelet and the scale.
Two types of the wavelet, namely Morlet (MORL) and Mexican hat (MEXH) are analyzed with an example profile (1500th column) from radio echograms of 20110329_02_020.The derived scalograms are shown in Figure 4a,b.Each row in the scalogram indicates the wavelet response to the radar amplitude variations at one scale.By screening through each row in the scalogram, the maximum responses to the wavelet can be found.The red line in Figure 4d shows radar returns below the bed reflection.Generally, there are no radar reflections below this area, so we use the maximum response from this part as a threshold, below which the local maxima can be discarded.The elimination is carried out on every scale and only the peaks above the background noise level are preserved.What we are interested in are the positions that possess high wavelet coefficients at every scale, thus, wavelet coefficients from all scales are summed at each position to produce a Coefficient Summation (CS), which is shown in Figure 4c,d.Both types of wavelet detect the real peaks of the radar signal, while the MEXH wavelet is more consistent in detecting the peaks at all scales compared with the MORL wavelet.Moreover, MEXH depresses the spikes that surround the real peaks, resulting in a higher CS of peaks.Thus, the peaks are more prominent when the MEXH wavelet is used.Starting from around a scale of 15, the MEXH wavelet creates more side spikes around the peaks between a depth of 1.0-1.5 km.The MORL wavelet is unable to suppress the noise effects at the small scale.If noise elimination is carried out at each scale, then the smallest scale need to be chosen cautiously, otherwise, the real peaks can be masked by a large number of false peaks.Consequently, we recommend using the MEXH wavelet in the peak detection.
To choose a proper set of scales for the CWT, the CS values on the peaks and the number of these peaks being detected in the scalogram are plotted in Figure 5a-c.The three plots present the CS values from column Nos.500, 1500 and 2500 of the radio echogram, 20110329_02_020, respectively.A high CS confirms with that either prominent peak signals or long-lasting mediate signals, such as the signals on column no.2500 shown in Figure 5c.The spikes are evaluated as low CS. Figure 5d-f shows the relationship between scales and the total number of peaks along the radar signal in Figure 5a-c, which are detected at each scale.These plots show that the total number of the detected peaks drops enormously around a scale of 3 and almost all peaks are detected within a scale of 15, which implies that most of the layering information is contained within scales from 3 to 15.Therefore, in the following study, we use a scale set of {3,4,5, . . .15} for detecting the peaks and calculating their CS values.

APP-Based Peak Detection
An alternative peak detection method is based on the Automatic Phase Picker (APP), which is commonly used in the seismological literature.Allen et al. proposed the use of the short-term average/long-term average (STA/LTA) method for peak detection [34][35][36][37].In this study, this peak detection method is improved by using the kurtosis, which is a higher order statistics for Gaussian distributions.The input of the APP-based peak detection is also the radar signal.First, two windows of different intervals are centered on all of the peaks, which are local maxima, of the radar signal.The ratio of the short-term average to the long-term average of the radar amplitude is calculated.The kurtosis of the radar signals within the long-term window is calculated by using Equation (4), where is the long-term window length and ̅ is the average value of the radar signal of the long-term window.
The STA/LTA ratio and the kurtosis are not only calculated with a single pair of short-term and long-term.Instead, we apply a set of short-term and long-term pairs to derive the STA/LTA ratios and the corresponding kurtosis.In addition, different weights are assigned to different STA/LTA ratios and the kurtosis calculated from different pairs of short-term and long-term lengths.The weights, which are between 0 and 1, are set to be decrease with increasing short-term and long-term lengths.The summation of the weights is 1.This weighting scheme allows the selection of narrow and high peaks.Finally, the summation of STA/LTA ratios is compared to the summation of kurtosis.The peaks with an overall STA/LTA ratio larger than 1/3 of its overall kurtosis are preserved and the others are discarded.The APP-based peak detection returns the summation of STA/LTA ratios as an output.Input parameters for CWT-based and APP-based peak detection are listed in Table 2.

APP-Based Peak Detection
An alternative peak detection method is based on the Automatic Phase Picker (APP), which is commonly used in the seismological literature.Allen et al. proposed the use of the short-term average/long-term average (STA/LTA) method for peak detection [34][35][36][37].In this study, this peak detection method is improved by using the kurtosis, which is a higher order statistics for Gaussian distributions.The input of the APP-based peak detection is also the radar signal.First, two windows of different intervals are centered on all of the peaks, which are local maxima, of the radar signal.The ratio of the short-term average to the long-term average of the radar amplitude is calculated.The kurtosis of the radar signals within the long-term window is calculated by using Equation (4), where n is the long-term window length and x is the average value of the radar signal of the long-term window.
The STA/LTA ratio and the kurtosis are not only calculated with a single pair of short-term and long-term.Instead, we apply a set of short-term and long-term pairs to derive the STA/LTA ratios and the corresponding kurtosis.In addition, different weights are assigned to different STA/LTA ratios and the kurtosis calculated from different pairs of short-term and long-term lengths.The weights, which are between 0 and 1, are set to be decrease with increasing short-term and long-term lengths.The summation of the weights is 1.This weighting scheme allows the selection of narrow and high peaks.Finally, the summation of STA/LTA ratios is compared to the summation of kurtosis.The peaks with an overall STA/LTA ratio larger than 1/3 of its overall kurtosis are preserved and the others are discarded.The APP-based peak detection returns the summation of STA/LTA ratios as an output.Input parameters for CWT-based and APP-based peak detection are listed in Table 2.

Selection of Seed Points
Applying the CWT to all profiles in a radio echogram results in a CS image that only shows the peaks with positive values of CS while values of other pixels are set to be zero.Points of higher CS value have a higher probability to be located on the englacial layers.
In manual layer-tracing, tracing englacial layers starts from points that possess a relatively large contrast of radar echo amplitude with the surrounding pixels.Besides, it is natural to trace the prominent layers before the faint ones.Comparatively, the automated layer-tracing should start from the points that obviously belong to layers.The CS image is a good evaluator for choosing seed points.Tracing all layers from a radio echogram should be ordered by the CS value of the selected seed points.
Figure 6 shows the Cumulative Distribution Functions (CDF) of the CS image produced by using the MCoRDS data (product IDs: 20110329_02_019 and 20110329_02_020).The CS values fit well with the lognormal distribution, which means that many points have low CS value.We only choose the points of CS larger than the expectation, E[X] of the lognormal distribution as seed points and arrange them with descending CS value, S n (n = 1, 2, 3, . . ., N).

CWT-Based Peak Detection APP-Based Peak Detection MEXH or MORL wavelet Long term windows Wavelet scales
Short term windows weights Pixels below bed to choose noisy signals

Selection of Seed Points
Applying the CWT to all profiles in a radio echogram results in a CS image that only shows the peaks with positive values of CS while values of other pixels are set to be zero.Points of higher CS value have a higher probability to be located on the englacial layers.
In manual layer-tracing, tracing englacial layers starts from points that possess a relatively large contrast of radar echo amplitude with the surrounding pixels.Besides, it is natural to trace the prominent layers before the faint ones.Comparatively, the automated layer-tracing should start from the points that obviously belong to layers.The CS image is a good evaluator for choosing seed points.Tracing all layers from a radio echogram should be ordered by the CS value of the selected seed points.
Figure 6 shows the Cumulative Distribution Functions (CDF) of the CS image produced by using the MCoRDS data (product IDs: 20110329_02_019 and 20110329_02_020).The CS values fit well with the lognormal distribution, which means that many points have low CS value.We only choose the points of CS larger than the expectation, of the lognormal distribution as seed points and arrange them with descending CS value, ( = 1, 2, 3, … , ).

Layer-Tracing Procedure
Taking advantage of the CS image and the selected seed points, tracing layers can be carried out in the radio echogram.For each seed point, one layer-tracing procedure is accomplished based on three assumptions, which are (1) the englacial layers are quasi-parallel; (2) there are no crossing layers; and (3) within a short distance, the englacial layers can be approximated by a straight line [15,17,18].
Figure 7 illustrates the layer-tracing procedure.Propagation from the seed point is separated into two directions, which are parallel or anti-parallel to the flight direction.Starting from the seed point and the CS image, a block is centered on the seed point and the layer follows straight lines

Layer-Tracing Procedure
Taking advantage of the CS image and the selected seed points, tracing layers can be carried out in the radio echogram.For each seed point, one layer-tracing procedure is accomplished based on three assumptions, which are (1) the englacial layers are quasi-parallel; (2) there are no crossing layers; and (3) within a short distance, the englacial layers can be approximated by a straight line [15,17,18].
Figure 7 illustrates the layer-tracing procedure.Propagation from the seed point is separated into two directions, which are parallel or anti-parallel to the flight direction.Starting from the seed point and the CS image, a block is centered on the seed point and the layer follows straight lines within the block as mentioned in assumption (3).For propagating the layer, the slope angle of the layer on the seed point needs to be known.The angle is detected by converting the block CS image into the Hough Transform domain.The equation of the HT transform is shown in Equation ( 5) [32,33].
where (i, j) corresponds to row and column in the block, ρ ∈ (− 2 w b ) is the distance from the image centre to the line feature, θ ∈ [−π/2, π/2) is the angle between the line feature and the horizontal axis.The w b is the block size which is listed in Table 3.The global maximum point in the HT domain corresponds to a dominant line feature in the block, and its slope angle is the corresponding angle, θ.When the straight line is successfully traced in the block, the location of the seed point is changed to the end of this line.When the dominant angle is calculated due to lines close to the centre but no layers are related to the seed point (Figure 7d), or there is no dominant angle (Figure 7e), the layer-tracing is stopped.For this purpose, a line is constructed at the center with the calculated slope angle and a buffer is then created around the line.The buffer distance to the central line is d t .The points outside the buffer would be eliminated to generate a new block image, which is input into the HT slope estimation again.If the CS points on the line are less than a percentage, p ht , to the block size, no slope angle is returned for further tracing.
Remote Sens. 2018, 10, 43 10 of 22 within the block as mentioned in assumption (3).For propagating the layer, the slope angle of the layer on the seed point needs to be known.The angle is detected by converting the block CS image into the Hough Transform domain.The equation of the HT transform is shown in Equation ( 5) [32,33].
where ( , ) corresponds to row and column in the block, ∈ (− √ , √ ) is the distance from the image centre to the line feature, ∈ − /2, /2) is the angle between the line feature and the horizontal axis.The is the block size which is listed in Table 3.The global maximum point in the HT domain corresponds to a dominant line feature in the block, and its slope angle is the corresponding angle, .When the straight line is successfully traced in the block, the location of the seed point is changed to the end of this line.When the dominant angle is calculated due to lines close to the centre but no layers are related to the seed point (Figure 7d), or there is no dominant angle (Figure 7e), the layer-tracing is stopped.For this purpose, a line is constructed at the center with the calculated slope angle and a buffer is then created around the line.The buffer distance to the central line is .The points outside the buffer would be eliminated to generate a new block image, which is input into the HT slope estimation again.If the CS points on the line are less than a percentage, , to the block size, no slope angle is returned for further tracing.In a summary, layer-tracing stops when it comes to anyone of the following conditions: (1) no slope angle is returned; (2) crossing a labelled (or successfully traced) line; (3) approaching a labelled line with less than a distance, and (4) exceeding a slope angle difference, ∆ .Table 3 lists all the parameters used in the layer tracing procedure.It is noteworthy that the numbers of labelled layers can be much less than the number of seed points because after one layer-tracing procedure, seed points which are close to a labelled layer less than a distance of , are removed from the set of seed points, .
There is no manual intervention in the layer-tracing procedure.However, several input parameters are critical to the derived results, which are , and ∆ (see Table 3).The block size, , controls the continuity of the layers since the larger the block, the more the surrounding layers are taken into account for estimating the slope angle, which would lead to a successful connection of discontinuous points appearing to be on one layer.A large block size leads to missing of the detailed curvature of layers due to the assumption of straight lines within the block.For example, if the block In a summary, layer-tracing stops when it comes to anyone of the following conditions: (1) no slope angle is returned; (2) crossing a labelled (or successfully traced) line; (3) approaching a labelled line with less than a distance, d t and (4) exceeding a slope angle difference, ∆θ s .Table 3 lists all the parameters used in the layer tracing procedure.It is noteworthy that the numbers of labelled layers can be much less than the number of seed points because after one layer-tracing procedure, seed points which are close to a labelled layer less than a distance of d t , are removed from the set of seed points, S n .
There is no manual intervention in the layer-tracing procedure.However, several input parameters are critical to the derived results, which are w b , d t and ∆θ s (see Table 3).The block size, w b , controls the continuity of the layers since the larger the block, the more the surrounding layers are taken into account for estimating the slope angle, which would lead to a successful connection of discontinuous points appearing to be on one layer.A large block size leads to missing of the detailed curvature of layers due to the assumption of straight lines within the block.For example, if the block size is set to be larger than the width of a fold, the curve of the fold may not be fully reconstructed.On the other hand, when the radio echogram suffers a loss of signal, the block can be adjusted to a larger size to improve layer continuity.The minimum distance, d t , controls the density of the traced layers.It is within this vertical distance that only one prominent layer can be traced, while any subsequent tracing is forbidden within this range.Though faint layers can be traced with decreasing values of d t , they may be short layers because the tracing may be stopped when they fade into the background.The slope angle difference, ∆θ s , is the angle difference between the current estimated slope angle with the last one from the previous block.This parameter controls how small the curves can be restored in the layer-tracing.The larger the ∆θ s is set, the smaller the folds that can be followed during tracing, otherwise, only smooth large folds can be restored.In some applications, such as tracing layers from ice drilling cores, prior knowledge about the location of layers may be known.In this case, manually selected seed points are easily imported into the layer-tracing module.

Post-Processing of Traced Englacial Layers
Since any complex radar sounding system can result in noise during data acquisition, a continuous ice sheet layer may be imaged as discontinuous layer segments in the radio echogram.When inspected manually, it is sometimes possible to connect the layer segments into the same layer.However, since the end condition has been set, these layer segments are usually coded with different labels during the tracing procedure.Therefore, it is necessary to post-process the traced englacial layers to connect the discontinuous layer segments that belong to the same layer.
Firstly, we go through all the traced layers and record their start points, end points and their initial labels.For one target layer segment, as illustrated by L t in Figure 8, all layer segments having their start point right to the end of L t are initially selected (purple lines in Figure 8) as the candidates for connection to the target layer segment L c .Then, for each of these candidates, a reference layer that is continuous between the end of L t and the start of the candidate layers is picked out (red continuous line in Figure 8).If the target layer segment and the candidate are located on different sides of the reference layer, no connection occurs, such as Figure 8a.Otherwise, if the difference of their vertical distance to the reference layer is below a threshold ∆d in Table 3, which here is |d 1 − d 2 | < ∆d, the candidate layer segment is connected to the target one (Figure 8b).The threshold for horizontal distance is not considered here because it is possible that line segments are disrupted for a long distance while they still may belong to the same layer.If two layer segments are close to a reference common layer, it is reasonable to think they are the same layer at least within one frame.This checking is repeated for the start point of the target layer segment for the left connection.
After the connection of the traced layer segments, a threshold for the minimum horizontal distance can be set up to eliminate any remaining short layers.This threshold can be tuned for different applications.Finally, the remaining layers are geocoded with the latitude, longitude and elevation from the meta information, which is available as part of the MCoRDS data.

Producing the CS Image and APP Ratio Image
To test this automated layer tracing method shown schematically in Figure 3, we first experiment it with one frame of a radio echogram as an example.We choose the 20th frame of the 2nd segment from the MCoRDS data acquired on 29 March 2011 as it recorded two ice folding areas.The echogram has 1839 rows and 3748 columns.One pixel in the image corresponds to a horizontal distance of ~13.6 m and a depth of ~2.8 m.The radar amplitude is shown in Figure 9. Figure 10a,b shows the MEXH CS image and seed points filtered by lognormal fitting.Black lines shown in Figure 9 represent the surface and bed boundaries.The CS image is produced by using a CWT-based peak detection with a MEXH wavelet and the scales of {3, 4, 5, …, s}.There are 210,487 peaks in the MEXH CS image.By using the lognormal distribution fitting, points of CS value larger than 44.5 are selected as seed points, which results in 59,402 points.Producing the CS image takes only ~60 s.Previously, it was shown in Section 3.2.1, that if small scales are screened, almost all pixels are selected by using the MORL wavelet.For comparison, the same radio echogram is processed by using a MORL wavelet with scales of 10-15.Figure 11 shows the CS image produced by using the MORL wavelet.Figure 12 shows the APP ratio image.The sets of short-term and long-term windows and weights are adjusted heuristically as {3, 5, 5, 7, 7, 9}, {5, 7, 9, 11, 13, 13} and {0.3, 0.25, 0.2, 0.15, 0.05, 0.05}.Producing the ratio image takes ~750 s as higher order statistics result in time complexity of

Producing the CS Image and APP Ratio Image
To test this automated layer tracing method shown schematically in Figure 3, we first experiment it with one frame of a radio echogram as an example.We choose the 20th frame of the 2nd segment from the MCoRDS data acquired on 29 March 2011 as it recorded two ice folding areas.The echogram has 1839 rows and 3748 columns.One pixel in the image corresponds to a horizontal distance of ~13.6 m and a depth of ~2.8 m.The radar amplitude is shown in Figure 9. Figure 10a,b shows the MEXH CS image and seed points filtered by lognormal fitting.Black lines shown in Figure 9 represent the surface and bed boundaries.The CS image is produced by using a CWT-based peak detection with a MEXH wavelet and the scales of {3, 4, 5, . . ., s}.There are 210,487 peaks in the MEXH CS image.By using the lognormal distribution fitting, points of CS value larger than 44.5 are selected as seed points, which results in 59,402 points.Producing the CS image takes only ~60 s.

Producing the CS Image and APP Ratio Image
To test this automated layer tracing method shown schematically in Figure 3, we first experiment it with one frame of a radio echogram as an example.We choose the 20th frame of the 2nd segment from the MCoRDS data acquired on 29 March 2011 as it recorded two ice folding areas.The echogram has 1839 rows and 3748 columns.One pixel in the image corresponds to a horizontal distance of ~13.6 m and a depth of ~2.8 m.The radar amplitude is shown in Figure 9. Figure 10a,b shows the MEXH CS image and seed points filtered by lognormal fitting.Black lines shown in Figure 9 represent the surface and bed boundaries.The CS image is produced by using a CWT-based peak detection with a MEXH wavelet and the scales of {3, 4, 5, …, s}.There are 210,487 peaks in the MEXH CS image.By using the lognormal distribution fitting, points of CS value larger than 44.5 are selected as seed points, which results in 59,402 points.Producing the CS image takes only ~60 s.Previously, it was shown in Section 3.2.1, that if small scales are screened, almost all pixels are selected by using the MORL wavelet.For comparison, the same radio echogram is processed by using a MORL wavelet with scales of 10-15.Figure 11 shows the CS image produced by using the MORL wavelet.Figure 12 shows the APP ratio image.The sets of short-term and long-term windows and weights are adjusted heuristically as {3, 5, 5, 7, 7, 9}, {5, 7, 9, 11, 13, 13} and {0.3, 0.25, 0.2, 0.15, 0.05, 0.05}.Producing the ratio image takes ~750 s as higher order statistics result in time complexity of Previously, it was shown in Section 3.2.1, that if small scales are screened, almost all pixels are selected by using the MORL wavelet.For comparison, the same radio echogram is processed by using a MORL wavelet with scales of 10-15.Figure 11 shows the CS image produced by using the MORL wavelet.Figure 12 shows the APP ratio image.The sets of short-term and long-term windows and weights are adjusted heuristically as {3, 5, 5, 7, 7, 9}, {5, 7, 9, 11, 13, 13} and {0.3, 0.25, 0.2, 0.15, 0.05, 0.05}.Producing the ratio image takes ~750 s as higher order statistics result in time complexity of NlogN.The numbers of peaks in the CS images and ratio image and their seed points selected by using a lognormal distribution fitting are listed in Table 4.
Remote Sens. 2018, 10, 43 13 of 22 . The numbers of peaks in the CS images and ratio image and their seed points selected by using a lognormal distribution fitting are listed in Table 4.   .The numbers of peaks in the CS images and ratio image and their seed points selected by using a lognormal distribution fitting are listed in Table 4.

Layer-Tracing from Single Frame
To compare the MEXH CS image, MORL CS image and the APP ratio image, the seed points are selected from each of these images and then imported into the following layer-tracing.The numbers of layers traced by using these three images are listed in Table 4.The three sets of traced layers are shown in Figure 13.As we can observe, more layers following along the ice folds are restored by using the MEXH CS image.By using the MEXH CS image, more layers are detected especially in the folding areas of the radio echograms.After post-processing, the number of layers traced by using the MEXH CS image remains the highest among the three results.This may be due to more layers detected in the folding area or more layer segments than the other two results.The fully traced layers from the three methods are shown in the Supplementary Material.In this paper, several close-up look images are shown to evaluate the difference between them.
Table 5 lists time consumption and the number of traced layers when using the different sets of block size and the minimum allowed distance, from which we can see that increasing the block size does not necessarily increase the number of traced layers, while it slightly increases the numbers of connected layers during post-processing.Decreasing , largely increases the number of traced layers.

Layer-Tracing from Single Frame
To compare the MEXH CS image, MORL CS image and the APP ratio image, the seed points are selected from each of these images and then imported into the following layer-tracing.The numbers of layers traced by using these three images are listed in Table 4.The three sets of traced layers are shown in Figure 13.As we can observe, more layers following along the ice folds are restored by using the MEXH CS image.By using the MEXH CS image, more layers are detected especially in the folding areas of the radio echograms.After post-processing, the number of layers traced by using the MEXH CS image remains the highest among the three results.This may be due to more layers detected in the folding area or more layer segments than the other two results.The fully traced layers from the three methods are shown in the Supplementary Material.In this paper, several close-up look images are shown to evaluate the difference between them.
Table 5 lists time consumption and the number of traced layers when using the different sets of block size and the minimum allowed distance, from which we can see that increasing the block size does not necessarily increase the number of traced layers, while it slightly increases the numbers of connected layers during post-processing.Decreasing d t , largely increases the number of traced layers.

Combination of Multi-Frame and Validation
Although the smallest unit of the released MCoRDS data is one frame, it is necessary to combine several frames to form a result that reflects subsurface features along longer traverses than 50 km.Therefore, we demonstrate the combination of traced englacial layers from several continuous frames.Figure 14b shows an example that is produced from a combination of three frames (19)(20)(21) of the 2nd segment acquired on 29 March 2011.The layers are traced frame by frame and then are subsequently combined.It is possible that some layers are traced in one frame while the corresponding ones are not traced in the adjacent frame.

Combination of Multi-Frame and Validation
Although the smallest unit of the released MCoRDS data is one frame, it is necessary to combine several frames to form a result that reflects subsurface features along longer traverses than 50 km.Therefore, we demonstrate the combination of traced englacial layers from several continuous frames.Figure 14b shows an example that is produced from a combination of three frames (19)(20)(21) of the 2nd segment acquired on 29 March 2011.The layers are traced frame by frame and then are subsequently combined.It is possible that some layers are traced in one frame while the corresponding ones are not traced in the adjacent frame.
To assess the effectiveness of the proposed layer-tracing method, the corresponding result from MacGregor et al.'s radiostratigraphy product is shown in Figure 14a, in which the colors denote ice age.In Figure 14b, the majority of layers have correspondences in the published radiostratigraphy, the average of their elevation differences (absolute value) is within 40 m, which corresponds to a vertical distance of about 15 pixels.The RRRAG product has no follow-on data after 2013, thus, seven out of twelve lines in Figure 1, which were acquired in 2011, 2012 and 2013 have the corresponding results in the RRRAG product.All of them are compared to the results from this study, which are provided in the Supplementary Materials.Comparisons of the RRRAG product with the proposed new method demonstrate that the traced englacial layers are comparable to those traced by MacGregor et al. [16].The proposed method works well on restoring the orientations of the ice layers, even the faint layers related the Eemian interglacial ice are restored successfully, which may not even be delineated in the RRRAG product, such as the ones denoted as black lines in Figure 14.When overlaid on the radio echogram, the black lines are located on faint layers that consist of amplitude peaks.The RRRAG dataset was To assess the effectiveness of the proposed layer-tracing method, the corresponding result from MacGregor et al.'s radiostratigraphy product is shown in Figure 14a, in which the colors denote ice age.In Figure 14b, the majority of layers have correspondences in the published radiostratigraphy, the average of their elevation differences (absolute value) is within 40 m, which corresponds to a vertical distance of about 15 pixels.
The RRRAG product has no follow-on data after 2013, thus, seven out of twelve lines in Figure 1, which were acquired in 2011, 2012 and 2013 have the corresponding results in the RRRAG product.All of them are compared to the results from this study, which are provided in the Supplementary Materials.Comparisons of the RRRAG product with the proposed new method demonstrate that the traced englacial layers are comparable to those traced by MacGregor et al. [16].The proposed method works well on restoring the orientations of the ice layers, even the faint layers related the Eemian interglacial ice are restored successfully, which may not even be delineated in the RRRAG product, such as the ones denoted as black lines in Figure 14.When overlaid on the radio echogram, the black lines are located on faint layers that consist of amplitude peaks.The RRRAG dataset was accomplished by manual picking of seed points, whose number is limited.Therefore, during the manual operation, only the prominent layers are traced rather than all possible layers.Figure 16a shows the isochrones from the RRRAG product and Figure 16b shows traced layers that are associated with the ice ages.There are 106 layers that are traced by the proposed method, among which 57 layers are confirmed with 43 isochrones from the RRRAG product, in which there are 49 isochrones in total.From Figure 16b we can see that the lower layers are not restored as continuous as those in the RRRAG product.This may be improved by setting a smaller block size for smaller increments in tracing procedure.Besides, there are jumps near a depth of 1.5 km around a distance of 0-30 km, and unrealistic slopes around a distance of 107 km and a depth of 2.2 km.Setting a threshold for the slope angle difference should be able to avoid these unrealistic slopes, however, to restore the sharp changes of the lower part in the radio echogram, the threshold is not used in this case.The jumps and unrealistic slopes are also attributed to the low signal to noise ratio, which can be partly demonstrated by the same difficulty in tracing continuous layers seen in Figure 16a.
Remote Sens. 2018, 10, 43 18 of 22 minimum distance allowed between two close layers is set to be 3 pixels.The red line is the bed layer, which comes with the product.
Figure 16a shows the isochrones from the RRRAG product and Figure 16b shows traced layers that are associated with the ice ages.There are 106 layers that are traced by the proposed method, among which 57 layers are confirmed with 43 isochrones from the RRRAG product, in which there are 49 isochrones in total.From Figure 16b we can see that the lower layers are not restored as continuous as those in the RRRAG product.This may be improved by setting a smaller block size for smaller increments in tracing procedure.Besides, there are jumps near a depth of 1.5 km around a distance of 0-30 km, and unrealistic slopes around a distance of 107 km and a depth of 2.2 km.Setting a threshold for the slope angle difference should be able to avoid these unrealistic slopes, however, to restore the sharp changes of the lower part in the radio echogram, the threshold is not used in this case.The jumps and unrealistic slopes are also attributed to the low signal to noise ratio, which can be partly demonstrated by the same difficulty in tracing continuous layers seen in Figure 16a.

Discussion
This study demonstrates a processing method for tracing englacial layers from MCoRDS data in two steps, the first of which is to enhance the layering features while suppressing the noise.Three different methods, namely CWT-based peak detection with MEXH and MORL wavelet and APP-

Discussion
This study demonstrates a processing method for tracing englacial layers from MCoRDS data in two steps, the first of which is to enhance the layering features while suppressing the noise.Three different methods, namely CWT-based peak detection with MEXH and MORL wavelet and the method.Since the tracing procedure is carried out block by block, it results in short layer segments when the radio echogram is noisy.Therefore, in the future, more effective algorithms are needed to fill the gaps.

Conclusions
Tracing englacial layers is important for studying ice sheet dynamics yet difficult to achieve with high efficiency and accuracy at the same time.In this work, we propose an automated method that combines a peak detection method (CWT or APP), with an HT-based layer-tracing procedure to trace englacial layers from MCoRDS radio echograms in Greenland.Following the initially traced layers, we also propose a method for connecting discontinuous layer segments and for connecting layers traced from multiple frames of radio echograms.In general, it is difficult to quantitatively assess the effectiveness of the layer tracing algorithm.In our approach, we compare the results to isochrones of the published RRRAG product.This comparison shows that layers traced by the proposed method are comparable to those traced using the phase gradient and manual inspection.Moreover, the proposed method traces layers with much higher densities and restores the layers near the folding areas.In addition, it needs no direct manual assistance, other than the setting of several parameters, and its results show much more detail variation of layers when compared with the RRRAG product in areas of no severe signal loss in the radio echograms.Given a sufficiently large training set, it is likely that machine learning techniques may be applied in future to eliminate the parameter setting stage.The MATLAB package developed for this proposed method is available to the scientific community to try to stimulate further improvements in the future.
It is noteworthy that the traced layers are not necessarily directly related to the ice age; thus, in the future, methods of relating ice age to elevation (or depth) can be investigated to apply to these traced layers after matching up ice cores with their dielectric properties.In the future, we plan to adapt this method to trace layers observed in the Martian Polar Layered Deposits (PLDs) by radar sounders, such as the Mars Advanced Radar for Subsurface and Ionosphere Sounding (MARSIS) and SHAllow RADar (SHARAD).

Supplementary Materials:
The following are available online at www.mdpi.com/2072-4292/10/1/43/s1, Figure S1: The outcome layers from layer-tracing procedures with inputs from the MEXH CS image, the MORL CS image and APP ratio image independently.Parameters for layer-tracing are a block size of 51 pixels, a maximum distance allowance of 7 pixels and a slope angle difference of 90 • (The MCoRDS radio echogram is 20110329_02_020), Figure S2: Comparison between the traced layers by using the proposed method and those in the RRRAG product.

Figure 1 .
Figure 1.The study site and the flight tracks of MCoRDS data, which are selected for testing the layer tracing method.The background image is the bed topography [23,24].The yellow transects were acquired in 2011, red 2012, black 2013, white 2014.

Figure 1 .
Figure 1.The study site and the flight tracks of MCoRDS data, which are selected for testing the layer tracing method.The background image is the bed topography [23,24].The yellow transects were acquired in 2011, red 2012, black 2013, white 2014.

Figure 2 .
Figure 2. (a) An example of a MCoRDS radio echogram (subset from a whole frame radio echogram of product ID: 20110329_02_020) demonstrates the typical geometry of the englacial layers; (b) the close-up look image showing the folding of ice; (c) Profile 1 and Profile 2 show typical variations of the radar amplitude, which consists of a drop from the surface reflection value, a slower decreasing signal until several hundred meters above bed reflection, and a low reflection zone of several hundred meters above the bed.

Figure 2 .
Figure 2. (a) An example of a MCoRDS radio echogram (subset from a whole frame radio echogram of product ID: 20110329_02_020) demonstrates the typical geometry of the englacial layers; (b) the close-up look image showing the folding of ice; (c) Profile 1 and Profile 2 show typical variations of the radar amplitude, which consists of a drop from the surface reflection value, a slower decreasing signal until several hundred meters above bed reflection, and a low reflection zone of several hundred meters above the bed.

Figure 3 .
Figure 3. Processing flowchart for automatically tracing englacial layers.The dashed lines shown are the alternative way of using APP-based peak detection compared to the CWT-based peak detection.

Figure 3 .
Figure 3. Processing flowchart for automatically tracing englacial layers.The dashed lines shown are the alternative way of using APP-based peak detection compared to the CWT-based peak detection.

Figure 4 .
Figure 4. (a) the scalogram produced by the Mexican Hat wavelet transform; (b) the scalogram produced by the Morlet wavelet transform; (c) the CWT coefficient summation with MEXH and (d) the CWT coefficient summation with MORL.

Figure 4 .
Figure 4. (a) the scalogram produced by the Mexican Hat wavelet transform; (b) the scalogram produced by the Morlet wavelet transform; (c) the CWT coefficient summation with MEXH and (d) the CWT coefficient summation with MORL.

Figure 5 .
Figure 5. (a-c) Detection of peaks on three radar signals.(d-f) total number of peaks detected at each scale.

Figure 5 .
Figure 5. (a-c) Detection of peaks on three radar signals.(d-f) total number of peaks detected at each scale.

Figure 6 .
Figure 6.The lognormal distribution fitting of the wavelet coefficients in the CS image, which is generated by using the MORL wavelet with scales of {3, 4, 5,…15}.Examples are CS images of (a) 20110329_02_019 and (b) 20110329_02_020.The red vertical lines are the value of the expectation of lognormal distribution.

Figure 6 .
Figure 6.The lognormal distribution fitting of the wavelet coefficients in the CS image, which is generated by using the MORL wavelet with scales of {3, 4, 5, . . .15}.Examples are CS images of (a) 20110329_02_019 and (b) 20110329_02_020.The red vertical lines are the value of the expectation of lognormal distribution.

Figure 7 .
Figure 7. Illustration demonstrating the procedure of tracing englacial layers, (a) the layer-tracing is propagated in two directions; (b) right propagation and the moving of a seed point (from red to blue); (c) the geometry of the HT transform; (d,e) the cases in which layer-tracing stops.

Figure 7 .
Figure 7. Illustration demonstrating the procedure of tracing englacial layers, (a) the layer-tracing is propagated in two directions; (b) right propagation and the moving of a seed point (from red to blue); (c) the geometry of the HT transform; (d,e) the cases in which layer-tracing stops.

Figure 8 .
Figure 8.The illustrations of layer connection during post-processing: (a) no connection to the right of the target layer segment; (b) connection to right of the target layer segment (blue line).The red line is the reference layers and the purple lines are the candidate layer segments.

Figure 9 .
Figure 9. MCoRDS amplitude data that are used for testing the proposed method.Product ID is 20110329_02_020.The axes show row numbers and column numbers within the image.The white rectangles, A1, A2 and A3 denote areas will be shown in detail in the following text.

Figure 8 .
Figure 8.The illustrations of layer connection during post-processing: (a) no connection to the right of the target layer segment; (b) connection to right of the target layer segment (blue line).The red line is the reference layers and the purple lines are the candidate layer segments.

Figure 8 .
Figure 8.The illustrations of layer connection during post-processing: (a) no connection to the right of the target layer segment; (b) connection to right of the target layer segment (blue line).The red line is the reference layers and the purple lines are the candidate layer segments.

Figure 9 .
Figure 9. MCoRDS amplitude data that are used for testing the proposed method.Product ID is 20110329_02_020.The axes show row numbers and column numbers within the image.The white rectangles, A1, A2 and A3 denote areas will be shown in detail in the following text.

Figure 9 .
Figure 9. MCoRDS amplitude data that are used for testing the proposed method.Product ID is 20110329_02_020.The axes show row numbers and column numbers within the image.The white rectangles, A1, A2 and A3 denote areas will be shown in detail in the following text.

Figure 10 .Figure 11 .
Figure 10.(a) The CS image produced by using the CWT-based peak detection with MEXH wavelet and scales of {3, 4, 5, … , 15}, 210,487 points are picked out; (b) 59,402 seed points filtered out by using the E[X] = 44.5 from the lognormal distribution fitting.The color indicates the CS values.

Figure 10 .
Figure 10.(a) The CS image produced by using the CWT-based peak detection with MEXH wavelet and scales of {3, 4, 5, . . . ,15}, 210,487 points are picked out; (b) 59,402 seed points filtered out by using the E[X] = 44.5 from the lognormal distribution fitting.The color indicates the CS values.

Figure 10 .Figure 11 .
Figure 10.(a) The CS image produced by using the CWT-based peak detection with MEXH wavelet and scales of {3, 4, 5, … , 15}, 210,487 points are picked out; (b) 59,402 seed points filtered out by using the E[X] = 44.5 from the lognormal distribution fitting.The color indicates the CS values.

Figure 11 .
Figure 11.(a) The CS image produced by using the CWT-based peak detection with a MORL wavelet and the scales of {10, 11, 12, . . . ,15}. (b) Seed points filtered by using the expectation of lognormal distribution fitting.The number of peaks in the image and seed points are listed Table 4.The color indicates the CS values.

Figure 12 .
Figure 12.(a) The CS image produced by using the APP-based peak detection.(b) Seed points filtered by using the expectation of lognormal distribution fitting.The number of peaks in the image and seed points is listed Table 4.The color indicates the STA/LTA ratios.

Figure 12 .
Figure 12.(a) The CS image produced by using the APP-based peak detection.(b) Seed points filtered by using the expectation of lognormal distribution fitting.The number of peaks in the image and seed points is listed Table 4.The color indicates the STA/LTA ratios.

Figure 13 .
Figure 13.(a) Englacial layers traced by the three methods on the slope of the folds, indicated by A1 in Figure 9; (b,c) englacial layers traced near the folding cores, which are A2 and A3 in Figure 9.The white lines show layers traced by the MEXH CS image, the red ones the MORL CS image, and the black ones the APP ratio image.

Figure 13 .
Figure 13.(a) Englacial layers traced by the three methods on the slope of the folds, indicated by A1 in Figure 9; (b,c) englacial layers traced near the folding cores, which are A2 and A3 in Figure 9.The white lines show layers traced by the MEXH CS image, the red ones the MORL CS image, and the black ones the APP ratio image.

Figure 14 .
Figure 14.Comparison of the traced englacial layers from radio echogram from MCoRDS data 20110329_02_020, with those derived from RRRAG isochrones [13]: (a) is the RRRAG isochrones; (b) the englacial layers traced by the proposed method after removing layer segments extending a distance of shorter than 10 km.The color lines are matched layer segments to the RRRAG isochrones.The black lines are layer segments that have no matched isochrones.

Figure 14 .
Figure 14.Comparison of the traced englacial layers from radio echogram from MCoRDS data 20110329_02_020, with those derived from RRRAG isochrones [13]: (a) is the RRRAG isochrones; (b) the englacial layers traced by the proposed method after removing layer segments extending a distance of shorter than 10 km.The color lines are matched layer segments to the RRRAG isochrones.The black lines are layer segments that have no matched isochrones.

Figure 16 .
Figure 16.(a) the isochrones from the RRRAG product; (b) the traced layers after being associated with different ice ages.The layers between depths of 1.0-2.0km are better restored than those in the upper and lower parts.

Figure 16 .
Figure 16.(a) the isochrones from the RRRAG product; (b) the traced layers after being associated with different ice ages.The layers between depths of 1.0-2.0km are better restored than those in the upper and lower parts.

Table 1 .
MCoRDS data that are investigated in this study.

Table 1 .
MCoRDS data that are investigated in this study.

Table 2 .
Input parameters of CWT-based and APP-based peak detection.

Table 2 .
Input parameters of CWT-based and APP-based peak detection.

Table 3 .
Parameters used in producing the CS image and tracing layers.

Table 4 .
The number of peaks detected by using CWT-based and APP-based peak detection, and numbers of layers that are traced by using the CS images and ratio image.Parameters are set to the default values in Table3.

Table 4 .
The number of peaks detected by using CWT-based and APP-based peak detection, and numbers of layers that are traced by using the CS images and ratio image.Parameters are set to the default values in Table3.

Table 4 .
The number of peaks detected by using CWT-based and APP-based peak detection, and numbers of layers that are traced by using the CS images and ratio image.Parameters are set to the default values in Table3.

Table 5 .
The number of layers traced by using different block sizes, and different thresholds for the indicated vertical distances.

Table 5 .
The number of layers traced by using different block sizes, and different thresholds for the indicated vertical distances.