Reflection full-waveform inversion (RFWI) can recover the low-wavenumber components of the velocity model along with the reflection wavepaths. However, this requires an expensive least-squares reverse time migration (LSRTM) to construct the perturbation image that can still suffer from cycle-skipping problems. As an inexpensive alternative to LSRTM, we use migration deconvolution (MD) with RFWI. To mitigate cycle-skipping problems, we develop a multiscale reflection phase inversion (MRPI) strategy that boosts the low-frequency data and should only explain the phase information in the recorded data, not its magnitude spectrum. We also use the rolling-offset strategy that gradually extends the offset range of data with an increasing number of iterations. Numerical results indicate that the MRPI + MD method can efficiently recover the low-wavenumber components of the velocity model and is less prone to getting stuck in local minima compared to conventional RFWI.

We use surface-wave scattering in ambient-noise cross-correlations to image near-surface scatterers and major faults under the populated area of Long Beach, California under a dense geophone array. Images are computed using empirical Green’s functions from ambient-noise cross-correlation, and therefore, we eliminate the need for prior velocity models and the costly modeling of surface-waves propagation. The scattered waves are inverted in the least-squares sense to map scatterers in the subsurface. Our results show a number of faults in the area including faults in the Newport-Inglewood fault zone (NIFZ).

Many explorationists think of surface waves as the most damaging noise in land seismic data. Thus, much effort is spent in designing geophone arrays and filtering methods that attenuate these noisy events. It is now becoming apparent that surface waves can be a valuable ally in characterizing the near-surface geology. This review aims to find out how the interpreter can exploit some of the many opportunities available in surface waves recorded in land seismic data. For example, the dispersion curves associated with surface waves can be inverted to give the S-wave velocity tomogram, the common-offset gathers can reveal the presence of near-surface faults or velocity anomalies, and back-scattered surface waves can be migrated to detect the location of near-surface faults. However, the main limitation of surface waves is that they are typically sensitive to S-wave velocity variations no deeper than approximately half to one-third the dominant wavelength. For many exploration surveys, this limits the depth of investigation to be no deeper than approximately 0.5–1.0 km.

We have developed an efficient approach for picking firstbreak wavefronts on coarsely sampled time slices of 3D shot gathers. Our objective was to compute a smooth initial velocity model for multiscale full-waveform inversion (FWI). Using interactive software, first-break wavefronts were geometrically modeled on time slices with a minimal number of picks. We picked sparse time slices, performed traveltime tomography, and then compared the predicted traveltimes with the data in-between the picked slices. The picking interval was refined with iterations until the errors in traveltime predictions fell within the limits necessary to avoid cycle skipping in early arrivals FWI. This approach was applied to a 3D ocean-bottom-station data set. Our results indicate that wavefront picking has 28% fewer data slices to pick compared with picking traveltimes in shot gathers. In addition, by using sparse time samples for picking, data storage is reduced by 88%, and therefore allows for a faster visualization and quality control of the picks. Our final traveltime tomogram is sufficient as a starting model for early arrival FWI.

S U M M A R Y Field experiments are used to unequivocally demonstrate seismic superresolution imaging of subwavelength objects in the near-field region of the source. The field test is for a conventional hammer source striking a metal plate near subwavelength scatterers and the seismic data are recorded by vertical-component geophones in the far-field locations of the sources. Timereversal mirrors (TRMs) are then used to refocus the scattered energy with subwavelength resolution to the position of the original source. A spatial resolution of λ/10, where λ is the dominant wavelength associated with the data, is seen in the field tests that exceeds the Abbe resolution limit of λ/2.

We present a migration method that does not require a velocity model to migrate backscattered surface waves to their projected locations on the surface. This migration method, denoted as natural migration, uses recorded Green's functions along the surface instead of simulated Green's functions. The key assumptions are that the scattering bodies are within the depth interrogated by the surface waves, and the Green's functions are recorded with dense receiver sampling along the free surface. This natural migration takes into account all orders of multiples, mode conversions and non-linear effects of surface waves in the data. The natural imaging formulae are derived for both active source and ambient-noise data, and computer simulations show that natural migration can effectively image near-surface heterogeneities with typical ambient-noise sources and geophone distributions.

