BAYESIAN ESTIMATION IN SEISMIC INVERSION. PART I: PRINCIPLES1

This paper gives a review of Bayesian parameter estimation. The Bayesian approach is fundamental and applicable to all kinds of inverse problems. Its basic formulation is probabilistic. Information from data is combined with a priori information on model parameters. The result is called the a posteriori probability density function and it is the solution to the inverse problem. In practice an estimate of the parameters is obtained by taking its maximum. Well-known estimation procedures like least-squares inversion or l₁ norm inversion result, depending on the type of noise and a priori information given. Due to the a priori information the maximum will be unique and the estimation procedures will be stable except (in theory) for the most pathological problems which are very unlikely to occur in practice. The approach of Tarantola and Valette can be derived within classical probability theory. The Bayesian approach allows a full resolution and uncertainty analysis which is discussed in Part II of the paper.     

Adjustment of regularization in ill-posed linear inverse problems by the empirical Bayes approach

Regularization is the most popular technique to overcome the null space of model parameters in geophysical inverse problems, and is implemented by including a constraint term as well as the data‐misfit term in the objective function being minimized. The weighting of the constraint term relative to the data‐fitting term is controlled by a regularization parameter, and its adjustment to obtain the best model has received much attention. The empirical Bayes approach discussed in this paper determines the optimum value of the regularization parameter from a given data set. The regularization term can be regarded as representing a priori information about the model parameters. The empirical Bayes approach and its more practical variant, Akaike's Bayesian Information Criterion, adjust the regularization parameter automatically in response to the level of data noise and to the suitability of the assumed a priori model information for the given data. When the noise level is high, the regularization parameter is made large, which means that the a priori information is emphasized. If the assumed a priori information is not suitable for the given data, the regularization parameter is made small. Both these behaviours are desirable characteristics for the regularized solutions of practical inverse problems. Four simple examples are presented to illustrate these characteristics for an underdetermined problem, a problem adopting an improper prior constraint and a problem having an unknown data variance, all frequently encountered geophysical inverse problems. Numerical experiments using Akaike's Bayesian Information Criterion for synthetic data provide results consistent with these characteristics. In addition, concerning the selection of an appropriate type of a priori model information, a comparison between four types of difference‐operator model – the zeroth‐, first‐, second‐ and third‐order difference‐operator models – suggests that the automatic determination of the optimum regularization parameter becomes more difficult with increasing order of the difference operators. Accordingly, taking the effect of data noise into account, it is better to employ the lower‐order difference‐operator models for inversions of noisy data.     

On the geostatistical approach to the inverse problem

The geostatistical approach to the inverse problem is discussed with emphasis on the importance of structural analysis. Although the geostatistical approach is occasionally misconstrued as mere cokriging, in fact it consists of two steps: estimation of statistical parameters (“structural analysis”) followed by estimation of the distributed parameter conditional on the observations (“cokriging” or “weighted least squares”). It is argued that in inverse problems, which are algebraically undetermined, the challenge is not so much to reproduce the data as to select an algorithm with the prospect of giving good estimates where there are no observations. The essence of the geostatistical approach is that instead of adjusting a grid-dependent and potentially large number of block conductivities (or other distributed parameters), a small number of structural parameters are fitted to the data. Once this fitting is accomplished, the estimation of block conductivities ensues in a predetermined fashion without fitting of additional parameters. Also, the methodology is compared with a straightforward maximum a posteriori probability estimation method. It is shown that the fundamental differences between the two approaches are: (a) they use different principles to separate the estimation of covariance parameters from the estimation of the spatial variable; (b) the method for covariance parameter estimation in the geostatistical approach produces statistically unbiased estimates of the parameters that are not strongly dependent on the discretization, while the other method is biased and its bias becomes worse by refining the discretization into zones with different conductivity.     

Use of VFSA for resolution, sensitivity and uncertainty analysis in 1D DC resistivity and IP inversion

