Recursive Enhancement of Weak Subsurface Boundaries and Its Application to SHARAD Data

Abstract

Sedimentary layers are composed of alternately deposited compositions in different periods, reflecting the geological evolution history of a planet. Orbital radar can detect sedimentary layers, but the radargram is contaminated by varying background noise levels. Traditional denoising methods, such as median filter, have difficulty dealing with such kinds of noise. We propose a recursive signal enhancement scheme to identify weak reflections from intense background noise. Numerical experiments with synthetic data and SHARAD radargrams illustrate that the proposed method can enhance the clarity of the radar echoes and reveal delicate sedimentary structures previously buried in the background noise. The denoising result presents better horizontal continuity and higher vertical resolution compared with those of the traditional methods.

1. Introduction

Mars polar caps comprise the majority of the planet’s known water ice reservoirs [1]. Their sedimentary layered structures, known as Martian polar layered deposits (PLDs), are believed to record the climate history of variant deposition rates and compositions [2]. During the exploration of the geology and stratigraphy of Martian PLDs, high-resolution imaging instruments, such as the MGS Mars Orbiter Camera (MOC) experiment [3] and the MRO High-Resolution Imaging Science Experiment (HiRISE) [4], returned abundant images of the polar regions of Mars. Subsequently, scientists have utilized these images at troughs and scarps of polar areas to analyze the morphology of the exposed sedimentary structures and inferred the internal structures (e.g., [5,6,7]).

Two radar sounder instruments, the Shallow Radar (SHARAD) [8] and the Mars Advanced Radar for Subsurface and Ionosphere Sounding (MARSIS) [9], are critical for detecting subsurface structures of Martian PLDs. Radar waves can propagate through the air and underground media. When encountering a specific change in dielectric permittivity, a portion of the energy is reflected back from the boundary of the layer [10]. After processing the backscattered power received by the antenna, we can delineate subsurface layers on a radargram—a two-dimensional image in which the vertical axis is the round-trip delay time and the horizontal axis is the along-track distance.

However, an unequivocal interpretation of the radargram is limited by the background noise, which unevenly increases the local mean grey level of a group of adjacent pixels. Studies attribute it to thermal and galactic background [8] and electromagnetic interference from other instruments or systems on the Mars Reconnaissance Orbiter [11]. The background noise of the SHARAD radargram is added to any noise that might be intrinsic to the instrument. For simplicity, the background noise is tackled as additive white Gaussian noise, which is stationary and uncorrelated among radargram cells.

In denoising synthetic aperture radar images, methods based on deep learning are prevalent [12,13,14]. However, these methods are applied mostly to the surface observations and rarely in subsurface sounding data because the noise-free SHARAD radargram is difficult to obtain. Other techniques involve using adaptive and non-adaptive filters on the signal processing; adaptive filters adapt their weightings to the noise level of the local area across the image, and non-adaptive filters utilize the same weightings uniformly across the whole image. Such filtering also eliminates subsurface reflections, especially for high-frequency components, and there are tradeoffs between the pertinency of filtering and the choosing of the filter type. Many forms of adaptive filtering have been studied extensively, such as the Lee filter [15,16], the Frost filter [17], and the Refined Gamma Maximum-A-Posteriori (RGMAP) filter [18]. Adaptive filtering is prominent in preserving edges and details. However, these all rely on specific mathematical assumptions [19] that can be harsh conditions for SHARAD radargram products and still produce visible artifacts, especially on sharp edges. Non-adaptive filtering is simple to implement and requires less capacity for calculation. There are two widespread applied forms of non-adaptive filtering: one is based on the moving average, and the other is based on the median (within a given window of pixels in the image). The former is better at enhancing continuity in a given direction, and the latter is better at preserving edges. Thus, we included the moving average as part of our method to improve lateral continuity.

Recently, the progressive image denoising (PID) method was proposed based on robust noise estimation [20]. It progressively reduces noise by deterministic annealing, which is heuristic and efficient at solving complex optimization problems in which non-single local extrema exist. We integrated the PID method and the moving-average method into our filtering scheme, the recursive subsurface boundary enhancement (RSE) method, aiming to filter out the background noise of the SHARAD radargram.

