Near-surface complexity is characteristic of rugged topographies and strong volume heterogeneities, which significantly affects seismic data recorded at the free surface. An efficient method is presented here to quantitatively analyze the near-surface seismological complexity. We use a boundary-volume element numerical method to discretize the generalized Lippmann-Schwinger integral equation formulated in terms of volume scattering and boundary scattering. The integral equation technique provided sufficient accuracy to model rough topographies, large-scale structures, and volume heterogeneities in a straightforward and efficient manner. Rather than solving the resultant global matrix equation at a very high computational cost, certain characteristic parameters of the scattering amplitude matrix were estimated and used as indicators to assess the near-surface seismological complexity. Applications to numerous examples validated the method's applicability, reliability, and efficiency, whereby the resultant matrix characteristic parameters varied consistently with the near-surface complexity in geometric structures, volume heterogeneities, and computational frequencies.