<%BANNER%>

Global Inversion Technique for Geotechnical Engineering Applications

Permanent Link: http://ufdc.ufl.edu/UFE0042099/00001

Material Information

Title: Global Inversion Technique for Geotechnical Engineering Applications
Physical Description: 1 online resource (135 p.)
Language: english
Creator: Tran, Khiem
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2010

Subjects

Subjects / Keywords: global, inversion, traveltimes, waveforms
Civil and Coastal Engineering -- Dissertations, Academic -- UF
Genre: Civil Engineering thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: Standard methods using seismic waves, which are routinely used for shallow subsurface investigation in engineering, have limitations in characterizing challenging profiles that include low-velocity layers and embedded cavities. The limitations are often due to insensitivity of data used for inversion and pitfalls of local inversion schemes employed in these methods. This research focuses on overcoming these limitations by developing two new methods using both sensitive data and a global inversion scheme, simulated annealing. The first method is an inversion technique to invert travel times for a wave velocity profile. The technique is based on an extremely fast finite-difference solution of the Eikonal equation to compute the first-arrival time through the velocity models by the multi-stencils fast marching method. The core of the simulated annealing, the Metropolis sampler, is applied in cascade with respect to shots to significantly reduce computer time. In addition, simulated annealing provides a suite of final models clustering around the global solution and having comparable least-squared error to allow determining uncertainties associated with inversion results. The capability of this inversion technique is tested with both synthetic and real experimental surface data sets. The inversion results show that this technique successfully maps 2-D velocity profiles with high variation. The inverted wave velocity from the real data appears to be consistent with cone penetration test (CPT), geotechnical borings, and standard penetration test (SPT) results. Employed for site characterization of deep foundation design, the developed technique is applied to combined surface and borehole data. Using the combined travel time data, the technique enables to provide credible information of material at the socket and partially detect anomalies near the socket. This becomes very important because the material at and near the socket often carries a majority of load from foundations. The inversion results of the combined data, including inverted profiles and associated uncertainties, enable to characterize spatial variability in geotechnical engineering physical parameters of subsurface formations useful in the design of deep foundations. This will be particularly useful in implementing the new LRFD design methodology that can explicitly account for spatial variability and uncertainty in design parameters. The second method is an inversion technique to invert full waveforms for a wave velocity profile. The full waveform inversion scheme is based on a finite-difference solution of 2-D elastic wave equation in the time distance domain. The strength of this approach is the ability to generate all possible wave types (body waves and surface waves, etc.) and thus to simulate and accurately model complex seismic wave fields that are then compared with observed data to infer complex subsurface properties. The capability of this inversion technique is also tested with both synthetic and real experimental data sets. The inversion results from synthetic data show the ability of detecting reverse models that are hardly detected by traditional inversion methods that use only dispersion property of Rayleigh waves. The inversion result from the real data is generally consistent with the cross-hole, SPT N-value, and material log results.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by Khiem Tran.
Thesis: Thesis (Ph.D.)--University of Florida, 2010.
Local: Adviser: Hiltunen, Dennis R.
Electronic Access: RESTRICTED TO UF STUDENTS, STAFF, FACULTY, AND ON-CAMPUS USE UNTIL 2011-08-31

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2010
System ID: UFE0042099:00001

Permanent Link: http://ufdc.ufl.edu/UFE0042099/00001

Material Information

Title: Global Inversion Technique for Geotechnical Engineering Applications
Physical Description: 1 online resource (135 p.)
Language: english
Creator: Tran, Khiem
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2010

Subjects

Subjects / Keywords: global, inversion, traveltimes, waveforms
Civil and Coastal Engineering -- Dissertations, Academic -- UF
Genre: Civil Engineering thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: Standard methods using seismic waves, which are routinely used for shallow subsurface investigation in engineering, have limitations in characterizing challenging profiles that include low-velocity layers and embedded cavities. The limitations are often due to insensitivity of data used for inversion and pitfalls of local inversion schemes employed in these methods. This research focuses on overcoming these limitations by developing two new methods using both sensitive data and a global inversion scheme, simulated annealing. The first method is an inversion technique to invert travel times for a wave velocity profile. The technique is based on an extremely fast finite-difference solution of the Eikonal equation to compute the first-arrival time through the velocity models by the multi-stencils fast marching method. The core of the simulated annealing, the Metropolis sampler, is applied in cascade with respect to shots to significantly reduce computer time. In addition, simulated annealing provides a suite of final models clustering around the global solution and having comparable least-squared error to allow determining uncertainties associated with inversion results. The capability of this inversion technique is tested with both synthetic and real experimental surface data sets. The inversion results show that this technique successfully maps 2-D velocity profiles with high variation. The inverted wave velocity from the real data appears to be consistent with cone penetration test (CPT), geotechnical borings, and standard penetration test (SPT) results. Employed for site characterization of deep foundation design, the developed technique is applied to combined surface and borehole data. Using the combined travel time data, the technique enables to provide credible information of material at the socket and partially detect anomalies near the socket. This becomes very important because the material at and near the socket often carries a majority of load from foundations. The inversion results of the combined data, including inverted profiles and associated uncertainties, enable to characterize spatial variability in geotechnical engineering physical parameters of subsurface formations useful in the design of deep foundations. This will be particularly useful in implementing the new LRFD design methodology that can explicitly account for spatial variability and uncertainty in design parameters. The second method is an inversion technique to invert full waveforms for a wave velocity profile. The full waveform inversion scheme is based on a finite-difference solution of 2-D elastic wave equation in the time distance domain. The strength of this approach is the ability to generate all possible wave types (body waves and surface waves, etc.) and thus to simulate and accurately model complex seismic wave fields that are then compared with observed data to infer complex subsurface properties. The capability of this inversion technique is also tested with both synthetic and real experimental data sets. The inversion results from synthetic data show the ability of detecting reverse models that are hardly detected by traditional inversion methods that use only dispersion property of Rayleigh waves. The inversion result from the real data is generally consistent with the cross-hole, SPT N-value, and material log results.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by Khiem Tran.
Thesis: Thesis (Ph.D.)--University of Florida, 2010.
Local: Adviser: Hiltunen, Dennis R.
Electronic Access: RESTRICTED TO UF STUDENTS, STAFF, FACULTY, AND ON-CAMPUS USE UNTIL 2011-08-31

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2010
System ID: UFE0042099:00001


This item has the following downloads:


Full Text

PAGE 1

GLOBAL INVERSION TECHNIQUE FOR GEOTECHNICAL ENGINEERING APPLICATIONS By KHIEM TAT TRAN A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2010 1

PAGE 2

2010 Khiem Tat Tran 2

PAGE 3

To my parents, whose lifetimes of hard work have made mine easier 3

PAGE 4

ACKNOWLEDGMENTS First of all, I would like to express my sin cere gratitude to my advisor, Dr. Dennis R. Hiltunen. His valuable support, encouragement, and guidance during my doctoral research were what made this possible. I have benefited tremendously from his experience and research methodology. He has provided for me not only his knowledge but also his kind consideration. I thank the other members of my committee, Dr Reynaldo Roque, Dr. Michael McVay, and Dr. Ashok V. Kumar for their valuable comments and suggestions. I would like to thank my parents for encouragi ng my studies. I deeply thank my wife who has supported my decisions and the results of thos e decisions for the past few years, and my son who has given me emotion to overcome all difficulti es of my studies. Lastly, I extend thanks to all of my friends who treated me like family. 4

PAGE 5

TABLE OF CONTENTS page ACKNOWLEDGMENTS...............................................................................................................4 LIST OF FIGURES.........................................................................................................................7 ABSTRACT...................................................................................................................................10 CHAPTER 1 INTRODUCTION................................................................................................................. .12 1.1 Problem Statement.......................................................................................................... ..12 1.2 Hypothesis........................................................................................................................14 1.3 Objectives.........................................................................................................................14 1.4 Scope.................................................................................................................................14 1.5 Organization of Dissertation.............................................................................................15 2 SITE CHARACTERIZATION MET HODS USING SEISMIC WAVES.............................17 2.1 Seismic Waves..................................................................................................................17 2.2 Site Characterization Methods Using Seismic Waves......................................................18 2.2.1 Goals of Seismic Wave Methods...........................................................................18 2.2.2 Borehole Methods..................................................................................................18 2.2.3 Surface Methods.....................................................................................................19 2.2.3.1 Methods using surface travel times..............................................................19 2.2.3.2 Methods using surface wave dispersion.......................................................20 2.3 Overview of the Surface Methods....................................................................................21 2.3.1 Overview of Methods Using Surface Travel Times...............................................21 2.3.2 Overview of Methods Using Surface Wave Dispersion.........................................23 2.4 Limitations of the Standard Methods Using Seismic Waves...........................................25 2.4.1 Limitations due to In sensitivity of Data.................................................................25 2.4.2 Limitations due to Local Inversion Techniques.....................................................26 2.5 Suggested Improvements..................................................................................................26 2.5.1 Improvements by Using Sensitive Data.................................................................26 2.5.2 Improvements by Employing Global Inversion Techniques..................................27 3 INVERSION OF FIRST-ARRIVAL TI MES USING SIMULATED ANNEALING...........33 3.1 Introduction............................................................................................................... ........33 3.2 Methodology.....................................................................................................................34 3.2.1 Forward Modeling..................................................................................................34 3.2.2 Optimization Method..............................................................................................35 3.3 Applications......................................................................................................................37 3.3.1 Applications on Synthetic Data..............................................................................37 5

PAGE 6

3.3.2 Applications on Real Test Data..............................................................................40 3.4 Conclusion for Chapter 3..................................................................................................4 7 4 INVERSION OF COMBINED BORE HOLE AND SURFACE TRAVEL TIMES..............56 4.1 Introduction............................................................................................................... ........56 4.2 Inversion of Synthetic Data..............................................................................................58 4.2.1 Synthetic Model 1...................................................................................................58 4.2.1 Synthetic Model 2...................................................................................................61 4.2.1 Synthetic Model 3...................................................................................................62 4.3 Inversion of Real Test Data..............................................................................................63 4.3.1 Newberry Test Site.................................................................................................63 4.3.2 Ocala Test Site........................................................................................................65 4.4 Conclusion for Chapter 4..................................................................................................6 6 5 TWO-DIMENSIONAL INVERSION OF FULL WAVEFORMS USING SIMULATED ANNEALING.................................................................................................80 5.1 Introduction............................................................................................................... ........80 5.2 Methodology.....................................................................................................................82 5.2.1 Forward Modeling..................................................................................................82 5.2.2 Optimization Method..............................................................................................86 5.3 Applications......................................................................................................................88 5.3.1 Applications on Synthetic Data..............................................................................88 5.3.1.1 One-dimensional synthetic models..............................................................88 5.3.1.2 Two-dimensional synthetic models..............................................................91 5.3.2 Applications on Real Test Data..............................................................................95 5.3.2.1 TAMU test site.............................................................................................95 5.3.2.2 Newberry test site.........................................................................................97 5.4 Conclusion for Chapter 5................................................................................................100 6 CLOSURE...................................................................................................................... ......128 6.1 Summary of Findings.....................................................................................................128 6.1.1 Inversion Technique Using Travel Times............................................................128 6.1.2 Inversion Technique Using Full Waveforms.......................................................129 6.2 Conclusions.....................................................................................................................130 6.3 Recommendations...........................................................................................................131 LIST OF REFERENCES.............................................................................................................132 BIOGRAPHICAL SKETCH.......................................................................................................135 6

PAGE 7

LIST OF FIGURES Figure page 2-1 Particle motions associated with body waves and surface waves (from Bolt 1976).........29 2-2 Bore-hole methods: a) cross-hole test and b) down-hole test (from Foti 2000)................30 2-3 Seismic refraction tomography: a) test set up, b) travel times, and c) inverted velocity tomogram...........................................................................................................................31 2-4 SASW: a) test setup, b) dispersion da ta, and c) inverted velocity profile.........................32 3-1 Synthetic model 1: a) true model, b) i nverted model, and c) uncertainty associated with the inverted model.....................................................................................................49 3-2 Synthetic model 2: a) true model, b) i nverted model, and c) uncertainty associated with the inverted model.....................................................................................................50 3-3 Synthetic model 3: a) true model, b) i nverted model, and c) uncertainty associated with the inverted model.....................................................................................................51 3-4 Comparison between the estimated and obser ved first-arrival time for the real data A_0-36.6............................................................................................................................52 3-5 Comparison between the estimated and obser ved first-arrival time for the real data A_36.6-73.2.......................................................................................................................52 3-6 Inversion result for the r eal data A_0-36.6: a) inverted model, and b) uncertainty associated with the inverted model....................................................................................53 3-7 Inversion result for the r eal data A_36.6-73.2: a) inverted model, and b) uncertainty associated with the inverted model....................................................................................54 3-8 Cone penetration tip resistance, line A..............................................................................55 4-1 Synthetic model 1: a) true model, b) inverted model of the surface data, c) uncertainty associated with the inverted m odel in figure b, d) inverted model of the combined data, and e) uncertainty associat ed with the inverted model in figure d...........69 4-2 Synthetic model 2: a) true model, b) inverted model c) uncertainty associated with the inverted model............................................................................................................. .70 4-3 Synthetic model 3: a) true model, b a nd c) inverted model and uncertainty using 2 shots, d and e) inverted model and uncertain ty using 4 shots, f and g) inverted model and uncertainty using 6 shots.............................................................................................72 4-4 Comparison between the observed and estim ated first-arrival times for Newberry surface data........................................................................................................................72 7

PAGE 8

4-5 Comparison between the observed and estim ated first-arrival times for Newberry combined data: a) surface data, and b) borehole data........................................................73 4-6 Newberry surface data: a) inve rted model, and b) uncertainty..........................................74 4-7 Newberry combined data: a) i nverted model, and b) uncertainty......................................75 4-8 Comparison between the observed and estim ated first-arrival time for Ocala surface data.....................................................................................................................................76 4-9 Comparison between the observed and es timated first-arrival time for Ocala combined data: a) surface data, and b) borehole data........................................................77 4-10 Ocala surface data : a) inverted model, and b) uncertainty.................................................78 4-11 Ocala combined data: a) inve rted model, and b) uncertainty............................................79 5-1 Discretization of the medium on a staggered grid...........................................................102 5-2 Comparison between wave fields genera ted by FDM and FEM: a) input model shear velocity Vs, and b) wave fields at 2 points x=30m and x=60m on the surface...............102 5-3 Shear wave velocity of 1-D models of Tokimatsu et. al. (1992): a) model 1 and b) model 2.............................................................................................................................103 5-4 Impulsive load: a) triangle wa velet and b) Ricker wavelet..............................................103 5-5 Synthetic model 1: particle velocity data from Plaxis2D observed at 10 points (A to F) located at 5 m spacing on the surface..........................................................................104 5-6 Synthetic model 1: comparison between estimated and observed wave fields...............105 5-7 Synthetic model 1: velocity hi stograms from all accepted models..................................106 5-8 Synthetic model 1: thickness histograms from all accepted models................................107 5-9 Synthetic model 1: comparison betw een the true and inverted models...........................108 5-10 Synthetic model 2: particle velocity data from Plaxis2D observed at 10 points (A to F) located at 5 m spacing on the surface..........................................................................109 5-11 Synthetic model 2: comparison between estimated and observed wave fields...............110 5-12 Synthetic model 2: velocity hi stograms from all accepted models..................................111 5-14 Synthetic model 2: comparison betw een the true and inverted models...........................113 5-15 Synthetic model 3: a) true model, and b) inverted model................................................114 8

PAGE 9

5-16 Synthetic model 3: comparison between estimated and observed wave fields...............115 5-17 Synthetic model 3: velocity hi stograms from all accepted models..................................116 5-18 Synthetic model 4: a) true model, and b) inverted model................................................117 5-19 Synthetic model 4: comparison between estimated and observed wave fields...............118 5-20 Synthetic model 4: velocity hi stograms from all accepted models..................................119 5-21 Synthetic model 5: a) true model, and b) inverted model................................................120 5-22 Synthetic model 5: comparison between observed and estimated wave fields...............121 5-23 Synthetic model 5: velocity hi stograms from all accepted models..................................121 5-24 Normalized spectral of the measured data at TAMU......................................................122 5-25 Comparison between observed wave field and estimated wave field for the real data set at TAMU....................................................................................................................122 5-26 Velocity histograms from all accepted models for the real data set at TAMU................123 5-27 Inversion result of shear wave velocity (m /s) profile for the real data set at TAMU......124 5-28 Composite soil profile of TAMU.....................................................................................124 5-29 Comparison between observed wave field and estimated wave field for the real data set at Newberry................................................................................................................ 125 5-30 Velocity histograms from all accepted models for the real data set at Newberry...........126 5-31 Newberry test site: a) S-wave tomogram from the full wave field, and b) P-wave tomogram from refraction travel times............................................................................127 9

PAGE 10

Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy GLOBAL INVERSION TECHNIQUE FOR GEOTECHNICAL ENGINEERING APPLICATIONS By Khiem Tat Tran August 2010 Chair: Dennis R. Hiltunen Major: Civil Engineering Standard methods using seismic waves, which are routinely used for shallow subsurface investigation in engineering, have limitations in characterizing challenging profiles that include low-velocity layers and embedded cavities. The limita tions are often due to insensitivity of data used for inversion and pitfalls of local inve rsion schemes employed in these methods. This research focuses on overcoming these limitati ons by developing two new methods using both sensitive data and a global invers ion scheme, simulated annealing. The first method is an inversion technique to invert travel times for a wave velocity profile. The technique is based on an extremely fast fin ite-difference solution of the Eikonal equation to compute the first-arrival time through the veloci ty models by the multi-stencils fast marching method. The core of the simulated annealing, the Metropolis sampler, is applied in cascade with respect to shots to significantly reduce computer time. In addi tion, simulated annealing provides a suite of final models clustering around the glob al solution and having comparable least-squared error to allow determining uncertainties associated with inversion results. The capability of this inversion technique is tested w ith both synthetic and real expe rimental surface data sets. The inversion results show that this technique successfully maps 2-D velocity profiles with high 10

PAGE 11

variation. The inverted wave velocity from the r eal data appears to be consistent with cone penetration test (CPT), geotech nical borings, and sta ndard penetration test (SPT) results. Employed for site characterization of deep f oundation design, the developed technique is applied to combined surface and borehole data. Using the combined travel time data, the technique enables to provide credible information of material at the sock et and partially detect anomalies near the socket. This becomes very important because the material at and near the socket often carries a majority of load from f oundations. The inversion re sults of the combined data, including inverted profiles and associated uncertainties, en able to characterize spatial variability in geotechnical engineering physical parameters of subsurface formations useful in the design of deep foundations. This will be par ticularly useful in implementing the new LRFD design methodology that can explicitly account for sp atial variability and uncertainty in design parameters. The second method is an inversion technique to invert full waveforms for a wave velocity profile. The full waveform inversion scheme is based on a finite-difference solution of 2-D elastic wave equation in the time distance domain. The strength of th is approach is the ability to generate all possible wave types (body waves and surface waves, etc.) and thus to simulate and accurately model complex seismic wave fields that are then compared with observed data to infer complex subsurface properties. The capability of this inversion technique is also tested with both synthetic and real experimental data sets. The inversion results from synthetic data show the ability of detecting reverse mode ls that are hardly detected by traditional inversion methods that use only dispersion property of Rayleigh waves. The inversion result from the real data is generally consistent with th e cross-hole, SPT N-value, and material log results. 11

PAGE 12

