Most conventional Reverse Time Migration (RTM) method for Ground Penetrating Radar (GPR) data under undulating surface condition employ Finite Difference Time Domain (FDTD) method for wavefield reverse time extrapolation. The staircase approximation inherent in Yee-grid FDTD discretization of undulating surfaces and complex interfaces generates geometric discretization errors, which consequently degrade the resolution of RTM imaging. To overcome this, this paper proposed a RTM imaging method for GPR data with undulating surfaces that utilizes the Finite Element Time Domain (FETD) with triangular mesh discretization. The FETD with triangular mesh discretization and non-split Complex Frequency Shifted Perfectly Matched Layer (CFS-PML) boundary condition is employed for GPR wavefield simulation and reverse time extrapolation. And the zero time imaging condition is applied to obtain the final imaging results. Numerical experiments of RTM imaging both of synthetic undulating surface data and field data demonstrated that compared with the RTM imaging results of FDTD with Yee mesh, the target and interface beneath the undulating surface in the RTM imaging results of FETD with triangular mesh are more convergent to the true positions, achieves higher imaging resolution, and requires fewer mesh elements and nodes. Furthermore, the proposed RTM method can effectively mitigate the adverse effects of undulating surface on migration quality, generating imaging results that are more reliable for geological interpretation.