Abstract: Numerical simulation is an important way to study the characteristics of seismic wavefield. In this paper, we established the unsplit complex frequency shifted convolutional perfectly matched layer (CFS-CPML) absorbing boundary condition, and staggered grid and high precision finite difference were used to simulate Love wavefield in viscoelastic medium. And then we analyzed the influence of the number of relaxation mechanisms on the Q value fitting and verified the absorbing effect of CFS-CPML boundary conditions on the waves with large incident angles. The numerical results show that the requirements of the present theoretical wavefield study can be met with higher calculation efficiency if five relaxation mechanisms and finite difference with spatial four-order accuracy are selected for this method. At the same time, the phase velocity of dispersion curve gradually deviates from that of the theoretical curve with the decrease of quality factor, and the high frequency energy of each pattern also weakens. The study provided certain theoretical basis for high precision surface wave inversion method.