PreviousNext No Access2015 Workshop: Depth Model Building: Full-waveform Inversion, Beijing, China, 18-19 June 2015Reflection Full-waveform Inversion for Inaccurate Starting ModelsAuthors: Abdullah AlTheyab*G. T. SchusterAbdullah AlTheyab*King Abdullah University of Science and Technology (KAUST)Search for more papers by this author and G. T. SchusterKing Abdullah University of Science and Technology (KAUST)Search for more papers by this authorhttps://doi.org/10.1190/FWI2015-005 SectionsSupplemental MaterialAboutPDF/ePub ToolsAdd to favoritesDownload CitationsTrack CitationsPermissions ShareFacebookTwitterLinked InRedditEmail Abstract We propose a seismic-reflection full-waveform inversion (FWI) workflow to obtain a velocity model when the starting velocity model contains low-wavenumber errors. This workflow is based on Gauss-Seidel iterations where we apply Gauss-Newton FWI over a small subset of the seismic data. The subset is selected such that each subset invert a narrow range of the high-wavenumber components of the model space. The different subsets are inverted in sequence where the final model from inverting a subset is the starting model for the next subset. An under-relaxation operator is applied to the cumulative update from one pass of Gauss-Seidel iterations to accelerate convergence. The main advantage of the approach is that it enhances the low-wavenumber updates. Applications to synthetic and field data show significant low-wavenumber updates and flattening of common-image gathers after many iterations. Permalink: https://doi.org/10.1190/FWI2015-005FiguresReferencesRelatedDetailsCited byTarget-oriented waveform inversion based on Marchenko redatumed dataYuzhao Lin and Huaishan Liu11 January 2023 | GEOPHYSICS, Vol. 88, No. 1Multiscale Data-Driven Seismic Full-Waveform Inversion With Field Data StudyIEEE Transactions on Geoscience and Remote Sensing, Vol. 60Seismic inversion by Newtonian machine learningYuqing Chen and Gerard T. Schuster10 June 2020 | GEOPHYSICS, Vol. 85, No. 4Multiscale reflection phase inversion with migration deconvolutionYuqing Chen, Zongcai Feng, Lei Fu, Abdullah AlTheyab, Shihang Feng, and Gerard Schuster19 December 2019 | GEOPHYSICS, Vol. 85, No. 1Full unwrapped phase inversion in the phase spaceYuzhao Lin, Kai Zhang, Zhenchun Li, Renwei Din, and Zhennan Yu10 August 2019References16 June 2017 2015 Workshop: Depth Model Building: Full-waveform Inversion, Beijing, China, 18-19 June 2015ISSN (online):2159-6832Copyright: 2015 Pages: 161 publication data© 2015 Published in electronic format with permission by the Society of Exploration GeophysicistsPublisher:Society of Exploration Geophysicists HistoryPublished Online: 19 Jun 2015 CITATION INFORMATION Abdullah AlTheyab* and G. T. Schuster, (2015), "Reflection Full-waveform Inversion for Inaccurate Starting Models," SEG Global Meeting Abstracts : 18-22. https://doi.org/10.1190/FWI2015-005 Plain-Language Summary PDF DownloadLoading ...

