Streaming prediction-error filters

Sergey Fomel; Jon Claerbout

doi:10.1190/geo2023-0646.1

Skip Nav Destination

Research Article| July 29, 2024 Early Publication

Streaming prediction-error filters

Sergey Fomel;

Sergey Fomel

The University of Texas at Austin

, Jackson School of Geosciences, Bureau of Economic Geology, Austin, Texas, USA. [email protected] (corresponding author)

https://orcid.org/0000-0002-9024-5137

Search for other works by this author on:

GSW

Google Scholar

Jon Claerbout

Stanford University

, Department of Geophysics, Stanford, California, USA. [email protected]

Search for other works by this author on:

GSW

Google Scholar

Author and Article Information

Sergey Fomel

https://orcid.org/0000-0002-9024-5137

The University of Texas at Austin

, Jackson School of Geosciences, Bureau of Economic Geology, Austin, Texas, USA. [email protected] (corresponding author)

Jon Claerbout

Stanford University

, Department of Geophysics, Stanford, California, USA. [email protected]

Publisher: Society of Exploration Geophysicists

Received: 03 Nov 2023

Revision Received: 22 Mar 2024

Accepted: 29 Apr 2024

First Online: 27 Jul 2024

Online ISSN: 1942-2156

Print ISSN: 0016-8033

© 2024 The Authors. Published by the Society of Exploration Geophysicists. All article content, except where otherwise noted (including republished material), is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License (CC BY-SA). See https://creativecommons.org/licenses/by-sa/4.0/. Distribution or reproduction of this work in whole or in part commercially or noncommercially requires full attribution of the original publication, including its digital object identifier (DOI). Derivatives of this work must carry the same license.

The Authors

Freely available online through the SEG open-access option.

Geophysics (2024) F89–F95.

https://doi.org/10.1190/geo2023-0646.1

ABSTRACT

Prediction-error filters (PEFs) are essential for seismic deconvolution and other geophysical estimation problems. We find that nonstationary multidimensional PEFs can be computed in a “streaming” manner, wherein the filter is updated incrementally by accepting one new data point at a time. The computational cost of estimating a streaming PEF is reduced to the cost of a single convolution. In other words, the cost of the PEF design while filtering is equivalent to the cost of applying the filter. Moreover, the nonlinear operation of finding and applying a streaming PEF is invertible at a similar cost, which enables a fast approach for missing data interpolation.

INTRODUCTION

Prediction-error filters (PEFs) play an essential role in geophysical data analysis. They can be used directly in different forms of seismic deconvolution (Webster, 1978; Robinson and Osman, 1996). More importantly, PEFs are a practical approximation for the covariance operators used in geophysical estimation problems (Claerbout, 2014). As such, they provide a way to absorb and characterize statistical patterns (Claerbout and Brown, 1999).

Because most natural patterns are nonstationary, extending the classic stationary formulation to nonstationarity is vital. This extension can be done either by patching, that is, breaking data into overlapping windows, or by a smoothly nonstationary estimation with regularization (Crawley et al., 1999; Curry, 2003). Shaping the regularization provides a particularly effective method for constraining the filter variability (Fomel, 2009; Liu and Fomel, 2011; Liu et al., 2012). Both approaches have an additional computational expense, particularly in storing variable filter coefficients (Ruan et al., 2015).

This paper presents an efficient approach for computing and applying nonstationary PEFs without excessive storage requirements. Using an adaptive-filtering approach (Widrow and Stearns, 1985; Haykin, 2002), we show that a nonstationary PEF can be updated on the fly by accepting one data point at a time. The cost of computing and applying such a PEF is reduced to the cost of a single convolution. In other words, the cost of the PEF design while filtering is equivalent to the cost of applying the filter. Moreover, the nonlinear operation of finding and applying a PEF has an exact inverse, which finds an immediate application in missing data interpolation problems. We extend the filter from one dimension to multiple dimensions using the helix transform (Claerbout, 1998) and show its application to simple benchmark problems.

