Accurate Sparse Recovery of Rayleigh Wave Characteristics Using Fast Analysis of Wave Speed (FAWS) Algorithm for Soft Soil Layers

Abstract

This paper presents a novel fast analysis of wave speed (FAWS) algorithm from the waveforms recorded by a random-spaced geophone array based on a compressive sensing (CS) platform. Rayleigh-type seismic surface wave testing is excited by a hammer source and conducted to develop the phase velocity characteristics of the subsoil layers in Shenyang Metro line 9. Data are filtered by a bandpass filter bank to pursue the dispersive profiles of phase velocity at various frequencies. The Rayleigh-type surface-wave dispersion curve for the soil layers at each frequency is conducted by the ℓ₁-norm minimization algorithm of CS theory. The traditional frequency-wavenumber transform technique and in-site downhole observation are employed as the comparison of the proposed technique. The experimental results indicate the proposed FAWS algorithm has a good agreement with both the results of conventional even-spaced geophone array and the in-site measurements, which provides an effective and efficient way for accurate non-destructive evaluation of the surface wave dispersion curve of the soil.

1. Introduction

In recent decades, with rapid urban population growth, underground constructions which extend human activities into underground spaces have attracted increasing attention [1,2]. A metro system is a kind of underground construction in an urban region to perfect alternative transportation in order to overcome traffic congestion, road accidents and environmental pollution caused by vehicles in the modern international metropolis [2,3]. A metro project faces significant challenges of construction for the complex conditions, e.g., soft soil layers, interaction caused by the urban infrastructures, vibration issues, accidents during the construction phase, etc. Thus, the safety level and durability of metro construction are of critical importance [4]. The constant need to improve the operational safety of construction has driven the development of non-destructive evaluation techniques and in-site monitoring methods aimed at the ongoing assessment of the metro rail corridors and near-surface soil layers [5,6].

The most important step of the underground construction in geotechnical evaluation is in situ identification of soil properties. The soil is a kind of multiphase, particulate medium with complex mechanical profiles. The analysis of wave propagation is an essential tool for the identification of the elastic and dissipative characteristics of soil layers in dynamic conditions. The mechanical properties of the soil are obtained by analyzing dynamic problems as important parameters for site response evaluation, vibration control, earthquake engineering, etc. [3,7]. The compressional and shear waves propagating in the subsoil has been widely used for the evaluation for various purposes with a great number of techniques: seismic reflection and refraction technique, downhole testing, seismic tomography, surface wave detection technique and so on [8].

In recent decades, as a non-destructive in situ evaluation technique, the surface wave method has been attracting increasing attention from both geophysicists and geotechnical engineers. In the 1960s, after being introduced by seismologists, Jone [9] developed a prototype for the surface wave measurement system. Multiple stations were employed to estimate the characterizations of the subsurface of the Earth. Jongmans and Demanet [10] studied the surface wave technique for the evaluation of the dynamic parameters of the subsoil. Foti et al. [8] discussed the applications of various wave modes (Rayleigh wave, Stoneley wave, and Love wave) for real engineering problems.

Because of the complex mechanical properties of the subsoil, the shallow shear wave velocity is a function witnessed from the characteristics of the soil layers and the frequency of ground motion. To obtain the wave velocity map, an on-site surface wave method has been proposed to involve multiple geophones to acquire the ground motion signals for the benefit of non-intrusive testing, effectiveness, and reliable data. The multiple channel data will be analyzed by interpretation algorithms (such as, the two-receiver spectral analysis of surface wave (SASW) method [11] and the multichannel analysis of surface waves (MASW) method [12,13], etc.) to achieve an experimental multimode dispersive image.