Objectives: To evaluate continuous positive airway pressure (CPAP) compliance and define predictors of CPAP compliance among Saudi patients with obstructive sleep apnea (OSA) after applying an educational program. Methods: This prospective cohort study included consecutive patients diagnosed to have OSA based on polysomnography between January 2012 and January 2014 in King Saud University, Riyadh, Kingdom of Saudi Arabia. All patients had educational sessions on OSA and CPAP therapy before sleep study, and formal hands-on training on CPAP machines on day one, day 7, and day 14 after diagnosis. The follow-up in the clinic was carried out at one, 4, and 10 months after initiating CPAP therapy. Continuous positive airway pressure compliance was assessed objectively. Logistic regression model was used to assess the predictors of CPAP adherence. Results: The study comprised 156 patients with a mean age of 51.9±12.1 years, body mass index of 38.4±10.6 kg/m2, and apnea hypopnea index of 63.7±39.3 events/hour. All patients were using CPAP at month one, 89.7% at month 4, and 83% at month 10. The persistence of CPAP-related side effects and comorbid bronchial asthma remained as independent predictors of CPAP compliance at the end of the study. Conclusion: With intensive education, support, and close monitoring, more than 80% of Saudi patients with OSA continued to use CPAP after 10 months of initiating CPAP therapy.

Summary We present a method for inverting seismic reflections using full-waveform inversion (FWI) with inaccurate starting models. For a layered medium, near-offset reflections (with zero angle of incidence) are unlikely to be cycle-skipped regardless of the low-wavenumber velocity error in the initial models. Therefore, we use them as a starting point for FWI, and the subsurface velocity model is then updated during the FWI iterations using reflection wavepaths from varying offsets that are not cycle-skipped. To enhance low-wavenumber updates and accelerate the convergence, we take several passes through the non-linear Gauss-Seidel iterations, where we invert traces from a narrow range of near offsets and finally end at the far offsets. Every pass is followed by applying smoothing to the cumulative slowness update. The smoothing is strong at the early stages and relaxed at later iterations to allow for a gradual reconstruction of the subsurface model in a multiscale manner. Applications to synthetic and field data, starting from inaccurate models, show significant low-wavenumber updates and flattening of common-image gathers after many iterations.

PreviousNext No AccessSEG Technical Program Expanded Abstracts 2015Controlled Noise SeismologyAuthors: Sherif M. Hanafy*Abdullah AlTheyabGerard T. SchusterSherif M. Hanafy*King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi ArabiaSearch for more papers by this author, Abdullah AlTheyabKing Abdullah University of Science and Technology (KAUST), Thuwal, Saudi ArabiaSearch for more papers by this author, and Gerard T. SchusterKing Abdullah University of Science and Technology (KAUST), Thuwal, Saudi ArabiaSearch for more papers by this authorhttps://doi.org/10.1190/segam2015-5906063.1 SectionsSupplemental MaterialAboutPDF/ePub ToolsAdd to favoritesDownload CitationsTrack CitationsPermissions ShareFacebookTwitterLinked InRedditEmail Abstract We use controlled noise seismology (CNS) to generate surface waves, where we continuously record seismic data while generating artificial noise along the profile line. To generate the CNS data we drove a vehicle around the geophone line and continuously recorded the generated noise. The recorded data set is then correlated over different time windows and the correlograms are stacked together to generate the surface waves. The virtual shot gathers reveal surface waves with moveout velocities that closely approximate those from active source shot gathers. Keywords: passive, 2D, surface wavePermalink: https://doi.org/10.1190/segam2015-5906063.1FiguresReferencesRelatedDetailsCited byAutoencoded Elastic Wave-Equation Traveltime Inversion: Toward Reliable Near-Surface TomogramIEEE Transactions on Geoscience and Remote Sensing, Vol. 61Wave-equation dispersion inversion of Love wavesJing Li, Sherif Hanafy, Zhaolun Liu, and Gerard T. Schuster12 August 2019 | GEOPHYSICS, Vol. 84, No. 5Wave-equation Rayleigh-wave dispersion inversion using fundamental and higher modesZhen-Dong Zhang and Tariq Alkhalifah27 May 2019 | GEOPHYSICS, Vol. 84, No. 4Two robust imaging methodologies for challenging environments: Wave-equation dispersion inversion of surface waves and guided waves and supervirtual interferometry + tomography for far-offset refractionsJing Li, Kai Lu, Sherif Hanafy, and Gerard Schuster25 October 2018 | Interpretation, Vol. 6, No. 4Dispersion inversion of guided P-waves in a waveguide of arbitrary geometryJing Li, Sherif Hanafy, and Gerard T. Schuster27 August 2018Tutorial for wave-equation inversion of skeletonized dataKai Lu, Jing Li, Bowen Guo, Lei Fu, and Gerard Schuster1 June 2017 | Interpretation, Vol. 5, No. 3Skeletonized wave-equation Qs tomography using surface wavesJing Li, Gaurav Dutta, and Gerard T. Schuster17 August 2017Mining and Geothermal and Near Surface Complete Session17 August 2017Wave-equation dispersion inversion10 December 2016 | Geophysical Journal International, Vol. 208, No. 3Opportunities and pitfalls in surface-wave interpretationGerard T. Schuster, Jing Li, Kai Lu, Ahmed Metwally, Abdullah AlTheyab, and Sherif Hanafy20 January 2017 | Interpretation, Vol. 5, No. 1Ray-map migration of transmitted surface wavesJing Li and Gerard T. Schuster25 August 2016 | Interpretation, Vol. 4, No. 4Skeletonized wave equation of surface wave dispersion inversionJing Li and Gerard Schuster1 September 2016Seismic Inversion Complete Session1 September 2016Skeletonized inversion of surface wave: Active source versus controlled noise comparisonJing Li and Sherif Hanafy12 July 2016 | Interpretation, Vol. 4, No. 3 SEG Technical Program Expanded Abstracts 2015ISSN (print):1052-3812 ISSN (online):1949-4645Copyright: 2015 Pages: 5634 publication data© 2015 Published in electronic format with permission by the Society of Exploration GeophysicistsPublisher:Society of Exploration Geophysicists HistoryPublished Online: 19 Aug 2015 CITATION INFORMATION Sherif M. Hanafy*, Abdullah AlTheyab, and Gerard T. Schuster, (2015), "Controlled Noise Seismology," SEG Technical Program Expanded Abstracts : 5102-5106. https://doi.org/10.1190/segam2015-5906063.1 Plain-Language Summary Keywordspassive2Dsurface wavePDF DownloadLoading ...

