Comparison of Earth Tides Analysis Programs

O. Dierks and J. Neumeyer

GeoForschungsZentrum Potsdam, Division 1,

Telegrafenberg, 14473 Potsdam, Germany

dierks@gfz-potsdam.de

1 Introduction

The increasing accuracy of gravimeters especially of those with superconducting technology and the comparison between many stations and working groups in projects like the “Global Geodynamic Project (GGP)” lead to the problem that possible differences in the results occur from differences in used analysis programs. This article describes three of the widely used analysis programs for earth tide analysis and shows the results of some comparing investigations made at the GFZ Potsdam for a diploma thesis.

2 Preprocessing

Before the tidal analysis can be performed the raw gravity data have to be repaired and filtered. This is done by using a pre-processing software. Two of these programs were also tested (PRETERNA 3.30 and TSOFT 1.1.4).

The preprocessing can be divided into several parts:

- calibration of the data

- calculation of the gravity residuals

- removing of spikes, steps and interfering signals

- filling of data gaps by theoretical values

- filtering and decimation to the required sampling rate for the tidal analysis software.

Central step is the correction of the data. The automated and hand made correction of steps and spikes is different for both programs. The automatically correction of large or multiple steps is not possible with both programs. In program PREGRED from the ETERNA-package it is difficult to overlook the data set with the small zoom functions. Corrected and uncorrected signals are visible in program TSOFT. Each correction can be rejected. Many mathematical and stochastically calculations can be applied on the data. A lot of datasets can be managed and handled in the program simultaneously. The user can program an automation of the correction without loosing visible control.

Within the other steps of the preprocessing the two programs behave equal. ETERNA-package is in some parts even more powerful by using the newest tidal potential catalogue and tuned digital numerical filters. The program TSOFT computes the filter coefficients each time they are needed according to the settings by the user. But this great choice possibly bears errors.

3 Description of the Analysis Programs

The following three widely used analysis programs were tested in this work:

- Program ANALYZE from the ETERNA-package, version 3.30, by H. - G. Wenzel [WEN-96a], [WEN-96b]

- Program BAYTAP-G in the version from 15.11.1999 by Y. Tamura [TAM-91]

- Program VAV in the version from April 2002 by A. P. Venedikov et. al [VEN-01]

The changing gravitational forces from sun, moon and the planets affect the earth. In the earth’s centre of mass this gravitational force is compensated by the centrifugal forces due to the motion of the earth around the sun and due to the motion of the earth around the barycentre of the moon-earth system respectively. Centrifugal acceleration is constant in every point of the earth but gravitational acceleration is different due to spatial extent of the earth. The small resulting acceleration is called tidal acceleration. The tidal acceleration in a fixed point of the earth changes with time because of earth rotation and movement of the participating bodies.

The changes of the tidal acceleration in a fixed point can be recorded by the use of e. g. a gravimeter. But the recorded data series y(t) does not only consists of the observed tidal gravity signal w(t) but also contains further information:

(1)

y(t) = w(t) + d(t) + a · a(t) + e(t).

Term d(t) describes the drift of the gravimeter. Term a(t) is a time series with meteorological or hydrological data. Coefficient a describes the influence of this additional parameter onto the gravity measurement. Further signals and measurement errors are combined in term e(t).

An analysis program corrects the observed signal y(t) by eliminating the drift series and the influence of the meteorological and hydrological signals. Using the coordinates of the station and a tidal potential catalogue a theoretical tidal gravity signal is computed. A comparison between this theoretical and the observed tidal gravity signal is used to estimate a set of tidal parameters (amplitude factor and phase lead) for the station. The tidal parameters amplitude factor and phase lead cannot be determined for each wave noted in the tidal potential catalogue. Following the Rayleigh-criterion the waves of the used tidal potential catalogue are stacked together to wave groups [VEN-61]. For each of these groups the tidal parameters amplitude factor and phase lead are estimated.

After some general comments on each program its methods are described to compute tidal parameters, influence of additional signals and accuracy of the results.

3.1 ANALYZE