CHAPTER 1 INTRODUCTION 1.1 Problem Statement Near surface soil conditions control the responses of f oundations and structures to earthquake and dynamic motions. To get the op timum engineering design, the modulus of underlying layers must be determined correc tly. A popular method to obtain modulus is nondestructive in situ testing. The measured elas tic wave fields carry substantial information about characteristics of media they propagate in, and wave evaluation techniques use this information to infer the properties of media. Th ere are many techniques that are routinely used for shallow subsurface investigation, such as techniques using wave velocity dispersion, including spectral analysis of surface waves (SASW) (Nazarian, 1984), multi -channel analysis of surface waves (MASW) (Park, et al., 1999), refrac tion microtremor (ReMi) (Louie, 2001), and passive-source frequency-wavenumber (f-k), and techniques us ing travel times, including conventional seismic refraction and seismic re fraction tomography (SRT). However, these standard techniques utilized in current engineering evaluations have limitations in dealing with challenging profiles, which include anomalies su ch as low-velocity z ones and cavities. The limitations of these methods are generally caused by: 1) insensitivity of data used for inversion, and 2) shortcomings of local inversion t echniques often employed in these methods. Regarding the insensitivity of data, standa rd techniques usually use only a limited portion of information carried by elastic wave fields, such as travel times or wave velocity dispersion. Synthetic tests in this research prove that surf ace travel times are not very sensitive to the anomalies. Thus, the anomalies are not well detected if using only surface tr avel times. For the techniques using wave velocity dispersion, ONeill (2003) found th at a low-velocity layer (LVL) of few meters thickness buried at depth of a few meters can not be inverted with confidence, due 12

PAGE 13

to loss of sensitivity with depth. In addition, the dispersion propert y is absolutely insensitive to lateral variation, thus techniques using wave ve locity dispersion are only for 1-D problems and can not be applied in regions with high lateral variation. Regarding the shortcomings of local inversi on techniques, geotechnical inversion problems are typically highly non-linear, non-unique, and influenced by inherent noise in measured data. Results from any local invers ion technique heavily depend on the initial model and priori information. In order to obtain a good inversion result, these local techniques require both a reasonable initial model and priori information th at are not always available. In addition, it does not guarantee that the model corresponding to the smallest error (best fitting) is closest to the true model because of the noise and non-uniqueness of inverted solutions. Any of many thousand models that have a similar degree of fitting can be the inversion result. The question here is how to decide which model is the best and the uncertainty associated with that model. This research focuses on overcoming the aforementioned limitations by developing new methods using both sensitive data and global inve rsion schemes. Firstly, data sensitive to anomalies, such as borehole data will be utilized to improve resolution of travel time inversion results. By using borehole data it is possible to increase depth of investig ation, to reduce the uncertainty associated with inversion result, a nd to improve resolution of inverted profiles. Secondly, full waveform data will be also utilized to solve 2-D reverse problems when the borehole data is not available. The reason to us e the full waveform data instead of the wave velocity dispersion is that it is more sensitive to the anomalies. Lastly, global optimization schemes based on simulated annealing will be employed in both travel time and full waveform inversions. These global technique s do not depend on the initial model and become important in regions where no prior information about subsurface profiles is available. 13

PAGE 14

1.2 Hypothesis Global inversion techniques using sensitive da ta such as borehole travel times or full waveforms can characterize subsurface profiles be tter than traditional te chniques that use only surface travel times or dispersion property. 1.3 Objectives The primary objective of the research is to develop inversion tec hniques that can well characterize challenging profiles with requirement of a little or no priori information, and can also evaluate the uncertainties associated with the inverted profiles. Sp ecific objectives of this research include the following: Develop a global inversion technique to invert travel times for a velocity structure. The technique is independent of initial model, and thus can be applied in regions where no priori information is available. Determine the uncertainty associated with th e inverted profile from travel times. This uncertainty represents the influence of noi se in measured data and non-uniqueness of inversion solution. Incorporate borehole data into a travel time inversion to improve the resolution of the inverted profile and reduce th e associated uncertainty. Develop another global inversion technique using a full waveform. This technique is primarily used when borehole data is not av ailable, particularly for cases of reverse velocity structures. 1.4 Scope Non-destructive tests including refraction test, full waveform test and/or invasive tests including con penetration test (CPT), geotechnical borings, sta ndard penetration test (SPT), and cross-hole test were conducted at three test sites: A National Geotechnical Experi ment site (NGES) at Texa s A & M University (TAMU). A Florida Department of Transportation (F DOT) storm water runoff retention basin in Alachua County off of state ro ad 26, Newberry, Florida. A test site near Ft. McCoy, Ocala, Florida. 14

PAGE 15

Two inversion techniques, which use travel times and full waveforms, have been developed and intensively tested on synthetic models before bei ng applied to real test data. For the technique using travel times, it is first tested on many synthe tic models that were previously tested on three commercial refraction tomography codes by Sheehan et al. (2005) to verify its capability, and then applied to many real data sets. The inverted profiles from non-destructive data are compared to CPT, SPT, and geotechnical borings results to appraise its feasibility for geotechnical engineering applications. Similarly, the technique using full waveforms is first tested on challenging 1-D and 2-D synthetic models that include lowand high-velocity layers to ve rify its capability of dealing with reverse profiles, and then applied to real te st data. The inverted resu lts from the real test data are also compared to SPT, cross-hole test, a nd refraction test results for partially appraising the accuracy of those inverted results. 1.5 Organization of Dissertation An overview of the following chapters follows. Chapter two provides an overview of site characterization methods using seismic waves, limitations of these standard methods, and suggested improvements. Chapter three presents a global inversion technique to invert travel times using simulated annealing. The technique is described in de tails including forward modeling, optimization method, and its applications to synthe tic and real test data sets. In this chapter, only travel times measured on the surface are used for inversion, and inverted results are compared to invasive test results to appraise the capability of the inversion technique. Chapter four demonstrates the coupling of so-called down-hole and refraction tomography technique, which uses combined bore-hole and su rface travel times. Comparisons of inversion results utilizing the combined data against results developed using just the surface data are 15

PAGE 16

presented to both quantitatively a nd qualitatively appraise the be nefit when adding data from a borehole. The focus here is to detect voids in rock that may be filled by air, water or soil. Chapter five presents another inversion technique to invert full waveforms using simulated annealing. This technique is also described in details including forward modeling, optimization method, and its applications to synt hetic and real test data sets. Th e primary applications of this technique are for reverse profiles when the borehole data is not available. 16

PAGE 17

CHAPTER 2 SITE CHARACTERIZATION MET HODS USING SEISMIC WAVES 2.1 Seismic Waves A wave is a perturbation that travels in a medium and carries energy (Doyle, 1995). Waves propagate in a medium can be divided into two categories: body waves and surface waves, whose particle motions are presented in Figure 21. Body waves that propagate in the interior of a medium include longitudinal waves (P-wave) and shear waves (S-wave). P-waves generate compression or dilatational particle motions paralle l to the direction of propagation, and S-waves generate particle motions perp endicular to the direction of propagation. Unlike body waves, surface waves only propagate in a shallow zone n ear a free surface. They are essentially of two different kinds: Love waves and Rayleigh waves. Love waves only exist when a medium consists of a soft superficial zone over a stiffe r zone, and the particle motions associated with them are transversal with respect to the direct ion of propagation. Rayleigh waves always exist near a free surface, and their part icle motions consist of vertical and horizontal components that are exactly 900 out of phase, and thus the particle motion path is an ellipse. Many techniques for soil characterization at very small strain levels are based on measurements of particle motions associated to wave propagation. This is made possible because of strong links existing between wave propagation characteristics and mech anical properties of the medium that need to be characterized. Thes e links are between wave velocities and elastic moduli as: (2-1) 2 SVG )2(2 2 S PVV (2-2) )23( E (2-3) 17

PAGE 18

where VP, VS are P-wave and S-wave velocities, are Lames parameters, G, E are shear and youngs moduli, and is density. 2.2 Site Characterization Methods Using Seismic Waves 2.2.1 Goals of Seismic Wave Methods The primary goal of site characterization methods using seismic waves is to determine elastic modulus of subsurface materials. The mo dulus is used either directly in design of engineering substructures or as an indication to delineate anomalies such as voids, low-velocity zones embedded in the subsurface profiles. 2.2.2 Borehole Methods Borehole methods include seismic down-hole and cross-hole methods (Figure 2-2). Down-hole tests can perform in a single borehole. In the down-hole test, an impulsive source is located on the surface, a single receiver that can move or a string of rece ivers are fixed to the walls of the borehole to measure the wave fields from the source. A plot of travel time versus depth can be generated, and the slope of the travel time curve represents the wave velocity at depth for a 1-D soil profile. Cross-hole tests use two (three ) boreholes. A source and one (t wo) receiver are located at the same depth in the boreholes, and the wave pr opagation velocity of the material between the boreholes at that depth is meas ured. By testing at various dept hs, a 1-D velocity profile is obtained. The borehole methods can be used to obtain cred ible 1-D velocity profiles. However, these methods are not popular as surface methods becau se of the requirement of boreholes and difficulties in measurements. 18

PAGE 19

2.2.3 Surface Methods Unlike borehole methods, surface methods are attractive to resear chers because the measurements on the surface are relatively easy, and the measured signals carry important information about the mechanical properties of the medium. These methods have been intensively developed and routin ely used for geotechnical site characterization in last few decades. They are roughly categorized into tw o groups: methods using travel times and methods using wave velocity dispersion. 2.2.3.1 Methods using surface travel times Methods using surface travel times are ofte n referred as seismic refraction methods including conventional methods and seismic refraction tomography (SRT). The seismic refraction methods consist of measuring data at a test site, picking travel times from the measured data, and doing inversion to invert th e travel times for a subsurface wave velocity profile. Figure 2-3 presents a test setup, picked tr avel times, and inverted velocity tomogram of a typical seismic refraction tomography. A data measurement is conducted at known points along the surface to obtain a seismic wave field generated by an energy source. The energy source is usually generated by a small explosive charge or accelerated weight drop, such as a sledgehammer. Energy that propagates out from the source, either travel s directly through the upper layer, known as a direct arrival, or travels down to and laterally along a lower high ve locity layer, known as a refracted arrival. This energy is detected on the surface using rece ivers arranged in a linear array, and travel times are derived by manually picki ng the first-arrival signals. The travel times are then used in an inversion to generate a velocity profile of the subsurface along the length of the refraction profile relying upon the assumption that th e velocity of the subsurface material increas es with depth. 19

PAGE 20

During the inversion, while the conventional me thods divide the model into continuous layers of constant velocity, SRT divides the model into a large numbe r of smaller constant velocity grid cells or nodes. The model is ad justed each iteration in attempt to match the calculated travel times with the measured as closel y as possible. This is repeated until the predefined stop criterion is achieved, and the inverted profile is then derived. 2.2.3.2 Methods using surface wave dispersion Two most popular methods using dispersion property of measured wave fields are spectral analysis of surface waves (SASW) (Nazarian, 1984) and multi-channel analysis of surface waves (MASW) (Park, et al., 1999). They both consist of three steps: measuring data from a test site, processing the measured data to obtain a disp ersion curve, and doing inversion to invert the dispersion curve for a subsurface shear wave veloc ity profile. Figure 2-4 presents a test setup, a typical dispersion curve, and an inverted ve locity profile of SASW. For MASW, multiple geophones are used instead of two; the dispersion data and inverted velocity profile are similar to those of SASW. For the first step, a wave field from an active source is measured by two geophones (SASW) or multiple geophones (MASW). The geophones are located far enough from the source in order to get dominant components of Raylei gh waves in the measured wave field. For the second step, a cross power spectrum analysis (SASW) or a multi-channel data processing method (MASW) is applied to the measured wave fiel d to obtain dispersion property. The measured wave field is here decomposed into components of different frequencies, and phase velocity of each component is calculated to establish a dispersi on curve that is used as experimental data for an inversion in the last step. Inversion of the dispersion curv e is a process for determining a shear wave velocity profile from frequency-phase velocity dispersion relations hip. This process consists of evaluation of a 20

PAGE 21

theoretical dispersion curve for an assumed pr ofile and comparison with the experimental dispersion curve. When the theoretical disper sion curve and the experimental dispersion relatively match, the assumed profile is the desired solution. The a ssumed profile is composed of horizontal layers that are homoge neous, isotropic and the shear velocity in each layer is constant and does not vary with depth. The theoretical di spersion curve calculation is often based on the matrix formulation of wave propagation in layered media given by Thomson (1950). Besides methods using the dispersion propert y of active wave fields, there are also methods using the dispersion property of pa ssive wave fields. Among those, refraction microtremor method (ReMi) (Louie, 2001) and passive-source frequency-wavenumber (f-k) method are more popular. Passive methods also re quire three steps as t hose in active methods. First, a wave field that is ambient vibrations is measured by a linear geophones array (ReMi) or a two-dimensional geophones array (f-k). S econd, a signal processing method, which is a slowness-frequency transformation (ReMi) or fr equency-wavenumber transformation (f-k), is applied to the measured signal to obtain an expe rimental dispersion curve. Lastly, an inversion technique, which is exactly the same as that of active methods, is used to invert the experimental dispersion curve for a subsurface shear wave velocity profile. 2.3 Overview of the Surface Methods 2.3.1 Overview of Methods Using Surface Travel Times Conventional refraction methods are useful fo r a simple profile consisting of a few constant velocity layers with linear interfaces (Sheehan et al., 2005). In this case of sharp contrasts in velocity at the interfaces, the conventional met hod is more suitable than SRT because of two reasons. First, SRT will always model a sharp contrast in velocity with a gradient in velocity, and second SRT requires that more shots and receivers be used than conventional techniques. This increase in field data coll ected significantly increa ses the amount of time 21

PAGE 22

necessary to perform the field testing. The time required to process this data, i.e., make firstarrival picks, and the time required for inversion is also si gnificantly increased. Seismic refraction tomography is able to reso lve velocity gradients and lateral velocity changes within the subsurface w ith greater ability than conventional methods. As such, SRT may be applied in settings where conventiona l refraction methods fail, such as areas of compaction, karst, and zone faults (Zhang and Toksoz, 1998). Carpenter et al. (2003) used SRT programs on sites known to be karstic in Illinois and Kentucky, to determine their accuracy by comp aring with known subsurface geology from outcrops, borehole logs, other geophysical met hods, and synthetic travel time data. The tomograms seemed to indicate consistency with known geology and borehole information. However, the authors mention problems with artif acts appearing such as pinnacles or cutters in the lower third of the model when too many veloc ity elements are used. These artifacts may be due to the insensitivity of the surface measured travel times to the deep and small velocity elements. These small elements can be assigned any velocity without ch anging the surface travel times, thus the artifacts are cr eated as a product of inversion. Sheehan et al. (2005) investigated the abil ity of SRT programs fo r karst terrains that frequently contain sinkholes, irregular and gradational bedrock interfaces, and voids in rocks that may be filled by air, water or soil. They found that SRT performs well in many situations where conventional methods fail, e.g., where vertical and lateral gradients compose a significant component of the velocity struct ure, and however SRT fails to detect embedded cavities. Again, the surface travel times are not sensitive to the embedded cavities, thus these cavities are hardly detected by surface-based SRT methods. 22

PAGE 23

Hiltunen and Cramer (2008) used SRT to characterize the subsurface at Pennsylvania bridge sites in mantled karst terrain. Refraction profiles ge nerated by commercial software SeisImager were compared with pile tip elevat ions at driving refusal and drilling data. They found that SRT appears to be able to characterize the soil/rock interface. However, there exists a zone where pile tip elevations are much deeper th an the soil/rock interfac e depth from SRT. This discrepancy is due to the presence of a low-velocity zone, which SR T fails to detect by just using the surface measurements. The failure of detect ing the low-velocity zone can be possibly explained as following. The invers ion technique of SeisImager is a deterministic method that depends on an initial model, which has veloc ity increasing with dept h. During the inversion, travel times and ray paths of each assumed model are calculated, and only cells having ray paths are updated. It is likely no ray paths go through the area where the low-velocity zone supposes to be, and thus this area is not fre quently updated and has the velocity similar to that of the initial model. 2.3.2 Overview of Methods Using Surface Wave Dispersion For active surface wave methods, SASW and MA SW have become very popular in recent years for shallow non-destructive testing of both layered natural (soil and rock) and artificial (pavement and concrete) structures. They typically perform well on normal profiles that have velocity increasing with depth. However, for reverse profiles consisti ng of both high-velocity and low-velocity layers, these methods n eed to be used with particular care. ONeill (2003) intensively investigated inversion using surface wave dispersion of multiple modes on a variety of models. He found th at the most difficult pa rameters to accurately interpret with the surface wave dispersion are layers below a hi gh-velocity layer (HVL), and a deep (a few meters) buried low-velocity layer (L VL) with a thickness of a few meters can not be inverted with confidence. He also found that for a case of real data from a well characterized test 23

PAGE 24

site, a 5 m thick stiff silt-sand layer at 20 m de pth, in a mostly soft, clay background cannot be interpreted with any surface wave dispersion method, due to loss of sensitivity with depth. The main reason behind these limitations is that the su rface wave dispersion is not very sensitive to thin (relatively to depth) embedded LVL or HVL. Jin et al. (2009) studied the role of forward model in surface-wave i nversion to delineate a HVL of 1.5 m thickness, 2 m depth, and 1500 m/s sh ear velocity (Vs). Th ey found that the HVL can be detected by an inversion technique c onsisting of simulated annealing followed by a linearized inversion, using only the fundamental mode. However, even though using simulated annealing to search over a larg e parameter space at the beginni ng of inversion, it still requires relatively strict constraints for the HVL such as thickness, 0 to 2 m; depth, 1 to 5 m; Vs, 1000 to 2000 m/s. That means the strictly priori informat ion about the HVL are needed to delineate it, and the technique would fail if the pr iori information is not available. For passive surface wave methods, ReMi and f-k methods have been routinely used recently for deep shear wave velocity profiles. The most important advantage of testing methods using passive waves is the ability to obtain deep depths of inve stigation with very little field effort. Desired Rayleigh waves from passive seis mic arrivals are relatively pure plane waves at low frequencies allow determining Vs profiles up to hundreds meter depth. Li (2008) applied ReMi and f-k methods on real data sets meas ured from eleven test sites distributed over a distance of about 180 km in the upper Mississippi Embayment in the central United States to obtain Vs prof iles up to depth of hundreds of meters. The Vs profiles are compared to those from SASW and MASW, which used a new developed low-frequency vibrator to generate wave fiel ds at very low frequencies down to 1 Hz. He found that the active surface wave methods are more reliable than passive wave methods for deep Vs profiling. 24

PAGE 25

2.4 Limitations of the Standard Methods Using Seismic Waves As reviewed in section 2.3, the limitations of the standard methods using seismic waves are generally caused by insensitivity of data used for inversion and the shortcomings of local inversion techniques. 2.4.1 Limitations due to Insensitivity of Data The standard seismic wave methods using either surface travel times or surface wave dispersion interpret subsurface velocity structures The travel time from a shot to a receiver on the surface is measured from the fastest ray that starts from the shot and travels through a medium to the receiver. This fast est ray tends to go through stiffer material (high velocity zones), and avoids softer material (low velocity zones). Thus the surface travel times are not sensitive to these low velocity zones, and c onsequently these zones can not be well characterized with only travel times. Surface wave dispersion data is more sensitive to a low velocity layer than surface travel times. However, the dispersion property is de veloped by taking accounts of whole material within a depth of approximately one wavelengt h for each frequency, and more sensitive to a shallow stiff layer than to a deep soft layer. When the depth of investigation increases, lower frequency or longer wavelength co mponents are required, larger volum es of material are utilized to derive the dispersion property, and the dispersion data becomes insensitive to a deep and thin low-velocity layer. Hence, the low-velocity layer is hardly detected by just using the surface wave dispersion. In addition, the dispersion property is absolutely insensitiv e to lateral variation, thus techniques using wave dispersion are only for 1-D problems and can not be applied in regions with high lateral variation. 25

PAGE 26

2.4.2 Limitations due to Local Inversion Techniques The main goal of inversion is to find earth models that explain th e seismic observations. The inversion involves finding an optimal value of a misfit function of several variables. The misfit function characterizes the differences betw een the observed and synthetic data calculated by using an assumed earth model. The earth mo del is described by physical parameters that characterize properties of soil/ro ck layers such as compressi on and shear wave velocities, density, etc. For many geotechnical applications, the misf it function is highly complicated, and its surface consists of multiple hills and valleys. T hus such a function will have many maxima and minima; the minimum of all the minima is calle d the global minimum and all other minima are called the local minima. Local inve rsion techniques such as grad ient descent methods typically attempt to find a local minimum in the close neighborhood of the starting solution. They use the local properties of the misfit f unction (gradients) to calculate the update to the current position and search in the downhill direction that only accepts a model having a smaller misfit value. Thus, these techniques will miss the global minimum if the starting solution is nearer to one of the local minima than the global minimum. 2.5 Suggested Improvements 2.5.1 Improvements by Using Sensitive Data The capability of the seismic wave methods will be significan tly improved by using sensitive data for inversion. The target of this research is to develop techniques that can characterize anomalies such as low-velocity laye rs, voids in rocks that may be filled by air, water, or soil. For this purpose, combined travel times measured both on the surface and in a borehole will be utilized to obtain high resolution inverted profiles. In addition, full waveforms will be employed to characterize reverse profiles when the borehole data is not available. The 26

PAGE 27

full waveform has been used by authors for the deep (> 1 km) subsurface investigation, in both the time domain (Shipp and Singh [2002], Sheen et al. [2006]) and the frequency domain (Pratt [1999]). In all of these applic ations, body waves are dominant components in wave fields used for inversion. Unlike these methods, this research uses Rayleigh waves as dominant components in wave fields. Because, Rayleigh waves are more sensitive to shear wave velocity and less sensitive to Poisons ratio or de nsity of a medium, they are usef ul for geotechnical applications where the shear wave velocity primarily controls properties of soil particles. The following chapters will prove that the bo rehole travel times and full waveforms are more sensitive to the anomalies than surface travel times and wave velocity dispersion, respectively. 2.5.2 Improvements by Employing Global Inversion Techniques The capability of the seismic wave methods will be also improve d by employing global inversion techniques. Although, these techniques re quire much more computer time than local techniques, they are feasible for geotechnical engineering problems with the advent of fast computers. Unlike local inversion techniques, global i nversion techniques such as simulated annealing, genetic algorithm, and other importa nce sampling approaches attempt to find the global minimum of the misfit function by searchi ng over a large parameter space. They have recently been applied in evaluation of various geophysical data sets (Sen and Stoffa [1991, 1995], Pullammanappallil and Louie [1994], Sharma and Kaikkonen [1998]). Most of the global inversion techniques are stocha stic and use more global info rmation to update the current position, thus they likely converg e to the global minimum. In th is study, simulated annealing will be employed for inversion problems in geotechnical engineering because it can be used in cases 27

PAGE 28

where the model-data relationship is highly non linear and produces multimodal misfit functions (Sambridge and Mosegaard, 2000). 28

PAGE 29

Figure 2-1. Particle motions associated with body waves and surface waves (from Bolt 1976) 29

PAGE 30

a) b) Figure 2-2. Bore-hole methods: a) cross-hole test and b) down-hole test (from Foti 2000) 30

