Effects of Random Heterogeneity on Seismic Reflection Images

The reflection seismic method, which is a technique to map geologic structure and stratigraphic features, has been adopted in a variety of applications such as oil and gas explorations, fundamental geological studies, engineering and hydrological studies. Furthermore, recently time-lapse seismic monitoring is considered as a promising technology to monitor changes in dynamic physical properties as a function of time by analyzing differences between seismic data sets from different epochs (e.g., Lumley, 2001). On the other hand, its widespread use has often revealed a weakness in seismic reflection methods when applied to complex structures. It is widely believed that highly dense spatial sampling increases the quality of final seismic reflection sections. However, the quality of final seismic sections obtained in real fields is often very poor for a variety of reasons such as ambient noise, heterogeneities in the rocks, surface waves, reverberations of direct waves within the near-surface, and seismic scattering, even if highly dense spatial sampling is adopted. In most of reflection seismic explorations, people implicitly assume that the subsurface target heterogeneities are sufficiently large and strong that other background heterogeneities only cause small fluctuations to the signals from the target heterogeneity. In this case, a clear distinction can be made between target structures and the small-scale background heterogeneities. However, if the small-scale heterogeneities are significantly strong and are of comparable size to the seismic wavelength, complicated waveforms often appear. This complication causes much difficulty when investigating subsurface structures by seismic reflection. In deep crustal studies (Brown et al., 1983) or geothermal studies (Matsushima et al., 2003), seismic data often have a poor signal-to-noise ratio. Complicated seismic waves are due to seismic wave scattering generated from the small-scale heterogeneities, which degrades seismic reflection data, resulting in attenuation and travel time fluctuations of reflected waves, and the masking of reflected waves by multiple scattering events. In this case, the conventional single-scattering assumption of migration may not be applicable; in other words, multiple scattering caused by strong heterogeneities may disturb the energy distribution in observed seismic traces (Emmerich et al., 1993).

