We develop a new automatic framework for the nondestructive multichannel analysis of surface waves that combines multimode dispersion spectrum matching and a finite-element method (FEM)-based inversion to enhance the accuracy of subsurface profiling in site investigation activities. This framework eliminates the need for manual identification of the Rayleigh wave energy component and multimode assignment, reducing dependence on operator experience and judgment. The dispersion spectrum is generated through an FEM model that simulates 2D seismic wave propagation, taking into account the actual acquisition layout and lateral variations in the subsurface. We introduce the Wasserstein distance (WD) for evaluating the difference between the observed and simulated spectra and incorporate Bayesian optimization for efficiently inverting S-wave velocity profiles. The effectiveness of our framework is demonstrated through synthetic data examples, and the superiority of the WD-based objective function is illustrated by comparing it with the conventional mean square error-based objective function. Subsequently, we conduct a field test on a reclaimed landfill to validate our framework. This test confirms the ability of the framework to retrieve multimode Rayleigh waves and demonstrates its effectiveness in providing high-resolution S-wave profiles of the shallow subsurface.