\righthead{Guerra and Biondi}
\lefthead{Pseudo-random-phase encoding}
\footer{SEP--136}
\title{New phase-encoding schemes 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 modeled areal data must be
combined into sets. However, crosstalk arises during migration due to the correlation of unrelated wavefields. 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 also introduce a method to compute random-phase functions by solving an optimization problem ... (more still to come). We show that adequately selecting the Gold codes can potentially shift the crosstalk out of the migration domain, 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
\par
 Conceptually, the synthesized areal data are computed by upward propagating source and receiver wavefields, using as initial conditions subsurface-offset-domain common-image gathers (SODCIG). 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, though, 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. For conciseness reason, from now on we use the term random instead of pseudo-random.
\par
The exploration geophysics community employs pseudo-random codes using intrinsic functions specific to the programing language. These pseudo-random codes present, generally, a uniform distribution and many of the algorithms are based on the ${\it multiplicative}$ ${\it congruential}$ ${\it algorithm}$, invented by \cite{lehmer}. 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, specially for systems using Code Division Multiple Access (CDMA), different pseudo-random codes have been widely used \cite[]{shi}. 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. These codes are binary sequences and have unique autocorrelation and cross-correlations properties which make them more suited to achieving 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[] {kasami} and Gold codes \cite[] {gold}. Here, we exploit the properties of the Gold codes to encode seismic data.
\par
In addition to using Gold codes, we introduce a novel method to compute phase functions by solving an optimization problem ...
\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 and the optimized random-phase functions. We illustrate the application of each phase-encoding scheme by migrating phase-encoded areal data synthesized by prestack exploding-reflector modeling.


\section{Prestack exploding-reflector modeling}
Starting from a prestack image obtained by wave-equation migration, areal source and receiver wavefields are modeled at the surface by
\begin{equation}
\begin{array}{cc}
S(x,\omega) = G(z_{\xi},x_{\xi}-h_{\xi};x,z=0,\omega)*I_s(z_{\xi},x_{\xi},h_{\xi}), \\
R(x,\omega) = G(z_{\xi},x_{\xi}+h_{\xi};x,z=0,\omega)*I_r(z_{\xi},x_{\xi},h_{\xi}), 
\end{array}
\label{eq:eq1}
\end{equation}
where $S(x,\omega)$ is the source wavefield at $z=0$; $R(x,\omega)$ is the receiver wavefield at $z=0$; $I_s(z_{\xi},x_{\xi},h_{\xi})$ and $I_r(z_{\xi},x_{\xi},h_{\xi})$ are the prestack images used as initial conditions for the source and receiver wavefield extrapolation, respectively, at a selected position, $x_{\xi}$,; $G(z_{\xi},x_{\xi} \pm h_{\xi};x,z=0,\omega)$ represents the operator that extrapolates the wavefields from the subsurface to the surface; $h_{\xi}$ is the subsurface offset; $z_{\xi}$ is depth; $\omega$ is the temporal frequency and $x$ is the spatial coordinate in the data space coinciding with $x_{\xi}$. The prestack images used as initial conditions for the source and receiver wavefield extrapolation should be dip-independent gathers computed by changing the dip along the offset direction according to the apparent geological dip \cite[]{biondi131}. Here we use one-way extrapolators for both modeling and migration of the areal data.
\par
By combining sets of individual modeling experiments into the same areal data, the amount of data input into migration can be significantly decreased. We achieve this by regularly selecting individual experiments and adding them up into their set, after being upward propagated, according to
\begin{eqnarray}
\widetilde{S_n}(x,\omega)& = &\sum_{n=1}^k \sum_{i=n,k,N} S_i(x,\omega), \nonumber \\
\widetilde{R_n}(x,\omega)& = &\sum_{n=1}^k \sum_{i=n,k,N} R_i(x,\omega), 
\label{eq:eq2}
\end{eqnarray}
where $\widetilde{S_n}(x,\omega)$, and  $\widetilde{R_n}(x,\omega)$ contain $k$ sets of summed areal sources and areal receivers, respectively; $N$ is the number of SODCIGs. Every $k$th areal datum is selected to compose one set. Pairs of $\widetilde{S_k}(x,\omega)$ and  $\widetilde{R_k}(x,\omega)$ are to be used as the areal source and the areal receiver wavefields, respectively, in areal shot migration.
\par
Migration of the combined areal data produces an image by cross-correlating the combined areal source and receiver wavefields,
\begin{equation}
\widetilde{I_m}(z_{\xi},x_{\xi},h_{\xi}) = \sum_{\omega} \widetilde{S_m}^{\star}(z_{\xi},x_{\xi}-h_{\xi},\omega) 
\widetilde{R_m}(z_{\xi},x_{\xi}+h_{\xi},\omega),
\label{eq:eq3}
\end{equation}
where $^{\star}$ represents complex conjugation. If $\widetilde{S_m}(x,\omega)$ and $\widetilde{R_m}(x,\omega)$ in equation~\ref{eq:eq2} are comprised by two summed areal shots, the image $\widetilde{I_m}(z_{\xi},x_{\xi},h_{\xi}) $ will be given by
\begin{eqnarray}
\lefteqn{\widetilde{I_m}(z_{\xi},x_{\xi},h_{\xi}) = I_1(z_{\xi},x_{\xi},h_{\xi}) + I_2(z_{\xi},x_{\xi},h_{\xi}) + } \nonumber \\
 &  &  \sum_{\omega}S_1^{\star}(z_{\xi},x_{\xi}-h_{\xi},\omega)R_2(z_{\xi},x_{\xi}+h_{\xi},\omega) + \nonumber \\
 &  &  \sum_{\omega}S_2^{\star}(z_{\xi},x_{\xi}-h_{\xi},\omega)R_1(z_{\xi},x_{\xi}+h_{\xi},\omega).
\label{eq:eq4}
\end{eqnarray}
In equation \ref{eq:eq4}, the last two summation terms represent the crosstalk. 
\par
\cite{guerra} show that in the prestack exploding-reflector method the crosstalk has two distinct origins. The crosstalk that is prominent in the zero-subsurface offset section, which herein we call type-1, is due to the cross-correlation of events in the source wavefield with that in the receiver wavefields modeled by the same SODCIG. Crosstalk in non-zero subsurface offsets, which herein we call type-2, is related to cross-correlation of events in the source wavefields with those in the receiver wavefields modeled by different SODCIGs. 
\par
\cite{guerra} introduce strategies to attenuate the crosstalk. Migration of $(x,z,\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 is further attenuated.

\section{Optimized random codes}
...

\section{Gold codes}
Before describing Gold codes it is useful to define maximum length sequences.
\par
Linear feedback shift registers (LFSR) are called state machine,s whose main operations are:to do the following

\begin{itemize}
\item compute the input bit according to the tap sequence;
\item shift the bit pattern;
\item register the output bit; and
\item insert the bit computed in the first step into the input bit position.  
\end{itemize}
Maximum length sequences (m-sequences) are composed by the output bits of a LFSR. They are, by definition, the largest codes that can be generated by a LFSR for a given tap sequence.
\par
The tap sequence defines which bits in the current state will be combined to determine the input bit for the next state, generally using module-2 addition ($\it exclusive~or$). Tap sequences can be represented by irreducible polynomials over ${\bf GF}(2)$ ($\it Galois~Field$ of order 2) -- polynomials with coefficients of either $0$ or $1$, which cannot be represented as the product of two or more polynomials. For example, $x^2+1$ is not irreducible over ${\bf GF}(2)$ because it can be factored into $x+1$:
\begin{equation}
(x+1)(x+1)=x^2+2x+1 \equiv x^2+1.
\end{equation}
Considering an irreducible polynomial, the corresponding tap sequence is given by the exponents whose polynomial coefficients are 1 \cite[] {dinan}. M-sequences are binary pseudo-random sequences of length $(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)$, with $gcd\{n,4\} = 1$, 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 autocorrelation of an m-sequence on the top and its cross-correlation with its preferred pair computed with $m=5$. 
\plot{prefpairs}{width=3in}{Correlation functions of preferred pairs of m-sequences. Top$:$ autocorrelation of an m-sequence. Bottom$:$ cross-correlation of a preferred pair of m-sequence.}

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. As 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 sequence, $\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 three 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 code.
\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 that generated with 84 circular shifts. 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. 
\plot{gold184}{width=\textwidth}{Correlation functions of Gold codes. a) Autocorrelation of Gold code generated with one circular shift of the preferred pair of m-sequence. b) Gold code generated with one circular shift of the preferred pair of m-sequence cross-correlated with that generated with 84 circular shifts. The peak of the cross-correlation occurs at lag 84.}
\plot{rand}{width=\textwidth}{Correlation functions of conventional random codes. a) Autocorrelation of the conventional random code. b) A cross-correlation.}

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 embedded 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 and the areal shot migration. Super-areal data are comprised of the collection of 10 modeling experiments initiated at every 10th CMP coordinate.
\par
Figure \ref{fig:perma} shows the areal shot migration of data generated by the prestack exploding-reflector modeling with no phase encoding applied. The panel on the left is the zero-subsurface-offset section, and the panel on the right is a SODCIG. 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 type-2 crosstalk. The crosstalk is periodic with half of the spacing of the modeling experiments in a super-areal shot. 
\plot{perma}{width=4.0in}{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.}

