Benefits of Coherent Large Beamwidth Processing of Radio-Echo Sounding Data

Coherent processing of radio echo sounding data of polar ice sheets is known to provide indication of bedrock properties and detection of internal layers. We investigate the benefits of coherent processing of a large azimuth beamwidth to retrieve and characterize the orientation and angular backscattering properties of the surface and subsurface features. MCoRDS data acquired over two distinct test areas in Greenland are used to demonstrate the specular backscattering properties of the ice surface and of the internal layers, as well as the much wider angular response of the bedrock. The coupling of internal layers’ 5 orientation with the bed topography is shown to increase with depth. Spectral filtering can be used to increase the SNR of the internal layers and for mitigating the surface multiple.


Introduction
Radio-echo sounding (RES) is a well established technique for remotely measuring the thickness of ice sheets.The use of synthetic aperture radar (SAR) focusing allows to improve gain and azimuth resolution of the echograms.Overall, state-of-theart SAR processing offers information about the spatial properties of the ice sheet and the strength of the response, which is used to determine ice thickness, internal layers' orientation and bedrock conditions, i.e. presence or absence of water.Several SAR algorithms were offered for focusing RES data, among them 1-D matched filtering (Legarsky et al., 2001), the ω − k migration (Gogineni et al., 2001), 2-D matched filtering (Heliere et al., 2007;Peters et al., 2007), and multilook time-domain back-projection (Mishra et al., 2016).Additionally, Holschuh et al. (2014) offer a method for improving SAR focusing of internal layers by introducing a correction of attenuation, migration and radial spreading for the returns from tilted internal layers.
Previous studies of angular backscattering properties of the ice sheet and bed are performed in (Jezek et al., 2009;Smith et al., 2010;Schroeder et al., 2013;MacGregor et al., 2015).Jezek et al. (2009) offer a technique for studying the backscattering properties of the ice sheet and bed using a special subaperture SAR approach.The authors study the dependency of the surface and bed return power on the incidence angle, the effect of the surface slope on the surface return power, they show that the response of the internal layers is specular, and propose incoherent presumming of subapertures to improve the SNR of internal layers.Smith et al. (2010) estimate an optimal value for the SAR beamwidth based on the bedrock signal-to-noise ratio (SNR).Schroeder et al. (2013) offer an approach for detecting the presence of subglacial water at the bed based on its angular  2016) introduce a novel approach for SAR processing of MCoRDS data, where the processing chain generates a number of SAR echograms, each corresponding to a particular incidence angle in along-track.A subset of the echograms with the highest SNR is then selected for further processing.The authors aim at improving the SNR of the weak bed echos in outlet glaciers and perform no further analysis of the backscattering properties of the ice-sheet and bed.
In this paper we introduce a new flexible technique to analyze the angular backscattering properties of the ice-sheet and bed, which can be applied to previously conventionally SAR focused echograms.Better understanding of those properties allows us to offers novel strategies for improving internal layers' and possibly bed SNR, to mitigate the surface multiple return, and to train sparsyfing dictionaries for model-based cross-track focusing methods such as (Wu et al., 2011;Heister and Scheiber, 2016).This paper begins with a description of the employed SAR focusing algorithm for RES data in Sect. 2. After that we introduce the technique for analyzing angular backscattering properties of the ice-sheet and bed in Sect.3. In Sect. 4 we analyze the results of processing two RES datatakes collected by the Center for Remote Control of Ice Sheets (CReSIS), Kansas, USA using their Multichannel Coherent Radar Depth Sounder (MCoRDS) during the Greenland campaign in 2008 (Gogineni, 2012).Based on the results of Sect.4, we discuss and demonstrate approaches for improving internal layers' visibility and for mitigating the surface multiple in Sect. 5. Finally, summary and conclusions are given.

