Data Processing of Gravity Base Network in Plateau Area: The Case of Qinghai Province, China

: The latest gravity survey of the gravity base network in Qinghai Province, China, was conducted with six Scintrex CG gravimeters and this gravity survey was tied to existed gravity reference stations. In this gravity network with long segments and very rugged topography, the calibration of scale factors is a time-consuming progress and its accuracy may be affected by many uncertainties, and the change in drift rates of the relative gravimeters are complex over time in this long-term survey. The reasonable calculation of scale factors and drift rates plays an important role in improving the gravity estimation accuracy. In this paper, based on the least squares, robust least squares, and Bayesian methods, various parameter calculation methods were employed to process this gravity network. The performance and practicality of each method were analyzed in terms of internal and external accuracy. The results indicated that the scale factors calibrated in the baseline ﬁeld had poor applicability due to insufﬁcient gravity difference, in this case, the scale factors estimated by the adjustment models were more accurate, which weakened the correlation between gravity differences and mutual differences. The drift rates estimated by the Bayesian method were relatively smooth over time, while drift rates estimated using symmetric observations were more practical for the gravimeter with highly variable drift. The weight constraints of observations can be optimized by the robust least squares method, the gravity values obtained by it were more consistent with absolute gravity values than those obtained by the least squares method, and the robust least squares method was recommended to process gravity data in plateau areas.


Introduction
A gravity base network forms the basis for gravity surveys and provides accurate gravity field information for mineral resources, geotectonic inversion, geoid refinement, and geodynamic research [1][2][3][4]. China's gravity base networks are established by relative gravity survey campaign based on a number of absolute gravity points [5]. The gravity base network of Qinghai Province, China, was established to fulfill the needs of infrastructure construction and geodesy research. The first phase of gravity measurements in this network was performed in the form of hybrid gravimetry, and five absolute gravity points and three relative gravity points were measured from 2017 to 2018. Later, in order to achieve province-wide coverage of this gravity network, the latest phase of gravity survey was conducted from September 2020 to October 2020 using two Scintrex CG-6 and four Scintrex CG-5 gravimeters.

Absolute Gravity Instrumentation
A077, A078, A079, A080, and A082 were measured by absolute gravimeters FG5X#246 and FG5#214. FG5X#246 belongs to the Institute of Geodesy and Geophysics, Chinese Academy of Sciences (IGG CAS), Wuhan, China. FG5#214 belongs to National Administration of Surveying, Mapping, and Geoinformation of China. In order to ensure the reliability of reference and metrology, ensure the long-term stability of all gravimeters, and trace the gravity measurement to the SI, Coordination Group for Absolute Gravity  Moreover, FG5#214 took part in the 10th International Comparison of Absolute Gravimeters (ICAG-2017), which was carried out on the Changping Campus of NIM, Beijing, China. The results showed that FG5#214 provided accurate measurements and was reliable in terms of metrology [14].
The effects of solid-earth tide, ocean tide loading, and barometric were calculated by the g9.0 software. ETGTAB model developed by Georg Wenzel was adopted to correct the effect of solid-earth tide, and FES2004 model, considered state-of-the-art, was used to correct the effect of ocean tide loading [5]. The barometric admittance factor of −0.3 µGal/hPa was adopted for atmosphere correction [15,16]. The effect of vertical gravity gradient was corrected based on the measured vertical gravity gradient at each gravity points. The correction of polar motion was estimated based on the pole positions published by the International Earth Rotation and Reference Systems Service (IERS) [15].

Hydrological Effects on Gravity
Seasonal effects of global hydrology on the first phase of gravity measurements were removed before adjustment. The effects were divided into loading and gravitational parts, the calculation was divided into several areas according to the spherical distance between the mass and the measurement point, and the contribution of each area was calculated using the formula provided by Farrell [18]. In this study, the GLDAS model was used to calculate the effect of global land water on reference points ( Table 2). The effects of hydrology are as follows: where ψ is the spherical distance between the mass and the measurement point, g E (ψ) is the loading effect per unit mass, g N (ψ) is the gravitational effect per unit mass, g represents the mean surface gravity, M represents the mass of the Earth, h n and k n represent the load Love numbers, and P n represent the Legendre polynomials.

