\def\figdir{./Fig}
\chapter{Pre-Stack Exploding-Reflector Model}
\label{chap:chap02}
\def\w{\omega}
\def\bX{\boldsymbol{x}}
\def\bxi{\boldsymbol{\xi}}

%\begin{abstract}
This chapter introduces the pre-stack exploding-reflector model (PERM). PERM uses $\it exploding$ reflectors as the initial conditions to synthesize data, in a manner similar to the well-known exploding-reflector model (ERM). However, PERM also considers reflectivity as a function of subsurface-offsets, as opposed to the zero-subsurface offset initial condition used by ERM. In this sense, PERM is a generalization of ERM since wave-equation migration of PERM data can generate a pre-stack image, which is not achievable with ERM. The modeling of PERM data can be performed using any wave propagation scheme; here I use the one-way wave-equation. As in the exploding-reflector model, the modeling experiments can be combined, potentially decreasing data size in orders of magnitude when compared to the original shot records.  Data size reduction is very appealing for migration velocity analysis, especially using wavefield extrapolation methods in 3D projects. After describing the exploding-reflector concept, I generalize it by describing the theory of PERM. The usefulness of PERM to migration velocity analysis, is demonstrated through migration examples.  
%\end{abstract} 
%
\section{Introduction} 
%
Migration is applied to seismic data to generate an image of the subsurface. For many years, migration was applied only in the post-stack domain, using the idea of exploding reflectors \cite[]{loewenthal01,iei}. The simple but powerful concept of the exploding-reflector model (ERM) states that a zero-offset time section can be considered as the recording at the surface of wavefields generated by simultaneous explosions of all points in the subsurface. The strength of the explosions is proportional to the reflection coefficient and, because the wavefield propagates from the subsurface to the surface, to obtain correct kinematics the medium velocity must be halved.
\par
Migration in the post-stack domain assumes that data is transformed to some approximation of the zero-offset acquisition geometry. Because of the required transformation to zero-offset, many algorithms have been developed like dip moveout (DMO) \cite[]{hale:741,black:1109,liner:182}, migration to zero-offset (MZO) \cite[]{tygel:557,bleistein:112}, and common-reflection surface \cite[]{gelchinsky,cruz:979}. DMO and MZO can be considered very particular cases of the more general azimuth moveout (AMO) \cite[]{biondi:574}, which must be applied to 3D data prior to common-azimuth migration \cite[]{biondi:1822}.
\par
In fairly simple geology, ERM perfectly matches the kinematics of the zero-offset data that would have been recorded with coincident source and receivers at the surface. This means that for ERM to hold, it is necessary that the downgoing path from the source to a point in the subsurface be the same as the upgoing path from the point in the subsurface back to the receiver. However, this assumption is often invalid in areas of geological complexity. Post-stack imaging below salt bodies, for instance, can give rise to numerous migration artifacts, in the presence of prismatic and strongly refracted waves \cite[]{peles:2455}. Therefore, post-stack migration does not produce reliable images in the presence of strong lateral velocity variation.
\par
In areas of complex geology pre-stack depth migration becomes mandatory not only for imaging purposes but also for velocity estimation. In such areas, migration by wavefield extrapolation has been widely used to produce the final image because it properly handles complex distortions of the wavefields. However, due to the high computational cost, wavefield extrapolation methods are rarely used to estimate the migration velocity model in 3D projects \cite[]{fei:2283}, where ray-based methods are the industry standards. In addition to the lower cost, ray-based methods are very flexible with respect to strategies for defining the velocity model \cite[]{stork:680,kosloff:1511,kosloff:1815,billette:1889}. But despite their advantages, ray-methods do not satisfactorily describe complex wave propagation in the presence of large lateral velocity contrasts. In this case, a more complete description of the wavefield complexities is needed, and therefore we face the challenge of decreasing the cost of migration velocity analysis by wavefield extrapolation while maintaining its robustness. 
\par
A typical way of decreasing the cost of wavefield extrapolation is to reduce the amount of input data. Data size can be reduced by selecting a smaller number of shots. Figure \ref{fig:perm01} shows pre-stack images computed with the multi-offset imaging condition \cite[]{james2002} using 380 shots 24 m apart (Figure \ref{fig:perm01}a) and 23 shots 384 m apart (Figure \ref{fig:perm01}b). The shot records are modeled with the one-way wave equation, using a smoothed version of the Marmousi velocity model. The same velocity model is used to migrate the data. The right panel is the zero-subsurface-offset section and the left panel is the subsurface-offset-domain common-image gather (SODCIG) selected at $x = 2760 m$ on the right panel. Notice the slanted lines in the SODCIG of Figure \ref{fig:perm01}b representing the angles according to which reflectors are illuminated. When using all the shots, the slanted lines constructively interfere and the energy is focused at zero-subsurface-offset (Figure \ref{fig:perm01}a).
\activeplot{perm01}{width=\textwidth}{}{a) Shot-profile migration of 380 off-end shots 24 m apart. b) Shot-profile migration of 23 off-end shots 384 m apart. Both images were computed with the correct velocity model. Notice slanted lines present in Figure \ref{fig:perm01}b.}
\par
By selecting only a few shots, the migration output shows poor angular coverage. An alternate way of reducing data size is by combining shot-profiles into areal shots based on the linearity of wavefield propagation. Data is multiplied by a comb function and stacked to originate one areal shot. The comb function is shifted until all shots are selected. The original angular coverage is maintained if the period of the sampling function is sufficiently big (Figure \ref{fig:perm02}a). However, if many shots are combined, crosstalk is generated since unrelated shots and receiver wavefields are cross-correlated during imaging (Figure \ref{fig:perm02}b). Compare Figures \ref{fig:perm02}a and \ref{fig:perm02}b with Figure \ref{fig:perm01}a. Crosstalk can completely overwhelm the reflectors, and the kinematic information for migration velocity analysis can be lost. 
\activeplot{perm02}{width=\textwidth}{}{a) Shot-profile migration of 95 areal shots resulting from the combination of 4 shot-profiles 2256 m apart. b) Shot-profile migration of 23 areal shots resulting from the combination of 16 shot-profiles 564 m apart. Both images were computed with the correct velocity model. Notice crosstalk occurring periodically in the SODCIG of Figure \ref{fig:perm02}b.}
%However, crosstalk is responsible for the poor image quality of Figure \ref{fig:perm02}b. In addition, it is clear the repetition pattern at the subsurface-offset gather. In this case, the crosstalk occurs at absolute subsurface-offsets approximately half the distance between the combined shots.
\par
The combination of wavefields is implicit in ERM. When reflectors are allowed to simultaneously explode, linearity of wavefield propagation is evoked to combine wavefields initiated at every point in the subsurface. However, crosstalk is not generated in post-stack migration based on ERM since no cross-correlation of wavefields is performed because the imaging condition is a simple summation over frequency, which extracts the image at zero time of wavefield propagation \cite[]{claerbout71}. Therefore, it is natural to consider the combination of wavefields in any generalization of ERM. 
\par
Combination of wavefields is exploited by the prestack-exploding reflector model (PERM) \cite[]{biondi:3056} to significantly reduce the data size used in migration velocity updates. Like ERM, PERM uses the concept of exploding-reflectors to propagate wavefields. However, instead of considering reflectivity only as a function of the spatial coordinates, reflectivity is also parameterized as a function of subsurface-offsets. Its elementary idea is to synthesize data necessary to correctly image a single reflector of an isolated SODCIG, so that migration of PERM data shows correct kinematics needed for performing velocity updates.
% Since modeling data for every single point in the subsurface can potentially increase data size, combination of modeling experiments is used to perform data reduction in a similar way as that for the previous shot-profile migration examples.
\par
An interesting feature of PERM is that it can be used in a target-oriented way, since the modeled data contains all the information  necessary to image a predefined region of the subsurface. This concept is present in different methods, such as controlled illumination \cite[]{Rietveld01,Rietveld02} and common-focus point (CFP) technology \cite[]{berkhout:937, berkhout:954, thorbecke:S273}. 
%One trace $n(t,x_s)$ of a CFP-gather results from the convolution in time between the shot gather at $x_s$ and a focusing operator, and horizontal stack which enhances the stationary phase contributions. The focusing operator is computed by modeling the seismic response of one point at the subsurface. In contrast with CFP, however, PERM does not operate on the original shot gathers. PERM data result from the wave propagation of the pre-stack image at one point at the subsurface. In this sense, PERM data is more closely related to the CFP operator itself than to the CFP data.
\par
In this chapter, I introduce and further develop PERM. It will be shown that data synthesized with PERM has the kinematic information necessary to perform migration velocity analysis. Although, in this thesis, migration velocity analysis is performed using wave-equation tomography, ray-based tomography can also be used in conjunction with PERM data.  Before introducing the theory of PERM, I will first describe the exploding-reflector concept of which PERM is a generalization. The usefulness of PERM for migration velocity analysis will be illustrated by comparing areal-shot migration of PERM wavefields with that of shot-profile migration for simple bidimensional models. The extension to a tridimensional (3D) medium is straightforward. In 3D, it will be shown that under the common-azimuth approximation \cite[]{biondi:1822}, PERM drastically decreases data size.

