Efficient computation on large regular grids of higher-order spatial statistics via fast Fourier transform

The complex spatial structures of natural variables are often caused by geological, physicochemical, meteorological, and biological processes that have shaped the emergence of the fields. The typical prediction of the spatial distributions of these phenomena is based on second-order geostatistical models. However, this approach has limitations, given the high complexity, non-Gaussian distributions, and nonlinear spatial connectivity models inherent in geological systems. Recently, researchers have suggested using higher-order spatial statistics, based on two- and three-point spatial statistics, to better capture spatial phenomena. Nevertheless, applying these methods requires intense numerical calculations, particularly in the case of extensive geostatistical models, and becomes especially intricate when utilized for conditioning realizations, such as in inverse problems. Spatial asymmetries and higher-order spatial cumulants, as well as their generalizations, are important higher-order statistics for characterizing non-Gaussian features. In this study, we focus on third-order statistics derived from two- and three-point spatial statistics. A MATLAB program has been developed to compute efficiently these spatial statistics using the FFT algorithm. The overall approach of these programs draws inspiration from the method successfully used for the fast calculation of variograms and cross-covariances using FFT. We recall the methodology associated with the computation of direct- and cross-variograms using FFT, as well as transiograms for categorical data. Codes are created to process regular grid data, whether it is complete or incomplete. Post-processing tools have been added to help geomodelers visualize the results. Using the FFT method is faster and delivers the same results as conventional spatial methods for this type of data. These programs are particularly valuable tools for geostatistical modeling and estimation when higher-order statistics are present in the spatial structures of natural variables, providing an efficient solution to the computational challenges associated with such applications.