Inversion of Geophysical Data

Presentation on theme: "Inversion of Geophysical Data"— Presentation transcript:

1 Inversion of Geophysical Data
강원대학교 지구물리학과 전기전자탐사 연구실 Inversion of Geophysical Data EM

2 Forward modeling Physical scale modeling Numerical modeling
강원대학교 지구물리학과 전기전자탐사 연구실 Forward modeling Physical scale modeling Numerical modeling Analytic solution Finite difference FDM Finite element FEM Integral equation IE EM

3 Inversion Inversion of geophysical data Non-linear Non-unique
강원대학교 지구물리학과 전기전자탐사 연구실 Inversion Inversion of geophysical data Non-linear Non-unique Iterative method Damped least-squares method Marquardt Smoothness constraint EM

4 The Data and Model Parameters
강원대학교 지구물리학과 전기전자탐사 연구실 The Data and Model Parameters If N measurements are performed in a particular experiment, we call it data vector of length N. If M parameters effect on data vector, we call it model parameter vector of length M. EM

5 Modeling and Inversion
강원대학교 지구물리학과 전기전자탐사 연구실 Modeling and Inversion The data and model parameters are in some way related and this relationship is modeling. Generally, the forward modeling (modeling) is to get a data vector for a given model parameter vector and the inverse modeling (inversion) is to find model parameters from a given data vector. model space model parameter vector data space data vector inversion modeling EM

6 강원대학교 지구물리학과 전기전자탐사 연구실
Linear problems EM

7 강원대학교 지구물리학과 전기전자탐사 연구실
Least-squares method The error vector is discrepancy between observed and calculated data . where, F represent the model responses for the model parameter p. Least-squares method is to find the model parameter vector p which minimize the square sum of error vector. EM

8 Least-squares method In the linear problem,
강원대학교 지구물리학과 전기전자탐사 연구실 Least-squares method In the linear problem, Then, the minimum of E is solved by setting the derivatives of E to zero which is the least-squares solution to the inverse problem Gp = d. EM

9 Non-linear inverse problems
강원대학교 지구물리학과 전기전자탐사 연구실 Non-linear inverse problems In non-linear problem, it is impossible to represent the system as a linear equation. Therefore, inversion starts from the linearization of the problem and find the model parameters which minimize the error in least-squares sense. Assuming that the p1 is an initial guess of model parameter and Dp is model perturbation which is very small. Then we can write Using the Taylor’s series expansion, the forward modeling F when p=p2 Is given by EM

10 Non-linear inverse problems
강원대학교 지구물리학과 전기전자탐사 연구실 Non-linear inverse problems Now, we can get the linear equation with respect to the perturbation vector Dp. The error vector is where J is Jacobian or sensitivity matrix. Then, the square sum of the error vector is EM

11 Non-linear inverse problems
강원대학교 지구물리학과 전기전자탐사 연구실 Non-linear inverse problems The minimum of E is solved by setting the derivatives of E with respect to model perturbation vector Dp to zero This result is very similar to that of linear inverse problem except that the system matrix G and data vector d are replaced with Jacobian matrix J and error vector e, respectively. But, most of geophysical inverse problem are usually ill-posed and the solution is non-unique and unstable. EM

12 강원대학교 지구물리학과 전기전자탐사 연구실
Optimization The well-known inversion method to deal with the ill-posedness is to minimize the object function which is composed of the data misfit and the model constraint functional. An iterative inversion is to find the model perturbation vector that minimizes model object function, fm, subject to minimizing the error (data misfit) E : where EM

13 Model Perturbation Vector
강원대학교 지구물리학과 전기전자탐사 연구실 Model Perturbation Vector The minimum of the object function S is solved by setting the derivatives of S with respect to model perturbation vector Dp to zero. EM

14 강원대학교 지구물리학과 전기전자탐사 연구실
Model object function Because of non-uniqueness, inversion of the geophysical data have a lot of solutions. Especially, under the condition that M >> N, it is Known from linear algebra that this system of equations has many (perhaps infinitely many ) solutions. How do we choose one that is useful ?  We need a quantitative way to distinguish between acceptable models. The solution is to find a solution that is "largest" or "smallest". But in order to define "smallest", we need a ruler. Norms are mathematical rulers to measure "length". EM

15 강원대학교 지구물리학과 전기전자탐사 연구실
Model object function We will define fm to be the norm of the model perturbation. fm will be called the model objective function. The procedure for selecting one model will be to find the solution with damping or smoothing constraint by minimizing Where, Wm is the model constraint or the model weighting matrix. Generally, damped least-squares inversions are divided according to what kind of the model weighting is imposed. EM

16 Data Weighting Matrix The data weighting matrix is generally given by
강원대학교 지구물리학과 전기전자탐사 연구실 Data Weighting Matrix The data weighting matrix is generally given by where sj is the standard deviation of each datum. However, this kind of data weighting has been usually not used since the standard deviations are seldom measured at the stage of the data acquisition. Thus, the data weighting matrix is regarded as the identity matrix in the usual inversion of geophysical data. EM

17 Model Weighting Matrix
강원대학교 지구물리학과 전기전자탐사 연구실 Model Weighting Matrix Inversion results are strongly dependent on the model weighting. The procedure for selecting one model weighting will be to find the solution with damping or smoothing constraint by minimizing. Generally, we can express the data weighting as follows: Where is the roughness of the model and L is the Lagrangian multiplier matrix. If n=0 and L=I, Marquardt-Levenburg method. If n=1 or 2 and L=I, smoothness constrained method (Occam) If n=1 or 2 and L is determined from the model resolution matrix and its spread function analysis, ACB method EM

