paper.tex

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

\end{equation}
It is easy to show algebraically that estimate~(\ref{eqn:dinv1}) 
is equivalent to estimate~(\ref{eqn:minv1}) under the condition
\begin{equation}
\mathbf{C} = \mathbf{P}\,\mathbf{P}^T = \left(\mathbf{D}^T\,\mathbf{D}\right)^{-1}\;,
\label{eqn:cond}
\end{equation}
where we need to assume the self-adjoint operator
$\mathbf{D}^T\,\mathbf{D}$ to be invertible.  

To prove the equivalence, consider the operator
\begin{equation}
\mathbf{G} = \mathbf{L}^T\,\mathbf{L}\,\mathbf{C}\,\mathbf{L}^T + \epsilon^2\,\mathbf{L}^T\;,
\label{eq:g}
\end{equation}
which is a mapping from the data space to the model space.
Grouping the multiplicative factors in two different ways, we can
obtain the equality
\begin{equation}
\mathbf{G} = \mathbf{L}^T\,\left(\mathbf{L}\,\mathbf{C}\,\mathbf{L}^T + 
\epsilon^2\,\mathbf{I}\right) = \left(\mathbf{L}^T\,\mathbf{L} + 
\epsilon^2 \mathbf{C}^{-1}\right)\,\mathbf{C}\,\mathbf{L}^T\;,
\label{eq:gmd}
\end{equation}
or, in another form,
\begin{equation}
\mathbf{C}\,\mathbf{L}^T\,\left(\mathbf{L}\,\mathbf{C}\,\mathbf{L}^T + 
\epsilon^2\,\mathbf{I}\right)^{-1} \equiv \left(\mathbf{L}^T\,\mathbf{L} + 
\epsilon^2\,\mathbf{C}^{-1}\right)^{-1}\,\mathbf{L}^T\;.
\label{eq:md}
\end{equation}
The left-hand side of equality (\ref{eq:md}) is exactly the projection
operator from formula (\ref{eqn:dinv1}), and the right-hand side is the
operator from formula (\ref{eqn:minv1}).

This proves the legitimacy
of the alternative data-space approach to data regularization: the
model estimation is reduced to a least-squares minimization of the
specially constructed compound model $\hat{\mathbf{p}}$ under the
constraint~(\ref{eqn:r}).

%\begin{latexonly}
We summarize the differences between the model-space and data-space
regularization in Table 1.

\tabl{table}{Comparison between model-space and data-space regularization}
{\centering
\begin{tabular}{|l||p{2in}|p{2in}|} \hline
  \multicolumn{3}{c}{\rule[-3mm]{0mm}{8mm}} \\
  \hline \rule[-3mm]{0mm}{8mm} 
  & \emph{Model-space} & \emph{Data-space} \\ \hline \hline
  \rule[-8mm]{0mm}{17mm} effective model & $\mathbf{m}$ & $\hat{\mathbf{p}} =
  \left[\begin{array}{c} \mathbf{p} \\ \mathbf{r} \end{array}\right]$ \\ \hline
  \rule[-8mm]{0mm}{17mm} effective data & $\hat{\mathbf{d}} =
  \left[\begin{array}{c} \mathbf{d} \\ \mathbf{0}
    \end{array}\right]$ & $\mathbf{d}$ \\ \hline \rule[-8mm]{0mm}{17mm}
  effective operator & $\mathbf{G_m} = \left[\begin{array}{c} \mathbf{L} \\
      \epsilon \mathbf{D} \end{array}\right]$ & $\mathbf{G_d} =
  \left[\begin{array}{cc} \mathbf{LP} & \epsilon \mathbf{I} \end{array}\right]$ \\
  \hline \rule[-3mm]{0mm}{8mm} optimization problem & minimize
  $\hat{\mathbf{r}}^T \hat{\mathbf{r}}$, \newline where \newline
  \rule[-3mm]{0mm}{8mm} $\hat{\mathbf{r}} = \hat{\mathbf{d}} - \mathbf{G_m m}$ & 
  minimize
  $\hat{\mathbf{p}}^T \hat{\mathbf{p}}$ \newline under the constraint
  \newline \rule[-3mm]{0mm}{8mm} $\mathbf{G_d} \hat{\mathbf{p}} = \mathbf{d}$ \\
  \hline \rule[-5mm]{0mm}{11mm} formal estimate for $\mathbf{m}$ &
  $\left(\mathbf{L}^T \mathbf{L} + 
    \epsilon^2 \mathbf{C}^{-1}\right) \mathbf{L}^T \mathbf{d}$, \newline
  \rule[-5mm]{0mm}{11mm} where $\mathbf{C}^{-1} = \mathbf{D}^T \mathbf{D}$ & 
  $\mathbf{C L}^T (\mathbf{L C L}^T
  + \epsilon^2 \mathbf{I})^{-1} \mathbf{d}$, 
  \newline \rule[-5mm]{0mm}{11mm} where $\mathbf{C}
  = \mathbf{P P}^T$. \\ \hline
\end{tabular}
}
%\end{latexonly}