We have implemented the proposed algorithm in the Julia programming language (Bezanson et al., 2017). For clarity, the essential code is included in the paper. All presented examples are reproducible using the Jupyter Notebook (Kluyver et al., 2016; Granger and Pérez, 2021) provided with the paper.

Theory

The basic equation defining a PEF is autoregression. For a given data stream

d

and a 1D PEF

a

⁠, the autoregression equation takes the form

[\begin{matrix} d_{n + 1} & d_{n} & d_{n - 1} & \dots & d_{n + 1 - k} \end{matrix}] [\begin{matrix} 1 \\ a_{1} \\ a_{2} \\ \dots \\ a_{k} \end{matrix}] \approx 0,

(1)

where it is assumed that

n \geq k

⁠. Equation represents the convolution of the data with the PEF.

Suppose the filter is updated with each new data point $d_{n + 1}$ ⁠. The additional constraint we can impose to control the variability of the filter is that the new filter $a$ stays close to the previous filter $\bar{a}$ ⁠.

The two conditions can be combined into one overdetermined linear system:

[\begin{matrix} d_{n} & d_{n - 1} & \dots & d_{n + 1 - k} \\ λ & 0 & 0 \\ 0 & λ & 0 \\ \dots & \dots & \dots \\ 0 & 0 & \dots & λ \end{matrix}] [\begin{matrix} a_{1} \\ a_{2} \\ \dots \\ a_{k} \end{matrix}] \approx [\begin{matrix} - d_{n + 1} \\ λ {\overline{a}}_{1} \\ λ {\overline{a}}_{2} \\ \dots \\ λ {\overline{a}}_{k} \end{matrix}],

(2)

where

λ

is the parameter that controls how much we allow

a

to deviate from

\bar{a}

⁠. In a shortened block-matrix notation, we can rewrite the linear system of equations as

[\begin{matrix} d^{T} \\ λ I \end{matrix}] a \approx [\begin{matrix} - d_{n + 1} \\ λ \overline{a} \end{matrix}],

(3)

where

d = [\begin{matrix} d_{n} \\ d_{n - 1} \\ \dots \\ d_{n + 1 - k} \end{matrix}], a = [\begin{matrix} a_{1} \\ a_{2} \\ \dots \\ a_{k} \end{matrix}],

(4)

and

I

is the identity matrix.

The least-squares solution of the overdetermined system is

a = (d d^{T} + λ^{2} I)^{- 1} (- d_{n + 1} d + λ^{2} \overline{a}) .

(5)

Next, we can use the Sherman-Morrison formula from linear algebra (Hager, 1989) to transform the inverse matrix as follows:

(d d^{T} + λ^{2} I)^{- 1} = \frac{1}{λ^{2}} (I - \frac{d d^{T}}{λ^{2} + d^{T} d}) .

(6)

Substituting the Sherman-Morrison equation and doing algebraic simplifications lead to the final result:

a = \frac{1}{λ^{2}} (I - \frac{d d^{T}}{λ^{2} + d^{T} d}) (- d_{n + 1} d + λ^{2} \bar{a}) = \overline{a} - (\frac{d_{n + 1} + d^{T} \overline{a}}{λ^{2} + d^{T} d}) d .

(7)

The Sherman-Morrison technique is a well-known approach in adaptive filtering by recursive least squares (Haykin, 2002).

This equation shows that the filter is updated by subtracting a scaled version of the data. The scale is proportional to the convolution residual computed using the previous version of the filter. Updating the filter requires only elementary algebraic operations (vector dot products) and no iterations. Moreover, computing the dot product

d^{T} d

in a streaming fashion requires at most two multiplications:

d^{T} d = {\overline{d}}^{T} \overline{d} + d_{n}^{2} - d_{n - k}^{2},

(8)

where

k

is the number of points in

a

and

d

⁠.

Inversion

The streaming residual

r_{n + 1}

from the left side of equation can be expressed as

r_{n + 1} = d_{n + 1} + d^{T} a,

(9)

or substituting the equation for

a

⁠,

