% This texfile was generated by executing Latify on file paper_act.tex
\documentclass[12pt]{sep2e}
\usepackage{html}
\begin{document}

\def\figdir{Fig}
\lefthead{Guerra and Biondi}
\righthead{Phase encoding with Gold codes}
\footer{SEP--136}
\title{Phase encoding with Gold codes for wave-equation migration }
\author{Claudio Guerra and Biondo Biondi}

%\maketitle

\begin{abstract}
Prestack exploding-reflector modeling aims to synthesize a small dataset comprised of areal shots, while preserving the correct kinematics to be used in iterations of migration velocity analysis. To achieve this goal, the amount of data is reduced by combining the modeled areal data into sets, we call super-areal data. However, crosstalk arises during migration due to the correlation of wavefields resulting from different modeling experiments. Phase encoding the modeling experiments can attenuate crosstalk during migration. In the geophysical community, the most used phase-encoding schemes are plane-wave-phase encoding and random-phase encoding. Here, we exploit the application of Gold codes commonly used in wireless communication, radar and medical imaging communities to phase encode data. We show that adequately selecting the Gold codes can potentially shift the crosstalk out of the migration domain, or the region of interest if a target-oriented approach is used, yielding an image free of crosstalk.
\end{abstract}

\section{Introduction}
\cite{Biondi.sep.124.biondo1,biondi131} introduced the concept of the prestack exploding-reflector modeling. This method synthesizes source and receiver wavefields along the entire survey at the surface, in the form of areal data, starting from a prestack migrated image cube computed with wave-equation migration. For migration velocity analysis, the aim is to generate a considerably smaller dataset than the one used in the initial migration, while maintaining the necessary kinematic information to update the velocity.
\par
 Conceptually, the synthesized areal data are computed by upward propagating source and receiver wavefields using subsurface-offset-domain common-image gathers (SODCIGs) as initial conditions. To decrease the number of experiments to migrate, we take advantage of the linearity of the wave propagation to combine several experiments into a set of composite records. Combining several experiments gives rise to crosstalk during imaging \cite[] {Biondi.sep.124.biondo1, guerra}. \cite {guerra} use pseudo-random-phase encoding \cite[]{romero} during the modeling step to attenuate crosstalk. 
\par
It is common, in the exploration geophysics community, to employ pseudo-random codes using intrinsic functions specific to the programing language. These pseudo-random codes present, generally, a uniform distribution. Their autocorrelation and cross-correlation functions have no special properties. The autocorrelation function presents nearly periodic side lobes with additive low-amplitude random variations. The peak-to-side lobe ratio is around 30. The cross-correlation function is pseudo-random, and its amplitudes are of the same order of magnitude as those of the non-zero lags of the autocorrelation function. Herein, theses codes are called conventional random codes.
\par
In wireless communication, especially for systems using Code Division Multiple Access (CDMA), a class of different pseudo-random codes have been widely used \cite[]{shi}. These codes are binary sequences and have unique autocorrelation and cross-correlations properties which make them more suited to achieve the above-mentioned objectives with minimal crosstalk. The autocorrelation function is represented by a large peak, whose amplitude equals the number of samples in the code, and the cross-correlation peaks, at non-zero lags, with the same amplitudes as that of the autocorrelation.  Examples of binary pseudo-random codes used by these communities are Golay \cite[]{golay, tseng}, Kasami \cite[] {kasami1} and Gold codes \cite[] {gold}. Medical imaging \cite[] {gran} and radar communities \cite[] {levanon} also exploit the statistical properties of these pseudo-random codes to increase bandwidth, signal-to-noise ratio and pulse compression. 
\par
\cite {quan} analyze orthogonal codes to encode simultaneous source signatures for cross-well surveys, and conclude that m-sequences and Gold codes provide the best results on the separation of the seismograms. Here, we exploit the properties of the Gold codes to encode the prestack exploding reflector modeling experiments.
\par
In the next section we give a brief description of the prestack exploding-reflector modeling. Then we discuss how to compute the Gold codes. To illustrate the effectiveness of phase encoding with Gold codes, we compare the migration of prestack-exploding-reflector modeled data encoded with conventional random codes and Gold codes.

