Crustal Imaging and Characterization Team 

GP Workbench Manual: Technical Manual, 

User’s Guide, and Software Guide 

By Charles P. Oden and Craig W. Moulton 

Open-File Report 2006–1365 

U.S. Department of the Interior 

U.S. Geological Survey 
U.S. Department of the Interior 

DIRK KEMPTHORNE, Secretary 

U.S. Geological Survey 

Mark D. Myers, Director 

U.S. Geological Survey, Reston, Virginia 2006 

Revised and reprinted: 2006 

For product and ordering information: 

World Wide Web: http://www.usgs.gov/pubprod 

Telephone: 1-888-ASK-USGS

For more information on the USGS—the Federal source for science about the Earth, 

its natural and living resources, natural hazards, and the environment: 

World Wide Web: http://www.usgs.gov 

Telephone: 1-888-ASK-USGS 

Suggested citation: 

Oden, Charles P., and Moulton, Craig W., 2006, GP Workbench Manual: Technical Manual, 

User’s Guide, and Software Guide: U.S. Geological Survey Open-File Report 2006-1365, 81 p. 

Any use of trade, product, or firm names is for descriptive purposes only and does not imply 

endorsement by the U.S. Government. 

Although this report is in the public domain, permission must be secured from the individual 

copyright owners to reproduce any copyrighted material contained within this report. 

ii 
Contents 

Introduction ...................................................................................................................................1

Disclaimer .....................................................................................................................................3

GP Workbench Views ...................................................................................................................4

Menu and Toolbar Functions.........................................................................................................6

File Menu ...................................................................................................................................6

Import......................................................................................................................................6

Export......................................................................................................................................8

ALLTEM Data..........................................................................................................................8

Process Menu ..........................................................................................................................11

Grid Parameters....................................................................................................................11

Position Shift .........................................................................................................................15

Filter......................................................................................................................................15

Variogram..............................................................................................................................17

View Header..........................................................................................................................18

GPR Forward Model .............................................................................................................18

Migrate ..................................................................................................................................18

UXO Forward Model .............................................................................................................20

ALLTEM Calibration ..............................................................................................................24

Estimate UXO Parameters....................................................................................................26

View Menu ...............................................................................................................................32

VDL Colors............................................................................................................................32

Surfer ....................................................................................................................................33

Make JPGs............................................................................................................................33

iii 
GP Wave Utilities ..................................................................................................................33

Window and Help Menus .........................................................................................................33

The Toolbar..............................................................................................................................34

GP Wave Utilities ........................................................................................................................38

The GP Wave Utilities View ........................................................................................................38

Menu and Toolbar Functions.......................................................................................................40

File Menu .................................................................................................................................41

Import and Export..................................................................................................................41

Process ....................................................................................................................................43

Wave Plot, Filters and Operators ..........................................................................................43

Wave Plot, Subtract Ramp....................................................................................................45

Wave Plot, Subtract Waveform .............................................................................................45

Wave Plot, Remove DC ........................................................................................................46

Wave Plot, Deconvolve .........................................................................................................46

Wave Plot, Convolve.............................................................................................................47

Wave Plot, Correlate .............................................................................................................47

Wave Plot, Envelope.............................................................................................................47

Wave Plot, Hilbert transform .................................................................................................48

Wave Plot, Path Length ........................................................................................................48

Wave Plot, Amplitude Response...........................................................................................48

Wave Plot, Amplitude Linearization.......................................................................................49

Wave Plot, Packet Time and Amplitude ................................................................................49

Wave Plot, Compare Waves .................................................................................................50

Image Plot, Extract Trace......................................................................................................50

iv 
Image Plot, 3D Depth Shift....................................................................................................50

Image Plot, 2D Depth Shift....................................................................................................52

Image Plot, System Response a(w)T(k,w)R(k,w)..................................................................52

Radiation Pattern ..................................................................................................................55

Ground Properties.................................................................................................................56

Ground Properties Model ......................................................................................................58

Compare Waveforms ............................................................................................................63

TM and TE Spectrum ............................................................................................................66

View, Window, and Help Menus...............................................................................................67

The Toolbar..............................................................................................................................67

Data Processing Methods ...........................................................................................................68

GP Workbench GPR Data Processing Procedures..................................................................68

GPR Section View Processing..............................................................................................69

GPR Plan View Processing...................................................................................................69

GP Workbench ALLTEM Data Processing Procedures ...........................................................70

Calibration Procedure ...........................................................................................................71

UXO Target Inversion Procedure ..........................................................................................72

GP Wave Utilities Data Processing Procedures.......................................................................73

Determining the Transfer Function........................................................................................73

Synchronize and Compare Experimental Data with Simulated Data.....................................74

References Cited ........................................................................................................................75

Appendix A..................................................................................................................................78

Appendix B..................................................................................................................................80

v 
Figures 

Figure 1. An example of the section view....................................................................................4

Figure 20. The left panel shows an input waveform (red) and an output waveform (blue). The 

Figure 21. Section view illustrating subsurface wave fronts and scan plane beneath a GPR 

Figure 2. An example of section view..........................................................................................5

Figure 3. The Grid Parameters Dialog box................................................................................12

Figure 4. The Filter Dialog box. .................................................................................................16

Figure 5. The Migration Dialog box. ..........................................................................................19

Figure 6. The UXO Forward Model Dialog box..........................................................................21

Figure 7. Simulated ALLTEM data for a horizontal rod-like UXO target. ...................................22

Figure 8. UXO Inversion Calibration Dialog box........................................................................24

Figure 9. The UXO Inversion Dialog box...................................................................................27

Figure 10. The progress dialog box shown during the UXO inversion process. ........................28

Figure 11. Survey data from the Yuma Proving Ground............................................................29

Figure 12. ALLTEM time domain waveform. .............................................................................30

Figure 13. The VDL Colors Dialog box......................................................................................32

Figure 14. The GP Workbench toolbar......................................................................................34

Figure 15. The GP Wave Utilities view showing two waveforms. ..............................................38

Figure 16. GP Wave Utilities view of wave field data. ...............................................................39

Figure 17. The Plot Parameters Dialog box...............................................................................40

Figure 18. Text File Import Parameters Dialog box. ..................................................................42

Figure 19. The Filter Dialog box. ...............................................................................................43

right panel shows the calculated response function for the signal-processing module. ..............49

antenna (adapted from Oden, 2006). ..........................................................................................51

Figure 22. The Depth Shift Dialog box. .....................................................................................52

Figure 23. Diagram of a GPR system response over a perfect horizontal planar reflector........54

Figure 24. Flow chart for deconvolving for reflector properties..................................................54

Figure 25. The Invert for Ground Properties Dialog box............................................................56

Figure 26. The model space grid...............................................................................................59

Figure 27. The Model Ground Properties Dialog box. ...............................................................60

Figure 28. The Compare Waveforms Dialog box. .....................................................................64

Figure 29. The GP Wave Utilities toolbar. .................................................................................67

Tables 

Table 1. A listing of file types that can be imported with GP Workbench. ....................................6

Table 2. Gridding methods.........................................................................................................13

Table 3. Grid file (*.xyz) entries..................................................................................................15

Table 4. File types imported by GP Wave Utilities. ....................................................................41

Table 5. Methods of extracting waveform attribute sets.............................................................62

vi 
GP Workbench Manual: Technical Manual, 

User’s Guide, and Software Guide 

By Charles P. Oden and Craig W. Moulton 

Introduction 

GP Workbench is a general-purpose geophysical data processing program written primarily for 

ground penetrating radar (GPR) data, and also includes support for several USGS prototype 

electromagnetic instruments. The two main programs in the package are GP Workbench and GP Wave 

Utilities. GP Workbench provides the standard two-dimensional (section view) and three-dimensional 

(plan view) processing for GPR data, as well as some research features. GP Workbench can produce 

high-quality graphics for reports when Surfer 8 or higher (Golden Software, Inc., Golden, CO, not 

included) is installed. GP Wave Utilities provides a wide range of processing algorithms for single 

waveforms. GP Wave Utilities is used primarily for calibrating radar systems and processing individual 

traces. Both programs are written to run on the Windows TM operating systems. This release of GP 

Workbench contains an install program called InstallGPWorkbench1_7_1_11.exe. Run this installation 

program to install GP Workbench and associated programs to a Windows operating system. The source 

code for the project can be found in the self-extracting zip-file GPWorkbench1_7_1_11.exe. See 

Appendix B for more information about compiling GP Workbench. 

GP Workbench can import GPR data file formats used by major commercial instrument 

manufacturers including Sensors and Software, GSSI, and Mala. GP Workbench also supports data file 

1
formats for some USGS systems such as the VETEM system (Wright and Chew, 2000), the ALLTEM 

system (Wright et al., 2006), and the RTDGPR system (Wright et al., 2005). The native file format is 

SU (Seismic Unix, Cohen and Stockwell, 2003). Files generated by GP Workbench can be read by SU 

as well as many other data processing packages. 

The intent of releasing this software is not to duplicate the commercial GPR processing 

packages, but to provide a framework for researchers to test and implement new processing algorithms. 

These programs provide a toolbox that is applicable to a wide range of equipment and problems. By 

adding to the source code provided here, scientists can investigate new processing algorithms. Both GP 

Workbench and GP Wave Utilities were written with Microsoft Visual C++ 6.0. The source code files 

extracted from GPWorkbench1_6_5_22.exe are organized into two main project directories (one for 

each program) in the standard manner for Visual C++ projects. An easily distributable installation 

program is generated using the freely available Inno Setup Compiler (Russel, 2004). A few functions 

rely on programs and code that were written by others. A GPR forward modeling program (Powers, 

1995) was modified and included. Some data structures and a program to summarize information in 

GPR file headers were adapted from Lucius and Powers (2002). The prjpt.exe program (Webring, 

1998) facilitates conversions from latitude and longitude to northing and easting (or vise versa). Some 

wavelet transform functions were adapted from Press et al. (1992). The MatClass library (Birchenhall, 

1993) is used for matrix calculations. GP Workbench uses a Windows TM automation interface to 

Golden Software’s Surfer 8 package (not included) to generate high-quality plan views of time-slice 

data. Finally, note that the structural organization of the software is not optimal, and that individual 

routines have not been optimized for speed. Although the software is object oriented, strict adherence 

to object-oriented programming methodologies has not been followed. The only conventions followed 

are to write in an active voice and to use ‘Hungarian notation’ in naming variables. 

2
GP Workbench uses a right-handed Cartesian coordinate system. In the instrument reference 

frame, the x-axis is directed to the right of the survey line, the y-axis is directed along the survey line, 

and the z-axis is directed upwards. In the map reference frame, the x-axis is in the direction of easting, 

the y-axis is in the direction of northing, and the z-axis is directed upwards. When the bearing of the 

survey instrument is zero degrees, the instrument reference frame and the map reference frame are 

identical. Unless otherwise noted, frequencies are in megahertz (MHz), time is in nanoseconds (ns), and 

position is in meters (m). 

Disclaimer 

This open-file report was prepared by an agency of the United States Government. Neither the 

United States Government nor any agency thereof nor any of their employees makes any warranty, 

expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or 

usefulness of any information, apparatus, product, or process disclosed in this report or represents that 

its use would not infringe privately owned rights. Reference therein to any specific commercial product, 

process, or service by trade name, trademark, manufacturer, or otherwise does not constitute or imply its 

endorsement, recommendation, or favoring by the United States Government or any agency thereof. 

Although all data and software in this open-file report have been used by the USGS, no 

warranty, expressed or implied, is made by the USGS as to the accuracy of the data and related 

materials and (or) the functioning of the software. The act of distribution shall not constitute any such 

warranty, and no responsibility is assumed by the USGS in the use of these data, software, or related 

materials. These programs are distributed "as is". There are certain to be "quirks" and "bugs" in some of 

the software. Registered trade names such as Microsoft, Windows, Surfer, etc., are used for reference in 

context and do not necessarily imply endorsement by the U.S. Geological Survey. 

3
GP Workbench Views 

An entire survey can be loaded and examined in section view or plan view. The section view 

provides hyperbola fitting for velocity analysis, zooming in time or position, and panning. Waveforms 

are displayed at mouse click position as well as time and depth to reflector. Figure 1 shows an example 

of the section view. The small vertical red lines at the top of the section indicate marked traces (i.e., 

traces with a known position). The cursor index values shown in the legend are the trace number and 

sample number at the current cross hair position. 

Figure 1. An example of the section view. 

4
An example of the plan view is shown in figure 2. The legend shows the time of the current 

time slice and the position at the cross hairs. The cross hairs are placed by clicking on the image with 

the mouse. Several plan view imaging methods are available with varying quality and speed. The 

ungridded method is the fastest but provides a low-quality image. The gridded method displays an 

image that is the result of a conical filter applied to the data (see Grid Parameters in the Menu and 

Toolbar Functions Section). The final method uses Surfer (Golden Software) to produce a high quality 

image. 

Figure 2. An example of section view. 

5
Menu and Toolbar Functions 

This section provides an overview of the GP Workbench menu commands and toolbar buttons. 

When the mouse pointer is placed over a menu command or a toolbar button, a help message appears in 

the status bar at the bottom of the screen that describes the function of the menu item or button. 

Additionally, a flying help message appears when the mouse pointer is placed over a toolbar button. 

File Menu 

The File menu contains the standard menu items offered by most Windows applications. A new 

document can be created or saved, and recent files can be opened. The page for printing can be set up, 

documents can be printed, and the printout can be previewed. Non-standard Windows menu items are 

described below. 

Import 

The Import and Export menu items offer additional functions. The Import, Survey Line (section 

view) and Import, 3D Survey (plan view) menu items allow the user to select and import data files with 

various formats as listed in table 1. 

Table 1. A listing of file types that can be imported with GP Workbench. 

File Type (extension) Format Description 