r_{n + 1} = d_{n + 1} + d^{T} [\bar{a} - (\frac{d_{n + 1} + d^{T} \bar{a}}{λ^{2} + d^{T} d}) d] = (d_{n + 1} + d^{T} \bar{a}) (1 - \frac{d^{T} d}{λ^{2} + d^{T} d}) = \frac{λ^{2} (d_{n + 1} + d^{T} \bar{a})}{λ^{2} + d^{T} d} .

(10)

Equation () can be directly inverted to reconstruct

d_{n + 1}

as follows:

d_{n + 1} = r_{n + 1} (1 + \frac{d^{T} d}{λ^{2}}) - d^{T} \bar{a} .

(11)

We see that streaming time-variable deconvolution is invertible: the input data and the time-varying filter can be reconstructed from the streaming residual. Note that because of the division by

λ^{2}

⁠, the inverse operation may become numerically unstable for small

λ

⁠. The forward and inverse filtering is implemented in Program 1.

Multiple dimensions

To extend the filter from one dimension to multiple dimensions, we follow the helix transform of Claerbout (1998). The change is minimal: the multidimensional data are arranged in a 1D stream on a helix and the definition of the data vector

d

changes to

d = [\begin{matrix} d_{n + 1 - l_{1}} \\ d_{n + 1 - l_{2}} \\ \dots \\ d_{n + 1 - l_{k}} \end{matrix}],

(12)

where

l_{1}, l_{2}, \dots, l_{k}

represent the lags of the helical filter. The computational cost and other benefits of streaming PEFs remain the same. The implementation is provided in Programs 2 and 3.

One caveat is that, in multiple dimensions, we may want the nonstationary filter to change smoothly along the helix and other dimensions. As shown in the following, this goal can be accomplished with only a minor change in the algorithm. Suppose that

{\overline{a}}_{1}

is the previous filter in the first (helical) dimension and

{\overline{a}}_{2}

is the previous filter in the second dimension. We can keep the newly estimated streaming PEF close to both filters by changing the matrix equation to

[\begin{matrix} d^{T} \\ λ_{1} I \\ λ_{2} I \end{matrix}] a \approx [\begin{matrix} - d_{n + 1} \\ λ_{1} {\overline{a}}_{1} \\ λ_{2} {\overline{a}}_{2} \end{matrix}],

(13)

where

λ_{1}

and

λ_{2}

control the filter variability in two directions. The least-squares solution of this system is

a = (d d^{T} + λ^{2} I)^{- 1} (- d_{n + 1} d + λ_{1}^{2} {\bar{a}}_{1} + λ_{2}^{2} {\bar{a}}_{2}),

(14)

where

λ^{2} = λ_{1}^{2} + λ_{2}^{2}

⁠. The new equation is equivalent to the previous equation with the simple substitution:

\overline{a} = \frac{λ_{1}^{2} {\overline{a}}_{1} + λ_{2}^{2} {\overline{a}}_{2}}{λ_{1}^{2} + λ_{2}^{2}} .

(15)

Thus, the only change in the algorithm is an increase in the storage to keep track of

{\overline{a}}_{1}

and

{\overline{a}}_{2}

⁠. The extension to more dimensions is analogous.

Cost comparison

Table 1 compares the computational cost of different approaches for computing PEFs. The cost of a streaming PEF is minimal and reduces to the cost of a single convolution. Streaming PEF can capture the input data’s nonstationary character without storing multiple copies of the filter.

Moreover, all streaming computations are local, which makes them particularly suitable for modern hardware accelerators such as general-purpose graphics processing units (Sanders and Kandrot, 2010) or Intel Xeon Phi (Jeffers and Reinders, 2013).

Missing data interpolation

The existence of the exact inversion for the application of a streaming PEF allows for a straightforward approach to missing data interpolation. When reading data in a streaming fashion, every time we meet a missing data point $d_{n + 1}$ ⁠, we can replace its value with a value computed from the residual. The residual value $r_{n + 1}$ ⁠, in this case, can be set to zero or a random number (white noise) with the expected variance of the residual. In the latter case, it is possible to generate multiple equiprobable distributions for the interpolation result (Clapp, 2000). The missing data interpolation is implemented in Program 4.