\section{Prestack exploding-reflector modeling}
Starting from a prestack image obtained by wave-equation migration represented by a single SODCIG, areal source and receiver wavefields are modeled at the surface by
\begin{equation}
\begin{array}{cc}
S(x,y,z=0,\omega;{\bf x}_m) = G({\bf x}_m-{\bf h};x,y,z=0,\omega)*I_s({\bf x},{\bf h};{\bf x}_m), \\
R(x,y,z=0,\omega;{\bf x}_m) = G({\bf x}_m+{\bf h};x,y,z=0,\omega)*I_r({\bf x},{\bf h};{\bf x}_m), 
\end{array}
\label{eq:eq1}
\end{equation}
where $S(x,y,z=0,\omega;{\bf x}_m)$ is the source wavefield and $R(x,y,z=0,\omega;{\bf x}_m)$ is the receiver wavefield. $I_s({\bf x},{\bf h};{\bf x}_m)$ and $I_r({\bf x},{\bf h};{\bf x}_m)$ are the prestack images used as initial conditions for the source and receiver wavefield extrapolation, respectively, at a selected position, ${\bf x}_m$. These prestack images should be dip-independent gathers. They are computed by re-mapping the dip along the offset direction according to the apparent geological dip \cite[]{biondi131}. $G({\bf x}_m\pm{\bf h};x,z=0,\omega)$ represents the operator that extrapolates the wavefields from the subsurface to the surface; ${\bf h}$ is the subsurface offset; $\omega$ is the temporal frequency; and ${\bf x}$ is the vector of spatial coordinates.
\par
In the case of using an one-way extrapolator, the source and receiver wavefields are upward continued according to the one-way wave equations
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 {s}^2({\bf x})-|{\bf k}| ^2}\right) S({\bf x},\omega;{\bf x}_m) = I_s({\bf x},{\bf h};{\bf x}_m) \\
S(x,y,z=z_{\rm max},\omega;{\bf x}_m) = 0 \end{array} \right.,
\label{eq:eq2}
\end{eqnarray}
and
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}-i\sqrt{\omega ^2 {s}^2({\bf x})-|{\bf k}| ^2}\right) R({\bf x},\omega;{\bf x}_m) = I_r({\bf x},{\bf h};{\bf x}_m) \\
R(x,y,z=z_{\rm max},\omega;{\bf x}_m) = 0 \end{array} \right.,
\label{eq:eq3}
\end{eqnarray}
where $s({\bf x})$ is the slowness at ${\bf x}$; ${\bf k}=(k_x,k_y)$ is the spatial wavenumber vector.
\par
Using the linearity of the wave propagation, sets of individual modeling experiments can be combined into the same areal data, such that the amount of data input into migration can be significantly decreased, reducing its cost. However, this procedure generates crosstalk when applying the imaging condition during migration.
\par
\cite{guerra} introduce strategies to attenuate the crosstalk. Migration of $({\bf x},\omega)$-random-phase encoded data disperses the crosstalk energy throughout the image as a pseudo-random background noise. By adding more realizations of random-phase encoded areal data, the speckled noise can be further attenuated. The encoded source wavefield, ${\widetilde S}({\bf x},{\bf p}_m,\omega)$, and the encoded receiver wavefield, ${\widetilde R}({\bf x},{\bf p}_m,\omega)$, are synthesized according to
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 {s}^2({\bf x})-|{\bf k}| ^2}\right) {\widetilde S}({\bf x},{\bf p}_m,\omega) 
= \widetilde{I}_s({\bf x},{\bf h},{\bf p}_m,\omega) \\
{\widetilde S}(x,y,z=z_{\rm max},{\bf p}_m,\omega) = 0 \end{array} \right.,
\label{eq:eq4}
\end{eqnarray}
and
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}-i\sqrt{\omega ^2 {s}^2({\bf x})-|{\bf k}| ^2}\right) {\widetilde R}({\bf x},{\bf p}_m,\omega) 
= \widetilde{I}_r({\bf x},{\bf h},{\bf p}_m,\omega) \\
{\widetilde R}(x,y,z=z_{\rm max},{\bf p}_m,\omega) = 0 \end{array} \right.,
\label{eq:eq5}
\end{eqnarray}
where $\widetilde{I}_s({\bf x},{\bf h},{\bf p}_m,\omega)$ and $\widetilde{I}_r({\bf x},{\bf h},{\bf p}_m,\omega)$ 
are the encoded SODCIGs after rotations. They are defined as follows:
\begin{eqnarray}
\widetilde{I}_s({\bf x},{\bf h},{\bf p}_m,\omega) &=& \sum_{{\bf x}_m}I_s({\bf x},{\bf h},{\bf x}_m)\beta({\bf x},{\bf x}_m,{\bf p}_m,\omega), \\
\widetilde{I}_r({\bf x},{\bf h},{\bf p}_m,\omega) &=& \sum_{{\bf x}_m}I_r({\bf x},{\bf h},{\bf x}_m)\beta({\bf x},{\bf x}_m,{\bf p}_m,\omega),
\end{eqnarray}
where $\beta({\bf x},{\bf x}_m,{\bf p}_m,\omega)=e^{i\gamma({\bf x},{\bf x}_m,{\bf p}_m,\omega)}$ is the phase-encoding function; the variable ${\bf p}_m$ is the index of different realizations of phase encoding. 
\par
The areal shot migration is performed by downward continuation of the areal source and receiver wavefields according to the following one-way wave equations
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}-i\sqrt{\omega ^2 {s}^2({\bf x})-|{\bf k}| ^2}\right) {\widehat S}({\bf x},{\bf p}_m,\omega) 
= 0 \\
{\widehat S}(x,y,z=0,{\bf p}_m,\omega) = {\widetilde S}(x,y,z=0,{\bf p}_m,\omega) \end{array} \right.,
\label{eq:eq6}
\end{eqnarray}
and
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 {s}^2({\bf x})-|{\bf k}| ^2}\right) {\widehat R}({\bf x},{\bf p}_m,\omega) 
= 0 \\
{\widehat R}(x,y,z=0,{\bf p}_m,\omega) = {\widetilde R}(x,y,z=0,{\bf p}_m,\omega) \end{array} \right.,
\label{eq:eq7}
\end{eqnarray}
where the encoded source wavefield, ${\widetilde S}({\bf x},{\bf p}_m,\omega)$, and the encoded receiver wavefield, ${\widetilde R}({\bf x},{\bf p}_m,\omega)$, are used as boundary conditions.
\par
The image, $\widehat{I}({\bf x},{\bf h})$, is obtained by cross-correlation of the source wavefield, ${\widehat S}({\bf x},{\bf p}_m,\omega)$, with the receiver wavefield, ${\widehat R}({\bf x},{\bf p}_m,\omega)$
\begin{equation}
\widehat{I}({\bf x},{\bf h}) = \sum_{\omega} \sum_{{\bf p}_m} \widehat{S}^{\star}({\bf x}-{\bf h},{\bf p}_m,\omega) 
\widehat{R}({\bf x}+{\bf h},{\bf p}_m,\omega),
\label{eq:eq8}
\end{equation}
where $^{\star}$ represents complex conjugation.

