The dispersion and attenuation characteristics of Rayleigh waves in a non-uniform viscoelastic half-space with a covering layer are studied in this paper. The half-space and the covering layer are modeled by the fractional-order Zener viscoelastic solid. Compared with the traditional integer-order viscoelastic model, the viscoelastic model with fractional-order derivative is more flexible in describing the complicated history-dependent mechanical behavior. Legendre and Laguerre orthogonal polynomials are used to approximate the displacement field of surface waves. The rectangular window function is used to merge the surface boundary conditions. The complicated problem of solving the complex wave number in the complex region is ultimately transformed into an eigenvalue problem. The correctness of the method is verified by comparing our outcomes with those in the literature. It is found that the present spectrum method avoids the complicated iterative process of root-finding and the pseudo-root and root-missing problem compared to the traditional root-finding method. Moreover, the uneven gradient distribution of the cover layer and half-space has opposite effects on the propagation of Rayleigh waves at low and high frequencies, and the same influence applies to fractional-order effects. The present study has important theoretical significance and application potential in the geophysical exploration and earthquake engineering.