Application of teleseismic tomography to the study of shallow structure beneath Shizigou in the western Qaidam basin

1.   Introduction

The studied area is located on the southern margin of Yingxiongling depression in the western Qaidam basin, and belongs to the middle portion of ShizigouYoushashan anticline, a third-order structure of western Qaidam depression (Figure 1). It is adjacent to the structures of Qigequan and Hongliuquan on its south, and to the slope of Altyn on its west. Moreover, it is next to the Youshashan oil field on its east, and to the Ganchaigou structure on its north.

Figure  1.  Location and tectonic settings of the studied area (denoted by a rectangle in the figure).

DownLoad: Full-Size Img PowerPoint

The previous studies (Li and Wang, 2001; Yin et al, 2006) showed that the Shizigou region of Qaidam basin is a potential site for the accumulation of oil and gas. This region has favorable hydrocarbon accumulation conditions and also belongs to a middle-depth oil and gas reservoir, therefore greatly potential for the exploration will result. Due to the fact of complication in structure, high undulation in topography and harsh geological surface condition, a satisfactory result has not been achieved by the near vertical controlled source seismic reflection profiling method during the past 50 years and other non-seismic methods could not meet the exploration requirement as well. So, the deep structure in this region is not well known, which seriously affect the progress in survey, exploration and assessment of the resource potential of oil and gas in this region. In order to overcome those difficulties, we try to gather information of deep structure in this region by the passive source seismic method for the first time. Compared with controlled source seismic methods, passive source seismic method has some merits such as lower exploration cost and greater exploration depth. This paper gives reconstruction of the three-dimensional seismic velocity structure by using the seismic tomography method to reveal how the structure extends to a great depth, thus it can provide some information for investigation and evaluation of oil and gas resources in the region.

Since medical tomography was successfully introduced to seismology by Aki et al(1976, 1977), it has made significant progress in application to seismological researches and lots of important results have been obtained in the study of deep structure (Shi et al, 1999; Zhao and Kayal, 2000; He et al, 2007). The technology of seismic tomography so-called can show underground structure visually on the image, so it has been highly appraised and widely applied to the researches on solid geophysics and geodynamics, energy exploration and exploitation, engineering and hazard prevention, and metal deposit exploration (Poupinet et al, 2003; Zhao et al, 2007; Fu et al, 2001; Qin et al, 2006; Wang et al, 2006; Shi et al, 2004). Among various methods of seismic tomography, seismic traveltime tomography (Rawlinson and Sambridge, 2003; Thurber, 2003) is developed earlier, and relatively mature. It is also the major method widely used in recent years. The present paper will try to obtain seismic images of the target region using teleseismic traveltime tomography with body wave.

2.   Inversion of teleseismic P wave tomography

2.1   Data acquisition and analysis

The structure beneath the Shizigou region is completely an anticline. The studied region (Figure 1) covers a range of ~200 km² in Shizigou. To delineate the deep structure of this region, a network consisted of 52 stations (Figure 2) equipped with portable high frequency sensors (Reftek 130+Mark product 2 Hz) was deployed in several batches. By redeploying the stations, we could reduce the station-to-station space and improve the resolution of the resultant tomograms.

Figure  2.  Location of the seismic stations in the studied area.

DownLoad: Full-Size Img PowerPoint

This study processed the dataset recorded by the stations during October of 2004 to August of 2007. Based on the earthquake catalogue released by the U.S. geological survey (USGS), all the data were extracted with a length of 60 s (with 20 s before and 40 s after the predicted arrival time computed with the IASP91 model) for each seismic event. Subsequently, all the data consisting of 565 teleseismic events were picked up event-by-event. The waveform correlation technique was used to pick up the P wave arrival time. To achieve a high quality of inversion, we selected the events recorded by the instruments with good time characteristic, and those with high signal to noise ratio and clear first arrival of P wave. At the same time, we removed the data recorded by less than 5 stations and with epicentral distance (Δ) less than 20° or approximate to 145° from the inversion data set (in order to decrease the influences from the surface effects and those from the different branches of phase PKP). Following this procedure, a total of 2 524 arrival times from 284 events were used in this seismic tomographic inversion.

The teleseismic events used in this tomographic inversion (Figure 3) distribute annularly around the studied region. Although most of the events locate in Japan and Southeast Asia area, most azimuths have been covered with events. The epicentral distance differs, and the seismic azimuth varies in the data set. In addition, the ray paths were densely crossed. As we known, all of these are premising to get good quality of results in the seismic tomographic inversion.

Figure  3.  Epicentral distribution of seismic events used in the inversion. The dots indicate the epicenter and the rectangle denotes the location of the target region.

DownLoad: Full-Size Img PowerPoint

