Search

Full-waveform inversion (FWI) is a nonlinear optimization problem that addresses the estimation of subsurface model parameters by matching the predicted to the observed seismograms. We formulate FWI as a constrained optimization problem in which the regularization term is minimized subject to the nonlinear data matching constraints. Unlike standard FWI, which solves this regularized problem with a penalty method, we use an augmented Lagrangian formulation in the framework of the method of multipliers (MM). This leads to a two-step recursive algorithm, which is called multiplier waveform inversion (MWI). First, the primal step solves an unconstrained minimization problem, which is similar to the standard FWI but with a modified data term. These modified data are obtained by adding the Lagrange multipliers to the data residuals of the classical FWI. Then, the dual step updates the Lagrange multipliers with basic gradient ascent steps. In this framework, they reduce to the running sum of the data residuals of the previous iterations. The performance of the overall algorithm is improved by considering that the MM does not require the exact solution of the primal step at each iteration. In fact, convergence is attained when the primal subproblem is solved with one (without an inner loop) iteration of a gradient-based method at each iteration of the MWI. The new algorithm can be easily implemented in existing FWI codes without computational overhead by adding the Lagrange multipliers to the data residuals at each iteration. We find with numerical examples that MWI has a faster convergence speed and much improved stability compared with FWI and can converge to an accurate solution of the inverse problem in the absence of low-frequency data, even with a constant step size.

Reverse time migration (RTM), as a state-of-the-art imaging technique, provides outstanding imaging capabilities due to its use of a full wave equation. Least-squares RTM (LSRTM) seeks the solution of a linearized wave equation via the minimization of a data misfit term; however, the quality of the results decreases when the assumptions of the method are not satisfied. This occurs, for example, when we use an erroneous velocity model or inadequate physics for inverting the data. In such cases, appropriate regularization is required to mitigate these shortcomings and stabilize the LSRTM solution. However, even for structurally simple earth models, the reflectivity images are complicated and may not be explained properly by particular regularization methods such as the Tikhonov, total variation (TV), or sparse regularization. Reflectivity images can be thought of as the difference between two structurally simpler components: a piecewise-constant component (the true squared slowness) and a smooth component (the background model). We have used a combined Tikhonov-TV regularizer to regularize these components separately, leading to an effective regularization for the reflectivity image. Because the background model is known in advance, this combined regularization reduces to a shifted TV regularization for which the associated optimization problem is solved efficiently using a new implementation of the Bregmanized operator splitting algorithm applied to the shifted TV method and the usual TV method. We determine the performance of our method with a set of numerical examples. The results confirm that our shifted regularization increases the robustness of LSRTM and allows us to estimate high-quality reflectivity images and properly update the background velocity model.

Efficient frequency-domain full-waveform inversion (FWI) of long-offset node data can be performed with a few frequencies. The seismic response of these frequencies can be computed with compact finite-difference stencils on regular Cartesian grid with direct or hybrid direct/iterative methods. Compactness, which is necessary to mitigate the fill-in induced by the lower-upper factorization, is implemented with second-order stencils, whereas accuracy is achieved by building a consistent mass matrix and a compound stiffness matrix with optimal weights so that the stencil covers the eight cells surrounding the central point, leading to a 27-point stencil. Classical approaches estimate constant weights by jointly minimizing numerical dispersion in homogeneous media for several numbers of grid points per wavelength ([Formula: see text]). Then, the impedance matrix is built by using these constant weights at each central point covering the heterogeneous subsurface model, leading to nonuniform wavefield accuracy. Instead, we estimate [Formula: see text]-dependent weights once and for all by minimizing dispersion separately for each value of [Formula: see text] in a range found in FWI applications. Then, we build the impedance matrix without computational overhead by selecting at each central point the weights corresponding to the local [Formula: see text], hence leading to a wavelength-adaptive stencil. This separate approach with adaptive weights makes wavefield accuracy uniform. We benchmark wavefields, which are computed in several 3D large-scale subsurface models with a sparse multifrontal direct solver and the nonadaptive/adaptive stencils, against analytical solutions when available and the highly accurate discretization-free convergent Born series method. Each benchmark reveals the higher accuracy of the adaptive stencil relative to the nonadaptive one. In the presence of sharp contrasts, the adaptive stencil also is more accurate than a classical [Formula: see text] finite-different time-domain staggered-grid stencil. The relevance of the adaptive stencil is finally illustrated with a 3D FWI case study in the 3.5–13 Hz frequency band.

