Compressive sensing for subsurface imaging using ground penetrating radar

 

Compressive sensing for subsurface imaging using ground penetrating radar

Ali C. Gurbuz^{a, , 1,} , James H. McClellan^b and Waymond R. Scott Jr.^b

^aDepartment of Electric and Electronics Engineering, TOBB University of Economics and Technology, Ankara 06560, Turkey

^bSchool of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA 30332-0250, USA

Received 27 March 2008; 

revised 17 February 2009; 

accepted 23 March 2009. 

Available online 1 April 2009.

Abstract

The theory of compressive sensing (CS) enables the reconstruction of sparse signals from a small set of non-adaptive linear measurements by solving a convex ℓ₁ minimization problem. This paper presents a novel data acquisition system for wideband synthetic aperture imaging based on CS by exploiting sparseness of point-like targets in the image space. Instead of measuring sensor returns by sampling at the Nyquist rate, linear projections of the returned signals with random vectors are used as measurements. Furthermore, random sampling along the synthetic aperture scan points can be incorporated into the data acquisition scheme. The required number of CS measurements can be an order of magnitude less than uniform sampling of the space–time data. For the application of underground imaging with ground penetrating radars (GPR), typical images contain only a few targets. Thus we show, using simulated and experimental GPR data, that sparser target space images are obtained which are also less cluttered when compared to standard imaging results.

Keywords: Compressive sensing; Synthetic aperture; Ground penetrating radar (GPR); ℓ₁ Minimization; Subsurface imaging; Sparsity; Source localization

Article Outline

1. Introduction

2. Theory: CS for GPR synthetic aperture imaging

2.1. Creating a dictionary for GPR data

2.2. Compressive sensing data acquisition

2.3. GPR imaging with compressive sensing

3. Selection of algorithm parameters

4. Discussion

5. Results

5.1. Performance in noise

5.2. Effect of the measurement matrix Φ

5.3. Random spatial sampling

5.4. Experimental results

6. Conclusions

References

1. Introduction

Imaging with a synthetic aperture typically requires the collection of a large amount of space–time data. The synthetic aperture is formed by scanning a sensor over the region of interest and recording the time signal returns for many spatial positions. When the sensor is a ground penetrating radar (GPR) the measurements are made by transmitting short electromagnetic pulses into the ground and recording the reflections [1]. The total subsurface response, which is a superposition of the responses from all reflectors within the medium, can be inverted using various imaging algorithms, e.g., time-domain standard backprojection (SBP) [2] X. Feng and M. Sato, Pre-stack migration applied to GPR for landmine detection, Inverse Problems 20 (2004), pp. 99–115. [2]. The SBP algorithm requires fine spatial sampling and Nyquist-rate time sampling of the received signals in order to perform matched filtering with the impulse response of the data acquisition process to form an image. Also, SPB does not use any prior knowledge about the image. In this paper, we approach imaging as an inverse problem based on the same model as used in the SBP algorithm, but we compute the image within a regularization framework, so we can take into account prior information.

We formulate the GPR subsurface imaging problem as a dictionary selection problem, where the dictionary entries are produced by discretizing the target space, and synthesizing the GPR model data for each discrete spatial position [3]. Recently, there has been considerable interest in representing signals as a superposition of elements from an overcomplete dictionary, and many methods have been developed to extract an “optimal” representation in terms of the given dictionary elements; this kind of search is called basis pursuit [4]. The assumptions made here are that the targets are point reflectors at discrete spatial positions, the targets do not interact so superposition is valid, and the wave propagation obeys ray theory. The reason we use a point target reflectivity model is simplicity in calculating the received data from the model. The point-target assumption is not crucial; other models can be used as long as the received data can be calculated for the assumed target model.

Generally, potential targets like buried landmines and metallic objects cover a small part of the total subsurface to be imaged. Thus, our basic a priori assumption about the GPR image is sparsity of targets. In other words, the number of targets is much less than the total number of discrete spatial positions, e.g., a image usually contains less than five targets, and rarely more than 10. Sparsity is not limited to the point-target model; other models can be used as long as the targets are sparse or compressible in some transform domain, e.g., smooth targets are sparse in the Fourier or wavelet domain.

Recent results in compressive sensing (CS) state that it is possible to reconstruct a K-sparse signal x=Ψs of length N from measurements [5], [6] and [7]. CS takes non-traditional linear measurements, y=Φx, in the form of randomized projections. The signal vector, x, which has a sparse representation in a transform domain Ψ, can be reconstructed exactly (with high probability) from M=C(μ²(Φ,Ψ)logN)K compressive measurements by solving a convex optimization problem of the following form:

(1) which can be solved efficiently with linear programming. The constant C is usually small and μ(Φ,Ψ) is the mutual coherence [5] between Φ and Ψ.

Instead of taking traditional time-samples, a very small set of informative measurements in the form of random linear projections will still allow us to obtain the subsurface image. This will reduce the load in the data acquisition part of the problem, but put more work on the signal processing. Note that we do not take the random projections of the target space; instead, we take random projections of the received signals at each scan point. Then, we formulate the CS reconstruction problem using our model(s) for the received signals as a function of target space positions.

