We have developed a 3D elastic full-waveform inversion (FWI) method for geotechnical site characterization. The method is based on a solution of 3D elastic-wave equations for forward modeling to simulate wave propagation and a local optimization approach based on the adjoint-state method to update the model parameters. The staggered-grid finite-difference technique is used to solve the wave equations together with implementation of the perfectly matched layer condition for boundary truncation. Seismic wavefields are acquired from geophysical testing using sensors and sources located in uniform 2D grids on the ground surface, and they are then inverted for the extraction of 3D subsurface wave velocity structures. The capability of the presented FWI method is tested on synthetic and field data sets. The inversion results from synthetic data indicate the ability of characterizing laterally variable low- and high-velocity layers. Field experimental data were collected using 96 receivers and a propelled energy generator to induce seismic wave energy. The field data result indicates that the waveform analysis was able to delineate variable subsurface soil layers. The seismic inversion results are generally consistent with invasive standard penetration test -values, including identification of a low-velocity zone.