Figure \ref{fig:perm0} shows the areal shot migration of data generated by the prestack exploding-reflector modeling with no combination of modeling experiments into super-areal shots. This result represents the ideal image we would 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. 
\plot{perm0}{width=4in}{Areal shot migration of synthesized data with no combination of the modeling experiments into super-areal data.}

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.
\plot{conv1r}{width=4in}{Areal shot migration of one realization of synthesized data with conventional random-phase encoding.}
Figure~\ref{fig:conv5r} shows the migration of 5 realizations of conventional random encoding modeling.  
\plot{conv5r}{width=4in}{Areal shot migration of five realizations of synthesized data with conventional random-phase encoding.}

Figure~\ref{fig:opti1} shows the migration of one realization of data encoded with codes computed by optimization, and Figure~\ref{fig:opti5} shows the migration of five realizations. No striking difference exists when comparing these results with those of the conventional random codes in Figures~\ref{fig:conv1r} and \ref{fig:conv5r}, respectively.
\plot{opti1}{width=4in}{Areal shot migration of one realization of synthesized data with optimized random-phase encoding.}
Both encoding methods disperse the crosstalk energy in a similar manner. (more to come ...)
\plot{opti5}{width=4in}{Areal shot migration of five realizations of synthesized data with optimized random-phase encoding.}