18 Flow chart of iterative non-linear inversion
강원대학교 지구물리학과 전기전자탐사 연구실 Flow chart of iterative non-linear inversion Initial model Forward modeling and error estimation Compute Jacobian matrix Solve system equations Update model parameters Final output yes no Converge? EM

19 Roughness matrix The model object function is given by
강원대학교 지구물리학과 전기전자탐사 연구실 Roughness matrix The model object function is given by For example, the amount of the roughness about the ith layer may be defined as EM

20 Roughness matrix for 1D model - 1st order
강원대학교 지구물리학과 전기전자탐사 연구실 Roughness matrix for 1D model - 1st order EM

21 Roughness matrix for 1D model – 2nd order
강원대학교 지구물리학과 전기전자탐사 연구실 Roughness matrix for 1D model – 2nd order EM

22 Roughness matrix for 2D model
강원대학교 지구물리학과 전기전자탐사 연구실 Roughness matrix for 2D model M N EM

23 강원대학교 지구물리학과 전기전자탐사 연구실
Jumping method The model constraint can be imposed on the model parameter or the model perturbation. When the model parameter is constrained, the object function is given by For iterative least-squares inversion Then, the objective function is expressed EM

24 강원대학교 지구물리학과 전기전자탐사 연구실
Jumping method The object function is differentiated with respect to the model perturbation vector This approach is regarded as a jumping method. EM

25 강원대학교 지구물리학과 전기전자탐사 연구실
Jumping method Therefore, the model parameter, not model perturbation vector, can be directly obtained in the jumping method as follows; EM

26 강원대학교 지구물리학과 전기전자탐사 연구실
Jumping method The alternative method to derive the formula for the model parameter EM

27 강원대학교 지구물리학과 전기전자탐사 연구실
Creeping vs Jumping In the creeping method, the final solution lies close to the initial guess, because the modified Jacobian keeps the changes small at each step of the process. Thus, the resultant model is strongly influenced by the initial choice. In the jumping method, on the other hand, the final answer should be independent of the starting guess since the scheme solves the minimization of the original non-linear functional (Constable, et al, 1987). Generally, the jumping method is popular in America and Europe. But the creeping in Korea. The creeping method is thought to be more stable than the jumping method. EM

28 ACB (Automatic Constraint Balancing) Inversion
강원대학교 지구물리학과 전기전자탐사 연구실 ACB (Automatic Constraint Balancing) Inversion In the inversion based on the least-squares method, the regularization parameter (Lagrangian multiplier, l ) controls the resolution and stability of inversion. But the selection of the optimum value of l is practically not easy. The ACB method is a smart scheme to provide the spatially varying l which is automatically obtained by the model resolution and its spread function analysis. This approach enhances the resolution and makes the inversion process more stable ( Yi, et al, 2003). regularization parameter, l EM

29 Model Resolution Matrix
강원대학교 지구물리학과 전기전자탐사 연구실 Model Resolution Matrix In the smoothness constraint (Occam’s) inversion, the objective function is given by The model perturbation to minimize the object function is Where J+ is called as the pseudo-inverse matrix and C is a roughness matrix. EM

30 Model Resolution Matrix
강원대학교 지구물리학과 전기전자탐사 연구실 Model Resolution Matrix Assuming there exists a true solution Dpg in the inversion problem, we have Thus, the expression for the estimated model parameter gives Now, the model resolution matrix R is defined as EM

31 Model Resolution Matrix
강원대학교 지구물리학과 전기전자탐사 연구실 Model Resolution Matrix The model resolution matrix characterizes whether the model parameter can be independently predicted or resolved. If R=I, then each model parameter is uniquely determined. If R is not identity matrix, the estimates of the model parameters are really weighted average of the true model parameters (Menke, 1984). Dp= Dpg EM

32 Backus-Gilbert Spread Function
강원대학교 지구물리학과 전기전자탐사 연구실 Backus-Gilbert Spread Function Because the resolution is best when the resolution matrix is an identity matrix, one possible measure of the resolution is the size, or spread of the off-diagonal elements. Rij i Index j Good Resolution Poor Resolution The Backus-Gilbert spread function is the weighted measure of spread that weights the (i, j) element of R according to its physical distance from the diagonal element. EM

33 Spread function The Backus-Gilbert spread function is defined as where
강원대학교 지구물리학과 전기전자탐사 연구실 Spread function The Backus-Gilbert spread function is defined as where wij : weighting factor, generally physical distance between model parameter i and j. dij : 1 if Cij =/ 0, otherwise 0 Rij : model resolution matrix Generally, the large spread is related to the poor resolution and the small spread to the high resolution. Especially, when R=I, the value of spread function is zero. EM

34 Regularization parameter matrix
강원대학교 지구물리학과 전기전자탐사 연구실 Regularization parameter matrix where l min and lmax are the lower and upper bounds of the regularization parameter, and SPmin and SPmax are the lower and upper bounds of the spread function, respectively. Finally, we can write the regularization parameter matrix as follows: EM

35 Regularization parameter matrix
강원대학교 지구물리학과 전기전자탐사 연구실 Regularization parameter matrix Finally, the object function in the ACB inversion is given by Then, the perturbation is written as In the ACB inversion, L is a spatially varying regularization parameter matrix which is automatically obtained from the model resolution matrix and its spread function. Consequently, this approach attains the enhanced resolution and the stability of inversion efficiently. EM

36 Jacobian matrix in 1D inversion
강원대학교 지구물리학과 전기전자탐사 연구실 Jacobian matrix in 1D inversion EM

37 Jacobian matrix in 1D inversion
강원대학교 지구물리학과 전기전자탐사 연구실 Jacobian matrix in 1D inversion EM

38 Jacobian matrix in 2D inversion
강원대학교 지구물리학과 전기전자탐사 연구실 Jacobian matrix in 2D inversion EM