Diffraction‐limited imaging and beyond – the concept of super resolution ‡

It is well‐known that experimental or numerical backpropagation of waves generated by a point‐source/‐scatterer will refocus on a diffraction‐limited spot with a size not smaller than half the wavelength. More recently, however, super‐resolution techniques have been introduced that apparently can overcome this fundamental physical limit. This paper provides a framework of understanding and analysing both diffraction‐limited imaging as well as super resolution.


I N T R O D U C T I O N
The diffraction limit in wave physics is well-known. De Rosny and  have demonstrated experimentally, using ultrasonic waves, that a time reversed wave backpropagates and refocuses at its initial source. In case of a point source, the spot size can not be smaller than half the wavelength. Both time-reversal and migration within the integral formulation (like Kirchhoff migration, Schneider 1978) can be represented by the same analytical expression (Esmersoy and Oristaglio 1988). This theoretical model is further analysed in this paper employing the general reciprocity theorems introduced by Bojarski (1983). Based on a simple 3D analysis the fundamentals of resolution are revisited and analysed further. In the literature there exist a vast amount of publications discussing the concept of resolution within the context of migration, diffraction tomography and imaging (Devaney 1982;Beylkin, Oristaglio and Miller 1985;Bleistein 1987;Hamran and Lecomte 1993;Tygel, Schleicher and Hubral 1994;Gelius 1995;Pratt, Shin and Hicks 1998;Schuster and Hu 2000;Gelius, Lecomte and Tabti 2002). A significant amount of published work also discusses backpropagation, migration and migration-inversion based on the use of the same anticausal Green's theorem discussed by Bojarski (1983) (Stolt 1978;Clayton and Stolt 1981;Stolt and Weglein 1985;Stolt and Benson 1985;Beylkin 1985;Weglein and Secrest 1990;Weglein et al. 2003).
The resolution analysis presented in the first part of this paper unifies the different ideas referred to above and provides an improved platform to understand the cause of diffractionlimited imaging. It is demonstrated that the monochromatic resolution function consists of both causal and non-causal parts even for ideal acquisition geometries. This is caused by the inherent properties of backpropagation not including the evanescent field contributions. As a consequence, only a diffraction-limited focus can be obtained unless there are ideal acquisition surfaces and an infinite source-frequency band. The diffraction-limited focus can further be visualized as superimposition of a converging and diverging wave as also discussed by Esmersoy and Oristaglio (1988) for a 2D case. This has also been verified experimentally by Rosny and Fink (2002). They introduced in their experiment an acoustic sink, in order to collapse the diverging wavefield and obtained super resolution with the time reversed wave focusing to a spot less than 1/14 th of the wavelength.
In the literature various attempts have been made to obtain images resolved beyond the classical diffraction limit, e.g., super resolution. The main direction of research has been to exploit the evanescent field components (Betzig et al. 1991;Schatzberg and Devaney 1992;Tabib-Azar et al. 1999;Stockman 2004;Grbic, Jiang and Merlin 2008). However, the evanescent waves can (and should) be neglected if the field measurements are performed at greater than two wavelengths from the scatterer(s). Alternatively, improvement of the image resolution of point like targets beyond the diffraction limit can apparently be obtained employing concepts adapted from conventional statistical multiple signal classification (MUSIC) (Schmidt 1986). The basis of this approach is the decomposition of the measurements into two orthogonal domains: signal and noise (nil) spaces. Such a decomposition can be obtained by a singular value decomposition (SVD) of the transfer matrix of the actual experiment. Potential applications of the MUSIC approach within radar imaging (Odendaal, Barnard and Pistorius 1994), microwave breast imaging (Scholz 2002) as well as within acoustic imaging (Lehman and Devaney 2003) have been published. Within an acoustic formulation, Lehman and Devaney (2003) employed the framework of time-reversal (Prada et al. 2002) to establish a super-resolution algorithm. In this paper a simpler and more direct approach is being used to analyse the concept of super resolution. It is demonstrated that by utilizing the null-space solutions of the wave problem, super resolution is apparently obtained since such solutions can give an extremely well localization of point-source targets compared to migration. The original work of Lehman and Devaney (2003) did not discuss the robustness of this super-resolution technique with respect to both acquisition aperture as well as noise. In this paper it is demonstrated that the method can be rather sensitive to noise for limited acquisition surfaces when monochromatic data are employed.

