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

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}