Large-scale stochastic linear inversion using hierarchical matrices

Geophysical Inverse Theory and Regularization Problems

A package of Fortran routines using linear programming techniques in linear inverse problems as gravimetry or geomagnetism is presented. These routines, with a complete description and examples of their applications, are available at the SEG Business Office. This paper presents a brief description of the package and the general features of the method; taking into account the nonuniqueness of the solution, the bounds of some parameters or some particular solutions, such as the “ideal body” [Parker (1975)], are computed. This package can be applied to inverse gravity or geomagnetic problems with a maximum data size on IBM 360/65 of the order such that the number of measurements ⩽ 250, number of unknowns ⩽ 600.

Trace inversion for reservoir parameters is affected by angle averaging of seismic data and wavelet distortion on the migration image. In an alternative approach to stochastic trace inversion, the data are inverted prestack before migration using 3D dynamic ray tracing. This choice makes it possible to interweave trace inversion with Kirchhoff migration. The new method, called ray-based stochastic inversion, is a generalization of current amplitude versus offset/amplitude versus angle (AVO/AVA) inversion techniques. The new method outperforms standard stochastic inversion techniques in cases of reservoir parameter estimation in a structurally complex subsurface with substantial lateral velocity variations and significant reflector dips. A simplification of the method inverts the normal-incidence response from reservoirs with approximately planar layering at the subsurface target locations selected for inversion. It operates along raypaths perpendicular to the reflectors, the direction that offers optimal resolution to discern layering in a reservoir. In a test on field data from the Gulf of Mexico, reservoir parameter estimates obtained with the simplified method, the estimates found by conventional stochastic inversion, and the actual values at a well drilled after the inversion are compared. Although the new method uses only 2% of the prestack data, the result indicates it improves accuracy on the dipping part of the reservoir, where conventional stochastic inversion suffers from wavelet stretch caused by migration.

The random microstructure of heterogeneous materials causes spatial variations in mechanical properties, significantly influencing material behaviour, which cannot be fully captured by deterministic models with fixed input parameters. Contemporary research has increasingly shifted toward stochastic modelling to address two distinct sources of non-determinism: epistemic uncertainty (reducible due to limited knowledge) and aleatory uncertainty (irreducible, reflecting inherent randomness). Stochastic models require a specialised approach to estimate their underlying deterministic material parameters, whereas deterministic models are calibrated by treating their parameters as random variables. Within a probabilistic framework, the identification of random parameters is formulated as a stochastic inverse problem. Its objective is to reproduce the variability observed in experimental data from an ensemble of specimens using a numerical model with random material parameters. To mitigate the computational cost of the numerous model simulations required, two widely used surrogate models polynomial chaos expansions and artificial neural networks are investigated. The preferred polynomial chaos provides analytical statistics of model responses enabling efficient variable transformations and sensitivity analysis, which are critical for uncertainty quantification. This thesis explores advanced methods for identifying aleatory uncertainty, with a particular focus on their application to the calibration of heterogeneous material models. Traditional methods, often based on Bayesian hierarchical modelling, typically assume a prescribed form for the parameter probability density function. The developed stochastic inversion framework derives the parameter distribution directly from the observed data distribution, constrained only by the feasible domain of material parameters to ensure model stability and validity. Both approaches are compared through two case studies involving heterogeneous material models. Additionally, this thesis addresses the calibration of the Lattice Discrete Particle Model for concrete, which necessitates specialised treatment due to its inherent stochasticity arising from the random generation of the material microstructure.