Relative Gravity Survey
The latest phase of the gravity survey lasted nearly 45 days. A total of 404 measurements were obtained using 6 relative gravimeters (2 Scintrex CG-6, and 4 Scintrex CG-5) on 35 segments, which included 10 closed loops (the maximum segment number of loops was 5), 1 annexed gravimetric line (B069-WUDL-YEGE-A082), and 1 branch line (BUQL-TUQL) ( Table 3). The 6 relative gravimeters were checked regularly, and performance tests were carried out on all gravimeters before gravity survey according to technology standard [19]. line covered a gravity difference of approximately 430 mGal ( Figure 1). The calibration line followed the road between B070 and B071 with a distance of 270 km, and the travel time of interval was less than 4 h, enabling the calibration to be completed within one day.
At each gravity point, the gravimeters were set for 3-5 min before data collection to obtain the data measured by gravimeters at steady state. The data collection required 5 series of 60 one-second measurements. Symmetrical survey scheme (A − B − C · · · · · · C − B − A) was applied, and each loop was closed within 48 h. The geodetic coordinates and altitude of the gravity points were surveyed based on the Qinghai CORS and the refined local geoid of Qinghai (root mean square (RMS) of 10 cm).

Relative Gravity Data Preprocess
Tide force causes periodic changes in the gravity value at each points, thus requiring the removal of its influence [20]. Longman formulas were used to correct the effect, which are expressed as follows [21]: and where δ gb (10 −8 ms −2 ) is the correction value for the solid-Earth tide; δ th is the tide factor, which was taken as 1.16; c m is the average distance from the geocenter to the selenocenter; r m is the distance from the selenocenter to the geocentric; c s is the average distance from the geocenter to the heliocenter; r s is the distance from the geocentric to the heliocentric; Z m is the geocentric zenith distance from the Moon to the gravity station; Z s is the geocentric zenith distance from the sun to the gravity station; ϕ is the latitude of the gravity station; and ϕ is the geocentric latitude of the station. The Earth's atmospheric density decreases with an increasing altitude, such that gravity is affected by changes in the atmospheric quality and by deformation of the Earth [22,23]. The study area is located in the plateau area with an average elevation exceeding 3000 m; the elevations of some gravity points exceed 5000 m, which requires an atmospheric correction [24]: and where δ ga (10 −8 ms −2 ) is the air pressure correction; p (hPa) is the measured air pressure at the gravity station, which ranges from 630 hPa to 817 hPa in this survey; p n (10 −8 ms −2 ) is the standard pressure at the gravity station; and H (m) is the elevation of the gravity station.
To obtain the gravity values on the ground, the correction of each point is as follows [25]: where θ is the vertical gravity gradient (10 −8 s −2 ) measured at each point and h is the height of the gravimeter sensor with an original fixed tripod.

Scale Factor and Drift Correction
The scale factors can be calibrated using the baseline field or obtained as unknown parameters in the adjustment process [26] as follows: where ∆G 12 is the difference in the known gravity values between two gravity points and ∆g 12 is the average gravity difference measured by the relative gravimeter between two gravity points [27]. The method of using the adjustment model to estimate the scale factor is described in the adjustment model. In addition to calibration using the baseline field, the scale factors were estimated by the adjustment model. The drift rate of the relative gravimeter is calculated as follows: where g (10 −8 ms −2 ) and g are the round-trip observations at the starting point of the survey line and t and t are the corresponding moments of the round-trip measurement at the starting point of the survey line. The drift correction for the gravity observation can be calculated as follows: where δg k is the drift correction value, k is the drift rate, and ∆t is the time interval between the observation station and the starting point. In the above model, the drift rate of each relative gravimeter was calculated using the symmetrical observations. In addition, we also adopted the Bayesian approach, whose prior assumption is that the variation in the drift rate is smooth throughout the gravity survey for any relative gravimeter in good condition. The Bayesian approach applies ABIC to estimate the optimal variance of each relative gravimeter and the variances for variations of each drift, thus the drift rate of each time period can be obtained [1]. The simplex method can speed up the ABIC minimization process and improve the computational efficiency [28,29].

