Next Article in Journal
Luminous Flux Utilization of Static Birefringent Fourier Transform Imaging Spectrometer with Zoomable Spectral Resolution
Next Article in Special Issue
Utilizing Multilevel Modeling to Measure Neighborhood Dynamics and Their Impact on House Prices
Previous Article in Journal
Synergic Effects of Nano Additives on Mechanical Performance and Microstructure of Lightweight Cement Mortar
Previous Article in Special Issue
Seismic Interferometry Method Based on Hierarchical Frequency Fusion and Its Application in Microtremor Survey
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Multi-Mode Imaging of Ambient Background Noise for Karst Detection in the Limestone Area Based on Frequency-Bessel Transform

1
Wuhan Center of China Geological Survey (Geosciences Innovation Center of Central South China), 69 Optics Valley Road, Wuhan 430205, China
2
Hubei Key Laboratory of Paleontology and Geological Environment Evolution, Wuhan Center of China Geological Survey, 69 Optics Valley Road, Wuhan 430205, China
3
Hubei Key Laboratory of Marine Geological Resources, China University of Geosciences, 388 Lumo Road, Wuhan 430074, China
*
Author to whom correspondence should be addressed.
Appl. Sci. 2023, 13(8), 5135; https://doi.org/10.3390/app13085135
Submission received: 28 February 2023 / Revised: 15 April 2023 / Accepted: 19 April 2023 / Published: 20 April 2023
(This article belongs to the Special Issue State-of-the-Art Earth Sciences and Geography in China)

Abstract

:
In response to the challenges of karst geophysical exploration in an environment with strong external interference, this paper proposes a new method, namely the frequency-Bessel transform method, for extracting multi-order dispersion curves of surface waves from background noise to characterize karst. The observation noise data of the Wuhan karst development area are used as an example, where the dolomitic limestone and limestone mixed with dolomite of the Jialing River Formation of the middle lower Triassic are widely developed in the observation area. The frequency-Bessel transform method involves performing a Bessel integral transformation on the cross-correlation coefficient of background noise in the frequency domain. Firstly, by synthesizing theoretical noise data and comparing it with the spatial autocorrelation method—which is currently the main method for extracting the fundamental dispersion curve of surface waves—it is verified that the frequency-Bessel transform method can extract the higher-mode dispersion curve. Then, by taking the actual measured single-point noise data as an example, the effect of applying the frequency-Bessel transform to the actual noise data is tested, and the inversion of the fine structure of the strata by the addition of higher-mode dispersion, the use of the damped least squares inversion method, and the joint inversion of fundamental and higher-mode dispersion curves are analyzed. The higher-mode dispersion curve of Rayleigh surface wave extracted by the frequency-Bessel transform is much clearer, and the 2D shear wave velocity structure profile obtained from inversion explains the karst development area, karst strip area, and thickness of the Quaternary overburden. The inferred results match with the actual borehole data. Multi-mode imaging of background noise based on the frequency-Bessel method can be applied to depict karst in complex backgrounds, and has significant potentiality in the field of ambient seismic noise tomography, providing a new idea and method for karst detection in near-surface engineering.

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.
In this paper, the frequency-Bessel function transform is used to extract the multi-order dispersion curves from background noise data from stations in the karst area, and then the damped least squares algorithm with mature methods, fast operation speed, and high inversion accuracy is used to invert the multi-order dispersion curves of the Rayleigh surface wave. With the introduction of the new method, the theoretical and practical effect of higher-mode dispersion in the extraction of background noise data by the frequency-Bessel function method is discussed. The damped least squares algorithm for the inversion of the shear wave velocity structure improves the ability to depict karst in complex near-surface environments, and provides a new idea and method for detecting karst in near-surface engineering.

2. Materials and Methods

2.1. Theory of the Frequency-Bessel Transform

The spatial autocorrelation theory has been derived in detail in many references [30,31,32], which is considered mature for extracting the fundamental dispersion cure; this section mainly introduces the derivation process of the frequency-Bessel transform. Previous research [33] proved that the Fourier transform of the cross-correlation function between the noise data records observed at any two stations x 1 and x 2 on the surface at the same time is proportional to the imaginary part of Green’s function; the difference between the two is only in amplitude, and the relationship between them can be expressed as follows:
C ( r , ω ) A · Im G z z ( r , z = 0 , ω )
where C ( r , ω ) is the Fourier transform of the cross-correlation function; A is a constant; r is the distance between the two stations; ω is the angular frequency; and G z z ( r , z = 0 , ω ) is the surface Green’s function of the vertical component. Therefore, the theoretical formula of the frequency-Bessel function transform is derived as follows:
0 + C ( r , ω ) J 0 ( k r ) r d r A · Im g z ( z = 0 , ω , k )
where J 0 ( k r ) is the first type of zero-order Bessel function; g z ( r , ω , k ) is the kernel function of Green’s function; and k is the wave number. According to the nature of the kernel function, when the wave number k is at the dispersion point, the dispersion point of the kernel function approaches infinity, the dispersion point is the singularity of the kernel function, and the value corresponds to the dispersion curve [34,35]. The integral operation on the left side of Equation (2) can be simplified as the sum of finite terms; C ( r , ω ) is expressed as:
C ( r , ω ) a i + b i r
where coefficient b i = C ( r i , ω ) C ( r i 1 , ω ) / ( r i r i 1 ) ; coefficient a i = C ( r i 1 , ω ) + C ( r i , ω ) C ( r i 1 , ω ) r i 1 / ( r i r i 1 ) ; and r i is the distance between stations of the cross-correlation spectrum C ( r i , ω ) in the frequency domain. According to the recursive property of Bessel function:
x J 0 ( x ) d x = x J 1 ( x ) x 2 J 0 ( x ) d x = x 2 J 1 ( x ) + x J 0 ( x ) J 0 ( x ) d x .
By using Equation (4), the integral term on the left side of the frequency-Bessel transform Equation (2) can be derived as [36]:
0 + C ( r , ω ) J 0 ( k r ) r d r i = 1 N a i c ω { J 1 ( ω c r ) } r i r i 1 + , b i c 2 ω 2 { k r 2 J 1 ( ω c r ) + r J 0 ( x ) } r i r i 1 b i c 2 ω 2 r i 1 r i J 0 ( ω c r ) d r .
According to Equation (5), the theoretical calculation of the frequency-Bessel function of the noise data is carried out, and the dispersion curves can be obtained by picking up the energy maximum in the dispersion diagram. The implementation process of this method mainly consists of the following aspects:
Firstly, for a given seismic observation array, the spectral cross-correlation functions for all possible station pairs are calculated after preprocessing the ambient seismic records. Secondly, the spectral CCFs are ordered as a function of the inter-station distance r i and subjected to the frequency-Bessel transform. The frequency-Bessel spectral map is calculated by the discrete summation formula for each given ω . Finally, with the help of image recognition algorithms, dispersion curves can be identified from images.
The frequency-Bessel transform is an effective method to extract multimode dispersion curves from ambient background noise. Several local and regional applications have demonstrated the applicability of the method approach [37].
With the high frequency dispersion information extracted by the frequency-Bessel method, we can obtain not only the short wavelength information of the fundamental mode, but also the long wavelength information of the higher mode of the surface wave [38]. Constrains on short and long wavelength information in the high frequency range can reduce the uncertainty of the inversion.

