/[escript]/trunk/doc/inversion/gravity.tex
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


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 % http://www.uq.edu.au
5 %
6 % Primary Business: Queensland, Australia
7 % Licensed under the Open Software License version 3.0
8 % http://www.opensource.org/licenses/osl-3.0.php
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 <p,q> = \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

  ViewVC Help
Powered by ViewVC 1.1.26