Classic Least Squares Adjustment
The error equation for the relative gravity observation can be written as: where g i and g j are the gravity adjustment values of points i and j, respectively; g RZi and g RZj are the preprocessed values of points i and j, respectively, after the solid-Earth tide correction, air pressure correction, instrument high correction, and drift correction; R i and R j are the readings of the gravimeter at points i and j, respectively, C K is the K-th grid value correction factor of the M-th degree polynomial lattice value function of the gravimeter; N is the number of period terms; X n and Y n are the periodic error parameters of the relative gravimeter; and T n is the period of the period error.
Remote Sens. 2022, 14, 1142 8 of 22 The weight of the relative gravity observations can be obtained using Helmert variance components estimation. If M denotes the number of relative gravimeters, the variance components, i.e., σ 2 = σ 2 1 · · · σ 2 M T can be calculated as follows [30]: and where P I (I = 1, 2, · · · , M) is the weight matrix for each relative gravimeter's observations, V I is the residual vector of each relative gravimeter calculated by the adjustment, N I is the normal equation of each relative gravimeter, N is the overall normal equation, and m I is the number of measurements from the I-th gravimeter. The steps in the Helmert variance components estimation are as follows. First, the relative gravity observations from each gravimeter were adjusted according to the equal weight treatment. The weights were then determined according to the posterior variance in the observations from each gravimeter, after which adjustment was performed again according to the posterior weights, repeating the above process until the unit weight medium error of the different observation types tends to be consistent [31]. The error equation for the a priori gravity value is as follows [8]: where V i is the residual of the gravity value at point i and g 0 i and g i are the priori gravity value with weight matrix P g and its adjustment value at point i, respectively. Combining Equations (13) and (16) yielded the following [8]: where L and V are the gravity measurements and the residuals vector, respectively, A is the design matrix, X is the unknown gravity vector at all stations. The unbiased estimation of the unknown parameters is: where P is the weight matrix, P = diag (P I , P g ). The calibrated scale factors of the baseline field are generally calculated before the adjustment of the gravity network. When calibrating scale factors in the survey area, certain phenomena may exist, such as insufficient gravity difference, excessively long survey line, and inconvenient traffic may exist. In this survey, reference gravity points were evenly distributed and there are a large number of redundant observations, which provides the feasibility to estimate reliable scale factors through adjustment model.

Robust Least Squares Adjustment
The robust least squares estimation can weaken the impact that gross errors have on the adjustment [32,33]. The error equation for the gravity observation with a weight of p i is as follows: The criteria for the robust estimation is [9]: Remote Sens. 2022, 14, 1142 , the equivalent weight of p i , which changes with the standardized residuals in the iterative calculation, is as follows [34]: The initial values ofX (0) , V (0) , andσ v i were obtained using the least squares method. The (k + 1)th solution of the robust estimation can be obtained after the determination of the (k + 1)th equivalent weight: The difference between two consecutive solutions of the same unknown parameter where ε is the threshold. Commonly used equivalent weight models include the Huber estimation [9], Danish method [35], IGG method [36], and Hampel three-tailed estimation [37]. The IGG III method was applied in this study, and the equivalent weight model was as follows [38]: where k 0 and k 1 are the thresholds and v i is the standardized residual. The IGG III method divides the observations into three parts: the observations with standardized residuals less than k 0 retain their original weights; the weights of observations with standardized residuals from k 0 to k 1 are reduced; and the observations with standardized residuals greater than k 1 are eliminated. The values of k 0 and k 1 are 1.5 and 4.5, respectively, in this study. The scale factors can be substituted into the error equation as unknown parameters for estimation.

