Statistical approaches to interpretation of seismic reflection data

Tectonophysics 329 (2000) 251±267

www.elsevier.com/locate/tecto

Statistical approaches to interpretation of seismic re¯ection data C.A. Hurich*, A. Kocurko Earth Science Department, Memorial University, St. John, Newfoundland, Canada A1B 3X5 Received 12 May 1999; accepted 2 December 1999

Abstract Statistical approaches offer a potentially powerful tool for interpretation of the complex seismic re¯ection patterns typically recorded from the crystalline crust. These methods focus on characterization of the complexity, scaling characteristics and spatial variability of the wave ®eld with the ultimate aim of recovering information on the geometric variation of the geology. In this paper we examine the potential and limitations of statistical approaches, focusing on the issues of quality of the statistical estimator, how noise and migration affect spatial statistics and the potential for both qualitative and quantitative approaches to statistics-based interpretation. The limited spatial bandwidth of typical seismic data presently inhibits complete quantitative recovery of the spatial properties of the geology. Migration, which should increase the spatial bandwidth, has to date, not proven particularly effective on complex wave ®elds. Mapping of relative variation in the spatial properties of the wave ®eld reveals subtle differences that are directly associated with the geometric character of the geology. Such mapping allows identi®cation and mapping of variations in the geologic fabric at scales signi®cantly larger than the seismic wavelength and objective access to scaling information that has previously been unavailable. q 2000 Elsevier Science B.V. All rights reserved. Keywords: re¯ection seismology; seismic interpretation; crustal seismology; fractal analysis; scaling

1. Introduction The geology of the crystalline crust is commonly characterized by a geometrically complex distribution of lithology and associated acoustic impedance resulting from multiple intrusion and deformation events. When illuminated by a seismic wave, the complex acoustic impedance ®eld results in complicated re¯ection patterns that are dif®cult to interpret using conventional interpretation procedures. In the absence of speci®c tectonic features, it is often only possible to generally categorize a volume of crust as re¯ective or non-re¯ective and perhaps to derive general infer* Corresponding author. Tel.: 11-709-737-2384; fax: 11-709737-2589. E-mail address: [email protected] (C.A. Hurich).

ences on tectonic fabric based on regions of common dip. Alternative methods of treating deep seismic data that provide descriptions of the complexity and spatial variability of the re¯ection wave ®eld offer the potential for extracting more of the information contained in crustal seismic data. An example of such an approach is presented by Vasudevan and Cook (1998). Fractal analysis of geologic maps, representative of various crustal environments, provides information on scaling characteristics and is a fruitful approach to describing the geometric complexity of the geology (Holliger and Levander, 1992, 1994a,b; Holliger et al., 1993; Levander et al., 1994a,b). Coupling velocity and density information to the geometric description provides a statistical description of the acoustic impedance ®eld. Theoretical arguments (Gibson, 1991;

0040-1951/00/$ - see front matter q 2000 Elsevier Science B.V. All rights reserved. PII: S 0040-195 1(00)00198-0

252

C.A. Hurich, A. Kocurko / Tectonophysics 329 (2000) 251±267

Pullammanappallil et al., 1997) and experimental results (Hurich, 1996) indicate that the spatial statistics of the acoustic impedance ®eld and the spatial statistics of the associated re¯ection wave ®eld are correlated. Correlation between the spatial statistics of the geology and the re¯ection wave ®eld suggests that determination of the spatial statistics of the re¯ection wave ®eld should lead to at least partial recovery of the spatial statistics of the geology. Mapping of variations in spatial statistics could provide another tool to aid interpretation of crustal seismic data. In this paper we discuss the potential of statistical mapping as a tool for enhanced interpretation of seismic data. We use both synthetic and real seismic data to illustrate both the value and limitations of statistical approaches. In particular, we focus on the issues of quality of the statistical estimator, how noise and migration affect spatial statistics and the potential for both qualitative and quantitative approaches to statistics-based interpretation. Our discussion is restricted to spatial (lateral) variations in the seismic wave ®eld due to the added complications of source signature removal and severe band-width restrictions in the temporal domain. 2. Approaches to determining spatial statistics There are several methods available for determining the spatial statistics of seismograms including analysis of the 2-D autocorrelation matrix or power spectrum (Holliger and Levander, 1994b; Hurich, 1996; Pullammanappallil et al., 1997), the analysis of variograms (Rea and Knight, 1998) and description of families of events derived from skeletonization through sample moments (Vasudevan and Cook., 1998). The autocorrelation, variogram and power spectrum operate directly on the wave ®eld and are mathematically related providing essentially the same information. Because these approaches operate directly on the wave ®eld they nominally incorporate noise as well as signal into the statistics but they retain amplitude and characteristic frequency information. Operating on skeletonized data requires reliance on a skeletonization algorithm to faithfully reproduce the wave ®eld. However, skeletonized data has the advantage of allowing analysis on subsets of the