The program is based on a method developed by Chojnicki [CHO-73] and improved by Schüller [SCH-76] and Wenzel [WEN-96a]. A least square adjustment is used to estimate the tidal parameters, the meteorological and hydrological regression parameters, the pole tide regression parameters and the TSCHEBYSCHEFF polynomial bias parameters for drift determination. The amount of data is nearly unlimited. Every kind of earth-tide data (gravity, strain, tilt and displacement) and up to eight channels with meteorological and hydrological data can be analysed. The user can determine the range of up to 85 wave groups. One tidal potential catalogue out of seven including the newest from HARTMANN and WENZEL [HAR-95] can be chosen to calculate the theoretical tidal signal. On the other hand the requirements on the format of the data and the parameter file are very stringent [WEN-96b].

The model used for least square adjustment is:

				(2)
where	ℓ(t)	=	Observed gravity signal
	v(t)	=	Improvements to the observations
	X_j, Y_j	=	Linear form of unknown tidal parameters H_j (amplitude factor) and DF_j (phase difference) for each wave group j: X_j=H_j·cosDF_j Y_j=H_j·sinDF_j
	CO_j, SI_j	=	Factor of theoretical tidal parameters A_j (amplitude) and F_j (phase) for each wave i in the wave group j, starting with wave a_i and ending with wave e_i: = Amplification factor from digital highpass filter (equal 1 if drift is approximated by polynomials)
	D_k, T_k	=	Coefficients (D_k) of TSCHEBYSCHEFF-polynomials T_k of degree k
	R_m, z_m	=	Regression coefficients (R_m) of additional channel number m (z_m)

A possible drift in the data is eliminated by highpass filtering or is approximated by TSCHEBYSCHEFF–polynomials (T_k) whose coefficients (D_k) are also estimated in the least square adjustment. The filter coefficients for different numerical digital filters are included in the ETERNA-package. But the method of high pass filtering can only be used when no long periodic waves shall be determined. Together with the analysis of long periodic waves the drift has to be eliminated by an approximation through TSCHEBYSCHEFF-polynomials.

The influence of the air pressure data (or other meteorological or hydrological signals z_m(t)) onto the gravity measurement is determined by a linear regression. In the case of highpass filtering the air pressure data are filtered too and the regression is computed with the filtered data.

The accuracy of each parameter is determined in the least square adjustment in the form of standard deviations. The standard deviations of the tidal parameters are too optimistic and therefore corrected. They are multiplied by a factor that is derived from the spectrum of the residuals [WEN-96b].

3.2 BAYTAP-G

This program is based on a method called Bayesian prediction, developed in 1976 by HARRISON and STEVENS [HAR-76]. The method has been modified in Japan since 1983 for the use with earth-tide data. All kinds of earth-tide data can be analysed, but only three additional channels with meteorological or hydrological data are possible. The requirements are very stringent to the format of data and parameter files. The arrangement of the tidal wave groups is done automatically depending on the length of the time series, but the user can change the wave group boundaries by editing the corresponding file [TAM-90], [TAM-91]. The tidal potential catalogues from TAMURA [TAM-87] or CARTWRIGHT-TAYLOR-EDDEN [CAR-73] can be used.

Tidal parameters, drift and meteorological parameters are estimated through an iterative method similar to least square adjustment by minimizing the term [TAM-90]

(3)

A_m and B_m are the linear expressions of the unknowns amplitude factor and phase lead for each m of the M wave groups at all. C_mj and S_mj are computed from the tidal potential catalogue using all j waves contained the m^th wave group. This tidal part is subtracted from each observation y_i (n datapoints in total) together with the drift-value d_i and the term describing the influence of additional channels x(t) onto the measurement (see equation (5)). D and Weight are called hyperparameters and can be defined in the parameter file.

The second line of equation (3) is used for drift computation. Within this program the drift is not approximated by low degree polynomials. Here the drift is computed separately in each datapoint. The drift behaviour is characterized by the formula:

(4)

d_i = 2d_i-1 - d_i-2 + u_i

