%% This file is automatically generated. Do not edit!
\documentclass[12pt,twoside,color]{geophysics}

\usepackage{amsmath,color}

\setfigdir{Fig}

\begin{document}
\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. PERM shares with ERM the ability to perform multiple modeling experiments simultaneously, greatly reducing data volume. 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, including 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, Because of this assumption, ERM does not model multiples and prismatic waves. Moreover, the coincident ray-path assumption is often invalid in areas of geological complexity. 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}. Instead Kirchoff migration and ray-based methods are still the industry standards. In addition to the lower cost, ray-based methods are very flexible with respect to strategies for defining the velocity model using selected reflectors \cite[]{stork:680,kosloff:1511,kosloff:1815,billette:1889}, which is not possible using the conventional migration-velocity analysis with wave-extrapolation methods. 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. In addition, it can be usefull to incorporate ray-based strategies, as the horizon-based or the grid-tomography approaches, in the framework of migration-velocity analysis by wavefield extrapolation.
\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}, in which the wavefields are spatially cross-correlated. The lags of the spatial cross-correlation are the so-called subsurface offsets. In Figure \ref{fig:perm01}a is the migration result of  380 shots 24 m apart, and in Figure \ref{fig:perm01}b is the migration result of 23 shots 384 m apart. 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. For the full migration using all the shots, the slanted lines constructively interfere and the energy is focused at zero-subsurface-offset (Figure \ref{fig:perm01}a). 
\plot{perm01}{width=0.8\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. This can be clearly seen in the angle-domain common-image gathers (top) and residual moveout panels (bottom) in Figure \ref{fig:perm03}a-c. They are selected at the same horizontal position as the SODCIG of Figures \ref{fig:perm01} and \ref{fig:perm02}. The left column shows angle gather and residual moveout panel for the full migration, the column in the middle corresponds to the migration of 23 sampled shots, and the right column corresponds to the migration of the combined shots. Whereas the sub-sampling of shots still yields a reasonable residual moveout information, the crosstalk resulting from the combination of many shots destroys the velocity information.
\plot{perm02}{width=0.8\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.}
\par
\plot{perm03}{width=0.8\textwidth}{Angle-domain common-image gathers (top) and residual moveout panels (bottom). Left: full migration. Middle: migration of 23 sampled shots. Right: migration of the combined shots. Crosstalk resulting from the combination of many shots destroys the velocity information.}
\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; rather 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 generate 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 the 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 essential 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 used in different methods, such as controlled illumination \cite[]{Rietveld01,Rietveld02} and common-focus point (CFP) \cite[]{berkhout:937, berkhout:954, thorbecke:S273}.
\par
%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. I will show 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 theoretically straightforward but challenging in practice because of computational cost and data handling issues. I will present the general 3D theory and show 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 the 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 strengths proportional to their reflection coefficients. 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 the 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 zero subsurface offset and that it is sufficient to parameterize the migrated image as a function of only 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
The generalization of ERM proposed in this thesis does not aim at computing pre-stack data $\it strictu~sensu$. Rather, it synthesizes data in the generalized source domain 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 the initial condition to propagate wavefields. The modeled data are, 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} 
Migration-velocity analysis using ray-based methods is very flexible with respect of strategies we can use to optimize the velocity model. For instance, regarding to the intermediate migration of the nonlinear iterations, Kirchhoff migration is the most commonly used. However, images obtained with wave-extrapolation can also be used \cite[]{bob00}, which yields more reliable moveout information in areas of geological complexity. 
\par
Also, in ray-based migration-velocity analysis reflectors must be identified. That is because traveltime perturbations can be caused by both velocity perturbations or perturbations in the reflector position. Two distinct approaches arise from the identification of reflectors: the horizon-based tomography \cite[]{kosloff:1511} and the grid tomography \cite[]{kosloff:1815,billette:1889}. The horizon-based tomography uses representative reflectors, which typically have good signal-to-noise ratio and characterize the main velocity changes, to define the macro-velocity structure. Grid tomography also selects coherent reflectors but without the need of being laterally continuous. Contrasting with the horizon-based approach, grid tomography uses several reflector segments to estimate the velocity perturbation. At the later stages of migration-velocity analysis, grid tomography can also be used to refine a velocity model obtained with horizon-based tomography.
\par
Migration-velocity analysis using wavefield-extrapolation methods does not have the same flexibility of its ray-based counterpart. This is partly because ray-based methods are a more mature technology. A characteristic of wave-extrapolation methods is the use of all events in the recorded wavefield to compute velocity perturbations. Ideally, the events should be consistent with the physics of wavefield-extrapolation method. When this is not the case, a preprocessing step is necessary to attenuate the events or the attenuation must be incorporated into the objective function  \cite[]{mulder:1184}. These authors show how multiples can strongly bias the velocity updates. Another strategy is to select only representative and coherent events for migration-velocity analysis using wavefield-extrapolation methods. However, differently from the ray-based migration-velocity analysis in which events are easily identified in the model space, event identification for migration-velocity analysis based on wavefield-extrapolation should be done in the data domain, which is a very difficult task.
\par
The use of PERM wavefields can enable us to incorporate ray-based strategies in wavefield-extrapolation methods. As we will see in this chapter, reflectors must be identified to model PERM wavefields. Hence, a horizon-based strategy is easily incorporated in migration-velocity analysis by wavefield-extrapolation. Moreover, due to the localized nature of the initial conditions for modeling PERM wavefields, migration of a single pair of PERM wavefields produces an image with a fairly good approximation of the correct kinematics. Therefore, using PERM wavefields potentially enable us to apply a grid-based strategy in migration-velocity analysis by wavefield-extrapolation.
\par
In addition to allowing more flexibility in migration-velocity analysis by wavefield-extrapolation, the size of PERM data set can be orders of magnitude smaller than the original data set enabling faster velocity updates.
\par
The fundamental idea of PERM is to model data that 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 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 PERM source $D_{\scriptscriptstyle P}$ and receiver $U_{\scriptscriptstyle P}$ 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)$. We do not consider 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; and 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 a 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.

