4010 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO.

7, JULY 2015

Efficient Algorithm Design for GPR

Imaging of Landmines
Kyle R. Krueger, James H. McClellan, Life Fellow, IEEE, and Waymond R. Scott Jr., Fellow, IEEE

Abstract—Ground-penetrating radar (GPR) is used to image is deals purely with the inversion algorithms, only the SFGPR
and detect subterranean objects, for example, in landmine de- will be examined and used in the laboratory experiments.
tection. Although full 3-D inversion of GPR measurements is Many GPR applications require high spatial resolution;
possible for simple algorithms such as backprojection, it is im-
practical when using more advanced algorithms that involve hence, a large antenna aperture would be required. When trying
1 -minimization. Many of the algorithms used for GPR imaging to locate very shallow buried targets, the desired resolution is
involve the storage, or online generation, of a huge dictionary on the order of centimeters. The GPR is usually a movable
matrix created from discretizing a high-dimensional nonlinear transmit–receive platform that is scanned over a region of
model. This parametric model includes all the target features that interest, e.g., some GPRs are handheld wands. Therefore, it is
need to be extracted, including 3-D location, object orientation,
and target type. As more parameters are added to the model, natural to use a synthetic-aperture approach to acquire coherent
the dimensionality increases. If uniform sampling is done over data over a large aperture [5]–[7]. Synthetic aperture radar
high-dimensional parameter space, the size of the dictionary and (SAR) has been shown to be effective in GPR applications
the complexity of the inversion algorithms rapidly grow, exceeding [2], [8]–[12]. A disadvantage with SAR is that it requires the
the capability of real-time processors. This paper shows that storage of the responses at each aperture position, increasing
strategic structuring of the dictionary, which takes advantage of
translational invariance in the model, can reduce the dictionary the memory requirements of the inversion system. A common
storage by several orders of magnitude and exploit the fast Fourier assumption in synthetic aperture is translational invariance of
transform for fast computation of previously highly impractical, the travel times from reflectors at a common depth, which
bordering on impossible, 3-D GPR imaging problems. facilitates coherent addition of the responses; this invariance
Index Terms—Compressive sensing (CS), fast Fourier trans- property is true in a layered-earth velocity model and can be
form (FFT), ground-penetrating radar (GPR). found in [8], [9], [13]–[16], as well as in [10]–[12].
The needs for extremely high resolution and data-acquisition
efficiency have led to investigations of many different inversion
I. I NTRODUCTION
techniques when dealing with the GPR problem. A few of these

G ROUND-penetrating radar (GPR) has been shown to be

an effective tool for imaging a large variety of sub-
terranean targets, as they are sensitive to changes in electri-
algorithms include backprojection (BP), orthogonal matching
pursuit (OMP), and 1 -minimization (1 M). All of these algo-
rithms fall into the dictionary-matching category. Dictionary
cal permittivity and conductivity of the subsurface [1]–[3]. matching is a technique where a synthetic model is built to
Development of GPR systems began in earnest in the early approximate the system response at different parameters and
1970s for applications ranging from detecting buried landmines then matched with the received measurements to determine the
to locating utilities [4]. Two common transmitter modes for intended environmental parameters. Each algorithm requires a
data acquisition are time-pulse GPR (TPGPR) and stepped- model dictionary, i.e., Ψ , which must be applied once, in the
frequency GPR (SFGPR). While these signals are Fourier trans- case of BP, or many times (on the order of 10 to 1000) for OMP
form pairs, the actual acquisition process of each method has and 1 M. Often, an explicit dictionary matrix, i.e., Ψ , is built by
its advantages and disadvantages. SFGPR systems are much enumerating the variables in the expected target response.
more robust to narrowband noise, but the receiver gain cannot The dictionary approach can become quite computationally
be adjusted for deep targets. Both the TPGPR and the SFGPR intensive since a 3-D imaging problem would require enumera-
systems are hindered by long data-acquisition times for many tion of a 2-D scan grid, enumeration of a 3-D simulated image
applications, but SFGPR is generally slower. Since this work space, and enumeration of the number of stepped frequencies,
or time samples. The matrix Ψ scales as O(N 6 ), if all variables
are generalized to equal discretizations of size N . If a reason-
Manuscript received November 1, 2013; revised April 15, 2014 and ably sized problem has N = 100, then 1012 complex elements
November 9, 2014; accepted November 25, 2014. This work was supported would have to be stored in computer memory. Each element
in part by the U.S. Army RDECOM CERDEC Night Vision and Elec- is going to be a complex double of size 64 bits, creating a
tronic Sensors Directorate, Science and Technology Division, Countermine
Branch, and in part by the U.S. Army Research Office under Grant W911NF- storage requirement of 8 TB. Since the application for these
11-1-0153. types of problems might require small mobile devices, 8 TB to
The authors are with the Department of Electrical and Computer Engi- store a part of the data used for the algorithm is unacceptable.
neering, Georgia Institute of Technology, Atlanta, GA 30332 USA (e-mail:
krkrueg@gmail.com). However, this paper will show that modifications to the struc-
Digital Object Identifier 10.1109/TGRS.2015.2388786 ture and application of Ψ with a functional representation can
0196-2892 © 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
KRUEGER et al.: EFFICIENT ALGORITHM DESIGN FOR GPR IMAGING OF LANDMINES 4011

