Introduction

The Intergovernmental Panel on Climate Change (IPCC) estimated that sea level would increase over a range of 18–59 cm over the next century. However, the IPCC also reported that models used to generate sea-level rise estimates did not include dynamic processes associated with rapid changes being observed in Greenland and Antarctica (Reference SolomonSolomon and others, 2007). Reference Pritchard, Arthern, Vaughan and EdwardsPritchard and others (2009) analyzed data from the Ice, Cloud and land Elevation Satellite (ICESat) and reported that both large ice sheets are losing mass, affirming the need for more realistic ice-sheet models. This need is well documented in a Scientific Committee on Antarctic Research (SCAR) report by Reference Van der VeenVan der Veen and others (2007). Finally, both bed topography and bed characteristics are identified as essential in developing improved models to predict ice-sheet response in a warming climate (Reference GogineniGogineni and others, 2007).

We conducted field experiments at Summit Camp (72.5783° N, 38.4596° W), Greenland, in July 2005 with a surface-based wideband synthetic aperture radar (SAR). The primary objective was to demonstrate the concept of wideband (120–300 MHz) SAR imaging of an ice-bed covered with very thick ice; our ultimate goal was to estimate the basal roughness and bed conditions from the backscattered signal. To accomplish this, a tracked vehicle with an eight-channel ice-penetrating radar was driven on a grid. The eight channels were then combined to generate left- and right-looking beams (Reference PadenPaden, 2006; Reference Allen, Paden, Dunson and GogineniAllen and others, 2008).

In this work, we apply tomographic techniques to data collected with our wideband and multiphase center cross-track array to generate fine-resolution bed topography from single-pass data. Application of tomographic techniques to produce three-dimensional (3-D) images of the subsurface (Reference Mast and JohanssonMast and Johansson, 1994; Reference Valle, Zanzi and RoccaValle and others, 1999; Reference Reigber and MoreiraReigber and Moreira, 2000) and to image vegetation and ice volume has been the subject of a few earlier studies (Reference Munson, O’Brien and JenkinsMunson and others, 1983; Jezek and others, 2008). However, none of these studies includes the application of parametric signal-processing algorithms (e.g. MUltiple SIgnal Classification (MUSIC)) to derive fine-resolution topography from single-pass data.

Data from closely spaced parallel tracks and perpendicular tracks (shown later in Fig. 7a) allow the results to be verified for self-consistency. Since the same spot on the ice–bed is imaged from multiple tracks, the elevation is independently measured multiple times. The measurements from two parallel tracks separated by 500 m and a track perpendicular to the two tracks are compared and found to be consistent, with a standard deviation of about 10 m.

Three-dimensional tomographic techniques produce a 3-D image of the ice and bed. For digital elevation model generation, we are interested in locating the ice–bed interface in the 3-D space. This work presents a simple, automated method for finding the ice–bed surface and removing point errors. Automation is critical for surface fitting due to the large number of data collected.

Equipment Description

The radar operates over the frequency range 120–300 MHz in the very high frequency (VHF) part of the spectrum. This frequency range is selected to obtain deep ice penetration with fine resolution. Radar configuration and pulse duration are designed to be fully programmable (Reference Paden, Allen, Gogineni, Jezek, Dahl-Jensen and LarsenPaden and others, 2005). The radar parameters and data collection configuration are given in Tables 1 and 2, respectively. The system transmits a linear chirp and pulse-compresses the received signals using a digital matched filter to achieve fine-range resolution. The received 10 μs signals from the eight antenna elements are multiplexed into two receivers using singlepole-quad-throw (SP4T) switches to map the base. One pulse of 1 μs duration is employed to map internal layers in the upper part of the ice, since the longer 10 μs pulse duration requires the receiver to be blanked when these close-in radar echoes are returning. Because of this, the effective pulse repetition frequency (PRF) is reduced from 6910 Hz to 1381 Hz. Only the 10 μs modes are used in this work. Also, the data are SAR-processed in two bands, 120–200 and 210–290 MHz, to facilitate a multispectral analysis of the backscattered signal. Only the 120–200 MHz band is used in this work due to its higher signal-to-noise ratio.