event population and the selective elimination of noise prior to analysis. Experimental comparison between an autocorrelation-based approach and a skeletonization-based approach (Hurich, 1996) suggests that the two methods measure different but related aspects of the wave ®eld but any such comparison is dependent on the details of the skeletonization procedure. All of the statistical analyses discussed in this paper are derived from analysis of the 2-D autocorrelation. Determination of spatial statistics based on the autocorrelation depends upon comparison of the autocorrelation of the sample to the autocorrelation of a model with appropriate statistical attributes. The von KaÂrmaÂn model has been used by a number of workers for describing spatial distributions in geological problems because it describes processes that are self-similar or self-af®ne and the Fourier transform, autocorrelation and power spectrum can be written conveniently (Frankel and Clayton, 1986). The 1-D von KaÂrmaÂn model can be described by three parameters, the correlation length (L) and the Hurst number (H) and the variance. For scale sizes smaller than L, the von KaÂrmaÂn model describes a power-law (fractal) process where H represents the power exponent. The fractal dimension (F) is related to H by the relation, F ˆ …E 1 1† 2 H where E is the Euclidian dimension. For scale sizes longer than L, the model represents a process that varies as white noise. For the analyses presented in this paper, we derive the autocorrelation matrix through Fourier transform-based estimation of the 2-D power spectrum and inverse transform of the power spectrum to obtain the autocorrelation. The observed autocorrelation is ®t to the von KaÂrmaÂn model through a least-squares optimization procedure to determine L and H (Pullammanappallil et al., 1997). In addition to the two von KaÂrmaÂn parameters, we use the autocorrelation to estimate the mean power, coherence and slowness (apparent dip) of the dominant energy in the analysis window. 3. Correlation between spatial statistics of geology and the seismic wave ®eld Interaction between the re¯ection wave ®eld and the geology involves a complex combination of

C.A. Hurich, A. Kocurko / Tectonophysics 329 (2000) 251±267

253

Fig. 1. Estimations of lateral correlation length (L) and Hurst number (H) for a suite of acoustic impedance models and the associated synthetic seismograms. The estimation space is the entire 20 £ 5 km impedance model and 20 km £ 1.8 s synthetic seismograms. (a) Estimated and theoretical L for the suite of acoustic impedance models. Solid line is the best linear ®t to the data, dashed line is the ®t for a perfect estimation. (b) Estimated H, theoretical H and theoretical L for the acoustic impedance models. (c) Estimated and theoretical L for the suite of synthetic seismograms. (d) Estimated H for the synthetic seismograms.

wave propagation phenomena including specular re¯ection, primary and higher order scattering and spatial and temporal interference. To use the spatial properties of the re¯ection wave ®eld to infer the spatial properties of the geology it must ®rst be established that the re¯ection wave ®eld retains relevant spatial information. It is also necessary to assess how the partitioning of the various modes of wave propagation affect the spatial properties of the wave ®eld. To this purpose we review the relevant theoretical and experimental arguments that demonstrate a correlation between the spatial properties of the acoustic impedance ®eld and the associated re¯ection wave ®eld. These arguments form the underpinning for statistical approaches to seismic interpretation.

3.1. Theoretical Pullammanappallil et al. (1997) relate the statistics of the re¯ection wave ®eld to that of a stochastic medium based on the approximate relations between the primary re¯ectivity series (PRS) and actual seismic data. The PRS is a conceptual model representing the ideal (perfectly imaged) seismic response. They demonstrate that, for both the t±x and t±k domains, normalized representations of the horizontal (spatial) characteristics of the stochastic medium and the associated PRS are equivalent. Assuming that seismic data recorded in the ®eld is a good representation of the PRS, the spatial statistics of the ®eld data should represent the spatial properties of the geology. The

254

C.A. Hurich, A. Kocurko / Tectonophysics 329 (2000) 251±267

validity of this assumption depends upon how well the real seismic data estimates the PRS in the presence of noise, multiple scattering and imperfect migration. 3.2. Experimental Hurich (1996) describes a set of numerical experiments in which the lateral correlation lengths of ensembles of self-similar acoustic impedance models are compared to the measured correlation lengths of the associated synthetic seismograms. The synthetic seismograms used in this analysis were generated using a ®nite difference solution to the acoustic wave equation that includes the effects of spatial interference and multiple scattering. These experiments demonstrate clear correlation between the characteristic L of the geologic models and the associated re¯ection wave ®elds. However, in these experiments, the estimator of L was chosen as the 50% down point on the autocorrelation and L derived from the wave ®elds is systematically smaller than the theoretical L of the model impedance ®elds. These experiments also demonstrate that the L of the wave ®eld depends upon the relative contributions of weak and strong scattering. Re¯ection wave ®elds dominated by strong scattering are characterized by a signi®cant decrease in L. In general, the relative contribution of weak and strong scattering is a function of the magnitude of impedance contrasts (re¯ection coef®cients) and the characteristic frequency of the seismic data. Consideration of the range of re¯ection coef®cients expected in the continental crust (Hurich and Smithson, 1987) and the frequencies typically recorded in deep seismic data suggests that weak scattering is dominant in the crystalline crust but in areas typi®ed by large re¯ection coef®cients ….0:1† strong scattering can play a signi®cant role. Because we do not have a priori knowledge of the scattering regime for a particular crustal volume and we cannot determine it directly from the seismic data, an ambiguity in inverting the spatial properties of the re¯ection wave ®eld for the spatial properties of the geology arises. Variation in the spatial properties of the seismic wave ®eld may be indicative of change in the dominant mode of scattering rather than variation in the spatial distribution of the acoustic impedance contrasts. The results of a suite of numerical experiments