Fig. 1. Detection flow for a GPR system with model-based inversion. Fig. 2. Uniform 2-D scan grid for GPR data acquisition.
be used to dramatically reduce the constraints in creating and
applying Ψ in a way that is increasingly more practical for the
GPR problem.
The key idea is that matrix Ψ can be structured by taking
advantage of the fact that Ψ can have a shift-invariance prop-
erty that creates a block-Toeplitz structure in Ψ . The block-
Toeplitz structure is associated with the dimensions where the
translational-invariance property is present. Taking advantage
of the translationally invariant structure of Ψ allows for the
matrix to be implicitly stored, and applied using a function,
i.e., gA , which reduces the storage requirements from O(N 6 ) to
Fig. 3. EM path through multiple media.
O(N 4 ). The computational expense to apply gA is much less,
O(N 4 ) instead of O(N 6 ) since gA will incorporate the fast A. Response Model
Fourier transform (FFT). If the example given in the previous The data-acquisition technique for a GPR includes the cre-
paragraph is reexamined using these data reduction techniques, ation of a synthetic aperture by moving a sensor to the positions
then the previous example that required 1012 elements would ls = (xs , ys , zs ). For a 3-D image, the scan positions, i.e., ls ,
only require 108 elements. This reduces the computer storage would be indexed over a grid in the 2-D plane as shown in
from approximately 8 TB to a much more manageable 800 MB Fig. 2, although almost impossible in practice, it is ideal for
before the use of a compressive sensing (CS) technique, which zs to be held constant.
could be used to further reduce the required storage and amount BP has been shown to be quite effective for imaging subsur-
of acquisition samples. face structures using a very simple scalar model for the system.
The remainder of this paper will be split into four sections. The imaging is effective even with the simplifications for the
Section II will describe the setup of the function gA with complete model, i.e.,
respect to the GPR landmine detection problem. Section III
will show simulated performance statistics for a variety of Aρ
R(ω, ls ; lt ) = e−jωτ (ls ,lt ,c,v) (1)
different inversion algorithms that can take advantage of the S(ls , lt , c, v)
new representation. Sections IV and V will provide results
such as the assumptions of a point scatter, a homogeneous
when processing measurements from a laboratory experiment
soil, and no interactions between the targets [2], [8]–[15]. The
and brief conclusions, respectively.
model in (1) describes a setup where a series of stepped-
frequency signals of strength A, at frequency ω, are sent from
II. M ODEL S ETUP AND I MPLEMENTATION
a transmitter located at ls and then reflected off a target at lo-
The GPR problem can be described as a parameter-detection cation lt = (xt , yt , zt ) with reflection coefficient ρ. The factor
problem that relies on model-based inversion. The data flow S(ls , lt , c, v) takes into account the electromagnetic spreading
and processing blocks for this detection system are shown in and the transmission through the boundary. The time delay
Fig. 1. There are two sources of sampling: acquisition sampling function τ is computed using the plane-wave approximation
and model discretization. Acquisition sampling is where a and an approximation of Snell’s law to calculate the wave path
sensor measures the environment at many locations. Model through air with velocity c and at the boundary of a medium
discretization is where a dictionary, i.e., Ψ , is constructed using with velocity v [16]. An example of the wave path from a
a sampling of a continuous model created from the knowledge bistatic GPR traveling into the ground and reflecting off a target
of the data-acquisition system and the parameters that need can be seen in Fig. 3, where τ = τ1 + τ2 + τ3 + τ4 .
to be extracted. The measurements collected with the data- Examples of a simulated time-domain and frequency-domain
acquisition system and the dictionary Ψ are combined during measurement can be seen in Fig. 4. There will also be a
the inversion step to obtain a 3-D image that estimates the reflection from the air–ground interface, and eliminating this
environment. ground response is an active research topic in its own right.
4012 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO. 7, JULY 2015

where the amplitude s(lt ) is a simple amplitude used to ap-

proximate the contribution of Aρ/S(ls , lt , c, v) from (1) and
(2), and s is a sparse vector that is only nonzero at the target lo-
cations. The indexing of vector r must follow the enumeration
of the triple (ω, xs , ys ) used for the measurement vectors ψ(lt ).
The indexing of s must follow that of the triple (xt , yt , zt ) used
for the target locations.
The inversion process is done to recover s from the measure-
ments r. However, the size of Ψ and computational complexity
Fig. 4. Simulated measurements for 2-D in the (a) time domain and the of the inversion algorithms make it problematic for real-world
(b) frequency domain showing the phase. applications. To address this issue, some structural changes
In particular, there is some work associated with ground re- are presented that can be made to simplify the way that the
moval for GPR geared toward a CS application introduced by dictionary Ψ is created, stored, and applied in the different
Tuncer et al.; however, ground removal will not be considered inversions.
in this paper [17]. That is, the simulations shown here use
a single medium with no velocity changes. The subsurface
laboratory experiments in Section IV have a fairly uniform B. Shift Invariance Property
ground response; hence, it is possible to subtract out the ground A simplification of Ψ is possible because the GPR acqui-
bounce from the measurements since the height of the sensor is sition system can have the extremely powerful property of
known in the controlled laboratory experiment. spatial shift invariance [18]–[20]. The key idea is that the model
The model dictionary, i.e., Ψ , that can be used with any response (1) at the horizontal aperture of sensors will shift
number of different imaging techniques can be created by in tandem with horizontal shifts in the target positions at a
enumerating (1) for all possible parameter discretizations. The fixed depth. This is true because the distance function d(ls , lt )
first step is to determine which parameters are associated with shown in Fig. 3 does not change with equal horizontal shifts of
the measurements and which will be extracted in finding the the sensor and target. There are a few important assumptions
targets. The highest priority variable in a landmine detection that are being made in order to make the shift invariance
system is target location, i.e., lt . For the remainder of this sec- true. The first assumption is that the soil model is reasonably
tion, only the SFGPR case will be studied; hence, the frequency homogeneous in the lateral plane within the ls grid used for a
ω and sensor locations ls are the measurement parameters. The particular image. The second assumption is that the height of
discussion could easily be switched to TPGPR. The SFGPR the sensor relative to the ground remains constant. The second
model for a single target can be rewritten as assumption is difficult to maintain in practice, particularly with
Aρ a handheld device. However, there are techniques that try to
R(ω, ls ; lt ) = e−jωτ (d(ls ,lt ),c,v) level out the ground response in an attempt to approximate
S(ls , lt , c, v)
Aρ the case when the radar is at a constant height [21]. This will
= ψ(ω, ls ; lt ) (2)
S(ls , lt ) slightly warp the data but was shown to be a rather effective
way of both removing ground bounce and laterally aligning
where d(ls , lt ) is a 3-D distance function. The measurements the measurements. The computation takes place with discrete
are made at a finite set of frequencies ω and a finite number grids for the positions ls and lt ; hence, the grids must also
of sensor locations ls = (xs , ys ). The velocities c and v are support the shift invariance. While it is not necessary in the
removed to make the notation more compact. general case for the sensor locations ls to be uniformly spaced,
One column vector of the dictionary matrix Ψ is the vector in this paper, it will be assumed that ls is uniformly spaced.
ψ(lt ) whose entries are all the measurements created by evalu- If ls is not physically uniformly spaced, it can be interpolated
ating the model R(ω, ls ; lt ) for a fixed lt while enumerating all onto a uniformly spaced grid. To show how the shift-invariance
possible triples of the three measurement space parameters, i.e., property simplifies the computation for the collection of SAR
ω, xs , and ys . The final dictionary is created by concatenating measurements, a 2-D example with the target in the air will
the target location response vectors as the column space, i.e., be used. Because the target and sensors are in the air, there
 
