Multichannel Analysis of Surface Waves (MASW) is one of the most widely used techniques in environmental and engineering geophysics to determine shear-wave velocities and dynamic properties, which is based on the elastic layered system theory. Wave propagation in the Earth, however, has been recognized as viscoelastic and the propagation of Rayleigh waves presents substantial differences in viscoelastic media as compared with elastic media. Therefore, it is necessary to carry out numerical simulation and dispersion analysis of Rayleigh waves in viscoelastic media to better understand Rayleigh-wave behaviors in the real world. We apply a pseudospectral method to the calculation of the spatial derivatives using a Chebyshev difference operator in the vertical direction and a Fourier difference operator in the horizontal direction based on the velocity-stress elastodynamic equations and relations of linear viscoelastic solids. This approach stretches the spatial discrete grid to have a minimum grid size near the free surface so that high accuracy and resolution are achieved at the free surface, which allows an effective incorporation of the free surface boundary conditions since the Chebyshev method is nonperiodic. We first use an elastic homogeneous half-space model to demonstrate the accuracy of the pseudospectral method comparing with the analytical solution, and verify the correctness of the numerical modeling results for a viscoelastic half-space comparing the phase velocities of Rayleigh wave between the theoretical values and the dispersive image generated by high-resolution linear Radon transform. We then simulate three types of two-layer models to analyze dispersive-energy characteristics for near-surface applications. Results demonstrate that the phase velocity of Rayleigh waves in viscoelastic media is relatively higher than in elastic media and the fundamental mode increases by 10-16% when the frequency is above 10. Hz due to the velocity dispersion of P and S waves. ?? 2011 Elsevier Ltd.