International audience; The augmented Lagrangian (AL) method provides a flexible and efficient framework for solving extended-space full-waveform inversion (FWI), a constrained nonlinear optimization problem whereby we seek model parameters and wavefields that minimize the data residuals and satisfy the wave-equation constraint. The AL-based wavefield reconstruction inversion, also known as iteratively refined wavefield reconstruction inversion, extends the search space of FWI in the source dimension and decreases the sensitivity of the inversion to the initial model accuracy. Furthermore, it benefits from the advantages of the alternating direction method of multipliers, such as generality and decomposability for dealing with nondifferentiable regularizers, e.g., total variation regularization, and large-scale problems, respectively. In practice, any extension of the method aiming at improving its convergence and decreasing the number of wave-equation solves would have great importance. To achieve this goal, we recast the method as a general fixed-point iteration problem, which enables us to apply sophisticated acceleration strategies such as Anderson acceleration. The accelerated algorithm stores a predefined number of previous iterates and uses their linear combination together with the current iteration to predict the next iteration. We investigate the performance of our accelerated algorithm on a simple checkerboard model and the benchmark Marmousi II and 2004 BP salt models through numerical examples. These numerical results confirm the effectiveness of our algorithm in terms of convergence rate and the quality of the final estimated model.

Different versions of the Radon transform (RT) are widely used in seismic data processing to focus the recorded seismic events. Multiple separation, data interpolation, and noise attenuation are some of the RT applications in seismic processing workflows. Unfortunately, conventional RT methods cannot focus events perfectly in the RT domain. This problem arises due to the blurring effects of the source wavelet and the nonstationary nature of the seismic data. Sometimes, the distortion results in a big difference between the original data and its inverse transform. We have developed a nonstationary deconvolutive RT (DecRT) to handle these two issues. Our algorithm takes advantage of a nonstationary convolution technique that builds on the concept of block convolution and the overlap method, in which the convolution operation is defined separately for overlapping blocks. Therefore, it allows the Radon basis function to take arbitrary shapes in the time and space directions. In addition, we introduce a nonstationary wavelet estimation method to determine time-space-varying wavelets. The wavelets and the Radon panel are estimated simultaneously and in an alternating way. Numerical examples demonstrate that our nonstationary DecRT method can significantly improve the sparsity of Radon panels. Hence, inverse RT does not suffer from the distortion caused by unfocused seismic events.

Efficient frequency-domain full-waveform inversion (FD-FWI) of wide-aperture data is designed by limiting inversion to a few frequencies and by solving the Helmholtz equation with a direct solver to process multiple sources efficiently. Some variants of FD-FWI, which process the wave equation as a weak constraint, have been proposed to increase the computational efficiency or extend the search space. Among them, the contrast-source reconstruction inversion (CSRI) reparameterizes FD-FWI in terms of contrast sources (CS) and contrasts and updates them in an alternating mode. This reparameterization allows for one lower-upper (LU) decomposition of the Helmholtz operator to be performed per frequency inversion hence further improving the computational efficiency of FD-FWI. However, iteratively refined wavefield reconstruction inversion (IR-WRI) relies on the alternating-direction method of multipliers to extend the search space by matching the data from the early iterations via an aggressive relaxation of the wave equation while satisfying it at the convergence point thanks to the defect correction performed by the Lagrange multipliers. In contrast to CSRI, IR-WRI requires redoing one LU decomposition when the subsurface model is updated. In both methods, the CSs or the wavefields are computed by solving in a least-squares sense an overdetermined linear system gathering an observation equation and a wave equation. A drawback of CSRI is that CSs are estimated approximately with one iteration of a conjugate gradient method, whereas the wavefields are reconstructed exactly by IR-WRI with a Gauss-Newton method. We have combined the benefits of CSRI and IR-WRI to decrease the number of LU decomposition during IR-WRI with a fixed-point (FP) algorithm while preserving its search space extension capability. Application on the 2D complex Marmousi and the BP salt models shows that our FP-based IR-WRI manages to reconstruct these models as accurately as the classic IR-WRI while reducing the number of LU factorizations considerably. A theoretical complexity analysis and a recent application of 3D FD-FWI based upon direct solver suggest that the FP algorithm should reduce the cost of IR-WRI by a factor of approximately 2 and 10 for 3D dense ocean bottom cable and sparse ocean bottom node acquisitions, respectively.