We designed the radar system for imaging the ice–bed with a cross-track array consisting of 16 phase centers. However, due to the failure of a transmit switch during installation in the field, we only collected data from eight phase-center positions. Reference PadenPaden (2006) contains a detailed analysis of the resolution and trade-offs for selection of the phase centers; only a review of the essential trade-offs is discussed here. A longer cross-track array typically yields higher-quality resolution. However, increasing the length of the cross-track array requires both a bigger sled than the 4 m wide sled used previously to collect data and more receive antennas. A wider sled increases the mechanical challenges, and more receive antennas increase the data storage rate and lower the signal-to-noise ratio, since the antennas must share the two receivers. For matched-filter processing and a flat surface, eight phase centers provide sufficient resolution for resolving a target about 500 m to the side so that the target on the opposite side (clutter) is suppressed. Having 16 phase centers allows the system to look even closer to nadir and still resolve the left and right sides of the track. One other trade-off that was considered was the use of an under-sampled array. The advantage of under-sampling is that fewer receive antennas are required, but this results in ambiguities in the data such that targets from different directions will produce similar measurements. Because of this disadvantage, we chose to fully sample the array with half-wavelength spacing to eliminate ambiguities.

The transmit antenna consists of two horn antennas arrayed in the H-plane. The feeds to this two-element array are phase-adjusted to form a nadir-looking beam. The two antennas allow two parallel transmit amplifiers to be used without the need for a lossy power combiner, and the antenna coupling positively affects the antenna return loss (Reference PadenPaden, 2006). The receive antennas consist of eight receive horn antennas arrayed in the E-plane. A diagram of the antenna network and the sled-mounted antenna networks is shown in Figure 1. The structure around the antennas (Fig. 1b) is designed to keep snow off the antennas to provide both consistent electrical performance and to avoid having too much snow weight on the sleds. The directional stability of the antennas was verified using the 3-D tomography algorithm to focus the internal layers. Because deep internal layers at Summit, Greenland, are flat to within less than one part per hundred (Reference PadenPaden, 2006), they should produce a specular response to a monostatic radar only in the nadir direction. The average measured scattering center from all strong internal layers along track 1 was less than half a resolution cell, corresponding to a 10 m horizontal displacement error.

Fig. 1. (a) Radar antenna array arrangement. (b) Transmit array, receive array, and radar system. Radar system housed in Compact PCI chassis is placed inside the cab of the tracked vehicle.

To reduce multiple reflections and improve radar sensitivity, low-noise amplifiers are integrated into the receive antennas. All other components in the receive chain are mounted on a printed-circuit board (PCB) with short electrical paths in the electomagnetic interference/radiofrequency interference (EMI/RFI)-shielded chassis shown in Figure 1b. The transmitter components, aside from the rack-mounted power amplifiers, are likewise mounted on a PCB in a separate compartment in the chassis. Attenuator pads are placed to suppress inter-device reflections where necessary. The transmitter and receiver designs are described in detail by Reference DunsonDunson (2006).

The receive channels are calibrated by measuring the time delay and phase through the entire system without the antennas using a calibration fixture. The calibration fixture is needed primarily to attenuate the high-power signal from the transmitter, so that the receiver will not saturate, and to terminate all unused channels. These signals, measured by the calibration fixture, are used after removing the effects of the calibration fixture to produce the matched filters for each receive channel used by the SAR processor.