\section{Gold codes}
Before describing Gold codes it is useful to define maximum length sequences.
\par
Linear feedback shift registers (LFSR) \cite[]{moon} are called state machines, whose components and functions are:
\begin{itemize}
\item the {\bf shift register} -- shifts the bit pattern and registers the output bit; and
\item the {\bf feedback function} -- computes the input bit according to the tap sequence and inserts the computed bit into the input bit position.
\end{itemize}
The output sequence of bits form pseudo-random binary sequences, which are completely controlled by the tap sequence. A tap sequence defines which bits in the current state will be combined to determine the input bit for the next state. The combination is generally performed using module-2 addition ($\it exclusive~or$ -- XOR). This means that adding the selected bit values defined by the tap sequence, if the sum is odd the output of the function is one; otherwise the output is zero. Table~\ref{tab:lfsr} shows the internal states and the output sequence of a 4-bit LFSR with tap sequence $[4,1]$. For the current state, the input bit (bit 1) is computed by the sum module--2 of the bits defined by the tap sequence(bits 1 and 4). The rest of the bits in the register (bits 2, 3 and 4) are obtained by shifting the bit values in the previous state to the right. For example, bit 2 of the current state is bit 1 of the previous state. The output bit of the current state is the last bit (bit 4) in the register from the previous state. 
\par
  \begin{table}[h]
   \begin{center}
    \begin{tabular}{|c|c|c|c|c|}
    \hline
    \multicolumn{4}{|c|}{{\bf Register States}} & \\ \hline
    bit1 & bit2 & bit3 & bit4 & {\bf Output} \\
    (Tap) &  &  & (Tap) & {\bf Sequence} \\ \hline
    1&1&0&1&   \\ \hline
    0&1&1&0&1  \\ \hline
    0&0&1&1&0  \\ \hline
    1&0&0&1&1  \\ \hline
    0&1&0&0&1  \\ \hline
    0&0&1&0&0  \\ \hline
    0&0&0&1&0  \\ \hline
    1&0&0&0&1  \\ \hline
    1&1&0&0&0  \\ \hline
    1&1&1&0&0  \\ \hline
    1&1&1&1&0  \\ \hline
    0&1&1&1&1  \\ \hline
    1&0&1&1&1  \\ \hline
    0&1&0&1&1  \\ \hline
    1&0&1&0&1  \\ \hline
    1&1&0&1&0  \\ \hline
    \end{tabular}
    \caption{Internal states and the output sequence of a 4-bit LFSR with tap sequence [4,1].}
    \label{tab:lfsr}
   \end{center}
  \end{table}