A simple 1D interpolation example using data from Figure 1 is shown in Figure 2. As evident from Figure 2, such interpolation, although capable of picking the nonstationary data pattern, remains one sided and may create discontinuities on the other side of the data gap. To help with this issue, we can rerun the streaming PEF using the opposite direction and stack the results.

Figure 3 shows a 2D missing data interpolation test using a synthetic pattern from Claerbout and Brown (1999). The input data for this test are shown in Figure 4b. The interpolation proceeds by running the streaming PEF twice in forward and backward directions (Figure 3a and 3b) and averaging the interpolation results (Figure 3c). The cost of such a procedure is the cost of two convolutions, as opposed to multiple iterations required in the conventional approach to missing data interpolation with PEFs (Naghizadeh and Sacchi, 2010; Liu and Fomel, 2011; Claerbout, 2014). A better interpolation result, shown in Figure 4c, is achieved by filling the residual inside the gap with small-variance white noise instead of zeros while reconstructing the data from the residual.

The random residual approach can generate multiple equiprobable realizations for the missing data interpolation problem in the spirit of geostatistical stochastic simulations (Clapp, 2000). Figure 5 shows three realizations created using different seeds for the pseudorandom number generator.

Figure 6 shows a similar interpolation test applied to 2D data with a nonstationary pattern. Although streaming PEF fails to achieve a perfect reconstruction in this case, it manages to capture most of the pattern and hide the location of the gap.

Figure 7 shows a similar test applied to 2D seismic data. The pattern made by reflection events with variable slopes is well captured.

Removing the dominant reflection events by PEF filtering leaves weaker hyperbolic diffraction events (Figure 8).

DISCUSSION

A streaming application, wherein the adaptive filter is updated one data point at a time, is appropriate in a continuous stream of a large amount of data, such as in passive monitoring of carbon storage (Alumbaugh et al., 2024). In more typical scenarios, wherein the data are immediately accessible, more accurate results can be achieved with other forms of regularization, such as regularizing adaptive filter coefficients with smoothing shaping operators (Fomel, 2009; Liu and Fomel, 2011; Liu et al., 2012). For greater efficiency, a hybrid regularization strategy can be developed. Such strategies are discussed by Geng et al. (2024), who also extend the streaming approach to other applications of seismic attributes.

In applications such as missing data reconstruction, the streaming approach represents an extreme point of the accuracy-efficiency trade-off. To achieve better accuracy, some efficiency can be sacrificed at the expense of performing extra iterations. The presented approach is also limited to the situations of missing values in regularly sampled data and will need to be modified for situations of irregular data.

In multidimensional applications, the helical boundary conditions allow for easy invertibility (Claerbout, 1998) but are not always appropriate. In this case, extra accuracy can also be bought at the cost of sacrificing some of the efficiency.

We hope that providing reproducible benchmarks with this paper will invite the reader to make direct comparisons with more advanced methods. Such comparisons are beyond the scope of this paper because they would require a different codebase.

CONCLUSION

We have presented an efficient approach for computing and applying nonstationary PEFs. Instead of storing multiple copies of varying filters, the streaming approach stores only one copy and updates the filter on the fly with every new data point. The cost of this procedure is equivalent to the cost of a single convolution and does not require multiple iterations. Moreover, the nonlinear operation of estimating and applying a streaming PEF has an exact inverse, which becomes helpful in missing data interpolation problems. A streaming approach for missing data interpolation does not require iterations and can be accomplished effectively at the cost of two convolutions. We anticipate many possible applications of the proposed technique in geophysical estimation problems.

All examples in this paper are reproducible using the attached Jupyter Notebook.

ACKNOWLEDGMENTS

We thank R. Clapp, S. Levin, Y. Liu, and R. Sarkar for inspiring discussions about prediction-error filters. We thank the sponsors of the Texas Consortium for Computational Seismology (TCCS) and the Stanford Exploration Project (SEP) for their financial support.

DATA AND MATERIALS AVAILABILITY

The data associated with this research are available and can be accessed via the following URL: https://doi.org/10.5281/zenodo.11099632.