Dispersion-curve inversion of Rayleigh waves to infer subsurface shear-wave velocity is a long-standing problem in seismology. Due to nonlinearity and ill-posedness, sophisticated regularization techniques are required to solve the problem for a stable velocity model. We have formulated the problem as a minimization problem with nonlinear operator constraint and then solve it by using an inexact augmented Lagrangian method, taking advantage of the Haney-Tsai Dix-type relation (a global linear approximation of the nonlinear forward operator). This replaces the original regularized nonlinear problem with iterative minimization of a more tractable regularized linear problem followed by a nonlinear update of the phase velocity (data) in which the update can be performed accurately with any forward modeling engine, for example, the finite-element method. The algorithm allows discretizing the medium with thin layers (for the finite-element method) and thus omitting the layer thicknesses from the unknowns and also allows incorporating arbitrary regularizations to shape the desired velocity model. In this research, we use total variation regularization to retrieve the shear-wave velocity model. We use two synthetic and two real data examples to illustrate the performance of the inversion algorithm with total variation regularization. We find that the method is fast and stable, and it converges to the solution of the original nonlinear problem.

Wavefield reconstruction inversion (WRI) extends the search space of Full Waveform Inversion (FWI) by allowing for wave equation errors during wavefield reconstruction to match the data from the first iteration. Then, the wavespeeds are updated from the wavefields by minimizing the source residuals. Performing these two tasks in alternating mode breaks down the nonlinear FWI as a sequence of two linear subproblems, relaying on the bilinearity of the wave equation. We solve this biconvex optimization with the alternating-direction method of multipliers (ADMM) to cancel out efficiently the data and source residuals in iterations and stabilize the parameter estimation with appropriate regularizations.Here, we extend WRI to viscoacoustic media for attenuation imaging. Attenuation reconstruction is challenging because of the small imprint of attenuation in the data and the crosstalks with velocities. To address these issues, we recast the multivariate viscoacoustic WRI as a triconvex optimization and update wavefields, squared slowness, and attenuation factor in alternating mode at each WRI iteration. This requires to linearize the attenuation-estimation subproblem via an approximated trilinear viscoacoustic wave equation. The iterative defect correction embedded in ADMM corrects the errors generated by this linearization, while the operator splitting allows us to tailor ℓ1 regularization to each parameter class. A toy numerical example shows that these strategies mitigate crosstalk artifacts and noise from the attenuation reconstruction. A more realistic example representative of the North sea further shows the ability of the method to jointly reconstruct the wavespeed and the attenuation starting from a very crude attenuation-free initial model, the moderate strength of the crosstalk artifacts and the sensitivity of the multi-parameter reconstruction to noise.

Conventional amplitude variation with offset (AVO) inversion analysis uses the Zoeppritz equations, which are based on a plane-wave approximation. However, because real seismic data are created by point sources, wave reflections are better modeled by spherical waves than by plane waves. Indeed, spherical reflection coefficients deviate from planar reflection coefficients near the critical and postcritical angles, which implies that the Zoeppritz equations are not applicable for angles close to critical reflection in AVO analysis. Elastic finite-difference simulations provide a solution to the limitations of the Zoeppritz approximation because they can handle near- and postcritical reflections. We have used a coupled acoustic-elastic local solver that approximates the wavefield with high accuracy within a locally perturbed elastic subdomain of the acoustic full domain. Using this acoustic-elastic local solver, the local wavefield generation and inversion are much faster than performing a full-domain elastic inversion. We use this technique to model wavefields and to demonstrate that the amplitude from within the local domain can be used as a constraint in the inversion to recover elastic material properties. Then, we focus on understanding how much the amplitude and phase contribute to the reconstruction accuracy of the elastic material parameters ([Formula: see text], [Formula: see text], and [Formula: see text]). Our results suggest that the combination of amplitude and phase in the inversion helps with the convergence. Finally, we analyze elastic parameter trade-offs in AVO inversion, from which we find that to recover accurate P-wave velocities we should invert for [Formula: see text] and [Formula: see text] simultaneously with fixed density.

