# Contents of /trunk/doc/inversion/gravity.tex

Revision 4045 - (show annotations)
Tue Oct 30 03:57:22 2012 UTC (6 years, 11 months ago) by caltinay
File MIME type: application/x-tex
File size: 10299 byte(s)
Imported inversion cookbook.


 1 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % Copyright (c) 2003-2012 by University of Queensland 4 5 % 6 % Primary Business: Queensland, Australia 7 % Licensed under the Open Software License version 3.0 8 9 % 10 % Development until 2012 by Earth Systems Science Computational Center (ESSCC) 11 % Development since 2012 by School of Earth Sciences 12 % 13 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 14 15 \chapter{Gravity Inversion} 16 We want to recover the density field $\rho$ from an acceleration field $g_i$ known on certain locations within the domain $\Omega$. It is 17 assumed that 18 \begin{equation} 19 \Omega = [x^{min}_0, x^{max}_0] \times 20 [x^{min}_1, x^{max}_1] \times 21 [x^{min}_2, x^{max}_2] 22 \end{equation} 23 If the density field $\rho$ is known the gravitational potential $\psi$ is given 24 as the solution of the PDE 25 \begin{equation}\label{GRAV:EQU:100} 26 \psi_{,ii} = 4\pi \cdot G \cdot \rho 27 \end{equation} 28 where $G=6.67300 \cdot 10^{-11} \frac{m^3}{kg \cdot s^2}$ is the gravitational constant. 29 The gravity force $g_i$ is given 30 as the negative of the gradient of the gravity potential $\psi$: 31 \begin{equation}\label{GRAV:EQU:101} 32 g_i = - \psi_{,i} 33 \end{equation} 34 The boundary conditions for PDE~\ref{GRAV:EQU:100} are set as follows: 35 On vertical faces of the domain we assume that the horizontal gravity forces are equal to zero: 36 \begin{equation}\label{GRAV:EQU:101a} 37 \begin{array}{ll} 38 g_0=0 & \mbox{ for } x_0=x^{min}_0 \mbox{ or } x_0=x^{max}_0 \\ 39 g_1=0 & \mbox{ for } x_1=x^{min}_1 \mbox{ or } x_1=x^{max}_1 \\ 40 \end{array} 41 \end{equation} 42 which translates to 43 \begin{equation}\label{GRAV:EQU:101aa} 44 n_i \cdot \psi_{,i} = 0 45 \end{equation} 46 on faces $x_0=x^{min}_0$, 47 $x_0=x^{max}_0$, 48 $x_1=x^{min}_1$ or 49 $x_1=x^{max}_1$. On the bottom surface we set 50 \begin{equation}\label{GRAV:EQU:101b} 51 \psi = 0 \mbox{ for } x_2=x^{min}_2 52 \end{equation} 53 which sets all horizontal gravity forces to zero $g_0=g_1=0$. 54 For the boundary conditions on the top surface we consider two options 55 \begin{itemize} 56 \item[(D)] 57 On the top surface we set 58 \begin{equation}\label{GRAV:EQU:101 BC D} 59 \psi = 0 \mbox{ for } x_2=x^{max}_2 60 \end{equation} 61 which sets all horizontal gravity forces to zero $g_0=g_1=0$. 62 \item[(N)] 63 On the top surface we set the vertical gravity forces to zero: 64 \begin{equation}\label{GRAV:EQU:101 BC N} 65 n_i \cdot \psi_{,i} = 0 \mbox{ for } x_2=x^{max}_2 66 \end{equation} 67 \end{itemize} 68 The PDE~\ref{GRAV:EQU:100} together with the boundary conditions~\ref{GRAV:EQU:101aa}, ~\ref{GRAV:EQU:101b} 69 and~\ref{GRAV:EQU:101 BC D} or~\ref{GRAV:EQU:101 BC N} 70 has a unique solution. For a given density $\rho$ we denote the unique solution $\psi$ by $\Psi[\rho]$. 71 72 73 The problem is to calculate the density distribution $\rho$ from the gravity force known at some parts of the region of interest 74 $\Omega$. In fact we want to minimize the value 75 \begin{equation}\label{GRAV:EQU:102} 76 J_{data}(\psi) = \frac{1}{2}\sum_{s} \int_{\Omega} \chi^{(s)}_i \cdot ( g_{i}- g^{(s)}_i)^2 dx 77 \footnote{Summation over $i$ is performed by Einstein notation.} 78 \end{equation} 79 where $g^{(s)}_i$ is a measurement of the gravity force $g_i=-\psi_{,i}$ and $\chi^{(s)}_i$ is a weighting factor. 80 $s$ indexes the surveys included in the inversion. 81 We assume that the gravity force $g^{(s)}_i$ is measured as deviation from a background gravity force $g^{bg}_i$, 82 i.e. the actually present gravity force is given as $g^{(s)}_i + g^{bg}$~\footnote{Notice that we use the same reference gravity force for all surveys.}. 83 This can be added to the 84 gravitational potential $\psi$ 85 \begin{equation}\label{GRAV:EQU:102x} 86 \psi^{total} = \psi + g^{bg}_i \cdot ( x_i - x^{min}_i ) 87 \end{equation} 88 The total gravitational potential $\psi^{total}$ still fulfills the PDE~\ref{GRAV:EQU:100} 89 with inhomogeneous versions of the boundary conditions. For the most common case of a purely vertical 90 background gravity force (i.e. $g^{bg}_0=g^{bg}_1=0$) the boundary conditions at the top of the domain 91 take the form $\psi^{total} =g^{bg}_2 \cdot ( x_2^{max} - x^{min}_2 )$ for case (D) or 92 $n_i \cdot \psi^{total}_{,i} =g^{bg}_2$ for case (N). Notice that the first case does not enforce 93 the background gravity force on the top surface but enforces the vertical gravity forces to be zero. 94 The second condition for case (D) enforces the vertical gravity force to the background gravity force but allows 95 for a non-zero vertical gravity force. In most cases the case (N) is the appropriate boundary condition 96 assuming the top boundary of the domain is sufficiently far away from any density variation. 97 98 99 In field observations $g^{(s)}_i$ 100 are known on a sub-domain of the domain only so one will set 101 \begin{equation}\label{GRAV:EQU:105} 102 \chi^{(s)}_i 103 = \left\{ 104 \begin{array}{lcl} 105 \frac{1}{(\sigma^{(s)}_i)^2} & & g^{(s)}_i \mbox{ is available} \\ 106 & \mbox{ where } & \\ 107 0 & & \mbox{ otherwise } \\ 108 \end{array} 109 \right. 110 \end{equation} 111 where $\sigma^{(s)}_i$ is the standard error deviation of measurement $g^{(s)}_i$, see~\cite{A}. 112 With this setting the data $g^{(s)}_i$ can be extended to any arbitrary value (typically to zero) on 113 all locations where no measurements of $g^{(s)}_i$ are available. 114 115 We need to make sure that the density $\rho$ takes physically meaningful values. 116 To achieve this we use a parametrization $C(m)$ for the $\rho$ where $C$ is a given function. We think as 117 $m$ being a dimensionless value between $-\infty$ to $\infty$. Possible settings for the function $C$ are: 118 \begin{itemize} 119 \item[(I)] We can set $\rho = C(m)= \rho^{scale} \cdot m$ where $\rho^{scale}$ is a reference density. 120 \item[(B)] To restrict the density $\rho$ between 121 a minimum density $\rho^{min}$ and a maximum density $\rho^{max}$ one sets 122 \begin{equation} 123 \rho = C(m) = \frac{\rho^{max} +\rho^{min}}{2} + \frac{\rho^{max} -\rho^{min}}{2} \cdot \tanh(m) 124 \end{equation} 125 \item[(L)] To use a logarithmic scale which ensures a positive density one sets 126 \begin{equation} 127 \rho = C(m) = \rho^{scale} \cdot e^{m} 128 \end{equation} 129 \end{itemize} 130 In the following the exact form of the parametrization $C$ is not relevant, however we will make use of the 131 fact that $C$ is sufficiently smooth and invertible, i.e. for any given valid density $\rho$ there is a unique 132 value $m$ in the parameter space such that $\rho=C(m)$. For the inverse function we use the notation $C^{-1}$. 133 In the minimization problem we will search for the value of $m$ providing the best fit to the data 134 rather than the density $\rho$. 135 136 137 The regularization term is added into the minimization problem so the problem has a unique answer. This 138 takes the form 139 \begin{equation}\label{GRAV:EQU:104} 140 J_{reg}(m) = \int_{\Omega} \omega \cdot (m-m^{ref})^2 + \omega_i \cdot L_i^2 \cdot [ (m -m^{ref})_{,i}] ^2 dx 141 \end{equation} 142 where $\omega$ and $\omega_i$ are weighting factors, $L_i$ is the smoothing length (e.g. $L_i= x^{max}_i-x^{min}_i$) 143 and $m^{ref}=C^{-1}(\rho^{ref})$ defines a reference parametrization for a given reference 144 density $\rho^{ref}$. 145 If the density (and therefore the parameter $m$) is known to be exactly $\rho^{ref}$ on a certain subregion $\Omega^{ref}$ of 146 the domain this information can be used 147 to directly constraint the problem, or alternatively one can choose the weights $\omega$. 148 149 The task is now to find $m$ which minimizes the cost function 150 \begin{equation}\label{GRAV:EQU:103} 151 f(m) = J_{data}(\Psi[C[m]]) + \mu \cdot J_{reg}(m) 152 \end{equation} 153 where $\mu>0$ is a fixed penalty factor subject to $m=m^{ref}$ on $\Omega^{ref}$. 154 155 \section{Solution methods} 156 To apply the Nonlinear Conjugate Gradient method (NLCG), see Appendix~\ref{sec:NLCG} or the L-BFGS method, see Appendix~\ref{sec:LBFGS} we need 157 to define an inner product $<.,.>$ and need to calculate the gradient of of $f$. 158 159 As inner product we use 160 \begin{equation}\label{GRAV:EQU:200} 161 = \int_{\Omega} p \cdot q \; dx 162 \end{equation} 163 With this notation the gravity potential $\phi$ is given as 164 \begin{equation}\label{GRAV:EQU:201} 165 < q_{,i},\psi_{,i} > = - 4\pi \cdot G \cdot < q , \rho > \mbox{ for all } q \mbox{ with } q=0 \mbox{ on } \Gamma_{z} 166 \end{equation} 167 where $\Gamma_{z}$ denotes the top and bottom surface of the domain for case (D) 168 and the top surface of the domain for case (N). 169 170 171 For any $q$ with $q=0$ on $\Omega^{ref}$ we get 172 \begin{equation}\label{GRAV:EQU:202} 173 < \nabla f(m),q> = < \nabla J_{data}(\Psi[C[m]]), q> + \mu \cdot < \nabla J_{reg}(m),q> 174 \end{equation} 175 with 176 \begin{equation}\label{GRAV:EQU:202a} 177 < \nabla J_{reg}(m),q> = 178 \int_{\Omega} 2 \omega q \cdot (m-m^{ref}) + 2 \omega_i \cdot L_i^2 \cdot q_{,i} (m-m^{ref})_{,i} dx 179 \end{equation} 180 and 181 \begin{equation}\label{GRAV:EQU:202b} 182 < \nabla J_{data}(\Psi[C[m]]), q> = 183 \sum_{s} \int_{\Omega} \chi^{(s)}_i \cdot ( \psi_{,i} + g^{(s)}_i) \cdot (\Psi[\frac{dC}{dm} \cdot q])_{,i} dx 184 \end{equation} 185 where $\rho=C(m)$ and $\psi=\Psi[\rho]$. We then set 186 \begin{equation}\label{GRAV:EQU:202c} 187 Z_i[\psi]= \sum_{s} \chi^{(s)}_i \cdot ( \psi_{,i} + g^{(s)}_i) 188 \end{equation} 189 and $Z^*[\psi]$ as the solution of the equation 190 \begin{equation}\label{GRAV:EQU:202d} 191 < p_{,i},Z^*[\psi]_{,i} > = < p_{,i} ,Z_i[\psi] > \mbox{ for all } p \mbox{ with } p=0 \mbox{ on } \Gamma_{z} 192 \end{equation} 193 with $Z^*[\psi]=0$ on $\Gamma_{z}$. By stetting $p=\Psi_0[\frac{dC}{dm} \cdot q]$ we get 194 \begin{equation}\label{GRAV:EQU:202e} 195 <(\Psi_0[\frac{dC}{dm} \cdot q]) _{,i},Z^*[\psi]_{,i} > = < (\Psi_0[\frac{dC}{dm} \cdot q])_{,i} ,Z_i[\psi] > 196 \end{equation} 197 and from~\ref{GRAV:EQU:201} with $q=Z^*[\psi]$ we get 198 \begin{equation}\label{GRAV:EQU:20e} 199 < (Z^*[\psi])_{,i},(\Psi_0[\frac{dC}{dm} \cdot q])_{,i} > = - 4\pi \cdot G \cdot < Z^*[\psi], \frac{dC}{dm} \cdot q > 200 \end{equation} 201 which leads to 202 \begin{equation}\label{GRAV:EQU:202f} 203 < \nabla J_{data}(\Psi[C[m]]), q> = < (\Psi_0[\frac{dC}{dm} \cdot q])_{,i} ,Z_i[\psi] > 204 = - 4\pi \cdot G \cdot < \frac{dC}{dm} Z^*[\psi] , q > 205 \end{equation} 206 Putting this all together we get 207 \begin{equation}\label{GRAV:EQU:202g} 208 < \nabla f(m),q> = 209 \mu \cdot < 2 \omega \cdot (m-m^{ref}), q> 210 + \mu \cdot < 2 \omega_i \cdot L_i^2 \cdot (m-m^{ref})_{,i}, q_{,i}> - 4\pi \cdot G \cdot < \frac{dC}{dm} Z^*[\psi] , q > 211 \end{equation} 212