  In this paper, we present a novel approach to geostatistical filtering which
tackles two challenges encountered when applying this method to complex spatial
datasets: modeling the non-stationarity of the data while still being able to
work with large datasets. The approach is based on a finite element
approximation of Gaussian random fields expressed as an expansion of the
eigenfunctions of a Laplace--Beltrami operator defined to account for local
anisotropies. The numerical approximation of the resulting random fields using
a finite element approach is then leveraged to solve the scalability issue
through a matrix-free approach. Finally, two cases of application of this
approach, on simulated and real seismic data are presented.
