Near-surface characterization has now gained significance among exploration geophysicists, and many methods are being proposed to retrieve the 2D structures of shallow soils. Because most of these methods are based on the modal inversion of the surface waves, they can only be applied to laterally homogeneous or smoothly heterogeneous soil models. We have developed a time-domain waveform inversion method for 2D near-surface exploration that offers an alternative approach to existing surface-wave techniques for layered soils with a flat surface. Our method directly fits the input Rayleigh waveforms to retrieve the 2D soil structure without need of any modal identification, allowing the inversion of soil models that can be challenging with modal-inversion-based approaches. In our method, the forward problem formulated in the time domain is based on a 2.5D staggered-grid finite-difference scheme to simulate the P-SV wavefield; soil modeling was achieved by dividing soil layers into specific number of blocks with discontinuous interfaces. The inversion strategy depends on attributing suitable values for the interface depth and S-wave velocity for each block to reconstruct a numerical soil model that fit the input waveforms. Because we cannot know the source signature during data acquisition, source deconvolution by a reference station is applied to observed and calculated waveforms to make a waveform inversion free of the source signature. Numerical experiments revealed that our method was able to sufficiently reconstruct soil structures with strong lateral velocity gradient or soils with a blind layer in noisy environments, using a single source and reasonable number of receivers. We also applied this method to real waveform data, and we succeeded in obtaining good correlation between the inverted 2D soil model and the existing borehole data.