Abstract: Boundary element method has unique advantages to solve elastic wave scattering problems, while the asymmetric and dense matrix equation greatly restrict the application of traditional boundary element in practical engineering. Therefore, based on the single-layer potential theory, combined with the fast multipole expansion technique, a new rapid indirect boundary element method is developed by using the Taylor series expansion of compressible and shear wave potential function for solving large scale elastic wave scattering problems efficiently and accurately. It is verified that this method greatly improves the computing speed and reduces the calculation storage. Then the elastic wave scattering problem for millions of DOFs can be solved effectively and accurately on an ordinary personal workstation. Finally, the scattering problem of P wave by ellipsoidal cavities group in elastic full-space is solved, and some important characteristics for the displacement and stress fields in steady state are discussed based on the numerical results. This method is especially suitable for solving the multiple scattering problems at low dimensionless frequency (ka<5.0).