Imaging the Fault Zone Structure of the Pearl River Estuary Fault in Guangzhou, China, from Waveform Inversion with an Active Source and Dense Linear Array

Abstract

Since high-resolution structure imaging of active faults within urban areas is vital for earthquake hazard mitigation, we perform a seismic survey line crossing the Pearl River Estuary Fault (PREF) in Guangzhou, China. First, ten shots of a new and environmentally friendly gas explosion source are excited with about 1 km spacing and recorded by 241 nodal short-period seismometers with an average spacing of 60 m. Then, based on these acquisition data, we adopt waveform inversion to explore the kinematic and dynamic information of early arrival wavefields to recover the subsurface structures. The inversion results indicate that while the low-velocity zone (LVZ) in depth surrounding the PREF is 2.5 km in width and extended to 0.7 km, another LVZ of 1.5 km in width and extended to 0.7 km in depth is surrounded by the Beiting–Nancun fault. We observe that the analysis of evolution and activities of the fault systems reveal no historical earthquakes in our study area; we interpret that the two LVZs controlled by the faults are probably attributed to the fluid dynamics, sediment source, and fault motion at different geological times, rather than fault-related damage zones. The results can provide significant basis for earthquake prevention and hazard assessment in Guangzhou. The finding also shows that the waveform inversion can effectively explore the fine structure of active faults in urban area with dense linear array and spare active source excitations.

1. Introduction

Since active faults, usually accompanied by observable movement or seismic activities during the last 10,000 years, are significant for earthquake hazard mitigation, fine imaging of the fault zone structure can help us understand the earthquake activity and potential earthquake hazard of the region near the fault better [1,2,3]. In particular, the potential threats could even be serious when the active faults are located within urban areas. Therefore, it is imperative to identify the fine structure of significant faults in urban areas to support the design of city infrastructure, earthquake hazard mitigation plans, and other factors that could affect urban safety [4].

Many geophysical prospective methods have been successfully developed in urban areas to detect the fault zone structure, including gravity, magnetic, and seismic methods [5]. By arranging dense linear arrays across the fault zone, we can conduct local and teleseismic earthquake travel time analysis [6,7,8], waveform modeling of fault zone trapped and head waves [2,9,10], and local earthquake travel time/ambient noise tomography [11,12,13]. However, while the above methods highly depend on the data quality of small-aperture arrays that can record the available phase and earthquakes within about 1–2 months, ambient noise tomography are challenging to robustly constrain depth [14].

Artificial seismic reflection with a large aperture has been adopted as an alternative method to identify fault characteristics [15,16,17]. Till date, many artificial sources, including explosive sources, vibrators, and methane gas explosion sources have been used for field detection [18], and to meet the increasing environmental protection requirement in urban areas, a new type of active source, the methane gas explosion source, is used [19]. As the methane gas explosion source uses gas-phase detonation reactions to generate seismic waves, the reaction products–water and carbon dioxide, cause little pollution. Therefore, the characteristic of methane gas has been analyzed in field experiments. Investigations revealed that since its excitation signal had enormous energy, a wide frequency band, and a long propagation distance (15–20 km), the first break and direct wave signals could be recorded with a high signal-to-noise ratio (SNR) [13]. Hence, combining the long propagation distance merits of methane gas explosion sources with several short-period instruments, deployment in a 2D line or 3D array with a receiver spacing of about 50–200 m, maintenance of sufficient coverage and wide aperture, and detection of the width and depth of fault zones are possible.

In previous research works, while the data have conventionally been exploited through travel time tomography, main discontinuities have been estimated by minimizing the misfit function of the observed and calculated first arrival travel times [20]. However, travel time tomography assumes that since the seismic data are infinite and of high frequency, the wave propagation would only be affected by the ray path, independent of any changes outside the infinitely narrow ray path, which violates the finite-frequency bandwidth of the seismic data as it uses only the travel time of the first arrival and discards abundant information in amplitude and phase [21,22]. Consequently, this method results in a suboptimal estimate of the earth’s velocity [23,24].