\par

 In number theory, Galois fields (${\bf GF}$) are finite fields in which all operations result in an element of the field \cite[]{lidl}. Addition, subtraction and multiplication of polynomials are defined in a finite field. Module 2 arithmetic forms the basis of ${\bf GF}(2)$ (Galois field of order 2). Addition and multiplication operations in GF(2) can be represented by bitwise operators XOR and AND, respectively. Table~\ref{tab:gf2} synthesizes the possible output values of addition and multiplication over ${\bf GF}(2)$.

  \begin{table}[h]
   \begin{center}
    \begin{tabular}{|c|c|c| c |c|c|c|}
    \hline
    + & 0 & 1 &  & $\times$ & 0 & 1 \\ \hline
    0 & 0 & 1 &  & 0 & 0 & 0 \\ \hline
    1 & 1 & 0 &  & 1 & 0 & 1 \\ \hline
    \end{tabular}
    \caption{${\bf GF}(2)$ addition and multiplication possible outcomes.}
    \label{tab:gf2}
   \end{center}
  \end{table}
\par

 If a polynomial can not be represented as the product of two or more polynomials, it is called an irreducible polynomial. For instance, $x^2+x+1$ is irreducible over ${\bf GF}(2)$ because it can not be factored. However, $x^2+1$ is not irreducible over ${\bf GF}(2)$ because, using normal algebra, $(x+1)(x+1)=x^2+2x+1$, and after reduction module--2, is $x^2+1$ (the term $2x$ is dropped). So, in ${\bf GF}(2)$, $x^2+1 \equiv x^2+2x+1$.
\par
The importance of studying irreducible polynomials ${\bf GF}(2)$ is that they are used to represent tap sequences. Considering an irreducible polynomial, the corresponding tap sequence is given by the exponents of the terms with coefficients of 1 \cite[] {dinan}. Special tap sequences can be used to generate particular  pseudo-random binary sequences. They are called maximum length sequences (m-sequences) and, by definition, are the largest codes that can be generated by a LFSR for a given tap sequence. Their length is $(b^n-1)$, where $n$ is the number of elements of the tap sequence, and $b=2,3~or~5$.
\par
The autocorrelation function of an m-sequence, $\Phi_{\it mls}(k)$, is given by 
\begin{equation}
\Phi_{\it mls}(k) = \left\{ \begin{array}{r}
b^n-1 \qquad\hbox{for}\quad k = 0, \cr
-1 \qquad\hbox{for}\quad k \neq 0,
\end{array} \right.
\label{eq:mseq}
\end{equation}
where $k$ is the lag of correlation.
In spite of the good autocorrelation properties, m-sequences, in general, are not immune to cross-correlation problems, and they may have large and unpredictable cross-correlation values. However, the so-called preferred pairs of m-sequences have cross-correlation functions which might assume the predicted values, $-1,~-1+p$, and $-1-p$, where $p=2^{(n+1)/2}$ for $n$ odd or $p=2^{(n+2)/2}$ for $n$ even. Given a $(2^{n}-1)$--length m-sequence, $a(k)$, and $gcd\{n,4\} = 1$ (greatest common divisor of $n$ and $4$), its preferred pair is the result of decimation computed by applying on $a(k)$ a circular shift of $q$ samples, where $q=2m+1$ and $gcd\{m,n\}=1$. Figure~\ref{fig:prefpairs} shows the cross-correlation of m-sequences that form a preferred pair computed with $m=5$. 
\activeplot{prefpairs}{width=3in}{CR}{Cross-correlation of a m-sequences that form a preferred pair.$\bf{[ER]}$}

The number of possible preferred pairs of m-sequences is limited, when compared to the requirements of practical applications of wireless communication. Preferred pairs of m-sequences, however, are used to generate Gold codes \cite[]{dinan}. 
\par
In CDMA, Gold codes are used as chipping sequences that allow several callers to use the same frequency, resulting in less interference and better utilization of the available bandwidth. Originally proposed by \cite {gold}, Gold codes can be computed by module-2 addition ($\it exclusive~or$) of circularly shifted preferred pairs of m-sequences of length $2^{n}-1$. The autocorrelation function of a Gold code, $\Phi_{\it gc}(k)$, is given by
\begin{equation}
\Phi_{\it gc}(k) = \left\{ \begin{array}{r}
\pm 2^{n}-1 \qquad\hbox{for}\quad k = 0, \cr
\pm 1 \qquad\hbox{for}\quad k \neq 0.
\end{array} \right.
\label{eq:gseq}
\end{equation}
More interestingly, the two valued cross-correlation function of Gold sequences, $\Psi_{\it gc}(k)$, is given by
\begin{equation}
\Psi_{\it gc}(k) = \left\{ \begin{array}{r}
\pm (2^{n}-1) \qquad\hbox{for}\quad k = \lambda, \cr
\mp 1 \qquad\hbox{for}\quad k \neq \lambda.
\end{array} \right.
\label{eq:xgold}
\end{equation}
where the correlation lag $\lambda$ is given by the difference between the number of circular shifts applied to the m-sequence to compute the Gold codes.
\par
Figure \ref{fig:gold184} illustrates the correlation properties of the Gold codes. The left part shows the autocorrelation of the Gold code generated with one circular shift of the preferred pair of m-sequence. The right part shows the cross-correlation of the Gold code generated with one circular shift of the preferred pair of m-sequence with the Gold code generated with 84 circular shifts. Note that the peak of the cross-correlation occurs at lag 84. For comparison, we show the correlation functions of conventional random codes in Figure~\ref{fig:rand}. The autocorrelation function presents nearly periodic side lobes with additive low-amplitude random variations. The peak-to-side lobe ratio is around 30. The cross-correlation function is pseudo-random, and its amplitudes are of the same order of magnitude as those of the non-zero lags of the autocorrelation function. 
\activeplot{gold184}{width=\textwidth}{CR}{Correlation functions of Gold codes. a) Autocorrelation of Gold code generated with one circular shift of the preferred pair of m-sequence. b) Cross-correlation of the Gold code generated with one circular shift of the preferred pair of m-sequence with that generated with 84 circular shifts. The peak of the cross-correlation occurs at lag -84.$\bf{[ER]}$}
\activeplot{rand}{width=\textwidth}{CR}{Correlation functions of conventional random codes. a) Autocorrelation of the conventional random code. b) A cross-correlation.$\bf{[ER]}$}