PAGE 31

a) b) c) Figure 2-3. Seismic refraction tomogr aphy: a) test setup, b) travel times, and c) inverted velocity tomogram 31

PAGE 32

a) b) 200 300 400 500 600 700 800 900 1000 1100 1200 010203040506070 Frequency (Hz)Phase Velocity (ft/s) c) Figure 2-4. SASW: a) test set up, b) dispersion data, and c) inverted velocity profile 32

PAGE 33

CHAPTER 3 INVERSION OF FIRST-ARRIVAL TI MES USING SIMULATED ANNEALING 3.1 Introduction This chapter presents an inversion technique to invert first-arrival time using simulated annealing. The scheme is based on an extremely fast finite-difference solution of the Eikonal equation to compute the first-arrival time through the velocity models by the multi-stencils fast marching method. The core of the simulated ann ealing, the Metropolis sampler, is applied in cascade with respect to shots to significantl y reduce computer time. In addition, simulated annealing provides a suite of final models clustering around the global solution and having comparable least-squared error to allow determ ining uncertainties associated with inversion results. Global inversion techniques such as simulated annealing, genetic algorithms, and other importance sampling approaches have recently been applied in evaluation of various geophysical data sets (Sen and Stoffa [1991, 1995], Pu llammanappallil and Louie [1994], Sharma and Kaikkonen [1998]). Unlike local inversion techni ques, global inversion techniques attempt to find the global minimum of the misfit function by searching over a large parameter space. Most of the global inversion techniques are stochastic and use more global information to update the current position, thus they likely converge to the global minimum. In this study, simulated annealing is employed for engineering problems because it can be used in cases where the model-data relationship is highly nonlinea r and produces multimodal misfit functions (Sambridge and Mosegaard [2000]). In addition, simulated annealing enables to provide a suite of final models clustering around the global solu tion and having comparable least-square error. This allows getting the inversion result by averagi ng all of these models to mitigate the influence of noise and the non-uniqueness of the inversion solutions. 33

PAGE 34

The most important advantage of global invers ion techniques is to avoid being trapped in local minima, and thus to allow final inversion results to be in dependent of the initial model. However, these global inversion techniques can require significant computer time, especially if the model contains a large number of model para meters. This disadvantage is reduced herein by using an extremely fast forwarding model solu tion, and by sampling in cascade with respect to shot position. First-arrival time inversions usually require many shots to obtain a high-resolution profile and each shot needs a forwarding m odel solution. Meanwhile, global optimization methods sample over a large number of trial mode ls, and only accept a small part of them. Using an acceptance rule in cascade, forwarding model solutions for only a few shots are often required to reject the biased models. For the cases presen ted herein, a saving in co mputer time of 70% is achieved when utilizing the acceptance rule in cascade. The technique is presented herein for syntheti c data sets created fr om models that were previously tested on three commercial refraction tomography codes by Sheehan et al. (2005), and two experimental data sets from a well documented test site. The inversion results show that this technique successfully maps 2-D velocity prof iles with high variati on. The inverted wave velocity from the real data appears to be c onsistent with cone pe netration test (CPT), geotechnical borings, and standard penetration test (SPT) results. 3.2 Methodology 3.2.1 Forward Modeling Two main approaches that have been routinely used to calculate first-arrival times are the shortest path method (SPM) and the Eikonal eq uation solution. SPM or iginated in network theory (Dijkstra [1959]) and was applied by Na kanishi and Yamaguchi (1986) and Moser (1991) to seismic ray tracing. The main advantage of SPM are its simplicity and capacity for simultaneous calculation of both th e first arrival time and the asso ciated ray path, but it takes 34

PAGE 35

more time for calculation and is not very effi cient to use for the global optimization methods. The Eikonal equation shown below has been solv ed by many authors, including (Vid ale [1988], Van Trier and Symes [1991], Nichols [1996] Sethian [1996, 1999], Ki m [1999], Chopp [2001], Hassouna and Farag [2007]). 2 2 2),( 1 zxVz u x u (3.1) In this equation, u=u(x,z) and V(x,z ) are the arrival time and veloci ty at point (x,z), respectively. Among these methods, the fast marching method (FMM) is typically considered the fastest and most stable and consistent me thod for solution of the Eikonal e quation. It was first presented by Sethian (1996) and has been improved by other au thors. In this study, the improved version of FMM by Hassouna and Farag (2007), the so called multistencils fast marching (MSFM), is utilized to compute first-arrival times for fo rward modeling. MSFM computes the solution at each point by solving the Eikonal equation along seve ral stencils, and then picks the solution that satisfies the upwind condition. 3.2.2 Optimization Method Simulated annealing was first proposed by Metropolis et al. ( 1953) and significantly improved by Kirkpatrick et al (1983). The basic concepts are borrowed from a physical annealing process. The annealing process occurs when a solid in a heat bath is initially heated by increasing the temperature such that all the part icles are distributed randomly in a liquid phase. This is followed by a slow cooling such that all particles arrange themselves in the low energy ground state where crystallization occurs. Simulated annealing has recently been applie d in evaluation of various geophysical data sets. Rather than discussing th e analogy of simulated annealing that has been well described by 35

PAGE 36

Sen and Stoffa (1995), a brief desc ription of the process used in this study is presented in the following: 1. Select an initial velocity model V and then calculate th e least-squared error Ei, for the ith shot, and average error Emean, as: N k ikik iVgd N E1 2 ,,)( 1, (3.2) M i i meanE M E11, (3.3) where dk,i and gk,i are the kth observed and computed travel times from the ith shot, respectively, N is the number of observati on points, and M is the number of shots. 2. Following the idea of Pullammanappallil and Lo uie (1994), perturb the velocity model by randomly selecting a box within the medium and randomly assign all of the cells within the box a new velocity between the minimum and maximum velocities chosen by the user. Next, to avoid unusual artifacts and high ve locity gradients that would violate the assumption of an Eikonal equation solution, the model is smoothened. The smoothing is accomplished by adjusting the velocity of each cell in contact with and just outside the boundary of the box. The adjusted velocity for e ach of these perimeter cells is established by averaging the velocities of the four cells in contact with the four edges of the perimeter cell. Following the process describe d, a new model is now obtained. 3. Apply the acceptance rule in cas cade with respect to shots for the new model. The travel times from one shot at a time are used to test the new model. For the ith shot, the leastsquared error Ei corresponding to the new model is calculated as equation (2). The new model unconditionally passes the ith shot if its least-squared erro r is smaller than that from the previous model, and conditionally passes the ith shot if its least-s quared error is bigger with the following probability: T E pi iexp, (3.4) where Ei is the positive difference between the leas t-squared errors of the two models. If the new model passes the first shot, then it is tested for the second shot and so on. It is immediately rejected if faili ng at a certain shot, i.e., the new model is only accepted if it passes all of the shots. 4. Repeat steps 2 and 3 until a desired number of accepted models are completed, then reduce the temperature to the next level according to a schedule: 2 )1( )( kT kT (3.5) where T(k) and T(k-1) are the temperatures at the kth level and the (k-1)th level. The criterion for convergence requires that the relative difference between the maximum and minimum average errors Emean from all accepted models at a particular temperature become 36

PAGE 37

very small (empirically say between 0.1% and 0. 5%). If the criterion for convergence is not satisfied, the temperature is reduced to the next level. At a sufficiently low temperature, the accepted models are in the vicinity of the global minimum and the optimization converges. The selection of a correct initial temperature is very important. Selecting a very high value leads to unnecessary computations, whereas star ting from a low value can lead to a local solution. Based on the idea of Sharma and Ka ikkonen (1998), in this study, the initial temperature is taken to be the nearest lower power of ten to the least-squared error; for example, if the least-squared error is 50.59, then the initial value is 10. In this way, the initial temperature is usually higher than the critical temperature, and any bias crea ted by the initial model will be destroyed. The cooling schedule is another important parameter in the global optimization. Many schedules have been investigated herein for i nversion problems up to 1500 cells, and it has been found that the number of accepted models at each temperature level should be at least 10,000, and the temperature needs to reach a very low va lue to obtain the global solution. In the studied cases, the final temperature was about 10-3 to 10-4 when the initial temperature was 10. When the number of cells increases, the pr ocedure needs a slower schedule, i.e., more accepted models at each temperature or a smaller change in temp erature, and hence the time of computation increases in order to get the gl obal solution. Therefore, the cooling schedule should be selected with care. 3.3 Applications 3.3.1 Applications on Synthetic Data Synthetic models refer to earth models whose velocity profile is assumed or known a priori. Using a known velocity structure, travel time data is calculated for an assumed test 37

PAGE 38

layout. This travel time data is then input to the inversion program as if it were data collected from a field test, and velocity st ructure is calculated from travel time data. Theoretically, the interpreted velocity profile should be the same as the model assumed at the start. Synthetic model studies were c onducted for at least two reasons. The first objective was to assess ability of the inversion t echnique to identify and delineate subsurface features that are routinely present in test sites. It was important to determine prior to actu al field testing whether the system is sensitive to these features, otherwise there was no point in conducting an actual field test. Second, in conducting model studies in which the answer is known, an analysis protocol can be developed to systematically and consistently analyze travel time data to develop believable velocity tomograms. In order to be confident of results interpreted from actual field data, one has to develop confidence that the pr ogram and analysis procedures are reputable. Synthetic model studies he lp build this confidence since inte rpreted velocity profiles can be compared to the known starting profile. The technique presented was intensively test ed on many synthetic models, and three of them with high lateral variation are presented here in in Figure 3-1. These three synthetic models (1, 2, and 3) are models 3, 5, and 6 from Sheehan et al. (2005), respectiv ely. They include a notch (model 1), a stair step (model 2), and broa d epikarst (model 3). For all three models, 25 shots into 48 geophones were used to create synthetic first-ar rival times. The geophone spacing was 1 m for models 1 and 2, and 2 m for model 3. For all models, all 25 shots were within the geophone spread and the shot sp acing was twice the geophone spacing. With the so arrangement, the number of shots, the number of geophones, a nd the geophone spacing are the same as those used in Sheehan et al. (2005). 38

PAGE 39

To run the inversion, the medium sizes were kept the same as the true models, and discretized with a grid size of 1 m for ease in picking the firstarrival times from the forward modeling. The velocity constraints were assigned by slightly increasing th e maximum true value and decreasing the minimum true value, and thus the velocity was allo wed to vary between 100 m/s and 4000 m/s for models 1, and 2, and betw een 100 m/s and 5000 m/s for model 3. All initial models were specified with a cons tant velocity of 2000 m/s, which leads to an ini tial temperature of 10. Velocities of cells were then pertur bed, models were updated, and temperature was reduced after each 10,000 accepted models until th e criterion of convergence was satisfied. Because of the non-uniqueness of the inversion pr oblem and to further a void unusual artifacts, in all three cases, 50,000 accepted models from the last five temperature levels and that have leastsquared errors within a few percent from the sm allest error are averaged to derive the final inverted model, and the associated uncertainty in the form of a coefficient of variation (COV). First, consider the results for model 1 pr esented in Figure 3-1. A primary feature of refraction tomography is delineati on of both vertical and horizontal changes in seismic wave velocity, e.g., detection of anomalies. Here it is observed that the inversion is able to reasonably characterize the presence, location, and general shape of the not ch. However, it is also observed that the inverted model is not identical to the tr ue model. Most notably, as discussed by Sheehan et al. (2005) and others, refraction tomography will always model a sharp contrast in velocity with a gradient in velocity (thi s is usually a result of some method of smoot hing of the velocity model such as utilized herein ). This reality largely accounts for the differences in the two models. In addition, refraction tomo graphy is not able to character ize zones of ma terial outside of the ray coverage. Here is where the uncertainty results (Figure 3-1c) are helpful. The uncertainty is low in zones with good ray coverage such as zones near th e surface, and high in 39

PAGE 40

zones with poor ray coverage, such as zones near the bottom corners of th e medium. In zones of high uncertainty, the inversion routine simply reports the velocity as the average value from all values randomly and uniformly withdrawn between the minimum and maximum velocity constraints selected by the user. Otherwise, in th e zones of low uncertaint y, the inverted velocity is independent of the constraints, and thus more reliable. From these results, it appears that the delineation between low and high uncertainty is reasonably established at a COV value of approximately 20%. Finally, the result from the presented technique seems better than that shown in Sheehan et al. (2005) for commercial codes GeoCT-II, SeisImager, and Rayfract because the results from these codes all show a high velocity artifact under the notch, and the notch is not well detected. Second, the results for model 2 are presente d in Figure 3-2. Similar to model 1, the anomaly (stair step) is reasonably recovered for both the location and general shape. The primary difference between the true and inverted models is again due to the gradient. The gradient is mostly inserted at the interface be tween layers, and the interface ca n be interpreted from areas of maximum gradient in the inverted model, such as the contour line of 2000 m/s shown in Figure 3-2b. In addition, the result from the presented te chnique seems better than that from GeoCT-II, and similar to those from SeisImager and Rayfract. Last, the results for model 3 are presented in Figure 3-3. It is observed that the inverted model is an excellent recovery of the true model. Both velociti es and interfaces between layers are very well inverted. The result from the pres ented technique is consis tent with that obtained from all three commercial codes. 3.3.2 Applications on Real Test Data The results from the synthetic models have s hown the capability of the presented technique in dealing with high lateral and vertical variat ion profiles. In the end, to gain acceptance and 40

PAGE 41

wide-spread application, it must be demonstrated that results from the technique compare well with ground truth information obtain ed from real test sites. Thus, it has been applied to a welldocumented test site. The University of North Florida (UNF) and the University of Florida (UF) have developed a Florida Department of Transportation (FDOT ) dry retention pond into a karstic limestone geophysical/ground proving test site in Alachua Count y, Florida. The site contains a number of survey lines and five PVC-cased boreholes extending to approximate ly 15 m. The test site is unique because the northern portion of the retention pond commonly experiences sinkhole activity, whereas the sout hern portion rarely expe riences sinkhole activity. The geophysical/ground proving test site is lo cated outside of Newberry, Florida on State Road 26 in Alachua County. The site is a pproximately 29 km from Gainesville, and approximately 150 km from Jacksonville. The test site is a dry reten tion pond, approximately 1.6 ha in size. The northern portion of the site ha s been susceptible to sinkhole formation and a number of large sinkholes have formed and been repaired. However, the southern portion has been relatively free of sinkholes and is an ideal location for characterizing karst limestone sites. The southern portion of the test site was subdivided into 26 nor th-south survey lines equally spaced a distance of 3.0 m apart. The lines were labeled A through Z from west to east across the site, and each line was 85.3 m long, with station 0 m located at the southern end of the site. Two 36.6-m long refraction tests were conducted end-to-end along lettere d site lines A, F, K, P, U, and Z, and beginning at stati on 0. Thus, the 12 refraction test s covered six of the site lines from station 0 to 73.2 m. Each 36.6-m long refraction test was conducte d with 4.5 Hz vertical geophones spaced equally at 0.61 m, for a total of 61 measurements. Seismic energy was created by vertically 41

PAGE 42

striking a metal ground plate with an 89 N sle dgehammer, thus producing compression wave (Pwave) first arrivals. Shot locations were spaced at 3.0-m intervals along the 36.6-m line and starting at 0 m, for a total of 13 shots. Since a 32-channel dynamic signal analyzer was used to collect time records, and each line required 61 measurements, each line was conducted in two stages. In stage one, 31 geophones were placed at 0.61-m intervals from station 0 m to station 18.3 m. In stage two, 31 geophones were placed at 0.61-m intervals fr om station 18.3 m to station 36.6 m. Since there was a designed overlap between the tw o stages at station 18.3 m, a total of 61 measurement locations were collected. For each stage, time records were collected at each of the 13 shot locations, and then the data we re combined to produce a complete shot gather for the survey line from station 0 m to station 36.6 m. Data from gridline A is chosen to present herein because these data displayed the most variable condition along the line. The measured travel time is presented in Figures 3-4 and 3-5. Cone penetration tests (CPT), geotechnical borings, and standard penetrati on tests (SPT) were also conduct ed on gridline A for partially verifying the inverted P-wave tomograms. To run the inversion, the veloc ity constraints and the medium depth need to be established. Even though the inversion techniqu e does not require strict velocity constraints for convergence, required computer time can be reduced significantly if they are available. To save computer time, an estimate of the minimum and maximum layer velocities must be determined. For these estimates, results from a simple two-layer velocity profile determined via traditional refraction analysis procedures (Burger [ 1992] and Redpath [1973]) were empl oyed. These simple velocity models were generated using trav el time curves from only the two outside source locations (e.g., 0 and 36.6 m for data set A_0-36.6) of the refrac tion surveys. From the simple profile thus determined, the first layer velocity was used to determine the minimum velocity constraint, and 42

PAGE 43

the second layer velocity was used to determine the maximum velocity constraint. For both data sets presented herein, the mini mum and maximum velocities we re determined as 300 m/s and 3000 m/s, respectively. It is very important to determine a reasona ble medium depth to get credible inversion results. If one assigns a very large medium, new models will be easily accepted because models are highly perturbed at less sensitive zones, or zones outside of ray coverage, and the inversion may not well characterize the velocity structure. For cases without any prior information, one can assign a large depth such as equa l to the geophone spread, and run a trial inversion. This run with a large medium may not give a high resolution of velocity structure, but the uncertainty result will provide rational information about the dept h of investigation. The delineation between low and high uncertainty zones can th en be used to determine the depth of the medium for an improved inversion run. Alternatively, the depth of the medium can also be estimated from experience from 1/3 to 1/2 of the geophone sprea d. For both data sets presented herein, the depth of the medium was simply taken as 1/3 of the geophone spread. In the inversion process, for convenience of calculating the firs t-arrival time in the forward modeling, the medium was discreti zed with a grid size equal to the geophone spacing of 0.61 m. The initial models were specified the same for both tests with a constant velocity of 1000 m/s, and the initial average least-squared erro rs were 63.0082 ms2 for data set A_0-36.6, and 30.1335 ms2 for data set A_36.6-73.2, thus the initial temperature was 10 for both. Velocities of cells were perturbed by randomly assigning values between the constraints, models were updated with the number of accepted models at each temperat ure level equal to 10,000, and the temperature was lowered until the criterion of convergence wa s reached at which the least-squared errors were 0.8652 ms2 for data set A_0-36.6, and 0.8927 ms2 for data set A_36.6-73. The criterion of 43

PAGE 44

convergence was not satisfied until the estimate d first-arrival times from the final accepted model were very close to the obs erved travel times (shown in Figur es 3-4 and 3-5). In order to reach the criteria of convergence, both inversion r uns were required to test more than one million trial models, and each run took a few hours on a standard laptop computer. Because of noise, subjective judgments in ma nually picking first arrival times, and nonuniqueness of inverted solutions, there is no guarantee that the model corresponding to the smallest error is closest to the true model. Th erefore, the inversion re sults should be inferred from many accepted models clustering around th e global minimum and having similar errors, instead of from only one model that has the smallest error. As with the synthetic models, the last 50,000 accepted models were used to derive the final inversion results and the associated uncertainties. Figures 3-6a and 3-7a present two-dime nsional compression (P) wave refraction tomograms determined from test data collected along the surface. These tomograms indicate that the P-wave velocity at this site generally increases with depth, and that the specific pattern of increase depends on lateral location across the s ite. Figures 3-6b and 3-7b show the uncertainties associated with the inverted P-wave profiles. Again, the uncertainty is consistent with expectations: low uncertainty in zones with good ray coverage (zones near the surface), and high uncertainty in zones with poor ray coverage (zones near the bottom corners of the medium). With the confidence gained from running synthetic models, the depth of inve stigation is about 10 m at the middle of the model where the COV from su rface to that depth is less than 20 percent. Following refraction data collection and analys is, invasive ground proving information was collected at the site to provide partial verificati on of the refraction test result interpretations, and these data are presented in the following paragraphs. 44

PAGE 45

Cone Penetration Tests (CPT): Ten CPT soundings were conducted at strategic locations across the site. Because line A displayed the most lateral variability along the line, four of the 10 soundings were conducted along line A. These four tests were located at the following horizontal stations: 19.8, 39.6, 44.2, and 65.5 m. The measured tip resistance results ar e shown in Figure 3-8, and the length of each test run is shown atop the tomogr ams in Figures 3-6a and 3-7a. These results are compared with the refraction tomograms from Figures 3-6a and 3-7a as follows: At station 19.8 m, the CPT tip resistance approa ched a large value of 30 MPa, and the test was terminated at a depth of about 9.2 m. St ation 19.8 m is near the middle of a valley feature on the tomograms (Figure 3-6a), and the CPT tip terminates at a P-wave velocity of m/s. The sounding at 39.6 m was terminated before the tip resistance appr oached a large value because the CPT rod system was bending se riously to the south as penetration was attempted. It is interesting to note on the tomograms (Figure 3-7a) that station 39.6 m is slightly to the left of a bloc k/pinnacle feature, and bending of the CPT rod to the south at the site (to the left or lowe r station number on tomogram) is consistent with this block feature. At stations 44.2 and 65.5 m the tests were te rminated at shallow depths less than 0.5 m because the CPT tip resistance approached a large value of 30 MPa. Stations 44.2 and 65.5 m are both located near the top of block/pinnacl e features on the tomograms (Figure 3-7a). However, in contrast to stat ion 19.8 m, the CPT tips terminated at P-wave velocities less than 1000 m/s at stations 44.2 and 65.5 m. It is reported that small, rock outcrops were visible near st ations 62.5-64 m and 66.1-67.1 m, which are on both sides of the CPT sounding at 65.5 m. Geotechnical Borings and Standard Penetration Tests (SPT): Eight geotechnical borings and SPT soundings were conducted at strategic locations across the site. Similar to above, the refraction tomogr ams and CPT results were used to select these locations. All of the borings in cluded drilling and recovery of rock cores through a minimum of 3 m of material, and in one co re, through 11 m of material. The following information is provided: 45

PAGE 46

Three of the eight borings were located along line A at the following stations: 19.8, 35.6, and 65.5 m. These borings coincided with CP T tests at 19.8 and 65.5, while the boring at 35.6 m was slightly to the left of the CPT sounding at 39.6 m. At station 19.8 m, the boring was advanced through predominantly sand overburden soil having SPT N-values less than 10 to a depth of 9.0 m. Below 9.0 m, coring was conducted to a depth of 12.5 m, and the material was reported to be tan limestone with fossils throughout. The recovery of this material was 100% throughout, and the rock quality designation (RQD) was reported as 100, except for the first 0.5 m which was broken at the top (RQD approximately 85). At station 35.6 m, the boring was advanced through sandy clay a nd sand overburden soil having SPT N-values between 7 and 10 to a depth of 2.3 m. Below 2.3 m, coring was conducted to a depth of 11.0 m, and the materi al was reported to be predominantly light tan to white limestone with fossils throughout. At station 65.5 m, the boring was advanced through sand overburden so il to a depth of only 0.3 m. Below 0.3 m, coring was conducted to a depth of 9.4 m, and the material was reported to be predominantly light tan to wh ite limestone with fo ssils throughout. The recovery of this material varied between 80-100%, with approximately 60% of the run reported at the 100% recovery le vel. The boring notes report that several zones appeared to be weak and broken, and the RQDs varied widely between 30 and 100. Nearly half the run indicated an RQD between 60 and 80, a s hort distance (10%) at RQD of 100, and the remaining 40% reported a RQD between 30 and 60. These results are compared with the refracti on tomograms from Figs. 6a and 7a, and the CPT results from Figure 3-8 as follows: The CPT and boring information at stati ons 19.8 and 65.5 m appear to be in good agreement. At 19.8, CPT testing was termin ated upon approaching a limiting large value at a depth of 9.2 m, while boring information reported the overburden soil/rock interface at a depth of 9.0 m. Similarly, at 65.5, CPT testing was terminated at a depth less than 0.5 m, while the top of rock was reported at 0.3 m. The undulating, valley/bowl to block/pinnacle features note d along the P-wave velocity tomograms appear to be the result of lateral variation in the overburden soil/rock interface along the length of the refracti on line. At 19.8, a valley/bowl feature appears, and the top of rock is found at 9.0 m, while at 65.5, a block/pinnacle feature appe ars, and the top of rock is found at 0.3 m. While the undulating features in the velocity tomograms ar e generally indicative of the soil/rock interface, the top of rock does not appear along a constant contour of P-wave velocity. Within a valley (station 19.8), the top of rock is found at approximately 2000 m/s, while at the top of a block (station 65.5), the top of rock is found at a velocity of approximately 500 m/s. However, the rock under station 19.8 m was reported competent and intact, with 100% recovery and RQD of 100 throughout all but a short length at top of 46

PAGE 47

core run, while the rock under st ation 65.5 m was reported to be of lower quality. Velocity differences between these two materials should be expected. Finally, it is interesting to note that the CP T appears to penetrate a particulate, sand material of higher velocity (station 19.8) easier than it will a rock of lower velocity (station 65.5). A possible explanation is as follows. Velocity is related to the small-strain modulus of the material, while CPT tip resistance is re lated to bearing capacity or strength of the material. Even though the rock at station 65.5 has a relatively low velocity as measured through a large volume of material the local strength beneath the cone tip is st ill large. A large, broken mass of this material under lo w confinement near the ground surface has low velocity. However, the local, broken pieces are still an intact, cemented material, and highly resistant to local CPT penetration. Alternatively, th e particulate, sand material under large confinement at 9 m is considerably stiffer, yet will undergo local shear failure under CPT penetration. Thus, these result s reinforce the premise that good site characterization practice should include measurement of multiple parameters to fully understand expected behavior. In summary, it would appear that the inva sive site characteri zation results provide excellent corroborating evidence for the refraction tomograms presen ted herein, and suggest that the refraction tomograms provide valuable info rmation regarding subsurface characteristics. 3.4 Conclusion for Chapter 3 A global optimization scheme based on simulate d annealing is presented to obtain nearsurface velocity profiles from travel times. Although the presented t echnique requires more computer time than the three mentioned commercial packages, it has the own superiorities. First, this inversion technique does not depend on the in itial model and becomes important in regions where prior information about subsurface profiles is not available. Second, simulated annealing provides a suite of final mode ls clustering around the global solution and having comparable least-squared error. This provides an inversion result by averaging all of these models to mitigate the influence of noise and the nonuniqueness of the inversion solutions. Last, the technique also enables to determine the uncertainties associated with inverted results. In cases, the inversion results of subsurface formations are used for the design of engineering structures such as foundations, the uncertainty will be particularly useful in implementing the new LRFD design 47

PAGE 48

methodology that can explicitly account for spatia l variability and uncertainty in design parameters. By using an extremely fast forwarding mode l solution and an acceptance rule in cascade to reduce the computer time to a few hours on a sta ndard laptop, the technique is feasible for practical engineering inversion problems. 48

PAGE 49

Velocity (m/s) Coefficient of Variation (%)a) c) 05101520253035404550 Distance (m) -15 -10 -5 0Depth (m) 800 1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 3000b) 05101520253035404550 Distance (m) -15 -10 -5 0Depth (m) 1000 m/s 3000 m/s 05101520253035404550 -15 -10 -5 0 05101520253035404550 Distance (m) -15 -10 -5 0Depth (m) 0 10 20 30 40 50 Figure 3-1. Synthetic model 1: a) true model, b) inverted model, and c) uncertainty associated with the inverted model 49