The stochastic geostatistical inversion approach is widely used in subsurface inverse problems to estimate unknown parameter fields and corresponding uncertainty from noisy observations. However, the approach requires a large number of forward model runs to determine the Jacobian or sensitivity matrix, thus the computational and storage costs become prohibitive when the number of unknowns, m, and the number of observations, n increase. To overcome this challenge in large‐scale geostatistical inversion, the Principal Component Geostatistical Approach (PCGA) has recently been developed as a “matrix‐free” geostatistical inversion strategy that avoids the direct evaluation of the Jacobian matrix through the principal components (low‐rank approximation) of the prior covariance and the drift matrix with a finite difference approximation. As a result, the proposed method requires about K runs of the forward problem in each iteration independently of m and n, where K is the number of principal components and can be much less than m and n for large‐scale inverse problems. Furthermore, the PCGA is easily adaptable to different forward simulation models and various data types for which the adjoint‐state method may not be implemented suitably. In this paper, we apply the PCGA to representative subsurface inverse problems to illustrate its efficiency and scalability. The low‐rank approximation of the large‐dimensional dense prior covariance matrix is computed through a randomized eigen decomposition. A hydraulic tomography problem in which the number of observations is typically large is investigated first to validate the accuracy of the PCGA compared with the conventional geostatistical approach. Then the method is applied to a large‐scale hydraulic tomography with 3 million unknowns and it is shown that underlying subsurface structures are characterized successfully through an inversion that involves an affordable number of forward simulation runs. Lastly, we present a joint inversion of head and tracer test data using MODFLOW and MT3DMS as coupled black‐box forward simulation solvers. These applications demonstrate the advantages of the PCGA, i.e., the scalability to high‐dimensional inverse problems and the ability to utilize multiple forward models as black boxes.

SummaryUncertainty quantification for linear inverse problems remains a challenging task, especially for problems with a very large number of unknown parameters (e.g., dynamic inverse problems) and for problems where computation of the square root and inverse of the prior covariance matrix are not feasible. This work exploits Krylov subspace methods to develop and analyze new techniques for large‐scale uncertainty quantification in inverse problems. In this work, we assume that generalized Golub‐Kahan‐based methods have been used to compute an estimate of the solution, and we describe efficient methods to explore the posterior distribution. In particular, we use the generalized Golub‐Kahan bidiagonalization to derive an approximation of the posterior covariance matrix, and we provide theoretical results that quantify the accuracy of the approximate posterior covariance matrix and of the resulting posterior distribution. Then, we describe efficient methods that use the approximation to compute measures of uncertainty, including the Kullback‐Liebler divergence. We present two methods that use the preconditioned Lanczos algorithm to efficiently generate samples from the posterior distribution. Numerical examples from dynamic photoacoustic tomography demonstrate the effectiveness of the described approaches.

A new stochastic inversion algorithm to invert three-dimensional (3D) electrical resistivity tomography (ERT) data has been implemented. The statistical information about the correlation among the model parameters and the correlation among observed data was employed to stabilize the ill-posed 3D ERT inverse problem. The data noise covariance matrix was diagonal based on the assumption of uncorrelated data noises. The model parameter covariance matrix, however, could be a full matrix depending on the correlation length between the model parameters. The resulting unsymmetric linear system was solved by bi-conjugate gradient (BICG) methods, such as BICGSTAB, a stabilized variant of BICG algorithm, and TFQMR, a transpose-free quasiminimal residual algorithm. An advantage of this algorithm is that any prior knowledge about the model parameters can be easily incorporated in the inversion process in the form of covariance model type and correlation length. This is the basis for our future implementation of a joint hydrologicalgeophysical inversion. This stochastic inverse technique will be used for characterizing, monitoring, and predicting fluid movement in heterogeneous vadose zone.

Understanding the subsurface is of prime importance for many geological and hydrogeological applications. Geophysical methods offer an economical alternative for investigating the subsurface compared to costly borehole investigations. However, geophysical results are commonly obtained through deterministic inversion of data whose solution is non-unique. Alternatively, stochastic inversions investigate the full uncertainty range of the obtained models, yet are computationally more expensive. In this research, we investigate the robustness of the recently introduced Bayesian evidential learning in one dimension (BEL1D) for the stochastic inversion of time-domain electromagnetic data (TDEM). First, we analyse the impact of the accuracy of the numerical forward solver on the posterior distribution, and derive a compromise between accuracy and computational time. We also introduce a threshold-rejection method based on the data misfit after the first iteration, circumventing the need for further BEL1D iterations. Moreover, we analyse the impact of the prior-model space on the results. We apply the new BEL1D with a threshold approach on field data collected in the Luy River catchment (Vietnam) to delineate saltwater intrusions. Our results show that the proper selection of time and space discretization is essential for limiting the computational cost while maintaining the accuracy of the posterior estimation. The selection of the prior distribution has a direct impact on fitting the observed data and is crucial for a realistic uncertainty quantification. The application of BEL1D for stochastic TDEM inversion is an efficient approach, as it allows us to estimate the uncertainty at a limited cost.

