Data Quality Assessment of Time-Variable Surface Microgravity Surveys in the Southeastern Tibetan Plateau

Abstract

Ground-based time-variable gravimetry with high accuracy is an important approach in monitoring geodynamic processes. The uncertainty of instruments including scale factor (SF) and drift rate are the primary factors affect the quality of observation data. Differing from the conventional gravity adjustment procedure, this study adopted the modified Bayesian gravity adjustment (MBGA) method, which accounts for the nonlinear drift rate, and where the SF is considered as one of the hyperparameters estimated using Akaike’s Bayesian information criterion. Based on the terrestrial time-variable gravity datasets (2018–2020) from the southeastern Tibetan Plateau, errors caused by nonlinear drift rate and SF were processed quantitatively through analysis of the gravity difference (GD) residuals and the mutual difference of the GD. Additionally, cross validation from absolute gravity (AG) values was also applied. Results suggest that: (1) the drift rate of relavive instruments show nonlinear characteristics, and owing to their different spring features, the drift rate of CG-5 is much larger than that of LCR-G gravimeters; (2) the average bias between the original and optimized SF of the CG-5 gravimeters is approximately 169 ppm, while that of the LCR-G is no more than 63 ppm; (3) comparison of the differences in gravity values (GV) suggests that the uncertainty caused by the nonlinear drift rate is smaller than that attributable to SF. Overall, the novel approach adopted in this study was found effective in removing errors, and shown to be adaptive and robust for large-scale hybrid surface gravity campaign which providing high accuracy gravity data for the geoscience research.

1. Introduction

Earth exhibits different forms of geodynamic phenomena stimulated by different physical processes such as surface fluid movement, earthquakes, material movement at internal boundaries, and interaction between earth layers. The high-precision time-variable gravity signal is the most basic and direct physical quantity reflecting the geodynamic characteristics of various environments in which the density of the medium changes [1,2]. Since the early part of the 21st century, earth observations from space by platforms such as the GRACE satellites have been developing continuously, and time-variable gravity signals acquired on the global scale for more than 10 years have been used widely in geoscience, hydrology, astronomy, and other related fields [3,4,5,6,7]. In comparison with satellite gravimetry, ground-based gravity measurements obtained at regular intervals at fixed stations on Earth’s surface are much closer to the field source, and the observed gravity signals can overcome the space resolution shortage (below 300 km) of gravity satellites. Terrestrial gravimetry is usually divided into relative gravimetry and absolute gravimetry [8,9,10,11,12,13]. For decades, with rapid improvement of different high-accuracy gravimeters and development of dense gravity survey networks, repeated surface gravity observations have been available to study geoscientific problems such as crustal deformation [14,15,16], groundwater change [17,18], oil, gas, and mineral exploration [19], glacier mass change [20], volcanoes [21,22,23], and earthquake precursor research [24,25,26,27,28,29,30].

A gravity signal obtained on Earth’s surface is a typical microgravity signal, and the uncertainties of such gravity data are complex, including both observational errors and field source factors of gravity variations. Observational errors require improved data adjustment theory that considers scale factor (SF), instrumental drift, and other factors. The field source factors of gravity variations require observation and elucidation of those factors that affect gravity signals, e.g., rainfall, groundwater, surface deformation, and air pressure. In fact, the former are the main factors, while the latter are complex and difficult to estimate quantitatively, and some errors are so small (<10 µGal) to be negligible [31]. Therefore, this study focused primarily on evaluating the quality of data in relation to the uncertainties of instruments.

Absolute gravimeters (e.g., the FG-5 and A-10 gravimeters) can provide direct gravity values with high accuracy; however, they tend to be too expensive, are generally cumbersome, and have rigorous requirements for the observation environment. Although relative instruments can only provide measurements of gravity variations, they are more economical, portable, and have minimal environmental requirements which making it suitable for large-scale field work. Moreover, depending on the control of absolute gravity (AG) values in the survey network, it is possible to obtain large-scale and high-precision terrestrial time-variable microgravity datasets. Currently, although various of high precision relative gravity instruments have emerged [8,9,10,11,12,13,32,33], portable spring-based relative gravimeters such as the Scintrex CG-5 gravimeter (CG-5) and the Lacoste & Romberg—G (LCR-G) gravimeter still play a key role in terrestrial gravity measurement. Generally, the zero-length spring is applied as the sensor for most relative gravimeters, and the measurement of the spring gravimeter is the feedback voltage, as it is proportional to the gravity change. The sensor of the CG-5 is a nonstatic fused quartz spring, which delivers accuracy better than 5 µGal (1 µGal = 1 × 10⁻⁸ m/s²) and reading resolution of 1 µGal [34]. In practice, nonlinear drift of a CG-5 relative gravimeter is inevitable due to the feature of quartz elastic, and it can reach up to 200 µGal/d [35,36], especially for a survey campaign with duration longer than 24 h [36]. Unlike the CG-5, the LCR-G gravimeter uses a metal spring as its sensor, which has less drift than a quartz system [35]. Additionally, another key factor is the main constant in relative gravimeters which called scale factor (SF) or calibration factor, and stress reduction in the spring of a gravimeter under the action of an alternating load will lead to a change of SF over time [37]. Generally, in China, SF calibration is performed every 3–5 years by the China Earthquake Administration on the long-baseline with known AG values. However, it is a time-consuming and expensive process, and the calibration period cannot meet the requirement for a twice-yearly observation survey. For some regional survey networks, the SF is calibrated on the short-baseline annually, which may be not accurate due to the insufficient gravity difference of the konwn AGs.

