Abstract: Elastic least-squares reverse time migration (ELSRTM) offers higher resolution, better amplitude balancing, less crosstalk noise, and broader bandwidth compared to conventional elastic reverse time migration (ERTM). However, most ELSRTM methods are implemented in the data domain, typically requiring multiple iterations, each involving significant computational cost. Moreover, the reflectivity perturbation model obtained from ELSRTM differs remarkably from the imaging result of ERTM. To address this issue, we propose an image-domain vector elastic least-squares reverse time migration (ID-VELSRTM) method based on the point spread function (PSF). This method yields imaging results that are consistent with the geological structures obtained from conventional ERTM, while simultaneously improving imaging quality and resolution. The inversion is performed in the model space of the image domain, offering significant computational 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 ID-VELSRTM. Within the least-squares inversion framework, space-variant deconvolution is combined with the fast iterative shrinkage-thresholding algorithm (FISTA) to refine the ERTM imaging results and finally produce high-precision PP- and PS-wave images. Numerical experiments using a graben model and the SEG/EAGE salt model demonstrate the effectiveness of the ID-VELSRTM method. Compared to ERTM, ID-VELSRTM achieves better amplitude balancing, less crosstalk noise, broader wavenumber coverage, and significantly enhanced imaging resolution.