174 seem incoherent and the small-scale heterogeneity is presumed to be randomly distributed, the statistical properties of seismic wave fluctuation relate to the statistical properties of this small-scale heterogeneity. Seismologists conclude that coda waves are one of the most convincing pieces of evidence for the presence of random heterogeneities in the Earth's interior. Seismic evidence suggests random heterogeneity on a scale ranging from tens of meters to tens of kilometers. In addition, geologic studies of exposed deep crustal rocks indicate petrologic variations in the lithosphere on a scale of meters to kilometers (Karson et al., 1984;Holliger and Levander, 1992). Well-logging data suggest that small-scale heterogeneities have a continuous spectrum (Shiomi et al., 1997).
From the viewpoint of seismic data processing, many authors have pointed out the disadvantages of the conventional CMP method proposed by Mayne (1962) when applied to complex structures. Based on a layered media assumption, the CMP stacking method does not provide adequate resolution for non-layered media. Since the 1970s, several prestack migration methods have been studied as improvements on CMP stacking. Sattlegger and Stiller (1974) described a method of prestack migration and demonstrated its advantages over poststack migration in complex areas. Prestack migration is divided into two types of techniques: prestack time migration (PSTM) and prestack depth migration (PSDM). PSTM is acceptable for imaging mild lateral velocity variations, while PSDM is required for imaging strong lateral velocity variations such as salt diapirism or overthrusting. A better image is obtained by PSDM when an accurate estimate of the velocity model exists; however, the advantage of PSTM is that it is robust and much faster than PSDM. From the viewpoint of the S/ N ratio, Matsushima et al. (2003) discussed the advantages of prestack migration over synthetic data containing random noise.
Wave phenomena in heterogeneous media are important for seismic data processing but have not been well recognized and investigated in the field of seismic exploration. There are only several studies which have taken into account the effect of scattering in the seismic reflection data processing. Numerical studies by Gibson and Levander (1988) indicate that different types of scattered noise can have different effects on the appearance of the final processed section. Gibson and Levander (1990) showed the apparent layering in CMP sections of heterogeneous targets. Emmerich et al. (1993) also concluded that the highly detailed interpretation, which is popular in crustal reflection seismology, is less reliable than believed, as far as the internal structure of scattering zones and scatterer orientations are concerned. Sick et al. (2003) proposed a method that compensates for the scattering attenuation effects from random isomorphic heterogeneities to obtain a more reliable estimation of reflection coefficients for AVO/ AVA analysis. It is important to understand how scattered waves caused by random heterogeneities affect data processing in seismic reflection studies and how these effects are compensated for. From the viewpoint of spatial sampling in time-lapse seismic survey, Matsushima and Nishizawa (2010a) reveal the effects of scattered waves on subsurface monitoring by using a numerical simulation of the seismic wave field and comparing the different responses of the final section by applying two different types of data processing: conventional CMP stacking and poststack migration. Matsushima and Nishizawa (2010a) demonstrate the existence of a small but significant difference by differentiating two sections with different spatial sampling. This small difference is attributed to the truncation artifact which is due to geometrical limitation and that cannot be practically prevented during data acquisition. Furthermore, Matsushima and Nishizawa (2010b) indicate that this small difference is also attributed to normal moveout (NMO)-stretch effect which cannot be practically prevented during data acquisition and processing.
A primary concern of this article is to study effects of random heterogeneity on seismic reflection images. We investigate the effect of spatial sampling on the images of seismic reflection, by comparing two set of images: one reproduced from simulated seismic data having a superimposed random noise in time series, and the other generated from numerically simulated wave fields in a same medium but containing random heterogeneity. We also investigate the relationship between the spatial sampling interval and the characteristic size of heterogeneities and also investigate from the viewpoint of spatial sampling how noise-like scattered wave fields that are produced from random isotropic heterogeneity influence the seismic section. We consider the adoption of highly dense spatial sampling with intervals smaller than the Nyquist interval to improve the final quality of a section. In this paper, three types of data processing, conventional CMP stacking, poststack migration and prestack migration are compared to examine different responses to the migration effect of different spatial sampling intervals. We generate 2-D finite-difference synthetic seismic data as input to this study. Our numerical models have a horizontal layered structure, upon which randomly distributed heterogeneities are imposed.

