Introduction
The topography and nature of the materials near the bottom of temperate glaciers are of interest to studies of mass balance, surging, erosion and sediment transport. Direct methods, such as drilling, and traditional glacial remote-sensing (geophysical) techniques are incapable of spatially mapping the character of temperate-glacier bottom structures. State-of-the-art geophysics data-acquisition and -processing techniques have the potential for delineating many of the basal ice characteristics of interest to glaciologists. In this paper we present three-dimensional (3-D) radar array-processing methods that can characterize the bottom of a temperate glacier with a much higher degree of accuracy than traditional two-dimensional (2-D) radar techniques.
High-resolution, short-pulse, ground-penetrating radar (GPR) surveys of temperate glaciers are usually interpreted for ice depth using either unprocessed single profiles or standard 2-D migration. For ice thicknesses of <100 m, GPR profiles are usually obtained with systems operating in a 30–500 MHz bandwidth. In this frequency range, temperate glaciers are inherently high-clutter radar environments with coherent but unwanted noise (Reference Watts and EnglandWatts and England, 1976). In many cases, this clutter noise severely limits the effectiveness of radar surveys. The clutter is generated by diffractions from small englacial water and debris inclusions (Reference Smith and EvansSmith and Evans, 1972) and laterally propagating airwaves reflected from surface discontinuities (Reference CloughClough, 1976; Reference ArconeArcone, 1984). When one uses short-pulse GPR antennas having low directivity (Reference ArconeArcone, 1995), these forms of clutter production are highly problematic and lead to difficulties in record-section interpretation. In cases where good signal penetration is achieved, despite clutter, interesting englacial and basal ice reflections and diffractions often suggest features relevant to the hydrological and dynamic characteristics of temperate ice sheets. These events have been interpreted as conduits or water bodies (Reference Blindow and ThyssenBlindow and Thyssen, 1986; Reference Arcone, Lawson and DelaneyArcone and others, 1995), thrust planes (Reference Lønne and LauritsenLønne and Lauritsen, 1996; Reference Murray, Gooch and StuartMurray and others, 1997) and basal ice layers (Reference Arcone, Lawson and DelaneyArcone and others, 1995; Reference Murray, Gooch and StuartMurray and others, 1997). However, some of the prominent diffractions in these profiles suggest that out-of-plane scattering may have coincidentally placed these reflected signals in arrival-time proximity with the basal ice reflections. A short-duration wavelet with a tightly confined, high-gain radiation pattern would contain less signal clutter (Reference Watts and EnglandWatts and England, 1976; Reference SkolnikSkolnik, 1962) and have good resolution capabilities.
Several approaches have been used to reduce radar clutter in temperate-glacier surveys. For example, electrically large antennas produce records with less clutter by achieving a high gain and narrow beamwidth. However, large antennas produce extended wavelets with many oscillations that compromise vertical resolution. As shown by Reference Smith and EvansSmith and Evans (1972) and Reference Watts and EnglandWatts and England (1976), long-wavelength impulse radar at center frequencies of 2–8 MHz can overcome scattering losses and bring out bottom reflections. Recent trends in temperate-glacier radar surveys indicate that this is currently the consensus approach for overcoming the clutter problem (Reference Trabant, Molnia and PostTrabant and others, 1991; Reference Harrison, Echelmeyer, Chacho, Raymond and BenedictHarrison and others, 1994; Reference Nolan, Motyka, Echelmeyer and TrabantNolan and others, 1995; Reference EchelmeyerEchelmeyer and others, 1996; Reference Fountain and JacobelFountain and Jacobel, 1997; Reference Welch, Pfeffer, Harper and HumphreyWelch and others, 1998). Despite the increased penetration, the large wavelengths involved leave little but a bottom profile to be interpreted in the record. Consequently, it would be desirable to exceed 100 m of signal penetration depth at much higher frequencies, and to be able to do so with good lateral and vertical resolution.
Our radar methods are modifications of 3-D seismic survey practices (Reference ClaerboutClaerbout, 1976; Reference SchneiderSchneider, 1978; Reference StoltStolt, 1978) that have been adapted for the radiation characteristics of GPR antennas, and to meet the practical constraints of GPR surveys, namely, the need to keep data-collection and data-processing costs within reasonable bounds. Reference SchneiderSchneider (1978) applied the 3-D Kirchhoff method using stacked common midpoint seismic data. There have since been many examples of 3-D migration in reflection seismology. Typically, seismic studies collect vast datasets with complexly arranged source and receiver arrays (Reference StoneStone, 1994). The cost of such studies is roughly two orders of magnitude larger than the cost of this study. In the comparatively restrictive operational context of GPR surveys, direct application of 3-D seismic methods is difficult to justify. Recent seismic migration innovations have concentrated on the difficult task of combining wave-speed model determination with migration. Fortunately, in glacier studies, radar wave propagation speed is essentially constant, and large quantities of monostatic data can be collected quickly.
Multidimensional GPR processing of glacier radar data has been performed in a limited number of cases. Reference FisherFisher and others (1989) processed aerial Scott Polar Research Institute radar data (Reference Evans and SmithEvans and Smith, 1969) with 3-D array techniques to define the bottom topography of the central Greenland ice sheet. However, the narrow-bandwidth/long-duration radar signal and aerial data-collection method are not practical for small-scale radar profiling, and migration using bottom return-time picks is inappropriate for determining the characteristics of the bottom reflections and the structure of the basal ice zone. On Storglaciaren, Sweden, Reference Walford and KennettWalford and Kennett (1989) used a bi-static, 6 MHz short-pulse transmitter/receiver (Tx/Rx) system to construct a single, three-component synthetic-aperture array. They used a low-population (100 elements) array with a comparatively small aperture (50 m by 50 m) and applied time-domain Kirchhoff processing to image in a 3-D volume. The results noted anomalies that suggested cavernous water inclusions at depths of 75 m. Reference Hamran and AarholtHamran and Aarholt (1993) used wave-number-domain synthetic-aperture radar (SAR) (a method also called diffraction tomography) on subpolar glaciers. Diffuse scattered events in their records from water inclusions were imaged to depths of 200 m with 300–370 MHz radar. Although an ice–rock interface was interpreted to 200 m with a 5–20 MHz radar, it is difficult to judge how well this interface was defined. Recently, Reference Welch, Pfeffer, Harper and HumphreyWelch and others (1998) processed 5 MHz GPR profiles collected on a tightly spaced orthogonal antenna grid. They used standard 2-D seismic migration methods to map a reflection surface directly under the data-collection surface. The 2-D migrated sections were then migrated again in the third dimension by reapplication of the standard 2-D wavenumber method. They also showed that arrays allowed resolution of bed features smaller than the Fresnel zone. Reference Fountain and JacobelFountain and Jacobel (1997) also collected multiple GPR survey lines and applied 2-D migration. They used widely separated survey lines, predominantly oriented transverse to the direction of ice flow. They were able to define a bed-reflection topography for most of South Cascade Glacier, Washington, U.S.A., by applying standard 2-D frequency-wavenumber seismic migration methods. However, both Reference Fountain and JacobelFountain and Jacobel (1997) and Reference Welch, Pfeffer, Harper and HumphreyWelch and others (1998) note that record reflections can originate from out-of-plane regions that are not localizable by 2-D array processing. As will be shown, depth estimates from 2-D migration can contain substantial errors when applied to a steeply dipping valley glacier bed.
In the context of characterizing a temperate glacier with radar, the goals of this study are to demonstrate that practical 3-D data acquisition and array-processing methods remediate the associated problems of high clutter and shallow signal penetration, and that these methods allow true spatial localization of reflections. We accomplish this by synthesizing a dipole array, which forms a steerable, narrow, flashlight-like beam pattern that preserves the short-duration wavelet properties inherent in each element. The theoretical resolution attributes of our approach are reported in Reference Moran, Greenfield, Arcone and DelaneyMoran and others (2000). That paper uses synthetic data with high noise levels to show that the signal-to-noise-ratio (SNR) improves dramatically with higher migration dimensionality, and that discrete scatters can be accurately localized outside the array coverage region. In the present work, these methods are applied to real data and the performance improvements are discussed and analyzed. The 50 MHz dataset was collected on a grid on the northeast branch of Gulkana Glacier, located in the central Alaska Range (Fig. 1). At this site the valley walls were steep and easily approached and suggested a complex bottom topography that was nearly ideal for showing the reflection-localization capabilities of 3-D arrays. We supplemented our grid survey with 400 MHz profiles to help define our noise field, and 12 MHz profiles to verify the character of the bed dip.
Aerial photograph of Gulkana Glacier with schematic of GPR survey lines. The 100 m by 340 m, 50 MHz GPR array is the rectangular box slightly to the left and below the U.S. Geological Survey hut.
Radar Equipment
50 and 400 MHz system
Our data were collected using a commercial GPR system (also known as impulse, short-pulse, transient, carrier-free, sinusoidal or monocycle radar). The operation of this system has been described previously (Reference ArconeArcone, 1995; Reference Arcone, Lawson and DelaneyArcone and others, 1995). The 16-bit control unit sets the pulse-repetition rate, trace parameters and time-variable gain function. The 400 MHz antenna transducers were commercially manufactured. They are backshielded, resistively loaded dipoles with 15 cm separation. The transmitter excites an 8 W peak power pulse on the antenna. The 50 MHz dipole antennas were developed in-house using a string of resistors. These antennas were unshielded and separated by 2.9 m. When collecting 50 MHz data we used a transmitter with 800 W peak power. Both antenna systems were towed behind a snow machine and were polarized with their antenna axes transverse to the profile direction.
12 MHz system
The 12 MHz system is similar to other non-commercial constructions based on digitizing oscilloscopes (Reference Jacobel, Anderson and RiouxJacobel and others, 1988; Reference Jones, Narod and ClarkeJones and others, 1989). The receiver is an eight-bit Tektronix model 524A digital oscilloscope. The commercially available transmitter (Reference Narod and ClarkeNarod and Clarke 1994) generates high-power (2.5 kW) pulses at 512 repetitions per second in synchronous operation with the oscilloscope. The antennas are 10 m long resistively loaded dipoles arranged (and polarized) in a co-linear mode for easy towing. Their center-to-center separation was 20 m. The receiver antenna was connected to the oscilloscope with fiber optics to avoid cable interference with the antenna radiation.
Antenna patterns
Theoretical far-field antenna radiation patterns in ice for resistively loaded dipoles of finite size are given by Reference ArconeArcone (1995) and Reference Arcone, Lawson and DelaneyArcone and others (1995). The patterns apply to the far field regardless of antenna orientation (parallel/opposing or parallel/co-linear) so long as the separation is a small fraction of the far-field survey depth. The patterns are irregular in the vicinity of the critical angles of 34° (sin−1(1/n), where n = 1.78 is the refractive index of ice and the angle is measured from vertical).
In ice, far-field beam patterns form within 1 m range at 400 MHz, by 20 m at 50 MHz, and by 80 m at 12 MHz. The patterns in both the E-plane (the electric field within this plane is parallel to the axis of the dipole) and the perpendicular H-plane (magnetic field) have wide angles of sensitivity, but the H-plane patterns have much greater power near, and at, the critical angles. Therefore, our 12 MHz co-linear antennas had greater sensitivity transverse to the direction of profiling, whereas our 50 and 400 MHz antennas had greater sensitivity along the profile direction. The directivity (F(θ, ϕ)) of an infinitesimal dipole (Reference Engheta, Papas and ElachiEngheta and others, 1982; Reference ArconeArcone, 1995) is incorporated into our migration model because it approximates well the pattern for finite-size dipoles and it has a closed-form solution. However, the directivity achieved by a large aperture array is far narrower than F(θ, ϕ) and precludes concern as to the exact form of F(θ, ϕ) beyond a few degrees off the center axis of the array beam.
Analytical far-field expressions for interfacial dipole radiation have a 1/range dependency for the beam’s geometric spreading loss, and a null pattern along the surface. However, it is well known (Reference BrekhovskikhBrekhovskikh, 1960; Reference BuddenBudden, 1960; Reference AnnanAnnan, 1973; Reference CloughClough, 1976) that near-field 1/range2 dependency terms exist along the surface and produce both ground waves traveling at the ice speed and lateral waves traveling at the air speed. Characteristic of surface wave propagation, these laterally propagating airwaves continually refract energy into the ice at the critical angle, which is then prone to scattering from near-surface inclusions.
Site Description
Gulkana Glacier is a small temperate valley glacier on the western slope of the central Alaska Range. It is easily accessible by snowmobile from the Richardson Highway. Gulkana Glacier has been well studied (Reference Péwé and RegerPéwé and Reger, 1983), especially as part of a long-term mass-balance investigation (Reference March and TrabantMarch and Trabant, 1996). We operated in the upper 2 km of the main branch of Gulkana Glacier (Fig. 1), which is within the ablation zone (no firn cover). Most of the accumulation occurs in and above several tributary icefalls feeding into the main branch. Gravity (Reference Ostenso, Sellmann and PéwéOstenso and others, 1965) and point radar soundings at 5 MHz (unpublished data from D. C. Trabant, Geophysical Institute, Fairbanks, Alaska) suggested that maximum depths in the main branch are approximately 300 m. Reference EchelmeyerEchelmeyer and others (1996) report that the elevation of Gulkana Glacier in this area has changed little in 50 years (an average decrease of about 11 m, despite severe recession of the terminus). Their data, along with those of Reference SellmannSellmann (1962) and Reference RutterRutter (1965), provided a check on the relative elevations measured in this study. The snow cover was generally 2 m deep and appeared dry at the time of our surveys (early April).
In choosing the array site we sought a location with a safe, smooth surface abutting a valley wall where the bottom topography deepened gradually in the cross-valley direction. These conditions are desirable for demonstrating improvements in clutter suppression and depth of penetration. Reconnaissance data collected at the base of a rock slope near the U.S. Geological Survey hut (Fig. 1) showed a gradually deepening bottom return that was undetectable by about 120–140 m depth. The surface of the snow was smooth and apparently safe. The selected data-collection area was groomed with the snowmobile. This facilitated towing of the radar-operator sled and antennas. We did not expect till layers or large englacial conduits in upper reaches of the glacier. Also, the smooth, relatively flat surface of the glacier gave no indication that significant 3-D bed topography was present.
Data Acquisition
We established a 100 m by 340 m regularly spaced array grid using a laser theodolite. The 100 m wide x dimension was oriented in the down-glacier, southerly direction. The 340 m y dimension was oriented transverse to ice flow, in a roughly easterly direction. The western end (y = 0 m) of the array was closest to the hut. We located station markers at 10 m intervals in x and 20 m intervals in y. The maximum surface-elevation difference within the array grid was 18 m, with the dominant elevation gradient occurring along the y axis of the grid.
The 50 MHz data have a trace duration of 3400 ns (16-bit, 1024 samples/trace). The 400 MHz data have a trace duration of 400 ns (16-bit, 2048 samples/trace). The 12 MHz records have durations of 8192 ns (eight-bit, 1024 samples/trace). These sample intervals give about 6, 13 and 10 samples per wavelet cycle at 50, 400 and 12 MHz, respectively. Although the 50 MHz data are mildly undersampled compared with typical intervals of ten or more points per wavelet, the data are not aliased in the frequency domain. We applied a time-variable gain to the 50 and 400 MHz data during recording. In the 50 MHz array data, all discernible signals of interest are contained in the first 2000 ns of the records. In the 12 MHz survey, we applied a 32-fold stack to provide a recording rate of 16 traces per second. The 400 and 50 MHz data were stacked only two- to threefold during the position-justification process.
Within the array grid we collected 51 parallel survey lines using the 50 MHz antennas. Each line was separated by approximately 2 m and oriented along the y axis of the grid (west to east). In addition, we collected 400 MHz data on ten y-oriented parallel lines spaced 10 m apart. The 12 MHz profiles discussed in this paper were obtained along a 100 m wide longitudinal corridor beginning near the maximum height of the main branch of the glacier (Fig. 1). Positions for the 12 MHz profiles were obtained with a GPS system having 30 m accuracy.
Image degradation in 3-D migration is a complicated function of the distribution of antenna-position errors, wavelet center frequency and the focus position. Though our time-domain approach allows us to use spatial sample intervals larger than a wavelength, a conservative bound for the mean error in antenna positions would ideally be less than one-fifth of an in situ ice wavelength, or about 70 cm at 50 MHz. Fortunately, these error tolerances are mean values and need only be maintained across a spatial extent corresponding with the largest array apertures to be processed. For example, a single 3-D migration result may use several overlapping, 40 m by 60 m sub-arrays that “tile” over an area of 40 m by 300 m. It is only within each sub-array that the sensor-position error tolerance needs to be maintained.
We used snowmobiles to tow the recording system and antennas at a nearly constant speed with a high trace-acquisition rate. Electronic position markers were entered in the record headers and then digitally justified to provide the same number of traces between recorded marker intervals. In the 50 MHz array grid, the y-position accuracy was better than 20 cm. Control over the x location of each line was problematic because the slightly down-sloping glacier surface caused antenna sleds to track as much as 1 m off the desired x position. We estimate our mean x position error to be <30 cm. A more desirable antenna-positioning method could be made by incorporation of differential global positioning system (GPS) coordinates directly into digital record headers instead of simple markers. In such a system, individual antenna positions could be specified with errors of <5 cm.
Data Processing
Initial data processing using commercial software included slight low-pass filtering, background removal, position normalization and elevation corrections. The range gain used in recording the data was maintained during array processing. The raw data were justified to give ten traces per meter for both the 50 and 400 MHz datasets. The 12 MHz data were justified to give five traces per meter. Elevation corrections for the 12 MHz profiles were obtained using single reading non-differential GPS with a position accuracy of roughly 30 m.
Deconvolution
Multiple cycles in the raw 50 MHz wavelets, attributed to weak multiple reflections between the surface and snow-ice interface (Reference Arcone, Lawson and DelaneyArcone and others, 1995), were removed prior to migration using an optimum Wiener deconvolution filter (Reference YilmazYilmaz, 1987). The deconvolved data improved the temporal resolution of bottom reflections and englacial diffractions (Fig. 2). The deconvolved bed-reflection wavelet is characterized by three large, half-cycle, oscillations (negative, positive, negative). This bottom-reflection attribute is consistent throughout our dataset. The parameters that produced the largest variations in deconvolved results were the shaping wavelet, the level of pre-whitening (random noise) and the segment of data used to design the filter coefficients. The deconvolution filter coefficients were calculated by forming an autocorrelation matrix from a raw bed-reflection wavelet (with a 106 ns period) and shaping this with a zero-lag boxcar wavelet having an 11 ns period. The resulting coefficients are applied to all the 50 MHz data. In several places the deconvolved data show closely spaced overlapping reflections. It is likely that these overlapping reflections result from multiple out-of-plane reflections arriving at nearly the same time.
Raw and deconvolved 50 MHz profiles along x = 11 m. The inset waveforms are taken from y = 150 m. The boxed wavelet is the bed reflection. The deconvolved data show significant broadening of wavelet bandwidth. This improved the migration results. The prominent bottom reflections show many subtle out-of-plane influences, such as intermittent layer-like overlapping Amplitude anomalies from out-of-plane reflections cause noise in migration results. The linear trends in the englacial clutter field are strong lateral waves.
Array processing
We imaged the subsurface using a slight modification of the time-domain Kirchhoff integral migration method given by Reference SchneiderSchneider (1978). The strength of the energy scattered from a subsurface point is given by
where
is the range-normalized single dipole radiation pattern (Reference Stutzman and ThieleStutzman and Thiele, 1981),
are the GPR surface observations,
is the desired image in the subsurface at R′, and
is the two-way travel time (retarded time) from the subsurface image point to the sensor point. We treated the exponents α and β as processing parameters because, currently, there is no theory that includes effects of transmitter and receiver radiation patterns or considers radiation pattern distortions ascribable to time-range gain. A systematic trial-and-error study on the 50 MHz data gave α = 2 and β = 1.5. It is likely that our values for these parameters are directly related to the time-range gain recorded with our data. Properties of this procedure for radar data are discussed by Reference Moran, Greenfield, Arcone and DelaneyMoran and others (2000).
Our implementation of Equation (1) follows a synthetic-aperture array approach. In a subsurface volume, a migration image is formed by evaluating Equation (1) at discrete voxels using data from a rectangular array of surface points, as shown in Figure 3. To image a y-z plane, a subset of the array data (sub-array) is first fixed with its center at some x and y position. The sub-array data are then used to image a portion of the y-z image plane. The sub-array is then moved, holding x fixed, to greater y, and an additional portion of the y-z plane is imaged. Generally, the two portions of the y-z plane have some overlap. In this overlap region the mean of the two results is taken. The center of the subarray is moved (along y) until the y-z image is completed. Sub-array apertures are generally about 40 m wide (x dimension) by 60 m long (y dimension), and overlapped by 40–60% for successive y positions. Figure 4 (explained in the Appendix) shows that these array dimensions result in a flashlight-like beamwidth of roughly 5° and that the resolution of the migration process is comparable to that of a delay and sum beamformer (Reference Capon, Greenfield and LacossCapon and others, 1969). Using overlap between sub-arrays reduces edge effects between imaged portions in the fully formed y-z image plane. To image an x-z plane the same procedure is used, with the y of the array center held constant as x is increased for successive sub-array positions (Fig. 3). Further discussion of our array processing is given in the Appendix.
Cartoon of our 3-D array and image-plane geometries. Data were collected continuously along a series of standard survey lines parallel to the y axis. A synthetic-aperture array is defined by a sequence of sub-arrays arranged parallel to the x or y axis. A migration depth image is formed in a set of vertically oriented planes.
2-D linear-array beam response vs focus position for an isotropic point reflector buried 100 m below the array center (d = 0). (a) Frequency-domain beam pattern and (b) migration reflector strength vs focus position. In both results the array lengths are given on each curve and the inter-element array spacing is 2 m. In all cases, the −6 dB beam width is <5 m.
Results
In cold-based glaciers, the ice–rock interface produces large radar reflections. However, there are several ways that reflections from the bottom can be lost. The most important are: (1) loss of energy to scattering from englacial inclusions, (2) masking by englacial clutter, and (3) a bottom topography that does not reflect rays to the data-coverage region. Although our data indicate that localized inclusions near the surface produced significant radar clutter, the 3-D array results show that loss of bottom reflections is principally caused by topographic changes in the bed surface. Our results also show significant improvement over 2-D migration methods.
Near-surface ice-clutter characterization
Near-surface ice inclusions and discontinuities, such as crevasses, produce radar clutter whose signature can extend deep into a record at arrival times similar to those of deep bed reflections. Aerial photographs (photo 87V2-186, Geodata Center, University of Alaska Fairbanks) from 1987 show that a complex system of narrow, crossing crevasses occurs in the area of our array grid. We use 400 MHz profiles to characterize the upper 33 m of the ice. They show a 2 m thick snow layer and clusters of hyperbolic point diffractions under the eastern two-thirds of our grid (Fig. 5); most of the hyperbolae reside in the upper 20 m of the ice. In general, these diffractors are associated with broad but modest depressions in the radar reflection horizon associated with the base of the snow layer. This would be expected when profiles cross narrow open crevasses at oblique angles. No snow-layer radar depressions are seen in the western end of the grid. At all depths, 2-D migration using the ice-propagation speed (ε r = 3.17), collapses the hyperbolae to points. This result implies that a significant fraction of our englacial noise originates from localized disturbances within the top 20 m of ice, and that the ice radar velocity is not impacted by the groomed snow surface. There is no evidence in the 400 MHz (or 50 MHz) data of interconnection between diffractors across the grid (i.e. linear patterns) or of phase changes from voids or conduits (Reference Arcone, Lawson and DelaneyArcone and others, 1995). Both the 400 and 50 MHz data, the ice underlying the western third of our grid, are relatively free of clutter (Figure 5). This is, in part, due to the lack of open crevasses in this region. The few weak diffractions originating in this zone also show hyperbola geometries characteristic of ε r = 3.17. This clear ice extends to the bedrock, at about 30–38 m in this zone.
Typical section of a 400 MHz record (a), which shows densely clustered point-diffraction hyperbolas, all of which migrate at the ice velocity (ε r =3.17). They are generally distributed along 20–40 m wide swaths (b) and are likely to be associated with the crossing pattern of crevasses (b), as interpreted from aerial photography. Prominent scattering is not observed in the western third of the grid.
In contrast with the far-field diffractions seen at 400 MHz, a dominant noise feature in the 50 MHz data (e.g. Fig. 2) is manifested by lateral waves. This noise is composed largely of numerous resonant events having a linear move-out of approximately 0.15 m ns−1. This is nearly the two-way propagation speed of air. These returns can be explained as laterally propagating airwaves, which continually refract into the ice (Reference CloughClough, 1976; Reference ArconeArcone, 1984). The refracted airwaves are scattered back to the receiver by near-surface discontinuities. Analogous phenomena are observed in seismic surveys (Reference Larner, Chambers, Yang, Lynn and WaiLarner and others, 1983). Lateral waves scattered from linear discontinuities (such as crevasses) at oblique angles will exhibit move-out velocities that are faster than airwave speeds (Clarke and Bentley, 1993). These are not seen in the 50 MHz data. Therefore, we believe our 50 MHz clutter was caused primarily by volumetric scattering from the localized shallow inclusions seen at 400 MHz.
Delineation of bottom topography
2-D array results
Traditional 2-D migration attempts to position all reflections in a vertical plane directly below the survey line. A 2-D y-z plane depth image from our data is shown in Figure 6. In the coordinate system of the array grid, the line position was along x = 11 m. The image shows many large- and small-amplitude circular “smears” in the depth image, which result from a strong noise field (Reference YilmazYilmaz, 1987, p. 274) and high-amplitude late arrivals in the time series (Fig. 2). Our 3-D results will show that these high-amplitude late arrivals are most likely out-of-plane reflections and are not properly localized by summation along a 2-D hyperbolic time–space curve. Also apparent in Figure 6 is a large discontinuity in the bed reflection in the interval, 75 m < y < 130 m. This results from the large bed dip along the survey line and (as will be shown later) a substantive component of bed dip perpendicular to the survey line. In such cases reflected rays intersect the surface outside the survey line, and 2-D processing is sufficient to resolve or correctly position the bed reflections.
2-D migration result along x = 11 m. Large- and small-amplitude anomalies from out-of-plane rejections have not been properly localized, and have also become poorly resolved, archingsmears of energy. This result used six linear sub-arrays; each has a 60 m aperture with 600 elements. Each sub-array overlapped the next by 17%.
Exploding-reflector-ray theory can partially explain the poorly resolved bed reflections from steep in-plane dips in 2-D migration results. Rays originating from regions of steep dip will generate low-density ray patterns (shadow zones) on the surface and have longer path lengths. The longer travel paths imply larger scattering losses, and the high angles of ray incidence at the surface may lie outside major lobes of the dipole beam pattern. These factors cause recorded bed-reflection amplitudes to be lower, and noise levels higher. This can account for many of the 2-D imaging degradations seen in Figure 6. However, if there are out-of-plane dip components (particularly in regions of steep dip) then 2-D migrations will fail to resolve the bed reflection, even though the raw-data records show high signal levels.
3-D array results
The 3-D array processing allows correct spatial positioning of reflections in x, y and z. As we shall show, in some circumstances reflectors can be correctly positioned outside the data-coverage area. Depth images for a sequence of four parallel y-z planes are given in Figure 7a–d. The images are generated from a single 3-D array that combined ten parallel survey lines. The array has a 6 dB beamwidth of <6° in x and 3° in y. In each image plane there are large-amplitude reflections from normally oriented facets of the bed. The relatively undistorted wavelets (compared to bed-reflection wavelets in Figure 2) associated with the prominent reflections indicate that the array has been properly focused. Within any given slice, the brightest reflections abruptly diminish as the array focus shifts along y or x (between images). The sharp bed termination in both the longitudinal and transverse directions demonstrates that the array processing is spatially suppressing highly coherent, large-amplitude reflected bed energy. This characteristic is similar to very narrow, flashlight-like, theoretical beamformer results (Reference Capon, Greenfield and LacossCapon and others, 1969). Figure 7e is a composite that superimposes each major reflection onto a single projection. It shows that a nearly continuous bed surface has been mapped out.
Four (a–d) parallel y-z plane image slices generated from a single 3-D array centered along x = 11 m. The prominent reflections occurring at various positions demonstrate localization of out-of-plane normal reflections. (e) Composite of all four reflections projected onto a single y-z plane. The results used a sub-array aperture of 20 m (in x) by 40 m (in y) with 671 elements, each separated by 2m and 1 m in x and y, respectively. Two sub-arrays covered the data interval between 0 m ≤ y ≤ 120 m.
Figure 8 shows four parallel x-z image planes produced from a single 3-D array result. Each image plane shows a coherent reflector between 15 m ≤ x ≤ 45 m with depths of 108–122 m. The image plane directly under the center line of the array (y = 175 m; Fig. 8d) shows an abrupt loss of focus and amplitude at roughly x = 40 m at 120 m depth. Later discussion will show that the bed slope abruptly increases (in the down ice-flow direction) in this region. In expectation of observing reflections from the up-dip (lower x) direction, we extended the image-plane boundaries to x = −30 m. In Figure 8a (y = 145 m) there is a distinct normal reflection at x = −20 m at 95 m depth. The loss of reflector continuity between this event and that imaged at x = 20 m must be caused by a steeply dipping slope in this region. The depth images in Figure 8 also contain low-amplitude, poorly focused reflection smears that result from poor sensor positioning in the x direction, and edge effects produced by limited array coverage. Relative sensor-position error within a sub-array introduces a slight misalignment of bed-reflection arrival times in the hyperbolic time–space summation performed in the Kirchhoff integration. This results in in complete signal cancellations as the focus point is shifted away from a normal reflection point on the bed. Edge effects will be seen when there is insufficient array coverage for positioning a normal reflection. In this example, extending the array coverage farther in x would improve image quality.
A sequence of four parallel x-z plane depth images generated from a single 3-D array run. All the images show well-resolved bottom reflections. The abrupt loss of the dominant reflector beyond x = 40 m is caused by an increase in the x component of the bed slope, (a) shows a normal reflection originating from x = −20 m at roughly 94 m depth. The results used two sub-arrays with apertures of 60 m (in x) by 40 m (in y) centered along y = 175 m. The overlap was 40%. Sub-array sensor spacing was 2 m (in both x and y), giving 600 elements.
We systematically applied 3-D migration over the entire data grid, and normal reflections were picked from the results. We fit a bi-cubic spline surface through these data to produce a topographic surface of the glacier bed (Fig. 9). The maximum resolved depth is roughly 140 m at y values of >250 m. The surface contains significant dip components in both the x-z and y-z planes. There is a blank region at x > 75 m where no reflections were observed. The open circles plotted on the bed are reflection points picked from several y-z depth slices for sub-arrays positioned along x = 35 m. The lines projecting from the bed-reflection points connect to the surface sub-array positions that received that particular normal reflection. Thus, these lines can be thought of as representing exploding reflector rays.
Perspective drawing of the bed topography as determined by 3-D array results. The open circles on the surface of the bed are points picked from several 3-D y-z depth slices for a sequence of array runs along x = 35 m. The lines projecting from the bed-reflection points connect with the position of the sub-arrays that observed the reflections. Thus, these lines can be thought of as the specular reflection ray paths.
Comparisons of 3-D, 2-D and raw-data depth profiles
In radar reflection environments such as Figure 9, there will be substantial error in depths determined from either raw data or 2-D migration. To quantify these errors, a topographic surface was constructed from the 2-D migration-bed profiles. This was subtracted from the 3-D topographic surface determined from our methods. The resulting difference surface is shown in Figure 10. The region of largest disagreement between the 2-D and 3-D surfaces follows a diagonal trend, where the steepest slopes (45°) occur in both surfaces. The maximum difference is roughly 14 m. In all locations the depth profiles from the raw data and the 2-D migrations are shallower than depths determined using the 3-D array. The 3D profiles are consistent with classic “U”-shaped cross-valley profiles, whereas the 2-D and raw-data profiles indicate more gently dipping surfaces. In the regions of steepest dip (around y = 100 m), the disagreement between the 3-D depth estimates and the 2-D and raw-data estimates is 14 and 35 m, respectively. Using these differences in depth, and assuming the 3-D result to be correct (with a depth of 95 m), implies that the 2-D depth estimate has a 15% error and the raw-data estimate has a 36% error.
In addition to correcting the depth errors, the bottom topography generated by the 3-D method extends outside the data-coverage surface in the up-dip direction and terminates well inside the data-coverage area in the down-dip direction. Specifically (in Fig. 10), at x positions in the interval −20 m ≤ x ≤ 1 m, the 3-D array has localized the reflecting surface, whereas the 2-D migration methods cannot position reflections in this region. In the interval 70 m ≤ x ≤ 100 m and y > 140 m, the 2-D migrations position reflections directly below the survey lines, whereas the 3-D array results show no resolvable reflections. Given a deepening reflector as x increases, it is obvious that the reflections seen in the 2-D image are really located up-dip at lower x.
Difference surface constructed from topographic surfaces generated by 2-D and 3-D array results. The maximum disagreement is roughly 14 m. There are regions where the 3-D array results show no returns (x > 70 m), whereas the 2-D array results erroneously place reflections in this region. At positions x < 1 m (shaded region on left), the 3-D results accurately position reflections outside the data-coverage area.
Independent quantitative support for the 3-D array results is given by 12 MHz survey lines that cross the eastern end of the grid. Figure 11a gives the geometry for two of these lines (lines 15 and 35) in the coordinate system of the array grid. Figure 11b shows the first 600 m of the 2-D migrated depth sections for lines 35 and 15. Line 35’s depth section passes diagonally across the grid between A and A′. In this interval there is a steepening of the bed slope. Line 15’s migrated section crosses the far eastern corner of the array grid in the interval labeled B–B′. It too indicates an abrupt increase in bed dip under the grid. The 12 MHz reflector depths are consistent with depth estimates determined from the 3-D array processing. Although these 2-D, 12 MHz results are processing out-of-plane reflections, the effect is small relative to the wavelength. Thus, we are justified in using the general character of the 12 MHz results to support our finding that steep dip in the down-flow direction resulted in the loss of bed reflections within the array grid.
(a) Layout of the 12 MHz survey lines relative to the 50 MHz array grid. The lines are generally parallel to the direction of ice flow, (b) 12 MHz, 2-D migrations along lines 35 and 15. Line 35 crosses the array grid diagonally between A and A′. The relative increase in bed slope at 190 m distance is roughly the rnidpoint of the array grid. The line 15 depth section crosses under the far eastern end of the array grid between B and B′. In this interval there is also a sudden increase in bed slope. This supports the 3D array results which show no bottom reflections originating in these regions of the array grid. Out-of-plane reflections are present in both sections, particularly before and after the steep dip sections.
False-layering suppression and SNR improvement
Suppression of reflections from outside the array beam increases SNR and reduces interpretation ambiguity. Figure 12 shows a side-by-side comparison of depth sections from 3-D and 2-D migration. Both images are y-z planes located at x = 21 m. The center line of the 3-D array was located along x = 71 m. This places the image plane 30 m up-dip (relative to x = 51 m, the northernmost line used in the 3-D array). Both the 2-D and 3-D results used comparable sub-array parameters. In the 2-D migration the bed profile is discernible only by the intermittent, weak, but spatially coherent, amplitude variation that trends counter to the dominant noise field. This migration image also shows an apparent layer sequence, which could easily be interpreted as evidence of glacial till or debris-filled basal ice. These features were common in many of our 2-D migrations. However, none of our image slices from the 3-D array processing shows such layer structures. A possible explanation for the apparent layering is the superposition of two (or more) slightly out-of-plane bed reflections.
Comparison of 3-D and 2-D y-z image planes along x = 21 m. In the 2-D migration the bed reflection is barely discernible above the noise field and there is a false-layering effect caused by nearly coincident out-of-plane arrivals. No corresponding layering is seen in the 3 -D depth section. The 3-D sub-arrays were positioned along x = 71 m and had aperture dimensions of 40 m (in x) by 60 m (in y) with sensors spaced at 2 m intervals, giving 600 elements. The 2-D result has a y aperture of 60 m with a sensor spacing of 0.5 m, giving 120 elements. Both results used three sub-arrays overlapped by 17%.
Relative to the 2-D migration image, the 3-D array image (Fig. 12) shows a bed profile with a much higher degree of spatial coherence and a higher SNR. Comparing the mean noise power to the mean bed-reflection power in the 2-D and 3-D results gives an SNR of 0 and 12 dB, respectively. In the 3-D result, the bed-reflection wavelet out to y = 160 m is relatively undistorted, indicating good focusing. Beyond y = 160 m, the reflector becomes less coherent and the wavelets show minor distortion. Repositioning of the image plane by as little as 5 m would bring this portion of the bed into better focus. Even with the slight focus problem, the 3-D array image convincingly contradicts the suggestion of layering seen in the 2-D migration image.
Conclusions
Using a 50 MHz monostatic short-pulse radar system, we have demonstrated that full-wavefield, 3-D SAR processing is practical and adds important capabilities for characterizing cluttered temperate glaciers. Using many parallel survey lines oriented transverse to the ice flow, we produced depth images oriented both parallel and transverse to the direction of ice flow. Our analysis methods offer several important improvements: 3-D localization of all radar reflections, improved interpretation of bottom reflections, increased SNR against most types of noise, and increased depth of signal penetration. These improvements result from the development of a very narrow, steerable array beam pattern, and preservation of the short-duration wavelets characteristic of our electrically small broadband dipoles. Large cost savings and better imaging resolution can undoubtedly be obtained by automatic incorporation of differential GPS data directly into the digital radar data records.
Classic beamformer array processing broadly defines “noise” as any energy arriving from directions outside the focus point of the array. By this definition, even reflections from the glacier bed are considered “noise” (and are suppressed) if they do not originate at the focus point in an imaged plane. It is this mechanism that allows high-resolution localization of normally oriented facets of the bed. In the temperate-glacier context, 3-D array processing reduces clutter because, in a sense, the array beam “sees” only the scattering inclusions that lie in the image plane. All other forms of clutter, including reflected airwaves, englacial inclusions and out-of-image-plane bed reflections, are suppressed.
The 3-D array was able to resolve 45° bed dips. The deepest reflections observed were at depths of 140 m. We conclude that changes in topography resulted in loss of signal, not signal scattering from ice inclusions. Near the steeply sloping wall of the glacier, we found that interpretation of raw-data profiles underestimated depths by 36%, and that 2-D migrations underestimated depths by 16%. Furthermore, 3-D array processing improved bed-reflection SNR by >12 dB, clearly bringing the reflector wavelet above the background noise level. Possible interpretation pitfalls such as mislocation of bottom geometry and false layering are clarified.
Though we have no direct evidence (such as drilling) to confirm our 3-D results, there is a high degree of internal consistency, and plausible, simple, physical explanations readily support our findings. For example: (1) reflection picks from many different sub-array positions formed a smooth continuous reflector surface; (2) the bottom topography is consistent with a classic “U”-shaped cross-valley profile, which is not seen in the 2-D results; (3) in the down-flow direction the bed surface has a reasonable dip profile; (4) we can only localize bottom reflections outside our arrays in the up-slope/up-flow direction; (5) with few exceptions, we cannot localize bottom reflections directly beneath our arrays, particularly at the southern end of the array grid (this is explained easily with simple ray theory and a dipping reflection surface); and (6) our resolved bottom signals are coherent wavelets with high SNR that have properties that fall well within resolution limits obtained from beamforming theory.
Acknowledgements
Support for this research was provided by the U.S. National Science Foundation (NSF grants OPP-9531800 and OPP-9531579), by the U.S. Army Cold Regions Research and Engineering Laboratory (CRREL) project 4A161102AT24, work unit AT24-SS-E10, and the CRREL project 96-ILIR, titled “Signal Penetration Depth and Noise Suppression Improvement Using A Time Domain Radar Array”.
We wish to express many thanks to L. Davis and P. Annan of Sensors and Software Inc., Mississauga, Ontario, Canada, for use of equipment and field assistance during data collection at Camp Borden. The Borden data were used extensively in developing the methods applied in this paper. Further thanks go to G. Kuehn of Arete Training and Solutions, Hobart, Tasmania, Australia, for his professional support in providing glacier safety, base-camp management and technical services while collecting glacier radar data. Thanks go to the organizers of and participants in the 12th Annual Arctic Man Ski and Sno-Go Classic, 1997, for their consideration during our experiments.
Appendix: Discussion of 3-D Array Properties
The Kirchhoff integral in Equation (1) is evaluated over the surface containing sensor observation points. In 3-D array processing, Reference SchneiderSchneider (1978) notes that this procedure works by summing over a 3-D hyperbolic time–space surface in data space. In 2-D applications, the integral is evaluated along a line, and the summation is along a hyperbolic time–space curve. When focusing at a spatial coordinate that contains a scattering inclusion, coherent energy will sum in phase and the signal level will be elevated. At focus locations where no scattering originates, the data summation is incoherent, resulting in a (relatively) reduced signal level.
The image formed by our migration procedure has much in common with the process of aiming direct sum beams at a subsurface point. For the common source–receiver configuration, the received energy has reflected from a specular point where the ray path of the signal is normal to the reflection surface. This idea is embodied in the exploding-reflector concept (Reference ClaerboutClaerbout, 1985). The location, Γ, is where this normally projecting ray intersects the recording surface. In carrying out the migration calculation of Equation (1), the main contributions to an imaged point on the reflector come from receivers surrounding Γ. Diffraction theory explains that the phase is stationary at Γ and slowly varying in the region surrounding Γ, allowing coherent summation of wave energy.
The array resolution in the subsurface depends on the array size. A single receiver has a horizontal resolution distance that is defined by the Fresnel zone radius
, where z is the depth and λ is the in situ wavelength. Reference Welch, Pfeffer, Harper and HumphreyWelch and others (1998) point out that migration using data over the complete surface gives a horizontal resolution distance of λ/2, which is approximately 2 m at 35 MHz. If data over the entire bounding surface are not available, or not used, the resolution distance becomes larger. Migration resolution distance is related to the concept of beamwidth from an array of sensors. Figure 4a shows the normal frequency-domain beam responses (Reference Capon, Greenfield and LacossCapon and others, 1969) for linear arrays of lengths 20, 40 and 80 m. The sensors are spaced at 2 m intervals. The target is a point reflector 100 m directly below the center of the array. The 6 dB beamwidth decreases with array length, giving resolution distances of 9.2, 5.0 and 2.6 m, respectively. This is equivalent to 5.3°, 2.9° and 1.5° in angular resolution. Since the array sizes are of the order of the depth, the far-field approximation to the beam pattern may not be valid (Reference Stutzman and ThieleStutzman and Thiele, 1981), so a numerical calculation was done for the linear-array response. Figure 4b shows the result of applying Equation (1) to the three linear arrays with isotropic sensor elements. The data use a 35 MHz short-pulse wavelet. The migration beamwidth is very similar to the frequency-domain beamwidth shown in Figure 4a. Extending these results to migration with rectangular 60 m by 40 m arrays gives approximately 4 m by 6 m horizontal resolutions for reflectors buried at 100 m depth.
Alternative migration methods, such as frequency-wave-number (Reference StoltStolt, 1978; Reference Welch, Pfeffer, Harper and HumphreyWelch and others, 1998) and finite difference (Reference ClaerboutClaerbout, 1976), require a spatial data sampling interval less than the Nyquist interval for the largest wave-number of interest. This roughly corresponds to one-quarter of the wavelength of the GPR center frequency. In contrast, time-domain Kirchhoff migration tolerates a spatial sample interval larger than one wavelength (about 3.4 m in ice at 50 MHz). This comparatively large sample interval is vital to keeping data-collection and -storage costs reasonable. The low spatial sample densities are also less computationally demanding, making the full 3-D waveform migration practical.
