Abstract
We developed a higher resolution method for the estimation of the three travel-time parameters that are used in the 2D zero-offset, Common-Reflection-Surface stack method. The underlying principle in this method is to replace the coherency measure performed using semblance with that of MUSIC (multiple signal classification) pseudospectrum that utilizes the eigenstructure of the data covariance matrix. The performance of the two parameter estimation techniques (i.e., semblance and MUSIC) was investigated using both synthetic seismic diffraction and reflection data corrupted with white Gaussian noise, as well as a multioffset ground penetrating radar (GPR) field data set. The estimated parameters employing MUSIC were shown to be superior of those from semblance.
1. Introduction
Many important tasks in seismic processing and imaging require the estimation of travel-time parameters. Such parameters include, among others, velocities (e.g., for stacking and time-migration purposes), travel-time slopes and curvatures (e.g., for slant, common-reflection-surface (CRS), multifocus (MF) stacks) and event picking for tomographical methods. As shared with many other areas of activity, a basic feature of seismic signals (referred to as events) is that they exhibit some sort of coherent or aligned energy. More specifically, seismic events (e.g., reflections or diffractions) align themselves along curves or surfaces (referred to as moveouts) within the data. The basic strategy for signal detection and information extraction is to express these moveouts as a function of a few, meaningful parameters and to estimate such parameters so as the moveout optimally approximates the events. In general, the search for parameters, sometimes referred to as wavefront shaping parameters, carry key information about the geological structure under investigation.
To assess how well a moveout, defined by some trial parameters, approximates a target signal, a number of quantifiers (or coherence measures) has been proposed. General discussions on coherency measures applied to seismic data can be found in the pioneering papers of [1–3] with a clear emphasis on the second-order coherence measure semblance. Semblance quantifies the likelihood between the trial moveout and the target event by stacking the data along that moveout and measuring the energy of the output.
Adopting the notation as in [4], for a given sample, , at a given (reference) trace, the so-called semblance coefficient, or simply semblance, , can be mathematically written in the form Here, the semblance coefficient is computed for samples taken from traces in a window centered about the trajectory defined by the moveout equation generated by the trial travel-time parameters (cf. Figure 1). In the following, the given sample, , and reference trace, as well as the number of samples, , and number of traces, , will be fixed throughout. As a consequence, we do not need to incorporate them into the semblance notation, which will be simply written as . To construct the window in Figure 1, proper interpolation is performed to select the appropriate samples. In the language of electrical engineers, the above-described windowing process steers the stacking along the trial moveout.
Semblance can be described in terms of the covariance matrix of the data. Following, for example, [5], within the selected time window along the chosen trial moveout, semblance can be written in the form where is a column vector of ones, which can be referred to as the unitary steering vector, and is the covariance of the data. Assuming that the different sources can be described by a zero-mean stochastic process, the data covariance matrix is given as where is the data matrix, in which is the recorded data at the th trace and th sample. As in usual notation, and represent the expected value and matrix trace, respectively. Moreover, superscripts and represent transpose and conjugate transpose, respectively. As pointed out by [2], Equation (2) provides the interpretation that semblance can be regarded, within the selected time window, as a normalized output/input energy ratio. The denominator, , is the normalization used by semblance in order to generate a maximum peak of unity at the “correct” moveout parameters (namely, the ones for which we have the optimal stack).
Even though semblance is a good measure of coherency, it can in many times provide insufficient resolution for the parameter estimation. That is the case, in particular, for interfering events. There is, thus, a motivation to look for alternatives to overcome these difficulties. Attempts have been made to further improve semblance by using only those parts of the data with higher resolving power [6] and also by introducing weights in the standard semblance formulation [7]. Statistical approaches have also been introduced to increase the resolution of the velocity analysis [8]. In this paper, an alternative to semblance-like techniques will be investigated.
As recognized in sonar and radar applications, methods exploiting the properties of the eigenstructure (namely, eigenvalues and eigenvectors) of the data covariance matrix can lead to far better resolution results than semblance [4, 9, 10]. The basic idea of the eigenstructure approach is to decompose the data covariance matrix into two orthogonal subspaces. The first is the signal subspace, which is generated by the eigenvectors associated to high eigenvalues. The second is the noise subspace, generated by the small or zero eigenvalues. In this paper, we use the eigenstructure method called multiple signal classification (MUSIC), introduced by [9]. MUSIC exploits the fact that the “correct” moveout, represented as a steering vector, must lie in the signal subspace and, therefore, is orthogonal to the noise subspace eigenvectors. As a consequence, the projection of the steering vector onto the noise subspace provides a nearly vanishing value. The inverse of such a projection (namely, the sum of the dot products of the steering vector with the noise eigenvectors) should peak when the steering vector represents a correct moveout.
This work can be seen as a followup of [10], in which the application of MUSIC to single-parameter velocity analysis and slant stacks is described. Here, we extend the application of MUSIC to common-reflection-surface (CRS) multi-parameter estimation. Besides the theoretical exposition of the technique, applications to first synthetic examples, consisting of dipping planar reflectors and point diffractors, are provided. Comparisons of the obtained results and conventional semblance confirm, at least for these initial examples, the expected far better resolution of MUSIC. To further support this conclusion a real multioffset GPR data set was also analysed. It was demonstrated that MUSIC, unlike semblance, was able to better resolve interfering events.
2. Classical Music: Narrowband and Uncorrelated Signals
In its original or classical form [9], MUSIC considers an array of receivers recording incoming reflected or diffracted signals, in an arbitrary background medium. In time domain, the data recorded by the th receiver can be modeled as where is the source pulse associated with event , and is the additive random noise at the th receiver. Finally, is the travel-time (or time delay) of the th incoming signal (or event) arriving at the th receiver. The superscript, , indicates that the moveouts depend on a set of one or more parameters, here denoted, by a so-called parameter vector, . The most popular trial-moveout example is the normal-moveout (NMO), applied for velocity analysis in the common-midpoint (CMP) configuration. In the 2D situation, the single parameter to be estimated is the NMO-velocity. An example of multiparameter moveout is the general hyperbolic moveout used by the common-reflection-surface (CRS) stacking method. As previously indicated, application of MUSIC to velocity analysis has been described by [10]. Here, we extend the analysis to CRS parameter estimation in 2D data. In this situation, three parameters are to be estimated. In order not to disturb the main flow, the description of the generalized hyperbolic or, more simply, the CRS travel-times, , is postponed to the appendix.
2.1. Narrowband Signals
For narrowband signals , the travel-times can be expressed as exponential phase shifts around the center angular frequency . For notation simplicity, that fixed frequency will be omitted. As a consequence, the data model in (4) can be recast as After time discretization, the above equation can be recast in matrix form as where and are, respectively, the data and additive noise matrices, and is the source matrix. Finally, is the array response matrix containing all the steering vectors MUSIC utilizes the eigenstructure of the data covariance matrix defined by (3). Substituting (6) into (3) and assuming uncorrelated noise with variance of , the covariance matrix can be recast as where and are, respectively, the source covariance and identity matrices. The MUSIC algorithm performs an eigen-decomposition of this covariance matrix where contains the eigenvalues satisfying , and is the matrix that consists of the corresponding (column) orthonormal eigenvectors of . The unitary matrix of eigenvectors can be decomposed further as , where the columns of comprise the eigenvectors corresponding to the largest eigenvalues of (the signal subspace), and with containing the remaining (noise) eigenvectors.
2.2. Uncorrelated Signals
For MUSIC to be applicable in our parameter search problem, the different source pulses, , should be uncorrelated resulting in a covariance matrix having full rank equal to the number of events recorded at the receivers. If the source vectors are linearly independent, then the matrix is positive definite which results in to be a positive semidefinite matrix with its rank spanning the steering vectors corresponding to the appropriate parameters we are searching. With the above condition satisfied and since the noise subspace is orthogonal to the signal subspace, the MUSIC pseudospectrum, , is given by where is the test steering vector and is the noise subspace projection matrix given by . Since the steering vectors are orthogonal to the eigenvectors spanning the noise subspace , it follows that the parameter estimates will occur at those parameter values for which we have This corresponds to large peaks in the MUSIC pseudospectrum as given by (11).
2.3. Wideband Uncorrelated Signals
As indicated above, the MUSIC algorithm was originally developed for narrowband and uncorrelated signal applications. If the condition of uncorrelated signals is maintained, an alternative to this situation is to decompose a wideband data into narrowband data components and then treat each narrowband separately [10]. The MUSIC pseudospectrum at the center angular frequency of the th narrowband now takes the form where and are respectively the test steering vector and the noise subspace projection matrix at the th center angular frequency . The strategy followed in this work is to Fourier transform the test data and select a narrowband close to the center frequency of the source pulse as input to MUSIC.
3. Seismic Music: Wideband and Correlated Signals
Seismic signals are highly correlated and require a special modification to be used by the original MUSIC algorithm. The consequence of having correlated sources is that there will be a rank deficiency in the source covariance matrix that will result in a mix of signal and noise subspaces. As a result, the MUSIC algorithm will loose its power to peak at the “right” parameters.
In order to handle correlated sources, spatial smoothing over the covariance matrix, can be employed [10]. The idea is to subdivide the array of sensors into identical overlapping subarrays of receivers (cf. Figure 2) and then compute the covariance for all the subarrays and average the result. If the covariance matrix for subarray is , the spatially smoothed covariance is given by
To be able to implement spatial smoothing within seismics, one has to taper the data within a window following the event(s) (cf. Figure 1). The purpose of this tapering is to make the delay times of the event linear (which is the basic requirement behind spatial smoothing) [10].
The other advantage of performing the analysis in a given window is to make the steering vectors, required for generating the MUSIC pseudospectrum, to be frequency independent. This allows us to handle wideband seismic data. This process of windowing the event can also be interpreted as steering of the correlation matrix before eigendecomposition and using unity steering vectors for generating the MUSIC pseudospectrum [4].
Ideally, when the window is “perfectly” matching the event, which will be the case of an optimal choice of the moveout parameters, the signal would be flattened and all traces will nearly have the same moveout. As a consequence, the steering vectors used in (11) will be simply replaced by a vector of ones making them frequency independent. In this situation, the MUSIC pseudospectrum generates a peak resulting in the identification of the optimal estimates of the parameters.
In practice, the windows are constructed by moveouts, defined by trial parameters. Peaking of the corresponding MUSIC pseudospectra identifies, thus, the “correct” parameters. Following this approach, [4] has shown that MUSIC can be applied for the single-parameter case of velocity analysis. The objective was, thus, to obtain a high-resolution velocity spectrum. In this work, we extend that strategy to the CRS multi-parameter estimating problem. In other words, our objective is to obtain high-resolution estimates of the CRS parameters, which are three in the present 2D situation.
4. Numerical Examples
In this section, we compare MUSIC and semblance for travel-time parameter estimation in the situations of classical MUSIC (narrowband uncorrelated signals) and seismic MUSIC (wideband correlated signals). For a simple model of a point diffractor and a dipping reflector with a homogeneous overburden, we analyzed the cases: (a) CMP configuration, which requires the determination of a single parameter, of CRS travel-time (cf. (A.1)) and (b) ZO configuration, which requires the determination of two parameters, namely, parameters and of the CRS travel-time (cf. (A.1)). The two events (diffraction and reflection) were chosen to be almost undistinguishable. All the test parameter points with RMS velocities within m/s (step size 7.5 m/s) were tested for parameter and points within (step size ) for parameter and points within (step size ) for parameter . As seen below, in all situations, MUSIC performed much better then semblance.
4.1. Classical Music
To illustrate the application of MUSIC for narrowband uncorrelated signals, we considered a point diffractor and a dipping reflector illuminated under a CMP configuration. For a given CMP gather, the data consists of (compare with (5)) where and are the sources and and are the travel-times for the diffractor and dipping reflector events, respectively. Moreover, is the additive noise. The travel-times for these two events are described by the ordinary NMO equations where the velocity coefficients for the diffractor, , and dipping reflector, , are given by Here, and represent the root mean square (RMS) velocity and the dip angle, respectively.
The sources, and are produced by a single narrowband source, , modified by two realizations of a random phase perturbation, and , so as to produce uncorrelated sources. In frequency domain, this process is generally described as A synthetic CMP gather was generated employing (15), (16), and (18) together with a Ricker zero-phase wavelet with a center frequency of 20 Hz (cf. Figure 3). The fold was 40 representing a half-offset range from 40 m to 820 m. The data was sampled with 2 ms and white Gaussian noise with a variance of of the maximum trace amplitude was added. The parameter estimation process was benchmarked using the classical semblance analysis of [2].
The output from MUSIC (cf. (13)) is shown in Figure 5 together with the result obtained using semblance. For both cases we used, a window size of 11 samples and in addition for MUSIC we considered two signal subspaces and the rest as noise subspaces. As a result, MUSIC is seen to outperform semblance and resolve the two parameters well. It is well known that the values output from MUSIC are arbitrary. To avoid this phenomenon, we have introduced a semblance balancing. This technique is discussed in detail in connection with the real-data example presented below.
To perform a two-parameter test, we have now simulated a zero-offset (ZO) section for the same previous point diffractor and dipping reflector (cf. Figure 4). The corresponding two ZO travel-times for diffraction and reflection are now given by with and (corresponding to a dip of and a homogeneous medium with constant velocity 2000 m/s). As seen from the Appendix, the above equations represent the generalized hyperbolic (CRS) travel-time of equation (A.1), in which the conditions have been implemented. As indicated in the Appendix, the far-left equation above represents the diffraction condition. The far-right equation is due to the fact that in this considered experiment, the N-wave is planar.
Based on (19), using the previous uncorrelated sources (18), synthetic ZO data were computed for midpoints between 40 m and 1040 m. The results from the two-parameter search ( and ) are shown for, respectively, MUSIC (cf. Figure 6) and semblance (cf. Figure 7). MUSIC gives well-resolved results, as opposed to the semblance, where the estimated parameters are more inaccurate.
4.2. Seismic Music
To examine the performance of MUSIC compared to semblance in case of wideband correlated signals, we generated synthetic data based on the travel-time (16) for a CMP gather and (19) for a ZO section. The parameter search was performed within a time window of 25 time samples with the ZO travel-time being the middle sample and following a hyperbolic delay trajectory defined by the travel-times. For the computation of the MUSIC pseudospectrum, the samples within the hyperbolic window were used to form the data covariance matrix and the associated eigendecomposition. In order to reduce the correlated source effect we performed spatial smoothing of the covariance matrix using 31 subarrays each consisting of 10 receivers for the CMP data and 37 subarrays each consisting of 15 receivers for the ZO data.
The results of the parameter search is shown in Figure 8 for the CMP data (i.e., determination of parameter ) and Figures 9 and 10 for the ZO section (i.e., determination of parameters and ). It is apparent that both semblance and MUSIC can resolve parameters and , but MUSIC shows a higher resolution in general. Moreover, for parameter only MUSIC is able to resolve the two events.
4.3. Real Data Example Using GPR Data
The first step of the CRS analysis determining the parameter can be regarded as a CMP-based velocity analysis. As indicated by our previous synthetic data example, MUSIC was seen to have a better potential than semblance for resolving interfering events (cf. Figures 5 and 8). We will now investigate whether that feature is confirmed in real data. Prior to our analysis, however, the following normalization issue has to be considered. As opposed to semblance, which produces normalized values between 0 and 1, MUSIC, despite its high-resolution capability, yields arbitrary amplitude values. Such behavior makes the simple replacement of semblance with MUSIC as a coherency measure, for example, in standard velocity analysis, not adequate.
In order to condition MUSIC to be a normalized quantity, we introduce a scaled version of it, denoted by semblance-balanced MUSIC or, more simply, SB-MUSIC. In the framework of velocity analysis, SB-MUSIC is defined as follows: for a given CMP location, as well as a selection of zero-offset time samples, and trial stacking velocities, , we let and represent the coherency values obtained from MUSIC and semblance, respectively. In other words, and represent velocity spectra associated with MUSIC and semblance coherency measures. Denoted by , SB-MUSIC is given by where Application of the above conditioning makes sure that those amplitude anomalies inherent to the original MUSIC velocity spectrum are balanced according to the energy level of semblance.
A real multioffset GPR data set was used to test out the feasibility of this approach. For an in depth description and discussion of these data, the readers are referred to [11]. Figure 11 shows an example of a typical velocity spectrum obtained from the GPR data using both semblance and SB-MUSIC. In these computations, we used a window size of eleven samples for both semblance and MUSIC. In addition, we performed spatial smoothing of subarray size 15 from a fold of 28 to ensure that MUSIC handles the correlated GPR signals properly. Figure 11 clearly demonstrates that interfering events are much better resolved in the SB-MUSIC spectrum (Figure 11(b)) than in its corresponding semblance spectrum (Figure 11(a)). In particular, as indicated by white arrows, it can be seen how two interfering events are unresolved by semblance (Figure 11(a)) and well resolved by SB-MUSIC (Figure 11(b)). To further validate the previous observation, the hyperbolic moveout curves corresponding to those two events were superimposed to the corresponding CMP-gather (cf. Figure 12(b)). These curves seem to correlate well with two interfering events. As a reference, the result obtained using semblance is also included (cf. Figure 12(a)). It can be regarded as a fit based on a mix between the two interfering events.
(a)
(b)
(a)
(b)
5. Conclusions
In this paper, we discussed the CRS travel-time parameters estimation problem in seismic signal processing. The conventional semblance algorithm was found to generate lower-resolution estimates of the parameters. For the purpose of obtaining higher-resolution parameter estimates, we replaced semblance with MUSIC algorithm. Such procedure allowed us to estimate the parameters within a resolution limit that is significantly better. This work can be seen as a followup of previous applications of MUSIC to single-parameter velocity analysis and slant stacks. Now, MUSIC has been extended to Common-Reflection-Surface (CRS) multiparameter estimation. Applications of the technique to first synthetic examples, consisting of dipping planar reflectors and point diffractors, and comparison to semblance, confirm, at least for these initial situations, the expected far better resolution of MUSIC. To further support this analysis, CMP velocity analysis has been applied to a real multioffset GPR data set. In this situation, better results were obtained upon the introduction of a scaled version of MUSIC, denoted semblance-balanced MUSIC. The new algorithm was seen to outperform semblance in resolving interfering events.
Appendix
General Hyperbolic Moveout
The CRS method uses the so-called generalized hyperbolic (normal) moveout, which is the natural generalization of the NMO, valid for CMP gathers, to CRS supergathers, in which source-receiver pairs are arbitrarily located around the (reference) central point, usually taken as a CMP. In 2D, the generalized hyperbolic moveout depends on three parameters, as opposed to conventional NMO, which depends on a single parameter (NMO velocity).
Mathematically, the generalized hyperbolic moveout, , associated with the event, , measured at receiver , is specified by the zero-offset (ZO) travel-time, , and (reference) trace location, , and given by (see Figure 13) where is the midpoint coordinate and is the half-offset coordinate for the th receiver. Here, is the CRS parameter vector, with three parameters, , and , to be estimated from the data. It is instructive to recall that these parameters are related to the angle and curvature quantities as follows [12]: where and are the curvatures of respectively the normal () and normal-incident-point (NIP) wavefronts, is the emergence angle and is the medium velocity. All these quantities are evaluated at the central point, . Still considering the CRS parameters, we make the following observations(a)In the CMP configuration of source-receiver pairs symmetrically located with respect to the central point, namely, , we have with the CMP, single parameter vector . Moreover, we have the relation with given by the lower-most (A.3).(b)In case the recorded data stems from a diffraction, the condition holds. This is because as the reflector shrinks to a point, the N-wave turns out to be identical to the NIP-wave [13]. As a consequence, the hyperbolic moveout of diffraction (or diffraction travel-time), reduces to with the diffraction, two-parameter vector .
Acknowledgments
The authors would like to thank Dr. Hervé Perroud for providing the GPR dataset. E. Asgedom has been funded by a PhD grant from the University of Oslo and the Norwegian Science Foundation. This work has been carried out partly while he was visiting State University of Campinas. M. Tygel acknowledges support of the Brazilian Council of Scientific and Technological Development (CNPq) and the sponsors of the Wave Inversion Technology (WIT) Consortium.