1 
gross 
4055 
\chapter{Gravity Inversion}\label{Chp:cook:gravity inversion} 
2 



3 
gross 
4185 
\section{Introduction} 
4 
azadeh 
4151 

5 
caltinay 
4205 
In this part of the documentation we give an introduction on how to use the 
6 


\downunder module and the inversion driver functions to perform the inversion 
7 


of gravity and magnetic data. 
8 


The driver functions enable geologists and geophysicists to apply the 
9 


\downunder module quickly and in an easy way without requiring detailed 
10 


knowledge of the theory behind inversion or programming skills. 
11 


However, users who are interested in specializing or extending the inversion 
12 
caltinay 
4266 
capacity are referred to Part~\ref{part2} of this manual. 
13 
caltinay 
4205 
It is beyond the intention of this manual to give a detailed introduction to 
14 


geophysical inversion, in particular to the appropriate preprocessing of data 
15 
gross 
4268 
sets. 
16 
azadeh 
4156 

17 
caltinay 
4205 
The \downunder module described here is designed to calculate estimations for 
18 


the 3D distribution of density and/or susceptibility from 2D gravity and 
19 


magnetic data measured in ground or airborne surveys. 
20 


This process is generally called inversion of geophysical data. 
21 


Following the standard assumption it is assumed that the data are measured as 
22 


perturbation of an expected gravity and/or magnetic field of the Earth. 
23 


In this context measured gravity and magnetic data are in fact describing 
24 


anomalies in the gravity and magnetic field. 
25 


As a consequence the inversion process provides corrections to an average 
26 


density (typically $2670 kg/m^3$) and susceptibility (typically $0$). 
27 


So in the following we will always assume that given data are anomalies and 
28 


therefore not in all cases explicitly use the terms gravity anomalies or 
29 


density corrections but just use the terms gravity and density. 
30 
azadeh 
4127 

31 
caltinay 
4205 
In this chapter we will give a detailed introduction into usage of the driver 
32 


functions for inversion of gravity data. 
33 


In the following chapters~\ref{Chp:cook:magnetic inversion} and~\ref{Chp:cook:joint inversion} 
34 


we will discuss the inversion of magnetic data, and the joint inversion of 
35 


gravity and magnetic data using \downunder. 
36 


Note, that the principles introduced in this chapter apply to magnetic and 
37 


joint inversion so the presentation for these problem classes is kept short 
38 


and users interested in magnetic data only should still work through this 
39 


chapter on gravity data. 
40 
azadeh 
4161 

41 
caltinay 
4205 
To run the examples discussed you need to have \escript (version 3.3.1 or 
42 


newer) installed on your computer. 
43 
caltinay 
4266 
Moreover, if you want to visualize the results you need to have access to a 
44 


data plotting software which is able to process \VTK input files, e.g. 
45 


\mayavi or \VisIt. 
46 


As \mayavi can be easily obtained and installed for most platforms the 
47 


tutorial includes commands to visualize output files using \mayavi . 
48 


However, it is pointed out that \VisIt is the preferred visualization tool for 
49 


\escript as it can deal with very large data sets more efficiently. 
50 
gross 
4185 

51 
gross 
4191 
\begin{figure} 
52 


\centering 
53 
gross 
4199 
\includegraphics[width=0.7\textwidth]{QLDWestGravityDataPlot.png} 
54 
caltinay 
4205 
\caption{Gravity Anomaly Data in $mgal$ from Western Queensland, Australia 
55 
caltinay 
4269 
(file \examplefile{data/QLDWestGravity.nc}). Data obtained from Geoscience Australia.} 
56 
caltinay 
4266 
%\AZADEH{Tracey R, Bacchin M, \& Wynne P. 2008. In preparation. AAGD07: A new absolute gravity datum for Australian gravity and new standards for the Australian National Gravity Database. Exploration Geophysics.}} 
57 
gross 
4191 
\label{FIG:P1:GRAV:0} 
58 


\end{figure} 
59 
azadeh 
4161 

60 
azadeh 
4127 
\begin{figure} 
61 
caltinay 
4205 
\centering% 
62 
caltinay 
4235 
\includegraphics[width=0.9\textwidth]{QLDGravContourMu10.png} 
63 
caltinay 
4266 
\caption{3D contour plot of the density distribution obtained by inversion of 
64 
caltinay 
4269 
file \file{data/QLDWestGravity.nc} (with $\mu=10$). 
65 
caltinay 
4205 
Colours represent values of density where high values are represented by 
66 
caltinay 
4235 
blue and low values are represented by red.} 
67 
gross 
4185 
\label{FIG:P1:GRAV:1} 
68 
azadeh 
4127 
\end{figure} 
69 



70 
gross 
4185 
\section{How does it work?} 
71 
caltinay 
4205 
The execution of the inversion is controlled by a script which, in essence, is 
72 


a text file and can be edited using any text editor. 
73 


The script contains a series of statements which are executed by an 
74 


interpreter which is an executable program reading the text file and executing 
75 


the statements linebyline. 
76 


In the case of \downunder the interpreter is \Python. 
77 


In order to be able to process the statements in each line of the script 
78 


certain rules (called syntax) need to be obeyed. 
79 


