International Laser Ranging Service

ILRS Normal Point Algorithm

Data Screening and Normal Point Formation

Re-Statement of Herstmonceux Normal Point Recommendation
A. T. Sinclair/RGO
June 1997

---

2012 Amendment:

The original Herstmonceux normal point definition specifies a standard normal point interval (SNPI) for each satellite based on altitude. This definition was adopted based on our then current firing rates of 5 – 10 Hz. Newer systems are now operating from 100 Hz to kHz rates and achieve sufficient full-rate data to accurately define normal points in time intervals far less than the SNPI. An amendment has been drafted to relax some of the normal point restrictions to allow the high repetition rate SLR systems to complete normal points in less than the SNPI and then move on to another satellite, thereby increasing data pass yield through pass interleaving.

---

Note:

The ILRS continues to strive for improvement in data reliability. At its meeting in Nice 2002, the ILRS Governing Board adopted a minimum data requirement for SLR normal points. This requirement, adopted at the recommendation of the Networks and Engineering Working Group is as follows:

Daytime normal points - minimum 6 data points

Night time normal points - minimum 3 data points

Fewer data points would be acceptable on lower satellites (5-second normal points) from those ranging systems with lower pulse repetition rates where these minimum requirements are not practical.

---

At the Laser Ranging Instrumentation Workshop held in Herstmonceux in 1984 a proposal was made for the method to be adopted for screening SLR data and for forming normal points. This proposal was published in the Proceedings of the Workshop, and later in the CSTG SLR Newsletter Vol 2, No 1 (October 1987). Unfortunately both of these descriptions contain some typographical errors, and there are some aspects of the wording which need clarification. Also some other small modifications are now needed in view of increased understanding of satellite signature effects and some signature effects caused by avalanche photo diode detectors. The purpose of this note is to carry out these revisions.

By about 1990 the improving precision of SLR stations enabled analysts to begin to see effects of satellite signature, which cause a non-symmetrical distribution of the ranging data from some satellites. This is particularly noticeable for stations operating at single photon return level. Also the characteristics of some SLR systems cause an intrinsic non-symmetrical distribution of the data, and this again is particularly marked for single photon systems. The method of data screening and the adopted mean reference point recommended in the Herstmonceux proposal are based on the assumption of an approximately Gaussian distribution of the data, and so need to be reconsidered for non-symmetrical data. An early note on this topic is by G.M. Appleby & A.T. Sinclair in CSTG SLR Newsletter, April 1991. In March 1991 a working group was set up to consider this matter, consisting of A.T. Sinclair, P.J. Dunn, R.J.Eanes, T.K. Varghese and G.M. Appleby. The report on its deliberations is presented and enlarged upon in a paper by A.T. Sinclair in the Proc. 8th Int. Laser Ranging Workshop, Annapolis 1992, entitled 'SLR Data Screening; location of peak of data distribution.' It points out that for a non-symmetrical data distribution the mean reference point proposed in the Herstmonceux algorithm will be displaced from the peak of the data. It considers that the peak might be the better reference point to adopt, and describes some ways of determining the peak of a distribution.

At the 9th International Laser Ranging Workshop in Canberra in 1994, papers were presented by J.J. Degnan and by R. Neubert which calculated theoretically the reflection signature of Lageos-type satellites, and calculated the correction from the mean reflection point to the centre of mass of the satellite. Subsequently further work was done by R. Neubert and by G.M. Appleby to calculate theoretically the correction from a reference point to the centre of mass, for a variety of ways of determining the reference point (e.g., mean, peak, leading-edge half-maximum), and for a variety of types of SLR system (e.g., single or multiple photon return level, and different detector types), and further work was done by A.T. Sinclair on ways of processing the data in order to determine these various reference points. These results were presented and discussed in depth at a meeting organised by Eurolas in Munich in March 1995, and the written papers are available in the Proceedings of that meeting.

This note re-states the basic Herstmonceux algorithm for forming normal points, with the typographical errors in earlier descriptions corrected (we hope!). It then describes the calculation of other quantities that define the distibution of the data, these being the skewness, kurtosis and the location of the peak relative to the mean of the data.

Re-statement of the Herstmonceux Algorithm

Formation of range residuals from a trend function, and data screening.

(1) Use high precision predictions to generate prediction residuals PR, where PR = observation - prediction. Include the best available estimate for a time bias for the predictions and possibly a correction to the UT1 value used to generate the predictions. The predicted range must include the refraction delay.

(2) Use a suitable range window to remove large outliers.

(3) Solve for a set of parameters (orbital parameters are preferable) to remove the systematic trends of the prediction residuals, while not introducing spurious high frequency signals (as can occur by fitting a high-order polynomial). The resulting fitted function is the trend function f, with value f(PR) at the epoch of the residual PR.

(4) Compute fit residuals FR, where FR = PR - f(PR). Compute the root-mean-square value RMS of the fit residuals. Identify outliers using a rejection level of nxRMS (where n is defined below), and omit these in step (3) of the next iteration.