The unique solvability of linear inverse coefficient problems for the evolutionary equation in a Banach space with the Caputo Fabrizio derivative is studied. An operator at the unknown function in the equation is assumed to be bounded, the equation is endowed with the Cauchy condition. For the inverse problem with a constant unknown coefficient and with an integral overdefinition condition in the sense of Riemann Stieltjes, which includes the condition of final overdefinition as a special case, a well-posedness criterion is obtained. Sufficient conditions for unique solvability and an estimate of the well-posedness for the solution are obtained for a linear inverse problem with a time-dependent unknown coefficient. The abstract results obtained are used in the study of inverse problems with an unknown coefficient depending only on spatial variables or only on time, for equations with polynomials of a self-adjoint elliptic differential operator with respect to spatial variables.

The geomagnetic components from eight low-latitude ground stations and the ATS 1 satellite are studied for substorms occurring on December 24 and 25, 1967, and on August 15, 1968. A simple, four-parameter, field-aligned line current model is employed, and the linear inverse problem is constructed to determine the set of model parameters that best fit the data. A straightforward review of linear inversion theory is included in the theory section. The four model parameters are the magnitude of the current, the L shell on which the current flows, and the azimuthal positions of the eastern and western portions of the circuit. Both growth phases and expansion phases are treated. The best fits to the ground and satellite data are shown to be independent of the L parameter for values of L ≳ 3.5 when only ground data are treated and for values of L ≳ 8.5 when satellite data are included. The four-parameter model, although it produces good qualitative fits to expansion phase data, is unable to fit the data within 2 standard deviations of the error in the observations. The model is able to fit the growth phase of one substorm within 1 standard deviation. In one substorm with an apparent growth phase a better fit of the data during the expansion phase was obtained by combining the line current model of the growth phase with an additional four-parameter line current model.

Over the last project six months, our project activities have concentrated on three areas: (1) performing a stochastic inversion of pattern 16 seismic data to deduce reservoir permeability, (2) development of the geochemical inversion strategy and implementation of associated software, and (3) completing the software implementation of TProGS and the geostatistical analysis that provides the information needed when using the software to produce realizations of the Midale reservoir. The report partially the following deliverables: D2: Model development: MCMC tool (synthetic fluid chemistry data); deliverable completed. D4: Model development/verification: MCMC tool (TProGS, field seismic/chemistry data) work product; deliverable requirements partially fulfilled. D5: Field-based single-pattern simulations work product; deliverable requirements partially fulfilled. When completed, our completed stochastic inversion tool will explicitly integrate reactive transport modeling, facies-based geostatistical methods, and a novel stochastic inversion technique to optimize agreement between observed and predicted storage performance. Such optimization will be accomplished through stepwise refinement of: (1) the reservoir model - principally its permeability magnitude and heterogeneity - and (2) geochemical parameters - primarily key mineral volume fractions and kinetic data. We anticipate that these refinements will facilitate significantly improved history matching and forward modeling of CO{sub 2} storage. Our tool uses the Markov Chain Monte Carlo (MCMC) methodology. Deliverable D1, previously submitted as a report titled ''Development of a Stochastic Inversion Tool To Optimize Agreement Between The Observed And Predicted Seismic Response To CO{sub 2} Injection/Migration in the Weyburn-Midale Project'' (Ramirez et al., 2009), described the stochastic inversion approach that will identify reservoir models that optimize agreement between the observed and predicted seismic response. The software that implements this approach has been completed, tested, and used to process seismic data from pattern 16. A previously submitted report titled ''Model verification: synthetic single pattern simulations using seismic reflection data'', Ramirez et al. 2010, partially fulfilled deliverable D3 by summarizing verification activities that evaluate the performance of the seismic software and its ability to recover reservoir model permeabilities using synthetic seismic reflection data. A future progress report will similarly describe summarizing verification activities of the geochemical inversion software, thereby completing deliverable D3. This document includes a chapter that shows and discusses permeability models produced by seismic inversion that used seismic data from pattern 16 in Phase 1A. It partially fulfills deliverable D5: Field-based single-pattern simulations work product. The D5 work product is supposed to summarize the results of applying NUFT/MCMC to refine the reservoir model and geochemical parameters by optimizing observation/prediction agreement for the seismic/geochemical response to CO{sub 2} injection/migration within a single pattern of Phase 1A/1B. A future progress report will show inversion results for the same pattern using geochemical data, thereby completing deliverable D5. This document also contains a chapter that fulfills deliverable D2: Model development: MCMC tool (synthetic fluid chemistry data). The chapter will summarize model development activities required to facilitate application of NUFT/MCMC to optimize agreement between the observed and predicted geochemical response to CO{sub 2} injection/migration. Lastly, this document also contains a chapter that partially fulfills deliverable D4: Model development/verification: MCMC tool (TProGS, field seismic/chemistry data) work product. This work product is supposed to summarize model development activities required for (1) application of TProGS to Weyburn, (2) use of TProGS within the MCMC tool, and (3) application of the MCMC tool to address field seismic and geochemical data. The chapter included here fulfills requirements 1 and 2. Requirement 3 will be addressed in a future progress report.