Ψ = ψ(lt 1 ) ψ(lt 2 ) · · · ψ(lt Nlt ) (3) is no change in medium for the wave; hence, the wave veloc-
ities will be dropped from the equation to keep the notation
where the Nlt values of lt are obtained by enumerating all simpler.
possible triples of the target location parameters, i.e., xt , yt , and First, rewrite the response vector from (4) as a sum of
zt . If the number of values for each parameter is N , then Ψ is an products
N 3 × N 3 matrix. For example, with N = 100, the dictionary
matrix is 106 × 106 with one trillion (1012 ) entries. 
Using the dictionary matrix in (3), the response vector can be r(ω, xs ) = s(xt , zt )e−jωτ (xs ,xt ,zt ) . (5)
zt xt
expressed as

r= s(lt )ψ(lt ) = Ψ s (4) Next, we discretize the x-dimension as follows: mΔx for m =
lt 1, 2, . . . , Nxs for xs and (h + α)Δx for h = 1, 2, . . . , Nxt for
KRUEGER et al.: EFFICIENT ALGORITHM DESIGN FOR GPR IMAGING OF LANDMINES 4013

xt , where α is a constant value 0  α < 1. The x discretization The convolution representation shown in (10), i.e., (shown
creates a new representation for the response vector, i.e., at the bottom of the page) uses the same dictionary as was
r(ω, mΔx) created in (2), but the discretization for xt has been replaced
by the constant α. Changing the discretization from xt to α
N xt
is significant because it drops the storage requirements from
= s ((h + α)Δx, zt ) e−jωτ (mΔx,hΔx,α,zt ) . (6)
zt h=1
Nxt to one. However, if there is a desire to upsample the image
locations, e.g., to get upsampling by a factor of 2, we could
If the time delay τ is examined for a monostatic system mea- use α1 = 0 and α2 = 0.5, then the process would need to be
suring in a single medium, we obtain repeated for each αi . With these steps outlined for the forward
τ (mΔx, hΔx, α, zt ) operator, the adjoint is fairly trivial.
 
1 Transitioning back to matrix notation, shift invariance leads
= (mΔx − (h + α)Δx)2 + zt2 (7)
c to a Toeplitz, or block-Toeplitz, structure in Ψ . For simplicity,
 
1 consider an example with α1 , α2 , and Nα = 2, which could
= ((m − h) − α)2 (Δx)2 + zt2 . be trivially expanded. The columns correspond to the sensor
c
positions, xs , the rows correspond to the simulated target
Thus, the time delay depends on the index difference (m − h),
locations, xt , and the entries in the D H matrix correspond to
and it can be shown that the inner sum in (6) is a discrete
the distance between xs and xt . The distance between xs and
convolution. The time delay function τ keeps horizontal shift
xt is going to be |m − h − αi |. A reasonable assumption made
invariance even when the system is not monostatic and when
here is that Nxt = Nxs . When they are not equal, the Toeplitz
the medium velocity is allowed to change with z, but not with
structure would just not be square. The difference matrix is
x or y. The horizontal shift invariance can be visualized with
shown in (10). The matrix in (10) is not Toeplitz, but with a
Fig. 3. If the target, T1 , and R1 are shifted by an equal value, d1
slight reorganization of the rows, it can become block-Toeplitz
and d2 will shift by the same amount; hence, τ will be the same
with Nα blocks, as shown in (11). The difference matrix, i.e.,
as before the shift. Therefore, with a fixed Δx, the exponential
D, is built directly into the representation matrix
in (6) is a function of the index difference (m − h), i.e.,
e−jωτ (mΔx,hΔx,zt ) = e−jωτ ((m−h)Δx,α,zt ) . (8) Ψ (ω1 , D, zt ) = ejωτ (D,tz1 ) (12)
Now, we can rewrite the inner sum of (6) as a convolution giving Ψ a block structure with Nα Nω Nzt blocks of size
because we have index shift invariance Nxs ×Nxt . It is well known that a block-Toeplitz matrix can
r(ω, mΔx) be stored with a single vector for each block.1 Toeplitz matrices
⎛ ⎞
have been shown to be an effective way to reduce computational
 
N xt

= ⎝ s ((h + α)Δx, zt ) e−jωτ ((m−h)Δx,α,zt ) ⎠ complexity in random sampling matrices for CS [22], [23]. A
zt h=1 graphical example of the structural changes that would need to
   be made to the storage and application of the Ψ matrix used for
 convolution w.r.t. m
2-D imaging in the time domain can be seen in Fig. 5. Fig. 5(a)
= s ((m + α)Δx, zt ) ∗m e−jωτ ((m)Δx,α,zt ) (9) shows the traditional dictionary, where every simulated target

zt
= s ((m + α)Δx, zt ) ∗m ψ(ω, mΔx, α, zt ). 1 For the examples given in this paper, it will be assumed that the sensor
zt spacing and the simulated target spacing are identical, i.e., Nα = 1, α1 = 0.

⎡ ⎤
α1 Δx − α1 ··· (Nxs − 1) Δx − α1 )
⎢ α2 Δx − α2 ··· (Nxs − 1) Δx − α2 ) ⎥
⎢ ⎥
⎢ Δx + α1 α1 ··· (Nxs − 2) Δx − α1 ) ⎥
⎢ ⎥
⎢ Δx + α2 α2 ··· (Nxs − 2) Δx − α2 ) ⎥
DH =⎢ ⎥ (10)
⎢ .. .. .. .. ⎥
⎢ . . . . ⎥
⎢ ⎥
⎣ (Nxt − 1) Δx + α1 (Nxt − 2) Δx + α1 ··· α1 ⎦
(N − 1) Δx + α2 (Nxt − 2) Δx + α2 ··· α2
⎡ xt ⎤
α1 Δx − α1 · · · (Nxs − 1) Δx − α1 )
⎢ Δx + α α1 · · · (Nxs − 2) Δx − α1 ) ⎥
⎢ 1 ⎥
⎢ .. .. .. .. ⎥
⎢ . . . . ⎥
⎢ ⎥
⎢ (Nxt − 1) Δx + α1 (Nxt − 2) Δx + α1 ··· α1 ⎥
DH =⎢
⎢
⎥ (11)
⎢ α2 Δx − α2 · · · (Nxs − 1) Δx − α2 ) ⎥
⎥
⎢ Δx + α α2 · · · (Nxs − 2) Δx − α2 ) ⎥
⎢ 2 ⎥
⎢ .. .. .. .. ⎥
⎣ . . . . ⎦
(Nxt − 1) Δx + α2 (Nxt − 2) Δx + α2 ··· α2
4014 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO. 7, JULY 2015

Fig. 5. Dictionary implementation. (a) Explicit enumeration with matrix multiplication. (b) Exploiting shift invariance by using a convolution operator.

vector multiplication to apply Ψ . The matrix-vector multipli-

