Study on the impact of calculation radius on IDW gravity modelling

Gravity modelling is a necessary process of applying gravity data to practice, and a vital research area in resource exploring and geodesy. Inverse distance weighted (IDW) interpolation is a widely used algorithm in gravity modelling and many other griding procedures. Based on the gravity anomaly data acquired from the airborne gravimetry over a region in the Northwest of China, IDW method is analysed and implemented with varying calculation radius in gravity modelling. The results show that over the study area, with a proper calculation radius, the IDW can give a more precise prediction than quadric surfaces fitting (QSF) method. The distances between sampled points and unsampled points are the only information taken into consideration in the IDW interpolation, so the analysis of the influence by changing the calculation radius of the unsampled points in this test provides a new research direction of error eliminating when using IDW.


Introduction
Last decades have seen large progresses made in the area of gravimetry, with different means of gravimetry developed rapidly, which makes it possible to acquire a huge amount of gravity data like gravity anomalies and disturbances. As a result, there is an urgent need to comprehensively use those gravity data to build gravity models both regionally and globally. Data used in the recovery of gravity field is usually discrete and randomly distributed and sometimes with different resolutions due to various of ways of measurements [1] . As a consequence, algorithms of interpolation are various, including inverse distance weighted, Kriging interpolation, quadric surfaces fitting and least squares collocation [2] . Different ways of modelling have different features. The advantage of IDW is that it only takes distance information into calculation when predicting the unknown points, which largely simplifies the calculation process and easy to implement. While the main weakness of this algorithm is the absence of relevant physical information over the study region which results in a low accuracy of the prediction [3] . Least squares collocation is a theory of function approximation, which has been developed gradually since the 1960s [4] . Least squares collocation is a complex and precise method in gravity modelling, especially when dealing with two or more kinds of gravity data, but the amount of calculation is usually large because of the computation of inversing large matrixes. The applicable scope of Kriging method is that the regional variables have spatial correlation, that is, if the results of variogram and structural analysis show that the regional variables have spatial correlation, then Kriging method can be used for interpolation, otherwise, it is not feasible. Therefore, considering the complexity and convenience of the calculation, IDW and QSF are selected.
The main objective of this study is to figure out the relation between the calculation radius and the prediction's accuracy, which is of vital importance to the further research on the determination of calculation radius. IDW and QSF are implemented to generate a regional gravity model based on a set of discrete gravity anomaly data over the survey region [5] . Datasets are extracted from the original data collected during the flight tests of airborne gravimetry conducted by National University of Defense Technology (NUDT), using a strapdown airborne gravimeter, which has an internal accuracy within 1.6mGal [6] . To see the impact of different calculation radiuses of IDW, the radius is set to grow step by step from 0.02° to 0.55° (considered in geodetic coordinate system, so the angular unit is used). This range of radius shows the comprehensive outcome according to the varying radius. With a certain properly selected radius in IDW, the error can be largely eliminated, performing better than QSF method.

Data source
In August, 2016, flight tests were conducted over a region in the Northwest of China, which is about 1600km 2 . In the datasets from the task of airborne gravimetry, the flying height data is averagely 1500m, and the range of the latitude covers about 0.35° and the range of the longitude covers about 0.40°.
There are 1075 points of gravity anomaly in the dataset extracted from the original dataset. For convenience, we assign a name for the dataset with 1075 points of gravity data, called DS, then we divided DS into two parts, DS1 and DS2. DS1 is employed as the experimental set, DS2 is used as the test set, which means the gravity model obtained from DS1, can be validated by the DS2. There are similar numbers of points in both DS1 (538 points) and DS2 (537 points). These two datasets are uniformly distributed in the study area.  Table 1 shows both DS1 and DS2 are similar to DS statistically, indicating three datasets have nearly the same statistical features, which means the idea of grouping is feasible.

