NG31B-0870 0800h
Nonlinear Joint Inversion to 3D Geological Models
We formulate an inverse problem to jointly incorporate all available observations from multiple geophysical surveys, geological data, and spatial statistics. The goal is to estimate the subsurface lithology distribution within a common earth model. For the $i$th particular geophysical surveying technique (e.g., gravity) write a discrete, linearized forward problem as the matrix equation ${\bf L}_i m_i = d_i$, where ${\bf L}_i$ maps a physical property $m_i$ (e.g., density distribution) into the corresponding observed data $d_i$ (e.g., $g_z$). While uniqueness theorems exist for some special cases, in practice finding $m_i$ from $d_i$ is ill-posed leading to a need for prior information to achieve a solution. Typically, priors are imposed via an assumed regularization of $m_i$. By merging and concatenating the $m_i$ as well as concatenating the $d_i$ into longer vectors $m$ and $d$, the joint forward operator now can be written as ${\bf L}$. The operator ${\bf L}$ is almost block diagonal with entries ${\bf L}_i$. However, merging operations on the $m$ vector lead to physical coupling terms appearing off diagonal for related phenomena (e.g., gravity and linearized seismic). Next, introduce a permutation operator ${\bf P}$ which re-orders rows of $m$ such that adjacent entries in the vector correspond to physical properties at the same region of space. (Call this the spatial grouping of the properties $m_v$.) From this we construct an operator ${\bf L}^v = {\bf L P}$ which relates the region grouped property vector $m_v$ to the data $d$. Next, consider a set of column vectors constructed as follows. The coordinate basis vector $e_i$ (zeros everywhere, except for a one in the $i$th row) represents lithology type $i$. Introduce a matrix ${\bf R}$ whose $i$th column represents the property values for lithology type $i$. Then the property values in a region are given by ${\bf R} v$ where $v$ is one of the basis vectors $e_i$ representing the composition of that region. We can represent the lithological parameterization of the model $m'$ as the concatenation of the $v$s corresponding to each region. By forming a block diagonal matrix ${\bf D}$ from repeated diagonal entries ${\bf R}$, one for each region, we construct the matrix relating the regional \textit{lithology} parameterization $m'$ to the region grouped \textit{property} vector $m_v$. Putting all of this together we obtain ${\bf L}' m' = d$ where ${\bf L}' = {\bf L P D}$. This is a constrained linear relationship between a lithological parameterization of the model and the data. This system restricts the degrees of freedom (DOFs) of $m'$ to that of single lithologies for each region. Geological sampling further restricts DOFs by directly constraining components of $m'$. Spatial statistics such as indicator variograms for lithology can also serve as further priors. The $m'$ are discrete vector data in each region that represent categorical variables. This makes the problem constrained, and hence nonlinear. Nevertheless, allowing linear mixing of lithologies during the solution process enables us to take advantage of linearity. If a linear mixing hypothesis is correct, such mixed lithology solutions might even be admissible. Nonlinear ${\bf L}$ could be treated in a similar fashion.
NG31B-0871 0800h
3D Inversion of magnetic total gradient in the presence of remanent magnetization
We present a general approach for inverting magnetic data in the presence of strong remanent magnetization. Inversion of magnetic data in the presence of remanence has been hampered by the need to specify the direction of magnetization. In such cases, the remanent magnetization can have a direction significantly different from that of the current inducing field and a magnitude large enough to alter the direction of total magnetization. As a result, the magnetization direction becomes an unknown quantity. Current interpretation techniques such as 3D inversion for susceptibility rely upon knowledge of this direction and their application becomes problematical. However, the total gradient (defined as the magnitude of the gradient vector of magnetic anomaly data) is independent of the magnetization direction when the causative bodies are two-dimensional and is weakly dependent upon the magnetization direction in three dimensions. Therefore, we can invert total gradient directly to recover the magnitude of magnetization without a precise knowledge of its direction. This is a nonlinear inverse problem because the total gradient depends nonlinearly upon the magnetization. It is also necessary to impose a bound constraint to ensure that the inverted magnitude of magnetization is nonnegative. This introduces further nonlinearity into the problem. We formulate the inversion by using Tikhonov regularization with the added bound constraints, and solve the resulting optimization using a primal logarithmic barrier method. We will present the algorithm and illustrate it with both synthetic and field examples. The ability to invert magnetic data with little information about the nature of remanent magnetization increases the areas in which we can apply 3D inversion of magnetic data. It opens the door to more quantitative interpretation of data in a variety of practical problems ranging from kimberlite imaging in diamond exploration, structural imaging in strongly magnetic terrain in mineral exploration, archeological investigations, and crustal studies.
NG31B-0872 0800h
Wavelet Domain Inversion Using Iteratively Reweighted Least Squares
We propose an approach for solving the discretized linear inverse problems in the wavelet domain. The solution minimizes the penalized least squares criterion with penalty being the level dependent $L_1$ norm of the Discreet Wavelet Transform of the discretized underlying signal. The solution of this mixed $L_2$ and $L_1$ norm minimization is non-linear. The non-linear solution in wavelet domain is sparse and is more appropriate for non-smooth signals. The numerical computation is carried out by iteratively reweighted least squares incorporating a line search acceleration. The proposed approach is compared with the classic regularization in extensive simulations. The tuning parameters are selected by data dependent methods including $L$-curve and cross validation.
NG31B-0873 0800h
Support Vector Machines for Non-linear Geophysical Inversion
Classical non-linear geophysical inversion can be simulated using computer learning via Support Vector Machines. Geophysical inverse problems are almost always ill-posed which means that many different models (i.e. descriptions of the earth) can be found to explain a given noisy or incomplete data set. Regularization and constraints encourage inversions to find physically realistic models. The set of preferred models needs to be defined a priori using as much geologic knowledge as is available. In inversion, it is assumed that data and a forward modeling process is known. The goal is to solve for a model. In the SVM paradigm, a series of models and associated data are known. The goal is to solve for a reverse modeling process. Starting with a series of initial models assembled using all available geologic information, synthetic data is created using the most realistic forward modeling program available. With the synthetic data as inputs and the known models as outputs, a Support Vector Machine is trained to approximate a local inverse to the forward modeling program. The advantages of this approach are that it is honest about the need to establish, a priori, the kinds of models that are reasonable in a particular field situation. There is no need to adjust the forward process to accommodate inversion, because SVMs can be easily modified to capture complicated, non-linear relationships. SVMs are transparent and require very little programming. If an SVM is trained using model/data pairs that are drawn from the same probability distribution that is implicit in the regularization of an inversion, then it will get very similar results to the inversion. Because SVMs can interpret as much data as desired so long as the conditions of an experiment do not change, they can be used to perform otherwise computationally expensive procedures. Support Vector Machines are trained to emulate non-linear seismic Amplitude Variation with Offset (AVO) inversions, gravity inversions and electromagnetic inversions. Training an SVM, including generating training data, is generally much faster than performing a non-linear inversion.
NG31B-0874 0800h
Optimal Survey Design Using Appraisal Analysis
Geophysical inversion produces images of the Earth model from finite noisy data. The scale of the model that can be resolved by inversion is affected by the sampling strategy of the data acquisition, associated data noise and a priori information that regularizes the ill-posed solution. Therefore an important questions is: How best to design surveys to obtain enhanced resolution in specific regions of the model? Thus the objective in survey design is to determine optimal survey parameters, such as position of sources/receivers and possibly frequencies in EM experiments, that would provide 'better' model resolution in a region of interest. We pose this as an inverse problem by maximising a resolution measure, the point spread function. The point spread function quantifies how an impulse in the true model is observed in the inversion result and hence the goal is to adjust the survey parameters so that the point spread function is as delta-like as possible. This is solved as a nonlinear optimisation problem with constraints on the parameters. Due to the highly nonlinear nature of the problem we examine two approaches for its solution. The first is a local, (Newton) strategy that use a primal interior point method to incorporate bounds on the parameters. We formulate it such that it can be applied to large scale problems. The second approach is a global method, simulated annealing. This proposed method is capable of enhancing resolution locally as well as globally. We examine various aspects of our method such as effects of data noise, data redundancy, multiregion survey design and compare with survey design techniques that enhances model resolution globally. We illustrate our methodology with a cross-well tomography example.
NG31B-0875 0800h
Fresnel Zone Georadar Attenuation Difference Tomography
Georadar attenuation difference tomography is useful tool for imaging temporal and spatial changes in bulk electrical conductivity due to fluid flow. The most common method of attenuation difference tomography employs the ray approximation where waves are assumed to propagate at infinite frequency. This approximation causes significant model error that generates artifacts and loss of resolution in the inverse estimates. In this paper we propose an efficient method of computing Fresnel volume sensitivities using scattering theory. These sensitivities account for finite frequency propagation and represent the physics of electromagnetic propagation more accurately than ray theory. Consequently the Fresnel theory sensitivity matrices predict data more accurately than the ray theory counterpart. We use singular value decomposition (SVD) analysis to show how and why this physical improvement allows Fresnel volume inversions to localize targets and resolve bulk conductivity changes better than ray-based inversions. By decomposing the sensitivity kernels using SVD we analyze the nature of the basis functions used to construct the inverse estimates for a synthetic case. The basis functions corresponding to the Fresnel and ray theory sensitivity kernels display a significantly different character. The Fresnel basis functions are more localized and are less oscillatory than the ray theory counterparts. The ray-based basis functions quickly become oscillatory and are dominated by the well know X-pattern. The slowly varying and localized composition of the Fresnel based basis function allows Fresnel volume inversions to resolve localized targets more accurately than ray based inversions. In addition,more higher order basis functions are required to appropriately fit the data in the ray based case. Thus we require fewer basis functions from the Fresnel operator to achieve the same misfit level. Computational efficiency of the scattering approach compared to the finite difference solution of Fresnel sensitivity will be presented in this work. Fresnel theory increases the utility of attenuation difference tomography in providing valuable flow and transport information for monitoring or estimating hydrogeologic parameters.
NG31B-0876 0800h
Statistical Inversion for Quantifying Uncertainties in Climate Prediction
The effort to estimate uncertainties in climate prediction, as similar to many problems in geophysical inversion, is defined by the dual challenges of 1) the non-linear dependencies between sources of modeling uncertainty and 2) the computational expense of any single model evaluation. Our analyses have focused on Greedy importance sampling based on multiple Very Fast Simulated Annealing (VFSA) as an ideal method for efficiently meeting the dual objectives of identifying problem solutions and sampling a multidimensional probably distribution for describing solution uncertainties. One question we have addressed is the extent to which problem characteristics affect the efficiency of multiple VFSA. We have found that certain algorithmic choices can be set based on broader characteristics of a problem such as its dimensionality and not the degree to which nonlinearities are important. We shall also present our approach to estimating observational uncertainties from a climate model, and how this knowledge may be incorporated into the analysis and expressed within estimates of the uncertainty in the parameter inversion.
http://www.ig.utexas.edu/people/staff/charles/uncertainties_in_model_predictio.htm
NG31B-0877 0800h
Phase Changes and State Estimation for Non-linear Systems
Prediction of atmospheric and oceanic states is difficult partly because of the lack of complete knowledge of the present states. An important goal of the data assimilation task is to provide such initial conditions for the atmospheric and oceanic dynamic equations used for forecasting, by updating the model results with some timely observational data. Data assimilation practiced today is mostly based on realization of linear estimation theory. In particular, a jointly-Gaussian distribution of the estimation errors and quasi-linear treatment of covariance dynamics have been assumed and practiced. This does not take account of the weather and climatic "regimes" (favored state values) that have been recognized both dynamically (e.g. by Lorenz 1963) and empirically. A fundamental problem is that a Gaussian distribution is uni-modal, while weather and climate often display multi-modal characteristics known also as "dynamic regimes". We have applied Monte Carlo numerical techniques to generic data assimilation problems without the assumptions of Gaussian distribution and quasi-linear dynamics. Specifically, we have extended the well-known assimilation technique of "ensemble Kalman filter" to a smoothing algorithm based on some theoretical developments by Kitagawa (1996) and others. The Monte-Carlo smoothing algorithm does yield more accurate state estimate than the corresponding filtering algorithm. A key component in both the filter and smoother is the generation of the random perturbations that can sample the solution space effectively and efficiently. When the size of ensemble is limited as in the case with a typical data assimilation practice, various empirically derived strategies to generate the stochastic forcing have shown potential to retain some filter/smoother performance.
NG31B-0878 0800h
On Numerical Recovery of the Cloud Moisture Structure Using Radiometric Measurements. Nonlinear Technique Solution of the Inverse Problem
Is it real to obtain a solution of a non-linear inverse problem without seeking an inverse operator but making only several simple iterations? Yes, it is possible, and the technique of [Streltsov] allows it without needs of hard work to find regularization parameters to get a stable solution. This method is also effective for remote sensing problems in which amount of measurements is less than amount of recovered parameters and additional a priory information is needed. In this work results and discussion of numerical simulation are shown of applying the technique to find solutions for a wide class of non-linear inverse problem. As an example, to formulate a problem, artificial UV radiometric data set formed by using two aircraft antennas with mutually transverse directions is constructed. For demonstration of numerical calculation results, parameters of m=40 and n=210 are chosen where m is dimension of sensor data and n is dimension of reconstructed physical profile. Reference Streltsov, Yuly. Nonlinear inverse problem. Recurrent technique. Natural and Technical Science, No.1. Moscow, 2004.
NG31B-0879 0800h
Cascadia Slow Earthquakes: Strategies for Time Independent Inversion of Displacement Fields
Continuous observations using Global Positioning System geodesy (CGPS) have revealed periodic slow or silent earthquakes along the Cascadia subduction zone with a spectrum of timing and periodicity. These creep events perturb time series of GPS observations and yield coherent displacement fields that relate to the extent and magnitude of fault displacement. In this study, time independent inversions of the surface displacement fields that accompany eight slow earthquakes characterize slip distributions along the plate interface for each event. The inversions employed in this study utilize Okada's elastic dislocation model and a non- negative least squares approach. Methodologies for optimizing the slip distribution smoothing parameter for a particular station distribution have also been investigated, significantly reducing the number of possible slip distributions and the range of estimates for total moment release for each event. The discretized slip distribution calculated for multiple creep events identifies areas of the Cascadia plate interface where slip persistently recurs. The current hypothesis, that slow earthquakes are modulated by forced fluid flow, leads to the possibility that some regions of the Cascadia plate interface may display fault patches preferentially exploited by fluid flow. Thus, the identification of regions of the plate interface that repeatedly slip during slow events may yield important information regarding the identification of these fluid pathways.
NG31B-0880 0800h
On Transverse Components of the Large-Scale Solar Magnetic Fields. Inverse Problem
Observed large-scale Solar magnetic fields data (Synoptic Data and daily magnetograms of the Wilcox Solar Observatory) were used to calculate statistical principal components (PC) of these fields. As it is known, all the PCs are mathematically orthogonal, i.e. transverse in every mathematical subspace. It is supposed that every of these components or some combination of them represents a physical property of the field, i.e. one of the PCs is responsible for the sight of view field component (or, in some cases, a radial one), another PCs bring properties of transverse components of the field (i.e. North-South and East-West). An inverse problem of reconstructing the field using the whole set of PCs leads to a linear algebra problem. Examples of such solution are presented. If one uses subsets of PCs to reconstruct the field, an ill-imposed inverse problem arises. As examples, within the inverse problem framework calculated and shown are fields using the first PC, second PC and the subset of these two components. It is proposed that such reconstructed maps could be used in boundary value problems, solar dynamo models, and to compare them with observed events on the Sun (filaments, flares, coronal mass ejections, radio emission sources etc).
NG31B-0881 0800h
Nonlinear complex principal component analysis of nearshore bathymetry
Spatially and temporally extensive nearshore bathymetric datasets have recently been analyzed with complex principal component analysis (CPCA) and a nonlinear generalisation of PCA known as NLPCA.cir to extract propagating spatial patterns that constitute most to the dominant lower-dimensional structure in the datasets. Both of these methods, however, do not fully capture the nonlinear features in the datasets. NLPCA.cir is restricted to extracting only nonlinear phase information (missing the amplitude variability), whereas CPCA captures both the amplitude and phase information but only linearly. A new Neural Network method called the nonlinear CPCA (NLCPCA) overcomes the deficiency of both the CPCA and NLPCA.cir, and captures both the phase and amplitude nonlinearly. In this study we examine the applicability of the NLCPCA, NLPCA.cir and CPCA methods to bathymetry data at Egmond (Netherlands), Hasaki Coast (Japan), and Duck (North Carolina, USA). All 3 sites are characterized by sandbars that have multiannual lifetimes and behave in an interannual quasi-periodic offshore directed manner. A cycle comprises bar birth in the inner nearshore, followed by up to several years of net offshore migration and final disappearance in the outer nearshore zone. At Duck, the underlying low-dimensional data structure was found to have only linear phase and amplitude variability and is well modelled by the CPCA. At Egmond, the data has a notable nonlinear phase variability (that is, rather peaked bars and wide shallow bar troughs) and is well modelled by the NLPCA.cir. At Hasaki, the data structure displays not only nonlinear phase variability but also amplitude variability, and the NLCPCA method is needed to model Hasaki well. In any propagating phenomena, it is difficult to know the structure of the data in advance as to which of one of the three mehthods should be used. In this study the simplest model representing well the data structure at Duck, Egmond and Hasaki is the CPCA, NLPCA.cir and NLCPCA, respectively. To avoid choosing a model, the generalized NLCPCA model can be used for all 3 cases.