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

3 
\section{Introduction} 
4 

5 
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 
capacity are referred to the Part~\ref{part2} of this manual. 
13 
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 
sets. We refer to~\cite{REF1, REF2, REF3}. 
16 

17 
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 

31 
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 

41 
To run the examples discussed you need to have \escript (version 3.3.1 or 
42 
newer) installed on your computer. 
43 
Moreover, if you want to visualize the results you must have access to a 
44 
plotting software which is able to process \VTK input files, e.g. \mayavi and \VisIt. 
45 
As \mayavi can easily be obtained and installed for most platforms the 
46 
tutorial is based on \mayavi as visualization. However, it is pointed out 
47 
that \VisIt is the preferred visualization tool for \escript as it can deal 
48 
with very large data sets more efficiently. 
49 

50 
\begin{figure} 
51 
\centering 
52 
\includegraphics[width=0.7\textwidth]{QLDWestGravityDataPlot.png} 
53 
\caption{Gravity Anomaly Data in $mgal$ from Western Queensland, Australia 
54 
(see file \examplefile{data/QLDWest_grav.nc}). \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.}} 
55 
\label{FIG:P1:GRAV:0} 
56 
\end{figure} 
57 

58 
\begin{figure} 
59 
\centering% 
60 
\includegraphics[width=0.7\textwidth]{QLDGravMu10Contour.png}% 
61 
\caption{3D contour plot of the density distribution obtained by inversion of 
62 
file \file{data/QLDWest_grav.nc} (with $\mu=10$). 
63 
Colours represent values of density where high values are represented by 
64 
blue and low values are represented by red. \AZADEH{mu=10? add colour bar?}} 
65 
\label{FIG:P1:GRAV:1} 
66 
\end{figure} 
67 