PAGE 50

Velocity (m/s) Coefficient of Variation (%)a) c) 05101520253035404550 Distance (m) -15 -10 -5 0Depth (m) 1000 1200 1400 1600 1800 2000 2200 2400 2600b) 05101520253035404550 Distance (m) -15 -10 -5 0Depth (m) 1000 m/s 3000 m/s 05101520253035404550 -15 -10 -5 0 05101520253035404550 Distance (m) -15 -10 -5 0Depth (m) 0 5 10 15 20 25 Figure 3-2. Synthetic model 2: a) true model, b) inverted model, and c) uncertainty associated with the inverted model 50

PAGE 51

0102030405060708090100Distance (m) -15 -10 -5 0Depth (m) 0102030405060708090100 -15 -10 -5 0 Velocit y ( m/s ) Coefficient of Variation ( % ) a) c) b)1000 m/s3000 m/s 0102030405060708090100Distance (m) -15 -10 -5 0Depth (m) 1000 1500 2000 2500 3000 3500 4000 0102030405060708090100Distance (m) -15 -10 -5 0Depth (m) 0 10 20 30 40 50 601000 m/s 2000 m/s 4000 m/sFigure 3-3. Synthetic model 3: a) true model, b) inverted model, and c) uncertainty associated with the inverted model 51

PAGE 52

0 5 10 15 20 25 30 35 0 5 10 15 20 25 30 35 40 Distance (m)Arrival time (ms) Observed Estimated Figure 3-4. Comparison between the estimated and obs erved first-arrival time for the real data A_0-36.6 35 40 45 50 55 60 65 70 75 0 5 10 15 20 25 30 35 Distance (m)Arrival time (ms) Observed Estimated Figure 3-5. Comparison between the estimated and obs erved first-arrival time for the real data A_36.6-73.2 52

PAGE 53

P wave velocit y ( m/s ) Coefficient of Variation (%)a) b) 0 5 101520253035 Distance (m) -10 -8 -6 -4 -2 0Depth (m) 0 2 4 6 8 10 12 14 16 18 20 22 24 26A-19.8 0 5 101520253035 Distance (m) -10 -8 -6 -4 -2 0Depth (m) 0 5 101520253035 -10 -8 -6 -4 -2 0 300 500 700 900 1100 1300 1500 1700 1900 2100 2300Figure 3-6. Inversion result for th e real data A_0-36.6: a) invert ed model, and b) uncertainty associated with the inverted model 53

PAGE 54

P wave velocit y ( m/s ) Coefficient of Variation (%) 40455055606570 Distance (m) -10 -8 -6 -4 -2 0Depth (m) 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34A-39.6 A-44.2 A-65.5a) b) 40455055606570 Distance (m) -10 -8 -6 -4 -2 0Depth (m) 400 600 800 1000 1200 1400 1600 1800 2000 2200 40455055606570 -10 -8 -6 -4 -2 0 Figure 3-7. Inversion result for th e real data A_36.6-73.2: a) invert ed model, and b) uncertainty associated with the inverted model 54

PAGE 55

Figure 3-8. Cone penetrati on tip resistance, line A 55

PAGE 56

CHAPTER 4 INVERSION OF COMBINED BOREHOLE AND SURFACE TRAVEL TIMES 4.1 Introduction Site characterization for design of deep foundati ons is very crucial, as unanticipated site conditions still represent the most common and mo st significant cause of problems and disputes that occur during construction. The problem is particularly acute in karst terrain where subsurface conditions are often hi ghly variable. The problem is further exaggerated by increased use of single, large-diameter, non-redundant drill ed shafts as foundation elements. Reliance on traditional subsurface characterization methods that invasively probe and sample a very small volume of material frequently results in problems and inefficiencies. For example, such techniques often produce poor assessment of the spat ial variation in subsurface conditions. This can lead to use of very conser vative design parameters and produ ce expensive design solutions. Perhaps even worse, poor assessment of spatial vari ation can also result in failure to identify a significant anomaly, which can lead to substr ucture failure. An additional concern with traditional techniques is that they are often unreliable in characterizing the engineering parameters of soft rock. For example, the standa rd penetration test (SPT) can fracture soft rock during driving, producing results that are too low, and design parameters that are too conservative or construction claims for differing site conditions. Characterizing the material to contain and surround a drilled shaft foundation could be conducted using a cross borehole tomography technique. A minimum of 3-4 boreholes would be required to develop an image. Wh ile this approach certainly has merit, it may be too costly in that multiple boreholes would be required. In this chapter, the coupling of so-called down-hole and refraction tomography techniqu e is presented as an alterna tive method, which uses only one borehole. Here a string of sensors is placed both horizontally along the ground surface and down 56

PAGE 57

a single borehole. The borehole sensors could be contained within another apparatus like a CPT or rock shear device. With the sensors so loca ted, energy is created and propagated from the ground surface and then detected via the sensor strings. The combined travel times are derived by picking first-arrival signals and used fo r inversion to obtain the subsurface profile. The addition of borehole data to surface data is expected to improve inversi on results, thus the need to know the capability of the combined surface and borehole data is critical. To both qualitatively and quantitatively appraise the capabili ty of the data, the gl obal inversion scheme presented in chapter 3 is used on many synthetic a nd real test data sets wi th or without boreholes to obtain both inverted profiles and associated quantitative uncertain ties. A comparison of tomograms utilizing the combined borehole and surface data against tomograms developed using just the surface data suggests that significant additional resolution of inverted profiles at depth are obtained with the addition of a borehole. Th e uncertainty estimates provide a quantitative assessment of the reliability of the interpreted pr ofiles. It is also found that the quantitative uncertainties associated with the inverted pr ofiles are significantly reduced when adding a borehole. Refraction tomography using combined data ca n be utilized for site characterization of deep foundation design. It enables to provide important physical pa rameters of material at the socket that often carries a majority of load fr om deep foundations. The in version results of the combined data enable to characterize spatial variability in geotechnical engineering physical parameters of subsurface formations useful in the design of deep foundations. This will be particularly useful in implementing the ne w LRFD design methodology that can explicitly account for spatial variability in design parameters. The refraction tomography using the combined data will improve the accuracy and re liability of engineering parameters determined 57

PAGE 58

for design of deep foundations, and will increase the volume of material actually evaluated in a subsurface investigation. The improved charac terization will lower uncertainty and risk, and reduce overall cost of subsurface construction. 4.2 Inversion of Synthetic Data 4.2.1 Synthetic Model 1 Synthetic model 1 (Figure 4-1a) was designed to illustrate the benefit of using the borehole data. It includes three layers increasing in velo city with depth and a lowvelocity zone buried in layers 2 and 3. Two sets of trav el time data were created from this model. The first set was calculated with 51 surface rece ivers at 1 m spacing and 17 shots at 3 m spacing. The second set was calculated with 51 surface receivers and 20 bor ehole receivers at 1 m spacing, and 16 shots at 3 m spacing. All the shots were on the su rface and within the surface geophone spread (no shot at borehole location). The borehole was pla ced at the middle of the medium, and it runs through the low-velocity zone to the bottom of the medium. This tr avel time data was then used to invert the velo city structures. To run the inversion, the medium sizes were kept the same as the true models, and discretized with a grid size of 1 m for ease in picking the fi rst-arrival times from forward modeling. The velocity constraints were assigned by slightly increasing the maximum true value and decreasing the minimum true value, and thus the velocity was allo wed to vary between 100 m/s and 5000 m/s. The initial model was specified with a constant velocity of 2000 m/s, which leads to an initial temperature of 10. Velocities of cells were then perturbed, models were updated, and temperature was reduced after each 10,000 accepted models until the criterion of convergence was satisfied. Because of the non-uniqueness of the inversion problem and to further avoid unusual artifacts, 50,000 accepted models from the last five temperature levels and 58

PAGE 59