Spatial sampling interval in seismic reflection
According to the Nyquist sampling theorem, sampling at two points per wavelength is the minimum requirement for sampling seismic data over the time and space domains; that is, the sampling interval in each domain must be equal to or above twice the highest frequency/ wavenumber of the continuous seismic signal being discretized. The phenomenon that occurs as a result of undersampling is known as aliasing. Aliasing occurs when recorded seismic data violate the criterion expressed in equation (1).
where x Δ is the spatial sampling interval which should be equal to or smaller than the spatial Nyquist sampling intervals is the minimum velocity, max f is the maximum frequency, and θ is the dip angle of the incident plane-wave direction.
On the other hand, in the case of zero-offset, the spatial sample interval should be equal to or smaller than a quarter-wavelength (Grasmueck et al., 2005). Aliasing occurs when recorded seismic data violate the criterion expressed in equation (2).
In the presence of structural dips or significant lateral velocity variations, adequate sampling becomes important for both vertical and lateral resolution. For the case of the maximum dip (θ=90), the spatial Nyquist sampling interval becomes a quarter-wavelength. Thus, quarter-wavelength spatial sampling is a minimum requirement for adequate recording. Vermeer (1990) defined the term " full-resolution recording" for unaliased shooting and recording of the seismic wave field at the basic signal-sampling interval. In practice, however, seismic data are often irregularly and/ or sparsely sampled in the space domain because of limitations such as those resulting from difficult topography or a lack of resources. In many cases, proper sampling is outright impossible. In order to avoid aliasing, standard seismic imaging methods discard some of the high frequency components of recorded signals. Valuable image resolution will be lost through processing seismic data (Biondi, 2001). Once seismic data are recorded, it is difficult to suppress aliasing artifacts without resurveying at a finer spatial sampling (Spitz, 1991).
In the case of migration processing, there are three types of aliasing (Biondi, 2001), associated with data, operator, and image spacing. Data space aliasing is the aliasing described above. Operator aliasing, which is common in Kirchhoff migration algorithms, occurs when the migration operator summation trajectory is too steep for a given input seismic trace spacing and frequency content. Kirchhoff migration approximates an integral with a summation and is subject to migration operator aliasing when trace spacings do not support the dip of the migration operator. In contrast, migration algorithms such as the f-k method or finite-difference methods only require that the input data volume be sampled well enough to avoid aliasing of the input volume (Abma et al., 1999). Adequate solution for operator aliasing is to control the frequency content (e.g., low-pass filtering at steep dips).
The anti-aliasing constraints to avoid operator aliasing can be easily derived from the Nyquist sampling theorem. The resulting anti-aliasing constraints are (Biondi, 1998): where data x Δ is the sampling rate of the data x-axis and op p is the operator dip.
Image space aliasing occurs when the spatial sampling of the image is too coarse to adequately represent the steeply dipping reflectors that the imaging operator attempts to build during the imaging process. Image space aliasing can be avoided simply by narrowing the image interval. But for a given spatial sampling of the image, to avoid image space aliasing we need to control the frequency content of the image. Similarly to the case of operator aliasing, the anti-aliasing constraints to avoid image space aliasing can be easily derived from the Nyquist sampling theorem. The resulting anti-aliasing constraints are (Biondi, 1998): where image x Δ is the image sampling rate of the x-axis and ref p is the reflector dip.
From the viewpoint of the S/ N ratio, dense spatial sampling increases the number of sources/ receiver pairs (i.e., stacking fold), which raises the effect of signal enhancement, that is, increases the S/ N ratio. The expected improvement in S/ N is proportional to the square root of the stacking fold under the assumption that it is purely random noise which has a flat power spectrum. Thus, highly dense spatial sampling improves the S/ N ratio of the section, even if the interval of spatial sampling becomes shorter than the Nyquist sampling interval.

Construction of synthetic data and seismic reflection imaging
We constructed two data sets. One is synthetic seismic data set generated from two-layer model where each layer has a constant velocity everywhere inside the layer. Random noise was added to the synthetic seismic data (random noise model=RN model). The other is synthetic seismic data set generated from two-dimensional random heterogeneous media where random velocity variation is superimposed on a layer above a reflector (random heterogeneous model=RH model). The second model will generate incoherent events by scattering of waves in the random heterogeneous media.

Random noise (RN) model
A numerical simulation model and source/ receiver arrangements are shown in Figure 1a. A reflector is placed at a depth of 2000 m, separating two layers having a constant velocity of 3800 m/ s and 4200 m/ s, respectively. Three different source-receiver intervals 80, 20, and 5 m were employed; each requiring 26, 101, and 401 sources and receivers, respectively. The reflected waves generated by a flat reflector were obtained by using the 2-D finite difference method as described below. In order to remove direct wavelets, the wavefield without the reflector was subtracted from the total wavefiled of the reflector model. We then obtain the wavefield containing only reflected waves. Random noise is added to the data containing only signal components (reflections) so that the S/ N ratio was 0.3. The S/ N ratio is defined as the following equation (5): where MAX S is the absolute value of the maximum amplitude of signal events in a stacked trace obtained from data consisting of only signal components, Noise (i) is the amplitude of the i-th sample in a stacked trace obtained from the random noise, and N is the total number of samples. The denominator of equation (5) equals the root-mean-square (rms) amplitude of the noise.