Summary Super-virtual refraction interferometry enhances the signal-to-noise ratio of far-offset refractions. However, when applied to 3D cases, traditional 2D SVI suffers because the stationary positions of the source-receiver pairs might be any place along the recording plane, not just along a receiver line. Moreover, the effect of enhancing the SNR can be limited because of the limitations in the number of survey lines, irregular line geometries, and azimuthal range of arrivals. We have developed a 3D SVI method to overcome these problems. By integrating along the source or receiver lines, the cross-correlation or the convolution result of a trace pair with the source or receiver at the stationary position can be calculated without the requirement of knowing the stationary locations. In addition, the amplitudes of the cross-correlation and convolution results are largely strengthened by integration, which is helpful to further enhance the SNR. In this paper, both synthetic and field data examples are presented, demonstrating that the super-virtual refractions generated by our method have accurate traveltimes and much improved SNR.

To reduce human labor in picking traveltimes in 3D data, first-arrivals are picked geographically on coarsely-sampled time-slices of the recorded shot gathers. Traveltime picks are interpolated in real-time in the polar coordinates to minimize the number of picks needed to track curved wavefronts. The picked firstarrival traveltimes are used in ray tomography to form an initial model for FWI. Results indicate an 80% reduction in human picking time compared to standard picking in the shot gather or common offset domains.

