## 1 Introduction

The Doppler Orbit Determination and Radiopositioning Integrated by Satellite (DORIS) system inaugurates a new means of scientific measurement based on the determination of orbit and radiopositioning [8,17]. It is a French civil system conceived and developed by CNES (French Space Agency), GRGS (French Research Group in Space Geodesy) and the IGN (French Mapping Agency) [3].

Since 1990, the beginning of the first DORIS mission, all of the measurements collected by the instruments embarked on SPOT2, SPOT 3, SPOT 4, TOPEX POSEIDON and JASON represent several tens of millions of precise observations placed at the disposal of the scientific community through the International DORIS Service [16].

The initial objectives of the Doris systems were Precise Orbit Determination of Low Earth Orbiting (LEO) satellites, gravity field determination and precise positioning [3] (expected to be at the 10 cm level before launch of SPOT2). In eight years, because of the restoration and installation of new stations and also by the improvement of models and strategies of computation, a significant improvement in the quality of network [7] was carried out, making it possible for the DORIS system to contribute significantly to the realization of the international reference terrestrial frame.

The richness of the measurements collected by the DORIS space system makes it possible today to represent the displacement of the stations on the ground in the form of station coordinate time series. Such series make it possible to detect, starting from analysis methods of observational, instrumental or geophysical parasitic phenomena origins, the local motion study of the stations; some methods are more sensitive to the different types of noises and their levels (autocorrelation function, spectral concentration of Allan variance, etc.) [12]; others are based on the search for periodicity, such as Fourier analysis, Singular Spectral Analysis (SA), Method of Maximum of entropy (MEM) [9], and also the wavelet method used in this article, which treats the non-stationary signals and makes it possible to extract information from the structures and the local irregularities.

The Fourier transform application provides only one piece of information on the total regularity of the signal and does not allow one to easily obtain information located in time. The following idea consists in representing our signal according to time and frequency. One of the first to have applied this principle to Fourier transforms is the physicist Dennis Gabor in 1940. One speaks then about Fourier transforms with windowing of the signal. The basic idea consists in cutting the signal in temporal sections and adapting, on each section, the Fourier transform. This analysis depends on the localisation of the section.

The main inconvenience of this technique is that the length of the bench (the scale) is fixed (Fig. 1) which prevents one from retaining the time-frequency representation and achieving a multiscale/multiresolution. In this case, we use the signal processing wavelet which takes into account the local regularity and makes a kind of zoom in/out on a portion of the signal by playing on the scale factor.

In the present article, we present an analysis methodology based on the wavelet technique using times series of DORIS station coordinates provided by the IGN/JPL center [18], and we apply the wavelet packet transform for denoising the signal. Our objective is to provide certain information useful to later geodynamic interpretations.

## 2 Discrete wavelet analysis

