# Contents of /branches/doubleplusgood/doc/inversion/ForwardMagnetic.tex

Revision 4345 - (show annotations)
Fri Mar 29 07:09:41 2013 UTC (5 years, 11 months ago) by jfenwick
File MIME type: application/x-tex
File size: 8623 byte(s)
Spelling fixes

 1 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % Copyright (c) 2003-2013 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 16 \section{Linear Magnetic Inversion}\label{sec:forward magnetic} 17 For the magnetic inversion we use the anomaly of the magnetic flux 18 density~\index{magnetic flux density} of the Earth. 19 The controlling material parameter is the susceptibility~\index{susceptibility} 20 $k$ of the rock. 21 With magnetization $M$ and inducing magnetic field anomaly $H^s$, the magnetic 22 flux density anomaly $B^s$ is given as 23 \begin{equation}\label{ref:MAG:EQU:1} 24 B_i = \mu_0 \cdot ( H^s_i + M_i ) 25 \end{equation} 26 where $\mu_0 = 4 \pi \cdot 10^{-7} \frac{Vs}{Am}$. 27 In this forward model we make the simplifying assumption that the magnetization 28 is proportional to the known geomagnetic flux density $B^b$: 29 \begin{equation}\label{ref:MAG:EQU:4} 30 \mu_0 \cdot M_i = k \cdot B^b_i \;. 31 \end{equation} 32 Values for the magnetic flux density can be obtained by the International 33 Geomagnetic Reference Field (IGRF)~\cite{IGRF} 34 (or the Australian Geomagnetic Reference Field (AGRF)~\cite{AGRF}). 35 In most cases it is reasonable to assume that that the background field is 36 constant across the domain. 37 38 The magnetic field anomaly $H^s$ can be represented by the gradient of a 39 magnetic scalar potential\index{scalar potential!magnetic} $\psi$. 40 We use the form 41 \begin{equation}\label{ref:MAG:EQU:6} 42 \mu_0 \cdot H^s_i = - \psi_{,i} 43 \end{equation} 44 With this notation one gets from Equations~(\ref{ref:MAG:EQU:1}) and~(\ref{ref:MAG:EQU:4}): 45 \begin{equation}\label{ref:MAG:EQU:7} 46 B_i = - \psi_{,i} + k \cdot B^b_i 47 \end{equation} 48 As the $B^s$ magnetic flux density anomaly we obtain the PDE 49 \begin{equation}\label{ref:MAG:EQU:8} 50 - \psi_{,ii} = - (k B^b_i)_{,i} 51 \end{equation} 52 which needs to be solved for a given susceptibility $k$. 53 The magnetic scalar potential is set to zero at the top of the domain 54 $\Gamma_{0}$. 55 On all other faces the normal component of the magnetic flux density anomaly 56 $B_i$ is set to zero, i.e. $n_i \psi_{,i} = k \cdot n_i B^b_i$ with outer 57 normal field $n_i$. 58 59 From the magnetic scalar potential we can calculate the magnetic flux density 60 anomaly via Equation~(\ref{ref:MAG:EQU:8}) to calculate the defect to the given 61 data. 62 If $B^{(s)}_i$ is a measurement of the magnetic flux density anomaly for 63 survey $s$ and $\omega^{(s)}_i$ is a weighting factor the data defect 64 $J^{mag}(k)$ in the notation of Chapter~\ref{chapter:ref:inversion cost function} is given as 65 \begin{equation}\label{ref:MAG:EQU:9} 66 J^{mag}(k) = \frac{1}{2}\sum_{s} \int_{\Omega} ( \omega^{(s)}_i \cdot (B_{i}- B^{(s)}_i) ) ^2 dx 67 \end{equation} 68 Summation over $i$ is performed. 69 The cost function kernel\index{cost function!kernel} is given as 70 \begin{equation}\label{ref:MAG:EQU:10} 71 K^{mag}(\psi_{,i},k) = \frac{1}{2}\sum_{s} ( \omega^{(s)}_i \cdot (k \cdot B^b_i - \psi_{,i} - B^{(s)}_i) ) ^2 72 \end{equation} 73 Notice that if magnetic flux density is measured in air one can ignore the 74 $k\cdot B^b_i$ as the susceptibility is zero. 75 76 In practice the magnetic flux density $b^{(s)}$ is measured along a certain 77 direction $d^{(s)}_i$ with a standard error deviation $\sigma^{(s)}$ at 78 certain locations in the domain. 79 In this case one sets $B^{(s)}_i=b^{(s)} \cdot d^{(s)}_i$ and the weighting 80 factors $\omega^{(s)}$ as 81 \begin{equation}\label{ref:MAG:EQU:11} 82 \omega^{(s)}_i 83 = \left\{ 84 \begin{array}{lcl} 85 f \cdot \frac{d^{(s)}_i}{\sigma^{(s)}} & & \mbox{data are available} \\ 86 & \mbox{ where } & \\ 87 0 & & \mbox{ otherwise } \\ 88 \end{array} 89 \right. 90 \end{equation} 91 where it is assumed that $d^{(s)}_i \cdot d^{(s)}_i =1$. With the objective to control the 92 gradient of the cost function the scaling factor $f$ is chosen in the way that 93 \begin{equation}\label{ref:MAG:EQU:12} 94 \sum_{s} \int_{\Omega} ( \omega^{(s)}_i B^{(s)}_i ) 95 \cdot ( \omega^{(s)}_j \frac{1}{L_j} ) \cdot L^2 \cdot 96 ( B^b_n \frac{1}{L_n} ) 97 \cdot k' \; 98 dx =\alpha 99 \end{equation} 100 where $\alpha$ defines a scaling factor which is typically set to one and $L$ is defined by equation~(\ref{ref:EQU:REG:6b}). 101 $k'$ is considering the 102 derivative of the density with respect to the level set function. 103 104 \subsection{Usage} 105 106 \LG{Add example} 107 108 \begin{classdesc}{MagneticModel}{domain, w, B, background_field, 109 \optional{, useSphericalCoordinates=False} 110 \optional{, fixPotentialAtBottom=False}, 111 \optional{, tol=1e-8}} 112 opens a magnetic forward model over the \Domain \member{domain} with 113 weighting factors \member{w} ($=\omega^{(s)}$) and measured magnetic flux 114 density anomalies \member{B} ($=B^{(s)}$). 115 The weighting factors and the measured magnetic flux density anomalies must be vectors. 116 \member{background_field} defines the background magnetic flux density $B^b$ 117 as a vector with north, east and vertical components. 118 If \member{useSphericalCoordinates} is \True spherical coordinates are used. 119 \member{tol} sets the tolerance for the solution of the PDE~(\ref{ref:MAG:EQU:8}). 120 If \member{fixPotentialAtBottom} is set to \True, the gravitational potential 121 at the bottom is set to zero in addition to the potential on the top. 122 \end{classdesc} 123 124 \begin{methoddesc}[MagneticModel]{rescaleWeights}{ 125 \optional{scale=1.} 126 \optional{k_scale=1.}} 127 rescale the weighting factors such condition~(\ref{ref:MAG:EQU:12}) holds where 128 \member{scale} sets the scale $\alpha$ 129 and \member{k_scale} sets $k'$. This method should be called before any inversion is started 130 in order to make sure that all components of the cost function are appropriately scaled. 131 \end{methoddesc} 132 133 134 \subsection{Gradient Calculation} 135 This section briefly explains how the gradient 136 $\frac{\partial J^{mag}}{\partial k}$ of the cost function $J^{mag}$ with 137 respect to the susceptibility $k$ is calculated. We follow the concept as outlined in section~\ref{chapter:ref:inversion cost function:gradient}. 138 139 The magnetic potential $\psi$ from PDE~(\ref{ref:MAG:EQU:8}) is solved in weak form: 140 \begin{equation}\label{ref:MAG:EQU:201} 141 \int_{\Omega} q_{,i} \psi_{,i} \; dx = \int_{\Omega} k \cdot q_{,i} B^b_i \; dx 142 \end{equation} 143 for all $q$ with $q=0$ on $\Gamma_{0}$. 144 In the following we set $\Psi[k]=\psi$ for a given susceptibility $k$ as 145 solution of the variational problem~(\ref{ref:MAG:EQU:201}). 146 If $\Gamma_{k}$ denotes the region of the domain where the susceptibility is 147 known and for a given direction $p$ with $p=0$ on $\Gamma_{k}$ one has 148 \begin{equation}\label{ref:MAG:EQU:201aa} 149 \int_{\Omega} \frac{\partial J^{mag}}{\partial k} \cdot p \; dx = \int_{\Omega} 150 \sum_{s} (\omega^{(s)}_j 151 ( B^{(s)}_j-B_{j})) \cdot ( \omega^{(s)}_i ( \Psi[p]_{,i} - p \cdot B^b_i ) ) \; dx 152 \end{equation} 153 with 154 \begin{equation}\label{ref:MAG:EQU:202c} 155 Y_i[\psi]= \sum_{s} (\omega^{(s)}_j 156 (B^{(s)}_j - B_{j}) ) \cdot \omega^{(s)}_i 157 \end{equation} 158 This is written as 159 \begin{equation}\label{ref:MAG:EQU:202cc} 160 \int_{\Omega} \frac{\partial J^{mag}}{\partial k} \cdot p \; dx = \int_{\Omega} 161 Y_i[\psi] \Psi[p]_{,i} - p \cdot Y_i[\psi]B^b_i \; dx 162 \end{equation} 163 We then set $Y^*[\psi]$ as the solution of the equation 164 \begin{equation}\label{ref:MAG:EQU:202d} 165 \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_{0} 166 \end{equation} 167 with $Y^*[\psi]=0$ on $\Gamma_{0}$. With $r=\Psi[p]$ we get 168 \begin{equation}\label{ref:MAG:EQU:202dd} 169 \int_{\Omega} \Psi[p]_{,i} Y^*[\psi]_{,i} \; dx = \int_{\Omega} \Psi[p]_{,i} ,Y_i[\psi] \; dx 170 \end{equation} 171 and from Equation~(\ref{ref:MAG:EQU:201}) with $q=Y^*[\psi]$ we get 172 \begin{equation}\label{ref:MAG:EQU:20e} 173 \int_{\Omega} Y^*[\psi]_{,i} \Psi[p]_{,i} \; dx = \int_{\Omega} p \cdot Y^*[\psi]_{,i} B^b_i \; dx 174 \end{equation} 175 which leads to 176 \begin{equation}\label{ref:MAG:EQU:20ee} 177 \int_{\Omega} \Psi[p]_{,i} ,Y_i[\psi] \; dx = \int_{\Omega} p \cdot Y^*[\psi]_{,i} B^b_i \; dx 178 \end{equation} 179 and finally 180 \begin{equation}\label{ref:MAG:EQU:201a} 181 \int_{\Omega} \frac{\partial J^{mag}}{\partial k} \cdot p \; dx = \int_{\Omega} 182 p \cdot (Y_i[\psi] - Y^*[\psi]_{,i}) B^b_i \; dx 183 \end{equation} 184 or 185 \begin{equation}\label{ref:MAG:EQU:201b} 186 \frac{\partial J^{mag}}{\partial k} = (Y^*[\psi]_{,i}-Y_i[\psi]) B^b_i 187 \end{equation} 188