The system also includes an internal calibration path without antenna feed cables. During field experiments, we observed <1° of phase change in this path. The antenna feed cables are 14 m long and are polyethylene foam cables with <10 ppm°C⁻¹ phase stability. The most extreme temperature change expected is 25°C. This corresponds to a maximum effective change of cable length of 3.5 mm, or <1°, at our operating frequency. Other considerations that help ensure minimal phase variation across receive channels are that 10 ppm °C⁻¹ is the worst case over the entire operating temperature range of the cable and that the same cable type is used for each receive path and should undergo phase changes similar to the other cables. Phase calibration is required since the eventual tomography algorithm needs to sense changes of phase on the order of 1/256 of a cycle, or 1.4°.

Sar Image Formation in Ice

The ice sheet is a layered medium of varying air and ice mixtures. At the top of the ice sheet, the newer snow has not compacted into ice and has an index of refraction around n _snow = 1.34 (Reference Bolzan and StrobelBolzan and Strobel, 1994). As the ice is compacted over time due to the overburden of snow and ice above it, the mixture turns to pure ice with an index of refraction of approximately n _ice = 1.78 (Reference Robin, Evans and BaileyRobin and others, 1969). Using geophysical measurements from the nearby boreholes (Greenland Icecore Project (GRIP) and Greenland Ice Sheet Project 2 (GISP2)), we approximate the layered media by discrete layers and estimate the permittivity, as described by Reference PadenPaden (2006). Only the real part is used in the propagation model and is plotted in Figure 2.

Fig. 2. Relative permittivity as a function of depth. Inset shows relative permittivity for the top 300 m of firn and ice.

Two assumptions are used here. First, the dielectric properties of the ice do not vary significantly away from the boreholes. This is reasonable, since the extent of the measurements is only tens of kilometers, and the time-averaged accumulation and weather patterns are unlikely to have varied much over this small area. The second assumption allows for thicker ice than that at the two boreholes by extrapolating the geophysical parameters for thicker ice using the bottom 100 m of each borehole. Since the thickest portion of the ice measured is only about 10% thicker (∼300 m) than where the boreholes (∼3000 m) were and the ice permittivity near the bottom is nearly uniform, this assumption likely produces minimal errors.

The permittivity profile shown in Figure 2 is used both in the SAR processing and the 3-D tomography. Frequency–wavenumber (f–k) migration using a layered-media approach is employed to focus the data (Reference Leuschen, Gogineni and TammanaLeuschen and others, 2000). Under the Born approximation, this algorithm is an exact matched filter for layered media, accounting for spherical spreading with refraction and propagation velocity changes. A 600 m synthetic aperture is used to focus the data, providing an along-track resolution of approximately 5 m at 3000 m depth. Data before and after migration correction are shown in Figure 3.

Fig. 3. Example of data before migration correction (a) and after migration correction (b). This example shows a section of bedrock with all eight channels focused toward nadir. The color bar indicates the relative signal power expressed in dB.

The permittivity profile in Figure 2 is also used to map the 3-D tomography’s direction-of-arrival estimate for each range bin to an off-nadir position. After SAR processing and 3-D tomography, the target’s along-track position, time delay and direction of arrival at the surface are known. To geocode the target, these parameters must be converted to along-track position, cross-track position and elevation, i.e. the right handed x-y-z coordinate system where x is along-track position, y is cross-track position pointing towards the left, and z is elevation.

A two-dimensional (2-D) table is constructed using a forward modeling approach. The forward model is supplied with an array of depths in ice and transmission angles into the ice, and the time delay and cross-track position are generated for each of these value pairs. The 2-D table is generated with fine enough spacing so that linear interpolation can be used to find a nearly exact inversion of this model. The forward model applies the ray-tracing assumption. This amounts to an application of Snell’s law at each boundary of the discrete layer model. Example results for 3000 m thick ice are shown in Figure 4.

Fig. 4. Time delay (a) and transmission angle (b) plotted versus off-nadir position.

3-D Tomography

