The underwater multichannel analysis of surface wave (UMASW) is becoming an essential tool for surveying subbottom shear wave velocity. Current practice is limited to interpretation based on the fundamental mode. This paper investigates the full dynamic response of the underwater-multilayered structure under two different types of sources (impact and explosion). Fundamental solutions in the transformed domain, after applying both Fourier transform and Fourier–Bessel transform, are derived utilizing the global stiffness matrix method, analogous to a 1-D finite element approach. Solutions in the physical domain are then obtained using the fast and accurate Fourier series and Fourier–Bessel series approach proposed in this paper. The most attractive feature of this novel approach is that the expansion coefficients (or Love numbers) can be pre-calculated, saved, and repeatedly used for other field points. Through numerical analyses, we quantitatively investigate the effect of the water depth and source/receiver types/locations on different wave features (i.e., Scholte, fast-guided, and acoustic modes) from different perspectives (i.e., dispersion curve, Green’s function, frequency-velocity spectrum (FVS), and waveform). We find that (1) unlike the impact source, stronger acoustic waves can be produced by the explosive source, which sometimes causes difficulty in identifying the Scholte waves. (2) The acoustic and interface-guided waves exhibit distinct behaviors depending on the source and receiver locations. Scholte waves and fast-guided waves are weakened when the source and receivers are elevated far from the water/soil interface. Moreover, (3) the response of the Scholte wave can be enhanced by muting the direct acoustic waves. However, (4) the fast-guided wave may also arise in underwater surveys when the Vs of the underlying half-space is much higher than those of the upper layers, exerting a significant influence on the shallow water scenarios and posing challenges for correct mode identification. With the ability to model the entire wavefield (or FVS) that takes into account all propagation modes and the actual survey configuration, inversion can be performed by fitting the entire FVS. This approach eliminates the need to pick dispersion curves and provides a more accurate and stronger model constraint.