We threw off some particular error data according to the residual curves for the observed and calculated P wave traveltimes (Figure 4). Figure 4 shows some examples of P wave traveltime residual curves, which demonstrates significant variation of travel-time and implies a significant change in velocity beneath the studied region accordingly. The positive relative residuals in Figure 4 indicate a low velocity anomaly corresponding to relatively lower altitude region in the southwest of the studied region, whereas the negative relative residuals illustrate a high velocity anomaly corresponding to higher altitude. This suggests that the observed residuals are mainly caused by the inhomogeneous velocity rather than the altitude of the stations. We can also infer that the velocity in the subsurface medium varies intensely based on the fluctuation of P wave relative residuals. The P wave relative residuals in Figure 4 also show some obvious variation with the epicentral distance and azimuth.

Figure  4.  The fitness of the observed and calculated traveltime residuals. The circle and the plus denote the observed and calculated traveltime residuals, respectively. The epicentral distance (Δ) and back azimuth (α) of the corresponding seismic event is shown on the top-left corner of each panel.

DownLoad: Full-Size Img PowerPoint

2.2   Inversion theory and method

The ACH method (Aki et al, 1977) is a mature one in teleseismic traveltime tomography for body waves, which is used to inverse velocity structure with P wave traveltime residuals recorded by seismographs, according to the relationship between seismic traveltime and velocity of underground medium. P wave traveltime residuals reflect velocity variation of the medium beneath the receiver array, which should be an integrated effects of all the residuals derived from the entire inhomogeneous medium along the ray path (Aki and Richards, 1980).

In a continuous medium, the traveltime of a seismic wave between the source and the receiver is given by the integral:

where L is the ray path from the source to the receiver, and v(x) is the velocity. According to the ACH inversion method, the parameterization of the medium is to divide it into N homogeneous blocks. Suppose L_n is the path length of the nth block that ray passed, and then the traveltime residual δt_n within the block is

where v_t is the true velocity, v_r is the referred velocity within a block of the model and δs_n is the slowness perturbation of the nth block. Thus, the traveltime residual of seismic wave path from the source to the receiver is given by:

If the total number of observed rays is K, given d^T=(δt₁, δt₂, …, δt_K), m^T=(δs₁, δs₂, …, δs_N) and G_ij=δL_n^ij(i=1, 2, …, K; j=1, 2, …, N), then equation (3) can be written in a matrix form Gm=d, where Gis a K×N matrix, m is a N×1 matrix and d is K×1 matrix. This equation can be solved by the damped least squares:

where θ is the damping factor, thus (4) is the basic formula in the ACH inversion.

2.3   Determination of the inversion parameters

Among the methods of seismic tomography, the ACH inversion is the most robust one. Owing to the turning points of the ray excluded from the model in this method, the initial model has little effect on the inversion results. However, the parameterization of the model, to a great degree, affects the results of the seismic inversion. To obtain a better parameterization of the model, Evans and Achauer (1993) used an offset and average method that made layers thinning vertically and smoothed the model horizontally. In the vertical direction, they thinned the model blocks to make the boundary of the blocks match the actual velocity interface. Besides, each layer of the original model was offset horizontally and a weighted average method was used for the inversion results to reconstruct the velocity anomaly.

By analyzing the logging data, we set up the referenced velocity model (Table 1) for the tomographic inversion. The block thicknesses on the top and at the bottom of the model are set as 250 m and 500 m, respectively. The horizontal parameterization of the model is restricted by the method of teleseismic tomography. Generally, the horizontal width of the blocks should not be much smaller than the station pace. Otherwise, the inversion results will become unstable. Therefore, the horizontal width of the blocks in the model is set as 2.0 km. In the teleseismic tomographic inversion, we applied a 3×3 offset and average method to smooth the inversions horizontally in the model, as the method proposed by Evans and Achauer (1993) to refine the resultant model.

Table  1.  The reference velocity model used in the ACH inversion

 | Show Table

DownLoad: CSV

The current inversion theory points out that the best result of a inversion is that one which can make the data misfit and the squared model length minimize simultaneously, rather than that one which can only make the data misfit smallest (Yang, 1997). We have to trade-off between minimizing the data misfit and the squared model length. Figure 5 shows the tradeoff curve of the data misfit and the squared model length used in the damped least squares when different damping coefficients are used in the inversion. Based on the tradeoff curve, we choose the optimal damping coefficient θ =8 s², which can reduce the data misfit from 0.011 177 s² to 0.001 922 s² with an improvement rate up to 82.8%.

Figure  5.  The tradeoff curve between the data misfit and the squared length of model perturbation as the damping coefficient θ changes. The numbers in the figure are the damping coefficients θ (unit in s²) used in the least squares inversion.

DownLoad: Full-Size Img PowerPoint

3.   Results and interpretation of the seismic tomograms