The classical gravity adjustment (CLS) method usually adopts the least squares method for the adjustment process [38,39,40], while the drift rate over the entire period of the campaign is considered linear. A number of different improvement methods have been proposed to remove errors attributable to SF bias [37,40,41] and instrumental drift [35,38,42]. However, the drift rate of the instruments is generally considered to be linear in the short period of loops (with few hours) or in the entire survey campaign, which is not accurate for the large-scale survey campaign with time span of months. Additionally, consideration of both SF calibration and nonlinear drift rate in the process of adjustment of relative gravimeters has rarely been discussed and applied in the large-scale continential gravity measurements. Different from the assumption of linear drift in the classical adjustment or piecewise zero drift based on the assumption of a smooth variation of the drift rate, Chen proposed the Bayesian gravity adjustment (BGA) method [34], in which Akaike’s Bayesian information criterion (ABIC) [43,44] was first applied to hybrid relative gravity survey adjustment. In this approach, the noise variance and drift rate variance are regarded as hyperparameters that can be estimated on the basis of the ABIC minimization criterion, and the unconditional marginal probability density function of absolute and relative gravity observations derived. Subsequently, Wang proposed the modified Bayesian gravity adjustment (MBGA) method [45], which adopts the SF as a third hyperparameter to be calculated.

In this study, using time-variable gravity data from the southeastern Tibetan Plateau acquired during 2018–2020, the effects of nonlinear drift rate and inaccurate SFs of different gravimeters were analyzed quantitatively using different gravity adjustment methods: the CLS method, BGA method, and MBGA method. Then, the GD (gravity difference between adjacent station pairs, which is the basic unit of relative gravity measurement and the basic element for the gravity adjustment) instead of the gravity measurements obtained at each station, were used as input information for the adjustment, and errors caused by nonlinear drift rate and SF were processed quantitatively through analysis of the GD residuals (difference between the observed GD and GD obtained from the adjustment method) and the mutual difference of GD (the difference of GD between every two gravimeter pairs). Moreover, Additionally, the cross-validation method was also proposed in the study to validate the optimized results.

2. Overview of the Study Area

The southeastern region of the Tibetan Plateau, located in the southern section of the North–South seismic belt of the Chinese mainland, is subject to intense compression and collision between the Indian plate and the Eurasian plate, making it an area of unique neotectonic movement with complex deformation and intense fracture activity [30,46,47]. Owing to the active faults in the area, a number of strong earthquakes have occurred since the beginning of the 20st century, e.g., the Ms7.8 Tonghai earthquake in 1970, Longling Ms7.4 earthquake in 1976, Gengma Ms7.6 earthquake in 1995, Lijiang Ms7.0 earthquake in 1996, Wenchuan Ms8.0 earthquake in 2008, Lushan Ms7.0 earthquake in 2013, and Jiuzhaigou Ms7.0 earthquake in 2017. Therefore, the study of terrestrial time-varying gravity data in this area is helpful for analyzing the variation characteristics of the density of the underground medium, reflecting its dynamic tectonic implications, and exploring the seismogenic structure.

