Near-surface Site Characterization Based on Joint Iterative Analysis of First-arrival and Surface-wave Data

Near-surface site characterization is of great significance in the fields of geotechnical engineering and resource exploration. In this paper, we propose a near-surface site characterization method based on the joint iterative analysis of first-arrival and surface-wave data (JIAFS). The proposed method combines the advantages of first-arrival traveltime tomography (FATT) and multichannel analysis of surface waves (MASW). First, the 1D S-wave velocity (vS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text{S}}}$$\end{document}) models obtained by MASW are interpolated to construct the pseudo-2D vS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text{S}}}$$\end{document} model. According to the available geological survey information and borehole data, the initial Poisson’s ratio (σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document}) model is estimated. Based on the estimated σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document} model, the pseudo-2D vS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text{S}}}$$\end{document} model is converted to a referenced P-wave velocity (vP\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text{P}}}$$\end{document}) model which is utilized to constrain the progress of FATT. This helps FATT overcome the inherent defect that it cannot effectively identify velocity-inversion interfaces and low-velocity zones. On the other hand, the vP\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text{P}}}$$\end{document} model obtained by FATT can provide a favorable priori information to improve the reliability of the results of MASW. Then, the vP\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text{P}}}$$\end{document} and vS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text{S}}}$$\end{document} models obtained by constrained FATT and MASW are used to update the σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document} model. In addition, the vP\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text{P}}}$$\end{document} and vS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text{S}}}$$\end{document} models are also used as initial models in the next iterative analysis. Finally, through the iteration of this process, the two inversion methods can make use of their own advantages to improve each other, so we can establish accurate near-surface vP\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text{P}}}$$\end{document}, vS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text{S}}}$$\end{document} and σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document} models under complex geological conditions. A velocity model including low-velocity zone is established for synthetic model test to analyze and verify the advantage of JIAFS. The vP\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text{P}}}$$\end{document}, vS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text{S}}}$$\end{document} and σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document} models obtained by JIAFS can accurately identify the low-velocity zone and match the true models well. In addition, the proposed method is applied to the field seismic data acquired for oil and gas exploration in Northwest China. Compared with the results of individual inversions and borehole data, JIAFS can establish more reliable 3D vP\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text{P}}}$$\end{document}, vS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text{S}}}$$\end{document} and σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document} models by interpolating the 2D inversion results, which reveals further details and enhances the geological interpretation significantly.


