Compared with other imaging algorithms (e.g., ray-based and one-way wave equation), reverse time migration (RTM) based on two-way wave equation exhibits great superiority, especially in dealing with steeply dipping structures. However, imaging with conventional single-component seismic data may become imperfect in some complicated structures (e.g., gas clouds). The elastic reverse time migration based on the elastodynamic equation using multi-component seismic data can extract PP and PS reflectivity containing subsurface information, thus it can be more consistent with characteristics of elastic wave propagation in the real earth's medium, and resulting seismic images can more accurately characterize the subsurface. To begin with, we employ the first-order velocity-stress equations to implement extrapolation of elastic vector wavefields, and separate P- and S-wavefields by computing the divergence and curl operators of the extrapolated velocity vector wavefields. Then, imaging data with pure wave modes can be computed by applying the source normalized cross-correlation imaging condition, thus avoiding the crosstalk between the unseparated wave modes. To address the polarity reversal problem of converted images, we propose an alternative method in the shot domain. We also develop an efficient method that reconstructs source wavefields in the reverse time direction to save storage in CPU and avoid large input/output for elastic reverse time migration. During the forward modeling, the method only saves velocity vector wavefields of all time interval within an efficient absorbing boundary and total wavefields in the last time interval. When we extrapolate the receiver wavefields in reverse time direction, simultaneously, we reconstruct the total source wavefields by the saved wavefields. Numerical examples with a graben and Marmousi2 models show that the polarity reversal correction method works well and that elastic reverse time migration can generate accurate images for complicated structures.