|COMPUTING THE SEISMIC
David Crossley, Saint
Step 1. Assemble in a directory all the 1 sec (5 sec or 10 sec) raw gravity and pressure files for each day of the sequence to be examined. It is important to have daily files for his analysis. It may be convenient to use 1 year at a time (about 365 days).
These files must not be corrected in any way for spikes, gaps earthquakes or other disturbances. You must have the volts to gravity calibration factor for the instrument.
Step 2. Fix the pressure files at the original sampling interval for spikes, gaps, and offsets so that problems in the pressure do not get transferred into the gravity data at the next step.
Step 3. Subtract from each day of gravity data the following:
(a) a synthetic elastic tide, based on a modern tidal potential (Tamura, Qinwen, or later) with recent values for the elastic tidal Love numbers. We chose the Qinwen potential with a cut-off of 0.0001, yielding 383 waves, based on the program GWAVE by Merriam. It does not make any difference whether the FCN correction is used or not.
Although in principal ocean tides can be added to the potential (e.g. through local tidal delta factors), this will be station dependent. Banka and Crossley found that the inclusion of ocean tides does not affect the noise levels in the seismic band (nor in fact does the use of a highly accurate tidal potential with many more waves than 383).
(b) atmospheric pressure with an admittance of -0.3 microgal per millibar. Again some better admittance may be chosen but will not affect the outcome significantly.
Step 4. Subtract from each day a 9th degree polynomial to reduce the amplitude of residual tides.
Step 5. For each of the days of reduced data, subtract the mean value and compute the RMS deviation of the data in microgal. Select the 5 days with the lowest RMS deviations as representing the quietest days seismically of that year.
Step 6. For each of the 5 quiet days, find the FFT computed as follows:
(a) use a 10% Cosine Bell window to equalize any jump between the beginning and end of the data (this would affect the wrap around in the FFT). Compensate for the loss of time domain power by multiplying all remaining values by a correction factor (1.06667 for the 10% cosine bell).
(b) pad the data with zeros to the (next+1) power of 2. For example, 86400 data points would be padded to 262114 or 218.
(c) compute the FFT and save the amplitude spectrum, unnormalized.
Step 7. Do a simple average of the 5 spectra files, i.e. for each frequency take the mean of the amplitudes from each file. This mean FFT file is to be plotted as a power spectrum as part of the SNM analysis. For comparisons with other spectra the axes should be in powers of 10 between 0.01 and 100 (mHz) in frequency and between 10-11 and 10-1 in power (microgal2 per mHz).
Step 8. Compute the mean Power Spectral Density from the mean FFT file in the period range 200-600 sec (this should avoid the sphere resonance of all the gravimeters we know about to date). This mean PSD can be converted to the final SNM value using:
SNM = log10(mean PSD) + 2.5
PART II - PROGRAMS
Some data and programs that have been used for the Boulder data and for the Cantley data are shown below and may be freely copied and used if found convenient.
1. Sample files are:
These are the 5 quietest days for 1990. The pressure files did not need any fixing for these days. Note these are not in GGP format because GGP failed to come to an agreement about standardizing the format of raw data.
2. Subtracting Tides and Pressure
I used the program g1sres.f to subtract tides and pressure and convert between volts and microgal. The tidal potential is that of Xi Qinwen, contained in the file qinwen.dat. The reduced data are in the *.tps files:
They are in a format called PLOTP (from a program called PLOTP for displaying gravity data that I wrote in 1990) that is much more compact than GGP or ETERNA.
3. Subtracting the Polynomial and Analyzing the RMS
I use a program gnoise.f that reads in the *.tps files, subtract the polynomial, saves the output and gives the ordered list of days with the lowest RMS. The files produced are *.dmp (see above). The summary file is gnoise.out (renamed to gnca90.out).
4. Computing the FFT
I use a program called period.f to do this step. Although the I/O is manual it goes quickly for only 5 days (it could obviously be modified to run in batch mode). The output files are (watch out, these are large files!)
5. Computing the mean FFT.
6. Producing the SNM.
The final step is to plot the mean .fft file as a power spectrum. I use PLOTP to do this (you can have this by emailing me at email@example.com if you want it) but almost any competent plotting program will do.
I have a small program snm.f that computes the mean PSD in the 200-600 s band and the SNM. For the above data the result should be
PART III - INTERPRETATION
ICET already received comments some years ago by publishing a 'Quality Factor' for the quality of data for a number of worldwide gravity stations. Nevertheless the idea of comparing the instrument+noise of different sites is obviously extremely useful and interesting to a number of people. The SNM (and the whole PSD which is actually much more useful) is only one method, but it is at least published.
I have recently re-analyzed all the raw data I have for Boulder and Cantley with the following results, to be used in conjunction with the referenced article.
The following results were obtained by Heikki Virtanen In January 2000, based on the programs and instructions listed above on this Web page.
I suggest that people at other GGP sites perform similar analyses. If you wish, I can do the analysis here if you will make available the appropriate data. In any case, if you do your own analysis (according to the above procedure) please send me the raw gravity and fixed pressure data for your quietest 5 days and I will confirm your results.
Banka, D. and Crossley, D.J., 1999. Noise levels of superconducting gravimeters at seismic frequencies, Geophys. J. Int., 139, 87-97.