Earth Impulse Response Research Articles

Two unequal shots in the same hole: statistical versus deterministic deconvolution

Determination of the source time function and earth impulse response for explosive source seismic reflection data is traditionally performed using predictive deconvolution, which assumes the reflection seismogram is the convolution of a minimum-phase wavelet with a white, random, stationary sequence of impulses. We compare this statistical method with a new deterministic method that uses two unequal shots close together in the same hole and assumes a seismogram is the convolution of a source time function with a Green's function, plus noise. The seismograms from the two shots share essentially the same Green's function to any receiver, because the shots are close together. The Green's function is estimated by deconvolving the two seismograms at a given receiver for the corresponding source time functions, modelled independently; source parameters are obtained by trial-and-error by minimising the difference between the two estimates of the Green's function. We test these competing theories using seismic reflection data acquired with two charges fired separately in the same hole into a receiver spread of 1178 bunched geophone groups at 10 m intervals: a 1 kg charge at 42 m depth and 2 kg at 40 m depth. The criterion for comparison of the theories is the similarity of the recovered impulse responses at any receiver. Each method recovers a wavelet, or source time function, for each shot and two estimates of the impulse response, or Green's function, at each receiver. The Green's functions recovered from the two shots using the deterministic method are almost identical and are measurably more similar than the impulse responses recovered using predictive deconvolution. That is, the deterministic method is better. The resolution of explosive source seismic reflection data can now be significantly improved with the new method and minimal additional field effort of loading one extra charge per shot hole.

Extracting Earth Impulse Response From Electromagnetic Data With Pseudorandom Current Using an Improved Correlation Identification

Extracting Earth Impulse Response From Electromagnetic Data With Pseudorandom Current Using an Improved Correlation Identification

Extended spectrum modulation method for MTEM harmonic noise filtering

Harmonic noise represented by power-line noise is the main interference of Multi-transient electromagnetic (MTEM), and traditional filtering methods inevitably suffer from noise residue and waveform distortion. The key step in MTEM is the extraction of the earth impulse response, which is equivalent to filtering the voltage response signal using an identification system. This paper modulates the spectrum of the identification system by extending the source signal to increase the attenuation at these target frequencies in order to filter out harmonic noise when extracting the earth impulse response, which we named as Extended Spectrum Modulation Method (ESMM). This paper first deduces the mathematical relationship between data length and target frequency. The filtering of power-line harmonics by ESMM is then demonstrated and compared with a conventional notch filter. This paper also proposes a rapid deployment method for ESMM to improve data utilization. Since ESMM filters out noise while extracting the earth impulse response, there is no additional hardware cost or time overhead. Finally, a field experiment demonstrates the superior performance and ultimate convenience of ESMM.

IMS observations of infrasound and acoustic-gravity waves produced by the January 2022 volcanic eruption of Hunga, Tonga: A global analysis

The 15 January 2022 Hunga, Tonga, volcano's explosive eruption produced the most powerful blast recorded in the last century, with an estimated equivalent TNT yield of 100–200 megatons. The blast energy was propagated through the atmosphere as various wave types. The most prominent wave was a long-period (>2000 s) surface-guided Lamb wave with energy comparable to that of the 1883 Krakatoa Lamb wave; both were clearly observed by pressure sensors (barometers) worldwide. Internal gravity, acoustic-gravity, and infrasound waves were captured in great detail by the entire infrasound component of the International Monitoring System (IMS). For instance, infrasound waves (<300 s period) were seen to circumnavigate Earth up to eight times. Atmospheric waves captured by the IMS infrasound network and selected barometers near the source provide insight on Earth's impulse response at planetary scales.

Electric field calculations for real-time space weather alerting systems