that have least-squared errors within a few percen t from the smallest error are used to derive the final inverted model and th e associated uncertainty. First, consider the inverted model (Figures 4-1b) using only surface travel time data. A primary feature of refraction tom ography is delineation of both vert ical and horizontal changes in seismic wave velocity, e.g., detection of anomalies. Here it is observed that the inversion is able to reasonably characterize the presence, location of the low-velocity zone. However, it fails to recover the general shape and true velocity of 1000 m/s of the lo w-velocity zone. This is the typical limitation of refraction tomography usin g only surface data, because the technique can only model structures with velocity increasing with depth. In addition, refraction tomography is not able to characterize zones of material out side of the ray coverage. Here is where the uncertainty results (Figures 4-1c) are helpful. The uncertainty is low in zones with good ray coverage, such as zones near the surface, and hi gh in zones with poor ray coverage, such as zones near the bottom corners of the medium. In zones of high uncertainty, the inversion routine simply reports the velocity as the average value from all values randomly and uniformly withdrawn between the minimum and maximum velocity constr aints selected by the user. Otherwise, in the zones of low un certainty, the inverted velocity is independent of the constraints and thus more reliable. From these results, it appears that the delinea tion between low and high uncertainty is reasonably established at a COV value of approximately 20%. Here the inversion fails to provide rational information of th e velocity structure below the depth of investigation of 15 m. This can be expl ained as follows. The half space (third layer) in this model has the highest velocity, thus the fastes t rays that travel through the half space can be detected on the surface if the geophone spread is large enough. However, these fastest rays only travel with in a few meters at the top of the half space rega rdless how very large the geophone 59

PAGE 60

spread is used. Therefore, the technique using onl y surface data fails to characterize the velocity structure below these few meters from the top of th e half space. In this case, the borehole data is necessary. The borehole location can be decided fr om the inversion results of the surface data (Figures 4-1b and 4-1c), which provide the cr edible information of the location of the lowvelocity zone. The borehole should be located at or near the low-velocity zone in order to well characterize this zone. Second, consider the inverted model (Figure 4-1d) using combined borehole and surface data. It is observed that the presence, location, gene ral shape, and true velocity of the lowvelocity zone are well characterized. Velocity structure at depth n ear the borehole is well recovered. However, even using th e combined data, the inverted model is not the same as the true model. Most notably, as discussed by Sh eehan et al. (2005) and others, refraction tomography will always model a sharp contrast in velocity with a gradient in velocity. This reality largely accounts for the differences in th e two models. In addition, from the uncertainty result (Figure 4-1e), the zone of investigation is again taken as the area having COV less than 20%. Finally, a comparison of the tomogram (Figures 4-1d) utilizing the combined borehole and surface data against the tomogram (Figures 4-1b) developed usi ng just the surface data suggests that significant additional resolu tion of inverted profiles at de pth are obtained with the addition of a borehole. The borehole data help well charact erize the low-velocity z one and also help to increase the depth of investig ation near the borehole. Simila rly, by comparing the uncertainty results (Figure 4-1c against Figure 4-1e), the unce rtainty is significantly re duced at depth for the case of using the borehole data. 60

PAGE 61

4.2.1 Synthetic Model 2 Synthetic model 2 (figure 2a) was designed to investigate whether a low-velocity zone can be characterize if the borehole does not run throu gh this zone. It also includes three layers increasing velocity with depth and a low-velocity zone of 6 m by 6 m buried in layers 2 and 3. This model is 10 m deeper than model 1 and the lo w-velocity zone is plac ed nearby the borehole, which runs from the top to the bottom at the middle of the medium. Two sets of travel time data were created from this model using the same test layouts as model 1, except using 30 borehole receivers at 1 m spacing in this case. As described for model 1 above, the inversion was conducted for both travel time data sets and the results are as follows. For the surface da ta only, the inversion results (not shown here) can not provide any information about the low-velo city zone. This can be explained that the lowvelocity zone is buried at a significant depth compared to the surface geophone spread, and also the surface refraction tomography can only model stru ctures with velocity increasing with depth. For the combined data, the inverted model is pr esented in Figure 4-2a. It is observed that the velocity structure near the borehole from the top to the bottom of the model is reasonably characterized. Particularly, the low-velocity z one is successfully recovered even though the borehole does not run through this zone. However, the inverted model and the true model are again not the same. The primary difference between the two models is du e to the gradient. Figure 4-2c presents the uncertainty associated with the inverted model. Similar to that of model 1, the uncertainty is low in zones with good ray coverage, such as z ones near the borehole, and high in zones with poor ray coverage, such as zones near the bottom corners of the medium. It is interesting to note that the location and general shape of low-velocity zone also can be recognized in the uncertainty image as a high uncerta inty area inside a low uncertainty zone. It is understood that the fastest rays from shots to ge ophones tend to avoid the lo w-velocity zone if no 61

PAGE 62

geophones are placed in this zone. Therefore, the low-velocity zone has lower ray coverage than high-velocity zones around, and it is less defined and has higher uncertainty. The indication of low-velocity zone in the uncertainty image si gnificantly increases the credence of a real lowvelocity zone instead of an artifact of the inversion product. 4.2.1 Synthetic Model 3 Model 3 (Figure 4-3a) was designed to optimize te st layouts in order to obtain a reasonable inverted profile near the borehole with the fewest possible number of shots. This model is similar to model 2, including 3 layers and a buried low-velo city zone of 4 m by 4 m, and the borehole is placed at the middle from the top to the bottom of the medium. The focus here is to delineate the velocity structure at de pth within 5 meters around the borehole and it is expected the goal can be achieved by using only a few shots on the ground surface. To optimize the number of shots, three cases were tested using three sets of co mbined data created with the number of surface shots as 2, 4, and 6. For the first case, 2 shots we re placed 5 m each side away from the borehole. And for other cases, shots were placed at both si des of the borehole at every 3m away from the borehole. Receivers were placed bo th on the surface within the shot spread and in the borehole at 1m spacing. The inversion results of the three cases are as follows: First, the inversion results us ing 2 shots are presented in Fi gures 4-3b and 4-3c. Here the inversion enables to characterize the velocity of the three layers but it fails to provide any information of the low-velocity zone. This ma y be explained that th ere are not enough rays through the low-velocity zone if using only one s hot each side of the borehole. Second, consider the inversion results using 4 shots (Figures 4-3d to 4-3e). It is observed that the velocity structure near the borehole is reasonably characterized. With respect to the low-velocity zone, the location and the true velocity of are we ll inverted. However its general sh ape is not fully recovered. Last, 62

PAGE 63

the inversion results using 6 shots are presented in Figures 4-3f and 4-3g. The velocity structure near the borehole is well characterized; especi ally the low-velocity zone is well recovered. From the inversion results for cases of model 3, it is possible to utilize just a few shot locations on the ground surface within a few meters around the borehole, surface geophones within this few meters, and a string of borehole geophones to reliably assess the subsurface properties at depth near the borehole. This is a significant reduction in te st effort in comparison to full surface arrays such as those used in models 1 and 2. 4.3 Inversion of Real Test Data 4.3.1 Newberry Test Site The description of the Newberry test site was provided in section 3.3.2. For the case presented here, the test was conducted on line K to measure both surface and borehole data. The surface data were measured with vertical geophones equally spaced at every 1 m, for a total of 37 measurements. The borehole data were measured in a borehole installed at the station 18 m with down-hole geophones at 1 m spacing, for a tota l of 17 measurements. Seismic energy was generated by vertically striking a metal ground plate with a sledgeha mmer, thus producing compression wave (P-wave) first arrivals. Shot locations were placed at 3 m spacing along 36 m and starting from 0 m, for a total of 13 shots. First consider the surface da ta only (Figure 4-4). To r un the inversion, the velocity constraints and the medium depth need to be established. Similar to surface data sets in chapter 3, the minimum and maximum velocities were de termined as 300 m/s and 3000 m/s, respectively and the depth of the medium was simply taken as 1/2 of the geophone spread. Second, for the case of the combined surface an d borehole data (Figure 4-5), the velocity constraints are assigned the same as those for th e case of the surface data only. The depth of the medium is simply taken the same as the depth of the borehole because it is expected that the 63

PAGE 64

measured data does not provide any information of the structure below the tip of the borehole. In addition, because of the influence of borehole cas ing, the data from the shot at the borehole location is excluded from the combined data. Thus, only 12 shots are used in this case instead of 13. In the inversion process of the two cases, the medium was discretized with a grid size equal to the geophone spacing of 1 m. The initial model was selected with a constant velocity of 1000 m/s, velocities of cells were perturbed by randomly assigning values between the constraints, models were updated, and the te mperature was reduced after every 10,000 accepted models until the criterion of convergence was reached. The criterion of convergence was not satisfied until the estimated first-arrival times fr om the final accepted model were very close to the observed travel times (show n in Figures 4-4 and 4-5). Because of noise, subjective judgments in ma nually picking first arrival times, and nonuniqueness of inverted solutions, there is no guarantee that the model corresponding to the smallest error is closest to the true model. Th erefore, the inversion re sults should be inferred from many accepted models clustering around th e global minimum and having similar errors, instead of from only one model that has the smallest error. Similar to synthetic model cases, last 50,000 accepted models were used to derive the inverted profiles and the associated uncertainties. Figure 4-6a presents the two-dimensional compression (P) wave refraction tomogram determined from test data collected along the su rface. This tomogram indicates that the P-wave velocity generally increases with depth within the investigation area. Figure 6b shows the uncertainty associated with the inverted P-wave profile. Again, the uncertainty is consistent to expectation, low uncertainty in zones with good ra y coverage and high uncer tainty in zones with 64

PAGE 65

poor ray coverage. With the confidence gained from running synthetic models, the depth of investigation is about 10 m wher e the COV from surface to that depth is less than 20 percent. Figure 4-7a presents the two-dimensional P-wave refraction tomogram determined from test combined data collected both along the surf ace and down the borehole. It is observed that a slightly low-velocity zone is found near the tip of the borehole. The position of this zone is also shown as the high value area in the uncertainty image (Figur e 4-7b). Also based on the low uncertainty zone, the inverted structure velocity is credible near the borehole at depth. Comparing the velocity profiles, Figure 4-6a against Figure 47a, the profiles are similar at shallow depth. This similarity is understandabl e because the shallow structure is mainly determined by the surface data, which is the same in both inversion runs At deeper depth, the profiles are very different because only the profile from the combined data is determined from the test measurement; the other is simply repor ted by inversion as an average of the minimum and maximum constraints. Similarly, a comparison of uncertainty results (Figure 4-6b against Figure 4-7b) indicates that the uncertainty is significantly redu ced at depth for the case of using a borehole. Thus, for deep foundation design, the inve rsion results of the co mbined data provide credible information of material at the socket, and also partially detect anomalies embedded near the socket. This is very important because the material at and near the socket will carry a majority of load from foundations. 4.3.2 Ocala Test Site The test site is located near Ft. McCoy, Florida. The test was conducted to measure both surface and borehole data. The surface data we re measured with ve rtical geophones equally spaced at every 1 m from 0 m to 30 m, for a to tal of 31 measurements. The borehole data were measured in a borehole installed at station 15 m with borehole geophones equally placed at every 1 m from 4 m to 18 m depth, fo r a total of 15 measurements. Seismic energy was generated by 65

PAGE 66

vertically striking a metal ground plate with a sledgehammer, thus producing compression wave (P-wave) first arrivals. Shot locations were placed from 1 m to 3 m spacing along 30 m on the surface and starting from 0 m, for a total of 20 shots. The shot spacing was 1 m for shots close to the borehole, 2 m for shots farther, and 3 m for shots at both ends. In a fashion similar to the Newberry data, tw o inversions were run for the surface data and the combined data. For both runs, the depth of me dium was selected as 18 m, equal to the depth of the deepest borehole geophone, and the minimu m and maximum velociti es were 200 m/s and 3000 m/s, respectively. The inversions converged when the estimated first-arrival times from the final accepted model were very close to the obser ved travel times (shown in Figures 4-8 and 49). Figures 4-10a and 4-10b present the two-di mensional P-wave refraction tomogram determined from test data collected along the su rface and the associated uncertainty. Here, it is determined that the depth of investigation is about 10 m, and the velocity structure in the characterized zone increases with depth. Figur es 4-11a and 4-11b show the P-wave tomogram and the associated uncertainty de termined from combined test data collected along the surface and in the borehole. It is obser ved that the addition of a bore hole significantly improves the resolution of inverted profile, reduces the associated uncertainty, and thus increases the depth of investigation. Again, the credible information of the bedrock is only obtained from the inversion of the combined data, and it is partic ularly useful for deep foundation design. 4.4 Conclusion for Chapter 4 The global inversion scheme presented in chapte r 3 is used on many synthetic and real test data sets with or without a borehole to obtain bo th inverted profiles and associated quantitative uncertainties. A comparison of inversion results utilizing the combined borehole and surface data 66

PAGE 67

against inversion results developed using just the surface data s uggests that significant additional resolution of inverted profiles at depth are obtai ned, and uncertainties ar e significantly reduced with the addition of a borehole. Employed for site characterization of deep foundation design, refr action tomography using combined surface and borehole data provides credible information of material at the socket and partially detects anomalies near the socket. This becomes very important because the material at and near the socket often carries a majority of load from foundations. The inversion results of the combined data, including inverted profiles and as sociated uncertainties, enable to characterize spatial variability in geotechnical engineering physical parameters of subsurface formations useful in the design of deep foundations. This wi ll be particularly useful in implementing the new LRFD design methodology that can explicitly acc ount for spatial variability and uncertainty in design parameters. 67

PAGE 68

05101520253035404550 -20 -15 -10 -5 0 Velocit y ( m/s ) Coefficient of Variation (%)a) c) b) 05101520253035404550Distance (m) -20 -15 -10 -5 0Depth (m) 4000 m/s 05101520253035404550Distance (m) -20 -15 -10 -5 0Depth (m) 0 10 20 30 40 501000 m/s 1000 m/s 2000 m/s 05101520253035404550Distance (m) -20 -15 -10 -5 0Depth (m) 1000 1500 2000 2500 3000 3500(Figure 4-1 is continued in the next page) 68

PAGE 69

05101520253035404550 -20 -15 -10 -5 0 Velocity (m/s) Coefficient of Va r iation (%)e) d) 05101520253035404550Distance (m) -20 -15 -10 -5 0Depth (m) 0 10 20 30 40 50 05101520253035404550Distance (m) -20 -15 -10 -5 0Depth (m) 1000 1500 2000 2500 3000 3500 Figure 4-1. Synthetic model 1: a) true model, b) inverted model of the surface data, c) uncertainty associated with the inverted m odel in figure b, d) inverted model of the combined data, and e) uncertainty associat ed with the inverted model in figure d. 69

PAGE 70

05101520253035404550 -30 -25 -20 -15 -10 -5 0 05101520253035404550 Distance (m) -30 -25 -20 -15 -10 -5 0Depth (m) Velocity (m/s) Coefficient of Variation (%)a) c) b) 05101520253035404550 Distance (m) -30 -25 -20 -15 -10 -5 0Depth (m) 4000 m/s 05101520253035404550 Distance (m) -30 -25 -20 -15 -10 -5 0Depth (m) 0 10 20 30 40 50 601000 m/s 1000 m/s 2000 m/s 1000 1500 2000 2500 3000 3500 4000 Figure 4-2. Synthetic model 2: a) true model, b) inverted model c) uncertainty associated with the inverted model 70

PAGE 71

51 01 5 -30 -25 -20 -15 -10 -5 0 51 01 5Distance (m) -30 -25 -20 -15 -10 -5 0Depth (m) Velocity (m/s) Coefficient of Va r iation (%)a) c) b) 0 5 10 15 20Distance (m) -30 -25 -20 -15 -10 -5 0Depth (m) 4000 m/s 51 01 5Distance (m) -30 -25 -20 -15 -10 -5 0Depth (m) 0 10 20 30 40 50 60 701000 m/s 1000 m/s 2000 m/s 1000 1500 2000 2500 3000 3500e) d) 51 01 5 -30 -25 -20 -15 -10 -5 0 51 01 5Distance (m ) -30 -25 -20 -15 -10 -5 0Depth (m) 51 01 5Distance (m) -30 -25 -20 -15 -10 -5 0Depth (m) 0 10 20 30 40 50 60 70 500 1000 1500 2000 2500 3000 3500 4000Velocity (m/s) Coefficient of Variation (%)(Figure 4-3 is continued in the next page) 71

PAGE 72

0 5 10 15 20 -30 -25 -20 -15 -10 -5 0 0 5 10 15 20Distance (m) -30 -25 -20 -15 -10 -5 0Depth (m) Velocity (m/s) Coefficient of Va r iation (%) 0 5 10 15 20Distance (m) -30 -25 -20 -15 -10 -5 0Depth (m) 0 10 20 30 40 50 60 70 80 1000 1500 2000 2500 3000 3500 4000f) g)Figure 4-3. Synthetic model 3: a) true model, b and c) inverted model and uncertainty using 2 shots, d and e) inverted model and uncertain ty using 4 shots, f and g) inverted model and uncertainty using 6 shots 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 45 Distance (m)Arrival time (ms) Observed Estimated Figure 4-4. Comparison between the observed and estimated first-arrival times for Newberry surface data 72

PAGE 73

a) 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 45 Distance (m)Arrival tim e (ms) Observed Estimated b) 5 10 15 20 25 30 35 -16 -14 -12 -10 -8 -6 -4 -2 0 Depth (m)Arrival time (ms) Observed Estimated Figure 4-5. Comparison between the observed and estimated first-arrival times for Newberry combined data: a) surface data, and b) borehole data 73

PAGE 74

P_wave velocit y ( m/s ) Coefficient of Variation (%)a) b) 0 5 101520253035Distance (m) -18 -16 -14 -12 -10 -8 -6 -4 -2 0Depth (m) 0 5 10 15 20 25 30 35 40 45 50 55 0 5 101520253035Distance (m) -18 -16 -14 -12 -10 -8 -6 -4 -2 0Depth (m) 0 5 101520253035 -18 0 400 600 800 1000 1200 1400 1600Figure 4-6. Newberry surface data: a) i nverted model, and b) uncertainty 74

PAGE 75

P wave velocit y ( m/s ) Coefficient of Variation (%)a) b) 0 5 101520253035Distance (m) -16 -14 -12 -10 -8 -6 -4 -2 0Depth (m) 5 10 15 20 25 30 35 40 45 50 0 5 101520253035Distance (m) -16 -14 -12 -10 -8 -6 -4 -2 0Depth (m) 0 35 -16 0 400 600 800 1000 1200 1400 1600 1800 2000Figure 4-7. Newberry combined data: a) inverted model, and b) uncertainty 75

PAGE 76

0 5 10 15 20 25 30 -10 0 10 20 30 40 50 60 Distance (m)Arrival time (ms) Observed Estimated Figure 4-8. Comparison between the observed and es timated first-arrival time for Ocala surface data 76

PAGE 77

a) 0 5 10 15 20 25 30 -10 0 10 20 30 40 50 60 Distance (m)Arrival time (ms) Observed Estimated b) 5 10 15 20 25 30 35 -18 -16 -14 -12 -10 -8 -6 -4 Depth (m)Arrival time (ms) Observed Estimated Figure 4-9. Comparison between the observed an d estimated first-arrival time for Ocala combined data: a) surface data, and b) borehole data. 77

PAGE 78

0 5 10 15 20 25 30Distance (m) -18 -16 -14 -12 -10 -8 -6 -4 -2 0Depth (m) 0 5 10 15 20 25 30 -18 0 P wave velocity (m/s) Coefficient of Variation (%)a) b) 0 5 10 15 20 25 30Distance (m) -18 -16 -14 -12 -10 -8 -6 -4 -2 0Depth (m) 0 5 10 15 20 25 30 35 40 45 50 55 400 600 800 1000 1200 1400Figure 4-10. Ocala surface data: a) inverted model, and b) uncertainty 78

PAGE 79

0 5 10 15 20 25 30Distance (m) -18 -16 -14 -12 -10 -8 -6 -4 -2 0Depth (m) 0 30 -16 0 P wave velocity (m/s) Coefficient of Variation (%)a) b) 0 5 10 15 20 25 30Distance (m) -18 -16 -14 -12 -10 -8 -6 -4 -2 0Depth (m) 0 5 10 15 20 25 30 35 40 45 50 55 400 600 800 1000 1200 1400 1600 1800Figure 4-11. Ocala combined data: a) inverted model, and b) uncertainty 79