We present results from the resolution and sensitivity analysis of 1D DC resistivity and IP sounding data using a non-linear inversion. The inversion scheme uses a theoretically correct Metropolis–Gibbs' sampling technique and an approximate method using numerous models sampled by a global optimization algorithm called very fast simulated annealing (VFSA). VFSA has recently been found to be computationally efficient in several geophysical parameter estimation problems. Unlike conventional simulated annealing (SA), in VFSA the perturbations are generated from the model parameters according to a Cauchy-like distribution whose shape changes with each iteration. This results in an algorithm that converges much faster than a standard SA. In the course of finding the optimal solution, VFSA samples several models from the search space. All these models can be used to obtain estimates of uncertainty in the derived solution. This method makes no assumptions about the shape of an a posteriori probability density function in the model space. Here, we carry out a VFSA-based sensitivity analysis with several synthetic and field sounding data sets for resistivity and IP. The resolution capability of the VFSA algorithm as seen from the sensitivity analysis is satisfactory. The interpretation of VES and IP sounding data by VFSA, incorporating resolution, sensitivity and uncertainty of layer parameters, would generally be more useful than the conventional best-fit techniques.     

Incorporating uncertainty in hydrological predictions for gauged and ungauged basins in southern Africa

Abstract

The increasing demand for water in southern Africa necessitates adequate quantification of current freshwater resources. Watershed models are the standard tool used to generate continuous estimates of streamflow and other hydrological variables. However, the accuracy of the results is often not quantified, and model assessment is hindered by a scarcity of historical observations. Quantifying the uncertainty in hydrological estimates would increase the value and credibility of predictions. A model-independent framework aimed at achieving consistency in incorporating and analysing uncertainty within water resources estimation tools in gauged and ungauged basins is presented. Uncertainty estimation in ungauged basins is achieved via two strategies: a local approach for a priori model parameter estimation from physical catchment characteristics, and a regional approach to regionalize signatures of catchment behaviour that can be used to constrain model outputs. We compare these two sources of information in the data-scarce region of South Africa. The results show that both approaches are capable of uncertainty reduction, but that their relative values vary.

Editor D. Koutsoyiannis

Citation Kapangaziwiri, E., Hughes, D.A., and Wagener, T., 2012. Incorporating uncertainty in hydrological predictions for gauged and ungauged basins in southern Africa. Hydrological Sciences Journal, 57 (5), 1000–1019.     

Bayesian inference of the Cole–Cole parameters from time- and frequency-domain induced polarization

The inversion of induced‐polarization parameters is important in the characterization of the frequency electrical response of porous rocks. A Bayesian approach is developed to invert these parameters assuming the electrical response is described by a Cole–Cole model in the time or frequency domain. We show that the Bayesian approach provides a better analysis of the uncertainty associated with the parameters of the Cole–Cole model compared with more conventional methods based on the minimization of a cost function using the least‐squares criterion. This is due to the strong non‐linearity of the inverse problem and non‐uniqueness of the solution in the time domain. The Bayesian approach consists of propagating the information provided by the measurements through the model and combining this information with a priori knowledge of the data. Our analysis demonstrates that the uncertainty in estimating the Cole–Cole model parameters from induced‐polarization data is much higher for measurements performed in the time domain than in the frequency domain. Our conclusion is that it is very difficult, if not impossible, to retrieve the correct value of the Cole–Cole parameters from time‐domain induced‐polarization data using standard least‐squares methods. In contrast, the Cole–Cole parameters can be more correctly inverted in the frequency domain. These results are also valid for other models describing the induced‐polarization spectral response, such as the Cole–Davidson or power law models.     

An approach to estimation problems containing uncertain parameters

We consider how to treat a finite-dimensional linear inverse problem when the form of the forward problem is known exactly, but is dependent upon some parameters whose exact value is uncertain and which enter the forward problem multiplicatively. We show one way to proceed when the uncertainty is treatable in a statistical manner. Predicting the secular variation ∂_tB(t) produced by a particular fluid flow V at the core-mantle boundary (when magnetic diffusion is ignored) is one such example, because the results depend on the main magnetic field B(t) originating in the core which is improperly known because of contamination by the crustal magnetic field. This infinite-dimensional inverse problem is often solved by projection on to a finite-dimensional basis, and the resulting parameters found by a maximum likelihood technique. If the main field is contaminated with errors from a Gaussian distribution, this paper describes how the maximum likelihood solution can take this into account, and we show the probability density function that must be maximised in this case. We give an example of the effects for a simple model system, and suggest possible areas of application.     

Uncertainty and multiple objective calibration in regional water balance modelling: case study in 320 Austrian catchments