Introduction
Accurate near-surface site characterization is essential for geotechnical engineering, oil and gas exploration, environmental study, mineral resources survey, soil mechanics, natural disaster prevention and urban underground space construction (Lakke et al. 2008;Zhu et al. 2008;Socco et al. 2010;Pegah and Liu 2016). The near-surface geophysical exploration methods mainly include gravity, electromagnetic and seismic methods. Due to the advantages of high accuracy, high resolution and large penetration depth, the seismic method has been widely used in near-surface site characterization.
Differences in the deposition time, structural morphology, lithology, degree of compaction, water content and surface undulations of the near-surface stratum cause near-surface anomalies and affect the propagation of seismic waves. Theoretically, compressional-and shear-wave velocities ( v P and v S ) are related to the mechanical parameters of rock (Foti et al. 2014). Compared with the soil skeleton, the compressibility of the pore fluid has a greater impact on v P . The pore fluid dominates the overall response of the saturated media; therefore, v P is a finite value in saturated soil. Since the pore fluid has no shearing resistance, the influence of the pore fluid on v S can be ignored (Foti et al. 2002). This avoids the water-masking effect in saturated media. Poisson's ratio ( ) is another important mechanical parameter of near-surface materials, which can be calculated from the ratio of v P and v S Li et al. 2022). In geotechnical engineering, the parameters ( v P , v S and ) are widely applied for risk assessment under specific site conditions, geological disaster prevention, assessing lithology and saturation of near-surface sediments, designing structural foundations and monitoring infrastructures (Ivanov et al. 2000;Park 2013;Sun et al. 2020).
The near surface is usually covered by loose sediments, which are called the weathered layers . For reflection seismic exploration, it is well known that static correction is a vital procedure in seismic data processing (Cox et al. 1999;Li et al. 2009). Static correction is the time-shifting processing of seismic traces to compensate for the effects of variations in elevation, weathering thickness, weathering velocity or reference to a datum (Su et al. 2021). So that all measurements are corrected on a flat plane without the influence of weathering or low-velocity material present (Sherriff and Geldart 1995). The calculation of statics correction requires several parameters relating to the near surface, such as elevation, thickness and velocity (Steeples et al. 1990). Therefore, obtaining accurate near-surface velocity models is the prerequisite and basic guarantee for static correction (Zhu et al. 1992;Boiero et al. 2013;Brodic et al. 2020). Besides, multicomponent converted-wave seismic exploration is a main method to study seismic anisotropy and has become a hot topic in the field of geophysics in recent years (Mari 1984;Pan et al. 2021).
For multicomponent converted-wave seismic data processing, accurate near-surface v P and v S models are both required. On the other hand, migration imaging is highly dependent on the accuracy of the shallow velocity structures. The shallow velocity error can be transmitted to the underlying stratum, which will affect the imaging accuracy of the deep geological structures. Therefore, obtaining an accurate near-surface velocity model is of great significance for reflection seismic exploration in complex survey areas.
At present, the main methods of near-surface seismic exploration include full-waveform inversion (FWI), first-arrival traveltime tomography (FATT) and multichannel analysis of surface waves (MASW). However, each method has its own advantages and defects. FWI is a cutting-edge technology to utilize the full content of seismic signal Zhang et al. 2021). By minimizing the difference between observed and modeled seismic waveforms, it can provide high-resolution models of subsurface velocity, absorption and reflectivity for seismic imaging and interpretation (Virieux and Operto 2009;Xing and Mazzotti 2019). Although FWI can get accurate near-surface parameters, it heavily relies on the reasonable initial models. Otherwise, it is a very time-consuming process because of the huge amount of computation (Tran and Hiltunen 2012;Pan et al. 2019).
The first-arrival waves are a type of seismic waves that start from the seismic source and first arrive at the receiver through the underground medium. As shown in Fig. 1, the firstarrival waves mainly include direct and refraction waves in homogeneous layered medium. They propagate near the surface and carry a large amount of near-surface information. FATT establishes the subsurface velocity structures through the inversion of the firstarrival traveltimes (Zhang and Toksoz 1998;Whiteley et al. 2020). The sensitivity matrix is computed by projecting the traveltime residuals along the raypaths. The tomographic system is usually constrained by smooth regularization and solved with a conjugate gradient algorithm. Because of the significant advantages of clear phases, reliable identification and universality, FATT has been widely used in various fields, such as earth structure studying, shallow-reservoir characterization, geophysical prospecting and geoengineering (Aki and Lee 1976;Nolet 1987;Hole 1992;Zhu et al. 2006).
From the geophysical point of view, the cementation and low compaction of shallow stratums may produce extremely complex geological structures. FATT suffers from inherent limitations associated with velocity-inversion interfaces, low-velocity zones and strong velocity contrasts (Ivanov et al. 2005). Because most raypaths of first-arrival waves Fig. 1 Schematic diagram of direct, refraction and Rayleigh waves in homogeneous layered medium 1 3 cannot pass through the velocity-inversion interfaces (the overlying velocity is higher than the velocity of underlying stratum), velocity-inversion interfaces might corrupt the tomographic solutions, and even result in inaccurate estimation of velocity and interface depth (Zhu et al. 2001;Liu et al. 2010). When a high-velocity bedrock is covered by the low-velocity layers with strong velocity contrasts, the seismic waves refracted by the highvelocity bedrock may reach the receivers earlier than those refracted by the low-velocity overburden. Therefore, strong velocity contrasts may lead to incorrect estimation of the depths of bedrock and low-velocity overburden. On the other hand, although the v S structure can be established by S-wave first-arrival traveltime tomography, it is limited by the specific equipment and the inadequate quality of S-wave seismic data. Besides, in multicomponent seismic data, S-wave first arrivals are sometimes mixed with the P-wave reflection signals and become difficult to identify and pick. And most existing automatic picking procedures are inapplicable for S-wave first arrivals. Therefore, accurate identification and picking of S-wave first-arrival traveltimes are another challenge for applying FATT to establish v S model. MASW based on the dispersion characteristic of surface waves has been rapidly developed and received more and more attention in recent years due to its noninvasive, costeffective and efficient features (Di Fiore et al. 2016;Barone et al. 2021;Rahimi et al. 2021). Rayleigh wave is a type of surface waves that propagates near the free surface with the characteristics of low frequency, low velocity and strong amplitude (Rayleigh 1885;Liu et al. 2020). It is formed by the superposition and interference of the multiple reflections of P and SV waves. As shown in Fig. 1, the particle motion of Rayleigh wave is elliptical and opposite to its propagation direction (Mi et al. 2020). In vertically heterogeneous formations, Rayleigh wave appears geometric dispersion phenomenon (phase velocity dependency from frequency) during propagation. Linear arrays are usually placed at the surface to record the surface-wave data generated by the sources. To obtain the distribution of surface-wave dispersion energy, the surface-wave data are then transformed to frequency-phase velocity domain by dispersion imaging methods Su et al. 2022). At present, the widely used dispersion imaging methods are τ-p transform method (McMechan and Yedlin 1981), F-K transform method (Yilmaz 1987), phase-shifting method (Park et al. 1998) and high-resolution linear Radon transform method (Luo et al. 2008).
Traditionally, multimode surface-wave dispersion curves are determined by extracting the trend of dispersion energy peak. In recent years, the automatic picking of dispersion curves has attracted the attention of researchers, which can significantly improve the processing efficiency and accuracy of surface-wave data (Rovetta et al. 2020;Dai et al. 2021;Song et al. 2021). At present, local linearized inversion method and global search inversion method are usually adopted for dispersion curve inversion (Wathelet et al. 2004;Cercato 2009). According to the dispersion curve inversion, the one-dimensional (1D) v S structure with high vertical resolution can be obtained . Multiple 1D v S structures can be constructed into two-dimensional (2D) and three-dimensional (3D) v S models through interpolation and smoothing Ikeda and Tsuji 2015).
In reflection seismic exploration, surface waves are traditionally considered as coherent noise interference and should be removed. Because surface-wave energy accounts for a large portion of seismic wave energy, removal processing results in a waste of near-surface information they carry (Socco et al. 2009). In recent years, the application of surface waves in reflected seismic data to characterize the near surface has attracted the attention of researchers (Ghanem et al. 2017;Wang et al. 2021b). In fact, Rayleigh wave contained 1 3 in the reflection seismic data has unique advantages for MASW. The high-fold data can offer an incredible statistical redundancy for the extraction of dispersion curves (Boiero et al. 2011). The receiving arrays in reflection seismic exploration are usually long enough to record high-mode surface waves (Strobbia et al. 2010). The use of multimode dispersion curves is very helpful to improve the depth and accuracy of inversion. However, MASW is also limited in some aspects. For example, surface-wave dispersion curve inversion can only get the layer thickness (h) and v S of the stratum, but cannot obtain v P . Because of the 1D approximation of MASW, the horizontal resolution of MASW is usually limited. In addition, although the MASW results have a high resolution at the shallow surface, its vertical resolution decreases with depth increases (Xia et al. 2012;Pan et al. 2018).
It is well known that the inversion problem is ill-posed, mix-determined and nonlinear. The joint inversion and analysis of different types of geophysical data is a hot research topic in the field of geophysics (Bergmann et al. 2014;Hellman et al. 2017;Jacob et al. 2018;Schwellenbach et al. 2020). If there is a link between the inversion parameters, different types of geophysical data can be directly combined Wagner et al. 2019). The combination of different exploration methods helps improve the interpretability and resolution of the results (Vozoff and Jupp 1975;Ronczka et al. 2018; Wagner and Uhlemann 2021). FATT and MASW utilize different physical principles and suffer from some intrinsic limitations. The advantages of two methods are complementary in some aspects (Donno et al. 2021). MASW can overcome some intrinsic limitations of FATT, such as velocity-inversion interfaces, low-velocity zones and strong velocity contrasts. On the other hand, FATT provides valuable information for MASW, such as the depths of bedrock and water table, the estimation of and the lateral variation of velocity.
Recently, 1D joint inversion of first-arrival traveltimes and surface-wave dispersion curves has attracted the attention of some researchers. Foti et al. (2003) used first-arrival and surface waves to jointly interpret field data and discussed the relative limitations and advantages of the two methods in terms of shallow water table, velocity-inversion interface and hidden layer. Dal Moro andPipan (2007, 2008) applied multi-and bi-objective evolutionary algorithm to 1D joint inversion of refraction traveltimes and dispersion curves. The results indicated that the joint inversion method can mitigate the inherent ambiguity associated with the hidden layer. A 1D least-squares joint inversion method of refraction traveltimes and surface-wave dispersion curves was proposed by Piatti et al. (2012). The result of synthetic model test showed that the joint inversion can obtain more accurate velocity structure than individual inversions. The joint inversion results of field data were used to estimate the porosity and mechanical properties of the soil sediments. For laterally varying layered media, Boiero and Socco (2014) proposed a 1D joint inversion workflow to construct v P and v S structures from P-wave refraction and surface-wave data.
The 2D joint inversion usually improves the effect of inversion by incorporating the normalized cross-gradients constraint based on the structural similarity of different geophysical datasets (Yari et al. 2021). Theoretically, the essence of FATT is global 2D inversion. The inversion result of FATT is 2D v P models. However, the essence of MASW is individual 1D inversion. The inversion result of MASW is 1D v S structure. The inherent mechanisms and objective functions of FATT and MASW are intrinsically different, so traditional 2D joint inversion of them is difficult to be performed. Combined inversion or mutually constrained inversion strategies have been widely applied for different type geophysical datasets to improve the results of the inversion and interpretation (Wisén and Christiansen 2005;Qin et al. 2020). In order to integrate the complementary advantages of FATT and MASW, Ivanov et al. (2006) presented the joint analysis method of refractions with surface waves, which used a reference v P models derived from MASW to constrain the process of 1 3 refraction traveltime tomography. Duret et al. (2016) proposed a combined inversion workflow based on complementary use of Rayleigh and refraction waves. The result of field data application proved that the combined inversion workflow can depict the near-surface velocity structure well, and significantly improve the quality of subsequent data processing.
In this paper, we propose a near-surface site characterization method based on the joint iterative analysis of first-arrival and surface-wave data (JIAFS) by adding the iterative process to the mutually constrained inversion. First, the 1D v S models obtained by MASW are interpolated to construct the pseudo-2D v S model, which can effectively identify the near-surface velocity-inversion interfaces and low-velocity zones. According to the available geological survey information and borehole data, the initial model is estimated. Based on the estimated model, the pseudo-2D v S model is converted to a referenced v P model which can be utilized to constrain the process of FATT. Then, the v P model obtained by FATT is used as a priori information to improve the reliability of the results of MASW. Thus, the results of the two inversion methods can be mutually improved, and a workflow of joint iterative analysis is established. Finally, through the iteration of this process, we can estimate accurate near-surface v P , v S and models under complex geological conditions. In synthetic model test, a velocity model including low-velocity zone is designed to analyze the advantage of the proposed JIAFS method in detail. Beside, JIAFS is also applied to the seismic data acquired for oil and gas exploration in Northwest China. The results of JIAFS are compared with the results of two individual inversion methods and borehole data to illustrate the practicality and accuracy of the proposed method.

