ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 4268 - (show annotations)
Mon Mar 4 00:41:08 2013 UTC (6 years, 1 month ago) by gross
File MIME type: application/x-tex
File size: 34481 byte(s)
some citation fixed.

1 \chapter{Gravity Inversion}\label{Chp:cook:gravity inversion}
3 \section{Introduction}
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 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.
17 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.
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.
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 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.
51 \begin{figure}
52 \centering
53 \includegraphics[width=0.7\textwidth]{QLDWestGravityDataPlot.png}
54 \caption{Gravity Anomaly Data in $mgal$ from Western Queensland, Australia
55 (see file \examplefile{data/QLDWest_grav.nc}). Data obtained from Geoscience Australia.}
56 %\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 \label{FIG:P1:GRAV:0}
58 \end{figure}
60 \begin{figure}
61 \centering%
62 \includegraphics[width=0.9\textwidth]{QLDGravContourMu10.png}
63 \caption{3-D contour plot of the density distribution obtained by inversion of
64 file \file{data/QLDWest_grav.nc} (with $\mu=10$).
65 Colours represent values of density where high values are represented by
66 blue and low values are represented by red.}
67 \label{FIG:P1:GRAV:1}
68 \end{figure}
70 \section{How does it work?}
71 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.
88 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:
92 \begin{pyc}\label{code: gravity1}
93 \
94 \begin{python}
95 # Header:
96 from esys.downunder import *
97 from esys.weipa import *
98 from esys.escript import unitsSI as U
100 # Step 1: set up domain
101 dom=DomainBuilder()
102 dom.setVerticalExtents(depth=40.*U.km, air_layer=6.*U.km, num_cells=25)
103 dom.setFractionalPadding(pad_x=0.2, pad_y=0.2)
104 dom.fixDensityBelow(depth=40.*U.km)
106 # Step 2: read gravity data
107 source0=NetCdfData(NetCdfData.GRAVITY, 'QLDWest_grav.nc', altitude=0.)
108 dom.addSource(source0)
110 # Step 3: set up inversion
111 inv=GravityInversion()
112 inv.setSolverTolerance(1e-4)
113 inv.setSolverMaxIterations(50)
114 inv.setup(dom)
116 # Step 4: run inversion
117 inv.getCostFunction().setTradeOffFactorsModels(10.)
118 rho = inv.run()
120 # Step 5: write reconstructed density to file
121 saveVTK("result.vtu", density=rho)
122 \end{python}
123 \end{pyc}
124 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.
132 Obviously the inversion needs to be fed with some gravity data. In this case
133 the data are loaded from the file \examplefile{QLDWest_grav.nc} which is
134 available in the \escript example file directory.
135 The data in this example are given in the \netcdf file format.
136 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 \begin{verbatim}
139 run-escript grav.py
140 \end{verbatim}
141 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}.
148 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 You can invoke the the visualization using the commands
155 \begin{verbatim}
156 mayavi2 -d result.vtu -m SurfaceMap
157 \end{verbatim}
158 from the command line.
159 Figure~\ref{FIG:P1:GRAV:1} shows the result of this inversion as a contour
160 plot\footnote{These plots were generated by \VisIt.}, while the gravity
161 anomaly data is shown in Figure~\ref{FIG:P1:GRAV:0}.
162 We will discuss data later in Section~\ref{SEC:P1:GRAV:DATA}.
164 Let us have a closer look at the script\footnote{In \Python lines starting
165 with `\#` are comments and are not processed further.}: Besides the header
166 section one can separate the script into five steps:
167 \begin{enumerate}
168 \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 \end{enumerate}
175 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.
184 \begin{figure}
185 \centering
186 \includegraphics[width=\textwidth]{dom2D.pdf}
187 \caption{2-D domain set-up for gravity inversion}
188 \label{FIG:P1:GRAV:2}
189 \end{figure}
191 \section{Creating the Inversion Domain}
192 First step in script~\ref{code: gravity1} is the definition of the domain over
193 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 $y$ (only $x$-direction is used in 2-D).
199 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.
215 \begin{figure}
216 \centering
217 \includegraphics[width=\textwidth]{dom2DFEM.pdf}
218 \caption{Cell distribution and boundary conditions for a 2-D domain}
219 \label{FIG:P1:GRAV:3}
220 \end{figure}
222 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.}.
243 In script~\ref{code: gravity1} the statement
244 \begin{verbatim}
245 dom=DomainBuilder()
246 \end{verbatim}
247 creates something like a factory to build a domain.
248 We then define the features of the domain we would like to create:
249 \begin{verbatim}
250 dom.setVerticalExtents(depth=40.*U.km, air_layer=6.*U.km, num_cells=25)
251 dom.setFractionalPadding(pad_x=0.2, pad_y=0.2)
252 \end{verbatim}
253 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.
258 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 \begin{verbatim}
263 dom.fixDensityBelow(depth=40.*U.km)
264 \end{verbatim}
265 introduces this constraint.
266 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.
269 \downunder uses the metre-kilogram-second based International System of Units
270 (SI)~\footnote{see \url{http://en.wikipedia.org/wiki/International_System_of_Units}}\index{SI}.
271 So all values must be converted to appropriate units.
272 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 latitude. In the script we have used the expression
275 \begin{verbatim}
276 depth=40.*U.km
277 \end{verbatim}
278 to define the depth of the domain to $40 km$.
279 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 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}.
288 \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 \begin{itemize}
294 \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 \end{itemize}
298 must have been applied to the data.
299 In general, data prepared in such a form are called Bouguer anomalies~\cite{Telford1990a}.
301 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 The grid uses a geographic coordinate system with latitudes and longitudes
307 assumed to be given in the Clarke 1866 geodetic system.
308 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 $-180^o$ to $180^o$ where negative signs refer to places in the western
315 hemisphere.
316 The colour at a location represents the value of the vertical Bouguer gravity
317 anomaly at this point at the surface of the Earth.
318 Values in this data set range from $-160 \; mgal$ to about $500 \; mgal$\footnote{The unit
319 $mgal$ means milli $gal$ (galileo) with $1 \; gal = 0.01 \frac{m}{sec^2}$.}
320 over a $121 \times 121$ grid.
322 In general, a patch of gravity data need to be defined over a plane
323 \verb|NX| $\times$ \verb|NY| where \verb|NX| and \verb|NY| define the number
324 of grid lines in the longitude (\verb|X|) and the latitude (\verb|Y|)
325 direction, respectively.
326 The grid is spanned from an origin with spacing \verb|DELTA_X| and
327 \verb|DELTA_Y| in the longitude and the latitude direction, respectively.
328 The gravity data for all grid points need to be given as an \verb|NX|
329 $\times$ \verb|NY| array.
330 If available, measurement errors can be associated with the gravity data.
331 The values are given as an \verb|NX| $\times$ \verb|NY| array matching the
332 shape of the gravity array.
333 Note that data need not be available on every single point of the grid, see
334 Section~\ref{SEC:P1:GRAV:REMARK:DATAHOLES} for more information on this.
336 Currently, two data file formats are supported, namely \emph{ER Mapper Raster}
337 \cite{ERMAPPER} files and \netcdf~\cite{NETCDF} files.
338 In the examples of this chapter we use \netcdf files and refer to
339 Section~\ref{SEC:P1:GRAV:REMARK:ERMAPPER} and
340 Section~\ref{sec:ref:DataSource:ERM} for more information on using ER Mapper
341 Raster files.
342 If you have data in any other format you have the option of writing a suitable
343 reader (for advanced users, see Chapter~\ref{Chp:ref:data sources}) or,
344 assuming you are able to read the data in \Python, refer to the example
345 script \examplefile{create_netcdf.py} which shows how to create a file
346 in the \netcdf file format~\cite{NETCDF} compatible with \downunder from
347 a data array.
349 \CIHAN{Use coarse data 10x10 or so?}
350 In script~\ref{code: gravity1} we use the statement
351 \begin{verbatim}
352 source0=NetCdfData(NetCdfData.GRAVITY, 'QLDWest_grav.nc', altitude=0.)
353 \end{verbatim}
354 to load the gravity data stored in \examplefile{QLDWest_grav.nc} in the \netcdf format.
355 Within the script the data set is now available under the name \verb|source0|.
356 We need to link the data set to the \verb|DomainBuilder| using
357 \begin{verbatim}
358 dom.addSource(source0)
359 \end{verbatim}
360 At the time the domain for the inversion is built the \verb|DomainBuilder|
361 will use the information about origin, extent and spacing along with the other
362 options provided to build an appropriate domain.
363 As at this point a flat Earth is assumed geographic coordinates used to
364 represent data in the input file are mapped to a (local) Cartesian coordinate
365 system. This is achieved by projecting the geographic coordinates into the
366 \emph{Universal Transverse Mercator} (UTM) coordinate system\footnote{See e.g.
367 \url{http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system}}.
369 It is possible to add several, possibly overlapping data sets to the same
370 domain:
371 \begin{verbatim}
372 source0=NetCdfData(NetCdfData.GRAVITY, 'data0.nc', altitude=0.)
373 source1=NetCdfData(NetCdfData.GRAVITY, 'data1.nc', altitude=0.)
374 dom.addSource(source0)
375 dom.addSource(source1)
376 \end{verbatim}
377 However, due to the coordinate transformation all data sets must be located in
378 the same UTM zone. If a single dataset crosses UTM zones only the zone of the
379 central longitude is used when projecting.
380 For example, if one data set lies mostly in zone 51 but contains areas of zone
381 52, it is transformed using zone 51.
382 In this case more data from zone 51 can be added, but not from any other zone.
384 It is important to keep an eye on the complexity of the inversion.
385 A good measure is the total number of cells being used.
386 Assume we have given a data set on a $20 \times 20$ grid and we add lateral
387 padding of, say, $20 \%$ to each side of the data, the lateral number of cells
388 becomes $(20\cdot 1.4)\times (20\cdot 1.4)=1.4^2\cdot 20^2\approx 2\cdot 10^2=800$.
389 If we use $20$ cells in the vertical direction we end up with a total number
390 of $800 \times 20 = 16,000$ cells.
391 This size can be easily handled by a modern desktop PC.
392 If we increase the grid size of the data to $40 \times 40$ points and use $40$
393 cells in the vertical extent we get a total of $(2\cdot 40^2)\cdot 40=128,000$
394 cells, a problem size which is considerably larger but can still be handled by
395 a desktop computer.
396 Taking this one step further, if the amount of data is increased to
397 $200\times 200$ points and we use $200$ cells in the vertical extent the
398 domain will contain $16,000,000$ ($16$ million) cells.
399 This scenario requires a computer with enough memory and (a) fast processor(s)
400 to run the inversion.
401 This estimate of complexity growth applies to the case where the increase
402 of data grid size is driven by an increase of resolution where it is
403 recommended to increase the vertical resolution in synch with the lateral
404 resolution. Note that if more than one data set is used the target resolution
405 will be the resolution of the finest data set.
406 In other cases the expansion of the region of interest drives an increase of
407 data grid size and the increase of total number of cells is less dramatic as
408 the vertical number of cells can remain constant while keeping a balanced
409 resolution in vertical and lateral direction.
411 \section{Setting up the Inversion and Running it}
412 We are now at step three of script~\ref{code: gravity1} in which the actual
413 inversion is set up.
414 First we create an empty inversion under the name \verb|inv|:
415 \begin{verbatim}
416 inv=GravityInversion()
417 \end{verbatim}
418 As indicated by the name we can use \verb|inv| to perform an inversion of
419 gravity data\footnote{\verb|GravityInversion| is a driver with a simplified
420 interface which is provided for convenience. See Part~\ref{part2} for more
421 details on how to write inversion scripts with more general functionality, e.g.
422 constraints.}. The inversion is an iterative process which sequentially
423 calculates updates to the density distribution in an attempt to improve the
424 match of the gravity field produced by the density distribution with the data.
425 Termination of the iteration is controlled by the tolerance which is set by
426 the user:
427 \begin{verbatim}
428 inv.setSolverTolerance(1e-4)
429 \end{verbatim}
430 Here we set the tolerance to $10^{-4}$, i.e. the iteration is terminated if
431 the maximum density correction is less than or equal to $10^{-4}$ relative to
432 the maximum value of estimated density anomaly.
433 In case the iteration does not converge a maximum number of iteration steps is
434 set:
435 \begin{verbatim}
436 inv.setSolverMaxIterations(50)
437 \end{verbatim}
438 If the maximum number of iteration steps (here $50$) is reached the iteration
439 process is aborted and an error message is printed.
440 In this case you can try to rerun the inversion with a larger value for the
441 maximum number of iteration steps.
442 If even for a very large number of iteration steps no convergence is achieved,
443 this is a strong indicator that the inversion has not been set up properly.
445 The statement
446 \begin{verbatim}
447 inv.setup(dom)
448 \end{verbatim}
449 links the inversion with the domain and the data.
450 At this step -- as we are solving a gravity inversion problem -- only
451 gravitational data attached to the domain builder \verb|dom| are considered.
452 Internally a cost function $J$ is created which is minimized during the
453 inversion iteration.
454 It is a combination of a measure of the data misfit of the gravity field from
455 the given density distribution and a measure of the smoothness of the density
456 distribution.
457 The latter is often called the regularization term.
458 By default the gradient of density is used as the regularization term, see
459 also Section~\ref{SEC:P1:GRAV:REMARK:REG}.
460 Obviously, the result of the inversion is sensitive to the weighting between
461 the misfit and the regularization.
462 This trade-off factor $\mu$ for the misfit function is set by the following
463 statement:
464 \begin{verbatim}
465 inv.getCostFunction().setTradeOffFactorsModels(0.1)
466 \end{verbatim}
467 Here we set $\mu=0.1$. The statement \verb|inv.setup| must appear in the
468 script before setting the trade-off factor.
469 A small value for the trade-off factor $\mu$ will give more emphasis to the
470 regularization component and create a smoother density distribution.
471 A large value of the trade-off factor $\mu$ will emphasize the misfit more
472 and typically creates a better fit to the data and a rougher density
473 distribution.
474 It is important to keep in mind that the regularization reduces noise in the
475 date and in fact gives the problem a unique solution.
476 Consequently, the trade-off factor $\mu$ may not be chosen too large in order
477 control the noise on the solution and ensure convergence in the iteration
478 process.
480 We can now actually run the inversion:
481 \begin{verbatim}
482 rho = inv.run()
483 \end{verbatim}
484 The answer as calculated during the inversion is returned and can be accessed
485 under the name \verb|rho|.
486 As pointed out earlier the iteration process may fail in which case the
487 execution of the script is aborted with an error message.
489 \section{Taking a Look}
490 In the final step of script~\ref{code: gravity1} the calculated density
491 distribution is written to an external file.
492 A popular file format used by several visualization packages such as
493 \VisIt~\cite{VISIT} and \mayavi~\cite{MAYAVI} is the \VTK file format.
494 The result of the inversion which has been named \verb|rho| can be written to
495 the file \file{result.vtu} by adding the statement
496 \begin{verbatim}
497 saveVTK("result.vtu", density=rho)
498 \end{verbatim}
499 at the end of script.
500 The inversion solution is tagged with the name \verb|density| in the result
501 file, however any other name for the tag could be used.
502 As the format is text-based (as opposed to binary) \VTK files tend to be very
503 large and take compute time to create, in particular when it comes to large
504 numbers of cells ($>10^6$).
505 For large problems it is more efficient to use the \SILO file format~\cite{SILO}.
506 \SILO files tend to be smaller and are faster generated and read.
507 It is the preferred format to import results into the visualization program
508 \VisIt~\cite{VISIT} which is particularly suited for the visualization of
509 large data sets.
510 Inversion results can directly exported into \SILO files using the statement
511 \begin{verbatim}
512 saveSilo("result.silo", density=rho)
513 \end{verbatim}
514 replacing the \verb|saveVTK(...)| statement.
515 Similar to \VTK files the result \verb|rho| is tagged with the name
516 \verb|density| so it can be identified in the visualization program.
518 \begin{figure}
519 \begin{center}
520 \subfigure[$\mu=0.1$]{%
521 \label{FIG:P1:GRAV:10 MU01}
522 \includegraphics[width=0.45\textwidth]{QLDGravContourMu01.png}
523 }%
524 \subfigure[$\mu=1.$]{%
525 \label{FIG:P1:GRAV:10 MU1}
526 \includegraphics[width=0.45\textwidth]{QLDGravContourMu1.png}
527 }\\ % ------- End of the first row ----------------------%
528 \subfigure[$\mu=10.$]{%
529 \label{FIG:P1:GRAV:10 MU10}
530 \includegraphics[width=0.45\textwidth]{QLDGravContourMu10.png}
531 }%
532 \subfigure[$\mu=100.$]{%
533 \label{FIG:P1:GRAV:10 MU100}
534 \includegraphics[width=0.45\textwidth]{QLDGravContourMu100.png}
535 }\\ % ------- End of the second row ----------------------%
536 \subfigure[$\mu=1000.$]{%
537 \label{FIG:P1:GRAV:10 MU1000}
538 \includegraphics[width=0.45\textwidth]{QLDGravContourMu1000.png}
539 }%
540 \end{center}
541 \caption{3-D contour plots of gravity inversion results with data from
542 Figure~\ref{FIG:P1:GRAV:0} for various values of the model trade-off
543 factor $\mu$. Visualization has been performed in \VisIt.
544 \AZADEH{check images.}}
545 \label{FIG:P1:GRAV:10}
546 \end{figure}
548 \begin{figure}
549 \begin{center}
550 \subfigure[$\mu=0.1$]{%
551 \label{FIG:P1:GRAV:11 MU01}
552 \includegraphics[width=0.45\textwidth]{QLDGravDepthMu01.png}
553 }%
554 \subfigure[$\mu=1.$]{%
555 \label{FIG:P1:GRAV:11 MU1}
556 \includegraphics[width=0.45\textwidth]{QLDGravDepthMu1.png}
557 }\\ % ------- End of the first row ----------------------%
558 \subfigure[$\mu=10.$]{%
559 \label{FIG:P1:GRAV:11 MU10}
560 \includegraphics[width=0.45\textwidth]{QLDGravDepthMu10.png}
561 }%
562 \subfigure[$\mu=100.$]{%
563 \label{FIG:P1:GRAV:11 MU100}
564 \includegraphics[width=0.45\textwidth]{QLDGravDepthMu100.png}
565 }\\ % ------- End of the second row ----------------------%
566 \subfigure[$\mu=1000.$]{%
567 \label{FIG:P1:GRAV:11 MU1000}
568 \includegraphics[width=0.45\textwidth]{QLDGravDepthMu1000.png}
569 }%
570 \end{center}
571 \caption{3-D slice plots of gravity inversion results with data from
572 Figure~\ref{FIG:P1:GRAV:0} for various values of the model trade-off
573 factor $\mu$. Visualization has been performed \VisIt.
574 \AZADEH{check images.}}
575 \label{FIG:P1:GRAV:11}
576 \end{figure}
578 Figures~\ref{FIG:P1:GRAV:10} and~\ref{FIG:P1:GRAV:11} show two different
579 styles of visualization generated in \VisIt using the result of the inversion
580 of the gravity anomalies shown in Figure~\ref{FIG:P1:GRAV:0}.
581 The inversions have been performed with different values for the model
582 trade-off factor $\mu$.
583 The visualization shows clearly the smoothing effect of lower values for the
584 trade-off factors.
585 For larger values of the trade-off factor the density distribution becomes
586 rougher showing larger details.
587 Computational costs are significantly higher for larger trade-off factors.
588 Moreover, noise in the data has an higher impact on the result.
589 Typically several runs are required to adjust the value for the trade-off
590 factor to the datasets used.
592 For some analysis tools it is useful to process the results in form of Comma
593 Separated Values (\CSV)~\footnote{see
594 \url{http://en.wikipedia.org/wiki/Comma-separated_values}.}.
595 Such a file can be created using the statement
596 \begin{verbatim}
597 saveDataCSV("result.csv", x=rho.getFunctionSpace().getX(), density=rho)
598 \end{verbatim}
599 in the script.
600 This will create a \file{result.csv} with columns separated by a comma.
601 Each row contains the value of the density distribution and the three
602 coordinates of the corresponding location in the domain.
603 There is a header specifying the meaning of the corresponding column.
604 \LG{Add example output CSV}.
605 Notice that rows are not written in a particular order and therefore, if
606 necessary, the user has to apply appropriate sorting of the rows.
607 Columns are written in alphabetic order of their corresponding tag names.
608 For the interested reader: the statement \verb|rho.getFunctionSpace()| returns
609 the type used to store the density data \verb|rho|.
610 The \verb|getX()| method returns the coordinates of the sampling points used
611 for the particular type of representation, see~\cite{ESCRIPT} for details.
614 \section{Remarks}
616 \subsection{Data With Holes}\label{SEC:P1:GRAV:REMARK:DATAHOLES}
617 \CIHAN{ADD some remarks holes in the data, example?}
620 \subsection{ER Mapper Raster File}\label{SEC:P1:GRAV:REMARK:ERMAPPER}
621 \CIHAN{ADD an example how to use ER Mapper raster files}
623 \subsection{Regularization Term}\label{SEC:P1:GRAV:REMARK:REG}
624 The \verb|GravityInversion| class supports the following form for the
625 regularization:
626 \begin{equation}
627 \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
628 \end{equation}
629 where the integral is calculated across the entire domain.
630 $\rho$ represents the density distribution where $\rho_{,0}$ $\rho_{,1}$ and
631 $\rho_{,2}$ are the spatial derivatives of $\rho$ with respect to the two
632 lateral and the vertical direction, respectively.
633 $w^{(0)}$, $w^{(1)}_0$, $w^{(1)}_1$ and $w^{(1)}_2$ are weighting
634 factors\footnote{A more general form, e.g. spatially variable values for the
635 weighting factors, is supported, see Part~\ref{part2}}.
636 By default these are $w^{(0)}=0$, $w^{(1)}_0=w^{(1)}_1=w^{(1)}_2=1$.
637 Other weighting factors can be set in the inversion set-up.
638 For instance to set $w^{(0)}=10$, $w^{(1)}_0=w^{(1)}_1=0$ and $w^{(1)}_2=100$
639 use the statement:
640 \begin{verbatim}
641 inv.setup(dom, w0=10, w1=[0,0,100])
642 \end{verbatim}
643 It is pointed out that the weighting factors are rescaled in order to improve
644 numerical stability. Therefore the relative size of the weighting factors is
645 relevant and using
646 \begin{verbatim}
647 inv.setup(dom, w0=0.1, w1=[0,0,1])
648 \end{verbatim}
649 would lead to the same regularization as the statement above.
651 %
652 % \section{Reference}
653 %
654 % There are three examples for 2D and 3D gravity inversions with artificial input data.
655 % 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.
656 %
657 % 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:
658 %
659 % \begin{verbatim}
660 % n_cells_in_data=100
661 % depth_offset=0.*U.km
662 % l_data = 100 * U.km
663 % l_pad=40*U.km
664 % THICKNESS=20.*U.km
665 % l_air=6*U.km
666 % \end{verbatim}
667 %
668 % The others assumptions comes with each examples.
669 % \begin{enumerate}
670 % \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})
671 %
672 % \begin{verbatim}
673 % n_cells_in_data=100
674 % n_humbs_h= 1
675 % n_humbs_v=1
676 % mu=100
677 % \end{verbatim}
678 %
679 % \begin{figure}
680 % \centering
681 % \includegraphics[width=\textwidth]{grav2D1.png}
682 % \caption{2D density model up) reference down) result}
683 % \label{fig:gravity2D1}
684 % \end{figure}
685 %
686 %
687 % \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})
688 %
689 % \begin{verbatim}
690 % n_cells_in_data=100
691 % n_humbs_h= 3
692 % n_humbs_v=1
693 % mu=100
694 % \end{verbatim}
695 %
696 % \begin{figure}
697 % \centering
698 % \includegraphics[width=\textwidth]{grav2D3.png}
699 % \caption{2D density model up) reference down) result}
700 % \label{fig:gravity2D3}
701 % \end{figure}
702 %
703 % \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})
704 %
705 % \begin{verbatim}
706 % n_cells_in_data=50
707 % n_humbs_h= 1
708 % n_humbs_v=1
709 % mu=10
710 % \end{verbatim}
711 %
712 % \begin{figure}
713 % \centering
714 % \includegraphics[width=\textwidth]{density3D-ref.png}
715 % \caption{3D density model of reference as synthetic data}
716 % \label{fig:gravity3D}
717 % \end{figure}
718 %
719 % \begin{figure}
720 % \centering
721 % \includegraphics[width=\textwidth]{gravity3D.png}
722 % \caption{3D density model of result}
723 % \label{fig:gravity3D1}
724 % \end{figure}
725 % \end{enumerate}

  ViewVC Help
Powered by ViewVC 1.1.26