Bayesian Adjustment
The Bayesian approach can estimate the variance of gravity observation noise and the variance of drift rate to determine the weights and drift rates of all gravimeters [1].
Assuming that the observation noise obeys a normal distribution and the variation in the drift rate is smooth, the relative gravity observation equation of gravimeter p, the absolute gravity observation equation, and the a priori assumption of drift rate of gravimeter p can be expressed by the following model: and where A p is the coefficient matrix corresponding to gravimeter p, D p is composed of the time interval between two adjacent gravity points for gravimeter p, v p is the drift rate vector for gravimeter p, y p is the vector derived from the observations of gravimeter p, W p is the weight matrix of gravimeter p, G is the coefficient matrix related to the absolute gravity observations, W g is the weight diagonal matrix of the absolute gravity, B p is a second-order smooth matrix, and W b,p is the weight matrix of the drift rate variations for gravimeter p [1].
Combining Equations (24) and (26) for all gravimeters with Equation (25), we obtained: where The ABIC serves as an additional constraint, which is expressed as follows [1]: where and H is the number of super parameters. The simplex optimization algorithm was used to increase the speed of the ABIC minimization in this study [28]. Each instrument observation noise variance and drift rate variance estimated by the ABIC minimization criterion are combined with the known absolute gravity observation variance to form the weight matrix, W [1]. The estimation of the unknown parameters can be expressed as follows: where x is the estimation of the gravity value, v is the estimation of the drift rate of all the gravimeters, and the scale factor of each gravimeter can be used as hyper-parameters to obtain its optimal value using the ABIC minimization criterion [39].

Analysis Procedures and Results
After the pre-processing of absolute and relative gravity data, the unknown parameters were estimated using the classic least squares method, robust least squares method and Bayesian approach. The characteristics of the different methods, used to determine the scale factors, drift rates, and weights, were analyzed. The applicability of the different methods for processing gravity conjunction data in plateau areas was analyzed.

Free Network Adjustment and Compatibility Test of Known Points
Free network adjustment evaluates the observation quality of the gravity network using the internal accord accuracy ( Figure 2). The results indicated that the gravity points with the highest accuracy were mostly distributed at the largest closed loop of the gravity network (B070, GUID, TODE, GCHA, and A080) or at the edge of the largest closed loop of the gravity network connected with multiple gravity points (QHTJ, B071, A078, XRID, B069, TREN, QHMY, JIAD, HNAN, QHDC, and BUQL). The gravity points located at the edge of the largest closed loop of the gravity network that merely connected with two gravity points indicated a slightly lower accuracy (A077, LENH, A079, JUZH, DARI, and MINH). Gravity points with the lowest accuracy were located at the annexed route or in the branch line of Qilian-Tuole (A082, YEGE, WUDL, and TUQL). F test was used to test the compatibility of known points by adjusting the gravity network into two groups, the first group is free network adjustment and the second group is the least squares adjustment involving known points. After adjustments of the two groups, the subsample variance of the first group was calculated as: (3) let significance level: α = 0.10, then F α = 1.14, F 0 < F α , and the null hypothesis holds, so the known gravity points are compatible.
first group is free network adjustment and the second group is the least squares adjustment involving known points. After adjustments of the two groups, the subsample variance of the first group was calculated as:

Scale Factor
In the adjustment where all of the reference points participate, the scale factors of each gravimeter were estimated by the classic least squares method, robust least squares method, and Bayesian approach. A slight difference existed in the scale factors for the same gravimeter when estimated by these methods, with a maximum difference of 1 × 10 −5 . However, the difference between the calibrated and estimated scale factors for gravimeters C259, C262, and C284 reached 5 × 10 −5 , the calibrated scale factor is less than one and the estimated scale factors were all greater than one for the C262 gravimeter ( Table 4). As such, further analysis was needed to discriminate the accuracy of the calibrated and estimated scale factors. The inaccuracy of the scale factors may lead to the existence of systematic errors in gravity differences because the gravity differences are multiplied by the corrected observations and the scale factor. For two gravimeters, if there is unequal error in their scale factors, mutual difference (i.e., the difference between the two gravity differences measured by them in the same section) will be correlated with the corresponding gravity difference. The mutual differences for two gravimeters with accurate scale factors follow a random distribution, which means they are independent of gravity differences.
When using the calibrated scale factors, a correlation was discovered between the mutual differences in the two gravimeters and gravity differences, and the slant linear fitting results indicated that the mutual differences relate to the value range of the gravity differences (Figure 3a,b). This phenomenon is related to the insufficient gravity differences in the baseline field. As the max gravity difference is about 600 mGal, and the gravity difference in the baseline field (430 mGal) could not cover the range of gravimeter readings in the survey area. The residual systematic errors may remain in the calibrated scale factors. When using the estimated scale factors, the correlation between the mutual differences and gravity differences was significantly lower; the range of the mutual differences became smaller, and the linear fitting result had an approximately zero slope (Figure 3a,b). Based on the random distribution of the mutual differences, it can be seen that the systematic errors were completely corrected.
The optimal accuracy of the estimated scale factors was obtained when all reference gravity points participated in the adjustment. Different numbers and distributions of the control points (i.e., reference gravity points used in the adjustment as initial-data) may affect the estimation of the scale factors (Table 5). When A080 (Xining, China) and A077 (Mangya, China), with the lowest gravity difference, were set as the control points, the scale factors estimated by the Bayesian approach were all less than one, significantly different than those controlled by all of the reference gravity points. The scale factors estimated by the classic least squares and robust least squares methods showed a difference of approximately 3.5 × 10 −4 compared to the scale factors controlled by all the reference gravity points. When A080 (Xining, China) and A082 (Qumalai, China), with the greatest gravity difference, were set as the control points, a relative higher difference occurred as compared with the Bayesian approach controlled by all of the reference gravity points; a scale factor difference of 2 × 10 −5 was observed between the least squares and robust least squares methods. When A080 (Xining, China), A082 (Qumalai, China), and A079 (Geermu, China) were set as the control points, the estimated scale factors showed a difference of approximately 1.5 × 10 −5 compared to those of the Bayesian approach controlled by all reference points; the differences were approximately 8 × 10 −6 and 1 × 10 −5 for the gravimeters using the classic least squares and robust least squares methods, respectively. The Bayesian approach was the most sensitive to the number of control points; the estimated parameters indicated a poorer accuracy when adjustments were controlled by a smaller number of control points. The classic least squares and robust least squares methods indicated no significant differences when the adjustment was controlled by the same points; an increase in either the number of control points or the gravity differences may improve the accuracy of scale factors.
The accuracy and reliability of the adjustment can be verified by some of the reference points not set as control points. Taking A080 (Xining, China) and A082 (Qumalai, China) as the control points, the network was adjusted by the least squares method using the scale factors listed in Table 5. The adjusted gravity values reached a discrepancy of 69.1 µGal on A079 using the scale factors calculated by A080 (Xining, China) and A077 (Mangya, China), while the discrepancy significantly decreased on A079, A078, and A077 when using the scale factors calculated by A080 (Xining) and A082 (Qumalai, China). The discrepancy decreased further when using the scale factors calculated by A080 (Xining), A082 (Qumalai, China), and A079 (Geermu, China). The adjustment had the optimal accuracy using the scale factors calculated at all of the reference points, with the largest discrepancy of 19.4 µGal ( Table 6). In summary, the gravity difference, quantity, and spatial distribution of the control points notably influenced the accuracy of the estimated scale factors, and the highest accuracy scale factors were estimated using the maximum number of control points, the greatest gravity differences, and the most widely distributed stations.