First-arrival Traveltime Tomography
The process of FATT consists of three steps: (1) picking first-arrival traveltimes, (2) ray tracing and (3) tomographic inversion. The shortest path algorithm is used to solve seismic raypaths tracing and shortest traveltime forward problem (Dijkstra 1959;Moser 1991). For a given 2D velocity model, the total traveltime is the integral over the raypaths. The underground model is divided a grid of cells, and each cell has a constant slowness. Slowness is the reciprocal of velocity. It should be noted that the size of the cell depends on the given velocity model and survey geometry. Too small cells result in the poor ray coverage and increase the non-uniqueness in FATT. On the contrary, too large cells cause inadequate resolution. Therefore, an adequate model parameterization is vital for FATT (Zhou 2003). In this paper, the triangular cells are used to divide the 2D underground model . When subsurface structural information such as the presence and orientation of lithological boundaries and abnormal velocity structures is available from MASW, the triangular cells allow flexible setting of velocity interfaces during FATT (Stenerud et al. 2009).
The integral process can be expressed as a matrix-vector product describing the traveltimes for all raypaths: where t is the first-arrival traveltime vector for the measured raypaths, s is the slowness vector containing the slowness values and the matrix L is the path matrix.
Like most geophysical inverse problem, Tikhonov regularization is utilized to solve illposed problems, and FATT inversion can be expressed as the following optimization problem: where m represents the model parameter vector. Φ(m) is the objective function. Φ d is the observation data fitting term. Φ m is the model regularization term. And r indicates the regularization parameter (Tikhonov and Arsenin 1977). It should be noted that the stability of inversion is enhanced by solving the logarithmic versions of parameters. For FATT, the model parameter during the inversion is the logarithm of velocity. The objective function can be further transformed into the following form: where ‖ * ‖ 2 2 represents the square of L2 norm. f (m) represents the forward response vector of model parameters. d denotes the observed first-arrival traveltime vector. m ref denotes the reference model parameter vector containing a priori information (Jordi et al. 2018). W and S are the weight and smoothness matrixes (Günther et al., 2006).
The optimization problem can be further expressed as the following linear equation to obtain the updated model parameters: k indicates the number of iterations. T represents the transpose operation. And J represents the Jacobian matrix including the partial derivatives of the model responses with respect to the model parameters. The LSQR algorithm presented by Paige and Saunders (1982) is used to solve Eq. (4) for iterative calculation. Finally, the iterative inversion stops until a previously established stopping criterion is satisfied, and the 2D v P inversion result is obtained.

