Band-limited Data Research Articles

Multichannel seismic attenuation compensation and interpolation with curvelet sparse constraint of frequency-wavenumber spectrum

Abstract High-resolution exploration is hampered by seismic attenuation, caused by the viscosity and heterogeneity of underground media. Conventional single-channel attenuation compensation methods can increase the corresponding vertical resolution partially. While the lateral continuity of compensated seismic data is always ignored and the noise resistance can be improved. Therefore, multichannel attenuation compensation methods were proposed, including the algorithm with sparse curvelet coefficient constraint of the time-space (t-x) data. However, nonstationary seismic data may be affected by irregular missing traces in field cases, which severely degrades its lateral continuity and negatively affects the performance of multichannel attenuation compensation. Therefore, we concentrate on simultaneous multichannel attenuation compensation and missing trace interpolation in a unified framework. Based on the inversion framework of sparsity promotion, we propose an approach for simultaneous multichannel compensation and interpolation, utilizing sparse curvelet coefficient constraint of the recovered principal frequency-wavenumber (f-k) spectrum. The size of the principal f-k spectrum is reduced by at least half compared with that of the corresponding t-x band-limited data. It significantly reduces the computational expense of curvelet transform-based processing. Synthetic and field data experiments validate the effectiveness of the proposed method in efficiency improvement with consistent performance when compared to the multichannel method with sparse curvelet coefficient constraint of the t-x data in improving the vertical resolution and lateral continuity. Furthermore, we have discussed a potential acceleration strategy based on random sampling in the normal compensation issue.

Diving-wave time-lapse delay for CO2 thin layer detection

SUMMARY We have derived an analytical approximate expression to estimate the delay in diving seismic waves due to thin layers of CO2. The expression is valid for high frequencies and can be used to estimate the delay in diving waves at seismic frequencies for large separations between the source and receiver (offset). The approximation may be used to assess CO2 detection limits using diving waves and to support survey planning for CO2 monitoring and full-waveform inversion (FWI) cycle skipping analysis. In this study, we analyse the diving-wave response to a thin layer of CO2 for band-limited data using acoustic finite-difference modelling, and compare the results against the analytical calculations. We find that the responses are offset-dependent and related to double- and single-leg interactions between the diving waves and the CO2. To test the methods, we created a synthetic representation of the 2010 subsurface conditions for the top CO2 layer at the Sleipner storage complex in the North Sea, by combining base and monitor post-stack seismic data with field velocity trends. Using the acoustic finite-difference method, we model pre-stack data that captures the complexity of field data and demonstrate the use of the diving-wave delay for CO2 migration monitoring and CO2 thin layer detection.

A tutorial of image-domain least-squares reverse time migration through point-spread functions

Least-squares reverse time migration (LSRTM) has shown great potential to improve the amplitude fidelity and spatial resolution of a reverse time migration (RTM) image. However, the main disadvantage is that it requires significant computational resources for the iterative solution. To ameliorate this problem, the image-domain LSRTM (IDLSRTM) approach through point-spread functions (PSFs) has been proven to be a viable and alternative technique to deconvolve the standard RTM image. In this paper, we develop a tutorial for the numerical implementation of IDLSRTM through PSFs. The Hessian matrix in the standard IDLSRTM approach is estimated through the spatial interpolation of precomputed PSFs on the fly, where the PSFs are computed by one-round forward modeling and migration. However, the resulting IDLSRTM scheme is highly ill-conditioned because of the incomplete acquisition condition, irregular subsurface illumination, and band-limited data. To stabilize inversion and improve the inverted image, we suggest that the PSFs and RTM images should be deblurred by applying a deblurring filter. The deblurring PSFs can reduce the condition number of the standard Hessian matrix and make the inverse problem more well conditioned, thus improving imaging quality and accelerating convergence. Furthermore, the regularization operators should be imposed to produce a reasonable inverted image and avoid the overfitting of the image-matching term. Based on several examples with synthetic and field data, we can conclude two significant points. The first point is that the standard IDLSRTM approach through the conventional PSFs can only recover a reflectivity image similar to the deblurring RTM and nonstationary matching filter (NMF) approaches. The second point is that the deblurring IDLSRTM approach through the deblurring PSFs can retrieve the reflectivity image with higher resolution and better amplitude fidelity than the standard IDLSRTM approach and the standard RTM approach with the deblurring filter or NMF.