*.su, *.txt Generic data formats 

*.dzt GSSI 

*.dt1 Sensors and Software 

*.rd3 Mala 

*.fna VETEM (USGS) 

*.aq1, *.aq2 RTDGPR (USGS) 

*.uxa, *.uxb, *.h1m, *.zzm ALLTEM (USGS) 

6
Data are saved in the *.SU format (see Seismic Unix, Cohen and Stockwell, 2003). Note that the trace 

identification code used in GP Workbench *.SU files ranges from 200 to 202, whereas most seismic 

programs use an identification code of 1. This difference may cause some difficulties when opening GP 

Workbench *.SU files with other programs. Also, note that *.SGY files, which are an industry standard 

in the seismic industry, have the same format as *.SU files with the exception of an added header at the 

beginning of the file. The Survey Line (section view) menu option imports data and displays it in 

section view, while the 3D Survey (plan view) option imports and displays data in plan view. In plan 

view, a data value recorded at a given time is displayed for each trace as a function of lateral position. 

The display time of the view can be adjusted using the up and down arrow keys, so that the user can 

observe the data in plan view while moving forward or backward in time. 

The position of each trace recorded by a prototype USGS instrument has an associated Global 

Positioning System (GPS) record. The GPS records contain the latitude and longitude of the instrument 

when each trace was recorded. The latitude and longitude of each trace is converted to northing and 

easting during the import process. When importing files containing GPS data, the import routines look 

for a file called Import.txt in the same directory as the data files being imported. This file contains 

information specifying how to convert from latitude and longitude to northing and easting. The 

specification of this conversion is location specific. If this file is not found, the assumed location is the 

Denver Federal Center in Lakewood, Colorado. A simple linear relation is used for the conversion. 

This linear relationship is usually sufficient when the lateral survey dimensions are less than a 

kilometer. The contents of the Import.txt file are shown below. 

Listing 1. Import.txt. 

GPRWorkbench Import File v. 1.0 


LocationName=YPG 


LatitudeReference=32.8531143 


LongitudeReference=114.2466241 


LatitudeFactor=108474.8000000 


LongitudeFactor=-96491.387000 


7
NorthingOffset=13.681 


EastingOffset=15.363 


The following relationships are used in the conversion: 

Northing = (Latitude − LatitudeReference)⋅ LatitudeFactor − NorthingOffset (1) 

Easting = (Longitude − LongitudeReference)⋅ LongitudeFactor − EastingOffset (2) 

The program prjpt.exe (Webring, 1998) can be used to convert latitude and longitude to UTM (universal 

transverse mercator) coordinates. 

Multiple files can be selected for importing. When selecting multiple files, note the order that 

the files are listed in the File Name edit box of the Open dialog box. The files will be imported and 

concatenated in the order they appear in this edit box. File concatenation is helpful for some data 

processing procedures, especially when combined with a filter operation (filters are described below). 

For example, the response to local geology that has been recorded in one file can be subtracted from the 

response of a target buried in similar geology that has been recorded in another file. 

Export 

The Export, Section Data menu item exports section view data to a *.csv text file that can be 

imported by other applications for further processing. The Export, Time Slice XYZ menu item exports 

data at the currently displayed time slice to a *.csv text file that can be imported by other applications 

for further processing. 

ALLTEM Data 

When working with ALLTEM data, some File Menu commands have enhanced functionality. 

Surveys with the ALLTEM system typically produce multiple data files with different combinations of 

coil polarizations. GP Workbench can open, close, save, and process these files simultaneously or 

8
individually. Files with extensions *.h1m or *.zzm are considered master ALLTEM files, and the 

others are called associated ALLTEM files. To work with these files simultaneously, open, close, save, 

or process a master file; GP Workbench will ask if you want to make the same operation on the 

associated files. The files associated with a given ALLTEM master file are defined in the file 

ALLTEM.def that is located in the same directory as the master file. A sample listing of this file is 

given below. Note that ALLTEM.def is a standard Windows section file, and therefore comments can 

be included as long as the first character in a comment line is a semicolon. 

Listing 2. ALLTEM.def. 

[ALLTEM_System_Definition] 


Number_of_polarization_combinations=7 


Master_file_extension=h1m 


StaticSample1=313 


StaticSample2=868 


PeakSample1=337 


PeakSample2=892 


[Combination_1] 


File_extension=h1m 


TX_polarity=z 


RX_polarity=z 


Cart_center_height=0.678 


RX_midpoint_position_x=0.0 


RX_midpoint_position_y=0.0 


RX_midpoint_position_z=0.0 


RX_offset_x=0.0 


RX_offset_y=0.0 


RX_offset_z=0.5 


[Combination_2] 


File_extension=hbl 


TX_polarity=z 


RX_polarity=z 


Cart_center_height=0.678 


RX_midpoint_position_x=-0.25 


RX_midpoint_position_y=-0.25 


RX_midpoint_position_z=0.0 


RX_offset_x=0.0 


RX_offset_y=0.0 


RX_offset_z=0.5 


[Combination_3] 


File_extension=hbr 


TX_polarity=z 


9
RX_polarity=z 


Cart_center_height=0.678 


RX_midpoint_position_x=0.25 


RX_midpoint_position_y=-0.25 


RX_midpoint_position_z=0.0 


RX_offset_x=0.0 


RX_offset_y=0.0 


RX_offset_z=0.5 


[Combination_4] 


File_extension=hfl 


TX_polarity=z 


RX_polarity=z 


Cart_center_height=0.678 


RX_midpoint_position_x=-0.25 


RX_midpoint_position_y=0.25 


RX_midpoint_position_z=0.0 


RX_offset_x=0.0 


RX_offset_y=0.0 


RX_offset_z=0.5 


[Combination_5] 


File_extension=hfr 


TX_polarity=z 


RX_polarity=z 


Cart_center_height=0.678 


RX_midpoint_position_x=0.25 


RX_midpoint_position_y=0.25 


RX_midpoint_position_z=0.0 


RX_offset_x=0.0 


RX_offset_y=0.0 


RX_offset_z=0.5 


[Combination_6] 


File_extension=vlr 


TX_polarity=y 


RX_polarity=y 


Cart_center_height=0.678 


RX_midpoint_position_x=0.0 


RX_midpoint_position_y=0.0 


RX_midpoint_position_z=0.28 


RX_offset_x=0.0 


RX_offset_y=0.5 


RX_offset_z=0.0 


[Combination_7] 


File_extension=vfb 


TX_polarity=x 


RX_polarity=x 


Cart_center_height=0.678 


RX_midpoint_position_x=0.0 


RX_midpoint_position_y=0.0 


RX_midpoint_position_z=0.28 


RX_offset_x=0.5 


RX_offset_y=0.0 


RX_offset_z=0.0 


 10 
This listing also shows the default settings that are used when this file does not exist in the same 

directory as the master file. The definitions in this file are used to calculate the coil locations and 

polarizations that were used when recording data files with the indicated extension. The location of 

each coil is specified with respect to the X-Y location (easting and northing respectively) of the center 

of the ALLTEM cube (Z is positive up); the transmitter coils are 1-meter-square loops centered about 

the center of the cube with the indicated polarization (TX_polarity). The center of the differential 

receiver coil pairs are specified by the RX_midpoint_position entry, and the location of each individual 

coil is found by adding or subtracting the RX_offset value from the midpoint of the coil pair. 

Process Menu 

The Process menu contains a variety of commands, which are used to make calculations on 

existing data and to generate new data. 

Grid Parameters 

The first Process Menu item, Grid Parameters, controls the processing steps used in processing 

the position information associated with each trace. In many surveys, the position of the instrument is 

not measured for each trace that is recorded. Rather, the positions of a few traces are recorded, and the 

remaining trace positions are determined using one of the methods described in table 2. This process is 

referred to as gridding. The Grid Parameters Dialog box shown in figure 3 controls how the positions 

of all traces are calculated. This dialog box writes an *.xyz file (same file name as imported data file) 

that describes how the trace positions are to be calculated when the grid button (see toolbar) is pressed. 

The *.xyz file is a text file that typically contains the positions of each marked trace or of each trace. 

Since this file is a text file, it can be edited as needed. When viewing imported or un-gridded data, the 

measured locations (or marked positions) are displayed as red marks at the top of each section view. 

11
Figure 3. The Grid Parameters Dialog box. 

After setting the gridding parameters with the above dialog box, press the grid button on the toolbar 

to grid the document. 

12
Table 2. Gridding methods. 

Method

Use marks on screen or in data file 

 Description 

An entry in the *.xyz file is made corresponding to each mark in the 

data file (and on screen). The marked XYZ positions are specified by 

the XYZ locations part of the dialog. The location of unmarked 

traces is determined by linear interpolation between marked traces. 

Mark each trace An entry in the *.xyz file is made corresponding to each trace in the 

data file (and on screen). The marked XYZ positions are specified by 

the XYZ locations part of the dialog. 

Use trace header GPS data and interpolate The location of unmarked traces is determined by linear interpolation 

between traces with known GPS positions. 

Use trace header position in data The location of each mark or trace is determined by the information 

in each trace header. Some data file formats provide a position for 

each trace, or for marked traces. If unmarked traces exist, the 

location of unmarked traces is determined by linear interpolation 

between marked traces. 

One fixed position Use the same position for all traces in the data file. 

In many cases, geophysical data are not sampled on a uniform grid. Some processing algorithms 

such as fast Fourier transforms require data that are sampled on a uniform grid. Check the resample on 

a uniform grid box to resample the data on a uniform grid. After an XYZ position has been calculated 

for each trace using one of the methods described above, a new trace for each uniformly spaced grid 

point is calculated by interpolating between the recorded traces at their recorded positions. This is a 

simple amplitude interpolation made independently at each time sample. This method of resampling 

can cause distortion in the waveforms. Interpolating the waveforms in the frequency domain yields 

better results. The frequency-domain interpolation will be implemented in a future release. The grid 

interval specifies the interval (in meters) between adjacent resampled traces. The bin width specifies the 

width (in meters) of a moving window that is used to calculate the re-sampled trace. For each uniformly 

spaced grid point, the window is centered on the grid point, and all original traces in the window are 

averaged with equal weighting to obtain the resampled trace. 

13
When gridding a master ALLTEM document (has an *.h1m or *.zzm extension), GP Workbench 

asks the user if the associated ALLTEM documents should also be gridded (see the File Menu Section 

for a description of associated ALLTEM documents). Note that a single grid file (*.xyz) is used for the 

master ALLTEM document and the associated documents. To use different gridding parameters with 

each ALLTEM document, each document must be gridded separately. In this case, select the document 

to be gridded (Window menu), make a grid file (*.xyz) for the current document (Process, Grid 

Parameters), and then press the grid button . Repeat for each document. 

The grid files (*.xyz) are ASCII text files that can be edited with a standard text file editor such 

as Notepad. A typical grid file is shown below. The first line of the file must be exactly as shown. 

Listing 3. Grid Parameters (*.xyz) file. 

GPR Workbench Mark File v.1.0 


G 


0.000 0.000 0.000 


The second line describes the data to follow. Valid entries on the second line are described in the table 

3. The lines following line two contain the X-Y-Z coordinates for each trace or marked trace. 

When gridding a plan view of a time slice, the position of unmarked traces is determined by 

interpolating between the known positions of marked traces. The plan view gridded image is result of a 

conical filter applied to all of the data in screen pixel space. The radius of the conical filter is specified 

with the Grid Parameters Dialog box. This radius is specified in virtual screen pixels, where the image 

is plotted on an area of 663 by 373 virtual screen pixels. The value of each screen pixel is the weighted 

average of all traces at the current time-slice value that fall within the radius of the filter’s conical 

window. The weighted value of each trace depends on where it lies under the conical window. The 

conical window has a maximum value of one at the center, and zero at the outside radius. Note that this 

conical filter is only used to display the data, and that no changes are made to the data with this gridding 

operation. 

14
Table 3. Grid file (*.xyz) entries. 

Line Two Entry Description 

M File contains the position of each marked trace. The position of each marked trace is listed on the 

following lines. 

MR As above, and the data will be resampled on a uniformly spaced grid. 

T File contains the positions for each trace. The position of each trace is listed on the following lines. 

TR As above, and the data will be resampled on a uniformly spaced grid. 

G Use the GPR positions in the data file for each marked trace. When importing data with GPR 

positions, a trace is marked when its position differs from that of the previous trace. Lines following 

line two are ignored. Note that the entries in the file Import.txt are used to convert from latitude and 

longitude to northing and easting. This option is used primarily for USGS prototype instruments. 

GR As above, and the data will be resampled on a uniformly spaced grid. 

H The position information in the file’s trace header is used to determine the position of each trace. 

Lines following line two are ignored. 

HR As above, and the data will be resampled on a uniformly spaced grid. 

F Line three contains a position that will be used for all traces. 

FH As above, and the data will be resampled on a uniformly spaced grid. 

Position Shift 

The Process, Position Shift menu item facilitates a positional shift of all traces in the current 

document. This operation is useful, for example, to account for an offset between the GPS position 

sensor and the instrument’s measurement point. After selecting this menu item, the user is presented 

with a dialog box in which the shifts along the x-axis (easting) and y-axis (northing) can be entered. 

When applying a position shift to a master ALLTEM document (has an *.h1m or *.zzm extension), GP 

Workbench asks the user if the associated ALLTEM documents should also be shifted (see the File 

Menu Section for a description of associated ALLTEM documents). 

Filter 

The Process, Filter menu item opens the dialog box shown in figure 4. The Filter Dialog box 

can be used to filter section view documents or plan view documents. Seven filter types are offered – 

15
three spatial filters, three temporal filters, and trace filtering. All filters operate on all traces in a 

document. 

To remove the entire spatial background, GP Workbench calculates the background value of all 

