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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4376 - (hide annotations)
Mon Apr 22 08:44:45 2013 UTC (6 years, 1 month ago) by gross
File MIME type: application/x-tex
File size: 10077 byte(s)
domumentation update on coordinate systems - part I
1 caltinay 4045
2     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 jfenwick 4154 % Copyright (c) 2003-2013 by University of Queensland
4 caltinay 4045 % 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 gross 4055 \section{Gravity Inversion}\label{sec:forward gravity}
16 gross 4373 For the gravity inversion we use the anomaly of the gravity acceleration~\index{gravity acceleration} of the Earth.
17 caltinay 4095 The controlling material parameter is the density~\index{density} $\rho$ of
18     the rock.
19     If the density field $\rho$ is known the gravitational potential $\psi$ is
20     given as the solution of the PDE
21 gross 4093 \begin{equation}\label{ref:GRAV:EQU:100}
22     -\psi_{,ii} = -4\pi G \cdot \rho
23 caltinay 4045 \end{equation}
24 caltinay 4095 where $G=6.6730 \cdot 10^{-11} \frac{m^3}{kg \cdot s^2}$ is the gravitational
25     constant.
26 gross 4125 The gravitational potential is set to zero at the top of the
27 gross 4121 domain $\Gamma_0$.
28 caltinay 4095 On all other faces the normal component of the gravity acceleration anomaly
29     $g_i$ is set to zero, i.e. $n_i \psi_{,i} = 0$ with outer normal field $n_i$.
30     The gravity force $g_i$ is given as the negative of the gradient of the gravity
31     potential $\psi$:
32 gross 4093 \begin{equation}\label{ref:GRAV:EQU:101}
33 caltinay 4045 g_i = - \psi_{,i}
34     \end{equation}
35 caltinay 4095 From the gravitational potential we can calculate the gravity acceleration
36     anomaly via Equation~(\ref{ref:GRAV:EQU:101}) to obtain the defect to the
37     given data.
38     If $g^{(s)}_i$ is a measurement of the gravity acceleration anomaly for
39     survey $s$ and $\omega^{(s)}_i$ is a weighting factor the data defect
40 gross 4142 $J^{grav}(k)$ in the notation of Chapter~\ref{chapter:ref:inversion cost function} is given as
41 gross 4093 \begin{equation}\label{ref:GRAV:EQU:9}
42 gross 4099 J^{grav}(k) = \frac{1}{2}\sum_{s} \int_{\Omega} ( \omega^{(s)}_i \cdot (g_{i}- g^{(s)}_i) ) ^2 dx
43 caltinay 4045 \end{equation}
44 gross 4099 Summation over $i$ is performed.
45 gross 4093 The cost function kernel\index{cost function!kernel} is given as
46     \begin{equation}\label{ref:GRAV:EQU:10}
47 gross 4099 K^{grav}(\psi_{,i},k) = \frac{1}{2}\sum_{s} ( \omega^{(s)}_i \cdot (\psi_{,i}+ g^{(s)}_i) ) ^2
48 caltinay 4045 \end{equation}
49 caltinay 4095 In practice the gravity acceleration $g^{(s)}$ is measured in vertical
50     direction $z$ with a standard error deviation $\sigma^{(s)}$ at certain
51     locations in the domain.
52     In this case one sets the weighting factors $\omega^{(s)}$ as
53 gross 4093 \begin{equation}\label{ref:GRAV:EQU:11}
54     \omega^{(s)}_i
55 caltinay 4045 = \left\{
56     \begin{array}{lcl}
57 gross 4093 f \cdot \frac{\delta_{iz}}{\sigma^{(s)}} & & \mbox{data are available} \\
58 caltinay 4045 & \mbox{ where } & \\
59     0 & & \mbox{ otherwise } \\
60     \end{array}
61     \right.
62     \end{equation}
63 gross 4102 With the objective to control the
64     gradient of the cost function
65     the scaling factor $f$ is chosen in the way that
66 gross 4093 \begin{equation}\label{ref:GRAV:EQU:12}
67 gross 4102 \sum_{s} \int_{\Omega} ( \omega^{(s)}_i g^{(s)}_i ) \cdot ( \omega^{(s)}_j \frac{1}{L_j} ) \cdot 4\pi G L^2 \cdot \rho' \; dx =\alpha
68 caltinay 4045 \end{equation}
69 gross 4102 where $\alpha$ defines a scaling factor which is typically set to one and $L$ is defined by equation~(\ref{ref:EQU:REG:6b}). $\rho'$ is considering the
70     derivative of the density with respect to the level set function.
71 caltinay 4045
72 gross 4102
73 gross 4093 \subsection{Usage}
74 caltinay 4045
75 gross 4093 \LG{Add example}
76 caltinay 4045
77 gross 4099 \begin{classdesc}{GravityModel}{domain,
78     w, g,
79 gross 4376 \optional{, coordinates=\None}
80 gross 4344 \optional{, fixPotentialAtBottom=False},
81 gross 4102 \optional{, tol=1e-8}
82     }
83 gross 4099 opens a gravity forward model over the \Domain \member{domain} with
84 gross 4373 weighting factors \member{w} ($=\omega^{(s)}$) and measured gravity acceleration anomalies \member{g} ($=g^{(s)}$).
85     The weighting factors and the measured gravity acceleration anomalies must be vectors
86     where components refer to the components
87     $(x_0,x_1,x_2)$ for the Cartesian coordinate system
88     and to $(\phi, \lambda, h)$ for the geodetic coordinate system.
89     If \member{reference} defines the reference coordinate system to be used, see Section~\ref{sec:cooridiate systems}.
90 gross 4099 \member{tol} set the tolerance for the solution of the PDE~(\ref{ref:GRAV:EQU:100}).
91 gross 4344 If \member{fixPotentialAtBottom} is set to \True, the gravitational potential
92     at the bottom is set to zero in addition to the potential on the top.
93 gross 4376 \member{coordinates} set the reference coordinate system to be used. By the default the
94     Cartesian coordinate system is used.
95 gross 4093 \end{classdesc}
96 caltinay 4045
97 gross 4102 \begin{methoddesc}[GravityModel]{rescaleWeights}{
98     \optional{scale=1.}
99     \optional{rho_scale=1.}}
100     rescale the weighting factors such condition~(\ref{ref:GRAV:EQU:12}) holds where
101     \member{scale} sets the scale $\alpha$
102     and \member{rho_scale} sets $\rho'$. This method should be called before any inversion is started
103     in order to make sure that all components of the cost function are appropriately scaled.
104     \end{methoddesc}
105    
106    
107 gross 4093 \subsection{Gradient Calculation}
108 caltinay 4095 This section briefly explains how the gradient
109     $\frac{\partial J^{grav}}{\partial \rho}$ of the cost function $J^{grav}$ with
110 gross 4142 respect to the density $\rho$ is calculated. We follow the concept as outlined in section~\ref{chapter:ref:inversion cost function:gradient}.
111 gross 4373 The gravity potential $\psi$ from PDE~(\ref{ref:GRAV:EQU:100}) is solved in
112 caltinay 4095 weak form:
113 gross 4093 \begin{equation}\label{ref:GRAV:EQU:201}
114     \int_{\Omega} q_{,i} \psi_{,i} \; dx = - \int_{\Omega} 4\pi G \cdot q \rho\; dx
115 caltinay 4045 \end{equation}
116 gross 4121 for all $q$ with $q=0$ on $\Gamma_0$.
117 caltinay 4095 In the following we set $\Psi[\cdot]=\psi$ for a given density $\cdot$ as
118     solution of the variational problem~(\ref{ref:GRAV:EQU:201}).
119     If $\Gamma_{\rho}$ denotes the region of the domain where the density is known
120 gross 4093 and for a given direction $p$ with $p=0$ on $\Gamma_{\rho}$ one has
121 gross 4121 \begin{equation}\label{ref:GRAV:EQU:201aa}
122 gross 4099 \int_{\Omega} \frac{\partial J^{grav}}{\partial \rho} \cdot p \; dx = \int_{\Omega}
123 gross 4093 \sum_{s} (\omega^{(s)}_j \cdot
124     (g^{(s)}_j-g_{j}) ) \cdot ( \omega^{(s)}_i \Psi[p]_{,i}) \; dx
125 caltinay 4045 \end{equation}
126 caltinay 4095 with
127 gross 4093 \begin{equation}\label{ref:GRAV:EQU:202c}
128 gross 4099 Y_i[\psi]= \sum_{s} (\omega^{(s)}_j \cdot
129 gross 4093 (g^{(s)}_j-g_{j}) ) \cdot \omega^{(s)}_i
130 caltinay 4045 \end{equation}
131 gross 4093 This is written as
132 gross 4121 \begin{equation}\label{ref:GRAV:EQU:202cc}
133 gross 4093 \int_{\Omega} \frac{\partial J^{grav}}{\partial \rho} \cdot p \; dx = \int_{\Omega}
134     Y_i[\psi] \Psi[p]_{,i} \; dx
135 caltinay 4045 \end{equation}
136 gross 4093 We then set $Y^*[\psi]$ as the solution of the equation
137     \begin{equation}\label{ref:GRAV:EQU:202d}
138     \int_{\Omega} r_{,i} Y^*[\psi]_{,i} \; dx = \int_{\Omega} r_{,i} ,Y_i[\psi] \; dx \mbox{ for all } p \mbox{ with } r=0 \mbox{ on } \Gamma_{top}
139 caltinay 4045 \end{equation}
140 gross 4121 with $Y^*[\psi]=0$ on $\Gamma_0$. With $r=\Psi[p]$ we get
141     \begin{equation}\label{ref:GRAV:EQU:202dd}
142 gross 4093 \int_{\Omega} \Psi[p]_{,i} Y^*[\psi]_{,i} \; dx = \int_{\Omega} \Psi[p]_{,i} ,Y_i[\psi] \; dx
143 caltinay 4045 \end{equation}
144 caltinay 4095 and from Equation~(\ref{ref:GRAV:EQU:201}) with $q=Y^*[\psi]$ we get
145 gross 4093 \begin{equation}\label{ref:GRAV:EQU:20e}
146     \int_{\Omega} Y^*[\psi]_{,i} \Psi[p]_{,i} \; dx = - \int_{\Omega} 4\pi G \cdot Y^*[\psi] \cdot p\; dx
147 caltinay 4045 \end{equation}
148     which leads to
149 gross 4093 \begin{equation}\label{ref:GRAV:EQU:20ee}
150     \int_{\Omega} \Psi[p]_{,i} ,Y_i[\psi] \; dx = - \int_{\Omega} 4\pi G \cdot Y^*[\psi] \cdot p \; dx
151     \end{equation}
152     and finally
153     \begin{equation}\label{ref:GRAV:EQU:201a}
154     \int_{\Omega} \frac{\partial J^{grav}}{\partial \rho} \cdot p \; dx = - \int_{\Omega}
155     4\pi G \cdot Y^*[\psi] \cdot p \; dx
156 caltinay 4045 \end{equation}
157 gross 4093 or
158     \begin{equation}\label{ref:GRAV:EQU:201b}
159     \frac{\partial J^{grav}}{\partial \rho} =- 4\pi G \cdot Y^*[\psi]
160     \end{equation}
161 caltinay 4045
162 gross 4373 \subsection{Geodetic Coordinates }
163 gross 4376 For geodetic coordinates $(\phi, \lambda, h)$, see Chapter~\ref{Chp:ref:coordinates}, the solution process needs to be slightly modified.
164 gross 4373 Observations are recorded along the geodetic coordinates axes $\alpha$ rather than the Cartesian axes $i$. In fact we
165     have in equation~\ref{ref:GRAV:EQU:9}:
166     \begin{equation}\label{ref:GRAV:EQU:300}
167     \omega^{(s)}_i \cdot (g_{i}- g^{(s)}_i) = \omega^{(s)}_{\alpha} \cdot (g_{{\alpha}}- g^{(s)}_{\alpha})
168     \end{equation}
169     where now $g^{(s)}_{\alpha}$ are the observational data with weighting factors $\omega^{(s)}_{\alpha}$. Using the
170 gross 4376 fact that $g_{{\alpha}} = - d_{\alpha \alpha} \psi_{,\alpha}$
171 gross 4373 equation~\ref{ref:GRAV:EQU:10} translates to
172     \begin{equation}\label{ref:GRAV:EQU:301}
173     J^{grav}(k) = \frac{1}{2}\sum_{s} \int_{\widehat{\Omega}}
174 gross 4376 ( \omega^{(s)}_{\alpha} \cdot (d_{\alpha \alpha} \psi_{,\alpha} + g^{(s)}_{\alpha} ) ) ^2 \; v \; d\widehat{x}
175 gross 4373 \end{equation}
176     where $\widehat{\Omega}$ and $d\widehat{x}$ refer to integration over the geodetic coordinates axes. This can be rearranged to
177     \begin{equation}\label{ref:GRAV:EQU:301}
178     J^{grav}(k) = \frac{1}{2}\sum_{s} \int_{\widehat{\Omega}}
179 gross 4376 ( \omega^{(s)}_{\alpha} v^{\frac{1}{2}} d_{\alpha \alpha} \cdot ( \psi_{,\alpha} + \frac{1}{d_{\alpha \alpha}} g^{(s)}_{\alpha} ) ) ^2 \; d\widehat{x}
180 gross 4373 =\frac{1}{2}\sum_{s} \int_{\widehat{\Omega}}
181     ( {\widehat{\omega}}^{(s)}_{\alpha}\cdot ( \psi_{,\alpha} + \widehat{g}^{(s)}_{\alpha} ) ) ^2 \; d\widehat{x}
182     \end{equation}
183     with
184     \begin{equation}\label{ref:GRAV:EQU:301b}
185 gross 4376 \widehat{\omega}^{(s)}_{\alpha} =\omega^{(s)}_{\alpha} v^{\frac{1}{2}} d_{\alpha \alpha} \mbox{ and }
186 gross 4373 \widehat{ g}^{(s)}_{\alpha}=
187 gross 4376 \frac{1}{d_{\alpha \alpha}} g^{(s)}_{\alpha}
188 gross 4373 \end{equation}
189     which means one can apply the Cartesian formulation to the geodetic coordinates using modified data.
190     The gravity potential is calculated from
191     \begin{equation}\label{ref:GRAV:EQU:302}
192 gross 4376 \int_{\widehat{\Omega}} v \; d_{\alpha \alpha}^2 q_{,\alpha} \psi_{,\alpha} \; d\widehat{x}
193 gross 4373 = - \int_{\widehat{\Omega}} (4\pi G v) \cdot q \rho\; d\widehat{x}
194     \end{equation}
195     see equation~\ref{ref:GRAV:EQU:201}, and the adjoint function $Y^*[\psi]$ for $Y_{\alpha}[\psi]$ is given from
196     \begin{equation}\label{ref:GRAV:EQU:303}
197 gross 4376 \int_{\widehat{\Omega}} v \; d_{\alpha \alpha}^2 q_{,\alpha} Y^*[\psi]_{,\alpha } \;d\widehat{x} =
198 gross 4373 \int_{\widehat{\Omega}} q_{,\alpha} Y_{\alpha}[\psi] \; d\widehat{x}
199     \end{equation}
200     and finally
201     \begin{equation}\label{ref:GRAV:EQU:201a}
202     \frac{\partial J^{grav}}{\partial \rho} = -
203     (4\pi G v) \cdot Y^*[\psi]
204     \end{equation}

  ViewVC Help
Powered by ViewVC 1.1.26