Utility-Locating Technology Development Using Multisensor Platforms (2014)

Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

12 C h a p t e r 3 Multichannel Gpr System Status, Findings, and applications UIT’s strategy for part of this project was to build on the multi- channel ground-penetrating radar (GPR) system, Terra Vision II, already constructed for enhanced utility detection. TerraVision is GSSI’s first truly one-pass, 3-D digital geophysical mapping system, and UIT has been using TerraVision II for util- ity mapping since 2009. TerraVision II collects a 5.12-ft wide rib- bon of data at a range of speeds and resolutions that adapt to the desired application. The internal survey wheel keeps track of data spacing while the GNSS or RTS keeps track of the relative sensor positions in project coordinate space. TerraVision II is designed to collect 1,200 ft of data at 3–5 mph, delivering the full resolution of one scan per inch, which allows for detection of smaller utility targets. From these data, 3-D images of the subsurface, along with 2-D map views may be constructed. Once the data set is collected or a project site has been cov- ered, the data files are transferred to the main system computer and archived to a storage media. The data are then transferred to processing software for analysis. The results of the analysis may be presented, edited, and saved in a 3-D interface that can be output to a CADD file. Appendices A and B are standard operating procedures (SOPs) for multichannel GPR digital geophysical mapping data collection and data processing. This document provides general procedures followed by UIT for geophysical data collection during subsurface utility field investigations with TerraVision II. Similar multichannel GPR systems are commercially available, but UIT has sustained efforts in making TerraVision II the state of the art through hardware improvements and through the development and maintenance of specialized proprietary data acquisition and data processing software. GPR is an integral part of a multisensor utility detection system. It is an extremely effective tool for mapping utilities 12-ft deep or more in favorable soils, specifically those that contain a high proportion of sand. With increasing clay con- tent, GPR’s penetration depth is decreased because higher soil conductivity values attenuate the imparted radar energy. Because the hardware behind the TerraVision II is well estab- lished, the R01B project focused on an approach for dealing with the substantial amount of data that it produces rather than on improvements to the hardware configuration. Processing and analysis of multichannel GPR data are complex and time-consuming exercises. The R01B project team worked on methods to improve the analysis process through the development of software algorithms designed for GPR feature extraction and user-defined advanced pro- cessing parameters. Several algorithms were explored, devel- oped, and tested for use in the system, providing a “toolkit” of methods for the user. The software improvements developed under the R01B project and implemented within SPADE are described in the next section. Improved applications for Geophysical Data analysis Software Software is the key to integrated systems. UIT has developed three software packages for operating, processing, and inter- preting multiple-sensor data. Data Acquisition Shell (DAS) is the software used to acquire geophysical data. DAS collects all of the field data from the various sensors while allowing the operator to track the system’s location and tracks over a site. Data Processing Engine (DPE) performs preprocessing for interpretation of the data, correctly applies position data to each data point, and prepares the data set to go in the visu- alization and interpretation software. The Semiautomated Process and Detect (SPADE) software package accepts final interpretation-ready data. The R01B project team and Sagentia, Ltd., have recently enhanced SPADE by making improvements to the semi- automated functions, making it easier for interpreters to Findings and Applications

13 manipulate GPR data for analysis. Other enhancements to SPADE include functions for 3-D migration of GPR data and automated target recognition routines. SPADE was designed to accept any kind of 2-D or 3-D spatial data into the workspace. This feature allows the interpreter to input CADD files of utility data, site photos, field notes, or any other spatially referenced data along with the geophysical data. UIT demonstrated this process to So-Deep, Inc., a subsurface engineering firm, during in-service testing. This process can be considered first-order data fusion and involves getting all data plotted together at the same scale for quick and accurate comparison in 3-D. For in-service testing, data from the R01B prototypes, the associ- ated interpretation layers, and all other relevant geo-referenced data were pulled together and compared against existing SUE records in one graphical and interactive window of SPADE. The foundational concepts behind SPADE are directly related to the SHRP 2 objectives focused on enhancing utility detection and location using multiple sensors and leveraging the presence of known utilities to enhance detection. Those concepts are • Integration of data from multiple sensors and existing map data; • Use of sensor-specific algorithms to process raw data and detect features in a largely automated way; and • Use of data fusion at the feature level, plus expert human interpretation to produce the best quality output from the available data. Historically, UIT’s workflow for picking features from 3-D GPR images involved working systematically through 2-D slices, a painstaking process that relies on the skill of the interpreter. In Phase 1 of Project R01B, Sagentia devel- oped algorithms for identifying extended structures, referred to as segments, in the GPR imagery and presenting them to the user in 3-D views. These segments have the potential to greatly reduce the time taken to reduce an image and to allow the user to focus on interpreting fea- tures in the context of other data (aerial photos, TDEMI images, etc.). The R01B research approach was built on the understanding of what works best in practice. The resultant enhancements to SPADE algorithms developed under this research included • Develop an initial set of segmentation algorithms and create prototype in MATLAB; • Test algorithms on archived UIT GPR data; • Create basic user interface for visualization and manipula- tion of segments and perform validation trials; • Test segmentation algorithms: 44 Peak finding, and 44 Peak joining; • Refine algorithms; • Review raw image quality; • View images with polarization effects; • Contribute to GPR migration; and • Port algorithms to VTK SPADE from MATLAB. Figure 3.1 illustrates how UIT’s approach to feature extrac- tion addresses the main challenges to the requirements when working with multisensor data. This principally relates to Good target coverage Wide range of target types, depths, soils High accuracy Low level of false positives (false targets) and false negatives (missed targets) Good speed Near-real-time detection Requirements Signal variability Signals from a target vary with depth and soil type Clutter GPR and seismic data are very complex Large data volume Clever but slow algorithms are of little use Challenges Data Segments Fast algorithms detect extended structures in the data Segments Features Target- and sensor-specific algorithms efficiently identify candidate features Fusion and Human Interpretation Human users review “most likely” candidate features, taking into account results from multiple sensors and knowledge of the environment Approach Figure 3.1. Illustration of UIT’s approach to algorithms.

14 GPR and potentially to seismic data. The system requires an optimal combination of machine and human intelligence. UIT’s approach involved the following concepts: 1. GPR data can be divided into substantial segments from adjacent data points with similar signal properties. The chal- lenge is in finding which specific signal properties are the results of individual utility target sources when the GPR data sets comprise numerous smaller segments, essentially due to clutter or noise. Furthermore, GPR signals vary with depth and soil type, which increases the challenge in matching adjacent data points of similar properties. 2. Once segments are extracted from the GPR data they can be joined to produce candidate features. Features are ana- lyzed and scored so that the user can be presented with a prioritized list of features for consideration. The algo- rithms do not determine which features are valid but assist the user in systematically evaluating the whole data set. 3. A human user can assess all of the features and segments identified by the software and decide which features cor- respond to real targets by making subtle interpretations of visual data, interpreting disparate types of data, and exploiting contextual awareness. This user’s task is signifi- cantly simplified because the software presents segments and candidate features in order of priority. Segments User Interface This section contains a brief introduction and explanation of the user guide for the segments user interface developed by subcontractor Sagentia, Ltd. The key functions provided by the user interface are as follows: • A single data object contains all segments, irrespective of their source. This allows segments from multiple images to be viewed together or combined. • The segments can be sorted by various attributes and orga- nized in folders, to help the user develop a workflow for identifying the segments that correspond to features of interest. • The segments can be shown in a 3-D view of SPADE, along with other types of image data. • Segments can be merged to form larger segments or split into smaller pieces. Points can also be picked from seg- ments. (Note: Two SPADE algorithms were developed for creating segments from GPR data: correlation segments and threshold segments. Both algorithms work by finding a small set of peaks in each GPR scan and then connecting neighboring peaks in adjacent scans into segments.) • The user interface for segmentation in SPADE allows the user to create segments from GPR data by selecting an impulse response width (samples), the minimum segment size (pixels), and the maximum number of segments (0 for no limit). The user may also view and edit segments. • The “spdSegments” object contains all the data for the seg- ments, irrespective of how they were created. (This is simi- lar to the “spdFeatures” object). Within spdSegments, all segments are stored in folders. The folders can be used to sort segments into different categories. • A dialog box enables other SPADE controls to be used at the same time. Segments functionalities include plot con- trols, folder controls, sort controls, segments listing, and manipulation controls used to select, merge, split, or delete segments. A 3-D picking function also enables picks to be made from segments. • Segments can be saved as Drawing eXchange Format (DXF) files; however, the export process can be slow because seg- ments can contain many thousands of polygons, and DXF is a verbose format. • Segments may also be imported into existing or user-created folders in SPADE, and their display properties can be adjusted by the user. Segmentation Algorithm Testing—Peak Finding Much of subcontractor Sagentia, Ltd.’s Phase 1 work focused on the segmentation algorithms testing. The objective of the peak finding algorithm was to identify all of the peaks within the scan that might belong to features. Using this algorithm, it is acceptable to pick an excess of peaks because only those peaks that are adjacent to similar peaks will be used to form segments. However, some simple criteria are required to avoid picking peaks that would confuse the peak joining algorithm. A common approach to identifying peaks that correspond to real features is to use a “spiking filter.” This approach attempts to fit a scan to a set of spikes convolved with the impulse response of the system. Various spiking fil- ters were tested, but with limited success. Wiener filtering attempts to exploit the different spectral characteristics of signal and noise. However, even noise and clutter in GPR scans tends to be bandwidth limited, similar to the signal. Orthogonal matching pursuit uses recently developed meth- ods in compressive sensing; and although its performance was reasonable, the algorithm ran very slowly. Another limi- tation of these algorithms was that they prioritize large- amplitude peaks over small-amplitude peaks, and the intent of segmentation algorithms was to prioritize peaks based on their connectivity to adjacent peaks, irrespective of their amplitude. Two simple algorithms were selected for detailed testing. In the first, the GPR scan is initially low-pass filtered to enhance signals in the bandwidth of the system and attenuate noise. All peaks are then selected, except for those closer to their

15 neighbors than the width of the impulse response, which were assumed to be residual noise. This algorithm method was termed correlation peaks. In the second, all peaks with amplitudes above a specified threshold were picked, except for those closer to their neighbors than the width of the impulse response, which were assumed to be residual noise. This provided the equivalent of an “isosurface” function, which is common in other image processing packages. This algorithm method was termed threshold peaks. Both algo- rithms were fast and they avoided picking noisy peaks, which was the most important criteria for the peak finding algo- rithm. The correlation peaks method, in particular, facilitated picking small-amplitude peaks. The correlation peaks algorithm involves three steps: 1. Filtering the scan with a broadband, low-pass filter; 2. Picking one peak from each assemblage (or train) of local extrema; and 3. Converting peak heights to impulse amplitudes. Step 1 was achieved by convolving the data with a narrow Gaussian, with half the width of the Ricker impulse response of the antenna. It conditions the scan by reducing noise, and it turns saturated parts of the scan back into proper local extrema. In Step 2, illustrated in Figure 3.2, the purpose was to ignore small local minima and maxima that result from noise. Rather than apply an amplitude threshold, which would risk rejecting small but genuine peaks, it uses a spatial criterion. If two peaks are closer together than the width of the Ricker impulse response, then they are unlikely to be due to separate peaks. The algorithm gathers assemblages of local minima and maxima that are closely spaced. If an assemblage has an odd number of extrema—as shown in Figure 3.2a— then the train is ripple on a genuine local maximum or mini- mum. The most extreme peak within the train is taken to be the real peak, and the others peaks are ignored. If a train has an even number of extrema as shown in Figure 3.2b, then the train is ripple on a monotonic section of signal and can be ignored altogether. After Step 2, suppose that n peaks have been identified, with the ith peak at position zi and having amplitude Ai. Step 3 is a deconvolution step. It assumes that the original set of impulses had the same locations but had amplitudes Ci. Then 1 A I z z Ci i j j j n∑ ( )= − = where I(z) is the impulse response of the antenna. Cj is the final impulse coefficient of the sequence. This linear set of equations can be inverted to give the original impulse coefficients, Ci. Step 3 is optional and omitting it also makes the peak find- ing step run more quickly. The threshold peaks algorithm involves a simple modifica- tion to the correlation peaks algorithms. Rather than starting with all local extrema in the GPR scan, it starts with only those that exceed a threshold, T. The data analyst can specify whether to consider only local maxima (with amplitude >T), local minima (with amplitude <-T), or both. This provides a crude method for ignoring small peaks due to noise. Segmentation Algorithm Testing—Peak Joining The peak finding algorithm developed by Sagentia in Phase 1 reduces each scan to a set of peaks, each with a certain depth location and amplitude. The algorithm connects peaks in adjacent scans to form segments. The peak joining algorithm connects peaks in adjacent scans to form segments. Suppose two peaks occur at depths z1, z2 and have amplitudes A1, A2. These are only considered adjacent if z1 - z2 is not large compared with the width of the antenna’s impulse, and A1, A2 have the same sign. These are moderately weak criteria and enable segments to be followed from scan to scan even if con- fronted by large changes in signal amplitude. For each peak in the first scan there may be two or three peaks in the second scan that meet these criteria. The connection is then made between the peaks that are closest together (i.e., for which z1 - z2 is smallest). Some additional topological criteria are applied that govern how segments grow and merge. There are two versions of this part of the algorithm: edge segments algorithm and square segments algorithm. The edge algorithm looks for neighboring peaks within the correlation distance and then joins peaks on a one-to-one basis, preferring to join peaks that are closer together given certain conditions: the adjacent peaks cannot change polarity, nor can a peak be joined to a segment if it would cause the segment to overlap in the time direction. The squares algorithm attempts to link peaks into squares (a mutu- ally connected set of four peaks). Out of potential groupings of points, squares are chosen such that the curvatures of both peak height and signal values are minimized and adjacent squares are joined to form segment sheets. Figure 3.2. Step 2 in the correlation peaks algorithm: (a) odd number of extrema, with the red dot as the most extreme, and (b) even number of extrema.

16 Figure 3.3a. Step 1 of four in peak joining. Figure 3.3b. Step 2 of four in peak joining. Figure 3.3c. Step 3 of four in peak joining. Figure 3.3d. Step 4 of four in peak joining. GPR images are typically too large to be read into memory, so the process of joining peaks into segments was constrained to work with one channel of the image at a time. In SPADE, each channel is a contiguous part of the image file, which enables it to be quickly read from disc. Figure 3.3 illustrates four steps that are used in the process of building segments and the data structures that are used. Figure 3.3a shows the state of the system having just found peaks within the current channel. These are shown as green circles. The peaks found in the previous channel are shown as red and blue circles. The connections between adjacent peaks in the previous channel are shown as red arrows. Every segment in the previous chan- nel is represented by a “tree” of peaks in which the blue circles represent the “root” of each tree. Two peaks are therefore part of the same segment if their connections can be followed to the same root peak. These trees are the first data structure used to keep track of segments. Each node peak also has the job of remembering the peaks from all the channels to date that belong to its segment. This information is stored in a “map” (a rectangular array over a range of scans and chan- nels). The maps for two segments are shown in Figure 3.3a as blue gridded rectangles. Non-zero entries in the array specify the depth location of a peak (shown as pale red circles), and zero values correspond to channel scan locations with no peak. These maps are the second data structure that is used to keep track of segments. Figure 3.3b shows the state of the system halfway through processing the current channel. Three peaks in the current channel are adjacent to peaks in the segment on the left. New connections have been made between peaks in the previous channel, in the current channel, and between the two. The tree structure for the segment has been reconfigured so that the root peak has moved into the current channel. The map for the segment has also expanded to include the peaks in the current channel. Figure 3.3c shows the state of the system just at the end of processing the current channel. Not only have the peaks on the right connected to the segment on the right, they have also connected to the peaks on the left in the current channel. This requires merging the two segments. There is now just one root peak (blue circle) and, associated with it, just one map, covering all of the peaks. Finally, Figure 3.3d shows the system about to advance to the next channel. The peaks in the previous channel are discarded; their existence is recorded in the map. Only the peaks in the current channel remain, with connections between them. If a segment has not connected to any peak in the current channel, it cannot grow any further. It is removed from the set of seg- ments being tracked and is added to the set of segments to be outputted. When the last channel has been processed, all the segments currently being tracked are transferred to the set of segments to be outputted.