Drift Correction
We used the symmetrical observation mode and conducted measurements between gravity points A and C as follows: A − B − C · · · · · · C − B − A, and the drift rates can be calculated based on the field data surveyed by frequently returning to the starting points. After corrections for the solid-Earth tide, air pressure, instrument height, and drift, the corrected relative gravity difference was used in the classic least squares and the robust least squares adjustments. The Bayesian approach assumes that the variation in the drift rate varies as a nonlinear smooth function with the elapsed time, and the optimal values of the drift rate variances were estimated based on ABIC. Figure 4 showed the estimated drift rate of each gravimeter adjusted by the Bayesian approach controlled by all the reference points and that calculated by the symmetric observations. The drift rate estimated by the Bayesian approach gently changed with time, whereas the change in the drift rate was irregularly based on the symmetric observations. C259 and C262 showed the lowest drift rate less than 0.04 mGal/d, and the estimated drift rate difference between the Bayesian approach and that the symmetric observations can be nearly neglected (Figure 4a,b). The estimated drift rate of C283 by the Bayesian approach significantly differed to the symmetric observations, which reached 0.13 mGal/d on some days (Figure 4c).  The drift rate estimated by Bayesian adjustment controlled by reference point A080 was compared to the drift rate controlled by all reference points ( Figure 5), it can be seen that the maximum difference was less than 0.015 mGal/d on C259, C262, and C298 (Figure 5a,b,e), while C283, C284, and C617 indicated similar drift rates (Figure 5c,d,f). For gravimeters C284, C298, and C627, an identical time trend for the estimated drift rate was obtained by the Bayesian approach and the symmetric observations and the drift rate calculated by the symmetric observations significantly undulated over several days (Figure 4d-f).
The drift rate estimated by Bayesian adjustment controlled by reference point A080 was compared to the drift rate controlled by all reference points ( Figure 5), it can be seen that the maximum difference was less than 0.015 mGal/d on C259, C262, and C298 (Figure 5a,b,e), while C283, C284, and C617 indicated similar drift rates (Figure 5c,d,f).
The Bayesian adjustment simultaneously estimated the gravity values of the network and the drift rates. The classic least squares method obtained the gravity values after the observations were corrected by the drift rate calculated using the symmetric observations. The estimations of hyper parameters using A079, A082, B070, and A077 as the control points were similar to those using all of the reference points. When A079, A082, B070, and A077 were used as the control points and all gravity observations participated in the adjustment, the two types of discrepancies obtained by two methods, respectively, have negligible difference within 3.0 µGal (Table 7). Generally, there was negligible difference between the adjustment results of the two methods when all gravity observations participated in the adjustment. C283 showed a difference of 0.14 mGal/d in the drift rates between the Bayesian approach and that using the symmetric observations on some days ( Figure 6). The classic least squares method showed better performance than the Bayesian approach because the unsmooth variations at C283 were significantly inconsistent with the smooth hypothesis of the Bayesian approach (Table 8). The drift rate estimated by Bayesian adjustment controlled by reference point A080 was compared to the drift rate controlled by all reference points ( Figure 5), it can be seen that the maximum difference was less than 0.015 mGal/d on C259, C262, and C298 (Figure 5a,b,e), while C283, C284, and C617 indicated similar drift rates (Figure 5c,d,f). The Bayesian adjustment simultaneously estimated the gravity values of the network and the drift rates. The classic least squares method obtained the gravity values after the observations were corrected by the drift rate calculated using the symmetric observations. The estimations of hyper parameters using A079, A082, B070, and A077 as the control points were similar to those using all of the reference points. When A079, A082, B070, and A077 were used as the control points and all gravity observations participated in the adjustment, the two types of discrepancies obtained by two methods, respectively, have negligible difference within 3.0 μGal (Table 7). Generally, there was negligible difference between the adjustment results of the two methods when all gravity observations participated in the adjustment. C283 showed a difference of 0.14 mGal/d in the drift rates between the Bayesian approach and that using the symmetric observations on some days ( Figure 6). The classic least squares method showed better performance than the Bayesian approach because the unsmooth variations at C283 were significantly inconsistent with the smooth hypothesis of the Bayesian approach (Table 8). Table 7. Discrepancies between the adjustment gravity values obtained using different drift rates and the check points (all observations participated in the adjustment).