We examine the value of additional information in multiple objective calibration in terms of model performance and parameter uncertainty. We calibrate and validate a semi‐distributed conceptual catchment model for two 11‐year periods in 320 Austrian catchments and test three approaches of parameter calibration: (a) traditional single objective calibration (SINGLE) on daily runoff; (b) multiple objective calibration (MULTI) using daily runoff and snow cover data; (c) multiple objective calibration (APRIORI) that incorporates an a priori expert guess about the parameter distribution as additional information to runoff and snow cover data. Results indicate that the MULTI approach performs slightly poorer than the SINGLE approach in terms of runoff simulations, but significantly better in terms of snow cover simulations. The APRIORI approach is essentially as good as the SINGLE approach in terms of runoff simulations but is slightly poorer than the MULTI approach in terms of snow cover simulations. An analysis of the parameter uncertainty indicates that the MULTI approach significantly decreases the uncertainty of the model parameters related to snow processes but does not decrease the uncertainty of other model parameters as compared to the SINGLE case. The APRIORI approach tends to decrease the uncertainty of all model parameters as compared to the SINGLE case. Copyright © 2006 John Wiley & Sons, Ltd.     

Petrophysical inversion of AVA data

We investigate the interactions between the elastic parameters, V_P, V_S and density, estimated by non-linear inversion of AVA data, and the petrophysical parameters, depth (pressure), porosity, clay content and fluid saturation, of an actual gas-bearing reservoir. In particular, we study how the ambiguous solutions derived from the non-uniqueness of the seismic inversion affect the estimates of relevant rock properties. It results that the physically admissible values of the rock properties greatly reduce the range of possible seismic solutions and this range contains the actual values given by the well. By means of a statistical inversion, we analyse how approximate a priori knowledge of the petrophysical properties and of their relationships with the seismic parameters can be of help in reducing the ambiguity of the inversion solutions and eventually in estimating the petrophysical properties of the specific target reservoir. This statistical inversion allows the determination of the most likely values of the sought rock properties along with their uncertainty ranges. The results show that the porosity is the best-resolved rock property, with its most likely value closely approaching the actual value found by the well, even when we insert somewhat erroneous a priori information. The hydrocarbon saturation is the second best-resolved parameter, but its most likely value does not match the well data. The depth of the target interface is the least-resolved parameter and its most likely value is strongly dependent on a priori information. Although no general conclusions can be drawn from the results of this exercise, we envisage that the proposed AVA–petrophysical inversion and its possible extensions may be of use in reservoir characterization.     

Estimation of soil hydraulic properties and their uncertainty: comparison between laboratory and field experiment

Often the soil hydraulic parameters are obtained by the inversion of measured data (e.g. soil moisture, pressure head, and cumulative infiltration, etc.). However, the inverse problem in unsaturated zone is ill‐posed due to various reasons, and hence the parameters become non‐unique. The presence of multiple soil layers brings the additional complexities in the inverse modelling. The generalized likelihood uncertainty estimate (GLUE) is a useful approach to estimate the parameters and their uncertainty when dealing with soil moisture dynamics which is a highly non‐linear problem. Because the estimated parameters depend on the modelling scale, inverse modelling carried out on laboratory data and field data may provide independent estimates. The objective of this paper is to compare the parameters and their uncertainty estimated through experiments in the laboratory and in the field and to assess which of the soil hydraulic parameters are independent of the experiment. The first two layers in the field site are characterized by Loamy sand and Loamy. The mean soil moisture and pressure head at three depths are measured with an interval of half hour for a period of 1 week using the evaporation method for the laboratory experiment, whereas soil moisture at three different depths (60, 110, and 200 cm) is measured with an interval of 1 h for 2 years for the field experiment. A one‐dimensional soil moisture model on the basis of the finite difference method was used. The calibration and validation are approximately for 1 year each. The model performance was found to be good with root mean square error (RMSE) varying from 2 to 4 cm³ cm^?3. It is found from the two experiments that mean and uncertainty in the saturated soil moisture (θ_s) and shape parameter (n) of van Genuchten equations are similar for both the soil types. Copyright © 2010 John Wiley & Sons, Ltd.     

Predicting return periods of hydrological droughts using the Pearson 3 distribution: a case from rivers in the Canadian prairies