Layered model overlapped with random heterogeneity
Random heterogeneous media are generally described by fluctuations of wave velocity and density, superposed on a homogeneous background. Their properties are given by an autocorrelation function parameterized by the correlation lengths and the standard deviation of the fluctuation. Random media with spatial variations of seismic velocity were generated by the same method as described in Frankel and Clayton (1986). The outline of the scheme is as follows: 1. Assign a velocity value v(x, z) to each grid point using a random number generator. 2. Fourier transform the velocity map into the wave number space. 3. Apply the desired filter in the wavenumber domain. The first two-layered random media model for two-dimensional acoustic wave simulation using the finite-difference method. The average velocity of the upper layer is 3800 m/ s with 3% standard deviation and correlation distance 10 m. (c) The second two-layered random media model with the same average velocity and standard deviation as for (b), except for a correlation distance of 50 m In this paper, the applied filter (Fourier transform of autocorrelation function, which is equal to the power spectral density function) has a von Karman probability distribution described by equation (6) where k is the wavenumber, β is the Hurst number that controls the components of small scale random heterogeneities, and a is the correlation distance indicating the characteristic heterogeneity size. The wavenumber k we use here is defined by equation (7): where λ is the wavelength. We use the above von Karman-type heterogeneous media with β =0.1. Saito et al. (2003) described that the value β =0.1 is nearly the same as the value for the power spectral density function of velocity fluctuation obtained from well-log data at depths shallower than 10 km (e.g., Shiomi et al., 1997;Goff and Holliger, 1999).
A homogeneous model and source/ receiver arrangements are the same as the case of the RN model. To estimate the relationship between the spatial sampling interval and the characteristic size of heterogeneities, two types of random heterogeneities were generated and implemented in the layered model as shown in Figures 1b and 1c. The velocity perturbations shown in Figure 1b were normalized to have a standard deviation 3% of the 3800 m/ s (upper) and 4200 m/ s (lower) layers on average and a characteristic heterogeneity size of 10 m (a=10 m). Figure 1c is the same as Figure 1b except for characteristic heterogeneity sizes of 50 m (a=50 m).
The level of scattering phenomena is a function of the wavelength and the average scale of heterogeneities. If the wavelength of a seismic wave is much longer than the scale length of heterogeneity, the system is considered a homogeneous material. Although scattering phenomenon is important only at wavelengths comparable to the scale length of heterogeneity, small-scale heterogeneities influence the seismic waveform with respect to the size of heterogeneities. Wu and Aki (1988) categorized the scattering phenomena into several domains. When ka < 0.01 (Quasi-homogeneous regime), the heterogeneous medium behaves like an effective homogeneous medium where scattering effects may be neglected. When 0.01 < ka < 0.1 (Rayleigh scattering regime), scattering effects may be characterized by Born approximation which is based on the single scattering assumption. When 0.1 < ka < 10 (Mie scattering regime), the sizes of the heterogeneities are comparable to the wavelength. The scattering effects are most significant. When ka > 10 (Forward scattering regime), the heterogeneous medium may be treated as a piecewise homogeneous medium where ray theory may be applicable.