The rest of the paper is organized as follows. We first introduce the SHARAD data and the details of the RSE method in Section 2. In Section 3, we perform numerical experiments with both the synthetic model and the SHARAD radargram. Finally, in Section 4, we discuss the performance of the RSE method, and we conclude in Section 5.

2. Data and Method

The SHARAD onboard the Mars Reconnaissance Orbiter operates at a center frequency of 20 MHz, with a vertical resolution of ~15 m in free space and ~8 m in water ice (with a dielectric permittivity value of 3.15) [8]. The SHARAD team produces two sets of standard data products: Experiment Data Records (EDRs) and Reduced Data Records (RDRs). The former includes unprocessed echoes, and the latter consists of radargrams that are processed with the range compression and synthetic aperture technique [8]. We used the radargrams of RDRs for processing in our study. Although SHARAD observations for closely repeating tracks across the same area could produce radargrams with different noise levels because of the SA-HGA positions and MRO roll [11], we do not necessarily use the best quality radargram among the closely repeated tracks for the RSE method. The RSE method is expected to enhance the quality and readability of the whole radargrams with respect to RDRs.

Before introducing the framework of the RSE method, we concisely review the PID method. For a 2D SHARAD radargram d, which is composed of a signal s and additive white Gaussian noise n with variance σ², the PID method iteratively removes the estimated noise instance using bilateral filtering in the time–space domain and frequency–wavenumber domain [20,21]; thus, it can adapt to various noise situations by maintaining the edge information and tiny structures. The mathematical expression of the radargram is:

$$d=s+n,$$

(1)

where $d$ is the SHARAD radargram with noise, $s$ is the ideal subsurface reflection, and $n$ is the noise. The PID method produces an estimated total noise instance integrated over iteration j:

$$\tilde{n}=\lambda {\displaystyle {\sum }_j{n}_j}$$

(2)

However, the direct application of the PID method to the SHARAD data cannot produce satisfactory results with clear subsurface reflections. As shown in Figure 1, the PID method with various noise sigma values blurs visible subsurface reflections or completely erases weak subsurface reflections whose SNRs are only slightly above the noise floor. In this case, it is difficult to precisely estimate the noise $n$ (usually the estimation is conservative), which generates an overestimated $\tilde{n}$. As a result, the visible signal is blurred, or the weak signal is erased from the denoised results. Here, we propose the recursive subsurface boundary enhancement (RSE) method to improve weak reflections away from intense noise.

First, we improve the lateral continuity of the subsurface reflections by horizontal moving average, in which a limited number of traces (e.g., 3–5) is suggested [22]. Stacking along 1–3 rows of the radargram could improve the clarity of the thick layers along the vertical direction, but this would potentially degrade the vertical resolution of the thin layers.

Generally, radar reflections from subsurface layers are laterally continuous; by contrast, the noise is diffusive and chaotic. On the basis of the differential attributes of signal and noise, we extract reflections with nearly horizontal reflections of the un-denoised radargram by calculating the clustered local maximum in the time domain of the denoised radargram. To compensate for the leaked signal within every step of PID, we overlap the scaled power of the extracted reflections to the original data, which can enhance the weak horizontal reflections. A recursive procedure of these steps helps to enhance the weak reflections further, as follows:

$$s_{i+1}=s_i+c{\overline{s}}_i,$$

(3)

where ${\overline{s}}_i$ is the local maximum of the denoised radargram in the i-th iteration, and c is defined as the reflection enhancement coefficient. The enhanced lateral reflections of the radargram are taken as the input of the PID step in the next iteration; meanwhile, we increase the noise variance of the PID method [20] to roughly meet the improved signal-to-noise ratio after each iteration. The flowchart of the RSE method is shown in Figure 2 and Algorithm 1.

Algorithm 1. Recursive subsurface boundary enhancement method *.

Input: SHARAD radargram d.

Output: Denoised radargram r.

Initialization: enhancement coefficient c = 10%

and noise variance estimation σi = 10, 50, 80, 100.

1: for k = 1 to 3 do

2:   Vertical stacking on d at each adjacent k pixels: dk

3:   for j = 1 to 3 do

4:     Lateral moving average on dk at adjacent j pixels: s0

5:     Denoise $s_0$ using PID method with σ0:r0

6:     for i = 1 to 3 do