SAR focusing
We perform SAR focusing of RES data using a modification of the range-Doppler algorithm.The processing is done in overlapping azimuthal blocks with each block processed as described in Algorithm 1.For each block we assume the platform to fly with a constant velocity v, the ice surface to have a constant along-track slope ψ, and the ice sheet to have a constant refractive index n ice = 1.78.We also assume that the electromagnetic wave propagation obeys Snell's law for a two-layer airice model shown in Fig. 1.The number of azimuthal samples in each block is selected to satisfy at least twice the desired SAR beamwidth of ∆θ = 30 • .We additionally assume that the azimuthal antenna pattern is broad enough so that its variation for incidence angles in the interval θ = [−15 • , 15 • ] can be safely ignored.
We now describe the inputs for Algorithm 1 using the notation where τ denotes range time, f τ denotes range frequency, η denotes azimuth time, and f η denotes azimuth frequency.The matched filter for range compression H RC is a complex conjugate of the Fourier transform of the transmitted signal weighted by the Hann function.The motion compensation filter H MOCO only corrects for a vertical component of the platform's deviation from a linear reference track in range frequency domain.where the geometric lengths that the electromagnetic wave travels in air and ice are (3) Both (Eq.2) and (Eq.3) depend on an unknown location of the refraction point s, which is a function of time η and depth d.The refraction point s can be found by solving a fourth-order polynomial equation (Heliere et al., 2007;Scheiber et al., 2008) or, more efficiently, by using Newton's optimization method, which iteratively finds s that minimizes (Eq. 1) with the following update rule at (i + 1)-th iteration where we initialize the refraction point with s 0 = 0.
Knowing s we calculate the phase shift of the received signal with respect to the time η = 0 when the platform crosses the origin of the x-axis where λ 0 is the wavelength of the transmitted wave in vacuum.
The Doppler frequency shift of the received signal is proportional to the derivative of Eq. ( 5) in time Knowing Eq. ( 6) for each depth d and azimuth position x, we compute the amount of range cell migration in range-Doppler Finally, we compute ∆φ(f η ) by interpolating Eq. ( 5) onto f η , and calculate the matched filter for SAR focusing H REF as We note that more precise and less restrictive SAR focusing algorithms for ice-sounder data exist, such as time-domain backprojection (Mishra et al., 2016); our choice of a particular approach described above is based on simplicity of implementation and its sufficiency for the subsequent analysis of the ice-sheet and bed angular backscattering properties.
Each subband is then accordingly zero-padded in azimuth so that all N subbands have the same size.After that an inverse azimuth Fourier transform is applied to each subband to get a set of N echograms I n , each containing returns coming predominantly from the corresponding incidence angle θ n .
The positions of ice-sheet features of interest such as surface, internal layers and bed are then manually selected from an echogram I incoh calculated as the incoherent sum We note that the ice sheet features can be tracked automatically, however, for a small amount of data we analyse in the paper, the manual selection is feasible.

Greenland MCoRDS data
We apply the approach presented in Sect. 3 to RES data collected by the CReSIS using their (MCoRDS) system (Gogineni, 2012).The main parameters of the radar and the acquisitions are summarized in Table 1.Two chirps with different duration were transmitted alternately on a pulse-to-pulse basis, with a 3 µs chirp intended to capture the surface and the shallow internal layer returns (shallow mode), and a 10 µs chirp intended to capture deeper internal layers and bed returns (deep mode).We employ the availability of multiple cross-track channels of MCoRDS to increase the SNR of nadir returns by combining SAR echograms of cross-track channels together using a conventional beamforming. We The full bandwidth echogram of track 1 is shown in Fig. 3  To study backscattering properties of internal layers we select a single internal layer with depth d ≈ 1870 m at azimuth x = 0 km.A deep layer is selected in order to avoid undesired contributions of the off-nadir surface returns.Figure 3(c) shows the internal layer's normalized power together with its slope, computed from (Eq. 10), drawn as a white line.We use bicubic interpolation to plot the figure.The internal layer response is specular, with θ max(I) corresponding to the layer's slope.A further insight into the behavior of internal layers' response is given in Fig. 3(d), where for each pixel of I incoh we colorcoded the incidence angle corresponding to the maximum intensity θ max(I) .Prior to plotting, we additionally applied a median filter of size (5, 5) and bicubic interpolation.The black lines on the figure correspond to the surface and bed return positions.The figure shows correlation between θ max(I) and the bed slope, with the blue and the red color appearing at azimuth positions with negative and positive bed slope correspondingly.Moreover, for a given azimuth position x 0 the absolute value of θ max(I) increases with depth, therefore, according to Fig. 3(c), the absolute value of internal layers' slope also increases with depth.This implies that the deeper the internal layer is located, the more its shape resembles the shape of the bed.
Figure 3(e) shows the normalized power of the bed response, where, prior to the normalization, we additionally compensate for the two-way propagation power loss of 2 dB/100 m.The incidence angle θ max (I) of the bed response varies in azimuth, overall the response is wide, meaning the bed is a rough surface for a radar with λ 0 = 2 m.
Figure 5 shows the dependency of the return power of the internal layer previously selected for Fig. 3(c) and bed for fixed azimuthal positions x = (3, 9, 46, 59) km.Those particular azimuth positions are selected to demonstrate the variety of shapes of reflective signatures for the bed and the persistent signature shape for the internal layer.We use quadratic interpolation to smoothen the signatures.
The full bandwidth echogram for track 2 is shown in Fig. 4(a), where we add a depth-dependent amplitude ramp of 2 dB/100 m.
Here the bed topography varies stronger as compared to track 1, the internal layers are visible close to the bed with gaps appearing at azimuth positions where the absolute value of bed slope is the highest; the surface multiple is also present in the    As expected, we observe larger color gradients for internal layers for track 2, whereas incidence angles of the surface multiple lie around θ n ≈ 0 • in white, corresponding to the ice surface.
The normalized power of the bed response for track 2 is shown in Fig. 4(e).
In Fig. 5 we compare the responses of the previously selected internal layer and the bed for fixed azimuth positions x = (3, 25, 32, 54) km.We again observe specular reflections from the internal layer and wider reflections from the bed, with 10 θ max (I) for the bed and the internal layer positively correlated for each selected position.We additionally compare distribution properties of the beamwidth of bed responses for both tracks in Fig. 7.According to the figure, the bed response of track 1 is overall wider than the one in track 2, suggesting that the bed in track 1 is in general rougher as compared to the bed in track 2.

