Dispersion analysis of spaced antenna scintillation measurement

We present a dispersion analysis of the phase of GPS signals received at high latitude. Basic theoretical aspects for spectral analysis of two-point measurement are given. To account for nonstationarity and statistical robustness a power distribution of the windowed Fourier transform cross-spectra as a function of frequency and phase is analysed using the Radon transform.


Introduction
The effects of interaction of electromagnetic waves with the plasma are used in various radio techniques to investigate the physical properties and dynamics of the ionosphere. In particular, satellite radio beacons are effectively applied to measure ionospheric total electron content (TEC). Beacon satellites are also used to study ionospheric irregularities. The amplitude, phase, and angle of arrival of the signals traversing the ionosphere are distorted by drifting small scale plasma density irregularities. Fluctuations of the wave parameters are called scintillations. Scintillations, when measured on the ground, provide information about the structure and motion of ionospheric irregularities (Yeh and Liu, 1982).
With increasing significance of satellite-based methods in precise positioning and timing, the scintillation studies become of growing importance. No wonder, scintillations can reduce the accuracy of positioning and even temporarily prevent the use of satellites in navigation and geodesy. Especially vulnerable, from this point of view, are equatorial and high-latitude regions. High-latitude ionosphere are often remarkably irregular and dynamic. For that reason, monitor-Correspondence to: M. Grzesiak (pajak@cbk.waw.pl) ing of the high-latitude scintillation is of outmost importance. With this in mind, we got involved in the MISTECS (Monitoring of Ionospheric Scintillation and Total Electron Content on Spitsbergen) project. Monitoring equipment has been installed in the Polish Polar Station in Hornsund on Spitsbergen.
One of the components of MISTECS is the spaced receivers experiment in which three GPS scintillation receivers (GISTM) are placed at the corners of the triangle (Fig. 1). The aim of the spaced receivers experiment is estimation of the diffraction pattern drift velocities and anisotropy of scintillation-producing electron density irregularities.
In general, two basic methods are used to derive drift velocity: (i) the cross-correlation analysis (Rino and Livingston, 1982;Kintner et al., 2004;Otsuka et al., 2006), and (ii) the dispersive (cross-spectral) analysis (Costa et al., 1988a;Briggs, 1968). The differences between two methods are discussed by Briggs and Golley (1968), and Wernik et al. (1983). Application of the dispersive method allows to study the frequency variations of the drift velocity and irregularity anisotropy. It was found that the velocity varies with frequency, while the velocity direction remains unchanged. The frequency variations of the velocity have been attributed to the velocity shear along the propagation path (Briggs and Golley, 1968;Wernik et al., 1983;. In this paper we propose a new method of analysis of spaced receivers data, which uses the windowed Fourier transform (WFT) to take care of possible nonstationarity of data, and Radon transform (RT) that can capture characteristic linear dependence of the cross-spectrum phase on frequency, provided the scintillation pattern on the ground is drifting. The advantage of the method is that it allows to obtain frequency variations of the drift velocity even for nonstationary data.

Data analyzed
Data analyzed in the paper come from the MISTECS project that was established during International Polar Year in 2007. The hardware part of the experiment consists of three GPS receivers placed around the Polish Polar Station (77 • N 15 • 33 E) in Hornsund, South Spitsbergen. The experiment aims in monitoring of ionospheric GPS signal scintillation making use of the amplitude and phase sampled with the frequency 50 Hz. Three spaced antennas (for setup and orientation see Fig. 1) gives also a possibility of determinig the ionospheric drifts based on two point measurement. More information on the MISTECS experiment can be found in Wernik et al. (2008). The data case under consideration was measured on 21 November 2007. It is a ten minutes long record of high resolution GPS phase (PRN 4) made arround noon magnetic time which indicates a conjunction to the magnetic cusp region under moderately active geomagnetic conditions (K p =4). Figure 2a shows the analyzed raw receiver GPS phase data before preprocessing. The substantial feature of the plot is a trend that mainly contributes to the measurements and is possibly due to the relative motion of the GPS satellite and receiver. This is confirmed by Fig. 2b -the data plotted are data from Fig. 2a   result is a sum of pieces of parabolas that is consistent with previous explanation for a global trend -the trajectory of the relative motion of GPS satellite and the receiver in known to be a second order to large extent. Figure 2c depicts the same data set after subtracting local trends of third order -the satellite contribution to phase measutement has been totaly removed. The removal of the trends was accomplished by removing coarse-scale coefficients of the Deslauriers-Dubuc biortogonal wavelet transform (Daubechies, 1992) of desired order. In Fig. 3 the data for three receivers are shown on the same plot -red for receiver numbered 05, green -06, blue -16. Remarkable features of time series are nonstationarity of individual receivers data and its mutual relation as well. This mutual relation between measurements taken at two distinct spatial position is a key property for deriving the velocity of a drifting field (Briggs, 1968).

Dispersive evolution
To have precise description of our experiment and justify data analysis method let us start with a description of a linear steady and homogenous evolution of a scalar field ψ which will be a model for the observed diffraction patern. We consider an autonomous system where an initial field ψ(r,0) evolve in time following the rule: Equation (1) should be thought of as the solution to an initial value problem of a partial differential equation and the K(r − r,t) kernel is closely related to the Green's function. Under some conditions (among others the translation invariance with respect to time and space coordinates, see Engel and Nagel, 2000) equation can be wtitten in the Fourier representation as: The function (k) describes physical properties of the system under consideration. Then averaging over an ensamble of initial conditions, and assuming that initial field is statistically homogenous with the power spectrum P (k) we get: is the expectation operator and C(ζ ,τ ) is the spatio-temporal autocorrelation function of the field. As we can see the translation invariance in space and time imply stationarity and homogenity of the random process. Now we can compute cross-spectrum for two spatially separated points available in our experiment, which is the basic quantity for spectral analysis: the ζ is vector that separates points of measurement. Unlike in Eq. (4) the phase of cross-spectrum will depend in general on (k), P (k) and the separation vector ζ . We illustrate this with a simple case of translation, where kernel K(r − r,t) takes the form: which is equivalent to the simple transport equation: and the formula for P (ζ ,ω) reads: Choosing the coordinate system such that v d ||k x we can write it as: dk x e ik x ζ x δ(k x v d − ω) dk y P (k)e ik y ζ y = 1 v d e i ωζx v d dk y P ω v d ,k y e ik y ζ y .
When the integral in Eq. (12) does not give contribution to the phase, what is the case for P (k) symmetric about v d (for example P (k) is isotropic), the phase of the cross-spectrum: contains direct information about drift velocity, where β is the angle between separation vector ζ and the drift velocity v d (see Fig. 4). If one takes two perpendicular separation vectors (three observation points) then writing v 1(2) = ζ 1(2) ω φ 1(2) = v d cosβ(sinβ) it can be shown that: The basic theory of two-point measurement presented applies, as we mentioned, to the diffraction pattern recorded by receivers on the ground. Its relation to the ionospheric electron concentration field is a subject of numerous papers (Wernik et al., 1983;Costa et al., 1988a). On the other hand there is a complementary method of deriving drift velocity of moving patterns -the cross-correlation method (Briggs and Golley, 1968;Costa et al., 1988a). Figure 5 represents phase of the Fourier cross-spectrum averaged over eight records as a function of frequency. One can see a region 0.01-1.5 Hz where the phase depends linearly on frequency but it looks very noisy and above 1.5 Hz it could be hardly seen if at all. We propose to use a technique described in Dudok de Wit et al. (1995) and Pincon et al. (1997) instead. The method takes advantage of time-frequency trans- forms (the wavelet or windowed Fourier transform for example) to identify nonstationary coherent features out of a stationary noisy part. For this purpose we take the windowed Fourier transform (WFT) (Daubechies, 1992) of received signal:

Dispersion analysis
which gives constant frequency resolution across a whole frequency range, then construct WFT cross-spectrum for signals f and g from two GPS receivers: W f g =f (τ,ω)g(τ,ω) = |f (τ,ω)g(τ,ω)|e i φ(τ,ω) . We build a sort of histograms of the power of W f g coefficients when its phases fall into a cartain binρ( φ x , ω) and call from now on the phase-frequency distribution. Such distribution normalized to its maximum for each frequency is shown in Fig. 6. The maxima organize along lines (the maxima lines) which show dominant, in terms of power, dependence of phase difference vs. frequency. In comparison with the FT the maxima lines of the WFT phase-frequency distribution follow linear dependence on frequency more notably and the frequency range of this dependence is broader. This result confirms the fact about the drift of the diffraction pattern observed on the ground and, indirectly, motion of ionospheric irregularities. Although the result above confirms that the diffraction pattern drifts, information we got is of qualitative nature. To get magnitude and direction of velocity vector one could compute for example an average where the weights w i (ω) denote total power of the WFT cross-spectrum coefficients whose phases are φ i (ω) at a given frequency ω. With this approach we have at least two problems: it is essentially the same except for lower frequency resolution as global FT method, and there is a risk to miss some features caused by nonstationarity -drift velocity can change resulting in more than one linear features over the same frequency range.
To overcome these difficulties we have used the Radon transform (RT) of the WFT phase-frequency distribution. The Radon transform of the two-dimenssional function f (x,y)is defined by Beylkin (1987): It is an integration along lines parametrized by their position (for example point of intersection with the line φ = 0) and direction (proportional in our case to the velocity). It is clear that we have just modified averaging over one frequency by introducing variable direction in the phase-frequency plane.
For simple linear structures averaging over maxima lines should give more significant result than suming allong neighbour (in position and direction) lines. Such an approach addresses also the problem of nonstationarity. As we pointed out changing drift velocity should give maxima lines with different direction over the same frequency range that can be captured by RT.
To compare our method to technique based on Fourier transform we applied both to a pair of one-dimensional test signals. First signal is just a coloured Gaussian noise with power-law power spectrum and second signal was obtained by translating first one by a time shift that changes at the middle of the time interval, which should change the drift velocity (we took time shifts so that resulting drift velocities change sign too). We added to the second signal a Gaussian white noise. Figure 7 shows relation φ(ω) using averaged Fourier transform -there is no clear evidence of drift in this analysis. In comparison Fig. 8 depicts the result obtained using our method. In upper panel phase-frequency distribution contains two families of parallel maxima lines with different slopes as one could expect for nonstationary drift. Highfrequency part of maxima lines diffuse in the noise which starts to dominate there since noise power spectrum is white. We also marked with arrows two maxima lines and shaw φ = 0 line. The intersection point of the maxima line with this line gives frequency value of maximum point in the RT domain shown in lower panel of Fig. 8. Figure 9 ilustrates the RT of phase-frequency distributions depicted in Fig. 6 for two pairs of GPS antennas -upper panel S-N pair and lower E-W, respectively. The isolated maxima correspond to maxima lines on the phase-frequency distribution but now we can determine a velocity that they define. Taking into account the relation between slope of maxima lines and velocity given by relations (14) and (15) it is possible to asses the velocity magnitude and direction, which is sumarized for arrow-marked maxima (Fig. 10) in Table 1. There is an interesting feature for both pairs of receivers pointed by white arrow in Fig. 6 and red one in Fig. 10-maxima line with opposite signs of slope within the same frequency range that suggest changing drift velocity i.e. nonstationary drift.

Conclusions
The topic of the paper is called "disspersion analysis" (Dudok de Wit et al., 1995) since it allows for different time response of a system subject to different space scale forcing. Our analysis shows that it is possible to state that small structures of diffraction pattern drift slower than large ones. Above 3 Hz it is hard to find one-to-one correspondence between RT for two pairs of receivers, which could be due to spatial aliasing (wavelenghts gets too small to be resolved over separation distance between receivers) or nonstationarity. While spatial aliasing can be reduced by decreasing the distance between receivers, nonstationarity requires to find forementioned correspondence i.e. which maximum on one picture correspond to which on the other one. Unfortunately information on time that could resolve this problem got lost during the analysis. Another solution could be to perform similar analysis on the joint phase-frequency distribution for two pairs of receivers ρ( φ x , φ y , ω), that requires (in our opinion) other than RT tools for analysis.