SUMMARY Space weather alerting systems for power systems require real-time calculations of the electric fields that drive geomagnetically induced currents. In this paper, we present a new method for calculating the Earth impulse response that can be convolved with real-time magnetic field data to give the required electric fields. We start with the Earth transfer function which can be expressed in two ways: as a relationship between the electric field, E, and the geomagnetic field, B, which has characteristics equivalent to that of a high-pass filter, and between E and the time derivative of the magnetic field, dB/dt, which has characteristics of a low-pass filter. An inverse Fourier transform of these transfer functions should then give the corresponding Earth impulse responses in the time domain. This works well for a uniform conductivity model for which the inverse Fourier transform of the transfer function has an analytic solution. However, the inverse Fourier transform of the transfer function for a non-uniform conductivity model requires numerical calculations and produces an acausal impulse response with oscillations because of the Gibbs phenomenon. To investigate the origin of the Gibbs oscillations, the real and imaginary parts of the transfer functions are transformed separately. This shows that the Gibbs oscillations arise from the inverse transform of the real and imaginary parts of the high-pass transfer function and the imaginary part of the low-pass transfer function. A new method is introduced that just transforms the real part of the low-pass transfer function and uses the requirement of causality to construct the full low-pass impulse response. From this, the derivative theorem of convolution is used to obtain the high-pass impulse response. Electric fields can then be calculated by convolution of the low-pass impulse response with the rate of change of the magnetic field, dB/dt, or by convolution of the high-pass impulse response with the magnetic field, B. Tests of the new method by comparison with analytic solutions for specified earth models and synthetic magnetic field data gave very high correlation coefficients, slopes near 1.0 and intercepts near zero, showing the accuracy of the new method. The method can be used with any Earth transfer function whether obtained from magnetotelluric measurements or from 1-D, 2-D or 3-D conductivity models. Thus, it provides a versatile technique that avoids the Gibbs phenomenon and produces a causal impulse response suitable for time domain calculations of electric fields.

Direct synthesis of time domain pseudo-random 3D electromagnetic response with a band-limited source

Methods based on pseudo-random coded emission technology can effectively suppress the noise in electromagnetic exploration and improve the depth and resolution of the detection. Additionally, the synthesis of electromagnetic response signals excited by pseudo-random code sources can provide not only theoretical data for data processing but also intermediate results for evaluations of field data. Currently, pseudo-random signal simulations ignore the presence of source effects generated by band-limited sources in the Earth's impulse response. Furthermore, the simulation methods are mainly completed based on one-dimensional and two-dimensional formats and thus cannot accurately reflect the pseudo-random response of three-dimensional geoelectric structures. In this paper, 3D pseudo-random electromagnetic signals with a band-limited source were synthesized according to the propagating routine of the electromagnetic field. First, Maxwell's Equations were solved using a three-dimensional time-domain finite volume method. To avoid the singularity of the source and incorporate source effect response, the Gaussian function was introduced to simulate the band-limited impulse sources. The Earth's impulse response was obtained, which contained source effects under three-dimensional geoelectric structures. Then, ideal pseudo-random electromagnetic signals could be composed by convoluting the Earth's impulse response with a pseudo-random source waveform. Since the transmitter and receiver instruments were band-limited, a low-pass filter with a Blackman window was introduced to simulate the instrument's system response. The simulated pseudo-random signals were obtained by adding noise and the ideal pseudo-random electromagnetic signals were filtered using a low-pass filter. By comparison with signals acquired in the field, the correctness of the pseudo-random signal synthesis with source effects was successfully verified.

Application of powerline noise cancellation method in correlation identification of controlled source electromagnetic method

