Beijing is an international metropolis, that is also an earthquake-prone city. The aims of this study are detailed quantifying and qualifying soil layer properties for an accurate seismic safety evaluation in the Beijing area. The time average shear-wave velocity in the first 30 m of subsoil, Vs30, is an important site parameter used in site response analysis, site classification, and seismic loss estimation. Mapping of Vs30over a city-scaled region is commonly done through proxy-based methods by correlating Vs30with geological or topographic information. In this paper, a geostatistical-based random field model is presented and applied to mapping Vs30over extended areas. This random field model is then coupled with Monte Carlo simulations to obtain an averaged Vs30map and its associated uncertainties. Unlike the traditional deterministic prediction model, this framework accounts for spatial variations of Vs30values and uncertainties, which makes the prediction more reliable. A total of 388 shear wave velocity measurements in the Beijing area are used to calculate Vs30values, from which the statistical and spatial properties for the random field realizations are inferred. New spatially correlated probabilistic Vs30maps for the Beijing area are then represented, and the effect of the maximum number of previously generated elements to correlate to in estimating Vs30maps is tested.