traces in the document and subtracts the background trace from each trace. The background trace can 

be the average value of all traces, the median value of all traces, or the RMS value of all the traces. The 

average, median, or RMS values are calculated for each sample time included in the traces. Spatial 

filters are a function of the trace number. The trace number (beginning with zero) is assigned 

Figure 4. The Filter Dialog box. 

16
sequentially to each trace when a data file is opened or imported. When removing a spatial window, a 

background trace is subtracted from each trace. The traces used to calculate the background trace can 

be taken from a moving window (specify filter width in traces) or a fixed window (check fixed window 

box and enter beginning and ending trace number). Alternatively, the trace number can be used to 

calculate an interpolated background trace between two fixed windows. The midpoint trace number of 

each fixed window is used as the endpoints in the interpolation. 

A set of temporal filters is also available. Specification of these filters is analogous to the spatial 

filters, except that the temporal filters are functions of time or sample number. Additionally, frequency 

filters are available. The frequency filters are 100-point FIR (finite impulse response) filters. The high 

pass and low pass cutoff frequencies can be specified to realize high pass, low pass, or band pass filters. 

No high (low) pass filter operation is performed if the high (low) pass cutoff frequency is zero. 

The final filter type can be used to remove traces from a document. This is useful, for instance, 

when isolating data containing the response to a single buried target to invert for target attributes. 

Specify the beginning and ending traces of a window, and delete the traces inside or outside this 

window. 

Variogram 

The next menu item, Process, Variogram, can be used to generate a horizontal semi-variogram 

image of the data. See Olhoeft (1994) for details on how the semi-variogram is calculated. The user 

will be prompted for the minimum and maximum range values to use in calculating the horizontal 

variogram. The document must be gridded to enable this menu item. 

17
View Header 

The Process, View Header (GPR_RHDR) command invokes the GPR_RHDR.exe utility (Lucius 

and Powers, 2002). This utility can be used to view the contents of GPR data files from most 

commercial manufacturers. File headers commonly have valuable information such as spatial sample 

interval, temporal sample interval, number of samples, and notes about the survey site. GPR_RHDR.exe 

does not return information about USGS prototype equipment. 

GPR Forward Model 

The Process, GPR Forward Model command invokes the 2.5 dimensional GPR modeling 

program called 2DRadar.exe (Powers, 1995). After choosing a *.2DP input parameter file for 

2DRadar, a status screen displays progress as the program executes. The *.2DP file specifies the model 

geometry and material properties of the subsurface. The simulation generates a synthetic dataset that is 

displayed in section view. 2DRadar can be used to simulate GPR data profiles for user-specified 

scenarios such as buried pipes in various media types, sand lenses, etc. After the simulation has 

completed, press the model view and section view buttons on the tool bar to change between 

views of the model used to generate the section data, and the section data. Consult Powers (1995) for 

more information on the *.2DP parameter files and the simulation method. Some example *.2DP files 

can be found in the directory where GP Workbench has been installed. 

Migrate 

The Process, Migrate command allows one of several frequency-domain migration techniques 

to be applied to the data. Note that the data must be gridded before migration is allowed. All of the 

migration techniques are based on the frequency-wavenumber method (Gazdag, 1978). The material 

properties of the background medium must be known for frequency-wavenumber migration. The 

18
dielectric permittivity and magnetic permeability are frequency dependent quantities, and are specified 

in terms of the Cole-Cole distribution (Olhoeft, 1985). Four variations on the frequency-wavenumber 

method are available (see fig. 5). The first method uses fixed values for the dielectric permittivity and 

magnetic permeability. The fixed values are taken as the DC or zero frequency values from the ColeCole distribution. According to the Kramers-Kronig relations (Jackson, 1999), if a medium has a 

frequency-dependent propagation velocity, then it must also have frequency-dependent attenuation. The 

second migration method accounts for the frequency-dependent velocity, but not the frequencydependent attenuation. The third migration method attempts to replace lost signal amplitude due to 

Figure 5. The Migration Dialog box. 

19
propagation in lossy media (Powers and Oden, 2004). Reversing attenuation during migration is an 

unstable process, and a gain cutoff must be used to prevent over-amplification of noise. Entering a 

smaller number for the gain cutoff reduces the amount of reverse attenuation that is applied. Using a 

gain cutoff that is too small results in no reversal of attenuation, while using a cutoff value that is too 

large results in an image dominated by noise. 

The final migration method uses knowledge of the source spectrum to limit the gain during the 

reverse attenuation operation (Oden, 2006; Oden et al., 2006). This gain-limiting scheme produces a 

higher resolution migrated image when working in lossy media. This scheme also prevents the overamplification of noise that would otherwise degrade the migrated image quality. To use this migration 

scheme, a document must be open that contains the system response waveform. The system response 

waveform is defined as the response of the system to an infinite horizontal conductive sheet that is 

buried in the ground just deep enough that it does not affect the current distribution on the GPR 

antennas. The system response can be determined through calibration, simulations, or experimental 

measurements (see Oden, 2006, for more details). During migration, the gain is limited by a frequencydomain limiting spectrum. This limiting spectrum is usually the system response (when the source 

interpolation is one). Choosing a source interpolation value closer to zero adjusts the limiting spectrum 

to be more like the received spectrum (i.e., the recorded data). Choosing a source spectrum cutoff factor 

greater than one (one is the usual setting) allows more amplification before gain limiting begins. This 

limited-gain migration algorithm is described in more detail in Oden (2006). 

UXO Forward Model 

The Process, UXO Forward Model menu item generates an ALLTEM synthetic dataset for a 

UXO (unexploded ordinance) target. The induced magneto-static fields in the UXO target can be 

described by three orthogonal dipoles. Similarly, the magnetic fields at a given moment in time 

20
Figure 6. The UXO Forward Model Dialog box. 

produced by the induced electro-dynamic response can be described by three orthogonal dipoles. The 

UXO Forward Model Dialog box is shown in figure 6. For the induced magneto-static response, the 

strengths of the three principal polarizability components are specified. The attitude of the first 

principal polarizability ( mˆ 1 ) is described in terms of azimuth and inclination from horizontal, the 

principal polarizability is horizontal (mˆ 2 = mˆ 1 × zˆ ), and the attitude of the third principal polarizability 

is the cross product of the first two ( mˆ 3 = mˆ 1 × mˆ 2 ). The location of the UXO can be specified, as well 

as the survey parameters. The survey parameters describe the spatial data density of the simulated data. 

21
2 2 2 2 

1.5 1.5 1.5 1.5 

1 1 1 1 

0.5 0.5 0.5 0.5 

0 0 0 0 

0 0.5 1 1.5 2 0 0.5 1 1.5 2 0 0.5 1 1.5 2 0 0.5 1 1.5 2 

2 2 2 

1.5 1.5 1.5 

1 1 1 

0.5 0.5 0.5 

0 0 0 

0 0.5 1 1.5 2 0 0.5 1 1.5 2 0 0.5 1 1.5 2 

2 2 2 2 

1.5 1.5 1.5 1.5 

1 1 1 1 

0.5 0.5 0.5 0.5 

0 0 0 0 

0 0.5 1 1.5 2 0 0.5 1 1.5 2 0 0.5 1 1.5 2 0 0.5 1 1.5 2 

2 2 2 

1.5 1.5 1.5 

1 1 1 

0.5 0.5 0.5 

0 0 0 

0 0.5 1 1.5 2 0 0.5 1 1.5 2 0 0.5 1 1.5 2 

Figure 7. Simulated ALLTEM data for a horizontal rod-like UXO target. The first seven plots 

show the magneto-static response, and the last seven show the electro-dynamic response. 

These figures show the amplitude response of the ALLTEM induction system as a function of 

the X (horizontal axis) and Y (vertical axis) position of the ALLTEM system over a stationary 

target. Each of the seven images is the response for a different coil set (see text). 

22
Typically, these parameters would be used to specify the line spacing and the distance between data 

points on a line. It is assumed that the survey starts at (x, y) = (0.0, 0.0). Next, a log file must be 

specified. The simulation results will be written to the log file in text file format. Finally, an output 

format needs to be specified. Plots (see fig. 7) can be generated using Golden Software’s Surfer 

graphing package (Surfer version 8.0 or greater must be installed). Alternatively, the user can generate 

GP Workbench documents in plan view or section view. Synthetic datasets can be generated in this 

manner for testing the UXO inversion algorithm (see the Estimate UXO Parameters section below). 

Since the ALLTEM system can collect data using many coil combinations, the forward operator looks 

for the file ALLTEM.def in the same directory as the log file. The ALLTEM.def file specifies which 

coil combinations to use in the simulations. A Surfer plot or a GP Workbench document is created for 

each coil combination listed in the ALLTEM.def file. 

The ALLTEM UXO forward operator approximates the induced field response of a subsurface 

UXO. This operator describes the induced magnetic fields in the UXO in terms of three orthogonal 

principal polarizabilities. A set of principal polarizabilities is used to describe the induced magnetostatic response, and another set is used to describe the induced electro-dynamic response. The forward 

operator A used in the inversion has the form 

A([rcart ,rtx ,rrx1,rrx2Ptx ,Prx ,t],[rs,uxo ,φs ,θ s ,ms,1,ms,2 ,ms,3,rt,uxo ,φt ,θt ,mt ,1,mt ,2 ,mt,3,τ 1,τ 2 ,τ 3 ]) (3) 

= y([rcart ,rtx ,rrx1,rrx2 ,Ptx ,Prx ,t]) 

where rcart is the location of the center of the ALLTEM cart, rtx, rrx1, and rrx2 are the locations of the 

transmitting and receiving loops, Ptx and Prx are the polarizations of the transmitter and receiver coils, t 

is time, and y are the simulated data. The UXO parameter set is listed in the second set of square 

brackets in the argument list of the forward operator, where rs,uxo is the location of the UXO, ms,1 

through ms,3 are the magnitude of the three orthogonal induced magneto-static principal 

polarizabilities,φs , andθ s are the azimuth and inclination of the m1 component, rt,uxo, mt,1 through mt,3, 

23
φt andθt are the analogous parameters for the transient electro-dynamic response, and τ1 through τ3 are 

the late time decay constants for the transient dipole moments. A separate set of UXO parameters are 

used for the static and transient response to provide more information to better facilitate discrimination. 

This parameter set will be estimated by the inverse algorithm (see the Inversion Algorithm section 

below). The ALLTEM UXO forward operator is described in more detail in Appendix A. 

ALLTEM Calibration 

The Process, ALLTEM Calibration menu item provides a calibration routine for the ALLTEM 

system. After the data files containing the response to the calibration standard have been imported, and 

these data have been pre-processed, then the calibration routine can be used. Loading and preprocessing the data are described in the Data Processing Methods Section of this Manual. The preFigure 8. UXO Inversion Calibration Dialog box. 

24
------ ---------------- 

------ ---------------- 

processed data contain the response to the calibration standard after removing the response of the 

system and of the geologic background. Make sure that the master ALLTEM document is selected (see 

the Window menu), bring up the calibration dialog box, then enter the location of the center of the 

calibration standard and the center of the ALLTEM coil cube. The default values are for the cube in the 

cart (center is 67.8 cm above the ground where z = 0.0), and a calibration standard positioned on the 

calibration standard holder located directly above the cart center. The calibration standard should be a 

ferrous, isotropic, conductive sphere. Enter the relative size of the induced magneto-static and electrodynamic principal polarizabilities. The inversion algorithm will estimate the relative strengths of the 

induced dipole moments in the target with respect to the response of the calibration standard. The 

calibration routine will write a *.cal file with the same name and in the same directory as the master 

ALLTEM document. The *.cal file is a formatted text file that can contains two matrices (column 

vectors) – one for the magnetic-static response, and another for the electro-dynamic response (see listing 

below). 

Listing 4. A sample ALLTEM calibration (*.cal) file. 

# Kmu matrix: 7, 1 


Row | Col: 1 


1 | -26.0269 


2 | -11.8141 


3 | -10.893 


4 | -11.7343 


5 | -10.0501 


6 | 0.0 


7 | 0.0 


# Ksig matrix: 7, 1 


Row | Col: 1 


1 | 73.1902 


2 | 50.394 


3 | 53.4475 


4 | 46.9004 


5 | 54.6051 


6 | 0.0 


7 | 0.0 


25
This file can be viewed or edited with a standard text file editor; however, one should only change the 

floating point numbers (numbers with a decimal point) and ensure these numbers remain right justified 

as shown in the listing. The numbers in each row are the calibration coefficient for the corresponding 

coil combination listed in the ALLTEM.def file. 

Estimate UXO Parameters 

The Process, Estimate UXO Parameters menu item inverts the UXO forward model (see 

Appendix A) to estimate the parameters of a UXO target from the ALLTEM data acquired over the 

target. The documents must contain a pre-processed data set for a single target that has had the 

background response of the system and of the ground removed. The master ALLTEM document must 

be selected (see the Window menu). The user must specify whether to use the magneto-static response, 

or the electro-dynamic response in the inversion (see fig. 9). Usually, an inversion will be run for both 

types of response. One or more initial models must be selected (see discussion below). Enter estimates 

of the noise level in the data, and the bearing of the survey lines. A calibration file for the current 

survey data must be selected. The inversion procedure will attempt to estimate the UXO parameters 

listed in Equation 3. During the inversion, the progress dialog box is shown (see fig. 10). When the 

inversion is complete, the text in the dialog can be selected with the mouse, copied to the clipboard, and 

pasted into another application such as a text editor where the results are being saved. The progress 

dialog box displays the parameter estimates and the normalized mean squared error at the current 

iteration. The number after the MSE label is the normalized mean squared error (i.e., the RMS 

difference between the predicted data and the actual data normalized by the RMS value of the data). 

The forward operator generates a data set that can be compared with survey data to invert for 