similar to those reported by Hurich (1996) except that L and H are determined by ®tting the von KaÂrmaÂn model are shown in Fig. 1a±d. In this set of experiments the acoustic impedance models are self-similar in the spatial dimensions. The horizontal correlation length (Lx) varies between 750 and 5000 m, the vertical correlation length (Lz) is constant at 100 m, the velocity ®eld is binary varying between 6.0 and 6.4 km/s, and H of the discrete ®eld is constant at 0.25 …F ˆ 2:75† (see discussion of parameterization of discrete and continuous von KaÂrmaÂn ®elds by Goff et al., 1994). The synthetic seismograms are determined using a ®nite-difference solution to the acoustic wave equation based on Kelly et al. (1976) with a frequency band between 2 and 35 Hz and a grid spacing of 25 m. Wave propagation occurs almost exclusively in the weak scattering regime. An example of one realization of the acoustic impedance ®eld and the resulting seismogram are shown in Fig. 2. Fig. 1a and b compare the estimated and theoretical values of L and H for the entire suite of acoustic impedance models. The results indicate that the estimator recovers the spatial statistics of the impedance ®eld but there is a slight tendency to underestimate L particularly for longer values of theoretical L. This underestimation of L must be, at least in part, the result of the limited size of the model space and points to the need for care when choosing sample sizes for analysis (Eneva, 1996). Estimation of H is relatively good (mean H ˆ 0:247; s ˆ 0:028† and notably the estimated H shows no dependence on L. This suite of experiments validates the von KaÂrmaÂn correlation function as an estimator for the spatial properties of the broadband acoustic impedance model. Fig. 1c and d shows the estimated and theoretical values of L and H for the synthetic seismograms associated with the suite of acoustic impedance models. Consistent with the experiments reported by Hurich (1996) there is a clear correlation …r2 ˆ 0:94† between L of the seismograms and the acoustic impedance ®elds but L of the seismic data is consistently signi®cantly less than that of the impedance ®eld (Fig. 1c). Hurich (1996) attributed the difference in L of the seismograms and the impedance models to the interaction of the various wave phenomena that contribute to the re¯ection wave ®eld. However, Bean et al. (1999) have recently suggested that the limited bandwidth of the seismic data also leads to a systematic

C.A. Hurich, A. Kocurko / Tectonophysics 329 (2000) 251±267

underestimation of L, a suggestion with which we concur. H of the seismograms signi®cantly overestimates the H of the impedance models (Fig. 1d). Overestimation of H implies that the wave ®eld is not as rough as the impedance model, a result consistent with limited temporal and spatial bandwidth and wavefront healing processes. The estimated H of the seismograms also shows a dependence on the theoretical L. For impedance models with theoretical L signi®cantly less than the Fresnel zone diameter (2500 m), H approaches unity indicating that the lateral scaling in the wave ®eld is essentially mono-

255

tonic, the wave ®eld is not self-af®ne. For impedance models with theoretical L larger than the Fresnel zone, although H is consistently overestimated, it is independent of L. 3.3. Summary Experimental results suggest that the re¯ection wave ®eld retains only a partial record of the spatial statistics of the impedance ®eld. Underestimation of L and overestimation of H based on the seismic data is a result of the band-limited nature of the seismic data and the interaction between various wave phenomena

20 km

5 km

a)

Vp = 6.0 km/s

Vp = 6.4 km/s

Hurst # - 0.25 Horizontal Correlation Length - 2000 m Vertical Correlation Length - 100 m

Reflection Time (s)

b)

1.7

2.2

2.7

3.2 Fresnel Zone Diameter

Fig. 2. Acoustic impedance model and associated synthetic seismogram. The seismogram is generated by a ®nite-difference solution for the acoustic wave equation and uses a plane wave to estimate the CMP stacked seismic section.

256

C.A. Hurich, A. Kocurko / Tectonophysics 329 (2000) 251±267

