2. SAR processing for layer enhancement

The proposed approach differs from methods based on the focusing of point scatters (see Fig. 2), in which azimuth processing is applied to reconstruct single points (also discussed as migration). SAR focusing techniques represent the optimal solution to improve the resolution of single points and rough surfaces. Instead,the proposed method matches azimuth history of a sloping specular reflector, which is the change in phase and range that a sloping specular reflector would produce in sequential traces within the processing aperture. For each range bin, this is obtained by modifying the along-track focusing to introduce a single phase shift Δφ _j to each trace within an aperture, which corresponds to a given specular slope. As shown in the block scheme of Figure 3, this is an iterative approach that guarantees the phase shift Δφ _j matches the layer-slope phase history (which is distinct from the phase history of a point scatterer).

Fig. 2. Illustration of the geometry for the proposed SAR (a), typical geometry of a point scatter focusing algorithm (b).

Fig. 3. Block scheme of the proposed method for LO-SAR processing and estimation of layers slope.

Our azimuth processing scheme consists of the following steps: (i) sidelobes weighting, (ii) range compression, (iii) phase-shifted along-track coherent summation, (iv) best phase-shift selection and (v) multilooking. For step (i), we use a Blackmann–Harris window, broadening the main lobe of the impulse response function (~ 80%) and reducing the peak-to-sidelobe ratio (− 50 dB). In step (ii), the received signal s _r is convolved with the complex conjugate of the time-reversed version of the transmitted signal s _t. It can be computed for each j-th trace as the following cross-correlation expression