Crosstalk attenuation is incomplete when using conventional or optimized random codes because their autocorrelations are not perfect spike, nor are their cross-correlations zero everywhere. 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, among the available codes, the ones 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 for which the Gold codes were selected sequentially. This means that when applying the multi-offset imaging condition, cross-correlation peaks can occur at a distance less than that enclosed by the SODCIG. The crosstalk is present everywhere in the SODCIG and potentially obscures the moveout information for migration velocity analysis. It appears also as `ghost' reflectors in the zero-subsurface-offset section.  It is interesting to notice that the crosstalk in Figure~\ref{fig:gold1x} still has its apexes at the same offsets as those present in Figure~\ref{fig:perm0}. However, the apexes are displaced in depth. In particular, the crosstalk that occurs close to the edge of the SODCIG in Figure~\ref{fig:perm0} is almost completely shifted out of the SODCIG. The amount of shift is defined by the lag where the cross-correlation peaks. In this example, the super-areal data  comprises Gold-encoded modeling experiments in which cross-correlation peaks are separated by multiples of 10 lags.
\plot{gold1x}{width=4in}{Areal shot migration of one realization of synthesized data with Gold phase encoding. Gold codes are selected sequentially such that the cross-correlation peak separation is multiple of 10.}

For the simple case of constant velocity, the depth shifts, $\delta z$, are given by
\begin{equation}
\delta z = \frac{0.5*v*n_{\lambda}}{n_{\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 choose 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.
\plot{gold1r}{width=4in}{Areal shot migration of one realization of synthesized data with Gold phase encoding. Gold codes are randomly selected.}

As before, migrating more realizations of randomly selected, Gold-encoded data further attenuates crosstalk. Figure~\ref{fig:gold5r} shows the migration of 5 realizations of randomly selected, Gold-encoded data. Comparison with Figure~\ref{fig:gold1r} shows that much of the remaining crosstalk energy has been attenuated. 
\plot{gold5r}{width=4in}{Areal shot migration of five realizations of synthesized data with Gold phase encoding. Gold codes are randomly selected.}

Considering that the crosstalk is shifted in depth, as can be seen in Figure~\ref{fig:gold1x}, and that the amount of shift is determined by the lag where the cross-correlation peaks (equation~\ref{eq:zshift}), one can choose the Gold codes according to a suitable interval that completely shifts the crosstalk out of the SODCIG. This strategy shares similar idea as the linear-phase encoding of \cite[]{romero} which aims to shift the crosstalk out of the migration domain by using a linear function of frequency.  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. The crosstalk is virtually eliminated. Compare with Figure~\ref{fig:perm0}.
\plot{gold1o}{width=4in}{Areal shot migration of one realization of synthesized data with Gold phase encoding. Gold codes are selected such that the crosstalk is shifted out of the migration domain.}

Actually, Gold phase encoding 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 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. So, 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}
Equation \ref{eq:linph} is the same as equation 26 of Romero's paper.
\par

This strategy can be less computationally demanding than using conventional random codes, because just one realization is necessary to achieve a good result, while even using more realizations of conventional random encoding does not produce an image with similar quality. Its use in more complex situations, where the shift in depth of the crosstalk is not easily predicted, however, deserves more research.


\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 image quality but also cost. While the elimination of the crosstalk is achieved with only one realization of Gold codes, when using conventional random codes more realizations are needed to obtain reasonable images, increasing the cost of modeling and migration. We also showed the equivalence of Gold phase encoding and linear phase encoding.
\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.

\bibliographystyle{seg}
\bibliography{SEP,ref134a}