that contribute to the re¯ection wave ®eld. The difference in the theoretical and experimental results appears to be primarily due to the difference between the PRS estimate of the wave ®eld which results in a spatially broadband seismogram and the unmigrated spatially band-limited seismograms used for the experimental study. 4. Utility and limitations of statistical approaches to seismic interpretation The utility of statistical approaches to seismic interpretation depends upon whether or not new information can be derived from the data and whether or not the new information provides insight into interpretation of the data. Aside from the issue of overcoming the limitations imposed by the limited bandwidth of seismic data and extracting quantitative estimates of the spatial properties of the associated geology, potential information also exists in mapping relative variations in the spatial statistics. Such a mapping is accomplished by applying the estimator to an appropriately sized sliding window of the seismic data (Vasudevan and Cook, 1998). In this section, we assess the potential for deriving new information from variations in the spatial statistics of seismic data by examining synthetic and real data examples for which the geology is well constrained. 4.1. Synthetic example The geologic model and associated synthetic seismogram for one realization of the suite of numerical experiments summarized in Fig. 1 are shown in Fig. 2. The model is 20 km long and includes a 5-km-thick heterogeneous zone overlain by a 5 km homogeneous layer. The heterogeneous zone is de®ned by a binary velocity function that varies between 6.0 and 6.4 km/s. The spatial distribution of the two lithologies is self-similar with H ˆ 0:25; L ˆ 2000 m and a vertical correlation length of 100 m. The synthetic seismogram contains frequencies in the band 2±35 Hz resulting in an average Fresnel zone diameter (Fd) of 2.5 km in the heterogeneous zone. Wave propagation in this example occurs dominantly in the weak scattering regime so the wave ®eld is dominated by specular re¯ections and small body

Table 1 Summary of estimates of the spatial statistics of the entire acoustic impedance model and various modi®cations of the synthetic seismogram Experiment

Correlation length (L) (m)

Hurst #(H)

Impedance model Synthetic seismogram Add noise Model noise Migration

1535 644 892 694 595

0.267 0.594 0.267 0.512 0.655

diffractions modi®ed by spatial and temporal interference (Hurich, 1996). Direct comparison between the synthetic seismogram, the statistical maps and geologic model provides an opportunity to evaluate the signi®cance of variations in the spatial statistics as well as identify the contribution of the various factors that make up the wave ®eld. Fig. 3 shows a mapping of L and H of the impedance model for the synthetic example. Results of analysis of the entire model space are shown in Table 1 for comparison. Analysis of the entire model correctly recovers the spatial statistics of the impedance ®eld within the variability expected for the individual realizations. Although the resolution of the statistical maps is limited by the size of the analysis window …x ˆ 5000; z ˆ 500 m†, local and rather subtle variations in the geometric character of the impedance model are clearly mapped in both L and H space implying signi®cant heterogeneity in the local spatial properties of the impedance model. The mapping of L (Fig. 3a and b) emphasizes local extremes in the correlation length that vary at scales signi®cantly larger than the seismic wavelength. Larger L values tend to be associated with regions of higher aspect ratio (x/z). The mapping of H (Fig. 3c and d) also indicates local variability but the range of variation is considerably less than L. Even though the variability of H is signi®cantly less than that of L, the two parameters are anticorrelated with regions of low L tending to be correlated with zones of higher H and vice versa. This anticorrelation between L and H is not consistent with the analyses of the entire model space shown in Fig. 1a and b but is consistent with the analyses of the associated suite of synthetic seismograms shown in Fig. 1c and d.

C.A. Hurich, A. Kocurko / Tectonophysics 329 (2000) 251±267 257

Fig. 3. Mapping of the (a and b) lateral correlation length (Lx) and (c and d) Hurst number (H) of the acoustic impedance model with and without the impedance model superimposed. The dimensions of the sliding analysis window …x ˆ 5000 m; y ˆ 500 m† are shown by the rectangle in the upper left corner of the impedance model. The window shifted by 1/4 of the horizontal or vertical dimension for each analysis in the map.

258

C.A. Hurich, A. Kocurko / Tectonophysics 329 (2000) 251±267

4.1.1. Selectivity of the seismic re¯ection system The seismic re¯ection technique preferentially images boundaries with high aspect ratios (length to thickness). The continuity of a seismic event is arguably the most important feature with regard to detecting or recognizing the event as signal, particularly if amplitude information is lost or modi®ed during processing. From this point of view, the geometric component of the impedance ®eld (geology) is at least as important as the magnitude of the re¯ection coef®cients for producing detectable re¯ections. Due to the contribution of the geometric component, it is likely that geology with a complex geometric distribution will have some regions that are favored for imaging over others. Such selective imaging, if it occurs, may weigh estimates of the spatial statistics for large regions of a seismogram toward the preferentially imaged geometries and away from less preferred geometries. On the other hand, understanding the origin of the selectivity and mapping of regions with preferred and less preferred geometries may provide signi®cant information on the fabric of the geology. To understand better the selectivity of the re¯ection seismic method and the potential for statistical bias, we map the variation in the mean power in the synthetic seismogram estimated by the zero-lag of the autocorrelation prior to normalization (Fig. 4a). The power map is superimposed onto the original impedance model in Fig. 4b. As expected, zones of high power in the seismogram are associated with the strongest and most continuous re¯ections, the events most likely to be detected in the presence of noise. It is more interesting however to identify the class of geometries that produces the high power re¯ections. Fig. 4b demonstrates that the high power re¯ections originate in zones typi®ed by relatively thin, continuous (high aspect ratio) impedance contrasts. In general, these re¯ective zones represent regions of the most intermixing of the contrasting impedances. The zones of lower power are dominantly associated with regions of lower aspect ratio and generally homogeneous impedance distribution. The power mapping demonstrates that the seismic re¯ection system is inherently selective of the geometries that are imaged. Re¯ection data from areas of complex geology image a skeleton of preferred re¯ectivity that is controlled as much or more by the local

geometry as by the magnitude of the impedance contrasts. In terms of interpretation, characterization of the scale and distribution of the holes in the re¯ective skeleton may provide another signi®cant source of geological information of equal importance as characterizing the re¯ective zones. The apparent selectivity of the seismic system casts some confusion on what is actually measured by the spatial statistics of a seismogram. Are the statistics dominated by the skeleton of re¯ections, the holes in the skeleton or some combination of both? We defer discussion of this important question until we present more analyses that help illuminate the problem. 4.1.2. Spatial statistics of the synthetic seismogram The synthetic seismogram shown in Fig. 2 contains a few approximately specular re¯ections 1±2 Fd in length that are associated with a few relatively smooth and continuous geologic boundaries. The majority of events are less than 1 Fd in length and are of two types. The ®rst type has amplitudes similar to the specular re¯ections and originates from a combination of small body (Fresnel) diffraction and constructive temporal interference (layer tuning). The second type of short event is of lower amplitude and dominantly represents constructive interference between the numerous low amplitude diffraction tails. Fig. 5 shows a mapping of the variation in L and H for the synthetic seismogram. The results of the whole window analysis are included in Table 1. We make three signi®cant observations from this mapping: (1) the L and H maps are correlated; (2) regions of homogeneous impedance that map with long L in the impedance ®eld, map with short L in the seismogram; and (3) L is sensitive to variations in the local dip of the seismic events. Correlation between L and H maps in the signal-only seismogram (long L mapping to small H and the inverse) is consistent with the experimental results presented in Fig. 1d. The correlation re¯ects the limited spatial bandwidth of the seismogram for L less than the Fresnel zone and the increased bandwidth for L greater than the Fresnel zone. Regions of relatively homogeneous impedance and low aspect ratio that are characterized by long L and low H in the maps of the impedance ®eld are represented by short L and high H in the seismogram. This difference re¯ects the predominance of short incoherent events associated with the homogeneous

C.A. Hurich, A. Kocurko / Tectonophysics 329 (2000) 251±267

259

Fig. 4. (a) Mapping of the average power for the synthetic seismogram. (b) Average power map superimposed on the acoustic impedance model. The power mapping demonstrates that the strongest (most detectable) re¯ections originate in regions of the impedance model characterized by high aspect ratio heterogeneities.

regions. Fluctuation of re¯ections about the horizontal and continuous local dip both result in low L estimates, but the geologic signi®cance of the two situations is quite different. Fluctuations about the horizontal re¯ect some form of discontinuity in the geology while continuous local dips generally re¯ect

true geologic dip. Continuous dipping events result in a skew of the long axis of the autocorrelation function away from the spatial axis. Determination of L along the long axis provides a better estimate of L. However, because the autocorrelation exists in x±t space, determination of L on any axis other than the spatial axis

260

C.A. Hurich, A. Kocurko / Tectonophysics 329 (2000) 251±267

Fig. 5. Mapping of the lateral correlation length (L) and Hurst number (H) for the synthetic seismogram. The dimensions of the sliding analysis window …x ˆ 5000 m; z ˆ 160 ms† are shown by the rectangle in the upper left corner of the impedance model. The window shifted by 1/4 of the horizontal or vertical dimension for each analysis in the map.

results in an estimate that does not have pure spatial dimensions. Our experience suggests that for moderate dips and qualitative interpretation purposes, the added discrimination between re¯ection discontinuity and geologic dip outweighs the problem of mixed dimensions.

The statistical mapping of the synthetic seismogram represents our best band-limited estimate of the statistics of the impedance model. It is clear that signi®cant differences associated with the selectivity and limited bandwidth of the seismic method make direct inversion of the seismic data for the spatial

C.A. Hurich, A. Kocurko / Tectonophysics 329 (2000) 251±267

statistics of the impedance ®eld using a von KaÂrmaÂn process based estimator dif®cult. 4.1.3. The effects of noise on spatial statistics Noise almost always represents a signi®cant component of the wave ®eld in crustal seismic data. The noise component certainly affects the spatial statistics of the wave ®eld and if the noise and signal cannot be adequately separated, noise represents a problem for extraction of both qualitative and quantitative information from the wave ®eld. Statistical analysis has the advantage that if the noise can be characterized statistically it can be included directly into the analysis. To investigate the effects of random noise on the statistics of the wave ®eld and the effectiveness of modeling the noise as part of the statistical analysis, we added band-limited Gaussian noise to the synthetic seismogram and remapped the variation in L and H. The noise population is designed so that a signal that is 40% of the maximum signal power is well de®ned in the seismic section. The results of mapping the noisy data without attempting to model the noise population are shown in Fig. 6 and the whole seismogram analysis is summarized in Table 1. The addition of noise signi®cantly lowers H in both the windowed and whole seismogram analyses although the pattern of relative variations in the windowed mapping remains essentially the same. A decrease in H suggesting an increase in the apparent roughness of the wave ®eld with the addition of noise is readily predictable and physically reasonable (compare the wave ®elds in Figs. 5 and 6). Interestingly, L increases signi®cantly in both the whole seismogram and windowed analysis, a result that at ®rst view is counter intuitive. However, we suggest that the decrease in H and increase in L result from the apparent broadening of the spatial bandwidth of the seismic data due to the addition of noise with a short correlation length. The increased spatial bandwidth results in a better ®t between the correlation function of the seismogram and the broad-band von KaÂrmaÂn model (as determined by least-squares error) and a shift in the model parameters. This result demonstrates a weakness in applying the von KaÂrmaÂn model to the band-limited seismogram and emphasizes the need to consider the noise population as part of the statistical analysis. The results of including the noise in the statistical analysis is shown in Fig. 7 and summarized in Table 1.