Particle Image Velocimetry (PIV) is a non-invasive measurement technique which resolves the flow velocity by taking instantaneous snapshots of tracer particle motion in the flow and uses digital image cross-correlation to estimate the particle shift up to subpixel accuracy. The measurement chain incorporates numerous sets of parameters, such as the particle displacements, the particle image size, the flow shear rate, the out-of-plane motion for planar PIV and image noise to name a few, and these parameters are interrelated and influence the final velocity estimate in a complicated way. In the last few decades, PIV has become widely popular by virtue of developments in both the hardware capabilities and correlation algorithms, especially with the scope of 3-component (3C) and 3-dimensional (3D) velocity measurements using stereo-PIV and tomographic-PIV techniques, respectively. The velocity field measurement not only leads to other quantities of interest such as Pressure, Reynold stresses, vorticity or even diffusion coefficient, but also provides a reference field for validating numerical simulations of complex flows. However, such a comparison with CFD or applicability of the measurement to industrial design requires one to quantify the uncertainty in the PIV estimated velocity field. Even though the PIV community had a strong impetus in minimizing the measurement error over the years, the problem of uncertainty estimation in local instantaneous PIV velocity vectors have been rather unnoticed. A typical norm had been to assign an uncertainty of 0.1 pixels for the whole field irrespective of local flow features and any variation in measurement noise. The first article on this subject was published in 2012 and since then there has been a concentrated effort to address this gap. The current dissertation is motivated by such a requirement and aims to compare the existing 2D PIV uncertainty methods, propose a new method to directly estimate the planar PIV uncertainty from the correlation plane and subsequently propose the first comprehensive methods to quantify the measurement uncertainty in stereo-PIV and 3D Particle Tracking Velocimetry (PTV) measurements.The uncertainty quantification in a PIV measurement is, however, non-trivial due to the presence of multitude of error sources and their non-linear coupling through the measurement chain transfer function. In addition, the advanced algorithms apply iterative correction process to minimize the residual which increases the complexity of the process and hence, a simple data-reduction equation for uncertainty propagation does not exist. Furthermore, the calibration or a reconstruction process in a stereo or volumetric measurement makes the uncertainty estimation more challenging. Thus, current uncertainty quantification methods develop a-posterior models utilizing the evaluated displacement information and combine it with either image information, correlation plane information or even calibration “disparity map” information to find the desired uncertainties in the velocity estimates.

