Dispersion properties of Rayleigh-type surface waves are widely used in environmental and engineering geophysics to image and characterize a shallow subsurface. In this paper, we numerically study the Rayleigh-type surface waves and their properties in 2D viscoelastic media. A finite difference method in a time-space domain is proposed, with an unsplit convolutional perfectly matched layer (C-PML) absorbing boundary condition. For two models that have analytical expressions of wave fields/dispersion curves, we calculate their wave fields and compare the analytical and numerical solutions to demonstrate the validity of this method. For the case where a medium has a high Poisson's ratio, say 0.49, traditional finite difference methods with a PML boundary condition are not stable when modeling Rayleigh waves but the proposed method is stable. For a laterally heterogeneous viscoelastic media model (Model 1) and a two-layer viscoelastic media model (Model 2) with a cavity, we use this method to obtain their corresponding Rayleigh waves. For several quality factors, the dispersion properties of these Rayleigh waves are analyzed. The results of Model 1 show that in a shallow subsurface, the phase velocity of a fundamental mode of the Rayleigh waves increases considerably with a quality factor $Q$ decreasing; the phase velocity increases with Poisson's ratio increasing. The results of Model 2 indicate that the energy of higher modes of the Rayleigh waves become strong when $Q$ decreases.