7:       Retain local peaks in each trace in ri-1:${\overline{s}}_{i-1}$

8:       Superpose $c{\overline{s}}_{i-1}$ to $s_{i-1}$:$s_i=s_{i-1}+c{\overline{s}}_{i-1}$

9:       Denoise $s_i$ using PID method with σi:ri

10:     end for

11:     return r

12:   end for

13: end for

* k, j and i are the indices of the iteration.

3. Numerical Examples

We performed numerical experiments based on synthetic data and SHARAD radargrams. The synthetic data contain various geological features, such as closely spaced layers, discontinues, and folds. The SHARAD radargrams contain rich subsurface reflectors and clutter. A comparison of the denoising results of the two models using different methods suggests the outperformance of the RSE method.

3.1. Synthetic Model

The size of the synthetic model is 500 pixels × 252 pixels, corresponding to a radargram with a 320 km × 9.45 μs size using a conversion of 460 km horizontally and 0.0375 μs vertically for one cell (Figure 3a). Two sections generate the model. One is the signal section (Figure 3b), and the other is the noise section. The signal section includes some laterally varying layers, a group of delicate layers, and some separated flat layers from top to bottom. A simulated radar pulse generates the layers. Note that we added some gaps in subsurface reflections to test the performance of the RSE method for recovering lateral continuity, with the gap width ranging from 1 to 5 columns from left to right, respectively (Figure 3b). The noise section is obtained from SHARAD radargram (orbit 224,401) using a window with the exact size of the model to intercept the background noise above the surface reflections. The simulated signal power is balanced with the average background noise level to test the performance of the RSE method when the SHARAD background noise is around the same level as the subsurface reflections.

Figure 4 shows the denoised results obtained by median filtering, the PID method, and the RSE method with different parameters. Median filtering produces the lowest-quality filtering results, with the delicate layers blurred (Figure 4a–c) and intense noise remaining (Figure 5a–c). With a low noise variance (e.g., 10) for the PID method, the background noise cannot be effectively reduced (Figure 4d). With a higher noise variance (e.g., 50 or 100), the PID method works better by presenting more continuous structures, but the subsurface reflections are blurred, with a poor vertical resolution of delicate layers. Therefore, for any noise variance listed, the PID method cannot recover the fine local layers in the model (Figure 3b).

By contrast, the results obtained by the proposed method have a much better vertical resolution of delicate layers while retaining the general shape of all layers. Additionally, horizontal smoothing among 3 or 5 pixels can improve the lateral continuity and partially compensate for the gaps. The possible drawback of the RSE method is that it introduces false reflections originating from random background noise. The false reflections become more evident as horizontal smoothing is performed with more pixels, while they typically have weaker power and less lateral continuity compared with real signals. Generally, the proposed method can better filter out SHARAD noise with a higher signal-to-noise ratio. It can obtain a higher vertical resolution by presenting all delicate local layers (Figure 4g–i) compared with the results obtained using the PID method (Figure 4d–f). As shown in Figure 5, the residual error of the proposed method is much smaller than that of the other two methods for all groups of parameters listed. The RSE method could enhance the reflection power while distinguishing the closely spaced layers. Comparing Figure 5a–c with Figure 5g–i, the gaps in the denoising results are relatively identical. This suggests that the median filtering method and the moving average have similar performance in recovering lateral continuity.

Figure 6 further shows the waves of the denoised results shown in Figure 4. The enhanced subsurface boundaries are discernible and highly match those of the theoretical models. This verifies the feasibility of the proposed method in extracting weak signals from a noisy background, even for models with gentle topography and lateral structure variations. Note that some artificial local peaks in the waveforms of the input data, denoted as “1” with arrows, are synthetic structures that can be suppressed well by the proposed method. For weak reflections in the waveforms of the input data (denoted as “2” with arrows), the proposed method can also enhance their energy correctly. While in some areas without evident reflections (marked as “3” with arrows), the proposed method’s results present some weak artifacts, but the energy is minimal.

3.2. SHARAD Radargrams