PreviousNext No AccessSEG Technical Program Expanded Abstracts 2013Hybrid linear and non-linear full-waveform inversion of Gulf of Mexico dataAuthors: Abdullah AltheyabXin WangAbdullah AltheyabKing Abdullah U of Science and TechnologySearch for more papers by this author and Xin WangKing Abdullah U of Science and TechnologySearch for more papers by this authorhttps://doi.org/10.1190/segam2013-0538.1 SectionsSupplemental MaterialAboutPDF/ePub ToolsAdd to favoritesDownload CitationsTrack CitationsPermissions ShareFacebookTwitterLinked InRedditEmail Abstract Standard Full-waveform inversion (FWI) often suffers from poor sensitivity to deep features of the subsurface model. To alleviate this problem, we propose a hybrid linear and non-linear optimization method to enhance the FWI results. In this method, iterative least-squares reverse-time migration (LSRTM) is used to estimate the model update at each nonlinear iteration, and the number of LSRTM iterations is progressively increased after each non-linear iteration. With this method, model updating along deep reflection wavepaths are automatically enhanced, which in turn improves imaging below the reach of diving-waves. This hybrid linear and non-linear FWI algorithm is implemented in the space-time domain to simultaneously invert the data over a range of frequencies. A multiscale approach is used where higher frequencies are iteratively incorporated into the inversion. Synthetic data are used to test the effectiveness of reconstructing both the high- and low-wavenumber features in the model without relying on diving waves in the inversion. We apply the method to Gulf of Mexico field data and illustrate the improvements after several iterations. Results show a significantly improved migration image in both the shallow and deep sections. Permalink: https://doi.org/10.1190/segam2013-0538.1FiguresReferencesRelatedDetailsCited byReverse time adjoint migrationHongwei Liu, Almomin Ali, and Yi Luo29 March 2019 | GEOPHYSICS, Vol. 84, No. 3References16 June 2017Far-field superresolution by imaging of resonance scatteringGerard T. Schuster* and Yunsong Huang5 August 2014 SEG Technical Program Expanded Abstracts 2013ISSN (print):1052-3812 ISSN (online):1949-4645Copyright: 2013 Pages: 5258 Publisher:Society of Exploration Geophysicists HistoryPublished Online: 19 Aug 2013 CITATION INFORMATION Abdullah Altheyab and Xin Wang, (2013), "Hybrid linear and non-linear full-waveform inversion of Gulf of Mexico data," SEG Technical Program Expanded Abstracts : 1003-1007. https://doi.org/10.1190/segam2013-0538.1 Plain-Language Summary PDF DownloadLoading ...

PreviousNext No AccessSEG Technical Program Expanded Abstracts 2013Time-domain incomplete Gauss-Newton full-waveform inversion of Gulf of Mexico dataAuthors: Abdullah AlTheyabXin WangGerard T. SchusterAbdullah AlTheyabKing Abdullah U of Science and TechnologySearch for more papers by this author, Xin WangKing Abdullah U of Science and TechnologySearch for more papers by this author, and Gerard T. SchusterKing Abdullah U of Science and TechnologySearch for more papers by this authorhttps://doi.org/10.1190/segam2013-1478.1 SectionsAboutPDF/ePub ToolsAdd to favoritesDownload CitationsTrack CitationsPermissions ShareFacebookTwitterLinked InRedditEmail Abstract We apply the incomplete Gauss-Newton full-waveform inversion (TDIGN-FWI) to Gulf of Mexico (GOM) data in the space-time domain. In our application, iterative least-squares reverse-time migration (LSRTM) is used to estimate the model update at each non-linear iteration, and the number of LSRTM iterations is progressively increased after each non-linear iteration. With this method, model updating along deep reflection wavepaths are automatically enhanced, which in turn improves imaging below the reach of diving-waves. The forward and adjoint operators are implemented in the space-time domain to simultaneously invert the data over a range of frequencies. A multiscale approach is used where higher frequencies are down-weighted significantly at early iterations, and gradually included in the inversion. Synthetic data results demonstrate the effectiveness of reconstructing both the high- and low-wavenumber features in the model without relying on diving waves in the inversion. Results with Gulf of Mexico field data show a significantly improved migration image in both the shallow and deep sections. Permalink: https://doi.org/10.1190/segam2013-1478.1FiguresReferencesRelatedDetailsCited ByAccelerating Hessian-free Gauss-Newton full-waveform inversion via l-BFGS preconditioned conjugate-gradient algorithmWenyong Pan, Kristopher A. Innanen, and Wenyuan Liao24 January 2017 | GEOPHYSICS, Vol. 82, No. 2References16 June 2017Wavefront picking for 3D tomography and full-waveform inversionGEOPHYSICS, Vol. 81, No. 6Adaptive overburden elimination with the multidimensional Marchenko equationGEOPHYSICS, Vol. 81, No. 5Full waveform inversion of diving & reflected waves for velocity model building with impedance inversion based on scale separation8 July 2015 | Geophysical Journal International, Vol. 202, No. 3 SEG Technical Program Expanded Abstracts 2013ISSN (print):1052-3812 ISSN (online):1949-4645Copyright: 2013 Pages: 5258 Publisher:Society of Exploration Geophysicists HistoryPublished: 19 Aug 2013 CITATION INFORMATION Abdullah AlTheyab, Xin Wang, and Gerard T. Schuster, (2013), "Time-domain incomplete Gauss-Newton full-waveform inversion of Gulf of Mexico data," SEG Technical Program Expanded Abstracts : 5175-5179. https://doi.org/10.1190/segam2013-1478.1 Plain-Language Summary PDF DownloadLoading ...