F U N D A M E N T A L S O F R E S O L U T I O N
Consider a closed surface S defining a volume V of space as shown in Fig. 1 and assume that receivers are distributed across the surface. A point source is located inside V at a position r 0 and illuminates the scatterers embedded in the (possible) non-uniform background medium. As shown in Appendix A (cf., equation (A6)), the monochromatic scattered wavefield associated with a compact perturbation is now given by the Lippman-Schwinger equation (Newton 1982): where p represents the total wavefield, G 0 is the background Green's function and k 0 and α are respectively the background wavenumber and the scattering potential defined as: with c and c 0 representing the space-variant velocities of respectively the scatterers and the background medium. Assume that the scattered waves are recorded by receivers placed on the surface S and backpropagated employing an integral formulation of migration, as described by the Kirchhoff integral (Schneider 1978;Wiggins 1984;Langenberg 1987;Esmersoy and Oristaglio 1988;Schleicher, Tygel and Hubral 2007). The backpropagated scattered wavefield can then be written explicitly as: Note that equation (3) represents an approximation of equation (A5) derived in Appendix A, where the second volume integral in equation (A5) has not been accounted for in equation (3). This leads to the non-causal behaviour of backpropagation as discussed below. For simplicity, often the high-frequency and far-field version of equation (3) is employed in imaging (i.e., migration) The surface integral defining the backpropagation of the scattered field can be replaced by an alternative volume integral formulation by combining equations (1) and (3) (Esmersoy and Oristaglio 1988): Comparison between equations (1) and (5) shows that the background Green's function G 0 ( r , r , ω) is now replaced by the backpropagation kernel B ( r , r , ω). This kernel has a simple physical interpretation: it represents backpropagation of the wavefield associated with a (secondary) point source located at r (cf., equation (3)) and is therefore closely related to the resolution function as discussed in the following. Ideally, the backpropagated field as given by equation (5) should resemble the scattered field as given by equation (1). However, due to the special properties of the backpropagation operation and consequently also kernel B this will only be fulfilled in case of an ideal acquisition aperture (e.g., receivers uniformly distributed over the closed surface S) and an infinite frequency band. To realize this, the introduction of the concept of a resolution function is very useful. Employing a U/D type of imaging condition (e.g., assuming a scatterer exists where the first-arrival of the downgoing wave is time-coincident with the upgoing wave; Claerbout 1971) an estimate of the scattering potential can be obtained from the equation (integrating over the available frequency band ω) where RF is the resolution function or point-spread function of the image when a 't = 0' image condition is applied (cf., equation (6)). Since kernel B is 'singular' when r = r , it is feasible to set p( r , ω)/ p( r, ω) ∼ = 1 in the expression for the resolution function, e.g.: In further analysis, it gives additional insight to also make use of the alternative expressions for RF ( r , r ) and B ( r , r , ω) as derived in Appendix B assuming a locally homogeneous background medium with k 0 = ω/c 0 (cf. equations (B11) and (B12) employing high-frequency and far-field assumptions): and In equation (8) the scattering wavenumber K spans out the Fourier transformed model space and the unity vectorυ defines its direction (cf., equation (B10)): For a given point r inside the scattering structure the unity vectorυ connects this point with a general receiver at the acquisition surface through a ray as shown in Fig. B1 in Appendix B. In case of a complete coverage of receivers, this unity vector spans the surface of a unity sphere surrounding r Gelius et al. 2002).

D I F F R A C T I O N -L I M I T E D I M A G I N G ( I D E A L A P E R T U R E A N D M O N O C H R O M A T I C )
According to equation (7) the governing part of the resolution function is given by the kernel B. For a monochromatic case, understanding the role of B also then automatically gives a good idea of the resolving power of imaging based on backpropagation as governed by equation (5) (and also migration based on integral formulations in general).
In case of an ideal aperture it is shown in Appendix A (cf., equation (A7)) that kernel B takes this special form: which represents a superposition of a time-advanced and timeretarded Green's function. This implies that kernel B has a causal and an anticausal part and as a consequence backpropagation of the recorded field associated with a point source as described by B will give a diffraction-limited focus. A similar result has been obtained by Esmersoy and Oristaglio (1988) by analysing 2D wave propagation. Considering a locally homogeneous background, it follows from equation (11): which gives a focus described by the sinc-function. In equation (12) k 0 is the wavenumber, c 0 is the velocity of the locally uniform background medium and λ 0 is the corresponding wavelength. Assuming that the size of the spherical focal spot is mainly defined by the main lobe of the radial sincfunction, its diameter d can be approximated from: The result in equation (13) is analogous to the focused beam size limit of imaging optics as determined by the diffraction of light. For an axial symmetric optical system the formula of the resolution limit can be obtained from the classical diffraction theory for electromagnetic waves, i.e., Rayleigh criterion (Born and Wolf 1999): where is the focused beam size. Moreover, NA is the socalled numerical aperture of a lens that characterizes the range of angles over which the system can accept or emit light. In the theoretical limit NA → 1 then → 0.61λ 0 . Within optics, the diffraction limit can also be explained by a simple diffraction theory for a single slit. When a parallel and monochromatic beam is incident on a slit the wave is diffracted and the amplitude of the electromagnetic wave far from the slit can again be described by a sinc-function (Born and Wolf 1999).
By introducing a finite frequency band, backpropagation of the kernel B as given by equation (5) can be better visualized. Consider now the simple example shown in Fig. 2. The wavefields generated by a point source embedded in a homogenous medium were recorded at receivers evenly distributed over a sphere (a total number of 1600 measurement points employed here and with the source placed in the centre). In the calculations the homogenous background medium was assigned a velocity of 2000 m/s and the source was represented by a band-limited pulse defined by a Ricker wavelet with a centre frequency of 20 Hz. Backpropagation in the time domain is now given by Fourier synthesis over the frequency band (here limited to be between 0-50 Hz): with B( r , r , ω) computed from equation (5). Moreover, the focus is represented by B( r , r , t = 0). Figure 3 shows snapshots of the backpropagated field associated with the kernel (note that each snapshot has been individually scaled so only a qualitative description is obtained). It can be easily seen from Fig. 3 the presence of non-causal contributions in snapshots (a) and (b). These contributions are non-physical since they appear at negative traveltimes where the physical (scattered) wavefield does not exist. The diffraction-limited focus given by snapshot (c) is caused by interference between converging and (polarity reversed) diverging waves as described by equation (11). The size of the focal spot in Fig. 3(c) is in the order of half of the centre wavelength as expected.
The non-causal nature of backpropagation kernel B can be further understood by assuming again a locally homogeneous background medium. The Green's function G 0 ( r , r , ω) appearing in equation (11) can then be written explicitly as: In equation (16) this Green's function is also symbolically decomposed into its homogeneous (G 0,H ) and inhomogeneous (G 0,I H ) parts. The homogenous part describes the propagating waves whereas the inhomogeneous part accounts for the evanescent field. In Appendix C it is demonstrated that the backpropagation kernel B in case of an ideal aperture and a uniform medium is closely related to the homogenous part G 0,H of the Green's function (cf., equations (C10)-(C13)), with the only difference being the limit of integration. A similar result has been obtained by Esmersoy and Oristaglio (1988) for the 2D case. This implies that backpropagation as given by equations (3) and (5) essentially only accounts for the propagating waves and does not account for the evanescent fields as also discussed by others (for example, Wapenaar and Haimé 1990). This causes again the backpropagation operation to be non-causal. As also discussed in Appendix C, the inhomogeneous part of the Green's function associated with a point-source is closely related to the real part of the same Green's function (cf., equation (C13)). The structure of the inhomogeneous Green's function will therefore be quite different: i.e., it will 'peak' (be singular) at the correct point-source location for every frequency. The evanescent field has therefore the potential of giving 'super resolution' (Betzig et al. 1991;Schatzberg and Devaney 1992;Tabib-Azar et al. 1999;Stockman 2004;Grbic et al. 2008).