Figure 7 shows the results obtained by the proposed method for SHARAD radargram (Figure 1a). The RSE method can extract weak reflections that are originally difficult to distinguish from the strongly noisy background using traditional methods (Figure 1). With the increases in the enhancement coefficient c and the moving average window with size H, more details of the subsurface reflections can be presented; however, some artifacts appear when using a larger c, as shown in Figure 7e, especially over the surface reflections. According to Equation (3), the signal superposed by each structural enhancement inevitably contains some background noise, so c should not be too large; otherwise, the artifacts would be boosted. Our tests show that an enhancement coefficient of 5–10% can produce acceptable results with a high SNR and low artifacts. Under a small c, a larger H can present better continuity horizontally. Therefore, we suggest using H = 5 and c = 10%, which deliver the best results among all listed, as shown in Figure 7c.

We further tested two sections of SHARAD radargram 2758501 (Figure 8) across Promethei Lingua, where closely spaced delicate subsurface structures exist. Section A is located in the west of Australe Sulci, with a high elevation and significant terrain variation. Subsurface layers interfere with strong surface clutter. In the upper left corner, there are about five blurred layers [23]. Each layer covers more than 2 pixels in the vertical direction (i.e., >2 range cells), and their thickness grows with depth. The bottom layer covers 6–7 range cells. The diffuse echo layers are visible but problematic to distinguish from the noisy background. Under the shallow layer and at the bottom of section A, there is a low reflection zone with an unclear boundary that is highly polluted by the background noise [23,24]. Section B is located at the overlap region of Promethei Lingua and Australe Sulci, with a slight topographic relief, a gentle slope, and subparallel but intermittent subsurface layers. According to the number of clustered layers and sedimentation pattern of this area, we can divide these layers into four units: units 1 to 3 contain several fine focused layers, and unit 4 contains two blurred layers. Two types of unconformities can be identified in this area [25,26]. In general, the subsurface echo power of this area is above the mean level of the background noise, but the subsurface reflections are not continuous laterally. This might be attributed to the operating mode of the instrument or the phase distortion of the signal, although after correction processing [11]. The vertical spacing of the focused layers is relatively close, and the nearest distinguishable range cells of two adjacent layers can only be about 1 to 2 pixels.

The denoised results of section A are shown in Figure 9. All blurred layers in the upper-left corner of the unit are recovered well. The most continuous layers are presented in Figure 9f. With an increasing number of pixels in the vertical stacking, the overall brightness of the results is continuously improved, rendering clear horizontal subsurface boundaries. This indicates that those subsurface reflections are unrelated to random scattering relevant to the permittivity and surface roughness [11]. The background noise is filtered out for the two low-reflection zones, and the boundaries are more prominent.

The denoised results of section B are shown in Figure 10. The adjacent subparallel layers are successfully separated, and each layer is more distinguishable from the other. According to Figure 10a, seven layers can be identified in unit 1, and 14 layers can be identified in unit 2. Although the moving average improves the lateral continuity of the reflections, it also introduces some apparent artifacts in areas dominated by background noise (see H1–H3 in Figure 10), which may cause misinterpretation. However, a comparison between Figure 10d–f and Figure 10a–c illustrates that vertical stacking can effectively eliminate the artifacts caused by the moving average, despite a slightly reduced vertical resolution. Four blurred layers in unit 4 can be identified in Figure 10e–h with vertical stacking of 2 or 3 pixels, which are unclear in Figure 10a–d; meanwhile, the vertical stacking reduces the vertical resolution of shallow delicate layers that can be well identified in Figure 10a–c. Therefore, vertical stacking is effective for deep reflections but harmful to superficial delicate layers; additionally, the horizontal moving average causes artifacts for large travel times but helps extract weak fine layers for short travel times.

Figure 11, Figure 12 and Figure 13 present the details of local segments in Figure 7, Figure 9 and Figure 10, respectively. Clear horizontal reflections are extracted from the noisy background for various reflection features. The waveforms further show locally dominant reflections consistent in lateral directions, as shown in the adjacent panels of the denoised results. Additionally, almost no signal power appears before the surface reflections, which indicates that the extracted clear reflections are signals rather than artifacts.

4. Discussion