Multichannel Analysis of Surface Waves
MASW is usually divided into three procedures: (1) seismic data acquisition, (2) dispersion information extraction and (3) dispersion curve inversion (Socco and Strobbia 2004). Surfacewave dispersion curve inversion refers to the process of finding a 1D velocity structure to optimally match its theoretical dispersion curves with the observed dispersion curves (Foti et al. 2011). For layered earth model, the phase velocity of Rayleigh wave is a function of four earth parameters ( v P , v S , is the objective function and h ) and frequency. Compared with v P and , the phase velocity is more sensitive to v S and h. Therefore, during the inversion, the parameters v P and are usually regarded as fixed known conditions which can be estimated from the results of FATT and a priori information Chai et al. 2012). The initial model usually consists of many uniform thin layers. The thickness of the thin layer depends on the minimum vertical resolution of the observed dispersion data, which is generally considered to be half of the minimum wavelength of the observed surface-wave data (Foti et al. 2018).
In this paper, we adopt the linearized dispersion curve inversion method of v S and h ). The problem of dispersion curve inversion method is similar to FATT introduced in the previous subsection and can also be represented in the form of Eq. (3), where f (m) represents the forward modeling of dispersion curves. The modeling dispersion curves are calculated by the fast-delta matrix method (Buchen and Ben-Hador 1996). d denotes the observed dispersion vector. m ref includes the a priori information obtained (4) Am k + 1 =d, from the results of FATT, geological survey or borehole data. This provides constrains on the individual inversion results of dispersion curves. W is the weight matrix of the reciprocals of data variances. And S is the smoothness, which is the discrete form of the Laplace operator (Constable et al. 1987).
The updates of the model parameters are obtained by iteratively solving the similar linear equation described in Eq. (4). Different from FATT, the Jacobian matrix J of dispersion curve inversion contains the partial derivative of phase velocity to v S and h (Wu et al. 2019b). Finally, the iterative inversion stops until a previously established stopping criterion is satisfied, and the 1D v S inversion result is obtained. Based on the 1D approximation of MASW, the inverted 1D v S structure reflects the underground structure below the linear receiver spread. Therefore, the inverted 1D v S structures are located at the midpoints of the according receiver spreads and can be interpolated to construct pseudo-2D v S model by the inverse distance weighting interpolation method.