Weights
Gravity differences usually cannot be regarded as equal-precision observations, the variance component estimation was used to obtain the variance of each gravimeter given the nominal accuracy of Scintrex CG-6 gravimeter was better than Scintrex CG-5, and the results indicated that gravimeters C259 and C262 had smaller posterior variances (Table 9). Several high-precision gravimeters were involved in this survey and the gravity observations showed a high precision and no observation was eliminated when the robust least squares adjustment method was used. The residuals of the gravity differences obtained by the least squares method and robust least squares method all obey normal distribution (Figure 7), and the adjusted gravity values obtained by the least squares and the robust least squares methods indicated negligible numerical difference (within 2 µGal). For the convenience of discussion, parts of the reference points were used as control points and the others were used as checkpoints to test the reliability of the robust least squares method.
The observations of C259 and C262 were taken as examples for adjustment, and gross error (threefold the mean square error) was added in some observations of C259 (B069-WUDL, B069-A079, B069-XRID, and GUID-B070). We designed four groups of schemes to test the performance of the least squares method and robust least squares method (Table 10). In terms of the discrepancies between the adjusted gravity values and check points, the result of the robust least squares method is better than the least squares method (Figure 8). When A080 was used as the control point, the discrepancies between the adjusted gravity values and the check points obtained by the least squares method ranged from 39 to 113 µGal. Compared to the least squares method, the discrepancies obtained by the robust least squares method ranged from 9 to 41 µGal. When A080 and A078 were used as the control points, the discrepancies obtained by the two methods were significantly reduced compared to the result obtained by only one control point; however, the result of the least squares method still had a discrepancy higher than 50 µGal at A082, A079, and A077 ( Figure 8b). When A080, A078, and A079 were used as control points, the discrepancies at A077 near A079 decreased significantly, by 31 µGal for the least squares method and 8 µGal for the robust least squares method (Figure 8c). When A080, A078, A079, and B071 are used as control points, the accuracy of the two methods was improved slightly (Figure 8d). The adjustment results obtained by the robust least squares method are more accurate in the case of gross error.

C284
Scintrex CG-5 0.000240 1.66 C298 Scintrex CG-5 0.000336 1.19 C617 Scintrex CG-5 0.000399 1.00 Several high-precision gravimeters were involved in this survey and the gravity observations showed a high precision and no observation was eliminated when the robust least squares adjustment method was used. The residuals of the gravity differences obtained by the least squares method and robust least squares method all obey normal distribution (Figure 7), and the adjusted gravity values obtained by the least squares and the robust least squares methods indicated negligible numerical difference (within 2 µGal). For the convenience of discussion, parts of the reference points were used as control points and the others were used as checkpoints to test the reliability of the robust least squares method.  The observations of C259 and C262 were taken as examples for adjustment, and gross error (threefold the mean square error) was added in some observations of C259 (B069-WUDL, B069-A079, B069-XRID, and GUID-B070). We designed four groups of schemes to test the performance of the least squares method and robust least squares method (Table  10). In terms of the discrepancies between the adjusted gravity values and check points, the result of the robust least squares method is better than the least squares method (Figure 8). When A080 was used as the control point, the discrepancies between the adjusted gravity values and the check points obtained by the least squares method ranged from 39 to 113 µGal. Compared to the least squares method, the discrepancies obtained by the robust least squares method ranged from 9 to 41 µGal. When A080 and A078 were used as the control points, the discrepancies obtained by the two methods were significantly reduced compared to the result obtained by only one control point; however, the result of the least squares method still had a discrepancy higher than 50 µGal at A082, A079, and A077 (Figure 8b). When A080, A078, and A079 were used as control points, the discrepancies at A077 near A079 decreased significantly, by 31 µGal for the least squares method and 8 µGal for the robust least squares method (Figure 8c). When A080, A078, A079, and B071 are used as control points, the accuracy of the two methods was improved slightly (Figure 8d). The adjustment results obtained by the robust least squares method are more accurate in the case of gross error.