Wave field calculation
We employed a second-order finite difference scheme for the constant density twodimensional acoustic wave equation described in the equation (8).
where P is the pressure in a medium and V(x,z) is the velocity as a function of x and z. The source wavelet was the Ricker wavelet with a dominant frequency of 20 Hz. The dominant frequency (20 Hz) and the average velocity (3800 m/ s) yielded the dominant wavelength (190 m). A uniform grid was employed in the x-z plane. To minimize grid dispersion in finite difference modeling, the grid size was set to be about one eighteenth of the shortest wavelength, which was calculated from the minimum velocity of 3600 m/ s, the maximum frequency of around 40 Hz ( max f =40), and a 5-m grid spacing. All edges of the finitedifference grid were set to be far from source/ receiver locations so that unnecessary events would not disturb the synthetic data. Source/ receivers were not located on the edge of the model, but within the model body. In this situation, scattered wave fields generated in the heterogeneous media above the source/ receiver locations would be included in the synthetic data. However, this does not affect the conclusions of this article.
The reflected waves generated by a flat reflector were obtained by using the 2-D finite difference method. In order to remove direct wavelets, the wavefield without the reflector was subtracted from the total wavefiled of the reflector model. We then obtain the wavefield containing only reflected waves. Figure 2a shows an example of the shot gather of reflected wavefield. In the case of the RN model, band-limited random noise (5-50 Hz) was added to the synthetic data containing only signal components (reflections) so that the S/ N ratio was 0.3 ( Figure 2b). In Figure 2b, reflected waves can hardly be detectable due to masking effect by random noise. In the case of the RH model, on the other hand, to compare results between random media of different characteristic lengths, wavelengths have to be described with reference to the characteristic lengths of random media. The product of the wavenumber k and the characteristic length a is used as an index for describing effects of random heterogeneity on seismic waves. In the present cases, the ka values at the dominant wavelengths are about 0.33 (a=10 m) and 1.65 (a=50 m), respectively. According to the classification by Wu and Aki (1988), our heterogeneous models are categorized as " Mie scattering regime" where strong scattering may occur and full waveform modeling is required. In order to remove direct wavelets, the total wave field calculated with the model shown in Figures 1b and 1c was subtracted from the wave field in a model with a constant velocity of 3800 m/ s to produce the wave field containing the reflected/ scattered wave field. Figure 3a shows an example of the shot gather from the scattered wave field in the case of a=10 m for source-receiver intervals of 5 m. Similarly, Figure 3b shows an example of the shot gather from the scattered wave field in the case of a=50 m for source-receiver intervals of 5 m. Although we can clearly see the reflection event in each shot gather shown in Figure 3, the shot gathers are full of chaotic diffraction patterns originating from random heterogeneities. Fig. 3. Examples of common-shot gather of a scattered wave field calculated for the model with different spatial sampling intervals and characteristic heterogeneity sizes: spatial sampling intervals of 5 m for the case of for a=10 m (Fig. 1b) and spatial sampling intervals of 5 m for the case of for a=50 m (Fig. 1c)   Fig. 4. Frequency-wavenumber (f-k) plots of the extracted shot gather of a scattered wave field with different spatial sampling intervals and characteristic heterogeneity sizes: spatial sampling intervals of (a) 80 m, (b) 20 m, and (c) 5 m for the case of for a=10 m The frequency-wavenumber (f-k) diagram is helpful for visualizing the sampling of a continuous wave field (Vermeer, 1990). The time window (from 0.65 to 1.05 s) including only scattered wave fields was extracted from each shot gather to calculate an f-k plot. Figures 4a through 4c show f-k plots of the extracted shot gather from the scattered wave field in the case of a=10 m for source-receiver intervals of 80, 20, and 5 m, respectively. Similarly, Figures 5a through 5c show f-k plots of the extracted shot gather from the scattered wave field in the case of a=50 m for source-receiver intervals of 80, 20, and 5 m, respectively. According to the spatial Nyquist sampling criterion defined in equation (1), θ=90). Thus, spatial sampling less than 45 m is sufficient to prevent spatial aliasing of the scattered wave field. In the case of the 80 m spatial sampling interval of Figure 4a and 5a, the sector of strong amplitudes in the f-k plot would be severely truncated, causing wrap-around effects. On the other hand, in the case of zero-offset defined by equation (2), spatial sampling of less than 22.5 m is sufficient to prevent spatial aliasing.

Results
Three types of data processing, conventional CMP stacking, poststack migration and prestack time migration (PSTM), were applied to the two types of model described above.