Here u_i is the stochastic part denoting a white noise sequence. d_i is the drift value at the current datapoint; d_i-1 and d_i-2 are the drift values in the two previous datapoints. The hyperparameter D can be used to fit the drift model to the data. A large value for parameter D causes an almost linear drift model; a small value leads to a drift model bending close to the data.

A similar possibility is given with hyperparameter Weight in the third line of equation (3). Here the variability of the tidal parameters can be chosen. But this option is only useful if too many tidal parameters shall be estimated from too poor data.

The influence onto gravity measurement is computed by regression for maximum three additional signals. But exceeding a simple linear regression the influence of more datapoints than the actual datapoint can be used:

(5)

Here parameter k gives the number of computation points for regression and Dt the time lag between the computation points if k > 0.

Within the iterative search the hyperparameter D is adjusted to get the best combination between parameters, measured data and tidal parameters. At the end of each turn an ABIC-value (ABIC = Akaike Bayesian Information Criterion) is computed. The solution with the smallest ABIC-value is the final one where data, parameter and drift fit each other best. This ABIC-value is also the only useful accuracy statement. A standard deviation is computed but following the author of this program it is simply derived from the ABIC-value. So this standard deviation is not comparable to standard deviations from the other programs.

3.3 VAV

Program VAV is based on a method called MV66 [VEN-66a], [VEN-66b] and an improvement of program NSV [VEN-97]. The data file can be adopted from program ANALYZE, but program VAV has its own format for data files and uses own input files for tidal wave grouping and parameter settings. The wave group arrangement is done automatically depending on the length of the data set. Also a grouping variant can be chosen from a special input file. The used tidal potential catalogue is from TAMURA [TAM-87].

The fundamental idea of the program NSV [VEN-97] is a filtration of the original data containing an elimination of the drift and the separation into several pairs of series (step 1). Each pair contains signals from one main tidal species (D, SD or TD). The unknown tidal parameters for each tidal species are determined simultaneously but separately (step 2). This leads also to a frequency dependent accuracy statement (AKAIKE Information Criterion (AIC-value) and standard deviation). Both steps are also contained in program VAV but the separation in step 1 is not restricted to main tidal species. Here a wide spectrum of frequencies can be chosen by the user. Step 2 is using all the separated tidal species in a single least squares adjustment. An improvement of program VAV is the possibility to use data with different sampling rates and with several gaps in the same run without the need for interpolation [VEN‑01].

The original data set Y is divided into N intervals y_i of equal length. Each set contains n data points (n · N = M = total number of data points). n differs between the intervals, if the data are unequally spaced. Tidal parameters and air pressure regression coefficient are determined in a least squares fit together with the drift polynomial coefficients by minimisation of the following expression:

AX + PZ + E = Y

(6)

where AX is the tidal model including terms for the air pressure correction. Vector X is the vector of unknowns. Matrix E is the noise of the measurement. The model of the drift PZ is explained next.

In each interval y_i the drift is approximated by a polynomial of low degree (0 ≤ k ≤ 3). The matrix-notation of this looks like

(7)

with P containing the known polynomials depending on the time in each interval and Z containing the unknown polynomial coefficients representing the drift.

For each interval m matrixes c are created. Each matrix represents one of the chosen frequencies for filtering and separation

(8)

With j running from 1 to m denoting the angular frequencies used for separation and i running from 1 to N denoting the intervals. t₁, t₂, … t_N is the time series in each interval. The matrixes c_ij are transformed to f_ij and than merged together for each frequency to

(9)

All m F-matrixes are merged together with matrix P to a matrix called D so that

(10)

size (D) = (M,M) and D^TD = DD^T = I

The resultant identity-matrix depends already on the transformation of the c‑matrixes.

Throughout filtering and separation we get

(11)

The least squares fit (6) can than be changed into one using the filtered values

(12)

GX + E’ = U

without changing the results as shown in [VEN-01], but with estimation of more realistic accuracy statements.

Beside this standard deviation from the least square adjustment an AIC-value is computed. This value is used to compare different solutions of the same dataset with different parameter settings. The solution with the smallest AIC-value is the best one.

4 Comparative analysis with synthetic data

