We present a semianalytical numerical method to simulate SH‐wave propagation in heterogeneous multiple layers with smooth property variations. Wave propagation in such piecewise heterogeneous media will dominate the reflection and transmission at interfaces, whereas the scattered waves by volume heterogeneities inside each geological layer are superimposed on the boundary waves to cause fluctuations in amplitude and phase. The incident, boundary‐scattering, and volume‐scattering waves are accurately superposed through the generalized Lippmann–Schwinger integral (GLSI) equation. The multilayered heterogeneous media with irregular interfaces are sliced into horizontal and heterogeneous slabs, with each discretized into a series of volume elements. The traditional generalized reflection–transmission matrices (GRTMs) method is extended to solve the GLSI equation for each heterogeneous slab in the frequency–wavenumber domain by incorporating the boundary–volume integral equation numerical techniques. The global GRTM is obtained by integrating all the slab’s GRTMs in terms of the boundary condition across adjacent slabs. The wavefield at an arbitrary depth can be computed in a recursive relation from a source term. The accuracy of the method is tested by comparing it with other full‐wave solutions to the SH problem. To show the applicability of the method, we calculate synthetic seismograms to demonstrate the significant impact of volume heterogeneities on seismic responses in layered heterogeneous media.