An optimality formulation for the design of mixed

An optimality formulation for the design of mixed-phase inverse filters
Milton J. Porsani∗ , CPGG/UFBA, and Bjørn Ursin, Dept. of Petroleum Engineering and Applied
Geophysics/Norwegian Univ. of Sc. and Tech.
Summary
We describe a new algorithm for the design of mixedphase filters. Solving the extended Yule-Walker
(EYW) equation, for lag α of the auto-correlation
function (ACF) on diagonal of the matrix, and
multiplying the Z-transform of the solution by the
minimum-delay wavelet, a finite polynomial of order α can be obtained. Its α roots correspond to
an estimate of roots for the minimum-delay wavelet.
By systematically combining and flipping these roots
from outside to inside the unit circle, one can generate mixed-delay wavelets and their corresponding
mixed-phase inverse filters. The optimal inverse filter is the one which gives the maximum value of
the Lp-norm of the filtered signal. Deconvolution of
synthetic and real data examples shows the performance of the new method.
Introduction
Inverse filtering or spiking deconvolution is routinely
applied to seismic data to remove the effect of the
wavelet (the source pulse convolved with earth filtering effects) in order to obtain the reflectivity series
(Robinson, 1957; Robinson and Treitel, 1980; Yilmaz, 1987). This produces good results in minimum
delay wavelets. However, in the case of a mixeddelay wavelet, the result is rather poor (Leinbach,
1995). If the wavelet is known, or can be reliably
estimated, an optimal pulse-shaping filter may be
applied to the data (Treitel and Robinson, 1964;
Berkhout, 1977). By delaying the desired output
pulse an additional phase-shift is introduced, resulting in a good filter performance.
Eisner and Hampson (1990) have presented an algorithm for decomposing a mixed-delay wavelet into its
minimum-delay and maximum-delay components.
This could be used in inverse filter design by convolving the causal and anti-causal inverse filters of
the two parts of the wavelet, resulting in an anticausal inverse filter.
In many cases the wavelet is not known and only
an estimate of its auto-correlation function (ACF)
is available. For minimum-delay wavelets, the YuleWalker (YW) equations can be solved to give the
minimum-phase inverse filter. By computing the
ACF of the inverse filter and solving the new
YW equations one can obtain an estimate of the
minimum-delay wavelet.
Porsani and Ursin (1996, 1998) solved the extended
Yule-Walker (EYW) equations, where the lag α of
the ACF is on the main diagonal of the filter equa-
tions, to decompose the minimum-delay wavelet into
a finite polynomial of order α multiplied by a polynomial of infinite order. Then they consider these
polynomials as the maximum- and minimum-delay
component of the wavelet, respectively. Its inverse
generates the causal and anti-causal components of
the mixed-phase inverse filter. For different values of
the lag α (α = 1, 2, . . .) the mixed-phase filter was
applied to seismic data, and the Lp-norm (p > 2)
of the output was computed. The filter which gave
the largest Lp-norm was selected as optimal. The
result was insensitive to the choice of norm, as long
as p > 2.
We present here a generalization of the previous procedure. For each value of α we compute the roots of
the Z-transform polynomial of the short filter. Each
root (a pair of complex conjugate roots is treated as
one root) may give rise to a maximum-phase component of the inverse filter. By systematically selecting
all possible combinations of roots and flipping them
to inside the unit circle, we generate, at most, 2α
different filters. One of these is the minimum-phase
inverse filter and one is the mixed-phase filter, described in our previous paper. Again we select the
filter which maximizes the Lp-norm of the filtered
output with p = 5. We only give a brief overview of
mixed-phase inverse filtering and the method. All
algorithmic details can be found in our previous paper (Porsani and Ursin, 1998).
In order to investigate the performance of the new algorithm, we first apply it to a synthetic mixed-delay
pulse and then to deconvolve a marine seismic section. Finally, we show its application to improve the
resolution of the semblance plots for velocity analysis.
Methodology
Porsani and Ursin (1998) have shown that mixedphase inverse filters, can be computed by solving the
EYW equation, for different lags of the ACF on diagonal of the matrix, and choosing the filter which
performs best when applied to seismic data. Here
we start with the same decomposition of the mixedphase inverse filter, but the phase of the all-pass
filter is varied by combining and flipping roots from
outside to inside the unit circle, in a systematic way,
providing a much larger number of possible inverse
filters. And again, we choose the filter which performs best when applied to seismic data. Starting
with Rp (τ ), an estimate of the ACF of the pulse, we
˜ t = c0 ,
compute the minimum-phase inverse filter, h
t
by solving the YW equations. This is the optimal
Mixed-phase inverse filter
filter in the sense of a minimum-delay pulse. Next
we compute Rh (τ ), the ACF of the inverse filter, and
compute the minimum-delay pulse p˜t by solving the
YW equations using Rh (τ ).
As in Porsani and Ursin (1998) we solve the EYW
equations with α = 1, 2, . . . to obtain the series cα
t.
Then we compute the finite, minimum-delay wavelet
bα
˜t ∗ cα
t =p
t.
(1)
Factoring its Z-transform results,
α α
B α (Z) = 1 + bα
1 Z + . . . + bα Z =
α
Y
j=1
(1 −
Z
) , (2)
rj
where the complex roots always occur in complex
conjugate pairs. The roots of B α (Z) are outside
the unit circle, so that |rj | > 1. In our previous
approach (Porsani and Ursin; 1998) B α (Z) was used
directly in the equation,
e
H(Z) = H(Z)
B α (Z)
.
Z α B α (Z −1 )
Once the optimal filter has been computed, we can
obtain an estimate of the wavelet form by taking the
inverse of equation (5).
Examples
The new inverse algorithm was tested on a synthetic
mixed-delay wavelet with two real and one complex
conjugate pair of roots inside the unit circle, and two
real and three complex-conjugate pair of roots outside the unit circle, as given in Table 1. The results
of applying the optimal inverse filters, for different
values of α, are shown in Figure 1. We remark that
both filters have same amplitude spectrum. It is seen
(also confirmed numerically) that α = 8 gives a perfect inverse filter, much better than the minimumphase inverse filter (α = 0). The estimated roots
which were used to generate the mixed-phase inverse
filters are given in Table 1. It is seen that the algorithm first selects the roots closest to the unit circle.
(3)
Now we vary the phase of the inverse filter systematically by only using some of the roots of B α (Z).
This is done by forming
B β (Z) =
β
Y
(1 −
j=1
Z
)
rj
(4)
and applying the filter
e
H(Z) = H(Z)
B β (Z)
Z β B β (Z −1 )
(5)
to the data. B β (Z) corresponds to a maximumdelay component of the pulse with roots rj−1 , j =
1, . . . , β. By varying the roots rj and the number of
roots β, we obtain, for a given value of α, at most 2α
different filters. These may be less than 2α filters,
because a pair of complex conjugate roots is treated
as one root (they are either both inside, or both outside the unit circle). β = 0 gives the minimum-phase
filter and β = α gives the mixed-phase filter of Porsani and Ursin (1998). For each filter we compute
the performance criterion,
P α,β 5
|e |
Φ(α, β) = Pt t0,0
(6)
5
t |et |
where eα,β
is the deconvolved trace for given values
t
of α and β. For each value of α only the best filter
is saved, and then the optimal filter is chosen by
selecting the α which gives the largest value of Φ.
Figure 1: Inverse filtering. Minimum-phase inversion
(a). Mixed-phase inversion for α = 1, . . . , 8, (b)-. . . -(i).
The next example deals with the prestack deconvolution of marine seismic data recorded with a 4 ms
sampling interval. A seismic section consisting of 50
CMP gathers with 24 traces each was deconvolved
and stacked using the same velocity function.
The optimal mixed- and minimum-phase filters,
shown in (b) and (a) in Figure 2, were obtained
from the averaged ACF of a CMP gather selected
inside the seismic section shown in Figure 3. These
filters were applied to all traces before stacking. The
stacked section of the original data without deconvolution is shown in Figure 3a. Figures 3b and 3c
shows the stacked sections with prestack minimumand mixed-phase deconvolution, using the filters in
Figures 2a and 2b, respectively. Clearly the mixedphase filter performs better than the minimumphase does.
Mixed-phase inverse filter
ESTIMATED ROOTS OF THE MINIMUM-PHASE WAVELET
α=1 α=2
α=3
α=4
α=5
α=6
α=7
α=8
1.197 1.053 1.111
1.111
1.111
1.111
1.111
1.111
-1.441 1.100±130◦ 1.100±130◦ 1.100±130◦ 1.100±130◦ 1. 100±130◦ 1.100±130◦
1.792
1.250±5.31◦ 1.251±5.4◦ 1.248±5.08◦ 1.248±5.02◦
1.463
1.313± 29.2◦ 1.337±30◦
1.479
ACTUAL
ROOTS
1/1.111
1.100± 130◦
1.25± 5◦
1/1.333 ± 30◦
1/1.428
1.3
-1.5
1.7 ± 40◦
Table 1: Estimated roots associated to the finite polinomial, B α (Z) for different values of α.
Acknowledgements
The authors wish to express their gratitude to the
geophysicist Adriano de Pinho Lima/PETROBRAS
for processing the marine seismic data presented.
To Petroleum Geo-Services and to the Conselho Nacional de Desenvolvimento Cient´ıfico e Tecnol´ogico,
(CNPq), Brazil, for their financial support.
References
Berkhout, A. J., 1977, Least squares inverse filtering
and wavelet deconvolution: Geophysics, 42, 13691383.
Figure 2:
Inverse filters:
mixed-phase (b).
Minimum-phase (a) and
Figures 4 and 5 show the resolution improvement
obtained in the semblance plot for the velocities by
using an optimal mixed-phase inverse filter. The
mixed-phase deconvolution, in Figure 5a, shows better continuity and resolution of the reflected events
as compared with the minimum-phase results in Figure 4a. The semblance plots using the data in Figures 4a and 5a are shown in Figures 4b and 5b. The
velocities are better defined in the semblance plot
for the mixed-phase deconvolved CMP.
Conclusions
We have described a new algorithm for mixed-phase
deconvolution. By using an estimate of the ACF
of the pulse and solving EYW equations, it is possible to estimate the roots of the minimum-delay
wavelet. These roots may be used to form all-pass
filters in order to change the phase characteristics of
the minimum-phase filter. The optimal inverse filter
was found by maximizing the Lp-norm of the filtered
seismic section. Synthetic and real data examples
demonstrated the effectiveness of the method.
Eisner, E., and Hampson, G., 1990, Decomposition into minimum and maximum-phase components: Geophysics, 55, 897-901.
Leinbach, J., 1995, Wiener spiking deconvolution
and minimum-phase wavelets: A tutorial: The
Leading Edge, 189-192.
Porsani, M. J., and Ursin, B., 1996, Mixed-phase
deconvolution of seismic and GPR data: 66th Ann.
Internat. Mtg., SEG, Expanded Abstracts, 1603–
1606.
Porsani, M. J., and Ursin, B., 1998, Mixed-phase
deconvolution, Geophysics, in press.
Robinson, E. A., and Treitel, S., 1980, Geophysical
signal analysis: Prentice-Hall, Englewood Cliffs.
Robinson, E. A., 1957, Predictive decomposition of
seismic traces: Geophysics, 22, 767-778.
Treitel, S., and Robinson, E. A., 1964, The stability
of digital filters: IEEE Trans. Geoscience Electronics, GE-2, 6-18.
¨ 1987, Seismic data processing: SEG,
Yilmaz, O.,
Tulsa.
Mixed-phase inverse filter
25
1.0
50
25
50
25
50
Time (s)
1.5
2.0
2.5
3.0
(a)
(b)
(c)
Figure 3: Stacked sections: without deconvolution (a), after minimum-phase deconvolution (b), and after
mixed-phase deconvolution (c).
1.3
10
20
Velocity (m/s)
1400
1500
1600
1300
1.3
1700
1.3
10
20
1300
1.3
1.5
1.5
1.5
1.5
1.7
1.7
1.7
1.7
Velocity (m/s)
1400
1500
1600
1700
0.6
0.5
1.9
1.9
2.1
2.1
0.5
0.5
0.5
2.3
2.3
0.6
Time (s)
Time (s)
0.6
0.6
0.5
1.9
1.9
2.1
2.1
2.3
2.3
2.5
2.5
0.5
0.5
0.5
0.5
2.5
2.5
2.7
0.6
0.5
0.5
0.5
2.7
2.7
2.7
2.9
2.9
3.1
3.1
0.5
2.9
2.9
3.1
3.1
0.6
0.5
0.5
0.5
0.6
0.5
0.6
0.5
0.5
(a)
0.8
0.7
(b)
0.6
0.5
(a)
0.6
(b)
Figure 4: Minimum-phase deconvolution of a marine
Figure 5: Mixed-phase deconvolution of a marine
CMP gather (a) and semblance plot of the velocity
analysis (b).
CMP gather (a) and semblance plot of the velocity
analysis (b).

Download Report