%, different strategies have been formulated to combine shot records based on the linearity of wavefield extrapolation. After the introduction of synthesis of plane waves \cite{Schultz}, several authors used this technique to perform imaging and velocity analysis using asymptotic methods \cite[]{akbar:1068} and wavefield extrapolation methods \cite[]{whitmore,mosher1998,zhang:E21}. \cite{liu01} provide a general framework to evaluate plane-wave decomposition in pre-stack source plane-wave migration and show the equivalence of shot profile and plane-wave migrations. Further generalizations of Schultz and Claerbout's method use wavefields of arbitrary shape in the subsurface to be upward propagated to the surface, where they define synthesis operators to be convolved with the shot profiles to generate areal shots \cite[]{Rietveld01,Rietveld02}. Migration of randomly phase-encoded shot-records \cite[]{romero} is a technique that enables migrating combined shots with minimal  cross-talk. Recently, random-phase encoding has been used to decrease the computational cost of waveform inversion \cite[]{ben-hadj-ali:2288}. The combination of phase-encoded and plane-wave areal shots can be used to drastically decrease the computational effort \cite[]{Sun01}. 
%\par
%Since the exploding-reflector assumes zero-offset propagating energy, any velocity estimation process which uses only the stacked section or near offset traces has to rely only on focusing and defocusing of the events \cite[]{harlan:1869} in the zero-subsurface-offset image. This procedure, however, is prone to problems caused by the ambiguity that exists between reflector curvature and velocity \cite[]{biondo09}.
%\par
%Fortunately, pre-stack images computed with wavefield extrapolation methods \cite[]{ottolini:237,prucha:824,rickett:883} contain the moveout information needed to update the velocity model. In shot-profile and areal shot migration, the pre-stack image is computed by time cross-correlation and extraction of the zero time lag of spatially shifted source and receiver wavefields. If migration velocity is not accurate, seismic energy spreads along curved trajectories at non-zero-subsurface-offsets instead of being focused at zero-subsurface-offset. % In addition, the pre-stack image can also be computed by time cross-correlation of time shifted source and receiver wavefields \cite[]{sava:S209}. In this case, if the velocity is inaccurate, seismic energy moves away from the zero time-lag.
\par

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{exploding-reflector model}
To describe the exploding-reflector model (ERM) let us analyze the modeling of seismic data under the acoustic approximation. Here, we consider the Born or single-scattering approximation. This consideration leads us to a linear operator whose adjoint is the migration operator.
\par
Let us start with the constant-density acoustic wave equation for a single temporal frequency $\omega$
\begin{eqnarray}
\left (\nabla^2+\omega^2 s^2({\bf x}) \right) P({\bf x},\omega) = 0,
\label{eq:acweq0} 
\end{eqnarray}
where $\nabla^2$ is the Laplacian operator, $s({\bf x})$ is the slowness field, and $P({\bf x},\omega)$ is the wavefield. Note that ${\bf x}=(x,y,z)$ is the vector of spatial coordinates. By introducing reflectivity as $r({\bf x}) = 1 - \frac{s^2({\bf x})}{s_0^2({\bf x})}$, where $s_0({\bf x})$ is a smooth background slowness, equation \ref{eq:acweq0} can be written as
\begin{eqnarray}
\left (\nabla^2+\omega^2 s_0^2({\bf x},\omega) \right) P({\bf x},\omega) \approx \omega^2s_0^2({\bf x}) r({\bf x}) P({\bf x},\omega).
\label{eq:acweq1}
\end{eqnarray}
%where $\nabla^2$ is the Laplacian operator, $s({\bf x})$ is the slowness field, and $D({\bf x},\omega)$ is the wavefield. Note that ${\bf x}=(x,y,z)$ is the vector of spatial coordinates. The slowness field can be described by $s({\bf x}) = s_0({\bf x}) + \Delta s({\bf x})$, where $s_0({\bf x})$ is a smooth background slowness and $\Delta s({\bf x})$ is a small slowness perturbation. Therefore, equation \ref{eq:acweq0} can be written as
\par
The total wavefield $P({\bf x},\omega)$ can be considered as the sum of an incident background wavefield $P_0({\bf x},\omega)$ with a perturbed or scattered wavefield $\Delta P({\bf x},\omega)$. Therefore, we can write
\begin{eqnarray}
\left (\nabla^2+\omega^2 s_0^2({\bf x}) \right) \Delta P({\bf x},\omega) \approx \omega^2s_0^2({\bf x}) r({\bf x}) P({\bf x},\omega),
\label{eq:acweq2}
\end{eqnarray}
given that the background wavefield $P_0({\bf x},\omega)$ is the solution of equation \ref{eq:acweq1} using the background slowness $s_0({\bf x})$. Notice that equation \ref{eq:acweq2} is a non-linear relationship between $r({\bf x})$ and $\Delta P({\bf x},\omega)$. A linear operator can be derived by using Green's function and the Born (weak scattering) approximation. The background Green's function $G_0({\bf x'},{\bf x},\omega)$ is the solution of equation \ref{eq:acweq0} using the background slowness in the presence of a point source $\delta({\bf x}-{\bf x'})$ at ${\bf x'}=(x',y',z')$:
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \nabla^2+\omega^2 s_0^2({\bf x}) \right) G_0({\bf x'},{\bf x},\omega) = 0 \\
G_0({\bf x'},{\bf x}={\bf x'},\omega) = \delta({\bf x}-{\bf x'}) \end{array} \right. .
\label{eq:acweq3}
\end{eqnarray}
Multiplying equation \ref{eq:acweq3} by $\left [ \omega^2s_0^2({\bf x}) r({\bf x}) P({\bf x},\omega) \right ]$, integrating with respect to ${\bf x'}$ over the volume $V$ in the subsurface, and comparing the result with equation \ref{eq:acweq2}, we see that
\begin{eqnarray}
\Delta P({\bf x},\omega) \approx {\displaystyle\int\nolimits_V} \omega^2s_0^2({\bf x'}) r({\bf x'}) G_0({\bf x'},{\bf x},\omega) P({\bf x'},\omega) d{\bf x'}.
\label{eq:acweq4}
\end{eqnarray}
Under the weak scattering assumption, the total wavefield $P({\bf x'},\omega)$ can be approximated by the background wavefield $P_0({\bf x'},\omega)$, and the linear relationship between the scattered wavefield $\Delta P({\bf x},\omega)$ and reflectivity $r({\bf x})$ reads
\begin{eqnarray}
\Delta P({\bf x},\omega) \approx {\displaystyle\int\nolimits_V} \omega^2s_0^2({\bf x'}) r({\bf x'}) G_0({\bf x'},{\bf x},\omega) P_0({\bf x'},\omega) d{\bf x'}.
\label{eq:acweq5}
\end{eqnarray}
\par
Using equation \ref{eq:acweq5} and assuming that the background slowness field and the reflectivity distribution are known, the scattered wavefields $\Delta P({\bf x_s},{\bf x_r},\omega)$ measured by receivers at ${\bf x_r}=(x_r,y_r,z_r)$ due to sources located at ${\bf x_s}=(x_s,y_s,z_s)$ are
\begin{eqnarray}
\Delta P({\bf x_s},{\bf x_r},\omega) \approx {\displaystyle\int\nolimits_V} \omega^2s_0^2({\bf x}) G_0({\bf x_s},{\bf x},\omega) r({\bf x}) G_0({\bf x},{\bf x_r},\omega) d{\bf x},
\label{eq:acweq6}
\end{eqnarray}
which amounts to convolving the source Green's function with the reflectivity, then with receiver Green's function, and summing the contributions of all points within $V$. By using one-way propagators, we downward propagate the source wavefield, convolve it with the reflectivity, and upward propagate this result up to the surface.
\par
Zero-offset data is obtained by selecting from $\Delta P({\bf x_s},{\bf x_r},\omega)$ traces with ${\bf x_s}={\bf x_r}$,
\begin{eqnarray}
\left. \Delta P({\bf x_s},\omega) \right|_{x_s=x_r} \approx {\displaystyle\int\nolimits} \omega^2s_0^2({\bf x}) r({\bf x}) \left. G_0({\bf x_s},{\bf x},\omega)\right|_{x_s=x_r} G_0({\bf x},{\bf x_s},\omega) d{\bf x}.
\label{eq:acweq7}
\end{eqnarray}
\par
ERM synthesizes zero-offset data by initiating virtual sources located on the reflectors. The fundamental consideration is that the downgoing and upgoing rays of a zero-offset source and receiver pair follow the same path. The exploding-reflector wavefield is propagated with half of the actual background velocity and the virtual sources explode with strength proportional to the reflection coefficient. Assuming that the background slowness and the reflectivity distribution are known, the wavefield generated by the exploding-reflector model $P_{ERM}({\bf x},\omega)$ propagates according to the following one-way wave-equation:
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 s^2({\bf x})-|{\bf k}| ^2}\right) P_{ERM}({\bf x},\omega) = r({\bf x}) \\
P_{ERM}(x,y,z=z_{\rm max},\omega) = 0 \end{array} \right.,
\label{erm01}
\end{eqnarray}
\par
ERM is kinematically correct only if source and receiver Green's functions are equal for ${\bf x_s}={\bf x_r}$ at a point $\bf x$ in the subsurface. Notice that this condition is unlikely to be fulfilled in the presence of strong lateral velocity variations \cite[]{peles:2455}. Moreover, ERM does not model multiples and prismatic reflections \cite[]{iei}. 
\par
In equation \ref{erm01}, the reflectivity $r({\bf x})$ acts as the initial condition for the wavefield propagation. Notice that a migrated image can replace the reflectivity in the initial condition. Taking into consideration the imaging principle \cite[]{claerbout71} and the computation of pre-stack images by wave-extrapolation methods \cite[]{iei,james2002}, ERM, as initially formulated, implicitly assumes that all the seismic energy is perfectly focused at the zero-subsurface offset and that it is sufficient to parameterize the migrated image only as a function of the position vector $\bf x$. This implies that the slowness field is accurate and illumination of the subsurface is sufficiently good. When this is not the case, the migrated image must be described not only as a function of ${\bf x}$ but also as a function of the subsurface offset or aperture angle. 
\par
%In addition, finite-offset data only can be modeled with ERM under the very strict assumption of an one-dimensional velocity structure within the region covered by the source and receiver ray-paths, as explored by the extended ERM of \cite{loewenthal02}. Therefore, ERM concept can only be used to image pre-stack data under the validity of time migration.

The generalization of ERM proposed in this thesis does not aim at computing pre-stack data. Rather, it synthesizes data whose migration with wavefield-extrapolation methods enables the computation of a pre-stack image. PERM uses the ERM concept with a poorly focused pre-stack image as initial condition to propagate wavefields. The modeled data is, potentially, orders of magnitude smaller than the original shot records and contains all necessary kinematic information to update the velocity model using ray-based tomography or, as used in this thesis, wave-equation-based tomography. This characterizes the pre-stack-exploding-reflector model, which will be described in the next section. 

%The generalization of ERM proposed in this thesis does not aim at computing pre-stack data. In contrast, it uses the ERM concept having as initial condition a pre-stack reflectivity to synthesize data whose migration with wave-extrapolation methods enables the computation of a pre-stack image. The synthesized data is, potentially, orders of magnitude smaller than the original shot records and contains all necessary kinematic information to update the velocity model using ray-based tomography or, as used in this thesis, wave-equation based tomography. This characterizes the pre-stack-exploding-reflector model, which will be described in the next section. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{pre-stack-exploding-reflector model} 
The fundamental idea of PERM is to model data which describes the correct kinematics of an isolated SODCIG. Many shot records contribute to form the image at a point in the subsurface. Therefore, to model data using the conventional one-way modeling, we would perform several modeling experiments consisting of downward continuing the source wavefield initiated by point sources at the surface, convolving the propagated source wavefield with the SODCIG, and upward continuing the convolution result up to the surface. As we do not know beforehand which shots contribute to forming the image at a point in the subsurface, we would have to model every shot originally present in the original dataset. % according to equation \ref{eq:perm01},
%\begin{eqnarray}
%\begin{array}{l}
%\Delta P({\bf x_s},{\bf x_r},\omega;{\bf x_m}) \approx {\displaystyle\int\nolimits_V} {\displaystyle\int\nolimits_{\bf h}} \omega^2s_0^2({\bf x}) G_0({\bf x_s},{\bf x-h},\omega;{\bf x_m}) r({\bf x,h;x_m}) \\
%~~~~~~~~~~~~\times G_0({\bf x+h},{\bf x_r},\omega;{\bf x_m}) d{\bf x} d{\bf h}, \end{array}
%\label{eq:perm01}
%\end{eqnarray}
%where ${\bf x_m}=(x_m,y_m)$ is the location of the isolated SODCIG; $G_0({\bf x_s},{\bf x},\omega;{\bf x_m})$ and $G_0({\bf x},{\bf x_r},\omega;{\bf x_m})$ are the source and receiver Green's functions, respectively.  Equation \ref{eq:perm01} is similar to equation \ref{eq:acweq6} except that reflectivity $r({\bf x,h};{\bf x_m})$ is also a function of subsurface-offset $\bf h$ and the Green's functions are shifted laterally by $\bf h$. This accommodates modeling data starting from a poorly focused image due to migration with an inaccurate velocity. Notice that the subsurface-offset $\bf h$ can be parameterized as ${\bf h}=(h_x,h_y,h_z)$. In this thesis, when describing 2D problems ${\bf h}=(h_x)$, and for the 3D case,  ${\bf h}=(h_x,h_y)$ without considering the computation of the vertical subsurface offset $h_z$ as introduced by \cite{biondi:1284}.
\par
Ideally, instead of performing many modeling experiments, we would like to synthesize a small amount of data with the condition that migration has the same kinematics as the initial SODCIG. This can be achieved by using a strategy similar to ERM. To compute SODCIGs using the new modeled data, we need to synthesize source and receiver wavefields.
\par
The modeling of source and receiver wavefields can be carried out by any wavefield-continuation scheme. Here, we use the following one-way wave equations:
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}-i\sqrt{\omega ^2 s_0^2({\bf x})-|{\bf k}| ^2}\right) D_{\scriptscriptstyle P}({\bf x},\omega;{\bf x_m}) = I_D({\bf x_m},{\bf h}) \\
D_{\scriptscriptstyle P}(x,y,z=z_{\rm max},\omega;{\bf x_m}) = 0 \end{array} \right.,
\label{eq:perm_bD}
\end{eqnarray}
and
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 s_0^2({\bf x})-|{\bf k}| ^2}\right) U_{\scriptscriptstyle P}({\bf x},\omega;{\bf x_m}) = I_U({\bf x_m},{\bf h}) \\
U_{\scriptscriptstyle P}(x,y,z=z_{\rm max},\omega;{\bf x_m}) = 0 \end{array} \right.,
\label{eq:perm_bU}
\end{eqnarray}
where $I_D({\bf x_m},{\bf h})$ and $I_U({\bf x_m},{\bf h})$ is the isolated SODCIG at the horizontal location ${\bf x_m}$ for a single reflector, suitable for the initial conditions for the source and receiver wavefields, respectively. The subsurface-offset $\bf h$ can be parameterized as ${\bf h}=(h_x,h_y,h_z)$. In this thesis, when describing 2D problems ${\bf h}=(h_x)$, and for the 3D case,  ${\bf h}=(h_x,h_y)$ without considering the computation of the vertical subsurface offset $h_z$ as introduced by \cite{biondi:1284}. The initial conditions are obtained by rotating the original unfocused SODCIGs according to the apparent geological dip of the reflector. For dipping reflectors, this rotation maintains the velocity information needed for migration velocity analysis. The rotation of SODCIGs will be described later in this chapter. If the initial condition has energy focused at the zero-subsurface offset or if no pre-stack information is available, the pre-stack image can be parameterized only by its spatial coordinates, and PERM is equivalent to ERM.

\par
%The pseudo-code describing PERM is
%\begin{algorithm}
%\caption{ Pre-stack exploding-reflector model of source and receiver wavefields (equations \ref{eq:perm_bD} and \ref{eq:perm_bU}). } 
%\centering
%\small
%\begin{algorithmic}
%\FORALL {spatial positions $n_x*n_y$}
%\STATE read pre-stack image I_{(D,U)}
%\FORALL {frequencies~$n_{\omega}$}
%\STATE wavefield=0.
%\FORALL {depths: z{max},2,-1 }
%\STATE wavefield=wavefield+FFT(I_{(D,U)})
%\STATE upward propagate wavefield
%\STATE wavefield=wavefield*taper
%\ENDFOR
%\ENDFOR
%\ENDFOR
%\end{algorithmic} 
%\normalsize 
%\label{tab:alg_perm}
%\end{algorithm}  
Let us now illustrate the generation of PERM data synthesized from a single SODCIG. We start with a pre-stack image computed with shot-profile migration of 401 split-spread shots at every 10 m and maximum offset of 2250 m, using a $10\%$ lower velocity (Figure \ref{fig:refpl01}). The model consists of one reflector at 750 m depth embedded in a medium with constant velocity of  1000 m/s. The pre-stack image has 81 subsurface-offsets ranging from -400 m to 400 m. Notice the poor focusing of energy around the zero-subsurface offset due to inaccurate velocity. 
\par
The SODCIG at $ x_m = 0~m$ was used as the initial condition for modeling the corresponding pair of PERM source and receiver wavefields using the same inaccurate velocity. The wavefields are upward propagated according to equations \ref{eq:perm_bD} and \ref{eq:perm_bU}. The PERM data is shown in Figure \ref{fig:refpl02}. Notice that the receiver wavefield (Figure \ref{fig:refpl02}a) occurs at positive times while the source wavefield (Figure \ref{fig:refpl02}b) occurs at negative times. According to the imaging principle \cite[]{claerbout71}, reflectors explode at time zero. This is the time at which the source wavefield impinges on the reflector. Because the receiver wavefield exists after the source wavefield has reached the reflector, the areal receiver data $U_{\scriptscriptstyle P}(x,y,z=0,\omega;{\bf x_m})$ is upward propagated forward in time. For the same reason, the areal source data $D_{\scriptscriptstyle P}(x,y,z=0,\omega;{\bf x_m})$ is upward propagated backward in time.

\activeplot{refpl01}{width=\textwidth}{}{Shot-profile migration of 401 split-spread shots 10 m apart with a $10\%$ slower velocity. The model consists of a horizontal reflector embedded in constant velocity of 1000 m/s.}

\activeplot{refpl02}{width=\textwidth}{}{Data synthesized by PERM having as the initial condition the SODCIG at $x_m=0$ m. a) The receiver wavefield. b) The source wavefield.}

