Abstract: Frequency behavior of Rayleigh waves in layered saturated poroelastic media is improtant for engineering geophysical prospecting and inversion of soil properties. The frequency equation of Rayleigh waves can be solved using the root searching algorithms in the conventional analysis method of the transfer or stiffness matrix. However, these algorithms are time-consuming in the complex domain and the root convergence is often poor. In this study, a rigid base is set at a certain depth of the half space which is estimated from the interested wavelength and the layers are discretized into a group of thin layers. The thin-layer model is then established to study the propagation behavior of Rayleigh waves in saturated poroelastic layered media under the plane strain conditions. In this model, the root searching problem is converted into the eigenvalue and eigenvector problem of matrices. The eigenvalues corresponding to the modes of Rayleigh waves can be sifted out from the calculated ones according to the attenuation of Rayleigh waves along depth. The frequency behavior, pore pressure and skeleton displacements of Rayleigh waves along depth can be calculated from the sifted eigenvalues and eigenvectors. The model developed is validated by comparing the numerical and the analytical results in the saturated poroelastic half space. The model is then used to study the frequency behavior of Rayleigh waves in several layered saturated porous half spaces. The method also provides some guidelines for investigating Rayleigh waves in layered media with more complex dynamic constitutive models.