We present the first results for a controlled source seismic experiment where we demonstrate that subwavelength scatterers can be imaged with a resolution of λ/8 using seismic data recorded in the far field of the scatterer. Using the Time Reverse Mirror (TRM) operation, we show that an image resolution Δx of 0.6 m can be obtained at the source location from data with seismic wavelengths ≤ 5 m. In other words, we show that it is possible to extract 220 Hz information from 55 Hz data using the TRM operation. These results also validate the theory of seismic superresolution associated with a seismic scanning tunneling macroscope.

In an effort to understand the transmission effects of localized heterogenities in the subsurface, we present the travel-time and amplitude distortions caused by localized variations in velocity and absorption. To examine the relative impact of velocity and absorption heterogeneities on seismic events, we conducted numerical experiments using visco-acoustic finite-difference modeling of the linearized waveequation for Newtonian fluids. We analyzed the distortions in the midpoint-offset domain. We find that the distortion caused by an anomaly that is both slow and absorptive is different from that an anomaly that is either slow or absorptive, but not both. Our results also indicate that amplitude distortion of highly absorptive anomalies (Q < 50) can be comparable to that of small velocity variation (less than 4%), and therefore absorption must be considered in seismic amplitude inversion and AVO analysis. INTRODUCTION Localized heterogeneities in the subsurface cause amplitude and travel-time distortions of seismic reflections from underlying reflectors. These distortions are problematic to imaging and AVO analysis. The distortions come in almost regular patterns and usually are recognizable by V-shaped trajectories in the midpoint-offset domain (X-shapes for split-spread acquisition geometry). These distortions can be used to find the locations of the anomalies, which can reveal valuable information for interpreters such as fault locations (Hatchell, 2000). Moreover, they can be used to invert for velocity and absorption anomalies. The analyses of several authors Vlad (2005), Hatchell (2000) and Harlan (1994) have considered mostly velocity anomalies, which cause focusing and defocusing effects. In this report, we stress that absorption must be considered in the analysis of these distortions. A seismic amplitude inversion that disregards absorption is likly to be biased, especially if velocity perturbations of interest are less than 4%. We examine the relative impact of localized velocity anomalies versus absorption anomalies on seismic amplitude.