After computing Gold codes, we use their phase information to encode the modeling experiments.

\section{Examples}
We illustrate the use of the  encoding methods on a simple model of a flat reflector, 0.5 km deep, embedded in a medium with a constant velocity of 2 km/s. The original data is migrated with a 5$\%$ slower velocity. We used the same slower velocity to perform the modeling. For the areal shot migration, we show examples of migrating Gold code phase-encoded data with the slower velocity. For comparison of the crosstalk behavior, when migrating with a velocity different from that used to model, we also show images migrated  with a 5$\%$ faster velocity. Super-areal data are comprised of the collection of 10 modeling experiments initiated at every 10th CMP coordinate. We computed the prestack migrated image with 61 subsurface offsets to observe how crosstalk is shifted when using Gold codes. It should have been reasonable to use a much smaller number of offsets, as the moveout information is  restricted to the 21 central offsets.
\par 
Figure \ref{fig:perm0} shows the areal shot migration of data generated by the prestack exploding-reflector modeling without combining the modeling experiments into super-areal shots. The panel on the left is the zero-subsurface-offset section, and the panel on the right is a SODCIG. This result represents the ideal image we would like to  obtain if the crosstalk could be eliminated. Our objective in phase encoding the modeling experiments is to achieve satisfactory crosstalk attenuation in such a way that the moveout information is not altered. 
\activeplot{perm0}{width=4in}{CR}{Areal shot migration of synthesized data with no combination of the modeling experiments into super-areal data.$\bf{[CR]}$}

Figure \ref{fig:perma} shows the areal shot migration of data generated by the prestack exploding-reflector modeling with no phase encoding applied. The super-areal data, input to areal shot migration, are comprised of modeling experiments initiated at every tenth SODCIG. The SODCIGs resulting from the areal shot migration, show strong crosstalk at subsurface offsets different from zero. The crosstalk is periodic with a period half of the spacing of the modeling experiments in a super-areal shot. 
\activeplot{perma}{width=4.0in}{CR}{Areal shot migration of synthesized data with no phase encoding applied. The super-areal data comprises 10 modeling experiments. Notice the crosstalk in the SODCIG.$\bf{[CR]}$}