Data-driven suppression of short-period multiples from laterally varying thin-layered overburden structures

Marchenko multiple elimination methods remove all orders of overburden-generated internal multiples in a data-driven way. In the presence of thin beds, however, these methods have been shown to underperform. This is because the underlying inverse problem requires the information about short-period internal multiple (SPIM) imprint on the inverse transmission to be correctly constrained. This has been addressed in 1.5D media with energy conservation and minimum-phase reconstruction. Extending the applications to two dimensions and, hence, making the step toward field data were believed to be hampered by the need for a multidimensional minimum-phase reconstruction which is (1) not unique and (2) no algorithm has been found to perform this in practice on band-limited data. Here, we address both of these problems with an approach that includes solving the Marchenko equation with a trivial constraint, evaluating the energy conservation condition of its solutions to find the spatially dependent error syndrome, using the 1.5D minimum-phase reconstruction for each shot gather to find the spatially dependent constraint, and finally using that inside another run of the Marchenko equation solver to find a much-improved result. We find that the method works because in 2D media the expression of SPIMs in the inverse transmission coda is approximately 1.5D. We then investigate a class of models and synthetic data sets to verify where the 1.5D approximation starts breaking down. Our analysis indicates that this approach could perform very well in settings with moderate lateral variations, which also is where the (short-period) internal multiples are most difficult to differentiate from primary reflections.

An Optimized Heading De-Noising Method of the Polarization Compass for Improving Orientation Accuracy

In recent years, the polarization compass (PC) has been extensively applied in the field of autonomous navigation due to its strong autonomy and anti-interference characteristics. Nevertheless, the existing research shows that the PC with high extinction ratio polarization sensor still offers poor orientation accuracy since there is a great quantity of noise in the heading data output from the compass. Here, we present a novel and promising scheme to denoise heading data to improve deficiency in heading data denoising of existing research. The purpose of this study aims to solve the problem that the orientation accuracy of PC caused by noise needs to be improved. Concretely, the major points of this research include an effective separation strategy of noise and useful signal based on the optimized variational mode decomposition (OVMD) and an effective heading data denoising method based on an adaptive singular spectrum analysis (AD-SSA). OVMD based on the center frequency is employed to decompose the original noisy heading data into a finite number of band-limited data intrinsic mode functions (DIMFs) so that the best number of decomposed layers for the noising data can be obtained. Subsequently, the above decomposed DIMFs are divided into low-frequency signal components, mixed components, and high-frequency noise components exploiting the maximum similar eigenvalue of the autocorrelation matrix corresponding to each DIMF as a classification method. An AD-SSA uses the correlation coefficient to further suppress the noise in classified DIMFs. These approaches are combined effectively to improve the orientation accuracy of the PC significantly and very contribute to the navigation application of unmanned mobile platforms such as unmanned aerial vehicles (UAVs). Compared to existing prior arts, the effectiveness of the proposed heading denoising scheme for the PC is demonstrated through the results obtained from the heading data sets on two representative sunny and dust days.

Graph learning from band-limited data by graph Fourier transform analysis

A graph provides an effective means to represent the statistical dependence or similarity among signals observed at different vertices. A critical challenge is to excavate graphs underlying observed signals, because of non-convex problem structure and associated high computational requirements. This paper presents a new graph learning technique that is able to efficiently infer the graph structure underlying observed graph signals. The key idea is that we reveal the intrinsic relation between the frequency-domain representation of general band-limited graph signals, and the graph Fourier transform (GFT) basis. Accordingly, we derive a new closed-form analytic expression for the GFT basis, which depends deterministically on the observed signals (as opposed to being solved numerically and approximately in the literature). Given the GFT basis, the estimation of the graph Laplacian, more explicitly, its eigenvalues, is convex and efficiently solved using the alternating direction method of multipliers (ADMM). Simulations based on synthetic data and experiments based on public weather and brain signal datasets show that the new technique outperforms the state of the art in accuracy and efficiency.

