Seismic Exploration and Signal Processing

Linearized inversion norms

In linearized inversion, the model \mathbf{m} is estimated by minimizing a measure of the residual \mathbf{r}=\mathbf{d}-\mathbf{G}\mathbf{m} between observed and predicted data. The choice of norm determines how strongly each residual influences the solution. The L_p norm is \|\mathbf{r}\|_p=(\sum_i |r_i|^p)^{1/p}, and the L_\infty norm is the largest absolute residual, \max_i |r_i|. The L_2 norm penalizes each residual by r^2, so large deviations are weighted quadratically, while the L_1 norm penalizes by |r|, which grows only linearly. For this reason the L_2 norm gives higher weight to outliers than the L_1 norm, and the L_1 norm is the more robust choice when the data contain spikes or bad points. The L_\infty norm is determined only by the single largest residual, so the point with the maximum deviation receives all the weight while the others are effectively ignored.

The type of system governs whether an exact solution exists. For a linearly independent system, an underdetermined system (fewer equations than unknowns) is always consistent and has infinitely many solutions, including the trivial one for homogeneous systems. An overdetermined system (more equations than unknowns) generally has no exact solution and is solved in the least-squares sense. A square system of independent equations is even-determined and has a unique solution. An ill-conditioned system, by contrast, amplifies small data errors, so its solution is unstable and cannot be trusted in the presence of noise.

Damped least squares minimizes \|\mathbf{d}-\mathbf{G}\mathbf{m}\|^2+\epsilon^2\|\mathbf{m}\|^2, which gives \mathbf{m}=(\mathbf{G}^T\mathbf{G}+\epsilon^2\mathbf{I})^{-1}\mathbf{G}^T\mathbf{d}. The damping term keeps the normal matrix well conditioned by suppressing the contribution of small singular values, so the damping parameter stabilizes the solution. The price is reduced resolution, because the estimate is biased toward smaller models. Damping does not remove the non-uniqueness that comes from the null space of \mathbf{G}. More generally, regularization adds prior information as a penalty term, which constrains an ill-posed problem and selects a physically plausible solution. It controls the trade-off between data misfit and model size, and it never converts a well-posed problem into an ill-posed one.

The generalized inverse \mathbf{G}^{-g} maps data to model estimates, \mathbf{m}_{est}=\mathbf{G}^{-g}\mathbf{d}. Substituting \mathbf{d}=\mathbf{G}\mathbf{m}_{true} gives the model resolution matrix \mathbf{R}_m=\mathbf{G}^{-g}\mathbf{G}, which is M\times M and describes how well each model parameter is resolved. The product \mathbf{G}\mathbf{G}^{-g} is the N\times N data resolution matrix, and it should not be confused with \mathbf{R}_m. For a full-rank \mathbf{G}=\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^T, the pseudoinverse is obtained by inverting the nonzero singular values and transposing the orthogonal factors, giving \mathbf{G}^{+}=\mathbf{V}\boldsymbol{\Sigma}^{-1}\mathbf{U}^T. For a square full-rank matrix \mathbf{U}\mathbf{U}^T=\mathbf{I}, so equivalent forms reduce to the same expression. For an underdetermined system with full row rank, the minimum-norm generalized inverse is \mathbf{G}^T(\mathbf{G}\mathbf{G}^T)^{-1}, whereas (\mathbf{G}^T\mathbf{G})^{-1}\mathbf{G}^T is the form for overdetermined systems.

The model covariance matrix has the variances of the estimated parameters on its diagonal. The standard deviation of a parameter is the square root of its diagonal element, \sigma_{m_i}=\sqrt{C_{ii}}, and it is expressed in the units of that parameter.

Filter transfer function and Hilbert transform

For a linear time-invariant filter, the transfer function is the ratio of output to input in the frequency domain, H(\omega)=Y(\omega)/X(\omega), and it equals the Fourier transform of the impulse response. Under the convention F(\omega)=\int f(t)e^{-i\omega t}\,dt, a causal decaying exponential e^{-\alpha t}u(t) has transform 1/(\alpha+i\omega). Taking the ratio of two such decays with constants \alpha and \beta gives H(\omega)=(\alpha+i\omega)/(\beta+i\omega). Always check the sign convention first, because the same transform written with e^{+i\omega t} produces the conjugate expression.

The poles of a transfer function are the roots of its denominator polynomial, and the zeros are the roots of its numerator. For example, s^2+3s+2=(s+1)(s+2) has poles at s=-1 and s=-2, and a zero in the numerator does not change these poles.

In digital filters, the terminology separates two classes. A non-recursive (FIR) filter has a transfer function that is a polynomial in z^{-1}, so it has no poles except trivial ones at the origin. A recursive (IIR) filter feeds back past outputs, so its transfer function may have poles alone or both poles and zeros. The transfer function of any digital filter is defined as H(z)=Y(z)/X(z), the ratio of the output to the input. The order of a digital filter counts both the past inputs and the past outputs that it uses, so it is not limited to inputs.