(1)$$ s_j^{\rm c} = \int_{0}^{+\infty} s(t)^{*}_{\rm t} s_{{\rm r}}(t+t')\, {\rm d}t', $$

where the superscript $^*$ denotes the complex conjugate.

After pulse compression, standard unfocused SAR processing coherently sums each trace recorded over the synthetic aperture length. The process is applied to each range bin in the same way as

(2)$$ s^{\rm s} = \sum_{\,j = 1}^N \vert s_j^{\rm c} \vert \cdot \exp{\left(i {\rm Arg}({s_j^{\rm c}}) \right)} $$

where $s_j^{\rm c}$ is the complex pulse compressed signal and N is the number of range lines included in the synthetic aperture, and where | · | denotes the modulus and Arg( · ) denotes the phase of a complex number.

For unfocused SAR processing, coherent summation increases the SNR of flat specular layers. The same approach applied to adjacent echoes received from sloping layers, which are common subsurface features, can reduce the SNR or cause layers to disappear (Holschuh and others, Reference Holschuh, Christianson and Anandakrishnan2014). To address this issue we modified the formulation in Eq. (2) to include a phase shift Δφ _j that compensates the path length difference by shifting each received signal phase before integration. Accounting for the phase shift, integration over the synthetic aperture now yields

(3)$$ s^{\rm s} = \sum_{\,j = 1}^N \vert s_j^{\rm c} \vert \cdot \exp{\left(i{\rm Arg}({s_j^{\rm c}}) + i\Delta \varphi_j\right)}, $$

where Δφ _j is the optimal phase shift associated to the j-th range line's aperture. Finally, multilooking is applied to the along-track processed data in order to reduce the speckle noise through incoherent averaging. We choose a rectangular windowing for multilooking (3 by 2 pixels, in along- and across-track) following Castelletti and others (Reference Castelletti, Schroeder, Hensley and Bruzzone2017).

3. Automatic layer slope estimation

In this section, we demonstrate how the local slope of englacial layers can be automatically estimated as a byproduct of the improved processing algorithm discussed above. To this end, consider a platform flying in the x direction at a constant height h above the ground, and a single internal layer at depth d(x) below the ground and locally inclined at an angle θ with respect to the horizontal (Fig. 4). Accounting for the difference in refraction index between ice and air, the two-way travel time of the radar signal emitted when the aircraft is at x = x ₀ reads

(4)$$t_2(x_0) = 2\displaystyle{{r_an_a} \over c} + 2\displaystyle{{r_i(x_0)n_i} \over c},$$

where c is the speed of light, n is the refraction index, and subscripts a,i indicate air and ice respectively. We note that the beam paths r _a, r _i depend only on the material properties and the geometry of the problem, and are related to the geometric quantities h, d and θ through

(5)$$r_i = \displaystyle{{d(x)} \over {\cos (\theta )}},\quad r_a = \displaystyle{h \over {\cos (\alpha _a)}},\quad \displaystyle{{n_a} \over {n_i}} = \displaystyle{{\sin (\alpha _i)} \over {\sin (\alpha _a)}},\quad \theta = \alpha _i,$$

where the latter relationship is Snell's law. It thus follows that the two-way travel time is also a known function of the problem geometry. The notion of phase shift arises when comparing the two-way travel time at different along-track coordinates x. Following Figure 4, we now look at the dashed beam path. For a given posting Δx, and assuming Δx small compared to the scale over which the layer exhibits spatial structure so as to be able to consider the layer as locally linear, the depth of the layer at the new location x + Δx will be

(6)$$ d(x_0 + \Delta x) = d_0 - \Delta x \sin(\theta)\cos(\theta). $$

where θ correspond to the layer slope. The beam path inside the ice then becomes

(7)$$ r_i(x_0 + \Delta_x) = r_i(x_0) - \Delta x \sin(\theta), $$

leading to a difference in two-way travel time between the two adjacent locations

(8)$$\Delta t_2 = t_2(x_0 + \Delta _x)-t_2(x_0) = -\displaystyle{{2\Delta x{\mkern 1mu} n_i\sin (\theta )} \over c},$$

which we note is purely a function of the local slope of the layer (for given material properties). Finally, the phase shift Δφ between two adjacent radar tracks can be expressed as a function of Δt ₂ as

(9)$$\Delta \varphi = -\displaystyle{{4\pi f\Delta x{\mkern 1mu} n_i\sin (\theta )} \over c}$$

where f is the center frequency of the radar. Equation (9) is the relationship between phase shift Δφ and layer slope θ; with Δφ for each pixel of the radargram provided from the processing, Eq. (9) can be used to estimate the layer slope as

(10)$$\theta _{\rm s} = \arcsin \left( {-\displaystyle{{c\Delta \varphi } \over {4\pi fn_i\Delta x}}} \right),$$

where θ_s is an estimated layer slope at the pixel level. This approach results in a single unambiguous slope value estimated for the pixel resulting from (and centered within) a processing aperture. As a result, the optimal Δφ is applied between each successive trace within the aperture so that larger – and opposite signed – total phase shifts are applied at the edges of the aperture.

The uncertainly in this slope estimate is determined by the Doppler resolution Δf _D of our azimuth processing, which is given by

(11)$$\Delta f_{\rm D} = \displaystyle{1 \over T} = \displaystyle{1 \over {Lv}}$$

where T is the dwell time, L is the processing aperture and v is the aircraft speed (Raney, Reference Raney1998). This allows us to resolve a sloping specular layer, with a Doppler frequency f _D of

(12)$$f_{\rm D} = \displaystyle{{2v} \over \lambda }\sin \alpha _a = \displaystyle{{2v} \over \lambda }n_i\sin \theta $$

(Peters and others, Reference Peters, Blankenship, Carter, Kempf, Young and Holt2007), with a slope resolution ΔΘ of

(13)$$\Delta \theta = \Delta f_{\rm D}\left( {\displaystyle{\delta \over {\delta \theta }}f_{\rm D}} \right) = \displaystyle{\lambda \over {2Ln_i\cos \theta }}$$

as shown in Figure 5.

Fig. 4. Schematic of the idealized setup used to derive the analytical relationship between layer slope and phase shift.

Fig. 5. Estimated slope resolution of our LO-SAR processor as a function of layer slope and processing aperture.