Harmonic Extrapolation of Seismic Reflectivity Spectrum for Resolution Enhancement: An Insight from Inas Field, Offshore Malay Basin

For decades, how to improve the resolution limit of seismic data has been a concern for seismologists. For this reason, several geoscientists have proposed various methods of improving the bandwidth of the data. Discovering an easy and scientifically reliable means of improving seismic data resolution would undoubtedly help geophysicists interpret more complex details of the subsurface geology. In this study, we transformed the bandlimited time-domain seismic data to the frequency domain using the Fourier analysis method, and a basis pursuit atomic algorithm was applied to decompose the real and imaginary parts of the spectrum into summations of cosines and sines. The resultant reflectivity spectrum (in the frequency domain) was deconvolved by a pre-estimated wavelet spectrum to obtain the true earth’s reflectivity data spectrum and was subsequently extrapolated to beyond the original band limit. The result shows an extended bandwidth from 68 Hz to 161 Hz and 80 Hz to 170 Hz for both synthetic trace model and the main seismic data, respectively. Consequently, this improved vertical resolution of sub-seismic geologic features, such as crevasse splay, levee, barrier bar complex, lagoon inlet channels, alluvial fans, and fluvial channels, and shows subtle facies variations in Inas field.

Geostatistical Inversion of Spectrally Broadened Seismic Data for Re-Evaluation of Oil Reservoir Continuity in Inas Field, Offshore Malay Basin

The application of geostatistics in seismic inversion techniques has been proven somewhat reliable in the delineation of reservoir properties and has recently attracted the attention of many geoscientists. However, there are cases where its prediction returned negative results after drilling. In this research, we re-evaluated a reservoir in Inas Field, whose geostatistical inversion result wrongly predicted sand continuity, resulting in the spudding of a dry hole. When a geostatistical seismic inversion is successfully applied, it provides an increase in seismic resolution and aids the prediction of sand continuity. Although this method relies more on the statistical data from a well because of the limitation of the seismic data in resolving thin geologic features, the spatial variation of reservoir parameters still depends on seismic data, which often have poor resolution quality. Therefore, to investigate the impact of bandlimited data on the geostatistical inversion, we harmonically extended the seismic bandwidth by applying a sparse-layer spectral inversion algorithm to the data. This algorithm increased the seismic data bandwidth from 80 Hz to 180 Hz, and its tuning thickness reduced from 32 m to 10 m at the reservoir interval. The resultant broadband (180 Hz), as well as the original seismic (narrowband of 80 Hz) data, were both used as input to build two separate geostatistical prediction models, respectively. Twenty (20) realizations of these models were generated, ranked into P10, P50, and P90, and the best case was selected for interpretation. These realizations were used to characterize the reservoir lithofacies distribution. When compared, the result of the broadband inversion, facies and sand distribution model showed that the reservoir facies changed towards the location of the dry well. The broadband geostatistical inversion efficiently improved the reservoir characterization process by not only producing an accurate estimation of the lateral extent of the reservoir heterogeneities but also generating outcomes that help us understand why other geostatistical inversion analyses of the target reservoir were misleading. Contrary to the popular assumption, it was discovered that the tuning effects of bandlimited data could affect the result of a geostatistical inversion and result in wrong facies predictions.

Deep Learning for Low-Frequency Extrapolation of Multicomponent Data in Elastic FWI