Abstract Powerline interference in the controlled source electromagnetic method has traditionally been one of the biggest conundrums plaguing geophysicists, and its conventional denoising methods primarily include filtering and noise estimation. The filter method leaches noise at specific frequency points, which might also filter useful signals; the noise estimation method significantly eliminates interference, whereas the premise is that the noise is stable after a short time and a recorder is necessary in the field. In the present study, using the periodicity and symmetry of powerline noise, we propose a subtraction and an addition method for cancellation of the powerline noise. First, the transmitted signal is optimized so that the equivalent transmitted signal is an m sequence; then the response signal is processed by using the cancellation method; subsequently, the correlation identification is applied and finally, we solve the earth impulse response by means of the Wiener filter deconvolution method. Simulation experiments and field data tests demonstrate that the powerline noise can be well suppressed by the cancellation method proposed in the present study, so that the system identification accuracy is greatly improved. The method is simple in principle and effective in removing powerline noise, which presents a novel perspective on noise elimination for system identification.

Inversion of explosive source land seismic data to determine source signature parameters

ABSTRACTIn seismic land exploration, where terrain conditions limit the use of vibrator sources, explosive sources are the preferred choice, and in situations where low‐frequency energy is essential, explosive is the preferred source, because it contains energy down to DC (0 Hz). I present a new method to estimate the source time functions of explosion source seismic data, using a modified acquisition method, data processing, modelling and inversion. Two shots of different sizes are fired into the same geophone spread with source–receiver geometry arranged such that the Green's functions are essentially the same. The spectral ratio of corresponding seismic traces is then the spectral ratio of the source time functions. A corresponding synthetic ratio filter is calculated from Blake's explosion source model for the sources within each pair. Each modelled source is defined by five parameters: the minimum radius a of the elastic zone, the internal pressure p0 at a, the density of the rock ρ and two elastic constants for an isotropic rock, the P‐wave velocity c and Poisson's ratio σ. A grid search finds the source parameters that minimise the difference between the measured and synthetic filters for each source pair, thus yielding the two source time functions, but not their absolute amplitudes. Deconvolution of the shot records within each pair for the minimum‐phase source time functions recovers the estimated earth impulse response gathers, convolved only with the response of the recording system. The improved knowledge of the source signature is then helpful in subsequent processing, especially for reverse time migration, full waveform inversion and amplitude‐versus‐offset analysis. I consider shots at different depths in the same shot hole and shots at the same depth in different shot holes. The former configuration minimises differences in geophone coupling, and, in operations, quality control procedures must be met. There is considerable variation in recovered source parameters from place to place, caused by variations in charge size, soil or rock properties and source coupling.

Analysis of Pseudo-Random Sequence Correlation Identification Parameters and Anti-Noise Performance

Using a pseudo-random sequence to encode the transmitted waveform can significantly improve the working efficiency and depth of detection of electromagnetic exploration. The selection of parameters of pseudo-random sequence plays an important role in correlation identification and noise suppression. A discrete cycle correlation identification method for extracting the earth impulse response is proposed. It can suppress the distortion in the early stage of the excitation field and the glitches of the cross correlation function by traditional method. This effectively improves the accuracy of correlation identification. The influence of the order and the cycles of m-series pseudo-random coding on its autocorrelation properties is studied. The numerical results show that, with the increase of the order of m-sequence, the maximum out-of-phase periodic autocorrelation function decreases rapidly. Therefore, it is very beneficial to achieve synchronization. The limited-cycle m-sequences have good autocorrelation properties. As the period of the m-sequence increases and the width of the symbol decreases, the overall autocorrelation becomes closer to the impact function. The discussion of the influence of symbol width and period of m-sequence on its frequency bandwidth and power spectral density shows that the narrower the symbol width, the wider its occupied band. The longer the period, the smaller the power spectral line spacing. The abilities of m-sequence to suppress DC (Direct-current) interference, Schumann frequency noise, and sine-wave noise are analyzed. Numerical results show that the m-sequence has excellent ability to suppress DC interference and Schumann frequency noise. However, for high-order harmonic noise, the correlation identification error appears severe oscillation in the middle and late stages of the impulse response. It indicates that the ability of m-sequence to suppress high-frequency sinusoidal noise is deteriorated. In practical applications, the parameters of the transmitted waveform should be reasonably selected in combination with factors including transmitter performance, hardware noise, and ambient noise level to achieve the best identification effect.