Typical SAR focusing provides time-delay resolution through pulse compression and along-track resolution through azimuth processing, but the cross-track dimension is ambiguous. If the cross-track beam pattern is not fine enough, then left-side and right-side discrimination is not possible. This is the situation typically found with the depth-sounder radars that have been used on the ice sheets to map the base for decades (Reference Gogineni, Chuah, Allen, Jezek and MooreGogineni and others, 1998). However, with multiple-phase centers in the cross-track dimension, beam-forming methods can be applied to provide resolution in this dimension. If the resolution is fine enough, an estimate of the scatterer’s cross-track position can be made. Figure 5 shows the single-antenna cross-track beamwidth in black and the improved antenna-array cross-track beam-width in red.

Fig. 5. Illustration of the direction-of-arrival problem for surface height estimation. After SAR processing, the only ambiguous dimension is the cross-track angle, θ, so that each range shell only includes two sources of scattering when no layover is present.

To choose a beam-forming algorithm that will provide sufficient resolution to estimate the topography, the following was considered. First, the eight channels of the radar system do not provide fine resolution when combined using the standard matched filter, since the extent of the antenna array is only a few wavelengths and the range is on the order of thousands of wavelengths. Second, unless there is layover, there are typically only two scattering sources, as shown in Figure 5: one from the left side and one from the right side. This means that the number of array elements is four times greater than the expected number of scatterers or sources. Finally, the signal-to-noise ratio is good for small incidence angles. The matched filter does not provide improved resolution with greater signal-to-noise ratio. In fact, a bright response is likely to hide nearby weaker responses.

The problem-at-hand of estimating the topography can be viewed as a direction-of-arrival estimation problem. While there are a number of algorithms that can be used, this work combines the eight receive channels using the MUSIC algorithm (Reference SchmidtSchmidt, 1986). Reference Proakis and ManolakisProakis and Manolakis (1996) compare direction-of-arrival estimation algorithms in a chapter on spectral estimation; this chapter shows that MUSIC achieves the best performance. MUSIC works best when the number of sources is small compared to the number of measurements. To ensure this, each channel is separately migrated (SAR-focused) so that for a given pixel, power should only come from two directions: the left and right sides. This assumes that there is very little or no volume scattering from beneath or above the ice–bed interface. In this work, the two largest eigenvalues are considered the signal space (one for the left side and one for the right side), and the six smallest eigenvalues are considered the noise space. Example results of the pseudo-spectrum are shown in Figure 6. As expected, near the nadir return the MUSIC algorithm breaks down to some degree because of layover caused by bed roughness. However, as the constant range circle extends outward (i.e. later range bins), the chance for layover (ambiguity) decreases and two very distinct lines form in the spectral estimation image.

Fig. 6. MUSIC spectral estimation example results.

To produce the covariance matrix for MUSIC from which the eigenvectors are derived, five consecutive along-track positions with 5 m separation are used. The assumption is that the bedrock topography along this 25 m stretch is sufficiently smooth to assume that the bed scattering is stationary. For a given along-track index, x, and range index, ρ, there are eight measurements, one taken from each of the eight receive antennas. Let the vector represent these eight measurements. We then create a matrix from five consecutive range lines: . The autocorrelation matrix is then estimated by . The eigenvectors, , and eigenvalues for the matrix are found. Let be the ith eigenvector and also let the eigenvectors be sorted according to the size of their corresponding eigenvalues so that has the largest eigenvalue and has the smallest eigenvalue. The pseudo-spectrum is then given by where p = 2 is the number of expected targets, M = 8 is the number of samples per measurement, and s = [1 ej ^{2 πF} ej ^{2.2 πF} … e ^{j(M −1)2πF} ] ^T is the steering vector for the normalized spatial frequency, F, that we are estimating a return for.

