High-precision measurement of satellite range and velocity using the EISCAT radar

This paper is a continuation of an earlier work by Nygrén et al.(2012), where the velocity of a hard target was determined from a set of echo pulses reflected by the target flying through the radar beam. Here the method is extended to include the determination of range at a high accuracy. The method is as follows. First, the flight time of the pulse from the transmitter to the target is determined at an accuracy essentially better than the accuracy given by the sampling interval. This method makes use of the fact that the receiver filtering creates slopes at the phase flips of the phase modulated echo pulse. A precise flight time is found by investigating the echo amplitude within this slope. A value of velocity is calculated from each echo pulse as explained in the earlier paper. Next, the ranges together with velocities from a single beam pass are combined to a measurement vector for a linear inversion problem. The solution of the inversion problem gives the time-dependent range and velocity from the time interval of satellite flight through the radar beam. The method is demonstrated using the EISCAT (European Incoherent Scatter) UHF radar and radio pulses reflected by a satellite. The achieved standard deviations of range are about 5–50 cm and those of velocity are about 3–25 mm s −1.

The goal of the present work is to use a radar to facilitate accurate satellite orbit determination.When a radar echo signal is sampled, the most straightforward approach for measuring the range would be to recognise the front end of the received echo pulse from the first sample rising above the noise level.Then the range can be calculated from the time difference between transmission and reception.This would lead to an accuracy determined by the sampling interval.A sampling interval of 1 µs, for instance, would lead to a range resolution of 150 m only, and therefore more sophisticated methods are needed for better resolutions.
An elementary way of measuring the velocity would be to observe two ranges from two subsequent pulses and to divide the range difference by the time between the two transmissions.A better approach is to calculate the Doppler shift of the echo pulse by fast Fourier transform.However, for very high accuracies, Fourier transforms with impractical lengths are required.Thus there is a need of finding a way of calculating the Doppler shift more accurately than Fourier transform can do in practice.
Published by Copernicus Publications on behalf of the European Geosciences Union.Nygrén et al. (2012; Paper I from now on) present a numerical method for determining the target velocity at a high accuracy.It consists of two parts.First, the Doppler shift of an individual pulse is found very accurately from the asymmetry of the amplitude peak given by a Fourier transform.The second stage is to adopt a time-dependent velocity model, which allows linear variation of the time derivative of velocity.The model is applied by fitting it to velocities given by subsequent echo pulses from the satellite passing the radar beam.When this method was applied to observations made using the EISCAT (European Incoherent Scatter) UHF radar, the best standard deviations of the velocity of a satellite passing the radar beam were of the order of a few millimetres per second.
The present paper puts forth a method of determining the range and beam-aligned velocity of a satellite from phasecoded radar pulses.The downconversion and scaling of the transmitted and received pulses is first explained and then a method is introduced to determine the time delays between respective phase shifts of the transmitted and received pulses.These delays are used to obtain range measurements with respective standard deviations, one observation for each transmitted pulse.
Each pulse also gives a velocity measurement by means of the method given in Paper I. Range and velocity observations from a single beam pass are collected together and used as measurements in an inversion problem.The inversion problem makes use of a time-dependent range model which is essentially a Taylor's expansion of range up to the third time derivative.The constants in this expansion are the unknowns in the inversion problem.The solution of this problem gives the best values of these constants with standard deviations.Once the constants are known, the time variation of range and velocity can be calculated.
The method is tested with EUMETSAT'S METOP-A (European Organisation for the Exploitation of Meteorological Satellites -Meteorological Operation A) satellite and the EISCAT UHF radar (Folkestad et al., 1983) using alternating codes (Lehtinen and Häggström, 1981).Some checking of the results is made by solving two separate smaller inversion programs.One of them uses only range observations and the other velocity observations as measurements.A good agreement is found between the results from these two independent sets of measurements.

