Dr. Mrinal Sen
Publications
Ray-Born inversion for fracture parameters
Liu, Y., and M. K. Sen, Scalar wave equation modeling with time-space domain dispersion-relation based staggered-grid finite-difference schemes, Bull. Seismol. Soc. Amer., 101, 141-159, 2011, 1 citation, doi:10.1785/0120100041, #2319 
The staggered-grid finite-difference (SFD) method is widely used in numerical modeling of wave equations. Conventional SFD stencils for spatial derivatives are usually designed in the space domain. However, when they are used to solve wave equations, it becomes difficult to satisfy the dispersion relations exactly. Liu and Sen (2009c) proposed a new SFD scheme for one-dimensional (1D) scalar wave equation based on the timeâââ‰â¬Åspace domain dispersion relation and plane wave theory, which is made to satisfy the exact dispersion relation. This new SFD scheme has greater accuracy and better stability than a conventional scheme under the same discretizations. In this paper, we develop this new SFD scheme further for numerical solution of 2D and 3D scalar wave equations. We demonstrate that the modeling accuracy is second order when the conventional 2M-th-order space-domain SFD and the second order time-domain finite-difference stencils are directly used to solve the scalar wave equation. However, under the same discretization, our 1D scheme can reach 2M-th-order accuracy and is always stable; 2D and 3D schemes can reach 2M-th-order accuracy along 8 and 48 directions, respectively, and have better stability. The advantages of the new schemes are also demonstrated with dispersion analysis, stability analysis, and numerical modeling.
Liu, Y., and M. K. Sen, Finite-difference modeling with adaptive variable-length spatial operators, Geophysics, 76, T79-T89, 2011, doi:10.1190/1.3587223, #2434
Most finite-difference simulation algorithms use fixed-length spatial operators to compute spatial derivatives. The choice of length is dictated by computing cost, stability, and dispersion criteria that are satisfied globally. We propose finite-difference schemes with adaptive variable-length spatial operators to decrease computing costs significantly without reducing accuracy. These schemes adopt long operators in regions of low velocity and short operators in regions of high velocity. Two methods automatically determine variable operator lengths. Dispersion analysis, along with 1D and 2D modeling, demonstrates the validity and efficiency of our schemes. In addition, a hybrid absorbing boundary condition helps reduce unwanted reflections from model boundaries. Our scheme is more efficient than those based on variable-grid methods for modeling, migration, and inversion of models with complex velocity structures because the latter require local grid refinement, which usually increases memory requirements and computing costs.
Liu, Y., and M. K. Sen, 3D acoustic wave equation modeling with time-space domain dispersion relation based staggered-grid finite-difference schemes and hybrid absorbing boundary conditions, Exploration Geophys., (in press), 2011, #2437
Liu, Y., and M. K. Sen, A hybrid absorbing boundary condition for elastic staggered grid modeling, Geophys. Prospecting, (in press), 2011, #2438
Sen, M. K., Seismic, Reflectivity method, Encyclopedia of Solid Earth Geophysics, Springer, edited by H. K. Gupta, 1269-1279, 2011, #2325
Sen, M. K., Inverse theory, Global optimization, Encyclopedia of Solid Earth Geophysics, Springer, edited by H. K. Gupta, 625-632, 2011, #2326
Sil, S., M. K. Sen, and B. Gurevich, Analysis of fluid substitution in a porous and fractured medium, Geophysics, 76, WA157-166, 2011, doi:10.1190/1.3564954, #2436
To improve quantitative interpretation of seismic data, we analyze the effect of fluid substitution in a porous and fractured medium on elastic properties and reflection coefficients. This analysis uses closed-form expressions suitable for fluid substitution in transversely isotropic media with a horizontal symmetry axis (HTI). For the HTI medium, the effect of changing porosity and water saturation on (1) P-wave moduli, (2) horizontal and vertical velocities, (3) anisotropic parameters, and (4) reflection coefficients are examined. The effects of fracture density on these four parameters are also studied. For the model used in this study, a 35% increase in porosity lowers the value of P-wave moduli by maximum of 45%. Consistent with the reduction in P-wave moduli, P-wave velocities also decrease by maximum of 17% with a similar increment in porosity. The reduction is always larger for the horizontal P-wave modulus than for the vertical one and is nearly independent of fracture density. The magnitude of the anisotropic parameters of the fractured medium also changes with increased porosity depending on the changes in the value of P-wave moduli. The reflection coefficients at an interface of the fractured medium with an isotropic medium change in accordance with the above observations and lead to an increase in anisotropic amplitude variation with offset (AVO) gradient with porosity. Additionally, we observe a maximum increase in P-wave modulus and velocity by 30% and 8%, respectively, with a 100% increase in water saturation. Water saturation also changes the anisotropic parameters and reflection coefficients. Increase in water saturation considerably increases the magnitude of the anisotropic AVO gradient irrespective of fracture density. From this study, we conclude that porosity and water saturation have a significant impact on the four studied parameters and the impacts are seismically detectable. ©2011 Society of Exploration Geophysicists
Tao, L., T. Hu, M. K. Sen, J. Yang, R. Wang, J. Wei, J. Wei, and S. Wang, A hybrid scheme for seismic modeling based on Galerkin method, Geophysical J. Int., 186, 1165-1178, 2011, doi:10.1111/j.1365-246X.2011.05094.x, #2435
Purely numerical methods for modelling seismic wave propagation are now fairly popular. In this paper, we report on the development of a hybrid finite element-finite difference method. Our method first employs semi-discretization of the finite element method in a part of the spatial domain (in the z-direction for 2-D situation) to obtain a wave equation in weak form, and then uses the finite difference method to solve it. This offers significant advantages in applying non-uniform grids with improved accuracy. To improve the computational efficiency, we introduce the spectral element method to replace the traditional finite element method, which leads to a diagonal mass matrix with sufficient accuracy. After that, we carry out detailed analyses of the dispersion behaviours for both uniform and non-uniform cases; the dispersion curves demonstrate high precision of our method in numerical modelling, especially in dealing with non-uniform grids. Some examples are presented to demonstrate performance of this method and confirm the analytic results.
Bansal, R., and M. K. Sen, Ray-Born inversion for fracture parameters, Geophys. J. Int., 180, 1274-1288, 2010, doi:10.1111/j.1365-246X.2009.04484.x, #2263
The purpose of this study is to investigate the possibility of inverting for the fracture parameters using the ray-Born approximation. The fracture parameters normal and tangential weaknesses ΔN and ΔT, derived using a low-frequency effective medium theory, are indicators of the crack density and fracture infill (liquid-filled or dry). In our implementation of ray-Born inversion, ΔN and ΔT are unknown model parameters. We treat the background medium as isotropic, and its elastic properties are assumed to be already determined using other conventional methods prior to inversion. In our current implementation, we assume that the fractures are vertical and the orientation of the symmetry axis of the effective anisotropic medium or the fracture normal is known.
We formulate the inverse problem in a least-squares sense and apply the conjugate-gradient method to minimize the error function. To overcome the problem of non-uniqueness in the inversion and to ensure that the inverted model parameters are physically meaningful and unique, we added an analytical function of model parameters to the error function. We show, using a numerical example, that the fracture parameters cannot be reliably obtained from single-component seismic data (vertical or horizontal). Simultaneous inversion of both components is needed for more robust estimates of ΔN and ΔT. The inversion scheme was implemented in the frequency domain, which allowed us to experiment with the number of frequencies used in the inversion. We found that a very small number of frequency samples covering the spectrum of the source wavelet could produce reasonable estimates of ΔN and ΔT. We tested our method on a synthetic data set generated using ray-Born forward modelling. For the experiment shown in the paper, the least-squares inversion scheme produced encouraging results. We also applied our method to a synthetic data set generated using a 3-D anisotropic finite-difference code. The results indicated that our method could predict the parameter trends, but with some inaccuracies in the values of ΔN and ΔT.
De Basabe, J. D., and M. K. Sen, Stability of the high-order finite elements for acoustic or elastic wave propagation with high-order time stepping, Geophys. J. Int., 181, 577-590, 2010, 3 citations, doi:10.1111/j.1365-246X.2010.04536.x, #2264
We investigate the stability of some high-order finite element methods, namely the spectral element method and the interior-penalty discontinuous Galerkin method (IP-DGM), for acoustic or elastic wave propagation that have become increasingly popular in the recent past. We consider the Lax-Wendroff method (LWM) for time stepping and show that it allows for a larger time step than the classical leap-frog finite difference method, with higher-order accuracy. In particular the fourth-order LWM allows for a time step 73 per cent larger than that of the leap-frog method; the computational cost is approximately double per time step, but the larger time step partially compensates for this additional cost. Necessary, but not sufficient, stability conditions are given for the mentioned methods for orders up to 10 in space and time. The stability conditions for IP-DGM are approximately 20 and 60 per cent more restrictive than those for SEM in the acoustic and elastic cases, respectively.
Jackson, C. S., M. K. Sen, P. L. Stoffa, and G. Huerta, Data directed importance sampling for climate model uncertainty estimation, in Advanced Computational Infrastructures for Parallel/Distributed Adaptive Applications, edited by M. Parashar, X. Li, and S. Chandra, editors, Wiley Publishers, 2010, #1912
Liu, Y., and M. K. Sen, Acoustic VTI modeling with a time-space domain dispersion-relation-based finite-difference scheme, Geophysics, 75, A11-A17, 2010, 3 citations, doi:10.1190/1.3374477, #2320
The finite-difference (FD) method is widely used in numerical modeling of wave equations. Conventional FD stencils for space derivatives are usually designed in the space domain. When they are used to solve wave equations, it is difficult to satisfy the dispersion relations exactly. We have designed a spatial FD stencil based on a time-space domain dispersion relation to simulate wave propagation in an acoustic vertically transversely isotropic (VTI) medium. Two-dimensional dispersion analysis and numerical modeling demonstrate that this stencil has greater precision than one used in a conventional FD, when the same number of grid points is used in the calculation of the spatial derivatives. Thus, the spatial FD stencil based on time-space domain dispersion relation can be used to replace a conventional one such that we can achieve greater accuracy with almost no increase in computational cost.
Liu, Y., and M. K. Sen, A hybrid scheme for absorbing edge reflections in numerical modeling of wave propagation, Geophysics, 75, A1-A6, 2010, 4 citations, doi:10.1190/1.3295447, #2324
We propose an efficient scheme to absorb reflections from the model boundaries in numerical solutions of wave equations. This scheme divides the computational domain into boundary, transition, and inner areas. The wavefields within the inner and boundary areas are computed by the wave equation and the one-way wave equation, respectively. The wavefields within the transition area are determined by a weighted combination of the wavefields computed by the wave equation and the one-way wave equation to obtain a smooth variation from the inner area to the boundary via the transition zone. The results from our finite-difference numerical modeling tests of the 2D acoustic wave equation show that the absorption enforced by this scheme gradually increases with increasing width of the transition area. We obtain equally good performance using pseudospectral and finite-element modeling with the same scheme. Our numerical experiments demonstrate that use of 10 grid points for absorbing edge reflections attains nearly perfect absorption.
Pecher, I. A., B. Milkereit, A. Sakai, M. K. Sen, N. L. Bangs, and J.-W. Huang, Vertical seismic profiles through gas-hydrate-bearing sediments, in Geophysical Characterization of Gas Hydrates, edited by M. Riedel, E. C. Wiloughby, and S. Chopra, SEG, 2010, doi:10.1190/1.9781560802197.ch8, #2291
Porsani, M. J., P. L. Stoffa, M. K. Sen, and R. K. Seif, Partitioned least squares operator for large-scale geophysical inversion, Geophysics, 75, R121-R128, 2010, 3 citations, doi:10.1190/1.3509780, #2313
Least-squares (LS) problems are encountered in many geophysical estimation and data analysis problems where a large number of observations (data) are combined to determine a model (some aspect of the earth structure). Examples of least squares in seismic exploration include several data processing algorithms, theoretically accurate LS migration, inversion for reservoir parameters, and background velocity estimation. A frequently encountered problem is that the volume of data in 3D is so large that the matrices required for the LS solution cannot be stored within the memory of a single computer. A new technique is described for parallel computation of the LS operator that is based on a partitioned-matrix algorithm. The classical LS method for solution of block-Toeplitz systems of normal equation (NE) to the general case of block-Hermitian and non-Toeplitz systems of NE. is generalized. Specifically, a solution of a block-Hermitian system of NE is shown that may be obtained recursively by linearly combining the solutions of lesser order that are related to the forward and backward subsystems of equations. This results in an efficient parallel algorithm in which each partitioned system can be evaluated independently. The application of the algorithm to the problem of 3D plane wave transformation is demonstrated.
Sil, S., R. P. Srivastava, and M. K. Sen, Observation of shear-wave splitting in the multi-component node data from Atlantis field, Gulf of Mexico, Geophys. Prospecting, 58, 953-964, 2010, 1 citation, doi:10.1111/J.1365-2478.2010.00871.x, #2323
In 2005, a multicomponent ocean bottom node data set was collected by BP and BHP Billiton in the Atlantis field in the Gulf of Mexico. Our results are based on data from a few sparse nodes with millions of shots that were analysed as common receiver azimuthal gathers. A first-order look at P-wave arrivals on a common receiver gather at a constant offset reveals variation of P-wave arrival time as a function of azimuth indicating the presence of azimuthal anisotropy at the top few layers. This prompted us to investigate shear arrivals on the horizontal component data. After preliminary processing, including a static correction, the data were optimally rotated to radial (R) and transverse (T) components. The R component shows azimuthal variation of traveltime indicating variation of velocity with azimuth; the corresponding T component shows azimuthal variation of amplitude and phase (polarity reversal). The observed shear-wave (S-wave) splitting, previously observed azimuthal P-wave velocity variation and azimuthal P-wave amplitude variation, all indicate the occurrence of anisotropy in the shallow (just below the seafloor) subsea sediment in the area. From the radial component azimuthal gather, we analysed the PP- and PS-wave amplitude variation for the first few layers and determined corresponding anisotropy parameter and VP/VS values. Since fracture at this depth is not likely to occur, we attribute the observed azimuthal anisotropy to the presence of microcracks and grain boundary orientation due to stress. The evidence of anisotropy is ubiquitous in this data set and thus it argues strongly in favour of considering anisotropy in depth imaging for obtaining realistic subsurface images, at the least.
Smith, D., M. K. Sen, and R. J. Ferguson, Data regularization and datuming by conjugate gradients, J. Seismic Explor., 19, 321-347, 2010, #2321
Srivastava, R. P., and M. K. Sen, Stochastic inversion of prestacked seismic data using fractal-based initial models, Geophysics, 75, R47-R59, 2010, 1 citation, doi:10.1190/1.3379322, #2322
In general, inversion algorithms rely on good starting models to produce realistic earth models. A new method, based on a fractional Gaussian distribution derived from the statistical parameters of available well logs to generate realistic initial models, uses fractal theory to generate these models. When such fractal-based initial models estimate P- and S-impedance profiles in a prestack stochastic inversion of seismic angle gathers, very fast simulated annealing ââ¬â a global optimization method ââ¬â finds the minimum of an objective function that minimizes data misfit and honors the statistics derived from well logs. The new stochastic inversion method addresses frequencies missing because of band limitation of the wavelet; it combines the low- and high-frequency variation from well logs with seismic data. This method has been implemented successfully using real prestack seismic data, and results have been compared with deterministic inversion. Models derived by a deterministic inversion are devoid of high-frequency variations in the well log; however, models derived by stochastic inversion reveal high-frequency variations that are consistent with seismic and well-log data.
Stoffa, P. L., M. K. Sen, R. K. Seifoullaev, and R. C. Pestana, Plane wave seismic data: Parallel and adaptive strategies for velocity analysis and imaging, in Advanced Computational Infrastructures for Parallel/Distributed Adaptive Applications, edited by M. Parashar, X. Li, and S. Chandra, Wiley Publishers, 2010, #1928
De Basabe, J. D., and M. K. Sen, New developments in the finite-element method for seismic propagation simulation, Leading Edge, 28, 562-567, 2009, doi:10.1190/1.3124931, #2092
Several versions of the finite-element method (FEM) have been proposed for seismic modeling since it was first applied in this field in the late 1960s and early 1970s. Here we bring forward some recently developed FEMs suitable for seismic modeling due to their high-order accuracy. These methods have been successfully applied to wave-propagation problems and improved accuracy, yielding an improved performance compared to classic FEM. ©2009 Society of Exploration Geophysicists
De Basabe, J. D., and M. K. Sen, Comment "Dispersion analysis of spectral element methods for elastic wave propagation," by G. Seriani and S. P. Oliveira, Wave Motion, 46, 92-93, 2009, 1 citation, doi:10.1016/j.wavemoti.2008.08.005, #2100
Dutta, U., M. K. Sen, N. N. Biswas, and Z. Yang, Investigation of shallow sedimentary structure of the Anchorage basin, Alaska, using simulated annealing inversion of site response, Bull. Seismol. Soc. Amer., 99, 326-339, 2009, 1 citation, doi:10.1785/0120070250, #2098
This study deals with shallow sedimentary structure of the Anchorage basin in Alaska. For this purpose, inversion of site response [SR(f)] data in the frequency range 0.5â11.0 Hz from various sites of the basin has been performed using the simulated annealing method to compute subsurface layer thickness, shear-wave velocity (β), density, and shear-wave quality factor. The one-dimensional (1D) models for the aforementioned parameters were obtained with preset bounds on the basis of available geological information such that the L-2 norm error between the observed and computed site response attained a global minimum. Next, the spatial distribution of the important parameter β was obtained by interpolating values yielded by the 1D models. The results indicate the presence of three distinct velocity zones as the source of spatial variation of SR(f) in the Anchorage basin. In the uppermost part of the basin, the β values of fine-grain Quaternary sediments mainly lie in the range of 180â500 m/sec with thickness varying from 15 to 50 m. This formation overlies relatively thick (80â200 m) coarse-grain Quaternary sediments with β values in the range of 600â900 m/sec. These two Quaternary units are, in turn, overlain on Tertiary sediments with β>1000 m/sec located at depths of 100 and 250 m, respectively, in the central and western side along the Knik Arm parts of the basin. The important implication of the result is that the sources of spatial variation of SR(f) in the Anchorage basin for the frequency band 0.5â11 Hz, besides in the uppermost 30 m, are found to be deeper than this depth. Thus, use of commonly considered geological formations in the depth intervals from 0 to 30 m for the ground-motion interpretation will likely yield erroneous results in the Anchorage basin.
Hong, T. C., and M. K. Sen, A new MCMC algorithm for seismic waveform inversion and corresponding uncertainty analysis, Geophys. J. Int., 177, 14-32, 2009, 3 citations, doi:10.1111/j.1365-246X.2008.04052.x, #2096
It is superior to formulate an inverse problem in a Bayesian framework and fully solve it by stochastically constructing the posterior probability density (PPD) distribution using Markov chain Monte Carlo (MCMC) algorithms. The estimated PPD can also be used to compute several measures of dispersion in the model space. However, for realistic application, MCMC methods can be computationally expensive and may lead to inaccurate PPD estimation as well as uncertainty analysis due to the strong non-linearity and high dimensionality. In this paper, to address the fundamental issues of efficiency and accuracy in parameter estimation and PPD sampling, we incorporate some new developments into a standard genetic algorithm (GA) to design more powerful algorithms for the practical geophysical inverse problems such as a non-linear pre-stack seismic waveform inversion.
First, a multiscale real-coded hybrid GA is developed to facilitate exploitation of the model space for optimal parameters at a fine scale. It is demonstrated that, by using real-coding and especially multiscaling to trade information between the model vectors defined at different resolutions, we attain a substantial speed-up in computation and obtain accurate parameter estimations. This new optimization method is further adapted to a new multiscale GA based MCMC method, in which multiple MCMC chains defined at different scales are run simultaneously in parallel. To gain the benefits of both the faster convergence of coarse scales and the greater detail of fine scales, realizations of chains at different scales are combined for intelligent proposals that facilitate exploration of the model space at the fine scale. In this study, the new MCMC is justified using an analytical example and its performance on PPD estimation, and uncertainty quantification is demonstrated using a non-linear seismic inverse problem. We find that incorporation of multiscaling in the Bayesian approach shows a great promise in solving inverse problems that involve computationally intensive forward modelling. We also find that the multiscaling is particularly attractive in addressing the model parametrization issue since the job of 'model selection' could be done based on the results of 'parameter estimation' at multiple scales, which circumvents separately solving a model selection problem in a Bayesian framework.
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.
Liu, Y., and M. K. Sen, A new time-space domain high-order finite-difference method for the acoustic wave equation, J. Computational Phys., 228, 8779-8806, 2009, 5 citations, doi:10.1016/j.jcp.2009.08.027, #2082
A new unified methodology was proposed in Finkelstein and Kastner (2007) [39] to derive spatial finite-difference (FD) coefficients in the joint timeâspace domain to reduce numerical dispersion. The key idea of this method is that the dispersion relation is completely satisfied at several designated frequencies. We develop this new timeâspace domain FD method further for 1D, 2D and 3D acoustic wave modeling using a plane wave theory and the Taylor series expansion. New spatial FD coefficients are frequency independent though they lead to a frequency dependent numerical solution. We prove that the modeling accuracy is 2nd-order when the conventional (2M)th-order space domain FD and the 2nd-order time domain FD stencils are directly used to solve the acoustic wave equation. However, under the same discretization, the new 1D method can reach (2M)th-order accuracy and is always stable. The 2D method can reach (2M)th-order accuracy along eight directions and has better stability. Similarly, the 3D method can reach (2M)th-order accuracy along 48 directions and also has better stability than the conventional FD method. The advantages of the new method are also demonstrated by the results of dispersion analysis and numerical modeling of acoustic wave equation for homogeneous and inhomogeneous acoustic models. In addition, we study the influence of the FD stencil length on numerical modeling for 1D inhomogeneous media, and derive an optimal FD stencil length required to balance the accuracy and efficiency of modeling. A new timeâspace domain high-order staggered-grid FD method for the 1D acoustic wave equation with variable densities is also developed, which has similar advantages demonstrated by dispersion analysis, stability analysis and modeling experiments. The methodology presented in this paper can be easily extended to solve similar partial difference equations arising in other fields of science and engineering.
Liu, Y., and M. K. Sen, Advanced finite-difference methods for seismic modeling, Geohorizons, 14 (2), 5-15, 2009, #2084
The finite-difference methods (FDMs) have been widely used in seismic modeling and migration. In this paper, we review the conventional arbitrary-order explicit FDMs and their recent developments, including arbitrary even-order implicit FDMs for standard grids and arbitrary even-order time-space domain FDMs for acoustic wave equations. For a given accuracy, an arbitrary even-order FDM can provide a trade-off between the order number and grid size. These explicit FDMs are the most popular in seismic modeling community. The finite-difference methods that are implicit in space are not very common because they generally require more computer resources than the explicit FDMs. Here we show that a new class of implicit FDMs can be derived that require solving a tri-diagonal system, which makes the resulting algorithm computationally very efficient. Therefore, some high-order explicit FDM may be replaced by an implicit FDM of some order to decrease the computation time without affecting the accuracy. We further demonstrate that we can also derive the FD coefficients in the joint time-space domain. These high-order spatial finite-difference stencils designed in joint time-space domain, when used in acoustic wave equation modeling, can provide even greater accuracy than those designed in the space domain alone under the same discretization. We demonstrate performance of these algorithms using some realistic 2D numerical examples.
Liu, Y., and M. K. Sen, An implicit staggered-grid finite-difference method for seismic modeling, Geophys. J. Int., 179, 459-474, 2009, 4 citations, doi:10.1111/j.1365-246X.2009.04305.x, #2086
We derive explicit and new implicit staggered-grid finite-difference (FD) formulas for derivatives of first order with any order of accuracy by a plane wave theory and Taylor's series expansion. Furthermore, we arrive at a practical algorithm such that the tridiagonal matrix equations are formed by the implicit FD formulas derived from the fractional expansion of derivatives. Our results demonstrate that the accuracy of a (2N+ 2)th-order implicit formula is nearly equivalent to or greater than that of a (4N)th-order explicit formula. The new implicit method only involves solving tridiagonal matrix equations. We also demonstrate that a (2N+ 2)th-order implicit formulation requires nearly the same amount of memory and computation as those of a (2N+ 4)th-order explicit formulation but attains the accuracy achieved by a (4N)th-order explicit formulation when additional cost of visiting arrays is not considered. Our analysis of efficiency and numerical modelling results for elastic wave propagation demonstrates that a high-order explicit staggered-grid method can be replaced by an implicit staggered-grid method of some order, which will increase the accuracy but not the computational cost.
Liu, Y., and M. K. Sen, A practical implicit finite-difference method: Examples from seismic modeling, J. Geophysics Engin., 6, 231-249, 2009, 4 citations, doi:10.1088/1742-2132/6/3/003, #2089
We derive explicit and new implicit finite-difference formulae for derivatives of arbitrary order with any order of accuracy by the plane wave theory where the finite-difference coefficients are obtained from the Taylor series expansion. The implicit finite-difference formulae are derived from fractional expansion of derivatives which form tridiagonal matrix equations. Our results demonstrate that the accuracy of a (2N + 2)th-order implicit formula is nearly equivalent to that of a (6N + 2)th-order explicit formula for the first-order derivative, and (2N + 2)th-order implicit formula is nearly equivalent to (4N + 2)th-order explicit formula for the second-order derivative. In general, an implicit method is computationally more expensive than an explicit method, due to the requirement of solving large matrix equations. However, the new implicit method only involves solving tridiagonal matrix equations, which is fairly inexpensive. Furthermore, taking advantage of the fact that many repeated calculations of derivatives are performed by the same difference formula, several parts can be precomputed resulting in a fast algorithm. We further demonstrate that a (2N + 2)th-order implicit formulation requires nearly the same memory and computation as a (2N + 4)th-order explicit formulation but attains the accuracy achieved by a (6N + 2)th-order explicit formulation for the first-order derivative and that of a (4N + 2)th-order explicit method for the second-order derivative when additional cost of visiting arrays is not considered. This means that a high-order explicit method may be replaced by an implicit method of the same order resulting in a much improved performance. Our analysis of efficiency and numerical modelling results for acoustic and elastic wave propagation validates the effectiveness and practicality of the implicit finite-difference method.
Liu, Y., and M. K. Sen, Numerical modeling of wave equation by truncated high-order finite difference method, Earthquake Science, 22, 205-213, 2009, #2094
Ma, J. T., M. K. Sen, and X. Chen, Free-surface multiple attenuation using inverse data processing in the coupled plane-wave domain, Geophysics, 74, V75-V81, 2009, 1 citation, doi:10.1190/1.3137059, #2090
ree-surface multiples contain a large amount of energy in the seismic data because of large reflectivity of the free surface. We propose a method for free-surface multiple attenuation by a simple muting in the inverse coupled plane-wave domain. Our method is based on inverse data processing and the well-known 2D invariant embedding technique. If the lateral variation in subsurface structure is smooth, the data are well compressed in the 2D coupled plane-wave domain, reducing computation costs and stabilizing the inversion procedure. Surface multiples and primaries are well separated in the inverse coupled plane-wave domain, and multiples can be eliminated by simple muting, which does not damage the primary energy. To reduce artifacts, wraparound, and noise introduced by the frequency-domain data matrix inversion, horizontal and vertical tapers are applied. A least-squares matrix inversion method is chosen to stabilize the inversion. Synthetic data examples show that plane-wave inverse data processing is stable and successful in attenuating free-surface multiples.
Ma, J. T., M. K. Sen, and X. Chen, Multiple attenuation using inverse data processing in the plane wave domain, (in Chinese with English abstract), Chinese J. Geophys., 52, 808-816, 2009, 1 citation, #2095
Sen, M. K., and A. Pal, A reflectivity method for laterally varying media: Homogeneous layers with curved interfaces, Geophys. J. Int., 178, 792-812, 2009, doi:10.1111/j.1365-246X.2009.04172.x, #2093
1-D reflectivity method in the frequency ray-parameter domain is the most popular method for computing synthetic seismograms. We explore an extension of this method to media consisting of homogeneous layers with curved interfaces. The 2-D extension of the method involves explicit evaluation of boundary conditions that results in a matrix formulation involving several matrix inversions in a coupled ray-parameter domain. First, we revisit this method by examining the intermediate results in the coupled shot-receiver ray-parameter domain for realistic models; numerous numerical artefacts are caused by matrix operations. Next we investigate an alternate asymptotic approach by examining a tangent plane or Kirchhoff formulation in the plane wave domain to replace the exact boundary condition evaluation. All other reflectivity operations such as the invariant imbedding or iterative computation of reflection, transmission, multiple and mode conversions can be readily applied even under this approximation. Our resulting algorithm computes noise free seismograms even with coarse sampling of the interfaces and ray parameters. Approximate calculation of reflection/transmission coefficients, however, does not include multiple interaction of a plane wave with an interface.
We show examples of synthetic seismograms for a suite of realistic models. Our computations are performed in the coupled slowness domain and thus, multiple shot-receiver data can be synthesized rapidly. Intermediate results in the coupled slowness space provide important insight into understanding wave propagation in heterogeneous media.
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.
Sil, S., and M. K. Sen, Seismic critical-angle anisotropy analysis in τ-p domain, Geophysics, 74, A53-A57, 2009, 1 citation, doi:10.1190/1.3137019, #2088
Seismic critical-angle reflectometry is a relatively new field for estimating seismic anisotropy parameters. The theory relates changes in the critical angle with azimuth of the seismic line to the principal axis and anisotropy parameters. Current implementation of the critical-angle reflectometry process has certain shortcomings in that the critical angle is determined from critical offset and the process is vulnerable to different approximation errors. Seismic critical-angle analysis in the plane-wave (-p) domain can handle these issues and has the potential to become an independent tool for estimating anisotropy parameters. The theory of seismic critical-angle reflectometry is modified to make it suitable for -p domain analysis. Then using full-wave synthetic seismograms at three different azimuths for a transversely isotropic medium with a horizontal axis of symmetry (HTI), the effectiveness of anisotropy parameter estimation is demonstrated. ©2009 Society of Exploration Geophysicists
Sil, S., and M. K. Sen, Azimuthal tau-p analysis in a weak orthorhombic medium, J. Seismic Explor., 18, 81-91, 2009, 3 citations, #2099
Srinivasan, S., and M. K. Sen, Stochastic modeling of facies distribution in a carbonate reservoir in the Gulf of Mexico, Geohorizons, 14 (2), 54-67, 2009, #2085
We examine the basis for classification of lithofacies in the Jurassic interval of a Gulf of Mexico reservoir. Previous attempts at facies classification based on the well logs alone yielded mixed results. While some facies exhibited distinct signatures, a majority of geological facies remained unresolved using well log attributes. The correspondence between seismic facies defined on the basis of multiple seismic attributes and those observed in cores was also poor. It was therefore deemed necessary to perform seismic inversion and derive acoustic impedance volumes. We analyzed the impedance volume to define three distinct types of flow facies which were subsequently utilized within a geostatistical model for interpolating facies between wells constrained to the impedance volume. The permanence of ratios hypothesis for data integration was implemented. Our results exhibit considerable refinement over the previous models defined solely on the basis of well attributes.
Srivastava, R. P., and M. K. Sen, Fractal-based stochastic inversion of poststack seismic data using very fast simulated annealing, J. Geophysics Engin., 6, 412-425, 2009, 1 citation, doi:10.1088/1742-2132/6/4/009, #2143
Seismic data do not contain some low- and high-frequency information because of the band-limited nature of the source wavelet. A deterministic inversion of such band-limited seismic data produces smooth models which are devoid of high-frequency variations observed in well logs. Stochastic inversion methods often based on random Gaussian priors can have a limitation of producing high frequencies in the desired model particularly the frequency band not constrained by the input seismic data. In this paper, we propose a new stochastic poststack inversion algorithm where fractal models constructed from statistical properties of well logs are used to generate a priori models. This provides a high-resolution model without injecting spurious high-frequency estimates in model space. Stacked seismic data are used in the inversion in which a suitable objective function is minimized using a nonlinear optimization method called 'very fast simulated annealing'. We demonstrate the effectiveness of our method for the estimation of acoustic impedance with the application to a field dataset.
Vedanti, N., and M. K. Sen, Seismic inversion tracks in situ combustion: A case study from Baloi oil field, India, Geophysics, 74, B103-B112, 2009, 4 citations, doi:10.1190/1.3129262, #2087
In situ combustion is one of the approaches used for secondary recovery of heavy oil in which time-lapse seismic data might be used to track the production process. We have analyzed three time-lapse seismic surveys carried out at a time interval of one year during pre- and postcombustion phases of in situ combustion at the Balol field in India. Interpretation of these land time-lapse data using standard seismic amplitude differences in terms of reservoir changes during production and injection can be erroneous. The Balol data, in particular, lacked calibration and had poor repeatability. We have addressed these issues by demonstrating the applicability of independent prestack inversion of baseline and monitor surveys to derive elastic attributes that help to produce clearer images of fluid movement. Independent estimates of wavelets from the baseline and two monitor surveys were used in prestack inversion to estimate acoustic impedance, shear impedance, and Poisson's ratio. The results show unexpected movement of the thermal front away from the injection wells toward the north and northwest of the injection wells, farther from the production wells. These results have been confirmed by well production data from the field north and northwest of Balol. ©2009 Society of Exploration Geophysicists
Bansal, R., and M. K. Sen, Finite difference modeling of S-wave splitting in anisotropic media, Geophys. Prospecting, 56, 293-312, 2008, 12 citations, doi:10.1111/j.1365-2478.2007.00693.x, #2108
We have implemented a 3D finite-difference scheme to simulate wave propagation in arbitrary anisotropic media. The anisotropic media up to orthorhombic symmetry were modelled using a standard staggered grid scheme and beyond (monoclinic and triclinic) using a rotated staggered grid scheme. The rationale of not using rotated staggered grid for all types of anisotropic media is that the rotated staggered grid schemes are more expensive than standard staggered grid schemes. For a 1D azimuthally anistropic medium, we show a comparison between the seismic data generated by our finite-difference code and by the reflectivity algorithm; they are in excellent agreement.
We conducted a study on zero-offset shear-wave splitting using the finite-difference modelling algorithm using the rotated staggered grid scheme. Our S-wave splitting study is mainly focused on fractured media. On the scale of seismic wavelenghts, small aligned fractures behave as an equivalent anisotropic medium. We computed the equivalent elastic properties of the fractures and the background in which the fractures were embedded, using low-frequency equivalent media theories. Wave propagation was simulated for both rotationally invariant and corrugated fractures embedded in an isotropic background for one, or more than one, set of fluid-filled and dry fractures. S-wave splitting was studied for dipping fractures, two vertical non-orthogonal fractures and corrugated fractures. Our modelling results confirm that S-wave splitting can reveal the fracture infill in the case of dipping fractures. S-wave splitting has the potential to reveal the angle between the two vertical fractures. We also notice that in the case of vertical corrugated fractures, S-wave splitting is sensitive to the fracture infill.
De Basabe, J. D., M. K. Sen, and M. F. Wheeler, The interior penalty discontinuous Galerkin method for elastic wave propagtion: Grid dispersion, Geophys. J. Int., 175, 83-93, 2008, 10 citations, doi:10.1111/j.1365-246X.2008.03915.x, #2106
Recently, there has been an increased interest in applying the discontinuous Galerkin method (DGM) to wave propagation. In this work, we investigate the applicability of the interior penalty DGM to elastic wave propagation by analysing it's grid dispersion properties, with particular attention to the effect that different basis functions have on the numerical dispersion. We consider different types of basis functions that naturally yield a diagonal mass matrix. This is relevant to seismology because a diagonal mass matrix is tantamount to an explicit and efficient time marching scheme. We find that the Legendre basis functions that are traditionally used in the DGM introduce numerical dispersion and anisotropy. Furthermore, we find that using Lagrange basis functions along with the Gauss nodes has attractive advantages for numerical wave propagation.
Dutta, U., M. K. Sen, N. N. Biswas, and Y. Zhang, Waveform modeling of earthquake data for delineating the sedimentary structure of the Anchorage basin, Alaska, 14th World Conference on Earthquake Engineering, Beijing, China, 2008, #2083
Filina, I. Y., D. D. Blankenship, M. Thoma, V. V. Lukin, V. N. Masolov, and M. K. Sen, New 3D bathymetry and sediment distribution in Lake Vostok: Implication for pre-glacial origin and numerical modeling of the internal processes within the lake, Earth Planet. Sci. Lett., 276, 106-114, 2008, 11 citations, doi:10.1016/j.epsl.2008.09.012 , #2103
A new distribution of water and unconsolidated sediments in subglacial Lake Vostok, East Antarctica was developed via inversion of airborne gravity data constrained by 60 seismic soundings. A model was developed for host rock with a density of 2550 kg/m3 that was inferred from prior 2D modeling. Our 3D bathymetry model of Lake Vostok corresponds better with seismic data (RMS of 125 m) than two previous models based on the same gravity dataset. The good match in both water and sediment thicknesses between the gravity model and seismic measurements confirms two major facts about Lake Vostok: (1) the lake is hosted by sedimentary rocks, and (2) the bottom of the lake is covered with a layer of unconsolidated sediments that does not exceed 300 m in the southern basin and thickens almost to 400 m in the northern basin. Our new bathymetry model suggests much shallower water thicknesses (up to twice the previous estimates) in the middle and northern parts of the lake, while the water layer is thicker in the southern basin. Numerical modeling of the internal processes in the lake reveals the relevance of our new bathymetry model to the basal mass balance. A significant decrease in transport is observed in the shallower northern basin, as well as a decrease of 33% in the turbulent kinetic energy. However, only minor differences were observed in the distribution of the calculated freezing and melting zones compared to previous models. Estimates for the sedimentation rates for six possible mechanisms were made. Possible sedimentation mechanisms are: (1) fluvial and periglacial, i.e. those that are active prior to the establishment of a large subglacial lake; (2) deposition due to overlying ice sheet, including melting out of the ice, as well as bulldozering by the overriding ice; and (3) suspended sediments from subglacial water flow including those deposited by periodical subglacial outbursts. The estimates for these mechanisms show that unconsolidated sediments of the observed thickness are most consistent with a lake that existed before glaciation.
Gangopadhyay, A., and M. K. Sen, A possible mechanism for the spatial distribution of seismicity in northern Gulf of Mexico, Geophys. J. Int., 175, 1141-1153, 2008, 2 citations, doi:10.1111/J.1365-246X.2008.03952.x, #2053
In an aseismic part of northern Gulf of Mexico (GOM), three earthquakes of Ms 5.2, Mw 4.6 and Mw 5.8 occurred between 2006 February and September. These earthquakes raised concerns because of their proximity to locations of active hydrocarbon production, and potential for tsunami generation, prompting investigations about their causes. Herein, we test a hypothesis using 2-D and 3-D models, where stress concentration due to contrast in mechanical properties between the salt deposits and surrounding sediments in the northern GOM, driven by background tectonic loading, may have triggered these earthquakes in particular, and cause earthquakes in the northern GOM in general. We perform the mechanical modelling by a 'distinct element method' implemented through commercially available computer programmes called 'Universal Distinct Element Code (UDEC)' and 'Three-Dimensional Distinct Element Code (3DEC)'. Our models consist of a simplified regional outline of the salt bodies surrounded by sediments. We assign mechanical properties to the salt and sediments based on known geology and seismic velocities. We load the models tectonically for a year, and observe the patterns of resulting shear stresses. The results of modelling suggest that some locations of relatively high shear stress correlate well with the spatial distribution of seismicity in the northern GOM, thereby suggesting a possible causal association.
Jackson, C. S., M. K. Sen, G. Huerta, Y. Deng, and K. P. Bowman, Error reduction and convergence in climate prediction, J. Climate, 21, 6698-6709, 2008, 8 citations, doi:10.1175/2008JCLI2112.1, #2029
Although climate models have steadily improved their ability to reproduce the observed climate, over the years there has been little change to the wide range of sensitivities exhibited by different models to a doubling of atmospheric CO2 concentrations. Stochastic optimization is used to mimic how six independent climate model development efforts might use the same atmospheric general circulation model, set of observational constraints, and model skill criteria to choose different settings for parameters thought to be important sources of uncertainty related to clouds and convection. Each optimized model improved its skill with respect to observations selected as targets of model development. Of particular note were the improvements seen in reproducing observed extreme rainfall rates over the tropical Pacific, which was not specifically targeted during the optimization process. As compared to the default model sensitivity of 2.4°C, the ensemble of optimized model configurations had a larger and narrower range of sensitivities around 3°C but with different regional responses related to the uncertain choice in optimized parameter settings. These results suggest current generation models, if similarly optimized, may become more convergent in their measure of global sensitivity to greenhouse gas forcing. However, this exploration of the possible sources of modeling and observational uncertainty is not exhaustive. The optimization process illustrates an objective means for selecting an ensemble of plausible climate model configurations that quantify a portion of the uncertainty in the climate model development process.
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
Kumar, C., M. K. Sen, and R. J. Ferguson, Depth migration anisotropy analysis in the time domain, Geophys. Prospecting, 56, 87-94, 2008, doi:10.1111/j.1365-2478.2007.00658.x, #1921
In areas of complex geology such as the Canadian Foothills, the effects of anisotropy are apparent in seismic data and estimation of anisotropic parameters for use in seismic imaging is not a trivial task. Here we explore the applicability of common-focus point (CFP)-based velocity analysis to estimate anisotropic parameters for the variably tilted shale thrust sheet in the Canadian Foothills model. To avoid the inherent velocity-depth ambiguity, we assume that the elastic properties of thrust-sheet with respect to transverse isotropy symmetry axis are homogeneous, the reflector below the thrust-sheet is flat, and that the anisotropy is weak. In our CFP approach to velocity analysis, for a poorly imaged reflection point, a traveltime residual is obtained as the time difference between the focusing operator for an assumed subsurface velocity model and the corresponding CFP response obtained from the reflection data. We assume that this residual is due to unknown values for anisotropy, and we perform an iterative linear inversion to obtain new model parameters that minimize the residuals. Migration of the data using parameters obtained from our inversion results in a correctly positioned and better focused reflector below the thrust sheet. For traveltime computation we use a brute force mapping scheme that takes into account weakly tilted transverse isotropy media. For inversion, the problem is set up as a generalized Newton's equation where traveltime error (differential time shift) is linearly dependent on the parameter updates. The iterative updates of parameters are obtained by a least-squares solution of Newton's equations. The significance of this work lies in its applicability to areas where transverse isotropy layers are heterogeneous laterally, and where transverse isotropy layers are overlain by complex structures that preclude a moveout curve fitting.
Sen, M. K., and N. Vedanti, Seismic inversion tracks in-situ combustion: 4D seismic at the Baloi field, India, Soc. Petroleum Engin., SPE-116600-PP, 2008, #2105
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.
Sil, S., and M. K. Sen, Azimuthal τ-p analysis in anisotropic media, Geophys. J. Int., 175, 587-597, 2008, 5 citations, doi:10.1111/j.1365-246X.2008.03882.x, #2104
For the purpose of transversely isotropic (TI) normal moveout (NMO) correction, we propose analysis of plane wave transformed azimuthal gathers, interactively using a single azimuth data at a time and a new delay time equation (developed in this paper), which is a function of two parameters at each azimuth. Results from independently estimated multi-azimuth gathers, then, can be combined to estimate stiffness or Thomsen coefficients. Azimuthal τâp analysis also avoids numerical ray tracing, resulting in a rapid algorithm. We demonstrate the applicability of our method using a set of P-wave synthetic seismograms from a multilayered medium, consisting of isotropic and HTI layers. Azimuth-dependent anisotropy parameters are derived by delay time fitting and NMO correction. The reflections from the bottom interface of an isotropic layer with an anisotropic overburden show apparent anisotropic traveltime behaviour, which is easily accounted for by our layer-stripping based azimuthal NMO analysis. Unlike the previous approximate HTI NMO correction equation, this equation performs better NMO correction for the HTI medium and is also applicable to the VTI medium. Presence of only two reduced parameters in the equation helps the anisotropic parameter estimation become less ambiguous.
Van Der Neut, J., M. K. Sen, and K. Wapenaar, Seismic reflection coefficients of faults at low frequencies: A model study, Geophys. Prospecting, 56, 287-292, 2008, 1 citation, doi:10.1111/j.1365-2478.2008.00701.x, #2109
We use linear slip theory to evaluate seismic reflections at non-welded interfaces, such as faults or fractures, sandwiched between general anisotropic media and show that at low frequencies the real parts of the reflection coefficients can be approximated by the responses of equivalent welded interfaces, whereas the imaginary parts can be related directly to the interface compliances. The imaginary parts of low frequency seismic reflection coefficients at fault zones can be used to estimate the interface compliances, which can be related to fault properties upon using a fault model. At normal incidence the expressions uncouple and the complex-valued P-wave reflection coefficient can be related linearly to the normal compliance. As the normal compliance is highly sensitive to the infill of the interface, it can be used for gas/fluid identification in the fault plane. Alternatively, the tangential compliance of a fault can be estimated from the complex-valued S-wave reflection coefficient. The tangential compliance can provide information on the crack density in a fault zone. Coupling compliances can be identified and quantified by the observation of PS conversion at normal incidence, with a comparable linear relationship.
Villagran, A., G. Huerta, C. S. Jackson, and M. K. Sen, Computational methods for parameter estimation in climate models, Bayesian Analysis, 3, 823-850, 2008, 3 citations, doi:10.1214/08-BA331, #2101
Intensive computational methods have been used by Earth scientists in a wide range of problems in data inversion and uncertainty quantification such as earthquake epicenter location and climate projections. To quantify the uncer- tainties resulting from a range of plausible model configurations it is necessary to estimate a multidimensional probability distribution. The computational cost of estimating these distributions for geoscience applications is impractical using traditional methods such as Metropolis/Gibbs algorithms as simulation costs limit the number of experiments that can be obtained reasonably. Several alternate sampling strategies have been proposed that could improve on the sampling effi- ciency including Multiple Very Fast Simulated Annealing (MVFSA) and Adaptive Metropolis algorithms. The performance of these proposed sampling strategies are evaluated with a surrogate climate model that is able to approximate the noise and response behavior of a realistic atmospheric general circulation model (AGCM). The surrogate model is fast enough that its evaluation can be embed- ded in these Monte Carlo algorithms. We show that adaptive methods can be superior to MVFSA to approximate the known posterior distribution with fewer forward evaluations. However the adaptive methods can also be limited by inad- equate sample mixing. The Single Component and Delayed Rejection Adaptive Metropolis algorithms were found to resolve these limitations, although challenges remain to approximating multi-modal distributions. The results show that these advanced methods of statistical inference can provide practical solutions to the cli- mate model calibration problem and challenges in quantifying climate projection uncertainties. The computational methods would also be useful to problems out- side climate prediction, particularly those where sampling is limited by availability of computational resources.
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.
De Basabe, J. D., and M. K. Sen, Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equations, Geophysics, 72, T81-T95, 2007, 23 citations, doi:10.1190/1.2785046, #1930
Purely numerical methods based on finite-element approximation of the acoustic or elastic wave equation are becoming increasingly popular for the generation of synthetic seismograms. We present formulas for the grid dispersion and stability criteria for some popular finite-element methods (FEM) for wave propagation, namely, classical and spectral FEM. We develop an approach based on a generalized eigenvalue formulation to analyze the dispersive behavior of these FEMs for acoustic or elastic wave propagation that overcomes difficulties caused by irregular node spacing within the element and the use of high-order polynomials, as is the case for spectral FEM. Analysis reveals that for spectral FEM of order four or greater, dispersion is less than 0.2% at four to five nodes per wavelength, and dispersion is not angle dependent. New results can be compared with grid-dispersion results of some classical finite-difference methods (FDM) used for acoustic or elastic wave propagation. Analysis reveals that FDM and classical FEM require a larger sampling ratio than a spectral FEM to obtain results with the same degree of accuracy. The staggered-grid FDM is an efficient scheme, but the dispersion is angle dependent with larger values along the grid axes. On the other hand, spectral FEM of order four or greater is isotropic with small dispersion, making it attractive for simulations with long propagation times.
Gangopadhyay, A., J. Pulliam, and M. K. Sen, Waveform modeling of teleseismic S, SP, SsPmP, and shear-coupled PL waves for crust- and upper-mantle velocity structure beneath Africa, Geophys. J. Int., 170, 1210-1226, 2007, 2 citations, doi:10.1111/j.1365-246X.2007.03470.x, #1911
We describe a waveform modelling technique and demonstrate its application to determine the crust- and upper-mantle velocity structure beneath Africa. Our technique uses a parallelized reflectivity method to compute synthetic seismograms and fits the observed waveforms by a global optimization technique based on a Very Fast Simulated Annealing (VFSA). We match the S, Sp, SsPmP and shear-coupled PL phases in seismograms of deep (200â800 km), moderate-to-large magnitude (5.5â7.0) earthquakes recorded teleseismically at permanent broad-band seismic stations in Africa. Using our technique we produce P- and S-wave velocity models of crust and upper mantle beneath Africa. Additionally, our use of the shear-coupled PL phase, wherever observed, improves the constraints for lower crust- and upper-mantle velocity structure beneath the corresponding seismic stations. Our technique retains the advantages of receiver function methods, uses a different part of the seismogram, is sensitive to both P- and S-wave velocities directly, and obtains helpful constraints in model parameters in the vicinity of the Moho. The resulting range of crustal thicknesses beneath Africa (21â46 km) indicates that the crust is thicker in south Africa, thinner in east Africa and intermediate in north and west Africa. Crustal P- (4.7â8 km s−1) and S-wave velocities (2.5â4.7 km s−1) obtained in this study show that in some parts of the models, these are slower in east Africa and faster in north, west and south Africa. Anomalous crustal low-velocity zones are also observed in the models for seismic stations in the cratonic regions of north, west and south Africa. Overall, the results of our study are consistent with earlier models and regional tectonics of Africa.
Gangopadhyay, A., J. Pulliam, and M. K. Sen, Modeling of teleseismic waveforms for crust and upper mantle velocity structure, in National Conference on Modern Trends in Geophysical Sciences and Techniques, Indian School of Mines, Jharkhand, India, 2007, #1977
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
Kumar, D., M. K. Sen, and N. L. Bangs, Gas hydrate concentration and characteristics within Hydrate Ridge inferred from multicomponent seismic reflection data, J. Geophys. Res., 112, B12306, 2007, 6 citations, doi:10.1029/2007JB004993, #1782
A seismic experiment composed of streamer and ocean bottom seismometer (OBS) surveys was conducted in the summer of 2002 at southern Hydrate Ridge, offshore Oregon, to map the gas hydrate distribution within the hydrate stability zone. Gas hydrate concentrations within the reservoir can be estimated with P wave velocity (V p ); however, we can further constrain gas hydrate concentrations using S wave velocity (V s ), and use V s through its relationship to V p (V p /V s ) to reveal additional details such as gas hydrate form within the matrix (i.e., hydrate cements the grains, becomes part of the matrix frame or floats in pore space). Both V p and V s can be derived simultaneously by inverting multicomponent seismic data. In this study, we use OBS data to estimate seismic velocities where both gas hydrate and free gas are present in the shallow sediments. Once V p and V s are estimated, they are simultaneously matched with modeled velocities to estimate the gas hydrate concentration. We model V p using an equation based on a modification of Wood's equation that incorporates an appropriate rock physics model and V s using an empirical relation. The gas hydrate concentration is estimated to be up to 7% of the rock volume, or 12% of the pore space. However, V p and V s do not always fit the model simultaneously. V p can vary substantially more than V s . Thus we conclude that a model, in which higher concentrations of hydrate do not affect shear stiffness, is more appropriate. Results suggest gas hydrates form within the pore space of the sediments and become part of the rock framework in our survey area.
Pulliam, J., M. K. Sen, and A. Gangopadhyay, Waveform modeling of teleseismic converted waves for crust and upper mantle structure beneath Canada and China, in Proc. 29th Monitoring Res. Rev: Ground-Based Nuclear Explosion Monitoring Technologies, Denver, CO, 1-22, 2007, #1975
Sen, M. K., F. D. Lane, and D. J. Foster, Anomalous reflection amplitudes from fractured reservoirs- Failure of AVOA?, Leading Edge, 26, 1148-1152, 2007, doi:10.1190/1.2780785, #1926
A fractured reservoir embedded in an isotropic background is equivalent to an anisotropic medium whose symmetry depends on fracture orientation. For this configuration, a variation of PP-reflection amplitude as a function of azimuth, at a given incidence angle is often observed. Here we report on an observation of unusual amplitude behavior for reflections from the top Cadotte and Falher formations in the Lynx Field, Western Canadian Sedimentary Basin, Alberta, Canada. Formation microscanner (FMS) images suggest considerable fracturing in the shallower Cadotte interval and resistivity measurements show fracturing in the Falher zone. Interestingly, reflection amplitudes from the Cadotte do not show the expected azimuthal variation consistent with a vertically fractured medium. Reflections from the deeper Falher interval, however, do exhibit large, unexpected azimuthal variation.
Shaw, R. K., A. Chatterjee, and M. K. Sen, Use of AVOA analysis to estimate dip of obliquely dipping fractures in hydrocarbon reservoirs, J. Seismic Explor., 15, 355-365, 2007, #1927
Bansal, R., and M. K. Sen, Two-point ray tracing in anisotropic media, J. Seismic Explor., 15, 25-44, 2006, #1855
Filina, I. Y., D. D. Blankenship, L. Roy, M. K. Sen, T. G. Richter, and J. W. Holt, Inversion of airborne gravity data acquired over subglacial lakes in East Antarctica, Antarctica: Contributions to Global Earth Sciences, edited by D. K. Futterer, D. Damaske, G. Kleinschmidt,H. Miller, and F. Tessensohn, Springer-Verlag, Heidelberg, Germany, 129-133, 2006, 5 citations, #1714
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.
Kumar, D., M. K. Sen, N. L. Bangs, C. Wang, and I. A. Pecher, Seismic anisotropy at Hydrate Ridge, Geophys. Res. Lett., 33, L01306, 2006, 9 citations, doi:10.1029/2005GL023945, #1843
P-wave velocity increases in the presence of gas hydrates and decreases in the presence of free gas in the sediments, making it an excellent means to investigate gas hydrate systems. However, seismic velocity is typically derived from surface seismic data without consideration of seismic anisotropy. The presence of anisotropy in the hydrate bearing sediments adds an additional complexity in data analysis; however anisotropy can help reveal the distribution of hydrates. Here we report on the evidence of seismic anisotropy at Hydrate Ridge along the Cascadia convergent margin. We find that the south summit is anisotropic, while the basin side (east of south summit) is isotropic. Anisotropy is likely caused by the hydrate veins. We interpret the anisotropy parameters in terms of the distribution and fabric of gas hydrates.
Kumar, D., M. K. Sen, and N. L. Bangs, Seismic characteristics of gas hydrates at Hydrate Ridge, offshore Oregon, Leading Edge, 25, 610-614, 2006, doi:10.1190/1.2202665, #1844
Gas hydrate is an ice-like substance that contains low molecular weight gases (mostly methane) in a lattice of water molecules (Sloan, 1998). Gas hydrates are widely present in permafrost and deep oceanic environments around the world. Hydrate Ridge, offshore Oregon, is one of the areas where hydrates are found in relatively high concentration. In marine environments, methane hydrates are usually stable at temperatures of 0ââ¬â20ð C, water depths greater than 500 m, and in sediments up to 300 m below the seafloor. Small amounts of free gas are often present below the gas-hydrate stability zone (GHSZ). The hydrate-to-free-gas contact gives a strong acoustic impedance contrast, which is evident on seismic sections as a bottom-simulating reflection (BSR). Offshore gas hydrate systems have received attention from the scientific community because of their potential to be
Kumar, C., M. K. Sen, and R. J. Ferguson, Traveltime sensitivity analysis for parameter estimation in TI media, J. Seismic Explor., 15, 119-131, 2006, #1854
Pulliam, J., M. K. Sen, and A. Gangopadhyay, Determination of crust and upper mantle structure beneath Africa using a global optimization based waveform modeling technique, Proc. 28th Seis. Res. Rev., 196-208, 2006, #1891
Sen, M. K., Seismic Inversion, Society of Petroleum Engineers, Richardson, TX, 120 pp., 2006, 4 citations, #1853
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.
Shaw, R. K., and M. K. Sen, Use of AVOA data to estimate fluid indicator in a vertically fractured medium, Geophysics, 71, C15-C24, 2006, 3 citations, doi:10.1190/1.2194896, #1856
Microstructural attributes of cracks and fractures, such as crack density, aspect ratio, and fluid infill, determine the elastic properties of a medium containing a set of parallel, vertical fractures. Although the tangential weakness T of the fractures does not vary with the fluid content, the normal weakness N exhibits significant dependence on fluid infill. Based on linear-slip theory, we used the ratio gN/T â termed the fluid indicator â as a quantitative measure of the fluid content in the fractures, with g representing the square of the ratio of S- and P-wave velocity in the unfractured medium. We used a Born formalism to derive the sensitivity to fracture weakness of PP- and PS-reflection coefficients for an interface separating an unfractured medium from a vertically fractured medium. Our formulae reveal that the PP-reflection coefficient does not depend on the 2D microcorrugation/surface roughness with ridges and valleys parallel to the fracture strike, whereas the PS-reflection coefficient is sensitive to this microstructural property of the fractures. Based on this formulation, we developed a method to compute the fluid indicator from wide-azimuth PP-AVOA data. Inversion of synthetic data corrupted with 10% random noise reliably estimates the normal and tangential fracture weaknesses and hence the fluid indicator can be determined accurately when the fractures are liquid-filled or partially saturated. As the gas saturation in the fractures increases, the quality of inversion becomes poorer. Errors of 15%â20% in g do not affect the estimation of fluid indicator significantly in case of liquid infill or partial saturation. However, for gas-saturated fractures, incorrect values of g may have a significant effect on fluid-indicator estimates.
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
Varela, O. J., C. Torres-Verdin, and M. K. Sen, Enforcing smoothness and assessing uncertainty in non-linear one-dimensional prestack seismic inversion, Geophys. Prospecting, 54, 239-259, 2006, 4 citations, doi:10.1111/j.1365-2478.2006.00531.x, #1857
Estimation of elastic properties of rock formations from surface seismic amplitude measurements remains a subject of interest for the exploration and development of hydrocarbon reservoirs. This paper develops a global inversion technique to estimate and appraise 1D distributions of compressional-wave velocity, shear-wave velocity and bulk density, from normal-moveout-corrected PP prestack surface seismic amplitude measurements. Specific objectives are: (a) to evaluate the efficiency of the minimization algorithm (b) to appraise the impact of various data misfit functions, and (c) to assess the effect of the degree and type of smoothness criterion enforced by the inversion. Numerical experiments show that very fast simulated annealing is the most efficient minimization technique among alternative approaches considered for global inversion. It is also found that an adequate choice of data misfit function is necessary for a reliable and efficient match of noisy and sparse seismic amplitude measurements. Several procedures are considered to enforce smoothness of the estimated 1D distributions of elastic parameters, including predefined quadratic measures of length, flatness and roughness.
Based on the general analysis of global inversion techniques, we introduce a new stochastic inversion algorithm that initializes the search for the minimum with constrained random distributions of elastic parameters and enforces predefined autocorrelation functions (semivariograms). This strategy readily lends itself to the assessment of model uncertainty. The new global inversion algorithm is successfully tested on noisy synthetic amplitude data. Moreover, we present a feasibility analysis of the resolution and uncertainty of prestack seismic amplitude data to infer 1D distributions of elastic parameters measured with wireline logs in the deepwater Gulf of Mexico. The new global inversion algorithm is computationally more efficient than the alternative global inversion procedures considered here.
Varela, O. J., C. Torres-Verdin, M. K. Sen, and I. G. Roy, Using time-lapse seismic data to detect variations of pore pressure and fluid saturation due to oil displacement by water: A numerical study based on one-dimensional pre-stack inversion, J. Geophysics Engin., 3, 177-193, 2006, doi:10.1088/1742-2132/3/2/009, #1858
We describe a numerical study performed to appraise the ability of seismic amplitude data to infer the time evolution of pore pressure and fluid saturation due to hydrocarbon production. To this end, we construct a synthetic, spatially heterogeneous hydrocarbon reservoir model that is subject to numerical simulation of multiphase fluid flow. Hydrocarbon production is assumed in the form of one water-injection well and four oil-producing wells. The synthetic reservoir model exhibits average porosities of 20% and poses significant vertical resolution constraints to the usage of seismic amplitude data to ascertain variations of pore pressure and fluid saturation. We assume the availability of migrated prestack seismic amplitude and use one-dimensional models to simulate trace by trace the seismic amplitude data before and after the onset of production. One-dimensional seismic amplitude data are simulated in time-lapse mode making use of a rock physics model that includes the effect of differential compaction between sands and shales as a function of depth of burial. The sensitivity study presented in this paper is based on one-dimensional inversion and hence sheds light on the vertical resolution properties of noisy seismic amplitude data. Multiphase fluid-flow parameters have a measurable impact on fluid saturation and pore pressure and hence on the spatial distribution and time evolution of elastic parameters. However, the inverted spatial distributions of elastic parameters at best correlate with smooth spatial averages of the actual distributions of pore pressure and fluid saturation. Because for the case under consideration time-lapse seismic amplitude variations are of the order of 5%, such a correlation would be difficult, if not impossible, to ascertain without the use of one-dimensional inversion. We show that the elastic parameters inverted from prestack seismic amplitude data provide more degrees of freedom to discriminate between time variations of pore pressure and fluid saturation in the reservoir compared to distributions of acoustic impedance inverted from poststack seismic amplitude data.
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
Pulliam, J., and M. K. Sen, Assessing uncertainties in waveform modeling of the crust and upper mantle, Proc. 27th Seis. Res. Rev.: Ground-Based Nuclear Explosion Monitoring Technologies, Palm Springs, CA, 152-160, 2005, #1852
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, I. G., M. K. Sen, and C. Torres-Verdin, Full waveform seismic inversion using a distributed system of computers, Concurrency and Computation: Practice and Experience, 17, 1365-1385, 2005, 1 citation, doi:10.1002/cpe.897, #1683
The aim of seismic waveform inversion is to estimate the elastic properties of the Earth's subsurface layers from recordings of seismic waveform data. This is usually accomplished by using constrained optimization often based on very simplistic assumptions. Full waveform inversion uses a more accurate wave propagation model but is extremely difficult to use for routine analysis and interpretation. This is because computational difficulties arise due to: (1) strong nonlinearity of the inverse problem; (2) extreme ill-posedness; and (3) large dimensions of data and model spaces. We show that some of these difficulties can be overcome by using: (1) an improved forward problem solver and efficient technique to generate sensitivity matrix; (2) an iteration adaptive regularized truncated Gauss-Newton technique; (3) an efficient technique for matrix-matrix and matrix-vector multiplication; and (4) a parallel programming implementation with a distributed system of processors. We use a message-passing interface in the parallel programming environment. We present inversion results for synthetic and field data, and a performance analysis of our parallel implementation.
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
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.
Ferguson, R. J., and M. K. Sen, Estimating the elastic parameters of anisotropic media using a joint inversion of P-wave and SV-wave traveltime error, Geophys. Prospecting, 52, 547-557, 2004, 3 citations, doi:10.1111/j.1365-2478.2004.00450.x, #1726
A method is presented to estimate the elastic parameters and thickness of media that are locally laterally homogeneous using P-wave and vertically polarized shear-wave (SV-wave) data. This method is a 'layer-stripping' technique, and it uses many aspects of common focal point (CFP) technology. For each layer, a focusing operator is computed using a model of the elastic parameters with which a CFP gather can be constructed using the seismic data. Assuming local homogeneity, the resulting differential time shifts (DTSs) represent error in the model due to anisotropy and error in thickness. In the (τ−p) domain, DTSs are traveltimes Δτ that connect error in layer thickness z, vertical slowness q, and ray parameter p. Series expansion is used to linearize Δτ with respect to error in the elastic parameters and thickness, and least-squares inversion is used to update the model.
For stability, joint inversion of P and SV data is employed and, as pure SV data are relatively rare, the use of mode-converted (PSV) data to represent SV in the joint inversion is proposed. Analytic and synthetic examples are used to demonstrate the utility and practicality of this inversion.
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.
Kumar, D., M. K. Sen, and R. J. Ferguson, Traveltime computation and prestack depth migration in tilted transversely isotropic media, Geophysics, 69, 37-44, 2004, 11 citations, doi:10.1190/1.16493733, #1656
The principal objective of our work is to develop a technique for prestack depth migration in tilted transversely isotropic (TTI) media in which the axis of symmetry is not vertical and may be spatially varying. Such models are required to image seismic data in geologically complex regions such as the Canadian Foothills. We have developed a 2D Kirchhoff integral-based migration algorithm in which the traveltime computation comprises the major task. Among the existing traveltime computation algorithms such as ray tracing with interpolation, ray bending, and eikonal solvers, a direct or a brute force approach of traveltime computation is generally highly robust. We have modified a direct method of first-arrival P-wave traveltime computation in 2D media that accounts for TTI. The algorithm requires that the group velocity be computed at each gridpoint, using either an analytic solution or by an approximate Fourier series expansion. The P-wave traveltime contours computed for complex geologic models show the pronounced effects from TTI media. Our results, using laboratory P-wave data collected over a physical model of an anisotropic thrust sheet, reveal that a 2D Kirchhoff migration based on our traveltime algorithm (TTI model) images the structure beneath the thrust sheet very well. The vertical transversely isotropic (assuming a vertical axis of symmetry) or isotropic imaging introduces false anticlinal structures. We compare our results with those obtained by a recursive-extrapolation method and find that our approach images the underside of one of the thrust sheets better.
Roy, I. G., M. K. Sen, C. Torres-Verdin, and O. J. Varela, Prestack inversion of a Gulf of Thailand ocean-bottom cable (OBC) data set, Geophysics, 69, 1470-1477, 2004, 2 citations, doi:10.1190/1.1836820, #1727
Shaw, R. K., and M. K. Sen, Born integral, stationary phase and linearized reflection coefficients in weak anisotropic media, Geophys. J. Int., 158, 225-238, 2004, 4 citations, doi:10.1111/j.1365-246X.2004.02283.x, #1684
We show that the linearized reflection coefficients for arbitrary anisotropic media embedded in an isotropic background can be derived directly from a Born formalism. Due to rapidly varying phases of the scattered waves from first-order perturbations in density and elastic parameters, the major contributions to the observed wavefield for any sourceâreceiver pair far from the volume of scatterers arise from the stationary points of a scattering integral, called the Born integral. For simple interface models, such integrals can be evaluated analytically using the method of stationary phase. The resulting scattering function relates linearly to the approximate (linearized) reflection coefficient through a scaling factor determined by the angle of incidence and the properties of the background medium. We consider a homogeneous isotropic background to express the approximate reflection coefficients as a sum of an isotropic and an anisotropic reflection coefficient. The isotropic coefficient is a weighted sum of density and isotropic perturbations about the background, whereas the anisotropic coefficient is a weighted sum of anisotropic perturbations where the weights depend on the angles of incidence, the properties of the background medium as well as the azimuth of the plane of reflection with respect to some symmetry plane of the weakly anisotropic medium.
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
Bhattacharya, B. B., Shalivahan, and M. K. Sen, Use of VFSA for resolution, sensitivity and uncertainty analysis in 1D DC resistivity and IP inversion, Geophys. Prospecting, 51, 393-408, 2003, 7 citations, doi:10.1046/j.1365-2478.2003.00379.x, #1657
We present results from the resolution and sensitivity analysis of 1D DC resistivity and IP sounding data using a non-linear inversion. The inversion scheme uses a theoretically correct MetropolisâGibbs' sampling technique and an approximate method using numerous models sampled by a global optimization algorithm called very fast simulated annealing (VFSA). VFSA has recently been found to be computationally efficient in several geophysical parameter estimation problems. Unlike conventional simulated annealing (SA), in VFSA the perturbations are generated from the model parameters according to a Cauchy-like distribution whose shape changes with each iteration. This results in an algorithm that converges much faster than a standard SA. In the course of finding the optimal solution, VFSA samples several models from the search space. All these models can be used to obtain estimates of uncertainty in the derived solution. This method makes no assumptions about the shape of an a posteriori probability density function in the model space. Here, we carry out a VFSA-based sensitivity analysis with several synthetic and field sounding data sets for resistivity and IP. The resolution capability of the VFSA algorithm as seen from the sensitivity analysis is satisfactory. The interpretation of VES and IP sounding data by VFSA, incorporating resolution, sensitivity and uncertainty of layer parameters, would generally be more useful than the conventional best-fit techniques.
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.
Pecher, I. A., W. S. Holbrook, M. K. Sen, D. Lizarralde, W. T. Wood, D. R. Hutchison, W. P. Dillon, H. Hoskins, and R. A. Stephen, Seismic anisotropy in gas-hydrate- and gas-bearing sediments on the Blake Ridge, from a walkaway vertical seismic profile, Geophys. Res. Lett., 30, 1733, 2003, 9 citations, doi:10.1029/2003GL017477, #1641
We present results from an analysis of anisotropy in marine sediments using walkaway vertical seismic profiles from the Blake Ridge, offshore South Carolina. We encountered transverse isotropy (TI) with a vertical symmetry axis in a gas-hydrate-bearing unit of clay and claystone with Thomsen parameters ε = 0.05 ± 0.02 and δ = 0.04 ± 0.06. TI increased to ε = 0.16 ± 0.04 and δ = 0.19 ± 0.12 in the underlying gas zone. Rock physics modeling suggests that the observed TI is caused by a partial alignment of clay particles rather than high-velocity gas-hydrate veins. Similarly, the increase of TI in the gas zone is not caused by thin low-velocity gas layers but rather, we speculate, by the sharp contrast between seismic properties of an anisotropic sediment frame and elongated gas-bearing pore voids. Our results underscore the significance of anisotropy for integrating near-vertical and wide-angle seismic data.
Pulliam, J., and M. K. Sen, Assessing uncertainties in waveform modeling of the crust and upper mantle, in Proc. 25th Seismic Res, Rev: Nuclear Explosion Monitoring: Building the Knowledge Base, Tucson, AZ, 2003, #1976
Sen, M. K., and A. Mukherjee, τ-p analysis in transversely isotropic media, Geophys. J. Int., 154, 647-658, 2003, 15 citations, doi:10.1046/j.1365-246X.2003.01997.x, #1638
Incorporating the effects of anisotropy in seismic processing is an active area of research. Here we describe a method for the analysis of P-wave seismic reflection data in the plane-wave domain for transversely isotropic media with a vertical axis of symmetry. The generalization of τâp methods to include transverse isotropy is straightforward in that we simply replace the isotropic vertical slowness term with one appropriate for anisotropy. Exact expressions for the vertical slowness do exist; P and SV require four parameters and SH require two. The P-wave traveltime trajectories are not sensitive to all four parameters and can be expressed as a function of two parameters only. We have developed a series representation for P-wave vertical slowness in terms of increasing powers of horizontal slowness; it is also convenient for interactive moveout analysis. The two-term expression derived from this series is adequate for the analysis of data from weakly anisotropic media and brings out the physics of wave propagation clearly. Thus our development is analogous to the non-hyperbolic normal moveout expression used in the offset-time domain. The τâp analysis enables us to derive interval anisotropy parameters directly without a two-point ray tracing, which can be used in further processing. We have developed a computer algorithm to iteratively and interactively fit the τâp trajectories to derive anisotropy parameters. The technique is demonstrated with applications to synthetic and field reflection seismic data from the Gulf of Mexico.
Sen, M. K., and I. G. Roy, Computation of differential seismograms and iteration adaptive regularization in prestack waveform inversion, Geophysics, 68, 2026-2039, 2003, 12 citations, doi:10.1190/1.1635056, #1639
Seismic waveform inversion is a highly challenging task. Nonlinearity, nonuniqueness, and robustness issues tend to make the problem computationally intractable. We have developed a simple regularized Gauss-Newtonâtype algorithm for the inversion of seismic data that addresses several of these issues. The salient features of our algorithm include an efficient approach to sensitivity computation, a strategy for band-limiting the Jacobian matrix, and a novel approach to computing regularization weight that is iteration adaptive. In this paper, we first review various forward modeling and differential seismogram computation algorithms and then evaluate different strategies for choosing the regularization weight. Under the assumption of locally 1D earth models, we design an efficient algorithm by rearranging recursion formula in the reflection matrix approach to compute plane wave seismograms and the Fréchet derivative matrix as a by-product of forward modeling. We then demonstrate that in a gradient-descentâtype optimization scheme, regularization is critical for obtaining stable and geologically realistic solutions. Although, in most applications, the regularization weight (relative importance between data and model misfit) is chosen in an ad-hoc manner; the robustness in model estimation and computational stability improve significantly by allowing adaptivity in the choice of the regularization weight in each iterative step. We evaluate performances of several methods, namely, an L-curve approach, generalized cross-validation technique, and methods based on a discrepancy principle with application to field ocean-bottom-cable data, and we propose a new hybrid approach in computing iteration adaptive regularization weight for prestack inversion.
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.
Pulliam, J., M. K. Sen, C. Frohlich, and S. P. Grand, Crustal Structure from waveform inversion of shear-coupled PL, in Proceedings of the 23rd Annual Seismic Research Symposium on Monitoring a Comprehensive Test Ban Treaty, Jackson, WY, 110-119, 2001, #1582
Sen, M. K., Pre-stack waveform inversion of plane wave seismograms: Isotropy to transverse isotropy, in Recent Advances in Exploration Geophysics, Canadian Society of Exploration Geophysicists, 85-96, 2001, #1578
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.
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., and M. K. Sen, Geophysical evidence for dewatering and deformation processes in the ODP Leg 170 area offshore Costa Rica, Earth Planet. Sci. Lett., 178, 125-138, 2000, 17 citations, doi:10.1016/S0012-821X(00)00069-8, #1508
We use a combination of borehole data from Ocean Drilling Program (ODP) Leg 170 and multichannel seismic reflection (MCS) data to quantify thickness changes in underthrust sediments away from the boreholes. Sediments thrust beneath the upper plate at convergent margins may be more rapidly loaded than in any other environment. Depending on the porosity and permeability of the available fluid pathways, these sediments can compact and dewater very rapidly, as observed in this area offshore the Nicoya Peninsula, Costa Rica. Rapid thinning and dewatering was previously interpreted in this area from MCS data, but the lack of velocity data in this deep-water environment caused ambiguity in the estimates of thickness change. We employ a non-linear inversion technique using detailed density data, primarily logs and some laboratory measurements and coincident MCS data to create 1D synthetic seismograms and detailed velocity functions at three ODP drill sites. Because only a small part of one hole was logged with a sonic tool and the laboratory measurements significantly underestimate in situ velocities, these results provide the most accurate estimate of the velocity profiles. We used these velocity functions to depth-migrate seven MCS lines in the vicinity of the trench and lower slope spanning a distance of 9 km along strike. Analysis of the depth-migrated images shows that there is significant variation along strike in how the underthrust section compacts, which appears to be related to the distribution of normal faults on the Cocos Plate. We interpret that preferentially rapid dewatering in the upper part of the underthrust section may lead to deformation below the original decollement and detachments at deeper stratigraphic levels.
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
Sen, M. K., (Book review): Seismic Inversion and Deconvolution: Dual-sensor Technology, by E. A. Robinson, Eos, Trans. Amer. Geophys. Un., 81, 368, 2000, doi:10.1029/00EO00275, #1521
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
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.
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
Pulliam, J., and M. K. Sen, Seismic anisotropy in the core-mantle transition zone, Geophys. J. Int., 135, 113-128, 1998, 32 citations, doi:10.1046/j.1365-246X.1998.00612.x, #1391
Split S waves observed at Hockley, Texas from events in the TongaâFiji region of the southwest Pacific show predominantly vertically polarized shear-wave (SV ) energy arriving earlier than horizontally polarized (SH) energy for rays propagating horizontally through D". After corrections are made for the effects of upper-mantle anisotropy beneath Hockley, a time lag of 1.5 to 2.0 s remains for the furthest events (93.9°â100.6° ), while the time lags of the nearer observations (90.5°â92.9° ) nearly disappear. At closer distances, the S waves from these same events do not penetrate as deeply into the lower mantle, and are not split. These observations suggest that a patch of D" beneath the central Pacific is anisotropic, while the mantle immediately above the patch is isotropic. The thickness of the anisotropic zone appears to be of the order of 100â200 km.
Observations of shear-wave splitting have previously been made for paths that traverse D" under the Caribbean and under Alaska. SH leads SV , the reverse of the Hockley observations, but in these areas the fact that SV leads SH in the HKT data shown here suggests a different sort of anisotropy under the central Pacific from that under Alaska and the Caribbean. The case of SH travelling faster than SV is consistent with transverse isotropy with a vertical axis of symmetry (VTI) and does not require variations with azimuth. The case of SV leading SH is consistent with transverse isotropy with a horizontal axis of symmetry (HTI), an azimuthally anisotropic medium, and with a VTI medium formed by a hexagonal crystal. Given that (Mg,Fe)SiO3 perovskite appears unlikely to form anisotropic fabrics on a large scale, the presence of anisotropy may point to chemical heterogeneity in the lowermost mantle, possibly due to mantleâcore interactions.
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, 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
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.
Sen, M. K., and G. Y. Xia, Parameter estimation in anisotropic media, , Proc., 67th SEG Annual Meeting, Dallas, TX, 1559-1562, 1997, #1364
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
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.
Matzel, E., M. K. Sen, and S. P. Grand, Evidence for anisotropy in the deep mantle beneath Alaska, Geophys. Res. Lett., 23, 2417-2420, 1996, 53 citations, #1227
We consider the possibility that the velocity structure of D″ is anisotropic. The data we examined consist of seismograms from 9 deep Japanese earthquakes recorded at WWSSN receiver stations in North America. The source‐receiver combinations span distances of 70°â106° with associated S waves passing through D″ beneath Alaska. Differential travel times of the S, Scd, ScS and SKS phases are used to constrain the velocity structure in D″. Shear waves refracted by D″ are observed beyond 72.2° and provide a sensitive measurement of the velocity structure in D″. Beyond 93°, the vertically polarized (SV) and horizontally polarized (SH) shear waves often appear distinctly split, although, at distances less than 89° the components are more nearly synchronous. Near 94°, SH occurs as a double arrival. SV in this range, however, remains a single arrival roughly synchronous with the second SH arrival. We have been unable to reproduce these effects in isotropic model synthetics. Synthetics for transversely isotropic models have been computed that do match these waveforms. The anisotropy was constrained to be only within D″, with a vertical symmetry axis. We conclude that these observations may be explained by an anisotropic D″ layer. The D″ discontinuity may be due to a transition to anisotropic mantle a few hundred kilometers above the core‐mantle boundary.
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.
Lindwall, D. A., M. K. Sen, and J. F. Gettrust, Detection of high shear wave velocities in marine sediment by inversion with simulated annealing, in Full Field Ocean Methods in Ocean and Seismo-acoustics, edited by O. Diastchok, A. Caiti, P. Gersoft, and H. Schmidt, Kluwer Academics, Netherlands, 383-388, 1995, #1181
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
Wood, W. T., and M. K. Sen, Non-linear inversion of seismic data by successive approximation of model parameter probability distribution functions, in Full Field Inversion Methods in Ocean and Seismo-acoustics, edited by O. Diastchok, A. Caiti, P. Gersoft, and H. Schmidt, Kluwer Academics, Netherlands, 147-152, 1995, #1180
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
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.
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.
Stoffa, P. L., and M. K. Sen, Seismic waveform inversion using global optimization, J. Seismic Explor., 1, 9-27, 1992, #920
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.
Sen, M. K., and L. N. Frazer, Multifold phase space path integral synthetic seismograms, Geophys. J. Int., 104, 479-487, 1991, 8 citations, #827
Sen, M. K., Modeling of wave propagation in northern Los Angeles Basin, Bull. Seismol. Soc. Amer., 81, 751-768, 1991, 2 citations, #856
Somerville, P., M. K. Sen, and B. Cohee, Simulation of strong ground motions recorded during the 1985 Michoacan, Mexico and Valparaiso, Chile earthquakes, Bull. Seismol. Soc. Amer., 81, 1-27, 1991, 51 citations, #912
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.
Sen, M. K., Deep structural complexity and site response in Los Angeles Basin, Proc., 4th U.S. National Symp. on Earthquake Engineering, 2, 545-554, 1990, #800
Tajima, F., and M. K. Sen, Implication from the aftershocks of the 1989 Loma Prieta earthquake, Geophys. Res. Lett., 17, 1421-1424, 1990, 2 citations, #802
We examined the P wave displacement spectra of nine aftershocks (M=2.5∼3.6) of the 1989 Loma Prieta earthquake recorded on PASSCAL/IRIS instruments to search for changes in the spectra that are indicative of heterogeneity in the source region. Temporal variations of spectral characteristics have been observed from the records of two events (#7 with M=2.7 and #8 with M=2.5) that occurred at an almost identical location with a time interval of 16 min and recorded at common stations. The different spectral characteristics can be attributed to the differences in the source time functions of these events. Also, the spectrum of event #2 (M=2.5) that occurred near San Francisco, is substantially richer in high frequency than those from events of a comparable size in the immediate aftershock area. This could be an indication of progressive stress concentration beyond the present aftershock area.