17 Performance and Refinement of the Segmentation Method Segments offer clear identification of features in data, provid- ing a global view of each structure in a 3-D interface. This can simplify the task of identifying pipe reflections, soil layers, and point objects by providing a more intuitive picture of the signal. In the course of reviewing the segmentation function- ality, some improvements were made to the algorithms. It was determined that segmentation finds most of the significant features in an image; however, it can still miss faint features which are disjointed by noise and, therefore, cannot entirely replace visual inspection of 2-D slices. Analysis of segments is a time-consuming process. A sub- stantial number of segments are produced; and although there are automated ranking systems to filter out the weakest seg- ments, dozens of segments still need to be analyzed by a data analyst for a relatively small image of 300 scans by 14 channels. Furthermore, the algorithm is not completely successful in identifying all features. Examples were discovered in which smaller objects were only partially segmented by the squares algorithm, and neither the edges nor squares algorithms were able to identify all of the pipes on some sample data sets. The majority of features were segmented, in whole or in part. How- ever, the numbers of false negatives are not significantly mini- mized; thus the effectiveness was determined to be comparable to that of traditional methods of using migrated time slices for GPR analysis and interpretation. On the one hand, the edges algorithm displayed a propensity to connect features to sheets of noise, and to add nonsmooth edges to segment sheets. This made features harder to pick out and identify and increased the analysis time, forcing the user to split and merge segments until a more representative feature became clear. On the other hand, the squares algorithm, which appeared to be only slightly less successful in segmenting features present in the data, did not exhibit such a greedy approach. It provided clearly defined (and more compact) shapes, making it easy to identify features, reducing operator error. This clarity of results made the squares algorithm marginally preferred during testing. Raw Image Quality Sagentia compiled a database of images from UIT TerraVision data sets to aid in the software and algorithms testing. The database identified a set of images, which span a range of fea- tures, identifying strength, feature type, feature direction (if linear), number of features, whether features were straight, whether they crossed, and whether they exhibited polarization effects. In all the images compiled, the ground reflection was very well aligned in time across each acquisition swath. The 14 GPR antennas of the TerraVision are array-staggered in the along-channel direction, with the biggest offset occurring between Channels 7 and 8. In most of the images reviewed, the channels were very well aligned in along-channel direction, which indicated that the antennas’ offsets were being correctly compensated for. However, in rare cases, there were visible misregistrations between adjacent swaths, and it was pre- sumed that such errors could be corrected after-the-fact by adjusting the map files. Figure 3.4 shows the impulse response of each of the 14 channels in the TerraVision system. These were calculated from an image for which the autocorrelation of the scans was averaged along each channel, and the autocorrelations were converted into impulse responses by assuming a zero-phase signal. The impulse responses were very similar from channel to channel, and they fit well to a Ricker wavelet. The most imperative issue observed in raw images was the presence of ringing in some channels. “Ringing” involves per- sistent oscillations that occur in a few specific channels. This is the most serious artifact that occurs in TerraVision data; it can easily obscure features, both in map view and in vertical slices. It is unlikely that this can be addressed by image pro- cessing. In many cases, the ringing increases with depth until the data saturates. Images Viewed with Polarization Artifacts Polarization effects in TerraVision data are the occurrence of a strong attenuation of GPR signal at every other channel of the array due to the angled orientation of utility target positions in the subsurface relative to the antennas’ direc- tion during acquisition. The effect of splitting channels based on polarization (and analyzing the two sets of images separately) on the ability to identify features from time slices of data was assessed for several pipe situations: diago- nal pipes, which show polarization effects; perpendicular Figure 3.4. Impulse response for each of 14 channels, derived from raw data by an autocorrelation method.

18 (to the swath) pipes; and longitudinal pipes (in the same direction as the swath). Most pipes that cut a swath at an angle around 45° exhibit polarization effects. Viewed from a time slice, the pipe appears banded, reducing its visibility. Splitting the image into sepa- rate polarizations produces one image in which diagonal pipes are clear and one in which such pipes have a highly attenuated response. (See Figure 3.5.) Splitting an image into two sepa- rate polarizations does, therefore, appear to aid pipe detection when using time slices; however, both polarizations must be checked so that pipes in orthogonal directions are not missed. Perpendicular pipes often show up clearly in time slices. For pipes that give strong signals, viewing only one polariza- tion can give a smoother image. However, the advantages of trying to pick out pipes with weak signals were less clear. When viewed in time slices, longitudinal pipes can be diffi- cult to distinguish from channel imbalances. Splitting polar- izations did not seem to reduce the visibility of longitudinal pipes in time slices, and in some cases resulted in a more clearly defined target. Polarization effects are also apparent in the detection of localized objects as well with the target dis- appearing every other channel in time-slice view. Separation of the data into separate polarizations reduced the ease of identification of these point-source signals in time-slice view, which was not a problem for linear target identification. In general, it was determined that polarization-sensitive meth- ods of viewing the raw multichannel GPR data do not neces- sarily provide the data analyst with an easy method to detect and identify features, particularly faint features. Contributions to GPR Migration Point features in the ground generate hyperbolic cones in raw GPR images; linear features generate hyperbolic sheets. Migra- tion algorithms stack images along hyperboloids, concentrat- ing the energy from hyperbolic cones and sheets back to points and lines respectively. This is not a mathematically perfect inversion, but it serves to increase contrast and sharpen fea- tures. UIT has historically had limited success in getting quality migrated images with multichannel GPR data. During the SHRP 2 R01B project, Sagentia released a bug fix for SPADE that addressed this, as well as a migration preview function that makes it easier for the data analyst to set the algorithm parameters. Three algorithms using a constant soil dielectric were studied: the Witten stack, Stolt migration, and Kirchoff migration. Each algorithm is mathematically similar and displayed similar image quality during testing. The Witten algorithm is simple and intuitive. It does not correspond to any mathematical model of wave propagation, unlike Stolt and Kirchoff migrations. However, it is also slow. The Witten 2-D code explicitly creates the weight image (as illustrated in Figure 3.6) and applies it to the whole raw image, ignoring the fact that the vast majority of weight values are zero. Stolt and Kirchoff migrations are based on a mathematical model of wave propagation in a medium with a constant speed of light (i.e., constant soil dielectric). Nonetheless, they can both be described in terms of image weights, as with the Witten algorithm. The Stolt algorithm decomposes the raw image into plane waves. The migration step involves altering the wave vector for each mode and recombining the modes to get the migrated image. The conversion to Fourier modes and back again is accomplished by fast Fourier transform (FFT), Figure 3.5. GPR time slice from diagonal pipe: (left) all channels, with polarization effects visible; (middle) Channel B, poor pipe visibility; (right) Channel A, clear pipe definition. Figure 3.6. Image weights for Witten migration.

19 which makes the algorithm relatively fast. The Kirchoff algo- rithm is a reverse time migration in which the raw image rep- resents energy collected at the surface as a function of time, and this is propagated backward to find the distribution of energy in the ground at time zero. The weight image of the Kirchoff migration (Figure 3.7) appeared less noisy than for the Stolt migration, resulting in a cleaner migrated image; but Kirchoff migration ran slower than Stolt migration because there is no FFT-style shortcut for its implementation. The effectiveness of the Stolt migration in SPADE was studied by migrating a small section of GPR data and tuning the “z-scale factor” until the migration was optimized. This factor was then used to migrate larger sets of data from the same site. These studies resulted in a new version of SPADE (delivered to UIT in March 2012) that included bug fixes for the migration algorithm contained. GPR Migration and Time Slices Comparisons of GPR time slices were made with and with- out migration applied. Linear features associated with sub- surface utilities were identified in the migrated GPR slice images and occasionally displayed slightly sharper definition than the unmigrated slices. Figure 3.8 shows a section of GPR data collected along three acquisition passes on a road with two subsurface pipes: one straight, one bent. The center swath is an unprocessed time slice; the outer swaths show the results of migration. It should be noted that although the migrated time slices offer tighter target resolution for several linear targets, some of the weak detection responses from subsurface utility sources were clearer without migration applied. Near the surface, features in migrated images stand out over a low noise background; but with deeper images, the noise proved more prominent, even for the migrated image. Weaker signals can align in deep time slices and be mistaken for linear features and potentially targets of interest. For this reason, it is important that the data analysts view both the migrated and unmigrated time slices images during data interpretation. Near-surface targets can be more precisely depicted on migrated GPR time slices, but deeper migrated time slices can display false targets of aligned noise. GPR Migration and Polarization Migration performance with polarization effects was also studied. In theory, the migration algorithm should be unaf- fected by polarization effects. Summing along a hyperbola in which every second channel has low signal strength will still result in a concentration of energy at the peak, albeit with a reduced signal strength compared with that from an image without polarization effects. To demonstrate this point, a series of synthetic images were constructed for analysis. Each image contained a simulated hyperbolic response from a pipe. Hyperbolic responses were constructed first with all channels responding and second with alternate channels highly attenuated (simulated polarization effects). Both hyperbolae images displayed a bright point after migration, although the image without polarization effects showed a sharper and stronger response, as expected. Indeed, migrated images of diagonal pipes, viewed in time slices, also showed pipes clearly, even when simulated polarization effects were present. When used on real project GPR data, the migration per- formed well in coping with images showing polarization effects. Pipes running diagonally across a swath (which do not always show polarization effects) were clearly imaged with no striping observed in the migrated image (Fig- ure 3.9). The contrast between noise and migrated pipe was not as good as for straight pipes without polarization effects (Figure 3.10). Figure 3.7. Weights applied to the raw image under Kirchoff migration (to 1 pixel at apex of hyperbola). Figure 3.8. Migrated-unmigrated time-slice comparison.

20 For images which show polarization effects, splitting the polarizations would be expected to increase the performance of the migration process, by removing channels with negligible response to maximize the signal in the mapped hyperbolic sheet. Sagentia used a three-step process to split channels: (1) separate the channels into A and B polarizations, (2) inter- polate between the remaining channels in each image, and (3) perform migration on both images separately. Using syn- thetic images, it was seen that migrating individual polarizations does offer some improvement in the strength and continuity of signal, over migrating both polarizations together. The migrated result from the polarization showing the strongest response was marginally stronger than the result without splitting polariza- tions. However, the added clarity is probably not worth the time cost of doubling the number of images to view. GPR Migration Assessment GPR migration within SPADE has the potential to be a power- ful tool for feature identification from GPR data. Migration col- lapses the pipe signal to a localized line, and this provides a quicker (and potentially more accurate) estimation of the depth of an object than achievable from the raw signal response. Fur- thermore, the localized signal could allow tracing algorithms to track the extent of a feature, after one point is identified. A comparison of migrated and unmigrated images reveals that, in general, pipes which show up in migrated images are also visible in unmigrated time slices. Migrated images tend to have lower noise, which can give sharper pipe definition. While some pipes are clearer after migration, others, par- ticularly those with weak signals, are not. The visibility of the thin signal line in migrated features, particularly when large swaths are viewed as a whole, is often poor when compared with the broad “bloom” of pipe features in unmi- grated time slices. Successful migration depends on an accurate choice for the z-scale factor (which takes account of the signal speed in local soil type). A number of methods are provided in SPADE for testing this scaling to include hyperbola fitting and migration parameter selection preview options. The benefits of using migration over assessing unmigrated time slices are not read- ily perceived; the improvements to workflow consist of a tran- sition to assessing time slices of large images (over vertical images of small subsections) and an opportunity for increased automation in feature tracing. high-Frequency Seismic Imaging proof-of-Concept prototype Status and Findings The seismic component of the R01B research was a first of its kind endeavor riddled with unknowns, highly technical con- cepts, and limitations posed to feasible configurations and computational electronics. So the project team focused on proof of concepts (specifically, shallow in situ soil seismic properties) in support of research objectives to design and con- struct the proposed high-frequency seismic imaging proto- type system. This moderate approach to development was due to several factors, including the slow process of finalizing project contracts and subcontracts, unforeseeable manufac- turing and vendor equipment delivery delays, limited experi- ence with and knowledge about new system electronics, and the lack of general support documentation. The seismic ele- ment of SHRP 2 R01B research work was performed by Owen Engineering Services with support from Psi-G, LLC and Bay Geophysical, Inc. a.) b.) Figure 3.9. GPR time slices of diagonal pipe: (a) unmodified image, (b) migrated image in which pipe is most visible in negative component of reflected signal. a.) b.) Figure 3.10. GPR time slices of perpendicular pipe (showing one perpendicular pipe and one pipe at an angle, with polarization effects): (a) unmodified image, (b) migrated image, with migration scaling factor tuned for perpendicular pipe.

21 Knowledge of seismic wave propagation in near-surface soil materials is critical to the design and optimization of an effec- tive seismic method for detecting and mapping underground utilities. Because of regional variability in soils and localized variations in utility backfills, in situ field measurements are essential to gaining a full understanding of underground util- ity environments for realistic high-resolution seismic system performance. Therefore, actual field testing and site-by-site characterizations of the various utility host media are manda- tory for accurate definition of the seismic properties. A regional survey of soil seismic properties was conducted emphasizing measurement of horizontally polarized shear (SH) wave propagation parameters along with similar measure- ments of compressional (P) waves. The primary seismic param- eters were the frequency-dependent viscoelastic attenuation versus depth and the propagation velocity versus depth—both measured along a common down-transmission vertical path. The methodology used in this task was specialized to high- resolution seismic attenuation and velocity measurements using a push-type soil probe system equipped with a three- component seismic CPT (cone-penetration test) toolhead attachment. The information gathered as a result of this test- ing aided in defining the general ranges of velocity and attenu- ation in near-surface materials and served to indicate the possibility of following through on the technical approach to develop an effective high-resolution seismic system for under- ground utility detection and mapping. Soil Seismic Properties Background Only a limited amount of quantitative information was avail- able on in situ seismic properties of shallow soils leading up to the SHRP 2 R01B research work. A noticeable limitation in nearly all of the documentation on seismic velocities in shallow soils was the relatively low frequency range for which practical measurements have been reported. Underground utility envi- ronments require defined seismic soils properties up to about 1,600 Hz. The available background of reported work at low frequencies has justifiably neglected the frequency dependence of velocity caused by dispersion effects inherent in viscoelastic attenuation; new measurements would require that velocity dispersion be accounted for in the soil seismic properties. Another characteristic of seismic propagation in shallow sur- face soils is the velocity gradient versus depth caused by com- paction in granular materials. Other factors such as grain size and shape and grain cementation also affect the elastic moduli, giving rise to differences between natural undisturbed soils and trench backfill soils. The soil modulus factors affect the P-wave and shear wave (S-wave) velocities; in general, soil materials having high mod- uli have higher seismic velocities and lower internal friction and, therefore, have lower attenuation (and higher attenuation coefficient, Q). Based on limited data researched, typical seis- mic interval velocities within the top 20–30 ft increase by about a factor of two relative to the velocity at the surface, depending on the soil type and compaction. Sedimentary layer variations in soil type may occur, but lateral variations in velocity are not significant at given localized test sites unless distinct areal changes in soil type are present. A seismic medium of this type is referred to as a transverse isotropic medium, characterized by only a vertical velocity gradient. In regard to R01B seismic soil properties testing, principal pertinent information was derived from documentation research and is revealed in a paper by Lew and Campbell (1985), which is illustrated in Figure 3.11. That information is as follows: 1. S-wave velocities at all of the test sites exhibited a significant increase versus depth in the top 30 ft of the ground; 2. Compressional wave velocities at all of the test sites exhib- ited a significant increase versus depth in the top 30 ft of the ground; 3. Statistical variations in S-wave velocity among the 270 test sites (95% confidence), for depths of 2–15 ft, are typically about ±17% for undisturbed soils and about ±11% for nonengineered backfills; 4. The highest velocity gradients occur in the top 10 ft com- pared with gradients in the underlying 10–30 ft, with firm soils and backfill materials having the greatest differences between these depth zones; and 5. The P-wave velocity to S-wave ratio ranged from 1.91 to 1.97 for all of the soil types in the top 30 ft, with the slightly greater effect in the soft natural soil. Although all of the test sites reported by Lew and Campbell were located in California, their specialized engineering soil classifications are directly relevant to similar soil conditions in other geographic regions. Equally important is that a large number of test sites were used in characterizing the velocities and vertical gradients, resulting in well-represented regression- analyzed nominal in situ velocity profiles for each soil-type classification. Shear wave velocity measurements versus depth in shallow soils were of direct interest in the planned R01B soil seismic properties survey, so several contributing factors were care- fully considered during test planning. Those factors centered on the implications of vertical velocity gradients and the implications of seismic attenuation in shallow soils. All of these considerations provided a clear rationale for conduct- ing new field tests to characterize soil seismic properties rel- evant to high-resolution underground utility detection and mapping. The scarcity of prior work in measuring attenua- tion and velocity dispersion in unconsolidated shallow soils indicated that many aspects of the new measurements would

22 demand close attention to devising proper field procedures. Thus, proven equipment technology capable of assuring acquisition of high-quality field data was used. In addition, the data processing techniques were expanded to include extracting velocity profiles in the transverse isotropic soils of interest when affected by frequency-dependent dispersion as well as extracting the amplitude attenuation coefficient, or Q, of the medium as a function of frequency over a range of high seismic frequencies not previously investigated. For accurate and consistent results, these tests required a dedicated mobile field system capable of measuring and recording seismic velocity and viscoelastic attenuation versus depth for both S-waves and P-waves at depths to 15–20 ft as functions of frequency over the range 50–1,600 Hz. A primary consideration in performing this study was to determine whether the transmission of high-frequency shear (S) waves or compressional (P) waves is practical given the attenuation that is likely to occur. To that end, the project team selected three separate regional sites to conduct the seis- mic soil properties testing: Manteno, Illinois; Houston, Texas; and Manassas Park, Virginia. Testing was performed during the fall of 2011. Before actual field testing, significant efforts Source: Adapted from Lew and Campbell (1985). Figure 3.11. Shear wave and compressional wave velocity versus depth in soft, intermediate, firm natural soils, and nonengineered backfills.

