prog.tex

国外免费地震资料处理软件包

\title{Guide to RSF programs}
\email{[e-mail]}
\author{Sergey Fomel}
\lefthead{Fomel}
\righthead{RSF programs}

\maketitle

\begin{abstract}

This guide introduces some of the most used RSF programs and illustrates their
usage with examples.

\end{abstract}

\section{Main programs}

The source files for these programs can be found under
\href{http://svn.sourceforge.net/viewcvs.cgi/rsf/trunk/filt/main/}{filt/main}
in the RSF distribution.

\noindent\doublebox{\parbox{\textwidth}{
\input{sfadd}
}}

\texttt{sfadd} is useful for combining (adding, dividing, or
multiplying) several datasets. What if you want to subtract two
datasets? Easy. Use the \texttt{scale} parameter as follows:
\begin{verbatim}
bash$ sfadd data1.rsf data2.rsf scale=1,-1 > diff.rsf
\end{verbatim}
or
\begin{verbatim}
bash$ sfadd < data1.rsf data2.rsf scale=1,-1 > diff.rsf
\end{verbatim}
The same task can be accomplished with the more general \texttt{sfmath} program:
\begin{verbatim}
bash$ sfmath one=data1.rsf two=data2.rsf output='one-two' > diff.rsf
\end{verbatim}
or
\begin{verbatim}
bash$ sfmath < data1.rsf two=data2.rsf output='input-two' > diff.rsf
\end{verbatim}
In both cases, the size and shape of \texttt{data1.rsf} and
\texttt{data2.rsf} hypercubes should be the same, and a warning
message is printed out if the the axis sampling parameters (such as
\texttt{o1} or \texttt{d1}) in these files are different.

\subsubsection{Implementation: \href{http://svn.sourceforge.net/viewcvs.cgi/rsf/trunk/filt/main/add.c?view=markup}{filt/main/add.c}}

The first input file is either in the list or in the standard input.
\moddex{add}{add}{101}{108}{filt/main}

Collect input files in the \texttt{in} array from all command-line
parameters that don't contain an ``\texttt{=}'' sign. The total number
of input files in \texttt{nin}.
\moddex{add}{add}{110}{116}{filt/main}

A helper function \texttt{check\_compat} checks the compatibility of
input files.
\moddex{add}{add}{412}{450}{filt/main}

Finally, we enter the main loop, where input data are getting read
buffer by buffer and combined in the total product depending on the
data type.
\moddex{add}{add}{173}{194}{filt/main}

The data combination program for floating point numbers is
\texttt{add\_float}.  
\moddex{add}{add}{253}{294}{filt/main}

\noindent\doublebox{\parbox{\textwidth}{
\input{sfattr}
}}

\texttt{sfattr} is a useful diagnostic program. It reports certain
statistical values for an RSF dataset: RMS (root-mean-square)
amplitude, mean value, norm value, variance, standard deviation,
maximum and minimum values, number of nonzero samples, and the total
number of samples.

If we denote data values as $d_i$ for $i=0,1,2,\ldots,n$, then the RMS
value is $\sqrt{\frac{1}{n}\,\sum\limits_{i=0}^n d_i^2}$, the mean
value is $\frac{1}{n}\,\sum\limits_{i=0}^n d_i$, the $L_2$-norm value
is $\sqrt{\sum\limits_{i=0}^n d_i^2}$, the variance is
$\frac{1}{n-1}\,\left[\sum\limits_{i=0}^n d_i^2 - \frac{1}{n}\left(\sum\limits_{i=0}^n d_i\right)^2\right]$, and the standard
deviation is the square root of the variance. Using \texttt{sfattr}
is a quick way to see the distribution of data values and check it for
anomalies.

\subsubsection{Implementation: \href{http://svn.sourceforge.net/viewcvs.cgi/rsf/trunk/filt/main/attr.c?view=markup}{filt/main/attr.c}}

Computations start by finding the input data (\texttt{in}) size
(\texttt{nsiz}) and dimensions (\texttt{dim}).
\moddex{attr}{attr}{72}{75}{filt/main}

In the main loop, we read the input data buffer by buffer.
\moddex{attr}{attr}{89}{108}{filt/main}

The data attributes are accumulated in corresponding double-precision
variables. 
\moddex{attr}{attr}{129}{131}{filt/main}

Finally, the attributes are reduced and printed out.
\moddex{attr}{attr}{161}{166}{filt/main}
\moddex{attr}{attr}{173}{182}{filt/main}

\noindent\doublebox{\parbox{\textwidth}{
\input{sfcat}
}}

\texttt{sfcat} and \texttt{sfmerge} concatenate two or more files
together along a particular axis. It is the same program, only
\texttt{sfcat} has the default \texttt{space=n} and \texttt{sfmerge}
has the default \texttt{space=y}.

Example of \texttt{sfcat}:
\begin{verbatim}
bash$ sfspike n1=2 n2=3 > one.rsf
bash$ sfin one.rsf
one.rsf:
    in="/tmp/one.rsf@"
    esize=4 type=float form=native
    n1=2           d1=0.004       o1=0          label1="Time" unit1="s"
    n2=3           d2=0.1         o2=0          label2="Distance" unit2="km"
        6 elements 24 bytes
bash$ sfcat one.rsf one.rsf axis=1 > two.rsf
bash$ sfin two.rsf
two.rsf:
    in="/tmp/two.rsf@"
    esize=4 type=float form=native
    n1=4           d1=0.004       o1=0          label1="Time" unit1="s"
    n2=3           d2=0.1         o2=0          label2="Distance" unit2="km"
        12 elements 48 bytes
\end{verbatim}

Example of \texttt{sfmerge}:
\begin{verbatim}
bash$ sfmerge one.rsf one.rsf axis=2 > two.rsf
bash$ sfin two.rsf
two.rsf:
    in="/tmp/two.rsf@"
    esize=4 type=float form=native
    n1=2           d1=0.004       o1=0          label1="Time" unit1="s"
    n2=7           d2=0.1         o2=0          label2="Distance" unit2="km"
        14 elements 56 bytes
\end{verbatim}
In this case, an extra empty trace is inserted between the two merged files.

The axes that are not being merged are checked for consistency:
\begin{verbatim}
bash$ sfcat one.rsf two.rsf > three.rsf
sfcat: n2 mismatch: need 3
\end{verbatim}

\noindent\doublebox{\parbox{\textwidth}{
\input{sfcmplx}
}}

\texttt{sfcmplx} simply creates a complex dataset from its real and
imaginary parts. The reverse operation can be accomplished with
\texttt{sfreal} and \texttt{sfimag}.

Example of \texttt{sfcmplx}:
\begin{verbatim}
bash$ sfspike n1=2 n2=3 > one.rsf
bash$ sfin one.rsf
one.rsf:
    in="/tmp/one.rsf@"
    esize=4 type=float form=native
    n1=2           d1=0.004       o1=0          label1="Time" unit1="s"
    n2=3           d2=0.1         o2=0          label2="Distance" unit2="km"
        6 elements 24 bytes
bash$ sfcmplx one.rsf one.rsf > cmplx.rsf
bash$ sfin cmplx.rsf
cmplx.rsf:
    in="/tmp/cmplx.rsf@"
    esize=8 type=complex form=native
    n1=2           d1=0.004       o1=0          label1="Time" unit1="s"
    n2=3           d2=0.1         o2=0          label2="Distance" unit2="km"
        6 elements 48 bytes
\end{verbatim}

\noindent\doublebox{\parbox{\textwidth}{
    \input{sfconjgrad}
}}

\texttt{sfconjgrad} is a generic program for least-squares linear
inversion with the conjugate-gradient method. Suppose you have an
executable program \texttt{<prog>} that takes an RSF file from the
standard input and produces an RSF file in the standard output. It may
take any number of additional parameters but one of them must be
\texttt{adj=} that sets the forward (\texttt{adj=0}) or adjoint
(\texttt{adj=1}) operations.  The program \texttt{<prog>} is typically
an RSF program but it could be anything (a script, a multiprocessor
MPI program, etc.) as long as it implements a linear operator
$\mathbf{L}$ and its adjoint. There are no restrictions on the data
size or shape. You can easily test the adjointness with
\texttt{sfdottest}. The \texttt{sfconjgrad} program searches for a
vector $\mathbf{m}$ that minimizes the least-square misfit 
$\|\mathbf{d - L\,m}\|^2$ for the given input data vector $\mathbf{d}$.

Here is an example. The \texttt{sfhelicon}
program implements Claerbout's multidimensional helical filtering
\cite[]{GEO63-05-15321541}. It requires a filter to be specified in
addition to the input and output vectors. We create a helical 
2-D filter using the Unix \texttt{echo} command.
\begin{verbatim}
bash$ echo 1 19 20 n1=3 n=20,20 data_format=ascii_int in=lag.rsf > lag.rsf
bash$ echo 1 1 1 a0=-3 n1=3 data_format=ascii_float in=flt.rsf > flt.rsf
\end{verbatim}
Next, we create an example 2-D model and data vector with \texttt{sfspike}.
\begin{verbatim}
bash$ sfspike n1=50 n2=50 > vec.rsf
\end{verbatim}
The \texttt{sfdottest} program can perform the dot product test to
check that the adjoint mode works correctly.
\begin{verbatim}
bash$ sfdottest sfhelicon filt=flt.rsf lag=lag.rsf \
> mod=vec.rsf dat=vec.rsf
sfdottest:  L[m]*d=5.28394
sfdottest: L'[d]*m=5.28394
\end{verbatim}
Your numbers may be different because \texttt{sfdottest} generates new
random input on each run.
Next, let us make some random data with \texttt{sfnoise}.
\begin{verbatim}
bash$ sfnoise seed=2005 rep=y < vec.rsf > dat.rsf
\end{verbatim}
and try to invert the filtering operation using \texttt{sfconjgrad}:
\begin{verbatim}
bash$ sfconjgrad sfhelicon filt=flt.rsf lag=lag.rsf \
mod=vec.rsf < dat.rsf > mod.rsf niter=10
sfconjgrad: iter 1 of 10
sfconjgrad: grad=3253.65
sfconjgrad: iter 2 of 10
sfconjgrad: grad=289.421
sfconjgrad: iter 3 of 10
sfconjgrad: grad=92.3481
sfconjgrad: iter 4 of 10
sfconjgrad: grad=36.9417
sfconjgrad: iter 5 of 10
sfconjgrad: grad=18.7228
sfconjgrad: iter 6 of 10
sfconjgrad: grad=11.1794
sfconjgrad: iter 7 of 10
sfconjgrad: grad=7.26941
sfconjgrad: iter 8 of 10
sfconjgrad: grad=5.15945
sfconjgrad: iter 9 of 10
sfconjgrad: grad=4.23055
sfconjgrad: iter 10 of 10
sfconjgrad: grad=3.57495
\end{verbatim}
The output shows that, in 10 iterations, the norm of the gradient vector decreases by almost 1000. 
We can check the residual misfit before
\begin{verbatim}
bash$ < dat.rsf sfattr want=norm
norm value = 49.7801
\end{verbatim}
and after
\begin{verbatim}
bash$ sfhelicon filt=flt.rsf lag=lag.rsf < mod.rsf | \
sfadd scale=1,-1 dat.rsf | sfattr want=norm
norm value = 5.73563
\end{verbatim}
In 10 iterations, the misfit decreased by an order of magnitude. The
result can be improved by running the program for more iterations.

\noindent\doublebox{\parbox{\textwidth}{
\input{sfcp}
}}

The \texttt{sfcp} and \texttt{sfmv} command imitate the Unix
\texttt{cp} and \texttt{mv} commands and serve for copying and moving
RSF files. Example:
\begin{verbatim}
bash$ sfspike n1=2 n2=3 > one.rsf
bash$ sfin one.rsf
one.rsf:
    in="/tmp/one.rsf@"
    esize=4 type=float form=native
    n1=2           d1=0.004       o1=0          label1="Time" unit1="s"
    n2=3           d2=0.1         o2=0          label2="Distance" unit2="km"
        6 elements 24 bytes
bash$ sfcp one.rsf two.rsf
bash$ sfin two.rsf
two.rsf:
    in="/tmp/two.rsf@"
    esize=4 type=float form=native
    n1=2           d1=0.004       o1=0          label1="Time" unit1="s"
    n2=3           d2=0.1         o2=0          label2="Distance" unit2="km"
        6 elements 24 bytes
\end{verbatim}