A special indirect boundary element method (IBEM) is proposed to investigate the waves scattering of plane P, SV, and SH waves by a 3D alluvial basin embedded in a multilayered half‐space. The new IBEM, which uses half‐space Green’s functions for uniformly distributed loads acting on an inclined plane as its fundamental solutions, has the merits of (1) excellent capability of dealing with the stratification of the basin and the external half‐space, (2) without the problem of singularity due to fictitious distributed loads being directly applied on the real boundaries, and (3) good adaptability to complex models with trapezoidal or triangular elements being used to discretize the boundaries. The validity and accuracy of the new method are verified by comparing its results with those in the literature. To illustrate the general applicability and efficiency of the new method further, 3D alluvial basins of varying shapes, depths, and sedimentary sequences embedded in a single layer overlying a homogeneous half‐space are numerically studied. Numerical results show that the basin’s shape, depth, and sedimentary sequence all have significant impact on the ground seismic responses; the incident angle also has noticeable effects on the surface motion, and these effects are more prominent at the observation points along the incident direction of the plane waves; for the case of layered model, the displacement spectral amplification is affected by the eigenmodes of the vibrations of the layers, both inside and outside the basin.