23 were made in preparing and configuring both the software and hardware components of the seismic data acquisition system for a full dry-run soil properties testing procedure. During this trial testing, in Traverse City, Michigan, the team determined the proper procedures for equipment setup, data acquisition, and initial data processing. A coordination plan was also developed to establish responsibilities, schedule suf- ficient time in obtaining one-call tickets, provide accurate locations of the testing sites, and allow sufficient up-front time to the host installations. Soil Seismic Properties Preparation, Equipment, and Data Acquisition The methodology of measuring soil seismic velocities and attenuation at high frequencies was given detailed consider- ation by subcontractor Owen Engineering Services. The in situ approach was selected because only such realistic tests could provide meaningful results relevant to actual field con- ditions. This fundamental criterion determined the planned equipment and measurement methodology. The system frequency response requirements were calculated on the basis of the project goal of detecting a 3-in.-diameter pipe at a depth of 12 ft. Assuming that the detection threshold occurs when the pipe radius is one-quarter wavelength at the pipe burial depth, the incident seismic wave must have a wavelength of 2DPλ= and a frequency of 2 f v v DP = λ = where v = speed of wave, l = wavelength, f = frequency, and DP = pipe diameter. Using this relationship, the lowest incident wave frequency for which a 3-in.-diameter pipe is detectable in soils having vs = 500 ft/s is 1,000 Hz. Larger-diameter pipes are detectable at proportionately lower frequencies. To produce the desired high-frequency seismic signals required for small-diameter pipe detection, vibrator sources generating continuous linear frequency sweeps were consid- ered, but separating independent frequencies for velocity and attenuation analysis would be difficult from a data processing perspective. Therefore, it was determined that the best way to measure attenuation of signals was to use monochromatic tonebursts with center frequencies of 100, 200, 300, 400, 500, 600, 700, 800, 900, 1,000, 1,100, 1,200, 1,300, 1,400, 1,500, and 1,600 Hz. Figure 3.12 shows the design features of these tonebursts, each having a time duration of 60 ms and sepa- rated in time by 120 ms. The amplitude envelopes of the tonebursts are cosine-squared in shape to confine the fre- quency spectrum to 100-Hz bandwidth at each toneburst frequency. The seismic microvibrator, MicroVib, was selected as the source for these tests, since it is capable of generating vibra- tions at the frequencies required. This proprietary source was designed by Owen Engineering Services and constructed under license by Bay Geophysical in the late 1990s. It has been used on numerous seismic surveys for generating SH-waves and P-waves at frequencies up to about 400–500 Hz. The MicroVib was fitted with a fabricated point-source ground- coupling base plate as shown in Figure 3.13. A new power amplifier was chosen as the device to provide current to the vibrator. This device must be capable of delivering constant current at the wide range of frequencies specified for the tonebursts. Accelerometers were chosen as the transducers to use for the measurements because of their uniform frequency response over the range of frequencies required. In particular, the selected seismic-grade accelerometers have exceptionally low self-noise in the frequency range of interest. A special tri- axial accelerometer probe housing was designed by Owen Engineering Services: the accelerometers were mounted as a set of three orthogonal transducers and pushed into a bore- hole for placement at various depths below surface from the source to measure the travel time and attenuation of the transmitted tonebursts (see Figure 3.14). Finally, to record the transmitted signals and generate con- trol signals for the vibrator, it was necessary to select a data logging and control system such as the Dewetron Model DEWE-3210 with high sensitivity, high sampling frequency, a wide digital dynamic range, a signal generator, and computer control capabilities. The selected equipment included the Dewetron data acquisition system, Wilcoxon Model 731-207 accelerometer, and the AE Techron Model LVC 5050 excita- tion power amplifier. After procurement of the equipment listed above, the com- ponents were laboratory tested to verify that the various com- ponents functioned correctly. During this shop check, several problems were encountered. Serious ground-loop interfer- ence was encountered with the power amplifier—whenever power to the data acquisition system was turned on or when- ever a triggered toneburst sequence was originated by the data acquisition system and power amplifier. The cause of this interference was found to be a configuration problem in the power amplifier and was later resolved through conversa- tions with the manufacturer. Also, the Dewetron software was not set up to generate and record causal measurements. That

24 Figure 3.13. Ground coupler base late for MicroVib. Figure 3.12. Design features of tonebursts: (upper graph) first six pulses of 16-tone- burst sequence; (lower graph) 16 tonebursts, 60 ms in duration, separated at 120-ms time intervals. Figure 3.14. Housing probe and triaxial accelerometers.

25 is, the software did not allow for synchronized generation of the toneburst excitation signal and the initiation of A-to-D sampling of the accelerometer signals. After some time, Dewetron software engineers were able to generate revisions to the software which allowed the synchronized seismic mea- surements to be made. Finally, several spurious resonance responses were observed in the operation of the MicroVib during shop testing. As the vibrator generated the stepped- frequency tonebursts, the vibration amplitudes were stronger at some frequencies than at others. Following the laboratory check of the system and resolu- tion of all non-vibrator-related problems, the system was set up at a soil probing site in Traverse City, Michigan, where the entire system could be tested in a mockup of the actual field measurements. A Model 540UD Geoprobe push-technology ground probe system was used to push the downhole acceler- ometer probe into a pilot hole previously pushed or cored by the Geoprobe system. The selected location had a soil profile that was entirely sand, at least in the upper 15 ft; and the Geo- probe system was unable to push the downhole accelerometer probe deeper than approximately 8 ft. Other difficulties were encountered during trial testing and lessons learned. Of sig- nificance was that the variable amplitude output of the vibra- tor versus frequency was confirmed to exist in the downhole accelerometer signals. This required additional shop testing and recalibration of the MicroVib before performing the soils properties testing at the regional U.S. locations. The purpose of the field measurements was to determine the S-wave and P-wave velocity and attenuation (or Q) depth profiles versus frequency for different types of soil. To make the desired velocity and attenuation measurements at dis- crete frequencies, the MicroVib source was used to generate a sequence of 16 monochromatic toneburst pulses that start at 100 Hz and step up in frequency at 100-Hz intervals to 1,600 Hz at the upper frequency limit. Each toneburst was shaped to have a cosine-squared envelope to minimize the spectral over- lap at adjacent toneburst pulses when analyzing the data. The 16 tonebursts were spaced 120 ms apart to produce a total sequence duration of 1.92 s. This toneburst sequence was trig- gered by the Dewetron data acquisition system at a short time delay of 10 ms after initiating the A-to-D sampling and record- ing of the six accelerometer sensor channels. The excitation sequence was fed to the Techron power amplifier and MicroVib to produce the ground vibrations. To measure attenuation of the propagated seismic signals, a reference accelerometer probe was used to record the radiated source waveforms (see Figure 3.15). This reference probe was placed at a depth of 18 in. in a shallow borehole located about 2–3 ft from the downhole accelerometer probe. Triaxial acceler- ometer signals from these probes were recorded as the downhole probe was pushed downward at 1-ft depth intervals. Simplisti- cally, the ratio of the individual toneburst signal amplitudes measured at the downhole probe divided by those measured at the reference probe can be used to derive the net seismic wave attenuation over the distance between the probes. At each borehole location where attenuation was to be measured, both reference and downhole transducer probes were calibrated simultaneously to ensure accuracy of measured amplitudes. The downhole probe was pushed into the ground at incre- ments of 1 ft. At each increase in depth, a toneburst signal was sent to the MicroVib, and accelerometer responses were recorded from the two sensor probes. This procedure was repeated until the downhole probe was unable to advance to further depths or until a maximum depth of 16 ft was reached. The downhole accelerometer probe was not rotated during the pushing proce- dures. During each such measurement, the MicroVib was main- tained at a specific orientation: vertical for P-waves or horizontal for SH-waves, with SH-wave polarization either in line with the direction of the line between the reference and downhole probe boreholes (radial orientation) or perpendicular to the line between the reference and downhole probes (transverse orienta- tion). These data acquisition procedures were applied in three separate nearby boreholes to permit downward-push measure- ments for the three vibrator source orientations. Soil Seismic Properties Data Analysis and Results from Houston Testing Site Seismic field data collected at the three regional locations was of mixed quality, resulting in usable data from potentially only five of the 11 total number of test boreholes. The quality Figure 3.15. Placement of downhole and reference probes during field testing.

26 of the field data acquired in all of the test boreholes was rela- tively low because of apparent spurious response effects in the measurement apparatus, variations in coupling of the seismic vibrator source and sensor probes at the ground sur- face and downhole, and possible anomalous propagation effects and seismic reflection interference from subsurface soil layering. These factors produced amplitude and phase variations in the recorded seismic waveforms, making both velocity and attenuation analyses difficult, or impossible. Additionally, several of the test boreholes were at locations where the Geoprobe push system was unable to drive the downhole sensor probe below depths of about 9–10 ft and/or were in close proximity to 60-Hz power lines that caused excessive inductive interference. Velocity Analysis A velocity analysis was performed by reviewing time-series records of S-wave or P-wave toneburst sequences traveling along near-vertical slant paths from the surface vibrator source to the downhole triaxial accelerometer probe. With the vibrator oriented horizontally and transverse to the borehole axis, SH-waves were generated at the base plate point coupler and transmitted downward to be detected by the y-axis accelerometer sensor in the downhole probe. With the vibrator oriented vertically, P-waves were generated at the base plate and transmitted downward to be detected by the z-axis acceleration sensor in the downhole probe. Separate boreholes located near one another in essentially the same soil structure were used for the S-wave and P-wave measurements. At each borehole location, the vibrator source was placed at an offset distance, typically 24–36 in., away from the borehole. The propagation paths are slant distances from the center of the vibrator source at the ground surface to the sensor location in the downhole probe, as determined by the measured depth and the source offset distance. The seismic wave velocity in near-surface soils increases signifi- cantly with depth, causing diverging refraction effects and curved-ray paths in the down-going wavefront. However, since the propagation paths were relatively short and nearly vertical, the curved-ray paths could be accurately approxi- mated by the geometrical slant-path distance associated with each measurement. The field procedures used to record the slant-path signals included the use of a second triaxial accelerometer probe placed at a fixed depth of 18 in. at a measured location between the source contact point and the borehole. The purpose of this second probe is to provide a reference seis- mic signal representing the radiated source wave. Detailed analysis techniques were developed to use this reference sig- nal in the velocity analysis. Unfortunately, the spurious responses and interference effects in the recorded downhole and reference probes caused significant errors in the derived slant-path velocities. To avoid these errors, the excitation current waveform applied to the vibrator source was used in place of the reference probe signal to determine the S-wave and P-wave travel times. Mean Slant-Path Velocity Preliminary analyses showed that phase distortions caused by anomalous interference effects yielded different travel times for each of the toneburst signals traveling along a given slant-path distance. This distortion precluded using the toneburst signals to measure the frequency-dependent velocity effects as planned as part of the soil attenuation analysis and made the velocity analysis ambiguous. As an alternative, the single toneburst signal at 200 Hz was selected for use in the velocity analysis. By using this relatively low- frequency waveform, the likelihood of encountering phase errors large enough to cause cycle skipping in the travel time analysis was avoided. Interval Velocity The velocity in each soil depth interval between successive downhole probe measurement positions was derived from the experimentally determined mean slant-path velocities. This was based on the assumptions that the downhole depth stations were uniformly spaced at 1-ft intervals and the prop- agation paths were accurately approximated as straight lines from the surface source to the downhole probe. The uni- formly spaced measurement depths allowed the slant-path distances to be divided into an equal number of uniform interval slant-length segments in which the interval velocity was assumed to be constant. Attenuation Analysis Seismic wave attenuation is caused by a combination of imper- fect elastic properties and scattering by inhomogeneities in the propagation medium. The predominant attenuation factor is viscous friction, which leads to a frequency-dependent expo- nential decay in amplitude, a factor that surpasses attenuation caused by geometric spreading (as when waves propagate through the distance of a particular medium). A practical method of determining the excess attenuation is to measure the amplitude at one position along a known propagation path and to compare it with the amplitude measured at a second posi- tion, usually closer to the source, along the same propagation path. By first compensating the measured amplitudes for the distance-dependent geometric spreading loss, the ratio of the amplitudes can be used to derive the exponential parameters of the attenuation factor. For this method to be effective, the

27 amplitudes must be free of any differing anomalous variations at the two measurement locations, and the measurements must be made at several frequencies to provide sufficient data for evaluating the frequency-dependent attenuation coefficient. The vibrator excitation signal used in these measure- ments was designed to generate toneburst signals having equal amplitudes and bandwidths to facilitate the measure- ments. However, in practice, the signals transmitted to the downhole and reference probes were found to be distorted by spurious frequency response and interference effects that caused amplitude variations in the signals detected at the two sensor probes. In many cases, the detected signals were sufficiently different to cause large errors when com- paring the amplitudes of the separate toneburst signals. As a result, only certain toneburst signals were of suitable quality for attenuation analysis, but these could not be selected in advance. Instead, all of the toneburst signals were processed under the assumption that they could be of use and any extreme outlier values would be identified and removed later. This turned out to be an extensive and time- consuming process. Effective Slant-Path Q The analysis objective was to determine the effective attenua- tion at as many toneburst frequencies as practical along the slant path between the surface source and the downhole probe. For this purpose, all of the toneburst signals (100–1,600 Hz) were processed in the frequency domain using spectral magni- tudes only to obtain the effective Q of the medium along the same slant paths used in the mean slant-path velocity analysis. The effective Q is inversely proportional to attenuation per wavelength in the medium. To minimize the variations in sig- nal amplitude, the downhole and reference signal frequency spectra were reduced to their mean values over the toneburst bandwidth before being normalized for further analysis. Interval Q Interval Q values were derived from the derived effective slant-path Q values. The assumptions and methodology used to determine the interval Q values were similar to those applied in deriving the interval velocities, with the exception that the cumulative attenuation along the slant path was the product of partial attenuation factors in each interval. The par- tial attenuation factors in each interval involve both the effec- tive slant-path Q and the mean slant-path velocity and, thus, were subject to combined error effects. In general, the Q of the soil medium tends to increase with depth, as with interval velocities. The increasing compaction and shear rigidity of the soil column tends to reduce the anelastic losses and, thus, increase the Q of the soil with depth. Analysis Results The data acquired at the three regional field locations were of mixed quality, requiring special attention to data processing methodology and procedures. The following illustrations are selected examples that show the typical nature of the recorded signals, the irregular responses, and the noise interference affecting the data quality. The vibrator excitation current waveform was described and illustrated earlier in Figure 3.12. Figure 3.16 shows the detailed SH-wave signals at depths of 1–8 ft at a Houston borehole location. Each trace contains toneburst pulses 100, 200, . . . , 1,600 Hz at 120-ms time intervals. These traces show the com- bined anomalous effects of vibrator response, ground coupling at the source and at the downhole probe, possible reflections from deeper soil layers, possible guided waves along the bore- hole, 60-Hz power generator interference, and noise. Figure 3.17 shows the amplitude frequency spectrum of the toneburst trace detected at 4-ft depth at the Houston bore- hole. The detected toneburst spectrum is scaled by a factor of 2,000 for convenient comparison. The 60-Hz power generator interference and its harmonics and fractional harmonics are evident. The relatively strong toneburst at 100 Hz followed by the weak toneburst at 200 Hz is a counterintuitive effect since the wavelengths are large compared with the borehole and downhole probe dimensions. Other higher-frequency tone- burst spectra also showed strong or weak amplitudes and their spectral content was broadened or reduced in comparison with that of the source excitation current toneburst spectra. Insufficient experimental data are available to identify the reasons for these differences. Each downhole-detected toneburst could be filtered and analyzed in comparison with the vibrator excitation current and/or reference probe signals. By comparing the quality of the filtered tonebursts at the three regional field sites, the 200-Hz toneburst was selected for use as the reference signal for analyz- ing the travel time delays and propagation velocities along the near-vertical depth paths. Figure 3.18 shows the 200-Hz time- domain tonebursts overlaid for propagation time-delay analy- sis at the Houston borehole. The clarity of the overlay match is evident in this 4-ft depth case. The time-lag shift required for this alignment was 7.92 ms, corresponding to a mean slant- path velocity of 543 ft/s. Although the amplitude envelope of the downhole-detected signal was idealized to overcome much of the distortion in the detected waveform, the phase of the downhole signal was not changed. This phase function contains both the propagation time-lag phase shift and any other phase distortions caused, for example, by variations in ground coupling. As shown in Figure 3.18, the envelope stretching caused by phase distor- tion of the 200-Hz time-domain toneburst signal is signifi- cant, making time-lag matching adjustments imperfect and/ or ambiguous as to which cycle in the overlay is correct

28 Figure 3.16. Downhole-detected signals at Houston borehole. Figure 3.17. Frequency spectrum of the toneburst trace at Houston borehole. Figure 3.18. Overlay of the idealized 200-Hz downhole toneburst and the 200-Hz current excitation toneburst for Houston borehole.

29 (introducing a possible cycle-skip error). In several cases, the interference and noise in the weaker signals detected at the greater depths could not be reliably analyzed in this manner to obtain accurate overlay time lag. In other cases, overlay cycle skipping was tested to determine the reasonable bounds for the derived velocity. Amplitude comparisons of downhole and reference probe toneburst signals were necessary to determine the excess attenuation caused by anelastic absorption and scattering in the soil medium. These comparisons were not productive when using the extracted time-domain toneburst waveforms, but they were improved when their spectral amplitudes were compared in the frequency domain. The spectral envelope displayed a common form of asymmetrical and bimodal dis- tortion caused by interference effects and resulting in inaccu- rate amplitude ratios. These variations were further suppressed by using only the ratio of the mean spectral amplitudes aver- aged over the toneburst filter bandwidth. Using this formula- tion, a useful number of toneburst signals could be compared to yield approximate exponential decay parameters defining the effective Q of the soil medium. Soil Seismic Properties at Regional Testing Sites The in situ experimental measurements of propagation veloc- ity and attenuation are the descriptive soil seismic parameters of interest in estimating the detectability of underground utility pipes. The analysis methods were determined adequate to yield useful approximate values of velocity and effective Q versus depth at the three regional sites. The soil seismic prop- erties differed for the three regions and were also different for multiple test locations at each regional site. For the sake of brevity, the detailed results presented here are only from data acquired at the Houston, Texas, test site location. Because of data quality limitations, the results show only approximate values of velocities and Qs but, more clearly, show their trends versus depth. The depth limits were caused either by probe-push penetration refusal or by excessive data analysis errors. Accuracy estimates of all the results are approximately 10% to 15%. Approximate values of velocity and effective Q versus depth at the Houston testing site are offered below. Mean Slant-Path S-Wave Velocity Versus Depth (Houston, Texas) Figure 3.19 shows the mean slant-path S-wave velocities obtained from measurements in natural soils at three bore- hole locations (D2, E1, and F2) in the Houston area. The mea- surements at Borehole D2 were physically limited in depth by the inability of the Geoprobe system to push the downhole accelerometer probe deeper than 9 ft. Similar limitations were encountered in Borehole E1 but with a maximum push depth of 14 ft. In all cases, unlike at Manteno, the quality of the recorded signals was sufficient to allow the mean slant-path S-wave velocity analysis to be productive at all depth positions down to maximum push penetration. Figure 3.19. Mean slant-path S-wave velocity, Houston Boreholes D2, E1, and F2.