CMP stacked sections
Conventional CMP stacking was applied to both random noise (RN) model and random heterogeneous (RH) model. The shot gathers were sorted into CMP gathers and corrected for NMO using constant velocity of 3800 m/ s, and finally stacked. Figures 6a through 6c show the CMP stacked sections for the RN model data shown in Figure 2b with different source/ receiver intervals at 80, 20 and 5 m, respectively. We can see that the S/ N ratio becomes larger with denser source/ receiver arrangements. The difference of the S/ N ratio among Figures 6a through 6c becomes larger with increasing the numbers of sources/ receivers. The low quality of both sides of CMP sections is due to the low fold in the CMP gathers at the margins of target area. 10, and 2.5 m, respectively. Although CMP stacking can act as a powerful mechanism for suppressing multiples and for the attenuation of many types of linear event noises such as airwaves and ground roll, we can see no significant differences among Figures 7a through 7c and among Figures 8a through 8c. However, a close examination of these sections reveals that image space aliasing occurs in the case of a CMP interval of more than 22.5 m (Figure 7a and 8a). Note that the effect of image space aliasing in the case of a=10 m is larger than the case of a=50 m. In each section of Figures 7 and 8, we can see a reflector at around 1.1 sec. and many discontinuously subhorizontal and dipping events that partly correlate with velocity heterogeneities of the model. Gibson and Levander (1988) mentioned that the limited bandwidth of the propagating seismic signal and spatial filtering attributable to CMP stacking cause these events, bearing no simple relation to the velocity anomalies of the model. While the reflector can be seen clearly from the chaotic background noise, we can see some arrival time fluctuations and amplitude variations in the observed reflector. These variations are attributed to the scattering effect of the heterogeneous media whose scale is smaller than the wavelength. In Figures 7 and 8, we can see no significant arrival time fluctuations but some amplitude variations in the observed reflector. These amplitude variations are attributed to the scattering attenuation (sometimes called apparent attenuation) in the heterogeneous media. When the heterogeneous scale is small, the amplitude is affected by the heterogeneity but the travel time is not strongly affected by the heterogeneity. In this situation, the assumptions of CMP stacking and simple hyperbolic reflection pattern can be fulfilled.  (Stolt, 1978) for a RN model with different source/ receiver intervals at 80, 20, and 5 m, respectively. The  Figure 9 are 40, 10, and 2.5 m, respectively. Although the resulting migrated sections suffer from the inadequate cancellation of migration smiles, we can see that the S/ N ratio becomes larger with denser source/ receiver arrangements. Figures 10a through 10c show poststack migrated sections using f-k migration (Stolt, 1978) with a random heterogeneous model for the case of a characteristic heterogeneity size of 10 m (a=10 m) with different source/ receiver intervals at 80, 20, and 5 m, respectively. Similarly, Figure 11 is the same as Figure 10 except for characteristic heterogeneity sizes of 50 m (a=50 m). We can see that numerous small segments are still detectable even after the poststack migration and that the results of poststack migration for the different heterogeneous models differ with different source/ receiver intervals. Although we can see  (Figure 10a and 11a). Note that the effect of image space aliasing in the case of a=10 m is larger than the case of a=50 m. In general, migration can improve lateral resolution by correcting the lateral mispositioning of dipping reflectors or collapsing diffraction patterns caused by a point scatterer. However, the application of poststack migration here does not improve seismic images in heterogeneous media. It is thought that the reason is that multiple-scattering effects in small-scale heterogeneities do not satisfy the assumption of migration theory based on single scattering. Although migration techniques assume that the seismic data to be migrated consists only of primary reflections and diffractions, these wave fields are attenuated and distorted by heterogeneities and multiple scattered wave fields are generated, producing apparent discontinuities in reflectors or diffractors.