The bed surface is generated in an automated way by taking the maximum range-bin return for each of the spatial frequency bins that the MUSIC algorithm produces. Considering Figure 6 as the visualization of the tomography output matrix, this is equivalent to taking the maximum for each column of the matrix. A total of 256 equally spaced spatial-frequency bins covering −0.5 to +0.5 are produced. 256 spatial-frequency bins sufficiently over-sample the result so that more bins only increase the computation time without improving results. With 256 spatial-frequency bins, the bin spacing corresponds to roughly 20–25 m basal resolution depending on range and incidence angle. The 20–25 m basal resolution is found using the ray-tracing method described in the previous section for each of the spatial frequency bins for 3000 m thick ice. To prevent occasional errors caused by the maximum finding routine, a 5 by 9 median filter is applied to the surface. If the difference between the median filter output and the surface is greater than a certain threshold (50 range bins), the value is considered to be in error and is replaced by the median filter output. About 0.6% of the pixels were replaced by this routine. After the point errors have been removed with the thresholding technique, a 3 by 3 median filter is applied to the data to reduce errors. The median filter further reduces the chance for point errors or outliers and has an averaging effect. The along-track resolution is 75 m after this step. The original cross-track pixel size was 20–25 m, but since the MUSIC algorithm provides resolution that is dependent on the signal-to-noise and clutter ratio, we can only state that the native resolution is degraded by 40–50 m after the median filter is applied. In other words, if the actual resolution was 50 m, it is closer to 100 m after the median filter is applied.

After the bed surface is found in terms of range bins and spatial-frequency bins, it is geocoded into a rectangular coordinate system. A 2-D table with the interpolation method described in the previous section is used. The range bin gives the primary index in the 2-D table, and the spatial frequency translated to surface transmission angle is used as the secondary index. The output from the table is depth (z coordinate) and cross-track position (y coordinate). The along-track, or x coordinate, is already known because the data have been SAR processed already.

The conversion from spatial wavenumber, ky, to surface transmission angle required to use the table is given by

where k = 2π/λ, λ = c/(n _snow f_c ), and c is the speed of light in a vacuum. Positive surface transmission angles will be assumed to mean the target is coming from the positive-y side of the platform (with a zero transmission angle being the nadir direction). The spatial wavenumber is related to the normalized spatial frequency (–0.5 to +0.5), F, that the MUSIC algorithm outputs by

where Δ _y is the antenna separation (24 in (∼61 cm)). Because the spacing is less than a half-wavelength at 160 MHz in snow, the bins for which k_y > k are not used.

Results and Discussion

Three tracks are considered in this work and are shown in Figure 7a: two parallel tracks spaced by approximately 500 m and a third perpendicular track. The primary track for comparison is 21 July 2005, sequence 6. The 22 July 2005 sequence-6 track nearly connects the GISP2 and GRIP boreholes; a no-vehicle research zone prevented the track from continuing all the way to GISP2. The 23 July 2005 sequence-20 track crosses the long east–west track near the GRIP borehole.

Fig. 7. (a) Tracks used in this work, where the circles indicate the start of each track. The origin is at Summit Camp, Greenland (72.5783° N, 38.4596° W). The GRIP (black circle) and GISP2 (black’x’) boreholes are also indicated. (b, c) Bed heights from two parallel tracks (1 and 2, respectively) separated by 500 m. (d) Relative magnitude height errors exceeding 10 m from two parallel tracks (1 and 2).

The 3-D tomography results for the two parallel tracks are shown in Figure 7b and c. Bed heights are with reference to the World Geodetic System 1984 (WGS84) ellipsoid. The along-track and cross-track dimensions are relative to the particular track displayed. Differences between these two registered results are given in Figure 7d, where along-track and cross-track dimensions are from track 1. The standard deviation of the difference is 9.1 m and the mean is 1.3 m. The 3-D tomography results for the two perpendicular tracks are shown in Figure 8a and b. Bed heights are also with reference to the WGS84 ellipsoid. Differences between these two registered results are given in Figure 8c. The standard deviation of the difference is 10.3 m and the mean is 0.8 m.

Fig. 8. (a, b) Bed heights from two perpendicular tracks (1 and 3, respectively). The GRIP borehole is marked by a circled asterisk. (c) Relative magnitude height errors exceeding 10 m for perpendicular tracks (1 and 3).