Abstract

The standardized series of monthly and weekly flow sequences, referred to as standardized hydrological index (SHI) series, from five rivers in the Canadian prairies were subjected to return period (T_r) analysis of drought length (L). The SHI series were truncated at drought probability levels q ranging from 0.5 to 0.05 with the intention of deducing drought events and corresponding drought lengths. The values of L were fitted to the Pearson 3, the gamma (2-parameter), the exponential (1-parameter), the Weibull 3 and the Weibull (2-parameter) probability density functions (pdfs). A priori assignment of one week or one month for the location parameter in the Pearson 3 pdf proved logical and also facilitated the rapid estimation of other parameters using either the method of moments or the method of maximum likelihood. The Pearson 3 turns out to be the most suitable pdf to describe and to estimate return periods of drought lengths. At the monthly and weekly time scales, it was inferred that the sample size (T, months or weeks) of SHI series could be treated equivalent to the return period of the largest recorded drought length. At the annual time scale, however, the sample size (T, years) should be modified using either the Hazen or the Gringorten plotting position formula to reflect the actual return period of the largest recorded drought length in years.

Editor D. Koutsoyiannis; Associate editor E. Gargouri     

Linearized amplitude variation with offset and azimuth and anisotropic poroelasticity

The normal-to-shear weakness ratio is commonly used as a fracture fluid indicator, but it depends not only on the fluid types but also on the fracture intensity and internal architecture. Amplitude variation with offset and azimuth is commonly used to perform the fluid identification and fracture characterization in fractured porous rocks. We demonstrate a direct inversion approach to utilize the observable azimuthal data to estimate the decoupled fluid (fluid/porosity term) and fracture (normal and shear weaknesses) parameters instead of the calculation of normal-to-shear weakness ratio to help reduce the uncertainties in fracture characterization and fluid identification of a gas-saturated porous medium permeated by a single set of parallel vertical fractures. Based on the anisotropic poroelasticity and perturbation theory, we first derive a linearized amplitude versus offset and azimuth approximation using the scattering function to decouple the fluid indicator and fracture parameters. Incorporating Bayes formula and convolution theory, we propose a feasible direct inversion approach in a Bayesian framework to obtain the direct estimations of model parameters, in which Cauchy and Gaussian distribution are used for the a priori information of model parameters and the likelihood function, respectively. We finally use the non-linear iteratively reweighted least squares to solve the maximum a posteriori solutions of model parameters. The synthetic examples containing a moderate noise demonstrate the feasibility of the proposed approach, and the real data illustrates the stabilities of estimated fluid indicator and dry fracture parameters in gas-saturated fractured porous rocks.     

A maximum likelihood approach to nonlinear inversion under constraints

Nonuniqueness in geophysical inverse problems is naturally resolved by incorporating prior information about unknown models into observed data. In practical estimation procedures, the prior information must be quantitatively expressed. We represent the prior information in the same form as observational equations, nonlinear equations with random errors in general, and treat as data. Then we may define a posterior probability density function of model parameters for given observed data and prior data, and use the maximum likelihood criterion to solve the problem. Supposing Gaussian errors both in observed data and prior data, we obtain a simple algorithm for iterative search to find the maximum likelihood estimates. We also obtain an asymptotic expression of covariance for estimation errors, which gives a good approximation to exact covariance when the estimated model is linearly close to a true model. We demonstrate that our approach is a general extension of various inverse methods dealing with Gaussian data. By way of example, we apply the new approach to a problem of inferring the final rupture state of the 1943 Tottori earthquake (M = 7.4) from coseismic geodetic data. The example shows that the use of sufficient prior information effectively suppresses both the nonuniqueness and the nonlinearity of the problem.     

GEOELECTRIC MODELS IN ENGINEERING GEOPHYSICS*

