/[escript]/trunk/doc/inversion/CookGravity.tex
ViewVC logotype

Annotation of /trunk/doc/inversion/CookGravity.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4286 - (hide annotations)
Thu Mar 7 04:28:11 2013 UTC (6 years, 6 months ago) by caltinay
File MIME type: application/x-tex
File size: 38035 byte(s)
Assorted spelling fixes.

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 3-D distribution of density and/or susceptibility from 2-D 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{3-D 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 line-by-line.
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(1e-4)
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     run-escript grav.py
140     \end{verbatim}
141 caltinay 4205 We are running \file{grav.py} through the \escript start-up 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/escript-finley}.
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 set-up 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{2-D domain set-up 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 set-up for a two-dimensional domain
196     (see also Figure~\ref{fig:domainBuilder} for 3-D).
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 2-D).
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 2-D 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 metre-kilogram-second 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 \verb|U.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 \verb|U| 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     pre-processed. In particular, corrections for
293 gross 4191 \begin{itemize}
294 caltinay 4205 \item free-air, 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 \verb|NX| $\times$ \verb|NY| where \verb|NX| and \verb|NY| define the number
323     of grid lines in the longitude (\verb|X|) and the latitude (\verb|Y|)
324     direction, respectively.
325     The grid is spanned from an origin with spacing \verb|DELTA_X| and
326     \verb|DELTA_Y| in the longitude and the latitude direction, respectively.
327     The gravity data for all grid points need to be given as an \verb|NX|
328     $\times$ \verb|NY| array.
329 caltinay 4266 If available, measurement errors can be associated with the gravity data.
330 caltinay 4205 The values are given as an \verb|NX| $\times$ \verb|NY| 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 \verb|source0|.
355     We need to link the data set to the \verb|DomainBuilder| 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 \verb|DomainBuilder|
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=1e-5)
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     \verb|errors|, 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 \verb|inv|:
433 gross 4199 \begin{verbatim}
434     inv=GravityInversion()
435     \end{verbatim}
436 caltinay 4205 As indicated by the name we can use \verb|inv| to perform an inversion of
437     gravity data\footnote{\verb|GravityInversion| 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(1e-4)
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 \verb|dom| 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 trade-off 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 \verb|inv.setup| must appear in the
486     script before setting the trade-off factor.
487     A small value for the trade-off factor $\mu$ will give more emphasis to the
488     regularization component and create a smoother density distribution.
489     A large value of the trade-off 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 trade-off 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 \verb|rho|.
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 \verb|rho| 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 \verb|density| in the result
519     file, however any other name for the tag could be used.
520     As the format is text-based (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 \verb|saveVTK(...)| statement.
533     Similar to \VTK files the result \verb|rho| is tagged with the name
534     \verb|density| 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{3-D contour plots of gravity inversion results with data from
560     Figure~\ref{FIG:P1:GRAV:0} for various values of the model trade-off
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{3-D slice plots of gravity inversion results with data from
590     Figure~\ref{FIG:P1:GRAV:0} for various values of the model trade-off
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     trade-off factor $\mu$.
601     The visualization shows clearly the smoothing effect of lower values for the
602     trade-off factors.
603     For larger values of the trade-off factor the density distribution becomes
604     rougher showing larger details.
605     Computational costs are significantly higher for larger trade-off 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 trade-off
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 Comma-separated Values (\CSV)\footnote{see
612     \url{http://en.wikipedia.org/wiki/Comma-separated_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 \verb|rho.getFunctionSpace()| returns
626     the type used to store the density data \verb|rho|.
627     The \verb|getX()| 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/escript-finley}.}.
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]{onshore-offshore.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 \verb|null_value| argument when constructing the
695     \verb|NetCdfData| 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} (not-a-number) 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 \verb|GravityInversion| 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 set-up.
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

  ViewVC Help
Powered by ViewVC 1.1.26