| 194 | | 68/2| RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S V G 202 4 GEODETSKI VESTNIK | letn. / Vol. 68 | št. / No. 2| ABSTRACT IZVLEČEK KLJUČNE BESEDE KEY WORDS Gravimetry, Kriging, Variogram, Anisotropy, Cross-Validation.gravimetrija, kriging, variogram, anizotropija, navzkrižna validacija UDK: 528.563:550.831(65) Klasifikacija prispevka po COBISS.SI: 1.02 Prispelo: 12. 2. 2024 Sprejeto: 8. 5. 2024 DOI: 10.15292/geodetski-vestnik.2024.02.194-210 REVIEW ARTICLES Received: 12. 2. 2024 Accepted: 8. 5. 2024 Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA The purpose of this work is to analyze and validate the gravimetric information in the study area in order to highlight the valid data for mapping the gravity anomaly from a grid generated by interpolation. In this context, we have performed a first step of data validation using the Kriging interpolation, this method was applied to a set of free-air anomalies reduced by the EGM08 geopotential model, from the test zone located in the north of Algeria, but the data were considered as values of an isotropic variable whose the variographic analysis depends only on the distances between points. In the second step the anisotropic character was taken into account, where the main tasks were to calculate the directional variograms for validating data and to perform interpolation by directional models. Finally, we have performed a computer application for interpolation and validation of gravimetric data by Kriging method, based on the algorithm for validation and grid generation of free-air anomalies, proposed in this article, where the results showed that the free-air anomaly is an isotropic variable which only depends on the distances between positions, and that the values resulting from interpolation are validated by the BGI grid, where the two maps generated one by interpolation and the other by the BGI grid are very similar to confirm the validity of the Kriging algorithm developed in this context. Namen tega dela je analiza in vrednotenje gravimetričnih podatkov na testnem območju severne Alžirije. S študijo želimo potrditi veljavne podatke za kartiranje anomalij težnosti iz grida, ustvarjene z interpolacijo. Tako smo izvedli prvi korak potrjevanja podatkov z uporabo interpolacije kriging; ta metoda je bila uporabljena za niz anomalij težnosti prostega zraka, reduciranih z uporabo geopotencialnega modela EGM08. Anomalije smo obravnavali kot izotropne spremenljivke, katerih variografska analiza je odvisna le od razdalj med točkami. V drugem koraku smo upoštevali anizotropni značaj, pri čemer sta bili glavni nalogi izračun smernih variogramov za vrednotenje podatkov in izvedba interpolacije s smernimi modeli. Nazadnje smo izvedli računalniško aplikacijo za interpolacijo in vrednotenje gravimetričnih podatkov z metodo kriging. Metoda temelji na algoritmu za vrednotenje in izdelavo grida anomalij težnosti prostega zraka, ki smo ga predlagali v naši raziskavi. Rezultati kažejo, da je anomalija prostega zraka izotropna spremenljivka, ki je odvisna le od razdalj med položaji, in da se vrednosti, pridobljene z interpolacijo, skladajo z vrednostmi anomalij, pridobljenimi iz grida ustanove BGI. Obe gravimetrični karti, naša, ustvarjena z interpolacijo, in karta ustanove BGI, sta si zelo podobni, kar potrjuje veljavnost algoritma kriging, ki smo ga predstavili v članku. Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | SI | E N GV_2024_2_Strokovni-del.indd 194 20. 06. 2024 14:32:17 | 195 | GEODETSKI VESTNIK | 68/2| RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S SI | E N Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | 1 INTRODUCTION Knowing the performances and properties of the Kriging interpolation method (Goovaerts, 1999, Chilès and Delfiner, 2012, Cressie, 2015) this work involves automating the validation process of the gravimet- ric data, and in the end to generate a regular grid from gravimetric data by this geostatistical method. The result of these calculations seems very interesting to the calculation of the local geoid considered as a reference of the altitudes. Often we find in different purviews that Kriging is applied considering that the variable studied and the variogram models are isotropic, which is not necessarily right in reality especially for our case where the variable studied is the gravimetric anomaly whose variogram models are calculated as a function of the distances between points. But the problem is that the variogram models are always sensitive to the anisotropy effect, caused by the variation of the normal gravity field relative to latitude, hence the needs for a new study of the problem by taking into account the anisotropy effect in order to measure the sensitivity of models and interpolation results to the longitude and latitude. In order to satisfy the needs of the problem introduced previously, the theme proposed in this paper aims to test the sensitivity of variogram models to variations in longitude and latitude, to test the sensitivity of the Kriging interpolation to anisotropy and better understanding the influence of anisotropy on data mapping and validation. The general methodology of this work is carried out following certain steps. First, the preliminary analysis that includes the spatial description of data, this first step allows us to clarify the choice of the stochastic model; this leads us to the second step which is the modeling, it is summarized by the variogram modeling which includes the calculation of the isotropic experimental variogram and the directional variograms and their fitting to the analytical models. Finally, interpolation by Kriging can be carried out and this is the subject of the third part, where interpolation is used as a tool for performing data validation and mapping the gravimetric anomaly. 2 VARIABLE DESCRIPTION This description aims at presenting the studied variable before proceeding to the study of the data distri- bution in the geographical space. In this paper it concerns the free-air gravity anomaly known as free-air anomaly. The difference between the measured value of the gravity acceleration g and the theoretical value γ at any point P(x, y, z) of the earth’s surface is caused by the density variation of masses in the earth’s interior; this difference is called “gravity anomaly”. In order to get this difference ∆g, some corrections and reductions must be performed (Heiskanen and Moritz, 1967, Moritz, 1989). Free-Air Anomaly (∆gfa) is defined as the difference between observed gravity on the physical surface (P) and gravity on the telluroid (Q), which is a surface where the normal gravity potential is equal to the actual gravity potential on the physical surface. In the context of this work, the geopotential model EGM08 (Pavlis et al., 2012), developed to the degree and order 2160, has been used to reduce the measured gravity anomalies. The development of the geopotential model (Rapp, 1998) is calculated in the spherical harmonic approximation by the next formula (Heiskanen and Moritz, 1967, Moritz, 1989): [ ]1 01 cos( ) cos( ) (cos ) , n n nm nm nmn m GM a V C m S m P r r λ λ θ+∞ = =   = + +       ∑ ∑ (1) GV_2024_2_Strokovni-del.indd 195 20. 06. 2024 14:32:17 | 196 | | 68/2| GEODETSKI VESTNIK RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S SI | E N Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | Where; (λ, θ, r) are the spherical coordinates of the measurement point, a is the semi-major axis of the reference ellipsoid, GM is the earth’s constant, Cnm and Snm are the normal coefficients of the geopotential model and Pnm are the first specie of Legendre’s functions, the infinite degree of the development can be replaced by Nmax which designates the maximum degree of the geopotential model. By adopting a spherical harmonic model, the gravity value g is given by: g = |grad V |, (2) By reference to the observed gravity value, we obtain: g = g0 + FC, (3) Where: g0 is the observed gravity value and FC is the free-air correction. Afterward, by comparing with the normal gravity γ0 on the reference ellipsoid, the free-air anomaly will be expressed by: ∆gfa = g0 − γ0 + FC, (4) On the basis of the following formula, known as the Somigliana formula (Heiskanen and Moritz, 1967, Moritz, 1989): ( ) 2 2 0 1 2 2 2 2 2 cos sin , cos sin e pa b a b γ ϕ γ ϕ γ ϕ ϕ + = + (5) Where: a and b are the semi-major axis and the semi-minor axis of the reference ellipsoid, γe and γp are the equatorial and polar normal gravity on the ellipsoid, ϕ is the geodetic (ellipsoidal) latitude. The Somigliana formula shows that the normal gravity on the ellipsoid depends on the latitude ϕ, which implies that the gravity anomaly is an anisotropic variable, which varies from one direction to the other, where the two principal directions are the polar and the equatorial directions. 3 INTERPOLATION PROCESS The first part of this application relates to the calculation of the empirical variogram and the second part concerns the fitting of the analytical variogram models (Cressie, 1985, Goovaerts, 1999, Gringarten and Deutsch, 2001) using the empirical variogram function. The third part, which is the most cumbersome in terms of software programming, is designed to the data validation and grid interpolation. 3.1 Kriging Model The choice of the model therefore imposes the type of Kriging to be used in order to perform the in- terpolation (Goovaerts, 1999, Chilès and Delfiner, 2012, Cressie, 2015). Thus, the structure of spatial correlation will be illuminated too and will be characterized in our case by the variogram. A variographic analysis will be carried out in order to better understand the selected variogram models. The basic Kriging model of any spatial variable z has the following form: z(s) = µ(s) + δ(s), (6) Where s is the position, µ is the trend function of the spatial variable that implies which Kriging inter- polation type should be applied, and δ is the stationary part of the spatial variable z. In this work we GV_2024_2_Strokovni-del.indd 196 20. 06. 2024 14:32:17 | 197 | GEODETSKI VESTNIK | 68/2| RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S SI | E N Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | have considered the trend µ constant, which leads us to perform an interpolation by ordinary Kriging. This type of Kriging is often advised by several authors, in particular Deutsch and Journel (1998), who claim that the essential algorithm remains the ordinary Kriging algorithm. 3.1.1 Experimental variogram Empirical variogram calculation is necessary for variographic analysis, as it is calculated directly from the data. The experimental variogram represents the link between data from field measurements and the model adopted to perform spatial interpolation. This calculation step is summarized by the decomposition of all the data into classes according to the distances between the positions of the data and the calculation of the variogram values for each class using the following formula: , 1 1 ( ) ( ) ( ) , 2 k i j N k i jk k h z s z s N γ =  = − ∑ (7) This formula expresses the dissimilarity between the data by a function of the distances between the positions si and sj of each points’ pair. For this purpose, an increment is used and the variogram values for the gravity anomaly are calculated as a function of the distances, to which respectively correspond the variances: , 1 1 ( ) ( ) ( ) , 2 k i j N k i jk k h g s g s N γ =  = ∆ − ∆ ∑ (8) Where ∆gi is the free-air anomaly in the position i, ki,j is the si, sj points’ pair’s number and Nk is the number of points’ pairs distant by a distance hki,j that satisfies: hk − ∆h < hki,j ≤ hk, (9) With; hk = k∆h for k = 1… K, where: K is the number of classes and ∆h is an increment less than the correlation distance. 3.1.2 Fitting variogram models Variographic analysis is the first step in modeling, it consists of finding the analytical variogram function that best suits the calculated experimental variogram. There are several functions in the literature which are used to be taken as a variogram function, we have chosen three functions which are widely used, and which are defined by their parameters; the range r, the sill c and the nugget c0. The calculation of these parameters is carried out by adjusting these parameters to the calculated values of the experimental variogram using an iterative method such as weighted least squares, knowing that the algorithm used in this work is the Gauss-Newton algorithm (Björck, 1996, Madsen et al., 2004, Gratton et al., 2007). However, it is good to always check the accuracy of the estimated parameters. Among the analytical models of variogram, existing in the literature, the chosen models (Marcotte, 2018) are: – Exponential Model 3 0 0( ) ( ) 1 , h rh c c c eγ −  = + − −    (10) GV_2024_2_Strokovni-del.indd 197 20. 06. 2024 14:32:17 | 198 | | 68/2| GEODETSKI VESTNIK RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S SI | E N Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | – Gauss model 2 3 0 0( ) ( ) 1 , h rh c c c eγ  −        = + − −     (11) – Spherical Model 3 0 0 3 1 ( ) if ( ) ,2 2 if h h c c c h r h r r c h r γ    + − − <    =      ≥ (12) Where: r is the range, c is the sill, c0 is the nugget and h is the distance between the positions si and sj of each points’ pair. 3.2 Kriging System For any set of spatial data, the ordinary Kriging interpolation of a point whose coordinates are known can be performed by solving the following system (Kitandis, 1997, Chilès and Delfiner, 2012): 10111 1 1 0 1 . 1 1 ... 1 0 1 γλγ γ γ γ λ γ µ             ⋅ =                      n n nn n n (13) The estimated value is calculated by the following linear equation: 1 * ,n i iiz zλ== ∑ (14) Where γij is the variogram function value of the distance between the two positions si and sj, z * is the estimated value of the variable z, λi is the Kriging weights of the estimated value, zi is the measured value at position si, and n is the number of measurement points. A minimum estimation variance is obtained by minimizing the following expression: 2 var( *),e z zσ = − (15) In the case of a linear estimate, the preceding formula becomes: 2 1 1 1 var( ) cov( , ) 2 cov( , ).n n ne i j i j i i ji j iz z z z zσ λ λ λ= = == + −∑ ∑ ∑ (16) Where: var(z) is the variance of the variable z, and cov(zi, zj) is the covariance value of a points’ pair (zi, zj). To have an unbiased estimate, it is necessary to add the following constraint: 1 1,n ii λ= =∑ (17) This constraint results from the fact that the local mean of the observations is considered constant throughout the study area; in this case the problem is solved by minimizing a quadratic function with an equality constraint. This problem is solved by the Lagrange optimization method to obtaining the Kriging system. GV_2024_2_Strokovni-del.indd 198 20. 06. 2024 14:32:18 | 199 | GEODETSKI VESTNIK | 68/2| RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S SI | E N Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | The estimation variance of this system, called “kriging variance” can be calculated by: 2 1 var( ) cov( , ) ,nk i iiz z zσ λ µ== − −∑ (18) Or through the variogram, as follows: 2 1 ( , ) ( , ) ,nk i ii z z z zσ λ γ γ µ== − −∑ (19) Where μ is the Lagrange factor from the Lagrange optimization method used to minimize the previous variance, and γ (zi, zj) is the variogram value of a points’ pair (zi, zj). 3.3 Interpolation Neighborhood There are many ways to define the geometry and size of a neighborhood, the simplest is to use a circular neighborhood centered at the point to be estimated z0 if the spatial variable is isotropic, and to use an elliptic neighborhood in the case of anisotropic variable. We consider hereafter that the estimation is carried out with a sliding neighborhood comprising k data taking the range r as limit value of the circular neighborhood’s radius for the isotropic interpolation and the ranges rmax and rmin as limit values of the semi-Major axis and semi-minor axis of an elliptic neighborhood in the case of anisotropic interpolation. 3.4 Cross Validation This step involves re-estimating the data for smoothing, the result of this operation is a set of estimated points, but which have the advantage of being linked by a stochastic model. Cross- validation (Dubrule, 1983, Myers, 1993, Zhang et al., 2010) consists to removing data one by one and calculating their estimated counterparts in the same positions by Kriging interpolation using the remaining data, cross validation errors are calculated by subtracting the data values from the estimated values. From these errors, one or more indices are calculated; the index of the root mean square error (RMSE) is often used. This index gives us a better idea on the estimation error that we want to minimize. In this phase cross-validation was used to perform data filtering for the exclusion of data subject to gross errors, and points grouped in narrow areas that can influence the generation of a gravity anomaly grid by Kriging interpolation. This type of filtering guarantees an optimal data quality via a process of automatic detection and rejection of inaccurate data, in order to ensure the reliability of predictions used in mapping the gravity anomaly. Thus, the rejection criterion carried out for this purpose is described by the following formula: 2 2observed predicted observed predicted ,g g k σ σ∆ − ∆ < + (20) In which the absolute difference between the data and its predicted value must not exceed a threshold calculated by the sum of the a-priori input noise and posteriori error estimate (Amin et al., 2005). This technique has been applied in the context of Geoid calculation in several countries, notably by Tscherning (1982) using the collocation method implemented in the GRAVSOFT software (Tschern- ing, 1994). GV_2024_2_Strokovni-del.indd 199 20. 06. 2024 14:32:18 | 200 | | 68/2| GEODETSKI VESTNIK RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S SI | E N Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | 4 DATA VALIDATION AND MAPPING ALGORITHM A software developed for data validation and grid generation by Kriging interpolation was designed as part of this work, based on the algorithm presented in Fig 1. By providing some essential parameters, it is possible to validate a dataset by this program or to calculate values of gravimetric anomaly in the nodes of a grid by Kriging interpolation and finally to mapping the gravimetric anomaly in the study area. The parameters of this program are: – Variogram parameters (range, sill and nugget), which are calculated by a subprogram developed for calculation and fitting the variogram. – The target file (set of points to be validated) – The source file (dataset used in the interpolation) – The reference ellipsoid parameters of the geopotential model (major-axis, minor-axis and flattening) Figure 1: Kriging interpolation flowchart GV_2024_2_Strokovni-del.indd 200 20. 06. 2024 14:32:18 | 201 | GEODETSKI VESTNIK | 68/2| RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S SI | E N Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | The advantage of Kriging method is that it guarantees a linear and unbiased interpolation tool with the property of minimal estimation variance. This variance is a very important for measuring the quality of the interpolation in the study area, since it is defined as the variance of the differences between the real values and the estimated values. In order to better understand the Kriging variance in our area, maps are made using grids generated by Kriging interpolation (Mueller et al., 2004); the calculation requires providing the mesh of grid and the other parameters of the validation program. The algorithm of this program is based on three steps; at the beginning we divide the study area into several cells, these cells are obtained by dividing the interval of the longitude into sub-intervals, and then dividing the interval of latitude in the same manner. The resulting nodes of this process will subsequently be interpolation points from the available dataset. Finally, it is essential in general to mapping the gravity anomaly in order to see better its attitude in the study area, the maps are made using grids calculated by the interpolation program. 5 RESULTS AND DISCUSSION 5.1 Study Area The data set used in this work is an extract of the international gravimetric network, precisely of the Alge- rian network (Fig 2). The study area is located in the northern region of Algeria, between the boundaries of longitude 0° to 8° and latitude 34° to 37°. A set of 1956 point with their free-air anomalies from a file provided by the International Gravimetric Bureau (BGI). Nevertheless, there are serious gaps in the gravimetric coverage of the country which needs to be filled for the different needs of geodesy. Figure 2: Dataset in the study area GV_2024_2_Strokovni-del.indd 201 20. 06. 2024 14:32:18 | 202 | | 68/2| GEODETSKI VESTNIK RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S SI | E N Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | 5.2 Isotropic Interpolation Results 5.2.1 Variogram models The comparison between the RMSE values and the analysis of the graphical representations of each vari- ogram model (Fig 3) show that the Gaussian model adjusts optimally the experimental variogram of the study area. However, it should be noted that the results obtained on the study area remain dependent on the geographical distribution and quality of the data used in this work. The results are represented in the following figure: Figure 3: Experimental variogram with the adjusted function models 5.2.2 Isotropic validation and mapping From the results of validation by an isotropic model (Table 1); the mean, absolute mean and standard deviation calculated from the data set are almost identical to their counterparts for the estimated values. The absolute mean of the differences close to 0.6 also confirms that the estimates of the total dataset are good. Therefore, the set is well validated but also specifying that the data whose estimates exceed the threshold of 20 mGal will be rejected, we found that 92.28% of the differences are below the threshold, which depends on the accuracy and density of the data and is used in the GRAVSOFT software. Table 1: Statistics of data and estimated values by the isotropic model Variables Min (Δgfa) [mGal] Max (Δgfa) [mGal] Mean (Δgfa) [mGal] S.d. (Δgfa) [mGal] Data -82.59 165.37 28.75 31.97 Estimates -65.81 137.15 29.32 31.95 Differences -71.43 83.63 0.57 2.28 In Fig 4 we perceive that Kriging variances decrease when there is a large presence of data and increases if there is a lack of data as in the case of two zones; the North where a part of the Mediterranean Sea appears and the South east, where there is a total absence of data, the study of the Kriging variances allowed us to even proclaim that the non-homogeneity of the gravimetric network in the study area much influ- GV_2024_2_Strokovni-del.indd 202 20. 06. 2024 14:32:18 | 203 | GEODETSKI VESTNIK | 68/2| RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S SI | E N Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | ences the interpolation. This implies that the results of the free-air gravimetric anomaly estimation and subsequently the calculation of the Geoid in this region remain dependent on the gravimetric network. Figure 4: Map of Kriging isotropic variance using a mesh grid 2.5’×2.5’ The results of the interpolation in the study area are represented by a map in Fig 5, these results are very similar and comparable to the values of the grid downloaded from the BGI website, which in turn are represented by a map in Fig 6, with a mesh grid 2.5’×2.5’. Figure 5: Map of Free-air anomaly by Kriging isotropic model using a mesh grid 2.5’×2.5’ In order to comment on the quality of free-air anomaly interpolation from all the gravimetric data in the study area, we have chosen to create a grid by interpolation from our available data, with the same mesh as that of the EGM08 model provided by the BGI, and to compare the two grids to note the advantages and disadvantages of free-air anomaly interpolation by Kriging. The statistics of the grid produced by Kriging and that of the EGM08 model represented in table 2, show us the presence of certain aberrant values due to the absence of data in certain regions, despite the fact that we have removed the majority of the values obtained in the lacking data regions. The averages of the two grids are almost identical and the average of the differences equal to 0.2 mGal is almost zero, this allows us to see that the difference between the two grids is stable and that the majority of the points have close values between the two grids. GV_2024_2_Strokovni-del.indd 203 20. 06. 2024 14:32:19 | 204 | | 68/2| GEODETSKI VESTNIK RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S SI | E N Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | Table 2: Statistics of the Free-air anomaly grid estimated by kriging and the EGM08 model grid of the BGI Variables Min (Δgfa) [mGal] Max (Δgfa) [mGal] Mean (Δgfa) [mGal] S.d. (Δgfa) [mGal] EGM08 Grid -97.29 158.18 23.75 41.94 Estimated Grid -99.79 169.98 23.77 42.12 Differences -205.88 246.71 0.02 4.35 Figure 6: Map of the EGM08 Free-air anomaly using a 2.5’×2.5’ mesh grid downloaded from the BGI website 5.3 Anisotropic Interpolation Results 5.3.1 Directional variograms In our case, if we return to the description of the variable, we can see that the normal gravity γ0 varies accord- ing to the latitude ϕ, due to the flattening of the reference ellipsoid used in the calculation of the terrestrial gravity field model. This implies that the gravity anomaly varies according to the latitude, in this case it can be considered as an anisotropic variable, knowing that the two principal directions are parallel to the merid- ians and to the equator respectively. In Fig 7 we can see that the two variograms have the same range, but with two different sills, it leads us to think that the gravity anomaly presents a zonal anisotropy in this area. Figure 7: Directional variograms with the Gaussian function models GV_2024_2_Strokovni-del.indd 204 20. 06. 2024 14:32:19 | 205 | GEODETSKI VESTNIK | 68/2| RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S SI | E N Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | 5.3.2 Geometric anisotropy In the following, we will use the anisotropic variogram models in order to get a final result il- lustrating the influence of anisotropy on interpolation and validation of gravimetric data. In order to apply an interpolation by an anisotropic model we made a projection of all the dataset on the tangent plane for each point and that to introduce the analytic expression of the angle θ used to generate the elliptical neighborhood of anisotropic interpolation. In this context, a polar stereographic projection (Snyder, 1987, Bugayevskiy and Snyder, 1995) was used which is characterized by: 0 0 2 tan sin 4 2 , 2 tan cos 4 2 X Rk Y Rk π ϕ λ π ϕ λ   = −        = − −    (21) Where (ϕ, λ) are the geographic coordinates of the projected point, R is the mean radius of the earth and k0 is a scale factor. Then, rθ can be calculated as a function of rmax and rmin by the following formula (Marcotte, 2018):  ( ) max min 1 2 2 2 2 2 min max * , cos * sin r r r r r θ θ θ = (22) The coordinates taken in this step are plane coordinates (X, Y ) from the stereographic projection. Inter- polation and validation were performed using a combined range rθ for the total dataset. 5.3.3 Anisotropic validation and mapping The results presented in table 3 show that the differences between the data values and their estimates by the anisotropic model are low, without taking into account some extreme values that influence the calculation of the minimum and the maximum. The absolute mean of the differences calculated by anisotropic model is close to 0.9 mGal, confirms the results of validation by the isotropic model, but it is superior than the absolute mean of the isotropic model, whose is almost equal to 0.6 mGal. A rate equal to 89.17% of data where the differences are below the threshold of 20 mGal shows that the majority of data are accepted, but remains less than the rate obtained by the isotropic validation which is equal to 92.28%. Table 3: Statistics of data and estimated values by the geometric anisotropy model Variables Min (Δgfa) [mGal] Max (Δgfa) [mGal] Mean (Δgfa) [mGal] S.d. (Δgfa) [mGal] Data -82.59 165.37 28.75 31.97 Estimates -79.14 153.60 21.00 32.04 Differences -113.01 149.55 0.89 2.46 To better express the quality of the estimated values by anisotropic Kriging interpolation, a map of anisotropic Kriging variance is made using a mesh grid 2.5’×2.5’. GV_2024_2_Strokovni-del.indd 205 20. 06. 2024 14:32:19 | 206 | | 68/2| GEODETSKI VESTNIK RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S SI | E N Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | Figure 8: Map of Kriging anisotropic variance using a mesh grid 2.5’×2.5’ In Fig 8 we notice a similarity between the map of Kriging variances resulting from the anisotropic model and that resulting from the isotropic model, of which we confirm the poor quality of the estimation in the northwest zone which coincides with the Mediterranean Sea and in the south-east zone where we consider a very low and insufficient presence of data. Nevertheless, we observe areas in the middle of the map where the Kriging variance is high, this is explained by the presence of outliers or aggregations of points close to each other, which influence the calculation. The following map (Fig 9) presents the gravity anomalies estimated by Kriging using the anisotropic model, which is very similar to the map made using the isotropic model and to the map of the BGI grid. Figure 9: Map of Free-air gravity anomaly by Kriging anisotropic model using a mesh grid 2.5’×2.5’ 6 COMPARATIVE STUDY In order to highlight the impact of the interpolation method on the prediction results obtained, and thereby to comment on our investigations, we deemed it useful to compare the values estimated by Krig- ing method with those obtained by other spatial interpolation methods. In this case, we used the Least Squares Collocation (LSC) and the Inverse Distance Weighting (IDW) methods. Least squares collocation (Moritz, 1978) is an optimal estimation method implemented in the GRAVSOFT software, applicable in GV_2024_2_Strokovni-del.indd 206 20. 06. 2024 14:32:19 | 207 | GEODETSKI VESTNIK | 68/2| RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S SI | E N Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | particular to gravimetric data. By this method, we carried out the estimation of a gravimetric signal from measurements taken at points scattered in the study area using a covariance function for expressing the covariance between the different measurements. The optimal model of the covariance function chosen in this study is the Exponential model with a spherical correlation distance ξ = 0.15° (~15km) and a variance C0 = 557.20 mgal 2. Among others, the inverse distance weighting method is part of interpolation methods class called «moving average methods», which predict the value of a regionalized variable at a point s0 by a weighted average of observations located in a neighborhood V(s0) defined a priori. The weights of this average are a function of the distance between the observed point and the estimation point, such that the closest points have more influence in the interpolation, therefore, this method give zero weight to the observations outside the interpolation neighborhood. By the inverse distance weighting method of a power d, we calculate the estimated value of a regionalized variable z at a point s0 by the weighted average: 0 0 0 0 ( ) 0( ) 1/ * ( ) ( ), 1/i i d i ids V s is V s s s z s z s s s∈ ∈ − = − ∑ ∑ (23) With: d is a strictly positive integer, si is the measurement point and z(si) is the observed value of the point si. Cross-validation on the data from the study area was carried out using respectively the inverse distance method weighting and the least squares collocation method in which the Exponential model with its adjusted parameters was adopted as a local covariance model to express the correlation between the observation values. The statistics of the prediction results obtained as well as their differences using the three interpolation methods are presented in table 4. Subsequently, the following histograms (Fig 10) show clearly that the LSC–Kriging difference is mainly within the amplitude interval ±5 mGal, this is due to the theoretical design of the two methods based on the principle of spatial correlation between the data, and that the IDW–Kriging difference is generally included in an amplitude interval ±20 mGal but the majority of differences are in the amplitude interval ±10 mGal, and this is due to the deterministic principle of the IDW method which is different to that of Kriging. Table 4: Statistics of estimation results by three interpolation methods Variables Interpolation method Min (Δgfa) [mGal] Max (Δgfa) [mGal] Average (Δgfa) [mGal] S.d. (Δgfa) [mGal] Data Free-air Anomalies -82.59 165.37 28.75 31.97 Estimates Kriging -65.81 137.15 29.32 31.95 LSC -54.81 140.37 29.18 31.74 IDW - 42.26 128.49 24.36 27.44 Differences LSC – Kriging -26.91 12.41 -0.01 0.59 IDW – Kriging - 36.66 29.27 0.19 1.34 The percentages of deviations below the threshold of 20 mGal are in the order of 92.1% and 87.5%, respectively, for the least squares collocation method and the inverse distance weighting method respec- tively. Therefore, the comparison of the estimation results shows that more than 90% of the doubtful points could be detected by the two interpolation methods, namely Kriging and LSC. This is due to the fact that these two methods are by definition linear estimates of the observations vector, and they are based on the principle of spatial correlation through the use of a covariance or variogram functions. GV_2024_2_Strokovni-del.indd 207 20. 06. 2024 14:32:20 | 208 | | 68/2| GEODETSKI VESTNIK RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S SI | E N Figure 10: Histograms of differences between methods 7 CONCLUSION The method used in this article is of major utility for spatial variables such as gravimetric measurements. All these measurements represent a data source used in general for the modeling of the terrestrial gravity field and its various components, particularly in the calculation of the geoid heights. First, interpolation and validation were performed assuming gravimetric anomaly as an isotropic variable, where variographic analysis depends only on the distances between points in all directions, but the Somigliana formula expresses the dependence of the normal gravity versus latitude, hence we consider that the interpolation of the gravity anomaly also depends on the latitude, from which the need arises to study the anisotropic nature of this variable in order to better understand its impact on our calculations. The comparison is done between the results of validation by the isotropic model and validation by directional models, the two principal directions of which are the direction of the meridians and of the equator. In this context, it results that the maximum range of the variogram is that calculated in the direction of the meridians. The results of validation using the anisotropic model are almost similar to the results obtained using the isotropic model; these results allow us to say that the estimation of the gravity anomaly by Kriging method as well as the validation does not depend on the direction. The algorithm developed in this work provided results of the gravity anomaly estimation as a regular grid which made it possible to produce maps with good quality, either by isotropic or anisotropic model, which allow an objective representation of the studied variable in the study area. These results are validated by the grid produced by the BGI and by the map obtained by this grid which is very similar to those obtained by applying the algorithm proposed in this work. Finally, the results allow us to say that the interpolation and validation of the free air anomalies by Kriging method in this area are not very influenced by the direction when calculat- ing the distances between points, and to even consider that the gravimetric anomaly is isotropic in our calculations but taking into account the sufficient coverage of the study area by the data and the zonal anisotropy, which requires other methodology and the division of the study area into several zones, based on detailed information of the mass density, before performing an interpolation. In order to exhibit the performance of the method used in this work a comparison with other gridding methods by spatial interpolation was applied, precisely the Least Squares Collocation method (LSC) used Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | GV_2024_2_Strokovni-del.indd 208 20. 06. 2024 14:32:20 | 209 | GEODETSKI VESTNIK | 68/2| RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S SI | E N in GRAVSOFT software and the Inverse Distance weighting method (IDW). The results of validation by Kriging method using the Gauss model as variogram function shows that 92.28% of the differences are below the threshold, and that allows us to confirm that the data provided by the BGI are of good quality and can be used in the geoid calculation of Algeria. In addition, the error rate calculated using the least squares collocation method remains close to that derived by the Kriging method and shows that 90% of the doubtful points could be detected by the two interpolation methods. However, unlike these two interpolation methods, the inverse distance weighting method has a slightly high error rate compared to the other methods and the non-existence of theoretical variance which constitutes crucial information on the precision of the estimated values obtained and consequently a means of deciding on the quality of the dataset used in this study. In this study we presented the Kriging method as a tool for free-air anomalies interpolation in order to create a grid for mapping this variable, but the disadvantage in this approach is the strong correlation of free-air anomaly to the topography, knowing that the objective of the work was to create an algorithm as well as a software to automate the methodology of this work. However, Kriging interpolation provides a smooth mesh which does not correlate with the topography in general, but the effect of the topography on free-air anomaly is often removed by the so-called residual terrain model method (RTM), knowing that we can also create a grid directly from the complete Bouguer anomalies (CBA) which are smoother than free-air anomaly (Čunderlik et al, 2021). The methodology as well as the algorithm of this work is feasible on all gravimetric data or on the Geoid undulation in order to carry out a mapping of these variables, thus the perspective of this work is to involve this methodology in the mapping of different levels of the gravity anomaly or in the interpolation of the residuals in the RTM method for further investigations in future studies. Literature and references: Amin, M.M., El-fatairy, S.M., Hassouna, R.M. (2005). A precise Geoidal map of the southern part of Egypt, by collocation: Toshka Geoid. FIG Working Week 2005 and GSDI-8, 16–21 April 2005, Cairo, Egypt. Björck, A. (1996). Numerical Methods for Least Squares Problems. Society for Industrial and Applied Mathematics, Philadelphia. Bugayevskiy, L. M., Snyder, J. P. (1995). Map Projections: A Reference Manual. Taylor and Francis, London. Chilès, J.P., Delfiner P. (2012). Geostatistics: Modeling Spatial Uncertainty, 2nd Edition. John Wiley & Sons, New York. Cressie, N. (1985). Fitting variogram models by weighted least square. Mathematical Geology, 17, 563–586. DOI: http://dx.doi.org/10.1007/BF01032109 Cressie, N. (2015). Statistics for Spatial Data, Revised Edition. John Wiley & Sons, New York. Čunderlik R., Tenzer, R., Macak, M, zahorec, P., Papčo, J., Ababio, A.N. (2023). A detailed quasigeoid model of the Hong Kong territories computed by applying a finite-element method of solving the oblique derivative boundary-value problem. Journal of geodetic science, 13, 1. 20220153. https://doi.org/10.1515/ jogs-2022-0153 Deutsch, C. V., Journel, A. G. (1998). GSLIB: Geostatistical Software Library and User’s Guide. Oxford University Press, New York. Dubrule, O. (1983). Cross validation of kriging in a unique neighborhood. Mathematical Geology, 15, 687–699. DOI: http://dx.doi.org/10.1007/BF01033232 Goovaerts, P. (1999). Geostatistics in soil science: state-of-the-art and perspectives. Geoderma, 89, 1–45. DOI: http://dx.doi.org/10.1016/S0016-7061(98)00078-0 Gratton, S., Lawless, A. S., Nichols, N. K. (2007). Approximate Gauss–Newton methods for nonlinear least squares problems. SIAM Journal on Optimization, 18 (1), 106–132. DOI: http://dx.doi.org/10.1137/050624935 Gringarten, E., Deutsch, C. V. (2001). Teacher’s aide: variogram interpretation and modeling. Mathematical Geology, 33, 507–534. DOI: http://dx.doi. org/10.1023/A:1011093014141 Heiskanen, W., Moritz, H. (1967). Physical Geodesy. W. H., Freeman and Company, San Francisco and London. International Gravimetric Bureau - BGI. http://bgi.obs-mip.fr/en , accessed 24. 11. 2023. Kitandis, P.K. (1997). Introduction to Geostatistics: Applications to Hydrogeology. Cambridge University Press, Cambridge. Madsen, K., Nielsen, H.B., Tingleff, O. (2004). Methods for non-linear least squares problems, 2nd edition, lecture notes. Informatics and Mathematical Modeling, Technical University of Denmark, Lyngby, 58p. http://www2.imm.dtu.dk/ Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | GV_2024_2_Strokovni-del.indd 209 20. 06. 2024 14:32:20 | 210 | | 68/2| GEODETSKI VESTNIK RE CE NZ IRA NI ČL AN KI | P EE R- RE VIE W ED AR TIC LE S SI | E N pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf , accessed 26. 6. 2023. Marcotte, D. (2018). Chapter 2: the variogram. In Course GML6402: geostatistics. Polytechnic School of Montreal, Montreal, 22p. https://cours.polymtl.ca/geo/ marcotte/gml6402.html , accessed on 26. 6. 2023. (In French) Moritz, H. (1978). Least-squares collocation. Reviews of Geophysics, American Geophysical Union, 16, 421–430. DOI: https://doi.org/10.1029/ RG016i003p00421 Moritz, H. (1989). Advanced Physical Geodesy, 2nd Revised Edition. Herbert Wichmann Verlag, Karlsruhe. Mueller, T.G., Pusuluri, N.B., Mathias, K.K., Cornelius, P.L., Barnhisel, R.I., Shearer, S.A. (2004). Map quality for ordinary kriging and inverse distance weighted interpolation. Soil Science Society of America journal. 68, 2042-2047. DOI: http://dx.doi.org/10.2136/sssaj2004.2042 Myers, D. (1993). Cross-validation and variogram estimation. Theory of Probability and Applications, 37, 345–347. DOI: http://dx.doi.org/10.1137/1137074 Pavlis, N.K., Holmes, S.A., Kenyon, S.C., Factor, J.K. (2012). The development and evaluation of the earth gravitational model 2008 (EGM2008). Journal of Geophysical Research, 117, 1–38. DOI: https://doi.org/10.1029/2011JB008916 Rapp R.H. (1998). Paste and future developments in geopotential modelling. In: Forsberg R, Feissel M. and Dietrich R., (eds), Geodesy on the Move. International Association of Geodesy Symposia, Springer, Heidelberg, Berlin, 119, 58–78. DOI: http://dx.doi.org/10.1007/978-3-642-72245-59 Snyder, J. P. (1987). Map Projection: A Working Manual. US Geological Survey Professional Paper 1395, US Government Printing Office, Washington, DC. Tscherning, C.C. (1982). Geoid determination for the Nordic countries using collocation. Proceedings of the General Meeting of the International Association of Geodesy, 7–15 May 1982, Tokyo, Japan. Tscherning, C. C. (1994). Geoid determination by least squares collocation using GRAVSOFT. Lecture notes for the International School for the Determination and Use of the Geoid, 10–15 October 1994, Milano, Italy. Zhang H. and Wang Y. (2010). Kriging and cross-validation for massive spatial data. Environmetrics, 21, 290–304. DOI: https://doi.org/10.1002/env.1023 Ilies Benikhlef Center of Space Techniques (CTS), Department of Spatial Geodesy, 01 Avenue de la Palestine B.P. 13, 31200 Arzew, Oran, Algeria. Oran 1 University Ahmed Ben Bella - Faculty of Exact and Applied Sciences - Department of mathematics, B.P. 1524, El M'Naouer - 31000 Oran, Algeria e-mail: benikhlef.ilies@gmail.com prof. dr. Abderrahmane Senoussaoui Oran 1 University Ahmed Ben Bella - Faculty of Exact and Applied Sciences - Department of mathematics, B.P. 1524, El M'Naouer - 31000 Oran, Algeria. e-mail: senoussaoui.abderrahmane@univ-oran.dz dr. Sid Ahmed Benahmed Daho, Research Director. Center of Space Techniques (CTS), Department of Spatial Geodesy, 01 Avenue de la Palestine B.P. 13, 31200 Arzew, Oran, Algeria. e-mail: sbenahmed@cts.asal.dz Benikhlef I., Senoussaoui A., Daho S. A. B. (2024). Kriging algorithm for mapping the gravimetric data in the north of Algeria. Geodetski vestnik, 68 (2), 194-210. DOI: https://doi.org/10.15292/geodetski-vestnik.2024.02.194-210 Ilies Benikhlef, Abderrahmane Senoussaoui, Sid Ahmed Benahmed Daho | ALGORITEM KRIGING ZA KARTIRANJE GRAVIMETRIČNIH PODATKOV NA SEVERU ALŽIRIJE | KRIGING ALGORITHM FOR MAPPING THE GRAVIMETRIC DATA IN THE NORTH OF ALGERIA | 194-210 | GV_2024_2_Strokovni-del.indd 210 20. 06. 2024 14:32:20