Figure~\ref{fig:conv1r} shows areal shot migration of one realization of phase encoding modeling with conventional random codes. The strong crosstalk observed in Figure~\ref{fig:perm0} is now dispersed throughout the image. The dispersed crosstalk can be further attenuated by migrating more random realizations, but this increases the cost of migration.
\activeplot{conv1r}{width=4in}{CR}{Areal shot migration of one realization of synthesized data with conventional random-phase encoding.$\bf{[CR]}$}
Figure~\ref{fig:conv5r} shows the migration of 4 realizations of conventional random encoding modeling.  
\activeplot{conv5r}{width=4in}{CR}{Areal shot migration of four realizations of synthesized data with conventional random-phase encoding.$\bf{[CR]}$}

Crosstalk attenuation is incomplete when using conventional random codes because their autocorrelations are not a perfect spike, nor are their cross-correlations zero. Gold codes partially satisfy these requirements: the autocorrelation is almost a perfect spike, except for -1's at lags different from zero, and, similarly, the cross-correlations are -1 everywhere except where they peak. Therefore, to obtain good results when using Gold codes, it is critical to select the codes which provide the best crosstalk attenuation. That is because the  cross-correlation functions have peaks with the same magnitude as those of the autocorrelation function. 
\par
In Figure~\ref{fig:gold1x}, the areal shot migration was performed on encoded data with Gold codes that have cross-correlation peaks at every 5th (Figure~\ref{fig:gold1x}a), 10th (Figure~\ref{fig:gold1x}b) and, 20th lags  (Figure~\ref{fig:gold1x}c). This means that when applying the multi-offset imaging condition, unrelated wavefields, delayed in time by the phase functions, will cross-correlate at depths different from that of the related wavefields, which were encoded with the same phase. The crosstalk of the Figure~\ref{fig:perma} is now shifted in depth according to the selected set of Gold codes. In Figure~\ref{fig:gold1x}a, the apexes of the crosstalk are displaced in depth at a constant spacing of, approximately, 0.1 km; in Figure~\ref{fig:gold1x}b the spacing is approximately 0.2 km, and in Figure~\ref{fig:gold1x}c, 0.4 km. Notice that in spite of the crosstalk is still present in Figure~\ref{fig:gold1x}c, its complete elimination can be achieved with a simple depth-windowing of the image.
\activesideplot{gold1x}{width=4in}{CR}{Areal shot migration of one realization of synthesized data with Gold code phase encoding. Gold codes have cross-correlation peaks at every 5th (a), 10th (b) and, 20th lag (c), and the depth shifts are, respectively, 0.1 km, 0.2 km, and 0.4 km.$\bf{[CR]}$}

The amount of depth shift of the crosstalk is defined by the lag where the cross-correlation of the Gold codes peaks. For the simple case of constant velocity, the depth shifts, $\delta z$, are given by
\begin{equation}
\delta z = \frac{n_{\lambda}v}{2n_{\omega}d_{\omega}}
\label{eq:zshift}
\end{equation}
where $n_{\omega}$ is the number of frequencies, $d_{\omega}$ is the frequency interval and $n_{\lambda}$ is the lag where the cross-correlation of the Gold codes peaks. For the present example, this amounts to $\delta z = 0.021*n_{\lambda}$~km.
\par
 One possibility to statistically attenuate the crosstalk is to randomly select the Gold codes. Figure~\ref{fig:gold1r} shows the areal shot migration of one realization of encoded data with randomly selected Gold codes. The crosstalk shows different patterns than that of the sequentially selected Gold codes.
\activeplot{gold1r}{width=4in}{CR}{Areal shot migration of one realization of synthesized data with Gold code phase encoding. Gold codes are randomly selected.$\bf{[CR]}$} 

As before, migrating more realizations of randomly selected Gold-encoded data further attenuates crosstalk. Figure~\ref{fig:gold5r} shows the migration of 4 realizations of randomly selected Gold-encoded data. Comparison with Figure~\ref{fig:gold1r} shows that much of the remaining crosstalk energy has been attenuated. However, this approach does not exploit the statistical properties of the Gold phase functions and the crosstalk shows up locally coherent, at random positions. This strategy is definitely not suited to provide a kinematically reliable image.
\activeplot{gold5r}{width=4in}{CR}{Areal shot migration of four realizations of synthesized data with Gold code phase encoding. Gold codes are randomly selected.$\bf{[CR]}$}