Schemes
Control Points Checkpoint a A080 A077, A078, A079, and A082 b A080 and A078 A077, A079, and A082 c A080, A078, and A079 A077, and A082 d A080, A078, A079, and B071 A077, and A082  The ABIC minimization can estimate the scale factor of each relative gravimeter with high accuracy, as well as estimating the optimal value of the observational noise variances to determine the weight of each gravimeter. Table 11 showed the observational noise variance of each gravimeter estimated by ABIC. Taking the weight estimated by ABIC as the initial weight, the robust adjustment weighted by ABIC yielded slightly more accurate adjustment values than those weighted by the variance component estimation (Table 12).

Discussion
Gravity networks located in plateau areas present challenging field survey conditions. In this study, six high-precision relative gravimeters were used for symmetric observations and eight reference gravity points distributed across the study area (Qinghai Province, China) were used as reference points. The data processing results using multiple methods indicate the sound accuracy of the network.
The scale factors played an important role in gravity data processing; multiple methods for the definition of each gravimeter's scale factor were applied in this study because the calibration was highly difficult owing to the large variations in elevation in this area.
Although the scale factor of each instrument was calibrated in the baseline prior to the gravity observation, the calibrated scale factors had poor suitability owing to the insufficient gravity difference in the calibration field. Systematic errors were found in the calibrated scale factors through a gravity difference in the simultaneous observation gravimeters and gravity difference in two adjacent stations. The estimated scale factors reduced the correlation between the mutual differences and gravity differences. The scale factors estimated by the three adjustment models were all more accurate than those calibrated by the baseline field with insufficient gravity differences.
The drift rates of the Bayesian approach and the least squares method using the symmetric observations showed similar accuracies when all observation stations were involved in the adjustment. This can be attributed to the high precision of the field survey and the nearly equal drift rates estimated by the Bayesian approach and by using symmetric observations. The variation in the drift rates at C283, which had poor stability, was not smooth owing to transportation, temperature, and other uncertainty factors arising from the plateau environment, which contradicted the prior assumptions of the Bayesian approach. When the observations of gravimeter C283 were merely involved in the classic least squares adjustment, the accuracy was significantly higher than that of the Bayesian approach, indicating that field observations are more accurate for complex drift undulation.
The equivalent weight of each observation calculated by the robust least squares method appeared to be more reasonable. This method weighted the observations with different precision by the same gravimeter, which reduced the effect of possible gross errors and yielded more reliable adjustment results.
The uncertainty of the terrestrial relative gravity observation mainly derived from the drift and scale factors of the gravimeters. All three adjustment models performed well in estimating the scale factors; a more reliable drift rate for the gravimeter was calculated using the symmetric observations. For gross error, the least squares adjustment result was significantly affected by gross error, whereas the robust least squares method performed best for gravity network adjustments in plateau areas.

Conclusions
This study applied the classic least squares, robust least squares, and Bayesian adjustment methods for gravity data processing in plateau areas, Qinghai Province, China. Using the different approaches, we estimated the scale factors, drift rates, and weights. The conclusions of this study can be summarized as follows:

1.
Scale factors calibrated by the baseline field have poor applicability in this gravity survey due to the insufficient gravity difference in the calibration field. When a large number of the control points with large gravity differences are across the study area, the scales factors estimated by the adjustment model also have high accuracy.

2.
The gravimeter drift rates estimated by the Bayesian approach showed a high accuracy for gravimeters with a smooth drift rate, while the drift rate calculated by the field symmetric observations indicated better performance for the gravimeters with low drift stability.

3.
The robust least squares method weighted by the ABIC with the scale factors as the unknown parameters showed the best performance, and the Bayesian and classic least squares methods showed similar adjustment accuracy.