The Hilbert transform is a linear, shift-invariant operation that can be written as convolution with 1/(\pi t), so \hat f(t)=\frac{1}{\pi}\,\mathrm{PV}\!\int f(\tau)/(t-\tau)\,d\tau. The Hilbert transform of the Dirac delta is therefore the kernel 1/(\pi t), with the sign reversed under the opposite convention, -1/(\pi t). This kernel has two branches of opposite sign on either side of the origin, it is singular at t=0, and its magnitude decays as 1/|t|. In the frequency domain the Hilbert transform multiplies the spectrum by -i\,\mathrm{sgn}(\omega), so applying it twice multiplies by (-i\,\mathrm{sgn}\,\omega)^2=-1 for \omega\neq 0. Hence the Hilbert transform of \hat f(t) equals -f(t), a sign change rather than the original signal.

The analytic (complex) signal is f_C=f+i\hat f, with the original signal as the real part. For \sin t, the Hilbert transform is -\cos t, so f_C=\sin t-i\cos t. For \cos\theta, the Hilbert transform is \sin\theta, so the analytic signal is A(\theta)=\cos\theta+i\sin\theta=e^{i\theta}.

2D seismic stacking fold and S/N

In a 2D survey, the fold is the number of traces that sample the same subsurface common midpoint (CMP). The fold is therefore defined by the maximum number of traces in a CMP gather, and these traces are later stacked. Common shot, common receiver and common offset gathers group traces by a different shared parameter, so their trace counts do not define the fold.

When all sources in a group are recorded by all receivers in another group, the number of contributing traces is the product of the two group sizes, N=N_s\times N_r. Random noise is uncorrelated between traces, while the reflection signal is coherent, so stacking N traces improves the signal-to-noise ratio by a factor of \sqrt{N}.

In a stacking chart, traces are plotted with the shot, receiver, mid-point and offset as coordinates. Traces that share a fixed value of one coordinate lie on a line parallel to the axis of the coordinate that varies along it. Lines of constant CMP run parallel to the offset axis, and lines of constant offset run parallel to the midpoint axis. Constant-shot traces lie on a line parallel to the receiver axis, and constant-receiver traces lie on a line parallel to the shot axis.

A reflected wavefront from a concave reflector focuses energy at a common point. The incident wavefront from a surface source is planar, with zero curvature, while a reflector with finite curvature has positive curvature, so the reflector curvature is greater than that of the incident wavefront. This curvature mismatch gives the focusing effect that concentrates reflected energy at the point.

The number of traces that can be assigned to a common midpoint depends on the geometry. The CMP interval is half the receiver interval, and each CMP bin corresponds to one zero-offset trace. The total number of zero-offset traces is the horizontal length covered by the survey divided by the CMP interval. That covered length is the span of the shot positions plus the near offset and the span of the receiver spread.

Seismic NMO and offset traveltime

For a reflector at depth z and offset x, the two-way reflection traveltime is t_R=\sqrt{4z^2+x^2}/v, and its normal moveout relative to the zero-offset time t_0=2z/v is \Delta t_R=(\sqrt{4z^2+x^2}-2z)/v. For small offsets this is approximately x^2/(4zv).

A diffraction from a point scatterer has a larger moveout than a reflection at the same offset. If the shot lies directly above the diffractor at depth z, the diffraction path has a vertical leg of length z and an inclined leg of length \sqrt{z^2+x^2}, so its moveout is \Delta t_D=(\sqrt{z^2+x^2}-z)/v\approx x^2/(2zv). The diffraction moveout is therefore twice the reflection moveout, and this ratio is useful when comparing the two events.

In a homogeneous medium the first arrival of a P-wave is the direct straight ray, so its traveltime is the source-to-receiver distance divided by the velocity, t=d/v. The difference in first-arrival time between two receivers is therefore the difference of their distances from the source divided by v. This calculation is made by writing each distance with the Pythagorean theorem from the horizontal and vertical separations.

In the X^2–T^2 formulation, the hyperbolic moveout of a horizontally layered medium is extended by a fourth-order term in X. The third term in the expansion, proportional to X^4, arises from the vertical variation of velocity, so it describes a heterogeneous medium. A homogeneous medium gives an exactly hyperbolic moveout with no X^4 term. The same fourth-order correction is relevant before AVO analysis, where the amplitudes must be corrected for geometrical spreading.

A multiple that reflects from a horizontal interface and the free surface follows a zigzag ray path. Each leg contributes its vertical extent, so the total vertical distance is the sum of the vertical legs. The traveltime is then obtained by unfolding the path into a single straight line with that total vertical extent and the horizontal offset, giving t=\sqrt{x^2+(\Sigma z)^2}/v.

Fourier transform pairs

The Fourier pair under the convention g(t)=\frac{1}{2\pi}\int G(\omega)e^{j\omega t}\,d\omega and G(\omega)=\int g(t)e^{-j\omega t}\,dt has a derivative property. Differentiating the inverse transform with respect to t gives \frac{dg}{dt}=\frac{1}{2\pi}\int j\omega\,G(\omega)e^{j\omega t}\,d\omega, so the transform of dg/dt is j\omega G(\omega). Differentiation in the time domain therefore multiplies the spectrum by j\omega, and the sign and the factor j must match the convention used.