For comparing of the three tidal analysis programs a theoretical benchmark series (limited to degree 3) has been calculated and kindly provided by Bernard Ducarme (Royal Observatory of Belgium). The ten-year dataset of synthetic tidal acceleration (1 hour sampling rate) was computed for the SG-station BE (Brussels, Belgium). The used series with disturbed data is also a benchmark series with added red noise from Bernard Ducarme. The real data analysed in the next chapter was measured at SG-station SU (Sutherland, South Africa). The benchmark series are analysed with the three programs ANALYZE, BAYTAP-G and VAV. The obtained results were compared to the theoretical tidal parameters included in the benchmark series (amplitude factor = 1.0, phase lag = 0.0).

4.1 Analysis of pure synthetic data

The first tidal analysis of each program was started with the standard parameter settings (default values) offered by the programs (see end of chapter 4.2 for finally used values). All analyses have been carried out with the Tamura catalogue [TAM-87].

The following figures show the differences of the resolved amplitude factors (DAF) against 1.0 and the resolved phase leads (or their difference against 0.0; DPL). Figures 4.1, 4.2, 4.5 and 4.6 show the results for the three programs and the 31 used tidal wave groups. 31 wave groups is the maximum number of wave groups to be used with program BAYTAP-G. The spaces between the wave groups on the horizontal scale are due to figures 4.3, 4.4, 4.7 and 4.8. Here additional a fine wave group splitting with 54 tidal wave groups is used.

Because of the absence of any perturbation in the input data (synthetic tidal acceleration) the programs should calculate the same tidal parameters as included in the input data. The difference between the amplitude factors (DAF) or phase lead (DPL), respectively are shown in figure 4.1 and figure 4.2 for the 31 selected wave groups.

a) ANALYZE

The difference of the amplitude factors DAF is smaller than 0.00015 for all the selected wave groups (maximum: 0.00014 for wave groups S1 and PSI1). For the phase lead the difference DPL is smaller than 0.016° (2N2). The application of the HANN-window or a change of the numerical highpass filter lead to nearly no changes.

b) BAYTAP-G

The deviations of the amplitude factors DAF are similar small (maximum: 0.00013). The biggest DAF are concentrated on the lower frequencies (wave groups SGMQ1, 2Q1, SIG1). The phase lead has a maximum deviation at M3 (0.051°). It is astonishing that all DAF are negative (amplitude factor smaller 1.0) and all DPL positive (phase lead greater 0.0).

Figure 4.1: Deviations of amplitude factors for analysis of pure benchmark data

Figure 4.2: Deviations of phase leads for analysis of pure benchmark data

Figure 4.3a: Deviations of amplitude factors after analysis of pure benchmark data with 31 wave group

Figure 4.3b: Deviations of amplitude factors after analysis of pure benchmark data with 54 wave groups

Figure 4.4a: Deviations of phase leads after analysis of pure benchmark data with 31 wave group

Figure 4.4b: Deviations of phase leads after analysis of pure benchmark data with 54 wave groups

c) VAV

The deviations of the amplitude factors DAF are smaller than 0.00010 (maximum: 0.00006 for wave group S1). The DPL are also comparable to the results of program ANAYLZE (maximum: 0.012° for 2N2 and 2Q1).

A further comparison was made by the use of a finer wave group splitting with 54 wave groups instead of 31. This grouping was used in a second run for the programs ANALYZE and VAV. Program BAYTAP-G is not able to use more than 31 wave groups. The comparison between programs ANALYZE and VAV is shown in the figures 4.3 (DAF) and 4.4 (DPL). In the upper part of both figures the results with 31 wave groups are shown. The wave groups marked with a cross are not analysed within this grouping. In the lower part of both figures the results of the finer grouping with 54 wave groups are shown.

The graphs are printed with the same scales to see directly any differences between the two grouping variants. As expected the differences between the two grouping variants are very small. The deviations DAF and DPL are of course high for the small, new created wave groups (e. g. 3MK1, ALF2, BET2, DEL2, 3MQ2). The results of some wave groups become a little bit better for both programs by using the finer wave grouping. For a few other wave groups they become worse.

4.2 Analysis of disturbed synthetic data

In a second test theoretical red noise was added to the synthetic tidal acceleration data. The analysis of this disturbed benchmark series leads to worse results for the three programs.

The analysis with programs ANALYZE and BAYTAP-G lead to similar results of DAF and in some parts also of DPL (figures 4.5 and 4.6). The maximal deviations are for DAF 0.009 (PSI1) and for DPL 0.6° (SIGMQ1).

The results from program VAV are different for many wave groups to both other programs. Maximum deviation for DAF is also on wave group PSI1 but with a value of 0.012. For DPL the maximum deviation is also on PSI1 (0.56°). But for many wave groups the values and also the sign of DAF and DPL are different to ANALYZE and BAYTAP-G.

In figures 4.7 (DAF) and 4.8 (DPL) the programs VAV and ANALYZE are compared again. The upper part of both figures shows again the results of the analysis with 31 wave groups, the lower parts the analysis with 54 wave groups. The results (DAF and DPL) of programs ANALYZE and VAV do nearly not change between both grouping variants except for the new created tidal wave groups and their direct neighbouring wave groups (figures 4.7 and 4.8). The small changes between both wave groupings are as uneven as in figures 4.3 and 4.4. Here DAF and DPL are very high for the new created wave groups, especially in the semi diurnal frequency band.

The following parameter settings were used within the three programs to analyse the benchmark series. The analysis of benchmark data with program ANALYZE leads to best results when using RIGIDEARTH=1, as was recommended by the author of the program. The HANN-window was not used (HANNWINDOW=0).

Figure 4.5: Deviations of amplitude factors for analysis of noisy benchmark

Figure 4.6: Deviations of phase leads for analysis of noisy benchmark

Figure 4.7a: Deviations of amplitude factors for analysis of noisy benchmark data with 31 wave groups

Figure 4.7b: Deviations of amplitude factors for analysis of noisy benchmark data with 54 wave groups

Figure 4.8a: Deviations of phase leads for analysis of noisy benchmark data with 31 wave groups

Figure 4.8b: Deviations of phase leads for analysis of noisy benchmark data with 54 wave groups

Parameter NUMFILNAME does not really have influence on the analysis (NUMFILNAME=n1h1h007.nlf and 002 were finally used). The results (DAF and DPL) get better for some wave groups for others they get worse. WAVEGROUPI and TIDALPOTEN were fixed to use the same wave grouping and tidal potential catalogue in all three programs.

Program BAYTAP-G leads to best results with LOVENM=0 (similar to RIGIDEARTH). Parameter DMIN, managing the search for the smallest ABIC-value, had to be set to something equal 0.1. Then the search for the smallest ABIC-value comes to an end. But the parameter DMIN has not much influence on the results as long as the solution with the smallest ABIC-value is detected. Parameter P4FLAG which includes the waves of potential of degree 4 is best set equal 1 (4^th order potential is contained) and parameter ORDER = 2 managing the variability of the drift (ORDER = 2 in the midterm of equation (3)). Parameters ITH=1 and IPOTEN=2 are fixed to use the tidal potential catalogue of TAMURA. SPAN and SHIFT were set to 0 to analyse the whole data set in just a single run.

The variety of parameters within program VAV is very big (three times as much as for programs ANALYZE and BAYTAP-G). But for the analysis of this benchmark series a lot of them are not used anyway. The time window (length of the filter intervals) had a length of 48 hours and the power of the drift polynomials was set to 3. In this drift time series the long periodic tidal waves are included. The long periodic signals cannot be eliminated by filters. The only other possibility is to analyse this long periodic waves too. The filter-frequencies were set to 15°/h, 30°/h, 45°/h. The wave grouping was the same as in the two other programs. The possibilities of separation or correction of the waves of third or fourth order potential were not used.

5 Comparative analysis of observed data

The monthly corrected and decimated datasets were connected to a long dataset with duration of ten month. After preparation of data in BAYTAP-format the wave group arrangement was adjusted so that all three programs use the same arrangement of 31 tidal wave groups together with the same tidal potential catalogue of TAMURA.

Program ANALYZE computes a regression coefficient between highpass filtered gravity data and highpass filtered air pressure data. Here no parameters can be set, the program just has to know which additional signal is supplied. Highpass filtering is applied to eliminate the drift. An approximation with TSCHEBYSCHEFF-polynomials is not used. This method only makes sense when analysing long period wave groups. Parameter RIGIDEARTH now has been set to 0 for real data and NUMFILNAME=n1h1h002.nlf and HANNWINDOW=1 led to better results.

The regression coefficient is estimated in program BAYTAP-G with the model shown in equation (5). In the analysis with the lowest ABIC-value the parameters were set to k_max=Dt = 0 = LAGP, IAUG = 1, LAGINT = 0 (arbitrary when LAGP = 0). As in the programs ANALYZE and VAV just a linear regression with only the actual air pressure value is computed than. Parameter LOVENM is set to 2, parameter DMIN is set as low as 0.01.

Figure 5.1: Deviations of amplitude factors for analysis of real data

Figure 5.2: Deviations of phase leads for analysis of real data

Figure 5.3: Amplitude spectrum of residuals computed by program ANALYZE

Figure 5.4: Amplitude spectrum of residuals computed by program VAV

Figure 5.5: Amplitude spectrum of residuals computed by program BAYTAP-G

The spectra from the residuals calculated by ANALYZE and VAV are very similar. While the decreasing amplitude below frequencies of 0.8 cpd results from highpass filtering there is no explanation for the very small amplitudes (different scale!) and the vanishing amplitudes already below frequencies of 1.5 cpd within the spectra of the residuals from program BAYTAP-G. Both other spectra show for single frequencies more energy than for most others. These are remaining signals that could not be assigned to one of the parts named in equation (1) (e. g. air pressure, tidal signal).

6 Conclusions

Combining the direct results of the different analyses for each program does not lead to an advantage of one of these three programs. The results (DAF and DPL) of one program are not better than for another throughout the different wave groups, grouping variants and analyses. But it has to be said that all these results were obtained with special data sets, a fixed wave grouping and just one tidal potential catalogue. By using further applications (e. g. catalogue of HARTMANN and WENZEL in analyses with ANALYZE) different and possibly better results may be obtained.

Program BAYTAP-G is in the used version not able to analyse the long periodic waves. But this could be very interesting when analysing longer data sets of a superconducting gravimeter with its high stability especially in this frequency band. The program BAYTAP-L is made for this frequency range but then there may be the problem that the two frequency parts of one data set (LP and D to TD) are analysed by two separate programs.

Program ANALYZE shows up with good results during the tests (except for the very small wave groups ALF2, BET2, DEL2), a good and complete documentation and a broad output without using special parameter settings. Only this program offers the possibility to use the newest and most accurately tidal potential catalogue of HARTMANN and WENZEL [HAR-95].

The numeric results from program VAV do nearly not differ to those from program ANALYZE. And when differences occur they are small and have changing sign. A problem with program VAV is the slight documentation. There are offered more than 70 parameters to manage the program but only just very little documentation about the effects of each parameter and how to handle them. Beside the advantage of many possibilities to control the analysis there is the danger of not knowing the influence onto the computation.

Acknowledgement

Many thanks to Bernard Ducarme (Royal Observatory of Belgium) for the computation of the benchmark series and the analyses, especially with ANALYZE and to Angel Venedikov (Geophysical Institute at the Bulgarian Academy of Science) for the analyses with VAV and the support for this program.

7 References

[CAR-73] Cartwright, D. E.; Edden, A. C.: Corrected tables of tidal harmonics. In: The Geophysical Journal. (1973), Vol. 3, No. 33

[CHO-73] Chojnicki, T.: Ein Verfahren zur Erdgezeitenanalyse in Anlehnung an das Prinzip der kleinsten Quadrate. In: Mitteilungen aus dem Institut für Theoretische Geodäsie der Universität Bonn, No. 15. Bonn: 1973

[HAR-76] Harrison, P. J.; Stevens, C. F.: Bayesian forecasting (with discussion). In: Journal of the Royal Statistical Society B, No. 38, pp. 205-247, 1976