Considering that the crosstalk is shifted in depth, as observed in Figure~\ref{fig:gold1x}, and that the amount of shift is determined by the lag where the cross-correlation of the Gold codes peaks (equation~\ref{eq:zshift}), one can choose the Gold codes according to a suitable interval that completely shifts the crosstalk away from the zone of interest. This strategy shares similar idea as the linear-phase encoding \cite{romero}, which aims to shift the crosstalk out of the migration domain by using a linear function of frequency. In the Appendix, we show that phase encoding with Gold codes is equivalent to linear phase encoding. Figure~\ref{fig:gold1o} shows the areal shot migration of data encoded by selecting every 50th Gold code, meaning that the depth shifts are multiples of 1 km, as predicted by equation~\ref{eq:zshift}, which should be adequate to completely push the crosstalk away from the SODCIG. Contrary to what was expected, the crosstalk is still present. Of course, its complete elimination can be achieved by windowing around the central traces (Figure~\ref{fig:goldwd}).
\activeplot{gold1o}{width=4in}{CR}{Areal shot migration of one realization of synthesized data with Gold phase encoding. Gold codes are selected such that the crosstalk is shifted away from the reflector. A simple windowing should completely eliminate the crosstalk.$\bf{[CR]}$}
\activeplot{goldwd}{width=4in}{CR}{Crosstalk of Figure~\ref{fig:gold1o} can be completely eliminated by windowing the SODCIG. Compare with Figure \ref{fig:perma}.$\bf{[CR]}$}

To understand the origin of the crosstalk in Figure~\ref{fig:gold1o}, we migrated the data without adding the results of individual migrations of the super areal data. Figure~\ref{fig:goldsp} shows three of these migrations. The three panels show the same spatial position. The top panel shows the migration of super areal data whose  modeling includes a SODCIG coincident with that displayed. Remember that we are using a ``comb'' function to select every 10th SODCIG to initiate the modeling experiments comprised by a super areal data. The central panel shows the migration of the super areal data modeled from the initial image shifted five CMP positions away from that of the uppermost panel, and the panel on the bottom shows the migration of the super areal data modeled from the initial image shifted nine CMP positions away from that of the uppermost panel. The modeling experiments which generated the super areal data for Figure~\ref{fig:goldsp}c are actually separated one CMP position from the ones for Figure~\ref{fig:goldsp}a, given that the modeling experiments in super areal data are separated every 10th CMP. Figure~\ref{fig:goldsp}a and Figure~\ref{fig:goldsp}c show that much of the reflector energy comes from the migration of modeling experiments initiated at CMP positions close to the considered image point. Figure~\ref{fig:goldsp}b shows that much of the crosstalk energy of Figure~\ref{fig:gold1o} is related to  migration of super areal shots whose modeling experiments are initiated at SODCIGs shifted five CMP positions with respect to the considered SODCIG. 
\activesideplot{goldsp}{width=4in}{CR}{Areal shot migration of three super areal data. The results were not stacked together. Notice that much of the crosstalk energy of Figure~\ref{fig:gold1o} is related to the migration of super areal shots ($b$) other than the one which modeling was initiated at that CMP position ($a$ and $c$).$\bf{[CR]}$}

%We migrated super areal data which modeling was initiated at CMP positions like the one in Figure~\ref{fig:goldsp}b, using Gold codes with cross-correlation peaks separated at every 5th (Figure~\ref{fig:gold3x}a), 10th (Figure~\ref{fig:gold3x}b), 20th (Figure~\ref{fig:gold3x}c), 30th (Figure~\ref{fig:gold3x}d), 40th (Figure~\ref{fig:gold3x}e), 45th (Figure~\ref{fig:gold3x}f), 47th (Figure~\ref{fig:gold3x}g), and 50th lags (Figure~\ref{fig:gold3x}h). Only the SODCIG at the same position as the one showed in Figure~\ref{fig:goldsp}b is displayed. Figure~\ref{fig:gold3x} shows that the crosstalk for positive subsurface offsets are shifted to greater depths while the crosstalk at negative subsurface offsets are shifted to shallower depths. In addition, it starts wrapping around in a complex manner to finally assume the configuration of the crosstalk in Figure~\ref{fig:gold1o} when the Gold codes with cross-correlation peaks at every 50th lags. The crosstalk in Figure~\ref{fig:gold3x}h is the cross-correlation of source and receiver wavefields encoded with Gold codes built by circularly shifting the preferred pair of m-sequences 250 positions to the left and 250 positions to the right.
%\sideplot{gold3x}{width=4in}{Areal shot migration of of data encoded with Gold codes with cross-correlation peaks at every: a) 5th  lags; b) 10th lags; c) 20th lags; d) 30th lags; e) 40th lags; f) 45th lags; g) 47th lags; and g) 50th lags.}