Compared to travel time tomography, waveform inversion utilizes much more characteristic information in seismogram records. By minimizing the misfit function that measures the difference between the waveforms of the observed and calculated seismograms, a very accurate estimate of both the coarse and finely detailed features in the actual model is presented [25,26]. Recently, field applicability of waveform inversion has also increased, attracting significant attention of researchers, for example, in retrieving imaging crustal structural details in southern Tibet [27] and time domain waveform inversion to reconstruct the shallow structures onshore in Argentina [28]. Though waveform inversion can provide high-resolution model parameters of a subsurface, its wide field applications have been impeded by some problems, including the strong nonlinearity between the wavefield and model parameters, the nonconvexity of misfit function, the complexity of the propagation medium, and so on [29,30]. Hence, some researchers have proposed promoting the field applications of waveform inversion, such as early arrival waveform inversion [21,31], adaptive waveform inversion [32], the construction method of the large-scale initial input model [33], and the envelope waveform inversion approach [34].

In June 2020, we deployed a temporarily linear dense array crossing the Pearl River Estuary Fault (PREF) in Guangzhou, China (Figure 1). Since urban environments pose many challenges to data acquisition: interference (electromagnetism and human activity noise) and limited space (dense buildings and hardened ground), ten methane gas explosion sources were excited along the array line at a maximum offset of 12.5 km. Then, based on this active source data, we attempted to reconstruct the high-resolution subsurface structure of PREF and describe a suitable workflow for fault zone detection in urban areas. This procedure used travel time tomography methods to obtain an initial model near the global minimum and applied waveform inversion using the early arrival wavefields of the observed and calculated data.

The rest of the article is organized as follows: Firstly, we introduce the field experiment as a linear survey line conducted at Guangzhou, China, after which the geological setting, acquisition geometry, and data pre-processing were presented. Next, we briefly cover the inversion theory in the time domain, introducing the overthrust velocity model to demonstrate the waveform inversion’s validity using the early arrival wavefields. Finally, the travel time tomography and field inversion results are presented, followed by the description and discussion, and finally, the conclusions.

2. Data

2.1. Geological Setting

The seismic data processed in this paper were collected from the alluvial plain of the Pearl River Delta (PRD). PRD is located in the tectonic faulted basin of the southeast coast, China, and has experienced the Hercynian–Indosinian, Yanshan, and Himalayan geological movements. This region also includes three main groups of faults: North-West (NW), East-West (EW), and North-East (NE), which intersect each other, cutting the delta basement into fault depressions and uplifts with different sizes and movement rates [35]. Notably, since the NW trending faults were mainly distributed in small coastal areas, generally extending within 100–200 km, they formed later than the other fault systems, mainly in the late Yanshanian and Himalayan periods. As a result, they cut off almost all NE and EW faults, thereby controlling coastal water systems and the formation/development of coastal harbors [36].

While its east and west sides are, respectively, the Pearl River Estuary Graben and Panyu fault uplift (PY), PRD is closely related to regional seismic activities and the formation/evolution of deltas [37,38]. Furthermore, with a total length of about 43 km, PREF starts from Jishan in the north to Lingdingyang in the south, cutting the Proterozoic metamorphic rock, Yanshanian granite, and Cretaceous clastic rock basement [35]. However, it is still weakly active, belonging to a brand-new active fault. Additionally, while a sub-grade NW trending small-fault Beiting–Nancun fault (BNF), which was active in the Late Pleistocene of the PY, showed no obvious signs of activity after entering the Holocene [39], its survey surface was covered by quaternary systems, showing no bedrock exposure. In terms of landform, however, while the west side was within a general elevation of 20–40 m, the east was a delta accumulation plain area, with flat terrain and surface elevation within 10 m.

2.2. Data Acquisition and Pre-Processing

In the mid-2020s, we performed a seismic acquisition survey across PREF and the south of the BNF in Guangzhou, China (Figure 1a). We observed that since the deviation was not too large, the dense linear array, involving 241 short-period instruments spaced at an average of about 60 m apart, was a straight line (blue triangle in Figure 1b). Furthermore, since the profile struck North-East–South-West and was almost perpendicular to the PREF, ten shots were fired with an average spacing of 1 km into the array by explosive housing charges in boreholes about 10 m deep (red circle in Figure 1b). Consequently, the maximum offset reached about 12.5 km, causing the average recording time intercepted from the continuous data of the short-period seismograph according to the excitation time to be 2.1 s and the sampling rate to be 1 ms. Evidently, the seismic data recorded by the short-period instruments showed high-quality first arrival signals (Figure 2a).