Related applications
In this section we offer two straightforward applications of the results provided in Sect. 4.

5
First, the fact that an internal layer's response is narrow means that for a given depth and azimuth it contributes only to a small azimuthal frequency range in a SAR echogram spectrum.For an azimuthal spectrum of a small azimuthal block of a SAR echogram we see that at each depth internal layer's contributions are clustered around the frequency that corresponds to the internal layer's slope.In order to improve the SNR of the internal layers, for each depth we select a frequency of a maximum intensity f max(I) (d), and fit it into a piecewise linear regression to make the estimation of the true internal layers' frequency   more robust, after that we nullify the spectrum at frequencies lying outside the interval f max(I) (d) ± 0.05∆B az .We apply this method to track 1 SAR echogram, the results are shown in Fig. 8, where subsets of the SAR echogram before and after the processing are shown at the top and bottom correspondingly.
The procedure results in a 21.8% sharpness improvement in terms of intensity squared metric (Fienup and Miller, 2003), with the mean intensity of the echograms normalized prior to the comparison.The corresponding SNR improvement is expected to be in the order of 10 dB.
Second, according to the Fig. 4(d), the contribution of the surface multiple return can be mitigated by identifying and filtering out its contribution in the azimuth frequency domain, therefore revealing previously masked internal layers.This approach however only works in areas where the θ max(I) for internal layers and the surface multiple differ.We demonstrate the results of the method applied to a part of track 2 in Fig. 9. particular incidence angle in along-track direction.We estimated and compared scattering characteristics of ice surface, internal layers and bed for two datatakes in Greenland.For those datatakes the surface and internal layers have narrow response, which corresponds to a smooth specular surface, while the bed response is wide, which corresponds to a rough surface.The scattering properties carry information which can be used to estimate characteristics of the bed roughness (Fung and Eom, 1983), with the specular bed response indicating the likely presence of subglacial water at the bed (Schroeder et al., 2013).
Based on the scattering characteristics of internal layers, we offered a post-processing technique to improve their visibility.
By taking a small azimuthal block of a SAR echogram, within which the orientation of internal layers varies slightly in alongtrack, we observe that internal layer's contribution to the block's azimuthal spectrum is sparse and is clustered around the frequency corresponding to the internal layer's slope.This observation directly suggests a way to improve internal layers' SNR by keeping only those spectral components where the internal layers contribution is present.This post-processing technique can improve spatial tracking and interpretation of both straight and sloped internal layers.As a subject for further studies we suggest that denoising of all ice-sheet features in a SAR echogram is possible by finding a sparse representation of the echogram given a sparsyfing dictionary learned on patches with high SNR.
We also demonstrated a way to reduce the undesired contribution of the multiple surface return, which masks internal layers on corresponding depths.The reduction is possible when the surface multiple and the masked layer contributions come from different incidence angles, in which case they are separable in azimuthal frequency domain.
backscattering characteristics.The authors estimate the specularity of the bed returns by comparing their power contributions 1 The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-61Manuscript under review for journal The Cryosphere Discussion started: 18 April 2018 c Author(s) 2018.CC BY 4.0 License. in two HiCARS 60 MHz (Peters et al., 2007) SAR echograms with synthetic apertures of 700 m and 2000 m.MacGregor et al. (2015) introduce two new methods for estimating the slope of internal layers, among them the Doppler centroid method, which uses the fact, that internal layers' returns are highly specular.The authors use azimuth Fourier transform of short overlaping range-compressed RES data blocks, and derive the slope of internal layers from the wavenumber of the corresponding Doppler centroids.Mishra et al. (
The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-61Manuscript under review for journal The Cryosphere Discussion started: 18 April 2018 c Author(s) 2018.CC BY 4.0 License.3Multiple subbands processingIn order to analyze the dependency of the backscattering properties of the ice sheet and bed on the incidence angle, we divide the azimuth spectrum of an echogram into N overlapping subbands of beamwidth ∆θ sub = 2 • and an overlap between two adjacent subbands of 1 • .The central frequency of n = (1, N )-th subband, f 0 (θ n ), corresponds to the incidence angle of interest select two tracks both flown over Greenland in summer 2008, both approximately 70 km long.The track flown from the inland towards Jakobshavn glacier is referred to as track 1, the track flown over Southeast Greenland in North-East direction is referred to as track 2. The regions of interest, their topography and flight trajectories are shown in Fig. 2. Those particular datatakes are chosen to demonstrate how different bed topography affects the reflective properties of the internal layers; the bed in track 1 has depth varying in the interval d bed ∈ [2170 m, 3030 m] and slopes varying in the interval ψ bed ∈ [−35 • , 33 • ], The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-61Manuscript under review for journal The Cryosphere Discussion started: 18 April 2018 c Author(s) 2018.CC BY 4.0 License.
(a).To produce the figure we combine echograms of the shallow and deep modes, rebin the echogram in azimuth by a factor of 8, and add a depth-dependent amplitude ramp of 1.5 dB/100 m 5 to improve visibility of deep internal layers and the bed.The internal layers are visible until d ≈ 2 km.The gaps in internal layers' visibility occur at azimuth positions where the bed slope is the steepest.First, we investigate reflective properties of the ice surface.Figure3(b) shows the normalized reflectivity power of the surface as a function of incidence angle.The surface response is specular, with the incidence angle corresponding to the maximum intensity θ max(I) varying slowly in azimuth. 10

