Acquisition geometry analysis with point-spread functions

Seismic data are traditionally acquired based on spatial sampling requirements, noise properties and budgetary constraints. However, designing a survey without taking into account the complexity of the subsurface may result in an image without the expected quality. Also, the subsequent preprocessing and processing steps may exploit or misuse the acquired data. The design should therefore incorporate the complexity of the subsurface and the (pre)processing steps that will be followed. We propose an analysis method that evaluates if the proposed combination of survey design, preprocessing and processing for a specific subsurface model fulfils a pre-defined quality criterion. With our method, we estimate a set of point-spread functions that correspond to the chosen combination,and we analyse their resolution and illumination-detection properties in the spatial and wavenumber domains, respectively. The estimated point-spread functions include the scattering and propagation effects generated by the sub-surface, including internal multiples. We show that in some cases, the use of internal multiples in imaging can improve amplitude and resolution compared with the use of primaries only. The proposed analysis method is also used to evaluate the effect of blending noise when blended acquisition is carried out.

Seismic surveys are traditionally designed based on properties such as fold and azimuth distribution (Vermeer, 2012).However, the illumination and detection properties of the survey design could be severely influenced by a complex subsurface.Today, in many parts of the world, subsurface models are available from legacy surveys.These can be of great value for designing future surveys as the acquisition parameters can be tuned to the particular subsurface structure to provide optimum imaging results.
After seismic data have been acquired, subsequent preprocessing, processing (imaging and velocity estimation) and interpretation (layer property estimation) are carried out.Traditionally, these steps are carried out sequentially and are considered as more or less independent.However, more and more it is recognized that they are strongly interrelated.For instance, the performance of certain data (pre)processing algorithms may strongly depend on the availability of certain wellsampled gathers, i.e. on the data acquisition, and the interpretation could require a certain minimum quality of the seismic image.For that reason, it does not make sense to design an acquisition geometry without taking into account what data (pre)processing steps will follow and what image quality is required.Although each step of the chain will influence the result in a different way, the final product will be affected by the complete sequence.
In this research, therefore, the whole seismic imaging method is considered, i.e. a 'holistic approach' is followed.The focus, however, is on data acquisition.It means that, in the ideal situation, the acquisition geometry is designed in such a way that -after data processing in a particular pre-defined way -the results have a certain pre-defined image quality.Here, the best design is the cheapest that leads to the goal.As this method is carried out in the pre-acquisition phase, it is fully carried out on synthetic data.
For benchmarking an imaging system, the imaging response of a collection of point scatterers is analysed.Ideally, for a point scatterer in the subsurface, the result of the imaging process should be a single, well-focused point.However, in practice, a coarse, irregular spatial sampling, a limited acquisition aperture, a complex subsurface and the limitations of the seismic (pre)processing cause the scattered energy to be imaged in non-perfect way: instead of a point, the resulting image has a distorted shape which resembles a band-limited wavelet.The result is called a point-spread function (PSF), and it is defined as the impulse response of the complete system of seismic acquisition, preprocessing and imaging at the location of the point scatterer.The collection of point scat-ters can be used to determine this impulse response at many locations.
The PSF, which is also called resolution function, has been studied extensively for prestack migration quality analysis (Berkhout, 1984;Beylkin, 1985).Several factors that influence the PSF are the acquisition geometry, frequency content, background velocity model and the migration algorithm.Vermeer (1999) studies the effect of frequency content, acquisition aperture and geometry on the resolution for a single scatterer in a homogeneous medium.Beylkin (1985) describes migration as a mapping from the acquisition coordinates space to the spatial frequency domain.The range of angles present will determine the spatial resolution.This domain can be obtained numerically via ray tracing for prestack migration.Ray tracing methods provide a computationally efficient way to model the illumination and detection properties of a survey in a complex subsurface (Lecomte, 2008).However, these properties can be either exploited or undermined by the chosen imaging algorithm.
In current imaging, the multiples (surface-related multiples as well as internal multiples) are not considered as noise, but instead they are considered as valuable signal.The handling of multiple reflections by the imaging algorithm may contribute to the imaging results, potentially relaxing the spatial sampling requirements.In Marchenko imaging for instance (Wapenaar et al., 2014), the estimated Green's functions include the internal multiple reflections that contribute to reflectivity estimation.Free-surface-related multiple reflections can also be taken into account (Singh et al., 2015(Singh et al., , 2017)).Least-squares and inversion approaches offer the same possibility (Brown and Guitton, 2005;Malcolm et al., 2009;Zhang and Schuster, 2014;Tu and Herrmann, 2015).In particular, the multiples may contribute if imaging by primaries is impossible, e.g. in the cases that primaries do not reach the particular area to be imaged with sufficient signal strength but multiples do.Consequently, the image quality depends on the subsurface reflectivity: only if sufficiently strong reflectors are present, the (internal) multiples are likely to contribute to the imaging; if not, their contribution can be ignored.In the first case, the use of an imaging algorithm such as full-wavefield migration (FWM) (Davydenko, 2016) that uses all the multiples can compensate for illumination in areas where the illumination by primaries only is not sufficient.This approach may reduce the sampling requirements in the acquisition surface, thereby reducing costs.In a similar manner, if a data set known for having strong multiples is processed with a primary-only migration, the signal corresponding to these multiples will not be processed adequately, resulting in undesired artefacts.In such a case, more data may be required in order to obtain an image with a sufficient signal-to-noise ratio.Therefore, when designing an acquisition geometry, the processing methodology should be taken into account for an optimum result.
The so-called focal beam analysis (Berkhout et al., 2001) provides a method to analyse the acquisition geometry of a survey based on the illumination and detection properties of the sources and receivers, respectively.The complexity of the subsurface can be taken into account (Veldhuizen et al., 2008) as well as the illumination by all multiples (Kumar, 2015).However, the process is only efficient for analysing a couple of target points in the subsurface.Repeating this process for all subsurface locations makes it too expensive.
Blended or simultaneous source acquisition (Beasley et al., 1998;Berkhout et al., 2008;Bouska, 2010) allows to use seismic sources whose responses overlap in time and space as well as in temporal and spatial frequency bands.Compared with conventional acquisition, a higher efficiency can be achieved.With algorithms such as FWM the blended data can be directly processed, i.e. without deblending, reducing the total number of shots and therefore increasing the computational efficiency.For other, more conventional algorithms, it is desired to deblend the records before further processing.In this case, the design of the blended acquisition has to be such that it allows for an optimal reconstruction of the unblended wavefields.As different deblending methods may pose different requirements on the spatial sampling of the input data and the allowed maximum blending factor, it is necessary to take these into account when designing the acquisition geometry.
In order to measure the performance of the seismic experiment, it is necessary to define an analysis method.The result of such an analysis would be a quantitative measure that allows one to judge whether or not the pre-defined quality criteria are met.In this paper, we focus on this analysis method, using the following approach.With the available velocity model we compute a set of PSFs by first modelling a seismic data set for the original reflectivity model and then modelling a second seismic data set for an almost identical model: the only difference is that a grid of point scatterers is added to the original reflectivity.Both data sets are preprocessed and imaged in the same way.The difference between the two final images is the sought set of PSFs.By computing the PSFs in this way, the effect of acquisition, preprocessing and imaging are all taken into account, and specially, if present, the contribution of multiples to the imaging.The PSFs are quantitatively analysed to assess the quality of the seismic experiment.