In migration velocity analysis, at every velocity iteration, data is migrated with an updated velocity. To verify how the crosstalk behaves when migrating with a different velocity, we used a 5$\%$ faster velocity to migrate the Gold encoded data modeled with the initial 5$\%$ slower velocity. We used sets of Gold codes with cross-correlation peaks at every 5th (Figure~\ref{fig:gold2x}a), 10th (Figure~\ref{fig:gold2x}b) and, 20th lag  (Figure~\ref{fig:gold2x}c) lags. The depth shift of the crosstalk is proportional to the increase in velocity. The behavior is similar to that of migration with a slower velocity, regarding that the interference with the region of interest decreases with the increase of the lag of the cross-correlation peaks.
\activesideplot{gold2x}{width=4in}{CR}{Areal shot migration of one realization of synthesized data with Gold phase encoding. Gold codes have cross-correlation peaks at every 5th (a), 10th (b) and, 20th lag (c).$\bf{[CR]}$}
 
The use  of Gold codes can be less costly than conventional random codes, given that just one realization is necessary to achieve an almost perfect result, while even using more realizations of conventional random encoding does not produce an image with similar quality. In addition, it is suited to be used in a horizon based strategy to migration velocity analysis, where a few important reflectors are chosen to do velocity update.
\par
However, an inadequate choice of the Gold codes is potentially more dangerous than using conventional random codes. Because the crosstalk may not be shifted out of the region of interest, coherent artifacts can coincide with reflectors and obliterate the moveout information. Conventional random codes, in turn, randomizes the crosstalk. This can lead to a noisy residual moveout scan or a noisy gradient. In a companion paper \cite{tang136} show that the prestack exploding-reflector modeled data, phase encoded with conventional random codes, is able to provide a reasonable direction for velocity update.

\section{Conclusion}
We showed that the use of Gold codes in phase encoding can virtually eliminate crosstalk if the codes are satisfactorily selected. Its advantage over conventional random codes is not only in the image quality but also cost. The elimination of the crosstalk can be achieved with only one realization of Gold codes, while for conventional random codes more realizations are needed to obtain reasonable images, but still of lower quality when compared to the ones obtained by using Gold codes. The need for more realizations increases the cost of modeling and migration. 
\par
We also showed that the phase of the cross-correlation of Gold codes is given by the lag in time where the cross-correlation peaks multiplied by frequency. This amounts to apply linear phase encoding to the encoded wavefields.
\par
We used a simple example with just one reflector to illustrate the encoding techniques. The performance of the Gold codes in the presence of more than one reflector needs to be tested. However, as just one realization is sufficient to drastically reduce the crosstalk, they can be affordable to use in a horizon-based approach to the modeling in which each reflector is modeled separately. Additional research is needed to analyze the performance of Gold codes in more complex situations, where the depth shift of the crosstalk may not be a simple function of velocity.

\activeplot{9002Dslowgold}{width=4in}{ER}{}
\activeplot{9006Dslowencz4r}{width=4in}{ER}{}

\section{Acknowledgements}
We thank Jerry Harris and Youli Quan for the suggestion of using codes applied in wireless communications.
\bibliographystyle{sep}
\bibliography{SEP,ref134a}
\newpage
\renewcommand{\theequation}{A-\arabic{equation}} % redefine the command that creates the equation no.
\setcounter{equation}{0}  % reset counter
\section*{APPENDIX}  % use *-form to suppress numbering
We show that phase encoding using Gold codes is equivalent to the linear phase encoding introduced by \cite[]{romero}. The cross-correlation function of Gold codes is given by equation~\ref{eq:xgold}. It corresponds to a spike of amplitude $\pm 2^n$ at lag $\lambda$, which depends on the difference between the number of circular shifts applied to the m-sequences to compute the Gold codes, plus a DC term, $\mp 1$. The phase of the cross-correlation function is given by the phase of the spike, $t_{\lambda} \omega$, and is equal to the phase difference of the input signals. If Gold codes have phases $\gamma_1(\omega)$ and $\gamma_2(\omega)$, the phase of their cross-correlation is 
\begin{equation}
\gamma_1(\omega) - \gamma_2(\omega) = t_{\lambda} \omega.
\label{eq:linph}
\end{equation}
According to equation \ref{eq:linph}, the phase of the cross-correlation of Gold codes is a linear function of the frequency.  Equation \ref{eq:linph} is equal to equation 26 of Romero's paper.
\par

\clearpage

\end{document}