UXO parameter estimates. For a given UXO target, a small data set containing the response to a single 

target is extracted from the field data set for use in the inversion. During a survey, the ALLTEM system 

26
records the response of various transmitting and receiving loop combinations, at successive positions 

along a survey line. A survey generally consists of a collection of survey lines. A small data set for an 

individual target is made by extracting received waveforms with locations inside a window over the 

target as shown in figure 11. For each waveform, two data values are extracted. One value is extracted 

for the static response (the mean value of the waveform values at static sample numbers specified in the 

ALLTEM.def file, see Listing 2), and another value is extracted for the transient response (see fig. 12). 

To obtain the transient response, the static response of each half-cycle (using the samples specified in 

the ALLTEM.def file) is subtracted from the recorded waveform. The transient value is taken at the 

sample time where the transient response is 25% of the peak transient value (the sample numbers for the 

peak value are specified in the ALLTEM.def file) so that late time decays are measured. When 

Figure 9. The UXO Inversion Dialog box. 

27
processing data from multiple ALLTEM coil combinations, the transient values for each document are 

taken at the same sample time. To find this sample time, the RMS transient response of all documents 

is used. The transient sample time is when the RMS value of the waveforms in all documents is 25% of 

the RMS peak in all documents. The resulting data set for a given target is made up of the two 

waveform values extracted from waveforms recorded at various ALLTEM cart locations using various 

Figure 10. The progress dialog box shown during the UXO inversion process. 

coil polarizations. For a standard survey, using a 50-cm line spacing, 20 cm between each record along 

each line, and seven different polarization combinations of the transmitting and receiving loops, there 

will be about 350 waveforms and 1,050 data points that can be used to estimate 19 parameters. 

During the early-time transient response, the induced dipole moments are oriented parallel to the 

primary magnetic field. The induced eddy currents will attenuate quickly if the cross-sectional area of 

the currents in the target is small. Conversely, the induced currents will decay slowly if the cross 

28
sectional area containing the currents is large. Therefore, the late-time transient data contain 

information about the orientation of the large area conductive planes in the target, and the targets aspect 

ratio (length to width ratio). For these reasons, late-time data are used in the forward model and 

inversion, and early-time data are not processed. 

Figure 11. Survey data from the Yuma Proving Ground. Data are for 1 meter square 

vertically polarized transmitter and receiver coils. A data set for an individual target is made by 

extracting received waveforms with locations inside a window over the target as shown. 

The inverse algorithm uses the Gauss-Newton method (Gill et al., 1986; Zhdanov, 2002) to 

reduce the RMS difference between the predicted data and the measured data. This RMS difference is 

referred to as the objective function. The predicted data are calculated using the forward operator with 

an initial set of UXO parameters, and the value of the objective function is calculated. The inversion 

29
algorithm then chooses a new set of UXO parameters (i.e., a new model) by moving in the downgradient direction of the objective function. The value of the objective function using the new 

parameter set will (usually) be less than the value using the initial parameter set. This process is 

repeated until the objective function is minimized. An acceptable solution (i.e., UXO parameter set) to 

ALLTEM Time Domain Waveform 

-15000 

-10000 

-5000 

0 

5000 

10000 

15000 

0 2000 4000 6000 8000 10000 12000 

Time (us) 

A

m

plit

u

d

e (

digitiz

er u

nits) 

Static response 

Transient response 

Peak response 

Figure 12. ALLTEM time domain waveform. One data point is extracted for the static 

response, and two data points are extracted for the transient response. 

the inverse problem is reached when the RMS difference between the predicted data and the actual data 

is less than the RMS noise level in the data. Since this is a non-linear problem, it is possible that the 

UXO parameter set determined by the inversion algorithm is located at a local minimum of the 

objective function rather than a global minimum. The parameter set at a local minimum may not be an 

acceptable solution to the inverse problem, and may have markedly different values than the parameter 

set at the global minimum. To solve this problem, the inversion can be performed with up to five 

different initial models (see fig. 9). An attempt to solve the inverse problem is made with each initial 

model, and the process stops when an acceptable solution is found. For each initial model, the inversion 

30
algorithm takes 100 steps (iterations) down gradient. If an acceptable solution is not found after 100 

iterations, the next initial model is started. While stepping down gradient, a reasonable estimate of the 

location and depth of the target is usually found in the early iterations. Therefore, the target location 

and depth do not change with subsequent initial models. For the first initial model, the target location is 

at the center of the data set, and the depth is 1 meter (z = -1). A similar inversion scheme has been used 

successfully for other non-linear problems (Oden, 2006; Oden et al. 2006). 

The noise level in the data can be estimated by dividing the amplitude of the signal variations 

measured over no target by the amplitude of the signal variations over the target. The RMS noise in the 

data has several components. They are system noise from the electronics and sensors, electromagnetic 

interference, the response from varying electromagnetic properties of the soil, the response to 

topography, and the response due to the changing attitude of the cart. The inversion procedure can more 

accurately estimate the UXO properties when the noise level is small. The uncertainty in a UXO 

parameter estimate can be calculated from the estimated value near the final solution according to 

pn+1 − pn ⋅σ RMS σ = , (4) p Σn+1 − Σn 2

where pn is the parameter estimate at the nth iteration, Σn is the normalized mean square error between 

the actual data and the predicted data at the nth iteration (obtained from the inversion progress dialog), 

σRMS is the RMS noise level in the data, and σp is the uncertainty in the estimate of parameter p. 

It should be noted that the inversion algorithm provides estimates of the UXO parameter set 

defined by the forward model. The inversion algorithm does not attempt to identify specific targets. A 

library should be compiled that maps specific UXO types to the UXO parameter set estimated by the 

inversion. The aspect ratio of the targets is determined from the three principal polarizabilities. If one 

component is significantly larger than the others, then the target is rod-like. If one component is 

significantly smaller than the others, then the target is disk-like. If all polarizabilities have similar 

31
magnitudes, then the target is ball-like. A separate location and attitude of the principal polarizabilities 

are determined for the magneto-static response and for the electro-dynamic response. Therefore, it may 

be possible to compare these locations and attitudes to better discriminate and identify targets when a 

priori information about the targets is known. 

View Menu 

The view menu contains commands for turning the toolbar and the status bar on or off. The nonstandard commands are described below. 

VDL Colors

 The View, VDL Colors menu item brings up the dialog shown in figure 13. This dialog can be 

used to edit the color map used to represent data values. The color map is defined by a collection of 

anchor color values. The colors that represent non-anchor values are determined by interpolating 

between anchor values. Anchor values can be added by typing in a new numerical value in the value 

list box and pressing the add/edit button to define a new color. Pressing the add/edit button invokes a 

color selection dialog box. The anchor value currently displayed by the list box can be deleted with the 

delete 

Figure 13. The VDL Colors Dialog box. 

32
current button. A set of anchor color values can be saved to a color map file (*.clr), or read from a 

color map file. 

Surfer 

Use this menu item to generate a rendered surface of the current time slice using Surfer (Golden 

Software). If the file crop.bln is found in the current document’s directory, then it is used to crop 

portions of the rendered surface. Usually, cropping is used to remove portions of the rendered surface 

that lie outside regions containing data. Consult the Surfer help file for more information on generating 

and using *.bln files. This menu item is only available in plan view. 

Make JPGs 

Use this menu item to generate a rendered surface for all time slices in the current document 

using Surfer (Golden software). Each rendered surface is saved in a *.jpg file. The entire set of *.jpg 

files can be used to make a movie, or viewed sequentially using a *.jpg viewer program (such as 

Windows Explorer). This menu item is only available in plan view. 

GP Wave Utilities 

If a trace is currently marked by the cross hairs in section view, then this menu item starts GP 

Wave Utilities, which opens with the currently marked trace. GP Workbench has many waveform 

utilities that may be helpful in processing and characterizing a given response. This menu item provides 

a method to open any data set that can be read by GP Workbench using GP Wave Utilities. 

Window and Help Menus 

The Window and Help menus contain standard windows menu items for selecting a view, 

arranging view windows, and obtaining information about GP Workbench and how to use it. 

33
The Toolbar 

The toolbar shown below provides many operations for viewing and processing data. Each 

button on the toolbar will be described in more detail below. 

Figure 14. The GP Workbench toolbar. 

The new document button generates a new empty document. A document must exist to be able 

to import data files. 

The Open-File dialog can be used to open GP Workbench document files (*.SU files). 

The save document saves the current document. If a file has not yet been associated with the 

document, the user is prompted for a file name. 

The print button prints the current document using the default printer. 

This button displays information about GP Workbench such as the version number. 

This button opens the help file (GP Workbench.pdf) that describes the various operations of GP 

Workbench. 

This button increases the display contrast. The current contrast setting is shown in the status bar 

(Disp Gain: 1.00). 

This button decreases the display contrast. The current contrast setting is shown in the status bar 

(Disp Gain: 1.00). 

This button turns the range gain used in the acquisition of the data on or off. 

This button increases the data multiplier that is applied to the data before plotting. The current

multiplier setting is shown in the status bar (Gain (A): 1.00). 

34
This button decreases the data multiplier that is applied to the data before plotting. The current 

multiplier setting is shown in the status bar (Gain (A): 1.00). 

This button increases the exponential gain that is applied to the data before plotting. The current 

gain setting is shown in the status bar (Gain (x): 1.00). 

This button decreases the exponential gain that is applied to the data before plotting. The current 

gain setting is shown in the status bar (Gain (x): 1.00). 

This button is used to set the zero time sample number. Many GPR systems begin recording 

GPR traces before the transmitter fires. On section view, click on the section image, then use the

arrow keys to move the cross hairs to the desired zero time sample. Press this button when the 

cross hairs are at the desired time sample. 

This button selects the x-coordinate (usually easting) of each trace position when plotting in 

section view. 

This button selects the y-coordinate (usually northing) of each trace position when plotting in 

section view. 

This button selects the radial position coordinate ( r =
 x

2 + y2 ) of each trace position when 

plotting in section view. 

This button toggles the mark indicator associated with each trace. It is assumed that the position 

of marked traces is known (see the Process, Grid Parameters menu item). Note that the screen 

resolution may not be sufficient to display all marks in sections with many traces that span a 

large distance. For these cases, the horizontal zoom ( button) should be used to select a 

horizontal plotting scale where individual marks are clearly visible. This prevents inadvertently 

toggling the marks of unintended traces. 

35
This button is used to enter the velocity hyperbola mode. The diffraction hyperbola of point-like 

scatterers can be used to estimate the velocity of the overlying material. Velocity hyperbola 

mode is only available in section view mode with data that have been gridded (see the Process, 

Grid Parameters menu item and the button). For valid results, the zero time must also be 

selected before using the velocity hyperbola mode (see button). To estimate the velocity of 

the medium, place the cross hairs at the apex of a diffraction hyperbola. Then press the 

button or ‘v’ key to enter the velocity hyperbola mode. Press the arrow keys until the model 

hyperbola fits the diffraction hyperbola in the data. Then press the button or ‘v’ key to exit 

the velocity hyperbola mode. After fitting the velocity, click on the section or use the arrow 

keys to position the cross hairs at another position. The status bar displays the index of 

refraction squared (Nsq: 10.0), the velocity in m/ns (Vel: 0.1), and depth of the target (Depth: 

6.3). 

This button increases the horizontal zoom (section view only). 

This button decreases the horizontal zoom (section view only). 

This pans to the right (section view only). 

This pans to the left (section view only). 

This button increases the vertical zoom (section view only). 

This button decreases the vertical zoom (section view only). 

This pans up (section view only). 

This pans down (section view only). 

36
This can be used to delete a trace at the cross hairs in section view or plan view. Noisy traces 

can be deleted in section view. Traces with bad position records can be eliminated in plan view. 

If the current document is a master ALLTEM document (*.h1m or *.zzm files), the user will be 

prompted to delete the same trace number on the associated documents that are currently open. 

This button toggles between gridded and ungridded data. See the Process, Grid Parameters

menu item for more details on the gridding algorithm. 

This button displays the model geometry used to generate synthetic GPR data. See the Process, 

GPR Forward Model menu item for more details on the forward model. 

This button performs and displays the wavelet (Debauchies) decomposition of the waveform at 

the current crosshairs. 

This button selects or indicates the section view. 

This button selects or indicates the plan view. 

37
GP Wave Utilities 

The GP Wave Utilities program has many processing algorithms designed for a single 

waveform. It also contains routines for working with synthetic GPR data generated using FDTD 

simulation programs such as GPRMax (Giannopoulos, 2004) and XFDTD (Remcom, Inc., State 

College, PA). Although some of the processing algorithms are specific to the RTDGPR (USGS) radar 

system, they can be modified for use with other systems or for more general use. 

The GP Wave Utilities View 

Two data types can be viewed with GP Wave Utilities. These are one-dimensional waveform 

data and three-dimensional wave field data. Both of these data types can be viewed in the time or 

Figure 15. The GP Wave Utilities view showing two waveforms. 

38
frequency domains. Figure 15 shows a typical view where two waveforms are being compared. Figure 

16 shows a typical view of wave field data. In any view, click on the plot to position the cross hairs. 

The time and amplitude at the cross hair is shown on the status bar at the bottom of the screen. The xaxis is time in seconds (samples), and the y-axis is amplitude (for waveform data) or position (for wave 

field data). Right click on the plot to bring up the Plot Parameters Dialog box (see fig. 17). 

Figure 16. GP Wave Utilities view of wave field data. The vertical axis is position in meters 

(position scales are on right side of plot) 

39
Figure 17. The Plot Parameters Dialog box.

 The Plot Parameters Dialog box is used to change the view between the time-domain or 

frequency-domain views. The plotting scales and limits can be set for both time and frequency 

domains. Multiple waveforms can be plotted on the same graph by selecting auxiliary traces from other 

documents that are open in GP Wave Utilities. When plotting three-dimensional wave field data, the 