In spite of the demonstration of accurate, high-resolution dispersion image, the application of a large number of geophones becomes a bottleneck to widespread use of the technique from the laboratory to large scalar engineering applications. Thus, in order to reduce the redundant geophones used in conventional measurements, this paper contributes to this objective by exploring a novel in-site technique based on compressive sensing (CS) to reduce the number of required geophone and measurements for the estimation of surface wave velocity of the subsoil. CS has been applied to recover the multimodal and dispersive properties of Lamb wave from observed data in the laboratory for analysis, reconstruction, and prediction of guided waves [14,15,16,17,18,19,20]. Jiang et al. applied CS to ultrasonic computerized tomography to reduce observation times for estimation of a steel tube slab structure [21,22]. This paper is organized as follows: Section 2 deals with the CS theoretical framework, wave field imaging, and the proposed fast analysis of the wave speed (FAWS) algorithm. Three actual experimental applications are implemented and analyzed in Section 3. Finally, conclusions and further developments are summarized in Section 4.

2. Materials and Methods

This section describes the process for constructing maps of local phase velocity estimates. The first part will introduce the compressive sensing methodology to apply a novel measurements strategy for evaluation. The second part will explain how the lower measurements approach to transfer the raw geophone data into local wave velocity estimates. The processing is performed with the optimized picking of the frequency-wavenumber spectral maxima by the ℓ₁-norm minimization algorithm, which is then transformed to the dispersion curve of surface wave velocity.

2.1. Compressive Sensing Framework

Compressive sampling is a new theory of information acquisition proposed by Donoho [23], Candès [24] and Tao et al. [25,26]. The fundamental mathematical theory of compressive sensing (CS) states that if a real-valued, finite-length, one-dimensional, discrete-time signal $x^{\prime }\in {ℝ}^N$ is S-sparse, i.e., only S ($S\ll N$) components of $x^{\prime }$ are non-zero, this signal can be exactly recovered from far fewer randomly chosen samples or linear measurements with an overwhelming probability. In addition, a compressible signal can also be treated as sparse if most of the transform coefficients are zeros or near zero in transform domains (e.g., Fourier, wavelet bases, spatial domain, etc.). CS is a theoretical framework that stands upon two pillars: sparsity and incoherence. Sparsity indicates that a signal’s information content can be represented in terms of a proper basis. The proper basis which transfers a compressible signal into sparse is called a “sparsifying basis”, $\Psi \in {ℝ}^{N\times N}$, which can be expressed as:

$$x=\Psi \mathit{\alpha }$$

(1)

where $x$ is the compressible signal, $\mathit{\alpha }$ is the transform coefficients vector with $S$ nonzero elements ($S$-sparse) and $\Psi$ is the sparsifying basis in a specific domain.

The general under-sampling measurements can be expressed as:

$$y=\Phi \Psi \mathit{\alpha }+\mathit{n}=\Theta \mathit{\alpha }+\mathit{n}$$

(2)

where $y\in {ℝ}^m$ is the measurement vector, $\Phi \in {ℝ}^{m\times N}$ is the under-sampling linear measurement matrix (m < N), $\mathit{n}$ represents the noise during the measuring process, and $\Theta =\Phi \Psi$ is defined as the final transfer matrix. The under-sampling measurement is the number of linear measurements $y$ smaller than the number of unknown variables in compressible signal $x$. It is an ill-posed problem to estimate the unknown variables in $x$ from few measurements $y$ [27,28]. As we mentioned above, because the unknown signal $x$ is compressible and sparse, the sparse vector $\mathit{\alpha }$ which contains information of signal $x$ could be reconstructed by solving an ${\ell }_1$-norm optimization problem.

In order to reconstruct the unknown compressible signal $x$, the second principle: incoherence is introduced due to its robustness in the presence of measurement noise [29]. Thus, the transfer matrix $\Theta$ is required to satisfy the following restricted isometry property (RIP) with an isometric constant $\delta$ smaller than unity [24]:

$$1-\delta \le \frac{\Vert \Theta v{\Vert }_2}{\Vert v{\Vert }_2}\le 1+\delta ,\delta >0$$

(3)

