Because of amplitude decay and phase dispersion of seismic waves, conventional migrations are insufficient to produce satisfactory images using data observed in highly attenuative geologic environments. We have developed a least-squares Gaussian beam migration method for viscoacoustic data imaging, which can not only compensate for amplitude decay and phase dispersion caused by attenuation, but it can also improve image resolution and amplitude fidelity through linearized least-squares inversion. We represent the viscoacoustic Green’s function by a summation of Gaussian beams, in which an attenuation traveltime is incorporated to simulate or compensate for attenuation effects. Based on the beam representation of the Green’s function, we construct the viscoacoustic Born forward modeling and adjoint migration operators, which can be effectively evaluated by a time-domain approach based on a filter-bank technique. With the constructed operators, we formulate a least-squares migration scheme to iteratively solve for the optimal image. Numerical tests on synthetic and field data sets demonstrate that our method can effectively compensate for the attenuation effects and produce images with higher resolution and more balanced amplitudes than images from acoustic least-squares Gaussian beam migration.