## 1 Introduction  
Karst is a cavity formed by carbonate under the long-term action of groundwater, and is still a technical problem in the field of geological and engineering exploration due to the irregularity and unevenness of its spatial development. Karst features, such as joints, fissures, and cavities, directly affect the mechanical properties of bedrock and are prone to triggering geological disasters, posing potential risks to human construction projects [1; 2; 3]. In near-surface karst exploration, geophysical methods have been proved to be effective in many karst surveys [4]. Ground penetrating radar (GPR) has high resolution and is usually used for the imaging of shallow caves [5; 6]. Three-dimensional geological modeling using resistivity tomography (ERT) data can describe the spatial characteristics of underground karst and has been successfully applied in the Lascaux area (France) [7]. In Algeria, the geological hazards of karst collapse in the Sahara Basin were investigated by using magnetotelluric three-dimensional inversion technology [8]. Seismic wave CT (computerized tomography) technology was used to detect the fine structure of karst in the Shenzhen (China) Metro Line 14 [9], but there is a high requirement for the spacing between boreholes in CT technology. In addition, seismic refraction tomography [10; 11; 12] is also a common geophysical method for karst detection, and microgravity detection has been used to detect shallow karst caves and filling karst caves in Norway [13]. All the above work is based on geophysical exploration technology using active source data. Due to the complex weathering of near-surface strata and rich environmental noise, the background noise imaging method based on spatial autocorrelation is widely used at present [14; 15; 16]. However, the theoretical calculation formula of the SPAC method usually assumes that the Rayleigh surface wave of the fundamental mode is dominant [17], and that the dispersion information of the higher mode is not available. However, researchers have found that the dispersion signal of the higher mode is more sensitive to the formation parameters. A large number of studies [18; 19; 20] have shown that the joint calculation of surface wave dispersion signals using fundamental and higher-mode modes can reduce the multiplicity of inversion, improve inversion accuracy, and obtain a more accurate S-wave velocity structure. The energy of the higher-mode surface wave is weak, the amplitude varies greatly with the period, and the energy distribution of different modes of the surface wave in different frequency ranges is also different, so it is difficult to extract higher-mode dispersion information from background noise data by traditional methods. The new algorithm of the frequency-Bessel function transform can theoretically calculate the fundamental and higher-mode dispersion curves of the Rayleigh surface wave [21]. By linking the frequency-Bessel transform to the spatial autocorrelation coefficients and extending the method to the multi-component spatial autocorrelation tensor, the phase velocity dispersion curve of the multi-mode Love wave and the ellipticity of the higher-mode Rayleigh surface wave can be extracted [22]. At present, the nonlinear algorithms for the inversion of the shear velocity structure of Rayleigh surface wave dispersion curves mainly include the simulated annealing method [23], genetic algorithm [24], wavelet transform method [25], particle swarm optimization [26], Bayesian inversion method [27], artificial neural network algorithm [28], and new intelligent algorithm [29]. The inversion of Rayleigh surface wave dispersion curves is nonlinear, multi-parameter, and multi-extremum, and the effect of the above algorithms is not very significant in practical application.