Shallow Structure Imaging Using Multi-mode Dispersion Curves Based on Multi-window HLRT in DAS Observation

Distributed acoustic sensing (DAS) is a high-density seismic observation system emerging in recent years. It can achieve low-cost and high-density observation and is one of the important methods for imaging city shallow structure. However, the high-mode dispersion spectrum information recorded by DAS system is difficult to utilize effectively. To solve this problem, we propose to extract multi-mode dispersion curves of active and passive sources through multi-window high-resolution linear Radon transform (HLRT) method and invert multi-mode dispersion curves of active and passive sources to construct shallow surface underground S-wave velocity structure. The research results show that the multi-window HLRT transform method can effectively extract multi-mode dispersion curves from DAS records by adding time windows; the improved objective function further enhances the stability and accuracy of inversion. The development of multi-mode dispersion curves extraction and inversion technology provides a new method for low-cost and high-resolution exploration of underground structures with DAS system.


Introduction
The Rayleigh surface wave exploration method has the characteristics of being non-destructive, noninvasive and high-resolution and is widely used in city geological exploration (Xia et al., 1999;Yoon, 2011). However, the Rayleigh surface wave exploration method is more dependent on high-density observation systems, and the high-density observation system is difficult to arrange on a large scale in cities because of instrument arrangement, instrument maintenance and high economic cost (Becker et al., 2017;Dou et al., 2017).
Distributed acoustic sensing (DAS) is a highdensity seismic observation system emerging in recent years. The distance between stations can reach the order of 1 * 10 m. Its system structure is relatively simple, the maintenance cost is very low, and it has the function of real-time data transmission (Becker et al., 2017;Daley et al., 2016;Dou et al., 2017). DAS technology monitors the dynamic strain of optical fiber by measuring the phase change of Rayleigh backscattered signals in optical fiber and then demodulates the seismic wave field records (Parker et al., 2014). DAS system has been widely and effectively used in shallow structure imaging (Tribaldos et al., 2021;Zeng et al., 2017), groundwater level monitoring (Becker et al., 2017;Mateeva et al., 2017), etc.
DAS system has the characteristics of long-distance, high-density, and long-term observation, and the observed value of the DAS system is the axial dynamic strain of optical fiber, which can be equivalent to the R component of the seismograph by demodulating the phase of light pulse. So, effective high-mode surface wave signals are easily obtained (Foti et al., 2000;Dal Socco et al., 2010;Moro et al., 2015;Hu et al., 2020). However, in the actual processing of DAS data, high-mode signals are difficult to use effectively. The reason is that most of the energy is concentrated in the fundamental dispersive spectrum, and the peak energy of the high-mode dispersive spectrum is lower, so high-mode dispersion curves are difficult to pick up. In addition, if the adjacent dispersive spectra are close to each other, the information of the dispersive spectrum will be polluted. As a result, multi-mode dispersion curves are very prone to mode mixing, resulting in misidentification of the dispersion curve mode, which brings a great challenge to the effective utilization of highmode dispersion curves (Luo et al., 2010;Maraschini et al., 2010;Tokimatsu et al., 1992;Zhang & Chan, 2003). Therefore, in the shallow surface structure imaging with DAS, obtaining a high resolution and high signal-to-noise dispersive spectrum image is very important for extracting the dispersion curve and subsequently suppressing the multi-solution of inversion results.
For the F-K transform, Tau-p transform, phase shift and slant stacking algorithm, which typically suffers from resolution loss and aliasing as a consequence of insufficient information, producing a highresolution dispersion image is difficult (Luo et al., 2009). High-resolution linear Radon transform can extract fundamental and high-mode dispersion curves (Luo et al., 2008(Luo et al., , 2009Pan et al., 2016). Wang et al. (2019) propose the frequency-Bessel (F-J) transform method, which can stably extract high-resolution fundamental and high-mode dispersion curves from ambient noise data. However, in actual exploration, it is a great challenge to accurately pick up the dispersion curve by directly transforming the seismic records to the f-v (frequency-phase velocity) domain. In the process of high-mode dispersion curves are extracted, the duration of fundamental-mode surface waves is longer than for other waveforms, and the amplitude difference is obvious (Visser et al., 2008;Yoshizawa & Ekström, 2010). The waveforms with small energy are masked when converting the waveforms into the f-v domain through fast Fourier transform (FFT). In addition, if the adjacent dispersive spectra are close to each other, the information of the dispersive spectrum will be polluted. Luo et al. (2009) propose to improve the resolution of highmode surface waves by mode separation, but it is necessary to intercept the target dispersive spectrum in the f-v domain, transfer it to the time domain and then transfer it to the f-v domain to extract the dispersion curve, which is complicated and tedious. Li and Chen (2020) add time windows to the F-J transform method, and fundamental and higher-mode dispersion curves are extracted by multi-window.
We propose to use the multi-window high-resolution linear Radon transform (HLRT) method to extract multi-mode dispersion curves, which enables the DAS system to extract richer dispersive spectrum information, to introduce uncertainty estimation into the inversion objective function and to improve the accuracy and stability of inversion. Our method provides a new method for near-surface high-resolution exploration in cities.

Introduction Multi-window HLRT Method
The high-resolution linear Radon transform (HLRT) extraction dispersion curve is relatively mature, and this section provides a brief overview of it. The forward LRT can be written in matrix form as: and the adjoint transformation, where L ¼ e ipf 2px is the forward LRT operator, d and m represent the shot gather and Radon panel, and m adj denotes a low-resolution Radon panel using the transpose or adjoint operator L T . It is difficult to achieve high-resolution images of dispersion energy from standard RT because they suffer from the loss of resolution and aliasing arising as a consequence of incomplete information, including limited aperture and discretization (Trad et al., 2003). For that purpose, they proposed a stochastic inversion technique that is able to retrieve a sparse Radon panel.
To invert Radon operators, a sparsity constraint should be considered when finding the model m that best fits the data while minimizing the number of model-space parameters. The model m can be defined by solving the following Eq. (3) (Trad et al., 2002(Trad et al., , 2003, where I denotes the identity matrix, W d is a matrix of data weights, and W m is a matrix of model weights that plays an extremely important role in the design of high-resolution Radon operators. k maintains balance between data misfit and model constraints. Detailed strategies for the solution of Eq. (3) can be found in many resources Ji, 2006;Luo et al., 2009;Sacchi & Ulrych, 1995;Schultz & Gu, 2013;Trad et al., 2002Trad et al., , 2003.
DAS system easily records multi-mode surface waves in city near-surface imaging applications, but there are some challenges in extracting high-mode dispersion curves. We give a velocity interval ([V1, V2]) to calculate time windows ([T1 = r/V2, T2 = r/ V1]) for each epicentral distance (r) to intercept waveforms. The intercepted waveforms are transformed into f-v domain by HLRT, and then dispersion curves are extracted. We call this method multi-window HLRT. The window function of the intercepted waveform is as follows: Fs is the sampling rate; out is the window function. To avoid the time window being too short, the time interval ([T1, T2]) should be extended appropriately. Generally, time window selection is flexible. As long as the fundamental waveform is separated, multi-window HLRT can be effectively performed.

Inversion Objective Function
In the process of multi-mode dispersion curve inversion, mode skip and mode mixing are often encountered. However, the DAS system very easily records high-mode dispersion curves. If the mode misidentification can be avoided, the advantages of DAS system in imaging near-surface can be further enhanced. Therefore, we use the following inversion objective function: where N is the number of data to be picked up, PV(f i ) is the phase velocity of the i-th frequency (f i ) to be picked up, and C(f i ) is the phase velocity (calculated multi-mode phase velocity) calculated by f i . r(f i ) is the uncertainty estimate of picking up PV(f i ). In the process of inversion, all data points are isolated from each other, so it is not necessary to identify the number of modes corresponding to the picked up data. Only the minimum square difference of C(f i ) of one mode with PV(f i ) is required. The introduction of r(f i ) increases the anti-error ability of inversion and also reflects that r(f i ) of extracting the dispersion curve is too high, which may cause the inversion results to be distorted. Therefore, it is necessary to increase the resolution of the dispersion spectrum and reduce the uncertainty of pick-up dispersion curves for the accuracy and reliability of inversion results.

Multi-window HLRT Method Test
A four-layer horizontal layered stratigraphic model (Table 1) Vol. 180, (2023) Shallow Structure Imaging Using Multi-mode Dispersion Curves 865 simulation. One hundred geophones received the signal. The channel spacing was set to 2 m, the minimum offset was set to 20 m, the sampling rate was set to 2000 Hz, the recording time was set to 2 s, and the simulation records were the radial component of the seismic signal.
Taking the 50th channel of synthetic seismogram as an example (Fig. 1), it can be seen that the fundamental and high-mode waveforms have obvious arrival boundaries. The duration of the fundamental waveform is longer, and the amplitudes are quite different. According to the waveform record in Fig. 1 The result of multi-window HLRT to extract the dispersion curves is shown in Fig. 2. When the dispersion spectrum is directly transformed without time windows, the mutual interference between the fundamental and higher-mode surface wave and the poor resolution of the dispersion spectrum between 5 and 20 Hz make it difficult to accurately extract the dispersion curve. The seismic waves intercepted by the Win1 are transformed to the f-v domain, because the high-mode surface waves are separated, and the 5-20 Hz fundamental dispersion information can be displayed. The fundamental dispersion spectrum is consistent with the theoretical phase velocity curve, which also shows the correctness of multi-window HLRT. The high-mode dispersion spectrum intercepted by Win2 shows that its resolution is significantly improved, and the pickup accuracy of the high-mode dispersion curve is improved because the fundamental and high-mode surface waves are separated. We add time windows to the frequency Bessel transform and transform the seismic records to the fv domain, as shown in Fig. 3, which also proves the effectiveness of the added time windows. Notably, in the formula of frequency Bessel transform, a general integral from 0 to positive infinity needs to be realized by numerical integration of observed data (Hu et al., 2020;Li & Chen, 2020;Wang et al., 2019). Due to our offset distance, the phase velocity energy peak superposition of the dispersion spectrum converges to the extreme value poorly. In the actual exploration of DAS system, the performance of the frequency Bessel transform will also be reduced due to the offset distance.

Objective Function Test
We use Bayesian trans-dimensional inversion (Bodin et al., 2012;Killingbeck et al., 2018) to test the objective function. The stratigraphic model is shown in Table 1. The prior distribution information on v S (defining upper and lower bounds) of each layer is shown in Table 2. The maximum number of floating nuclei was set to 15, the number of confined nuclei was set to 3, the burn-in chain length was set to 10,000 and 1,000,000 iterations, r change -v S (standard deviation of Gaussian proposal on change v S value) was set to 20 m/s, r move -depth (standard deviation of Gaussian proposal on depth move) was set to 1 m, and r birth/death -v S (standard deviation of Gaussian proposal on birth and death) was set to 400 m/s. The inversion results of the objective function with uncertainty estimation r(f i ) are shown in Fig. 4a-d. When the picked data are not added with noise, the posterior probability density distribution of v S is very concentrated, and the fitting degree between the best fitting data and the picked data is also very high, which shows the correctness of the objective function. When 10% white noise is added to the picked data, the posterior probability density distribution of v S is concentrated near the real model, and the best fitting data are close to the picked data with 10% white noise. The inversion results of the objective function without uncertainty estimation r(f i ) are shown in Fig. 4e, f. Inversion results pursue over-fitting, the posterior probability density distribution of v S is very different from the actual model, and many abnormal thin layers appear. The inversion results show that the introduction of uncertainty estimation r(f i ) ensures the anti-interference capability of the algorithm and improves the stability of inversion of actual data.

Multi-window HLRT Extract DAS Passive Source Dispersion Curve
The experimental site is located at the North Square of Geological Palace of Jilin University, and about 400 m optical fiber was arranged by shallow burying. Our data were from 53th to 78th channel recordings (Fig. 5). In the data collection process, the channel interval was set to 5 m, the sampling rate was  set to 2000 Hz, and the total time for continuous recording was 3 h (12:00-15:00). The 1-min continuous recordings were sampled down to 100 Hz; 1-20 Hz band-pass filtered recordings were selected as normalization factor of time domain normalization. Spectrum whitening was performed in the frequency band, and the NCFs were calculated. Finally, the 1-min NCFs were stacked to 3 h. Figure 6 shows the NCFs of the 80th channel as the virtual source and 78th-60th channels as the receiving channels. After 3-h NCFs had been stacked, the surface wave signals were more significant and had a higher signal-to-noise ratio. The obtained NCFs are converted into dispersion spectra as shown in Fig. 7.
The DAS records are directly transformed to the f-v domain without time windows, the dispersion spectrum has poor continuity and low resolution, and it is difficult to accurately pick up the dispersion curve (Fig. 7a). We set the large time window Win(150-300 m/s) to intercept DAS records and enhance surface waves. Because noise data are eliminated, the dispersion spectrum is continuous (Fig. 7b). With Win(150-180 m/s), the fundamental surface waves are separated well, and the fundamental dispersion curve at higher frequencies can be extracted (Fig. 7c). With Win(180-300 m/s), because the actual data are more complex, setting the time windows can only weaken the energy of the fundamental waveforms in the f-v domain, which is difficult to eliminate. Therefore, the information of the fundamental dispersive spectrum is included in Fig. 7d. Fortunately, it is consistent with the information on the dispersive spectrum in Fig. 7c and is not distorted. Figure 7 shows that the actual dispersion spectrum is complex, and the mode identification of the dispersion curve faces great challenges. Therefore, it is necessary to use our improved objective function for inversion.

Multi-window HLRT Extract DAS Active Source Dispersion Curve
In Fig. 5, when the excitation point is in the same straight line as optical fiber, optical fiber is equivalent to the R component seismograph that can record Rayleigh surface wave signals. To remove the interference of other signals to extract dispersion curves, a 5-40-Hz band-pass filter was applied to the raw data of the DAS system. The velocity range of surface waves is between 100 and 500 m/s (Fig. 8b). Then, a F-K filter was applied in the range of 100-600 m/s. Figure 8c shows that the noise recorded by DAS system is suppressed well by band-pass filtering and F-K filtering. The final DAS records were obtained by picking up the arrival time of seismic (Fig. 8d).
The DAS records are transformed to the dispersion spectrum, and the active source can obtain the information of higher frequency and higher mode than the passive source (Fig. 9). When there are no time windows, the resolution of the dispersion spectrum is low (Fig. 9a). With Win (110-200 m/s), although the surface wave signal is enhanced, highmode dispersion spectra of the active source have stronger energy and interfere with each other, so the continuity of the dispersion spectrum is poor (Fig. 9b). With Win (110-150 m/s), the fundamental and part of the high-mode spectra are clearly displayed (Fig. 9c). With Win (150-200 m/s), the resolution of the high-mode dispersion spectrum becomes higher, which is complementary to the information of dispersion spectrum intercepted by Win (110-150 m/s). By comprehensively comparing Fig. 9, the dispersion curve can be accurately picked up.

Trans-dimensional Inversion of Active and Passive Source
In theory, the passive source can provide lower frequency band information and increase exploration depth, while active source can provide higher frequency band and higher mode dispersion information and increase exploration resolution (Luo et al., 2007;Xia et al., 2003).
We invert the dispersion curves extracted from the above active and passive sources. The prior distributions on v S are shown in Table 3. The maximum number of floating nuclei was set at 15, the number of confined nuclei was set at 3, the burnin chain length was set to 10,000, with 1 million iterations, r change -v S was set to 20 m/s, r move -depth was set to 1 m, and r birth/death -v S was set to 300 m/s. The inversion results are shown in Fig. 10. The best fitting data fit well with the picked data, and are near 10-15 Hz for passive source inversion. Our objective function automatically determines the mode of dispersion curves. In addition, the passive source can detect more deeply, and active source shows a high velocity layer near 10 m. Combined with the analysis of active and passive source results, within the 95% credible intervals, the posterior probability density distribution of v S is concentrated and has high credibility at depths of about 1-15 m. At depths of about 20-80 m, the posteriori probability density distribution of v S in active source is very ambiguous, and it is difficult to reflect effective information.
Because the dispersion spectrum becomes wider in low frequency band, the uncertainty of the pick-up dispersion curve increases, and the sensitivity of v S to stratum changes decreases with the increase of depth.  Therefore, the interval of posterior probability density distribution of v S becomes wider in passive source, but it can be seen that the posteriori probability density of the v S of 25-50 m is concentrated at about 380 m/s; under 50 m, v S is concentrated at about 650 m/s. Each time the excitation point was moved to the west by 10 m, the composite dispersion curves were extracted as described above. The obtained dispersion curves were inverted, and the v S -h model was taken as the velocity structure below the midpoint of receiving channels; the S-wave velocity structure of experimental field site shown in Fig. 11. The underground structure is divided into four layers. The first layer, with S-wave velocity of about 200 m/s and thickness of about 5 m, is a loose soil layer. The second layer, with the S-wave velocity of about 200-400 m/s and depth of about 5-20 m, has a mixture of soil and gravel. The third layer, with S-wave velocity of 450-600 m/s and depth of 20-50 m, is presumed to be a weathered layer. The fourth layer, below 50 m, with the S-wave velocity [ 650 m/s, is presumed to be bedrock.

Discussion
In Fig. 2c, interference occurs on adjacent highmode dispersion curves near 25-30 Hz, which reduces the confidence of dispersion information. Through further analysis in Fig. 1, there are obvious arrival boundaries and amplitude differences between the fundamental and high-mode waveforms, but there is no such difference between the high-mode waveforms. It is difficult to separate the high-mode waveforms through time windows because of the differences between waveforms. Therefore, we use a sliding time window strategy to extract high-mode dispersion curves (Fig. 12).
Slide Win2 V1 to 350 m/s and 450 m/s to intercept the waveforms, respectively; the waveforms Inversion of stratum profile along optical fiber: a inversion of stratum profile for active source, b inversion of stratum profile for passive source Vol. 180, (2023) Shallow Structure Imaging Using Multi-mode Dispersion Curves 875 intercepted by Win (350-500 m/s) are converted into f-v domain. The resolution of the high-mode dispersion spectrum at 30-40 Hz is significantly improved, which is consistent with the theoretical dispersion curve. Win (450-500 m/s) dispersion spectrum is distorted because too many high-mode waveforms are cut off. Therefore, in the actual data processing, the setup of sliding windows depends on experience and many attempts. We suggest that the dispersion information with low certainty should be ignored, and only the high certainty be extracted by comparing the dispersion spectrum of different time windows.

Conclusions
We extract the high-resolution multi-mode dispersion curves from active and passive source signals through the multi-window HLRT method and invert the dispersion curves with improved objective function. Our methods are also used for data processing of actual DAS systems. We reach the following conclusions.
DAS observation system has the characteristics of long distance and high density, and it is easy to collect high-mode surface wave signals. However, in the actual data processing, it is difficult to effectively use the high-mode signals. The multi-window HLRT method can effectively extract the high-resolution Rayleigh surface spectrum by adding time windows to make full use of the seismic signals collected by DAS system.
In addition, we use the objective function independent of the mode number of dispersion curves, which largely avoids human-caused mode identification errors. The introduction of uncertainty estimation as a constraint enhances the anti-noise ability of the algorithm. These strategies further enhance the advantage of DAS system in urban nearsurface imaging.