# Contents of /trunk/doc/inversion/CookMagnetic.tex

Revision 4433 - (show annotations)
Fri May 31 12:09:58 2013 UTC (6 years, 2 months ago) by gross
File MIME type: application/x-tex
File size: 8601 byte(s)
some clarifications on geodetic coordinates.
order of background magnetic flux density component has been corrected: input is now B_east, B_north, B_vertical.


 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 = [2201.*U.Nano*U.Tesla, 31232.*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 = [2201.*U.Nano*U.Tesla, 31232.*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_E, B_N, 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