\plot{refpl01}{width=0.8\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.}

\plot{refpl02}{width=0.8\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 those of the original shot-profile migration (compare Figures \ref{fig:refpl03} and \ref{fig:refpl01}). However, for farther subsurface offsets energy is not adequately imaged. This can be easily explained by analyzing how the image $I_{\scriptscriptstyle P}({\bf x},{\bf h};{\bf x_m})$ is formed when applying the multi-offset imaging condition to the downward propagated PERM wavefields modeled from a SODCIG located at $\bf x_m$. For one particular frequency, it reads
\begin{eqnarray}
I_{\scriptscriptstyle P}({\bf x},{\bf h};{\bf x_m}) = D^{*}_{\scriptscriptstyle P}({\bf x-h};{\bf x_m})U_{\scriptscriptstyle P}({\bf x+h};{\bf x_m}),
\label{eq:I_perm}
\end{eqnarray}
where `$^{*}$' stands for the complex conjugate. Notice that the maximum absolute distance at which wavefields still correlate is $\left |2h_{max} \right |$ for symmetric SODCIGs with respect to subsurface offset or, more generally, twice the subsurface offset range. Therefore, to ensure that the areal-shot migrated image has kinematics at all subsurface offsets similar to those in the original isolated SODCIG, we need to model PERM data from a set of SODCIGs within that neighborhood around the central SODCIG.
\plot{refpl03}{width=0.8\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. 

\plot{refpl04}{width=0.8\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.} 
\plot{refpl05}{width=0.8\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 are 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.

\plot{refpl06}{width=0.8\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 and the image computed with PERM wavefields resembles the original shot-profile migrated image. However, by carefully examining these images we can see that the former has stronger side lobes due to the squaring of the wavelet, which is mathematically explained by the modeling and migration of PERM wavefields initiated at a single SODCIG selected at $\bf x_m$. For one particular frequency, and considering, 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, PERM modeling can be described by
\begin{eqnarray}
D_{\scriptscriptstyle P}(\boldsymbol\xi;{\bf x_m}) = \sum_{\bf x}\sum_{\bf h}G_0(\boldsymbol\xi,{\bf x-h})I({\bf x_m,h}),
\label{eq:perm_opD}
\end{eqnarray}
and
\begin{eqnarray}
U_{\scriptscriptstyle P}(\boldsymbol\xi;{\bf x_m}) = \sum_{\bf x}\sum_{\bf h}G_0(\boldsymbol\xi,{\bf x+h})I({\bf x_m,h}).
\label{eq:perm_opU}
\end{eqnarray}
The pre-stack image is injected to the modeling by projecting the subsurface offsets $\bf h$ on the spatial axis $\bf x$. The Green's function $G_0$ upward propagates the wavefields from the subsurface (represented by the the coordinates $\bf x$) up to the depth level where the wavefields are collected (represented by the coordinates $\boldsymbol\xi$). Its subscript 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 wavefields are recursively downward propagated  in depth according to 
\begin{eqnarray}
D_{\scriptscriptstyle P}({\bf x};{\bf x_m}) = \sum_{\boldsymbol\xi}G^*_1({\boldsymbol\xi},{\bf x}) D_{\scriptscriptstyle P}({\boldsymbol\xi};{\bf x_m}),
\label{eq:perm_opDz} 
\end{eqnarray}
and
\begin{eqnarray}
U_{\scriptscriptstyle P}({\bf x};{\bf x_m}) = \sum_{\boldsymbol\xi}G^*_1({\boldsymbol\xi},{\bf x}) U_{\scriptscriptstyle P}({\boldsymbol\xi};{\bf x_m}).
\label{eq:perm_opUz}
\end{eqnarray}
Note that in equations \ref{eq:perm_opDz} and \ref{eq:perm_opUz} the subscript of the Green's function 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. 
\par
The lateral shifts of the wavefields for the multi-offset imaging condition are represented by
\begin{eqnarray}
D_{\scriptscriptstyle P}({\bf x-h};{\bf x_m}) = \sum_{\boldsymbol\xi}G^*_1({\boldsymbol\xi},{\bf x-h}) \sum_{\bf x'}\sum_{\bf h'}G_0(\boldsymbol\xi,{\bf x'-h'})I({\bf x_m,h'}),
\label{eq:perm_lshD} 
\end{eqnarray}
and
\begin{eqnarray}
U_{\scriptscriptstyle P}({\bf x+h};{\bf x_m}) = \sum_{\boldsymbol\xi}G^*_1({\boldsymbol\xi},{\bf x+h}) \sum_{\bf x'}\sum_{\bf h'}G_0(\boldsymbol\xi,{\bf x'+h'})I({\bf x_m,h'}).
\label{eq:perm_lshU}
\end{eqnarray}
\par
The PERM image is obtained by inserting equations \ref{eq:perm_lshD} and  \ref{eq:perm_lshU} into equation \ref{eq:I_perm}
\begin{eqnarray}
\begin{array}{lll}
I_{\scriptscriptstyle P}({\bf x},{\bf h};{\bf x_m}) & = & \displaystyle{\sum_{\boldsymbol\xi}\sum_{\bf x'}\sum_{\bf h'}}G_0(\boldsymbol\xi,{\bf x'-h'}) G^*_1({\boldsymbol\xi},{\bf x-h}) I({\bf x_m,h'})  \\
&\times& \displaystyle{\sum_{\boldsymbol\xi}\sum_{\bf x'}\sum_{\bf h'}}G^*_1({\boldsymbol\xi},{\bf x+h}) G_0(\boldsymbol\xi,{\bf x'+h'})I({\bf x_m,h'}),
\end{array}
\label{eq:I_perm1}
\end{eqnarray} 
which finally gives
\begin{eqnarray}
\begin{array}{lll}
I_{\scriptscriptstyle P}({\bf x},{\bf h};{\bf x_m}) & = & \displaystyle{\sum_{\boldsymbol\xi'}\sum_{\bf x'}\sum_{\bf h'}}\displaystyle{\sum_{\boldsymbol\xi''}\sum_{\bf x''}\sum_{\bf h''}}G_0(\boldsymbol\xi',{\bf x'-h'}) G^*_1({\boldsymbol\xi'},{\bf x-h}) \\
&\times& G^*_1({\boldsymbol\xi''},{\bf x+h}) G_0(\boldsymbol\xi'',{\bf x''+h''})I({\bf x_m,h'})I({\bf x_m,h''}),
\end{array}
\label{eq:I_perm2}
\end{eqnarray} 
\par  
From equation \ref{eq:I_perm2}, we can see that, if the velocity used in the downward propagation is the same as that in the upward propagation, the pre-stack image $I_{\scriptscriptstyle P}$ at $x_m$ 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 and the amplitude variation in the PERM image will be more pronounced than that in the original image. The last amplitude effect can be mitigated by taking the square root of the absolute value of initial conditions while keeping is polarity.
\par
Migration of PERM wavefields generated at all SODCIGs is obtained by summing over all $\bf x_m$
\begin{eqnarray}
\begin{array}{lll}
I_{\scriptscriptstyle P}({\bf x},{\bf h}) & = & \displaystyle{\sum_{\bf x_m}\sum_{\boldsymbol\xi'}\sum_{\bf x'}\sum_{\bf h'}}\displaystyle{\sum_{\boldsymbol\xi''}\sum_{\bf x''}\sum_{\bf h''}}G_0(\boldsymbol\xi',{\bf x'-h'}) G^*_1({\boldsymbol\xi'},{\bf x-h}) \\
&\times& G^*_1({\boldsymbol\xi''},{\bf x+h}) G_0(\boldsymbol\xi'',{\bf x''+h''})I({\bf x_m,h'})I({\bf x_m,h''}).
\end{array}
\label{eq:I_perm3}
\end{eqnarray} 
\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 compute 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 inaccuracies in the migration velocity. The image-point dispersal causes events with different reflection angles from the same reflection point in the subsurface to 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 velocity slower (Figure \ref{fig:pdisp}a) and faster (Figure \ref{fig:pdisp}b) than the true velocity. For simplicity, let us consider constant velocity in the vicinity of the image point, so source and receiver rays are straight.
\par
 When the migration velocity is too low, the reflector is imaged at a shallower 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 that is  twice the apparent reflection angle $\gamma$. Source and receiver rays cross deeper than the image points, causing events to curve downward in the SODCIG. 
\par
When the 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 to curve upward in the SODCIG. 
\plot{pdisp}{width=0.9\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) Underestimated velocity, and b) overestimated 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. Furthermore, accurate dip information is difficult to obtain, especially when events cross because of velocity inaccuracy.
\par
\cite{biondi:1283} point out that, at least to 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, so that all the energy of a point in the subsurface is contained by the corresponding SODCIG injected into the 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 directions of the wavefields locally at the image point. In 2D, $\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$ in equation \ref{eq:talpha} 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 performed 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}}$, with $x_w$ being 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}, reduce 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
%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. 
In 3D, the cross-line offsets also must be rotated according to the apparent geological dip in the cross-line direction in addition to the in-line rotation. By assuming that source and receiver rays are coplanar such that they cross, the 3D transformation to the reflection-angle domain \cite[]{tisserant}  is given by 
\begin{eqnarray}
k_{h_x}=-k_{z} \sec\alpha_{y'}\tan\gamma,
\label{eq:3Dkhx0}
\end{eqnarray} 
and
\begin{eqnarray}
k_{h_y}=-k_{y'} \tan\gamma \tan \alpha_{x'},
\label{eq:3Dkhy0}
\end{eqnarray}
where $\alpha_{x'}$ and $\alpha_{y'}$ are the apparent geological dips in the in-line and cross-line  directions, respectively. Notice that the 3D transformation is dependent on the apparent geological dip in contrast with the 2D case. In the spatial domain, equations \ref{eq:3Dkhx0} and \ref{eq:3Dkhy0} define slant-stack transformations along the paths
\begin{eqnarray}
z=z_{h_{x}}+\frac{h_x}{\cos\alpha_{y'}}\tan\gamma,
\label{eq:3Dkhx1}
\end{eqnarray}
and
\begin{eqnarray}
z=z_{h_{y}}+h_y \tan\alpha_{x'}\tan\alpha_{y'}\tan\gamma,
\label{eq:3Dkhy1}
\end{eqnarray}
respectively. The term $\cos\alpha_{y'}$ in equation \ref{eq:3Dkhx1} stretches the inline-subsurface-offset axis, while the combination of terms $\tan\alpha_{x'}$ and  $\tan\alpha_{y'}$ in equation \ref{eq:3Dkhy1} can stretch or shrink the crossline-subsurface-offset axis. 
\par
Under the common-azimuth approximation, the 3D rotation is similar to the 2D rotation (equation \ref{eq:sltstk1}) except for the subsurface-offset stretching factor (equation \ref{eq:3Dkhx1}). Later in this chapter, we will see areal-shot migration results of 3D-PERM wavefields modeled from the initial conditions computed with and without considering the stretching term.
\par
%Instead of performing a true 3D slant-stack, which can require excessive computer memory, we can take advantage of the fact 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$.
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 a velocity of 1000 m/s and migrated with velocity underestimated by $10\%$ (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}.
\plot{dip01}{width=0.8\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 a depth of 2500 m embedded in a medium with a constant velocity of 1000 m/s.}
\plot{dip02}{width=0.8\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 incorrect 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
\plot{dip03}{width=0.8\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.}
\plot{dip04}{width=0.8\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 in the right panels.
\plot{dip05}{width=0.8\textwidth}{Initial conditions for modeling a) source and b) receiver wavefields. The dipping reflector is oriented in opposite directions in the SODCIG. Rotation affects neither the horizontal reflector nor the-zero subsurface offset, as can be seen in 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 since each reflector was injected separately into the modeling.
\plot{dip06}{width=0.8\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. 
\plot{dip07}{width=0.8\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. The SODCIG on the left is selected at the horizontal position where the dipping reflector was laterally shifted to. 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 the image from wavefields computed with the proposed rotation is much flatter (Figure \ref{fig:dip09}b).
\plot{dip08}{width=0.8\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.}  
\plot{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 using a combination of modeling experiments can decrease the size of PERM data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Combination of modeling experiments}
In the previous examples, if we were to fully image the reflectors, there might be twice as many areal shots as in the original shot-profiles. In equation \ref{eq:I_perm3}, there is a summation over $\bf x_m$ since the initial conditions are isolated SODCIGs. To decrease the number of modeling experiments, we can apply the concept of generalized sources and use the linearity of wavefield propagation 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 {I_D}({\bf x-h}) \\
\widehat{D_{\scriptscriptstyle P}}(x,y,z=z_{\rm max},\omega) = 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 {I_U}({\bf x+h}) \\
\widehat{U_{\scriptscriptstyle P}}(x,y,z=z_{\rm max},\omega) = 0 \end{array} \right.,
\label{eq:comb_bU}
\end{eqnarray}
where $\widehat{I_D}({\bf x-h})$ and $\widehat{I_U}({\bf x+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)$ and $\widehat{U_{\scriptscriptstyle P}}({\bf x},\omega)$, the source and receiver wavefields, respectively. The selection of SODCIGs can be thought of as the multiplication of the pre-stack image by a 2D--comb function, which 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
Since the wavefields are initiated on the reflectors using combinations of SODCIGs, the idea of generalized sources applies, characterizing the image-space generalized source domain.
\par
Again, for one particular frequency and a plane reflector so that the initial conditions are the same for modeling source $\widehat D_{\scriptscriptstyle P}$  and receiver  $\widehat U_{\scriptscriptstyle P}$ wavefields, the modeling of a pair of PERM wavefields starting from combined SODCIGs can be described by
\begin{eqnarray}
\widehat{D_{\scriptscriptstyle P}}(\boldsymbol\xi;\Delta {\bf x}) = \sum_{\bf x}\sum_{\bf h}\sum_m G_0(\boldsymbol\xi,{\bf x-h}) \delta({\widehat{\bf x}-m\Delta{\bf x}})I(\widehat{\bf x}{\bf ,h}),
\label{eq:perm_opDdx}
\end{eqnarray}
and
\begin{eqnarray}
\widehat{U_{\scriptscriptstyle P}}(\boldsymbol\xi;\Delta {\bf x}) = \sum_{\bf x}\sum_{\bf h}\sum_m G_0(\boldsymbol\xi,{\bf x+h}) \delta({\widehat{\bf x}-m\Delta{\bf x}})I(\widehat{\bf x}{\bf ,h}),
\label{eq:perm_opUdx}
\end{eqnarray}
where  $\displaystyle{\sum_m\delta({\widehat{\bf x}-m\Delta{\bf x}})}$ is the 2D--sampling function. 
\par
The wavefields are recursively downward propagated in depth according to 
\begin{eqnarray}
\widehat{D_{\scriptscriptstyle P}}({\bf x};\Delta {\bf x}) = \sum_{\boldsymbol\xi}G^*_1({\boldsymbol\xi},{\bf x}) \widehat{D_{\scriptscriptstyle P}}({\boldsymbol\xi}),
\label{eq:perm_opDzdx} 
\end{eqnarray}
and
\begin{eqnarray}
\widehat{U_{\scriptscriptstyle P}}({\bf x};\Delta {\bf x}) = \sum_{\boldsymbol\xi}G^*_1({\boldsymbol\xi},{\bf x}) \widehat{U_{\scriptscriptstyle P}}({\boldsymbol\xi}).
\label{eq:perm_opUzdx}
\end{eqnarray}
\par
The lateral shifts of the wavefields for the multi-offset imaging condition are represented by
\begin{eqnarray}
\begin{array}{lll}
\widehat{D_{\scriptscriptstyle P}}({\bf x-h};\Delta {\bf x})&=&\displaystyle{\sum_{\boldsymbol\xi}G^*_1({\boldsymbol\xi},{\bf x-h})} \\
&\times& \displaystyle{\sum_{\bf x'}\sum_{\bf h'}\sum_m G_0(\boldsymbol\xi,{\bf x'-h'}) \delta({\widehat{\bf x}-m\Delta{\bf x}})I(\widehat{\bf x}{\bf ,h'})},
\end{array}
\label{eq:perm_lshDdx} 
\end{eqnarray}
and
\begin{eqnarray}
\begin{array}{lll}
\widehat{U_{\scriptscriptstyle P}}({\bf x+h};\Delta {\bf x})&=&\displaystyle{\sum_{\boldsymbol\xi}G^*_1({\boldsymbol\xi},{\bf x+h})} \\
&\times& \displaystyle{\sum_{\bf x'}\sum_{\bf h'}\sum_m G_0(\boldsymbol\xi,{\bf x'+h'}) \delta({\widehat{\bf x}-m\Delta{\bf x}})I(\widehat{\bf x}{\bf ,h'})}.
\end{array}
\label{eq:perm_lshUdx}
\end{eqnarray}
\par
Applying the cross-correlation imaging condition to the wavefields of equations \ref{eq:perm_lshDdx} and \ref{eq:perm_lshUdx} gives
\begin{eqnarray}
\begin{array}{lll}
\widehat{I_{\scriptscriptstyle P}}({\bf x},{\bf h};\Delta {\bf x}) & = & \displaystyle{\sum_{\boldsymbol\xi'}\sum_{\bf x'}\sum_{\bf h'}\sum_m}\displaystyle{\sum_{\boldsymbol\xi''}\sum_{\bf x''}\sum_{\bf h''}\sum_n}G_0(\boldsymbol\xi',{\bf x'-h'}) G^*_1({\boldsymbol\xi'},{\bf x-h}) \\
&\times& G^*_1({\boldsymbol\xi''},{\bf x+h}) G_0(\boldsymbol\xi'',{\bf x''+h''}) \delta({\widehat{\bf x}}-m\Delta{\bf x}) \delta({\widehat{\bf x}}-n\Delta{\bf x}) \\
&\times& I(\widehat{\bf x}{\bf ,h'})I(\widehat{\bf x}{\bf ,h''}),
\end{array}
\label{eq:I_perm2dx}
\end{eqnarray} 
\par  
which can be recast as
\begin{eqnarray}
\begin{array}{lll}
\widehat{I_{\scriptscriptstyle P}}({\bf x},{\bf h};\Delta {\bf x}) & = & I_{\scriptscriptstyle P}({\bf x},{\bf h};\Delta {\bf x}) \\
&+& \displaystyle{\sum_{\boldsymbol\xi'}\sum_{\bf x'}\sum_{\bf h'}}\displaystyle{\sum_{\boldsymbol\xi''}\sum_{\bf x''}\sum_{\bf h''}\sum_{n \ne m}}G_0(\boldsymbol\xi',{\bf x'-h'}) G^*_1({\boldsymbol\xi'},{\bf x-h}) \\
&\times& G^*_1({\boldsymbol\xi''},{\bf x+h}) G_0(\boldsymbol\xi'',{\bf x''+h''}) \delta({\widehat{\bf x}}-m\Delta{\bf x}) \delta({\widehat{\bf x}}-n\Delta{\bf x}) \\
&\times& I(\widehat{\bf x}{\bf ,h'})I(\widehat{\bf x}{\bf ,h''}).
\end{array}
\label{eq:I_perm2dx1}
\end{eqnarray} 
The first term in the right-hand side of equation \ref{eq:I_perm2dx1} is the desired image we would obtain by independently modeling and migrating PERM wavefields. The second term represents crosstalk. To obtain a crosstalk-free image, the sampling period must be large enough that wavefields initiated at different SODCIGs do not correlate. As previously shown when discussing how to compute PERM images with kinematics similar to those of the original shot-profile migration, PERM wavefields generated from SODCIGs within an interval equals to twice the subsurface-offset range $\boldsymbol \eta$ still contribute to the image at the central SODCIG. Crosstalk will occur if the sampling period is shorter than that interval. This is easily seen by realizing that the terms $\delta({\widehat{\bf x'}}-m\Delta{\bf x}) I(\widehat{\bf x'}{\bf ,h'})$ and $\delta({\widehat{\bf x''}}-n\Delta{\bf x}) I(\widehat{\bf x''}{\bf ,h''})$ are 2D--periodic-rectangular functions with period $\Delta{\bf x}$ and width $\boldsymbol \eta$. We want the spatial 2D--correlation of these functions to be a 2D--periodic triangular function, without interference between the individual 2D--triangles. This is achieved by setting $\Delta{\bf x} > {\boldsymbol \eta}$. 
\par
The combination of modeling experiments is illustrated in Figure \ref{fig:comb01b}. We model PERM data starting with the rotated images from the previous section and combine SODCIGs into sets using the sampling period of 163 SODCIGs. 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. 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. 
\plot{comb01b}{width=0.85\textwidth}{Areal-shot migration of PERM data synthesized from sets of SODCIGs selected with sampling period of 163 SODCIGs. Notice that no crosstalk is generated when the sampling period is larger than twice the subsurface-offset range.}
\par
We saw that carefully combining the modeling experiments decreases the data size while maintaining the correct kinematics, 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, at least in 2D. 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, the number of independent experiments is only less than one half as many as in the original dataset. Data reduction techniques like plane-wave decomposition, for instance, could lower data size by a factor of ten. 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 to inject more than one reflector in the modeling.
\par
Although PERM theory was developed in 3D, all the examples I have shown so far have been 2D. Next, I discuss a 3D example under the common-azimuth approximation \cite[]{biondi:1822} and show that in this case the SODCIGs in the $y$ direction can be continuously sampled, and the number of modeling experiments will depend only on the 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 the $x$ and $y$ directions. Considering the usual parameters, the number of modeling experiments may be as low as several hundreds. 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}. CAM reduce the dimensionality of the pre-stack wavefields, and therefore the cost of migration, by assuming zero cross-line offset. That does not mean that the cross-line offset wavenumber is zero. Rather, its asymptotic approximation is a 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 in $\bf x$ and $h_x$. 
\par
Because of the zero-cross-line offset assumption, when using  CAM images 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. 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. Contrast this case with the five-dimensional initial conditions for the full azimuth case of Figure \ref{fig:cam01}a. 
\par
The continuous sampling of SODCIGs in the cross-line direction yields one more order of magnitude of data reduction. Therefore, under the common-azimuth approximation, 3D-PERM data size can be two orders of magnitude smaller than 3D-plane wave data.
\plot{cam01}{width=1.1\textwidth}{The initial conditions for synthesizing PERM data from CAM images can be specified as in b) because no pre-stack information exists in the cross-line direction, in contrast with the full azimuth situation in a).}
\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. There are 96 in-lines and cross-lines spaced 25 m apart. The offset interval is 25 m.  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 zero subsurface offset, with $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 zero subsurface offset, with $x = 750$ m.
\plot{cam02}{width=0.8\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 stretching the inline-subsurface-offset axis by $\sec\alpha_{y'}$ and rotating the SODCIGs in the in-line direction, since no cross-line offset is computed in migration. 
\par
%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.
To illustrate the effect of not stretching the in-line-subsurface-offset axis, source and receiver wavefields were modeled from two different initial conditions: one after simple rotation, and the other after stretching and rotation. Both used continuous sampling along the cross-line direction and sampling period of 48 in the in-line direction. This period is sufficient to avoid crosstalk during the areal-shot migration, given that the number of subsurface-offsets of the pre-stack image is 21. Only the synthesized 3D receiver wavefield for the non stretched case 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.
\plot{cam03}{width=0.8\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.}
The 3D migrations of the 48 areal shots with the velocity underestimated by 5$\%$ are shown in Figures \ref{fig:cam04}a-b for the non stretched case and Figures \ref{fig:cam05}a-b for the stretched case. SODCIGs are positioned at $(x = 750$ m, $y = 600$ m$)$ in Figures \ref{fig:cam04}a and \ref{fig:cam05}a, and at $(x = 750$ m, $y = 1000$ m$)$ in Figures \ref{fig:cam04}a and \ref{fig:cam05}a. To facilitate the comparison with the CAM images of Figures \ref{fig:cam02}a-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 overall matches those of the SODCIGs computed with CAM. However, notice how the images computed with PERM wavefields from the subsurface-offset stretched SODCIGs show slightly better amplitudes at farther subsurface offsets.
\plot{cam04}{width=0.8\textwidth}{3D-areal-shot migration of PERM data from non stretched subsurface offset SODCIGs. 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} and \ref{fig:cam05}.}
\plot{cam05}{width=0.8\textwidth}{3D-areal-shot migration of PERM data from stretched subsurface offset SODCIGs. 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} and \ref{fig:cam04}.}
\par
The correct kinematics shown in Figures \ref{fig:cam04} and \ref{fig:cam05} enables the use of 3D PERM wavefields computed from CAM images in the optimization of migration velocity. As the initial conditions for the modeling are continuously sampled in the cross-line direction, data size is drastically reduced. Chapter \ref{chap:chap06} will use these wavefields to optimize the migration velocity for a 3D survey from the North Sea.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
\section{Conclusions}
In this chapter image-space generalized sources were obtained by combining PERM wavefields. We saw that wavefields synthesized by PERM 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 that will be computed during areal-shot migration of PERM data. Recall that SODCIGs in the initial conditions must be separated by at least twice the maximum absolute subsurface-offset value to prevent crosstalk. Implicit to PERM is the requirement that reflectors must be identified in order to avoid reflector crosstalk during migration. 3D Pre-stack interpretation can be cumbersome, but it allows avoiding the use of reflectors with low signal-to-noise ratios in migration-velocity estimation. Moreover, in commercial software for migration velocity estimation reflector picking is a standard and almost entirely automated procedure.
\par
Whereas in 2D PERM data size is comparable to that of the plane-wave decomposition, in 3D it is one order of magnitude smaller when computing cross-line subsurface-offsets. Further data size reduction by another order of magnitude is achieved if the initial conditions are computed with common-azimuth migration.
%
\par
\cleardoublepage
%\bibliographystyle{seg}
%\bibliography{SEP,SEG,MISC,GEOTLE,allrefs}

\end{document}