2.2. Test with Synthetic Background Noise Data

In order to test the effect of the frequency-Bessel transform in extracting multi-order dispersion curves from background noise data, we designed a theoretical model of four layers of horizontally layered media to synthesize noise data. The P-wave and S-wave velocity, density, and depth of each stratum are based on the borehole data collected in the paper. The specific parameters of physical properties are shown in Table 1.
By the simulation method proposed by Wang et al., in the paper, the center frequency of 6–10 Hz Reyker wavelet was selected as the source time function, and the discrete wavenumber method was used to randomly generate 1000 vertical noise signals. The source intensity was between 0.001 and 1, with a signal sampling interval of 6 ms, and acquisition time of 60 s. The observation system (Figure 1) adopted a linear arrangement of 60 stations for reception, represented by red dots with a spacing of 2.5 m. The randomly distributed noise is located between the inner and outer circles, represented by blue dots, with an inner circle radius of 500 m and an outer circle radius of 1500 m.
The synthesized noise data were segmented by 5 s, and the dispersion energy superposition was calculated by the spatial autocorrelation method and the frequency-Bessel transform algorithm, respectively (Figure 2). Among them, the dispersion spectrum calculated by the spatial autocorrelation (a) had a well-developed fundamental dispersion, mainly concentrated at 7–20 Hz, but there was no higher-mode dispersion information in the dispersion spectrum. The dispersion spectrum calculated by the frequency-Bessel transform (b) developed multi-order mode dispersion, with the dispersion of fundamental mode being concentrated at 7–20 Hz, and the higher-mode dispersion also being clear and continuous. It can be seen that the information of the dispersion spectrum provided by the frequency-Bessel transform is more abundant. The synthetic noise data also confirms the theoretical basis and computational effectiveness of the frequency-Bessel transform in extracting multi-order dispersion information.

3. Application to Field Data

3.1. Geological Background

The geotectonic structure of Wuhan City straddles two major primary tectonic units: the Yangtze landmass and the Qinling-Dabie orogenic belt. It is bounded by the Xiangfan-Guangji Fault, the southern side is the southeastern Eurasian Fold Belt of the lower Yangzte landmass, and the northern side is the Qinling arc-basin system of the Qinling-Dabie orogenic belt. Wuhan is dominated by a denudation landform and an alluvial lacustrine plain landform, the surface is widely distributed within a Quaternary loose accumulation layer, and the underlying bedrock is dominated by four sets of carbonate rocks.
The karst area in Wuhan is mainly distributed in strips in the Daye stratigraphic community of the lower Yangzte stratigraphic division in the southern Yangzte stratigraphic region, with hidden karst as the main form. The karst forms are mainly karst gaps, karst pores, and small-scale karst caves, with filling karst caves as the main form.
The main soluble rocks are dolomitic limestone and limestone intercalated with dolomite of the Jialing River Formation in the middle and lower Triassic, and limestone of Daye Formation in the lower Triassic. Meanwhile, there is also chert nodule limestone of the Qixia Formation in the lower Permian, thick layered dolomite and limestone of the Chuanshan Formation and Huanglong Formation in the upper Carboniferous, and claystone intercalated with limestone of the Hezhou Formation in the lower Permian. Among them, the highest degree of karst development is found in the Huanglong Formation and Chuanshan Formation in the upper Carboniferous.

3.2. Background Noise Data Acquisition

The point of noise data acquisition is located in the karst development area of Jinshui Street in Wuhan, China, with roads and villages evenly distributed around the acquisition point, and rich noise sources (Figure 3). The instrument used in this research was the three-component McSEIS-AT seismic station produced by Japan’s OYO company, with a main geophone frequency of 2 Hz. There are no specific requirements of the frequency-Bessel calculation for the observation system, which can be either randomly distributed or linearly arranged. A total of 7 three-component seismic stations were arranged on the site to observe in a linear arrangement (Figure 4), with a station spacing of 5 m, sampling rate of 4 ms, and acquisition time of 30 min. After the completion of each observation duration, the overall arrangement was moved forward by 1 station spacing and the collection continued for 30 min.
A total of 45 sets of background noise data were collected this time. The observation system was designed to facilitate the comparison of the actual dispersion curves calculated by the spatial autocorrelation method and the frequency-Bessel transform.