F R A M E WO R K O F S U RV E Y A NA LY S I S
We describe 3D seismic data using the matrix notation introduced by Berkhout (1982).For each frequency component, a data set can be formulated as where P is the seismic signal recorded at depth level z d , generated by the sources at level z s .Each element contains the amplitude and phase information of one-shot record, recorded by one receiver, for the frequency component under consideration.S and D are the source and receiver matrix, respectively.Together they describe the survey geometry (number of sources and receivers and their locations), as well as the source and receiver properties (directivity, spectral properties etc.).Matrix X represents the transfer operator of the medium and contains propagation and reflection effects.It can be considered to be the ideal-sampled seismic data.
The modelling is carried out with full-wavefield modelling (FWMod) (Berkhout, 2014).In FWMod, the modelling is recursive in depth.First, the downgoing wavefield is computed and then the upgoing wavefield.These two steps are called one round trip.The implementation requires multiple round trips, each round trip adding one additional order of multiple reflections: source wavefield.Matrices U − and U + are the upgoing and downgoing full-wavefield propagators, respectively.They are computed as follows: Here W + (z k+1 , z k ) is the downward propagation operator from depth level z k to depth level z k+1 ; W − (z k−1 , z k ) is the upward propagation operator from depth level z k to depth level z k−1 ; T + (z k ) is the transmission operator of the downgoing wavefield crossing depth level z k from above; T − (z k ) is the transmission operator of the upgoing wavefield crossing depth level z k from below.
Therefore, operator U + (z n , z m ) includes all propagation and transmission effects in the downgoing direction from depth level z m to depth level z n , and similarly for operator U − (z n , z m ).The use of this type of modelling allows to specify reflectivity and transmissivity independently from propagation velocity.The outcome of the modelling process after N round trips (corresponding to the order of multiples being modelled) is the seismic signal p − N at the recording surface z d .The modelling is repeated for multiple sources at the surface z s , resulting in the data matrix P(z d ; z s ) in equation (1).Note that the number of round trips N is also indicative for the cost ratio between using FWMod and using primaries-only modelling.
To get a blended data set P (z d ; z s ), we multiply the unblended data with blending operator (z s ) (Berkhout et al., 2008): (4) Matrix (z s ) contains the blending code, i.e. time delays or phase shifts for simultaneous sources.Each column of (z s ) represents a blending experiment and each row the specific source to be blended.To simplify the notation, we assume that the sources and the receivers are located at the same acquisition surface, so that z d and z s can be omitted.