image plot consists of multiple traces (i.e., like raster lines on a television screen). The extents of a 

wave field image plot can be adjusted by selecting the range of traces to display and the minimum and 

maximum plotting time (or frequency). 

Menu and Toolbar Functions 

This section provides an overview of the GP Wave Utilities menu commands and toolbar 

buttons. When the mouse pointer is placed over a menu command or a toolbar button, a help message 

appears in the status bar at the bottom of the screen that describes the function of the menu item or 

40
button. Additionally, a flying help message appears when the mouse pointer is placed over a toolbar 

button. 

File Menu 

The File menu contains the standard menu items offered by most Windows applications. A new 

document can be created or saved, and recent files can be opened. The page for printing can be set up, 

documents can be printed, and the printout can be previewed. Non-standard Windows menu items are 

described below. 

Import and Export 

The Import and Export menu items offer additional functions. The Import menu items allow the 

user to select and import data files with various formats as listed below. 

Table 4. File types imported by GP Wave Utilities. 

File Type (extension) Format Description 

*.txt ASCII Text file 

*.aq1 RTDGPR (USGS) 

*.out GPRMax Field files (fields, current, etc.) 

*.src; *.txt GPRMax Source files 

*.g0 XFDTD Wave files 

*.fss XFDTD Field files 

*.src XFDTD Source files 

*.ao XFDTD Impedance files 

GP Workbench (see GP Workbench, View, GP Wave Any file imported by GP Workbench (see table 1) Utilities command) 

When importing text files, the dialog shown in Figure 18 is displayed. This dialog displays the 

first few lines of the text file. Text files containing data often have a few descriptive lines before the 

41 
Figure 18. Text File Import Parameters Dialog box. 

data begins. These must be ignored. The first column of numbers must be the abscissa (time or 

frequency). The user can specify which column to begin importing the range data (amplitude, 

magnitude and phase, or real and imaginary). Data can be imported in the time-domain or the 

frequency-domain. The fast Fourier transform (FFT) is an integral part of the GP Wave Utilities 

program. Therefore, a few data samples with the final sample value are concatenated to the imported 

data so that the proper number of data points are available for FFT operation. This usually involves 

adding a small number of samples (~10) because the prime factor FFT is used. Data are saved in the 

GP Wave Utilities *.fdw format. 

When exporting data, an ASCII text file is written. The contents of this file depend on whether a 

time-domain or frequency-domain plot is currently displayed. For time-domain data, two columns of 

data are written: time and amplitude. For frequency-domain data, five columns of data are written: 

frequency, real component, imaginary component, amplitude, and phase. 

42
Process 

The menu items in the Process menu section are where most of the GP Wave Utilities functions 

can be found. 

Wave Plot, Filters and Operators 

This menu item brings up the dialog box shown in figure 19. Select one of the operations on the 

left side of the dialog and enter the corresponding parameters in the relevant subsections. Two 

frequency filters are available – a finite impulse response (FIR) filter and a fast Fourier transform (FFT) 

based filter. The FIR filter convolves a 100-point filter function with the time domain waveform. The 

FFT filter applies the filter in the frequency-domain. Note that since FFTs require a circularly 

continuous function, if the first and last samples have a significantly different level, then a step is 

Figure 19. The Filter Dialog box. 

43
present in the data. If this is the case, the frequency-domain representation data will have energy across 

the band because a step function contains all frequencies. To avoid this effect, use the Process, Wave 

Plot, Subtract Ramp menu item. If any FFT operations are anticipated for the data in the document, 

then FFT filters should be used rather than the FIR filters. The cosine-taper time window weights the 

data with a cosine taper in time. This window has a zero value at times before the begin time and after 

the end time. The value is one at times between the begin time and end time. The taper frequency is 

used to specify the cosine-taper frequency and thus specifies the rise and fall time of the window. The 

resample operation resamples the data over the specified time window and sample rate. If the time 

window extends beyond the range of the current data, then the data are padded using the nearest existing 

data point. Data are removed if the time window is smaller than the data extent. New resampled data 

points are determined by interpolating between existing data points. The linear operator operates on the 

x-axis values (e.g., time or frequency) or the y-axis values (e.g., amplitude). The operation is given in 

equation 5, 

newx = (multiplier ⋅ x) + offset + (squareroot ⋅ x) + (log10 ⋅ log10(x)) . (5) 

Equation 5 is applied to every data point in the document. Note that the newy and y are used in place of 

newx and x in equation 5 when operating on the y-axis values. The frequency-dependent transmission 

coefficient (T) can be calculated from the frequency-dependent reflection coefficient (R) that is loaded 

into the current document. The intrinsic impedance Z of layers one (contains incident and reflected 

waves) and two (contains transmitted wave) must be specified. The transmission and reflection 

coefficient calculations are for waves normally incident on a planar interface. Alternatively, the 

frequency-dependent reflection coefficient (R) can be calculated from the frequency-dependent 

transmission coefficient (T) loaded into the current document. A waveform can also be rotated or 

shifted in time. When rotating a waveform, the last part of the waveform is rotated to the beginning part 

44
of the waveform (or vise-versa depending on the direction of rotation – positive or negative). When 

shifting a waveform, the waveform is padded on one side using the nearest sample value, and truncated 

on the other side. The final operation in this dialog box, the moving average, convolves the waveform 

with a boxcar filter function having the specified number of points. In all cases, the current waveform 

in the current document will be replaced by the resulting waveform. 

Wave Plot, Subtract Ramp 

Use this menu item to subtract a ramp from the waveform so that the first and last sample both 

have the same value. If these samples have the same value, the broadly distributed frequency 

components of a step function will be added to the frequency-domain spectrum of the waveform, which 

may cause problems with some waveform processing operations. This happens because an FFT 

assumes that a waveform is periodic. Ramp subtraction is a routine part of processing many 

waveforms. The current waveform in the current document will be replaced by the resulting waveform. 

Wave Plot, Subtract Waveform 

Use this menu option to subtract a waveform (sample by sample) in another document from the 

waveform in the current document. Both waveforms must have the same sample interval and the same 

number of samples. If necessary, use the re-sample function (Process, Wave Plot, Filters and 

Operators menu item) to re-sample the waveforms so that they have the same sample interval and the 

same number of samples. The user will be asked to select a currently opened document containing the 

waveform to subtract from the current waveform. This is a time-domain operation. The current 

waveform in the current document will be replaced by the resulting waveform. 

45
Wave Plot, Remove DC 

The menu item subtracts the average or DC value from the current waveform. This is a timedomain operation. The current waveform in the current document will be replaced by the resulting 

waveform. 

Wave Plot, Deconvolve 

Select this menu item to deconvolve (in the time-domain) the specified waveform from the 

current waveform. The user will be asked to select a currently opened document containing the 

waveform to deconvolve from the current waveform. Deconvolution can be an unstable procedure 

when parts of the frequency-domain spectrum of the waveform to deconvolve have low amplitudes in 

comparison to the other frequency-domain components. Dhaene et al. (1994) describe an adaptive filter 

to reduce many of the problems associated with frequency-domain deconvolution of time-domain data. 

The realization of their filter used in this work is given by 

F(ω,γ , λ) = 

| x(ω)

| 

| 

x 

2

(ω

+γ 

) 

ω 

|

2 

2 + λ 

, (6) 

where λ is a pre-whitening or peak reduction parameter, and γ is a low pass filter or smoothing 

parameter. In the resulting deconvolution, 

h(ω) = F(ω,γ ,λ) y(ω) , (7)

x(ω) 

the peak reduction parameter reduces the effect of noise at frequencies where x(ω) is small, and the 

smoothing parameter reduces the effect at high frequencies where most of the power in x(ω) is noise. In 

this work, the parameters λ and γ are selected to maximize the power in the resulting deconvolution at 

frequencies below 200 MHz, and minimize the power above 400 MHz. The assumption is that the most 

power in x(ω) and y(ω) is distributed over the lower frequencies, and that the power at higher 

46
frequencies is mostly noise. This assumption is based on the fact that most time-domain 

electromagnetic equipment produces a negligible amount of energy above a certain frequency (about 

250 MHz for the RTDGPR built by the USGS). The code can be easily modified to use other 

frequencies or to ask the user to provide these frequencies. The current waveform in the current 

document will be replaced by the resulting waveform. 

Wave Plot, Convolve 

Use this menu item to convolve two traces in the time-domain. The user will be asked to select 

a currently opened document containing the waveform to convolve with the current waveform. The 

current waveform in the current document will be replaced by the resulting waveform. 

Wave Plot, Correlate 

Use this menu item to correlate two traces in the time-domain. The user will be asked to select a 

currently opened document containing the waveform to correlate with the current waveform. The 

current waveform in the current document will be replaced by the resulting waveform. 

Wave Plot, Envelope 

This menu item is used to calculate the envelope of the current waveform. Two methods are 

available. The first method, the peak method, simply calculates the average value of the absolute value 

of the waveform using a moving boxcar window of the user-specified length. For the second method, 

the Hilbert envelope method, the waveform in the current document and its Hilbert transform are 

combined in quadrature in the frequency domain. This method generally produces a more compact 

envelope than the peak method. The current waveform in the current document will be replaced by the 

resulting waveform. 

47
Wave Plot, Hilbert transform 

Use this menu item to calculate the Hilbert transform of the waveform in the current document. 

The current waveform in the current document will be replaced by the resulting waveform. 

Wave Plot, Path Length 

Use this menu item to calculate the RMS path length of the waveform. This is done by summing 

the quadrature path lengths of all line segments making up the waveform. The quadrature path length is 

found by combining the vertical and horizontal lengths of a line segment in quadrature. The current 

waveform in the current document will be replaced by the resulting waveform. 

Wave Plot, Amplitude Response 

Use this menu item to build a non-linear function that maps one waveform (e.g., an input 

waveform to an amplifier) to another (e.g., the waveform output from the amplifier). This function can 

be used to build a non-linear amplitude response function for a signal-processing module such as a 

logarithmic amplifier. The currently selected document must contain the output of the signal-processing 

module. The user is prompted to select another opened document containing the input to the signalprocessing module. The time scale and sample interval must be the same in both waveforms. The 

sample times of the output waveform should be adjusted so that the waveform value at a given time 

represents the signal-processing module output for the input at the same time. By comparing the input 

amplitude to the output amplitude at each sample, an amplitude response function is built. After 

calculating the amplitude response function, the response function is displayed. In the response 

function plot, the input amplitude is plotted on the vertical axis, and the output amplitude is plotted on 

the horizontal axis. The resulting response function can be used to linearize the non-linear response of a 

signal-processing module (see the Process, Wave Plot, Amplitude Linearization menu item). For 

48
example, after calculating the amplitude response of a logarithmic amplifier, the logarithmic 

compression imposed by the amplifier can be removed from the signal. Figure 20 shows the input and 

output waveforms for a signal-processing module, and the calculated amplitude response function. 

Figure 20. The left panel shows an input waveform (red) and an output waveform (blue). The 

right panel shows the calculated response function for the signal-processing module. 

Wave Plot, Amplitude Linearization 

Use this menu item to calculate the input to a non-linear signal-processing module (e.g., a 

logarithmic amplifier) given the module’s output waveform and the non-linear response function. The 

output waveform should be loaded in the current document. Then the user will be prompted to select 

the opened document containing the non-linear response function that was calculated using the Process, 

Wave Plot, Amplitude Response menu item, as shown in the right side of figure 20. 

Wave Plot, Packet Time and Amplitude 

Use this menu item to calculate the RMS value of the waveform, and the temporal centroid of 

the RMS value (i.e., the time in which half of the contribution to the total RMS value is due to earlier 

parts of the waveform). These calculated values can be used in subsequent data processing, where, for 

instance, waveforms need to be compared or arrival times are needed. 

49
Wave Plot, Compare Waves 

This menu item was designed to compare various early time waveforms produced by the 

RTDGPR (USGS). The user will be prompted to select an opened document containing a waveform to 

compare with the waveform in the current document. This menu item makes a temporary copy of each 

waveform, and temporally aligns them at their largest positive peak. Then waveform data after 40 ns 

are set to zero by a 125 MHz cosine filter. Next, the selected waveform is scaled so that it has the same 

RMS value as the waveform in the current document. Finally, the procedure reports the following 

values: the time shift needed to align the first positive peak of each waveform, the RMS difference 

between the waveforms, and the scaling factor used. This menu item only makes a comparison between 

the waveforms. No data are modified. 

Image Plot, Extract Trace 

Use this menu item to extract a trace from an image plot (a.k.a. wave field plot; see fig. 16), so 

that it can be viewed and manipulated like a waveform (see fig. 15). An image is comprised of many 

horizontal traces that are vertically adjacent to each other. The user is prompted for the trace number to 

extract. The trace numbering begins at the bottom of the screen. The currently displayed trace numbers 

can be selected using the Plot Parameters Dialog box (see fig. 17). The plot parameters dialog box is 

invoked by right-clicking on the plot. 

Image Plot, 3D Depth Shift 

Use this menu item to mathematically propagate a digital representation of an electromagnetic 

wave field through a user-specified medium. The current document must contain wave field data (see 

the Process, Calc TE and TM Spectrum menu item). Depending on the size of the wave field, this depth 

shift operation may require a PC with 1 GB or more of RAM so that efficient calculations can be made 

50
(approximately 20 seconds on a 3 GHz Pentium). The calculation may be exceedingly lengthy if 

insufficient RAM is available. The required RAM is a function of the spatial and temporal sample 

density of the wave field. It is assumed that the wave field in the current document represents an 

electromagnetic wave field on a horizontal scan plane as shown in figure 21. The scan plane is an X-Y 

plane for a constant Z (or depth) value. It is further assumed that the waves are propagating downward. 