Joint Iterative Analysis of First-arrival and Surface-wave Data
Conventionally, FATT and MASW are two individual inversion processes, as shown in is usually limited. FATT is an effective method for determining wave velocities in the presence of lateral velocity variations. This provides a favorable priori information constrain for MASW and improves the reliability of the results of MASW. On the contrary, MASW can effectively identify the near-surface velocity-inversion interfaces and low-velocity hidden layers, which is complementary to FATT.
In order to integrate the complementary advantages of the two methods, is used as a bridge to connect the two methods. The joint iterative analysis method of first-arrival and surface-wave data is proposed in this paper, and the corresponding workflow is shown in Fig. 2b. First, the 1D v S models obtained by MASW are interpolated to construct the pseudo-2D v S model. According to the available geological survey information and borehole data, the initial model is estimated. Based on the estimated model, the pseudo-2D v S model is converted to a referenced 2D v P model which can be utilized to constrain the process of FATT. Then, the 2D v P model obtained by FATT provides a priori information to constrain the individual inversions of dispersion curves. The v P and v S models obtained by FATT and MASW are used to update the model and set as the initial models in the next iterative analysis. Finally, through the iteration of this process, the two inversion methods can make use of their own advantages to improve the effects of each other, and the workflow of joint iterative analysis is established. By introducing multiple iterative operations in the process of mutually constrained inversion, accurate and reliable near-surface v P , v S and structures can be obtained.

