% copyright (c) 1998 Jon Claerbout\title{Spatial aliasing and scale invariance}\author{Jon Claerbout}\maketitle\label{paper:lal}\sx{alias}\sx{spatial aliasing}\sx{scale invariance}\parLandforms are not especially predictable.Therefore, crude PEF approximations are often satisfactory.Wavefields are another matter.Consider the ``shape'' of the acoustic wavefrontsat this moment in the room you are in.The acoustic wavefield has statistical order in many senses.If the 3-D volume is filled with waves emitted from a few point sources,then (with some simplifications)what could be a volume of information is actually a few 1-D signals.When we work with wavefronts we can hope for more dramatic,even astounding, results from estimating properly.%\par%In its simplest form,%the \bx{Nyquist} condition says that%we can have no frequencies higher than two points per wavelength.%In migration, this is a strong constraint on data collection.%It seems there is no escape.%Yet, in applications dealing with a CMP gather%% (such as in Figure~\ref{lal/conj-stretch} or \ref{lal/conj-stack}),%we see data with spatial frequencies that exceed Nyquist%and we are not bothered, because after \bx{NMO},%these frequencies are OK.%Nevertheless, such data is troubling because%it breaks many of our conventional programs,%such as downward continuation with finite differences%or with Fourier transforms.%(No one uses focusing for \bx{stacking}.)%Since NMO overcomes%the limitation imposed%by the simple statement of the Nyquist condition,%we revise the condition to say that%the real limitation is on the spectral bandwidth,%not on the maximum frequency.%Mr.~Nyquist does not tell us where that bandwidth must be located.%Further, it seems that precious bandwidth {\it need not be contiguous.}%The signal's spectral band can be split into pieces%and those pieces positioned in different places.%Fundamentally, the issue is whether the total bandwidth exceeds Nyquist.%%Noncontiguous Nyquist bands are depicted in Figure~\FIG{nytutor}.%%%\ACTIVESIDEPLOT{nytutor}{width=3in}{nytutor}{%%      Hypothetical spatial frequency bands.%%      Top is typical.%%      Middle for data skewed with $\tau=t-px$.%%      Bottom depicts data with wave arrivals%%      from three directions.%%      }%%%%\par%Noncontiguous bandwidth arises naturally with%two-dimensional data where there are several plane-waves present.%There the familiar spatial Nyquist limitation %oversimplifies real life because the plane-waves link time and space.\par\noindent\boxit{         The plane-wave model links an axis that is not aliased (time)        with axes (space) that often are.        }\par%%\subsection{Aliasing a plane-wave chirp}%Spatial aliasing is often described as being%like temporal aliasing but on another axis.%Mathematically,%aliasing on the two axes is independent,%however, the plane-wave model links them%in our day-to-day experience.%The plane wave model links an axis (time) that is not aliased%with axes (space) that might be.%\par%A function generally useful with seismic vibrator sources%is the ``gated-\bx{chirp} function''.%Within a specified time gate%this function is a sinusoid with a time-variable frequency.%Here we will take this signal to be defined by%\begin{equation}%f(\tau ) \quad = \quad%       \left\{%               \begin{array}{ll}%                       \sin( \omega_{\max} \tau^2 / 2\tau_{\max} )%                               & \quad 0 \le \tau \le \tau_{\max}  \\%                       0       & \quad \mbox{otherwise}%               \end{array}%       \right.%\end{equation}%The phase $\phi$ and instantaneous frequency $\omega$%of this chirp signal is%\begin{eqnarray}%\phi   &=& 0.5\ \omega_{\max} \ \tau^2 / \tau_{\max}   \\%\omega &=& {d\phi \over d\tau} \quad = \quad \omega_{\max} \ \tau / \tau_{\max}%\end{eqnarray}%which shows that as%$\tau$ ranges from $0$ to $\tau_{\max}$,%the frequency $\omega$ ranges from $0$ to $\omega_{\max}$.%\par%Now consider the chirp signal seen on a wave, say $f(t-px)$.%For constant $x$ we see a chirp signal,%possibly shifted on the time axis.%For constant $t$ we see a chirp function%along the $x$ axis,%but the $-p$ causes the chirp to%be reversed and stretched.%Figure~\FIG{alias2d} shows%the chirp wave function $f(t-px)$ for $p=0.45$%and for $\omega_{\max}=2\pi$.%\activeplot{alias2d}{width=6in,height=3in}{}{%       A coarsely sampled chirp function on a wavefront.%       View figure at a grazing angle from various directions.%       }%Since $\omega_{\max}=2\pi$ is double the Nyquist frequency,%you see aliasing on the time axis in Figure~\FIG{alias2d}.%In other words, along the left edge, the frequency first%increases with time until it reaches the Nyquist frequency and then,%as frequency goes beyond Nyquist,%the apparent frequency decreases.%Something else happens on the space axis%which is best seen along the bottom edge of the plot.%There the frequency sweeps steadily from zero to Nyquist.%\par%Thus, we see that the continuum is different from the sampled world.%In a continuum, the function $f(t-px)$%at constant $t$ looks like the same function at constant $x$%except for shifting, reversal, and stretching.%In the sampled world we see different functions on each axis.%The sampled world matches the continuum only at low frequencies.%\par%Real life would be more like Figure~\FIG{alias2d}%if we would interchange the time and space axis%because it is easy to sample densely enough in time,%or to high-cut filter the data before sampling it%but it is hard to get enough observing locations in space.%The result of insufficient density of stations%is that seismic data is often aliased in space.%Spatial aliasing is when a spatial frequency%looks lower than it would%if the data were more densely sampled in space.%%\section{ORTHOGONALITY OF CROSSING PLANE WAVES}%%\par%%As an alternative to the customary approach%%of defining and analyzing aliasing by Fourier analysis,%%I suggest that spatial aliasing may be defined and analyzed%%with reference to plane waves.%%Consider two sinusoidal signals of different frequencies.%%Since they have different frequencies%%they should be orthogonal,%%but they need not be orthogonal if there is aliasing,%%because aliasing can make a high frequency look like a low one.%%Now let us see if we can substitute%%the word ``dip''%%for the word ``frequency.''%%Can we say that plane waves with different dips should be orthogonal?%%%%\par%%Normally, waves do not contain zero frequency.%%Thus the time integral of a waveform normally vanishes.%%Likewise,%%for a dipping plane wave, the time integral vanishes.%%Likewise,%%a line integral%%across the $(t,x)$-plane%%along a straight line that crosses a plane wave%%or a dipping plane wave vanishes.%%Likewise,%%two plane waves with different slopes should be orthogonal%%if one of them has zero mean.%%The other wave, however, need not have zero mean.%%This ``other wave'' need not be a wave at all,%%but could be an impulse function,%%say $\delta(t-px)$,%%spread over a mesh.%%Once it is spread out, it will look like an impulsive plane wave.%%The purpose of this impulsive plane wave is to multiply onto data,%%thereby enabling us to do a line integral of the data.%%What is important is how we represent this impulse function on the mesh.%%Theoretically, we have line integrals that should vanish%%when they cross a zero-mean plane wave.%%When in practice they do not,%%it means that we have not done%%a good job of representing the line integral.%%%%%%\par%%Fourier analysis and sinc functions%%play an important role in integrals along straight lines.%%Production reflection seismic data processing%%involves much effort with line integrals%%along hyperbolic curves and other shapes---rarely straight lines.%%Hopefully, the theoretical concepts above will be suggestive%%of further theoretical and practical developments.%%%%%\subsection{Rambling thoughts}%%%\par%%%Next I assume considerable facility with Fourier analysis.%%%Now imagine the function for Figure \FIG{alias2d} had been %%%a sinc function with its width parameter adjusted so%%%to match aliasing on the coarse axis (axis 1).%%%The FT of a line is another line.%%%I have a sinc convolved with a line.%%%So in $(\omega, k)$-space it seems like a line%%%convolved (multiplied?) with a rectangle.%%%%%%%%%%%\subsection{Relation of missing data to inversion}%We take {\it data space} to be a uniform mesh%on which some values are given and some are missing.%We rarely have missing values on a time axis,%but commonly have missing values on a space axis,%i.e.,~missing signals.%Missing signals (traces) happen occasionally for miscellaneous reasons,%and they happen systematically because of \bx{alias}ing and \bx{truncation}.%The aliasing arises for economic reasons---saving instrumentation%by running receivers far apart.%Truncation arises at the ends of any survey,%which, like any human activity, must be finite.%Beyond the survey lies more hypothetical data.%The traces we will find for the \bx{missing data}%are not as good as real observations,%but they are closer to reality%than supposing unmeasured data is zero valued.%Making an image with a single application of an adjoint modeling operator%amounts to assuming that data vanishes beyond its given locations.%\bxbx{Migration}{migration}%is an example of an economically important process%that makes this assumption.%Dealing with missing data is a step beyond this.%In \bx{inversion}, restoring \bx{missing data}%reduces the need for arbitrary model filtering.%\subsection{The world is filled with local plane waves.}%\sx{plane waves}%In your ears now are sounds from various directions.%From moment to moment the directions change.%Momentarily, a single direction (or two) dominates.%Your ears sample only two points in x-space.%Earthquake data is a little better.%Exploration data is much better and%sometimes seems to satisfy the Nyquist requirement,%especially when we forget that the world is 3-D.\parWe often characterize data fromany region of $(t,x)$-space as ``good'' or ``noisy''when we really mean it contains ``few'' or ``many'' plane-wave eventsin that region.Where regions are noisy,there is no escaping the simple form of the Nyquist limitation.Where regions are good we may escape it.Real data typically contains both kinds of regions.Undersampled data with a broad distribution of plane waves is nearly hopeless.Undersampled data with a sparse distribution of plane wavesoffer us the opportunity to resample without aliasing.Consider data containing a spherical wave.The angular bandwidth in a plane-wave decomposition appears huge{\it until we restrict attention to a small region} of the data.(Actually a spherical wave contains very little informationcompared to an arbitrary wave field.)It can be very helpful in reducing the local angular bandwidthif we can deal effectively with tiny pieces of data.%as we did in chapter \UNI .If we can deal with tiny pieces of data,then we can adapt to rapid spatial and temporal variations.This chapter shows such tiny windows of data.\section{INTERPOLATION BEYOND ALIASING}%I have long marveled at the ability of humans to interpolate%seismic data containing mixtures of dips where spatial frequencies%exceed the Nyquist limits.%These limits are hard limits on migration programs.%Costly field-data-acquisition activities%are designed with these limits in mind.%I feared this human skill of going beyond the limit was deeply nonlinear%and beyond reliable programming.%Now, however, I have obtained results%comparable in quality to those of S.~\bx{Spitz},%and I am doing so in a way that seems reliable---using two-stage,%linear least squares.%\par%Here we attack the problem of missing data arising in a regular pattern%such as alternate traces being missing.%This defeats earlier methods%because every output requires missing data as inputs.%\subsection{Traditional 2-D interpolation before aliasing}A traditional method of data interpolation on a regular meshis a four-step procedure:(1) Set zero values at the points to be interpolated;(2) \bx{Fourier transform};(3) Set to zero the high frequencies;and(4) Inverse transform.This is a fine method and is suitable for many applicationsin both one dimension and higher dimensions.However,this method fails to take advantage of our prior knowledgethat seismic data has abundant fragments of plane wavesthat link an axis that is not aliased (time)to axes that often are (space).%Where the method falls down is where more is needed than simple%interlacing---for example,%when signal values are required beyond the ends of%the data sample.%The simple Fourier method of interlacing also loses its applicability%when known data is irregularly distributed.%An example of an application in two dimensions of the methodology%of this section is given in the section on tomography%beginning on page~\pageref{lal/tomography'}.\subsection{Interlacing a filter}\sx{interlacing a filter}\sx{filter ! interlaced}The filter below canbe designed despite alternate missing traces.This filter destroys plane waves.If the plane wave should happen to pass halfway betweenthe ``d'' and the ``e'',those two points could interpolate the halfway point,at least for well-sampled temporal frequencies,and the time axis should always be well sampled.For example, $d=e=-.5$ would almost destroy the plane waveand it is an aliased planewave for its higher frequencies.\begin{equation}   \begin{array}{ccccccccc}      a     &\cdot &b     &\cdot &c     &\cdot &d     &\cdot &e     \\      \cdot &\cdot &\cdot &\cdot &\cdot &\cdot &\cdot &\cdot &\cdot \\      \cdot &\cdot &\cdot &\cdot &1     &\cdot &\cdot &\cdot &\cdot \end{array}\label{eqn:ilacefil}\end{equation}We coulduse module \texttt{pef} \vpageref{lst:pef}to find the filter (\ref{eqn:ilacefil}),if we set up the %constraint array {\tt afre($\ast$,$\ast$)} lag table \texttt{lag}appropriately.Then we could throw away alternate zeroed rows and columns (rescale the lag) to get the filter\begin{equation}   \begin{array}{ccccc}      a     &b     &c     &d     &e     \\      \cdot &\cdot &1     &\cdot &\cdot \end{array}\label{eqn:noilace}\end{equation}which could be used with subroutine \texttt{mis1()} \vpageref{lst:mis1},to find the interleaved databecause boththe filters (\ref{eqn:ilacefil}) and (\ref{eqn:noilace})have the same dip characteristics.\inputdir{lace}\parFigure~\ref{fig:lace3} shows three plane waves recorded on five channelsand the interpolated data.\plot{lace3}{width=6in}{  Left is five signals, each showing three arrivals.  With the data shown on the left  (and no more),  the signals have been interpolated.  Three new traces appear between each given trace,  as shown on the right.}Both the original data and the interpolated data can be describedas ``beyond \bx{alias}ing,'' because on the input datathe signal shifts exceed the signal duration.The calculation requires only a few secondsof a two-stage least-squares method, in which the first stageestimates a PEF (inverse spectrum) of the known data,and the second uses the PEF to estimate the missing traces.Figure~\ref{fig:lace3} comes from PVIwhich introduces the clever method described above.We will review how that was done and examine the F90 codesthat generalize it to $N$-dimensions.Then we'll go on to more general methodsthat allow missing data in any location.Before the methods of this sectionare applied to field data for migration,data must be broken into many overlapping tilesof size about like those shown hereand the results from each tile pieced together.That is described later in chapter \ref{pch/paper:pch}.\parA PEF is like a differential equation.The more plane-wave solutions you expect,the more lags you need on the data.Returning to Figure~\ref{fig:lace3},the filter must cover four traces (or more)to enable it to predict three plane waves.In this case,\texttt{na=(9,4)}.As usual, the spike on the 2-D PEF is at\texttt{center=(5,1)}.We see the filter is expanded by a factor of\texttt{jump=4}.The data size is\texttt{nd=(75,5)}and \texttt{gap=0}.Before looking at the code\texttt{lace} \vpageref{lst:lace}for estimating the PEF,it might be helpful to recall the basic utilities\texttt{line2cart} and\texttt{cart2line}\vpageref{lst:cartesian}for conversion between a multidimensional space andthe helix filter lag.We need to sweep across the whole filterand ``stretch'' its lags on the 1-axis.We do not need to stretch its lags on the 2-axisbecause the data has not yet been interlaced by zero traces.\moddex{lace}{fill missing traces by rescaling PEF}{31}{69}{user/gee}The line \texttt{ii[0] *= jump}means we interlace the 1-axis but not the 2-axis becausethe data has not yet been interlaced with zero traces.\begin{comment}For a 3-D filter\texttt{aa(na1,na2,na3)},the somewhat obtuse expression\texttt{(/na(1)*jump, na(2:)/)}is a three componentvector containing\texttt{(na1*jump, na2, na3)}.\end{comment}\parAfter the PEF has been found, we can get missing data inthe usual way with with module