With limited recording apertures, finite-frequency source functions, and irregular subsurface illuminations, traditional imaging methods have been insufficient to produce satisfactory reflectivity images with high resolution and amplitude fidelity. This is because most traditional imaging approaches are commonly formulated as the adjoint instead of the inverse operator with respect to the forward-modeling operator. In addition, intrinsic attenuation introduces amplitude loss and phase dispersion during wave propagation. Without considering these effects, migrated images might be kinematically and dynamically incorrect. We have developed a viscoacoustic least-squares reverse time migration (LSRTM) method based on a time-domain complex-valued wave equation. According to the Born approximation, we first linearized the viscoacoustic wave equation and derived a demigration operator. Then, using the complex-valued Lagrange multiplier method, we derived the adjoint viscoacoustic wave equation and corresponding sensitivity kernel. With the forward and adjoint operators, a linear inverse problem is formulated to estimate the subsurface reflectivity model. A total-variation regularization scheme is introduced to enhance the robustness of our viscoacoustic LSRTM, and a diagonal Hessian is used as the preconditioner to accelerate the convergence. Three synthetic examples are used to demonstrate that our approach enables us to compensate attenuation effects, improve imaging resolution, and enhance amplitude fidelity in comparison with the adjoint imaging method.