PAGE 80

CHAPTER 5 TWO-DIMENSIONAL INVERSION OF FU LL WAVEFORMS US ING SIMULATED ANNEALING 5.1 Introduction The technique using combined travel times presented in chapter 4 demonstrates a great capability of characterizing anomalies in subsur face profiles. However, it requires the borehole data, which is not always avai lable. The question here is how do we solve problems of 2-D reverse models with only surface data? This chapter suggests a possible solution by using an inversion of full waveforms. The full waveform ha s been used by authors for the deep (> 1 km) subsurface investigation, in both the time doma in (Shipp and Singh [2002], Sheen et al. [2006]) and the frequency domain (Pratt [1999]). In all of these applications, body waves are dominant components in wave fields used for inversi on. Unlike these methods, the technique of full waveform inversion presented in this chapter uses Rayleigh waves as dominant components in wave fields. Because, Rayleigh waves are more sensitive to shear wave velocity and less sensitive to Poisons ratio or de nsity of a medium, they are usef ul for geotechnical applications where the shear wave velocity primarily controls properties of soil particles. The presented technique is based on a finite -difference solution of 2-D elastic wave equation in the time distance domain. The strength of this approach is the ability to generate all possible wave types (body waves and surface waves) and thus to simulate and accurately model complex seismic wave fields that are then compared with observed data to infer complex subsurface properties. The technique uses full in formation of elastic wave fields to increase resolution of inversion results, es pecially dealing with reverse m odels. It also employs a global inversion technique, simulated anne aling, to invert the full wave fi elds for near-surface velocity profiles. 80

PAGE 81

Unlike local inversion techniques, simulate d annealing attempts to find the global minimum of the misfit function by searching over a large parameter space and uses more global information to update the current position, thus it likely converges to the global minimum. In addition, simulated annealing provides a suite of final models clustering around the global solution and having comparable le ast-square error. This allows getting the inversion result by averaging all of these models to mitigate the influence of noise and the non-uniqueness of the inversion solutions. The most important advantage of simulated annealing is to avoid be ing trapped in local minima, and thus to allow final inversion results to be independent of the initial model. However, this technique can require significant computer time, especially if the model contains a large number of model parameters. Th is disadvantage is reduced herein by using a fast forwarding model solution, and by assuming that interfaces w ithin a medium are multi-linear in order to limit the number of model parameters. The technique is first tested on many different synthetic data sets cr eated from challenging reverse models with high-velocity and low-velocity layers at different depths that are hardly inverted by traditional inversion methods usi ng only dispersion property of Rayleigh waves. Then it is applied to real expe rimental data sets, and the inversion results are compared to invasive test results including those of cross-hole, SPT N-value, a nd material log, or compared to results of independent refraction tests. Inversion results from both synthetic and experimental data show superiority of this technique ove r the traditional met hods via its accuracy. 81

PAGE 82

5.2 Methodology 5.2.1 Forward Modeling Virieux (1986) presented a finite-difference sc heme to solve the first-order linear partial differential equation describing el astic wave propagation on a staggered grid. For 2-D case, the following two sets of equations are solved: Equations governing particle velocity: zxt v zxt vzz zx z zx xx x 1 1 (5.1) Equations governing stress tensor: x v z v t x v z v t z v x v tz x zx x z zz z x xx 2 2 (5.2) In these equations, (vx, vz) is the particle velocity vector, ( xx, zz, xz) is the stress tensor, is the density and are Lame coefficients. Because of the need of a fast forwardi ng model solution for the global optimization method, some special boundary conditions are empl oyed in this study to allow limiting medium to the smallest possible size. Thes e boundary conditions are as follows: Free-surface condition on top: 0 0zz zx (5.3) or equivalently from equation (5.2) and (5.3): 82

PAGE 83

x v z v x v z vx z z x 2 (5.4) Symmetrical condition at le ft hand side (load line) 0 0x zxv (5.5) Absorbing condition A1 (Clayton and Engquist, 1977) at bottom 0 0 z v Vp t v z v Vs t vz z x x (5.6) and at right hand side 0 0 x v Vp t v x v Vs t vz z x x (5.7) where Vs, Vp are shear and P-wave velocities. Derivatives are discretized by using center finite-differences. A ssuming equations are verified at nodes, discretization leads to a unique staggered grid, as shown at Figure 5.1. The explicit numerical scheme (Virieux, 1986), equivalent to the systems (5.1) and (5.2) is: k ji k ji ji k ji k ji ji k ji k jiTxz Txz z t B Txx Txx x t BUU2/1, 2/1, ,,2/1 ,2/1 2/1 2/1 k ji k ji ji k ji k ji ji k ji k jiTxz Tzz z t B Txz Txz x t B V V,2/1 1,2/1 2/1,2/1 2/1, 2/1,1 2/1,2/1 2/1 2/1,2/1 2/1 2/1,2/1 83

PAGE 84

2/1 2/1,2/1 2/1 2/1,2/1 ,2/1 2/1 2/1 ,1 ,2/1 ,2/1 1 ,2/12 k ji k jiji k ji k ji ji k ji k jiU V z t L UU x t ML Txx Txx (5.8) 2/1 2/1 ,1 ,2/1 2/1 2/1,2/1 2/1 2/1,2/1 ,2/1 ,2/1 1 ,2/12 k ji k ji ji k ji k ji ji k ji k jiUU x t L V V z t ML Tzz Tzz 2/1 2/1,2/1 2/1 2/1,2/1 2/1, 2/1 2/1 1, 2/1,2/1, 1 2/1, k ji k ji ji k ji k ji ji k ji k jiV V x t MUU z t M Txz Txz In these equations, k, i, and j are the indices for time, x-axis, and z-axis discretization. t, x, and z are grid steps for time, x-axis, z-axis, respectively. Numerical velocity (U,V)=(vx,vz) at time (k+1/2) t, and numerical stress (Txx, Tzz, Txz)= ( xx, zz, xz) at time (k+1) t, are computed explicitly from velocity at time (k-1/2) t and stress at time k t. B is 1/ and M, L are Lame coefficients ( ). For free-surface condition, the geometry is c hosen so that the free surface is located exactly through the upper part of the staggered grid-p oints at z=0, that is, xx, zz, and vx are located on the free surface. To update the wave fi eld in the proximity of the free surface, both vertical velocity vz and shear stress xz of one row above the surface are required. By imaging xz as an odd function around the free surface, it is assured that xz =0 there. From equation (5.4), the vz above surface can be obtained as: 2/1 1, 2/1 1,1 1,2/1 1,2/1 2/1 2/3,2/1 2/1 2/1,2/1)2( k i k i i i k i k iUU x z ML L VV (5.9) Similarly, assuming xz and vx are odd functions around the sy mmetrical line to update the wave field in the proximity of left hand side boundary. For absorbing boundary conditions (5.6) and (5.7), backward finite-difference scheme is used in both spatial and time directions to update 84

PAGE 85

velocity components vx, vz at vertical gridlines Nx, Nx-1/2 for right hand side boundary, and at horizontal gridlines Nz, Nz-1/2 for the bottom boundary. The initial conditions are set to satisfy equilibri um at time t=0, i.e., stress and velocity are zero everywhere in the medium. Then the medi um is perturbed by changing vertical stress zz at source that is often modeled as a triangle wavelet (Figure 5-4a) or a Ricker wavelet R(t) (Figure 5-4b): 2 0 22 2 0 22exp 21)( ttf ttf tRc c (5.10) where fc is the center of frequency band and t0 is time shift. For homogeneous media, the numerical stabi lity condition for this explicit scheme: 1 122 zx tVp, (5.11) and for heterogeneous cases, Vp value in the equation (5.11) should be the maximum P-wave velocity in the media. x and z need to be selected to satisf y at least ten points per minimum expected wavelength to avoid num erical dispersion, and then t is selected from the stability condition. In inversion process, t is supposed to change due to different input Vp values, thus the wave fields are interpolated to a fixed time interval for compar ing to measured wave fields. The code is developed in Matlab in which all stresses and particle velocities are calculated in matrix forms at each time step and then advanced along time direction. The accuracy of this finite-difference method (FDM) is verified by co mparing its wave field to that generated by a finite element method (FEM) that is available in commercial software Plaxis2D. Figure 5.2 shows the results from two methods are almost id entical. By using FDM, computer time can be saved significantly. On the same standard com puter, FDM needs less than 1 second, while FEM needs a few minutes to generate the wave fiel d shown in the Figure 5.2. This fast finite 85

PAGE 86

difference solution brings possibility of bei ng applied to the global optimization method for inversion presented in next section. 5.2.2 Optimization Method Rather than discussing the analogy of simulate d annealing that has been well described by Sen and Stoffa (1995), a brief desc ription of the process used in this technique is presented herein. 1. Select an initial model P so that each component Pi to satisfy maxand then calculate the least-squared error E as: min iiiPPP N k kkPgd N E1 2)( 1, (5.12) where dk and gk are the kth observed and computed responses, respectively, and N is the number of observation points, equal to the number of receivers times the number of time steps. 2. Following the idea of Ingber (1989, 1993), randomly select one parameter (say ith parameter) and perturb it according to the Cauchy probability distribution. The updating factor yi for the ith parameter is computed from the equation: 1 1)5.0sgn(12 0 0iu i iT T T T uy (5.13) In the above equation, yi varies between -1 and +1, ui is a random number between 0 and 1, T is the current temperature, and T0 is the initial temperature. Thus, the parameter Pi is updated to from its previous value by: 1 j iPj iP ) (min max 1 i ii j i j iPPyPP (5.14) The random number ui sometimes needs to be reselected to guarantee that the value from equation (5.14) sa tisfies the constraint Following the process described, the new model is now obtained and then the least-squared error corresponding to this new model is calculated. 1 j iPmax 1 min i j i iPPP 3. The new model is unconditionally accepted if its least-squared error is smaller than that from the previous model, and conditionally accepte d if its least-squared error is bigger with the following probability: T E p exp, (5.15) where E is the positive difference between the l east-squared errors of the two models. 4. Repeat steps 2 and 3 until a desired number of accepted models is completed, then reduce the temperature to the next level according to a schedule: 86

PAGE 87

2 )1( )( kT kT (5.16) where T(k) is the temperature at the kth level. The criterion for convergence requires that the relative difference between the maximum and minimum least-squared errors from all accepte d models at a particular temperature to become very small (empirically say between 0. 5% and 1%). If the cr iterion for convergence is not satisfied, the temperature is reduced to the ne xt level. At a sufficiently low temperature, the accepted models are in the vicinity of the global minimum and the optimization converges. The selection of a correct initial temperature is very important. Selecting a very high value leads to unnecessary computations, whereas starting from a low value leads to a local solution. In this study, the initial temperature is taken as the nearest lower pow er of ten to the least-squared error; for example, if the least-s quared error is 0.35 then the initial value is 0.1. In this way, the initial temperature is usually higher than the criti cal temperature and any bias in the initial model will be destroyed. The cooling schedule is a nother important parameter in the global optimization using simulated annealing. Many schedules have been tested for inversion problems with the number of model parameters up to 40, and it is found that the number of accepted models at each temperature level should be at least 200 times of the number of model parameters, and the temperature needs to reach to a very low value to obtain the global solution. In this study, the final temperature is about 10-4 to 10-5 when the initial temperature is 0.1. When the number of model parameters increases, the procedure ne eds a slower schedule and hence the time of computation increases in order to get the global soluti on. Therefore, the cooling schedule should be selected with care, particul arly when dealing with a larg e number of model parameters. 87

PAGE 88

5.3 Applications 5.3.1 Applications on Synthetic Data Synthetic models refer to earth models w hose velocity profile is assumed or known a priori. Using a known velocity structure, wave field data is calculated for an assumed test layout. This wave field data is then input to the inversion program as if it were data collected from a field test (observed data), and velocity structure (inverted model) is calculated from the observed data. By comparison of the inverted mode l against the true model, the capability of the inversion technique is appraised. 5.3.1.1 One-dimensional synthetic models Two 1-D synthetic challenging models (models 1 and 2) that include highand lowvelocity layers are presented in this section. They are from Tokimatsu et al. (1992) and were previously used for studies of inversion using surface wave velocity dispersion by ONeill (2003). Synthetic model 1 Model 1 consists of a buried low-velocity layer (second laye r) between two high-velocity layers, followed by a half space (Figure 5-3a). The observed wave fi eld was generated in Plaxis2D using 10 receivers at 4 m spacing on th e surface, and an active source placed 4 m away from the first receiver. The active source was mode led as a triangle impulsive load (Figure 5-4a), and such a source produced a wave field with dominant components between 10 Hz and 30 Hz. The observed wave field (at 10 points from A to J) is presented in Figure 5-5. It shows that the magnitude of the wave field is decayed significantly from point A (closest to the source), to point J (farthest from the source) due to damping. In or der to increase the contribution of the far field signals, the receivers were treated equally by no rmalizing the maximum ma gnitude of measured 88

PAGE 89

signals at each receiver to unity before being used for inversion, thus removing all kinds of damping (material and radiat ion) in wave propagation. Before running inversion, constrai nts, medium sizes, and the grid size need to be specified. For velocity constraints, S-wave velocities were allowed to va ry between 50 m/s and 300 m/s for the 3 layers, and between 200 m/s and 500 m/s fo r the half space. For thickness constraints, thicknesses were allowed to vary between 1 m and 5 m for layers 1 and 2, and between 5 m and 10 m for layer 3. For the medium sizes, the medium was selected as 20 m depth and 60 m width. The width was extended 20 m horiz ontally from the last receiver to limit the boundary affect during forward modeling. Finally, based on the r eceiver spacing and the minimum true S-wave velocity, the grid size was specified as 0.5 m, which can accurately model the maximum frequency of about 30 Hz at whic h the criterion of 10 mesh point s per wavelength is satisfied. To run inversion, Poissons ratio and the medium density were kept the same as those used to generate the observed data to limit the number of model parameters to 7 (3 thicknesses and 4 velocities). The inversio n began with the initial model of a constant S-wave velocity of 200 m/s and all thicknesses of 5 m, model parameters were perturbed within the c onstraints, models were updated, and the temperature was reduced afte r every 1400 (200 times of the number of model parameters) accepted models until the criterion of convergence was satisfied. The convergence was only found when the estimated data was very close to the observed (Figure 5-6), and the least-squared error was reduced from the initial value of 0.2087 to the final value of 0.0240. Figure 5-6 presents a comparison between th e observed and the estimated data. A good match is achieved at every channel along the w hole range of the offset. Figure 5-7 presents velocity histograms of all accepted models at all temperature levels. At high temperatures, the accepted models are located over a large mode l space, and thus accepted velocities are any 89

PAGE 90

random values between the constraints. Otherwis e, at low temperatures, the accepted models cluster around the global minimum, and thus the mode values of the histograms are close to those of the true model. Similarly, Figure 58 presents thickness histograms of all accepted models, and again the mode valu es are observed to be close to those of the true model. Figure 5-9 compares the inverted S-wave profile against the true model. It is found that the true model is well recovered. The differences of velocities and thicknesses between the true and the inverted models are less than 15%. These differences may be due to the incomparability of numerical methods used to gene rate the observed data (FEM) a nd the estimated data (FDM). ONeill (2003) did inversion studies on this model using dispersion data and a local inversion scheme. Even though the same syntheti c modeling methods were used for the observed and estimated dispersion data, he found that the true model wa s hardly recovered. It was relatively recovered if the true depth of the half space was fixe d during the inversion, and the initial model was close enough to the true model, such as the initial model having the true depth interfaces and a reasonable shear wave velocity. It is clear that the technique presented here requires much less prior information to invert the true model. The superiority of the technique is attributed to the full waveform data used for inversion that is more sensitive to a low-veloci ty layer (LVL) than the dispersion data, and the advantage of the global inversion scheme th at is independent of the initial model. Synthetic model 2 Model 2 consists of a buried high-velocity layer (second laye r) between two low-velocity layers, followed by a half space (Figure 5-3b). Similar to model 1, the observed wave field (Figure 5-10) was generated in Pl axis2D using 10 receivers at 4 m spacing on the surface, and an impulsive source (Figure 5-4a) placed 4 m away from the first receiver. 90

PAGE 91

In the same manner of model 1, an invers ion was conducted using the observed data of model 2. During the inversion, th e least-squared error decreased from 0.2025 for an initial model of a constant velocity of 200 m/s to 0.0378 for the final accepted model. Again, the final inverted model (Figure 5-14) is a good recove ry of the true model, correct ly inverting both interfaces and velocities. The good recovery of the velocity model leads to a good match between the observed and estimated wave fields (Figure 5-11). The velocity and thickness histograms of all accepted models are also shown in Figures 5-12 and 5-13, whose mode values are close to those of the true model. Similar to model 1, Oneil (2003 ) found that this true model wa s also hardly recovered. It was relatively recovered if the tr ue depth of the half space was fixed during the inversion, and the initial model had the true depth interf aces and a reasonable shear wave velocity. 5.3.1.2 Two-dimensional synthetic models Three 2-D synthetic models (models 3, 4, and 5) are presented here. The first two models only consist of linear interfaces, and they can be inverted by da ta from a few receivers on the surface. The two synthetic data sets from thes e simple configurations were generated by the commercial software Plaxis2D. The last model c onsists of multi-linear interfaces that require data from many receivers along the surface to rec over. For convenience, the synthetic data set from this more complicated configuration wa s generated by the finite-difference solution. Synthetic model 3 Model 3 consists of a buried high-velocity layer (second laye r) between two low-velocity layers, followed by a half space (Figure 5-15a). The medium is extended 20 m from the last receiver to limit the boundary affect. The observed wave field was generated in Plaxis2D using 6 receivers at 10 m spacing on the surface, and an active source placed 10 m away from the first receiver. The active source was 91

PAGE 92

modeled as a triangle impulsive lo ad (Figure 5-4a), and such a s ource produced a wave field with dominant components between 10 Hz and 30 Hz. To run inversion, the medium sizes were kept the same as those used in Plaxis2D. Poissons ratio and density were al so kept the same as values used to generate the observed wave field, thus the number of model parameters is 10 (4 shear veloci ties, 3 left thicknesses, 3 right thicknesses). The velocities of all layers were allowed to vary betw een 200 m/s and 1000 m/s, and the thickness was allowed to vary between 2 m and 10 m. The medium was discretized with a grid size of 1 m, and the material property at each node was assigned based on its tributary area, i.e., the material property was interpolated when interfaces lied between grid lines. Based on the grid size and the minimum true velocity the maximum frequency that can be model accurately in forward modeling is 30 Hz, thus both observed and pr edicted wave fields were low band-pass filtered to remove all signa ls above 30 Hz before comparing. The inversion began with an initial model selected with a constant sh ear velocity of 500 m/s, and left and right thicknesses of 7.5 m. With such the initial model, the initial leastsquared error was 0.1480, thus the initial temperature was taken as 0.1. During inversion, velocities and thicknesses were perturbed, models were updated, and the temperature was reduced after every 2000 accepted models until the criterion of convergen ce was satisfied at which the final smallest error was 0.0048. The criterion of convergence was satisfied when the difference between the observed and predicted wave fields became very small, as shown in Figure 5-16 where an excellent match is found. Figure 5-17 presents histograms of velocities from all accepted models for all temperature levels. It demonstrates that at high temperatur es, the accepted models are located over a large model space and at low temperatures, the acc epted models cluster around the true model. 92

PAGE 93

Therefore, the inversion result should be inferred from a suite of models that have similar leastsquared error near the global minimum. The la st 2000 accepted models at the final temperature are used to get an average model for the invers ion result (Figure 5-15b), and it is observed that the true model is excellently recovered. Synthetic model 4 Synthetic model 4 consists of a buried low-ve locity layer (second la yer) between two highvelocity layers followed by a half space as shown in Figure 5-18a. Synthetic data was generated in Plaxis2D exactly the same as model 3. Similar to model 3, during inve rsion, the least-squared erro r decreased from 0.1088 for an initial model of constant velocity of 500 m/s to 0.0066 for the final accepted model. Again, the final inversion result (Figure 518b) is an excellent recovery of the true model, correctly inverting both interfaces and velocities. The excellent recovery of the velocity model leads to an excellent match between the observed and estimate d wave fields (Figure 5-19). The histogram of all accepted models is also shown in Figure 5-20, whose mode values are close to the velocities of the true model. Synthetic model 5 Model 5 consists of three lowand high-veloc ity layers followed by a half space, with six linear segments in each interf ace (Figure 5-21a). The observed wave field was generated by the finite-difference solution using 51 receivers at 1 m spacing on th e surface, and an active source placed 10 m away from the first receiver. The ac tive source was modeled as a Ricker wavelet (Figure 5-4b) having a central fre quency of 15 Hz. The Ricker wavele t is typically better than the triangle wavelet for modeling an impulsive load, because a wave field generated by the Ricker wavelet has less numerical noise, and thus less filtering is required. 93

PAGE 94

To run inversion, the source, Poissons ratio, density, and the number of linear segments in each interface were kept the same as values used to generate the observed wave field, thus the number of model parameters was 25 (7 thickn esses in each layer of the 3 layers, and 4 velocities). The velocity was allowed to vary between 200 m/ s and 900 m/s, and the thickness was allowed to vary between 2 m and 10 m. The inversion began with the initial model of a constant velocity of 400 m/s, parameters were perturbed, and the temper ature was reduced after 5000 accepted models until convergence. The run tested more than 150,000 trial models in about half a day on a standard laptop. Figure 5-21b presents the inverted model, wh ich is taken as the average model from the last 5000 accepted models at the fi nal temperature. One more time, excellent recovery of the true model is found. Both layer veloci ties and interfaces are accurately inverted with differences less than 10% from those of the true model. The observed wave field, the estimated wave fi eld associated with th e last accepted model, and residuals are shown alongside in Figure 5-22. It is observed that the data has been well fit across the whole range of offsets, and the two wave fields are al most identical. The histogram of all accepted models is also shown in Figure 5-23. Again, the mode values are close to those of the true model. The inversion result of model 5 demonstrates that it is possibl e using a full wave field from only one shot to characterize 2-D profiles with in terfaces of linear segments of a few meters. However, one should not expect too much from ju st one shot, for example, dividing the medium into too many layers with interf aces of multiple small segments, or dividing the medium into many small cells. In such cases, full wave fields from multiple shots are required for characterization. 94

PAGE 95

5.3.2 Applications on Real Test Data The results from the synthetic models have s hown the capability of the presented technique in dealing with reverse and high variation profiles. In the end, to gain acceptance and widespread application, the technique must be dem onstrated that tomography results compare well with ground truth information obtained from real test sites. Here it is applied to well-documented testing sites. 5.3.2.1 TAMU test site Data were collected at the National Geotechni cal Experiment site (NGES) on the campus of Texas A & M University (TAMU). The TAMU s ite is well documented, and consists of an upper layer of approximately 10 m of medium dens e, fine, silty sand followed by hard clay. The multi-channel test were conducted using 31 receiv ers at a spacing of 1.22 m, giving a total receiver spread of 36.6 m, and an active source 6 m away from the first receiver. In addition, the cross-hole test was conducted at a set of nearby cased boreholes spaced approximately 3 m apart. Before running inversion, the velocity cons traints and the medium depth need to be established. Although the techni que does not require strict co nstraints for convergence, the computer time can be reduced significantly if th ey are available. To save computer time, an estimate of velocity constraints and medium si zes must be established. For the estimate of velocity constraints, a spectral analysis was em ployed. By applying the spectral analysis to the measured data set, the measured wave field was decomposed into components of different frequencies, and velocity of each component was determined. Figure 5-24 presents the normalized spectral obtained by us ing the cylindrical beamformer technique. It is observed that most energy of the measured wave field con centrates in a narrow band, and Rayleigh wave velocities are determined to va ry from 120 to 500 m/s. Thus, the constraint of shear wave velocity, which is a little bit bi gger than Rayleigh wave velocity was taken as a range of 100 to 95

PAGE 96

600 m/s. For the estimate of medium sizes, the depth of investigati on was assumed as one half of the geophone spread. Based on the geophone spacing and the dispersion property from the spectral an alysis, the grid size was selected as 0.61 m that is sustainable to the maximum frequency of 30 Hz at which the criterion of 10 mesh poi nts per wavelength is satisfied. Ther efore, the measured data was low band-pass filtered through a bandw idth of frequencies below 30 Hz to remove high frequency signals and background noise before using for invers ion. In addition, the filtered data from every receiver was treated equally by normalizing the maximum magnitude to un ity, thus removing all damping in wave propagation. The central frequenc y of the filtered data was about 15 Hz, and it was selected to be the central frequency of a Ricker wavelet for an active source in forward modeling. Because most of the signals in the full wave field measured on the surface are Rayleigh waves, which are not very sensitive to Poisson s ratio and density, and because the number of model parameters needs to be limited to redu ce computer time, the Poissons ratio and the density were kept constant as 0.3 and 1800 kg/m3, respectively for the entire medium during inversion. It was assumed that the medium consists of 3 layers and a half space, and each interface was divided into five horizontally equa l linear segments. Thus, six thicknesses for each layer and 4 velocities were searched; th e number of model parameters was 22. To run inversion, the shear velocity was allowed to vary between the constraints, the thickness was allowed to vary between 2 m and 8 m, and the temperature was reduced after every 4400 accepted models until the criterion of convergence was satisfied. The least-squared error decreased from 0.1560 for the initial model of constant velocity of 400 m/s, to 0.0160 for the final accepted model after searching over 120,000 trial models in about a half day on a 96

PAGE 97

standard laptop. The criterion of convergence was only satisfied when the estimated full wave field from the final accepted model was very close to the observed wave field as shown in Figure 5-25, where a good match is found. Because of noise and non-uniqueness of invers ion solutions, there is no guarantee that the model corresponding to the smallest error is closest to the true m odel. Therefore, the inversion results should be inferred from many accepted models clustering around the global minimum and having similar errors, instead of from only one m odel that has the smallest error. Similar to synthetic data sets, the last 4400 accepted models at the final temperature were averaged for the inversion result. Figure 5-26 presents the velocity histograms from all accepted models. It is observed that the first three layers are charac terized better than the half spa ce below. Figure 5-27 presents the 2-dimensional shear wave velocity within the ge ophone spread (0 36.6 m) from the full wave field inversion. The tomogram indicat es that the shear wave velocity structure at the test site is a lightly reverse profile with a burie d low-velocity layer (layer 3). For comparison, the Vs profile from the inversio n result at distance 10 m that is near to the boreholes used for the cross-hole test is plotted together with cross-hole measurements, SPT Nvalues, and material logs in Fi gure 5-28. It is observed that th e full wave field inverted shear wave profile compares well with the cross-hole results, and the presented technique successfully detects the thin low-velocity layer that is ha rdly detected by traditional surface wave methods that use only dispersion property of Rayleigh waves. There also appears to be reasonable consistency between the shear wave velocity re sults and the SPT N-valu es and material log. 5.3.2.2 Newberry test site The description of the Newberry test site is provided in sectio n 3.3.2. For the case presented here, the full wave field measurement was conducted at line F from station 36.6 m to 97

PAGE 98

73.2 m, using 31 vertical geophones at a spacing of 1.22 m, giving a to tal geophone spread of 36.6 m, and the wave field was generated by an active source 6 m away from the first geophone. To run the inversion of the measured full wa ve field, the velocity constraints and the medium depth need to be established. For the ve locity constraints, the cylindrical beam-former technique was again employed for spectral analysis of the measured data, and the minimum and maximum shear velocity constraints were determ ined as 100 m/s and 1200 m/s, respectively. For the medium depth, it was simply ta ken as 1/3 of the geophone spread. The medium was discretized with a grid si ze of 0.61 m, which can accurately model the maximum frequency of 40 Hz, and thus the measur ed data was low band-pass filtered to remove all components above 40 Hz before using for inve rsion. The center frequency at 18 Hz of the filtered data was selected to be the central frequency of Ricker wavelet for the active source in forward modeling. Similar to TAMU data, during inversion, the Poissons ratio and the density were kept constant as 0.3 and 1800 kg/m3, respectively, for the entire me dium. It was assumed that the medium consists of 3 layers and a half sp ace, and each interface was divided into five horizontally equal linear segments Six thicknesses for each of three layers and 4 velocities were searched, thus the total number of model parameters was 22. The inversion began with the in itial model of a constant velo city of 500 m/s, which led to the initial least-squared error of 0.225. The model was then perturbed and updated by allowing the shear velocity to vary between the constraint s, the thickness to vary between 1 m and 6 m, and the temperature was reduced after ever y 4400 accepted models until the criterion of convergence was satisfied. The convergence was f ound at the final least-squared error of 0.0360 after searching over 120,000 trial models in about a half day on a standard laptop. 98

PAGE 99

The observed wave field, estimat ed wave field associated with the last accepted model, and residuals are set alongside in Figur e 5-29. It is observed that the data has been fit relatively well across the whole range of offsets except the last few channels at both ends of the geophone array. The misfit at these channels may be explained th at the laterally variab le subsurface velocity structure can not be completely modeled by a 2-D multi-linear interface profile. To find a better fit wave field, more complexity of velocity profile needs to be added into the medium, such as a medium consisting of small cells with different velocities. Figure 5-30 presents velocity histograms of all accepted models during the inversion. Again, it demonstrates that at high temperatur es, the accepted models are located over a large model space and thus the velocities randomly vary between the c onstraints. Otherwise, at low temperatures, the accepted models cluster ar ound the global minimum model, whose layer velocities are the mode values in these histograms. Figure 5-31a presents the 2-dimensional shea r wave (S-wave) velocity with in the geophone spread (36.6-m length) from the full wave field inversion, which is the average model from last 4400 accepted models. The tomogram indi cates that the shear wa ve velocity structure at the test site is increased with depth. It is noted that this increased-velocity-with-depth profile can be well characterized by the seismic refraction tomography. For comparison, the result of the seismic refraction tomography (SRT) developed by Quigley (2006) using commercial software Rayf ract is also presented on Figure 5-31b. The refraction test data was measured on the same lo cation of the full wave field test by vertical geophones spaced equally at 0.61 m, for a total of 61 measurements, and shot locations were spaced at 3.0-m intervals along the 36.6 -m line, for a total of 13 shots. 99

PAGE 100

The comparison of the S-wave velocity profile (Figure 5-31a) generated from the full wave field against the P-wave velocity profile (Fig ure 5-31b) generated from the refraction travel times is as following: As expected, the P-wave veloci ty is about twice of the S-wa ve velocity generally over the medium. There is a good agreement in depth of the top bedrock, which is about 9 m depth in both tomograms. The interface between layer 2 a nd layer 3 in Figure 5-31a is similar to the interface (contour of 1500 m/s) in Figure 5-31b. Slopes in interfaces at the left of th e medium are found in both tomograms. There are some differences in velocity at shallow depth because the full wave field inversion assumes the medium consisting of laye rs, not consisting of cells, and the velocity of each layer is an average value of many cells. The relatively good agreement between the inve rsion results of the full wave field and travel times suggests that the full waveform tec hnique should be even ap plied to normal profiles that can be well characterized by the SRT for two reasons. First, the full waveform technique can characterize sharp contrasts in velocity, which SRT always models by gradient in velocity. Second, the full waveform technique requires much less effort in both field testing and manually data processing than that of SR T. The SRT usually requires many shots in data measurement, and days for data processing, e.g., ma nually picking the first-arrival signals. 5.4 Conclusion for Chapter 5 A global inversion technique ba sed on simulated annealing is presented to obtain nearsurface velocity profiles from full wave fields. This inversion technique does not depend on the initial model and becomes important in regions where the prior information about subsurface profiles is not available. Simulated annealing pr ovides a suite of fina l models clustering around the global solution and having comparable least-squa re error. This allows getting the inversion 100

PAGE 101

result by averaging all of these models to mitigate the influence of noise and the non-uniqueness of the inversion solutions. The capability of the technique to accuratel y invert 2-D reverse models demonstrates superiority over 1-D traditional methods that us e only dispersion property of Rayleigh waves. In addition, the technique though re quires significant computer tim e for inversion; it possibly replaces the seismic refraction tomography to characterize normal profiles (increasing velocity with depth) because it requires much less effort in both site testing and manually data processing. 101

PAGE 102

Figure 5-1. Discretization of the medium on a staggered grid. a) b) X, mY, m 0 20 40 60 80 0 5 10 15 20 25 30 35 300 m/s 200 m/s 400 m/s 0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 Receiver @ 30mNormalized magnitude FDM FEM 0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 Receiver @ 60m time in secNormalized magnitude FDM FEM Figure 5-2. Comparison between wa ve fields generated by FDM and FEM: a) input model shear velocity Vs, and b) wave fields at 2 points x=30m and x=60m on the surface 102