Full waveform inversion (FWI) strongly depends on an accurate starting model to succeed. This is particularly true in the elastic regime: The cycle-skipping phenomenon is more severe in elastic FWI compared to acoustic FWI, due to the short S-wave wavelength. In this paper, we extend our work on extrapolated FWI (EFWI) by proposing to synthesize the low frequencies of multi-component elastic seismic records, and use those "artificial" low frequencies to seed the frequency sweep of elastic FWI. Our solution involves deep learning: we separately train the same convolutional neural network (CNN) on two training datasets, one with vertical components and one with horizontal components of particle velocities, to extrapolate the low frequencies of elastic data. The architecture of this CNN is designed with a large receptive field, by either large convolutional kernels or dilated convolution. Numerical examples on the Marmousi2 model show that the 2-4Hz low frequency data extrapolated from band-limited data above 4Hz provide good starting models for elastic FWI of P-wave and S-wave velocities. Additionally, we study the generalization ability of the proposed neural network over different physical models. For elastic test data, collecting the training dataset by elastic simulation shows better extrapolation accuracy than acoustic simulation, i.e., a smaller generalization gap.

Self-Supervised Deep Learning for Nonlinear Seismic Full Waveform Inversion

Seismic full waveform inversion (FWI) is able to build high-resolution velocity model based on the full information carried by seismic wave. However, FWI requires an accurate enough initial model to ensure convergence. In this paper, we propose a new nonlinear FWI method to mitigate the initial model dependence problem. Specifically, we firstly propose a nonlinear operator within the hybrid model- and data-driven framework based on the frequency controllable envelope operator (FCEO) and a deep learning architecture U-Net. FCEO is used to obtain the envelope of a band-limited data and U-Net realizes the mapping from this envelope to that corresponding to a lower frequency band. The U-Net is trained in a self-supervised manner that avoids the reliance on labeled data and benefits the generalization ability. Based on the nonlinear operator, a nonlinear FWI method is proposed by defining a new misfit function. In addition, the calculation of gradient is derived using the adjoint-state method. Using numerical examples, we investigate the performance of the proposed nonlinear operator and the new nonlinear FWI method. The results clearly demonstrate that the proposed nonlinear operator is effective in obtaining low-frequency envelope data, and the new nonlinear FWI method has advantages over common method in mitigating cycle-skipping and building an initial model for conventional FWI.

Time-Domain Multidimensional Deconvolution: A Physically Reliable and Stable Preconditioned Implementation

Multidimensional deconvolution constitutes an essential operation in a variety of geophysical scenarios at different scales ranging from reservoir to crustal, as it appears in applications such as surface multiple elimination, target-oriented redatuming, and interferometric body-wave retrieval just to name a few. Depending on the use case, active, microseismic, or teleseismic signals are used to reconstruct the broadband response that would have been recorded between two observation points as if one were a virtual source. Reconstructing such a response relies on the the solution of an ill-conditioned linear inverse problem sensitive to noise and artifacts due to incomplete acquisition, limited sources, and band-limited data. Typically, this inversion is performed in the Fourier domain where the inverse problem is solved per frequency via direct or iterative solvers. While this inversion is in theory meant to remove spurious events from cross-correlation gathers and to correct amplitudes, difficulties arise in the estimation of optimal regularization parameters, which are worsened by the fact they must be estimated at each frequency independently. Here we show the benefits of formulating the problem in the time domain and introduce a number of physical constraints that naturally drive the inversion towards a reduced set of stable, meaningful solutions. By exploiting reciprocity, time causality, and frequency-wavenumber locality a set of preconditioners are included at minimal additional cost as a way to alleviate the dependency on an optimal damping parameter to stabilize the inversion. With an interferometric redatuming example, we demonstrate how our time domain implementation successfully reconstructs the overburden-free reflection response beneath a complex salt body from noise-contaminated up- and down-going transmission responses at the target level.

New algorithms for spherical harmonic analysis of area mean values over blocks delineated by equiangular and Gaussian grids

