The purpose of this work is to generalize the frozen Gaussian approximation (FGA) theory to solve the 3-D elastic wave equation and use it as the forward modelling tool for seismic tomography with high-frequency data. FGA has been previously developed and verified as an efficient solver for high-frequency acoustic wave propagation (P wave). The main contribution of this paper consists of three aspects: (1) We derive the FGA formulation for the 3-D elastic wave equation. Rather than standard ray-based methods (e.g. geometric optics and Gaussian beam method), the derivation requires to do asymptotic expansion in the week sense (integral form) so that one is able to perform integration by parts. Compared to the FGA theory for acoustic wave equation, the calculations in the derivation are much more technically involved due to the existence of both P and S waves, and the coupling of the polarized directions for SH and SV waves. In particular, we obtain the diabatic coupling terms for SH and SV waves, with the form closely connecting to the concept of Berry phase that is intensively studied in quantum mechanics and topology (Chern number). The accuracy and parallelizability of the FGA algorithm is illustrated by comparing to the spectral element method for 3-D elastic wave equation in homogeneous media. (2) We derive the interface conditions of FGA for 3-D elastic wave equation based on an Eulerian formulation and the Snell’s law. We verify these conditions by simulating high-frequency elastic wave propagation in a 1-D layered Earth model. In this example, we also show that it is natural to apply the FGA algorithm to geometries with non-Cartesian coordinates. (3) We apply the developed FGA algorithm for 3-D seismic wave-equation-based traveltime tomography and full waveform inversion, respectively.