I M A G I N G W I T H A F R E Q U E N C Y B A N D ( I D E A L A P E R T U R E )
As discussed in the previous section, the backpropagation kernel B plays the role of a monochromatic resolution function. From the previous analysis of B the phenomenon of diffraction-limited imaging could be easily understood. In case of imaging with a frequency band the backpropagation kernel B needs to be recomputed for the relevant range of frequencies before the corresponding resolution function can be obtained from equation (7). By increasing the frequency band of the illuminating source pulse an improved point-spread function (focus) can be obtained as shown in Fig. 4.
In the limit of infinite frequencies the focus will be perfectly resolved (i.e., ideal point-spread or resolution function). This can be seen by combining equations (7) and (11) assuming a locally homogeneous background (also introducing spherical In order to arrive at the final result in equation (17) the relationship between a delta-function, expressed by respectively  spherical and Cartesian coordinates, has been employed as given by Bracewell (1999). Combination of equations (6) and (17) now gives α( r) = α( r ), e.g., perfectly resolved scatterers as expected.
The same result follows from use of the alternative expression of the resolution function as given by equation (8). In case of a complete coverage of receivers, the unity vectorυ in equation (10) spans the surface of a unity sphere surrounding r . Moreover, it follows from the same equation that the length of the scattering wavenumber K is proportional to the frequency and that its direction is defined byυ. Hence introducing an infinite frequency band in addition to an ideal aperture implies that K samples the complete Fourier space of the model locally around r , e.g., the resolution function takes the limit:

T H E E F F E C T O F A P E R T U R E O N I M A G E R E S O L U T I O N
It has been demonstrated that in case of an ideal aperture and a monochromatic signal, the corresponding resolution function is diffraction-limited. Only by introducing an infinite band of frequencies can an ideal point-spread function be obtained.
In case of a limited aperture, the image resolution will be even more distorted and the backpropagation kernel B can be decomposed into two parts: where B H corresponds to an ideal aperture (cf. equation (11)) and B A represents the additional effect of a possible limited acquisition geometry. By combining equations (5) and (11) the aperture part of backpropagation kernel B can be expressed as (also making use of equation (B12)): In equation (20) the vector k 0υ scaled by ω −1 represents the slowness vector associated with the ray from a diffractor at r to the receiver. The (half) angle that is made by this slowness vector and the slowness vector from the source at the diffractor (here being a fixed source) is called the diffraction angle (Gelius 2002) (cf., Fig. B1 where ϑ represents symbolically the scattering angle in case of 2D). Correspondingly, this diffraction angle represents a common definition of the aperture.
To visualize the additional effect of a limited aperture, the experiment described in Fig. 2 is repeated but this time receivers are only placed along the right part of a ring in the vertical slice through the surface of the sphere. The corresponding snapshots of kernel B show a more complex wave focus (caused by treating a 3D problem as a two-dimensional one). Again the focus (Fig. 5b) is formed by the interaction of converging and diverging waves and is more distorted than the focus obtained in the ideal aperture case (cf., Fig. 3c).
The effect of both a limited aperture as well as a finite frequency band will distort the final point-spread function. Assuming a locally homogenous background model, the resolution function for the scattering example considered here is approximated well by equation (8).

T H E C O N C E P T O F S U P E R R E S O L U T I O N
Various attempts exist in the literature aiming to remove the resolution artefacts caused by the limited bandwidth as well as the finite acquisition aperture (e.g., deconvolving the space-variant non-ideal resolution function). One approach is to compute the space-variant resolution functions employing ray-based techniques and then correct for it by deconvolution (Gelius et al. 2002;Sjoeberg, Gelius and Lecomte 2003). Alternatively, a least-squares representation of migration or imaging (LSM) can be introduced (Nemeth, Wu and Schuster 1999;Duquet, Marfurt and Dellinger 2000;Kuhl and Sacchi 2003). Within this formulation the Hessian matrix represents the resolution function and is inverted for in an iterative manner. However, this type of algorithm is computationally intensive and also needs careful stabilization and conditioning. An alternative implementation of LSM denoted migration deconvolution (MD) has therefore been introduced, where the exact Hessian is replaced by a less expensive approximation (Hu, Schuster and Valasek 2001;Yu et al. 2006). Recently, full-waveform inversion (FWI) has also been given much attention again because of the rapid developments in computing power (Virieux and Operto 2009). FWI is based on a fullwaveform forward-modelling engine and a local differential approach in which the gradient and the Hessian operators are efficiently estimated.
Additional improvements in resolution can be obtained if also the evanescent energy can be made use of. If that is the case super-resolved imaging can be feasible (Betzig et al. 1991;Schatzberg and Devaney 1992;Tabib-Azar et al. 1999;Stockman 2004;Grbic et al. 2008). However, this is not a practical approach in case of seismic data since the actual field measurements are performed at greater than two wavelengths from the scatterer(s).
Recently, there has been a growing interest in utilizing the diffracted wave energy to improve velocity estimates as well as enhancing the resolution of finer details associated with faults and pinch outs (Fomel, Landa and Taner 2007;Moser and Howard 2008). Diffractions associated with pointlike features can apparently be imaged beyond the classical resolution limit by employing a MUSIC (MUltiple SIgnal Classification) type of algorithm (Lehman and Devaney 2003;Gelius 2009). The key to such super resolution is to decompose the measurements into a signal space and a complementary noise (nil) space. Such a decomposition is obtained by employing a singular value decomposition (SVD) analysis of the data set. The orthogonality between the two spaces can be used to construct 'images' based on the nilspace solutions. Such nil-space based 'images' have the potential of giving extreme precise localizations of the scatterers and show characteristics similar to the imaging of evanescent energy.
The original work of Lehman and Devaney (2003) derived the main results within the frame work of time reversal. In this paper a more direct approach is employed that provides a simpler analysis of the problem. Moreover, the work of Lehman and Devaney (2003) did not discuss the robustness of this super-resolution technique with respect to both acquisition aperture as well as noise. In this paper it is shown that in case of random noise the super-resolution power may break down in case of limited apertures when employing monochromatic data.

Figure 6
Multi-source multi-receiver experiment involving a cluster of D scatterers.

Basic algorithm of super resolution
In this section the basic concept of super resolution is introduced. The derivation presented here is somewhat simpler than in the original work of Lehman and Devaney (2003) by avoiding the time-reversal analysis. Consider a multi-source multi-receiver experiment (cf., Fig. 6) with a total of D scatterers (scattering strength d i , i = 1,2, . . D) embedded in a possible inhomogeneous background. In case of possible transient signals it is assumed that a temporal Fourier transform has been applied. Hence, the following formulation is derived for a monochromatic case. Let N s and N r represent respectively the total number of sources and receivers and the monochromatic signalS (vector with dimension N s ) represent the transmitted signal from the source array. Introducing the complex transfer N r xN s matrix K, gives the following relationship (assuming a noise-free case): whereR represents the data measured at the receiver array. Assume that the rank of the matrix is D < min(N r, N s ) with D representing the total number of scatterers. This implies a physical experiment with D linear scattering contributions, e.g., Born type with no interactions between the scatterers.
Introducing a singular value decomposition (SVD) of the system matrix K in equation (21) now gives (Vaccaro 1994): where subscript s means signal space and ⊥ denotes the complementary space (nil space). Moreover, superscript H means complex conjugated transposed. Due to the orthogonal signal and nil spaces in equation (22) the following holds: Also due to inherent normalization: Assume now that the scatterers are fully resolved (e.g., ideal array point-spread functions with respect to both source and receiver side), which mathematically can be stated as the conditions whereḠ 0s ( r i ) andḠ 0r ( r i )represent a background Greens function (column) vector with respect to respectively the source and the receiver array that focuses at a point scatterer located at the position r i . Due to equation (25) the same Greens function vectors (normalized versions) will form the singular functions associated with the scatterers. Hence, the left-and right-singular matrices take the form: It is known that the singular vectors output from a SVD analysis of a complex-valued matrix system will be non-unique by an arbitrary phase. The phase angles θ i in equation (26) represent symbolically this non-uniqueness. The diagonal matrix s in equation (22) now contains the corresponding (amplitude weighted) scattering strengths along its diagonal. Based on these observations a monochromatic MUSIC type of pseudo-spectrum operator can be introduced that peaks at the true target locations (the asterisk implies complex conjugating and the superscript T means transposed): Introduce now the projection matrices for the nil space with respect to both the receiver and source array side: These matrices have the following characteristics: Projection matrices for the signal space can be introduced correspondingly: which show the complementary characteristics of those in equation (28). Finally, by combining equations (27) and (28) the alternative form of the MUSIC pseudo-spectrum operator can be established: With reference to the previous discussion about evanescent waves, the use of the noise-space in this way gives a result that is similar to the effect of including evanescent energy: e.g., it peaks at the location of the scatterer(s). However, this does not imply that the nil-space singular functions represent the true evanescent waves; they only give an overall effect in the 'image' that resembles well the one obtained from employing the true evanescent fields. In case of strong scatterers there will in general no longer exist a one to one mapping between the non-zero eigenvalues and the scatterers. Instead, a given eigenvalue (and corresponding singular functions) will now represent a combination of several scatterers. Even if a SVD analysis can be employed the transition between signal and noise space will no longer be so sharp and also the extreme localization capability of the nil-solutions will degrade.
To illustrate the power of the super-resolution technique the simple scattering example shown in Fig. 7(a) was considered. It involves a fairly ideal acquisition geometry represented by a ring surrounding two scatterers. Coincident sources and receivers (transceiver system) were distributed evenly along the ring (half centre-wavelength separation). Figure 7(b) shows the array resolution function with respect to one of the scatterers (combined effect of source and receiver array and a single frequency of 20 Hz). Such an array resolution function is computed from an expression analogue to the one in equation (25) but this time employing a fixed position r j corresponding to the true location of the scatterer considered and with the other position vector r i running through all relevant image points. This array point spread function (PSF) is seen to peak well at the target location even if it is not completely fulfilling the condition in equation (25), e.g., it is not being zero at the location of the other scatterer. A similar result is obtained for the PSF corresponding to the other point target.
Data were generated in the time domain by employing a Ricker wavelet with a centre frequency of 20 Hz and the Foldy-Lax theory (Green and Lumme 2005) to take into account interactions between the two scatterers. The separation between the two scatterers was one quarter of the centre wavelength. The actual super-resolution computation was carried out in the frequency domain after carrying out a Fourier transform. Figure 7(c) shows the super resolved result with the two scatterers being well resolved. Note that this result is obtained employing a single frequency of 20 Hz. As a comparison, Fig. 7(d) shows the result obtained employing a Kirchhoff type of migration using the same single-frequency data. In this case the two scatterers can not be distinguished because of the diffraction limit.
In the next example a rather limited acquisition geometry was employed as shown in Fig. 8(a), e.g., involving only four sources and three receivers organized in a random way. The corresponding monochromatic (20 Hz) array resolution function is shown in Fig. 8(b) with respect to one of the scatterers. Again a similar result is obtained for the array PSF associated with the second point target. Comparison between the array resolution function in Fig. 8(b) and the ideal criterion given by equation (25) indicates clearly that this latter assumption is violated. However, it turns out that the superresolution technique is more robust when used in practice as long as the array PSF of a given scatterer peaks at its true location and shows significantly lower values at the location of the other scatterer(s) as shown in Fig. 8 Fig. 8(c), which shows the result obtained employing the superresolution technique (monochromatic data corresponding to the centre frequency of 20 Hz). The quality of this 'image' is comparable to the one in Fig. 7(c) where a rather ideal acquisition geometry was applied. Finally, employing the same monochromatic data as input to the Kirchhoff migration the image in Fig. 8(d) was arrived upon. As expected, due to the sparse acquisition geometry and also the use of monochromatic data only, the reconstruction of the two scatterers is heavily distorted.

Including random noise
In the original paper of Lehman and Devaney (2003) the issue of noise was not addressed. As demonstrated in Appendix D (cf., equation (D24)), in case of random noise the expected value of the MUSIC pseudo-spectrum operator now takes the form (σ 2 being the variance of the noise): To obtain an improved understanding of equation (32) assume for simplicity the same singular value η for all scatterers. This assumption leads to the simplified form: where It follows from equation (33) that the effect of noise is to degrade the projection matrices of the nil space caused by leakage from the signal-space projection matrices. Revisit now the controlled-data example shown in Fig. 7, which represents a fairly ideal acquisition case. Adding random noise (corresponding to S/N of 20 dB) to the data gave still a super-resolved result comparable to that in Fig. 7(c) (not shown here). This is due to the fact that the nil-space projection matrices have a much larger magnitude than the corresponding signal-space projection matrices. This is illustrated in Fig. 9 showing the magnitude of the receiver side projection matrices with respect to both signal-and nil-space (source side projection matrices are here the same because of the acquisition symmetry). Note that because of equations (28)-(30) the following relations hold: P ⊥,r + P s,r = I Nr xNr , and P ⊥,s + P s,s = I Ns xNs .
Since the projection matrix of the signal space has a much smaller magnitude than the corresponding nil-space projection matrix in Fig. 9, the latter is close to a diagonal matrix (cf. equation (35)). Next, the controlled-data experiment from Fig. 8(a) was revisited. Again random noise was added to the time-domain data corresponding to a S/N of 20 dB. Figure 10(a) shows the output from the super-resolution technique (monochromatic input at 20 Hz) demonstrating its large sensitiveness to noise. The 'image' is no longer able to resolve the positions of the two scatterers. As a reference, the corresponding result employing Kirchhoff migration to the same monochromatic noisy data is shown in Fig. 10(b). As expected, migration is much more robust with respect to random noise (apparently no significant differences between Figs 8d and 10b). Figure 10(c,d) shows plots of respectively the nil-and signalspace projection matrices (both receiver and source side).
Compared to the case in Fig. 9, the nil-and signal-space matrices are now much more similar in both pattern and magnitude level. This seems to be an effect of distributing the sources and receivers in a more random way. Because of these similarities one may assume that the projection matrix of the nil space is less distorted in case of noise. Consider now a band of frequencies (covering the most energetic part of the signal spectrum) and make the following assumptions: i) by varying the frequency the false events observed for a single frequency (cf., Fig. 10a) become non-stationary and can add more or less out if a proper bandwidth is employed and ii) it is reasonable that for one or more frequencies within the band considered the random noise is weak and the super-resolution capability has been more or less preserved. This leads to the introduction of a multi-frequency pseudo-spectrum operator: with P ω,MUSIC given by equation (31). By using frequencies between 15-45 Hz and the algorithm in equation (37), the result shown in Fig. 10(e) was obtained. Compared to the image in Fig. 10(a) a significant improvement has been achieved, where the two scatterers are now quite well resolved. Correspondingly, Fig. 10(d) shows the result obtained employing migration and the same frequency band. The latter reconstruction is seen to be much worse and the scatterers are not resolved.
In the last controlled data example, again a limited aperture was considered (cf., Fig. 11a) but unlike the geometry used in Fig. 10 the sources and receivers are now distributed in a well organized pattern (e.g., line geometry). In case of noise-free data super resolution can be easily achieved also for this set up (not shown here). However, when adding random noise to the time-domain data corresponding to a S/N of 20 dB the degraded result shown in Fig. 11(c) was obtained (monochromatic input at 20 Hz). From Fig.11(c) it now seems that only one scatterer exists for this model. The result of Kirchhoff migration applied to the same monochromatic noisy data is shown in Fig. 11(d). As expected, migration can not resolve the features at all. Figure 11(b) shows plots of both the nil-and signal-space projection matrices with respect to the receiver side (for this particular example due to its symmetry the source side projection matrices are just the same). Comparison between Figs 11(b) and 9 shows that the overall magnitude of the signal-space projection matrix is now larger relative to the corresponding nil-space projection matrix. This implies an in-creased sensitivity to noise that has also been reflected in the degraded result shown in Fig. 11(c). Moreover, since the projection matrices in Fig. 11(b) are still much closer in similarity to those in Fig. 9 than those in Fig. 10(c,d), the use of a frequency band will most likely not improve the super-resolution capability. This is also demonstrated in the numerical experiment where Fig. 11(e) shows the result obtained using frequencies between 15-45 Hz. Comparison with Fig. 11(c) shows that essentially no improvement has been obtained. Finally, Fig. 11(f) shows the result obtained employing migration and the same frequency band. As expected Kirchhoff migration can not resolve the two scatterers even in the case of multifrequency data.