Biographies and photographs of the authors are not available.

Figure 1.

View large Download slide

From top to bottom: the input nonstationary data, streaming deconvolution residual, and data reconstructed from the residual by the inverse streaming operation.

Figure 2.

View large Download slide

The 1D missing data interpolation example. From top to bottom: the ideal data, input data with missing gaps, and data reconstructed by streaming PEF.

Figure 3.

View large Download slide

The 2D missing data interpolation example. (a) Data reconstructed using forward streaming, (b) data reconstructed using backward streaming, and (c) sum of the forward and backward reconstructions. The residual is set to zero inside the gap.

Figure 4.

View large Download slide

The 2D missing data interpolation example. (a) The ideal data, (b) input data with a missing gap, and (c) data reconstructed by the streaming PEF. The residual is set to small-variance white noise inside the gap.

Figure 5.

View large Download slide

Three different realizations of the missing data interpolation. The images coincide at the locations of the known data but differ inside the gap.

Figure 6.

View large Download slide

(a–c) The 2D missing data interpolation example with a nonstationary data pattern. (a) The ideal data, (b) input data with a missing gap, and (c) data reconstructed by the streaming PEF.

Figure 7.

View large Download slide

(a–c) The 2D missing data interpolation with seismic data. (a) The ideal data, (b) input data with a missing gap, and (c) data reconstructed by the streaming PEF.

Figure 8.

View large Download slide

(a) The seismic data, (b) the residual after filtering streaming PEF, and (c) the data reconstruction from the residual by the inverse streaming operation. The reconstruction of the data from the residual is implemented in Program 5.

Program 1.

Julia code for streaming PEF.

function stream!(inv::Bool, d::Vector{T}, r::Vector{T}, na::Int, λ::Real) where T <: Real

# inv: inverse flag (d → r or r → d)

# d: data

# r: residual

# na: size of PEF

# λ: regularization parameter

a = zeros(T, na)#streaming PEF

dd = da = zero(T)#d (dot) d, d (dot) a

for ia in 1:na

if (inv)

d[ia] = r[ia]

else

r[ia] = d[ia]

end

dd += d[ia]*d[ia]

end

for i1 in na + 1:n1

if (inv)# from r to d

rn = r[i1] / λ

dn = rn * (λ + dd) - da

d[i1] = dn

else# from d to r

dn = d[i1]

rn = (dn + da)/(λ + dd)

r[i1] = λ * rn

end

# update PEF

for ia in 1:na

a[ia] -= rn * d[i1-ia]

end

# update dd and da

dd += dn*dn - d[i1-na] * d[i1-na]

da = dn * a[1]

for ia in 2:na

da += a[ia] * d[i1-ia + 1]

end

View large

Table 1.

Rough cost comparison among different PEF estimation methods.

Method	Storage	Cost
Stationary PEF	$O (N_{a})$	$O (N_{a}^{2} N_{d})$
Nonstationary PEF	$O (N_{a} N_{d})$	$O (N_{a}^{2} N_{d})$
Streaming PEF	$O (N_{a})$	$O (N_{a} N_{d})$

Here, $N_{a}$ is the filter size and $N_{d}$ is the data size.

View large

Program 2.

Julia code for mapping helical filter lags.

function helix(lag::Vector{Tuple{T, T}}, ci::CartesianIndices) where T <: Integer

# lag: filter lags

# ci: grid dimensions

# output: lags along the helix

# middle of the grid

mid = CartesianIndex(Tuple(last(ci)). ÷ 2)

# helix index of middle

hmid = LinearIndices(ci)[mid]

# from Cartesian shift to helix shift

return LinearIndices(ci)[map(x -> CartesianIndex(x) + mid, lag)].- hmid

end

View large

Program 3.

Julia code for multidimensional missing data interpolation using helical streaming PEF.

function stream_missing_helix!(d, k,

lag::Vector{Tuple{I, I}},

λ::Real, std=0) where I <: Integer

# d: data (modified by the program)

# k: mask for known data

# lag: PEF lags

# λ: regularization parameter