Development of a full-waveform voltage and current recording device for multichannel transient electromagnetic transmitters

Abstract. Due to the pressing demand for metallic ore exploration technology in China, several new technologies are being employed in the relevant exploration instruments. In addition to possessing the high resolution of the traditional transient electromagnetic method, high-efficiency measurements, and a short measurement time, the multichannel transient electromagnetic method (MTEM) technology can also sensitively determine the characteristics of a low-resistivity geologic body, without being affected by the terrain. Besides, the MTEM technology also solves the critical, existing interference problem in electrical exploration technology. This study develops a full-waveform voltage and current recording device for MTEM transmitters. After continuous acquisition and storage of the large, pseudo-random current signals emitted by the MTEM transmitter, these signals are then convoluted with the signals collected by the receiver to obtain the earth's impulse response. In this paper, the overall design of the full-waveform recording apparatus, including the hardware and upper-computer software designs, the software interface display, and the results of field test, is discussed in detail.

Qualitative analysis of MTEM response using instantaneous attributes

This paper introduces new technique for qualitative analysis of multi-transient electromagnetic (MTEM) earth impulse response over complex geological structures. Instantaneous phase and frequency attributes were used in place of the conventional common offset section for improved qualitative interpretation of MTEM data by obtaining more detailed information from the earth impulse response. The instantaneous attributes were used to describe the lateral variation in subsurface resistivity and the visible geological structure with respect to given offsets. Instantaneous phase attribute was obtained by converting the impulse response into a complex form using the Hilbert transform. Conversely, the polynomial phase difference (PPD) estimator was favored over the center finite difference (CFD) approximation method in calculating the instantaneous frequency attribute because it is computationally efficient and has the ability to give a smooth variation of the instantaneous frequency over a common offset section. The observed results from the instantaneous attributes were in good agreement with both the subsurface model used and the apparent resistivity section obtained from the MTEM earth impulse response. Hence, this study confirms the capability of both instantaneous phase and frequency attributes as highly effective tools for MTEM qualitative analysis.

Improved data segmentation method for EM excited by m-sequence: A new approach in powerline noise reduction

Electromagnetic (EM) data obtained by using an m-sequence source is typically segmented with an m-sequence cycle, and the random noise is suppressed by stacking the earth impulse responses recovered from each data segment. However, the stack cannot suppress the powerline noise, which is normally several orders of magnitude greater than the random noise. In this paper, the powerline noise reduction with different data lengths were compared by numerical simulation; the results show that there is always an optimal data length, not a single cycle but several cycles that can achieve maximum powerline noise suppression. The concept of an identification system was proposed and the identified earth impulse response is obtained by measuring the system output with the input of the voltage response. The optimal data length was summarized in two cases by analyzing the identification system with different data lengths. Coding frequency optimization was proposed to shorten the optimal data length to increase the stacking times. Finally, the validity of the theory described in this paper was verified by field data processing. Different from the traditional subtraction method and notch filter, the powerline suppression method presented here does not need a noise record or filter design; the powerline noise is suppressed accompanied by the identification processing without any additional workload; the simple and practical method provides a new insight into the powerline noise reduction. With this method, the utilization of field data is improved and the anti-noise performance of the m-sequence is also fully exploited.

2D Multitransient Electromagnetic Response Modeling of South China Shale Gas Earth Model Using an Approximation of Finite Difference Time Domain with Uniaxial Perfectly Matched Layer

