PROCESSING OF RETURN BACKSCATTER ENERGY

Methods used for the analysis of waveforms acquired during the second flight of the Shuttle Laser Altimeter (SLA) are summarized here. They represent modifications made to codes developed for SLA-01. SLA waveform processing is implemented in the Interactive Data Language (IDL) environment. The methods are dynamic and continue to be modified as experience is gained. Therefore, the following procedures reflect the current SLA ‘state-of-the-art’ processing.

Parameterizing the Return Signal

Waveform processing algorithms parameterize the return signal resulting from the interaction of the transmitted laser pulse with the intercepted surface, and identify the response from the multiple targets encountered within the footprint. The SLA detector output voltage is continuously sampled by a high speed, 8-bit digitizer. Upon detection of a backscatter return by the ranging electronics, the digitizer time series is sampled and stored, yielding a waveform record of received laser backscatter energy. The digitizer memory is sampled so as to record detector output voltage beginning slightly before the ranging electronic’s detection of the backscatter return and extending in time to include the maximum range of within-footprint heights expected for land surfaces. Thus a time series of the complete backscatter return for land surfaces is recorded. Returns are modeled as a single Gaussian function, or as the combination of several Gaussian peaks when measurement of multiple ranges from a single return is required. In this manner, the vertical extent and approximate height distribution of intercepted surfaces can be derived from the return signal. Most of the waveforms are single peaked and can be fit by a single Gaussian function, characterized by its maximum amplitude, location of this maximum amplitude in time with respect to the ranging electronics detection time, and its half width. When complex surfaces are intercepted within the footprint (as with the presence of complex surface topography, clouds, vegetation, buildings), multiple returns are present in the waveforms and multi-Gaussian functions are used to model these more complex waveforms.

In brief, the waveform processing steps consist of:

  1. Identifying and processing only shots that are classified as valid surface returns form land and ocean based on a comparison of the laser bounce point elevation (orthometric height) to a reference surface (5 minute resolution Terrain Base Digital Elevation Model for land returns, and mean sea level for ocean returns).
  2. Determining the noise baseline and calculating noise mean and standard deviation, establishing a waveform threshold level for signal above noise.
  3. Identifying start and end of signal above waveform threshold.
  4. Identifying saturated returns.
  5. Subtracting mean noise level from the signal.
  6. Characterizing the basic properties of the signal.
  7. Scaling waveform engineering units to physical units (detector output voltage vs. time) based on scaling factors and calibration constants.
  8. Identifying returns to be excluded from processing based on anomalous characteristics.
  9. Smoothing the signal.
  10. Establishing initial estimates for peak positions, amplitudes and half-widths based on the first and second derivatives of the signal, with exception handling for saturated returns.
  11. Applying constrained function fitting to obtain first estimates of peak amplitudes, with exception handling for saturated returns. Editing peaks based on their amplitude and proximity (zero amplitude peaks are eliminated), with exception handling for saturated shots.
  12. Re-evaluating peak’s significance if necessary, and solving for peak amplitude, location, and half-width.
  13. Deriving distances from the start of the waveform signal to: 1) centroid of all peaks, 2) centroid of the last peak.

For our purpose, a Gaussian peak is defined as follows:

, (2)

where . (3)

F is the analytical function representing the model, and the parameters to be solved are:

A(0)= Gaussian peak maximum amplitude

A(1)= location in the time axis

A(2)= 1-sigma deviation from its mean

Xi is the independent variable, Yi is the observations (independent variables), and ERRi are the 1-sigma uncertainties in the observations. The residuals are calculated in the following manner:

Residuals = (Yi - F(Xi)) / ERRI (4)

If ERR are the 1-sigma uncertainties in Y, then the total chi-squared value will be:

(5)

SLA uses the IDL routine MPFIT to fit Gaussian distributions to the waveform. This is a recently-added, user-supplied IDL routine authored by Craig B. Markwardt, NASA/GSFC Code 622, which uses the The Levenberg-Marquardt technique as a particular strategy to iteratively search for the best fit in the parameter space (Bevington and Robinson, 1992). A Gaussian distribution is used because SLA does not digitize the shape of the transmitted pulse. Lacking information on the shape of the transmitted pulse on a per shot basis, a Gaussian distribution is used as a reasonable approximation. A non-Gaussian, user-supplied fitting function can be input to MPFIT. Updated versions can be found on http://astrog.physics.wisc.edu/~craigm/idl/fitting.html#MPFIT. MPFIT uses the Levenberg-Marquardt technique to solve the least-squares problem. Within certain constraints, these routine will find the set of parameters which best fits the data in the least-squares sense; that is, the sum of the weighted squared differences between the model and data is minimized by minimizing the Chi-square value, calculating derivatives numerically via a finite difference approximation. A set of starting parameters are specified, and the user can apply constraints to individual parameters by setting boundaries on the lower and/or upper side. Otherwise, it is assumed that all parameters are free and unconstrained. The step size to be used in calculating the numerical derivatives of the model with respect to the parameters can be defined by the user or it is computed automatically, and the covariance matrix can also be computed.