The user will be presented with the Depth Shift Dialog box shown in figure 22. Enter the Cole-Cole 

material properties of the medium, and depth shift value to use. A new wave field will be calculated for 

the horizontal plane whose depth is the depth of the original scan plan plus the depth shift. 

Ground Surface 

Transmitting Antenna 

Lateral Wave 

Air Wave 

Scan Plane 

Spherical Wave 

Figure 21. Section view illustrating subsurface wave fronts and scan plane beneath a GPR 

antenna (adapted from Oden, 2006). 

51
Figure 22. The Depth Shift Dialog box. 

Image Plot, 2D Depth Shift 

This menu item is essentially identical to the Process, Image Plot, 3D Depth Shift menu item, 

except that the calculations are made in two dimensions – X and Z. Wave spreading and attenuation are 

not accounted for in the Y direction. This calculation is faster than the three-dimensional calculations, 

but they are not accurate. Despite this, the two-dimensional calculations can provide useful 

information. The two-dimensional results are especially inaccurate when the scan plan or the depthshifted plane are in the near field zone of the antenna (see fig. 21). See Oden (2006) for a discussion on 

the near field regions of GPR antennas. 

Image Plot, System Response a(w)T(k,w)R(k,w) 

This menu item can be used to calculate the frequency-domain system response of a GPR system 

due to a perfect planar reflector at an arbitrary depth as shown in figure 23. Here, the system response is 

defined as the received waveform produced by a given GPR system over a perfect planar reflector at a 

specified depth, and for specified soil properties. This response includes the signal sent to the 

52
transmitting antenna a(ω), the transmitting antenna response function T(k, ω), and the receiving antenna 

response function R(k, ω), where ω is radian frequency, and k is spatial frequency. Spatial frequency is 

accounted for in the X and Y directions, but is shown as k here for simplicity. The product a(ω)T(k, ω) 

is the frequency-domain wave field at the reflector. The product a(ω)T(k, ω)R(k, ω), the system 

response, is the frequency-domain waveform at the receiving antenna. To be concise, these response 

functions and their relations are written symbolically and many mathematical and physical details have 

been omitted. Consult Oden (2006) for a more detailed treatment. The system response for a given 

GPR system over a perfect planar reflector at a specified depth, and for specified soil properties can be 

calculated using the following steps. 

1) Build a scan plane wave field (see figure 21 and the Process, TE and TM Spectrum menu item). 

This scan plane wave field is the product a(ω)T(k, ω) for the waves at the scan plane. 

2) Use the Process, Image Plot, 3D Depth Shift menu item to shift the wave field to the desired 

depth. The wave field now represents the product a(ω)T(k, ω) for the waves at the reflecting 

plane. 

3) Use the Process, Image Plot, System Response a(w)T(k,w)R(k,w) menu item to calculate the 

product a(ω)T(k, ω)a(ω)R(k, ω)/Z0, where Z0 is the characteristic impedance of the transmission 

line connected to the receiver antenna terminals (usually 200 ohms for GPR antennas using 

banana jack style connectors). It is assumed that the receiving antenna is identical to the 

transmitting antenna. 

4) Deconvolve a(ω) and multiply by Z0 to obtain the received waveform a(ω)T(k, ω)R(k, ω). 

5) If needed, convolve by the response of the receiver electronics b(ω) to obtain a(ω)T(k, ω)R(k, ω) 

b(ω). 

53
Ground 

surface 

at z=0 

Antennas 

Scan 

plane at 

depth z0 

Horizontal 

planar 

reflector at 

depth z 

Figure 23. Diagram of a GPR system response over a perfect horizontal planar reflector. 

This procedure can be used to deconvolve for the electrical properties of a horizontal reflector. 

This deconvolution procedure is outlined in figure 24. For this procedure, the fields on the scan plane 

are determined by numerical simulations of the GPR antennas system, and laboratory measurements of 

the signal sent to the transmitting antenna. 

Determine fields on 

scan plane below 

surface from a priori 

numerical simulations 

Calculate a(ω)T(k, ω) 

R(k, ω)b(ω) for 

perfect reflector at 

estimated depth z 

Calculate correlation 

between recorded 

waveform and 

a(ω)T(k, ω)R(k, ω) 

b(ω) 

Is correlation 

minimum phase? 

Deconvolve a(ω)T(k, 

ω)R(k, ω) b(ω) with 

extracted wavelet to 

estimate reflection 

coefficient and 

reflector properties 

Extract wavelet 

reflected subsurface 

plane from recorded 

data 

Update estimated 

depth z based on 

phase of correlation N Y

Add time offset to 

recorded data so that 

transmitter fires at 

fire time used in 

simulations 

Figure 24. Flow chart for deconvolving for reflector properties. 

54
Radiation Pattern 

This menu item invokes an experimental routine to estimate a two-dimensional radiation pattern 

from a GPR antenna. Note that this routine is experimental, and the results have not been carefully 

verified. The routine reads a simulated wave field that has been recorded on a scan plane below a GPR 

antenna (see fig. 21). Accurate simulations of actual GPR antennas can be made using the finitedifference time-domain methods (Oden, 2006; Giannopoulos, 2004). The routine imports a twodimensional wave field (X versus time) from data files recorded during the simulation. The imported 

wave field is converted to the frequency-domain representation. The frequency-domain wave field is 

comprised of a collection of plane waves, each with a certain frequency and traveling in a certain 

direction. The routine adds up the amplitude of all plane waves traveling in each direction to obtain the 

radiation pattern. 

The GPRMax simulations must be set up to record the fields on the scan plane using snapshot 

files (*.snp). The simulations calculate the wave field at a series of times steps (e.g., 0 to 100 ns in 20 

ps intervals). The wave field values on the scan plane for each time step are recorded in a snapshot file. 

A different file is written for each time step with the name *_nnnn.snp, where nnnn is the index of the 

time step. If 1,000 time steps are used in the simulation, the first file will be named *_0000.snp and the 

last *_0999.snp. The scan plane, as defined by the field values recorded in the snapshot files, must be 

centered beneath the transmitting antenna, and have the same spatial extent in the X and Y directions. 

To begin the calculation, select the *_0000.snp file. The process may require about an hour to run 

depending on the number of time steps in the simulation, the size of the scan plane, and the spatial 

sample interval. After calculating the radiation pattern, the pattern is displayed as a waveform. A 

radiation pattern with poor resolution will result from using a scan plane that is too small. 

55
Ground Properties 

Use this menu item to estimate the shallow soil properties under ground-coupled GPR antennas 

using the early-arriving portions of the received waveform. This procedure requires an existing ground 

properties model that specifies the response of the antennas as a function of the soil properties (see the 

Process, Ground Properties Model menu item). The user will be presented with the Invert for Ground 

Properties Dialog box shown in figure 25. The first step is to specify the waveform(s) to use in 

estimating the soil properties. The waveform in the currently loaded document can be used. 

Alternatively, the waveforms from a number of GPRMax (Giannopoulos, 2004) simulations (*.out files) 

Figure 25. The Invert for Ground Properties Dialog box. 

56
can be used. The first waveform listed in each selected *.out file will be used. Multiple *.out files can 

be selected. Finally, the waveform(s) in a RTDGPR (USGS) data file can also be used. The user can 

either specify a single trace number, or specify that all of the waveforms in the file should be used. 

When multiple waveforms have been specified, a set of soil property estimates will be made for each 

waveform. 

Next, the preprocessing parameters must be specified. Preprocessing may be needed to make 

the waveforms in the specified input files represent the same situations as the processed waveforms that 

were used to generate the ground properties model. For instance, if the input file is the simulated 

response of the antenna system, then it may be necessary to convolve the simulated waveform with the 

response of the receiver electronics. An amplitude scale factor can also be specified. Actual data may 

need to be time shifted before they can be used with the selected ground properties model. For the case 

of the RTDGPR (USGS), a fiducial marker is placed on the recorded waveform to indicate when the 

transmitter fired. An amplitude threshold must be specified so that the fiducial can be found. The 

fiducial is interpreted to occur at the time of the first peak that is higher than the specified amplitude 

threshold. 

Next, a ground properties model file must be selected. The ground properties model files 

contain information specifying the waveform attributes to extract, and how these attributes vary with 

changing soil permittivity, conductivity, and antenna height above the ground. The parameters to find 

in the inversion must be specified. Sometimes, one of the parameters (soil permittivity, conductivity, 

and antenna height above the ground) is already known (perhaps through independent measurement). 

This helps to constrain the inversion. The uncertainty between the measured waveform attributes and 

the predicted waveform attributes must be specified. This uncertainty includes noise in the data, the 

difference between the GPRMax simulated attributes and the actual attributes, and the descretization 

57
error in the forward operator (i.e., the ground properties model). Specifying each of these uncertainty 

components as an RMS difference, and adding them in quadrature leads to the RMS difference between 

the predicted waveform attributes and the recorded waveform attributes. Normalizing this RMS 

difference by the RMS attribute amplitude of the recorded waveform gives the normalized uncertainty 

to enter in the dialog box. When the normalized uncertainty value is unknown, a value of 0.1 (10%) can 

be used. A provision for selecting an uncertainty components file has been provided for future 

development. For a more detailed discussion about uncertainty and the process of inverting for soil 

property estimates, consult Oden (2006). Finally, specify the information to record in the output file 

(output type) and press OK to begin the inversion. If Internal inversion details is selected, then the 

output file will contain the convergence information at each iteration of the inversion process. If 

Resolution details is selected, then the output file will contain the parameter estimates for each input 

waveform. If only one input waveform is selected, then the results will be presented on screen. 

Ground Properties Model 

Use this menu item to calculate a response model that can be used to estimate the soil properties 

just beneath a pair of ground-coupled transmitting and receiving GPR antennas. A suite of simulations 

must have been run using GPRMax (Giannopoulos, 2004) to model the antenna pair response to a 

variety of soil properties. Each GPRMax simulation models the response of a transmitting-receiving 

antenna pair, and the results of each simulation are written to a *.out file. The *.out files can contain 

electromagnetic field values for many spatial locations in the simulation grid. When using these 

simulation files to calculate ground response, the first waveform in each *.out file must be the electric 

field waveform at the receiving antenna feed port. A suite of simulations are needed that model the 

antenna response for various ground properties (permittivity and conductivity), and the antenna height 

above the ground. Figure 26 graphically depicts the range over which the ground properties and the 

58
antenna standoff (i.e., antenna height above the ground) can vary for a given ground properties model. 

This is known as the model space grid. For the soil properties and antenna standoff represented at each 

grid corner, a simulation has been run and the antenna response is known. This collection of 

simulations is used to create a response model that can be used to estimate ground properties beneath 

GPR antennas in actual surveys. Once the shallow soil properties are known beneath GPR antennas in a 

given location, the antenna response can be calculated. This in turn enables other data processing 

techniques such as deterministic deconvolution for subsurface reflector properties. Consult Oden 

(2006) for a detailed discussion on estimating shallow soil properties beneath ground-coupled GPR 

antennas, determining the antenna response at a given location, and the data processing techniques that 

these steps enable. 

Figure 26. The model space grid. The simulated response is known at the corners of each 

grid cell. 

59
The Process, Ground Properties Model menu item invokes the Model Ground Properties Dialog box 

shown in figure 27. The first step is to choose a suite of simulations files. The simulation files should 

cover a uniform model space such as the one shown in figure 26. The file names containing GPRMax 

simulation results that will be used to create the ground properties model must follow a prescribed 

naming convention. GP Wave Utilities requires the soil properties and antenna standoff used 

Figure 27. The Model Ground Properties Dialog box. 

60
in the simulation to be encoded in the file name. The beginning of the file name can be anything, but 

the end must be written as follows 

*_p9_c0_02_s7.out, 

where the number following the p is the RDP, the number following the c is the conductivity (in mS/m), 

and the number following the s is standoff (in cm). An underscore is used in place of the decimal point 

in the conductivity value, and to separate fields. The file extension must be *.out. Next, a variety of 

processing steps are needed to convert the simulated antenna response to the response of the entire GPR 

system. An amplitude scaling factor can be set. A response waveform can be convolved with the 

simulated response. For example, the response of the receiver electronics and/or the signal presented at 

the transmitting antenna feed port can be convolved with the simulated antenna response. A signal can 

also be deconvolved. This may be necessary when the transmitting antenna excitation used in the 

simulation is different than the actual excitation (i.e., the transmitting antenna feed port signal). Next, a 

user-specified finite-impulse-response band-pass filter can be applied to the simulated waveforms. 

Finally, a 125 MHz cosine taper window can be used to eliminate energy outside of the specified time 

window. For example, one may eliminate energy at late times when it is undesirable to use the late time 

energy to predict ground properties. After pre-processing the simulated waveforms, attributes of the 

waveforms will be extracted for the ground properties model. The ground properties model specifies 

the waveform attributes as a function of the soil properties and antenna standoff. This ground properties 

model can be inverted to estimate the soil properties under GPR antennas from survey data (see the 

Process, Ground Properties menu item). A ground properties model can be made using different 

attributes sets (i.e., attributes of a given type) of the simulated waveforms. Table 5 lists the various 

methods for extracting waveform attributes. The FD amplitudes and the 20 data points on envelope 

methods provide the best results. 

61
After specifying the input files and processing parameters, press the save button to calculate a 

ground properties model and save the model to disk. A file save dialog box will be presented so that the 

user can specify the file name (a *.txt file) in which to store the model. Pressing the open button allows 

an existing ground properties model to be opened. A two-dimensional contour map of each waveform 

attribute can be made using the view button. This button uses Surfer (Golden Software) to generate the 

contour maps. Each waveform attribute varies as a function of soil permittivity, conductivity, and 

antenna standoff (i.e., three dimensions). A contour map of the attributes can be made as a function of 

any two of these parameters by adjusting the view parameters. For example, selecting Eps vs Sig plots 

