The problem of reconstructing underground obstacles from near-field, surface seismic measurements is investigated within the framework of a linear sampling method. Although the latter approach has been the subject of mounting attention in inverse acoustics dealing with far-field wave patterns in infinite domains, there have apparently not been any attempts to apply this new method to the interpretation of near-field elastic wave forms such as those relevant to the detection of subterranean objects. Aimed at closing this gap, a three-dimensional inverse analysis of elastic waves scattered by an obstacle (or a system thereof), manifest in the surface ground motion patterns, is formulated as a linear integral equation of the first kind whose solution becomes unbounded in the exterior of the hidden scatterer. To provide a comprehensive theoretical foundation for this class of imaging solutions, generalization of the linear sampling method to near-field elastodynamics and semi-infinite domains is highlighted in terms of its key aspects. A set of numerical examples is included to illustrate the performance of the method. On replacing the featured elastodynamic half-space Green function by its free-space counterpart, the proposed study is directly applicable to infinite media as well.