In this study, we introduced multitransient electromagnetic (MTEM) method as an effective tool for shale gas exploration. We combined the uniaxial perfectly matched layer (UPML) equation with the first derivative diffusion equation to solve for a finite difference time domain (FDTD) UPML equation, which was discretized to form an algorithm for 3D modeling of earth impulse response and used in modeling MTEM response over 2D South China shale gas model. We started with stepwise demonstration of the UPML and the FDTD algorithm as an effective tool. Subsequently, quantitative study on the convergence of MTEM earth impulse response was performed using different grid setup over a uniform earth material. This illustrates that accurate results can be obtained for specified range of offset. Furthermore, synthetic responses were generated for a set of geological scenarios. Lastly, the FDTD algorithm was used to model the MTEM response over a 2D shale gas earth model from South China using a PRBS source. The obtained apparent resistivity section from the MTEM response showed a similar geological setup with the modeled 2D South China shale gas section. This study confirmed the competence of MTEM method as an effective tool for unconventional shale gas prospecting and exploitation.

Wiener estimation of the Green’s function

I consider the problem of finding the impulse response, or Green’s function, from a measured response including noise, given an estimate of the source time function. This process is usually known as signature deconvolution. Classical signature deconvolution provides no measure of the quality of the result and does not separate signal from noise. Recovery of the earth impulse response is here formulated as the calculation of a Wiener filter in which the estimated source signature is the input and the measured response is the desired output. Convolution of this filter with the estimated source signature is the part of the measured response that is correlated with the estimated signature. Subtraction of the correlated part from the measured response yields the estimated noise, or the uncorrelated part. The fraction of energy not contained in this uncorrelated component is defined as the quality of the filter. If the estimated source signature contains errors, the estimated earth impulse response is incomplete, and the estimated noise contains signal, recognizable as trace-to-trace correlation. The method can be applied to many types of geophysical data, including earthquake seismic data, exploration seismic data, and controlled source electromagnetic data; it is illustrated here with examples of marine seismic and marine transient electromagnetic data.

Comparison of pseudo‐random binary sequence and square‐wave transient controlled‐source electromagnetic data over the Peon gas discovery, Norway

ABSTRACTWe discuss the problem of source control in controlled‐source electromagnetic (CSEM) surveying and compare and contrast equal energy transient square‐wave and transient pseudo‐random binary sequence source signatures for the same towed‐streamer electromagnetic survey line over the Peon gasfield in the Norwegian sector of the North Sea. The received response of the transient square‐wave data was 11 dB greater than that of the pseudo‐random binary sequence data, due to diffusive attenuation of higher frequencies present in the more broadband pseudo‐random binary sequence signature. Deconvolution of the pseudo‐random binary sequence data recovers the total impulse response function, increases the signal‐to‐noise ratio by 32.6 dB and separates most of the air wave from the earth impulse response by the causality principle. The recovered impulse responses have more detailed information in the frequency domain than the transient square‐wave data. The pseudo‐random binary sequence data were acquired with a 10 Hz source bit rate but contain no information about the Peon gasfield at frequencies above 2 Hz. The bit rate could have been reduced to 4 Hz, increasing the signal energy below 2 Hz by 150% and thus, potentially, increasing the signal‐to‐noise ratio by a further 4 dB.Because the total earth impulse response can be recovered from the broad‐bandwidth pseudo‐random binary sequence data, further time‐domain processing may be applied, including correlated noise removal, which can increase the signal‐to‐noise ratio by as much as 20 dB, and air wave removal using the causality principle. The information in the arrival time of the peak of the earth response provides the potential for traveltime to resistivity mapping to provide a starting model for inversion.

Lagrangian and energy forms for retrieving the impulse response of the Earth due to random electromagnetic forcing

We distinguish between trivial and nontrivial differences in retrieving the real or imaginary parts of the Green's function. Trivial differences come from different Green's function definitions. The energy and lagrangian forms constitute nontrivial differences. Magnetic noise sources suffice to extract the quasistatic electromagnetic-field Earth impulse response in the lagrangian form. This is of interest for Earth subsurface imaging. A numerical example demonstrates that all source vector components are necessary to extract a single-field vector component.

Analytic minimum information deconvolution and its application to ocean bottom seismometer data