The next section describes the CS formulation of the subsurface imaging problem and discusses its advantages and how possible problems are handled. Section 3 explains the choice of algorithm parameters and Section 4 discusses handling non-ideal conditions. Section 5 presents results for simulated and experimental GPR data, along with a performance analysis of the proposed algorithm for different parameters.

2. Theory: CS for GPR synthetic aperture imaging

Ground penetrating radar uses electromagnetic (EM) radar pulses to image the subsurface. When the transmitted EM wave reflects from a buried object or a boundary with different dielectric constants, the receiving antenna records the reflected return signal [8] and [9]. Although the response of possible targets like landmines could be more complex, we use the common point-target model [2], [8] and [9] due to its simplicity. It is important to note that the point target model is not crucial to the method developed in this paper; a more sophisticated target model could be used. We use the model to write the received signal reflected from a point target as a time delayed and scaled version of the transmitted signal s(t),

(2) where τ_i(p) is the total round-trip delay between the antenna at the ith scan point and the target at p, σ_p is the reflection coefficient of the target and A_i,p is a scaling factor used to account for any attenuation and spreading losses.

For imaging, the GPR sensor is moved in space to collect reflected time data at different scan positions. The collection of these scan points constitute a synthetic aperture that forms the space–time data. Thus it is essential to know the space–time response from a point target. We typically use a bistatic GPR with separate transmitting and receiving antennas separated by a fixed transmitter–receiver (T–R) distance d_tr, as shown in Fig. 1(a). Calculating the travel-time τ_i requires knowledge of the path the wave travels from the transmitter antenna to the target and then to the receiver, as well as the wave velocities in the different media.

Full-size image (24K)

Fig. 1. (a) Bistatic GPR measurement setup. T denotes the transmit antenna, R the receiving antenna, and d_tr the fixed offset between T and R. The midpoint, x_i, is chosen as the location of the sensor at ith scan point. (b) Space–time impulse response of the GPR data acquisition process for a single point target scanned from -70 to by a bistatic GPR with and .

View Within Article

At the boundary between two different media with different dielectric properties, e.g., air and soil, the propagation direction changes according to Snell's law. Taking wave refraction into account is especially critical for near-field imaging problems, such as mine detection, since the distance between the antennas and the target is relatively small. Calculation of the exact refraction point requires the solution of a 4th degree polynomial. Several approximations to compute the refraction point are available [10]. After finding the refraction points, the distances d_1:4,i of the four segments of the path shown in Fig. 1(a) can be found and the τ_i times can be calculated as

(3) To collect data the T–R antenna pair is moved along the scan direction.² Fig. 1(b) shows the space–time data collected over a single point target at when scanning along the x-axis from -70 to at a fixed y slice. The scan location is referenced to the middle of the T–R pair. The response from a single point target follows a spatially variant curve in the space–time domain as shown in Fig. 1(b). When using data collected over this region of interest, the processing is called synthetic aperture imaging.

A model which is a simple linear superposition of targets is obtained if we ignore the interaction between the targets:

(4) where d(u_x,u_y,t) is the measured GPR space–time data at scan point (u_x,u_y), τ(u_x,u_y;x,y,z) is the total travel time from the transmitting antenna to the target space point (x,y,z) and on to the receiver, σ(x,y,z) is the reflection coefficient at the target space point (x,y,z), and A_L(u_x,u_y,x,y,z) accounts for spreading and losses. In (4) (x,y,z) constitutes the target space π_T which lies in the product space [x_i,x_f]×[y_i,y_f]×[z_i,z_f]. Here (x_i,y_i,z_i) and (x_f,y_f,z_f) denote the initial and final positions of the target space along each axis.³ The goal is to generate the reflectivity profile or in other words an image of the target space using the measurements d(u_x,u_y,t).

Standard imaging algorithms do matched filtering⁴ with the measured data with the impulse response of the data acquisition process (Fig. 1(b)) for each spatial location. Time-domain standard backprojection [2], [10] and [11] can be expressed as

(5) where f(x,y,z) is the created subsurface image, w(u_x,u_y) is the weighting function on scan positions typically used to reduce sidelobes in images, and the shifted is the impulse response.

Normally we work with the discrete versions of (4) and (5). If the measurements and the image are concatenated into vectors, a sampled f is related to the observed discrete (noisy) space–time domain data d through a matrix. This matrix embodies the assumed point-target model, and it can also be used to generate a dictionary of GPR responses.

2.1. Creating a dictionary for GPR data

A discrete inverse operator can be created by discretizing the spatial domain target space π_T and synthesizing the GPR model data for each discrete spatial position. Discretization generates a finite set of target points , where N determines the resolution and each π_j is a 3D vector [x_j,y_j,z_j]. Using (2) and (3) the signal at the GPR receiver can be calculated for a given element of using σ_p=1 in (2). For the ith scan point, the jth column of Ψ_i is the received signal that corresponds to a target at π_j as shown in Fig. 2. The nth index of the jth column can be written as

(6)