Viscous seismic wave propagation simulation using the fractional order equation has attracted much recent attention. However, conventional finite-difference (FD) methods of the fractional partial difference equation adopt a global difference operator to approximate the fractional derivatives, which reduces the computational efficiency dramatically. To improve the efficiency of the FD method, we have developed a reasonable truncated stencil pattern by strict mathematical derivation and adopted an asymptotic local FD (ALFD) method. Theoretical analysis and numerical results indicate that the ALFD method is accurate and efficient. In fact, our numerical results illustrate that the numerical solution solved by the ALFD method has a maximum relative error not exceeding 0.014% compared to the reference solution (applied to a finely meshed computational domain). The computation speed of ALFD is also significantly faster than that of the original FD method. The computational time of the three ALFD methods satisfying a different preset accuracy is only approximately 2.71%, 1.26%, and 0.78% of that of the original fractional wave equation FD method. The ALFD method provides a useful tool for viscoelastic seismic wavefield propagation simulation.