# std: standard deviation for the residual

n1, na = length(d), length(lag)

hlag = helix(lag, CartesianIndices(d))

maxlag = maximum(hlag)

T = eltype(d)

a = zeros(T, na)# streaming PEF

dd = da = zero(T)# d (dot) d, d (dot) a

for ia in 1:na

dd += d[maxlag + 1-hlag[ia]]^2

end

for i1 in maxlag + 1:n1

if (k[i1])

dn = d[i1]

rn = (dn + da)/(λ + dd)

else# assume r_n is random

rn = std * randn()/λ

dn = rn * (λ + dd) - da

d[i1] = dn

end

# update PEF

for ia in 1:na

a[ia] -= rn * d[i1-hlag[ia]]

end

# update dd and da

dd += dn * dn - d[i1-maxlag] * d[i1-maxlag]

da = dn * a[1]

for ia in 2:na

da += a[ia] * d[i1 + 1-hlag[ia]]

end

View large

Program 4.

Julia code for missing data interpolation using streaming PEF.

function stream_missing!(d::Vector{T}, k::Vector{Bool}, na::Int, λ::Real) where T <: Real

# d: data (modified by the program)

# k: mask for known data

# na: size of PEF

# λ: regularization parameter

a = zeros(T, na)# streaming PEF

dd = da = zero(T)# d (dot) d, d (dot) a

for ia in 1:na

dd += d[ia]*d[ia]

end

for i1 in na + 1:n1

if (k[i1])# from d to r

dn = d[i1]

rn = (dn + da)/(λ + dd)

else# assume r = 0

dn = - da

rn = zero(T)

d[i1] = dn

end

# update PEF

for ia in 1:na

a[ia] -= rn * d[i1-ia]

end

# update dd and da

dd += dn*dn - d[i1-na] * d[i1-na]

da = dn * a[1]

for ia in 2:na

da += a[ia] * d[i1-ia + 1]

end

View large

Program 5.

Julia code for multidimensional helical streaming PEF.

function stream_helix!(inv::Bool, d, r, lag::Vector{Tuple{I, I}}, λ::Real) where I <: Integer

# inv: inverse flag (d → r or r → d)

# d: data

# r: residual

# lag: PEF lags

# λ: regularization parameter

n1, na = length(d), length(lag)

hlag = helix(lag, CartesianIndices(d))

maxlag = maximum(hlag)

T = eltype(d)

a = zeros(T, na)# streaming PEF

for i1 in 1:maxlag

if (inv)

d[i1] = r[i1]

else

r[i1] = d[i1]

end

dd = da = zero(T)# d (dot) d, d (dot) a

for ia in 1:na

dd += d[maxlag + 1-hlag[ia]]^2

end

for i1 in maxlag + 1:n1

if (inv)

rn = r[i1]/λ

dn = rn * (λ + dd) - da

d[i1] = dn

else

dn = d[i1]

rn = (dn + da)/(λ + dd)

r[i1] = λ * rn

end

# update PEF

for ia in 1:na

a[ia] -= rn * d[i1-hlag[ia]]

end

# update dd and da

dd += dn * dn - d[i1-maxlag] * d[i1-maxlag]

da = dn * a[1]

for ia in 2:na

da += a[ia] * d[i1 + 1-hlag[ia]]

end

View large

Streaming prediction-error filters

ABSTRACT

INTRODUCTION

Theory

Inversion

Multiple dimensions

Cost comparison

Missing data interpolation

DISCUSSION

CONCLUSION

ACKNOWLEDGMENTS

DATA AND MATERIALS AVAILABILITY

Supplementary data

Data & Figures

About GSW

Resources

Explore

Connect

Streaming prediction-error filters

ABSTRACT

INTRODUCTION

Theory

Inversion

Multiple dimensions

Cost comparison

Missing data interpolation

DISCUSSION

CONCLUSION

ACKNOWLEDGMENTS

DATA AND MATERIALS AVAILABILITY

Supplementary data

Data & Figures

About GSW

Resources

Explore

Connect

This Feature Is Available To Subscribers Only