We have developed a simple method to halve the memory required to store the forward wavefield for the imaging condition in adjoint-state methods such as reverse time migration and full-waveform inversion. It stores the wavefield at only half the Nyquist rate, and it uses the wave equation to calculate the second time derivative, allowing approximate reconstruction of the forward wavefield at the required Nyquist rate. We have determined that the method produces an image with the Marmousi model that is not visibly different compared to a traditional (full Nyquist) implementation.

The traditional hyperbolic Radon transform (RT) decomposes seismic data into a sum of constant amplitude basis functions. This limits the performance of the transform when dealing with real data in which the reflection amplitudes include the amplitude variation with offset (AVO) variations. We adopted the Shuey-Radon transform as a combination of the RT and Shuey’s approximation of reflectivity to accurately model reflections including AVO effects. The new transform splits the seismic gather into three Radon panels: The first models the reflections at zero offset, and the other two panels add capability to model the AVO gradient and curvature. There are two main advantages of the Shuey-Radon transform over similar algorithms, which are based on a polynomial expansion of the AVO response. (1) It is able to model reflections more accurately. This leads to more focused coefficients in the transform domain and hence provides more accurate processing results. (2) Unlike polynomial-based approaches, the coefficients of the Shuey-Radon transform are directly connected to the classic AVO parameters (intercept, gradient, and curvature). Therefore, the resulting coefficients can further be used for interpretation purposes. The solution of the new transform is defined via an underdetermined linear system of equations. It is formulated as a sparsity-promoting optimization, and it is solved efficiently using an orthogonal matching pursuit algorithm. Applications to different numerical experiments indicate that the Shuey-Radon transform outperforms the polynomial and conventional RTs.

Given the ill-conditioned nature of Dix inversion, the resultant Dix interval-velocity field is often unrealistic, noisy, and highly dependent on the quality of the provided root-mean-square velocities. While the classic least-squares regularization techniques, e.g., various forms of Tikhonov regularization, lead to somewhat suboptimal stability, we formulated the Dix inversion as a new constrained optimization problem. This enables one to incorporate prior knowledge as soft and/or hard bounds for the optimization, effectively treating it as a denoising problem. The solution to the problem is achieved by a bound-constrained total variation (TV) regularization. TV regularization has the advantage of being able to recover the discontinuities in the model, but it often comes with a large memory and compute requirements. Therefore, we have developed a simple and memory-efficient algorithm using iterative refinement strategy. The quality of the new algorithm is also cross-examined against different strategies, which are currently used in practice. Overall, we observe that the proposed method outperforms classic Dix inversion methods on the synthetic and real data examples.

We have developed an adaptive singular spectrum analysis (ASSA) method for seismic data denoising and interpolation purposes. Our algorithm iteratively updates the singular-value decomposition (SVD) of current spatial patches using the most recently added spatial sample. The method reduces the computational cost of classic singular spectrum analysis (SSA) by requiring QR decompositions on smaller matrices rather than the factorization of the entire Hankel matrix of the data. A comparison between results obtained by the ASSA and SSA methods, in which the SVD applies to all of the traces at once, proves that the ASSA method is a valid way to cope with spatially varying dips. In addition, a comparison of the ASSA method with the windowed SSA method indicates gains in efficiency and accuracy. Synthetic and real data examples illustrate the effectiveness of our method.