cation has to be replaced with a specifically designed func-
tion to perform the equivalent of the forward, i.e., gA , and
H
adjoint, i.e., gA , operators of a reduced matrix Ψ α to the
sparse vector s [19]. It is important to note that Ψ α is built
by enumerating (2), where lt is only enumerated for z and
Fig. 6. Storage requirements for the GPR system. not x or y.
The first step is to make sure that the discretization of ls and
position lt is enumerated into the columns of Ψ , and the dic- lt provide the translational invariance discussed in Section II-B.
tionary is applied using standard matrix-vector multiplication. Next, identify how the convolution operation in Fig. 5(b) is
Fig. 5(b) shows a reduced size Ψ , where Nα = 1, that does not going to be performed to take advantage of the translational
enumerate in the horizontal dimensions and is applied using a invariance. The FFT can be used to perform circular convolu-
convolution operator along the horizontal position, instead of a tion efficiently, and with the use of a zero-padding operator,
simple matrix-vector multiplication. Z, linear convolution. The FFT can only be performed in the
The particular Ψ structure shown in Fig. 5(b) has the horizontal dimensions if the time and frequency samples are
added bonus that the computational operations required are the same for each ls . Since there is translational invariance in
O(N log(N )), by using the FFT, instead of O(N 2 ) for each both dimensions of ls , zero padding must take place in both
dimension where the shift invariance can be exploited. A tradi- dimensions. The simplest thing to do is to add Nxs /2 samples
tional (explicit) Ψ used to image three dimensions can be stored to the beginning and the end of the xs dimension and do the
and applied in O(N 6 ), assuming all measurements and parame- same for y. The zero pad will allow for shifts to take place
ters are equally discretized. On the other hand, when Ψ has the within the desired range of ls without wraparound effects from
Toeplitz structure in both horizontal lt dimensions equivalent the circular convolution. A slightly more efficient way would
to ls , taking advantage of the FFT would reduce the storage be to zero-pad with a number based on the maximal beamwidth
to O(N 4 ) and the computation to O(N 4 ). A flowchart for the of the antenna in space. Then, the FFT needs to be taken in
special properties of the GPR and the resultant dimensionality both x and y of ls . The same operations must be performed to
reductions can be seen in Fig. 6. The addition of using a CS the corresponding x and y of lt in the sparse vector s. These
inversion would allow for a further reduction in the frequency are then combined through a series of index multiplications
domain by using M  N frequencies. specifically shown in Algorithm 1. The summation over zt
should be interpreted as a for loop over depth that adds the
C. Implementation Specifics for Structure Change
results of FFT (horizontal) convolution done for each depth.
Now that the structure within the dictionary has been iden- The final step is to then compress the measurements with Φ if
tified, the inversion algorithms no longer use a simple matrix- a compression algorithm is being used. In the case where no
KRUEGER et al.: EFFICIENT ALGORITHM DESIGN FOR GPR IMAGING OF LANDMINES 4015

compression is required, Φ can be removed or equivalently set for faster computation when using MATLAB. For example,
H
to the identity matrix. running BP using gA on the 3-D simulation from Section III-C
takes 5.28 s if the multiplications and additions are done in a
loop and 0.13 by using vectorized operations to perform the
same task.
The algorithms described in Algorithm 1 and Algorithm 2
H
describe gA and gA if a general Φ were to be used in a
compression algorithm. However, to use a general Φ requires
complete sampling during data acquisition and eliminates the
usefulness of compressive algorithms in terms of data acquisi-
tion. Making sure that Φ is a random sampling matrix where
the random frequencies, or random time samples, at each ls
are equal, would allow for Φ to be applied to Ψ α before the
zero pad and FFT take place. Applying Φ before the zero pad
and FFT allows for the compressed sampling to be performed
during data acquisition. In other words, if ΦΨ has the same
shift-invariant property as Ψ , then Φ can be applied before the
zero pad and FFT operations.

III. S IMULATED P ERFORMANCE

H
The forward gA and adjoint operators gA can be used in
different ways to perform inversion of the GPR measurements.
In some cases, the operators are applied only once, or a few
Now that a function to perform the forward operation has times, in other algorithms thousands of times. This section
been created, a similar function for the adjoint operation can be is split into four subsections. Section III-A will introduce the
described. The adjoint matrix Ψ H is the conjugate transpose; three different inversion algorithms that will be examined
hence, it possesses Toeplitz subblocks, and FFT convolution in the remaining subsections: BP, OMP, and compressed
can be done. In addition, the adjoint can be created by working 1 -minimization (C1 M). All of the simulated examples in this
the forward algorithm backward. The specifics for how to section use point targets with 379 equally spaced frequencies
H
perform the adjoint operator gA can be found in Algorithm 2. from 0.50 to 8.06 GHz. The measurements are created by using
the forward model from (1) and adding noise. There is assumed
to be no interaction between the target responses. Section III-B
directly compares the computation times of the the three
algorithms with and without exploiting translational invariance.
Only 2-D examples can be run in this case because the 12-GB
computer used to run the simulations was unable to run the
3-D inversions using explicit matrices. Section III-C shows the
different algorithms successfully being inverted in 3-D while
comparing their capabilities in dealing with compressed mea-
surements. The 3-D compressed-algorithm simulation helps
show the advantages of using an algorithm that looks for spar-
sity over those, such as BP, that do not. The last simulation in
Section III-D is set up to compare solving a full 3-D inversion
exploiting translational invariance and creating a 3-D image
out of 2-D slices using the explicit matrix method. The 2-D
slice inversion is unable to take advantage of the full synthetic
aperture and will introduce imaging artifacts that are not
present in the full 3-D inversion.

A. Inversions
The most common inversion method for GPR is BP. In the
The description of the algorithms above can be modified to notation of (4), the standard BP image is created simply by
get a more time-efficient implementation. Many of the steps can multiplying the measurements, i.e., r, by the adjoint Ψ H , i.e.,
Ψα and taking
be performed offline, for instance, calculating Ψ̃
the transform of b. Moreover, the for loops are easy to vectorize sbp = Ψ H r = (Ψ
ΨH Ψ )s = gA
H
(r). (13)
4016 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO. 7, JULY 2015

roughly equivalent to the sparsity of the image. When compared