30 The higher velocities in Borehole F2 indicate a change in the soil seismic properties, characterized by an increase in shear modulus, at depths below about 10 ft. The uniform increase in slant-path S-wave velocity below 9 ft in Borehole E1 is indicative of a uniform soil exhibiting only the effects of soil column compaction. Interval S-Wave Velocity Profile (Houston, Texas) Figure 3.20a–c shows the derived interval S-wave velocity profiles in natural soil at Houston Boreholes D2, E1, and F2. The power-law regression curves fitted to these interval velocity profiles tend to smooth the numerical scatter in the derived values. Since interval velocity calculations were not productive at depths deeper than 10 ft in Borehole F2, the regression curve for that borehole does not show a higher velocity gradient at the deeper depths. Mean Slant-Path P-Wave Velocity Versus Depth (Houston, Texas) Figure 3.21 shows the mean slant-path P-wave velocities obtained from measurements in natural soil at three test bore- holes (D2, E1, and F2) very near the boreholes used for S-wave measurements in the Houston area. The maximum Geoprobe system push depths in these boreholes were 8 ft, 13 ft, and 16 ft, respectively. The slant-path P-wave velocities in Borehole D2, although limited to a depth of 8 ft, are uniform and consistent with the S-wave velocities in the same soil column. The slant-path P-wave velocities in Boreholes E1 and F2 show wide varia- tions at several depth points below about 4 ft that are incon- sistent with the S-wave velocities. These variations are attributed to possible differences in the separate borehole measuring conditions (ground coupling, etc.) or artifacts in the data analysis. Interval P-Wave Velocity Profile (Houston, Texas) Figure 3.22a–c shows the derived interval P-wave velocity profiles in natural soils at the three Houston Boreholes D2, E1, and F2. The power-law regression curves fitted to the derived interval velocities bring the scattered values into better agree- ment with the S-wave velocities and, lacking better infor- mation, at these P-wave test boreholes, provide reasonable characterizations of P-wave velocity profiles at these sites. Effective Slant-Path S-Wave Q Versus Depth (Houston, Texas) Figure 3.23a–c shows the effective slant-path S-wave Q and related attenuation rate in natural soils at three boreholes (D2, E1, and F2) in the Houston area. a.) Figure 3.20. Interval S-wave velocity, Houston boreholes: (a) D2, (b) E1, and (c) F2. (Continued on next page.)

31 c.) b.) Figure 3.20. Interval S-wave velocity, Houston boreholes: (a) D2, (b) E1, and (c) F2. (Continued from previous page.)

32 Figure 3.21. Mean slant-path P-wave velocity, Houston Boreholes D2, E1, and F2. a.) Figure 3.22. Interval P-wave velocity, Houston boreholes: (a) D2, (b) E1, and (c) F2. (Continued on next page.)

33 c.) b.) Figure 3.22. Interval P-wave velocity, Houston boreholes: (a) D2, (b) E1, and (c) F2. (Continued from previous page.)

34 Sufficient data obtained at depths to 9 ft in Borehole D2 and deeper in Boreholes E1 and F2 allowed reasonable determination of the effective Q values of soils at the three locations. The calculated values for Boreholes D2 and E1 are shown to illustrate the variations in the data. However, the data in Borehole F2 were widely scattered and the effec- tive slant-path S-wave Q values are shown represented by a smoothed curve. The resulting Q values are in the ranges 22–32 for Borehole D2, 20–26 for Borehole E1, and 9–22 for Borehole F2. Interval S-Wave Q Profile (Houston, TX) Figure 3.24a–c shows the derived interval S-wave Q and related attenuation rate profiles in natural soils at the three boreholes (D2, E1, and F2) in the Houston area. Variations in the interval S-wave Q values at Boreholes E1 and F2 were widely scattered and are shown only by their smoothed values. The interval Q values in Borehole D2 are acceptable within the top 5 ft of the ground and are repre- sented at deeper depths by the mean value of 28.6 in the top 5 ft. The interval Q values in Boreholes D2 and E1 are consis- tent with the effective slant-path S-wave Q values. However, the interval Q values in Borehole F2 exhibit a much wider range with increasing depth, largely because of the apparently higher S-wave velocity at depths below 10 ft and the data smoothing process. Effective Slant-Path P-Wave Q Versus Depth (Houston, Texas) Figure 3.25a–c shows the effective slant-path P-wave Q and related attenuation rate derived in natural soils at the three boreholes (D2, E1, and F2) in the Houston area. The effective slant-path P-wave Q values showed more scat- ter than the corresponding S-wave Q values. These variations were smoothed to avoid extreme values, resulting in slightly increasing Q values versus depth. The calculated P-wave Q val- ues are in the range 13–21 in Borehole D2, 15–27 in Borehole E1, and 31–47 in Borehole F2. Interval P-Wave Q Profile (Houston, Texas) Figure 3.26a–c shows the derived interval P-wave Q and related attenuation rate profiles in natural soils at the three boreholes (D2, E1, and F2) in the Houston area. The derived interval Q values are widely variable in the three boreholes. The interval P-wave Q values in Boreholes D2 and E1 are represented by their mean values of 13 and 18, respectively. The values in Borehole F2 increase with depth, represented by a fractional power-law curve over the range 31–85. a.) b.) c.) Figure 3.23. Effective slant-path S-wave Q, Houston boreholes: (a) D2, (b) E1, and (c) F2.

35 a.) b.) c.) Figure 3.24. Interval S-wave Q profile, Houston boreholes: (a) D2, (b) E1, and (c) F2. a.) b.) c.) Figure 3.25. Effective slant-path P-wave Q, Houston boreholes: (a) D2, (b) E1, and (c) F2.

36 Soil Seismic Properties Collective Results The soil seismic properties measurements were a qualified success in determining S-wave and P-wave velocities and attenuation rates at the three regional field sites. These results are considered to be representative of the seismic parameters of interest and are empirically estimated to have an accuracy tolerance of approximately 10% to 15%. The measurement objectives were demanding in that new custom-designed downhole instrumentation was required for depth-profile S-wave and P-wave measurements at seismic frequencies up to 1,600 Hz in soils typical of locations where underground utilities are installed. Programming, assembly, and checkout of the Dewetron digital data acquisition system required several revisions and adjustments in both software and hardware to implement the required toneburst excitation and data recording process. The surface seismic vibrator (OES MicroVib) was equipped with a custom-designed ground- coupling device for the purpose of achieving point-source S-wave and P-wave radiation. This small-contact-area ground- coupling technique was successful for S-wave operation but ineffectual for P-wave operation, resulting in a larger ground- coupling base plate being used for the P-wave measurements. Soil conditions at the Traverse City, Michigan, trial testing location selected for system assembly and checkout were unsuitable (because of their loose and unconsolidated granular nature) for use in quantitative tests at downhole depths greater than just a few feet below surface. This condition limited the full-scale testing and perfection of the planned Geoprobe sys- tem push-probe method and, therefore, limited the seismic measurement checkout tests to basic equipment operability. As a result, only minimal seismic measurements were possible, and they were not sufficient to provide preliminary downhole response measurements and sample data for evaluating the complete system performance. At the first regional field test site at Manteno, Illinois, the combination of vibrator source frequency response, ground- coupling effects related to field conditions, and possible anom- alous propagation responses associated with the intended high-frequency downhole sensor measurements was observed to cause significant interference effects in the recorded data. In particular, the borehole diameter and downhole probe dimen- sions are comparable with the seismic wavelengths at fre- quencies of about 1,200 Hz and higher. Therefore, operating the system at these higher frequencies may potentially excite guided waves along the borehole and spurious responses in the downhole sensor, resulting in accelerometer interference effects. These effects were noticeable in the recorded data, although not separable into their individual causes as related to the surface and downhole conditions. In general, the interfer- ence effects limited the data analysis to toneburst frequencies in the range of 100–1,100 Hz. a.) b.) c.) Figure 3.26. Interval P-wave Q profile, Houston boreholes: (a) D2, (b) E1, and (c) F2.

37 Irregular system responses at 1,100 Hz and lower caused amplitude and phase variations that introduced errors in the derived propagation velocity and attenuation. Unconventional analysis methods were devised and adapted to minimize these errors, with trade-offs between obtaining approximate values of S-wave and P-wave velocity and Q and their depth trends in the tested soils and not obtaining any useful information from the data. The measured S-wave and P-wave velocities and Qs for the three regional field sites are plotted in Figure 3.27, Figure 3.28, Figure 3.29, and Figure 3.30 for comparison. These compari- sons show that the S-waves and P-waves at the three regional field sites have the same general trends in soil seismic proper- ties versus depth. The mean combined velocity and Q profiles for S-waves and P-waves at the three Houston borehole sites are shown as bold dash-dot lines. The respective velocity and Q magnitudes at the three sites are also relatively close in comparison, suggesting that the seismic properties of near-surface natural soils are governed, in a first-order sense, by their shallow depth and natural com- paction in combination with moisture content. Their simi- larities, as illustrated above, suggest that they may be merged to provide useful composite profiles relevant to these particu- lar sites and perhaps other similar sites. These curves follow a fraction power-law depth dependence with exponents of 0.521 and 0.586 for S-waves and P-waves, respectively. Estimated Seismic Reflection System Performance A cross section of near-surface seismic propagation and reflection applicable to subsurface pipe detection is shown in Figure 3.31. The seismic source and sensor units generate and receive horizontally polarized shear (SH) waves and are located relatively close together to provide near-vertical transmission and reflection paths in a soil medium having a Figure 3.27. Mean slant-path velocity and Q profiles at three regional sites, S-waves. Figure 3.28. Interval velocity and Q profiles at three regional sites, S-waves.

38 positive velocity gradient versus depth. The features indicated in the cross section were divided into categories used to formu- late the two-way transmission-reflection loss based on (1) the noise threshold established by the proposed seismic sensor array, and (2) the practical detection limit established by the dynamic range of the seismic data recording system. Those categories, which were characterized completely by the subcontractor Owen Engineering Services, include • Seismic transmission and reflection paths dependent on the soil velocity and attenuation depth profiles; • Pipe target reflection; and • Seismic operating system, including ground coupling. Pipe Detection Performance (Houston, Texas) Figure 3.32 shows the detection depth of different size pipes using an SH-wave seismic reflection system at the Houston Borehole D2 test site. SH-wave pipe detection versus depth studies were completed for all borehole locations of the seis- mic soils properties testing sites where data quality was suffi- cient to do so. Those results show trends similar to the Houston boreholes represented here. The graphs shown here represent three distinct locations within the Houston metro area sepa- rated by several miles. The practical detection threshold indi- cates the probable maximum detection depth for pipes with diameters in the range of 3–6 in. Larger-diameter pipes are detectable at deeper depths. The slopes of the curves for each pipe size are governed primarily by the shear wave attenuation in the Houston D2 soil column and the target pipe diameter. The measurements at Borehole D2 were limited to the indi- cated depth of 9 ft by the Geoprobe system push refusal. The starting sweep frequency is high in this case because of the higher S-wave velocity; and the upper frequency limit is arbi- trarily selected as being adequate for resolving the detection depth of the pipes. Conservative estimates of the pipe detection Figure 3.29. Mean slant-path velocity and Q at three regional sites, P-waves. Figure 3.30. Interval velocity and Q profiles at three regional sites, P-waves.

39 Figure 3.31. Pipe reflection cross section in near-surface soil medium. Figure 3.32. SH-wave reflection system performance, Houston Borehole D2.

40 depths for the different size pipes are indicated by extending the slopes of the curves down to the practical threshold limit (dashed lines). Pipes having a diameter of 3 in. are detectable at a depth of about 10.5 ft, and larger pipes are detectable at deeper depths. Pipes 6 in. in diameter are detectable at depths of about 17 ft. Because of the dominant effects of the soil properties, as long as the system source excitation frequency sweep signal starts at a time when a low frequency limit equals the frequency range for which distinct pipe reflections occur [or (fL = fr)], the indicated pipe detection depths (based on the 120-dB practical system dynamic range threshold) will not change significantly when modest changes are introduced in the other operating system parameters. However, as a cautionary note, because the velocity and Q values were not measured at depths greater than 9 ft, the soil conditions at greater depths are actually unknown. Figure 3.33 shows the detection depth of different size pipes using an SH-wave seismic reflection system at the Houston Borehole E1 test site. The practical detection threshold indi- cates the probable maximum detection depth for pipes in the 3-in.-to-6-in.-diameter range. Larger-diameter pipes are detectable at deeper depths. The slopes of the curves for each pipe size are governed primarily by the shear wave attenuation in the Houston E1 soil column and the pipe diameter. The velocity and attenuation analyses for Borehole E1 were not practical at depths below 14 ft because of Geoprobe system push refusal. The shear wave velocity and Q at this borehole site were lower than at the Borehole D2 site. The selected sweep frequency range of 400–1,600 Hz is a nominal range compat- ible with detecting pipes in the 3-in.-to-24-in.-diameter range in Houston E1 soil. Conservative estimates of the pipe detec- tion depths for the different size pipes are indicated by extend- ing the slopes of the curves down to the practical threshold limit (dashed lines). Pipes having a diameter of 3 in. are detect- able at a depth of about 5.5 ft, and shallower. Pipes 6 in. in diameter are detectable at depths of about 11.5 ft. Because of the dominant effects of the soil properties, as long as the system source excitation frequency sweep signal starts at fL = fr, the indicated pipe detection depths (based on the 120-dB practical system dynamic range threshold) will not change significantly when modest changes are introduced in the other operating system parameters. Figure 3.34 shows the detection depth of different size pipes using an SH-wave seismic reflection system having the design and operating parameters (described above) at the Houston Borehole F2 test site. The practical detection threshold indi- cates the probable maximum detection depth for pipes in the 3-in.-to-6-in.-diameter range. Larger-diameter pipes are detectable at deeper depths. The slopes of the curves for each pipe size are governed primarily by the shear wave attenuation in the Houston F2 soil column and the target pipe diameter. The velocity and attenuation analyses for Borehole F2 were not practical at depths below 9 ft because of limited data quality. The shear wave velocity and Q at this borehole site were approximately the same as those at Borehole E1 but lower than at Borehole D2. The selected sweep frequency range of 400–1,600 Hz is a nominal range compatible with detecting pipes in the 3-in.-to-24-in.-diameter range in Houston F2 soil. Conservative estimates of the pipe detection Figure 3.33. SH-wave reflection system performance, Houston Borehole E1.

41 depths for the different size pipes are indicated by extending the slopes of the curves down to the practical threshold limit (dashed lines). Pipes having a diameter of 3 in. are detectable at a depth of about 5.5 ft and shallower. Pipes 6 in. in diam- eter are detectable at depths of about 11.5 ft. Pipes larger than 8-in. diameter are detectable at depths of about 15 ft and deeper. Because of the dominant effects of the soil properties, as long as the system source excitation frequency sweep signal starts at fL = fr, the indicated pipe detection depths (based on the 120-dB practical system dynamic range threshold) will not change significantly when modest changes are introduced in the other operating system parameters. Seismic Soils Properties Testing Assessment The project team performed the soil testing analysis on the data gathered by pushing an instrumented cone into the ground and measuring seismic parameters at 1-ft intervals on the way down. The testing was conducted at four regional sites: Traverse City, Michigan; Manteno, Illinois; Houston, Texas; and Manassas, Virginia. The Geoprobe with instrumented cone was not always in perfect contact with the soil surround- ing the tip during the tests, leading to gaps in the data coverage in many of the test holes. Owen Engineering Services and Psi-G struggled to find at least one completely satisfactory data set from the multiple tests that were done at each site. This diffi- culty was predictable given past experience in making in situ measurements; several test holes were made at each site to attempt to make up for the predicted ground-coupling prob- lems. The seismic team felt that enough good data were obtained at each site for the planned analysis to be done. The research that has been done so far has been largely aimed at determining the basic physical properties (seismic velocity and wave attenuation) of shear waves in soils when using frequencies in the range of 100–1,600 Hz. This work was performed because no literature exists in this area of study; without knowing these parameters, the project team could not establish the specifications for seismic measure- ment and imaging systems. Three basic results came out of this work. First, velocity and attenuation of shear waves in a wide range of soils are within ranges that make it possible to specify a measurement system that can be operated within the capabilities of modern electronic components, such as analog-to-digital converters and amplifiers. Second, the tests demonstrated that shear waves could be generated and prop- agated in the subsurface soils within the frequency range of interest, most of the time. In some cases the linearity of soil behavior is in question and further testing to track down this variable will be done in the final work phase as described below. Third, subsurface soil environments are even more complex and heterogeneous than the R01B team expected. The last point has the effect of requiring more work than expected in structuring and operating sources, constructing receiver arrays, and performing data processing. Figure 3.34. SH-wave reflection system performance, Houston Borehole F2.

