Hostname: page-component-cd9895bd7-gxg78 Total loading time: 0 Render date: 2024-12-22T14:58:40.508Z Has data issue: false hasContentIssue false

Interface Tracking in Digitally Recorded Glaciological Data

Published online by Cambridge University Press:  09 May 2017

A.P.R. Cooper*
Affiliation:
Scott Polar Research Institute, University of Cambridge, Cambridge CB2 1 ER, U.K.
Rights & Permissions [Opens in a new window]

Abstract

As more data are recorded in digital form the importance of automatically extracting parameters of glaciological significance increases. This paper addresses the problem, with particular reference to tracking bedrock or internal reflectors in digitally recorded Radio Echo Sounding (RES) data. It has been found that the simplest solution to this problem is a “supervised” system, where operator decisions may be added interactively, either on operator command or upon loss of track. Increasing internal decision making within the program may reduce the number of operator interventions required, but is unlikely to eliminate them completely.

Algorithms are presented and discussed for determining the position of the interface, for predicting the position of the interface in successive records, for determining the loss of track condition, and for re-acquiring track after loss of track.

Type
Research Article
Copyright
Copyright © International Glaciological Society 1987

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.

Fig. 1 General flow chart of the tracking process

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.

Fig. 2 Examples of typical returned waveforms, illustrating the required tracking 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

(1)

(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

(2)

(3)

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.

Fig. 3 Flow chart of RES tracking program used for data gathered in Svalbard, 1985.

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:

(A1)

(A2)

therefore

(A3)

(A4)

therefore

(A5)

(A6)

therefore

(A7)

therefore

(A8)

therefore

(A9)

In order to correct for the curvature of the Earth, we have to add the further correction SF:

(A10)

but

(A11)

therefore

(A12)

and therefore

(A13)

The total correction, DE + SF, is given by:

(A14)

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.

References

Albert, D G 1986 FORTRAN subroutines for zero-phase digital frequency filters. CRREL Special Report: 864 Google Scholar
Drewry, D J, Liestøl, O 1985 Glaciological investigations of surging ice caps in Nordaustlandet, Svalbard, 1983. Polar Record 22(139): 359378 Google Scholar
Gorman, M R, Cooper, A P R 1987 A digital radio echo-sounding and navigation recording system. Annals of Glaciology 9: 8184 Google Scholar
Harrison, C H 1970 Reconstruction of subglacial relief from radio echo sounding records. Geophysics 35(6): 10991115 Google Scholar
Harrison, C H 1971 Radio echo sounding: focussing effects in wavy strata. Geophysical Journal of the Royal Astronomical Society 24: 383400 Google Scholar
McIntyre, N F, Drewry, D J 1984 Modelling ice-sheet surfaces for ERS-1’s radar altimeter. ESA Journal 8(3): 261274 Google Scholar
Rapley, C G and 22 others 1983 A study of satellite radar altimeter operation over ice-covered surfaces. Paris, European Space Agency (ESA Contract Report 5182/82/F/CG(SC))Google Scholar
Rapley, C G and 16 others 1985 Applications and scientific uses of ERS-1 radar altimeter data: final report. Noordwijk, European Space Agency (ESA Contract Report 5684/83/NL/BI)Google Scholar
Sy, A 1985 An alternative editing technique for oceanographic data. Deep Sea Research 32A(12): 15911599 CrossRefGoogle Scholar
Wingham, D J, Rapley, C G, Griffiths, H 1986 New techniques in satellite altimeter tracking systems. In IGARSS ‘86 Symposium, Zurich, September 1986. Paris, European Space Agency: 13391344 (ESA Special Publication 254)Google Scholar
Figure 0

Fig. 1 General flow chart of the tracking process

Figure 1

Fig. 2 Examples of typical returned waveforms, illustrating the required tracking points.

Figure 2

Fig. 3 Flow chart of RES tracking program used for data gathered in Svalbard, 1985.