There is a large number of online tutorials for \Python available\footnote{e.g. 
80 


\url{http://www.tutorialspoint.com/python} and \url{http://doc.pyschools.com}}. 
81 


We also refer to the \escript cook book \cite{ESCRIPTCOOKBOOK} and user's 
82 


guide \cite{ESCRIPT} which is in particular useful for users who like to dive 
83 


deeper into \downunder. 
84 


For this part of the manual no \Python knowledge is required but it is 
85 


recommended that users acquire some basic knowledge on \Python as they 
86 


progress in their work with \downunder. 
87 
gross 
4185 

88 
caltinay 
4205 
The following script~\ref{code: gravity1}\footnote{The script is similar to 
89 


\examplefile{grav_netcdf.py} within the \escript example file directory.} is a 
90 


simple example to run an inversion for gravity data: 
91 



92 


\begin{pyc}\label{code: gravity1} 
93 
gross 
4185 
\ 
94 
caltinay 
4205 
\begin{python} 
95 
gross 
4185 
# Header: 
96 


from esys.downunder import * 
97 


from esys.weipa import * 
98 
caltinay 
4205 
from esys.escript import unitsSI as U 
99 
gross 
4185 

100 


# Step 1: set up domain 
101 


dom=DomainBuilder() 
102 
caltinay 
4205 
dom.setVerticalExtents(depth=40.*U.km, air_layer=6.*U.km, num_cells=25) 
103 
gross 
4185 
dom.setFractionalPadding(pad_x=0.2, pad_y=0.2) 
104 
caltinay 
4205 
dom.fixDensityBelow(depth=40.*U.km) 
105 
gross 
4185 

106 


# Step 2: read gravity data 
107 
caltinay 
4269 
source0=NetCdfData(NetCdfData.GRAVITY, 'GravitySmall.nc') 
108 
gross 
4185 
dom.addSource(source0) 
109 



110 


# Step 3: set up inversion 
111 


inv=GravityInversion() 
112 


inv.setSolverTolerance(1e4) 
113 


inv.setSolverMaxIterations(50) 
114 


inv.setup(dom) 
115 



116 


# Step 4: run inversion 
117 


inv.getCostFunction().setTradeOffFactorsModels(10.) 
118 


rho = inv.run() 
119 



120 


# Step 5: write reconstructed density to file 
121 
caltinay 
4205 
saveVTK("result.vtu", density=rho) 
122 


\end{python} 
123 
gross 
4185 
\end{pyc} 
124 
caltinay 
4205 
The result, in this case the density distribution, is written to an external 
125 


file for further processing. You can copy and paste the text of the script 
126 


into a file of any name, let's say for further reference we use the file 
127 


name \file{grav.py}. 
128 


It is recommended to use the extension \file{.py} to identify the file as a 
129 


\Python script. 
130 


We will discuss the statements of the script later in this chapter. 
131 
gross 
4185 

132 
caltinay 
4269 
Obviously the inversion needs to be fed with some gravity data. You can find 
133 


example data from western Queensland, Australia in two resolutions in the 
134 


\escript example directory. In this case the data are loaded from the file 
135 


\examplefile{GravitySmall.nc} which is given in the \netcdf file format. 
136 
caltinay 
4205 
After you have copied this file into the directory in which you have saved the 
137 


script \file{grav.py} you can run the program using the command line 
138 
gross 
4185 
\begin{verbatim} 
139 


runescript grav.py 
140 


\end{verbatim} 
141 
caltinay 
4205 
We are running \file{grav.py} through the \escript startup command since 
142 


\escript is used as a back end for the inversion algorithm\footnote{Please see 
143 


the \escript user's guide~\cite{ESCRIPT} on how to run your script in parallel 
144 


using threading and/or MPI.}. 
145 


Obviously it is assumed that you have an installation of \escript available on 
146 


your computer, see \url{https://launchpad.net/escriptfinley}. 
147 
gross 
4185 

148 
caltinay 
4205 
After the execution has successfully completed you will find the result file 
149 


\file{result.vtu} in the directory from where you have started the execution 
150 


of the script. 
151 


The file has the \VTK format and can be imported easily into many 
152 


visualization tools. 
153 


One option is the \mayavi package which is available on most platforms. 
154 
caltinay 
4285 
You can invoke the visualization using the commands 
155 
gross 
4185 
\begin{verbatim} 
156 


mayavi2 d result.vtu m SurfaceMap 
157 


\end{verbatim} 
158 
caltinay 
4205 
from the command line. 
159 


Figure~\ref{FIG:P1:GRAV:1} shows the result of this inversion as a contour 
160 
caltinay 
4269 
plot\footnote{These plots were generated by \VisIt using the higher resolution 
161 


data.}, while the gravity anomaly data is shown in Figure~\ref{FIG:P1:GRAV:0}. 
162 
caltinay 
4205 
We will discuss data later in Section~\ref{SEC:P1:GRAV:DATA}. 
163 
gross 
4185 

164 
caltinay 
4269 
Let us take a closer look at the script\footnote{In \Python lines starting 
165 


with `\#` are comments and are not processed further.}. Besides the header 
166 
caltinay 
4205 
section one can separate the script into five steps: 
167 
gross 
4185 
\begin{enumerate} 
168 
caltinay 
4205 
\item set up domain on which the inversion problem is solved 
169 


\item load the data 
170 


\item setup the inversion problem 
171 


\item run the inversion 
172 


\item further processing of the result. Here we write the reconstructed 
173 


density distribution to a file. 
174 
gross 
4185 
\end{enumerate} 
175 
caltinay 
4205 
In the following we will discuss the steps of the scripts in more detail. 
176 


Before we do this it is pointed out that the header section, following 
177 


\Python conventions, makes all required packages available to access within 
178 


the script. 
179 


At this point we will not discuss this in more details but emphasize that the 
180 


header section is a vital part of the script. 
181 


It is is required in each \downunder inversion script and should not be 
182 


altered except if additional modules are needed. 
183 
gross 
4185 

184 
azadeh 
4156 
\begin{figure} 
185 


\centering 
186 
caltinay 
4205 
\includegraphics[width=\textwidth]{dom2D.pdf} 
187 


\caption{2D domain setup for gravity inversion} 
188 
gross 
4185 
\label{FIG:P1:GRAV:2} 
189 


\end{figure} 
190 



191 
caltinay 
4205 
\section{Creating the Inversion Domain} 
192 
caltinay 
4269 
First step in Script~\ref{code: gravity1} is the definition of the domain over 
193 
caltinay 
4205 
which the inversion is performed. 
194 


We think of the domain as a block with orthogonal, plain faces. 
195 


Figure~\ref{FIG:P1:GRAV:2} shows the setup for a twodimensional domain 
196 


(see also Figure~\ref{fig:domainBuilder} for 3D). 
197 


The lateral coordinates along the surface of the Earth are denoted by $x$ and 
198 
caltinay 
4266 
$y$ (only $x$direction is used in 2D). 
199 
caltinay 
4205 
The $z$ direction defines the vertical direction where the part above the 
200 


surface has positive coordinates and the subsurface negative coordinates. 
201 


The height of the section above the surface, which is assumed to be filled 
202 


with air, needs to be set by the user. 
203 


The inversion assumes that the density in the section is known to be 
204 


zero\footnote{Always keeping in mind that these are not absolute values but 
205 


anomalies.}. 
206 


The density below the surface is unknown and is calculated through the 
207 


inversion. The user needs to specify the depth below the surface in which the 
208 


density is to be calculated. 
209 


The lateral extension of the domain is defined by the data sets fed into the 
210 


inversion. 
211 


It is chosen large enough to cover all data sets (in case more than one is 
212 


used). In order to reduce the impact of the boundary a padding zone around the 
213 


data sets can be introduced. 
214 



215 
gross 
4185 
\begin{figure} 
216 


\centering 
217 
caltinay 
4205 
\includegraphics[width=\textwidth]{dom2DFEM.pdf} 
218 


\caption{Cell distribution and boundary conditions for a 2D domain} 
219 
gross 
4191 
\label{FIG:P1:GRAV:3} 
220 
gross 
4185 
\end{figure} 
221 



222 
caltinay 
4205 
The reconstruction of the gravity field from an estimated density distribution 
223 


is the key component of the inversion process. 
224 


To do this \downunder uses the finite element method (FEM). 
225 


We need to introduce a grid over the domain, see Figure~\ref{FIG:P1:GRAV:3}. 
226 


The number of vertical cells is set by the user while the number of horizontal 
227 


cells is derived from the grid spacing of the gravity data set(s). 
228 


It is assumed that gravity field data given are constant across a cell. 
229 


To be able to reconstruct the gravity field some assumptions on the values of 
230 


the gravity field on the domain boundary have to be made. 
231 


\downunder assumes that on all faces the lateral gravity field component 
232 


equals zero. No assumptions on the horizontal components are 
233 


made\footnote{It is assumed that the gravity potential equals zero on the top 
234 


and bottom surface, see Section~\ref{sec:forward gravity} for details}% 
235 


\footnote{Most inversion codes use Green's functions over an unbounded domain 
236 


to reconstruct the gravity field. This approach makes the assumption that the 
237 


gravity field (or potential) is converging to zero when moving away from the 
238 


region of interest. The boundary conditions used here are stronger in the 
239 


sense that the lateral gravity component is enforced to be zero in a defined 
240 


distance of the region of interest but weaker in the sense that no constraint 
241 


on the horizontal component is applied.}. 
242 
gross 
4185 

243 
caltinay 
4205 
In script~\ref{code: gravity1} the statement 
244 
gross 
4185 
\begin{verbatim} 
245 


dom=DomainBuilder() 
246 


\end{verbatim} 
247 
caltinay 
4205 
creates something like a factory to build a domain. 
248 


We then define the features of the domain we would like to create: 
249 
gross 
4185 
\begin{verbatim} 
250 
caltinay 
4205 
dom.setVerticalExtents(depth=40.*U.km, air_layer=6.*U.km, num_cells=25) 
251 
gross 
4185 
dom.setFractionalPadding(pad_x=0.2, pad_y=0.2) 
252 


\end{verbatim} 
253 
caltinay 
4205 
Here we specify the depth of the domain to $40 km$, the thickness of the air 
254 


layer above the surface to $6 km$ and the number of vertical cells to $25$. 
255 


We also introduce a lateral padding of $20 \%$ of the expansion of the gravity 
256 


data on each side of the data and in both lateral directions. 
257 
gross 
4185 

258 
caltinay 
4205 
In some cases it can be appropriate to assume that the density below a certain 
259 


depth is zero\footnote{As we are in fact calculating density corrections this 
260 


means that the density is assumed to be the average density.}. 
261 


The statement 
262 
gross 
4185 
\begin{verbatim} 
263 
caltinay 
4205 
dom.fixDensityBelow(depth=40.*U.km) 
264 
gross 
4185 
\end{verbatim} 
265 
caltinay 
4205 
introduces this constraint. 
266 
caltinay 
4214 
As in the case discussed here if the depth for zero density is not less than 
267 


the depth of the domain no constraint at depth is applied to the density. 
268 
gross 
4185 

269 
caltinay 
4205 
\downunder uses the metrekilogramsecond based International System of Units 
270 
caltinay 
4269 
(SI)\footnote{see \url{http://en.wikipedia.org/wiki/International_System_of_Units}}\index{SI}. 
271 
caltinay 
4205 
So all values must be converted to appropriate units. 
272 
caltinay 
4214 
This does not apply to geographic coordinates which in \downunder are given in 
273 


fractional degrees (as a floating point number) to represent longitude and 
274 
caltinay 
4205 
latitude. In the script we have used the expression 
275 
gross 
4191 
\begin{verbatim} 
276 
caltinay 
4205 
depth=40.*U.km 
277 
gross 
4191 
\end{verbatim} 
278 
caltinay 
4205 
to define the depth of the domain to $40 km$. 
279 
caltinay 
4266 
The expression \verbU.km denotes the unit $km$ (kilometer) and ensures 
280 


appropriate conversion of the value $40$ into the base unit $m$ (meter). 
281 
caltinay 
4205 
It is recommended to add units to values (where present) in order to make sure 
282 


that the final values handed into \downunder is given with the appropriate 
283 


units. 
284 


The physical units module of \escript, which we have imported here under the 
285 


name \verbU in the script header, defines a large number of physical units 
286 


and constants, please see~\cite{ESCRIPT} and~\cite{ESCRIPTONLINE}. 
287 
gross 
4185 

288 
caltinay 
4205 
\section{Loading Gravitational Data}\label{SEC:P1:GRAV:DATA} 
289 


In practice gravity acceleration is measured in various ways, for instance by 
290 


airborne surveys~\cite{Telford1990a}. 
291 


\downunder assumes that all data supplied as input are already appropriately 
292 


preprocessed. In particular, corrections for 
293 
gross 
4191 
\begin{itemize} 
294 
caltinay 
4205 
\item freeair, to remove effects from altitude above ground; 
295 


\item latitude, to remove effects from ellipsoidicity of the Earth; 
296 


\item terrain, to remove effects from topography 
297 
gross 
4191 
\end{itemize} 
298 
caltinay 
4205 
must have been applied to the data. 
299 


In general, data prepared in such a form are called Bouguer anomalies~\cite{Telford1990a}. 
300 
gross 
4185 

301 
caltinay 
4205 
To load gravity data into \downunder the data are given on a plane parallel 
302 


to the surface of the Earth at a constant altitude, see 
303 


diagram~\ref{FIG:P1:GRAV:2}. 
304 


The data need to be defined over a rectangular grid covering a subregion of 
305 


the Earth surface. 
306 
caltinay 
4214 
The grid uses a geographic coordinate system with latitudes and longitudes 
307 


assumed to be given in the Clarke 1866 geodetic system. 
308 
caltinay 
4205 
Figure~\ref{FIG:P1:GRAV:0} shows an example of such a data set from the 
309 


western parts of Queensland, Australia. 
310 


The data set covers a rectangular region between $140^o$ and $141^o$ east 
311 


and between $20^o$ and $21^o$ south. 
312 


Notice that latitude varies between $90^o$ to $90^o$ where negative signs 
313 


refer to places in the southern hemisphere and longitude varies between 
314 
caltinay 
4269 
$180^o$ to $180^o$ where negative signs refer to places west of Greenwich. 
315 
caltinay 
4205 
The colour at a location represents the value of the vertical Bouguer gravity 
316 


anomaly at this point at the surface of the Earth. 
317 
caltinay 
4266 
Values in this data set range from $160 \; mgal$ to about $500 \; mgal$\footnote{The unit 
318 
caltinay 
4205 
$mgal$ means milli $gal$ (galileo) with $1 \; gal = 0.01 \frac{m}{sec^2}$.} 
319 


over a $121 \times 121$ grid. 
320 
gross 
4191 

321 
caltinay 
4271 
In general, a patch of gravity data needs to be defined over a plane 
322 
caltinay 
4205 
\verbNX $\times$ \verbNY where \verbNX and \verbNY define the number 
323 


of grid lines in the longitude (\verbX) and the latitude (\verbY) 
324 


direction, respectively. 
325 


The grid is spanned from an origin with spacing \verbDELTA_X and 
326 


\verbDELTA_Y in the longitude and the latitude direction, respectively. 
327 


The gravity data for all grid points need to be given as an \verbNX 
328 


$\times$ \verbNY array. 
329 
caltinay 
4266 
If available, measurement errors can be associated with the gravity data. 
330 
caltinay 
4205 
The values are given as an \verbNX $\times$ \verbNY array matching the 
331 


shape of the gravity array. 
332 
caltinay 
4266 
Note that data need not be available on every single point of the grid, see 
333 


Section~\ref{SEC:P1:GRAV:REMARK:DATAHOLES} for more information on this. 
334 



335 


Currently, two data file formats are supported, namely \emph{ER Mapper Raster} 
336 


\cite{ERMAPPER} files and \netcdf~\cite{NETCDF} files. 
337 


In the examples of this chapter we use \netcdf files and refer to 
338 


Section~\ref{SEC:P1:GRAV:REMARK:ERMAPPER} and 
339 


Section~\ref{sec:ref:DataSource:ERM} for more information on using ER Mapper 
340 


Raster files. 
341 


If you have data in any other format you have the option of writing a suitable 
342 


reader (for advanced users, see Chapter~\ref{Chp:ref:data sources}) or, 
343 


assuming you are able to read the data in \Python, refer to the example 
344 


script \examplefile{create_netcdf.py} which shows how to create a file 
345 
caltinay 
4205 
in the \netcdf file format~\cite{NETCDF} compatible with \downunder from 
346 


a data array. 
347 
gross 
4191 

348 
caltinay 
4205 
In script~\ref{code: gravity1} we use the statement 
349 
gross 
4199 
\begin{verbatim} 
350 
caltinay 
4269 
source0=NetCdfData(NetCdfData.GRAVITY, 'GravitySmall.nc') 
351 
gross 
4199 
\end{verbatim} 
352 
caltinay 
4269 
to load the gravity data stored in \examplefile{GravitySmall.nc} in the 
353 


\netcdf format. 
354 
caltinay 
4205 
Within the script the data set is now available under the name \verbsource0. 
355 


We need to link the data set to the \verbDomainBuilder using 
356 
gross 
4199 
\begin{verbatim} 
357 


dom.addSource(source0) 
358 


\end{verbatim} 
359 
caltinay 
4205 
At the time the domain for the inversion is built the \verbDomainBuilder 
360 


will use the information about origin, extent and spacing along with the other 
361 


options provided to build an appropriate domain. 
362 


As at this point a flat Earth is assumed geographic coordinates used to 
363 


represent data in the input file are mapped to a (local) Cartesian coordinate 
364 


system. This is achieved by projecting the geographic coordinates into the 
365 


\emph{Universal Transverse Mercator} (UTM) coordinate system\footnote{See e.g. 
366 


\url{http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system}}. 
367 



368 
caltinay 
4271 
There are a few optional arguments you can add when constructing a data source. 
369 


While Section~\ref{sec:ref:DataSource} has a detailed description of all 
370 


arguments it is worth noting a few. 
371 


Firstly, it is important to know that data slices are assumed to be at altitude 
372 


$0$ by default. This can be easily changed though: 
373 


\begin{verbatim} 
374 


source0=NetCdfData(NetCdfData.GRAVITY, 'GravitySmall.nc', altitude=2.5*U.km) 
375 


\end{verbatim} 
376 


Another important setting is the scale or unit of the measurements. 
377 


The default is dependent on the data type and for gravity anomalies a scale 
378 


of $\frac{\mu m}{sec^2}$ (or $0.1 \; mgal$) is assumed. 
379 


To change the default scale to $mgal$ (which is $10^{5} \frac{m}{sec^2}$), 
380 


you could use: 
381 


\begin{verbatim} 
382 


source0=NetCdfData(NetCdfData.GRAVITY, 'GravitySmall.nc', scale_factor=1e5) 
383 


\end{verbatim} 
384 


Finally, it is possible to specify measurement errors (i.e. uncertainties) 
385 


alongside the data. 
386 


Since these can never be zero, a value of $2$ units is used if nothing else 
387 


has been specified. 
388 


The error value is assumed to be given in the same units as the data so the 
389 


default value translates to an error of $0.2 \; mgal$. 
390 


There are two possibilities to specify the error, namely by providing a 
391 


constant value which is applied to all data points: 
392 


\begin{verbatim} 
393 


source0=NetCdfData(NetCdfData.GRAVITY, 'GravitySmall.nc', error=1.7) 
394 


\end{verbatim} 
395 


or, if the information is available in the same \netcdf file under the name 
396 


\verberrors, provide \downunder with the variable name: 
397 


\begin{verbatim} 
398 


source0=NetCdfData(NetCdfData.GRAVITY, 'GravitySmall.nc', error="errors") 
399 


\end{verbatim} 
400 



401 
caltinay 
4205 
It is important to keep an eye on the complexity of the inversion. 
402 


A good measure is the total number of cells being used. 
403 


Assume we have given a data set on a $20 \times 20$ grid and we add lateral 
404 


padding of, say, $20 \%$ to each side of the data, the lateral number of cells 
405 


becomes $(20\cdot 1.4)\times (20\cdot 1.4)=1.4^2\cdot 20^2\approx 2\cdot 10^2=800$. 
406 


If we use $20$ cells in the vertical direction we end up with a total number 
407 


of $800 \times 20 = 16,000$ cells. 
408 


This size can be easily handled by a modern desktop PC. 
409 


If we increase the grid size of the data to $40 \times 40$ points and use $40$ 
410 


cells in the vertical extent we get a total of $(2\cdot 40^2)\cdot 40=128,000$ 
411 


cells, a problem size which is considerably larger but can still be handled by 
412 


a desktop computer. 
413 


Taking this one step further, if the amount of data is increased to 
414 


$200\times 200$ points and we use $200$ cells in the vertical extent the 
415 


domain will contain $16,000,000$ ($16$ million) cells. 
416 


This scenario requires a computer with enough memory and (a) fast processor(s) 
417 


to run the inversion. 
418 


This estimate of complexity growth applies to the case where the increase 
419 


of data grid size is driven by an increase of resolution where it is 
420 


recommended to increase the vertical resolution in synch with the lateral 
421 


resolution. Note that if more than one data set is used the target resolution 
422 
caltinay 
4279 
will be the resolution of the finest data set (see also 
423 


Section~\ref{SEC:P1:GRAV:REMARK:MULTIDATA}). 
424 
caltinay 
4205 
In other cases the expansion of the region of interest drives an increase of 
425 


data grid size and the increase of total number of cells is less dramatic as 
426 


the vertical number of cells can remain constant while keeping a balanced 
427 


resolution in vertical and lateral direction. 
428 
gross 
4191 

429 
caltinay 
4205 
\section{Setting up the Inversion and Running it} 
430 


We are now at step three of script~\ref{code: gravity1} in which the actual 
431 


inversion is set up. 
432 


First we create an empty inversion under the name \verbinv: 
433 
gross 
4199 
\begin{verbatim} 
434 


inv=GravityInversion() 
435 


\end{verbatim} 
436 
caltinay 
4205 
As indicated by the name we can use \verbinv to perform an inversion of 
437 


gravity data\footnote{\verbGravityInversion is a driver with a simplified 
438 


interface which is provided for convenience. See Part~\ref{part2} for more 
439 


details on how to write inversion scripts with more general functionality, e.g. 
440 


constraints.}. The inversion is an iterative process which sequentially 
441 


calculates updates to the density distribution in an attempt to improve the 
442 


match of the gravity field produced by the density distribution with the data. 
443 


Termination of the iteration is controlled by the tolerance which is set by 
444 


the user: 
445 
gross 
4199 
\begin{verbatim} 
446 


inv.setSolverTolerance(1e4) 
447 


\end{verbatim} 
448 
caltinay 
4205 
Here we set the tolerance to $10^{4}$, i.e. the iteration is terminated if 
449 


the maximum density correction is less than or equal to $10^{4}$ relative to 
450 


the maximum value of estimated density anomaly. 
451 


In case the iteration does not converge a maximum number of iteration steps is 
452 


set: 
453 
gross 
4199 
\begin{verbatim} 
454 


inv.setSolverMaxIterations(50) 
455 


\end{verbatim} 
456 
caltinay 
4205 
If the maximum number of iteration steps (here $50$) is reached the iteration 
457 


process is aborted and an error message is printed. 
458 


In this case you can try to rerun the inversion with a larger value for the 
459 


maximum number of iteration steps. 
460 


If even for a very large number of iteration steps no convergence is achieved, 
461 
caltinay 
4269 
it is very likely that the inversion has not been set up properly. 
462 
gross 
4191 

463 
gross 
4199 
The statement 
464 


\begin{verbatim} 
465 


inv.setup(dom) 
466 


\end{verbatim} 
467 
caltinay 
4205 
links the inversion with the domain and the data. 
468 


At this step  as we are solving a gravity inversion problem  only 
469 


gravitational data attached to the domain builder \verbdom are considered. 
470 


Internally a cost function $J$ is created which is minimized during the 
471 


inversion iteration. 
472 


It is a combination of a measure of the data misfit of the gravity field from 
473 


the given density distribution and a measure of the smoothness of the density 
474 


distribution. 
475 


The latter is often called the regularization term. 
476 


By default the gradient of density is used as the regularization term, see 
477 


also Section~\ref{SEC:P1:GRAV:REMARK:REG}. 
478 


Obviously, the result of the inversion is sensitive to the weighting between 
479 


the misfit and the regularization. 
480 


This tradeoff factor $\mu$ for the misfit function is set by the following 
481 


statement: 
482 
gross 
4199 
\begin{verbatim} 
483 


inv.getCostFunction().setTradeOffFactorsModels(0.1) 
484 


\end{verbatim} 
485 
caltinay 
4205 
Here we set $\mu=0.1$. The statement \verbinv.setup must appear in the 
486 


script before setting the tradeoff factor. 
487 


A small value for the tradeoff factor $\mu$ will give more emphasis to the 
488 


regularization component and create a smoother density distribution. 
489 


A large value of the tradeoff factor $\mu$ will emphasize the misfit more 
490 


and typically creates a better fit to the data and a rougher density 
491 


distribution. 
492 


It is important to keep in mind that the regularization reduces noise in the 
493 


date and in fact gives the problem a unique solution. 
494 


Consequently, the tradeoff factor $\mu$ may not be chosen too large in order 
495 


control the noise on the solution and ensure convergence in the iteration 
496 


process. 
497 
gross 
4185 

498 
gross 
4199 
We can now actually run the inversion: 
499 


\begin{verbatim} 
500 


rho = inv.run() 
501 


\end{verbatim} 
502 
caltinay 
4205 
The answer as calculated during the inversion is returned and can be accessed 
503 


under the name \verbrho. 
504 


As pointed out earlier the iteration process may fail in which case the 
505 


execution of the script is aborted with an error message. 
506 
gross 
4199 

507 
caltinay 
4205 
\section{Taking a Look} 
508 


In the final step of script~\ref{code: gravity1} the calculated density 
509 


distribution is written to an external file. 
510 


A popular file format used by several visualization packages such as 
511 


\VisIt~\cite{VISIT} and \mayavi~\cite{MAYAVI} is the \VTK file format. 
512 


The result of the inversion which has been named \verbrho can be written to 
513 


the file \file{result.vtu} by adding the statement 
514 
gross 
4203 
\begin{verbatim} 
515 


saveVTK("result.vtu", density=rho) 
516 


\end{verbatim} 
517 
caltinay 
4205 
at the end of script. 
518 


The inversion solution is tagged with the name \verbdensity in the result 
519 


file, however any other name for the tag could be used. 
520 


As the format is textbased (as opposed to binary) \VTK files tend to be very 
521 


large and take compute time to create, in particular when it comes to large 
522 


numbers of cells ($>10^6$). 
523 


For large problems it is more efficient to use the \SILO file format~\cite{SILO}. 
524 


\SILO files tend to be smaller and are faster generated and read. 
525 


It is the preferred format to import results into the visualization program 
526 


\VisIt~\cite{VISIT} which is particularly suited for the visualization of 
527 


large data sets. 
528 


Inversion results can directly exported into \SILO files using the statement 
529 
gross 
4203 
\begin{verbatim} 
530 


saveSilo("result.silo", density=rho) 
531 


\end{verbatim} 
532 
caltinay 
4205 
replacing the \verbsaveVTK(...) statement. 
533 


Similar to \VTK files the result \verbrho is tagged with the name 
534 


\verbdensity so it can be identified in the visualization program. 
535 
gross 
4185 

536 
gross 
4203 
\begin{figure} 
537 
caltinay 
4205 
\begin{center} 
538 
gross 
4203 
\subfigure[$\mu=0.1$]{% 
539 
caltinay 
4205 
\label{FIG:P1:GRAV:10 MU01} 
540 
azadeh 
4229 
\includegraphics[width=0.45\textwidth]{QLDGravContourMu01.png} 
541 
caltinay 
4205 
}% 
542 


\subfigure[$\mu=1.$]{% 
543 


\label{FIG:P1:GRAV:10 MU1} 
544 
azadeh 
4229 
\includegraphics[width=0.45\textwidth]{QLDGravContourMu1.png} 
545 
caltinay 
4205 
}\\ %  End of the first row % 
546 


\subfigure[$\mu=10.$]{% 
547 


\label{FIG:P1:GRAV:10 MU10} 
548 
azadeh 
4229 
\includegraphics[width=0.45\textwidth]{QLDGravContourMu10.png} 
549 
caltinay 
4205 
}% 
550 


\subfigure[$\mu=100.$]{% 
551 


\label{FIG:P1:GRAV:10 MU100} 
552 
azadeh 
4229 
\includegraphics[width=0.45\textwidth]{QLDGravContourMu100.png} 
553 


}\\ %  End of the second row % 
554 


\subfigure[$\mu=1000.$]{% 
555 


\label{FIG:P1:GRAV:10 MU1000} 
556 


\includegraphics[width=0.45\textwidth]{QLDGravContourMu1000.png} 
557 
caltinay 
4205 
}% 
558 


\end{center} 
559 


\caption{3D contour plots of gravity inversion results with data from 
560 


Figure~\ref{FIG:P1:GRAV:0} for various values of the model tradeoff 
561 


factor $\mu$. Visualization has been performed in \VisIt. 
562 


\AZADEH{check images.}} 
563 


\label{FIG:P1:GRAV:10} 
564 
gross 
4203 
\end{figure} 
565 



566 
gross 
4204 
\begin{figure} 
567 
caltinay 
4205 
\begin{center} 
568 
gross 
4204 
\subfigure[$\mu=0.1$]{% 
569 
caltinay 
4205 
\label{FIG:P1:GRAV:11 MU01} 
570 
azadeh 
4229 
\includegraphics[width=0.45\textwidth]{QLDGravDepthMu01.png} 
571 
caltinay 
4205 
}% 
572 


\subfigure[$\mu=1.$]{% 
573 


\label{FIG:P1:GRAV:11 MU1} 
574 
azadeh 
4229 
\includegraphics[width=0.45\textwidth]{QLDGravDepthMu1.png} 
575 
caltinay 
4205 
}\\ %  End of the first row % 
576 


\subfigure[$\mu=10.$]{% 
577 


\label{FIG:P1:GRAV:11 MU10} 
578 
azadeh 
4229 
\includegraphics[width=0.45\textwidth]{QLDGravDepthMu10.png} 
579 
caltinay 
4205 
}% 
580 


\subfigure[$\mu=100.$]{% 
581 


\label{FIG:P1:GRAV:11 MU100} 
582 
azadeh 
4229 
\includegraphics[width=0.45\textwidth]{QLDGravDepthMu100.png} 
583 


}\\ %  End of the second row % 
584 


\subfigure[$\mu=1000.$]{% 
585 


\label{FIG:P1:GRAV:11 MU1000} 
586 


\includegraphics[width=0.45\textwidth]{QLDGravDepthMu1000.png} 
587 
caltinay 
4205 
}% 
588 


\end{center} 
589 


\caption{3D slice plots of gravity inversion results with data from 
590 


Figure~\ref{FIG:P1:GRAV:0} for various values of the model tradeoff 
591 


factor $\mu$. Visualization has been performed \VisIt. 
592 


\AZADEH{check images.}} 
593 


\label{FIG:P1:GRAV:11} 
594 
gross 
4204 
\end{figure} 
595 



596 
caltinay 
4205 
Figures~\ref{FIG:P1:GRAV:10} and~\ref{FIG:P1:GRAV:11} show two different 
597 


styles of visualization generated in \VisIt using the result of the inversion 
598 


of the gravity anomalies shown in Figure~\ref{FIG:P1:GRAV:0}. 
599 


The inversions have been performed with different values for the model 
600 


tradeoff factor $\mu$. 
601 


The visualization shows clearly the smoothing effect of lower values for the 
602 


tradeoff factors. 
603 


For larger values of the tradeoff factor the density distribution becomes 
604 


rougher showing larger details. 
605 


Computational costs are significantly higher for larger tradeoff factors. 
606 
caltinay 
4269 
Moreover, noise in the data has a higher impact on the result. 
607 
caltinay 
4205 
Typically several runs are required to adjust the value for the tradeoff 
608 


factor to the datasets used. 
609 
gross 
4204 

610 
caltinay 
4269 
For some analysis tools it is useful to process the results in form of 
611 
caltinay 
4271 
Commaseparated Values (\CSV)\footnote{see 
612 


\url{http://en.wikipedia.org/wiki/Commaseparated_values}}. 
613 
caltinay 
4205 
Such a file can be created using the statement 
614 
gross 
4203 
\begin{verbatim} 
615 
gross 
4234 
saveDataCSV("result.csv", x=rho.getFunctionSpace().getX(), density=rho) 
616 
gross 
4203 
\end{verbatim} 
617 
caltinay 
4205 
in the script. 
618 


This will create a \file{result.csv} with columns separated by a comma. 
619 


Each row contains the value of the density distribution and the three 
620 


coordinates of the corresponding location in the domain. 
621 


There is a header specifying the meaning of the corresponding column. 
622 


Notice that rows are not written in a particular order and therefore, if 
623 


necessary, the user has to apply appropriate sorting of the rows. 
624 


Columns are written in alphabetic order of their corresponding tag names. 
625 


For the interested reader: the statement \verbrho.getFunctionSpace() returns 
626 


the type used to store the density data \verbrho. 
627 


The \verbgetX() method returns the coordinates of the sampling points used 
628 


for the particular type of representation, see~\cite{ESCRIPT} for details. 
629 
gross 
4203 

630 
gross 
4185 
\section{Remarks} 
631 



632 
caltinay 
4279 
\subsection{ER Mapper Raster Files}\label{SEC:P1:GRAV:REMARK:ERMAPPER} 
633 


The \downunder module can read data stored in ER Mapper Raster files. A data 
634 


set in this format consists of two files, a header file whose name usually ends 
635 


in \verb.ers and the actual data file which has the same filename as the 
636 


header file but without any file extension. 
637 


These files are usually produced by a commercial software package and the 
638 


contents can be quite diverse. 
639 


Therefore, it is not guaranteed that every data set is supported by \downunder 
640 


but the most common types of raster data should work\footnote{If your data 
641 


does not load please contact us through \url{https://launchpad.net/escriptfinley}.}. 
642 



643 


The interface for loading ER Mapper files is very similar to the \netcdf 
644 


interface described in Section~\ref{SEC:P1:GRAV:DATA}. 
645 


To load gravity data stored in the file pair\footnote{These files are available 
646 


in the example directory.} \examplefile{GravitySmall.ers} (the 
647 


header) and \examplefile{GravitySmall} (the data) without changing any of the 
648 


defaults use: 
649 


\begin{verbatim} 
650 


source0=ErMapperData(ErMapperData.GRAVITY, 'GravitySmall.ers') 
651 


\end{verbatim} 
652 


If your data set does not follow the default naming convention you can specify 
653 


the name of the data file explicitly: 
654 


\begin{verbatim} 
655 


source0=ErMapperData(ErMapperData.GRAVITY, 'GravitySmall.ers', 
656 


datafile='GravityData') 
657 


\end{verbatim} 
658 


Please note that there is no way for the reader to determine if the two files 
659 


really form a pair so make sure to pass the correct filenames when constructing 
660 


the reader object. 
661 


The same optional arguments explained in sections \ref{SEC:P1:GRAV:DATA} and 
662 


\ref{SEC:P1:GRAV:REMARK:DATAHOLES} are available for ER Mapper data sets. 
663 


However, due to the limitation of the file format only a constant error value 
664 


is supported. 
665 



666 
gross 
4199 
\subsection{Data With Holes}\label{SEC:P1:GRAV:REMARK:DATAHOLES} 
667 
caltinay 
4279 
As described previously in this chapter input data is always given in the form 
668 


of a rectangular grid with constant cell size in each dimension. 
669 
caltinay 
4271 
However, there are cases when this is not necessarily the case. 
670 


Consider an onshore data set which includes parts of the offshore region as in 
671 


Figure~\ref{FIG:P1:GRAV:onshoreoffshore}. 
672 


The valid data in this example has a value range of about $600$ to $600$ and 
673 


the inversion is to be run based on these values only, disregarding the 
674 


offshore region. 
675 


In order to achieve that, the offshore region is \emph{masked} by using a 
676 


constant value which is not found within the onshore area. 
677 


Figure~\ref{FIG:P1:GRAV:onshoreoffshore} clearly shows this masked area in 
678 


dark blue since a mask value of $1000$ was used. 
679 
gross 
4185 

680 
caltinay 
4271 
\begin{figure}[ht] 
681 


\centering 
682 


\includegraphics[width=0.45\textwidth]{onshoreoffshore.png} 
683 


\caption{Plot of a rectangular gridded onshore data set that includes 
684 


offshore regions which have a value (here $1000$) not found within the 
685 


real data (Bouguer anomalies in Tasmania, courtesy Geoscience Australia)} 
686 


\label{FIG:P1:GRAV:onshoreoffshore} 
687 


\end{figure} 
688 
gross 
4199 

689 
caltinay 
4271 
The \netcdf conventions supported in \downunder include a standard way of 
690 


specifying such a mask value. 
691 


The example script \examplefile{create_netcdf.py} demonstrates how this is 
692 


accomplished in an easy way with any data. 
693 


If, for any reason, the mask value in the input file is invalid it can be 
694 


overridden via the \verbnull_value argument when constructing the 
695 


\verbNetCdfData object: 
696 


\begin{verbatim} 
697 


source0=NetCdfData(NetCdfData.GRAVITY, 'data0.nc', null_value=1000) 
698 


\end{verbatim} 
699 


In this example, all data points that have a value of $1000$ are ignored and 
700 


not used in the inversion. 
701 


Please note that the special value \emph{NaN} (notanumber) is sometimes used 
702 


for the purposes of masking in data sets. 
703 


Areas marked with this value are always disregarded in \downunder. 
704 



705 
caltinay 
4279 
\subsection{Multiple Data Sets}\label{SEC:P1:GRAV:REMARK:MULTIDATA} 
706 


It is possible to run a single inversion using more than one input data set, 
707 


possibly in different file formats. 
708 


To do so, simply create the data sources and add them to the domain builder: 
709 
caltinay 
4272 
\begin{verbatim} 
710 
caltinay 
4279 
source0=NetCdfData(NetCdfData.GRAVITY, 'data0.nc') 
711 


source1=ErMapperData(ErMapperData.GRAVITY, 'data1.ers') 
712 


dom.addSource(source0) 
713 


dom.addSource(source1) 
714 
caltinay 
4272 
\end{verbatim} 
715 
caltinay 
4279 
However, there are some restrictions when combining data sets: 
716 


\begin{itemize} 
717 


\item Due to the coordinate transformation all data sets must be located in 
718 


the same UTM zone. If a single dataset crosses UTM zones only the zone 
719 


of the central longitude is used when projecting. 
720 


For example, if one data set lies mostly in zone 51 but contains areas 
721 


of zone 52, it is transformed using zone 51. In this case more data 
722 


from zone 51 can be added, but not from any other zone. 
723 


\item All data sets should have the same spatial resolution but this is not 
724 


enforced. Combining data with different resolution is currently 
725 


considered experimental but works best when the resolutions are 
726 


multiples of each other. For example if the first data set has a 
727 


resolution (or cell size) of $100$ metres and the second has a cell 
728 


size of $50$ metres then the target domain will have a cell size of 
729 


$50$ metres (the finer resolution) and each point of the coarse data 
730 


will occupy two cells (in the respective dimension). 
731 


\end{itemize} 
732 
caltinay 
4272 

733 
gross 
4199 
\subsection{Regularization Term}\label{SEC:P1:GRAV:REMARK:REG} 
734 
caltinay 
4205 
The \verbGravityInversion class supports the following form for the 
735 


regularization: 
736 
gross 
4199 
\begin{equation} 
737 


\int w^{(0)} \cdot \rho^2 + w^{(1)}_0 \rho_{,0}^2 + w^{(1)}_1 \rho_{,1}^2 + w^{(1)}_2 \rho_{,2}^2\; dx 
738 


\end{equation} 
739 
caltinay 
4205 
where the integral is calculated across the entire domain. 
740 


$\rho$ represents the density distribution where $\rho_{,0}$ $\rho_{,1}$ and 
741 


$\rho_{,2}$ are the spatial derivatives of $\rho$ with respect to the two 
742 


lateral and the vertical direction, respectively. 
743 


$w^{(0)}$, $w^{(1)}_0$, $w^{(1)}_1$ and $w^{(1)}_2$ are weighting 
744 


factors\footnote{A more general form, e.g. spatially variable values for the 
745 


weighting factors, is supported, see Part~\ref{part2}}. 
746 


By default these are $w^{(0)}=0$, $w^{(1)}_0=w^{(1)}_1=w^{(1)}_2=1$. 
747 


Other weighting factors can be set in the inversion setup. 
748 


For instance to set $w^{(0)}=10$, $w^{(1)}_0=w^{(1)}_1=0$ and $w^{(1)}_2=100$ 
749 


use the statement: 
750 
gross 
4199 
\begin{verbatim} 
751 


inv.setup(dom, w0=10, w1=[0,0,100]) 
752 


\end{verbatim} 
753 
caltinay 
4205 
It is pointed out that the weighting factors are rescaled in order to improve 
754 


numerical stability. Therefore the relative size of the weighting factors is 
755 


relevant and using 
756 
gross 
4199 
\begin{verbatim} 
757 


inv.setup(dom, w0=0.1, w1=[0,0,1]) 
758 


\end{verbatim} 
759 
caltinay 
4205 
would lead to the same regularization as the statement above. 
760 
gross 
4185 