where $\delta$ is the RIP constant which is defined as the smallest value that meets the requirement of Equation (3), $v\in {ℝ}^N$ is all S-sparse vector. RIP constant $\delta$ is a parameter representing the character of the nearly orthogonal matrices operating on sparse vectors. $\delta \approx 0$ indicates that the transfer matrix is a nearly orthonormal matrix, while $\delta \approx 1$ indicates that the vectors in the transfer matrix $\Theta$ are redundant. Random selected matrices have been proven to meet the requirement of RIP [24], which guarantees that the CS problem can be reconstructed with an overwhelming probability. In addition, the sampling number should be larger than the required number $m>\mu ·S·\log \left(N/S\right)$ where $\mu$ is a constant (generally, $\mu =4$) dependent on the basis. Thus, as a random matrix (e.g., random Gaussian, Bernoulli and the partial Fourier matrix), it has the benefits of accurate reconstruction, good compatible with data, easy to apply, etc.

Considering the noise, the sparsest solution is solved by the $l_1$-norm optimization algorithm using a toolbox embedded in MATLAB (developed by M. Grant and S. P. Boyd, available at: http://cvxr.com/cvx/download/):

$$\widehat{\mathit{\alpha }}=\arg \min (\Vert \Theta \alpha -\mathit{y}{\Vert }_2+\lambda \Vert \mathit{\alpha }{\Vert }_1),\mathrm{where}\Theta =\Phi \Psi$$

(4)

where, $\widehat{\mathit{\alpha }}$ is the reconstructed transform coefficients vector, $\lambda$ is the Lagrangian multiplier, $\Vert ·{\Vert }_1$ is the ${\ell }_1$-norm and $\Vert ·{\Vert }_2$ is the ${\ell }_2$-norm.

2.2. Wave Field Imaging

The speed of a Rayleigh wave propagating in the elastic Earth’s surface is dispersive which means the velocity of the Rayleigh wave depends on its frequency. In addition, multiple modes of Rayleigh wave are involved at any frequency. Thus, it is difficult to analyze the real-world dataset with the interference of a multi-mode Rayleigh wave to obtain an accurate, high-resolution dispersion image. A geophone array containing hundreds of geophones was used to pursue the wave speed and to distinguish the involved modes of Rayleigh wave for technicians and engineers (see in Figure 1).

In order to rapidly estimate the velocity of a surface wave from multiple channel data, a CS platform is proposed to improve the conventional geophone array assignment and analysis method. A randomly allocated geophone array is assigned in a line instead of the traditional uniformly spacing geophone linear array in this study. The signal measured by each geophone in the frequency domain is a sum of frequency dispersive modes of Rayleigh waves:

$$y\left(r,\omega \right)={\displaystyle \sum }_{j=1}^N\frac{1}{\sqrt{r}}{G}_j\left(\omega \right)e^{-i{k}_j\left(\omega \right)r}$$

(5)

where $G_j\left(\omega \right)$ is the complex-valued amplitude of each wave mode, $k_j\left(\omega \right)$ represents the real-valued frequency dependent wavenumber of each mode, $\omega$ represents the angular frequency, and $r$ is the distance from the source to the received sensor.

Based on the character of dispersion, for each mode j = 1, 2, 3,$\cdots$, N, the complex-valued amplitude $G_j\left(\omega \right)$ and wavenumber $k_j\left(\omega \right)$ in Equation (5) are continuous functions of frequency $\omega$ and wavenumber $\kappa$ in the Rayleigh wave equation:

$$\mathit{G}\left(\omega ,k\right)=\{\begin{array}{cc}{G}_j\left(\omega \right)& \mathrm{if} k_j\left(\omega \right)=\kappa \mathrm{for}\mathrm{some} j\\ 0& \mathrm{otherwise}\end{array}$$

(6)

$$k_j\left(\omega \right)=\{\begin{array}{cc}\kappa & \mathrm{for}\mathrm{some} j\mathrm{at}\mathrm{special}\mathrm{frequencies}\\ 0& \mathrm{otherwise}\end{array}$$

(7)

where $\kappa$ indicates the actual wavenumber in the Earth’s surface, which contains the roots of total S modal Rayleigh wave of the wave equation.

Considering the randomly allocated geophone array of m geophones, the signal received by the geophone array is the superposition of all modal components of the Rayleigh waves:

$$\mathit{y}\left(\omega \right)=\mathit{A}·\mathit{\alpha }\left(\omega \right)$$

(8)

$$\mathit{y}(\omega )=\left[\begin{array}{c}{y}_1(\omega )\\ y_2(\omega )\\ ⋮\\ y_m(\omega )\end{array}\right]=\left[\begin{array}{cccc}\frac{1}{\sqrt{r_1}}{e}^{-i{k}_1(\omega )r_1}& \frac{1}{\sqrt{r_1}}{e}^{-i{k}_2(\omega )r_1}& \cdots & \frac{1}{\sqrt{r_1}}{e}^{-i{k}_N(\omega )r_1}\\ \frac{1}{\sqrt{r_2}}{e}^{-i{k}_1(\omega )r_2}& \frac{1}{\sqrt{r_2}}{e}^{-i{k}_2(\omega )r_2}& \cdots & \frac{1}{\sqrt{r_2}}{e}^{-i{k}_N(\omega )r_2}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{1}{\sqrt{r_m}}{e}^{-i{k}_1(\omega )r_m}& \frac{1}{\sqrt{r_m}}{e}^{-i{k}_2(\omega )r_m}& \cdots & \frac{1}{\sqrt{r_m}}{e}^{-i{k}_N(\omega )r_m}\end{array}\right]\left[\begin{array}{c}{G}_1(\omega )\\ G_2(\omega )\\ ⋮\\ G_N(\omega )\end{array}\right]$$

(9)

where $\mathit{A}$ is a $m\times N$ transfer matrix. In addition, the condition for the RIP is satisfied due to the random allocation of the geophone array. $\mathit{\alpha }\left(\omega \right)$ represents the solution to a discrete inverse problem and a discretized approximation of $\mathit{G}\left(\omega \right)$. For most frequency $\omega$ and wavenumber $k\left(\omega \right)$, $\mathit{G}\left(\omega ,k\right)$ is zero based on the theoretical dispersion image and wave equation. Thus, its discretization $\mathit{\alpha }\left(\omega \right)$ is sparse as well. These imply that the acquired signal $\mathit{y}\left(\omega \right)$ could be recovered by chasing the sparse solution for $\mathit{\alpha }\left(\omega \right)$ based on the CS theory from the underdetermined inverse problem:

$$\widehat{\mathit{\alpha }}\left(\omega \right)=\arg \min \left(\Vert \mathit{A}·\mathit{\alpha }\left(\omega \right)-\mathit{y}\left(\omega \right){\Vert }_2+\lambda \Vert \mathit{\alpha }\left(\omega \right){\Vert }_1\right)$$

(10)

It is noted that $\mathit{G}\left(\omega \right)$ is a two-dimensional function across frequency and wavenumber, thus, a filter bank is proposed to estimate the dispersion image in the various frequency bands. The benefit of the filter bank is to solve Equation (10) in each certain narrow band to focus on the sparsest representation of the wavenumber (velocity) from measurements. The dispersion curves of the wavenumber (phase velocity) are obtained from one frequency range to another.

The filter bank is established by a family of bandpass filters which are used to separate the information according to their frequency components. Second-order bandpass filters with tuned frequencies are employed to assemble the bandpass filter bank to separate the information at various interested frequencies (see in Figure 2). The central frequency and the spacing are the most important parameters for the parallel bandpass filter bank. The basic bandpass transfer function of the filter is expressed as Equation (11):

$$H_{BP}\left(s\right)=\frac{\frac{{\omega }_p}{Q_p}}{s^2+\frac{{\omega }_p}{Q_p}s+{\omega }_p^2}$$

(11)

where $s$ is the complex variable, $\omega$ is the angular frequency in rads⁻¹, that is, $\omega =2\pi f$, ${\omega }_p$ is the central angular frequency, and $Q_p$ is the the quantity Q-factor of the filter; higher $Q_p$ produces narrower bandwidth.

The 3dB bandwidth of the bandpass filter is defined as:

$${\omega }_{BW}={\omega }_{+3dB}-{\omega }_{-3dB}$$

(12)

where ${\omega }_{BW}$ is the bandwidth of the bandpass filter, ${\omega }_{+3dB}$ and ${\omega }_{-3dB}$ are the half-power frequencies, which are defined as:

$$\left|H_{BP}\left(j·{\omega }_{+3dB}\right)\right|=\frac{1}{\sqrt{2}}, {\omega }_{+3dB}={\omega }_p\sqrt{1+\frac{1}{2Q_p^2}+\frac{1}{Q_p}\sqrt{1+{\left(\frac{1}{2Q_p}\right)}^2}}$$

(13)

$$\left|H_{BP}\left(j·{\omega }_{-3dB}\right)\right|=\frac{1}{\sqrt{2}}, {\omega }_{-3dB}={\omega }_p\sqrt{1+\frac{1}{2Q_p^2}-\frac{1}{Q_p}\sqrt{1+{\left(\frac{1}{2Q_p}\right)}^2}}$$

(14)

In summary, the proposed FAWS is divided into five steps to carry out for the on-site experimental testings: (1) preparing step, determining the detection range and the measurement number; (2) arrangement step, generating random Bernoulli matrix and allocating the geophones based on the random matrix; (3) acquiring step, generating surface wave and recording the raw data of the ground-motion waveforms by the ununiformed geophone array; (4) processing step, filtering the raw data by a narrow-bandpass filter bank to decompose and extract the information of the surface wave from the recorded data at various interested frequencies; (5) inversion step, computing the wavenumber or velocity of surface wave at each frequency and assembling the inversed information as the final dispersive mapping of the surface wave of the subsoil to evaluate its mechanical properties. Figure 3 shows a summary flow diagram that describes the transformation process in detail.

3. Results and Discussion

A real-world experiment was employed to study the effectiveness of the proposed FAWS algorithm.

The experiment was carried out on a station of the metro system for the city of Shenyang, Liaoning province, China. Olympic Metro Station is a transfer station of Shenyang metro line 9 and line 2 in the downtown of Shenyang city. The property of the super shallow buried soil layer should be tested on site to make sure the load-carrying capacity of the underground construction. The transfer station is oriented in a west–east direction as shown in Figure 4a. For metro line 2, the station was built in 2016, thus, a tunnel is constructed to connect the previous station (line 2) to the new one of metro line 9. To estimate the mechanical properties of the soil layer, a linear array of 50 geophones with the spacing of 0.5 m was used to acquire the site response of ground under the impact. The testing site is flat, free of obstacles, and of one-dimensional geometry. The total length of 25 m test array ran through the west–east direction of the metro line 9 before the construction of the new metro station, as shown in Figure 4a. The in-site testing were conducted in summer with daily temperature from 20 to 30 °C.

A sledgehammer with a 4 kg weight was employed to excite the vertical ground motion. The impact source caused by the hammer blow is located at the offset distance of 0.5 m. The locations of the first and last geophones are 0.5 m and 25 m, respectively. The data were recorded by multiple geophones with 5 Hz resonant frequency and a digital geometrics seismography to save as a .dat file with a sampling frequency of 1000 Hz. The record length of the multiple channel data is 600 ms. A schematic of the experimental setup is shown in Figure 4b.

3.1. Conventional Even-Spaced Array Results

As a comparison, the data were sampled by regular geophones and stored without any filtering or filter bank. The raw data acquired with the impulsive sledgehammer source and the 50 vertical geophones are demonstrated in Figure 5a. In addition, the time-domain signal recorded by the number 10 geophone with a 5-m distance from the impact source was transformed to a time-frequency spectrogram using Gabor wavelet transform, as shown in Figure 5b. It is noted that the energy in the measured wave field concentrates in the frequency band of 40–80 Hz.

To address the dispersion curve of phase velocity of the near-surface earth, waveform analysis is processing using frequency-wavenumber transform (f-k transform) by two-dimensional fast Fourier transform (2D-FFT). This technique is extensively used for many near-surface applications as a full-waveform inversion approach. All the time-domain signals recorded by the even-spaced array were transformed to a wavenumber-frequency spectrogram and a dispersion image using 2D-FFT as a function of two variables: time and offset.

The reference wavenumber-frequency spectrogram of the conventional technique is obtained by applying 2D FFT on the even-space signals transferring them from the time domain into the frequency (f-) domain, from spatial domain (offset) into wavenumber (k-) domain. The transform of the multi-sensed surface wave data is aimed at identifying the wavenumber with the energy propagating at each frequency. The f-k transform provides an obvious image of the multiple modal propagations, as shown in Figure 6. The fundamental mode and the first higher mode could be clearly separated at the different dominant ranges. The resulting image represents the energy density as functions of the wavenumber: one wave mode is focused on 50 Hz with a wavenumber of 1 ${\mathrm{m}}^{-1}$, and the other one is around 60 Hz with a wavenumber of 2 ${\mathrm{m}}^{-1}$(Figure 6a). Similarly, the phase velocity-frequency image of the soil layers yields two main dispersion curves: one mode ranges between 150 to 200 m/s with a central frequency of 60 Hz, the other curve is about 400 m/s with 50 Hz frequency, as shown in Figure 6b.

3.2. Results of the Proposed Fast Analysis of Wave Speed (FAWS) Method

To validate the capability of the proposed FAWS method in the CS framework, geophones were randomly chosen to acquire the data based on the Bernoulli matrix. Figure 7a shows the random distribution of the active geophones. The total number of the geophones in the non-even-spaced array is 20 instead of 50 geophones in the whole 25-m measurement length. The locations of the activated geophones are randomly selected based on Bernoulli matrix as mentioned in the previous section. The black block represents a geophone at the corresponding location; by contrast, the white block means no geophone is arranged there. The bandwidth and the filter spacing of the bandpass filter bank is 0.5 Hz and 1 Hz, respectively. Thus, there are 98 filters in the desired range of 3 to 100 Hz.

Considering the impulsive wave propagation admits only a few modal wave number at each frequency in the dispersion image, the sparse f-k spectrogram is reconstructed based on the ℓ₁-norm minimization algorithm of CS theory (Equation (10)). The reconstruction results are shown in Figure 7. The two modes can be obviously separated in the wavenumber-frequency spectrum. The reconstructed results are demonstrated in an image whose amplitudes is mapped by a white-red colortable: the higher magnitude, the darker red color. It can be seen the wavenumber-frequency spectrum have two peaks: one peak is observed at 50 Hz with the wavenumber of 1 ${\mathrm{m}}^{-1}$, the other is in the range of 50 Hz to 60 Hz with the wavenumber around 2 ${\mathrm{m}}^{-1}$, which are in good agreement with the conventional 2D-FFT algorithm.

Meanwhile, as shown in Figure 7c, the two dispersion curves of the phase velocity of surface wave are accurately identical in the whole range, which means the proposed method transform can be used to identify the properties of the two modes. Thus, in the application of Shenyang metro line 9, a total number of 20 geophones are random selected from 50 geophone space array for measuring the seismic wave to evaluate the mechanical properties of the soil. This proposed technique reduces the acquisition time to achieve high speed in characterization of the shallow subsurface by removing 60% of the Nyquist sampling grid.

3.3. In-Site Downhole Testing Results

In addition, the conventional in-site downhole testing method was also employed to directly obtain the measurements of compression (P-) and shear (S-) wave velocity as a comparison. Figure 8 shows the prepared borehole with a diameter of 127 mm for downhole testing according to the ASTM standard test method [30]. Four boreholes were drilled 15 m on the ground to put the downhole receivers. A 14-bit ADC measurement equipment (XG-1 type 3-axis receiver manufactured by Langfang Dadi corporation, Ltd. Langfang, China) was employed to acquire the wave velocities through a sampling frequency of 1000 Hz. After the 3-axis receiver were placed at the desired test locations in the downhole, the energy source is activated to generate ground motion, the waves were recorded three times by the receivers to improve the signal-to-noise ratio (SNR). The experimental results are listed in

Table 1. Experimental in-site downhole results at Olympic station of Shenyang metro line 9.

Bore Hole Number	Soil Type	Depth (m)	Shear Wave Velocity (m/s)	Compression Wave Velocity (m/s)
ZDBC-01	Topsoil	2.9	170	445
	Silt and clay	5.3	200	529
	Sand	7.8	269	709
ZDBC-02	Topsoil	2.2	154	424
	Silt and clay	4.3	195	532
	Sand	8.5	244	634
ZDBC-03	Topsoil	3	161	431
	Silt and clay	5.2	199	523
	Sand	6.8	261	681
ZDBC-04	Topsoil	3.4	172	439
	Silt and clay	6	187	488
	Sand	5.6	257	678

. There are 3 kinds of soil layers: topsoil layer, silt and clay layer, and sand layer with various depths at different locations. The S-wave velocities of the top two layers are identified in the ranges between 150 to 200 m/s, and the P-wave velocity is in the range of 400–500 m/s, as shown in Table 1. The results of in-site observation provide sufficient evidence to the surface wave identification methods for both the proposed FAWS method and the conventional even-spaced geophone array. As we can see, the results are in good agreement with the previous estimation results obtained by the proposed method. Considering the time-consuming drilling, complicated implement, and redundant observations, the proposed FAWS technique is an effective and efficient game-changing technique with accurate estimation of the wave velocities of soil properties in actual engineering applications.

4. Conclusions

In this study, an accurate FAWS algorithm is proposed based on compressive sensing theory. Fewer measurements of a random-spaced geophone array are needed to provide the desired dispersion curve of surface wave velocity for the non-destructive evaluation of the soil layers. In-site experimental testing is deployed in Shenyang metro line 9 to validate the capability of the proposed technique. A bandpass filter bank is introduced to reconstruct the multi-modal dispersion characters of wavenumber and surface wave velocities at various frequencies. According to the sparse features of the Rayleigh-type seismic surface wave of the Earth, an $l_1$-norm optimization algorithm is employed to reconstruct the dispersion curve of the velocity and wavenumber-frequency spectrum.

As comparisons, both a conventional even-spaced geophone array and in-site downhole testing are carried out to evaluate the soils in Shenyang metro line 9 for the future construction. The dispersion curve identified by the FAWS method was picked up with a good agreement with the trend of the dispersion image of surface wave processed by the comparison technique. In addition, the experimental results indicate the proposed technique can identify both the wavenumber-frequency distribution and velocity-frequency image of the surface wave of subsoil with very good accuracy. Considering the measurements required are fewer than the traditional method, the proposed approach is a game-changing technique with the advantage of being low cost, rapid, and safe for the identification of the soil properties for both geotechnical research and engineering construction.

It should be emphasized that the FAWS technique has its own specific requirements and limitations, which might need comparisons and supplements from other in-situ testing for a complete, accurate evaluation of deeper soil properties.

In future work, the proposed technique requires additional study for the optimization of the arrangement of the geophones and the filter bank. In addition, the proposed in-situ FAWS technique will be applied to evaluate the soil properties on a larger scale for both engineering applications and geophysical research. More efficient and powerful algorithms, e.g., orthogonal matching pursuit (OMP), CoSaMP, and subspace-pursuit (SP), will be investigated in future work. Also, the characteristics of deep soil layer properties will also be studied to pursue the tomography of the subsoil.