Abstract: Elastic least-squares reverse time migration (ELSRTM) offers higher resolution, better amplitude balancing, fewer crosstalk noises, and broader bandwidth compared to traditional elastic reverse time migration (ERTM). However, most current ELSRTM methods are implemented in the data domain, typically requiring multiple iterations, each iteration involving significant computational cost. Additionally, the perturbation model of the reflection coefficients obtained from conventional ELSRTM differs apparently from the imaging results of ERTM. Therefore, we propose an image-domain vector elastic least-squares reverse time migration (ID-VELSRTM) method based on the point spread function (PSF). This method aims to achieve imaging results that are consistent with the geological structures obtained via ERTM while improving image quality and resolution. The inversion is performed in the model space of the image domain, offering significant computational efficiency advantages over data-domain methods. Using the PSF from optics as a localized approximation of the Hessian matrix, we construct an image-domain objective function for the ID-VELSRTM method. In the least-squares inversion framework, space-variant deconvolution combined with the fast iterative shrinkage-thresholding algorithm (FISTA) is applied to the ERTM imaging results to obtain the final PP and PS images. Numerical experiments with layered, graben, and SEG/EAGE salt models demonstrate the effectiveness of the ID-VELSRTM method. Compared to ERTM, ID-VELSRTM achieves better amplitude balancing, fewer crosstalk noises, a broader wavenumber range, and higher imaging resolution.