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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4207 - (show annotations)
Mon Feb 18 03:17:11 2013 UTC (6 years, 6 months ago) by azadeh
File MIME type: application/x-tex
File size: 33693 byte(s)
adding bar for gravity image

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 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
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 line-by-line.
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(1e-4)
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 run-escript grav.py
138 \end{verbatim}
139 We are running \file{grav.py} through the \escript start-up 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/escript-finley}.
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 set-up 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{2-D domain set-up 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 set-up for a two-dimensional domain
194 (see also Figure~\ref{fig:domainBuilder} for 3-D).
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 2-D 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 metre-kilogram-second 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 \verb|U.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 \verb|U| 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 pre-processed. In particular, corrections for
291 \begin{itemize}
292 \item free-air, 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 \verb|NX| $\times$ \verb|NY| where \verb|NX| and \verb|NY| define the number
321 of grid lines in the longitude (\verb|X|) and the latitude (\verb|Y|)
322 direction, respectively.
323 The grid is spanned from an origin with spacing \verb|DELTA_X| and
324 \verb|DELTA_Y| in the longitude and the latitude direction, respectively.
325 The gravity data for all grid points need to be given as an \verb|NX|
326 $\times$ \verb|NY| array.
327 If available errors can be associated with the gravity data.
328 The values are given as an \verb|NX| $\times$ \verb|NY| 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 \verb|source0|.
355 We need to link the data set to the \verb|DomainBuilder| using
356 \begin{verbatim}
357 dom.addSource(source0)
358 \end{verbatim}
359 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 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 \verb|inv|:
414 \begin{verbatim}
415 inv=GravityInversion()
416 \end{verbatim}
417 As indicated by the name we can use \verb|inv| to perform an inversion of
418 gravity data\footnote{\verb|GravityInversion| 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(1e-4)
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 \verb|dom| 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 trade-off 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 \verb|inv.setup| must appear in the
467 script before setting the trade-off factor.
468 A small value for the trade-off factor $\mu$ will give more emphasis to the
469 regularization component and create a smoother density distribution.
470 A large value of the trade-off 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 trade-off 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 \verb|rho|.
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 \verb|rho| 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 \verb|density| in the result
500 file, however any other name for the tag could be used.
501 As the format is text-based (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 \verb|saveVTK(...)| statement.
514 Similar to \VTK files the result \verb|rho| is tagged with the name
515 \verb|density| 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{3-D contour plots of gravity inversion results with data from
537 Figure~\ref{FIG:P1:GRAV:0} for various values of the model trade-off
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{3-D slice plots of gravity inversion results with data from
563 Figure~\ref{FIG:P1:GRAV:0} for various values of the model trade-off
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 trade-off factor $\mu$.
574 The visualization shows clearly the smoothing effect of lower values for the
575 trade-off factors.
576 For larger values of the trade-off factor the density distribution becomes
577 rougher showing larger details.
578 Computational costs are significantly higher for larger trade-off 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 trade-off
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/Comma-separated_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 \verb|rho.getFunctionSpace()| returns
600 the type used to store the density data \verb|rho|.
601 The \verb|getX()| 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 \verb|GravityInversion| 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 set-up.
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]{density3D-ref.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

  ViewVC Help
Powered by ViewVC 1.1.26