Earth Sci. Res. J. Vol. 11, No. 2 (December 2007): 108-116 AccurAte grAvity AnomAly interpolAtion: A cAse-study in cAmeroon, centrAl AfricA J. Kamguia1, c. t. tabod2, J. m. tadjou1, e. manguelle-dicoum2, r. nouayou2 and l. H. Kande1 1 National Institute of Cartography (NIC) – Cameroon, PO Box 157, Yaounde – Cameroon 2 Department of Physics, Faculty of Science, University of Yaounde – Cameroon Corresponding author: Joseph Kamguia: kjerryfr@yahoo.fr; j.kamguia@geomaticsystem.com; 0023799440426 ABstrAct Many treatments in geodesy and geophysics require regularly gridded gravity anomalies. The gridding of gravity data needs interpolation. For the predicted data to be accurate, the smoothest type of gravity anomaly should be used along with the most indicated prediction method. This paper presents the comparison of various prediction methods applied on different types of gravity anomalies and considering the relative geological complexity of the study area. Many algorithms are tested and the suitability of each type of anomaly and each prediction method discussed in a case-study in Cameroon (Central Africa), using a set of 43,000 gravity data points to determine the must accurate prediction technique. Key words: Prediction, Gravity anomaly, Residual free-air anomalies, Global geopotential model, Complex geology, Geophysics. resumen Muchos tratamientos en geodesia y geofísica requieren que las anomalías gravimétricas se encuentren en grillas regulares. Esta disposición de los datos de gravedad requiere interpolación. Para que los datos predichos sean exactos, la anomalía gravimétrica más fina debe ser utilizada junto con el método más indicado para predicción. Este artículo presenta la comparación de varios métodos de predicción aplicados a diversos tipos de anomalías gravimétricas y considerando la relativa complejidad geológica del área del estudio. Se probaron una gran cantidad de algoritmos y la conveniencia de cada tipo de de anomalía y método de predicción fueron discutidos en una aplicación en Camerún (África central), usando un conjunto de 43.000 datos gravimétricos para determinar la técnica de predicción más adecuada.. 108 eArtH sciences reseArcH JournAl Manuscript received July 15 2007. Accepted for publication December 5 2007. Palabras claves: Predicción, Anomalía Gravimétrivca, Anomalías Residuales Libres de Aire, Modelo Geopotencial Global, Geología Compleja, Geofísica. introduction Interpolation is a process widely used in earth science. It estimates the value of a parameter at a point from neighbourings. The process may be extrapolation in some areas (Heiskanen and Moritz, 1967). In Geodesy, gravity anomalies are interpolated when computing geoid undulations or the quasi-geoid height anomalies. In Geophysics, prediction is used to compute regular grids of gravity anomalies required for presentation purposes or further treatments. The accuracy of the predicted gravity anomalies depends mainly on their smoothness and on the method used. It also depends the data density, their geographical distribution and the grid intervals. In general, the smoother the gravity anomalies, the more accurate the predicted values. Many interpolation methods are now proposed (El Abbas et al., 1990), with some algorithms available in open sources such as the Generic Mapping Tools or GMT (Wessel and Smith, 1995). Comparisons have already been made between surface fitting algorithms, in order to check their reliability in predicting Bouguer gravity anomalies (El Abbas et al., 1990). However, many types of gravity anomalies are routinely used in geodesy and in geophysics. It may sometimes be difficult to determine which type is a priori suited for prediction. A constant crustal density value generally used in different computations might be too far from the reality below the earth topographic surface in complex geological areas. Residual anomalies that are deduced from free-air (FAA), simple Bouguer (BA), complete Bouguer (CBA) and isostatic (IA) anomalies might also be smoother than the original anomalies in some areas. These anomalies are now easily computed using the spherical harmonic coefficients of a global geopotential model (GGM) (Kamguia et al., 2007). This paper presents the results of comparison of various types of gravity anomaly predicted with some methods commonly and in different areas of Cameroon. These anomalies were computed from a homogeneous database containing 43,000 data points. The aim is to give some hints as to the most appropriate gravity anomaly prediction strategy according to the density of the gravity net, their geographical distribution and the relative complexity of the geology of the study area. To achieve this, four areas with different geological characteristics are chosen. In these areas of different extension and gravity net density, six prediction algorithms are tested: the Minimum Curvature Splines in Tension (Smith and Wessel, 1990); the Least Square Polynomial Fitting (two variants); Krigging (Krige, 1978) and the inverse distance to a power method (two variants). 2. Gravity anomaly prediction 2.1. Characteristics of gravity anomalies Different types of gravity anomaly are used in earth sciences. Gravity data recorded on land must be adjusted for elevation above or below sea level to yield the free-air anomalies (FAA). They are nearly equivalent to what would be observed if all the topographic masses were condensed onto this reference surface (Blakely, 1996). Simple and complete Bouguer gravity anomalies (BA and CBA) are obtained by completely removing the masses that exist between the level of observation and sea level. They are in theory smoother and then more ACCuRATE GRAVITy ANoMALy INTERPoLATIoN: A CASE-STuDy IN CAMERooN, CENTRAL AFRICA 109 Kamguia et al. ESRJ Vol. 11, No. 2. December 2007 110 suitable for prediction than the FAA (Heiskanen and Moritz, 1967). However, in complex geological contexts, the gravity field may have special characteristics. Therefore, BA and CBA may be rougher than FAA. The isostatic reduction regularizes the earth’s crust according to a model of isostasy. Depending on the complexity of the geology of the study area, these anomalies may be less or more suited in prediction processes than FAA, BA and CBA. From spherical harmonic coefficients of global geopotential models (GGM), the long wavelength components of FAA, BA, CBA and IA are easily computed. After subtraction from respective total field gravity anomalies, the residual anomalies obtained may be smoother in some areas. The long wavelength anomalies ( )GGMg∆ are computed from the harmonic coefficients of the GGM using equation (1) and subtracted from the total field gravity anomalies )( g∆ (Heiskanen and Moritz, 1967): (1) where GM is the geocentric gravitational constant; ),,( qlr the spherical coordinates of the computation point; g the normal gravity on the reference ellipsoid; a the equatorial radius of the earth; (sin )n mP q the fully normalized associated Legendre functions for degree n and order m ; mn C∆ and mn S∆ the normalized EGM-GGM harmonic coefficients, reduced for the even zonal harmonic for the ellipsoid and complete to degree and order 360max =N . 2.2. Prediction of gravity anomalies and statistical analysis The gravity is a regionalized field, since it varies from place to place in a continuous manner but with no possibility of associating a specified mathematical function. The most accurate gravity anomaly prediction technique is the one that maintains the broad and fine features of the original gravity data processed, without introducing undue distortions. Therefore, the predicted and the measured values must at the same point should be mathematically and/or physically related. To achieve this, the criterion of Crain and Bhattacharyya (1967) is important. Data point Grid node with a reference point of gridded data (i) Input data in gridding algorithms (ii) Gridded cells waiting for gridded data at nodes   (i) (ii) figure 1. Choice of reference points of gridded data   max 2 2 0 ( 1) cos sin (sin ) nN n GGM n m n m n m n m GM a g n C m S m P r r                To apply the criterion, some data points must be selected to serve as reference points (Fig. 1). A reference point is a data point with coordinates nearly coinciding with those of a grid node. The predicted and known values are therefore comparable in a simple statistical analysis of their differences. 3. A case-study in Cameroon (Central Africa) 3.1. Characteristics of the gravity field in Cameroon Free-air gravity anomalies are smoother than BA and CBA in Cameroon, as is indicated in Table 1. In this table, the standard deviation (STD) of the BA and CBA are nearly 10 mGal greater than the FAA’s. This means that the gravity field in Cameroon possesses special features. one can therefore conclude that: ● A constant density )67.2( =r used in the Bouguer reduction does not effectively suppress the effect of the topography above the reference surface in Cameroon; ● The disturbing geological masses responsible of the observed gravity anomalies cannot be easily determined here; ● The geology of Cameroon is more complex than it looks like. In this context, interpolation of gravity anomalies may be subject to many errors if the more indicated method and especially the suitable type of gravity anomaly are not used. Table 1: Statistical comparison of Bouguer and free-air anomalies in Cameroon (unit: mGal) Four test areas ),,,( 4321 ZZZZ (Fig. 2) were selected in order to precise the special features of the gravity field and its behaviour in prediction processes. Their selection is based on the differences in geological characteristics, the density and distribution of the data points, and the relative roughness of the topography relief of the area. The principal characteristics of these areas are summarized in Table 2. Table 2: Characteristics of the four test areas The first test area 1Z is located in the northern sedimentary region of Cameroon. The area is characterized by a long wavelength relatively positive Bouguer anomaly attributed by Adighije (1981), and Fairhead and Okereke (1987) to an uplift of the Moho. The topographic relief has a mean orthometric height of about 1000 m. Gravity data are more denser in this test area. The second test area 2Z is situated on the North eastern extension of the Cameroon volcanic line (CVL). It is characterised by large negative Bouguer anomalies oriented E-W from western Cameroon to western C. A. R. These are attributed to a light mass anomaly structure in the upper mantle beneath the Adamawa uplift (Fairhead and Okereke, 1987; Poudjom et al., 1992). The mean elevation is greater than 1000 m. The gravity data density is the smallest in this test area. 3Z is located on the continental part of the CVL. This line is associated with a series of NE- SW trending gravity lows not well defined by the reconnaissance survey of Collignon (1968). The mean height of the area is nearly 2000 m, with a maximum of 4010 m for Mount-Cameroon. The gravity data density is intermediate to those of 1Z and 2Z . The area 4Z is situated in southern Cameroon. The mean height is lower than in the 111 Anomaly Max. Min. Mean STD Variability Free-air 315,5 –85,4 –0,03 25,0 Smoother Simple Bouguer 253,4 –125,5 –26,5 34,4 Rougher Complete Bouguer 254,2 –121,4 –25,5 33,6 Rougher Number ofTest area points used control points Data mean density 1Z  EE  1613  ;  NN  118  2179 36 1 point per 25 km 2Z  EE  1513  ;  NN  85  672 16 1 point per 2100 km 3Z  EE  119  ;  NN  74  1465 9 1 point per 250 km 4Z  EE  1412  ;  NN  52  987 16 1 point per 2100 km ACCuRATE GRAVITy ANoMALy INTERPoLATIoN: A CASE-STuDy IN CAMERooN, CENTRAL AFRICA 112 previous test areas and the gravity data density is comparable to what is observed in 3Z . The distributions of data in each test area are shown in figures 3 to 6. The reference points are well distributed in 1Z . The numbers of these points in the test areas are enough for an informative analysis of predicted gravity anomalies. Table 3 shows the statistics of the anomalies to be interpolated in the test areas. From the analysis based on the STD, one can see that the total field Bouguer and isostatic anomalies are not always the smoothest. The gravity field has special characteristics in Cameroon and its surroundings. The high frequency components of this field, created by the topography relief, are not entirely removed after the Bouguer and isostatic reductions. By smoothing FAA, BA, CBA and IA with subtraction of the longer wavelength components computed using equation (1), four residual types of anomaly specified by (res) in table 3 are computed. Table 3: Statistics based on STD of different types of anomaly in the test areas (unit: mGal) The differences between the measured and interpolated anomaly at the reference points are computed in each area and for each type of anomaly. A statistical comparison is made. The analysis is made on the standard deviations (STD) of the differences between anomalies. Where the STD are small, the interpolation technique is said to be most accurate. The GGM used in residual anomaly computations is EGM- GGM, the most representative of gravity data in the central African subregion (Kamguia et al., 2007). The residual free-air anomalies FAA(res) are smoother than all the other types of anomaly, except in area 2Z where the total gravity Bouguer anomalies (BA and CBA) are the smoothest (STD = 13 mGal). The difference between the STD of the residual anomalies is nearly 1 mGal. The STD of residual Bouguer and isostatic anomalies are greater than the STD of total gravity Bouguer and isostatic anomalies and than the STD of residual free-air anomalies in some areas. In conclusion, CBA and BA are the smoothest anomalies in 2Z while FAA(res) is the smoothest in 1Z , 3Z and 4Z . FAA(res) is therefore smoothest type of anomaly in Cameroon. It is the only residual anomaly used in the analysis made below. 3.2. Comparison of prediction results The FAA(res) were treated with the total gravity free-air, Bouguer (simple and complete) and isostatic anomalies with the above interpolation methods (and their variants). Table 4 shows the results of statically analysis based on STD of the differences between predicted and known data at the reference points. Table 4: Summary of the STD of interpolation results in test areas (unit = mGal) From Table 4, one can conclude that the accuracy of predicted anomalies depends both Maximum Minimum Mean Standard deviation Anomaly Z1 Z2 Z3 Z4 Z1 Z2 Z3 Z4 Z1 Z2 Z3 Z4 Z1 Z2 Z3 Z4 FAA 71 89 155 67 -30 -60 -40 -30 2,4 26 25 9,9 13 22 32 20 BA 3,7 -46 105 -12 -73 -126 -119 -105 -41 -90 -37 -62 13 13 42 22 CBA 4,3 -43 106 -9 -72 -121 -95 -102 -39 -83 -30 -59 12 13 41 20 IA 68 85 153 65 -30 -63 -44 -32 0,9 23 22 7,7 13 22 39 20 FAA (res) 56 57 87 27 -34 -84 -84 -37 -1,6 5,8 -6 -1,2 11 19 29 11 BA (res) 3 -55 34 -40 -94 -163 -187 -109 -44 -110 -68 -73 15 20 58 13 CBA (res) 4 -52 34 -37 -86 -153 -172 -106 -43 -102 -60 -70 14 18 54 12 IA (res) 53 53 85 26 -36 -87 -87 -39 -3,1 2,4 -8,4 -3,5 12 20 30 12 Kamguia et al. ESRJ Vol. 11, No. 2. December 2007 FAA BA CBA IA FAA (res) Methods Z1 Z2 Z3 Z4 Z1 Z2 Z3 Z4 Z1 Z2 Z3 Z4 Z1 Z2 Z3 Z4 Z1 Z2 Z3 Z4 Minimum Curvature 4 6 7 2 3 4 3 3 4 4 3 3 4 6 7 2 3 6 3 1 Polynomial quadratic 17 27 29 16 12 12 13 19 12 13 16 19 17 26 29 16 10 26 11 13 Polynomial cubic 16 27 27 15 11 12 13 18 11 13 15 19 17 27 27 15 10 26 11 12 Krigging 3 3 6 2 3 3 2 4 3 3 2 4 3 3 6 2 2 3 2 1 Inverse distance (1/r2) 5 5 7 2 4 4 4 4 4 4 4 4 5 5 7 3 3 5 2 1 Inverse distance (1/r3) 4 4 7 2 3 4 3 4 3 3 3 4 4 4 7 2 3 4 2 1 113 Longitude La tit ud e Z1 NIGER CAMEROON NIGERIA CHAD C. A. R. GABON P. R. CONGO E. GUINEA ATLANTIC OCEAN Z2 Z4 Z3 Study area Lac Chad figure 2. Location of the four test areas 13 14 15 16 8 9 10 11 figure 3. Data and control point distributions in + = data point ; ● = reference point ACCuRATE GRAVITy ANoMALy INTERPoLATIoN: A CASE-STuDy IN CAMERooN, CENTRAL AFRICA 114 Kamguia et al. ESRJ Vol. 11, No. 2. December 2007 13 14 15 5 6 7 8 figure 4. Data and control point distributions in + = data point ; ● = reference point 9 10 11 4 5 6 7 figure 5. Data and control point distributions in + = data point ; ● = reference point 115 ACCuRATE GRAVITy ANoMALy INTERPoLATIoN: A CASE-STuDy IN CAMERooN, CENTRAL AFRICA on the interpolation method, the variability of the anomaly used and the test area considered. However, the choice of the type of anomaly to interpolate is more pertinent than the method and the area. The smoother the anomaly interpolated, the less are the interpolation errors, represented here by the STD. For the same type of anomaly interpolated, the krigging and the inverse distance to power 3 methods gave the most accurate results, with approximately the same precision. The best precisions are 2 mGal, 3 mGal, 2 mGal and 1 mGal in test areas ),,,( 4321 ZZZZ respectively. Hence, when the best gridding method is used in conjunction with the smoothest anomalies, a precision of about 2 mGal is obtained after interpolation. 4- Discussions and conclusion Accurate interpolated gravity anomalies are needed in geodetical and geophysical analyses. These anomalies are obtained if the interpolation method and the type of anomaly used are carefully selected. Theoretically, Bouguer gravity anomalies are smoother and therefore suited for prediction processes than free-air anomalies. However, this is only true in an area with simple geology where a constant crustal density of 67.2 is nearest to the real value. For specific geological, geophysical and geodetical studies in such areas, if Bouguer anomalies are interpolated, the information deduced from the gravity interpretations (such as subsurface modelling, basin exploration, geoid computation …) will be biased. Residual free-air anomalies can be smoother than total gravity free-air, Bouguer and isostatic anomalies in many areas. They are easily computed using the spherical harmonic coefficients of a global geopotential model. After reconstruction of the total gravity free-air anomalies from interpolated residual anomalies, any type of gravity anomaly may be deduced. These residual free-air anomalies are associated 11 12 13 14 2 3 4 5 figure 6. Data and control point distributions in + = data point ; ● = reference point 116 Kamguia et al. ESRJ Vol. 11, No. 2. December 2007 to smaller errors in gravity anomaly prediction in Cameroon. The krigging method and the inverse distance to power 3 performed approximately equally as best prediction methods. The geology of Cameroon seems more complicated and only a digital density model may introduce adjusted crustal density values into Bouguer gravity reductions. The number of interpolation methods tested is very limited, just for clarity reasons. If the data were very dense in Cameroon, many grid intervals would have been tested. AcKnowledgments The authors would like to thank the following individuals and organizations: Mr. Sylvain Bonvalot of IRD (Institut de Recherche pour le Developpement), for giving the IRD gravity data during the computation of the first geoid model of Cameroon, and Mr. Bernard Langellier of BGI (Bureau de Gravimetrie International) for making these data accessible. references Adighije, C. I., 1981. A gravity interpretation of the Benue Trough, Nigeria. Tectonophysics, 79, pp. 109-128. Blakely, J. R., 1996. Potential theory in gravity and magnetic applications. Cambridge university Press. 441 pages. Collignon, F., 1968. Gravimetrie de reconnaissance, Cameroun, Rapport oRSToM, 35 p. Crain, I. K. and Bhattacharyya, B. K., 1967. Treatment of non-equispaced two dimensional data with a digital computer. Geoexploration, Vol 5, pp. 173-194. El Abbas, T., Jallouli, C., Albouy, y. and Diament, M., 1990. A comparison of surface fitting algorithms for geophysical data. Terra-Nova, Vol. 2, N° 5, pp. 467-475. Fairhead, J. D. and Okereke, C. S., 1987. A regional gravity study of the West African rift system in Nigeria and Cameroon and its tectonic interpretation. Tectonophysics, 143, pp. 141-159. Heiskanen, W. A. and Moritz, H., 1967. Physical Geodesy. Freeman, Sanfrancisco, 364 pages. Kamguia, J., Manguelle-Dicoum, E., Tadjou, J. M., Tabod, C. T. and Kande, H. L., 2007. The first geoid model of Cameroon (Central Africa). Accepted for publication in the Nordic Journal of Surveying and Real Estate Research (NJSR). Krige, D. G., 1978. Geostatistics for ore evaluation. South African Institute of Mining and Metallurgy, Johannesburg, South Africa. 50 pp. Poudjom-Djomani, Y. H., Diament, M. and Albouy, Y., 1992. Mechanical behaviour of the lithosphere beneath the Adamawa uplift (Cameroon, West Africa) based on gravity data. Journal of African Earth Sciences, 15 (1), pp. 81-90. Smith, W. H. F. and Wessel, P., 1990. Gridding with continuous curvature splines in tension. Geophysics, Vol. 55, N°3, pp. 293-305. Wessel, P. and Smith W. H. F., 1995. New version of GMT released: Transactions of the American Geophysical union 72 (441), pp. 445-446.