3.3. Dispersion Curve Calculation

The collected noise data needed to be preprocessed, and the current preprocessing methods [39] include de-instrument response, mean removal, linear trend removal, band-pass filtering, spectral whitening, etc. In the paper, we mainly used mean removal, linear trend removal, and one-bit normalization. The specific processing is as follows: the noise record was segmented as a small window signal, with a general value of 5 s to 30 s, and the specific window value requires parameter testing; the spatial cross-correlation coefficients between each two sets of station pairs within the hour window were calculated, and the cross-correlation spectrum was averaged, which is the process of improving imaging resolution; by integrating the superimposed correlation spectrum, the frequency phase velocity energy spectrum of the noise signal was calculated, and the fundamental and higher-mode dispersion curves in the dispersion spectrum were further extracted.
The spatial autocorrelation method and frequency-Bessel transform method were used to calculate the dispersion spectrum. The noise was divided according to the 20 s time window, and the frequency range was set between 1 Hz and 20 Hz, with a frequency interval of 0.1 Hz. The phase velocity calculation range was 10–1000 m/s, with a phase velocity interval of 5 m/s.

3.4. Multi-Mode Dispersion Curve Comparison

Figure 5 shows the dispersion energy spectrum calculated by the working stations, including the spatial autocorrelation method and the frequency-Bessel function transform. There are some differences in dispersion modes between the two methods. The first higher-mode curve calculated by the frequency-Bessel transform method for partial frequencies above 10 Hz is more obvious, while by the spatial autocorrelation method the higher-mode data are not detected and a phenomenon of mode crossover with fundamental dispersion is exhibited. The frequency-Bessel transform method has a larger dispersion distribution range of the fundamental mode in frequencies below 10 Hz than the spectral range calculated by the spatial autocorrelation method. The signal-to-noise ratio of dispersion energy is also significantly high, and the dispersion energy is continuous, which makes it easier to identify and extract the dispersion points. When the frequency is higher than 10 Hz, the energy carried by the fundamental mode is much lower than that of the higher-mode mode. When the fundamental signal is disturbed by the higher-mode dispersion or body wave, it is necessary to use the higher-mode mode information to calculate the shear wave velocity structure.
Therefore, the comparison diagram of the dispersion energy spectrum also indicates that the higher-mode dispersion curve cannot be extracted from background noise only by the spatial autocorrelation method, resulting in multiple solutions and instability of inversion by using only the fundamental mode.

4. Joint Inversion

4.1. Damped Least Squares Inversion

In recent years, there have been many studies on the inversion methods of Rayleigh surface wave fundamental dispersion curves, among which the damped least squares algorithm has been widely used due to its fast operation speed and high inversion accuracy. In the paper, the damped least squares method is used to invert the fundamental and higher-mode dispersion curves of the surface wave extracted from background noise, and to calculate the shear wave velocity structure of the strata. The multi-order joint inversion by the damped least squares algorithm is similar to the fundamental inversion method. The difference is that the higher-mode mode data are added in the inversion process, and the partial derivatives of the fundamental mode and the higher-mode mode to the shear velocity need to be calculated separately when calculating the Jacobi matrix [40]. Due to the low sensitivity of the Rayleigh wave phase velocity to density and P-wave velocity, only two parameters, namely S-wave velocity and layer thickness, are calculated in the inversion of the theoretical model in the paper.
The objective function of the damped least squares inversion [41] is defined as:
ϕ = J Δ x Δ b 2 W J Δ x Δ b 2 + α Δ x 2 2
where x = v s 1 , v s 2 , v s 3 , , v s n T is the n-dimensional length vector of shear wave velocity, b = b 1 , b 2 , b 3 , , b m T is the m-dimensional length vector of Rayleigh wave phase velocity, Δ x is the correction of shear wave velocity, Δ b is the difference between the initial value model response and the observed data, J is the Jacobi matrix with m rows and n columns, α is the damping factor, W is the weighting matrix, W = L T L , and L is the diagonal matrix, which represents the partial derivative of phase velocity with respect to shear wave velocity. According to the singular value decomposition method, the correction of shear wave velocity can be obtained as follows:
Δ x = V Λ 2 + α I 1 Λ U T d
where I is the identity matrix, d = L b , and d is the phase velocity matrix. The damping factor α controls the direction and convergence rate of the shear wave velocity correction, and the stability and convergence of the inversion process can be controlled by adjusting the α value. According to Equation (7), without considering the number of damping factor changes, the appropriate S-wave velocity correction can be calculated once by singular value decomposition, and the inversion process can be realized.

4.2. Multi-Order Joint Inversion of S-Wave Velocity Structure