261

The noise was modeled as a von KaÂrmaÂn process with L ˆ 50 m and H ˆ 0:1 with the noise representing 25% of the total power in each window. The properties of the noise were determined by inspection of the autocorrelation for the whole seismogram. Comparison between Figs. 5 and 7 and the whole seismogram analyses summarized in Table 1 demonstrate that inclusion of the noise in the statistical analysis results in reasonably good recovery of the spatial statistics of the signal-only wave ®eld. Modeling noise with small characteristic L is an effective and generally robust procedure because the signal-only wave ®eld is de®cient in the short correlation lengths due to the limited spatial bandwidth of the seismic system. Statistical characterization of the short L noise population in real seismic data is more dif®cult than our well controlled synthetic example because many seismic processing procedures tend to increase the characteristic L and the intersection between noise and signal populations may be dif®cult to objectively determine. Objective discrimination between the small L noise and re¯ection signal populations is the target of ongoing research. It is clear that our approach to small L noise is not effective for noise with characteristic L of the same order as the re¯ection signal. 4.1.4. Statistical analysis of unmigrated and migrated data In the preceding discussion we have argued that the limited spatial bandwidth of unmigrated seismic data imposes limitations on the spatial information that can be derived for the impedance ®eld. Migration is the obvious approach to increasing the spatial bandwidth, if the migration algorithm is successful in properly migrating complex wave ®elds with modest signal-tonoise ratios. Traditional migration schemes do not properly handle multiply scattered events so the presence of signi®cant multiply scattered energy further complicates the problem. Pullammanappallil et al. (1997) compared L and H estimates for CMP stacked and migrated synthetic seismograms and reported only minor differences in estimates for the two data domains. They choose to carry out analyses on real data in the common shot and CMP stack domains. Our experiences with statistical analyses on unmigrated and migrated data are similar to Pullammanappallil et al. (1997). As an example, a mapping of L and H for a migrated version of the synthetic seismogram

262

C.A. Hurich, A. Kocurko / Tectonophysics 329 (2000) 251±267

Fig. 6. Mapping of the lateral correlation length (L) and Hurst number (H) for the synthetic seismogram in the presence of band-limited Gaussian noise. The analysis window is the same as in Fig. 5. The presence of noise increases the estimate of L and decreases the estimate of H by arti®cially broadening the spatial bandwidth of the seismogram.

is shown in Fig. 8 and the whole seismogram analysis is summarized in Table 1. The migration was carried out with a phase-shift algorithm and 1-D linear velocity function. The results of both the mapping and the whole window analysis are similar to the simulated stack for both L and H.

Comparisons between analyses of unmigrated and migrated seismic data indicate migration even of noiseless synthetic data do not signi®cantly improve estimation of spatial statistics. It is permissible and indeed may be preferable to carry out statistical analysis of real seismic data in the shot or stack domain.

C.A. Hurich, A. Kocurko / Tectonophysics 329 (2000) 251±267

263

Fig. 7. Mapping of the lateral correlation length (L) and Hurst number (H) for the synthetic seismogram by modeling the noise as a distinct subpopulation of the wave. The analysis window is the same as in Fig. 5. Modeling the noise effectively recovers the spatial statistics of the original seismogram.

The limited quality of a migrated image in the presence of modest signal-to-noise ratios and multiple scattering may degrade rather than enhance statistical analyses. It is clear however that in terms of the proper placement of dipping re¯ections in seismic data, migration is required.

4.2. A real data example Maps of the spatial statistics of seismic data over the Bjerkreim-Sokndal layered intrusion in southern Norway (Deemer and Hurich, 1997) are shown in Fig. 9. The data provide an example against which we can

264

C.A. Hurich, A. Kocurko / Tectonophysics 329 (2000) 251±267

Fig. 8. Mapping of the lateral correlation length (L) and Hurst number (H) for a migrated version of the synthetic seismogram. The analysis window is the same as in Fig. 5.

test some of the inferences derived from the synthetic experiments. These data also provide a particularly useful example of the utility of statistical mapping as an aid to seismic interpretation because downplunge projection of surface geology allows con®dent identi®cation of the origin of the many re¯ections and comparison between re¯ection patterns and the spatial properties.

In general, the L map (Fig. 9c) indicates uniformly small correlation lengths with the exception of a few distinct zones with strong, coherent re¯ections and longer correlation lengths. The Fresnel zone diameter varies from about 300 to 1500 m from the top to the bottom of the seismic section. Assuming that statics corrections are reasonably complete, we infer that with the exception of the few zones of specular or

C.A. Hurich, A. Kocurko / Tectonophysics 329 (2000) 251±267 Fig. 9. (a) Conventional interpretation of migrated and coherency ®ltered seismic data across the Bjerkreim-Sokndal layered intrusion based on down-plunge projection. (b) Mapping of the Hurst number (H) and (c) lateral correlation length (L) of the unmigrated version of the seismic data. The spatial statistics, particularly the Hurst number, provide additional constraints on the actual location of the base of the intrusion, the signi®cance of chaotic re¯ections within an apparently homogeneous leuconorite (MCU1A) and the location and geometry of the subvertical contact between underlying anorthosites and migmatitic country rock.

265

266

C.A. Hurich, A. Kocurko / Tectonophysics 329 (2000) 251±267

near-specular re¯ections, the majority of the re¯ected wave ®eld must be associated with sub-Fresnel zone scaled geologic variations. The mapping of H (Fig. 9b) shows considerably more variation than L and is more revealing about variations in the scattering dominated wave ®eld. Visual comparison between the seismic data and the H mapping indicates that low (but not the lowest) H values are associated with the least coherent (roughest) zones of the wave ®eld and increasing H values are associated with increasingly coherent (smoother) regions in the wave ®eld (for example compare zones labeled A and C in Fig. 9b). From a scaling point of view, low H values imply similar power for all scale lengths less than L and higher values of H imply that power is partitioned preferentially into some scale lengths over others. We suggest that these variations in H re¯ect variations in the organization and lithologic mix of the geology at a sub-Fresnel zone scale. The zones of specular or near specular re¯ection and highest L values have the lowest H values mapped. This observation is consistent with the observed anticorrelation between L and H in the synthetic experiments but is apparently inconsistent with the suggested correlation between H and wave ®eld smoothness for scattering-dominated wave ®elds. This inconsistency suggests a fundamental difference in the response of the broad band von KaÂrmaÂn based estimator for wave ®elds dominated by scattering and specular re¯ection. This inconsistency is the topic of ongoing research. The majority of the coherent re¯ections in the Bjerkreim-Sokndal data mark the boundaries of megacyclic cumulate sequences that represent new in¯uxes of magma into the chamber and a variation in pyroxene abundance. Three important interpretation issues remain after correlation with the downplunge projection: (1) the exact location of the base of the intrusion; (2) the signi®cance of relatively short chaotic re¯ections from within large zones of apparently homogeneous leuconorite (eg. megacycle MCU1A); and (3) the location of the contact between the anorthosites which crop out west of the intrusion and most likely underlie part of the intrusion and the migmatitic country rock that outcrops east of the intrusion. Addressing these questions requires an objective assessment of subtle variations in the re¯ection wave ®eld.

The Bjerkreim-Sokndal intrusion represents a latestage intrusive body emplaced along the eastern edge of the much larger Egersund Anorthosite complex. The base of the intrusion is marked by a transition between anorthosite and leuconorite but is not marked by a speci®c distinctive set of re¯ections. The base of the intrusion as shown in Fig. 9a, a migrated and coherency ®ltered version of the seismic data, is generalized from a down-plunge projection and the synformal nature of the intrusion. Mapping the variation in H on the unmigrated seismic data (Fig. 9b) con®rms a general and distinctive change in the geometric character of the re¯ection wave ®eld that correlates with the projected base of the intrusion. This correlation con®rms the geologic signi®cance of a rather subtle change in the wave ®eld. The correlation between variation in H and the base of the intrusion also con®rms the signi®cance of the short chaotic re¯ections within the seemingly massive leuconorite by establishing that they have geologic signi®cance and are not simply noise or processing artifacts. The chaotic re¯ections most likely represent isolated anorthosite pods within the leuconorite (J. Scoates, pers. commun.). Beneath the intrusion, the zone of lower H in the western 3/4 of the pro®le (A in Fig. 9) probably represents anorthosites of the Egersund suite that outcrop west of the intrusion but not to the east. The higher H zone (C) in the eastern 1/3 of the pro®le likely represents the migmatitic country rock that outcrops to the east of the Bjerkreim-Sokndal intrusion. The variation in H between the anorthosites and the migmatites probably is consistent with an increase in the organization of moderate aspect ratio lithologic variations. 5. Discussion Statistical approaches offer an alternative method of seismic interpretation that focuses on the complexity and spatial variability of the wave ®eld. The results of our experimental work and recent work reported by Bean et al. (1999) suggest that the re¯ection wave ®eld retain only a partial record of the spatial statistics of the impedance ®eld. The limited spatial bandwidth of the seismic data results in underestimation of the L and overestimation of H of the impedance ®eld and distinct limitations for directly deriving the spatial

C.A. Hurich, A. Kocurko / Tectonophysics 329 (2000) 251±267