with BP, the computation time is greater. The computation
time of OMP is approximately twice the number of iterations
multiplied by the time it takes to run the BP algorithm. On
the other hand, the sparse solutions should be much easier to
interpret, particularly when compressive sampling is employed
to use fewer (random) frequencies.
An increasingly popular technique is CS, which has been
extensively studied in the imaging community [27]–[29]. The
ability of CS to take a sparse representation of a signal and
Fig. 7. Standard BP of a simulated target with lt resolution of (a) 2 cm and project it onto a much lower-dimensional measurement space
(b) 1 cm.
can decrease the data-acquisition requirements dramatically.
CS also produces sparse images, which some have claimed pro-
The output vector sbp must be reshaped to display the 2-D vide an increase to the detection resolution above that of other
images seen in Fig. 7. Fig. 7(a) shows a single target image with commonly used detection algorithms for GPR [12]. A common
2-cm discretization in lt , and Fig. 7(b) shows the same target inversion used in CS is the compressed 1 -minimization tech-
imaged with 1-cm discretization. The 1-cm target discretization nique (C1 M).
case is an example of when there are two α values, α1 = 0, The images created using BP and C1 M with the same
α2 = 0.5, and Δx = 2 cm. Both examples have the same 2-cm two randomly selected frequencies at each scan position and
discretization for the sensor positions ls . an SNR of 10 dB are shown in Fig. 9. The advantage of
Although there is only one target, it is easy to see that there C1 M is fairly easy to see because the locations of the three
are many nonzero image pixels containing energy, not just one targets are obvious in Fig. 9(b). On the other hand, the BP
at the true target location. This is because Ψ is not an orthogonal image in Fig. 9(a) is very difficult to interpret. Reducing the
matrix, and the correlation between two nonidentical columns number of measurements taken has degraded the effectiveness
is not always zero. The coherence of the dictionary, defined by of BP dramatically while providing no degradation to the C1 M
the maximum correlation over column pairs, will increase for solution. The algorithm used to solve the 1 -minimization in
a finer grid where the spacing of lt decreases. The increase in C1 M is SPGL1 basis pursuit denoising (BPDN), which works
coherence can be graphically observed in Fig. 7 where the finer well because of its ability to handle model error and noise, its
grid, Fig. 7(b), has higher valued off-target pixels close to the ability to handle functional inputs for the dictionary, and its
actual target. The coherence is important when choosing the relatively fast computation times [30]. The error parameters are
discretization of lt for the dictionary because it can decrease set using CV as in OMP.
the sparsity of the problem and make it difficult to differentiate The major disadvantage with C1 M is that it is an iterative
between two close target positions. method, requiring many more iterations than that of OMP. The
OMP is an iterative BP algorithm with an application- BPDN algorithm for C1 M performs the forward and adjoint
dependent stopping condition, which is typically set with cross operations many more times—sometimes on the order of hun-
validation (CV) or with a residual constraint [24]–[26]. The dreds to thousands of times. This sort of computational com-
stopping conditions used in this paper are set by CV and a plexity makes 3-D C1 M impractical with explicit dictionary
sparsity bound that does not allow the resulting images to enumeration, but the computational and storage improvements
have a sparsity below a certain level. OMP encourages sparse gained by exploiting shift invariance and using the FFT make it
answers; however, it does not guarantee a minimum-norm practical [12].
solution. The results will show that for the current application,
the performance is satisfactory. Having sparse solutions is an
B. Two-Dimensional Comparison to Previous Methods
advantage because it can create an image that does not require
any interpretation. For example, Fig. 7, if taken literally, would The easiest way to compare the new functional methods with
have many more than one target, so the user is responsible for older explicit methods is through small 2-D simulations. The 2-
interpreting the image blur as an artifact and not additional D simulations are used because 3-D simulations with explicit
targets. On the other hand, sparse solutions are sensitive to enumeration are impossible; they quickly scale beyond the
noise in the measurements. A three-target simulation is run storage capabilities of the computer used for the experiments.
and imaged with the OMP algorithm at different signal-to-noise For example, a 3-D problem, solved using the explicit matrix
ratio (SNR) levels. The three targets correspond to s([4, 8]) = method, with equal N discretizations, and a computer with
5, s([10, 12]) = 8, and s([22, 24]) = 15. Fig. 8 shows the mea- 4 GB of dedicated data memory could have a maximum N of
surements and results from the three-target setup. Fig. 8(c) is a approximately 28. Increasing the data memory size to 12 GB
perfect reconstruction of the original image, but once the noise only brings this value up to 33, and this assumes that the entire
overtakes the signal, it is difficult to get a perfect reconstruction, data memory is used only for Ψ . However, using the functional
as is seen by the missing target and false alarm at (40,40) cm in method enables a maximum N of 150 in the 4-GB computer
Fig. 8(d). and 196 in the 12-GB computer.
H
OMP is an iterative algorithm that applies gA and gA at each Fig. 10 compares the measured times for a 2-D GPR problem
step. For sparse solutions, the number of iterations is going to be run using BP and OMP, using different size spatial dimensions,
KRUEGER et al.: EFFICIENT ALGORITHM DESIGN FOR GPR IMAGING OF LANDMINES 4017

Fig. 8. Three-target simulation. (a) Time measurements. (b) Actual image. (c) OMP solution with 10-dB SNR. (d) OMP solution with −5-dB SNR.

Fig. 9. Images created with two random frequencies using (a) BP and
(b) C1 M.
H
Fig. 11. Approximately the ratio of storage and time saved by using gA and gA
instead of explicit matrix multiplication for different values of N .

Fig. 10. Comparing timing of BP with an explicit matrix to the functional

implementation (fBP) and OMP using an explicit matrix to the functional
implementation (fOMP). Fig. 12. Timing comparison of functional and explicit C1 M for different
numbers of compressed measurements.


N = 5 Nls Nlt , and fixed Nω = 401 frequencies, to show how The functional and explicit models are also compared using
the algorithms scale as the dimensionality of the image space a C1 M algorithm in Fig. 12, and again, the functional method
increases. Functional BP (fBP) is around 15 times faster than is much more time efficient by about a factor of 5 at 0.2×106
BP for N = 70. Functional OMP (fOMP) is also around 15 measurements and a factor of 15 at 106 measurements.
times faster than OMP and is actually about four times faster
than BP for N = 70. These values are again just for the
C. Three-Dimensional Compressed Simulations
2-D inversion where there is only one dimension containing the
redundancy; in 3-D, there are two. The slow data-acquisition time of SFGPR has led to height-
To see what the gains would be for the 3-D problem, the ened interest in algorithms that use compressed measurements.
approximate savings ratios for time and storage are calculated The following simulated example will only use compressed
for the 2-D and 3-D solvers in Fig. 11. To check how close algorithms. Compressed algorithms are those in which a sam-
the approximation is to reality, the ratios for experimental pling scheme is applied to the measurements and the dictionary
processing for N = 70 came out to be about 15 times, Fig. 11 before use in the algorithm. For this application, the same ran-
shows about 17 times for the 2-D processing at N = 70, which dom selection of stepped frequencies is taken at each individual
is acceptably close. This means that for a 3-D problem with scan position, i.e., ls , along the synthetic aperture.
size N = 70 that is solved with one of the discussed algorithms, The simulation run was on a problem size that would be rea-
there would be a reduction by a factor of approximately 600 to sonable in a real application. The image grid has a resolution of
the time and storage by using the functional algorithms. 2 cm and has 60 equal discretizations for xt , yt , and zt , giving
4018 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO. 7, JULY 2015