Figure 2 .
Figure 2. MCoRDS datatakes on a map.The map of Greenland is plotted using a stereographic projection with a central meridian of 41 • W and a central parallel of 72 • N. Isolines on tracks' maps correspond to a surface elevation change of 250 m.

Figure 5 .
Figure 5. Internal layer and bed returns for track 1.

Figure 4 Figure 4
Figure4(a) shows the normalized reflectivity power of the surface.The surface response is similar to the one for track 1, with higher variation of θ max(I) occurring starting from azimuth x = 65 km.Reflective properties of a single internal layer with depth d ≈ 440 m at azimuth x = 0 km are shown in Fig.4(b).Here we select a shallow layer because θ max(I) for deeper layers would lie outside the interval θ n ∈ [−14 • , 14 • ] previously selected in Sect.3. The incidence angle θ max(I) in Fig.4(c) varies stronger and more frequently as compared to that in Fig.3(c).5

Figure 6 .
Figure 6.Internal layer and bed returns for track 2.

Figure 8 .
Figure 8. Improvement of internal layers' visibility for track 1.

Figure 9 .
Figure 9. Mitigation of the surface multiple for track 2.

106
Summary and conclusionsIn this paper we offered a new approach to study scattering characteristics of ice-sheets which is based on division of a conventionally focused large beamwidth ice sounder SAR echogram into a set of subbands each of which corresponds to a 12 The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-61Manuscript under review for journal The Cryosphere Discussion started: 18 April 2018 c Author(s) 2018.CC BY 4.0 License.
the corresponding intervals for track 2 are d bed ∈ [640 m, 1970 m] and ψ bed ∈ [−62 • , 65 • ].We calculate slopes of the bed and internal layers as

Table 2 .
Beamwidth of the responses.