Next: Interpolation example
Up: Fomel: Splines in tension
Previous: Mathematical theory of splines
In the onedimensional case, a finitedifference representation of the
squared Laplacian can be defined as a centered 5point filter with
coefficients (1,4,6,4,1). On the same grid, the Laplacian operator
can be approximated to the same order of accuracy with the filter
(1/12,4/3,5/2,4/3,1/12). Combining the two filters in accordance
with equation (3) and performing a spectral factorization
with one of the standard methods
Claerbout (1976, 1992), we can obtain a
3point minimumphase filter, suitable for inverse filtering.
Figure 1 shows a family of onedimensional minimumphase
filters for different values of the parameter t.
Figure 2 demonstrates the interpolation results obtained
with these filters on a simple onedimensional synthetic. As expected,
a small tension value (t=0.01) produces a smooth interpolation, but
creates artificial oscillations in the unconstrained regions around
sharp changes in the gradient. The value of t=1 leads to linear
interpolation with no extraneous inflections, but with discontinuous
derivatives. Intermediate values of t allow us to achieve a
compromise: a smooth surface with constrained oscillations.
otens
Figure 1 Onedimensional minimumphase
filters for different values of the tension parameter t. The
filters range from the second derivative for t=0 to the first
derivative for t=1.

 
int
Figure 2 Interpolating a simple
onedimensional synthetic with recursive filter preconditioning for
different values of the tension parameter t. The input data is
shown on the top. The interpolation results range from a natural
cubic spline interpolation for t=0 to linear interpolation for
t=1.
To design the corresponding filters in two dimensions, I define the
finitedifference representation of operator (3) on a
5by5 stencil. The filters coefficients are chosen with the help of
the Taylor expansion to match the desired spectrum of the operator
around the zero spatial frequency. The matching conditions lead to
the following set of coefficients for the squared Laplacian:
1/60 
2/5 
7/30 
2/5 
1/60 
2/5 
14/15 
44/15 
14/15 
2/5 
7/30 
44/15 
57/5 
44/15 
7/30 
2/5 
14/15 
44/15 
14/15 
2/5 
1/60 
2/5 
7/30 
2/5 
1/60 
= 1/60
1 
24 
14 
24 
1 
24 
56 
176 
56 
24 
14 
176 
684 
176 
14 
24 
56 
176 
56 
24 
1 
24 
14 
24 
1 
Laplacian representation with the same order of accuracy has the
coefficients
1/360 
2/45 

2/45 
1/360 
2/45 
14/45 
4/5 
14/45 
2/45 

4/5 
41/10 
4/5 

2/45 
14/45 
4/5 
14/45 
2/45 
1/360 
2/45 

2/45 
1/360 
=
1/360
1 
16 

16 
1 
16 
112 
288 
112 
16 

288 
1476 
288 

16 
112 
288 
112 
16 
1 
16 

16 
1 
For the sake of simplicity, I assumed an equal physical spacing in x
and y directions. The coefficients can be easily adjusted for
anisotropic spacing. Figures 3 and 4 show
the spectra of the finitedifference representations of
operator (3) for the different values of the tension
parameter. The finitedifferent spectra appear as fairly isotropic.
They match the exact expressions at small frequencies.
specc
Figure 3 Spectra of the finitedifference
splinesintension schemes for different values of the tension
parameter (contour plots).
specp
Figure 4 Spectra of the finitedifference
splinesintension schemes for different values of the tension
parameter (crosssection plots). Dashed lines show the exact spectra
for continuous operators.
Regarding the finitedifference operators as twodimensional
autocorrelations and applying the efficient WilsonBurg method of
spectral factorization Claerbout (1999); Sava et al. (1998), I obtain
twodimensional minimumphase filters suitable for inverse filtering.
The exact filters contain many coefficients, which rapidly decrease in
magnitude at a distance from the first coefficient. For reasons of
efficiency, it is advisable to restrict the shape of the filter so
that it contains only the valuable coefficients. Keeping all the
coefficients that are 1000 times smaller in magnitude than the
leading coefficient creates a 53point filter for t=0 and a 35point
filter for t=1, with intermediate filter lengths for intermediate
values of t. When the ratio is changed to 200, we obtain 25 and
16point filters, respectively. The restricted filters don't factor
the autocorrelation exactly, but provide an effective approximation of
the exact factors. As outputs of the WilsonBurg spectral
factorization process, they obey the minimumphase condition.
splin
Figure 5 Inverse filtering with the tension
filters. The left plots show the inputs composed of filters and
spikes. Inverse filtering turns filters into impulses and turns
spikes into inverse filter responses (middle plots). Adjoint
filtering creates smooth isotropic shapes (right plots). The tension
parameter takes values 0, 0.3, 0.7, and 1 (from top to bottom).
Figure 5 shows the twodimensional filters for different
values of t and illustrates inverse recursive filtering, which is
the essence of the helix method
Claerbout (1998a,b, 1999). The case of t=1
leads to the filter known as helix derivative
Claerbout (1999); Zhao (1999). The filter values are spread mostly on
two columns. The other boundary case of t=0 leads to a threecolumn
filter, which serves as the minimumphase version of the Laplacian. As
expected from the theory, the inverse impulse response of this filter
is noticeably smoother and wider than the inverse response of the
helix derivative. Filters corresponding to intermediate values of t
exhibit intermediate properties. Theoretically, the inverse impulse
response of the filter corresponds to the Green function of
equation (3). The theoretical Green function for the case
of t=1 is
where r is the distance from the impulse:
. In the case of
t=0, the Green function is smoother at the origin:
 
(6) 
The theoretical Green function expression for an arbitrary value of
t is not known, but we can assume that its smoothness lies between
the two boundary conditions.
In the next section, I illustrate an application of helical inverse
filtering to a twodimensional interpolation problem.
Next: Interpolation example
Up: Fomel: Splines in tension
Previous: Mathematical theory of splines
Stanford Exploration Project
4/27/2000