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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4051 - (hide annotations)
Wed Oct 31 03:40:47 2012 UTC (6 years, 11 months ago) by caltinay
File MIME type: application/x-tex
File size: 10327 byte(s)
Split inversion cookbook into two parts...

1 caltinay 4045
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 caltinay 4051 \chapter{Gravity Inversion}\label{chp:gravityinversion}
16 caltinay 4045 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