Processing the radar pulses
Data given by the EISCAT UHF radar (Folkestad et al., 1983) are used in this paper to observe a satellite flying through the radar beam.The pulse length is 1920 µs and the signal is sampled at 1 µs intervals as explained in Paper I. A 32-bit alternating code is used in the experiment, and the bit length is 60 µs.Each IPP (inter pulse period) is sampled continuously starting from the very beginning.Hence the transmitted wave form is also recorded.
The transmission and echo signals are used for determining the time delay from transmission to reception at a high accuracy.In principle, this could be done using the time differences of the front or back ends of the transmitted and the received pulse.It turns out, however, that the ends of the transmitted and received pulse are polluted by small disturbances.This problem was detected by our analysis and it was tracked down to experiment design, where the initial and final phase flips where performed too near the pulse-on and pulseoff times.The problem can be corrected and it should not appear in the future.The available data can still be used, since alternating codes contain several phase flips inside the modulation patterns when no power switch on or off takes place.These phase flips can be used instead of front or rear ends of the pulse.Times of each phase flip in the transmission and the corresponding flip in the radar echo give a single range estimate.Thus several range estimates are obtained from a single phase coded transmitted pulse, and this also improves the accuracy of the range determination.
In order to measure the time differences between the phase flips of the echo pulses and the transmitted ones, one has to convert both of them to zero frequency and to scale the amplitudes to the same values.The transmission also contains notable distortions of the ideal alternating code modulation envelopes, and they should be removed.There is an approximately linear amplitude droop of a few per cent during the length of the pulse.In addition, there is a phase droop which is equivalent to a Doppler shift.To correct these distortions, the transmitted waveform after downconversion to zero frequency is modelled by (1) Here A is a constant complex amplitude, β is a real constant modelling the amplitude droop, ν p is a Doppler frequency modelling the phase droop, E(t) is the ideal filtered modulation envelope of the alternating code and ε(t) is noise.The values of the modulation envelope are +1 or −1 according to the alternating code pattern, except close to the phase flips where the filtered envelope varies smoothly.
After the transmitted signal is downconverted to zero frequency and sampled at constant sampling intervals, the model parameters A, β and ν p are determined by nonlinear least squares fitting.A few sample points at the smooth transient near each phase flip are omitted in this fit.When the parameter values are known, the real-valued ideal modulation envelope E(t) is obtained from Eq. (1) by division.An example of the result is shown in the top panel in Fig. 1.
The received pulse consists of the same modulation envelope as the transmitted one.It also contains the same amplitude and phase droops, but its complex amplitude is different.In addition, it has experienced a Doppler shift.The Doppler frequency is determined in the manner explained in Paper I. Using this frequency, the received signal is first downconverted to zero frequency.Then the amplitude and phase droops are removed using the parameter values determined from the transmitted pulse.Finally a real-valued echo pulse is obtained by dividing the echo signal by a complex amplitude determined from the echo pulse itself.This is obtained by averaging 40 samples (the bit length is 60 samples), taken from a region where the modulation envelope does not change its sign and which lies outside of any region of phase flip.An example of the resulting echo pulse is shown in the bottom panel in Fig. 1.One should notice here that the amplitude and phase droop corrections are determined very accurately from the transmitted pulse with a high SNR (signal to noise ratio).When the same parameters are used in correcting the echo pulse, a similar accuracy is achieved, although the SNR of the echo pulse is essentially lower.Another point to notice is that the length of the echo pulse is also changed due to the Doppler effect.There is no need to investigate this, however, because the analysis will deal with the time differences of the corresponding phase flips in the echo pulses and transmitted pulses.
The receiver works as a filter and its filtering properties can be expressed in terms of a single impulse response.Precise measurements of the receiver impulse response are not available but the filter bandwidth gives the length of the impulse response.Here a model is used which consists of a boxcar impulse response of the analogue receiver with a length of 0.3 µs and a triangular impulse response of the digital receiver with a length of 1.2 µs.The total impulse response is the convolution of these two responses and it is shown in the top panel in Fig. 2.
Filtering the radar signal means that the output signal of the receiver is the convolution of the impulse response and the unfiltered signal.When this is applied to the phase flips, the resulting phase shifts are not sharp, but they consist of slopes.These slopes are convolutions of the impulse response and a step function, i.e. they are step responses of the filter.The applied sampling interval is such that it is slightly shorter than the length of the impulse response in such a way that at least one but never more than two of the samples lie within the slope.These are called slope points here.The slope points are shown by dots in Fig. 1.
The theoretical shapes of the filtered phase shifts are shown in the bottom panel in Fig. 2; the panel on the left indicates amplitude shift from −1 to +1 and on the right a shift back to −1.When the sampling interval is 1 µs and the time between subsequent phase flips is a multiple of 60 µs, the slope points at these two types of phase shifts should lie symmetrically around zero.Slope points can be located anywhere within the slopes.If the slope point for one phase shift happens to be the black dot on the left-hand panel, the corresponding slope point for the opposite phase shift should be the black dot on the right-hand panel.The same applies to the pair of open dots.
The top panel in Fig. 1 reveals that the slope points are not actually located symmetrically as expected.Inspection of the signal path at the radar site with an oscilloscope has shown that the asymmetry is most probably due to an  actual (previously ignored) hardware timing difference between generation of the up and down phase shifts.Even the shapes of the two slopes appear to be slightly different.In this paper, we nevertheless assume the ideal shapes shown in Fig. 2, but acknowledge that this is a possible source of a small second-order error in range (first-order timing distortions cancel out, because approximately the same error would be present both in transmission and in reception).

