\documentclass[manuscript]{geophysics}
\newcommand{\rs}[1]{\mathstrut\mbox{\scriptsize\rm #1}}
\newcommand{\rr}[1]{\mbox{\rm #1}} 
\usepackage{algorithm}	
\usepackage{algorithmic}	
\usepackage{amsmath}	
\usepackage{multicol}
\usepackage{placeins}
\newcommand{\remark}[1]{{ \bf [ \footnotesize #1 ]}}  
\begin{document} 
%: 
\title{Migration-velocity analysis using image-space generalized wavefields}
\renewcommand{\thefootnote}{\fnsymbol{footnote}}  

\ms{GEOPHYSICS} % manuscript number 

\address{ 
\footnotemark[1]Stanford Exploration Project, \\ 
397 Panama Mall, Stanford, CA, 94395 \\ 
} 
\author{Claudio Guerra and Biondo Biondi}

%\def\figdir{./Fig}
\righthead{MVAing - Norne}
\lefthead{Ayeni}
\footer{Norne 4D}
  
\begin{abstract}
	Seismic matching (or cross-equalization) attenuates unwanted or non-production related artifacts from time-lapse seismic data sets.
	Two robust post-imaging matching methods are considered. 
	First, an efficient cyclic 1D correlation method is used to estimate vertical and lateral displacement components between baseline and monitor images. 
	Next, a cyclic 1D match-filtering method that estimates the optimal filtering parameters attenuates residual artifacts.
	These two methods provide a robust time-lapse cross-equalization scheme.
	Using the Norne time-lapse seismic data sets, we show that in order to obtain accurate estimates of production-related changes, all displacement components must be taken into account.
\end{abstract}

\section {Introduction} 
	Time-lapse (4D) seismic provides important information about production-related changes in reservoir properties.
	Over the past decade, time-lapse seismic technology has developed into an integral part of most hydrocarbon reservoir management projects.
	Several authors have discussed important technical considerations and the business impact of time-lapse seismic \cite[]{johnston:fulmar,rickett:2001,whitcombe:2251,zou:764,ebaid:holstein}. 
\par
	Imperfect survey geometry repetition, non-repeatable ambient noise and other discrepancies generate artifacts that contaminate production-related differences between seismic data sets.
	Before reservoir property changes can be extracted from time-lapse seismic data, these artifacts must be attenuated. 
	The process of attenuating artifacts from time-lapse data sets is generally called seismic cross-equalization \cite[]{rickett:2001}.
	In this paper, two robust post-imaging cross-equalization steps are considered. 
\par
	Because reservoir depletion changes the stress-states in and around a producing reservoir, seismic travel times and path lengths are different between surveys.
	However, in practice, the baseline and monitor data sets are migrated with a single (baseline) velocity, leading to \emph{apparent} vertical and lateral displacements between seismic images.
	Because it is expensive to estimate all components of the displacement vector field, it is common practice to only consider vertical displacements.
	Instead of this \emph{vertical-only} approximation, displacement components are efficiently computed with a cyclic 1D correlation method \cite[]{hale:V99}.
	It is shown that in order to obtain accurate estimates of reservoir property changes, all displacement components must be taken into account.
\par
	To ensure that only production-related changes within the reservoir are interpreted, differences in seismic signal in non-reservoir regions should be minimal.
	After aligning all monitor images to the baseline (by removing the estimated apparent displacements), an optimal cyclic 1D match-filtering method is used to attenuate differences in the non-reservoir regions.
	This cyclic filtering method is an extension of the optimal matching method described by \cite{gboyega:seg3}.
	Together with the cyclic correlation method described above, this matching method forms a robust and effective seismic cross-equalization scheme. 
\par
	In this paper, I first briefly discuss the cyclic shifting and match-filtering methods.
	Next, I apply these to four time-lapse seismic data sets from the Norne field. 
	The Norne study shows that neglecting lateral displacements can cause spurious time-lapse signals.
	Finally, I show that these cross-equalization methods significantly improve the quality of the time-lapse images and their derivatives. 
%
\section {Apparent displacement vectors by cyclic 1D search}
	Components of the apparent displacement vectors between time-lapse seismic data sets can be efficiently computed using a cyclic 1D search approach \cite[]{hale:V99}.
	This method involves a sequential cyclic 1D search for cross-correlation peaks in different directions.
	This removes the need for a computationally expensive 3D search for cross-correlation peaks or the inaccurate assumption that displacements are only in the vertical direction.
	To improve the resolution of displacement estimates in all directions, the input data are first band-passed to remove very low frequencies. 
	This preprocessing step serves a similar purpose as the prediction-error filtering step of \cite{hale:V99}.
\par
	This implementation utilizes correlation gates and lags of varying sizes and an acceptance criterion based on the relative values of the cross-correlation peaks.
	To first capture the low-frequency, high-amplitude displacements, this iterative procedure starts with large correlation gates and lags.
	The correlation gates and lags are then decreased as a function of iteration.
	Displacement estimates are filtered based on the displacement size and cross-correlation peak.
	For each iteration, the results (shifted image(s) and accumulated displacements) from one direction (e.g. vertical) are the inputs to the next direction (e.g. in-line).
	For each direction, shifts are accumulated after interpolating the image and displacement components from the previous step.
	Algoirthm~\ref{alg1} summarizes this procedure. 
\par 
	The estimated displacements are used to align the baseline and monitor data sets prior to time-lapse image computation.
	They are also useful for estimating velocity and geomechanical changes.
%
\subsection {Velocity, path-length and travel time changes}
	Depletion-related strain and velocity changes will perturb the travel path and length of a normally-incident ray.
	The distance $dz$ travelled by the ray in time $dt$ at velocity $v$ is given by
		\begin{equation} 
			dz = v dt.
		\end{equation}
	For small perturbations $\delta z\ll dz$, $\delta v\ll dv$ and $\delta t\ll dt$, we can write 
		\begin{equation}
			 \frac {d(\delta z)}{dz}  = \frac{\delta v}{v} + \frac {d(\delta t)}{dt},
		\end{equation}
	where, $d(\delta z)/dz=\epsilon_{zz}$ is the vertical ``depth'' strain and $d(\delta t)/dt=\epsilon_{tt}$ is the apparent ``time'' strain.
	The $R$-factor \cite[]{hatchell:05} relates the fractional change in velocity ${\delta v}/{v}$ to the depth strain $\epsilon_{zz}$ and time strain $\epsilon_{tt}$: 
		\begin{equation}
			\frac{\delta v}{v} = -R \epsilon_{zz} = -\frac{R}{1+R}  \epsilon_{tt}.
			\label{eqn:dzt}	 
		\end{equation}
	The $R$-factor can be measured from core samples in the laboratory \cite[]{hatchell:05,bathija:dilation}, derived from the estimated time-shifts \cite[]{hawkins:stress}, or derived from other production-related data \cite[]{hodgson:pressure}.
	For most reservoir rocks, $R$ values of $4$ to $6$ have been shown to be reasonable \cite[]{hatchell:compaction,hodgson:pressure}.
	In this paper, a constant $R$ value of $5$ is assumed.
%
%
\begin{figure}
\begin{algorithm}[H]                      % enter the algorithm environment
\caption{Apparent displacements by cyclic 1D search}          % give the algorithm a caption
\label{alg1}                           	% and a label for \ref{} commands later in the document
\begin{algorithmic}                    	% enter the algorithmic environment 
\STATE $b \Leftarrow preprocess(Base)$
\STATE $m \Leftarrow preprocess(Monitor)$  
	\WHILE{$s^{ia=1:3}_{max} \geq s^{ia=1:3}_{stop}; j=1:niter$}
		\FOR{$axis=1$:$3$} 
		\IF{$axis=1$}
			\STATE ${ s^{1}} = s[xcorr_{peak}(m_{j} | b)]$\   {!compute displacements btw $m_{j}$ \& $b$ along axis 1}
			\STATE $ z^{1} \Leftarrow { s^{1}} <z^{1}>+ {s^{1}}$ \ \ \ \  {!interpolate and accum. displacements along axis 1} 
			\STATE $z^{2} \Leftarrow { s^{1}} <z^{2}>$  \ \ \ \ \ \ \ \ \ \   {!interpolate displacements along axis 2} 
			\STATE $z^{3} \Leftarrow { s^{1}} <z^{3}>$ \ \ \ \ \ \ \ \ \ \    {!interpolate displacements along axis 3} 
			\STATE $m^{1}_{j} \Leftarrow { s^{1}} <m_{j}>$ \ \ \ \ \ \ \ \   {!interpolate monitor along axis 1}
		\ELSIF{$axis=2$}
			\STATE ${ s^{2}} = s[xcorr_{peak}(m^{1}_{j} | b)]$
			\STATE $z^{1} \Leftarrow { s^{2}} <z^{1}>$ 
			\STATE $ z^{2} \Leftarrow { s^{2}} <z^{2}>+  {s^{2}}$
			\STATE $z^{3} \Leftarrow { s^{2}} <z^{3}>$
			\STATE $m^{2}_{j} \Leftarrow { s^{2}} <m^{1}_{j}>$
		\ELSIF{$axis=3$}
			\STATE ${ s^{3}} = s[xcorr_{peak}(m^{2}_{j} | b)]$ 
			\STATE $z^{1} \Leftarrow { s^{3}} <z^{1}>$ 
			\STATE $ z^{2} \Leftarrow { s^{3}} <z^{2}>$
			\STATE $z^{3} \Leftarrow { s^{3}} <z^{3}>+  {s^{3}}$
			\STATE $m^{3}_{j} \Leftarrow { s^{3}} <m^{2}_{j}>$
			\STATE $m_{j+1} \Leftarrow m^{3}_{j}$
		\ENDIF
		\ENDFOR
	\ENDWHILE
	\STATE $Monitor \Leftarrow z^{1:3}<Monitor>$ 
\end{algorithmic} 
$s$: displacements \\ 
$z$: cumulative displacements \\ 
$b$: preprocessed base \\ 
$m$: preprocessed monitor \\
$xcorr_{peak}$: correlation peak \\
${ a} < b>$: interpolation of $b$ with displacements $a$ \\
\end{algorithm} 
\end{figure}
%
\section {Cyclic 1D match-filtering by Evolutionary Programming}
	The match-filtering problem between two data sets, one baseline and one monitor, can be posed in the frequency domain as follows:
\begin{equation}
	\tensor{F}(\omega)=\frac{\tensor{\overline M}(\omega)\tensor{B}(\omega)}{\tensor{\overline M}(\omega)\tensor{M}(\omega)+\epsilon^2}.
\label{eqn:filt6}
\end{equation}
	where $\tensor{F}(\omega)$ is the match filter, $\tensor{B}(\omega)$ is the baseline data, and $\tensor{M}(\omega)$ is the monitor data at frequency $\omega$.
	The damping factor ($\epsilon$) in the denominator prevents zero-division.
	The overline denotes complex conjugation.
	Match filters can be estimated globally (for all trace locations), locally (on a trace-by-trace basis) or both.
	In this paper, only local 1D match filters are considered. 
\par	
	At each trace location, local match filters are non-linear functions of the estimation parameters.
	Usually, these parameters are selected by a manual trial-and-error approach.
	This can be a tedious or impossible challenge due to the strong non-linearity of the problem and the large size of seismic data sets. 
	Furthermore, in most cases a single set of parameters is used for all traces locations.
	While a particular parameter set may result in satisfactory filter estimates in some parts of the data volume, it may result in very poor filters in other parts.
\par 
	Here, an evolutionary programming (EP) method is used estimate the optimal filter parameters at each trace location.
	EP belongs to a class of global optimization methods called evolutionary algorithms.
	These algorithms solve optimization problems using Darwinian evolutionary principles of natural selection.
	Although any global optimization method may be used, the EP method was chosen because of its easy implementation, easy parameterization and flexibility.
	The optimized estimation parameters include estimation window size and taper, filter-length and the damping factor. 
	\cite{gboyega:seg3} provide a more detailed description and an application of evolutionary programming to the match filtering problem.
	An extended version of this filter estimation approach that is used in this paper is presented in Algorithm~\ref{alg2}.
%
\begin{figure}
\begin{algorithm}[H]                      % enter the algorithm environment
\caption{Optimal cyclic 1D matching by evolutionary programming}          % give the algorithm a caption 
\label{alg2}                           	% and a label for \ref{} commands later in the document
\begin{algorithmic}                    	% enter the algorithmic environment
\FOR{$irepeat=1$:$nrepeat$}
	\FOR{$itrace=1$:$ntrace$}
		\STATE  {(Re)Initialize}: $estimation\ window, filter\ length, damping, etc.$
		\STATE {Evaluate}:  $fitness\ (nrms)$ 
		\WHILE{$nrms_{iter} \geq nrms_{stop}; iter=1:niter$}
	  	 	\STATE  {Select}: ${fittest}\ individuals$ 
	 	 	\STATE  {Mutate}: $new\ offsprings$ 
	  	 	\STATE  {Evaluate}: $fitness\ (nrms)$
		\ENDWHILE
	\ENDFOR
 	\STATE Apply constraints  
\ENDFOR 
\end{algorithmic} 
\end{algorithm}
\end{figure}
%
\section {Case study}
	The Norne field is located in the Norwegian North sea.
	The main field is a 9 km x 3 km horst block comprised of high-porosity, high-permeability, high net-to-gross lower and middle Jurassic sandstones.
	The field was discovered in 1991 and production started in 1997.
	For this study, four high-resolution Q-marine time-lapse seismic data sets acquired in 2001, 2003, 2004 and 2006 are considered.
	The 2001 data is taken as the baseline.
	\cite{osdal:1134} and \cite{aarre:3174} provide more detailed description of the geology, production history, time-lapse seismic acquisition and processing, and some interpretation of the time-lapse seismic response at the Norne. 
\par
	The apparent vertical displacements (time-shifts) for three cycles between the 2006 and the 2001 data are shown in Figure~\ref{fig:z-shift-1,z-shift-2,z-shift-3}.
	In this figure, and in similar ones throughout this paper, the orthogonal lines show slices through a 3D volume at the positions indicated by each line.
	Figure~\ref{fig:z-shift-3-TopNot} shows the 2006 time-shifts together with the mapped top of the reservoir.
	The corresponding lateral displacements in the inline and cross-line directions (after three cycles) are shown in Figure~\ref{fig:x-shift-3,y-shift-3}.
	Figures~\ref{fig:shift-compare1}~to~\ref{fig:shift-compare3} show the absolute vertical and lateral displacements overlain by arrows indicating the displacement direction.
	Estimated velocity changes at different production stages are shown in Figure~\ref{fig:dvv-shift-2003,dvv-shift-2004,dvv-shift-2006}.
	It is assumed that the time-shift estimates obtained by considering displacements in all directions are accurate.
	Differences between time-shifts, velocity changes and vertical strains obtained by considering all displacement components and those obtained by considering only vertical displacements are shown in Figure~\ref{fig:z-shift-3-error,dvv-shift-2006-error,ez-shift-2006-error}.
%
	The optimal match-filter parameters obtained by the evolutionary programming method are shown in Figure~\ref{fig:optimal-1}.
	Figures~\ref{fig:amp-shift-2006,amp-shift-2006-shifted,amp-shift-2006-matched}~and~\ref{fig:amp-shift-2006-b,amp-shift-2006-b-shifted,amp-shift-2006-b-matched} show time-lapse images obtained after different processing steps.
	Amplitude errors resulting from considering only vertical displacements are shown in Figures~\ref{fig:amp-shift-2006-b-error}.
%
\section{Discussion}
	The sequential cyclic 1D search method provides a robust alternative to a full 3D search for cross-correlation peaks. 
	By iteratively estimating displacements in a cyclic manner using variable correlation gates and lags, it is possible to obtain high-resolution estimates of displacement components (Figures~\ref{fig:z-shift-1,z-shift-2,z-shift-3} to \ref{fig:x-shift-3,y-shift-3}).
	These displacement components contain different information about changes within and around the reservoir and may be interpreted qualitatively (Figures~\ref{fig:shift-compare1}~to~\ref{fig:shift-compare3} ).
	Analysis of the estimated displacements show that in general, velocity within the reservoir decreses with time (Figure~\ref{fig:dvv-shift-2003,dvv-shift-2004,dvv-shift-2006}).
	A velocity decrease of up to 5 percent is observed between 2001 and 2006 (Figure~\ref{fig:dvv-shift-2006}).
\par 
	There are inherent errors caused by considering only vertical displacements in time-lapse seismic data sets (Figure~\ref{fig:z-shift-3-error}).
	In the current example, this approximation leads to time-shift overestimation below the reservoir, and time-shift underestimation above the reservoir.
	Errors in time-shift estimates, which are more pronounced around dipping reflectors, will map into errors in velocity change and vertical-strain estimates (Figures~\ref{fig:dvv-shift-2006-error} and \ref{fig:ez-shift-2006-error}).
\par
	Filter parameters derived from the proposed method (Figure~\ref{fig:optimal-1}) enable computation of the optimal filters at each trace location.
	Before computing time-lapse images, the monitor data sets are shifted into alignment with the baseline by interpolating them with the estimated displacements.
	After data alignment (Figures~\ref{fig:amp-shift-2006-shifted} and \ref{fig:amp-shift-2006-b-shifted}), optimal match-filtering removes residual artifacts, thereby improving the time-lapse data-quality (Figures~\ref{fig:amp-shift-2006-matched} and \ref{fig:amp-shift-2006-b-matched}).
	Finally, especially at steeply dipping flanks, considering only vertical displacements can lead to erroneous time-lapse amplitude changes (Figure~\ref{fig:amp-shift-2006-b-error}).
%
\section{Conclusions}
	Time-lapse seismic data sets can be efficiently cross-equalized by cyclic 1D matching methods.
	The cyclic 1D shifting and optimized match-filtering methods form an efficient and robust cross-equalization scheme.
	From the case study, it is apparent that in order to obtain accurate estimates of reservoir property changes, all displacement components must be taken into account.
	Where estimated parameters are important for decision making, considering only vertical displacements may result in serious undesirable consequences.
	Future work will consider ways to fully integrate all displacement components in time-lapse seismic imaging.
%	
\section{Acknowledgements} 
	Thanks to Statoil for the release of the Norne data. 
	Thanks to the NTNU/IO CENTRE for their cooperation and coordination of the Norne Cases. 
	Special thanks to Trine Alsos (Staoil), Mohsen Dadashpour (NTNU) and David Echeverria (Stanford Smart Fields) for their continued assistance throughout this study.
%%
%\multiplot{3}{z-shift-1,z-shift-2,z-shift-3}{height=0.255\textheight}{Apparent vertical (time-) shifts estimates between the baseline and the 2006 data sets after one (top), tow (middle) and three (bottom) cycles. 
%Note the improvement in resolution with the number of cycles. [{\bf CR}].} 
%%
%\plot{z-shift-3-TopNot}{height=0.255\textheight}{Time-shifts between the baseline and the 2006 data and the top of the reservoir. 
%Note the similarity between the displacements correspond to known producing segments in the filed. [{\bf CR}].}
%%
%\multiplot{2}{x-shift-3,y-shift-3}{height=0.255\textheight}{Apparent lateral displacements between the baseline and the 2006 data. [{\bf CR}].}
%%
%\plot{shift-compare1}{height=0.7\textheight}{Absolute vertical and lateral displacements at the top of the reservoir between the baseline and the 2003 monitor.
%The arrows indicate the displacement direction. [{\bf CR}].}
%%
%\plot{shift-compare2}{height=0.7\textheight}{Absolute vertical and lateral displacements at the top of the reservoir between the baseline and the 2004 monitor.
%The arrows indicate the displacement direction. [{\bf CR}].}
%%
%\plot{shift-compare3}{height=0.7\textheight}{Absolute vertical and lateral displacements at the top of the reservoir between the baseline and the 2006 monitor.
%The arrows indicate the displacement direction. [{\bf CR}].}
%%
%\multiplot{3}{dvv-shift-2003,dvv-shift-2004,dvv-shift-2006}{height=0.255\textheight}{Fractional velocity change between the baseline and the 2003 (a), 2004 (b) and 2006 (c) monitor data. 
%Note that the velocity within the reservoir decreases with time. [{\bf CR}].}
%%
%\multiplot{3}{z-shift-3-error,dvv-shift-2006-error,ez-shift-2006-error}{height=0.255\textheight}{Time-shift (a), velocity change (b) and vertical strain (c) errors between the baseline and the 2006 data as a result of assuming on vertical displacements. 
%In each case, the result obtained by considering only vertical displacements were subtracted from those obtained by considering all directions with the cyclic search method. 
%The red line indicates the reservoir top. [{\bf CR}].}
%%
%\plot{optimal-1}{height=0.93\textheight}{Optimal filter estimation parameters obtained by the evolutionary programming method. [{\bf CR}].}
%%
%\multiplot{3}{amp-shift-2006,amp-shift-2006-shifted,amp-shift-2006-matched}{height=0.255\textheight}{Time-lapse images after different processing steps. 
%The top panel is a time-slice through the reservoir.
%Note the improvements in the quality of time-lapse amplitudes within the reservoir from top to bottom. [{\bf CR}].}
%%
%\multiplot{3}{amp-shift-2006-b,amp-shift-2006-b-shifted,amp-shift-2006-b-matched}{height=0.255\textheight}{Time-lapse images after different processing steps. 
%The top panel is a time-slice below the reservoir region where no time-lapse amplitude changes are expected. 
%Note the improvement in the data quality from top to bottom. [{\bf CR}].}
%%
%\plot{amp-shift-2006-b-error}{height=0.255\textheight}{Errors in time-lapse amplitude computed as a difference between the time-lapse image obtained by considering only vertical displacements and from considering all displacement components. [{\bf CR}].}
%%
%%\newpage
\bibliographystyle{seg}
\bibliography{migration}

%\newpage
%\fbox{\includegraphics[width=\textwidth]{/net/koko/wrk1/sep140/gayeni1/Fig/filler5.pdf}}
\end{document}
