This study uses an ensemble Kalman method for near-surface seismic site characterization of 154 network earthquake monitoring stations in California to improve the resolution of S-wave velocity (V_S) and P-wave velocity (V_P) profiles—up to the resolution depth—coupled with better quantification of uncertainties compared to previous site characterization studies at this network. These stations were part of the Yong et al. site characterization project, with selected stations based on future recordings of ground motions that are expected to exceed 10 per cent peak ground acceleration in 50 yr. To estimate V_S and V_P from experimental dispersion data, Yong et al. investigated these stations using linearized (local search and iteration) routines, and Yong et al. later studied a subset of these stations using nonlinear (global search and optimization) routines. In both studies, the selection of model parameters—that is, discretization of the V_S and V_P profiles with only five fixed thickness layers—was mainly based on trial and error. In contrast, this paper uses an approximate Bayesian method to assimilate experimental dispersion data and sequentially update an ensemble of particle estimates that span the V_S and V_P parameter spaces. Doing so, we systematically determine the most probable profiles conditioned on the experimental dispersion data, the introduced noise levels, and a priori knowledge in the form of physical constraints. We consider two configurations to discretize the soil depth from the surface to half of the maximum discernible wavelength obtained from the experimental dispersion data, namely refined and coarse models, and two initial models for each configuration to study solution multiplicity. Our results suggest that using the refined model for the top surface layers improves the resolution of near-surface site characteristics and the model’s success rate in capturing dispersion data at high frequencies. All models result in similar V_S but distinct V_P profiles, with increasing uncertainty at deeper layers, suggesting that the fundamental mode of Rayleigh wave dispersion data is not adequate to constrain the P-wave velocity profile and the S-wave velocity close to the resolution depth.