C O N C L U D I N G R E M A R K S
This paper has provided a framework of understanding and analysing both diffraction-limited imaging as well as super resolution. Extrapolation of measurement data employing both the principles of migration and time-reversal give the same estimate of the scattered wavefields. In case of a point scatterer it is shown explicitly that the band-limited character of the focus is caused by a superimposition of converging and diverging waves. Hence, the backpropagation operation is non-causal by nature since it is based on the anti-causal Green's theorem. Imaging beyond the classical diffraction limit has been proposed in the literature by making use of the evanescent field contributions. However, this approach is not practical in case of seismic imaging in general since the evanescent waves are so weak because of attenuation, they are masked by the noise. However, point-diffracted data can apparently be super resolved by making use of the null-space solutions. This concept applies equally to a single or a collection of scatterers as long as ideally there are no strong interactions between them (Born type). In this paper, it was demonstrated that two scatterers with strong interaction still could be super resolved. Hence, the Born assumption can be somewhat relaxed if there are a small number of scatterers. However, the apparently super-resolved focus obtained by this technique has to be interpreted as an extreme localization rather than a quantitative image (the value of the peak is arbitrary). A possible extension of the idea to image segments (i.e., line of scatterers) should be further investigated. The super-resolution technique is obviously not able to form complete and reliable images of the subsurface. However, by separating diffractions from reflections and also by further identifying a local target region (for  example, diffractions associated with local faults or alternatively local changes computed from residual time-lapse data) then super-resolution methods have the potential of adding more details to the big picture.

