Chu, C., P. L. Stoffa, and R. K. Seif, High-order rotated staggered finite difference modeling of 3D elastic wave propagation in general anisotropic media, SEG Expanded Abstr., 28, 291, 2009, doi:10.1190/1.3255458, #2359 
We analyze the dispersion properties and stability conditions of the high-order convolutional finite difference operators and compare them with the conventional finite difference schemes. We observe that the convolutional finite difference method has better dispersion properties and becomes more efficient than the conventional finite difference method with the increasing order of accuracy. This makes the high-order convolutional operator a good choice for anisotropic elastic wave simulations on rotated staggered grids since its enhanced dispersion properties can help to suppress the numerical dispersion error that is inherent in the rotated staggered grid structure and its efficiency can help us tackle 3D problems cost-effectively.
Chu, C., P. L. Stoffa, and R. K. Seif, 3D Seismic modeling and reverse-time migration with the parallel Fourier method using non-blocking collective communications, SEG Expanded Abstr., 28, 2677, 2009, doi:10.1190/1.3255403, #2360
The major performance bottleneck of the parallel Fourier method on distributed memory systems is the network communication cost. In this study, we investigate the potential of using non-blocking all-to-all communications to solve this problem by overlapping computation and communication. We present the runtime comparison of a 3D seismic modeling problem with the Fourier method using non-blocking and blocking calls, respectively, on a Linux cluster. The data demonstrate that a performance improvement of up to 40% can be achieved by simply changing blocking all-to-all communication calls to non-blocking ones to introduce the overlapping capability. A 3D reverse-time migration result is also presented as an extension to the modeling work based on non-blocking collective communications.
Chu, C., P. L. Stoffa, and R. K. Seif, 3D Elastic wave modeling using modified high-order time stepping schemes with improved stability conditions, SEG Expanded Abstr., 28, 2662, 2009, doi:10.1190/1.3255400, #2361
We present two Lax-Wendroff type high-order time stepping schemes and apply them to solving the 3D elastic wave equation. The proposed schemes have the same format as the Taylor series expansion based schemes, only with modified temporal extrapolation coefficients. We demonstrate by both theoretical analysis and numerical examples that the modified schemes significantly improve the stability conditions.
Hu, C. S., and P. L. Stoffa, Slowness-driven Gaussian-beam prestack depth migration for low-fold seismic data, Geophysics, 74, WCA35-WCA45, 2009, doi:10.1190/1.3250268, #2114
Subsurface images based on low-fold seismic reflection data or data with geometry acquisition limitations, such as obtained from ocean-bottom seismography (OBS), are often corrupted by migration swing artifacts. Incorporating prestack instantaneous slowness information into the imaging condition can significantly reduce these migration swing artifacts and improve image quality, especially for areas with poor illumination. We combine the horizontal surface slowness information of observed seismic data with Gaussian-beam depth migration to implement a new slowness-driven Gaussian-beam prestack depth migration whereby Fresnel weighting is combined naturally with beam summation. The prestack instantaneous slowness information is extracted from the original OBS or shot gathers using local slant stacks and is combined with a local semblance analysis. During migration, we propagate the seismic energy downward, knowing its instantaneous slowness information. At each image location, the beam summation is localized in a resolution-dependent Fresnel zone; the instantaneous slowness information controls the beam summation. Synthetic and real data examples confirm that slowness-driven Gaussian-beam migration can suppress most noise from inadequate stacking and give a clearer migration result.
Jin, L., M. K. Sen, and P. L. Stoffa, Optimal model parameterization in stochastic inversion for reservoir properties using time-lapse seismic and production data, SEG Ann. Meeting, San Antonio, TX, 2009, #1954
Jin, L., M. K. Sen, and P. L. Stoffa, Fusion based classification method and its application, J. Seismic Explor., 18, 103-117, 2009, #2091
Jin, L., P. L. Stoffa, M. K. Sen, R. K. Seif, and A. R. Sena, Pilot point parameterization in stochastic inversion for reservoir properties using time-lapse seismic and production data, J. Seismic Explor., 18, 1-20, 2009, 1 citation, #2097
Jin, L., M. K. Sen, and P. L. Stoffa, One-dimensional prestack seismic waveform inversion using ensemble Kalman filter, SEG Expanded Abstr., 28, 1810, 2009, doi:10.1190/1.3255204, #2371
Estimation of reservoir parameters is important for prediction and optimization of oil production. Time-lapse seismic data, commonly acquired in mature oil fields to identify profitable bypassed zones, can be used to estimate reservoir parameters by posing the history matching process as a joint inversion problem. Due to the size of the problem and existence of multiple minima (of the error function), only global optimization techniques, such as Very Fast Simulated Annealing (VFSA) have an opportunity to find the global minimum. However, even for the VFSA technique, the computational cost of the forward modeling (particularly, the reservoir simulation) can make the inversion process impractical in real applications if we consider that many inversions are required to estimate the solution uncertainty. To improve the convergence characteristics of the VFSA technique, we have developed a modification of this technique where a regulation of the local temperature associated with the model variables is driven by the local changes of the error function. We call this modification Local Thermal Regulation (LTR) and show that the time (or number of iterations) spent for convergence can be reduced by 25 to 50%, while maintaining the same degree of ergodicity of the conventional VFSA technique.
Pestana, R. C., and P. L. Stoffa, Rapid expansion method (REM) for time-stepping in reverse time migration (RTM), SEG Expanded Abstr., 28, 2819, 2009, doi:10.1190/1.3255434, #2368
We show that the wave equation solution using a conventional finite-difference scheme, derived commonly by the Taylor series approach, can be derived directly from the rapid expansion method (REM). After some mathematical manipulation we consider an analytical approximation for the Bessel function where we assume that the time step is sufficiently small. From this derivation we find that if we consider only the first two Chebyshev polynomials terms in the rapid expansion method we can obtain the second order time finite-difference scheme that is frequently used in more conventional finite-difference implementations. We then show that if we use more terms from the REM we can obtain a more accurate time integration of the wave field. Consequently, we have demonstrated that the REM is more accurate than the usual finite-difference schemes and it provides a wave equation solution which allows us to march in large time steps without numerical dispersion and is numerically stable. We illustrate the method with post and pre stack migration results.
Sena, A. R., P. L. Stoffa, M. K. Sen, and R. K. Seif, Assessing the value of time-lapse seismic data in joint inversion for reservoir parameter estimation in an oil reservoir subjected to water flooding recovery: A synthetic example, SEG Expanded Abstr., 28, 3904, 2009, doi:10.1190/1.3255683, #2366
Joint inversion of time-lapse seismic and production data for reservoir parameter estimation is an extension of automatic history matching where seismic data introduce spatial constraints. In joint inversion, the objective function (global error) is a weighted sum of both the seismic and production errors, and the seismic term includes as many time-lapses as available. An inherent question in this problem is âwill the availability of more seismic time-lapses help to appreciably improve the inversion results?â In this work we attempt to answer this question in a stringent case: a clastic oil reservoir subjected to water flooding recovery. This case is challenging because of the very small changes of elastic impedance due to substitution of oil by water. In this study, the reservoir under consideration is a sub-volume of the SPE 10 model, and we use Very Fast Simulated Annealing with Local Thermal Regulation as the optimization technique. The results suggest that the interval between time-lapse seismic surveys could be planned around the time when the water break is expected to occur in the production well. However, time-lapse seismic surveys between these two times could be important to improve the inversion results if rock-physics and initial fluid composition in the reservoir are not very well constrained.
Sena, A. R., M. K. Sen, P. L. Stoffa, R. K. Seif, and L. Jin, Joint inversion of time-lapse seismic and production data using VFSA with local thermal regulation and pilot point parameterization, SEG Expanded Abstr., 29, 1810, 2009, doi:10.1190/1.3255204, #2367
Estimation of reservoir parameters is important for prediction and optimization of oil production. Time-lapse seismic data, commonly acquired in mature oil fields to identify profitable bypassed zones, can be used to estimate reservoir parameters by posing the history matching process as a joint inversion problem. Due to the size of the problem and existence of multiple minima (of the error function), only global optimization techniques, such as Very Fast Simulated Annealing (VFSA) have an opportunity to find the global minimum. However, even for the VFSA technique, the computational cost of the forward modeling (particularly, the reservoir simulation) can make the inversion process impractical in real applications if we consider that many inversions are required to estimate the solution uncertainty. To improve the convergence characteristics of the VFSA technique, we have developed a modification of this technique where a regulation of the local temperature associated with the model variables is driven by the local changes of the error function. We call this modification Local Thermal Regulation (LTR) and show that the time (or number of iterations) spent for convergence can be reduced by 25 to 50%, while maintaining the same degree of ergodicity of the conventional VFSA technique.
Shahin, A., P. L. Stoffa, R. H. Tatham, and D. Sava, Multicomponent seismic time-lapse cross-plot and its applications, SEG Expanded Abstr., 28, 1227, 2009, doi:10.1190/1.3255073, #2364
A synthetic poorly consolidated sandstone reservoir of varying thickness, porosity, and clay content, embedded in shale is considered for this study. A reservoir segment which has constant properties equal to average reservoir properties is used for calculating seismic attributes including travel times and reflection coefficients (RCs). For a range of water saturation and pore pressure, Multi-component (MC) seismic, i.e., P-P, P-SV, and SH-SH, travel times trough and RCs at top of the reservoir segment have been calculated. Then, we computed changes in travel times and RCs with respect to a reference travel time and RC calculated at a reference saturation and pressure. We plotted changes in RCs versus changes in travel times. The corresponding time-lapse cross-plot shows interesting patterns for saturation and pressure changes and has a potential for quantitatively discriminating pressure and saturation changes. Next, we deployed a statistical method to determine the efficacy of MC seismic in detecting production-induced time-lapse changes. The significant and representative data in time-lapse cross-plot allows us to statistically analyze the detectability of a known scenario of saturation and pressure changes using MC seismic attributes. Applying different thresholds for travel times and RCs, we constructed joint probability detectors (JPDs) that help us to compare the likelihood of detection of a known change in dynamic reservoir properties using different component of seismic data. Repeating the above mentioned procedures for reservoir segments with various porosities, thicknesses, and clay contents gives insight into a classification of various reservoir types or different parts of a laterally varying reservoir in terms of saturation and pressure discrimination and statistical detectability.
Stoffa, P. L., and R. C. Pestana, Numerical solution of the acoustic wave equation by the rapid expansion method (REM)- A one step time evolution algorithm, SEG Expanded Abstr., 28, 2672, 2009, doi:10.1190/1.3255402, #2369
We revisit a numerical solution for the time evolution of the acoustic wave equation where the wave field can be advanced in time in one large step, the rapid expansion method, or REM. This method has the advantage of being numerically stable even for large time steps and when combined with a pseudo spectral Fourier method for the spatial derivatives results in a highly accurate solution for the propagating wave fields. We show that the number of terms in the integration that produces the time snap shots can be significantly reduced making the method computationally attractive for high resolution simulations.
Chu, C., and P. L. Stoffa, A pseudospectral-finite difference hybrid approach for large-scale seismic modeling and RTM on parallel computers, SEG Expanded Abstr., 27, 2087, 2008, doi:10.1190/1.3059302, #2362
We present numerical comparisons between the finite difference method and the Fourier pseudospectral method on the 3D SEG/EAGE salt model to investigate the potentials of the two methods for large scale seismic modeling and reverse-time migration problems on distributed-memory parallel systems. We compare the accuracy of the seismic modeling results by the two methods. We discuss the parallel implementations and present the runtime and speedup rate data on a Linux cluster. Based on these comparison results, we propose a pseudospectral-finite difference hybrid method for largescale seismic modeling and RTM problems on massively distributed-memory parallel systems.
Hu, C. S., P. L. Stoffa, and K. D. McIntosh, First arrival stochastic tomography: Automatic background velocity estimation using beam semblances and VFSA, Geophys. Res. Lett., 35, L23307, 2008, 1 citation, doi:10.1029/2008GL034776, #2113
We present a new tomography method based on the local beam semblance and the very fast simulated annealing (VFSA) global optimization method. The data space is the local beam semblance calculated using local slant stacks for overlapping offset windows, i.e. beam windows, of the original common-shot or common-receiver gathers. On each beam semblance panel, the first coherency peak can be identified with a particular ray parameter, first-arrival traveltime and beam center position. The forward problem can be solved with any ray tracer to find arrivals matching the identified peaks. Our inversion scheme uses VFSA to find the maximum-a-posteriori (MAP) solution and estimates the uncertainty by applying Bayesian analysis of all the sampled models for a specified model parameterization. This integration of automatic local semblance evaluation instead of first-arrival picking and a fast forward modeling method combined with VFSA to determine the optimal model makes our method robust, efficient and accurate.
Hu, C. S., P. L. Stoffa, and K. D. McIntosh, First arrival stochastic tomography: Automatic background velocity estimation using beam semblances and VFSA, SEG Expanded Abstr., 27, 3305, 2008, doi:10.1190/1.3064031, #2370
We present a new stochastic tomography method for automatic background velocity estimation based on the local beam semblance of common-shot or common-receiver gathers and the very fast simulated annealing (VFSA) global optimization method. The data space is the local beam semblance which is calculated using the local slant stacks for overlapping offset windows, i.e. beam windows, of the original common-shot or common-receiver gathers. On each beam semblance panel, the first coherency peak represents the first-arrival energy which can now be located with a particular ray parameter and traveltime and associated with the offset at the center of each beam window. The forward problem can be solved with an eikonal solver with/without ray tracing or beam method to find these first coherency peaks. Our inversion scheme uses VFSA to find the maximum a posteriori (MAP) solution and estimate the uncertainty by the application of Bayesian analysis to all the samples which are based on an optimal model parameterization. This integration of automatic local semblance evaluations instead of first-arrival picking, optimal model parameterization estimation and fast ray tracing or beam forward modeling method through VFSA makes our stochastic tomography robust, efficient and accurate.
Jin, L., M. K. Sen, P. L. Stoffa, and R. K. Seif, Time-lapse seismic attribute analysis for a water-flooded reservoir, J. Geophysics Engin., 5, 210-220, 2008, 2 citations, doi:10.1088/1742-2132/5/2/008, #1988
One of the goals of time-lapse seismic monitoring is the direct detection of the fluid front and two-phase contact area. However, several factors affect the quality of time-lapse seismic difference data and decrease detectability. One of these factors is random noise. In this paper, we propose five different methods aimed at improving the quality and detectability of noisy time-lapse seismic difference data. Common to these methods is the transform of the differences to a domain where the time-lapse signal and random noise are well separated. Our proposed methods include direct Fourier transform based spectral decomposition, bispectra, wavelet transform, singular value decomposition and hybrid methods. We also propose a method that combines multiple time-lapse difference data and gives a final difference which enhances the common part and attenuates the differences of the multiple difference images resulting in a better detectability than the original images. A synthetic time-lapse model is used to demonstrate the feasibility of our proposed methods.
Jin, L., P. L. Stoffa, and M. K. Sen, Stochastic inversion for reservoir properties using parallel learning based VFSA and pilot point parameterization, Soc. Petroleum Engin., SPE-118818-PP, 2008, #2102
Sena, A. R., M. K. Sen, and P. L. Stoffa, Modeling of ground penetrating radar data in stratified media using the reflectivity technique, J. Geophysics Engin., 5, 129-146, 2008, doi:10.1088/1742-2132/5/2/001, #2110
Horizontally layered media are often encountered in shallow exploration geophysics. Ground penetrating radar (GPR) data in these environments can be modelled by techniques that are more efficient than finite difference (FD) or finite element (FE) schemes because the lateral homogeneity of the media allows us to reduce the dependence on the horizontal spatial variables through Fourier transforms on these coordinates. We adapt and implement the invariant embedding or reflectivity technique used to model elastic waves in layered media to model GPR data. The results obtained with the reflectivity and FDTD modelling techniques are in excellent agreement and the effects of the airâsoil interface on the radiation pattern are correctly taken into account by the reflectivity technique. Comparison with real wide-angle GPR data shows that the reflectivity technique can satisfactorily reproduce the real GPR data. These results and the computationally efficient characteristics of the reflectivity technique (compared to FD or FE) demonstrate its usefulness in interpretation and possible model-based inversion schemes of GPR data in stratified media.
Shahin, A., P. L. Stoffa, R. H. Tatham, and D. Sava, Sensitivity analysis of multicomponent seismic attributes to fluid content and pore pressure, SEG Expanded Abstr., 27, 247, 2008, doi:10.1190/1.3054798, #2365
We analyzed the sensitivity of multicomponent (MC) seismic (P-P, P-SV, and SH-SH) reflection coefficients (RCs) and traveltimes to water saturation and pore pressure by taking the appropriate partial derivatives. Applying this approach to a poorly consolidated sandstone reservoir partially saturated with light oil and brine, demonstrates that P-P traveltimes have the largest sensitivity to water saturation, but the least sensitivity to pore pressure. In contrast, SH-SH traveltimes have the least sensitivity to water saturation, and the most sensitivity to pressure. P-SV traveltimes have intermediate sensitivities to pressure and saturation, but are more affected by saturation than pressure. By analyzing the sensitivity of MC seismic RCs to water saturation at the reservoir top, i.e., shale over oil-saturated sandstone, and at the oil-water contact (OWC), i.e., oil-saturated sandstone over brine-saturated sandstone, we found that the absolute value of amplitudes at all angles is greatest for P-P, smallest for P-SV, and intermediate for SH-SH. In addition, the absolute values of AVO (amplitude variation vs. offset) gradients at the reservoir top and OWC can be organized in descending order as (P-SV, SH-SH, P-P) and (P-P, P-SV, SH-SH), respectively. For the sensitivity of both interfaces to pressure (the reservoir top and OWC), angle-dependent relations are extracted.
Wood, W. T., W. S. Holbrook, M. K. Sen, and P. L. Stoffa, Full waveform inversion of reflection seismic data for ocean temperature profiles, Geophys. Res. Lett., 35, L040608, 2008, 11 citations, doi:10.1029/2007GL032359, #2111
We show that ocean temperature profiles can be accurately recovered using only acoustic methods employed at the sea surface. Using a towed air gun array and a hydrophone streamer, thermohaline boundaries are ensonified at a suite of frequencies and angles, yielding travel time trajectories and reflectivities. These data are inverted via full waveform inversion to estimate sound speed and, subsequently, a temperature profile. The high lateral data density of the seismic technique offers the potential of acoustically derived temperature profiles to be used to constrain models of ocean mixing and internal waves. Results on realistic synthetic data show that sound speed can be recovered with arbitrary accuracy when using broadband data, with known source function and recording geometry. Application to field seismic data (corroborated by expendable bathythermograph) shows that even with a seismic acquisition system not specifically calibrated for seismic oceanography, temperature contrasts within the ocean can be recovered to within one degree Celsius.
Hong, T. C., M. K. Sen, P. L. Stoffa, H. Klie, S. G. Thomas, A. Rodriguez, and M. F. Wheeler, Integrated time-lapse seismic inversion for reservoir petrophysics and fluid flow imaging, SEG Ann. Meeting, San Antonio, TX, 2007, #1955
Hu, C. S., M. K. Sen, P. L. Stoffa, and K. D. McIntosh, Plane wave gaussian beam prestack depth migration, SEG Ann. Meeting, San Antonio, TX, 2200-2204, 2007, #1960
Jin, L., M. K. Sen, and P. L. Stoffa, Fusion based classification method and its application, SEG Ann. Meeting, San Antonio, TX, 2007, #1952
Jin, L., M. K. Sen, T. C. Hong, and P. L. Stoffa, Joint estimation of porosity and saturation by combining a rock physics model and constrained pre-stack seismic waveform inversion, SEG Ann. Meeting, San Antonio, TX, 1957-1961, 2007, #1961
Hu, C. S., K. D. McIntosh, and P. L. Stoffa, Voronoi cell based staggered grid SH wave numerical simulation, SEG Ann. Meeting, New Orleans, LA, 2270-2274, 2006, #1956
Hu, C. S., K. D. McIntosh, H. J. A. Van Avendonk, and P. L. Stoffa, Hybrid ray tracer and amplitude calculation with finite difference, graph theory and ray bending
, SEG Ann. Meeting, New Orleans, LA, 3408-3412, 2006, #1957
Klie, H., W. Bangerth, X. Gai, M. F. Wheeler, P. L. Stoffa, M. K. Sen, M. Parashar, U. Catalyurek, J. Saltz, and T. Kurc, Models, methods, and middleware for grid-enabled multiphysics oil resovoir management, Engineering with Computers, Engineering with Computers, 22, 349-370, 2006, 5 citations, doi:10.1007/s00366-006-0035-9, #1840
Meeting the demands for energy entails a better understanding and characterization of the fundamental processes of reservoirs and of how human made objects affect these systems. The need to perform extensive reservoir studies for either uncertainty assessment or optimal exploitation plans brings up demands of computing power and data management in a more pervasive way. This work focuses on high performance numerical methods, tools and grid-enabled middleware systems for scalable and data-driven computations for multiphysics simulation and decision-making processes in integrated multiphase flow applications. The proposed suite of tools and systems consists of (1) a scalable reservoir simulator, (2) novel stochastic optimization algorithms, (3) decentralized autonomic grid middleware tools, and (4) middleware systems for large-scale data storage, querying, and retrieval. The aforementioned components offer enormous potential for performing data-driven studies and efficient execution of complex, large-scale reservoir models in a collaborative environment.
Sena, A. R., P. L. Stoffa, and M. K. Sen, Split-step Fourier migration of GPR data in lossy media, Geophysics, 71, K77-K91, 2006, doi:10.1190/1.2217157, #1797
Migration of ground-penetrating radar (GPR) data traditionally has been implemented assuming uniform, nonconductive, and nondispersive media. However, in many real applications, the effects of heterogeneities, conductivity, or dispersion can be important, so it is necessary to consider these effects to image the data correctly. We implement the split-step Fourier technique for migration of GPR data in 2D media transverse-electric (TE) or transverse-magnetic (TM) propagation modes and demonstrate how it takes into account, naturally and efficiently, the effects of dispersion, attenuation, and heterogeneities. Compensation for attenuation during migration implies applying a gain that can make the numerical algorithm unstable. We introduce a homogeneous plane-wave approximation that gives greater stability to the migration technique, allowing a stable and satisfactory migration of the GPR data up to depths equivalent to three times the skin depth computed at the dominant frequency of the radar signal. Multiple slowness references of split-step migration for lossy media are introduced and compared with the single-reference slowness technique. Finally, we propose two limited gain modifications of the migration algorithm and study their usefulness when highly lossy zones are in the media.
Stoffa, P. L., M. K. Sen, R. K. Seifoullaev, R. C. Pestana, and J. T. Fokkema, Plane-wave depth migration, Geophysics, 71, S261-S272, 2006, 4 citations, doi:10.1190/1.2357832, #1839
We present fast and efficient plane-wave migration methods for densely sampled seismic data in both the source and receiver domains. The methods are based on slant stacking over both shot and receiver positions (or offsets) for all the recorded data. If the data-acquisition geometry permits, both inline and crossline source and receiver positions can be incorporated into a multidimensional phase-velocity space, which is regular even for randomly positioned input data. By noting the maximum time dips present in the shot and receiver gathers and constant-offset sections, the number of plane waves required can be estimated, and this generally results in a reduction of the data volume used for migration. The required traveltime computations for depth imaging are independent for each particular plane-wave component. It thus can be used for either the source or the receiver plane waves during extrapolation in phase space, reducing considerably the computational burden. Since only vertical delay times are required, many traveltime techniques can be employed, and the problems with multipathing and first arrivals are either reduced or eliminated. Further, the plane-wave integrals can be pruned to concentrate the image on selected targets. In this way, the computation time can be further reduced, and the technique lends itself naturally to a velocity-modeling scheme where, for example, horizontal and then steeply dipping events are gradually introduced into the velocity analysis. The migration method also lends itself to imaging in anisotropic media because phase space is the natural domain for such an analysis
Zhang, X., B. Rutt, U. Catalyurek, T. Kurc, P. L. Stoffa, M. K. Sen, and J. Saltz, Supporting scalable and distributed data subsetting and aggregation in large-scale seismic data analysis, Int. J. High Performance Computing Applications, 20, 423-438, 2006, doi:10.1177/1094342006067471, #1841
The ability to query and process very large, terabyte-scale datasets has become a key step in many scientific and engineering applications. In this paper, we describe the application of two middleware frameworks in an integrated fashion to provide a scalable and efficient system for execution of seismic data analysis on large datasets in a distributed environment. We investigate different strategies for efficient querying of large datasets and parallel implementations of a seismic image reconstruction algorithm. Our results on a state-of-the-art mass storage system coupled with a high-end compute cluster show that our implementation is scalable and can achieve about 2.9 Gigabytes per second data processing rate â about 70% of the maximum 4.2GB/s application-level raw I/O bandwidth of the storage platform.
Gai, X., J. Rungamornrat, H. Klie, W. Bangerth, M. F. Wheeler, P. L. Stoffa, M. K. Sen, and R. K. Seifoullaev, Fully integrated reservoir flow, geomechanics and seismic modelling: A tool for better reservoir characterization and geomechanical prediction using 4D seismic
, SEG Ann. Meeting, Houston, TX, RD 3.2, 1367-1370, 2005, #1958
Kurc, T., U. Catalyurek, X. Zhang, J. Saltz, R. Martino, M. F. Wheeler, M. Peszynska, A. Sussman, C. Hansen, M. K. Sen, R. K. Seifoullaev, P. L. Stoffa, C. Torres-Verdin, and M. Parashar, A simulation and data analysis system for large scale data driven, oil reservoir simulation studies, Concurrency and Computation: Practice and Experience, 17, 1441-1467, 2005, 5 citations, doi:10.1002/cpe.898, #1728
The main goal of oil reservoir management is to provide more efficient, cost-effective and environmentally safer production of oil from reservoirs. Numerical simulations can aid in the design and implementation of optimal production strategies. However, traditional simulation-based approaches to optimizing reservoir management are rapidly overwhelmed by data volume when large numbers of realizations are sought using detailed geologic descriptions. In this paper, we describe a software architecture to facilitate large-scale simulation studies, involving ensembles of long-running simulations and analysis of vast volumes of output data.
Matossian, V., V. Bhat, M. Parashar, M. Peszynska, M. K. Sen, P. L. Stoffa, and M. F. Wheeler, Autonomic oil reservoir optimization on the grid, Concurrency and Computation: Practice and Experience, 17, 1-26, 2005, 10 citations, doi:10.1002/cpe.871, #1721
The emerging Grid infrastructure and its support for seamless and secure interactions is enabling a new generation of autonomic applications where the application components, Grid services, resources, and data interact as peers to manage, adapt and optimize themselves and the overall application. In this paper we describe the design, development and operation of a prototype of such an application that uses peer-to-peer interactions between distributed services and data on the Grid to enable the autonomic optimization of an oil reservoir. Copyright © 2005 John Wiley & Sons, Ltd.
Mukherjee, A., M. K. Sen, and P. L. Stoffa, Traveltime computation and pre-stack time migration in transversely isotropic media, J. Seismic Explor., 13, 201-225, 2005, 1 citation, #1793
Roy, L., M. K. Sen, D. D. Blankenship, P. L. Stoffa, and T. G. Richter, Inversion and uncertainty estimation of gravity data using simulated annealing: An application over Lake Vostok, East Antarctica, Geophysics, 70, J1-J12, 2005, 12 citations, doi:10.1190/1.1852777, #1664
Interpretation of gravity data warrants uncertainty estimation because of its inherent nonuniqueness. Although the uncertainties in model parameters cannot be completely reduced, they can aid in the meaningful interpretation of results. Here we have employed a simulated annealing (SA)âbased technique in the inversion of gravity data to derive multilayered earth models consisting of two and three dimensional bodies. In our approach, we assume that the density contrast is known, and we solve for the coordinates or shapes of the causative bodies, resulting in a nonlinear inverse problem. We attempt to sample the model space extensively so as to estimate several equally likely models. We then use all the models sampled by SA to construct an approximate, marginal posterior probability density function (PPD) in model space and several orders of moments. The correlation matrix clearly shows the interdependence of different model parameters and the corresponding trade-offs. Such correlation plots are used to study the effect of a priori information in reducing the uncertainty in the solutions. We also investigate the use of derivative information to obtain better depth resolution and to reduce underlying uncertainties. We applied the technique on two synthetic data sets and an airborne-gravity data set collected over Lake Vostok, East Antarctica, for which a priori constraints were derived from available seismic and radar profiles. The inversion results produced depths of the lake in the survey area along with the thickness of sediments. The resulting uncertainties are interpreted in terms of the experimental geometry and data error.
Roy, L., M. K. Sen, K. D. McIntosh, P. L. Stoffa, and Y. Nakamura, Joint inversion of first arrival seismic trave- time and gravity data, J. Geophysics Engin., 2, 277-289, 2005, 6 citations, doi:10.1088/1742-2132/2/3/011, #1764
oint interpretation of disparate geophysical datasets can be a powerful tool in deriving meaningful subsurface geologic models. Here we describe a method for joint inversion of first arrival travel-time and gravity data with application to field data from a geologically complex subduction zone. As both tomographic and gravity inversions are prone to non-uniqueness, incorporation of prior information in the model description is crucial to the success of such algorithms. We employ a layer-based model description, in which interfaces (which may also be called iso-velocity lines) are defined by a summation of arc-tangent functions. Arc-tangent functions are highly flexible in mapping smooth interfaces as well as the nearly discontinuous changes in depth of an interface. Within each layer, the velocity is assumed to vary linearly with depth at each surface location. Because of the non-uniqueness of the gravity inversion, we use prior knowledge to relate the velocity to density values. In our application here, the density is related to the velocity using a fourth-order polynomial whose coefficients are assumed to be known. The nonlinear optimization problem is solved by a very fast simulated annealing (VFSA) technique. At each iteration, travel times are generated by the solution of the Eikonal equation while the gravity anomalies are computed using a standard formula. The objective function consists of two parts: one measures the misfit in travel time and the other measures the misfit of gravity anomalies. We applied our technique to field data collected over the Ryukyu subduction zone offshore Taiwan during an ocean bottom seismometer (OBS) experiment (called TAICRUST) conducted in the year 1995. Application to one NS and one EW trending line resulted in acceptable fit of both travel-time and gravity data. The resulting models are helpful in the interpretation of the local geology.
Sena, A. R., P. L. Stoffa, and M. K. Sen, Migration of ground penetrating radar data in heterogeneous and dispersive media, in New Strategies for European Remote Sensing, Proc. 24th Symp. Eur. Assoc. Remote Sensing Laboratories, Dubrovnik, Croatia, 25-27 May, 2004, 711-719, 2005, #1913
Stoffa, P. L., M. K. Sen, R. K. Seifoullaev, H. Klie, X. Gai, W. Bangerth, J. Rungamornrat, and M. F. Wheeler, An analysis of flow-simulation scales and seismic response, SEG Ann. Meeting, Houston, TX, RC P2.5, 1461-1465, 2005, #1959
Xia, Y. L., G. Kocurek [], P. L. Stoffa, and M. K. Sen, Using different hydrological variables to assess the impacts of atmospheric forcing errors on optimization and uncertainty analysis of the CHASM surface model at a cold catchment, J. Geophys. Res., 110, D01101, 2005, 2 citations, doi:10.1029/2004JD005130, #1735
Estimation of parameters for land-surface models, along with their corresponding uncertainties, relies on the input data for the atmospheric forcing variables including atmospheric pressure, temperature, humidity, wind speed, precipitation, and incoming shortwave and longwave radiation. Most studies assume that forcing data are accurate and contain no random or systematic observational errors. In practice, there are indeed systematic errors in precipitation measurements, especially for snowfall, due to wind-caused undercatch. Incoming shortwave and longwave radiation fluxes are often not directly measured, but estimated from empirical formulations. Uncertainties in these forcing data may substantially affect optimization and uncertainty estimates of land surface models. In this study, we used 18-year forcing and calibration data as well as information about the uncertainties in the forcing variables at Valdai, Russia, to study the impacts of forcing errors on selection of optimal model parameters and their uncertainty estimates when three different hydrological variables were used for calibration. The results show that forcing errors have few effects on the selection of optimal model parameter sets when monthly evapotranspiration and runoff are calibrated. However, forcing errors do introduce significant effects on the selection of optimal model parameters when daily snow water equivalent is calibrated. Forcing errors also significantly affect uncertainty estimates of the land surface model parameters. In addition, constraints of forcing errors are different when different hydrological variables are calibrated. All three hydrological variables constrain the incoming longwave radiation error well, and the snow water equivalent and runoff constrain winter snowfall errors well. However, all three hydrological variables cannot constrain the incoming solar radiation error well. We highlight in this study that runoff is shown to be a good observable to use for calibration, the reason being that it integrates multiple hydrological processes; and the results support the theory that typical rain/snow gauges have 10â20% undercatch.
Xia, Y. L., G. Kocurek [], P. L. Stoffa, and M. K. Sen, Optimal parameter and uncertainty estimation of a surface model: Sensitivity to parameter ranges and model complexities, Adv. Atmos. Sci., 22, 142-157, 2005, 1 citation, #1794
Aldunate, G. C., R. C. Pestana, and P. L. Stoffa, Migracao sismica 2-D pre-empilhamento em profundidade com operadores de extrapolacao "split-step", (in portuguese), Revista Brasileira de Geofisica, 22(2), 153-161, 2004, doi:10.1590/S0102-261X2004000200005, #1796
Three 2D prestack depth migration techniques using split-step extrapolation operators were developed and tested on seismic data sorted into common shot gathers. In the first method, which we name simultaneous split-step migration (SS-S), the migration procedure is carried out simultaneously for the sources and receivers. The recorded receiver data are depropagated in depth and the source wavefield is downward propagated using the split-step operators for both. The final depth section is achieved by summing all the frequencies of interest after the correlation of the propagated and depropagated wavefields, for each depth level and by the sum of all migrated shot gathers. To decrease the computational time of the SS-S method, we can calculate the source wavefield's through a finite difference solution of the eikonal equation. This second method we call the hybrid split-step migration method (SS-H). In the third migration method, we combine the SS-S with the PSPI method. In this case the wavefields are depropagated using split-step operators for different velocities and then interpolated as in the PSPI method. We called this method PSPI-SS. The choice of the split-step operator for migration is mainly due to its easy implementation, high accuracy and robustness even in situations with very strong lateral velocity variation. The results we present in this work were obtained using the Marmousi data and also the SEG-EAGE salt model, which present very high geologic complexity. The results obtained with the three differentes methods were compared and all show satisfactory images.
Bangerth, W., H. Klie, M. F. Wheeler, P. L. Stoffa, and M. K. Sen, On optimization algorithms for the reservoir oil well placement problem, Computational Geosciences, 10, 303-319, 2004, 7 citations, doi:10.1007/s10596-006-9025-7, #1750
Determining optimal locations and operation parameters for wells in oil and gas reservoirs has a potentially high economic impact. Finding these optima depends on a complex combination of geological, petrophysical, flow regimen, and economical parameters that are hard to grasp intuitively. On the other hand, automatic approaches have in the past been hampered by the overwhelming computational cost of running thousands of potential cases using reservoir simulators, given that each of these runs can take on the order of hours. Therefore, the key issue to such automatic optimization is the development of algorithms that find good solutions with a minimum number of function evaluations. In this work, we compare and analyze the efficiency, effectiveness, and reliability of several optimization algorithms for the well placement problem. In particular, we consider the simultaneous perturbation stochastic approximation (SPSA), finite difference gradient (FDG), and very fast simulated annealing (VFSA) algorithms. None of these algorithms guarantees to find the optimal solution, but we show that both SPSA and VFSA are very efficient in finding nearly optimal solutions with a high probability. We illustrate this with a set of numerical experiments based on real data for single and multiple well placement problems.
Jackson, C. S., M. K. Sen, and P. L. Stoffa, An efficient stochastic Bayesian approach to optimal parameter and uncertainty estimation for climate model predictions, J. Climate, 17, 2828-2841, 2004, 21 citations, #1626
One source of uncertainty for climate model predictions arises from the fact that climate models have been optimized to reproduce observational means. To quantify the uncertainty resulting from a realistic range of model configurations, it is necessary to estimate a multidimensional probability distribution that quantifies how likely different model parameter combinations are, given knowledge of the uncertainties in the observations. The computational cost of mapping a multidimensional probability distribution for a climate model using traditional means (e.g., Monte Carlo or Metropolis/Gibbs sampling) is impractical, requiring 104â106 model evaluations for problems involving less than 10 parameters. This paper examines whether such a calculation is more feasible using a particularly efficient but approximate algorithm called Bayesian stochastic inversion, based on multiple very fast simulated annealing (VFSA). Investigated here is how the number of model parameters, natural variability, and the degree of nonlinearity affect the computational cost and accuracy of estimating parameter uncertainties within a surrogate climate model that is able to approximate the noise and response behavior of a realistic atmospheric GCM. In general, multiple VFSA is one to two orders of magnitude more efficient than the Metropolis/Gibbs sampler, depending primarily on dimensionality of the parameter space analysis. The average cost of estimating parameter uncertainties is only moderately affected by noise within the model as long as the signal-to-noise ratio is greater than 5. Also the average cost of estimating parameter uncertainties nearly doubles for problems in which parameters are nonlinearly related.
Mu, Q., C. S. Jackson, and P. L. Stoffa, A multivariate empirical-orthogonal-function-based measure of climate model performance, J. Geophys. Res., 109, D15101, 2004, 3 citations, doi:10.1029/2004JD004584, #1734
A measure of the average distance between climate model predictions of multiple fields and observations has been developed that is based on the use of empirical orthogonal functions (EOFs). The application of EOFs provides a means to use information about spatial correlations in natural variability to provide a more balanced view of the significance of changes in model predictions across multiple fields, seasons, and regions. A comparison is made between the EOF-based measure and measures that are normalized by grid point variance and spatial variance for changes in the National Center for Atmospheric Research Community Climate Model, Version 3.10 (CCM3.10), parameter controlling initial cloud downdraft mass flux (ALFA), an important parameter within the Zhang and McFarlane [1995] convection scheme. All measures present consistent views that increasing ALFA from its default value creates significant improvements in precipitation, shortwave radiation reaching the surface, and surface latent heat fluxes at the expense of degrading predictions of total cloud cover, near-surface air temperature, net shortwave radiation at the top of the atmosphere, and relative humidity. However, the relative importance of each of these changes, and therefore the average view of the change in model performance, is significantly impacted by the details of how each measure of model performance handles regions with little or no internal variability. In general, the EOF-based measure emphasizes regions where modeled-observational differences are large, excluding those regions where internal variability is small.
Tsoflias, G., J.-P. Van Gestel, P. L. Stoffa, D. D. Blankenship, and M. K. Sen, Vertical fracture detection by exploiting the polarization properties of ground-penetrating radar signals, Geophysics, 69, 803-810, 2004, 14 citations, doi:10.1190/1.1759466, #1622
Vertically oriented thin fractures are not always detected by conventional single-polarization reflection profiling ground-penetrating radar (GPR) techniques. We study the polarization properties of EM wavefields and suggest multipolarization acquisition surveying to detect the location and azimuth of vertically oriented fractures. We employ analytical solutions, 3D finite-difference time-domain modeling, and field measurements of multipolarization GPR data to investigate EM wave transmission through fractured geologic formations. For surface-based multipolarization GPR measurements across vertical fractures, we observe a phase lead when the incident electric-field component is oriented perpendicular to the plane of the fracture. This observation is consistent for nonmagnetic geologic environments and allows the determination of vertical fracture location and azimuth based on the presence of a phase difference and a phase lead relationship between varying polarization GPR data.
Xia, Y. L., P. L. Stoffa, C. S. Jackson, and M. K. Sen, Effect of forcing data errors on calibration and uncertainty estimates of the CHASM model: A multi-dataset study, in Observations, Theory, and Modeling of Atmospheric and Oceanic Variability, World Sci. Ser. Meteor. Asia, Singapore, 3, 340-355, 2004, #1669
Xia, Y. L., G. Kocurek [], C. S. Jackson, P. L. Stoffa, and M. K. Sen, Impacts of data length on optimal parameter and uncertainty estimation of a land surface model, J. Geophys. Res., 109, D07101, 2004, 12 citations, doi:10.1029/2003JD004419, #1719
The optimal parameters and uncertainty estimation of land surface models require that appropriate length of forcing and calibration data be selected for computing error functions. Most of the previous studies used less than two years of data to optimize land surface models. In this study, 18-year hydrometeorological data at Valdai, Russia, were used to run the Chameleon Surface Model (CHASM). The optimal parameters were obtained by employing a global optimization technique called very fast simulated annealing. The uncertainties of model parameters were estimated by the Bayesian stochastic inversion technique. Forty-four experiments were conducted by using different lengths of data from the 18-year record, and a total of about 3 million parameter sets were produced. This study found that different calibration variables require different lengths of data to obtain optimal parameters and uncertainty estimates which are insensitive to the period selected. In the case of optimal parameters, monthly root-zone soil moisture, runoff, and evapotranspiration require 8, 3, and 1 years of data, respectively. In the case of uncertainty estimates, monthly root-zone soil moisture, runoff, and evapotranspiration require 8, 8, and 3 years of data, respectively. Spin-up has little impact on the selection of optimal parameters and uncertainty estimates when evapotranspiration and runoff were calibrated. However, spin-up affects the selection of optimal parameters when soil moisture was calibrated.
Xia, Y. L., M. K. Sen, C. S. Jackson, and P. L. Stoffa, Multidataset study of optimal parameter and uncertainty estimation of a land surface model with Bayesian stochastic inversion and multicriteria method, J. Applied Meteor., 43, 1477-1497, 2004, 4 citations, #1720
This study evaluates the ability of Bayesian stochastic inversion (BSI) and multicriteria (MC) methods to search for the optimal parameter sets of the Chameleon Surface Model (CHASM) using prescribed forcing to simulate observed sensible and latent heat fluxes from seven measurement sites representative of six biomes including temperate coniferous forests, tropical forests, temperate and tropical grasslands, temperate crops, and semiarid grasslands. Calibration results with the BSI and MC show that estimated optimal values are very similar for the important parameters that are specific to the CHASM model. The model simulations based on estimated optimal parameter sets perform much better than the default parameter sets. Cross-validations for two tropical forest sites show that the calibrated parameters for one site can be transferred to another site within the same biome. The uncertainties of optimal parameters are obtained through BSI, which estimates a multidimensional posterior probability density function (PPD). Marginal PPD analyses show that nonoptimal choices of stomatal resistance would contribute most to model simulation errors at all sites, followed by ground and vegetation roughness length at six of seven sites. The impact of initial root-zone soil moisture and nonmosaic approach on estimation of optimal parameters and their uncertainties is discussed.
Ahmed, I., P. L. Stoffa, and M. K. Sen, Residual migration velocity analysis in the offset-depth domain, J. Seismic Explor., 12, 237-257, 2003, #1722
Jackson, C. S., Y. L. Xia, M. K. Sen, and P. L. Stoffa, Optimal parameter and uncertainty estimation of a land surface model: A case example using data from Cabauw, Netherlands, J. Geophys. Res., 108, 4853, 2003, 18 citations, doi:10.1029/2002JD002991, #1627
Land surface models involve a large number of interdependent parameters that affect the physics of how surface energy fluxes are partitioned between latent heat, sensible heat, net radiative, and ground heat fluxes. The goal of an optimal parameter and uncertainty analysis of a land surface model is to identify a range of parameter sets that enable model predictions to be bounded within observational uncertainties. Here we apply Bayesian stochastic inversion (BSI) using very fast simulated annealing (VFSA) to identify parameter sets of the Chameleon surface model (CHASM) land surface model that are consistent with the uncertainty limits ascribed to a high-quality data set collected from Cabauw, Netherlands. These results are compared to the parameter sets obtained through the multicriteria (MC) approach. All analyses evaluate model performance against daily and monthly mean observations of sensible, latent, and ground heat fluxes. BSI and MC identify similar âbest fitâ model parameter sets that improve CHASM performance over default parameter settings. The three most important CHASM parameters at Cabauw are minimum stomatal resistance, vegetation roughness length, and vegetation fraction cover. BSI is based on a Bayesian inference model such that that it expresses uncertainty in terms of a posterior probability density function, different moments of which provide information about parameter means and covariances. Although MC gives a range of possible optimal parameters through the concept of a Pareto set, we found that these ranges did not provide a consistent or representative view of the uncertainty within the observational data. The BSI algorithm in the current study is particularly efficient in that it only requires about double the number of model evaluations than the MC algorithm. This is a substantial saving over other more accurate methods to evaluate uncertainty such as the Metropolis/Gibbs' sampler that requires at least 40 times more computations than the BSI algorithm to obtain similar results.
Sen, M. K., P. L. Stoffa, R. K. Seifoullaev, and J. T. Fokkema, Numerical and field investigations of GPR: Toward an airborne GPR, Subsurface Sensing Technologies and Applications, 4, 41-60, 2003, doi:10.1023/A:1023011413969, #1672
The main objective of our work is to explore the feasibility of designing an airborne RADAR system much like the University of Texas Institute for Geophysics airborne ice-penetrating radar systems for use in the detection of underground structures. Here we summarize our initial efforts toward achieving this goal. The studies conducted so far include two major tasks: (1) development of hybrid numerical-analytical techniques for simulating the radar response of underground targets, and (2) a ground penetrating radar field experiment conducted at a known site that has underground tunnels. In our numerical simulation, we used analytic solutions to propagate the EM wave field from an airplane to the ground surface and then used a finite difference scheme to compute the response of complex geology including some underground structure. During the field investigation, we collected several surface GPR profiles with varying source-receiving antenna separations and configurations. The data were processed using algorithms similar to those used in seismic data processing. Despite the high moisture content in the field area during data collection, we were able to identify features in the data that can be modeled as diffractions from underground targets. The results are encouraging. Our future work will include revisiting the site during the summer months since we expect to be able to image the underground facilities better when the ground is dry. Based on the success of the ground based GPR field surveys, we will investigate the possibilities of designing and testing an airborne system.
Jiao, J., P. L. Stoffa, M. K. Sen, and R. K. Seifoullaev, Residual migration-velocity analysis in the plane-wave domain, Geophysics, 67, 1258-1269, 2002, 8 citations, doi:10.1190/1.1500388, #1569
Over the last few years, migration-velocity analysis methods have been developed for 2-D and 3-D models by extending the assumptions and approximations used for rms velocity models. Computational requirements for these analyses have increased dramatically because top-down layer-stripping migration is needed to derive interval velocities directly instead of using rms velocities and then converting into interval velocities. We establish exact equations for 1-D and 2-D residual velocity analysis in the depth-plane-wave domain and use these in an iterative and interactive migration velocity analysis program. The new method updates interval velocities directly in a top-down residual-difference correction for all layers after prestack depth migration instead of top-down layer-stripping migration followed by residual analysis. This makes the new method a suitable tool for migration velocity analysis, especially for 3-D surveys. We test the method on synthetic and field data. The field data results show that a reasonable velocity model is obtained and most common image gathers are correctly imaged using no more than four iterations.
Muszala, S. P., P. L. Stoffa, and L. A. Lawver, An application for removing cultural nose from aeromagnetic data, Geophysics, 66, 213-219, 2001, 2 citations, doi:10:1190/1.1444897, #1515
A high-resolution aeromagnetic survey collected over the North Slope of Alaska by World Geoscience was acquired in response to a need for highly detailed data in an area where traditional geophysical techniques are expensive and prohibitive (McConnell, 1995) (Figure 1). These data were recently released to the University of Texas' Institute for Geophysics and provide a unique opportunity to investigate the problem of cultural noise suppression in aeromagnetic data. The data contain isolated magnetic anomalies that are presumably from the many drill platforms and their accompanying cultural objects such as buildings and pipe repositories (Figure 2). We present a new, automated method to reduce the amplitude of these cultural anomalies without affecting the magnetic signal from the surrounding geology.
Pestana, R. C., and P. L. Stoffa, Plane wave prestack time migration, J. Seismic Explor., 9, 211-222, 2001, #1517
Van Gestel, J.-P., and P. L. Stoffa, Application of Alford rotation to ground-penetrating radar data, Geophysics, 66, 1781-1792, 2001, 9 citations, doi:10:1190/1.1487120, #1510
We investigate the application of Alford rotation to ground-penetrating radar (GPR) data. By recording the reflected field amplitudes using four different configurations, we extract information about the orientation of buried objects that have angle-dependent reflectivity. In theory this method can be successfully applied to find the orientation of dipping layers, cylinders, and vertical fractures. Modeling results show angle-dependent reflections in all three cases; as a result, we can exactly determine the orientation of these targets. Analysis of a field survey at a controlled GPR test site in which reflections were collected from an elongate cylinder buried in a homogeneous soil show good prediction of the angle of orientation of the cylinder and confirm the expected theoretical and modeling results. The Alford rotation method requires accurate data acquisition for effective practical implementation. Improved results will require exact knowledge of the radiation pattern of the GPR antennas under different circumstances.
Calderon-Macias, C., M. K. Sen, and P. L. Stoffa, Artificial neural networks for parameter estimation in geophysics, Geophys. Prospecting, 48, 21-47, 2000, 26 citations, doi:10.1046/j.1365-2478.2000.00171.x, #1461
Artificial neural systems have been used in a variety of problems in the fields of science and engineering. Here we describe a study of the applicability of neural networks to solving some geophysical inverse problems. In particular, we study the problem of obtaining formation resistivities and layer thicknesses from vertical electrical sounding (VES) data and that of obtaining 1D velocity models from seismic waveform data. We use a two-layer feedforward neural network (FNN) that is trained to predict earth models from measured data. Part of the interest in using FNNs for geophysical inversion is that they are adaptive systems that perform a non-linear mapping between two sets of data from a given domain. In both of our applications, we train FNNs using synthetic data as input to the networks and a layer parametrization of the models as the network output. The earth models used for network training are drawn from an ensemble of random models within some prespecified parameter limits. For network training we use the back-propagation algorithm and a hybrid back-propagationâsimulated-annealing method for the VES and seismic inverse problems, respectively. Other fundamental issues for obtaining accurate model parameter estimates using trained FNNs are the size of the training data, the network configuration, the description of the data and the model parametrization. Our simulations indicate that FNNs, if adequately trained, produce reasonably accurate earth models when observed data are input to the FNNs.
Eldholm, O., M. F. Coffin, P. L. Stoffa, and J. A. Austin, Academic/industrial cooperation in marine seismology: An international model, Proc., Offshore Tech. Conf., 32, OTC 11944, 2000, #1498
Jiao, J., P. L. Stoffa, M. K. Sen, and R. K. Seifoullaev, Residual migration velocity analysis, 70th Ann. Int. Meeting SEG, Expanded Abstracts, Calgary, Canada, MIG7.8, 2000, #1570
Liu, F., M. K. Sen, and P. L. Stoffa, Dip selective 2-D multiple attenuation in the plane-wave domain, Geophysics, 65, 264-274, 2000, 5 citations, doi:10:1190/1.1444717, #1460
In many geological settings, strong reflections at the air-water interface contribute to most of the multiple energy in the recorded seismograms. Here, we describe a method for free-surface multiple attenuation using a reflection operator model of a seismic record, derived using the well-known invariant embedding technique. We implement this method in the 2-D plane-wave domain, where lateral variation of the geological structure of the earth is taken into account by the coupling of different ray parameters. In situations where the lateral variations are smooth, the data are well compressed in the 2-D plane-wave domain and the resultant bandlimited matrices significantly reduce the computation cost. One important feature of the proposed method is its flexibility, which allows for the removal of multiples from selected reflections. To generate multiple free data, wave-theory-based multiple attenuation methods attempt to estimate either the source function or the subsurface reflectivity. Our method takes advantage of both approaches, such that we initially predict multiple traveltime using the reflectivity approach and then seek a source function to predict the amplitudes. Synthetic and real data examples show that this method is stable and successful in attenuating the surface multiples.
McIntosh, K. D., F. E. Akbar, C. Calderon-Macias, P. L. Stoffa, S. Operto, G. L. Christeson, Y. Nakamura, T. H. Shipley, E. R. Flueh, A. U. Stavenhagen, and G. Leandro, Large aperture seismic imaging at a convergent margin: Techniques and results from the Costa Rica seismogenic zone, Marine Geophysical Researches, 21, 451-474, 2000, 5 citations, doi:10.1023/A:1026597927732, #1514
In March and April 1995 a cooperative German, Costa Rican, and United States research team recorded onshore-offshore seismic data sets along the Pacific margin of Costa Rica using the R/V Ewing. Off the Nicoya Peninsula we used a linear array of ocean bottom seismometers and hydrophones (OBS/H) with onshore seismometers extending across much of the isthmus. In the central area we deployed an OBS/H areal array consisting of 30 instruments over a 9 km by 35-km area and had land stations on the Nicoya Peninsula adjacent to this marine array and also extending northeast on the main Costa Rican landmass. Our goal in these experiments was to determine the crustal velocity structure along different portions of this convergent margin and to use the dense instrument deployments to create migrated reflection images of the plate boundary zone and the subducting Cocos Plate. Our specific goal in the central area was to determine whether a subducted seamount is present at the location of the 1990, M 7 earthquake off the Nicoya Peninsula and can thus be linked to its nucleation. Subsequently we have processed the data to improve reflection signals, used the data to calculate crustal velocity models, and developed several wide-aperture migration techniques, based on a Kirchhoff algorithm, to produce reflection images. Along the northern transect we used the ocean bottom data to construct a detailed crustal velocity model, but reflections from the plate boundary and top and bottom of the subducting Cocos plate are difficult to identify and have so far produced poor images. In contrast, the land stations along this same transect recorded clear reflections from the top of the subducting plate or plate boundary, within the seismogenic zone, and we have constructed a clear image from this reflector beneath the Nicoya shelf. Data from the 3-D seismic experiment suffer from high-amplitude, coherent noise (arrivals other than reflections), and we have tried many techniques to enhance the signal to noise ratio of reflected arrivals. Due to the noise, an apparent lack of strong reflections from the plate boundary zone, and probable structural complexity, the resulting 3-D images only poorly resolve the top of the subducting Cocos Plate. The images are not able to provide compelling evidence of whether there is a subducting seamount at the 1990 earthquake hypocenter. Our results do show that OBS surveys are capable of creating images of the plate boundary zone and the subducting plate well into the seismogenic zone if coherent reflections are recorded at 1.8 km instrument spacing (2-D) and 5 km inline by 1 km crossline spacing for 3-D acquisition. However, due to typical high amplitude coherent noise, imaging results may be poorer than expected, especially in unfavorable geologic settings such as our 3-D survey area. More effective noise reduction in acquisition, possibly with the use of vertical hydrophone arrays, and in processing, with advanced multiple removal and possibly depth filtering, is required to achieve the desired detailed images of the seismogenic plate boundary zone.
Pestana, R. C., P. L. Stoffa, and J. R. Santos, Plane wave pre-stack migration, 70th Ann. Int. Meeting SEG, Expanded Abstracts, Calgary, Canada, MIG3.7, 2000, #1575
Porsani, M. J., P. L. Stoffa, M. K. Sen, and R. K. Chunduru, Fitness functions, genetic algorithms and hybrid optimization in seismic waveform inversion, J. Seismic Explor., 9, 143-164, 2000, 6 citations, #1518
Van Gestel, J.-P., and P. L. Stoffa, Migration using multi-configuration GPR data, Proc. Eighth Int. Conf. Ground Penetrating Radar, GPR2000, Gold Coast, Australia, 448-452, 2000, #1526
Xia, G. Y., M. K. Sen, and P. L. Stoffa, Mapping of elastic properties of gas hydrates in the Carolina Trough by waveform inversion, Proc. 68th Ann. Meeting, Soc. Explor. Geophys, New Orleans LA, 1146-1149, 2000, #1465
Xia, G. Y., M. K. Sen, and P. L. Stoffa, Mapping of elastic properties of gas hydrates in the Carolina trough by waveform inversion, Geophysics, 65, 735-744, 2000, 1 citation, doi:10:1190/1.1444772, #1469
Gas hydrates are frozen methane gas that forms at appropriate pressure and temperature conditions. They are found in the marine sediments along continental margins worldwide. They have the economic potential of being tapped as a fuel source and also have the potential as a "greenhouse" agent after being freed into the atmosphere. In seismic sections, the occurrence of the base of gas hydrates, in some areas is often marked by a bright amplitude reflection. Such reflections follow the sea floor topography and are called bottom-simulating reflectors (BSR). The BSRs have negative polarity with respect to the sea-floor reflection and, in a common shot or a CDP gather, the amplitude increases with offset. The negative impedance contrast causing BSRs may be due to negative velocity contrast between hydrated sediments and normal sediment below or due to the presence of free gas at the base of the hydrates. In this paper, we carry out a prestack seismic waveform inversion of multichannel seismic data collected in the offshore of South Carolina to investigate the origin of the BSRs. We apply a multistage seismic waveform inversion for this purpose. A nonlinear optimization method is applied to estimate the low-frequency component of the velocity, whereas an amplitude-variation-with-offset inversion is applied to determine high-frequency components of the velocity field. Our detailed seismic waveform inversion along the seismic line results in at least three low-velocity zones where the velocity is well below the velocity of the normal sediments. Such low-velocity zones correlate very well with negative fluid factors indicating the presence of free gas. Thus we conclude that the BSR is caused by free gas at the base of the hydrates in this region. We identify at least two other layers of free gas beneath the hydrates. The thickness of each of these layers is below the resolution of the source wavelet. Our results confirm similar findings reported from Ocean Drilling Program drilling and vertical-seismic-profiling analysis in the same general area.
Liu, F., M. K. Sen, and P. L. Stoffa, Surface multiple attenuation for multi-component ocean bottom seismometer data, Proc. 69th Ann. Meeting, Soc. Explor. Geophys., Houston, TX, 1999, #1572
Pestana, R. C., P. L. Stoffa, and M. K. Sen, Multiple attenuation in the plane-wave domain by match filtering, J. Seismic Explor., 8, 167-179, 1999, #1462
Sen, V., M. K. Sen, and P. L. Stoffa, PVM based 3-D Kirchhoff depth migration using dynamically computed travel-times: An application in seismic data processing, Parallel Computing, 25, 231-248, 1999, 2 citations, doi:10.1016/S0167-8191(98)00113-6, #1411
Seismic depth migration is an image reconstruction technique used to generate realistic images of the Earth's interior from surface recordings of sound waves reflected by buried geological structures. Migration algorithms that have been developed for digitally recorded seismic data are computationally intensive and are typically implemented on massively parallel architecture. In this paper we describe a highly accurate method of computing seismic travel-times and study the feasibility of a PVM-based SPMD parallelization of a migration algorithm based on this method. The parallel algorithm developed in this study is coarse grained and utilizes a simple geometrical decomposition of the input problem. We show that near linear speedup can be achieved on a homogeneous cluster and the algorithm can be ported to a variety of platforms, i.e. clusters of Sun workstations, Cray Y-MP, nCUBE-2, Cray T3E and SGI Origin 2000. The performance of the scheme on different computational platforms ultimately depends on a number of different factors, e.g., machine architecture, machine speed, network configuration etc. For a simple test case that was studied in this research, the performances of modest PVM clusters were found to be encouraging. This underscores the advantages of adopting a strategy based on parallelizing computationally intensive tasks and also demonstrates the suitability of the travel-time scheme for developing coarsely parallel applications.
Akbar, F. E., P. L. Stoffa, and M. K. Sen, Three-dimensional plane-wave Kirchoff depth migration, Proc. 68th Ann. Meeting, Soc. Explor. Geophys., New Orleans, LA, 1708-1711, 1998, #1466
Akbar, F. E., C. Calderon-Macias, V. Sen, M. K. Sen, and P. L. Stoffa, A comparative study of first arrival time computation for 3-D inhomogeneous isotropic velocity models, Proc. 68th Ann. Meeting, Soc. Explor. Geophys., New Orleans, LA, 1728-1731, 1998, #1467
Anandakrishnan, S., D. D. Blankenship, R. B. Alley, and P. L. Stoffa, Influence of subgalcial geology on the position of a West Antarctic ice stream from seismic observations, Nature, 394, 62-65, 1998, 99 citations, doi:10.1038/27889, #1400
Ice streams drain much of the interior West Antarctic Ice Sheet and buffer the main ice reservoir from oceanic influences1,2. The slow-flowing interior feeds the floating Ross Ice Shelf with ice via fast-flowing ice streams3 that are believed to modulate sea-level change through their control of inland ice storage. Understanding ice-stream behaviour, and predicting the response to climate change4, requires a better knowledge of the subglacial geology5,6. It is known that a thawed ice-bed and high-pressure basal water are necessary, but not sufficient, conditions to cause ice streaming7,8. Moreover, it has been hypothesized that a soft sedimentary bed is also required, because of its intrinsic low frictional resistance to flow9, and owing to its high erodibility so as to generate till that can deform and lubricate ice motion10,11, or to bury rough features and smooth the bed for sliding. Here we use seismic observations to provide evidence that one margin of the upglacier part of an ice stream is directly above the boundary of a basin with such sedimentary fill. The ice stream is within the basin and the ice outside the basin is slow-flowing. The basin fill presents an order-of-magnitude lower frictional resistance to ice flow than the subglacial material outside the basin. We conclude that the ice stream position is dependent on subglacial geology.
Calderon-Macias, C., M. K. Sen, and P. L. Stoffa, Automatic NMO correction and velocity estimation by a feedforward neural network, Geophysics, 63, 1696-1707, 1998, 11 citations, doi:10:1190/1.444465, #1401
We describe a new method of automatic normal moveout (NMO) correction and velocity analysis that combines a feedforward neural network (FNN) with a simulated annealing technique known as very fast simulated annealing (VFSA). The task of the FNN is to map common midpoint (CMP) gathers at control locations along a 2-D seismic line into seismic velocities within predefined velocity search limits. The network is trained while the velocity analysis is performed at the selected control locations. The method minimizes a cost function defined in terms of the NMO-corrected data. Network weights are updated at each iteration of the optimization process using VFSA. Once the control CMP gathers have been properly NMO corrected, the derived weights are used to interpolate results at the intermediate CMP locations. In practical situations in which lateral velocity variations are expected, the method is applied in spatial data windows, each window being defined by a separate FNN. The method is illustrated with synthetic data and a real marine data set from the Carolina Trough area.
Coffin, M. F., O. Eldholm, P. L. Stoffa, and J. A. Austin, Looking ahead to the future of marine reflection seismology, Eos, Trans. Amer. Geophys. Un., 79, 614-615, 1998, doi:10.1029/98EO00442, #1418
Liu, F., M. K. Sen, and P. L. Stoffa, 2-D multiple attenuation operators in τ-P domain, Proc. 68th Ann. Meeting, Soc. Explor. Geophys, New Orleans, LA, 1256-1259, 1998, #1464
Luhurbudi, E. C., J. Pulliam, J. A. Austin, S. Saustrup, and P. L. Stoffa, Removal of diurnal tidal effects from an ultra-high resolution 3-D marine seismic survey on the continental shelf offshore New Jersey, Geophysics, 63, 1036-1040, 1998, 2 citations, doi:10:1190/1.1444381, #1277
An ultra-high-resolution 3-D, single-channel seismic survey was performed off the coast of New Jersey in 1993 to study the late Quaternary history of sedimentation on the northwest Atlantic continental margin (see Davies et al., 1992) as a part of the Office of Naval Research STRATAFORM initiative (Nittrouer and Kravitz, 1995). Three different sets of profiles were acquired (Figure 1), but only the set with highest spatial density is discussed here. A single ten-element receiver recorded 300 ms of data for every shot during the survey, which covers a total area of 0.6 km (north-south) Ã 7.75 km (east-west) (see Table 1). The deep-towed HuntecTM source (deployed at ~30 m depth) produced frequencies of 500 to 3500 Hz; a band-pass filter with corner frequencies at 1000 and 3500 Hz was applied during preprocessing
Sen, M. K., F. Liu, P. L. Stoffa, and J. T. Fokkema, A unified treatment of free-surface multiple-elimination methods, J. Seismic Explor., 7, 129-143, 1998, 4 citations, #1410
Sen, V., P. L. Stoffa, I. W. D. Dalziel, D. D. Blankenship, A. M. Smith, and S. Anandakrishnan, Seismic surveys of central West Antarctica: Data and processing examples from the ANTALITH field tests (1994-1995), Terra Antartica, 5, 761-772, 1998, #1447
Sen, M. K., P. L. Stoffa, and G. Y. Xia, AVO and seismic waveform inversion in the plane wave domain: Application to gas hydrate data, Geohorizons, 3 (2), 5-14, 1998, #1463
Tanis, M. C., P. L. Stoffa, and R. C. Pestana, Prestack depth migration in the source-offset domain, Proc. 68th Ann. Meeting, Soc. Explor. Geophys., New Orleans LA, 1839-1842, 1998, #1470
Varela, C. L., P. L. Stoffa, and M. K. Sen, Background velocity estimation using non-linear optimization for reflection tomography and migration misfit, Geophys. Prospecting, 46, 51-78, 1998, 12 citations, doi:10.1046/j.1365-2478.1998.720308.x, #1188
We show that it is possible to estimate the background velocity for prestack depth migration in 2D laterally varying media using a non-linear optimization technique called very fast simulated annealing (VFSA). We use cubic splines in the velocity model parametrization and make use of either successive pairs of shot gathers or several constant-offset sections as input data for the inversion. A Kirchhoff summation scheme based on first-arrival traveltimes is used to migrate/model the input data during the velocity analysis. We evaluate and compare two different measures of error. The first is defined in the recorded data or (x,t) domain and is based on a reflection-tomography criterion. The second is defined in the migrated data or (x,z) domain and is based on a migration-misfit criterion. Depth relaxation is used to improve the convergence and quality of the velocity analysis while simultaneously reducing the computational cost. Further, we show that by coarse sampling in the offset domain the method is still robust.
Our non-linear optimization approach to migration velocity analysis is evaluated for both synthetic and real seismic data. For the velocity-analysis method based on the reflection-tomography criterion, traveltimes do not have to be picked. Similarly, the migration-misfit criterion does not require that depth images be manually compared. Interpreter intervention is required only to restrict the search space used in the velocity-analysis problem. Extension of the proposed schemes to 3D models is straightforward but practical only for the fastest available computers.
Xia, G. Y., M. K. Sen, and P. L. Stoffa, 1-D elastic waveform inversion: A divide-and-conquer approach, Geophysics, 63, 1670-1684, 1998, 7 citations, doi:10:1190/1.444463, #1412
Subsurface rock properties are manifested in seismic records as variations in traveltimes, amplitudes, and waveforms. It is commonly acknowledged that traveltimes are sensitive to the long wavelength part of the velocity, whereas amplitudes are sensitive to the short wavelength part of the velocity. The inherent sensitivity of seismic velocity at different wavelengths suggests an approach that decomposes the waveform data into traveltime and amplitude components. Therefore we propose a divide-and-conquer approach to the elastic waveform inversion problem. We first estimate the smoothly varying background velocity from the traveltime and the rapidly changing perturbations from the amplitude by amplitude variation with offset (AVO) inversion based on linearized reflection coefficient. Then we combine the perturbation with the background to obtain a starting model to be used in the final waveform inversion that models all converted waves and internal multiples assuming a 1-D earth model. For estimating the background velocity, we use the flatness of events as the objective criterion, and simulated annealing as a search tool. Three different model parameterization schemes (constant velocity blocks, splines, and arctangent models) are compared, with the arctangent having the most flexibility and least artifacts. Having obtained the background velocities, we analyze the AVO effects to estimate the perturbations to the background, for which we use a linearized inversion method. The combination of the perturbation and background should be sufficiently close to the true model so that the inverse problem becomes quasi-linear. A full elastic waveform inversion is used to fine-tune the combined model to obtain P-wave and S-wave velocity and density, again using either a nonlinear optimization method or an iterative linearized solution. Application of the inversion algorithm to synthetic data from an 84-layer model was able to predict the full reflectivity data and recover the true model parameters. Application to one seismic line in the Carolina Trough area found a thin gas zone which produces strong Bottom Simulating Reflectors (BSRs).
Calderon-Macias, C., M. K. Sen, and P. L. Stoffa, Hopfield neural networks and mean field annealing for seismic deconvolution and multiple attenuation, Geophysics, 62, 992-1002, 1997, 8 citations, doi:10:1190/1.1444205, #1216
We describe a global optimization method called mean field annealing (MFA) and its application to two basic problems in seismic data processing: Seismic deconvolution and surface related multiple attenuation. MFA replaces the stochastic nature of the simulated annealing method with a set of deterministic update rules that act on the average value of the variables rather than on the variables themselves, based on the mean field approximation. As the temperature is lowered, the MFA rules update the variables in terms of their values at a previous temperature. By minimizing averages, it is possible to converge to an equilibrium state considerably faster than a standard simulated annealing method. The update rules are dependent on the form of the cost function and are obtained easily when the cost function resembles the energy function of a Hopfield network. The mapping of a problem onto a Hopfield network is not a precondition for using MFA, but it makes analytic calculations simpler. The seismic deconvolution problem can be mapped onto a Hopfield network by parameterizing the source and the reflectivity in terms of binary neurons. In this context, the solution of the problem is obtained when the neurons of the network reach their stable states. By minimizing the cost function of the network with MFA and using an appropriate cooling schedule, it is possible to escape local minima. A similar idea can also be applied to design an operator that attenuates surface related multiple reflections from plane-wave transformed seismograms assuming a 1-D earth. The cost function for the multiple elimination problem is based on the criterion of minimum energy of the multiple suppressed data.
Chunduru, R. K., M. K. Sen, and P. L. Stoffa, Hybrid optimization methods for geophysical inversion, Geophysics, 62, 1196-1207, 1997, 33 citations, doi:10:1190/1.1444220, #1299
Local and global optimization algorithms are used commonly in geophysical data inversion. Each type of algorithm has unique advantages and disadvantages. Here we propose several methods of combining the two algorithms such that we can overcome their drawbacks and make use of the salient features of the two methods. In particular, we combined a local conjugate gradient (CG) method with a global very fast simulated annealing (VFSA) approach to solve problems of geophysical interests. We conducted a systematic study to find an efficient strategy to combine CG and VFSA optimization schemes and recommend a couple of ways for future implementations. Seven different hybrid algorithms were first tested on a set of field 1-D Schlumberger resistivity sounding data and their performances were compared with stand-alone genetic algorithm (GA), simulated annealing, and local search algorithms. Almost all of the proposed hybrid algorithms were found to be computationally more efficient than the conventional global optimization approaches. Having found the most efficient of the hybrid approaches we apply them to the problem of seismic velocity analysis using seismograms recorded in the offset-time domain. Finally, we applied the hybrid algorithm to a 2-D field resistivity profiling data collected over a disseminated sulfide zone at Safford Arizona and compared our hybrid inversion results with the previously published results.
Tanis, M. C., and P. L. Stoffa, Parallel implementation of 3-D split-step Fourier depth migration algorithm on T3E, Proc., SEG 66th Ann. Meeting, Denver, CO, 1433-1436, 1997, #1366
Xia, G. Y., M. K. Sen, and P. L. Stoffa, AVO analysis of Mobil offshore data by a linearized inversion in the τ-p domain, in AVO Inversion of Mobil Data Set, SEG Press, edited by R. G. Keys, Chapter, 9, 167-185, 1997, #1165
Akbar, F. E., M. K. Sen, and P. L. Stoffa, Prestack plane-wave Kirchhoff migration in laterally varying media, Geophysics, 61, 1068-1079, 1996, 10 citations, doi:10:1190/1.1444028, #1164
Conventional Kirchhoff migration methods are applied to seismograms in the offset-time (x â t) domain. We describe the theory, numerical details, and examples of a prestack depth migration method in the plane-wave domain that is valid in laterally inhomogeneous media. The theory is based on a Kirchhoff-Helmholtz formulation of the wavefield and uses plane-wave-transformed shot gathers for migration. We use geometric ray theory for the source wavefield continuation operator and a plane-wave expansion of the receiver wavefield in the integrand of the Kirchhoff-Helmholtz integral. The source and receiver plane-wave traveltimes are computed at each grid point in the subsurface using a finite-difference approximation of the eikonal equation with appropriate initial and boundary conditions. We developed an efficient technique to compute imaging time by a combination of these two times. The technique allows us to design algorithms for migrating shot gather or constant ray-parameter sections efficiently. We evaluate the efficiency and accuracy of the algorithm with three sets of synthetic data examples with varying degrees of complexity. We also compare the performance of the parallel algorithms using Parallel Virtual Machine (PVM). Migration of a marine data set results in excellent images of a mud volcano and the top of an accretionary prism of an active continental margin on the Nicoya Peninsula of Costa Rica.
Akbar, F. E., P. L. Stoffa, M. K. Sen, and C. L. Varela, Automated background velocity estimation in the plane wave domain using VFSA, Proc., SEG 66th Ann. Meeting, Denver, CO, 731-734, 1996, #1369
Christeson, G. L., Y. Nakamura, K. D. McIntosh, and P. L. Stoffa, Effect of shot interval on ocean bottom seismograph and hydrophone data, Geophys. Res. Lett., 23, 3783-3786, 1996, 10 citations, #1246
Data collected by 18 ocean bottom receivers for a seismic line shot at both 50‐m (∼24 s shot interval) and 125‐m (∼58 s shot interval) shot spacing provide a direct field comparison of the effect of shot interval on marine wide‐angle seismic data. Our results indicate that both shot spacings produce high‐quality refraction data in shallow water (<1000 m) on hydrophone and vertical channel data. In deeper water, the data quality of the 50‐m line is adequate for the vertical channel, but it is often poor at large offsets for the hydrophone channel in comparison to the 125‐m shot spacing data. A theoretical model to explain these observations provides further information useful for designing an experiment using ocean‐bottom receivers.
Chunduru, R. K., M. K. Sen, and P. L. Stoffa, 2-D resistivity inversion using spline parameterization and simulated annealing, Geophysics, 61, 151-161, 1996, 16 citations, doi:10:1190/1.1443935, #1060
Successful inversion of geophysical data depends on prior information, proper choice of inversion scheme, and on effective parameterization of the model space such that the model representation is appropriate and efficient. Inversion of resistivity data has long been recognized as a nonlinear or quasi-linear problem. Traditionally, 2-D resistivity inversion has been performed by trial and error methods and with linear and iterative linear methods. The linear and iterative linear methods are limited because of the requirement of good prior knowledge of the subsurface. Unlike linear and iterative linear methods, most nonlinear inversion schemes do not depend strongly on the starting solution, but prior information helps to reduce the computational cost and to obtain geologically meaningful results. In the present study, we have applied a nonlinear optimization scheme called very fast simulated annealing (VFSA) in the inversion of 2-D dipole-dipole resistivity data to image the subsurface. Unlike Metropolis simulated annealing (SA) in which each new model is drawn from a uniform distribution, VFSA draws a model from a Cauchy-like distribution, which is also a function of a control parameter called temperature. The advantage of using such a scheme is that at high temperatures, the algorithm allows for searches far beyond the current position, while at low temperatures, it looks for improvement in the close vicinity of the current model. We have used the mean square error between the synthetics and original data as the error function to be minimized. The synthetic response for 2-D models was obtained by finite-difference modeling, and cubic splines were used to parameterize the model space to get smooth images of the subsurface and to reduce computational cost. VFSA was used to estimate the conductivity at each spline node location. The inversion was applied to various synthetic data to study the influence of the starting solution and the location of the spline nodes. Finally, we applied it to real data collected over a disseminated sulfide zone at Safford, Arizona, and compared the results with those obtained from a linearized inversion and from a model based on geologic and well-log data. The VFSA results are in good agreement with the previously published results.
Chunduru, R. K., M. K. Sen, and P. L. Stoffa, Development of efficient hybrid optimization for geophysical inversion, Proc., SEG 66th Ann. Meeting, Denver, CO, 1130-1133, 1996, #1370
Jervis, M., M. K. Sen, and P. L. Stoffa, Prestack migration velocity estimation using nonlinear methods, Geophysics, 61, 138-150, 1996, 20 citations, doi:10:1190/1.1443934, #1053
We describe here methods of estimating interval velocities based on two nonlinear optimization methods; very fast simulated annealing (VFSA) and a genetic algorithm (GA). The objective function is defined using prestack seismic data after depth migration. This inverse problem involves optimizing the lateral consistency of reflectors between adjacent migrated shot records. In effect, the normal moveout correction in velocity analysis is replaced by prestack depth migration. When the least-squared difference between each pair of migrated shots is at a minimum, the true velocity model has been found. Our model is parameterized using cubic-B splines distributed on a rectangular grid. The main advantages of using migrated data are that they do not require traveltime picking, knowledge of the source wavelet, and expensive computation of synthetic waveform data to assess the degree of data-model fit. Nonlinear methods allow automated determination of the global minimum without relying on estimates of the gradient of the objective function, the starting model, or making assumptions about the nature of the objective function itself. For the velocity estimation problem, the VFSA converges 4 to 5 times faster than the GA for both a 2-D synthetic example and a structurally complex real data example from the Gulf of Mexico. Though computationally intensive, this problem requires few model parameters, and use of a fast traveltime code for Kirchhoff migration makes the algorithm tractable for real earth problems.
Pulliam, J., J. A. Austin, E. C. Luhurbudi, S. Saustrup, and P. L. Stoffa, An ultrahigh resolution 3-D survey of the shallow subsurface on the continental shelf of New Jersey, Leading Edge, 15, 839-845, 1996, doi:10.1190/1.1437378, #1218
Regional seismic surveys have identified a wedge of late Quaternary sediment extending 150 km south from the Hudson River apron along the edge of the continental shelf off New Jersey. The bottom of the sediment wedge is defined by a prominent reflector that is assumed to be an erosional surface carved during a lowstand of sea level, probably corresponding to the Wisconsin Maximum glaciation (about 18 000 years ago). An extremely high resolution 3-D seismic reflection survey of the southern part of the wedge (Figure 1) was carried out in 1993 as part of the STRATAFORM (STRATA FORMation on Margins) initiative, funded by the U.S. Office of Naval Research.
Pulliam, J., P. L. Stoffa, E. C. Luhurbudi, S. Saustrup, and J. A. Austin, 3-D depth migration of an ultra high resolution seismic survey on New Jersey's continental shelf, Proc., SEG 66th Ann. Meeting, Denver, CO, 847-850, 1996, #1367
Sen, M. K., and P. L. Stoffa, Bayesian inference, Gibbs' sampler and uncertainty estimation in geophysical inversion, Geophys. Prospecting, 44, 313-350, 1996, 85 citations, doi:10.1111/j.1365-2478.1996.tb00152.x, #1066
The posterior probability density function (PPD), σ(m|dobs), of earth model m, where dobs are the measured data, describes the solution of a geophysical inverse problem, when a Bayesian inference model is used to describe the problem. In many applications, the PPD is neither analytically tractable nor easily approximated and simple analytic expressions for the mean and variance of the PPD are not available. Since the complete description of the PPD is impossible in the highly multi-dimensional model space of many geophysical applications, several measures such as the highest posterior density regions, marginal PPD and several orders of moments are often used to describe the solutions. Calculation of such quantities requires evaluation of multidimensional integrals. A faster alternative to enumeration and blind Monte-Carlo integration is importance sampling which may be useful in several applications. Thus how to draw samples of m from the PPD becomes an important aspect of geophysical inversion such that importance sampling can be used in the evaluation of these multi-dimensional integrals. Importance sampling can be carried out most efficiently by a Gibbs' sampler (GS). We also introduce a method which we called parallel Gibbs' sampler (PGS) based on genetic algorithms (GA) and show numerically that the results from the two samplers are nearly identical.
We first investigate the performance of enumeration and several sampling based techniques such as a GS, PGS and several multiple maximum a posteriori (MAP) algorithms for a simple geophysical problem of inversion of resistivity sounding data. Several non-linear optimization methods based on simulated annealing (SA), GA and some of their variants can be devised which can be made to reach very close to the maximum of the PPD. Such MAP estimation algorithms also sample different points in the model space. By repeating these MAP inversions several times, it is possible to sample adequately the most significant portion(s) of the PPD and all these models can be used to construct the marginal PPD, mean) covariance, etc. We observe that the GS and PGS results are identical and indistinguishable from the enumeration scheme. Multiple MAP algorithms slightly underestimate the posterior variances although the correlation values obtained by all the methods agree very well. Multiple MAP estimation required 0.3% of the computational effort of enumeration and 40% of the effort of a GS or PGS for this problem. Next, we apply GS to the inversion of a marine seismic data set to quantify uncertainties in the derived model, given the prior distribution determined from several common midpoint gathers.
Sen, M. K., and P. L. Stoffa, (Book review): Wavelets in Geophysics, edited by E. Foufoula-Georgiou and P. Kumar, Mathematics of Computation, 65, 213, 1996, #1158
Xia, G. Y., M. K. Sen, and P. L. Stoffa, Two-step velocity estimation and AVO inversion in the τ-p domain,, Proc., SEG 66th Ann. Meeting, Denver, CO, 1727-1730, 1996, #1368
Zhao, L.-S., M. K. Sen, P. L. Stoffa, and C. Frohlich, Application of very fast simulated annealing to the determination of the crustal structure beneath Tibet, Geophys. J. Int., 125, 355-370, 1996, 33 citations, doi:10.1111/j.1365-246X.1996.tb00004.x, #1157
The simulated annealing method has recently been applied to several multiparameter optimization problems, including those of geophysical inversion. A new variant of simulated annealing, called very fast simulated annealing (VFSA), overcomes some of the drawbacks of a conventional simulated annealing; it has been found to be a practical tool for geophysical inversion (Sen & Stoffa 1995). the method is particularly useful for non-linear problems with multiple parameters where a grid-search method is prohibitively expensive. Here, we have applied VFSA to retrieve the crustal structure beneath seismic stations in Tibet using teleseismic body-waveform data. Our approach is to compare the radial-component records with generalized ray synthetics directly in the time domain. For any given crustal structure, we have formulated a direct relationship between the radial- and the vertical-component signal, so that we can generate synthetic radial-component data from recorded vertical-component seismograms (Zhao & Frohlich 1996). We have tested the VFSA inversion algorithm using synthetics and confirmed that it works very well, i.e. it finds solutions very close to the true solution. Using the algorithm, we have determined the crustal structures beneath 11 stations in Tibet. From south to north, the total crustal thickness is quite uniform. Our tectonic interpretation of these crustal models is that they may represent crust from the Eurasian plate injected beneath the crust of the converging Indian plate. Certain features of the models are consistent with the presence of a small convection cell or plume beneath north-central Tibet, as suggested by Molnar (1990).
Chunduru, R. K., M. K. Sen, P. L. Stoffa, and R. Nagendra, Non-linear inversion of resistivity profiling data for some regular geometrical bodies, Geophys. Prospecting, 43, 979-1003, 1995, 18 citations, doi:10.1111/j.1365-2478.1995.tb00292.x, #1041
The inversion of resistivity profiling data involves estimation of the spatial distribution of resistivities and thicknesses of rock layers from the apparent resistivity data values measured in the field as a function of electrode separation. The drawbacks of using traditional curve-matching techniques to solve this inverse problem have been overcome by iterative linear techniques but these require good starting models even if the shape of the causative body is asssumed known. In spite of the recent developments in inversion techniques, no robust method exists for the inversion of resistivity profiling data for the simple model of dikes and spheres which are the classical models of geophysical prospecting.
We apply three different non-linear inversion schemes to invert synthetic resistivity profiling data for the classical models embedded in a uniform matrix of contrasting resistivity. The three non-linear algorithms used are called the Metropolis simulated annealing (SA), very fast simulated annealing (VFSA) and a genetic algorithm (GA). We compare the performance of the three algorithms using synthetic data for an outcropping vertical dike model. Although all three methods were successful in obtaining optimal solutions for arbitrary starting models, VFSA proved to be computationally the most efficient.
Oh, J., J. A. Austin, J. D. Phillips, M. F. Coffin, and P. L. Stoffa, Seaward-dipping reflectors offshore the southeastern United States: Seismic evidence for extensive volcanism accompanying sequential formation of the Carolina Trough and Blake Plateau Basin, Geology, 23, 9-12, 1995, 14 citations, doi:10.1130/0091-7613(1995)023<0009:SDROTS>2.3.CO;2, #1070
Deep-penetration multichannel seismic reflection profiles off the southeastern United States reveal widespread seaward-dipping reflectors (SDRs). Similar features have been imaged and sampled on other North Atlantic rifted margins, where voluminous volcanism has accompanied continental breakup. Beneath the Carolina trough are two sets of SDRs, one along a basement hinge zone and another seaward of the East Coast magnetic anomaly axis. The hinge SDRs lie beneath, and apparently developed prior to, a flood basalt that erupted at 184 ± 3 Ma and is marked by a prominent reflector, J. Beneath the northern Blake Plateau basin, only the hinge SDRs are observed, but they developed after J. We suggest that the inferred north-to-south age difference of SDR emplacement implies a heretofore unrecognized time-transgressive breakup of northwest Africa and North America during the early Middle Jurassic.
Sen, M. K., and P. L. Stoffa, Global Optimization Methods in Geophysical Inversion, Advances in Exploration Geophysics 4, Elsevier, 281 pp., 1995, 219 citations, #1075
Zhou, R., F. Tajima, and P. L. Stoffa, Application of genetic algorithms to constrain near-source velocity structure for the 1989 Sichuan earthquakes, Bull. Seismol. Soc. Amer., 85, 590-605, 1995, 19 citations, #1043
Zhou, R., F. Tajima, and P. L. Stoffa, Earthquake source parameter determination using genetic algorithms, Geophys. Res. Lett., 22, 517-520, 1995, 10 citations, #1072
Determination of earthquake source parameters is a non‐linear and multi‐dimensional problem in which trade‐off problems between the parameters can be significant and complicated. Accordingly global optimization is desirable to find the best fit parameters in a large model parameter space. We apply Genetic Algorithms (GA) to search for optimal source parameters in a multiparameter space. Results show that the GA method can rapidly explore a large multiparameter model space and find optimal solutions for the source parameters with a relatively small number of function evaluations. Results with the data of the 1988 Nepal‐India border earthquake indicate that GA determined source parameters better explain the observed waveforms.
Faria, E. L., and P. L. Stoffa, Traveltime computation in transversely isotropic media, Geophysics, 59, 272-281, 1994, 27 citations, doi:10:1190/1.1443539, #945
An approach for calculating first-arrival traveltimes in a transversely isotropic medium is developed and has the advantage of avoiding shadow zones while still being computationally fast. Also, it works with an arbitrary velocity grid that may have discontinuities. The method is based on Fermat's principle. The traveltime for each point in the grid is calculated several times using previously calculated traveltimes at surrounding grid points until the minimum time is found. Different ranges of propagation angle are covered in each traveltime calculation such that at the end of the process all propagation angles are covered. This guarantees that the first-arrival traveltime for a specific grid point is correctly calculated. The resulting algorithm is fully vectorizable. The method is robust and can accurately determine first-arrival traveltimes in heterogeneous media. Traveltimes are compared to finite-difference modeling of transversely isotropic media and are found to be in excellent agreement. An application to prestack migration is used to illustrate the usefulness of the method.
Faria, E. L., and P. L. Stoffa, Finite-difference modeling in transversely isotropic media, Geophysics, 59, 282-289, 1994, 29 citations, doi:10:1190/1.1443590, #983
We developed a modeling algorithm for transversely isotropic media that uses finite-difference operators in a staggered grid. Staggered grid schemes are more stable than the conventional finite-difference methods because the differences are actually based on half the grid spacing. This modeling algorithm uses the full elastic wave equation that makes possible the modeling of all kinds of waves propagating in transversely isotropic media. The spatial derivatives are represented by fourth-order, finite-difference operators while the time derivative is represented by a secondorder, finite-difference operator. The algorithm has no limitation on the acquisition geometry or on the heterogeneity of the media. The program is currently formulated to work in a 2-D transversely isotropic medium but can readily be extended to 3-D. Snapshots can be obtained at any time with no additional computational cost. A four-layer model is used to show the usefulness of the method. Horizontal and vertical component seismograms are modeled in transversely isotropic media and compared with seismograms modeled in the corresponding isotropic media.
Holbrook, W. S., E. C. Reiter, G. M. Purdy, D. S. Sawyer, P. L. Stoffa, J. A. Austin, J. Oh, and J. Makris, Deep structure of the U.S. Atlantic continental margin, offshore South Carolina, from coincident ocean bottom and multichannel seismic data, J. Geophys. Res., 99, 9155-9178, 1994, 78 citations, #971
We present the results of a combined multichannel seismic reflection (MCS) and wide-angle, ocean bottom seismic profile collected in 1988 across the Carolina Trough on the U.S. Atlantic continental margin. Inversion of vertical-incidence and wide-angle travel time data has produced a velocity model of the entire crust across the continent-ocean transition. The margin consists of three structural elements: (1) rifted continental crust, comprising 1â4 km of post-rift sedimentary rocks overlying a 30â34 km thick subsedimentary crust, (2) transitional crust, a 70- to 80-km-wide zone comprising up to 12 km of postrift sedimentary rocks overlying a 10- to 24-km-thick subsedimentary crust, and (3) oceanic crust, comprising 8 km of sedimentary rocks overlying an 8-km-thick crystalline crust. The boundary between rifted continental and transitional crust, marked by the Brunswick magnetic anomaly, represents an abrupt change in physical properties, with strong lateral increases in seismic velocity, density, and magnetic susceptibility. The transitional crust contains mid-crustal seaward-dipping reflections observed on the MCS section and has seismic velocities of 6.5â6.9 km/s in the midcrust and 7.2â7.5 km/s in the lower crust. Modeling of potential field data shows that transitional crust also produces the prominent, margin-parallel gravity anomaly and the Brunswick and East Coast magnetic anomalies. These observations support the interpretation that the transitional crust was formed by magmatism during continental breakup. The prodigious thickness (up to 24 km) of igneous material rivals that interpreted on continental margins of the North Atlantic (e.g., Hatton Bank and Vøring Plateau), which formed in the vicinity of the Iceland hotspot. These observations, when combined with other transects across the margin, confirm previous suggestions that the U.S. Atlantic margin is strongly volcanic and further imply that the magmatism was not the result of a long-lived mantle plume.
Sen, M. K., A. Datta-Gupta, P. L. Stoffa, L. Lake, and G. Pope, Stochastic reservoir modeling using simulated annealing and genetic algorithm, Proc., Ann. Meet., Soc. of Petrol. Engineers, Paper SPE 24754, 939-950, 1994, #941
Shipley, T. H., G. F. Moore, N. L. Bangs, J. C. Moore, and P. L. Stoffa, Seismically inferred dilatancy distribution, northern Barbados Ridge decollement: Implications for fluid migration and fault strength, Geology, 22, 411-414, 1994, 92 citations, doi:10.1130/0091-7613(1994)022<0411:SIDDNB>2.3.CO;2, #1025
A 5 x 25 km, three-dimensional seismic survey of the lower part of the northern Barbados Ridge accretionary prism creates a three-dimensional image of a major active decollement fault. The fault is usually a compound negative-polarity reflection modeled as a low-velocity, high-porosity zone less than ∼14 m thick. This thickness is significantly less than that defined by drilling of a >40 m zone of deformation at Ocean Drilling Program (ODP) Site 671B, located within the surveyed area. We infer that the seismically defined fault is a thin, high-porosity zone and is thus an undercompacted, high-fluid-pressure dilatant section. If these inferences are correct, then map-view variations in seismic-reflection waveform and amplitude illustrate complex patterns of fault-zone fluid content and fluid migration paths. The amplitude map suggests kilometre-wide channels of locally high porosity and thus focused fluid flow. These paths are only subparallel to the expected minimum head, as inferred from the shape of the overlying sediment wedge; other factors must modify fluid concentrations and ultimately migration. Several areas of positive-polarity fault reflections define square-kilometre-sized regions inferred to be lower porosity sections producing strong asperities in an otherwise weak fault. One, coincident with Site 671B, may explain the success of drilling through the fault here. All other holes drilled in the area were within the negative-polarity regions and were unsuccessful in penetrating through the entire fault zone, possibly because of instability associated with high fluid pressures and a weak fault. ODP Leg 156 planned for 1994 will test inferences related to fault permeability and fluid pressures.
Squires, L. J., P. L. Stoffa, and G. Cambois, Borehole transmission tomography for velocity plus statics, Geophysics, 59, 1028-1036, 1994, 9 citations, doi:10:1190/1.1443659, #1012
The accuracy of velocity tomograms reconstructed from borehole transmission traveltime data is highly sensitive to traveltime statics. We present a least-squares tomography algorithm that includes a traveltime static term. The algorithm solves for both the velocity field and the traveltime statics simultaneously. This enables us to separate traveltime signal from traveltime noise, reducing the tomographic velocity artifacts caused by the statics. The incorporation of a priori constraints on the poorly determined spectral components of the velocity field further improves accuracy by reducing velocity artifacts as a result of uneven ray coverage. Application of the algorithm to numerical crosswell data results in velocity and statics' estimates that are accurate to within 1 percent. Application of the algorithm to Exxon's Friendswood tomography data results in velocity and statics' estimates that correlate with independent data.
Wood, W. T., P. L. Stoffa, and T. H. Shipley, Quantitative detection of methane hydrate through high-resolution seismic velocity analysis, J. Geophys. Res., 99, 9681-9695, 1994, 56 citations, #950
A laterally extensive, high-resolution travel time velocity analysis and acoustic wave form, inversion were used to quantitatively determine methane hydrate content in deep water sediments of the Blake Ridge off the southeast U.S. coast. The interval acoustic velocity (Vp ) analyses were performed in the τ-p domain by interactively picking the τ-p trajectories of prominent reflections in each of 50 plane wave-decomposed common midpoint gathers. The reflections correspond to seismic stratigraphic boundaries so that lateral Vp changes due to lithology changes are mitigated, and Vp changes due to changing hydrate content are enhanced. Two separate interval Vp analyses were performed, one with thick (∼0.4 km) layers which yielded lower uncertainty but also lower resolution, and one with thinner layers (∼0.1 km), yielding higher resolution but slightly larger uncertainties. Results show no correlation between low-sediment reflectivity and Vp . However, in the areas exhibiting a bottom simulating reflector (BSR) a high Vp interval (∼2.0 km/s and 0.15 km thick) is seen immediately above the BSR. Where the BSR is strongest a 256-layer, least squares acoustic wave form inversion reveals the BSR to be caused by a Vp decrease from ∼2.0 to ∼1.5 km/s, with little or no change in density. The inversion also reveals a thin (0.025 km) layer of anomalously low Vp lying immediately below the BSR. Two models of methane hydrate distribution are tested, each indicating that the volume of methane hydrate in the intervals of elevated Vp is up to ∼25% of the total volume.
Cambois, G., and P. L. Stoffa, Surface-consistent phase decomposition in the log/Fourier domain, Geophysics, 58, 1099-1111, 1993, 2 citations, doi:10:1190/1.1443494, #926
In the log/Fourier domain, decomposing the amplitude spectra of seismic data into surface-consistent terms is a linear problem that can be solved, very efficiently, one frequency at a time. However, the nonunique definition of the complex logarithm makes it much more difficult to decompose the phase spectra. The instability of phase unwrapping has previously prevented any attempt to decompose phase spectra in the log/Fourier domain. We develop a fast and robust partial unwrapping algorithm, which makes it possible to efficiently decompose the phase spectra of normal moveout-corrected (NMO-) data into surface-consistent terms, in the log/Fourier domain. The dual recovery of amplitude and phase spectra yields a surface-consistent deconvolution technique where only the average reflectivity is assumed to be white, and only the average wavelet is required to be minimum-phase. Each individual deconvolution operator may be mixed-phase, depending on its estimated phase spectra. For example, surface-consistent time shifts and phase rotations, as well as any other surface-consistent phase effects, are included in the phase spectra of the surface-consistent deconvolution operators. Consequently, static shifts are estimated and removed without ever picking horizons or crosscorrelations.
Jervis, M., P. L. Stoffa, and M. K. Sen, 2-D migration velocity estimation using a genetic algorithm, Geophys. Res. Lett., 20, 1495-1498, 1993, 7 citations, #989
We address the problem of velocity estimation in heterogeneous media using a combination of nonlinear inversion and migration velocity analysis. In velocity estimation, the travel time information in seismic reflection data are nonlinearly related to the velocity perturbations in the subsurface. By taking a functional of seismic traces, the migrated data themselves, we define a misfit criterion which greatly reduces the oscillatory nature of the objective function. Migration is inherently a smoothing process; it collapses diffractions, focuses reflected energy and suppresses random noise. We use the lateral consistency of reflectors after migration as a measure of model misfit. If we compare one migrated shot record with an adjacent record, the misfit will be only slightly affected by large velocity variations even though reflectors may show large errors in depth positioning. The closer the spacing between shot records the smaller the variations in misfit; hence the smoother the objective function.
We search for global minima of the objective function thus defined, using a genetic algorithm (GA) and a linearized inversion scheme. We illustrate the techniques and results from the algorithms by applying them to a realistic scale synthetic data set. As expected, the success of a linearized scheme depends strongly on the starting model while GA does not depend on the choice of the initial population of models. GA is computationally intensive but our choice of spline parameterization reduces the number of model parameters significantly and makes the algorithm computationally tractable for real earth problems.
Sen, M. K., B. B. Bhattacharya, and P. L. Stoffa, Nonlinear inversion of resistivity sounding data, Geophysics, 58, 496-507, 1993, 41 citations, doi:10:1190/1.1443432, #899
The resistivity interpretation problem involves the estimation of resistivity as a function of depth from the apparent resistivity values measured in the field as a function of electrode separation. This is commonly done either by curve matching using master curves or by more formal linearized inversion methods. The problems with linearized inversion schemes are fairly well known; they require that the starting model be close to the true solution. In this paper, we report the results from the application of a nonlinear global optimization method known as simulated annealing (SA) in the direct interpretation of resistivity sounding data. This method does not require a good starting model but is computationally more expensive. We used the heat bath algorithm of simulated annealing in which the mean square error (difference between observed and synthetic data) is used as the energy function that we attempt to minimize. Samples are drawn from the Gibbs probability distribution while the control parameter the temperature is slowly lowered, finally resulting in models that are very close to the globally optimal solutions. This method is also described in the framework of Bayesian statistics in which the Gibbs distribution is identified as the a posteriori probability density function in model space. Computation of the true posterior distribution requires computation of the energy function at each point in model space. However, a fairly good estimate of the most significant portion(s) of the function can be obtained from simulated annealing run in a reasonable computation time. This can be achieved by making several repeat runs of SA, each time starting with a new random number seed so that the most significant portion of the model space is adequately sampled. Once the posterior density function is known, many measures of dispersion can be made. In particular, we compute a mean model and the a posteriori covariance matrix. We have applied this method successfully to synthetic and field data. The resulting correlation covariance matrices indicate how the model parameters affect one another and are very useful in relating geology to the resulting resisitivity values.
Cambois, G., and P. L. Stoffa, Surface-consistent deconvolution in the log/Fourier domain, Geophysics, 57, 823-840, 1992, 7 citations, doi:10:1190/1.1443296, #884
In the surface-consistent hypothesis, a seismic trace is the convolution of a source operator, a receiver operator, a reflectivity operator (representing the subsurface structure) and an offset-related operator. In the log/Fourier domain, convolutions become sums and the log of signal amplitude at a given frequency is the sum of source, receiver, structural, and offset-related terms. Recovering the amplitude of the reflectivity for a given frequency is then a linear problem (very similar to a surface-consistent static correction problem). However, this linear system is underconstrained. Thus, among the infinite number of possible solutions, a particular one must be selected. Studies with real data support the choice of a spatially band-limited solution. The surface-consistent operators can then be calculated very efficiently using an inverse Hessian method. Applications to real seismic data show improvement compared with previous techniques. Surface-consistent deconvolution is robust and fast in the log/Fourier domain. It allows the use of long operators, improves statics estimation, and removes the amplitude variations due to source or receiver coupling.
Sen, M. K., and P. L. Stoffa, Rapid sampling of model space using genetic algorithms: Examples from seismic waveform inversion, Geophys. J. Int., 108, 281-292, 1992, 82 citations, doi:10.1111/j.1365-246X.1992.tb00857.x, #871
The seismic waveform inversion problem is usually cast into the framework of Bayesian statistics in which prior information on the model parameters is combined with the data and physics of the forward problem to estimate the a posteriori probability density (PPD) in model space. The PPD is a function of an objective or fitness function computed from the observed and synthetic data. In general, the PPD or the fitness function is multimodal and its shape is unknown. Global optimization methods such as simulated annealing (SA) and genetic algorithms (GA) do not require that the shape of the fitness function be known. In this paper, we investigate GA to rapidly sample the most significant portion or portions of the PPD, when very little prior information is available. First, we use a simple three operator (selection, crossover and mutation) GA acting on a randomly chosen finite population of haploid binary coded models. We use plane wave transformed synthetic seismic data and a normalized cross-correlation function [E(m)] in the frequency domain as a fitness function. A moderate value of crossover probability, a low value of mutation probability, a high value of update probability and a proper population size are required to reach very close to the global maximum of the fitness function. Next, with an attempt to accelerate convergence we show that the concepts from simulated annealing can be used in stretching of the fitness function, i.e., we use exp [E(m)/T] rather than E(m) as the fitness function, where T is a control parameter analogous to temperature in simulated annealing. By a schemata analysis, we show that at low temperatures, schemata with above average fitness values are reproduced in large numbers causing a much more rapid convergence of the algorithm. A high value of temperature T assigns nearly equal selection probability to most of the schemata and thus retains diversity among the members of the population. Thus a GA with a step function type cooling schedule (very high temperature in the beginning followed by rapid cooling to a very low temperature) improves the performance dramatically: high values of the fitness function are obtained rapidly using only half as many models as would be required by a conventional GA. Similar performance could also be achieved by first using a high mutation probability and then decreasing the mutation probability to a very low value, while retaining the same low temperature throughout.
We also address the problem of 'genetic drift' which causes the finite GAs to converge to one peak or the other when the algorithm is applied to a highly multimodal fitness function with several peaks of nearly the same height. A parallel genetic algorithm based on the concept of 'punctuated equilibria' is implemented to circumvent the problem. We run several GAs each with a finite subpopulation in parallel and collect many good models from each one of these runs. These are then used to grasp the most significant portion(s) of the PPD in model space. We then compute the weighted mean model and use the derived good models to estimate uncertainty in the derived model parameters.
Sen, M. K., and P. L. Stoffa, Genetic inversion of AVO, Leading Edge, 11, 27-29, 1992, doi:10.1190/1.1436845, #887
Dan Hampson's article AVO inversion, theory and practice (June 1991 TLE) describes some recent developments, particularly linearized inversion techniques and Monte Carlo methods that can be applied to the inversion of AVO data to determine rock properties. He also described the problem of nonuniqueness associated with the inversion process and pointed out that we have to think about inversion as a process that gives many reasonable answers.
Shipley, T. H., K. D. McIntosh, E. A. Silver, and P. L. Stoffa, Three-dimensional seismic imaging of the Costa Rica accretionary prism: Structural diversity in a small volume of the lower slope, J. Geophys. Res., 97, 4439-4459, 1992, 75 citations, #901
Conventional two-dimensional seismic reflection investigations have been generally relied upon to provide images of large to medium scale structural features in accretionary prisms. We undertook a three-dimensional seismic reflection survey of a small part of a prism arcward of the Middle America Trench off Costa Rica to more correctly image structure and to use the improved structural information to examine the processes of accretion. This survey reveals small features, with dimensions of hundreds of meters, while also defining features thousands of meters in lateral extent, both of which were underappreciated in conventional two-dimensional data from the same area. We have imaged active off scraping at the trench and both duplexing and out-of-sequence faulting a few kilometers arcward of the trench. Fault spacing and reflector geometry vary dramatically over a space of several hundred meters. Some of these variations are related to visible changes in morphology of the underlying oceanic basement, but others are not so easily documented. Fault surface reflections define an architecture which may control gross fluid motion through the prism. This architecture is apparently formed by duplexing and out-of-sequence faulting and has been maintained by periodic motion on some of the out-of-sequence faults. The slope sediment apron records multiple phases of deformation. Abundant small offset reverse faults break the seafloor and indicate recent shortening of a broad region of the underlying prism. A primary result of this survey is appreciation of the structural diversity across a small width of an accretionary prism.
Squires, L. J., S. N. Blakeslee, and P. L. Stoffa, The effects of statics on tomographic velocity reconstructions, Geophysics, 57, 353-362, 1992, 16 citations, doi:10:1190/1.1443249, #921
Seismic first arrival times from crosshole, VSP, and reversed VSP (RVSP) experiments are collectively inverted by least-squares for the velocity distribution between two boreholes. The tomographic reconstruction exhibits a large lateral velocity contrast that is not supported by the surface reflection data from the same location. After examining the traveltime residuals from the three tomographic datasets separately, we conclude that the velocity contrast is due primarily to static delays in the RVSP first arrival times. When a first-order correction is made for the statics, tomographic inversion results in a velocity reconstruction that is more consistent with the surface reflection data. To isolate the velocity errors produced by the RVSP statics, we compute a residual tomogram by subtracting the statics adjusted tomogram from the original. The residual tomogram shows that the statics introduce errors not only in the region sampled by the RVSP rays, but they indirectly contaminate other regions of the tomogram as well. We reproduce this velocity error distribution as part of a model study designed to simulate the effects of statics on tomographic velocity reconstructions. Results indicate that traveltime errors on the order of 2 percent can result in tomographic velocity errors of up to 7 percent.
Stoffa, P. L., W. T. Wood, T. H. Shipley, G. F. Moore, E. Nishiyama, M. A. B. Botelho, A. Taira, H. Tokuyama, and K. Suyehiro, Deepwater high-resolution expanding spread and split spread seismic profiles in the Nankai Trough, J. Geophys. Res., 97, 1687-1713, 1992, 18 citations, #803
In deep water the source-receiver offsets that are required for accurate velocity determinations cannot be achieved with single-ship multichannel seismic methods. Two ships, one equipped with a multichannel receiving array and the other with a seismic source, have previously been employed to acquire common midpoint, expanding spread profiles, principally to determine deep crustal velocity structures. We extend this method to higher resolution in space and time to determine the velocities of sedimentary layers in deep water offshore Japan in the Nankai Trough. This high-resolution two-ship data acquisition method used a 13.1-L water gun source array; a 1.6-km, 96-channel receiving array with 0.0166-km active group; and shore-based navigation. Analysis of the data was performed in the τ-p domain by successive downward continuation of the plane wave data. Interactive velocity analysis methods for both one-dimensional and two-dimensional Earth models are described for both common source/receiver and common midpoint profiles. Results in one of two areas surveyed show a low-velocity zone below the subduction decollement which is consistent with models of low wedge taper, high pore fluid pressure, and reflection polarity reversal described previously by other researchers. The velocity profiles show the expected landward increase in velocity assumed to be due to lateral strain and porosity decrease, but the effect is small, only slightly greater than would be expected in an area of no lateral strain.
Stoffa, P. L., and M. K. Sen, Seismic waveform inversion using global optimization, J. Seismic Explor., 1, 9-27, 1992, #920
Diebold, J. B., and P. L. Stoffa, The traveltime equation, tau-p mapping, and inversion of common midpoint data, Slant-stack processing, edited by G. H. Gardner and L. Liu, Soc. Expl. Geophys., 1991, #2227
Loewenthal, D., and P. L. Stoffa, Synthetic acoustic seismograms by deverberating sources, J. Acoustical Soc. Amer., 90, 1101-1105, 1991, 4 citations, doi:10.1121/1.402299, #897
Using the dereverberated-source concept, the synthetic seismograms for stratified equal travel-time layering for source and receiver at the surface as well as for a buried receiver as in the VSP data are derived. Here the simpler acoustic representation of seismograms is treated. These results are extended to data of equal vertical delay time, , and horizontal ray parameter, p. Since only time domain operators are used, it is not necessary to use Fourier transforms to compute the seismic response for any depth point, or any ray parameter. In addition to extending Goupillaud's equal travel-time layer synthetic seismogram method to a buried receiver and to the case of non-normal incidence, the use of the dereverberated-source concept intuitively illustrates the nature of the seismic impulse response for a one-dimensional earth.
Moore, G. F., D. E. Karig, T. H. Shipley, A. Taira, P. L. Stoffa, and W. T. Wood, Structural framework of the ODP Leg 131 area, Nankai Trough, Proc. Ocean Drilling Prog., Init. Rept., 131, 15-20, 1991, #890
Oh, J., J. D. Phillips, J. A. Austin, and P. L. Stoffa, Deep-penetration seismic reflection images across the southeastern United States continental margin, in Continental Lithosphere: Deep Seismic Reflections, edited by H.-J. Durbaum, Amer. Geophys. Union Geodynamics Ser., 22, 225-240, 1991, #865
Sen, M. K., and P. L. Stoffa, Nonlinear one-dimensional seismic waveform inversion using simulated annealing, Geophysics, 56, 1624-1638, 1991, 161 citations, doi:10:1190/1.1442973, #821
The seismic inverse problem involves finding a model m that either minimizes the error energy between the data and theoretical seismograms or maximizes the cross-correlation between the synthetics and the observations. We are, however, faced with two problems: (1) the model space is very large, typically of the order of 5050; and, (2) the error energy function is multimodal. Existing calculus-based methods are local in scope and easily get trapped in local minima of the energy function. Other methods such as "simulated annealing" and "genetic algorithms" can be applied to such global optimization problems and they do not depend on the starting model. Both of these methods bear analogy to natural systems and are robust in nature. For example, simulated annealing is the analog to a physical process in which a solid in a "heat bath" is heated by increasing the temperature, followed by slow cooling until it reaches the global minimum energy state where it forms a crystal. To use simulated annealing efficiently for 1-D seismic waveform inversion, we require a modeling method that rapidly performs the forward modeling calculation and a cooling schedule that will enable us to find the global minimum of the energy function rapidly. With the advent of vector computers, the reflectivity method has proved successful and the time of the calculation can be reduced substantially if only plane-wave seismograms are required. Thus, the principal problem with simulated annealing is to find the critical temperature, i.e., the temperature at which crystallization occurs. By initiating the simulated annealing process with different starting temperatures for a fixed number of iterations with a very slow cooling, we noticed that by starting very near but just above the critical temperature, we reach very close to the global minimum energy state very rapidly. We have applied this technique successfully to band-limited synthetic data in the presence of random noise. In most cases we find that we are able to obtain very good solutions using only a few plane wave seismograms.
Stoffa, P. L., T. H. Shipley, W. Kessinger, D. F. Dean, R. Elde, E. A. Silver, D. L. Reed, and A. Aguilar, Three-dimensional seismic imaging of the Costa Rica accretionary prism: Field program and migration examples, J. Geophys. Res., 96, 21693-21712, 1991, 18 citations, #849
Seismic reflection techniques are the primary geophysical tool used to examine the structure of continental margins. On convergent plate margins, widely separated seismic reflection profiles often do not image complex structural features because out-of-plane reflections and diffractions obscure the seismic images. Yet the structural features are often critical to understanding the tectonic processes. Thus, the resulting geologic interpretations are based on images that may not accurately represent the true subsurface structure. To image the complex geologic structures of an active continental margin, a three-dimensional (3-D) seismic survey was conducted off of the Nicoya Peninsula of Costa Rica in April 1987. Over 60,000 shot points of 96-trace multichannel data were collected in a 9 Ã 22 km grid. This detailed survey was located over the accretionary prism and designed to investigate the active tectonic processes and evolution of this continental margin. In this paper we outline the data acquisition program, including the navigation processing critical to successful imaging, and the seismic processing methodology and then show several examples of images improved by 3-D surveying. The resulting 3-D seismic data were migrated both as independent lines and as a 3-D volume. Seismic images were produced as conventional depth and time sections and as cross-sectional slices in depth and time. We present comparisons of the migration results that show significant improvements of the subsurface image using 3-D techniques. To illustrate these improvements, we present images of two complex geologic areas: a mud volcano and reflections near the top of the accretionary prism, and the top of the subducting oceanic basement beneath the prism. Comparisons between the original stacked data, the data of the two-dimensional migration, and the results of the 3-D migration illustrate the value of 3-D techniques in studies of active margins. In addition to the advantages of the 3-D imaging which include better structural delineation and improved signal-to-noise ratio in the final image, we show that the ability to view the 3-D data volume in both section and plan view offers significant interpretational advantages.
Stoffa, P. L., and M. K. Sen, Nonlinear multiparameter optimization using genetic algorithms: Inversion of plane-wave seismograms, Geophysics, 56, 1794-1810, 1991, 155 citations, doi:10:1190/1.1442992, #855
Seismic waveform inversion is one of many geophysical problems which can be identified as a nonlinear multiparameter optimization problem. Methods based on local linearization fail if the starting model is too far from the true model. We have investigated the applicability of "Genetic Algorithms" (GA) to the inversion of plane-wave seismograms. Like simulated annealing, genetic algorithms use a random walk in model space and a transition probability rule to help guide their search. However, unlike a single simulated annealing run, the genetic algorithms search from a randomly chosen population of models (strings) and work with a binary coding of the model parameter set. Unlike a pure random search, such as in a "Monte Carlo" method, the search used in genetic algorithms is not directionless. Genetic algorithms essentially consist of three operations, selection, crossover, and mutation, which involve random number generation, string copies, and some partial string exchanges. The choice of the initial population, the probabilities of crossover and mutation are crucial for the practical implementation of the algorithm. We investigated the effects of these parameters in the inversion of plane-wave seismograms in which a normalized crosscorrelation function was used as the objective or fitness function (E). We also introduce the concept of "update" probability to control the influence of past generations. The combination of a low value of mutation probability (~0.01), a moderate value of the crossover probability (~0.6) and a high value of update probability (~0.9) are found to be optimal for the convergence of the algorithm. Further, we show that concepts from simulated annealing can be used effectively for the stretching of the fitness function which helps in the convergence of the algorithm. Thus, we propose to use exp (E/T) rather than E as the fitness function, where T (analogous to temperature in simulated annealing) is a properly chosen parameter which can change slowly with each generation. Also, by repeating the GA optimization procedure several times with different randomly chosen initial model populations, we derive "a very good subset" of models from the entire model space and calculate the a posteriori probability density (m) exp (E(m)/T). The (m) 's are then used to calculate a "mean" model, which is found to be close to the true model.
Stoffa, P. L., P. Buhl, J. B. Diebold, and F. Wenzel, Direct mapping of seismic data to the domain of intercept time and ray parameter - A plane-wave decomposition, Slant-stack processing, edited by G. H. Gardner and L. Liu, Soc. Expl. Geophys., 1991, #2228
Austin, J. A., P. L. Stoffa, J. D. Phillips, J. Oh, D. S. Sawyer, G. M. Purdy, E. C. Reiter, and J. Makris, Crustal structure of the southeast Georgia embayment-Carolina trough: Preliminary results of a composite seismic image of a continental suture(?) and a volcanic passive margin, Geology, 18, 1023-1027, 1990, 35 citations, doi:10.1130/0091-7613(1990)018<1023:CSOTSG>2.3.CO;2, #819
New deep-penetration multichannel seismic reflection data, combined with refraction results and magnetics modeling, support a hypothesis that the Carolina trough is a Mesozoic volcanic passive margin exhibiting a seaward-dipping wedge and associated underplating. The structure of Carolina platform continental crust is consistent with the late Paleozoic continental collision that produced the Appalachians, but imbrication has had no obvious effect on shallower structures produced by Mesozoic extension and volcanism. The origin of prominent magnetic anomalies crossing the Southeast Georgia embayment can be explained by processes attending Mesozoic separation of Africa and North America, and is not related to a Paleozoic continental suture, as previously postulated.
Moore, G. F., T. H. Shipley, P. L. Stoffa, D. E. Karig, A. Taira, S. Kuramoto, H. Tokuyama, and K. Suyehiro, Structure of the Nankai Trough accretionary zone from multichannel seismic reflection data, J. Geophys. Res., 95, 8753-8765, 1990, 147 citations, #846
New multichannel seismic reflection data collected over the Nankai Trough image the accretionary complex in two areas: the International Program of Ocean Drilling leg 87 transect area (western area) and the region of upcoming Ocean Drilling Program leg 131 (eastern area). The incoming Shikoku Basin sedimentary section consists of hemipelagic clays and thin terrigenous turbidites. The basin section is overlain by a trench wedge that is 12â16 km wide and 350â750 m thick at the thrust front. Accretionary deformation begins in a protothrust zone that is characterized by thickening and seaward tilting of the trench wedge. The zone in the western area is 4.5 km wide and is characterized by âkinkâ folds; the zone in the eastern area is only 2.5 km wide and does not exhibit such folds. The frontal thrusts in each area are imaged as fault plane reflections and ramp upward from within the basin hemipelagic section. The overthrusting sediments form fault-bend folds over these ramps. Thrust spacing at the toe of the slope is 1.5â2.5 km. The second thrust cuts up from an inferred décollement within the Shikoku Basin sedimentary section. In the eastern area, a reflection marking the top of the basin pelagic sediment section changes from normal to reversed polarity about 6.3 km seaward of the thrust front and underlies the entire protothrust zone. This reflector continues with reversed polarity under the accretionary complex and is at the level of the basal décollement. The underlying basin pelagic section is apparently thrust undisturbed beneath the accretionary prism. The reversal of polarity indicates a change in reflection coefficient that is due to a combination of decreasing seismic velocity and density across the interface. This decrease in velocity and density may indicate that the décollement is a zone of high porosity due to fluid expulsion from deeper within the accretionary prism. The reflections from the first and second thrusts are also reversed polarity, possibly indicating that they also are pathways of fluid expulsion. The critical wedge taper of the western area is greater than that of the eastern area, an observation that is consistent with the existence of an overpressured décollement in the eastern area.
Shipley, T. H., P. L. Stoffa, and D. F. Dean, Underthrust sediments, fluid migration paths, and mud volcanoes associated with the accretionary wedge off Costa Rica: Middle American Trench, J. Geophys. Res., 95, 8743-8752, 1990, 75 citations, #925
Styles of deformation and tectonic responses resulting from the convergence of oceanic and continental plates are strongly dependent on fluids in the sediments. We estimate the volume of fluids and sediment underthrust beneath the toe of an accretionary wedge and describe evidence of fluid migration observed in seismic data in the Middle America Trench near Costa Rica. Reduction in normal incident travel time of the underthrust oceanic plate section arcward of the deformation front may be related to fluid expulsion. Within 4 km of the deformation front, where overburden is <800 m, more than half of the total water in the pore spaces may have been expelled. The rate of fluid transfer from the underthrust section is highest in the zone 0ââ¬â4 km arcward of the deformation front. The clarity of the reflection associated with the décollement is exceptional and shows a phase reversal, relative to the seafloor. This high-amplitude phase-reversed reflection ends about 4 km arcward of the deformation front, suggesting a rapid reduction in the density and/or velocity contrast between the dewatering underthrust section and overlying offscraped sediments at this position. Fluid migration within the wedge produces mud volcanoes or ridges in the midslope region having relief of about 50 m and extent of perhaps 500 m Ãâ 1000 m. A reflection traced beneath a seafloor mud volcano crosses through the accretionary wedge/slope cover boundary and may mark a fluid conduit. Locally, high-amplitude reflections associated with the base of slope cover are coincident with closure of structure and indicate local traps to upward organized or diffusive fluid flow. Some fluids appear to migrate along pathways marked by reflections which extend deep within the accretionary wedge but not always extending into the underthrust section. Near the trench where fluid expulsion rates are highest, water may exit to the seafloor seaward of the downslope pinch-out of the slope section.
Stoffa, P. L., J. T. Fokkema, R. De Luna Freire, and W. Kessinger, Split-step Fourier migration, Geophysics, 55, 410-421, 1990, 169 citations, doi:10:1190/1.1442850, #1104
The split-step Fourier method is developed and applied to migrating stacked seismic data in two and three dimensions. This migration method, which is implemented in both the frequency-wavenumber and frequency-space domains, takes into account laterally varying velocity by defining a reference slowness (reciprocal of velocity) as the mean slowness in the migration interval and a perturbation term that is spatially varying. The mean slowness defines a reference vertical wavenumber which is used in the frequency-wavenumber domain to downward continue the data across a depth interval as in constant-velocity phase-shift migration. The perturbation term is used to define a "source" contribution that is taken into account by the application of a second phase shift in the frequency-space domain. Since the method does not include the effects of second and higher order spatial derivatives of the slowness field, the method theoretically is accurate only when there are no rapid lateral slowness variations combined with steep angles of propagation. However, synthetic and real examples indicate that good results are obtained for realistic geologic structures.
Stoffa, P. L., Analysis and processing of wide-angle reflection and refraction seismic data in the τ-P domain, Adv. Geophys. Data Processing, 2, 81-117, 1989, #692
Stoffa, P. L., (Editor), τ-p: A Plane Wave Approach to the Analysis of Seismic Data, Modern Approaches in Geophysics 4, Kluwer, 176 pp., 1989, #924
Diebold, J. B., P. L. Stoffa, and the LASE Study Group, A large aperture seismic experiment in the Baltimore Canyon Trough, in Geology of North America, Vol. I-2: Atlantic Continental Margin, edited by R. E. Sheridan, and J. S. Grow, Geol. Soc. Amer., 388-398, 1988, #773
Johansen, B., O. Eldholm, M. Talwani, P. L. Stoffa, and P. Buhl, Expanding spread profile at the northern Jan Mayen Ridge, Polar Res., 6, 95-104, 1988, #771
Stoffa, P. L., Acquisition and analysis of wide-angle reflection and refraction seismic data, in Petroleum Resources of China and Related Subjects, edited by H. C. Wagner, L. C. Wagner, F. F. H. Wang, and F. L. Wong, Circum-Pacific Council for Energy & Mineral Resources Earth Sci. Series, 10, 695-718, 1988, #790
Loewenthal, D., P. L. Stoffa, and E. L. Faria, Suppressing the unwanted reflections of the full wave equation, Geophysics, 52, 1007-1012, 1987, 19 citations, doi:10:1190/1.1442352, #923
In many instances in exploration geophysics we are interested in the so-called one-way wave equation. This equation allows the wave fields to propagate in the positive depth direction, but not in the reverse (âZ) direction. Some modeling and migration methods, such as the f-k method (Stolt, 1978) and the phase-shift method (Gazdag, 1978), produce in a natural way the one-way wave equation. The main advantage of the one-way wave equation is that it does not give rise to multiples or interlayer reverberations and enables the observation of primary events only.
Stoffa, P. L., Analysis of seismic data in the τ-p domain, Revista Brasileira de Geofisica, 4 (2), 31-43, 1986, #772
Gamboa, L. A. P., M. Truchan, and P. L. Stoffa, Middle and Upper Jurassic depositional environments at outer shelf and slope of Baltimore Canyon Trough, AAPG Bull., 69, 610-621, 1985, 8 citations, #624
New CDP data acquired in the Baltimore Canyon Trough during project LASE made it possible to map a continuous Jurassic sedimentary sequence from the continental margin to the abyssal plain without interruption by basement structures. Intense carbonate sedimentation is inferred at the outer shelf during the Middle and Late Jurassic. Carbonate sedimentation probably started during the Middle Jurassic with a platform that prograded seaward with the development of ramps. By the Late Jurassic, a major reef complex had developed at the outer continental shelf. The onset of reef growth can be tentatively dated as 138 Ma by using the J1 reflector dated by the Deep Sea Drilling Project. A well-developed reef-talus deposit can be identified overlying the interface that generates th J1 reflector. A detailed analysis of semblance-derived interval velocities in the reef-talus sequence indicates a compressional velocity of 4.3-4.5 km/sec (14,100-14,800 ft/sec) for that interval, which was part of a major barrier reef along the United States eastern margin. After the reef formed, the deep oceanic basin was mostly starved from shelf-derived sediments until the reef died and was buried by clastic sediments. By correlation of our seismic data and COST well information, we infer that in the Baltimore Canyon Trough this reef had terminated by about the end of the Jurassic Period.
McCowan, D. W., P. L. Stoffa, and J. B. Diebold, Fan filters for data with variable spatial sampling, IEEE Trans. Acoustics, Speech and Signal Processing, 32, 1154-1159, 1984, 1 citation, #691
We develop a version of the velocity or "fan filter" which is applicable to seismic data with variable spatial sampling. Our formulation of the fan filter involves convolving a multichannel difference operator with the Hilbert transform of the data. We present results of a pre-stack application of this filter to wide aperture variable spaced CDP data. The results show that, in our example, the fanpass filter is able to attenuate multiples because of their residual moveout after NMO.
Mutter, J. C., M. Talwani, and P. L. Stoffa, Evidence for a thick oceanic crust adjacent to the Norwegian margin, J. Geophys. Res., 89, 483-502, 1984, 97 citations, #585
The oceanic crust created during this first few million years of accretion in the Norwegian-Greenland Sea lies at an unusually shallow depth for its age, has a smooth upper surface, and in many places the results of multichannel seismic reflection profiling reveal that its upper layers comprise a remarkable sequence of arcuate, seaward-dipping reflectors. These have been attributed to lava flows generated during a brief period of subaerial seafloor spreading. We describe the results of inversions of digitally recorded sonobuoy measurements and two-ship expanded spread profiles collected over the oceanic crust adjacent to the Norwegian passive margin. We find that the crust of the deep Lofoten Basin is indistinguishable from normal oceanic crust in thickness and structure. Closer to the margin we observe up to a four times expansion in thickness of layers with velocities equal to those of oceanic layer 2, while the layer 3 region retains approximately the same thickness. The area over which the seaward-dipping reflectors can be observed on reflection profiles corresponds to the region of greatest expansion in âLayer 2â thickness. In the very oldest crust immediately adjacent to an escarpment that probably marks the continent-ocean boundary, we see evidence for a low velocity zone overlying an indistinct reflector that may mark the dyke-lava interface in the thick crust. Comparing the structure of the thick crust to that of eastern Iceland, we find a strong resemblance, especially in the expansion in thickness of material with layer 2 velocities. These results support the suggestion that during the earliest stages of spreading extrusive volcanism at the ridge crest was unusually voluminous, building a thick pile of lavas erupted from a subaerial spreading center.