Fig. 13. Images created from a simulated environment with (a) CBP with 5% threshold, (b) CBP with 10% threshold, (c) COMP, and (d) C1 M.
TABLE I a suboptimal approach was taken by solving small 2-D slices
T IMING S TATISTICS FOR 3-D S IMULATION
and concatenating them together to form a 2.5-D image [12]. If
we compare the two approaches, it is easy to show why solving
the large 3-D problem is important. It is important to note that
the functional method can be done in 2-D as well and would
be more time efficient than using the explicit method, as was
a total number of image pixels equal to Nlt = 603 = 216, 000. shown in Section III-B.
As previously mentioned, the scan positions ls = (xs , ys , zs ) = This experiment consists of randomly placing two targets
(xt , yt , 0), or α = 0 and Δx = 2 cm, allow for the translational in a 3-D volume and using C1 M to invert the measure-
invariance to be present. The possible frequencies, i.e., ω, ments and create an image. The frequency range used is from
used in the data collection were 379 frequencies in the range 0 MHz to 5.02 GHz with Nω = 158 equally spaced frequencies.
2π(500 MHz) to 2π(8.06 GHz). The compression takes place The scan positions correspond to a colocated transmitter and
in the selection of a random set of the same ten frequencies at receiver. The transmitter locations are uniformly spaced in a
each ls . The total explicit dictionary size for this problem would 2-D square from ys = −96 cm to 90 cm at 6-cm intervals,
be Nls Nω ×Nlt = 602 (10)×603 , which turns out to be ap- giving Nys = 32, and, similarly, from xs = −96 cm to 90 cm
proximately 60 GB. A typical personal computer has less than at 6-cm intervals, giving Nxs = 32. The horizontal target lo-
12 GB of memory; hence, directly solving this problem with cations are in the range yt = −90 cm to 84 cm at 6-cm
an explicit dictionary matrix is impossible without calculating intervals, giving Nyt = 30; likewise for xt . The fact that the
the responses during run time in a loop. However, using the scan positions and the target locations have the same horizontal
translational invariance exploit removes the need to store xt and discretization leads to the translational-invariance simplifica-
yt , and the size of the dictionary becomes O(Nls )Nω ×Nzt = tion. The depth locations need not be constrained, but for this
experiment, z ranges from 270 cm to 420 cm deep at 6-cm
602 (10)×60, which is approximately 17 MB, a reduction by a
intervals, giving Nzt = 26. The number of compressed mea-
factor of more than 103 ! The simulation environment consisted
surements, i.e., M , is selected between 0.1% and 2.4% of the
of three point targets with an SNR of 15 dB. Fig. 13 shows
total measurements N = Nω Nls = 158×322 ≈ 105 . The num-
the results of the imaging with three different compressed
ber of discretizations in the target location parameter space is
algorithms. Two thresholds were used in the compressed BP
P = Nlt = 302 ×26 ≈ 104 . Without using CS or the functional
(CBP) case in Fig. 13(a) and (b), and it can be seen that
representation, the explicit matrix would be of size N ×P ≈
changing the user-defined threshold can dramatically affect
109 . However, employing CS and the functional representation,
the output images. The compressed OMP (COMP) and C1 M
N becomes M , P = 26, and the total storage requirements for a
algorithms in Fig. 13(c) and (d), respectively, achieve exact
matrix using 2.4% of the total measurements becomes M ×P =
reconstruction.
2400×26 ≈ 105 , a reduction of four orders of magnitude.
The measured timing statistics for this simulation can be
Fig. 14 shows the detection accuracy advantages that are
seen in Table I. A large portion of the time it takes to run the
available when solving the full 3-D imaging problem as op-
algorithms is associated with the time it takes to run the forward
posed to imaging the volume with many 2-D slices. When
and adjoint operators. The ratio of the amount of time per
using 2-D slices and not solving the full 3-D problem, the
function call is shown in the third column of Table I. The time
synthetic aperture in the xt dimension is not being directly
it takes to run COMP is almost exclusively dependent on the
H exploited. The sparsity in the 2-D slice image in Fig. 14(b)
speed of gA and gA , whereas a more complicated algorithm,
shows the resolution issues in xt . The full 3-D inversion that
i.e., C1 M, has a little more time devoted to overhead. For
is made possible by using the functional representations gA
algorithms such as C1 M, where the dictionary is going to H
and gA is shown in Fig. 14(c). The full 3-D inversion provides
be applied hundreds of times, any efficiency increase to the
a much higher resolution image in xt than using 2-D slices.
dictionary is going to linearly scale.
The higher resolution is a byproduct of the fact that the full
3-D inversion takes advantage of the full 2-D array aperture.
D. Comparing Full 3-D to Sliced 2-D
The accuracy metric used to analyze the solutions is the Earth
The computational inefficiency of previous methods has Mover’s Distance (EMD) [31]. The EMD is used because,
prevented the direct implementation of 3-D inversion; hence, unlike mean square error, EMD takes the support into account
KRUEGER et al.: EFFICIENT ALGORITHM DESIGN FOR GPR IMAGING OF LANDMINES 4019

Fig. 14. (a) Full 3-D C1 M solution using the FFT method with exact reconstruction. (b) Solution using 2-D slices. (c) EMD comparison of 2-D slice and full
3-D solutions.

