The frozen Gaussian approximation (FGA) is an efficient solver for high frequency wave propagation. This work is to generalize the FGA to solve the 3-D elastic wave equation and use it as the forward modeling tool for seismic tomography with high-frequency initial datum. The evolution equation is derived by weak asymptotic analysis in conjunction with projecting onto an orthonormal frame; this is numerically verified and analytically proven to have same asymptotic error as the eigenfunction decomposition. Examples from seismology are given by forward modeling, solving an inverse problem and generating data sets to train neural networks.