A C K N O W L E D G E M E N T S
This work has been partly funded by the Norwegian Science Foundation (NFR). The authors are thankful to two anonymous reviewers and an associate editor for many helpful suggestions to improve the manuscript.

A P P E N D I X A : B A S I C F O R M U L A T I O N A N D C O N C E P T S
Consider a closed surface S defining a volume V of space as shown in Fig. 1 and assume that receivers are evenly distributed across the surface. A point source is located inside V at a position r 0 and illuminates the scatterers embedded in the (possible) non-uniform background medium. The total wavefield p can be described by the scalar inhomogeneous time-reduced Helmholtz wave equation: In equation (A1) k 0 and α are respectively the background wavenumber and the scattering potential defined as k 0 ( r , ω) = ω/c 0 ( r ) and α( r ) = [(c 2 0 ( r )/c 2 ( r )) − 1]. The sound velocities c and c 0 represent the space-variant velocities of respectively the scatterers and the background medium.
In the following, rigorous solutions of equation (A1) will be established employing the technique of Bojarski (1983), which involves Kirchhoff integration of the considered waveequation over volume V with the time-retarded and timeadvanced Green's function, respectively.
The governing equation for the time-retarded Green's function G 0 of the background medium is given by Taking the complex conjugate of both sides of equation (A2) yields the corresponding equation for the time-advanced Note that equation (A3) is only valid within a medium with no attenuation. If the medium is visco-acoustic then the wavenumber k 0 must also be complex conjugated. Making use of the reciprocity theorems of Bojarski (1983) gives two alternative solutions for the total field (an ideal aperture): and In equations (A4) and (A5) the notation ∂/∂n implies normal derivative. The contribution from the surface integral in equation (A4) will vanish if the surface surrounding the volume of interest is assumed to be at infinity and radiation conditions are employed (Bojarski 1983;Weglein and Secrest 1990;Ramirez and Weglein 2009). Equation (A4) is also only true if the observation point r is explicitly within the volume V (Ramirez and Weglein 2009). However, the surface integral in equation (A5) will not cancel out under similar circumstances. Based on the assumptions discussed above, an expression for the scattered field p s can now be constructed from equation (A4): which is the so-called Lippman-Schwinger equation (Newton 1982). This equation represents an all-time all-space solution.
Hence, the volume represents the whole space and no surface integral exists. If the volume is considered to be finite without including the surface integral, then the assumption is that the perturbation is compact and contained within V and that there are no sources or scatterers outside the volume (Weglein and Secrest 1990). Moreover, a combination of equations (A4) and (A5) can also give a very useful expression. Since the left-hand side of these equations are equal, by introducing α = 0 (e.g., no scatterer) in the same two equations and setting them equal gives the result (also replacing r 0 with r without losing generality): This latter equation can be interpreted as backpropagation of the waves generated from a point source at r and measured at the closed surface S. It can also be easily shown that equation (A7) can serve as a link both to time-reversal and seismic interferometry. By employing the principle of reciprocity of the Green's function, equation (A7) can be written on the alternative form: which represents the time-reversal analogy (Fink 1999) of equation (A7) (i.e., now with a point-source placed at r ). Equation (A8) can be interpreted as follows: the wavefield from a point source at r is measured along surface S. These measurements are then time-reversed (which corresponds to the operation of complex conjugating in the frequency domain) and forward propagated to obtain a focus. Comparison between equations (A7) and (A8) shows that the basic operation of backpropagation within a migration formulation or a time-reversal formulation gives the same diffraction-limited image. Also by employing reciprocity, equation (A7) can be transformed to another alternative form with a new physical interpretation: Equation (A9) serves as the starting point of seismic interferometry (Wapenaar 2004) and the closely-related concept of virtual-source imaging (Bakulin and Calvert 2006). This can be more easily seen if the high-frequency and far-field version of equation (A9) is employed: Realize first that the operation under the integral represents a cross-correlation in the frequency domain. The physical interpretation is now as follows: the contributions from active sources distributed over the surface S are measured at two receiver positions r and r . Cross-correlation of these data (left-hand side) transforms one of the receivers to a virtual source (cf., right-hand side) as discussed by Bakulin and Calvert (2006).
Alternatively, by assuming random underground sources (associated with reservoir production as an example), equation (A10) can be used to transform passive seismic surface recordings to active experiments (Wapenaar 2004). This is known as seismic interferometry in the literature.