PAGE 103

a)-20.0 -18.0 -16.0 -14.0 -12.0 -10.0 -8.0 -6.0 -4.0 -2.0 0.0 0100200300400 Shear wave velocity (m/s)Depth (m) b)-20.0 -18.0 -16.0 -14.0 -12.0 -10.0 -8.0 -6.0 -4.0 -2.0 0.0 0100200300400 Shear wave velocity (m/s)Depth (m) Figure 5-3. Shear wave velocity of 1-D models of Tokimatsu et. al. (1992): a) model 1 and b) model 2 0 0.2 0.4 0.6 0.8 1 -50 0 50 100 Time (s)Force (N) 0 0.2 0.4 0.6 0.8 1 -50 0 50 100 Time (s)Force (N) Figure 5-4. Impulsive load: a) triangle wavelet and b) Ricker wavelet 103

PAGE 104

0 0.1 0.2 0.3 0.4 0.5 0.6 -0.02 -0.01 0 0.01 0.02 Dynamic time [s] Vy [m/s] Chart 1 Point A Point B Point C Point D Point E Point F Point G Point H Point I Point J Figure 5-5. Synthetic model 1: par ticle velocity data from Plaxis2D observed at 10 points (A to F) spaced equally at 5 m on the surface 104

PAGE 105

0 0.2 0.4 0.6 0.8 -1 0 1 receiver 1 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 2 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 3 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 4 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 5Normalized magnitude 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 6Normalized magnitude 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 7 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 7 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 9 time (s) 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 10 time (s) Estimated Observed Figure 5-6. Synthetic model 1: comparison be tween estimated and observed wave fields. 105

PAGE 106

50 100 150 200 250 300 0 5000 10000 Layer#1 50 100 150 200 250 300 0 5000 10000 Layer#2Number of accepted model 50 100 150 200 250 300 0 2000 4000 Layer#3 200 250 300 350 400 450 500 0 500 1000 1500 Layer#4 Velocity (m/s) Figure 5-7. Synthetic model 1: velocity histograms from all accepted models 106

PAGE 107

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 2000 4000 6000 8000 Layer#1 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 2000 4000 6000 Layer#2Number of accepted model 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5 10 0 500 1000 1500 2000 Layer#3 Figure 5-8. Synthetic model 1: thickn ess histograms from all accepted models 107

PAGE 108

-20.0 -18.0 -16.0 -14.0 -12.0 -10.0 -8.0 -6.0 -4.0 -2.0 0.0 050100150200250300350400 Shear wave velocity (m/s)Depth (m) True model Inverted model Figure 5-9. Synthetic model 1: comparison between the true and inverted models 108

PAGE 109

0 0.2 0.4 0.6 0.8 -0.06 -0.03 0 0.03 0.06 0.09 Dynamic time [s] Vy [m/s] Chart 1 Point A Point B Point C Point D Point E Point F Point G Point H Point I Point J Figure 5-10. Synthetic model 2: par ticle velocity data from Plaxis2D observed at 10 points (A to F) spaced equally at 5 m on the surface 109

PAGE 110

0 0.2 0.4 0.6 0.8 -1 0 1 receiver 1 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 2 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 3 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 4 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 5 time (s)Normalized magnitude 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 6 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 7 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 8 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 9 0 0.2 0.4 0.6 0.8 -1 0 1 receiver 10 time (s) Estimated Observed Figure 5-11. Synthetic model 2: comparison be tween estimated and observed wave fields. 110

PAGE 111

50 100 150 200 250 300 0 5000 10000 Layer#1 50 100 150 200 250 300 0 2000 4000 Layer#2Number of accepted model 50 100 150 200 250 300 0 5000 10000 Layer#3 200 250 300 350 400 450 500 0 500 1000 Layer#4 Velocity (m/s) Figure 5-12. Synthetic model 2: velocity histograms from all accepted models 111

PAGE 112

1 1.5 2 2.5 3 3.5 4 4.5 5 0 1000 2000 Layer#1 1 1.5 2 2.5 3 3.5 4 4.5 5 0 2000 4000 Layer#2Number of accepted model 5 6 7 8 9 10 0 500 1000 1500 Layer#3 Velocity (m/s) Figure 5-13. Synthetic model 2: thickne ss histograms from all accepted models 112

PAGE 113

-20.0 -18.0 -16.0 -14.0 -12.0 -10.0 -8.0 -6.0 -4.0 -2.0 0.0 050100150200250300350400450 Shear wave velocity (m/s)Depth (m) True model Inverted model Figure 5-14. Synthetic model 2: comparison between the true and inverted models 113

PAGE 114

a) b) 01020304050607080 Distance (m) -30 -25 -20 -15 -10 -5 0Depth (m) 400 500 600 700Shear wave velocity (m/s)Figure 5-15. Synthetic model 3: a) tr ue model, and b) inverted model 114

PAGE 115

0 0.1 0.2 0.3 0.4 0.5 -1 -0.5 0 0.5 1 receiver 1 Estimated Observed 0 0.1 0.2 0.3 0.4 0.5 -1 -0.5 0 0.5 1 receiver 2 0 0.1 0.2 0.3 0.4 0.5 -1 -0.5 0 0.5 1 receiver 3Normalized magnitude 0 0.1 0.2 0.3 0.4 0.5 -1 -0.5 0 0.5 1 receiver 4 0 0.1 0.2 0.3 0.4 0.5 -1 -0.5 0 0.5 1 receiver 5 time (s) 0 0.1 0.2 0.3 0.4 0.5 -1 -0.5 0 0.5 1 receiver 6 time (s) Figure 5-16. Synthetic model 3: comparison be tween estimated and observed wave fields 115

PAGE 116

200 400 600 800 1000 0 1 2 3 4 5 x 104 Number of accepted modelsVelocity (m/s) Layer#1 200 400 600 800 1000 0 2000 4000 6000 8000 Number of accepted modelsVelocity (m/s) Layer#2 200 400 600 800 1000 0 2000 4000 6000 8000 10000 Number of accepted modelsVelocity (m/s) Layer#3 200 400 600 800 1000 0 1000 2000 3000 4000 5000 Number of accepted modelsVelocity (m/s) Layer#4Figure 5-17. Synthetic model 3: velocity histograms from all accepted models 116

PAGE 117