and specifying a view standoff of 7 would generate a contour plot of each attribute where the horizontal 

axes of the graph are soil permittivity (epsilon) and conductivity, and the standoff is held to a constant 

value of 7. By changing the view parameters, the character of the changing attributes can be viewed on 

Table 5. Methods of extracting waveform attribute sets. 

Attribute Set Description Results 

FD amplitudes FFT spectral amplitudes from 10-200 MHz. If the 

waveform is 100 ns long, an amplitude value will be 

available at 10 MHz intervals. 

Good 

TD amplitudes Peak amplitudes of waveform after applying a band pass 

filter bank. The filter bandwidths were 10 MHz and 

Marginal 

TD times 

TD amplitudes and times 

Time and amplitude of 

first 4 peaks 

RMS amplitude and time 

centroid 

centered at frequencies between 10-200 MHz at 

approximately 10 MHz intervals. 

Time of energy peak of waveform after applying a band 

pass filter bank. The filter bandwidths were 10 MHz and 

centered at frequencies between 10-200 MHz at 

approximately 10 MHz intervals. 

All parameters from the previous two methods 

Time and amplitude of first four peaks. 

RMS amplitude and time of energy peak of waveform. 

Poor 

Poor 

Poor 

Poor 

20 data points 20 samples evenly spaced throughout the selected time 

window. 

Poor 

20 data points on 

envelope 

20 samples of the Hilbert envelope of the waveform evenly 

spaced throughout selected time window. 

Good 

62
any plane (as long as one parameter remains fixed) that cuts through the model space shown in figure 

26. This is useful for investigating the topography of a ground properties model (i.e., a forward model), 

and assessing how an inversion based on this forward model would perform. 

Compare Waveforms 

Use this menu item to compare waveforms measured in the field or laboratory with simulated 

waveforms. This is useful, for instance, when testing several GPR antenna models to find the one that 

produces simulation results similar to the actual GPR antenna response. This menu item produces the 

Compare Waveforms Dialog box (see fig. 28). This is a modeless dialog box, meaning that it does not 

have to be dismissed to access the other menu items or toolbar buttons. Click and drag on the blue 

header at the top of the dialog box to move the dialog box as needed. The dialog box will compare a 

single waveform from a field or laboratory test with a number of simulated waveforms. Use the 

corresponding select buttons to invoke the file selection dialogs. After selecting the files to be 

compared, check the invert boxes as necessary (the polarity of the simulated waveform may be different 

than that of the recorded waveform). The only supported file format for the field or laboratory data is 

the *.aq1 file recorded by the RTDGPR (USGS). The only simulation files supported at this time are 

*.out files created by GPRMax. The *.out files contain waveform data for the electromagnetic fields at 

many locations in the simulated model space. Enter the grid indices of the simulation grid 

corresponding to the physical location for the electric field waveform of interest (e.g., the receiving 

antenna port). Further operations to the simulated data may be needed before they can be compared 

with the recorded data. For example, the simulated data may only account for the antenna response, and 

the response of the receiving electronics must be included before a proper comparison to the recorded 

data can be made. Specify *.txt files containing waveforms to deconvolve and convolve if necessary. 

These waveforms should be written as ASCII text files where the first column contains the time values, 

63
Figure 28. The Compare Waveforms Dialog box. 

64
and the second column contains the data values. Note that any waveform that has been imported or 

created by GP Workbench can be exported in the ASCII text file format (see the File, Export menu 

item). The next step is to synchronize the waveforms using either the positive peak method or by 

correlation. The positive peak method adds a time shift to the simulated data so that the largest positive 

peak in the simulated waveform is aligned with the largest positive peak in the recorded data. The 

correlation method adds a time shift to the simulated data until a zero lag correlation is obtained 

between the simulated and recorded data. Next, a cosine taper can be applied to remove portions of the 

waveforms that should not be compared. For instance, a cosine taper could be used to remove a late 

time reflection from the recorded waveform that is not present in the simulated data. Finally, the 

simulated data can be amplitude scaled before being compared with the recorded data. If a value of zero 

is entered, an amplitude scaling factor will be calculated that causes the simulated data to most closely 

match the recorded data. 

After processing the data sets, the comparison can be made. Four methods of comparing the 

waveforms are available. The time-domain method simply calculates the RMS amplitude difference 

between the waveforms on a sample-by-sample basis. The time-domain envelope method calculates the 

Hilbert envelope of each trace before calculating the RMS amplitude difference. The frequency-domain 

method takes the fast Fourier transform of each trace, and then calculates the RMS amplitude difference 

on a frequency-by-frequency basis. Finally, the trace length method determines the trace length (i.e., 

the sum of the length of each line segment making up the waveform) of each waveform, and compares 

these lengths. When all of the parameters for processing the waveforms and making the comparisons 

are set, press the process button to make the comparison. The amplitude-scaling factor and the 

difference between the simulated data and the recorded data will be displayed in the results pane. If the 

create waveform documents box is checked, then a GP Workbench document will be created for the 

65
recorded data file and each simulated data file specified in the dialog box. The waveforms will be 

processed according to the parameters in the Compare Waveforms Dialog box before a document is 

created. A visual comparison between waveforms can be made by creating waveform documents, and 

then plotting multiple waveforms on the same graph (see the Plot Parameters Dialog box in Figure 17). 

TM and TE Spectrum 

Use this menu item to import a scan plane wave field (see fig. 21) from a set of numerical 

simulations. The simulations are needed to model the electromagnetic fields on a scan plane beneath 

the transmitting antenna of a GPR system. Accurate simulations of actual GPR antennas can be made 

using the finite-difference time-domain methods (Oden, 2006; Giannopoulos, 2004). These simulations 

can only directly calculate the fields near the antenna. If the fields further from the antenna are desired, 

then a field propagator can be used (see the Process, Image Plot 3D Depth Shift menu item). To import 

a scan plane wave field, GP Wave Utilities reads in the time- and space-domain wave field values from 

the files produced in the simulation. During import, these time-domain and space-domain wave field 

values are converted to frequency-domain field values. An electromagnetic wave field can be 

efficiently described in terms of the transverse magnetic (TM) spectrum, and the transverse electric (TE) 

spectrum. The TM and TE spectra are a subset of the frequency-domain wave field values. The salient 

point is that this import process converts the wave field into a more compact and readily usable format, 

which in turn facilitates wave field calculations such as those made with the Process, Image Plot, 3D 

Depth Shift menu item. 

The GPRMax simulations must be set up to record the fields on the scan plane using snapshot 

files (*.snp). The simulations calculate the wave field at a series of times steps (e.g., 0 to 100 ns in 20 

ps intervals). The wave field values on the scan plane for each time step are recorded in a snapshot file. 

A different file is written for each time step with the name *_nnnn.snp, where nnnn is the index of the 

66
snapshot. If 1000 time steps are used in the simulation, the first file will be named *_0000.snp and the 

last *_0999.snp. The scan plane, as defined by the locations of the field values recorded in the snapshot 

files, must be centered beneath the transmitting antenna and have the same spatial extent in the X and Y 

directions. To begin the import process, select the *_0000.snp file. The import process may require a 

few hours to run depending on the number of time steps in the simulation, the size of the scan plane, and 

the spatial sample interval. After importing the wave field, it can be saved as a GP Wave Utilities file 

(*.fdw) for storage to disk and much quicker retrieval. For a detailed description of this process, refer to 

Oden (2006). 

View, Window, and Help Menus 

These menus contain the standard items for a Windows application. 

The Toolbar 

The toolbar shown below provides many operations for viewing and processing data. Each 

button on the toolbar will be described in more detail below. 

Figure 29. The GP Wave Utilities toolbar. 

The new document button generates a new empty document. A document must exist in which to 

import data files. 

The Open-File dialog can be used to open GP Wave Utilities document files (*.fdw files). 

The save document saves the current document. If a file has not yet been associated with the 

document, the user is prompted for a file name. 

The print button prints the current document using the default printer. 

67
This button displays information about GP Wave Utilities such as the version number. 

This button shows the wave field in the x-plane in either the time domain (x, t) of the frequency 

domain (kx, ω). This button only functions with the wave field data. 

This button shows the wave field in the y-plane in either the time domain (y, t) of the frequency 

domain (ky, ω). This button only functions with the wave field data. 

This button shows the magnitude of the wave field in the frequency domain. This button only 

functions with the wave field data. 

This button shows the real part of the wave field in the frequency domain. This button only 

functions with the wave field data. 

This button shows the imaginary part of the wave field in the frequency domain. This button 

only functions with the wave field data. 

Data Processing Methods 

The processing steps used to process each geophysical data set are different. There are, 

however, some basic processing outlines that are broadly applicable to many data sets. In this section, 

processing overviews are given for GPR and ALLTEM datasets. 

GP Workbench GPR Data Processing Procedures 

The goals of data processing vary. Even so, there are basic processing steps that are commonly 

followed when processing GPR data. An outline for section view processing and for plan view 

processing are presented below. 

68
GPR Section View Processing 

The following list provides a good start for processing a single file in section view. 

1) Import the data file in section view. 

2) After importing the data, it needs to be gridded. The exact procedure will vary depending on the 

type of data and the gridding algorithm used. If marked positions are used, the user should make 

sure that all marks are present, and no extra marked positions are included. See the toggle marks 

button on the toolbar section. Specify the grid parameters using the Grid Parameters Dialog 

box (see fig. 3). Then press the grid button on the toolbar. 

3) Click on the plot and use the up and down arrow keys to position the cross hairs at time zero 

(when the transmitter fired). Press the zero time button on the toolbar. 

4) Filter the image to improve the contrast of reflections with background and to remove noise. 

Use the Process, Filter menu item. Commonly applied filters include background removal and 

frequency filters. 

5) Perform a velocity analysis. Locate a reflector on the section that exhibits hyperbolic move out. 

Click on the plot and use the arrow keys to position the cross hairs on the apex of the hyperbola. 

Next, click the velocity hyperbola button. Use the arrow keys to adjust the shape of the fitting 

hyperbola to match the diffraction hyperbola pattern. The velocity and index of refraction will 

be displayed in the status bar on the bottom of the main window. Click the velocity hyperbola 

button to exit the hyperbola fitting mode. Now, click on any point on the cross section and read 

the depth of the cross hairs from the status bar. 

GPR Plan View Processing 

The following steps outline a basic procedure for processing plan view GPR data. 

69
1) Import the data files in plan view. If specific positional processing must be applied to each file, 

then each file should be opened in section view, gridded, and saved. These processed files can 

then be imported together in plan view. 

2) After importing the data, it may need to be gridded. Usually, plan view data are not re-sampled 

on a uniform grid. It is important to have interpolated positions calculated between marked 

traces before continuing. This can be done on a line-by-line basis in section view before 

importing the files in plan view, or it can be done in plan view. Specify the grid parameters 

using the Grid Parameters Dialog box (see fig. 3). Then press the grid button on the toolbar. 

3) Filter the image to improve the contrast of reflections with background and to remove noise. 

Use the Process, Filter menu item. Commonly applied filters include background removal and 

frequency filters. 

4) Use the up and down arrows to select the time-slice of interest, then select the View, Surfer menu 

item. This will generate a high quality Surfer (Golden Software) image. See the View, Surfer 

menu item for more details. 

GP Workbench ALLTEM Data Processing Procedures 

The GP Workbench program defines a UXO target in terms of three orthogonal principal 

polarizabilities. When specifying these dipole moments, the attitude of the first (largest) principal 

polarizability is specified. The direction of second dipole is determined from the cross product of the 

first polarizability with the vertical axis. The direction of third polarizability is the cross product of the 

first two polarizabilities. Only one time constant is used to simulate the transient response. This is the 

effective time constant at given late time, and the dipole moments are the magnitude of the response at 

this time. When inverting actual data, data for the transient response for all polarizations are extracted 

70
at the same moment in late-time. This moment is when the RMS amplitude of all channels is equal to 

25% of the RMS peak value of all the channels. 

The data processing procedures outlined below require that calibration data be recorded using 

the steel ball after the system has warmed up and the response has stabilized. 

Calibration Procedure 

1. Import data files recorded at a fixed location to obtain the system response. Import 

*.h1m and associated documents. Run a 5-point boxcar spatial filter. Select a trace to 

use in the calibration and run a fixed window spatial filter using this trace (i.e., replace 

all traces with selected trace). Remove all traces except the first 10 (use the Filter 

Dialog box). Use Save As to save files as *.h1m.SU. 

2. Import data files recorded at a fixed location (same location that the system response was 

measured) where the response of the steel reference ball was measured. Import *.h1m 

and associated documents. Run a 5-point boxcar spatial filter. Select a trace to use in the 

calibration and run a fixed window spatial filter using this trace (i.e., replace all traces 

with selected trace). Remove all traces except the first 10 (use the Filter Dialog box). 

Use Save As to save files as *.h1m.SU. 

3. Concatenate system response and steel shot response files by importing the *.h1m.SU 

files saved in steps 1 and 2 above. The files will be concatenated in the order they appear 

in the import file select list. Make sure that system response file is the first in the list. 

4. Subtract background response. Use the Filter Dialog box, remove spatial window, fixed 

window 2-3. Process associated files concurrently. 

5. Remove all records in each file except the last 5 (use the Filter Dialog box). 

71
6. Select Process, ALLTEM Calibration from the menu. This will write the calibration file 

Cal.h1m.cal. 

UXO Target Inversion Procedure 

1. Import data files in section view. Remove noisy traces (using toolbar button). Process 

associated files concurrently. Use Save As to save as (*.h1m.SU, and associated files). 

2. Import data files in plan view. Remove traces with bad GPR readings (using toolbar 

button). Process associated files concurrently. Use Save As to save as (*.h1m.SU, and 

associated files). 

3. Make an import.txt file for the directory containing the data. This file contains the 

