Abstract: Dispersion curves are widely used to constrain lithospheric structures. Existing methods are capable of extracting overtone dispersion curves from both ambient noise and earthquake data to some extent. However, when applied to earthquake records with broad azimuthal station coverage, the azimuthal-dependent radiation pattern of seismic sources may lead to inaccurate measurements of dispersion curves. To theoretically address this azimuthal effect introduced by the earthquake source radiation pattern, this study introduces the concept of surface integration into the frequency-Bessel (F-J) framework. Starting from the semi-analytical expressions of seismic wavefields, we derive a full set of dispersion spectrograms that include contributions from Rayleigh waves, Love waves, and leaking modes. Synthetic tests demonstrate that surface integration can improve the effects from source radiation pattern. To accommodate realistic array deployments with irregular station distributions, we implement a discretized version of surface integration based on triangular subarrays. Applications to both synthetic and real earthquake data validate our method’s ability in extracting comprehensive dispersion spectrograms from earthquake event records. Applications to both synthetic and real data validate the method's effectiveness in extracting the full set of dispersion spectrograms from earthquake records. By incorporating multiple spectrograms and combining with a mode separation strategy, the proposed method enables robust extraction of multimode dispersion curves. This extended F-J framework offers a flexible tool for analyzing diverse dispersion characteristics and holds promise for improving subsurface imaging by providing more complete and reliable structural constraints.