Amplitude energy attenuation and phase distortion of seismic waves caused by formation viscoelasticity reduce the resolution of reverse time migration (RTM) images. Q-RTM often is used to compensate the attenuation effects and improve the resolution of seismic imaging. However, serious high-frequency noise and tremendous amplitude will be produced during the wavefield extrapolation of Q-RTM, resulting in its inability to be imaged. Many Q-RTM algorithms solve the problem of instability through low-pass filtering in the wavenumber domain, but the method is less efficient in computation and has a truncation effect in the wavefield. We have developed a stable and efficient Q-RTM method, in which a regularization term was introduced into the viscoacoustic wave equation to suppress the high-frequency noise, and the finite-difference method was used to solve the viscoacoustic wave equation with a regularization term. We used the model example to visually demonstrate the instability of wavefield extrapolation in Q-RTM and compared the effect and computational efficiency of the two stabilization processing methods, low-pass filtering and regularization. Meanwhile, our method is not involved in solving the fractional derivatives by using the pseudo-spectral method, the computational efficiency also can be improved. We tested the Q-RTM approach on a simple layered model, Marmousi model, and real seismic data. The results of numerical examples demonstrated that the Q-RTM method can solve the problem of instability effectively and obtain a higher resolution image with lower computational cost.