While the 2-D representations provide a consistent medium for comparison of point differences, a 3-D representation provides better visualization and interpretation of surface features. The full swaths for both parallel tracks, divided into two parts to keep the aspect ratio reasonable, are shown in Figure 9a–d, and the whole overlapped region from the perpendicular tracks is shown in Figure 10a and b for this purpose.

Fig. 9. Bed heights from parallel tracks (1 and 2) shown in two subsections to preserve aspect ratio. Information for one of the parallel tracks can be obtained by combining the images in (a) and (c); (a) is a continuation of (c). Information for the other parallel track can be obtained by combining the images in (b) and (d); (b) is a continuation of (d).

Fig. 10. Bed heights from two perpendicular tracks (1 and 3).

The black dot in Figure 10 is the location of the GRIP drill site. The results show a strong self-consistency in the bed elevation estimates, implying that the 3-D tomography, surface-finding routine, and geocoding are working well. The mean errors are plotted versus cross-track position for tracks 1 and 2 in Figure 11. The errors appear to be lower for small cross-track positions. This makes sense, because the beam-forming routines have the best resolution at nadir. Spatial frequency resolution degrades away from nadir, and the ground projection of the beam increases away from nadir. The errors also have a Gaussian-like distribution, as shown in the inset in Figure 11. Further spatial averaging in each individual result reduces the difference error between the two results, suggesting that a trade-off can be made between spatial resolution and error.

Fig. 11. Errors relative to the cross-track position of tracks 1 and 2. RMS: root-mean-square.

The only ground truth in the imaged area is the GRIP borehole, which is located near the perpendicular track crossing. GRIP ice-core drilling was stopped in silty ice close to bedrock at 3028.8 m (GRIP/GISP2, 1997). In Figure 8a and b, the borehole is indicated by a circled asterisk. The borehole appears to be located above a hillside. The slope is over 15° (60 m rise over 240 m run). The borehole estimate from track 1 is 3050 m, and the borehole estimate from track 3 is 3037 m. Since the borehole is directly beneath track 3, its resolution of the borehole is finer, so the error should be smaller. The error in track 3’s measurement may be due to the ice dielectric assumed. An error of 10 m at 3000 m depth translates to an error of 0.3% in the index of refraction, which is within the error bars of laboratory measurements of ice. The estimate from track 1 may be off due to errors in the surface-finding routine or a geocoding error: the hillside is shifted to the northeast as compared to track 3.

Conclusions

The use of a direction-of-arrival estimation routine (MUSIC) to provide left and right isolation and a bed height estimate has been shown to work successfully when coupled with an ice-sheet dielectric model. By comparing parallel and perpendicular tracks and showing that the difference is small, we have shown that the results are self-consistent. By comparing the results to the GRIP borehole (our only ground truth), we have established that the results also provide absolute accuracy. The spatial resolution in the along-track direction is 75 m. While the cross-track resolution is not theoretically determined, as it can be with the SAR resolution, analysis of the misregistration errors suggests that it may be similar to the along-track resolution. Finer cross-track resolution can be achieved by adding additional cross-track phase centers, as was intended by the use of a ping-pong transmit mode that failed during field installation.

We are developing a 16-element thinned antenna array for installation on the NASA P-3 aircraft for sounding and imaging of the Greenland ice sheet. The array can be used in ping-pong mode to collect data with 32 phase centers. The technique described in this paper can be applied to data to be collected for generating bed topography for many areas, including some of the fast-flowing outlet glaciers. SAR images produced from data collected with the large array can also be analyzed to determine basal conditions. Information about the bed topography and basal conditions is essential for developing realistic ice-sheet models to improve predictions of ice sheets’ contribution to sea-level rise in a warming climate. The technique presented in this paper, coupled with an advanced radar consisting of a large cross-track array, has the potential to obtain much-needed data for the development of next-generation ice-sheet models.