/[escript]/branches/doubleplusgood/doc/inversion/CookMagnetic.tex
ViewVC logotype

Contents of /branches/doubleplusgood/doc/inversion/CookMagnetic.tex

Parent Directory Parent Directory | Revision Log Revision Log


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: 8599 byte(s)
Spelling fixes
1 \chapter{Magnetic Inversion}\label{Chp:cook:magnetic inversion}
2
3 \begin{figure}
4 \centering
5 \includegraphics[width=0.7\textwidth]{QLDWestMagneticDataPlot.png}
6 \caption{Magnetic anomaly data in $nT$ from Western Queensland, Australia
7 (file \examplefile{data/QLDWestMagnetic.nc}). Data obtained from Geoscience Australia.}
8 \label{FIG:P1:MAG:0}
9 \end{figure}
10
11 Magnetic data report the observed magnetic flux density over a region above the
12 surface of the Earth.
13 Similar to the gravity case the data are given as deviation from an expected
14 background magnetic flux density $B^b$ of the Earth.
15 Example data in units of $nT$ (nano Tesla) are shown in Figure~\ref{FIG:P1:MAG:0}.
16 It is the task of the inversion to recover the susceptibility distribution $k$
17 from the magnetic data collected. The approach for inverting magnetic data is
18 almost identical to the one used for gravity data.
19 In fact the \downunder script~\ref{code: magnetic1} used for the magnetic
20 inversion is very similar to the script~\ref{code: gravity1} for gravity inversion.
21
22 \begin{pyc}\label{code: magnetic1}
23 \
24 \begin{python}
25 # Header:
26 from esys.downunder import *
27 from esys.weipa import *
28 from esys.escript import unitsSI as U
29
30
31 # Step 1: set up domain
32 dom=DomainBuilder()
33 dom.setVerticalExtents(depth=40.*U.km, air_layer=6.*U.km, num_cells=25)
34 dom.setFractionalPadding(pad_x=0.2, pad_y=0.2)
35 B_b = [31232.*U.Nano*U.Tesla, 2201.*U.Nano*U.Tesla, -41405.*U.Nano*U.Tesla]
36 dom.setBackgroundMagneticFluxDensity(B_b)
37 dom.fixSusceptibilityBelow(depth=40.*U.km)
38
39 # Step 2: read magnetic data
40 source0=NetCdfData(NetCdfData.MAGNETIC, 'MagneticSmall.nc', scale_factor=U.Nano * U.Tesla)
41 dom.addSource(source0)
42
43 # Step 3: set up inversion
44 inv=MagneticInversion()
45 inv.setSolverTolerance(1e-4)
46 inv.setSolverMaxIterations(50)
47 inv.fixMagneticPotentialAtBottom(False)
48 inv.setup(dom)
49
50 # Step 4: run inversion
51 inv.getCostFunction().setTradeOffFactorsModels(0.1)
52 k = inv.run()
53
54 # Step 5: write reconstructed susceptibility to file
55 saveVTK("result.vtu", susceptibility=k)
56 \end{python}
57 \end{pyc}
58
59 \begin{figure}
60 \centering
61 \includegraphics[width=0.7\textwidth]{QLDMagContourMu01.png}
62 \caption{Contour plot of the susceptibility from a three-dimensional magnetic inversion (with $\mu=0.1$).
63 Colours represent values of susceptibility where high values are represented by
64 blue and low values are represented by red.}
65 \label{FIG:P1:MAG:1}
66 \end{figure}
67
68 The structure of the script is identical to the gravity case.
69 Following the header section importing the necessary modules the domain of the
70 inversion is defined in step one.
71 In step two the data are read and added to the domain builder.
72 Step three sets up the inversion and step four runs it.
73 Finally in step five the result is written to the result file, here
74 \file{result.vtu} in the \VTK format.
75 Results are shown in Figure~\ref{FIG:P1:MAG:1}.
76
77 Although scripts for magnetic and gravity inversion are largely identical there
78 are a few small differences which we are going to highlight now.
79 The magnetic inversion requires data about the background magnetic flux density
80 over the region of interest which is added to the domain by the statements
81 \begin{verbatim}
82 B_b = [31232.*U.Nano*U.Tesla, 2201.*U.Nano*U.Tesla, -41405.*U.Nano*U.Tesla]
83 dom.setBackgroundMagneticFluxDensity(B_b)
84 \end{verbatim}
85 Here it is assumed that the background magnetic flux density is constant across
86 the domain and is given as the list
87 \begin{verbatim}
88 B_b= [ B_N, B_E, B_V ]
89 \end{verbatim}
90 in units of Tesla (T) where
91 \member{B_N}, \member{B_E} and \member{B_V} refer to the north, east and
92 vertical component of the magnetic flux density, respectively.
93 Values for the magnetic flux density can be obtained by the International
94 Geomagnetic Reference Field (IGRF)~\cite{IGRF} (or the Australian Geomagnetic
95 Reference Field (AGRF)~\cite{AGRF} via \url{http://www.ga.gov.au/oracle/geomag/agrfform.jsp}).
96 Similar to the gravity case susceptibility below a certain depth can be set to
97 zero via the statement
98 \begin{verbatim}
99 dom.fixSusceptibilityBelow(depth=40.*U.km)
100 \end{verbatim}
101 where here the susceptibility below $40km$ is prescribed (this has no effect as
102 the depth of the domain is $40km$)\footnote{Notice that the method called is
103 different from the one in the case of gravity inversion.}.
104
105 Magnetic data are read and added to the domain with the following statements:
106 \begin{verbatim}
107 source0=NetCdfData(NetCdfData.MAGNETIC, 'MagneticSmall.nc', \
108 scale_factor=U.Nano * U.Tesla)
109 dom.addSource(source0)
110 \end{verbatim}
111 The first argument \member{NetCdfData.MAGNETIC} identifies the data read from
112 file \file{MagneticSmall.nc} (second argument) as magnetic data.The argument
113 \file{scale_factor} specifies the units (here $nT$) of the magnetic flux
114 density data in the file.
115 If scalar data are given it is assumed that the magnetic flux density anomalies
116 are measured in direction of the background magnetic flux density\footnote{The
117 default for \file{scale_factor} for magnetic data is $nT$.}.
118
119 Finally the inversion is created and run:
120 \begin{verbatim}
121 inv=MagneticInversion()
122 inv.fixMagneticPotentialAtBottom(False)
123 k = inv.run()
124 \end{verbatim}
125 The result for the susceptibility is named \member{k}. In this case the magnetic potential is
126 not fixed at the bottom of the domain. The magnetic potential is still set zero at the top of the domain.
127
128 We then write the result
129 to a \VTK file using
130 \begin{verbatim}
131 saveVTK("result.vtu", susceptibility=k)
132 \end{verbatim}
133 where the result of the inversion is tagged with the name \member{susceptibility}
134 as an identifier for the visualization software.
135
136 \begin{figure}
137 \begin{center}
138 \subfigure[$\mu=0.001$]{%
139 \label{FIG:P1:MAG:10 MU0001}
140 \includegraphics[width=0.45\textwidth]{QLDMagContourMu0001.png}
141 }%
142 \subfigure[$\mu=0.01$]{%
143 \label{FIG:P1:MAG:10 MU001}
144 \includegraphics[width=0.45\textwidth]{QLDMagContourMu001.png}
145 }\\ % ------- End of the first row ----------------------%
146 \subfigure[$\mu=0.1$]{%
147 \label{FIG:P1:MAG:10 MU01}
148 \includegraphics[width=0.45\textwidth]{QLDMagContourMu01.png}
149 }%
150 \subfigure[$\mu=1.$]{%
151 \label{FIG:P1:MAG:10 MU1}
152 \includegraphics[width=0.45\textwidth]{QLDMagContourMu1.png}
153 }\\ % ------- End of the second row ----------------------%
154 \subfigure[$\mu=10.$]{%
155 \label{FIG:P1:MAG:10 MU10}
156 \includegraphics[width=0.45\textwidth]{QLDMagContourMu10.png}
157 }%
158 \end{center}
159 \caption{3-D contour plots of magnetic inversion results with data from
160 Figure~\ref{FIG:P1:MAG:0} for various values of the model trade-off
161 factor $\mu$. Visualization has been performed in \VisIt.}
162 \label{FIG:P1:MAG:10}
163 \end{figure}
164
165 \begin{figure}
166 \begin{center}
167 \subfigure[$\mu=0.001$]{%
168 \label{FIG:P1:MAG:11 MU0001}
169 \includegraphics[width=0.45\textwidth]{QLDMagDepthMu0001.png}
170 }%
171 \subfigure[$\mu=0.01$]{%
172 \label{FIG:P1:MAG:11 MU001}
173 \includegraphics[width=0.45\textwidth]{QLDMagDepthMu001.png}
174 }\\ % ------- End of the first row ----------------------%
175 \subfigure[$\mu=0.1$]{%
176 \label{FIG:P1:MAG:11 MU01}
177 \includegraphics[width=0.45\textwidth]{QLDMagDepthMu01.png}
178 }%
179 \subfigure[$\mu=1.$]{%
180 \label{FIG:P1:MAG:11 MU1}
181 \includegraphics[width=0.45\textwidth]{QLDMagDepthMu1.png}
182 }\\ % ------- End of the second row ----------------------%
183 \subfigure[$\mu=10.$]{%
184 \label{FIG:P1:MAG:11 MU10}
185 \includegraphics[width=0.45\textwidth]{QLDMagDepthMu10.png}
186 }%
187 \end{center}
188 \caption{3-D slice plots of magnetic inversion results with data from
189 Figure~\ref{FIG:P1:MAG:0} for various values of the model trade-off
190 factor $\mu$. Visualization has been performed \VisIt.}
191 \label{FIG:P1:MAG:11}
192 \end{figure}
193
194 Figures~\ref{FIG:P1:MAG:10} and~\ref{FIG:P1:MAG:11} show results from the
195 inversion of the magnetic data shown in Figure~\ref{FIG:P1:MAG:0}.
196 In Figure~\ref{FIG:P1:MAG:10} surface contours are used to represent the
197 susceptibility while Figure~\ref{FIG:P1:MAG:11} uses contour lines
198 on a lateral plane intercept and two vertical plane intercepts.
199 The images show the strong impact of the trade-off factor $\mu$ on the result.
200 Larger values give more emphasis to the misfit term in the cost function
201 leading to rougher susceptibility distributions.
202 The result for $\mu=0.1$ seems to be the most realistic.
203

  ViewVC Help
Powered by ViewVC 1.1.26