[HAR-95] Hartmann, T.; Wenzel, H.-G.: Catalogue HW95 of the Tide Generating Potential. In: Marees Terrestres Bulletin d’Informations, No. 123, pp. 9278-9301, Paul Melchior (Edit.). Bruxelles: Association Internationale de Géodésie, 1995

[SCH-76] Schüller, K.: Ein Beitrag zur Auswertung von Erdgezeitenregistrierungen. In: Deutsche Geodätische Kommission bei der Bayerischen Akademie der Wissenschaften, No. 227. München: 1976

[TAM-87] Tamura, Y.: A harmonic development of tide-generating potential. In: Marees Terrestres Bulletin d’Informations, No. 99, pp. 6813-6855, Paul Melchior (Edit.). Bruxelles: Association Internationale de Géodésie, 1987

[TAM-90] Tamura, Y.: BAYTAP-G Users Manual. Mizusawa: National Astronomical Observatory

[TAM-91] Tamura, Y.; Sato, T.; Ooe, M.; Ishiguro, M.: A procedure for tidal analysis with a Bayesian information criterion. In: Geophysical Journal International, No. 104, pp. 507-516. Oxford: 1991

[VEN-61] Venedikov, A. P.: Application à l’Analyse Harmonique des Observations des Marées Terrestres de la Méthode des Moindres Carrées. In: C. R. Acad. Bulg. Sci. (1961), No. 14, pp. 671-674

[VEN-66a] Venedikov, A. P.: Une méthode d’analyse des Marées terrestres à partir d’enregistrements des longueur arbitraires. In: Bull. Cl. Sci., 5^e S., t. LIII, fasc 3, Communications de l’Observatoire Royal de Belgique, Série géoph., No. 71, pp. 463-485. Acad. Royal de Belgique, 1966

[VEN-66b] Venedikov, A. P.: Sur la constitution de filters numériques pour le traitement des enregistrements des Marées terrestres. In: Bull. Cl. Sci., 5^e S., t. LIII, fasc 6, Communications de l’Observatoire Royal de Belgique, Série géoph., No. 76, pp. 827-845. Acad. Royal de Belgique, 1966

[VEN-97] Venedikov, A. P., Vieira, R.; de Toro, C.; Arnoso, J.: A new program developed in Madrid for tidal data processing. In: Marees Terrestres Bulletin d’Informations, No. 126, pp. 9669-9704, Paul Melchior (Edit.). Bruxelles: Association Internationale de Géodésie, 1997

[VEN-98] Venedikov, A. P.: Output when calling the help function in NSV98. Extracted from source code <NCONI.FOR> in the version from 01. Sept. 1998

[VEN-01] Venedikov, A. P.; Arnoso, J.; Vieira, R.: Program VAV/2000 for tidal analysis of unevenly spaced data with irregular drift and coloured noise. In: Journal of the Geodetic Society of Japan, Vol. 47, No. 1, pp. 281-286, 2001

[WEN-94] Wenzel, H.-G.: Erdgezeitenanalyse. In: Bericht zum DGG-Seminar Gezeiten in Oberwolfach / Schwarzwald Oktober 1994, Sonderband II / 1995, Deutsche Geophysikalische Gesellschaft e.V. Münster 1995, pp. 19 - 38

[WEN-96a] Wenzel, H.-G.: The nanogal software: Earth tide data processing package ETERNA 3.30. In: Marees Terrestres Bulletin d’Informations, No. 124, pp. 9425-9439, Paul Melchior (Edit.). Bruxelles: Association Internationale de Géodésie, 1996

[WEN-96b] Wenzel, H.-G.: ETERNA Manual. Version 3.30, Karlsruhe: Black Forest Observatory

		Coefficient in nms^-2/hPa	Standard deviation in nms^-2/hPa
ANALYZE		-2,82	0,02
BAYTAP-G		-3,16	0,02
VAV		-2,81	0,09
VAV	D	-2,94	0,11
	SD	-2,75	0,20
	TD	-2,44	0,11