a) b) 01020304050607080 Distance (m) -30 -25 -20 -15 -10 -5 0Depth (m) 400 500 600 700 800Shear wave velocit y ( m/s ) Figure 5-18. Synthetic model 4: a) tr ue model, and b) inverted model 117

PAGE 118

0 0.1 0.2 0.3 0.4 0.5 -1 0 1 receiver 1 Estimated Observed 0 0.1 0.2 0.3 0.4 0.5 -1 0 1 receiver 2 0 0.1 0.2 0.3 0.4 0.5 -1 0 1 receiver 3Normalized magnitude 0 0.1 0.2 0.3 0.4 0.5 -1 0 1 receiver 4 0 0.1 0.2 0.3 0.4 0.5 -1 -0.5 0 0.5 1 receiver 5 time (s) 0 0.1 0.2 0.3 0.4 0.5 -1 -0.5 0 0.5 1 receiver 6 time (s) Figure 5-19. Synthetic model 4: comparison be tween estimated and observed wave fields 118

PAGE 119

200 400 600 800 1000 0 1000 2000 3000 4000 Number of accepted modelsVelocity (m/s) Layer#1 200 400 600 800 1000 0 2000 4000 6000 8000 Number of accepted modelsVelocity (m/s) Layer#2 200 400 600 800 1000 0 1000 2000 3000 4000 Number of accepted modelsVelocity (m/s) Layer#3 200 400 600 800 1000 0 500 1000 1500 Number of accepted modelsVelocity (m/s) Layer#4 Figure 5-20. Synthetic model 4: velocity histograms from all accepted models 119

PAGE 120

051015202530354045505560 Distance (m) -20 -15 -10 -5 0Depth (m) 300 400 500 600 700 800Shear wave velocit y ( m/s ) 051015202530354045505560 Distance (m) -20 -15 -10 -5 0Depth (m) 300 400 500 600 700 800Shear wave velocit y ( m/s ) a) b) Figure 5-21. Synthetic model 5: a) tr ue model, and b) inverted model 120

PAGE 121

0 20 40 60 0 0.1 0.2 0.3 0.4 0.5 Observed datatime (s)Receiver position (m) 0 20 40 60 0 0.1 0.2 0.3 0.4 0.5 Estimated datatime (s)Receiver position (m) 0 20 40 60 0 0.1 0.2 0.3 0.4 0.5 Residualtime (s)Receiver position (m) Figure 5-22. Synthetic model 5: comparison be tween observed and estimated wave fields 200 300 400 500 600 700 800 900 0 5 x 104 Layer#1 200 300 400 500 600 700 800 900 0 1 2 x 104 Layer#2Number of accepted model 200 300 400 500 600 700 800 900 0 2 4 x 104 Layer#3 200 300 400 500 600 700 800 900 0 1 2 x 104 Layer#4 Velocity (m/s) Figure 5-23. Synthetic model 5: velocity histograms from all accepted models 121

PAGE 122

1 0 20 40 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Observed datatime (s)Receiver position (m) 0 20 40 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Estimated datatime (s)Receiver position (m) 0 20 40 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Residualtime (s)Receiver position (m) Figure 5-24. Normalized spectral of the measured data at TAMU Figure 5-25. Comparison between observed wave fiel d and estimated wave field for the real data set at TAMU. 510152025303540 Frequency (Hz) 0 100 200 300 400Rayleig Wave Velocity (m/s500 600) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 122

PAGE 123

100 200 300 400 500 600 0 2 4 x 104 Layer#1 100 200 300 400 500 600 0 1 2 x 104 Layer#2 100 200 300 400 500 600 0 5000 10000 Layer#3Number of accepted models 100 200 300 400 500 600 0 2000 4000 Layer#4 Velocity (m/s) Figure 5-26. Velocity histograms from all accepted models for the real data set at TAMU. 123

PAGE 124

124 Figure 5-27. Inversion result of sh ear wave velocity (m/s) profile for the real data set at TAMU. SV-Wave Velocity (m/s) 0 2 4 6 8 10 12 14 16 18 20 22 24 0100200300400500Depth (m) Crosshole test Full waveform Standard Penetration Resistance (SPT N-Value)12 15 18 15 14 13 17 16 20 60 70 62 520102030405060708090 Soil Profile0 2 4 6 8 10 12 14 16 18 20 22 24 Dark Gray Clay EOB 15.2mTan Sand Tan Sandy Clay Gravel Dark Gray Clay w/Gravel Tan Silty Fine Sand Tan Silty Fine Sand Figure 5-28. Composite soil profile of TAMU.

PAGE 125

0 20 40 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Observed datatime (s)Receiver position (m) 0 20 40 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Estimated datatime (s)Receiver position (m) 0 20 40 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Residualtime (s)Receiver position (m) Figure 5-29. Comparison between observed wave fiel d and estimated wave field for the real data set at Newberry. 125

PAGE 126

0 200 400 600 800 1000 1200 0 2 4 x 104 Layer#1 0 200 400 600 800 1000 1200 0 1 2 x 104 Layer#2 0 200 400 600 800 1000 1200 0 1 2 x 104 Layer#3Number of accepted models 0 200 400 600 800 1000 1200 0 1 2 x 104 Layer#4 Velocity (m/s) Figure 5-30. Velocity histograms from all accepted models for the real data set at Newberry. 126

PAGE 127

127 0 5 101520253035 Distance (m) -12 -10 -8 -6 -4 -2 0Depth (m) 300 400 500 600 700 800 900 1000 1100 0 5101520253035 Distance (m) -10 -8 -6 -4 -2 0Depth (m) 300 500 700 900 1100 1300 1500 1700 1900 2100 2300 2500 2700 2900S-wave velocity (m/s) P-wave velocit y ( m/s ) a) b) Figure 5-31. Newberry test site: a) S-wave tomogram from the full wave field, and b) P-wave tomogram from refraction travel times

PAGE 128

CHAPTER 6 CLOSURE 6.1 Summary of Findings Two inversion techniques using simulated annealing have been developed to invert travel times and full waveforms for wave velocity profil es. They are first tested on synthetic data and then applied to real test data. The inverted profil es from the real test data are compared to other independent test results such as those of CPT, SPT, geotechnical borings, and cross-hole test to verify the feasibilities of the techniques for geotechnical engineering applications. 6.1.1 Inversion Technique Using Travel Times The technique using travel times is based on an extremely fast finite -difference solution of the Eikonal equation to compute the first-arriva l time through the velocity models by the multistencils fast marching method. The core of the simulated annealing, the Metropolis sampler, is applied in cascade with respect to shots to si gnificantly reduce computer time. This technique has been applied to surface data, and to combined surface and borehole data. For the cases of only surface data, the following findings have been derived: Tested on synthetic data, the presented techni que produces inverted re sults consistent to those from commercial codes including Ge oCT-II, SeisImager, and Rayfract, and it performs better than those codes in ca ses of sharp contrasts in velocity. Applied to real test data, the presented tec hnique produces inverted results that appear consistent to those from i nvasive tests including CPT, SPT, and geotechnical borings. Using surface travel times, the technique can well characterize normal profiles that increase in velocity with depth. Although the presented technique requires more computer time than the three mentioned commercial packages, it has several advantages First, this inversion technique does not depend on the initial model and becomes impor tant in regions where prior information about subsurface profiles is not available. Sec ond, simulated annealing provides a suite of final models clustering around the global solu tion and having comparable least-squared error. This provides an inversion result by av eraging all of these models to mitigate the influence of noise and the non-uniqueness of the inversion so lutions. Last, the technique also enables to determine the uncertaintie s associated with inverted results. 128

PAGE 129

As a typical limitation of surface-based refractio n methods, the presented technique fails to detect anomalies such as low-velocity zone s and embedded cavities. It also fails to characterize velocity structure below a few meters from the top of the half space (bottom layer), because the fastest rays only travel w ithin these few meters at the top of the half space regardless of how very large the geophone spread is used. For the cases of combined surface and borehol e data, the following findings have been derived: Tested on synthetic combined data, the pres ented technique well characterizes velocity structures at depth near the borehole. Presence, location, ge neral shape, and tr ue velocity of embedded low-velocity zone s are well characterized. Also tested on synthetic combined data, the presented technique requires just a few shot locations on the ground surface within a few meters around a borehole, surface geophones within this few meters, and a string of borehole geophones to reliably assess the subsurface properties at depth near the borehole. Applied to real test combined data, the pr esented technique produces credible inverted velocity profiles. A comparison of tomograms utilizing the combined data against tomograms developed using just the surface data suggests that significant additional resolution of inverted profiles at depth are obt ained with the addition of a borehole. It is also found that the quantitativ e uncertainties associated w ith the inverted profiles are significantly reduced when adding a borehole. Employed for site characterization of deep foundation design, refr action tomography using combined surface and borehole data provides credib le information of material at the socket and partially detect anomalies near the sock et. This becomes very important because the material at and near the sock et often carries a majority of load from foundations. The inversion results of the combined data, in cluding inverted prof iles and associated uncertainties, enable to characterize spatial va riability in geotechnica l engineering physical parameters of subsurface formations useful in the design of deep foundations. This will be particularly useful in implementing the ne w LRFD design methodology that can explicitly account for spatial variability and uncertainty in design parameters. 6.1.2 Inversion Technique Using Full Waveforms The technique using full waveforms is based on a finite-difference so lution of 2-D elastic wave equation in the time distance domain. It uses full information of elastic wave fields to increase resolution of inversion results, especially dealing with reverse models. It also employs a global inversion technique, simulated annealing, to invert the full wave fields for near-surface velocity profiles. The technique is first tested on many different synthetic data sets created from 129

PAGE 130

challenging reverse models with hi gh-velocity and low-velocity laye rs at different depths. Then it is applied to experimental data sets, and the inversi on results are compared to invasive test results including cross-hole, SPT N-value, and material log, or compared to results of independent refraction tests. Th e following findings about the t echnique have been derived: Tested on 1-D synthetic models, the presente d technique performs well on reverse models and shows superiority over the traditional techniques us ing Rayleigh wave velocity dispersion. Tested on 2-D synthetic models, the presente d technique shows that 2-D profiles with multi-linear interfaces of a few meter segments can be characterized using a full wave field from only one source. Applied to real full waveform data, the presen ted technique produces inverted results that appear consistent to those from invasive tests including SPT, cross-hole test, and geotechnical borings. The inverted results fr om the full waveform data also appear consistent to the inverted results from the refraction travel times. Besides its superiority over th e traditional methods in dealin g with reverse profiles, the presented technique possibly replaces the se ismic refraction tomography to characterize normal profiles (increasing veloci ty with depth) because it requi res much less effort in both site testing and manua l data processing. The presented technique does not depend on th e initial model and becomes important in regions where prior information about subsurface profiles is not available. Simulated annealing provides a suite of final models clustering around the global solution and having comparable least-squared error. This provides an inversion re sult by averaging all of these models to mitigate the influence of noi se and the non-uniqueness of the inversion solutions. 6.2 Conclusions Based on the findings outlined above, the following conclusions are made. First, the technique using travel times can well characterize normal profiles (increasing in velocity with depth) using only surface data and it also well characterize pr ofiles with embedded anomalies using combined surface and borehole data. The i nverted results from travel times appear consistent to those from inva sive tests including CPT, SPT, and geotechnical borings. Second, the technique using full waveforms can accurate ly characterize challe nging reverse models 130

PAGE 131

including highand low-velocity layers with only surface data. The inverted results from full waveforms are consistent with those from invasive tests incl uding SPT, cross-hole test, and geotechnical borings. The inverted re sults also appear consistent to the inverted results from the refraction travel times. Last, both inversion techni ques do not depend on the initial model and become important in regions where prior informa tion about subsurface profiles is not available. Simulated annealing provides a suite of fina l models clustering around the global solution and having comparable least-squared error. This provides an inversion result by averaging all of these models to mitigate the influence of noise and the non-uniqueness of the inversion solutions. 6.3 Recommendations The following recommendations are suggested after reviewing all of the findings and conclusions previ ously discussed: The developed inversion tec hniques should be applied in commercial software. Full waveform inversion technique should be investigated more. For example, using data from multiple shots to increase the resolution of inverted profiles, in such a case, the medium can be divided into many small cells and velocity of each cell is determined independently. Guidelines for array design and depth of inve stigation of the full waveform technique are needed. The global inversion techniques could be follo wed by a local inversion technique in order to save computer time. 131

PAGE 132

LIST OF REFERENCES Burger, H. R. (1992). Explorati on geophysics of the shallow subs urface. Prentice-Hall, Inc., Englewood Cliffs, New Jersey. Carpenter, P. J., Higuera-Diaz, I. C., Thompson, M. D., Atre, S., and Mandell, W. (2003). Accuracy of seismic refraction tomography codes at karst sites. Geophysical Site Characterization: Seeing Beneath the Sur face, Proceedings of a Symposium on the Application of Geophysics to Engineeri ng and Environmental Problems, San Antonio, Texas, April 6-10 832-840. Chopp, D.L. (2001). Some improvements on fast marching method. SIAM J. Scientific Computing 23(1), 230-244. Clayton, R. and Engquist, B. (1977). Absorbing boundary condition for acoustic and elastic waves. Bull. Seimol. Soc. Am. 67(6), 1529-1540. Dijkstra, E.W. (1959). A note on two pr oblems in connection with graphs. Numer. Math. J. 1, 269-271. Doyle H. (1995). Seismology. J. Wiley & sons, Chichester. Foti S. (2000). Multistation methods for geot echnical characterizati on using surface waves. Ph.D. dissertation, Politecnico di Torino. Hassouna, M.S. and Farag, A.A. (2007). Multisten cils fast marching methods: a highly accurate solution to the eikonal equa tion on cartesian domains, IEEE Trans. Patt. Anal. Mach. Int. 29(9), 1-12. Hiltunen, D. R. and Cramer, B. J. (2008). Appli cation of seismic refraction tomography in karst terrane. Journal of Geotechnical and Ge oenvironmental Engineering, 134(7), 938-948. Ingber, L. (1989). Very fast simulated reannealing. Math. Comput. Modeling 12(8), 967-993. Ingber, L. (1993). Simulated ann ealing: practice versus theory. Math. Comput. Modeling 18(11), 29-57. Jin, X., Luke B., and Calderon-Macias C. (2009). Role of forward model in surface-wave studies to delineate a buri ed high-velocity layer. Journal of Environmental and Engineering Geophysics, 14(1), 1-14. Kim, S. (1999). ENO-DNO-PS: A stable, second-order accuracy eikonal solver. Soc. Exploration Geophysicists 1747-1750. Li J. (2008). Study of surface wave methods for deep shear wave velocity profiling applied in the upper mississippi embayment. Ph.D. dissertation, University of Missouri Columbia. 132

PAGE 133

Louie, J. N. (2001). Faster, better, shear-wave velocity to 10 0 meters depth from refraction microtremor arrays. Bulletin of Seismological Society of America 91(2), 347-364. Metropolis, N., Rosenbluth A., Rosenbluth M., Te ller A., and Teller E., (1953). Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 21(6), 1087-1092. Moser, T.J. (1991). Shortest path calculation of seismic rays. Geophysics, 56, 59-67. Nakanishi, J. and Yamaguchi, K. (1986). A numerical experiment on non-linear image reconstruction from first arrival times for two-dimensional island structure. J. Phys. Earth 34, 195-201. Nazarian, S. (1984). In situ determination of elastic moduli of soil deposits and pavement systems by spectral-analysis-of-surface-wa ves method. Ph.D. Dissertation, The University of Texas at Austin. Nichols, D.E. (1996). Maximum energy travel times calculated in seismic frequency band, Geophysics, 61, 253-263. Oneil, A. (2003). Full-waveform reflectivity for modeling, inversion, and appraisal of seismic surface wave dispersion in shallo w site investigations. Ph.D. dissertation, The University of Western Australia. Park, C. B., Miller, R. D., and Xia, J. ( 1999), Multi-channel anal ysis of surface wave (MASW), Geophysics, 64(3), 800-808. Pratt, R. G. (1998). Seismic waveform inversion in the frequency domain, part 1: theory and verification in a physical scale model. Geophysics, 64(3), 888-901. Pullammanappallil, S.K. and Louie, J.N. (1994). A generalized simulated annealing optimization for inversion of first arrival time. Bull. Seimol. Soc. Am ., 84(5), 1397-1409. Quigley, T. P. (2006). Ground proving seismi c refraction tomography (SRT) in laterally variable karstic limestone terrain. Masters Thesis, University of Florida. Redpath, B. B. (1973). Seismic refraction exploration for engineering site investigations. Technical Report TR E -4 U. S. Army Engineer Waterways Experiment Station Explosive Excavation Research Laboratory, Livermore, California, May. Sambridge M. and Mosegaard K. (2000). Mont e carlo methods in geophysical inverse problems, Reviews of Geophysics 40(3). Sharma, S.P. and Kaikkonen, P. (1998). Two-d imensional non-linear inve rsion of VLF-R data using simulated annealing. Geophysics J, Int. 133, 649-668. Sheehan, J.R., Doll, W.E., and Mandell, W.A. ( 2005). An evaluation of methods and available software for seismic refraction tomography analysis. Journal of Environmental and Engineering Geophysics 10, 21-34. 133

PAGE 134

Sheen, D. H., Tuncay, K., Baag, C. E., and Ortoleva1, P. J. (2006). Time domain Gauss Newton seismic waveform inve rsion in elastic media. Geophysical Journal International 167, 1373-1384. Shipp, R. M. and Singh, S. C. (2002). Two-d imensional full wavefiel d inversion of wideaperture marine seismic streamer data. Geophysical Journa l International 151, 325-344. Sen, M.K. and Stoffa, P.L. (1991). Nonlinear one-dimensional seismic waveform inversion using simulated annealing, Geophysics, 56, 1624-1638. Sen, M.K. and Stoffa, P.L. (1995). Global o ptimization methods in geophysical inversion, Adv. Explor. Geophysics, 4, Elsevier Sci., New York. Sethian, J.A. (1996a), Fast marching set level method for monotonically advancing fronts, Pro. Nat. Acad. Sci. 93. Sethian, J.A. (1999). Level set methods and fast marching methods, second ed. Campridge Univ. Press. Tokimatsu, K., Tamura, S., and Kojima, H. (1992). Effects of multiple modes on Rayleigh wave dispersion characteristics. J. Geotech. Eng. 118 (2), 1529. Thomson W.T. (1950). Transmission of elastic waves through a stratif ied solid medium. J. Applied Physics, 21(1), 89-93. Van Trier, J. and Symes, W. (1991). Upwind fi nite-difference calculati on of travel times. Geophysics, 56, 812-821. Vidale, J.E. (1988). Finite-diffe rence travel time calculation. Bull. Seis. Soc. Am. 78, 20622076. Virieux, J. (1986). P-SV wave propagation in heterogeneous me dia: velocity-stress finitedifference method. Geophysics, 51(4), 889-901. Zhang, J. and Toksoz, M.N., (1998). Non linear refraction traveltime tomography. Geophysics, 63(5), 1726-1737. 134

PAGE 135

BIOGRAPHICAL SKETCH Khiem Tat Tran was born in 1978 in Thanh Ho a, Vietnam, and remained in Thanh Hoa until he graduated from Lam Son High School in 1996. He enrolled in Hanoi University of Civil Engineering, and graduated with a Bachelor of Science in civil engineer ing in spring 2001. He decided that it would be most beneficial to gain a few year s of work experience before continuing on with graduate studies so he worked for five years in Vietnam until he moved to United States for studying. He enrolled at the Univ ersity of Florida in Ga inesville, FL in August of 2006 where he worked as a graduate research assistant under Dr. Dennis Hiltunen. He completed his Master of Scie nce in May of 2008 and Doctor of Philosophy in August of 2010. 135