Uncertainty based multi-step seismic analysis for near-surface imaging

1 Introduction

Imaging of near-surface structures is often performed rapidly using a single method. In most commercial applications, this comprises the analysis of surface wave dispersion curves (MASW) [1], or first breaks travel time tomography alone. The inversion of seismic data in the near-surface is a problem and is a non-unique procedure, without other geophysical measurements [2, 3, 4, 5, 6, 7, 8, 9]. To limit this shortcoming, it is possible to use more data or to use interpretation techniques that utilise the whole of the recorded data [10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. These techniques have different resolution capabilities. One interpretation method is to focus on the most precise technique for a given problem. We are convinced that a combination of several techniques will result in better recovery of data for a studied structure, but it is necessary to understand the limitations of the techniques used, especially its uncertainty [20, 21]. Here, using simple field data, we present a complex approach to combine three seismic methods in a strict uncertainty guided procedure, in order to interpret near-surface structure. The results obtained are optimally fitted to the given data and are presented with the uncertainty analysis at each processing step.

2 Motivation

Near-surface geophysical studies have to be simple and cost-effective. In most cases, studied geological structures are rather simple. Our motivation for this work is to develop a simple but effective method that will use shallow seismic data more completely (Figure 1). The application of similar procedures to those used for deep experiments [22, 23] but in a simplified form, would allow us to overcome multiple processing problems and obtain a final result with an estimated uncertainty. The use of additional analyses which are effective in near-surface studies (e.g. MASW), will provide additional information for the correct interpretation of the subsurface.

Figure 1

Example of real data with specified different types of seismic waves. The utilisation of Refracted Waves (A), Surface Waves (B), Reflected Waves (C) and wide-angle Reflected Waves (D) increases the certainty of the final result.

3 Geology

Data acquisition in this case study was conducted at the test site of the Central Geophysical Observatory in Poland. The geology of the region consists of thick Quaternary sediments (definitely thicker than 100 m), which are post-glacial and were created in the Pleistocene. According to the Central Geological Database (http://bazagis.pgi.gov.pl) in the vicinity of the test site, three shallow boreholes were drilled. All of them show quaternary deposits (Mała Wieś borehole 97.8m depth, Belsk 1/6 borehole 15m depth, Pgr Stara Wieś borehole 53m depth). Under the Cainozoic deposits, the Maastrichtian deposits are present [24]. However simple, the geology of the research area is difficult for the near-surface seismic methods because of the high wave attenuation. The well-mixed post-glacial sediments, with small velocity differences, are characterised by the small reflection coefficient differences between the layers. For that reason, the proposed research area is well suited as a test site for the proposed approach to data processing, that will be mainly applied in similar geological conditions. The most interesting geological boundary between Cainozoic and Mesozoic is lying on the depth that exceeds the penetration range of the author’s seismic equipment.

4 Fieldwork

During fieldwork, a 400 m seismic line (Figure 2) covered by 100 stations (60 Data-Cube type (light-blue pointer) and 40 Reftek 130 (orange pointer) 4.5 Hz type receivers (Table 1) was set. The spacing between geophones was 4 m (Figure 2). As a seismic source, the authors used a car-mounted PEG 40 accelerated weight drop in the shooting scheme, assuming signal excitation between geophones (Figure 2 blue stars). Terrain coverage in the fieldwork area changed with distance from the east side of the profile to the forest zone in the west. Because of this, the application of a car-mounted PEG 40 accelerated hammer provided the best possible signal to noise (S/N) ratio for the first 332 m of the profile, while the rest of the shots were executed with a sledgehammer (Figure 2 white stars). At every shot point, a total of six signal excitations was conducted for further vertical stacking. This was an unconventional data acquisition system for a geoengineering purpose (with a single channel, stand-alone stations) and required very careful time measurement for the duration of the experiment. The timing device used was based on a GPS and allowed us to measure each source with a precision better than 2 ms. Because of negligible elevation changes of less than 1m, geodesy measurements were not needed.

Figure 2

The geometry of acquisition profile from project GPB Belsk. The total range of 400 m long line was regularly covered with 100 seismic stations. Sources were regularly spaced at 4 m intervals. After 300m along the profile, the application of a car-mounted, accelerated weight drop was impossible due to rough terrain. The last 72 m of the acquisition line was shot by sledgehammer.

5 MASW data analysis

The application of multi-method, seismic data analysis requires a careful approach to the data quality and an estimation of the uncertainty of the results. Careful quality control was conducted to eliminate errors in the data due to possible timing problems. During that process, some outlier records with time deviations greater than 2 ms, were deleted. The verified data was vertically stacked, with the application of the seismic line geometry. To recover information about the layers nearest to the ground surface, classical MASW analysis was conducted [1, 25, 26, 28]. Dispersion curves for 100 regularly spaced shots were obtained and inverted by applying the genetic algorithm [29] (Eiben et al. 1994). During the MASW processing procedure, the raw data was filtered using the band-pass filter in the range 3 to 60 Hz. Moreover, automated gain control correction and the manual muting of refraction was applied to improve the S/N ratio. For correct tracing of surface waves, 50 channels for every shot point were selected, as an optimal range where energy was visible. 1D V_S models were merged by kriging interpolation and with the assumption that each location in the 2D profile is determined to be in the centre of the spread of the single shot modelling (Figure 3A). Uncertainty estimation (Figure 3B) for the results obtained was based on a statistical comparison between inversion results. However 1D and 2D uncertainties are not the same [30, 31]. In simple cases, a proposed approximation is sufficient. To achieve 2D V_S uncertainty decay in the subsurface model, a similar kriging interpolation method was applied to the results. In the final model (Figure 3A), the authors were able to recognise a dipping layer, with a high V_S velocity contrast. The results obtained by the MASW method were recalculated to V_P by the genetic algorithm (with corresponding uncertainty information) and used in the seismic tomography stage of processing for further analysis [32]. The MASW modelling results provide useful information about low-velocity zones despite their significant uncertainty. Because classical First Arrival Seismic Tomography is blind to the presence of such structures [33], solving the methodological problem is essential to near-surface research. Moreover, the best possible recognition of the shallowest layers is essential for proper static and residual corrections in seismic reflection imaging. In the presented case study, however, MASW analysis proved that there was no low-velocity zone in the research area.

Figure 3

Results of the multichannel analysis of surface waves. To achieve the best possible data utilisation, 101 1D models were created and then interpolated to produce a 2D model (C). For the dataset under consideration: at each shot point, the fundamental and first higher-mode of Rayleigh wave dispersion was visible (A, B). To obtain a 2D model, single 1D models were interpolated. The final 2D model (C) revealed a high-velocity change, visible with a slight dip at 30 m under the surface. Panel (D) presents uncertainty analysis based on the misfit of curves during the inversion process, which revealed a fast degradation of the model precision with depth.

6 Seismic travel-time tomography analysis

Application of the MASW method as the first processing procedure, allowed the authors to receive information about the nearest subsurface. However, because of the high levels of uncertainty in the method, particularly the fast degradation of results with depth, the modelling of layers deeper than 20 to 30 m is ineffective in the presented case. Information about uncertainty in the deeper parts of the subsurface can be treated as boundary conditions in a more accurate method. Such approximations are very useful as a priori information to seismic travel time tomography. The results obtained by this method strongly depend on the starting model. This study presents a novel approach to utilising this information. By assuming that deviations in the results received from MASW limit the intervals of possible velocities, the authors created multiple (40 in total) starting models restricted by the results of the surface analysis method. For seismic tomography modelling, the authors used JIVE3D code [34], based on ray approximation. The starting model, based on MASW results, was further smoothed by including additional information about near-surface velocities from shallow refractions. As a result of seismic tomography modelling, the authors received a smooth velocity image, with the ray coverage shown in Figure 4, and statistically estimated the uncertainty (Figure 5) of the results to be similar to the method used by [35] Melendez et al. in Tomo3D code.

Figure 4

Ray coverage for the last iteration of tomographic inversion. Model is well resolved down to 30 m, and the maximum penetration depth is about 70 m.

Figure 5

(A), Travel time tomography result, obtained with JIVE3D code, with uncertainty analysis (B) based on a modification of initial model introduced in Tomo3D [35]. Obtained results reveal a smooth gradient of velocity with depth and decrease of velocity towards the west. Uncertainty analysis displays high reliability of received model down to 50 m.

In the method proposed by [35] Melendez et al., multiple initial models are created and analysed using basic statistics. In this study, the best initial tomography model (obtained from inversion) was randomly modified in the range of uncertainty from previous MASW analysis. As mentioned above, 40 models were created with different velocities and gradients. For every obtained model, the same inversion path was applied. For all of the results obtained, a mean value (Figure 5A) and standard deviation (Figure 5B) was calculated. The travel-time tomography results obtained show higher accuracy down to a depth of 40 m (Figure 5B) when compared to the MASW method (Figure 3D). Moreover, possible dipping of the layers is also visible in the image (Figure 5A). Information received about horizontal velocities, and the uncertainty of the results was used in the velocity analysis of the seismic imaging, which is the final step of the multi-step seismic analysis. Even though only 70 m of the model were covered by the rays (Figure 4), the deeper velocities were obtained by extrapolation of the velocity field. When compared to seismic tomography, MASW results are more certain between a depth range of 20 to 50 m, with an average difference of 120 m/s. This almost corresponds to a 50% increase in model accuracy for that depth range, assuming that the V_P to V_S ratio is a factor of 1.7. However large this may appear to be, such an approximation for near surface data is effective. Moreover, the combination of results obtained from two different, but connected, physical phenomena increases the reliability of the velocity model. The results obtained from surface wave modelling are more certain in the utmost 10 m, with an estimated uncertainty of 20 m/s of VS, which is about 35 m/s of Vp assuming the Vp/Vs factor of 1.7.

In comparison, the uncertainty of seismic V_P tomography is much larger (150 m/s) for that part of the model. However, when the results are compared directly, the average V_P velocity recalculated from Vs for the first 10m of the model (with a ratio factor of 1.7) was 680 m/s, where the tomographic velocity average from the first breaks for that area was 750 m/s. If we assume an extreme value of V_P/Vs ratio (of, say, 2.0), with an uncertainty of 35 m/s of Vs, then the results are comparable. For that reason, the authors decided to apply tomographic velocity to the reflection imaging.

7 Seismic imaging

There is a difficulty in applying seismic imaging to shallow results obtained from velocity analysis as in many cases no external data is available and so gathering any information about vertical velocities is problematic. The problem results from the presence of strong surface waves and refractions mixing with observed reflections and, as a result, the S/N ratio is low enough to disturb classical, semblance-

based velocity analysis. In the presented work,the authors used a classical image processing scheme [36] (Table 2).

At every step of the processing, multiple parameter settings were tested. During the data quality control step, multiple “death” traces were removed, and the corrected polarity of the traces was recorded by the Reftek 130 Texan stations. Moreover, panel comparison of vertical stacking algorithms was conducted in order to receive the best possible signal correlation. Because there were two types of receiver station, a geometry database was created manually, to obtain residual and static corrections. To improve the S/N ratio, simple noise suppression based on 1D and 2D F-K filters was conducted, with later deconvolution. Instead of applying stand-alone semblance velocity analysis, a smooth tomographic velocity model with estimated uncertainties was used. It is worth noting that the semblance-based analysis applied to the GPB Belsk data-set before adding the information from external methods, was not working properly. In the case of difficult data (e.g. where it was difficult to distinguish between noise and possible reflections), the application of an external velocity field was the only solution in order to obtain useful results. Due to anisotropy, refraction velocities are much faster than stacking ones. In order to fit refraction velocities to the reflection seismic, Dix’ equation [37], and semblance analysis was used for correcting the data fit. As a quality factor of the velocity field, horizontal straightening was observed and analysed. The best results received from the depth-based, estimated velocity field (Table 3) presented a visible NMO correction (Figure 6).

Figure 6

Application of the velocity model received from tomographic modelling to the reflection seismic analysis (NMO), allowing recognition of possible reflectors. Shot gathers in front of the profile (A), and the rear (B) reveal a straightening effect after application of the velocity model, which was impossible to obtain by the classic semblance analysis method.

Next, to the main velocity model, two disturbed velocity fields (faster and slower) were calculated, based on tomographic uncertainty. After applying an identical fitting procedure for all velocity fields, horizontal stacking, post-stack F-K filtering, deconvolution and time to depth conversion were carried out. As a result, the authors received the main seismic stack (Figure 5A) with two other deviation-based results. Comparison between images and manual picking of selected reflecting horizons, allowed us to estimate the difference in depth of the given reflections (Figure 7B). The authors were able to distinguish the main reflectors in the structure and estimate the uncertainty of those horizons in the depth domain (Figure 7A). Application of the tomographic velocity field, not only allowed the authors to solve the problem of correct velocity field estimation, but also resulted in additional information about results uncertainty in the depth domain that could easily be related to the real rock structures. In comparison with the classical approach to velocity analysis for NMO and TTD migration, the approach presented here allowed for fast recognition of the velocity field. However, only extrapolated values from the tomographic field were used for imaging depths greater than 50 m because it was properly recognised that the shallow subsurface allowed for correct correlations of the velocity field with possible reflectors at greater depths.

Figure 7

(A) Seismic stack obtained after data imaging of the GPB Belsk profile. Application of the velocity model received from external methods allowed correct stacking of the data, which is visible for the shallow structures, down to 50 m depth. Since tomographic methods allow for uncertainty estimation, stacks with extreme models allowed us to estimate the uncertainty of recognised reflectors (B).

8 Discussion of the results

In comparison to the standard processing flow that uses only a single data type, the methodology presented in this paper provides a better and more complete result. The main idea of every data assimilation procedure is to utilise every bit of additional information. This in all cases will lead to increases in the final result quality. As proved in the article [22], propagation of uncertainty between methods allows obtaining detailed information about the in-depth structures, even fora difficult sub-salt structures. Even in the industrial seismic profiling, where both high-resolution data, and geophysical information from borehole soundings are available, estimation of the results uncertainty is extremely useful for a correct correlation of the geophysical and geological information. The migration of the information between methods allows tightening the possible space of solutions and limits the interpretation errors. However, in the near-surface experiments, obtaining a certain geological model of the subsurface is more complicated.

In most cases, shallow experiments are conducted in places where no geophysical (e.g. VSP profiling) or geological information from boreholes is available. Even if they exist, in most cases usability of such data is limited to the final geological interpretation of the upper-most part of the seismic image. Moreover, the biggest problem with the estimation of the velocity field for seismic migration is unresolved. Adding the fact, that due to the scale of survey targets, all topographic and residual shifts have a much higher impact on the final result, the situation seems to be unresolvable. The methodology presented in work is a solution to those problems. Because of relatively shallow geological structures, estimation of migration velocity by the refraction methods is effective. To overcome the strong dependence of the final result on the initial model, the MASW and LVL based model was used. In addition to limiting the space of solution, such approach considers nearest surface effects, that are beyond the resolution of travel time tomography. In this paper, the results obtained by the MASW method correlate well with the tomographic solution. As presented, the application of such methodology requires information about result uncertainty, on each step of the processing. That information, not only improve the overall result’s quality but also allow for correct migration of the data between the methods. The seismic stack (Figure 7), where the increase of the horizon uncertainty with depth is visible, explicitly presents the need for such estimations. In the presented case study, the Quaternary post-glacial Pleistocene sediments were the main imaging target. All possible reflections are connected with structures formed during that period. The correlation of our results, with the data from shallow boreholes, especially the deepest ones in “Mała Wieś” and“ Stara Wieś” validated, that the thickness of Quaternary sediments is much higher than 100 meters because there are no strong reflections from larger depths on the seismic stack. Any deeper investigation requires more powerful seismic source and additional acquisition channels for better recording of the post-critical reflections. In the presented methodology, many assumptions which cannot be used for more complex geology was applied. In cases of rough geology, the joint 2D inversion of MASW results should be applied, instead of the interpolation of single 1D models.

Presented case study proves that the applied methodology was effective, as should be easily applied in similar near-surface studies. As a cost and time effective method, presented procedure greatly simplifies and improves the results of near-surface surveys without the need for additional data. In addition to modern high-resolution seismic surveys [39, 40], presented methodology increases the overall quality of the final result. Additionally, the presented methodology can be used in other scales, where layer striping approach is used, for example in industrial scale [41], or even lithospheric research [42]. For further improvement of the results of the presented methodology, the information from additional geophysical methods, based on the different physical parameter (e.g. ERT, GPR) can be used [43]. The scalability of the technique, to other scales, makes it cost and time efficient, and allows to solve multiple processing problems. Additional uncertainty parameter of the final seismic image will greatly simplify the interpretation, which is extremely useful not only in small scale experiments but also industrial-scale projects.

9 Conclusions

This work presents an uncertainty-based, multi-step seismic analysis that allows us to obtain final results with estimated reliability. Direct use of the uncertainty information between the methods allowed the authors to restrict the area of possible solutions for the consecutive methods. That approach leads to final results which showed detailed information about the subsurface through seismic imaging. Moreover, the application of external methods to the reflection seismic, allowed us to find solutions where the classical approach to velocity analysis was ineffective. The case study presented by the authors reveals that propagation of the information about model uncertainty should be treated equally in the results. The assumption that deviations in the results received from previous methods limit the possible solution range of more accurate consecutive ones leads to more certain results. Moreover, that approach allowed simplified processing of the data because the final model from the previous analysis is treated as a starting point for subsequent ones. The presented methodology allows us to overcome the problems that arise in near-surface surveys that are unsolvable using classical processing tools. In the presented study, an example of finding the velocity field for seismic imaging allowed for recovery information about reflection horizons in the dataset. The methodology presented by the authors is correct only if information about uncertainty is estimated and treated with caution. This approach allows correct, partial results transfer between methods,which results in less processing mistakes and a correct final solution. A similar methodology could be used for solving other special geophysical problems.