Although the two approaches lead to similar theoretical results, they
behave quite differently in the process of iterative optimization.  In
the next section, we illustrate this fact with many examples and show
that in the case of incomplete optimization, the second
(preconditioning) approach is generally preferable.

\section{One-dimensional synthetic examples}

\inputdir{oned}

\plot{data}{width=6in}{ The input data (right) are
  irregularly spaced samples of a sinusoid (left).}

In the first example, the input data $\mathbf{d}$ were randomly
subsampled (with decreasing density) from a sinusoid (Figure
\ref{fig:data}). The forward operator $\mathbf{L}$ in this case is
linear interpolation. In other words, we seek a regularly sampled
model $\mathbf{m}$ on 200 grid points that could predict the data with a
forward linear interpolation. Sparse irregular distribution of the
input data makes the regularization necessary.  We applied convolution
with the simple $(1,-1)$ difference filter as the operator $\mathbf{D}$
that forces model continuity (the first-order spline). An appropriate
preconditioner $\mathbf{P}$ in this case is the inverse of $\mathbf{D}$,
which amounts to recursive causal integration \cite[]{gee}. Figures
\ref{fig:im1} and \ref{fig:fm1} show the results of inverse
interpolation after exhaustive 300 iterations of the
conjugate-direction method. The results from the model-space and
data-space regularization look similar except for the boundary
conditions outside the data range.  As a result of using the causal
integration for preconditioning, the rightmost part of the model in
the data-space case stays at a constant level instead of decreasing to
zero. If we specifically wanted a zero-value boundary condition, we
could easily implement it by adding a zero-value data point at the
boundary.

\plot{im1}{width=3in}{ Estimation of a continuous
  function by the model-space regularization. The difference operator
  $\mathbf{D}$ is the derivative operator (convolution with $(1,-1)$).}

\plot{fm1}{width=3in}{Estimation of a continuous
  function by the data-space regularization. The preconditioning
  operator $\mathbf{P}$ is causal integration.}

The model preconditioning approach provides a much faster rate of convergence.
We measured the rate of convergence using the power of the model residual,
defined as the least-squares distance from the current model to the final
solution.  Figure \ref{fig:schwab1} shows that the preconditioning (data
regularization) method converged to the final solution in about 6 times fewer
iterations than the model regularization. Since the cost of each iteration for
both methods is roughly equal, the computational economy is evident.  Figure
\ref{fig:early1} shows the final solution, and the estimates from model- and
data-space regularization after only 5 iterations of conjugate directions. The
data-space estimate appears much closer to the final solution.

\plot{schwab1}{width=3in}{Convergence of the iterative
  optimization, measured in terms of the model residual.  The ``d''
  points stand for data-space regularization; the ``m'' points for
  model-space regularization.  }

\plot{early1}{width=4in}{ The top figure is the exact
  solution found in 250 iterations. The middle is with data-space
  regularization after 5 iterations. The bottom is with model-space
  regularization after 5 iterations.  }

Changing the preconditioning operator changes the regularization
result. Figure \ref{fig:fm6} shows the result of data-space
regularization after a triangle smoother is applied as the model
preconditioner. Triangle smoother is a filter with the $Z$-transform
$\frac{\left(1-Z^N\right)\left(1-Z^{-N}\right)}{(1-Z)\left(1-Z^{-1}\right)}$
\cite[]{Claerbout.blackwell.92}. We chose the filter length $N=6$.

\sideplot{fm6}{width=3in}{Estimation of a smooth function
  by the data-space regularization. The preconditioning operator
  $\mathbf{P}$ is a triangle smoother.}

If, instead of looking for a smooth interpolation, we want to limit the number
of frequency components, then the optimal choice for the model-space
regularization operator $\mathbf{D}$ is a prediction-error filter (PEF). To
obtain a mono-frequency output, we can use a three-point PEF, which has the
$Z$-transform representation $D (Z) = 1 + a_1 Z + a_2 Z^2$. In this case, the
corresponding preconditioner $\mathbf{P}$ can be the three-point
\emph{recursive} filter $P (Z) = 1 / (1 + a_1 Z + a_2 Z^2)$. To test this
idea, we estimated the PEF $D (Z)$ from the output of inverse linear
interpolation (Figure \ref{fig:fm1}), and ran the data-space regularized
estimation again, substituting the recursive filter $P (Z) = 1/ D(Z)$ in place
of the causal integration.  We repeated this two-step procedure three times to
get a better estimate for the PEF. The result, shown in Figure \ref{fig:pm1},
exhibits the desired mono-frequency output. One can accommodate more frequency
components in the model using a longer filter.