Synthetic Model Test
For illustrating the influence of low-velocity zones on FATT, two three-layer velocity models with horst structures are established for synthetic model test. As shown in Fig. 3a, the velocity of Model A increases from shallow to deep. The blue dots at the surface of the models represent the locations of the shots and geophones. The synthetic first-arrival traveltime data are generated assuming 75 shots and geophones spaced by 2 m. The shortest path method is used to calculate the raypaths with three additional nodes on each edge of the triangular mesh to increase the calculation accuracy (Giroux and Larouche 2013). The raypath coverage and tomographic inversion result are illustrated in Fig. 3b and c, and the white lines represent the raypaths. It can be seen that the raypath coverage can uniformly cover the underground geological structures. The velocity structure interfaces of FATT result are clear, and the horst structure in the second layer can be accurately identified. The horst structure in Model B is a low-velocity zones, as shown in Fig. 3d. From the raypath coverage and FATT result in Fig. 3e and f, it can be seen that the raypaths are difficult to pass through the velocity-inversion interfaces above the low-velocity structure. This causes a shielded region for the raypaths. The velocity field in the shielded region strongly depends on the given initial velocity model and cannot be updated effectively, as the black arrow indicated in Fig. 3f. This is due to the limitations of the basic assumptions of FATT. When there is a velocity-inversion interface, the dependence of FATT on the initial velocity model can be alleviated by adding constrains of the inversion results of MASW. For comparison, Model B is also used to assess the performance of JIAFS. The corresponding true v P , v S and models are demonstrated in Fig. 4a-c. In Fig. 4d, the velocity of initial v P model increases uniformly from shallow to deep. It should be noted that the actual initial v S models are 30-layer 1D velocity models. The 2D initial v S model in Fig. 4e are obtained by interpolating all 75 1D initial v S models which are arranged at the positions of X = 2, 4, …, 148, 150 m. The fast-delta matrix method is used to calculate the synthetic dispersion curves (Buchen and Ben-Hador 1996). A total of 75 sets of fundamental-mode dispersion curves are calculated as observed data for inversion. The initial model is set to be a constant ( = 0.33), as shown in Fig. 4f. Figure 5 demonstrates the four iterations results of JIAFS. The 2D v S models in Fig. 5 are obtained through interpolating and smoothing all 75 1D inversion results. It can be seen that as the number of iterations increases, the influence of velocity-inversion interfaces on the underlying stratum is gradually mitigated. Comparing the v P results of FATT and JIAFS in Figs. 3f and 5j, JIAFS can more accurately identify the low-velocity horst structure and velocity interfaces. In addition, the v S and results of JIAFS can also match the true models well, as shown in Fig. 5k and l. To highlight the quantitative differences, the 1D v P , v S and models at X = 24, 50, 74, 100 and 126 m are extracted for comparison in Fig. 6. After four iterations, the quantitative consistency between the inversion results and the true models is significantly improved. The fit between the inverted and observed data is shown in Fig. 7. It can be seen that the inverted first-arrival traveltimes and dispersion curves can match the observed data well.
In addition, for quantitatively analyze the quality of inversion results, we introduced the definition of relative error (RE): where D true and D true represent the true and estimated data of a cell, the operator "||" stands for absolute value operation. Therefore, a smaller RE represents a more accurate inversion  Fig. 9. It can be seen that the inversion results of JIAFS can better fit the true models, especially for the hidden low-velocity structure.

