Abstract
One of the most important operating procedures after the installation of a superconducting gravimeter (SG) is its calibration. The calibration process can identify and evaluate possible time variability in the scale factor and in the hardware anti-aliasing filter response. The SG installed in Cantley, Canada is calibrated using two absolute gravimeters and the data are analysed in the time and frequency domains to estimate the SG scale factor. In the time domain, we use the weighted linear regression method whereas in the frequency domain we use the least squares response method. Rigorous statistical procedures are applied to define data disturbances, outliers, and realistic data noise levels. Using data from JILA-2 and FG5-236 separately, the scale factor is estimated in the time and frequency domains as μGal/V and μGal/V, respectively. The relative accuracy in the time domain is 0.015%. We cannot identify any significant periodicity in the scale factor. The hardware anti-aliasing filter response is tested by injecting known waves into the control electronics of the system. Results show that the anti-aliasing filter response is stable and conforms to the global geodynamics project standards.
1. Introduction and Background
Nowadays, most systems and devices are controlled by digital computers where digital output is the norm. The superconducting gravimeter (SG) is no exception; its output is digitized by an A/D converter. The SG is a relative gravimeter and must be calibrated accurately to estimate its scale factor before further use of its records and data interpretation. The voltage output from the SG has no physical meaning until it is converted to the proper units of Gal by using specific calibration techniques. Therefore, physical experiment or externally induced acceleration can be used to represent or simulate the input or the physics (gravity changes) behind this voltage variation.
Reference [1] reported that a highly accurate scale factor is essential to constrain mantle anelasticity. In [2], it is mentioned that an accurate scale factor is needed, for ocean tide loading, removal of synthetic tide, and amplitude determination of the Earth's normal modes in the sub-seismic band. To meet the above-mentioned application requirements, it is very important to estimate the scale factor with relative accuracy of less than 0.1% [3]. However, it is really a challenge to calibrate a device using another one, which is much less precise: the absolute gravimeter (AG) precision is at the Gal level while the SG is at the nGal level. In other words, the SG precision is about three orders of magnitude higher than that of the AG.
2. Methods of SG Calibration
Various methods have been developed and implemented to calibrate the SG, specifically the amplitude calibration. One such method is called “relative method.” It uses another gravimeter, for example an absolute gravimeter or a relative gravimeter with a well-known scale factor. According to this method, both gravimeters have to run in parallel that is, at the same time and as close as possible to one another so they can sense the same gravity variations. Then, the least squares regression method is applied directly to both records to estimate the scale factor. The accuracy of this method is dependent on the data length and the precision of the gravimeters used in the experiment.
More often than not, the SG is calibrated using the relative method. Reference [3] ran the SG model TT70 and JILA-5 at the Strasbourg station for one day and estimated the calibration parameter with relative precision of less than 1.0% using linear regression. Later on, a compact SG (GWR C026) was installed at Strasbourg station. Reference [4] used over three years of AG and SG data to check the stability of the scale factor in the time domain. Reference [5] calibrated the SG at Membach, Belgium, using 47 days of data involving a large number of AG drops. A precision of about 0.1% was obtained through linear regression and tidal amplitude factors comparison. In addition to this high precision, the scale factor was checked in the frequency domain. Instead of running only one AG parallel to the SG, [6] used four FG5 AGs operating simultaneously for an average of a one-week period. The calibration factor was estimated from this experiment within 0.1%. Reference [7] also proved that this precision can be achieved with FG5 AG, recording for at least five days. The Japanese stations at Esashi, Canberra, and Matsushiro have been calibrated using the relative calibration method (e.g., [8–11]).
An alternative calibration method is known as the “absolute calibration” method. This method either uses a moving mass close to the SG or fixes the SG itself on a moving platform subjected to a known acceleration. The induced acceleration from the moving mass can be determined using Newton's law while the second method, which fixes the SG on a moving platform, uses the second time derivative of the vertical displacement of the platform. The accuracy of this method is limited to the accuracy of the gravitational constant, the homogeneity of the moving mass, and the geometry of the position of the mass with respect to the SG sphere [2]. Furthermore, the time span of the experiment is short; therefore, the regression of the induced acceleration and the output voltage from the SG has to be extrapolated to estimate the scale factor at zero frequency “DC” [12].
The SG absolute calibration method was first introduced by [13] by moving a sphere filled with mercury vertically under the SG; they achieved a precision of 0.2%. Reference [14] applied the same concept by moving a 273 kg annular mass around the SG and then used Newton's law to calculate the Newtonian attraction, which can easily be linked to the SG output; the scale factor was estimated with 0.3% relative accuracy. Reference [15] developed and applied the acceleration platform, which was used at the Frankfurt station. So far, the reported results from the platform provide a calibration precision of 0.01% [12]. This method is difficult to use because it requires special installations that are not easy nor are they available at Cantley station. Instead we used the relative method.
The SG at Cantley, Canada, was calibrated twice since the first installation by different researchers. Reference [16] estimated the scale factor with 0.4% precision by applying the theoretical tide method. Reference [17] estimated the scale factor to be Gal/V using JILA-2 gravimeter, with a relative precision of 0.13%. Hence, the SG at any station must continuously be calibrated whenever data are available using different methods. In order to contribute to the previous calibration effort of the SG at Cantley station, we recalculated the SG scale factor based on long series from the superconducting and absolute gravimeters from 1997 to 2009. The calibration was carried out in the time and frequency domains. To the best of our knowledge, the Least Squares Response Method (LSRM) has never been used to calibrate the SG in the frequency domain.
The least squares response method (LSRM) developed in [18], which will be revisited again in the methodology section, is used to estimate the scale factor in the frequency domain. The LS spectra of both SG and AG are used to define common peaks. These common peaks are identified through the product spectrum of the two factor spectra. Then, the frequencies of the statistically significant peaks at 95% CL are used to estimate the amplitude and phase by suppressing its period in the time series (AG and SG). Finally, at each frequency the amplitude and phase difference of the scale factor is estimated, and then the least squares fitting method is used to estimate the scale factor in different frequency bands.
In addition to the frequency domain analysis mentioned above, we estimate the scale factor in the time domain by separating the gravity data into blocks based on the duration of the experiment. Then, the scale factor of each block is estimated using the weighted least squares regression and its value corresponds to the middle of the block. The weighted average of all individual scale factors and its standard error are finally calculated.
3. Methodology
3.1. Least Squares Spectral Analysis
We use the Least Squares Spectral Analysis (LSSA) to estimate the spectra of the absolute and superconducting gravimeter series and subsequently produce their product spectrum. The LSSA has its roots in [19, 20]. Its advantages over classical Fourier analysis have already been presented in the literature and need not be repeated here (see, e.g., [21–23]). For the sake of comprehensiveness, we only present the fundamental formulas and emphasize the statistical properties of this method.
Consider a time series observed at discrete times , , not necessarily evenly spaced, which is essentially equivalent to the presence of gaps in the time series. This time series also has an associated covariance matrix that carries information on the uncertainty of the observed values as well as on their correlation (usually not available). The LSSA spectrum is described by the percentage variance of the spectral content of a specific cyclic frequency , which is the ratio of the quadratic norm of the specific signal to the total quadratic norm of the series: where is the projection of onto the model space spanned by the trigonometric base functions of cyclic frequency ; . For cyclic frequency , where and are estimated coefficients from the least squares fitting. For more details about the geometrical interpretation of the least squares transform and functional analysis see, for example, [24, 25]. Reference [22] showed that the probability density function (PDF) of the LS spectrum is a distribution defined by two parameters and , where and depends on the number of data points (length of the series) and the number of unknown parameters estimated by the LS estimation (degrees of freedom of the LS system).
3.2. Least Squares Product Spectrum
The LS product spectrum has been used as an effective and rigorous approach to identify common spectral peaks among the time series obtained from superconducting gravimeters for the purpose of detecting very weak signals buried in noise, such as the Slichter triplet [26]. In our case, two different time series (AG and SG) are used to obtain their product spectrum. It should be pointed out that both series do not necessarily need to have the same sampling rate, only their corresponding spectra must be of the same band and resolution (number of spectral values in the band) to make possible the mutual multiplication. This procedure can be applied to two or more time series to define the common features (multichannel system or multi-input single-output system).
The PDF of the product spectrum of two series can be derived using standard statistical approaches (e.g., [27]). This PDF is very important for the definition of the significance of the common peaks at different confidence levels. Moreover, the mean and the variance of the product spectrum can easily be defined. Let a random variable denote the LS spectrum, which is distributed as , with and dependent upon the degrees of freedom of the linear LS system as follows [22, Theorem 3]: where is the number of data points in the time series and is the number of unknown parameters that are estimated simultaneously with the LS spectrum. This distribution has finite first and second moments [22]:
For the derivation of the PDF of the product LS spectrum () of two random variables and , that is, the derivation of the PDF of the random variable , we follow the same procedure as in [26] by deriving the PDF of the new random variable . However, the derivation of the PDF in [26] was applied for random variables, where was larger than 10, which resulted in a PDF approximately following the Gaussian distribution. In this paper, we deal with two random variables both having the same distribution but different parameters. These two random variables are the percentage variances of both AG and SG series. First, we take its natural logarithm so the new random variable becomes
Before deriving the PDF of , we must project the PDF of to . If is a random variable of continuous type, such that [26] in the set , we must find the PDF of the new random variable . We introduce the transformation , which gives . We now have the mapping of into , and the Jacobian of transformation is . The PDF of is then , or , with and , which can be obtained by applying the well-known definition of the mean and the variance of PDF, see, for example, [27]: where is Euler's constant, is the logarithmic factorial function, and is the th derivative of [28]. Now, the mean and the variance of the new random variable can be defined as follows: the mean is the sum of the two individual means and the variance is the sum of the two variances (e.g., [27]):
We can now obtain the PDF of the sum of the natural logarithms of the two spectra , where and are in the set and is in the set , by where is the integration variable. Since is the sum of two independent random variables and , both following , distribution, then follows the same distribution but with different parameters [27]. So, the general formula of is . From the PDF characteristics (e.g., and ), two out of the four unknown parameters ( and ) can be defined analytically giving the final PDF: where and are two parameters dependent on and , both of which are smaller than unity and can be estimated from the numerical integration of (7). Notice that is the function. The above PDF underlines the product LS spectrum and can be used to identify the significance of product peaks of both series at any confidence level (e.g., 90%, 95%, or 99%).
3.3. Least Squares Response Method
After producing the product LS spectrum of AG and SG records and using its PDF as a guide, we can identify statistically significant peaks at any confidence level. In this paper, we use the 95% confidence level to identify significant peaks. These peaks (actually their periods) are suppressed separately in the AG and SG series to estimate their amplitudes and phases using weighted least squares as follows: where and are the offsets or averages in the AG and SG records, respectively, and are the parameters of AG and SG linear trends, if any, is the number of common significant peaks, and are the amplitudes of AG and SG constituents for the th frequency, respectively, and , are their associated phases. All the above parameters have an associated covariance matrix estimated from the LS estimation process. Only the statistically significant amplitudes and phases are used to estimate the scale factor. The magnitude and phase difference of the scale factor is then estimated for each th frequency from
By applying the covariance law to (11) and (12), we can calculate the standard deviations and of the scale factor and phase difference, respectively:
We note here that the scale factor is calculated only at specific cyclic frequencies at which both AG and SG data have statistically significant amplitudes.
4. SG and AG Data at Cantley Station
The SG at Cantley station started recording gravity data on November 7, 1989. Shortly after the installation the SG was drifting excessively (0.5 Gal per day), and several years later a new gravimeter sensor was installed and started recording in January 1997 and contributed to GGP during phase I and II. The SG has been operating on the pier of the Canadian Absolute Gravity Site (CAGS) close to JILA-2 and most recently FG5-236 Absolute Gravimeter. The JILA-2 gravimeter contributed calibration data until January 2007, and then it was replaced by FG5-236 in May 2007. In this work, we use AG data from both JILA-2 and FG5-236 for the calibration. Having the AG and SG running in parallel provides a good opportunity to examine the drift of the SG and check the signals coming out of both [17].
Over a period of 13 years (1997–2009), there have been 180 AG calibration experiments. The experiment durations range from a few to 200 hours with 7200 h of total length. The observation periods were most likely synchronised with the maximum Earth tide at Cantley station. Also, it is worth mentioning that the AG observations are unevenly spaced but the average sampling interval is very close to 60 s, which is not a problem when applying the LSSA method. Figure 1 shows the data availability and their duration per calibration experiment. Over 0.52 million sets are available from the AG. These records are distributed along the entire 13-year period and contaminated by all kinds of noises, for example, earthquakes, snow storms, and others.
5. Data Preprocessing and Filtering
Data preprocessing is one of the most important steps for obtaining a high precision scale factor. Three major factors or parameters have to be checked and tested before we estimate the scale factor namely, data offsets, data outliers, and the variance of the original data. The SG data are not free from small or large offsets, which need to be identified before the calibration process. Large offsets can be defined easily even by visual inspection but small offsets are not. To identify small offsets, all known signals must be removed first and then check the significance of the existing offsets as well as the residual outliers. Also, with this approach, the realistic noise level of the residuals can be estimated. The procedure below is followed when dealing with the SG data.(1)For each AG block or experiment, the corresponding SG data segment is retrieved but extended 15 s beyond each end of the AG record to accommodate the filtering process. However, if there are no SG data corresponding to the AG records for any reason, the AG record is discarded.(2)A nominal scale factor (from previous calibration) is used to transform the SG data from volts to Gal. Then, the Earth tide, ocean loading, pressure, and polar motion are removed from the raw SG data. Earth and ocean tides are removed based on the synthetic tide generated from “GWAVE” [29] while the pole tide is computed using the daily polar motion data from the International Earth Rotation and Reference Systems Service (IERS). A constant atmospheric admittance (−0.3 Gal/mbar) is used to model the pressure effect. Hence, the resulting residual series is practically free from most of the known signals.(3)A Parzen window function with lag = 15 s (cutoff frequency Hz and bandwidth Hz) is used to remove the high frequency component and provide a measure of the noise level. This is achieved by propagating the variance of the SG residual series segment within the window to the filtered value [23, 26]. Each filtered value is estimated exactly at the time AG measurements are available. In other words, the Parzen window is centered at the same time marks as the AG measurements.(4)For each SG segment corresponding to the AG experiment, the LSSA v5.02 software [22, 23] is used for the following purposes: (a) offset detection and identification of the time they occur, (b) residual outlier detection and flagging, (c) detection of any periodicities or deterministic signals followed by their suppression, and (d) test on the variance factor using Chi-square test [24].(5)Based on the previous step, (a) the times of statistically significant offsets are identified for later analysis, (b) outliers more than 4.215 (95.0% CL for in-context test with 2000 data points [24]) are rejected. Their corresponding values in the AG record are excluded as well, and (c) the standard deviation estimated in step (3), which is associated with the filtered value, is rescaled (adjusted) using the a posteriori variance factor when the Chi-Square test fails in step (4d).(6)Finally, all removed signals from steps (2) and (4c) above are restored in the residual series and translated back into volts using the same scale factor as in step (2) as well the standard deviation estimated in step (5c).
The AG data are provided to us as follows: time stamp, gravity set, and standard deviation. The standard deviation of each gravity value is a formal estimate. The standard deviation is between 1.0 and 3.0 Gal, which is rather too optimistic. Hence, the procedure below is followed to calibrate the standard deviation and outlier detection.(1)Similar to step (2) in the SG data preprocessing, all known signals, such as Earth and ocean tide, atmospheric pressure, and polar motion are removed from the AG records.(2)At each data point, we retrieve all AG data that extend 5.0 min at either side to estimate the standard deviation of the gravity value at the selected point. Then, the mean and the trend of this small data segment are subtracted and the estimated standard deviation of the residual values in the segment replaces the original value.(3)Residual outliers are rejected along with their corresponding values in the SG record. The Chi-Square tests on the variance factor for all segments passed. This indicates that the standard deviation estimated in step 2 represents a realistic noise level of the AG records. (4)All removed signals are restored in the residual AG series, as with the SG.
Figure 2 shows a segment from AG and SG before and after calibrating the standard deviation. The top-left panel (a) shows the AG with an uncalibrated error (provided with the original data), which is very optimistic, while the top-right panel (b) shows the AG after removing all signals with the calibrated standard deviation (gray error bars) based on the realistic noise level. The lower panels of Figure 2 represent the SG data in volts.
(a)
(b)
(c)
(d)
6. Scale Factor Determination
6.1. Scale Factor in the Frequency Domain
Due to the large number of AG data points (~0.52 million), a lowpass filter is applied to decimate the data (AG and SG) and reduce the high frequency noise. A Parzen window with lag = 450 s ( mHz and mHz) is applied to produce data points with an average sampling interval of 15 min. The LS spectrum in the band 6.00000–0.00025 cpd (0.25–4015.00 d) of the filtered data is estimated using the LSSA v5.02 software. The identification of all significant peaks in the product LS spectrum at 95% CL is done rigorously using the probability density function (PDF) of the product LS spectrum defined by (8) (dashed lines in Figure 3). In fact there are many peaks that are significant, but they are close to each other. When these significant peaks are suppressed, we find that most of the estimated amplitudes and phases are highly correlated (from the variance-covariance matrix). Due to the high noise level in the AG data and the huge number of tidal constituents especially in the diurnal, semidiurnal and terdiurnal bands, it becomes difficult to resolve these peaks. Thus, we filter out these peaks and we use only the ones that satisfy the following two criteria: (1) highest peak and (2) smallest correlation with others.
(a)
(b)
(c)
(d)
(e)
(f)
Next, all significant peaks identified in the previous step are suppressed using (9) and (10) to estimate their amplitudes and phases simultaneously with other systematic noise (e.g., trend and offsets). We use (11)–(14) to estimate the magnitude and phase difference of the scale factor associated with the covariances. We note here that the scale factor is calculated only at specific frequencies at which both AG and SG time series have significant amplitudes. Then, the weighted LS regression is used to define the best fit to all admittance points to obtain a smooth response function [18].
The estimated amplitude of the scale factor in the frequency domain is Gal/V, while the phase difference is . The relative accuracy is 0.096%, which is higher than that in the time domain. The phase difference (negative sign for all the SGs in GGP network) indicates how the feedback system in SG responds to gravity variation. Our estimation of the phase difference indicates that SG lags behind AG by , which is not significant. To check the scale variability we select five bands (ter-diurnal (TD), semidiurnal (SD), diurnal (D), long period (LP), which includes Mf, Mm, and Mtm, and other seasonal (Y), which includes Ssa, Sa, and others) to estimate the scale factor amplitude and phase difference at the midpoint of these bands. The red diamonds in Figure 4 and Table 1 represent these values.
6.2. Scale Factor in the Time Domain
All AG and their corresponding SG calibration experiments are used to estimate the scale factor in the time domain. For each individual experiment block the scale factor (SF) can be estimated by fitting a straight line through AG and SG residuals by minimizing the square of the weighted residuals (in the least squares sense) as follows: where AG and SG are, respectively, the AG and the SG data after subtracting the mean value. The estimated scale factor and its associated standard error correspond to the middle of the block (see Figure 5(c)). Figure 5(a) shows one of this regression lines for one calibration experiment. Then, the least squares weighted average is applied to determine the scale factor based on 11 years of data from JILA-2 and two years from FG5-236. Using only JILA-2 data, the SG scale factor is estimated in the time domain as Gal/V. The FG5-236 short time series (May 2007–February 2009) is used to reestimate the scale factor in the time domain to be Gal/V, which is statistically compatible with the one estimated from JILA-2 data. Combining both datasets the final scale factor becomes Gal/V with a relative accuracy of 0.015%, which is comparable with the absolute method of calibration of SG when using a platform. This is the first time to achieve such high precision for the scale factor due to the fact that there are very long AG records for calibration. Further analysis of the scale factor in the time domain reveals that there is no statistically significant drift. However, the LSSA is applied to investigate the existence of any periodicity in the scale factor. The frequency domain representation of this series shows no evidence of any periodicities in the scale factor as there is no single peak greater than the 95% CL (Figure 5(b)).
(a)
(b)
(c)
7. Calibration of SG Electronics
The SG electronic package includes electronics that sense and control all subsystems (e.g., gravimeter, tilt, and temperature). An important feature of the gravimeter control is the gravity circuit card, which is mounted in the analog control chassis. On this card, there is an eight-pole Bessel hardware filter with a corner frequency of Hz. This filter, known as GGP1, is used as antialiasing hardware filter for digitizing the continuous gravity signal to 1 s sampling interval with a Nyquist frequency Hz. The GGP1 filter conforms with recommendations put forth by the GGP. The red line in Figure 6 shows the theoretical filter response in the frequency domain represented by the gain and phase response.
(a)
(b)
This analog Bessel filter is characterized by almost constant group delay across the entire passband (0–50 mHz), thus preserving the wave shape of filtered signals. In addition, the filter does not attenuate the signals up to the frequency of 2.8 mHz (as seen from Figure 6). However, this has to be tested frequently by injecting a series of sine and/or square waves to the sensor electronics and study its actual response [30].
In January 22, 2008 the original filter (provided with the instrument) at Cantley station was replaced by the GGP1 filter. Six months later, the transfer function (gain and phase) of this filter was tested by injecting 14 sine and square waves with periods ranging from 3 s to 200 s and amplitudes ranging from 0.013 V to 0.6386 V, where the lowest amplitude is equivalent to 1 Gal. The black dots in Figures 6(a) and 6(b), show the gain and phase response, respectively.
It is obvious from Figure 6(a), that any noise at frequencies higher than 0.160 Hz (<6 s) is filtered out totally while the real cutoff frequency from this experiment is 0.070 Hz (14.3 s), which is close to the theoretical Hz. Also, the filter preserves the signals at frequencies lower than 0.005 Hz (200 s). Based on this test (as seen in Figure 6(b)), the estimated phase delay at Hz is 14.4 s, which is close to the specified filter parameters required by GGP.
8. Summary and Concluding Remarks
In this paper, the scale factor of the Cantley SG was determined in both time and frequency domains. The time domain estimation gave better accuracy ( Gal/V) than the frequency domain ( Gal/V). A couple of reasons may have contributed to this result: (1) the frequency domain estimation considers the noise level of both gravimeters SG and AG, while the time domain calculation deals only with AG noise; (2) using the whole dataset may hide real existing offsets and reduce the accuracy level in the frequency domain, while in the time domain the data are already segmented into blocks. However, these two estimates are statistically compatible (the compatibility of two means was statistically tested using -test at 95% CL).
In the frequency domain, the scale factor in diurnal and semidiurnal bands is more accurate than in the other bands as these two bands have the strongest Earth tide signals. In addition, the scale factor determined in the time domain is closer to diurnal and semidiurnal bands estimation although not commensurate with the accuracy levels.
Finally, we could not identify any drift or periodicities in the scale factor. Also, the new installed GGP filter was tested to conform to the published specifications. Last but not least, we recommend the scale factor posted at the GGP website be changed to the new and accurate one estimated in this research from the time domain.
Acknowledgments
This research was financially supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada, Ryerson University Post-Doctoral Fellowship, and the Start-Up Grant from University of Cape Town for the first author. The authors thank the manager of Cantley, Canada station for making the SG and AG data available for this research.