The following constraints are applied during the fitting search:

1) lower bounding constraint that all return position, amplitude and width parameters be non-negative; 2) upper and lower bounding constraints on each peak position such that its final position can not be outside the estimated position bounded by its half-width (stopping peaks from migrating away into larger adjacent peaks); 3) if detector saturation occurs, only that part of the signal up to the saturation point is fit, the estimated peak position is seeded earlier in time to compensate for the falsely broadened return, and the area under the saturated peak is preserved in the Gaussian fit; 4) if the digitizer 8-bit dynamic range is exceeded causing clipping of the signal, then the signal before and after clipping is fit making up for the clipped part of the signal with reasonable accuracy; 5) shots with more than 10 peaks were excluded as they were difficult to fit, usually lacking convergence after 40 iterations (these are probably low-lying cloud returns mis-classified as ocean or land surface returns).

Deriving Elevations from the Waveform

The SLA ranging electronics provide the range from the start of the transmit pulse to the start of the backscatter return. It is this range that is used in the geolocation processing and, thus, the location of the geolocated bounce point (latitude, longitude, and elevation) refers to the highest detected feature within the 100 m diameter laser footprint. The distances from the start of the waveform signal to the centroid of all the peaks, the centroid of the last peak, and the end of signal, all provided in the SLA-02 data structure, can be used to correct the bounce point elevation depending on the character of the return signal (e.g., single or multi-peaked), the assumed nature of the surface type, and the intent of the end user.

For a measure of the mean elevation of illuminated surfaces within the laser footprint (assuming uniform reflectance of all the surfaces at the 1064 nm laser wavelength) the appropriate range would be from the centroid of the transmit pulse to the centroid of the backscatter return. An approximation of this mean elevation is obtained by:

[geolocated bounce point elevation] + [transmit pulse centroid] – [centroid distance for all return peaks] (6)

The transmit pulse centroid corresponds to the distance from the start to the centroid of the transmit pulse. SLA does not provide a digitized record of the transmit pulse. However, the transmit pulse impulse response can be obtained from returns from flat, smooth surfaces such as water, defining the narrowest possible returns which have not been broadened in time by surface relief. The impulse response pulse width is a function of peak amplitude. Examination of a plot of waveform centroid versus peak amplitude for single-peak returns defines an envelope of data points, with the minimum boundary defining the impulse response centroid which varies from 2 meters for low amplitude returns to 4 meters for high amplitude returns.

For multiple-peaked waveforms where it is assumed that the last waveform peak is due to laser energy backscattered from the ground and that preceding peaks are energy returned from higher surfaces such as vegetation or buildings, the mean elevation of the ground surface can be approximated by:

[geolocated bounce point elevation] + [transmit pulse centroid] – [centroid distance to last peak] (7)

Inference that the last peak corresponds to the ground surface within a 100 m diameter footprint requires that the height distribution is simple, as for example due to an open vegetation canopy or building above a flat ground surface. Ground slope across the footprint can cause convolution of the ground return with returns from overlying surfaces.

For a measure of the lowest detected surface within a footprint, the appropriate range is the distance from the start of the transmit-pulse to the end of the waveform signal, with a correction for the full width of the transmit pulse impulse response. This elevation can be approximated by:

[geolocated bounce point elevation] + 2*[transmit pulse centroid]– [distance to end of waveform signal] (8)

These calculations assume the laser vector is at nadir. For off-nadir pulses, the distances along the laser vector to the centroids and end of signal should be modified by multiplying by the cosine of the off-nadir pointing angle, although this modification is very small for the near-nadir SLA observations.

The waveform processing and resulting products thus provide a means to correct the SLA first return elevation to a mean elevation for illuminated surfaces, a mean elevation of the last return, and the elevation of the lowest return. The appropriate use of these elevations will depend on the assumed character of the surface within the laser footprint and the intent of the user.
 

Bibliography:

Bevington, P.R. and D.K. Robinson, 1992. Data Reduction and Error Analysis for the Physical Sciences.  Mc Graw Hill, Inc., 2nd Edition, pp. 141-167.
 

Documentation extracted from: "Processing of Shuttle Laser Altimeter Range and Return Pulse Data in suppor of SLA-02", by Claudia C. Carabajal, David H. Harding, Scott B. Luthcke, Waipang Fong, Shelley C. Rowton and J.J. Frawley, to be published in the Proceedings of the ISPRS Workshop on Mapping Surface Structure and Topography by Airborn and Spaceborn Lasers, 1999.

Responsible NASA official: David Harding

Technical Contact: Claudia C. Carabajal

Webpage Developer: Claudia Carabajal

Email claudia@stokes.gsfc.nasa.gov with questions, comments or suggestions.

Last modified January 12, 2000.