grmdoc.tex

su 的源代码库

{\documentstyle[authordate1-4,psfig]{project96}
\newtheorem{algorithm}{Algorithm}
\begin{document}
\title{GRM - Generalized Reciprocal Method \\
 with automated optimum parameters and depth estimation}\shorttitle{GRM}
\author[Steven D. Sheaffer]{Steven D. Sheaffer\\
Center for Wave Phenomena\\
Colorado School of Mines}

\shortauthor{Sheaffer}

\maketitle

\begin{abstract}

\end{abstract}

\begin{keywords}
\end{keywords}

\section*{Introduction}

Refraction methods for determination of the depth to shallow interfaces
are utilized in both near-surface investigations and deep reflection surveys.
For reflection data processing,
the depth to the base of the low velocity surface layer (the "weathering"
layer) is useful for the determination of static corrections.  
One class of methods involves the use of first arrival times of 
refracted "head" waves generated by an increase in velocity across the
interface.

For the simple case of a plane dipping refractor with constant velocity
above and below, with or without surface topography, the well known plus-minus 
method (Hagedoorn, 1959) will yield the depth normal to the interface
at all surface points.  The method requires a forward and reverse shot,
producing a travel time from each at every geophone location between them.
Figure 1. shows the raypath for fixed shots A and B and a particular midpoint
$i$.  If the function
	\begin{center}		
	\( T^{-}(i) =  T_{ai}-T_{bi}+T_{ab} \)
	\end{center}
is evaluated at each midpoint $i$ for
fixed A and B, the result will be a line with slope equal to the reciprocal 
of the apparent refractor velocity $ 1/v'_{r}$.  A second function, 
	\begin{center}
	\( T^{+}(i) = T_{ai}+T_{bi}-T_{ab} \)
	\end{center}
can be seen to be equivalent to a solution over the inner triangle, which is
evaluated to be,
	\begin{center}
	\( T^{+}(i) = T_{a'i}+T_{b'i}-T_{a'b'} =
	\ \displaystyle{\frac{2d\cos\theta_{cr}}{v_{0}}} \)
	\end{center}
where $v_{0}$ is the known near surface velocity, 
and $\theta_{cr}$ is the critical angle, which can be estimated from
$v_{0}$ and $v'_{r}$ 
(since the dip is not known but is assumed small, the apparent velocity is 
taken to be the true velocity.) 
The normal depth \textit{d} is found by multiplying $T^{+}(i)$ with a 
velocity function that accoutns for the cosine, 
	\begin{center}
	\(d(i) = \frac{1}{2} T^{+}(i)\displaystyle{\frac{v_{0}}{\cos\theta_{cr}}} \)
	\end{center}
In reality, however, we do not want to be limited to plane refractors, but are primarily 
interested in imaging irregular refractor surfaces.  It is easlily seen
that the plus-minus method will break down for a non-planar surface.  Since the 
refractor surface 
at $a'$ and $b'$ will most likely have different dips, triangle geometry will
not form, and the method will produce erroneous results.  
The problem of dip conflict can be overcome by adopting the ray geometry shown
in Figure 2., where both rays leave the irregular surface at the same point,
and are subject to the same dip.  The geometry of the inner triangle is
now oblique, but still suggests solution.   It is this "optimal" geometry that
is the basis of the Generalized Reciprocal Method (GRM) developed by 
Palmer (1982).

To use the GRM geometry, some modification of the previous theory is required,
not only to account for the new geometry but also to address the fact that 
we do not know 
what separation XY on the surface corresponds to the rays leaving the
refractor at a common point (the optimium XY).  The geophone location \textit{i} 
is defined as halfway between X and Y, as measured on the surface.
The Velocity Analysis
Function is analogous to a one-way version of $T^{-}$,
	\begin{center}
	\(T_{v}(i) = \frac{1}{2}(T_{ay}-T_{bx}+T_{ab}) \)
	\end{center}
The reciprocal of its slope will give $V'_{r}$, as before.  Observing
Figure 2., it should 
be clear that if the rays leave the refractor at different points, $T_{v}(i)$
will be perturbed by the additional path between them.  Since the 
refractor is irregular, this additional path will vary in length along
the refractor, causing the function $T_{v}(i)$ to be rough.  Therefore, if it is evaluated
for various values of XY, the smoothest will give the optimum XY.  This
is the value of the offset that, on average, will give rays which come
closest to a common take off point for each midpoint \textit{i}.

Once XY is determined, the Time-Depth function is a modified $T^{+}$:
	\begin{center}
	\(T_{g}(i) = \frac{1}{2}(T_{ay}+T_{bx}-T_{ab}-XY/V'_{r})  \)
	\end{center}
For zero dip, this is analogous to a 1-way version of $T^{+}$, but in the 
presence of dip, its geometry is oblique, and does not give the exact 
normal depth.  Again, GRM converts $T_{g}(i)$ into its estimate of depth 
by multiplication with the cosine scaled velocity function as in the 
Plus-minus calculation, so that
	\begin{center}
	\(d_{grm}(i) = T_{g}(i) \displaystyle{\frac{v_{0}}{\cos\theta_{cr}}} \)
	\end{center}
To understand the relationship between $d_{grm}$ and the actual depth normal
to the refractor, $d$, a reference line $d'$ can be drawn on Figure 2.
connecting the midpoint $i$, halfway between X and Y and the take-off 
point on the refractor.  Calculation of $d_{grm}$ over the oblique triangle
gives the result that 
	\begin{center}
	\( d_{grm} = d' \left[ \displaystyle{\frac{\cos^{2}\alpha - \sin^{2}\theta_{cr}\cos\alpha} 
                       {\cos^{2}\alpha - \sin^{2}\theta_{cr}} }\right]   \) 
	\end{center}
where $\alpha$ is the local dip. The pole occurs at 
	\begin{center}
	\(\theta_{cr} + \alpha = \displaystyle{\frac{\pi}{2}} \)
	\end{center}
which is the configuration where the refracted wave is parallel to the surface,
and the triangle becomes infinitely long.

It is clear from this that GRM will only yield accurate results for small dips,
since:
	\begin{center}
	\( d_{grm}\approx d' \)
	\end{center}
and,
	\begin{center}
	\( d' \approx d \)
	\end{center}
only when $\alpha$ is small.

\section*{SU Implementation of GRM}

The current version of GRM will define the depths to
a single refractive interface, given a set of previously determined first
arrival times.  The method can be extended to multiple refractors,
given a set of first arrivals for each, where each interface is 
modeled as a single interface overlain by a layer of constant velocity
approximately equal to the average velocity of the actual structure
above (see Palmer, 1982).  So, the code can be applied once for each
refractor from the surface down, and previous results used to determine
the average velocity for the next calculation. 

The GRM method is normally implemented interactively. Optimum XY values
are chosen by inspection of plots of $T_{v}(i)$ using various XY values,
and visually determining the smoothest function.
GRM depths are estimates of normal depths, and might be thought of as 
being, in essense, "unmigrated", since they are normal to an unknown
point on an interface with unknown local dip.
These depths are often not converted to true depths, under the assumption
of small dips, where they would be approximately equal to vertical depths.
Alternatively, they can also be converted to vertical depths interactively
by plotting arcs with radii
equal to the GRM depth at each geophone ($i$), and then interpolating an
envelope along the tangents of the arcs. 

GRM attempts to remove
the interactive steps, and automatically determine these values.
Optimum XY is determined by investigating $T_{v}(i)$ for a range of test XY values,
and summing the finite-difference approximation of the Laplacian 
	\begin{center}
	\(\nabla^{2} T(i)\approx \displaystyle{\frac{T_{v}(n+1) - 2T_{v}(n) + T_{v}(n-1)}
	       {dx^{2}}} \)
	\end{center}
over the length of the data.  The function with the smallest total Laplacian
is assumed to be the smoothest, giving XY optimum.  The maximum length of $T_{v}(i)$
for some XY, though, is dependent upon XY.
The function must reach back to $x(i)-XY/2$ and forward to $x(i)+XY/2$ for times 
required by the calculation, cutting off the ends of the available data.
To keep the comparison consistent, all of the $T_{v}(i)$
used in the comparison are calculated to the length allowed by the largest
test XY.  The user can choose the maximum XY used in searching for XY optimum
by setting the optional parameter $\texttt{xymax}$.  Note, though, that
this means setting large values of $\texttt{xymax}$ will cause the summation
to occur over smaller regions of $T_{v}(i)$ to determine optimality.  A user-specified
XY can be be defined by setting the optional parameter
$\texttt{xy}$, which 
will then be used as XY in the calculation.  

Once the optimum XY has been defined, the average slope of the resulting optimum 
$T_{v}(i)$ is found by summing the forward difference
	\begin{center}
	\(\nabla T(i) \approx \displaystyle{\frac{T_{v}(n+1) - T_{v}(n)}{dx}}    \)
	\end{center}
and dividing by the length.  The reciprocal of this average gives $V'_{r}$.
The reciprocal of the average slope is used instead of the average of the 
reciprocals of the slope at each point to reduce the effect of small local
slopes. 

Now, the depths can be found using the time-depth function and the velocity
function, both of which are fully defined with the addition of the optimum XY
and $V'_{r}$, detemined previously.  Again, the range of midpoints available for
the calculation depend on XY, since they function reaches back to $x(i)-XY/2$
and forward to $x(i)+XY/2$ for times.

Once the GRM depths are determined, the program attempts to convert them
to vertical depths.  Since the refraction points and dips are not known,
it will attempt to 
find a surface, tangent to arcs with radii $d_{grm}(i)$ around each point $i$.
The procedure is to take each pair of adjacent midpoints, $i$ and $i+1$, and 
assume a half-circle around each with radius equal to the GRM depth at 
each location.  Then, determine the line that is tangent to both circles.

Consider the two points $x_{i},y_{i}$ on circle the circle surrounding $i$, 
and $x_{i+1},y_{i+1}$
on the circle surrounding $i+1$, where the slope of the tangent at both points is equal.
Then if the slope of 
the line $x_{i},y_{i}$ to $x_{i+1},y_{i+1}$ between them has a slope that
is equal
to their mutual tangent slope, the line is tangent to both circles at the 
points specified.  

This line is found by moving along the half-circle centered at $i$, calculating the slope
of the tangent at each point $x_{i},y_{i}$ on the half-circle.  For each of these
points, the point on the adjacent circle with the same tangent slope is the
one that is the same angular distance around the arc.  The slope of the line
between these two points is compared to the tangent slope at $x_{i},y_{i}$.
 Since
the calculation occurs at discretized points along circle $i$, the slopes
will not be exactly equal, so a residual of the difference in the two slopes
is recorded for the entire half-circle , and the smallest residual is chosen
as the slope of tangency to both.  The distance below geophone $i$ of this 
line is found and reported as the depth to the refractor.

This discritization is controlled by the optional parameter \texttt{depthres},
which
allows the user to set the size of the increment in $x$ that is used when finding
the tangent.  The value of the residual in the slope is converted to a 
error in depth and reported in the output file.  If these errors
are deceided to be too large by the user, a smaller value of \texttt{depthres}
be used to try to "zero-in" on a better match.

This method has the possiblity of overestimating the depth where the difference
in radii at $i$ is much larger than that at $i+1$, but since the radii are
GRM determined depths, these
values will be fairly smooth along the refractor.  The GRM procedure tends 
to smooth out large dips, and therefore the tangent method should never 
encounter adjacent variations in radii that are large enough to cause a 
significant error of this type. 

\section*{Input and Output Specifications}

Input files for GRM are 4-column ASCII files, where the first two columns
are the X and Y coordinates of the receivers, where X has been defined 
with the forward shotpoint at X=0.  The last two columns are forward
and reverse times, respectively.  The method can only calculate depths
for midpoints with both times, so any locations with only one time will
yield a result that will be erroneous and not necessarily zero, so do
not include them.  The number of data records in the file must be 
specified with the required parameter $\texttt{nt}$.

There are two variations of the 1D array that can be handled.  First, an
array where there are receivers spanning the entire distance
between the two shotpoints, with both forward and reverse times at 
all of them.  Here the first and last input data records will be 
equal and give the shotpoint-to-shotpoint time, $t_{ab}$.  This 
configuration is specified by setting the required parameter
$\texttt{abtime}$ to zero.

The second variation is that where receivers with 
both arrival times only span a portion of the distance between the
shotpoints.  In this case, the time $t_{ab}$ must be specified
explicitly, by setting the required parameter $\texttt{abtime}$ to
$t_{ab}$.
This is perhaps the more common situation.

Currently, the receiver spacing is set with the required parameter
$\texttt{dx}$, and must be constant.  A future version will 
include variable station spacing.
Velocity in the surface layer must be specified with the required 
parameter $\texttt{v0}$.

As previously mentioned, the automatic determination of the 
optimum XY can be overridden by setting the optional parameter
$\texttt{xy}$ to the desired offset value.  GRM will still calculate
an optimum and will provide it in the output for reference.

Also mentioned previously is the optional parameter $\texttt{xymax}$
which limits the range of offsets used in the determination
of the optimum XY.  If not specified, the default is 20$dx$.

Finally, the optional parameter $\texttt{depthres}$ sets the increment
of the search for the tangent line in the depth determination step,
as discussed.

The output file generated by GRM is ASCII, and is partitioned into a
header and two
lists of depths.  The first lines of the file will give information
on XY and the average apparent refractor velocity determined by the 
program.  Next will be the list of depths determined by GRM, with 
columns giving the location of the receiver $(x(i),y(i))$, the GRM depth
below the surface $(d_{grm}(i))$, and
an absolute depth $(y(i)-d_{grm}(i))$.
The last section will be a list of depths determined by the 
vertical depth estimator, with $(x(i),y(i))$, the vertical depth estimate
below the receiver $(z(i))$, and absolute depth $(y(i)-z(i))$.

\section*{References}

$\textit{   }$

Hagedoorn, J. G., 1959, The plus-minus method of interpreting seismic
refraction sections: Geophys.
Prosp., 7, no. 2, 158-182. 

$\textit{   }$

Palmer, D., 1982, The Generalized Reciprocal Method of Seismic Refraction
Interpretation: SEG, Tulsa. 

\bibliography{references1}

\bibliographystyle{authordate1}

\end{document}