statistics of the impedance ®eld from the seismic data. Our conclusions however differ from other workers (eg. Pullammanappallil et al., 1997) that report more successful recovery of the statistics of the impedance ®eld. Recovery of the spatial statistics of the impedance ®eld remains a major objective of our research in this area because knowledge of the spatial statistics opens up the entire realm of forward stochastic modeling. We suggest that adaptive revision of the von KaÂrmaÂn model to match the bandwidth of the seismic data and insure proper recovery of the wave ®eld statistics followed by extrapolation based on the wide-band von KaÂrmaÂn model may be a fruitful approach to establishing the parameters for stochastic modeling. Despite the dif®culties of completely recovering the spatial properties of the impedance ®eld, mapping of relative variations in the spatial properties of the re¯ection wave ®eld represents a potentially powerful tool for interpretation. Mapping of subtle variations in the wave ®eld allows identi®cation and mapping of variations in the geologic fabric at scales signi®cantly larger than the seismic wavelength as well as identi®cation of subtle variations in scattering associated with sub-Fresnel zone geologic variations. An excellent example of the utility of this approach is the application of statistical mapping to aid resolution of outstanding interpretation questions in the Bjerkriem-Sokndal data. The mapping of the synthetic example also clearly identi®es regions of more or less homogeneity in an otherwise uninterpretably complex seismogram. Preliminary work in this area also suggests that multiple scale domains are present in some seismic data suggesting overprinting of scale domains in the geology. Acknowledgements We thank S. Deemer, S. Pullammanappallil, K. Vasudevan and an anonymous reviewer for their considerable aid in clarifying portions of this manuscript. This research is supported by an NSERC individual research grant to CAH. References Bean, C.J., Marsan, D., Martini, F., 1999. Statistical measures of crustal heterogeneity from seismic data: the role of seismic bandwidth. Geophys. Res. Lett. 26, 3241±3244.

267

Deemer, S.J., Hurich, C.A., 1997. A seismic image of the basal portion of the Bjerkreim±Sokndal layered intrusion. Geology 25, 1107±1110. Eneva, M., 1996. Effect of limited data sets in evaluating the scaling properties of spatially distributed data: an example from mining-induced seismic activity. Geophys. J. Int. 124, 773±786. Frankel, A., Clayton, R.W., 1986. Finite-difference simulations of seismic scattering: implications for propagation of short-period seismic waves in the crust and models of crustal heterogeneity. J. Geophys. Res. 91, 6465±6489. Gibson, B.S., 1991. Analysis of lateral coherency in wide-angle seismic images of heterogeneous targets. J. Geophys. Res. 96, 10,261±10,273. Goff, J.A., Holliger, K., Levander, A., 1994. Modal ®elds: a new method for characterization of random seismic velocity heterogeneity. Geophys. Res. Lett. 21, 493±496. Holliger, K., Levander, A.R., 1992. A stochastic view of lower crustal fabric based on evidence from the Ivrea Zone. Geophys. Res. Lett. 19, 1153±1156. Holliger, K., Levander, A.R., 1994a. Structure and seismic response of extended continental crust: stochastic analysis of the StronaCeneri and Ivrea zones, Italy. Geology 22, 79±82. Holliger, K., Levander, A.R., 1994b. Seismic structure of gneissic/ granitic upper crust: geological and petrophysical evidence from the Strona-Ceneri zone (northern Italy) and implications for crustal seismic exploration. Geophys. J. Int. 119, 497±510. Holliger, K., Levander, A.R., Goff, J.A., 1993. Stochastic modeling of the re¯ective lower crust: petrophysical and geological evidence from the Ivrea zone (northern Italy). J. Geophys. Res. 98, 11,967±11,980. Hurich, C.A., 1996. Statistical description of seismic re¯ection wave®elds: a step towards quantitative interpretation of deep seismic re¯ection pro®les. Geophys. J. Int. 125, 719±728. Hurich, C.A., Smithson, S.B., 1987. Compositional variation and the origin of deep crustal re¯ections. Earth Planet. Sci. Lett. 85, 416±426. Kelly, K.R., Ward, R.W., Treitel, S., Alford, R.M., 1976. Synthetic seismograms: a ®nite-difference approach. Geophysics 41, 2±27. Levander, A.R., Hobbs, R.W., Smith, S.K., England, R.W., Snyder, D.B., Holliger, K., 1994a. The crust as a heterogeneous `optical' medium or `crocodiles in the mist'. Tectonophysics 232, 281±297. Levander, A.R., England, R.W., Smith, S.K., Hobbs, R.W., Goff, J.A., Holliger, K., 1994b. Stochastic characterization and seismic response of upper and middle crustal rocks based on the Lewisian gneiss complex, Scotland. Geophys. J. Int. 119, 243±259. Pullammanappallil, S., Levander, A., Larkin, S., 1997. Estimation of crustal stochastic parameters from seismic exploration data. J. Geophys. Res. 102, 15,269±15,286. Rea, J., Knight, R., 1998. Geostatistical analysis of ground-penetrating radar data: a means of describing spatial variation in the subsurface. Water Resour. Res. 34, 329±339. Vasudevan, K., Cook, F.A., 1998. Skeletons and fractals Ð a statistical approach to deep crustal seismic data processing and interpretation. Tectonophysics 286, 93±109.