Numerical modelling
In a numerical experiment, the calculation of a point-spread function (PSF) involves the modelling of seismic data with a chosen acquisition geometry and subsurface model.As mentioned in the introduction, to obtain the PSFs corresponding to a complex subsurface, two data sets are modelled.The first, P 1 , corresponds to the original reflectivity model, and the second, P 2 , corresponds to the second reflectivity model which is the sum of the original reflectivity model and an additional grid of point scatters.The original velocity model is used for modelling both data sets.Although it is not physically possible to decouple velocity from reflectivity, FWMod offers this option and we use it as a means to simulate the presence of the point scatterers while taking into account the complexity of the subsurface.The obtained data sets are Here P is the data related to the original subsurface, and P + P is the data related to the subsurface plus point scatterers.These scatterers have a low reflectivity.This is to make sure that the multiples in P 2 that contribute to the imaging can be considered to be present in P only, and not in P. In the case of blended acquisition, these data sets become P 1 and P 2 .

Imaging
We choose full-wavefield migration (FWM) to illustrate the possible advantages that the use of the full wavefield in imaging could bring to survey design.FWM is an inversion-based algorithm that aims to estimate the reflectivity for a given data set.We assume the velocity model is known.However, it is possible to combine reflectivity estimation with velocity estimation in the so-called joint migration-inversion (JMI) (Staal, 2015).Here we focus on the imaging only, i.e. on the reflectivity estimation.In our examples, we simplify the reflectivity to be angle-independent (i.e. the reflection operators become reflection coefficients and, consequently, R ∩ and R ∪ become diagonal matrices).In FWM, the reflectivity model is iteratively updated by computing the data residual, i.e. the difference between the observed data and the data forward-modelled through FWMod.This data residual is mapped to the model space via adjoint modelling.After proper scaling, the update is applied to the current reflectivity model.A complete mathematical description is offered by Davydenko (2016).
Two images are obtained after processing P 1 and P 2 in exactly the same way.These are I 1 and I 2 , respectively: Here, I is the reflectivity image of the subsurface, i.e. the collection of estimated Rvalues for all depth levels obtained after preprocessing and processing of data set P 1 .The term I + δI on the right-hand side of equation ( 8) is the reflectivity image of the original subsurface model, that differs from I by the term δI due to the possible influence of the scattering P from equation ( 6) on the imaging.The term I is the reflectivity image of the point scatterers, generated by P. As we chose the reflectivity of the point scatterers to be small compared with the reflectors in the original model, we neglect the possible multiples generated by the scatterers.Therefore, for FWM, the image contribution δI can be neglected as well.By subtracting the two reflectivity images I 2 and I 1 , the reflectivity image of the point scatterers, I, is obtained.The latter is the collection of PSFs mentioned earlier.The theory presented in this and the previous section is valid for 2D and 3D seismic experiments.
However, for convenience, we present 2D examples in the remainder of this paper.To illustrate the concept of PSFs as discussed in this section, we selected a part of the Marmousi model (Fig. 1) to model seismic data using FWMod.One data set is modelled using the original reflectivity model (Fig. 1b).We then added a collection of 81 point scatterers to the same model, distributed regularly in the subsurface (Fig. 1c), and modelled a second seismic data set.The original velocity model (Fig. 1a) was used in both data sets.The spatial sampling of the sources and receivers was ideal in the two data sets, and both sets were (pre)processed in exactly the same way.Finally, the two obtained reflectivity images were subtracted, resulting in the set of PSFs shown in Figure 2. It is these PSFs that can be further analysed.This analysis is the topic of the next section.