(5) Iterate steps (3) and (4) until the process converges. Previously rejected data should be reconsidered in step (4) at each iteration.

Appropriate values for n are :

n = 3.0 for systems using multiple photon detection systems

n = 2.5 for systems that detect the first photo-electron.

The reason for using a tighter rejection limit for systems that detect the first photo-electron (e.g., PMT or APD detectors) is that some of these detectors give a data distribution that is skewed towards long ranges. Also, if the system is operating at single photon return levels, then the effect of satellite signature can also cause skewness towards long ranges. The use of the tighter rejection limit reduces the effect of this skew tail of the distribution.

Normal point formation

(6) Subdivide the accepted fit residuals FR into fixed intervals (bins) starting from 0^h UTC. The bin sizes for various satellites are listed at the end of this note.

(7) Compute the mean value and the mean epoch of the accepted fit residuals within each bin i. Let n_i be the number of accepted fit residuals within this bin.

(8) Locate the particular observation O_i with its fit residual FR_i , whose observation epoch t_i is nearest to the mean epoch of the accepted fit residuals in bin i.

(9) The normal point is computed as

(10) Compute the RMS of the accepted fit residuals in each bin i from their mean value :

where the summation over j is over the accepted points within the bin.

If there is only one point in the bin then take RMS_i = pass RMS.

(11) Report in the normal point data record :

ti	epoch of the normal point
NPi	the normal point range value
ni	the number of observations in bin i
RMSi	the RMS of the observations within the bin from their mean

Comments

1. The algebraic form in which this algorithm is expressed rather hides the underlying principle, which is very simple. The quantities FR are the variations of the individual observations from the trend function, and if it has fitted well then these are essentially just the random errors of the observations. At step (9) the random error of a particular observation O_i is subtracted from it, and the mean of the random errors of all the observations within its bin is added to it. The result is that the random error of the normal point is reduced to that of the mean of the bin rather than that of an individual observation. Also, because the trend function at the observation is subtracted and the mean trend function over the bin is added back, the normal point is to a large extent independent of the particular trend function selected, provided the fit is sufficiently good that there is not much variation of the trend function from the observation residuals over the bin.

2. This re-statement of the algorithm differs from the original in a few details. The rejection level specified is either 3.0xRMS or 2.5xRMS, whereas in the original only 3.0x was specified. RMS is now used instead of as it is better defined for a distribution of data that may not be Gaussian. In step (10) the RMS of the residuals within each bin is calculated. The original algorithm used a divisor (n_i - 1) instead of n_i within the square root, in order to give the 'standard error of a single observation'. This caused some confusion of meaning, and so the RMS is now preferred.

Analysis of the distribution of the data

The parameters skewness and kurtosis are useful for describing some of the features of the distribution of the data. Use only the fit residuals that have been retained within the nxRMS screening in steps (3) and (4).

Let x_i ( i =1, n ) be the screened fit residuals, with mean x_m

Form 2nd, 3rd, 4th moments of the data :
for j=2, 3, 4.

(The 2nd moment ₂ is just the RMS of the data, formed at steps (3) and (4)).

Then,

Skewness = ( = 0 for a symmetrical distribution )

Kurtosis = ( = 3 for a Gaussian distribution )

In some texts the expression used for skewness is the square of this expression, but that loses information on the sign, which for this formula is positive if the data are skewed towards long ranges. If the rejection level used is 2.5xRMS then the value for the kurtosis will probably be about 2.6 rather than 3.0.

If the data have insignificant skewness, then the mean value of the retained fit residuals will be located at the peak of the distribution of the data. However if the data are significantly skew, then the mean will be offset from the peak, in the direction of the skewness. The difference of the mean from the peak is another useful indicator of the distribution of the data. For the Lageos satellites the offset of the mean from the peak due to satellite signature is about 4mm.

There are several possible methods of determining an estimate of the location of the peak. A simple graphical method is to plot a histogram of the distribution of the data. However the results of this are rather dependent on the histogram bin width chosen, and the method cannot give a resolution much better than the bin width. A simple numerical method is to determine a further iterated mean with a very tight rejection level, such as 1.0xRMS, but it is important that the RMS should be held fixed at the value determined in step (4), otherwise this process sometimes fails to converge. A better method is described in a paper 'Data screening and peak location' by A.T. Sinclair in Proceedings of Annual Eurolas Meeting, March 1995. A standard subroutine, DISTRIB, for determining the mean and peak of a data distribution, is available from the Eurolas Data Center.

At present only the mean and the RMS from the mean are reported in the normal point data record, but at some stage these other statistics will be also included, and meanwhile it will be advantageous for a station to compute these parameters, skewness, kurtosis and mean minus peak, for its own information.

Recommended Normal Point Bin Widths

(Editor's Note: These values can be found in separate documentation.)

Related information:

• Overview of Normal Point Data