Elastic least-squares reverse time migration (ELSRTM) is a powerful tool to retrieve high-resolution subsurface images of the earth’s interior. By minimizing the differences between synthetic and observed data, ELSRTM can improve spatial resolution and reduce migration artifacts. However, conventional ELSRTM methods usually assume constant density models, which cause inaccurate amplitude performance in resulting images. To partially remedy this problem, we have developed a new ELSRTM method that considers P- and S-wave velocity and density variations. Our method can simultaneously obtain P- and S-wave velocity and density images with enhanced amplitude fidelity and suppressed parameter crosstalk. In addition, it can provide subsurface elastic impedance images by summing the inverted velocity images with the density image. Compared with the conventional ELSRTM method, our method can improve the quality of final images and provide more accurate reflectivity estimates. Numerical experiments on a horizontal reflector model and a Marmousi-II model demonstrate the effectiveness of this method.