The basic idea of the RSE method is to gradually enhance weak subsurface boundaries using a recursive scheme of denoising, moving average, and partial stacking. Here, we used the PID as the core of the whole framework to obtain robust noise estimation. Note that the PID method, like other common image filtering methods, assumes an additive white Gaussian noise (AWGN) model. By contrast, common noise in the SAR image usually does not follow a Gaussian distribution but rather a Rayleigh distribution. The clearly enhanced quality of the SHARAD radargram using the RSE method indicates that the SHARAD background noise potentially follows a Gaussian distribution. With the development of the denoising technique, new methods with better accuracy or faster computational efficiency can be directly applied to the RSE framework to further increase the resolution capabilities. The proposed method can be applied to various noisy data-containing structures, including general images and scientific observations. The primary deficiency of the RSE method is that the denoised results rely on parameter selection according to professional experience. The enhancement coefficient and pixel numbers of the moving average should not be too large in case of evident artifacts. The possible cause of the artifacts is that the moving average enhances the lateral continuity of both the signals and the background noise. Thus, the PID method could mistakenly estimate the noise because of its enhanced correlation. The artifacts shown in Figure 10 are consistent with those present in the synthetic model results (Figure 4g–i); by contrast, the absence of artifacts in Figure 9 section A, with the same denoising parameter, suggests that the RSE method is sensitive to the background noise level, with more evident artifacts under a higher background noise level.

SHARAD radargram 2758501 across the Promethei Lingula area (Figure 8 section B) most closely recorded the spaced reflections. The four different geological units present in Figure 8 box B give a misleading interpretation of the penetration depth of the SHARAD signal in the Promethei Lingula region, which seems deeper than that of the adjacent regions. The possible reasons for this include the tilting of the layers, the surface roughness, and the geologic structure or physical state of the materials [27]. This section has dense subsurface reflections with short spacing. In order to exclude possible clutters that can generate weak boundary features and confuse the interpretation of subsurface reflectors, we used MOLA data [28] to generate simulated cluttergrams [29] and compared them with SHARAD radargrams. The clutter simulation result shown in Figure 14 illustrates that these subsurface reflections are not caused by surface clutter. In addition, we assumed that the reflections in the Promethei Lingula region (Section B in Figure 8) are not attributed to multiple reflections. Studies have correlated SHARAD radargrams with optical images in which outcrops are observed in exposed troughs [25,26]. The layers in optical images are more closely spaced than those in radargrams. If the packets of layers were caused by multiple reflections, they would be more intensively spaced and have weaker echo power. Apart from SHARAD and MARSIS, two ground-penetrating radar payloads operate on Mars’s surface [30,31,32]. The observations of these two payloads will continuously increase our understanding of the local details of the subsurface structures of Mars. However, noise is unavoidable, according to the experiences of similar payloads on the near side and far side of the Moon (e.g., [33,34,35,36]). Thus, weak signal enhancement is also necessary to improve the signal-to-noise ratio and identify local fine sedimentary structures under the rover’s track.

5. Conclusions

Noise interferes with the interpretation of subsurface reflections of layers in the SHARAD radargram. To best eliminate noise while maintaining the most prominent reflection features, we propose the RSE method to enhance the weak subsurface reflections recursively, particularly for regions with nearly flat topography and subsurface layers. Taking the PID method as a possible implementation, the RSE method consists of several steps. First, we enhanced the lateral subsurface reflections by the moving average within a given window (3 to 5 pixels) along the horizontal direction. Second, we reduced the outlier pixels that have variations that are too large in amplitude compared with the surrounding pixels using the PID method with a slight estimated noise variance (a controlling parameter of the PID method) [19]. Third, we enhanced the subsurface structures by adding a small portion of the local peaks of the lateral reflections. Then, we denoised using the PID method with a larger noise variance. Next, we applied the last two steps recursively (usually 3–5 times) to further increase the signal-to-noise ratio. Finally, we obtained the denoised results with nearly horizontal subsurface reflections evidently enhanced. The numerical experiments on synthetic and SHARAD radargrams illustrate that the RSE method can effectively remove the noise and retain the layered structures with a high vertical resolution. Although the RSE method relies on selecting input parameters, it could enhance the weak reflections for both the shallow visible layers and the deep blurred layers, which are difficult to distinguish from the original SHARAD radargram using existing methods. The denoising results are almost noise-free, which might shed light on reliable geological interpretations of the evolutionary history of sedimentary structures on Mars. The proposed method can be applied to a wide range of radar data processing situations and many other scientific observations.