Objectives: Focus of the Research experiment is to evaluate the uncertainty in Ethephon determination. Validation of analytical method and identify the factors effecting the uncertainty in measurement. Methods/Statistical Analysis: A Titrimetric Method was used for determination of Ethephon content. A certified commercially formulated solution of Ethephon 39.00 % (w/w) S.L has been taken as a sample for analysis. An estimation of uncertainty in quantification of Ethephon was determined by identifying the following parameters, uncertainty influencing factors and its contribution to measurement, Type A Uncertainty (UA), Type B Uncertainty (UB), Combined Uncertainty (UC), Expanded Uncertainty (UE) and uncertainty budget. An Experimentally calculated value of combined uncertainty in preparation of 0.1 N NaOH solutions is 0.132%. A value of relative standard uncertainty in repeatability test of Ethephon is 0.084%. Determined % bias value of test method is 0.06%. An estimated value of combined uncertainty in Ethephon determination is 0.373 %. Obtained value of expanded uncertainty with coverage factor of 2.26 with 95% confidence is 0.843%. Obtained % RSD of Ethephon test result is 0.26 %. LOQ value of analytical test method is 0.36%. A Measured uncertainty value in Ethephon determination is 39.06 ± 0.33%. Findings: This paper produces Technique for analytical method validation and uncertainty measurements in quantitative analysis of Ethephon. Application/Improvements: As per regulatory for laboratory accreditation bodies, it is mandatory to evaluate the uncertainty in measurement. Presented method provides extensive details regarding uncertainty measurements.

The resolution of many large-scale inverse problems using MCMC methods\nrequires a step of drawing samples from a high dimensional Gaussian\ndistribution. While direct Gaussian sampling techniques, such as those based on\nCholesky factorization, induce an excessive numerical complexity and memory\nrequirement, sequential coordinate sampling methods present a low rate of\nconvergence. Based on the reversible jump Markov chain framework, this paper\nproposes an efficient Gaussian sampling algorithm having a reduced computation\ncost and memory usage. The main feature of the algorithm is to perform an\napproximate resolution of a linear system with a truncation level adjusted\nusing a self-tuning adaptive scheme allowing to achieve the minimal computation\ncost. The connection between this algorithm and some existing strategies is\ndiscussed and its efficiency is illustrated on a linear inverse problem of\nimage resolution enhancement.\n

Millions of digital images are captured by imaging devices on a daily basis. The way imaging devices operate follows an integral process from which the information of the original scene needs to be estimated. The estimation is done by inverting the integral process of the imaging device with the use of optimization techniques. This linear inverse problem, the inversion of the integral acquisition process, is at the heart of several image processing applications such as denoising, deblurring, inpainting, and super-resolution. We describe in detail the use of linear inverse problems in these applications. We review and compare several state-of-the-art optimization algorithms that invert this integral process. Linear inverse problems are usually very difficult to solve. Therefore, additional prior assumptions need to be introduced to successfully estimate the output signal. Several priors have been suggested in the research literature, with the Total Variation (TV) being one of the most prominent. In this thesis, we review another prior, the l0 pseudo-norm over the gradient domain. This prior allows full control over how many non-zero gradients are retained to approximate prominent structures of the image. We show the superiority of the l0 gradient prior over the TV prior in recovering genuinely piece-wise constant signals. The l0 gradient prior has shown to produce state-of-the-art results in edge-preserving image smoothing. Moreover, this general prior can be applied to several other applications, such as edge extraction, clip-art JPEG artifact removal, non-photorealistic image rendering, detail magnification, and tone mapping. We review and evaluate several state-of-the-art algorithms that solve the optimization problem based on the l0 gradient prior. Subsequently we apply the l0 gradient prior to two applications where we show superior results as compared to the current state-of-the-art. The first application is that of single-image reflection removal. Existing solutions to this problem have shown limited success because of the highly ill-posed nature of the problem. We show that the standard l0 gradient prior with a modified data-fidelity term based on the Laplacian operator is able to sufficiently remove unwanted reflections from images in many realistic scenarios. We conduct extensive experiments and show that our method outperforms the state-of-the-art. In the second application of haze removal from visible-NIR image pairs we propose a novel optimization framework, where the prior term penalizes the number of non-zero gradients of the difference between the output and the NIR image. Due to the longer wavelengths of NIR, an image taken in the NIR spectrum suffers significantly less from haze artifacts. Using this prior term, we are able to transfer details from the haze-free NIR image to the final result. We show that our formulation provides state-of-the-art results compared to haze removal methods that use a single image and also to those that are based on visible-NIR image pairs.