By analyzing the results of teleseismic tomography in Figure 7, we may obtain:

Figure  7.  The average seismic tomographic results for 9 inversions with 3×3 offset models where x and y separately denote distance in EW and NS direction. (a)-(f) is the corresponding velocity perturbation at the layers 1, 2, 4, 6, 8 and 10 of the model in different depths (z equal to -0.125, -0.375, -1.25, -2.25, -3.25 and -4.25 km, respectively); the black bold line in sub-figure (a) indicates the fault zone, the dashed line A-B indicates the location of the vertical cross-section (g); sub-figure (g) shows the velocity perturbation of vertical cross-section A-B in (a)

DownLoad: Full-Size Img PowerPoint

1) There is a low velocity anomaly (with a velocity perturbation about 4%) in the southwest of the studied area. The low velocity anomaly is surrounded by high velocity anomalies in the north, northeast and east directions. The intense variation of the near surface velocity is discernible in the shallow horizontal slices with a difference in perturbation up to 16% between the high velocity and the low velocity.

2) The southwest side of the studied area covers the border region between the basin and the mountain, where the Huatugou fault zone crops, and the altitude is relatively lower (Figure 6). The results of the seismic tomographic inversion reveal that the Huatugou fault zone is characterized by a low velocity anomaly, which is approximately marked by a black line in Figure 7a. The low velocity anomaly is characterized by a zone of 10 km×3 km in size with a strike approximate in northwest-southeast orientation. The traveltime residuals (Figure 4) also display they are positive in value in this area. It should just result from the low velocity anomaly.

Figure  6.  Cross-section diagram of Shizigou structure.

DownLoad: Full-Size Img PowerPoint

3) A high velocity anomaly lies in the north, north east and east part of the studied area, encompassing the low velocity zone. From investigation on superficial geology, the Shizigou anticline structure is just located in this region, where the altitude is higher and the topography is terribly undulated. The tomogram slices of our results illustrate the strike and the coverage area of this structure, which is well consistent with the actual geometry of the structure (Figure 1). As shown in Figure 4, the traveltime residuals observed at the stations with higher altitude are negative. Thus, there must be a high velocity geological body in the deep part of this region.

4) Figure 7g is a vertical transect (with its position indicated by a dashed line A-B in Figure 7a) nearly perpendicular to the strike of the anomaly. From this figure, we can obviously see the low velocity anomaly zone with a velocity perturbation up to 4%. This low velocity zone is shown in a high dip angle toward the northeast, which is similar to the Huatugou fault in appearance (Figure 6), indicating that the Huatugou fault must dip northeastward. Besides, the high velocity anomaly on the northeast of the low velocity anomaly becomes smaller in velocity perturbation and extends almost vertically.

5) As shown in Figures 7a-7c and g, there is a trend that the anomalies of velocity perturbation distribute in several different layers. At the depth of z = -0.5 km or around, the low velocity anomaly is divided into two portions by a local high velocity anomaly body, whereas the high velocity anomaly is also divided by a visible low velocity zone within the high velocity anticline structure. This feature should be reliable because it is shown clearly in the tomographic results.

Therefore, the velocity inhomogeneity is displayed in Shizigou region, which is characterized by a low velocity anomaly in the southwest region of the studied area surrounded by high velocity anomalies. The low velocity anomaly dives in a high dipping angle on the southwest side of the high velocity anomaly dipping northeastward.

4.   Discussion and conclusions

In this teleseismic tomography research, we reconstruct the three-dimensional velocity structure beneath the Shizigou in western Qaidam basin. The satisfactory tomographic results reveal the shallow structure and its character in the studied area. Nevertheless, limited by the resolution of the method itself, it is still far from mature to be used in oil and gas exploration. The resolution of the seismic tomography is affected by some factors, such as the distance between stations, the parameterization of initial model and the denseness of seismic rays. In this study, the average of the station-to-station distance is about 1 km. The maximum block size of the model in parameterization is 2 km, which results in a sufficient ray path density for inversion. Thus, the inversion results can well resolve the anomaly bigger than 2 km.

In a word, this research has revealed to some degree the velocity structure of the oil field in the west Qaidam, which provide some information for the deep-seated geological structure of the exploration-difficult area. These results provide us a basis for the research on the deep-seated three-dimensional geological structure. The inversion results obtained show a strong velocity perturbation in the shallow crust with some obviously low and high velocity anomalies distributing in belts, which can be used to outline the anticline and the fault zone.

Therefore, the teleseismic tomography can be used to delineate the strike and the extents of the structures in the studied area. We tested the teleseismic tomography method to resolve small scale structures, which may be beneficial to the practical use in future. The earthquake observation work is still going on in the studied area. We believe that we can improve our results to achieve a better velocity model for this region when more and more seismic data become available soon.