A wavelet is a function $\psi \left(x\right)$ of L_{2}(R) (L_{2}(R): space continuous functions square integrable having properties of localization and oscillation [6]. Just as the Fourier transform can be defined as being a projection on the basis of the complex exponential, we introduce the wavelet transform as projection on the base of wavelets functions:

$$\text{TO}(a,b)=\frac{1}{\sqrt{\text{a}}}{\int}_{-\infty}^{+\infty}x(t)\psi \left(\frac{\text{t}-\text{b}}{\text{a}}\right)\text{dt}$$ | (1) |

^{2}

Starting from a mother wavelet [15], we can create a Hilbert base of L_{2}(R) by translation and dilation of mother wavelet $\psi $:

$$\left\{{\psi}_{a,b}(t)=\frac{1}{\sqrt{a}}\psi \left(\frac{t-b}{a}\right),\text{\hspace{1em}}(a,b)\in {R}^{2}\right\}$$ |

Thus, we break up a signal x(t) in a basis of the Hilbert space L_{2}(R) defined by the family of functions ${\psi}_{a,b}(t)$ and the expression of the wavelet coefficients becomes:

$$\text{TO}(a,b)={\int}_{-\infty}^{\text{+}\infty}x(t){\psi}_{a,b}(t)\text{dt}$$ |

Here, a and b represent the scale factor and the translation factor, respectively.

The functions ${\psi}_{a,b}\left(t\right)$ are obtained starting from dilation and translation of the mother function $\psi \left(t\right)$.

A. Grossmann and J. Morlet showed that if $\psi (t)$ is real valued, all of these wavelets can be considered as an orthonormal base, which means that any signal of energy can be written as a linear combination of wavelets ${\psi}_{a,b}(t)$ [11].

The wavelet coefficients (1) depend on the real a and b. We can decide on an infinite number of values for these two parameters which vary continuously, i.e.: a∈R^{+}\{0}, b∈R; we speak therefore about a continuous wavelet transform (CWT), which makes it possible to give us a better description of the pace of the signal. When we have a discretisation of the resolution parameter a and position b, we speak thus about the discrete wavelet transform (a = 2^{j}, b = k2^{j}, j, k∈Z^{2}).

The wavelet transform can be also considered as a process of decomposition of the signal into approximations and details; we speak then about the multiresolution analysis.

The concept of the multiresolution is to use several times the wavelet transformation (discrete case) in order to break up the starting signal f into several signals (scale of signal) containing different information. We define thus the family of sub-spaces {Vj} of L_{2} (R) which verifies:

$${V}_{j}\subset {V}_{j+1}\subset {L}_{2}\text{\hspace{1em}}\forall j\in Z;\text{\hspace{1em}}\cap {V}_{j}=\left\{0\right\}\text{\hspace{1em}}\cup {V}_{j}={L}_{2}$$ |

The approximation with the resolution j is defined by following projection:

$${A}_{n}^{j}=<f/{\Phi}_{j,n}>$$ |

Φ_{j,n} constitutes a basis of V_{j} and verifies:

$${\Phi}_{j,n}(t)={(\sqrt{2})}^{-j}\Phi ({2}^{-j}t-n)\text{\hspace{1em}}\text{with}\text{\hspace{1em}}j,n\in Z$$ |

Φ (t) ∈ V_{0} called the mother function or the scale function [10] defined by:

$$\Phi (t)=\sum _{n}h(n){\Phi}_{-1,n}={2}^{1/2}\sum _{n}h(n)\Phi (2t-n)$$ |

_{-1,n}> and Φ

_{-1,n}= 2

^{1/2}Φ (2t − n).

The loss information during this projection corresponds to f projections on a detail space of Wj, which satisfies: V_{j} ⊕ W_{j} = V_{j+1,} the detail coefficient is defined:

$${d}_{nj}=<f/{Y}_{j,n}>\text{\hspace{1em}}\text{such}\text{\hspace{0.17em}}\text{as}:{\psi}_{j,n}(t)={(\sqrt{2})}^{-j}\psi ({2}^{-j}t-n)$$ |

All this can be also translated like a filtering of our signal by two filters, a high-pass filter and a low-pass filter such as:

$${A}_{n}^{j}=\sum _{k}\tilde{h}(2n-l){A}_{n}^{j-1}\text{\hspace{1em}}\text{and}\text{\hspace{1em}}{d}_{n}^{j}=\sum _{k}\tilde{g}(2n-l){A}_{n}^{j-1}$$ |

The original signal passes through two complementary filters (Fig. 2): a high-pass filter gives a signal of details and a low-pass filter, a signal from the approximations. In its formulation, the wavelet transform can be interpreted as an analysis band of constant overpressure filter. In such a bench, each filter (band pass) can result from a single gauge by a dilation or compression frequency.

## 3 Wavelet packet transform

The analysis of the signal by the wavelet transform makes it possible to extract the noise from the original signal and to obtain a filtered signal, which causes a loss of information on the high frequencies of the original signal. The wavelet packet analysis [4] offers a rich analysis of the signal, thanks to its best model of decomposition in the tree form. The wavelet packet allows one to break up the signal on many bases and to choose within the meaning of certain criteria which represent the signal as closely as possible. In fact, in other words, not only the low-pass filtered versions of the signal are broken up, but also the high-pass filtered versions. The high frequencies are also cut out in sub-bands and the decomposition tree becomes symmetrical (Fig. 3). The resulting coefficients from this decomposition are characterised by three parameters: the level of decomposition, frequency index and temporal index. Each wavelet packet is carrying triple information (f: frequency, s: scale and p: position). This brings different techniques to code the signal and offers a richer range of possibilities for analysing it. We use in our application the wavelet packet of Meyer [14].

## 4 Denoising signal methods by wavelets

Denoising is a problem of signal treatment that is well-known and largely studied.

Let y(t) = x(t) + n(t), be the signal noise observed. The problem is to obtain an estimate of the original signal x(t). According to the model employed for these random variables x and n, several principles of denoising have been proposed. Some methods are based on the calculation of the Hodler exhibitors [13], and others use the thresholding method of denoising proposed by Donoho and Johnstone [5,6]. This work treats the second criterion based on the thresholding method. In the general case two types of thresholdings are used:

- • a hard thresholding:
$${S}_{\text{hard}}=\left\{\begin{array}{ll}X\hfill & \text{if}\left|X\right|>T\hfill \\ 0\hfill & \text{if}\left|X\right|<T\hfill \end{array}\right\}$$ If the observation signal coefficient is lower than a certain threshold T, it is considered as being pure noise and it is replaced by zero; if not, it is kept such at its value.

- • a soft thresholding:
$${S}_{\text{soft}}=\left\{\begin{array}{ll}sign(X)(|X|-T)\hfill & \text{if}\left|X\right|>T\hfill \\ 0\hfill & \text{if}\left|X\right|<T\hfill \end{array}\right\}$$

If the observation signal coefficient is lower than a certain threshold T, it is regarded as being pure noise and is replaced by zero; if not it is set to a narrowed value of the threshold.

In the case of a Gaussian additive white vibration, Donoho and Johnstone also proposed the universal thresholding (1994) defined as follows [5,6]: $T=\sigma \sqrt{2\mathrm{log}(n)}$, n is the number of signal points and σ is the estimated noise level: σ = M/0.6745 and M represents the absolute median estimated on the first scale.

In the case of correlated noise, Johnstone and Silverman proposed the thresholding (level depend) defines as follows [10]:

$$T={\sigma}_{t}\sqrt{2\mathrm{log}({n}_{t})}$$ |

_{i}points and σ = M

_{i}/0.6745.

Within the framework of this article, we noted that the Gaussian white vibration is more appropriate for our signal. We thus studied only the white vibration case and we used the thresholding minimax suggested by Donoho and Johnstone which uses the risk minimum principle [5,6].

The noise reduction method employing the wavelet transform consists in applying the signal decomposition concept in a wavelet adaptive base. This method is summarised by the four following stages:

- • the signal decomposition by choosing a wavelet and a resolution level N;
- • the choice of the best tree for this decomposition;
- • the selection of a threshold to the wavelet coefficients;
- • the signal reconstruction by wavelet starting from the coefficients of the resolution level N and the details coefficients of the level 1 to N.

## 5 Time series analyzes by wavelets

When we observe a time series of station coordinates, we want to detect the station displacement, therefore we search a slope, a tendency (change in the short or the long term of the series), periodic terms or irregular movements. This article proposes an analysis by wavelet transform of time series of DORIS station coordinates.

DORIS data used during the treatment make up of a set of weekly residuals time series coordinates East (dE), North (dN) and Vertical (dH) from the eight DORIS stations as we can see in Table 1.

**Table 1**

Representation of the eight DORIS stations selected.

Tableau 1 Représentation des huit stations DORIS choisies.

Stations | Geographical situation | Latitude | Longitude | Observation period (year) |

DIOA | Greece | 38°05′ | 23°56 | 1216 |

COLA | Sri Lanka | 6°53′ | 79°52′ | 11.76 |

FAIB | Alaska | 64°58′ | −147°31′ | 5.42 |

KRAB | Krasnoyarsk | 56° | 92°48′ | 7.47 |

SAKA | Yuzhno-Sakha linsk | 47°1′ | 142°43′ | 11.5 |

SODB | Scomo | 18°43′ | −110°57′ | 6.82 |

THUB | Greenland | 76°32′ | −68°49′ | 2.29 |

SYPB | Japan | −69° | 39°35′ | 5.94 |

They were selected on the following criterion basis [20]:

- • stations presenting discontinuities: KRAB and SAKA;
- • stations affected by the earthquakes: COLA, FAIB, SAKA;
- • stations affected by volcanic movements: SODB;
- • stations presenting data gaps and an antenna instability: DIOA;
**•**station representing better quality of geodetic data with high latitudes: SYPB and THUB.

### 5.1 Temporal representation of the Doppler Orbit Determination and Radiopositioning Integrated by Satellite (DORIS) time series

This stage consists in representing on a graph the residual coordinates (dE, dN and dH) of the DORIS stations according to time in order to carry out an examination a priori on the series distributions for the observation period. Fig. 4 represents the series distributions of the eight DORIS stations:

### 5.2 Discrete wavelet transform – study of tendencies and periods

#### 5.2.1 Tendencies

Knowing that the tendency represents the change in the short or the long term signal, wavelet analysis makes it possible to extract the tendency and the change in tendency of the series starting from the last resolution of the approximations (the 8th approximation) by using the discrete wavelet treatment of Meyer [14] which is very appropriate to our signal (Fig. 5).

The examination of the signal series of the eight DORIS stations (Fig. 5) enabled us to detect:

- • for DIOA station (Fig. 6): a descendent tendency from the beginning of the signal to 1997 and a data gap between 1997 and 1998 (null tendency). Another change of tendency (ascendent tendency) starting from 2001, and that for the three components dE, dH and dN;
- • for COLA station (Fig. 7): an overall descendent tendency for dE component. The descendent tendency becomes ascendent from the end of 2001 for both components dH and dN;
- • for FAIB station (Fig. 8): an ascendent tendency on all the signals for the two components dE and dH and a descendent tendency for the dN component;
- • for KRAB station (Fig. 9): an overall descendent tendency on all the signal for the dE components and ascendent tendency for the two other components, dH and dN;
- • for SAKA station (Fig. 10): an overall ascendent tendency on all the signal for the dE component, an overall descendent tendency for the two other components dH and dN. From the beginning of 2002 the tendency becomes ascendent for the three components;
- • for SODB station (Fig. 11): an overall descendent tendency on all the signal for the dE and dH components, and overall ascendent tendency for the dN component;
- • for THUB station (Fig. 12): an ascendent tendency on all the signal for the dE component and a descendent tendency for the two other components dH and dN;
- • for SYPB station (Fig. 13): an ascendent tendency on all the signal for the two components dE and dN and an overall descendent tendency for the dH component.

Finally, we notice that tendencies are almost identical for:

- • three components dE, dH and dN for DIOA station;
- • two components dH and dN for COLA, KRAB, SAKA and THUB stations;
- • two components dE and dH for FAIB and SODB stations;
- • two components dE and dN for SYPB station.

The signal evolution of the vertical component (dH) is similar to the North component (dN) which is due probably to the fact that the SPOTS satellites are heliosynchronous, the traces are north-south to the equator and consequently the East component (dE) of stations is not as well determined as the other components.

#### 5.2.2 Periods

The approximations are used to detect a tendency, a slope and irregularities of a signal, the details are employed to identify cycles (periodic signals). Fig. 6–13 represent the detailed signals starting from the 4th resolution of the three components for the eight stations.

According to these figures, we observe that the first details represent very important frequential extensions, i.e. small periods, but which are vague to distinguish. More the level of resolution increases, less the components become rich in frequencies (component periodicals are important). We can thus identify various periodic components on the resolutions and various levels of the details for the three components (dE, dH and dN) of the eight stations. Their periods increase with the resolution level. By making a zoom in the graph of the 4th detail of the vertical component of station SYPB, we notice a four-months period (= 118 days) (Fig. 14) which can result in a possible period error on the Z component of the geocentre [12,19].

We also note, on the vertical component level of DIOA, FAIB, KRAB and SYPB stations, a periodic signal of one year which is detached clearly on the level from the 5th detail D5. Concerning THUB station, we identified a six-months period on the level of the 4th detail of the East component, but we could not obtain the one-year periodic component owing to the fact that the observation period of this station is very short. For the other stations COLA, SAKA which are affected by earthquakes [19], the periodic components appear with difficulty from the graphs. Table 2 summarizes the various cycles detected on the level of the three components of the eight DORIS stations.

**Table 2**

Representation of the periods starting from the detailed coefficients of the stations.

Tableau 2 Représentation des périodes à partir des coefficients des détails des stations.

Components | Details | |||

D4 | D5 | D6 | D7 | |

dE, dH and dN | Period of 4 to 6 months | Period of 6 to 12 months | Period of 12 months with more than 24 months | Period of 24 months with more than 36 months |

The examination of these characteristics of the time series, revealed by the wavelets, will be able to help identify the natural geophysical phenomena related to the systematical error in the spatial geodesic computations (earthquakes, imperfections of the ITRF2000 model [2] which is replaced by the ITRF2005) [1].

### 5.3 Continuous wavelet transform (CWT): characterization of the series in time frequency

The CWT is the sum during the time of the signal multiplied by the function-shifted wavelet along the signal. This process produces the wavelet coefficients which depend on the scale and time; it allows us to characterize the time series in time and frequency, to give us a qualitative idea of the time series contents and to emphasize unexpected components (cycles, transitory problems), and to carry out a denoising of the series.

The wavelet coefficients obtained during our analysis are represented in coloured mode, in 3D (Fig. 15). The marked color represents the importance of the amplitude coefficients: dark shows small coefficients and light expresses high values.

Fig. 15 represents the periodic components along the signal; more the scale increases, more the signal has low frequencies and high values of coefficients. The periodic terms are variable in frequency or in amplitude.

A second analysis was applied (Fig. 16) to the denoised signal. We obtain new wavelet coefficients with better periodic components. We notice thus, for the eight vertical components of the eight stations, we detect the one-year periods for a low scale and six-months periods for large scale.

## 6 Wavelet packet transform

An interface of denoising signal using the wavelet packet transform with the soft thresholding applied to the vertical components (dH) of the eight DORIS stations was developed under MATLAB. The results obtained enabled us to notice that the analysed signal (filtered) of the eight stations (Fig. 17) became smoother, except that there remains some fluctuations on the level of the following stations:

- • DIOA: at the beginning of the signal;
- • COLA: towards the end of 1994 and also at the end of 2001 and at the beginning of 2002;
- • FAIB: in the middle of the year 2000;
- • KRAB: towards the beginning of 1998 and also the end of 1999;
- • SAKA: towards the beginning of the year 1993et 1995 and also beginning and at the end of 1997;
- • SODB: in the middle of the year 1998 and 1999 and also towards the end of 2001 and the beginning of 2002;
- • SYPB: towards the beginning of the year 2000 and the end of 2002.

On the other hand, at the THUB station, we note that the signal is smooth, which is certainly due to the fact that the observation period is relatively important.

All these fluctuations can be caused by a geophysical phenomenon (an earthquake), or by an antenna change or bad knowledge of one of the parameters being used with the estimate of the residues as coordinates. We also notice that all these fluctuations represent the high frequencies, which translates the analysis effectiveness of the wavelet packet transform.

## 7 Conclusion

The results obtained during the analysis of the residual time series of DORIS station coordinates showed that the wavelet method makes it possible to better determine and to quantify the systematic signals such as periodicity and tendency. Indeed, the tests carried out highlighted, on the one hand, the tendencies of the components (dE), (dN) and (dH) of each DORIS station. We notice almost identical tendencies of the vertical and the North component which is due probably to the fact that the East component less well determined as the other components because the SPOTS satellites are heliosynchrones. In addition, the annual and semi-annual periodic terms appear clearly on the detailed level of the wavelet coefficient for the eight station coordinates. Moreover, the 118 days period which arises from the graph of the 4th detail (D4) of the SYPB station provides a possible explanation error on the Z component of the geocentre. The extraction of the original signal noise by the wavelet transform can cause an information loss at high frequencies. The wavelet packet transform offers a richer analysis of the signal, thanks to its best model of decomposition in tree form and makes it possible to store more information on the high frequencies. Indeed, the wavelet packet transform application to the time series of DORIS station coordinates has allowed us, on the one hand, to obtain a smoother filtered signal, and on the other hand, to detect fluctuations representing the great frequencies which can be the cause of parasitic phenomena of observational, instrumental or geophysical origin.

The analysis of the geodetic time series is a field of research which is currently in a phase of development of the statistical studies adapted to the criteria of this new geodesy with four dimensions, while being often based on methods already developed in other disciplines. The work presented in this article is a contribution to these methodological developments based on the technique of the wavelets.

## Acknowledgements

The authors are very grateful to the referees for their suggestions, comments, constructive criticisms and hints which lead to substantial improvement of the paper. They also would like to thank the University of Oran Es-senia for supporting this research topics under CNEPRUS B01820060148 and B3101/01/05.