Abstract: Seismic wave propagation in viscoacoustic media results in amplitude attenuation and phase distortion, which degrade energy focusing and final resolution of imaging. Current reverse time migration (RTM) for viscoacoustic media is predominantly based on fractional-order differential equations. These approaches operate in the frequency-wavenumber domain to decouple amplitude attenuation from phase dispersion, thereby compensating for amplitude loss while correcting phase distortion. The pseudospectral method is routinely employed for numerical computation; however, its efficiency in a 3D setting is insufficient for field data application. To address this limitation, this paper adopts the generalized standard linear solid (GSLS) model to achieve a multi-parameter, high-accuracy characterization of the quality factor (Q) within a defined frequency range and thereby enable the implementation using an integer-order differential operator. The computation is performed using a high-order finite-difference scheme, which offers flexibility in computational domain decomposition and is well-suited for GPU acceleration, thereby achieving an effective balance between imaging quality and computational efficiency. Application to synthetic and field data demonstrates that this method successfully implements viscoacoustic RTM in the time-space domain with enhanced imaging accuracy and resolution via effective amplitude compensation and phase correction. The proposed method exhibits the capability for industrial-scale application to field data.