42 Seismic Source Evaluation A major issue in evaluating the soil testing data was that many data sets contained data that didn’t meet quality criteria for pro- cessing and computing of the desired parameters. With further study, the project team determined that many of the data sets were experiencing nonlinear results because larger than neces- sary signals were being generated by the microvibrator source. The MicroVib was developed for study of deeper geologic tar- gets, perhaps greater than 500 ft. The source energy necessary to do that is so large as to affect the structure of the surface soils on which the source sits and through which signals must pass. While the MicroVib’s operating system was more than capable of reproducing the high-frequency signals required for this work, its output could not be “dialed down” to levels that would not drive the soils nonlinearly. That resulted in higher source levels than were appropriate for shallower depths and uncon- solidated soils. A number of issues surrounding the source operations needed to be altered to solve this problem. By the time the project team figured out the issues with this system, the project was almost complete. As a way to move as far forward as allowed by the time and money available, the team decided to design a completely new source generator that would account for the issues discovered in the soil testing. The prototype new shear wave seismic source is shown in Figure 3.34A. It is over-engineered to solve a range of problems and suspected problems and is clearly not ready for more than very rudimentary prototype testing. As the prototype goes through further development, it will be refined to get to a commercially workable system. This prototype and a receiver prototype have now been developed. They work well enough on the bench to say that the team is very close to proving feasibility, but the system has not yet been tested in the field. Seismic Modeling Software Dr. Nevin Simicevic and the team at Louisiana Tech Univer- sity (LTU) have produced the basic code for 3-D modeling of seismic signals. Models were created using soil parameters derived from the seismic field testing. This work progression has been significant and will become increasingly valuable in the interpretation of test results moving forward toward a seismic prototype system. Dr. Mark Baker of Geomedia Research & Development is already beginning to use the results for his seismic prototype development work on SHRP 2 Renewal Project R01C, Innovations to Improve the Extent of Locatable Zone. Wave Propagation Program Modeling WPP Modeling—Soil Properties Soil densities, the velocities of the waves in the soil, and the coefficients needed to describe the wave attenuation are the soil properties needed to describe the elastic and acoustic wave prop- agation using the finite difference time-domain (FDTD) method. Depending on the type of the soil used in the simula- tion, they should be measured or, if the data exist, taken from the literature. An example of the parameters used by LTU is taken from the work of Michael L. Oelze, William D. O’Brien, and Robert G. Darmody—Measurement of Attenuation and Speed of Sound in Soils (Oelze et al. 2002). Table 3.1 describes the type of soil used, and Table 3.2 describes the soil mean bulk densities. Table 3.3 describes the mean acoustic propagation speed in the soil, and Table 3.4 shows the mean attenuation coefficients valid in the frequency range between 2 kHz and 6 kHz. The attenuation coefficients in Table 3.4 are obtained from the relation 20log 0 x p x p ( ) α = − where p(x) is the measured pressure as a function of the thickness, x, traveled by the acoustic wave. The measured pressure is normalized to the pressure at the sound source. The WPP code uses quality factors, Q, instead of the attenu- ation coefficients. The relation between the Q factors and the attenuation coefficient is Q f v = • • pi α where v is the speed of the wave, and f is the frequency. Figure 3.34A. Newly designed seismic source preprototype.

43 Reprinted by permission, ASA, CSSA, SSSA. Table 3.1. Chemical Composition of a Soil Material (Oelze et al. 2002) Table 3.2. Mean Bulk Densities of a Soil Material (Oelze et al. 2002) Reprinted by permission, ASA, CSSA, SSSA. Reprinted by permission, ASA, CSSA, SSSA. Table 3.3. Mean Propagation Speed of Acoustic Waves in Soil (Oelze et al. 2002) Reprinted by permission, ASA, CSSA, SSSA. Table 3.4. Mean Attenuation Coefficients of Acoustic Waves in Soil (Oelze et al. 2002) While some of the soil properties described in this report apply only to the propagation of acoustic or P-waves, the same formalism applies to S-wave propagation, which will be studied in the rest of this report. The parameters for S-wave propagation in soils are less available in the literature than those for the propagation of P-waves. WPP Modeling—Sources In the WPP code, the soil can be deformed by a time-dependent force inducing the propagation of elastic or viscoelastic waves. Several shapes of time-dependent sources are possible; but, generally, they can be separated into three groups: a unipolar pulse with a power spectrum in the frequency domain peak- ing at 0 Hz (Figure 3.35), a multipolar pulse for which the power density at 0 Hz is zero (Figure 3.36), and the Gaussian window function used to represent a quasi-monochromatic wave (Figure 3.37). For the pulsed sources, the time durations of the source were chosen to conform to the frequency range requirement of the R01B project, 30–4,000 Hz. In addition,

44 because of the increase of wave attenuation with the increase in wave frequency, the selected frequency range was as low as possible, as long as it satisfied imposed requirements on the resolution. WPP Modeling—Attenuation To calculate the attenuation of the waves in the material instead of the attenuation coefficients, the WPP code uses the Q factors of the material. The relation between the Q factors and the attenuation coefficient is Q f v i i = pi α where v is the speed of wave, and f is the frequency. WPP uses the linear viscoelastic material model to simulate the wave attenuation. Using this model results in a requirement for more computer memory and more processing time. In a compromise between the computational requirements and physical accuracy, the number of relaxation mechanisms was set to three, resulting in solving 12 differential equations simul- taneously. In this case, the modeling was restricted to the fre- quency band of [w, 100w]; in that range, the quality factors (Q) for both P- and S-waves are assumed to be constant for a particular material. Quality factors used in the simulations were provided by the subcontractor Owen Engineering Services. An example Figure 3.35. Unipolar pulse example (in the shape of a very smooth bump and its power spectrum; the frequency parameter for this pulse was 1,000 Hz). Figure 3.36. Bipolar pulse example (in the shape of the time integral of the Ricker function, proportional to the time derivative of the Gaussian function; the frequency parameter for this pulse was 1,000 Hz). Figure 3.37. Gaussian window of monochromatic wave example (with central frequency of 200 Hz; the number of windowed cycles was 12).

45 of provided parameters used in the simulation is shown in Figures 3.24 through 3.26. WPP Modeling—Computational Requirements The factors contributing to the computational requirements (RAM and number of CPUs) depend on the size of the physical volume, the frequencies and velocities of waves, and the atten- uation parameters. For a realistic and useful simulation, an ideal computational volume would be of the size 400 cm × 400 cm × 300 cm. The computation volume of 400 cm × 400 cm × 300 cm, if required to be discretized by cells of the size 5 mm × 5 mm × 5 mm, results in a total of 3.84 × 108 computational cells, requiring large RAM and a large number of CPUs. Imposing attenuation may result in even smaller discretization cells, more memory, and more CPUs. Such a computation can only take place in a high-performance com- puting (HPC) environment. The WPP code is written to take advantage of the parallelization in a HPC environment; and Louisiana Tech University has access to the Louisiana Optical Network Initiative (LONI), providing a powerful system of high-performance computers. Testing of the speed of WPP was performed at the LONI system of supercomputers; the speed of computation was measured as a function of the number of CPUs involved. Using the same volume size and discretization, the model was run for a total of 9,468 time steps. Identical calculations were carried out on 8, 20, 40, 80, and 100 processors separately. Figure 3.38 shows the time in minutes it took the computa- tion to finish, as a function of the number of processors used. Calculation that takes 309 minutes using 8 processors takes only 24 minutes using 100 processors (~13 times faster). It can also be observed in Figure 3.38 that the increase of computation speed is not linearly proportional to the num- ber of processors. WPP Modeling—Monitoring of Quantitative Values Many physical situations were simulated using the WPP com- puter code, including the propagation of S- and P-waves through a soil with or without a reflecting target, with or with- out attenuation, and using pulsed or quasi-continuous waves. The output from the simulation consisted of physical quanti- ties in the form of a time series of the values of the displace- ment at specific points in space—underground or on the ground. Particularly, to estimate possible measured values, a network of imaginary motion sensors was positioned on the ground. They were used to study the difference in the strength and the time of arrival of the reflected signal. The schematic of the sensor positions is shown in Figure 3.39. The amplitude of the displacement and the time of arrival of the reflected pulse were recorded at those positions throughout the entire simulation run. The results from the simulation of pulse propagation and the reflection from buried pipes were then used to test possible signal processing, analysis methods, and target recognition. The goal was to determine the positions of buried pipes from the time series of the values of the displace- ment, or its time derivatives, recorded at the position of sensors. WPP Modeling—Simulation of Propagation of Shear Waves Without Attenuation A series of simulations of the signal propagation was performed assuming no physical attenuation of the waves. Quantitative values of the displacement were recorded as a time series at Figure 3.38. Total computation time to complete as function of the number of processors, in minutes. Figure 3.39. Schematic of the position of a source and motion sensors at the top of the soil.