Group velocity is an important characteristic of surface wave that is defined as the velocity of an envelope of frequencies. Although many studies have shown the promises of analyzing the group velocity to obtain subsurface S-wave velocity, the estimation of the group velocity is not straightforward due to the uncertainties of selecting an optimum envelope of frequencies. Conventional transformations or filtering algorithms used to define an optimum envelope usually give reasonable results just for a narrow frequency or velocity range. We introduced a new approach for the estimation of the group velocity using the sparse S transform (SST) and sparse linear Radon transform (SLRT). In SST, the width of the Gaussian window is optimally calculated by energy concentration to eliminate energy smearing in the time-frequency (TF) domain, and then the sparsity is applied to enhance the TF resolution. Compared with conventional methods for the estimation of the group velocity based on the generalized S transforms, SST does not require any adjustment to the Gaussian window and yields accurate estimates of the group velocity. We apply SST to each seismic trace of a seismic shot record to obtain a 3D cube of frequency, time, and offset. For any frequency, we obtain a common frequency gather of time and offset to which we apply SLRT to obtain the group velocity of the surface wave. Our approach is robust at calculating high-resolution distinguishable dispersion curves of the group velocity in particular when data are extremely sparse.

Seismic deconvolution used for improving the bandwidth of data is inherently nonstationary, mixed phase, and blind. Due to some restricting assumptions imposed by conventional deconvolution methods, they are either stationary or semiblind. A fully nonstationary blind deconvolution method is proposed that is able to simultaneously take into account different sources of nonstationarity and to improve the bandwidth of highly nonstationary seismic data in a fully blind manner. Based on the concept of block convolution and the overlap method, the convolutional model of seismic data is generalized to consider nonstationary cases and to model nonstationary data. This generalized convolutional model is then used for nonstationary blind deconvolution, in which the statistical characteristics of the wavelets are allowed to arbitrarily change in the vertical and horizontal directions. Given a nonstationary seismic record, several time-space-varying wavelets are simultaneously determined with the reflectivity model in an alternating direction algorithm using a variational approach. Numerical tests are presented showing the high performance of our nonstationary blind deconvolution for improving the temporal resolution of data in comparison with their stationary counterparts. The results indicate that in comparison with patched deconvolution, our nonstationary method is more robust and stable for different window sizes and it produces better results with a higher signal-to-noise ratio.

Seismic velocity analysis is one of the most crucial and, at the same time, the most laborious tasks during seismic data processing. This becomes even more difficult and time-consuming when nonhyperbolicity has to be considered in the velocity analysis. Nonhyperbolic velocity analysis provides very useful information during the processing and interpretation of seismic data. The most common approach for considering anisotropy during velocity analysis is to describe the moveout based on a nonhyperbolic equation. The nonhyperbolic moveout equation in vertically transverse isotropic (VTI) media is defined by two parameters: normal moveout (NMO) velocity [Formula: see text] and anellipticity [Formula: see text] (or horizontal velocity [Formula: see text]). We have developed a new approach based on polynomial chaos (PC) expansion for automating nonhyperbolic velocity analysis of common-midpoint (CMP) data in VTI media. For this purpose, we use the PC expansion to approximate the nonhyperbolic semblance function with a very fast-to-simulate function in terms of [Formula: see text] and [Formula: see text]. Then, using particle swarm optimization, we stochastically look for the optimum NMO and horizontal velocities that provide the maximum semblance. In contrary to common approaches for nonhyperbolic velocity analysis in which the two parameters are estimated iteratively in an alternating fashion, we find [Formula: see text] and [Formula: see text] simultaneously. This approach is tested on various data including a simple convolutional model, an anisotropic benchmark model, and a real data set. In all cases, the new method provided acceptable results. Reflections in the CMP corrected using the optimum velocities are properly flattened, and almost no residual moveout is observed.

