ARTICLE IN PRESS
Computers & Geosciences 31 (2005) 989–999 www.elsevier.com/locate/cageo
Differential reduction to the pole$ G.R.J. Coopera,, D.R. Cowanb a
School of Geosciences, University of the Witwatersrand, Johannesburg, South Africa b Cowan Geodata Services, 12 Edna Road, Dalkeith, Western Australia
Received 2 December 2004; received in revised form 31 January 2005; accepted 22 February 2005
Abstract Due to the dipolar nature of the geomagnetic field, magnetic anomalies located anywhere other than at the magnetic poles are asymmetric even when the magnetic source distribution is symmetrical. This complicates interpretation. Pole reduction (RTP) takes the anomaly, as measured at any latitude, and transforms it into that which would have been measured if the body had been laid at the magnetic pole i.e. the area where the field inclination is vertical and the anomalies from symmetrical bodies are symmetrical. The algorithm is usually applied in the frequency domain which has the advantage of being computationally fast, but the disadvantage of restricting the application of the algorithm to regions possessed of constant geomagnetic inclination and declination. Various frequency domain algorithms have been devised to solve this problem, but all have problems of one form or another. A simple algorithm for RTP is suggested here that employs a Taylor series expansion in the space domain (other work has used the Taylor series expansion approach in the frequency domain). The algorithm is demonstrated on synthetic data and on aeromagnetic data from the Northern Territory, Australia. r 2005 Elsevier Ltd. All rights reserved. Keywords: Aeromagnetics; Pseudogravity; Geophysics
1. Introduction Pole reduction is an operator which takes magnetic anomalies and changes their asymmetric form to the symmetric form which would have been observed had the causative magnetic bodies lain at the magnetic poles. The frequency domain operator is (Baranov, 1957) Aðu; vÞ A0 ðu; vÞ ¼ (1) ðsin y þ i cos y sinðf þ aÞÞ2 $ Code on server at http://www.iamg.org/CGEditor/index.htm. Corresponding author. Tel.: +27 117176608; fax: +27 117176579. E-mail addresses:
[email protected] (G.R.J. Cooper),
[email protected] (D.R. Cowan). URL: http://www.wits.ac.za/science/geophysics/gc.htm.
where A(u,v) is the amplitude at frequencies ðu; vÞ, y and f are the geomagnetic inclination and declination, respectively, and a is tan1(v/u). The method has several (well-known) problems when implemented in this fashion. It is unstable at low magnetic latitudes, it gives incorrect results if the causative magnetic bodies possess unknown remanent magnetisation, and lastly, because of the frequency domain implementation of the algorithm, y and f must remain constant throughout the area of application of the filter. There are several approaches that can be taken to solve this latter problem. Lu et al. (2003) used a parallel computer to reduce the dataset to the pole nxm times, where the dataset contains nxm datapoints. The inclination and declination could be different at each grid point as required, and only the response centred on the current point was retained from each RTP operation.
0098-3004/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2005.02.005
ARTICLE IN PRESS 990
G.R.J. Cooper, D.R. Cowan / Computers & Geosciences 31 (2005) 989–999
The method is effective but requires considerable computer power. The equivalent layer method can also be used to apply RTP when the field parameters vary (Von Frese et al., 1981; Silva, 1986). In this case, the inversion stage that determines the layer susceptibilities uses sources with inclination and declination that vary over the dataset. Then the forward model is calculated with all layer dipole inclinations set to 901. The computational effort required to perform the inversion is again the main problem with the method. If the variations of the field parameters are small, then they can be considered as perturbations about the average field values of the region. Arkani-Hamed (1988) used this idea in the frequency domain, allowing the crustal magnetisation to vary continuously over a plane, using an iterative algorithm. However, as pointed out by Swain (2000), the method is rarely used in practice because of the large data storage requirements of the algorithm and because it is unstable at low magnetic latitudes. Additionally, the iterative nature of the algorithm makes its computational requirements yet more demanding.
2. Differential RTP as a space domain perturbation An alternative method is to apply the perturbations in the space domain rather than the frequency domain. So, using a Taylor series expansion qRTP q2 RTP þ 0:5Dinc2 qinc qinc2 2 qRTP q RTP þ Ddec þ ð2Þ þ 0:5Ddec2 qdec qdec2
RTPvar ¼ RTPmean þ Dinc
RTPmean is the dataset reduced to the pole using the average field inclination and declination of the area. Dinc is the difference between the inclination at a given point and the average inclination, and Ddec is computed similarly. The derivatives are computed in the space domain by differencing. Fig. 1a shows a synthetic dataset data consisting of three dipoles, and the field inclination varies from 301 to 901. The standard RTP image (shown in Fig. 1b) used an inclination of 601 for the whole dataset, whereas the differential RTP image (Fig. 1c) used the space domain Taylor series expansion of equation 1. Figs. 1d–f plot NS profiles through the anomaly maxima of the data in Figs. 1a–c. The anomaly shapes in Fig. 1e are almost as asymmetrical as those of the original data, whereas those in Fig. 1f show a noticeable improvement. The question arises as to how many terms are required in Eq. 2 to achieve sufficient accuracy. This obviously depends on how ‘‘sufficient accuracy’’ is defined, but an estimate of the error at each point can be obtained from the remainder term of the Taylor series after n terms,
given by Wolfram:1 Rn ¼
ðx x0 Þnþ1 f nþ1 x , ðn þ 1Þ!
(3)
where x lies between x and x0. In terms of the pole reduction problem this becomes Rn ¼
Dincnþ1 Ddecnþ1 nþ1 nþ1 RTPMEAN RTPMEAN inc þ dec , ðn þ 1Þ! ðn þ 1Þ! (4)
where inc* lies between 0 and Dinc, and dec* lies between 0 and Ddec. In practice, the number of terms necessary depends on the manner in which the geomagnetic field varies over the area of interest. In the synthetic dipole model of Fig. 1, the inclination varied in a linear manner from south to north and almost no improvement (compared to the theoretical result) was achieved by using more than the first two terms of Eq. 2. If the geomagnetic field varied in a more complex manner then more terms would be required. Hence, prior examination of the variation of the field parameters (perhaps by polynomial surface fitting) is useful in determining the number of terms to use.
3. Differential pseudogravity Once the differential reduction to the pole dataset has been computed it is relatively simple to convert to pseudogravity by vertical integration and the application of a scale factor. Pseudogravity converts the magnetic field into the gravity field that would be observed if the magnetisation distribution were to be replaced with an identical density distribution (Blakely, 1995, p. 344). It is a useful technique for the interpretation of major magneto-tectonic provinces as it simplifies anomaly patterns and focuses on large-scale features rather than local detail. In the wavenumber domain, the nth order vertical integral is given by (Blakely, 1995, p. 326) A0 ðkÞ ¼ AðkÞjkjn ,
(5)
where AðkÞ is the amplitude at a wavenumber k, k ¼ 2p=l and l is wavelength. To convert to pseudogravity A0 (k) is multiplied by Gr=M (assuming that r=M is a constant, where M is magnetisation and r is density) before performing an inverse Fourier transform.
4. Application to aeromagnetic data from the Northern Territories, Australia The new differential reduction to the pole algorithm was applied to the regional aeromagnetic grid of the 1 Wolfram, S., 2005. http://mathworld.wolfram.com/TaylorSeries.html.
ARTICLE IN PRESS G.R.J. Cooper, D.R. Cowan / Computers & Geosciences 31 (2005) 989–999
991
Fig. 1. Differential RTP of synthetic data: (a) Synthetic aeromagnetic dataset consisting of three dipoles, with field inclination varying from 301 to 901, and declination being constant throughout; (b) the dataset in Fig. 1a after RTP using a fixed inclination of 601; (c) the dataset in Fig 1a after RTP using space domain Taylor series expansion of Eq. (1); (d) NS profile through centre of data in Fig. 1a; (e) NS profile through centre of data in Fig. 1b and (f) NS profile through centre of data in Fig. 1c.
Northern Territory, Australia. The magnetic field inclination ranges from 571 in the south to 331 in the north. The simplified geology of the Northern Territory is shown in Fig. 2a. Proterozoic Orogens of the North Australian Craton and the Central Australian Mobile Belts form the basement. The North Australian Craton is interpreted as complex accreted terranes (Myers et al., 1996) with Archaean granite–gneiss and granitoids forming the basement to Palaeoproterozoic geosynclinal strata which were deformed and metamorphosed during the Barramundi Orogeny (1860–1850 Ma). The North Australian Craton Proterozoic Orogens include the Pine Creek Orogen, Tanami Region, Tennant Region, Murphy and Arnhem Inliers. The Proterozoic Orogens of the Central Australian Mobile Belts include part of the Arunta Region and the Mesoproterozoic Musgrave Block. The Arunta region is a complex basement inlier which has been active from around 1880 Ma to the Alice Springs Orogeny (400–300 Ma). Largely unmetamorphosed Palaeo- to Mesoproterozoic Platform sediments unconformably overlie the Proterozoic Orogens including the McArthur, South Nicholson and Victoria–Birrindudu Basins and the Ashburton and Davenport Provinces of the Tennant Inlier. Extensive Neoproterozoic–Palaeozoic intracratonic basins cover much of the Northern Territory. Neoproterozoic sediments in the Georgina, Wiso, Amadeus and Ngalia Basins are overlain by Cambrian Basalts and Cambrian and Ordovician
clastics and carbonate sediments. Finally there are Mesozoic Basins including the Eromanga-Pedirka and Dunmarra Basins. Previous RTP of this aeromagnetic dataset was performed by Geoscience Australia by breaking the dataset into a series of 300 300 km overlapping panels, using the mean geomagnetic inclination and declination of each panel, retaining only the central 100 100 km portion, and then stitching the panels together (Milligan, 2004, Pers. Comm.). According to Milligan, this worked well for the vertical magnetic gradient but was less successful for total magnetic intensity. As well as being inelegant, this method has the disadvantage that long wavelength features (comparable in size or greater than the panel size) are not correctly reduced to the pole. The new algorithm adjusts the field inclination progressively rather than in fixed increments. A total magnetic intensity colour shade image for the Northern Territory is shown in Fig. 2b along with a standard reduction to the pole image (Fig. 2c) using a fixed field inclination corresponding to the centre of the area. Comparison of the total magnetic intensity data in Fig. 2b with the differential reduction to the pole in Fig. 2d shows that reduction to the pole has produced significant improvement in image clarity by simplifying magnetic patterns, especially in the north of the image. At these lower magnetic latitudes, magnetic anomaly shapes are more difficult to relate to source locations and reduction to the pole makes qualitative magnetic
ARTICLE IN PRESS 992
G.R.J. Cooper, D.R. Cowan / Computers & Geosciences 31 (2005) 989–999
Fig. 2. Application to regional aeromagnetic data: (a) Simplified geology of Northern Territory, Australia; (b) Total intensity aeromagnetic data of Northern Territory, Australia; (c) Data from Fig. 2b after application of standard RTP; (d) Data from Fig. 2b after application of differential RTP and (e) Difference between data in Fig. 2d and c.
ARTICLE IN PRESS G.R.J. Cooper, D.R. Cowan / Computers & Geosciences 31 (2005) 989–999
993
Fig. 2. (Continued)
interpretation simpler. The image contrasts the highfrequency active magnetics of the Proterozoic Orogens with the subdued, long-wavelength response of the basins and allows the interpreter to follow basement magnetic trends below cover.
Comparison of standard (Fig. 2c) and differential RTP (Fig. 2d) shows differences both in magnetic relief and magnetic texture. Fig. 2e shows the difference between the differential reduction to the pole and standard reduction to the pole using the field inclination
ARTICLE IN PRESS 994
G.R.J. Cooper, D.R. Cowan / Computers & Geosciences 31 (2005) 989–999
Fig. 2. (Continued)
at the centre of the area. Although the differences are not large in amplitude (due to the maximum inclination difference from the mean value being only 121) significant changes in anomaly shapes are obvious and the differential reduction to the pole is more coherent and has correct DC levels.
Fig. 3 shows the differential pseudogravity produced by vertical integration of the differential reduction to the pole data. The image provides a much clearer picture of major basement magneto-tectonic provinces, providing a balanced view of major structural units with major basement magnetic anomaly trends below the Georgina
ARTICLE IN PRESS G.R.J. Cooper, D.R. Cowan / Computers & Geosciences 31 (2005) 989–999
995
Fig. 2. (Continued)
and McArthur Basins being much clearer than in Fig. 2d. To illustrate the use of the reduction to the pole data we have selected a subarea shown in Fig. 2d. Fig. 4a is a simplified geological map of the subarea which covers part of the large Georgina Basin and the adjacent
Arunta Province and Davenport Fold Belt basement. Fig. 4b shows a differential reduction to the pole image of the subarea and Fig. 4c the geology image blended with the reduction to the pole data. The images highlight the abrupt contrast in anomaly wavelengths between outcropping Arunta Province basement and the deep
ARTICLE IN PRESS 996
G.R.J. Cooper, D.R. Cowan / Computers & Geosciences 31 (2005) 989–999
Fig. 2. (Continued)
basement of the Altjawarra Domain basement which is covered by 1–2 km of Neoproterozoic and Palaeozoic sediments (Tompkins et al., 2003). The Altjawarra Domain comprises a central, geophysically distinct zone known as the Altjawarra Cratonic Nucleus, and adjacent, related crustal elements. The Altjawarra Cratonic Nucleus extends for approximately 145 km in an east–west orientation and 125 km in a north–south
direction, incorporating large areas of low-amplitude and rugose magnetic texture joints, fractures and minor dykes, all of which are characteristic of a granitoid– gneiss terrane. These areas are separated by narrow, high-amplitude linear and curvilinear anomalies, interpreted to represent complexly folded and faulted mafic and ultramafic dominated terranes (or greenstone belts) traversed by prominent northeast to north-northeast,
ARTICLE IN PRESS G.R.J. Cooper, D.R. Cowan / Computers & Geosciences 31 (2005) 989–999
997
Fig. 3. Differential pseudogravity.
west-northwest and east-west structures. The ratio of interpreted granitoids and gneisses to greenstones approximates 80:20, similar to Archaean terranes elsewhere. The conspicuous magnetic high (or in Fig. 4b) is the approximately circular Ooratippra anomaly, with magnetic amplitudes up to 2000 nT and a coincident
gravity high. There are a large number of plug-like magnetic anomalies varying in size from a few hundred metres to several kilometres. The modeled depth to the magnetic source indicates the platform cover varies from 1 to 2 km. The central region of the Altjawarra Domain is characterized by an E–W oriented magnetic fabric and
ARTICLE IN PRESS 998
G.R.J. Cooper, D.R. Cowan / Computers & Geosciences 31 (2005) 989–999
Fig. 4. Application to detail aeromagnetic data: (a) Simplified geology of Georgina Basin and adjacent areas; (b) Differential reduction to pole for Georgina Basin and adjacent areas and (c) Blend of reduction to pole and simplified geology.
ARTICLE IN PRESS G.R.J. Cooper, D.R. Cowan / Computers & Geosciences 31 (2005) 989–999
is bounded to the north and south by palaeo-rift features of approximately E–W orientation. The Davenport Fold Belt shows a complex pattern of folding and faulting with source depth close to surface. Collision of the Altjawarra Domain with the Davenport Province to the northwest is likely to have occurred during consolidation of the North Australian Craton, which was complete by approximately 1880 Ma. A postcollision northeast trending suture between these blocks is postulated, possibly resulting in intrusion of Mesoproterozoic granitoids, including the Ooratippra basement feature.
4. Conclusions A simple algorithm for the pole reduction of magnetic datasets with variable geomagnetic inclination and declination was presented. The algorithm is computationally simple and gave good results both on synthetic data and on aeromagnetic data from Australia.
Acknowledgements We thank the Northern Territory Geological Survey for providing the original total magnetic intensity grid and the geological overview maps.
999
References Arkani-Hamed, J., 1988. Differential reduction to the pole of regional magnetic anomalies. Geophysics 53 (12), 1592–1600. Baranov, V., 1957. A new method for interpretation of aeromagnetic maps:pseudo-gravimetric anomalies. Geophysics 22 (2), 359–383. Blakely, R.J., 1995. Potential Theory in Gravity and Magnetic Applications. Cambridge University Press, New York, NY 435pp. Lu, R.S., Mariano, J., Willen, D.E., 2003. Differential reduction of magnetic anomalies to the pole on a massively parallel computer. Geophysics 68 (6), 1945–1951. Myers, J.S., Shaw, R.D., Tyler, I.M., 1996. Tectonic evolution of Proterozoic Australia. Tectonics 15, 1431–1446. Silva, J.C.B., 1986. Reduction to the pole as an inverse problem and its application to low latitude anomalies. Geophysics 51 (2), 369–382. Swain, C.J., 2000. Reduction to the pole of regional magnetic data with variable field direction, and its stabilisation at low inclinations. Exploration Geophysics 31, 78–83. Tompkins, L.A., Taylor, W., Cowan, D.R., 2003. Diamond prospectivity of the Altjawarra Craton, Australia. Eighth International Kimberlite Conference, Vancouver Island, BC, Canada (Extended Abstract). Von Frese, R.R.B., Hinze, W.J., Braile, L.W., 1981. Spherical earth gravity and magnetic analysis by equivalent point source inversion. Earth and Planetary Science Letters 53, 69–83.