Field Data Application
The field seismic data used in this paper was acquired for oil and gas exploration in Northwest China. The acquisition parameters were not originally designed to record first-arrival and surface waves. The acquisition geometry included eight receiving arrays, as shown in Fig. 10. The red pentacles and blue triangles represent the locations of sources and receivers. Each array included 35 sources and 80 receivers. The natural frequency of the geophones was 4 Hz. The seismic waves were excited by the explosive sources. According to the available information, the surface of the survey area is flat and covered by loose Quaternary sediments. The depth of the water table varies from 10 to 60 m. Due to the complexity of the geological conditions in the shallow underground, especially the existence of velocity-inversion interfaces, the seismic imaging of some deep targets is challenging. The proposed JIAFS method is applied to field seismic record to characterize the near surface by using the first-arrival and surface-wave data. In Fig. 11a, the seismic record shows obvious Rayleigh and first-arrival waves. The recording length was 7 s in time domain. The time sampling interval was 2 ms. The picked first arrivals are denoted by the red lines. In the spatial windows demonstrated by the green dashed frames A and B in Fig. 11a, the seismic records according to the linear  Fig. 11b and c. The adopted spatial window contains 35 seismic traces, covering the main distribution range of Rayleigh waves. The multimode dispersion curves represented by the white lines are automatically picked by the clustering algorithm (Wang et al. 2021a). The picked first arrivals and dispersion curves are used as observed data in subsequent inversions. Based on the 1D approximation of MASW, every inverted 1D v S structure reflects the underground structure below the linear receiver spread. Each seismic shot can obtain two 1D v S structures corresponding to the positive and negative offsets. The two 1D v S structures are located at the midpoints of the spatial windows. In Fig. 12, the orange and purple points represent the spatial distributions of all inverted 1D v S structures obtained from positive and negative offsets.
The inversion results of FATT,MASW and JIAFS are compared in Figs. 13,14,15,16,17 and 18. Due to the insufficient ray coverage on both sides, the X coordinate ranges of the 2D inversion results correspond to the green dashed frame in Fig. 10. The 2D inversion models of JIAFS are the results after three iterations, because the subsequent iterations no longer significantly improve the inversion results. From the inversion results, it can be seen that the overall trends of the velocity models obtained by the three methods are consistent. From shallow to deep, the 2D velocity models can be divided into low-velocity, medium-and high-velocity layers. The joint analysis of FATT and MASW reduces the non-uniqueness of the inverse problem and improves the reliability of the inversion results.
Since the v P of shallow stratums is greatly affected by the water saturation, the v P is very high when the water is saturated below the water table. However, the v S is mainly related to lithology and less affected by water saturation (Gomes and Lopes 2014). Therefore, has a more intuitive reflection on the porosity and water saturation of stratums. We demonstrate 2D models obtained by individual inversions and JIAFS in Figs. 17 and 18. Obviously, the lateral variation of is significant. The maximum is about 0.47 due to the shallow water table and loose surface sediment. In a certain depth range below the water table, the water saturation of the stratum is very high, and the is greater than 0.465. As the depth gradually increases, the stratum pressure increases, the porosity and water saturation of the stratum gradually decrease, and the also starts to decrease from about 0.465.
To verify the reliability and accuracy of the inversion results, the borehole data at the X coordinates of 4900 and 7700 m are embedded into the 2D inversion results. It can be seen that, compared with FATT and MASW, the velocity interfaces of the inversion results of JIAFS are more consistent with the lithological stratification of the borehole data. The 1D velocity models in Fig. 19 are extracted from the positions where the X coordinates are 4900 and 7700 m in Line 1. The blue and red solid curves indicate the v S models obtained by MASW and JIAFS. The blue and red dotted curves indicate the v P models obtained by FATT and JIAFS. In contrast, the results of JIAFS match the borehole data better than individual inversions, and the thin gravel layer can be accurately identified in 1D v P and v S result, as shown in Fig. 19b. In addition, we interpolate the 2D Fig. 6 a, b and c Extracted 1D v P v S and models of the first iteration; d, e and f Extracted 1D v P , v S and models of the second iteration; g, h and i Extracted 1D v P , v S and models of the third iteration; j, k and l Extracted 1D v P , v S and models of the fourth iteration 1 3 inversion results in the Y direction to obtain the 3D v P , v S and models in Figs 20 and 21. The sections in any direction can be obtained by slicing the 3D models, which is of great significance for the near-surface site characterization and subsequent processing of seismic data. Compared with FATT and MASW, the 3D results of JIAFS have higher resolution for geological structures and are more consistent with a prior information. This further verifies the practicability of JIAFS.