A seismic trace corresponding to an anelastic layered earth model can be sparsely represented by a structured dictionary of properly attenuated wavelets, through a nonstationary sparse deconvolution. The sparseness of the coefficients (reflectivity series), however, considerably decreases when the wavelets are incorrectly modeled due to an incorrect [Formula: see text] model. Mathematically, the wavelets, as the elements of the dictionary, are nonlinearly related to the earth quality factor [Formula: see text] (the inverse of the attenuation coefficient). A parametric dictionary-learning strategy enables interval-[Formula: see text] estimation and compensation by training the dictionary atoms from the input trace adaptively to provide a sparse representation of it. We assumed a piecewise [Formula: see text] model by dividing the dictionary elements into several groups, each containing several wavelets whose temporal supports are close to each other and can be described by a single [Formula: see text]-value. The dictionary is learned iteratively where at each iteration only one group of the wavelets is optimized by searching for the corresponding optimum [Formula: see text]-value, leading to an iterative construction of the [Formula: see text] model. The main advantages of our method for interval-[Formula: see text] estimation are its stability because it performs in a forward-modeling manner and its accuracy because of the resolved interferences by the sparsity constraint. Our method is tested on synthetic and field data sets, and the results that we obtained demonstrate the stability and accuracy of the method for interval-[Formula: see text] estimation and compensation.

The inversion of prestack seismic data using amplitude variation with offset (AVO) has received increased attention in the past few decades because of its key role in estimating reservoir properties. AVO is mainly governed by the Zoeppritz equations, but traditional inversion techniques are based on various linear or quasilinear approximations to these nonlinear equations. We have developed an efficient algorithm for nonlinear AVO inversion of precritical reflections using the exact Zoeppritz equations in multichannel and multi-interface form for simultaneous estimation of the P-wave velocity, S-wave velocity, and density. The total variation constraint is used to overcome the ill-posedness while solving the forward nonlinear model and to preserve the sharpness of the interfaces in the parameter space. The optimization is based on a combination of Levenberg’s algorithm and the split Bregman iterative scheme, in which we have to refine the data and model parameters at each iteration. We refine the data via the original nonlinear equations, but we use the traditional cost-effective linearized AVO inversion to construct the Jacobian matrix and update the model. Numerical experiments show that this new iterative procedure is convergent and converges to a solution of the nonlinear problem. We determine the performance and optimality of our nonlinear inversion algorithm with various simulated and field seismic data sets.

Various effects cause a seismic trace to consist of inhomogeneous (phase-shifted) wavelets, which can be modeled by constant phase rotations of a base wavelet. Analytic deconvolution, an analytic signal formulation of the conventional statistical deconvolution, has been proposed for deconvolution of such traces, but it seems to be ineffective in the presence of wavelet interferences and noise because it has no control over the estimated phase. Morphological deconvolution (MD) allows full control over the estimated phase and reflectivity; hence, it is more effective than analytic deconvolution in deconvolving traces with nonstationary phase variations. MD replaces the convolution model of a seismic trace with the sum of a set of convolutions with phase-rotated wavelets. The solution to the resulting under-determined system of equations provides the reflectivity as a function of the time and phase shifts, and it is found by a sparse plus low-rank optimization. The sparsity constraint forces decomposition of the trace into the least number of wavelets, whereas the low-rank constraint allows for regularizing the phase function. MD is easily applicable in multichannel form (MMD), allowing regularization of the estimated phase in the time and space directions. The efficacy and robustness of the MD and MMD methods are determined on several synthetic and field data sets in the presence of interference and noise.

There are lots of geophysical problems that include computationally expensive functions (forward models). Polynomial chaos (PC) expansion aims to approximate such an expensive equation or system with a polynomial expansion on the basis of orthogonal polynomials. Evaluation of this expansion is extremely fast because it is a polynomial function. This property of the PC expansion is of great importance for stochastic problems, in which an expensive function needs to be evaluated thousands of times. We have developed PC expansion as a novel technique to solve nonlinear geophysical problems. To better evaluate the methodology, we use PC expansion for automating the velocity analysis. For this purpose, we define the optimally picked velocity model as an optimizer of a variational integral in a semblance field. However, because computation of a variational integral with respect to a given velocity model is rather expensive, it is impossible to use stochastic methods to search for the optimal velocity model. Thus, we replace the variational integral with its PC expansion, in which computation of the new function is extremely faster than the original one. This makes it possible to perturb thousands of velocity models in a matter of seconds. We use particle swarm optimization as the stochastic optimization method to find the optimum velocity model. The methodology is tested on synthetic and field data, and in both cases, reasonable results are achieved in a rather short time.