QUA N T I TAT I V E A NA LY S I S Resolution and illumination detection
To quantify the quality of the final image, the individual PSFs are analysed.Each of them can be considered separately in order to obtain the resolution in the image domain and involved angle distribution in the wavenumber domain at the location of each scatterer.Figure 3 shows the results for the ideal (band-limited) image of a point scatterer in the spatial and wavenumber domains.The PSF in the spatial domain shows a well-focused event whose sharpness is limited by the bandwidth of the seismic data.The wavenumber spectrum on the right shows a full range of illumination and detection angles; see the yellow area.Again, its limited radial extent is the result of the band limitation.Figure 4 shows a set of one hundred PSFs.In this example, the acquisition geometry consists of a uniform spatial sampling of sources and receivers and the model velocity is constant.One of the PSFs has been selected for a zoomed view in Figure 4(b) and its corresponding spectrum is shown in Figure 4(c).Ideally, the point scatterer should be present at its true position only, being (x, z) = (1813, 2488) m.However, the energy is spread in an area tens of metres around that position.In addition, the PSF is no longer circularly symmetric; it has instead a shape whose cross-sections resemble band-limited wavelets.For a Ricker wavelet as source pulse, a Ricker wavelet shape is obtained in the vertical direction.A Gaussian spatial wavelet is obtained in the horizontal direction (von Seggern, 1991).These cross-sections are not strictly vertical or horizontal but have an orientation related to the illumination and detection angles at the position of the point scatterer.Each PSF may have a different orientation which we call the inclination angle.This inclination is more evident for the PSFs closer to the edges, as the acquisition aperture seen by the corresponding point scatterers is relatively small compared with the aperture seen by point scatterers in the central section.The wavenumber spectrum shows a bowtie shape.Compared with Figure 3, the higher angles are missing.Also, it can be seen that the very low frequencies are no longer present.The angles and frequencies present in the wavenumber spectrum correspond to the combined effect of angles of illumination by the sources and angles of detection by the receivers and their combined spectral properties.
To compute a quantitative measure of resolution and to identify the inclination angle of the PSF, we analyse the two cross-sections that pass across the point of maximum amplitude of the PSF.As mentioned, two axes can be distinguished in a PSF.First, a lower resolution axis, along which the PSF has a Gaussian shape.Second, usually perpendicular to the first, a higher-resolution axis.For the latter, we compute its envelope by taking the absolute value of its Hilbert transform.Then, we measure the distance between the points with half of the peak amplitude.The low-resolution axis is defined as the cross-section with the longest distance between these two points.The high-resolution axis has the shortest distance between those two points.These distances are our measures of resolution (see Fig. 5).
The resolution measurements can be computed for all PSFs in the subsurface.Figure 6 shows the measure of resolution for the set of PSFs in Figure 2. The high-and lowresolution measures are displayed in Figure 6(a) and (b), respectively.The arrows indicate the orientation of the respective axes.Their length corresponds to the respective measures of resolution, i.e. the longer the arrow, the lower the resolution.In addition, this resolution has been computed for the whole subsurface via interpolation between PSFs and colour coded in the background.The measurements of high resolution (Fig. 6a) show a tendency of decreasing with depth.The main reason is the increase of seismic velocities with depth.The yellow part at the lower-right section coincides with the higher velocities observed in Figure 1(a).The measurements of low-resolution (Fig. 6b) show a similar trend.As the acquisition aperture seen by the point scatterers decreases with depth, the range of illumination and detection angles decreases and therefore so does the resolution.