Discussions and Conclusions
Accurate near-surface site characterization and velocity modeling are critical in many fields, such as geotechnical engineering, soil mechanics and oil and gas exploration. At present, FATT and MASW methods are widely used for near-surface velocity modeling and have respective advantages and intrinsic limitations. The 1D assumption of MASW limits its horizontal resolution. FATT cannot effectively identify the low-velocity hidden layers. However, because of the differences in the physical mechanisms of inversions, 2D joint inversion is not applicable for the two methods. In order to combine the complementary advantages of the FATT and MASW, a near-surface site characterization method based on the joint iterative analysis of first-arrival and surface-wave data is proposed in this paper. The synthetic model test and field data application were used to verify the performance of the proposed JIAFS method. When there is a low-velocity zone, the limitations and advantages of FATT and MASW were demonstrated through the synthetic model test. The result of synthetic model test demonstrated that the strategy of JIAFS can combine the advantages of FATT and MASW to obtain more accuratev P and v S models. Besides, is an important mechanical parameter of near-surface materials. At present, the method to obtain the near-surface models is mainly through the petrophysical experiments. But this is expensive for large area site characterization. According to the joint analysis of first-arrival and surface-wave data, accurate near-surface models can also be established. JIAFS was also applied to field seismic data acquired for oil and gas prospecting in Northwest China. The comparison between the results of JIAFS, FATT and MASW illustrated the superiority of JIAFS. The reliability and accuracy of the JIAFS results were verified by the borehole data. The 3D v P , v S and models were obtained by interpolating the 2D inversion results. This is of great significance for geologists to accurately characterize near-surface sites.
For subsequent research, the 3D FATT could be considered to be introduced into the process of JIAFS. It should be noted that, since we introduced the iterative process into mutually constrained inversion, this increases the extra computational cost compared to the individual inversions. Although the increased computational cost in the 2D joint iterative analysis is acceptable, for 3D joint iterative analysis, the computational efficiency and cost of the algorithm must be considered. On the other hand, the dispersion curve inversions are constrained separately in this paper. For geological environments where the lateral variations are not too pronounced, laterally constrained inversion strategy could also be considered in JIAFS to help geological interpretation (Socco et al. 2009). Besides, the v P and v S models obtained by JIAFS can be applied to the subsequent processing of seismic data to further verify its advantage, e.g., statics correction, stacking and migration imaging.