Prestack time migrated sections
In this paper, we obtained PSTM sections using a diffraction stacking method proposed by Matsushima et al. (2003). Figures 12a through 12c show PSTM sections for a RN model with different source/ receiver intervals at 80, 20, and 5 m, respectively. We can see that the S/ N ratio becomes larger with denser source/ receiver arrangements. Figures 13a through 13c show the PSTM sections using a diffraction stacking method (Matsushima et al., 2003) for a random heterogeneous model with a characteristic heterogeneity size of 10 m (a=10 m), as shown in Figure 1b with different source/ receiver intervals at 80, 20, and 5 m, respectively. Similarly, Figure 14 is the same as Figure 13 except for characteristic heterogeneity sizes of 50 m (a=50 m). Each PSTM section is full of migration smiles, producing the appearance that the section is heavily over-migrated, thus reducing the quality of the image. A possible explanation of this phenomenon is that the wave field is distorted by heterogeneities, which in turn produce apparent discontinuities in reflectors or diffractors. These discontinuities do not have associated diffraction hyperbolae, so that the migration, instead of collapsing the absent hyperbolae, propagates the noise represented by the discontinuity along wavefronts. As a result, the seismic section is full of migration smiles that are heavily over-migrated. Warner (1987) pointed out that deep continental data are often best migrated at velocities that are up to 50 % less than appropriate interval velocities from crustal refraction experiments or directly from stacking velocities. His explanation for this behavior is that near surface features distort and attenuate the seismic wave field and produce apparent discontinuities in deep reflections. During the process of migration, reflections are invented in order to cancel out the missing diffractions thereby producing a smiley section that appears over-migrated. Although PSTM is expected to provide more realistic images compared to conventional poststack migration (Gibson and Levander, 1988), we can see no significant differences among Figures 13a  through 13c, and also among Figures 14a through 14c. Similar to the case of poststack migration, the reason is thought to be that multiple-scattering effects in small-scale heterogeneities do not satisfy the assumption of migration theory based on single scattering.  Figures 15a thorough 15c show the center trace of the corresponding section in the case of the RN model with three different spatial sampling intervals. We can see that there is little difference of the S/ N ratio among the data processing variants when the spatial sampling Fig. 15. Comparison of the center trace of the corresponding section in the case of the RN model with three different data processing and three different spatial sampling intervals Fig. 16. Comparison of the center trace of the corresponding section in the case of the RH model (a=10 m) with three different data processing and three different spatial sampling intervals interval is 80 m (i.e., the number of sources/ receivers is small). However, the difference of the S/ N ratio becomes larger with shortening the spatial sampling interval (i.e., increasing numbers of sources/ receivers), and the PSTM does a much better job of imaging the reflector. Huygens' principle explains this mechanism as follows. A reflector is presumed to consist of Huygens' secondary sources, in which case imaging a reflector is considered to be equivalent to imaging each point scatterer separately and summing the imaged point scatterers at the end (Matsushima et al., 1998). A point scatterer can be delineated more appropriately by PSTM than by CMP stacking or poststack migration. In this case, an adequate zero-offset section cannot be obtained by CMP stacking without dip moveout (DMO) corrections.  Figures 17a thorough 17c show the center trace of the corresponding section with three different spatial sampling intervals in the case of the RH model (a=10) and RH model (a=50), respectively. We can see that there is little difference of the S/ N ratio between different spatial sampling intervals except the shallow part of each section (less than 0.2 sec.) in each data processing. However, the difference of the S/ N ratio among the data processing variants is obvious, that is, the PSTM does a much better job of imaging the reflector in the randomly heterogeneous media. The reason can be explained by the Huygens' principle as described above.

