Clutter distributions for tomographic image standardization in ground-penetrating radar
Brian M. Worthmann, David H. Chambers, David S. Perlmutter, Jeffrey E. Mast, David W. Paglieroni, Christian T. Pechard, Garrett A. Stevenson, Steven W. Bond
JJOURNAL OF L A TEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015 1
Clutter distributions for tomographic imagestandardization in ground-penetrating radar
Brian M. Worthmann, David H. Chambers,
Senior Member, IEEE,
David S. Perlmutter, Jeffrey E. Mast,David W. Paglieroni,
Senior Member, IEEE,
Christian T. Pechard, Garrett A. Stevenson,and Steven W. Bond,
Member, IEEE
Abstract —Multistatic ground-penetrating radar (GPR) signalscan be imaged tomographically to produce three-dimensionaldistributions of image intensities. In absence of objects ofinterest, these intensities can be considered to be estimates ofclutter. These clutter intensities spatially vary over several ordersof magnitude, and vary across different arrays, which makesdirect comparison of these raw intensities difficult. However, bygathering statistics on these intensities and their spatial variation,a variety of metrics can be determined. In this study, the clutterdistribution is found to fit better to a two-parameter Weibulldistribution than Gaussian or lognormal distributions. Basedupon the spatial variation of the two Weibull parameters, scaleand shape, more information may be gleaned from these data.How well the GPR array is illuminating various parts of theground, in depth and cross-track, may be determined from thespatial variation of the Weibull scale parameter, which may inturn be used to estimate an effective attenuation coefficient inthe soil. The transition in depth from clutter-limited to noise-limited conditions (which is one possible definition of GPRpenetration depth) can be estimated from the spatial variationof the Weibull shape parameter. Finally, the underlying clutterdistributions also provide an opportunity to standardize imageintensities to determine when a statistically significant deviationfrom background (clutter) has occurred, which is convenient forburied threat detection algorithm development which needs tobe robust across multiple different arrays.
Index Terms —ground penetrating radar, tomography, clutter,landmine detection, Weibull distribution.
I. I
NTRODUCTION G ROUND-penetrating radar (GPR) allows for buried ob-ject detection, such as landmines or utility pipes. GPRcan be used to detect spatial changes in dielectric permittivity,meaning both metallic and nonmetallic objects are detectable.In the system used for this study, a vehicle-mounted mul-tistatic radar array takes time series measurements, whichare used to form three dimensional images in real-time. Inthe absence of buried objects, the images are not empty, butrather contain clutter, where the word “clutter” is used here
The authors would like to acknowledge support from the Office of NavalResearch. This work was performed under the auspices of the U.S. Departmentof Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.B.M. Worthmann, D.H. Chambers, D.S. Perlmutter, D.W. Paglieroni,C.T. Pechard, G.A. Stevenson, and S.W. Bond are with Lawrence LivermoreNational Laboratory, Livermore, CA 94550 USA. (e-mails, respectively:[email protected], [email protected], [email protected],[email protected], [email protected], [email protected] [email protected]).Jeffrey E. Mast is with Teres Technologies, Inc., Loveland, CO 80537 USA(e-mail: [email protected]).Manuscript submitted 3 Aug 2020 to denote random volumetric scattering associated with spatialvariations in soil permittivity – in other words, homogeneousscattering associated with soil composition and structure, asopposed to scattering associated with discrete buried objectssuch as rocks or man-made objects. The amplitudes of thoseclutter voxels are then statistically analyzed to determinetheir deviation from the background. This study discusses thestatistics of those clutter voxels, including empirical fits toprobability distribution functions (PDFs), spatial variations ofthe distributions’ parameters, and how standardization can beperformed on image voxels to improve buried object visibility.
A. Clutter statistics
Most literature involving the fitting of clutter amplitudesor intensities to probability distribution functions (PDFs) isassociated with monostatic radars looking at weather-, land- orsea-based clutter, with the latter receiving the most attention.A variety of clutter distributions are considered. These includeone-parameter distributions, such as the zero-mean Gaussianmodel [1], [2] which was an adequate model for thermalnoise, and the Rayleigh distribution for amplitudes [3] (or chi-squared distribution for intensities [4]) which was found tobe suitable for high-grazing angle radar turns. A number oftwo-parameter distributions have been proposed for handlinglow-grazing angle returns, including the lognormal distribution[3], the Weibull distribution [5], [6], the Rice distributionfor amplitudes [7] (or noncentral chi-squared distribution forintensities [4]), the heavy-tailed Rayleigh distribution [8], andthe K-distribution [9]–[11] and similar generalizations [12].For additional flexibility, more parameters may be added, suchas with the six-parameter generalized compound distribution[13], which contains within it numerous special cases, includ-ing the generalized gamma distribution, the hypergeometricdistribution, and all of the other one- and two-parametersystems named here. The overview provided here is far fromcomplete, as there have been several decades of research intoradar clutter modeling [14]. It is also noted that some of thesedistributions have physical interpretations for their parameters,while others are simply empirical fits. In this study, onlythe two-parameter lognormal and Weibull distributions areconsidered, for simplicity. Some physical interpretation of theshape and scale parameters is provided, but the parametersare fundamentally an empirical fit. This is in agreement withthe claim by Shnidman that the Weibull and log-normaldistributions are used because “they produce a reasonable fit to a r X i v : . [ phy s i c s . g e o - ph ] J a n OURNAL OF L A TEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015 2 measured clutter densities, rather than any underlying physicalargument” [4].There are other investigations into clutter associated withground penetrating radar, including physical modeling of soilscattering [15], a graphical method to determine goodness-of-fit for different clutter probability distributions [16], and astudy on the relationships between complex permittivity andthe observed signal-to-clutter ratios [17]. These studies allutilize statistics gathered on time-series data directly, whereasthe study here utilizes ground penetrating radar imaging data,formed from an application of computed tomography on themeasured time-series data.While the GPR system utilized for this study is downward-looking, many characteristics are shared with forward-lookingGPR systems. Computed tomography images associated withforward-looking GPR are detailed in [18], [19] and [20].Several studies (e.g. [21], [22] and [23]) use electromagneticsimulations to study the effects of clutter originating fromrough surface scattering, and present methods to reduce thattype of clutter. A detailed study by Liao and Dogaru [24]showed the quantitative effect of rough surface scattering, anddetermined that scattering amplitudes are well-described by K-distributions, or Rayleigh distributions at low surface rough-ness. Another forward-looking GPR imaging study [19] foundclutter amplitudes follow a truncated Rayleigh distribution,and targets follow a three-component Gaussian mixture model.In this study, clutter intensities are empirically found to well-described by lognormal and Weibull distributions instead (seeSection II-A).
B. Multistatic tomographic imaging for GPR
Many commercially available GPR systems operate in abistatic or multi-monostatic configuration, where only onetransmitter and one receiver are active at a time. In the vehicle-based system [25] used for this study, a multistatic array isused, where only one transmitter (Tx) is active at a time, butall receivers (Rx) are taking measurements, as illustrated inFigure 1. The array used here has N = 16 transmitters andan equivalent number of receivers located at approximatelythe same position in cross-track and are separated by asmall distance in down-track. After cycling through all 16transmitters, = 256 time series measurements are stored,representing one frame. This is an increase in the amount ofreceived data by a factor of when compared to a multi-monostatic configuration. These additional measurements, andtheir associated spatial diversity, permit improved signal-to-noise ratios compared to a multi-monostatic system.Each frame’s time series are pre-processed to remove thecoupling pulse, radio interference, and the reflection fromthe ground, leaving only the radar returns from subsurfacefeatures, such as buried objects or volumetric scattering (clut-ter). Note that imperfect removal of the surface reflectionwould also lead to clutter, but this study does not distinguishbetween origins of clutter (i.e. surface roughness or volumetricscattering). The pre-processed time series serve as inputsto a plane-to-plane backpropagation algorithm [26], whichleverages temporal and spatial Fourier transforms to focus the … … … Tx
Fig. 1. Schematic illustrating the multistatic transmitter and receiver archi-tecture. There are N transmitters (in blue) and N receivers (in green). Theblack horizontal lines represent the ground, and only the arrows represent thereflections from the ground, though for clarity, other reflections, such as thosefrom subsurface objects or clutter are omitted for clarity of the figure. Cross Track: 𝑥 Depth: 𝑧 Along-track: 𝑦 Frame (a) (b)Fig. 2. Left: Schematic for the spatial convention used here, with x serving asthe cross-track variable, y serving as the along-track variable, and z servingas the depth variable (positive oriented into the ground). Right: A photographof the vehicle with the GPR array mounted to the front of the vehicle. Thearray is nominally 30 cm above the surface. received signal amplitudes to a uniformly-spaced spatial grid.The focused tomographic image serves as the data utilized forthis study, as well as the input to detection and classificationalgorithms [25], [27]. The spatial convention used here isillustrated in Figure 2, along with a photograph of the vehiclewith the GPR array mounted to the front.In the construction of the images used in this study, not allof the time series are utilized. Visualizing the time seriesas a × matrix of time series, if only the diagonal wereused, that would be comparable to a multi-monostatic array.We found empirically that utilizing matrix elements withinfour rows or columns of the diagonal gave better imagingperformance. In other words, using the language of reference[27], the multistatic degree is set to four (where a multistaticdegree of zero corresponds to the multi-monostatic case). Thischoice is primarily limited by the beamwidth of the antenna,with a wider beamwidth supporting a higher multistatic degree.Given the receivers are spaced at 15 cm, and the nominal arrayheight above the ground is 30 cm, a multistatic degree offour is equivalent to using surface reflections of 45 degrees orless. The antenna used in this study has an approximate half-width half-max beamwidth of 20–25 degrees, which means OURNAL OF L A TEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015 3
Fig. 3. Plot of the ultra-wideband pulse shape, amplitude (in arbitrary units)versus time, in nanoseconds. The inset at the bottom-right shows a photographof the RLVD antenna used. the far off-diagonal elements receive approximately 12 dB lesspower (or 25% of the amplitude) as the mononstatic (diagonal)elements receive.Since the backpropagation algorithm operates in the fre-quency domain, only certain frequencies (100 MHz to 3.5GHz) are imaged in this study. The synthetic aperture radar(SAR) algorithm detailed in Chambers et al. [25] was turnedoff for the purposes of this study, so that magnitude trends indepth could be assessed without the additional complicationof along-track summations that vary in depth. The images thatare created are 2.4 meters (m) (128 pixels) in cross-track ( x ),90cm (46 pixels) in depth ( z ) and a variable number of framesalong-track ( y ), with a frame spacing of 2 centimeters (cm).In depth, the first 10 cm correspond to the air gap with arefractive index of 1. The remaining 80 cm correspond to thesoil, with an assumed refractive index of 2.0. C. Dataset overview
In this study, the data from only one array is used. Specif-ically, a resistively-loaded V-dipole (RLVD) antenna is used(pictured in Figure 3), along with its monostatic pulse shapefor ground reflections. The center frequency of the pulseis approximately 500 MHz, though all frequencies between100 MHz and 3.5 GHz are utilized in the plane-to-planebackpropagation imaging algorithm. The same antenna isused for transmission and receiving, meaning there are 32antennas inside the array pictured in Fig. 2b. The 16 receivingantennas are spaced 15 cm apart in the cross-track directionfor a physical array aperture of 2.25 m. The 16 transmittingantennas have the same spacing, and are offset from the rowof receivers by 26.5 cm in the along-track direction. Due tothe antenna beamwidth, the imaging aperture of the array isfound empirically to be approximately 2.4 m.The data used here were taken across two campaigns(separated by one year) to a southwestern United States GPRtesting location. A total of 12 lanes at this location wereimaged by this array, representing approximately 9.3 km (or5.5 acres) of unique earth. Each run (a single pass over asingle lane) is comprised of thousands of frames—on average,40,000 frames per run. Each frame represents 2 cm (1 voxel)in the along-track direction, 2.4 m (128 voxels) in the cross-track direction, and 0.9 m (46 voxels) in the depth direction.A total of 189 runs were used in this study, representing 45billion individual voxels, or approximately 150 km (or 90 acres) of imagery. Due to the spatial sparseness of objectsthat were intentionally buried in these lanes, the overwhelmingmajority of this tomographic image data can be consideredto be measurements of the background, or clutter. Thus, forpurposes of gathering statistics on the clutter, the first andsecond moments of the distributions are not expected to beinfluenced by the presence of buried objects, as these representless than 1% of the image volume. While these buried objectshave little effect on the first and second moments of thedistribution, they may cause deviations from the anticipateddistribution for image amplitudes above the 99th percentile –see discussion in Sec II-D.
D. Methodology Overview
For clarity, we summarize below the processing steps ap-plied to the recorded data which serve as the input data forthe analyses given in the remaining sections of this paper.1) Fully multistatic ( × ) time series are recorded every2 cm.2) The coupling pulse is removed from the time series.3) The surface reflection is found and subtracted from thetime series.4) The time series are re-aligned to a nominal height of10cm above the surface.5) Plane-to-plane backpropagation is performed on thesetime series.6) The image is formed through the summation across arange of frequencies (given by the bandwidth) and arange of transmitter/receiver pairs (given by the multi-static degree)7) Image magnitudes are stored as a function of depth,cross-track, and along-track.8) Repeat steps (1) through (7) for each of 189 runs.In Section II, these image magnitudes are fit to probabilitydistributions. In Section III, physical interpretations of theprobability distribution parameters are provided. In Section IV,these image magnitudes are standardized to suppress clutterand accentuate buried targets, and Section V provides theconclusions drawn from this study.II. P ROBABILITY D ISTRIBUTIONS
A. PDF and CDF definitions
The data described in the previous section are fit to twodistributions, the lognormal and Weibull distribution. Theirdefinitions are given below: P Lognormal ( X = x | µ, σ ) = 1 √ πσ exp (cid:32) − (ln x − µ ) σ (cid:33) (1) P Weibull ( X = x | λ, κ ) = κλ (cid:16) xλ (cid:17) κ − exp (cid:16) − (cid:16) xλ (cid:17) κ (cid:17) (2)Thus, µ and λ can be seen as “scale” parameters, while σ and κ can be seen as “shape” parameters. X is therandom variable containing the clutter amplitudes, and x is the OURNAL OF L A TEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015 4 corresponding dummy variable. The cumulative distributionfunction of each function is given by:CDF
Lognormal ( X = x | µ, σ ) = 12 (cid:18) erf (cid:18) ln x − µ √ σ (cid:19)(cid:19) (3)CDF Weibull ( X = x | λ, κ ) = 1 − exp (cid:16) − (cid:16) xλ (cid:17) κ (cid:17) (4) B. Distribution parameter estimates
To estimate lognormal distribution parameters, the followingequations are used, where E ( · ) represents an expected value: µ = E (ln X ) (5) σ = (cid:113) E ((ln X ) ) − ( E (ln X )) (6)To estimate the Weibull distribution parameters, the firstand second moments of the clutter amplitude data, E ( X ) and E ( X ) can be used: Γ (cid:0) κ (cid:1) Γ (cid:0) κ (cid:1) = E ( X ) E ( X ) (7) λ = E ( X )Γ (cid:0) κ (cid:1) (8)In Eqn. 7, it is not possible to analytically evaluate κ directly, however the function on the left is well behaved,so a simple numerical interpolation scheme is sufficient todetermine κ , given the right-hand side ratio of moments.For some intuition for the dataset, a histogram of the voxelmagnitudes for all 45 billion voxels is given in Figure 4.Note that this is a log-log plot. On a linear vertical axis,the distribution would simply appear bell-curve shaped, andcould be interpreted to be a Gaussian distribution (that is,a lognormal distribution, since the horizontal axis is log-arithmic). However, a lognormal distribution would give asymmetric parabola on a log-log plot, and it is clear to seein Figure 4 that the distribution is skewed right, with the lefttail decidedly linear, not parabolic. This is in good agreementwith a Weibull distribution, where in the X (cid:28) λ limit,the distribution should appear linear on a log-log plot. Theright tail is where buried objects exist, and it is not obviousfrom this plot where the transition from buried object toclutter occurs. Of course, as with any detection problem, thesedistributions overlap. However, due to the low prevalence ofburied objects, only the extreme right edge of the plot wouldbe contaminated with buried objects. The additional structureto the plot beyond a magnitude of ∼ is likely due to thestrong spatial relationship the scale parameter has with space,or in particular, depth.For the purposes of evaluating the two fit parameters forboth distributions, each run, cross-track, and depth positionare treated as statistically independent of each other. In otherwords, there are 1.1 million different shape and scale pa-rameters for both distributions (189 runs ×
128 cross-trackpositions ×
46 depth positions = ∼ Fig. 4. Log-log plot of the histogram for all 45 billion voxels (with histogrambins spaced evenly in log space).
C. Spatial variation of distribution parameters
The histogram parameter plots in Fig. 5 hide the spatialvariation of the scale and shape parameters. To illustrate those,Figure 6 shows the cross-track and depth variation of the fourdistribution parameters’ values averaged across the 189 runs.In all six plots, the values are weighted by the number offrames in each run (so that longer runs are weighted morehighly than shorter runs). The possible physical interpretationsfor the results in Fig. 6 are reserved for Section III.
D. Goodness of fit
The final consideration that should be made here is thedetermination of which distribution fits the data best. Toaddress that, a modified version of a Q-Q (quantile-quantile)plot is created, as shown in Figure 7. To generate Fig. 7,a set of linearly spaced Z -score values are chosen, whichcan be converted into their CDF equivalents by way of erf (cid:0) √ Z (cid:1) = 2 · CDF, which is the relationship between Z -score and CDF for a normal distribution. Then, using eachrun and ( x, z ) position’s best fit distribution parameters, thoseCDFs are mapped to a clutter magnitude. Next, the empiricaldistribution is estimated for each clutter magnitude—or inother words, what percent of the data for that run and ( x, z ) position is less than or equal to the given clutter magnitude.This empirical CDF is evaluated for every clutter magnitude(one for each Z -score being sampled). Then a parametric plotof the empirical CDF vs the theoretical CDF is created foreach run and ( x, z ) position, which is effectively a Q-Q plot.If the horizontal and vertical axes scaled linearly with CDF,then the rare (infrequent) values that help determine the fit ofthe distribution would be pushed into the extreme edges ofthe plot, with most of the plot showing strong correspondencebetween median values of distribution. However, by scalinglinearly with Z -score, as these plots do, these rare values aregiven equal significance in the plot as the common values nearthe median. In other words, the middle 80% is given the samespace on the horizontal (and vertical) axis as the lower 10%and the upper 10%. Similarly, the 99% to 99.9% (or 1% to0.1%) of the data takes up approximately the same space onthe plot as one of the middle quartiles ( e.g.
25% to 50%).
OURNAL OF L A TEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015 5 (a)(b)(c)(d)Fig. 5. Histogram of best-fit parameters for the lognormal distribution ( µ and σ are Fig. 5a and 5b respectively), and the Weibull distribution ( σ and κ areFig. 5c and 5d respectively) This plot can be generated for all 1.1 million data sets[one for each of the runs, and ( x, z ) positions]. To avoidplotting 1.1 million lines, instead the mean (solid coloredlines), the ± standard deviation (darker-shaded regions) andthe maximum/minimum (lighter-shaded regions) are shown.The closer the lines and shaded regions are to the blackdiagonal line, the better the fit. If the shaded regions fallabove the black line, this means the observed distribution hasa fatter tail than the proposed distribution. Similarly, fallingbelow the black line means the observed distribution has athinner tail than the proposed distribution. Additionally, toquantify the goodness of fit, the Kolmogorov-Smirnov (KS)statistic is generated, which is effectively a measure of thelargest difference between the empirical CDF and the testdistribution’s CDF for any quantile [28], where a smaller KSstatistic implies a better fit to the proposed distribution.These Q-Q plots support the following conclusions:1) The observed clutter distributions have consistently fat-ter tails on the left (low clutter magnitudes) than wouldbe expected for a lognormal distribution.2) The observed clutter distribution has thinner tails on theleft (low clutter magnitudes) than would be expected fora Weibull distribution, but the match is quantitativelybetter than that of the lognormal distribution.3) Both plots show significant deviations above the 99%CDF. However, this is where buried objects are expectedto reside, and therefore a mismatch in that region isnot surprising, as image magnitudes for buried objectswould be expected to followed a different distributionthan for clutter, and would likely be strongly dependenton the buried objects’ scattering physics.4) The average KS statistic (averaged over each run andeach cross-track and depth position) for the Weibulldistribution (0.0312) is about half the value of theKS statistic using the Lognormal distribution (0.0555),which implies that a Weibull distribution is a better fitto the data than Lognormal.Additionally, it can be seen qualitatively that the Weibullplot in Fig. 7a is closer to linear (particularly below the99th percentile) than the Lognormal plot in Fig. 7b, whichhas a consistent curve away from the black line. For allof these reasons, it is claimed that the Weibull PDF is abetter fit to the observed distribution than the Lognormal PDF.Though, it is also important to point out, that for the middle(which represents ∼ % of the observed distribution), bothdistributions fit reasonably well. Only at the more extremevalues in the left tail are deviations appreciable.III. P ARAMETER I NTERPRETATIONS
This section considers the spatial variation of the distribu-tion parameters given in the previous section. More specifi-cally, it delineates what the spatial variation means in physicalterms about either the array or the environment.
A. Attenuation coefficient
A plot of the mean linear amplitude as a function of depthis given in Figure 8. The linear amplitude is averaged across
OURNAL OF L A TEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015 6 (a) (b)(c) (d)Fig. 6. Cross-track and depth variation of various moments and parameters. The left column shows scale parameters, with the top being the fitted µ ( x, z ) parameter for the lognormal distribution, and the bottom being the logarithm of the fitted λ ( x, z ) parameter for the Weibull distribution. The right columngives the shape parameters, with the top being σ ( x, z ) for the lognormal distribution, and the bottom being κ ( x, z ) for the Weibull distribution. All fourplots have a horizontal axis of 2.4m and a vertical axis of 0.9m, with a black horizontal line indicating the ground surface. The color axes are indicated tothe right of each plot. all runs (weighted by the number of frames in a run), and thenaveraged in cross-track.Theoretically, this plot of power vs. depth would giveamplitudes near the noise floor above the surface, and thenwould decay linearly with increasing depth until returningto the noise floor. Roughly, that trend does hold, but onlyapproximately so. Above the surface, the amplitudes here arelikely above the would-be noise floor, as this part of the imagecorresponds to early arrival times, which is where the couplingpulse dominates the signal. Imperfect removal of the couplingpulse can lead to higher amplitudes in this air gap region.The first 20 cm of soil shows an unexpected peak inamplitude. The most likely culprit for this is imperfect removalof the ground bounce, the removal of which is known tobe difficult on non-flat ground, as well as the imperfectlyimpulsive source waveform (see Fig. 3, which shows someringing after the main pulse). However, from about 20 – 50 cmin depth, the plot is roughly linear, with a best-fit slope of 24dB/m, which can be considered to be an indirect measurementof the attenuation coefficient. This estimate ignores otherforms of losses, such as spherical spreading, or the beampattern of the antennas, so 24 dB/m could be considered anupper bound on the attenuation coefficient. Nonetheless, thisvalue is reasonable, given the the frequency content of thepulse and the geology of the test location (a hot desert regionwhose soil type is primarily loamy sand, containing 2% orless of silt or clay). B. Cross-track beam pattern
In Figure 9, the mean amplitude is given versus cross-trackposition, where the mean is a weighted average across runs, and averaged in depth. The small stair-steps visible in thisplot are reflected in the vertical striations seen in the leftcolumn of Fig. 6, and are associated with a discrete changein the number of transmitters/receivers contributing to a givenvoxel ( i.e. imposing the multistatic degree with a hard cut-off, rather than a smooth transition). Despite these small stair-steps, cross-track trends can be seen, which can be attributed tothe cross-track beam pattern of the antennas. This is again anindirect measurement, but it shows that voxels on the extremeedges are “illuminated” at about half the intensity (6 dB ofpower) as an equivalent voxel at the center. There are likelytwo contributions to the cross-track amplitude variations: (1)the beam pattern of the individual antennas, and (2) the useof multi-static data with a finite aperture.
C. Transition from clutter-limited to noise-limited
Consider a plot of the average Weibull shape parameteras a function of depth, as shown in Figure 10. A value of κ = 1 corresponds to an exponential distribution, whereasa value of κ = 2 corresponds to a Rayleigh distribution.A Rayleigh distribution is interesting because that would bethe distribution one would anticipate for the magnitude of acomplex number whose real and imaginary parts are zero-mean Gaussian distributed. The expectation is that noise, i.e. thermal noise, would produce this zero-mean Gaussiandistribution for the real and imaginary parts, leading to κ = 2 for the best-fit Weibull distribution. Thus, for deeper voxels,the best-fit Weibull distribution is getting closer to a Rayleighdistribution, which would be the noise-limited case. However,it is not possible to claim the exact depth at which the voxelmagnitudes are no longer clutter-limited, but instead noise- OURNAL OF L A TEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015 7 (a)(b)Fig. 7. A modified Q-Q plot that allows the determination of whichdistribution fits the data best, with the Lognormal distribution in 7a, andWeibull distribution in 7b. The black diagonal line is the line of perfect matchbetween the empirical and theoretical distributions. The blue/red lines indicatethe average Q-Q plot across all runs and ( x, z ) positions. The darker-shadedregions represent ± standard deviation across all runs and ( x, z ) positions.The lighter-shaded regions represent the maximum and minimum valuesobserved. Note that the plots are shown with linear sampling in the Z -scoreequivalent axes, labeled across the top and right, though the correspondingCDF percentages are labeled on the left and bottom axes. Fig. 8. The average power, in decibels (dB), averaged cross-track, versusdepth, in centimeters. The decibel scale is relative – no milliwatt reference isimplied in this figure. The red dashed line is a line of best fit between 20 –50 cm of depth, which has a slope of approximately 24 dB/m, which can beinterpreted as an indirect measure of the attenuation coefficient for this soil.Fig. 9. Average power versus cross-track, averaged across runs and in depth limited. For example, κ = 1 . might truly still be clutter-limited, though it cannot be confirmed without conductingsimulations or more detailed analyses of this data. But if oneuses a nominal cut-off of κ = 1 . (being halfway between theextreme cases of a Rayleigh distribution and an exponentialdistribution), then a penetration depth of approximately 35 cmcan be estimated. Note: this does not mean targets cannot bedetected deeper than this depth. What it does suggest is thata more powerful transmitter (or lower-noise receiver) wouldlikely push this transition point deeper. Or, alternatively, if theregion of interest is no deeper than 35 cm, a more powerfultransmitter would not be expected to improve performance.IV. GPR I MAGE S TANDARDIZATION
Given knowledge of the underlying clutter distributionsand their spatial variation, this information can be used totransform voxel magnitudes from their raw values into anothervalue which more accurately reflects their deviation from thebackground. To do this, the CDF for each voxel can beestimated based on the underlying distribution. Plotting CDFsdirectly, however, is inconvenient, as the most interesting datapoints would occur in the range between 0.9 and 1.0. However,
OURNAL OF L A TEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015 8
Fig. 10. Plot of the average Weibull shape parameter versus depth, wherethe shape parameter is a weighted average across runs, and an average incross-track. κ = 1 corresponds to an exponential distribution, and κ = 2 corresponds to a Rayleigh distribution. calculating a Z -score on logarithmic data allows brightnessvalues to be compared more intuitively, on a ‘linear’ axis,much like the rescaled axes in Fig. 7. This transformed imageis termed as the “standardized logarithmic intensity” image,or SLI image. This is the motivation behind image standard-ization – the exact implementation of which is considered inthe next subsection. A. Standardized Logarithmic Intensity (SLI) definition
Consider the raw tomographic image output R ( x, y, z ) . Thestandardized logarithmic intensity (SLI) image output is givenby: SLI ( x, y, z ) = log R ( x, y, z ) − µ ( x, z ) σ avg + f airgap ( z ) (9)The logarithm shown can be whichever base, as long as the µ and σ used are calculated on that same base (the plots in Fig.6 are generated with base e ). µ ( x, z ) is the lognormal mean,which is a function of cross-track and depth, as given in Fig.6a. σ avg is in fact a scalar, and is the mean of Fig. 6b (averagedin quadrature in depth and cross-track). Note that, for real-time systems, µ and σ can be updated for each down-trackposition, and thus could be considered functions of y as well.Finally, f airgap ( z ) is defined as − | z | / (10 cm ) for z < , andzero for z > . This ad-hoc function was designed to suppressair gap imaging artifacts arising from pre-processing that areunintentionally accentuated in image standardization. Whilethis function is designed to suppress voxels in the airgap, itis also designed to have no effect on the subsurface voxels,which is where the most interesting image features occur.The choice in the above equation to use σ avg rather than σ ( x, z ) is somewhat arbitrary, but has its advantages. First, ina real-time implementation of this system, it is possible that σ ( x, z ) may take many frames to converge to an appropriatevalue, whereas σ avg would be expected to converge somewhatfaster. Additionally, σ ( x, z ) was found to be a weak functionof spatial position (see Fig. 6b), so the approximation that σ ( x, z ) ≈ σ avg is reasonable. Finally, it also preserves theinterpretation of SLI ( x, y, z ) being approximately equal toa normalization of R ( x, y, z ) instead of a standardization(where normalization refers to dividing by a constant, andstandardization refers to subtracting a mean and dividing bya standard deviation). Owing to the nature of logarithms,and ignoring the f airgap ( z ) function above (which is zero forsubsurface voxels anyway), SLI ( x, y, z ) can be re-written as:SLI ( x, y, z ) = log (cid:32)(cid:20) R ( x, y, z )exp( µ ( x, z )) (cid:21) /σ avg (cid:33) (10)In this form, it is more easily seen that SLI images aresimply logarithms of the normalized image, which is thenrescaled by an exponent of /σ avg . Thus, even though SLIis named as a standardized image, and its definition certainlyresembles a standardization process, it can also be thought ofas substantively equivalent to a normalized image, thanks tothe use of a constant σ avg , as opposed to a spatially varyingone. Additionally, keeping σ constant means only µ ( x, z ) canamplify or attenuate data, not σ as well.Also, it should be noted that the SLI definition utilizesmagnitudes, not intensities, though its name contains theword “intensity”. However, due to the logarithm, using trueintensities (proportional to magnitude-squared) would producea µ and σ that are a factor of two larger, and therefore thestandardization equation would cancel out this effect. B. Standardized images of targets
Four example targets are selected to illustrate the utility ofimage standardization. In Figure 11, four pairs of images areshown. On the left are raw images, shown after a logarithm isapplied, i.e. log ( R ( x, y, z )) , and on the right are SLI ( x, y, z ) images, standardized according to Section IV-A. In all fourraw image plots, the raw color scale is defined to be between13 (dark blue) and 16 (dark red). In all four standardized imageplots, the SLI color scale is set to be between 0.5 (dark blue)and 2.5 (dark red). Since the images are fundamentally three-dimensional volumes, a max projection across each of the threeaxes is used to display these cubes. Recall that x (cross-track)is the 240-cm-long dimension, y (along-track) is the 300-cm-long dimension, and z (depth) is the 90-cm-long dimension.The upper dark blue region seen in the x – z and y – z plotsshow the air gap, where the region is made darker blue thanexpected due to the use of f airgap ( z ) in the SLI definition.For these targets, the standardized images are generallysuperior when detecting buried objects, as the target cantechnically be seen in the raw images in all four cases, butthe presence of shallow clutter tends to obscure them (here,‘shallow’ refers to depths of a few wavelengths). By standard-izing with respect to the underlying clutter distribution, onlyvoxels with anomalously high magnitudes are accentuated,with the rest being suppressed. These four images illustratethat standardization is generally positive for improving the OURNAL OF L A TEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015 9 (a) (b)(c) (d)(e) (f)(g) (h)Fig. 11. Raw images (left column) and standardized images (right column) for four different targets: a circular plastic-encased landmine (first row), a large,square, plastic landmine seen at the edge of the array (second row), a deep large, metallic target (third row), and a shallow plastic-encased circular landminesurrounded by bright clutter (fourth row).
OURNAL OF L A TEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015 10 contrast between clutter and targets, which is convenient forthe development of an automated detection algorithm. Thetarget in the first row is fairly typical, in that the targetcan be seen in unstandardized image, but the standardizedimage (Fig 11b) makes it clearer. The other three targetexamples were chosen to accentuate a different aspect of thestandardization: compensating for cross-track variation (Fig11d), compensating for depth variation (Fig 11f), and thesuppression of shallow clutter to reveal a slightly deeper target(Fig 11h).While these figures show an improvement in contrast withstandardization, note that this apparent improvement in con-trast is primarily associated with compensating for the spatialvariation of the clutter distribution. Standardization by itselfdoes not increase intrinsic contrast, it simply measures contrastin a way that can be easily compared with other regions ofthe image, or with other systems.Standardizing in this way relies upon the underlying clutterdistribution being unimodal for a given ( x, z ) row, with well-defined means and variances. The unimodal requirement istypically satisfied, unless there are significant along-trackvariations in the environment that could lead to multimodaldistributions. The means and variances of multimodal distribu-tions are less meaningful for standardization, and could presentchallenges. However, if multimodal distribution occurred dueto the background environment changing discontinuously (ie.the vehicle abruptly transitioned from soil to concrete), therunning mean and variance calculation could be reset based onsome threshold. Other than this, only pathological clutter dis-tributions, such as the Cauchy-Lorentz distribution whose vari-ance is undefined, would present difficulties for this method ofstandardization. In the authors’ experience, this methodologyto standardize GPR images was found to be robust across avariety of testing locations and GPR imaging settings.V. C ONCLUSIONS
This study supports the following three conclusions:1) The underlying clutter magnitude distribution in GPRimagery is modeled well by a Weibull distribution ratherthan a log-normal distribution. This suggests that cluttermagnitudes in GPR imagery follows similar distributionsas sea-clutter gathered from synthetic aperture radars.2) The Weibull shape parameter can be used to estimatethe transition in depth from clutter-limited to noise-limited. In this study, objects at a depth of 35 cm orless are unlikely to be improved by a more powerfultransmitter, whereas deeper objects could benefit from amore powerful transmitter.3) Image standardization that appropriately compensatesfor attenuation and finite antenna beam-widths leadsto images with significantly improved clutter-to-targetcontrast, and allows for more fair comparison betweentargets seen at the edge of the array, deep targets, andshallow targets surrounded by bright clutter.A
CKNOWLEDGMENT
The authors would like to acknowledge support from theOffice of Naval Research. This work was performed under the auspices of the U.S. Department of Energy by LawrenceLivermore National Laboratory under Contract DE-AC52-07NA27344. R
EFERENCES[1] J. Marcum, “A statistical theory of target detection by pulsed radar,”
IEEE Transactions on Information Theory , vol. 6, no. 2, pp. 59–267,Apr. 1960. [Online]. Available: https://doi.org/10.1109/tit.1960.1057560[2] P. Swerling, “Probability of detection for fluctuating targets,”
IEEETransactions on Information Theory , vol. 6, no. 2, pp. 269–308, Apr.1960. [Online]. Available: https://doi.org/10.1109/tit.1960.1057561[3] G. V. Trunk, “Radar properties of non-rayleigh sea clutter,”
IEEE Transactions on Aerospace and Electronic Systems , vol.AES-8, no. 2, pp. 196–204, Mar. 1972. [Online]. Available:https://doi.org/10.1109/taes.1972.309490[4] D. Shnidman, “Generalized radar clutter model,”
IEEE Transactions onAerospace and Electronic Systems , vol. 35, no. 3, pp. 857–865, Jul.1999. [Online]. Available: https://doi.org/10.1109/7.784056[5] M. Sekine and Y. Mao,
Weibull radar clutter , 1st ed., ser. 3. London,UK: The Institution of Engineering and Technology, 1 1990, vol. 1.[6] D. Schleher, “Radar detection in weibull clutter,”
IEEE Transactionson Aerospace and Electronic Systems , vol. AES-12, no. 6, pp.736–743, Nov. 1976. [Online]. Available: https://doi.org/10.1109/taes.1976.308352[7] D. Drumheller, “General expressions for rician density and distributionfunctions,”
IEEE Transactions on Aerospace and Electronic Systems ,vol. 29, no. 2, pp. 580–588, Apr. 1993. [Online]. Available:https://doi.org/10.1109/7.210098[8] E. Kuruoglu and J. Zerubia, “Modeling SAR images with ageneralization of the rayleigh distribution,”
IEEE Transactions onImage Processing , vol. 13, no. 4, pp. 527–533, Apr. 2004. [Online].Available: https://doi.org/10.1109/tip.2003.818017[9] D. Iskander and A. Zoubir, “Estimation of the parameters of thek-distribution using higher order and fractional moments [radarclutter],”
IEEE Transactions on Aerospace and Electronic Systems ,vol. 35, no. 4, pp. 1453–1457, 1999. [Online]. Available: https://doi.org/10.1109/7.805463[10] J. Jao, “Amplitude distribution of composite terrain radar clutter andthe k-distribution,”
IEEE Transactions on Antennas and Propagation ,vol. 32, no. 10, pp. 1049–1062, Oct. 1984. [Online]. Available:https://doi.org/10.1109/tap.1984.1143200[11] K. D. Ward, R. J. A. Tough, and S. Watts,
Sea Clutter: Scattering, theK Distribution and Radar Performance , 1st ed., ser. 20. London, UK:The Institution of Engineering and Technology, 1 2006, vol. 1.[12] M. A. Ritchie, K. Woodbridge, and A. G. Stove, “Analysis of sea clutterdistribution variation with doppler using the compound k-distribution,”in . IEEE, 2010. [Online]. Available:https://doi.org/10.1109/radar.2010.5494569[13] Anastassopoulos, G. Lampropoulos, A. Drosopoulos, and N. Rey, “Highresolution radar clutter statistics,”
IEEE Transactions on Aerospaceand Electronic Systems , vol. 35, no. 1, pp. 43–60, 1999. [Online].Available: https://doi.org/10.1109/7.745679[14] S. Watts, “Radar sea clutter: Recent progress and future challenges,” in . IEEE, Sep. 2008. [Online].Available: https://doi.org/10.1109/radar.2008.4653882[15] K. Takahashi, J. Igel, and H. Preetz, “Clutter modeling for ground-penetrating radar measurements in heterogeneous soils,”
IEEE Journalof Selected Topics in Applied Earth Observations and Remote Sensing ,vol. 4, no. 4, pp. 739–747, 2011.[16] A. C. Gurbuz, “Determination of background distribution for ground-penetrating radar data,”
IEEE Geoscience and Remote Sensing Letters ,vol. 9, no. 4, pp. 544–548, 2012.[17] I. Giannakis, A. Giannopoulos, and A. Yarovoy, “Model-based eval-uation of signal-to-clutter ratio for landmine detection using ground-penetrating radar,”
IEEE Transactions on Geoscience and Remote Sens-ing , vol. 54, no. 6, pp. 3564–3573, 2016.[18] D. Comite, F. Ahmad, T. Dogaru, and M. G. Amin, “Adaptive detectionof low-signature targets in forward-looking gpr imagery,”
IEEE Geo-science and Remote Sensing Letters , vol. 15, no. 10, pp. 1520–1524,2018.[19] D. Comite, F. Ahmad, D. Liao, T. Dogaru, and M. G. Amin, “Multi-view imaging for low-signature target detection in rough-surface clutterenvironment,”
IEEE Transactions on Geoscience and Remote Sensing ,vol. 55, no. 9, pp. 5220–5229, 2017.
OURNAL OF L A TEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015 11 [20] D. Comite, F. Ahmad, T. Dogaru, and M. Amin, “Coherence-factor-based rough surface clutter suppression for forward-looking GPRimaging,”
Remote Sensing , vol. 12, no. 5, p. 857, Mar. 2020. [Online].Available: https://doi.org/10.3390/rs12050857[21] V. Kovalenko, A. G. Yarovoy, and L. P. Ligthart, “A novel clutter sup-pression algorithm for landmine detection with gpr,”
IEEE Transactionson Geoscience and Remote Sensing , vol. 45, no. 11, pp. 3740–3751,2007.[22] A. van der Merwe and I. J. Gupta, “A novel signal processing techniquefor clutter reduction in gpr measurements of small, shallow land mines,”
IEEE Transactions on Geoscience and Remote Sensing , vol. 38, no. 6,pp. 2627–2637, 2000.[23] T. Dogaru and L. Carin, “Time-domain sensing of targets buried undera rough air-ground interface,”
IEEE Transactions on Antennas andPropagation , vol. 46, no. 3, pp. 360–372, 1998.[24] D. Liao and T. Dogaru, “Full-wave characterization of rough terrainsurface scattering for forward-looking radar applications,”
IEEE Trans-actions on Antennas and Propagation , vol. 60, no. 8, pp. 3853–3866,2012.[25] D. H. Chambers, D. W. Paglieroni, J. E. Mast, and N. R. Beer,“Real-time vehicle-mounted multistatic ground penetrating radar imag-ing system for buried object detection,” Lawrence Livermore NationalLaboratory, Livermore, CA, USA, Tech. Rep. 1, 6 2011, LLNL-TR-615452.[26] C. Leuschen and R. Plumb, “A matched-filter-based reverse-timemigration algorithm for ground-penetrating radar data,”
IEEETransactions on Geoscience and Remote Sensing , vol. 39, no. 5, pp. 929–936, May 2001. [Online]. Available: https://doi.org/10.1109/36.921410[27] D. W. Paglieroni, D. H. Chambers, J. E. Mast, S. W. Bond,and N. R. Beer, “Imaging modes for ground penetrating radarand their relation to detection performance,”
IEEE Journal ofSelected Topics in Applied Earth Observations and Remote Sensing ,vol. 8, no. 3, pp. 1132–1144, Mar. 2015. [Online]. Available:https://doi.org/10.1109/jstars.2014.2357718[28] J. H. Drew, A. G. Glen, and L. M. Leemis, “Computing thecumulative distribution function of the kolmogorov–smirnov statistic,”
Computational Statistics & Data Analysis , vol. 34, no. 1, pp. 1–15,Jul. 2000. [Online]. Available: https://doi.org/10.1016/s0167-9473(99)00069-9
Brian M. Worthmann is a postdoctoral researcherat LLNL with a background in wave-based remotesensing and signal processing. He received his PhDat the University of Michigan in Applied Physics,where he worked on model-based source localizationin ocean acoustics. At LLNL, he has worked withground penetrating radar signal and image process-ing, as well as statistical data analysis for hardwareand performance characterization.
David H. Chambers
David H. Chambers(M’97–SM’04) is originally from southeast Kansas.He received the Bachelor’s degrees in both physicsand mechanical engineering, and the Master’sdegree in physics from Washington Universityin St. Louis, St Louis, MO, USA, from 1976 to1982, and the Ph.D. degree in theoretical andapplied mechanics from the University of Illinoisat Urbana-Champaign, in 1987. Afterward, he tooka position as a Physicist with Lawrence LivermoreNational Laboratory, Livermore, CA, USA, workingon the propagation of high-energy laser beams through the atmosphere andradar imaging of the ocean surface. His research interests include groundpenetrating radar, gravity gradiometry, statistical theory of fission chains,and applications of time reversal symmetry to target characterization andcommunications for both acoustic and radar imaging. Dr. Chambers is aFellow of the Acoustical Society of America and a member of the AmericanPhysical Society and Society of Industrial and Applied Mathematics. He haspublished papers in the IEEE Transactions for Antennas and Propagation,Journal of the Acoustical Society of America, Physical Review, and others.
David S. Perlmutter
David S. Perlmutter is a staffengineer at Lawrence Livermore National Labora-tory (LLNL) interested in applied signal processingand machine learning. He received a M.S. in elec-trical engineering at the University of Washington(’15), where he researched techniques for low-doseCT image reconstruction. Previously he was a staffengineer at MIT Lincoln Laboratory, working onLiDAR imaging technologies. At LLNL, his expe-rience includes modeling, control and data analysisof ultrafast optical and RF sensing systems.
Jeffrey E. Mast
Jeffrey E. Mast received the B.S.,M.S., and Ph.D. degrees in electrical and com-puter engineering from the University of Illinois atUrbana-Champaign, Champaign, IL, USA, in 1987,1989, and 1993, respectively. He was a General Mo-tors Scholar as an undergraduate and graduated withhonors. As a graduate student, he was a Fellow forthe U.S. Army Advanced Construction TechnologyCenter specializing in radar imaging systems fornondestructive evaluation of civil structures. From1993 to 1998, he held a position at Lawrence Liv-ermore National Laboratory (LLNL) in the Imaging and Detection Programdeveloping radar imaging systems for submarine detection, bridge deck androadway inspection, and unexploded ordnance detection. He is currentlyPresident of Teres Technologies, Inc., Loveland, CO, USA, and serves as aConsultant for LLNL and several private and public companies. His researchinterests include radar imaging and remote sensing, wide area motion imagery,high performance computing, and real-time systems.
David W. Paglieroni
David W. Paglieroni (S’80–M’87–SM’02) received the B.S. (summa cum laude),M.S., and Ph.D. degrees in electrical and computerengineering from the University of California, Davis,CA, USA, in 1982, 1984, and 1986. In 1987, hejoined Lockheed Martin, San Jose, CA, USA, wherehe served as a Manager of the Imagery SoftwareSection. Since 1999, he has been with the LawrenceLivermore National Laboratory (LLNL), Livermore,CA, USA, where he is currently a senior memberof the technical staff in the Engineering Directorate.Dr. Paglieroni provides subject matter expertise on a variety of nationalsecurity projects. He holds numerous patents and has published numerouspeer-reviewed papers in various journals and conferences. His current areasof interest and expertise include computer vision and AI-based data fusionfor risk assessment / situational awareness.
OURNAL OF L A TEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015 12
Christian T. Pechard
Christian Pechard received hisM.S degree in computer engineering from ESIEE-Paris, France, in 1995. From 1996 to 2007, heworked in the field of networking, telecommunica-tion, real-time embedded systems, network manage-ment at Cisco Systems, San Jose, California. From2007 to 2010, he worked in the field of enterprisewireless, real-time embedded system for TrapezeNetworks, Pleasanton, California. Since 2010, hejoined Lawrence Livermore National Laboratory(LLNL) to research and develop applications in thefield of telecommunication, wireless protocols, real-time embedded systems.He is currently responsible for software pipeline development on a large multi-years ground penetrating radar program.
Garrett A. Stevenson is a signal and image pro-cessing engineer at LLNL with a background in em-bedded deep learning, computer vision, nondestruc-tive characterization, and navigation/localization. Hereceived his M.S. in Computer Science at C.S.U.East Bay with emphases in Computer Vision andEmbedded Systems. He currently works in drugdiscovery and explosives detection using CT imagesin Ground Penetrating Radar and Checked Baggage.