When working with this field data, it is necessary to improve the data quality before inversion [40]. Therefore, we designed a pre-processing procedure for this field data, including quality control procedures to remove the low signal traces and to correct the reversal polarity, band-pass filtering (20–100 Hz) and curvelet denoising [41] to effectively suppress the noise, amplitude normalization to partially compensating for the source size variability and receiver-ground coupling/attenuation effects. By analyzing the seismogram records, it was evident that the early arrivals, defined as those events that arrive within a few periods of the first arrivals, were within high SNR. Hence, we designed time windows to truncate the early arrival wavefields and implement waveform, with most of the ten shots gathered being good in quality. Figure 2 shows the field shot gathered for a single shot located 9 km before (Figure 2a) and after pre-processing (Figure 2b).

3. Method

3.1. Waveform Inversion Theory

Generally, waveform inversion is a nonlinear optimization problem that seeks to minimize misfit functions. It begins from an initial input model, calculates the search direction and step length, then updates the target model iteratively, finally enabling the derivation of the high-resolution subsurface model by exploiting full information content from the observed seismic data [29]. Hence, considering the propagation characteristics of this field data and detection aims, we used only the early arrival wavefields to construct a misfit function (Equation (1)), utilizing the kinematic and dynamic information of early arrival wavefields to recover the velocity model accurately. With the early arrival, inversion processes reduced nonlinearity and complexity, presenting a more reliable convergence than conventional waveform inversion [21,42].

The goal of waveform inversion is to minimize the misfit function $f$, which measures the observed and calculated data. For the field data, we used the L1-norm misfit function that is weakly sensitive to noise in the inversion process and can provide more reliable results, calculated as follows [43]:

$$f=\left|{\displaystyle \sum _s{\displaystyle \sum _g\left[{{\displaystyle p}}^{cal}\left(x,z,t\right)-{{\displaystyle p}}^{obs}\left(x,z,t\right)\right]m\left(x,z,t\right)}}\right|$$

(1)

where ${{\displaystyle p}}^{cal}\left(x,z,t\right)$ and ${{\displaystyle p}}^{obs}\left(x,z,t\right)$ represent the calculated and observed data, respectively, and $m\left(x,z,t\right)$ represents a mute function that silences all energy, except for the early arrival wavefields; $s$ and $g$ represent the source and geophone, respectively. Subsequently, the calculated data ${{\displaystyle p}}^{cal}\left(x,z,t\right)$ could be simulated by implementing forward modeling in the time domain using the following acoustic wave equation:

$$\frac{1}{{{\displaystyle v}}^2\left(x,z\right)}\frac{{{\displaystyle \partial }}^{2}p\left(x,z,t\right)}{{{\displaystyle \partial t}}^2}={{\displaystyle \nabla }}^{2}p(x,z)+s\left(x,z,t\right)$$

(2)

Next, a time domain finite-difference scheme with tenth-order in space and second-order in time was used to solve Equation (2), after which the Perfectly Matched Layer absorbing boundary condition [44] was implemented to prevent reflections and diffractions at the boundary effectively. For simulated experiments, while the source time function $s$ in Equation (2) was definite, in field data, the source time function was always unknown. Therefore, a source time function assumption procedure is inevitable to make the inversion more reliable and improve the results [45]. For this field data inversion, the source was estimated using a method similar to that of Pratt [46]:

$$s=\frac{{\left(p^{cal}\right)}^t{\left(p^{obs}\right)}^{*}}{{\left(p^{cal}\right)}^t{\left(p^{cal}\right)}^{*}}$$

(3)

where superscript $t$ represents the transpose, and superscript $*$ represents the complex conjugate.

While we implemented the waveform inversion in the time domain, we also needed to estimate the source wavelet in the frequency domain and then inverse it for the time domain [28]. To this end, the misfit function gradient was efficiently calculated via the adjoint-state method, as adopted in most cases [47], and is expressed as shown below:

$$g=\frac{2}{v^3}{\displaystyle {\int }_0^T\lambda \frac{{\partial }^{2}p\left(x,z,t\right)}{\partial {{\displaystyle t}}^2}}dt$$

(4)

where $T$ is the maximum time recorded on the dataset, $p$ is the forward-propagated source wavefield, and $\lambda$ represents the back-propagated residual wavefield calculated by Equation (2).

We then estimated the back-propagated residual of the L1-norm misfit function, as shown below:

$$backc\left(r\right)=sign(r)=\left\{\begin{array}{l}+1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}r>0\\ 0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}r=0\\ -1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{r}<0\end{array}\right\}$$

(5)

where $\hspace{0.17em}r$ represents the wavefield residual error, which is set as the source in Equation (2) to calculate the back-propagated residual wavefield.

In this process, the parameters were updated iteratively to achieve calculated data from the velocity model that can match the observed data using the expression shown below:

$${{\displaystyle v}}_{k+1}={{\displaystyle v}}_k+{{\displaystyle \alpha }}_k{{\displaystyle h}}_k$$

(6)

${{\displaystyle v}}_{k+1}$ and ${{\displaystyle v}}_k$ represent the (k+1)th and kth iteration model parameters, ${{\displaystyle \alpha }}_k$ represents the step length at the kth iteration, and ${{\displaystyle h}}_k$ represents the search direction at the kth iteration.

Finally, while step length is usually calculated using line search [48], the search direction was calculated using gradient methods, such as the conjugate gradient and quasi-Newton methods, such as the limited memory BFGS (L-BFGS) method, with the inversion process stopping when the misfit function decrease no longer exceeds the threshold, or when the maximum iteration number is achieved.

3.2. Synthetic Example: The Overthrust Model

A synthetic example was subsequently adopted to verify the validity of the waveform inversion using the early arrival wavefields. The overthrust velocity model (Figure 3a), comprising a linearly increasing velocity model as the initial input model (Figure 3b), is shown in Figure 3. Investigations revealed that while the synthetic data generated from the 30 shots were separated by 500 m, the receivers were regularly distributed in the horizontal direction at 25 m intervals and a maximum offset of 15 km. Hence, we used a Ricker wavelet with a 10-Hz peak frequency as the source function, with the time step set to 2 ms and the total recording time set to 6 s.

Here, the early arrival wavefields were defined as those events that arrived within a few periods of the first arrivals. The truncated early arrival waveform of the synthetic data generated from the overthrust velocity model (Figure 3a) is shown in Figure 4. It was evident that the truncated early arrival wavefields (Figure 4b) mainly contained the direct and refracted waves.

Though we only adopted early arrival wavefields in this process, the misfit was more linear for the velocity than conventional full-waveform inversions, indicating that the high-frequency data used could cause the misfit function to be highly nonlinear and inversion to suffer from the local minimum problem [22]. Therefore, we used a multiscale inversion strategy for time domain waveform inversion to further promote waveform inversion converging to the global minimum. Then, we adopted the Wiener filter function to filter the data and obtain different frequency band datasets for the multiscale inversion [29], which sequentially applied waveform inversion to three frequency band-pass ranges with the following cutoff frequencies: [0.5–1.6], [1.5–4.9], [4.8–16.3]. The final inversion result is shown in Figure 5. It was evident that the inversion adopted early arrival wavefields to reconstruct the layer structures and the top thrust fault of the overthrust model. However, we also observed that the left and right deep parts were slightly and poorly inverted, considered to be due to an inaccurate initial velocity model and insufficient data coverage.

For more details, we extracted the profiles of velocity variation profiles with a depth at positions of 5 km (Figure 6a) and 10 km (Figure 6b). Evidently, the waveform inversion recovered the high accuracy structure of the overthrust velocity model within a 1.5 km depth (Figure 6), fully showing that waveform inversion using early arrival wavefields could recover the subsurface structure with high resolution.

4. Application of the Waveform Inversion

4.1. First Arrival Travel Time Tomography