\par
Areal-shot migration of data from Figure \ref{fig:refpl02} with the same inaccurate velocity shows kinematics at near subsurface-offsets similar to that of the original shot-profile migration (compare Figures \ref{fig:refpl03} and \ref{fig:refpl01}). For farther subsurface-offsets, however, energy is not adequately imaged. This can be easily explained by analyzing how the image $I({\bf x},{\bf h})$ is formed when applying the multi-offset imaging condition 
\begin{eqnarray}
I({\bf x},{\bf h}) = \sum_{\omega} \sum_{{\bf x}_s} D^{*}({\bf x}-{\bf h},{\bf x}_s,\omega) U({\bf x}+{\bf h},{\bf x}_s,\omega). 
\label{eq:ic_bI}
\end{eqnarray}
In equation \ref{eq:ic_bI}, $D({\bf x}-{\bf h},{\bf x}_s,\omega)$ and  $U({\bf x}+{\bf h},{\bf x}_s,\omega)$ are the source and receiver wavefields, respectively, shifted laterally by ${\bf h}$. The image is formed by summing the contributions of all shot positions  ${\bf x}_s$ and all frequencies. Notice that `$^{*}$' stands for the complex conjugate. From equation \ref{eq:ic_bI} we can see that farther subsurface-offsets get contributions from the cross-correlation of wavefields within a neighborhood around ${\bf x}$. Specifically for the PERM case, farther subsurface-offsets of the areal-shot migrated image at ${\bf x_m}$ (the location of the isolated SODCIG) also get contributions from the cross-correlation of wavefields generated at SODCIGs within a neighborhood around ${\bf x_m}$. Therefore, to ensure that the areal-shot migrated image has kinematics at all subsurface-offsets similar to that in the original isolated SODCIG, we need to model PERM data from a set of SODCIGs around ${\bf x_m}$. From equation \ref{eq:ic_bI}, it is clear that the maximum distance at which the wavefields still contribute to form the image at the maximum subsurface-offset is given by $2(max\left | h_x \right |,max\left | h_y \right |)$.
\par

\activeplot{refpl03}{width=\textwidth}{}{Areal-shot migration of PERM data shown in Figure \ref{fig:refpl02} with a $10\%$ slower velocity. By comparing with Figure \ref{fig:refpl01} we see that far subsurface-offsets are not properly imaged. }

The areal-shot migration of PERM data synthesized by isolated SODCIGs within the interval $(-2max\left | h_x \right |,2max\left | h_x \right |)$ is shown in Figure \ref{fig:refpl04}. By using more data, energy is adequately imaged at far subsurface-offsets (compare with Figures \ref{fig:refpl01} and \ref{fig:refpl03}). To further understand the behavior of the areal-shot migrated image, let us examine the reflection  angle-domain common-image gathers (ADCIGs) \cite[]{sava:1065}. ADCIGs computed from the SODCIGs of Figures \ref{fig:refpl01}, \ref{fig:refpl03} and \ref{fig:refpl04} attest to the more accurate imaging when migrating data modeled from the set of SODCIGs around  ${\bf x_m}$ (Figure \ref{fig:refpl05}). Although the ADCIG from the the image computed with a single pair of PERM data (Figure \ref{fig:refpl05}b) shows reasonable kinematics, the amplitude of wide-aperture angles is weaker than that of the original ADCIG (Figure \ref{fig:refpl05}a). Notice that the amplitude behavior of the ADCIG computed with several areal shots from SODCIGs within the neighborhood of ${\bf x_m}$ (Figure \ref{fig:refpl05}c) better matches that of the original isolated SODCIG. 

\activeplot{refpl04}{width=\textwidth}{}{Areal-shot migration of PERM data having a set of isolated SODCIGs around $x_m=0$ m as the initial condition with a $10\%$ slower velocity. By comparing with Figure \ref{fig:refpl01} we see that the kinematics of far subsurface-offsets is properly recovered.} 
\activeplot{refpl05}{width=\textwidth}{}{ADCIGs selected at $x_m=0$ m. a) Computed from the shot-profile migration; b) computed from the areal-shot migration of one pair of PERM data modeled from the SODCIG at $x_m=0$ m; and c) computed from the areal-shot migration of pairs of PERM data modeled from a set of SODCIGs around $x_m=0$ m. Notice that although the kinematics is similar, the amplitudes in c) better match those of a).}
\par
%The localized nature of ... grid tomography
%\par
If, instead of using the incorrect migration velocity, we input the correct migration velocity to the areal-shot migration, energy nicely focuses at zero-subsurface offset (Figure \ref{fig:refpl06}). This property will be used to perform migration velocity updates in Chapter 3.

\activeplot{refpl06}{width=\textwidth}{}{Areal-shot migration of PERM data having a set of isolated SODCIGs around $x_m=0$ m as the initial condition with the correct velocity. Energy nicely focuses at zero-subsurface offset.}
\par
In the previous examples we saw that PERM data contains all the kinematic information needed to perform migration velocity analysis. However, by carefully examining the migration results of the original data and PERM data we can see that the latter has stronger side lobes. 
% \sideplot{refpl07}{width=4in}{Comparison of the waveforms for the shot-profile migration of the original data (dashed blue) and the areal-shot migration of PERM data (continuous red). Note stronger side lobes in the areal-shot migration of PERM data.}
\par
The stronger side lobes are due to the squaring of the wavelet. It occurs because PERM source and receiver wavefields contain the waveform of the original image. Let us consider, for simplicity, a plane reflector so that the initial conditions are the same for modeling source $D_{\scriptscriptstyle P}$  and receiver  $U_{\scriptscriptstyle P}$ wavefields. For only one frequency, PERM can be described by
\begin{eqnarray}
D_{\scriptscriptstyle P}(\boldsymbol\xi) = {\mathcal L}_0(\boldsymbol\xi,{\bf x-h})I({\bf x},{\bf h}),
\label{eq:perm_opD}
\end{eqnarray}
and
\begin{eqnarray}
U_{\scriptscriptstyle P}(\boldsymbol\xi) = {\mathcal L}_0(\boldsymbol\xi,{\bf x+h})I({\bf x},{\bf h}).
\label{eq:perm_opU}
\end{eqnarray}
%${\bf M}$ is the operator which injects the original image $I$ at every depth level and
The operator ${\mathcal L}$ is the recursive upward propagation of the wavefields from the subsurface (represented by the the coordinates $\bf x$) up to the surface (represented by the coordinates $\boldsymbol\xi$). This operator is also a function of the subsurface-offset because of the way the pre-stack image is injected into the modeling. The subscript of $\mathcal L$ denotes that the wavefield propagation is performed with the background slowness $s_0({\bf x})$ used to migrate the original shot-profiles. The source wavefield is upward propagated backwards in time and the receiver wavefield is upward propagated forward in time.
\par
The pre-stack image $I_{\scriptscriptstyle P}$ is computed with PERM data, at every depth level, using the multi-offset imaging condition
\begin{eqnarray}
I_{\scriptscriptstyle P}({\bf x},{\bf h}) = D^{*}_{\scriptscriptstyle P}({\bf x-h})U_{\scriptscriptstyle P}({\bf x+h}).
\label{eq:I_perm}
\end{eqnarray}
The wavefields in depth are recursively downward propagated according to 
\begin{eqnarray}
D^*_{\scriptscriptstyle P}({\bf x}) = {\mathcal L}'_1({\boldsymbol\xi},{\bf x}) D_{\scriptscriptstyle P}({\boldsymbol\xi}),
\label{eq:perm_opDz} 
\end{eqnarray}
and
\begin{eqnarray}
U_{\scriptscriptstyle P}({\bf x}) = {\mathcal L}'_1({\boldsymbol\xi},{\bf x}) U_{\scriptscriptstyle P}({\boldsymbol\xi}).
\label{eq:perm_opUz}
\end{eqnarray}
where $'$ stands for the complex conjugate. Note that in equations \ref{eq:perm_opDz} and \ref{eq:perm_opUz} the subscript of the operator ${\mathcal L}$ indicates the use of a different migration velocity. The source wavefield is downward propagated forward in time, and the receiver wavefield is downward propagated backwards in time. Inserting equation \ref{eq:perm_opD} into equation \ref{eq:perm_opDz} and equation \ref{eq:perm_opU} into equation \ref{eq:perm_opUz}, and substituting the expressions for the wavefields in depth into equation \ref{eq:I_perm} we get
\begin{eqnarray}
\begin{array}{lll}
I_{\scriptscriptstyle P}({\bf x},{\bf h}) & = &\left \{ {\mathcal L}'_0({\boldsymbol\xi},{\bf x-h}){\mathcal L}_1(\boldsymbol\xi,{\bf x-h})I({\bf x},{\bf h}) \right \} \\
&\times& \left \{ {\mathcal L}'_1({\boldsymbol\xi},{\bf x+h}){\mathcal L}_0(\boldsymbol\xi,{\bf x+h})I({\bf x},{\bf h}) \right \}.
\end{array}
\label{eq:I_perm1}
\end{eqnarray} 
From equation \ref{eq:I_perm1}, we can see that the pre-stack image $I_{\scriptscriptstyle P}$ is approximately a squared version of the original image. This means that in addition to the stronger side lobes, reflectors in the PERM image will always have positive polarity.
\par
For the simple case of a horizontal reflector in a constant velocity medium, we have shown that migration of PERM data produces images with the same kinematics as the shot profile migration. Now, we introduce a dipping reflector in the same constant background velocity medium. In the presence of a non-zero geological dip, a pre-processing of the initial conditions is necessary to obtain correct kinematics. This pre-processing step is represented by a rotation of the pre-stack image according to the apparent geological dip.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Dip-independent initial conditions}
Shot-profile and areal-shot migrations by wavefield extrapolation computes pre-stack images by means of the multi-offset imaging condition \cite[]{james2002}, in which source and receiver wavefields are laterally shifted prior to time correlation. However, the shift between wavefields might not be restricted to the horizontal direction. For instance, vertical shifts of the wavefields produce the vertical-subsurface offset gathers which provide reliable velocity information in the presence of steep dips \cite[]{biondi:1284}.
\par
Ideally, wavefields should be shifted along the geological dip direction. According to \cite{biondi:1283}, SODCIGs computed this way do not suffer from image-point dispersal in the presence of dip and migration velocity errors. The image-point dispersal makes events with different reflection angle from the same reflection point in the subsurface be imaged at different locations.
\par 
The image-point dispersal in 2D is illustrated in Figure \ref{fig:pdisp} for the case of migrating with a slower velocity (Figure \ref{fig:pdisp}a) and higher velocity (Figure \ref{fig:pdisp}b) than the true velocity. Let us consider constant velocity in the vicinity of the image point, so source and receiver rays are straight for simplicity.
\par
 When migration velocity is too low, the reflector is imaged at a smaller depth. The image point computed with horizontal shifts of the wavefields $\mathsf {I_{hx}}$ is shifted down-dip with respect to the image point computed with shifts along the apparent geological dip $\mathsf {I_{hg}}$. The geological dip is called apparent because of the migration velocity error. The point $\mathsf I$ is where source and receiver rays cross at an angle which is  twice the apparent reflection angle $\gamma$. Source and receiver rays cross deeper than the image points, causing events curve down in the SODCIG. 
\par
When migration velocity is too high, the reflector is imaged at a greater depth. The image point $\mathsf {I_{hx}}$ is shifted up-dip with respect to $\mathsf {I_{hg}}$. Source and receiver rays cross shallower than the image points, causing events curve up in the SODCIG. 
\activeplot{pdisp}{width=\textwidth}{}{Geometry for the computation of SODCIGs. Source, receiver and image points are labeled with $\mathsf S$, $\mathsf R$ and $\mathsf I$, respectively. The subscript $\mathsf {hx}$ corresponds to subsurface offsets computed with horizontal shift. The subscript $\mathsf {hg}$ corresponds to subsurface offsets computed by shifting along the apparent geological dip $\alpha$. a) too low velocity, and b) too high velocity. Modified from \cite{biondi:1283}.}
%\par
%Using simple trigonometry, from Figure \ref{fig:pdisp} we can compute the image-point dispersal for the in-line subsurface offset as
%\begin{eqnarray}
%\Delta {\bf I}_{hx} = {\bf h}_g \tan \gamma \tan \alpha,
%\label{eq:i_disp}
%\end{eqnarray}
%where ${\bf h}_g$ is the vector representing the geological-subsurface offset aligned along the apparent geological dip. The geological-subsurface offset, in turn, is related to the horizontal-subsurface offset $h_x$ by
%\begin{eqnarray}
%h_g = h_x \cos \alpha.
%\label{eq:h0_hx}
%\end{eqnarray}
%Notice that the image-point dispersal depends indirectly on migration velocity through ${\bf h}_g$, $\gamma$ and $\alpha$.
\par
Generating SODCIGs along the geological-dip direction overcomes the problem of the image-point dispersal. However, it is computationally demanding since wavefields must be stored at various depths and, in addition, accurate dip information is difficult to obtain, especially when events cross because of velocity inaccuracy.
\par
\cite{biondi:1283} point out that, at least in the first order, the reflection-angle domain is immune to image-point dispersal. This is because the SODCIG to ADCIG transformation shifts events to the line connecting $\mathsf {I}$ and $\mathsf {I_{hg}}$ in Figure \ref{fig:pdisp} at the same image point shared by all the reflection angles.
\par
In the presence of dip, to accurately model PERM data it is crucial that the initial conditions are free of image-point dispersal to guarantee that all the energy of a point in the subsurface is contained by the correspondent SODCIG injected to modeling. Since SODCIGs along the geological dip are not easily computed, can we pre-process the SODCIGs computed with horizontal shifts of the wavefields such that they are transformed into a good approximation of the SODCIGs along the geological dip?
\par
To answer this question, let us first examine the angle relationships in Figure \ref{fig:pdisp}. The angles $\gamma + \alpha$ and $\gamma - \alpha$ are the source and receiver ray angles, respectively. They are the propagation direction of the wavefields locally at the image point. Notice that $\alpha$ and $\gamma$ are related to slopes in the pre-stack image according to
\begin{eqnarray}
\tan \alpha = -\frac{dz_m}{dx_m}
\label{eq:talpha}
\end{eqnarray}
and
\begin{eqnarray}
\tan \gamma = -\frac{dz}{dh_x},
\label{eq:tgamma}
\end{eqnarray}
where the subscript $m$ refers to the local nature of the relationship. The solutions of the differential equations \ref{eq:talpha} and \ref{eq:tgamma} define slant-stack paths, which allow us to transform the 2D pre-stack image $I(x,z,h_x)$ into $I(x,z,\alpha,\gamma)$ by angle decomposition according to the following integrals:
\begin{eqnarray}
\left . I(x,z,\alpha,\gamma) = \int_{x_{m_{i}}}^{x_{m_{f}}} \int_{-h_x}^{h_x} W(x_m-x) \frac{dI(x,z,h_x)}{dz} dx_{m} dh_{x} \right |_{\substack{z=z_h+h_x\tan\gamma\\~z=z_m+x_m\tan\alpha}}
\label{eq:sltstk1}
\end{eqnarray}
where the derivative with respect to $z$ is to recover the correct phase. The  local window $W(x_m-x)$ is used in the local slant-stack integral on $x_m$, being defined as
\begin{eqnarray}
\left\{ \begin{array}{l}
1 , x_{m_{i}} \leq x \leq x_{m_{f}} \nonumber , \\
0,  elsewhere \end{array} \right .
\label{eq:window}
\end{eqnarray}
where $x_{m_{i}}={x_m-\frac{x_w}{2}}$ and $x_{m_{f}}={x_m+\frac{x_w}{2}}$, $x_w$ is the width of the local window.
\par
Again, using simple trigonometry, we have 
\begin{eqnarray}
\tan (\gamma + \alpha) = \frac{\tan \gamma + \tan \alpha}{1 - \tan \gamma \tan \alpha},
\label{eq:theta0a}
\end{eqnarray}
\begin{eqnarray}
\tan (\gamma - \alpha) = \frac{\tan \gamma - \tan \alpha}{1 + \tan \gamma \tan \alpha}.
\label{eq:theta0b}
\end{eqnarray}
\par
To align the initial conditions with the geological dip, we need to change the dip along the subsurface-offset axis according to the apparent geological dip, yielding the new subsurface offset $\widetilde{h_{x_s}}$ and $\widetilde{h_{x_r}}$ for the initial conditions of the modeling of source and receiver wavefield, respectively. This is accomplished by solving the following differential equations:
\begin{eqnarray}
\tan (\gamma + \alpha) = -\frac{dz}{d\widetilde {h_{x_s}}}, 
\label{eq:theta1}
\end{eqnarray}
\begin{eqnarray}
\tan (\gamma - \alpha) = \frac{dz}{d\widetilde {h_{x_r}}}.
\label{eq:theta2}
\end{eqnarray}
The solutions of equations \ref{eq:theta1} and \ref{eq:theta2} define new slant-stack operations which, in combination with equations \ref{eq:theta0a} and \ref{eq:theta0b}, reduces the dimensionality of the decomposed pre-stack image (equation \ref{eq:sltstk1}) by transforming $I(x,z,\alpha,\gamma)$ into $I_D(x,z,\widetilde{h_{x_s}})$ and $I_U(x,z,\widetilde{h_{x_r}})$.
\par
In 3D, we must also consider rotation of the cross-line offsets according to the apparent geological dip in the cross-line direction in addition to the in-line rotation. To apply the rotation in 3D, we need to assume that source and receiver rays are coplanar, such that they cross. In 3D, this assumption might not be valid when migrating with an inaccurate velocity in complex overburden. The co-planarity assumption is also present in the formulation of the common-azimuth migration \cite[]{biondi:1822}, which turned out to provide sufficiently accurate images in areas of complex geology.
\par
Instead of performing a true 3D slant-stack, which can be computer memory intensive, we can take advantage that a 3D slant-stack can be split into two 2D slant-stacks on the in-line and cross-line directions. This can be easily seen by using Fourier transforms to conveniently express a generic 3D slant-stack integral 
\begin{eqnarray}
\begin{array}{lll}
U(\phi_1,\phi_2,\phi_3) &=& \int \int u(z,x,y) dx~dy \left . \right |_{t=\phi_1+\phi_2 x+\phi_3 y} \\
 &=& \int \int u(\phi_1+\phi_2 x+\phi_3 y,x,y) dx~dy
\end{array}
\label{eq:sltstk1}
\end{eqnarray}
as
\begin{eqnarray}
U(\phi_1,\phi_2,\phi_2)=\int d\phi_1 e^{-ik\phi_1} \int dx e^{-i\phi_2 x}\int dy e^{-i\phi_3 y}\int dz e^{ikz}u(z,x,y).
\label{eq:sltstk2}
\end{eqnarray}
where $\phi_1$, $\phi_2$ and $\phi_3$ are, respectively, the output $z$, the slanting parameter corresponding to $x$, and the slanting parameter corresponding to $y$.
\par
To illustrate the generation of dip-independent initial conditions 801 split-spread shots 10 m apart with maximum offset of 3250 m were modeled with velocity of 1000 m/s and migrated with a $10\%$ slower velocity (Figure \ref{fig:dip01}). The model has a $20^{\circ}$ dipping reflector and a horizontal reflector at a depth of 2500 m. The SODCIG located at 0 m was used as the initial condition for the modeling of PERM data without applying the pre-processing described above. Since PERM models one event of one isolated SODCIG, the dipping reflector and the horizontal reflector originate two different pairs of PERM wavefields. This means that reflectors used in the modeling need to be interpreted in the pre-stack volume. The pair of source and receiver wavefields for the dipping reflector are shown in Figure \ref{fig:dip02}.
\activeplot{dip01}{width=\textwidth}{}{Shot-profile migration of 801 split-spread shots 10 m apart with velocity $10\%$ slower than the true velocity. The model is represented by a $20^{\circ}$ dipping reflector and a horizontal reflector at depth of 2500 m embedded in constant velocity of 1000 m/s.}
\activeplot{dip02}{width=\textwidth}{}{Data synthesized by PERM having as the initial condition the dipping reflector in the SODCIG at $x_m=0$ m. a) The receiver wavefield. b) The source wavefield.}
\par
Since the modeling of a single non-rotated SODCIG carries no dip information, migration of the corresponding PERM data using the correct velocity does not shift events laterally as can be seen in Figure \ref{fig:dip03}. As expected, the horizontal reflector focuses at zero-subsurface offset. However, notice how the dipping reflector still presents a residual curvature. Migration of PERM data from SODCIGs within a neighborhood around $x_m = 0$ m is shown in Figure \ref{fig:dip04}. Again, the residual curvature is present and any migration velocity analysis using this result will lead to wrong velocity updates. This residual curvature is a result of not having corrected the image-point dispersal. Unless stated, SODCIGs in the figures are selected at $x = 0$ m.
\par
\activeplot{dip03}{width=\textwidth}{}{Areal-shot migration of PERM data shown in Figure \ref{fig:dip02} using the correct velocity. The horizontal reflector is focused at zero-subsurface offset, but the dipping reflector shows residual curvature.}
\activeplot{dip04}{width=\textwidth}{}{Areal-shot migration with correct velocity of PERM data having a set of isolated SODCIGs around $x_m = 0$ m as the initial conditions. As in Figure \ref{fig:dip03}, the horizontal reflector is focused at zero-subsurface offset, but the dipping reflector shows residual curvature.}
\par
The rotation was applied to the image in Figure \ref{fig:dip01}, and the new PERM data was modeled using the initial conditions shown in Figure \ref{fig:dip05}. Notice that the initial condition for modeling the source wavefield (Figure \ref{fig:dip05}a) and the initial condition for modeling the receiver wavefield (Figure \ref{fig:dip05}b) have the dipping event oriented in opposite directions in the SODCIG. The rotation changes neither the horizontal reflector nor the-zero subsurface offset, as can be seen on the right panels.
\activeplot{dip05}{width=\textwidth}{}{Initial conditions for modeling a) source and b) receiver wavefields. The dipping reflector is oriented in opposite directions in the SODCIG. Rotation does no affect the horizontal reflector nor the-zero subsurface offset, as can be seen on the right panels.}
\par
 The source and receiver wavefields for the dipping reflector after rotation are shown in Figure \ref{fig:dip06}. The events in Figures \ref{fig:dip06}a and \ref{fig:dip06}b are shown in the same areal shot for illustration only. Actually, they pertain to different areal shots.
\activeplot{dip06}{width=\textwidth}{}{Dip-independent PERM data for the dipping reflector from the rotated SODCIG at $x_m = 0$ m. a) The receiver wavefield. b) The source wavefield.}
\par
Areal-shot migration of dip-independent PERM data is shown in Figure \ref{fig:dip07}. Notice that the segment of the dipping reflector is shifted laterally with respect to that of the horizontal reflector. Since the dip-independent wavefields carry information about the dip of the reflector, the observed reflector movement is now consistent with migration with a higher velocity. 
\activeplot{dip07}{width=\textwidth}{}{Areal-shot migration with the correct velocity of dip-independent PERM data having the rotated the SODCIGs at $x_m = 0$ m as the initial condition. Compare with Figure \ref{fig:dip03}. The dipping reflector is now focused, in contrast to the image in Figure \ref{fig:dip03}, where it shows residual curvature.}
\par
Migration with the correct velocity of dip-independent PERM data modeled from a set of SODCIGs in a neighborhood around $x_m = 0$ m confirms the correctness of the rotation (Figure \ref{fig:dip08}). The focusing of the dipping reflector  around zero-subsurface offset is greatly improved when compared with Figure \ref{fig:dip04}. The corresponding ADCIGs confirm the more consistent move-out after rotation (Figure \ref{fig:dip09}). Note the residual move-out in the angle gather corresponding to the image computed with wavefields with non-rotated initial conditions (Figure \ref{fig:dip09}a) and how much flatter is the image from wavefields computed with the proposed rotation (Figure \ref{fig:dip09}b).
\activeplot{dip08}{width=\textwidth}{}{Areal-shot migration  with correct velocity of dip-independent PERM data having a set of rotated SODCIGs around $x_m=0$ m as the initial conditions. Compare with Figure \ref{fig:dip04}. The focusing of the dipping reflector is greatly improved when using the rotated initial conditions.}  
\activeplot{dip09}{width=0.66\textwidth}{}{ADCIGs of images computed with the correct migration velocity using PERM data having: a) non-rotated initial conditions, and b) rotated initial conditions. Note the residual move-out in a) and the flatter response in b).} 
\par
In the example with two reflectors, we initiated each modeling experiment from one isolated SODCIG of one single reflector to compute a pre-stack image restricted to a certain region in the output space. However, depending on the number of reflectors and the size of the prestack image, this procedure can generate a dataset even larger than the original shot-profiles, defeating the original purpose of PERM, which is to synthesize a smaller dataset to be used in migration velocity analysis. We see next, that combination of modeling experiments can be used to decrease the size of PERM data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Combination of modeling experiments}
In the previous examples, if we were to fully image the reflectors, the number of areal shots might be twice as many as the original shot-profiles. To decrease the number of modeling experiments, linearity of wave propagation can be used to combine isolated SODCIGs and inject them simultaneously into one single model experiment as
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}-i\sqrt{\omega ^2 s_0^2({\bf x})-|{\bf k}| ^2}\right) \widehat {D_{\scriptscriptstyle P}}({\bf x},\omega;\widehat{{\bf x}}) = \widehat {I_D}(\widehat{{\bf x}},{\bf h}) \\
\widehat{D_{\scriptscriptstyle P}}(x,y,z=z_{\rm max},\omega;\widehat{{\bf x}}) = 0 \end{array} \right.,
\label{eq:comb_bD}
\end{eqnarray}
and
\begin{eqnarray}
\left\{ \begin{array}{l}
\left( \frac{\partial}{\partial z}+i\sqrt{\omega ^2 s_0^2({\bf x})-|{\bf k}| ^2}\right) \widehat{U_{\scriptscriptstyle P}}({\bf x},\omega;\widehat{{\bf x}}) = \widehat {I_U}(\widehat{{\bf x}},{\bf h}) \\
\widehat{U_{\scriptscriptstyle P}}(x,y,z=z_{\rm max},\omega;\widehat{{\bf x}}) = 0 \end{array} \right.,
\label{eq:comb_bU}
\end{eqnarray}
where $\widehat{I_D}(\widehat{{\bf x}},{\bf h})$ and $\widehat{I_U}(\widehat{{\bf x}},{\bf h})$ are the combination of SODCIGs for a single reflector to be used as the initial conditions for the modeling of combined wavefields, $\widehat{D_{\scriptscriptstyle P}}({\bf x},\omega;\widehat{{\bf x}})$ and $\widehat{U_{\scriptscriptstyle P}}({\bf x},\omega;\widehat{{\bf x}})$, the source and receiver wavefields, respectively. Notice that $\widehat {\bf x}$ is a series of evenly spaced points on the reflector, where the SODCIGs are selected. In 3D, the selection of SODCIGs can be thought of as the multiplication of the pre-stack image by a 2D comb function. This sampling function is shifted laterally to select new set of SODCIGs to initiate the modeling of another pair of combined wavefields. After shifting along one period of the sampling function in the $x$ and $y$ directions, all the points on the reflector are used in the modeling. Consequently, the number of modeling experiments equals the number of lateral shifts of the sampling function.
\par
The choice of the sampling period determines how much crosstalk will be present in the migrated image, in a manner similar to the shot-profile migration illustrated in the introduction of this chapter. To obtain a crosstalk-free image, the sampling period needs to be set to a distance such that wavefields initiated at different SODCIGs do not correlate. As previously shown when discussing how to compute PERM images with kinematics similar to that of the original shot-profile migration, PERM wavefields generated from SODCIGs within the interval $2(max\left | h_x \right |,max\left | h_y \right |)$ still contribute to the image at the central SODCIG. For the same reason, no crosstalk is generated during migration if the period of the sampling function is larger than that interval. Since the focusing of energy in the SODCIG is velocity error dependent, so is the the period of the sampling function. Therefore, for small velocity errors a small sampling period can be used and, consequently, a smaller number of combined modeling experiments is needed.
\par
The combination of modeling experiments is illustrated in Figure \ref{fig:comb01}. We model PERM data starting with the rotated images from the previous section and combine SODCIGs into sets using the sampling period of 41, 81 and 163 SODCIGs. As can be seen in Figures \ref{fig:comb01}a and \ref{fig:comb01}b, crosstalk occurs in the SODCIG according to a period which is half of the sampling period. Recalling that the number of subsurface-offsets is 81, it is expected that no crosstalk will occur when migrating the set with sampling period of 163 SODCIGs (Figure \ref{fig:comb01}c). In this case, since each reflector is used separately in the modeling, the total number of areal shots is 326, which is less than half of the original shot profiles. 
\activeplot{comb01}{width=0.85\textwidth}{}{Areal-shot migration of PERM data synthesized from sets of SODCIGs selected with sampling period of: a) 41, b) 81, and c) 163 SODCIGs. Notice that no crosstalk is generated when the sampling period is larger than twice the maximum absolute subsurface-offset.}
\par
To further decrease the number of modeling experiments, we can combine the reflectors in the initial condition. By doing so, two events are present both in the source wavefield and in the receiver wavefield, as in Figure \ref{fig:dip06}. During migration, in addition to the expected cross-correlations, cross-correlation of wavefields from different reflectors will occur, generating crosstalk between reflectors. Areal-shot migration of PERM data synthesized from a set of SODCIGs with sampling period of 163 is shown in Figure \ref{fig:comb02}. Both reflectors are simultaneously injected in the initial condition. Notice the presence of reflector crosstalk labeled with `C'.
\activeplot{comb02}{width=0.85\textwidth}{}{Areal-shot migration of PERM data synthesized from a set of SODCIGs selected with sampling period of 163. The two reflectors are simultaneously injected to the model. Notice the reflector crosstalk resulting from the cross-correlation of the wavefields from the horizontal reflector with  that from the dipping reflector.}
\par
To illustrate the reflector crosstalk generation, snapshots of the wavefield propagation are shown in (Figure \ref{fig:comb03}). Both wavefields are represented in the same panel and signaled with `D' (from downgoing) for the source wavefield and `U' (from upgoing) for the receiver wavefield. The panels on the left show the wavefields in the time-depth domain selected at horizontal positions where reflector crosstalk occurs in Figure \ref{fig:comb02}. The panels on the right are taken at the propagation time when the wavefields cross on the left panels. The crossing times vary from 0.14 s to 1.24 s, for the SODCIGs from right to left, respectively. The crossing times are a function of the distance between reflectors and the propagation velocity. In Figure \ref{fig:comb03}a notice that for propagation times less than 0.14 s minus the period of the wavelet in time no crosstalk will occur.
%For the simple case of constant velocity, the time at which reflector crosstalk occurs is  
\activeplot{comb03}{width=0.80\textwidth}{}{Snapshots of the wavefield propagation of the Figure \ref{fig:comb02}. The panels on the left are selected at a horizontal position where the crosstalk occurs in Figure \ref{fig:comb02}. The panels on the right are taken at the propagation time when the wavefields cross on the left panel.}
\par
Crosstalk can be mitigated if we take advantage of the imaging principle, which states that the reflector image is formed at zero time of wavefield propagation. This is especially true when using migration velocity is accurate. When the migration velocity is inaccurate, the focusing of the image departs from the zero time, and this can potentially be used to update the velocity model \cite[]{sava:S209,yang:3715}. As Figure \ref{fig:comb03} shows, crosstalk is formed at times different from zero. Therefore, if the wavefields are cross-correlated within a time window centered at zero, reflector crosstalk might not occur depending on the distance in depth between the reflectors \cite[]{biondi07}. The time-windowed imaging condition for a single areal shot reads
\begin{eqnarray}
I{\scriptscriptstyle {P_w}}({\bf x},{\bf h}) = \sum_{\frac{-t_{\scriptscriptstyle w}}{2} \leq t \leq \frac{t_{\scriptscriptstyle w}}{2}} {\mathcal F}^{-1} \left [ D^{*}({\bf x}-{\bf h}) \right ] {\mathcal F}^{-1}\left [ U({\bf x}+{\bf h}) \right ],
\label{eq:wic_bI}
\end{eqnarray}
where $t_{\scriptscriptstyle w}$ is the length of the time window. When using one-way propagators, the wavefields are inverse Fourier transformed to time by ${\mathcal F}^{-1}$. Migration using the time-windowed imaging condition of  equation \ref{eq:wic_bI}  is shown in Figure \ref{fig:comb04}. The length of the time window is 0.06 s which avoids the cross-correlation of events from different reflectors. The image is completely free of reflector crosstalk. 
\activeplot{comb04}{width=0.9\textwidth}{}{Areal-shot migration using the time-windowed imaging condition (equation \ref{eq:wic_bI}). The reflector crosstalk is completely eliminated.}
\par
We saw that by carefully combining the modeling experiments the data size is decreased while the correct kinematics is maintained, which is important for migration velocity analysis. However, combining the modeling experiments using a decorrelation distance between events does not achieve a significant data reduction. In this case, data reduction depends on the number of subsurface-offsets which are necessary to capture all the relevant velocity information. In the example using the time-windowed imaging condition, the number of independent experiments was reduced five times if compared with the original dataset. Data reduction techniques like plane-wave decomposition, for instance, could lower that number ten times. As we will see in Chapter \ref{chap:chap03}, further data reduction can be achieved by using the phase-encoding technique \cite[]{romero} to linearly combine the modeling experiments. This will enable us to use a shorter sampling period of SODCIGs, and also injecting more than one reflector in the modeling.
\par
Although PERM theory was developed in 3D, all the examples so far were 2D. Next, it is discussed an 3D example under the common-azimuth approximation \cite[]{biondi:1822}. It will be shown that in this case the SODCIGs in the $y$ direction can be continuously sampled and the number of modeling experiments will be determined only by sampling period in the $x$ direction, drastically decreasing data size.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{3D-{PERM} from common-azimuth migrated images}
In the way PERM is formulated there is no restriction on the dimensionality of the pre-stack image used as the initial condition for the modeling, which means that if the original data have sufficient cross-line offsets as in the acquisition geometries with wide-range of azimuths \cite[]{regone:SM231,kapoor:926,moldoveanu:908}, the initial conditions are a five-dimensional hypercube on $\bf x$, $h_x$ and $h_y$. 
\par
To synthesize PERM data starting with the five-dimensional initial conditions such that no crosstalk is generated during migration, the minimum number of modeling experiments is $4n_{h_{x}}n_{h_{y}}$, where $n_{h_{x}}$ and $n_{h_{y}}$ are the number of subsurface offsets in $x$ and $y$ directions. Considering usual parameters, the number of modeling experiments would be as low as some hundreds or few thousands. This data reduction is very substantial if we compare, for instance, with data reduction achieved by 3D-plane-wave migration. Using plane waves, to obtain artifact-free SODCIGs due to the lack of illumination from some propagation directions we need to migrate roughly 2000 plane waves. This means that 3D-PERM data size can be one order of magnitude smaller than 3D-plane wave data.
\par
Despite the recent good migration results obtained in geologically complex areas using wide-azimuth data, narrow-azimuth acquisition is still the industry standard. Narrow-azimuth data can be efficiently imaged by common-azimuth wave-equation migration (CAM) \cite[]{biondi:1822}. That is because  the dimensionality of the pre-stack wavefields is decreased by assuming zero-cross-line offset in contrast with the full-azimuth migration. That does not mean that the cross-line offset wavenumber is zero. Rather, its asymptotic approximation is function of the the in-line midpoint and in-line offset wavenumbers. Therefore, instead of a five-dimensional hypercube CAM images are  four-dimensional hypercubes on $\bf x$ and $h_x$. 
\par
Because of the lower dimensionality of CAM images, when using them as the initial conditions to synthesize PERM data, the SODCIGs in the cross-line direction can be sampled continuously as depicted in Figure \ref{fig:cam01}b in contrast with the five-dimensional initial conditions for the full azimuth case of Figure \ref{fig:cam01}a. Recall that PERM is equivalent to ERM if energy is focused at the zero-subsurface offset as well as if this is the only available subsurface offset. 
\par
The continuous sampling of SODCIGs in the cross-line direction yields one more order of magnitude of data reduction. So, under the common-azimuth approximation, 3D-PERM data size can be two orders of magnitude smaller than 3D-plane wave data.
\activeplot{cam01}{width=1.1\textwidth}{}{The initial conditions for synthesizing PERM data from CAM images can be specified as in a) because no pre-stack information exists in the cross-line direction, in contrast with the full azimuth situation in b).}
\par
To illustrate the validity of the above assumptions, a split-spread data with maximum offset of 1587.5 m was computed using 3D-Born modeling \cite[]{james96} on a 30$^{\circ}$ dipping reflector with 45$^{\circ}$ azimuth with respect to the acquisition direction, which is aligned with the in-line direction. The number of in-lines and cross-lines is 96 each spaced at every 25 m. The offset interval is 25m.  The velocity used in the modeling is the 1D function $v(z) = (1500 + 0.5z)$ m/s.
\par
The Born data was input to CAM with a 5$\%$ slower velocity. Migration results can be seen in Figures \ref{fig:cam02}a and \ref{fig:cam02}b for SODCIGs positioned at $(x = 750$ m,$y = 600$ m$)$ and  $(x = 750$ m,$y = 1000$ m$)$, respectively. The panel on the left is the SODCIG which contains 21 subsurface-offsets ranging from $-250$ m to $250$ m. The panel in the middle is the in-line at the zero-subsurface offset, and $y = 600$ m (Figure \ref{fig:cam02}a) and $y = 1000$ m (Figure \ref{fig:cam02}b). The panel on the right is the cross-line  at the zero-subsurface offset, and $x = 750$ m.
\activeplot{cam02}{width=0.9\textwidth}{}{Common-azimuth migration of 3D-Born data modeled from a 30$^{\circ}$ dipping reflector with 45$^{\circ}$ azimuth with respect to the acquisition direction. The panel in the middle is the in-line at the zero-subsurface offset, and $y = 600$ m (Figure \ref{fig:cam02}a) and $y = 1000$ m (Figure \ref{fig:cam02}b). The panel on the right is the cross-line  at the zero-subsurface offset, and $x = 750$ m.}
\par
In the common-azimuth regime, the computation of the dip-independent initial conditions is performed by simply rotating the SODCIGs  in the in-line direction, since no cross-line offset is computed in migration. The underlying assumption that source and receiver rays cross, which makes the rotation valid in 3D, implies that they are coplanar. This assumption is in agreement with the co-planarity condition present in the formulation of CAM.
\par
Source and receiver wavefields were modeled using continuous sampling along the cross-line direction and sampling period of 48 in the in-line direction. This period is sufficient to not generate crosstalk during the areal-shot migration, given that the number of subsurface-offsets of the pre-stack image is 21. The synthesized 3D receiver wavefield is shown in Figure \ref{fig:cam03}. The left panel is the in-line at $y = 1200$ m, the right panel is the cross-line at $x = 1400$ m, and the top panel is the time-slice at $t = 0.5$ s.
\activeplot{cam03}{width=0.9\textwidth}{}{3D-PERM receiver wavefield. The left panel is the in-line at $y = 1200$ m, the right panel is the cross-line at $x = 1400$ m, and the top panel is the time-slice at $t = 0.5$ s.}
\par
The 3D migration of the 48 areal shots with the 5$\%$ slower velocity is shown in Figures \ref{fig:cam04}a and and \ref{fig:cam04}b for SODCIGs positioned at $(x = 750$ m,$y = 600$ m$)$ and  $(x = 750$ m,$y = 1000$ m$)$, respectively. For facilitating the comparison with the CAM images of Figures \ref{fig:cam02}a and \ref{fig:cam02}b, the polarity of the areal-shot migrated image is inverted due to the squaring of the wavelet. The kinematics of the SODCIGs computed with PERM wavefields matches that of the SODCIGs computed with CAM, especially between $-200$ m and $200$ m. This enables the use of 3D PERM wavefields computed from CAM images in optimization of migration velocity, as will be shown in Chapter \ref{chap:chap06} for a 3D survey from the North-Sea.
\activeplot{cam04}{width=0.9\textwidth}{}{3D-areal-shot migration of PERM data. The panel in the middle is the in-line at the zero-subsurface offset, and $y = 600$ m (Figure \ref{fig:cam02}a) and $y = 1000$ m (Figure \ref{fig:cam02}b). The panel on the right is the cross-line  at the zero-subsurface offset, and $x = 750$ m. Compare with Figure \ref{fig:cam02}.}
\newpage 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
\section{Conclusions}
In this chapter we saw that wavefields synthesized by PERM are able to provide migrated images with correct kinematics while decreasing data size. Data reduction is achieved by combining the modeling experiments and is controlled by the number of subsurface offsets which will be computed during areal-shot migration of PERM data. Recall that SODCIGs in the initial conditions must be separated by a distance at least twice the maximum absolute subsurface-offset value to prevent crosstalk. Implicit to PERM is the fact that reflectors must be identified such that reflector crosstalk is avoided during migration. 3D Pre-stack interpretation can be cumbersome, but it allows, for instance, avoiding the use of reflectors with low signal-to-noise ratio in migration velocity estimation. Moreover, in commercial softwares for migration velocity estimation reflector picking is a usual and almost entirely automated procedure.
\par
While in 2D data reduction with PERM is comparable to that of the plane-wave decomposition, in 3D it can be of one order of magnitude less if cross-line subsurface-offsets are to be computed. If the initial conditions are computed with common-azimuth migration, PERM data size can potentially be two orders of magnitude less than 3D plane-wave data.
%
\par
\cleardoublepage
%\bibliographystyle{seg}
%\bibliography{SEP,SEG,MISC,GEOTLE,allrefs}

