\def\w{\omega}
\def\bX{\boldsymbol{x}}
\def\bxi{\boldsymbol{\xi}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{pre-stack-exploding-reflector model} 
The fundamental idea of PERM is to model data that describes the correct kinematics of an isolated SODCIG computed by a wavefield-extrapolation method. In conventional Born modeling, as we do not know beforehand which shot records contribute to forming the image at a point in the subsurface, we would have to perform several modeling experiments to synthesize these shot records.  Ideally, instead of several modeling experiments, we would like to synthesize a small amount of data with the condition that migration produces the same kinematics as the initial SODCIG. This can be achieved by using the prestack exploding-reflector model (PERM). A prestack image is used as the initial conditions for the modeling of upward and downward PERM wavefields. The modeling is performed in a way similar to the exploding-reflector model (ERM). Since migration of PERM data can produce a prestack image, PERM can be considered a generalization of ERM.
\par
The modeling of PERM  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} \pm i\sqrt{\omega ^2 s_0^2({\bf x})-|{\bf k}| ^2}\right) P_{\scriptscriptstyle u,d}({\bf x},\omega;{\bf x_m}) = I_{\scriptscriptstyle u,d}({\bf x_m},{\bf h}) \\
P_{\scriptscriptstyle u,d}(x,y,z=z_{\rm max},\omega;{\bf x_m}) = 0 \end{array} \right.,
\label{eq:perm_bUD}
\end{eqnarray}
where $P_{\scriptscriptstyle u,d}$ represent upward and downward PERM wavefields, $I_{\scriptscriptstyle u,d}({\bf x_m},{\bf h})$ are subsets of the prestack image centered around the SODCIG under study at ${\bf x_m}$ for a single reflector, suitable for the initial conditions for the upward and downward PERM wavefields, respectively. The sign in front of the square root is positive for the upward PERM wavefield and negative for the downward PERM wavefield. The subsurface-offset $\bf h$ is parameterized as ${\bf h}=(h_x,h_y)$, where $h_x$ and $h_y$ are the inline- and the crossline-subsurface offsets, respectively. If the prestack image has energy focused at zero subsurface offset, the initial conditions can be parameterized only by its spatial coordinates, and PERM is equivalent to ERM.  
\par
The initial conditions are obtained by rotating the original unfocused SODCIGs according to the apparent geological dip of the reflector.  This corrects for the image-point dispersal due to velocity inaccuracy, causing events with different reflection angles from the same reflection point in the subsurface to be imaged at different locations. The geometry of the rotation is shown in Figures \ref{fig:pdisp}a and \ref{fig:pdisp}b for velocities which are slower and faster than the correct velocity, respectively.  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 in-line-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 cross-line-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. Appendix 1 shows the complete derivation of equations *** and ***.



\begin{eqnarray}
\tan (\gamma - \alpha) = \frac{dz}{d\widetilde {h_{x_u}}}, 
\label{eq:theta1}
\end{eqnarray}
\begin{eqnarray}
\tan (\gamma + \alpha) = -\frac{dz}{d\widetilde {h_{x_d}}}.
\label{eq:theta2}
\end{eqnarray}
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:


It is represented the situation of a single subsurface offset Notice that when the migration velocity is slower, the image point  . 
\par
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 in-line-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 cross-line-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
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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 in-line-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}

