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

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 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 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}