Application of Hierarchical Matrices to Linear Inverse Problems in Geostatistics

Application des matrices hiérarchiques aux problèmes d’inversion linéaire en géostatistique

¹ Institute for Computational and Mathematical Engineering, Stanford University - USA
² Department of Civil and Environmental Engineering, Stanford University - USA
³ Department of Mechanical Engineering, Stanford University - USA
e-mail: arvindks@stanford.edu - sivaambi@stanford.edu - yuel@stanford.edu - peterk@stanford.edu - darve@stanford.edu

^* Corresponding author

Abstract

Characterizing the uncertainty in the subsurface is an important step for exploration and extraction of natural resources, the storage of nuclear material and gasses such as natural gas or CO₂. Imaging the subsurface can be posed as an inverse problem and can be solved using the geostatistical approach [Kitanidis P.K. (2007) Geophys. Monogr. Ser. 171, 19-30, doi:10.1029/171GM04; Kitanidis (2011) doi: 10.1002/9780470685853. ch4, pp. 71-85] which is one of the many prevalent approaches. We briefly describe the geostatistical approach in the context of linear inverse problems and discuss some of the challenges in the large-scale implementation of this approach. Using the hierarchical matrix approach, we show how to reduce matrix vector products involving the dense covariance matrix from 𝒪(m²) to 𝒪(m log m), where m is the number of unknowns. Combined with a matrix-free Krylov subspace solver, this results in a much faster algorithm for solving the system of equations that arise from the geostatistical approach. We illustrate the performance of our algorithm on an application, for monitoring CO₂ concentrations using crosswell seismic tomography.

Résumé

La caractérisation de l’incertitude en subsurface est une étape importante, tant, par exemple, dans le cadre de la recherche et de l’extraction de ressources naturelles, que dans celui du stockage de matériaux radioactifs ou de gaz tels que le gaz naturel ou le CO₂. L’imagerie en subsurface peut être posée comme un problème inverse et être résolue à l’aide d’une approche géostatistique [Kitanidis P.K. (2007) Geophys. Monogr. Ser. 171, 19-30, doi :10.1029/171GM04; Kitanidis (2011) doi : 10.1002/9780470685853. ch4, pp. 71-85] qui est l’une des méthodes importantes de ce domaine. Nous commençons par décrire brièvement l’approche géostatistique dans le contexte d’un problème inverse linéaire, puis discutons les difficultés rencontrées dans le cadre d’une implémentation à grande échelle. Ensuite, en utilisant une approche basée sur les matrices hiérarchiques, nous montrons comment réduire le coût de calcul d’un produit matrice-vecteur de 𝒪(m²) à 𝒪(m log m) dans le cas de matrices de covariance denses; m désignant ici le nombre d’inconnues. Combinée avec un solveur de Krylov, qui ne requiert pas la construction explicite de la matrice, cette méthode conduit à un algorithme beaucoup plus rapide pour résoudre le système d’équations associé à l’approche géostatistique. Nous illustrons enfin la performance de notre algorithme dans un situation spécifique, surveiller la concentration de CO₂ à l’aide de la tomographie sismique entre puits.