Discussion
It is important to discriminate between two different types of noise: a random noise in time series and a noise-like wave field produced from random heterogeneity. One may regard the scattered waves generated from heterogeneous media as a random noise in field seismic data. Some authors (e.g., Matsushima et al., 2003) have added random noise to their synthetic data for simulating field seismic data. However, the noise is a consequence of the wave phenomena in heterogeneous media, and is not same as the noise that randomly appears in the time-series (Levander and Gibson, 1991). Scales and Snieder (1998) concluded that the noise in a seismic wave is not merely a time-series which is independent from the original seismic wave but a signal-induced wave mostly consisting of scattered waves. This is important for seismic data processing but not well recognized in the field of seismic explorations. To generate the signal induced noise, the noise should be calculated from the interaction between the small-scale random heterogeneity and the original seismic wave. However, we should also note that the small-scale random heterogeneities are not known and should be estimated by other methods like numerical experiments.
We demonstrate that one can obtain better final section in terms of its S/ N ratio as the intervals of spatial sampling becomes shorter (with increasing the numbers of sources/ receivers) for the case of random noises model where the added random noise is a completely independent time-series against seismic traces. Thus, this type of random noise cancels each other by applying CMP stacking, poststack migration, and PSTM. On the other hand, scattering waves generated from random media is now recognized as a mutually dependent noise among the seismic traces, which indicates the interaction between the short-wavelength heterogeneity and the source and reflected wavelet. Although these scattered waves appear as random noises, they are thought to be an accumulation of many scattered waves which themselves partially coherent. Thus, this type of scattering noise should be categorized into coherent noise if we classify noise types. In general, coherent noise can not be reduced after processing the data, merely by increasing the source strength or shortening the sampling interval.
It is widely believed that highly dense spatial sampling increases the quality of final seismic sections. There are two aspects to the improvement of the quality. One is that a shorter spatial sampling interval can reduce the migration noise caused by spatial aliasing. The other is that the increase in the number of sources/ receivers raises the effect of signal enhancement to increase the S/ N (signal to noise) ratio.
In random heterogeneous media, three types of data processing, conventional CMP stacking, poststack migration, and PSTM, were applied and compared to examine different responses to different sampling intervals. Each data process without data space aliasing achieves very similar final sections for different sampling intervals. Safar (1985) studied the effects of spatial sampling on the lateral resolution of a surface seismic reflection survey when carrying out scatterer point imaging by applying migration, and found almost no effect of spatial sampling on lateral resolution. Safar (1985) also demonstrated the generation of migration noise caused by a large sampling interval. Migration noise is a consequence of spatial aliasing that is related to frequency, velocity, and dip of a seismic event. A shorter sampling interval cannot improve spatial resolution very much, even if there is no noise. The same conclusion was obtained by Vermeer (1999). The results we have obtained correlate well with those of these previous studies.
Our numerical experiments indicate that the highly dense spatial sampling does not improve resolution of the section except the shallow part of the section when the subsurface structure contains random heterogeneity, even if the interval of spatial sampling becomes shorter than the Nyquist sampling interval. However, we found the existence of a significant difference among the data processing variants. We demonstrate that the prestack migration method has the advantage of imaging reflectors with higher S/ N ratios than typically obtained with the conventional CMP stacking method with/ without the poststack migration. We explained the possible mechanism by the Huygens' principle. A point scatterer can be delineated more appropriately by PSTM than by CMP stacking or poststack migration.
In our numerical experiments for RH models, two different heterogeneity sizes (a=10, 50 m) with three different spatial sampling (5, 20, 80 m) were applied. Our numerical experiments show that the effect of image space aliasing depends on the relationship between the heterogeneity size and the spatial sampling interval. Frequency components of scattering waves generated from random media depend on the heterogeneity size. When spatial sampling is too coarse, steeper-dip events are relatively aliased. To avoid spatial aliasing in heterogeneous media, it is important to know how dense the source/ receiver arrangements should be in data acquisition. Narrower interval in spatial sampling can provide a clearer image of heterogeneous media. Qualitatively, spatial sampling should be smaller than the size of heterogeneities. Further consideration on quantifying the relationship between spatial sampling and the size of heterogeneities is needed. We also note that the small-scale random heterogeneities are not known and cannot be effectively estimated prior to data acquisition.

Conclusions
We have shown from the viewpoint of spatial sampling how the two different types noise, a random noise in time series and a noise-like wavefield produced from random isotropic heterogeneity, influence the final section. We use a 2-D finite difference method for numerically modeling acoustic wave propagation. In the presence of the time-series random noise, a final section can be obtained with a higher S/ N ratio with shortening the interval of spatial sampling, that is, the increasing the numbers of sources/ receivers improve the reflection image. On the other hand, in the case of random heterogeneous model, a final section is influenced by the interval of spatial sampling in different way as that of timeseries random noise. Highly dense spatial sampling does not seem to improve the final quality of a section regardless of the relationship between the spatial sampling interval and the characteristic size of heterogeneities, even when the interval of spatial sampling is smaller than the Nyquist interval. We have pointed out the importance of discrimination between two different types of noise: a random noise in time series and a noise-like wave field produced from random heterogeneity. We have also demonstrated that the prestack migration method has the advantage of imaging reflectors with higher S/ N ratios than typically obtained with the conventional CMP stacking method with/ without the poststack migration in both RN and RH model, which can be explained by the Huygens' principle.

Acknowledgments
The author greatly acknowledges the thorough reviews and constructive comments of an anonymous reviewer, which helped increase the quality of the manuscript. This study was