Inverse distance weighted (IDW)
IDW is a widely used interpolation algorithm in the research of gravity field modelling and the computation of GPS height anomaly. IDW is applied to calculate the unsampled points according to other known points over the survey region. The basic principle is that the nearer a known point is, the bigger influence it has on the unsampled point, while the larger the distance is, the less impact it has on the unknown point. First, distances between the unsampled point and the rest of the points need to be calculated within the twodimensional plane.
In equation (1) i Dis is the distance between the th i point and the unsampled point. When applying IDW in gravity modelling, all distances between the sampled data and the unknown points should be calculated. When there is not a large number of points of gravity data, it would be difficult to figure out the drawbacks of the IDW method. When the amount of data is large, the computation amount of the algorithm grows rapidly, which makes it a low efficiency algorithm. Mostly, the amount of gravity data is large enough that it takes a long time to accomplish the computation work. As a result, selecting a limitation for the calculation radius can help improve the efficiency.

Quadric surfaces fitting (QSF)
The mathematical model of quadric surfaces fitting can be described as the formula below:   In formula (2),   , x y is the position of points in DS1, and i a stands for the 6 coefficients in the model, so there should be no less than 6 known points in the DS1. In fact, the quantity of the data in DS1 is much larger than 6, so when the number of known points is n , which is larger than 6, the error equation is as follows: x y x y x y x y x y x y x y x y x y According to the principle of least square algorithm, we can solve the equation and give the solution of the coefficients vector X .
In formula (4), L is the vector of the magnitude of gravity anomaly in DS1.

Methods validation and results
There are several steps in the whole process of gravity modelling and precision analysis.
1. Implementing the traditional IDW method without setting a radius limitation. Calculating of the gravity anomaly model by using quadric method.
2. Analysing the precision of the results. 3. Conducting the IDW with different limitations on the calculating radius. Comparing the results with the results of quadric method acquired in the second step.    The root mean square error (RMSE) of IDW without limitation of calculation radius is 7.98mGal. As is shown in table 2, standard deviation, mean value, minimum and maximum are 7.98, -0.38, -24.31 and 17.43 respectively.

Implementation of IDW and QSF
The RMSE of QSF is 4.20mGal, which is smaller than the result of IDW without radius limitation. Table 3 shows the standard deviation, mean value, minimum and maximum are 4.20, -0.034, -13.16 and 13.61 respectively.
By comparing the results from the tests of IDW and quadric method, it can be figured out that IDW algorithm results in a greater RMSE than the latter does. Besides, QSF method performs better in terms of the mean value and max-min value of modelling error.

IDW with limitation on calculation radius
In the algorithm of IDW, the distance between the sampled points and the unsampled points is the only information taken into calculation. Results show the accuracy is improved by controlling the calculating radius at a certain low level.   Fig. 7. Gravity anomaly over the study region after griding by IDW without limitation on calculation radius. Fig. 8. Gravity anomaly over the study region after griding by radius limited IDW. Figure 5 shows a minimum RMSE of the IDW method is attainable when the calculation radius is selected properly, in this context, the radius should be set at 0.031° and the minimum RMSE is 2.58mGal, which is much lower than that of QSF. When the limitation of the radius grows more than 0.031°, error of this method is increased gradually. Figure 7 and 8 show that the whole picture by IDW with and without limitation on calculation radius. The comparison between these two figures indicates that the improved method gives more details than the normal IDW over the survey region.
Error of the result from radius constrained IDW performs ideal statistical properties that mean value is -0.05mGal, standard deviation is 2.58mGal, min-max value is 20.45mGal, indicating that with a constrained calculating radius, results of IDW can be stable and ideal.

Conclusion
The results from the calculation based on the airborne gravity anomaly data show that the calculating distance has a significant impact on the precision of the unsampled points' prediction. According to this study, when the radius rises to a level higher than 0.031°, the error increases with the distance growing. On the contrary, with a distance less than 0.031°, IDW performs worse or even unable to complete the calculation due to the lack of sampled points. With a shorter radius computed in the IDW, amount of the whole calculation in the algorithm is decreased due to fewer points of data are taken into calculation.
As a result, with a properly chosen calculation radius when applying IDW, the prediction's accuracy can be improved to 2.58mGal, which is better than the accuracy of QSF method with 4.20mGal of RMSE. Selecting a best radius for the algorithm is a nearly impossible task when constructing a model for the gravity data in a real scenario without knowing the data to be estimated, while in this study, these two datasets from the same original dataset perform similar statistical properties, so it is possible to figure out a proper radius by assigning one of the datasets a predicting dataset and another a validating dataset to give a similar result of precision as in a real predictive calculation. The quantitative relation between the radius constraint and the RMSE of the IDW will be the focus in future research.