Real strata are often not the ideal models of increasing velocity. When there is a low velocity or high velocity layer sandwiched between a layer in the middle, the energy distribution of each model of the surface wave will change, and the velocity curves of each model will also exhibit the phenomenon of mode crossover. Researchers [42] compared the dispersion curves extracted from the numerical simulation of the Rayleigh wave field with those calculated from the dispersion equation, and analyzed the energy distribution characteristics of the multi-order mode in the Rayleigh wave field. When dispersion of higher-mode energy dominates, the development of the fundamental surface wave cannot be observed in the dispersion energy spectrum. It is very necessary to introduce the joint inversion of fundamental and higher-mode surface wave dispersion. At present, the joint inversion of fundamental and higher-mode dispersion curves is mainly applied in the direction of the active source of the Rayleigh surface wave, and the inversion effect has been theoretically and practically verified. However, in the past, background noise dispersion curves were mostly inverted using the fundamental dispersion curve, and the joint inversion effect of measured data needs to be verified by more application cases.
In the paper, we take the dispersion curves of measured noise data in Figure 4 as an example. An initial number of 10 layers was adopted for the inversion model of the fundamental and higher-mode dispersion curves, and the values of S-wave velocity, P-wave velocity, density, layer thickness, and depth of the model layers were given. S-wave velocity and P-wave velocity were converted according to a certain Poisson’s ratio, and the layer thickness increment factor was 1.1. The minimum S-wave velocity value was 100 m/s, the maximum S-wave velocity value was 3000 m/s, the maximum number of iterations was 20, and the root mean square error limit was 0.10 m/s.
Figure 6a shows a monotonically increasing trend of shear wave velocity in the stratum of the fundamental mode dispersion inversion, which does not match the stratum structure of the test area. With the addition of the first higher-mode dispersion information, Figure 7a shows the S-wave velocity inversion in both shallow 12.4 m and deep 67.4 m, which is closer to the real stratum situation. Figure 6b shows the dispersion curve of measured noise data picked up in the fundamental mode and the theoretical dispersion curve calculated after the initial model iteration, and Figure 7b shows the dispersion curve of the inversion in the fundamental mode, and the first higher-mode fitting the measured data dispersion curve. It can be seen that the joint inversion of the multi-order dispersion with the addition of higher-mode dispersion has a higher resolution of the stratum and can distinguish more detailed layers. In the dispersion curve below the 3 Hz band, the multi-order mode fits the measured dispersion curve better than the fundamental mode does. In the dispersion curve above 15 Hz, the fundamental and higher-mode dispersion shows a trend of mode crossover, which is close to the occurrence of a missing stratum, or inversion in real strata.
The joint inversion of the fundamental and higher-mode dispersion of measured noise data through the frequency-Bessel transform results in a S-wave velocity structure that is more consistent with the sedimentary characteristics of real strata, laying a reliable theoretical basis for the calculation and interpretation of two-dimensional profiles, and at the same time providing a new method for the detection of karst development under the Quaternary overburden.

5. Results

The above section discusses the use of the frequency-Bessel transform algorithm to calculate the multi-mode dispersion energy spectrum of noise data, and this section mainly verifies the feasibility and practical application effect of this method by taking 45 sets of noise data measured in the Wuhan karst development area as an example.
Figure 8 shows the dispersion energy spectrum of nine positions, from point 0 to point 40 calculated at equal intervals, with seven seismic stations as a group. In the dispersion spectrum, there is a significant development of fundamental and higher-mode energy, with a wide range of dispersion energy distribution where the fundamental dispersion is mainly developed in the frequency band below 10 Hz, and the frequency band above 10 Hz is mainly of a higher-mode energy type. On the basis of the dispersion energy spectrum, the multi-mode dispersion curves can be extracted for inversion. In the paper, the known geological data were used to establish a constrained parameter model for stratum P-wave velocity, S-wave velocity, and density, and the damped least squares algorithm was used to iteratively invert the extracted multi-mode dispersion curves, obtaining 45 sets of data curves of stratum S-wave velocity changing with depth. Finally, a S-wave velocity profile was produced.
The survey line was located in the karst area, where karst is relatively developed, and the dissolution of the rock is relatively serious. The pores of karst after dissolution will be filled with a large amount of water, which is shown in the lower S-wave velocity value in the dissolved area, and not in the intact and undissolved area. As shown in Figure 9, the black curve marks the bottom interface of the Quaternary, with a maximum burial depth of 32 m and a minimum burial depth of 16 m. The S-wave velocity of the bottom interface is about 260 m/s. By tracking the value of S-wave velocity, the undulating shape of the bottom of the Quaternary can be sketched. The inversion profile marks four areas of low shear wave velocity anomalies, namely anomaly a at 25 m, anomaly b at 50 m, anomaly c at 100 m, and anomaly d at 100 m to 140 m. It can be inferred that the karst fissures in the underlying rock layers are relatively developed. At the same time, low-speed strips appear near the sections 155 m and 185 m, showing the characteristics of large burial depth in the longitudinal direction and small width in the transverse direction, which are different from the above-mentioned four anomalies. This can be assumed to be the location of karst strip development, and karst fissures are very developed.
According to the borehole data at the 110 m position of the S-wave velocity profile, it was revealed that within 4 m of the stratum, there is a miscellaneous fill layer, 4 m to 10 m is a silty clay layer, 10 m to 27 m a silty sand layer without water content, 27 m to 28 m lime dolomite (strongly weathered, mixed with purplish red gray dolomite, with a developed karst structure and numerous karst fissures forming pores in the rock layer with a pore diameter of 2–7 mm), and 28 m to 45 m gray-yellow lime dolomite with a developed karst structure and numerous dissolution fissures. It can be seen that the borehole data and the inverted shear wave velocity anomaly d in Figure 9 correspond well with each other, verifying the correctness of the interpretation.
By analyzing the distribution profile of the S-wave velocity in Figure 9, it can also better correspond to four velocity change layers: 80 m/s to 135 m/s is a clay layer within 10 m of the surface layer, 135 m/s to 260 m/s a silty sand layer, 260 m/s to 360 m/s a dissolved lime dolomite layer, and 350 m/s to 560 m/s underlying bedrock.