A P P E N D I X B : T H E R E S O L U T I O N F U N C T I O N I N C A S E O F A G E N E R A L A P E R T U R E A N D A L I M I T E D F R E Q U E N C Y B A N D
The resolution function is defined in equation (7) and valid for an available frequency band ω: A limited aperture implies that the inner integration is only taken over a part of the closed surface S (e.g., limited coverage of receivers). Introduce the high-frequency approximation and represent the Green's functions by: In equation (B2) a Taylor expansion of the phase term has been introduced since points r are assumed close to r . Moreover, the amplitude factor is assumed to be real, e.g., a non-viscous medium. Next, assume far-field configuration and let S represent the wavefront surface at infinity due to a source located at r . Hence, the normal derivative terms in equation (B1) can be evaluated on the wavefront surface S by using (c 0 being the velocity along the surface): Combination of equations (B1)-(B3) gives the far-field and high-frequency approximation of the kernel: Assume that a unity sphere is surrounding the point r (e.g., |υ| = 1 in Fig. B1). Conservation of energy within the ray tube gives the relationship (Keller 1987): with the spherical coordinate system defined bŷ υ = sin υ 1 cos υ 2î + sin υ 1 sin υ 2ĵ + cos υ 1k , Assuming a homogenous medium within the infinitesimal sphere it follows that (spherical wave): and also that Combination of equations (B4), (B5), (B7) and (B8) gives (also making use of the fact that ∇τ ( r , r ) = − 1 In this derivation the following assumptions have been employed: r A locally homogeneous volume surrounds the point-source located at r .
r All points r considered are close to r (in accordance with equation (B2)). Introduce now the coordinate transformation (υ, ω) → K where: is the scattering wavenumber vector spanning out the Fourier space of the model (cf., equation (B11)). After this transformation, equation (B9) takes the form For a given point r inside the scattering structure the unity vectorυ connects this point with a general receiver at the acquisition surface through a ray as shown in Fig. B1. Moreover, it follows from equation (B10) that the length of the scattering wavenumber K is proportional to the frequency and that its direction is defined byυ. For a limited aperture and a limited frequency band K will therefore span out a subspace K. Comparison between equations (B9) and (7) also gives this expression for the backpropagation kernel:

A P P E N D I X C : T H E H O M O G E N E O U S P A R T O F T H E G R E E N ' S F U N C T I O N A S S O C I A T E D W I T H A P O I N T -S O U R C E ( S P H E R I C A L W A V E A N D U N I F O R M B A C K G R O U N D )
In case of a uniform medium the Green's function associated with a point-source is given by a spherical wave: Angular-spectrum or plane-wave decomposition of this Green's function can be carried out by employing the Weyl expansion (Mandel and Wolf 1995): Introduce now a transformation from wavenumber coordinates (k x , k y ) to spherical coordinates (υ 1 , υ 2 ) (assuming a sphere with fixed radius k 0 ): k x = k 0 sin υ 1 cos υ 2 , k y = k 0 sin υ 1 sin υ 2 , (k z = k 0 cos υ 1 ), which corresponds to the Jacobian dk x dk y = k 0 k z sin υ 1 dυ 1 dυ 2 .
Also introduce the position vector Combination of equations (C6), (C8) and (C9) gives the compact form (considering upper half-space where z is positive, lower half-space just symmetrical) In case of ideal aperture, it follows from equation (B12) in Appendix B that (replacingr −r withR without any loss of generality): Direct comparison between equations (C10) and (C11) now reveals that the backpropagation kernel B is closely related to the homogenous part of the Green's function. The only difference in their expressions is the integration limit of the angle υ 1 . Evaluating equations (C10) and (C11) in the plane (x,y) (e.g., assuming z = 0) gives the simplified result: B(R, ω) = G 0 (R, ω) − G * 0 (R, ω) = 2iIm G 0 (R, ω) = 2G 0,H (R, ω).

