To complement velocity distributions, seismic attenuation provides additional important information on fluid properties of hydrocarbon reservoirs in exploration seismology, as well as temperature distributions, partial melting, and water content within the crust and mantle in earthquake seismology. Full waveform inversion (FWI), as one of the state-of-the-art seismic imaging techniques, can produce high-resolution constraints for subsurface (an)elastic parameters by minimizing the difference between observed and predicted seismograms. Traditional waveform inversion for attenuation is commonly based on the standard-linear-solid (SLS) wave equation, in which case the quality factor (Q) has to be converted to stress and strain relaxation times. When using multiple attenuation mechanisms in the SLS method, it is difficult to directly estimate these relaxation time parameters. Based on a time domain complex-valued viscoacoustic wave equation, we present an FWI framework for simultaneously estimating subsurface P wave velocity and attenuation distributions. Because Q is explicitly incorporated into the viscoacoustic wave equation, we directly derive P wave velocity and Q sensitivity kernels using the adjoint-state method and simultaneously estimate their subsurface distributions. By analyzing the Gauss-Newton Hessian, we observe strong interparameter crosstalk, especially the leakage from velocity to Q. We approximate the Hessian inverse using a preconditioned L-BFGS method in viscoacoustic FWI, which enables us to successfully reduce interparameter crosstalk and produce accurate velocity and attenuation models. Numerical examples demonstrate the feasibility and robustness of the proposed method for simultaneously mapping complex velocity and Q distributions in the subsurface.