6. Discussion

After the frequency-Bessel method adopted in the paper was proposed, exciting results have been achieved in recent years in the large-scale imaging of the shear wave velocity structure of the Kanto Basin in Japan [43] and the three-dimensional shear wave velocity structure in the continental United States [44,45]. Although some scholars have made some attempts to study the frequency-Bessel method in detecting shallow geological structures, there have been no reports of combining the frequency-Bessel method with karst exploration. Compared to the background noise surface wave method, the use of the multi-channel surface wave method in detecting railway hazards [46,47] has achieved certain results. Due to the significant noise interference in the urban environment of the study area in the paper, it is difficult to collect high-quality raw data by conventional geophysical methods, including multi-channel surface wave technology. Therefore, surface wave analysis based on background noise will be very meaningful for karst detection.
The karst development area studied in the paper was previously explored by the high-density resistivity method, and the interpreted anomaly location of karst [48] is very consistent with the results of the joint inversion using multi-mode dispersion in the paper. However, the use of multi-mode dispersion information has a higher lateral resolution of karst and is more effective in vertical stratification, which verifies the feasibility and correctness of karst imaging based on multi-mode dispersion.
Due to the low detection efficiency and relatively high construction difficulty of the high-density resistivity method, the data are easily disturbed by external interference. In comparison, using the background noise as the signal source, the fundamental and higher-mode dispersion curves can be extracted through the frequency-Bessel transform, and then the S-wave velocity structure of strata can be inverted, making it easier to achieve in data acquisition and processing.
In addition, according to the interpretation results [49] of ground penetrating radar carried out in the work area, there are discontinuous stripes in the profile, accompanied by small arc-shaped reflections in some areas, and multiple reflections in the karst area, indicating that ground penetrating radar has a certain sensitivity to karst. However, compared to the 90 m depth of inversion through the frequency-Bessel transform in the paper, the abnormal location reflected by ground penetrating radar is shallow, and the detection effect is greatly affected by surface conditions [50]. Therefore, the method has strong applicability, does not require the arrangement of active sources like other geophysical methods, and can achieve good results in noise-rich urban environments.

7. Conclusions

In this study, the dispersion curves of multi-mode Rayleigh surface wave were extracted from an ambient background noise signal, and through the dispersion curves the shear velocity structure of karst in the Wuhan area was inverted. The conclusions are as follows:
(1)
Adding higher-mode dispersion information can improve the inversion fitting effect of the low-frequency band of the fundamental dispersion curve, and the joint inversion of the fundamental and higher-mode of the measured data improves the ability to identify the fine structure of the strata.
(2)
The shear wave velocity structure profile obtained from the inversion of background noise data in the limestone area of Wuhan depicts the range of the karst low-speed area, and the abnormal positions are in good agreement with the borehole data.
(3)
The thickness of the Quaternary overburden calculated by the frequency-Bessel algorithm is more accurate, and the undulation of the bottom interface of the Quaternary is closer to the real weathering state of the stratum.
To sum up, the frequency-Bessel transform method provides a new idea and method for accurately detecting overburden thickness and karst in near-surface engineering, and we can try to utilize this method extensively whether in surface wave research or in solving practical geological problems.

Author Contributions

Conceptualization, S.C. and F.C.; methodology, S.C.; software, S.C.; validation, D.L.; formal analysis, S.C.; investigation, D.L.; resources, S.C.; data curation, S.C. and J.X.; writing—original draft preparation, S.C.; writing—review and editing, F.C.; visualization, S.C.; supervision, D.L.; project administration, S.C.; funding acquisition, D.L. and J.X. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Hubei Key Laboratory Foundation from the Hubei Key Laboratory of Paleontology and Geological Environment Evolution, Wuhan Center of China Geological Survey, grant number PEL202212; the National Key Research and Development Program of China, grant number 2020YFC1512401; and the China Geological Survey Project, grant numbers DD20221729 and DD20221734.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

