1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% Copyright (c) 20032013 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/osl3.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=1e8}, 
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)}_jB_{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 