Effect of imaging
In traditional prestack depth, migration multiples are considered as noise.Therefore, they need to be eliminated before imaging.In full-wavefield imaging methods, multiples are not eliminated but used to estimate reflectivity.As pointed out in the introduction, multiple reflections may reach areas in the subsurface that primaries do not.In particular, internal multiples, generated deeper in the subsurface than surface-related multiples, may provide additional illumination from below.In this section, we use the proposed analysis method to evaluate the effect of using internal multiples for imaging.For synthetic models and different acquisition scenarios, we compute sets of point-spread functions (PSFs) using full-wavefield migration (FWM) as well as primaries-only migration.For the latter, we use a least-squares migration with a correlation-type imaging condition, which only considers the primary wavefield in its modelling engine (FWMod with N = 1).Therefore, multiple reflections do not contribute to the reflectivity estimation in that case.Subsequently, we measure resolution, analyse illumination and detection and evaluate the differences.First, we consider a set-up with two synthetic models of a salt dome embedded in a layered medium (Fig. 7).The velocity of the deepest layer, i.e. below z = 1000 m, is 2000 m/s for model 1 (Fig. 7a) and 3300 m/s for model 2 (Fig. 7b).The density is kept constant.The reflectivity of the interface at z = 1000 m is 0.05 and 0.34 for models 1 and 2, respectively.A different reflectivity at this interface means that the amplitude of the reflected upgoing wavefield, and hence the strength of the (internal) multiples, is different.The amplitude of these internal multiples is proportional to the reflectivity of that interface, i.e. weak internal multiples in the first case and strong internal multiples in the second.These internal multiples may be used for imaging depending on the chosen imaging method.It is of interest to have a good image of the area beneath the flank of the salt, as we assume that a target of economic importance is located there.
We consider an acquisition scenario in which the receivers are placed at the complete acquisition surface, while the sources are located at the right side only, i.e. from x = 1400 m to x = 2500 m.Seismic data are modelled with FWMod and migrated using primaries-only migration (Fig. 8). Figure 8 Figure 8 shows that with primaries-only migration, the reflector at z = 1000 m in models 1 and 2 is clearly imaged.It is much stronger in Figure 8(b) than in Figure 8(a), as expected.As the sources are only located at the right side, given the geometry of the model, there will not be primary reflections generated by the flank of the salt that are recorded by the receivers.Therefore, it is not possible to obtain an image of this section using primaries-only migration.None of the images contains the flank of the salt dome.Additionally, strong artefacts appear inside the salt dome.The multiple reflections between the top of the salt dome and the water bottom cannot be explained by this type of imaging.Hence, they appear as spurious events beneath the top of the salt.
Figure 9 shows the imaging results for both models using FWM.Even for the case of weak internal multiples, in Figure 9(a) a part of the flank of the salt is imaged.The artefacts in Figure 8(a, b) are no longer present.The flank is better imaged when stronger internal multiples are present (Fig. 9b).It is important to note that the synthetic data generated for these examples are noise-free.In the presence of noise, it is possible that the weak internal multiples are not strong enough to obtain a signal-to-noise ratio (SNR) suitable for imaging the target zone.Therefore, stronger internal multiples may be needed, as those used for the imaging in Figure 9(b).
To analyse the difference in resolution and illumination, we calculate a set of PSFs for model 2. Figure 10 shows a comparison between the PSFs obtained with primaries-only migration and FWM.The PSFs in the right-hand side for both models are almost identical.The point scatterers in this part of the model are well illuminated by the primary wavefield.Therefore, primaries-only migration suffices to obtain a good image.This is not the case for the point scatterers beneath the flank of the salt dome.There, the PSFs from the primaries-only migration (Fig. 10a) have a lower resolution and amplitude than the PSFs from the FWM (Fig. 10b).This is because the internal multiples generated between the salt dome and the interface at z = 1000 m are used for imaging.
Figure 11 shows a comparison between the PSFs centred at (x, z) = (1000, 640) m (lower red boxes in Fig. 10a and b). Figure 11(b) shows the PSF obtained through FWM.It has a higher amplitude than the one obtained through primariesonly migration (Fig. 11a).The energy from the internal multiples is now used for imaging; therefore, it contributes to the reflectivity estimation.The red lines in both figures correspond to the low-resolution axes.In order to compare their resolution, they are normalized and plotted in Figure 11(c).The lowresolution axis corresponding to FWM has a better resolution than the one corresponding to primaries-only migration.The wavenumber spectrum of the FWM PSF (Fig. 11e) shows a wider range of illumination and detection angles than the spectrum of the primaries-only PSF (Fig. 11d).The improve-ment in resolution and amplitude is more evident in the region above the salt dome.Figure 12 shows the illumination and detection angles for the PSF at (x, z) = (1000, 240) m (upper red boxes in Fig. 10a and b).Here there are strong internal multiples reflected between the top of the salt and the water bottom.As there are no sources at the acquisition surface directly above this region, the additional illumination by the internal multiples helps to improve the image of the zone considerably.
To illustrate the procedure for a more complex subsurface model and a more challenging acquisition scenario, we calculate a set of PSFs for the Marmousi model (Fig. 1).The acquisition geometry consists of receivers at the complete surface and a gap in the source distribution from x = 1000 m to x = 3000 m, assuming there is an obstacle in this region.As in the previous examples, our seismic data were modelled without surface-related multiples and we focus on the effects generated by the internal multiples.Figure 13 shows a comparison between the PSFs obtained with primaries-only migration and the ones obtained with FWM.It is visible that for both cases the acquisition gap has a considerable effect.When comparing with the PSFs in Figure 4, which correspond to a full source sampling, artefacts appear in both images in the region beneath the source gap, at the centre of the model.The primaries-only migration result is more affected as it relies on the primary reflections generated by the wavefield of the adjacent sources.Therefore, the shallower PSFs have a decreased resolution.The difference plot shows that the PSFs in this part have more energy in the FWM case.This difference can be attributed to the internal multiples used for the imaging, as the acquisition geometry is kept constant.
The additional illumination provided by the internal multiples could help to reduce the spatial sampling requirements of the survey.To illustrate this we compare the PSFs obtained from two acquisition geometries, see Figure 14.The source interval is 125 m for the first geometry and 500 m for the second.The spatial sampling of the receivers is 12.5 m for both geometries.Figure 14(a) shows the set of PSFs for the first geometry and Figure 14(b) for the second.Figure 14(c) shows their difference.The differences are barely noticeable.This means that even with the sparse source distribution, a good result was obtained.Note that in general the 'interpolating effect' of the use of multiples is expected to be largest when surface multiples are used.However, we did not use these as our focus is on the contribution of internal multiples.
The examples in Figure 14 allowed to analyse the effect of different acquisition geometries while the processing algorithm remained the same.The analysis method shows the joint effect of illumination by the sources and detection by the receivers for all locations in the subsurface.If the separate effects of the sources and the receivers are needed, for instance, in survey design, our method could be complemented with the focal beam analysis method for a couple of target locations.