For the field data set, combining waveform inversion with travel time tomography was a reliable choice that could generate a sufficiently accurate initial velocity model. Hence, using this large-scale velocity model as the input, we further mitigated the nonlinearity of this complex inverse problem to reduce the risk of falling into a local minimum [49]. Travel time tomography usually requires tedious travel time-picking, the initial model of the one-dimensional parallel layer velocity was approximately derived based on an analysis of the time-offset relationship for the dataset and according to the picking first arrival travel time. As a result, we obtained the tomography results by minimizing the misfit function of picking the observed and calculated travel times, as shown in Figure 7a, including the travel time-fitting of the corresponding shot, as shown in Figure 7b and the P-wave velocity disturbance results after five iterations, the travel time distributions of the corresponding seismograms, as Figure 7c shows. While the inversion travel time of each shot fitted well with the picking travel time, the travel time tomography recovered the long-wavenumber background velocity structure. We also observed that the detailed features, such as the short- and medium-wavelength velocity perturbations, were blurred. In any case, this tomography result was essential as the initial velocity model for waveform inversion.

4.2. Waveform Inversion Results

To balance the computation cost with the inversion accuracy, this inverse problem was solved by adopting the L-BFGS optimization algorithm [50] until the calculated early arrival wavefields satisfactorily match those of the field data. However, while the inexact line search method has also been adopted to compute the step length [48], the layer-stripping approach that includes the frequency, offset, and time-dapping is usually adopted to hierarchically mitigate the nonlinearity in inversion processes [51]. Hence, our study used different frequency bands to invert the velocity model and sequentially applied waveform inversion to three frequency bands with the following cutoff frequencies: [1.4–4.9], [4.3–14.7], and [14.4–48.9]. Then, while the frequency bands were calculated using the method based on the Wiener filter function [29], the output velocity model at the lower frequency band was set as the input model for the higher frequency band. Subsequently, the background model obtained from travel time tomography (Figure 7a) was set as the initial input model. Investigations revealed that when the multiscale frequency strategy with an accurate initial model decreased the risk of cycle skipping efficiently, the maximum iteration number was seven at each frequency band. The inversion flowchart for this field data inversion is shown in Figure 8.

Multiscale velocity models of this survey line obtained from early arrival waveform inversion with three frequency band selection strategies are shown in Figure 9. From top to bottom, the results of the first (Figure 9a), second (Figure 9b), and third (Figure 9c) frequency bands are presented, respectively. We discovered that the lower frequency band stage tends to recover the large-scale structure because of its more linear relationship with the low wavenumber component of the velocity model. Furthermore, with an increase in frequency, gradual and clear characterization of the increasing number of apparent features was also possible, yielding high-resolution results at the final sequential inversion strategy stage.

In comparison to the velocity model obtained from travel time tomography, waveform inversion significantly improved the spatial resolution of the subsurface structure, recovering a thin layer with low velocity (2000–2500 m/s) at the top of the model. Beneath this near-surface low-velocity layer, the structures were cross-cut by multiple faults, with velocities ranging from 3000 m/s to 6000 m/s. Consequently, some detailed features were recovered that had not been constrained strongly by the travel times, such as the distinctly high-velocity anomalies at 2 and 9 km (6000 m/s). Meanwhile, note that the most noticeable structures were the two low-velocity zones (3800–4500 m/s), whose widths were marked. The deepest layer of the model revealed high velocity (5400–6200 m/s), as detected in the granitic bedrock.

5. Discussion

Our investigations showed that the results from the quaternary sediment bottom boundary delineated by a thin layer of waveform inversion agreed well with the sediment thickness results calculated during a previous Horizontal Vertical Spectral Ratio research [13]. They also discovered that the sedimentary thickness had gradually thickened from west to east since the Late Pleistocene. Conversely, after the Mid-Holocene, investigations showed that the coastline retreated southward, and the fault fell below the sea level and accepted sedimentation, causing the modern sediments near the PREF to be superimposed on the fault block [35,39,52]. This finding was highly consistent with a previous geological survey and borehole measurement, indicating that while the fault was generally concealed and hidden under quaternary sediments, only a small part was exposed [53].