We greatly appreciate the comments and suggestions from two anonymous reviewers that significantly improved the quality of the manuscript. We are grateful to Shichuan Yuan for sharing the synthetic data used in the synthetic analysis and Chao Shen for the discussion about the effect of extracting multi-mode dispersion from field data. We also thank Lei Wang for his help in drawing the geological map of Wuhan.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Gómez-Ortiz, D.; Martín-Crespo, T. Assessing the risk of subsidence of a sinkhole collapse using ground penetrating radar and electrical resistivity tomography. Eng. Geol. 2012, 149–150, 1–12. [Google Scholar] [CrossRef]
  2. De Waele, J.; Gutiérrez, F.; Parise, M.; Plan, L. Geomorphology and natural hazards in karst areas: A review. Geomorphology 2011, 134, 1–8. [Google Scholar] [CrossRef]
  3. Epting, J.; Huggenberger, P.; Glur, L. Integrated investigations of karst phenomena in urban environments. Eng. Geol. 2009, 109, 273–289. [Google Scholar] [CrossRef]
  4. Amanatidou, E.; Vargemezis, G.; Tsourlos, P. Combined application of seismic and electrical geophysical methods for karst cavities detection: A case study at the campus of the new University of Western Macedonia, Kozani, Greece. J. Appl. Geophys. 2022, 196, 104499. [Google Scholar] [CrossRef]
  5. Fernandes, A.L.; Medeiros, W.E.; Bezerra, F.H.R.; Oliveira, J.G.; Cazarin, C.L. GPR investigation of karst guided by comparison with outcrop and unmanned aerial vehicle imagery. J. Appl. Geophys. 2015, 112, 268–278. [Google Scholar] [CrossRef]
  6. Martínez-Moreno, F.J.; Galindo-Zaldívar, J.; Pedrera, A.; Teixido, T.; Ruano, P.; Peña, J.A.; González-Castillo, L.; Ruiz-Constán, A.; López-Chicano, M.; Martín-Rosales, W. Integrated geophysical methods for studying the karst system of Gruta de las Maravillas (Aracena, Southwest Spain). J. Appl. Geophys. 2014, 107, 149–162. [Google Scholar] [CrossRef]
  7. Verdet, C.; Sirieix, C.; Marache, A.; Riss, J.; Portais, J.-C. Detection of undercover karst features by geophysics (ERT) Lascaux cave hill. Geomorphology 2020, 360, 107177. [Google Scholar] [CrossRef]
  8. Foudili, D.; Bouzid, A.; Berguig, M.C.; Bougchiche, S.S.; Abtout, A.; Guemache, M.A. Investigating karst collapse geohazards using magnetotellurics: A case study of M’rara basin, Algerian Sahara. J. Appl. Geophys. 2019, 160, 144–156. [Google Scholar] [CrossRef]
  9. Lin, S.; Wang, W.; Jin, C.; Deng, X.; Liu, Z. Application and discussion of seismic CT in detailed karst detection: A case of Shenzhen metro line 14. Sci. Technol. Eng. 2019, 19, 18–23. [Google Scholar]
  10. Gan, F.; Han, K.; Lan, F.; Chen, Y.; Zhang, W. Multi-geophysical approaches to detect karst channels underground—A case study in Mengzi of Yunnan Province, China. J. Appl. Geophys. 2017, 136, 91–98. [Google Scholar] [CrossRef]
  11. QADY, G.E.; Hafez, M.; Abdalla, M.A.; Ushijima, K. Imaging subsurface cavities using geoelectric tomography and ground-penetrating radar. J. Cave Karst Stud. 2005, 67, 174–181. [Google Scholar]
  12. Leucci, G. Evaluation of karstic cave stability using integrated geophysical methods_2003. GeoActa 2003, 2, 75–88. [Google Scholar]
  13. Solbakk, T.; Fichler, C.; Wheeler, W.; Lauritzen, S.-E.; Ringrose, P. Detecting multiscale karst features including hidden caves using microgravimetry in a Caledonian nappe setting: Mefjell massif, Norway. Nor. J. Geol. 2018, 98, 359–378. [Google Scholar] [CrossRef]
  14. Wang, K.; Qian, J.; Zhang, H.; Gao, J.; Bi, D.; Gu, N. Seismic imaging of mine tunnels by ambient noise along linear arrays. J. Appl. Geophys. 2022, 203, 104718. [Google Scholar] [CrossRef]
  15. Jin, C.; Lin, S.; Wang, J.; Zhou, H.; Cheng, M. Estimation of Shallow Shear Velocity Structure in a Site with Weak Interlayer Based on Microtremor Array. Appl. Sci. 2022, 13, 185. [Google Scholar] [CrossRef]
  16. Tün, M.; Pekkan, E.; Özel, O.; Guney, Y. An investigation into the bedrock depth in the Eskisehir Quaternary Basin (Turkey) using the microtremor method. Geophys. J. Int. 2016, 207, 589–607. [Google Scholar] [CrossRef]
  17. Zhou, X.; Lin, J.; Zhang, H.; Jiao, J. Mapping extraction dispersion curves of multi-mode Rayleigh waves in microtremors. Chin. J. Geophys. 2014, 57, 2631–2643. [Google Scholar]
  18. Pan, L.; Chen, X.; Wang, J.; Yang, Z.; Zhang, D. Sensitivity analysis of dispersion curves of Rayleigh waves with fundamental and higher modes. Geophys. J. Int. 2019, 216, 1276–1303. [Google Scholar] [CrossRef]
  19. Luo, Y.; Xia, J.; Liu, J.; Liu, Q.; Xu, S. Joint inversion of high-frequency surface waves with fundamental and higher modes. J. Appl. Geophys. 2007, 62, 375–384. [Google Scholar] [CrossRef]
  20. Xia, J.; Miller, R.D.; Park, C.B.; Tian, G. Inversion of high frequency surface waves with fundamental and higher modes. J. Appl. Geophys. 2003, 52, 45–57. [Google Scholar] [CrossRef]
  21. Wang, J.; Wu, G.; Chen, X. Frequency-Bessel Transform Method for Effective Imaging of Higher-Mode Rayleigh Dispersion Curves From Ambient Seismic Noise Data. J. Geophys. Res. Solid Earth 2019, 124, 3708–3723. [Google Scholar] [CrossRef]
  22. Hu, S.; Luo, S.; Yao, H. The Frequency-Bessel Spectrograms of Multicomponent Cross-Correlation Functions From Seismic Ambient Noise. J. Geophys. Res. Solid Earth 2020, 125, e2020JB019630. [Google Scholar] [CrossRef]
  23. Beaty, K.S.; Schmitt, D.R.; Sacchi, M. Simulated annealing inversion of multimode Rayleigh wave dispersion curves for geological structure. Geophys. J. Int. 2002, 151, 622–631. [Google Scholar] [CrossRef]
  24. Feng, S.; Takeshi, S.; Hiroaki, Y. Effectiveness of multi-mode surface wave inversion in shallow engineering site investigations. Explor. Geophys. 2005, 36, 26–33. [Google Scholar] [CrossRef]
  25. Axe, T. An Unsupervised Wavelet Transform Method for Simultaneous Inversion of Multimode surface waves. J. Environ. Eng. 2005, 10, 287–294. [Google Scholar]
  26. Cai, W.; Song, X.; Yuan, S.; Hu, Y. Fast and stable Rayleigh-wave dispersion-curve inversion based on particle swarm optimization. Oil Geophys. Prospect. 2018, 53, 25–34. [Google Scholar] [CrossRef]
  27. Fu, D.; Liu, J.; Zhou, L.; Xu, H.; Liao, J.; Chen, S.; Guo, D. Inversion of multimode Rayleigh-wave dispersion curves of soft interlayer based on Bayesian theory. Chin. J. Geotech. Eng. 2015, 37, 321–329. [Google Scholar]
  28. Wang, Y.; Song, X.; Zhang, X. Research on nonlinear inversion of seismic surface waves based on artificial neural network algorithm. Oil Geophys. Prospect. 2021, 56, 979–991. [Google Scholar] [CrossRef]
  29. Yu, D.; Song, X.; Zhang, X.; Zhao, S.; Cai, W. Rayleigh wave dispersion inversion based on grasshopper optimization algorithm. Oil Geophys. Prospect. 2019, 54, 288–301. [Google Scholar] [CrossRef]
  30. Aki, K. Space and time spectra of stationary stochastic waves with special reference to microtremors. Bull. Earthq. Res. Inst. 1957, 35, 415–456. [Google Scholar]
  31. Okada, H. Theory of efficient array observations of microtremors with special reference to the SPAC method. Explor. Geophys. 2006, 37, 73–85. [Google Scholar] [CrossRef]
  32. You, Z.; Xu, P.; Ling, S.; Du, Y.; Zhang, R.; Yao, J.; Zhang, H. Estimation of shallow subsurface S-wave velocity structure in urban area based on inversion of apparent dispersion curve. J. Geophys. Eng. 2020, 17, 940–955. [Google Scholar] [CrossRef]
  33. Sanchez-Sesma, F.J. Retrieval of the Green’s Function from Cross Correlation: The Canonical Elastic Problem. Bull. Seismol. Soc. Am. 2006, 96, 1182–1191. [Google Scholar] [CrossRef]
  34. Dai, w.; Pan, L.; Wang, J.; Yang, Z.; Chen, X. Application of frequency-Bessel transform method in shallow exploration of the beach of Chaohu lake. Comput. Tech. Geophys. Geochem. Explor. 2021, 43, 290–295. [Google Scholar]
  35. Yang, Z.; Chen, X.; Pan, L.; Wang, J.; XU, J.; Zhang, D. Multi-channel analysis of Rayleigh waves based on the Vector Wavenumber Tansform Method (VWTM). Chin. J. Geophys. 2019, 62, 298–305. [Google Scholar]
  36. Li, X.; Chen, X.; Yang, Z.; Wang, B.; Yang, B. Application of hingher-order surface waves in shallow exploration: An example of the Suzhou river. Chin. J. Geophys. 2020, 63, 247–255. [Google Scholar]
  37. Zhan, W.; Pan, L.; Chen, X. A widespread mid-crustal low-velocity layer beneath Northeast China revealed by the multimodal inversion of Rayleigh waves from ambient seismic noise. J. Asian Earth Sci. 2020, 196, 104372. [Google Scholar] [CrossRef]
  38. Xi, C.; Xia, J.; Mi, B.; Dai, T.; Liu, Y.; Ning, L. Modified frequency–Bessel transform method for dispersion imaging of Rayleigh waves from ambient seismic noise. Geophys. J. Int. 2021, 225, 1271–1280. [Google Scholar] [CrossRef]
  39. Bensen, G.D.; Ritzwoller, M.H.; Barmin, M.P.; Levshin, A.L.; Lin, F.; Moschetti, M.P.; Shapiro, N.M.; Yang, Y. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements. Geophys. J. Int. 2007, 169, 1239–1260. [Google Scholar] [CrossRef]
  40. Luo, Y.; Xia, J.; Liu, J.; Liu, Q. Joint inversion of fundamental and higher mode Rayleigh waves. Chin. J. Geophys. 2008, 51, 242–249. [Google Scholar]
  41. Xia, J.; Miller, R.D.; Park, C.B. Estimation of near-surface shear-wave velocity by inversion of Rayleigh waves. Geophysics 1999, 56, 691–700. [Google Scholar] [CrossRef]
  42. Shao, G.; Li, Q.; Wu, H. Dispersion curves and mode energy distribution of Rayleigh wave based on wavefield numerical simulation. Oil Geophys. Prospect. 2015, 50, 306–315. [Google Scholar] [CrossRef]
  43. Wu, H.; Chen, X.; Pan, L. S-wave velocity imaging of the Kanto basin in Japan using the frequency-Bessel transformation method. Chin. J. Geophys. 2019, 62, 3400–3407. [Google Scholar]
  44. Wu, G.x.; Pan, L.; Wang, J.n.; Chen, X. Shear Velocity Inversion Using Multimodal Dispersion Curves From Ambient Seismic Noise Data of USArray Transportable Array. J. Geophys. Res. Solid Earth 2020, 125, e2019JB018213. [Google Scholar] [CrossRef]
  45. Li, Z.; Zhou, J.; Wu, G.; Wang, J.; Zhang, G.; Dong, S.; Pan, L.; Yang, Z.; Gao, L.; Ma, Q.; et al. CC-FJpy: A Python Package for Extracting Overtone Surface-Wave Dispersion from Seismic Ambient-Noise Cross Correlation. Seismol. Res. Lett. 2021, 92, 3179–3186. [Google Scholar] [CrossRef]
  46. Tang, T.; Xiao, X.; Li, Q. The application of transient Rayleigh wave to detecting Railway karst disaster. Chin. J. Eng. Geophys. 2019, 16, 339–350. [Google Scholar]
  47. Yang, Y.; Zhu, D. Application of multi-channel surface wave method based on CMPCC processing technology in karst exploration. Chin. J. Eng. Geophys. 2020, 17, 559–566. [Google Scholar]
  48. Liu, D.; Xu, J.; Liu, L.; He, J.; Qi, X.; Chen, S. Application of the integrated geophysical methods in the fine exploration of karst collapses: A case study of Wuhan City. Geol. Explor. 2022, 58, 865–874. [Google Scholar]
  49. He, J.; Liu, L.; Li, Q.; Liu, D.; Chen, B.; Zhang, A.; Zhao, Y. Techniques for detecting underground space in hidden karst region: Taking Wuhan as an example. Hydrogeol. Eng. Geol. 2020, 47, 47–56. [Google Scholar] [CrossRef]
  50. Zhao, Y.; Wang, J.; Zhang, W. Applied analysis of ground penetrating radar in Wuhan karst prospecting. Coal Geol. China 2019, 31, 108–112. [Google Scholar]
Figure 1. Observation system of synthetic background noise.
Figure 1. Observation system of synthetic background noise.
Applsci 13 05135 g001
Figure 2. Dispersion spectrograms calculated by different methods. (a) Dispersion spectrum calculated by the spatial autocorrelation (SPAC) method. (b) Dispersion spectrum calculated by the frequency-Bessel method.
Figure 2. Dispersion spectrograms calculated by different methods. (a) Dispersion spectrum calculated by the spatial autocorrelation (SPAC) method. (b) Dispersion spectrum calculated by the frequency-Bessel method.
Applsci 13 05135 g002
Figure 3. The geological map and location of survey lines in the study area.
Figure 3. The geological map and location of survey lines in the study area.
Applsci 13 05135 g003
Figure 4. Schematic layout of field data acquisition. (a) Diagram of station observation system. (b) Wave diagram of data collected by 7 seismic stations.
Figure 4. Schematic layout of field data acquisition. (a) Diagram of station observation system. (b) Wave diagram of data collected by 7 seismic stations.
Applsci 13 05135 g004
Figure 5. Comparison diagram of dispersion energy spectrum of measured data. (a) Dispersion spectrum of the spatial autocorrelation method. (b) Dispersion spectrum of the frequency-Bessel transform method.
Figure 5. Comparison diagram of dispersion energy spectrum of measured data. (a) Dispersion spectrum of the spatial autocorrelation method. (b) Dispersion spectrum of the frequency-Bessel transform method.
Applsci 13 05135 g005
Figure 6. Inversion results of fundamental mode. (a) S-wave velocity structure of fundamental mode inversion. (b) Dispersion curve of measured noise data picked up in fundamental mode and the theoretical dispersion curve calculated after the initial model iteration.
Figure 6. Inversion results of fundamental mode. (a) S-wave velocity structure of fundamental mode inversion. (b) Dispersion curve of measured noise data picked up in fundamental mode and the theoretical dispersion curve calculated after the initial model iteration.
Applsci 13 05135 g006
Figure 7. Joint inversion results of fundamental and higher-mode curves. (a) S-wave velocity structure of fundamental and higher-mode inversion. (b) Dispersion curve of the inversion in the fundamental and first higher-mode fitting the measured data dispersion curve.
Figure 7. Joint inversion results of fundamental and higher-mode curves. (a) S-wave velocity structure of fundamental and higher-mode inversion. (b) Dispersion curve of the inversion in the fundamental and first higher-mode fitting the measured data dispersion curve.
Applsci 13 05135 g007
Figure 8. Multi-mode dispersion energy spectrum of measured data.
Figure 8. Multi-mode dispersion energy spectrum of measured data.
Applsci 13 05135 g008
Figure 9. Shear-wave velocity profile of strata by multi-mode inversion.
Figure 9. Shear-wave velocity profile of strata by multi-mode inversion.
Applsci 13 05135 g009
Table 1. Model parameters.
Table 1. Model parameters.
Layer NoDepth (m)P-Wave Velocity (m/s)S-Wave Velocity (m/s)Density
(g/cm3)
102251301.72
244502601.91
3106583801.97
42715939202.87
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Chen, S.; Liu, D.; Cheng, F.; Xu, J. Multi-Mode Imaging of Ambient Background Noise for Karst Detection in the Limestone Area Based on Frequency-Bessel Transform. Appl. Sci. 2023, 13, 5135. https://doi.org/10.3390/app13085135

AMA Style

Chen S, Liu D, Cheng F, Xu J. Multi-Mode Imaging of Ambient Background Noise for Karst Detection in the Limestone Area Based on Frequency-Bessel Transform. Applied Sciences. 2023; 13(8):5135. https://doi.org/10.3390/app13085135

Chicago/Turabian Style

Chen, Song, Daohan Liu, Fei Cheng, and Junjie Xu. 2023. "Multi-Mode Imaging of Ambient Background Noise for Karst Detection in the Limestone Area Based on Frequency-Bessel Transform" Applied Sciences 13, no. 8: 5135. https://doi.org/10.3390/app13085135

APA Style

Chen, S., Liu, D., Cheng, F., & Xu, J. (2023). Multi-Mode Imaging of Ambient Background Noise for Karst Detection in the Limestone Area Based on Frequency-Bessel Transform. Applied Sciences, 13(8), 5135. https://doi.org/10.3390/app13085135

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop