## Abstract

We demonstrate for the first time that seismic resonant multiples, usually considered as noise, can be used for super-resolution imaging in the far-field region of sources and receivers. Tests with both synthetic data and field data show that resonant multiples can image reflector boundaries with resolutions more than twice the classical resolution limit. Resolution increases with the order of the resonant multiples. This procedure has important applications in earthquake and exploration seismology, radar, sonar, LIDAR (light detection and ranging), and ultrasound imaging, where the multiples can be used to make high-resolution images.

- seismic resonant multiples
- super-resolution imaging

## INTRODUCTION

In 1873, Abbe recognized a fundamental diffraction limit of optics: whenever an object is imaged by a classical optical lens, features of the object smaller than λ/2 become unrecoverable. Such fine, detailed information is lost because light emerges with these fine features as evanescent energy, which decays exponentially away from the object and is not carried by the propagating waves. Evanescent energy can be used for imaging to achieve a resolution that is higher than the Abbe limit (*1*, *2*). As an example, Lerosey *et al*. (*3*) demonstrated that subwavelength (that is, super-resolution) imaging of the source location could be achieved by refocusing scattered microwave energy from near-field scatterers. They denoted this refocusing operation as a time-reversal mirror (*4*, *5*), which is similar to reverse time migration by exploration geophysicists (*6*–*8*) except that no velocity model is required because the Green’s functions for refocusing are recorded in the data.

Super-resolution imaging is important for many applications in optical imaging, including nanolithography and imaging of live biological specimens (*9*–*11*). For nanolithography, subwavelength focusing is used with ultraviolet light to etch nanoscale transistors, but such light cannot be used with live biological specimens without harming the subject. Recent developments show that super-resolution can be achieved in the far field of both sources and receivers. For example, Huang and Zheludev (*12*) demonstrated that the construction of a mask at the aperture of an optical or plasmonic device will create super-oscillatory constructive interference in the object plane. The super-oscillatory characteristic means that bandlimited functions can oscillate arbitrarily faster than the highest Fourier components they contain (*13*). Numerical simulations show that optical super-resolution imaging can be achieved with a resolution higher than 1/10 λ in the far-field region. However, the imaged spot is characterized by a significantly diminished amplitude and is surrounded by severe sidebands. The energy of the focused spot is often several orders of magnitude less than that of the side lobes, making it difficult to apply this technology to seismic waves.

Analogs of practical super-resolution imaging in seismology are currently limited to a few applications where the objects of interest are in the near field of the sources or receivers (*14*). However, most of the scattering bodies of interest in seismology are in the far field of both the sources and receivers; hence, near-field scattering methods are not applicable. We now show both theoretically and experimentally that subwavelength imaging of reflectors with arbitrary dip angles is possible with resonant multiples. This largely overcomes the limitation of resonant super-resolution to mostly horizontal layers (*15*).

### Resolution limits of resonant multiples

A resonant multiple is defined as a multiple where the upgoing and downgoing ray paths coincide (see Fig. 1). For example, the offset between the source and the receiver in Fig. 1A is zero, such that the *n*th-order free-surface multiple has a travel time described by(1)where *v* is the velocity of the topmost layer and *N* = 0 for primary reflection. The vertical resolution in locating the depth of the interface can be estimated by first taking the derivative of Eq. 1 with respect to *z* to get(2)

Defining the dominant period by *T*_{0}, substituting δ*T* → *T*_{0}/2 and λ = *vT*_{0} into Eq. 1, and rearranging give the vertical resolution limit δ*z* → Δ*z* (*16*)(3)

The limit states that Δ*z* is the minimum separation in depth between two reflectors that are distinguishable from one another in the data. The greater the order *N* of the resonance, the higher the resolution in depth. If the bed is tilted with respect to the horizontal recording surface, then this resolution limit is calculated along the direction perpendicular to the interface in Fig. 1B.

The resolution limit for seismically imaging a dipping reflector is now derived for a resonant multiple with arbitrary order *N*. An example of a third-order (*N* = 3) resonant multiple is depicted in Fig. 1B, where the resonant multiple makes *N* bounces on the free surface, for instance, at *A*_{1} and *A*_{2}. Let the distance from *S* perpendicular to the reflector be . Then, we have and .

The total two-way travel time *T* is given by(4)

Here, *v* is the homogeneous velocity in the upper layer and allows Eq. 4 to subsume the case of a primary reflection (*N* = 0), where *T* = 2*d*/*v*. Perturbing the travel time δ*T* with respect to the perturbation in δ*d* gives(5)

Setting *v*δ*T* → λ/2 gives the *n*th-order resolution limit along the direction perpendicular to the interface(6)where is the resolution coefficient. Setting α = 0 and *N* = 1 gives the vertical resolution limit Δ*d* → Δ*z* = λ/8 for a horizontal reflector, which agrees with Eq. 3 for *N* = 1. Equation 6 also reveals that the resolution limit becomes higher with decreasing interface dip angle.

The resolution limit in Eq. 6 is that for bulk resolution where the entire layer interface in Fig. 1B is shifted by δ*d* along the perpendicular to the interface. For a localized perturbation at *B*_{1} where only the point *B*_{1} is perturbed by the amount δ*d*, the depth resolution limit is inversely proportional to the number of bounces at *B*_{1} that is, Δ*d* = 0.25λ/(*M* cos γ), where γ is the angle of incidence at *B*_{1} and *M* is the number of bounces at *B*_{1}. In Fig. 1B, for example, *M* = 2, even though the ray diagram is that for a third-order multiple. If the rays at the source point *S* and at *A*_{2} are normal to the free surface as in Fig. 1C, then this forms a pure resonance system. Unlike the illustration in Fig. 1B, the higher-order pure resonance rays in Fig. 1C will completely coincide with and *M* will grow with the order *N* of the multiple.

### Imaging resonant free-surface multiples

Relocating a reflection event in a trace to its reflecting point on the interface is known in the geophysics community as migration. For example, if the zero-offset trace *d*(**g**, *t*) only contains a diffraction from a point scatterer at **x**_{o}, then it can be represented by,(7)where *w*(*t*) is the impulsive-like source wavelet (with dominant period *T*_{0}) associated with the point source and geophone at **g** and is the travel time for waves to propagate from **g** to **x**_{o}. Migration answers the following question: Given the trace *d*(**g**, *t*), what are the possible locations of point scatterers that can explain an event that arrives at **g** with travel time ? For a homogeneous medium, the answer is that the scatterers could have been located anywhere along the annulus in Fig. 2A. Here, the radius of the inner annulus centered at **g** is equal to , and the outer annulus radius is . Thus, the thickness of the annulus is Δ*r*^{prim} = λ/2 = 0.5*T*_{o}*v*, which is the difference between primary reflections that reflect off the inner and outer borders of the annulus. This answer, also known as a single-trace migration image, is rather vague because the scatterer could have been located anywhere along the thick annulus. To reduce uncertainty in the answer, the migration images from other traces at **g** ∈ *B* can be added together to obtain the stacked migration image *m*(**x**),(8)where *B* is the set of recording positions along the horizontal recording line; is a preconditioning factor that compensates for geometrical spreading, obliquity, and wavelet stretch; is the Fourier transform of the recorded data *d*(**g**,*t*); and **x** is the trial image point.

Similarly, the single-trace migration image for a first-order resonant event is shown in Fig. 2B. For example, the thickness of the migration annulus for a flat reflector directly below the geophone position is Δ*r*^{mult} = λ/4 = 0.25*T*_{o}*c*, where Δ*r*^{mult}/*c* is the travel time difference between the first-order resonant reflections that scatter from the outer and inner parts of the resonant annulus. Hence, the spatial resolution associated with the resonance migration operator is twice that of the primary migration operator.

The general migration operator for migrating first-order resonant multiples associated with the Fig. 2B model is given by(9)where is the Fourier transform of the source wavelet. Here, τ_{gx} and τ_{xA} are the propagation times from the trial image point **x** to the geophone at **g** and the normal incidence point at **A** on the free surface, respectively (*17*).

The workflow for implementing Eq. 9 is to (i) select any trial image point at **x** in Fig. 2B and compute the imaging time 2(τ_{gx} + τ_{xA}) in Eq. 9 and (ii) use Eq. 9 to migrate the resonant multiple. The above procedure is valid for heterogeneous velocity models where the travel times are computed by ray tracing. Another way of implementing Eq. 9 is to use wavepath migration (*18*, *19*). Its workflow is similar to the previous one, except that the trial image point at **x** has to be on the wavepath from **g** to the image.

## RESULTS

We now present the results of resonant imaging for both synthetic data and field data from two different experiments. The first field data experiment is sonar imaging of the seafloor with 120- and 200-kHz ultrasonic sources. The second experiment is a conventional exploration survey where the airgun source has a frequency content from 5 to 60 Hz.

### Synthetic test

The synthetic experiment is designed to illustrate super-resolution imaging achieved by migrating resonant multiples. Eight point scatterers are placed in a homogeneous background model, and diffractions are generated by Kirchhoff modeling (*20*). The data are recorded on the surface and then migrated by wavepath migration to obtain the images in Fig. 3. The resonant multiples are migrated using Eq. 9 to produce the image in Fig. 3B, which has a finer resolution than the primary image in Fig. 3A. The top four point scatterers are both vertically and horizontally distinguishable from one another in the resonant multiple migration image compared to the blurred primary image.

### Field data tests

The resonant multiples in two field data sets will be tested in their ability to achieve super-resolution imaging. The first data set consists of 120- and 200-kHz echo sounder data that can be used to map out the topography of the seafloor and to detect organisms in the water column such as schools of fish. The experiments were carried out in shallow waters close to the King Abdullah University of Science and Technology (KAUST) by the survey ship shown in Fig. 4A. The second data set is a low-frequency marine survey used to map out deep oil deposits in a sedimentary basin.

**Sonar survey**. Figure 4 (B and C) shows the recorded 120- and 200-kHz sonar profiles. The seafloor mapped out by the 200-kHz echo sounder is used as the ground truth for the 120-kHz echo sounder. That is, the 120-kHz image computed from the resonant multiples in the data in Fig. 4B should achieve about the same resolution as the 200-kHz data in Fig. 4C. Sonar readings were taken 10 times per second as the boat moved at a speed of approximately 1.5 m/s over coral reefs. The primary reflections and first-order resonant multiples were migrated to produce the reflectivity images in Fig. 4D. As predicted from our theory, the image obtained from the 120-kHz resonant multiples has about the same resolution as the image computed from 200-kHz primaries. Also, the resonant multiple image has almost twice the resolution of the 120-kHz primary image. This verifies the theoretical prediction that migrating first-order resonant multiples achieves almost twice the resolution that imaging the primaries does. This assumes no significant loss of high frequencies, which is true in this example.

**Exploration seismic survey**. Figure 5 (A and B) shows two migration images, one from primary reflections and the other from resonant multiples. The data set is collected by a geophysical survey company. The primary reflections are migrated using a conventional Kirchhoff method (*20*), whereas the resonant multiples are moveout-corrected and stacked for a higher signal-to-noise ratio (*17*) and then migrated using Eq. 9. Here, moveout correction aligns the non–zero-offset multiples with zero-offset resonant multiples. As expected, the resonant multiple image in Fig. 5B has about twice the resolution of the primary image in Fig. 5A.

## DISCUSSION

Super-resolution imaging with resonant multiples has possible applications in many fields of science. Figure 6 depicts four possible applications. For 35-GHz radar imaging of cloud structure (www.arm.gov/instruments/instrument.php?id=mmcr), the internal multiples in Fig. 6A might provide higher resolution of cloud structure than the primary reflections. Figure 6B suggests that resonant reflections in full-waveform (*21*) bathymetric LIDAR (light detection and ranging) data (www.nauticalcharts.noaa.gov/hsd/lidar.html) could be used to get better estimates of the water-floor topography or underwater objects than possible with standard bathymetric LIDAR. This might also be true for topographic LIDAR where resonant reflections between the ground and trees can be used to better characterize the underside features of the tree canopy (www.youtube.com/watch?v=HfV7jJgrw4Q). In an urban environment, it might be possible to image the underside of a bridge with resonant multiples. If successful, this opens the possibility of achieving both ground-based and topography LIDAR using just the aerial LIDAR unit. If a mirror (see the pink mirror in Fig. 6B) is attached in front of a ground-based LIDAR, then at least twice the resolution can be achieved with full-waveform LIDAR imaging. This assumes that the topography of the land surface is known. In the ultrasound imaging of the fetus in Fig. 6C, the resonant multiples should provide higher resolution than the primaries if they can be separated from one another. For nondestructive testing and engineering, the surface features of an object can be imaged with super-resolution, as depicted in Fig. 6D.

Increasing the resolution of reflectivity images by migrating multiples is similar to the principle of a Fabry-Perot interferometer, where a light beam is allowed to bounce multiple times between two mirrors in an arm of the interferometer before it interferes with the beam from the other arm. If the arms are slightly different in length, then multiple bounces amplify the slight phase difference between the two different beams, as markedly illustrated with the recent detection of gravitational waves by a Fabry-Perot interferometer (*22*). Increasing the number of resonant bounces in a medium increases the sensitivity of the data to small changes in the medium.

## MATERIALS AND METHODS

### Sonar survey

Acoustic measurements were made using a Simrad EK60 echo sounder (7° beam width) operating at two frequencies (120 and 200 kHz). The pulse duration was 64 μs. Raw data (as shown in fig. S1, A and B) were demeaned and band-pass–filtered such that the peak frequency of each trace in the 200-kHz data was approximately twice that of the 120-kHz data (peak frequencies are shown in fig. S1C). After filtering, the primary reflections and resonant multiples were separated on the basis of their arrival time and normalized individually (as shown in Fig. 4, B and C). This normalization was for the enhancement of the multiple signal. Last, the primary reflections and resonant multiples were imaged using a conventional Kirchhoff method (*20*) and based on Eq. 9 to produce the images of the seafloor in Fig. 4D.

### Exploration seismic survey

The marine data set was collected by a geophysical survey company. There are 334 shots with a shot interval of about 80 m. Each shot has 1126 hydrophones as receivers. The primaries and multiples were separated on the basis of their different arrival times. Then, moveout corrections were applied to primaries and multiples, respectively, to align the non–zero-offset data with the zero-offset data (as shown in fig. S2). After the alignment, stacking across different offsets yields a higher signal-to-noise ratio, zero-offset primary reflections, and resonant multiples (*17*, *20*). Then, the data were filtered by a low-pass filter with a high cutoff frequency of 22 Hz (as shown in fig. S3). Finally, the stacked primaries and resonant multiples were imaged separately to obtain images shown in Fig. 5 (A and B).

## SUMMARY

We have presented a theory for transforming the coherent noise of resonant multiples into signals that can be used for super-resolution imaging. The theory is validated with both synthetic data and field data tests. The limitations to this method are that the resonant multiples must be separated from the primaries, and there must be no significant attenuation in the model; otherwise, the multiples will lose high frequencies compared to the primaries. This procedure has implications for exploration and earthquake seismology, especially if resolution can be more than doubled compared to previous methods. We also believe that it has important applications in radar, sonar, LIDAR, nondestructive testing, and ultrasound imaging, where the noise from multiples can be transformed into high-resolution images.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/2/5/e1501439/DC1

fig. S1. Raw sonar data and their peak frequencies after processing.

fig. S2. Moveout corrections in the common midpoint gathers align the non–zero-offset multiples (covered by yellow shadows) with the zero-offset resonant multiples.

fig. S3. Near-offset marine data before and after moveout correction and stacking.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**

**Funding:**This work was supported by funding from KAUST.

**Author contributions:**G.S., B.G., and Y.H. conceived the theory. A.R. and B.G. collected the sonar data. B.G. implemented the synthetic and field data results. G.S., B.G., and Y.H. wrote the paper.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are available at https://figshare.com/s/7a1d308d49470d430864 or are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2016, The Authors