The convolution of two rectangular (box-car) pulses is a piecewise linear function. When the two pulses have equal durations the result is a triangle. When the durations are unequal and finite, both the rising and falling edges are linear ramps joined by a flat top, so the result is a trapezoid.

Convolution in time corresponds to multiplication of spectra, so the phase spectra add. Cross-correlation corresponds to X_1(f)X_2^*(f), whose phase is \phi_1(f)-\phi_2(f), so the phase spectra are subtracted. Deconvolution in the time domain corresponds to dividing the spectra, X_1(f)/X_2(f), which also subtracts the phase spectra. Subtraction of two time series subtracts the signals linearly and does not subtract their phase spectra.

2D seismic processing sequence

A 2D reflection dataset is conditioned in a fixed order before it is converted from time to depth. Deconvolution is applied first to shape the source wavelet and suppress multiples. Filtering removes unwanted frequency content, and equalization balances the amplitudes. Migration then repositions dipping reflectors to their true subsurface locations, and coherency analysis is used last to assess the migrated data before time–depth conversion. Migration placed before the amplitude and wavelet conditioning steps reverses this workflow.

In a 2D stacked section, the vertical axis is two-way travel time. The horizontal axis is the common midpoint (CMP) location, because the stacking process combines all traces that share a midpoint into one trace. Receiver locations, source locations and offsets are removed or combined during stacking, so they cannot serve as the horizontal axis.

Weathered layer static correction

The weathered layer is a near-surface zone of low velocity, and it delays the arrivals. A static correction removes these delays by referencing all traces to a datum, so each trace is shifted by the one-way traveltime through the weathered layer between its source or receiver and the datum. For a homogeneous weathered layer of thickness h and velocity V_w, this one-way time is t=h/V_w.

If a source or receiver lies at the base of the weathered layer or at the datum, the path to the datum has no weathered-layer leg, so its static correction is zero. When a receiver is above the datum, the path within the weathered layer contributes the static, and its magnitude is the thickness divided by the weathered-layer velocity.

The first Fresnel zone of a reflection is the area on the reflector from which energy arrives within half a wavelength of the zero-offset path. For a zero-offset reflection from depth z in a medium with wavelength \lambda=V/f, the radius is r_F=\sqrt{\lambda z/2}, which is equivalent to r_F=(V/2)\sqrt{t/f} with t the two-way time. The datum does not change the radius, so it is evaluated at the reflector depth below the datum.

Head wave moveout in layered earth

A head wave is a refracted arrival that travels along the top of a higher-velocity layer below the surface and returns to the surface at the critical angle. In horizontal, isotropic and homogeneous layers, its traveltime is t=x/v_2+t_i, where x is the offset, v_2 is the velocity of the refracting layer and t_i is the intercept time. This relation is linear in offset, so the head-wave moveout appears as a straight line on a time–distance plot.

For a two-layer model with velocities v_1<v_2, the critical angle satisfies \sin\theta_c=v_1/v_2. The critical distance is the surface distance at which the refracted arrival first emerges, and for a first layer of thickness h it is x_c=2h\tan\theta_c. The critical distance is computed from the critical angle and the layer thickness, so it depends on both velocities and the depth of the interface.

Normal incidence reflection coefficient

The acoustic impedance of a layer is Z = \rho V, the product of density and P-wave velocity. For a wave hitting an interface at normal incidence, the reflection coefficient is R = (Z_2 - Z_1)/(Z_2 + Z_1), where Z_1 is the impedance of the incident layer and Z_2 that of the layer below. The sign of R tells you the polarity of the reflection: a positive R keeps the polarity of the incident wave, and a negative R reverses it. Only the impedance contrast matters, so R is dimensionless and unchanged if both impedances are scaled by the same factor. A common mistake is to use V or \rho alone instead of their product. Example: if Z_1 = 4.0\times10^6 and Z_2 = 6.0\times10^6 kg/(m^2s), then R = 2/10 = 0.20.

Z-transform of discrete sequences

The Z-transform of a sequence x[n] is X(z) = \sum_n x[n]\,z^{n} in the form used by some GATE items, or \sum_n x[n]\,z^{-n} in the standard textbook form. Always read the convention given in the question, because the same sequence gives mirrored powers under the two definitions. To transform a finite sequence, write one power of z for each nonzero sample at its index, and drop the zero-valued samples. A delay of one sample multiplies the transform by z^{-1} in the standard form. For a causal recursive filter y[n] = a\,x[n] + b\,x[n-1] + c\,y[n-1] with the system at rest before the first sample, the outputs follow by iteration from n=0. Example: for y_n = x_n + 0.5\,y_{n-1} with x_0=2, x_1=4, x_2=0, we get y_0=2, y_1=5, and y_2=2.5.