\sideplot{pm1}{width=3in}{Estimation of a mono-frequency
  function by the data-space regularization. The preconditioning
  operator $\mathbf{P}$ is a recursive filter (the inverse of PEF).}

\subsection{Regularization after binning: missing data interpolation}

\inputdir{mis1}

One of the factors affecting the convergence of iterative data
regularization is clustering of data points in the output bins. When
least-squares optimization assigns equal weight to each data point, it
may apply inadequate effort to fit a cluster of data points with
similar values in a particular output bin. To avoid this problem, we
can replace the regularized optimization with a less accurate but more
efficient two-step approach: data binning followed by missing data
interpolation.

Missing data interpolation is a particular case of data
regularization, where the input data are already given on a regular
grid, and we need to reconstruct only the missing values in empty
bins.  The basic principle of missing data interpolation is formulated
as follows \cite[]{Claerbout.blackwell.92}:
\begin{quote}
  A method for restoring missing data is to ensure that the restored
  data, after specified filtering, has minimum energy.
\end{quote}
Mathematically, this principle can be expressed by the simple
equation
\begin{equation}
\mathbf{D m \approx 0}\;,
\label{eqn:basic}
\end{equation}
where $\mathbf{m}$ is the model vector and $\mathbf{D}$ is the specified
filter.  Equation~(\ref{eqn:basic}) is completely equivalent to
equation~(\ref{eqn:mreg}).  The approximate equality sign means that
equation (\ref{eqn:basic}) is solved by minimizing the squared norm
(the power) of its left side.  Additionally, the known data values
must be preserved in the optimization scheme. Introducing the mask
operator $\mathbf{K}$, which is a diagonal matrix with
zeros at the missing data locations and ones elsewhere, we can rewrite
equation (\ref{eqn:basic}) in the extended form
\begin{equation}
\mathbf{D (I-K) m} \approx - \mathbf{D K m} = - \mathbf{D m}_k\;,
\label{eqn:rigor}
\end{equation}
where $\mathbf{I}$ is the identity operator and $\mathbf{m}_k$
represents the known portion of the data. It is important to note that
equation (\ref{eqn:rigor}) corresponds to the
regularized linear system
\begin{equation}
\left\{\begin{array}{l}
\mathbf{K m} = \mathbf{m}_k\;, \\
\epsilon \mathbf{D m} \approx \mathbf{0}
\end{array}\right.
\label{eqn:reg}
\end{equation}
in the limit of the scaling coefficient $\epsilon$ approaching zero.
System~(\ref{eqn:reg}) is equivalent to
system~(\ref{eqn:main}-\ref{eqn:mreg}) with the masking
operator~$\mathbf{K}$ playing the role of the forward interpolation
operator~$\mathbf{L}$. Setting $\epsilon$ to zero implies associating
more weight on the first equation in (\ref{eqn:reg}) and using the
second equation only to constrain the null space of the solution.
Applying the general theory of data-space regularization from
the previous section, we can immediately transform system
(\ref{eqn:reg}) to the equation
\begin{equation}
\mathbf{K P p} \approx  \mathbf{m}_k\;,
\label{eqn:prec2}
\end{equation}
where $\mathbf{P}$ is a preconditioning operator and $\mathbf{p}$ is the
preconditioning variable, connected with $\mathbf{m}$ by the simple
relationship~(\ref{eqn:prec}).  According to equations~(\ref{eqn:dinv1})
and~(\ref{eqn:minv1}) from the previous section, equations
(\ref{eqn:prec2}) and (\ref{eqn:rigor}) have exactly the same
solutions if condition~(\ref{eqn:cond}) is satisfied.  If $\mathbf{D}$ is
represented by a discrete convolution, the natural choice for
$\mathbf{P}$ is the corresponding deconvolution (inverse recursive
filtering) operator:
\begin{equation}
\mathbf{P}  = \mathbf{D}^{-1}\;.
\label{eqn:decon}
\end{equation}

We illustrate the missing data problem with a simple 1-D synthetic data
test taken from \cite{gee}.  Figure~\ref{fig:mall} shows the
interpolation results of the unpreconditioned technique with two
different filters $\mathbf{D}$. For comparison with the preconditioned scheme, we
changed the boundary convolution conditions from internal to truncated
transient convolution. As in the previous example, the system was
solved with a conjugate-gradient iterative optimization.

\plot{mall}{width=6in,height=4in}{Unpreconditioned
  interpolation with two different regularization filters. Left plot:
  the top shows the input data; the middle, the result of
  interpolation; the bottom, the filter impulse response. The right
  plot shows the convergence process for the first four iterations.}

As depicted on the right side of the figures, the interpolation
process starts with a ``complicated'' model and slowly extends
it until the final result is achieved. 

Preconditioned interpolation behaves differently (Figure
\ref{fig:sall}). At the early iterations, the model is simple. As the
iteration proceeds, new details are added into the model.  After a