Spherical harmonic analysis is widely used in all aspects of geoscience. Exact quadrature methods are available for the spherical harmonic analysis of band-limited point values at the grid points of equiangular and Gaussian grids. However, no similarly exact quadrature methods are available for the spherical harmonic analysis of area mean values over the blocks delineated by these grids. In this study, new algorithms appropriate for the exact spherical harmonic analysis of the band-limited area mean values over the blocks delineated by equiangular and Gaussian grids are proposed. For band-limited data, precision that is between that of the least-squares estimation method and of the approximate quadrature methods can be achieved by using the new algorithms. Regarding the computational complexity, fewer operations are needed by the new methods as compared to those needed by the least-squares estimation method and the approximate quadrature methods in the preparation stage when the maximum degree of the spherical harmonic analysis is very large. Simulation experiments are performed to compare the ability to recover the spherical harmonic coefficients by using the least-squares estimation method, the approximate quadrature methods and these new algorithms from aliased data with aliasing components of realistic magnitudes. The results suggest that these new algorithms, with time complexity one order less than that of the least-squares estimation method in the solving stage, perform roughly the same as the least-squares estimation method in recovering spherical harmonic coefficients from the aliased data.

Intrinsic non-uniqueness of the acoustic full waveform inverse problem

SUMMARYIn the context of seismic imaging, full waveform inversion (FWI) is increasingly popular. Because of its lower numerical cost, the acoustic approximation is often used, especially at the exploration geophysics scale, both for tests and for real data. Moreover, some research domains such as helioseismology face true acoustic media for which FWI can be useful. In this work, an argument that combines particle relabelling and homogenization is used to show that the general acoustic inverse problem based on band-limited data is intrinsically non-unique. It follows that the results of such inversions should be interpreted with caution. To illustrate these ideas, we consider 2-D numerical FWI examples based on a Gauss–Newton iterative inversion scheme and demonstrate effects of this non-uniqueness in the local optimization context.

Extrapolated full-waveform inversion with deep learning

The lack of low-frequency information and a good initial model can seriously affect the success of full-waveform inversion (FWI), due to the inherent cycle skipping problem. Computational low-frequency extrapolation is in principle the most direct way to address this issue. By considering bandwidth extension as a regression problem in machine learning, we have adopted an architecture of convolutional neural network (CNN) to automatically extrapolate the missing low frequencies. The band-limited recordings are the inputs of the CNN, and, in our numerical experiments, a neural network trained from enough samples can predict a reasonable approximation to the seismograms in the unobserved low-frequency band, in phase and in amplitude. The numerical experiments considered are set up on simulated P-wave data. In extrapolated FWI (EFWI), the low-wavenumber components of the model are determined from the extrapolated low frequencies, before proceeding with a frequency sweep of the band-limited data. The introduced deep-learning method of low-frequency extrapolation shows adequate generalizability for the initialization step of EFWI. Numerical examples show that the neural network trained on several submodels of the Marmousi model is able to predict the low frequencies for the BP 2004 benchmark model. Additionally, the neural network can robustly process seismic data with uncertainties due to the existence of random noise, a poorly known source wavelet, and a different finite-difference scheme in the forward modeling operator. Finally, this approach is not subject to strong assumptions on signals or velocity models of other methods for bandwidth extension and seems to offer a tantalizing solution to the problem of properly initializing FWI.

Goal-Oriented Inversion-Based NMO Correction Using a Convex $l_{2,1}$ -Norm

Normal moveout (NMO) correction is a routine step in seismic data processing, which has an important impact on other seismic processing procedures, seismic inversion, and interpretation. We propose a goal-oriented inversion-based NMO correction method using a convex $l_{2,1}$ -norm. The proposed method corrects the superresolution multichannel offset-dependent reflectivity rather than the bandlimited data itself sample by sample, block by block, or wavelet by wavelet. Therefore, the proposed method can essentially reduce the amplitude and even phase distortion introduced by data-based NMO correction methods in the presence of strong wavelet interference. We impose two goal-oriented constraints including both the temporal sparsity and the horizontal continuity of reflectivity, which are approximately represented by a convex $l_{2,1}$ -norm, on the geometric moveout relationship from offset-dependent trajectories to zero offsets to build a new objective function for NMO correction. The goal-desired temporal sparsity of reflectivity can induce the superresolution solution; meanwhile, the goal-desired horizontal continuity introduces a reasonable intrinsic structure to further limit the solution space and is particularly suitable to processing interfering reflections. Attributing to these two additional constraints, the new NMO correction method can flatten the interfering events and the intersecting events with favorable offset-dependent amplitude and phase variations even in the presence of noise. Synthetic and real data examples are adopted to verify the performance of our method. The results show that goal-oriented inversion-based NMO correction using the $l_{2,1}$ -norm is a potentially effective, stable, and high-quality NMO correction tool, especially for strong wavelet interference and at far offsets.

