Elastic wave reverse time migration (RTM) and least squares reverse time migration (LSRTM) are reliable methods for imaging subsurface geological structures. However, existing methods based on RTM and LSRTM do not apply to linear slip fault interfaces. Therefore, developing a mathematical approach for imaging linear slip fault interfaces is necessary. Linear slip interfaces consider the shape of the interface and can be used to simulate fault interfaces. In this study, we combined the linear slip theory with RTM techniques and proposed an improved RTM and LSRTM imaging method for fault media. Initially, based on the explicit wave equation for fault media, we used a rotated staggered grid to obtain the forward propagation expression for fault media. This can provide a more accurate elastic wave numerical simulation. Using the adjoint state method, we derived an adjoint equation for fault media and a gradient expression of the Lamé coefficient and density parameters. This allowed the use of RTM and LSRTM techniques for obtaining images of the fault media. The simulation data test results show that the improved elastic wave numerical simulation scheme is effective and stable for linear slip fault media. The RTM method of using source-normalised cross-correlation imaging conditions can accurately image linear slip fault media. Analysing single-channel imaging and error curves of the compliance parameter showed that the LSRTM method can accurately determine the reflection coefficient of fault interfaces. Our research provides a new approach for fault media imaging studies.