Maximum-likelihood estimation of the Matérn covariance structure of isotropic spatial random fields on finite, sampled grids

TL;DR

The paper develops a statistically rigorous and computationally efficient spectral-domain maximum-likelihood estimator for isotropic Matérn Gaussian fields on finite, sampled grids. It advances the debiased Whittle likelihood by exactly accounting for discretization and edge effects through a blurred spectral density, and it provides full uncertainty quantification including the parameter covariance. The framework delivers practical tools for estimating $(\sigma^2,\nu,\rho)$, validating model adequacy with a dedicated residual test, and applying to real geophysical data with interpretable results. The authors also supply open-source Matlab and Python resources to enable replication and broader use across 1D–3D and multi-variate settings.

Abstract

We present a statistically and computationally efficient spectral-domain maximum-likelihood procedure to solve for the structure of Gaussian spatial random fields within the Matern covariance hyperclass. For univariate, stationary, and isotropic fields, the three controlling parameters are the process variance, smoothness, and range. The debiased Whittle likelihood maximization explicitly treats discretization and edge effects for finite sampled regions in parameter estimation and uncertainty quantification. As even the best parameter estimate may not be good enough, we provide a test for whether the model specification itself warrants rejection. Our results are practical and relevant for the study of a variety of geophysical fields, and for spatial interpolation, out-of-sample extension, kriging, machine learning, and feature detection of geological data. We present procedural details and high-level results on real-world examples.

Figures (16)

• Figure 1: Random fields generated from stationary isotropic Matérn models with variances ${\sigma^2}$, differentiabilities $\nu$, and correlation lengths $\rho$, as indicated. (Top:) Normalized spectral densities, ${\mathcal{S}}_{\boldsymbol{\theta}}(k)/{\mathcal{S}}_{\boldsymbol{\theta}}(0)$, from eq. (\ref{['matern2']}). The vertical black lines identify the wavenumbers $k_\alpha$ at which the power reaches $100\times\alpha$ per cent of the variance, from eq. (\ref{['kasolve']}), in wavelengths $\lambda_{100\alpha}=2\pi/k_\alpha$, as labeled. (Middle:) Correlations, the normalized spatial covariances, ${\mathcal{C}}_{\boldsymbol{\theta}}(r)/{\sigma^2}$, from eq. (\ref{['Kr']}). The vertical blue lines are drawn at the values $\pi\rho$, the distances at which the correlations die down to approximately one third of the variance. (Bottom:) Field realizations. The blue circles have radii $\pi\rho$, drawn for visual guidance. In the titles, $\mathsf{m}$ and $\mathsf{s}$ identify the sample means and standard deviations.
• Figure 2: The Fejér blurred spectral density $\bar{{\mathcal{S}}}_{\boldsymbol{\theta}}(\mathbf{k})$ approximates the expectation of the periodogram $|H(\mathbf{k})|^2$, of gridded and cosine-tapered data generated from the population density ${\mathcal{S}}_{\boldsymbol{\theta}}(\mathbf{k})$. (a) A single realization ${\mathcal{H}}(\mathbf{x})$, (b) its modified periodogram $|H(\mathbf{k})|^2$, and (c) the blurred spectrum $\bar{{\mathcal{S}}}_{\boldsymbol{\theta}}(\mathbf{k})$. (d) The cosine-tapered unit window, (e) the ratio of the average periodogram to the blurred spectral density, and (f) the average periodogram, over 100 realizations.
• Figure 3: The sample variance $\mathrm{s^2}$ systematically underestimates the true process variance ${\sigma^2}$. It is negatively biased by the presence of spatial correlation embodied by the Matérn parameters $\nu$ and $\rho$, listed in the titles. Black open circles ('mean') are averages of the sample variances $\mathrm{s^2}$ for data patches of different sizes, as observed over 40 lattice simulations, normalized by the actual variance ${\sigma^2}$. Solid blue lines ('full-covariance') predict the average behavior by incorporating the bias according to eq. (\ref{['svar2']}), evaluating eq. (\ref{['Kr']}) on all of the pairwise distances available in the grids. Essentially hidden underneath the blue ones are solid black lines ('blurred-likelihood') resulting from calculations that use eq. (\ref{['svar3']}). Solid red lines ('full-likelihood') are from eq. (\ref{['svar6']}). As detailed in the text, the quality of the various approximations is to be interpreted in terms of the Matérn correlation parameters $\nu$ and $\rho$, in relation to the sampling spacings $(\Delta x,\Delta y)$, which were kept constant at 10 km, and the field sizes $(M,N)$, which increased from left to right, as shown. The vertical black lines are drawn at the values $2\pi\rho$, a distance beyond which the bias in the sample variance estimator decreases to about a third of the true value, speaking empirically.
• Figure 4: The sample variance $\mathrm{s^2}$ is a biased, inconsistent, and inefficient estimator for the true process variance ${\sigma^2}$. The maximum-likelihood estimator is asymptotically unbiased, consistent and efficient. Conducting 40 lattice simulations on differently sized data patches, with Matérn parameters $(\sigma^2,\nu,\rho)$ as listed in the titles, the left panels show the behavior of the sample variance $\mathrm{s^2}$, and the right panels that of the maximum-likelihood variance estimator ('MLE'), both normalized by the actual variance ${\sigma^2}$. The gray bars span the 5th to 95th percentiles of the estimates at the quoted patch sizes, the black open circles are the mean estimates, and the solid blue lines their predictions from eq. (\ref{['svar2']}), as in Fig. \ref{['varbias']}. The magenta curves are the scaled spatial correlation functions, with the vertical black lines at $2\pi\rho$. The means of the MLE for field sizes smaller than $2\pi\rho$ were calculated over the 80th percentile of the estimates.
• Figure 5: Behavior of the relative error of the Whittle maximum-likelihood estimators of the Matérn parameters $(\sigma^2,\nu,\rho)$ for the two sets of true values listed in the titles, conducted on square lattices growing in size from $M=N=2$ to $M=N=128$ pixels of size $\Delta x=\Delta y=10~\mathrm{km}$, quoting $M\Delta x$ in multiples of $\pi\rho$. Gray bars cover the 5th through 95th percentiles of the estimates for each set of up to 40 simulations. Black filled circles are the means of the estimates, computed over the 80th percentile of the sets for fields whose linear dimension $M\Delta x<2\pi\rho$, but over the full set of up to 40 estimates beyond that size. With growing grid size, the estimates reveal themselves to be less biased with shrinking variance. The blue line tracks the absolute value of the relative error multiplied by $\sqrt{MN}$, offset by $-0.5$ for visual clarity.