Performance of Vector Fitting Algorithm Applied to Bandpass and Baseband Systems

This article presents the performance evaluation of Vector Fitting Algorithm (VFA) from a system identification perspective. In this paper, VFA has been first applied to known baseband and bandpass systems such as Butterworth lowpass and bandpass filters to analyze the algorithm’s pole-residue extraction ability for band-limited noisy data. The poles identified by the algorithm for different bandwidths and noise powers are compared with the actual system poles of the baseband and bandpass systems. It is concluded that the algorithm is capable of identifying the actual system poles even if the capture bandwidth is less than the 3 dB bandwidth, which is a significant observation of this paper. It is also seen that the system identification performance with noisy data is better for baseband systems when compared to bandpass systems. Further, a practical investigation has been done to evaluate VFA performance for modeling a microstrip coupled line filter in the presence of noise.

Squeezing Kappa (κ) Out of the Transportable Array: A Strategy for Using Bandlimited Data in Regions of Sparse Seismicity

The κ parameter (Anderson and Hough, 1984), and namely its site‐specific component ( κ ), is important for predicting and simulating high‐frequency ground motion. We develop a framework for estimating κ and addressing uncertainties under the challenging conditions often imposed in practice: (1) low seismicity (limited, poor‐quality, distant records); (2) limited‐bandwidth data from the Transportable Array (TA; maximum usable frequency 16 Hz); and (3) low magnitudes ( M L 1.2–3.4) and large uncertainty in stress drop (corner frequency). We cannot resolve stress drop within the bandwidth, so we propose an approach that only requires upper and lower bounds on its regional values to estimate κ . To address uncertainties, we combine three measurement approaches (acceleration spectrum slope [AS]; displacement spectrum slope [DS]; and broadband spectral fit). We also examine the effect of crustal amplification, and find that neglecting it can affect κ by up to 35%. DS estimates greatly exceed AS estimates. We propose a reason behind this bias, related to the residual effect of the corner frequency on κ AS and κ DS. For our region, we estimate a frequency‐independent mean S ‐wave Q of 900±300 at 9–16 Hz, and an ensemble mean κ over all sites of 0.033±0.014 s. This value is similar to the native κ of the Next Generation Attenuation‐West 2 ground‐motion prediction equations, indicating that these do not need to be adjusted for κ for use in southern Arizona. We find that stress‐drop values in this region may be higher compared to estimates of previous studies, possibly due to trade‐offs between stress drop and κ . For this dataset, the within‐approach uncertainty is much larger than the between‐approach uncertainty, and it cannot be reduced if the data quality is not improved. The challenges discussed here will be relevant in studies of κ for other regions with bandlimited data; for example, any region where data come primarily from the TA.

Using instantaneous frequency and colored inversion attributes to distinguish and determine the sandstones facies of the Late Ordovician Mamuniyat reservoir, R-field in Murzuq Basin, Libya

The Mamuniyat petroleum reservoir in southwestern Libya is comprised of clean sandstones and intercalated shale and sand facies that are characterized by spatial porosity variations. Seismic reflection data from the field exhibit relatively low vertical seismic resolution, side lobes of reflection wavelets, reflection interference, and low acoustic impedance contrast between the reservoir and the units underneath the reservoir, which make mapping those facies a difficult task. In the absence of broadband seismic data, optimizing frequency bands of bandlimited data can be used to suppress pseudoreflectors resulting from side-lobe effects and help to separate the clean sandstone facies of the reservoir. We have optimized the data based on our investigation of seismic frequency bands and used instantaneous frequency analysis to reveal the reflection discontinuity that is mainly associated with the reservoir boundary of the sandstone facies of the clean Mamuniyat reservoir. We also preformed rock-physics diagnostic modeling and inverted the seismic data using spectral-based colored inversion into relative acoustic impedance. The inverted impedance matches the up-scaled impedance from the well data and the inversion of relative acoustic impedance confirms the conclusion that was drawn from the instantaneous frequency results. The interpretation of facies distributions based on the instantaneous frequency was supported by the inversion results and the rock-physics models.