68 
\section{How does it work?} 
69 
The execution of the inversion is controlled by a script which, in essence, is 
70 
a text file and can be edited using any text editor. 
71 
The script contains a series of statements which are executed by an 
72 
interpreter which is an executable program reading the text file and executing 
73 
the statements linebyline. 
74 
In the case of \downunder the interpreter is \Python. 
75 
In order to be able to process the statements in each line of the script 
76 
certain rules (called syntax) need to be obeyed. 
77 
There is a large number of online tutorials for \Python available\footnote{e.g. 
78 
\url{http://www.tutorialspoint.com/python} and \url{http://doc.pyschools.com}}. 
79 
We also refer to the \escript cook book \cite{ESCRIPTCOOKBOOK} and user's 
80 
guide \cite{ESCRIPT} which is in particular useful for users who like to dive 
81 
deeper into \downunder. 
82 
For this part of the manual no \Python knowledge is required but it is 
83 
recommended that users acquire some basic knowledge on \Python as they 
84 
progress in their work with \downunder. 
85 

86 
The following script~\ref{code: gravity1}\footnote{The script is similar to 
87 
\examplefile{grav_netcdf.py} within the \escript example file directory.} is a 
88 
simple example to run an inversion for gravity data: 
89 

90 
\begin{pyc}\label{code: gravity1} 
91 
\ 
92 
\begin{python} 
93 
# Header: 
94 
from esys.downunder import * 
95 
from esys.weipa import * 
96 
from esys.escript import unitsSI as U 
97 

98 
# Step 1: set up domain 
99 
dom=DomainBuilder() 
100 
dom.setVerticalExtents(depth=40.*U.km, air_layer=6.*U.km, num_cells=25) 
101 
dom.setFractionalPadding(pad_x=0.2, pad_y=0.2) 
102 
dom.fixDensityBelow(depth=40.*U.km) 
103 

104 
# Step 2: read gravity data 
105 
source0=NetCdfData(NetCdfData.GRAVITY, 'data/QLDWest_grav.nc', altitude=0.) 
106 
dom.addSource(source0) 
107 

108 
# Step 3: set up inversion 
109 
inv=GravityInversion() 
110 
inv.setSolverTolerance(1e4) 
111 
inv.setSolverMaxIterations(50) 
112 
inv.setup(dom) 
113 

114 
# Step 4: run inversion 
115 
inv.getCostFunction().setTradeOffFactorsModels(10.) 
116 
rho = inv.run() 
117 

118 
# Step 5: write reconstructed density to file 
119 
saveVTK("result.vtu", density=rho) 
120 
\end{python} 
121 
\end{pyc} 
122 
The result, in this case the density distribution, is written to an external 
123 
file for further processing. You can copy and paste the text of the script 
124 
into a file of any name, let's say for further reference we use the file 
125 
name \file{grav.py}. 
126 
It is recommended to use the extension \file{.py} to identify the file as a 
127 
\Python script. 
128 
We will discuss the statements of the script later in this chapter. 
129 

130 
Obviously the inversion needs to be fed with some gravity data. In this case 
131 
the data are loaded from the file \examplefile{data/QLDWest_grav.nc} which is 
132 
available in the \escript example file directory. 
133 
The data are given in the \netcdf file format for gridded data. 
134 
After you have copied this file into the directory in which you have saved the 
135 
script \file{grav.py} you can run the program using the command line 
136 
\begin{verbatim} 
137 
runescript grav.py 
138 
\end{verbatim} 
139 
We are running \file{grav.py} through the \escript startup command since 
140 
\escript is used as a back end for the inversion algorithm\footnote{Please see 
141 
the \escript user's guide~\cite{ESCRIPT} on how to run your script in parallel 
142 
using threading and/or MPI.}. 
143 
Obviously it is assumed that you have an installation of \escript available on 
144 
your computer, see \url{https://launchpad.net/escriptfinley}. 
145 

146 
After the execution has successfully completed you will find the result file 
147 
\file{result.vtu} in the directory from where you have started the execution 
148 
of the script. 
149 
The file has the \VTK format and can be imported easily into many 
150 
visualization tools. 
151 
One option is the \mayavi package which is available on most platforms. 
152 
You can invoke the the visualization using the commands 
153 
\begin{verbatim} 
154 
mayavi2 d result.vtu m SurfaceMap 
155 
\end{verbatim} 
156 
from the command line. 
157 
Figure~\ref{FIG:P1:GRAV:1} shows the result of this inversion as a contour 
158 
plot\footnote{These plots were generated by \VisIt.}, while the gravity 
159 
anomaly data is shown in Figure~\ref{FIG:P1:GRAV:0}. 
160 
We will discuss data later in Section~\ref{SEC:P1:GRAV:DATA}. 
161 

162 
Let us have a closer look at the script\footnote{In \Python lines starting 
163 
with `\#` are comments and are not processed further.}: Besides the header 
164 
section one can separate the script into five steps: 
165 
\begin{enumerate} 
166 
\item set up domain on which the inversion problem is solved 
167 
\item load the data 
168 
\item setup the inversion problem 
169 
\item run the inversion 
170 
\item further processing of the result. Here we write the reconstructed 
171 
density distribution to a file. 
172 
\end{enumerate} 
173 
In the following we will discuss the steps of the scripts in more detail. 
174 
Before we do this it is pointed out that the header section, following 
175 
\Python conventions, makes all required packages available to access within 
176 
the script. 
177 
At this point we will not discuss this in more details but emphasize that the 
178 
header section is a vital part of the script. 
179 
It is is required in each \downunder inversion script and should not be 
180 
altered except if additional modules are needed. 
181 

182 
\begin{figure} 
183 
\centering 
184 
\includegraphics[width=\textwidth]{dom2D.pdf} 
185 
\caption{2D domain setup for gravity inversion} 
186 
\label{FIG:P1:GRAV:2} 
187 
\end{figure} 
188 

189 
\section{Creating the Inversion Domain} 
190 
First step in script~\ref{code: gravity1} is the definition of the domain over 
191 
which the inversion is performed. 
192 
We think of the domain as a block with orthogonal, plain faces. 
193 
Figure~\ref{FIG:P1:GRAV:2} shows the setup for a twodimensional domain 
194 
(see also Figure~\ref{fig:domainBuilder} for 3D). 
195 
The lateral coordinates along the surface of the Earth are denoted by $x$ and 
196 
$y$ (only $x$direction is shown). 
197 
The $z$ direction defines the vertical direction where the part above the 
198 
surface has positive coordinates and the subsurface negative coordinates. 
199 
The height of the section above the surface, which is assumed to be filled 
200 
with air, needs to be set by the user. 
201 
The inversion assumes that the density in the section is known to be 
202 
zero\footnote{Always keeping in mind that these are not absolute values but 
203 
anomalies.}. 
204 
The density below the surface is unknown and is calculated through the 
205 
inversion. The user needs to specify the depth below the surface in which the 
206 
density is to be calculated. 
207 
The lateral extension of the domain is defined by the data sets fed into the 
208 
inversion. 
209 
It is chosen large enough to cover all data sets (in case more than one is 
210 
used). In order to reduce the impact of the boundary a padding zone around the 
211 
data sets can be introduced. 
212 

213 
\begin{figure} 
214 
\centering 
215 
\includegraphics[width=\textwidth]{dom2DFEM.pdf} 
216 
\caption{Cell distribution and boundary conditions for a 2D domain} 
217 
\label{FIG:P1:GRAV:3} 
218 
\end{figure} 
219 

220 
The reconstruction of the gravity field from an estimated density distribution 
221 
is the key component of the inversion process. 
222 
To do this \downunder uses the finite element method (FEM). 
223 
We need to introduce a grid over the domain, see Figure~\ref{FIG:P1:GRAV:3}. 
224 
The number of vertical cells is set by the user while the number of horizontal 
225 
cells is derived from the grid spacing of the gravity data set(s). 
226 
It is assumed that gravity field data given are constant across a cell. 
227 
To be able to reconstruct the gravity field some assumptions on the values of 
228 
the gravity field on the domain boundary have to be made. 
229 
\downunder assumes that on all faces the lateral gravity field component 
230 
equals zero. No assumptions on the horizontal components are 
231 
made\footnote{It is assumed that the gravity potential equals zero on the top 
232 
and bottom surface, see Section~\ref{sec:forward gravity} for details}% 
233 
\footnote{Most inversion codes use Green's functions over an unbounded domain 
234 
to reconstruct the gravity field. This approach makes the assumption that the 
235 
gravity field (or potential) is converging to zero when moving away from the 
236 
region of interest. The boundary conditions used here are stronger in the 
237 
sense that the lateral gravity component is enforced to be zero in a defined 
238 
distance of the region of interest but weaker in the sense that no constraint 
239 
on the horizontal component is applied.}. 
240 

241 
In script~\ref{code: gravity1} the statement 
242 
\begin{verbatim} 
243 
dom=DomainBuilder() 
244 
\end{verbatim} 
245 
creates something like a factory to build a domain. 
246 
We then define the features of the domain we would like to create: 
247 
\begin{verbatim} 
248 
dom.setVerticalExtents(depth=40.*U.km, air_layer=6.*U.km, num_cells=25) 
249 
dom.setFractionalPadding(pad_x=0.2, pad_y=0.2) 
250 
\end{verbatim} 
251 
Here we specify the depth of the domain to $40 km$, the thickness of the air 
252 
layer above the surface to $6 km$ and the number of vertical cells to $25$. 
253 
We also introduce a lateral padding of $20 \%$ of the expansion of the gravity 
254 
data on each side of the data and in both lateral directions. 
255 

256 
In some cases it can be appropriate to assume that the density below a certain 
257 
depth is zero\footnote{As we are in fact calculating density corrections this 
258 
means that the density is assumed to be the average density.}. 
259 
The statement 
260 
\begin{verbatim} 
261 
dom.fixDensityBelow(depth=40.*U.km) 
262 
\end{verbatim} 
263 
introduces this constraint. 
264 
As in the case discussed here the depth for zero density is not less than the 
265 
depth of the domain no constraint at depth is applied to the density. 
266 

267 
\downunder uses the metrekilogramsecond based International System of Units 
268 
(SI)~\footnote{see \url{http://en.wikipedia.org/wiki/International_System_of_Units}}\index{SI}. 
269 
So all values must be converted to appropriate units. 
270 
This does not apply to geographic coordinates where \downunder uses the degree 
271 
with fractional degree (as a floating point number) to represent longitude and 
272 
latitude. In the script we have used the expression 
273 
\begin{verbatim} 
274 
depth=40.*U.km 
275 
\end{verbatim} 
276 
to define the depth of the domain to $40 km$. 
277 
The expression \verbU.km is the unit $km$ (kilometer) and an appropriate 
278 
conversion of $40$ as length in $km$ into the base unit $m$ (meter). 
279 
It is recommended to add units to values (where present) in order to make sure 
280 
that the final values handed into \downunder is given with the appropriate 
281 
units. 
282 
The physical units module of \escript, which we have imported here under the 
283 
name \verbU in the script header, defines a large number of physical units 
284 
and constants, please see~\cite{ESCRIPT} and~\cite{ESCRIPTONLINE}. 
285 

286 
\section{Loading Gravitational Data}\label{SEC:P1:GRAV:DATA} 
287 
In practice gravity acceleration is measured in various ways, for instance by 
288 
airborne surveys~\cite{Telford1990a}. 
289 
\downunder assumes that all data supplied as input are already appropriately 
290 
preprocessed. In particular, corrections for 
291 
\begin{itemize} 
292 
\item freeair, to remove effects from altitude above ground; 
293 
\item latitude, to remove effects from ellipsoidicity of the Earth; 
294 
\item terrain, to remove effects from topography 
295 
\end{itemize} 
296 
must have been applied to the data. 
297 
In general, data prepared in such a form are called Bouguer anomalies~\cite{Telford1990a}. 
298 

299 
To load gravity data into \downunder the data are given on a plane parallel 
300 
to the surface of the Earth at a constant altitude, see 
301 
diagram~\ref{FIG:P1:GRAV:2}. 
302 
The data need to be defined over a rectangular grid covering a subregion of 
303 
the Earth surface. 
304 
The grid uses the geographic coordinate system. 
305 
Figure~\ref{FIG:P1:GRAV:0} shows an example of such a data set from the 
306 
western parts of Queensland, Australia. 
307 
The data set covers a rectangular region between $140^o$ and $141^o$ east 
308 
and between $20^o$ and $21^o$ south. 
309 
Notice that latitude varies between $90^o$ to $90^o$ where negative signs 
310 
refer to places in the southern hemisphere and longitude varies between 
311 
$180^o$ to $180^o$ where negative signs refer to places in the western 
312 
hemisphere. 
313 
The colour at a location represents the value of the vertical Bouguer gravity 
314 
anomaly at this point at the surface of the Earth. 
315 
Values range from $160 \; mgal$ to about $500 \; mgal$\footnote{The unit 
316 
$mgal$ means milli $gal$ (galileo) with $1 \; gal = 0.01 \frac{m}{sec^2}$.} 
317 
over a $121 \times 121$ grid. 
318 

319 
In general, a patch of gravity data need to be defined over a plane 
320 
\verbNX $\times$ \verbNY where \verbNX and \verbNY define the number 
321 
of grid lines in the longitude (\verbX) and the latitude (\verbY) 
322 
direction, respectively. 
323 
The grid is spanned from an origin with spacing \verbDELTA_X and 
324 
\verbDELTA_Y in the longitude and the latitude direction, respectively. 
325 
The gravity data for all grid points need to be given as an \verbNX 
326 
$\times$ \verbNY array. 
327 
If available errors can be associated with the gravity data. 
328 
The values are given as an \verbNX $\times$ \verbNY array matching the 
329 
shape of the gravity array. 
330 
The \Python script \examplefile{create_netcdf.py} shows how to create a file 
331 
in the \netcdf file format~\cite{NETCDF} compatible with \downunder from 
332 
a data array. 
333 
see Section~\ref{SEC:P1:GRAV:REMARK:DATAHOLES} 
334 
BLA 
335 

336 
\CIHAN{ADD an example how ro use ER Mapper raster files} 
337 

338 
In addition to the \netcdf file format 
339 

340 

341 
Two data formats are supported, namely ER Mapper Raster File and \netcdf. 
342 
Here we will focus on the \netcdf file format~\cite{NETCDF} and refer to 
343 
subsection~\ref{SEC:P1:GRAV:REMARK:ERMAPPER} and 
344 
Section~\ref{sec:ref:DataSource:ERM} for more information on the ER Mapper 
345 
Raster format. 
346 

347 
BLA 
348 
In script~\ref{code: gravity1} we use the statement 
349 
\begin{verbatim} 
350 
source0=NetCdfData(NetCdfData.GRAVITY, 'data/QLDWest_grav.nc', altitude=0.) 
351 
\end{verbatim} 
352 
\CIHAN{REWISE} \AZADEH{ADD coarse data 10x10 or so?} 
353 
to load the gravity data stored in \examplefile{data/QLDWest_grav.nc} in the \netcdf format. 
354 
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 
\begin{verbatim} 
357 
dom.addSource(source0) 
358 
\end{verbatim} 
359 
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 
It is possible to add several, possibly overlapping data sets to the same 
369 
domain: 
370 
\begin{verbatim} 
371 
source0=NetCdfData(NetCdfData.GRAVITY, 'data0.nc', altitude=0.) 
372 
source1=NetCdfData(NetCdfData.GRAVITY, 'data1.nc', altitude=0.) 
373 
dom.addSource(source0) 
374 
dom.addSource(source1) 
375 
\end{verbatim} 
376 
However, due to the coordinate transformation all data sets must be located in 
377 
the same UTM zone. If a single dataset crosses UTM zones only the zone of the 
378 
central longitude is used when projecting. 
379 
For example, if one data set lies mostly in zone 51 but contains areas of zone 
380 
52, it is transformed using zone 51. 
381 
In this case more data from zone 51 can be added, but not from any other zone. 
382 

383 
It is important to keep an eye on the complexity of the inversion. 
384 
A good measure is the total number of cells being used. 
385 
Assume we have given a data set on a $20 \times 20$ grid and we add lateral 
386 
padding of, say, $20 \%$ to each side of the data, the lateral number of cells 
387 
becomes $(20\cdot 1.4)\times (20\cdot 1.4)=1.4^2\cdot 20^2\approx 2\cdot 10^2=800$. 
388 
If we use $20$ cells in the vertical direction we end up with a total number 
389 
of $800 \times 20 = 16,000$ cells. 
390 
This size can be easily handled by a modern desktop PC. 
391 
If we increase the grid size of the data to $40 \times 40$ points and use $40$ 
392 
cells in the vertical extent we get a total of $(2\cdot 40^2)\cdot 40=128,000$ 
393 
cells, a problem size which is considerably larger but can still be handled by 
394 
a desktop computer. 
395 
Taking this one step further, if the amount of data is increased to 
396 
$200\times 200$ points and we use $200$ cells in the vertical extent the 
397 
domain will contain $16,000,000$ ($16$ million) cells. 
398 
This scenario requires a computer with enough memory and (a) fast processor(s) 
399 
to run the inversion. 
400 
This estimate of complexity growth applies to the case where the increase 
401 
of data grid size is driven by an increase of resolution where it is 
402 
recommended to increase the vertical resolution in synch with the lateral 
403 
resolution. Note that if more than one data set is used the target resolution 
404 
will be the resolution of the finest data set. 
405 
In other cases the expansion of the region of interest drives an increase of 
406 
data grid size and the increase of total number of cells is less dramatic as 
407 
the vertical number of cells can remain constant while keeping a balanced 
408 
resolution in vertical and lateral direction. 
409 

410 
\section{Setting up the Inversion and Running it} 
411 
We are now at step three of script~\ref{code: gravity1} in which the actual 
412 
inversion is set up. 
413 
First we create an empty inversion under the name \verbinv: 
414 
\begin{verbatim} 
415 
inv=GravityInversion() 
416 
\end{verbatim} 
417 
As indicated by the name we can use \verbinv to perform an inversion of 
418 
gravity data\footnote{\verbGravityInversion is a driver with a simplified 
419 
interface which is provided for convenience. See Part~\ref{part2} for more 
420 
details on how to write inversion scripts with more general functionality, e.g. 
421 
constraints.}. The inversion is an iterative process which sequentially 
422 
calculates updates to the density distribution in an attempt to improve the 
423 
match of the gravity field produced by the density distribution with the data. 
424 
Termination of the iteration is controlled by the tolerance which is set by 
425 
the user: 
426 
\begin{verbatim} 
427 
inv.setSolverTolerance(1e4) 
428 
\end{verbatim} 
429 
Here we set the tolerance to $10^{4}$, i.e. the iteration is terminated if 
430 
the maximum density correction is less than or equal to $10^{4}$ relative to 
431 
the maximum value of estimated density anomaly. 
432 
In case the iteration does not converge a maximum number of iteration steps is 
433 
set: 
434 
\begin{verbatim} 
435 
inv.setSolverMaxIterations(50) 
436 
\end{verbatim} 
437 
If the maximum number of iteration steps (here $50$) is reached the iteration 
438 
process is aborted and an error message is printed. 
439 
In this case you can try to rerun the inversion with a larger value for the 
440 
maximum number of iteration steps. 
441 
If even for a very large number of iteration steps no convergence is achieved, 
442 
this is a strong indicator that the inversion has not been set up properly. 
443 

444 
The statement 
445 
\begin{verbatim} 
446 
inv.setup(dom) 
447 
\end{verbatim} 
448 
links the inversion with the domain and the data. 
449 
At this step  as we are solving a gravity inversion problem  only 
450 
gravitational data attached to the domain builder \verbdom are considered. 
451 
Internally a cost function $J$ is created which is minimized during the 
452 
inversion iteration. 
453 
It is a combination of a measure of the data misfit of the gravity field from 
454 
the given density distribution and a measure of the smoothness of the density 
455 
distribution. 
456 
The latter is often called the regularization term. 
457 
By default the gradient of density is used as the regularization term, see 
458 
also Section~\ref{SEC:P1:GRAV:REMARK:REG}. 
459 
Obviously, the result of the inversion is sensitive to the weighting between 
460 
the misfit and the regularization. 
461 
This tradeoff factor $\mu$ for the misfit function is set by the following 
462 
statement: 
463 
\begin{verbatim} 
464 
inv.getCostFunction().setTradeOffFactorsModels(0.1) 
465 
\end{verbatim} 
466 
Here we set $\mu=0.1$. The statement \verbinv.setup must appear in the 
467 
script before setting the tradeoff factor. 
468 
A small value for the tradeoff factor $\mu$ will give more emphasis to the 
469 
regularization component and create a smoother density distribution. 
470 
A large value of the tradeoff factor $\mu$ will emphasize the misfit more 
471 
and typically creates a better fit to the data and a rougher density 
472 
distribution. 
473 
It is important to keep in mind that the regularization reduces noise in the 
474 
date and in fact gives the problem a unique solution. 
475 
Consequently, the tradeoff factor $\mu$ may not be chosen too large in order 
476 
control the noise on the solution and ensure convergence in the iteration 
477 
process. 
478 

479 
We can now actually run the inversion: 
480 
\begin{verbatim} 
481 
rho = inv.run() 
482 
\end{verbatim} 
483 
The answer as calculated during the inversion is returned and can be accessed 
484 
under the name \verbrho. 
485 
As pointed out earlier the iteration process may fail in which case the 
486 
execution of the script is aborted with an error message. 
487 

488 
\section{Taking a Look} 
489 
In the final step of script~\ref{code: gravity1} the calculated density 
490 
distribution is written to an external file. 
491 
A popular file format used by several visualization packages such as 
492 
\VisIt~\cite{VISIT} and \mayavi~\cite{MAYAVI} is the \VTK file format. 
493 
The result of the inversion which has been named \verbrho can be written to 
494 
the file \file{result.vtu} by adding the statement 
495 
\begin{verbatim} 
496 
saveVTK("result.vtu", density=rho) 
497 
\end{verbatim} 
498 
at the end of script. 
499 
The inversion solution is tagged with the name \verbdensity in the result 
500 
file, however any other name for the tag could be used. 
501 
As the format is textbased (as opposed to binary) \VTK files tend to be very 
502 
large and take compute time to create, in particular when it comes to large 
503 
numbers of cells ($>10^6$). 
504 
For large problems it is more efficient to use the \SILO file format~\cite{SILO}. 
505 
\SILO files tend to be smaller and are faster generated and read. 
506 
It is the preferred format to import results into the visualization program 
507 
\VisIt~\cite{VISIT} which is particularly suited for the visualization of 
508 
large data sets. 
509 
Inversion results can directly exported into \SILO files using the statement 
510 
\begin{verbatim} 
511 
saveSilo("result.silo", density=rho) 
512 
\end{verbatim} 
513 
replacing the \verbsaveVTK(...) statement. 
514 
Similar to \VTK files the result \verbrho is tagged with the name 
515 
\verbdensity so it can be identified in the visualization program. 
516 

517 
\begin{figure} 
518 
\begin{center} 
519 
\subfigure[$\mu=0.1$]{% 
520 
\label{FIG:P1:GRAV:10 MU01} 
521 
\includegraphics[width=0.45\textwidth]{QLDGravMu01Contour.png} 
522 
}% 
523 
\subfigure[$\mu=1.$]{% 
524 
\label{FIG:P1:GRAV:10 MU1} 
525 
\includegraphics[width=0.45\textwidth]{QLDGravMu1Contour.png} 
526 
}\\ %  End of the first row % 
527 
\subfigure[$\mu=10.$]{% 
528 
\label{FIG:P1:GRAV:10 MU10} 
529 
\includegraphics[width=0.45\textwidth]{QLDGravMu10Contour.png} 
530 
}% 
531 
\subfigure[$\mu=100.$]{% 
532 
\label{FIG:P1:GRAV:10 MU100} 
533 
\includegraphics[width=0.45\textwidth]{QLDGravMu100Contour.png} 
534 
}% 
535 
\end{center} 
536 
\caption{3D contour plots of gravity inversion results with data from 
537 
Figure~\ref{FIG:P1:GRAV:0} for various values of the model tradeoff 
538 
factor $\mu$. Visualization has been performed in \VisIt. 
539 
\AZADEH{check images.}} 
540 
\label{FIG:P1:GRAV:10} 
541 
\end{figure} 
542 

543 
\begin{figure} 
544 
\begin{center} 
545 
\subfigure[$\mu=0.1$]{% 
546 
\label{FIG:P1:GRAV:11 MU01} 
547 
\includegraphics[width=0.45\textwidth]{QLDGravMu01Slice.png} 
548 
}% 
549 
\subfigure[$\mu=1.$]{% 
550 
\label{FIG:P1:GRAV:11 MU1} 
551 
\includegraphics[width=0.45\textwidth]{QLDGravMu1Slice.png} 
552 
}\\ %  End of the first row % 
553 
\subfigure[$\mu=10.$]{% 
554 
\label{FIG:P1:GRAV:11 MU10} 
555 
\includegraphics[width=0.45\textwidth]{QLDGravMu10Slice.png} 
556 
}% 
557 
\subfigure[$\mu=100.$]{% 
558 
\label{FIG:P1:GRAV:11 MU100} 
559 
\includegraphics[width=0.45\textwidth]{QLDGravMu100Slice.png} 
560 
}% 
561 
\end{center} 
562 
\caption{3D slice plots of gravity inversion results with data from 
563 
Figure~\ref{FIG:P1:GRAV:0} for various values of the model tradeoff 
564 
factor $\mu$. Visualization has been performed \VisIt. 
565 
\AZADEH{check images.}} 
566 
\label{FIG:P1:GRAV:11} 
567 
\end{figure} 
568 

569 
Figures~\ref{FIG:P1:GRAV:10} and~\ref{FIG:P1:GRAV:11} show two different 
570 
styles of visualization generated in \VisIt using the result of the inversion 
571 
of the gravity anomalies shown in Figure~\ref{FIG:P1:GRAV:0}. 
572 
The inversions have been performed with different values for the model 
573 
tradeoff factor $\mu$. 
574 
The visualization shows clearly the smoothing effect of lower values for the 
575 
tradeoff factors. 
576 
For larger values of the tradeoff factor the density distribution becomes 
577 
rougher showing larger details. 
578 
Computational costs are significantly higher for larger tradeoff factors. 
579 
Moreover, noise in the data has an higher impact on the result. 
580 
Typically several runs are required to adjust the value for the tradeoff 
581 
factor to the datasets used. 
582 

583 
For some analysis tools it is useful to process the results in form of Comma 
584 
Separated Values (\CSV)~\footnote{see 
585 
\url{http://en.wikipedia.org/wiki/Commaseparated_values}.}. 
586 
Such a file can be created using the statement 
587 
\begin{verbatim} 
588 
saveCSV("result.csv", x=rho.getFunctionSpace().getX(), density=rho) 
589 
\end{verbatim} 
590 
in the script. 
591 
This will create a \file{result.csv} with columns separated by a comma. 
592 
Each row contains the value of the density distribution and the three 
593 
coordinates of the corresponding location in the domain. 
594 
There is a header specifying the meaning of the corresponding column. 
595 
\LG{Add example output CSV}. 
596 
Notice that rows are not written in a particular order and therefore, if 
597 
necessary, the user has to apply appropriate sorting of the rows. 
598 
Columns are written in alphabetic order of their corresponding tag names. 
599 
For the interested reader: the statement \verbrho.getFunctionSpace() returns 
600 
the type used to store the density data \verbrho. 
601 
The \verbgetX() method returns the coordinates of the sampling points used 
602 
for the particular type of representation, see~\cite{ESCRIPT} for details. 
603 

604 

605 
\section{Remarks} 
606 

607 
\subsection{Data With Holes}\label{SEC:P1:GRAV:REMARK:DATAHOLES} 
608 
\CIHAN{ADD some remarks holes in the data, example?} 
609 

610 

611 
\subsection{ER Mapper Raster File}\label{SEC:P1:GRAV:REMARK:ERMAPPER} 
612 
\CIHAN{ADD an example how to use ER Mapper raster files} 
613 

614 
\subsection{Regularization Term}\label{SEC:P1:GRAV:REMARK:REG} 
615 
The \verbGravityInversion class supports the following form for the 
616 
regularization: 
617 
\begin{equation} 
618 
\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 
619 
\end{equation} 
620 
where the integral is calculated across the entire domain. 
621 
$\rho$ represents the density distribution where $\rho_{,0}$ $\rho_{,1}$ and 
622 
$\rho_{,2}$ are the spatial derivatives of $\rho$ with respect to the two 
623 
lateral and the vertical direction, respectively. 
624 
$w^{(0)}$, $w^{(1)}_0$, $w^{(1)}_1$ and $w^{(1)}_2$ are weighting 
625 
factors\footnote{A more general form, e.g. spatially variable values for the 
626 
weighting factors, is supported, see Part~\ref{part2}}. 
627 
By default these are $w^{(0)}=0$, $w^{(1)}_0=w^{(1)}_1=w^{(1)}_2=1$. 
628 
Other weighting factors can be set in the inversion setup. 
629 
For instance to set $w^{(0)}=10$, $w^{(1)}_0=w^{(1)}_1=0$ and $w^{(1)}_2=100$ 
630 
use the statement: 
631 
\begin{verbatim} 
632 
inv.setup(dom, w0=10, w1=[0,0,100]) 
633 
\end{verbatim} 
634 
It is pointed out that the weighting factors are rescaled in order to improve 
635 
numerical stability. Therefore the relative size of the weighting factors is 
636 
relevant and using 
637 
\begin{verbatim} 
638 
inv.setup(dom, w0=0.1, w1=[0,0,1]) 
639 
\end{verbatim} 
640 
would lead to the same regularization as the statement above. 
641 

642 
% 
643 
% \section{Reference} 
644 
% 
645 
% There are three examples for 2D and 3D gravity inversions with artificial input data. 
646 
% In first step, an area with synthetic density section was suggested. Then based on forward modeling its gravitational data was collected. Afterwards with generated gravity data, escripts simulated a volume of inverted density. Whilst new density mass could be compared with the synthetic density section to verify the inversion. 
647 
% 
648 
% Some of the presumptions were the same for all of the examples to simplify the situation to make a logical comparison between synthetic input and output. which is as followed: 
649 
% 
650 
% \begin{verbatim} 
651 
% n_cells_in_data=100 
652 
% depth_offset=0.*U.km 
653 
% l_data = 100 * U.km 
654 
% l_pad=40*U.km 
655 
% THICKNESS=20.*U.km 
656 
% l_air=6*U.km 
657 
% \end{verbatim} 
658 
% 
659 
% The others assumptions comes with each examples. 
660 
% \begin{enumerate} 
661 
% \item A 2D density section with a maximum in center was assumed. The reference density and inverted will be shown. The padding area is excluded. (\ref{fig:gravity2D1}) 
662 
% 
663 
% \begin{verbatim} 
664 
% n_cells_in_data=100 
665 
% n_humbs_h= 1 
666 
% n_humbs_v=1 
667 
% mu=100 
668 
% \end{verbatim} 
669 
% 
670 
% \begin{figure} 
671 
% \centering 
672 
% \includegraphics[width=\textwidth]{grav2D1.png} 
673 
% \caption{2D density model up) reference down) result} 
674 
% \label{fig:gravity2D1} 
675 
% \end{figure} 
676 
% 
677 
% 
678 
% \item A 2D density properties with two maximum in corners and one minimum in the center was inverted. The result have eliminated the effects in padding area.(\ref{fig:gravity2D3}) 
679 
% 
680 
% \begin{verbatim} 
681 
% n_cells_in_data=100 
682 
% n_humbs_h= 3 
683 
% n_humbs_v=1 
684 
% mu=100 
685 
% \end{verbatim} 
686 
% 
687 
% \begin{figure} 
688 
% \centering 
689 
% \includegraphics[width=\textwidth]{grav2D3.png} 
690 
% \caption{2D density model up) reference down) result} 
691 
% \label{fig:gravity2D3} 
692 
% \end{figure} 
693 
% 
694 
% \item A 3D model with a maximum in the center was used as input data and the result after simulation in shown in the next image which determined a good distribution of density through the model in the main area.(\ref{fig:gravity3D} and \ref{fig:gravity3D1}) 
695 
% 
696 
% \begin{verbatim} 
697 
% n_cells_in_data=50 
698 
% n_humbs_h= 1 
699 
% n_humbs_v=1 
700 
% mu=10 
701 
% \end{verbatim} 
702 
% 
703 
% \begin{figure} 
704 
% \centering 
705 
% \includegraphics[width=\textwidth]{density3Dref.png} 
706 
% \caption{3D density model of reference as synthetic data} 
707 
% \label{fig:gravity3D} 
708 
% \end{figure} 
709 
% 
710 
% \begin{figure} 
711 
% \centering 
712 
% \includegraphics[width=\textwidth]{gravity3D.png} 
713 
% \caption{3D density model of result} 
714 
% \label{fig:gravity3D1} 
715 
% \end{figure} 
716 
% \end{enumerate} 
717 