Following the Wenchuan Ms8.0 earthquake in 2008, the China Earthquake Administration advanced the construction of its terrestrial gravity observation network to obtain high-precision and large-scale seismic precursor gravity variations before the occurrence of strong earthquakes. Thus, a hybrid gravity survey network in the southeastern Tibetan Plateau covering all major fault structures and under the control of AG stations has been established (Figure 1). The distribution of the survey network is not uniform owing to the topographic limitations. And there are approximately 19 absolute gravimetric stations and 537 relative gravimetric stations. No fewer than 500 gravity segments are repeatedly observed in each survey period. Each campaign has duration of approximately 90 d and spans approximately 12° of longitude and 13° of latitude. Generally, no fewer than two relative gravimeters are used simultaneously in a survey campaign: two CG-5 gravimeters (#1169 and #1170; hereafter, CG-169 and CG-170, respectively) are usually used in Yunnan and two LCR-G gravimeters (#149 and #150; hereafter, LCR-149, and LCR-150, respectively) are used in Sichuan.

The GD distribution of adjacent station pairs is shown in Figure 2. The maximum GD in Sichuan is approximately 390 mGal (10⁻⁵ m/s²), the maximum GD in Yunnan is approximately 305 mGal, the minimum GD of the entire survey network is <100 µGal, and the average GD is approximately 60 mGal. Generally, the distance between two adjacent stations is of the order of several tens of kilometers, and road vehicles are used to travel between stations to perform observations that typically have a time interval of several hours. To remove the influence of seasonal rainfall changes, survey campaigns are accomplished semiannually (generally starting in March and July each year).

3. Methods of Modified Bayesian Gravity Adjustment (MBGA)

Generally, more than two relative gravimeters are used in a gravity survey campaign that will contain redundant observations. To resolve this problem, the CLS method adopts the least squares method to establish an adjustment model in which the drift rate is considered linear. However, in an actual field gravity campaign, drift rate always shows nonlinear characteristics that reflect the duration of the survey period and the complexity of the geological structures, especially for CG-5 relative gravimeters. Moreover, the practical solution to the problem is attributed to an overdetermined system of equations, which leads to a nonunique solution. In this study, on the basis of BGA method [36], the modified Bayesian gravity network adjustment optimization model (MBGA) [45] was applied, in which the ABIC criterion was adopted to optimize the solutions, and it was assumed that the effects of ocean tides and pole shift were removed before the adjustment. Additionally, gravity changes caused by crustal deformation and precipitation were not considered during the same observation period.

Assuming that the observation error, absolute gravity error, and drift error of the gravimeter p obey the normal distribution, then:

$${{\displaystyle A}}_p\tilde{x}+{{\displaystyle D}}_p{{\displaystyle \tilde{v}}}_p-{{\displaystyle \tilde{y}}}_p~Normal(0,{{\displaystyle W}}_p^{-1})$$

(1)

$$G\tilde{x}-\tilde{g}~Normal(0,{{\displaystyle W}}_g^{-1})$$

(2)

$${{\displaystyle B}}_p{{\displaystyle \tilde{v}}}_p~Normal(0,{{\displaystyle W}}_{b,p}^{-1})$$

(3)

Equations (1) and (2) are the assumption used by the least-squares criterion [48] which are the observation equation for relative measurements and observation equation for absolute measuremts, respectively. Equation (3) is the calculation model of nonlinear drift which added in the BGA method, in which the drift rate of the relative gravimeter is assumed temporally smoothed.

Here, A_p, D_p, G, and B_p are the observation order matrix of GD, observation time matrix of GD, absolute point matrix, and second-order smoothing operator, respectively. Similarly, $\tilde{x}$, ${\tilde{y}}_p$, ${\tilde{v}}_p$, $\tilde{g}$ represent the gravity values at all stations, vector of the observed GD values, drift rate of all the gravimeters, and the observed AG values, respectively. In the above, Normal(.) represents the normal distribution with a mean value of 0 and variance of $W^{-1}$, and the weight matrix is expressed as: $W_p=diag\left({\sigma }_p^{-2},{\sigma }_p^{-2},\cdots ,{\sigma }_p^{-2}\right)$, $W_g=diag\left({\sigma }_g^{-2},{\sigma }_g^{-2},\cdots ,{\sigma }_g^{-2}\right)$, where ${\delta }_p^2$, and ${\delta }_g^2$ represent the observed noise variance and absolute gravity variance, respectively. $W$ and $W_b$ are regarded as the weights of the relative gravity observations and weights of the drift rate variations, respectively [45].

In Equation (1), GD of gravimeter p can be represented as:

$${{\displaystyle \tilde{y}}}_p={{\displaystyle \Delta y}}_p{{\displaystyle l}}_p-{{\displaystyle \Delta T}}_p\alpha -{{\displaystyle \Delta P}}_p\beta$$

(4)

where $\Delta y_p$ is the gravity reading data of GD, $l_p$ is the SF of gravimeter p, $\Delta T_p$ is the theoretical Earth tide, $\alpha$ is the earth tide factor, $\Delta P_p$ is the gravitational load difference of atmospheric pressure, and $\beta$ is the pressure admittance coefficient, respectively. In the process of gravity adjustment, the Earth tide factor [49] and pressure admittance coefficient [38,50] were considered as approximations of 1.16 and −0.3 µGal/hPa, respectively. Here, the $l_p$ is considered as a constant in the BGA method, while it is introduced as a new hyperparameter to be estimated in the MBGA method, together with the two hyperparameters (observed noise variance and drift rate variance).

Equations (1)–(3) can be combined as the full Bayesian modal for the adjustment:

$$SX-Y~Normal(0,{\tilde{W}}^{-1})$$

(5)

Here,

$$S=\left[\begin{array}{cc}A& D\\ G& 0\\ 0& B\end{array}\right],X=\left[\begin{array}{l}\tilde{x}\\ \tilde{v}\end{array}\right],Y=\left[\begin{array}{c}\tilde{y}\\ \tilde{g}\\ 0\end{array}\right],\tilde{W}=\left[\begin{array}{ccc}W& 0& 0\\ 0& {{\displaystyle W}}_g& 0\\ 0& 0& {{\displaystyle W}}_b\end{array}\right]$$

According to the ABIC principle, the posterior probability density function of the drift rate and gravity station values can be expressed as follows:

$$Posterior=p(v,x|g,y)=\frac{p(g,y|v,x)p(v)}{{\displaystyle \iint p(g,y|v,x)p(v)dvdx}}=\frac{p(y|v,x)p(g|x)p(v)}{{\displaystyle \iint p(y|v,x)p(g|x)dvdx}}$$

(6)

Here, because the solution is not unique, estimation of the model parameters is a typical ill-posed problem; therefore, we introduced the ABIC minimization criterion [44] as an additional constraint, and the optimal values of the parameters were estimated by fitting analysis using the Bayesian adjustment optimization model. Then, the ABIC formula can be written as presented in Equation (7) [36,45], which balances model parameter complexity and the ability of the model to describe the dataset (i.e., the likelihood function) by calculating the minimum value of the ABIC:

$$\begin{array}{ll}ABIC& =-2\times \log (\mathrm{maximum}\mathrm{likelihood})+2\times (\mathrm{No}.\mathrm{of}\mathrm{hyperparameters})\\ & =-\log \left|\begin{array}{cc}W& 0\\ 0& W_g\end{array}\right|+\log \left|S^T\tilde{W}S\right|-\log {\left|2\pi B^T{W}_{b}B\right|}_{+}+\min U(X)+2\times \dim (\theta )\end{array}$$

(7)

Here, $U\left(X\right)={(A\tilde{x}+D\tilde{v}-\tilde{y})}^{T}W\left(A\tilde{x}+D\tilde{v}-\tilde{y}\right)+{(Bv)}^T{W}_b\left(Bv\right)$, and $\theta =\left\{{\sigma }_1^2,{\sigma }_2^2,\cdots {\sigma }_p^2,{\sigma }_{b,1}^2,{\sigma }_{b,2}^2,\cdots {\sigma }_{b,p}^2,l_1,l_2,\cdots l_p\right\}$, where $\dim (\theta )$ is the dimension of $\theta$. The number of hyperparameters depends on the number of relative gravimeters, so hyperparameters is 3p in this study. And the hyperparemeters can be estimated through minimizing the ABIC by employing Newton method in this study.

4. Quality Assessment of Actual Gravity Datasets

The nonlinear drift characteristics of relative gravimeters and the variation of SFs are the two main errors in the processing of terrestrial time-varying microgravity data. The CLS, BGA, and MBGA methods were used separately to process the time-varying microgravity data obtained by the southeastern Tibetan Plateau gravity survey network during 2018–2020. The CLS method used the original SF and the linear drift rate, the BGA method used the original SF and considered the nonlinear drift rate, and the MBGA method considered both the nonlinear drift rate and the calibrated SF.

4.1. Drift Estimation of the Different Relative Gravimeters

In static and field observations, instruments exhibit temporal variation in the display of the zero position, which is called instrumental drift and represents one of the inherent features of change of relative gravimeters with time. On the basis of July 2020 gravity data, the daily drift rate of the different spring-based relative gravimeters is shown in Figure 3, fitted using the BGA method. Here, the drift rate is an average value calculated daily and expressed as the drift rate change per hour.

It is evident that all the instruments show nonlinear drift characteristics, although the amplitude of the drift rate between each gravimeter presents large deviation. The drift rate of the CG-5 gravimeters (up to approximately 37 µGal/h) is much higher than that of the LCR-G instruments (broadly below 1.5 µGal/h). It is also evident that the drift rate of the CG-5 instruments shows asynchronous difference, indicating that the CG-169 gravimeter is more sensitive to the observation route or the observation environment, and that the drift rate exhibits substantial nonlinear characteristics. The deviation between the maximum and minimum values of drift rate of CG-169 is approximately 10 µGal/h, i.e., the maximum drift difference will reach up to 240 µGal in a single day, which far exceeds the observational error in the measurement campaign.

Figure 4 shows the variation of GD residuals of the survey time series obtained using the CLS (Figure 4a) and BGA (Figure 4b) methods for different relative gravimeters. The amplitude of the GD residuals of the CG-5 gravimeters varies markedly, fluctuating in the range of ±30 and ±10 µGal for CG-169 and CG-170, respectively. The two CG-5 instruments also present obvious nonrandom variation characteristics when using the CLS method, and the variation trend is consistent with the trend of change of drift rate shown in Figure 3a. However, the GD residuals obtained using the BGA method present as a random white noise signal, which means that the variation of GD residuals has no obvious correlation with drift rate. In Figure 4e–h, the GD residuals of the LCR-G instruments are constrained to ±20 µGal, and because the drift rate amplitude of LCR-G instruments is negligibly small, there is no obvious difference between the GD residuals of the two LCR-G gravimeters. Nevertheless, we have found that the number of GD residuals (<5 µGal) from BGA method is slightly larger than that from CLS method in Figure 5. In addition, the SD of GD residuals decreased from 8.0 to 7.8 µGal and from 7.7 to 6.9 µGal for CG-5 and LCR-G gravimeter respectively, which means that the BGA method is more effective than the CLS method in removing errors caused by nonlinear drift.

The nonlinear instrument drift rate of the different relative gravimeters obtained using the BGA method during 2018–2020 is shown in Figure 6. All instruments were obviously affected by irregular nonlinear drift, and the general trends of the drift rate variation appear broadly consistent for the two instrument pairs in the same observation period. The nonlinear drift rates of the two LCR-G gravimeters change only slightly with low amplitude of approximately 1 µGal/h. However, the drift rate of the CG-5 instruments increases with observation time during the survey campaigns, while on the long term, the drift rate decreases gradually with increasing age and use of the instruments (Figure 6a,b). The mean drift rate of CG-169 and CG-170 stabilized at approximately 27 and 33 µGal, respectively, whereas the nonlinear drift rate of CG-169 is unstable and its amplitude changed drastically in comparison with that of CG-170.

Despite the decrease in drift rate of the two CG-5 m with age and use, the drift rate in the second half of the year is markedly higher than that in the first half of the same year. The nonlinear drift is possibly affected by spring age, temperature, transportation, and tilt [51], and relative gravimeters are likely to provide incorrect measurements of gravity change (of the order of tens of microgals) in conditions under which large temperature variations might occur [48]. Given that the difference in temperature between the two campaigns conducted in the same year might reach up to 20 °C, we attribute this phenomenon mainly to the notable temperature difference.

4.2. Uncertainty in Estimation of Scale Factor (SF)

The SF of gravimeters mainly varies with the range of the gravity readings of field-derived measurements; therefore, inaccuracy in the SF will lead to systematic error. Usually, instruments will be given the initial SF by the manufacturer, although in practice, the SF should be calibrated before each survey campaign using a baseline with known AG values, and where the AG interval on the baseline should not be less than the maximum GD of the network. In this study, the SFs were calibrated using the long-baseline at the beginning of 2018 of approximately 1.4 Gal, but were not recalibrated in the following years. The MBGA method adopted here can calibrate the SF without need of baseline calibration. In the next section, the influences of different methods and different AG values on SF calibration are evaluated quantitatively.

4.2.1. Calibration of Scale Factors

The MBGA method can use AG values to calibrate the SF, but it can also use one meter with a known SF to calibrate another meter used synchronously, which is especially useful for a survey network with only a few AG values.

• Scale factor calibrated by another gravimeter

For the example of survey data from March 2018, the CG-170 gravimeter was relatively stable, and its SF was calibrated using the long-baseline; the original SF of CG-169 and CG-170 was 1.000253 and 1.000066, respectively (ratio: 1.000187) (

Table 1. One gravimeter used to calibrate the scale factor of another.

Relative Gravimeters	Scale Factors
	Original 1	Calibrated by CG-170	Deviation (ppm 3)	Original 2	Calibrated by CG-170	Deviation (ppm)
CG-169	1.000253	1.000259	6	0.999956	0.999949	7
CG-170	1.000066	-	-	0.999752	-	-

¹ calibrated using the long-baseline; ² calibrated using the MBGA method; ³ ppm represents 10⁻⁶.

). Testing the SF of CG-169 calibrated by CG-170 revealed that the SF deviation of CG-169 was just 6 ppm, which could lead to a difference of <2 µGal in the maximum GD (300 mGal). Then, the SF of CG-170 optimized using the MBGA method was used to calibrate the SF of CG-169. It was found that even taking account of the different SF of CG-170 in the calibration, the SF ratio between the two instruments was almost the same. Therefore, it was considered that if the SF of one instrument is known exactly, the SF of another instrument could be calibrated using the known one, and in the same survey campaign, there is a fixed proportional relationship between the SFs of instrument pairs. However, if the SF of the first instrument is not accurate, it will introduce inaccuracy into the calibration result of a second instrument.

2.

Scale Factor Calibrated Using Different AG Values

Generally, calibration of SF was carried out in the baseline field before the gravity survey, and the AG values on the baseline are observed quasi-synchronously. However, AG measurements are usually completed only once a year, and owing to various uncontrollable reasons, the actual number of AG values might be less. Therefore, although relatively stable stations in a survey network can be selected as starting points for SF calibration, the different interval of AG values will affect the calibration results. In the following, taking the observation data of July 2020 in Yunnan as an example, the influences of different AG values on SF calibration are analyzed in combination with ABIC values.

Based on Equation (7), the noise variance, drift variance, and SF are taken as hyperparameters, and the estimation of hyperparameters is optimized by calculating the ABIC minimum value.

Table 2. Scale factors and corresponding Akaike’s Bayesian information criterion (ABIC) values by using the modified Bayesian gravity adjustment (MBGA) method and the Bayesian gravity adjustment (BGA) method with different absolute gravity (AG) values.

AG Stations Applied in the Calibration	Gravity Interval of AGs (mGal)	Scale Factors
		BGA			MBGA
		CG-169	CG-170	Values of ABIC	CG-169	CG-170	Values of ABIC
Case 1: All 10 stations	403	1.000335	0.999987	−7120	1.000147	0.999813	−7209
Case 2: Lijiang, Kunming Funing, Gengma, Ruili	364	1.000335	0.999987	−7094	1.000146	0.999813	−7174
Case 3: Lijiang, Funing	302	1.000335	0.999987	−7071	1.000122	0.999788	−7142
Case 4: Lijiang, Gengma	150	1.000335	0.999987	−7067	1.000087	0.999754	−7137
Case 5: Lijiang, Kunming	56	1.000335	0.999987	−7030	0.999881	0.999548	−7106

lists the calibrated SFs of two CG-5 instruments in the YN Survey network obtained using AG values from different stations. Obviously, the ABIC values obtained using the MBGA method are all less than those derived using the BGA method. When SF calibration is performed using the AG values of all the 10 stations (Xiangcheng, Lijiang, Dali, Panzhihua, Kunming, Ruili, Gengma, Simao, Mengzi, Funing) in the survey network (Case 1), because the GD of the AG values can cover the maximum GD (approximately 305 mGal) of the entire network of YN and because there is a sufficient number of stations participating in the adjustment, the ABIC value is sufficiently small, which means that the calibrated SF is accurate. When taking only Lijiang, Kunming, Funing, Gengma, and Ruili into the adjustment to calibrate the SF (Case 2), because the number of AG values involved in the adjustment is reduced, the ABIC value is increased slightly; however, the calibrated SFs of two instruments remain almost the same as in Case 1 as the AG value can cover the maximum GD of the network. The influence on SF calibration of using only two AG values with different gravity intervals was tested. Case 3: When Lijiang and Funing are used as the control stations to optimize the SFs, the gravity interval of the two AG stations is approximately 302 mGal, which broadly covers the maximum GD. Although the ABIC value increased owing to the reduced number of AG stations, the SF deviation compared with Case 1 is only 25 ppm, and the deviation of the maximum GD is only 7.5 µGal, which is within the uncertainty of measurement error. Case 4: Further reduction of the difference between the two AG values to approximately half that of the maximum GD in the survey area, the deviation is approximately 60 ppm, the effect on the maximum GD reaches to 18 µGal, and the corresponding ABIC value keeps increasing. Therefore, we think that such calibration results are not reliable. Case 5: Finally, when the GD of the AG values is 56 mGal (Case 5), which is only one fifth that of the maximum GD, the ABIC value increases notably, and the deviation of the SF reaches 260 ppm, making the average error up to 16 µGal.

The above results lead to the following conclusions. (1) When under the same conditon, the ABIC values obtained using the MBGA method are all less than those derived using the BGA method, which means that the SFs calibrated by MBGA method are more appropriate. (2) The optimal SFs were obtained when all the AG stations participated in the calibration, while when few AG values are available in a measurement survey campaign, no fewer than two AG stations are needed for SF calibration. (3) With increase of the gravity interval among the AG values, the ABIC value decreases, and SFs that are more accurate can be obtained. Therefore, the gravity interval of the AGs participated in the adjustment should be more larger to cover the maximum of the GD in the survey area. (4) When the number of AG stations involved in calibration in the same campaign remains the same, the ABIC value can be considered the criterion on which to judge the calibration results. (5) Irrespective of how the number of AG values changes in a survey campaign, the calibrated SF changes almost synchronously, and the ratio between the two instruments remains broadly the same. (6) The difference between the two AG stations involved in the calibration process should cover no less than half that of the maximum GD in the network, and the difference of the SFs are no more than 100 ppm.

4.2.2. Correlation Analysis of GD Mutual Differences

In practice, the accuracy of SF calibration can be analyzed from the mutual difference of the GD between two instruments. Normally, if the calibrated SF of two instruments is accurate, the mutual difference of the GD between each instrument pair will not increase with increase of the GD, that is, there is a weak correlation between GD and mutual difference of GD. Figure 7 shows the mutual difference of the GD of the two CG-5 relative gravimeters and the two LCR-G gravimeters in July 2020.

As shown in Figure 7a, the range of the mutual difference between the two CG-5 instruments is constrained to within ±30 µGal. Because the SF of this measurement period is the result of the long baseline calibration, the fitting line obtained using the MBGA method is just slightly better than that of the BGA method, and the slope of both fitting lines is broadly horizontal, i.e., the mutual differences do not increase notably with increase of the GD. In contrast, the range of the mutual difference between the LCR-G instruments is relatively small (Figure 7b), i.e., constrained to within ±15 µGal. However, it is evident that there is significant correlation between the GD and the mutual difference of the GD when using the BGA method, but that the correlation is significantly weaker when using the MBGA method. The SFs calibrated using the MBGA method are improved obviously in the gravity segment with a GD of >100 mGal, in which the deviation of the mutual difference obtained using the BGA method and the MBGA method can reach 15 µGal, and the range of the mutual difference has decreased from ±15 µGal to ±10 µGal. However, for the segments with small GD, the optimization effects of the MBGA method are not obvious.

Similarly, we also analyzed the GD residuals obtained using the BGA method (Figure 4b,d,f,h) and the MBGA method (Figure 8). It can be seen that because the nonlinear drift characteristics of the instruments have been taken into account in both methods, the GD residuals of the two results all show nonrandom variation characteristics, i.e., the changes in GD with time are broadly independent of nonlinear drift rate variations. For some gravity segments, the results obtained using the MBGA method can decrease the GD residuals. Owing to the small improvement in the correlation between the GD and mutual differences of the GD of the CG-5 instruments when using the MBGA method, the reduction of the GD residuals is not significant. However, for the LCR-G gravimeters, because the correlation has reduced obviously, the GD residuals have decreased.

Consideration of the histograms of the mutual differences (Figure 7c,d) reveals that the performance of the MBGA method is slightly better than that of the BGA method, especially for the LCR-G gravimeters. Therefore, from the perspective of the correlation between the mutual difference and the GD, the SFs calibrated using the MBGA method are considered relatively accurate (under the condition of a sufficient number of AG values participating in the calibration).

4.2.3. Variations of Calibrated Scale Factors

After testing the SF calibration methods, we adopted the MBGA method to recalibrate the SFs of each survey campaign conducted in the southeastern Tibetan Plateau during 2018–2020. Because the AG stations in the survey area are always observed only once a year, the SF calibration for March 2018 and July 2018 was based on the AG values observed in the first half of 2018, the SF calibration for March 2019, July 2019, and March 2020 was based on the AG values observed in the first half of 2019, and the SF calibration for July 2020 was based on the AG values observed in the second half of 2020.

The SFs estimated using the MBGA method were compared with those obtained using the BGA method. It should be noted that the original SFs obtained using the BGA method are usually actual calibrated or calibrated on the baseline. The SFs of the four instruments calibrated using the MBGA method show different deviations (

Table 3. Differences in scale factors obtained using the MBGA and BGA methods.

Date of Survey Campaign	Scale Factor of Different Gravimeters
	CG-169		CG-170		LCR-149		LCR-150
	Optimized (MBGA)	Difference (ppm)	Optimized (MBGA)	Difference (ppm)	Optimized (MBGA)	Difference (ppm)	Optimized (MBGA)	Difference (ppm)
March 2018	0.999956	297	0.999752	314	1.000047	51	0.988529	51
July 2018	0.999991	301	0.999738	328	0.999917	113	0.988404	109
March 2019	1.000126	37	0.999899	29	0.999959	62	0.988443	36
July 2019	1.000097	55	0.999835	29	0.999902	43	0.988383	31
March 2020	1.000094	139	0.999829	121	1.000060	14	0.989273	44
July 2020	1.000147	188	0.999813	174	0.999894	97	0.989136	92
Average		169		166		63		60

). The average difference of the CG-5 instruments is >160 ppm, and the maximum difference reaches 328 ppm, which could cause deviation of approximately 100 µGal on the maximum gravity segment with a value of 300 mGal. Thus, in comparison with the nonlinear drift variation of the instruments, it is evident that errors caused by SF bias might have greater impact on the accuracy of the adjustment results of an entire survey network. In comparison, the average deviation of the LCR-G instruments is no more than 65 ppm, and the difference on the maximum gravity segment (390 mGal) caused by the SF is no more than 30 µGal. As mentioned above, the differences in SF variations are related to the inherent features of the different spring-based gravimeters. In principle, the SF depends primarily on the range of the gravity reading and not on time; however, for the CG-5 instrument, the reading range changes slightly with time owing to long-term instrumental drift [38]. Therefore, the high drift rate variations of CG-5 gravimeters will cause larger deviation in the reading ranges, thereby affecting the difference of the calibrated SFs. It is also revealed that even though there is considerable difference in the SF deviation among different survey campaigns, the deviations between each gravimeter pairs are almost synchronous in the same campaign (Figure 9).

4.3. Cross-Validation by Using Absolute Gravity (AG) Values

Many factors can affect the drift rate of gravimeters, such as spring age, temperature inside the instrument, and transportation; therefore, using repeated observations is a practical approach when estimating the drift rate [34]. Additionally, the correlation between the mutual differences and the GD of segments can be seen as an indicator of whether the SFs of two pairs of instruments used simultaneously are preliminarily in accord. However, when the optimization of the mutual difference is not obvious, although the deviation of the SFs might exceed 150 ppm, the deviation is almost synchronous, which is not sufficient to evaluate the reliability of the SFs calibrated on the basis of the mutual differences of the GD. Therefore, high-precision AG values obtained quasi-synchronously by independent absolute gravimetry in the survey network were adopted to further verify the accuracy of the SFs and to revisit the drift rate of the gravimeters. The cross-validation method used some AG values for the gravity adjustment, while the remainder were used as validation points. Then, the calculated AG values were compared with the corresponding values obtained from the actual absolute observations to evaluate the quality of the adjustment results.

The survey data of July 2020 were analyzed and the gravity deviations of the AG values were calculated separately using the CLS, BGA, and MBGA methods. As shown in Figure 10, it clearly show that the deviation of the AG values calculated using the MBGA method are smaller than those derived using either the CLS method or the BGA method. For the case in which the AG stations of Kunming and Ya’an are introduced as the control stations participating in the adjustment (Figure 10a), and from the results of the CLS method, most deviations are no greater than 60 μGal, except for the stations of Ya’an, Songpan, Wanzhou, and Luzhou, which might reflect substantial influence by nonlinear drift, and there is no significant difference between the results of the BGA and MBGA methods when testing smaller numbers of starting AG stations. With further increase in the number of control stations (i.e., selecting stations Xiangcheng, Lijiang, Kunming, Panzhihua, Songpan, Ya’an, and Luzhou as the known stations), the cross-validation results (Figure 10b) show that the deviations obtained using the MBGA method decreases obviously; the average value is approximately 13 µGal, which means that the optimized AG deviations can reduce the gravity changes by approximately 20 µGal with the increasing number of known AG stations. Therefore, using the cross-validation of limited number of AG stations, reliable time-varying gravity data can be obtained by minimizing the uncertainty of the SFs and nonlinear drift.

To further investigate the spatial difference of gravity values (GVs) caused by drift rate and SF, the spatial distribution of the difference of GVs obtained using the MBGA, BGA, and CLS methods in the southeastern Tibetan Plateau are visualized in Figure 11. The spatial difference of the adjustment results from the CLS and BGA methods, displayed in Figure 11a, shows that the difference caused by nonlinear drift is small, i.e., within the range of ±12 µGal. In comparison, it is evident that the SF is the main factor affecting the gravity difference of GVs, with the maximum difference of approximately 60 µGal. The difference caused by drift rate is obvious in the larger loops or at the edge of the survey network, in which the survey time is long or where stations are not connected in a loop. The difference caused by SF is most obvious where the GD is large (e.g., in western parts of Pixian and Ya’an or to the south of Mengzi) or far from the constraint of AG stations, and this phenomenon is consistent with expectation.

Finally, through adjustment of the measurement data of the six campaign periods during 2018–2020 using the CLS and MBGA methods separately, the overall data quality was evaluated comprehensively. As shown in

Table 4. Adjustment results obtained using the CLS and MBGA methods for each campaign during 2018–2020.

Survey Date	Relative Gravimeters	Accuracy of GV (µGal)		SD of GV (µGal)		SD of GD Residuals (µGal)
		CLS	MBGA	CLS	MBGA	CLS	MBGA
March 2018	CG-5 (#1169, #1170) LCR-G (#149, #150, #132, #843)	6.8	5.9	2.2	1.3	9.3	8.2
July 2018	CG-5 (#1169, #1170, #834, #845) LCR-G (#149, #150)	9.0	7.8	3.2	2.8	12.5	8.7
March 2019	CG-5 (#1169, #1170, #834, #845) LCR-G (#149, #150)	9.3	7.0	2.8	1.8	11.3	8.5
July 2019	CG-5 (#1169, #1170, #834, #845) LCR-G (#149, #150)	8.4	8	2.1	1.8	10.5	9.6
March 2020	CG-5 (#1169, #1170) LCR-G (#149, #150)	7.8	6.9	1.7	2.3	10.7	9.2
July 2020	CG-5 (#1169, #1170) LCR-G (#149, #150)	6.1	5.7	1.6	1.3	9.3	7.7

GV: gravity values calculated by adjustment; SD: standard deviation; GD: gravity difference between adjacent station pairs. Even some other gravimeters (CG-834, CG-845, LCR-132, LCR-843) are used, however, owing to the survey time is about several days, we do not analyze them in this study.

, the accuracy of the GVs obtained using the MBGA method is better than that of the GVs derived using the CLS method, which is no more than 10 μGal, and the range of the GD residuals is 6.8–11.3 μGal.

5. Conclusions

Most geophysical phenomena are associated with mass transport, and gravity signals can be observed in ground-based gravity measurements. Among the various errors associated with spring-based gravimeters, the nonlinear drift rate and inaccurate SF are two of the most prominent. In practical, these can be estimated using redundant gravity observations, including increasing the number of AG values, increasing the frequency of the baseline calibration, and adopting optimized adjustment methods. In this study, we adopted the MBGA method to quantified the errors caused by nonlinear drift and SF bias. Through processing gravity data obtained by the southeastern Tibetan Plateau survey network during 2018–2020, the following conclusions were derived.

In comparison with the CLS method, which considers a linear drift rate, the drift rate of CG-5 and LCR-G relative gravimeters, used in the study, showed nonlinear drift characteristics when using either the MBGA method or the BGA method. Additionally, the amplitude of the drift rate showed obvious differences among the different instruments, i.e., the CG-5 gravimeter amplitude was much larger than that of the LCR-G gravimeters. Moreover, the nonlinear drift rate variations of the two CG-5 gravimeters were not synchronous. With the increase of service life, the drift rate of the CG-5 instruments decreased gradually and stabilized at 30–40 µGal/h. In contrast, the drift rate of the LCR-G was relatively small, i.e., generally no greater than 1.5 µGal/h.

By contrasting the GD residuals obtained using the CLS and MBGA methods, it was found that the GD residuals of the CLS method have obvious nonrandom signals, and that this change is synchronized with the variation of the drift rate. Conversely, the GD residuals obtained using the BGA and MBGA methods present as random white noise signals that have no significant correlation with the variation of drift rate. Thus, the gravity adjustment based on the Bayesian estimation principle can remove errors caused by the nonlinear drift of the instrument.

No fewer than two AG stations are needed for SF calibration, further, the gravity difference between the two AG stations should cover no less than half that of the maximum GD in the network. And when the number of AG stations involved in calibration remains the same, the ABIC value can be considered as the criterion to judge the calibration results.

Analysis of the correlation between the GD and mutual difference of the GD from the two gravimeter pairs revealed an obvious correlation from inaccurate SF, while the correlation obtained using the MBGA method was significantly weakened. Therefore, the optimized SFs from the MBGA method are more effective in reducing the system errors caused by inaccurate SFs. The adjustment results indicated that the average difference of the SFs of the LCR-G gravimeters was <63 ppm, while the average difference was >150 ppm for the CG-5 instruments owing to the different characteristics of the spring. Moreover, the variation of SF is also associated with instrumental drift.

Cross validation of the AG values revealed that the optimized SFs calculated using the MBGA method are better than those derived using the BGA method. Comparison of the difference of GVs among the different adjustment methods indicated that the uncertainty caused by the nonlinear drift rate (±12 µGal) is smaller than that associated with accurate SFs (−40 to +60 µGal).

In summary, under the constraint of AG values, the MBGA method can reduce the GV uncertainty caused by nonlinear instrumental drift and inaccurate SFs without calibration using the long-baseline, thereby allowing high-precision terrestrial microgravity measurements to be obtained, and lay a solid foundation for the application of time-varying gravity data in seismic monitoring and geoscience research.