Feasibility of approximate broadband estimation of acoustic impedance profile from first principles and band-limited reflection data

Because seismic reflection data are band limited, acoustic impedance profiles derived from them are nonunique. The conventional inversion methods counter the nonuniqueness either by stabilizing the answer with respect to an initial model or by imposing mathematical constraints such as sparsity of the reflection coefficients. By making a nominal assumption of an earth model locally consisting of a stack of homogeneous and horizontal layers, we have formulated a set of linear equations in which the reflection coefficients are the unknowns and the recursively integrated seismic trace constitute the data. Drawing only on first principles, the Zoeppritz equation in this case, the approach makes a frontal assault on the problem of reconstructing reflection coefficients from band-limited data. The local layer-cake assumption and the strategy of seeking a singular value decomposition solution of the linear equations counter the nonuniqueness, provided that the objective is to reconstruct a smooth version of the impedance profile that includes only its crude structures. Tests on synthetic data generated from elementary models and from measured logs of acoustic impedance demonstrated the efficacy of the method, even when a significant amount of noise was added to the data. The emergence of consistent estimates of impedance, approximating the original impedance, from synthetic data generated for several frequency bands has inspired our confidence in the method. The other attractive outputs of the method are as follows: (1) an accurate estimate of the impedance mean, (2) an accurate reconstruction of the direct-current (DC) frequency of the reflectivity, and (3) an acceptable reconstruction of the broad outline of the original impedance profile. These outputs can serve as constraints for either more refined inversions or geologic interpretations. Beginning from the restriction of band-limited data, we have devised a method that neither requires a starting input model nor imposes mathematical constraints on the earth reflectivity and still yielded significant and relevant geologic information.

Seismic sparse-spike deconvolution via Toeplitz-sparse matrix factorization

We have developed a new sparse-spike deconvolution (SSD) method based on Toeplitz-sparse matrix factorization (TSMF), a bilinear decomposition of a matrix into the product of a Toeplitz matrix and a sparse matrix, to address the problems of lateral continuity, effects of noise, and wavelet estimation error in SSD. Assuming the convolution model, a constant source wavelet, and the sparse reflectivity, a seismic profile can be considered as a matrix that is the product of a Toeplitz wavelet matrix and a sparse reflectivity matrix. Thus, we have developed an algorithm of TSMF to simultaneously deconvolve the seismic matrix into a wavelet matrix and a reflectivity matrix by alternatively solving two inversion subproblems related to the Toeplitz wavelet matrix and sparse reflectivity matrix, respectively. Because the seismic wavelet is usually compact and smooth, the fused Lasso was used to constrain the elements in the Toeplitz wavelet matrix. Moreover, due to the limitations of computer memory, large seismic data sets were divided into blocks, and the average of the source wavelets deconvolved from these blocks via TSMF-based SSD was used as the final estimation of the source wavelet for all blocks to deconvolve the reflectivity; thus, the lateral continuity of the seismic data can be maintained. The advantages of the proposed deconvolution method include using multiple traces to reduce the effect of random noise, tolerance to errors in the initial wavelet estimation, and the ability to preserve the complex structure of the seismic data without using any lateral constraints. Our tests on the synthetic seismic data from the Marmousi2 model and a section of field seismic data demonstrate that the proposed method can effectively derive the wavelet and reflectivity simultaneously from band-limited data with appropriate lateral coherence, even when the seismic data are contaminated by noise and the initial wavelet estimation is inaccurate.