Introduction
Tracking is a process whereby the position of an interface is determined by an automatic or semi-automatic method. It is applicable to many glaciological data, especially altimeter data, or Radio Echo Sounding (RES) data. The problem of tracking may be divided into several parts, as follows (Fig.1). First, an algorithm must be provided to enable the program to locate the interface in the digitized trace. Secondly, the program must be provided with a means of predicting where in the trace the interface is likely to be found on the following trace, in order to eliminate unnecessary searching, to eliminate glaciologically impossible points and, in the case of instruments with a narrow window (e.g. satellite altimeters), to position the window so that the surface lies within it. Thirdly, the program must be able to recognize when no point in the area of the trace defined above fulfils the criteria specified by the first algorithm. Finally, a means of re-acquiring track after its loss must be provided. In addition to these four processes that are specific to the tracking problem, it will almost certainly be necessary to pre-process the signal to reduce the noise level as far as possible, and to convert the data from instrumental units to geophysical units.
Pre-Processing
Many techniques are available for reducing the noise level in a signal, and the final choice of method will be dependent on the characteristics of the data being processed. If the data contain isolated spikes, a spike-editing procedure such as the median filter (Reference SySy 1985) should be used to remove them. Otherwise, a very useful technique is that of low-pass filtering (e.g. Reference AlbertAlbert 1986), provided it is known that there is a limiting frequency above which there is unlikely to be useful information in the signal, as is the case for pulsed radar systems. There are many other useful techniques such as deconvolution, Wiener filtering and maximum entropy deconvolution which are of use, especially when the signal being processed may be regarded as the convolution of a series of spikes with a wavelet, as in seismic processing.
It may also be useful at this point to convert the data from the instrument units to geophysical units. This is often of benefit if the relationship between the two is non-linear, as, for example, in the case where an amplifier with a logarithmic response is being used.
Locating the Track Point
Various algorithms are in use for specifying the point in the recorded trace that is to be regarded as representing the interface to be tracked. These depend on the characteristics of the surface being tracked, and also on the signal returned by the instrument in question. For example, the surface of the sea is tracked by satellite radar altimeters using a simple model of the expected radar return, which has the advantage of allowing various useful properties of the pulse to be extracted directly from the tracking algorithm. This approach is very successful for returns from ocean surfaces, which are statistically homogeneous and geometrically simple. Unfortunately, for ice surfaces the situation is not so clear cut, because the surface is geometrically complex and because the scattering mechanism from the ice surface is not well understood. This prevents the use of simple models of the return pulse shape, and other approaches have to be considered. The most useful point to track in the return pulse is the first arrival time, as this is the simplest point to interpret in terms of the overall shape of the surface being studied (Reference HarrisonHarrison 1970). This point can be detected either by setting a threshold level and testing for the first point in the pulse to exceed this level, or by differentiating the return pulse and testing for positive-going zero crossings, which represent minima in the returned signal. In either case, noise in the signal will produce spurious results, and secondary criteria such as testing the next maximum in the signal must be used to reject invalid points. The first method is suitable for signals where there is only one interface, or for signals where the interfaces are separated by intervals where the signal level returns to the noise level, and the second is more suitable where the returns from successive interfaces overlap, as is the case in RES when the bedrock return arrives during the wide-angle reflections from the surface. Fig.2 illustrates the difference between the two types of signal described above, showing possible alternate track points.
Predicting the Next Track Point
Prediction of the position of the interface in the record succeeding the current one is the heart of the tracking problem. It is dependent on the geophysical characteristics of the interface being tracked, and also on instrument characteristics, especially the beamwidth in the case of the ranging instruments which are primarily being considered here. A further constraint on the type of algorithm to be used is the time available for the processing of the signal, which is clearly more limited in the case of a satellite altimeter which is tracking in real time than in the case of post-processing RES data.
At this point, it is necessary to discuss the exact meaning of the terms “broad” and “narrow” beamwidths. In this paper, the beamwidth of the instrument is considered to be broad if the nearest point of the surface lies within the antenna beam pattern, and narrow if it lies outside. If we consider a uniformly sloping surface, with undulations superimposed on it, then the beamwidth is broad if the height of the undulations above the sloping surface does not exceed
(This equation is derived in the appendix)
where hmax is the maximum height of an undulation above the sloping surface, R is the range of the instrument from the intersection of the slope with the nadir, Φ is the angle of slope of the surface, ⊖ is the beamwidth of the instrument, and Re is the radius from the centre of the Earth to the nadir point. Note that the last two terms of the equation correct for the curvature of the Earth, and are negligible for airborne instruments.
In the case of broad beam instruments, the maximum rate of change of range to the surface with distance along track is 1.0 (Reference HarrisonHarrison 1970). This means that it is only necessary to scan a portion of the digitized signal equivalent to ± the along track distance between successive returns from the instrument. In the case of airborne instruments where the signal is being digitized fairly frequently, this limits to manageable proportions the volume of data to be inspected for a particular interface. This constraint is not true for narrow beam instruments, but in these cases the digitized pulses returned by the instrument will be difficult to intepret, as the first arrivals will be coming from a region on the extremity of the antenna beam pattern.
For instruments which need a window to be positioned prior to the digitization of the returned signal, predictive tracking methods are commonly used, of which the most widely known is the alpha-beta tracker implemented in the Seasat altimeter (Reference RapleyRapley and others 1983). In this method, the range is predicted using the equations
where Hn is the range estimate at the nth, α is a constant, ΔHn is the error in the range estimate at the nth point, H is the estimate of the rate of change of range at the nth point, β is a constant and Δt is the time interval between pulses. The time taken by the tracker to respond to changes in the surface gradient, or to step changes in the surface elevation, may be adjusted by varying a and β. For Seasat, α was set to 0.25 and β was set to 0.015625, giving a response time of 0.28 s.
In these methods, assumptions are made concerning the form of the surface being tracked, and these assumptions are used to predict from past statistics the future position of the interface. Unfortunately, the major assumption of this method is that the differential of the surface observed varies smoothly, without discontinuities. Reference HarrisonHarrison (1971) showed that this is not a valid assumption for ice surfaces, or indeed for any undulating reflecting surface where the radius of curvature of the surface is less than the range from the sensor. Therefore such methods cannot be depended on over ice or other complex surfaces.
As the surface slopes of more than 90% of the Antarctic continent are less than the half beamwidth of altimeters like the Seasat altimeter (Reference McIntyre and DrewryMcIntyre and Drewry 1984), such altimeters are effectively operating in a broad beam mode and therefore returns may not change their range at a higher rate than 1.0. This gives rise to an interesting possibility for tracking such surfaces without toss of lock. To outline one possible scheme, the position of the returned waveform within the window would be determined on a per pulse basis, and the position of the window would be changed so that the leading edge of the returned pulse is in the centre. Then, as long as the rate of advance per pulse of the satellite does not exceed the half width of the range window, the leading edge of the next pulse must appear within the window so placed, and the process may then be repeated. This method would require fast computing, as only about 1 ms would be available for the positioning algorithm. However, in the case of Seasat the rate of advance per pulse was ~7 m, and the range window was 30 m wide. Thus the range window was sufficiently wide to allow for more than twice the maximum possible shift in the position of the return, provided that the window was centred on the return from the previous pulse. The offset centre of gravity tracker developed by Reference Wingham, Rapley and GriffithsWingham and others (1986) would be a very suitable algorithm for finding the position of the leading edge, as it is insensitive to noise and changes in the pulse shape, and operates linearly over a wide range of displacements from the centre of the instrument window. It should also be borne in mind that modern “off the shelf” microprocessors are capable of running at rates in excess of 1 million instructions per second, and that special purpose processors are much faster. In order to provide averaged pulses for the determination of significant wave height and back-scatter coefficient, it would still be necessary to maintain a value of the range rate, but this could be computed from a smoothed value of the corrections applied to the position of the range window.
Determining loss of Lock
This can be a considerable problem for real-time algorithms, especially over complex surfaces such as ice sheets. Off-line, it is considerably simpler, as more time is available, and it is also possible to include human intervention in the loop. Indeed, for RES, human intervention is the only wholly reliable method of determining the loss-of-lock condition, as it is possible for tracking software to be misled by strong diffraction hyperbolæ, which are indistinguishable from a very steeply sloping interface and are often stronger reflections. In the case of satellite altimeters, criteria based on threshold values of the power difference between the early part of the return and the late part of the return are used. These can easily be invalidated by pulses of unusual shape, so that actual loss of lock occurs much earlier than it is detected by the instrument, as is seen in the case of Seasat data over Antarctica, or over sea ice (Reference RapleyRapley and others 1985).
Acquisition of Lock
The acquisition or re-acquisition of the track point after loss of lock is an essential part of any complete tracking system, and, in real-time systems, must be as fast as practicable, in order to minimize the volume of data lost during this phase. It must also be robust, in order to avoid spurious results. As with the loss-of-lock algorithm, the only method that is wholly reliable is human intervention, which can only be used in post-processing of data gathered in the field. In the case of satellite altimeters, it is usual to operate the altimeter in a completely different mode from the normal tracking mode (Reference RapleyRapley and others, 1983).
Implementation of a RES Tracking Scheme
Following the acquisition of a large volume of digitally recorded RES data in 1983, with further data from 1985 and 1986 (Reference Drewry and LiestølDrewry and Liestøl 1985; Reference Gorman and CooperGorman and Cooper 1987), various tracking programs have been implemented by the author, both on a large mainframe computer and on Z80 based microcomputers. The mainframe has the advantage of processing the data much faster, so that relatively complex data processing can be carried out, but has the serious disadvantage of not allowing manual interaction with the data processing, so that direct human intervention in the case of loss of lock and re-acquisition of lock cannot take place. In many ways, processing this data on a microcomputer has proved more successful, as an operator can intervene in the processing whenever appropriate, and hence simpler tracking algorithms can be used. The main disadvantage is the slow speed of the microcomputer, especially when floating-point arithmetic is to be carried out. Therefore little pre-processing has been carried out, as the rate of processing then drops to unacceptable levels. With more advanced microprocessors than a Z80, there is less of a problem, and experiments will be carried out in the near future using an Intel 80286 fitted with an 80287 maths co-processor.
The scheme followed on microprocessors has been that outlined in Fig.3. For each record, the data are first displayed. This has been one of the least satisfactory parts of the microcomputer implementation, as the resolution of the graphical display has not been adequate to allow the data to be displayed as a profile along the track, with power modulating the intensity of the display as in a Z-scope display. The alternative adopted has been to display successive records as plots of delay time against power, as in an A-scope display. This allows the full resolution of the data to be displayed, and the location of the interfaces can be seen quite clearly. The disadvantage is that the continuity of interfaces cannot readily be visualized. Then the program checks whether an interface was found on the last record (the “track lost flag” referred to in Fig.3), and also whether there has been any keyboard input from the operator. In either case, the operator is given the opportunity to position the track point manually, or, alternatively, to set the “track lost flag” manually, thus preventing mis-tracking by the automatic routines. If the “track lost flag” is set, the results are written to an output file, with identifying data. Otherwise, the interface is located precisely by searching within a window, whose width is constrained as described above, for minima in the signal. In this particular implementation, an interface was defined as occurring when the returned power had risen by one twentieth of the difference in power between the minimum and the next maximum, in order to avoid difficulties with instrumental noise near the minimum. This approach was required because there was instrumental noise near the minimum, and could have been improved on if more pre-processing had been possible. The results were consistent, with a small systematic offset which was determined by inspection of the data. Various constraints were added to the simple test described here, such as a minimum difference in power between the minimum and the maximum, which improved the selectivity of the system. If no suitable interface was found, then the “track lost flag” was set. In either case, the result was written out to file.
Acknowledgements
Much of this study was funded by NERC grant GR3/4463 to D J Drewry. M R Gorman, J A Dowdeswell and W G Rees all contributed much helpful discussion during the development of tracking programs for RES.
Appxendix A
The equation describing the broad beam criterion is derived in this appendix.
If A is the sensor location, C is the nadir point at the surface, and CE is the mean surface, with a slope of Φ, then for a sensor beamwidth of Θ, DE represents the maximum elevation of a point above the slope for the closest approach of the point to be within the beam-limited footprint of the sensor. We can calculate DE as follows:
therefore
therefore
therefore
therefore
therefore
In order to correct for the curvature of the Earth, we have to add the further correction SF:
but
therefore
and therefore
The total correction, DE + SF, is given by:
This equation can now be converted from the geometric notation used here into the algebraic notation used in the main text (Equation 1), as follows:
where hmax(= DE + SF) is the maximum height of an undulation above the sloping surface, R (= AC) is the range of the instrument from the intersection of the slope with the nadir, φ is the angle of slope of the mean surface, θ is the beamwidth of the instrument, and Re (= OC) is the radius from the centre of the Earth to the nadir point.