The Radon transform (RT) plays an important role in seismic data processing for its ability to focus seismic events in the transform domain. The band-limited nature of seismic events due to the blurring effects of the source wavelet, however, causes a decrease in the temporal resolution of the transform. We have developed the deconvolutive RT (DecRT) as a generalization of conventional RT and to increase the temporal resolution. Unlike the conventional counterpart, the new basis functions can take an arbitrary shape in the time direction. This method is thus proposed to adaptively infer the temporal wave shape from the input data while finding a sparse representation of it. The new transform significantly improves the sparsity and thus the temporal resolution of the resulting seismic data. The applicability of the hyperbolic DecRT in seismic data processing is demonstrated for random noise attenuation, primary and multiple separation, high-quality stacking, and automatic velocity model building. The results obtained on synthetic and field data sets confirm the effectiveness of the method in improving the time and slowness/curvature resolutions compared with conventional transforms, which leads to improved seismic processing results in the deconvolutive Radon domains.

The inversion of seismic reflection data for acoustic impedance (AI) is a common and accepted method for the interpretation of poststack seismic data. The original mathematical problem is nonlinear due to the nonlinearity of relation between the AI and the earth reflectivity series and also the insufficiency of information about the source wavelet. Furthermore, the problem is ill posed due to the lack of low and high frequencies in the data. We have developed a multichannel blind inversion and solved it for obtaining the AI model and the wavelet directly from seismic reflection data. We found a solution to the overall problem by alternating between two subproblems, corresponding to the AI and wavelet recovery. Having an estimation of the wavelet, the algorithm directly inverts multichannel data for a high-resolution AI model, having blocky structures in the sense of total variation (TV) regularization, while satisfying a priori low-frequency information. The advantages of the split Bregman technique and the discrete cosine transform are used to build a fast and efficient algorithm for solving the nonlinear impedance inversion with TV regularization. Having an estimate of the AI model, the wavelet is updated by restricting it to have a sparse representation in a wavelet transformation domain while predicting the observed data. Numerical tests using simulated 2D data obtained from the benchmark Marmousi model and also 2D field data confirm that the proposed algorithm stably generates accurate estimates of complex AI models and complicated wavelets, simultaneously.

Gabor deconvolution consists of approximating a nonstationary problem as several stationary subproblems via the Gabor frame, solving each subproblem independently, and then recombining/projecting the subsolutions into an approximate solution to the original nonstationary problem. The approximations, however, cause inherent instability due to systematic errors that prevent the algorithm from converging to the true solution even for noise-free cases. Furthermore, the method will be time consuming when a nonlinear optimization is used for shaping the reflectivity structure. An alternative projected Gabor deconvolution is considered, which is based on the Gabor expansion. However, in contrast to the conventional form, in the new method, first the problem is projected into the time domain, and then an inversion is performed to obtain the final reflectivity. Compared with the Gabor deconvolution, the projected alternative (1) exhibits an improved convergence property, (2) is more efficient for sparse deconvolution because only a single optimization is required to be solved, and (3) is more flexible for incorporating prior information about noise and reflectivity structure via a least-squares method. Numerical tests using simulated and field data are presented showing that the new method generates more accurate and stable reflectivity models compared with the Gabor deconvolution.

The analysis of acoustic impedance (AI) allows for the mapping of seismic reflection data to lithology, and hence it plays an important role in the interpretation of poststack seismic data. The AI is obtainable from the inversion of the earth reflectivity series. Efficient deconvolution methods have been developed for recovering the reflectivity series from band-limited poststack data, which are multiple free, zero offset, and migrated. However, calculation of the AI from the reflectivity, when considering the spatial correlation of the impedance parameters, demands the solution of a constrained nonlinear inverse problem. Two efficient algorithms are proposed for solving the nonlinear impedance problem in multichannel form with the total-variation (TV) constraint to recover impedance maps with blocky structures. The first uses the continuous earth model for reflectivity, which allows linearizing the problem in the logarithm domain. The second uses the discrete (layered) earth model for reflectivity. In this case, a nonlinear problem is solved in the original domain. The main properties of the algorithms, which are built based on the split-Bregman technique and the discrete cosine transform are as follows: (1) They are multichannel inversion procedures, so they respect the spatial and temporal correlation of the data, (2) they impose the TV constraint on the parameters to generate blocky structures, (3) they are fast such that a large 3D impedance model with a blocky structure can be recovered easily on a personal computer, and (4) the computational complexities of both algorithms are the same. Numerical tests using simulated 2D data obtained from a benchmark Marmousi model and also 2D and 3D field data confirmed that the proposed algorithm generates more accurate with higher resolution impedance models compared with the conventional methods.

