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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4376 - (show annotations)
Mon Apr 22 08:44:45 2013 UTC (6 years ago) by gross
File MIME type: application/x-tex
File size: 11773 byte(s)
domumentation update on coordinate systems - part I
1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % Copyright (c) 2003-2013 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
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^r_i)_{,i}
51 \end{equation}
52 with $B^r_i=B^b_i$ 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
105 \subsection{Usage}
106
107 \LG{Add example}
108
109 \begin{classdesc}{MagneticModel}{domain, w, B, background_field,
110 \optional{, coordinates=\None}
111 \optional{, fixPotentialAtBottom=False},
112 \optional{, tol=1e-8},
113 }
114 opens a magnetic forward model over the \Domain \member{domain} with
115 weighting factors \member{w} ($=\omega^{(s)}$) and measured magnetic flux
116 density anomalies \member{B} ($=B^{(s)}$).
117 The weighting factors and the measured magnetic flux density anomalies must be vectors.
118 \member{background_field} defines the background magnetic flux density $B^b$
119 as a vector with north, east and vertical components.
120 \member{tol} sets the tolerance for the solution of the PDE~(\ref{ref:MAG:EQU:8}).
121 If \member{fixPotentialAtBottom} is set to \True, the gravitational potential
122 at the bottom is set to zero in addition to the potential on the top.
123 \member{coordinates} set the reference coordinate system to be used. By the default the
124 Cartesian coordinate system is used.
125 \end{classdesc}
126
127 \begin{methoddesc}[MagneticModel]{rescaleWeights}{
128 \optional{scale=1.}
129 \optional{k_scale=1.}}
130 rescale the weighting factors such condition~(\ref{ref:MAG:EQU:12}) holds where
131 \member{scale} sets the scale $\alpha$
132 and \member{k_scale} sets $k'$. This method should be called before any inversion is started
133 in order to make sure that all components of the cost function are appropriately scaled.
134 \end{methoddesc}
135
136
137 \subsection{Gradient Calculation}
138 This section briefly explains how the gradient
139 $\frac{\partial J^{mag}}{\partial k}$ of the cost function $J^{mag}$ with
140 respect to the susceptibility $k$ is calculated. We follow the concept as outlined in section~\ref{chapter:ref:inversion cost function:gradient}.
141
142 The magnetic potential $\psi$ from PDE~(\ref{ref:MAG:EQU:8}) is solved in weak form:
143 \begin{equation}\label{ref:MAG:EQU:201}
144 \int_{\Omega} q_{,i} \psi_{,i} \; dx = \int_{\Omega} k \cdot q_{,i} B^r_i \; dx
145 \end{equation}
146 for all $q$ with $q=0$ on $\Gamma_{0}$.
147 In the following we set $\Psi[k]=\psi$ for a given susceptibility $k$ as
148 solution of the variational problem~(\ref{ref:MAG:EQU:201}).
149 If $\Gamma_{k}$ denotes the region of the domain where the susceptibility is
150 known and for a given direction $p$ with $p=0$ on $\Gamma_{k}$ one has
151 \begin{equation}\label{ref:MAG:EQU:201aa}
152 \int_{\Omega} \frac{\partial J^{mag}}{\partial k} \cdot p \; dx = \int_{\Omega}
153 \sum_{s} (\omega^{(s)}_j
154 ( B^{(s)}_j-B_{j})) \cdot ( \omega^{(s)}_i ( \Psi[p]_{,i} - p \cdot B^b_i ) ) \; dx
155 \end{equation}
156 with
157 \begin{equation}\label{ref:MAG:EQU:202c}
158 Y_i[\psi]= \sum_{s} (\omega^{(s)}_j
159 (B^{(s)}_j - B_{j}) ) \cdot \omega^{(s)}_i
160 \end{equation}
161 This is written as
162 \begin{equation}\label{ref:MAG:EQU:202cc}
163 \int_{\Omega} \frac{\partial J^{mag}}{\partial k} \cdot p \; dx = \int_{\Omega}
164 Y_i[\psi] \Psi[p]_{,i} - p \cdot Y_i[\psi]B^b_i \; dx
165 \end{equation}
166 We then set adjoint function $Y^*[\psi]$ as the solution of the equation
167 \begin{equation}\label{ref:MAG:EQU:202d}
168 \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}
169 \end{equation}
170 with $Y^*[\psi]=0$ on $\Gamma_{0}$. With $r=\Psi[p]$ we get
171 \begin{equation}\label{ref:MAG:EQU:202dd}
172 \int_{\Omega} \Psi[p]_{,i} Y^*[\psi]_{,i} \; dx = \int_{\Omega} \Psi[p]_{,i} Y_i[\psi] \; dx
173 \end{equation}
174 and from Equation~(\ref{ref:MAG:EQU:201}) with $q=Y^*[\psi]$ we get
175 \begin{equation}\label{ref:MAG:EQU:20e}
176 \int_{\Omega} Y^*[\psi]_{,i} \Psi[p]_{,i} \; dx = \int_{\Omega} p \cdot Y^*[\psi]_{,i} B^r_i \; dx
177 \end{equation}
178 which leads to
179 \begin{equation}\label{ref:MAG:EQU:20ee}
180 \int_{\Omega} \Psi[p]_{,i} ,Y_i[\psi] \; dx = \int_{\Omega} p \cdot Y^*[\psi]_{,i} B^r_i \; dx
181 \end{equation}
182 and finally
183 \begin{equation}\label{ref:MAG:EQU:201a}
184 \int_{\Omega} \frac{\partial J^{mag}}{\partial k} \cdot p \; dx = \int_{\Omega}
185 p \cdot (Y^*[\psi]_{,i} B^r_i - Y_i[\psi] B^b_i) \; dx
186 \end{equation}
187 or
188 \begin{equation}\label{ref:MAG:EQU:201b}
189 \frac{\partial J^{mag}}{\partial k} = Y^*[\psi]_{,i} B^r_i - Y_i[\psi] B^b_i
190 \end{equation}
191
192 \subsection{Geodetic Coordinates }
193 For geodetic coordinates $(\phi, \lambda, h)$, see Appendix~\ref{sec:geodetic}, the solution process needs to be slightly modified.
194 Observations are recorded along the geodetic coordinates axes $\alpha$ rather than the Cartesian axes $i$. In fact we
195 have in equation~\ref{ref:MAG:EQU:9}:
196 \begin{equation}\label{ref:MAG:EQU:300}
197 \omega^{(s)}_i \cdot (B_{i}- B^{(s)}_i) = \omega^{(s)}_{\alpha} \cdot (B_{{\alpha}}- B^{(s)}_{\alpha})
198 \end{equation}
199 where now $B^{(s)}_{\alpha}$ are the observational data with weighting factors $\omega^{(s)}_{\alpha}$. Using the
200 fact that $B_{{\alpha}} = k \cdot B^b_{{\alpha}} -d_{\alpha \alpha} \psi_{,\alpha}$
201 equation~\ref{ref:MAG:EQU:10} translates to
202 \begin{equation}\label{ref:MAG:EQU:301}
203 J^{mag}(k) = \frac{1}{2}\sum_{s} \int_{\widehat{\Omega}}
204 ( \omega^{(s)}_{\alpha} \cdot (d_{\alpha \alpha} \psi_{,\alpha} - k \cdot B^b_{{\alpha}} + B^{(s)}_{\alpha} ) ) ^2 \; v \; d\widehat{x}
205 \end{equation}
206 where $\widehat{\Omega}$ and $d\widehat{x}$ refer to integration over the geodetic coordinates axes. This can be rearranged to
207 \begin{equation}\label{ref:MAG:EQU:301}
208 J^{mag}(k) = \frac{1}{2}\sum_{s} \int_{\widehat{\Omega}}
209 ( \omega^{(s)}_{\alpha} v^{\frac{1}{2}} d_{\alpha \alpha} \cdot (
210 \psi_{,\alpha} - k \cdot v_{\alpha \alpha} B^b_{{\alpha}} + v_{\alpha \alpha} B^{(s)}_{\alpha} ) ) ^2 \; d\widehat{x}
211 =\frac{1}{2}\sum_{s} \int_{\widehat{\Omega}}
212 ( {\widehat{\omega}}^{(s)}_{\alpha}\cdot ( \psi_{,\alpha} - k \cdot \widehat{B}^b_{{\alpha}}+ \widehat{B}^{(s)}_{\alpha} ) ) ^2 \; d\widehat{x}
213 \end{equation}
214 with
215 \begin{equation}\label{ref:MAG:EQU:301b}
216 \widehat{\omega}^{(s)}_{\alpha} = \omega^{(s)}_{\alpha} v^{\frac{1}{2}} d_{\alpha \alpha} \mbox{ , }
217 \widehat{B}^{(s)}_{\alpha}=
218 \frac{1}{d_{\alpha \alpha}} B^{(s)}_{\alpha} \mbox{ and } \widehat{B}^b_{{\alpha}} = \frac{1}{d_{\alpha \alpha}} B^b_{{\alpha}}
219 \end{equation}
220 which means one can apply the Cartesian formulation to the geodetic coordinates using modified data.
221
222
223 The magnetic potential is calculated from
224 \begin{equation}\label{ref:MAG:EQU:302}
225 \int_{\widehat{\Omega}} v \; d_{\alpha \alpha}^2 q_{,\alpha} \psi_{,\alpha} \; d\widehat{x}
226 = \int_{\widehat{\Omega}} v \; d_{\alpha \alpha} k \cdot q_{,\alpha} B^r_{\alpha} \; d\widehat{x}
227 = \int_{\widehat{\Omega}} k \cdot q_{,\alpha} \widehat{B}^r_{\alpha} \; d\widehat{x}
228 \end{equation}
229 with
230 \begin{equation}\label{ref:MAG:EQU:302b}
231 \widehat{B}^r_{\alpha} =v \; d_{\alpha \alpha} \widehat{B}^r_{\alpha}
232 \end{equation}
233 see equation~\ref{ref:MAG:EQU:201}, and the adjoint function $Y^*[\psi]$ for $Y_{\alpha}[\psi]$ is given from
234 \begin{equation}\label{ref:MAG:EQU:303}
235 \int_{\widehat{\Omega}} v \; d_{\alpha \alpha}^2 q_{,\alpha} Y^*[\psi]_{,\alpha } \;d\widehat{x} =
236 \int_{\widehat{\Omega}} r_{,{\alpha}} ,Y_{\alpha}[\psi] \;d\widehat{x}
237 \end{equation}
238 and finally
239 \begin{equation}\label{ref:MAG:EQU:201a}
240 \frac{\partial J^{mag}}{\partial k} = Y^*[\psi]_{,{\alpha}} B^r_{\alpha} - Y_i[\psi] B^b_{\alpha}
241 \end{equation}
242
243

  ViewVC Help
Powered by ViewVC 1.1.26