# User:Ageary/Chapter 9

Series Problems-in-Exploration-Seismology-and-their-Solutions.jpg Geophysical References Series Problems in Exploration Seismology and their Solutions Lloyd P. Geldart and Robert E. Sheriff http://dx.doi.org/10.1190/1.9781560801733 ISBN 9781560801153 SEG Online Store

## 9.1 Fourier series

9.1a Show that the Fourier series coefficients, ${\displaystyle a_{n}}$ and $\displaystyle b_{n}$ in equation (9.1a), are given by equations (9.1d,e).

Background

If ${\displaystyle g(t)}$ is a periodic function, that is, one that repeats exactly at intervals ${\displaystyle T}$ (called the period) and has a finite number of maxima, minima, and discontinuities in each interval ${\displaystyle T}$, then it can be expressed in a Fourier series in any one of the three following equivalent forms, ${\displaystyle n}$ being an integer:

 {\displaystyle {\begin{aligned}g(t)={\frac {1}{2}}a_{0}+\sum \limits _{n=1}^{+\infty }(a_{n}\cos \omega _{n}t+b_{n}\sin \omega _{n}t)\end{aligned}}} (9.1a)
 {\displaystyle {\begin{aligned}={\frac {1}{2}}c_{0}+\sum \limits _{n=1}^{+\infty }c_{n}\cos(\omega _{n}t-\gamma _{n})\end{aligned}}} (9.1b)
 \displaystyle \begin{align} =\sum\limits_{n=-\infty}^{+\infty} \alpha _{n} e^{\hbox {j}\omega _{n} t}, \end{align} (9.1c)