Effect of preprocessing
The preprocessing used to deblend seismic records can leave undesired (blending) noise.In this section, we examine the blending noise from the PSFs. Figure 15 shows a comparison of the PSFs obtained for the model in Figure 1 (upper-left section), for the cases of deblended data (Fig. 15a), pseudodeblended data (Fig. 15b) and their differences with respect to the unblended case (Fig. 15c and d, respectively).The blended data sets correspond to the sources located in the range from x = 0 m to x = 2000 m being blended with those in the range from x = 2000 m to x = 4000 m.Each blended experiment consists of two sources separated by 2000 m and with random time delays ranging from 0 ms to 200 ms between the shots.The deblending algorithm uses an iterative process where coherency filters and thresholding are applied in the commonreceiver domain to remove the blending noise (Mahdad et al., 2011).Figure 15(c) shows the effect of the blending noise once observed in the image domain when using FWM.The measured signal-to-noise ratio (SNR) was 7.52 dB.When using pseudo-deblending only (Fig. 15b), the blending noise in the image domain is considerably larger (Fig. 15d), leading to an SNR of −11.33 dB.
This result shows the combined effect of the acquisition geometry, the chosen preprocessing algorithm and the chosen imaging algorithm.In this case, carpet shooting and detection were used, so there is no further improvement possible for the acquisition geometry.It means that if the blending noise in the final image is considered to be too large, a better-quality deblending algorithm has to be used.Or, alternatively, the blending code (blending factor and time delays) should be chosen such that better deblending results are obtained.Here, we assume that the choice of imaging algorithm would be less relevant and that a lower-quality deblending would always result in a lower-quality final image.A similar type of analysis could be carried out to assess the performance of denoising algorithms in the preprocessing stage.For instance, two different denoising algorithms could be applied to the same data set-with-noise to get two denoised data sets.Subsequently, the corresponding sets of PSFs could be computed, and the corresponding SNR ratios estimated with respect to the PFSs of the noise-free case.This would provide a quantitative criterion for comparing the performance of the different denoising algorithms after the complete chain of acquisition, preprocessing and imaging.