TABLE II issues that can arise from that type of inversion. However,
S UMMARY OF L ANDMINE TARGET P ROPERTIES
the 2-D slice method could also be improved by exploiting
the translational invariance if it was desired to still be used.
It would be faster and more storage efficient than using the
explicit method to calculate each slice because each slice would
garner a computational time improvement to what was shown in
Section III-B. The purpose behind inverting these data again
is to show that a real-world experiment that was previously
too computationally intense to invert with full 3-D C1 M can
be completed with the exploitation of shift invariance. The
discretizations used for this problem are as follows: xs and ys
were both taken from −60 cm to 60 cm at 2-cm increments, xt
and yt taken from −58 cm to 56 cm at 2-cm increments, and
zt taken from 1 cm to 20 cm with 0.5-cm increments.
Fig. 15(c) and (d) shows the depth-summed top view of the
created images using uncompressed BP with 379 frequencies
when calculating the error. Targets with a small error in location and C1 M with 100 frequencies. Both of these were com-
have a lower EMD than they would with MSE. For example, puted with the functional implementation of the forward and
the MSE between two vectors with a single spike where the adjoint operators. It is very difficult to distinguish where some
spike location differs by one or by 100 is identical. However, of the weaker targets are in the BP image, but it becomes
the EMD will provide a higher error when the support is off by much clearer in the C1 M image. For example, the mines
100. An SNR analysis showing the EMD reduction produced at (xt , yt ) = (−45, 5) and (xt , yt ) = (0, 50), and the cylinder
by solving the full 3-D inversion instead of the 2-D slices is at (xt , yt ) = (45, 50). The imaging improvement of C1 M
found in Fig. 14(d). over BP combined with the fact that it can be done with
much fewer measurements is a substantial reason as to why
C1 M would be preferably used over BP. In the work by
IV. L ABORATORY E XPERIMENT
Gürbüz et al., full 3-D inversion for this experiment was
The final experiment uses data that are collected from the impossible with C1 M, but now, with the modification to the
system described by Counts et al. in [10]. The system consists storage and application of Ψ , it can be done. A problem with
of a linear array of resistive-vee antennas, microwave switches, this imaging in a CS application is the relatively low sparsity of
a vector network analyzer, and a 3-D positioner, under com- the resulting image, caused by using a point-target model. The
puter control. Data from a single transmit–receive pair spaced point-target model images each reflection as its own individual
12 cm apart as shown in Fig. 15(a) and measured at 379 equally target, which is why certain targets, such as that at (xt , yt ) =
spaced frequencies from 0.50 to 8.06 GHz are used in this work. (−50, 50), are larger than the others. The low sparsity in a
The targets are buried in damp compact sand, and the ground- CS problem can cause some very small targets to be excluded
truth location of the targets can be seen in Fig. 15(b). The because the accuracy of the inversion is based on the sparsity
depth of the targets is shown in the parenthesis on the graph. and the number of measurements used. This can be seen by
The landmines are inert, and their relative sizes are indicated the exclusion of the target at (xt , yt ) = (−45, 40). If a more
by their outlines. The properties of the buried objects for this accurate and complex model is used, the sparsity could be
experiment can be seen in Table II from Counts et al. [10]. dramatically increased, and even fewer measurements would
A 2-D slice C1 M algorithm was employed for inver- be required. Even with a more complex model, most likely, the
sion to avoid the computational inefficiencies in the work of translational invariance would remain, and the efficient model-
Gürbüz et al. [12], but we saw in the previous section the ing technique described in this paper could still be leveraged.
4020 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO. 7, JULY 2015

Fig. 15. Underground experiment for multiple objects. (a) Sensor setup. (b) Target locations. The values in parentheses correspond to the depths of the individual
targets. (c) Solution using BP. (d) Full 3-D CS solution using functional gA method.

With that being said, the point of the paper is not to validate R EFERENCES
the use of CS in this application. The point of the paper is to [1] S. R. J. Axelsson, “Analysis of random step frequency radar and compari-
show that an increase in the modeling efficiency can make it son with experiments,” IEEE Trans. Geosci. Remote Sens., vol. 45, no. 4,
more tractable, along with the other dictionary-based inversion pp. 890–904, Apr. 2007.
[2] L. P. Peters, Jr., J. J. Daniels, and J. D. Young, “Ground penetrating radar
algorithms. as a subsurface environmental sensing tool,” Proc. IEEE, vol. 82, no. 12,
pp. 1802–1822, Dec. 1994.
[3] D. J. Daniels, “Surface-penetrating radar,” Electron. Commun. Eng. J.,
V. C ONCLUSION vol. 8, no. 4, pp. 165–182, Aug. 1996.
[4] R. M. Lerner, “ Ground radar system,” U.S. Patent 3 831 173, Aug. 20,
An efficient restructuring of the storage and application of 1974.
the dictionary used in GPR inversion algorithms have been [5] L. J. Cutrona and E. N. Leith, “On the application of coherent optical
processing techniques to synthetic-aperture radar,” Proc. IEEE, vol. 54,
discussed. The dramatic reduction in computation times and no. 8, pp. 1026–1032, Aug. 1966.
storage constraints allow for previously unsolvable problems [6] W. M. Brown and L. J. Porcello, “An introduction to synthetic-aperture
to be solved. The model must be set up with a laterally homo- radar,” IEEE Spectr., vol. 6, no. 9, pp. 52–62, Sep. 1969.
[7] J. C. Kirk, “A discussion of digital processing in synthetic aperture radar,”
geneous medium, or the medium must be closely approximated IEEE Trans. Aerosp. Electron. Syst., vol. AES-11, no. 3, pp. 326–337,
as laterally homogeneous. The measurements need to be taken, May 1975.
or normalized, so that the radar height remains, or appears to [8] P. Millot and A. Berges, “Ground based SAR imaging tool for the design
of buried mine detectors,” in Proc. Detect. Abandoned Land Mines (Conf.
remain, relatively constant with respect to the surface. In this Publ. No. 431), 1996, pp. 7–9.
paper, only the point-target model is addressed; however, a [9] G. F. Stickley, “Synthetic aperture radar for the detection of shallow
direct improvement to this method, specifically for CS recovery buried objects,” Detection of Abandoned Land Mines, no. 431, pp. 7–9,
1996.
algorithms, would be to use a more complex target model to [10] T. Counts, A. C. Gurbuz, W. R. Scott, Jr., J. H. McClellan, and K. Kim,
more sparsely represent the buried objects. As long as the trans- “Multistatic ground-penetrating radar experiments,” IEEE Trans. Geosci.
lational invariance holds with a particular model, any target Remote Sens., vol. 45, no. 8, pp. 2544–2553, Aug. 2007.
[11] A. C. Gurbuz, J. H. McClellan, and W. R. Scott Jr., “Compressive sensing
model could be directly plugged into this model representation for subsurface imaging using ground penetrating radar,” Signal Process.,
and used efficiently. vol. 89, no. 10, pp. 1959–1972, Oct. 2009.
KRUEGER et al.: EFFICIENT ALGORITHM DESIGN FOR GPR IMAGING OF LANDMINES 4021

