manuscript.doc

<h1 id="euler-deconvolution-of-potential-field-data">Euler deconvolution of potential field data</h1>
<p>Leonardo Uieda, Vanderlei C. Oliveira Jr, and Valéria C. F. Barbosa. Observatório Nacional, Rio de Janeiro, RJ, Brazil.</p>
<p>In this tutorial we'll talk about a widely used method of interpretation for potential-field data called Euler deconvolution. Our goal is to demonstrate its usefulness and, most importantly, call attention to some pitfalls encountered in the interpretation of the results. The code and synthetic data required to reproduce our results and figures can be found in the accompanying IPython notebooks (<strong>ipython.org/notebook</strong>) at <strong>dx.doi.org/10.6084/m9.figshare.923450</strong> or <strong>github.com/pinga-lab/paper-tle-euler-tutorial</strong>. The notebooks also expand the analysis presented here. We encourage you to download the data and try it on your software of choice. For this tutorial we'll use the implementation in the open-source Python package Fatiando a Terra (<strong>fatiando.org</strong>).</p>
<h3 id="brief-overview">Brief overview</h3>
<p>Euler deconvolution was first developed by Thompson <span class="citation">(1982)</span> and later extended by Reid et al. <span class="citation">(1990)</span>. Since then, it has been adapted and improved upon by Keating <span class="citation">(1998)</span>, Mushayandebvu et al. <span class="citation">(2004)</span>, and many others. This popularity is largely due to its great simplicity of implementation and use, making it the tool of choice for a quick initial interpretation.</p>
<p>In many cases, maps of gravity and magnetic data (and transformations thereof) provide good constraints on the horizontal location of an anomaly source. Euler deconvolution adds an extra dimension to the interpretation. It estimates a set of (x, y, z) points that, ideally, fall inside the source of the anomaly. Euler deconvolution requires the x-, y-, and z-derivatives of the data and a parameter called the <em>structural index</em> (SI). The SI is an integer number that is related to the homogeneity of the potential field and varies for different fields and source types <span class="citation">(<em>Stavrev and Reid</em>, 2007)</span>. For example, in the case of total field magnetic anomaly data, a sphere is represented by an SI of 3, whereas a dike is represented by an SI of 1. There are methods that can estimate the SI and we refer the reader to Barbosa et al. <span class="citation">(1999)</span> and Melo et al. <span class="citation">(2013)</span> for more information.</p>
<h3 id="an-example-with-synthetic-data">An example with synthetic data</h3>
<p>The best way to fully grasp the usefulness and limitations of any geophysical method is to apply it in a controlled setting using synthetic data. To make things more interesting, we'll apply the Euler deconvolution on data generated by a model that differs slightly from ideal sources (Figure 1).</p>
<p>First, we need to calculate the 3 directional derivatives of the data. For gridded data, this can be done with the Fourier transform using module <code>fourier</code> from Fatiando a Terra. Assuming that the data has been loaded into Numpy arrays <code>x</code>, <code>y</code>, <code>z</code> with the respective x-, y-, z-coordinates of each data point and the array <code>data</code> with our magnetic data (see the accompanying notebook):</p>
<pre><code>from fatiando.gravmag import fourier
from fatiando import utils
# Need data in SI units
data_si = utils.nt2si(data)
# shape = (ny, nx) number of points in y and x
dx = fourier.derivx(x, y, data_si, shape)
dy = fourier.derivy(x, y, data_si, shape)
dz = fourier.derivz(x, y, data_si, shape)</code></pre>
<p>Now that we have calculated the derivatives, we have to create a solver to run the deconvolution. Fatiando provides the <code>Classic</code> solver which uses the whole dataset to estimate a single point (source location). The usual strategy for multiple sources, introduced by Reid et al. <span class="citation">(1990)</span>, is to use a <em>moving window</em> over the data. A point is estimated for each window using the data that fall inside it, producing a large number of solutions. The challenge then becomes separating the good solutions from the spurious ones. FitzGerald et al. <span class="citation">(2004)</span> provide an overview of different selection criteria. Fatiando provides the <code>MovingWindow</code> solver that can run a given Euler solver on a moving window (see the accompanying notebook for more details on the selection criterion used):</p>
<pre><code>from fatiando.gravmag.euler import Classic, MovingWindow
classic = Classic(x, y, z, data_si, dx, dy, dz, struct_index)
solver = MovingWindow(classic, windows=(50, 50), size=(1000, 1000))
solver.fit()</code></pre>
<p>Calling <code>solver.fit()</code> runs the deconvolution using the specified <span class="math">50 × 50</span> windows of size 1 km. The estimated points are stored in the <code>solver.estimate_</code> variable.</p>
<p>For comparison, we ran the deconvolution using various different combinations of structural index and window size. Figure 2 shows scatter plots of the solutions (the typical way Euler solutions are presented). Notice how changing the window size influences the spread of the solutions and can affect our interpretation. For example, for a 1 km window and structural index 3, we cannot tell apart the two model bodies in the left. Furthermore, increasing the structural index will increase the depths of the solutions.</p>
<p>We can view these solutions in 3D using the Fatiando module <code>myv</code>, a wrapper for the Mayavi visualization software (<strong>code.enthought.com/projects/mayavi</strong>):</p>
<pre><code>from fatiando.vis import myv
myv.figure()
myv.points(solver.estimate_)
myv.axes(myv.outline(model_bounds))
myv.show()</code></pre>
<p>This will create something similar to Figure 3. Notice that the Euler solutions do fall approximately inside the sources. However, it is curious how the solutions make arch-like structures that have no relation to the actual forms of the sources. In fact, Silva et al. <span class="citation">(2001)</span> have shown that Euler solutions are biased and that the bias depends on the moving-window scheme and selection criterion used.</p>
<h3 id="last-few-words-of-caution">Last few words of caution</h3>
<p>We repeat here for emphasis: Euler deconvolution solutions <strong>do not</strong> represent the three-dimensional outline of the sources. They indicate, at most, an approximate source location. It is crucial to remember that Euler solutions are subject to non-uniqueness and bias just like any other geophysical inverse problem. Tests using synthetic data are indispensable to assert the plausibility of our interpretations.</p>
<p><em>Corresponding author: leouieda@gmail.com</em></p>
<div class="figure">
<img src="fig/data-model-low.png" alt="Our model (a) and synthetic total field anomaly data (b). The model simulates a dike (blue), a sill (green), and an intrusive body (red). The data is corrupted with 5 nT pseudo-random Gaussian noise." /><p class="caption">Our model (a) and synthetic total field anomaly data (b). The model simulates a dike (blue), a sill (green), and an intrusive body (red). The data is corrupted with 5 nT pseudo-random Gaussian noise.</p>
</div>
<div class="figure">
<img src="fig/euler-solutions-low.png" alt="Euler deconvolution solutions for varying structural index (SI) and moving window size." /><p class="caption">Euler deconvolution solutions for varying structural index (SI) and moving window size.</p>
</div>
<div class="figure">
<img src="fig/euler-solutions-3d-composite-low.png" alt="3D view of the Euler deconvolution solutions (black dots) for a moving window of 3 km and structural index 3." /><p class="caption">3D view of the Euler deconvolution solutions (black dots) for a moving window of 3 km and structural index 3.</p>
</div>
<h2 id="references">References</h2>
<p>Barbosa, V. C. F., J. B. C. Silva, and W. E. Medeiros (1999), Stability analysis and improvement of structural index estimation in Euler deconvolution, <em>Geophysics</em>, <em>64</em>, 48–60, doi:10.1190/1.1444529.</p>
<p>FitzGerald, D., A. Reid, and P. McInerney (2004), New discrimination techniques for Euler deconvolution, <em>Computers &amp; Geosciences</em>, <em>30</em>, 461–469, doi:10.1016/j.cageo.2004.03.006.</p>
<p>Keating, P. B. (1998), Weighted Euler deconvolution of gravity data, <em>Geophysics</em>, <em>63</em>, 1595–1603.</p>
<p>Melo, F. F., V. C. F. Barbosa, L. Uieda, V. C. Oliveira Jr, and J. B. C. Silva (2013), Estimating the nature and the horizontal and vertical positions of 3D magnetic sources using Euler deconvolution, <em>Geophysics</em>, <em>78</em>, 87, doi:10.1190/geo2012-0515.1.</p>
<p>Mushayandebvu, M. F., V. Lesur, A. B. Reid, and J. D. Fairhead (2004), Grid Euler deconvolution with constraints for 2D structures, <em>Geophysics</em>, <em>69</em>, 489–496, doi:10.1190/1.1707069.</p>
<p>Reid, A. B., J. M. Allsop, H. Granser, A. J. Millett, and I. W. Somerton (1990), Magnetic interpretation in three dimensions using Euler deconvolution, <em>Geophysics</em>, <em>55</em>, 80–91, doi:10.1190/1.1442774.</p>
<p>Silva, J. B. C., V. C. F. Barbosa, and W. E. Medeiros (2001), Scattering, symmetry, and bias analysis of source-position estimates in Euler deconvolution and its practical implications, <em>Geophysics</em>, <em>66</em>, 1149–1156, doi:10.1190/1.1487062.</p>
<p>Stavrev, P., and A. Reid (2007), Degrees of homogeneity of potential fields and structural indices of Euler deconvolution, <em>Geophysics</em>, <em>72</em>, 1, doi:10.1190/1.2400010.</p>
<p>Thompson, D. T. (1982), EULDPH: A new technique for making computer-assisted depth estimates from magnetic data, <em>Geophysics</em>, <em>47</em>, 31–37, doi:10.1190/1.1441278.</p>