where $\displaystyle \omega _{0} =2\pi f_{0} =2\pi /T=$ fundamental angular frequency, ${\displaystyle \omega _{n}=n\omega _{0}=}$ harmonic frequency, and the amplitudes and phases are given by the following equations where ${\displaystyle 0, except for equation (9.1h) where ${\displaystyle -\infty \leq n\leq +\infty }$:

 {\displaystyle {\begin{aligned}a_{n}={\frac {2}{T}}\int _{-T/2}^{+T/2}g(t)\cos \,\omega _{n}t\,\mathrm {d} t,\end{aligned}}} (9.1d)
 {\displaystyle {\begin{aligned}b_{n}={\frac {2}{T}}\int _{-T/2}^{+T/2}g(t)\sin \,\omega _{n}t\,\mathrm {d} t,\end{aligned}}} (9.1e)
 {\displaystyle {\begin{aligned}c_{n}={\frac {2}{T}}\int _{-T/2}^{+T/2}g(t)\cos \,(\omega _{n}t-\gamma _{n})\,\mathrm {d} t,\end{aligned}}} (9.1f)
 {\displaystyle {\begin{aligned}\gamma _{n}=\tan ^{-1}(b_{n}/a_{n})=\mathrm {phase} ,\gamma _{0}=0,\end{aligned}}} (9.1g)
 \displaystyle \begin{align} \alpha _{n} =\frac{1}{T} \int_{-T /2}^{+T /2}g(t)e^{-{\hbox {j}}\omega _{n}t} \mathrm{d}t. \end{align} (9.1h)

Solution

We multiply both sides of equation (9.1a) by $\displaystyle {\rm \; cos\;}\omega _{r} t$ , ${\displaystyle r}$ being any nonzero positive integer, and integrate between the limits $\displaystyle -T/2$ and ${\displaystyle +T/2}$. (We omit the limits on the integrals in the following proofs since they are all the same.) Equation (9.1a) now becomes

{\displaystyle {\begin{aligned}\int g(t)\cos \omega _{r}t\mathrm {d} t&=\sum \limits _{n=1}^{+\infty }a_{n}\int \cos \,\omega _{r}t\cos \omega _{n}t\,\mathrm {d} t\\&+\sum \limits _{n=1}^{+\infty }b_{n}\int \cos \omega _{r}t\sin \omega _{n}t\,\mathrm {d} t.\end{aligned}}}

By using the identities

{\displaystyle {\begin{aligned}\cos x\cos y={\frac {1}{2}}[\cos(x+y)+\cos(x-y)]\end{aligned}}}

{\displaystyle {\begin{aligned}{\hbox{and}}\quad \quad \quad \cos x\sin y={\frac {1}{2}}[\sin(x+y)-\sin(x-y)],\end{aligned}}}

the first integrand reduces to the sum of two cosines of the form ${\displaystyle \cos(\omega _{r}\pm \omega _{n})t}$ while the second integrand becomes the difference between two sines of the form ${\displaystyle \sin(\omega _{r}\pm \omega _{n})t.}$ If $\displaystyle r\ne n$ , the integrals are zero, because at the limits the arguments of the cosines and sines are equal to $\displaystyle 2p\omega _{0} T/2=2p\pi$ , ${\displaystyle p}$ integral. When ${\displaystyle r=n}$, the first integrand becomes $\displaystyle \cos^{2} \omega _{n} t=\frac{1}{2} (1+\cos\,2\omega _{n} t)$ and the right-hand side equals $\displaystyle a_{n} T/2$ . When ${\displaystyle r=n}$, the second integrand becomes $\displaystyle \cos \omega _{n} t\sin\omega _{n} t\,\mathrm{d}t=\frac{1}{2} \sin\,2\omega _{n}t\, \mathrm{d}t$ , which again gives zero upon integration.

Thus we are left with

{\displaystyle {\begin{aligned}\int g(t)\cos \omega _{n}t\,\mathrm {d} t=a_{n}T/2,\quad n=0,1,2,\ldots .\end{aligned}}}

Solving for ${\displaystyle a_{n}}$ gives equation (9.1d).

To verify equation (9.1e), we multiply both sides of equation (9.1a) by ${\displaystyle \sin \omega _{r}t}$ and integrate, all integrals on the right side of the equation giving zero except for the one with the integrand ${\displaystyle b_{n}\sin ^{2}\omega _{n}t=(b_{n}/2)(1-\cos \,2\omega _{n}t)}$. Proceeding as before, we arrive at equation (9.1e).

When we set $\displaystyle n=0$ in equations (9.1d,e) we get ${\displaystyle b_{0}=0}$ and

{\displaystyle {\begin{aligned}a_{0}/2={\frac {1}{T}}\int g(t)\,\mathrm {d} t=\mathrm {average\ value\ of} \,g(t).\end{aligned}}}

9.1b Verify equations (9.1f,g) for the Fourier series coefficients $\displaystyle c_{n}$ and $\displaystyle \gamma _{n}$ .

Solution

To verify equation (9.1f) for ${\displaystyle c_{n}}$, we multiply both sides of equation (9.1b) by $\displaystyle \cos(\omega_{r} t-\gamma _{r})$ , ${\displaystyle r\neq 0}$, and proceed as in part (a). All integrals vanish except when ${\displaystyle r=n}$, and the result is equation (9.1f).

When ${\displaystyle r=0}$, we have

\displaystyle \begin{align} \int g(t)\,\mathrm{d}t=\frac{c_{0}}{2} \int \mathrm{d}t=\frac{c_{0} T}{2} ,\\ \hbox{so}\quad \quad \quad \frac{c_{0}}{2} =\frac{1}{T} \int g(t)\,\mathrm{d}t = \mathrm{average\ value\ of}\ g(t). \end{align}

To verify equation (9.1g) for ${\displaystyle \gamma _{n}}$, we compare terms in equations (9.1a) and (9.1b) that involve ${\displaystyle \omega _{n}t}$, ${\displaystyle n\neq 0}$. The result is

\displaystyle \begin{align} a_{n}\cos \omega_{n} t+b_{n} \sin\omega _{n} t&=c_{n} \cos(\omega _{n} t-\gamma _{n})\\ &=c_{n}(\cos \omega_{n}t\,\cos\,\gamma_{n}+\sin \omega_{n}t\,\sin\, \gamma_{n}). \end{align}

Equating coefficients of ${\displaystyle \cos \omega _{n}t}$ and ${\displaystyle \sin \omega _{n}t}$, we get

{\displaystyle {\begin{aligned}a_{n}=c_{n}\cos \gamma _{n},\quad b_{n}=c_{n}\sin \gamma _{n},\quad c_{n}^{2}=a_{n}^{2}+b_{n}^{2},\quad \tan \gamma _{n}=b_{n}/a_{n}.\end{aligned}}}

When ${\displaystyle n=0}$, ${\displaystyle a_{0}=c_{0}}$, ${\displaystyle b_{0}=0}$, ${\displaystyle \gamma _{0}=0}$.

9.1c Verify the coefficients for the exponential form, equation (9.1h).

Solution

To derive equation (9.1h), we use Euler’s equations (Sheriff and Geldart, 1995, problem 15.12a), namely,

{\displaystyle {\begin{aligned}\cos x={\frac {1}{2}}(e^{{\hbox{j}}x}+e^{-{\hbox{j}}x}),\quad \sin \,x={\frac {1}{2{\hbox{j}}}}(e^{{\hbox{j}}x}-e^{-{\hbox{j}}x}),\\{\hbox{or}}\quad \quad \quad e^{{\hbox{j}}x}=(\cos x+{\hbox{j}}\sin x),\quad e^{-{\hbox{j}}x}=(\cos x-{\hbox{j}}\sin x).\end{aligned}}}

Multiplying both sides of equation (9.1c) by $\displaystyle e^{-\hbox {j}\omega _{r} t}$ , ${\displaystyle r\neq 0}$, and integrating, we get a series of integrals with integrands such as ${\displaystyle e^{{\hbox{j}}(\omega _{n}-\omega _{r})t}}$. Using Euler’s formula, ${\displaystyle e^{{\hbox{j}}(\omega _{n}-\omega _{r})t}={\hbox{cos}}(\omega _{n}-\omega _{r})t+{\hbox{j}}\ {\hbox{sin}}(\omega _{n}-\omega _{r})t}$, so the integrals vanish for all values of ${\displaystyle n}$ except when ${\displaystyle n=r}$; for this value we get

{\displaystyle {\begin{aligned}\int g(t)e^{-{\hbox{j}}\omega _{n}t}\,\mathrm {\rm {d}} t=\alpha _{n}\int \mathrm {d} t=\alpha _{n}T,\end{aligned}}}

which is equation (9.1h), ${\displaystyle n\neq 0}$. For ${\displaystyle n=0}$,

{\displaystyle {\begin{aligned}\alpha _{0}={\frac {1}{T}}\int g(t)\,\mathrm {d} t=\mathrm {average\ value\ of} \,g(t).\end{aligned}}}

## 9.2 Space-domain convolution; Vibroseis acquisition

The techniques and concepts of convolution, aliasing, ${\displaystyle z}$-transforms, and so on can be applied to other domains than the time-domain. Express the source and group patterns of Figure 9.2a as functions of ${\displaystyle x}$ (horizontal coordinate) and convolve the two to verify the effective pattern shown.

Background

The convolution of two functions ${\displaystyle g(t)}$ and $\displaystyle h(t)$ , written ${\displaystyle g(t)*h(t)}$, is defined as

 {\displaystyle {\begin{aligned}g(t)*h(t)&=\int \limits _{-\infty }^{+\infty }g(\tau )h(t-\tau )\,\mathrm {d} {\tau }\\&=\int \limits _{-\infty }^{+\infty }h(\tau )g(t-\tau )\,\mathrm {d} {\tau }=h(t)*g(t).\end{aligned}}} (9.2a)

For digital functions (see problem 9.4), the integrals become sums:

 {\displaystyle {\begin{aligned}g_{t}*h_{t}=\sum \limits _{k}^{}g_{k}\ h_{t-k}=\sum \limits _{k}^{}h_{k}\ g_{t-k}=g_{t}*g_{t}.\end{aligned}}} (9.2b)

These equations state that convolution is equivalent to superposition in which each element of one function is replaced by the other function multiplied (weighted) by the element being replaced. The sum of all values at a given time is the value of ${\displaystyle g(t)*h(t)}$ at that instant.

Figure 9.2a.  Vibroseis acquisition. (i) Field locations for six sweeps (shown displaced vertically); (ii), array pattern for four vibroseis trucks spaced 33 m apart moving forward 16.5 m for each sweep. Each receiver array consists of geophones spaced uniformly over 100 m. For the next location, an additional group of phones is added on one side and dropped on the other. The numbers in (ii) give the multiplicity.
Figure 9.2a.  Boxcar and its transform.
Figure 9.2b.  Two boxcars.

Another explanation of convolution is the following: We note that ${\displaystyle g(-t)}$ is $\displaystyle g(t)$ reflected in the ${\displaystyle t}$-axis; in other words, the curve of ${\displaystyle g(-t)}$ is the same as that of ${\displaystyle g(t)}$ except that it is reversed in direction [keeping ${\displaystyle g(0)}$ fixed]. The function ${\displaystyle g(t-\tau )}$ is $\displaystyle g(t)$ displaced ${\displaystyle \tau }$ units to the right. Thus, in equation (9.2a) the value of ${\displaystyle g(t)*h(t)}$ at the time ${\displaystyle t=\tau }$ is obtained by moving ${\displaystyle h(t)}$ to the right ${\displaystyle \tau }$ units, reflecting it in the ${\displaystyle t}$-axis, then summing (integrating) the products of corresponding abscissas. The result is the same whichever function is displaced and reflected.

Arrays are discussed in problem 8.6, aliasing in problem 9.4, and ${\displaystyle z}$-transforms in Sheriff and Geldart, 1995, section 15.5.3.

Solution

Taking the spatial sampling interval as 16.5 m, the source pattern consists of [1,1,2,2,3,3,3, 3,2,2,1,1] for a total array length of 200 m. The geophone group is an array 100 m long; to convolve the source and receiver arrays, they should have the same spatial intervals, so we take six receivers spaced 16.5 m apart. To convolve the two arrays, we replace each element of the receiver array with the source array (taking the weights as unity). The result is that shown at the bottom of Figure 9.2a.

## 9.3 Fourier transforms of the unit impulse and boxcar

9.3a Because a unit impulse ${\displaystyle \delta ({\textit {t}})}$ is zero except at ${\displaystyle t=0}$, where it equals $\displaystyle +1$ , we can apply the Fourier transform equation (9.3c) and find that ${\displaystyle \delta ({\textit {t}})\leftrightarrow }$ ${\displaystyle +1}$. Show that

{\displaystyle {\begin{aligned}\delta (t-t_{0})\leftrightarrow e^{-{\hbox{j}}\varpi t_{0}}.\end{aligned}}}

Background

Equations (9.1a,b,c) apply to a harmonic function which repeats after every interval $\displaystyle T$ . If we let ${\displaystyle T}$ increase, the repetitions occur at longer intervals and in the limit when ${\displaystyle T=\infty ,}$ the function becomes aperiodic, that is, it no longer repeats. To see the effect of this on equation (9.1c), we replace $\displaystyle \alpha _{n}$ in equation (9.1c) with the right-hand side of equation (9.1h). The result is

 \displaystyle \begin{align} g(t)=\sum\limits_{n=-\infty}^{+\infty} \left(\frac{1}{T}\int_{-T /2}^{+T /2} g(t)e^{-\hbox {j}\omega _{n}t} \mathrm{d}t\right)e^{\hbox {j}\omega_{n} t}. \end{align} (9.3a)

We now let $\displaystyle T\to +\infty$ so that $\displaystyle \omega _{0} \to 0$ , which causes the two adjacent frequencies, ${\displaystyle \omega _{n}}$ and ${\displaystyle \omega _{n+1}}$, differing by ${\displaystyle \omega _{0}}$, to approach each other. In the limit ${\displaystyle \omega _{n}}$ becomes a continuous variable ${\displaystyle \omega }$, $\displaystyle 1/T$ becomes ${\displaystyle {\hbox{d}}\omega /2\pi }$, the summation becomes an integral, and equation (9.3a) becomes

 {\displaystyle {\begin{aligned}g(t)={\frac {1}{2\pi }}\int \limits _{-\infty }^{+\infty }e^{+{\hbox{j}}\omega t}\left[\int \limits _{-\infty }^{+\infty }g(t)e^{-{\hbox{j}}\omega t}{\hbox{d}}t\right]\mathrm {d} \omega .\end{aligned}}} (9.3b)

It is convenient to represent the inner integral by a symbol and write equation (9.3b) as two equations:

 {\displaystyle {\begin{aligned}G(f)=\int \limits _{-\infty }^{+\infty }g(t)e^{-{\hbox{j}}2\pi ft}\,\mathrm {d} t\quad \mathrm {or} \quad G(\omega )=\int \limits _{-\infty }^{+\infty }g(t)e^{-{\hbox{j}}\omega t}\mathrm {d} t,\end{aligned}}} (9.3c)

 {\displaystyle {\begin{aligned}g(t)=\int \limits _{-\infty }^{+\infty }G(f)e^{{\hbox{j}}2\pi ft}\mathrm {d} f\quad \mathrm {or} \quad g(t)={\frac {1}{2\pi }}\int \limits _{-\infty }^{+\infty }G(\omega )e^{{\hbox{j}}\omega t}\mathrm {d} \omega .\end{aligned}}} (9.3d)

The function ${\displaystyle G(\omega )}$ [or ${\displaystyle G(f)}$] is the Fourier transform of ${\displaystyle g(t)}$ while ${\displaystyle g(t)}$ is the inverse Fourier transform of ${\displaystyle G(\omega )}$ [or $\displaystyle G(f)$ ]. The relation between ${\displaystyle g(t)}$ and ${\displaystyle G(\omega )}$ can be indicated by a double arrow:

 \displaystyle \begin{align} g(t)\leftrightarrow G(\omega)\quad \mathrm{or}\quad g(t)\leftrightarrow G(f). \end{align} (9.3e)

The unit impulse $\displaystyle \delta(t)$ , also called the Dirac delta, is by definition zero everywhere except at ${\displaystyle t=0}$ where it equals $\displaystyle +1$ ; similarly, $\displaystyle \delta(t-t_{0})$ is zero except when ${\displaystyle (t-t_{0})=0}$, where it equals ${\displaystyle +1}$. In digital notation (see problem 9.4) we write ${\displaystyle \delta _{t}}$, $\displaystyle \delta_{(t-t_{0})}$ , but the meaning is the same.

The convolution ${\displaystyle g(t)*h(t)}$ is discussed in problem 9.2. The convolution theorem [see Sheriff and Geldart, 1995, equation (15.145)] states that

 {\displaystyle {\begin{aligned}g(t)*h(t)\leftrightarrow G(\omega )H(\omega ),\end{aligned}}} (9.3f)

where ${\displaystyle G(\omega )}$ and $\displaystyle H(\omega)$ are the Fourier transforms of ${\displaystyle g(t)}$ and ${\displaystyle h(t)}$. For digital functions we use ${\displaystyle z}$-transforms (see Sheriff and Geldart, 1995, section 15.5.3), the equivalent of equation (9.3f) being

 \displaystyle \begin{align} g_{t} *h_{t} \leftrightarrow G(z)H(z). \end{align} (9.3g)

The boxcar, written ${\displaystyle {\hbox{box}}_{2t_{0}}(t)}$, is a function whose value is everywhere zero except in the interval from $\displaystyle -t_{0}$ to ${\displaystyle +t_{0}}$, where its value is ${\displaystyle +1}$. A boxcar in the frequency domain is written ${\displaystyle {\hbox{box}}_{2f_{0}}(f)}$ or ${\displaystyle {\hbox{box}}_{2\omega _{0}}(\omega )}$.

Solution

The transform of ${\displaystyle \delta (t-t_{0})}$ is given by equation (9.3c):

 {\displaystyle {\begin{aligned}\delta (t-t_{0})\leftrightarrow \int \limits _{-\infty }^{+\infty }\delta (t-t_{0})e^{-{\hbox{j}}\omega t}\mathrm {d} t=e^{-{\hbox{j}}\omega t_{0}}.\end{aligned}}} (9.3h)

9.3b Show that

 {\displaystyle {\begin{aligned}\delta (t)*g(t)=g(t),\quad \delta _{t}*g_{t}=g_{t},\end{aligned}}} (9.3i)
 {\displaystyle {\begin{aligned}\delta (t-\tau )*g(t)=g(t-\tau ),\quad \delta _{t-\tau }*g_{t}=g_{t-\tau }.\end{aligned}}} (9.3j)

Solution

Convolution involves replacing each element of one function with the other function and since ${\displaystyle \delta (t)}$ involves only one nonzero element at zero time, replacing an element at ${\displaystyle t=0}$ with ${\displaystyle g(t)}$ gives $\displaystyle g(t)$ , thus proving equation (9.3i). Likewise, replacing an element at time ${\displaystyle \tau }$ with $\displaystyle g(t)$ gives ${\displaystyle g(t-\tau )}$, thus proving equation (9.3j).

An alternative proof for both cases above can be obtained by using transforms. Because the transform of both ${\displaystyle \delta (t)}$ and $\displaystyle \delta_{t}$ is ${\displaystyle +1}$, equations (9.3f, g) show that the transforms of $\displaystyle \delta(t)*g(t)$ and ${\displaystyle \delta _{t}*g_{t}}$ give ${\displaystyle G(\omega )}$ and ${\displaystyle G(z)}$, respectively.

To prove equations (9.3j), we first use equation (9.3h) and get

{\displaystyle {\begin{aligned}\delta (t-\tau )*g(t)\leftrightarrow e^{-{\hbox{j}}\omega \tau }G(\omega )\leftrightarrow g(t-\tau ),\end{aligned}}}

using Sheriff and Geldart, 1995, equation (15.136). Next we use Sheriff and Geldart, 1995, problem (15.27) to write

\displaystyle \begin{align} \delta_{t-\tau} *g_{t} \leftrightarrow z^{\tau} G(z)\leftrightarrow g_{t-\tau}. \end{align}

9.3c Show that a boxcar of height ${\displaystyle h}$ extending from ${\displaystyle -f_{0}}$ to $\displaystyle +f_{0}$ in the frequency domain has the transform

{\displaystyle {\begin{aligned}h\ \mathrm {box} _{2f_{0}}(f)\leftrightarrow A{\frac {\sin(2\pi f_{0}t)}{2\pi f_{0}t}}=A\sin {\hbox{c}}(2\pi f_{0}t),\end{aligned}}}

where ${\displaystyle A=2hf_{0}=}$ area of the boxcar.

Solution

Using equation (9.3d) a boxcar in the frequency domain becomes ${\displaystyle g(t)}$ in the time domain:

 \displaystyle \begin{align} g(t)&=\int_{-\infty}^{+\infty} h \mathrm{box}_{2f_0} (f)e^{\hbox {j}2\pi ft} \mathrm{d}f=h\int_{-f_0}^{+f_0} e^{\hbox {j}2\pi ft}\, \mathrm{d}f\\ &=\left (\frac{h}{2\pi jt}\right) e^{\hbox {j}2\pi ft} |_{+f_0}^{-f_0}=\frac{h}{2\pi jt} (e^{\hbox {j}2\pi f_0t} -e^{-{\hbox {j}}2\pi f_0t}).\\ &=\frac{h}{\pi t} \sin(2\pi f_{0} t)=\frac{2f_{0} h}{(2\pi f_{0} t)} \sin\,(2\pi f_{0} t)=A \sin c (2\pi f_{0} t), \end{align} (9.3k)

where ${\displaystyle A=2hf_{0}}$ and sinc $\displaystyle x=(\sin x)/x.$

9.3d Calculate the transform of the pair of displaced boxcars in Figure 9.3b. Discuss the relation between your result and equation (9.3k) for a single boxcar centered at the origin.

Solution

Equation (9.3d) gives for the transform of the pair of boxcars:

 {\displaystyle {\begin{aligned}f(t)&={\frac {h}{2\pi }}\left[\int \limits _{-\omega _{2}}^{-\omega _{1}}e^{{\hbox{j}}\omega t}d\omega +\int \limits _{+\omega _{1}}^{+\omega _{2}}e^{{\hbox{j}}\omega t}d\omega \right]={\frac {h}{2\pi jt}}\left[e^{{\hbox{j}}\omega t}|_{-\omega _{2}}^{-\omega _{1}}+e^{{\hbox{j}}\omega t}|_{+\omega _{1}}^{+\omega _{2}}\right]\\&={\frac {h}{2\pi jt}}[(e^{-{\hbox{j}}\omega _{1}t}-e^{-{\hbox{j}}\omega _{2}t})+(e^{{\hbox{j}}\omega _{2}t}-e^{-{\hbox{j}}\omega _{1}t})]={\frac {h}{\pi t}}(\sin \omega _{2}t-\sin \omega _{1}t)\\&=(h/\pi )(\omega _{2}\sin \mathrm {c} \ \omega _{2}t-\omega _{1}\sin \mathrm {c} \ \omega _{1}t)\\&=A_{2}\sin \mathrm {c} \ \omega _{2}t-A_{1}\sin \mathrm {c} \ \omega _{1}t,\end{aligned}}} (9.3l)

where ${\displaystyle A_{2}=2hf_{2}}$, $\displaystyle A_{1} =2hf_{1}$ , and we have used Euler’s formulas (see Sheriff and Geldart, 1995, problem 15.12a) to get the sines. Thus, the transform is the difference between two sinc functions corresponding to the upper and lower limiting frequencies of the boxcars.

Note that equation (9.3$\displaystyle \ell$ ) can be regarded as giving the result of two boxcars, one extending from ${\displaystyle -\omega _{2}}$ to $\displaystyle +\omega _{2}$ , the other extending from ${\displaystyle -\omega _{1}}$ to $\displaystyle +\omega _{1}$ , the effect of the latter being subtracted from that of the former.

To compare equation (9.3l) with equation (9.3k), we set $\displaystyle \omega _{1} =0$ . Then sinc $\displaystyle \omega _{1} t= \sin\,0/0=1$ , so ${\displaystyle \omega _{1}}$ sinc $\displaystyle \omega _{1} t=0$ and ${\displaystyle f(t)}$ becomes

\displaystyle \begin{align} f(t)=2hf_{2} \sin \mathrm{c}\ \omega _{2} t. \end{align}

Setting $\displaystyle f_{2} =f_{0}$ gives equation (9.3k).

## 9.4 Extension of the sampling theorem

Explain why sampling at 4 ms is sufficient to reproduce exactly a signal whose spectrum is within the range $\displaystyle 0$ to $\displaystyle 125\ {\hbox {Hz}} =f_{N}$ , but not if the bandwidth is shifted upward?

Background

To convert an analog (continuous) signal to the digital form, values of the signal are measured at a fixed interval ${\displaystyle \Delta }$, called the sampling interval. The result is a series of numbers representing the amplitudes and polarities of the signal at the times $\displaystyle t=n\Delta$ , where $\displaystyle n$ is an integer (often we omit $\displaystyle \Delta$ and specify the time by giving ${\displaystyle n}$ only). We write ${\displaystyle g_{t}}$ for the digital function corresponding to ${\displaystyle g(t)}$.

Sampling is equivalent to multiplying the analog signal by a comb, a function consisting of an infinite series of unit impulses spaced at a fixed interval ${\displaystyle \Delta }$. The equation of comb(t) is

 {\displaystyle {\begin{aligned}\mathrm {comb} (t)=\sum \limits _{r=-\infty }^{+\infty }\delta (t+r\Delta ),\\{\mbox{so}}\qquad \qquad g_{t}=g(t)\sum \limits _{r=-\infty }^{+\infty }\delta (t+r\Delta ).\end{aligned}}} (9.4a)

The transform of comb(t) is another comb with impulses in the frequency domain at intervals ${\displaystyle 1/\Delta }$:

 {\displaystyle {\begin{aligned}\mathrm {comb} (t)\leftrightarrow k\,\mathrm {comb} (f),\end{aligned}}} (9.4b)

where $\displaystyle k=2\pi /\Delta$ [see Sheriff and Geldart, 1995, equation (15.155)].

Figure 9.4a.  Sampling and recovering a waveform. The left column is time domain, the right is frequency domain.

The Nyquist frequency, $\displaystyle f_{N}$ , is half the frequency of sampling $\displaystyle 1/\Delta$ , that is,

 \displaystyle \begin{align} f_{N} =1/2\Delta. \end{align} (9.4c)

While a signal that does not contain frequencies higher than $\displaystyle f_{N}$ will be recorded accurately, a frequency higher than $\displaystyle f_{N}$ by the amount ${\displaystyle \Delta f}$ (i.e., a frequency that is not sampled at least twice per cycle) will produce an alias frequency equal to ${\displaystyle (f_{N}-\Delta f)}$ (Sheriff and Geldart, 1995, section 9.2.2c).

The convolution theorem, equation (9.3f), has a converse [see Sheriff and Geldart, 1995, equation (15.146)]:

 \displaystyle \begin{align} g(t)h(t)\leftrightarrow \frac{1}{2\pi} G(\omega)*H(\omega). \end{align} (9.4d)

Solution

We examine first the case where the signal frequency spectrum extends from ${\displaystyle 0}$ to 125 Hz, next consider modifications when the spectrum is shifted upward 125 Hz $\displaystyle (=f_{N}$ for 4-ms sampling), and finally, the case where the signal is shifted upward by an amount ${\displaystyle f\neq 125n}$.

Part (i) of Figure 9.4a shows the signal and its frequency spectrum, the latter being confined to the range ${\displaystyle -125\ {\hbox{to}}+125}$ Hz [the negative frequencies arise when we use Euler’s equations (Sheriff and Geldart, 1995, problem 15.12a) to transform equation (9.1a) into (9.1c). Part (ii) shows comb(t), with elements 4 ms apart, and its transform, a comb with elements spaced ${\displaystyle 1/\Delta =250}$ Hz apart. Digitizing the signal is equivalent to multiplying it by comb(t) [part (iii)], which is equivalent to convolving the spectrum ${\displaystyle G(f)}$ with comb(f) [see equation (9.4b)]; this results in a repetition of the spectrum for each impulse in comb(f) [see equation (9.3j)]. To eliminate the repeated spectra, we multiply by a boxcar [part (iv)], thus getting back the original spectrum [part (v)]. Multiplication in the frequency domain is equivalent to convolution in the time domain, so we convolve with a sinc function in part (V) and recover the original time-domain signal.

Figure 9.4b.  Repeat of Figure 9.4a for a 125–250 Hz signal.
Figure 9.4c.  Repeat of Figure 9.4a for spectrum overlapping the Nyquist frequency.

If the signal is shifted upward by ${\displaystyle 125n}$, where $\displaystyle n$ is an integer, then the sampling interval is larger than 1/2 cycle for all frequencies and the resulting spectrum is that of the alias function rather than of the signal, as illustrated in Figure 9.4b. Because we no longer have the signal spectrum after sampling, we cannot recover the signal.

If the signal is shifted upward so that it overlaps the Nyquist frequency (Figure 9.4c), the sampling causes the frequencies above ${\displaystyle f_{\rm {N}}}$ to overlap and distort the signal spectrum. Now we cannot separate the signal spectrum from the distorted spectrum and hence cannot recover the signal.

## 9.5 Alias filters

9.5a The standard alias filter such as that shown in Figure 9.5a has a 3-dB point at about half the Nyquist frequency ${\displaystyle f_{N}}$ and a very steep slope, so that noise above ${\displaystyle f_{N}}$ is highly attenuated relative to the passband of the system. Assuming an original flat spectrum, alias filtering with a 125-Hz, 72-dB/octave filter, and subsequent resampling from 2 to 4 ms (without additional alias filtering), graph the resulting alias noise versus frequency.

Background

The passband of a system is the range of frequencies that are unattenuated by passage through the system. The limits are usually taken as the frequencies for which the attenuation is 3 dB.

Figure 9.5a.  Seismic filter responses.

Solution

With the 3-dB point on the high-cut filter at 125 Hz and a 72 dB/octave slope, the attenuation at ${\displaystyle f_{N}=250}$ Hz is 75 dB. Frequencies ${\displaystyle \Delta {f}}$ higher than the Nyquist frequency, $\displaystyle f_{N} +\Delta f$ , alias to appear as the frequencies ${\displaystyle f_{N}-\Delta {f}}$, i.e., they fold back about ${\displaystyle f_{N}}$, as shown in Figure 9.5a.

A 65-Hz 72-dB/octave filter is usually applied before resampling to prevent aliasing. If resampling to 4 ms is done without this additional filtering, ${\displaystyle f_{N}}$ for the resampled data is $\displaystyle 1/(2\times 0.004)=125$ Hz; the aliased noise is shown by the foldback curves in Figure 9.5a.

9.5b Some believe that standard alias filters may be unnecessarily restrictive. The standard alias filter for 4-ms sampling is about 65 Hz, 72-dB/octave. Graph the alias noise versus frequency for a 90-Hz, 72-dB/octave filter for 4-ms sampling and draw conclusions.

Solution

The 90-Hz filter is shown in Figure 9.5a. Because the sampling interval is 4 ms, $\displaystyle f_{N}$ is 125 Hz and the alias noise is given by the foldback curves shown. The 90-Hz filter gives a broader passband than the standard 65-Hz alias filter if signal frequencies above about 80 Hz and signals attenuated more than 60 dB are not important.

## 9.6 The convolutional model

9.6a The earth acts as a filter to seismic energy to change the waveshape so that what we record is a distorted version of the waveshape we send into the earth. Assume that the signal generated by a seismic source is so short that it can be represented adequately by an impulse; what aspects of passage through the earth will distort the impulse and what will be the effects?

Solution

The signal from a seismic source will be distorted: (a) by the region near the source where the stresses are so high as to produce nonlinear effects $\displaystyle s_{i}$ , (b) by the near-surface region where absorption is especially large ${\displaystyle n_{i}}$, (c) by absorption in passage through the main body of the earth $\displaystyle p_{i}$ , and (d) perhaps by other factors. The first of these will change an impulsive signal to a waveform ${\displaystyle s_{i}}$, which we can think of as a series of impulses with amplitudes that give the shape ${\displaystyle s_{i}}$. Passage through the near surface will replace each of these impulses by ${\displaystyle n_{i}}$ multiplied by the scaled value ${\displaystyle s_{i}}$, that is, ${\displaystyle s_{i}*n_{i}}$. Absorption will change this series to $\displaystyle s_{i} *n_{i} *p_{i}$ (where the absorption ${\displaystyle p_{i}}$ is frequency dependent). Reflection from interfaces in the earth can also be thought of as additional convolutions changing the waveform to ${\displaystyle s_{i}*n_{i}*p_{i}*r_{i}}$, etc. Thus the effect on the waveform can be thought of as a sequence of convolutions. In ordinary reflection seismic work the objective is to determine the reflectivity ${\displaystyle r_{i}}$ so we combine all the other waveshape-changing factors into an equivalent wavelet, generally called the embedded wavelet ${\displaystyle w_{i}}$:

{\displaystyle {\begin{aligned}s_{i}*n_{i}*p_{i}*r_{i}=w_{i}*r_{i},\end{aligned}}}

where ${\displaystyle w_{i}=s_{i}*n_{i}*p_{i}}$. This concept is called the convolutional model.

9.6b What is the effect if the source is not impulsive?

Solution

A nonimpulsive source ${\displaystyle v_{i}}$, such as the vibroseis, has the effect of replacing ${\displaystyle s_{i}}$ with ${\displaystyle v_{i}}$, but otherwise the reasoning in part (a) is unchanged.

## 9.7 Water reverberation filter

9.7a Verify equation (9.7d) for the inverse filter for water reverberation by convolving it with equation (9.7a) for the water reverberation.

Background

Multiple reflections within a water layer can cause severe reverberation (ringing) when the reflection coefficient at the sea floor is sufficiently high. We take the reflection coefficients at the top and bottom of the water layer as ${\displaystyle -1}$ and $\displaystyle R$ , and ${\displaystyle n\Delta }$ as the round-trip traveltime in the layer, ${\displaystyle n}$ being an integer and $\displaystyle \Delta$ the sampling interval. Bounces can occur before or after a wave is reflected at deeper reflectors. We assign unit amplitude to a primary reflection with traveltime ${\displaystyle t}$. A reflection from the same horizon that has suffered one bounce in the water layer has traveltime $\displaystyle t+n\Delta$ and amplitude $\displaystyle -R$ . Since the bounce may take place either before or after reflection at depth, there are two waves arriving at ${\displaystyle t+n\Delta }$, the combined amplitude being ${\displaystyle -2R}$. There are three waves that suffer two bounces: two bounces before the deep reflection, two after, and one before plus one after, so the multiple that arrives at $\displaystyle t+2n\Delta$ has amplitude $\displaystyle +3R^{2}$ . Reflections suffering three bounces arrive at $\displaystyle t+3n\Delta$ with amplitude ${\displaystyle -4R^{3}}$, and so on. The impulse response (see Sheriff and Geldart, 1955, 279) of the water layer is thus the infinite series

 \displaystyle \begin{align} w_{t} =[1,\; -2R,\; 3R^{2} ,\; -4R^{3} ,\; 5R^{4} ,\ldots\ldots.]\quad \mathrm{when}\ n=1; \end{align} (9.7a)

 {\displaystyle {\begin{aligned}{\mbox{or}}\qquad \qquad w_{t}=[1,\;0,\;-2R,\;0,\;3R^{2},\;0,\;-4R^{3},\ldots .]\qquad \mathrm {when} \ n=2;\end{aligned}}} (9.7b)

and so on for other values of ${\displaystyle n}$.

A filter (see problem 7.11) that removes the effect of another filter is an inverse filter, and the inverse filter in this case is $\displaystyle i_{t}$ ; it satisfies the equation

 \displaystyle \begin{align} w_{t} *i_{t} =\delta_{t}; \end{align} (9.7c)

 \displaystyle \begin{align} i_{t} &=[1,\; 2R,\; R^{2} ,\; 0,\; 0,\; .\; .\; .]\qquad \mathrm{when}\ n=1 \end{align} (9.7d)

 {\displaystyle {\begin{aligned}&=[1,\;0,\;2R,\;0,\;R^{2},\;0,\;.\;.\;.]\quad \mathrm {when} \ n=2.\end{aligned}}} (9.7e)

[See Sheriff and Geldart, 1995, section 15.5.3 for a discussion of $\displaystyle z$ -transforms.]

Solution

For ${\displaystyle n=1}$, $\displaystyle i_{t} =\left[1,\; 2R,\; R^{2} \right]$ , and ${\displaystyle w_{t}=[1,-2R,3R^{2},-4R^{3},5R^{4},\ldots ]}$. We wish to verify that ${\displaystyle w_{t}*i_{t}=\delta _{t}}$. Using the procedure described in problem 9.2, we replace each element of $\displaystyle i_{t}$ by a scaled version of ${\displaystyle w_{t}}$; thus

 1, –$\displaystyle 2R$ , +$\displaystyle 3R^{2}$ , –${\displaystyle 4R^{3}}$, +$\displaystyle 5R^{4},\ldots$ +${\displaystyle 2R}$, –$\displaystyle 4R^{2}$ , +${\displaystyle 6R^{3}}$, –$\displaystyle 8R^{4}, \ldots$ +${\displaystyle R^{2}}$, –$\displaystyle 2R^{3}$ , +${\displaystyle 3R^{4},\ldots }$ = 1, 0, 0, 0, 0, … ${\displaystyle =\delta _{t}}$.

We can proceed in the same way for other values of $\displaystyle n$ .

9.7b The spectrum of the water-layer filter is shown in Figure 9.7a for $\displaystyle n=1$ ; the large peaks occur at the “singing frequency” $\displaystyle nf_{1}$ , where ${\displaystyle n}$ is any odd number. Sketch the amplitude spectrum of the inverse filter.

Solution

In the time domain, $\displaystyle w_{i} *i_{t} =\delta_{t}$ ; in the frequency domain this becomes $\displaystyle W(f)\times I(f)=1$ . Thus the amplitude spectrum of ${\displaystyle I(f)}$ is the inverse of that of ${\displaystyle W(f)}$, that is, for a given frequency,

{\displaystyle {\begin{aligned}|I(f)|=1/|W(f)|,\end{aligned}}}

as shown in Figure 9.7b.

Figure 9.7a.  Spectrum of water-layer filter; ${\displaystyle f_{1}=V_{W}/4h,h=}$ water depth, $\displaystyle R = 0.5$ .
Figure 9.7b.  Inverse filter spectrum (dashed).

9.7c Verify your sketch of the water-reverberation inverse filter by transforming

 \displaystyle \begin{align} \left[1, \; 2R,\; R^{2}\right]\leftrightarrow [1+2Re^{-{\hbox {j}}2\pi f{\Delta}} + R^{2} e^{-{\hbox {j}}4\pi f{\Delta}}, \end{align} (9.7f)

and calculating the value of the spectrum for ${\displaystyle R=0.5}$, ${\displaystyle f=0}$, ${\displaystyle 1/2\Delta }$, and ${\displaystyle 1/\Delta }$.

Solution

We start with the ${\displaystyle z}$-transform of $\displaystyle i_{t}$ , that is, ${\displaystyle I(f)=(1+2Rz+R^{2}z^{2})}$. This is a complex quantity and we multiply by the conjugate complex (see Sheriff and Geldart, 1995, section 15.1.5) to get the magnitude squared. Since $\displaystyle z=e^{\hbox {j}2\pi\Delta}$ , the conjugate complex is ${\displaystyle e^{+{\hbox{j}}2\pi \Delta }}$ and $\displaystyle zz^{-1} =1$ ; thus,

\displaystyle \begin{align} |I(f)|^{2} &=(1+2Rz+R^{2} z^{2})(1+2Rz^{-1} +R^{2} z^{-2})\\ &=(1+Rz)^{2} (1+Rz^{-1})^{2} =[(1+Rz)(1+Rz^{-1})]^{2},\end{align} \displaystyle \begin{align} \mbox {so} \qquad \qquad |I(f)|=(1+Rz)(1+Rz^{-1})=[1+R^{2} +R(e^{+\hbox {j}2\pi f\Delta} +e^{-{\hbox {j}}2\pi f\Delta})]. \end{align}

Substituting ${\displaystyle R=0.5}$ and using Euler’s formula (Sheriff and Geldart, 1995, problem 15.12a), we get

\displaystyle \begin{align} |I(f)|=[1.25+\cos\,(2\pi f\Delta)]. \end{align}

Substituting ${\displaystyle f=0}$, ${\displaystyle 1/2\Delta }$, and ${\displaystyle 1/\Delta }$, we get the following values

{\displaystyle {\begin{aligned}|I(f)|=2.25,\ 0.25,\ 2.25,\ \quad \mathrm {that\ is} ,9/4,\ 1/4,\ 9/4.\end{aligned}}}

The peaks and troughs in Figure 9.7a have the relative values 9, 1, 9; since these values are normalized, the agreement is exact.

## 9.8 Calculating crosscorrelation and autocorrelation

Four causal wavelets are given by ${\displaystyle a_{t}=[2,\;-1]}$, ${\displaystyle b_{t}=[4,1]}$, $\displaystyle c_{t} =[6,\; -7,\; 2]$ , ${\displaystyle d_{t}=[4,9,2]}$. Calculate the crosscorrelations and autocorrelations ${\displaystyle \phi _{ab}}$, $\displaystyle \phi _{ac}$ , ${\displaystyle \phi _{ca}}$, $\displaystyle \phi _{aa}$ , and $\displaystyle \phi _{cc}$ in both the time and frequency domains.

Background

A causal wavelet as defined in problem 5.21 has zero values when ${\displaystyle t<0}$.

The crosscorrelation ${\displaystyle \phi _{gh}(\tau )}$ of ${\displaystyle g_{t}}$ and ${\displaystyle h_{t}}$ tells us how similar the two functions are when $\displaystyle h_{t}$ is shifted the amount ${\displaystyle \tau }$ relative to $\displaystyle g_{t}$ . The crosscorrelation is given by the equation

 {\displaystyle {\begin{aligned}\phi _{gh}(\tau )=\sum \limits _{k}^{}g_{k}h_{k+\tau }.\end{aligned}}} (9.8a)

This equation means that $\displaystyle h_{t}$ is displaced ${\displaystyle \tau }$ units to the left relative to ${\displaystyle g_{t}}$, corresponding values multiplied, and the products summed to give $\displaystyle \phi_{gh} (\tau)$ . The result is the same if we move $\displaystyle g_{t}$ to the right ${\displaystyle \tau }$ units, that is,

 \displaystyle \begin{align} \phi _{gh} (\tau)=\phi_{hg} (-\tau). \end{align} (9.8b)

The functions $\displaystyle \phi_{gh} (\tau)$ and $\displaystyle g_{t} *h_{t}$ are closely related. Using two simple curves, it is easily shown that we can crosscorrelate by reversing $\displaystyle g_{t}$ and convolving the reversed function with ${\displaystyle h_{t}}$, that is,

 {\displaystyle {\begin{aligned}\phi _{gh}(\tau )=g_{-\tau }*h_{\tau }=g_{\tau }*h_{-\tau }\end{aligned}}} (9.8c)

(since convolution is commutative). Reversing a function in time changes the sign of $\displaystyle t=n\Delta$ , so the exponent of $\displaystyle z$ changes sign also, ${\displaystyle z}$ becoming ${\displaystyle z^{-1}}$ and $\displaystyle G(z)$ becoming the conjugate complex ${\displaystyle {\bar {G}}(z)}$. The convolution theorem (equation (9.3f)) now becomes the crosscorrelation theorem:

 {\displaystyle {\begin{aligned}\phi _{gh}(\tau )\leftrightarrow {\bar {G}}(z)H(z).\end{aligned}}} (9.8d)

If ${\displaystyle h_{t}}$ is the same as ${\displaystyle g_{t}}$, we get the autocorrelation of ${\displaystyle g_{t}}$ and equations (9.8a,d) become

 {\displaystyle {\begin{aligned}\phi _{gg}(\tau )=\sum \limits _{k}^{}g_{k}\ g_{k+\tau }\leftrightarrow |G(z)|^{2},\end{aligned}}} (9.8e)

the transform relation being the autocorrelation theorem. Since the two functions are the same, $\displaystyle \phi_{gg} (\tau)$ does not depend upon the direction of displacement. The autocorrelation for ${\displaystyle \tau =0}$ equals the sum of the data elements squared, hence is called the energy of the trace;

 {\displaystyle {\begin{aligned}\phi _{gg}(0)=\sum \limits _{k}^{}g_{k}^{2}.\end{aligned}}} (9.8f)

Both the autocorrelation and the crosscorrelation are often normalized; in the case of the autocorrelation, equation (9.8e) becomes

 {\displaystyle {\begin{aligned}\phi _{gg}(\tau )_{\mathrm {norm} }=\left(\sum \limits _{k}^{}g_{k}\ g_{k+\tau }\right)\left(\sum \limits _{k}^{}g_{k}^{2}\right)^{-1}=\phi _{gg}(\tau )/\phi _{gg}(0).\end{aligned}}} (9.8g)

The normalized crosscorrelation is

 \displaystyle \begin{align} \phi_{gh} (\tau)_{\mathrm{norm}} =\phi_{gh} (\tau)/[\phi_{gg} (0)\phi_{hh} (0)]^{1/2}. \end{align} (9.8h)

Solution

Time-domain calculations:

Using equation (9.8a) to calculate ${\displaystyle \phi _{ab}}$, i.e., ${\displaystyle b_{t}}$ is first shifted to the left; we have:

{\displaystyle {\begin{aligned}\phi _{ab}(+1)=2\times 1=2;\quad \phi _{ab}(0)=2\times 4-1\times 1=7;\quad \phi _{ab}(-1)=-1\times 4=-4.\end{aligned}}}

So

{\displaystyle {\begin{aligned}\phi _{ab}(\tau )=[-4,\mathop {7} \limits ^{\downarrow },\;2],\end{aligned}}}

where ${\displaystyle \downarrow }$ marks the value at $\displaystyle \tau=0$ . Proceeding in the same way, we find that

{\displaystyle {\begin{aligned}&\phi _{ac}(+2)=4,\quad \phi _{ac}(+1)=-16,\quad \phi _{ac}(0)=19,\quad \phi _{ac}(-1)=-6;\\&\phi _{ac}(\tau )=[-6,\mathop {19} \limits ^{\downarrow },-16,4].\end{aligned}}}

To find ${\displaystyle \phi _{ca}(\tau )}$ we displace ${\displaystyle a_{t}}$ instead of ${\displaystyle c_{t}}$ and obtain ${\displaystyle \phi _{ca}=[4,\;-16,\;\mathop {19} \limits ^{\downarrow },\;-6]}$, which equals ${\displaystyle \phi _{ac}(-\tau )}$.

Autocorrelations are found in the same way:

{\displaystyle {\begin{aligned}\phi _{aa}(\tau )=[-2,\;\mathop {5} \limits ^{\downarrow },\;-2];\quad \phi _{cc}=[12,\;-56,\;\mathop {89} \limits ^{\downarrow },\;-56,\;12].\end{aligned}}}

Frequency-domain calculations

The ${\displaystyle z}$-transforms of ${\displaystyle a_{t}}$, ${\displaystyle b_{t}}$, and ${\displaystyle c_{t}}$ are ${\displaystyle A(z)=(2-z)}$, ${\displaystyle B(z)=(4+z)}$, $\displaystyle C(z)=(6-7z+2z^{2})$ . The conjugate complexes of these transforms are are ${\displaystyle A(z)=(2-z^{-1})}$, ${\displaystyle B(z)=(4+z^{-1})}$, ${\displaystyle C(z)=(6-7z^{-1}+2z^{-2})}$. Using these transforms, we get

{\displaystyle {\begin{aligned}&\phi _{ab}(\tau )\leftrightarrow {\bar {A}}(z)B(z)=(2-z^{-1}(4+z)\\&=-4z^{-1}+7+2z\leftrightarrow [-4,\;\mathop {7} \limits ^{\downarrow },\;2];\\&\phi _{ac}(\tau )\leftrightarrow {\bar {A}}(z)C(z)=(2-z^{-1})(6-7z+z^{2})\\&=-6z^{-1}+19-16z+4z^{2}\leftrightarrow [-6,\;\mathop {19} \limits ^{\downarrow },\;-16+4];\\&\phi _{ca}(\tau )\leftrightarrow {\bar {C}}(z)A(z)=(2z^{-2}-7z^{-1}+6)(2-z)\\&=4z^{-2}-16z^{-1}+19-6z\leftrightarrow [4,\;-16,\;\mathop {19} \limits ^{\downarrow },\;-6];\\&\phi _{aa}(\tau )\leftrightarrow (2-z^{-1})(2-z)=-2z^{-1}+5-2z\leftrightarrow [-2,\;\mathop {5} \limits ^{\downarrow },\;-2];\\&\phi _{cc}(\tau )\leftrightarrow (6-7z^{-1}+2z^{-2})(6-7z+2z^{2})\\&=12z^{-2}-56z^{-1}+89-56z+12z^{2}\leftrightarrow [12,\;-56,\;\mathop {89} \limits ^{\downarrow },\;-56,\;12].\end{aligned}}}

## 9.9 Digital calculations

Fill in the values in Table 9.9a.

Solution

Table 9.9b shows Table 9.9a completed.

Table 9.9a. Digital wavelets and operations.
${\displaystyle t=-3}$ $\displaystyle t=-2$ ${\displaystyle t=-1}$ ${\displaystyle t=0}$ ${\displaystyle t=+1}$ ${\displaystyle t=+2}$ ${\displaystyle t=+3}$ $\displaystyle t=+4$ $\displaystyle t=+5$
${\displaystyle a_{t}=\left[2,\;1\;,\;-2,\;1\right]}$
${\displaystyle b_{t}=-2a_{t}}$
${\displaystyle c_{t}{=3a_{t-2}}}$
${\displaystyle {\hbox{d}}_{t}{=a_{-t}}/2}$
${\displaystyle e_{t}=\pi a_{3-t}}$
${\displaystyle f_{t}=\left[-1,\;1\right]}$
${\displaystyle g_{t}{=a_{t}*f_{t}}}$
${\displaystyle \delta _{t+2}}$
${\displaystyle \delta _{2-t}}$
${\displaystyle \phi _{ff}(t)}$
${\displaystyle \phi _{fa}(t)}$
Table 9.9b. Completed table.
${\displaystyle t=-3}$ ${\displaystyle t=-2}$ $\displaystyle t=-1$ ${\displaystyle t=0}$ ${\displaystyle t=+1}$ $\displaystyle t=+2$ ${\displaystyle t=+3}$ $\displaystyle t=+4$ ${\displaystyle t=+5}$
${\displaystyle a_{t}=\left[2,\;1\;,\;-2,\;1\right]}$ 0 0 0 2 1 ${\displaystyle -2}$ 1 0 0
${\displaystyle b_{t}=-2a_{t}}$ 0 0 0 $\displaystyle -4$ ${\displaystyle -2}$ 4 ${\displaystyle -2}$ 0 0
$\displaystyle c_t{=3a_{t-2}}$ 0 0 0 0 0 6 3 ${\displaystyle -6}$ 3
${\displaystyle {\hbox{d}}_{t}{=a_{-t}}/2}$ 1/2 ${\displaystyle -1}$ 1/2 1 0 0 0 0 0
${\displaystyle e_{t}=\pi a_{3-t}}$ 0 0 0 ${\displaystyle \pi }$ ${\displaystyle -2\pi }$ ${\displaystyle \pi }$ $\displaystyle 2\pi$ 0 0
$\displaystyle f_{t} =\left[-1,\; 1\right]$ 0 0 0 ${\displaystyle -1}$ 1 0 0 0 0
${\displaystyle g_{t}{=a_{t}*f_{t}}}$ 0 0 0 ${\displaystyle -2}$ 1 3 ${\displaystyle -3}$ 1 0
${\displaystyle \delta _{t+2}}$ 0 1 0 0 0 0 0 0 0
${\displaystyle \delta _{2-t}}$ 0 0 0 0 0 1 0 0 0
$\displaystyle \phi_{ff}(t)$ 0 0 ${\displaystyle -1}$ 2 ${\displaystyle -1}$ 0 0 0 0
${\displaystyle \phi _{fa}(t)}$ 0 0 0 ${\displaystyle -1}$ $\displaystyle -1$ ${\displaystyle -3}$ 3 ${\displaystyle -1}$ 0

## 9.10 Semblance

Semblance ${\displaystyle S_{T}}$ is the ratio of the energy of a stack of ${\displaystyle N}$ traces to ${\displaystyle N}$ times the sum of the energies of the ${\displaystyle N}$ component traces, all summed over some time interval. Show that ${\displaystyle S_{T}=1}$ when the traces are identical.

Background

While crosscorrelation is a quantitative measure of the similarity of two traces, semblance ${\displaystyle S_{T}}$ is a measure of the similarity of a number of traces. Assuming ${\displaystyle N}$ traces, we sum (stack) the traces at time ${\displaystyle t}$, square the sum to get the total energy, and sum the values over a time interval ${\displaystyle {m}\Delta }$. The equation for ${\displaystyle S_{T}}$ is

 {\displaystyle {\begin{aligned}S_{T}=\sum \limits _{t=t_{0}}^{t_{0}+m\Delta }\left(\sum \limits _{i=1}^{N}g_{ti}\right)^{2}/N\left[\sum \limits _{t=t_{0}}^{t_{0}+m\Delta }\sum \limits _{i=1}^{N}(g_{ti})^{2}\right].\end{aligned}}} (9.10a)

Solution

The sum of ${\displaystyle N}$ traces is ${\displaystyle (\sum {_{i=1}^{N}g_{ti}})}$ and its energy is ${\displaystyle (\sum {_{i=1}^{N}}g_{ti})^{2}}$. If the traces are identical, the numerator of equation (9.10a) becomes ${\displaystyle N^{2}[\sum {^{t_{0}+m\Delta }_{t=t_{0}}}g_{t}^{2}]}$. The denominator becomes ${\displaystyle N\sum {^{t_{0}+m\Delta }_{t=t_{0}}}(\sum {^{N}_{i=1}}g_{t}^{2})=N^{2}\sum {}_{t=t_{0}}^{t_{0}+m\Delta }g_{t}^{2}}$. Since numerator and denominator are equal, ${\displaystyle S_{T}=1}$.

## 9.11 Convolution and correlation calculations

9.11a Convolve ${\displaystyle g_{t}=[2,\;5,\;-2,\;1]}$ with ${\displaystyle f_{t}=[6,\;-1,\;-1]}$.

Background

Convolution is discussed in problem 9.2, crosscorrelation and autocorrelation in problem 9.8, minimum-phase wavelets in Sheriff and Geldart, 1995, section 9.4 and section 15.5.6a, ${\displaystyle z}$-transforms in Sheriff and Geldart, 1995, section 15.5.3. One of the definitions of a minimum-phase function is that, when the function is expressed as the product of its factors in the form ${\displaystyle \prod {_{n}}{(a_{n}-z)}}$, the magnitudes of the roots ${\displaystyle |a_{n}|}$ must all be greater than unity.

Solution

To evaluate the convolution, we replace each element of ${\displaystyle f_{t}=[6,\;-1,\;-1]}$ with a scaled version of ${\displaystyle g_{t}=[2,\;5,\;-2,\;1]}$ and then sum values for the same times.

 6${\displaystyle g_{t}}$: 12, 30, –12, 6 –${\displaystyle g_{t-1}}$: –2, –5, 2, –1 –${\displaystyle g_{t-2}}$: –2, –5, 2, –1 Sum: 12, 28, –19, 3, 1, –1

The same result can be obtained by reversing one function and sliding it past the other, multiplying values at the same times and summing the results:

 ${\displaystyle g_{t}\to }$ 2 5 –2 1 Sum t ${\displaystyle f_{-t}\to }$ –1 –1 6 0 –1 –1 6 +12 1 –1 –1 6 +28 2 –1 –1 6 –19 3 –1 –1 6 +3 4 –1 –1 +1 5 –1 –1

Still another method is to use ${\displaystyle z}$-transforms; thus,

{\displaystyle {\begin{aligned}\left[2,\;5,\;-2,\;1\right]*[6,\;-1,\;-1]\leftrightarrow (2+5z-2z^{2}+z^{3})(6-z-z^{2})\\\leftrightarrow (12+28z-19z^{2}+3z^{3}+z^{4}-z^{5}).\end{aligned}}}

Transforming back to the time domain, we have

{\displaystyle {\begin{aligned}f_{t}*g_{t}=[12,\;28,\;-19,\;+3,\;+1,\;-1].\end{aligned}}}

9.11b Crosscorrelate ${\displaystyle g_{t}=[2,\;5,\;-2,\;1]}$ with ${\displaystyle f_{t}=[6,\;-1,\;-1]}$. For what shift are these functions most nearly alike?

Solution

Using equation (9.8a) for ${\displaystyle \phi _{g_{f}}(\tau )}$, we write

 ${\displaystyle \tau }$ ${\displaystyle g_{t}\;\to }$ 2 5 -2 1 ${\displaystyle \phi _{gf}(\tau )}$ +2 ${\displaystyle f_{t+\tau }\;\to }$ 6 -1 -1 -2 +1 6 -1 -1 -7 0 6 -1 -1 +9 -1 6 -1 -1 +31 -2 6 -1 -13 -3 6 +6

The functions are most alike for a shift of ${\displaystyle -1}$.

9.11c Convolve ${\displaystyle g_{t}=[2,\;5,\;-2,\;1]}$ with ${\displaystyle f_{-t}=[-1,\;-1,\;6]}$. Compare with the answer in part (b) and explain.

Solution

Using the convolution method of replacing the first function by a scaled version of itself, the scale factors being the elements of the second function, we obtain

 –2 5 –2 1 –2 –5 2 –1 12 30 –12 6 –2, –7, 9, 31, –13, 6

Thus, the convolution is ${\displaystyle [-2,\;-7,\;9,\;31,\;-13,\;6]}$, the same result achieved by the correlation in part (b). Thus, ${\displaystyle g_{t}*f_{-t}=\phi _{fg}(\tau )}$.

9.11d Autocorrelate ${\displaystyle f_{t}=[6,\;-1,\;-1]}$ and ${\displaystyle h_{t}=[3,\;-5,\;-2]}$. The autocorrelation of a function is not unique to that function (see Sheriff and Geldart, 1995, 552). Other wavelets having the same autocorrelation as the preceding are ${\displaystyle [-2,\;-5,\;3]}$ and ${\displaystyle [-1,\;-1,\;6]}$. Which of the four is the minimum-phase wavelet?

Solution

We use ${\displaystyle z}$-transforms to get the autocorrelations:

{\displaystyle {\begin{aligned}f_{t}=[6,\;-1,\;-1]\leftrightarrow (6-z-z^{2});\\\phi _{ff}(\tau )\leftrightarrow (6-z-z^{2})(6-z^{-1}-z^{-2})\\\leftrightarrow -6z^{-2}-5z^{-1}+38-5z-6z^{2};\\\phi _{ff}(\tau )=-6,\;-5,\;\mathop {38} \limits ^{\downarrow },\;-5,\;-6;\\h_{t}=[3,\;-5,\;-2]\leftrightarrow (3-5z-2z^{2});\\\phi _{hh}(\tau )\leftrightarrow (3-5z-2z^{2})(3-5z^{-1}-2z^{-2})\\\leftrightarrow -6z^{-2}-5z^{-1}+38-5z-6z^{2};\\\phi _{hh}(\tau )=[-6,\;-5,\mathop {38} \limits ^{\downarrow },\;-5,\;-6].\end{aligned}}}

To determine which of the four wavelets are minimum-phase, we factor each of the transforms:

 Roots ${\displaystyle (6-z-z^{2})=(3+z)(2-z);}$ ${\displaystyle -3,2}$ ${\displaystyle (3-5z-2z^{2})=(3+z)(1-2z);}$ ${\displaystyle -3,1/2}$ ${\displaystyle (-2-5z+3z^{2})=(1+3z)(-2+z);}$ ${\displaystyle -1/3,2}$ ${\displaystyle (-1-z+6z^{2})=(1+3z)(-1+2z).}$ ${\displaystyle -1/3,1/2}$

Since only the first wavelet has all roots of magnitude greater than 1, it is the only one that is minimum-phase.

9.11e What is the normalized autocorrelation of ${\displaystyle f_{t}=[6,\;-1,\;-1]}$? What is the normalized crosscorrelation in part (b)? What do you conclude from the magnitude of the largest value of this normalized crosscorrelation?

Solution

To normalize ${\displaystyle \phi _{ff}(\tau )}$, we divide by ${\displaystyle \phi _{ff}(0)=38}$ [see equation (9.8g)] from part (d), giving

{\displaystyle {\begin{aligned}\phi _{ff}(\tau )_{\mathrm {norm} }=\left[{\frac {-3}{19}},\;{\frac {-5}{38}},1,\;{\frac {-5}{38}},\;{\frac {-3}{19}}\right].\end{aligned}}}

To normalize ${\displaystyle \phi _{gf}(\tau )}$, we divide by ${\displaystyle [\phi _{ff}(0)\phi _{gg}(0)]^{1/2}}$ [see equation (9.8h)]. We use equation (9.8f) to get ${\displaystyle \phi _{gg}(0):\phi _{gg}(0)=2^{2}+5^{2}+2^{2}+1^{2}=34}$. Using the values of ${\displaystyle \phi _{gf}(\tau )}$ from part (b), we get

{\displaystyle {\begin{aligned}\phi _{gf}(\tau )_{\mathrm {norm} }=\phi _{gf}(\tau )/[34\times 38]^{1/2}=\phi _{gf}(\tau )/35.94\approx \phi _{gf}(\tau )/36,\\{\mbox{so}}\qquad \qquad \phi _{gf}(\tau )_{\mathrm {norm} }\approx \left[{\frac {-1}{18}},\;{\frac {-7}{36}},\;{\frac {\mathop {1} \limits ^{\downarrow }}{4}},\;{\frac {31}{36}},\;{\frac {-13}{36}},\;{\frac {1}{6}}\right]\\\approx [-0.056,\;-0.194,\;\mathop {0.250} \limits ^{\downarrow },\;0.861\;,\;0.361\;,\;0.167].\end{aligned}}}

The largest value comes when ${\displaystyle \tau =-1}$, so the two wavelets are most similar when ${\displaystyle f_{t}}$ is shifted one time unit to the right with respect to ${\displaystyle g_{t}}$. The maximum value is ${\displaystyle 31/36=0.86}$, which means that ${\displaystyle f_{t}}$ and ${\displaystyle g_{t}}$ are quite similar, but not the same.

## 9.12 Properties of minimum-phase wavelets

9.12a Of the four wavelets given in problem 9.8, which are minimum-phase?

Background

Minimum-phase wavelets are discussed briefly in problem 9.11 and in more detail in Sheriff and Geldart, 1995, section 9.4 and section 15.5.6); ${\displaystyle z}$-transforms are discussed in Sheriff and Geldart, 1995, section 15.5.3.

Solution

The four wavelets are: ${\displaystyle a_{t}=[2,\;-1]}$, ${\displaystyle b_{t}=[4,1]}$, ${\displaystyle c_{t}=[6,\;-7,\;2]}$, and ${\displaystyle d_{t}=[4,9,2]}$. The roots of the wavelets ${\displaystyle a_{t}}$ and ${\displaystyle b_{t}}$ are 2 and ${\displaystyle -4}$, so both are minimum-phase because the magnitudes of the roots are greater than unity. Since ${\displaystyle c_{t}\leftrightarrow (6-7z+2z^{2})=(3-2z)(2-z)}$, the roots are 3/2 and 2; ${\displaystyle c_{t}}$ is therefore also minimum-phase. Finally, ${\displaystyle d_{t}\leftrightarrow (4+9z+2z^{2})=(4+z)(1+2z)}$, the roots being ${\displaystyle -4}$ and ${\displaystyle -1/2}$. Because ${\displaystyle |-1/2|<1}$, ${\displaystyle d_{t}}$ is mixed-phase.

9.12b Find ${\displaystyle a_{t}*b_{t}}$ and ${\displaystyle a_{t}*c_{t}}$ by calculating in the time domain.

Solution

{\displaystyle {\begin{aligned}a_{t}*b_{t}=[2,\;-1]*[4,\;1]\;=8,\;-4\\2,\;-1\\=[8,\;-2,\;-1];\end{aligned}}} {\displaystyle {\begin{aligned}a_{t}*c_{t}=[2,\;-1]*[6,\;-7,\;2]\;=12,\;-6\\-14,\;7\\4,\;-2\\=[12,\;-20,\;11,\;2].\end{aligned}}}

9.12c Repeat part (b) except using transforms.

Solution

{\displaystyle {\begin{aligned}A(z)=(2-z)\;,\;B(z)=(4+z)\;,\;C(z)=(6-7z+2z^{2});\\a_{t}*b_{t}\leftrightarrow (2-z)(4+z)=8-2z-z^{2}\leftrightarrow [8,\;-2,\;-1];\\a_{t}*c_{t}\leftrightarrow (2-z)(6-7z+2z^{2})=12-20z+11z^{2}-2z^{3}\\=[12,\;-20,\;11,\;-2].\end{aligned}}}

9.12d Find ${\displaystyle a_{t}*b_{t}*c_{t}}$.

Solution

{\displaystyle {\begin{aligned}a_{t}*b_{t}*c_{t}=(a_{t}*b_{t})*c_{t}=[8,\;-2,\;-1]*[6,\;-7,\;2]\\=48,\;-12,\;-6\\-56,\;14,\;7\\16,\;-4,\;-2\\=[48,\;-68,\;24,\;3,\;-2].\\\end{aligned}}}

9.12e Does the largest value of a minimum-phase wavelet have to come at ${\displaystyle t=0}$?

Solution

The wavelet ${\displaystyle (a+z)(b+z)=[ab+(a+b)z+z^{2}]}$ is minimum-phase if ${\displaystyle |a|>1}$ and ${\displaystyle |b|>1}$. The ratio of the first two terms is ${\displaystyle ab/(a+b)=a/(1+a/b)}$ and it has its minimum absolute value when ${\displaystyle a}$ and $\displaystyle b$ have the same signs. When ${\displaystyle a}$ and ${\displaystyle b}$ have the same sign and are both slightly larger than unity, the ratio is close to 1/2 and the second term is larger than the first. As $\displaystyle a$ and/or ${\displaystyle b}$ increase, the ratio increases; the first and second terms are equal when $\displaystyle a=b=2$ . If $\displaystyle a$ and $\displaystyle b$ differ significantly in magnitude, the second term can be larger than the first for large values of $\displaystyle a$ or ${\displaystyle b}$; e.g., if ${\displaystyle b\approx 1}$, the ratio is ${\displaystyle \approx 0.9}$ when $\displaystyle a=9$ .

If ${\displaystyle a}$ and $\displaystyle b$ have opposite signs, the ratio cannot be smaller than 1 since the two terms in the denominator have opposite signs and the denominator cannot exceed the numerator.

When the wavelet has three factors, the ratio of the first to second term takes the form ${\displaystyle abc/(ab+ac+bc)}$. When ${\displaystyle a}$, ${\displaystyle b}$, and ${\displaystyle c}$ are all close to unity and of the same polarity, the magnitude of the ratio is ${\displaystyle \approx 1/3}$ and the second term is larger than the first. Generalizing to $\displaystyle n$ terms the ratio can be ${\displaystyle \approx 1/n}$.

9.12f Can a minimum-phase wavelet be zero at ${\displaystyle t=0}$?

Solution

If a wavelet is zero at ${\displaystyle t=0}$, it is of the form (${\displaystyle 0,a,b,c,\ldots }$}, so

{\displaystyle {\begin{aligned}\left[0,\;a,\;b,\;c,\;.\;.\;.\right]\;\leftrightarrow (0+az+bz^{2}+cz^{3}+\cdots \cdot )=z(a+bz+cz^{2}+\cdots ).\end{aligned}}}

Since one of the roots is ${\displaystyle z=0<1}$, the wavelet is not minimum-phase.

When we deal with an individual wavelet, we avoid the root ${\displaystyle z=0}$ by taking $\displaystyle t=0$ as the time when the first nonzero value occurs.

## 9.13 Phase of composite wavelets

Using the wavelets

{\displaystyle {\begin{aligned}W_{1}(z)=(2-z)^{2}(3-z)^{2};\quad W_{2}(z)=(4-z^{2})(9-z^{2}),\end{aligned}}}

calculate the composite wavelets:

\displaystyle \begin{align} W_{1}(z)+W_{2}(z) ;\quad W_{1}(z)+zW_{2}(z) ;\quad zW_{1}(z)+W_{2}(z) ;\quad z^{-1} W_{1}(z)+W_{2}(z). \end{align}

Plot the composite wavelets in the time domain. All of these composite wavelets have the same frequency spectrum but different phase spectra because multiplication by ${\displaystyle z}$ shifts the phase. The results illustrate the effects of phase.

Background

It is shown in Sheriff and Geldart, 1995, section 15.1.5 that the phase of a complex quantity, ${\displaystyle z=a+jb}$, is ${\displaystyle \gamma ={\rm {tan}}^{-1}(b/a)}$, that is, ${\displaystyle {\rm {\;tan\;}}\gamma =}$ (imaginary part)/(real part). Since ${\displaystyle z^{n}=e^{-{\hbox{j}}\omega n\Delta }={\rm {\;cos\;}}(\omega n\Delta )+j{\rm {\;sin\;}}(\omega n\Delta )}$, ${\displaystyle {\rm {\;tan\;}}\gamma ={\rm {\;tan\;}}(\omega n\Delta )}$. If a wavelet has several elements, the imaginary part of the wavelet will be the sum of several sines and the real part will be the sum of several cosines. When the wavelet is multiplied by $\displaystyle z^{n}$ , both sums will change, hence the phase changes.

Solution

{\displaystyle {\begin{aligned}W_{1}(z)=(2-z)^{2}(3-z)^{2}=(4-4z+z^{2})(9-6z+z^{2})\\=36-60z+37z^{2}-10z^{3}+z^{4}\leftrightarrow [36,\;-60,\;37.\;-10,\;1],\\W_{2}(z)=(4-z^{2})(9-z^{2})=36-13z^{2}+z^{4}\leftrightarrow [36,\;0,\;-13,\;0,\;1],\\W_{3}(z)=W_{1}(z)+W_{2}(z)=72-60z+24z^{2}-10z^{3}+2z^{4}\\\leftrightarrow [72,\;-60,\;24,\;-10,\;2],\\W_{4}(z)=W_{1}(z)+zW_{2}(z)=36-24z+37z^{2}-23z^{3}+z^{4}+z^{5}\\\leftrightarrow [36,\;-24,\;37,\;-23,\;1,\;1],\\W_{5}(z)=zW_{1}(z)+W_{2}(z)=36+36z-73z^{2}+37z^{3}-9z^{4}+z^{5}\\\leftrightarrow [36,\;36,\;-73,\;37,\;-9,\;1],\\W_{6}(z)=z^{-1}W_{1}(z)+W_{2}(z)=36z^{-1}-24+37z-23z^{2}+z^{3}+z^{4}\\\leftrightarrow [36,\;\mathop {24} \limits ^{\downarrow },\;-37,\;-23,\;-1,\;1].\\\end{aligned}}}

The convention is that (unless otherwise specified) the first nonzero element fixes the wavelet origin ${\displaystyle t=0}$ (see problem 9.12f). We start our plots in Figure 9.13a at ${\displaystyle t=-1}$ [except for ${\displaystyle W_{6}(z)]}$ in order to show the complete waveforms.

The two wavelets $\displaystyle W_{1}$ and ${\displaystyle W_{2}}$ are plotted first in Figure 9.13a and then the four sums of the wavelets. Note that two of the four factors that make up each of $\displaystyle W_{1}$ and ${\displaystyle W_{2}}$ are the same and the other two differ only in signs:

{\displaystyle {\begin{aligned}W_{1}(z)=(2-z)(2-z)(3-z)(3-z);\\W_{2}(z)=(2-z)(2+z)(3-z)(3+z).\\\end{aligned}}}

Figure 9.13a.  The composite wavelets.

However, the waveshapes differ significantly, especially in apparent frequency. Clearly relatively small changes in the equation of a wavelet can produce significant changes in the waveshape.

Because delaying the second wavelet or advancing the first produces the same waveshape, composite wavelets ${\displaystyle W_{4}}$ and ${\displaystyle W_{6}}$ are the same except for a time shift. Wavelets ${\displaystyle W_{3}}$ and ${\displaystyle W_{5}}$ have distinctly different waveshapes from those of ${\displaystyle W_{4}}$. We conclude that a shift of one component relative to the other has a significant effect on the location of peaks and troughs and hence on the waveshape.

## 9.14 Tuning and waveshape

9.14a The following wavelet (upper curve in Figure 9.14a) is approximately minimum-phase: ${\displaystyle \left[11,\;14,\;5,\;-10,\;-12,\;-6,\;3,\;5,\;2,\;0,\;-1,\;-1,\;0\right]}$, the sampling interval being 2 ms. Use ${\displaystyle V_{sd}=2.0}$ km/s for the velocity in sand, ${\displaystyle V_{sh}=1.5}$ km/s for the velocity in shale, and the shale-to-sand reflection coefficient ${\displaystyle =0.1}$ (value scaled up and rounded off). Determine the reflected waveshape for sands 0, 2, 4, 6, 8, and 10 m thick encased in shale. (A thickness of 6 m is approximately a quarter wavelength.)

Background

Reflection and transmission coefficients are discussed in problem 3.6, minimum-phase and zero-phase wavelets in Sheriff and Geldart, 1995, section 15.5.6a,d, respectively.

Tuning, the build up because of constructive interference that occurs when the travelpath difference between the waves reflected at top and bottom of a bed produces a half-cycle shift, which occurs at a thickness of ${\displaystyle \lambda /4}$ when reflection coefficients are of opposite polarity. The opposite effect of lowering the amplitude (detuning) occurs at a thickness of ${\displaystyle \lambda /4}$ if the reflections from top and base of the bed have the same polarity.

Solution

Because ${\displaystyle R=+0.1}$ for a shale-to-sand interface, the reflection at the sand-to-shale interface at the base of the sand is reversed in polarity compared to that from the top. The reflection energy is ${\displaystyle R^{2}=1\%}$ so 99% of the incident energy is transmitted into the sand; therefore we ignore transmission losses. The wave reflected at the base of the bed is delayed by the two-way traveltime, ${\displaystyle \tau =2z/V_{sd}=2z/2000=z}$ ms. Thus the reflection from the base is delayed by one time sample for each 2 m of sand thickness, i.e., ${\displaystyle \tau =z/2}$ time samples. Since the reflection coefficient is ${\displaystyle \pm 0.1}$ for both interfaces, we can take 0.1 as a scale factor and omit it in our calculations.

Figure 9.14a.  Wavelets used: (i) Minimum-phase wavelet [parts (a) to (c)]; (ii) low-frequency minimum-phase wavelet (part d); (iii) zerophase wavelet (part e).

To get the composite wave, that is, the sum of the wavelets ${\displaystyle g_{t}}$ and ${\displaystyle -g_{t-\tau }}$, we displace ${\displaystyle -g_{t}}$ by the amount ${\displaystyle \tau }$ and add it to ${\displaystyle g_{t}}$. When ${\displaystyle z=0}$, there is no reflection. The embedded wavelet is wavelet (i) shown in Figure 9.14a. Calculations for different sand thicknesses are shown in Table 9.14a. The embedded wavelet ${\displaystyle g_{t}}$ and the composite reflections are plotted in Figure 9.14b. The curves for ${\displaystyle \tau =1}$, 2, 3 are roughly displaced copies of each other and have successively larger amplitudes. The ${\displaystyle \tau =3}$ curve at the tuning thickness has the maximum amplitude, the largest amplitudes of the ${\displaystyle \tau =2}$ and ${\displaystyle \tau =4}$ curves being slightly smaller. If measuring arrival time by timing the first peak, one would pick too early when the bed is thinner than the tuning thickness. The first trough of the curve for ${\displaystyle \tau =4}$ is broadened and a change of phase is clear in the first trough of the ${\displaystyle \tau =5}$ curve, indicating that the bed is thicker than the resolvable limit of ${\displaystyle \lambda /4}$.

 For ${\displaystyle \tau =1(2m)}$ ${\displaystyle g_{t}}$ 11, 14, 5, –10, –12, –6, 3, 5, 2, 0, –1, –1, 0 ${\displaystyle -g_{t-1}}$ –11, –14, –5, 10, 12, 6, –3, –5, –2, 0, 1, 1, 0 Sum 11, 3, –9, –15, –2, 6, 9, 2, –3, –2 –1, 0, 1, 0 For ${\displaystyle \tau =2(4m)}$ ${\displaystyle g_{t}}$ 11, 14, 5, –10, –12, –6, 3, 5, 2, 0, –1, –1, 0 ${\displaystyle -g_{t-2}}$ –11, –14, –5, 10, 12, 6, –3, –5, –2, 0, 1, 1, 0 Sum 11, 14, –6, –24, –17, 4, 15, 11, –1, –5, –3, –1, 1, 1, 0 For ${\displaystyle \tau =3(6m\approx \lambda /4)}$ ${\displaystyle g_{t}}$ 11, 14, 5, –10, –12, –6, 3, 5, 2, 0, –1, –1, 0 ${\displaystyle -g_{t-3}}$ –11, –14, –5, 10, 12, 6, –3, –5, –2, 0, 1, 1, 0 Sum 11, 14, 5, –21, –26, –11, 13, 17, 8, –3, –6, –3, 0 1, 1, 0 For ${\displaystyle \tau =4(8m)}$ ${\displaystyle g_{t}}$ 11, 14, 5, –10, –12, –6, 3, 5, 2, 0, –1, –1, 0 ${\displaystyle -g_{t-4}}$ –11, –14, –5, 10, 12, 6, –3, –5, –2, 0, 1, 1, 0 Sum 11, 14, 5, –10, –23, –20, –2, 15, 14, 6, –4, –6, –2 0, 1, 1 0 For ${\displaystyle \tau =5(10m)}$ ${\displaystyle g_{t}}$ 11, 14, 5, –10, –12, –6, 3, 5, 2, 0, –1, –1, 0 ${\displaystyle -g_{t-4}}$ –11, –14, –5, 10, 12, 6, –3, –5, –2, 0, 1, 1, 0 Sum 11, 14, 5, –10, –12, –17, –11, 0, 12, 12, 5, –4, –5, –2, 0, 1, 1, 0

9.14b Repeat for the sand overlain by shale and underlain by limestone. Assume a sand-limestone reflection coefficient of ${\displaystyle +0.1}$.

Solution

The calculations (see Table 9.14b) are the same as in part (a) except that polarity at the sand/limestone interface is the same as at the shale/sand interface. The top curve in Figure 9.14b is for zero sand thickness so that the contact is shale/limestone with a reflection coefficient of ${\displaystyle +0.2}$. Timing the first peak would give erroneous depths for the ${\displaystyle \tau =1}$, 2, and 3 curves and the ${\displaystyle \tau =2,3}$ curves show rather clear phasing, indicating that more than one reflector is involved. The amplitude of the first trough decreases as the sand thickens for the ${\displaystyle \tau =1}$, 2, 3 curves and its character varies considerably for curves for ${\displaystyle \tau =3}$, 4, 5.

Figure 9.14b.  Shale/sand/shale reflections.

9.14c Determine the waveshape for two sands, each 2 m thick and separated by 1.5 m of shale, the sequence being encased in shale; this illustrates a “tuned” situation. Compare the results with those for 4 m and 6 m of sand in part (a), that is, for the same net and gross thicknesses.

Solution

The composite reflected wave is the sum ${\displaystyle g_{t}-g_{t-1}+g_{t-2}-g_{t-3}}$, calculated in Table 9.14c.

 ${\displaystyle g_{t}}$ 11, 14, 5, −10, −12, −6, 3, 5, 2, 0, −1, −1, 0 ${\displaystyle -g_{t-1}}$ −11, −14, −5, 10, 12, 6, −3, −5, −2, 0, 1, 1, 0 ${\displaystyle g_{t-2}}$ 11, 14, 5, −10, −12, −6, 3, 5, 2, 0, −1, −1, 0 ${\displaystyle -g_{t-3}}$ −11, 14, −5, 10, 12, 6, −3, −5, −2, 0, 1, 1, 0 Sum 11, 3, 2, −12, −10, −9, 7, 8, 6, 0, −4, −2, 0, 0, 1, 0
 For ${\displaystyle \tau =1}$ (2 m) ${\displaystyle g_{t}}$ 11, 14, 5, −10, −12, −6, 3, 5, 2, 0, −1, −1, 0 ${\displaystyle g_{t-1}}$ 11, 14, 5, −10, −12, −6, 3, 5, 2, 0, −1, −1, 0 Sum 11, 25, 19, −5, −22, −18, −3, 8, 7, 2 −1, −2, −1, 0 For ${\displaystyle \tau =2}$ (4 m) ${\displaystyle g_{t}}$ 11, 14, 5, −10, −12, −6, 3, 5, 2, 0, −1, −1, 0 ${\displaystyle g_{t-2}}$ 11, 14, 5, −10, −12, −6, 3, 5, 2, 0, −1, −1, 0 Sum 11, 14, 16, 4, −7, −16, −9, −1, 5, 5, 1, −1, −1, −1, 0 For ${\displaystyle \tau =3}$ (6 m ${\displaystyle \approx \lambda /4}$) ${\displaystyle g_{t}}$ 11, 14, 5, −10, −12, −6, 3, 5, 2, 0, −1, −1, 0 ${\displaystyle g_{t-3}}$ 11, 14, 5, −10, −12, −6, 3, 5, 2, 0, −1, −1, 0 Sum 11, 14, 5, 1, 2, −1, −7, −7, −4, 3, 4, 1, 0 −1, −1, 0 For ${\displaystyle \tau =4}$ (8 m) ${\displaystyle g_{t}}$ 11, 14, 5, −10, −12, −6, 3, 5, 2, 0, −1, −1, 0 ${\displaystyle g_{t-4}}$ 11, 14, 5, −10, −12, −6, 3, 5, 2, 0, −1, −1, 0 Sum 11, 14, 5, −10, −1, 8, 8, −5, −10, −6, 2, 4, 2 0, −1, −1, 0 For ${\displaystyle \tau =5}$ (10 m) ${\displaystyle g_{t}}$ 11, 14, 5, −10, −12, −6, 3, 5, 2, 0, −1, −1, 0 ${\displaystyle g_{t-5}}$ 11, 14, 5, −10, −12, −6, 3, 5, 2, 0, −1, −1, 0 Sum 11, 14, 5, −10, −12, 5, 17, 10, −8, −12, −7, 2, 5, 2, 0, −1, −1, 0

Single sands 4 and 6 m thick from part (a):

 6 m 11, 14, 5, -21, -26, -11, 13, 17, 8, -3, -6, -3, 0, 1, 1, 0; 4 m 11, 14, -6, -24, -17, 4, 15, 11, -1, -5, -3, -1, 1, 1, 0.

The composite wavelet for two 2-m sands is the upper curve in Figure 9.14d and those for single 4-m and 6-m sands are the lower curves, i.e., for the same net and same gross sand thicknesses. The composite curve for the two thin sands shows phasing where the reflections from the top and base of each sand interfere. The gross thickness of 6 m is above the tuning thickness (note the peak-to-trough time difference between the 4-m and 6-m sands, evidence that 6 m is larger than the tuning thickness). Where the gross thickness of an interval is smaller than a quarter wavelength, information as to the thickness of the different lithology (sand, in this case) is contained in amplitude rather than in time measurements (or waveshape changes).

Figure 9.14c.  Shale/sand/limestone reflections.

9.14d Repeat part (a) with the low-frequency wavelet [6,11,14,14,10,5,–2,–10,–11,–12,–10,–6,0,3,4,5,4,3,1,0] [Figure 9.14a (ii)], which is the minimum-phase wavelet in Figure 9.14a(i), but stretched out so that it has about half the dominant frequency. Compare with the results of part (a) to illustrate the effect of frequency on resolution.

Solution

We proceed as in part (a) after replacing the former embedded wavelet with the new one. The results for the shale/sand/shale are shown in Table 9.14d and plotted in Figure 9.14e.

We compare Figures 9.14b and 9.14e which differ only in the frequencies, both having the same embedded wavelet waveshape. In Figure 9.14e the curves all have nearly the same shape but differ in amplitude except for the ${\displaystyle \tau =5}$ curve where the trough and second peak are broadened because ${\displaystyle \tau =5}$ is about equal to ${\displaystyle \lambda /4}$. The resolution is poorer than with the higher-frequency wavelet in Figure 9.14a(i).

Figure 9.14d.  Effect of thin sands.
 For ${\displaystyle \tau =1}$ (2 m) ${\displaystyle g_{t}}$ 6, 11, 14, 14, 10, 5, –2, –10, –11, –12, –10, –6, 0, 3, 4, 5, 4, 3, 1, 0 ${\displaystyle -g_{t-1}}$ –6, –11, –14, –14, –10, –5, 2, 10, 11, 12, 10, 6, 0, –3, –4 –5, –4, –3, –1, 0 Sum 6, 5, 3, 0, 4, 5, 7, 8, 1, 1, 2, 4, 6, 3, 1, 1, 1, 1, 2, 1, 0 For ${\displaystyle \tau =2}$ (4 m) ${\displaystyle g_{t}}$ 6, 11, 14, 14, 10, 5, –2, –10, –11, –12, –10, –6, 0, 3, 4, 5, 4, 3, 1, 0 ${\displaystyle -g_{t-3}}$ –6, –11, –14, –14, –10, –5, 2, 10, 11, 12, 10, 6, 0, –3, –4 –5, –4, –3, –1, 0 Sum 6, 11, 8, 3, 4, 9, 12, 15, 9, 2, 1, 6, 10, 9, 4, 2, 0, 2, 3, 3, 1, 0 For ${\displaystyle \tau =3}$ (6 m) ${\displaystyle g_{t}}$ 6, 11, 14, 14, 10, 5, –2, –10, –11, –12, –10, –6, 0, 3, 4, 5, 4, 3, 1, 0 ${\displaystyle -g_{t-3}}$ –6, –11, –14, –14, –10, –5, 2, 10, 11, 12, 10, 6, 0, –3, –4, –5, –4, –3, –1, 0 Sum 6, 11, 14, 8, 1, 9, 16, 20, 16, 10, 0, 5, 12, 13, 10, 5, 1, 1, 4, 4, 3, 1, 0 For ${\displaystyle \tau =4}$ (8 m) ${\displaystyle g_{t}}$ 6, 11, 14, 14, 10, 5, –2, –10, –11, –12, –10, –6, 0, 3, 4, 5, 4, 3, 1, 0 ${\displaystyle -g{t-4}}$ –6, –11, –14, –14, –10, –5, 2, 10, 11, 12, 10, 6, 0, –3, –4, –5, –4, –3, –1, 0 Sum 6, 11, 14, 14, 4, –6, –16, –24, –21, –17, –8, 4, 11, 15, 14, 11, 4, 0, –3, –5, –4, –3, –1, 0 For ${\displaystyle \tau =5}$ (10 m) ${\displaystyle g_{t}}$ 6, 11, 14, 14, 10, 5, –2, –10, –11, –12, –10, –6, 0, 3, 4, 5, 4, 3, 1, 0 ${\displaystyle -g_{t-5}}$ –6, –11, –14, –14, –10, –5, 2, 10, 11, 12, 10, 6, 0, –3, –4, –5, –4, –3, –1, 0 Sum 6, 11, 14, 14, 10, –1, –13, –24, –25, –22, –15, –4, 10, 14, 16, 15, 10, 3, –2, –4, –5, –4, –3, –1, 0

9.14e Repeat parts (a) and (b) using the zero-phase wavelet [1, 1, ${\displaystyle -1}$, ${\displaystyle -4}$, ${\displaystyle -6}$, ${\displaystyle -4}$, 10, 17, 10, ${\displaystyle -4}$, ${\displaystyle -6}$, ${\displaystyle -4}$, ${\displaystyle -1}$, 1, 1] [Figure 9.14a(iii)], which has nearly the same spectrum as wavelet (i).

Solution

The calculations are shown in Tables 9.14e and 9.14f and the reflection waveshapes in Figures 9.14f and 9.14g.

These curves have been plotted about their points of symmetry or asymmetry. The top curve in Figure 9.14f is the embedded wavelet, in Figure 9.14g it is for zero sand thickness. As with the minimum-phase wavelet, timing would be in error where the thickness is less than ${\displaystyle \lambda /4}$ and the maximum amplitude in Figure 9.14f occurs where the thickness is ${\displaystyle \lambda /4}$. The peak-trough times increase for ${\displaystyle \tau =4}$ and ${\displaystyle \tau =5}$ in Figure 9.14f, and there is distinct phasing for $\displaystyle \tau=5$ . In Figure 9.14g, there is a distinct change of shape at $\displaystyle \tau=3$ and there is clear resolution and an amplitude minimum at ${\displaystyle \tau =3}$.

Figure 9.14e.  Shale/sand/shale reflection with low-frequency wavelet in Figure 9.14a(ii).
Figure 9.14f.  Shale/sand/shale reflection, zero-phase wavelet.
Figure 9.14g.  Shale/sand/limestone reflection, zero-phase wavelet.
 For ${\displaystyle \tau =1}$ (2 m) ${\displaystyle g_{t}}$ 1, 1, –1, –4, –6. –4, 10, 17, 10, –4, –6, –4, –1, 1, 1, 0 ${\displaystyle -g_{t-1}}$ 1, 1, 1, 4, 6. 4, –10, –17, –10, 4, 6, 4, 1, –1, –1, 0 Sum 1, 0, –2, –3, –2, 2, 14, 7, –7, –14, –2, 2, 3, 2, 0, –1, 0 For ${\displaystyle \tau =2}$ (4 m) ${\displaystyle g_{t}}$ 1, 1, 1, –4, –6. –4, 10, 17, 10, –4, –6, –4, 1, 1, 1, 0 ${\displaystyle -g_{t-2}}$ –1, –1, 1, 4, 6. 4, –10, –17, –10, 4, 6, 4, 1, –1, –1, 0 Sum 1, 1, –2, –5, –5, 0, 16, 21, 0, –21, –16, 0, 5, 5, 2, –1, –1, 0 For ${\displaystyle \tau =3}$ (6 m) $\displaystyle g_{t}$ 1, 1, –1, –4, –6. –4, 10, 17, 10, –4, –6, –4, –1, 1, 1, 0 ${\displaystyle -g_{t-3}}$ –1, –1, 1 4, 6. 4, –10, –17, –10, 4, 6, 4, 1, –1, –1, 0 Sum 1, 1, –1, –5, –7, –3, 14, 23, 14, –14, –23, –14, 3, 7, 5, 1, –1, –1, 0 For ${\displaystyle \tau =4}$ (8 m) ${\displaystyle g_{t}}$ 1, 1, –1, –4, –6. –4, 10, 17, 10, –4, –6, –4, –1, 1, 1, 0 ${\displaystyle -g_{t-4}}$ –1, –1, 1 4, 6. 4, –10, –17, –10, 4, 6, 4, 1, –1, –1, 0 Sum 1, 1, –1, –4, –7, –5, 11, 21, 16, 0, –16, –21, –11, 5, 7, 4, 1, –1, –1, 0 For ${\displaystyle \tau =5}$ (10 m) ${\displaystyle g_{t}}$ 1, 1, –1, –4, –6. –4, 10, 17, 10, –4, –6, –4, –1, 1, 1, 0 $\displaystyle -g_{t-4}$ –1, –1, 1 4, 6. 4, –10, –17, –10, 4, 6, 4, 1, –1, –1, 0 Sum 1, 1, –1, –4, –6, –5, 9, 18, 14, 2, –2, –14, –18, –9, 4, 6, 4, 1, –1, –1, 0
 For ${\displaystyle \tau =1}$ (2 m) ${\displaystyle g_{t}}$ 1, 1, -1, -4, -6. -4, 10, 17, 10, -4, -6, -4, -1, 1, 1, 0 $\displaystyle -g_{t-1}$ 1, 1, -1, -4, -6. -4, 10, 17, 10, -4, -6, -4, -1, 1, 1, 0 Sum 1, 2, 0, -5, -10, -10, 6, 27, 27, 6, -10, -10, -5, 0, 2, 1, 0 For ${\displaystyle \tau =2}$ (4 m) ${\displaystyle g_{t}}$ 1, 1, -1, -4, -6. -4, 10, 17, 10, -4, -6, -4, -1, 1, 1, 0 ${\displaystyle -g_{t-2}}$ 1, 1, -1, -4, -6, -4, 10, 17, 10, -4, -6, -4, -1, 1, 1, 0 Sum 1, 1, 0, -3, -7, -8, 4, 13, 20, 13, 4, -8, -7, -3, 0, 1, 1, 0 For $\displaystyle \tau = 3$ (6 m) ${\displaystyle g_{t}}$ 1, 1, -1, -4, -6. -4, 10, 17, 10, -4, -6, -4, -1, 1, 1, 0 $\displaystyle -g_{t-3}$ 1, 1, -1, -4, -6. -4, 10, 17, 10, -4, -6, -4, -1, 1, 1, 0 Sum 1, 1, -1, -3, -5, -5, 6, 11, 6, 6, 11, 6, -5, -5, -3, -1, 1, 1, 0, For ${\displaystyle \tau =4}$ (8 m) ${\displaystyle g_{t}}$ 1, 1, -1, -4, -6. -4, 10, 17, 10, -4, -6, -4, -1, 1, 1, 0 ${\displaystyle -g_{t-4}}$ 1, 1, -1, -4, -6. -4, 10, 17, 10, -4, -6, -4, -1, 1, 1, 0 Sum 1, 1, 1, 4, 5, 3, 9, 13, 4, 0, 4, 13, 9, 3, 5, 4, 1, 1, 1, 0 For ${\displaystyle \tau =5}$ (10 m) $\displaystyle g_{t}$ 1, 1, -1, -4, -6. -4, 10, 17, 10, -4, -6, -4, -1, 1, 1, 0 ${\displaystyle -g_{t-5}}$ 1, 1, -1, -4, -6. -4, 10, 17, 10, -4, -6, -4, -1, 1, 1, 0 Sum 1, 1, -1, -4, -6, -3, 11, 16, 6, -10, -10, 6, 16, 11, -3, -6, -4, -1, 1, 1, 0

## 9.15 Making a wavelet minimum-phase

The wavelet ${\displaystyle [-0.9505,\;-0.0120,\;0.9915]}$ is not minimum-phase. How would you change it to make it mimimum-phase without changing it any more than necessary? Give two methods.

Background

Minimum-phase wavelets are discussed in Sheriff and Geldart, 1995, section 9.4 and section 15.5.6a.

Solution

The transform of the wavelet is ${\displaystyle W(z)=-0.9505-0.0120z+0.9915z^{2}}$. The roots of ${\displaystyle W(z)=0}$ are 0.9652 and ${\displaystyle -0.9531}$, so to make the wavelet minimum-phase, we must change it so that the magnitudes of both roots are greater than unity. To increase both roots by the same factor, we multiply both roots by a number slightly greater than ${\displaystyle 1/0.9531=1.0492}$. If we use the multiplier 1.0500 the roots become 1.0135 and ${\displaystyle -1.00076}$ and the revised equation is

{\displaystyle {\begin{aligned}(z-1.0135)(z+1.0008)=-1.0143-0.0127z+z^{2}.\end{aligned}}}

To obtain the same value at $\displaystyle t=0$ as we had before, we multiply by ${\displaystyle 0.9505/1.0143=0.9371}$ and the transform becomes ${\displaystyle 0.9505-0.0119z+0.9371z^{2}}$ which gives the wavelet [${\displaystyle -0.9505}$, ${\displaystyle -0.0119}$, 0.9371]. The modified wavelet is nearly the same as the original except that its “tail” has been decreased from 0.9915 to 0.9371, a 5.5% decrease.

A second method is to use a taper to reduce the tail of the wavelet. A commonly used taper is ${\displaystyle k^{t}}$ where ${\displaystyle t}$ is in milliseconds and ${\displaystyle k}$ is slightly less than unity, for example, 0.9950. Using this value for ${\displaystyle k}$ and assuming that the sampling interval is 2 ms, the values of the taper are 1.0000, ${\displaystyle 0.9950^{2}}$, ${\displaystyle 0.9950^{4}}$, that is, 1.0000, 0.9900, 0.9801 and we get ${\displaystyle W(z)=-0.9505-0.9900\times 0.0120z+0.9801\times 0.9915z^{2}=-0.9505-0.0119z+0.9718z^{2}}$. The roots of ${\displaystyle W(z)=0}$ are now 0.9904, ${\displaystyle -0.9783}$, both less than unity so we did not apply enough taper. We next try using ${\displaystyle 0.9850^{t}}$, which reduces the wavelet to ${\displaystyle [-0.9505,\;-0.0116,\;0.9333]}$. The roots are now 1.0154 and ${\displaystyle -1.0030}$ and ${\displaystyle W(z)=(z-1.0154)(z+1.0030)=-1.0184-0.0124z+z^{2}}$. Multiplying by $\displaystyle 0.9505/1.0184=0.9333$ makes the wavelet [0.9505, $\displaystyle -0.0116$ , 0.9333].

Comparing the two methods, we see that the second method has changed the wavelet slightly more than the first. In practice, a wavelet will have many more than three elements and the second method will generally be more practical. To verify that the taper is large enough to achieve its objective, we have to find all of the roots.

## 9.16 Zero-phase filtering of a minimum-phase wavelet

Show that the result of passing a minimum-phase signal through a zero-phase filter is mixed phase.

Background

Zero-phase signals are discussed in Sheriff and Geldart, 1995, section 15.5.6d, where it is shown that the spectrum of a zero-phase signal comprises products of pairs of factors of the form ${\displaystyle (az-1)(az^{-1}-1)=(1+a^{2})-a(z+z^{-1})=(1+a^{2})-2a{\rm {\;cos\;}}(\omega \Delta )}$, where ${\displaystyle a}$ can be complex. Because the imaginary part is zero, the phase is zero.

Solution

Let ${\displaystyle \omega _{t}}$ be the minimum-phase signal and ${\displaystyle f_{t}}$ the zero-phase filter. Time-domain filtering is accomplished by convolution, $\displaystyle \omega _{t} *f_{t}$ ; in the frequency domain the result is $\displaystyle W(z)F(z)$ . The factors of ${\displaystyle F(z)}$ occur in pairs of the form ${\displaystyle (az-1)(az^{-1}-1)}$, and each pair has roots ${\displaystyle z=a}$, ${\displaystyle 1/a}$. If ${\displaystyle |a|>1}$, then ${\displaystyle |1/a|<1}$. Thus, one member of each pair of roots is not minimum-phase and consequently the filtered signal is mixed-phase.

## 9.17 Deconvolution methods

Deconvolution, whose objective is undoing the results of a prior convolution, is a somewhat open-ended collection of methods, a number of which are briefly described in the background. List the assumptions of different methods, such as invariant wavelet, randomness of the reflectivity or of the noise, whether a source wavelet is the same as the recorded wavelet, and so on.

Background

Various deconvolution methods and their characteristics are described in the following text. Instead of giving a solution, we leave it up to the reader to list the assumptions.

A change in signal waveshape can be regarded as filtering, the effect being equivalent to convolution of the signal ${\displaystyle g_{t}}$ with a filter, ${\displaystyle h_{t}=f_{t}*g_{t}}$. Deconvolution attempts to remove the effects of ${\displaystyle f_{t}}$ and obtain the original signal ${\displaystyle g_{t}}$. Deconvolution in the time domain consists of convolution with an inverse filter ${\displaystyle i_{t}}$, i.e., $\displaystyle h_{t} *i_{t}$ (see problems 9.18 and 9.19); in the frequency doman its equivalent is multiplication of the transforms, ${\displaystyle H(\omega )I(\omega )}$ (see Sheriff and Geldart, 1995, section 9.5 for more details). An inverse filter convolved with the filter yields an impulse, $\displaystyle i_{t} *f_{t} \leftrightarrow I(\omega)F(\omega)=1\leftrightarrow \delta_{t}$ . The objective of deconvolution often is to determine the shape of the embedded wavelet (see problem 9.6).

Most deconvolution methods are based on autocorrelations of individual traces. An autocorrelation measures the repetition in a time series. Presumably the embedded wavelet is repeated for every reflection and thus the early part of the autocorrelation is largely determined by the shape of the embedded wavelet (see Figure 9.17a). An embedded wavelet that is long will contribute significantly to several half-cycles of the autocorrelation. Because we generally want a short embedded wavelet, we want the autocorrelation’s amplitude to die off rapidly.

The wavelet changes shape as it travels through the earth. Nevertheless, most deconvolution methods assume constant waveshape, sometimes called stationarity. Obtaining a representative autocorrelation requires that an appreciable number of time samples be included in the calculation process, so usually 500 or more samples (1 s at 0.002 s sampling) are included. Because there is no phase information in an autocorrelation, the phase spectrum of the embedded wavelet cannot be recovered from it. It is often assumed that the embedded wavelet is minimum-phase.

Changes in the wavelet shape with time are generally accommodated by subdividing a trace into time segments, e.g., finding a deconvolution operator for the early portion of a trace and another for a later portion, assuming that the wavelet shape is constant during each portion. The portions analyzed often overlap, e.g., one autocorrelation is calculated for the data between 0.5 and 2.0 s, a second for the data between 1.5 and 3.0 s. Each deconvolution operator is then applied over its respective range, giving two outputs in the overlap region. These two outputs are then added together in different proportions in the overlap region, the output at 1.5 s being 100% of that given by the early operator and that at 2.0 s being 100% of that given by the latter operator, the proportions changing linearly during the overlap time. This is often called adaptive deconvolution.

While many deconvolution methods exist, they can be roughly divided into deterministic and statistical methods, but the division is often unclear because deterministic methods may employ statistics and statistical methods may utilize knowledege about the nature of the convolution to be undone. Deterministic deconvolution requires that we know or can reasonably assume the mechanism or properties of the convolution to be undone.

One type of deconvolution is dereverberation or deringing, whose objective is to remove the effects of energy bouncing repeatedly in a near-surface layer—usually the water layer with marine data. This requires a knowledge of the repetition time of the near-surface bounces and the relative amplitudes of successive bounces. In the marine case the repetition time is assumed to be given by the water depth and the amplitudes by the sea-floor reflectivity, which sometimes can be obtained by trials. One marine deghosting method used with ocean-bottom recording employs velocity geophones and hydrophones (see problem 7.10).

Deghosting is a type of deterministic deconvolution where we assume that the ghost is a replica of the original signal ${\displaystyle G(z)}$ with amplitude reduced by the factor ${\displaystyle R}$ and delayed by ${\displaystyle n\Delta }$, where ${\displaystyle R}$ is the reflection coefficient at the ghosting interface (note that ${\displaystyle R}$ is usually negative) for a wave approaching from below and ${\displaystyle n\Delta }$ is the two-way traveltime between the source and the ghosting interface. The transform of the signal plus ghost is then

 {\displaystyle {\begin{aligned}H(z)=G(z)(1+Rz^{n}).\end{aligned}}} (9.17a)

On land we usually can get ${\displaystyle n}$ fairly accurately from the depth of the source and the thickness and velocity of the LVL, but we must determine (or guess) the value of ${\displaystyle R}$. In marine work, we get $\displaystyle n$ from the depth of the water. Deghosting is often done in a recursive manner, called recursive deconvolution.

System deconvolution is sometimes deterministically carried out to remove the filtering action of recording instrumentation or data processing.

In spiking deconvolution we assume that the impulse response of the earth $\displaystyle r_{t}$ , whose elements are the reflection coefficients at the various interfaces, are randomly distributed, and, hence, the autocorrelation of ${\displaystyle r_{t}}$ has a nonzero value (a “spike”) only at zero time shift, that is,

 {\displaystyle {\begin{aligned}\phi _{rr}(t)\approx k\delta _{t}.\end{aligned}}} (9.17b)

[“Random” here means unpredictable, i.e., one cannot predict the arrival time of a primary reflection based on the arrival times of earlier reflections.] If ${\displaystyle \omega _{t}}$ is the embedded wavelet, the geophone output is ${\displaystyle g_{t}=e_{t}*\omega _{t}}$. We can use the method of least squares (least-squares filtering is discussed in problem 9.22) to find the optimum inverse filter that will give a result that has the properties we assume for ${\displaystyle e_{t}}$ (see Sheriff and Geldart, 1995, 296). Spiking deconvolution can be carried out in either the time or frequency domain. Because we consider ${\displaystyle e_{t}}$ as random, its spectrum contains all frequencies in equal abundance, that is, its spectrum should be flat; techniques for achieving this are called spectral flattening or flattening deconvolution. Flattening is usually done only over the passband where the signal is assumed to be dominant. Since the autocorrelation of ${\displaystyle e_{t}}$ is small except for zero shift, the inverse filter ${\displaystyle I(z)}$ is

 {\displaystyle {\begin{aligned}I(z)=1/H(z),\end{aligned}}} (9.17c)

where ${\displaystyle H(z)}$ is the transform of the observed seismic trace. We can get the amplitude spectrum from the autocorrelation, but we also need to know the phase of ${\displaystyle H(z)}$ in order to solve the problem. There is no phase information in the autocorrelation, and so we have to assume the phase to get a solution. We usually assume that ${\displaystyle H(z)}$ is minimum-phase. Spiking deconvolution is often done to shorten the duration of the embedded wavelet, but it may make the data noisy if too much noise is included.

Predictive deconvolution uses information about the primary reflections to predict multiples produced by the same reflectors. Long-path multiples cause systematic repetition of a trace which produces significant values of the autocorrelation following the time delay associated with the multiple repetition time (see Figure 9.17a). Hence deconvolution to remove multiples is generally based on the portion of an autocorrelation trace after a lag time $\displaystyle L$ . A predictive deconvolution filter does not begin to act until the time $\displaystyle L$ . Autocorrelation elements that occur after ${\displaystyle L}$ are regarded as produced by multiples and a least-squares method can be used to find the filter that will predict the multiples. The predicted multiples can then be subtracted to get rid of the recorded multiples.

Entropy is a measure of the disorder or unpredictability of a system, the entropy increasing as the disorder increases. The autocorrelation of a time series is not unique, and a spectrum corresponds to a number of different time series. Maximum-entropy filtering attempts to select the time series that has the maximum entropy (maximum disorder).

We sometimes try to represent a seismic trace as a sequence minimizing the number of reflection events, i.e., involving only a few large reflections. This can be accomplished, for example, by minimizing the error in a least-fourth power sense rather than in a least-squares sense (as is done in Wiener filtering). This type of operation is sometimes called minimum-entropy (sparse-spike) deconvolution.

Homomorphic deconvolution (Sheriff and Geldart, 1995, 298) involves transformation into the cepstral domain, the cepstrum ${\displaystyle g^{*}(\zeta )}$ being defined by the relations

{\displaystyle {\begin{aligned}{\rm {\;ln\;}}[G(z)]=G^{*}(z)\leftrightarrow g^{*}(\zeta ),\end{aligned}}}

Where $\displaystyle \leftrightarrow$ denotes the Fourier transformation. In the cepstral domain the geophone input $\displaystyle g_{t} =e_{t} *\omega _{t}$ becomes ${\displaystyle g^{*}(\zeta )=e^{*}(\zeta )+\omega ^{*}(\zeta )}$. If ${\displaystyle e_{t}}$ varies more rapidly than $\displaystyle \omega _{t}$ , the two may be separable by frequency filtering.

## 9.18 Calculation of inverse filters

Assuming that the signature of an air-gun array is a unit impulse and that the recorded wavelet after transmission through the earth is ${\displaystyle [-12,\;-4,\;+3,\;+1]}$, find the inverse filter that will remove the earth filtering. How many terms should the filter include?

Solution

The inverse filter ${\displaystyle i_{t}}$ (see problem 9.7) is a filter that will restore the source impulse, i.e.,

{\displaystyle {\begin{aligned}g_{t}*i_{t}=\delta _{t},\end{aligned}}}

or in the frequency domain where ${\displaystyle i_{t}\leftrightarrow I(z)}$,

\displaystyle \begin{align} G(z)I(z)=1. \end{align}

Thus,

 {\displaystyle {\begin{aligned}I(z)=1/G(z)=1/[-12-4z+3z^{2}+z^{3}]\\=-{\frac {1}{12}}\left[1+{\frac {z}{3}}-{\frac {z^{2}}{4}}-{\frac {z^{3}}{12}}\right]^{-1}=-{\frac {1}{12}}(1+A)^{-1},\end{aligned}}} (9.18a)

where ${\displaystyle A={\frac {z}{3}}-{\frac {z^{2}}{4}}-{\frac {z^{3}}{12}}}$.

Since ${\displaystyle z}$ has the magnitude ${\displaystyle |z|=1}$, the magnitude of ${\displaystyle A<1}$ for all values of ${\displaystyle z}$, and we can expand equation (9.18a) using the binomial theorem [see equation (4.1b)] and Sheriff and Geldart, 1995, equation (15.43):

 {\displaystyle {\begin{aligned}I(z)=-{\frac {1}{12}}(1-A+A^{2}-A^{3}+A^{4}-\ldots ).\end{aligned}}} (9.18b)

We first find ${\displaystyle (1-A+A^{2}-A^{3}+A^{4})}$ neglecting powers higher than ${\displaystyle z^{4}}$:

\displaystyle \begin{align} -A&=-0.3333z+0.2500z^{2} +0.0833z^{3} ,\\ A^{2} &= +0.1111z^{2} -0.1667z^{3} +0.0069z^{4} ,\\ -A^{3} &= -0.0370z^{3} +0.0833z^{4} ,\\ A^{4} &= +0.0123z^{4} \\ \mathrm{Sum} &=-0.3333z+0.3611z^{2} -0.1204z^{3} +0.1025z^{4} \\ I(z)&=-\frac{1}{12} (1-A+A^{2} -A^{3} +A^{4})\\ &=-\frac{1}{12} (1-0.3333z+0.3611z^{2} -0,1204z^{3} +0.1025z^{4}). \end{align}

We can verify the accuracy of ${\displaystyle I(z)}$ by multiplying ${\displaystyle G(z)}$ by ${\displaystyle I(z)}$. We have

\displaystyle \begin{align} I(z):1-0.3333z+0.3611z^{2} -0, 1204z^{3} +0.1025z^{4} ,\\ G(z):1+0.3333z-0.2500z^{2} -0.0833z^{3} \\ 1-0.3333z+0.3611z^{2} -0.1204z^{3} +0.1025z^{4} \\ +0.3333z-0.1111z^{2} +0.1204z^{3} -0.0401z^{4} +0.0342z^{5} \\ -0.2500z^{2} +0.0833z^{3} -0.0903z^{4} +0.0301z^{5} -0.0256z^{6}\\ {-0.0833z^{3} +0.0278z^{4} -0.0301z^{5} +0.0100z^{6} -0.0085z^{7}}\\ {1+0+0+0-0.0001z^{4} +0.0342z^{5} -0.0156z^{6} -0.0085z^{7}}. \end{align}

We see that the inverse filter is exact as far as the term $\displaystyle z^{3}$ and terms for higher powers are small. The overall effect is to create a small tail whose energy is 0.00149 or 0.1%.

To determine how the accuracy depends on the number of terms used in ${\displaystyle I(z)}$, we observe the effect on the pro