The anisotropic properties of subsurface media cause waveform distortions in seismic wave propagation, resulting in a negative influence on seismic imaging. In addition, wavefields simulated by the conventional coupled pseudo-acoustic equation are not only affected by SV-wave artifacts but are also limited by anisotropic parameters. We propose a least-squares reverse time migration (LSRTM) method based on the pure qP-wave equation in vertically transverse isotropic media. A finite difference and fast Fourier transform method, which can improve the efficiency of the numerical simulation compared to a pseudo-spectral method, is used to solve the pure qP-wave equation. We derive the corresponding demigration operator, migration operator, and gradient updating formula to implement the LSRTM. Numerical tests on the Hess model and field data confirm that the proposed method has a good correction effect for the travel time deviation caused by underground anisotropic media. Further, it significantly suppresses the migration noise, balances the imaging amplitude, and improves the imaging resolution.