The most advanced data processing/imaging algorithms aim at optimizing a suitably defined cost function consisting of a data misfit term, posed in the data domain, and a regularization term. For strong phase errors in the data, however, such algorithms fail to achieve an acceptable solution because they force us to match the incorrect phase information. To remedy this deficiency, we replaced the time domain misfit term by a Fourier magnitude fidelity term to recover the desired signal from just the amplitude spectrum of the observations. Such a problem is known as phase retrieval, and it is a subject of interest for those who cannot measure the phase information about the system under study. The proposed phase retrieval problem is solved by a fast iterative shrinkage/thresholding algorithm, and it is used to tackle two common phase problems arising in seismic processing: data recovery in the presence of residual static shifts and deconvolution with a missing wavelet phase. Under sparsity constraints, both problems can be solved using just the amplitude spectrum of the observed data. In both cases, using computer-simulated data and field data, the regularized phase retrieval algorithm was able to obtain better results than the conventional methods in which the misfit term is posed in the time domain.

The variation of frequency content of a seismic trace with time carries information about the properties of the subsurface reflectivity sequence. Time-frequency (TF) analysis is a significant tool to extract such information for seismostratigraphic interpretation purposes. However, several TF transforms have been reported in the literature; higher resolution and sensitivity to local changes of the signal have always mattered. We have developed an adaptive high-resolution TF transform that is performed in two sequential steps: First, the window length is adaptively determined for each sample of the signal such that it leads to maximum compactness of energy in the resulting TF plane. Second, the generated nonstationary windows are used to inversely decompose the signal under study via a convex constrained sparse optimization, where a mixed norm of the TF coefficients is minimized subject to invertibility of the transform. Later on, the optimized transform is used as an efficient tool for seismic data analysis such as thin-bed characterization and thin-bedded gas reservoir detection. In the case of gas reservoir detection, based on amplitude versus offset analysis in the TF domain, a simple new method called the difference section was evaluated. The results of various numerical examples from synthetic and field data revealed a remarkable performance of the proposed method compared with the state-of-the-art TF transforms.

Residual statics estimation in complex areas is one of the main challenging problems in seismic data processing. It is well known that the result of this processing step has a profound effect on the quality of final reconstructed image. A novel method is presented to compensate for surface-consistent residual static corrections based on sparsity maximization, which has proved to be a powerful tool in the analysis and processing of signals and related problems. The method is based on the hypothesis that residual static time shift represents itself by noise-like features in the Fourier or curvelet domain. Residual time shift corrections are then retrieved by optimizing the sparsity in these domains. Here, the statics model is considered as a maximizer of [Formula: see text]-norm ([Formula: see text]) of the data coefficients in the sparse domain, and a fast and efficient algorithm is presented to iteratively solve the corresponding nonlinear optimization problem. Applications on synthetic and real data show very high performance of the presented algorithm.

Peeter Akerberg Recent Ph.D. graduates are invited to submit abstracts of their dissertations by completing and submitting the submission form 〈http://www.seg.org/SEGportalWEBproject/prod/SEG-Publications/Pub-Geophysics/Documents/dissertationform.pdf〉 at http://seg.org/dissertationabstracts. Abstracts are reviewed and accepted or rejected based on their relevance to the readers of Geophysics. Abstracts must be written in English and defended in 2009 or later. Month and year defended: March, 2010 http://library.ut.ac.ir/english/ (not posted yet) Ill-posed inverse problems play an important role in geophysical studies. Basic information about geophysical models is needed to find …