C O N C L U D I N G R E M A R K S
Our proposed analysis method makes use of a priori knowledge of the subsurface to evaluate how a target point in the subsurface is illuminated and how its seismic response is detected.It takes into account the complete chain of seismic acquisition, preprocessing and imaging.The result is a set of point-spread functions (PSFs) that are analysed to determine the local resolution and amplitude fidelity in the final image.In our study, we chose to study the effect of the internal multiples.To this end, the subsurface model should contain the major reflectors which control the amplitude of the internal multiples.If strong internal multiples are present, they may help to illuminate exploration targets in areas where illumination by primaries is insufficient.Subsequently, the signal-to-noise ratio can be improved.It was shown that the use of internal multiples can provide additional imaging energy and improve resolution.To compute PSFs for the case of an imaging algorithm that makes use of multiples, the imaging process has to be carried out twice, with and without point scatterers being present.Their difference section contains the PSFs.In this way, illumination by multiples is included while still ending up with PSFs only.
The design of an acquisition geometry is always a compromise, being limited by the acquisition budget and physical constraints.Therefore, it must be evaluated beforehand to check whether or not the proposed seismic experimental setup of acquisition and (pre)processing meets the pre-defined quality criteria.For example, to achieve the objective, a fullwavefield imaging method may be necessary.On the other hand, if such method is not available, the acquisition geometry has to be modified, probably by adding more sources and/or receivers in order to achieve the required target illumination.The use of blended acquisition could also be evaluated at this stage.