Nvidia’s graphics processing units (GPU) powered with Compute Unified Device Architecture (CUDA), the supporting API, have allowed a significant speedup to finite difference time domain (FDTD) seismic modeling and, consequently, to reverse time migration (RTM). To utilize the power of GPUs for velocity analysis, we implemented kernels for generating offset-domain common image gathers (ODCIGs). With 4GB of memory, a single Tesla 10 series GPU can perform the 2D RTM with generation of ODCIGs. Computing the ODCIGs takes the majority of the algorithm execution time because of the large volume of output. We examine the performance of a general 2D RTM algorithm with ODCIGs computation fully offloaded to a single GPU device. We optimized the imaging kernel utilizing the available shared memory on the GPU to double the throughput of the kernel. INTRODUCTION Reverse time migration (RTM) is a full wave equation imaging technique that constructs an image that best represents the subsurface structure. Seismic data are migrated using an estimate of the wave propagation velocity in the subsurface. Estimates of the velocity field can be inaccurate at the first imaging attempt, and subsurface offset gathers can provide a measure of the errors in velocity estimation (Biondi and Symes, 2004). In addition, they can give amplitude versus offset (AVO) information if amplitudes are handled properly. Algorithms for RTM and generation of ODCIGs are known to be computationally expensive and sometimes unaffordable. Therefore, an efficient implementation of RTM with ODCIGs generation algorithm is vital for minimizing the time required for creating a complete image. Reverse time migration falls into the computational class of convolution with a stencil. The stencil computation workload can be divided among many processing units in an embarrassingly parallel fashion. However, the main performance limitation on modern computer architectures is the memory latency. Cache-aware algorithms minimize data traffic by taking advantage of spatial and temporal locality and/or the data prefetch capabilities of modern CPUs. Another way of hiding memory latency is to have more threads than cores to execute some threads while other threads are waiting for memory access. The performance gain given by this technique is not