coefficients defining the linear mapping of latitude and longitude to northing and easting. 

4. Pick a target to process. Identify the lines containing the response from this target. For 

each line, do steps 3-4. It is assumed that each line has leading and trailing data with no 

targets from which the system response and the soil response can be determined. 

5. Subtract background response. Determine trace numbers for leading data and trailing 

data using viewer. Use the Filter Dialog box, remove spatial window, interpolate 

between windows 2-3 and 101-102, where the numbers are the trace numbers determined 

above. Process associated files concurrently. 

6. Grid the data without resampling on a uniform grid. 

7. Remove all records in each file except those containing the response from the target (use 

the Filter Dialog box). Use Save As to save files as TargetX.h1m.SU (or similar) and 

associated files. 

8. With TargetX.h1m.SU and associated files opened in plan view (this interpolates 

between GPS marks without resampling), select Process, Estimate UXO Parameters 

72
from the menu. Several initial models may need to be run to obtain a valid solution. 

After inversion, text listing UXO parameter estimates can be copied from the progress 

dialog box. 

GP Wave Utilities Data Processing Procedures 

The various procedures offered by GP Wave Utilities can be combined in many ways to solve a 

myriad of signal processing problems. Two examples are given here. For the first, the transfer function 

of a signal processing electronics module is determined. For the second, simulated data and measured 

data are synchronized and compared. 

Determining the Transfer Function 

In this processing example, the transfer function of an amplifying analog to digital (A-D) 

converter module is determined. 

1. Inject a signal from a signal generator to the input terminals of the A-D module. Record the 

signal at the input terminals with a digitizing oscilloscope. Most oscilloscopes record data in a 

text file format. Also, record the digital output of the module. 

2. Import the data recorded with the oscilloscope using the File, Import, ASCII menu item. Create 

a new document (File, New) and import the digitized output of the module. 

3. Switch to the oscilloscope data, and resample the data at the sample rate and time interval used 

when recording the module output. 

4. Switch to the module data, then use the Process, Wave Plot, Deconvolve menu item to 

deconvolve the oscilloscope data from the module output data. Next, right click on the resulting 

plot and select the frequency-domain representation. The plot now shows the frequency-domain 

transfer function of the amplifier/A-D module. 

73
Synchronize and Compare Experimental Data with Simulated Data 

1. Import the simulated data. For example, use the File, Import, GPRMax, Fields menu item. 

Import the experimental data (e.g., that from an oscilloscope) using the File, Import, ASCII 

menu item. 

2. Switch to the oscilloscope data, and resample the data at the sample rate and time interval used 

when creating the simulated data. 

3. Use the Process, Wave Plot, Compare Waves menu item to determine the time shift to apply to 

the oscilloscope data to synchronize it with the simulated data. 

4. Use the Process, Wave Plot, Filters and Operators menu item to shift the oscilloscope data by 

the amount of time determined from the previous step. 

5. Now the data are synchronized and can be visually compared. Right click on the plot and select 

the simulated data as an auxiliary trace. Now both traces will be plotted on the same graph. Use 

the Process, Wave Plot, Compare Waves menu item to determine the RMS error between the 

waveforms and the amplitude scale factor between them. 

74
References Cited 

Birchenhall, C.R., 1993, A draft guide to MatClass, a matrix class for C++, version 1.0d: Dept. of 

Econometrics and Social Statistics, U. Manchester, UK, ftp://ftp.mcc.ac.uk/pub/matclass/pc. 

Cohen, J.K., and Stockwell, J.W., Jr., 2003, CWP/SU: Seismic Unix release 37: a free package for 

seismic research and processing: Center for Wave Phenomena, Colorado School of Mines, Golden, 

Colorado, 141 p. 

Dhaene, T., Martens, L., and Zutter, D., 1994, Extended Bennia-Riad criterion for iterative frequencydomain deconvolution: IEEE Trans. Instrumentation and Measurement, v. 42, no. 2, p. 176-180. 

Gazdag, J., 1978, Wave equation migration with the phase shift method: Geophysics, v. 43, p. 1342-

1351. 

Giannopoulos, A., 2004, GPRMax – A ground penetrating radar simulation tool (v. 1.5): available 

online at http://www.gprmax.org. 

Gill, P.E., Murry, W., and Wright, M.H., 1986, Practical optimization: Amsterdam, Elsevier, 401 p. 

Jackson, J.D., 1999, Classical electrodynamics (3rd ed.): Hoboken , NJ, Wiley and Sons, 808 p. 

Lucius, J.E., and Powers, M.H., 2002, GPR data processing computer software for the PC: U.S. 

Geological Survey Open-File Report 02-166, 98 p. 

Olhoeft, G.R., 1985, Electrical properties of rocks and minerals, short course notes: available online at 

http://mars.mines.edu/pub/85ShortCourseNotes.pdf, 218 p. 

Olhoeft, G.R., 1994, Geophysical observations of geological, hydrological, and geochemical 

heterogeneity, in Proceedings of the Symposium on the Application of Geophysics to Engineering and 

Environmental Problems: Boston, MA, p. 129-141. 

75
Oden, C.P., 2006, Calibration and data processing techniques for ground penetrating radar systems with 

applications in dispersive ground: Dept. of Geophysics, Colorado School of Mines, Golden, Colo., 

PhD Dissertation, 249 p. 

Oden, C.P., Powers, M.H., Wright, D.L., and Olhoeft, G.R., 2006, Improving GPR image resolution in 

lossy ground using dispersive migration, IEEE Trans. Geoscience and Remote Sensing, submitted. 

Press, W.H., Teukolsky, S.A., Vetterling, W.T., and Flannery, B.P., 1992, Numerical recipes in C – the 

art of scientific computing (2nd ed.): New York, Cambridge University Press, 994 p. 

Powers, M. H., 1995, Dispersive Ground Penetrating Radar Modeling in 2D: Department of 

Geophysics, Colorado School of Mines, Golden, Colo., PhD Dissertation, 198 p. 

Powers, M.H., and Oden, C.P., 2004, Migration of dispersive GPR data, in Proceedings of the 10th Int'l. 

Conf. on Ground Penetrating Radar, Delft, The Netherlands, p. 333-336. 

Russel, J., 2004, Inno Setup Compiler (v. 5.0.6): available online at http://www.jrsoftware.org. 

Webring, M.W., 1998, Unpublished program based on the projection drivers from the ODDF subroutine 

package: Webring, M.W., 1998, ODDF: A file I/O subroutine package implementing NASA PDS 

data description and USGS map projections, version 1.6: U.S. Geological Survey Open-File Report 

98-765. 

Wright, D.L., and Chew, W.C., 2000, Enhancements to and characterization of the very early time 

electromagnetic (VETEM) prototype instrument and applications to shallow subsurface imaging at 

sites in the DOE complex, final report: http://www.osti.gov/em52/final_reports/60162.pdf. 

Wright, D. L., Oden, C. P., Powers, M. H., Moulton, C. W., Hutton, S. R., Kibler, J. D., Olhoeft, G. R., 

and Woodruff, W. F., 2005, A ground penetrating radar system for high loss environments: Proc. 

Symposium on the Application of Geophysics to Engineering and Environmental Problems, April 3-7, 

Atlanta, GA, p. 1-10. 

76
Wright, D.L., Moulton, C.W., Asch, T.H., Brown, P.J., Hutton, S.R., Nabighian, M.N., Li, Y., 2006, 

ALLTEM for UXO applications – First field tests: Proc. Symposium on the Application of 

Geophysics to Engineering and Environmental Problems, April 2-6, Seattle, WA, p. 1761-1775. 

Zhdanov, M.S., 2002, Geophysical inverse theory and regularization problems, Amsterdam, Elsevier, 

628 p. 

77
Appendix A 

The ALLTEM UXO forward operator approximates the induced field response of a subsurface 

UXO. The forward model includes both the induced magneto-static response, and the electro-dynamic 

response from induced transient eddy currents at a fixed instant in time. The magneto-static UXO 

response is modeled as three orthogonal magnetic dipoles, and the electro-dynamic response is modeled 

as three orthogonal exponentially-decaying magnetic dipoles. It is assumed that the target and the 

ALLTEM cart are in a non-magnetic, non-conducting whole space. The modeled magneto-static 

induction (T) at a receiver coil Bs(r) is calculated using 

3Rˆ (ms ⋅ Rˆ ) − ms 

R3 

t ⎡ ⎤

B (r) k (A1) ⎢

⎣ ⎥⎦

=
 ⋅
 s 

t where k is a calibration matrix, r is the location of the receiver, r′ is the UXO location, R = r − r′ , 

R = R , Rˆ = R / R , and ms is the static induced dipole moment (A-m2). The static induced dipole 

moment is given by 

ms,1mˆ s,1 ms,1mˆ s,2 ms,1mˆ s,3 ⎡ ⎤ ⎢

⎢

⎢

⎣

⎥

⎥

⎥

⎦

mˆ mˆ mˆ ⋅
Hp m =
 m m m (r′) (A2) s s,2 s ,1 s ,2 s ,2 s ,2 s,3

ms,3mˆ s,1 ms,3mˆ s,2 ms,3mˆ s,3 

where Hp (r′) is the primary magnetic field (A/m), and the matrix is the polarizability tensor (m3

) and the 

three induced magnetic moments are related by 

mˆ 2 = mˆ 1 × zˆ and mˆ 3 = mˆ 1 × mˆ 2 . (A3) 

The primary field Hp (r′) at the UXO location is calculated using the Biot-Savart Law (Jackson, 1999) 

for the 1-meter-square loop transmitting coils. Similar relationships are used to model the electrodynamic induction at a receiver coil, 

78
3Rˆ (mt (t) ⋅ Rˆ ) − mt t ⎡ (t)⎤

Bt (r) k ⎢ , (A4)

⎣ ⎥⎦

=
 ⋅
 R3 

mt ,1mˆ t ,1e−t /τ 1 mt ,1mˆ t ,2e−t /τ 1 mt ,1mˆ t ,3e−t /τ 1 ⎡ ⎤

⎥

⎥

⎥

⎦

mˆ ,1 

−t /τ 2 m mˆ −t /τ 2 m mˆ −t /τ m (t) =
 H ( ′), (A5) m 2

mt,3mˆ t ,1e−t /τ 3 mt,3mˆ t ,2e−t /τ 3 mt,3mˆ t,3e−t /τ 3 

e e e ⋅
 r t t ,2 t t ,2 t ,2 t ,2 t ,3 p 

⎢

⎢

⎢

⎣

where τ1-τ3 are the late time decay constants. In practice, the decay constant are ignored, and the 

electro-dynamic response is calculated for a given instant in time. 

To calculate the ALLTEM response to a UXO target using a given loop polarization, the 

primary field Hp (r′) is calculated at the UXO location for the relevant transmitting loop, and the B fields 

at the relevant receiver loop locations are calculated. Each transmitter coil is modeled as a one-metersquare loop. For a given measurement, the ALLTEM uses two receiver coils to make a differential 

measurement. The B fields at each coil location are subtracted to model the differential measurement. 

Each receiver coil is approximated as an ideal small dipole located at the center of the actual coil. The 

error introduced by this approximation increases as the distance between a receiver coil and the UXO 

decreases, and as the area of the receiver coil increases. 

79
Appendix B 

This Appendix describes the steps for setting up the GP Workbench development project so that 

changes can be made to the code as needed, and new executables can be compiled. GP Workbench was 

written using Microsoft Visual C++ v. 6. This development system or a compatible development 

system must be installed to make changes and recompile GP Workbench. 

1. Make a directory for the GP Workbench code. Extract the contents of the self-extracting zip-file 

GPWorkbench1_7_1_11.exe to this directory. 

2. Open the GP Wave Utils directory, and then open the MatClass directory. Double click on the 

MatClass.dsw file to bring up the MatClass project. Select Build, Set Active Configuration from 

the menu. Select the debug configuration. Select Build, Rebuild All to build a debug version of 

the project. Repeat for the release version of the project. Copy the file \FDTDProc – 

MatClass\MatClass\Release\MatClass.lib to the \FDTDProc – MatClass\Release directory and 

to the \GP Workbench\Release directory. Copy the file \FDTDProc – 

MatClass\MatClass\Debug\MatClass.lib to the \FDTDProc – MatClass\Debug directory and to 

the \GP Workbench\Debug directory. 

3. Open the GP Wave Utils directory, and then open the FDTDProc - MatClass directory. Double 

click on the FDTDProc.dsw file to bring up the GP Wave Utilities project. Select Build, Set 

Active Configuration from the menu. Select the debug configuration. Select Build, Rebuild All 

to build a debug version of the project. Repeat for the release version of the project. Look for a 

new FDTDProc.exe in the debug and release directories. Copy the release version of the 

executable to the program files directory where GP Workbench was installed (typically 

C:\Program Files\GP Workbench\). 

80
4. Open the GP Workbench directory. Double click on the Radar.dsw file to bring up the GP 

Workbench project. Select Build, Set Active Configuration from the menu. Select the debug 

configuration. Select Build, Rebuild All to build a debug version of the project. Repeat for the 

release version of the project. Look for a new GPWorkbench.exe in the debug and release 

directories. Copy the release version of the executable to the program files directory where GP 

Workbench was installed (typically C:\Program Files\GP Workbench\). Note that when running 

GP Workbench from the \GP Workbench\Release\ or the \GP Workbench\Debug\ directories, 

some GP Workbench functions may not operate correctly because the program may need certain 

files that exist in the C:\Program Files\GP Workbench\ directory but not in the development 

directories. 

5. The Inno Setup Compiler (Russel, 2004) can be used to make a single file install program for 

GP Workbench. Make the necessary changes to the file GP 

Workbench\GPWorkbenchInstaller.iss to make the install program. Consult the documentation 

for the Inno Setup Compiler for more information. 

81