Man's engineering activities are concentrated on the uppermost part of the earth's crust which is called engineering-geologic zone. This zone is characterized by a significant spatialtemporal variation of the physical properties status of rocks, and saturating waters. This variation determines the specificity of geophysical and, particularly, geoelectrical investigations. Planning of geoelectric investigations in the engineering-geologic zone and their subsequent interpretation requires a priori) geologic-geophysical information on the main peculiarities of the engineering-geologic and hydrogeologic conditions in the region under investigation. This information serves as a basis for the creation of an initial geoelectric model of the section. Following field investigations the model is used in interpretation. Formalization of this a priori) model can be achieved by the solution of direct geoelectric problems. An additional geologic-geophysical information realized in the model of the medium allows to diminish the effect of the “principle of equivalence” by introducing flexible limitations in the section's parameters. Further geophysical observations as well as the correlations between geophysical and engineering-geologic parameters of the section permit the following step in the specification of the geolectric model and its approximation to the real medium. Next correction of this model is made upon accumulation of additional information. The solution of inverse problems with the utilization of computer programs permits specification of the model in the general iterational cycle of interpretation.     

Recent progress in estimating uncertainty in geomagnetic field modeling

Abstract

Recent work pertaining to estimating error and accuracies in geomagnetic field modeling is reviewed from a unified viewpoint and illustrated with examples. The formulation of a finite dimensional approximation to the underlying infinite dimensional problem is developed. Central to the formulation is an inner product and norm in the solution space through which a priori information can be brought to bear on the problem. Such information is crucial to estimation of the effects of higher degree fields at the Core-Mantle boundary (CMB) because the behavior of higher degree fields is masked in our measurements by the presence of the field from the Earth's crust. Contributions to the errors in predicting geophysical quantities based on the approximate model are separated into three categories: (1) the usual error from the measurement noise; (2) the error from unmodeled fields, i.e. from sources in the crust, ionosphere, etc.; and (3) the error from truncating to a finite dimensioned solution and prediction space. The combination of the first two is termed low degree error while the third is referred to as truncation error.

The error analysis problem consists of “characterizing” the difference δz = z—z, where z is some quantity depending on the magnetic field and z is the estimate of z resulting from our model. Two approaches are discussed. The method of Confidence Set Inference (CSI) seeks to find an upper bound for |z—?|. Statistical methods, i.e. Bayesian or Stochastic Estimation, seek to estimate E(δz² ), where E is the expectation value. Estimation of both the truncation error and low degree error is discussed for both approaches. Expressions are found for an upper bound for |δz| and for E(δz² ). Of particular interest is the computation of the radial field, B., at the CMB for which error estimates are made as examples of the methods. Estimated accuracies of the Gauss coefficients are given for the various methods. In general, the lowest error estimates result when the greatest amount of a priori information is available and, indeed, the estimates for truncation error are completely dependent upon the nature of the a priori information assumed. For the most conservative approach, the error in computing point values of B_r at the CMB is unbounded and one must be content with, e.g., averages over some large area. The various assumptions about a priori information are reviewed. Work is needed to extend and develop this information. In particular, information regarding the truncated fields is needed to determine if the pessimistic bounds presently available are realistic or if there is a real physical basis for lower error estimates. Characterization of crustal fields for degree greater than 50 is needed as is more rigorous characterization of the external fields.     

Pitfalls in Nonlinear Inversion

— We discuss and illustrate graphically with simple 2-D problems, four common pitfalls in geophysical nonlinear inversion. The first one establishes that the Lagrange multiplier, used to incorporate a priori information in the geophysical inverse problem, should be the largest value still compatible with a reasonable data fitting. This procedure should be used only when the interpreter is sure about the importance of the a priori information used to stabilize the inverse problem relative to the geophysical observations. Because this is rarely the case, the user should use the smallest Lagrange multiplier still producing stable solutions. The second pitfall is an attempt to automatically estimate the Lagrange multiplier by decreasing it along the iterative process used to solve the nonlinear optimization problem. Consequently, at the last iteration, the Lagrange multiplier may be so small that the problem may become ill-posed and any computed solution in this case is meaningless. The third pitfall is related to the incorporation of a priori information by a technique known as “Jumping.” This formulation, from the viewpoint of the class of Acceptable Gradient Methods, is incomplete and may lead to a premature halt in the iteration, and, consequently, to solutions far from the true one. Finally, the fourth pitfall is an inadequate convergence criterion which stops the iteration when the data misfit drops just below the noise level, irrespective of the fact that the functional to be minimized may not have attained its minimum. This means that the a priori information has not been completely incorporated, so that this stopping criterion partially neutralizes the effect of the stabilizing functional, and opens the possibility of obtaining unstable, meaningless estimates.