The small wavelength content of this inverted velocity model provides new insights into the structure of the PREF and BNF. It was evident from Figure 10b,c that the PREF comprised two main branches: the Hualong fault (HLF) to the east and the Wenchong fault (WCF) to the west. At the right part of the velocity model, a clear LVZ (3800–4500 m/s) was observed with about 2.5 km in width and 0.7 km in depth under the quaternary surrounding the PREF. Investigations also showed that at a horizontal distance of 3.0 km, another LVZ (3800–4500 m/s) with about 1.5 km in width and 0.7 km in depth surrounded the BNF, with the velocity difference reaching 25%–40% inside and outside the two LVZs. Generally, the widths of the fault damage zones after earthquakes are on the order of 100 m, such as the ruptured structure of the Longmenshan Fault after the 2008 Wenchuan earthquake and the fault damage zones of the 1992 Mw7.3 Landers earthquake [54,55]. In addition, there is no historical earthquake record around the area that investigated the dense linear array, it is an earthquake gap in this region [56]. Our study showed that such LVZ width surrounding the fault was probably attributed to the fluid dynamics, sediment source, and fault motion at different geological times. Since the Late Pleistocene, while PRD has experienced substantial subsidence, the formation and development of the block were strongly related to the basement faults due to the great differences in fault activity, showing that while the faults controlled the sedimentation process of different structural units, changes in the sedimentary layers and thickness were affected by the differential uplift and subsidence of the underlying fault blocks. Accordingly, a study previously reported that the high-velocity structures and anomalies in the velocity profile were mainly granite bodies intruded by basic dikes and related to the granitic bedrock layer formed by tertiary magmatism [38].

Figure 10c denotes the imaging results obtained using ambient noise double beamforming within short-period dense arrays along this seismic profile. Wu et al. [38] discovered that while the east side of HLF had an obvious LVZ (dotted line), a velocity reduction surrounding the BNF was also present, though it may be unclear (the arrow indicates). These ambient noise imaging results confirm the reliability of our waveform inversion velocity model.

Hence, in the future, where the maximum earthquake magnitude in Guangzhou has been estimated to be 6.0 ≤ M < 6.5 [36], the discovered LVZ associated with the BNF and PREF may be a weakening layer, exerting significant influences on the amplification of ground motion during future earthquakes [57,58]. Our high-resolution results, therefore, provide an important factor for earthquake hazard assessment, disaster reduction, and tectonic evolution research of the delta.

Though this study also established that the inversion result using the early arrival wavefields (Figure 10b) indicated a finer velocity structure than travel time tomography, the resolution may not be satisfied because of the sparse excitation (about 1.0 km) in urban areas, causing only the use of the early arrival wavefields. To this end, integrating travel time tomography with waveform inversion using the early arrival wavefields of methane gas explosion source was considered an effective approach and workflow for resolving high-resolution fault zone structure.

6. Conclusions

This study successfully derived a high-resolution subsurface structure by waveform inversion using the early arrival wavefield data excited by a methane gas explosion source, followed by recording with a dense linear array deployed across the northern portion of the PREF in Guangzhou, China. The waveform inversion results showed an obvious velocity difference between the PREF concealed under the quaternary system and the clear LVZ of the fault that was 2.5 km in width, extending to 0.7 km in depth. Outstandingly, at about 3.0 km in the horizontal direction, another LVZ of 1.5 km in width that extended to 0.7 km in depth, surrounding the BNF, was also indicated clearly, with the velocity difference being able to reach 25%–40% for inside and outside the LVZs. From our findings, we considered that these two LVZs likely represented nonuniform sediment distributions as a controlling factor during the formation and evolution of PRD, PREF and BNF. Hence, since the LVZs associated with fault fracture systems dramatically amplify ground motion during earthquakes, they could be important in understanding earthquake physics and hazard assessments.

Summarily, the inversion results presented significant insights into this application’s fault structure, effectively improving our understanding of the LVZ surrounding the fault system in Guangzhou. Meanwhile, combining methane source sparse firing with short-period dense array observation has been confirmed as an effective method for detecting the fault zone structure in urban areas. Therefore, we have presented successful waveform inversion processes to image the LVZ structure, using field data of this acquisition, considering that this acquisition and inversion methods should have broad applications in other cities.