Seismic refraction/wide angle reflection data for the study of the earth's crust generally lack resolution due to the low frequency wavelets from seismic sources. Standard deconvolution techniques to improve resolution of the data are limited by phase‐shifted waveforms such as supercritical reflections from the Moho as well as various refraction phases from the crust. We extended the minimum information deconvolution technique to the analytic formulation to deconvolve seismic data that contain phase‐shifted wavelets. When applied to a seismogram, analytic deconvolution yields the amplitude of the earth's impulse response function in the form of the signal envelope. An application example with ocean bottom seismometer data acquired along the Canadian continental margin of the Atlantic shows a remarkable improvement in resolution by compressing the wavelets effectively.

Comment on ‘Aspects of 1D seismic modelling using the Goupillaud principle’ by Evert Slob and Anton Ziolkowski1

AbstractThe paper by Slob and Ziolkowski (1993) is apparently a comment on my paper (Szaraniec 1984) on odd‐depth structure. In fact the basic understanding of a seismogram is in question. The fundamental equation for an odd‐depth model and its subsequent deconvolution is correct with no additional geological constraints. This is the essence of my reply which is contained in the following points. The discussion by Slob and Ziolkowski suffers from incoherence. On page 142 the Goupillaud (1961) paper is quoted:“… we must use a sampling rate at least double that… minimum interval…”. In the following analysis of such a postulated model Slob and Ziolkowski say that “… two constants are used in the model: Δt as sampling rate and 2Δt as two‐way traveltime”. By reversing the Goupillaud postulation all the subsequent criticism becomes unreliable for the real Goupillaud postulation as well as the odd‐depth model. Slob and Ziolkowski take into consideration what they call the total impulse response. This is over and above the demands of the fundamental property of an odd‐depth model.Following a similar approach I take truncated data in the form of a source function, S(z), convolved with a synthetic seismogram (earth impulse response), R̃(z), the free surface being included. The problem of data modelling is a crucial one and will be discussed in more detail below. By my reasoning, however, the function image may be considered as a mathematical construction introduced purely to work out the fundamental property. In this connection there is no question of this construction having a physical meaning.It is implicit that in terms of system theory, K(z) stands for what is known as input impedance. Our understandings of data are divergent but Slob and Ziolkowski state erroneously that:“Szaraniec (1984) gives (21) as the total impulse response…”. This point was not made.This inappropriate statement is repeated and echoed throughout the paper making the discussion by Slob and Ziolkowski, as well as the corrections proposed in their Appendix A, ineffective.Thus, my equation (2) is quoted in the form image which is in terms of the reflection response Gsc and holds true at least in mathematical terms. No wonder that “this identity is not valid for the total impulse response” (sic), which is denoted as G(z). None the less a substitution of G for Gsc is made in Appendix A, equation (A3). The equation numbers in my paper and in Appendix A are irrelevant, but (A3) is substituted for (32) (both numbers of equations from the authors’ paper). Afterwards, the mathematical incorrectness of the resulting equation is proved (which was already evident) and the final result (A16) is quite obviously different from my equation (2). However, the substitution in question is not my invention. With regard to the problem of data modelling, I consider a bi‐directional ID seismic source located just below the earth's surface. The downgoing unit impulse response is accompanied by a reflected upgoing unit impulse and the earth response is now doubled. The total impulse response for this model is thus given by image where (—r0) =— 1 stands for the surface reflection coefficient in an upward direction. Thus image that is to say, the total response to a unit excitation is identical with the input impedance as it must be in system theory.The one‐directional 1D seismic source model is in question. There must be a reaction to every action. When only the downgoing unit impulse of energy is considered, what about the compensation? In more realistic modelling, an early part of a total seismogram is unknown (absent) and the seismogram is seen in segments or through the windows. That is why in the usual approach, especially in dynamic deconvolution problems, synthetic data in the presence of the free surface are considered as an equivalent of the global reflection coefficient. It is implicit that model arises from a truncated total seismogram represented as a source function convolved with a truncated global reflection coefficient. Validation or invalidation of the truncation procedure for a numerically specified model may be attempted in the frame of the odd‐depth assumption. My equations (22) and (23) have been designed for investigating the absence or presence of truncated energy.The odd‐depth formalism allows the possibility of reconstructing an earlier part of a seismogram (Szaraniec 1984), that is to say, a numerical recovery of unknown moments which are unlikely designed by Slob and Ziolkowski for the data.