A P P E N D I X D : M U S I C P S E U D O -S P E C T R U M I N C A S E O F R A N D O M N O I S E
Vaccaro (1994) introduced a second-order perturbation expansion for the SVD. The derivation presented here follows his outline but specializes to a first-order analysis for simplicity. Limiting to first-order perturbations should be satisfactorily enough to describe the case of weaker noise. Start now with the noise corrupted version of the system matrix (Vaccaro 1994): and its corresponding SVD-form: K will generally have full rank but is 'close' to a matrix of rank D if the noise N is small. In the following this is assumed. Introduce now for the perturbed left-hand and right-hand singular vector matrices (linear or first-order perturbation assumed) Due to orthogonality between the signal and nil space the following condition holds: Combination between equations (D3) and (D4) then gives In order to arrive at the result in equation (D7) one has to make use of equations (23) and (24). Moreover, assuming a firstorder perturbation analysis the second-order term Q H Q is neglected. The same check can be made on the other perturbed singular matrices giving the same result. By analogy with equation (31) the monochromatic MUSIC operator now takes the form where the nil-space projection matrices are defined as and The issue is now to find explicit expressions for the first-order perturbation matrices Q and L. Consider first the matrix Q. The starting point is equation (D1) and multiply now this equation with the factorŨ H ⊥ on both sides, e.g.: Employing equations (22) Hence, the expected value of the MUSIC pseudo-spectrum operator in equation (D8) is now: