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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4185 - (show annotations)
Fri Feb 8 09:06:21 2013 UTC (6 years, 8 months ago) by gross
File MIME type: application/x-tex
File size: 24515 byte(s)
modifications to cookbook
1 \chapter{Gravity Inversion}\label{Chp:cook:gravity inversion}
2
3 \section{Introduction}
4
5 In this part the documentation we give an introduction on how to us the \downunder module using the
6 inversion driver functions to perform the inversion of gravity and magnetic data. The driver functions
7 enables geologists and geophysicists apply the \downunder module quickly and in easy way
8 with out a detailed knowledge of the theory of inversion or programming skills. However, users
9 how are interested in specializing or extending the inversion capacity are referred to the second part~\ref{part2}
10 of this manual. It is beyond the intention of this manual to give a in-detailed introduction to
11 geophysical inversion in particular to the appropriate preprocessing of data sets. We refer to~\cite{REF1, REF2, REF3}.
12
13
14 The \downunder function described here are designed for calculate estimations for the
15 3-D distribution of density and/or susceptibility from 2-D gravity and magnetic data measured in ground
16 or airborne surveys. This process is generally called inversion of geophysical data.
17 Following the standard assumption it is assumed that the
18 data are measured as perturbation of an expected gravity and/or magnetic
19 field of the Earth. In this context measured gravity and magnetic data are in fact describing anomalies in
20 gravity and magnetic field. As a consequence the inversion process
21 provides corrections to an average density (typically $2670 kg/m^3$) and susceptibility (typically $0$).
22 So in the following we will always assume that given data for anomalies are given and therefore
23 not in all cases explicitly use the terms such gravity anomalies or density corrections but just use the terms
24 gravity and density.
25
26 In this chapter we will give a detailed introduction into usage of the driver functions
27 for inversion of gravity data. In the following chapters~\ref{Chp:cook:magnetic inversion} and~\ref{Chp:cook:joint inversion}
28 we will discuss the inversion of magnetic and the joint inversion of gravity and magnetic data using
29 \downunder. As the same principles as for gravity data apply the presentation for these problem classes is kept short and
30 users interested in magnetic data only should still work through this chapter on gravity data.
31
32 To run the examples discussed you need to have \escript (version 3.3.1 or newer) installed on your computer.
33 Moreover, if you want to visualize result, you must have access to plotting software which is able to process
34 \VTK input files, e.g. \mayavi and \VisIt. As \mayavi can easily be obtained and installed for most platforms
35 the tutorial is based on \mayavi as visualization. However, it pointed out
36 that \VisIt is the preferred visualization tool for \escript as it can deal with very large data sets more effectively.
37
38
39
40 \begin{figure}
41 \centering
42 \includegraphics[width=\textwidth]{density10.png}
43 \caption{3D contour plot of density distribution of inversion of data \file{bouguer_anomaly.nc}. Colors
44 represent values of density where high values are represented by red and low values are represented by blue.}
45 \label{FIG:P1:GRAV:1}
46 \end{figure}
47
48 \section{How does it work?}
49 The execution of the inversion is controlled by script which is in essence is a text file. This can be edited using any
50 text editor. The script contains a serious of statements which are executed from an interpreter which
51 is an executable program reading the text file and executing the statements line-by-line. In the case
52 of \downunder the interpreter is \python. In order to be able to process the statements in each line of the script
53 certain rules (called syntax) need to be obeyed. There is large number of online tutorials for \python
54 available\footnote{e.g. \url{http://www.tutorialspoint.com/python} and \url{http://doc.pyschools.com}}. We also
55 refer to the \escript cook book \cite{ESCRIPTCOOKBOOK} and users guide \cite{ESCRIPT} which is in particular useful for users which who to dive deeper
56 \downunder. For this part of the manual no \python knowledge is required but it is recommended that
57 users acquire some basic knowledge on \python as they progressing in their work with \downunder.
58
59 The following script~\ref{code: garvity1}\footnote{The script is similar to the \examplefile{gravity_netcdf.py} script
60 find in the \escript example file directory.} is a simple example to run an inversion for gravity data:
61 \begin{pyc}
62 \
63 \begin{verbatim*}
64 # Header:
65 from esys.downunder import *
66 from esys.weipa import *
67 import esys.escript.unitsSI as U
68
69 # Step 1: set up domain
70 dom=DomainBuilder()
71 dom.setVerticalExtents(depth= 40.*U.km, air_layer=6.*U.km, num_cells=25)
72 dom.setFractionalPadding(pad_x=0.2, pad_y=0.2)
73 dom.fixDensityBelow(depth= 40. * U.km)
74
75 # Step 2: read gravity data
76 source0=NetCdfData(NetCdfData.GRAVITY, 'bouguer_anomaly.nc',altitude=0.)
77 dom.addSource(source0)
78
79 # Step 3: set up inversion
80 inv=GravityInversion()
81 inv.setSolverTolerance(1e-4)
82 inv.setSolverMaxIterations(50)
83 inv.setup(dom)
84
85 # Step 4: run inversion
86 inv.getCostFunction().setTradeOffFactorsModels(10.)
87 rho = inv.run()
88
89 # Step 5: write reconstructed density to file
90 saveVTK("result.vtk", density=rho)
91 \end{verbatim*}\label{code: garvity1}
92 \end{pyc}
93 The result, in this
94 case the density distribution, is written to an internal file for further processing. You can copy and past the text of the
95 script into a file of any name, let's say for further reference we use the file name \file{grav.py}. It is recommendable
96 to use the extention \file{py} to identify the file as \python script. We will discuss the statements
97 of the script later in this chapter.
98
99 Obviously the inversion needs to be fed with some gravity data. In this case
100 the data are loaded from the file \file{bouguer_anomaly.nc} which available in the \escript example file directory. The data are given
101 in the \netcdf data format for gridded data. After you have copied this file into the directory in which you have
102 saved the script \file{grav.py} you can run the program using the command line
103 \begin{verbatim}
104 run-escript grav.py
105 \end{verbatim}
106 We are running \file{grav.py} through the \escript start-up command as \escript is used as a back end for the inversion
107 algorithm\footnote{Please see the \escript users guide~\cite{ESCRIPT} on how to run
108 your script in parallel using threading and/or MPI.}. Obviously it is assumed that you have an
109 installation of \escript available on your computer, see \url{https://launchpad.net/escript-finley}.
110
111 After the execution has successfully been completed you will find the result file \file{result.vtk} in the directory
112 from where you have started the execution of the script. the file is using the \VTK format and can be important
113 easily into many visualization tools. One option is the \mayavi package which is available on most platforms. You can invoke the
114 the visualization using the commands
115 \begin{verbatim}
116 mayavi2 -d result.vtu -m SurfaceMap
117 \end{verbatim}
118 from the command line. Figure~\ref{FIG:P1:GRAV:1} shows the result using contour plotting\footnote{Plots
119 have been created with \VisIt.}.
120
121 Let us have a closer look at the script\footnote{In \python lines starting with `\#` are commends and are not processed further.}: Besides the header section one can separate the script into five steps:
122 \begin{enumerate}
123 \item set up domain on which the inversion problem is solved
124 \item load the data
125 \item set-up the inversion problem
126 \item run the inversion
127 \item further processing of the result, here we write the reconstructed density distribution to a file.
128 \end{enumerate}
129 In the following we will discuss the steps of the scripts in more details. Before we do this it is pointed out that
130 the header section which, by following \python standards, makes all the package we need available to access within the script. At this
131 point we will not discuss this in more details but emphasize that the header section is a vital part of the script. It is is required
132 in each \downunder inversion script and should not be altered accept if additional modules are needed.
133
134 \begin{figure}
135 \centering
136 \includegraphics[width=\textwidth]{dom2D.png}
137 \caption{2D domain set-up for gravity inversion.}
138 \label{FIG:P1:GRAV:2}
139 \end{figure}
140
141 \section{How to create the inversion Domain}
142 First step in the script~\ref{code: garvity1} is the definition of the domain over which the inversion is performed.
143 We think of the domain as a block with orthogonal, plain faces. Figure~\ref{FIG:P1:GRAV:2} shows the set-up for
144 a two dimensional domain (see also Figure~\ref{fig:domainBuilder} for 3D). The lateral coordinates along the
145 surface of the Earth are denoted by $x$ and $y$ (on $x$-direction is shown). The $z$ direction defines the
146 vertical direction where the part above the surface has positive coordinates and the
147 subsurface negative coordinates. The height of the section above the surface which is assumed to be filled to with air
148 needs to be set by the user. The inversion assumes that the density in the section is known to be zero. The density
149 below the surface is unknown and is calculated through the inversion. The user needs to specify the depth below
150 the surface in which the density is to be calculated.
151 The lateral extension of the domain is defined by the data sets fed into the inversion. It is chosen large enough
152 to cover all data sets. In order to reduce the impact of boundary a padding zone around the data sets can be introduced.
153
154 \begin{figure}
155 \centering
156 \includegraphics[width=\textwidth]{dom2DFEM.png}
157 \caption{Cell distribution and boundary conditions for 2D domain.}
158 \label{FIG:P1:GRAV:2}
159 \end{figure}
160
161 The reconstruction of the gravity field from an estimated density distribution is a the key component of the inversion process.
162 To do this \downunder uses the finite element method (FEM) to do this. We need to introduce a grid over the domain, see Figure~\ref{FIG:P1:GRAV:2}.
163 The number of vertical cells is set by the user while the number of vertical cells is derived from the grid spacing of the
164 gravity data set(s). It is assumed that gravity field data given are constant across a cell.
165 To be able to reconstruct the gravity field some assumptions on the values of the gravity field on the domain boundary
166 have to made. \downunder assumes that on all faces the lateral gravity field component equals
167 zero. No assumptions on the
168 horizontal components are made\footnote{It is assumed that the gravity potential equals zero on the top and bottom
169 surface, see Section~\ref{sec:forward gravity} for details}\footnote{Most inversion codes are using Green's function over an unbounded domain to reconstruct the gravity field. This approach
170 makes the assumption the gravity field (or potential) is converging to zero when moving away from the region of interest. The
171 boundary conditions used here are stronger in the sense that the lateral gravity component is enforced to zero in
172 a defined distance of the region of interest but weaker in the sense that no constraint on the horizontal component is applied.}.
173
174 In script~\ref{code: garvity1} the statement
175 \begin{verbatim}
176 dom=DomainBuilder()
177 \end{verbatim}
178 creates same form of a factory to build a domain.
179 We define the features of the domain we would like to create:
180 \begin{verbatim}
181 dom.setVerticalExtents(depth= 40.*U.km, air_layer=6.*U.km, num_cells=25)
182 dom.setFractionalPadding(pad_x=0.2, pad_y=0.2)
183 \end{verbatim}
184 Here we define the depth of the domain to $40 km$, the thickness of the air layer above the surface to $6km$ and
185 the number of vertical cells to $25$. We also introduce a lateral padding of $20 \%$ of the expansion of
186 the gravity data on each side of the data and in both lateral directions.
187
188 In some cases it can be appropriate to assume that the density below a certain depth is
189 zero\footnote{As we are in fact calculate density corrections this means that the density is assumed to be
190 the average density.}. The statement
191 \begin{verbatim}
192 dom.fixDensityBelow(depth= 40. * U.km)
193 \end{verbatim}
194 introduces this constrain. As in the case discussed here the depth for zero density is not less than the
195 depth of the domain no constraint at depth is applied to the density.
196
197 ADD U
198
199 \section{Getting the Data in}
200
201 \section{Set-up Inversion and run the show}
202
203 \section{Take a look}
204
205 \section{Remarks}
206
207 =====
208
209
210 Before we do this a few words about the header section of the script:
211 \begin{verbatim*}
212 from esys.downunder import *
213 from esys.weipa import *
214 import esys.escript.unitsSI as U
215 \end{verbatim*}
216 The \verb|import| statement makes certain functions directly available for usage in the script. After \downunder
217 has been imported in the first import statement all functions in the \downunder package cane directly been used,
218 eg. we can call \verb|DomainBuilder| to create a domain for the inversion. Similar
219 the second import statement makes all the \weipa functions ready for usage in the script.
220 Here we will in particular use the \verb|saveVTK| function to write the inversion result to
221 an external file\footnote{See also \escript users guide~\cite{ESCRIPT}.}.
222
223
224 For instance we can
225 directly call which is part Here we
226 load all the functions of \downunder which is the inversion tool kit and \weipa which supports the output
227 of the inversion into various file formats. Notice that both package use \escript so this is also imported behind the scenes.
228 The third \verb|import| statement makes the
229 $40 km$
230 \verb|40 * U.km|
231
232
233
234
235
236
237 ==============================
238
239 will not make a difference between gravity and gravity anomalies
240 and density and density corrections. We will loosely use the terms density and gravity when referring to corrections and anomalies.
241
242 When using \downunder driver function it is also assumed that the necessary corrections (e.g. for height and topography) are applied
243 before data are fed into the inversion.
244
245 To solve the inversion problem
246 the earth is modeled using a large number of points - so called nodes - in a rectangular grid.
247 Density and susceptibility are estimated on these nodes.
248 There distribution is obtained by minimizing a cost function to fit the observed data and to achieve
249 a desirable smoothness of the distribution of density and susceptibility.
250
251
252
253 Density and susceptibility as well as the resulting gravity and magnetic field
254 are calculated as a distribution across a hexahedron which contains a subsurface region of the Earth filled
255 with solid material of unknown density and susceptibility and a surface area filled with air
256 in which the gravity and magnetic field are measured and density and susceptibility are in fact known.
257
258 The inversion problem is solved in an iterative manner
259 where an estimation of the density and susceptibility distribution
260 is improved through a series of update steps in order to fit the result with the input data.
261 This iteration process is terminated after a stopping criteria is met, e.g. when the estimated
262 termination error reaches a given tolerance.
263
264
265
266
267
268 \begin{figure}
269 \centering
270 \includegraphics[width=\textwidth]{density11.png}
271 \caption{Depth image across previous 3D gravity inversion which presents discrepancy in density. Diversity in density are detected with colors and contours.}
272 \end{figure}
273
274
275
276
277
278 \section{Gravity Data}
279 Theoretically weighs depends on the force of gravity at
280 that position and the force of gravity varies with elevation, rock densities, latitude and topography.
281 Mass, however, does not depend on gravity but is a constitutive quantity throughout the earth. So the
282 spring stretch in mass suspension is related to the gravity force. In addition with a constant mass,
283 difference in spring stretch illustrate the changes in acceleration of gravity. The principal of gravity exploration
284 is based on the topography of the basement, thickness of the sedimentary section with various density, porosity of
285 the layer and elevation.
286
287 The amount of $g$ at sea level is about 980 $kg/s^2$ or 980,000 $mgal$. Gravity acceleration
288 is measured in two types. The first corresponds to specify the absolute magnitude of gravity
289 at any place and the second refers to the alteration in gravity from one place to another.
290 In gravity study variation of this value which is caused by underground structure, is plotted as a residual gravity.
291 Closed variation in residual gravity indicates discrepancy in subsurface geological structure.
292
293 The Gravity data are taken from onshore and offshore observation recorded at many gravity stations with high precision
294 of determination in elevation and position (latitude and longitude). All raw reading gravity observations require
295 to process with many corrections.
296
297 The sun and moon gravitational forces make curvature in Earth's shape. These tide effects change figure of oceans,
298 atmosphere and even solid body of the Earth, which impress gravity measurement and it is necessary to compensate it,
299 which vary with location, date and time of the day.
300
301 The surface of the Earth is lumpy on land and water. However for Geophysical and Geological study, a smooth closed
302 surface is assumed. The main one is a spheroid flattened at the poles which is called ellipsoid. The new
303 data are used to defined a best-fitting obtained ellipsoid. The second suggestion is geoid which
304 is really mathematical convenience. There is a uniform mass between gravity stations and ellipsoid, that's
305 effect must be removed with corrections.
306
307 Also the level of topography for hilly and valley measurements is important. The gravity amount which is
308 made up by that equal the mass of hill or valley must be added as a terrain correction to have a measurement on a level surface.
309
310 Because gravity descend towards the poles the latitude correction must be added to the observed gravity.
311
312 Free-air correction that must be added to observation, ignores the effects of material between the measurement
313 and reference level which is positive for above sea-level station and is negative for station below sea-level.
314
315 The gravitational acceleration or Bouguer correction is calculated for known thickness (between measurements
316 station and ellipsoid) and density (depends on local rock) which must be subtracted from the measurements
317 gravity if station is above sea-level. and the station is below sea-level this must be added.
318
319 For this adjustment first of all Tidal correction apply based on tidal table or calculated tidal effect
320 for given time and location of gravity data. The second one is terrain correction. Then Latitude, Free air and
321 Bouguer correction have applied. The standard Bouguer density is 2670 $kg/m^3$.
322
323 (in the input file, bouguer anomaly is used as residual gravity which have to be gridded.)
324
325
326
327
328 \section{Input File}
329
330 For starting up the inversion 2 files are needed. Each of the
331 two files contains a series of parameters which must appear in the correct order, as described in
332 the next paragraph. Each parameter is marked by a keyword, which is followed either on the same line
333 or the next line by one or more parameter values.
334
335 The first file is contained the gravity anomaly, including number of points, the accurate location
336 (latitude and longitude) of the observed position and the value of the anomaly after all gravity corrections.
337 The main part in collecting the data is spacing of the observation point which means that all data must be gridded and
338 its spacing affected the simulation resolution wisely.
339
340 The second is run_gravity with py extension which is related to the codes though
341 it needs some constraints to have a good results in inversion.
342
343 A small part of sample of run_gravity:
344
345 \begin{verbatim}
346 mu=100
347 THICKNESS=20.*U.km
348 DATASET='QLD_west.nc'
349 PAD_X = 0.2
350 PAD_Y = 0.2
351 l_air = 6. * U.km
352 n_cells_v = 25
353 \end{verbatim}
354
355 Run_gravity file consist many options to implement which control how inversion is performed such as padding area, depth , MU factor,\ldots.
356
357 \begin{description}
358 \item[MAX\_ITER]
359 Specifying maximum iteration depends on model, the area which have been selected to have an inversion in it, your hard capacity and the level of users satisfaction. However the best result were built with 200 iterations. In addition all steps of inversion could be traced and the suited one selected.
360
361 \item[l_pad or PAD\_X, PAD\_Y] To implement the boundary constraints in this file, padding area in $x$ and $y$ direction should be determined. In a rectangular area, same padding for both direction is preferable. If directional area will be fixed with elongation in one orientation it does not matter to change the padding area. In addition for 2D inversion padding just add in one direction.
362
363 \item[DATASET] In run_gravity file just the location and the name of the data file have to be fixed. Both source of data and run_gravity files must be in a folder.
364
365 \item[THICKNESS] Depth of the model should be assigned to have dipper or shallower inversion also it is assigned to the layer where shows inversion in. The parameter values must be real numbers, and they represent depths in km.
366
367 \item[l_air] Length of air is height of the model above the sea level.
368
369 \item[n_cells_v] The last important part of the inversion property is the number of elements in data of the model which shows the finer or coarser cells in the model so its delimitation have affected on resolution of the inversion.
370
371 \item[mu]It is defined in accordance with the noise of data and it has a wide range to select from 0.0001 to 100 for each inversion.
372
373 \end{description}
374
375 \section{Output File}
376
377 At the end of each inversion iteration, package produces a new inversion file with 'gravin.silo' name which is stored in the file. This silo file shows the inversion result which does not have the number of iterations stage.
378
379 In terminal indicates the specifications of inversion iterations. In this descriptions of the paths of all inversion during stages of modeling are cleared. The format of the file will be described here as it is designed to be used directly for analysis or debugging.
380
381 After final iteration the silo file, the result of inversion, is visible with some software which shows silo file format. It illuminate density distribution in the area (not in the padding) which create the gravity input data.
382
383 \section{Reference}
384
385 There are three examples for 2D and 3D gravity inversions with artificial input data.
386 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.
387
388 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:
389
390 \begin{verbatim}
391 n_cells_in_data=100
392 depth_offset=0.*U.km
393 l_data = 100 * U.km
394 l_pad=40*U.km
395 THICKNESS=20.*U.km
396 l_air=6*U.km
397 \end{verbatim}
398
399 The others assumptions comes with each examples.
400 \begin{enumerate}
401 \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})
402
403 \begin{verbatim}
404 n_cells_in_data=100
405 n_humbs_h= 1
406 n_humbs_v=1
407 mu=100
408 \end{verbatim}
409
410 \begin{figure}
411 \centering
412 \includegraphics[width=\textwidth]{grav2D1.png}
413 \caption{2D density model up) reference down) result}
414 \label{fig:gravity2D1}
415 \end{figure}
416
417
418 \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})
419
420 \begin{verbatim}
421 n_cells_in_data=100
422 n_humbs_h= 3
423 n_humbs_v=1
424 mu=100
425 \end{verbatim}
426
427 \begin{figure}
428 \centering
429 \includegraphics[width=\textwidth]{grav2D3.png}
430 \caption{2D density model up) reference down) result}
431 \label{fig:gravity2D3}
432 \end{figure}
433
434 \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})
435
436 \begin{verbatim}
437 n_cells_in_data=50
438 n_humbs_h= 1
439 n_humbs_v=1
440 mu=10
441 \end{verbatim}
442
443 \begin{figure}
444 \centering
445 \includegraphics[width=\textwidth]{density3D-ref.png}
446 \caption{3D density model of reference as synthetic data}
447 \label{fig:gravity3D}
448 \end{figure}
449
450 \begin{figure}
451 \centering
452 \includegraphics[width=\textwidth]{gravity3D.png}
453 \caption{3D density model of result}
454 \label{fig:gravity3D1}
455 \end{figure}
456 \end{enumerate}

  ViewVC Help
Powered by ViewVC 1.1.26