Ranges from individual radar pulses
The times of phase flips can now be determined using the step responses.Assuming that the time dependence of the step response is p(t) and the signal value at a given observed slope point at a time t s is s s , the true time of phase flip is t s − τ s , where τ s obeys s s = p(τ s ).This can be verified using Fig. 2.
If the slope point of a given phase flip has an amplitude equal to that of the filled dot in the bottom left panel in Fig. 2, the true time of the phase flip is 0.5 µs earlier than the time of the observed slope point.
In addition to possible systematic errors caused, e.g. by an inaccuracy in the slope-shape model, there are also statistical errors which are caused by noise.In order to calculate the statistical error of τ s , both the delay τ and the signal value s are treated as random variables.Then where ε is due to Gaussian noise with zero mean.In terms of the SNR of the complex signal, the standard deviation of ε is When τ is fixed, s obeys the conditional probability density D(s|τ ).This is proportional to the conditional probability density D(τ |s), and therefore The Taylor expansion of p(τ ) around τ s is At high values of SNR, σ ε is small so that D(τ |s) makes a narrow peak.Then one can approximate the density function in Eq. ( 4) using the two first terms in the Taylor expansion.Since p is the convolution of the impulse response h and a step function with a step size of 2, Using this result and noticing that p(τ s ) = s s , the conditional density of τ at s = s s is  where the standard deviation of τ is The last expression in this equation is obtained using Eq. ( 3).This result gives an estimate for the the standard deviation of the time of the phase flip, when the signal value at the slope point is measured.The standard deviation is determined by two factors, the SNR and the value of the impulse response at the time corresponding to the time at the slope point.Equation ( 8) shows that, for a given SNR, the determination of τ is most accurate if the time of the slope point corresponds to the maximum value of the impulse response.This, on the other hand, corresponds to the time of the steepest slope of the step response.This is also otherwise an obvious result, because varying s causes smallest variation in τ at the point where s(τ ) has its steepest slope.In this work the amplitudes of the applied slope points always lie within the limits ±0.7.
Alternating codes consist of a cycle of pulses with different modulation patterns, and the number of phase flips in individual pulses is variable.When the pulse front and rear ends are excluded, the number of phase flips in the applied 32-bit code varies between 10 and 21.The times of the phase flips are determined both for the transmitter envelope and for the echo pulse in the manner described above.
If the calculated time of phase flip n is t (TX) n for the transmitted pulse and t (RX) n for the echo pulse, the time delay of the received phase flip from the transmitted one is and the range of the satellite is The time of the reflection of the phase flip at the satellite is Due to the high satellite velocity, the time difference between the received and transmitted phase flips varies from flip to flip.For instance, if the beam-aligned satellite velocity is 4000 m s −1 (see Paper I), the satellite range varies nearly by 8 m during the period of pulse reflection.This value is so large that it should be taken into account in the analysis.Still the pulse is so short that it is sufficient to assume the satellite velocity to be constant during the period of pulse reflection.
Then one can make a linear fit to the ranges calculated from the phase flips of a single echo.
Fitting the calculated ranges is demonstrated in the top panel in Fig. 3.Here the modulation consists of 20 phase flips and the dots indicate the ranges calculated using Eq. ( 9).In order to show the variation, the plotted quantity is R n − r 10 .Here R n is the range calculated from the individual phase flips according to Eq. ( 10 is equal to t 10 = 5433.801µs.The error bars in the top panel in Fig. 3 are standard deviations calculated using Eq. ( 8).The dots without error bars are points which do not obey the criterion −0.7 < s < 0.7 and they are not used in the fitting.The continuous line is the result of a weighted least squares linear fit to the other points.The bottom panel in Fig. 3 shows the standard deviations of the fitted ranges.It turns out that all standard deviations are below 1 m, but the the smallest one is slightly higher than 50 cm and it is obtained for n = 10.The corresponding range and its standard deviation will be used in calculating the range variation of the satellite during the beam pass.

Formulation of the inverse problem
When the target passes through the radar beam, a few hundred echo pulses are typically received.Range and velocity with respective standard deviations can be determined from each pulse separately; range as explained in Sect. 3 and velocity as explained in Paper I. In Paper I, a time-dependent velocity model was fitted to the observed velocities.Here the principle is expanded to a model of time varying range.In order to achieve the best possible accuracy, the range is fitted simultaneously to the observed ranges and velocities.This can be done using the formalism described below.
The time behaviour of range is taken to obey a model This model implies that the beam-aligned velocity has a time dependence and its time derivative is Although a(t) is not the beam-aligned component of the target acceleration, but only the second time derivative of range, it will be called acceleration here.
The task is to determine the best values of the parameters r 0 , v 0 , a 0 and ȧ and their standard deviations.
The number of pulses used in range determination is denoted by n r .Results like those in Fig. 3 are used for obtaining ranges at certain instants of time.The range with the smallest standard deviation and the time of the corresponding phase flip are taken from each echo pulse.The range and time for pulse number i are denoted by r i and t (r)i , respectively.The ranges are collected into a column vector r = (r 1 , r 2 , . . ., r nr ) T , where T indicates transpose.The four unknown parameters in Eq. ( 12) make an unknown vector With this notation, the relation between the unknowns and the measured ranges can, according to Eq. ( 12), written as where The number of pulses used in velocity measurements is not necessarily the same as that for range measurements and even the times of range and velocity measurements from the same pulse are not the same.Therefore the number of velocity measurements is denoted by n v , and the velocities and the times by v i and t (v)i , i = 1, 2 . . ., n v .In analogy with Eq. ( 15), the observed velocities are collected into a column vector The three unknown parameters in Eq. ( 13) are collected to a column vector which contains the three last components of the unknown vector x.Then Eq. ( 13) gives where Equations ( 17) and ( 21) represent two linear inversion problems which could be solved separately.Actually, Eq. ( 21) is the same problem which was presented in Paper I. Instead of solving the two problems separately, a more compact way is to combine them into a single one.This has benefits which will be discussed later.The two vectors r and v are combined into a single measurement vector and the two matrices B (r) and B (v) into a single theory matrix where 0 (v) is a n v × 1 zero matrix.Then the so-called direct theory of the combined linear problems gets a simple form and Then the standard deviations of m can be arranged into a single vector and the diagonal covariance matrix of the measurement vector is The inversion problem in Eq. ( 25) can now be solved following the procedure presented in Paper I. Hence the best value of the unknown vector is and its error covariance matrix is Since the vector x has 4 components, x is a 4 × 4 matrix.
The best values of the unknown parameters r 0 , v 0 , a 0 and ȧ, i.e. r0 , v0 , â0 and â, are the components of x and their variances are obtained from the diagonal components of x .One can now calculate the time behaviour of range from Eq. ( 17) and that of velocity from Eq. ( 21).
The standard deviation of range at each instant of time is given by the square roots of the diagonal elements of the matrix Correspondingly, the standard deviation of velocity at each instant of time is obtained from the diagonal components of the matrix where xv is the 3 × 3 matrix in the bottom right corner of x .
According to Eq. ( 21), the time variation of acceleration can be written as where x(a) = ( â0 , â) T and B (a) is a matrix similar to the first two columns of B (v) .Then the standard deviation for each time is obtained as a square root of the respective diagonal element of the covariance matrix where xa is a 2×2 matrix in the bottom right corner of x .
of flip n = 10 is t 10 = 5433.801µs and the corresponding fitted range is r 10 = 1 Standard deviations of fitted ranges.

Inversion results
The method is demonstrated with data from the EIS-CAT UHF radar, obtained from beam passes of EUMET-SAT's METOP-A satellite.The satellite is 3-axis stabilised and, including its solar panel, it has a size of 17.6 m × 6.7 m × 5.4 m.The method described in Sect. 4 was applied to the same 12 beam passes as those in Paper I. The middle panel in Fig. 4 shows the corresponding standard deviations from individual pulses.The smallest values of these standard deviations are encountered close to the centre of the radar beam and their magnitude is of the order of 0.5 m.The plot also reveals a set of outliers making descending structures.
The standard deviations of the inversion results, calculated according to Eq. ( 32), are plotted in the bottom panel in Fig. 4.These standard deviations are of the order of 5 cm and remain rather constant during the beam pass.The result shows that the inversion improves the standard deviation roughly by an order of magnitude.Figure 5 shows the corresponding results for velocity.The dots in the top panel are the same as the dots in the top panel in Fig. 4 of Paper I, but the red curve is obtained by fitting the model simultaneously to both ranges and velocities.The middle and bottom panels show the standard deviations of measured and fitted velocities, respectively.Also in this case the standard deviations are improved roughly by an order of magnitude in fitting, leading to minimum values of a few millimetres per second.Acceleration is not shown here, since its behaviour is similar to that shown in Paper I.
The analysis made in Paper I revealed that the measured velocities showed semiperiodic variations around the bestfit curve during some beam passes.These variations were clearly not of statistical origin.Therefore, even in this case it is of interest to investigate how the measured values are scattered around the fitted curves.
Figures 6 and 7 show deviations of individual range and velocity measurements from the inversion results.Deviations of ranges are portrayed in the top panels and the corresponding deviations of velocities in the bottom panels.The red bars are errors in terms of standard deviations.
Figure 6 contains results from beam pass 6.The range observations in the top panel contain semiperiodic ascending structures with deviations of the order of a few metres.These deviations are too large in comparison with the standard deviations of the observation.The gaps in time come from pulses which violate the amplitude criterion ±0.7 explained in Sect.3. Also in the bottom panel, the scatter of individual velocity observations around the fitted value looks too broad in comparison with the standard deviations.This suggests that both range and velocity observations contain errors which are not of statistical origin.The bottom panel shows no variations similar to those in the top panel.Still, the observed velocities may follow some semiperiodic variation which is partly invisible due to the scatter of the individual points.In some other beam passes clearer periodic variations are seen, as shown in Paper I. fitted curves seems too broad in comparison with their standard deviations.
There are two possible reasons for the observed broad distributions.One is that the true slopes at the phase shifts are different from the applied model.Another point to consider is the size of the satellite and its geometrical structure.The signal model assumes a point target, whereas the satellite is several metres in size and it carries different instruments and a large solar panel.When the whole satellite is illuminated by the incident wave, reflections from different parts of the satellite are summed up in the received signal.Then the exact meaning of the calculated range is not clear at all.In Paper I it was assumed that the periodic velocity variations would be due to the satellite rotation but, since the satellite is 3-axis stabilised, this explanation does not seem to be valid.Whatever the reasons are, the investigation of plots like those in Figs. 6 and 7 shows that the true errors or range are not likely to be higher than 1 m, which is clearly smaller than the satellite size.Also, the errors of velocities must usually be smaller than some centimetres per second.
Results from analyses made for all 12 beam passes are displayed in Table 1.The table contains ranges, velocities and accelerations at a single instant of time for each beam pass, together with their standard deviations.The reference times are chosen at points where the standard deviation of range is smallest.The standard deviations of velocity and acceleration are typically of the order of a few mm s −1 and a few mm s −2 .

Separate inversions
Equations ( 17) and ( 21) define two inversion problems which can be solved separately.In the first one, the measurement vector consists of ranges and in the second one the measurements are velocities.The error covariance matrices of these inversion results can be calculated in the same way as for the complete inversion problem, i.e. in the way similar to Eq. ( 31).The two sets of measurements in the separate inversion problems are completely independent, since ranges are determined from delay times and velocities from Doppler shifts.When results from the two separate inversion problems are compared with those given by the joint inversion problem, one can evaluate the benefit of the method applied in this paper.In addition, both methods are able to give fitted velocities and accelerations.Since they come from independent measurements, the reliability of the method can be checked by comparing the results.
Figure 8 shows the difference of range r r given by the range inversion problem in Eq. ( 17) and range r given by the complete inversion problem in Eq. ( 25).The top panel comes from beam pass 1 and the bottom panel from beam www.ann-geophys.net/31/859/2013/Ann.Geophys., 31, 859-870, 2013 Table 1.Results from 12 beam passes on 1 December 2010, after 16:00 UTC (Universal Time Coordinated).The numbers in the first column refer to the beam passes in Table 1 of Paper I. The next three columns show the beam azimuths and elevations as well as the maximum SNR.
The following column gives the start times of the IPPs in minutes and seconds, t gives the times of the results in the next columns (after the start of the IPPs), and the last six columns contain the ranges, the line-of-sight velocities and the accelerations with their standard deviations.
No. az el SNR mm:ss In the middle of the radar beam the agreement is quite good.
In the case of beam pass 1, for instance, r r − r = 0.2 m.At the edges of the beam the difference may be of the order of a few metres.It would be expected that the difference always departs from zero by less than the given error limits.Indeed, this mostly happens but not at all times.Figure 4 shows that the standard deviations of r from beam pass 1 are of the order of 50 mm and nearly constant during the whole beam pass.
Close to the centre of the beam, the corresponding standard deviations of r r (not shown) are only about 26 % higher than those of r at the minimum of standard deviations, but at the edges they are about 10-fold.Results for beam pass 6 are quite similar.Hence including velocity measurements in the inversion improves the accuracy of range only slightly close to the centre of the beam, but the improvement during the whole beam pass may be substantial.Velocities obtained by different kind of inversions are compared in Fig. 9.All these results are from beam pass 1.The top panel shows the difference of velocities obtained from velocity observations (v v ) and those obtained by combining both velocity and range observations (v).The difference is maximally of the order of a couple of millimetres per second.
The ratio of the standard deviations given by the two methods is plotted in the middle panel in Fig. 9.The results shows that including the range measurements in the inversion only very slightly improves the accuracy.Thus the improvement in velocity gained by the inversion problem in Eq. ( 25) in comparison with that in Eq. ( 21) is nearly negligible.
The bottom panel in Fig. 9 shows the difference of velocity v r obtained from range measurements and v, i.e. that given by combining both velocity and range observations.The former is obtained from Eq. ( 13) with the parameter values of v 0 , a 0 and ȧ given by the inversion problem in Eq. ( 17).The grey region indicates the statistical error of v r in terms of 3 σ limits.It is found that the velocities given by the two methods are well in agreement.Close to the beam centre the results are nearly identical.For instance, at the time of minimum standard deviation of range, v r − v = −9.3mm s −1 and this departs from zero by less than its 3 σ limit.In most of the 12 beam passes this difference is smaller than 0.5 m s −1 .
Figure 10 shows a similar comparison of acceleration given by different inversions, also from beam pass 1.The top panel gives the difference of accelerations from velocity measurements (a v ) and all measurements (a).The results are nearly identical, the difference being always smaller than 1 mm s −2 .The second panel indicates that the ratio of the two standard deviations is almost precisely equal to unity, so that same standard deviations are given by both methods.The bottom panel demonstrates the acceleration a r obtained by range measurements and the grey region displays the 3 σ limits.At the edges of the beam, the difference may be of the order of 2 m s −2 and, even at the time of minimum standard deviation of range, a r − a = 0.52 m s −2 .At this instant of time, the respective standard deviations are σ ar = 107 mm s −2 and σ a = 2.7 mm s −2 .Considering that a ≈ 18 m s −2 , the disagreement between a and a r is rather large in comparison with the corresponding disagreements of range and velocity.

Summary and discussion
This paper, together with the earlier work by Nygrén et al. (2012), presents a compound algorithm for determination of satellite velocity and range from radar pulses observed during a beam pass of the satellite.The method first calculates the Doppler shifts of individual echo pulses at accuracies beyond what can be obtained by Fourier transform of any practical length.The beam-aligned velocities are obtained from the Doppler shifts.This method is described by Nygrén et al. (2012).Next, the time differences of corresponding phase flips in the transmitted and received phase coded pulses are determined.These time differences give range values at times when the phase flips hit the target.The method is based on the filtering effect on the pulse shape.The receiver filter creates slopes at the phase flips and the amplitudes of the samples at the slopes allow the determination of the times of the phase flips at a high accuracy.
The ranges and velocities from a single beam pass of the satellite are combined into a stochastic inversion problem.Solving the problem gives a Taylor's expansion of range as a function of time up to the third time derivative.Thus the result gives the time variation of range and velocity as well as the linear time variation of the second time derivative of range, together with their standard deviations.
The method is demonstrated using data from the EIS-CAT UHF radar, obtained from beam passes of EUMET-SAT's METOP-A satellite.The smallest standard deviations of range from the analysed 12 beam passes vary from 5 to 50 cm and those of velocity from 3 to 24 mm s −1 .These are the statistical error limits given by the applied model assuming a point target.The best results shown in this paper compare well with those obtained by lasers, altimeters and GPS receivers.For instance, Tapley et al. (1994), Nouël et al. (1994) andVisser et al. (2009) report position errors of the order of 5 cm or somewhat below.
The target satellite is 3-axis stabilised and, including its solar panel, it has a size of 17.6 m × 6.7 m × 5.4 m.Thus the standard deviation of range is much smaller than the satellite size.Since the whole satellite is illuminated by the incident wave, it is impossible to know what part of the satellite the inversion result refers to.The analysis also shows that the observed deviations of individual range observations from the fitted range are often larger than what one would expect from their standard deviations.As discussed in Sect.5, this is possibly caused by the analysis method creating non-statistical errors.The results indicate that the true error of the fitted results is not likely to be much higher than 1 m.This is still smaller than the satellite size.
The determination of the true range is affected by a factor which is not taken into account in this paper.The point of transmission should be defined at an accuracy which compares to the standard deviation of the calculated range.This is a problem itself and it is not considered here.
Another fact to be taken into account is that the direction of the satellite position cannot be determined precisely by the monostatic radar experiment used in this paper.This is because the results of the EISCAT UHF radar come from the beam with a width of the order of 0.5 • .The variation of SNR during the satellite pass gives some clue of a more precise direction but, because the satellite hardly ever flies through the centre of the beam, even this information cannot be very precise.A tristatic radar would offer a possibility of more accurate determination of satellite position.By directing all three antennas to the satellite, one could find the ranges from the three antennas, and then the position and velocity would be unambiguously determined.
Since the radar pulse travels in the ionosphere, the radio wave may be sightly affected by refraction and then the measured distance does not exactly correspond to a straight line between the radar and the satellite.When extremely precise measurements are made, the ionospheric effects are not necessarily negligible (Hoque and Jakowski, 2011).These effects are largest at low elevations of the radar beam, but they are not considered in the present work.
When considering the velocities, their true errors must also be larger than those shown by the standard deviations.This is seen by comparing the individual standard deviations and the scatter of observations around the fitted velocities.Also, as shown in the figures of the earlier paper (Nygrén et al., 2012), the velocities sometimes show periodic variations around the fitted values.Whatever the reason is, the true errors of velocity are usually not higher than a few centimetres per second.It is possible to divide the inversion problem into two separate parts, one of them using only range observations and the other only velocity measurements.Comparison of results from these inversion gives a view on the benefits gained by the full inversion.The results show that the full inversion problem only slightly improves the accuracy in the middle of the radar beam, where smallest standard deviations are obtained.However, a clear improvement is achieved at the edges of the beam, and therefore the results given by the full inversion problem are useful when applied in determining the orbit parameters, for instance.
The separate inversion problems allow the determination of velocity separately from measured ranges and from Doppler shifts, which are completely independent observations.The standard deviations given by the former method are larger than those given by the latter one, but the agreement between the two results is good.Closest to the centre of the radar beam, the difference is typically some tens of centimetres per second and rarely more than a few metres per second.The difference usually lies within a 3 σ limit of the least accurate result.This agreement between the two results greatly supports their reliability.
Ultimately, the accuracy of our method must be verified against the operational orbit of the observed satellite.Before this can be meaningfully done with the inherently very high accuracy provided by the present method, the systematic offsets and distortions due to the receiver and the antenna need to be carefully quantified.This will require further inspection and measurements of the radar system itself.

Fig. 1 .Fig. 1 .
Fig. 1.Top panel: An example of a sampled transmitted pulse after downconversion to ze and conversion to a real signal.Bottom panel: Same for an echo signal.The open dots in phase flips.

Fig. 3 .
Fig. 3. Top panel: Linear fit of ranges from the phase flips of a single echo pulse.The d calculated from individual phase flips and the error bars are their standard deviations.Th of flip n = 10 is t 10 = 5433.801µs and the corresponding fitted range is r 10 = 1315.4058Standard deviations of fitted ranges.

Fig. 3 .
Fig. 3. Top panel: linear fit of ranges from the phase flips of a single echo pulse.The dots indicate ranges calculated from individual phase flips and the error bars are their standard deviations.The time of reflection of flip n = 10 is t 10 = 5433.801µs and the corresponding fitted range is r 10 = 1315.4058km.Bottom panel: standard deviations of fitted ranges.

Fig. 4 .Fig. 4 .
Fig. 4. Ranges and their standard deviations from satellite pass 1. Zero time tim measurement for this specific beam pass.Top panel: The red curve is the fitted r trend r L and the dots are ranges obtained from individual pulses.Middle panel: from individual pulses.Bottom panel: Standard deviations of fitted ranges.
Figure 4 portrays ranges from satellite pass 1.The red curve in the top panel indicates the ranges given by the inversion, after subtracting a linear trend r L .The dots indicate ranges from individual pulses, calculated as described in Sect.3.

Fig. 5 .
Fig. 5. Velocities and their standard deviations from satellite pass 1. Zero time time refers to the start time of the measurement for this specific beam pass.Top panel: The red curve is the fitted velocity after subtracting a linear trend v L and the dots are velocities obtained from individual pulses.Middle panel: Standard deviations of velocities from individual pulses.Bottom panel: Standard deviations of fitted velocities.

Fig. 6 .Fig. 5 .
Fig. 6.Deviations of individual measured ranges r m (top) and velocities v m (bottom) from the inversion results r and v.The standard deviations are indicated by red bars.The observations are from beam pass 6.

Fig. 6 .Fig. 6 .Fig. 7 .Fig. 7 .
Fig. 6.Deviations of individual measured ranges r m (top) and velocities v m (bott r and v.The standard deviations are indicated by red bars.The observations are

Fig. 8 .Fig. 8 .
Fig. 8. Differences of ranges obtained by solving the complete inversion problem and the inversion problem for mere range measurements.The results for the two methods are denoted by r and r r , respectively.The grey areas indicate the 3σ limits of r r .The top panel is from beam pass 1 and the bottom panel from beam pass 6.

Fig. 9 .Fig. 9 .
Fig. 9. Comparison of velocities from beam pass 1, obtained by three differen ferences of velocities given by the complete inversion problem and the inversi measurements.The results for the two methods are denoted by v and v v , respect of the standard deviations of v v and v, denoted by σ vv and σ v , respectively.velocities given by the complete inversion problem and the inversion problem f The latter results are denoted by v r .The grey areas indicate the 3σ limits of v r .

Fig. 10 .Fig. 10 .
Fig. 10.Comparison of accelerations from beam pass 1, obtained by three different inversions.Top panel: Differences of accelerations given by the complete inversion problem and the inversion problem for mere velocity measurements.The results for the two methods are denoted by a and av, respectively.Middle panel: The ratio of the standard deviations of av and a, denoted by σav and σa, respectively.Bottom panel: Differences of velocities given by the complete inversion problem and the inversion problem for mere range measurements.The latter results are denoted by ar.The grey areas indicate the 3σ limits of ar.

www.ann-geophys.net/31/859/2013/ Ann. Geophys., 31, 859-870, 2013
The standard deviations of the measurements are needed in solving the linear inversion problem.The standard deviations of ranges and velocities are collected, respectively, into column vectors ) Ann. Geophys., 31, 859-870, 2013 www.ann-geophys.net/31/859/2013/ The difference is plotted by heavy lines and the grey regions indicate statistical errors of r r in terms of 3 σ limits.