Wavelets, well logs and Wiener filters

The deconvolution of source wavelets from seismic traces can provide useful estimates of the earth &apos;s impulse response and thereby aid in geological interpretation. This signal analysis tool has been especially useful in the resolution of thin rock layers and has received widespread application since its development in the 1950s. The deconvolution problem often involves the estimation of a wavelet before removing it by digital filtering. Even in the noiseless case, this task is difficult since the seismic trace can be viewed as the convolution of an unknown wavelet coefficient sequence with an unknown earth impulse response. In cases where multiples have been eliminated, this impulse response becomes transformed into a sequence of reftection coefficients. As pointed out by Ziolkowski (1982), this deconvolution will yield non-unique solutions since it essentially involves solving one equation for two unknown sequences: the wavelet coefficients and the reflection coefficients. However, this inverse problem is not hopeless since the two sequences generally have different statistical properties. Moreover, well log data may be used to sort out the non-uniqueness. The Earth&apos;s impulse response derived from well log information and the statistical estimate of the source wavelet can be modified together by least squares inversion in order to model the seismic trace and provide reliable wavelet estimates. Recent wavelet estimation research has shown a similarity between wavelet estimates obtained by least squares inversion and those found by deconvolution of the earth&apos;s impulse response from the seismic trace. Upon examination of these procedures, it turns out that least squares inversion with a constrained earth&apos;s impulse response is equivalent to deriving the Wiener filter which shapes the earth&apos;s impulse response to the recorded seismic trace. When reflection coefficients derived from a well log are simply used as an initial guess for model, parameters, the unconstrained inversion procedure provides wavelet estimates which are similar to wavelet estimates obtained by deconvolution of the earth&apos;s impulse response. The similarities are demonstrated for synthetic and real data cases, and explained in heuristic terms.

DEBUBBLING: A GENERALIZED LINEAR INVERSE APPROACH*

AbstractOn seismograms recorded at sea bubble pulse oscillations can present a serious problem to an interpreter. We propose a new approach, based on generalized linear inverse theory, to the solution of the debubbling problem. Under the usual assumption that a seismogram can be modelled as the convolution of the earth's impulse response and a source wavelet we show that estimation of either the wavelet or the impulse response can be formulated as a generalized linear inverse problem. This parametric approach involves solution of a system of equations by minimizing the error vector (ΔX = Xobs– Xcal) in a least squares sense. One of the most significant results is that the method enables us to control the accuracy of the solution so that it is consistent with the observational errors and/or known noise levels.The complete debubbling procedure can be described in four steps: (1) apply minimum entropy deconvolution to the observed data to obtain a deconvolved spike trace, a first approximation to the earth's response function; (2) use this trace and the observed data as input for the generalized linear inverse procedure to compute an estimated basic bubble pulse wavelet; (3) use the results of steps 1 and 2 to construct the compound source signature consisting of the primary pulse plus appropriate bubble oscillations; and (4) use the compound source signature and the observed data as input for the generalized linear inverse method to determine the estimated earth impulse response—a debubbled, deconvolved seismogram. We illustrate the applicability of the new approach with a set of synthetic seismic traces and with a set of field seismograms.A disadvantage of the procedure is that it is computationally expensive. Thus it may be more appropriate to apply the technique in cases where standard analysis techniques do not give acceptable results. In such cases the inherent advantages of the method may be exploited to provide better quality seismograms.