[12] A. C. Gurbuz, J. H. McClellan, and W. R. Scott Jr., “A compressive sens- Kyle R. Krueger received the B.S. degree in elec-
ing data acquisition and imaging method for stepped frequency GPRs,” trical engineering from the University of Texas at
IEEE Trans. Signal Process., vol. 57, no. 7, pp. 2640–2650, Jul. 2009. Austin, Austin, TX, USA, in 2008 and the M.S.
[13] X. Feng and M. Sato, “Pre-stack migration applied to GPR for landmine and Ph.D. degrees from the Georgia Institute of
detection,” Inverse Probl., vol. 20, no. 6, pp. S99–S115, Dec. 2004. Technology, Atlanta, GA, USA, in 2010 and 2013,
[14] C. J. Leuschen and R. G. Plumb, “A matched-filter-based reverse-time mi- respectively.
gration algorithm for ground-penetrating radar data,” IEEE Trans. Geosci. His research involves increasing the algorithmic
Remote Sens., vol. 39, no. 5, pp. 929–936, May 2001. efficiency of dictionary detection problems with ap-
[15] L. Song and Q. H. Liu, “Ground-penetrating radar land mine imag- plications to landmines.
ing: Two-dimensional seismic migration and three-dimensional inverse
scattering in layered media,” Radio Sci., vol. 40, no. 1, pp. RS1S90-
1–RS1S90-15 Feb. 2005.
[16] E. M. Johansson and J. E. Mast, “Three-dimensional ground-penetrating
radar imaging using synthetic aperture time-domain focusing,” in Proc. James H. McClellan (S’69–M’74–SM’79–
SPIE, Sep. 1994, vol. 2275, pp. 205–214. F’85–LF’12) received the B.S. degree in electrical
[17] M. Tuncer and A. Gurbuz, “Ground reflection removal in compressive engineering from Louisiana State University, Baton
sensing ground penetrating radars,” IEEE Geosci. Remote Sens. Lett., Rouge, LA, USA, in 1969 and the M.S. and Ph.D.
vol. 9, no. 1, pp. 23–27, Jan. 2012. degrees from Rice University, Houston, TX, USA,
[18] K. R. Krueger, J. H. McClellan, and W. R. Scott, Jr., “Dictionary reduction in 1972 and 1973, respectively.
technique for 3-D stepped-frequency GPR imaging using compressive From 1973 to 1982, he was a member of the Re-
sensing and the FFT,” in Proc. SPIE, Apr. 2012, p. 83650Q. search Staff at Lincoln Laboratory and then a Profes-
[19] K. R. Krueger, J. H. McClellan, and W. R. Scott, Jr., “3-D imaging for sor with the Massachusetts Institute of Technology,
ground penetrating radar using compressive sensing with block-toeplitz Cambridge, MA, USA. From 1982 to 1987, he was
structures,” in Proc. IEEE 7th SAM Signal Process. Workshop, Jun. 2012, with Schlumberger Well Services. Since 1987, he has
pp. 229–232. been a Professor with the School of Electrical and Computer Engineering,
[20] K. R. Krueger, J. H. McClellan, and W. R. Scott, Jr., “Sampling techniques Georgia Institute of Technology (Georgia Tech), Atlanta, GA, USA, where he
for improved algorithmic efficiency in electromagnetic sensing,” in Proc. presently holds the John and Marilu McCarty Chair. He is a coauthor of the texts
Int. Conf. Sampling Theory Appl., Jul. 2013, pp. 61–64. Number Theory in Digital Signal Processing, Computer Exercises for Signal
[21] P. A. Torrione, L. M. Collins, F. Clodfelter, S. Frasier, and I. Starnes, Processing, DSP First: A Multimedia Approach, and Signal Processing First,
“Application of the LMC algorithm to anomaly detection using the Wich- which received the IEEE McGraw-Hill Jacob Millman Award as outstanding
mann/NIITEK ground-penetrating radar,” in Proc. SPIE, 2003, vol. 5089, innovative textbook in 2003.
pp. 1127–1136. Dr. McClellan received the W. Howard Ector Outstanding Teacher Award
[22] J. Romberg, “Compressive sensing by random convolution,” SIAM J. at Georgia Tech in 1998 and the Education Award from the IEEE Signal
Imag. Sci., vol. 2, no. 4, pp. 1098–1128, Oct. 2009. Processing Society in 2001. In 1987, he received the Technical Achievement
[23] W. U. Bajwa, J. D. Haupt, G. M. Raz, S. J. Wright, and R. D. Nowak, Award for his work on finite-impulse-response filter design and, in 1996, the
“Toeplitz-structured compressed sensing matrices,” in Proc. Workshop Society Award, both from the IEEE Signal Processing Society. In 2004, he
Stat. Signal Process., Aug. 2007, pp. 294–298. was a corecipient of the IEEE Jack S. Kilby Signal Processing medal. He is
[24] J. A. Tropp and A. C. Gilbert, “Signal recovery from random measure- a member of Tau Beta Pi and Eta Kappa Nu.
ments via orthogonal matching pursuit,” IEEE Trans. Inf. Theory, vol. 53,
no. 12, pp. 4655–4666, Dec. 2007.
[25] A. C. Gurbuz, “Sparsity enhanced fast subsurface imaging with GPR,” in Waymond R. Scott Jr. (S’81–M’82–SM’03–F’08)
Proc. Int. Conf. GPR, Jun. 2010, pp. 1–5. received the B.E.E., M.S.E.E., and Ph.D. degrees
[26] P. Boufounos, M. Duarte, and R. Baraniuk, “Sparse signal reconstruction from the Georgia Institute of Technology, Atlanta,
from noisy compressive measurements using cross validation,” in Proc. GA, USA, in 1980, 1982, and 1985, respectively.
Workshop Stat. Signal Process., 2007, pp. 299–303. In 1986, he joined the School of Electrical and
[27] D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, vol. 52, Computer Engineering, Georgia Institute of Tech-
no. 4, pp. 1289–1306, Apr. 2006. nology, as an Assistant Professor, where he was
[28] R. G. Baraniuk, “Compressive sensing [lecture notes],” IEEE Signal Pro- subsequently promoted to Professor. His research in-
cess. Mag., vol. 24, no. 2, pp. 118–121, Jul. 2007. volves the interaction of electromagnetic and elastic
[29] E. J. Candès and M. B. Wakin, “An introduction to compressive sam- fields with materials. This research spans a broad
pling,” IEEE Signal Process. Mag., vol. 25, no. 2, pp. 21–30, Mar. 2008. range of topics, including the measurement of the
[30] E. van den Berg and M. P. Friedlander, “SPGL1: A Solver for Large- properties of materials, experimental and numerical modeling, and systems
Scale Sparse Reconstruction,” Jun. 2007, [Online]. Available: for the detection of buried objects. Currently, his research is concentrated on
http://www.cs.ubc.ca/labs/scl/spgl1. investigating techniques for detecting objects buried in the earth. This work has
[31] Y. Rubner, C. Tomasi, and L. J. Guibas, “A metric for distributions with many practical applications, for example, the detection of underground utilities,
applications to image databases,” in Proc. Int. Conf. Comput. Vis., 1998, buried hazardous waste, buried structures, unexploded ordnance, and buried
pp. 59–66. landmines.