Light propagating in a water-filled pool is perturbed by the water surface, creating patterns on the pool floor. In this report, I use ray tracing to compute an approximation of the light intensity field on the pool floor using point source and exploding surface models. The ultimate goal is to infer the water surface from the intensity field. In a seismic imaging, it is similar to imaging using amplitudes rather than travel-times. I present a geometric approach to invert for a discretized surface from a ray-count representation of the intensity field. With this formulation, the inversion becomes a combinatorial problem, which can be solved using non-deterministic search techniques. The formulation has a large inherent null space. The low cost of the technique allows a large number of iterations to be applied. INTRODUCTION Light propagation is similar in both nature and theory to seismic wave propagation in the subsurface. It is not surprising to see a large number of common problems between the fields of seismology and optics. Claerbout (2007) poses a question about the relationship between light patterns under water and seismology. One prominent problem in reflection seismology is estimating seismic velocities in the heterogeneous subsurface. Velocity is a measure that depends on travel time; i.e. it is determined by the moveout of seismic events. It is always costly to estimate subsurface velocities, whether using conventional velocity analysis or inversion methodologies. Unfortunately, imaging algorithms rely on the accuracy of velocity models. Therefore, a curious scientist might ask: Can we use amplitudes alone for imaging? For the sake of simplicity, let us consider a simple imaging problem involving a single wavefield caustic: a single reflector including a syncline. Using the exploding reflector concept, we can simulate a zero-offset section (Claerbout, 1985). What we will observe in the section is a bowtie that is caused by the syncline. The syncline causes many ray paths between the exploding reflector to the receiver to converge at some point before reaching the surface, forming a caustic. Given the correct velocity model, migration algorithms like Kirchhoff migration will resolve the bowtie into a syncline (Yilmaz, 2001). Al Theyab 2 Light under water Now, let us consider a swimming-pool experiment where the goal is to infer the water surface from the light patterns on the floor of the pool like the ones shown in Figure 1. We might encounter similar difficulties as in seismic imaging. It is hard to estimate the exact speed of light in the pool. In this case the question is as follows: Can we use the light intensity field at the bottom of the pool to infer the surface of water? The two questions are closely related, and the pool experiment is easier to comprehend intuitively. In this paper, I implement the forward simulation using ray theory to build a light intensity field at the pool floor that can be used for inferring the surface. I will not follow the exact physics as it is done in optical simulations: I assume infinite frequencies for ray theory and representing light with a finite number of rays. The forward simulation is done using ray tracing and beam tracing. I finally discuss a Monte Carlo ray tracing inversion that is based on ray counting. Figure 1: Light patterns at the bottom of a pool (Claerbout, 2007). [NR] FORWARD SIMULATION I simulate the light intensity field at the pool floor using ray tracing. For our purposes, we assume that light travels with an infinite velocity; i.e. any ray refracting through the water surface will project instantly onto another point on the floor. If multiple rays for any reason converge at a point before reaching the floor, they form a refracted caustic (Shah et al., 2007). This simplistic view of the experiment has received considerable attention in the field of computer graphics. Shah et al. (2007), for example, introduce a real-time technique for rendering caustics from reflective and refractive surfaces. Ray tracing techniques have a major deficiency, in that we need an infinite number of rays to simulate a realistic distribution of light. This shortcoming is discussed by Watt (1990) who suggests using beam tracing as an alternative. The vast majority of the techniques developed in computer graphics, however, are developed for a 3D space having point light source(s), polygonal objects, and a single observation point. Al Theyab 3 Light under water Our simplified case inherits only a minor subset of the wide range of techniques used in computer graphics. To start, we will not have an observer point, and therefore do not need to ray trace from the pool floor to the observer point. Also, we have only a single smooth surface of moderate relief, and the incident light direction is limited to rays arriving from above the surface. In other words, we have an unique refraction from each point on the surface. Because we have moderate topology in the water surface, we will assume that the refracted rays do not re-intersect the water surface. Surface Representation In (x, y, z)-space, the surface of the water is given as a discretized function on a uniformly sampled grid in one or two dimensions. The function has points that represent the depth of the surface from the xy-plane. If we have a surface S = f(x, y) in 3D space, the normal vector to that surface is ~n = − ∂x î− ∂S ∂y ĵ + k̂ √ ( ∂x )2 + ( ∂y )2 + 1 , (1) and for a surface S = f(x), in 2D space, it is ~n = − ∂x î + k̂ √ ( ∂x )2 + 1 . (2) where î, ĵ, and k̂ are unit vectors along the x-, y-, z-axes, respectively. To compute the first partial derivatives, we can use finite-divided-difference approximations (Chapra and Canale, 2002). A design decision must be made at this point regarding where the refracted rays will start with respect to the points on the surface function. For a 2D simulation with S(x), we can have the ray coming out of the points. In that case, the preferred scheme is centered finite-divided-difference: ∂ ∂x S(xi) ≈ S(xi+1)− S(xi−1) 2∆x . (3) The same differencing scheme can be used for the partial derivatives with respect to x and y for the 3D simulation. Using this scheme, it is not possible to compute the partial derivatives for the edges of the surface, since more points are needed for the computation. The second design option is to have the rays launching from segments of the surface between the points in the 2D simulation. For the 3D case, we divide every quad – i.e. square between four adjacent points– using two diagonals. For every possible triangle within a quad, we compute the normal and have a ray starting from the center of the triangle. The advantage of this approach is that we have many more rays than the number of surface samples. For computational stability, the surface must be smooth locally with respect to a segment or a quad. Al Theyab 4 Light under water An easy way to construct a synthetic water surface with a sinusoidal wave with a decaying factor as a function of the distance from the source: S(x, y) = c0 cos ( c2|~ d| ) e−c1| ~ d| , (4) where ~ d is the distance from the source, and c0, c1, and c2 are arbitrary positive constants to customize amplitudes, decay rate, and angular frequency of the ripples, respectively. Waves originating from several source points can be summed to form a more complex surface. This method of construction is for a time-invariant surface without reflection at the surface boundaries. A more realistic water surface can be obtained using finite differencing of the 2D wave equation. For the sake of simplicity, I use the explicit differencing scheme in the (x, y, t)-domain. Figure 2 is the result of finite differencing with three sources that are Gaussian wavepackets of different amplitudes. Figure 2: Finite-difference simulation of a water surface with one strong and two weaker sources disturbing the water surface. [ER]