Figure 1
Figure 1 Section of the Marmousi model.(a) Velocity model.(b) Reflectivity model.(c) Original reflectivity plus grid of point scatterers.The latter were amplified 10 times for display.

Figure 2
Figure 2 Set of PSFs corresponding to the model in Figure 1.The acquisition geometry consists of a uniform spatial sampling of sources and receivers.Note the PSFs are different, depending on their location in the model.

Figure 3
Figure 3 (a) Image of a well-illuminated point scatterer.The wavenumber domain transformation (b) contains a complete range of illuminationdetection angles.The resolution is only limited by the seismic bandwidth.

Figure 4
Figure 4 Set of PSFs (a), zoomed image of the PSF in the red box (b), and its wavenumber spectrum (c).The PSF in (b) is no longer symmetrical as in Figure 3(a).The wavenumber spectrum in (c) shows the illumination and detection deficiencies.

Figure 5
Figure 5 (a) Example of a PSF with its low (b) and high (c) resolution axes across the point of maximum amplitude.The amplitudes in (b) and (c) are normalized and treated separately.

Figure 6
Figure 6 (a) The measurements of the high-resolution indicate a decreasing trend with depth.Resolution is specially low in the deeper part of the model.(b) The low-resolution measurement indicates a similar trend of decreasing resolution with depth and structurally complex zones.

©Figure 7
Figure 7 Salt dome models.The reflectivity of the interface at z = 1000 is lower for model 1 (a) than for model 2 (b).The latter generates stronger internal multiples.The red stars on the top of both models indicate the position of the sources.The blue triangles indicate the position of the receivers.

Figure 8
Figure 8 Images of models 1 (a) and 2 (b) using primaries-only migration.With this set-up of acquisition geometry and imaging, it is not possible to image the flank of the salt dome.There is an artefact (false structure) inside the salt dome.

Figure 9
Figure 9 Images of models 1 (a) and 2 (b) using FWM.When using the internal multiples, the flank of the salt dome is imaged.With stronger internal multiples the image further improves.The artefacts inside the salt dome in Figure 8 are no longer present.

Figure 10
Figure 10 Set of PSF obtained with primaries-only migration (a) and FWM (b).

Figure 11 Figure 12 Figure 13
Figure 11 Zoom at the PSF centred at (x, z) = (1000, 640) m for primaries-only migration (a) and FWM (b).The low-resolution axis shows a better resolution for the FWM PSF than for primaries-only migration (c).The wavenumber spectrum of primaries-only migration PSF (d) shows a reduced angle coverage compared with the FWM spectrum (e).

Figure 14
Figure 14 Set of PSFs for the model in Figure 1 obtained with FWM using different sampling intervals.(a) Sources every 125 m.(b) Sources every 500 m.(c) Difference between (a) and (b) amplified 20 times.The computed SNR is 12.12 dB.

Figure 15
Figure 15 PSFs for deblended data (a) (SNR = 7.52 dB), pseudo-deblended data (b) (SNR = −11.33dB) and their differences, (c) and (d), respectively, with respect to the unblended case in Figure 2 (upper-left section).The noise present in Figure (c) and (d) can be attributed to the imaging response of the blending noise.
(a) corresponds to the reflectivity image of model 1, while Figure 8(b) is the reflectivity image of model 2.