46 chosen points under or on the ground. The unipolar case was represented by the shear pulse excitation in the form of the very smooth bump along the x-axis on the surface of the ground. The frequency parameter, w, was the same as in the case of the bipolar excitation—1,000 Hz. The shape of this pulse and its power spectrum are shown in Figure 3.35. The pulse was reflected from a plastic pipe with the diameter of 1 ft, positioned 115 cm under the surface of the soil. The soil and the plastic pipe were modeled in the same way as in the case of the bipolar excitation. The size of the computational volume was also the same. For the case of unipolar excitation, the change in the shape of the x-component of the S-wave dis- placement in the x–z plane as the pulse propagates through the soil and reflects from the plastic pipe is shown in Figure 3.40. Figure 3.40. The propagation and reflection of the x-component of the unipolar shear pulse displacement in the x–z plane. (The animation is at http://www.phys.latech.edu/~neven/uit/final.)

47 Figure 3.41. The propagation and reflection of the z-component of the unipolar shear pulse displacement in the x–z plane. (The animation is at http://www.phys.latech.edu/~neven/uit/final.) The y-component of the same S-wave displacement is small in the x–z plane. The propagation of the z-component is shown in Figure 3.41. Also in this case, the quantitative values of the displacement were shown at chosen points under or on the ground. The values of the x-component of the S-pulse displacement for each time step of the computation are shown in Figure 3.42 at positions of 25, 50, and 75 cm under the ground, between the source of the pulse and the plastic pipe. The same study was done by replacing the plastic pipe with water. Furthermore, to test the target recognition software, simulations were per- formed in which the pulse propagated through the soil and was reflected off multiple plastic pipes, each with diameters of 1 ft and positioned at different depths.

48 Figure 3.42. Values of the x-component of the unipolar S-wave displacement for each time step, at positions of 25 (top), 50 (middle), and 75 (bottom) cm underground (the form of excitation was a smooth bump; the transmitted and reflected pulses are shown). Figure 3.43. The experimental vibrator pulse and its Fourier power spectrum. While still not introducing the attenuation, LTU attempted to simulate the propagation using the pulse as measured in the trial seismic soil test in Traverse City, Michigan. LTU downloaded the vibrator data and performed the Fourier analysis. The shape of the pulse and its Fourier spectrum are shown in Figure 3.43. The Fourier analysis shows that the power spectrum peaks at a frequency of ~10 kHz, much higher than the frequency used in previous simulations. In the simulation, an experimentally obtained pulse shape was modified into the pulse shape compatible with one of the input modes of the Wave Propagation Program (WPP). The physical properties of such a pulse differ very little from the physical properties of an experimental pulse. The shape of the inputted pulse is shown in Figure 3.44. The simulation of the propagation of the Traverse City shear wave was performed in a volume equaling 200 cm × 200 cm × 170 cm, which was due to the power spectrum peaking at a fre- quency of ~10 kHz, discretized into 1,001 cubes × 1,001 cubes × 851 cubes for an overall total of 8.53 × 108 grid points. The length of a cell side was 2 mm, and the velocity of the shear wave was 214 m/s, resulting in ~10 samples per wave length at a frequency of 10 kHz. The stepping time was 5.41 × 10-7 s. The pulse propagated through the soil and was reflected by a plastic pipe that was the diameter of 1 ft, positioned at a depth of ~140 cm vertically under the position of the source. The P-wave velocity of the soil was 432 m/s, and the S-wave veloc- ity was 214 m/s. For the plastic pipe, the P-wave velocity was

49 Figure 3.44. The shape of the vibrator pulse as input to WPP. 2,458 m/s and the S-wave velocity was 1,164 m/s. The pro- gram ran on more than 120 processors over a time period of 30 hours, making 30,000 time steps. The execution of such a program is at the limit of LONI capabilities, assuming that the resources are shared with other users. At this stage, no physical pulse attenuation was computed. Assuming geometrical, but not physical attenuation, it was shown that such pulse can detect a buried pipe; but due to the frequency content and width of the pulse, the detection was not as straight forward as with simple pulses used in previous simulations. To facilitate the use of time-of-flight techniques to reconstruct targets positions, the Hilbert transformation was used. WPP Modeling—Simulation of Propagation of Shear Waves with Attenuation Simulations that include physical attenuation require more computational resources than simulations without attenuation. Still, several simulations, which included physical attenuation of the propagation of shear waves, were performed. Depending on the frequency of the wave and the allowed time of the exe- cution, the conditions of the simulations, including the size of the physical volume and the time of propagation, varied. In the first set of examples, signals from the vibrator’s monochromatic tonebursts, as shown in Figure 3.12, were used. Since the time duration of a toneburst was 60 ms, larger computational volumes were required. In the first simulation, a shear wave in the form of a toneburst of 200 Hz was used. The pulse shape and its frequency power spectrum are shown in Figure 3.37. Because of its low-frequency content, the computation could be performed relatively fast in a larger volume. The shear wave was excited along the x-axis on the surface of the ground. The propagation was simulated in the physical vol- ume equaling 400 cm × 400 cm × 366 cm. The parameters of the soil—including its density, P-wave and S-wave velocities, and the Q factors—varied with the depth of the soil. They are shown in Table 3.5. The simulated pulse propagation is shown in Figure 3.45. The attenuation of the pulse is shown in Figure 3.46. While the simulation of the propagation of a shear wave in the form of a toneburst of 200 Hz produced results when the attenuation was incorporated in the computation, several steps were necessary to fully understand the physical meaning of those results. The contribution to the attenuation also comes from the dependence of the amplitude on the distance from the source; for spherical waves, this dependence is 1/r, on numerical precision, on reflection from different type of soils, and so on. To fully understand the physical meaning of the simulation results and to properly incorporate the attenuation, a systematic study of the propagation of a shear wave in differ- ent conditions was done. In the process, the sizes of the ampli- tudes were determined using the Hilbert transformation. To test the geometrical 1/r amplitude attenuation, a shear wave isotropic point source of the 200-Hz toneburst was excited along the x-axis on the surface of the ground. The propagation was simulated in the physical volume Table 3.5. Input Parameters for WPP Modeling Simulation Block vp vs rho z1 z2 qp qs 193 88 1155 0.000 0.305 17.4 14.5 253 127 1236 0.305 0.610 17.5 15.2 296 157 1286 0.610 0.914 17.6 16.2 331 183 1322 0.914 1.219 17.6 17.2 362 206 1352 1.219 1.524 17.7 17.8 388 226 1375 1.524 1.829 17.7 18.5 412 245 1396 1.829 2.134 17.8 18.9 434 255 1415 2.134 2.438 17.9 17.6 454 280 1430 2.438 2.743 18 21.7 474 296 1446 2.743 3.048 18.1 23.9 492 311 1460 3.048 3.35 18.2 24.7 508 326 1472 3.35 3.66 18.2 28.6 Note: vp and vs are depth-dependent P-wave and S-wave velocities, in m/s; rho is depth-dependent density, in kg/m3, and z1 and z2 are the limits of the soil depth, in m; qp and qs are depth-dependent Q factors.

50 Figure 3.45. The propagation of the simulated pulse in the x–z plane. (The animation is at http:// www.phys.latech.edu/~neven/uit/final.)

51 Figure 3.46. Time series of the strength of the x-component of displacement at depths of, from top to bottom, 50, 100, 150, and 200 cm, showing the signal attenuation (the source was at a depth of 6 cm).

52 Figure 3.47. Values of the 200-Hz toneburst’s amplitudes at different depths (normalized to the amplitude at the source, in agreement with the predicted 1/r attenuation). Figure 3.48. Values of the 200-Hz toneburst’s amplitudes at different depths with attenuation (normalized to the amplitude at the source, in agreement with the predicted 1/r attenuation). equaling 400 cm × 400 cm × 366 cm. The parameters of the soil included the density of 1,155 kg/m3 and P-wave and S-wave velocities of 193 m/s and 88 m/s, respectively. The amplitudes of the wave were measured at depths incre- mented by 50 cm. The results of the values of the ampli- tudes at different depths, normalized to the amplitude at the source, and their agreement with the predicted 1/r attenuation, is shown in Figure 3.47. The perfect agreement between the values of the amplitudes in the simulated wave propagation and the theoretical prediction was also used as a test of the simulation software. In the second test, to test the effects of the soil properties, the same shear wave isotropic point source of the 200-Hz toneburst was excited along the x-axis on the surface of the ground. Again, the propagation was simulated in the physical volume equaling 400 cm × 400 cm × 366 cm, but the param- eters of the soil were taken from Table 3.5. The amplitudes of the wave were measured at the depths incremented by 50 cm. The results of the values of the amplitudes at different depths, normalized to the amplitude at the source, are shown in Fig- ure 3.48. Also shown in Figure 3.48, the attenuation is greater than the 1/r geometrical attenuation. As expected, the transi- tion of the wave from one soil density and propagation veloc- ity to another one increased the attenuation. In the third test, to test the effects of the nonisotropic source distribution, several isotropic point sources of the 200-Hz toneburst were distributed at the surface of the ground and excited along the x-axis. The array of sources’ distribution mimicked the soft-soil vibrator used by sub- contractor Owen Engineering Services. The propagation was simulated in the physical volume equaling 400 cm × 400 cm × 366 cm, and the parameters of the soil were taken from Table 3.5. The amplitudes of the wave were measured at the depths incremented by 50 cm. The results of the values of the ampli- tudes at different depths, normalized to the amplitude at the source, are shown in Figure 3.49. Again, as expected, the attenuation is slower than the 1/r geometrical attenuation since the shape of the wave changed from the spherical-type wave to the plane-wave-type wave. Finally, a full simulation was put together with the physical attenuation. Densities of the soil, velocities, and attenuation Q factors were taken from Table 3.5. An array of several iso- tropic point sources of the 200-Hz toneburst mimicking the soft-soil vibrator used by Owen Engineering Services was posi- tioned at the surface of the ground and excited along the x-axis. The propagation was simulated in the physical volume of the size 400 cm × 400 cm × 366 cm. The amplitudes of the wave were measured at depths incremented by 50 cm. The results of the values of the amplitudes at different depths, normalized to the amplitude at the source, are shown in Figure 3.50. They are compared, in the same figure, with the results without the physical attenuation. The overall attenuation increased, but

53 Figure 3.49. Values of the 200-Hz toneburst’s amplitudes at different depths, MicroVib source (normalized to the amplitude at the source compared with the geometrical 1/r attenuation; the source of the wave was the soft-soil Owen vibrator). Figure 3.50. Values of the 200-Hz toneburst’s amplitudes at different depths, MicroVib source and physical attenuation (normalized to the amplitude at the source compared with the geometrical 1/r attenuation; the source of the wave was the soft-soil Owen vibrator). since the Q values are large, the attenuation did not increase as much as expected. To test the attenuation dependence on the wave frequency, the propagation of 400-Hz and 800-Hz tonebursts was sim- ulated. As in the case of the propagation of the 200-Hz tone- burst, densities of the soil, velocities, and attenuation Q factors were taken from Table 3.5. An array of several isotro- pic point sources was positioned at the surface of the ground and excited along the x-axis. The propagation was simulated in the physical volume of the size 400 cm × 400 cm × 366 cm. The amplitudes of the wave were measured at depths incre- mented by 50 cm; and the results of the values of the ampli- tudes at different depths, normalized to the amplitude at the source, were reviewed for 400-Hz tonebursts and 800-Hz tonebursts. As expected, since the attenuation of the wave increases with frequency, a significant decrease in the wave amplitudes was observed when the wave frequency increased. The 800-Hz tone- burst propagation is shown in Figure 3.51. When compared with the propagation of the 400-Hz toneburst, the attenua- tion is unexpectedly lower. Since the wave is produced by a source array, this may be attributed to a narrower beam width for higher frequency. In the final test, the attenuation of the 1,200-Hz toneburst was compared with the attenuation of the 800-Hz toneburst; and, as expected, a significant decrease in the wave amplitudes was observed for a higher frequency. It was also observed that the beam width got narrower with higher frequency. The simulation of the signal attenuation for the pulsed source was done using the shear wave excited in the form of a very smooth bump along the x-axis on the surface of the ground. All the other parameters were the same as in the sim- ulation of the tonebursts. The pulse shape and power spec- trum of this pulse are shown in Figure 3.52. The amplitudes of the pulse in the form of a very smooth bump were mea- sured at depths incremented by 50 cm. The attenuation was comparable to the attenuation of the 200-Hz to 300-Hz tonebursts. Finally, LTU simulated the signal attenuation after reflec- tion from a target. The soil parameters and the geometry were kept the same as in the case of no target. The full 3-D simula- tion of the propagation of shear waves was performed in a volume equaling 400 cm × 400 cm × 367 cm. The shear wave was excited in the form of the very smooth bump shown in Figure 3.52. The source was the soft-soil vibrator. The pulse was reflected from plastic pipe that had a diameter of 2 ft and was positioned 250 cm under the surface of the soil. For the plastic pipe, the P-wave velocity was 2,458 m/s, and the S-wave velocity was 1,164 m/s, with no signal attenuation inside of the pipe. The propagation of the pulse is shown in Figure 3.53.

54 Figure 3.51. The propagation, in x–z plane, of the x-component of the 800-Hz toneburst. (The animation is at http://www.phys.latech.edu/~neven/uit/final.)

55 Figure 3.53 and the accompanying animation show that the reflected signal is clearly visible. They also show that the strength of the reflected signal is on the order of the magni- tude of physical and numerical background. The attenuation of the reflected signal was studied by per- forming two simulations: with the target and without the target. The amplitudes of the pulse were recorded at depths incremented by 50 cm. The physical and numerical back- ground was subtracted from the signal. The difference between the amplitudes with no target is shown in Figure 3.54. The reflected amplitude is clearly visible at the depth of 150 cm (at the time 23 ms); its value is comparable to the back- ground at the depth of 100 cm (at the time 26 ms) and becomes less than the background at the depths of 50 cm and 0 cm. While the reflected signal could be resolved with multiple sen- sors and sophisticated signal processing, the goal here was to estimate signal attenuation. To avoid complicating signal processing, the attenuation of the reflected signal was estimated in three steps: 1. The attenuation from the surface of the ground to the depth of 150 cm is calculated from the values of the ampli- tudes of the 500-Hz smooth bump at different depths when physical attenuation was included and found to be ~35 dB. 2. The attenuation due to the reflection is calculated by comparing amplitudes of the 500-Hz smooth bump at different depths when physical attenuation was included with the information from Figure 3.54 and found to be ~7 dB. 3. The attenuation from the depth of 150 cm back to the surface of the ground is again ~35 dB. The total attenuation of the signal for the reflection of a pulse shown in Figure 3.52 from a 2-ft-diameter plastic pipe Figure 3.52. Shape and power spectrum of a very smooth bump (the frequency parameter, , for this pulse was 500 Hz). Figure 3.53. The propagation and reflection of the x-component of the S-wave excited in form of a smooth bump with attenuation. (The animation is at http://www.phys.latech .edu/~neven/uit/final.) (Continued on next page.)

56 2 m under the ground is estimated to be ~77 dB. The simula- tion for a 1-ft-diameter plastic pipe 2 m under the ground showed that the total attenuation was ~87 dB. While the simulations of the propagation of the pulses in the form of a very smooth bump and the reflection from 2-ft- or 1-ft-diameter plastic pipes buried 2 m under the ground are only some of the possible cases, they were a testing foundation for all the tools needed in the simulation. The tools required to carry out the simulations include the software, a proper description of the source, a proper description of the properties of the soil and the target, and a proper description of attenuation. This case also showed that the simulation requires access to state-of-the-art supercomputers, but the computation can be performed in a reasonable amount of time. Signal Processing and Analysis Method The results from the simulation of pulse propagation were also used for testing possible signal processing and analysis and target recognition. The position of the buried pipes could be determined from the time series of the values of the dis- placement, or its time derivatives, at chosen points in space under or on the ground. In the case tested here, the imaginary motion sensors were positioned on the ground, as shown in Figure 3.53. The propagation and reflection of the x-component of the S-wave excited in form of a smooth bump with attenuation. (The animation is at http://www.phys.latech .edu/~neven/uit/final.) (Continued from previous page.)

57 Figure 3.39. The amplitude of the displacement and the time of arrival of the reflected pulse were recorded at those posi- tions through the entire simulation run. In the case of a pulsed source, the position of a pipe could be determined from the values of the displacements and their time of arrival. The possible methods used in electro- magnetics to estimate the position of a target, but applied here for the reflection of the elastic pulses, are shown in Figure 3.55. The method used to determine the positions of the buried pipes was based on the difference of the time propagation of the transmitted and reflected signals from different sensors. Because the distance between the target and a particular sen- sor is different, the round-trip time of flight of the signal is also different; as a consequence, the time of signal arrival can be used to locate the target. There are two possibilities: round-trip time of flight (RTOF) or time difference of arrival (TDOA). In the case of RTOF, the position of the target is obtained as a point of the intersection of sphere sizes that are determined from the round-trip time of flight and the speed of propaga- tion of the elastic pulse. In the case of TDOA, the position of the target is obtained by the intersection of hyperboloids; the size of each hyperboloid is determined from the difference of the time of the signal arrival to different receivers. The method of choice described in this report is the time difference of arrival. The advantage of using the TDOA method instead of RTOF is the less-strict requirements on the transceivers’ clock synchronization. The position of targets using the TDOA method was obtained by determining the time difference in the arrival of the reflected signal between different receivers. The points in space corresponding to the same TDOA between two receiv- ers define a hyperboloid, the size of which is determined by the difference in the time of signal arrival. Different pairs of a.) c.) b.) d.) Figure 3.54. Values of the background-subtracted amplitudes of the 500-Hz smooth bump reflected from target [target is 2-ft-diameter plastic pipe; depths are (a) 0 cm, (b) 50 cm, (c) 100 cm, and (d) 150 cm; MicroVib source]. Measurement Principles Propagation-Time Based Time of Arrival (TOA) Roundtrip Time of Flight (RTOF) Time Difference (TDOA) Angle of Arrival (AOA) Received Signal Strength (RSS) Figure 3.55. Localization methods used with ultra-wideband electromagnetic pulses (Vossiek et al. 2003).

58 receivers define different hyperboloids, and the point of their intersection defines the position of the target. Verification of TDOA Method Many simulations were performed to demonstrate the adequacy of the TDOA method. The results of some of the performed simulations, starting from a simple one-pipe case to a more complicated two-pipes case, are described in this chapter. The simulations were performed without the inclusion of the physical attenuation of the signal in the soil. The TDOA of the reflected signal between different receivers was found by first calculating the pair-wise cross-correlation of the signals. The TDOAs corresponded to the maxima in the cross-correlation. In a 2-D space constant, TDOA defines a hyperbola. The position of the target is obtained by the inter- section of hyperbolas corresponding to different TDOAs between different receivers. In the first simulation, the pulse propagated through the soil and was reflected from a plastic pipe that had a 1-ft diameter and was positioned at a depth of 125 cm. The P-wave velocity of the soil was 432 m/s, and the S-wave velocity was 214 m/s. For the plastic pipe, the P-wave velocity was 2,458 m/s, and the S-wave velocity was 1,164 m/s. The density of the soil was 1,650 kg/m3, and the pipe was 1,400 kg/m3. The computed x- and z-components of the elastic wave displacement propa- gation for this case are shown in Figures 3.56 and 3.57. The position of the buried pipe was obtained using the TDOA method. It is shown as color contours in Figure 3.58 and is superimposed to the exact pipe position shown as a black cir- cle. To determine the position, the data from five receivers were used. The x- and the z-components of the elastic wave displace- ment were tested; and, for this situation, the z-component showed a better position determination. The case of two buried pipes is more complicated. The pipes were positioned at depths of 2.0 m and 1.5 m—1.0 m left and 1.0 m right of the transmitter. The positions of buried pipes, obtained using the TDOA method, are shown as color contours in Figure 3.59 and are superimposed to the exact pipes positions shown as black circles. The simulation described in this report demonstrated that the TDOA method was able to resolve two buried pipes positioned at different depths under the ground. Modeling Summary WPP is written to take advantage of parallelization in a high- performance computing (HPC) environment, and the Loui- siana Optical Network Initiative (LONI) provides a powerful HPC environment for using WPP. WPP comes with a reper- toire of test scenarios that both demonstrate how the WPP may be used and provide checks as to whether WPP is running correctly. To date, the WPP scenarios that have run on LONI have been the default Lamb’s problem WPP test sce- nario and an adaptation of a scenario that modeled a pipe. In the case of the pipe scenario, a large amount of time was devoted to devising tools for the postprocessing of the output files from WPP scenarios. While WPP does supply a collec- tion of basic MATLAB scripts for visualizing its output, these scripts do not form a complete tool by themselves and are weakly documented. In the interests of performance and scal- ability to large data sets, a C++ program was written to per- form the postprocessing tasks. The postprocessing program addressed two major visual- izations of the data generated from a single WPP run. Using WPP’s output image files, the program can make movies and images of cross sections of the computational domain in dif- ferent scenarios. WPP also outputs sac files, which are files that contain the entire history of the simulation from the per- spective of a single point in the computational domain. Using these files, the postprocessing program can plot the histories for a given component of the computational domain (x, y, and z) at a particular location in the computational domain, or make movies of a series of such histories sampled along a path (presently a line) in the computational domain. The postprocessing program depended on two third-party libraries. MathGL was used to perform plotting and graphics, and some of the libraries in the Boost C++ libraries were used on several computational tasks that show up during the post- processing. The postprocessing program depended on a sin- gle third-party software, FFmpeg, which is used via system calls to produce movies. advanced tDeMI System Status and Findings The SHRP 2 R01B time-domain electromagnetic induction (TDEMI) component of research is also an innovative approach to digital geophysical mapping. The primary objec- tive was to develop a functional, advanced, sensor coil array prototype that continuously and dynamically records accu- rately positioned TDEMI geophysical data based on user set parameters and multiple coil configuration geometries. As mentioned before, the prototype TDEMI system was mod- eled after a first prototype version developed with other fund- ing by SAIC and the U.S. Naval Research Laboratory. This 5-by-5 coil sensor array has proven to be very successful for unexploded ordnance (UXO) classification under static data acquisition conditions. These successful contributions to an alternate application increased demand for the technology in the munitions response community, thus limiting the resources available for developing the R01B system during the project life cycle. The plan was to leverage the technology in two ways: (1) by collecting test data using the existing

59 Figure 3.56. The propagation through the soil and reflection from the plastic pipe of the x-component of the shear pulse. (The animation is at http://www.phys.latech.edu/ ~neven/uit/final.)

60 Figure 3.57. The propagation through the soil and reflection from the plastic pipe of the z-component of the shear pulse. (The animation is at http://www.phys.latech .edu/~neven/uit/final.)

61 5-by-5 array in the first prototype, and (2) by asking the developer to build a custom 5-by-1 array for integration for the SHRP 2 R01B EMI prototype. The TDEMI element of the R01B research work was performed by SAIC, G&G Sciences, and UIT, with support from Geosoft. After some difficulties in construction and testing by the subcontractors, UIT received the system in November 2011, and initial testing has been performed. The instrumentation received is shown in Figure 3.60, including Figure 3.58. Position of the buried pipe obtained using the TDOA method. Figure 3.59. Position of two buried pipes obtained using the TDOA method. a. Five coils plus a spare; b. Five preamps plus a spare; c. Approximately 20 ft of cabling for five transmit-receive coils; d. National Instruments electronics box with embedded computer and analog-to-digital converter (A/D) card; e. G&G electronics box attached to National Instruments box; f. Power cables for National Instruments box and G&G box; and g. Package of National Instruments software and CDs. Figure 3.60. TDEMI system instrumentation.

62 The system arrangement is shown in Figure 3.61. Equip- ment needed to run the system included the following: a. AC power or 18-30V DC for the National Instruments box; b. +12 and -12 volts for G&G box; c. Monitor, mouse, and keyboard (These were needed initially to set up networking on the embedded computer. Once the networking is configured, a laptop computer and a network cable can be used to “remote desktop” onto embedded com- puter.); and d. Data copied from the embedded computer using a USB thumb drive. EM3D is the software program developed by G&G Sciences for data acquisition. It has the following characteristics: a. Multiple windows for configuration, mapping, and GPS features; b. EM3D window used to configure i. Base file name; ii. Transmit pulse length, pulse repetition, and pulse stacking; iii. Single-shot or continuous operation; and iv. Single transmit coil or all five transmit coils sequentially; c. EM3DPlot window used to plot collected data and to export binary *.tem files to comma-delimited *.csv files. Screenshots of running program are represented in Figure 3.62. TDEMI Prototype Bench Testing Before delivery of the system, G&G Sciences performed static bench testing to ensure system stability and expected func- tionality. Testing involved the analysis of decay curves, so cali- bration results are offered here. Plots of test data collected are shown in Figures 3.63 through 3.72. a. Data were collected with only the center transmitter firing and with all five transmitters firing sequentially. b. Data collection was performed as follows: i. Single-shot with pulse on/off duration of 25 ms, nine repetitions of bipolar pulses, and the data collected and averaged three times. Total collection time for a single transmit coil is 2.7 s (0.025 × 2 × 2 × 9 × 3). For all five transmit coils, the collection time is 13.5 s. Figure 3.61. Block diagram of the advanced TDEMI system.

63 Figure 3.62. Screenshots of EM3D data acquisition software. Figure 3.63. Single-shot, measured transmit currents. ii. Continuous with pulse on/off duration of 8.33 ms, a single bipolar pulse, and no averaging. The time for a single transmit should be 33.33 ms (8.33 × 2 × 2). For five transmits sequentially, the collection time is 0.16667 s. c. Figure 3.63 plots the single-shot transmit currents measured as coils 0 through 4 (left to right) sequentially fire. The cur- rents vary slightly with coil and peak between 9 and 10 amps. d. Figure 3.64 plots the sequentially firing transmit currents as the array fires continuously. The transmit currents change slightly from shot to shot. e. To check that the five receivers were responding equally to the same signal, a 2-in.-diameter steel sphere was placed at a fixed distance above the center of each coil. Figure 3.65 plots the current normalized and background-subtracted

64 receive signals for this sphere’s data. The top plot shows the monostatic responses (same transmit-receive coil directly under sphere). The bottom plot shows the adja- cent bistatic response (transmit under the sphere, receive coil to the left and/or right). The amplitudes are equal, and the array was well calibrated. Figure 3.66 plots the same sphere data in separate plots. The solid lines are positive signal, and the dotted lines are negative. The consistent polarities indicate that all of the transmit-receive coils are wound in the same direction. f. Figure 3.67 plots the background signals measured on all 25 Tx-Rx pairs (array in the air and no metallic objects present). Note that there is significant ringing in the receive coil response to the transmit current turn-off in time gates Figure 3.64. Continuous shots, measured transmit currents. Figure 3.65. Current normalized and background-subtracted receive signals for 2-in. calibration sphere. (Continued on next page.)

65 Figure 3.65. Current normalized and background-subtracted receive signals for 2-in. calibration sphere. (Continued from previous page.) Figure 3.66. Two-inch calibration sphere above each transmit coil, mono- and bistatic terms.

66 earlier than 0.100 ms (from millivolts to volts in ampli- tude). The monostatic pairs have measurable ringing even after this (tenths of a millivolt). The sensor noise levels are on the order of 0.01 millivolts in the later gates. g. Figure 3.68 plots the result of subtracting one background from one measured several minutes later. For early gates, the ringing varied slightly, and there was significant spu- rious signal present. After 0.100 ms, the background dif- ferencing process removed most of the sensor ringing and results in sensor noise levels on the order of 0.01 millivolts. h. Plots were generated for a collection of continuous back- ground signals over periods of several minutes. The ringing response was the same as the single-shot data; the later noise levels were somewhat higher due to the lack of averaging. i. Figure 3.69 plots the result of subtracting the average background from continuous backgrounds. The noise in the later gates was closer to 0.1 millivolts. j. Figure 3.70 plots the recorded time between data samples for continuously collected data. The red symbols are for a single transmit coil firing constantly, and the black sym- bols are for all five transmit coils firing continuously in a sequence. The dotted lines are the expected data rates based on the pulse settings. It appears that there is some delay between each data collection. k. Figures 3.71 and 3.72 show the array signals from a 1-in.- diameter black steel pipe centered under coil 2 but orthog- onal to the array. The pipe was 30 cm and 60 cm below the array. The data were collected in the single-shot mode with repeated pulses and averaging. Figure 3.67. Background signals measured on all 25 Tx-Rx pairs.

67 These tests indicated that the system was performing com- parably to the original NRL TEMTADS 5-by-5 EMI array. SAIC sent a representative to assist UIT in the setup of the sensor platform and to conduct basic training on the R01B TDEMI system prototype. The sensor array was prepared for data acquisition transport on UIT’s existing metal-free trailer, as shown in Figure 3.73. The arrangement of sensors on the trailer offers a detection swath footprint of approximately 190 cm; the coils sit approximately 26 cm from the ground surface; and the GNSS antenna is centered on the sensor array and offset backward from the direction of acquisition travel approximately 35 cm. After assembling the power source and electronic connections, UIT conducted a series of static tests to become more familiar with the data acquisition software, the output data, and the general hardware configuration. TDEMI Prototype Data Acquisition Parameters The parameters selected for initial testing of the TDEMI pro- totype with the EM3D program are presented in Table 3.6 and offer an explanation of user-defined characteristics of the acquisition software. The list of parameters is copied from the *.AcqParams file. Delt is the data acquisition sampling rate in seconds. BlockT is the total measurement time in seconds. nRepeats is how many times the transmit cycle is repeated within BlockT. In this case, if one repeats nine times in 0.9 s, the transmit cycle is 0.1 s. Each transmit cycle consists of a bipolar pulse. With the Duty Cycle set at 50% (DtyCyc = 0.5), the transmit cycle consists of ¼th positive pulse, ¼th off, ¼th negative pulse, and ¼th off. So, in this case, the transmit pulse on time is 0.1/4 = 0.025 s, and the measurements made Figure 3.68. Signal shot measured on all 25 Tx-Rx pairs with differenced background.

68 while the pulse is off can only go out to 0.025 s. As one changes BlockT and nRepeats, one gets different pulse on/off times, which limits how late the time gates can go. A setting of BlockT = 0.1 and nRepeats = 1 would result in the same 0.025 on/off time and the same number of possible gates. A setting of BlockT = 0.03333 and nRepeats = 1 would only have a 0.00833-s pulse and gates out to 0.00833 s. The parameter nStk just repeats the measurement that many times and “stacks” the result to provide averaging and reduce general noise. If BlockT is 0.9 s and nStk = 3, the total time to take a measurement is 2.7 s. The AcqMode determines how the measurements made every 2.00E-06 over the total time block are averaged together. If AcqMode = FullWaveform, the entire set of measurements over BlockT is dumped. This setting is good for diagnosing problems, but it generates very large data files. In some cases, it may lock up the program. If set to DecaysOnly, the off time period after the negative pulse is subtracted from the off time period after the positive pulse (this cancels 60-Hz noise sources but adds signal). If there are multiple repeats, these are averaged together; but the final averaged data are still dumped at the 2.00E-6 sampling. The files are still quite large. Again, this is good for diagnosing noise sources and equip- ment problems. The third AcqMode is the one mostly used: DecaysDeci- mated. After dealing with the plus/minus pulse data and aver- aging over the nRepeats, the data are averaged over windows. The width of these windows is determined as a percentage by the parameter GateWid. A setting of 0.15 means the width of the window will be 15% of the time gate. This makes earlier Figure 3.69. Continuous data collection with background averages subtracted.

69 Figure 3.70. Measured data rate for single transmitter fired (red) and all five fired (black) (dotted line indicates expected rate). gates averaged over fewer of the 2.00E-06 samples and later gates averaged over more. The higher the percentage is set, the fewer the time gates will result. Averaging over wider gates reduces high-frequency noise. SAIC has worked with settings from 0.05 (5%) to 0.20 (20%) and is unsure what the limits are. The last setting is GateHOff, which is the “hold-off time” in seconds after the transmit pulse is off before data are kept. UIT was uncertain on the details of how this feature was supposed to work and so did not try changing it. SAIC suggested keeping it at the manufacturer’s setting—around 3.0E-5 to 5.0E-5. SAIC has used several of the G&G TDEMI sensors for the purpose of identifying buried unexploded ordnance and for stationary measurements; in those cases it has used BlockT = 0.9, nRepeats = 9, nStk = 3, and GateWid = 0.05. These set- tings allow multiple transmits (one system has 25) to be fired and data collected on the minute(s) time scale and seems to beat down noise sufficiently to detect items of interest at depths of interest. For smaller targets, SAIC has tried increas- ing the stacking to 10–20. The noise goes down with the square root of the stacks. However, at greater numbers of stacks, SAIC was concerned that instrument drift might be a factor. SAIC studied, and is still studying, the noise issues in col- lecting continuous, moving data with the TDEMI prototype. The sensor should not be moving very far on the time scale of the measurement. This results in increased noise because one cannot average over repeats and stacks. It was determined that the fastest measurement the TDEMI prototype can make is 1/30th of a second. At a speed of 1 m/s, the sensor moves 3 cm during the measurement. With BlockT = 0.0333, nRepeats = 1, and nStk = 1, the pulse time becomes 0.00833; so, later time gates are not possible. Greater gate widths reduce some of the noise, so SAIC suggested setting GateWid to 0.20 or higher. With BlockT = 0.0333, nRepeat could also be set higher—at the cost of shorter pulses and losing late time gates. TDEMI Prototype Dynamic Testing In preparation for in-service testing, UIT conducted a series of dry-run tests while operating the TDEMI prototype system in continual mode and integrated with GNSS. The objectives were to obtain further familiarity with the system setup, directly compare the TDEMI prototype data and detection capabilities with UIT’s current Geonics EM61 MKII towed sensor array

70 Figure 3.71. TDEMI measurements with 1-in.-diameter pipe 30 cm below sensor array. geophysical system, and gather information on optimal acquisi- tion settings for detecting underground metallic utilities. In addition to static testing, dynamic detection and positioning testing of the prototype TDEMI system was performed by lay- ing out metal chains of various sizes straight across the ground surface at various spacing and orientations. Operators of the system traversed the chains at differing speeds, using different acquisition parameter settings. Figure 3.74 is a general depic- tion of the type of field tests conducted for the dynamic testing. This initial attempt to view and evaluate the TDEMI data was performed with Geosoft’s Oasis Montaj software. By March 2012, the project team was able to put together a brief outline summary of the data collection methods and preliminary results and observations: 1. Stationary data collection a. Sensor array stationary test of metallic objects placed beneath. b. Data acquisition set to collect data continuously (including a single-shot mode), with i. Acquisition parameters: nRepeats = 1, BlockT = 0.0333 s. ii. This setting results in bipolar pulse (+on, off, -on, off), with a pulse on time of 8.33 ms and a total bipolar pulse length of 4 × 8.33 = 33.33 ms. iii. Transmitters were fired sequentially. Total time to fire all five was 5 × BlockT value plus some extra time gap between firings. Total time recorded in data files was roughly 0.21 s (5 × 0.0333 + 0.043 extra). This total time and array speed determine data sampling density along track when the array is moving. iv. Time gate width was set to 15%, resulting in 36 loga- rithmically spaced time gates out to 8.33 ms. Larger gate widths (20% to 40%) would average over more of the time decay and reduce high-frequency noise at the cost of fewer time gates over the same decay

71 Figure 3.72. TDEMI measurements with 1-in.-diameter pipe 60 cm below sensor array. Figure 3.73. Configuration and testing of the 5-by-1 TDEMI sensor array on UIT’s metal-free trailer. Table 3.6. Parameters for the Initial Testing of the TDEMI Prototype with the EM3D Program Parameter Setting Delt (s) 2.00E-06 BlockT (s) 0.9 nRepeats (n) 9 DtyCyc 0.5 (50%) nStk (n) 3 AcqMode DecaysDecimated GateWid 0.15 (15%) GateHOff (s) 3.00E-05

72 Figure 3.74. TDEMI prototype dynamic field test design.

73 range. At a gate width of 30%, there would be 16 time gates. c. Three test objects: large metal object (LMO), medium metal object (MMO), and small metal object (SMO). They are all 4-in.-long steel pipe nipples: 1.25-in., 1.0-in., and 0.75-in. diameters, respectively. d. Each object was placed under an array coil, and a data file was collected for 60–70 s with roughly 300 samples collected. e. Data files were collected for all three objects under each of the five coils in two different orientations. The orientations had the pipe lying flat with the long pipe axis pointing along the direction of travel (called par- allel) and the long pipe axis pointing orthogonal to the direction of travel (called perpendicular). f. One background file was collected for each object. 2. Stationary data processing and observations a. All data were normalized by the measured transmit current. b. For each file, the mean and standard deviation were calculated for each data channel (transmit-receive coil combination). c. The mean was also calculated for the background files by channel. d. The signal from a given object is found by subtracting the average background from the average measured signal. e. Figure 3.75 plots the average/background-subtracted array time-decay response as a function of time gate to the LMO in the two orientations directly under the center coil. i. The top five plots are for the LMO perpendicular, and the bottom five parallel. ii. The five plots across are for the five receive coils across the array. iii. Each plot has five colored curves. The color indi- cates which transmit is firing: 1 = black, 2 = red, 3 = green, 4 = blue, 5 = cyan. For five receives and five transmits, there are 25 data channels total. iv. The largest response is for the center receive coil when the center transmit coil is fired (green curve in center plot), since the object is directly below this transmit-receive pair. Coaxial transmit- receive pairs are termed monostatic. v. Solid curves are positive responses. Dotted curves are negative. Some of the separated transmit- receive pairs have negative responses. The sepa- rated pairs are termed bistatic. vi. The vertical black dotted line indicates roughly were the transmit-receive coil response has stabi- lized from the “ringing” of the transmit current turn off. vii. The slanted dotted line indicates a rough noise floor of the stationary array for the given data acquisition settings (see Figure 3.75). f. Figure 3.76 plots the standard deviation of the mea- sured time-decay responses for two of the stationary data files. This is the stationary root mean square (RMS) noise level of the array for the given continuous data and data acquisition parameters. (Note that Fig- ure 3.75 data are averaged over 300 samples and have a lower noise floor.) The five plots across indicate the receive coil, and the five colors indicate the transmit coil—as in Figure 3.75. i. The noise levels vary by receive coil. ii. The noise levels vary by data file. iii. It is probable that a time varying noise source is close to the array. iv. For later time gates, noise should fall off as t-½ as indicated by slanted, dotted lines. g. Figure 3.77 summarizes the response of the three objects (by row) and the two orientations (by column) as the objects are placed under each coil. i. The upper five curves are the monostatic responses with the object directly underneath. The color indicates which coil the object is under. ii. The lower curves are bistatic responses from adja- cent coils. iii. For the parallel cases, the monostatic and bistatic responses have the same curve shape. This is because the field from the transmit coils is only directed orthogonally to the long axis of the pipe nipple. iv. For the perpendicular cases, the bistatic response curves have a different shape. This is because some of the transmit field is directed along the long axis of the pipe nipple. v. The array transmit fields never intersect the test objects in three orthogonal directions. Because the objects are axisymmetric, only one or two unique directions are covered by the transmit field. h. Because of the limited “illumination” by the array transmit fields, a complete inversion of all dipole model parameters was not possible with the test object data. i. By constraining the object location to be directly centered under a given coil (in this case the center coil), a limited inversion of the polarization terms is possible. ii. Figure 3.78 plots the decay polarizations terms as a function of time gate for the three objects (LMO = black, MMO = red, and SMO = green). iii. For the parallel cases, the long axis is not illumi- nated; and only the equal, transverse polarizations are inverted. For the perpendicular cases, the long axis and transverse polarizations are both found.

74 Figure 3.75. Plot showing the average/background-subtracted array time-decay response as a function of time gate(s) to the LMO. Figure 3.76. Plot showing the rough noise floor of the stationary array.

75 Figure 3.77. Plot summing the array response of the three objects (by row) and the two orientations (by column) as the objects are placed under each coil. iv. For this relatively high signal-to-noise ratio (SNR), stationary data, the inverted polarizations were sufficient to distinguish the three pipe sizes. 3. Dynamic data processing and observations a. The baseline data were collected with the 8.33-ms pulse, 33.33 ms for the total bipolar pulse, and 0.21 s for all five transmitters to fire sequentially. b. At an array speed of 1 m/s, this means each set of trans- mit data is sampled roughly every 0.21 m along track. c. There are recorded time stamps only for each set of all five transmitters firing. It is assumed that the data for each transmit is uniformly spaced over this interval. d. The data acquisition program records and averages all of the GNSS data that come in while the five transmit- ters are firing. Because the GNSS data rate during field testing was slower than this time, many of the “5 trans- mit” data records are not updated in position at all (i.e., GNSS every 0.5 s and data record every 0.21 s).

76 e. To map the data, only new GNSS updates were kept. The GNSS antenna was centered behind the array. A course-over-ground had to be calculated from the GNSS trajectory and the GNSS position projected for- ward 0.35 m to the array. The recorded GNSS times with an arbitrary time shift were used to interpolate GNSS positions at GNSS times to transmit firing times. This shift was determined by looking at dynamic data collected in back-and-forth paths over two chains sep- arated at different spatial intervals and positions at a variety of orientations. f. Figure 3.79 plots the RMS noise measured with the array moving over a reasonable stretch of ground just before the placed chain. The dynamic noise levels are higher than the stationary levels and are worse on the array edges. The source of this noise has yet to be determined. g. Figure 3.80 plots rasters of the monostatic channels as the array drives back and forth over parallel chains laid on the ground surface. i. The outer sensors 1, 4, and 5 show the increased noise. Figure 3.78. Plot showing the decay polarizations terms as function of time gate(s) for the three metal objects. Figure 3.79. Plot displaying the RMS noise measured with the array moving over reasonable stretch of ground before chain emplacement.

77 ii. From top to bottom, the chains are spaced at 0-ft, 1-ft, and 2-ft apart. The array can resolve two chains at the 2-ft chain separation distance. iii. Figure 3.81 shows similar data for the 3-ft and 4-ft separations. h. Different array speeds and different acquisition param- eters were also tested. i. In addition to the basic fast rate (nRepeats = 1, BlockT = 0.0333), a medium rate of (nRepeats = 3, BlockT = 0.1) and a slow rate of (nRepeats = 9 and BlockT = 0.3) were used. These settings provide more averaging to reduce noise. ii. The medium rate records a new transmit sample every 0.56 s. At a speed of 1 m/s, this means a sam- ple every half meter or so. At the slow rate, there is a sample every 1.6 s with a 1 m/s distance of 1.6 m between samples. The slow data rate is not viable at reasonable acquisition speeds. Figure 3.80. Plot displaying rasters of the monostatic channels as the array drives back and forth over parallel chains laid on the ground surface at separation distances of 0 ft, 1 ft, and 2 ft.

78 iii. Speeds were varied from less than 0.5 m/s up to roughly 1.5 m/s. Figure 3.82 plots the responses over chains placed in an X pattern. At the fast data rate, the pattern is well resolved at a slow speed (top plot), but not as well at the fast speed (bottom plot). iv. Figure 3.83 plots the same X-pattern data at the medium data rate and the two speeds. At a slow speed, the pattern is still resolved. At the fast speed, the pattern dissipates. v. At the slowest data sampling rate, all of the chain patterns were completely washed out. vi. At the price of later time gates, two other data acquisition settings were considered to reduce noise: (1) One is nRepeats = 3 and BlockT = 0.0333. This setting produces an “on” pulse of 2.77 ms. The time for all five transmits to fire would still be 0.21 s, and at 1 m/s a sample every 0.21 m. (2) The other is nRepeats = 9 and BlockT = 0.1. Again, this is a 2.77-ms pulse; but with a time of 0.56 s per sample and a slow speed of 0.5 m/s, Figure 3.81. Plot displaying rasters of the monostatic channels as the array drives back and forth over parallel chains laid on the ground surface at separation dis- tances of 3 ft and 4 ft. there would still be a measurement roughly every 0.25 m. (3) In addition, the gate width parameter could be increased to 20–40% at the cost of fewer gates spaced over the same time-decay period. TDEMI Prototype Data Manipulation Although preliminary field testing was performed with the R01B TDEMI system, additional testing was needed before the in-service testing process to determine optimal acquisi- tion parameter settings; to compare against currently used TDEMI systems; and to determine effective ways to orga- nize, view, and process the expansive data sets generated by the prototype. UIT elected to conduct data evaluations with Geosoft Oasis Montaj software. Initial data management included two methods. The first involved viewing the data on a point-sample-by-point-sample basis. Each line of the Geosoft database represents a specific data point, and the array channels of each line represent the various Tx-Rx pairs’ decay curves, as well as the five signal transmit strengths.

79 The second involved viewing the data on a sensor-by-sensor basis. Each line of the Geosoft database represents a specific sensor, and the channels of each sensor line represent the exact numerical values for each point of the transmit-decay arrays. Using either method, the process proved to be quite time consuming, so UIT continued searching for optimal data acquisition parameter settings and data workflows to account for the relatively large volume of data being generated. UIT’s collaboration with Geosoft programmers resulted in the development of an import template with which TDEMI prototype data could be imported into a Geosoft database; each line of the database represents a specific sensor of the array. Individual channels of the database can carry values— such as Easting, Northing, Elevation, Heading, Pitch, Roll, Point Number, GPS Coordinated Universal Time (UTC), GateTime, and the various decay curve arrays of the Tx-Rx pairs. UIT data analysts beta-tested these new user-interface Figure 3.82. Plot showing the array responses over chains placed in an X pattern at speeds of ~0.5 m/s and ~1.5 m/s (fast data sampling rate).

80 functions and offered feedback to developers on a continual basis. Figure 3.84 illustrates the graphical user interface from this tool development of data import to a geophysical analysis software package. The primary notion was to evaluate the system’s performance as a digital geophysical mapping tool with which data are collected in continued and dynamic mode rather than in a static “cued interrogation” mode as is used for munitions response applications with similar TDEMI technologies. For the R01B TDEMI prototype, the raw TEM data files are recorded and stored on the field acquisition computer. The operator is allowed to set a file name prefix within the setup display window, and the EM3D acquisition software auto- matically names individual data sets sequentially and numer- ically thereafter. These raw files are generated in binary format and come with an extension of “.TEM.” To change raw files into a usable format, the operator is required to convert the .TEM data files into an ASCII format (.CSV). This procedure Figure 3.83. Plot showing the array responses over chains placed in an X pattern at speeds of ~0.5 m/s and ~1.5 m/s (medium data sampling rate).

81 is performed within the EM3D environment and is typically conducted at the end of a production day. This action results in the generation of two distinct files for each data set col- lected: the actual numerical data file (shown in Figure 3.85), and the user-defined acquisition parameters file (shown in Figure 3.86). On import of TEM data CSV files, a separate Geosoft data- base is created for each data session. The geophysical data are split and organized for each data set imported, such that each line of the Geosoft database corresponds to one sensor of the sensor array. A monostatic (Rx coil = Tx coil) data channel is also created within the database so that data can be viewed in a similar fashion to existing TDEMI digital geophysical map- ping systems in which transmit coils are the same (and equally positioned) as receiving coils. Given the vast amount of data acquired (compared with current digital geophysical map- ping EMI systems), sensor values are stored as arrays within each cell of the Geosoft database as shown in Figure 3.84, rather than individual numerical digits. Once the TDEMI prototype data are imported to Geosoft Oasis Montaj, data analysts are able to perform all data pro- cessing, data interpretation, and mapping procedures. Data points may be translated into the site-specific coordinate sys- tem and plotted against existing geo-referenced records to check for accuracy and completeness; data can be corrected for instrument latency and drift; data can be reviewed in pro- file view for multiple data channels (see Figure 3.87), and/or single channels can be gridded and color contoured for col- lective views within map space (see Figure 3.88). The software involves dynamically linked databases, maps, and 3-D views to give the analyst all the spatial and detection factors needed to perform and devise a well-informed subsurface assess- ment. Figure 3.89 shows a 3-D view in Geosoft Oasis Montaj; each time gate measurement is combined and gridded as a full block of data (voxel). The time gate intervals of the TDEMI prototype act like a pseudo-depth component; the earliest time gate measurements are depicted at the actual ground surface elevations, and later time gate measurements are depicted at deeper intervals. This is done by flipping the gate time about the 0-time axis and essentially making the gate times correspond to equivalent depth values below the ground surface. Although this does not provide any sort of true depth-to-target estimations, it does allow the analyst to view, manipulate, and rotate all of the EMI data within one graphical window. Finally, TDEMI data interpretations of suspect utility locations can be synthesized down, graphically drafted, and exported from the software. Mapping deliver- ables to engineering firms can be produced in a variety of compatible computer-aided design (CAD) and GIS formats. TDEMI Prototype Testing Summary During the course of SHRP 2 R01B research work, a practical prototype was constructed and tested using advanced TDEMI Figure 3.84. Screenshot of Oasis Montaj, representing customized data import function and resultant Geosoft database array channel for each Tx-Rx coil pair.

82 Figure 3.85. Sample raw excerpt (one data point) of TDEMI prototype data file. Figure 3.86. Sample raw TDEMI prototype acquisition parameters file.

83 Figure 3.87. Screenshot of Oasis Montaj profile view of monostatic channel for Sensor 1 (colored profile) and individual array view (smaller window) of one individual data point. Figure 3.88. Screenshots of various TDEMI prototype grids in Oasis Montaj; images are at fixed color scale and represent the same production zone at different time gates (TDEMI decay rates).

84 technology first spawned by SAIC and NRL’s 5-by-5 TDEMI sensor array. The construction and development of the R01B prototype was incremental, as several first efforts were thwarted due to resource constraints and unforeseeable hard- ware malfunctions early in the project. Ultimately, a 5-by-1 TDEMI sensor array was assembled on a metal-free cart and thoroughly tested at UIT’s New York office location as well as the two in-service testing sites. Testing involved the use of a variety of acquisition parameters and careful analysis of the acquired TDEMI data sets. UIT successfully worked with Geosoft to design and create an import template to manage the input and organization of the rich TDEMI prototype data sets produced. Initial results suggest that the system provides improved sensitivity in detecting targets but no significant improvements to horizontal target resolution compared with current TDEMI systems. Advanced processing of the TDEMI prototype data is ongo- ing in an attempt to make better sense of the acquired data. This effort involves building on existing interpretative pro- cesses, possibly by modifying and implementing target utility characterization and depth estimation techniques that use a combination of modeling and inversion routines created by SAIC for the munitions response industry. First attempts have been marginally successful with the sample data sets tested. Although the data clearly indicate the presence of buried utili- ties, it is not yet possible to accurately invert for target param- eters. This is primarily due to two factors. First, the SNR is low; it is high enough for detection purposes, but larger signals are needed for stable inversions. Second, the prototype’s timing is not disciplined enough to support continuous data collections. This can be observed in Figure 3.90. The time step between data points—defined as BlockT × nTransmitters × nStk—is 1/30 × 5 × 1 = 0.16667, assuming continuous data collection with all five transmitters firing. If a bit extra is added for calculations, the time step is expected to be ~0.2 s. When the CPUms (i.e., milliseconds as counted by the CPU) in the data samples is observed, however, the time step is seen to vary—even within a single data file. This is problematic for inversions, because the data measurements must be positioned accurately relative to each other, and that requires accurately matching up the time stamp of the TDEMI data with a time stamp from the GNSS. Through discussion with the designer and fabricator of the unit, the project team learned that this system has hardware limitations that result in a variable time step. These limitations Figure 3.89. Viewing TDEMI prototype data with 3-D visualization and in conjunction with geo-referenced utility records data.

85 were not an issue with the units used for munitions classifica- tion because those systems collect data in static mode. Only during the past few years has the TDEMI system’s fabricator learned to improve the hardware (and controlling software) such that the units can be used in continuous data collection mode. This problem has been corrected by the system fabrica- tor and will not be a problem in future units. In-Service testing Two R01B prototypes (TDEMI system and software enhance- ments for multichannel GPR system) were demonstrated and evaluated during the in-service testing in July 2012. These tools were operated by the research team on actual, current highway projects; and their deployment was observed by an authorized Subsurface Utility Engineering (SUE) firm, So- Deep, Inc., that had previously performed utility mapping at the sites. The data collected by the research team could not be processed in real time, and utilities could not be visualized in the field. All that could be done in the field was confirm that raw data were collected. The data processing and analysis were done in the office, observed by the same SUE firm. The goal of the in-service testing was to determine if the developed tools can image utilities at accurate depths and elevations, detect utilities that cannot be detected with current methods, reduce the effort needed to find utilities, and confirm the detection capabilities of current methods. On completion of the in- service testing, the SUE firm provided an evaluation report compiling its observations. This report, Field Evaluation of Tools Developed in the SHRP 2 R01B and R01C Projects (Ans- pach and Skahn 2014), is available at http://www.trb.org/ main/blurbs/171470.aspx. The selected project sites and SUE firm evaluators for in- service testing were • Virginia Department of Transportation (DOT) at Stringfel- low Road, Fairfax, Virginia; So-Deep, Inc. (see Figure 3.91); and • Georgia DOT at Talbotton Road, Columbus, Georgia; So-Deep, Inc. (see Figure 3.92). In-Service Testing Objectives and Methods The SHRP 2 R01B and R01C project teams completed the overall in-service testing project at multiple sectors along production routes designated at both the Virginia and Geor- gia sites. For the Virginia site, data acquisition began at the intersection of Stringfellow Road and Lee Jackson Memorial Figure 3.90. Example TDEMI prototype data time steps.

86 Highway (US-50) and moved south along the production route to specific areas representative of diverse underground utility conditions as determined by the project team, DOT, and SUE firm. For the Georgia site, data acquisition was con- ducted in a similar manner, beginning at the intersection of Talbotton Road and Hamilton Road and moving eastward and westward along the production route. At each location, Day 1 required the establishment of sur- vey control (Figure 3.93), a visual site assessment, and thor- ough review of the existing records of utility quality level validation data. The goal of this activity was to define data collection specifications; identify project logistical support requirements; discuss and validate project staging, access, and scheduling dependencies; and determine data acquisition strategies at the project areas of interest. The R01B project team deployed two geophysical systems: the multichannel GPR TerraVision II equipment (Figure 3.94) and the advanced TDEMI system prototype (Figure 3.95). These systems, working together, were considered to offer the greatest likelihood of detecting the subsurface utilities in these localized areas. The primary factors that affected the ability of these geophysical methods to detect subsurface objects include object size, mass, orientation, distance of the object from the sensors, and cultural surface conditions. For GPR detection of utilities, soil conditions need to allow for GPR pen- etration; for metallic utilities, there needs to be sufficient Figure 3.91. Stringfellow Road, Fairfax, Virginia, in-service testing site. Figure 3.92. Talbotton Road, Columbus, Georgia, in-service testing site.

87 Figure 3.93. Establishment of survey control before geophysical investigations during in-service testing. Figure 3.94. TerraVision II data acquisition during in-service testing. Figure 3.95. TDEMI prototype data acquisition during in-service testing. contrast between the metallic utility objects and surrounding materials to distinguish targets within the TDEMI data sets. The data acquisition software was unique to each prototype evaluated, and the data collected from these geophysical sys- tems were processed in separate specialized analysis software packages. For the TerraVision II, data were acquired with UIT’s Data Acquisition Shell (DAS) software. The multichan- nel GPR data were then processed, analyzed, and interpreted using UIT’s SPADE software. For the TDEMI prototype, data were acquired using EM3D software (developed by system fabricator, G&G Sciences). The TDEMI data sets were then processed, analyzed, and interpreted using Geosoft Oasis Montaj software. Both analysis software packages, SPADE and Oasis Montaj, had been recently adapted to improve and enhance the geophysical data processing routines specific to each system. Data processing and quality control of the production data followed predefined workflows, as defined in the standard operating procedures of Appendices A through D and in UIT’s in-service testing report. Daily equipment function tests were performed on all equipment to document that data acquisi- tion was performing as designed. The test regimen included (1) functional checks to ensure the position and geophysical sensor instrumentation was functioning properly before and at the end of each data acquisition session, (2) processing checks to ensure the data collected were of sufficient quality and quantity to meet the project objectives, and (3) interpre- tation checks to ensure the processed data were representative of the site conditions. Field notes were also reviewed for obser- vations such as cultural features and equipment malfunctions that could influence data quality and interpretation. In-Service Testing Evaluation of Process and Reporting The project team and SUE firm monitored and documented observations related to each system tested. Each definable fea- ture of work was evaluated for efficiency, quality, usability, and overall functionality. Data collection was the primary focus from July 9, 2012, to July 20, 2012. Geophysical data processing and interpretation for each system followed in the subsequent weeks. Once data processing was complete, target geophysical anomaly selections were made at predicted utility locations. For TDEMI data, horizontal target selections were made at precise peak instrument response locations (relative to surrounding background) that were believed to be associ- ated with subsurface utilities. For multichannel GPR data, horizontal target selections were made, along with depth esti- mates, at the interpreted target utility locations. All interpre- tations were made within the context of supplied resources (labor, existing records, etc.) that a typical SUE contractor might have available for any given project. A qualified UIT

88 geophysicist reviewed all of the data sets with regard to data quality, coverage, and validity of the target selections. All the field notes, processing logs, and quality control (QC) test result tables were preserved and made available to the evaluating SUE firm and project sponsor. Daily data were packaged in *.zip files by the date the data were collected, using the following folder structure: • Raw data; • All raw data from the field; • All raw ASCII files created from the raw data; • Processed data/daily; • All QC tests and production data *.xyz files created after processing in SPADE or Oasis Montaj; and • All packed QC tests and production data *.map files cre- ated in SPADE or Oasis Montaj. All raw field data [real-time kinematic (RTK)–GNSS] con- tained an associated time stamp. Corrections were applied for positioning offsets, instrument bias (including instrument latency), and instrument drift; in addition, final processed data were filtered or normalized. Each data file was named logically and sequentially so that the file name could be easily correlated with the project-specific naming conventions. Data within the files could be delineated into individual fields for each value reported. Values reported in data files included local, geographic, and/or projected coordinates for each mea- surement (one or more values which were the data associated with each data channel measurement) and the time stamp for each measurement. The primary goal of the R01B and R01C in-service testing was to have the SUE firm observe the application of the pro- totype tools (hardware and software) that had been devel- oped to meet actual field conditions. This observation covered not only the deployment of the tools but also the methods and means by which the data gathered were processed and depicted. After this observation, the SUE firm was asked to review interpretation results, compare results with existing records, and generate a report detailing the assessment. The SHRP 2 project team solicited input from the User Group Panel to develop relevant questions and thoughts for the SUE firm’s report. The preferred approach for the in-service test- ing was the independent SUE provider’s active participation in the assessment of the prototype tool(s), working in close coordination with UIT. In-Service Testing SUE Firm Summary The SUE firm report—Observation, Evaluation, and Report on the SHRP 2 R01B and R01C Tools (Anspach and Skahn 2014)—states that the TerraVision II and TDEMI prototype appear to be good tools for certain projects and may enhance, but not replace, traditional utility mapping methods. It was determined that significant further testing was warranted: a comprehensive test hole program would be beneficial going forward, especially for determination of unknowns, reliabil- ity of depths, and identification of areas of anomalies. Also, a concern arose that DOTs may be reluctant to increase their utility mapping budgets to accommodate the costs asso- ciated with the new tools. One solution was proposed: some of the utility mapping costs of the new tools could be offset with funds budgeted from other DOT departments that would benefit from the additional data the prototype tools provide. These might include but are not limited to paving and maintenance functions, archeological surveys, environ- mental surveys, geotechnical base-lining and bore hole place- ment development, arborist (for historical tree root determinations), septic field mapping, limits of cemeteries, and reduction of unknown or differing site conditions for construction departments. Tables 3.7 and 3.8 summarize the results of the R01B tools compared with the existing SUE records for each of the in- service testing sites. In-Service Testing Results and Experience The in-service testing for this project was done on a “stand alone” basis: the two main system components (on which development occurred and which were far enough along for testing) were tested on two sites where other work has been carried out by So-Deep. So-Deep used standard pipe and cable locators and test holes to map the utilities at the sites, Table 3.7. So-Deep QL-B Versus R01B Tool Utility Footage—Stringfellow Road, Virginia Utility Results (Anspach and Skahn 2014) Utility Water Gas Telephone Electric Cable TV Sanitary Storm Unknowna So-Deep QL-B footage 1,860 ft 2,795 ft 1,695 ft 75 ft NA NA NA NA R01B footage 1,740 ft 825 ft 205 ft 0 NA NA NA 210 ft Percent footage found by R01B 93% 33% 12% 0% NA NA NA NA Note: NA = not available. a Out of So-Deep scope or newly installed lines.

89 Table 3.8. So-Deep QL-B Versus R01B Tool Utility Footage—Talbotton Road, Georgia Utility Results (Anspach and Skahn 2014) Utility Water Gas Telephone Electric Cable TV Sanitary Storm Unknowna So-Deep QL-B footage 2,635 ft 1,655 ft 1,405 ft 90 ft 85 ft 1,815 ft 220 ft 880 ft R01B footage 1,855 ft 860 ft 135 ft 0 ft 0 ft 610 ft 290 ft 775 ft Percent footage found by R01B 70% 51% 9% 0% 0% 33% 100%+ na Note: na = not applicable. a Data for Unknown (unknown utility or instrument response) is shown for comparison purposes. The unknown lines found by So-Deep do not necessarily coincide with the unknown lines found by the new tools; therefore, percentages of lines found do not apply. then the new systems were brought in for a comparison study. Those comparisons are shown in Tables 3.6 and 3.7. The main outcomes from the testing showed that in conditions that were optimal for each system, each system performed reasonably well. In less than optimal conditions, such as clays soils for GPR, the system didn’t necessarily perform well. This was the most straightforward method of performing an hon- est test, but it was not really a test of the “system” approach to performing multisensor utility mapping. Essentially one set of system components were tested against a separate set of system components, which ultimately is nonsensical. The intention of multisensor systems is to perform the mapping by using all available system components in concert with each other, not competitively. In the normal course of conducting projects, the 3-D GPR system would not have been deployed at the Virginia site because of poor GPR soil conditions, but the other system components would have been deployed. If the costs of the So-Deep part and the R01B parts of these projects were added up, the resulting project costs would mostly likely be unpalatable to the DOTs involved (this was not done in this case). The reason that multisensor mapping is not more deeply in the market now is that both consultants and DOTs look at the components as separate items. In most cases, the “extra” cost of doing the added geophysical work is not included. If projects are scoped to include the system approach, the benefits of the multisensor geophysical data become evident. The benefits come from finding unknown utilities (which happens much more often in practice than on the two projects studied in R01B), limited but targeted use of test holes, better 3-D information, enhanced project and public safety from having better mapping, less public inconvenience due to road closures, and environmental benefits from fewer test holes and errant digs during construction. Multisensor mapping has been done on numerous projects by UIT before and since the R01B project and has measurable benefits. See Young and Keaton (2014) for a project example as well as a description of how enhanced utility mapping with multisensor systems fits well into sustainable engineering practice.